Homec4science
Diffusion Homework_sp4e_Peruzzo_SáezUribe (master) Browse

R9490/exercise_2.py
e18f3edf424emaster

Last Change
Raw File

exercise_2.py
View Options

"""
EXERCISE 2 - part 2 and 3
Compare the Conjugate Gradient that has been implemented with the minimizer provided by SCIPY
"""
import argparse
import numpy as np
from conjugate_gradient import conjugate_gradient
from optimizer import optimizer
import plot
# Initiate and describe parser
parser = argparse.ArgumentParser(description='Find the minimizer of symmetric definite systems of linear equations by using the conjugate gradient method')
# Define type of arguments to be received
# Receive a list of elements to build matrix A and vector b
parser.add_argument('elements', metavar='N', type=float, nargs='+',
help='enter 6 numbers: first 4 for matrix A in the order (1,1) (1,2) (2,1) (2,2), and last 2 for vector b. NOTE: matrix A in exercise 2 is equal to 2 * matrix of size 2x2 defined in exercise 1')
# Receive the type of minimizer
parser.add_argument('--method', default = 'scipy',
help='choose method to perform minimization, options are: "ourCG" or "scipy", if you do not define method then by default is "scipy"')
# Receive instruction whether plor or not
parser.add_argument('--plot', default = 'not',
help='indicate if you want to plot, options are: "yes" or "not", if you do not indicate any option then by default is "not"')
# Get arguments
args = parser.parse_args()
# Build matrix A and vector b
A = np.array([[args.elements[0], args.elements[1]], [args.elements[2], args.elements[3]]])
b = np.array([args.elements[4], args.elements[5]])
# NOTE: matrix A here (in exercise 2) is equal to 2 * matrix os size 2x2 defined in exercise 1
# Construct the quadratic function S(x,y)
def Sx(xy):
A = np.array([[args.elements[0], args.elements[1]], [args.elements[2], args.elements[3]]])/2
b = np.array([args.elements[4], args.elements[5]])
x = np.asarray(xy[0])
y = np.asarray(xy[1])
xy = np.array([x, y])
res = xy.T @ A @ xy - xy.T @ b
return res
# Initial guess, tolerance and maximum number of iterations for our conjugate gradient method
x0 = np.array([2, 1])
tol = 1e-9
maxIter = 100;
# Call corresponding optimizer function
if args.method == 'ourCG':
res = conjugate_gradient(A,b,x0,tol,maxIter)
minimizer = res[0]
points = np.asarray([res[1],res[2]]).transpose()
method="CG"
else:
method = "BFGS"
minimizer, resultINFO, points = optimizer(Sx, method, x0, tol)
# Plot (in case)
if args.plot == 'yes':
fig = plot.plotminimizer(Sx, points,method)
else:
print('not plot requested')
# Print result
print('_________________________________________________________')
print('The minimizer is: x = ', minimizer[0], ', y = ', minimizer[1], 'and S(x,y) = ',
Sx(minimizer))
# Exact solution
exact = np.linalg.inv(A)@b
print('The exact solution is: x = ', exact[0], 'and y = ', exact[1],
'and S(x,y) = ',Sx([exact[0], exact[1]]))
c4science · Help