##################################
#                                #
# Last modified 2017/11/28       # 
#                                #
# Georgi Marinov                 #
#                                # 
##################################

import sys
import string
import scipy.stats
import numpy
import math
import random
from sets import Set
import time

def run():

    if len(sys.argv) < 8:
        print 'usage: python %s genome_length GC_fraction AT->AT AT->GC GC->AT GC->GC N_Iterations outfilename' % sys.argv[0]
        sys.exit(1)

    GL = int(sys.argv[1])
    GC = float(sys.argv[2])
    ATtoAT = int(sys.argv[3])
    ATtoGC = int(sys.argv[4])
    GCtoAT = int(sys.argv[5])
    GCtoGC = int(sys.argv[6])
    N = int(sys.argv[7])
    outfilename = sys.argv[8]

    GC = int(GC*GL)
    AT = int(GL-GC)

    p_mut_AT = (ATtoAT + ATtoGC + 0.0)/AT
    p_mut_GC = (GCtoAT + GCtoGC + 0.0)/GC

    p_mut_ATtoAT = ATtoAT/(ATtoAT + ATtoGC + 0.0)
    p_mut_ATtoGC = ATtoGC/(ATtoAT + ATtoGC + 0.0)

    p_mut_GCtoAT = GCtoAT/(GCtoAT + GCtoGC + 0.0)
    p_mut_GCtoGC = GCtoGC/(GCtoAT + GCtoGC + 0.0)

    print p_mut_AT, p_mut_GC, p_mut_ATtoAT, p_mut_ATtoGC, p_mut_GCtoAT, p_mut_GCtoGC

    outfile = open(outfilename, 'w')

    outline = '#Iteration\tAT->AT\tAT->GC\tGC->AT\tGC->GC\tu/(u+v)\tv/u\t'
    outfile.write(outline + '\n')

    for i in range(N):
        if i % 1000 == 0:
            print i
        simATtoAT = 0.0
        simATtoGC = 0.0
        simGCtoAT = 0.0
        simGCtoGC = 0.0
        for j in range(AT):
            PM = random.random()
            if PM < p_mut_AT:
                PMAT = random.random()
                if PMAT < p_mut_ATtoAT:
                    simATtoAT += 1
                else:
                    simATtoGC += 1
        for j in range(GC):
            PM = random.random()
            if PM < p_mut_GC:
                PMGC = random.random()
                if PMGC < p_mut_GCtoAT:
                    simGCtoAT += 1
                else:
                    simGCtoGC += 1
#        print simATtoAT, simATtoGC, simGCtoAT, simGCtoGC
        u = simATtoGC/AT
        v = simGCtoAT/GC
        GCeq = u/(u+v)
        if u == 0:
            MutBias = 'nan'
        else:
            MutBias = v/u
        outline = str(i+1) + '\t' + str(simATtoAT) + '\t' + str(simATtoGC) + '\t' + str(simGCtoAT) + '\t' + str(simGCtoGC) + '\t' + str(GCeq) + '\t' + str(MutBias)
        outfile.write(outline + '\n')

    outfile.close()

run()