Multiple improvements to pymcmodels

[jchodera]

• Add support for specified F_PL and its uncertainty dF_PL
• Introduce LogNormalWrapper to make lognormal priors more convenient (Add LogNormalWrapper for pymc2 #84)
• Eliminate normal prior for concentrations (leaving lognormal)
• Convert strictly positive quantities to lognormal priors for robustness (Convert all strictly positive quantities to log-space for pymc2 models #85)
• Separate handling of [L] = 0 wells (Allow handling of [L] = 0 data. #81)
• Add rudimentary nosetest to make sure I haven't broken things completely
• Decouple buffer background fluorescence for +/- protein rows

I didn't have time to update any examples to try specifying the F_PL and its uncertainty dF_PL as optional arguments to pymcmodels.make_model(), but feel free to add to this branch/PR directly as you test it out, @sonyahanson!

[jchodera]

@sonyahanson : I'm having trouble figuring out how to test any of these improvements.

What can I run to actually test things? I can't find anything in examples/direct-fluorescence-assay that will actually generate some output to test, and the Jupyter notebooks I've found (like 3a Bayesian fit xml file - SrcGefitinib) don't seem to actually use make_model from pymcmodels.py.

[jchodera]

@sonyahanson : Regarding [L]=0 support for +/- protein: Do we want to support only a single well for [L]=0 with protein and a single well without protein, or do you want to be able to add an arbitrary number of [L]=0 control wells with or without protein?

[sonyahanson]

You need to make sure the xml file path in inputs_p38_singlet actually points to your xml files. I can make the behavior here, more intuitive. (e.g. throw an error saying it didn't find the xml files).

We definitely need more/any tests, and better documentation on the examples. This is on my to do list to work on in the next week anyway. Let me know of suggestions.

Regarding [L]=0, whichever is easiest for now. A single well is perfectly fine for now, but maybe in the future we would want an arbitrary number, but maybe that future will see a pretty large overhaul of the API anyway...

[jchodera]

You need to make sure the xml file path in inputs_p38_singlet actually points to your xml files. I can make the behavior here, more intuitive. (e.g. throw an error saying it didn't find the xml files).

That would be super useful!

We definitely need more/any tests, and better documentation on the examples. This is on my to do list to work on in the next week anyway. Let me know of suggestions.

Priorities for tests are probably:

• Does it run at all with python 2/3?
• Does it produce correct results for synthetic data?
• Can it process a variety of XML files (with short runs) without dying?

Regarding [L]=0, whichever is easiest for now. A single well is perfectly fine for now, but maybe in the future we would want an arbitrary number, but maybe that future will see a pretty large overhaul of the API anyway...

The difference in API between one and multiple wells is very little, but it does change the implementation under the hood.

[jchodera]

OK, things seem to work locally now, so I'm going to finish the [L]=0 case and then add a very rudimentary test just to make sure things run.

[jchodera]

I think I was going about the [L]=0 case the wrong way by trying to break it out as a separate set of five things to pass into the API. Instead, I think I can just detect which elements of Lstated are zero and handle those as special cases internally, which doesn't require any changes to the API at all. Trying this now.

[jchodera]

I've implemented the scheme I suggested above for [L]=0. We autodetect when either Lstated or Pstated have zero concentration entries and deal with those appropriately.

I've also added a quickmodel run as a travis test.

This PR should now be complete.

[jchodera]

@sonyahanson : Review and merge when ready!

[gregoryross]

It's great that quickmodel has been added to travis. However, quickmodel runs the default number of MCMC steps, currently set as 20000 PyMC moves, which may take a while on travis. How about we augment the argparser on quickmodel's so that we can specify far fewer moves on travis? Parsing the number of MCMC moves to quickmodel will make it easier to use as well.

[gregoryross]

This all looks good to me! My comment on quickmodel is minor, and could be added in another PR if we want to use these changes asap.

[sonyahanson]

I'm going to work through this now. One thing that pops out immediately: why were the defaults for run_mcmc changed from nthin=50, nburn=500, niter=1000 to nthin=50, nburn=0, niter=20000?

I agree with Greg that travis should be run with fewer iterations.

[jchodera]

One thing that pops out immediately: why were the defaults for run_mcmc changed from nthin=50, nburn=500, niter=1000 to nthin=50, nburn=0, niter=20000?

The output of quickmodel was absolute garbage otherwise. Runtime was ~1 min/dataset on my machine, and total processing time for the quickmodel example was < 10 min, which seems reasonable for travis.

I'm not certain why travis isn't running, however. We might have to fix that in a different branch.

I like the idea of adding the command-line option to quickmodel, but think the new defaults are sensible. I'll add the command-line argument idea as an issue.

[sonyahanson]

https://github.com/choderalab/assaytools/blob/master/scripts/quickmodel.py#L301

[sonyahanson]

I will be careful to use the same nthin,niter,nburn to compare the master branch to this updated 'improved' branch just to make sure the comparisons are fair.

[jchodera]

Looks like you're discarding the initial non-equilibrated DeltaG values, but not discarding the initial non-equilibrated traces when plotting:
https://github.com/choderalab/assaytools/blob/master/scripts/quickmodel.py#L129-L134

Can I fix that too?

[sonyahanson]

Yes, that's correct, feel free to fix.

[sonyahanson]

Maybe include but colored differently?

[jchodera]

I will be careful to use the same nthin,niter,nburn to compare the master branch to this updated 'improved' branch just to make sure the comparisons are fair.

I totally understand. How about this: The current settings are configured to do the same amount of sampling. Maybe run the data and compare it to what you have printed out?

[sonyahanson]

on it

[jchodera]

Here's the Bosutinib-CD example from above:

[jchodera]

The pre-equilibration traces are still shown, just very lightly shaded. They're mostly lying underneath the darker traces, but you can see them a little bit in the salmon traces here:

[jchodera]

I'm moving on to other projects now, so I'll let you folks take the PR over from here.

I can tackle the outlier detection in a separate PR when I have time to cycle back.

[sonyahanson]

yeah, I think this is fine for now, but don't think we should merge quite yet.

[jchodera]

I've just fixed a bug in the code where the same ligand concentration was being used for both the +protein and -protein rows. This prevented the true ligand concentrations from adjusting to deviations from the expected binding curve even when dispensing errors occurred.

I've also temporarily disabled the Metropolis step methods, and am testing how things work without them.

[sonyahanson]

So it seems like the big problem was allowing the metropolis tuning throughout, since we've now taken out the burnin step.

But making this changed improved our overlap between repeats from this:

To this:

Merge away. We can address further problems in a future PR.

[sonyahanson]

This change fcde1cf doesn't quite match the commit message, what was the motivation? I'm curious if this also effected the change in our results...

[jchodera]

Whoops---there were two changes in rapid succession:

• Tune throughout (this was the major remedy)
• Also adjust the scale of the AdaptiveMetropolis method that makes joint proposals between F_PL and DeltaG. This was not being tuned before, but once tuning was working, the exact choice became less relevant. I changed it to 0.1 to match all the other Metropolis step methods

[jchodera]

Can we merge this and forge ahead?

[jchodera]

Merging this so I can forge ahead.