import numpy as np
import matplotlib.pyplot as plt
import random
from math import *
import operator as op
from sympy import FallingFactorial




def ncr(n,k):
    return reduce(lambda acc,m:acc*(n-k+m)/m, range(1,k+1),1) 

#this is a helper function to calclate the combinitorial function "N choose k"
#def ncr(n,r):
#    if n < r:
#        return 0
#    elif n == r:
#        return 1
#    else:
#        return int(FallingFactorial(n,r)/gamma(r+1))



def strings_to_floats(rate_list):
    l = [float(item) for item in rate_list]
    return l

# This function takes our data from the parsed .psc file and returns the
# reactions found in the data as a list
def ReactionList(data):
    R = []
    for item in data:
       if type(item) is list:
           R.append(item)
    return R

# This is a placeholder for our initial population vector
def PopulationVector(data):
    d = dict()
    for item in data:
        if type(item) is dict:
            d.update(item.items())

    return d
    
# This function retrieves the reaction rates from our dictionary and returns
# them in a list that will become our reaction rate vector. It takes our
# reaction rate dictionary as a parameter.
def Cmu(RList):
    c_mu = []
    for reactants,products,rate in RList:
        c_mu.append(rate)
    return c_mu


def calcM(RList):
    """This function returns the number of reactions in our dictionary"""
    M = len(RList)
    return M


def calcN(Rlist):
    species = sum([reactants + products 
                   for (reactants,products,rate) in Rlist],[])
    return len(set(species))

# this function is used to caclulate the discrete RV mu in step 2 of this
# algorithim. mu is the index of the reaction that occurs which leads to the
# next system state. The function takes as parameters: the dictionary a_0, a
# uniform RV r_2, and the sum of all values in a_0 (which is calculated in
# Step 1)
def calcMu(a_0,r_2,a_sum):
    s = 0
    v = int()
    while s < a_sum*r_2:
        s = s + a_0[v]
        v = v + 1
    return v

# This function calculates h_i for every species of reactants and stores them
# to a list which the function returns. The function takes our reaction list
# and our population dictionary as parameters 
def calc_h(RList,X_i):
    h_mu=[]
    for reactants,products,rate in RList:
        h_i = [ncr(int(X_i[species]),reactants.count(species))
                                    for species in set(reactants)]
        h_mu.append(reduce(op.mul,h_i,1))
    return h_mu


def stochsimm(PSCmodel):

    # Step 0: initilize values
    # call our reaction list
    RList = ReactionList(PSCmodel)
    RList = RList[0]
    # c_mu: the reaction constans for each defined reaction
    #       value is the average probability reaction mu will occur in the next
    #       infitesimal time period dt
    c_mu = Cmu(RList)
    c_mu = strings_to_floats(c_mu)
    # X_i: population of each reactant species
    X_i = PopulationVector(PSCmodel)
    # M: number of possible reactions in system
    M = calcM(RList)
    # N: number of molecule species in system
    N = calcN(RList) 
    # pop: the development of molecule X population over time
    pop = list()
    # time: the list of times when a defined reaction occurs
    time = list()
    # initial start time and reaction counter
    t = 0
    counter = 0
    #h_mu = calc_h(RList,X_i)

    # define conditions for when algorithim is run
    while counter <= 10000:
        # Step 1: parameter calculation
        # update value of h
        # h: the number of distinct reactant combinations availible in the state
        h_mu = calc_h(RList,X_i)
        # calculate each a_mu and take thier sum
        # a_mu:average probability reaction mu will occur given the curent state
        #      of the system
        a_0 = [c_mu[i]*h_mu[i] for i in range(M)]
        a_sum = sum(a_0)
        # When a_sum is 0 there is no posible reaction that can occour in the
        # System. When this occours, end the algorithim here
        if a_sum == 0:
            break
        # Step 2: Monte-Carlo simulation
        # generate r_1 and r_2
        r_1 = random.uniform(0,1)
        r_2 = random.uniform(0,1)
        # calculate the continuous RV tau
        # tau: time till reaction mu that leads to the next system state
        tau = (1/a_sum)*log(1/r_1)
        # calculate the discrete RV mu
        # mu: the index number that identifies the reaction that occurs
        mu = calcMu(a_0,r_2,a_sum)
        # Step 3: System update and data record
        # update the system time and add to the list
        t = t + tau
        time.append(t)
        # update the population levels in the system and add them to the list
        reaction = RList[mu-1]
        for reactants in reaction[0]:
            X_i[reactants] -= 1
        
        for products in reaction[1]:
            X_i[products] += 1
    
        
        pop.append(dict(X_i))
    
        #reaction counter increases
        counter += 1


    return pop
    #test code here
    
   
def getATP(population):
    ATP = [X_i['ATP'] for X_i in population]
    return ATP


    


    
print('loaded gill_alg')


7ad7d43e68bed0655670d2568d6a384a92a9b182
<mask token>


def ncr(n, k):
    return reduce(lambda acc, m: acc * (n - k + m) / m, range(1, k + 1), 1)


def strings_to_floats(rate_list):
    l = [float(item) for item in rate_list]
    return l


def ReactionList(data):
    R = []
    for item in data:
        if type(item) is list:
            R.append(item)
    return R


<mask token>


def calcM(RList):
    """This function returns the number of reactions in our dictionary"""
    M = len(RList)
    return M


def calcN(Rlist):
    species = sum([(reactants + products) for reactants, products, rate in
        Rlist], [])
    return len(set(species))


def calcMu(a_0, r_2, a_sum):
    s = 0
    v = int()
    while s < a_sum * r_2:
        s = s + a_0[v]
        v = v + 1
    return v


<mask token>

<mask token>


def ncr(n, k):
    return reduce(lambda acc, m: acc * (n - k + m) / m, range(1, k + 1), 1)


def strings_to_floats(rate_list):
    l = [float(item) for item in rate_list]
    return l


def ReactionList(data):
    R = []
    for item in data:
        if type(item) is list:
            R.append(item)
    return R


def PopulationVector(data):
    d = dict()
    for item in data:
        if type(item) is dict:
            d.update(item.items())
    return d


def Cmu(RList):
    c_mu = []
    for reactants, products, rate in RList:
        c_mu.append(rate)
    return c_mu


def calcM(RList):
    """This function returns the number of reactions in our dictionary"""
    M = len(RList)
    return M


def calcN(Rlist):
    species = sum([(reactants + products) for reactants, products, rate in
        Rlist], [])
    return len(set(species))


def calcMu(a_0, r_2, a_sum):
    s = 0
    v = int()
    while s < a_sum * r_2:
        s = s + a_0[v]
        v = v + 1
    return v


def calc_h(RList, X_i):
    h_mu = []
    for reactants, products, rate in RList:
        h_i = [ncr(int(X_i[species]), reactants.count(species)) for species in
            set(reactants)]
        h_mu.append(reduce(op.mul, h_i, 1))
    return h_mu


def stochsimm(PSCmodel):
    RList = ReactionList(PSCmodel)
    RList = RList[0]
    c_mu = Cmu(RList)
    c_mu = strings_to_floats(c_mu)
    X_i = PopulationVector(PSCmodel)
    M = calcM(RList)
    N = calcN(RList)
    pop = list()
    time = list()
    t = 0
    counter = 0
    while counter <= 10000:
        h_mu = calc_h(RList, X_i)
        a_0 = [(c_mu[i] * h_mu[i]) for i in range(M)]
        a_sum = sum(a_0)
        if a_sum == 0:
            break
        r_1 = random.uniform(0, 1)
        r_2 = random.uniform(0, 1)
        tau = 1 / a_sum * log(1 / r_1)
        mu = calcMu(a_0, r_2, a_sum)
        t = t + tau
        time.append(t)
        reaction = RList[mu - 1]
        for reactants in reaction[0]:
            X_i[reactants] -= 1
        for products in reaction[1]:
            X_i[products] += 1
        pop.append(dict(X_i))
        counter += 1
    return pop


<mask token>

<mask token>


def ncr(n, k):
    return reduce(lambda acc, m: acc * (n - k + m) / m, range(1, k + 1), 1)


def strings_to_floats(rate_list):
    l = [float(item) for item in rate_list]
    return l


def ReactionList(data):
    R = []
    for item in data:
        if type(item) is list:
            R.append(item)
    return R


def PopulationVector(data):
    d = dict()
    for item in data:
        if type(item) is dict:
            d.update(item.items())
    return d


def Cmu(RList):
    c_mu = []
    for reactants, products, rate in RList:
        c_mu.append(rate)
    return c_mu


def calcM(RList):
    """This function returns the number of reactions in our dictionary"""
    M = len(RList)
    return M


def calcN(Rlist):
    species = sum([(reactants + products) for reactants, products, rate in
        Rlist], [])
    return len(set(species))


def calcMu(a_0, r_2, a_sum):
    s = 0
    v = int()
    while s < a_sum * r_2:
        s = s + a_0[v]
        v = v + 1
    return v


def calc_h(RList, X_i):
    h_mu = []
    for reactants, products, rate in RList:
        h_i = [ncr(int(X_i[species]), reactants.count(species)) for species in
            set(reactants)]
        h_mu.append(reduce(op.mul, h_i, 1))
    return h_mu


def stochsimm(PSCmodel):
    RList = ReactionList(PSCmodel)
    RList = RList[0]
    c_mu = Cmu(RList)
    c_mu = strings_to_floats(c_mu)
    X_i = PopulationVector(PSCmodel)
    M = calcM(RList)
    N = calcN(RList)
    pop = list()
    time = list()
    t = 0
    counter = 0
    while counter <= 10000:
        h_mu = calc_h(RList, X_i)
        a_0 = [(c_mu[i] * h_mu[i]) for i in range(M)]
        a_sum = sum(a_0)
        if a_sum == 0:
            break
        r_1 = random.uniform(0, 1)
        r_2 = random.uniform(0, 1)
        tau = 1 / a_sum * log(1 / r_1)
        mu = calcMu(a_0, r_2, a_sum)
        t = t + tau
        time.append(t)
        reaction = RList[mu - 1]
        for reactants in reaction[0]:
            X_i[reactants] -= 1
        for products in reaction[1]:
            X_i[products] += 1
        pop.append(dict(X_i))
        counter += 1
    return pop


def getATP(population):
    ATP = [X_i['ATP'] for X_i in population]
    return ATP


<mask token>

import numpy as np
import matplotlib.pyplot as plt
import random
from math import *
import operator as op
from sympy import FallingFactorial


def ncr(n, k):
    return reduce(lambda acc, m: acc * (n - k + m) / m, range(1, k + 1), 1)


def strings_to_floats(rate_list):
    l = [float(item) for item in rate_list]
    return l


def ReactionList(data):
    R = []
    for item in data:
        if type(item) is list:
            R.append(item)
    return R


def PopulationVector(data):
    d = dict()
    for item in data:
        if type(item) is dict:
            d.update(item.items())
    return d


def Cmu(RList):
    c_mu = []
    for reactants, products, rate in RList:
        c_mu.append(rate)
    return c_mu


def calcM(RList):
    """This function returns the number of reactions in our dictionary"""
    M = len(RList)
    return M


def calcN(Rlist):
    species = sum([(reactants + products) for reactants, products, rate in
        Rlist], [])
    return len(set(species))


def calcMu(a_0, r_2, a_sum):
    s = 0
    v = int()
    while s < a_sum * r_2:
        s = s + a_0[v]
        v = v + 1
    return v


def calc_h(RList, X_i):
    h_mu = []
    for reactants, products, rate in RList:
        h_i = [ncr(int(X_i[species]), reactants.count(species)) for species in
            set(reactants)]
        h_mu.append(reduce(op.mul, h_i, 1))
    return h_mu


def stochsimm(PSCmodel):
    RList = ReactionList(PSCmodel)
    RList = RList[0]
    c_mu = Cmu(RList)
    c_mu = strings_to_floats(c_mu)
    X_i = PopulationVector(PSCmodel)
    M = calcM(RList)
    N = calcN(RList)
    pop = list()
    time = list()
    t = 0
    counter = 0
    while counter <= 10000:
        h_mu = calc_h(RList, X_i)
        a_0 = [(c_mu[i] * h_mu[i]) for i in range(M)]
        a_sum = sum(a_0)
        if a_sum == 0:
            break
        r_1 = random.uniform(0, 1)
        r_2 = random.uniform(0, 1)
        tau = 1 / a_sum * log(1 / r_1)
        mu = calcMu(a_0, r_2, a_sum)
        t = t + tau
        time.append(t)
        reaction = RList[mu - 1]
        for reactants in reaction[0]:
            X_i[reactants] -= 1
        for products in reaction[1]:
            X_i[products] += 1
        pop.append(dict(X_i))
        counter += 1
    return pop


def getATP(population):
    ATP = [X_i['ATP'] for X_i in population]
    return ATP


print('loaded gill_alg')

import numpy as np
import matplotlib.pyplot as plt
import random
from math import *
import operator as op
from sympy import FallingFactorial




def ncr(n,k):
    return reduce(lambda acc,m:acc*(n-k+m)/m, range(1,k+1),1) 

#this is a helper function to calclate the combinitorial function "N choose k"
#def ncr(n,r):
#    if n < r:
#        return 0
#    elif n == r:
#        return 1
#    else:
#        return int(FallingFactorial(n,r)/gamma(r+1))



def strings_to_floats(rate_list):
    l = [float(item) for item in rate_list]
    return l

# This function takes our data from the parsed .psc file and returns the
# reactions found in the data as a list
def ReactionList(data):
    R = []
    for item in data:
       if type(item) is list:
           R.append(item)
    return R

# This is a placeholder for our initial population vector
def PopulationVector(data):
    d = dict()
    for item in data:
        if type(item) is dict:
            d.update(item.items())

    return d
    
# This function retrieves the reaction rates from our dictionary and returns
# them in a list that will become our reaction rate vector. It takes our
# reaction rate dictionary as a parameter.
def Cmu(RList):
    c_mu = []
    for reactants,products,rate in RList:
        c_mu.append(rate)
    return c_mu


def calcM(RList):
    """This function returns the number of reactions in our dictionary"""
    M = len(RList)
    return M


def calcN(Rlist):
    species = sum([reactants + products 
                   for (reactants,products,rate) in Rlist],[])
    return len(set(species))

# this function is used to caclulate the discrete RV mu in step 2 of this
# algorithim. mu is the index of the reaction that occurs which leads to the
# next system state. The function takes as parameters: the dictionary a_0, a
# uniform RV r_2, and the sum of all values in a_0 (which is calculated in
# Step 1)
def calcMu(a_0,r_2,a_sum):
    s = 0
    v = int()
    while s < a_sum*r_2:
        s = s + a_0[v]
        v = v + 1
    return v

# This function calculates h_i for every species of reactants and stores them
# to a list which the function returns. The function takes our reaction list
# and our population dictionary as parameters 
def calc_h(RList,X_i):
    h_mu=[]
    for reactants,products,rate in RList:
        h_i = [ncr(int(X_i[species]),reactants.count(species))
                                    for species in set(reactants)]
        h_mu.append(reduce(op.mul,h_i,1))
    return h_mu


def stochsimm(PSCmodel):

    # Step 0: initilize values
    # call our reaction list
    RList = ReactionList(PSCmodel)
    RList = RList[0]
    # c_mu: the reaction constans for each defined reaction
    #       value is the average probability reaction mu will occur in the next
    #       infitesimal time period dt
    c_mu = Cmu(RList)
    c_mu = strings_to_floats(c_mu)
    # X_i: population of each reactant species
    X_i = PopulationVector(PSCmodel)
    # M: number of possible reactions in system
    M = calcM(RList)
    # N: number of molecule species in system
    N = calcN(RList) 
    # pop: the development of molecule X population over time
    pop = list()
    # time: the list of times when a defined reaction occurs
    time = list()
    # initial start time and reaction counter
    t = 0
    counter = 0
    #h_mu = calc_h(RList,X_i)

    # define conditions for when algorithim is run
    while counter <= 10000:
        # Step 1: parameter calculation
        # update value of h
        # h: the number of distinct reactant combinations availible in the state
        h_mu = calc_h(RList,X_i)
        # calculate each a_mu and take thier sum
        # a_mu:average probability reaction mu will occur given the curent state
        #      of the system
        a_0 = [c_mu[i]*h_mu[i] for i in range(M)]
        a_sum = sum(a_0)
        # When a_sum is 0 there is no posible reaction that can occour in the
        # System. When this occours, end the algorithim here
        if a_sum == 0:
            break
        # Step 2: Monte-Carlo simulation
        # generate r_1 and r_2
        r_1 = random.uniform(0,1)
        r_2 = random.uniform(0,1)
        # calculate the continuous RV tau
        # tau: time till reaction mu that leads to the next system state
        tau = (1/a_sum)*log(1/r_1)
        # calculate the discrete RV mu
        # mu: the index number that identifies the reaction that occurs
        mu = calcMu(a_0,r_2,a_sum)
        # Step 3: System update and data record
        # update the system time and add to the list
        t = t + tau
        time.append(t)
        # update the population levels in the system and add them to the list
        reaction = RList[mu-1]
        for reactants in reaction[0]:
            X_i[reactants] -= 1
        
        for products in reaction[1]:
            X_i[products] += 1
    
        
        pop.append(dict(X_i))
    
        #reaction counter increases
        counter += 1


    return pop
    #test code here
    
   
def getATP(population):
    ATP = [X_i['ATP'] for X_i in population]
    return ATP


    


    
print('loaded gill_alg')

