dA/dt = -k1*A
dB/dt = k1*A-k2*B-k3*B
dC/dt = k3*B-k4*C

#include <stdio.h>
#include "lmmin.h"
#include "libalglib/diffequations.h"

using namespace alglib;

/* Model for alglib ODE solver */
void
ode_model(const real_1d_array &x, double t, real_1d_array &dxdt, void *ptr)
{
    double* k = (double*) ptr;
    // Debug code:
    /*printf("t = %.5f\n", t);
    printf("k = %.5f\n", k[0]);
    printf("--\n");*/
    dxdt[0] = -k[0]*x[0];
    dxdt[1] = k[0]*x[0];
}

/* data structure to transmit arrays and fit model */
typedef
struct {
    double *t; // Array of timepoints
    double *x; // Array of datapoints (experimental)
    double *sx; // Array of known start values
} data_struct;

/*
  Evaluation of parameters.
  Calls alglibs ODE solver to generate data with estimated parameters.
  Returns the difference of experimental and simulated timepoints.
*/
void
evaluate_model(
  const double *par, int m_dat, const void *data,
  double *fvec, int *info)
{
    // ODE solver routine
    /* for readability, explicit type conversion */
    data_struct *D;
    D = (data_struct*)data;

    // Initial Values for y
    real_1d_array y;
    y.setcontent(2,D->sx);

    // Time points
    real_1d_array t;
    t.setcontent(m_dat/2,D->t);
    double k[] = {par[0]};

    double eps = 0.000001; // Error tolerance
    double h = 0.001; // Step length

    ae_int_t m; // Number of result values
    real_1d_array ttbl; // Array of time values
    real_2d_array ytbl; // Matrix of result data

    odesolverreport rep;
    odesolverstate stt;
    odesolverrkck(y, t, eps, h, stt);

    odesolversolve(stt, ode_model, &k);
    odesolverresults(stt, m, ttbl, ytbl, rep);


    // calculate difference to experimental data
    for (int i = 0; i < m_dat/2; i++ ) {
        fvec[i] = D->x[i] - ytbl[i][0];
        fvec[i+m_dat/2] = D->x[i+m_dat/2] - ytbl[i][1];
    }
}

int
main()
{
    // Array of initial parameter guess
    double par[] = {0.1}; // Use 3 as initial guess for strange behaviour
    // Number of parameters to estimate
    int n_par = (int)(sizeof(par)/sizeof(*par));

    // Array of starting values
    double sx[2] = {10,0};
    // Array of time values
    double t[11]  = {0,1,2,3,4,5,6,7,8,9,10};

    // Number of datapoints
    int m_dat = (int)(sizeof(t)/sizeof(*t))*(sizeof(sx)/sizeof(*sx));

    // Array of testdata
    double x[22] = {10,8.24,6.09,4.95,3.95,2.65,2.08,1.91,1.26,1.03,0.91,0,2.26,4.05,5.49,6.19,7.44,8.46,8.43,9.05,9.13,8.41};

    // Data struct to send data to lmmin
    data_struct data = { t, x, sx };

    // Initialization of lmmin parameters
    lm_status_struct status;
    lm_control_struct control = lm_control_double;

    control.verbosity = 0; // No output during parameter estimation

    // perform the fit
    lmmin( n_par, par, m_dat, (const void*) &data,
           evaluate_model, &control, &status );

    /* print results */
    printf( "\nResults:\n" );
    printf( "status after %d function evaluations:\n  %i = %s\n",
            status.nfev, status.outcome, lm_infmsg[status.outcome] );

    printf("obtained parameters:\n");
    int i;
    for ( i=0; i<n_par; ++i )
    printf("  par[%i] = %12g\n", i, par[i]);
    printf("obtained norm:\n  %12g\n", status.fnorm );

    return 0;
}