Source file lbfgs.ml

# 1 "lbfgs.pp.ml"
(* File: lbfgs.ml

   Copyright (C) 2011

     Christophe Troestler <Christophe.Troestler@umons.ac.be>
     WWW: http://math.umons.ac.be/an/software/

   This library is free software; you can redistribute it and/or modify
   it under the terms of the GNU General Public License version 3 or
   later as published by the Free Software Foundation, with the special
   exception on linking described in the file LICENSE.

   This library is distributed in the hope that it will be useful, but
   WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the file
   LICENSE for more details. *)


open Bigarray
open Printf

type 'l vec = (float, float64_elt, 'l) Array1.t
type wvec = fortran_layout vec (* working vectors *)
(* FORTRAN 77 "integer" is mandated to be half the size of DOUBLE PRECISION  *)
type 'l int_vec = (int32, int32_elt, 'l) Array1.t
type wint_vec = fortran_layout int_vec (* working int vectors *)

external setulb :
  n:int ->
  m:int ->
  (* BEWARE: C-style offsets *)
  c_ofsx:int -> x:'l vec ->
  c_ofsl:int -> l:'l vec ->
  c_ofsu:int -> u:'l vec -> nbd:'l int_vec ->
  f:float -> g:'l vec -> factr:float -> pgtol:float ->
  wa:wvec ->      (* dim: depends on FORTRAN version; see [wa_min_size] *)
  iwa:wint_vec -> (* dim: 3n *)
  task:Bytes.t ->  (* length: 60 *)
  iprint:int ->
  csave:Bytes.t -> (* length: 60 (working string) *)
  lsave:wint_vec -> (* logical working array of dimension 4 *)
  isave:wint_vec -> (* dim: 44 *)
  dsave:wvec ->   (* dim: 29 *)
  float
    = "ocaml_lbfgs_setulb_bc" "ocaml_lbfgs_setulb"
(* Return the value of the function 'f'. *)

let max i j = if (i: int) > j then i else j (* specialized version *)
let min i j = if (i: int) < j then i else j

type work = {
  n: int;   (* dimension of the problem used to create this work *)
  wa: wvec;
  iwa: wint_vec;
  task: Bytes.t;
  csave: Bytes.t;
  lsave: wint_vec;
  isave: wint_vec;
  dsave: wvec;
}

let wvec ty n = Array1.create ty fortran_layout n

(* The work space [wa] changes with the L-BFGS-B version. *)

# 73 "lbfgs.pp.ml"
let wa_min_size n m = (
# 73 "lbfgs.pp.ml"
                        (2 * m + 5) 
# 73 "lbfgs.pp.ml"
                              ) * n + (
# 73 "lbfgs.pp.ml"
                                        m * (11 * m + 8) 
# 73 "lbfgs.pp.ml"
                                              )
let wa_n_of_size s m = max 0 ((s - (
# 74 "lbfgs.pp.ml"
                                     m * (11 * m + 8) 
# 74 "lbfgs.pp.ml"
                                           )) / (
# 74 "lbfgs.pp.ml"
                                                  (2 * m + 5) 
# 74 "lbfgs.pp.ml"
                                                        ))

let unsafe_work n m =
  { n = n;
    wa = wvec float64 (wa_min_size n m);
    iwa = wvec int32 (3 * n);
    (* FORTRAN requires the strings to be initialized with spaces: *)
    task = Bytes.make 60 ' ';
    csave = Bytes.make 60 ' ';
    lsave = wvec int32 4;
    isave = wvec int32 44;
    dsave = wvec float64 29;
  }

let work ?(corrections=10) n =
  if corrections <= 0 then
    failwith "Lbfgs.work: corrections must be > 0";
  if n <= 0 then
    failwith "Lbfgs.work: n must be > 0";
  unsafe_work n corrections

(* Check that the work is large enough for the current problem. *)
let check_work n m work =
  if Array1.dim work.wa < wa_min_size n m
    || Array1.dim work.iwa < 3 * n then
    let n_min =
      min (wa_n_of_size (Array1.dim work.wa) m) (Array1.dim work.iwa / 3) in
    failwith(sprintf
               "Lbfgs.min: dim of work too small for problem size n = %i, \
                valid n <= %i" n n_min)

let set_start s =
  (* No final '\000' for FORTRAN *)
  Bytes.set s 0 'S';  Bytes.set s 1 'T';  Bytes.set s 2 'A';
  Bytes.set s 3 'R';  Bytes.set s 4 'T'

exception Abnormal of float * string;;

let is_space c =
  c = ' ' || c = '\t' || c = '\n'

let rec strip_final_spaces s i =
  if i <= 0 then ""
  else if is_space (Bytes.get s i) then strip_final_spaces s (i - 1)
  else Bytes.sub_string s 0 i

let extract_c_string s =
  try strip_final_spaces s (Bytes.index s '\000')
  with Not_found -> strip_final_spaces s (Bytes.length s - 1)

type print =
| No
| Last
| Every of int
| Details
| All
| Full

let int_of_print = function
| No -> -1
| Last -> 0
| Every i ->
  if i <= 0 then -1
  else if i >= 98 then 98
  else i
| Details -> 99
| All -> 100
| Full -> 101

type state = work
(* Distinguish it from the first to avoid questionning a workspace not
   being used.  This information is only available when task=NEW_X. *)

let is_constrained w = w.lsave.{2} <> 0l
let nintervals w = Int32.to_int w.isave.{22}
let nskipped_updates w = Int32.to_int w.isave.{26}
let iter w = Int32.to_int w.isave.{30}
let nupdates w = Int32.to_int w.isave.{31}
let nintervals_current w = Int32.to_int w.isave.{33}
let neval w = Int32.to_int w.isave.{34}
let neval_current w = Int32.to_int w.isave.{36}

let previous_f w = w.dsave.{2}
let norm_dir w = w.dsave.{4}
let eps w = w.dsave.{5}
let time_cauchy w = w.dsave.{7}
let time_subspace_min w = w.dsave.{8}
let time_line_search w = w.dsave.{9}
let slope w = w.dsave.{11}
let normi_grad w = w.dsave.{13}
let slope_init w = w.dsave.{15}

let is_nan x = (x: float) <> x

module F =
struct
  type vec = (float, float64_elt, fortran_layout) Array1.t
  let layout = fortran_layout
  ;;
# 1 "lbfgs_FC.pp.ml"
(* File: lbfgs_FC.ml

   Copyright (C) 2011

     Christophe Troestler <Christophe.Troestler@umons.ac.be>
     WWW: http://math.umons.ac.be/an/software/

   This library is free software; you can redistribute it and/or modify
   it under the terms of the GNU General Public License version 3 or
   later as published by the Free Software Foundation, with the special
   exception on linking described in the file LICENSE.

   This library is distributed in the hope that it will be useful, but
   WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the file
   LICENSE for more details. *)

(* FORTRAN/C code *)

let empty_vec = Array1.create float64 layout 0

let nbd_of_lu n c_ofsl (lopt: vec option) c_ofsu (uopt: vec option) =
  let nbd = Array1.create int32 layout n in
  match lopt, uopt with
  | None, None -> Array1.fill nbd 0l; (empty_vec, empty_vec, nbd)
  | Some l, None ->
    let min_diml = c_ofsl + n in
    if Array1.dim l < min_diml then
      invalid_arg(sprintf "%s.min: dim(l): valid=[%i..[ got=%i"
                    
# 30 "lbfgs_FC.pp.ml"
                     "Lbfgs.F"  
# 30 "lbfgs_FC.pp.ml"
                        min_diml (Array1.dim l));
    for i = 
# 31 "lbfgs_FC.pp.ml"
             1  
# 31 "lbfgs_FC.pp.ml"
                  to 
# 31 "lbfgs_FC.pp.ml"
                      n  
# 31 "lbfgs_FC.pp.ml"
                             do
      if is_nan l.{i} || l.{i} = infinity then
        invalid_arg(sprintf "%s.min: l.{%i} = %f => empty domain" 
# 33 "lbfgs_FC.pp.ml"
                                                                   "Lbfgs.F"  
# 33 "lbfgs_FC.pp.ml"
                                                                      i l.{i});
      nbd.{i} <- if l.{i} = neg_infinity then 0l else 1l
    done;
    (l, empty_vec, nbd)
  | None, Some u ->
    let min_dimu = c_ofsu + n in
    if Array1.dim u < min_dimu then
      invalid_arg(sprintf "%s.min: dim(u): valid=[%i..[ got=%i"
                    
# 41 "lbfgs_FC.pp.ml"
                     "Lbfgs.F"  
# 41 "lbfgs_FC.pp.ml"
                        min_dimu (Array1.dim u));
    for i = 
# 42 "lbfgs_FC.pp.ml"
             1  
# 42 "lbfgs_FC.pp.ml"
                  to 
# 42 "lbfgs_FC.pp.ml"
                      n  
# 42 "lbfgs_FC.pp.ml"
                             do
      if is_nan u.{i} || u.{i} = neg_infinity then
        invalid_arg(sprintf "%s.min: u.{%i} = %f => empty domain" 
# 44 "lbfgs_FC.pp.ml"
                                                                   "Lbfgs.F"  
# 44 "lbfgs_FC.pp.ml"
                                                                      i u.{i});
      nbd.{i} <- if u.{i} = infinity then 0l else 3l
    done;
    (empty_vec, u, nbd)
  | Some l, Some u ->
    let min_diml = c_ofsl + n in
    if Array1.dim l < min_diml then
      invalid_arg(sprintf "%s.min: dim(l): valid=[%i..[ got=%i"
                    
# 52 "lbfgs_FC.pp.ml"
                     "Lbfgs.F"  
# 52 "lbfgs_FC.pp.ml"
                        min_diml (Array1.dim l));
    let min_dimu = c_ofsu + n in
    if Array1.dim u < min_dimu then
      invalid_arg(sprintf "%s.min: dim(u): valid=[%i..[ got=%i"
                    
# 56 "lbfgs_FC.pp.ml"
                     "Lbfgs.F"  
# 56 "lbfgs_FC.pp.ml"
                        min_dimu (Array1.dim u));
    for i = 
# 57 "lbfgs_FC.pp.ml"
             1  
# 57 "lbfgs_FC.pp.ml"
                  to 
# 57 "lbfgs_FC.pp.ml"
                      n  
# 57 "lbfgs_FC.pp.ml"
                             do
      if is_nan l.{i} || l.{i} = infinity then
        invalid_arg(sprintf "%s.min: l.{%i} = %f => empty domain" 
# 59 "lbfgs_FC.pp.ml"
                                                                   "Lbfgs.F"  
# 59 "lbfgs_FC.pp.ml"
                                                                      i l.{i});
      if is_nan u.{i} || u.{i} = neg_infinity then
        invalid_arg(sprintf "%s.min: u.{%i} = %f => empty domain" 
# 61 "lbfgs_FC.pp.ml"
                                                                   "Lbfgs.F"  
# 61 "lbfgs_FC.pp.ml"
                                                                      i u.{i});
      nbd.{i} <-
        if l.{i} = neg_infinity then (if u.{i} = infinity then 0l else 3l)
        else (if u.{i} = infinity then 1l else 2l)
    done;
    (l, u, nbd)

let min ?(print=No) ?work ?nsteps ?stop
    ?(corrections=10) ?(factr=1e7) ?(pgtol=1e-5)
    ?n ?(ofsl=
# 70 "lbfgs_FC.pp.ml"
               1 
# 70 "lbfgs_FC.pp.ml"
                   ) ?l ?(ofsu=
# 70 "lbfgs_FC.pp.ml"
                                1 
# 70 "lbfgs_FC.pp.ml"
                                    ) ?u f_df ?(ofsx=
# 70 "lbfgs_FC.pp.ml"
                                                      1 
# 70 "lbfgs_FC.pp.ml"
                                                          ) (x: vec) =
  if ofsl < 
# 71 "lbfgs_FC.pp.ml"
             1  
# 71 "lbfgs_FC.pp.ml"
                  then invalid_arg(sprintf "%s.min: ofsl < %i" 
# 71 "lbfgs_FC.pp.ml"
                                                                "Lbfgs.F"  
# 71 "lbfgs_FC.pp.ml"
                                                                    1 
# 71 "lbfgs_FC.pp.ml"
                       );
  if ofsu < 
# 72 "lbfgs_FC.pp.ml"
             1  
# 72 "lbfgs_FC.pp.ml"
                  then invalid_arg(sprintf "%s.min: ofsu < %i" 
# 72 "lbfgs_FC.pp.ml"
                                                                "Lbfgs.F"  
# 72 "lbfgs_FC.pp.ml"
                                                                    1 
# 72 "lbfgs_FC.pp.ml"
                       );
  if ofsx < 
# 73 "lbfgs_FC.pp.ml"
             1  
# 73 "lbfgs_FC.pp.ml"
                  then invalid_arg(sprintf "%s.min: ofsx < %i" 
# 73 "lbfgs_FC.pp.ml"
                                                                "Lbfgs.F"  
# 73 "lbfgs_FC.pp.ml"
                                                                    1 
# 73 "lbfgs_FC.pp.ml"
                       );
  (* Convert to C-style offsets *)
  let c_ofsx = ofsx - 
# 75 "lbfgs_FC.pp.ml"
                       1 
  
# 76 "lbfgs_FC.pp.ml"
  and c_ofsl = ofsl - 
# 76 "lbfgs_FC.pp.ml"
                       1 
  
# 77 "lbfgs_FC.pp.ml"
  and c_ofsu = ofsu - 
# 77 "lbfgs_FC.pp.ml"
                       1  
# 77 "lbfgs_FC.pp.ml"
                            in
  let n = match n with
    | None ->
      let max_ofsx = 
# 80 "lbfgs_FC.pp.ml"
                      Array1.dim x  
# 80 "lbfgs_FC.pp.ml"
                                        in
      if ofsx > max_ofsx then
        invalid_arg(sprintf "%s.min: ofsx: valid=[%i..%i] got=%i"
                      
# 83 "lbfgs_FC.pp.ml"
                       "Lbfgs.F"  
# 83 "lbfgs_FC.pp.ml"
                           1  
# 83 "lbfgs_FC.pp.ml"
                                max_ofsx ofsx);
      Array1.dim x - c_ofsx
    | Some n ->
      if n <= 0 then invalid_arg(
# 86 "lbfgs_FC.pp.ml"
                                  "Lbfgs.F"  
# 86 "lbfgs_FC.pp.ml"
                                     ^ ".min: n <= 0");
      let min_dim = c_ofsx + n in
      if Array1.dim x < min_dim then
        invalid_arg(sprintf "%s.min: dim(x): valid=[%i..[ got=%i"
                      
# 90 "lbfgs_FC.pp.ml"
                       "Lbfgs.F"  
# 90 "lbfgs_FC.pp.ml"
                          min_dim (Array1.dim x));
      n in
  if corrections <= 0 then failwith(
# 92 "lbfgs_FC.pp.ml"
                                     "Lbfgs.F"  
# 92 "lbfgs_FC.pp.ml"
                                        ^ ".min: corrections <= 0");
  let l, u, nbd = nbd_of_lu n c_ofsl l c_ofsu u in
  let w = match work with
    | None -> unsafe_work n corrections
    | Some w -> check_work n corrections w; w in
  set_start w.task; (* task = "START" *)
  let continue = ref true in
  let f = ref nan
  and g = Array1.create float64 layout n in
  let stop_at_x = match nsteps, stop with
    | None, None -> (fun _w -> false)
    | Some n, None -> (fun w -> Int32.to_int w.isave.{30} > n)
    | None, Some f -> f
    | Some n, Some f -> (fun w -> Int32.to_int w.isave.{30} > n || f w) in
  while !continue do
    f := setulb ~n ~m:corrections ~c_ofsx ~x ~c_ofsl ~l ~c_ofsu ~u ~nbd
      ~f:!f ~g ~factr ~pgtol
      ~wa:w.wa ~iwa:w.iwa ~task:w.task ~iprint:(int_of_print print)
      ~csave:w.csave ~lsave:w.lsave ~isave:w.isave ~dsave:w.dsave;
    match Bytes.get w.task 0 with
    | 'F' (* FG *) -> f := f_df x g
    | 'C' (* CONV *) ->
      (* the termination test in L-BFGS-B has been satisfied. *)
      continue := false
    | 'A' (* ABNO *) -> raise(Abnormal(!f, extract_c_string w.task))
    | 'E' (* ERROR *) -> invalid_arg (extract_c_string w.task)
    | 'N' (* NEW_X *) -> if stop_at_x w then continue := false
    | _ -> assert false
  done;
  1. *. !f (* unbox f *)


let max ?print ?work ?nsteps ?stop ?corrections ?factr ?pgtol
    ?n ?ofsl ?l ?ofsu ?u f_df ?ofsx (x: vec) =
  (* Play with -f *)
  let neg_f_df (x: vec) (g: vec) =
    let v = f_df x g in
    for i = 
# 129 "lbfgs_FC.pp.ml"
             1  
# 129 "lbfgs_FC.pp.ml"
                  to 
# 129 "lbfgs_FC.pp.ml"
                      Array1.dim g  
# 129 "lbfgs_FC.pp.ml"
                                        do
      g.{i} <- -. g.{i}
    done;
    -. v in
  let v = min ?print ?work ?nsteps ?stop ?corrections ?factr ?pgtol
    ?n ?ofsl ?l ?ofsu ?u neg_f_df ?ofsx x in
  -. v


(* Local Variables: *)
(* compile-command: "make -k -C .." *)
(* End: *)
# 180 "lbfgs.pp.ml"
end

module C =
struct
  type vec = (float, float64_elt, c_layout) Array1.t
  let layout = c_layout
  ;;
# 1 "lbfgs_FC.pp.ml"
(* File: lbfgs_FC.ml

   Copyright (C) 2011

     Christophe Troestler <Christophe.Troestler@umons.ac.be>
     WWW: http://math.umons.ac.be/an/software/

   This library is free software; you can redistribute it and/or modify
   it under the terms of the GNU General Public License version 3 or
   later as published by the Free Software Foundation, with the special
   exception on linking described in the file LICENSE.

   This library is distributed in the hope that it will be useful, but
   WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the file
   LICENSE for more details. *)

(* FORTRAN/C code *)

let empty_vec = Array1.create float64 layout 0

let nbd_of_lu n c_ofsl (lopt: vec option) c_ofsu (uopt: vec option) =
  let nbd = Array1.create int32 layout n in
  match lopt, uopt with
  | None, None -> Array1.fill nbd 0l; (empty_vec, empty_vec, nbd)
  | Some l, None ->
    let min_diml = c_ofsl + n in
    if Array1.dim l < min_diml then
      invalid_arg(sprintf "%s.min: dim(l): valid=[%i..[ got=%i"
                    
# 30 "lbfgs_FC.pp.ml"
                     "Lbfgs.C"  
# 30 "lbfgs_FC.pp.ml"
                        min_diml (Array1.dim l));
    for i = 
# 31 "lbfgs_FC.pp.ml"
             0  
# 31 "lbfgs_FC.pp.ml"
                  to 
# 31 "lbfgs_FC.pp.ml"
                      n - 1  
# 31 "lbfgs_FC.pp.ml"
                             do
      if is_nan l.{i} || l.{i} = infinity then
        invalid_arg(sprintf "%s.min: l.{%i} = %f => empty domain" 
# 33 "lbfgs_FC.pp.ml"
                                                                   "Lbfgs.C"  
# 33 "lbfgs_FC.pp.ml"
                                                                      i l.{i});
      nbd.{i} <- if l.{i} = neg_infinity then 0l else 1l
    done;
    (l, empty_vec, nbd)
  | None, Some u ->
    let min_dimu = c_ofsu + n in
    if Array1.dim u < min_dimu then
      invalid_arg(sprintf "%s.min: dim(u): valid=[%i..[ got=%i"
                    
# 41 "lbfgs_FC.pp.ml"
                     "Lbfgs.C"  
# 41 "lbfgs_FC.pp.ml"
                        min_dimu (Array1.dim u));
    for i = 
# 42 "lbfgs_FC.pp.ml"
             0  
# 42 "lbfgs_FC.pp.ml"
                  to 
# 42 "lbfgs_FC.pp.ml"
                      n - 1  
# 42 "lbfgs_FC.pp.ml"
                             do
      if is_nan u.{i} || u.{i} = neg_infinity then
        invalid_arg(sprintf "%s.min: u.{%i} = %f => empty domain" 
# 44 "lbfgs_FC.pp.ml"
                                                                   "Lbfgs.C"  
# 44 "lbfgs_FC.pp.ml"
                                                                      i u.{i});
      nbd.{i} <- if u.{i} = infinity then 0l else 3l
    done;
    (empty_vec, u, nbd)
  | Some l, Some u ->
    let min_diml = c_ofsl + n in
    if Array1.dim l < min_diml then
      invalid_arg(sprintf "%s.min: dim(l): valid=[%i..[ got=%i"
                    
# 52 "lbfgs_FC.pp.ml"
                     "Lbfgs.C"  
# 52 "lbfgs_FC.pp.ml"
                        min_diml (Array1.dim l));
    let min_dimu = c_ofsu + n in
    if Array1.dim u < min_dimu then
      invalid_arg(sprintf "%s.min: dim(u): valid=[%i..[ got=%i"
                    
# 56 "lbfgs_FC.pp.ml"
                     "Lbfgs.C"  
# 56 "lbfgs_FC.pp.ml"
                        min_dimu (Array1.dim u));
    for i = 
# 57 "lbfgs_FC.pp.ml"
             0  
# 57 "lbfgs_FC.pp.ml"
                  to 
# 57 "lbfgs_FC.pp.ml"
                      n - 1  
# 57 "lbfgs_FC.pp.ml"
                             do
      if is_nan l.{i} || l.{i} = infinity then
        invalid_arg(sprintf "%s.min: l.{%i} = %f => empty domain" 
# 59 "lbfgs_FC.pp.ml"
                                                                   "Lbfgs.C"  
# 59 "lbfgs_FC.pp.ml"
                                                                      i l.{i});
      if is_nan u.{i} || u.{i} = neg_infinity then
        invalid_arg(sprintf "%s.min: u.{%i} = %f => empty domain" 
# 61 "lbfgs_FC.pp.ml"
                                                                   "Lbfgs.C"  
# 61 "lbfgs_FC.pp.ml"
                                                                      i u.{i});
      nbd.{i} <-
        if l.{i} = neg_infinity then (if u.{i} = infinity then 0l else 3l)
        else (if u.{i} = infinity then 1l else 2l)
    done;
    (l, u, nbd)

let min ?(print=No) ?work ?nsteps ?stop
    ?(corrections=10) ?(factr=1e7) ?(pgtol=1e-5)
    ?n ?(ofsl=
# 70 "lbfgs_FC.pp.ml"
               0 
# 70 "lbfgs_FC.pp.ml"
                   ) ?l ?(ofsu=
# 70 "lbfgs_FC.pp.ml"
                                0 
# 70 "lbfgs_FC.pp.ml"
                                    ) ?u f_df ?(ofsx=
# 70 "lbfgs_FC.pp.ml"
                                                      0 
# 70 "lbfgs_FC.pp.ml"
                                                          ) (x: vec) =
  if ofsl < 
# 71 "lbfgs_FC.pp.ml"
             0  
# 71 "lbfgs_FC.pp.ml"
                  then invalid_arg(sprintf "%s.min: ofsl < %i" 
# 71 "lbfgs_FC.pp.ml"
                                                                "Lbfgs.C"  
# 71 "lbfgs_FC.pp.ml"
                                                                    0 
# 71 "lbfgs_FC.pp.ml"
                       );
  if ofsu < 
# 72 "lbfgs_FC.pp.ml"
             0  
# 72 "lbfgs_FC.pp.ml"
                  then invalid_arg(sprintf "%s.min: ofsu < %i" 
# 72 "lbfgs_FC.pp.ml"
                                                                "Lbfgs.C"  
# 72 "lbfgs_FC.pp.ml"
                                                                    0 
# 72 "lbfgs_FC.pp.ml"
                       );
  if ofsx < 
# 73 "lbfgs_FC.pp.ml"
             0  
# 73 "lbfgs_FC.pp.ml"
                  then invalid_arg(sprintf "%s.min: ofsx < %i" 
# 73 "lbfgs_FC.pp.ml"
                                                                "Lbfgs.C"  
# 73 "lbfgs_FC.pp.ml"
                                                                    0 
# 73 "lbfgs_FC.pp.ml"
                       );
  (* Convert to C-style offsets *)
  let c_ofsx = ofsx - 
# 75 "lbfgs_FC.pp.ml"
                       0 
  
# 76 "lbfgs_FC.pp.ml"
  and c_ofsl = ofsl - 
# 76 "lbfgs_FC.pp.ml"
                       0 
  
# 77 "lbfgs_FC.pp.ml"
  and c_ofsu = ofsu - 
# 77 "lbfgs_FC.pp.ml"
                       0  
# 77 "lbfgs_FC.pp.ml"
                            in
  let n = match n with
    | None ->
      let max_ofsx = 
# 80 "lbfgs_FC.pp.ml"
                      Array1.dim x - 1  
# 80 "lbfgs_FC.pp.ml"
                                        in
      if ofsx > max_ofsx then
        invalid_arg(sprintf "%s.min: ofsx: valid=[%i..%i] got=%i"
                      
# 83 "lbfgs_FC.pp.ml"
                       "Lbfgs.C"  
# 83 "lbfgs_FC.pp.ml"
                           0  
# 83 "lbfgs_FC.pp.ml"
                                max_ofsx ofsx);
      Array1.dim x - c_ofsx
    | Some n ->
      if n <= 0 then invalid_arg(
# 86 "lbfgs_FC.pp.ml"
                                  "Lbfgs.C"  
# 86 "lbfgs_FC.pp.ml"
                                     ^ ".min: n <= 0");
      let min_dim = c_ofsx + n in
      if Array1.dim x < min_dim then
        invalid_arg(sprintf "%s.min: dim(x): valid=[%i..[ got=%i"
                      
# 90 "lbfgs_FC.pp.ml"
                       "Lbfgs.C"  
# 90 "lbfgs_FC.pp.ml"
                          min_dim (Array1.dim x));
      n in
  if corrections <= 0 then failwith(
# 92 "lbfgs_FC.pp.ml"
                                     "Lbfgs.C"  
# 92 "lbfgs_FC.pp.ml"
                                        ^ ".min: corrections <= 0");
  let l, u, nbd = nbd_of_lu n c_ofsl l c_ofsu u in
  let w = match work with
    | None -> unsafe_work n corrections
    | Some w -> check_work n corrections w; w in
  set_start w.task; (* task = "START" *)
  let continue = ref true in
  let f = ref nan
  and g = Array1.create float64 layout n in
  let stop_at_x = match nsteps, stop with
    | None, None -> (fun _w -> false)
    | Some n, None -> (fun w -> Int32.to_int w.isave.{30} > n)
    | None, Some f -> f
    | Some n, Some f -> (fun w -> Int32.to_int w.isave.{30} > n || f w) in
  while !continue do
    f := setulb ~n ~m:corrections ~c_ofsx ~x ~c_ofsl ~l ~c_ofsu ~u ~nbd
      ~f:!f ~g ~factr ~pgtol
      ~wa:w.wa ~iwa:w.iwa ~task:w.task ~iprint:(int_of_print print)
      ~csave:w.csave ~lsave:w.lsave ~isave:w.isave ~dsave:w.dsave;
    match Bytes.get w.task 0 with
    | 'F' (* FG *) -> f := f_df x g
    | 'C' (* CONV *) ->
      (* the termination test in L-BFGS-B has been satisfied. *)
      continue := false
    | 'A' (* ABNO *) -> raise(Abnormal(!f, extract_c_string w.task))
    | 'E' (* ERROR *) -> invalid_arg (extract_c_string w.task)
    | 'N' (* NEW_X *) -> if stop_at_x w then continue := false
    | _ -> assert false
  done;
  1. *. !f (* unbox f *)


let max ?print ?work ?nsteps ?stop ?corrections ?factr ?pgtol
    ?n ?ofsl ?l ?ofsu ?u f_df ?ofsx (x: vec) =
  (* Play with -f *)
  let neg_f_df (x: vec) (g: vec) =
    let v = f_df x g in
    for i = 
# 129 "lbfgs_FC.pp.ml"
             0  
# 129 "lbfgs_FC.pp.ml"
                  to 
# 129 "lbfgs_FC.pp.ml"
                      Array1.dim g - 1  
# 129 "lbfgs_FC.pp.ml"
                                        do
      g.{i} <- -. g.{i}
    done;
    -. v in
  let v = min ?print ?work ?nsteps ?stop ?corrections ?factr ?pgtol
    ?n ?ofsl ?l ?ofsu ?u neg_f_df ?ofsx x in
  -. v


(* Local Variables: *)
(* compile-command: "make -k -C .." *)
(* End: *)
# 191 "lbfgs.pp.ml"
end

(* Local Variables: *)
(* mode: tuareg *)
(* compile-command: "make -k -C .." *)
(* End: *)