• API
• FAQ
• Tools
• Archive
SHARE
TWEET # Untitled a guest May 2nd, 2017 66 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
1. import java.util.Random;
3.
4. public class EColiSim {
5.
6.
7.     /* E. Coli K-12 has about 4.5 million base pair genome,
8.     and a mutation rate of about 1e-3 per new individual,
9.     for a total chance of about 2.22e-10 mutation rate per base-pair per individual.
10.
11.     Our sim is going to use a simplified genome: 5 base pairs total.
12.     The mutation rate is going to be about the same: 2.22e-10 mutation rate.
13.     */
14.     final static double mutationRate = Double.parseDouble("2.22e-10");
15.
16.     final static int[] fiveChooseX = new int[]{ 1, 5, 10, 10, 5, 1 };
17.
18.     final static double countPerMilliLiter = Double.parseDouble("1e+9");
19.     final static double olympicSwimmingPoolSizeInLiter = Double.parseDouble("2.5e+9");
20.     final static double initialPopulationCount = countPerMilliLiter * olympicSwimmingPoolSizeInLiter * 1000.0;
21.
22.     final static int numBases = 4; //A C G T
23.     final static int numBasePairs = 5;
24.     final static int numPopulations = 1024; //4 pow 5;
25.
26.     //When 1, math is fully accurate (more or less).
27.     //When greater than 1, it multiplies the mutation rate by this number, and the generation increment count is also this number.
28.     //Putting this number greater than 1 decreases sim time, at the cost of accuracy.
29.     //However, for our purposes here, it should /overestimate/ the time it takes to reach the solution,
30.     //because it doesn't compound the numbers properly.
31.     //Ex: there's a difference between an interest rate of 12% per year compounded annually, and 1% per month compounded monthly.
32.     //In particular, a positive interest return rate of 12% per year compounded annually is substantially less than 1% per month compounded monthly.
33.     final static int skippingFactor = 1000;
34.
35.     public static void main(final String[] args) {
36.
37.         //initialize first gen
38.         final Generation firstGeneration = new Generation();
39.         firstGeneration.setPopulationSize(0, initialPopulationCount);
40.         for (int ident = 1; ident < numPopulations; ++ident) {
41.             firstGeneration.setPopulationSize(ident, 0.0);
42.         }
43.
44.         //start the sim
45.         long prevGenerationCount = 0;
46.         Generation prevGeneration = firstGeneration;
47.         for (;;) {
48.             final Generation nextGeneration = new Generation();
49.             for (int nextIdent = 0; nextIdent < numPopulations; ++nextIdent ) {
50.                 double sum = 0.0;
51.                 for (int prevIdent = 0; prevIdent < numPopulations; ++prevIdent ) {
52.                     final int distance = getDistance(nextIdent, prevIdent);
53.                     final double chance = getChance(distance);
54.                     final double prevPopulationSize = prevGeneration.getPopulationSize(prevIdent);
55.                     double newPopulationContributingFactor = chance * prevPopulationSize;
56.                     newPopulationContributingFactor = randomlyRoundUpOrDown(newPopulationContributingFactor);
57.
58.                     sum += newPopulationContributingFactor;
59.                 }
60.                 nextGeneration.setPopulationSize(nextIdent, sum);
61.             }
62.
63.             //loop increment
64.             prevGeneration = nextGeneration;
65.             prevGenerationCount += skippingFactor;
66.
67. //            //for debugging: print every population value for this generation
68. //            System.out.println("Generation count " + prevGenerationCount);
69. //            for (int ident = 0; ident < numPopulations; ++ident ) {
70. //                System.out.println("" + Generation.getPopulationName(ident) + ": " + nextGeneration.getPopulationSize(ident));
71. //            }
72. //            System.out.println("----");
73.
74.             //for normal use, print summary information for this generation
75.             System.out.println(""
76.                     + prevGenerationCount + ", " + nextGeneration.getPopulationSize(0)
77.                     + " ... "
78.                     + prevGenerationCount + ", " + nextGeneration.getPopulationSize(0 + 1)
79.                     + " ... "
80.                     + prevGenerationCount + ", " + nextGeneration.getPopulationSize(0 + 1 + 4)
81.                     + " ... "
82.                     + prevGenerationCount + ", " + nextGeneration.getPopulationSize(0 + 1 + 4 + 16)
83.                     + " ... "
84.                     + prevGenerationCount + ", " + nextGeneration.getPopulationSize(0 + 1 + 4 + 16 + 64)
85.                     + " ... "
86.                     + prevGenerationCount + ", " + nextGeneration.getPopulationSize(numPopulations - 1)
87.                     );
88.
89.             //sanity
90.             double totalCount = getTotalCount(nextGeneration);
91.             if (totalCount * 1.1 < initialPopulationCount || totalCount / 1.1 > initialPopulationCount) {
92.                 throw new RuntimeException("the math of the sim is wrong, or errors compounded sufficiently to get huge errors");
93.             }
94.
95.             //loop break test
96.             if (nextGeneration.getPopulationSize(numPopulations - 1) >= 1.0) {
97.                 break;
98.             }
99.         }
100.         System.out.println("Generation count " + prevGenerationCount);
101.         for (int ident = 0; ident < numPopulations; ++ident ) {
102.             System.out.println("" + Generation.getPopulationName(ident) + ": " + prevGeneration.getPopulationSize(ident));
103.         }
104.         System.out.println("----");
105.     }
106.
107.     static double randomlyRoundUpOrDown(double x) {
108.         if (ThreadLocalRandom.current().nextDouble(0, 1) <= x - Math.floor(x)) {
109.             x = Math.floor(x) + 1;
110.         }else{
111.             x = Math.floor(x) ;
112.         }
113.         return x;
114.     }
115.
116.     static int getDistance(int i, int j) {
117.         int d = 0;
118.         for (int index = 0; index < numBasePairs; ++index) {
119.             int filter = (1+2) << (2*index);
120.             if ((i & filter) != (j & filter)) {
121.                 ++d;
122.             }
123.         }
124.         return d;
125.     }
126.
127.     static double getChance(int distance) {
128.         double mutationRateForTarget = mutationRate * (1.0 / (double)numBases); //There are 4 bases, ACGT, so only 1/4 for getting our target base
129.         mutationRateForTarget *= (double)skippingFactor;
130.         return ((double)fiveChooseX[distance])
131.                 * Math.pow(mutationRateForTarget, distance)
132.                 * Math.pow(1.0 - mutationRateForTarget, 5 - distance);
133.     }
134.
135.     static double getTotalCount(Generation gen) {
136.         double sum = 0.0;
137.         for (int ident = 0; ident < numPopulations; ++ident) {
138.             sum += gen.getPopulationSize(ident);
139.         }
140.         return sum;
141.     }
142.
143.     static final class Generation {
144.         public double getPopulationSize(int ident) {
145.             return populations[ident];
146.         }
147.         public void setPopulationSize(int ident, double val) {
148.             populations[ident] = val;
149.         }
150.         public static String getPopulationName(int ident) {
151.             final StringBuilder name = new StringBuilder();
152.             name.append("Population [");
153.             for (int index = 0; index < numBasePairs; ++index) {
154.                 int i = (ident >> (2 * index)) & (1+2);
155.                 name.append("" + i);
156.             }
157.             name.append("]");
158.             return name.toString();
159.         }
160.         private final double[] populations = new double[numPopulations];
161.     }
162.
163. }
RAW Paste Data
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.

Top