Working on binary input

ocaml/Block.ml

@@ -62,19 +62,60 @@ let to_bytes b =
          Weight ;
          ...  ] ]  *)
   let l =
-    [ Property.to_bytes      b.property ;
+    [  Property.to_bytes      b.property ;
-      Sample.to_bytes        b.value  ;
+       Sample.to_bytes        b.value  ;
-      Weight.to_bytes        b.weight ;
+       Weight.to_bytes        b.weight ;
-      bytes_of_int           b.pid ;
+       bytes_of_int           b.pid ;
-      Block_id.to_bytes      b.block_id ;
+       Block_id.to_bytes      b.block_id ;
-      Compute_node.to_bytes  b.compute_node ]
+       Compute_node.to_bytes  b.compute_node ]
     |> List.map (fun x -> [ bytes_of_int (Bytes.length x) ;  x ] )
     |> List.concat
   in
   let result =
     Bytes.concat Bytes.empty (zero :: l)
   in
-  Bytes.set_int64_le result 8 (Int64.of_int (Bytes.length result));
+  Bytes.set_int64_le result 0 (Int64.of_int ((Bytes.length result) - 8));
+  result
+
+
+let read_bytes b =
+  (* Reads m, the first 8 bytes as an int64 containing the number of bytes to read.
+     Then, read the next m bytes and return a tuple containing the decoded data and the rest.
+  *)
+  let l = Bytes.length b in
+  let m =
+    Bytes.get_int64_le b 0
+    |> Int64.to_int
+  in
+  let nl = l-m-8 in
+  if nl > 0 then
+    (Bytes.sub b 8 m, Some (Bytes.sub b (8+m) nl))
+  else
+    (Bytes.sub b 8 m, None)
+
+
+
+let of_bytes b =
+  let b, _rest =
+    read_bytes b
+  in
+  let rec loop accu s =
+    match read_bytes s with
+    | data, None -> List.rev (data :: accu)
+    | data, (Some rest) -> loop (data :: accu) rest
+  in
+  let result =
+    match loop [] b with
+    | value :: weight :: property :: compute_node :: pid :: block_id :: [] ->
+          { property = Property.of_bytes property;
+            value    = Sample.of_bytes value;
+            weight   = Weight.of_bytes weight;
+            pid          = int_of_bytes pid;
+            block_id     = Block_id.of_bytes block_id;
+            compute_node = Compute_node.of_bytes compute_node;
+          }
+    | _ -> assert false
+  in
   result
 
 
@@ -87,11 +128,6 @@ let to_string b =
   (string_of_int          b.pid)
   (Block_id.to_int        b.block_id)
 
-(*
-let of_string s =
-  Bytes.of_string s
-  |> of_bytes
-*)
 
 let dir_name = lazy(
   let ezfio_filename =

ocaml/Qmcchem_dataserver.ml

@@ -228,13 +228,11 @@ let run ?(daemon=true) ezfio_filename =
       begin
         Ezfio.set_electrons_elec_coord_pool_size walk_num_tot ;
         let walkers_list =
-          Array.map Array.to_list walkers_array
+          Array.to_list walkers_array
-          |> Array.to_list
+          |> Array.concat
-          |> List.concat
+          |> Array.map float_of_string
-          |> List.rev_map float_of_string
-          |> List.rev
         in
-        Ezfio.set_electrons_elec_coord_pool (Ezfio.ezfio_array_of_list
+        Ezfio.set_electrons_elec_coord_pool (Ezfio.ezfio_array_of_array
           ~rank:3 ~dim:[| elec_num+1 ; 3 ; walk_num_tot |] ~data:walkers_list);
         let t0 =
           Unix.gettimeofday ()

ocaml/Random_variable.ml

@@ -194,6 +194,12 @@ let sum { property ; data } =
 (** Calculation of the average and error bar *)
 let ave_error { property ; data } =
 
+  (* sum:   \sum_k  x_k *. w_k
+     ansum: \sum_k  w_k
+     avsum: \sum_k  x_k *. w_k
+     avcu0: avsum / ansum
+     avsq:  \sum_k (1. -. (w_k /. ansum_k)) *. (x_k -. avcu0)^2 *. w_k)
+  *)
   let rec loop ~sum ~avsq ~ansum ~avsum ~n ?idx = function
     | [] ->
       begin
@@ -205,12 +211,8 @@ let ave_error { property ; data } =
       end
     | (x,w) :: tail ->
       begin
-        let avcu0 =
+        let avcu0 = avsum /. ansum in
-          avsum /. ansum
+        let xw = x *. w in
-        in
-        let xw =
-          x *. w
-        in
         let ansum, avsum, sum =
           ansum +. w ,
           avsum +. xw ,

ocaml/Sample.ml

@@ -43,7 +43,6 @@ let to_string = function
       Array.map string_of_float x
       |> Array.to_list
       |> String.concat " "
-      |> Printf.sprintf "%s"
 
 let to_bytes = function
   | One_dimensional  x -> Qptypes.bytes_of_float x
@@ -54,3 +53,12 @@ let to_bytes = function
         |> Bytes.set_int64_le b (i*8) ) x;
       b
 
+let of_bytes b =
+  match Bytes.length b with
+  | 8 -> let x = Qptypes.float_of_bytes b in
+         One_dimensional x
+  | l -> let len = l/8 in
+         Multidimensional ( Array.init len (fun i ->
+                               Bytes.get_int64_le b (i*8)
+                               |> Int64.float_of_bits ),
+                           len )

ocaml/Sample.mli

@@ -5,5 +5,6 @@ val of_float        : float -> t
 val of_float_array  : dim:int -> float array -> t
 val to_string       : t -> string
 val to_bytes        : t -> bytes
+val of_bytes        : bytes -> t
 val dimension       : t -> int
 

ocaml/qptypes_generator.ml

@@ -7,9 +7,12 @@ let global_replace x =
   |> Str.global_replace (Str.regexp "Int.to_bytes") "bytes_of_int"
   |> Str.global_replace (Str.regexp "Int64.to_bytes") "bytes_of_int64"
   |> Str.global_replace (Str.regexp "Float.to_bytes") "bytes_of_float"
+  |> Str.global_replace (Str.regexp "Float.of_bytes") "float_of_bytes"
   |> Str.global_replace (Str.regexp "Int.of_bytes") "int_of_bytes"
+  |> Str.global_replace (Str.regexp "Int64.of_bytes") "int64_of_bytes"
   |> Str.global_replace (Str.regexp "String.\\(to\\|of\\)_string") ""
   |> Str.global_replace (Str.regexp "String.to_bytes") "Bytes.of_string"
+  |> Str.global_replace (Str.regexp "String.of_bytes") "Bytes.to_string"
 
 let input_data = "
 * Positive_float : float
@@ -182,8 +185,22 @@ let bytes_of_int i =
   |> bytes_of_int64
 
 
+let int64_of_bytes b =
+  Bytes.get_int64_le b 0
+
+
+let int_of_bytes b =
+  int64_of_bytes b
+  |> Int64.to_int
+
+
+let float_of_bytes b =
+  int64_of_bytes b
+  |> Int64.float_of_bits
+
+
 let bytes_of_float f =
-  Int64.of_float f
+  Int64.bits_of_float f
   |> bytes_of_int64
 
 "
@@ -195,12 +212,14 @@ module %s : sig
   val of_%s : %s %s -> t
   val to_string : t -> string
   val to_bytes  : t -> bytes
+  val of_bytes  : bytes -> t
 end = struct
   type t = %s [@@deriving sexp]
   let to_%s x = x
   let of_%s %s x = ( %s x )
   let to_string x = %s.to_string x
   let to_bytes  x = %s.to_bytes  x
+  let of_bytes  b = %s.of_bytes  b
 end
 
 "
@@ -224,7 +243,7 @@ let parse_input input=
         and name = String_ext.strip name in
         let typ_cap = String.capitalize_ascii typ in
         let newstring = Printf.sprintf template name typ typ typ params_val typ typ
-          typ typ params ( String_ext.strip text ) typ_cap typ_cap
+          typ typ params ( String_ext.strip text ) typ_cap typ_cap typ_cap
         in
         List.rev (parse (newstring::result) tail )
   in
@@ -274,6 +293,10 @@ end = struct
 end
 "
 
+(*
+  val of_bytes  : bytes -> t
+  let of_bytes  x = %s.of_bytes  x
+*)
 
 let parse_input_ezfio input=
   let parse s =
@@ -320,7 +343,8 @@ let input_lines filename =
 let create_ezfio_handler () =
   let lines =
     input_lines "ezfio.ml"
-    |> List.mapi (fun i l -> if i > 417 then Some l else None)
+    (* /!\ Change when ezfio.ml changes *)
+    |> List.mapi (fun i l -> if i > 444 then Some l else None)
     |> List.filter (fun x -> x <> None)
     |> List.map (fun x ->
         match x with

src/MAIN/admc.py

@@ -0,0 +1,125 @@
+#!/usr/bin/env python3
+
+from mpi4py import MPI
+import sys
+import gzip
+import random
+import math
+import subprocess
+
+admc_exec = "/home/scemama/qmcchem/src/MAIN/admc"
+n_walk_per_proc = 10
+
+def start():
+    return subprocess.Popen(
+        [ admc_exec, sys.argv[1] ],
+        stdin=subprocess.PIPE,
+        stdout=subprocess.PIPE,
+        stderr=subprocess.PIPE)
+
+
+def read(process,len_walk):
+    line = process.stdout.readline().decode("utf-8").strip()
+    walk_num = int(line)
+    walkers = []
+    print(walk_num)
+    for k in range(walk_num):
+        w = []
+        for i in range(len_walk):
+            line = process.stdout.readline().decode("utf-8").strip()
+            w.append( line )
+        w = '\n'.join(w)
+        walkers.append(w)
+
+    _, E, W = process.stdout.readline().decode("utf-8").split()
+    return walkers, float(E), float(W)
+
+
+def write(process, message):
+    process.stdin.write(f"{message}\n".encode("utf-8"))
+    process.stdin.flush()
+
+
+def terminate(process):
+    process.stdin.close()
+    process.terminate()
+    process.wait(timeout=0.2)
+
+def print_energy(EnergyWeight, Energy2Weight, Weight, N):
+    e  = EnergyWeight / Weight
+    e2 = Energy2Weight / Weight
+    err = math.sqrt(abs(e*e - e2) / max(1,(N-1)) )
+    print("%f +/- %f"%(e, err))
+    return err
+
+def main():
+    try:
+      input_dir = sys.argv[1]
+    except:
+        print("syntax: argv[0] [FILE]")
+        sys.exit(-1)
+
+    # Pool of electron coordinates
+    with gzip.open(input_dir+"/electrons/elec_coord_pool.gz","r") as f:
+        data = f.read().decode("utf-8").split()
+
+    len_walk = int(data[1])*int(data[2])
+    icount = 0
+    buffer = []
+    walkers = []
+    for d in data[4:]:
+        buffer.append(d)
+        icount += 1
+        if (icount == len_walk):
+            walkers.append(buffer)
+            buffer = []
+            icount = 0
+
+    walkers = [ '\n'.join(x) for x in walkers ]
+    do_loop = True
+
+    EnergyWeight  = 0.
+    Energy2Weight = 0.
+    Weight        = 0.
+    NSamples      = 0.
+
+    # Start processes
+    proc = start()
+    while do_loop:
+
+        # Once every 1000, shuffle the list of walkers
+        if random.random() < 0.01:
+            print("SHUFFLE")
+            random.shuffle(walkers)
+
+        # Pick new walkers
+        new_coords = walkers[:n_walk_per_proc]
+        walkers = walkers[n_walk_per_proc:]
+
+        # Send new walkers to the process
+        write(proc, '\n'.join(new_coords))
+
+        # Fetch new walkers from the process
+        new_coords, e_new, w_new = read(proc, len_walk)
+        walkers += new_coords
+
+        # Print energy
+        ew = e_new * w_new
+        EnergyWeight  += ew
+        Energy2Weight += e_new * ew
+        Weight        += w_new
+        NSamples      += 1.
+        print (len(walkers))
+        err = print_energy(EnergyWeight, Energy2Weight, Weight, NSamples)
+
+        if err < 1.e-3:
+                do_loop = False
+
+    terminate(proc)
+    return
+
+
+
+
+if __name__ == "__main__":
+    main()