fork download 
  1. using System;
  2. using System.Collections.Generic;
  3. using System.Collections.ObjectModel;
  4. using System.Linq;
  5.  
  6. public class CoalescentTreeNode
  7. {
  8.    public string Name { get; set; }
  9.  
  10.    private CoalescentTree tree;
  11.    public CoalescentTree Tree
  12.    {
  13.       get { return this.tree; }
  14.       internal set { this.tree = value; }
  15.    }
  16.  
  17.    private CoalescentTreeNode left;
  18.    public CoalescentTreeNode DescendantLeft
  19.    {
  20.       get { return this.left; }
  21.       internal set { this.left = value; }
  22.    }
  23.  
  24.    private CoalescentTreeNode right;
  25.    public CoalescentTreeNode DescendantRight
  26.    {
  27.       get { return this.right; }
  28.       internal set { this.right = value; }
  29.    }
  30.  
  31.    private CoalescentTreeNode ancestor;
  32.    public CoalescentTreeNode Ancestor
  33.    {
  34.       get { return this.ancestor; }
  35.       internal set { this.ancestor = value; }
  36.    }
  37.  
  38.    private double time;
  39.    public double WaitingTime
  40.    {
  41.       get { return this.time; }
  42.       internal set { this.time = value; }
  43.    }
  44.  
  45.    public bool IsLeaf
  46.    {
  47.       get
  48.       {
  49.          return DescendantLeft == null && DescendantRight == null;
  50.       }
  51.    }
  52.  
  53.    internal bool HasDescendant
  54.    {
  55.       get
  56.       {
  57.          return DescendantLeft != null && DescendantRight != null;
  58.       }
  59.    }
  60.  
  61.    public CoalescentTreeNode()
  62.    {
  63.    }
  64.  
  65.    internal double GetTotalTime()
  66.    {
  67.       double sum = 0.0;
  68.       CoalescentTreeNode current = DescendantLeft;
  69.       while (current != null)
  70.       {
  71.          sum += current.WaitingTime;
  72.          current = current.DescendantLeft;
  73.       }
  74.       return sum;
  75.    }
  76. }
  77.  
  78. public class CoalescentTree : IFormattable
  79. {
  80.    private readonly List<CoalescentTreeNode> nodes;
  81.  
  82.    public CoalescentTreeNode Root
  83.    {
  84.       get { return IsIncomplete ? null : this.nodes[0]; }
  85.    }
  86.  
  87.    public bool IsIncomplete
  88.    {
  89.       get { return this.nodes.Count != 1; }
  90.    }
  91.  
  92.    public CoalescentTree()
  93.    {
  94.       this.nodes = new List<CoalescentTreeNode>();
  95.    }
  96.  
  97.    internal CoalescentTree(IEnumerable<CoalescentTreeNode> nodes)
  98.    {
  99.       this.nodes = new List<CoalescentTreeNode>(nodes);
  100.    }
  101.  
  102.    public CoalescentTreeNode Coalesce(CoalescentTreeNode left, CoalescentTreeNode right, double time)
  103.    {
  104.       if (!Contains(left) || left.Ancestor != null)
  105.       {
  106.          throw new ArgumentException("left");
  107.       }
  108.       if (!Contains(right) || right.Ancestor != null)
  109.       {
  110.          throw new ArgumentException("right");
  111.       }
  112.       double leftTime = left.GetTotalTime();
  113.       if (leftTime > time)
  114.       {
  115.          throw new ArgumentException("left");
  116.       }
  117.       double rightTime = right.GetTotalTime();
  118.       if (rightTime > time)
  119.       {
  120.          throw new ArgumentException("right");
  121.       }
  122.  
  123.       left.WaitingTime = time - leftTime;
  124.       right.WaitingTime = time - rightTime;
  125.  
  126.       CoalescentTreeNode ancestor = new CoalescentTreeNode()
  127.       {
  128.          DescendantLeft = left,
  129.          DescendantRight = right,
  130.       };
  131.       left.Ancestor = ancestor;
  132.       right.Ancestor = ancestor;
  133.  
  134.       Remove(left);
  135.       Remove(right);
  136.       Add(ancestor);
  137.       return ancestor;
  138.    }
  139.  
  140.    public void Add(CoalescentTreeNode node)
  141.    {
  142.       if (node == null)
  143.       {
  144.          throw new ArgumentNullException("node");
  145.       }
  146.       if (node.Tree != null)
  147.       {
  148.          throw new InvalidOperationException();
  149.       }
  150.       this.nodes.Add(node);
  151.       node.Tree = this;
  152.    }
  153.  
  154.    private bool Remove(CoalescentTreeNode node)
  155.    {
  156.       return this.nodes.Remove(node);
  157.    }
  158.  
  159.    private bool Contains(CoalescentTreeNode node)
  160.    {
  161.       return this.nodes.Contains(node);
  162.    }
  163.  
  164.    internal IList<CoalescentTreeNode> GetNodes()
  165.    {
  166.       return new ReadOnlyCollection<CoalescentTreeNode>(nodes);
  167.    }
  168.  
  169.    public override string ToString()
  170.    {
  171.       return ToString(null, null);
  172.    }
  173.  
  174.    public string ToString(string format, IFormatProvider formatProvider)
  175.    {
  176.       if (Root.HasDescendant)
  177.       {
  178.          return string.Format("({0},{1});", ToString(Root.DescendantLeft, format, formatProvider), ToString(Root.DescendantRight, format, formatProvider));
  179.       }
  180.       return string.Format("{0}:{1};", Root.Name, Root.WaitingTime);
  181.    }
  182.  
  183.    private string ToString(CoalescentTreeNode node, string format, IFormatProvider formatProvider)
  184.    {
  185.       if (node == null)
  186.       {
  187.          return string.Empty;
  188.       }
  189.       if (node.HasDescendant)
  190.       {
  191.          return string.Format("({0},{1}):{2}", ToString(node.DescendantLeft, format, formatProvider), ToString(node.DescendantRight, format, formatProvider), node.WaitingTime.ToString(format, formatProvider));
  192.  
  193.       }
  194.       return string.Format("{0}:{1}", node.Name, node.WaitingTime.ToString(format, formatProvider));
  195.    }
  196. }
  197.  
  198. public class CoalescentSimulator
  199. {
  200.    private readonly Random random;
  201.  
  202.    public CoalescentSimulator()
  203.    {
  204.       this.random = new Random();
  205.    }
  206.  
  207.    public CoalescentTree PerformSimulation(IEnumerable<CoalescentTreeNode> samples)
  208.    {
  209.       CoalescentTree tree = new CoalescentTree(samples);
  210.  
  211.       double time = 0.0;
  212.       while (tree.IsIncomplete)
  213.       {
  214.          IList<CoalescentTreeNode> nodes = tree.GetNodes();
  215.          int count = nodes.Count;
  216.          time = GenerateNextCoalescentEventTime(time, count);
  217.          int leftIndex;
  218.          int rightIndex;
  219.          GetRandomTwoIndexes(count, out leftIndex, out rightIndex);
  220.          tree.Coalesce(nodes[leftIndex], nodes[rightIndex], time);
  221.       }
  222.  
  223.       return tree;
  224.    }
  225.  
  226.    private void GetRandomTwoIndexes(int count, out int index1, out int index2)
  227.    {
  228.       index1 = random.Next(count);
  229.       index2 = random.Next(count - 1);
  230.       if (index2 >= index1)
  231.       {
  232.          ++index2;
  233.       }
  234.    }
  235.  
  236.    private double GenerateNextCoalescentEventTime(double currentTime, int sampleCount)
  237.    {
  238.       double rate = (double)sampleCount * (sampleCount - 1);
  239.       double waitingTime = random.NextExponential(rate);
  240.       return currentTime + waitingTime;
  241.    }
  242. }
  243.  
  244. public static class RandomExtension
  245. {
  246.    public static double NextExponential(this Random random, double rate)
  247.    {
  248.       double r;
  249.       do
  250.       {
  251.          r = random.NextDouble();
  252.       } while (r == 0.0);
  253.       return -Math.Log(r) / rate;
  254.    }
  255. }
  256.  
  257. class Program
  258. {
  259.    static void Main()
  260.    {
  261.       CoalescentSimulator simulator = new CoalescentSimulator();
  262.       for (int loop = 0; loop < 10; loop++)
  263.       {
  264.          var nodes = Enumerable.Range(0, 6).Select(i => new CoalescentTreeNode() { Name = (i + 1).ToString() }).ToArray();
  265.          CoalescentTree tree = simulator.PerformSimulation(nodes);
  266.          Console.WriteLine(tree.ToString("0.0000", null));
  267.       }
  268.    }
  269. }
Success #stdin #stdout 0.02s 38080KB
comments (?)
stdin 
Standard input is empty
stdout 
(4:0.5722,(5:0.3363,((6:0.0023,2:0.0023):0.3334,(3:0.2688,1:0.2688):0.0670):0.0005):0.2359);
((4:0.0970,6:0.0970):0.8588,(2:0.4736,(3:0.4464,(5:0.0237,1:0.0237):0.4226):0.0272):0.4823);
((((2:0.0222,3:0.0222):0.0224,5:0.0447):0.0245,1:0.0692):0.5238,(4:0.4575,6:0.4575):0.1355);
((3:0.1441,(1:0.0062,4:0.0062):0.1379):0.2228,((5:0.1537,2:0.1537):0.0695,6:0.2232):0.1437);
((5:0.3036,1:0.3036):0.0008,((3:0.0320,4:0.0320):0.2176,(6:0.0449,2:0.0449):0.2048):0.0547);
((5:0.1800,(3:0.0472,6:0.0472):0.1328):0.2682,(2:0.0253,(4:0.0160,1:0.0160):0.0092):0.4229);
((4:1.2542,((6:0.0677,2:0.0677):0.2818,5:0.3495):0.9047):0.1914,(1:0.0146,3:0.0146):1.4310);
(((4:0.1030,1:0.1030):0.0090,3:0.1119):1.4891,(2:0.7473,(5:0.0837,6:0.0837):0.6636):0.8537);
((2:0.0430,3:0.0430):0.5369,((1:0.0965,6:0.0965):0.1061,(5:0.0570,4:0.0570):0.1456):0.3773);
(((1:0.0390,(6:0.0021,3:0.0021):0.0369):0.2990,(5:0.1011,2:0.1011):0.2369):0.3773,4:0.7153);
https://ideone.com/t1fkB
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!