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;
  2. import java.util.concurrent.ThreadLocalRandom;
  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. OK, I Understand
 
Top