type
  PFitData=^TFitData;
  TFitData=record
    t:Pdouble;
    y:Pdouble;
    n:Integer;
  end;
 
function market(alpha,kappa, t:double):double;
begin
  result:=alpha*exp(-kappa*t);
end;
 
function func_f(x: Pgsl_vector; params: Pointer; f: Pgsl_vector): Integer;cdecl;
var
  data:PFitData;
  alpha:double;
  kappa:double;
  i:Integer;
  t,yi,y:double;
begin
  data:= PFitData(params);
  alpha := gsl_vector_get(x, 0);
  kappa := gsl_vector_get(x, 1);
 
  for I := 0 to data^.n-1 do
  begin
    t := data^.t[i];
    yi := data^.y[i];
    y := market(alpha,kappa, t);
 
    gsl_vector_set(f, i, yi - y);
  end;
 
 
  Result:= Ord(GSL_SUCCESS);
end;
 
 
procedure TFormMain.btnASClick(Sender: TObject);
var
  fdf:gsl_multifit_nlinear_fdf;
 
  n:Integer;
 
  p:Integer;
 
  x:pgsl_vector;
  fdf_params:gsl_multifit_nlinear_parameters;
  r:pgsl_rng;
 
  i:Integer;
 
  t:double;
  fit_data:TFitData;
 
  alpha,kappa:double;
 
  tp:pgsl_multifit_nlinear_type;
  gsl_multifit_nlinear_trs_lmaccel2:function():Pgsl_multifit_nlinear_trs;cdecl;
 
  gsl_vector_set2:procedure (v: Pgsl_vector; i: NativeUInt; x: Double); cdecl;
begin
    n:=300; //number of data points to fit
    p:=2; //parameters count
 
 
  x := gsl_vector_alloc(p);
 
  fdf_params :=gsl_multifit_nlinear_default_parameters();
 
  GetMem(fit_data.t, n * sizeof(double));
  GetMem(fit_data.y, n * sizeof(double));
  fit_data.n := n;
 
  //generate simulate data
  for I := 0 to n-1 do
  begin
    t:=i/n;
    fit_data.t[i]:=t;
    fit_data.y[i]:=market(3,4, t);
  end;
 
  //  /* define function to be minimized */
  FillChar(fdf, sizeof(gsl_multifit_nlinear_fdf),0);
  fdf.f := @func_f;
  fdf.df := nil;//func_df;
  fdf.fvv := nil;//func_fvv;
  fdf.n := n;
  fdf.p := p;
  fdf.params := @fit_data;
 
//  /* starting point */
  gsl_vector_set(x, 0, 1.0);
  gsl_vector_set(x, 1, 1.0);
 
  fdf_params.trs := gsl_multifit_nlinear_trs_lmaccel;
 
 
//  solve_system(x, &fdf, &fdf_params);
 
  tp:=gsl_multifit_nlinear_trust;
 
  var max_iter:Integer := 200;
  var xtol:double := 1.0e-8;
  var gtol:double := 1.0e-8;
  var ftol:double := 1.0e-8;
 
  //fdf_params.trs^.alloc(@fdf,n,p);
  var work:pgsl_multifit_nlinear_workspace := gsl_multifit_nlinear_alloc(tp, @fdf_params, n, p);
 
  var f2:pgsl_vector := gsl_multifit_nlinear_residual(work);
  var y:pgsl_vector := gsl_multifit_nlinear_position(work);
 
  var info:integer;
  var chisq0, chisq, rcond:double;
 
  ///* initialize solver */
  gsl_multifit_nlinear_init(x, @fdf, work);
 
  //* store initial cost */
  gsl_blas_ddot(f2, f2, @chisq0);
 
  //* iterate until convergence */
  gsl_multifit_nlinear_driver(max_iter, xtol, gtol, ftol,
                              nil, nil, @info, work);
 
  //* store final cost */
  gsl_blas_ddot(f2, f2, @chisq);
 
  //* store cond(J(x)) */
  gsl_multifit_nlinear_rcond(@rcond, work);
 
  gsl_vector_memcpy(x, y);
 
  gsl_multifit_nlinear_free(work);
 
    alpha := gsl_vector_get(x, 0);
    kappa := gsl_vector_get(x, 1);
 
 
  gsl_vector_free(x);
 
end;