##########----------Import packages
import os, sys, time, argparse
import numpy as np
import pandas as pd
import random
from multiprocessing import Pool

##########----------calc fractional occupany for haploid
def calcFO_Haploid(conc1, m1, k):
    numer = conc1 / (k ** m1)
    denom = 1 + conc1 / (k ** m1)
    return (numer / denom)

#########----------calc fractional occupancy
def calcFO(conc1, conc2, m1, m2, k):
    numer = conc1 / (k ** m1)
    denom = 1 + conc1 / (k ** m1) + conc2 / (k ** m2)
    return (numer / denom)

##########----------calc expression of B based on allele at A
##########  if A is 0 it reduce expression, if A is 1 it enhances expression
def calcExpB(A, conc1, conc2, m1, m2, k):
    FO = calcFO(conc1, conc2, m1, m2, k)
    if A==0:
        return (-1 * FO)
    else:
        return (FO)

#########----------return gt given two alleles
#########   assume values are 0 or 1 only
def calcGt(allele1, allele2):
        if allele1 == 0 and allele2 == 0:
            return 0
        elif allele1 == 1 and allele2 == 1:
            return 2
        else:
            return 1

##########----------return binding affinity for alpha based on sex
def calcMAlpha(SEX, ALPHA):
    if SEX=='m':
        return ALPHA
    else:
        return round(1-ALPHA, 1)

##########----------return binding affinity for beta based allele of A
def calcMBeta(A, BETA):
    if A==0:
        return BETA
    else:
        return round(1-BETA, 1)

##########----------calc absolute fitness based on sex, standardize concentration of B and selection regime
##########  SR = 0 for sexual antagonism
##########  SR = 1 for sexually concordant selection using 'female' fitness eq
##########  SR = 2 for sexually concordant selection using 'male' fitness eq

def calcAbsFitness(SEX, STCONC, SR, S_F, S_M, GAMMA_F, GAMMA_M):
    if SR == 0:
        if SEX == 'm':
            return (1 - (S_M * STCONC) ** GAMMA_M)
        else:
            return (1 - (S_F * (1 - STCONC) ** GAMMA_F))
    elif SR == 1:
        return (1 - (S_F * (1 - STCONC) ** GAMMA_F))
    elif SR == 2:
        return (1 - (S_M * STCONC) ** GAMMA_M)
'''
def calcAbsFitness(SEX, STCONC, SR, S_F, S_M, GAMMA_F, GAMMA_M):
    if SR == 0:
        if SEX == 'm':
            return 0.9
        else:
            return 0.1
    elif SR == 1:
        return 0.25
    elif SR == 2:
        return 0.75
'''
##########----------mutate binding values (alpha, beta)
##########   increments of 0.1 and value bound between 0 and 1
def mutateBindValue(mu, bValue):
    rBinom = np.random.binomial(1, mu)
    if rBinom == 1:
        newBValue = -1
        while (newBValue < 0 or newBValue > 1):
            epsilon = random.choices([-0.1, 0.1], k=1)
            newBValue = round(bValue + epsilon[0], 1)
    else:
        newBValue = bValue
    return (newBValue)

###########----------generate gamete through selection on parents
def selectionGamete(listIndiv, maxFit, popSize):
    #####-----while loop stop when get enough gametes
    listID, nGam = [], 0
    lGamAlpha, lGamA, lGamBeta = [], [], []
    while nGam < popSize:
        #####-----choose random parent and get relative fitness
        pID = int(random.randint(0, len(listIndiv) - 1))
        pFit = listIndiv[pID].getAbsFitness()
        relFit = pFit/maxFit
        #####-----selection by testing relFit against random uniform variate
        rVar = np.random.uniform(0, 1, 1)
        if relFit >= rVar:
            nGam += 1
            pGam = listIndiv[pID].getGamete()
            lGamAlpha.append(pGam[0])
            lGamA.append(pGam[1])
            lGamBeta.append(pGam[2])
    #####-----return list of gametes
    return ([lGamAlpha, lGamA, lGamBeta])

##########----------get highest fitness
def getMaxAbsFit(listIndiv):
    listAbsFit = [x.getAbsFitness() for x in listIndiv]
    maxAbsFit = np.max(listAbsFit)
    return(maxAbsFit)

##########----------get average fitness
def getAvgAbsFit(listIndiv):
    listStat = [x.getAbsFitness() for x in listIndiv]
    avgStat = np.mean(listStat)
    return (avgStat)

##########----------get average fitness by genotype of protein A
def getAvgAbsFit_GtA(listIndiv):
    listFit = [x.getAbsFitness() for x in listIndiv]
    listGt = [x.getGtA() for x in listIndiv]
    n0, n1, n2, fitA0, fitA1, fitA2 = 0, 0, 0, [], [], []
    for i in range(len(listIndiv)):
        gt, fit = listGt[i], listFit[i]
        if gt == 0:
            n0 += 1
            fitA0 = fitA0 + [fit]
        elif gt == 1:
            n1 += 1
            fitA1 = fitA1 + [fit]
        elif gt == 2:
            n2 += 1
            fitA2 = fitA2 + [fit]
    avgFitA0 = np.mean(fitA0) if n0 > 0 else 0
    avgFitA1 = np.mean(fitA1) if n1 > 0 else 0
    avgFitA2 = np.mean(fitA2) if n2 > 0 else 0
    return [avgFitA0, avgFitA1, avgFitA2]

##########----------get mean for total expression of A
def getAvgExpA(listIndiv):
    listStat = [sum(x.getExpA()) for x in listIndiv]
    avgStat = np.mean(listStat)
    return (avgStat)

##########----------get variance for total expression of A
def getVarExpA(listIndiv):
    listStat = [sum(x.getExpA()) for x in listIndiv]
    varStat = np.var(listStat)
    return (varStat)

##########----------get expression of allele A0 and A1 from heterozygotes averaged across individuals
def getAvgAlleleExpA_Het(listIndiv):
    hetConcA0, hetConcA1 = [], []
    for INDIV in listIndiv:
        gt = INDIV.getGtA()
        a1, a2 = INDIV.getA()
        concA1, concA2 = INDIV.getExpA()
        #####-----get conc of allele A0 and A1 in heterozygote
        if gt == 1:
            if a1 == 0:
                hetConcA0 += [concA1]
                hetConcA1 += [concA2]
            else:
                hetConcA0 += [concA2]
                hetConcA1 += [concA1]
    if len(hetConcA0) > 0:
        avgConcC0, avgConcA1 = np.mean(hetConcA0), np.mean(hetConcA1)
        return ([avgConcC0, avgConcA1])
    else:
        return ([-1, -1])

##########----------get average standardized experssion of B
def getAvgStExpB(listIndiv):
    listStat = [x.getStExpB() for x in listIndiv]
    avgStat = np.mean(listStat)
    return (avgStat)

##########----------get variance for standardized experssion of B
def getVarStExpB(listIndiv):
    listStat = [x.getStExpB() for x in listIndiv]
    varStat = np.var(listStat)
    return (varStat)


##########----------get frequency of A=1 allele
def getFreqA(listIndiv):
    listStat = [x.getGtA() for x in listIndiv]
    freqA = np.sum(listStat)/(2*len(listIndiv))
    return (freqA)

##########----------get genotype frequency at coding region of gene A
def getGtFreqA(listIndiv):
    listStat = [x.getGtA() for x in listIndiv]
    gtFreqA = [listStat.count(x)/len(listIndiv) for x in [0,1,2]]
    return (gtFreqA)

##########----------get average value of ALPHA
def getAvgAlpha(listIndiv):
    listStat = []
    [listStat.extend(x.getAlpha()) for x in listIndiv]
    avgStat = np.mean(listStat)
    return (avgStat)

##########----------get average value of BETA
def getAvgBeta(listIndiv):
    listStat = []
    [listStat.extend(x.getBeta()) for x in listIndiv]
    avgStat = np.mean(listStat)
    return (avgStat)

#########----------count occurence of each binding value
def countBindValue(listStat):
    lValues = [0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1]
    lCounts = [listStat.count(x) for x in lValues]
    return (lCounts)

#########----------count occurence of values for ALPHA
def getDistAlpha(listIndiv):
    listStat = []
    [listStat.extend(x.getAlpha()) for x in listIndiv]
    listCounts = countBindValue(listStat)
    return (listCounts)

#########----------count occurence of values for BETA
def getDistBeta(listIndiv):
    listStat = []
    [listStat.extend(x.getBeta()) for x in listIndiv]
    listCounts = countBindValue(listStat)
    return (listCounts)


##########----------count num individuals with exp B in discrete bins
##########  bins are MIN <= x < MAX
##########      last bin uses 0.9 <= x < 1.1 because maximum expression of B is 1
def countBinnedExpB(listStat):
    lMin = [0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9]
    lMax = [0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.1]
    lCounts = []
    for MIN, MAX in zip(lMin, lMax):
        lCounts += [sum([1 for x in listStat if MIN <= x < MAX])]
    return (lCounts)

def getDistStExpB(listIndiv):
    listStat = [x.getStExpB() for x in listIndiv]
    listCounts = countBinnedExpB(listStat)
    return (listCounts)

##########----------calculate correlation between loci
def getCorrel(listIndiv):
    #####get lists of value from each indiv
    N = len(listIndiv)
    alpha1, alpha2, a1, a2, beta1, beta2 = [-1]*N, [-1]*N, [-1]*N, [-1]*N, [-1]*N, [-1]*N
    for i, x in enumerate(listIndiv):
        alpha1[i], alpha2[i] = x.getAlpha()
        a1[i], a2[i] = x.getA()
        beta1[i], beta2[i] = x.getBeta()
    #####combine list between homologous Cm for haplotype correlations
    alphaAll, aAll, betaAll = alpha1 + alpha2, a1 + a2, beta1 + beta2
    #####haplotype correlations (across loci)
    cAlphaA, cAlphaBeta, cABeta = float("nan"), float("nan"), float("nan")
    if np.var(alphaAll) > 0 and np.var(aAll) > 0:
        cAlphaA = np.corrcoef(alphaAll, aAll)[0, 1]
    if np.var(alphaAll) > 0 and np.var(betaAll) > 0:
        cAlphaBeta = np.corrcoef(alphaAll, betaAll)[0, 1]
    if np.var(aAll) > 0 and np.var(betaAll) > 0:
        cABeta = np.corrcoef(aAll, betaAll)[0, 1]
    #####return
    return [cAlphaA, cAlphaBeta, cABeta]



####################--------------------Object to represent individual
class indiv():
    #####-----Initizlize individual
    def __init__(self, SEX, CONC_D, GA, ALPHA1, ALPHA2, A1, A2, BETA1, BETA2, KA, KB, ASAT, S_F, S_M, GAMMA_F, GAMMA_M, UA, UB):
        
        #####-----assign sex and allele
        #####   'alpha' is binding site allele for gene A (real number between 0 and 1)
        #####   'a' is protein coding allele for gene A (binary 0 or 1)
        #####   'beta' is binding site allele for gene B (real number between 0 and 1)
        self.sex = SEX
        self.uAlpha, self.uBeta = UA, UB
        self.alpha1, self.alpha2 = ALPHA1, ALPHA2
        self.a1, self.a2 = A1, A2
        self.beta1, self.beta2 = BETA1, BETA2
        self.gtA = calcGt(self.a1, self.a2)
        self.absFitness = 1

        #####-----assign binding affinity between sex TF and gene A binding site
        self.mSexAlpha1 = calcMAlpha(self.sex, self.alpha1)
        self.mSexAlpha2 = calcMAlpha(self.sex, self.alpha2)

        #####-----assign binding affinity between gene A product and gene B binding site
        #####   depends on value of 'a' (0 or 1) and 'beta'
        self.mA1Beta1 = calcMBeta(self.a1, self.beta1)
        self.mA1Beta2 = calcMBeta(self.a1, self.beta2)
        self.mA2Beta1 = calcMBeta(self.a2, self.beta1)
        self.mA2Beta2 = calcMBeta(self.a2, self.beta2)

        #####-----calc concentration of product from gene A (TF for gene B)
        self.concA1 = calcFO_Haploid(CONC_D, self.mSexAlpha1, KA)*GA
        self.concA2 = calcFO_Haploid(CONC_D, self.mSexAlpha2, KA)*GA

        #####-----calc expression of gene B
        #####   depends on concentration of gene A products (concA1, concA2) ,their identity (a1, a2) and binding sites affinity (mABeta)
        self.concB1 = 0.5 * (calcExpB(self.a1, self.concA1, self.concA2, self.mA1Beta1, self.mA2Beta1, KB) + calcExpB(self.a2, self.concA2, self.concA1, self.mA2Beta1, self.mA1Beta1, KB))
        self.concB2 = 0.5 * (calcExpB(self.a1, self.concA1, self.concA2, self.mA1Beta2, self.mA2Beta2, KB) + calcExpB(self.a2, self.concA2, self.concA1, self.mA2Beta2, self.mA1Beta2, KB))
        self.concB = self.concB1 + self.concB2

        #####-----calc max and min expression of B using saturating conc of A
        self.concMinB = -1 * calcFO_Haploid(ASAT, 0, KB)
        self.concMaxB = calcFO_Haploid(ASAT, 0, KB)

        #####-----standardize expression of B
        self.stConcB = (self.concB-self.concMinB)/( self.concMaxB-self.concMinB)

        #####-----calc absolute fitness using stConcB and selection parameters
        #self.absFitness = calcAbsFitness(SEX, self.stConcB, S_F, S_M, GAMMA_F, GAMMA_M)

    #####-----Propoerties of individual
    def getSex(self):
        return self.sex
    def getGtA(self):
        return self.gtA
    def getA(self):
        return [self.a1, self.a2]
    def getAlpha(self):
        return [self.alpha1, self.alpha2]
    def getBeta(self):
        return [self.beta1, self.beta2]
    def getMSAlpha(self):
        return [self.mSexAlpha1, self.mSexAlpha2]
    def getMABeta(self):
        return [self.mA1Beta1, self.mA1Beta2, self.mA2Beta1, self.mA2Beta2]
    def getExpA(self):
        return [self.concA1, self.concA2]
    def getExpB(self):
        return [self.concB1, self.concB2, self.concB]
    def getSatExpB(self):
        return [self.concMinB, self.concMaxB]
    def getStExpB(self):
        return self.stConcB
    def getAbsFitness(self):
        return self.absFitness
    def getSmy(self):
        SMY = [self.sex, self.gtA, self.a1, self.a2, self.alpha1, self.alpha2, self.beta1, self.beta2, self.concA1, self.concA2, self.stConcB, self.absFitness]
        return SMY

    #####-----randomly pass on one allele as gamete
    #####   assume complete linkage for gene A and free recomb between A and B
    def getGamete(self):
        #####-----get random 0 or 1
        rBinom1 = np.random.binomial(1, 0.5)
        rBinom2 = np.random.binomial(1, 0.5)
        #####-----rBinom1, rBinom2 controls inheritance of gene A and B, respectively
        if rBinom1 == 0:
            if rBinom2 == 0:
                return [mutateBindValue(self.uAlpha, self.alpha1), self.a1, mutateBindValue(self.uBeta, self.beta1)]
            else:
                return [mutateBindValue(self.uAlpha, self.alpha1), self.a1, mutateBindValue(self.uBeta, self.beta2)]
        else:
            if rBinom2 == 0:
                return [mutateBindValue(self.uAlpha, self.alpha2), self.a2, mutateBindValue(self.uBeta, self.beta1)]
            else:
                return [mutateBindValue(self.uAlpha, self.alpha2), self.a2, mutateBindValue(self.uBeta, self.beta2)]


####################--------------------Simulation type 1
####################  random assigment of ALPHA, A and BETA values
def runSim01(CONC_D, GA, KA, KB, ASAT, SR, S_F, S_M, GAMMA_F, GAMMA_M, UA, UB, N, PROPF, NUMGEN):
    ##########----------initialize dfs to collect data
    SIMTYPE = 1
    COLNAME = [
        'SIMTYPE','CONC_D','GA','KA','KB','ASAT','SR','S_F','S_M','GAMMA_F','GAMMA_M','UA','UB','N','PROPF','NUMGEN','GEN',
        'maxFit_F','maxFit_M','maxFit','avgFit_F','avgFit_M','avgFit',
        'avgFit_F00','avgFit_F01','avgFit_F11','avgFit_M00','avgFit_M01','avgFit_M11','avgFit_00','avgFit_01','avgFit_11',
        'avgExpA_F', 'avgExpA_M', 'avgExpA', 'varExpA_F', 'varExpA_M', 'varExpA',
        'avgStExpB_F', 'avgStExpB_M', 'avgStExpB', 'varStExpB_F', 'varStExpB_M', 'varStExpB',
        'avgHetExpA0_F', 'avgHetExpA1_F', 'avgHetExpA0_M', 'avgHetExpA1_M',
        'freqA_F','freqA_M','freqA','gtA00_F','gtA01_F','gtA11_F','gtA00_M','gtA01_M','gtA11_M','gtA00','gtA01','gtA11',
        'avgAlpha_F','avgAlpha_M','avgAlpha','avgBeta_F','avgBeta_M','avgBeta',
        'cAlphaA_F', 'cAlphaBeta_F', 'cABeta_F', 'cAlphaA_M', 'cAlphaBeta_M', 'cABeta_M'
    ]
    dfGeneral = pd.DataFrame(columns=COLNAME)

    COLNAME = [
        'SIMTYPE', 'CONC_D', 'GA', 'KA', 'KB', 'ASAT', 'SR', 'S_F', 'S_M', 'GAMMA_F', 'GAMMA_M', 'UA', 'UB', 'N', 'PROPF', 'NUMGEN', 'GEN',
        'alF0', 'alF1', 'alF2', 'alF3', 'alF4', 'alF5', 'alF6', 'alF7', 'alF8', 'alF9', 'alF10',
        'alM0', 'alM1', 'alM2', 'alM3', 'alM4', 'alM5', 'alM6', 'alM7', 'alM8', 'alM9', 'alM10',
        'beF0', 'beF1', 'beF2', 'beF3', 'beF4', 'beF5', 'beF6', 'beF7', 'beF8', 'beF9', 'beF10',
        'beM0', 'beM1', 'beM2', 'beM3', 'beM4', 'beM5', 'beM6', 'beM7', 'beM8', 'beM9', 'beM10',
        'al0', 'al1', 'al2', 'al3', 'al4', 'al5', 'al6', 'al7', 'al8', 'al9', 'al10',
        'be0', 'be1', 'be2', 'be3', 'be4', 'be5', 'be6', 'be7', 'be8', 'be9', 'be10',
        'BF0', 'BF1', 'BF2', 'BF3', 'BF4', 'BF5', 'BF6', 'BF7', 'BF8', 'BF9',
        'BM0', 'BM1', 'BM2', 'BM3', 'BM4', 'BM5', 'BM6', 'BM7', 'BM8', 'BM9',
    ]
    dfDist = pd.DataFrame(columns=COLNAME)

    ##########----------calculate num of male and females
    N_F= int(N*(1-PROPF))
    N_M = N - N_F

    ##########----------Initialize alleles
    ##########  equal prob for each allele of each locus
    listBind = [0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0]
    #listBind = [1.0]
    lAlpha1 = random.choices(listBind, k=N)
    lAlpha2 = random.choices(listBind, k=N)
    lA1 = np.random.binomial(1, 0.5, N)
    lA2 = np.random.binomial(1, 0.5, N)
    lBeta1 = random.choices(listBind, k=N)
    lBeta2 = random.choices(listBind, k=N)
    random.shuffle(lAlpha1)
    random.shuffle(lAlpha2)
    random.shuffle(lA1)
    random.shuffle(lA2)
    random.shuffle(lBeta1)
    random.shuffle(lBeta2)

    ##########----------Initialize sexes based on N_F, N_M
    popSex = ['f'] * int(N_F) + ['m'] * int(N_M)
    random.shuffle(popSex)

    ##########----------Initilize N individuals in population
    POPUL_F, POPUL_M = [], []
    for i in range(N):
        SEX = popSex[i]
        if SEX == 'f':
            POPUL_F.append(
                indiv('f', CONC_D, GA, lAlpha1[i], lAlpha2[i], lA1[i], lA2[i], lBeta1[i], lBeta2[i], KA, KB, ASAT, S_F,
                      S_M, GAMMA_F, GAMMA_M, UA, UB))
        else:
            POPUL_M.append(
                indiv('m', CONC_D, GA, lAlpha1[i], lAlpha2[i], lA1[i], lA2[i], lBeta1[i], lBeta2[i], KA, KB, ASAT, S_F,
                      S_M, GAMMA_F, GAMMA_M, UA, UB))
    ##########----------Initial calculation of absolute fitness
    for INDIV in POPUL_F + POPUL_M:
        INDIV.absFitness = calcAbsFitness(INDIV.getSex(), INDIV.getStExpB(), SR, S_F, S_M, GAMMA_F, GAMMA_M)

    ##########----------Run for NUMGEN generations
    ##########  if 'A' gets fixed (freq = 0 or 1), then sim stops after MaxExtraGEN generations
    extraGen, MaxExtraGen = 0, 500
    for gen in range(NUMGEN):
        #####-----determine highest fitness
        maxAbsFit_F = getMaxAbsFit(POPUL_F)
        maxAbsFit_M = getMaxAbsFit(POPUL_M)
        #####-----get gametes using relative fitness
        listFemGam = selectionGamete(POPUL_F, maxAbsFit_F, N)
        listMalGam = selectionGamete(POPUL_M, maxAbsFit_M, N)
        lAlpha1, lA1, lBeta1 = listFemGam[0], listFemGam[1], listFemGam[2]
        lAlpha2, lA2, lBeta2 = listMalGam[0], listMalGam[1], listMalGam[2]
        #####-----create offspring from gametes
        offspringSex = ['f'] * int(N_F) + ['m'] * int(N_M)
        random.shuffle(offspringSex)
        OFFSPRING_F, OFFSPRING_M = [], []
        for i in range(N):
            SEX = offspringSex[i]
            if SEX == 'f':
                tempInd = indiv('f', CONC_D, GA, lAlpha1[i], lAlpha2[i], lA1[i], lA2[i], lBeta1[i], lBeta2[i], KA, KB,
                                ASAT, S_F, S_M, GAMMA_F, GAMMA_M, UA, UB)
                OFFSPRING_F.append(tempInd)
            else:
                tempInd = indiv('m', CONC_D, GA, lAlpha1[i], lAlpha2[i], lA1[i], lA2[i], lBeta1[i], lBeta2[i], KA, KB,
                                ASAT, S_F, S_M, GAMMA_F, GAMMA_M, UA, UB)
                OFFSPRING_M.append(tempInd)
        #####-----new population created from offspring
        POPUL_F = OFFSPRING_F.copy()
        POPUL_M = OFFSPRING_M.copy()

        #####-----Calculate fitness
        for INDIV in POPUL_F + POPUL_M:
            INDIV.absFitness = calcAbsFitness(INDIV.getSex(), INDIV.getStExpB(), SR, S_F, S_M, GAMMA_F, GAMMA_M)

        #####-----Collect data every 10 generations
        if (gen + 1) % 10 == 0:
            vParam = [SIMTYPE, CONC_D, GA, KA, KB, ASAT, SR, S_F, S_M, GAMMA_F, GAMMA_M, UA, UB, N, PROPF, NUMGEN, gen+1]
            vMaxFit = [getMaxAbsFit(POPUL_F), getMaxAbsFit(POPUL_M), getMaxAbsFit(POPUL_F+POPUL_M)]
            vAvgFit = [getAvgAbsFit(POPUL_F), getAvgAbsFit(POPUL_M), getAvgAbsFit(POPUL_F+POPUL_M)]
            vAvgFit_GtA = getAvgAbsFit_GtA(POPUL_F) + getAvgAbsFit_GtA(POPUL_M) + getAvgAbsFit_GtA(POPUL_F+POPUL_M)
            vAvgStExpA = [getAvgExpA(POPUL_F), getAvgExpA(POPUL_M), getAvgExpA(POPUL_F + POPUL_M)]
            vVarStExpA = [getVarExpA(POPUL_F), getVarExpA(POPUL_M), getVarExpA(POPUL_F + POPUL_M)]
            vAvgStExpB = [getAvgStExpB(POPUL_F), getAvgStExpB(POPUL_M), getAvgStExpB(POPUL_F+POPUL_M)]
            vVarStExpB = [getVarStExpB(POPUL_F), getVarStExpB(POPUL_M), getVarStExpB(POPUL_F+POPUL_M)]
            vAvgAlleleExpA_het = getAvgAlleleExpA_Het(POPUL_F) + getAvgAlleleExpA_Het(POPUL_M)
            vFreqA = [getFreqA(POPUL_F), getFreqA(POPUL_M), getFreqA(POPUL_F+POPUL_M)]
            vGtA = getGtFreqA(POPUL_F) + getGtFreqA(POPUL_M) + getGtFreqA(POPUL_F+POPUL_M)
            vAvgAlpha = [getAvgAlpha(POPUL_F), getAvgAlpha(POPUL_M), getAvgAlpha(POPUL_F+POPUL_M)]
            vAvgBeta = [getAvgBeta(POPUL_F), getAvgBeta(POPUL_M), getAvgBeta(POPUL_F+POPUL_M)]
            vCorrel_F, vCorrel_M = getCorrel(POPUL_F), getCorrel(POPUL_M)

            vDistAlpha_F, vDistAlpha_M = getDistAlpha(POPUL_F), getDistAlpha(POPUL_M)
            vDistBeta_F, vDistBeta_M = getDistBeta(POPUL_F), getDistBeta(POPUL_M)
            vDistAlpha, vDistBeta = getDistAlpha(POPUL_F+POPUL_M), getDistBeta(POPUL_M+POPUL_M)
            vDistStExpB_F, vDistStExpB_M = getDistStExpB(POPUL_F), getDistStExpB(POPUL_M)

            vAllData = vParam + vMaxFit + vAvgFit + vAvgFit_GtA + \
                       vAvgStExpA + vVarStExpA + vAvgStExpB + vVarStExpB + vAvgAlleleExpA_het + \
                       vFreqA + vGtA + vAvgAlpha + vAvgBeta + vCorrel_F + vCorrel_M
            dfGeneral.loc[len(dfGeneral)] = vAllData

            vDistData = vParam + vDistAlpha_F + vDistAlpha_M + vDistBeta_F + vDistBeta_M + vDistAlpha + vDistBeta + vDistStExpB_F + vDistStExpB_M
            dfDist.loc[len(dfDist)] = vDistData

        '''
        #####-----Check freq of A, if fixed then stops after MaxExtraGen generations
        freqA = getFreqA(POPUL_F+POPUL_M)
        if freqA == 0 or freqA == 1:
            extraGen += 1
        if extraGen == MaxExtraGen:
            break
        '''
        if gen == NUMGEN-1:
            COLNAME = ['ID', 'SEX', 'GTA', 'A1', 'A2', 'ALPHA1', 'ALPHA2', 'BETA1', 'BETA2', 'concA1', 'concA2', 'stconB', 'absFitness']
            dfIndiv = pd.DataFrame(columns=COLNAME)
            #[self.sex, self.gtA, self.a1, self.a2, self.alpha1, self.alpha2, self.beta1, self.beta2, self.concA1, self.concA2, self.stConB, self.absFitness]
            for i, INDIV in enumerate(POPUL_F+POPUL_M):
                SMY = INDIV.getSmy()
                dfIndiv.loc[len(dfIndiv)] = [i]+SMY

    #####-----return data
    return dfGeneral, dfDist, dfIndiv



####################--------------------Function that runs multiple replicates of simulation
####################    each replicate output file has a different name (increment by an integer)
def runReps(DIR, NAME, REPS, CONC_D, GA, KA, KB, ASAT, SR, S_F, S_M, GAMMA_F, GAMMA_M, UA, UB, N, PROPF, NUMGEN):
    ##########----------Run simulations for REPS number of replicates
    for r in range(REPS):
        print ('Running [{}] rep [{}]'.format(NAME, r + 1))
        ##########----------run simulation
        dfGeneral, dfDist, dfIndiv = runSim01(CONC_D, GA, KA, KB, ASAT, SR, S_F, S_M, GAMMA_F, GAMMA_M, UA, UB, N, PROPF, NUMGEN)
        ##########----------output data
        outGeneral = DIR + NAME + '.' + str(r + 1) + '.general.csv'
        outDist = DIR + NAME + '.' + str(r + 1) + '.dist.csv'
        outIndiv = DIR + NAME + '.' + str(r + 1) + '.indiv.csv'
        dfGeneral.to_csv(path_or_buf=outGeneral, index=False)
        dfDist.to_csv(path_or_buf=outDist, index=False)
        dfIndiv.to_csv(path_or_buf=outIndiv, index=False)


####################--------------------Function that reads parameters from csv file and runs simulations
def runParams(NPROC, PARFILE):
    ##########----------read parameter file
    try:
        PARAM = pd.read_csv(PARFILE)
    except OSError:
        print('ERROR: could not open [{}]'.format(PARFILE))
        sys.exit()
    ##########----------Check that all NAMES for outputing files are unique
    nRow = len(PARAM)
    nUniqName = len(pd.unique(PARAM['NAME']))
    if nUniqName != nRow:
        print('ERROR: one or more NAME are repeated: {}'.format(PARAM['NAME'].tolist()))
        sys.exit()
    ##########----------Run simulation for each set of parameters
    pool = Pool(processes=NPROC)
    for i in range(len(PARAM)):
        #####-----get params for row
        parRow = PARAM.loc[i]
        DIR, NAME, REPS, N, PROPF, NUMGEN = str(parRow['DIR']), str(parRow['NAME']), int(parRow['REPS']), int(
            parRow['N']), float(parRow['PROPF']), int(parRow['NUMGEN'])
        CONC_D, GA, KA, KB, ASAT = float(parRow['CONC_D']), float(parRow['GA']), float(parRow['KA']), float(
            parRow['KB']), float(parRow['ASAT'])
        SR, S_F, S_M, GAMMA_F, GAMMA_M, UA, UB = int(parRow['SR']), float(parRow['S_F']), float(parRow['S_M']), float(
            parRow['GAMMA_F']), float(parRow['GAMMA_M']), float(parRow['UA']), float(parRow['UB'])
        #print(DIR, NAME, REPS, CONC_D, GA, KA, KB, ASAT, S_F, S_M, GAMMA_F, GAMMA_M, UA, UB, N, PROPF, NUMGEN)
        #####-----run simlation if directory exists
        if os.path.isdir(DIR):
            print('Running simulations for [{}]'.format(NAME))
            pool.apply_async(runReps, args=(DIR, NAME, REPS, CONC_D, GA, KA, KB, ASAT, SR, S_F, S_M, GAMMA_F, GAMMA_M, UA, UB, N, PROPF, NUMGEN))
        else:
            print('DIR [{}] does not exist. Simulations for [{}] were not performed.'.format(DIR, NAME))
    #####-----makes sure all processes finish before continung script
    pool.close()
    pool.join()

####################--------------------Read arguments for MAIN
def parseArg():
    ##########----------Initialize parser
    parser = argparse.ArgumentParser()
    ##########----------Read arguments
    parser.add_argument("-n", "--NPROC", help="Number of processes to run simultaneously", type=int, required=True)
    parser.add_argument("-p", "--PARFILE", help="CSV file containing parameters for simulations", type=str, required=True)
    args = parser.parse_args()
    ##########----------Return arguments
    return args.NPROC, args.PARFILE

####################--------------------MAIN
def main():
    ##########----------get arguments from command line
    NPROC, PARFILE = parseArg()
    print("Num of processes to run: [{}]".format(NPROC))
    print("Reading parameters from: [{}]\n".format(PARFILE))
    ##########----------start program
    starttime = time.time()
    #PARFILE = 'D:/Users/eddie/DomSimData/TestMulti/PARAM.csv'
    runParams(NPROC, PARFILE)
    print('That took {} seconds'.format(time.time() - starttime))

if __name__ == "__main__":
    main()