Python scipy.optimize.newton_krylov() Examples

def solve(self):
    """Solve the selfconsistent problem"""
    return
    self.iterate() # first iteration
    def r2z(v): # convert real to complex
      n = v.shape[0]//2
      return v[0:n] + 1j*v[n:2*n]
    def z2r(v): # convert complex to real
      return np.concatenate([v.real,v.imag])

    def fopt(cij): # function to return the error
      self.cij = r2z(cij) # store this vector
      self.cij2v()  # update the v vectors
      self.iterate() # do an iteration
      print(self.cij)
      self.v2cij() # convert the v to cij
      return z2r(self.cij)-cij
    cij = optimize.broyden1(fopt,z2r(self.cij),f_tol=1e-8,max_rank=10)
#    cij = optimize.fsolve(fopt,z2r(self.cij),xtol=1e-8)
#    cij = optimize.anderson(fopt,z2r(self.cij),f_tol=1e-6,w0=0.1)
#    cij = optimize.newton_krylov(fopt,z2r(self.cij),f_tol=1e-6,outer_k=8)
    self.cij = cij
    self.cij2v() # update 

def newton_krylov_solve(simulation, Estart=None, conv_threshold=1e-10, max_num_iter=50,
				 averaging=True, solver=DEFAULT_SOLVER, jac_solver='c2r',
				 matrix_format=DEFAULT_MATRIX_FORMAT):
	# THIS DOESNT WORK YET! 

	# Stores convergence parameters
	conv_array = np.zeros((max_num_iter, 1))

	# num. columns and rows of A
	Nbig = simulation.Nx*simulation.Ny

	# Defne the starting field for the simulation
	if Estart is None:
		if simulation.fields['Ez'] is None:
			(_, _, E0) = simulation.solve_fields()
		else:
			E0 = deepcopy(simulation.fields['Ez'])
	else:
		E0 = Estart

	E0 = np.reshape(E0, (-1,))

	# funtion for roots
	def _f(E):
		E = np.reshape(E, simulation.eps_r.shape)
		fx = nl_eq_and_jac(simulation, Ez=E, compute_jac=False, matrix_format=matrix_format)
		return np.reshape(fx, (-1,))

	print(_f(E0))

	E_nl = newton_krylov(_f, E0, verbose=True)

	# I'm returning these in place of Hx and Hy to not break things
	return (E_nl, E_nl, E_nl, conv_array) 

def attractive_hubbard(h0,mf=None,mix=0.9,g=0.0,nk=8,solver="plain",
        maxerror=1e-5,**kwargs):
    """Perform the SCF mean field"""
    if not h0.check_mode("spinless"): raise # sanity check
    h = h0.copy() # initial Hamiltonian
    if mf is None:
      try: dold = inout.load(mf_file) # load the file
      except: dold = np.random.random(h.intra.shape[0]) # random guess
    else: dold = mf # initial guess
    ii = 0
    os.system("rm -f STOP") # remove stop file
    def f(dold):
      """Function to minimize"""
#      print("Iteration #",ii) # Iteration
      h = h0.copy() # copy Hamiltonian
      if os.path.exists("STOP"): return dold
      h.add_swave(dold*g) # add the pairing to the Hamiltonian
      t0 = time.time()
      d = onsite_delta_vev(h,nk=nk,**kwargs) # compute the pairing
      t1 = time.time()
      print("Time in this iteration = ",t1-t0) # Difference
      diff = np.max(np.abs(d-dold)) # compute the difference
      print("Error = ",diff) # Difference
      print("Average Pairing = ",np.mean(np.abs(d))) # Pairing
      print("Maximum Pairing = ",np.max(np.abs(d))) # Pairing
      print()
#      ii += 1
      return d
    if solver=="plain":
      do_scf = True
      while do_scf:
        d = f(dold) # new vector
        dold = mix*d + (1-mix)*dold # redefine
        diff = np.max(np.abs(d-dold)) # compute the difference
        if diff<maxerror: 
          do_scf = False
    else:
        print("Solver used:",solver)
        import scipy.optimize as optimize 
        if solver=="newton": fsolver = optimize.newton_krylov
        elif solver=="anderson": fsolver = optimize.anderson
        elif solver=="broyden": fsolver = optimize.broyden2
        elif solver=="linear": fsolver = optimize.linearmixing
        else: raise
        def fsol(x): return x - f(x) # function to solve
        dold = fsolver(fsol,dold,f_tol=maxerror)
    h = h0.copy() # copy Hamiltonian
    h.add_swave(dold*g) # add the pairing to the Hamiltonian
    inout.save(dold,mf_file) # save the mean field
    scf = SCF() # create object
    scf.hamiltonian = h # store
    return scf # return SCF object