ReactionSystem.py 文件源码

def run(self, t,seed=None,plot=False):
            """
            Gillespie Implementation
            """
            l = self.reactions[:,0]
            r = self.reactions[:,1]
            change=r-l
            time =0
            times =[]
            output =[]
            times.append(time)
            output.append(self.population)
            while time<t:           
                lam = np.prod(comb(self.population,l)*factorial(l),axis=-1)*self.rates
                lam_sum = np.sum(lam)
                dt = np.log(1.0/np.random.uniform(0,1))*(1.0/(lam_sum))
                react = np.where(np.random.multinomial(1,lam/lam_sum)==1)[0][0]
                self.population = self.population + change[react]
                time += dt
                output.append(self.population)
                times.append(time)
            output = np.array(output)
            times = np.array(times)
            if plot:
                for i in xrange(output.shape[1]):
                    label = self.species[i]
                    plt.plot(times, output[:, i], label=label)
                plt.legend(loc='best')
                plt.xlabel('t')
                plt.title('Population Profile')
                plt.grid()
                plt.show()
            return times,output


# def main():
#   reactions = [[[1,0],[0,1]], [[0,1],[1,0]]]
#   rates = [1,1]
#   system = MassActionSystem(reactions,rates)
#   system.set_concentrations([2.,4.])
#   system.run()
#   print system.current_concentrations()
#   print system.dydt()
#   system = StochasticSystem(reactions,rates)
#   population = [10,0]
#   system.set_population(population)
#   t,o = system.run(10,plot=True)
#   return