Commit 3d706ee4 authored by Philippe Veber's avatar Philippe Veber
Browse files

tk/Mutsel_simulator_cpg: cosmetic changes

parent 6d715ccb
open Core_kernel
open Phylogenetics
module Mutsel_cpg = struct
module NSCodon = Codon.Universal_genetic_code.NS
module NSCodon = Codon.Universal_genetic_code.NS
module Evolution_model = struct
type context_CpG = {
lh: Nucleotide.repr option;
rh: Nucleotide.repr option;
......@@ -28,25 +27,26 @@ module Mutsel_cpg = struct
let context_CpG_of_codons p ?lh ?rh () =
let open Option.Monad_infix in
{
lh = lh >>| (fun lh -> inspect_codon lh |> trd3);
rh = rh >>| (fun rh -> inspect_codon rh |> fst3);
codon = inspect_codon p;
lh = lh >>| inspect_codon >>| trd3 ;
rh = rh >>| inspect_codon >>| fst3 ;
codon = inspect_codon p ;
}
let fitness_of_profile ?(beta = 1.) =
Amino_acid.Vector.map ~f:(fun x -> beta *. Float.log x)
let flat_fitness () =
Amino_acid.Vector.init (fun _ -> 1. /. float Amino_acid.card)
|> fitness_of_profile
let make_param state pos _branch =
let fitness_of_profile ?(beta = 1.) =
Amino_acid.Vector.map ~f:(fun x -> beta *. Float.log x)
in
let flat_fitness () =
Amino_acid.Vector.init (fun _ -> 1. /. float Amino_acid.card)
|> fitness_of_profile
in
let context_CpG =
if pos = 0
then context_CpG_of_codons state.(pos) ~rh:state.(pos+1) ()
else if pos = (Array.length state - 1 )
then context_CpG_of_codons state.(pos) ~lh:state.(pos-1) ()
else context_CpG_of_codons state.(pos) ~lh:state.(pos -1) ~rh:state.(pos + 1) ()
if pos = 0 then
context_CpG_of_codons state.(pos) ~rh:state.(pos+1) ()
else if pos = (Array.length state - 1 ) then
context_CpG_of_codons state.(pos) ~lh:state.(pos-1) ()
else
context_CpG_of_codons state.(pos) ~lh:state.(pos -1) ~rh:state.(pos + 1) ()
in
let pi = Nucleotide.flat_profile () in
let rho = Linear_algebra.Lacaml.Vector.init 6 ~f:(fun _ -> 1. /. 6.) in
......@@ -63,12 +63,12 @@ module Mutsel_cpg = struct
let rate_increase_CpG (pos, x_a, x_b) context =
match Nucleotide.(pos, inspect x_a, inspect x_b), context with
| (0, C, T), { codon=(_,G, _) ; _ }
| (1, C, T), { codon=(_,_, G) ; _ }
| (2, C, T), { rh=Some G ; _ }
| (0, G, A), { lh=Some C ; _ }
| (1, G, A), { codon=(C, _, _) ; _ }
| (2, G, A), { codon=(_,C,_) ; _ } -> true
| (0, C, T), { codon = (_,G, _) ; _ }
| (1, C, T), { codon = (_,_, G) ; _ }
| (2, C, T), { rh = Some G ; _ }
| (0, G, A), { lh = Some C ; _ }
| (1, G, A), { codon = (C, _, _) ; _ }
| (2, G, A), { codon = (_,C,_) ; _ } -> true
| _ -> false
let fixation_probability delta =
......@@ -118,6 +118,4 @@ module Mutsel_cpg = struct
|> NSCodon.Vector.normalize
end
module M(BI: Simulator.Branch_info) = struct
include Simulator.Make(Mutsel_cpg.NSCodon)(Mutsel_cpg)(BI)
end
module Make(BI: Simulator.Branch_info) = Simulator.Make(NSCodon)(Evolution_model)(BI)
open Phylogenetics
module NSCodon = Codon.Universal_genetic_code.NS
module Mutsel_cpg : sig
module NSCodon = Codon.Universal_genetic_code.NS
module Evolution_model : sig
type context_CpG = {
lh: Nucleotide.repr option;
rh: Nucleotide.repr option;
......@@ -26,9 +24,8 @@ module Mutsel_cpg : sig
val stationary_distribution: param -> NSCodon.vector
val make_param: NSCodon.t array -> int -> 'a -> param
end
module M(BI: Simulator.Branch_info) : sig
include module type of Simulator.Make(Mutsel_cpg.NSCodon)(Mutsel_cpg)(BI)
module Make(BI: Simulator.Branch_info) : sig
include module type of Simulator.Make(NSCodon)(Evolution_model)(BI)
end
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment