Added set_mo_class ocaml script

ocaml/.gitignore

@@ -1,4 +1,5 @@
 ezfio.ml
+qptypes.ml
 *.byte
 *.native
 _build

ocaml/Makefile

@@ -9,16 +9,21 @@ $(info  -----------------------------------------------)
 $(error )
 endif
 
+
+EXECUTABLES=set_mo_class
+
 LIBS=
 PKGS=
 OCAMLCFLAGS=-g
 OCAMLBUILD=ocamlbuild -cflags $(OCAMLCFLAGS) -lflags -g
-MLFILES=$(wildcard *.ml) ezfio.ml
+MLFILES=$(wildcard *.ml) ezfio.ml qptypes.ml
 MLIFILES=$(wildcard *.mli) 
 ALL_TESTS=$(patsubst %.ml,%.byte,$(wildcard test_*.ml))
 
 
-default: $(ALL_TESTS)
+default: $(ALL_TESTS) executables
+
+executables: $(patsubst %, %.native, $(EXECUTABLES))
 
 %.inferred.mli: $(MLFILES)
 	$(OCAMLBUILD) $*.inferred.mli -cflags -i -use-ocamlfind  $(PKGS)
@@ -32,6 +37,13 @@ default: $(ALL_TESTS)
 ezfio.ml: ${QPACKAGE_ROOT}/EZFIO/Ocaml/ezfio.ml
 	cp ${QPACKAGE_ROOT}/EZFIO/Ocaml/ezfio.ml .
 
+qptypes_generator.byte: qptypes_generator.ml
+	$(OCAMLBUILD) qptypes_generator.byte  -use-ocamlfind 
+
+qptypes.ml: qptypes_generator.byte
+	./qptypes_generator.byte > qptypes.ml
+	rm qptypes_generator.byte
+
 ${QPACKAGE_ROOT}/EZFIO/Ocaml/ezfio.ml:
 	$(MAKE) -C ${QPACKAGE_ROOT}/src ezfio
 

ocaml/bit.ml

@@ -0,0 +1,44 @@
+
+(*
+Type for bits 
+==============
+
+Zero | One
+
+*)
+
+type bit      = 
+   | One 
+   | Zero 
+
+let to_string = function
+  | Zero -> "0"
+  | One  -> "1"
+;;
+
+let and_operator a b = 
+  match a, b with
+  | Zero, _ -> Zero
+  | _, Zero -> Zero
+  | _, _ -> One
+;;
+
+let or_operator a b = 
+  match a, b with
+  | One, _ -> One
+  | _, One -> One
+  | _, _ -> Zero
+;;
+
+let xor_operator a b =
+  match a, b with
+  | One, Zero -> One
+  | Zero, One -> One
+  | _, _ -> Zero
+;;
+
+let not_operator = function
+  | One -> Zero
+  | Zero -> One
+;;
+

ocaml/bitlist.ml

@@ -0,0 +1,148 @@
+open Qptypes;;
+
+(*
+Type for bits strings
+=====================
+
+list of Bits
+*)
+
+type bit_list = Bit.bit list
+
+(* String representation *)
+let to_string b = 
+  let rec do_work accu = function
+    | [] -> accu
+    | head :: tail ->
+      let new_accu = (Bit.to_string head) ^ accu
+      in do_work new_accu tail
+  in
+  do_work "" b
+;;
+
+
+(* Create a bit list from an int64 *)
+let of_int64 i = 
+  let rec do_work = function
+  | 0L -> [ Bit.Zero ]
+  | 1L -> [ Bit.One ]
+  | i -> let b =
+    match (Int64.logand i 1L ) with
+    | 0L -> Bit.Zero
+    | 1L -> Bit.One 
+    | _ -> raise (Failure "i land 1 not in (0,1)")  
+    in b:: ( do_work (Int64.shift_right_logical i 1) )
+  in
+  let adjust_length result = 
+    let rec do_work accu = function
+    | 64 -> accu 
+    | i when i>64 -> raise (Failure "Error in of_int64 > 64")
+    | i when i<0 -> raise (Failure "Error in of_int64 < 0")
+    | i -> do_work (accu@[Bit.Zero])  (i+1)
+    in
+    do_work result (List.length result)
+  in
+  adjust_length (do_work i)
+;;
+
+(* Create an int64 from a bit list *)
+let to_int64 l =
+  assert ( (List.length l) <= 64) ;
+  let rec do_work accu = function
+    | [] -> accu
+    | Bit.Zero::tail -> do_work Int64.(shift_left accu 1) tail
+    | Bit.One::tail  -> do_work Int64.(logor one (shift_left accu 1)) tail
+  in do_work Int64.zero (List.rev l)
+;;
+
+(* Create a bit list from a list of int64 *)
+let of_int64_list l = 
+  let list_of_lists = List.map of_int64 l in
+  let result = List.rev list_of_lists in
+  List.flatten result
+;;
+
+(* Create an int64 list from a bit list *)
+let to_int64_list l =
+  let rec do_work accu buf counter = function
+    | [] -> 
+        begin
+          match buf with
+          | [] -> accu
+          | _  -> (List.rev buf)::accu
+        end
+    | i::tail -> 
+        if (counter < 64) then
+          do_work accu (i::buf) (counter+1) tail
+        else
+          do_work ( (List.rev (i::buf))::accu) [] 1 tail
+  in
+  let l = do_work [] [] 1 l
+  in
+  List.map to_int64 l
+;;
+
+(* Create a bit list from a list of MO indices *)
+let of_mo_number_list n_int l = 
+  let n_int = N_int_number.to_int n_int in
+  let length = n_int*64 in
+  let a = Array.make length (Bit.Zero) in
+  List.iter (fun i-> a.((MO_number.to_int i)-1) <- Bit.One) l;
+  Array.to_list a
+;;
+
+
+
+(* logical operations on bit_list *)
+let logical_operator2 op a b =
+  let rec do_work_binary result a b = 
+  match a, b with
+  | [], [] -> result
+  | [], _  | _ , [] -> raise (Failure "Lists should have same length")
+  | (ha::ta), (hb::tb) -> 
+    let newbit = op ha hb
+    in  do_work_binary (newbit::result) ta tb
+  in
+  List.rev (do_work_binary [] a b)
+;;
+
+let logical_operator1 op b =
+  let rec do_work_unary result b = 
+  match b with
+  | [] -> result
+  | (hb::tb) -> 
+    let newbit = op hb
+    in  do_work_unary (newbit::result) tb
+  in
+  List.rev (do_work_unary [] b)
+;;
+
+let and_operator a b = logical_operator2 Bit.and_operator a b;;
+let xor_operator a b = logical_operator2 Bit.xor_operator a b;;
+let  or_operator a b = logical_operator2  Bit.or_operator a b;;
+let not_operator   b = logical_operator1 Bit.not_operator   b ;;
+
+let test_module () = 
+  let test = of_int64_list ([-1231L;255L]) in
+  print_string (to_string test);
+  print_newline ();
+  print_string (string_of_int (String.length (to_string test)));
+  print_newline ();
+  print_string ( Bit.to_string Bit.One );
+
+  let a = of_int64_list ([-1L;0L]) 
+  and b = of_int64_list ([128L;127L])
+  in begin
+   print_newline ();
+   print_newline ();
+   print_string (to_string a);
+   print_newline ();
+   print_string (to_string b);
+   print_newline ();
+   print_string (to_string (and_operator a b));
+   print_newline ();
+   print_string (to_string (or_operator a b));
+   print_newline ();
+   print_string (to_string (xor_operator a b));
+  end
+;;

ocaml/qptypes.ml

@@ -1,105 +0,0 @@
-module Positive_float : sig
-  type t
-  val to_float : t -> float
-  val of_float : float -> t
-end = struct
-  type t = float
-  let to_float x = x
-  let of_float x = ( assert (x >= 0.) ; x )
-end
-
-
-module Strictly_positive_float : sig
-  type t
-  val to_float : t ->float 
-  val of_float : float -> t
-end = struct
-  type t =float 
-  let to_float x = x
-  let of_float x = ( assert (x > 0.) ; x )
-end
-
-
-module Positive_int : sig
-  type t
-  val to_int : t -> int
-  val of_int : int -> t
-end = struct
-  type t = int
-  let to_int x = x
-  let of_int x = ( assert (x >= 0) ; x )
-end
-
-
-module Strictly_positive_int : sig
-  type t
-  val to_int : t -> int
-  val of_int : int -> t
-end = struct
-  type t = int
-  let to_int x = x
-  let of_int x = ( assert (x > 0) ; x )
-end
-
-
-module Non_empty_string : sig
-  type t
-  val to_string : t -> string
-  val of_string : string -> t
-end = struct
-  type t = string
-  let to_string x = x
-  let of_string x = ( assert (x <> "") ; x )
-end
-
-(*
-module MO_number : sig
-  type t
-  val to_int : t -> int
-  val of_int : int -> t
-end = struct
-  type t = int
-  let to_int x = x
-  let of_int x = ( assert (x > 0) ;
-  if (Ezfio.has_mo_basis_mo_tot_num ()) then
-    assert (x <= (Ezfio.get_mo_basis_mo_tot_num ())); x )
-end
-
-
-module AO_number : sig
-  type t
-  val to_int : t -> int
-  val of_int : int -> t
-end = struct
-  type t = int
-  let to_int x = x
-  let of_int x = ( assert (x > 0) ;
-  if (Ezfio.has_ao_basis_ao_num ()) then
-    assert (x <= (Ezfio.get_ao_basis_ao_num ())); x )
-end
-
-
-module N_int_number : sig
-  type t
-  val to_int : t -> int
-  val of_int : int -> t
-end = struct
-  type t = int
-  let to_int x = x
-  let of_int x = ( assert (x > 0) ;
-  if (Ezfio.has_determinants_n_int ()) then
-    assert (x == (Ezfio.get_determinants_n_int ())); x )
-end
-
-module Det_number : sig
-  type t
-  val to_int : t -> int
-  val of_int : int -> t
-end = struct
-  type t = int
-  let to_int x = x
-  let of_int x = ( assert (x > 0) ;
-  if (Ezfio.has_determinants_det_num ()) then
-    assert (x <= (Ezfio.get_determinants_det_num ())); x )
-end
-*)

ocaml/qptypes_generator.ml

@@ -0,0 +1,92 @@
+open Core.Std;;
+
+let input_data = "
+* Positive_float : float  
+  assert (x >= 0.) ; 
+
+* Strictly_positive_float : float  
+  assert (x > 0.) ; 
+
+* Negative_float : float  
+  assert (x <= 0.) ; 
+
+* Strictly_negative_float : float  
+  assert (x < 0.) ; 
+
+* Positive_int : int  
+  assert (x >= 0) ; 
+
+* Strictly_positive_int : int  
+  assert (x > 0) ; 
+
+* Negative_int : int  
+  assert (x <= 0) ; 
+
+* Strictly_negative_int : int  
+  assert (x < 0) ; 
+
+* Non_empty_string : string
+  assert (x <> \"\") ;
+
+* MO_number : int  
+  assert (x > 0) ; 
+  if (Ezfio.has_mo_basis_mo_tot_num ()) then
+    assert (x <= (Ezfio.get_mo_basis_mo_tot_num ()));
+
+* AO_number : int  
+  assert (x > 0) ; 
+  if (Ezfio.has_ao_basis_ao_num ()) then
+    assert (x <= (Ezfio.get_ao_basis_ao_num ()));
+
+* N_int_number : int 
+  assert (x > 0) ; 
+  if (Ezfio.has_determinants_n_int ()) then
+    assert (x == (Ezfio.get_determinants_n_int ()));
+
+* Det_number : int 
+  assert (x > 0) ; 
+  if (Ezfio.has_determinants_det_num ()) then
+    assert (x <= (Ezfio.get_determinants_det_num ()));
+"
+;;
+
+
+let template = format_of_string "
+module %s : sig
+  type t
+  val to_%s : t -> %s
+  val of_%s : %s -> t
+end = struct
+  type t = %s
+  let to_%s x = x
+  let of_%s x = ( %s x )
+end
+
+"
+;;
+
+let parse_input input=
+  let rec parse result = function
+    | [] -> result
+    | ( "" , ""   )::tail -> parse result tail
+    | ( t  , text )::tail -> 
+        let  name , typ  = String.lsplit2_exn ~on:':' t
+        in
+        let typ  = String.strip typ
+        and name = String.strip name
+        in
+        let newstring = Printf.sprintf template name typ typ typ typ typ typ typ
+        ( String.strip text )
+        in
+        List.rev (parse (newstring::result) tail )
+  in
+     String.split ~on:'*' input
+  |> List.map ~f:(String.lsplit2_exn ~on:'\n') 
+  |> parse []
+  |> String.concat 
+  |> print_string
+;;
+
+let () = parse_input input_data;;
+
+

ocaml/qputils.ml

@@ -0,0 +1 @@
+let (/) (a:string) (b:string) = a^"/"^b;;

ocaml/range.ml

@@ -0,0 +1,72 @@
+open Core.Std;;
+
+(* A range is a string of the type:
+ * 
+ * "[36-53,72-107,126-131]"
+ *
+ * that should represent the list of integers
+ * [ 37 ; 37 ; 38 ; ... ; 52 ; 53 ; 72 ; 73 ; ... ; 106 ; 107 ; 126 ; 127 ; ...
+ * ; 130 ; 131 ]
+*)
+
+
+type t = int list ;;
+
+let expand_range r =
+  match String.lsplit2 ~on:'-' r with
+  | Some (s, f) ->
+      begin
+        let start = Int.of_string s
+        and finish =  Int.of_string f
+        in
+        assert (start <= finish) ;
+        let rec do_work = function
+          | i when i=finish -> [ i ]
+          | i     -> i::(do_work (i+1))
+        in do_work start
+      end
+  | None -> 
+      begin
+        match r with
+          | "" -> []
+          | _  -> [Int.of_string r] 
+      end
+;;
+
+let of_string s =
+  assert (s.[0] = '[') ;
+  assert (s.[(String.length s)-1] = ']') ;
+  let s = String.sub s 1 ((String.length s) - 2) in
+  let l = String.split ~on:',' s in
+  let l = List.map ~f:expand_range l in
+  List.concat l |> List.dedup ~compare:Int.compare |> List.sort ~cmp:Int.compare
+;;
+
+let to_string l =
+  let rec do_work buf symbol = function
+    | [] -> buf
+    | a::([] as t) -> 
+          do_work (buf^symbol^(Int.to_string a)) "" t
+    | a::(b::q as t) ->
+        if (b-a = 1) then
+          do_work buf "-" t
+        else
+          do_work (buf^symbol^(Int.to_string a)^","^(Int.to_string b)) "" t
+  in
+  let result = 
+    match l with
+    | [] ->
+        "[]"
+    | h::t  ->
+        do_work ("["^(Int.to_string h)) "" l in
+  (String.sub result 0 ((String.length result)))^"]"
+;;
+
+let test_module () =
+  let s = "[72-107,36-53,126-131]" in
+  let l = of_string s in
+  print_string s ; print_newline () ;
+  List.iter ~f:(fun x -> Printf.printf "%d, " x) l ; print_newline () ;
+  to_string l |> print_string ;  print_newline () ;
+;;
+

ocaml/set_mo_class.ml

@@ -0,0 +1,170 @@
+open Qputils;;
+open Qptypes;;
+open Core.Std;;
+
+(*
+ * Command-line arguments
+ * ----------------------
+ *)
+
+let build_mask from upto n_int =
+  let from  = MO_number.to_int from
+  and upto  = MO_number.to_int upto
+  and n_int = N_int_number.to_int n_int
+  in
+  let rec build_mask bit = function
+    | 0 -> []
+    | i -> 
+        if ( i = upto ) then
+          Bit.One::(build_mask Bit.One (i-1))
+        else if ( i = from ) then
+          Bit.One::(build_mask Bit.Zero (i-1))
+        else
+          bit::(build_mask bit (i-1))
+  in
+  let starting_bit = 
+    if ( (upto >= n_int*64) || (upto < 0) ) then Bit.One
+    else Bit.Zero
+  in
+  build_mask starting_bit (n_int*64)
+  |> List.rev
+;;
+
+let mo_tot_num = ref 0;;
+let n_int      = ref (N_int_number.of_int 1);;
+
+let apply_mask mask = 
+  let full_mask = build_mask (MO_number.of_int 1) (MO_number.of_int !mo_tot_num) !n_int
+  in
+  let newmask = Bitlist.and_operator full_mask mask
+  in
+  (* TODO *)
+  newmask |> Bitlist.to_string |> print_endline;
+  string_of_int !mo_tot_num |> print_endline;
+;;
+
+
+
+let failure s = 
+  raise (Failure s)
+;;
+
+let run ?active ?(inactive="[]") ezfio_filename =
+
+  Ezfio.set_file ezfio_filename ;
+  if not (Ezfio.has_mo_basis_mo_tot_num ()) then
+    failure "mo_basis/mo_tot_num not found" ;
+
+  mo_tot_num := Ezfio.get_mo_basis_mo_tot_num () ;
+  n_int := N_int_number.of_int (Ezfio.get_determinants_n_int ()) ;
+
+  let inactive_mask = Range.of_string inactive
+    |> List.map ~f:MO_number.of_int
+    |> Bitlist.of_mo_number_list !n_int
+  and active_mask = 
+    let s = 
+    match active with
+    | Some range -> Range.of_string range
+    | None       -> Range.of_string ("[1-"^(Int.to_string !mo_tot_num)^"]")
+    in
+      List.map ~f:MO_number.of_int s
+      |> Bitlist.of_mo_number_list !n_int
+  in
+  let mask = 
+    Bitlist.not_operator inactive_mask
+    |> Bitlist.and_operator active_mask
+  in apply_mask mask
+;;
+
+let ezfio_file =
+  let failure filename = 
+        eprintf "'%s' is not an EZFIO file.\n%!" filename;
+        exit 1
+  in
+  Command.Spec.Arg_type.create
+  (fun filename ->
+    match Sys.is_directory filename with
+    | `Yes -> 
+        begin
+          match Sys.is_file (filename / ".version") with
+          | `Yes -> filename
+          | _ -> failure filename
+        end
+    | _ -> failure filename
+  )
+;;
+
+let default range =
+  let failure filename = 
+        eprintf "'%s' is not a regular file.\n%!" filename;
+        exit 1
+  in
+  Command.Spec.Arg_type.create
+  (fun filename ->
+    match Sys.is_directory filename with
+    | `Yes -> 
+        begin
+          match Sys.is_file (filename / ".version") with
+          | `Yes -> filename
+          | _ -> failure filename
+        end
+    | _ -> failure filename
+  )
+;;
+
+let spec =
+  let open Command.Spec in
+  empty 
+  +> flag "inactive"  (optional string) ~doc:"range Range of inactive orbitals"
+  +> flag "active"  (optional string)   ~doc:"range Range of active orbitals"
+  +> anon ("ezfio_filename" %: ezfio_file)
+;;
+
+let command = 
+    Command.basic 
+    ~summary: "Set the active/inactive orbitals in an EZFIO directory"
+    ~readme:(fun () ->
+      "The range of MOs has the form : \"[36-53,72-107,126-131]\"
+        ")
+    spec
+    (fun inactive active ezfio_filename () -> run ?inactive
+    ?active ezfio_filename )
+;;
+
+let () =
+    Command.run command
+
+
+(*
+let test_module () = 
+  let { Ezfio.rank ; Ezfio.dim ; Ezfio.data } = Ezfio.get_bitmasks_generators () in
+  let test = 
+    Ezfio.flattened_ezfio_data data 
+    |> Array.to_list 
+    |> List.map Int64.of_int 
+    |> Bitlist.of_int64_list 
+  in
+  print_string (Bitlist.to_string test);
+  print_newline ();
+  print_string (string_of_int (String.length (Bitlist.to_string test)));
+  print_newline ();
+
+  let a = Bitlist.of_int64_list ([-1L;0L]) 
+  and b = Bitlist.of_int64_list ([128L;127L])
+  in begin
+    print_newline ();
+    print_newline ();
+    Bitlist.to_string a |> print_string;
+    print_newline ();
+    Bitlist.to_string b |> print_string;
+    print_newline ();
+    Bitlist.and_operator a b |> Bitlist.to_string |> print_string;
+    print_newline ();
+    Bitlist.or_operator a b  |> Bitlist.to_string |> print_string;
+    print_newline ();
+    Bitlist.xor_operator a b |> Bitlist.to_string |> print_string;
+  end
+;;
+*)
+
+(*test_module ();;*)