+1 vote

I was implementing gillespie algorithm and when i ran the code it gave me a value error saying can't broadcast input array from shape(2) into shape(5). Can anyone please help me out?

import multiprocessing
import numpy as np
import scipy.stats as st
import numba
import matplotlib.pyplot as plt
import seaborn as sns

simple_update = np.array([[-1,1,0,0,0],
[1,-1,0,0,0],
[0,-1,1,0,0],
[0,1,-1,0,0],
[0,0,-1,1,0],
[0,0,1,-1,0],
[0,0,0,-1,1],
[0,0,0,1,-1]], dtype = np.int)

def simple_propensity(params, population):

alpha_1,alpha_2,alpha_3,alpha_4,beta_1,beta_2,beta_3,beta_4 = params

C0,C1,C2,C3,O = population

return np.array([alpha_1*C0,
beta_1*C1,
alpha_2*C1,
beta_2*C2,
alpha_3*C2,
beta_3*C3,
alpha_4*C3,
beta_4*O])

def sample_discrete(probs):
q = np.random.randn()

# Find index
i = 0
p_sum = 0.0
while p_sum < q:
p_sum += probs[i]
i += 1
return i - 1

def gillespie_draw(params, propensity_func, population):

props = propensity_func(params, population)

# Sum of propensities
props_sum = props.sum()

# Compute time
time = np.random.exponential(1.0 / props_sum)

# Compute discrete probabilities of each reaction
rxn_probs = props / props_sum

# Draw reaction from this distribution
rxn = sample_discrete(rxn_probs)

return rxn, time
def gillespie_ssa(params, propensity_func, update, population_0,
time_points):

pop_out = np.empty((len(time_points), update.shape[1]), dtype=np.int)

# Initialize and perform simulation
i_time = 1
i = 0
t = time_points[0]
population = population_0.copy()
pop_out[0,:] = population
while i < len(time_points):
while t < time_points[i_time]:
# draw the event and time step
event, dt = gillespie_draw(params, propensity_func, population)

# Update the population
population_previous = population.copy()
population += update[event,:]

# Increment time
t += dt

# Update the index
i = np.searchsorted(time_points > t, True)

# Update the population
pop_out[i_time:min(i,len(time_points))] = population_previous

# Increment index
i_time = i

return pop_out

# Specify parameters for calculation
params = np.array([10, 10, 0.4])
time_points = np.linspace(0, 50, 101)
population_0 = np.array([0, 0])
n_simulations = 100

# Seed random number generator for reproducibility
np.random.seed(42)

# Initialize output array
pops = np.empty((n_simulations, len(time_points), 2))

# Run the calculations
for i in range(n_simulations):
pops[i,:,:] = gillespie_ssa(params, simple_propensity, simple_update,population_0, time_points)

fig, ax = plt.subplots(1, 2, figsize=(14, 5))

for i in range(n_simulations):
ax[0].plot(time_points, pops[i,:,0], '-', lw=0.3, alpha=0.2,
color=sns.color_palette()[0])
# Set up subplots
fig, ax = plt.subplots(1, 2, figsize=(14, 5))

for i in range(n_simulations):
ax[0].plot(time_points, pops[i,:,0], '-', lw=0.3, alpha=0.2,
color=sns.color_palette()[0])
# Plot protein trajectories
for i in range(n_simulations):
ax[1].plot(time_points, pops[i,:,1], 'k-', lw=0.3, alpha=0.2,
color=sns.color_palette()[0])

ax[0].set_xlabel('time')
ax[1].set_xlabel(' time')
ax[0].set_ylabel('number of subunits')
ax[1].set_ylabel('number of calcium')
plt.tight_layout()
closed with the note: other
Dec 23, 2019 in Python
reshown Jan 13, 2020 2,665 views
Hello

Which line is giving you this error?
in the gillespie_ssa function

pop_out[0,:] = population
Hi @Satyam, did you solve the issue?

## Python Error ""ValueError: could not broadcast input array from shape (4) into shape (1000)""

There are better ways of achieving the ...READ MORE

## ValueError: could not broadcast input array from shape (4,1) into shape (4)

Hey @Giorgio, You can try this hope this ...READ MORE

## ValueError: could not broadcast input array from shape (360,270,3) into shape (360,280,3)

Hi@akhtar, In the above error it shows could not ...READ MORE

## ValueError: could not broadcast input array from shape (224,224,9) into shape (224,224)

img_shape = 224 test_data = [] test_labels = [] for ...READ MORE

+1 vote

## ValueError: operands could not be broadcast together with shapes (3,) (1000,)

This is the part of my code, why ...READ MORE

## how to print array integer without [] bracket in python like result = 1,2,3,4,5

Hey @abhijmr.143, you can print array integers ...READ MORE

## How do I obtain the index list in a NumPy Array of all the NaN values present using Python?

Hi, it is pretty simple, to be ...READ MORE