RE: [Caml-list] nonlinear fit function binding

[<RABIH.ELCHAAR@sgam.com>]

Bigarrays, in my understanding are essentially c structure, with a pointer to data, reference of usages, ...
It is very similar to the gsl vector matrix structure.
You can look at the ocamlgsl binding to see how easily bigarrays are transformed into gsl vector and matrices.
To illustrate my answier, I will try to be more specific, and give some c code snapshots.
Let's says you have your library function called minimizer(...,array1, ...)

array1 is a c double array.
Your goal is to create a c function called caml_wrap_minimizer, which will be called from CAML.
This function will take into argument the function f to minimize, the initial starting point, constraints, target point, ...
View from caml, the minimize function have a signature like (including levmar definitions)


(************************************ wrapper.ml ********************************************)
open Bigarray
type float_vector = (float, float64_elt, c_layout) Array1.t 
type float_matrix = (float, float64_elt, c_layout) Array2.t 

type stopping=
  |Small_Grandient
  |Small_Dp
  |Stop_Itmax
  |Matrix_Singular
  |Error_Reduction_Saturated
  |Small_e_2

type calibrationInfo=
    {initial_e_2:float;
     stopping_e_2:float;
     stopping_J_T_e_inf:float;
     stopping_Dp_2:float;
     stopping_max_J_T_J:float;
     iterations:int;
     stopping_reason:stopping;
     function_evaluation:int;
     jacobian_evalutaion:int;}

external workspace_size: ?m:int -> ?n:int -> ?analyticJacobian:bool -> int =
  "workspace_size"

(** gneral minimization *)
external dlevmar_der: ?tau:float -> ?epsilon1:float -> ?epsilon2:float ->
  ?epsilon3:float ->
    ?f:(float_vector -> float_vector -> unit) ->
      ?j:(float_vector -> float_matrix ->unit) ->
        ?p:float_vector -> ?x:float_vector -> ?itmax:int ->
          ?work:float_vector -> ?covar:float_matrix -> calibrationInfo =
            "dlevmar_der_bytecode" "dlevmar_der_native"


external dlevmar_dif: ?tau:float -> ?epsilon1:float -> ?epsilon2:float ->
  ?epsilon3:float -> ?delta:float ->
    ?f:(float_vector -> float_vector -> unit) ->
      ?p:float_vector -> ?x:float_vector -> ?itmax:int ->
        ?work:float_vector -> ?covar:float_matrix ->
          ?tempp:float_vector -> ?tempx:float_vector -> calibrationInfo =
          "dlevmar_dif_bytecode" "dlevmar_dif_native"

(************************************ wrapper.ml ********************************************)

f is the function that computes x = f(p) and our goal is
starting from p, find the point p* verifying
f(p*)=x

In case we know the jacobian analytically, the jacobian function j is passed.
If not, we have to use numerical jacobian computation.

What you have to keep in mind is that the vector p will be modified in memory by the c function.
I will emphasize a bit more the c code.

(************************************ wrapper.c ********************************************)

static double *Vdata(value  v)
{
    return Bigarray_val(v)->data;
}

static value fg; /* : vector -> vector -> unit */
static double * ctempp, * ctempx;
static value tempp;
static value tempx;

static void wrap_f_dif1(double *p, double *hx, int m, int n, void *adata)
{
  memcpy(Vdata(tempp),p,m*sizeof(double));
  callback2(fg,tempp,tempx);
  memcpy(hx,Vdata(tempx),n*sizeof(double));
}

CAMLprim value dlevmar_dif_native(value tau, value epsilon1, value epsilon2, value epsilon3, value delta,
                                  value f, value p, value x, value itmax, 
                                  value work, value covar, value temp_p, value temp_x)
{
  double info[LM_INFO_SZ];
  double opts[LM_OPTS_SZ];
  CAMLparam5(tau,epsilon1,epsilon2,epsilon3,f);
  CAMLxparam5(p,x,itmax,covar,work);
  CAMLxparam3(temp_p,temp_x,delta);
  CAMLlocal1(res);
  int n;
  int m;
  
  double *cp;
  double *cfvec;
  double *ccovar;
  double *cwork;
  int citmax;
  int ret;


  fg = Field(f,0);
  caml_register_global_root(&fg);

  tempp = Field(temp_p,0);
  ctempp = Vdata(tempp);

  tempx = Field(temp_x,0);
  ctempx = Vdata(tempx);


  caml_register_global_root(&tempp);
  caml_register_global_root(&tempx);

  ...

  ret = dlevmar_dif(
              &wrap_f_dif1,...
   
  ...
  caml_remove_global_root(&tempp);
  caml_remove_global_root(&tempx); 
  caml_remove_global_root(&fg);

}
(************************************ wrapper.c ********************************************)


When you look thoroughly to the wrap_f_dif1 function, you get a clear understanding of what I am trying to do.
After having registered f, temp and tempw passed from caml as global_roots, when I have to call my caml function for minimization, i copy the c pointers of my c function to the bigarray structure behind temp and tempx, and do a callback2 to have the function computed in the caml world.

Hope this helps.

Rabih
________________________________________
De : caml-list-bounces@yquem.inria.fr [mailto:caml-list-bounces@yquem.inria.fr] De la part de Matthieu Dubuget
Envoyé : lundi 4 juin 2007 09:38
À : caml-list@inria.fr
Objet : Re: [Caml-list] nonlinear fit function binding

Thanks for your reply, Rabih.
2007/5/31, RABIH.ELCHAAR@sgam.com <RABIH.ELCHAAR@sgam.com>: 
I did a wrapping of the levmar library.

Unfortunately, since I am to use this in a commercial product, and I was not able (yet ?) to convince my boss to release the sources, I can't use with levmar, wich is GPLed... I may switch to cminpack. But i think this is not related to my current question. 

The approach I chose is to allocate bigarrays from Caml, and pass them to the main C function, along with the function f.
Does it mean that each time f function is called by the minimization, the C wrappers has to copy the params C array into your Bigarray pre-allocated storage place? That's why I choosed to leave the C array untouched and wrap a Bigarray around it  (should not be a big amount of memory). Am I wrong in going this way? 

In the main function, I register f and all the bigarrays as global roots, and then call as many times the c sub function to compute the minimization.

I do not understand this part. f is passed to the main (C) function from OCaml: why do you need to register it as a global root? Same question for the Bigarrays? Does it mean that else, nothing in OCaml world would point to them?  
Feel free to ask for more precisions if needed
Hope this helps.
Surely. I really need to acquire a better understanding of GC and C interface in OCaml. 
Really, I do not understand what happens to those C allocated Bigarrays...

Salutations

Matthieu
Ce message et toutes les pieces jointes (ci-apres le "message") sont confidentiels et etablis a l'intention exclusive de ses destinataires. 
Toute utilisation ou diffusion non autorisee est interdite. 
Tout message electronique est susceptible d'alteration. 
Societe Generale Asset Management et ses filiales declinent toute responsabilite au titre de ce message s'il a ete altere, deforme ou falsifie. 
  
Decouvrez l'offre et les services de Societe Generale Asset Management sur le site www.sgam.fr 
  
                                ******** 
  
This message and any attachments (the "message") are confidential and intended solely for the addressees. 
Any unauthorised use or dissemination is prohibited. 
E-mails are susceptible to alteration. 
Neither Societe Generale Asset Management nor any of its subsidiaries or affiliates shall be liable for the message if altered, changed or falsified. 
 
Find out more about Societe Generale Asset Management's proposal on www.sgam.com