Commits

Anonymous committed b761506

Initial import.

Comments (0)

Files changed (3)

+
+
+(*** Type definitions ***)
+type expr = 
+     Factor of string
+   | Metab  of string
+   | Value  of int
+   | AND    of expr * expr
+   |  OR    of expr * expr
+   | NOT    of expr
+
+type rule = Rule of string * expr
+
+(*** Parser ***)
+let gram = Grammar.gcreate (Plexer.gmake ())
+
+let expr = Grammar.Entry.create gram "expr"
+let rule = Grammar.Entry.create gram "rule"
+
+let make_xor x y = ( AND(OR(x, y), NOT(AND(x, y))) )
+
+EXTEND
+   rule:
+      [ [ r = STRING; ":"; x = expr -> Rule(r, x) ] ];
+
+   expr:
+      [ 
+        [ x = expr; "AND"; y = expr -> AND(x, y)
+        | x = expr; "OR";  y = expr ->  OR(x, y)
+        | x = expr; "XOR"; y = expr -> make_xor x y ] 
+      | [ "NOT"; x = expr -> NOT(x) ]
+      | [ x = INT    -> Value(int_of_string x)
+        | x = STRING -> Metab(x) ]
+      | [ "("; x = expr; ")" -> x ]
+      ];
+END
+
+
+let get_expr r = match r with
+   Rule( n, e ) -> e
+
+let get_name r = match r with
+   Rule( n, e ) -> n
+
+
+let apply_to_exprs f rule_list = 
+   let apply_rule r = 
+      match r with
+         Rule( n, e ) -> Rule( n, f e )
+   in
+      List.map apply_rule rule_list
+
+let iter_apply_to_exprs f rule_list =
+   let n = List.length rule_list in
+   let old_array = Array.of_list rule_list in
+   let new_array = Array.create n (Rule("", Value 0)) in
+   begin
+      for i = 0 to n-1 do
+         new_array.(i) <- Rule( get_name old_array.(i), 
+                                f (get_expr old_array.(i)) );
+      done;
+      Array.to_list new_array
+   end
+
+let map_exprs f rule_list = 
+   let apply_rule r = 
+      match r with
+         Rule( n, e ) -> f e
+   in
+      List.map apply_rule rule_list
+
+(*** File manipulations ***)
+let file_to_list file = 
+   let channel = open_in file in
+      let rec file_to_list_aux l =
+         let s = 
+            try
+               input_line channel
+            with
+               End_of_file -> "" 
+         in
+            match s with
+                 "" -> List.rev l
+               | _  -> file_to_list_aux (s :: l)
+      in
+         file_to_list_aux []
+
+
+let parse_rule r = 
+   Grammar.Entry.parse rule (Stream.of_string r)
+
+let parse_rules rules = List.map parse_rule rules
+
+exception Parser_error of int
+
+let parse_file file = 
+   let raw_rules = file_to_list file in
+   let rules = ref [] in
+      begin
+         for i = 0 to List.length raw_rules - 1 do
+            try
+               rules := !rules @ [(parse_rule (List.nth raw_rules i))]
+            with
+               _ -> raise (Parser_error (i+1))
+         done;
+         !rules
+      end
+
+
+(*** Rule manipulation -- substitution ***)
+
+let sub_metab metab_names rule = 
+   let rec sub_aux expr = 
+      match expr with
+           Metab f -> if (List.mem f metab_names) then
+                          Metab f
+                       else
+                          Factor f
+         | AND(x, y) -> AND( sub_aux x, sub_aux y )
+         |  OR(x, y) ->  OR( sub_aux x, sub_aux y )
+         | NOT(x)    -> NOT( sub_aux x )
+         | x -> x
+   in
+      Rule( get_name rule, sub_aux (get_expr rule) )
+
+(* given a list of metabolite names, substitute change Factors to Metabs *)
+let sub_metabs metab_names rule_list = 
+   List.map (sub_metab metab_names) rule_list
+
+let find_uniques li = 
+   let rec unique_aux u l = 
+      match l with 
+           hd :: tl -> if (List.mem hd u) then
+                          unique_aux u tl
+                       else
+                          unique_aux (u @ [hd]) tl
+         | [] -> u
+   in
+      unique_aux [] li
+
+let get_rule_names rule_list = 
+   let aux x = (get_name x)
+   in
+      find_uniques (List.map aux rule_list)
+
+let get_metab_names rule_list = 
+   let rule_names = get_rule_names rule_list
+   in
+      let rec get_aux expr = 
+         match expr with
+              AND(x, y) -> (get_aux x) @ (get_aux y)
+            |  OR(x, y) -> (get_aux x) @ (get_aux y)
+            | NOT( x )  -> get_aux x
+            | Metab(x) -> if (List.mem x rule_names) then
+                              []
+                           else
+                              [x]
+            | _ -> []
+      in
+         let exprs = List.map get_expr rule_list
+         in
+            find_uniques (List.flatten (List.map get_aux exprs))
+
+let get_factor_names rule_list = 
+   let rule_names = get_rule_names rule_list
+   in
+      let rec get_aux expr = 
+         match expr with
+              AND(x, y) -> (get_aux x) @ (get_aux y)
+            |  OR(x, y) -> (get_aux x) @ (get_aux y)
+            | NOT( x )  -> get_aux x
+            | Metab(x) -> if (List.mem x rule_names) then
+                              [x]
+                           else
+                              []
+            | Factor(x) -> [x]
+            | _ -> []
+      in
+         let exprs = List.map get_expr rule_list
+         in
+            find_uniques (List.flatten (List.map get_aux exprs))
+
+let get_all_names rlist = 
+   (get_factor_names rlist) @ (get_metab_names rlist)
+
+
+
+
+(*** searching ***)
+let get_rule name rule_list = 
+   List.find (function x -> name = (get_name x) ) rule_list
+
+(*** Rule manipulations -- evaluation and simplification ***)
+
+let get_factor_rule factor rule_list = 
+   let aux x = (get_name x = factor)
+   in
+      List.find aux rule_list
+
+let is_solved_factor factor rule_list = 
+   let is_aux r =
+      match (get_expr r) with
+           Value x -> true
+         | _ -> false
+   in
+     is_aux (get_factor_rule factor rule_list)
+
+let reduce_rule solved_list rule = 
+   let rec is_solved x = 
+      List.mem_assoc x solved_list
+   and get_value x = 
+      Value (List.assoc x solved_list)
+   and reduce_expr r = 
+      match r with
+           Factor f -> if ( is_solved f ) then
+                          get_value f
+                       else
+                          Factor f
+         | Metab m  -> if ( is_solved m ) then
+                          get_value m
+                       else
+                          Metab m
+         | AND( Value x, Value y ) -> if ( x*y = 1 ) then
+                                         Value 1
+                                      else
+                                         Value 0
+         | AND( Value x, e )       -> if ( x = 0 ) then
+                                         Value 0
+                                      else 
+                                         reduce_expr e
+         | AND( e, Value x )       -> if ( x = 0 ) then
+                                         Value 0
+                                      else
+                                         reduce_expr e
+         | AND( e1, e2 ) -> AND( reduce_expr e1, reduce_expr e2 )
+
+         |  OR( Value x, Value y ) -> if ( x + y > 0 ) then
+                                         Value 1
+                                      else
+                                         Value 0
+         |  OR( Value x, e )       -> if ( x = 1 ) then
+                                         Value 1
+                                      else
+                                         reduce_expr e
+         |  OR( e, Value x )       -> if ( x = 1 ) then
+                                         Value 1
+                                      else
+                                         reduce_expr e
+         |  OR( e1, e2 ) -> OR( reduce_expr e1, reduce_expr e2 )
+
+         | NOT( Value x ) -> if ( x = 1 ) then
+                                Value 0
+                             else
+                                Value 1
+         | NOT( e ) -> NOT( reduce_expr e )
+
+         | x -> x
+
+      in
+         Rule( get_name rule, reduce_expr (get_expr rule) )
+
+
+let get_solved_factors rules = 
+   let is_solved r = 
+      match (get_expr r) with
+           Value x -> true
+         | _ -> false
+   and make_pair r = 
+      match (get_expr r) with
+        Value x -> (get_name r, x)
+      | _ -> ("", -1)
+   in
+      List.map make_pair (List.filter is_solved rules)
+
+let rec reduce_rules solved_list rules = 
+   let new_rules = List.map (reduce_rule solved_list) rules
+   in
+      if ( new_rules = rules ) then
+         rules
+      else
+         reduce_rules (solved_list @ (get_solved_factors new_rules)) new_rules
+
+let send_down_not rule = 
+   let rec send_aux expr = 
+      match expr with
+           NOT( AND(x, y) ) ->  OR( send_aux (NOT(x)), send_aux (NOT(y)) )
+         | NOT(  OR(x, y) ) -> AND( send_aux (NOT(x)), send_aux (NOT(y)) )
+         | NOT( NOT( x ) )  -> send_aux x
+         | NOT( Value 1 ) -> Value 0
+         | NOT( Value 0 ) -> Value 1
+
+         | AND(x, y) -> AND( send_aux x, send_aux y )
+         |  OR(x, y) ->  OR( send_aux x, send_aux y )
+         | NOT( x )  -> NOT( send_aux x )
+
+         | x -> x
+   in
+      Rule( (get_name rule), send_aux (get_expr rule) )
+
+let send_down_nots rule_list = 
+   List.map send_down_not rule_list
+
+
+(*** Importing procedure ***)
+
+(* Use this function to convert a file (named "file") to a list of rules *)
+let import_rules_from_file file = 
+   let rule_list = parse_file file in
+   let metab_names = get_metab_names rule_list in
+   let subed_rules = sub_metabs metab_names rule_list in
+   let reduced_rules = reduce_rules [] subed_rules in
+   let downed_rules = send_down_nots reduced_rules in
+      downed_rules
+
+
+(*** factoring routines ***)
+
+let rec has_or expr = 
+   match expr with
+         OR(x, y) -> true
+      | AND(x, y) -> (has_or x) or (has_or y)
+      | NOT( x )  -> has_or x
+      | _ -> false
+
+type direction = Left | Right
+
+let copy_by_or dir expr = 
+   let rec copy_aux e = 
+      match e with
+            OR(x, y) -> ( match dir with
+                             Left  -> x
+                           | Right -> y   )
+         | AND(x, y) -> if ( has_or x ) then
+                           AND( copy_aux x, y )
+                        else
+                           AND( x, copy_aux y )
+         | NOT( x ) -> NOT( copy_aux x )
+         | x -> x
+   in
+      copy_aux expr
+
+let copy_left  = copy_by_or Left
+let copy_right = copy_by_or Right 
+
+let split_once_by_or expr = 
+   if ( has_or expr ) then
+      [copy_left expr; copy_right expr]
+   else
+      [expr]
+
+           
+let rec swap_and_or expr = 
+   match expr with
+        AND(x, y) ->  OR( swap_and_or x, swap_and_or y )
+      |  OR(x, y) -> AND( swap_and_or x, swap_and_or y )
+      | NOT( x )  -> NOT( swap_and_or x )
+      | x -> x
+
+let split_expr_by_ors expr = 
+   let rec split_aux elist = 
+      let new_elist = List.flatten (List.map split_once_by_or elist)
+      in
+         if ( new_elist = elist ) then
+            elist
+         else
+            split_aux new_elist
+   in
+      split_aux [expr]
+
+let split_expr_by_ands expr = 
+   List.map swap_and_or (split_expr_by_ors (swap_and_or expr))
+
+
+type factored_rule = FactRule of string * expr list
+
+let factor_rule_by_ands r = 
+   FactRule ( get_name r, split_expr_by_ands (get_expr r) )
+
+let factor_rules_by_ands rule_list = 
+   List.map factor_rule_by_ands rule_list
+
+let factor_rule_by_ors r = 
+   FactRule ( get_name r, split_expr_by_ors (get_expr r) )
+
+let factor_rules_by_ors rule_list = 
+   List.map factor_rule_by_ors rule_list
+
+
+(*** Scale free analysis ***)
+
+let count_factors_in_rule rule = 
+   let rec count_aux e = 
+      match e with 
+           AND( x, y ) -> (count_aux x) @ (count_aux y)
+         |  OR( x, y ) -> (count_aux x) @ (count_aux y)
+         | NOT( x )    -> count_aux x
+         
+         | Factor( x ) -> [x]
+         | _ -> []
+   in
+      List.length (find_uniques (count_aux (get_expr rule)))
+
+let count_factors_in_rules rule_list = 
+   List.map count_factors_in_rule rule_list
+
+let count_metabs_in_rule rule = 
+   let rec count_aux e = 
+      match e with 
+           AND( x, y ) -> (count_aux x) @ (count_aux y)
+         |  OR( x, y ) -> (count_aux x) @ (count_aux y)
+         | NOT( x )    -> count_aux x
+         
+         | Metab( x ) -> [x]
+         | _ -> []
+   in
+      List.length (find_uniques (count_aux (get_expr rule)))
+
+let count_metabs_in_rules rule_list = 
+   List.map count_metabs_in_rule rule_list
+
+
+
+(*** correlation matrix routines ***)
+
+let not_mem  l a = not (List.mem a l)
+let all_true l   = List.for_all ( (=) true ) l
+
+let has_no_factors rule = not ((count_factors_in_rule rule) > 0)
+
+let has_been_seen name_list seen = 
+   List.mem true (List.map (function x -> List.mem x seen) name_list)
+
+let get_factors_in_rule rule = 
+   let rec get_aux e = 
+      match e with
+           AND(x, y) -> (get_aux x) @ (get_aux y)
+         |  OR(x, y) -> (get_aux x) @ (get_aux y)
+         | NOT( x )  -> get_aux x
+      
+         | Factor x -> [x]
+         | _ -> []
+   in
+      get_aux (get_expr rule)
+
+let is_embeddable gene_name rule_list = 
+   let rec is_embed_aux seen name = 
+      let rule = get_rule name rule_list in
+         if ( has_no_factors rule ) then
+            true
+         else
+            let factors = get_factors_in_rule rule in
+               if ( has_been_seen factors seen ) then
+                  false
+               else
+                  let recurse x = is_embed_aux (x :: seen) x in
+                     all_true (List.map (recurse) factors)
+   in
+      is_embed_aux [gene_name] gene_name
+
+
+(* this should only be used after the rule has been shown to be embeddable *)
+let embed_rule name rules = 
+   let rec embed_aux e = 
+      match e with
+           AND(x, y) -> AND( embed_aux x, embed_aux y )
+         |  OR(x, y) ->  OR( embed_aux x, embed_aux y )
+         | NOT( x )  -> NOT( embed_aux x )
+         
+         | Factor x -> embed_aux (get_expr (get_rule x rules))
+         | _ -> e
+   in
+      Rule( name, embed_aux (get_expr (get_rule name rules)) )
+
+let embed_rules names rules = 
+   let rec embed_aux rlist = 
+      match rlist with
+           hd :: tl -> if (List.mem (get_name hd) names) then
+                           (embed_rule (get_name hd) rules) :: (embed_aux tl)
+                       else
+                          hd :: (embed_aux tl)
+         | [] -> []
+   in
+      embed_aux rules
+
+let get_embeddable_rule_names rules = 
+   let rule_names = get_rule_names rules in
+   let pred = function name -> is_embeddable name rules in
+      List.filter pred rule_names
+
+let get_nonembeddable_rules rules = 
+   let rule_names = get_rule_names rules in
+   let embeddable_names = get_embeddable_rule_names rules in
+   let pred = function x -> not (has_been_seen [x] embeddable_names)
+   in
+      List.filter (pred) rule_names
+
+let partition_rules rules = 
+   let pred = function x -> is_embeddable (get_name x) rules
+   in
+      List.partition pred rules
+
+let get_embeddable_rules rules = 
+   fst (partition_rules rules)
+
+let get_nonembeddable_rules rules = 
+   snd (partition_rules rules)
+
+(* calculating c-matrix coefficients *)
+
+let rec s11 expr i = 
+   match expr with
+        Factor(x) -> 1
+      | Metab(x)  -> 1
+      | NOT( Metab(x)  ) -> if ( x = i ) then 0 else 1
+      | NOT( Factor(x) ) -> if ( x = i ) then 0 else 1
+      | NOT(_) -> 1  (* should not happen; send_down_nots *)
+
+      | Value(x)  -> if ( x = 1 ) then 1 else 0
+      
+      |  OR(x, y) -> (s11 x i)*((s11 y i) + (s01 y i))
+                      + (s11 y i)*(s01 x i)
+      | AND(x, y) -> (s11 x i)*(s11 y i)
+and s01 expr i = 
+   match expr with
+        Factor(x) -> if ( x = i ) then 0 else 1
+      | Metab(x)  -> if ( x = i ) then 0 else 1
+      | NOT( Metab(x)  ) -> 1
+      | NOT( Factor(x) ) -> 1
+      | NOT(_) -> 1  (* should not happen; send_down_nots *)
+      
+      | Value(x)  -> if ( x = 0 ) then 1 else 0
+      
+      |  OR(x, y) -> (s01 x i)*(s01 y i)
+      | AND(x, y) -> (s01 x i)*((s01 y i) + (s11 y i))
+                      + (s01 y i)*(s11 x i)
+and s10 expr i = 
+      match expr with
+        Factor(x) -> if ( x = i ) then 0 else 1
+      | Metab(x)  -> if ( x = i ) then 0 else 1
+      | NOT( Metab(x)  ) -> 1
+      | NOT( Factor(x) ) -> 1
+      | NOT(_) -> 1  (* should not happen; send_down_nots *)
+      
+      | Value(x)  -> if ( x = 1 ) then 1 else 0
+      
+      |  OR(x, y) -> (s10 x i)*((s10 y i) + (s00 y i))
+                      + (s10 y i)*(s00 x i)
+      | AND(x, y) -> (s10 x i)*(s10 y i)
+and s00 expr i = 
+   match expr with
+        Factor(x) -> 1
+      | Metab(x)  -> 1
+      | NOT( Metab(x)  ) -> if ( x = i ) then 0 else 1
+      | NOT( Factor(x) ) -> if ( x = i ) then 0 else 1
+      | NOT(_) -> 1  (* should not happen; send_down_nots *)
+      
+      | Value(x)  -> if ( x = 0 ) then 1 else 0
+      
+      |  OR(x, y) -> (s00 x i)*(s00 y i)
+      | AND(x, y) -> (s00 x i)*((s00 y i) + (s10 y i))
+                      + (s00 y i)*(s10 x i)
+
+let rec n_states expr = 
+   match expr with
+        NOT(x)    -> n_states x
+      |  OR(x, y) -> (n_states x)*(n_states y)
+      | AND(x, y) -> (n_states x)*(n_states y)
+      | _ -> 2
+
+let rec c_coeff expr i = 
+   (float_of_int ((s11 expr i) + (s00 expr i) - (s01 expr i) - (s10 expr i)))
+    /. (float_of_int (n_states expr))
+
+let get_coeff_cons r i = (i, (c_coeff r i))
+
+let get_coeffs mlist r = 
+   List.map (get_coeff_cons (get_expr r)) mlist
+
+let numbered_hash_of_list l = 
+   let hash = Hashtbl.create (List.length l) in
+   let arr  = Array.of_list l in
+      begin
+         for i = 0 to Array.length arr - 1 do
+            (* uses Matlab indexing *)
+            Hashtbl.add hash arr.(i) (i+1)
+         done;
+
+         hash
+      end
+
+let get_rule_hash rule_list = 
+   numbered_hash_of_list (get_rule_names rule_list)
+
+let get_atom_hash rule_list = 
+   numbered_hash_of_list (get_all_names rule_list)
+
+let build_c_matrix coeff_list_list rule_list =
+   let atoms = get_atom_hash rule_list in
+   let n_rows = List.length (List.hd coeff_list_list)
+   and n_cols = List.length coeff_list_list in
+   let matrix = Array.make_matrix n_rows n_cols 0.0 in
+      let add_column j coeff_list = 
+         let coeff_array = Array.of_list coeff_list in
+            for i = 0 to Array.length coeff_array - 1 do
+               let index = (Hashtbl.find atoms (fst coeff_array.(i))) - 1 in
+                  matrix.(index).(j) <- snd coeff_array.(i)
+            done
+      in
+         let column_array = Array.of_list coeff_list_list in
+            begin
+               for j = 0 to Array.length column_array - 1 do
+                  add_column j column_array.(j)
+               done;
+
+               matrix
+            end
+
+
+let make_c_matrix rules = 
+   let names = get_all_names rules in
+      build_c_matrix (List.map (get_coeffs names) (send_down_nots rules)) rules
+
+let export_c_matrix infile outfile embed = 
+   let raw_rules = import_rules_from_file infile in
+   let embeddable = get_embeddable_rule_names raw_rules in
+   let rules = if ( embed ) then
+                  embed_rules embeddable raw_rules
+               else
+                  raw_rules
+   in
+   let c_matrix = make_c_matrix rules in
+   let rule_names = get_rule_names rules in
+   let names = get_all_names rules in
+      begin
+         open_file outfile;
+         print (string_list_to_cell rule_names "column_names");
+         print (string_list_to_cell names "row_names");
+         if ( embed ) then
+            print (string_list_to_cell embeddable "embeddable_names");
+         print (make_array_of_float_matrix c_matrix "c_matrix");
+         close_file ()
+      end
+
+
+(*** Rule export -- algebraic format ***)
+
+let rules_to_alg rule_list = 
+   let atom_nums = get_atom_hash rule_list 
+   and rule_nums = get_rule_hash rule_list in
+
+   let rec expr_to_alg expr = 
+      match expr with
+           Metab  m -> indexed_array "x" (Hashtbl.find atom_nums m)
+         | Factor f -> indexed_array "x" (Hashtbl.find atom_nums f)
+         | Value  x -> string_of_int x
+      
+         | AND(x, y) -> "(" ^ (expr_to_alg x) ^ "*" ^ (expr_to_alg y) ^ ")"
+         |  OR(x, y) -> "(" ^ (expr_to_alg x) ^ " + " ^ (expr_to_alg y) 
+                        ^ " - " ^ (expr_to_alg (AND(x, y))) ^ ")"
+         | NOT( x )  -> "(1 - " ^ (expr_to_alg x) ^ ")"
+   and rule_to_alg rule = 
+      (indexed_array "g" (Hashtbl.find rule_nums (get_name rule)))
+         ^ " = " ^ (expr_to_alg (get_expr rule)) ^ ";\n"
+   in
+      List.map rule_to_alg rule_list
+
+
+let export_rules_to_alg infile outfile name_file = 
+   let rules = import_rules_from_file infile in
+   let rule_names = get_rule_names rules in
+   let names = get_all_names rules in
+      begin
+         open_file name_file;
+         print (string_list_to_cell rule_names "gene_names");
+         print (string_list_to_cell names "metab_names");
+         close_file ();
+
+         open_file (outfile ^ ".m");
+         print ("function g = " ^ outfile ^ "(x)\n\n");
+         let algs_array = Array.of_list (rules_to_alg rules) in
+            for i = 0 to (Array.length algs_array) - 1 do
+               print algs_array.(i)
+            done;
+         print "\n";
+         close_file ()
+      end
+
+
+(** correlation calculation exporting **)
+let rule_pair_to_alg rule_list = 
+   let atom_nums = get_atom_hash rule_list 
+   and rule_nums = get_rule_hash rule_list in
+
+   let rec expr_to_alg expr = 
+      match expr with
+           Metab  m -> indexed_array "x" (Hashtbl.find atom_nums m)
+         | Factor f -> indexed_array "x" (Hashtbl.find atom_nums f)
+         | Value  x -> string_of_int x
+      
+         | AND(x, y) -> "(" ^ (expr_to_alg x) ^ "*" ^ (expr_to_alg y) ^ ")"
+         |  OR(x, y) -> "(" ^ (expr_to_alg x) ^ " + " ^ (expr_to_alg y) 
+                        ^ " - " ^ (expr_to_alg (AND(x, y))) ^ ")"
+         | NOT( x )  -> "(1 - " ^ (expr_to_alg x) ^ ")"
+   in
+      List.map (function x -> expr_to_alg (get_expr x)) rule_list
+
+let export_all_embeddable_pairs pair_name_list pair_list outfile = 
+   let pairs = Array.of_list pair_list in
+   let pair_names = Array.of_list pair_name_list in
+      begin
+         open_file outfile;
+
+         for i = 1 to Array.length pairs do
+            begin
+               let n = List.length (get_all_names pairs.(i-1)) in
+                  print ("n = " ^ (string_of_int n) ^ ";\n");
+               print ("results{" ^ (string_of_int i) ^ ",1} = "
+                     ^ string_of_string (fst pair_names.(i-1)) ^ ";\n");
+               print ("results{" ^ (string_of_int i) ^ ",2} = "
+                     ^ string_of_string (snd pair_names.(i-1)) ^ ";\n");
+               let algs = rule_pair_to_alg pairs.(i-1) in
+                  begin
+                     print ("\nrule = @(x) [" ^ (List.nth algs 0) ^ ", " 
+                        ^ (List.nth algs 1) ^ "];\n");
+                     print ("results{" ^ (string_of_int i) 
+                        ^ ",3} = sim_pair(rule, n);\n");
+                  end;
+               if (i mod 1000 = 0) then
+                  begin
+                     print ("\n" ^ string_of_int i ^ "\n");
+                     print "\nsave results_sim.mat results;\n";
+                  end;
+            end
+         done;
+         close_file ();
+      end
+
+
+let make_pairs_of_rules infile outfile = 
+   let raw_rules = import_rules_from_file infile in
+   let embeddable = get_embeddable_rule_names raw_rules in
+   let rules = embed_rules embeddable raw_rules in
+
+   let rarr = Array.of_list rules in
+      let pairs = ref []
+      and name_pairs = ref [] in
+         begin
+            for i = 0 to Array.length rarr - 1 do
+               for j = 0 to i - 1 do
+                  begin
+                     pairs := [[rarr.(i); rarr.(j)]] @ !pairs;
+                     name_pairs := [(get_name rarr.(i), get_name rarr.(j))] 
+                                    @ !name_pairs
+                  end
+               done
+            done;
+         
+            export_all_embeddable_pairs !name_pairs !pairs outfile
+         end
+
+
+(*** generating random rules ***)
+let zero_one () = Random.int 2
+let true_false () = ( zero_one () = 0 )
+
+let is_even n = ( n mod 2 = 0 )
+let is_odd  n = not (is_even n)
+
+let half_floor n = n / 2
+let half_ceil  n = n - (half_floor n)
+
+let generate_random_metab n_total = 
+   Metab( "m" ^ (string_of_int ((Random.int n_total)+1)) )
+
+
+let rec generate_random_expr n_metabs n_total = 
+   match n_metabs with
+        1 -> if ( true_false () ) then
+               NOT( generate_random_metab n_total )
+             else
+                generate_random_metab n_total
+      | _ -> let x = (generate_random_expr (half_floor n_metabs) n_total)
+             and y = (generate_random_expr (half_ceil  n_metabs) n_total)
+             in
+               if ( true_false () ) then
+                  AND( x, y )
+               else
+                   OR( x, y )
+
+let rec generate_random_rule n_metabs n_total = 
+   Rule( "i" ^ (string_of_int (Random.int 999999999)), 
+         generate_random_expr n_metabs n_total )
+
+let rec generate_random_rule_list n_rules total_metabs metab_bound = 
+   let random_array = Array.create n_rules 0 in
+   begin
+      for i = 0 to Array.length random_array - 1 do
+         random_array.(i) <- Random.int metab_bound
+      done;
+
+      let random_list = Array.to_list random_array in
+         let f x = generate_random_rule x total_metabs 
+         in
+            List.map f random_list
+   end
+
+
+let rec expr_to_string expr = 
+   match expr with
+        Metab  x  -> dblqte ^ x ^ dblqte
+      | Factor x  -> dblqte ^ x ^ dblqte
+      | Value  x  -> string_of_int x
+
+      | NOT( x ) -> "NOT(" ^ (expr_to_string x) ^ ")"
+      | AND( x, y ) -> "(" ^ (expr_to_string x) ^ ") AND (" ^
+                         (expr_to_string y) ^ ")"
+      |  OR( x, y ) -> "(" ^ (expr_to_string x) ^ ") OR (" ^
+                         (expr_to_string y) ^ ")"
+
+let rule_to_string rule = 
+   match rule with
+      Rule( r, expr ) -> dblqte ^ r ^ dblqte ^ " : " ^ (expr_to_string expr)
+
+
+let create_rule_file n_rules metab_bound total_metabs outfile = 
+   for j = 1 to (n_rules / 10) do
+      open_append outfile;
+      for i = 1 to 10 do
+         Gc.full_major ();
+         let n_metabs = Random.int metab_bound in
+            show (rule_to_string (generate_random_rule n_metabs total_metabs))
+      done;
+      close_file ();
+   done
+
+(*** Final assembly ***)
+(* ! must have sent down nots before using these functions ! *)
+
+let not_string = "NOT__"
+
+let substitute_nots rules = 
+   let nots_subs = ref [] in
+   let add_sub name = 
+         nots_subs := !nots_subs @ [(name, not_string ^ name)]
+   in let make_sub e = 
+      match e with
+           Metab  x -> begin
+                          add_sub x;
+                          Metab  (not_string ^ x)
+                       end
+         | Factor x -> begin
+                          add_sub x;
+                          Factor (not_string ^ x)
+                       end
+         | Value  x -> if ( x = 0 ) then
+                          Value 1
+                       else
+                          Value 0
+         | _ -> e  (* should never happen; only called under NOTs *)
+   in
+      let rec sub_aux e = 
+         match e with
+              AND(x,y) -> AND( sub_aux x, sub_aux y )
+            |  OR(x,y) ->  OR( sub_aux x, sub_aux y )
+            | NOT( x ) -> make_sub x
+            | x -> x
+      in
+         let subed_rules = apply_to_exprs (sub_aux) rules
+         in
+            ( subed_rules, find_uniques (!nots_subs) )
+      
+
+
+
+type p_expr = 
+     P_Value  of int
+   | P_Factor of string
+   | P_Metab  of string
+
+   | Split_AND of p_expr * p_expr
+   | Subed_AND of p_expr * p_expr
+   | Split_OR  of p_expr * p_expr
+   | Subed_OR  of p_expr * p_expr
+
+   | P_NOT    of p_expr
+
+type p_rule = 
+   P_Rule of string * p_expr
+
+type split_or_subed = 
+     Split
+   | Subed
+
+
+let rec p_expr_to_expr p_expr = 
+   match p_expr with
+        P_Value(x)  -> Value x
+      | P_Metab(x)  -> Metab x
+      | P_Factor(x) -> Factor x
+
+      | Split_AND(x,y) -> AND( p_expr_to_expr x, p_expr_to_expr y )
+      | Subed_AND(x,y) -> AND( p_expr_to_expr x, p_expr_to_expr y )
+      | Split_OR (x,y) ->  OR( p_expr_to_expr x, p_expr_to_expr y )
+      | Subed_OR (x,y) ->  OR( p_expr_to_expr x, p_expr_to_expr y )
+
+      | P_NOT(x) -> NOT( p_expr_to_expr x )
+
+let p_rule_to_rule p_rule = 
+   match p_rule with
+      P_Rule(name, expr) -> Rule( name, p_expr_to_expr expr )
+
+let p_rules_to_rules p_rules = 
+   List.map (p_rule_to_rule) p_rules
+
+let get_p_name p_rule = 
+   get_name (p_rule_to_rule p_rule)
+
+let get_p_expr p_rule = 
+   match p_rule with
+      P_Rule(name, p_expr) -> p_expr
+
+let rec split_or_subed expr = 
+   let n_a = match expr with
+                  AND(x,y) -> n_expr x
+                |  OR(x,y) -> n_expr x
+                | _ -> 1
+   and n_b = match expr with
+                  AND(x,y) -> n_expr y
+                |  OR(x,y) -> n_expr y
+                | _ -> 1
+   in
+      if ((n_b > n_a) && (n_a >= 3)) || ((n_a > n_b) && (n_b >= 3)) then
+         Subed
+      else if (n_a = 3) && (n_b = 3) then
+         (* this is one of two arbitrary cases *)
+         Subed
+      else
+         (*Split*)
+         (* always sub *)
+         Subed
+   
+and n_gt expr = 
+   match expr with
+        Value  _ -> 1
+      | Factor _ -> 1
+      | Metab  _ -> 1
+      | NOT( x ) -> 1
+
+      | AND(x, y) -> if (( split_or_subed (AND(x,y)) ) = Subed) then
+                        1 + (n_gt x) + (n_gt y)
+                     else
+                        (n_gt x) * (n_gt y)
+
+      |  OR(x, y) -> if (( split_or_subed (OR(x,y)) ) = Subed ) then
+                        2 + (n_gt x) + (n_gt y)
+                     else
+                        (n_gt x) + (n_gt y)
+and n_lt expr = 
+   match expr with
+        Value  _ -> 1
+      | Factor _ -> 1
+      | Metab  _ -> 1
+      | NOT( x ) -> 1
+
+      | AND(x, y) -> if (( split_or_subed (AND(x,y)) ) = Subed ) then
+                        2 + (n_lt x) + (n_lt y)
+                     else
+                        (n_lt x) + (n_lt y)
+
+      |  OR(x, y) -> if (( split_or_subed (OR(x,y)) ) = Subed ) then
+                        1 + (n_lt x) + (n_lt y)
+                     else
+                        (n_lt x) * (n_lt y)
+
+and n_expr expr = (n_gt expr) + (n_lt expr)
+
+
+let rec prepare_expr expr = 
+   match expr with
+        Value x  -> P_Value x
+      | Metab x  -> P_Metab x
+      | Factor x -> P_Factor x
+      
+      | AND( x, y ) -> if ( split_or_subed expr = Subed ) then
+                          Subed_AND( prepare_expr x, prepare_expr y )
+                       else
+                          Split_AND( prepare_expr x, prepare_expr y )
+      |  OR( x, y ) -> if ( split_or_subed expr = Subed ) then
+                           Subed_OR( prepare_expr x, prepare_expr y )
+                       else
+                           Split_OR( prepare_expr x, prepare_expr y )
+      | NOT( x ) -> P_NOT( prepare_expr x )
+
+let rule_to_p_rule rule = 
+   match rule with
+      Rule( name, expr ) -> P_Rule( name, prepare_expr expr )
+
+let rules_to_p_rules rules = 
+   List.map (rule_to_p_rule) rules
+
+let sub_string = "sID__"
+
+let counter = ref 0
+let generate_unique_gene_id () = 
+   begin
+      counter := !counter + 1;
+      sub_string ^ (string_of_int !counter)
+   end
+
+let stored_rule_list = ref []
+
+let set_rule_list rules = 
+   stored_rule_list := rules
+
+let get_rule_list () = 
+   !stored_rule_list
+
+let append_rule rule = 
+   stored_rule_list := !stored_rule_list @ [rule]
+
+(* appends a new rule with expr; returns the rule name
+ * this function is used to make substitutions *)
+let rec add_expr p_expr = 
+   let add_sub pe = 
+      let uid = generate_unique_gene_id () in
+      let new_rule = P_Rule (uid, pe) in
+         begin
+            append_rule (make_substitution new_rule);
+            uid
+         end
+   in
+      match p_expr with
+           Subed_AND(x,y) -> P_Factor (add_sub p_expr)
+         | Subed_OR (x,y) -> P_Factor (add_sub p_expr)
+         | _ -> p_expr
+     
+
+and make_substitution rule = 
+   let rec sub_expr expr = 
+      match expr with
+(*  these cases should never be used; substitute only  *)
+(*           Split_AND (x,y) -> Split_AND ( sub_expr x, sub_expr y )   *)
+(*         | Split_OR  (x,y) -> Split_OR  ( sub_expr x, sub_expr y )   *)
+         
+           Subed_AND (x,y) -> Split_AND( add_expr x, add_expr y )
+         | Subed_OR  (x,y) -> Split_OR ( add_expr x, add_expr y )
+         | _ -> expr
+   in
+      P_Rule( get_p_name rule, sub_expr (get_p_expr rule) )
+                                 
+let make_substitutions rule_list = 
+   begin
+      set_rule_list [];
+      let subed_rules = List.map (make_substitution) rule_list
+      in
+         subed_rules @ (get_rule_list ())
+   end
+
+
+type in_eq_type = 
+     LEQ
+   | GEQ
+
+type in_eq = 
+     GT of (string * int) list * int
+   | LT of (string * int) list * int
+
+let is_gt_ineq ineq = 
+   match ineq with
+        GT(_, _) -> true
+      | _ -> false
+
+let is_lt_ineq ineq = 
+   not (is_gt_ineq ineq)
+
+let get_gt_ineqs ineq_list = 
+   List.filter (is_gt_ineq) ineq_list
+
+let get_lt_ineqs ineq_list = 
+   List.filter (is_lt_ineq) ineq_list
+
+(*
+let get_var ineq = 
+   match ineq with
+        GT( n, _, _ ) -> n
+      | LT( n, _, _ ) -> n
+let get_list ineq = 
+   match ineq with
+        GT( _, l, _ ) -> l
+      | LT( _, l, _ ) -> l
+let get_int ineq = 
+   match ineq with
+        GT( _, _, i ) -> i
+      | LT( _, _, i ) -> i
+*)
+
+(********************************************************)
+
+(*
+let combine_ineqs typ l1 l2 = 
+   let al1 = Array.of_list l1
+   and al2 = Array.of_list l2 in
+   let n_l1 = Array.length al1
+   and n_l2 = Array.length al2 in
+   let carr = Array.create (n_l1 * n_l2) (GT("", [], 0)) in
+   begin
+      for i = 0 to n_l1 - 1 do
+         for j = 0 to n_l2 - 1 do
+            carr.(i+n_l1*j) <- if (typ = GEQ) then
+                                  GT( "", 
+                                    (get_list al1.(i)) @ (get_list al2.(j)),
+                                    (get_int al1.(i)) + (get_int al2.(j)) - 1 )
+                               else
+                                  LT( "", 
+                                    (get_list al1.(i)) @ (get_list al2.(j)),
+                                    (get_int al1.(i)) + (get_int al2.(j)) - 1 )
+         done;
+      done;
+      Array.to_list carr
+   end
+*)
+
+(*******************************************************)
+
+(* this function is only used when true splits (everything not subed) *)
+(*
+let rec p_expr_to_ineqs_with_true_splits p_expr = 
+   match p_expr with
+        Split_AND(x,y) ->   (get_lt_ineqs (p_expr_to_ineqs x))
+                          @ (get_lt_ineqs (p_expr_to_ineqs y))
+                          @ (combine_ineqs GEQ 
+                                           (get_gt_ineqs (p_expr_to_ineqs x))
+                                           (get_gt_ineqs (p_expr_to_ineqs y)))
+      | Split_OR(x,y)  ->   (get_gt_ineqs (p_expr_to_ineqs x))
+                          @ (get_gt_ineqs (p_expr_to_ineqs y))
+                          @ (combine_ineqs LEQ
+                                           (get_lt_ineqs (p_expr_to_ineqs x))
+                                           (get_gt_ineqs (p_expr_to_ineqs y)))
+
+      | P_Metab  x -> [ GT("", [x], 0); LT("", [x], 0) ]
+      | P_Factor x -> [ GT("", [x], 0); LT("", [x], 0) ]
+      | P_Value  x -> [ GT("",  [], x); LT("",  [], x) ]
+      | _ -> []
+*)
+
+let p_expr_to_ineqs p_expr name = 
+   let name_of pe = match pe with
+        P_Metab  x -> x
+      | P_Factor x -> x
+      | _ -> "ERROR"
+   and atom_to_ineq atom_name = 
+      [ GT([ (name, 1); (atom_name, -1) ], 0);
+        LT([ (name, 1); (atom_name, -1) ], 0)  ]
+   in
+      match p_expr with
+           Split_AND(x,y) -> [ GT([ (name,     -4); 
+                                    (name_of x,  2); 
+                                    (name_of y,  2) ], -1);
+                               LT([ (name,      -4); 
+                                    (name_of x,  2); 
+                                    (name_of y,  2) ],  3) ] 
+         | Split_OR (x,y) -> [ GT([ (name,       3);
+                                    (name_of x, -1);
+                                    (name_of y, -1) ],  0);
+                               LT([ (name,       3);
+                                    (name_of x, -1);
+                                    (name_of y, -1) ],  2) ]
+         
+         | P_Factor x -> atom_to_ineq x
+         | P_Metab  x -> atom_to_ineq x
+
+(*
+let name_ineq name ineq = 
+   if (is_gt_ineq ineq) then
+      GT( name, (get_list ineq), (get_int ineq) )
+   else
+      LT( name, (get_list ineq), (get_int ineq) )
+*)
+
+let p_rule_to_ineqs p_rule = 
+   match p_rule with
+      P_Rule( name, p_expr ) -> p_expr_to_ineqs p_expr name
+
+let p_rules_to_ineqs p_rule_list = 
+   List.flatten (List.map (p_rule_to_ineqs) p_rule_list)
+
+let iter_p_rules_to_ineqs p_rule_list = 
+   let rule_array = Array.of_list p_rule_list in
+   let ineq_array = Array.create (Array.length rule_array) [] in
+      begin
+         for i = 0 to Array.length rule_array - 1 do
+            ineq_array.(i) <- p_rule_to_ineqs rule_array.(i)
+         done;
+         List.flatten (Array.to_list ineq_array)
+      end
+
+
+let name_list_to_index_list name_list name_hash = 
+   List.map (Hashtbl.find name_hash) name_list
+
+let add_inv index_list = 
+   List.map ( ( * ) (-1) ) index_list
+
+
+(* this function takes a pair of a list of rules 
+ * and a list of not substitutions; this is the 
+ * output of the substitute_nots function *)
+let get_var_names p_rules not_subs = 
+   let rules = p_rules_to_rules p_rules 
+   in
+      let split = List.split not_subs in
+         find_uniques (   (get_rule_names rules) 
+                        @ (get_metab_names rules) 
+                        @ (fst split) 
+                        @ (snd split) )
+
+(* rules should not have the nots replaced *)         
+let get_metab_indexes rules var_names = 
+   let metabs = (get_metab_names rules) in
+      let hsh = numbered_hash_of_list var_names in
+         List.map (Hashtbl.find hsh) metabs
+
+let get_not_eq_index_lists not_subs var_names = 
+   let hsh = numbered_hash_of_list var_names in
+   let f x = [Hashtbl.find hsh (fst x); Hashtbl.find hsh (snd x)]
+   in
+      List.map (f) not_subs
+
+let rec weave l1 l2 = 
+   match l1 with
+        hd :: tl -> [hd] @ [(List.hd l2)] @ (weave tl (List.tl l2))
+      | [] -> []
+
+let ineq_to_array name_hash ineq = 
+   let transform l coeff = 
+      let sl = List.split l
+      in
+         let coeffs = if ( coeff < 0 ) then
+                         add_inv (snd sl)
+                      else
+                         (snd sl)
+         in
+            weave (name_list_to_index_list (fst sl) name_hash) coeffs
+   in
+      match ineq with
+           LT( l, i ) ->      [i] @ (transform l ( 1))
+         | GT( l, i ) ->   [-1*i] @ (transform l (-1))
+
+
+let ineqs_to_arrays ineqs name_hash = 
+   List.map (ineq_to_array name_hash) ineqs
+
+
+(* finding rules with duplicates *)
+
+let rec get_strings e = 
+   match e with
+        AND(x,y) -> (get_strings x) @ (get_strings y)
+      |  OR(x,y) -> (get_strings x) @ (get_strings y)
+      | NOT( x ) -> get_strings x
+      | Metab  x -> [x]
+      | Factor x -> [x]
+      | Value  x -> [(string_of_int x)]
+
+let get_strings_of_rule r = 
+   get_strings (get_expr r)
+
+let rule_contains_duplicates r = 
+   let strs = get_strings_of_rule r
+   in
+      (List.length (find_uniques strs)) < (List.length strs)
+
+let get_rules_with_duplicates rule_list = 
+   List.filter (rule_contains_duplicates) rule_list
+
+
+(*** Final conversion routines ***)
+
+let convert_rules_to_ineqs raw_rules embed = 
+   let metab_names = get_metab_names raw_rules in
+   let rules = if ( embed ) then 
+                  let embeddable_rules = get_embeddable_rule_names raw_rules
+                  in
+                     get_embeddable_rules (embed_rules embeddable_rules
+                                                       raw_rules)
+               else
+                  raw_rules
+   in
+   let results = substitute_nots rules in
+   let not_subs = snd results in
+
+   let p_rules = make_substitutions (rules_to_p_rules (fst results)) in
+   let ineqs = iter_p_rules_to_ineqs p_rules in
+      (ineqs, (p_rules, not_subs))
+
+let export_rules infile embed = 
+   let rules = import_rules_from_file infile in
+   let results = convert_rules_to_ineqs rules embed
+   in
+   let ineqs    = fst results
+   and p_rules  = fst (snd results)
+   and not_subs = snd (snd results) in
+
+   let metab_names = get_metab_names rules in
+   let var_names = get_var_names p_rules not_subs in
+   let var_hash  = numbered_hash_of_list var_names in
+
+   let list_to_string l = delimiter_list (string_of_int) "\t" l
+   in
+      begin
+         open_file "rule_names.m";
+         show (string_list_to_cell var_names "cifba.var_names");
+         show (string_list_to_cell metab_names "cifba.metab_names");
+         show (int_list_to_array (get_metab_indexes rules var_names)
+                                 "cifba.metab_indexes");
+         close_file ();
+
+         open_file "rule_ineqs.txt";
+         List.iter (function x -> show (list_to_string x))
+                   (ineqs_to_arrays ineqs var_hash);
+         close_file ();
+
+         open_file "rule_not_eqs.txt";
+         List.iter (function x -> show (list_to_string x)) 
+                   (get_not_eq_index_lists not_subs var_names);
+         close_file ()
+      end
+
+(*** GPR conversion ***)
+let change_name_in_expr old_name new_name expr = 
+   let rec change e = 
+      match e with
+           Metab x  -> if ( x = old_name ) then 
+                          Metab new_name
+                       else
+                          Metab x
+         | Factor x -> if ( x = old_name ) then
+                          Factor new_name
+                       else
+                          Factor x
+         | AND(x,y) -> AND( change x, change y )
+         |  OR(x,y) ->  OR( change x, change y )
+         | NOT( x ) -> NOT( change x )
+        
+   in
+      change expr
+
+
+let change_name_in_rules rules old_name new_name = 
+   iter_apply_to_exprs (change_name_in_expr old_name new_name) rules
+
+let change_gpr_names rules = 
+   let names = get_metab_names rules in
+   let sub_array = Array.create (List.length names) "" in
+   let sub_list = 
+      begin
+         for i = 1 to (List.length names) do
+            sub_array.(i-1) <- "x(" ^ (string_of_int i) ^ ")"
+         done;
+         Array.to_list sub_array
+      end
+   in
+      let rs = ref rules in
+         begin
+            for i = 0 to (List.length names) - 1 do
+               rs := change_name_in_rules !rs 
+                                          (List.nth names i)
+                                          (List.nth sub_list i)
+            done;
+            !rs
+         end
+
+let rec expr_to_gpr_string expr = 
+   match expr with
+        Metab  x -> x
+      | Factor x -> x
+
+      | NOT( x ) -> "1 - (" ^ (expr_to_gpr_string x) ^ ")"
+      | AND( x, y ) -> "(" ^ (expr_to_gpr_string x) ^ ") & (" ^
+                         (expr_to_gpr_string y) ^ ")"
+      |  OR( x, y ) -> "(" ^ (expr_to_gpr_string x) ^ ") | (" ^
+                         (expr_to_gpr_string y) ^ ")" 
+      | _ -> ""
+
+let rules_to_gpr_strings rules = 
+   map_exprs (expr_to_gpr_string) rules
+
+
+let export_gpr infile = 
+   let rules = import_rules_from_file infile in
+   let gene_names = get_metab_names rules in
+   let rxn_names  = get_rule_names rules in
+   let rule_strings = rules_to_gpr_strings rules in
+   let subed_rule_strings = rules_to_gpr_strings (change_gpr_names rules) in
+
+   begin
+      open_file "gpr_names.m";
+      show (string_list_to_cell gene_names "gene_names");
+      show (string_list_to_cell rxn_names "rxn_names");
+      show (string_list_to_cell rule_strings "grRules");
+      show (string_list_to_cell subed_rule_strings "rules");
+      close_file ()
+   end
+
+(*** File IO methods ***)
+
+let channel = ref (open_out "temp.temp")
+
+(* returns a function to write strings to the channel *)
+let open_file filename = 
+   channel := open_out filename
+
+let open_append filename = 
+   channel := open_out_gen [Open_append] 777 filename
+
+let close_file () = 
+   close_out !channel
+
+let print str = output_string !channel str
+
+let show str = print (str ^ "\n")
+
+(*** export methods ***)
+
+let dblqte  = String.make 1 '\x22'
+let snglqte = String.make 1 '\x27'
+
+let bool_to_string b = 
+   match b with
+        true  -> "1"
+      | false -> "0"
+
+let string_of_string str = 
+   snglqte ^ str ^ snglqte
+
+let no_format x = x
+
+let delimiter_list string_of_t delim li = 
+   let rec delim_aux l = 
+      match l with
+           hd :: [] -> string_of_t hd
+         | hd :: tl -> (string_of_t hd) ^ delim ^ (delim_aux tl)
+         | [] -> ""
+   in
+      delim_aux li
+
+let comma_list string_of_t = 
+   delimiter_list string_of_t ", "
+let semicolon_list string_of_t = 
+   delimiter_list string_of_t "; "
+
+let make_array_of_list string_of_t = 
+   comma_list string_of_t
+
+let t_list_to_array t_list_to_string t_list name = 
+   name ^ " = [" ^ (t_list_to_string t_list) ^ "];\n"
+
+let t_list_to_cell t_list_to_string t_list name = 
+   name ^ " = {" ^ (t_list_to_string t_list) ^ "};\n"
+
+(* returns a string "[1, 2, 3]" *)
+let make_array_of_int_list    = make_array_of_list string_of_int
+let make_array_of_float_list  = make_array_of_list string_of_float
+let make_array_of_bool_list   = make_array_of_list bool_to_string
+let make_array_of_string_list = make_array_of_list string_of_string
+
+(* returns a string "name = [1, 2, 3];\n" *)
+let    int_list_to_array = t_list_to_array make_array_of_int_list
+let  float_list_to_array = t_list_to_array make_array_of_float_list
+let   bool_list_to_array = t_list_to_array make_array_of_bool_list
+let string_list_to_array = t_list_to_array make_array_of_string_list
+
+(* returns a string "name = {1, 2, 3};\n" *)
+let    int_list_to_cell = t_list_to_cell make_array_of_int_list
+let  float_list_to_cell = t_list_to_cell make_array_of_float_list
+let   bool_list_to_cell = t_list_to_cell make_array_of_bool_list
+let string_list_to_cell = t_list_to_cell make_array_of_string_list
+
+let make_array_of_float_matrix m name = 
+   let lists = Array.create (Array.length m) "" in
+   begin
+      for i = 0 to Array.length m - 1 do
+         lists.(i) <- comma_list (string_of_float) (Array.to_list m.(i))
+      done;
+
+      name ^ " = [" ^ (semicolon_list (no_format) (Array.to_list lists))
+         ^ "];\n"
+   end
+
+let indexed_array name index = 
+   name ^ "(" ^ (string_of_int index) ^ ")"
+
+#load "camlp4o.cma";;
+#load "pa_extend.cmo";;
+
+#use "MatlabIO.ml";;
+#use "Boolean.ml";;
+