10
1
mirror of https://gitlab.com/scemama/QCaml.git synced 2024-10-30 10:48:12 +01:00
QCaml/Notebooks/SpVec.ipynb

50 KiB

None <html> <head> </head>
In [1]:
#use "topfind";;
#require "jupyter.notebook";;
#require "lacaml.top" ;;
- : unit = ()
Findlib has been successfully loaded. Additional directives:
  #require "package";;      to load a package
  #list;;                   to list the available packages
  #camlp4o;;                to load camlp4 (standard syntax)
  #camlp4r;;                to load camlp4 (revised syntax)
  #predicates "p,q,...";;   to set these predicates
  Topfind.reset();;         to force that packages will be reloaded
  #thread;;                 to enable threads

- : unit = ()
/home/scemama/qp2/external/opam/default/lib/bytes: added to search path
/home/scemama/qp2/external/opam/default/lib/base64: added to search path
/home/scemama/qp2/external/opam/default/lib/base64/base64.cma: loaded
/home/scemama/qp2/external/opam/default/lib/ocaml/compiler-libs: added to search path
/home/scemama/qp2/external/opam/default/lib/ocaml/compiler-libs/ocamlcommon.cma: loaded
/home/scemama/qp2/external/opam/default/lib/result: added to search path
/home/scemama/qp2/external/opam/default/lib/result/result.cma: loaded
/home/scemama/qp2/external/opam/default/lib/ppx_deriving/runtime: added to search path
/home/scemama/qp2/external/opam/default/lib/ppx_deriving/runtime/ppx_deriving_runtime.cma: loaded
/home/scemama/qp2/external/opam/default/lib/ppx_deriving_yojson/runtime: added to search path
/home/scemama/qp2/external/opam/default/lib/ppx_deriving_yojson/runtime/ppx_deriving_yojson_runtime.cma: loaded
/home/scemama/qp2/external/opam/default/lib/ocaml/unix.cma: loaded
/home/scemama/qp2/external/opam/default/lib/uuidm: added to search path
/home/scemama/qp2/external/opam/default/lib/uuidm/uuidm.cma: loaded
/home/scemama/qp2/external/opam/default/lib/easy-format: added to search path
/home/scemama/qp2/external/opam/default/lib/easy-format/easy_format.cma: loaded
/home/scemama/qp2/external/opam/default/lib/biniou: added to search path
/home/scemama/qp2/external/opam/default/lib/biniou/biniou.cma: loaded
/home/scemama/qp2/external/opam/default/lib/yojson: added to search path
/home/scemama/qp2/external/opam/default/lib/yojson/yojson.cma: loaded
/home/scemama/qp2/external/opam/default/lib/jupyter: added to search path
/home/scemama/qp2/external/opam/default/lib/jupyter/jupyter.cma: loaded
/home/scemama/qp2/external/opam/default/lib/jupyter/notebook: added to search path
/home/scemama/qp2/external/opam/default/lib/jupyter/notebook/jupyter_notebook.cma: loaded
/home/scemama/qp2/external/opam/default/lib/ocaml/compiler-libs/ocamlbytecomp.cma: loaded
/home/scemama/qp2/external/opam/default/lib/ocaml/compiler-libs/ocamltoplevel.cma: loaded
/home/scemama/qp2/external/opam/default/lib/ocaml/bigarray.cma: loaded
/home/scemama/qp2/external/opam/default/lib/lacaml: added to search path
/home/scemama/qp2/external/opam/default/lib/lacaml/lacaml.cma: loaded
/home/scemama/qp2/external/opam/default/lib/lacaml/top: added to search path
/home/scemama/qp2/external/opam/default/lib/lacaml/top/lacaml_top.cma: loaded

Sparse Vector module


Types

In [2]:
module L = Lacaml.D

type t =
  {
    dim: int ;
    nnz: int ;
    indices: (int, Bigarray.int_elt, Bigarray.fortran_layout) Bigarray.Array1.t ; (* Indices *)
    values: L.Vec.t
  }
Out[2]:
module L = Lacaml.D
Out[2]:
type t = {
  dim : int;
  nnz : int;
  indices :
    (int, Bigarray.int_elt, Bigarray.fortran_layout) Bigarray.Array1.t;
  values : L.Vec.t;
}

Printers

In [3]:
let pp ppf t =
  let pp_data ppf t =
    for i=1 to t.nnz do
      Format.fprintf ppf "@[(%d,@ %g)@]@;" t.indices.{i} t.values.{i}
    done
  in
  Format.fprintf ppf "@[{@[dim:@ %d@]@;@[%a@]}@]" t.dim pp_data t
Out[3]:
val pp : Format.formatter -> t -> unit = <fun>
In [33]:
#install_printer pp

Creators

In [4]:
let make0 dim =
  { dim ;
    nnz = 0;
    indices = Bigarray.(Array1.create int fortran_layout) 32 ;
    values = L.Vec.create 32;
  }
Out[4]:
val make0 : int -> t = <fun>
In [5]:
let of_vec ?(threshold=0.) v =
  let dim = L.Vec.dim v in
  let buffer_idx = Bigarray.(Array1.create int     fortran_layout) dim in
  let buffer_val = Bigarray.(Array1.create float64 fortran_layout) dim in
  let check = 
    if threshold = 0. then
      fun x -> x <> 0.
    else
      fun x -> (abs_float x) > 0.
  in
  let rec aux k i =
    if i > dim then
      k-1
    else if check v.{i} then
      ( buffer_idx.{k} <- i ;
        buffer_val.{k} <- v.{i} ;
        aux (k+1) (i+1)
      )
    else
      aux k (i+1)
  in
  let nnz = aux 1 1 in
  let indices = Bigarray.(Array1.create int fortran_layout) nnz in
  let values  = L.Vec.create nnz in
  for i=1 to nnz do
    indices.{i} <- buffer_idx.{i};
    values.{i} <- buffer_val.{i};
  done ;
  { dim ; nnz ; indices ; values }
Out[5]:
val of_vec : ?threshold:float -> Lacaml.D.vec -> t = <fun>
In [6]:
let of_array ?(threshold=0.) a =
  L.Vec.of_array a
  |> of_vec ~threshold
Out[6]:
val of_array : ?threshold:float -> float array -> t = <fun>
In [7]:
let copy t =
  let indices = 
    Bigarray.(Array1.create int fortran_layout) t.nnz
  in
  Bigarray.Array1.blit t.indices indices ;
  let values = L.copy t.values in
  { dim = t.dim ;
    nnz = t.nnz ;
    indices ; values }
Out[7]:
val copy : t -> t = <fun>
In [8]:
let map ?(threshold=0.) f t =
  let indices = 
    Bigarray.(Array1.create int fortran_layout) t.nnz
  in
  let check =
    if threshold = 0. then
      fun x -> x <> 0.
    else
      fun x -> abs_float x > threshold
  in
  let values = L.Vec.create t.nnz in
  let nnz = ref 0 in
  for i=1 to t.nnz do
    let w = f t.values.{i} in
    if check w then (
      incr nnz;
      values.{!nnz} <- w ;
      indices.{!nnz} <- t.indices.{i}
    )
  done;
  { dim = t.dim ;
    nnz = !nnz ;
    indices ; values }
Out[8]:
val map : ?threshold:float -> (float -> float) -> t -> t = <fun>

Test

In [9]:
let x = make0 10 

let dense_a = 
   L.Vec.of_array [| 1. ; -2. ; 0. ; 0. ; 0.5 ; 1.e-8 ; 0. ; 3. ; 0. |]
   
let sparse_a = of_vec dense_a 

let _ =
    copy sparse_a = sparse_a
    
let _ =
    copy sparse_a == sparse_a
Out[9]:
val x : t =
  {dim = 10; nnz = 0; indices = <abstr>;
   values =
              R1           R2 R3              R30          R31          R32
    6.91705E-310 6.91705E-310  0 ... 6.91705E-310 6.91705E-310 6.91705E-310}
Out[9]:
val dense_a : Lacaml.D.vec = R1 R2 R3     R7 R8 R9
                              1 -2  0 ...  0  3  0
Out[9]:
val sparse_a : t =
  {dim = 9; nnz = 5; indices = <abstr>;
   values = R1 R2  R3    R4 R5
             1 -2 0.5 1E-08  3}
Out[9]:
- : bool = true
Out[9]:
- : bool = false

Accessors

In [10]:
let dim t = t.dim
let nnz t = t.nnz
let indices t = t.indices
let values  t = t.values
let density t = float_of_int t.nnz /. float_of_int t.dim
Out[10]:
val dim : t -> int = <fun>
Out[10]:
val nnz : t -> int = <fun>
Out[10]:
val indices :
  t -> (int, Bigarray.int_elt, Bigarray.fortran_layout) Bigarray.Array1.t =
  <fun>
Out[10]:
val values : t -> L.Vec.t = <fun>
Out[10]:
val density : t -> float = <fun>
In [11]:
let get t i =
  if i < 1 || i > dim t then invalid_arg "index out of bounds";
                                                                                              
  let rec binary_search index value low high =
    if high = low then
      if index.{low} = value then
        low
      else
        raise Not_found
    else let mid = (low + high) / 2 in
      if index.{mid} > value then
        binary_search index value low (mid - 1)
      else if index.{mid} < value then
        binary_search index value (mid + 1) high
      else
        mid
  in
  try
    let k = 
      let id = indices t in
      binary_search id i id.{1} (nnz t)
    in
    t.values.{k}
  with Not_found -> 0.
Out[11]:
val get : t -> int -> float = <fun>
In [12]:
let iter f t =
  for k=1 to nnz t do
    f t.indices.{k} t.values.{k}
  done
Out[12]:
val iter : (int -> float -> 'a) -> t -> unit = <fun>

Test

In [13]:
dense_a = L.Vec.init (dim sparse_a) (get sparse_a);;
iter (fun i v -> Printf.printf "%d %f\n%!" i v) sparse_a;;
Out[13]:
- : bool = true
1 1.000000
2 -2.000000
5 0.500000
6 0.000000
8 3.000000
Out[13]:
- : unit = ()

Converters

In [14]:
let to_assoc_list t = 
  let rec aux k accu =
    if k = 0 then
      accu
    else
      aux (k-1) ( (t.indices.{k}, t.values.{k})::accu )
  in
  aux (nnz t) []
  
  
let to_vec t =
  let result = L.Vec.make0 (dim t) in
  iter (fun k v -> result.{k} <- v) t;
  result
Out[14]:
val to_assoc_list : t -> (int * float) list = <fun>
Out[14]:
val to_vec : t -> Lacaml.D.vec = <fun>

Test

In [44]:
to_assoc_list sparse_a;;
to_vec sparse_a = dense_a;;
Out[44]:
- : (int * float) list = [(1, 1.); (2, -2.); (5, 0.5); (6, 1e-08); (8, 3.)]
Out[44]:
- : bool = true

Operations

One-vector operations

In [45]:
let immutable f t =
  let result = copy t in
  f result;
  result
  
let scale_mut x t = 
  L.scal x t.values 

let scale x = immutable @@ scale_mut x
  
let neg t =
  { t with values = L.Vec.neg t.values }
  
let norm t =
  L.nrm2 t.values
Out[45]:
val immutable : (t -> 'a) -> t -> t = <fun>
Out[45]:
val scale_mut : float -> t -> unit = <fun>
Out[45]:
val scale : float -> t -> t = <fun>
Out[45]:
val neg : t -> t = <fun>
Out[45]:
val norm : t -> float = <fun>

Test

In [47]:
let sparse_b = copy sparse_a;;
scale_mut 0.5 sparse_b;;

equal sparse_b @@ map (fun x -> x *. 0.5) sparse_a;;


let dense_b = L.copy dense_a;;
L.scal 0.5 dense_b;;

let sparse_c = scale 0.5 sparse_a;;

equal sparse_b sparse_c;;
    
let _ =
    let n1 = neg sparse_a in
    neg n1 = sparse_a
Out[47]:
val sparse_b : t = {dim: 9 (1, 1) (2, -2) (5, 0.5) (6, 1e-08) (8, 3) }
Out[47]:
- : unit = ()
Out[47]:
- : bool = true
Out[47]:
val dense_b : L.vec = R1 R2 R3     R7 R8 R9
                       1 -2  0 ...  0  3  0
Out[47]:
- : unit = ()
Out[47]:
val sparse_c : t = {dim: 9 (1, 0.5) (2, -1) (5, 0.25) (6, 5e-09) (8, 1.5) }
Out[47]:
- : bool = true
Out[47]:
- : bool = true

Two-vector operations

In [19]:
let equal x y =
  (x.nnz + y.nnz = 0) || (
  (x.dim = y.dim) &&
  (x.nnz  = y.nnz) &&
  Bigarray.Array1.(sub x.indices 1 x.nnz = sub y.indices 1 y.nnz) &&
  (Array.sub (L.Vec.to_array x.values) 0 (x.nnz-1) = Array.sub (L.Vec.to_array y.values) 0 (y.nnz-1) )
  )
Out[19]:
val equal : t -> t -> bool = <fun>

aX + Y

Run along all the entries of X and Y simultaneously with indices k and l. m is the index of the new array.

if k<l, update using a*x[k].

if k>l, update using y[l].

if k=l, update using a*x[k] + y[l].

In [20]:
let axpy ?(threshold=0.) ?(alpha=1.) x y =

  if dim x <> dim y then                                                                      
    invalid_arg "Inconsistent dimensions";

  let check = (* Test if value should be added wrt threshold *)
    if threshold = 0. then
      fun x -> x <> 0.
    else
      fun x -> abs_float x > threshold
    in
    
  let f = (* if a=1 in ax+y, then do x+y. If a=0 then do y *)
    if alpha = 1. then
        fun x y -> x +. y
    else if alpha = 0. then
        fun _ y -> y
    else
        fun x y -> alpha *. x +. y
  in
  
  let dim = dim x in
  let nnz = x.nnz + y.nnz in
  let new_indices =  Bigarray.(Array1.create int fortran_layout) nnz in
  let new_values = L.Vec.create nnz in
  
  let rec aux k l m =
    match k <= x.nnz, l <= y.nnz with
    | true , true  -> (* Both arrays are running *)
      begin
        if x.indices.{k} < y.indices.{l} then (
          let w = f x.values.{k} 0. in
          if check w then (
              new_indices.{m} <- x.indices.{k};
              new_values.{m} <- w ;
              (aux [@tailcall]) (k+1) l (m+1)
          ) else (aux [@tailcall]) (k+1) l m
        )
        else if x.indices.{k} > y.indices.{l} then (
          let w = y.values.{l} in
          if check w then (
              new_indices.{m} <- y.indices.{l};
              new_values.{m} <- w;
              (aux [@tailcall]) k (l+1) (m+1)
          ) else (aux [@tailcall]) k (l+1) m
        )
        else (
          let w = f x.values.{k} y.values.{l} in
          if check w then (
              new_indices.{m} <- x.indices.{k};
              new_values.{m} <- w;
              (aux [@tailcall]) (k+1) (l+1) (m+1)
          ) else (aux [@tailcall]) (k+1) (l+1) m
        )
      end
    | false, true  -> (* Array x is done running *)
      begin
        let m = ref m in
        for i=l to y.nnz do
          let w = y.values.{i} in
          if check w then (
              new_indices.{!m} <- y.indices.{i};
              new_values.{!m}  <- w;
              incr m;
          )
        done; !m
      end
    | true, false  -> (* Array y is done running *)
      begin
        let m = ref m in
        for i=k to x.nnz do
          let w = alpha *. x.values.{i} in
          if check w then (
              new_indices.{!m} <- x.indices.{i};
              new_values.{!m}  <- w;
              incr m;
          )
        done; !m
      end
    | false, false -> (* Both arrays are done    *)
      m
  in
  let nnz = (aux 1 1 1) - 1 in
  { dim ; nnz ;
    indices = new_indices ; values = new_values }
Out[20]:
val axpy : ?threshold:float -> ?alpha:float -> t -> t -> t = <fun>
In [21]:
let add ?(threshold=0.) x y = axpy ~threshold ~alpha:1. x y

let sub ?(threshold=0.) x y = axpy ~threshold ~alpha:(-1.) y x
Out[21]:
val add : ?threshold:float -> t -> t -> t = <fun>
Out[21]:
val sub : ?threshold:float -> t -> t -> t = <fun>

Test

In [22]:
let m_A = L.Vec.of_array [| 1. ; 2. ; 0. ; 0. ; 0.01 ; -2. ; 0. ; -1.e-3 ; 0. ; 0.|]
let m_B = L.Vec.of_array [| 0. ; 1. ; 2. ; 0. ; 0. ; 0.01 ; -2. ; 0. ; -1.e-3 ; 2. |]

let m_As = of_vec m_A
let m_Bs = of_vec m_B

let m_C = L.copy m_B ;;
L.axpy ~alpha:2. m_A m_C;;
m_C;;

let m_D = L.copy m_A;;
L.axpy ~alpha:2. m_B m_D;;
m_D;;

let m_Cs = axpy ~alpha:2. m_As m_Bs


let _ =  equal (of_vec m_C) m_Cs;;

let m_Ds = axpy ~alpha:2. m_Bs m_As

let _ = equal (of_vec m_D) m_Ds;;
Out[22]:
val m_A : Lacaml.D.vec = R1 R2 R3         R8 R9 R10
                          1  2  0 ... -0.001  0   0
Out[22]:
val m_B : Lacaml.D.vec = R1 R2 R3     R8     R9 R10
                          0  1  2 ...  0 -0.001   2
Out[22]:
val m_As : t = {dim: 10 (1, 1 (2, 2 (5, 0.01 (6, -2 (8, -0.001 }
Out[22]:
val m_Bs : t = {dim: 10 (2, 1 (3, 2 (6, 0.01 (7, -2 (9, -0.001 (10, 2 }
Out[22]:
val m_C : L.vec = R1 R2 R3     R8     R9 R10
                   0  1  2 ...  0 -0.001   2
Out[22]:
- : unit = ()
Out[22]:
- : L.vec = R1 R2 R3         R8     R9 R10
             2  5  2 ... -0.002 -0.001   2
Out[22]:
val m_D : L.vec = R1 R2 R3         R8 R9 R10
                   1  2  0 ... -0.001  0   0
Out[22]:
- : unit = ()
Out[22]:
- : L.vec = R1 R2 R3         R8     R9 R10
             1  4  4 ... -0.001 -0.002   4
Out[22]:
val m_Cs : t =
  {dim: 10
  (1, 2 (2, 5 (3, 2 (5, 0.02 (6, -3.99 (7, -2 (8, -0.002 (9, -0.001 (10, 2 }
Out[22]:
- : bool = true
Out[22]:
val m_Ds : t =
  {dim: 10
  (1, 1 (2, 4 (3, 4 (5, 0.01 (6, -1.98 (7, -4 (8, -0.001 (9, -0.002 (10, 4 }
Out[22]:
- : bool = true

Dot product

In [23]:
let dot x y =

  let rec aux accu k l =
    if k <= x.nnz && l <= y.nnz then
      begin
        if x.indices.{k} < y.indices.{l} then 
          (aux [@tailcall]) accu (k+1) l
        else if x.indices.{k} > y.indices.{l} then (
          (aux [@tailcall]) accu k (l+1)
        )
        else (
          (aux [@tailcall]) (accu +. x.values.{k} *. y.values.{l}) (k+1) (l+1)
        )
      end
    else
      accu
  in
  aux 0. 1 1
Out[23]:
val dot : t -> t -> float = <fun>

Test

In [24]:
dot sparse_a sparse_b = L.dot dense_a dense_b;;
Out[24]:
- : bool = true

Tests

In [25]:
#require "alcotest";;
/home/scemama/qp2/external/opam/default/lib/astring: added to search path
/home/scemama/qp2/external/opam/default/lib/astring/astring.cma: loaded
/home/scemama/qp2/external/opam/default/lib/cmdliner: added to search path
/home/scemama/qp2/external/opam/default/lib/cmdliner/cmdliner.cma: loaded
/home/scemama/qp2/external/opam/default/lib/seq: added to search path
/home/scemama/qp2/external/opam/default/lib/stdlib-shims: added to search path
/home/scemama/qp2/external/opam/default/lib/stdlib-shims/stdlib_shims.cma: loaded
/home/scemama/qp2/external/opam/default/lib/fmt: added to search path
/home/scemama/qp2/external/opam/default/lib/fmt/fmt.cma: loaded
/home/scemama/qp2/external/opam/default/lib/fmt/fmt_cli.cma: loaded
/home/scemama/qp2/external/opam/default/lib/fmt/fmt_tty.cma: loaded
/home/scemama/qp2/external/opam/default/lib/alcotest: added to search path
/home/scemama/qp2/external/opam/default/lib/alcotest/alcotest.cma: loaded
In [26]:
let x1 = L.Vec.map (fun x -> if abs_float x < 0.6 then 0. else x) (L.Vec.random 100)
and x2 = L.Vec.map (fun x -> if abs_float x < 0.3 then 0. else x) (L.Vec.random 100)                                                                                                            

let x3 = L.Vec.map (fun x -> 2. *. x) x1
and x4 = L.Vec.add x1 x2
and x5 = L.Vec.sub x1 x2
and x6 = 
  let v = L.copy x2 in
  L.axpy ~alpha:3. x1 v;
  v 

  
let v1 = x1
and v2 = x2
and v3 = x3 
and v4 = x4 
and v5 = x5 
and v6 = x6 

  
let v1_s = of_vec x1
and v2_s = of_vec x2
and v3_s = of_vec x3
and v4_s = of_vec x4
and v5_s = of_vec x5
and v6_s = of_vec x6

    
let zero   = L.Vec.make0 100
and zero_s = of_vec (L.Vec.make0 100)
Out[26]:
val x1 : Lacaml.D.vec =
        R1 R2 R3     R98      R99 R100
  0.852078  0  0 ...   0 0.656419    0
val x2 : Lacaml.D.vec =
         R1       R2 R3           R98 R99   R100
  -0.606197 0.411059  0 ... -0.368989   0 0.9001
Out[26]:
val x3 : Lacaml.D.vec =
       R1 R2 R3     R98     R99 R100
  1.70416  0  0 ...   0 1.31284    0
val x4 : Lacaml.D.vec =
        R1       R2 R3           R98      R99   R100
  0.245881 0.411059  0 ... -0.368989 0.656419 0.9001
val x5 : Lacaml.D.vec =
       R1        R2 R3          R98      R99    R100
  1.45827 -0.411059  0 ... 0.368989 0.656419 -0.9001
val x6 : L.vec =
       R1       R2 R3           R98     R99   R100
  1.95004 0.411059  0 ... -0.368989 1.96926 0.9001
Out[26]:
val v1 : Lacaml.D.vec =
        R1 R2 R3     R98      R99 R100
  0.852078  0  0 ...   0 0.656419    0
val v2 : Lacaml.D.vec =
         R1       R2 R3           R98 R99   R100
  -0.606197 0.411059  0 ... -0.368989   0 0.9001
val v3 : Lacaml.D.vec =
       R1 R2 R3     R98     R99 R100
  1.70416  0  0 ...   0 1.31284    0
val v4 : Lacaml.D.vec =
        R1       R2 R3           R98      R99   R100
  0.245881 0.411059  0 ... -0.368989 0.656419 0.9001
val v5 : Lacaml.D.vec =
       R1        R2 R3          R98      R99    R100
  1.45827 -0.411059  0 ... 0.368989 0.656419 -0.9001
val v6 : L.vec =
       R1       R2 R3           R98     R99   R100
  1.95004 0.411059  0 ... -0.368989 1.96926 0.9001
Out[26]:
val v1_s : t =
  {dim: 100
  (1, 0.852078 (4, 0.733763 (7, -0.818022 (10, -0.692304 (11, -0.980468
  (15, 0.974355 (16, -0.665109 (24, 0.889086 (27, 0.736393 (30, -0.843748
  (31, -0.750514 (34, 0.778919 (38, -0.790597 (40, 0.720996 (41, 0.978861
  (42, -0.679388 (43, -0.853596 (46, 0.977185 (47, -0.811648 (49, -0.702925
  (50, 0.89789 (54, -0.675687 (56, 0.752085 (57, -0.626277 (66, 0.664662
  (67, 0.855765 (69, 0.673572 (71, -0.764248 (72, 0.61481 (73, -0.849851
  (74, -0.983965 (84, -0.764701 (87, -0.834976 (92, -0.851806 (94, 0.696202
  (95, 0.844479 (99, 0.656419 }
val v2_s : t =
  {dim: 100
  (1, -0.606197 (2, 0.411059 (4, 0.329122 (5, -0.740514 (7, 0.387312
  (8, 0.977894 (9, -0.43281 (11, -0.971917 (15, -0.471397 (16, 0.863036
  (17, -0.497091 (19, 0.8659 (20, 0.572763 (21, 0.963705 (23, -0.568959
  (24, 0.643522 (25, 0.697603 (26, -0.706857 (29, -0.874151 (30, 0.559791
  (32, -0.392582 (34, -0.488787 (36, -0.988593 (38, -0.817274 (39, -0.78555
  (42, 0.505945 (43, 0.370838 (44, -0.734224 (45, -0.634506 (46, 0.304503
  (47, -0.330526 (48, -0.800156 (49, -0.570671 (50, -0.514998 (51, -0.47082
  (52, 0.731985 (53, -0.501758 (54, -0.856375 (55, 0.422431 (56, 0.819669
  (57, 0.847078 (58, 0.301306 (59, -0.917181 (60, 0.391153 (61, -0.936169
  (62, 0.855439 (63, 0.730839 (64, 0.763087 (65, 0.777606 (66, 0.945573
  (67, -0.584689 (68, 0.312286 (69, 0.840451 (72, 0.980958 (73, 0.383844
  (77, 0.509676 (79, 0.86618 (80, -0.767171 (82, -0.955884 (83, -0.326832
  (88, 0.846553 (89, -0.578475 (90, -0.522823 (91, 0.802248 (92, 0.975437
  (95, -0.816165 (96, -0.999559 (97, 0.738769 (98, -0.368989 (100, 0.9001 }
val v3_s : t =
  {dim: 100
  (1, 1.70416 (4, 1.46753 (7, -1.63604 (10, -1.38461 (11, -1.96094
  (15, 1.94871 (16, -1.33022 (24, 1.77817 (27, 1.47279 (30, -1.6875
  (31, -1.50103 (34, 1.55784 (38, -1.58119 (40, 1.44199 (41, 1.95772
  (42, -1.35878 (43, -1.70719 (46, 1.95437 (47, -1.6233 (49, -1.40585
  (50, 1.79578 (54, -1.35137 (56, 1.50417 (57, -1.25255 (66, 1.32932
  (67, 1.71153 (69, 1.34714 (71, -1.5285 (72, 1.22962 (73, -1.6997
  (74, -1.96793 (84, -1.5294 (87, -1.66995 (92, -1.70361 (94, 1.3924
  (95, 1.68896 (99, 1.31284 }
val v4_s : t =
  {dim: 100
  (1, 0.245881 (2, 0.411059 (4, 1.06289 (5, -0.740514 (7, -0.43071
  (8, 0.977894 (9, -0.43281 (10, -0.692304 (11, -1.95238 (15, 0.502957
  (16, 0.197927 (17, -0.497091 (19, 0.8659 (20, 0.572763 (21, 0.963705
  (23, -0.568959 (24, 1.53261 (25, 0.697603 (26, -0.706857 (27, 0.736393
  (29, -0.874151 (30, -0.283957 (31, -0.750514 (32, -0.392582 (34, 0.290132
  (36, -0.988593 (38, -1.60787 (39, -0.78555 (40, 0.720996 (41, 0.978861
  (42, -0.173443 (43, -0.482758 (44, -0.734224 (45, -0.634506 (46, 1.28169
  (47, -1.14217 (48, -0.800156 (49, -1.2736 (50, 0.382892 (51, -0.47082
  (52, 0.731985 (53, -0.501758 (54, -1.53206 (55, 0.422431 (56, 1.57175
  (57, 0.220801 (58, 0.301306 (59, -0.917181 (60, 0.391153 (61, -0.936169
  (62, 0.855439 (63, 0.730839 (64, 0.763087 (65, 0.777606 (66, 1.61024
  (67, 0.271076 (68, 0.312286 (69, 1.51402 (71, -0.764248 (72, 1.59577
  (73, -0.466006 (74, -0.983965 (77, 0.509676 (79, 0.86618 (80, -0.767171
  (82, -0.955884 (83, -0.326832 (84, -0.764701 (87, -0.834976 (88, 0.846553
  (89, -0.578475 (90, -0.522823 (91, 0.802248 (92, 0.123631 (94, 0.696202
  (95, 0.0283147 (96, -0.999559 (97, 0.738769 (98, -0.368989 (99, 0.656419
  (100, 0.9001 }
val v5_s : t =
  {dim: 100
  (1, 1.45827 (2, -0.411059 (4, 0.404641 (5, 0.740514 (7, -1.20533
  (8, -0.977894 (9, 0.43281 (10, -0.692304 (11, -0.00855015 (15, 1.44575
  (16, -1.52814 (17, 0.497091 (19, -0.8659 (20, -0.572763 (21, -0.963705
  (23, 0.568959 (24, 0.245564 (25, -0.697603 (26, 0.706857 (27, 0.736393
  (29, 0.874151 (30, -1.40354 (31, -0.750514 (32, 0.392582 (34, 1.26771
  (36, 0.988593 (38, 0.0266769 (39, 0.78555 (40, 0.720996 (41, 0.978861
  (42, -1.18533 (43, -1.22443 (44, 0.734224 (45, 0.634506 (46, 0.672682
  (47, -0.481123 (48, 0.800156 (49, -0.132254 (50, 1.41289 (51, 0.47082
  (52, -0.731985 (53, 0.501758 (54, 0.180688 (55, -0.422431 (56, -0.067584
  (57, -1.47335 (58, -0.301306 (59, 0.917181 (60, -0.391153 (61, 0.936169
  (62, -0.855439 (63, -0.730839 (64, -0.763087 (65, -0.777606 (66, -0.280911
  (67, 1.44045 (68, -0.312286 (69, -0.166879 (71, -0.764248 (72, -0.366148
  (73, -1.2337 (74, -0.983965 (77, -0.509676 (79, -0.86618 (80, 0.767171
  (82, 0.955884 (83, 0.326832 (84, -0.764701 (87, -0.834976 (88, -0.846553
  (89, 0.578475 (90, 0.522823 (91, -0.802248 (92, -1.82724 (94, 0.696202
  (95, 1.66064 (96, 0.999559 (97, -0.738769 (98, 0.368989 (99, 0.656419
  (100, -0.9001 }
val v6_s : t =
  {dim: 100
  (1, 1.95004 (2, 0.411059 (4, 2.53041 (5, -0.740514 (7, -2.06675
  (8, 0.977894 (9, -0.43281 (10, -2.07691 (11, -3.91332 (15, 2.45167
  (16, -1.13229 (17, -0.497091 (19, 0.8659 (20, 0.572763 (21, 0.963705
  (23, -0.568959 (24, 3.31078 (25, 0.697603 (26, -0.706857 (27, 2.20918
  (29, -0.874151 (30, -1.97145 (31, -2.25154 (32, -0.392582 (34, 1.84797
  (36, -0.988593 (38, -3.18907 (39, -0.78555 (40, 2.16299 (41, 2.93658
  (42, -1.53222 (43, -2.18995 (44, -0.734224 (45, -0.634506 (46, 3.23606
  (47, -2.76547 (48, -0.800156 (49, -2.67945 (50, 2.17867 (51, -0.47082
  (52, 0.731985 (53, -0.501758 (54, -2.88344 (55, 0.422431 (56, 3.07592
  (57, -1.03175 (58, 0.301306 (59, -0.917181 (60, 0.391153 (61, -0.936169
  (62, 0.855439 (63, 0.730839 (64, 0.763087 (65, 0.777606 (66, 2.93956
  (67, 1.98261 (68, 0.312286 (69, 2.86117 (71, -2.29274 (72, 2.82539
  (73, -2.16571 (74, -2.9519 (77, 0.509676 (79, 0.86618 (80, -0.767171
  (82, -0.955884 (83, -0.326832 (84, -2.2941 (87, -2.50493 (88, 0.846553
  (89, -0.578475 (90, -0.522823 (91, 0.802248 (92, -1.57998 (94, 2.0886
  (95, 1.71727 (96, -0.999559 (97, 0.738769 (98, -0.368989 (99, 1.96926
  (100, 0.9001 }
Out[26]:
val zero : Lacaml.D.vec = R1 R2 R3     R98 R99 R100
                           0  0  0 ...   0   0    0
val zero_s : t = {dim: 100 }
In [27]:
let test_conversion () =
    Alcotest.(check bool) "sparse -> dense  1" true (to_vec v1_s = v1  );
    Alcotest.(check bool) "sparse -> dense  2" true (to_vec v2_s = v2  );
    Alcotest.(check bool) "dense  -> sparse 1" true (of_vec v1   = v1_s);
    Alcotest.(check bool) "dense  -> sparse 2" true (of_vec v2   = v2_s)
 
  
  let test_operations () =
    Alcotest.(check  bool)  "dense   scale"        true  (let w = L.copy v1 in L.scal 2. w ; w  =  v3);
    Alcotest.(check  bool)  "sparse  scale"        true  (scale  2.    v1_s  =  v3_s);
    
    Alcotest.(check  bool)  "dense   dense   add"  true  (L.Vec.add  v1  v2    =  v4);
    (*
    Alcotest.(check  bool)  "dense   sparse  add"  true  (add    v1    v2_s  =  v4_s);
    Alcotest.(check  bool)  "sparse  dense   add"  true  (add    v1_s  v2    =  v4_s);
    Alcotest.(check  bool)  "sparse  dense   add"  true  (add    v1    v2_s  =  v4_s);
    *)
    Alcotest.(check  bool)  "sparse  sparse  add"  true  (equal (add    v1_s  v2_s)  v4_s);
    
    Alcotest.(check  bool)  "dense   dense   sub"  true  (L.Vec.sub    v1    v2    =  v5);
    (*
    Alcotest.(check  bool)  "dense   sparse  sub"  true  (sub    v1    v2_s  =  v5_s);
    Alcotest.(check  bool)  "sparse  dense   sub"  true  (sub    v1_s  v2    =  v5_s);
    Alcotest.(check  bool)  "sparse  dense   sub"  true  (sub    v1    v2_s  =  v5_s);
    *)
    Alcotest.(check  bool)  "sparse  sparse  sub"  true  (equal (sub v1_s v2_s) v5_s);
    
    Alcotest.(check  bool)  "dense   dense   sub"  true  (L.Vec.sub  v1 v1   =  zero);
    (*
    Alcotest.(check  bool)  "dense   sparse  sub"  true  (sub    v1    v1_s  =  zero_s);
    Alcotest.(check  bool)  "sparse  dense   sub"  true  (sub    v1_s  v1    =  zero_s);
    *)
    Alcotest.(check  bool)  "sparse  sparse  sub"  true  (equal (sub v1_s v1_s) zero_s);
    
    Alcotest.(check  bool)  "dense   dense   axpy"  true  (let w = L.copy v2 in L.axpy ~alpha:3. v1  w ;  w  =  v6);
    (*
    Alcotest.(check  bool)  "dense   sparse  axpy"  true  (sub ~threshold:1.e-12 (axpy ~alpha:3.  v1    v2_s) v6_s = zero_s);
    Alcotest.(check  bool)  "sparse  dense   axpy"  true  (sub ~threshold:1.e-12 (axpy ~alpha:3.  v1_s  v2) v6_s = zero_s);
    *)
    Alcotest.(check  bool)  "sparse  sparse  axpy"  true  (equal (sub ~threshold:1.e-12 (axpy ~alpha:3.  v1_s  v2_s) v6_s) zero_s)
  
let test_dot () =  
    let d1d2 = L.dot x1 x2
    and d1d1 = L.dot x1 x1
    and d2d2 = L.dot x2 x2 
    in
    (*
    Alcotest.(check (float 1.e-10)) "sparse x dense  1" (dot v1_s v2  )  d1d2;
    Alcotest.(check (float 1.e-10)) "sparse x dense  2" (dot v1_s v1  )  d1d1;
    Alcotest.(check (float 1.e-10)) "sparse x dense  3" (dot v2_s v2  )  d2d2;
    Alcotest.(check (float 1.e-10)) "dense  x sparse 1" (dot v1   v2_s)  d1d2;
    Alcotest.(check (float 1.e-10)) "dense  x sparse 2" (dot v1   v1_s)  d1d1;
    Alcotest.(check (float 1.e-10)) "dense  x sparse 3" (dot v2   v2_s)  d2d2;
    *)
    Alcotest.(check (float 1.e-10)) "sparse x sparse 1" (dot v1_s v2_s)  d1d2;
    Alcotest.(check (float 1.e-10)) "sparse x sparse 2" (dot v1_s v1_s)  d1d1;
    Alcotest.(check (float 1.e-10)) "sparse x sparse 3" (dot v2_s v2_s)  d2d2
    
    
let test_case () = 
  [ 
    "Conversion", `Quick, test_conversion;
    "Operations", `Quick, test_operations;
    "Dot product", `Quick, test_dot;
  ]
Out[27]:
val test_conversion : unit -> unit = <fun>
Out[27]:
val test_operations : unit -> unit = <fun>
Out[27]:
val test_dot : unit -> unit = <fun>
Out[27]:
val test_case : unit -> (string * [> `Quick ] * (unit -> unit)) list = <fun>
In [28]:
(*)
Alcotest.run ~argv:[|"ignored"|] "Unit tests" [
    "SpVec", test_case ();
];;
Printf.printf "%!";;
*)
File "[28]", line 1, characters 0-3:
Warning 1: this `(*' is the start of a comment.
Hint: Did you forget spaces when writing the infix operator `( * )'?
</html>