fork download
copy 
  1. using System;
  2. using System.Collections.Generic;
  3. using System.Collections.ObjectModel;
  4. using System.Linq;
  5. using System.Text;
  6.  
  7. public class Locus : IEnumerable<double>
  8. {
  9.    private readonly HashSet<double> mutations;
  10.  
  11.    public Locus()
  12.    {
  13.       this.mutations = new HashSet<double>();
  14.    }
  15.  
  16.    public Locus(IEnumerable<double> mutations)
  17.    {
  18.       this.mutations = new HashSet<double>(mutations);
  19.    }
  20.  
  21.    public static Locus operator |(Locus c1, Locus c2)
  22.    {
  23.       if (c1 == null)
  24.       {
  25.          throw new ArgumentNullException("c1");
  26.       }
  27.       if (c2 == null)
  28.       {
  29.          throw new ArgumentNullException("c2");
  30.       }
  31.       return new Locus(c1.mutations.Union(c2.mutations));
  32.    }
  33.  
  34.    public void AddMutation(double site)
  35.    {
  36.       this.mutations.Add(site);
  37.    }
  38.  
  39.    public bool IsMutant(double site)
  40.    {
  41.       return this.mutations.Contains(site);
  42.    }
  43.  
  44.    public IEnumerator<double> GetEnumerator()
  45.    {
  46.       return this.mutations.GetEnumerator();
  47.    }
  48.    System.Collections.IEnumerator System.Collections.IEnumerable.GetEnumerator()
  49.    {
  50.       return GetEnumerator();
  51.    }
  52. }
  53.  
  54. public class CoalescentTreeNode
  55. {
  56.    public string Name { get; set; }
  57.  
  58.    private CoalescentTree tree;
  59.    public CoalescentTree Tree
  60.    {
  61.       get { return this.tree; }
  62.       internal set { this.tree = value; }
  63.    }
  64.  
  65.    private CoalescentTreeNode left;
  66.    public CoalescentTreeNode DescendantLeft
  67.    {
  68.       get { return this.left; }
  69.       internal set { this.left = value; }
  70.    }
  71.  
  72.    private CoalescentTreeNode right;
  73.    public CoalescentTreeNode DescendantRight
  74.    {
  75.       get { return this.right; }
  76.       internal set { this.right = value; }
  77.    }
  78.  
  79.    private CoalescentTreeNode ancestor;
  80.    public CoalescentTreeNode Ancestor
  81.    {
  82.       get { return this.ancestor; }
  83.       internal set { this.ancestor = value; }
  84.    }
  85.  
  86.    private double time;
  87.    public double WaitingTime
  88.    {
  89.       get { return this.time; }
  90.       internal set { this.time = value; }
  91.    }
  92.  
  93.    private Locus locus;
  94.    public Locus Locus
  95.    {
  96.       get
  97.       {
  98.          if (this.locus == null)
  99.          {
  100.             this.locus = new Locus();
  101.          }
  102.          return this.locus;
  103.       }
  104.       internal set { this.locus = value; }
  105.    }
  106.  
  107.    public bool IsLeaf
  108.    {
  109.       get
  110.       {
  111.          return DescendantLeft == null && DescendantRight == null;
  112.       }
  113.    }
  114.  
  115.    internal bool HasDescendant
  116.    {
  117.       get
  118.       {
  119.          return DescendantLeft != null && DescendantRight != null;
  120.       }
  121.    }
  122.  
  123.    public CoalescentTreeNode()
  124.    {
  125.    }
  126.  
  127.    public void AddMutation(double site)
  128.    {
  129.       foreach (var node in EnumerateAll(this))
  130.       {
  131.          node.Locus.AddMutation(site);
  132.       }
  133.    }
  134.  
  135.    internal static IEnumerable<CoalescentTreeNode> EnumerateAll(CoalescentTreeNode baseNode)
  136.    {
  137.       if (baseNode == null)
  138.       {
  139.          throw new ArgumentNullException("baseNode");
  140.       }
  141.  
  142.       yield return baseNode;
  143.       if (baseNode.DescendantLeft != null)
  144.       {
  145.          foreach (var descendant in EnumerateAll(baseNode.DescendantLeft))
  146.          {
  147.             yield return descendant;
  148.          }
  149.       }
  150.       if (baseNode.DescendantRight != null)
  151.       {
  152.          foreach (var descendant in EnumerateAll(baseNode.DescendantRight))
  153.          {
  154.             yield return descendant;
  155.          }
  156.       }
  157.    }
  158.  
  159.    internal static IEnumerable<CoalescentTreeNode> EnumerateLeaves(CoalescentTreeNode baseNode)
  160.    {
  161.       if (baseNode == null)
  162.       {
  163.          throw new ArgumentNullException("baseNode");
  164.       }
  165.  
  166.       CoalescentTreeNode node = baseNode;
  167.       while (node.DescendantLeft != null)
  168.       {
  169.          node = node.DescendantLeft;
  170.       }
  171.       yield return node;
  172.       while (node != baseNode)
  173.       {
  174.          node = node.Ancestor;
  175.          if (node.DescendantRight != null)
  176.          {
  177.             foreach (var value in EnumerateLeaves(node.DescendantRight))
  178.             {
  179.                yield return value;
  180.             }
  181.          }
  182.       }
  183.    }
  184.  
  185.    internal double GetTotalTime()
  186.    {
  187.       double sum = 0.0;
  188.       CoalescentTreeNode current = DescendantLeft;
  189.       while (current != null)
  190.       {
  191.          sum += current.WaitingTime;
  192.          current = current.DescendantLeft;
  193.       }
  194.       return sum;
  195.    }
  196. }
  197.  
  198. public class CoalescentTree : IEnumerable<CoalescentTreeNode>, IFormattable
  199. {
  200.    private readonly List<CoalescentTreeNode> nodes;
  201.  
  202.    public CoalescentTreeNode Root
  203.    {
  204.       get { return IsIncomplete ? null : this.nodes[0]; }
  205.    }
  206.  
  207.    public bool IsIncomplete
  208.    {
  209.       get { return this.nodes.Count != 1; }
  210.    }
  211.  
  212.    public CoalescentTree()
  213.    {
  214.       this.nodes = new List<CoalescentTreeNode>();
  215.    }
  216.  
  217.    internal CoalescentTree(IEnumerable<CoalescentTreeNode> nodes)
  218.    {
  219.       this.nodes = new List<CoalescentTreeNode>(nodes);
  220.    }
  221.  
  222.    public CoalescentTreeNode Coalesce(CoalescentTreeNode left, CoalescentTreeNode right, double time)
  223.    {
  224.       if (!Contains(left) || left.Ancestor != null)
  225.       {
  226.          throw new ArgumentException("left");
  227.       }
  228.       if (!Contains(right) || right.Ancestor != null)
  229.       {
  230.          throw new ArgumentException("right");
  231.       }
  232.       double leftTime = left.GetTotalTime();
  233.       if (leftTime > time)
  234.       {
  235.          throw new ArgumentException("left");
  236.       }
  237.       double rightTime = right.GetTotalTime();
  238.       if (rightTime > time)
  239.       {
  240.          throw new ArgumentException("right");
  241.       }
  242.  
  243.       left.WaitingTime = time - leftTime;
  244.       right.WaitingTime = time - rightTime;
  245.  
  246.       CoalescentTreeNode ancestor = new CoalescentTreeNode()
  247.       {
  248.          DescendantLeft = left,
  249.          DescendantRight = right,
  250.          Locus = left.Locus | right.Locus,
  251.       };
  252.       left.Ancestor = ancestor;
  253.       right.Ancestor = ancestor;
  254.  
  255.       Remove(left);
  256.       Remove(right);
  257.       Add(ancestor);
  258.       return ancestor;
  259.    }
  260.  
  261.    public void Add(CoalescentTreeNode node)
  262.    {
  263.       if (node == null)
  264.       {
  265.          throw new ArgumentNullException("node");
  266.       }
  267.       if (node.Tree != null)
  268.       {
  269.          throw new InvalidOperationException();
  270.       }
  271.       this.nodes.Add(node);
  272.       node.Tree = this;
  273.    }
  274.  
  275.    private bool Remove(CoalescentTreeNode node)
  276.    {
  277.       return this.nodes.Remove(node);
  278.    }
  279.  
  280.    private bool Contains(CoalescentTreeNode node)
  281.    {
  282.       return this.nodes.Contains(node);
  283.    }
  284.  
  285.    internal IList<CoalescentTreeNode> GetNodes()
  286.    {
  287.       return new ReadOnlyCollection<CoalescentTreeNode>(nodes);
  288.    }
  289.  
  290.    public override string ToString()
  291.    {
  292.       return ToString(null, null);
  293.    }
  294.  
  295.    public string ToString(string format, IFormatProvider formatProvider)
  296.    {
  297.       if (format != null && Char.ToUpper(format[0]) == 'A')
  298.       {
  299.          return ToAlignment();
  300.       }
  301.       else
  302.       {
  303.          if (Root.HasDescendant)
  304.          {
  305.             return string.Format("({0},{1});", ToString(Root.DescendantLeft, format, formatProvider), ToString(Root.DescendantRight, format, formatProvider));
  306.          }
  307.          return string.Format("{0}:{1};", Root.Name, Root.WaitingTime);
  308.       }
  309.    }
  310.  
  311.    private string ToString(CoalescentTreeNode node, string format, IFormatProvider formatProvider)
  312.    {
  313.       if (node == null)
  314.       {
  315.          return string.Empty;
  316.       }
  317.       if (node.HasDescendant)
  318.       {
  319.          return string.Format("({0},{1}):{2}", ToString(node.DescendantLeft, format, formatProvider), ToString(node.DescendantRight, format, formatProvider), node.WaitingTime.ToString(format, formatProvider));
  320.       }
  321.       return string.Format("{0}:{1}", node.Name, node.WaitingTime.ToString(format, formatProvider));
  322.    }
  323.  
  324.    private string ToAlignment()
  325.    {
  326.       StringBuilder alignment = new StringBuilder();
  327.       HashSet<double> allMutations = new HashSet<double>(this.SelectMany(node => node.Locus));
  328.       foreach (var node in this)
  329.       {
  330.          if (alignment.Length != 0)
  331.          {
  332.             alignment.AppendLine();
  333.          }
  334.          StringBuilder sequence = new StringBuilder();
  335.          foreach (var site in allMutations)
  336.          {
  337.             sequence.Append(node.Locus.IsMutant(site) ? '1' : '0');
  338.          }
  339.          alignment.Append(node.Name).Append(": ").Append(sequence.ToString());
  340.       }
  341.       return alignment.ToString();
  342.    }
  343.  
  344.    public IEnumerator<CoalescentTreeNode> GetEnumerator()
  345.    {
  346.       if (IsIncomplete)
  347.       {
  348.          throw new InvalidOperationException();
  349.       }
  350.       foreach (var node in CoalescentTreeNode.EnumerateLeaves(Root))
  351.       {
  352.          yield return node;
  353.       }
  354.    }
  355.    System.Collections.IEnumerator System.Collections.IEnumerable.GetEnumerator()
  356.    {
  357.       return GetEnumerator();
  358.    }
  359. }
  360.  
  361. public class CoalescentSimulator
  362. {
  363.    private readonly Random random;
  364.  
  365.    public CoalescentSimulator()
  366.    {
  367.       this.random = new Random();
  368.    }
  369.  
  370.    public CoalescentTree PerformSimulation(IEnumerable<CoalescentTreeNode> samples, double theta = 0.0)
  371.    {
  372.       if (theta < 0)
  373.       {
  374.          throw new ArgumentOutOfRangeException("theta");
  375.       }
  376.       CoalescentTree tree = new CoalescentTree(samples);
  377.  
  378.       double time = 0.0;
  379.       while (tree.IsIncomplete)
  380.       {
  381.          IList<CoalescentTreeNode> nodes = tree.GetNodes();
  382.          int count = nodes.Count;
  383.          time = GenerateNextCoalescentEventTime(time, count);
  384.          int leftIndex;
  385.          int rightIndex;
  386.          GetRandomTwoIndexes(count, out leftIndex, out rightIndex);
  387.          CoalescentTreeNode leftNode = nodes[leftIndex];
  388.          CoalescentTreeNode rightNode = nodes[rightIndex];
  389.          tree.Coalesce(leftNode, rightNode, time);
  390.          AddMutation(leftNode, theta);
  391.          AddMutation(rightNode, theta);
  392.       }
  393.  
  394.       return tree;
  395.    }
  396.  
  397.    private void GetRandomTwoIndexes(int count, out int index1, out int index2)
  398.    {
  399.       index1 = random.Next(count);
  400.       index2 = random.Next(count - 1);
  401.       if (index2 >= index1)
  402.       {
  403.          ++index2;
  404.       }
  405.    }
  406.  
  407.    private double GenerateNextCoalescentEventTime(double currentTime, int sampleCount)
  408.    {
  409.       double rate = (double)sampleCount * (sampleCount - 1);
  410.       double waitingTime = random.NextExponential(rate);
  411.       return currentTime + waitingTime;
  412.    }
  413.  
  414.    private void AddMutation(CoalescentTreeNode node, double theta)
  415.    {
  416.       int count = random.NextPoisson(node.WaitingTime * theta);
  417.       while (--count >= 0)
  418.       {
  419.          double site = random.NextDouble();
  420.          node.AddMutation(site);
  421.       }
  422.    }
  423. }
  424.  
  425. public static class RandomExtension
  426. {
  427.    public static double NextExponential(this Random random, double rate)
  428.    {
  429.       double r;
  430.       do
  431.       {
  432.          r = random.NextDouble();
  433.       } while (r == 0.0);
  434.       return -Math.Log(r) / rate;
  435.    }
  436.  
  437.    public static int NextPoisson(this Random random, double mean)
  438.    {
  439.       double lambda = Math.Exp(mean);
  440.       int k = 0;
  441.       if (!double.IsInfinity(lambda))
  442.       {
  443.          while ((lambda *= random.NextDouble()) > 1.0)
  444.          {
  445.             k++;
  446.          }
  447.       }
  448.       else
  449.       {
  450.          lambda = mean;
  451.          while ((lambda += Math.Log(random.NextDouble())) > 0.0)
  452.          {
  453.             k++;
  454.          }
  455.       }
  456.       return k; 
  457.    }
  458. }
  459.  
  460. class Program
  461. {
  462.    static void Main()
  463.    {
  464.       CoalescentSimulator simulator = new CoalescentSimulator();
  465.       var nodes = Enumerable.Range(0, 6).Select(i => new CoalescentTreeNode() { Name = (i + 1).ToString() }).ToArray();
  466.       CoalescentTree tree = simulator.PerformSimulation(nodes, 5.0);
  467.       Console.WriteLine(tree.ToString("0.000", null));
  468.       Console.WriteLine(tree.ToString("A", null));
  469.    }
  470. }
Success #stdin #stdout 0.03s 38176KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
((3:0.015,(6:0.011,5:0.011):0.003):0.545,(4:0.035,(1:0.019,2:0.019):0.017):0.524);
3: 1111100
6: 1111100
5: 1111100
4: 0000010
1: 0000010
2: 0000011
https://ideone.com/aHskQ
language:
C# (gmcs 5.20.1)
created:
13 years ago
visibility:
public

Share or Embed source code

Discover > Sphere Engine API

The brand new service which powers Ideone!

Discover > IDE Widget

Widget for compiling and running the source code in a web browser!