removed merge debris [ci skip]

This commit is contained in:
edoapra 2020-07-01 11:12:46 -07:00
parent 471a9083ad
commit c9099c8623
No known key found for this signature in database
GPG key ID: 472C4EBB4F3AEDFA

View file

@ -152,17 +152,10 @@ def quadratic_step(trust, g0, h0):
if h0 > 0:
delta2 = -g0/h0
if abs(delta2) > trust:
<<<<<<< HEAD
print(" Step restriction: %f %f " % (delta2, trust))
delta2 = abs(trust*delta2)/delta2
else:
print(" Negative curvature ")
=======
print (" Step restriction: %12.5f %12.5f " % (delta2, trust))
delta2 = abs(trust*delta2)/delta2
else:
print (" Negative curvature ")
>>>>>>> b1a8eff3c66eb281ed1c2ac4f342fd72cd2a7e46
delta2 = -abs(trust*g0)/g0
return delta2
@ -176,42 +169,25 @@ def linesearch(func, x0, s, lsgrad, eps):
# bracketed the minimum or gone downhil with enough
# energy difference to start fitting
<<<<<<< HEAD
print(" Line search: step alpha grad hess value")
print(" ---- --------- -------- -------- ----------------")
=======
print (" Line search: step alpha grad hess value")
print (" ---- --------- -------- -------- ----------------")
>>>>>>> b1a8eff3c66eb281ed1c2ac4f342fd72cd2a7e46
trust = 0.2
alpha0 = 0.0
f0 = func(x0)
<<<<<<< HEAD
print(" %9.2e %8.1e %16.8f" % (alpha0, lsgrad, f0))
=======
print (" %9.2e %8.1e %16.8f" % (alpha0 , lsgrad , f0))
>>>>>>> b1a8eff3c66eb281ed1c2ac4f342fd72cd2a7e46
if lsgrad < 0:
alpha1 = alpha0 + trust
else:
alpha1 = alpha0 - trust
f1 = func(takestep(x0,s,alpha1))
<<<<<<< HEAD
print(" %9.2e %16.8f" % (alpha1, f1))
while f1 > f0:
if trust < 0.00125:
print(" system is too badly conditioned for initial step")
=======
print (" %9.2e %16.8f" % \
(alpha1, f1))
while f1 > f0:
if trust < 0.00125:
print (" system is too badly conditioned for initial step")
>>>>>>> b1a8eff3c66eb281ed1c2ac4f342fd72cd2a7e46
return (alpha0,f0) # Cannot seem to find my way
trust = trust * 0.5
if lsgrad < 0:
@ -219,12 +195,8 @@ def linesearch(func, x0, s, lsgrad, eps):
else:
alpha1 = alpha0 - trust
f1 = func(takestep(x0,s,alpha1))
<<<<<<< HEAD
print(" %9.2e %16.8f" % (alpha1, f1))
=======
print (" %9.2e %16.8f" % \
(alpha1, f1))
>>>>>>> b1a8eff3c66eb281ed1c2ac4f342fd72cd2a7e46
g0 = lsgrad
h0 = (f1-f0-alpha1*g0)/alpha1**2
@ -244,14 +216,6 @@ def linesearch(func, x0, s, lsgrad, eps):
if iter == 1:
f2prev = f2
<<<<<<< HEAD
print(" %9.2e %16.8f" % (alpha2, f2))
# Check for convergence or insufficient precision to proceed further
if (abs(f0-f1)<(10*eps)) and (abs(f1-f2)<(10*eps)):
print(" ")
print(" Insufficient precision ... terminating LS")
=======
print (" %9.2e %16.8f" % \
(alpha2, f2))
@ -259,20 +223,14 @@ def linesearch(func, x0, s, lsgrad, eps):
if (abs(f0-f1)<(10*eps)) and (abs(f1-f2)<(10*eps)):
print (" "),
print (" Insufficient precision ... terminating LS")
>>>>>>> b1a8eff3c66eb281ed1c2ac4f342fd72cd2a7e46
break
if (f2-f2prev) > 0:
# New point is higher than previous worst
if nbackstep < 3:
nbackstep = nbackstep + 1
<<<<<<< HEAD
print(" ")
print(" Back stepping due to uphill step")
=======
print (" "),
print (" Back stepping due to uphill step")
>>>>>>> b1a8eff3c66eb281ed1c2ac4f342fd72cd2a7e46
trust = max(0.01,0.2*abs(alpha2 - alpha0)) # Reduce trust radius
alpha2 = alpha0 + 0.2*(alpha2 - alpha0)
continue
@ -292,20 +250,12 @@ def linesearch(func, x0, s, lsgrad, eps):
alpha1, alpha2, f1, f2 = alpha2, alpha1, f2, f1
(f0, g0, h0) = quadfit(alpha0, f0, alpha1, f1, alpha2, f2)
<<<<<<< HEAD
print(" %4i %9.2e %8.1e %8.1e %16.8f" %(iter, alpha0, g0, h0, f0))
if (h0>0.0) and (abs(g0) < 0.03*abs(lsgrad)):
print(" ")
print(" gradient reduced 30-fold ... terminating LS")
=======
print (" %4i %9.2e %8.1e %8.1e %16.8f" % \
(iter, alpha0, g0, h0, f0))
if (h0>0.0) and (abs(g0) < 0.03*abs(lsgrad)):
print (" "),
print (" gradient reduced 30-fold ... terminating LS")
>>>>>>> b1a8eff3c66eb281ed1c2ac4f342fd72cd2a7e46
break
# Determine the next step
@ -313,13 +263,8 @@ def linesearch(func, x0, s, lsgrad, eps):
alpha2 = alpha0 + delta
df = g0*delta + 0.5*h0*delta*delta
if abs(df) < 10.0*eps:
<<<<<<< HEAD
print(" ")
print(" projected energy reduction < 10*eps ... terminating LS")
=======
print (" "),
print (" projected energy reduction < 10*eps ... terminating LS")
>>>>>>> b1a8eff3c66eb281ed1c2ac4f342fd72cd2a7e46
break
@ -458,11 +403,7 @@ def hessian_update_bfgs(hp, dx, g, gp):
for j in range(n):
h[i][j] = h[i][j] + dg[i]*dg[j]/dxdg - hdx[i]*hdx[j]/dxhdx
else:
<<<<<<< HEAD
print(' BFGS not updating dxdg (%e), dgdg (%e), dxhdx (%f), dxdx(%e)' % (dxdg, dgdg, dxhdx, dxdx))
=======
print (' BFGS not updating dxdg (%e), dgdg (%e), dxhdx (%f), dxdx(%e)' % (dxdg, dgdg, dxhdx, dxdx))
>>>>>>> b1a8eff3c66eb281ed1c2ac4f342fd72cd2a7e46
return h
@ -500,27 +441,16 @@ def quasinr(func, guess, tol, eps, printvar=None):
(value,g,h) = numderiv(func, x, step, eps)
gmax = max(map(abs,g))
<<<<<<< HEAD
print(' ')
print(' iter gmax value ')
print(' ---- --------- ----------------')
print("%4i %9.2e %16.8f" % (iter,gmax,value))
=======
print (' ')
print (' iter gmax value ')
print (' ---- --------- ----------------')
print ("%4i %9.2e %16.8f" % (iter,gmax,value))
>>>>>>> b1a8eff3c66eb281ed1c2ac4f342fd72cd2a7e46
if (printvar):
printvar(x)
if gmax < tol:
<<<<<<< HEAD
print(' Converged!')
=======
print (' Converged!')
>>>>>>> b1a8eff3c66eb281ed1c2ac4f342fd72cd2a7e46
break
if iter == 0:
@ -532,22 +462,14 @@ def quasinr(func, guess, tol, eps, printvar=None):
(v,e) = jacobi(hessian)
emax = max(map(abs,e))
emin = emax*1e-4 # Control noise in small eigenvalues
<<<<<<< HEAD
print('\n Eigenvalues of the Hessian:')
=======
print ('\n Eigenvalues of the Hessian:')
>>>>>>> b1a8eff3c66eb281ed1c2ac4f342fd72cd2a7e46
printvector(e)
# Transform to spectral form, take step, transform back
gs = mxv(v,g)
for i in range(n):
if e[i] < emin:
<<<<<<< HEAD
print(' Mode %d: small/negative eigenvalue (%f).' % (i, e[i]))
=======
print (' Mode %d: small/negative eigenvalue (%f).' % (i, e[i]))
>>>>>>> b1a8eff3c66eb281ed1c2ac4f342fd72cd2a7e46
s[i] = -gs[i]/emin
else:
s[i] = -gs[i]/e[i]
@ -558,12 +480,8 @@ def quasinr(func, guess, tol, eps, printvar=None):
for i in range(n):
trust = max(abs(x[i]),abs(x[i]/sqrt(max(1e-4,abs(hessian[i][i])))))
if abs(s[i]) > trust:
<<<<<<< HEAD
print(' restricting ', i, trust, abs(x[i]), abs(x[i]/sqrt(abs(hessian[i][i]))), s[i])
=======
print (' restricting ', i, trust, abs(x[i]), \
abs(x[i]/sqrt(abs(hessian[i][i]))), s[i])
>>>>>>> b1a8eff3c66eb281ed1c2ac4f342fd72cd2a7e46
scale = min(scale,trust/abs(s[i]))
if scale != 1.0:
for i in range(n):
@ -572,11 +490,7 @@ def quasinr(func, guess, tol, eps, printvar=None):
(alpha,value) = linesearch(func, x, s, dot(s,g), eps)
if alpha == 0.0:
<<<<<<< HEAD
print(' Insufficient precision to proceed further')
=======
print (' Insufficient precision to proceed further')
>>>>>>> b1a8eff3c66eb281ed1c2ac4f342fd72cd2a7e46
break
for i in range(n):
@ -601,15 +515,6 @@ def cgminold(func, dfunc, guess, tol):
g = dfunc(x)
gmax = max(map(abs,g))
<<<<<<< HEAD
print(' ')
print(' iter gmax value ')
print(' ---- --------- ----------------')
print("%4i %9.2e %16.8f" % (iter,gmax,value))
if gmax < tol:
print(' Converged!')
=======
print (' ')
print (' iter gmax value ')
print (' ---- --------- ----------------')
@ -617,7 +522,6 @@ def cgminold(func, dfunc, guess, tol):
if gmax < tol:
print (' Converged!')
>>>>>>> b1a8eff3c66eb281ed1c2ac4f342fd72cd2a7e46
break
if (iter == 0) or ((iter%20) == 0):
@ -658,15 +562,6 @@ def cgmin(func, dfunc, guess, tol, precond=None, reset=None):
g = dfunc(x)
gmax = max(map(abs,g))
<<<<<<< HEAD
print(' ')
print(' iter gmax value ')
print(' ---- --------- ----------------')
print("%4i %9.2e %16.8f" % (iter,gmax,value))
if gmax < tol:
print(' Converged!')
=======
print (' ')
print (' iter gmax value ')
print (' ---- --------- ----------------')
@ -674,7 +569,6 @@ def cgmin(func, dfunc, guess, tol, precond=None, reset=None):
if gmax < tol:
print (' Converged!')
>>>>>>> b1a8eff3c66eb281ed1c2ac4f342fd72cd2a7e46
break
if precond:
@ -742,27 +636,16 @@ def cgmin2(func, guess, tol, eps, printvar=None,reset=None):
(value,g,hh) = numderiv(func, x, step, eps)
gmax = max(map(abs,g))
<<<<<<< HEAD
print(' ')
print(' iter gmax value ')
print(' ---- --------- ----------------')
print("%4i %9.2e %16.8f" % (iter,gmax,value))
=======
print (' ')
print (' iter gmax value ')
print (' ---- --------- ----------------')
print ("%4i %9.2e %16.8f" % (iter,gmax,value))
>>>>>>> b1a8eff3c66eb281ed1c2ac4f342fd72cd2a7e46
if (printvar):
printvar(x)
if gmax < tol:
<<<<<<< HEAD
print(' Converged!')
=======
print (' Converged!')
>>>>>>> b1a8eff3c66eb281ed1c2ac4f342fd72cd2a7e46
break
if (iter % reset) == 0:
@ -793,20 +676,12 @@ def cgmin2(func, guess, tol, eps, printvar=None,reset=None):
if alpha == 0.0:
# LS failed, probably due to lack of precision.
if beta != 0.0:
<<<<<<< HEAD
print("LS failed - trying preconditioned steepest descent direction")
=======
print ("LS failed - trying preconditioned steepest descent direction")
>>>>>>> b1a8eff3c66eb281ed1c2ac4f342fd72cd2a7e46
for i in range(n):
s[i] = -g[i]
(alpha,value) = linesearch(func, x, s, dot(s,g), eps)
if alpha == 0.0:
<<<<<<< HEAD
print(" Insufficient precision to proceed further")
=======
print (" Insufficient precision to proceed further")
>>>>>>> b1a8eff3c66eb281ed1c2ac4f342fd72cd2a7e46
break
for i in range(n):
@ -843,24 +718,6 @@ if __name__ == '__main__':
return sum
<<<<<<< HEAD
print('\n\n TESTING QUASI-NR SOLVER \n\n')
quasinr(f, [1.,0.5,0.3,-0.4], 1e-4, 1e-10)
print('\n\n TESTING GC WITH NUM. GRAD. AND DIAG. PRECOND.\n\n')
cgmin2(f, [1.,0.5,0.3,-0.4], 1e-4, 1e-10, reset=20)
print('\n\n TESTING GC WITH ANAL. GRAD. AND WITHOUT OPTIONAL PRECOND.\n\n')
cgmin(f, df, [1.,0.5,0.3,-0.4], 1e-4)
print('\n\n TESTING GC WITH ANAL. GRAD. AND WITH OPTIONAL PRECOND.\n\n')
cgmin(f, df, [1.,0.5,0.3,-0.4], 1e-4, precond=precond)
print('\n\n TESTING GC WITH ANAL. GRAD. AND NO PRECOND.\n\n')
cgminold(f, df, [1.,0.5,0.3,-0.4], 1e-4)
print('\n\n TESTING JACOBI EIGENSOLVER\n\n')
=======
print ('\n\n TESTING QUASI-NR SOLVER \n\n')
quasinr(f, [1.,0.5,0.3,-0.4], 1e-4, 1e-10)
@ -877,7 +734,6 @@ if __name__ == '__main__':
cgminold(f, df, [1.,0.5,0.3,-0.4], 1e-4)
print ('\n\n TESTING JACOBI EIGENSOLVER\n\n')
>>>>>>> b1a8eff3c66eb281ed1c2ac4f342fd72cd2a7e46
n = 50
a = zeromatrix(n,n)
for i in range(n):
@ -886,11 +742,7 @@ if __name__ == '__main__':
(v,e)= jacobi(a)
<<<<<<< HEAD
print(' eigenvalues')
=======
print (' eigenvalues')
>>>>>>> b1a8eff3c66eb281ed1c2ac4f342fd72cd2a7e46
printvector(e)
#print ' v '
#printmatrix(v)
@ -901,8 +753,4 @@ if __name__ == '__main__':
err = max(err,abs(ev[i][j] - e[i]*v[i][j]))
err = err/(n*max(1.0,abs(e[i])))
if err > 1e-12:
<<<<<<< HEAD
print(' Error in eigenvector ', i, err)
=======
print (' Error in eigenvector ', i, err)
>>>>>>> b1a8eff3c66eb281ed1c2ac4f342fd72cd2a7e46