c

20 files changed, 391 insertions, 129 deletions

diff --git a/algorithms/dune b/algorithms/dune
index a98c160..90f35b8 100644
--- a/algorithms/dune
+++ b/algorithms/dune
@@ -1,6 +1,6 @@
 (library
  (inline_tests)
  (name algorithms_lib)
- (libraries core owl genetics_lib utils_lib structures_lib)
+ (libraries core owl progress genetics_lib utils_lib structures_lib)
  (preprocess
   (pps ppx_inline_test ppx_jane)))
diff --git a/algorithms/pbil.ml b/algorithms/pbil.ml
index 63d8ee8..e4ea6af 100644
--- a/algorithms/pbil.ml
+++ b/algorithms/pbil.ml
@@ -4,29 +4,42 @@ open Utils_lib
 
 module Input = struct
   type t =
-    { gene_length : int
+    { chromosome_length : int
     ; population_size : int
     ; learning_coef : float
     ; mutation_p : float
     ; mutation_disturbance_coef : float
-    ; evaluate : Individual.t -> float
+    ; evaluate : Individual.t -> int
     ; termination_condition : Population.t -> iteration:int -> bool
+    ; gather_verbose_info : bool
     }
   [@@deriving sexp_of]
 end
 
 module Output = struct
+  module Iteration = struct
+    type t =
+      { distribution : float array
+      ; best_individual : Individual.t
+      ; mean_individual : Individual.t
+      ; worst_individual : Individual.t
+      }
+    [@@deriving sexp]
+  end
+
   type t =
     { distribution : float array
     ; population : Population.t
+    ; iterations_info : Iteration.t list option [@sexp.option]
     }
-  [@@deriving sexp_of]
+  [@@deriving sexp]
 end
 
 let adjust_distribution distribution individual learning_coef =
-  Gene.iteri (Individual.gene individual) ~f:(fun i value ->
+  Chromosome.iteri (Individual.chromosome individual) ~f:(fun i value ->
     let p = distribution.(i) in
-    distribution.(i) <- Float.((p * (1. - learning_coef)) + (value * learning_coef)))
+    distribution.(i)
+    <- Float.((p * (1. - learning_coef)) + (of_int value * learning_coef)))
 ;;
 
 let mutate_distribution distribution mutation_p mutation_disturbance_coef =
@@ -36,24 +49,43 @@ let mutate_distribution distribution mutation_p mutation_disturbance_coef =
     else
       Float.(
         (p * (1. - mutation_disturbance_coef))
-        + (Random.coin_flip () * mutation_disturbance_coef)))
+        + (of_int (Random.coin_flip ()) * mutation_disturbance_coef)))
+;;
+
+let get_info distribution population evaluate =
+  { Output.Iteration.distribution = Array.copy distribution
+  ; best_individual = Population.best_individual population ~f:evaluate
+  ; mean_individual = Population.mean_individual population ~f:evaluate
+  ; worst_individual = Population.worst_individual population ~f:evaluate
+  }
 ;;
 
 let pbil
-  { Input.gene_length
+  { Input.chromosome_length
   ; population_size
   ; learning_coef
   ; mutation_p
   ; mutation_disturbance_coef
   ; evaluate
   ; termination_condition
+  ; gather_verbose_info
   }
   =
-  let distribution = Array.init gene_length ~f:(fun _ -> 0.5) in
-  let population = Population.random population_size ~gene_length in
-  let rec pbil_aux iteration distribution population =
+  Log.info_s
+    [%message
+      "Running PBIL"
+        (chromosome_length : int)
+        (population_size : int)
+        (learning_coef : float)
+        (mutation_p : float)
+        (mutation_disturbance_coef : float)
+        (gather_verbose_info : bool)];
+  let t0 = Time_ns.now () in
+  let distribution = Array.init chromosome_length ~f:(fun _ -> 0.5) in
+  let population = Population.random population_size ~chromosome_length in
+  let rec pbil_aux iteration ~distribution ~population ~iterations_info =
     match termination_condition population ~iteration with
-    | true -> distribution, population
+    | true -> { Output.distribution; population; iterations_info }
     | false ->
       let best_individual = Population.best_individual population ~f:evaluate in
       adjust_distribution distribution best_individual learning_coef;
@@ -62,8 +94,28 @@ let pbil
         List.init population_size ~f:(fun _ -> Individual.derive distribution)
         |> Population.create_exn
       in
-      pbil_aux (iteration + 1) distribution population
+      let iterations_info =
+        match iterations_info with
+        | None -> None
+        | Some info -> Some (get_info distribution population evaluate :: info)
+      in
+      pbil_aux (iteration + 1) ~distribution ~population ~iterations_info
+  in
+  let output =
+    pbil_aux
+      0
+      ~distribution
+      ~population
+      ~iterations_info:
+        (if gather_verbose_info
+         then Some [ get_info distribution population evaluate ]
+         else None)
+  in
+  let t1 = Time_ns.now () in
+  let output =
+    { output with iterations_info = Option.map output.iterations_info ~f:List.rev }
   in
-  let distribution, population = pbil_aux 0 distribution population in
-  { Output.distribution; population }
+  Log.info_s
+    [%message "Pbil finished. Elapsed time" (Time_ns.diff t1 t0 : Time_ns.Span.t)];
+  output
 ;;
diff --git a/algorithms/pbil.mli b/algorithms/pbil.mli
index cea76bb..85fc6e4 100644
--- a/algorithms/pbil.mli
+++ b/algorithms/pbil.mli
@@ -3,23 +3,35 @@ open Genetics_lib
 
 module Input : sig
   type t =
-    { gene_length : int
+    { chromosome_length : int
     ; population_size : int
     ; learning_coef : float
     ; mutation_p : float
     ; mutation_disturbance_coef : float
-    ; evaluate : Individual.t -> float
+    ; evaluate : Individual.t -> int
     ; termination_condition : Population.t -> iteration:int -> bool
+    ; gather_verbose_info : bool
     }
   [@@deriving sexp_of]
 end
 
 module Output : sig
+  module Iteration : sig
+    type t =
+      { distribution : float array
+      ; best_individual : Individual.t
+      ; mean_individual : Individual.t
+      ; worst_individual : Individual.t
+      }
+    [@@deriving sexp]
+  end
+
   type t =
     { distribution : float array
     ; population : Population.t
+    ; iterations_info : Iteration.t list option
     }
-  [@@deriving sexp_of]
+  [@@deriving sexp]
 end
 
 val pbil : Input.t -> Output.t
diff --git a/app/pbil_benchmark/src/evaluators.ml b/app/pbil_benchmark/src/evaluators.ml
index 4c39e25..d918f89 100644
--- a/app/pbil_benchmark/src/evaluators.ml
+++ b/app/pbil_benchmark/src/evaluators.ml
@@ -1,30 +1,32 @@
 open! Core
 open Genetics_lib
 
-let one_max individual = Individual.gene individual |> Gene.fold ~init:0. ~f:Float.( + )
+let one_max individual =
+  Individual.chromosome individual |> Chromosome.fold ~init:0 ~f:( + )
+;;
 
 let deceptive_one_max individual =
   let result = one_max individual in
-  if Float.(result = 0.) then Individual.length individual + 1 |> Float.of_int else result
+  if result = 0 then Individual.length individual + 1 else result
 ;;
 
 let k_deceptive_one_max k individual =
-  Gene.foldi (Individual.gene individual) ~init:(0., 0) ~f:(fun i (sum, zero_cnt) value ->
-    let zero_cnt = if Float.(value = 0.) then zero_cnt + 1 else zero_cnt in
-    let sum =
-      if i % k = k - 1 && zero_cnt = k then sum +. Float.of_int (k + 1) else sum +. value
-    in
-    let zero_cnt = if i % k = k - 1 then 0 else zero_cnt in
-    sum, zero_cnt)
+  Chromosome.foldi
+    (Individual.chromosome individual)
+    ~init:(0, 0)
+    ~f:(fun i (sum, zero_cnt) value ->
+      let zero_cnt = if value = 0 then zero_cnt + 1 else zero_cnt in
+      let sum = if i % k = k - 1 && zero_cnt = k then sum + k + 1 else sum + value in
+      let zero_cnt = if i % k = k - 1 then 0 else zero_cnt in
+      sum, zero_cnt)
   |> Tuple2.get1
 ;;
 
 let%expect_test "evaluators" =
   let individuals =
     [ [ 1; 1; 1; 1; 1; 1 ]; [ 0; 0; 0; 0; 0; 0 ]; [ 0; 0; 1; 1; 0; 0 ] ]
-    |> List.map ~f:(List.map ~f:Float.of_int)
     |> List.map ~f:Array.of_list
-    |> List.map ~f:Gene.create
+    |> List.map ~f:Chromosome.create_exn
     |> List.map ~f:Individual.create
   in
   let evaluators =
@@ -35,8 +37,9 @@ let%expect_test "evaluators" =
   in
   List.iter individuals ~f:(fun individual ->
     List.iter evaluators ~f:(fun (name, evaluator) ->
-      print_s [%message name (individual : Individual.t) (evaluator individual : float)]));
-  [%expect{|
+      print_s [%message name (individual : Individual.t) (evaluator individual : int)]));
+  [%expect
+    {|
     (one_max (individual (1 1 1 1 1 1)) ("evaluator individual" 6))
     (deceptive_one_max (individual (1 1 1 1 1 1)) ("evaluator individual" 6))
     (2_deceptive_one_max (individual (1 1 1 1 1 1)) ("evaluator individual" 6))
diff --git a/app/pbil_benchmark/src/evaluators.mli b/app/pbil_benchmark/src/evaluators.mli
index 4889e9a..80a8dbf 100644
--- a/app/pbil_benchmark/src/evaluators.mli
+++ b/app/pbil_benchmark/src/evaluators.mli
@@ -1,6 +1,6 @@
 open! Core
 open Genetics_lib
 
-val one_max : Individual.t -> float
-val deceptive_one_max : Individual.t -> float
-val k_deceptive_one_max : int -> Individual.t -> float
+val one_max : Individual.t -> int
+val deceptive_one_max : Individual.t -> int
+val k_deceptive_one_max : int -> Individual.t -> int
diff --git a/app/pbil_benchmark/src/main.ml b/app/pbil_benchmark/src/main.ml
index d6053d9..82dd1e4 100644
--- a/app/pbil_benchmark/src/main.ml
+++ b/app/pbil_benchmark/src/main.ml
@@ -1,7 +1,53 @@
 open! Core
 open Algorithms_lib
+open Genetics_lib
 open Utils_lib
 
+let plot_best_mean_max
+  h
+  (input : Pbil.Input.t)
+  (iterations_info : Pbil.Output.Iteration.t list)
+  trials
+  =
+  Plotting.plot_multiple_int64
+    ~h
+    ~name:
+      [%string
+        "pbil_best_mean_worst_%{input.chromosome_length#Int}_%{input.population_size#Int}_%{trials#Int}_%{input.learning_coef#Float}_%{input.mutation_p#Float}_%{input.mutation_disturbance_coef#Float}"]
+    [ List.map iterations_info ~f:(fun iteration ->
+        input.evaluate iteration.best_individual |> Int64.of_int)
+      |> Array.of_list
+    ; List.map iterations_info ~f:(fun iteration ->
+        input.evaluate iteration.mean_individual |> Int64.of_int)
+      |> Array.of_list
+    ; List.map iterations_info ~f:(fun iteration ->
+        input.evaluate iteration.worst_individual |> Int64.of_int)
+      |> Array.of_list
+    ]
+    ~pen_size:0.5
+;;
+
+let plot_distribution
+  h
+  (input : Pbil.Input.t)
+  (iterations_info : Pbil.Output.Iteration.t list)
+  trials
+  =
+  let get_dim dim (iterations_info : Pbil.Output.Iteration.t list) =
+    List.map iterations_info ~f:(fun iteration -> iteration.distribution.(dim))
+    |> Array.of_list
+  in
+  let first_iteration = List.hd_exn iterations_info in
+  let length = Array.length first_iteration.distribution in
+  Plotting.plot_multiple_float
+    ~h
+    ~name:
+      [%string
+        "pbil_distributions_%{input.chromosome_length#Int}_%{input.population_size#Int}_%{trials#Int}_%{input.learning_coef#Float}_%{input.mutation_p#Float}_%{input.mutation_disturbance_coef#Float}"]
+    (List.init length ~f:(fun i -> get_dim i iterations_info))
+    ~pen_size:0.5
+;;
+
 let create_command evaluator =
   Command.basic_or_error
     ~summary:"one_max evaluation"
@@ -18,21 +64,43 @@ let create_command evaluator =
          "mutation-disturbance-coef"
          (required float)
          ~doc:"FLOAT mutation disturbance coefficient"
+     and gather_verbose_info =
+       flag
+         "gather-verbose-info"
+         no_arg
+         ~doc:"Gather verbose information about iterations."
      in
      fun () ->
-       let output =
-         Pbil.pbil
-           { Pbil.Input.gene_length = n
-           ; population_size
-           ; learning_coef
-           ; mutation_p
-           ; mutation_disturbance_coef
-           ; evaluate = evaluator
-           ; termination_condition = (fun _ ~iteration -> iteration >= trials)
-           }
+       let input =
+         { Pbil.Input.chromosome_length = n
+         ; population_size
+         ; learning_coef
+         ; mutation_p
+         ; mutation_disturbance_coef
+         ; evaluate = evaluator
+         ; termination_condition = (fun _ ~iteration -> iteration >= trials)
+         ; gather_verbose_info
+         }
        in
-       Log.info_s [%sexp (output : Pbil.Output.t)];
-       Ok ())
+       let output = Pbil.pbil input in
+       let best_individual = Population.best_individual output.population ~f:evaluator in
+       let best_result = evaluator best_individual in
+       Log.info_s
+         [%message
+           "Results"
+             (output.distribution : float array)
+             (best_individual : Individual.t)
+             (best_result : int)];
+       match output.iterations_info with
+       | None -> Ok ()
+       | Some iterations_info ->
+         let h = Owl_plplot.Plot.create ~m:1 ~n:2 "pbil.jpg" in
+         Owl_plplot.Plot.subplot h 0 0;
+         let%bind.Or_error h = plot_best_mean_max h input iterations_info trials in
+         Owl_plplot.Plot.subplot h 0 1;
+         let%bind.Or_error h = plot_distribution h input iterations_info trials in
+         Owl_plplot.Plot.output h;
+         Ok ())
 ;;
 
 let command =
diff --git a/app/qap/src/qap.ml b/app/qap/src/qap.ml
index c52288d..c749c22 100644
--- a/app/qap/src/qap.ml
+++ b/app/qap/src/qap.ml
@@ -87,11 +87,12 @@ let random_sampling ~n ~trials =
   done;
   let t1 = Time_ns.now () in
   print_endline [%string "Elapsed time: %{Time_ns.diff t1 t0#Time_ns.Span}"];
-  Plotting.histogram_int64
-    ~pen_size:3.
-    ~bin:100
-    ~name:[%string "qap_random_sampling_%{n#Int}_trials_%{trials#Int}"]
-    costs;
+  Owl_plplot.Plot.output
+    (Plotting.histogram_int64
+       ~pen_size:3.
+       ~bin:100
+       ~name:[%string "qap_random_sampling_%{n#Int}_trials_%{trials#Int}"]
+       costs);
   costs
 ;;
 
@@ -104,11 +105,12 @@ let simulated_annealing ~n ~trials ~radius ~alpha =
     Algorithms_lib.Simulated_annealing.simulated_annealing
       { trials; radius; alpha; objective_function; initialize_sample; random_neighbour }
   in
-  Plotting.simple_plot_int64
-    ~name:
-      [%string
-        "qap_simulated_annealing_%{n#Int}_trials_%{trials#Int}_alpha_%{alpha#Float}_radius_%{radius#Int}"]
-    result.costs;
+  Owl_plplot.Plot.output
+    (Plotting.simple_plot_int64
+       ~name:
+         [%string
+           "qap_simulated_annealing_%{n#Int}_trials_%{trials#Int}_alpha_%{alpha#Float}_radius_%{radius#Int}"]
+       result.costs);
   Log.info_s
     [%message
       (Array.min_elt result.costs ~compare:Int64.compare : Int64.t option)
diff --git a/lib/genetics/chromosome.ml b/lib/genetics/chromosome.ml
new file mode 100644
index 0000000..e9dd1ce
--- /dev/null
+++ b/lib/genetics/chromosome.ml
@@ -0,0 +1,22 @@
+open! Core
+open Owl
+
+type t = int array [@@deriving sexp]
+
+let create arr =
+  match Array.for_all arr ~f:(fun i -> i = 1 || i = 0) with
+  | false ->
+    Or_error.error_s
+      [%message
+        "Cannot create a chromosome, not all values in array are 0 or 1" (arr : int array)]
+  | true -> Ok arr
+;;
+
+let create_exn arr = create arr |> Or_error.ok_exn
+let of_list = Array.of_list
+let length = Array.length
+let fold = Array.fold
+let foldi = Array.foldi
+let init length ~value = Array.init length ~f:(fun _ -> value)
+let random length = Array.init length ~f:(fun _ -> Stats.uniform_int_rvs ~a:0 ~b:1)
+let iteri = Array.iteri
diff --git a/lib/genetics/chromosome.mli b/lib/genetics/chromosome.mli
new file mode 100644
index 0000000..44a1fc6
--- /dev/null
+++ b/lib/genetics/chromosome.mli
@@ -0,0 +1,16 @@
+open! Core
+
+type t [@@deriving sexp]
+
+(* TODO: don't use int explicitely here: write this as a module type with an abstract array type. *)
+val create_exn : int array -> t
+val create : int array -> t Or_error.t
+val of_list : int list -> t
+val length : t -> int
+val fold : t -> init:'a -> f:('a -> int -> 'a) -> 'a
+val foldi : t -> init:'a -> f:(int -> 'a -> int -> 'a) -> 'a
+val init : int -> value:int -> t
+val iteri : t -> f:(int -> int -> unit) -> unit
+
+(** Create a random gene. *)
+val random : int -> t
diff --git a/lib/genetics/gene.ml b/lib/genetics/gene.ml
deleted file mode 100644
index d940e69..0000000
--- a/lib/genetics/gene.ml
+++ /dev/null
@@ -1,13 +0,0 @@
-open! Core
-open Owl
-
-type t = float array [@@deriving sexp]
-
-let create = Fn.id
-let of_list = Array.of_list
-let length = Array.length
-let fold = Array.fold
-let foldi = Array.foldi
-let init length ~value = Array.init length ~f:(fun _ -> value)
-let random length = Array.init length ~f:(fun _ -> Stats.uniform_rvs ~a:0. ~b:1.)
-let iteri = Array.iteri
diff --git a/lib/genetics/gene.mli b/lib/genetics/gene.mli
deleted file mode 100644
index 70ea69e..0000000
--- a/lib/genetics/gene.mli
+++ /dev/null
@@ -1,14 +0,0 @@
-open! Core
-
-type t [@@deriving sexp]
-
-val create : float array -> t
-val of_list : float list -> t
-val length : t -> int
-val fold : t -> init:'a -> f:('a -> float -> 'a) -> 'a
-val foldi : t -> init:'a -> f:(int -> 'a -> float -> 'a) -> 'a
-val init : int -> value:float -> t
-val iteri : t -> f:(int -> float -> unit) -> unit
-
-(** Create a random gene. *)
-val random : int -> t
diff --git a/lib/genetics/genetics_lib.ml b/lib/genetics/genetics_lib.ml
index 4e3f80b..8042989 100644
--- a/lib/genetics/genetics_lib.ml
+++ b/lib/genetics/genetics_lib.ml
@@ -1,3 +1,3 @@
-module Gene = Gene
+module Chromosome = Chromosome
 module Population = Population
 module Individual = Individual
diff --git a/lib/genetics/individual.ml b/lib/genetics/individual.ml
index bf08647..88f0655 100644
--- a/lib/genetics/individual.ml
+++ b/lib/genetics/individual.ml
@@ -1,13 +1,13 @@
 open! Core
 open Utils_lib
 
-type t = Gene.t [@@deriving sexp_of]
+type t = Chromosome.t [@@deriving sexp]
 
 let create = Fn.id
-let length = Gene.length
-let gene = Fn.id
-let random = Gene.random
+let length = Chromosome.length
+let chromosome = Fn.id
+let random = Chromosome.random
 
 let derive distribution =
-  Array.map distribution ~f:(fun p -> Random.binary p) |> Gene.create |> create
+  Array.map distribution ~f:(fun p -> Random.binary p) |> Chromosome.create_exn |> create
 ;;
diff --git a/lib/genetics/individual.mli b/lib/genetics/individual.mli
index ba4ac1c..aeb3c94 100644
--- a/lib/genetics/individual.mli
+++ b/lib/genetics/individual.mli
@@ -1,11 +1,11 @@
 open! Core
 
-type t [@@deriving sexp_of]
+type t [@@deriving sexp]
 
 (* For now it's an identity function, but maybe later there will be need for some additional information. *)
-val create : Gene.t -> t
+val create : Chromosome.t -> t
 val length : t -> int
-val gene : t -> Gene.t
+val chromosome : t -> Chromosome.t
 val random : int -> t
 
 (** Creates the individual based on the genes probabilities. Returns error if the array has wrong length. *)
diff --git a/lib/genetics/population.ml b/lib/genetics/population.ml
index 76017ca..668e712 100644
--- a/lib/genetics/population.ml
+++ b/lib/genetics/population.ml
@@ -1,6 +1,6 @@
 open! Core
 
-type t = Individual.t list [@@deriving sexp_of]
+type t = Individual.t list [@@deriving sexp]
 
 let create individuals =
   match individuals with
@@ -23,29 +23,41 @@ let best_individual t ~f =
     | individual :: rest ->
       let value = f individual in
       let best_individual, best_value =
-        if Float.(value > best_value)
-        then individual, value
-        else best_individual, best_value
+        if value > best_value then individual, value else best_individual, best_value
       in
       best_individual_aux best_individual best_value rest
   in
-  best_individual_aux (List.hd_exn t) 0. t
+  let fst = List.hd_exn t in
+  best_individual_aux fst (f fst) t
 ;;
 
-let random length ~gene_length =
-  List.init length ~f:(fun _ -> Individual.create (Gene.random gene_length))
+let worst_individual t ~f = best_individual t ~f:(fun ind -> -f ind)
+
+let mean_individual t ~f =
+  let individuals =
+    List.map t ~f:(fun ind -> f ind, ind)
+    |> List.sort ~compare:(fun (v1, _) (v2, _) -> Int.compare v1 v2)
+  in
+  List.split_n individuals (List.length t / 2)
+  |> Tuple2.get2
+  |> List.hd_exn
+  |> Tuple2.get2
+;;
+
+let random length ~chromosome_length =
+  List.init length ~f:(fun _ -> Individual.create (Chromosome.random chromosome_length))
 ;;
 
 let%expect_test "best_individual" =
   let population =
-    [ [| 0.; 0.1; 0.2; 0.3 |]; [| 1.; 1.; 1.; 1. |]; [| 0.; 0.; 0.; 0. |] ]
-    |> List.map ~f:Gene.create
+    [ [| 0; 0; 1; 1 |]; [| 1; 1; 1; 1 |]; [| 0; 0; 0; 0 |] ]
+    |> List.map ~f:Chromosome.create_exn
     |> List.map ~f:Individual.create
     |> create
   and empty_population = create []
   and different_lengths_population =
-    [ [| 1. |]; [| 1.; 2. |] ]
-    |> List.map ~f:Gene.create
+    [ [| 1 |]; [| 1; 0 |] ]
+    |> List.map ~f:Chromosome.create_exn
     |> List.map ~f:Individual.create
     |> create
   in
@@ -54,15 +66,20 @@ let%expect_test "best_individual" =
   print_s [%sexp (different_lengths_population : t Or_error.t)];
   [%expect
     {|
-    (Ok ((0 0.1 0.2 0.3) (1 1 1 1) (0 0 0 0)))
+    (Ok ((0 0 1 1) (1 1 1 1) (0 0 0 0)))
     (Error "Cannot create an empty population.")
     (Error
-     ("Not all individuals have the same length." (individuals ((1) (1 2))))) |}];
+     ("Not all individuals have the same length." (individuals ((1) (1 0))))) |}];
   let population = Or_error.ok_exn population in
-  let best_ind =
-    best_individual population ~f:(fun ind ->
-      Gene.fold (Individual.gene ind) ~init:0. ~f:( +. ))
-  in
+  let evaluation ind = Chromosome.fold (Individual.chromosome ind) ~init:0 ~f:( + ) in
+  let best_ind = best_individual population ~f:evaluation
+  and mean_ind = mean_individual population ~f:evaluation
+  and worst_ind = worst_individual population ~f:evaluation in
   print_s [%sexp (best_ind : Individual.t)];
-  [%expect {| (1 1 1 1) |}]
+  print_s [%sexp (mean_ind : Individual.t)];
+  print_s [%sexp (worst_ind : Individual.t)];
+  [%expect {|
+    (1 1 1 1)
+    (0 0 1 1)
+    (0 0 0 0) |}]
 ;;
diff --git a/lib/genetics/population.mli b/lib/genetics/population.mli
index 950984a..2877037 100644
--- a/lib/genetics/population.mli
+++ b/lib/genetics/population.mli
@@ -1,10 +1,14 @@
 open! Core
 
-type t [@@deriving sexp_of]
+type t [@@deriving sexp]
 
 (** Returns error if the list is empty. *)
 val create : Individual.t list -> t Or_error.t
 
 val create_exn : Individual.t list -> t
-val best_individual : t -> f:(Individual.t -> float) -> Individual.t
-val random : int -> gene_length:int -> t
+val random : int -> chromosome_length:int -> t
+
+(* TODO: this type Individual.t -> int should be a type of some module "Evaluator" *)
+val best_individual : t -> f:(Individual.t -> int) -> Individual.t
+val mean_individual : t -> f:(Individual.t -> int) -> Individual.t
+val worst_individual : t -> f:(Individual.t -> int) -> Individual.t
diff --git a/lib/utils/plotting.ml b/lib/utils/plotting.ml
index d4a8af3..c5dcb31 100644
--- a/lib/utils/plotting.ml
+++ b/lib/utils/plotting.ml
@@ -2,22 +2,27 @@ open! Core
 open Owl_plplot
 
 let initialize_plot
+  ?h
   ?(xlabel = "values")
   ?(ylabel = "")
-  ?(font_size = 8.)
+  ?(font_size = 4.)
   ?(pen_size = 1.)
   ~name
   ()
   =
-  let h = Plot.create [%string "%{name}.jpg"] in
+  let h =
+    Option.value_or_thunk h ~default:(fun () -> Plot.create [%string "%{name}.jpg"])
+  in
   Plot.set_xlabel h xlabel;
   Plot.set_ylabel h ylabel;
   Plot.set_font_size h font_size;
   Plot.set_pen_size h pen_size;
+  Plot.set_title h name;
   h
 ;;
 
 let simple_plot_floats
+  ?h
   ?xlabel
   ?ylabel
   ?font_size
@@ -27,7 +32,7 @@ let simple_plot_floats
   values
   =
   Log.info_s [%message "Plotting" name];
-  let h = initialize_plot ?xlabel ?ylabel ?font_size ?pen_size ~name () in
+  let h = initialize_plot ?h ?xlabel ?ylabel ?font_size ?pen_size ~name () in
   let length = Array.length values in
   let indices = Array.init length ~f:Float.of_int in
   Plot.plot
@@ -35,11 +40,12 @@ let simple_plot_floats
     ~spec
     (Owl_dense_matrix.D.of_array indices 1 length)
     (Owl_dense_matrix.D.of_array values 1 length);
-  Plot.output h
+  h
 ;;
 
-let simple_plot_int64 ?xlabel ?ylabel ?font_size ?pen_size ?spec ~name values =
+let simple_plot_int64 ?h ?xlabel ?ylabel ?font_size ?pen_size ?spec ~name values =
   simple_plot_floats
+    ?h
     ?xlabel
     ?ylabel
     ?font_size
@@ -50,6 +56,7 @@ let simple_plot_int64 ?xlabel ?ylabel ?font_size ?pen_size ?spec ~name values =
 ;;
 
 let histogram_int64
+  ?h
   ?xlabel
   ?ylabel
   ?font_size
@@ -59,7 +66,7 @@ let histogram_int64
   ~name
   values
   =
-  let h = initialize_plot ?xlabel ?ylabel ?font_size ?pen_size ~name () in
+  let h = initialize_plot ?h ?xlabel ?ylabel ?font_size ?pen_size ~name () in
   Plot.histogram
     ~h
     ~bin
@@ -68,5 +75,56 @@ let histogram_int64
        (Array.map ~f:Int64.to_float values)
        1
        (Array.length values));
-  Plot.output h
+  h
+;;
+
+let plot_multiple_float ?h ?xlabel ?ylabel ?font_size ?pen_size ?specs ~name values =
+  let get () = Owl.Stats.uniform_int_rvs ~a:0 ~b:255 in
+  let specs =
+    Option.value_or_thunk specs ~default:(fun () ->
+      List.init (List.length values) ~f:(fun _ -> [ Plot.RGB (get (), get (), get ()) ]))
+  in
+  let%bind.Or_error length =
+    match values with
+    | [] ->
+      Or_error.error_s [%message "There is nothing to plot, values is empty list" name]
+    | values :: _ -> Ok (Array.length values)
+  in
+  let%bind.Or_error () =
+    match List.for_all values ~f:(fun values -> Array.length values = length) with
+    | false ->
+      Or_error.error_s
+        [%message
+          "Not all values have the same length, cannot plot." (values : float array list)]
+    | true -> Ok ()
+  in
+  let indices = Array.init length ~f:Float.of_int in
+  let h = initialize_plot ?h ?xlabel ?ylabel ?font_size ?pen_size ~name () in
+  match
+    List.iter2 specs values ~f:(fun spec values ->
+      Plot.plot
+        ~h
+        ~spec
+        (Owl_dense_matrix.D.of_array indices 1 length)
+        (Owl_dense_matrix.D.of_array values 1 length))
+  with
+  | Unequal_lengths ->
+    Or_error.error_s
+      [%message
+        "Specs have different size than the values"
+          (List.length specs : int)
+          (List.length values : int)]
+  | Ok _ -> Ok h
+;;
+
+let plot_multiple_int64 ?h ?xlabel ?ylabel ?font_size ?pen_size ?specs ~name values =
+  plot_multiple_float
+    ?h
+    ?xlabel
+    ?ylabel
+    ?font_size
+    ?pen_size
+    ?specs
+    ~name
+    (List.map values ~f:(Array.map ~f:Float.of_int64))
 ;;
diff --git a/lib/utils/plotting.mli b/lib/utils/plotting.mli
index 8af5afa..5a94832 100644
--- a/lib/utils/plotting.mli
+++ b/lib/utils/plotting.mli
@@ -1,28 +1,41 @@
 open! Core
 open Owl_plplot
 
+val initialize_plot
+  :  ?h:Plot.handle
+  -> ?xlabel:string
+  -> ?ylabel:string
+  -> ?font_size:float
+  -> ?pen_size:float
+  -> name:string
+  -> unit
+  -> Plot.handle
+
 val simple_plot_floats
-  :  ?xlabel:string
+  :  ?h:Plot.handle
+  -> ?xlabel:string
   -> ?ylabel:string
   -> ?font_size:float
   -> ?pen_size:float
   -> ?spec:Plot.spec list
   -> name:string
   -> float array
-  -> unit
+  -> Plot.handle
 
 val simple_plot_int64
-  :  ?xlabel:string
+  :  ?h:Plot.handle
+  -> ?xlabel:string
   -> ?ylabel:string
   -> ?font_size:float
   -> ?pen_size:float
   -> ?spec:Plot.spec list
   -> name:string
   -> int64 array
-  -> unit
+  -> Plot.handle
 
 val histogram_int64
-  :  ?xlabel:string
+  :  ?h:Plot.handle
+  -> ?xlabel:string
   -> ?ylabel:string
   -> ?font_size:float
   -> ?pen_size:float
@@ -30,4 +43,26 @@ val histogram_int64
   -> ?spec:Plot.spec list
   -> name:string
   -> int64 array
-  -> unit
+  -> Plot.handle
+
+val plot_multiple_float
+  :  ?h:Plot.handle
+  -> ?xlabel:string
+  -> ?ylabel:string
+  -> ?font_size:float
+  -> ?pen_size:float
+  -> ?specs:Plot.spec list list
+  -> name:string
+  -> float array list
+  -> Plot.handle Or_error.t
+
+val plot_multiple_int64
+  :  ?h:Plot.handle
+  -> ?xlabel:string
+  -> ?ylabel:string
+  -> ?font_size:float
+  -> ?pen_size:float
+  -> ?specs:Plot.spec list list
+  -> name:string
+  -> int64 array list
+  -> Plot.handle Or_error.t
diff --git a/lib/utils/random.ml b/lib/utils/random.ml
index b3fd05f..a1362df 100644
--- a/lib/utils/random.ml
+++ b/lib/utils/random.ml
@@ -2,5 +2,5 @@ open! Core
 open Owl
 
 let uniform_01 () = Stats.uniform_rvs ~a:0. ~b:1.
-let binary p = if Float.(uniform_01 () < p) then 1. else 0.
+let binary p = if Float.(uniform_01 () < p) then 1 else 0
 let coin_flip () = binary 0.5
diff --git a/lib/utils/random.mli b/lib/utils/random.mli
index 7eb339c..17361aa 100644
--- a/lib/utils/random.mli
+++ b/lib/utils/random.mli
@@ -3,6 +3,6 @@ open! Core
 val uniform_01 : unit -> float
 
 (** Return 1 with probability p, 0 with probability 1 - p. *)
-val binary : float -> float
+val binary : float -> int
 
-val coin_flip : unit -> float
+val coin_flip : unit -> int