Multiple improvements to pymcmodels

[jchodera]

• [x] Add support for specified F_PL and its uncertainty dF_PL
• [x] Introduce LogNormalWrapper to make lognormal priors more convenient (#84)
• [x] Eliminate normal prior for concentrations (leaving lognormal)
• [x] Convert strictly positive quantities to lognormal priors for robustness (#85)
• [x] Separate handling of [L] = 0 wells (#81)
• [x] Add rudimentary nosetest to make sure I haven't broken things completely
• [x] 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!

[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.

[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?

More specifically:

• nthin seemed fine, so I didn't change that; it could possibly be reduced
• nburn > 0 is unnecessary now that we use automatic equilibration detection
• niter=1000 gave so few uncorrelated samples that the results were not meaningful

[sonyahanson]

Maybe I'm misunderstanding? quickmodel already is a command-line tool that uses argparse? https://github.com/choderalab/assaytools/blob/master/scripts/quickmodel.py#L424

[sonyahanson]

I have been using those values (nburn=500, niter=1000 to nthin=50) for a while, and the output was not garbage, it is difficult for me to compare your other changes to my current results if these things are changed unnecessarily. The idea is to see if these changes improve our results...

[jchodera]

I can certainly change the defaults back, but I honestly wasn't able to verify that things were working with the old parameters on any of the datasets because the output was so problematic.

Here's a concrete example on bostutinib. Compare the old parameters

nthin=50, nburn=500, niter=1000

with the new parameters

nthin=50, nburn=0, niter=20000

[jchodera]

@sonyahanson : Mystery solved! The defaults were changed in this PR, and this change specifically.

@gregoryross had moved those settings into optional arguments, and probably reduced them for testing.

I'll have quickmodel use the recommended defaults.

[jchodera]

• Default quickmodel parameters (prior to being accidentally changed in #79) have been restored
• quickmodel now supports --nsamples and --nthin command-line arguments
• Travis now uses --nsamples 20 to reduce test runtime

Anything else we should put in this PR?

[sonyahanson]

Hmmm, I think this is still not back to what it was. I think the way this original commit was done was good: https://github.com/choderalab/assaytools/commit/b662469912886bcda96b50876197d269bca2c444

[jchodera]

What specifically are you looking for in terms of the default settings here?

[jchodera]

(I can't tell if you're suggesting to restore the behavior before or after that particular PR)

[sonyahanson]

I think it's fine, I just wanted to reference that commit since I found it, and since it's what I thought was happening. I can change the quickmodel parameters to match the previous setting no problem.

[jchodera]

Let's keep nburn=0 because it really is unnecessary (and discards effort) given that you are already using automated equilibration detection.

[sonyahanson]

I hadn't realized that Greg had made a second commit that changed this further, the original behavior I have been using is in fact this:

   nthin = 20
   nburn = nthin*10000
   niter = nthin*10000

[jchodera]

I hadn't realized that Greg had made a second commit that changed this further, the original behavior I have been using is in fact this:

The current behavior does the same amount of sampling, but does not arbitrarily discard the first 10000 samples. Instead, it lets automatic equilibration decide how much to discard.

[jchodera]

Whoops, no, I'm wrong---I confused the number of samples with the number of iterations. Fixing!

[sonyahanson]

I'm still on the fence about the automated equilibration detection, sometimes it makes things worse, and I have yet to see a case where it makes things better, but perhaps it is because we actually have had quite an extensive nburn phase.

[jchodera]

I'm still on the fence about the automated equilibration detection, sometimes it makes things worse, and I have yet to see a case where it makes things better, but perhaps it is because we actually have had quite an extensive nburn phase.

Let's keep nburn = 0 and stick with automatic equilibration detection until we have good evidence that it is causing harm.

[jchodera]

Oh, but are you running equilibration detection on the log deviance or the DeltaG?

[sonyahanson]

DeltaG

[jchodera]

I suspect it should be used on the log deviance instead, but let's run with the current "restored" behavior for now and see how it goes.

[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 https://github.com/choderalab/assaytools/pull/83/commits/fcde1cfbef1bbf85d245a6a1706b280297e75c1a 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.