Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.
# for i in range(numboots):
# boot = np.random.randint(L,size=L)
# Rhoboot.append(calc_rho(None,**{'r_':Rhos[boot]}))
# Rhoboot = np.array(Rhoboot)
# Rho_err = np.std(Rhoboot)
if FDCheck:
Sep = printcool("Numerical Derivative:")
GRho1 = property_derivatives(mvals, h, pdb, FF, Xyzs, Settings, Sim, kT, calc_rho, {'r_':Rhos}, Boxes)
FF.print_map(vals=GRho1)
Sep = printcool("Difference (Absolute, Fractional):")
absfrac = ["% .4e % .4e" % (i-j, (i-j)/j) for i,j in zip(GRho, GRho1)]
FF.print_map(vals=absfrac)
print "Box energy:", np.mean(Energies)
print "Monomer energy:", np.mean(mEnergies)
Sep = printcool("Enthalpy of Vaporization: % .4f +- %.4f kJ/mol, Derivatives below" % (Hvap_avg, Hvap_err))
FF.print_map(vals=GHvap)
print Sep
# Define some things to make the analytic derivatives easier.
Gbar = np.mean(G,axis=1)
def deprod(vec):
return flat(np.mat(G)*col(vec))/N
def covde(vec):
return flat(np.mat(G)*col(vec))/N - Gbar*np.mean(vec)
def avg(vec):
return np.mean(vec)
## Thermal expansion coefficient and bootstrap error estimation
def calc_alpha(b = None, **kwargs):
if b == None: b = np.ones(L,dtype=float)
if 'h_' in kwargs:
# When all initial parameter values are zero, it could be a problem...
if maxval == 0:
maxval += 1
rsfactors[termtype] = maxval
rsfac_list.append(termtype)
# Physically motivated overrides
rs_override(rsfactors,termtype)
# Overrides from input file
for termtype in self.priors:
rsfac_list.append(termtype)
rsfactors[termtype] = self.priors[termtype]
# for line in os.popen("awk '/rsfactor/ {print $2,$3}' %s" % pkg.options).readlines():
# rsfactors[line.split()[0]] = float(line.split()[1])
if printfacs:
bar = printcool("Rescaling Factors (Lower Takes Precedence):",color=1)
logger.info(''.join([" %-35s : %.5e\n" % (i, rsfactors[i]) for i in rsfac_list]))
logger.info(bar)
## The array of rescaling factors
self.rs = ones(len(self.pvals0))
for pnum in range(len(self.pvals0)):
for termtype in rsfac_list:
if termtype in self.plist[pnum]:
self.rs[pnum] = rsfactors[termtype]
if X > (X_prev + self.err_tol):
Best_Step = 0
# Toggle switch for rejection (experimenting with no rejection)
Rejects = True
GOODSTEP = 0
Reevaluate = True
trust = max(ndx*(1./(1+a)), self.mintrust)
logger.info("Rejecting step and reducing trust radius to % .4e\n" % trust)
if Rejects:
xk = xk_prev.copy()
if Reevaluate:
restep = True
color = "\x1b[91m"
logger.info("%6s%12s%12s%12s%14s%12s%12s\n" % ("Step", " |k| "," |dk| "," |grad| "," -=X2=- ","Delta(X2)", "StepQual"))
logger.info("%6i%12.3e%12.3e%12.3e%s%14.5e\x1b[0m%12.3e% 11.3f\n\n" % (ITERATION_NUMBER, nxk, ndx, ngr, color, X, stdfront, Quality))
printcool("Objective function rises!\nRe-evaluating at the previous point..",color=1)
ITERATION_NUMBER += 1
data = self.Objective.Full(xk,Ord,verbose=True)
GOODSTEP = 1
X, G, H = data['X'], data['G'], data['H']
X_prev = X
dx *= 0
ndx = norm(dx)
nxk = norm(xk)
ngr = norm(G)
Quality = 0.0
color = "\x1b[0m"
else:
color = "\x1b[91m"
G = G_prev.copy()
H = H_stor.copy()
data = deepcopy(datastor)
logger.info(str(vals_in) + '\n')
logger.info('These parameter values will be used:\n')
logger.info(str(scanvals) + '\n')
minvals = None
minobj = 1e100
for pidx in idx:
if MathPhys:
logger.info("Scanning parameter %i (%s) in the mathematical space\n" % (pidx,self.FF.plist[pidx]))
vals = self.mvals0.copy()
else:
logger.info("Scanning parameter %i (%s) in the physical space\n" % (pidx,self.FF.plist[pidx]))
self.FF.use_pvals = True
vals = self.FF.pvals0.copy()
counter = 1
for i in scanvals:
printcool("Parameter %i (%s) Value is now % .4e ; Step %i/%i" % (pidx, self.FF.plist[pidx],i,counter,len(scanvals)), color=1,sym="@")
vals[pidx] = i
data = self.Objective.Full(vals,Order=0,verbose=True)
if data['X'] < minobj:
minobj = data['X']
minvals = vals.copy()
logger.info("Value = % .4e Objective = % .4e\n" % (i, data['X']))
ITERATION += 1
counter += 1
return minvals
xmin, Jmin, T, feval, iters, accept, status = optimize.anneal(xwrap(self.Objective.Full), self.mvals0, lower=self.mvals0-1*self.trust0*np.ones(self.np),
upper=self.mvals0+self.trust0*np.ones(self.np),schedule='boltzmann', full_output=True)
scodes = {0 : "Points no longer changing.",
1 : "Cooled to final temperature.",
2 : "Maximum function evaluations reached.",
3 : "Maximum cooling iterations reached.",
4 : "Maximum accepted query locations reached.",
5 : "Final point not the minimum amongst encountered points."}
logger.info("Simulated annealing info:\n")
logger.info("Status: %s \n" % scodes[status])
logger.info("Function evaluations: %i" % feval)
logger.info("Cooling iterations: %i" % iters)
logger.info("Tests accepted: %i" % iters)
return xmin
elif Algorithm == "basinhopping":
printcool("Minimizing Objective Function using Basin Hopping Method" , ansi=1, bold=1)
T = xwrap(self.Objective.Full)(self.mvals0)
Result = optimize.basinhopping(xwrap(self.Objective.Full), self.mvals0, niter=self.maxstep, T=T, stepsize=self.trust0, interval=20,
minimizer_kwargs={'method':'nelder-mead','options':{'xtol': self.convergence_step,'ftol':self.convergence_objective}}, disp=True)
logger.info(Result.message + "\n")
return Result.x
elif Algorithm == "cg":
printcool("Minimizing Objective Function using\nPolak-Ribiere Conjugate Gradient Method" , ansi=1, bold=1)
return optimize.fmin_cg(xwrap(self.Objective.Full,callback=False),self.mvals0,fprime=gwrap(self.Objective.Full),gtol=self.convergence_gradient)
elif Algorithm == "tnc":
printcool("Minimizing Objective Function using\nTruncated Newton Algorithm (Unconfirmed)" , ansi=1, bold=1)
Result = optimize.fmin_tnc(fgwrap(self.Objective.Full,callback=False),self.mvals0,
maxfun=self.maxstep,ftol=self.convergence_objective,pgtol=self.convergence_gradient,xtol=self.convergence_objective)
return Result.x
elif Algorithm == "ncg":
printcool("Minimizing Objective Function using\nNewton-CG Algorithm" , ansi=1, bold=1)
Result = optimize.fmin_ncg(xwrap(self.Objective.Full,callback=False),self.mvals0,fprime=gwrap(self.Objective.Full,callback=False),
color = "\x1b[92m"
my_gfunc.x_best = Objective
else:
color = "\x1b[91m"
if verbose:
if self.print_vals:
logger.info("k=" + ' '.join(["% .4f" % i for i in mvals]) + '\n')
logger.info("|Grad|= %12.3e X2= %s%12.3e\x1b[0m d(X2)= %12.3e\n\n" % (norm(Answer),color,Objective,dx))
return Answer
my_gfunc.x_best = None
return my_gfunc
if Algorithm == "powell":
printcool("Minimizing Objective Function using Powell's Method" , ansi=1, bold=1)
return optimize.fmin_powell(xwrap(self.Objective.Full),self.mvals0,ftol=self.conv_obj,xtol=self.conv_stp,maxiter=self.maxstep)
elif Algorithm == "simplex":
printcool("Minimizing Objective Function using Simplex Method" , ansi=1, bold=1)
return optimize.fmin(xwrap(self.Objective.Full),self.mvals0,ftol=self.conv_obj,xtol=self.conv_stp,maxiter=self.maxstep,maxfun=self.maxstep*10)
elif Algorithm == "anneal":
printcool("Minimizing Objective Function using Simulated Annealing" , ansi=1, bold=1)
return optimize.anneal(xwrap(self.Objective.Full),self.mvals0,lower=-1*self.trust0*np.ones(self.np),upper=self.trust0*np.ones(self.np),schedule='boltzmann')
elif Algorithm == "cg":
printcool("Minimizing Objective Function using Conjugate Gradient" , ansi=1, bold=1)
return optimize.fmin_cg(xwrap(self.Objective.Full,verbose=False),self.mvals0,fprime=gwrap(self.Objective.Full),gtol=self.conv_grd)
# logger.info("The force field has been written to the '%s' directory.\n" % self.resdir)
outfnm = self.save_mvals_to_input(xk)
# logger.info("Input file with optimization parameters saved to %s.\n" % outfnm)
printcool("The force field has been written to the %s directory.\n"
"Input file with optimization parameters saved to %s." % (self.resdir, outfnm), color=0)
# "To reload these parameters, use %s as the input\n"
# "file without changing the '%s' directory." %
# (outfnm, outfnm, self.FF.ffdir), color=0, center=False, sym2='-')
## Write out stuff to checkpoint file
self.writechk()
## Print out final message
logger.info("Wall time since calculation start: %.1f seconds\n" % (time.time() - t0))
if self.failmsg:
bar = printcool("It is possible to commit no errors and still lose.\nThat is not a weakness. That is life.",ansi="40;93")
else:
bar = printcool("Calculation Finished.\n---==( May the Force be with you! )==---",ansi="1;44;93")
return xk
self.Objective.ObjDict_Last[key] = val
restep = False
dx, dX_expect, bump = self.step(xk, data, trust)
old_pk = self.FF.create_pvals(xk)
old_xk = xk.copy()
# Increment the iteration counter.
ITERATION_NUMBER += 1
# Take a step in the parameter space.
xk += dx
if self.print_vals:
pk = self.FF.create_pvals(xk)
dp = pk - old_pk
bar = printcool("Mathematical Parameters (Current + Step = Next)",color=5)
self.FF.print_map(vals=["% .4e %s %.4e = % .4e" % (old_xk[i], '+' if dx[i] >= 0 else '-', abs(dx[i]), xk[i]) for i in range(len(xk))])
logger.info(bar)
bar = printcool("Physical Parameters (Current + Step = Next)",color=5)
self.FF.print_map(vals=["% .4e %s %.4e = % .4e" % (old_pk[i], '+' if dp[i] >= 0 else '-', abs(dp[i]), pk[i]) for i in range(len(pk))])
logger.info(bar)
# Evaluate the objective function and its derivatives.
data = self.Objective.Full(xk,Ord,verbose=True)
X, G, H = data['X'], data['G'], data['H']
ndx = norm(dx)
nxk = norm(xk)
ngr = norm(G)
drc = abs(flat(dx)).argmax()
dX_actual = X - X_prev
try:
Quality = dX_actual / dX_expect
except:
logger.warning("Warning: Step size of zero detected (i.e. wrong direction). Try reducing the finite_difference_h parameter\n")
Quality = 1.0 # This is a step length of zero.
# Increment the parameters.
xk += dx
ndx = np.linalg.norm(dx)
#================================#
#| Increase iteration number. |#
#================================#
ITERATION += 1
self.iteration += 1
# The search code benefits from knowing the step size here.
if self.trust0 < 0:
trust = ndx
# Print parameter values.
if self.print_vals:
pk = self.FF.create_pvals(xk)
dp = pk - pk_prev
bar = printcool("Mathematical Parameters (Current + Step = Next)",color=5)
self.FF.print_map(vals=["% .4e %s %.4e = % .4e" % (xk_prev[i], '+' if dx[i] >= 0 else '-', abs(dx[i]), xk[i]) for i in range(len(xk))])
logger.info(bar)
bar = printcool("Physical Parameters (Current + Step = Next)",color=5)
self.FF.print_map(vals=["% .4e %s %.4e = % .4e" % (pk_prev[i], '+' if dp[i] >= 0 else '-', abs(dp[i]), pk[i]) for i in range(len(pk))])
logger.info(bar)
# Write checkpoint file.
self.chk = {'xk': xk, 'X' : X, 'G' : G, 'H': H, 'X_hist': X_hist, 'trust': trust}
if self.wchk_step:
self.writechk()
outfnm = self.save_mvals_to_input(xk)
logger.info("Input file with saved parameters: %s\n" % outfnm)
# Check for whether the maximum number of optimization cycles is reached.
if ITERATION == self.maxstep:
logger.info("Maximum number of optimization steps reached (%i)\n" % ITERATION)
break
# Check for whether the step size is too small to continue.
For each element in the Hessian, use a five-point stencil in both
parameter indices to compute a finite-difference derivative, and
compare it to the analytic result.
This is meant to be a foolproof checker, so it is pretty slow. We
could write a faster checker if we assumed we had accurate first
derivatives, but it's better to not make that assumption.
The second derivative is computed by double-wrapping the objective
function via the 'wrap2' function.
"""
Adata = self.Objective.Full(self.mvals0,Order=2)['H']
Fdata = np.zeros((self.FF.np,self.FF.np),dtype=float)
printcool("Checking second derivatives by finite difference!\n%-8s%-20s%-20s%13s%13s%13s%13s" \
% ("Index", "Parameter1 ID", "Parameter2 ID", "Analytic","Numerical","Difference","Fractional"),bold=1,color=5)
# Whee, our double-wrapped finite difference second derivative!
def wrap2(mvals0,pidxi,pidxj):
def func1(arg):
mvals = list(mvals0)
mvals[pidxj] += arg
return f1d5p(fdwrap(self.Objective.Full,mvals,pidxi,'X',Order=0),self.h)
return func1
for i in range(self.FF.np):
for j in range(i,self.FF.np):
Fdata[i,j] = f1d5p(wrap2(self.mvals0,i,j),self.h)
Denom = max(abs(Adata[i,j]),abs(Fdata[i,j]))
Denom = Denom > 1e-8 and Denom or 1e-8
D = Adata[i,j] - Fdata[i,j]