nonlinear fit function binding

["matthieu.dubuget" <matthieu.dubuget@laposte.net>]

Hello,

I'm want to bind one non linear fitting function of MINPACK
(Fortran).

My first attempt seems to work, but...

1/ I used one GNU C extension (nested function definition):
I'd like
to know if any more portable solution exists, with a license
suitable
for a commercial product

2/ I'm still not sure about my use of the CAML* macros (see in
C part,
later)





/* Seen from C, the fitting FORTRAN function I first bound is: */
extern void lmdif1_(void fcn(int *m, int *n, double *x, double
*fvec, int *iflag),
	     int * m, int * n, double * x, double * fvec,
	     double * tol, int * info, int * iwa, double * wa, int *
lwa);

/* where fcn() is a user supplied C function pointer. */





(* From OCaml, I wanted to have something like this: *)
type float1darray =
    (float, Bigarray.float64_elt, Bigarray.c_layout)
Bigarray.Genarray.t

(* Here is the kind of function I want to fit.
   The first argument contains the parameters to adjust.
   The second one is the abscisse *)
type fitable = float1darray -> float -> float

(* Since I also want to be able to use directly C functions: *)
type fitfunc =
    OFitfunc of (string * fitable) (* OCaml function. The first
				   parameter is one string. It
                                   allows to name a registered

                                   callback, and pass it to C
*)

  | CFitfunc of nativeint         (* This is directly a C
pointer to a
			          function written in C. The type of C
			          function is ffittable (see in C
			          part). almost same thing as type
			          fitable *)

(* Here is the interface of the C function *)
external fit :
  fitfunc -> float1darray -> float1darray -> float1darray ->
float -> int
  = "minpack_fit"

(* The fitting function I intend to use:
   x contains absisses of experimental datas
   y contains experimental values
   p0 contains an initial guess of the parameters

   The function returns the info code given by MINPACK function.

   Side effect: p0 datas are modified and contain the adjusted
   parameters.
*)
val minimize : fcn:fitfunc -> x:float1darray -> y:float1darray ->
     p0:float1darray -> ?tol:float -> unit -> int

(* implementation: *)
let minimize ~fcn ~x ~y ~p0 ?(tol=1e-6) () =
  let () = match fcn with
      OFitfunc (nom, fc)  -> Callback.register nom fc;
    | CFitfunc _ -> ()
  in
    fit fcn x y p0 tol


/* And know, the C world: */

/* This is the kind of function pointer I can use. */
typedef double (*ffittable)(double, double *);

CAMLprim value minpack_fit(value func, value x, value y,
                           value params0, value tol){
  CAMLparam5(func, x, y, params0, tol);

  double * my_x, * my_y, *my_params, *fvec, *wa;
  double my_tol;
  int m, n, lwa, info, *iwa;
  ffittable f = NULL;
  value * ofunc = NULL;   /* I'm not sure here. Should I have
used a
  			     CAMLlocal macro */

  /* .

     ...Erased sanity checks for shortening the mail...
     .
  */

  /* my_* are C converted versions of '*' parameters */
  my_x = Data_bigarray_val(x);
  my_y = Data_bigarray_val(y);
  my_params = Data_bigarray_val(params0);

  /* Number of experimental data points */
  m = Bigarray_val(x)->dim[0];

  /* Number of parameters */
  n = Bigarray_val(params0)->dim[0];

  my_tol=Double_val(tol);

  /*.
    ... Allocations of fvec, iwa, wa....
    ... Initialisation lwa....
    .
  */

  if (Tag_val(func) == 0){	/* OFitfunc of (string * fitable) */

    /* Nested function definition: GCC EXTENSION */
    void my_fcn(int *m, int *n, double *params, double *fvec,
int *iflag){

      /* Should any CAML* macro be used here? */
      value baparams;
      int i;

      /* And here? */
      ofunc = caml_named_value(
String_val(Field(Field(func,0), 0)) );
      assert(ofunc);

      baparams =
        alloc_bigarray_dims(BIGARRAY_FLOAT64
|BIGARRAY_C_LAYOUT,1, params, *n);

      /* For each data point, fvec is filled with the distance
between
      computed value and experimental values */

      for (i=0; i<(*m); i++){
	fvec[i] =
          Double_val(
             caml_callback2(*ofunc, baparams,
caml_copy_double(my_x[i]))
          ) - my_y[i];
      }

      /* Should baparams be freed there? And if yes, how? */
    };

   /* MINPACK FORTRAN SUBROUTINE CALL */
   lmdif1_(*my_fcn, &m, &n,
	   my_params, fvec, &my_tol, &info, iwa, wa, &lwa);


  } else {
    if (Tag_val(func) == 1) {	/* CFitfunc of nativeint */
      f = (ffittable) Nativeint_val(Field(func, 0));

      /* Nested function definition: GCC EXTENSION */
      void my_fcn(int *m, int *n, double *params, double
*fvec, int *iflag){
	int i;
	for (i=0; i<(*m); i++){
	  fvec[i] = (*f)(my_x[i], params) - my_y[i];
	}
      };

   /* MINPACK FORTRAN SUBROUTINE CALL */
   lmdif1_(*my_fcn, &m, &n, my_params, fvec, &my_tol, &info,
iwa, wa, &lwa);


    } else {
      assert(0);
    }
  }


  /* free what can be freed...*/
   CAMLreturn(Val_int(info));
}

Salutations

Matthieu Dubuget

Créez votre adresse électronique prenom.nom@laposte.net
1 Go d'espace de stockage, anti-spam et anti-virus intégrés.