mirror of https://gitlab.com/scemama/QCaml.git
50 KiB
50 KiB
None
<html>
<head>
</head>
</html>
In [1]:
#use "topfind";;
#require "jupyter.notebook";;
#require "lacaml.top" ;;
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]:
Out[2]:
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]:
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]:
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]:
In [6]:
let of_array ?(threshold=0.) a =
L.Vec.of_array a
|> of_vec ~threshold
Out[6]:
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]:
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]:
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]:
Out[9]:
Out[9]:
Out[9]:
Out[9]:
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]:
Out[10]:
Out[10]:
Out[10]:
Out[10]:
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]:
In [12]:
let iter f t =
for k=1 to nnz t do
f t.indices.{k} t.values.{k}
done
Out[12]:
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]:
Out[13]:
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]:
Out[14]:
Test¶
In [44]:
to_assoc_list sparse_a;;
to_vec sparse_a = dense_a;;
Out[44]:
Out[44]:
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]:
Out[45]:
Out[45]:
Out[45]:
Out[45]:
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]:
Out[47]:
Out[47]:
Out[47]:
Out[47]:
Out[47]:
Out[47]:
Out[47]:
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]:
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]:
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]:
Out[21]:
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]:
Out[22]:
Out[22]:
Out[22]:
Out[22]:
Out[22]:
Out[22]:
Out[22]:
Out[22]:
Out[22]:
Out[22]:
Out[22]:
Out[22]:
Out[22]:
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]:
Test¶
In [24]:
dot sparse_a sparse_b = L.dot dense_a dense_b;;
Out[24]:
Tests¶
In [25]:
#require "alcotest";;
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]:
Out[26]:
Out[26]:
Out[26]:
Out[26]:
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]:
Out[27]:
Out[27]:
Out[27]:
In [28]:
(*)
Alcotest.run ~argv:[|"ignored"|] "Unit tests" [
"SpVec", test_case ();
];;
Printf.printf "%!";;
*)