Page MenuHomec4science
Files F66393754

matrix_greedy.py
No OneTemporary
Actions

File Metadata

Created
Mon, Jun 10, 06:09

matrix_greedy.py
View Options

import numpy as np
import scipy.sparse as sp
from matplotlib import pyplot as plt
from rrompy.hfengines.base import MatrixEngineBase as MEB
from rrompy.reduction_methods.distributed_greedy import \
RationalInterpolantGreedy as Pade
from rrompy.reduction_methods.distributed_greedy import \
RBDistributedGreedy as RB
test = 1
timed = False
method = "Pade"
method = "RB"
verb = 200
errorEstimatorKind = "BARE"
#errorEstimatorKind = "BASIC"
#errorEstimatorKind = "EXACT"
N = 100
solver = MEB(verbosity = verb)
solver.npar = 1
solver.nAs = 2
if test == 1:
solver.As = [sp.spdiags([np.arange(1, 1 + N)], [0], N, N),
- sp.eye(N)]
elif test == 2:
solver.setSolver("SOLVE")
fftB = np.fft.fft(np.eye(N)) * N**-.5
solver.As = [fftB.dot(np.multiply(np.arange(1, 1 + N), fftB.conj()).T),
- np.eye(N)]
np.random.seed(420)
solver.nbs = 1
solver.bs = [np.random.randn(N) + 1.j * np.random.randn(N)]
mu0 = 10.25
murange = [1.25, 19.25]
mutars = np.linspace(murange[0], murange[1], 500)
if method == "Pade":
params = {'muBounds':murange, 'nTestPoints':200, 'Delta':0, 'S':5,
'greedyTol':1e-2, 'polybasis':"CHEBYSHEV",
'errorEstimatorKind':errorEstimatorKind}
approx = Pade(solver, mu0 = mu0, approxParameters = params,
verbosity = verb)
elif method == "RB":
params = {'muBounds':murange, 'nTestPoints':500, 'greedyTol':1e-2, 'S':5}
approx = RB(solver, mu0 = mu0, approxParameters = params,
verbosity = verb)
if timed:
from time import clock
start_time = clock()
approx.greedy()
print("Time: ", clock() - start_time)
else:
approx.greedy(True)
approx.samplingEngine.verbosity = 0
approx.trainedModel.verbosity = 0
approx.verbosity = 0
normApp = np.zeros(len(mutars))
norm = np.zeros_like(normApp)
res = np.zeros_like(normApp)
err = np.zeros_like(normApp)
for j in range(len(mutars)):
normApp[j] = approx.normApprox(mutars[j])
norm[j] = approx.normHF(mutars[j])
res[j] = (approx.estimatorNormEngine.norm(approx.getRes(mutars[j]))
/ approx.estimatorNormEngine.norm(approx.getRHS(mutars[j])))
err[j] = approx.normErr(mutars[j]) / approx.normHF(mutars[j])
resApp = approx.errorEstimator(mutars)
plt.figure()
plt.semilogy(mutars, norm)
plt.semilogy(mutars, normApp, '--')
plt.semilogy(np.real(approx.mus(0)),
1.25*np.max(norm)*np.ones(approx.mus.shape, dtype = float), 'rx')
plt.xlim(murange)
plt.grid()
plt.show()
plt.close()
plt.figure()
plt.semilogy(mutars, res)
plt.semilogy(mutars, resApp, '--')
plt.semilogy(np.real(approx.mus(0)),
4.*np.max(resApp)*np.ones(approx.mus.shape, dtype = float), 'rx')
plt.xlim(murange)
plt.grid()
plt.show()
plt.close()
plt.figure()
plt.semilogy(mutars, err)
plt.xlim(murange)
plt.grid()
plt.show()
plt.close()
polesTrue = np.arange(1, 1 + N)
polesTrue = polesTrue[polesTrue >= murange[0]]
polesTrue = polesTrue[polesTrue <= murange[1]]
polesApp = approx.getPoles()
mask = (np.real(polesApp) < murange[0]) | (np.real(polesApp) > murange[1])
print("Outliers:", polesApp[mask])
polesApp = polesApp[~mask]
plt.figure()
plt.plot(np.real(polesApp), np.imag(polesApp), 'kx')
plt.plot(polesTrue, np.zeros_like(polesTrue), 'm.')
plt.axis('equal')
plt.grid()
plt.show()
plt.close()

Event Timeline

c4science · Help