diff --git a/contrib/python/mathutil.py b/contrib/python/mathutil.py index 4e82c957b6..dd80cfed35 100644 --- a/contrib/python/mathutil.py +++ b/contrib/python/mathutil.py @@ -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