nnjmvmv

(* wrapper around the fortan mvndst library from Genz. must be in $HOME/lib as
 * a shared library *)

open Bigarray
open Lacaml.D
open Ctypes
open Util
open Containers

[@@@landmark "auto"]
(*let integration = Landmark.register "integration"*)

let integer = Ctypes_static.int32_t   (* fortran integer type representation *)

let os =
  assert (Sys.os_type = "Unix");
  let ic = Unix.open_process_in "uname" in
  let uname = input_line ic in
  let () = close_in ic in
  if uname = "Darwin" then `MacOs else `Linux

(* get the library *)
let mvnd =
  let filename = Sys.getenv "HOME" ^ "/lib/mvndst." ^
                 (if os = `MacOs then "dylib" else "so") in
  Dl.dlopen ~flags:[Dl.RTLD_LAZY] ~filename

(* the integration routine which gives the integral in value *)
let integrate = Foreign.foreign ~from:mvnd "mvndst_"
    (ptr integer @-> (* N *)
     ptr double  @-> (* LOWER *)
     ptr double  @-> (* UPPER *)
     ptr integer @-> (* INFIN *)
     ptr double  @-> (* CORREL *)
     ptr integer @-> (* MAXPTS *)
     ptr double  @-> (* ABSEPS *)
     ptr double  @-> (* RELEPS *)
     ptr double  @-> (* ERROR *)
     ptr double  @-> (* VALUE *)
     ptr integer @-> (* INFORM *)
     returning void);;

(* c_layout arrays to get the start address *)
let c_make0 n =
  let a = Array1.create float64 c_layout n in
  Array1.fill a 0.; a

let v2c fa =
  fa |> genarray_of_array1
  |> (fun fga -> Genarray.change_layout fga c_layout)
  |> array1_of_genarray

(* arrange a symmetric positive matrix in the way the library expects *)
let lower_vec mat =
  let d = Mat.dim1 mat in
  let v = Vec.create (d * (d - 1) / 2) in
  for i = 1 to d do for j = 1 to i - 1 do
      v.{j + (i-2)*(i-1)/2} <- mat.{i,j} done done;
  v

(* correlation and variances *)
let corr_vars ci =
  let c = U.po_inv ci in                                (* covariance matrix *)
  let vars = Mat.copy_diag c in
  let istds =vars |> Vec.sqrt |> Vec.reci in
  Mat.scal_cols c istds;
  Mat.scal_rows istds c;
  lower_vec c, vars

exception Integration_error of string

let log_gauss_integral
    ?lower ?upper
    ?infin
    ?maxpts
    ?releps ?abseps
    ci =
  let default_or default = Option.get_or ~default in
  let n = Mat.dim1 ci in
  let lower, upper =
    default_or (Vec.make0 n) lower, default_or (Vec.make0 n) upper in
  let () =
    let test vec = if not (Vec.dim vec = n)
      then raise (Invalid_argument "need full starting vector") in
    test lower; test upper; () in
  let infin = default_or
      (let a = Lacaml.Common.create_int32_vec n in
       Array1.fill a 1l; a)                (* default from lower to infinity *)
      infin in
  let maxpts = default_or (100*n) maxpts in
  (*Printf.printf "maxpts %d%!" maxpts;*)
  let releps = default_or (1e-2) releps in
  let abseps = default_or (1e-2) abseps in
  let correl, vars = corr_vars ci in
  let jac_log_det = 0.5 *. (vars |> Vec.log |> Vec.sum) in
  (* pointers *)
  let ip k = allocate integer (Int32.of_int k) in
  let fp f = allocate double f in
  let ap v = bigarray_start array1 (v2c v) in
  let _N = ip n in
  let _LOWER, _UPPER, _INFIN = ap lower, ap upper, ap infin in
  let _CORREL = ap correl in
  let _MAXPTS, _ABSEPS, _RELEPS = ip maxpts, fp abseps, fp releps in
  let _ERROR,  _VALUE, _INFORM = fp 0., fp 0., ip 0 in
  (*Landmark.enter integration;*)
  let integral = integrate                                         (* do it! *)
      _N _LOWER _UPPER _INFIN _CORREL _MAXPTS _ABSEPS _RELEPS
      _ERROR _VALUE _INFORM;
    if (!@_INFORM <> 0l)
    then raise (Integration_error
                  (Printf.sprintf "integration problem %d, error estimate %f"
                     (Int32.to_int !@_INFORM) !@_ERROR));
    !@_VALUE in
  (*Landmark.exit integration;*)
  jac_log_det +. log integral