choderalab/pymbar — pull request #83 by mrshirts (2014-07-09 23:18:52+00:00)
Title: proposed conversion from u_kln->u_kn formalism

includes a slightly modified harmonic_oscillator_update.py, which gives identical behavior with the new code to pymbar-examples/harmonic_oscillator.py run with the previous code up to numerical precision.

This new code automatically converts from old-style u_kln matrices internally.

It supports all existing functionality, including updating the weave code and the _pymbar.c code.

The one major difference is in handling computeExpectations. This change hopefully provides a much more general solution which can possibly reduce the total amount of repetitive code as well.

Previously, computeExpectations recognized whether observables were state_dependent or state_independent by whether the observable array was KxN_max, or KxLxN_max. Now, this is done by giving an explicit state_dependent flag. If state_dependent is False, computeExpectations handles the code like the 2D version before (although now the array is 1D, the observable at each state).
The big difference is if state_dependent is True. In this case, computeExpectations expects K different observables (the K different observables at each state - for example, the potential energies in each of the K states), as well as K different energies (which defaults to just the u_kn matrix MBAR was generated with). It then calculates the expectations of the first observable and the first state, 2nd observable and the 2nd state, and so forth, which is what the previous code did.

In order to make this happen, I've created a computeGeneralExpectations routine, that takes in I observable arrays and K state observables. In theory, this would mean that there are I*K observables, which means that there are (I*K)^2 terms in the covariance matrix, which is just silly. Instead, this computeGeneralExpectation matrix expects a Sx2 array state_list which describes how many average observables should be computed.

I have demonstrated this by allowing computeExpectations to perform it's functions by calling computeGeneralExpectation for all the heavy lifting. For the state_dependent=False case, the state_list would be [[0,1],[0,2],[0,3]. . . [0,K]], which means that the 0th observable is evaluated at the K different state, for S=K total states. For the state_dependent=True case, the state_list would look like [[0,0],[1,1],[2,2]. . . [K,K]], so there are S=K average observables, with the 0th observable averaged at the 0th state, 1st observable averaged at the 1st state, etc.

Note that this computeGeneralExpectations can easily be used to do the heavy lifting with computeMultipleExpectations, with state list [1,0],[2,0],[3,0]. . . [K,0]], but I haven't done that until people look it everything. It should work just fine for computePerturbedExpectations as well, and possibly computeEntropyandEnthalpy with some tweaking. I'd propose using some version of computeGeneralExpectations to replace a number of the existing repetitive codes eventually.

computeGeneralExpectations can be improved (right now, the u_ln that you add are recomputed even if they are already solved for in the u_kn matrix), but should be a good baseline.

I'm up for lots of other improvements; but I though it would be easier to start by solving the biggest problem first independently of other issues.