Commit 4487046a authored by Guillaume Anciaux's avatar Guillaume Anciaux
Browse files

Merge branch 'master' of gitlab.epfl.ch:anciaux/sp4e

parents 60c35af2 52ab38fc
#!/usr/bin/env python3
import numpy as np
def solve_GMRES(A, b, x0, max_iter=30, tol=1e-15, callback=None):
# Initialisation
r = b - np.einsum('ij, j->i', A, x0)
h = np.zeros((max_iter+1, max_iter))
q = np.array([0 for i in range(max_iter)], dtype=object)
q[0] = r/np.linalg.norm(r)
n = min(max_iter, len(b))
for k in range(n):
# Arnoldi iteration
Q, h = arnoldi_iteration(A, q, h, k, max_iter)
# least square solution
x, sol, dependent_variable = least_square_solution(
r, max_iter, h, q, x0, callback=callback)
# Calcul of the residual
residual = np.linalg.norm(np.einsum('ij, j->i', h, sol)-np.linalg.norm(dependent_variable))/np.linalg.norm(dependent_variable)
if residual < tol:
return x, 0
return x, 1
def arnoldi_iteration(A, q, h, k, max_iter):
r_new = np.einsum('ij, j -> i', A, q[k])
for _l in range(k+1):
h[_l, k] = np.einsum('i, i', q[_l], r_new)
r_new = r_new - q[_l]*h[_l, k]
h[k+1, k] = np.linalg.norm(r_new)
if (h[k + 1, k] != 0 and k != max_iter - 1):
q[k + 1] = r_new/h[k + 1, k]
return q, h
def least_square_solution(r, max_iter, h, q, x0, callback=None):
dependent_variable = np.zeros(max_iter + 1)
dependent_variable[0] = np.linalg.norm(r)
sol = np.linalg.lstsq(h, dependent_variable, rcond=None)[0]
x = np.dot(np.transpose(q), sol) + x0
if callback is not None:
callback(x)
return x, sol, dependent_variable
# Instructions to launch the program
- In order to launch the program, simply launch:
```python
python3 ./main.py a_11 a_12 a_21 a_22 b_1 b_2
```
With the coefficients $a_ij$ and $b_i$ the terms of the matrix $A$ and vector $b$ in the quadratic form $1/2 x^T A x - x^T b$, as advertised by the help:
```python
python3 ./main.py -h
```
This will compute the minimum of the quadratic function using one optimizer.
- In order to change the method you can add an additional argument:
```python
python3 ./main.py a_11 a_12 a_21 a_22 b_1 b_2 -m BFGS
python3 ./main.py a_11 a_12 a_21 a_22 b_1 b_2 -m GMRES
python3 ./main.py a_11 a_12 a_21 a_22 b_1 b_2 -m my_GMRES
```
- In order to make a plot, use the following option:
```python
python3 ./main.py a_11 a_12 a_21 a_22 b_1 b_2 -p 1
```
Example of launch:
```python
python3 ./main.py -m my_GMRES 8 1 1 3 2 4 -p 1
```
\ No newline at end of file
#!/usr/bin/env python3
from optimizer import optimize
import numpy as np
import argparse
################################################################
# main function: parsing of program arguments
################################################################
if __name__ == '__main__':
parser = argparse.ArgumentParser(
description='Optimization exercise in python')
parser.add_argument(
'coeffs', type=float,
nargs='+',
help=('Specify the coefficients of the quadratic function,'
' in the following order: a_11, a_12, a_21, a_22, b_1, b_2.'
' The coefficients are the terms of the matrix A'
' and the vector b in the quadratic form 1/2 x^T A x - x^T b.'))
parser.add_argument(
'-m', '--method', type=str,
help=('Specify the method to use for the optimization.'
' Can be BFGS/GMRES/my_GMRES'), default='BFGS')
parser.add_argument(
'-p', '--plot', type=bool,
help=('Choose to plot or not the solution.'
' True/False'), default=False)
args = parser.parse_args()
A = np.array(args.coeffs[:4]).reshape(2, 2)
b = np.array(args.coeffs[4:])
# transform to dictionary
kwargs = vars(args)
del kwargs['coeffs']
# definitions of the quadratic form
def func(X, A, b):
"functor calculation of the quadratic form"
X = X.reshape(X.flatten().shape[0]//2, 2)
# computes the quadratic form with einsum
res = 0.5*np.einsum('ai,ij,aj->a', X, A, X)
res -= np.einsum('i,ai->a', b, X)
return res
optimize(A, b, func=lambda x: func(x, A, b),
jac=lambda x: A@x-b,
hess=lambda x: A,
**kwargs
)
#!/usr/bin/env python3
import matplotlib.pyplot as plt
import numpy as np
from mpl_toolkits.mplot3d import Axes3D
from scipy.optimize import minimize as scipySolve
from scipy.sparse.linalg import lgmres as gmresSolve
from GMRES import solve_GMRES
def plotFunction(func, iterations, title):
'Plotting function for plane and iterations'
x = np.linspace(-3, 3, 100)
y = np.linspace(-3, 3, 100)
x, y = np.meshgrid(x, y)
X = np.stack([x.flatten(), y.flatten()], axis=1)
z = func(X).reshape((100, 100))
xy = np.array(iterations)
zs = func(xy)
fig = plt.figure()
ax = Axes3D(fig)
ax.plot_surface(x, y, z, linewidth=0, antialiased=True,
cmap=plt.cm.viridis, alpha=0.2)
ax.contour(x, y, z, 10, colors="k", linestyles="solid")
ax.plot(xy[:, 0], xy[:, 1], zs, 'ro--')
ax.grid(False)
ax.view_init(elev=62., azim=143)
ax.set_title(title)
plt.xlabel('$x$')
plt.ylabel('$y$')
plt.show()
################################################################
# generic optimization program
################################################################
def optimize(A, b, method='BFGS', plot=False, func=None, **kwargs):
if func is None:
raise RuntimeError('Need a function to optimize')
tol = 1e-6
x0 = np.array([0, 0])
if method == 'BFGS':
optimizer = scipySolve
# avoids a warning
if 'hess' in kwargs:
del kwargs['hess']
elif method == 'GMRES':
optimizer = gmresSolve
elif method == 'my_GMRES':
optimizer = solve_GMRES
else:
raise RuntimeError('unknown solve method: ' + str(method))
if method == 'BFGS':
iterations = []
iterations.append(x0)
res = optimizer(func, x0, tol=tol,
callback=lambda Xi: iterations.append(Xi), **kwargs)
if method == 'GMRES' or method == 'my_GMRES':
iterations = []
iterations.append(x0)
res, info = optimizer(A, b, x0, tol=tol,
callback=lambda Xi: iterations.append(Xi))
print('Method: ', method, ' Solution :', res)
if plot is True:
plotFunction(func, iterations, method)
from __future__ import print_function
from series import Series
################################################################
class ComputeArithmetic(Series):
def __init__(self):
Series.__init__(self)
def __del__(self):
pass
def getAnalyticPrediction(self):
return 1.*self.current_term*(self.current_term+1)/2.
def computeTerm(self, k):
return 1.*k
from __future__ import print_function
import math
from series import Series
################################################################
class ComputePi(Series):
def __init__(self,):
Series.__init__(self)
self.current_value = 0
def __del__(self):
pass
def compute(self, N):
Series.compute(self, N)
return math.sqrt(6.*self.current_value)
def computeTerm(self, k):
return 1./(1.*k*k)
def getAnalyticPrediction(self):
return math.pi
import sys
from abc import ABC, abstractmethod
class DumperSeries(ABC):
def __init__(self, series):
self.series = series
self.precision = 4
self._file = sys.stdout
self.outfile = ''
def __del__(self):
pass
@abstractmethod
def dump(self):
pass
def setPrecision(self, precision):
self.precision = precision
def setOutFile(self, outfile):
self.outfile = outfile
if outfile == '':
return
self._file = open(outfile, 'w')
def __str__(self):
return self.dump()
################################################################
from __future__ import print_function
import sys
import argparse
from compute_pi import ComputePi
from compute_arithmetic import ComputeArithmetic
from print_series import PrintSeries
from plot_series import PlotSeries
################################################################
def main():
parser = argparse.ArgumentParser(
description='Series exercise in python')
parser.add_argument('-s', '--series_type', type=str,
help=('pi or arithmetic'),
required=True)
parser.add_argument('-d', '--dumper_type', type=str,
help=('print or plot'),
required=True)
parser.add_argument('-n', '--maxiter', type=int, required=True,
help='number of loop iteration to compute the series')
parser.add_argument('-f', '--frequency', type=int, default=1,
help='frequency at which dumps/plots are made')
parser.add_argument('-o', '--outfile', type=str, default="",
help='file where to store the output')
args = parser.parse_args()
series_type = args.series_type
dumper_type = args.dumper_type
maxiter = args.maxiter
freq = args.frequency
outfile = args.outfile
#################################################################
# creation of ojects
#################################################################
if series_type == "pi":
serie_object = ComputePi()
elif series_type == "arithmetic":
serie_object = ComputeArithmetic()
else:
raise RuntimeError("unknown series type: " + series_type)
if dumper_type == "print":
dumper_object = PrintSeries(serie_object, maxiter, freq)
elif dumper_type == "plot":
dumper_object = PlotSeries(serie_object, maxiter, freq)
else:
print("unknown dumper type: " + dumper_type)
sys.exit(-1)
dumper_object.setPrecision(20)
dumper_object.setOutFile(outfile)
#################################################################
# start of the run
#################################################################
dumper_object.dump()
if __name__ == "__main__":
main()
from __future__ import print_function
from dumper_series import DumperSeries
import matplotlib.pyplot as plt
################################################################
class PlotSeries(DumperSeries):
def __init__(self, series, maxiter, freq):
DumperSeries.__init__(self, series)
self.maxiter = maxiter
self.freq = freq
def __del__(self):
pass
def dump(self):
fig = plt.figure()
axe = fig.add_subplot(1, 1, 1)
x = []
numerical = []
analytic = []
for i in range(1, int(self.maxiter/self.freq)):
res = self.series.compute(i*self.freq-1)
x.append(i*self.freq)
numerical.append(res)
analytic.append(self.series.getAnalyticPrediction())
axe.plot(x, numerical, marker='o', label='Numerical')
axe.plot(x, analytic, label='Analytical')
axe.set_xlabel(r'$k$')
axe.set_ylabel(r'Series')
axe.legend()
if self.outfile == '':
plt.show()
else:
plt.savefig(self.outfile)
return "Plotted"
from __future__ import print_function
import math
from dumper_series import DumperSeries
################################################################
class PrintSeries(DumperSeries):
def __init__(self, series, maxiter, freq):
DumperSeries.__init__(self, series)
self.maxiter = maxiter
self.freq = freq
def __del__(self):
pass
def dump(self):
out = ""
for i in range(1, int(self.maxiter/self.freq)):
res = self.series.compute(i*self.freq-1)
res2 = self.series.compute(i*self.freq)
my_format = '{:.' + str(self.precision) + 'e}'
out += str(i*self.freq) + " "
out += my_format.format(res) + " "
out += my_format.format(res2 - res)
try:
out += " " + my_format.format(
math.fabs(res2 - self.series.getAnalyticPrediction()))
except Exception:
pass
out += "\n"
self._file.write(out)
from abc import ABC, abstractmethod
class Series(ABC):
def __init__(self):
self.current_term = 0
self.current_value = 0
def __del__(self):
pass
def compute(self, N):
if self.current_term <= N:
N -= self.current_term
else:
self.current_value = 0.
self.current_term = 0
for k in range(0, N):
self.addTerm()
return self.current_value
@abstractmethod
def getAnalyticPrediction(self):
pass
def addTerm(self):
self.current_term += 1
self.current_value += self.computeTerm(self.current_term)
def computeTerm(self, k):
pass
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment