#!/usr/bin/python

# Example illustrating the use of MBAR for computing the hydration free energy of OPLS 3-methylindole
# in TIP3P water through alchemical free energy simulations.

#===================================================================================================
# IMPORTS
#===================================================================================================

import numpy
import pymbar # multistate Bennett acceptance ratio estimator.
import timeseries # for timeseries analysis 
import commands
import re
import pdb

#===================================================================================================
# CONSTANTS
#===================================================================================================

convert_atmnm3_to_kJmol = 1.01325e5*(1e-09)**3 * 6.02214 * (1e23) / 1000 # Convert pV from atmospheres*nm^3 into kJ/mol
kB = 1.381*6.02214/1000.0  # Boltzmann's constant (kJ/mol/K)
relative_tolerance = 1e-10
verbose = False
optional = False  # set this to false if we only want the results for the paper
NR = 1
Nboot = 10
plotspline = False

#methods = ['TI','TI-CUBIC','FEXP','REXP','BAR','MBAR'];
#methods = ['TI','TI-CUBIC','FEXP','REXP','BAR'];
#methods = ['TI','TI-CUBIC'];
methods = ['FEXP','REXP','BAR','UBAR','MBAR','UMBAR'];

############vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv#######################
#####THIS SECTION OF METHODS CAN BE TAKEN OUT OF THE DISTRIBUTION########

# expanded list of methods
#methods = ['TI','FEXP','REXP','UBAR','RBAR','BAR','PMBAR','MBAR']; 
#methods = ['TI','FEXP','REXP','UBAR','RBAR','BAR','PMBAR','MBAR','UMBAR']; 

#####THIS SECTION OF METHODS CAN BE TAKEN OUT OF THE DISTRIBUTION########
############^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^#######################


#===================================================================================================
# PARAMETERS
#===================================================================================================

temperature = 298.0 # temperature (K)
pressure = 1.0 # pressure (atm)
nequil = 100; # number of samples assumed to be equilibration
#datafile_directory = '../../large-datasets/trp50ns/' # directory in which datafiles are store
#datafile_directory = '../../large-datasets/trp38/' # directory in which datafiles are stored
#datafile_directory = '../../trunk/examples/alchemical-free-energy/data/3-methylindole-38steps/' # directory in which datafiles are stored
#datafile_directory = '../../trunk/examples/alchemical-free-energy/data/3-methylindole-11steps/' # directory in which datafiles are stored
datafile_directory = '../../manuscripts/comment-on-entropy-estimator/dgdata/'
datafile_prefix  = 'UAM_short_mgibbs.dhdl' # prefixes for datafile sets 

#===================================================================================================
# HELPER FUNCTIONS
#===================================================================================================
def sortbynum(item):
   vals = item.split('.')
   for v in reversed(vals):
      if v.isdigit():
         return int(v)
   print "Error: No digits found in filename, can't sort ", item

class naturalcubicspline:

   def __init__(self, x):
      
      L = len(x)
      
      self.x = x.copy()

      # define some space
      H = numpy.zeros([L,L],float)
      M = numpy.zeros([L,L],float)
      BW = numpy.zeros([L,L],float)
      AW = numpy.zeros([L,L],float)
      DW = numpy.zeros([L,L],float)

      self.h = self.x[1:L]-x[0:L-1]
      ih = 1.0/self.h
      h = self.h

      # define the H and M matrix, from p. 371 "applied numerical methods with matlab"
      H[0,0] = 1
      H[L-1,L-1] = 1
      for i in range(1,L-1):
         H[i,i] = 2*(h[i-1]+h[i])
         H[i,i-1] = h[i-1]
         H[i,i+1] = h[i]

      for i in range(1,L-1):
         M[i,i] = -3*(ih[i-1]+ih[i])
         M[i,i-1] = 3*(ih[i-1])
         M[i,i+1] = 3*(ih[i])
            
      CW = numpy.dot(numpy.linalg.inv(H),M)  # this is the matrix translating c to weights in f.
                                                   # each row corresponds to the weights for each c.
            
      # from CW, define the other coefficient matrices
      for i in range(0,L-1):
         BW[i,:]    = -(h[i]/3)*(2*CW[i,:]+CW[i+1,:])
         BW[i,i]   += -ih[i]   
         BW[i,i+1] += ih[i]
         DW[i,:]    = (ih[i]/3)*(CW[i+1,:]-CW[i,:])
         AW[i,i]    = 1

      self.AW = AW.copy()
      self.BW = BW.copy()
      self.CW = CW.copy()
      self.DW = DW.copy()
     
      # find the integrating weights

      self.wsum = numpy.zeros([L],float)
      self.wk = numpy.zeros([L-1,L],float)
      dlamw = numpy.zeros([4],float)
      for k in range(0,L-1):
         dlamw[0] = (self.h[k])
         dlamw[1] = (self.h[k]**2)/2.0
         dlamw[2] = (self.h[k]**3)/3.0
         dlamw[3] = (self.h[k]**4)/4.0
         w = DW[k,:]*dlamw[3] + CW[k,:]*dlamw[2] + BW[k,:]*dlamw[1] + AW[k,:]*dlamw[0]
         self.wk[k,: ] = w 
         self.wsum += w

   def interpolate(self,y,xnew):

      L = len(self.x)
      if len(y) != L:
         print "Error"
          
      # get the array of actual coefficients by multiplying the coefficient matrix by the values  
      a = numpy.dot(self.AW,y)
      b = numpy.dot(self.BW,y)
      c = numpy.dot(self.CW,y)
      d = numpy.dot(self.DW,y)

      N = len(xnew)
      ynew = numpy.zeros([N],float)
      delta = min(xnew[1:N]-xnew[0:N-1])/2
      ind = numpy.searchsorted(self.x, xnew+delta,side='left')-1
      h = xnew-self.x[ind]
      ynew = d[ind]*h*h*h + c[ind]*h*h + b[ind]*h + a[ind];

      return ynew
      
#===================================================================================================
# MAIN
#===================================================================================================

# Compute inverse temperature in 1/(kJ/mol)
beta = 1./(kB*temperature)

#===================================================================================================
# Read all snapshot data
#===================================================================================================

# Generate a list of all datafiles with this prefix.
# NOTE: This uses the file ordering provided by 'ls' to determine which order the states are in.
filenames = commands.getoutput('ls %(datafile_directory)s/%(datafile_prefix)s.*' % vars()).split()

# Determine number of alchemical intermediates.
K = len(filenames)
# Determine length of files.
filenames = commands.getoutput('ls %(datafile_directory)s/%(datafile_prefix)s*' % vars()).split()   

# sort the files numerically.
filenames.sort(key=sortbynum)
   
nsnapshots = numpy.zeros(K, int) # nsnapshots[k] is the number of snapshots for state k; 
# no larger than number of lines starting out.
for k in range(K):
   # Temporari1ly read the file into memory.
   infile = open(filenames[k], 'r')
   lines = infile.readlines()
   infile.close()
   
   # Determine maxnumber of snapshots from quickly parsing file and ignoring header lines.
   nsnapshots[k] = 0
   for line in lines:
      if ((line[0] == '#') or (line[0] == '@')):
         continue
      nsnapshots[k] += 1
maxn = max(nsnapshots)   

# Load all of the data
u_klt = numpy.zeros([K,K,maxn], numpy.float64) # u_klt[k,m,t] is the reduced potential energy of snapshot t of state k evaluated at state m
for k in range(K):
   # File to be read
   filename = filenames[k]   
   
      # Read contents of file into memory.
   print "Reading %s..." % filename
   infile = open(filename, 'r')
   lines = infile.readlines()
   infile.close()

   t = 0
   n_components = 0
   n_states = 0
   bPV = False

   # Parse the file.
   pe = numpy.zeros(K,numpy.float64); # temporary storage for energies
   for line in lines:

      # split into elements
      elements = line.split()
      if ((line[0] == '#') or (line[0] == '@')):
         if (line[0] == '@'):
            # it's an xvg legend entry -- load in the information
            if (line[2] == 's'):  
               # it's a xvg entry, and a lambda component, so note it down
               if (re.search('-lambda',line)):     
                  #it's a listing of the lambdas
                  n_components +=1
                  lv = numpy.zeros([K,n_components],float)
               elif (re.search("\\\\8D\\\\4H \\\\8\\l\\\\4",line)): 
                  lambda_string = elements[5]
                  lambda_list = re.sub('[()\"]','',lambda_string)
                  lambdas = lambda_list.split(',');
                  for i in range(n_components):
                     lv[n_states,i] = lambdas[i]
                  n_states+=1;   
               #elif (re.search("pv",line)):     
               #    bPV = 1;   # for testing now, eliminated PV term
      else:                           
         if ((t==0) and (k==0)):     # we don't know the number of components until here.
            dhdlt = numpy.zeros([K,n_components,maxn],float) 

         time = float(elements.pop(0))

         #
         #  If we print the energy in the dhdl file; if not, delete this line.
         #
         energy = float(elements.pop(0))            
            
         # now record the derivative with respect to lambda
      
         for nl in range(n_components):
            dhdlt[k,nl,t] = float(elements.pop(0))
         # now record the potential energy differences.   
         for l in range(K):
            pe[l] = float(elements.pop(0)) + energy
                  
         # pressure-volume contribution
         if (bPV):
            pv = float(elements.pop(0))
         else:
            pv = 0
         # compute and store reduced potential energy at each state
         for l in range(K):
            u_klt[k,l,t] = beta * (pe[l] + pv)
               
         t += 1   
   nsnapshots[k] = t

#===================================================================================================
# Preliminaries: Subsample data to obtain uncorrelated samples
#===================================================================================================   
Nequilibrated = max(nsnapshots) - nequil
ur_kln = numpy.zeros([K,K,Nequilibrated], numpy.float64) # ur_kln[k,m,n] is the reduced potential energy of uncorrelated sample index n from state k evaluated at state m
dhdlr = numpy.zeros([K,n_components,Nequilibrated], float) #dhdl is value for dhdl for each component in the file at each time.
NR_k = numpy.zeros(K, int) # NR_k[k] is the number of uncorrelated samples from state k over all replicates
g = numpy.zeros(K,float) # autocorrelation times for the data
for k in range(K):
   # Determine indices of uncorrelated samples from potential autocorrelation analysis at state k.
   # alternatively, could use the energy differences -- here, we will use total dhdl
   dhdl_sum = numpy.sum(dhdlt[k,:,:],axis=0)
   g[k] = timeseries.statisticalInefficiency(dhdl_sum[nequil:nsnapshots[k]])
   indices = numpy.array(timeseries.subsampleCorrelatedData(dhdl_sum[nequil:nsnapshots[k]],g=g[k])) # indices of uncorrelated samples
   N = len(indices) # number of uncorrelated samples
   NR_k[k] = N      
   indices += nequil
   for n in range(n_components):
      dhdlr[k,n,0:N] = dhdlt[k,n,indices]
   for l in range(K):         
      ur_kln[k,l,0:N] = u_klt[k,l,indices]

print "Correlation times:"
print g
print ""
print "number of uncorrelated samples:"
print NR_k
print ""

#Now, figure out how many replicates we have:
SN_k = numpy.zeros([NR+1,K],int)
for k in range(K):
   for r in range(NR):
      SN_k[r+1,k] += SN_k[r,k] + NR_k[k]/NR
   # last one might have a few more based on divisibilities
   SN_k[NR,k] = NR_k[k]
 
#===================================================================================================
# Preliminaries: Calculate average TI values
#===================================================================================================   


replicate_dF = list()
replicate_ddF = list()
replicate_ddFboot = list()

for r in range(NR):
   print "Replicate: ",
   print r+1

   # for bootstrap: first is the normal arrangement, then the rest are random indices.

   N_k = SN_k[r+1,:] - SN_k[r,:]

   dhdl_original = numpy.zeros([K,n_components,maxn], float)
   u_kln_original = numpy.zeros([K,K,maxn],numpy.float64)

   for k in range(K):
      u_kln_original[k,:,0:N_k[k]] = ur_kln[k,:,SN_k[r,k]:SN_k[r+1,k]] 
      dhdl_original[k,:,0:N_k[k]] = dhdlr[k,:,SN_k[r,k]:SN_k[r+1,k]]        

   u_kln = numpy.zeros([K,K,max(N_k)],numpy.float64)
   dhdl = numpy.zeros([K,n_components,max(N_k)],float)
   random_indices = numpy.zeros([K,Nboot+1,maxn],int);
   for k in range(K):
      random_indices[k,0,0:N_k[k]] = range(N_k[k]);
      
   for b in range(1,Nboot+1):
      numpy.random.seed()
      for k in range(K):
         random_indices[k,b,0:N_k[k]]=numpy.random.random_integers(0,high=N_k[k]-1,size=N_k[k])

   lchange = numpy.zeros([K,n_components],bool)   # booleans for which lambdas are changing 
   for j in range(n_components):    
      # need to identify range over which lambda doesn't change, and not interpolate over that range.
      for k in range(K-1):
         if (lv[k+1,j]-lv[k,j] > 0):
            lchange[k,j] = True 
            lchange[k+1,j] = True 

   # construct a map back to the original components, and forward

   mapk = numpy.zeros([K,n_components],int)   # map back to the original k from the components
   mapl = numpy.zeros([K,n_components],int)
   mapc = 0

   for j in range(n_components):    
      for k in range(sum(lchange[:,j])-1):
         mapk[k,j] = k + mapc
         mapl[k+mapc,j] = k
      mapc += (k+1)

   # put together the spline weights for the different components
   cubspl = list()  
   for j in range(n_components):    
      spl = naturalcubicspline(lv[lchange[:,j],j])
      cubspl.append(spl)

   bootstrap_list = list()

   for b in range(Nboot+1):
      print "Bootstrap: ",
      print b
      for k in range(K):
         for l in range(K):
            u_kln[k,l,0:N_k[k]] = u_kln_original[k,l,random_indices[k,b,0:N_k[k]]]
         for n in range(n_components):
            dhdl[k,n,0:N_k[k]] = dhdl_original[k,n,random_indices[k,b,0:N_k[k]]]

   #===================================================================================================
   # Calculate free energies with different methods
   #===================================================================================================    
            
      df_allk = list()
      ddf_allk = list()

      ave_dhdl = numpy.zeros([K,n_components],float)
      std_dhdl = numpy.zeros([K,n_components],float)

      # compute TI data

      for k in range(K):
         ave_dhdl[k,:] = beta*numpy.average(dhdl[k,:,0:N_k[k]],axis=1)        
         std_dhdl[k,:] = beta*numpy.std(dhdl[k,:,0:N_k[l]],axis=1)/numpy.sqrt(N_k[k]-1)

      # separate interpolation for each of the components.  Only interpolate range over which lambda changes.
      # lets look at how good the splines are:

         
      if (plotspline):   
         import scipy
         import scipy.interpolate
         import matplotlib.pylab as plt
         xin = numpy.arange(0,1,0.0005)
         yint1 = numpy.zeros([len(xin),n_components],float)
         yint2 = numpy.zeros([len(xin),n_components],float)
         for j in range(n_components):
            yint1[:,j] = cubspl[j].interpolate(ave_dhdl[lchange[:,j],j],xin)
            rep3 = scipy.interpolate.splrep(lv[lchange[:,j],j],ave_dhdl[lchange[:,j],j],s=0,k=3)
            yint2[:,j] = scipy.interpolate.splev(xin,rep3)
            plt.figure(j+1)
            plt.plot(xin,yint1[:,j],'r')
            plt.plot(xin,yint2[:,j],'b')
            plt.savefig('interpolate-' + str(j) + '-' + str(len(lv[:,j])) + '.pdf')

      for k in range(K-1):
         df = dict()
         ddf = dict()

         dlam = lv[k+1,:]-lv[k,:]

         for name in methods:
            if name == 'TI':
               #===================================================================================================
               # Estimate free energy difference with TI.
               #===================================================================================================   
               df['TI'] = numpy.dot(dlam/2,ave_dhdl[k]+ave_dhdl[k+1])        
               ddf['TI'] = numpy.sqrt(numpy.dot((dlam/2)**2,(std_dhdl[k])**2+(std_dhdl[k+1])**2)) 
               # above gives the correct error for 2 states, but consecutive parts are correlated, so we need to redo at the end.
                                      
            if name == 'TI-CUBIC':
               df['TI-CUBIC'] = 0 
               ddf['TI-CUBIC'] = 0 

               for j in range(n_components):
                  if dlam[j] > 0:
                     lj = lchange[:,j]
                     df['TI-CUBIC'] += numpy.dot(cubspl[j].wk[mapl[k,j]],ave_dhdl[lj,j])
                        #print k, df['TI-CUBIC'], df_j
                     ddf[name] += numpy.dot(cubspl[j].wk[mapl[k,j]]**2,std_dhdl[lj,j]**2)
               ddf['TI-CUBIC'] = numpy.sqrt(ddf['TI-CUBIC'])
      
            if name == 'FEXP':
               #===================================================================================================
               # Estimate free energy difference with Forward-direction EXP.
               #===================================================================================================   
                           
               w_F = u_kln[k,k+1,0:N_k[k]] - u_kln[k,k,0:N_k[k]] 
               (df['FEXP'], ddf['FEXP']) = pymbar.EXP(w_F)

            if name == 'REXP':
               #===================================================================================================
               # Estimate free energy difference with Reverse-direction EXP.
               #===================================================================================================   
               w_R = u_kln[k+1,k,0:N_k[k+1]] - u_kln[k+1,k+1,0:N_k[k+1]] 
               (rdf,rddf) = pymbar.EXP(w_R)
               (df['REXP'], ddf['REXP']) = (-rdf,rddf)

            if name == 'BAR':
               #===================================================================================================
               # Estimate free energy difference with BAR.
               #===================================================================================================   

               # w_F and w_R computed above (if order is ok . . . )                     
               (df['BAR'], ddf['BAR']) = pymbar.BAR(w_F, w_R, relative_tolerance=relative_tolerance, verbose = verbose)      

            if name == 'PMBAR':
               #===================================================================================================
               # Estimate free energy difference with Pairwise MBAR -- should be the same as BAR
               #===================================================================================================   
        
               MBAR = pymbar.MBAR(u_kln[k:(k+2),k:(k+2),:], N_k[k:(k+2)], relative_tolerance=relative_tolerance, verbose = verbose, method='Newton-Raphson')
               (Deltaf_ij, dDeltaf_ij) = MBAR.getFreeEnergyDifferences()	
               df['PMBAR'] = Deltaf_ij[0,1]
               ddf['PMBAR'] = dDeltaf_ij[0,1]     

            if name == 'UBAR':
                pdb.set_trace()
                #===================================================================================================
                # Estimate free energy difference with unoptimized BAR -- assume dF is zero, and just do one evaluation
                #===================================================================================================   
                # w_F and w_R computed above                      

                (df['UBAR'], ddf['UBAR']) = pymbar.BAR(w_F, w_R, optimized = no, verbose = verbose)      

            if name == 'RBAR':
                #===================================================================================================
                # Estimate free energy difference with Unoptimized BAR over range of free energy values, and choose the one 
                # That is self consistently best.
                #===================================================================================================   

                # w_F and w_R computed above                      
                min_diff = 1E6;
                best_udf = 0
                for trial_udf in range (-20,20,1):
                    (udf, uddf) = pymbar.BAR(w_F, w_R, DeltaF = trial_udf, optimized = no,verbose = verbose)      
                    diff = numpy.abs(udf - trial_udf);
                    if (diff < min_diff):
                       best_udf = udf
                       best_uddf = uddf
                       min_diff = diff
                (df['RBAR'], ddf['RBAR']) = (best_udf,best_uddf)

         df_allk = numpy.append(df_allk,df)
         ddf_allk = numpy.append(ddf_allk,ddf)

      for name in methods:
         if name == 'MBAR':

            #===================================================================================================
            # Estimate free energy difference with MBAR -- all states at once
            #===================================================================================================   

            # Initialize MBAR (computing free energy estimates, which may take a while)
            #print "Computing free energy differences..."
            MBAR = pymbar.MBAR(u_kln, N_k, verbose = verbose, method = 'Newton-Raphson', relative_tolerance = relative_tolerance) # use faster Newton-Raphson solver

            if (verbose):
               # Get matrix of dimensionless free energy differences and uncertainty estimate.
               print "Computing covariance matrix..."

            (MBAR_Deltaf_ij, MBAR_dDeltaf_ij) = MBAR.getFreeEnergyDifferences(uncertainty_method='svd-ew')

            if (verbose):
               # Matrix of free energy differences
               print "Deltaf_ij:"		
               print MBAR_Deltaf_ij
               # Matrix of uncertainties in free energy difference (expectations standard deviations of the estimator about the true free energy)
               print "dDeltaf_ij:"
               print MBAR_dDeltaf_ij
     	
            for k in range(K-1):
               df_allk[k]['MBAR'] = MBAR_Deltaf_ij[k,k+1]                      
               ddf_allk[k]['MBAR'] = MBAR_dDeltaf_ij[k,k+1]

         if name == 'UMBAR':
            #===================================================================================================
            # Estimate free energy difference with MBAR -- one step, matrix inversion method
            #===================================================================================================   
            # Initialize MBAR (computing free energy estimates, which may take a while)
            print "Computing free energy differences..."

            UMBAR = pymbar.MBAR(u_kln, N_k, verbose = verbose, method='Newton-Raphson',maximum_iterations = 1, relative_tolerance = relative_tolerance,initialize='zeros') 

            UMBAR_Deltaf_ij = UMBAR.getFreeEnergyDifferences(compute_uncertainty=False)
            UMBAR_dDeltaf_ij = numpy.zeros([K,K],float)
            cov_ij = numpy.zeros([K,K],float)
            # need to get variances another way
            # cov(fi-fj) = cov(fi,fi) + cov(fj,fj) - 2*cov(fi,fj)
            # cov(fi,fj) = cov(ln c_i \sum_k \sum_n_k W_n_k,i, ln c_j \sum_l \sum_n_l W_n_l,j) 
            #            = cov(ln c_i \sum_k \sum_n_k W_n_k,i, ln c_j \sum_l \sum_n_lW_n_l,j) 
            #      \approx cov(c_i \sum_k \sum_n_k W_n_k,i , c_j \sum_l \sum_n_l W_n_l,j) / (c_i \sum_k \sum_n_k W_n_k,i)(c_j \sum_l \sum_n_l W_n_l,j)  
            #            = cov(c_i \sum_k \sum_n_k W_n_k,i, c_j \sum_l \sum_n_l W_n_l,j) / (c_i)(c_j)  
            #            = cov(\sum_k \sum_n_k W_n_k,i, \sum_l \sum_n_l W_n_l,j)  
            #            = \sum_k \sum_l cov(\sum_n_k W_n_k,i, \sum_n_l W_n_l,j)  
            #            = \sum_k cov(\sum_n_k W_n_k,i, \sum_n_k W_n_k,j)  
            #            = \sum_k N_k**2 cov(W_n_k,i, W_n_k,j)  
            W_nk = UMBAR.W_nk
            for i in range(K):
               for j in range(K):
                  SumN = 0
                  for k in range(K):
                     m = numpy.vstack((W_nk[SumN:SumN+N_k[k],i],W_nk[SumN:SumN+N_k[k],j])) 
                     cov_ij[i,j] += (N_k[k]**2)*numpy.cov(m)[0,1]
                     SumN += N_k[k]
            for i in range(K):
               for j in range(K):
                  UMBAR_dDeltaf_ij[i,j] = numpy.sqrt(cov_ij[i,i]+cov_ij[j,j] - 2*cov_ij[i,j])

            for k in range(K-1):
                df_allk[k]['UMBAR'] = UMBAR_Deltaf_ij[k,k+1]                      
                ddf_allk[k]['UMBAR'] = UMBAR_dDeltaf_ij[k,k+1]


   #All done with calculations, now summarize and print stats
                

      if (b==0):           
         dF = dict()
         ddF = dict()
         ddFboot = dict()
         dU = dict()
         dS = dict()

         for name in methods:

            if name == 'MBAR': 
               dF['MBAR'] = MBAR_Deltaf_ij[0,K-1]                            
               ddF['MBAR'] = MBAR_dDeltaf_ij[0,K-1]
            elif name == 'UMBAR':
               dF['UMBAR'] = UMBAR_Deltaf_ij[0,K-1]                            
               ddF['UMBAR'] = UMBAR_dDeltaf_ij[0,K-1]
            elif name[0:2] == 'TI':
               dF[name] = 0   
               ddF[name] = 0
               for k in range(K-1):
                   dF[name] += df_allk[k][name]

               for j in range(n_components):      
                  lj = lchange[:,j]                     

                  if name == 'TI-CUBIC':
                     ddF[name] += numpy.dot((cubspl[j].wsum)**2,std_dhdl[lj,j]**2)
                        
                  if name == 'TI':
                     lam = lv[lj,j]
                     L = len(lam)
                     h = lam[1:L]-lam[0:L-1]
                     wsum = 0.5*(numpy.append(h,0) + numpy.append(0,h))
                     ddF[name] += numpy.dot(wsum**2,std_dhdl[lj,j]**2)

               ddF[name] = numpy.sqrt(ddF[name])

            else:
               dF[name] = 0
               ddF[name] = 0
               for k in range(K-1):
                  dF[name] += df_allk[k][name]
                  ddF[name] += (ddf_allk[k][name])**2
               ddF[name] = numpy.sqrt(ddF[name]);
      else:             
         bootlist = dict();
         for name in methods:
            if name == 'MBAR': 
               bootlist['MBAR'] = MBAR_Deltaf_ij[0,K-1]                            
            elif name == 'UMBAR':
               bootlist['UMBAR'] = UMBAR_Deltaf_ij[0,K-1]                            
            else:
               bootval = 0
               for k in range(K-1):
                  bootval += df_allk[k][name]
               bootlist[name] = bootval
         bootstrap_list.append(bootlist)
     
      if (b==Nboot):
         for name in methods:
            bvals = numpy.zeros([Nboot],float)
            for i in range(Nboot):
               bvals[i] = bootstrap_list[i][name]
            ddFboot[name] = numpy.std(bvals)

   replicate_dF = numpy.append(replicate_dF,dF)
   replicate_ddF = numpy.append(replicate_ddF,ddF)
   replicate_ddFboot = numpy.append(replicate_ddFboot,ddFboot)

for name in methods:
   print '%10s (kJ/mol):' % name,

   rep_dF = numpy.zeros(NR,float)   
   rep_ddF = numpy.zeros(NR,float);
   rep_ddFboot = numpy.zeros(NR,float);

   for r in range(NR):
      rep_dF[r] = replicate_dF[r][name]/beta
      rep_ddF[r] = replicate_ddF[r][name]/beta
      rep_ddFboot[r] = replicate_ddFboot[r][name]/beta

   ave_dF = numpy.average(rep_dF)
   std_dF = numpy.std(rep_dF,ddof=1)
   ave_ddF = numpy.average(rep_ddF)
   std_ddF = numpy.std(rep_ddF,ddof=1)

   boot_dF = numpy.average(rep_ddFboot)
   std_boot_dF = numpy.std(rep_ddFboot,ddof=1)

   print '%8.3f +/- %6.3f (predicted %6.3f +/- %6.3f, bootstrap %6.3f +/- %6.3f)' % (ave_dF,std_dF,ave_ddF,std_ddF,boot_dF,std_boot_dF),
   print rep_dF,
   print rep_ddF,
   print rep_ddFboot