Question. Why not return a fresh task every run?

[sglyon]

I am working on writing a generic Gibbs sampler that shares design principles with the code here. (the WIP is on the gibbs branch of my fork)

I was wondering why we don't always return a fresh task after running a chain? Specifically, I noticed that a task is stopped before being returned in parallel (run_serialmc_exit), but not if you run in serial. Unless I'm mistaken this is just a way to tell Julia that it is ok to serialize the task so it can be returned from another process.

When I wrote a simple random walk metropolis routine I found that it was easier to always spin a new task at the end of running a chain (in serial or in parallel) because that way the task is always ready to resume (again, in serial or in parallel). I found that I could get the task to resume at the same place without losing any samples simply by consuming the task one more time to find a new model.init value, spinning a new task with that updated model, then making sure that the first produce in the sampler was model.init.

Is there a reason we don't do that here?

[papamarkou]

The reason I had set the task to stop in the parallel case is because, at least at that time, it was not possible to serialize tasks (I haven't checked Julia's base since then, so it may be possible by now to serialize them).

I am not sure whether you refer to the returned chain, which holds a task field; if so, the task field will be dropped altogether from the returned chain, see my fork of MCMC.jl, where the whole type system for MCMC has been altered quite a lot. Although the core changes have been achieved in the updated implementation of my fork, I haven't merged them because I haven't ported all the samplers yet.

Perhaps I could possibly disable some existing samplers in order to merge, so that everyone can be on the same page with the new MCMC interface and code (for ex there is not much point for you to put this effort on the basis of an old-to be redundant type system, you'd rather use the new one to write up the Gibbs sampler you want to contribute).

@fredo-dedup, would it be a sensible senario to merge the changes I have made disabling temporarily some samplers (not the most common ones that are ready pretty much)?

[sglyon]

@scidom wow, thanks for the heads up. I just had a look at your fork and looks like you have done a lot of work refactoring.

I'd love to base my Gibbs sampler off the codebase in your brach. Should I pull your master and rebase my work off it? Or, should we wait for @fredo-dedup to see if we can get your version merged before I work off of it?

[papamarkou]

Yes, from my side feel free to pull off my master to work on Gibbs, in fact I am happy if you find the new (inter)face of MCMC cleaner and more practical. In any case we will merge the changes in due time, so I think you are safe to base your work on my fork.

Let's see what @fredo-dedup thinks about a "premature" merge of my fork to make sure all involved Julian parties are happy :) I should finish the porting very soon, was held a bit behind due to a small health issue with eyes due to excessive computer use.

[papamarkou]

With a bit of an extra self-push, I ported the acceptance rejection sampler, so I got 8 out of the 11 samplers ported; I now have 3 left, namely the important NUTS sampler and the less important in widely-used applications RMHMC and MMALA (I am not sure we really need these 2 right now but for the sake of completeness I will port them too). I promise I will get everything done by the coming Friday so that we can merge.

[fredo-dedup]

Hi @scidom,
it looks like you've ported the most important samplers so I have no issue with a merge with master (with or without a new package version as you see fit).
I am not clear about the old task vs fresh task question though( but don't hold anything because of me). The task had the advantage of being a container for the state of the sampler when there was some adaptation algo involved (RAM, NUTS, ....). But there are probably better solutions for that. Clearly, I need to catch up on your new work !

[papamarkou]

Yeap, tasks are very helpful as you introduced them in our package more than a year ago, and I would never drop them! Basically, what I had done on that respect was to specify Job types, such as tasked and plain jobs; this gives fine-grained control allowing the user to choose whether he wants to run his MCMC with or without tasks.

Since I have your approval, I will then merge the code; I will do this by tomorrow the latest, in order to make some more small additions. After this, I will write up elaborate "readthedocs" documentation, so that we can all be on the same page with the package updates.

[goedman]

Hi @scidom I'd been trying if I could get scidom/MCMC.jl to run for the slice sampler, but that doesn't work yet on my system:

julia> include("/Users/rob/.julia/v0.3/MCMC/test/test_slice_sampler.jl")
ERROR: unrecognized keyword argument "steps"
 in include at /Applications/Julia-0.3.2.app/Contents/Resources/julia/lib/julia/sys.dylib
 in include_from_node1 at /Applications/Julia-0.3.2.app/Contents/Resources/julia/lib/julia/sys.dylib
while loading /Users/rob/.julia/v0.3/MCMC/test/test_slice_sampler.jl, in expression starting on line 19

This is on Julia v0.3.2. I'll try again after you merge the changes.

As of this morning ARS1 and ARS2 indeed work fine, you beat me to that! But still, a good way to study your improvements.

julia> Pkg.test("MCMC")
INFO: Testing MCMC
Running tests:
  * test_empmctuner.jl *
    Testing basic EmpiricalMCTuner constructors...
    Testing that EmpiricalMCTuner tuners works with all samplers...
Burnin iteration 100 of 1000: 77.0 % acceptance rate
Burnin iteration 200 of 1000: 83.0 % acceptance rate
...
Burnin iteration 1000 of 1000: 86.0 % acceptance rate
  * test_ARS1.jl *

M: [10.0 10.0]

Acceptance rate: 30.307743584046218
Parameter 1
Min        -2.7485353323597135
Mean       -0.07148964097150934
Max        3.1297758086060483
MC Error   0.0576240046952458
ESS        369.6677390517674
AC Time    24.348892394798675

Acceptance rate: 31.407621375402734
Parameter 1
Min        -2.882154115906877
Mean       0.04760503880756451
Max        3.260958543198918
MC Error   0.05031067651408723
ESS        445.5541614003709
AC Time    20.201808847907454

  * test_ARS2.jl *

M: 120.0

Acceptance rate g0: 11.209865570492168

Parameter 1
Min        -1.1490068050198476
Mean       -0.012927468168042612
Max        1.3299385490286677
MC Error   0.021511701242054978
ESS        362.0319829549099
AC Time    24.862444269519276

Parameter 2
Min        0.8513403757360513
Mean       2.0289615677895063
Max        3.134343916078227
MC Error   0.01904927156690838
ESS        415.7464965061753
AC Time    21.650212510849876

INFO: MCMC tests passed
INFO: No packages to install, update or remove

[papamarkou]

Cheers @goedman, yes, I know, I haven't updated most of the tests to reflect the syntax updates. I updated the slice sampler test and pushed the relevant changes to my fork, you may try it now, it should work.

I will do a little more work tonight before I merge everything to our JuliaStats MCMC package.

[papamarkou]

Just to keep you up to speed informally, I dropped the model specification syntax based on * operator overloading. There are 3 alternatives

• the old run(model, sampler, runner) (in fact there are a couple of hidden optional arguments there, the whole command is run(model, sampler, runner, tuner, job)),
• the new run(MCSystem(model, sampler, runner)) and
• the simplest (new) interface, which you can find in src/ui.jl in my fork. This hides the creation of model, sampler and runner under the hook, as there has been an interest to provide a simpler way of interacting with the package.

It is still possible to specify a bunch of MCMC simulations; this can be done by using the MCSystem syntax instead of the older one based on the * opearator. I will update the README to shortly explain all these.

[goedman]

Thanks @scidom . Works fine now!

[papamarkou]

You're welcome @goedman.

[papamarkou]

I merged the major updates. The test failure is due to the separate recent issue with the ArrayViews.jl package. I still need to port more code, but all the basic functionality is available. I didn't want to hinder progress of any Julia user or developer who wants to contribute to MCMC, this is why I pushed this forward without waiting for the whole codebase to be ported. @spencerlyon2 you should be able now to proceed with the Gibbs sampler coding.

[sglyon]

@scidom

This is great. Thanks for pushing this through so quickly. I will base my Gibbs code on the current master.

I will probably submit a rfc soon. I have made some design decisions that makes the Gibbs model feel/work different than the likelihood model we currently have and am open to any suggestions that could make my code fit more harmoniously with the likelihood based code.

Thanks again and you'll hear from me shortly.

[sglyon]

Closed as no longer as relevant with new api. If we feel this issue should be addressed again under new code base at can reopen

[papamarkou]

Sure, @spencerlyon2, I would be very happy to see your PR for the Gibbs sampler and how you decided to software-engineer it. I agree that the existing likelihood-based model is somewhat orthogonal to the needs of a Gibbs sampling scheme, so we are certainly open to new ideas, hoping that their fusion will amplify the potential of the package.

[papamarkou]

@spencerlyon2 i am sorry about this, but I needed to make yet a minor change in the code base. The MCSystem and MCJob types were highly overlapping. To avoid any redundancy, I dropped the MCSystem by collapsing these two into a single job type (MCJob). This has some minor benefits:

• the MCStash subtypes are now a field in the job types. This does not change the way the stash is handled at a lower level (it is still allocated on the heap), but it is much cleaner coding to have explicit access to it as a job type field rather than manipulating it indirectly,
• people are accustomed to running jobs rather than systems (a very usual example here is job schedulers), so the latest update fits existing terminology.

I promise that the dust has settled and you won't be seeing any other major code changes for a while.

[sglyon]

@scidom, no problem. I was about to submit a PR, but I'll adjust my code before submitting.

When I do open the PR it will be at a 70% done level and I'm hoping we can iterate together on converging to the optimal way to describe a Gibbs sampler.

[sglyon]

Also I have some questions. Is there an easier way to chat about some design things? Perhaps google hangouts or a skype chat. I'm thinking IM rather than video call, but it would be helpful to me to talk in a bit more interactive way (that also doesn't generate so many notifications for others).

If you have some time and access to a technology like that I think it would make my process for learning the new api a bit smoother.

[papamarkou]

Thanks very much for your patience with the ongoing updates the last few days, it is appreciated. I am sure that we will manage to get the Gibbs sampler in a good shape after discussing over your forthcoming PR.

Sure, you may reach me at theodore dot papamarkou at gmail dot com if you would like to use Google hangout. I have up to 30 mins free now as we will be spending time with my wife in a bit; if you are free now for a short chat, that would be great, otherwise we could arrange it tomorrow evening.