Merge pull request #319 from Lnaden/revert_returns

Make all outputs default to 3.0.3 behavior

README.md

@@ -45,14 +45,15 @@ To analyze this data, we first initialize the `MBAR` object:
 ```
 Estimating dimensionless free energy differences between the sampled thermodynamic states and their associated uncertainties (standard errors) simply requires a call to `getFreeEnergyDifferences()`:
 ```python
->>>  results = mbar.getFreeEnergyDifferences()
+>>>  results = mbar.getFreeEnergyDifferences(return_dict=True)
 ```
 Here `results` is a dictionary with keys `Deltaf_ij`, `dDeltaf`, and `Theta`. `Deltaf_ij[i,j]` is the matrix of dimensionless free energy differences `f_j - f_i`, `dDeltaf_ij[i,j]` is the matrix of standard errors in this matrices estimate, and `Theta` is a covariance matrix that can be used to propagate error into quantities derived from the free energies.
+By default, `return_dict` is `False` and the return will be a tuple.
 
 Expectations and associated uncertainties can easily be estimated for observables `A(x)` for all states:
 ```python
 >>> A_kn = x_kn # use position of harmonic oscillator as observable
->>> results = mbar.computeExpectations(A_kn)
+>>> results = mbar.computeExpectations(A_kn, return_dict=True)
 ```
 where `results` is a dictionary with keys `mu`, `sigma`, and `Theta`, where `mu[i]` is the array of the estimate for the average of the observable for in state i, `sigma[i]` is the estimated standard deviation of the `mu` estimates,  and `Theta[i,j]` is the covarinace matrix of the log weights. 
 

examples/constant-force-optical-trap/force-bias-optical-trap.py

@@ -211,9 +211,6 @@ def construct_nonuniform_bins(x_n, nbins):
 # Compute PMF in unbiased potential (in units of kT).
 print("Computing PMF...")
 (f_i, df_i) = mbar.computePMF(u_kn, bin_kn, nbins)
-results = mbar.computePMF(u_kn, bin_kn, nbins)
-f_i = results['f_i']
-df_i = results['df_i']
 # compute estimate of PMF including Jacobian term
 pmf_i = f_i + numpy.log(bin_width_i)
 # Write out unbiased estimate of PMF
@@ -237,7 +234,7 @@ def construct_nonuniform_bins(x_n, nbins):
 # compute observed and expected histograms at each state
 for l in range(0,K):
     # compute PMF at state l
-    results = mbar.computePMF(u_kln[:,l,:], bin_kn, nbins)
+    results = mbar.computePMF(u_kln[:,l,:], bin_kn, nbins, return_dict=True)
     f_i = results['f_i']
     df_i = results['df_i']
     # compute estimate of PMF including Jacobian term

examples/harmonic-oscillators/harmonic-oscillators-distributions.py

@@ -135,7 +135,7 @@
 
   # Initialize the MBAR class, determining the free energies.
   mbar = MBAR(u_kln, N_k, relative_tolerance=1.0e-10,verbose=False) # use fast Newton-Raphson solver
-  results = mbar.getFreeEnergyDifferences()
+  results = mbar.getFreeEnergyDifferences(return_dict=True)
   Deltaf_ij_estimated = results['Delta_f']
   dDeltaf_ij_estimated = results['dDelta_f']
 

examples/harmonic-oscillators/harmonic-oscillators.py

@@ -25,7 +25,8 @@
 from __future__ import print_function
 import sys
 import numpy
-from pymbar import testsystems, EXP, EXPGauss, BAR, MBAR, utils
+from pymbar import testsystems, EXP, EXPGauss, BAR, MBAR
+from pymbar.utils import ParameterError
 
 #=============================================================================================
 # HELPER FUNCTIONS
@@ -142,7 +143,7 @@ def GetAnalytical(beta,K,O,observables):
 print("      Testing getFreeEnergyDifferences       ")
 print("=============================================")
 
-results = mbar.getFreeEnergyDifferences()
+results = mbar.getFreeEnergyDifferences(return_dict=True)
 Delta_f_ij_estimated = results['Delta_f']
 dDelta_f_ij_estimated = results['dDelta_f']
 
@@ -174,7 +175,7 @@ def GetAnalytical(beta,K,O,observables):
   k1 = nonzero_indices[i+1]
   w_F = u_kln[k, k1, 0:N_k[k]]   - u_kln[k, k, 0:N_k[k]]       # forward work                                  
   w_R = u_kln[k1, k, 0:N_k[k1]] - u_kln[k1, k1, 0:N_k[k1]]    # reverse work                                  
-  results = BAR(w_F, w_R)
+  results = BAR(w_F, w_R, return_dict=True)
   df_bar = results['Delta_f']
   ddf_bar = results['dDelta_f']
   bar_analytical = (f_k_analytical[k1]-f_k_analytical[k]) 
@@ -190,7 +191,7 @@ def GetAnalytical(beta,K,O,observables):
 for k in range(K-1):
   if N_k[k] != 0:
     w_F = u_kln[k, k+1, 0:N_k[k]]   - u_kln[k, k, 0:N_k[k]]       # forward work
-    results = EXP(w_F)
+    results = EXP(w_F, return_dict=True)
     df_exp = results['Delta_f']
     ddf_exp = results['dDelta_f']
     exp_analytical = (f_k_analytical[k+1]-f_k_analytical[k]) 
@@ -202,7 +203,7 @@ def GetAnalytical(beta,K,O,observables):
 for k in range(1,K):
   if N_k[k] != 0:
     w_R = u_kln[k, k-1, 0:N_k[k]] - u_kln[k, k, 0:N_k[k]]         # reverse work                                  
-    (df_exp,ddf_exp) = EXP(w_R)
+    (df_exp,ddf_exp) = EXP(w_R, return_dict=True)
     df_exp = -results['Delta_f']
     ddf_exp = results['dDelta_f']
     exp_analytical = (f_k_analytical[k]-f_k_analytical[k-1]) 
@@ -218,7 +219,7 @@ def GetAnalytical(beta,K,O,observables):
 for k in range(K-1):
   if N_k[k] != 0:
     w_F = u_kln[k, k+1, 0:N_k[k]]   - u_kln[k, k, 0:N_k[k]]       # forward work                                  
-    results = EXPGauss(w_F)
+    results = EXPGauss(w_F, return_dict=True)
     df_gauss = results['Delta_f']
     ddf_gauss = results['dDelta_f']
     gauss_analytical = (f_k_analytical[k+1]-f_k_analytical[k]) 
@@ -230,7 +231,7 @@ def GetAnalytical(beta,K,O,observables):
 for k in range(1,K):
   if N_k[k] != 0:
     w_R = u_kln[k, k-1, 0:N_k[k]] - u_kln[k, k, 0:N_k[k]]         # reverse work                                  
-    results = EXPGauss(w_R)
+    results = EXPGauss(w_R, return_dict=True)
     df_gauss = results['Delta_f']
     ddf_gauss = results['dDelta_f']
     gauss_analytical = (f_k_analytical[k]-f_k_analytical[k-1]) 
@@ -286,7 +287,7 @@ def GetAnalytical(beta,K,O,observables):
     for k in range(0,K):
       A_kn[k,0:N_k[k]] = x_kn[k,0:N_k[k]]**2
 
-  results = mbar.computeExpectations(A_kn, state_dependent = state_dependent)
+  results = mbar.computeExpectations(A_kn, state_dependent = state_dependent, return_dict=True)
   A_k_estimated = results['mu']
   dA_k_estimated = results['sigma']
 
@@ -340,7 +341,7 @@ def GetAnalytical(beta,K,O,observables):
   stdevs = numpy.abs(As_k_error[Nk_ne_zero]/dAs_k_estimated[Nk_ne_zero])
   print(stdevs)
 
-  results = mbar.computeExpectations(A_kn, state_dependent = state_dependent, output = 'differences')
+  results = mbar.computeExpectations(A_kn, state_dependent = state_dependent, output = 'differences', return_dict=True)
   A_kl_estimated = results['mu']  
   dA_kl_estimated = results['sigma']
 
@@ -383,7 +384,7 @@ def GetAnalytical(beta,K,O,observables):
 for i,observe in enumerate(observables_single):
   A_ikn[i,:,:] = A_kn_all[observe]
 for i in range(K):
-  results = mbar.computeMultipleExpectations(A_ikn, u_kln[:,i,:], compute_covariance=True)
+  results = mbar.computeMultipleExpectations(A_ikn, u_kln[:,i,:], compute_covariance=True, return_dict=True)
   A_i = results['mu']
   dA_ij = results['sigma']
   Ca_ij = results['covariances']
@@ -398,7 +399,7 @@ def GetAnalytical(beta,K,O,observables):
 print("      Testing computeEntropyAndEnthalpy")
 print("============================================")
 
-results = mbar.computeEntropyAndEnthalpy(u_kn = u_kln, verbose = True)
+results = mbar.computeEntropyAndEnthalpy(u_kn = u_kln, verbose = True, return_dict=True)
 Delta_f_ij = results['Delta_f']
 dDelta_f_ij = results['dDelta_f']
 Delta_u_ij = results['Delta_u']
@@ -467,7 +468,7 @@ def GetAnalytical(beta,K,O,observables):
     for l in range(L):
       unew_kln[k,l,0:N_k[k]] = (K_extra[l]/2.0) * (x_kn[k,0:N_k[k]]-O_extra[l])**2
 
-results = mbar.computePerturbedFreeEnergies(unew_kln)
+results = mbar.computePerturbedFreeEnergies(unew_kln, return_dict=True)
 Delta_f_ij_estimated = results['Delta_f']
 dDelta_f_ij_estimated = results['dDelta_f']
 
@@ -519,7 +520,7 @@ def GetAnalytical(beta,K,O,observables):
     A_kn = A_kn_all['position^2']
 
   A_k_estimated, dA_k_estimated 
-  results = mbar.computeExpectations(A_kn,unew_kln[:,[nth],:],state_dependent=state_dependent)
+  results = mbar.computeExpectations(A_kn,unew_kln[:,[nth],:],state_dependent=state_dependent, return_dict=True)
   A_k_estimated = results['mu'] 
   dA_k_estimated = results['sigma']
   # need to additionally transform to get the square root
@@ -543,7 +544,7 @@ def GetAnalytical(beta,K,O,observables):
 print("      Testing computeOverlap   ")
 print("============================================")
 
-results = mbar.computeOverlap()
+results = mbar.computeOverlap(return_dict=True)
 O = results['scalar']
 O_i = results['eigenvalues']
 O_ij = results['matrix']
@@ -577,7 +578,7 @@ def GetAnalytical(beta,K,O,observables):
 print("so standard (scaled) results should be very close to MBAR results.")
 print("No standard estimate exists for states that are not sampled.")
 A_kn = x_kn
-results = mbar.computeExpectations(A_kn)
+results = mbar.computeExpectations(A_kn, return_dict=True)
 val_mbar = results['mu']
 err_mbar = results['sigma']
 err_standard = numpy.zeros([K],dtype = numpy.float64)
@@ -744,7 +745,7 @@ def generate_pmf_data(ndim=1, nbinsperdim=15, nsamples = 1000, K0=20.0, Ku = 100
 
 # Compute PMF.
 print("Computing PMF ...")
-results = mbar.computePMF(u_n, bin_n, nbins, uncertainties = 'from-specified', pmf_reference = zeroindex)
+results = mbar.computePMF(u_n, bin_n, nbins, uncertainties = 'from-specified', pmf_reference = zeroindex, return_dict=True)
 f_i = results['f_i']
 df_i = results['df_i']
 
@@ -833,7 +834,7 @@ def generate_pmf_data(ndim=1, nbinsperdim=15, nsamples = 1000, K0=20.0, Ku = 100
 # Compute PMF.          
 print("Computing PMF ...")
 [f_i, df_i] 
-results = mbar.computePMF(u_n, bin_n, nbins, uncertainties = 'from-specified', pmf_reference = zeroindex)
+results = mbar.computePMF(u_n, bin_n, nbins, uncertainties = 'from-specified', pmf_reference = zeroindex, return_dict=True)
 f_i = results['f_i']
 df_i = results['df_i']
 

examples/heat-capacity/heat-capacity.py

@@ -281,23 +281,23 @@ def PrintResults(string,E,dE,Cv,dCv,types):
     E_kln = u_kln  # not a copy, we are going to write over it, but we don't need it any more.
     for k in range(K):
         E_kln[:,k,:]*=beta_k[k]**(-1)  # get the 'unreduced' potential -- we can't take differences of reduced potentials because the beta is different; math is much more confusing with derivatives of the reduced potentials.
-    results = mbar.computeExpectations(E_kln, state_dependent = True)
+    results = mbar.computeExpectations(E_kln, state_dependent = True, return_dict=True)
     E_expect = results['mu']
     dE_expect = results['sigma']
     allE_expect[:,n] = E_expect[:]
 
     # expectations for the differences, which we need for numerical derivatives
-    results = mbar.computeExpectations(E_kln,output='differences', state_dependent = True)
+    results = mbar.computeExpectations(E_kln,output='differences', state_dependent = True, return_dict=True)
     DeltaE_expect = results['mu']
     dDeltaE_expect = results['sigma']
     print("Computing Expectations for E^2...")
 
-    results = mbar.computeExpectations(E_kln**2, state_dependent = True)
+    results = mbar.computeExpectations(E_kln**2, state_dependent = True, return_dict=True)
     E2_expect = results['mu']
     dE2_expect = results['sigma']
     allE2_expect[:,n] = E2_expect[:]
 
-    results = mbar.getFreeEnergyDifferences()
+    results = mbar.getFreeEnergyDifferences(return_dict=True)
     df_ij = results['Delta_f']
     ddf_ij = results['dDelta_f']
 

examples/parallel-tempering-2dpmf/parallel-tempering-2dpmf.py

@@ -293,7 +293,7 @@ def read_file(filename):
 # f_i[i] is the dimensionless free energy of bin i (in kT) at the temperature of interest
 # df_i[i,j] is an estimate of the covariance in the estimate of (f_i[i] - f_j[j], with reference
 # the lowest free energy state.
-results = mbar.computePMF(u_kn, bin_kn, nbins, uncertainties='from-lowest')
+results = mbar.computePMF(u_kn, bin_kn, nbins, uncertainties='from-lowest', return_dict=True)
 f_i = results['f_i']
 df_i = results['df_i']
 

examples/umbrella-sampling-pmf/umbrella-sampling.py

@@ -147,7 +147,7 @@
 mbar = pymbar.MBAR(u_kln, N_k, verbose = True)
 
 # Compute PMF in unbiased potential (in units of kT).
-results = mbar.computePMF(u_kn, bin_kn, nbins)
+results = mbar.computePMF(u_kn, bin_kn, nbins, return_dict=True)
 f_i = results['f_i']
 df_i = results['df_i']
 

pymbar/bar.py

@@ -147,7 +147,7 @@ def BARzero(w_F, w_R, DeltaF):
     return fzero
 
 
-def BAR(w_F, w_R, DeltaF=0.0, compute_uncertainty=True, uncertainty_method='BAR',maximum_iterations=500, relative_tolerance=1.0e-12, verbose=False, method='false-position', iterated_solution=True):
+def BAR(w_F, w_R, DeltaF=0.0, compute_uncertainty=True, uncertainty_method='BAR',maximum_iterations=500, relative_tolerance=1.0e-12, verbose=False, method='false-position', iterated_solution=True, return_dict=False):
     """Compute free energy difference using the Bennett acceptance ratio (BAR) method.
 
     Parameters
@@ -176,16 +176,18 @@ def BAR(w_F, w_R, DeltaF=0.0, compute_uncertainty=True, uncertainty_method='BAR'
     iterated_solution : bool, optional, default=True
         whether to fully solve the optimized BAR equation to consistency, or to stop after one step, to be 
         equivalent to transition matrix sampling.
+    return_dict : bool, default False
+        If true, returns are a dict, else they are a tuple
 
     Returns
     -------
-    result_vals : dictionary
-    
-    Possible keys in the result_vals dictionary
     'Delta_f' : float
         Free energy difference
+        If return_dict, key is 'Delta_f'
     'dDelta_f': float
         Estimated standard deviation of free energy difference
+        If return_dict, key is 'dDelta_f'
+
 
     References
     ----------
@@ -203,15 +205,15 @@ def BAR(w_F, w_R, DeltaF=0.0, compute_uncertainty=True, uncertainty_method='BAR'
 
     >>> from pymbar import testsystems
     >>> [w_F, w_R] = testsystems.gaussian_work_example(mu_F=None, DeltaF=1.0, seed=0)
-    >>> results = BAR(w_F, w_R)
+    >>> results = BAR(w_F, w_R, return_dict=True)
     >>> print('Free energy difference is {:.3f} +- {:.3f} kT'.format(results['Delta_f'], results['dDelta_f']))
     Free energy difference is 1.088 +- 0.050 kT
 
     Test completion of various other schemes.
 
-    >>> results = BAR(w_F, w_R, method='self-consistent-iteration')
-    >>> results = BAR(w_F, w_R, method='false-position')
-    >>> results = BAR(w_F, w_R, method='bisection')
+    >>> results = BAR(w_F, w_R, method='self-consistent-iteration', return_dict=True)
+    >>> results = BAR(w_F, w_R, method='false-position', return_dict=True)
+    >>> results = BAR(w_F, w_R, method='bisection', return_dict=True)
 
     """
 
@@ -228,8 +230,8 @@ def BAR(w_F, w_R, DeltaF=0.0, compute_uncertainty=True, uncertainty_method='BAR'
         nfunc = 0
 
     if method == 'bisection' or method == 'false-position':
-        UpperB = EXP(w_F)['Delta_f']
-        LowerB = -EXP(w_R)['Delta_f']
+        UpperB = EXP(w_F, return_dict=True)['Delta_f']
+        LowerB = -EXP(w_R, return_dict=True)['Delta_f']
 
         FUpperB = BARzero(w_F, w_R, UpperB)
         FLowerB = BARzero(w_F, w_R, LowerB)
@@ -242,10 +244,14 @@ def BAR(w_F, w_R, DeltaF=0.0, compute_uncertainty=True, uncertainty_method='BAR'
             if compute_uncertainty:
                 result_vals['Delta_f'] = 0.0 
                 result_vals['dDelta_f'] = 0.0
-                return result_vals
+                if return_dict:
+                    return result_vals
+                return 0.0, 0.0
             else:
-                result_vals['Delta_f'] = 0.0 
-                return result_vals
+                result_vals['Delta_f'] = 0.0
+                if return_dict:
+                    return result_vals
+                return 0.0
 
         while FUpperB * FLowerB > 0:
             # if they have the same sign, they do not bracket.  Widen the bracket until they have opposite signs.
@@ -487,12 +493,16 @@ def BAR(w_F, w_R, DeltaF=0.0, compute_uncertainty=True, uncertainty_method='BAR'
             print("DeltaF = {:8.3f} +- {:8.3f}".format(DeltaF, dDeltaF))
         result_vals['Delta_f'] = DeltaF
         result_vals['dDelta_f'] = dDeltaF
+        if return_dict:
+            return result_vals
+        return DeltaF, dDeltaF
     else:
         if verbose:
             print("DeltaF = {:8.3f}".format(DeltaF))
         result_vals['Delta_f'] = DeltaF
-
-    return result_vals
+        if return_dict:
+            return result_vals
+        return DeltaF
 
 #=============================================================================================
 # For compatibility with 2.0.1-beta