LinAlgError raised when running SMC ABC

[BryanRumsey]

Error

Traceback (most recent call last):
  File "/stochss/stochss/handlers/util/scripts/start_job.py", line 120, in <module>
    job.run(verbose=args.verbose)
  File "/stochss/stochss/handlers/util/model_inference.py", line 958, in run
    results = smc_abc_inference.infer(*infer_args, **infer_kwargs)
  File "/opt/conda/lib/python3.8/site-packages/sciope/inference/smc_abc.py", line 220, in infer
    kweights = self.perturbation_kernel.pdf(population, new_samples)
  File "/opt/conda/lib/python3.8/site-packages/sciope/utilities/perturbationkernels/multivariate_normal.py", line 45, in pdf
    pdfs.append(multivariate_normal.pdf(x, x0[i], self.cov))
  File "/opt/conda/lib/python3.8/site-packages/scipy/stats/_multivariate.py", line 580, in pdf
    params = self._process_parameters(mean, cov, allow_singular)
  File "/opt/conda/lib/python3.8/site-packages/scipy/stats/_multivariate.py", line 417, in _process_parameters
    psd = _PSD(cov, allow_singular=allow_singular)
  File "/opt/conda/lib/python3.8/site-packages/scipy/stats/_multivariate.py", line 172, in __init__
    raise np.linalg.LinAlgError(msg)
numpy.linalg.LinAlgError: When `allow_singular is False`, the input matrix must be symmetric positive definite.

Model

def create_genetic_toggle_switch(parameter_values=None):
    model = gillespy2.Model(name='Genetic_Toggle_Switch')
    model.volume = 1

    # Variables
    U = gillespy2.Species(name='U', initial_value=10, mode='discrete')
    V = gillespy2.Species(name='V', initial_value=10, mode='discrete')
    model.add_species([
        U, V
    ])

    # Parameters
    alpha1 = gillespy2.Parameter(name='alpha1', expression='1')
    alpha2 = gillespy2.Parameter(name='alpha2', expression='1')
    beta = gillespy2.Parameter(name='beta', expression='2')
    gamma = gillespy2.Parameter(name='gamma', expression='2')
    mu = gillespy2.Parameter(name='mu', expression='1')
    model.add_parameter([
        alpha1, alpha2, beta, gamma, mu
    ])

    # Reactions
    cu = gillespy2.Reaction(
        name='cu',
        reactants={}, products={'U': 1},
        propensity_function='alpha1/(1+pow(V,beta))',
        ode_propensity_function='alpha1/(1+pow(V,beta))'
    )
    cv = gillespy2.Reaction(
        name='cv',
        reactants={}, products={'V': 1},
        propensity_function='alpha2/(1+pow(U,gamma))',
        ode_propensity_function='alpha2/(1+pow(U,gamma))'
    )
    du = gillespy2.Reaction(
        name='du', rate='mu',
        reactants={'U': 1}, products={}
    )
    dv = gillespy2.Reaction(
        name='dv', rate='mu',
        reactants={'V': 1}, products={}
    )
    model.add_reaction([
        cu, cv, du, dv
    ])

    # Timespan
    tspan = gillespy2.TimeSpan.arange(1, t=100)
    model.timespan(tspan)
    return model

Inference

def configure_simulation(trajectories=1):
    solver = gillespy2.SSACSolver(model=model, delete_directory=False)
    kwargs = {
        'number_of_trajectories': trajectories,
        'seed': None,
        'solver': solver,
    }
    return kwargs

kwargs = configure_simulation(100)
results = model.run(**kwargs)

unshaped_obs_data = results.to_array()
shaped_obs_data = unshaped_obs_data.swapaxes(1, 2)
obs_data = shaped_obs_data[:,1:, :]
print(obs_data.shape)

def process(raw_results):
    return raw_results.to_array().swapaxes(1, 2)[:,1:, :]

def simulator(parameter_point):
    model = create_genetic_toggle_switch()

    labels = [
        'alpha1', 'alpha2', 'beta', 'gamma', 'mu'
    ]
    for ndx, parameter in enumerate(parameter_point):
        model.listOfParameters[labels[ndx]].expression = str(parameter)

    kwargs = configure_simulation()
    raw_results = model.run(**kwargs)

    return process(raw_results)

values = numpy.array([
    1.0, 1.0, 2.0, 2.0, 1.0
])
parameter_names = [
    'alpha1', 'alpha2', 'beta', 'gamma', 'mu'
]

lower_bounds = numpy.array([
    0.5, 0.5, 1, 1, 0.5
])
upper_bounds = numpy.array([
    1.5, 1.5, 3, 3, 1.5
])

uni_prior = uniform_prior.UniformPrior(lower_bounds, upper_bounds)

summ_func = identity.Identity()

eps_selector = RelativeEpsilonSelector(20, max_rounds=1)

c = Client()

smc_abc_inference = smc_abc.SMCABC(
    obs_data, sim=simulator, prior_function=uni_prior, summaries_function=summ_func.compute
)

smc_abc_results = smc_abc_inference.infer(
    num_samples=10, batch_size=10, eps_selector=eps_selector, chunk_size=10
)