CODEWITHSUNDEEP
pythonnumpylambda
Patrickens
Patrickens

Reputation: 323

Lambda function does not return correct value

Im trying to make a variant of the Gillespie algorithm, and to determine the reaction propensities Im trying to automatically generate the propensity vector using lambda expressions. However when creating SSA.P all goes wrong. The last loop in the block of code, PROPLOOP, returns two propensities, where the one generated using P_alternative is the correct one. The question is: how do I get the same values for SSA.P as for SSA.P_alternative?

import numpy as np
from numpy.random import uniform
class Markov:
  def __init__(self,z0,t0,tf,rates,stoich):
    self.S=stoich

    self.z0=z0
    self.rates=rates

    self.P=self.propensities()
    self.P_alternative=[
      lambda z,rate:(0.5*rate[0]*z[0]*(z[0]-1)),
      lambda z,rate:rate[1]*np.prod(z[0]),
      lambda z,rate:rate[2]*np.prod(z[1]),
      lambda z,rate:rate[3]*np.prod(z[1]),
      lambda z,rate:rate[4]*np.prod(z[np.array([0,1])]),
      lambda z,rate:rate[5]]

    self.t0=t0
    self.tf=tf


  def propensities(self):
    prop=[]
    for i,reac in enumerate(self.S.T):
      if all(z>=0 for z in reac):
        prop.append(lambda z,rate:rate[i])

      if any(z==-1 for z in reac):
        j=np.where(reac==-1)[0]
        prop.append(lambda z,rate:rate[i]*np.prod(z[j]))

      if any(z==-2 for z in reac):
        j=np.where(reac==-2)[0][0]
        prop.append(lambda z,rate:(0.5*rate[i]*z[j]*(z[j]-1))[0])

    return prop


stoich=np.array([
        [-2, -1,  2,  0, -1,  0],
        [ 1,  0, -1, -1, -1,  1],
        [ 0,  0,  0,  1,  1,  0]])

rates=np.array([1.0,0.02,200.0,0.0004,0.9,0.9])

z0=np.array([540,730,0])

SSA=Markov(z0=z0,t0=0,tf=100,rates=rates,stoich=stoich)

#PROPLOOP; the values should be equal for both SSA.P and SSA.P_alternative, where SSA.P_alternative is the correct one
for i in xrange(len(SSA.P)):
  print "Inexplicably wrong",SSA.P[i](z0,rates)
  print "Correct answer",SSA.P_alternative[i](z0,rates), "\n"

output is:

Inexplicably wrong 130977.0
Correct answer 145530.0 

Inexplicably wrong 354780.0
Correct answer 10.8 

Inexplicably wrong 354780.0
Correct answer 146000.0 

Inexplicably wrong 354780.0
Correct answer 0.292 

Correct answer 354780.0
Correct answer 354780.0 

Inexplicably wrong 0.9
Correct answer 0.9

Upvotes: 4

Views: 896

Answers (1)

Blckknght
Blckknght

Reputation: 104722

The issue is that you're creating your lambda functions in a loop, and they refer to the variables i and j that may change as the loop goes on.

The lambda doesn't copy the values of i or j when it is created, it just keeps a reference to the namespace that they are defined in. When it uses the variables when it is called later, it looks them up in that namespace. Since your lambdas get called after the loop (and indeed, the whole function) has ended, they all see the final values the variables were given, which is not what you intended. This explains why the two versions of your code give the same output on the last iteration. The final value of i and j is the expected one for the last lambda function.

You can work around this issue by making the lambda keep a copy of the current value of i and j when it is defined. The easiest way to do this is with a default argument:

for i,reac in enumerate(self.S.T):
  if all(z>=0 for z in reac):
    prop.append(lambda z, rate, i=i: rate[i]) # add i=i here and further down

  if any(z==-1 for z in reac):
    j=np.where(reac==-1)[0]
    prop.append(lambda z, rate, i=i, j=j: rate[i]*np.prod(z[j]))

  if any(z==-2 for z in reac):
    j=np.where(reac==-2)[0][0]
    prop.append(lambda z, rate, i=i, j=j: (0.5*rate[i]*z[j]*(z[j]-1))[0])

The i=i (and j=j where necessary) in the lambda definitions makes the variables arguments of the lambda function with a default value that is the current value of i (and j) in the outer namespace. Since you only pass two arguments when you call the lambda function, the saved default values will be used.

Upvotes: 6

Related Questions