Backends other than NDArray don't work with sequential sampling

[delwaterman]

Toy example of Backends does not work with default chains

While trying to use a toy example to test out Backends for PyMC3 3.3, I found that none of the backends (except the default in memory one) work with multiple chains:

Normal works:

17:20 $ python
Python 3.6.1 (default, Nov  8 2017, 14:29:33)
[GCC 4.2.1 Compatible Apple LLVM 7.0.2 (clang-700.1.81)] on darwin
Type "help", "copyright", "credits" or "license" for more information.
>>> import pymc3 as pm
>>> model = pm.Model()
>>> with model:
...     a = pm.Normal('a', mu=0, sd=1)
...     trace = pm.sample(1000, n_init=1000, cores=1, njobs=1)
...
Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
Sequential sampling (2 chains in 1 job)
NUTS: [a]
100%|████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 1500/1500 [00:01<00:00, 1445.02it/s]
100%|████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 1500/1500 [00:00<00:00, 2644.67it/s]

Text Backend Fails:

17:11 $ python
Python 3.6.1 (default, Nov  8 2017, 14:29:33)
[GCC 4.2.1 Compatible Apple LLVM 7.0.2 (clang-700.1.81)] on darwin
Type "help", "copyright", "credits" or "license" for more information.
>>> import pymc3 as pm
>>> model = pm.Model()
>>> with model:
...     a = pm.Normal('a', mu=0, sd=1)
...     db_text = pm.backends.Text("text-test-42")
...     trace = pm.sample(1000, n_init=1000, trace=db_text, cores=1, njobs=1)
...
Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
Sequential sampling (2 chains in 1 job)
NUTS: [a]
100%|████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 1500/1500 [00:00<00:00, 2404.39it/s]
100%|████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 1500/1500 [00:00<00:00, 2485.21it/s]
Traceback (most recent call last):
  File "<stdin>", line 3, in <module>
  File "/Users/orion.delwaterman/.pyenv/versions/3.6.1/envs/hiring-horizons/lib/python3.6/site-packages/pymc3/sampling.py", line 439, in sample
    trace = _sample_many(**sample_args)
  File "/Users/orion.delwaterman/.pyenv/versions/3.6.1/envs/hiring-horizons/lib/python3.6/site-packages/pymc3/sampling.py", line 494, in _sample_many
    return MultiTrace(traces)
  File "/Users/orion.delwaterman/.pyenv/versions/3.6.1/envs/hiring-horizons/lib/python3.6/site-packages/pymc3/backends/base.py", line 265, in __init__
    raise ValueError("Chains are not unique.")
ValueError: Chains are not unique.

SQLite Fails:

17:16 $ python
Python 3.6.1 (default, Nov  8 2017, 14:29:33)
[GCC 4.2.1 Compatible Apple LLVM 7.0.2 (clang-700.1.81)] on darwin
Type "help", "copyright", "credits" or "license" for more information.
>>> import pymc3 as pm
>>> model = pm.Model()
>>> with model:
...     a = pm.Normal('a', mu=0, sd=1)
...     db_sqllite = pm.backends.SQLite("test-sqllite")
...     trace = pm.sample(1000, n_init=1000, trace=db_sqllite, cores=1, njobs=1)
...
Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
Sequential sampling (2 chains in 1 job)
NUTS: [a]
100%|████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 1500/1500 [00:00<00:00, 2328.49it/s]
100%|████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 1500/1500 [00:00<00:00, 2453.27it/s]
Traceback (most recent call last):
  File "<stdin>", line 4, in <module>
  File "/Users/orion.delwaterman/.pyenv/versions/3.6.1/envs/hiring-horizons/lib/python3.6/site-packages/pymc3/sampling.py", line 439, in sample
    trace = _sample_many(**sample_args)
  File "/Users/orion.delwaterman/.pyenv/versions/3.6.1/envs/hiring-horizons/lib/python3.6/site-packages/pymc3/sampling.py", line 494, in _sample_many
    return MultiTrace(traces)
  File "/Users/orion.delwaterman/.pyenv/versions/3.6.1/envs/hiring-horizons/lib/python3.6/site-packages/pymc3/backends/base.py", line 265, in __init__
    raise ValueError("Chains are not unique.")
ValueError: Chains are not unique.

HDF5 Fails:

16:55 $ python
Python 3.6.1 (default, Nov  8 2017, 14:29:33)
[GCC 4.2.1 Compatible Apple LLVM 7.0.2 (clang-700.1.81)] on darwin
Type "help", "copyright", "credits" or "license" for more information.
>>> import pymc3 as pm
>>> model = pm.Model()
>>> with model:
...     a = pm.Normal('a', mu=0, sd=1)
...     db = pm.backends.HDF5('test-hdf5-3')
...     trace = pm.sample(1000, n_init=1000, trace=db, cores=1, njobs=1)
...
Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
Sequential sampling (2 chains in 1 job)
NUTS: [a]
100%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 1500/1500 [00:07<00:00, 193.66it/s]
  0%|                                                                                                                                                     | 0/1500 [00:00<?, ?it/s]
Traceback (most recent call last):
  File "<stdin>", line 4, in <module>
  File "/Users/orion.delwaterman/.pyenv/versions/3.6.1/envs/hiring-horizons/lib/python3.6/site-packages/pymc3/sampling.py", line 439, in sample
    trace = _sample_many(**sample_args)
  File "/Users/orion.delwaterman/.pyenv/versions/3.6.1/envs/hiring-horizons/lib/python3.6/site-packages/pymc3/sampling.py", line 482, in _sample_many
    step=step, random_seed=random_seed[i], **kwargs)
  File "/Users/orion.delwaterman/.pyenv/versions/3.6.1/envs/hiring-horizons/lib/python3.6/site-packages/pymc3/sampling.py", line 526, in _sample
    for it, strace in enumerate(sampling):
  File "/Users/orion.delwaterman/.pyenv/versions/3.6.1/envs/hiring-horizons/lib/python3.6/site-packages/tqdm/_tqdm.py", line 862, in __iter__
    for obj in iterable:
  File "/Users/orion.delwaterman/.pyenv/versions/3.6.1/envs/hiring-horizons/lib/python3.6/site-packages/pymc3/sampling.py", line 614, in _iter_sample
    strace.setup(draws, chain, step.stats_dtypes)
  File "/Users/orion.delwaterman/.pyenv/versions/3.6.1/envs/hiring-horizons/lib/python3.6/site-packages/pymc3/backends/hdf5.py", line 154, in setup
    self._set_sampler_vars(sampler_vars)
  File "/Users/orion.delwaterman/.pyenv/versions/3.6.1/envs/hiring-horizons/lib/python3.6/site-packages/pymc3/backends/base.py", line 80, in _set_sampler_vars
    raise ValueError("Can't change sampler_vars")
ValueError: Can't change sampler_vars
>>> model.unobserved_RVs
[a]

When I tried using a single chain it worked correctly for all three.

Versions and main components

• PyMC3 Version: 3.3
• Theano Version: 1.0.1
• Python Version: 3.6.1
• Operating system: OS X
• How did you install PyMC3: (conda/pip) pip

[junpenglao]

Thanks for reporting @delwaterman.
Change title as it happens with cores=1. Setting random_seed does not help.

[shkr]

@junpenglao Yeah. Same as the other issue #2801 copying the trace when running on 1 core for every new chain will fix this error.

The same trace is being reused (which can be confirmed by print([id(c) for c in traces])) for both the "chains" when only 1 job is used. So the a duplicate chain error is being thrown.

Again just like the other issue, this will fix the error.

https://github.com/pymc-devs/pymc3/blob/master/pymc3/sampling.py#L501

def _sample_many(draws, chain, chains, start, random_seed, step, **kwargs):
    from copy import deepcopy
    traces = []
    for i in range(chains):
        kwargs['trace'] = deepcopy(kwargs['trace'])
        trace = _sample(draws=draws, chain=chain + i, start=start[i],
                        step=step, random_seed=random_seed[i], **kwargs)
        if trace is None:
            if len(traces) == 0:
...

[shkr]

In this case, since it is hard to see a workaround just creating a new copy of the trace-backend is made, given you do want the chains even when not created in parallel to be independent of each other.

[junpenglao]

Yeah, deepcopy will fix the issue but it will also give the performance a big hit.
Would it be appropriate to do: if sequential sampling is called, we keep sampling but reshape at the end to make it into multichain? And we do deepcopy for SGFS sampler in sequential sampling only.

[shkr]

Okay. I will ponder a bit more on alternates to deepcopy if possible for SGFS and subsequently update the other #2801 PR

One question about this issue. For folding one trace backend into two or more, would require modification to the BaseTrace class. Is that what you are referring to ? Just want to think about the expected change in the code.

[junpenglao]

Actually, probably not a good idea to folding one long trace into more, see this discussion.

[shkr]

@junpenglao Isn't the initialize point the same for all chains in pymc3 ?

The reason you want multiple chains started from diverse points is to check if your chain is behaving well. (This is what Rhat does) The (very small) efficiency loss is worth it to check you’re not computing rubbish.

that's the reason for not doing one long chain in the thread

[junpenglao]

Isn't the initialize point the same for all chains in pymc3 ?

The default 'jitter+adapt_diag' gives different starting point, but I think in general the starting point is always the model.test_point. It would be a good idea to move the jitter part into pm.sample actually.

[shkr]

Okay, makes sense. I think I can start a PR with a 'jitter' by default in pm.sample, and then follow up with another PR that fixes this issue and Issue handled in PR #2801 for njobs=1 and chains > 1 situation. Let me know your thoughts, and whether that plan makes sense.

[junpenglao]

That sounds good to me!

[junpenglao]

Also, just some related ideas:
I think the PR should work towards the idea of creating a function that does "forward random sample from prior" (something like pm.sampling.sample_prior?). Currently, we have most of the components:

• in smc there is a part of code sample from the prior:
  https://github.com/pymc-devs/pymc3/blob/801accb5f236ab9daa89a8fcd9d09a3ba4ed0a39/pymc3/step_methods/smc.py#L186-L193
  But there are some cases that this doesn't work: DensityDist RV with no random method, or a distribution without random method. In those cases, use the part from below:
• the jitter in https://github.com/pymc-devs/pymc3/blob/801accb5f236ab9daa89a8fcd9d09a3ba4ed0a39/pymc3/sampling.py#L1291-L1295
• the jitter failed sometimes as well that the uniform random actually takes some parameter outside of the support (it is generally fine in NUTS, as the default support is (-inf, inf)), there should be some failback options for these cases.