hdiff output

r27871/histogram_reweighting1d.py 2016-07-06 15:37:36.744863093 +0100 r27870/histogram_reweighting1d.py 2016-07-06 15:37:39.828904781 +0100
  1: import numpy as np #to access np.exp() not built int exp  1: import numpy as np #to access np.exp() not built int exp
   2: #import timeseries # for timeseries analysis 
   3: #import commands
   4: #import pdb;
   5: #import pickle
   6: from wham_potential import WhamPotential
   7: #import matplotlib.pyplot as plt
   8: #from matplotlib.pyplot import *
   9: import wham_utils
  2:  10: 
  3: from histogram_reweighting.wham_potential import WhamPotential 
  4: from histogram_reweighting import wham_utils 
  5:  11: 
  6: class Wham1d(object): 
  7:     """Combine 1d histograms of energy at multiple temperatures into one density of states 
  8:  12: 
  9:     Parameters 
 10:     ---------- 
 11:     Tlist : list of floats 
 12:         list of temperatures 
 13:     binenergy : list of floats 
 14:         the lower edges of the energy bins 
 15:     visits1d : ndarray 
 16:         two dimensional array where visits1d[r,e] is the histogram value 
 17:         for replica r at energy bin e. 
 18:     k_B : float 
 19:         the Boltzmann constant 
 20:  13: 
  14: 
  15: 
  16: class wham1d(object):
  17:     """ class to combine 1d histograms of energy E at
  18:     multiple temperatures into one best estimate for the density of states
  19: 
  20:     input will be:
  21:     filenames : list of filenames where the data can be found
  22:     Tlist # Tlist[k] is the temperature of the simulation in filenames[k] 
  23: 
  24:     binenergy = zeros(nebins, float64) #lower edge of bin energy
  25:     visits1d =  zeros([nrep,nebins], integer) #1d histograms of data
 21:     """ 26:     """
 22:     def __init__(self, Tlist, binenergy, visits1d, k_B=1., verbose=True): 27:     def __init__(self, Tlist, binenergy, visits1d):
 23:         if visits1d.shape != (len(Tlist), len(binenergy)): 28:         if visits1d.shape != (len(Tlist), len(binenergy)):
 24:             raise ValueError("visits1d has the wrong shape") 29:             raise ValueError("visits1d has the wrong shape")
 25:         #define some parameters 30:         #define some parameters
 26:         self.k_B = 1. 31:         self.k_B = 1.
 27:  32: 
 28:         self.nrep = len(Tlist) 33:         self.nrep = len(Tlist)
 29:         self.nebins = len(binenergy) 34:         self.nebins = len(binenergy)
 30:         self.Tlist = np.asarray(Tlist) 35:         self.Tlist = np.array(Tlist)
 31:         self.binenergy = np.asarray(binenergy) 36:         self.binenergy = np.array(binenergy)
 32:         self.visits1d = np.asarray(visits1d) 37:         self.visits1d = np.array(visits1d)
 33:          
 34:         self.verbose = verbose 
 35:          
 36:  38: 
 37:     def minimize(self): 39:     def minimize(self):
 38:         """compute the best estimate for the density of states""" 
 39:         nreps = self.nrep 40:         nreps = self.nrep
 40:         nbins = self.nebins 41:         nbins = self.nebins
 41:         visitsT = (self.visits1d) 42:         visitsT = (self.visits1d)
 42:         #print "min vis", np.min(visitsT) 43:         #print "min vis", np.min(visitsT)
  44:         self.logP = np.where( visitsT != 0, np.log( visitsT ), 0 )
 43:         #print "minlogp", np.min(self.logP) 45:         #print "minlogp", np.min(self.logP)
 44:         self.reduced_energy = self.binenergy[np.newaxis,:] / (self.Tlist[:,np.newaxis] * self.k_B) 46:         self.reduced_energy = self.binenergy[np.newaxis,:] / (self.Tlist[:,np.newaxis] * self.k_B)
 45:          47:         
 46:         self.whampot = WhamPotential(visitsT, self.reduced_energy) 48:         self.whampot = WhamPotential(visitsT, self.reduced_energy)
 47:          49:         
 48:         if False: 
 49:             X = np.random.rand( nreps + nbins ) 
 50:         else: 
 51:             # estimate an initial guess for the offsets and density of states 
 52:             # so the minimizer converges more rapidly 
 53:             offsets_estimate, log_dos_estimate = wham_utils.estimate_dos(self.visits1d, 
 54:                                                                          self.reduced_energy) 
 55:             X = np.concatenate((offsets_estimate, log_dos_estimate)) 
 56:  
 57:         E0, grad = self.whampot.getEnergyGradient(X) 
 58:         rms0 = np.linalg.norm(grad) / np.sqrt(grad.size) 
 59:          
 60:         try: 
 61:             from pele.optimize import lbfgs_cpp as quench 
 62:             if self.verbose: 
 63:                 print "minimizing with pele lbfgs" 
 64:             ret = quench(X, self.whampot, tol=1e-3, maxstep=1e4, nsteps=10000) 
 65:         except ImportError: 
 66:             from wham_utils import lbfgs_scipy 
 67:             if self.verbose: 
 68:                 print "minimizing with scipy lbfgs" 
 69:             ret = lbfgs_scipy(X, self.whampot, tol=1e-3, nsteps=10000) 
 70:         #print "quench energy", ret.energy 
 71:          
 72:         if self.verbose: 
 73:             print "chi^2 went from %g (rms %g) to %g (rms %g) in %d iterations" % ( 
 74:                 E0, rms0, ret.energy, ret.rms, ret.nfev) 
 75:          50:         
  51:         X = np.random.rand( nreps + nbins )
  52:         E = self.whampot.getEnergy(X)
  53:         #print "energy", E 
  54: #        self.whampot.test_potential(X)
  55:         
  56:         #print "quenching"
  57:         from wham_utils import lbfgs_scipy
  58:         ret = lbfgs_scipy(X, self.whampot)
  59: #        try: 
  60: #            from pele.optimize import mylbfgs as quench
  61: #            ret = quench(X, self.whampot, iprint=-1, maxstep=1e4)
  62: #        except ImportError:
  63: #            from pele.optimize import lbfgs_scipy as quench
  64: #            ret = quench(X, self.whampot)            
  65:         #print "quench energy", ret.energy
 76:         X = ret.coords 66:         X = ret.coords
  67:         
 77:         self.logn_E = X[nreps:] 68:         self.logn_E = X[nreps:]
 78:         self.w_i_final = X[:nreps] 69:         self.w_i_final = X[:nreps]
 79:          70:         
 80:  71: 
 81:     def calc_Cv(self, ndof, Tlist=None, ntemp=100): 72: #    def calc_Cv_new(self, NDOF, TRANGE=None, NTEMP=100):
 82:         if Tlist is None: 73: #        from pele.thermodynamics import dos_to_cv
 83:             Tlist = np.linspace(self.Tlist[0], self.Tlist[-1], ntemp) 74: #        dT = (self.Tlist[-1] - self.Tlist[0]) / NTEMP
 84:         have_data = np.where(self.visits1d.sum(0) > 0)[0] 75: #        Tlist = np.arange(self.Tlist[0], self.Tlist[-1], dT)
 85:         return wham_utils.calc_Cv(Tlist, self.binenergy, self.logn_E, ndof, 76: ##        print self.logn_E
 86:                                   have_data=have_data, k_B=self.k_B) 77: #        lZ, U, U2, Cv = dos_to_cv(self.binenergy, self.logn_E, Tlist, K=NDOF)
 87:                                    78: #        cvdata = np.zeros([len(Tlist), 6])
  79: #        cvdata[:,0] = Tlist
  80: #        cvdata[:,1] = lZ
  81: #        cvdata[:,2] = U # average potential energy
  82: #        cvdata[:,3] = U2
  83: #        cvdata[:,5] = Cv
  84: #        
  85: #        eavg = U + float(NDOF) / 2 * Tlist # average energy including the kinetic degrees of freedom
  86: #        cvdata[:,4] = eavg
  87: #        return cvdata
  88:         
  89: 
  90:     def calc_Cv(self, NDOF, TRANGE=None, NTEMP=100, use_log_sum=None):
  91: #        return self.calc_Cv_new(NDOF, TRANGE, NTEMP)
  92:         return wham_utils.calc_Cv(self.logn_E, self.visits1d, self.binenergy,
  93:                 NDOF, self.Tlist, self.k_B, TRANGE, NTEMP, use_log_sum=use_log_sum)
  94: 
  95: 
 88:  96: 


r27871/histogram_reweighting2d.py 2016-07-06 15:37:37.096867855 +0100 r27870/histogram_reweighting2d.py 2016-07-06 15:37:40.152909170 +0100
  1: import numpy as np 
  2: import scipy 
  3: from numpy import exp, log  1: from numpy import exp, log
  4:   2: import numpy as np #to access np.exp() not built int exp
   3: #from math import *
   4: #import scipy.optimize
   5: #import timeseries # for timeseries analysis 
  5: import wham_utils  6: import wham_utils
   7: from wham_utils import logSum, logSum1
   8: #import matplotlib.pyplot as plt
   9: #from matplotlib.pyplot import *
  10: #mbar = pickle.load(open("mbar.pickle","rb"))
  6:  11: 
  7:  12: 
  8:  13: 
  9:  14: 
 10:  15: 
 11:  16: 
 12: class Wham2d(object): 17: class wham2d(object):
 13:     """ class to combine 2d histograms of energy E and order parameter q at 18:     """ class to combine 2d histograms of energy E and order parameter q at
 14:     multiple temperatures into one best estimate for the histogram 19:     multiple temperatures into one best estimate for the histogram
 15:    20:   
 16:     input will be: 21:     input will be:
 17:     filenames : list of filenames where the data can be found 22:     filenames : list of filenames where the data can be found
 18:     Tlist # Tlist[k] is the temperature of the simulation in filenames[k]  23:     Tlist # Tlist[k] is the temperature of the simulation in filenames[k] 
 19:    24:   
 20:     binenergy = zeros(nebins, float64) #lower edge of bin energy 25:     binenergy = zeros(nebins, float64) #lower edge of bin energy
 21:     binq =      zeros(nqbins, float64) #lower edge of q bins 26:     binq =      zeros(nqbins, float64) #lower edge of q bins
 22:     visits2d =  zeros([nrep,nebins,nqbins], integer) #2d histogram of data 27:     visits2d =  zeros([nrep,nebins,nqbins], integer) #2d histogram of data
 23:     """ 28:     """
 24:     #============================================================================================= 29:     #=============================================================================================
 25:     # Constructor. 30:     # Constructor.
 26:     #============================================================================================= 31:     #=============================================================================================
 27:     def __init__(self, Tlist, binenergy, binq, visits2d, k_B=1., verbose=True): 32:     def __init__(self, Tlist, binenergy, binq, visits2d):
 28:         if visits2d.shape != (len(Tlist), len(binenergy), len(binq)): 
 29:             raise ValueError("Visits2d has the wrong shape") 
 30:      33:     
 31:         #define some parameters 34:         #define some parameters
 32:         self.k_B = k_B 35:         self.k_B=1.
 33:         self.LOGMIN = -1e200 36:         self.LOGMIN = -1e200
 34:         self.verbose = verbose 
 35:      37:     
 36:         self.nrep = len(Tlist) 38:         self.nrep = len(Tlist)
 37:         self.nebins = len(binenergy) 39:         self.nebins = len(binenergy)
 38:         self.nqbins = len(binq) 40:         self.nqbins = len(binq)
 39:         self.Tlist = np.array(Tlist) 41:         self.Tlist = np.array(Tlist)
 40:         self.binenergy = np.array(binenergy) 42:         self.binenergy = np.array(binenergy)
 41:         self.binq = np.array(binq) 43:         self.binq = np.array(binq)
 42:         self.visits2d = np.array(visits2d, dtype = np.integer) 44:         self.visits2d = np.array(visits2d, dtype = np.integer)
 43:          45:         
 44:     def minimize(self): 46:     def minimize(self):
 45:         #shape(visits2d) is now (nqbins, nebins, nreps) 47:         #shape(visits2d) is now (nqbins, nebins, nreps)
 46:         #we need it to be (nreps, nqbins*nebins) 48:         #we need it to be (nreps, nqbins*nebins)
 47:         #first reorder indices 49:         #first reorder indices
 48:         nreps = self.nrep 50:         nreps = self.nrep
 49:         nebins = self.nebins 51:         nebins = self.nebins
 50:         nqbins = self.nqbins 52:         nqbins = self.nqbins
 51:         nbins = self.nebins * self.nqbins 53:         nbins = self.nebins * self.nqbins
 52:          54:         #visits = np.zeros([nreps, nebins, nqbins], np.integer)
 53:         reduced_energy = np.zeros([nreps, nebins, nqbins]) 55:         reduced_energy = np.zeros([nreps, nebins, nqbins])
  56: #        for k in range(self.nrep):
  57: #            for j in range(self.nqbins):
  58: #                for i in range(self.nebins):
  59: #                    #visits[k,i,j] = self.visits2d[i,j,k]
  60: #                    reduced_energy[k,i,j] = self.binenergy[i] / (self.Tlist[k]*self.k_B)
 54:         for j in range(self.nqbins): 61:         for j in range(self.nqbins):
 55:             reduced_energy[:,:,j] = self.binenergy[np.newaxis,:] / (self.Tlist[:,np.newaxis]*self.k_B) 62:             reduced_energy[:,:,j] = self.binenergy[np.newaxis,:] / (self.Tlist[:,np.newaxis]*self.k_B)
 56:                      63:                     
 57:         visits = self.visits2d 64:         visits = self.visits2d
 58:         visits = np.reshape( visits, [nreps, nbins ])  65:         visits = np.reshape( visits, [nreps, nbins ]) 
 59:         reduced_energy = np.reshape( reduced_energy, [nreps, nbins])            66:         reduced_energy = np.reshape( reduced_energy, [nreps, nbins])           
 60:         self.logP = np.where( visits != 0, np.log( visits.astype(float) ), 0 ) 67:         self.logP = np.where( visits != 0, np.log( visits.astype(float) ), 0 )
 61:  68: 
 62:          69:         
 63:         from wham_potential import WhamPotential 70:         from wham_potential import WhamPotential
 64:         whampot = WhamPotential(visits, reduced_energy ) 71:         whampot = WhamPotential(visits, reduced_energy )
 65:          72:         
 66:         nvar = nbins + nreps 73:         nvar = nbins + nreps
 67:         if False: 74:         X = np.random.rand(nvar)
 68:             X = np.random.rand(nvar) 
 69:         else: 
 70:             # estimate an initial guess for the offsets and density of states 
 71:             # so the minimizer converges more rapidly 
 72:             offsets_estimate, log_dos_estimate = wham_utils.estimate_dos(visits, 
 73:                                                                          reduced_energy) 
 74:             X = np.concatenate((offsets_estimate, log_dos_estimate)) 
 75:             assert X.size == nvar 
 76:  
 77:         E0, grad = whampot.getEnergyGradient(X) 
 78:         rms0 = np.linalg.norm(grad) / np.sqrt(grad.size) 
 79:  
 80:  
 81:         from wham_utils import lbfgs_scipy 75:         from wham_utils import lbfgs_scipy
 82:         ret = lbfgs_scipy(X, whampot) 76:         ret = lbfgs_scipy(X, whampot)
 83: #        print "initial energy", whampot.getEnergy(X) 77: #        print "initial energy", whampot.getEnergy(X)
 84: #        try:  78: #        try: 
 85: #            from pele.optimize import mylbfgs as quench 79: #            from pele.optimize import mylbfgs as quench
 86: #            ret = quench(X, whampot, iprint=10, maxstep = 1e4) 80: #            ret = quench(X, whampot, iprint=10, maxstep = 1e4)
 87: #        except ImportError: 81: #        except ImportError:
 88: #            from pele.optimize import lbfgs_scipy as quench 82: #            from pele.optimize import lbfgs_scipy as quench
 89: #            ret = quench(X, whampot)             83: #            ret = quench(X, whampot)            
 90:  84: 
 91:         if self.verbose: 85:         print "quenched energy", ret.energy
 92:             print "chi^2 went from %g (rms %g) to %g (rms %g) in %d iterations" % ( 86:         
 93:                 E0, rms0, ret.energy, ret.rms, ret.nfev) 
 94:  
 95:          87:         
 96:         #self.logn_Eq = zeros([nebins,nqbins], float64) 88:         #self.logn_Eq = zeros([nebins,nqbins], float64)
 97:         X = ret.coords 89:         X = ret.coords
 98:         self.logn_Eq = X[nreps:] 90:         self.logn_Eq = X[nreps:]
 99:         self.w_i_final = X[:nreps] 91:         self.w_i_final = X[:nreps]
100:          92:         
101:         self.logn_Eq = np.reshape(self.logn_Eq, [nebins, nqbins]) 93:         self.logn_Eq = np.reshape(self.logn_Eq, [nebins, nqbins])
102:         self.logn_Eq = np.where( self.visits2d.sum(0) == 0, self.LOGMIN, self.logn_Eq ) 94:         self.logn_Eq = np.where( self.visits2d.sum(0) == 0, self.LOGMIN, self.logn_Eq )
103:          95:         
104:         #renormalize logn_Eq 96:         #renormalize logn_Eq
153:             for i in range(nebins):145:             for i in range(nebins):
154:                 logP_Eq[i,:] = logn_Eq[i,:]-(binenergy[i] - EREF)/(self.k_B*T)146:                 logP_Eq[i,:] = logn_Eq[i,:]-(binenergy[i] - EREF)/(self.k_B*T)
155:       147:       
156:             logP_Eq[self.allzero2dind[0], self.allzero2dind[1]] = self.LOGMIN148:             logP_Eq[self.allzero2dind[0], self.allzero2dind[1]] = self.LOGMIN
157:             expoffset = np.nanmax(logP_Eq)149:             expoffset = np.nanmax(logP_Eq)
158:             #print "T expoffset ", T, expoffset150:             #print "T expoffset ", T, expoffset
159:             logP_Eq -= expoffset151:             logP_Eq -= expoffset
160:             #P_q = np.exp(logP_Eq).sum(0)152:             #P_q = np.exp(logP_Eq).sum(0)
161:             # sum over the energy153:             # sum over the energy
162:             for j in range(nqbins):154:             for j in range(nqbins):
163:                 logP_q[j] = scipy.misc.logsumexp( logP_Eq[:,j] )155:                 logP_q[j] = logSum( logP_Eq[:,j] )
164:             logP_q[self.nodataq] = np.NaN156:             logP_q[self.nodataq] = np.NaN
165:             F_q[:,n] = -self.k_B*T*logP_q[:]157:             F_q[:,n] = -self.k_B*T*logP_q[:]
166:             fmin = np.nanmin(F_q[:,n])158:             fmin = np.nanmin(F_q[:,n])
167:             F_q[:,n] -= fmin159:             F_q[:,n] -= fmin
168:     160:     
169:         return TRANGE,F_q161:         return TRANGE,F_q
170:   162:   
171:     def calc_qavg(self, TRANGE = []):163:     def calc_qavg(self, TRANGE = []):
172:         """calculate the average q as a function of temperature"""164:         """calculate the average q as a function of temperature"""
173:         #put some variables in this namespace165:         #put some variables in this namespace
208:             for i in range(nebins):200:             for i in range(nebins):
209:                 logP_Eq[i,:] = logn_Eq[i,:]-(binenergy[i] - EREF)/(self.k_B*T)201:                 logP_Eq[i,:] = logn_Eq[i,:]-(binenergy[i] - EREF)/(self.k_B*T)
210:       202:       
211:             logP_Eq[self.allzero2dind[0], self.allzero2dind[1]] = self.LOGMIN203:             logP_Eq[self.allzero2dind[0], self.allzero2dind[1]] = self.LOGMIN
212:             expoffset = logP_Eq.max()204:             expoffset = logP_Eq.max()
213:             #print "T expoffset ", T, expoffset205:             #print "T expoffset ", T, expoffset
214:             logP_Eq -= expoffset206:             logP_Eq -= expoffset
215:             #P_q = np.exp(logP_Eq).sum(0)207:             #P_q = np.exp(logP_Eq).sum(0)
216:             # sum over the energy208:             # sum over the energy
217:             for j in range(nqbins):209:             for j in range(nqbins):
218:                 logP_q[j] = scipy.misc.logsumexp(logP_Eq[:,j])210:                 logP_q[j] = wham_utils.logSum( logP_Eq[:,j] )
219:             logP_q[self.nodataq] = np.NaN211:             logP_q[self.nodataq] = np.NaN
220:       212:       
221:             #find mean q213:             #find mean q
222:             qmin = min(binq)214:             qmin = min(binq)
223:             qmin -= 0.1215:             qmin -= 0.1
224:             lqavg = -1.0e30216:             lqavg = -1.0e30
225:             lnorm = -1.0e30217:             lnorm = -1.0e30
226:             for i in range(0,nqbins): 218:             for i in range(0,nqbins): 
227:                 if not np.isnan(logP_q[i]):219:                 if not np.isnan(logP_q[i]):
228:                     lnorm = np.logaddexp(lnorm, logP_q[i]) 220:                     lnorm = wham_utils.logSum1( lnorm, logP_q[i] ) 
229:                     lqavg = np.logaddexp(lqavg, logP_q[i] + log(binq[i] - qmin))221:                     lqavg = wham_utils.logSum1( lqavg, logP_q[i] + log(binq[i] - qmin) )
230:             self.qavg[n] = exp(lqavg - lnorm) + qmin222:             self.qavg[n] = exp(lqavg - lnorm) + qmin
231:             #print lqavg223:             #print lqavg
232:     224:     
233:         return TRANGE, self.qavg225:         return TRANGE,self.qavg
234:   226:   
235:   227:   
236:     def calc_Cv(self, ndof, Tlist=None, ntemp=100):228:     def calc_Cv(self, NDOF):
 229:         nebins = self.nebins
237:         visits1d = self.visits2d.sum(2)230:         visits1d = self.visits2d.sum(2)
238:         have_data = np.where(visits1d.sum(0) > 0)[0]231:         logn_E = np.zeros(nebins)
239:         logn_E = scipy.misc.logsumexp(self.logn_Eq, axis=1)232:         for i in range(nebins):
240: #        logn_E = np.zeros(nebins)233:             logn_E[i] = wham_utils.logSum(self.logn_Eq[i,:])
241: #        for i in range(nebins):234: 
242: #            logn_E[i] = scipy.misc.logsumexp(self.logn_Eq[i,:])235:         return wham_utils.calc_Cv(logn_E, visits1d, self.binenergy, \
243: 236:                                   NDOF, self.Tlist, self.k_B)
244:         if Tlist is None:237: 
245:             Tlist = np.linspace(self.Tlist[0], self.Tlist[-1], ntemp)238:   
246:         return wham_utils.calc_Cv(Tlist, self.binenergy, logn_E, ndof, 
247:                                   have_data=have_data, k_B=self.k_B) 
248: 239: 


r27871/makecvplots.py 2016-07-06 15:37:39.488900195 +0100 r27870/makecvplots.py 2016-07-06 15:37:42.692943384 +0100
369:     time1 = Vall[0][ind].time369:     time1 = Vall[0][ind].time
370:     # main loop over indices data370:     # main loop over indices data
371:     for isub in indices[1:]:371:     for isub in indices[1:]:
372:         # compute the cv for all data collected between time2 and time1372:         # compute the cv for all data collected between time2 and time1
373:         time2 = Vall[0][isub].time373:         time2 = Vall[0][isub].time
374:         # print "time 2 = ", time2374:         # print "time 2 = ", time2
375:         print "computing Cv for the time period", time2, "to", time1 375:         print "computing Cv for the time period", time2, "to", time1 
376: 376: 
377:         # visits corresponding to isub will be subtracted from ind  377:         # visits corresponding to isub will be subtracted from ind  
378:         binenergy, visits = makeHistogram(Vall, ind, isub, neglect=neglect, verbose=verbose)378:         binenergy, visits = makeHistogram(Vall, ind, isub, neglect=neglect, verbose=verbose)
379:         379: 
380:         assert(not np.any(np.isnan(visits))) 
381:         if useWHAM:380:         if useWHAM:
382:             wham = WHAM.Wham1d(Tlist, binenergy, visits)381:             wham = WHAM.wham1d(Tlist, binenergy, visits)
383:             wham.minimize()382:             wham.minimize()
384:             cvdata = wham.calc_Cv(NDOF, ntemp=500)383:             cvdata = wham.calc_Cv(NDOF, NTEMP=500, use_log_sum=True)
385:         else:384:         else:
386:             cvdata = calcCvNoWHAM(binenergy, visits, Tlist, NDOF)385:             cvdata = calcCvNoWHAM(binenergy, visits, Tlist, NDOF)
387:         assert(not np.any(np.isnan(cvdata))) 
388:         #checkWHAM(binenergy, visits, cvdata, Tlist)386:         #checkWHAM(binenergy, visits, cvdata, Tlist)
389:         datalist.append( (time1, time2, cvdata) )387:         datalist.append( (time1, time2, cvdata) )
390: 388: 
391:     if len(indices) < ncurves:389:     if len(indices) < ncurves:
392:         # plot Cv for the last Visits file without any subtraction.390:         # plot Cv for the last Visits file without any subtraction.
393:         # Set time2 to None to indicate that no subtraction was done391:         # Set time2 to None to indicate that no subtraction was done
394:         print "computing Cv from all existing data" 392:         print "computing Cv from all existing data" 
395:         binenergy, visits = makeHistogram(Vall, ind, None)393:         binenergy, visits = makeHistogram(Vall, ind, None)
396:         if useWHAM:394:         if useWHAM:
397:             wham = WHAM.Wham1d(Tlist, binenergy, visits)395:             wham = WHAM.wham1d(Tlist, binenergy, visits)
398:             wham.minimize()396:             wham.minimize()
399:             cvdata = wham.calc_Cv(NDOF, ntemp=500)397:             cvdata = wham.calc_Cv(NDOF, NTEMP=500, use_log_sum=True)
400:         else:398:         else:
401:             cvdata = calcCvNoWHAM(binenergy, visits, Tlist, NDOF)399:             cvdata = calcCvNoWHAM(binenergy, visits, Tlist, NDOF)
402:         datalist.append( (time1, None, cvdata) )400:         datalist.append( (time1, None, cvdata) )
403: 401: 
404: 402: 
405:     if True:403:     if True:
406:         # save the data to a file404:         # save the data to a file
407:         for time1, time2, cvdata in datalist:405:         for time1, time2, cvdata in datalist:
408:             if time2 is None:406:             if time2 is None:
409:                 fname = "%s%d" % (outprefix, time1)407:                 fname = "%s%d" % (outprefix, time1)


r27871/test_reweighting_2d.py 2016-07-06 15:37:37.804877426 +0100 r27870/test_reweighting_2d.py 2016-07-06 15:37:40.824918264 +0100
  1: import unittest  1: svn: E195012: Unable to find repository location for 'svn+ssh://svn.ch.private.cam.ac.uk/groups/wales/trunk/SCRIPTS/python/ptmc/histogram_reweighting/test_reweighting_2d.py' in revision 27870
  2: import numpy as np 
  3:  
  4: from matplotlib import pyplot as plt 
  5:  
  6: import test_reweighting 
  7: from histogram_reweighting.histogram_reweighting2d import Wham2d 
  8:  
  9: class HarmonicOscillatorQ(test_reweighting.HarmonicOscillator): 
 10:     """add routines to build a 2d visits histogram""" 
 11:      
 12:      
 13:     def random_histogram_2d(self, N, binenergy, binq, T): 
 14:         binenergy = np.append(binenergy, binenergy[-1] + (binenergy[-1] - binenergy[-2])) 
 15:         binq = np.append(binq, binq[-1] + (binq[-1] - binq[-2])) 
 16:  
 17:         x = np.random.normal(size=self.d * N) 
 18:         x = x.reshape([N,self.d]) 
 19:         x2 = np.sum(x**2, axis=1) 
 20:         assert x2.size == N 
 21:         Elist = x2 * T / 2 
 22:         xcoord = x[:,0] # the first coordinate 
 23:          
 24:         evar = np.var(Elist) 
 25:         print "Cv computed from energy variance", evar / T**2 + float(self.d)/2 
 26:          
 27:         counts, bins1, bins2 = np.histogram2d(Elist, xcoord, bins=[binenergy, binq]) 
 28:          
 29:         if False: 
 30:             c = np.where(counts==0, np.nan, counts) 
 31:             plt.imshow(c) 
 32:             plt.show() 
 33:              
 34:             v = counts.sum(1) 
 35:             plt.clf() 
 36:             plt.plot(binenergy[:-1], v) 
 37:             plt.show() 
 38:  
 39:         return counts 
 40:  
 41:     def make_visits_2d(self, Tlist, binenergy, binq, N, random=True): 
 42:         visits = [] 
 43:         for i, T in enumerate(Tlist): 
 44:             if random: 
 45:                 counts = self.random_histogram_2d(N, binenergy, binq, T) 
 46:             else: 
 47:                 counts = self.make_analytic_histogram_2d(N, binenergy, binq, T) 
 48:             visits.append(counts) 
 49:  
 50:         if False: 
 51:             plt.clf() 
 52:             plt.title("before numpy array") 
 53:             for T, counts in izip(Tlist, visits): 
 54:                 log_nE = np.log(counts) + binenergy / T 
 55:                 plt.plot(binenergy, log_nE) 
 56:             plt.show() 
 57:  
 58:         visits = np.array(visits) 
 59:         print visits.shape,  (len(Tlist), len(binenergy), len(binq)) 
 60:         assert visits.shape == (len(Tlist), len(binenergy), len(binq)) 
 61:  
 62:  
 63:         if False: 
 64:             plt.clf() 
 65:             print visits.shape, binenergy.shape, Tlist.shape 
 66:             log_nET = np.log(visits) + binenergy[np.newaxis, :] / Tlist[:,np.newaxis] #+ wham.w_i_final[:,np.newaxis] 
 67:             plt.plot(binenergy, np.transpose(log_nET)) 
 68:             plt.show() 
 69: #        raise Exception 
 70:  
 71:          
 72:         return visits 
 73:          
 74:          
 75: class TestReweighting2d(unittest.TestCase): 
 76:     def setUp(self): 
 77:         self.d = 3 
 78:         self.N = 100000 
 79:         self.Tlist = test_reweighting.get_temps(.2, 1.6, 8) 
 80:         self.Tlist = np.array(self.Tlist) 
 81:         self.binenergy = np.linspace(0, 20, 100) 
 82:         self.binq = np.linspace(0, 7, 20) 
 83:  
 84:     def test(self): 
 85:         Tlist = self.Tlist 
 86:         binenergy = self.binenergy 
 87:         binq = self.binq 
 88:         ho = HarmonicOscillatorQ(self.d) 
 89:         visits2d = ho.make_visits_2d(self.Tlist, self.binenergy, self.binq, self.N, random=True) 
 90:          
 91:         wham = Wham2d(Tlist, binenergy, binq, visits2d) 
 92:         wham.minimize() 
 93:         cvdata = wham.calc_Cv(self.d, Tlist=self.Tlist) 
 94:         print cvdata[:,5] 
 95:         for cv in cvdata[:,5]: 
 96:             self.assertAlmostEqual(cv, self.d, delta=.7) 
 97:          
 98: #    def test_qavg(self): 
 99: #        ho = HarmonicOscillatorQ(self.d) 
100: #        visits2d = ho.make_visits_2d(self.Tlist, self.binenergy, self.binq, self.N, random=True) 
101: #        wham = Wham2d(self.Tlist, self.binenergy, self.binq, visits2d) 
102: #        wham.minimize() 
103: #         
104: #        qavg = wham.calc_qavg(self.Tlist) 
105:  
106:              
107:          
108:          
109: if __name__ == "__main__": 
110:     unittest.main() 


r27871/test_reweighting.py 2016-07-06 15:37:37.464872828 +0100 r27870/test_reweighting.py 2016-07-06 15:37:40.476913594 +0100
  1: import unittest  1: svn: E195012: Unable to find repository location for 'svn+ssh://svn.ch.private.cam.ac.uk/groups/wales/trunk/SCRIPTS/python/ptmc/histogram_reweighting/test_reweighting.py' in revision 27870
  2: from itertools import izip 
  3: import numpy as np 
  4:  
  5: import matplotlib.pyplot as plt 
  6:  
  7: from histogram_reweighting1d import Wham1d 
  8: from histogram_reweighting import wham_utils 
  9:  
 10: def get_temps(Tmin, Tmax, nreplicas): 
 11:     """ 
 12:     set up the temperatures 
 13:     distribute them exponentially 
 14:     """ 
 15:     #dT = (Tmax - Tmin) / (nreplicas-1) 
 16:     CTE = np.exp( np.log( Tmax / Tmin ) / (nreplicas-1) ) 
 17:     Tlist = [Tmin* CTE**i for i in range(nreplicas)] 
 18:     return Tlist 
 19:  
 20:  
 21: class HarmonicOscillator(object): 
 22:     def __init__(self, d): 
 23:         self.d = d 
 24:      
 25:     def random_energy(self, T): 
 26:         x = np.random.normal(size=self.d) 
 27:         x2 = x.dot(x) 
 28:         E = T * x2 / 2 
 29:         return E 
 30:  
 31:     def random_energies(self, N, T): 
 32:         x = np.random.normal(size=self.d * N) 
 33:         x = x.reshape([N,self.d]) 
 34:         x2 = np.sum(x**2, axis=1) 
 35:         assert x2.size == N 
 36:         Elist = x2 * T / 2 
 37:         return Elist 
 38:  
 39:     def make_analytic_histogram(self, N, binenergy, T): 
 40:         assert self.d > 2 
 41:         binenergy = np.append(binenergy, binenergy[-1] + (binenergy[-1] - binenergy[-2])) 
 42:         ecenters = (binenergy[1:] + binenergy[:-1]) / 2 
 43:         assert np.all(ecenters > 0) 
 44:         logP = float(self.d-2) / 2 * np.log(ecenters) - ecenters / T 
 45:          
 46:               
 47:          
 48:         # now normalize it to N 
 49:         import scipy 
 50:         logsumP = scipy.misc.logsumexp(logP) 
 51:         logP = logP + np.log(N) - logsumP 
 52:         visits = np.exp(logP) 
 53: #        visits = np.round(visits) 
 54:         if True: 
 55:             e1 = np.average(ecenters, weights=visits)  
 56:             e2 = np.average(ecenters**2, weights=visits) 
 57:             print "Cv from analytic histogram", (e2 - e1**2) / T**2 + float(self.d)/2 
 58:          
 59: #        if True: 
 60: #            plt.clf() 
 61: #            plt.plot(binenergy[:-1], visits.transpose()) 
 62: #            plt.show() 
 63:          
 64:         return visits  
 65:          
 66:  
 67:     def random_histogram(self, N, binenergy, T): 
 68:         Elist = self.random_energies(N, T) 
 69:      
 70:         evar = np.var(Elist) 
 71:         print "Cv computed from energy variance", evar / T**2 + float(self.d)/2 
 72:          
 73:         binenergy = np.append(binenergy, binenergy[-1]+(binenergy[-1] - binenergy[-2])) 
 74:          
 75:         counts, bins = np.histogram(Elist, binenergy) 
 76:          
 77:         if True: 
 78:             e = binenergy[:-1] 
 79: #            plt.plot(e, counts) 
 80:             plt.plot(e, np.log(counts) + e / T ) 
 81:         return counts 
 82:      
 83:     def make_visits(self, Tlist, binenergy, N, random=True): 
 84: #        binenergy = np.linspace(0, 40, 1000) 
 85: #        visits = np.zeros([len(Tlist), len(binenergy)]) 
 86:         visits = [] 
 87:         for i, T in enumerate(Tlist): 
 88:             if random: 
 89:                 counts = self.random_histogram(N, binenergy, T) 
 90:             else: 
 91:                 counts = self.make_analytic_histogram(N, binenergy, T) 
 92:             visits.append(counts) 
 93: #        plt.show() 
 94:  
 95:         if False: 
 96:             plt.clf() 
 97:             plt.title("before numpy array") 
 98:             for T, counts in izip(Tlist, visits): 
 99:                 log_nE = np.log(counts) + binenergy / T 
100:                 plt.plot(binenergy, log_nE) 
101:             plt.show() 
102:  
103:         visits = np.array(visits) 
104:  
105:  
106:         if False: 
107:             plt.clf() 
108:             print visits.shape, binenergy.shape, Tlist.shape 
109:             log_nET = np.log(visits) + binenergy[np.newaxis, :] / Tlist[:,np.newaxis] #+ wham.w_i_final[:,np.newaxis] 
110:             plt.plot(binenergy, np.transpose(log_nET)) 
111:             plt.show() 
112: #        raise Exception 
113:  
114:          
115:         return visits 
116:      
117:      
118:     def make_analytic_dos(self, binenergy): 
119:         """make a binned density of states 
120:          
121:             dos = E**(d/2) 
122:         """ 
123:         assert self.d > 2 
124:         binenergy = np.append(binenergy, binenergy[-1] + (binenergy[-1] - binenergy[-2])) 
125:         ecenters = (binenergy[1:] + binenergy[:-1]) / 2 
126:         assert np.all(ecenters > 0) 
127:         log_dos = float(self.d-2) / 2 * np.log(ecenters) 
128:         return log_dos 
129:  
130:  
131: class TestHistogramReweighting(unittest.TestCase): 
132:     def setUp(self): 
133:         self.d = 3 
134:         self.N = 100000 
135:         self.Tlist = get_temps(.2, 1.6, 8) 
136: #        self.Tlist = [2.000000000000000111e-01, 
137: #                2.691800385264712658e-01, 
138: #                3.622894657055626966e-01, 
139: #                4.876054616817901977e-01, 
140: #                6.562682848061104357e-01, 
141: #                8.832716109390498227e-01, 
142: #                1.188795431309558781e+00, 
143: #                1.600000000000000089e+00] 
144:         self.Tlist = np.array(self.Tlist) 
145:         self.binenergy = np.linspace(0, 20, 1000) 
146:  
147:     def test1(self): 
148:         np.random.seed(0) 
149:         ho = HarmonicOscillator(self.d) 
150:         self.visits = ho.make_visits(self.Tlist, self.binenergy, self.N, random=True) 
151:         assert self.visits.shape == (len(self.Tlist), len(self.binenergy)) 
152:         visits = self.visits 
153:         binenergy = self.binenergy 
154:         Tlist = self.Tlist 
155:         if False: 
156:             plt.clf() 
157:             print visits.shape, binenergy.shape, Tlist.shape 
158:             log_nET = np.log(visits) + binenergy[np.newaxis, :] / Tlist[:,np.newaxis] #+ wham.w_i_final[:,np.newaxis] 
159:             plt.plot(binenergy, np.transpose(log_nET)) 
160: #            for T, counts in izip(Tlist, visits): 
161: #                log_nE = np.log(counts) + binenergy / T 
162: #                plt.plot(binenergy, log_nE) 
163:             plt.show() 
164:          
165:         wham = Wham1d(Tlist, binenergy, visits.copy()) 
166:         wham.minimize() 
167:         cvdata = wham.calc_Cv(3, Tlist=Tlist) 
168: #        print cvdata.shape 
169: #        print cvdata 
170: #        print cvdata 
171: #        print "Cv values", cvdata[:,5] 
172:          
173:         for cv in cvdata[:,5]: 
174:             self.assertAlmostEqual(cv, self.d, delta=.3) 
175:  
176:         if False: 
177:             plt.clf() 
178:             plt.plot(Tlist, cvdata[:,5]) 
179:             plt.show() 
180:          
181:          
182:          
183:         if False: 
184:             plt.clf() 
185:             plt.plot(binenergy, wham.logn_E) 
186:             plt.show() 
187:         if False: 
188:             plt.clf() 
189:             for r, T in enumerate(Tlist): 
190:                 v = visits[r,:] 
191:     #            plot(bin,v) 
192:             plt.plot(binenergy, np.log(np.transpose(visits))) 
193:             plt.show() 
194:         if False: 
195:             plt.clf() 
196:             log_nET = np.log(visits) + binenergy[np.newaxis, :]  / Tlist[:,np.newaxis] + wham.w_i_final[:,np.newaxis] 
197:             plt.plot(binenergy, np.transpose(log_nET)) 
198:             plt.plot(binenergy, wham.logn_E) 
199:             plt.show() 
200:         if False: 
201:             plt.clf() 
202:             plt.plot(binenergy, wham.logn_E, 'k', lw=2) 
203:             newlogn_E = wham_utils.dos_from_offsets(Tlist, binenergy, visits, 
204:                                                     wham.w_i_final) 
205:             plt.plot(binenergy, newlogn_E, '--r', lw=.5) 
206:             plt.show() 
207:      
208:     def test_analytic(self): 
209:         ho = HarmonicOscillator(self.d) 
210:         visits = ho.make_visits(self.Tlist, self.binenergy, self.N, random=False) 
211:          
212:         wham = Wham1d(self.Tlist, self.binenergy, visits) 
213:         wham.minimize() 
214:         cvdata = wham.calc_Cv(3, Tlist=self.Tlist) 
215:          
216:         for cv in cvdata[:,5]: 
217:             self.assertAlmostEqual(cv, self.d, delta=.01) 
218:  
219:      
220:     def test2(self): 
221:         ho = HarmonicOscillator(self.d) 
222:         self.visits = ho.make_visits(self.Tlist, self.binenergy, self.N, random=True) 
223:         assert self.visits.shape == (len(self.Tlist), len(self.binenergy)) 
224:         self.reduced_energy = self.binenergy[np.newaxis,:] / (self.Tlist[:,np.newaxis]) 
225:  
226:         wham_utils.estimate_dos(self.visits, self.visits) 
227:              
228:  
229: if __name__ == "__main__": 
230:     unittest.main() 


r27871/test_utils.py 2016-07-06 15:37:38.140881969 +0100 r27870/test_utils.py 2016-07-06 15:37:41.200923315 +0100
  1: import unittest  1: svn: E195012: Unable to find repository location for 'svn+ssh://svn.ch.private.cam.ac.uk/groups/wales/trunk/SCRIPTS/python/ptmc/histogram_reweighting/test_utils.py' in revision 27870
  2:  
  3: import numpy as np 
  4:  
  5: from histogram_reweighting import wham_utils 
  6: import test_reweighting 
  7:  
  8:  
  9: class TestCvCalc(unittest.TestCase): 
 10:     def test(self): 
 11:         d = 3 
 12:         ho = test_reweighting.HarmonicOscillator(d) 
 13:         binenergy = np.linspace(0, 40, 200) 
 14:         Tlist = np.linspace(.2, 1.6, 8) 
 15:         log_dos = ho.make_analytic_dos(binenergy) 
 16:          
 17:         cvdata = wham_utils.calc_Cv(Tlist, binenergy, log_dos, d) 
 18:         cv = cvdata[:,5] 
 19:         self.assertLess(np.max(np.abs(cv - d)), 0.05) 
 20:          
 21:     def test_with_zeros(self): 
 22:         d = 3 
 23:         ho = test_reweighting.HarmonicOscillator(d) 
 24:         binenergy = np.linspace(0, 40, 200) 
 25:         Tlist = np.linspace(.2, 1.6, 8) 
 26:         log_dos = ho.make_analytic_dos(binenergy) 
 27:         have_data = np.array(range(len(binenergy)-100)) 
 28: #        log_dos[0] = 0. 
 29:         log_dos[-100:] = 0. 
 30:          
 31:         cvdata = wham_utils.calc_Cv(Tlist, binenergy, log_dos, d, have_data=have_data) 
 32:         cv = cvdata[:,5] 
 33:         self.assertLess(np.max(np.abs(cv - d)), 0.05) 
 34:          
 35: if __name__ == "__main__": 
 36:     unittest.main() 


r27871/wham_potential.py 2016-07-06 15:37:38.488886596 +0100 r27870/wham_potential.py 2016-07-06 15:37:41.720930391 +0100
  1: import numpy as np  1: import numpy as np
   2: #from pele.potentials.potential import BasePotential
  2:   3: 
  3: from pele.potentials import BasePotential  4: from pele.potentials import BasePotential
  4:  
  5: class WhamPotential(BasePotential):  5: class WhamPotential(BasePotential):
  6:     """  6:     """
  7:     the idea behind this minimization procedure is as follows   7:     the idea behind this minimization procedure is as follows 
  8:       8:     
  9:     from a simulation at temperature T you find the probability of finding energy  9:     from a simulation at temperature T you find the probability of finding energy
 10:     E is P(E,T).  We know this can be compared to the density of states n(E) as 10:     E is P(E,T).  We know this can be compared to the density of states n(E) as
 11:      11:     
 12:         P(E,T) = n(E) exp(-E/T_i) / w_i 12:         P(E,T) = n(E) exp(-E/T_i) / w_i
 13:      13:     
 14:     Where w_i is a constant that is not known.  The density of  14:     Where w_i is a constant that is not known.  The density of 
 41:         reduced_energy:  E/T_i  a 2d array of shape( nreps, nbins) giving the 41:         reduced_energy:  E/T_i  a 2d array of shape( nreps, nbins) giving the
 42:             reduced energy of each bin 42:             reduced energy of each bin
 43:  43: 
 44:         note: this works perfectly well for 2d histograms as well.  In this case the 2d  44:         note: this works perfectly well for 2d histograms as well.  In this case the 2d 
 45:             histograms should be linearized 45:             histograms should be linearized
 46:          46:         
 47:         """ 47:         """
 48:         self.nreps, self.nbins = P.shape 48:         self.nreps, self.nbins = P.shape
 49:         assert P.shape == reduced_energy.shape 49:         assert P.shape == reduced_energy.shape
 50:         self.P = P 50:         self.P = P
  51:         #self.reduced_energy = reduced_energy
 51:          52:         
 52:         if np.any(self.P < 0): 53:         if np.any(self.P < 0):
 53:             raise ValueError("P has negative values") 54:             raise ValueError("P has negative values")
 54:          55:         
 55:         SMALL = 0. # this number is irrelevant, as long as it's not NaN 56:         SMALL = 1e-100
 56:         self.log_n_rE = np.where(self.P==0, SMALL,  57:         self.logP = np.where(self.P==0, SMALL, np.log(self.P))
 57:                                  np.log(self.P) + reduced_energy) 58:         self.n_rE = self.logP + reduced_energy
 58:          59:         
  60:         if False: #see how much effort is wasted
  61:             nallzero = len( np.where(np.abs(self.logP.sum(1)) < 1e-10 )[0])
  62:             print "WHAM: number of degrees of freedom", self.nreps + self.nbins
  63:             print "WHAM: number of irrelevant d.o.f. ", nallzero
  64: 
 59:     def getEnergy(self, X): 65:     def getEnergy(self, X):
 60:         """ 66:         """
 61:         X: is the array of unknowns of length nrep + nbins 67:         X: is the array of unknowns of length nrep + nbins
 62:             X[0:nrep] = {w_i}         : the replica unknowns 68:             X[0:nrep] = {w_i}         : the replica unknowns
 63:             X[nrep:] = {log(n_F(E))}  : the bin unknowns 69:             X[nrep:] = {log(n_F(E))}  : the bin unknowns
 64:  70: 
 65:         R(E,T_i) = log(n_F(E)) - log(w_i) - log( P(E,T_i)*exp(E/T_i) ) 71:         R(E,T_i) = log(n_F(E)) - log(w_i) - log( P(E,T_i)*exp(E/T_i) )
 66:      72:     
 67:         energy = sum_E sum_i P(E,T_i)*|R(E,T_i)|^2 73:         energy = sum_E sum_i P(E,T_i)*|R(E,T_i)|^2
 68:         """ 74:         """
 69:         wi = X[:self.nreps] 75:         wi = X[:self.nreps]
 70:         lognF = X[self.nreps:] 76:         lognF = X[self.nreps:]
 71:         energy = np.sum( self.P * (lognF[np.newaxis,:] - wi[:,np.newaxis] - self.log_n_rE)**2 ) 77:         """
  78:         energy = 0.
  79:         for irep in range(self.nreps):
  80:             for ibin in range(self.nbins):
  81:                 R = lognF[ibin] - wi[irep] -( logP[irep, ibin]  + reduced_energy[irep, ibin])
  82:                 energy += logP[irep, ibin] * R**2
  83:         """
  84:         energy = np.sum( self.P * (lognF[np.newaxis,:] - wi[:,np.newaxis] - self.n_rE)**2 )
 72:         return energy 85:         return energy
 73:  86: 
 74:     def getEnergyGradient(self, X): 87:     def getEnergyGradient(self, X):
 75:         """ 88:         """
 76:         X: is the array of unknowns of length nrep + nbins 89:         X: is the array of unknowns of length nrep + nbins
 77:             X[0:nrep] = {w_i}         : the replica unknowns 90:             X[0:nrep] = {w_i}         : the replica unknowns
 78:             X[nrep:] = {log(n_F(E))}  : the bin unknowns 91:             X[nrep:] = {log(n_F(E))}  : the bin unknowns
 79:  92: 
 80:         R(E,T_i) = log(n_F(E)) - log(w_i) - log( P(E,T_i)*exp(E/T_i) ) 93:         R(E,T_i) = log(n_F(E)) - log(w_i) - log( P(E,T_i)*exp(E/T_i) )
 81:      94:     
 82:         energy = sum_E sum_i P(E,T_i)*|R(E,T_i)|^2 95:         energy = sum_E sum_i P(E,T_i)*|R(E,T_i)|^2
 83:         """ 96:         """
 84:         wi = X[:self.nreps] 97:         wi = X[:self.nreps]
 85:         lognF = X[self.nreps:] 98:         lognF = X[self.nreps:]
 86:         R = lognF[np.newaxis,:] - wi[:,np.newaxis] - self.log_n_rE 99:         R = lognF[np.newaxis,:] - wi[:,np.newaxis] - self.n_rE
 87:         100:         
 88:         energy = np.sum( self.P * (R)**2 )101:         energy = np.sum( self.P * (R)**2 )
 89:         102:         
 90:         gradient = np.zeros(len(X))103:         gradient = np.zeros(len(X))
 91:         gradient[:self.nreps] = -2. * np.sum( self.P * R, axis=1 )104:         gradient[:self.nreps] = -2. * np.sum( self.P * R, axis=1 )
 92:         gradient[self.nreps:] = 2. * np.sum( self.P * R, axis=0 )105:         gradient[self.nreps:] = 2. * np.sum( self.P * R, axis=0 )
 93:         #print np.shape(gradient)106:         #print np.shape(gradient)
 94:         #print gradient107:         #print gradient
 95:     108:     
 96:         return energy, gradient109:         return energy, gradient


r27871/wham_utils.py 2016-07-06 15:37:39.144895535 +0100 r27870/wham_utils.py 2016-07-06 15:37:42.192936994 +0100
  1: import numpy as np  1: import numpy as np
  2: from numpy import log, exp  2: from numpy import log, exp
  3: from scipy.misc import logsumexp 
  4: from scipy.optimize import fmin_l_bfgs_b 
  5: try: 
  6:     from scipy.optimize import Result 
  7: except ImportError: 
  8:     # they changed the name from Result to OptimizeResult at some point 
  9:     from scipy.optimize import OptimizeResult as Result 
 10:   3: 
 11: def dos_from_offsets(visits, log_dos_all, offsets, nodata_value=0.): 
 12:     log_dos_all = log_dos_all + offsets[:,np.newaxis] 
 13:      
 14:     ldos = np.sum(log_dos_all * visits, axis=0) 
 15:     norm = visits.sum(0) 
 16:     ldos = np.where(norm > 0, ldos / norm, nodata_value) 
 17:     return ldos 
 18:   4: 
 19: def estimate_dos(visits, reduced_energy):  5: 
 20:     """estimate the density of states from the histograms of bins  6: def logSum1(a, b):
 21:      
 22:     Notes 
 23:     ----- 
 24:     The density of states is proportional to visits * exp(E/T).  Therefore 
 25:     the log density of states for each replica is related to one another by an 
 26:     additive constant.  This function will find those offsets and produce a 
 27:     guess for the total density of states. 
 28:      
 29:     """  7:     """
 30:     SMALL = 0.  8:     return log( exp(a) + exp(b) )
 31:     if len(visits.shape) != 2:  9:     """
 32:         raise ValueError("visits must have two dimensions") 10:     if a > b:
 33:     if visits.shape != reduced_energy.shape: 11:         return a + log(1.0 + exp(-a + b) )
 34:         raise ValueError("visits has a different shape from reduced_energy") 12:     else:
 35:     log_dos = np.where(visits==0, SMALL,  13:         return b + log(1.0 + exp(a - b) )
 36:                        np.log(visits) + reduced_energy)  
 37:      
 38:     offsets = [0.] 
 39:     nreps = visits.shape[0] 
 40:     for i in xrange(1,nreps): 
 41:         # find the difference in the log density of states 
 42:         ldos_diff = log_dos[i-1,:] - log_dos[i,:] 
 43:         # weight the difference by the minimum visits in each bin 
 44:         weights = np.where(visits[i-1,:] < visits[i,:], visits[i-1,:], visits[i,:]) 
 45:         new_offset = np.average(ldos_diff, weights=weights) 
 46:         offsets.append( offsets[-1] + new_offset) 
 47:     offsets = np.array(offsets) 
 48:     ldos = dos_from_offsets(visits, log_dos, offsets) 
 49:  14: 
 50:     return offsets, ldos 
 51:  15: 
  16: def logSum(log_terms):
  17:     """
  18:     Compute the log of a sum of terms whose logarithms are provided.
 52:  19: 
 53: def calc_Cv(Tlist, binenergy, logn_E, ndof, have_data=None, k_B=1.): 20:     REQUIRED ARGUMENTS  
 54:     """compute the heat capacity and other data from the density of states 21:       log_terms is the array (possibly multidimensional) containing the logs of the terms to be summed.
 55:      22: 
 56:     Parameters 23:     RETURN VALUES
 57:     ---------- 24:       log_sum is the log of the sum of the terms.
 58:     Tlist : list 
 59:         list of temperatures at which to do the computations 
 60:     binenergy : array 
 61:         The lower edge of the bins for the density of states data 
 62:     logn_E : numpy array 
 63:         the binned density of states. 
 64:     ndof : int 
 65:         number of degrees of freedom 
 66:     have_data : numpy array 
 67:         array of indices where data is available.  If have_data is None  
 68:         then all data is assumed good 
 69:     k_B : float 
 70:         Boltzmann's constant.  Conversion factor between temperature and energy 
 71:      
 72:     """ 25:     """
 73:     nz = have_data 26:     try:
  27:         return logSumFast(log_terms)
  28:     except ImportError:
  29:         #print "using slow logsum, install scipy weave inline blitz to speed up"
  30:         return logSumSlow(log_terms)
  31: 
  32: def logSumSlow(log_terms):
  33:     log_sum = log_terms[0] 
  34:     for lt in log_terms[1:]:
  35:         if log_sum > lt:
  36:             log_sum = log_sum + log(1.0 + exp(-log_sum + lt) )
  37:         else:
  38:             log_sum = lt + log(1.0 + exp(log_sum - lt) )
  39:     # return the log sum
  40:     return log_sum
  41: 
  42: #def logSumFort(log_terms):
  43: #    import _wham_utils
  44: #    return _wham_utils.logsum(log_terms)
  45: 
  46: #def logSumFast(log_terms):
  47: #    """
  48: #    there are two options available in case the user doesn't have package scipy (as
  49: #    is the case on some clusters) or doesn't have f2py
  50: #    """
  51: #    try:
  52: #        return logSumFort(log_terms)
  53: #    except ImportError:
  54: #        return logSumWeave(log_terms)
  55: 
  56: def logSumFast(log_terms):
  57:     lmax = np.max(log_terms)
  58: #    lsub = log_terms - lmax
  59:     result = np.log(np.sum(np.exp(log_terms - lmax))) + lmax
  60:     return result
  61:     
  62: 
  63: def calc_Cv(logn_E, visits1d, binenergy, NDOF, Treplica, k_B, TRANGE=None, NTEMP=100, use_log_sum = None):
  64:     if use_log_sum == None:
  65:         try:
  66:             logSumFast( np.array([.1, .2, .3]) )
  67:             use_log_sum = True
  68:             #print "using logsum"
  69:         except ImportError:
  70:             #dont use log sum unless the fast logsum is working
  71:             #print "not using logsum because it's too slow.  Install scipy weave to use the fast version of logsum"
  72:             use_log_sum = False
  73:     else:
  74:         #print "use_log_sum = ", use_log_sum
  75:         pass
 74:      76:     
 75:     assert(not np.any(np.isnan(logn_E[nz]))) 
 76:     assert(not np.any(np.isnan(binenergy[nz]))) 
 77:      77:     
 78:     Ebin_min = binenergy[nz].min() - 1. 
 79:     nebins = len(binenergy) 
 80:  78: 
 81:     # now calculate partition functions, energy expectation values, and Cv 79:     #put some variables in this namespace
 82:     dataout = np.zeros([len(Tlist), 6]) 80:     nrep, nebins = np.shape(visits1d)
  81:     #print "nreps, nebins", nrep, nebins
  82: 
  83:     nz = np.where( visits1d.sum(0) != 0)[0]
  84: 
  85:     #find the ocupied bin with the minimum energy
  86:     EREF = np.min(binenergy[nz]) - 1.
  87: 
  88:     #now calculate partition functions, energy expectation values, and Cv
  89:     if TRANGE is None:
  90:         #NTEMP = 100 # number of temperatures to calculate expectation values
  91:         TMAX = Treplica[-1]
  92:         TMIN = Treplica[0]
  93:         TRANGE = np.linspace(TMIN, TMAX, NTEMP)
  94: #        TINT=(TMAX-TMIN)/(NTEMP-1)
  95: #        TRANGE = [ TMIN + i*TINT for i in range(NTEMP) ]
  96:     NTEMP = len(TRANGE)
  97: 
  98:     dataout = np.zeros( [NTEMP, 6] )
 83:     if abs((binenergy[-1] - binenergy[-2]) - (binenergy[-2] - binenergy[-3]) ) > 1e-7: 99:     if abs((binenergy[-1] - binenergy[-2]) - (binenergy[-2] - binenergy[-3]) ) > 1e-7:
 84:         print "calc_Cv: WARNING: dE is not treated correctly for uneven energy bins"100:         print "calc_Cv: WARNING: dE is not treated correctly for exponential energy bins"
 85: 101: 
 86:     for count, T in enumerate(Tlist):102:     for count,T in enumerate(TRANGE):
 87:         kBT = k_B * T103:         kBT = k_B*T
 88:         # find expoffset so the exponentials don't blow up104:         #find expoffset so the exponentials don't blow up
 89:         dummy = logn_E[nz] - binenergy[nz] / kBT105:         dummy = logn_E[nz] - (binenergy[nz] - EREF)/kBT
 90: 106:         expoffset = np.max(dummy)
 91:         lZ0 = logsumexp(dummy)107:         if not use_log_sum:
 92:         lZ1 = logsumexp(dummy + log(binenergy[nz] - Ebin_min))108:             dummy = np.exp(dummy)
 93:         lZ2 = logsumexp(dummy + 2. * log(binenergy[nz] - Ebin_min))109:             Z0 = np.sum(dummy)
 110:             Z1 = np.sum( dummy *  (binenergy[nz] - EREF) )
 111:             Z2 = np.sum( dummy *  (binenergy[nz] - EREF)**2 )
 112:             lZ0 = np.log(Z0)
 113:             lZ1 = np.log(Z1)
 114:             lZ2 = np.log(Z2)
 115:         else:
 116:             lZ0 = logSum( dummy )
 117:             lZ1 = logSum( dummy + log(binenergy[nz] - EREF) )
 118:             lZ2 = logSum( dummy + 2.*log(binenergy[nz] - EREF) )
 119:             Z0 = np.exp( lZ0 )
 120:             Z1 = np.exp( lZ1 )
 121:             Z2 = np.exp( lZ2 )
 94: 122: 
 95:         i = nebins-1123:         i = nebins-1
 96:         #if abs(binenergy[-1] - binenergy[-2]) > 1e-7:124:         #if abs(binenergy[-1] - binenergy[-2]) > 1e-7:
 97:         #    print "calc_Cv: WARNING: dE is not treated correctly for exponential bins"125:         #    print "calc_Cv: WARNING: dE is not treated correctly for exponential bins"
 98:         if i == nebins-1:126:         if i == nebins-1:
 99:             dE = binenergy[i] - binenergy[i-1]127:             dE = binenergy[i]-binenergy[i-1]
100:         else:128:         else:
101:             dE = binenergy[i+1] - binenergy[i]129:             dE = binenergy[i+1]-binenergy[i]
102:             130:             
103:         if dE/kBT < 1e-7:131:         if dE/kBT < 1.0E-7:
104:             OneMExp = -dE/kBT132:             ONEMEXP=-dE/kBT
105:         else:133:         else:
106:             OneMExp = 1.0 - np.exp(dE/kBT)134:             ONEMEXP= 1.0-np.exp(dE/kBT)
107: 135: 
108:         Eavg = ndof * kBT / 2.0 + kBT + dE / OneMExp + exp(lZ1-lZ0) + Ebin_min136:         Eavg = NDOF*kBT/2.0 + 1.0*(kBT + dE/ONEMEXP) + exp(lZ1-lZ0) + EREF
109:         137:         
110:         # correction to the Cv for finite size of bins138:         Cv = NDOF/2. + 1.0*(1.0 - dE**2 * exp(dE/kBT)/(ONEMEXP**2*kBT**2)) \
111:         bin_correction = 1.0 - dE**2 * exp(dE / kBT) / (OneMExp**2 * kBT**2)139:             - exp(lZ1 - lZ0)**2 / kBT**2 + exp(lZ2-lZ0) / kBT**2
112:         Cv = k_B * (ndof / 2. + (exp(lZ2-lZ0) - exp(lZ1 - lZ0)**2) / kBT**2140:             #- (Z1/(Z0*kBT))**2 + Z2/(Z0*kBT**2)
113:                     + bin_correction 
114:                     ) 
115:         141:         
116:         dataout[count,0] = T142:         dataout[count,0] = T
117:         dataout[count,1] = lZ0143:         dataout[count,1] = lZ0 + expoffset
118:         dataout[count,2] = lZ1144:         dataout[count,2] = lZ1 + expoffset
119:         dataout[count,3] = lZ2145:         dataout[count,3] = lZ2 + expoffset
120:         dataout[count,4] = Eavg146:         dataout[count,4] = Eavg
121:         dataout[count,5] = Cv147:         dataout[count,5] = Cv
 148: 
 149:         
 150:         #np.array([kBT, Z0*np.exp(expoffset), Z1*np.exp(expoffset), Z2*np.exp(expoffset), Eavg, Cv, np.log(Z0)+expoffset, np.log(Z1)+expoffset, np.log(Z2)+expoffset]).tofile(fout," ")
122:         151:         
 152:         #fout.write("\n")
 153: 
123:     return dataout154:     return dataout
124: 155: 
125: 156: 
126: def lbfgs_scipy(coords, pot, iprint=-1, tol=1e-3, nsteps=10000):157: def lbfgs_scipy(coords, pot, iprint=-1, tol=1e-3, nsteps=1500):
127:     """158:     """
128:     a wrapper function for lbfgs routine in scipy159:     a wrapper function for lbfgs routine in scipy
129:     160:     
130:     .. warn::161:     .. warn::
131:         the scipy version of lbfgs uses linesearch based only on energy162:         the scipy version of lbfgs uses linesearch based only on energy
132:         which can make the minimization stop early.  When the step size163:         which can make the minimization stop early.  When the step size
133:         is so small that the energy doesn't change to within machine precision (times the164:         is so small that the energy doesn't change to within machine precision (times the
134:         parameter `factr`) the routine declares success and stops.  This sounds fine, but165:         parameter `factr`) the routine declares success and stops.  This sounds fine, but
135:         if the gradient is analytical the gradient can still be not converged.  This is166:         if the gradient is analytical the gradient can still be not converged.  This is
136:         because in the vicinity of the minimum the gradient changes much more rapidly then167:         because in the vicinity of the minimum the gradient changes much more rapidly then
137:         the energy.  Thus we want to make factr as small as possible.  Unfortunately,168:         the energy.  Thus we want to make factr as small as possible.  Unfortunately,
138:         if we make it too small the routine realizes that the linesearch routine169:         if we make it too small the routine realizes that the linesearch routine
139:         isn't working and declares failure and exits.170:         isn't working and declares failure and exits.
140:         171:         
141:         So long story short, if your tolerance is very small (< 1e-6) this routine172:         So long story short, if your tolerance is very small (< 1e-6) this routine
142:         will probably stop before truly reaching that tolerance.  If you reduce `factr` 173:         will probably stop before truly reaching that tolerance.  If you reduce `factr` 
143:         too much to mitigate this lbfgs will stop anyway, but declare failure misleadingly.  174:         too much to mitigate this lbfgs will stop anyway, but declare failure misleadingly.  
144:     """175:     """
145:     assert hasattr(pot, "getEnergyGradient")176:     assert hasattr(pot, "getEnergyGradient")
 177:     from scipy.optimize import Result, fmin_l_bfgs_b
146:     res = Result()178:     res = Result()
147:     res.coords, res.energy, dictionary = fmin_l_bfgs_b(pot.getEnergyGradient, 179:     res.coords, res.energy, dictionary = fmin_l_bfgs_b(pot.getEnergyGradient, 
148:             coords, iprint=iprint, pgtol=tol, maxfun=nsteps, factr=10.)180:             coords, iprint=iprint, pgtol=tol, maxfun=nsteps, factr=10.)
149:     res.grad = dictionary["grad"]181:     res.grad = dictionary["grad"]
150:     res.nfev = dictionary["funcalls"]182:     res.nfev = dictionary["funcalls"]
151:     warnflag = dictionary['warnflag']183:     warnflag = dictionary['warnflag']
152:     #res.nsteps = dictionary['nit'] #  new in scipy version 0.12184:     #res.nsteps = dictionary['nit'] #  new in scipy version 0.12
153:     res.nsteps = res.nfev185:     res.nsteps = res.nfev
154:     res.message = dictionary['task']186:     res.message = dictionary['task']
155:     res.success = True187:     res.success = True
156:     if warnflag > 0:188:     if warnflag > 0:
157:         print "warning: problem with quench: ",189:         print "warning: problem with quench: ",
158:         res.success = False190:         res.success = False
159:         if warnflag == 1:191:         if warnflag == 1:
160:             res.message = "too many function evaluations"192:             res.message = "too many function evaluations"
161:         else:193:         else:
162:             res.message = str(dictionary['task'])194:             res.message = str(dictionary['task'])
163:         print res.message195:         print res.message
164:         print "    the energy is", res.energy, "the rms gradient is", np.linalg.norm(res.grad) / np.sqrt(res.grad.size), "nfev", res.nfev196:         print "    the energy is", res.energy, "the rms gradient is", np.linalg.norm(res.grad) / np.sqrt(res.grad.size), "nfev", res.nfev
165:         print "    X: ", res.coords197:         print res.coords
166:     #note: if the linesearch fails the lbfgs may fail without setting warnflag.  Check198:     #note: if the linesearch fails the lbfgs may fail without setting warnflag.  Check
167:     #tolerance exactly199:     #tolerance exactly
168:     if False:200:     if False:
169:         if res.success:201:         if res.success:
170:             maxV = np.max( np.abs(res.grad) )202:             maxV = np.max( np.abs(res.grad) )
171:             if maxV > tol:203:             if maxV > tol:
172:                 print "warning: gradient seems too large", maxV, "tol =", tol, ". This is a known, but not understood issue of scipy_lbfgs"204:                 print "warning: gradient seems too large", maxV, "tol =", tol, ". This is a known, but not understood issue of scipy_lbfgs"
173:                 print res.message205:                 print res.message
174:     res.rms = res.grad.std()206:     res.rms = res.grad.std()
175:     return res207:     return res
176: 208: 
 209: 
 210: if __name__ == "__main__":
 211:     pass
177: 212: 


legend
Lines Added 
Lines changed
 Lines Removed

hdiff - version: 2.1.0