let expbF (expData : float[]) (x : float[]) (fi : float[]) myobj =
    for i in 0..fi.Length-1 do
        fi.[i] <- ((x.[0] * System.Math.Exp(-x.[1] * float i) + x.[2]) - expData.[i]) / 0.1

let expbDf (x : float[]) (fi : double[]) (jac : double[,]) myobj =
    let a = x.[0]
    let lambda = x.[1]
    for i in 0..fi.Length-1 do
        let t = float i   
        let e = System.Math.Exp(-lambda * t)
        jac.[i, 0] <- e / 0.1
        jac.[i, 1] <- -t * a * e / 0.1
        jac.[i, 2] <- 1.0 / 0.1

let testDV rawExpData =
    let mutable fitExp : float[] = [|1.0; 0.0; 0.0; |]
    let mutable state : alglib.minlmstate  =  (alglib.minlmstate())
    let mutable rep : alglib.minlmreport = alglib.minlmreport()

    alglib.minlmcreatevj(Array.length rawExpData, fitExp, &state)
    alglib.minlmsetcond(state, 0.0, 0.0, 0.0, 0)
    alglib.minlmoptimize(state, alglib.ndimensional_fvec(expbF rawExpData), expbDf, null, null )
    alglib.minlmresults(state, &fitExp, &rep)
    rep.terminationtype, alglib.ap.format(fitExp, 2)

namespace GSL {
#include <gsl/gsl_blas.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_multifit_nlin.h>
#include <gsl/gsl_fit.h>

}

struct data {
   size_t n;
   const double * y;
};

int expb_f(const GSL::gsl_vector * x, void *data, GSL::gsl_vector * f)
{
   size_t n = ((struct data *)data)->n;
   const double *y = ((struct data *)data)->y;

   double A = gsl_vector_get(x, 0);
   double lambda = gsl_vector_get(x, 1);
   double b = gsl_vector_get(x, 2);

   size_t i;

   for (i = 0; i < n; i++)
   {
      /* Model Yi = A * exp(-lambda * i) + b */
      double t = i;
      double Yi = A * exp(-lambda * t) + b;
      gsl_vector_set(f, i, (Yi - y[i]) / 0.1);
   }

   return GSL::GSL_SUCCESS;
}
int expb_df(const GSL::gsl_vector * x, void *data, GSL::gsl_matrix * J)
{
   size_t n = ((struct data *)data)->n;

   double A = GSL::gsl_vector_get(x, 0);
   double lambda = GSL::gsl_vector_get(x, 1);

   size_t i;

   for (i = 0; i < n; i++)
   {
      /* Jacobian matrix J(i,j) = dfi / dxj, */
      /* where fi = (Yi - yi)/sigma[i],      */
      /*       Yi = A * exp(-lambda * i) + b  */
      /* and the xj are the parameters (A,lambda,b) */
      double t = i;
      double s = 0.1;
      double e = exp(-lambda * t);
      gsl_matrix_set(J, i, 0, e / s);
      gsl_matrix_set(J, i, 1, -t * A * e / s);
      gsl_matrix_set(J, i, 2, 1 / s);
   }
   return GSL::GSL_SUCCESS;
}
int expb_fdf(const GSL::gsl_vector * x, void *data, GSL::gsl_vector * f, GSL::gsl_matrix * J)
{
   expb_f(x, data, f);
   expb_df(x, data, J);

   return GSL::GSL_SUCCESS;
}
double* getTestData() {
   const size_t n = 40;
   auto y = new double[n];

   GSL::gsl_rng_env_setup();

   auto type = GSL::gsl_rng_default;
   GSL::gsl_rng * r;
   r = GSL::gsl_rng_alloc(type);
   for (size_t i = 0; i < n; i++)
   {
      double t = i;
      y[i] = 1.0 + 5 * exp(-0.1 * t) + gsl_ran_gaussian(r, 0.1);
   };
   return y;
}
int main(void)
{
   double* const y = ExpFit::getTestData();

   const GSL::gsl_multifit_fdfsolver_type * T = GSL::gsl_multifit_fdfsolver_lmsder;
   GSL::gsl_multifit_fdfsolver * s = gsl_multifit_fdfsolver_alloc(T, 40, 3);
   struct data d = { 40, y };
   GSL::gsl_multifit_function_fdf f;
   f.f = &expb_f;
   f.df = &expb_df;
   f.fdf = &expb_fdf;
   f.n = 40;
   f.p = 3;
   f.params = &d;

   double x_init[3] = { 1.0, 0.0, 0.0 };
   GSL::gsl_vector_view x = GSL::gsl_vector_view_array(x_init, 3);
   gsl_multifit_fdfsolver_set(s, &f, &x.vector);

   int status;
   unsigned int iter = 0;
   do
   {
      iter++;
      status = gsl_multifit_fdfsolver_iterate(s);

      if (status)
      {
         break;
      }

      status = gsl_multifit_test_delta(s->dx, s->x, 1e-4, 1e-4);
   } while (status == GSL::GSL_CONTINUE && iter < 500);
   gsl_multifit_fdfsolver_free(s);
        delete[](y)
   return 0;
}

let expbDf (expData : float[]) (x : float[]) (fi : double[]) (jac : double[,]) myobj =
    let a = x.[0]
    let lambda = x.[1]
    for i in 0..fi.Length-1 do
        fi.[i] <- ((x.[0] * System.Math.Exp(-x.[1] * float i) + x.[2]) - expData.[i]) / 0.1
    for i in 0..fi.Length-1 do
        let t = float i   
        let e = System.Math.Exp(-lambda * t)
        jac.[i, 0] <- e / 0.1
        jac.[i, 1] <- -t * a * e / 0.1
        jac.[i, 2] <- 1.0 / 0.1