Singular Matrix Error

[BryanRumsey]

devils_v0_month.ipynb.zip

Data File: DFDT Population Devil Population

We are consistently getting the following error using the above model:

---------------------------------------------------------------------------
LinAlgError                               Traceback (most recent call last)
<ipython-input-15-a92b38b04305> in <module>
      7     # First compute the fixed(observed) mean
      8     abc.compute_fixed_mean(chunk_size=2)
----> 9     res = abc.infer(num_samples=100, batch_size=10, chunk_size=2)

/opt/conda/lib/python3.8/site-packages/sciope/inference/rep_smc_abc.py in infer(self, num_samples, alpha, R_trial, c, p_min, batch_size, chunk_size)
    208                 for i in range(n_cull, num_samples):
    209                     perturb_tasks.append(self._perturb_resample(population[i, :], distances[i], R_trial, tol))
--> 210                 res, = dask.compute(perturb_tasks)
    211 
    212                 # Update the population with the perturbed population

/opt/conda/lib/python3.8/site-packages/dask/base.py in compute(*args, **kwargs)
    566         postcomputes.append(x.__dask_postcompute__())
    567 
--> 568     results = schedule(dsk, keys, **kwargs)
    569     return repack([f(r, *a) for r, (f, a) in zip(results, postcomputes)])
    570 

/opt/conda/lib/python3.8/site-packages/distributed/client.py in get(self, dsk, keys, workers, allow_other_workers, resources, sync, asynchronous, direct, retries, priority, fifo_timeout, actors, **kwargs)
   2669                     should_rejoin = False
   2670             try:
-> 2671                 results = self.gather(packed, asynchronous=asynchronous, direct=direct)
   2672             finally:
   2673                 for f in futures.values():

/opt/conda/lib/python3.8/site-packages/distributed/client.py in gather(self, futures, errors, direct, asynchronous)
   1946             else:
   1947                 local_worker = None
-> 1948             return self.sync(
   1949                 self._gather,
   1950                 futures,

/opt/conda/lib/python3.8/site-packages/distributed/client.py in sync(self, func, asynchronous, callback_timeout, *args, **kwargs)
    843             return future
    844         else:
--> 845             return sync(
    846                 self.loop, func, *args, callback_timeout=callback_timeout, **kwargs
    847             )

/opt/conda/lib/python3.8/site-packages/distributed/utils.py in sync(loop, func, callback_timeout, *args, **kwargs)
    323     if error[0]:
    324         typ, exc, tb = error[0]
--> 325         raise exc.with_traceback(tb)
    326     else:
    327         return result[0]

/opt/conda/lib/python3.8/site-packages/distributed/utils.py in f()
    306             if callback_timeout is not None:
    307                 future = asyncio.wait_for(future, callback_timeout)
--> 308             result[0] = yield future
    309         except Exception:
    310             error[0] = sys.exc_info()

/opt/conda/lib/python3.8/site-packages/tornado/gen.py in run(self)
    733 
    734                     try:
--> 735                         value = future.result()
    736                     except Exception:
    737                         exc_info = sys.exc_info()

/opt/conda/lib/python3.8/site-packages/distributed/client.py in _gather(self, futures, errors, direct, local_worker)
   1811                             exc = CancelledError(key)
   1812                         else:
-> 1813                             raise exception.with_traceback(traceback)
   1814                         raise exc
   1815                     if errors == "skip":

/opt/conda/lib/python3.8/site-packages/sciope/inference/rep_smc_abc.py in _perturb_resample()
    113             # Metropolis Perturbation
    114             if nume > 0:
--> 115                 nume = nume * self.perturbation_kernel.pdf(param.reshape(1, -1), proposal.reshape(1, -1))[0, 0,]
    116                 dnm = self.prior_function.pdf(param) * \
    117                       self.perturbation_kernel.pdf(proposal.reshape(1, -1), param.reshape(1, -1))[0, 0]

/opt/conda/lib/python3.8/site-packages/sciope/utilities/perturbationkernels/multivariate_normal.py in pdf()
     43                 pdfs.append(multivariate_normal.logpdf(x, x0[i], self.cov))
     44             else:
---> 45                 pdfs.append(multivariate_normal.pdf(x, x0[i], self.cov))
     46         return np.vstack(pdfs)
     47 

/opt/conda/lib/python3.8/site-packages/scipy/stats/_multivariate.py in pdf()
    514         dim, mean, cov = self._process_parameters(None, mean, cov)
    515         x = self._process_quantiles(x, dim)
--> 516         psd = _PSD(cov, allow_singular=allow_singular)
    517         out = np.exp(self._logpdf(x, mean, psd.U, psd.log_pdet, psd.rank))
    518         return _squeeze_output(out)

/opt/conda/lib/python3.8/site-packages/scipy/stats/_multivariate.py in __init__()
    163         d = s[s > eps]
    164         if len(d) < len(s) and not allow_singular:
--> 165             raise np.linalg.LinAlgError('singular matrix')
    166         s_pinv = _pinv_1d(s, eps)
    167         U = np.multiply(u, np.sqrt(s_pinv))

LinAlgError: singular matrix

[prasi372]

Hi, thanks for reporting the crash. It is due to the covariance matrix being ill-conditioned, while applying the perturbation during SMC ABC. The ideal way to solve this would be to rescale the range of the prior. I performed a quick test by changing the prior and simulator as,

from sciope.utilities.priors import uniform_prior

# take default from mode 1 as reference
default_param = np.array(list(model.listOfParameters.items()))[:, 1]

parameter_names = []
bound = []
for exp in default_param:
    bound.append(float(exp.expression))
    parameter_names.append(exp.name)

bound = np.log(bound)

# Set the bounds
bound = np.array(bound)
dmin = bound * 0.1
dmax = bound * 2.0

# Here we use uniform prior
uni_prior = uniform_prior.UniformPrior(dmin, dmax)

Here we have just added a np.log transformation to the bounds. Then we do the opposite while simulating as,

# Here we use the GillesPy2 Solver
def simulator(params, model):
    params = np.exp(params)
    model_update = set_model_parameters(params, model)

    res = model_update.run(**kwargs)
    devils = res['Devils']
    infected = res['I']

    return np.vstack([devils, infected]).reshape(1, 2, -1)

Note that only the first line (np.exp) has been added. This should work, or you could also rescale to [0,1] for example. Do let us know if this consistently fixes the issue.

[BryanRumsey]

We attempted this solution and are consistently get this error:

/opt/conda/lib/python3.8/site-packages/sciope/inference/rep_smc_abc.py:224: RuntimeWarning:

divide by zero encountered in double_scalars

---------------------------------------------------------------------------
OverflowError                             Traceback (most recent call last)
<ipython-input-18-74032fd839d8> in <module>
      6     # First compute the fixed(observed) mean
      7     abc.compute_fixed_mean(chunk_size=2)
----> 8     res = abc.infer(num_samples=100, batch_size=10, chunk_size=2)

/opt/conda/lib/python3.8/site-packages/sciope/inference/rep_smc_abc.py in infer(self, num_samples, alpha, R_trial, c, p_min, batch_size, chunk_size)
    222                 N_acc = np.sum(N_accs)
    223 
--> 224                 R = int(round(np.log(c) / np.log(1 - p_acc)))
    225 
    226                 # Perturb again with better estimate

OverflowError: cannot convert float infinity to integer

[prasi372]

The latest commit on sciope master now takes care of this scenario. I have also updated your notebook to use vanilla SMC ABC, can you please check if the attached notebook works as expected for you, with the latest code on sciope master? Please note that the replenishment variant might not work and needs more investigation and setup. The vanialla SMC should work as expected. devils_test.zip