ValueError could not broadcast input array from shape 2 into shape 5 [closed]

Question

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()

Comments

Hello

Which line is giving you this error?

---

in the gillespie_ssa function

pop_out[0,:] = population

---

Hi @Satyam, did you solve the issue?