From 54f2bae5f6405ebdd41144848629ee5d7228a1ae Mon Sep 17 00:00:00 2001
From: Anthony Scemama <scemama@irsamc.ups-tlse.fr>
Date: Mon, 2 May 2016 21:51:09 +0200
Subject: [PATCH] Added operators on random variables

---
 ocaml/Qmcchem_result.ml                |  8 ++---
 ocaml/Random_variable.ml               | 45 ++++++++++++++++++++++++--
 src/PROPERTIES/properties_energy.irp.f |  2 +-
 src/SAMPLING/pdmc_step.irp.f           |  2 +-
 4 files changed, 48 insertions(+), 9 deletions(-)

diff --git a/ocaml/Qmcchem_result.ml b/ocaml/Qmcchem_result.ml
index c8cbc04..c1737d6 100644
--- a/ocaml/Qmcchem_result.ml
+++ b/ocaml/Qmcchem_result.ml
@@ -145,15 +145,15 @@ let display_summary ~range =
   in
   List.iter properties ~f:print_property ;
 
-(** TODO *)
+(*
   let open Random_variable in
-    let p = (of_raw_data ~range Property.E_loc_zv)
-         +! (of_raw_data ~range Property.E_loc)
+    let p =  (of_raw_data ~range Property.E_loc)
+          +! (of_raw_data ~range Property.E_loc_zv) 
     in
     Printf.printf "%20s : %s\n"
      ("E_loc_zv(+)")
      (Random_variable.to_string p);
-(** TODO *)
+ *)
 
   let cpu =
     Random_variable.of_raw_data ~range Property.Cpu
diff --git a/ocaml/Random_variable.ml b/ocaml/Random_variable.ml
index 94c46b2..c1ae8ad 100644
--- a/ocaml/Random_variable.ml
+++ b/ocaml/Random_variable.ml
@@ -457,20 +457,59 @@ let merge_per_compute_node_and_block_id =
        (Block_id.to_int block.Block.block_id) )
 
 
-(** Create a new random variable which is a sum of 2 *)
-let (+!) p1 p2 = 
+(** Create float, variable operators *)
+let one_variable_operator ~update_value p f = 
   merge
+  ~update_value
+  ~hashable:String.hashable
+  ~create_key:(fun block ->
+      Printf.sprintf "%s %10.10d %10.10d" 
+       (Compute_node.to_string block.Block.compute_node)
+       (Block_id.to_int block.Block.block_id) 
+       (Pid.to_int block.Block.pid) )
+  ~update_weight:(fun wc wb -> wc )
+  p
+
+let ( +@ ) p f = one_variable_operator p f
+  ~update_value: (fun  wc vc wb vb sw -> f +. (wc *. vc +. wb *. vb) /. sw) 
+
+let ( *@ ) p f = one_variable_operator p f 
+  ~update_value: (fun  wc vc wb vb sw -> f *. (wc *. vc +. wb *. vb) /. sw) 
+
+let ( -@ ) p f = one_variable_operator p f 
+  ~update_value: (fun  wc vc wb vb sw -> (wc *. vc +. wb *. vb) /. sw -. f)
+
+let ( /@ ) p f = one_variable_operator p f 
+  ~update_value: (fun  wc vc wb vb sw -> (wc *. vc +. wb *. vb) /. sw /. f)
+
+
+(** Create two variable operators *)
+let two_variable_operator ~update_value p1 p2 = 
+  merge
+  ~update_value
   ~hashable:String.hashable
   ~create_key:(fun block ->
       Printf.sprintf "%s %10.10d %10.10d" 
        (Compute_node.to_string block.Block.compute_node)
        (Block_id.to_int block.Block.block_id) 
        (Pid.to_int block.Block.pid) )
-  ~update_value: (fun wc vc wb vb sw -> (vc +. vb) )
   ~update_weight:(fun wc wb -> wc )
   { property = p1.property ;
     data = List.concat [ p1.data ; p2.data ] } 
 
+let ( +! ) = two_variable_operator  
+  ~update_value: (fun wc vc wb vb sw -> (vc +. vb) )
+
+let ( *! ) = two_variable_operator  
+  ~update_value: (fun wc vc wb vb sw -> (vc *. vb) )
+
+let ( -! ) = two_variable_operator  
+  ~update_value: (fun wc vc wb vb sw -> (vc -. vb) )
+
+let ( /! ) = two_variable_operator  
+  ~update_value: (fun wc vc wb vb sw -> (vc /. vb) )
+
+
 
 
 (** Merge two consecutive blocks *)
diff --git a/src/PROPERTIES/properties_energy.irp.f b/src/PROPERTIES/properties_energy.irp.f
index 727feec..1924d49 100644
--- a/src/PROPERTIES/properties_energy.irp.f
+++ b/src/PROPERTIES/properties_energy.irp.f
@@ -253,7 +253,7 @@ BEGIN_PROVIDER [ double precision, E_loc_zv ]
  BEGIN_DOC
  ! Zero-variance parameter on E_loc
  END_DOC
- E_loc_zv = 0.d0
+ E_loc_zv = E_loc
 END_PROVIDER
 
 
diff --git a/src/SAMPLING/pdmc_step.irp.f b/src/SAMPLING/pdmc_step.irp.f
index 8dd011b..9b3758b 100644
--- a/src/SAMPLING/pdmc_step.irp.f
+++ b/src/SAMPLING/pdmc_step.irp.f
@@ -174,7 +174,7 @@ END_SHELL
 
     if (calc_E_loc_zv) then
       if (dabs(pdmc_weight(i_walk)*pop_weight_mult) > 1.d-6) then
-        E_loc_zv = (E_ref-E_loc)/(pdmc_weight(i_walk)*pop_weight_mult)
+        E_loc_zv  = E_loc + (E_ref-E_loc)/(pdmc_weight(i_walk)*pop_weight_mult)
       else
         E_loc_zv = 0.d0
       endif