-1

I'm trying to simulate Gillespies algorithm for cellular systems in python 3.8 using a simple system with the following entities Enzyme, Substrate, Enzyme-Substrate complex, Product.

I have the following code that calculates propensity functions for a series of reactions which are represented as rows of an array: 

propensity = np.zeros(len(LHS))

def propensity_calc(LHS, popul_num, stoch_rate):
    for row in range(len(LHS)):     
            a = stoch_rate[row]        
            for i in range(len(popul_num)):      
                if (popul_num[i] >= LHS[row, i]):             
                    binom_rxn = (binom(popul_num[i], LHS[row, i]))    
                    a = a*binom_rxn            
                else: 
                    a = 0   
                    break 
            propensity[row] = a
    return propensity.astype(float)

The input arrays are as follows: 

popul_num = np.array([200, 100, 0, 0])
LHS = np.array([[1,1,0,0], [0,0,1,0], [0,0,1,0]])
stoch_rate = np.array([0.0016, 0.0001, 0.1000])

The function works as expected, until I try to call it from within the following while loop: 

while tao < tmax: 
propensity_calc(LHS, popul_num, stoch_rate)
a0 = sum(propensity)
if a0 <= 0:
    break
else:
    t = np.random.exponential(a0)   
    print(t)   
    # sample time system stays in given state from exponential distribution of the propensity sum. 
    if tao + t > tmax: 
        tao = tmax
        break
j = stats.rv_discrete(name="Reaction index", values=(num_rxn, rxn_probability)).rvs()   # array of reactions increasing by 1 until they get to the same length/size as rxn_probability
print(j)
tao = tao + t       
popul_num = popul_num + state_change_matrix[j]      # update state of system/ popul_num

The other variables in the while loop are as follows:

a0 = sum(propensity)

def prob_rxn_fires(propensity, a0):
prob = propensity/a0   
return prob 
rxn_probability = (prob_rxn_fires(propensity, a0))

num_rxn = np.arange(1, rxn_probability.size + 1).reshape(rxn_probability.shape)

When I run the code calling the calc_propensity function from within the while loop, it makes it through the first iteration of the while loop with the following error: 

ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()

The error is first thrown on the following line of the calc_propensity function:

if (popul_num[i] >= LHS[row, i]):

but for some reason the code keeps running until it reaches the same line in the calc_propensity function but in the second function call (the while loop) and I dont understand why?

cheers

hpaulj
  • 221,503
  • 14
  • 230
  • 353
Mike5298
  • 377
  • 1
  • 13
  • You shouldn't use `len` with numpy arrays, try `myarray.shape[dimension]` instead – Jakob Filser May 05 '20 at 14:35
  • 1
    Please [edit] your question and format your code properly. Python is indentation sensitive, several parts have incorrect syntax at the moment. – MisterMiyagi May 05 '20 at 14:35
  • Does this answer your question? [Use a.any() or a.all()](https://stackoverflow.com/questions/34472814/use-a-any-or-a-all) – MisterMiyagi May 05 '20 at 14:36
  • 1
    It looks to my like the `popul_num` and/or `LHS` object is not what you think it is. Printing it right before the line where it crashes will give you insight into what you are actually trying to compare. – Jakob Filser May 05 '20 at 14:45
  • It runs ok for me (commenting out the stuff that uses undefined `binom`), You have to do the debugging, starting with identifying the problem variables just prior to the error. Error's like this arise when numpy arrays are used in a comparison (like `>=`). The result is a boolean array, which can't be used in an `if` statement. Something must be changing the dimensions of the arrays being indexed in that statement. – hpaulj May 05 '20 at 15:57

1 Answers1

2

It appears that the values in the if statement are not to be used in the if conditional in the first place. You should use the .any() or .all() functions. Here is a link to illustrate what the ValueError stands for:

https://sopython.com/canon/119/the-truth-value-of-an-array-with-more-than-one-element-is-ambiguous/

Nachiket Kulkarni
  • 121
  • 1
  • 7
  • From what it looks like, `popul_num` is *supposed* to be a rank 1 array and `LHS` is *supposed* to be a rank 2 array. In that case the comparison would be perfectly fine, so I guess the problem is more that the objects aren't actually what they should be – Jakob Filser May 05 '20 at 14:56
  • Thanks a lot for the suggestions! I'm just not sure why the calc_propensity function works fine when its not called from inside the while loop though but then crashes when it is? – Mike5298 May 06 '20 at 07:42