choderalab/pymbar — pull request #102 by kyleabeauchamp (2014-08-09 21:59:53+00:00)
Title: Fast numexpr / scipy MBAR solver.

Lots of stuff going on here, but it's all logically related so a "big" PR is mostly unavoidable.

1. Moved current version of MBAR to `old_mbar.py` as a reference implementation for testing against.
2. Cleared all MBAR solver code into `mbar_solvers.py`
3. Added fast NumExpr based `logsumexp()` function to use as foundation of fast-but-precise MBAR calculations. If NumExpr is unavailable, a pure-numpy implementation is used.
4. Implemented matrix multiplication-based faster-but-imprecise MBAR calculations.
5. Replaced hand-written solvers with a single interface to `scipy.optimize.minimize` and `scipy.optimize.root`, which are scipy's unified interfaces to minimization and root-finding algorithms, respectively.
6. Because balancing precision and performance is pretty challenging, we need to make multiple passes through the solver--first get a rough guess with BFGS, then refine it with the MINPACK implementation of NR / HYBR. `mbar_solvers.solve_mbar()` uses a "protocol" to call BFGS and NR code several times in a way that achieves similar precision while matching the performance of the old MBAR code.
7. Added `pymbar/testsystems/pymbar_datasets.py` to automatically load the large datasets from compressed HDF5 files, which I will add to the `pymbar-datasets` repository soon.
8. Added `tests/test_mbar_solvers.py` to test the accuracy of the new MBAR code on several challenging systems, and also compare to the old MBAR code.
9. Added scripts/benchmark_pymbar.py to do performance and accuracy benchmarks against the old mbar code on several challenging datasets, including the `gas` and `8proteins` datasets.

Overall, the speed is somewhat improved for the challenging `gas` and `8proteins` datasets. However, this restructuring helps us eliminate the `weave` dependence (needed for Py3 compliance) and the hand-written C code. I also think that `mbar_solvers.py` is an improvement in readability. If we need more speed in the future, we can plug in faster versions of the gradient / hessian calculation.