New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

SMC #1923

Closed
wants to merge 14 commits into
base: master
from

Conversation

Projects
None yet
9 participants
@hvasbath
Copy link
Contributor

hvasbath commented Mar 20, 2017

Sequential Monte Carlo - renamed from ATMCMC here:
#1569

Sampling speed improvements, rebased on master ...

This was referenced Mar 20, 2017

@hvasbath

This comment has been minimized.

Copy link
Contributor

hvasbath commented Mar 20, 2017

One test is failing I have no clue why. What is different there to the others ...?

return x[(n_steps - 1)::n_steps]


def two_gaussians(x):

This comment has been minimized.

@twiecki

twiecki Mar 20, 2017

Member

This can be simplified to:

pm.NormalMixture('two_gaussians', mu=[mu1, mu2], sigma=[dsigma, isigma], w=[w1, w2])

This comment has been minimized.

@hvasbath

hvasbath Mar 20, 2017

Contributor

Will try to do that, last time I tried it was absolutely not straight forward, as mu is n-dimensional .

This comment has been minimized.

@twiecki

twiecki Mar 20, 2017

Member

Oh, I missed that. So mu1 and mu2 are not scalars? In that case it might be too much trouble. @AustinRochford Do you think that should be possible?

This comment has been minimized.

@AustinRochford

AustinRochford Mar 20, 2017

Member

@twiecki yes, right now Mixture only supports 1d mixtures. There are some subtleties in supporting multidimensional distributions and weights that vary with observed data points (necessary for DDR, for example) that I am thinking about but have not resolved.

@twiecki

This comment has been minimized.

Copy link
Member

twiecki commented Mar 20, 2017

@hvasbath Can you convert the example to an IPython NB, with description of the what and why?

@twiecki

This comment has been minimized.

Copy link
Member

twiecki commented Mar 20, 2017

======================================================================
ERROR: test_sample (pymc3.tests.test_smc.TestSMC)
----------------------------------------------------------------------
Traceback (most recent call last):
  File "/home/travis/build/pymc-devs/pymc3/pymc3/tests/test_smc.py", line 72, in test_sample
    rm_flag=False)
  File "/home/travis/build/pymc-devs/pymc3/pymc3/step_methods/smc.py", line 691, in ATMIP_sample
    step.select_end_points(mtrace)
  File "/home/travis/build/pymc-devs/pymc3/pymc3/step_methods/smc.py", line 360, in select_end_points
    n_steps = len(mtrace)
  File "/home/travis/build/pymc-devs/pymc3/pymc3/backends/base.py", line 316, in __len__
    chain = self.chains[-1]
IndexError: list index out of range

Looks like a python 3 issue. Perhaps a iterator is used where you expect a list?

@hvasbath

This comment has been minimized.

Copy link
Contributor

hvasbath commented Mar 20, 2017

It looks like the sampling hasnt been executed, as the mtrace is empty.
So there is maybe a problem with the atext.paripool function regarding python3? Can you/someone look at that function if you see a python 3 problem there?

I have absolutely no experience with notebooks. They cannot be opened with a texteditor.
What do I have to do to create a notebook? If it is difficult with system crashing installations to python3 etc I will likely have no time for that...

@hvasbath

This comment has been minimized.

Copy link
Contributor

hvasbath commented Mar 20, 2017

It can also be that the outfolder that is defined in the example causes the problem, because it is being created in the execution path? Maybe that crashes Travis? What do you think? Shall I also use temporal directories here? @twiecki

import theano.tensor as tt
from matplotlib import pyplot as plt

test_folder = ('ATMIP_TEST')

This comment has been minimized.

@hvasbath

hvasbath Mar 20, 2017

Contributor

This test folder here ...

@twiecki

This comment has been minimized.

Copy link
Member

twiecki commented Mar 20, 2017

Shall I also use temporal directories here?

Definitely, there's nice python API for that, I think.

@hvasbath

This comment has been minimized.

Copy link
Contributor

hvasbath commented Mar 20, 2017

Ok created the notebook and removed the example.

@twiecki

This comment has been minimized.

Copy link
Member

twiecki commented Mar 20, 2017

Seems like NB did not run to completion.

@hvasbath

This comment has been minimized.

Copy link
Contributor

hvasbath commented Mar 20, 2017

I dont know why the figure is not being displayed in the notebook. You used exactly the same syntax in the getting started case... Must be some local problem on my computer I guess.

K figured it out needs %matplotlib inline command

@hvasbath

This comment has been minimized.

Copy link
Contributor

hvasbath commented Mar 20, 2017

Ok it definitely was not the temporal directory ... So we need a python3 expert here ...

@twiecki

This comment has been minimized.

Copy link
Member

twiecki commented Mar 21, 2017

@hvasbath The trick is to use %matplotlib inline as the first cell. Then you don't need the plt.shows.

"pm.traceplot(mtrace, transform=last_sample, combined=True);\n"
"axs = pm.traceplot(mtrace, transform=last_sample, combined=True)\n",
"\n",
"plt.show()"

This comment has been minimized.

@twiecki

twiecki Mar 21, 2017

Member

shouldn't need this.

@hvasbath

This comment has been minimized.

Copy link
Contributor

hvasbath commented Mar 21, 2017

Ok did a python3 installation on another computer and was playing around. Somehow the work list that is being created in smc.py line. 857 -860 has a length of zero. I have no clue why, probably because of some changes with zip or maybe the map that has an influence on creating the progressbar lists.
Then in the smc_text backend in the paripool function there is a usage of map and apparently this function has been changed a lot from python2 to 3.
To conclude it will be absolutely not straight forward to get it running for both python versions simultaneously with the same code. For that I am lacking too much knowledge about python3 and I wont have time looking into it for the next monthes ...

@twiecki

This comment has been minimized.

Copy link
Member

twiecki commented Mar 21, 2017

@hvasbath

This comment has been minimized.

Copy link
Contributor

hvasbath commented Mar 21, 2017

cannot import name 'zip'

@twiecki

This comment has been minimized.

Copy link
Member

twiecki commented Mar 21, 2017

Sorry, from six.moves import zip

I think you can do the same with map.

@hvasbath

This comment has been minimized.

Copy link
Contributor

hvasbath commented Mar 21, 2017

Thanks for the hints @twiecki !Just tried it, also putting list() around the map. Now it runs with njobs=1, but still not for n>1
Somehow it doesnt run through the generator. Only does the first iteration and then leaves the loop.

@hvasbath

This comment has been minimized.

Copy link
Contributor

hvasbath commented Apr 6, 2017

Sorry for being unresponsive. I just returned from a vacation. I am glad it also runs for you @ColCarroll within reasonable times! Putting an EXPERIMENTAL label on it is perfectly fine for me.

I am sorry for the amount of code, I would have been happy if it was less, probably it could have been. But I am programming in python only since 1 1/2 years so I guess there might be a lot of room for improvement in terms of code style and factorization. Anywhere a "proper" programmer by profession?

I am keen to get your remarks @aloctavodia !

@twiecki

This comment has been minimized.

Copy link
Member

twiecki commented Apr 6, 2017

We should also test how well SMC does. It should do pretty poorly in high dimensions, but perhaps for mixture models it works well?

@aloctavodia

This comment has been minimized.

Copy link
Member

aloctavodia commented Apr 6, 2017

@twiecki I guess it should perform the same as Metropolis, except that it should do better for multimodal posteriors, right?. One possible extension of SMC could be to use NUTS instead of Metropolis.

@hvasbath I am not a "proper" programmer, but I would like to help, from improving the code to testing the performance/suitability to different scenarios.

How should I proceed to help? Should I wait for the merge under "experimental conditions" ?

@ColCarroll

This comment has been minimized.

Copy link
Member

ColCarroll commented Apr 6, 2017

I can make a PR today on @hvasbath's repo adding warnings (other option is to develop on a branch, which might make sense here since it is quite decoupled).

Any good references to read up on this?

@hvasbath

This comment has been minimized.

Copy link
Contributor

hvasbath commented Apr 6, 2017

@twiecki No it does very well so far under any circumstances, so far I had no rpoblem with any model. It does especially well in high dimensions thats the whole point of it. It is just a matter of number of chains you use to sample. The higher the dimensions the higher the number of chains should be defined.
There are two references given in the code, another one probably the most original one:
http://www.stats.ox.ac.uk/~doucet/delmoral_doucet_jasra_sequentialmontecarlosamplersJRSSB.pdf

@aseyboldt

This comment has been minimized.

Copy link
Member

aseyboldt commented Apr 6, 2017

You could also use the fixtures in sampler_fixtures.py to test if posterior samples for a couple of distributions are consistent with what we expect. This helped me a lot when I was working on NUTS. The tests seemed to detect even relatively small biases. For this you need to add a new class in sampler_fixtures like BaseSampler and add a couple of tests that use this in test_posteriors.py

@fonnesbeck

This comment has been minimized.

Copy link
Member

fonnesbeck commented Apr 9, 2017

Perhaps rather than a EXPERIMENTAL tag, we could adopt a pymc3.sandbox submodule as a staging area for algorithms that we would like merged but are not yet ready for use in production.

@twiecki

This comment has been minimized.

Copy link
Member

twiecki commented Apr 10, 2017

@fonnesbeck I like the idea in principle but worry about name-space issues. This PR, for example, adds code in various parts, including the backends, would those live in the sandbox too?

@ferrine

This comment has been minimized.

Copy link
Member

ferrine commented Apr 10, 2017

I like sandbox idea, all changes outside it should be just reviewed and approved by pymc-team. That points on our not perfect enough inner API or architecture that is not extensible sometimes.

@fonnesbeck

This comment has been minimized.

Copy link
Member

fonnesbeck commented Apr 10, 2017

The sandbox would be more for methodologically experimental code, rather than for infrastructure like backends, I should think. More complicated changes would have to remain in their own PR until they are ready.

@hvasbath

This comment has been minimized.

Copy link
Contributor

hvasbath commented Apr 13, 2017

I like the sandbox idea! Somehow I am lost, is there anything that I should do next?

@aloctavodia

This comment has been minimized.

Copy link
Member

aloctavodia commented Apr 13, 2017

I have found another reference that could be useful to read. I just gave it a quick look, but it seems like an improved version of the Del Moral's paper.
https://arxiv.org/abs/1504.05753, http://ieeexplore.ieee.org/document/7339702/

BTW, I have been confused all this time thinking ATMCMC/SMC was a variant of Replica Exchange/Parallel Tempering, but they are really different algorithms!

@twiecki

This comment has been minimized.

Copy link
Member

twiecki commented Apr 13, 2017

@hvasbath Seems like there are still failing tests, have you looked at those?

@hvasbath

This comment has been minimized.

Copy link
Contributor

hvasbath commented Apr 13, 2017

They are all running. As @ColCarroll stated somehow in travis the tests take much longer then when run locally. So thats a problem in the travis setup, which cannot take care of. @aloctavodia thanks for the reference I will check it out!

@hvasbath

This comment has been minimized.

Copy link
Contributor

hvasbath commented Apr 13, 2017

Oh that paper seems really cool. I always thought what a waste to keep only the last samples! They mainly improved on that. Thats great! However, I wont be able to implement that soon as I have to focus on writing my articles. But I will do it at some point- if no one else did take care of it in the meanwhile ...

@twiecki

This comment has been minimized.

Copy link
Member

twiecki commented Apr 13, 2017

I think this PR is close to being merged. I would add an experimental warning (we can always move it to a sandbox separately before 3.1). However, the timeouts in the unittest are a blocker, not sure how we can resolve that. Notably, it only happens in python 2.

@ColCarroll

This comment has been minimized.

Copy link
Member

ColCarroll commented Apr 13, 2017

@ColCarroll

This comment has been minimized.

Copy link
Member

ColCarroll commented Apr 13, 2017

Super strange. I ran the file run_test.py:

import pytest

if __name__ == '__main__':
    pytest.main(['-vx', 'pymc3/tests/test_smc.py'])

Both python run_test.py and python -m cProfile -o time.out run_test.py take ~6 minutes on my machine, but cprofilev run_test.py takes ~80seconds (cprofilev). Luckily, both give a profile output, and the big difference is:

  • normal run:
ncalls  tottime  percall  cumtime  percall filename:lineno(function)
11577  300.409    0.026  300.409    0.026 {method 'acquire' of 'thread.lock' objects}
  • cprofilev run:
   ncalls  tottime  percall  cumtime  percall filename:lineno(function)
     131    0.000    0.000    0.000    0.000 {method 'acquire' of 'thread.lock' objects}

I don't know what this means yet, but it accounts for all of the timing difference. Will check some more after work.

@hvasbath

This comment has been minimized.

Copy link
Contributor

hvasbath commented Apr 13, 2017

Wow! Thanks a lot for investigating that! I thought we use no threading as it is mutliprocessing it should be forking? Then I wouldnt be surprised. Is that the GIL he's wainting for, because he uses threading instead of forking?
In my pre-restructuring version anything to python3 it also has this short runtime ... Somehow joblib does something we dont want ;) ...

@junpenglao

This comment has been minimized.

Copy link
Member

junpenglao commented Apr 16, 2017

Great work @hvasbath ! Just wondering what is the status on this? I have a model which doesnt work very well in NUTS (the geometry is not continuous) and using Metropolis returns broadcast error #1983. I am curious to try SMC on it.

@ColCarroll

This comment has been minimized.

Copy link
Member

ColCarroll commented Apr 16, 2017

@ColCarroll ColCarroll referenced this pull request Apr 16, 2017

Merged

SMC experimental #2045

@ColCarroll

This comment has been minimized.

Copy link
Member

ColCarroll commented Apr 17, 2017

#2045 was merged with these changes -- thanks again @hvasbath !

@ColCarroll ColCarroll closed this Apr 17, 2017

@twiecki

This comment has been minimized.

Copy link
Member

twiecki commented Apr 17, 2017

@hvasbath

This comment has been minimized.

Copy link
Contributor

hvasbath commented Apr 17, 2017

Great finally done! Thanks for helping to finish including this to pymc3! Although it is kind of sad that this didnt increase my contribution stats and it went all to @ColCarrol. I know this is kind of nerdy but I like these stats ;) .

@twiecki

This comment has been minimized.

Copy link
Member

twiecki commented Apr 18, 2017

@twiecki

This comment has been minimized.

Copy link
Member

twiecki commented Apr 18, 2017

@hvasbath OK, the cleanest solution was to revert and then re-commit as you: f4803ce Sorry about that.

@hvasbath

This comment has been minimized.

Copy link
Contributor

hvasbath commented Apr 18, 2017

Wow thanks @twiecki I didnt know such a thing is possible. Cool!

@twiecki

This comment has been minimized.

Copy link
Member

twiecki commented Apr 18, 2017

git commit --amend --author='...' ;)

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment