import java.util.Random;
import java.util.concurrent.ThreadLocalRandom;

public class EColiSim {


    /* E. Coli K-12 has about 4.5 million base pair genome,
    and a mutation rate of about 1e-3 per new individual,
    for a total chance of about 2.22e-10 mutation rate per base-pair per individual.

    Our sim is going to use a simplified genome: 5 base pairs total.
    The mutation rate is going to be about the same: 2.22e-10 mutation rate.
    */
    final static double mutationRate = Double.parseDouble("2.22e-10");

    final static int[] fiveChooseX = new int[]{ 1, 5, 10, 10, 5, 1 };

    final static double countPerMilliLiter = Double.parseDouble("1e+9");
    final static double olympicSwimmingPoolSizeInLiter = Double.parseDouble("2.5e+9");
    final static double initialPopulationCount = countPerMilliLiter * olympicSwimmingPoolSizeInLiter * 1000.0;

    final static int numBases = 4; //A C G T
    final static int numBasePairs = 5;
    final static int numPopulations = 1024; //4 pow 5;

    //When 1, math is fully accurate (more or less).
    //When greater than 1, it multiplies the mutation rate by this number, and the generation increment count is also this number.
    //Putting this number greater than 1 decreases sim time, at the cost of accuracy.
    //However, for our purposes here, it should /overestimate/ the time it takes to reach the solution,
    //because it doesn't compound the numbers properly.
    //Ex: there's a difference between an interest rate of 12% per year compounded annually, and 1% per month compounded monthly.
    //In particular, a positive interest return rate of 12% per year compounded annually is substantially less than 1% per month compounded monthly.
    final static int skippingFactor = 1000;

    public static void main(final String[] args) {

        //initialize first gen
        final Generation firstGeneration = new Generation();
        firstGeneration.setPopulationSize(0, initialPopulationCount);
        for (int ident = 1; ident < numPopulations; ++ident) {
            firstGeneration.setPopulationSize(ident, 0.0);
        }

        //start the sim
        long prevGenerationCount = 0;
        Generation prevGeneration = firstGeneration;
        for (;;) {
            final Generation nextGeneration = new Generation();
            for (int nextIdent = 0; nextIdent < numPopulations; ++nextIdent ) {
                double sum = 0.0;
                for (int prevIdent = 0; prevIdent < numPopulations; ++prevIdent ) {
                    final int distance = getDistance(nextIdent, prevIdent);
                    final double chance = getChance(distance);
                    final double prevPopulationSize = prevGeneration.getPopulationSize(prevIdent);
                    double newPopulationContributingFactor = chance * prevPopulationSize;
                    newPopulationContributingFactor = randomlyRoundUpOrDown(newPopulationContributingFactor);

                    sum += newPopulationContributingFactor;
                }
                nextGeneration.setPopulationSize(nextIdent, sum);
            }

            //loop increment
            prevGeneration = nextGeneration;
            prevGenerationCount += skippingFactor;

//            //for debugging: print every population value for this generation
//            System.out.println("Generation count " + prevGenerationCount);
//            for (int ident = 0; ident < numPopulations; ++ident ) {
//                System.out.println("" + Generation.getPopulationName(ident) + ": " + nextGeneration.getPopulationSize(ident));
//            }
//            System.out.println("----");

            //for normal use, print summary information for this generation
            System.out.println(""
                    + prevGenerationCount + ", " + nextGeneration.getPopulationSize(0)
                    + " ... "
                    + prevGenerationCount + ", " + nextGeneration.getPopulationSize(0 + 1)
                    + " ... "
                    + prevGenerationCount + ", " + nextGeneration.getPopulationSize(0 + 1 + 4)
                    + " ... "
                    + prevGenerationCount + ", " + nextGeneration.getPopulationSize(0 + 1 + 4 + 16)
                    + " ... "
                    + prevGenerationCount + ", " + nextGeneration.getPopulationSize(0 + 1 + 4 + 16 + 64)
                    + " ... "
                    + prevGenerationCount + ", " + nextGeneration.getPopulationSize(numPopulations - 1)
                    );

            //sanity
            double totalCount = getTotalCount(nextGeneration);
            if (totalCount * 1.1 < initialPopulationCount || totalCount / 1.1 > initialPopulationCount) {
                throw new RuntimeException("the math of the sim is wrong, or errors compounded sufficiently to get huge errors");
            }

            //loop break test
            if (nextGeneration.getPopulationSize(numPopulations - 1) >= 1.0) {
                break;
            }
        }
        System.out.println("Generation count " + prevGenerationCount);
        for (int ident = 0; ident < numPopulations; ++ident ) {
            System.out.println("" + Generation.getPopulationName(ident) + ": " + prevGeneration.getPopulationSize(ident));
        }
        System.out.println("----");
    }

    static double randomlyRoundUpOrDown(double x) {
        if (ThreadLocalRandom.current().nextDouble(0, 1) <= x - Math.floor(x)) {
            x = Math.floor(x) + 1;
        }else{
            x = Math.floor(x) ;
        }
        return x;
    }

    static int getDistance(int i, int j) {
        int d = 0;
        for (int index = 0; index < numBasePairs; ++index) {
            int filter = (1+2) << (2*index);
            if ((i & filter) != (j & filter)) {
                ++d;
            }
        }
        return d;
    }

    static double getChance(int distance) {
        double mutationRateForTarget = mutationRate * (1.0 / (double)numBases); //There are 4 bases, ACGT, so only 1/4 for getting our target base
        mutationRateForTarget *= (double)skippingFactor;
        return ((double)fiveChooseX[distance])
                * Math.pow(mutationRateForTarget, distance)
                * Math.pow(1.0 - mutationRateForTarget, 5 - distance);
    }

    static double getTotalCount(Generation gen) {
        double sum = 0.0;
        for (int ident = 0; ident < numPopulations; ++ident) {
            sum += gen.getPopulationSize(ident);
        }
        return sum;
    }

    static final class Generation {
        public double getPopulationSize(int ident) {
            return populations[ident];
        }
        public void setPopulationSize(int ident, double val) {
            populations[ident] = val;
        }
        public static String getPopulationName(int ident) {
            final StringBuilder name = new StringBuilder();
            name.append("Population [");
            for (int index = 0; index < numBasePairs; ++index) {
                int i = (ident >> (2 * index)) & (1+2);
                name.append("" + i);
            }
            name.append("]");
            return name.toString();
        }
        private final double[] populations = new double[numPopulations];
    }

}