diff --git a/README.md b/README.md
index 2ae0ba7..4998fbb 100644
--- a/README.md
+++ b/README.md
@@ -1,31 +1,36 @@
RROMPy
============
Module for the solution and rational model order reduction of parametric PDE-based problem. Coded in Python 3.6.
## Prerequisites
-**RBniCS** requires
-* **numpy**.
-* **scipy**.
+**RROMPy** requires
+* **numpy**;
+* **scipy**;
+* **fenics**.
+
+### Fenics
+Most of the PDE engines already provided rely on [FEniCS](http://fenicsproject.org/). If you do not have FEniCS installed, you may want to create an [Anaconda3/Miniconda3](http://anaconda.org/) environment using the provided [file](conda-fenics.yml) by running the command
+```
+conda env create --file conda-fenics.yml
+```
+This will create an environment where Fenics can be used. In order to use FEniCS, the environment must be activated through
+```
+source activate fenicsenv
+```
## Installing
Clone the repository
```
git clone https://c4science.ch/source/RROMPy.git
```
and install the package by typing
```
python3 setup.py install
```
-### Fenics
-Most of the PDE engines already provided rely on [FEniCS](https://fenicsproject.org/). If you do not have FEniCS installed, you may want to create an [Anaconda3/Miniconda3](http://anaconda.org/) environment using the provided [file](./conda-fenics.yml) by running the command
-```
-conda env create --file conda-fenics.yml
-```
-This will create an environment where Fenics can be used. In order to use FEniCS, the environment must be activated through
-```
-source activate fenicsenv
-```
## License
-This project is licensed under the GNU GENERAL PUBLIC LICENSE license - see the [LICENSE](./LICENSE) file for details.
+This project is licensed under the GNU GENERAL PUBLIC LICENSE license - see the [LICENSE](LICENSE) file for details.
+
+## Acknowledgments
+Part of the funding that made this module possible has been provided by the Swiss National Science Foundation through the FNS Research Project No. 200021.
diff --git a/examples/scipy/solver.py b/examples/base/solver.py
similarity index 100%
rename from examples/scipy/solver.py
rename to examples/base/solver.py
diff --git a/examples/greedy/squareBubbleHomog.py b/examples/greedy/squareBubbleHomog.py
new file mode 100644
index 0000000..53abf14
--- /dev/null
+++ b/examples/greedy/squareBubbleHomog.py
@@ -0,0 +1,95 @@
+import numpy as np
+from rrompy.hfengines.scipy import HelmholtzSquareBubbleProblemEngine as HSBPE
+from rrompy.reduction_methods.lagrange_greedy import ApproximantLagrangePadeGreedy as Pade
+from rrompy.reduction_methods.lagrange_greedy import ApproximantLagrangeRBGreedy as RB
+from rrompy.utilities.base import squareResonances
+
+verb = 2
+algo = "Pade"
+#algo = "RB"
+
+k0s = np.power(np.linspace(95, 127, 25), .5)
+k0 = np.mean(np.power(k0s, 2.)) ** .5
+kl, kr = min(k0s), max(k0s)
+
+polesexact = np.unique(np.power(squareResonances(kl**2., kr**2., False), .5))
+
+rescaling = lambda x: np.power(x, 2.)
+rescalingInv = lambda x: np.power(x, .5)
+params = {'muBounds':[kl, kr], 'nTrainingPoints':500,
+ 'greedyTol':1e-2, 'nTestPoints':5}
+
+solver = HSBPE(kappa = 12 ** .5, theta = np.pi / 3, n = 30,
+ verbosity = verb)
+solver.omega = np.real(k0)
+if algo == "Pade":
+ approx = Pade(solver, mu0 = k0, approxParameters = params,
+ verbosity = verb)
+else:
+ approx = RB(solver, mu0 = k0, approxParameters = params, verbosity = verb)
+
+#from time import clock
+#start_time = clock()
+approx.greedy(True)
+#print("Time: ", clock() - start_time)
+
+approx.samplingEngine.verbosity = 0
+approx.verbosity = 0
+from matplotlib import pyplot as plt
+MassInv = approx.massInv
+normApp = np.zeros_like(k0s)
+norm = np.zeros_like(k0s)
+res = np.zeros_like(k0s)
+err = np.zeros_like(k0s)
+for j in range(len(k0s)):
+ normApp[j] = approx.normApp(k0s[j])
+ norm[j] = approx.normHF(k0s[j])
+ res[j] = (solver.norm(MassInv.solve(approx.getRes(k0s[j])))
+ / solver.norm(MassInv.solve(approx.getRHS(k0s[j]))))
+ err[j] = approx.normErr(k0s[j]) / approx.normHF(k0s[j])
+resApp = approx.errorEstimator(k0s)
+
+plt.figure()
+plt.semilogy(k0s, norm)
+plt.semilogy(k0s, normApp, '--')
+plt.semilogy(polesexact,
+ 2.*np.max(norm)*np.ones_like(polesexact, dtype = float), 'm.')
+plt.semilogy(np.real(approx.mus),
+ 4.*np.max(norm)*np.ones_like(approx.mus, dtype = float), 'rx')
+plt.xlim([kl, kr])
+plt.grid()
+plt.show()
+plt.close()
+
+plt.figure()
+plt.semilogy(k0s, res)
+plt.semilogy(k0s, resApp, '--')
+plt.semilogy(polesexact,
+ 2.*np.max(resApp)*np.ones_like(polesexact, dtype = float), 'm.')
+plt.semilogy(np.real(approx.mus),
+ 4.*np.max(resApp)*np.ones_like(approx.mus, dtype = float), 'rx')
+plt.xlim([kl, kr])
+plt.grid()
+plt.show()
+plt.close()
+
+plt.figure()
+plt.semilogy(k0s, err)
+plt.semilogy(polesexact,
+ 2.*np.max(err)*np.ones_like(polesexact, dtype = float), 'm.')
+plt.xlim([kl, kr])
+plt.grid()
+plt.show()
+plt.close()
+
+polesApp = approx.getPoles()
+plt.figure()
+plt.plot(np.real(polesApp), np.imag(polesApp), 'kx')
+plt.plot(np.real(polesexact), np.imag(polesexact), 'm.')
+plt.xlim(kl, kr)
+plt.ylim(- (kr - kl) * .5, + (kr - kl) * .5)
+#plt.axis('scaled')
+plt.grid()
+plt.show()
+plt.close()
+
diff --git a/examples/greedy/squareTransmissionNonHomog.py b/examples/greedy/squareTransmissionNonHomog.py
new file mode 100644
index 0000000..51a00f0
--- /dev/null
+++ b/examples/greedy/squareTransmissionNonHomog.py
@@ -0,0 +1,77 @@
+import numpy as np
+from rrompy.hfengines.scipy import HelmholtzSquareTransmissionProblemEngine as HSTPE
+from rrompy.reduction_methods.lagrange_greedy import ApproximantLagrangePadeGreedy as Pade
+
+verb = 7
+homog = True
+#homog = False
+
+k0s = np.power(np.linspace(4, 15, 200), .5)
+k0 = np.mean(np.power(k0s, 2.)) ** .5
+kl, kr = min(k0s), max(k0s)
+
+rescaling = lambda x: np.power(x, 2.)
+rescalingInv = lambda x: np.power(x, .5)
+params = {'muBounds':[kl, kr], 'nTrainingPoints':500,
+ 'greedyTol':1e-3, 'nTestPoints':10}
+
+solver = HSTPE(nT = 1, nB = 2, theta = np.pi * 20 / 180, kappa = 4.,
+ n = 20, verbosity = verb)
+solver.omega = np.real(k0)
+approx = Pade(solver, mu0 = k0, approxParameters = params,
+ verbosity = verb, homogeneized = homog)
+approx.greedy(True)
+
+approx.samplingEngine.verbosity = 0
+approx.verbosity = 0
+from matplotlib import pyplot as plt
+MassInv = approx.massInv
+normApp = np.zeros_like(k0s)
+norm = np.zeros_like(k0s)
+res = np.zeros_like(k0s)
+err = np.zeros_like(k0s)
+for j in range(len(k0s)):
+ normApp[j] = approx.normApp(k0s[j])
+ norm[j] = approx.normHF(k0s[j])
+ res[j] = (solver.norm(MassInv.solve(approx.getRes(k0s[j],
+ homogeneized = homog)))
+ / solver.norm(MassInv.solve(approx.getRHS(k0s[j],
+ homogeneized = homog))))
+ err[j] = (approx.normErr(k0s[j], homogeneized = homog)
+ / approx.normHF(k0s[j], homogeneized = homog))
+resApp = approx.errorEstimator(k0s, homogeneized = homog)
+
+polesApp = approx.getPoles()
+polesApp = polesApp[np.abs(np.imag(polesApp)) < 1e-3]
+plt.figure()
+plt.semilogy(k0s, norm)
+plt.semilogy(k0s, normApp, '--')
+plt.semilogy(np.real(approx.mus),
+ 4.*np.max(norm)*np.ones_like(approx.mus, dtype = float), 'rx')
+plt.semilogy(np.real(polesApp),
+ 2.*np.max(norm)*np.ones_like(polesApp, dtype = float), 'k.')
+plt.xlim([kl, kr])
+plt.grid()
+plt.show()
+plt.close()
+
+plt.figure()
+plt.semilogy(k0s, res)
+plt.semilogy(k0s, resApp, '--')
+plt.semilogy(np.real(approx.mus),
+ 4.*np.max(resApp)*np.ones_like(approx.mus, dtype = float), 'rx')
+plt.semilogy(np.real(polesApp),
+ 2.*np.max(resApp)*np.ones_like(polesApp, dtype = float), 'k.')
+plt.xlim([kl, kr])
+plt.grid()
+plt.show()
+plt.close()
+
+plt.figure()
+plt.semilogy(k0s, err)
+plt.semilogy(np.real(polesApp),
+ 2.*np.max(err)*np.ones_like(polesApp, dtype = float), 'k.')
+plt.xlim([kl, kr])
+plt.grid()
+plt.show()
+plt.close()
diff --git a/examples/pod/LagrangePoles.py b/examples/pod/LagrangePoles.py
new file mode 100644
index 0000000..141a5f6
--- /dev/null
+++ b/examples/pod/LagrangePoles.py
@@ -0,0 +1,48 @@
+from matplotlib import pyplot as plt
+import numpy as np
+from rrompy.hfengines.scipy import HelmholtzSquareBubbleProblemEngine as HSBPE
+from rrompy.reduction_methods.lagrange import ApproximantLagrangePade as Pade
+from rrompy.utilities.parameter_sampling import QuadratureSampler as QS
+from rrompy.utilities.base import squareResonances
+
+verb = 0
+sampling = "Uniform"
+
+ks = [1, 46 ** .5]
+solver = HSBPE(kappa = 12 ** .5, theta = np.pi / 3, n = 20,
+ verbosity = verb)
+k0 = np.mean(np.power(ks, 2.)) ** .5
+k0 = 3.46104724
+solver.omega = np.real(k0)
+
+rescaling = lambda x: np.power(x, 2.)
+rescalingInv = lambda x: np.power(x, .5)
+if sampling == "Uniform":
+ sampler = QS(ks, "UNIFORM", rescaling, rescalingInv)
+
+nsets = 15
+paramsPade = {'S':2, 'POD':True, 'sampler':sampler}
+approx = Pade(solver, mu0 = k0, approxParameters = paramsPade,
+ verbosity = verb)
+
+poles = [None] * (nsets - 1)
+polesexact = np.unique(np.power(squareResonances(ks[0]**2., ks[1]**2., False),
+ .5))
+for i in range(1, nsets):
+ print("N = {}".format(4 * i))
+ approx.approxParameters = {'N': 4 * i, 'M': 4 * i, 'S': 4 * i + 1}
+ approx.setupApprox()
+ poles[i - 1] = approx.getPoles()
+
+
+for i in range(1, nsets):
+ plt.figure()
+ plt.plot(np.real(poles[i - 1]), np.imag(poles[i - 1]), 'kx')
+ plt.plot(polesexact, np.zeros_like(polesexact), 'm.')
+ plt.plot(k0, 0, 'r*')
+ plt.xlim(ks)
+ plt.ylim((ks[0] - ks[1]) / 2., (ks[1] - ks[0]) / 2.)
+ plt.title("N = {}, Neff = {}".format(4 * i, len(poles[i - 1])))
+ plt.grid()
+ plt.show()
+ plt.close()
\ No newline at end of file
diff --git a/examples/scipy/LagrangeSweep.py b/examples/pod/LagrangeSweep.py
similarity index 91%
rename from examples/scipy/LagrangeSweep.py
rename to examples/pod/LagrangeSweep.py
index d1d09c9..45be186 100644
--- a/examples/scipy/LagrangeSweep.py
+++ b/examples/pod/LagrangeSweep.py
@@ -1,102 +1,102 @@
from copy import copy
import numpy as np
#from rrompy.hfengines.scipy import HelmholtzBoxScatteringProblemEngine as HBSPE
from rrompy.hfengines.scipy import HelmholtzCavityScatteringProblemEngine as HBSPE
from rrompy.reduction_methods.lagrange import ApproximantLagrangePade as Pade
from rrompy.reduction_methods.lagrange import ApproximantLagrangeRB as RB
from rrompy.utilities.parameter_sweeper import ParameterSweeper as Sweeper
from rrompy.utilities.parameter_sampling import QuadratureSampler as QS
from rrompy.utilities.parameter_sampling import WarpingFunction as WF
from rrompy.utilities.parameter_sampling import WarpingSampler as WS
verb = 0
npoints = 100
homog = True
homog = False
LSratio = 2. / 3.
sampling = "Uniform"
#sampling = "Cheby"
#sampling = "SinC"
assert LSratio <= 1. + np.finfo(float).eps
ks = [10 + 0.j, 14 + 0.j]
solver = HBSPE(kappa = 3, n = 15)
solver.omega = np.real(np.mean(ks))
mutars = np.linspace(9, 15, npoints)
homogMSG = "Homog"
if not homog: homogMSG = "Non" + homogMSG
-#filenamebase = '../data/output/HelmholtzBoxLSweep'
-filenamebase = '../data/plots/LagrangeScatteringSquare1p5'
-filenamebase = filenamebase + sampling + "/HelmholtzBoxLSweep" + homogMSG
+filenamebase = '../data/output/ScatteringSquareLSweep'
+#filenamebase = '../data/plots/LagrangeScatteringSquare1p5'
+#filenamebase = filenamebase + sampling + "/HelmholtzBoxLSweep" + homogMSG
k0 = np.mean(ks)
shift = 7
shift = np.int(8 / LSratio - 1)
nsets = 3
stride = np.int(8 / LSratio)
Smax = stride * (nsets - 1) + shift + 2
rescaling = lambda x: np.power(x, 2.)
rescalingInv = lambda x: np.power(x, .5)
if sampling == "Uniform":
sampler = QS(ks, "UNIFORM", rescaling, rescalingInv)
if sampling == "Cheby":
sampler = QS(ks, "CHEBYSHEV", rescaling, rescalingInv)
if sampling == "SinC":
warping = WF(call = lambda x: x - 2. * (1. - LSratio) / np.pi * np.sin(np.pi * x),
repr = "x-{}*sin(pi*x)".format(2. * (1. - LSratio) / np.pi))
sampler = WS(ks, warping, rescaling, rescalingInv)
paramsPade = {'S':Smax, 'POD':True, 'sampler':sampler}
paramsRB = copy(paramsPade)
paramsPoly = copy(paramsPade)
paramsSetsPade = [None] * nsets
paramsSetsRB = [None] * nsets
paramsSetsPoly = [None] * nsets
for i in range(nsets):
paramsSetsPade[i] = {'N': np.int(LSratio * (stride * i + shift + 1)),
'M': np.int(LSratio * (stride * i + shift + 1)),
'S': stride * i + shift + 2}
paramsSetsRB[i] = {'R': np.int(LSratio * (stride * i + shift + 1)),
'S': stride * i + shift + 2}
paramsSetsPoly[i] = {'N': 0,
'M': np.int(LSratio * (stride * i + shift + 1)),
'S': stride * i + shift + 2}
appPade = Pade(solver, mu0 = k0, approxParameters = paramsPade,
- verbosity = verb, homogeneize = homog)
+ verbosity = verb, homogeneized = homog)
appRB = RB(solver, mu0 = k0, approxParameters = paramsRB,
- verbosity = verb, homogeneize = homog)
+ verbosity = verb, homogeneized = homog)
appPoly = Pade(solver, mu0 = k0, approxParameters = paramsPoly,
- verbosity = verb, homogeneize = homog)
+ verbosity = verb, homogeneized = homog)
sweeper = Sweeper(mutars = mutars, mostExpensive = 'Approx')
sweeper.ROMEngine = appPade
sweeper.params = paramsSetsPade
filenamePade = sweeper.sweep(filenamebase + 'Pade.dat')
sweeper.ROMEngine = appRB
sweeper.params = paramsSetsRB
filenameRB = sweeper.sweep(filenamebase + 'RB.dat')
sweeper.ROMEngine = appPoly
sweeper.params = paramsSetsPoly
filenamePoly = sweeper.sweep(filenamebase + 'Poly.dat')
sweeper.plotCompare([filenamePade, filenameRB, filenamePoly], ['muRe'],
['normHF', 'normApp'], ['S'], onePlot = True,
save = filenamebase + 'Norm', ylims = {'top' : 1e1},
saveFormat = "png", labels = ["Pade'", "RB", "Poly"],
# figsize = (5, 3.75))
figsize = (10, 7.5))
sweeper.plotCompare([filenamePade, filenameRB, filenamePoly], ['muRe'],
['normResRel'], ['S'], save = filenamebase + 'Res',
ylims = {'top' : 1e1}, saveFormat = "png",
labels = ["Pade'", "RB", "Poly"],
# figsize = (5, 3.75))
figsize = (10, 7.5))
sweeper.plotCompare([filenamePade, filenameRB, filenamePoly], ['muRe'],
['normErrRel'], ['S'], save = filenamebase + 'Err',
ylims = {'top' : 1e1}, saveFormat = "png",
labels = ["Pade'", "RB", "Poly"],
# figsize = (5, 3.75))
figsize = (10, 7.5))
diff --git a/examples/scipy/LagrangeSweepUnstable.py b/examples/pod/LagrangeSweepUnstable.py
similarity index 100%
rename from examples/scipy/LagrangeSweepUnstable.py
rename to examples/pod/LagrangeSweepUnstable.py
diff --git a/examples/scipy/PadeLagrange.py b/examples/pod/PadeLagrange.py
similarity index 96%
rename from examples/scipy/PadeLagrange.py
rename to examples/pod/PadeLagrange.py
index e3f330f..300165c 100644
--- a/examples/scipy/PadeLagrange.py
+++ b/examples/pod/PadeLagrange.py
@@ -1,107 +1,107 @@
import numpy as np
from rrompy.hfengines.scipy import HelmholtzSquareBubbleProblemEngine as HSBPE
from rrompy.hfengines.scipy import HelmholtzSquareTransmissionProblemEngine as HSTPE
from rrompy.hfengines.scipy import HelmholtzBoxScatteringProblemEngine as HBSPE
from rrompy.reduction_methods.lagrange import ApproximantLagrangePade as Pade
from rrompy.utilities.parameter_sampling import QuadratureSampler as QS
-testNo = 2
+testNo = 1
verb = 0
homog = True
homog = False
if testNo == 1:
k0s = np.power([10 + 0.j, 14 + 0.j], .5)
k0 = np.mean(k0s)
ktar = (11 + .5j) ** .5
rescaling = lambda x: np.power(x, 2.)
rescalingInv = lambda x: np.power(x, .5)
params = {'N':4, 'M':3, 'S':5, 'POD':True,
'sampler':QS(k0s, "CHEBYSHEV", rescaling, rescalingInv)}
solver = HSBPE(kappa = 12 ** .5, theta = np.pi / 3, n = 40,
verbosity = verb)
solver.omega = np.real(k0)
approx = Pade(solver, mu0 = k0, approxParameters = params,
verbosity = verb)
approx.setupApprox()
# approx.plotSamples()
approx.plotApp(ktar, name = 'u_Pade''')
approx.plotHF(ktar, name = 'u_HF')
approx.plotErr(ktar, name = 'err')
approx.plotRes(ktar, name = 'res')
appErr, solNorm = approx.normErr(ktar), approx.normHF(ktar)
resNorm, RHSNorm = approx.normRes(ktar), approx.normRHS(ktar)
print(('SolNorm:\t{}\nErr:\t{}\nErrRel:\t{}').format(solNorm, appErr,
np.divide(appErr, solNorm)))
print(('RHSNorm:\t{}\nRes:\t{}\nResRel:\t{}').format(RHSNorm, resNorm,
np.divide(resNorm, RHSNorm)))
print('\nPoles Pade'':')
print(approx.getPoles())
############
elif testNo == 2:
k0s = [3.85 + 0.j, 4.15 + 0.j]
k0 = np.mean(k0s)
ktar = 4 + .15j
rescaling = lambda x: np.power(x, 2.)
rescalingInv = lambda x: np.power(x, .5)
params = {'N':8, 'M':9, 'S':10, 'POD':True,
'sampler':QS(k0s, "CHEBYSHEV", rescaling, rescalingInv)}
solver = HSTPE(nT = 2, nB = 1, theta = np.pi * 45/180, kappa = 4., n = 50,
verbosity = verb)
solver.omega = np.real(k0)
approx = Pade(solver, mu0 = k0, approxParameters = params,
- verbosity = verb, homogeneize = homog)
+ verbosity = verb, homogeneized = homog)
approx.setupApprox()
# approx.plotSamples()
approx.plotApp(ktar, name = 'u_Pade''')
approx.plotHF(ktar, name = 'u_HF')
approx.plotErr(ktar, name = 'err')
approx.plotRes(ktar, name = 'res')
appErr, solNorm = approx.normErr(ktar), approx.normHF(ktar)
resNorm, RHSNorm = approx.normRes(ktar), approx.normRHS(ktar)
print(('SolNorm:\t{}\nErr:\t{}\nErrRel:\t{}').format(solNorm, appErr,
np.divide(appErr, solNorm)))
print(('RHSNorm:\t{}\nRes:\t{}\nResRel:\t{}').format(RHSNorm, resNorm,
np.divide(resNorm, RHSNorm)))
print('\nPoles Pade'':')
print(approx.getPoles())
############
elif testNo == 3:
k0s = [2, 5]
k0 = np.mean(k0s)
ktar = 4.5 - 0.j
params = {'N':10, 'M':11, 'S':15, 'POD':True,
'sampler':QS(k0s, "CHEBYSHEV")}
solver = HBSPE(R = 7, kappa = 3, theta = - np.pi * 75 / 180, n = 40,
verbosity = verb)
solver.omega = np.real(k0)
approx = Pade(solver, mu0 = k0, approxParameters = params,
- verbosity = verb, homogeneize = homog)
+ verbosity = verb, homogeneized = homog)
approx.setupApprox()
approx.plotSamples()
approx.plotApp(ktar, name = 'u_Pade''')
approx.plotHF(ktar, name = 'u_HF')
approx.plotErr(ktar, name = 'err')
approx.plotRes(ktar, name = 'res')
appErr, solNorm = approx.normErr(ktar), approx.normHF(ktar)
resNorm, RHSNorm = approx.normRes(ktar), approx.normRHS(ktar)
print(('SolNorm:\t{}\nErr:\t{}\nErrRel:\t{}').format(solNorm, appErr,
np.divide(appErr, solNorm)))
print(('RHSNorm:\t{}\nRes:\t{}\nResRel:\t{}').format(RHSNorm, resNorm,
np.divide(resNorm, RHSNorm)))
print('\nPoles Pade'':')
print(approx.getPoles())
diff --git a/examples/scipy/PadeTaylor.py b/examples/pod/PadeTaylor.py
similarity index 96%
rename from examples/scipy/PadeTaylor.py
rename to examples/pod/PadeTaylor.py
index 4fe5a30..5461db0 100644
--- a/examples/scipy/PadeTaylor.py
+++ b/examples/pod/PadeTaylor.py
@@ -1,97 +1,97 @@
import numpy as np
from rrompy.hfengines.scipy import HelmholtzSquareBubbleProblemEngine as HSBPE
from rrompy.hfengines.scipy import (
HelmholtzSquareTransmissionProblemEngine as HSTPE)
from rrompy.hfengines.scipy import HelmholtzBoxScatteringProblemEngine as HBSPE
from rrompy.reduction_methods.taylor import ApproximantTaylorPade as Pade
testNo = 2
verb = 0
homog = True
#homog = False
if testNo == 1:
params = {'N':4, 'M':3, 'E':4, 'sampleType':'Arnoldi', 'POD':True}
k0 = 12 ** .5
ktar = 10.5 ** .5
solver = HSBPE(kappa = 12 ** .5, theta = np.pi / 3, n = 40,
verbosity = verb)
solver.omega = np.real(k0)
approx = Pade(solver, mu0 = k0, approxParameters = params,
verbosity = verb)
approx.setupApprox()
# approx.plotSamples()
approx.plotApp(ktar, name = 'u_Pade''')
approx.plotHF(ktar, name = 'u_HF')
approx.plotErr(ktar, name = 'err')
approx.plotRes(ktar, name = 'res')
appErr, solNorm = approx.normErr(ktar), approx.normHF(ktar)
resNorm, RHSNorm = approx.normRes(ktar), approx.normRHS(ktar)
print(('SolNorm:\t{}\nErr:\t{}\nErrRel:\t{}').format(solNorm, appErr,
np.divide(appErr, solNorm)))
print(('RHSNorm:\t{}\nRes:\t{}\nResRel:\t{}').format(RHSNorm, resNorm,
np.divide(resNorm, RHSNorm)))
print('\nPoles Pade'':')
print(approx.getPoles())
############
elif testNo == 2:
params = {'N':6, 'M':7, 'E':7, 'sampleType':'Arnoldi', 'POD':True}
k0 = 16 ** .5
ktar = 15 ** .5
solver = HSTPE(nT = 2, nB = 1, theta = np.pi * 45/180, kappa = 4., n = 50,
verbosity = verb)
solver.omega = np.real(k0)
approx = Pade(solver, mu0 = k0, approxParameters = params,
- verbosity = verb, homogeneize = homog)
+ verbosity = verb, homogeneized = homog)
approx.setupApprox()
# approx.plotSamples()
approx.plotApp(ktar, name = 'u_Pade''')
approx.plotHF(ktar, name = 'u_HF')
approx.plotErr(ktar, name = 'err')
approx.plotRes(ktar, name = 'res')
appErr, solNorm = approx.normErr(ktar), approx.normHF(ktar)
resNorm, RHSNorm = approx.normRes(ktar), approx.normRHS(ktar)
print(('SolNorm:\t{}\nErr:\t{}\nErrRel:\t{}').format(solNorm, appErr,
np.divide(appErr, solNorm)))
print(('RHSNorm:\t{}\nRes:\t{}\nResRel:\t{}').format(RHSNorm, resNorm,
np.divide(resNorm, RHSNorm)))
print('\nPoles Pade'':')
print(approx.getPoles())
############
elif testNo in [3, 4]:
if testNo == 3:
params = {'N':7, 'M':8, 'E':8, 'sampleType':'Krylov', 'POD':True}
else:
params = {'N':7, 'M':8, 'E':8, 'sampleType':'Arnoldi', 'POD':True}
k0 = 3
ktar = 4.+0.j
solver = HBSPE(R = 5, kappa = 3, theta = - np.pi * 75 / 180, n = 30,
verbosity = verb)
solver.omega = np.real(k0)
approx = Pade(solver, mu0 = k0, approxParameters = params,
- verbosity = verb, homogeneize = homog)
+ verbosity = verb, homogeneized = homog)
approx.setupApprox()
approx.plotSamples()
approx.plotApp(ktar, name = 'u_Pade''')
approx.plotHF(ktar, name = 'u_HF')
approx.plotErr(ktar, name = 'err')
approx.plotRes(ktar, name = 'res')
appErr, solNorm = approx.normErr(ktar), approx.normHF(ktar)
resNorm, RHSNorm = approx.normRes(ktar), approx.normRHS(ktar)
print(('SolNorm:\t{}\nErr:\t{}\nErrRel:\t{}').format(solNorm, appErr,
np.divide(appErr, solNorm)))
print(('RHSNorm:\t{}\nRes:\t{}\nResRel:\t{}').format(RHSNorm, resNorm,
np.divide(resNorm, RHSNorm)))
print('\nPoles Pade'':')
print(approx.getPoles())
diff --git a/examples/scipy/RBLagrange.py b/examples/pod/RBLagrange.py
similarity index 97%
rename from examples/scipy/RBLagrange.py
rename to examples/pod/RBLagrange.py
index bbfe2bc..ac5474f 100644
--- a/examples/scipy/RBLagrange.py
+++ b/examples/pod/RBLagrange.py
@@ -1,99 +1,99 @@
import numpy as np
from rrompy.hfengines.scipy import HelmholtzSquareBubbleProblemEngine as HSBPE
from rrompy.hfengines.scipy import HelmholtzSquareTransmissionProblemEngine as HSTPE
from rrompy.hfengines.scipy import HelmholtzBoxScatteringProblemEngine as HBSPE
from rrompy.reduction_methods.lagrange import ApproximantLagrangeRB as RB
from rrompy.utilities.parameter_sampling import QuadratureSampler as QS
testNo = 3
verb = 0
homog = True
homog = False
if testNo == 1:
k0s = np.power([10 + 0.j, 14 + 0.j], .5)
k0 = np.mean(k0s)
ktar = (11 + .5j) ** .5
rescaling = lambda x: np.power(x, 2.)
rescalingInv = lambda x: np.power(x, .5)
params = {'S':5, 'R':4, 'POD':True,
'sampler':QS(k0s, "CHEBYSHEV", rescaling, rescalingInv)}
solver = HSBPE(kappa = 12 ** .5, theta = np.pi / 3, n = 40,
verbosity = verb)
solver.omega = np.real(k0)
approx = RB(solver, mu0 = k0, approxParameters = params,
verbosity = verb)
approx.setupApprox()
# approx.plotSamples()
approx.plotApp(ktar, name = 'u_RB')
approx.plotHF(ktar, name = 'u_HF')
approx.plotErr(ktar, name = 'err')
approx.plotRes(ktar, name = 'res')
appErr, solNorm = approx.normErr(ktar), approx.normHF(ktar)
resNorm, RHSNorm = approx.normRes(ktar), approx.normRHS(ktar)
print(('SolNorm:\t{}\nErr:\t{}\nErrRel:\t{}').format(solNorm, appErr,
np.divide(appErr, solNorm)))
print(('RHSNorm:\t{}\nRes:\t{}\nResRel:\t{}').format(RHSNorm, resNorm,
np.divide(resNorm, RHSNorm)))
############
elif testNo == 2:
k0s = [3.85 + 0.j, 4.15 + 0.j]
k0 = np.mean(k0s)
ktar = 4 + .15j
rescaling = lambda x: np.power(x, 2.)
rescalingInv = lambda x: np.power(x, .5)
params = {'S':10, 'R':9, 'POD':True,
'sampler':QS(k0s, "CHEBYSHEV", rescaling, rescalingInv)}
solver = HSTPE(nT = 2, nB = 1, theta = np.pi * 45/180, kappa = 4., n = 50,
verbosity = verb)
solver.omega = np.real(k0)
approx = RB(solver, mu0 = k0, approxParameters = params,
- verbosity = verb, homogeneize = homog)
+ verbosity = verb, homogeneized = homog)
approx.setupApprox()
# approx.plotSamples()
approx.plotApp(ktar, name = 'u_RB')
approx.plotHF(ktar, name = 'u_HF')
approx.plotErr(ktar, name = 'err')
approx.plotRes(ktar, name = 'res')
appErr, solNorm = approx.normErr(ktar), approx.normHF(ktar)
resNorm, RHSNorm = approx.normRes(ktar), approx.normRHS(ktar)
print(('SolNorm:\t{}\nErr:\t{}\nErrRel:\t{}').format(solNorm, appErr,
np.divide(appErr, solNorm)))
print(('RHSNorm:\t{}\nRes:\t{}\nResRel:\t{}').format(RHSNorm, resNorm,
np.divide(resNorm, RHSNorm)))
############
elif testNo == 3:
k0s = [2, 5]
k0 = np.mean(k0s)
ktar = 4.5 - 0.j
params = {'S':15, 'R':10, 'POD':True, 'sampler':QS(k0s, "CHEBYSHEV")}
solver = HBSPE(R = 7, kappa = 3, theta = - np.pi * 75 / 180, n = 40,
verbosity = verb)
solver.omega = np.real(k0)
approx = RB(solver, mu0 = k0, approxParameters = params,
- verbosity = verb, homogeneize = homog)
+ verbosity = verb, homogeneized = homog)
approx.setupApprox()
# approx.plotSamples()
approx.plotApp(ktar, name = 'u_RB')
approx.plotHF(ktar, name = 'u_HF')
approx.plotErr(ktar, name = 'err')
approx.plotRes(ktar, name = 'res')
appErr, solNorm = approx.normErr(ktar), approx.normHF(ktar)
resNorm, RHSNorm = approx.normRes(ktar), approx.normRHS(ktar)
print(('SolNorm:\t{}\nErr:\t{}\nErrRel:\t{}').format(solNorm, appErr,
np.divide(appErr, solNorm)))
print(('RHSNorm:\t{}\nRes:\t{}\nResRel:\t{}').format(RHSNorm, resNorm,
np.divide(resNorm, RHSNorm)))
diff --git a/examples/scipy/RBTaylor.py b/examples/pod/RBTaylor.py
similarity index 96%
rename from examples/scipy/RBTaylor.py
rename to examples/pod/RBTaylor.py
index b95241c..00b729a 100644
--- a/examples/scipy/RBTaylor.py
+++ b/examples/pod/RBTaylor.py
@@ -1,90 +1,90 @@
import numpy as np
from rrompy.hfengines.scipy import HelmholtzSquareBubbleProblemEngine as HSBPE
from rrompy.hfengines.scipy import HelmholtzSquareTransmissionProblemEngine as HSTPE
from rrompy.hfengines.scipy import HelmholtzBoxScatteringProblemEngine as HBSPE
from rrompy.reduction_methods.taylor import ApproximantTaylorRB as RB
testNo = 4
verb = 0
homog = True
#homog = False
if testNo == 1:
params = {'E':4, 'R':4, 'sampleType':'Arnoldi', 'POD':True}
k0 = 12 ** .5
ktar = 10.5 ** .5
solver = HSBPE(kappa = 12 ** .5, theta = np.pi / 3, n = 40,
verbosity = verb)
solver.omega = np.real(k0)
approx = RB(solver, mu0 = k0, approxParameters = params,
verbosity = verb)
approx.setupApprox()
# approx.plotSamples()
approx.plotApp(ktar, name = 'u_RB')
approx.plotHF(ktar, name = 'u_HF')
approx.plotErr(ktar, name = 'err')
approx.plotRes(ktar, name = 'res')
appErr, solNorm = approx.normErr(ktar), approx.normHF(ktar)
resNorm, RHSNorm = approx.normRes(ktar), approx.normRHS(ktar)
print(('SolNorm:\t{}\nErr:\t{}\nErrRel:\t{}').format(solNorm, appErr,
np.divide(appErr, solNorm)))
print(('RHSNorm:\t{}\nRes:\t{}\nResRel:\t{}').format(RHSNorm, resNorm,
np.divide(resNorm, RHSNorm)))
############
elif testNo == 2:
params = {'E':7, 'R':7, 'sampleType':'Arnoldi', 'POD':True}
k0 = 16**.5
ktar = 15**.5
solver = HSTPE(nT = 2, nB = 1, theta = np.pi * 45/180, kappa = 3.,
n = 50, verbosity = verb)
solver.omega = np.real(k0)
approx = RB(solver, mu0 = k0, approxParameters = params,
- verbosity = verb, homogeneize = homog)
+ verbosity = verb, homogeneized = homog)
approx.setupApprox()
# approx.plotSamples()
approx.plotApp(ktar, name = 'u_RB')
approx.plotHF(ktar, name = 'u_HF')
approx.plotErr(ktar, name = 'err')
approx.plotRes(ktar, name = 'res')
appErr, solNorm = approx.normErr(ktar), approx.normHF(ktar)
resNorm, RHSNorm = approx.normRes(ktar), approx.normRHS(ktar)
print(('SolNorm:\t{}\nErr:\t{}\nErrRel:\t{}').format(solNorm, appErr,
np.divide(appErr, solNorm)))
print(('RHSNorm:\t{}\nRes:\t{}\nResRel:\t{}').format(RHSNorm, resNorm,
np.divide(resNorm, RHSNorm)))
############
elif testNo in [3, 4]:
if testNo == 3:
params = {'E':8, 'sampleType':'Krylov', 'POD':True}
else:
params = {'E':8, 'sampleType':'Arnoldi', 'POD':True}
k0 = 3
ktar = 4.25+.5j
solver = HBSPE(R = 5, kappa = 3, theta = - np.pi * 75 / 180,
n = 30, verbosity = verb)
solver.omega = np.real(k0)
approx = RB(solver, mu0 = k0, approxParameters = params,
- verbosity = verb, homogeneize = homog)
+ verbosity = verb, homogeneized = homog)
approx.setupApprox()
# approx.plotSamples()
approx.plotApp(ktar, name = 'u_RB')
approx.plotHF(ktar, name = 'u_HF')
approx.plotErr(ktar, name = 'err')
approx.plotRes(ktar, name = 'res')
appErr, solNorm = approx.normErr(ktar), approx.normHF(ktar)
resNorm, RHSNorm = approx.normRes(ktar), approx.normRHS(ktar)
print(('SolNorm:\t{}\nErr:\t{}\nErrRel:\t{}').format(solNorm, appErr,
np.divide(appErr, solNorm)))
print(('RHSNorm:\t{}\nRes:\t{}\nResRel:\t{}').format(RHSNorm, resNorm,
np.divide(resNorm, RHSNorm)))
diff --git a/examples/scipy/TaylorPoles.py b/examples/pod/TaylorPoles.py
similarity index 100%
rename from examples/scipy/TaylorPoles.py
rename to examples/pod/TaylorPoles.py
diff --git a/examples/scipy/TaylorSweep.py b/examples/pod/TaylorSweep.py
similarity index 94%
rename from examples/scipy/TaylorSweep.py
rename to examples/pod/TaylorSweep.py
index 9e9c8a3..3cf37e6 100644
--- a/examples/scipy/TaylorSweep.py
+++ b/examples/pod/TaylorSweep.py
@@ -1,77 +1,77 @@
import numpy as np
from rrompy.hfengines.scipy import HelmholtzSquareBubbleProblemEngine as HSBPE
from rrompy.hfengines.scipy import HelmholtzBoxScatteringProblemEngine as HBSPE
from rrompy.reduction_methods.taylor import ApproximantTaylorPade as Pade
from rrompy.reduction_methods.taylor import ApproximantTaylorRB as RB
from rrompy.utilities.parameter_sweeper import ParameterSweeper as Sweeper
testNo = 1
verb = 0
homog = True
homog = False
k0 = (12 + 0.j) ** .5
npoints = 100
shift = 5
nsets = 2
stride = 2
Emax = stride * (nsets - 1) + shift + 1
if testNo == 1:
solver = HSBPE(kappa = 12 ** .5, theta = np.pi / 3, n = 15,
verbosity = verb)
solver.omega = np.real(k0)
params = {'Emax':Emax, 'sampleType':'ARNOLDI', 'POD':True}
ktars = np.linspace(7**.5, 16**.5, npoints)
filenamebase = '../data/output/HelmholtzBubbleTSweep'
homog = False
elif testNo == 2:
solver = HBSPE(R = 5, kappa = 3, theta = - np.pi * 75 / 180, n = 10,
verbosity = verb)
solver.omega = np.real(k0)
params = {'Emax':Emax, 'sampleType':'Arnoldi', 'POD':True}
ktars = np.linspace(7**.5, 17**.5, npoints)
homogMSG = ""
if not homog: homogMSG = "Non"
filenamebase = '../data/output/HelmholtzBoxTSweep' + homogMSG + "Homog"
paramsSetsPade = [None] * nsets
paramsSetsRB = [None] * nsets
paramsSetsPoly = [None] * nsets
for i in range(nsets):
paramsSetsPade[i] = {'N':stride*i+shift+1, 'M':stride*i+shift+1,
'E':stride*i+shift+1}
paramsSetsRB[i] = {'E':stride*i+shift+1,'R':stride*i+shift+2}
paramsSetsPoly[i] = {'N':0, 'M':stride*i+shift+1,
'E':stride*i+shift+1}
appPade = Pade(solver, mu0 = k0, approxParameters = params,
- verbosity = verb, homogeneize = homog)
+ verbosity = verb, homogeneized = homog)
appRB = RB(solver, mu0 = k0, approxParameters = params,
- verbosity = verb, homogeneize = homog)
+ verbosity = verb, homogeneized = homog)
appPoly = Pade(solver, mu0 = k0, approxParameters = params,
- verbosity = verb, homogeneize = homog)
+ verbosity = verb, homogeneized = homog)
sweeper = Sweeper(mutars = ktars, mostExpensive = 'Approx')
sweeper.ROMEngine = appPade
sweeper.params = paramsSetsPade
filenamePade = sweeper.sweep(filenamebase + 'Pade.dat')
sweeper.ROMEngine = appRB
sweeper.params = paramsSetsRB
filenameRB = sweeper.sweep(filenamebase + 'RB.dat')
sweeper.ROMEngine = appPoly
sweeper.params = paramsSetsPoly
filenamePoly = sweeper.sweep(filenamebase + 'Poly.dat')
sweeper.plotCompare([filenamePade, filenameRB, filenamePoly], ['muRe'],
['normHF', 'normApp'], ['E'], onePlot = True,
save = filenamebase + 'Norm',
saveFormat = "png", labels = ["Pade'", "RB", "Poly"])
sweeper.plotCompare([filenamePade, filenameRB, filenamePoly], ['muRe'],
['normResRel'], ['E'], save = filenamebase + 'Res',
saveFormat = "png", labels = ["Pade'", "RB", "Poly"])
sweeper.plotCompare([filenamePade, filenameRB, filenamePoly], ['muRe'],
['normErrRel'], ['E'], save = filenamebase + 'Err',
saveFormat = "png", labels = ["Pade'", "RB", "Poly"])
diff --git a/examples/scipy/airfoil.py b/examples/pod/airfoil.py
similarity index 94%
rename from examples/scipy/airfoil.py
rename to examples/pod/airfoil.py
index 498aa50..2121fff 100644
--- a/examples/scipy/airfoil.py
+++ b/examples/pod/airfoil.py
@@ -1,212 +1,212 @@
from copy import deepcopy as copy
import numpy as np
import fenics as fen
import ufl
from rrompy.hfengines.scipy import HelmholtzBoxScatteringProblemEngine as HSP
from rrompy.reduction_methods.taylor import ApproximantTaylorPade as TP
from rrompy.reduction_methods.taylor import ApproximantTaylorRB as TRB
from rrompy.reduction_methods.lagrange import ApproximantLagrangePade as LP
from rrompy.reduction_methods.lagrange import ApproximantLagrangeRB as LRB
from rrompy.utilities.parameter_sweeper import ParameterSweeper as Sweeper
from rrompy.utilities.parameter_sampling import QuadratureSampler as QS
from rrompy.utilities.base.fenics import fenONE
from operator import itemgetter
def subdict(d, ks):
return dict(zip(ks, itemgetter(*ks)(d)))
verb = 0
####################
homog = True
#homog = False
####################
test = "solve"
test = "Taylor"
test = "Lagrange"
test = "TaylorSweep"
test = "LagrangeSweep"
plotSamples = True
k0 = 10
kLeft, kRight = 8 + 0.j, 12 + 0.j
ktar = 11
ktars = np.linspace(8, 12, 21) + 0.j
PI = np.pi
R = 2
def Dboundary(x, on_boundary):
return on_boundary and (x[0]**2+x[1]**2)**.5 < .95 * R
kappa = 10
theta = PI * - 45 / 180.
mu = 1.1
epsilon = .1
mesh = fen.Mesh('../data/mesh/airfoil.xml')
c, s = np.cos(theta), np.sin(theta)
x, y = fen.SpatialCoordinate(mesh)[:]
u0R = - fen.cos(kappa * (c * x + s * y))
u0I = - fen.sin(kappa * (c * x + s * y))
checkReal = x**2-x+y**2
rhop5 = ((x**2+y**2)/((x-1)**2+y**2))**.25
phiroot1 = fen.atan(-y/(x**2-x+y**2)) / 2
phiroot2 = fen.atan(-y/(x**2-x+y**2)) / 2 - PI * ufl.sign(-y/(x**2-x+y**2)) / 2
kappam1 = (((rhop5*fen.cos(phiroot1)+.5)**2.+(rhop5*fen.sin(phiroot1))**2.)/
((rhop5*fen.cos(phiroot1)-1.)**2.+(rhop5*fen.sin(phiroot1))**2.)
)**.5 - mu
kappam2 = (((rhop5*fen.cos(phiroot2)+.5)**2.+(rhop5*fen.sin(phiroot2))**2.)/
((rhop5*fen.cos(phiroot2)-1.)**2.+(rhop5*fen.sin(phiroot2))**2.)
)**.5 - mu
Heps1 = .9 * .5 * (1 + kappam1/epsilon + fen.sin(PI*kappam1/epsilon) / PI) + .1
Heps2 = .9 * .5 * (1 + kappam2/epsilon + fen.sin(PI*kappam2/epsilon) / PI) + .1
cTT = ufl.conditional(ufl.le(kappam1, epsilon), Heps1, fenONE)
c_F = fen.Constant(.1)
cFT = ufl.conditional(ufl.le(kappam2, epsilon), Heps2, fenONE)
c_F = fen.Constant(.1)
cT = ufl.conditional(ufl.ge(kappam1, - epsilon), cTT, c_F)
cF = ufl.conditional(ufl.ge(kappam2, - epsilon), cFT, c_F)
a = ufl.conditional(ufl.ge(checkReal, 0.), cT, cF)
###
solver = HSP(R, np.abs(k0), theta, n = 1, verbosity = verb,
degree_threshold = 8)
solver.V = fen.FunctionSpace(mesh, "P", 3)
solver.diffusivity = a
solver.DirichletBoundary = Dboundary
solver.RobinBoundary = "REST"
solver.DirichletDatum = [u0R, u0I]
###
if test == "solve":
uinc = solver.liftDirichletData(k0)
if homog:
uhtot = solver.solve(k0, homogeneized = homog)
uh = uhtot + uinc
else:
uh = solver.solve(k0, homogeneized = homog)
uhtot = uh - uinc
print(solver.norm(uh))
print(solver.norm(uhtot))
solver.plot(fen.project(a, solver.V).vector(), what = 'Real',
name = 'a')
solver.plot(uinc, what = 'Real', name = 'u_inc')
solver.plot(uh, what = 'ABS')
solver.plot(uhtot, what = 'ABS', name = 'u + u_inc')
elif test in ["Taylor", "Lagrange"]:
if test == "Taylor":
params = {'N':8, 'M':8, 'R':8, 'E':8,
'sampleType':'Arnoldi', 'POD':True}
parPade = subdict(params, ['N', 'M', 'E', 'sampleType', 'POD'])
parRB = subdict(params, ['R', 'E', 'sampleType', 'POD'])
approxPade = TP(solver, mu0 = k0, approxParameters = parPade,
- verbosity = verb, homogeneize = homog)
+ verbosity = verb, homogeneized = homog)
approxRB = TRB(solver, mu0 = k0, approxParameters = parRB,
- verbosity = verb, homogeneize = homog)
+ verbosity = verb, homogeneized = homog)
else:
params = {'N':8, 'M':8, 'R':9, 'S':9, 'POD':True,
'sampler':QS([kLeft, kRight], "CHEBYSHEV")}
parPade = subdict(params, ['N', 'M', 'S', 'POD', 'sampler'])
parRB = subdict(params, ['R', 'S', 'POD', 'sampler'])
approxPade = LP(solver, mu0 = np.mean([kLeft, kRight]),
approxParameters = parPade, verbosity = verb,
- homogeneize = homog)
+ homogeneized = homog)
approxRB = LRB(solver, mu0 = np.mean([kLeft, kRight]),
approxParameters = parRB, verbosity = verb,
- homogeneize = homog)
+ homogeneized = homog)
approxPade.setupApprox()
approxRB.setupApprox()
if plotSamples:
approxPade.plotSamples()
approxPade.plotHF(ktar, name = 'u_HF')
approxPade.plotApp(ktar, name = 'u_Pade''')
approxPade.plotErr(ktar, name = 'err_Pade''')
approxPade.plotRes(ktar, name = 'res_Pade''')
approxRB.plotApp(ktar, name = 'u_RB')
approxRB.plotErr(ktar, name = 'err_RB')
approxRB.plotRes(ktar, name = 'res_RB')
HFNorm, RHSNorm = approxPade.normHF(ktar), approxPade.normRHS(ktar)
PadeRes, PadeErr = approxPade.normRes(ktar), approxPade.normErr(ktar)
RBRes, RBErr = approxRB.normRes(ktar), approxRB.normErr(ktar)
print('HFNorm:\t{}\nRHSNorm:\t{}'.format(HFNorm, RHSNorm))
print('PadeRes:\t{}\nPadeErr:\t{}'.format(PadeRes, PadeErr))
print('RBRes:\t{}\nRBErr:\t{}'.format(RBRes, RBErr))
print('\nPoles Pade'':')
print(approxPade.getPoles())
elif test in ["TaylorSweep", "LagrangeSweep"]:
if test == "TaylorSweep":
shift = 10
nsets = 2
stride = 3
Emax = stride * (nsets - 1) + shift + 1
params = {'Emax':Emax, 'sampleType':'Arnoldi', 'POD':True}
paramsSetsPade = [None] * nsets
paramsSetsRB = [None] * nsets
for i in range(nsets):
paramsSetsPade[i] = {'N':stride*i+shift + 1, 'M':stride*i+shift+1,
'E':stride*i+shift + 1}
paramsSetsRB[i] = {'E':stride*i+shift + 1,'R':stride*i+shift + 2}
approxPade = TP(solver, mu0 = kappa,approxParameters = params,
- verbosity = verb, homogeneize = homog)
+ verbosity = verb, homogeneized = homog)
approxRB = TRB(solver, mu0 = kappa, approxParameters = params,
- verbosity = verb, homogeneize = homog)
+ verbosity = verb, homogeneized = homog)
else:
shift = 10
nsets = 2
stride = 3
Smax = stride * (nsets - 1) + shift + 2
paramsPade = {'S':Smax, 'POD':True,
'sampler':QS([kLeft, kRight], "CHEBYSHEV")}
paramsRB = copy(paramsPade)
paramsSetsPade = [None] * nsets
paramsSetsRB = [None] * nsets
for i in range(nsets):
paramsSetsPade[i] = {'N': stride*i+shift+1, 'M': stride*i+shift+1,
'S': stride*i+shift+2}
paramsSetsRB[i] = {'R': stride*i+shift+2, 'S': stride*i+shift+2}
approxPade = LP(solver, mu0 = np.mean([kLeft, kRight]),
approxParameters = paramsPade, verbosity = verb,
- homogeneize = homog)
+ homogeneized = homog)
approxRB = LRB(solver, mu0 = np.mean([kLeft, kRight]),
approxParameters = paramsRB, verbosity = verb,
- homogeneize = homog)
+ homogeneized = homog)
sweeper = Sweeper(mutars = ktars, mostExpensive = 'Approx')
homogMSG = ""
if not homog: homogMSG = "Non"
filenamebase = '../data/output/airfoil' + test[:-5] + homogMSG + "Homog"
sweeper.ROMEngine = approxPade
sweeper.params = paramsSetsPade
filenamePade = sweeper.sweep(filenamebase +'Pade.dat')
sweeper.ROMEngine = approxRB
sweeper.params = paramsSetsRB
filenameRB = sweeper.sweep(filenamebase +'RB.dat')
if test == "TaylorSweep":
constr = ['E']
else:
constr = ['S']
sweeper.plotCompare([filenamePade, filenameRB], ['muRe'],
['normHF', 'normApp'], constr, onePlot = True,
save = filenamebase + 'Norm',
saveFormat = "png", labels = ["Pade'", "RB"])
sweeper.plotCompare([filenamePade, filenameRB], ['muRe'], ['normResRel'],
constr, save = filenamebase + 'Res',
saveFormat = "png", labels = ["Pade'", "RB"])
sweeper.plotCompare([filenamePade, filenameRB], ['muRe'], ['normErrRel'],
constr, save = filenamebase + 'Err',
saveFormat = "png", labels = ["Pade'", "RB"])
diff --git a/examples/scipy/laplaceGaussianTaylor.py b/examples/pod/laplaceGaussianTaylor.py
similarity index 100%
rename from examples/scipy/laplaceGaussianTaylor.py
rename to examples/pod/laplaceGaussianTaylor.py
diff --git a/examples/scipy/membraneTaylor.py b/examples/pod/membraneTaylor.py
similarity index 100%
rename from examples/scipy/membraneTaylor.py
rename to examples/pod/membraneTaylor.py
diff --git a/examples/scipy/parametricDomain.py b/examples/pod/parametricDomain.py
similarity index 100%
rename from examples/scipy/parametricDomain.py
rename to examples/pod/parametricDomain.py
diff --git a/examples/scipy/scatteringSquare.py b/examples/pod/scatteringSquare.py
similarity index 100%
rename from examples/scipy/scatteringSquare.py
rename to examples/pod/scatteringSquare.py
diff --git a/rrompy/reduction_methods/base/__init__.py b/rrompy/reduction_methods/base/__init__.py
index fdc4ada..993f983 100644
--- a/rrompy/reduction_methods/base/__init__.py
+++ b/rrompy/reduction_methods/base/__init__.py
@@ -1,25 +1,28 @@
# Copyright (C) 2018 by the RROMPy authors
#
# This file is part of RROMPy.
#
# RROMPy is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# RROMPy is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with RROMPy. If not, see .
#
from .generic_approximant import GenericApproximant
+from .pade_utils import checkRobustTolerance, checkPolyfitRank
__all__ = [
- 'GenericApproximant'
+ 'GenericApproximant',
+ 'checkRobustTolerance',
+ 'checkPolyfitRank'
]
diff --git a/rrompy/reduction_methods/base/generic_approximant.py b/rrompy/reduction_methods/base/generic_approximant.py
index 804e5f6..9b25cae 100644
--- a/rrompy/reduction_methods/base/generic_approximant.py
+++ b/rrompy/reduction_methods/base/generic_approximant.py
@@ -1,530 +1,530 @@
# Copyright (C) 2018 by the RROMPy authors
#
# This file is part of RROMPy.
#
# RROMPy is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# RROMPy is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with RROMPy. If not, see .
#
from abc import abstractmethod
import numpy as np
from copy import copy
from rrompy.sampling.base.sampling_engine_base import SamplingEngineBase
from rrompy.utilities.base.types import Np1D, DictAny, HFEng, sampleEng, strLst
from rrompy.utilities.base import purgeDict, verbosityDepth
__all__ = ['GenericApproximant']
class GenericApproximant:
"""
ABSTRACT
ROM approximant computation for parametric problems.
Args:
HFEngine: HF problem solver.
mu0(optional): Default parameter. Defaults to 0.
approxParameters(optional): Dictionary containing values for main
parameters of approximant. Recognized keys are:
- 'POD': whether to compute POD of snapshots; defaults to True.
Defaults to empty dict.
verbosity(optional): Verbosity level. Defaults to 10.
Attributes:
HFEngine: HF problem solver.
mu0: Default parameter.
w: Weight for norm computation.
approxParameters: Dictionary containing values for main parameters of
approximant. Recognized keys are in parameterList.
parameterList: Recognized keys of approximant parameters:
- 'POD': whether to compute POD of snapshots.
verbosity: Verbosity level.
POD: Whether to compute POD of snapshots.
samplingEngine: Sampling engine.
uHF: High fidelity solution with wavenumber lastSolvedHF as numpy
complex vector.
lastSolvedHF: Wavenumber corresponding to last computed high fidelity
solution.
uApp: Last evaluated approximant as numpy complex vector.
lastApproxParameters: List of parameters corresponding to last
computed approximant.
"""
def __init__(self, HFEngine:HFEng, mu0 : complex = 0,
- approxParameters : DictAny = {}, homogeneize : bool = False,
+ approxParameters : DictAny = {}, homogeneized : bool = False,
verbosity : int = 10):
self._preInit()
self.verbosity = verbosity
if self.verbosity >= 10:
verbosityDepth("INIT", ("Initializing approximant engine of "
"type {}.").format(self.name()))
self.HFEngine = HFEngine
self._HFEngine0 = copy(HFEngine)
if not hasattr(self, "parameterList"):
self.parameterList = []
self.parameterList += ["POD"]
self.mu0 = mu0
- self.homogeneize = homogeneize
+ self.homogeneized = homogeneized
self.approxParameters = approxParameters
self._postInit()
def _preInit(self):
if not hasattr(self, "depth"): self.depth = 0
else: self.depth += 1
def _postInit(self):
if self.depth == 0:
if self.verbosity >= 10:
verbosityDepth("DEL", "Done initializing.\n")
del self.depth
else: self.depth -= 1
def name(self) -> str:
return self.__class__.__name__
def __str__(self) -> str:
return self.name()
def __repr__(self) -> str:
return self.__str__() + " at " + hex(id(self))
def setupSampling(self, SamplingEngine : sampleEng = SamplingEngineBase):
"""Setup sampling engine."""
self.samplingEngine = SamplingEngine(self.HFEngine,
verbosity = self.verbosity)
@property
def approxParameters(self):
"""Value of approximant parameters."""
return self._approxParameters
@approxParameters.setter
def approxParameters(self, approxParams):
if not hasattr(self, "approxParameters"):
self._approxParameters = {}
approxParameters = purgeDict(approxParams, self.parameterList,
dictname = self.name() + ".approxParameters",
baselevel = 1)
keyList = list(approxParameters.keys())
if "POD" in keyList:
self.POD = approxParameters["POD"]
elif hasattr(self, "POD"):
self.POD = self.POD
else:
self.POD = True
@property
def POD(self):
"""Value of POD."""
return self._POD
@POD.setter
def POD(self, POD):
if hasattr(self, "POD"): PODold = self.POD
else: PODold = -1
self._POD = POD
self._approxParameters["POD"] = self.POD
if PODold != self.POD:
self.setupSampling()
self.resetSamples()
@property
- def homogeneize(self):
- """Value of homogeneize."""
- return self._homogeneize
- @homogeneize.setter
- def homogeneize(self, homogeneize):
- if not hasattr(self, "_homogeneize"):
- self._homogeneize = None
- if homogeneize != self.homogeneize:
- self._homogeneize = homogeneize
+ def homogeneized(self):
+ """Value of homogeneized."""
+ return self._homogeneized
+ @homogeneized.setter
+ def homogeneized(self, homogeneized):
+ if not hasattr(self, "_homogeneized"):
+ self._homogeneized = None
+ if homogeneized != self.homogeneized:
+ self._homogeneized = homogeneized
self.resetSamples()
def solveHF(self, mu : complex = None):
"""
Find high fidelity solution with original parameters and arbitrary
parameter.
Args:
mu: Target parameter.
"""
if mu is None: mu = self.mu0
if (not hasattr(self, "lastSolvedHF")
or not np.isclose(self.lastSolvedHF, mu)):
self.uHF = self.samplingEngine.solveLS(mu,
- homogeneized = self.homogeneize)
+ homogeneized = self.homogeneized)
self.lastSolvedHF = mu
def resetSamples(self):
"""Reset samples."""
if hasattr(self, "samplingEngine"):
self.samplingEngine.resetHistory()
else:
self.setupSampling()
def plotSamples(self, name : str = "u", save : str = None,
what : strLst = 'all', saveFormat : str = "eps",
saveDPI : int = 100, **figspecs):
"""
Do some nice plots of the samples.
Args:
name(optional): Name to be shown as title of the plots. Defaults to
'u'.
what(optional): Which plots to do. If list, can contain 'ABS',
'PHASE', 'REAL', 'IMAG'. If str, same plus wildcard 'ALL'.
Defaults to 'ALL'.
save(optional): Where to save plot(s). Defaults to None, i.e. no
saving.
saveFormat(optional): Format for saved plot(s). Defaults to "eps".
saveDPI(optional): DPI for saved plot(s). Defaults to 100.
figspecs(optional key args): Optional arguments for matplotlib
figure creation.
"""
self.samplingEngine.plotSamples(name = name, save = save, what = what,
saveFormat = saveFormat,
saveDPI = saveDPI,
**figspecs)
@abstractmethod
def setupApprox(self):
"""
Setup approximant. (ABSTRACT)
Any specialization should include something like
self.computeDerivatives()
if not self.checkComputedApprox():
...
self.lastApproxParameters = copy(self.approxParameters)
"""
pass
def checkComputedApprox(self) -> bool:
"""
Check if setup of new approximant is not needed.
Returns:
True if new setup is not needed. False otherwise.
"""
return (hasattr(self, "lastApproxParameters")
and self.approxParameters == self.lastApproxParameters)
@abstractmethod
def evalApproxReduced(self, mu:complex):
"""
Evaluate reduced representation of approximant at arbitrary parameter.
(ABSTRACT)
Any specialization should include something like
self.setupApprox()
self.uAppReduced = ...
Args:
mu: Target parameter.
"""
pass
@abstractmethod
def evalApprox(self, mu:complex):
"""
Evaluate approximant at arbitrary parameter. (ABSTRACT)
Any specialization should include something like
self.evalApproxReduced(mu)
self.uApp = ...
Args:
mu: Target parameter.
"""
pass
def getHF(self, mu:complex, homogeneized : bool = False) -> Np1D:
"""
Get HF solution at arbitrary parameter.
Args:
mu: Target parameter.
homogeneized(optional): Whether to remove Dirichlet BC. Defaults to
False.
Returns:
HFsolution.
"""
self.solveHF(mu)
- if self.homogeneize and not homogeneized:
+ if self.homogeneized and not homogeneized:
return self.uHF + self.HFEngine.liftDirichletData(mu)
- if not self.homogeneize and homogeneized:
+ if not self.homogeneized and homogeneized:
return self.uHF - self.HFEngine.liftDirichletData(mu)
return self.uHF
def getRHS(self, mu:complex, homogeneized : bool = False) -> Np1D:
"""
Get linear system RHS at arbitrary parameter.
Args:
mu: Target parameter.
homogeneized(optional): Whether to remove Dirichlet BC. Defaults to
False.
Returns:
Linear system RHS.
"""
return self.HFEngine.residual(None, mu, homogeneized = homogeneized)
def getAppReduced(self, mu:complex) -> Np1D:
"""
Get approximant at arbitrary parameter.
Args:
mu: Target parameter.
Returns:
Reduced approximant.
"""
self.evalApproxReduced(mu)
return self.uAppReduced
def getApp(self, mu:complex, homogeneized : bool = False) -> Np1D:
"""
Get approximant at arbitrary parameter.
Args:
mu: Target parameter.
homogeneized(optional): Whether to remove Dirichlet BC. Defaults to
False.
Returns:
Approximant.
"""
self.evalApprox(mu)
- if self.homogeneize and not homogeneized:
+ if self.homogeneized and not homogeneized:
return self.uApp + self.HFEngine.liftDirichletData(mu)
- if not self.homogeneize and homogeneized:
+ if not self.homogeneized and homogeneized:
return self.uApp - self.HFEngine.liftDirichletData(mu)
return self.uApp
def getRes(self, mu:complex, homogeneized : bool = False) -> Np1D:
"""
Get residual at arbitrary parameter.
Args:
mu: Target parameter.
homogeneized(optional): Whether to remove Dirichlet BC. Defaults to
False.
Returns:
Approximant residual.
"""
return self.HFEngine.residual(self.getApp(mu, homogeneized), mu,
homogeneized = homogeneized)
def getErr(self, mu:complex, homogeneized : bool = False) -> Np1D:
"""
Get error at arbitrary parameter.
Args:
mu: Target parameter.
homogeneized(optional): Whether to remove Dirichlet BC. Defaults to
False.
Returns:
Approximant error.
"""
return self.getApp(mu, homogeneized) - self.getHF(mu, homogeneized)
@abstractmethod
def getPoles(self) -> Np1D:
"""
Obtain approximant poles.
Returns:
Numpy complex vector of poles.
"""
pass
def normHF(self, mu:complex, homogeneized : bool = False) -> float:
"""
Compute norm of HF solution at arbitrary parameter.
Args:
mu: Target parameter.
homogeneized(optional): Whether to remove Dirichlet BC. Defaults to
False.
Returns:
Target norm of HFsolution.
"""
return self.HFEngine.norm(self.getHF(mu, homogeneized))
def normRHS(self, mu:complex, homogeneized : bool = False) -> Np1D:
"""
Compute norm of linear system RHS at arbitrary parameter.
Args:
mu: Target parameter.
homogeneized(optional): Whether to remove Dirichlet BC. Defaults to
False.
Returns:
Norm of linear system RHS.
"""
return self.HFEngine.norm(self.getRHS(mu, homogeneized))
def normApp(self, mu:complex, homogeneized : bool = False) -> float:
"""
Compute norm of approximant at arbitrary parameter.
Args:
mu: Target parameter.
homogeneized(optional): Whether to remove Dirichlet BC. Defaults to
False.
Returns:
Target norm of approximant.
"""
return self.HFEngine.norm(self.getApp(mu, homogeneized))
def normRes(self, mu:complex, homogeneized : bool = False) -> float:
"""
Compute norm of approximant residual at arbitrary parameter.
Args:
mu: Target parameter.
homogeneized(optional): Whether to remove Dirichlet BC. Defaults to
False.
Returns:
Target norm of (A(mu)app(mu) - f(mu)).
"""
return self.HFEngine.norm(self.getRes(mu, homogeneized))
def normErr(self, mu:complex, homogeneized : bool = False) -> float:
"""
Compute norm of approximant error at arbitrary parameter.
Args:
mu: Target parameter.
homogeneized(optional): Whether to remove Dirichlet BC. Defaults to
False.
Returns:
Target norm of (approximant - HFsolution).
"""
return self.HFEngine.norm(self.getErr(mu, homogeneized))
def plotHF(self, mu:complex, name : str = "uHF", save : str = None,
what : strLst = 'all', saveFormat : str = "eps",
saveDPI : int = 100, homogeneized : bool = False, **figspecs):
"""
Do some nice plots of the HF solution at arbitrary parameter.
Args:
mu: Target parameter.
name(optional): Name to be shown as title of the plots. Defaults to
'u'.
what(optional): Which plots to do. If list, can contain 'ABS',
'PHASE', 'REAL', 'IMAG'. If str, same plus wildcard 'ALL'.
Defaults to 'ALL'.
save(optional): Where to save plot(s). Defaults to None, i.e. no
saving.
saveFormat(optional): Format for saved plot(s). Defaults to "eps".
saveDPI(optional): DPI for saved plot(s). Defaults to 100.
homogeneized(optional): Whether to remove Dirichlet BC. Defaults to
False.
figspecs(optional key args): Optional arguments for matplotlib
figure creation.
"""
uHF = self.getHF(mu, homogeneized)
self.HFEngine.plot(uHF, name = name, save = save, what = what,
saveFormat = saveFormat, saveDPI = saveDPI,
**figspecs)
def plotApp(self, mu:complex, name : str = "uApp", save : str = None,
what : strLst = 'all', saveFormat : str = "eps",
saveDPI : int = 100, homogeneized : bool = False, **figspecs):
"""
Do some nice plots of approximant at arbitrary parameter.
Args:
mu: Target parameter.
name(optional): Name to be shown as title of the plots. Defaults to
'u'.
what(optional): Which plots to do. If list, can contain 'ABS',
'PHASE', 'REAL', 'IMAG'. If str, same plus wildcard 'ALL'.
Defaults to 'ALL'.
save(optional): Where to save plot(s). Defaults to None, i.e. no
saving.
saveFormat(optional): Format for saved plot(s). Defaults to "eps".
saveDPI(optional): DPI for saved plot(s). Defaults to 100.
homogeneized(optional): Whether to remove Dirichlet BC. Defaults to
False.
figspecs(optional key args): Optional arguments for matplotlib
figure creation.
"""
uApp = self.getApp(mu, homogeneized)
self.HFEngine.plot(uApp, name = name, save = save, what = what,
saveFormat = saveFormat, saveDPI = saveDPI,
**figspecs)
def plotRes(self, mu:complex, name : str = "res", save : str = None,
what : strLst = 'all', saveFormat : str = "eps",
saveDPI : int = 100, homogeneized : bool = False, **figspecs):
"""
Do some nice plots of approximation residual at arbitrary parameter.
Args:
mu: Target parameter.
name(optional): Name to be shown as title of the plots. Defaults to
'u'.
what(optional): Which plots to do. If list, can contain 'ABS',
'PHASE', 'REAL', 'IMAG'. If str, same plus wildcard 'ALL'.
Defaults to 'ALL'.
save(optional): Where to save plot(s). Defaults to None, i.e. no
saving.
saveFormat(optional): Format for saved plot(s). Defaults to "eps".
saveDPI(optional): DPI for saved plot(s). Defaults to 100.
homogeneized(optional): Whether to remove Dirichlet BC. Defaults to
False.
figspecs(optional key args): Optional arguments for matplotlib
figure creation.
"""
uRes = self.getRes(mu, homogeneized)
self.HFEngine.plot(uRes, name = name, save = save, what = what,
saveFormat = saveFormat, saveDPI = saveDPI,
**figspecs)
def plotErr(self, mu:complex, name : str = "err", save : str = None,
what : strLst = 'all', saveFormat : str = "eps",
saveDPI : int = 100, homogeneized : bool = False, **figspecs):
"""
Do some nice plots of approximation error at arbitrary parameter.
Args:
mu: Target parameter.
name(optional): Name to be shown as title of the plots. Defaults to
'u'.
what(optional): Which plots to do. If list, can contain 'ABS',
'PHASE', 'REAL', 'IMAG'. If str, same plus wildcard 'ALL'.
Defaults to 'ALL'.
save(optional): Where to save plot(s). Defaults to None, i.e. no
saving.
saveFormat(optional): Format for saved plot(s). Defaults to "eps".
saveDPI(optional): DPI for saved plot(s). Defaults to 100.
homogeneized(optional): Whether to remove Dirichlet BC. Defaults to
False.
figspecs(optional key args): Optional arguments for matplotlib
figure creation.
"""
uErr = self.getErr(mu, homogeneized)
self.HFEngine.plot(uErr, name = name, save = save, what = what,
saveFormat = saveFormat, saveDPI = saveDPI,
**figspecs)
diff --git a/rrompy/reduction_methods/base/pade_utils.py b/rrompy/reduction_methods/base/pade_utils.py
new file mode 100644
index 0000000..87744e9
--- /dev/null
+++ b/rrompy/reduction_methods/base/pade_utils.py
@@ -0,0 +1,46 @@
+# Copyright (C) 2018 by the RROMPy authors
+#
+# This file is part of RROMPy.
+#
+# RROMPy is free software: you can redistribute it and/or modify
+# it under the terms of the GNU Lesser General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# RROMPy is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public License
+# along with RROMPy. If not, see .
+#
+
+import numpy as np
+from rrompy.utilities.base.types import Np1D, Np2D, Tuple
+from rrompy.utilities.warning_manager import warn
+
+machineEps = np.finfo(float).eps
+
+def checkRobustTolerance(ev:Np1D, E:int, tol:float) -> dict:
+ N = len(ev) - 1
+ ts = tol * np.linalg.norm(ev)
+ diff = N - np.sum(np.abs(ev) >= ts)
+ if diff <= 0: return {}
+ Enew = E - diff
+ Nnew = min(N, Enew)
+ if Nnew == N:
+ strN = ""
+ else:
+ strN = "N from {} to {} and ".format(N, Nnew)
+ warn(("Smallest {} eigenvalues below tolerance.\nReducing {}E from {} to "
+ "{}.").format(diff + 1, strN, E, Enew))
+ newParameters = {"N" : Nnew, "E" : Enew}
+ return newParameters
+
+def checkPolyfitRank(xs:Np1D, ys:Np2D, E:int, rcond:float,
+ w : Np1D = None) -> Tuple[bool, Np2D]:
+ G = np.polyfit(xs, ys, deg = E, w = w, rcond = rcond, full = True)
+ if G[2] < E + 1: return False, G[2] - 1
+ return True, G[0]
+
diff --git a/rrompy/reduction_methods/lagrange/approximant_lagrange_pade.py b/rrompy/reduction_methods/lagrange/approximant_lagrange_pade.py
index 132a1d5..f77b827 100644
--- a/rrompy/reduction_methods/lagrange/approximant_lagrange_pade.py
+++ b/rrompy/reduction_methods/lagrange/approximant_lagrange_pade.py
@@ -1,425 +1,455 @@
# Copyright (C) 2018 by the RROMPy authors
#
# This file is part of RROMPy.
#
# RROMPy is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# RROMPy is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with RROMPy. If not, see .
#
from copy import copy
import numpy as np
+from rrompy.reduction_methods.base import (checkRobustTolerance,
+ checkPolyfitRank)
from .generic_approximant_lagrange import GenericApproximantLagrange
from rrompy.utilities.base.types import Np1D, DictAny, List, HFEng
from rrompy.utilities.base import purgeDict, verbosityDepth
from rrompy.utilities.warning_manager import warn
__all__ = ['ApproximantLagrangePade']
class ApproximantLagrangePade(GenericApproximantLagrange):
"""
ROM Lagrange Pade' interpolant computation for parametric problems.
Args:
HFEngine: HF problem solver.
mu0(optional): Default parameter. Defaults to 0.
approxParameters(optional): Dictionary containing values for main
parameters of approximant. Recognized keys are:
- 'POD': whether to compute POD of snapshots; defaults to True;
- 'S': total number of samples current approximant relies upon;
defaults to 2;
- 'sampler': sample point generator; defaults to uniform sampler on
[0, 1];
- 'E': coefficient of interpolant to be minimized; defaults to
min(S, M + 1);
- 'M': degree of Pade' interpolant numerator; defaults to 0;
- - 'N': degree of Pade' interpolant denominator; defaults to 0.
+ - 'N': degree of Pade' interpolant denominator; defaults to 0;
+ - 'interpRcond': tolerance for interpolation via numpy.polyfit;
+ defaults to None;
+ - 'robustTol': tolerance for robust Pade' denominator management;
+ defaults to 0.
Defaults to empty dict.
Attributes:
HFEngine: HF problem solver.
mu0: Default parameter.
mus: Array of snapshot parameters.
ws: Array of snapshot weigths.
approxParameters: Dictionary containing values for main parameters of
approximant. Recognized keys are in parameterList.
parameterList: Recognized keys of approximant parameters:
- 'POD': whether to compute POD of snapshots;
- 'S': total number of samples current approximant relies upon;
- 'sampler': sample point generator;
- 'E': coefficient of interpolant to be minimized;
- 'M': degree of Pade' interpolant numerator;
- 'N': degree of Pade' interpolant denominator;
+ - 'interpRcond': tolerance for interpolation via numpy.polyfit;
- 'robustTol': tolerance for robust Pade' denominator management.
extraApproxParameters: List of approxParameters keys in addition to
mother class's.
S: Number of solution snapshots over which current approximant is
based upon.
sampler: Sample point generator.
M: Numerator degree of approximant.
N: Denominator degree of approximant.
POD: Whether to compute POD of snapshots.
+ interpRcond: Tolerance for interpolation via numpy.polyfit.
robustTol: Tolerance for robust Pade' denominator management.
samplingEngine: Sampling engine.
uHF: High fidelity solution with wavenumber lastSolvedHF as numpy
complex vector.
lastSolvedHF: Wavenumber corresponding to last computed high fidelity
solution.
Q: Numpy 1D vector containing complex coefficients of approximant
denominator.
P: Numpy 2D vector whose columns are FE dofs of coefficients of
approximant numerator.
uApp: Last evaluated approximant as numpy complex vector.
lastApproxParameters: List of parameters corresponding to last
computed approximant.
"""
def __init__(self, HFEngine:HFEng, mu0 : complex = 0.,
- approxParameters : DictAny = {}, homogeneize : bool = False,
+ approxParameters : DictAny = {}, homogeneized : bool = False,
verbosity : int = 10):
self._preInit()
if not hasattr(self, "parameterList"):
self.parameterList = []
- self.parameterList += ["E", "M", "N", "robustTol"]
+ self.parameterList += ["E", "M", "N", "interpRcond", "robustTol"]
super().__init__(HFEngine = HFEngine, mu0 = mu0,
approxParameters = approxParameters,
- homogeneize = homogeneize,
+ homogeneized = homogeneized,
verbosity = verbosity)
self._postInit()
@property
def approxParameters(self):
"""
- Value of approximant parameters. Its assignment may change E, M, N and
- S.
+ Value of approximant parameters. Its assignment may change E, M, N,
+ robustTol and S.
"""
return self._approxParameters
@approxParameters.setter
def approxParameters(self, approxParams):
approxParameters = purgeDict(approxParams, self.parameterList,
dictname = self.name() + ".approxParameters",
baselevel = 1)
- approxParametersCopy = purgeDict(approxParameters, ["E", "M", "N"],
+ approxParametersCopy = purgeDict(approxParameters, ["E", "M", "N",
+ "interpRcond",
+ "robustTol"],
True, True, baselevel = 1)
if hasattr(self, "M"):
Mold = self.M
self._M = 0
if hasattr(self, "N"):
Nold = self.N
self._N = 0
if hasattr(self, "E"):
self._E = 0
GenericApproximantLagrange.approxParameters.fset(self,
approxParametersCopy)
keyList = list(approxParameters.keys())
+ if "interpRcond" in keyList:
+ self.interpRcond = approxParameters["interpRcond"]
+ elif hasattr(self, "interpRcond"):
+ self.interpRcond = self.interpRcond
+ else:
+ self.interpRcond = None
if "robustTol" in keyList:
self.robustTol = approxParameters["robustTol"]
elif hasattr(self, "robustTol"):
self.robustTol = self.robustTol
else:
self.robustTol = 0
if "M" in keyList:
self.M = approxParameters["M"]
elif hasattr(self, "M"):
self.M = Mold
else:
self.M = 0
if "N" in keyList:
self.N = approxParameters["N"]
elif hasattr(self, "N"):
self.N = Nold
else:
self.N = 0
if "E" in keyList:
self.E = approxParameters["E"]
else:
self.E = min(self.S - 1, self.M + 1)
@property
def M(self):
"""Value of M. Its assignment may change S."""
return self._M
@M.setter
def M(self, M):
if M < 0: raise ArithmeticError("M must be non-negative.")
self._M = M
self._approxParameters["M"] = self.M
if hasattr(self, "S") and self.S < self.M + 1:
warn("Prescribed S is too small. Updating S to M + 1.")
self.S = self.M + 1
@property
def N(self):
"""Value of N. Its assignment may change S."""
return self._N
@N.setter
def N(self, N):
if N < 0: raise ArithmeticError("N must be non-negative.")
self._N = N
self._approxParameters["N"] = self.N
if hasattr(self, "S") and self.S < self.N + 1:
warn("Prescribed S is too small. Updating S to N + 1.")
self.S = self.N + 1
@property
def E(self):
"""Value of E. Its assignment may change S."""
return self._E
@E.setter
def E(self, E):
if E < 0: raise ArithmeticError("E must be non-negative.")
self._E = E
self._approxParameters["E"] = self.E
if hasattr(self, "S") and self.S < self.E + 1:
warn("Prescribed S is too small. Updating S to E + 1.")
self.S = self.E + 1
@property
def robustTol(self):
"""Value of tolerance for robust Pade' denominator management."""
return self._robustTol
@robustTol.setter
def robustTol(self, robustTol):
if robustTol < 0.:
warn("Overriding prescribed negative robustness tolerance to 0.")
robustTol = 0.
self._robustTol = robustTol
self._approxParameters["robustTol"] = self.robustTol
@property
def S(self):
"""Value of S."""
return self._S
@S.setter
def S(self, S):
if S <= 0: raise ArithmeticError("S must be positive.")
if hasattr(self, "S"): Sold = self.S
else: Sold = -1
vals, label = [0] * 3, {0:"M", 1:"N", 2:"E"}
if hasattr(self, "M"): vals[0] = self.M
if hasattr(self, "N"): vals[1] = self.N
if hasattr(self, "E"): vals[2] = self.E
idxmax = np.argmax(vals)
if vals[idxmax] + 1 > S:
warn("Prescribed S is too small. Updating S to {} + 1."\
.format(label[idxmax]))
self.S = vals[idxmax] + 1
else:
self._S = S
self._approxParameters["S"] = self.S
if Sold != self.S:
self.resetSamples()
def getTargetFunctionalOptimum(self):
"""Compute optimum of target LS functional."""
self.setupApprox()
- if self._funcOptimumWarn:
- warn("Target functional optimum numerically zero.")
- return np.abs(self._funcOptimum)
+ if self._funcOptimumIllCond:
+ if self.POD:
+ jOpt = np.linalg.norm(self.P[:, 0])
+ else:
+ jOpt = self.HFEngine.norm(self.samplingEngine.samples.dot(
+ self.P[:, 0]))
+ else:
+ jOpt = np.abs(self._funcOptimum)
+ return jOpt
def setupApprox(self):
"""
Compute Pade' interpolant.
SVD-based robust eigenvalue management.
"""
if not self.checkComputedApprox():
if self.verbosity >= 5:
verbosityDepth("INIT", "Setting up {}.". format(self.name()))
self.computeSnapshots()
if self.N > 0:
if self.verbosity >= 7:
verbosityDepth("INIT", ("Starting computation of "
"denominator."))
while self.N > 0:
TN = np.vander(self.radiusPade(self.mus), N = self.N + 1)
if self.POD:
data = self.samplingEngine.RPOD
else:
data = self.samplingEngine.samples
RHSFull = np.empty((self.S, data.shape[0] * (self.N + 1)),
dtype = np.complex)
for j in range(self.S):
RHSFull[j, :] = np.kron(TN[j, :], data[:, j])
- G = np.polyfit(self.radiusPade(self.mus), RHSFull,
- deg = self.E, w = self.ws)
+ fullR, G = checkPolyfitRank(self.radiusPade(self.mus),
+ RHSFull, self.E, w = self.ws,
+ rcond = self.interpRcond)
+ if not fullR:
+ Enew = G
+ Nnew = min(self.N, Enew)
+ Mnew = min(self.M, Enew)
+ if Nnew == self.N:
+ strN = ""
+ else:
+ strN = "N from {} to {} and ".format(self.N, Nnew)
+ if Mnew == self.M:
+ strM = ""
+ else:
+ strM = "M from {} to {} and ".format(self.M, Mnew)
+ warn(("Polyfit is poorly conditioned.\nReducing {}{}E "
+ "from {} to {}.").format(strN, strM,
+ self.E, Enew))
+ newParameters = {"N" : Nnew, "M" : Mnew, "E" : Enew}
+ self.approxParameters = newParameters
+ continue
G = G[0, :].reshape((self.N + 1, data.shape[0])).T
if self.POD:
if self.verbosity >= 7:
verbosityDepth("INIT", ("Solving svd for square "
"root of gramian matrix."),
end = "")
_, ev, eV = np.linalg.svd(G, full_matrices = False)
ev = ev[::-1]
eV = eV[::-1, :].conj().T
self._funcOptimum = ev[0]
else:
if self.verbosity >= 10:
verbosityDepth("INIT", "Building gramian matrix.",
end = "")
G2 = self.HFEngine.innerProduct(G, G)
if self.verbosity >= 10:
verbosityDepth("DEL", "Done building gramian.",
inline = True)
if self.verbosity >= 7:
verbosityDepth("INIT", ("Solving eigenvalue "
"problem for gramian "
"matrix."), end = "")
ev, eV = np.linalg.eigh(G2)
self._funcOptimum = ev[0] ** .5
if ev[0] <= 0 or np.isclose(np.abs(ev[0] / ev[-1]), 0.):
- self._funcOptimumWarn = True
+ self._funcOptimumIllCond = True
else:
- self._funcOptimumWarn = False
+ self._funcOptimumIllCond = False
if self.verbosity >= 7:
verbosityDepth("DEL", " Done.", inline = True)
- ts = self.robustTol * np.linalg.norm(ev)
- Nnew = np.sum(np.abs(ev) >= ts)
- diff = self.N - Nnew
- if diff <= 0: break
- Enew = self.E - diff
- Mnew = min(self.M, Enew - 1)
- if Mnew == self.M:
- strM = ""
- else:
- strM = ", M from {0} to {1},".format(self.M, Mnew)
- print(("Smallest {0} eigenvalues below tolerance.\n"
- "Reducing N from {1} to {2}{5} and E from {3} to "
- "{4}.").format(diff + 1, self.N, Nnew,
- self.E, Enew, strM))
- newParameters = {"N" : Nnew, "M" : Mnew, "E" : Enew}
+ newParameters = checkRobustTolerance(ev, self.E,
+ self.robustTol)
+ if not newParameters:
+ break
self.approxParameters = newParameters
- if self.N <= 0:
- eV = np.ones((1, 1))
+ if self.N <= 0:
+ eV = np.ones((1, 1))
self.Q = np.poly1d(eV[:, 0])
if self.verbosity >= 7:
verbosityDepth("DEL", "Done computing denominator.")
else:
- if self.POD:
- data = self.samplingEngine.RPOD
- else:
- data = self.samplingEngine.samples
- g = np.polyfit(self.radiusPade(self.mus), data.T,
- deg = self.E, w = self.ws)[0, :].flatten()
- if self.POD:
- self._funcOptimum = np.linalg.norm(g)
- else:
- self._funcOptimum = self.HFEngine.norm(g)
- self._funcOptimumWarn = False
+ self._funcOptimumIllCond = True
self.Q = np.poly1d([1])
if self.verbosity >= 7:
verbosityDepth("INIT", "Starting computation of numerator.")
Qevaldiag = np.diag(self.Q(self.radiusPade(self.mus)))
- P = np.polyfit(self.radiusPade(self.mus), Qevaldiag, deg = self.M,
- w = self.ws).T
+
+ while self.M >= 0:
+ fullR, P = checkPolyfitRank(self.radiusPade(self.mus),
+ Qevaldiag, self.M, w = self.ws,
+ rcond = self.interpRcond)
+ if fullR:
+ P = P.T
+ break
+ warn(("Polyfit is poorly conditioned. Reducing M from {} to "
+ "{}. Exact snapshot interpolation not guaranteed.")\
+ .format(self.M, P))
+ self.M = P
self.P = np.atleast_2d(P)
if self.POD:
self.P = self.samplingEngine.RPOD.dot(self.P)
if self.verbosity >= 7:
verbosityDepth("DEL", "Done computing numerator.")
self.lastApproxParameters = copy(self.approxParameters)
if hasattr(self, "lastSolvedApp"): del self.lastSolvedApp
if self.verbosity >= 5:
verbosityDepth("DEL", "Done setting up approximant.\n")
def radiusPade(self, mu:Np1D, mu0 : float = None) -> float:
"""
Compute translated radius to be plugged into Pade' approximant.
Args:
mu: Parameter(s) 1.
mu0: Parameter(s) 2. If None, set to self.mu0.
Returns:
Translated radius to be plugged into Pade' approximant.
"""
if mu0 is None: mu0 = self.mu0
return self.HFEngine.rescaling(mu) - self.HFEngine.rescaling(mu0)
def getPVal(self, mu:List[complex]):
"""
Evaluate Pade' numerator at arbitrary parameter.
Args:
mu: Target parameter.
"""
self.setupApprox()
if self.verbosity >= 10:
verbosityDepth("INIT",
"Evaluating numerator at mu = {}.".format(mu),
end = "")
try:
len(mu)
except:
mu = [mu]
- powerlist = np.vander(self.radiusPade(mu), self.M + 1).T
+ powerlist = np.vander(self.radiusPade(mu), self.P.shape[1]).T
p = self.P.dot(powerlist)
if len(mu) == 1:
p = p.flatten()
if self.verbosity >= 10:
verbosityDepth("DEL", " Done.", inline = True)
return p
def getQVal(self, mu:List[complex]):
"""
Evaluate Pade' denominator at arbitrary parameter.
Args:
mu: Target parameter.
"""
self.setupApprox()
if self.verbosity >= 10:
verbosityDepth("INIT",
"Evaluating denominator at mu = {}.".format(mu),
end = "")
q = self.Q(self.radiusPade(mu))
if self.verbosity >= 10:
verbosityDepth("DEL", " Done.", inline = True)
return q
def evalApproxReduced(self, mu:complex):
"""
Evaluate Pade' approximant at arbitrary parameter.
Args:
mu: Target parameter.
"""
self.setupApprox()
if (not hasattr(self, "lastSolvedApp")
or not np.isclose(self.lastSolvedApp, mu)):
if self.verbosity >= 5:
verbosityDepth("INIT",
"Evaluating approximant at mu = {}.".format(mu))
self.uAppReduced = self.getPVal(mu) / self.getQVal(mu)
self.lastSolvedApp = mu
if self.verbosity >= 5:
verbosityDepth("DEL", "Done evaluating approximant.")
def evalApprox(self, mu:complex):
"""
Evaluate approximant at arbitrary parameter.
Args:
mu: Target parameter.
"""
self.evalApproxReduced(mu)
self.uApp = self.samplingEngine.samples.dot(self.uAppReduced)
def getPoles(self) -> Np1D:
"""
Obtain approximant poles.
Returns:
Numpy complex vector of poles.
"""
self.setupApprox()
return self.HFEngine.rescalingInv(self.Q.r
+ self.HFEngine.rescaling(self.mu0))
diff --git a/rrompy/reduction_methods/lagrange/approximant_lagrange_rb.py b/rrompy/reduction_methods/lagrange/approximant_lagrange_rb.py
index 8d0751f..d9a725b 100644
--- a/rrompy/reduction_methods/lagrange/approximant_lagrange_rb.py
+++ b/rrompy/reduction_methods/lagrange/approximant_lagrange_rb.py
@@ -1,292 +1,291 @@
# Copyright (C) 2018 by the RROMPy authors
#
# This file is part of RROMPy.
#
# RROMPy is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# RROMPy is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with RROMPy. If not, see .
#
from copy import copy
import numpy as np
import scipy as sp
from .generic_approximant_lagrange import GenericApproximantLagrange
from rrompy.utilities.base.types import Np1D, DictAny, HFEng
from rrompy.utilities.base import purgeDict, verbosityDepth
from rrompy.utilities.warning_manager import warn
__all__ = ['ApproximantLagrangeRB']
class ApproximantLagrangeRB(GenericApproximantLagrange):
"""
ROM RB approximant computation for parametric problems.
Args:
HFEngine: HF problem solver.
mu0(optional): Default parameter. Defaults to 0.
approxParameters(optional): Dictionary containing values for main
parameters of approximant. Recognized keys are:
- 'POD': whether to compute POD of snapshots; defaults to True;
- 'S': total number of samples current approximant relies upon;
defaults to 2;
- 'sampler': sample point generator; defaults to uniform sampler on
[0, 1];
- 'R': rank for Galerkin projection; defaults to S.
Defaults to empty dict.
Attributes:
HFEngine: HF problem solver.
mu0: Default parameter.
mus: Array of snapshot parameters.
ws: Array of snapshot weigths (unused).
approxRadius: Dummy radius of approximant (i.e. distance from mu0 to
farthest sample point).
approxParameters: Dictionary containing values for main parameters of
approximant. Recognized keys are in parameterList.
parameterList: Recognized keys of approximant parameters:
- 'POD': whether to compute POD of snapshots;
- 'S': total number of samples current approximant relies upon;
- 'sampler': sample point generator;
- 'R': rank for Galerkin projection.
extraApproxParameters: List of approxParameters keys in addition to
mother class's.
S: Number of solution snapshots over which current approximant is
based upon.
sampler: Sample point generator.
R: Rank for Galerkin projection.
POD: Whether to compute POD of snapshots.
samplingEngine: Sampling engine.
projMat: Projection matrix for RB system assembly.
uHF: High fidelity solution with wavenumber lastSolvedHF as numpy
complex vector.
lastSolvedHF: Wavenumber corresponding to last computed high fidelity
solution.
uApp: Last evaluated approximant as numpy complex vector.
lastApproxParameters: List of parameters corresponding to last
computed approximant.
As: List of sparse matrices (in CSC format) representing coefficients
of linear system matrix wrt theta(mu).
bs: List of numpy vectors representing coefficients of linear system
RHS wrt theta(mu).
thetaAs: List of callables representing coefficients of linear system
matrix wrt mu.
thetabs: List of callables representing coefficients of linear system
RHS wrt mu.
ARBs: List of sparse matrices (in CSC format) representing coefficients
of compressed linear system matrix wrt theta(mu).
bRBs: List of numpy vectors representing coefficients of compressed
linear system RHS wrt theta(mu).
"""
def __init__(self, HFEngine:HFEng, mu0 : complex = 0.,
- approxParameters : DictAny = {}, homogeneize : bool = False,
+ approxParameters : DictAny = {}, homogeneized : bool = False,
verbosity : int = 10):
self._preInit()
if not hasattr(self, "parameterList"):
self.parameterList = []
self.parameterList += ["R"]
super().__init__(HFEngine = HFEngine, mu0 = mu0,
approxParameters = approxParameters,
- homogeneize = homogeneize,
+ homogeneized = homogeneized,
verbosity = verbosity)
if self.verbosity >= 10:
verbosityDepth("INIT", "Computing affine blocks of system.")
self.As, self.thetaAs = self.HFEngine.affineBlocksA(self.mu0)
self.bs, self.thetabs = self.HFEngine.affineBlocksb(self.mu0,
- self.homogeneize)
+ self.homogeneized)
if self.verbosity >= 10:
verbosityDepth("DEL", "Done computing affine blocks.")
self._postInit()
def resetSamples(self):
"""Reset samples."""
super().resetSamples()
self.projMat = None
@property
def approxParameters(self):
"""
Value of approximant parameters. Its assignment may change M, N and S.
"""
return self._approxParameters
@approxParameters.setter
def approxParameters(self, approxParams):
approxParameters = purgeDict(approxParams, self.parameterList,
dictname = self.name() + ".approxParameters",
baselevel = 1)
approxParametersCopy = purgeDict(approxParameters, ["R"], True, True,
baselevel = 1)
GenericApproximantLagrange.approxParameters.fset(self,
approxParametersCopy)
keyList = list(approxParameters.keys())
if "R" in keyList:
self.R = approxParameters["R"]
elif hasattr(self, "R"):
self.R = self.R
else:
self.R = self.S
@property
def R(self):
"""Value of R. Its assignment may change S."""
return self._R
@R.setter
def R(self, R):
if R < 0: raise ArithmeticError("R must be non-negative.")
self._R = R
self._approxParameters["R"] = self.R
if hasattr(self, "S") and self.S < self.R:
warn("Prescribed S is too small. Updating S to R.")
self.S = self.R
def checkComputedApprox(self) -> bool:
"""
Check if setup of new approximant is not needed.
Returns:
True if new setup is not needed. False otherwise.
"""
return (self.projMat is not None and super().checkComputedApprox())
def setupApprox(self):
"""Compute RB projection matrix."""
if not self.checkComputedApprox():
if self.verbosity >= 5:
verbosityDepth("INIT", "Setting up {}.". format(self.name()))
self.computeSnapshots()
if self.verbosity >= 7:
verbosityDepth("INIT", "Computing projection matrix.",
end = "")
if self.POD:
U, _, _ = np.linalg.svd(self.samplingEngine.RPOD,
full_matrices = False)
self.projMat = self.samplingEngine.samples.dot(U[:, : self.R])
else:
self.projMat = self.samplingEngine.samples[:, : self.R]
if self.verbosity >= 7:
verbosityDepth("DEL", " Done.", inline = True)
self.lastApproxParameters = copy(self.approxParameters)
if hasattr(self, "lastSolvedApp"): del self.lastSolvedApp
- self.assembleReducedSystem()
if self.verbosity >= 5:
verbosityDepth("DEL", "Done setting up approximant.\n")
def assembleReducedSystem(self):
"""Build affine blocks of RB linear system through projections."""
self.setupApprox()
if self.verbosity >= 10:
verbosityDepth("INIT", "Projecting affine terms of HF model.",
end = "")
- projMatH = self.projMat.T.conjugate()
+ projMatH = self.projMat.T.conj()
self.ARBs = [None] * len(self.As)
self.bRBs = [None] * len(self.bs)
if self.verbosity >= 10:
verbosityDepth("MAIN", ".", end = "", inline = True)
for j in range(len(self.As)):
self.ARBs[j] = projMatH.dot(self.As[j].dot(self.projMat))
if self.verbosity >= 10:
verbosityDepth("MAIN", ".", end = "", inline = True)
for j in range(len(self.bs)):
self.bRBs[j] = projMatH.dot(self.bs[j])
if self.verbosity >= 10:
verbosityDepth("DEL", "Done.", inline = True)
def solveReducedSystem(self, mu:complex) -> Np1D:
"""
Solve RB linear system.
Args:
mu: Target parameter.
Returns:
Solution of RB linear system.
"""
- self.setupApprox()
+ self.assembleReducedSystem()
if self.verbosity >= 10:
verbosityDepth("INIT",
"Assembling reduced model for mu = {}.".format(mu),
end = "")
ARBmu = self.thetaAs(mu, 0) * self.ARBs[0][:self.R,:self.R]
bRBmu = self.thetabs(mu, 0) * self.bRBs[0][:self.R]
if self.verbosity >= 10:
verbosityDepth("MAIN", ".", end = "", inline = True)
for j in range(1, len(self.ARBs)):
ARBmu += self.thetaAs(mu, j) * self.ARBs[j][:self.R, :self.R]
if self.verbosity >= 10:
verbosityDepth("MAIN", ".", end = "", inline = True)
for j in range(1, len(self.bRBs)):
bRBmu += self.thetabs(mu, j) * self.bRBs[j][:self.R]
if self.verbosity >= 10:
verbosityDepth("DEL", "Done.", inline = True)
if self.verbosity >= 5:
verbosityDepth("INIT",
"Solving reduced model for mu = {}.".format(mu),
end = "")
uRB = np.linalg.solve(ARBmu, bRBmu)
if self.verbosity >= 5:
verbosityDepth("DEL", " Done.", inline = True)
return uRB
def evalApproxReduced(self, mu:complex):
"""
Evaluate RB approximant at arbitrary wavenumber.
Args:
mu: Target parameter.
"""
- self.setupApprox()
+ self.assembleReducedSystem()
if (not hasattr(self, "lastSolvedApp")
or not np.isclose(self.lastSolvedApp, mu)):
if self.verbosity >= 5:
verbosityDepth("INIT",
"Computing RB solution at mu = {}.".format(mu))
self.uAppReduced = self.solveReducedSystem(mu)
self.lastSolvedApp = mu
if self.verbosity >= 5:
verbosityDepth("DEL", "Done computing RB solution.")
def evalApprox(self, mu:complex):
"""
Evaluate approximant at arbitrary parameter.
Args:
mu: Target parameter.
"""
self.evalApproxReduced(mu)
self.uApp = self.projMat[:, :self.R].dot(self.uAppReduced)
def getPoles(self) -> Np1D:
"""
Obtain approximant poles.
Returns:
Numpy complex vector of poles.
"""
warn(("Impossible to compute poles in general affine parameter "
"dependence. Results subject to interpretation/rescaling, or "
"possibly completely wrong."))
- self.setupApprox()
+ self.assembleReducedSystem()
if len(self.ARBs) < 2:
return
A = np.eye(self.ARBs[0].shape[0] * (len(self.ARBs) - 1),
dtype = np.complex)
B = np.zeros_like(A)
A[: self.ARBs[0].shape[0], : self.ARBs[0].shape[0]] = - self.ARBs[0]
for j in range(len(self.ARBs) - 1):
Aj = self.ARBs[j + 1]
B[: Aj.shape[0], j * Aj.shape[0] : (j + 1) * Aj.shape[0]] = Aj
II = np.arange(self.ARBs[0].shape[0],
self.ARBs[0].shape[0] * (len(self.ARBs) - 1))
B[II, II - self.ARBs[0].shape[0]] = 1.
return self.HFEngine.rescalingInv(sp.linalg.eigvals(A, B)
+ self.HFEngine.rescaling(self.mu0))
diff --git a/rrompy/reduction_methods/lagrange/generic_approximant_lagrange.py b/rrompy/reduction_methods/lagrange/generic_approximant_lagrange.py
index 8c3c94b..00756f3 100644
--- a/rrompy/reduction_methods/lagrange/generic_approximant_lagrange.py
+++ b/rrompy/reduction_methods/lagrange/generic_approximant_lagrange.py
@@ -1,217 +1,202 @@
# Copyright (C) 2018 by the RROMPy authors
#
# This file is part of RROMPy.
#
# RROMPy is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# RROMPy is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with RROMPy. If not, see .
#
import numpy as np
from rrompy.reduction_methods.base.generic_approximant import (
GenericApproximant)
from rrompy.utilities.base.types import DictAny, HFEng
from rrompy.utilities.base import purgeDict, verbosityDepth
__all__ = ['GenericApproximantLagrange']
class GenericApproximantLagrange(GenericApproximant):
"""
ROM Lagrange interpolant computation for parametric problems (ABSTRACT).
Args:
HFEngine: HF problem solver.
mu0(optional): Default parameter. Defaults to 0.
approxParameters(optional): Dictionary containing values for main
parameters of approximant. Recognized keys are:
- 'POD': whether to compute POD of snapshots; defaults to True;
- 'S': total number of samples current approximant relies upon;
- 'sampler': sample point generator; defaults to uniform sampler on
[0, 1].
Defaults to empty dict.
Attributes:
HFEngine: HF problem solver.
mu0: Default parameter.
mus: Array of snapshot parameters.
ws: Array of snapshot weigths.
approxParameters: Dictionary containing values for main parameters of
approximant. Recognized keys are in parameterList.
parameterList: Recognized keys of approximant parameters:
- 'POD': whether to compute POD of snapshots;
- 'S': total number of snapshots current approximant relies upon;
- 'sampler': sample point generator.
extraApproxParameters: List of approxParameters keys in addition to
mother class's.
S: Number of solution snapshots over which current approximant is
based upon.
sampler: Sample point generator.
POD: Whether to compute POD of snapshots.
samplingEngine: Sampling engine.
uHF: High fidelity solution with wavenumber lastSolvedHF as numpy
complex vector.
lastSolvedHF: Wavenumber corresponding to last computed high fidelity
solution.
uApp: Last evaluated approximant as numpy complex vector.
lastApproxParameters: List of parameters corresponding to last
computed approximant.
"""
def __init__(self, HFEngine:HFEng, mu0 : complex = 0.,
- approxParameters : DictAny = {}, homogeneize : bool = False,
+ approxParameters : DictAny = {}, homogeneized : bool = False,
verbosity : int = 10):
self._preInit()
if not hasattr(self, "parameterList"):
self.parameterList = []
self.parameterList += ["S", "sampler"]
super().__init__(HFEngine = HFEngine, mu0 = mu0,
approxParameters = approxParameters,
- homogeneize = homogeneize,
+ homogeneized = homogeneized,
verbosity = verbosity)
self._postInit()
def setupSampling(self):
"""Setup sampling engine."""
if not hasattr(self, "POD"): return
if self.POD:
from rrompy.sampling.scipy.sampling_engine_lagrange_pod import (
SamplingEngineLagrangePOD)
super().setupSampling(SamplingEngineLagrangePOD)
else:
from rrompy.sampling.scipy.sampling_engine_lagrange import (
SamplingEngineLagrange)
super().setupSampling(SamplingEngineLagrange)
@property
def mus(self):
"""Value of mus. Its assignment may reset snapshots."""
return self._mus
@mus.setter
def mus(self, mus):
musOld = self.mus if hasattr(self, 'mus') else None
self._mus = np.array(mus)
_, musCounts = np.unique(self._mus, return_counts = True)
if len(np.where(musCounts > 1)[0]) > 0:
raise Exception("Repeated sample points not allowed.")
if (musOld is None or len(self.mus) != len(musOld)
or not np.allclose(self.mus, musOld, 1e-14)):
self.resetSamples()
self.autoNode = None
@property
def approxParameters(self):
"""Value of approximant parameters. Its assignment may change S."""
return self._approxParameters
@approxParameters.setter
def approxParameters(self, approxParams):
approxParameters = purgeDict(approxParams, self.parameterList,
dictname = self.name() + ".approxParameters",
baselevel = 1)
approxParametersCopy = purgeDict(approxParameters, ["S", "sampler"],
True, True, baselevel = 1)
GenericApproximant.approxParameters.fset(self, approxParametersCopy)
keyList = list(approxParameters.keys())
if "S" in keyList:
self.S = approxParameters["S"]
elif hasattr(self, "S"):
self.S = self.S
else:
self.S = 2
if "sampler" in keyList:
self.sampler = approxParameters["sampler"]
elif not hasattr(self, "S"):
from rrompy.utilities.parameter_sampling import QuadratureSampler
self.sampler = QuadratureSampler([0., 1.], "UNIFORM")
del QuadratureSampler
-
- @property
- def POD(self):
- """Value of POD."""
- return self._POD
- @POD.setter
- def POD(self, POD):
- if hasattr(self, "POD"): PODold = self.POD
- else: PODold = -1
- self._POD = POD
- self._approxParameters["POD"] = self.POD
- if PODold != self.POD:
- self.setupSampling()
- self.resetSamples()
@property
def S(self):
"""Value of S."""
return self._S
@S.setter
def S(self, S):
if S <= 0: raise ArithmeticError("S must be positive.")
if hasattr(self, "S"): Sold = self.S
else: Sold = -1
self._S = S
self._approxParameters["S"] = self.S
if Sold != self.S:
self.resetSamples()
@property
def sampler(self):
"""Value of sampler."""
return self._sampler
@sampler.setter
def sampler(self, sampler):
if 'generatePoints' not in dir(sampler):
raise Exception("Sampler type not recognized.")
if hasattr(self, '_sampler'):
samplerOld = self.sampler
self._sampler = sampler
self._approxParameters["sampler"] = self.sampler
if not 'samplerOld' in locals() or samplerOld != self.sampler:
self.resetSamples()
def computeSnapshots(self):
- """
- Compute snapshots of solution map.
- """
+ """Compute snapshots of solution map."""
if self.samplingEngine.samples is None:
if self.verbosity >= 5:
verbosityDepth("INIT", "Starting computation of snapshots.")
self.mus, self.ws = self.sampler.generatePoints(self.S)
self.mus = np.array([x[0] for x in self.mus])
self.samplingEngine.iterSample(self.mus,
- homogeneize = self.homogeneize)
+ homogeneized = self.homogeneized)
if self.verbosity >= 5:
verbosityDepth("DEL", "Done computing snapshots.")
def checkComputedApprox(self) -> bool:
"""
Check if setup of new approximant is not needed.
Returns:
True if new setup is not needed. False otherwise.
"""
return (self.samplingEngine.samples is not None
and super().checkComputedApprox())
def normApp(self, mu:complex, homogeneized : bool = False) -> float:
"""
Compute norm of approximant at arbitrary parameter.
Args:
mu: Target parameter.
homogeneized(optional): Whether to remove Dirichlet BC. Defaults to
False.
Returns:
Target norm of approximant.
"""
- if not self.POD or self.homogeneize != homogeneized:
+ if not self.POD or self.homogeneized != homogeneized:
return super().normApp(mu, homogeneized)
return np.linalg.norm(self.getAppReduced(mu))
+
diff --git a/rrompy/reduction_methods/hermite/__init__.py b/rrompy/reduction_methods/lagrange_greedy/__init__.py
similarity index 62%
rename from rrompy/reduction_methods/hermite/__init__.py
rename to rrompy/reduction_methods/lagrange_greedy/__init__.py
index ed60590..f926ff6 100644
--- a/rrompy/reduction_methods/hermite/__init__.py
+++ b/rrompy/reduction_methods/lagrange_greedy/__init__.py
@@ -1,18 +1,30 @@
# Copyright (C) 2018 by the RROMPy authors
#
# This file is part of RROMPy.
#
# RROMPy is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# RROMPy is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with RROMPy. If not, see .
#
+from .generic_approximant_lagrange_greedy import (
+ GenericApproximantLagrangeGreedy)
+from .approximant_lagrange_greedy_pade import ApproximantLagrangePadeGreedy
+from .approximant_lagrange_greedy_rb import ApproximantLagrangeRBGreedy
+
+__all__ = [
+ 'GenericApproximantLagrangeGreedy',
+ 'ApproximantLagrangePadeGreedy',
+ 'ApproximantLagrangeRBGreedy'
+ ]
+
+
diff --git a/rrompy/reduction_methods/lagrange_greedy/approximant_lagrange_greedy_pade.py b/rrompy/reduction_methods/lagrange_greedy/approximant_lagrange_greedy_pade.py
new file mode 100644
index 0000000..7b5cdbd
--- /dev/null
+++ b/rrompy/reduction_methods/lagrange_greedy/approximant_lagrange_greedy_pade.py
@@ -0,0 +1,294 @@
+# Copyright (C) 2018 by the RROMPy authors
+#
+# This file is part of RROMPy.
+#
+# RROMPy is free software: you can redistribute it and/or modify
+# it under the terms of the GNU Lesser General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# RROMPy is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public License
+# along with RROMPy. If not, see .
+#
+
+from copy import copy
+import numpy as np
+from rrompy.reduction_methods.base import (checkRobustTolerance,
+ checkPolyfitRank)
+from .generic_approximant_lagrange_greedy import (
+ GenericApproximantLagrangeGreedy)
+from rrompy.reduction_methods.lagrange import ApproximantLagrangePade
+from rrompy.utilities.base.types import DictAny, List, HFEng
+from rrompy.utilities.base import purgeDict, verbosityDepth
+from rrompy.utilities.warning_manager import warn
+
+__all__ = ['ApproximantLagrangePadeGreedy']
+
+class ApproximantLagrangePadeGreedy(GenericApproximantLagrangeGreedy,
+ ApproximantLagrangePade):
+ """
+ ROM greedy Lagrange Pade' interpolant computation for parametric problems.
+
+ Args:
+ HFEngine: HF problem solver.
+ mu0(optional): Default parameter. Defaults to 0.
+ approxParameters(optional): Dictionary containing values for main
+ parameters of approximant. Recognized keys are:
+ - 'POD': whether to compute POD of snapshots; defaults to True;
+ - 'muBounds': list of bounds for parameter values; defaults to
+ [[0], [1]];
+ - 'greedyTol': uniform error tolerance for greedy algorithm;
+ defaults to 1e-2;
+ - 'maxIter': maximum number of greedy steps; defaults to 1e2;
+ - 'refinementRatio': ratio of training points to be exhausted
+ before training set refinement; defaults to 0.2;
+ - 'nTrainingPoints': number of training points; defaults to
+ maxIter / refinementRatio;
+ - 'nTestPoints': number of starting test points; defaults to 1;
+ - 'trainingSetGenerator': training sample points generator;
+ defaults to uniform sampler within muBounds;
+ - 'interpRcond': tolerance for interpolation via numpy.polyfit;
+ defaults to None;
+ - 'robustTol': tolerance for robust Pade' denominator management;
+ defaults to 0.
+ Defaults to empty dict.
+
+ Attributes:
+ HFEngine: HF problem solver.
+ mu0: Default parameter.
+ mus: Array of snapshot parameters.
+ approxParameters: Dictionary containing values for main parameters of
+ approximant. Recognized keys are in parameterList.
+ parameterList: Recognized keys of approximant parameters:
+ - 'POD': whether to compute POD of snapshots;
+ - 'muBounds': list of bounds for parameter values;
+ - 'greedyTol': uniform error tolerance for greedy algorithm;
+ - 'maxIter': maximum number of greedy steps;
+ - 'refinementRatio': ratio of training points to be exhausted
+ before training set refinement;
+ - 'nTrainingPoints': number of training points;
+ - 'nTestPoints': number of starting test points;
+ - 'trainingSetGenerator': training sample points generator;
+ - 'interpRcond': tolerance for interpolation via numpy.polyfit;
+ - 'robustTol': tolerance for robust Pade' denominator management.
+ extraApproxParameters: List of approxParameters keys in addition to
+ mother class's.
+ POD: whether to compute POD of snapshots.
+ muBounds: list of bounds for parameter values.
+ greedyTol: uniform error tolerance for greedy algorithm.
+ maxIter: maximum number of greedy steps.
+ refinementRatio: ratio of training points to be exhausted before
+ training set refinement.
+ nTrainingPoints: number of training points.
+ nTestPoints: number of starting test points.
+ trainingSetGenerator: training sample points generator.
+ interpRcond: Tolerance for interpolation via numpy.polyfit.
+ robustTol: Tolerance for robust Pade' denominator management.
+ samplingEngine: Sampling engine.
+ uHF: High fidelity solution with wavenumber lastSolvedHF as numpy
+ complex vector.
+ lastSolvedHF: Wavenumber corresponding to last computed high fidelity
+ solution.
+ Q: Numpy 1D vector containing complex coefficients of approximant
+ denominator.
+ P: Numpy 2D vector whose columns are FE dofs of coefficients of
+ approximant numerator.
+ uApp: Last evaluated approximant as numpy complex vector.
+ lastApproxParameters: List of parameters corresponding to last
+ computed approximant.
+ """
+
+ def __init__(self, HFEngine:HFEng, mu0 : complex = 0.,
+ approxParameters : DictAny = {}, homogeneized : bool = False,
+ verbosity : int = 10):
+ self._preInit()
+ if not hasattr(self, "parameterList"):
+ self.parameterList = []
+ self.parameterList += ["interpRcond", "robustTol"]
+ super().__init__(HFEngine = HFEngine, mu0 = mu0,
+ approxParameters = approxParameters,
+ homogeneized = homogeneized,
+ verbosity = verbosity)
+ self._postInit()
+
+ @property
+ def approxParameters(self):
+ """
+ Value of approximant parameters. Its assignment may change robustTol.
+ """
+ return self._approxParameters
+ @approxParameters.setter
+ def approxParameters(self, approxParams):
+ approxParameters = purgeDict(approxParams, self.parameterList,
+ dictname = self.name() + ".approxParameters",
+ baselevel = 1)
+ approxParametersCopy = purgeDict(approxParameters, ["interpRcond",
+ "robustTol"],
+ True, True, baselevel = 1)
+ GenericApproximantLagrangeGreedy.approxParameters.fset(self,
+ approxParametersCopy)
+ keyList = list(approxParameters.keys())
+ if "interpRcond" in keyList:
+ self.interpRcond = approxParameters["interpRcond"]
+ elif hasattr(self, "interpRcond"):
+ self.interpRcond = self.interpRcond
+ else:
+ self.interpRcond = None
+ if "robustTol" in keyList:
+ self.robustTol = approxParameters["robustTol"]
+ elif hasattr(self, "robustTol"):
+ self.robustTol = self.robustTol
+ else:
+ self.robustTol = 0
+
+ def errorEstimator(self, mus:List[np.complex],
+ homogeneized : bool = None) -> List[np.complex]:
+ """
+ Standard residual-based error estimator. Unreliable for unstable
+ problems.
+ """
+ if homogeneized is None:
+ homogeneized = self.homogeneized
+ self.setupApprox()
+ self.initMassMatrixInverse()
+ nmus = len(mus)
+ jOpt = self.getTargetFunctionalOptimum()
+ if np.isnan(jOpt):
+ err = np.empty_like(mus)
+ err[:] = np.inf
+ return err
+ musTile = np.tile(self.HFEngine.rescaling(mus).reshape(-1, 1),
+ [1, len(self.mus)])
+ smusCol = self.HFEngine.rescaling(self.mus).reshape(1, -1)
+ num = np.abs(np.prod(musTile - smusCol, axis = 1))
+ den = np.abs(self.getQVal(mus))
+ if self.HFEngine.nbs == 1:
+ RHS = self.getRHS(mus[0], homogeneized = homogeneized)
+ err = (jOpt * num / den
+ / self.HFEngine.norm(self.massInv.solve(RHS)))
+ else:
+ err = [0.] * nmus
+ for j in range(nmus):
+ RHS = self.getRHS(mus[j], homogeneized = homogeneized)
+ err[j] = (jOpt * num[j] / den[j]
+ / self.HFEngine.norm(self.massInv.solve(RHS)))
+ return err
+
+ def setupApprox(self):
+ """
+ Compute Pade' interpolant.
+ SVD-based robust eigenvalue management.
+ """
+ if not self.checkComputedApprox():
+ if self.verbosity >= 5:
+ verbosityDepth("INIT", "Setting up {}.". format(self.name()))
+ self.S = len(self.mus)
+ M = self.S - 1
+ N = self.S - 1
+ self.greedy()
+ if N > 0:
+ if self.verbosity >= 7:
+ verbosityDepth("INIT", ("Starting computation of "
+ "denominator."))
+ while N > 0:
+ TN = np.vander(self.radiusPade(self.mus), N = N + 1)
+ if self.POD:
+ data = self.samplingEngine.RPOD
+ else:
+ data = self.samplingEngine.samples
+ RHSFull = np.empty((self.S, data.shape[0] * (N + 1)),
+ dtype = np.complex)
+ for j in range(self.S):
+ RHSFull[j, :] = np.kron(data[:, j], TN[j, :])
+ fullR, G = checkPolyfitRank(self.radiusPade(self.mus),
+ RHSFull, M,
+ rcond = self.interpRcond)
+ if not fullR:
+ warn(("Polyfit is poorly conditioned. Starting "
+ "preemptive termination of computation of "
+ "approximant."))
+ self._funcOptimum = np.nan
+ self._funcOptimumIllCond = False
+ self.Q = np.poly1d([1.])
+ self.P = np.diag([np.nan] * len(self.mus))
+ self.lastApproxParameters = copy(self.approxParameters)
+ if hasattr(self, "lastSolvedApp"):
+ del self.lastSolvedApp
+ if self.verbosity >= 7:
+ verbosityDepth("DEL", ("Aborting computation of "
+ "denominator."))
+ if self.verbosity >= 5:
+ verbosityDepth("DEL", ("Aborting computation of "
+ "approximant.\n"))
+ return
+ G = G[0, :].reshape((data.shape[0], N + 1))
+ if self.POD:
+ if self.verbosity >= 7:
+ verbosityDepth("INIT", ("Solving svd for square "
+ "root of gramian matrix."),
+ end = "")
+ _, ev, eV = np.linalg.svd(G, full_matrices = False)
+ ev = ev[::-1]
+ eV = eV[::-1, :].conj().T
+ self._funcOptimum = ev[0]
+ else:
+ if self.verbosity >= 10:
+ verbosityDepth("INIT", "Building gramian matrix.",
+ end = "")
+ G2 = self.HFEngine.innerProduct(G, G)
+ if self.verbosity >= 10:
+ verbosityDepth("DEL", "Done building gramian.",
+ inline = True)
+ if self.verbosity >= 7:
+ verbosityDepth("INIT", ("Solving eigenvalue "
+ "problem for gramian "
+ "matrix."), end = "")
+ ev, eV = np.linalg.eigh(G2)
+ self._funcOptimum = ev[0] ** .5
+ if ev[0] <= 0 or np.isclose(np.abs(ev[0] / ev[-1]), 0.):
+ self._funcOptimumIllCond = True
+ else:
+ self._funcOptimumIllCond = False
+ if self.verbosity >= 7:
+ verbosityDepth("DEL", " Done.", inline = True)
+ newParameters = checkRobustTolerance(ev, M, self.robustTol)
+ if not newParameters:
+ break
+ N = newParameters["N"]
+ M = newParameters["E"]
+ if N <= 0:
+ eV = np.ones((1, 1))
+ self.Q = np.poly1d(eV[:, 0])
+ if self.verbosity >= 7:
+ verbosityDepth("DEL", "Done computing denominator.")
+ else:
+ self._funcOptimumIllCond = True
+ self.Q = np.poly1d([1])
+ if self.verbosity >= 7:
+ verbosityDepth("INIT", "Starting computation of numerator.")
+ Qevaldiag = np.diag(self.Q(self.radiusPade(self.mus)))
+ while M >= 0:
+ fullR, P = checkPolyfitRank(self.radiusPade(self.mus),
+ Qevaldiag, M,
+ rcond = self.interpRcond)
+ if fullR:
+ P = P.T
+ break
+ warn(("Polyfit is poorly conditioned. Reducing M from {} to "
+ "{}. Exact snapshot interpolation not guaranteed.")\
+ .format(M, P))
+ M = P
+ self.P = np.atleast_2d(P)
+ if self.POD:
+ self.P = self.samplingEngine.RPOD.dot(self.P)
+ if self.verbosity >= 7:
+ verbosityDepth("DEL", "Done computing numerator.")
+ self.lastApproxParameters = copy(self.approxParameters)
+ if hasattr(self, "lastSolvedApp"): del self.lastSolvedApp
+ if self.verbosity >= 5:
+ verbosityDepth("DEL", "Done setting up approximant.\n")
diff --git a/rrompy/reduction_methods/lagrange_greedy/approximant_lagrange_greedy_rb.py b/rrompy/reduction_methods/lagrange_greedy/approximant_lagrange_greedy_rb.py
new file mode 100644
index 0000000..7781de6
--- /dev/null
+++ b/rrompy/reduction_methods/lagrange_greedy/approximant_lagrange_greedy_rb.py
@@ -0,0 +1,319 @@
+# Copyright (C) 2018 by the RROMPy authors
+#
+# This file is part of RROMPy.
+#
+# RROMPy is free software: you can redistribute it and/or modify
+# it under the terms of the GNU Lesser General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# RROMPy is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public License
+# along with RROMPy. If not, see .
+#
+
+import numpy as np
+from copy import copy
+from .generic_approximant_lagrange_greedy import (
+ GenericApproximantLagrangeGreedy)
+from rrompy.reduction_methods.lagrange import ApproximantLagrangeRB
+from rrompy.utilities.base.types import DictAny, HFEng, List
+from rrompy.utilities.base import verbosityDepth
+from rrompy.utilities.warning_manager import warn
+
+__all__ = ['ApproximantLagrangeRBGreedy']
+
+class ApproximantLagrangeRBGreedy(GenericApproximantLagrangeGreedy,
+ ApproximantLagrangeRB):
+ """
+ ROM greedy RB approximant computation for parametric problems.
+
+ Args:
+ HFEngine: HF problem solver.
+ mu0(optional): Default parameter. Defaults to 0.
+ approxParameters(optional): Dictionary containing values for main
+ parameters of approximant. Recognized keys are:
+ - 'POD': whether to compute POD of snapshots; defaults to True;
+ - 'muBounds': list of bounds for parameter values; defaults to
+ [[0], [1]];
+ - 'greedyTol': uniform error tolerance for greedy algorithm;
+ defaults to 1e-2;
+ - 'maxIter': maximum number of greedy steps; defaults to 1e2;
+ - 'refinementRatio': ratio of training points to be exhausted
+ before training set refinement; defaults to 0.2;
+ - 'nTrainingPoints': number of training points; defaults to
+ maxIter / refinementRatio;
+ - 'nTestPoints': number of starting test points; defaults to 1;
+ - 'trainingSetGenerator': training sample points generator;
+ defaults to uniform sampler within muBounds.
+ Defaults to empty dict.
+
+ Attributes:
+ HFEngine: HF problem solver.
+ mu0: Default parameter.
+ mus: Array of snapshot parameters.
+ approxParameters: Dictionary containing values for main parameters of
+ approximant. Recognized keys are in parameterList.
+ parameterList: Recognized keys of approximant parameters:
+ - 'POD': whether to compute POD of snapshots;
+ - 'muBounds': list of bounds for parameter values;
+ - 'greedyTol': uniform error tolerance for greedy algorithm;
+ - 'maxIter': maximum number of greedy steps;
+ - 'refinementRatio': ratio of training points to be exhausted
+ before training set refinement;
+ - 'nTrainingPoints': number of training points;
+ - 'nTestPoints': number of starting test points;
+ - 'trainingSetGenerator': training sample points generator;
+ - 'robustTol': tolerance for robust Pade' denominator management.
+ extraApproxParameters: List of approxParameters keys in addition to
+ mother class's.
+ POD: whether to compute POD of snapshots.
+ muBounds: list of bounds for parameter values.
+ greedyTol: uniform error tolerance for greedy algorithm.
+ maxIter: maximum number of greedy steps.
+ refinementRatio: ratio of training points to be exhausted before
+ training set refinement.
+ nTrainingPoints: number of training points.
+ nTestPoints: number of starting test points.
+ trainingSetGenerator: training sample points generator.
+ samplingEngine: Sampling engine.
+ projMat: Projection matrix for RB system assembly.
+ uHF: High fidelity solution with wavenumber lastSolvedHF as numpy
+ complex vector.
+ lastSolvedHF: Wavenumber corresponding to last computed high fidelity
+ solution.
+ uApp: Last evaluated approximant as numpy complex vector.
+ lastApproxParameters: List of parameters corresponding to last
+ computed approximant.
+ As: List of sparse matrices (in CSC format) representing coefficients
+ of linear system matrix wrt theta(mu).
+ bs: List of numpy vectors representing coefficients of linear system
+ RHS wrt theta(mu).
+ thetaAs: List of callables representing coefficients of linear system
+ matrix wrt mu.
+ thetabs: List of callables representing coefficients of linear system
+ RHS wrt mu.
+ ARBs: List of sparse matrices (in CSC format) representing coefficients
+ of compressed linear system matrix wrt theta(mu).
+ bRBs: List of numpy vectors representing coefficients of compressed
+ linear system RHS wrt theta(mu).
+ """
+
+ def __init__(self, HFEngine:HFEng, mu0 : complex = 0.,
+ approxParameters : DictAny = {}, homogeneized : bool = False,
+ verbosity : int = 10):
+ self._preInit()
+ super().__init__(HFEngine = HFEngine, mu0 = mu0,
+ approxParameters = approxParameters,
+ homogeneized = homogeneized,
+ verbosity = verbosity)
+ if self.verbosity >= 10:
+ verbosityDepth("INIT", "Computing affine blocks of system.")
+ self.As, self.thetaAs = self.HFEngine.affineBlocksA(self.mu0)
+ self.bs, self.thetabs = self.HFEngine.affineBlocksb(self.mu0,
+ self.homogeneized)
+ if self.verbosity >= 10:
+ verbosityDepth("DEL", "Done computing affine blocks.")
+ self._postInit()
+
+ @property
+ def S(self):
+ """Value of S."""
+ return self._S
+ @S.setter
+ def S(self, S):
+ self._S = S
+
+ @property
+ def R(self):
+ """Value of R."""
+ return self._S
+ @R.setter
+ def R(self, R):
+ warn(("R is used just to simplify inheritance, and its value cannot "
+ "be changed from that of S."))
+
+ def resetSamples(self):
+ """Reset samples."""
+ super().resetSamples()
+ self.projMat = None
+ self.resbb = None
+ self.resAb = None
+ self.resAA = None
+
+ def checkComputedApprox(self) -> bool:
+ """
+ Check if setup of new approximant is not needed.
+
+ Returns:
+ True if new setup is not needed. False otherwise.
+ """
+ return (self.projMat is not None and super().checkComputedApprox())
+
+ def errorEstimator(self, mus:List[np.complex],
+ homogeneized : bool = None) -> List[np.complex]:
+ """
+ Standard residual-based error estimator. Unreliable for unstable
+ problems (inf-sup constant is missing).
+ """
+ if homogeneized is None:
+ homogeneized = self.homogeneized
+ nmus = len(mus)
+ err = [0.] * nmus
+ self.assembleReducedSystem()
+ self.assembleReducedResidualBlocks(homogeneized)
+ assert homogeneized == self.homogeneized
+ nAs = len(self.As)
+ nbs = len(self.bs)
+ for j in range(nmus):
+ mu = mus[j]
+ uAppReduced = self.getAppReduced(mu)
+ prodbb = 0.
+ prodAb = 0.
+ prodAA = 0.
+ for i1 in range(nbs):
+ rhobi1 = self.thetabs(mu, i1)
+ for i2 in range(nbs):
+ rhobi2 = self.thetabs(mu, i2).conj()
+ prodbb += (rhobi1 * rhobi2
+ * self.resbb[homogeneized][i1][i2])
+ for i1 in range(nAs):
+ rhoAi1 = self.thetaAs(mu, i1)
+ for i2 in range(nbs):
+ rhobi2 = self.thetabs(mu, i2).conj()
+ prodAb += (rhoAi1 * rhobi2
+ * self.resAb[homogeneized][i1][i2])
+ for i1 in range(nAs):
+ rhoAi1 = self.thetaAs(mu, i1)
+ for i2 in range(nAs):
+ rhoAi2 = self.thetaAs(mu, i2).conj()
+ prodAA += (rhoAi1 * rhoAi2
+ * self.resAA[homogeneized][i1][i2])
+ err[j] = np.abs(((uAppReduced.T.conj().dot(prodAA.dot(uAppReduced))
+ - 2. * np.real(prodAb.dot(uAppReduced))) / prodbb
+ + 1.)[0]) ** .5
+ return err
+
+ def setupApprox(self):
+ """Compute RB projection matrix."""
+ if not self.checkComputedApprox():
+ if self.verbosity >= 5:
+ verbosityDepth("INIT", "Setting up {}.". format(self.name()))
+ self.greedy()
+ self.S = len(self.mus)
+ if self.verbosity >= 7:
+ verbosityDepth("INIT", "Computing projection matrix.",
+ end = "")
+ self.projMat = self.samplingEngine.samples
+ if self.verbosity >= 7:
+ verbosityDepth("DEL", " Done.", inline = True)
+ self.lastApproxParameters = copy(self.approxParameters)
+ if hasattr(self, "lastSolvedApp"): del self.lastSolvedApp
+ if self.verbosity >= 5:
+ verbosityDepth("DEL", "Done setting up approximant.\n")
+
+ def assembleReducedResidualBlocks(self, homogeneized:bool):
+ """Build affine blocks of RB linear system through projections."""
+ self.initMassMatrixInverse()
+ self.assembleReducedSystem()
+ computeResbb = (self.resbb is None
+ or homogeneized not in self.resbb.keys())
+ computeResAb = (self.resAb is None
+ or homogeneized not in self.resAb.keys()
+ or self.resAb[homogeneized][0][0].shape[0] != self.S)
+ computeResAA = (self.resAA is None
+ or homogeneized not in self.resAA.keys()
+ or self.resAA[homogeneized][0][0].shape[0] != self.S)
+ if computeResbb or computeResAb or computeResAA:
+ if self.verbosity >= 7:
+ verbosityDepth("INIT", "Projecting affine terms of residual.",
+ end = "")
+ nAs = len(self.As)
+ nbs = len(self.bs)
+ if computeResbb:
+ if self.resbb is None: self.resbb = {}
+ self.resbb[homogeneized] = []
+ for i in range(nbs):
+ resbbi = []
+ Mbi = self.massInv.solve(self.bs[i])
+ for j in range(i):
+ Mbj = self.massInv.solve(self.bs[j])
+ resbbi += [self.HFEngine.innerProduct(Mbi, Mbj)]
+ resbbi += [self.HFEngine.innerProduct(Mbi, Mbi)]
+ self.resbb[homogeneized] += [resbbi]
+ for i in range(nbs):
+ for j in range(i + 1, nbs):
+ self.resbb[homogeneized][i] += [
+ self.resbb[homogeneized][j][i].conj()]
+ if computeResAb:
+ if self.resAb is None: self.resAb = {}
+ if homogeneized not in self.resAb.keys():
+ self.resAb[homogeneized] = []
+ for i in range(nAs):
+ resAbi = []
+ MAi = self.massInv.solve(self.As[i].dot(self.projMat))
+ for j in range(nbs):
+ Mbj = self.massInv.solve(self.bs[j])
+ resAbi += [self.HFEngine.innerProduct(MAi, Mbj)]
+ self.resAb[homogeneized] += [resAbi]
+ else:
+ Sold = self.resAb[homogeneized][0][0].shape[0]
+ for i in range(nAs):
+ MAi = self.massInv.solve(self.As[i].dot(
+ self.projMat[:, Sold :]))
+ for j in range(nbs):
+ Mbj = self.massInv.solve(self.bs[j])
+ self.resAb[homogeneized][i][j] = np.append(
+ self.resAb[homogeneized][i][j],
+ self.HFEngine.innerProduct(MAi, Mbj))
+ if computeResAA:
+ if self.resAA is None: self.resAA = {}
+ if homogeneized not in self.resAA.keys():
+ self.resAA[homogeneized] = []
+ for i in range(nAs):
+ resAAi = []
+ MAi = self.massInv.solve(self.As[i].dot(self.projMat))
+ for j in range(i):
+ MAj = self.massInv.solve(self.As[j].dot(
+ self.projMat))
+ resAAi += [self.HFEngine.innerProduct(MAi, MAj)]
+ resAAi += [self.HFEngine.innerProduct(MAi, MAi)]
+ self.resAA[homogeneized] += [resAAi]
+ for i in range(nAs):
+ for j in range(i + 1, nAs):
+ self.resAA[homogeneized][i] += [
+ self.resAA[homogeneized][j][i].conj()]
+ else:
+ Sold = self.resAA[homogeneized][0][0].shape[0]
+ for i in range(nAs):
+ MAi = self.massInv.solve(self.As[i].dot(self.projMat))
+ for j in range(i):
+ MAj = self.massInv.solve(self.As[j].dot(
+ self.projMat))
+ resAAcross1 = self.HFEngine.innerProduct(
+ MAi[:, Sold :], MAj[:, : Sold])
+ resAAcross2 = self.HFEngine.innerProduct(
+ MAi[:, : Sold], MAj[:, Sold :])
+ resAAdiag = self.HFEngine.innerProduct(
+ MAi[:, Sold :], MAj[:, Sold :])
+ self.resAA[homogeneized][i][j] = np.block(
+ [[self.resAA[homogeneized][i][j], resAAcross1],
+ [ resAAcross2, resAAdiag]])
+ resAAcross = self.HFEngine.innerProduct(MAi[:, Sold :],
+ MAi[:, : Sold])
+ resAAdiag = self.HFEngine.innerProduct(MAi[:, Sold :],
+ MAi[:, Sold :])
+ self.resAA[homogeneized][i][i] = np.block(
+ [[self.resAA[homogeneized][i][i], resAAcross],
+ [ resAAcross.T.conj(), resAAdiag]])
+ for i in range(nAs):
+ for j in range(i + 1, nAs):
+ self.resAA[homogeneized][i][j] = (
+ self.resAA[homogeneized][j][i].conj())
+ if self.verbosity >= 7:
+ verbosityDepth("DEL", ("Done setting up affine "
+ "decomposition of residual."))
diff --git a/rrompy/reduction_methods/lagrange_greedy/generic_approximant_lagrange_greedy.py b/rrompy/reduction_methods/lagrange_greedy/generic_approximant_lagrange_greedy.py
new file mode 100644
index 0000000..7b4d4ae
--- /dev/null
+++ b/rrompy/reduction_methods/lagrange_greedy/generic_approximant_lagrange_greedy.py
@@ -0,0 +1,452 @@
+# Copyright (C) 2018 by the RROMPy authors
+#
+# This file is part of RROMPy.
+#
+# RROMPy is free software: you can redistribute it and/or modify
+# it under the terms of the GNU Lesser General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# RROMPy is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public License
+# along with RROMPy. If not, see .
+#
+
+import numpy as np
+from rrompy.reduction_methods.base.generic_approximant import (
+ GenericApproximant)
+from rrompy.reduction_methods.lagrange.generic_approximant_lagrange import (
+ GenericApproximantLagrange)
+from rrompy.utilities.base.types import DictAny, HFEng, Tuple, List
+from rrompy.utilities.base import purgeDict, verbosityDepth
+from rrompy.utilities.warning_manager import warn
+
+__all__ = ['GenericApproximantLagrangeGreedy']
+
+class GenericApproximantLagrangeGreedy(GenericApproximantLagrange):
+ """
+ ROM greedy Lagrange interpolant computation for parametric problems
+ (ABSTRACT).
+
+ Args:
+ HFEngine: HF problem solver.
+ mu0(optional): Default parameter. Defaults to 0.
+ approxParameters(optional): Dictionary containing values for main
+ parameters of approximant. Recognized keys are:
+ - 'POD': whether to compute POD of snapshots; defaults to True;
+ - 'muBounds': list of bounds for parameter values; defaults to
+ [[0], [1]];
+ - 'greedyTol': uniform error tolerance for greedy algorithm;
+ defaults to 1e-2;
+ - 'maxIter': maximum number of greedy steps; defaults to 1e2;
+ - 'refinementRatio': ratio of training points to be exhausted
+ before training set refinement; defaults to 0.2;
+ - 'nTrainingPoints': number of training points; defaults to
+ maxIter / refinementRatio;
+ - 'nTestPoints': number of starting test points; defaults to 1;
+ - 'trainingSetGenerator': training sample points generator;
+ defaults to uniform sampler within muBounds;
+ - 'robustTol': tolerance for robust Pade' denominator management;
+ defaults to 0.
+ Defaults to empty dict.
+
+ Attributes:
+ HFEngine: HF problem solver.
+ mu0: Default parameter.
+ mus: Array of snapshot parameters.
+ approxParameters: Dictionary containing values for main parameters of
+ approximant. Recognized keys are in parameterList.
+ parameterList: Recognized keys of approximant parameters:
+ - 'POD': whether to compute POD of snapshots;
+ - 'muBounds': list of bounds for parameter values;
+ - 'greedyTol': uniform error tolerance for greedy algorithm;
+ - 'maxIter': maximum number of greedy steps;
+ - 'refinementRatio': ratio of training points to be exhausted
+ before training set refinement;
+ - 'nTrainingPoints': number of training points;
+ - 'nTestPoints': number of starting test points;
+ - 'trainingSetGenerator': training sample points generator.
+ - 'robustTol': tolerance for robust Pade' denominator management.
+ extraApproxParameters: List of approxParameters keys in addition to
+ mother class's.
+ POD: whether to compute POD of snapshots.
+ muBounds: list of bounds for parameter values.
+ greedyTol: uniform error tolerance for greedy algorithm.
+ maxIter: maximum number of greedy steps.
+ refinementRatio: ratio of training points to be exhausted before
+ training set refinement.
+ nTrainingPoints: number of training points.
+ nTestPoints: number of starting test points.
+ trainingSetGenerator: training sample points generator.
+ robustTol: tolerance for robust Pade' denominator management.
+ samplingEngine: Sampling engine.
+ uHF: High fidelity solution with wavenumber lastSolvedHF as numpy
+ complex vector.
+ lastSolvedHF: Wavenumber corresponding to last computed high fidelity
+ solution.
+ uApp: Last evaluated approximant as numpy complex vector.
+ lastApproxParameters: List of parameters corresponding to last
+ computed approximant.
+ """
+
+ TOL_INSTABILITY = 1e-6
+
+ def __init__(self, HFEngine:HFEng, mu0 : complex = 0.,
+ approxParameters : DictAny = {}, homogeneized : bool = False,
+ verbosity : int = 10):
+ self._preInit()
+ if not hasattr(self, "parameterList"):
+ self.parameterList = []
+ self.parameterList += ["muBounds","greedyTol", "maxIter",
+ "refinementRatio", "nTrainingPoints",
+ "nTestPoints", "trainingSetGenerator"]
+ super(GenericApproximantLagrange, self).__init__(
+ HFEngine = HFEngine, mu0 = mu0,
+ approxParameters = approxParameters,
+ homogeneized = homogeneized,
+ verbosity = verbosity)
+ self._postInit()
+
+ @property
+ def approxParameters(self):
+ """Value of approximant parameters. Its assignment may change S."""
+ return self._approxParameters
+ @approxParameters.setter
+ def approxParameters(self, approxParams):
+ approxParameters = purgeDict(approxParams, self.parameterList,
+ dictname = self.name() + ".approxParameters",
+ baselevel = 1)
+ approxParametersCopy = purgeDict(approxParameters,
+ ["muBounds","greedyTol", "maxIter",
+ "refinementRatio", "nTrainingPoints",
+ "nTestPoints",
+ "trainingSetGenerator"],
+ True, True, baselevel = 1)
+ GenericApproximant.approxParameters.fset(self, approxParametersCopy)
+ keyList = list(approxParameters.keys())
+ if "muBounds" in keyList:
+ self.muBounds = approxParameters["muBounds"]
+ elif hasattr(self, "muBounds"):
+ self.muBounds = self.muBounds
+ else:
+ self.muBounds = [[0.], [1.]]
+ if "greedyTol" in keyList:
+ self.greedyTol = approxParameters["greedyTol"]
+ elif hasattr(self, "greedyTol"):
+ self.greedyTol = self.greedyTol
+ else:
+ self.greedyTol = 1e-2
+ if "maxIter" in keyList:
+ self.maxIter = approxParameters["maxIter"]
+ elif hasattr(self, "maxIter"):
+ self.maxIter = self.maxIter
+ else:
+ self.maxIter = 1e2
+ if "refinementRatio" in keyList:
+ self.refinementRatio = approxParameters["refinementRatio"]
+ elif hasattr(self, "refinementRatio"):
+ self.refinementRatio = self.refinementRatio
+ else:
+ self.refinementRatio = 0.2
+ if "nTrainingPoints" in keyList:
+ self.nTrainingPoints = approxParameters["nTrainingPoints"]
+ elif hasattr(self, "nTrainingPoints"):
+ self.nTrainingPoints = self.nTrainingPoints
+ else:
+ self.nTrainingPoints = np.int(np.ceil(self.maxIter
+ / self.refinementRatio))
+ if "nTestPoints" in keyList:
+ self.nTestPoints = approxParameters["nTestPoints"]
+ elif hasattr(self, "nTestPoints"):
+ self.nTestPoints = self.nTestPoints
+ else:
+ self.nTestPoints = 1
+ if "trainingSetGenerator" in keyList:
+ self.trainingSetGenerator = (
+ approxParameters["trainingSetGenerator"])
+ elif hasattr(self, "trainingSetGenerator"):
+ self.trainingSetGenerator = self.trainingSetGenerator
+ else:
+ from rrompy.utilities.parameter_sampling import QuadratureSampler
+ self.trainingSetGenerator = QuadratureSampler(self.muBounds,
+ "UNIFORM")
+ del QuadratureSampler
+
+ @property
+ def S(self):
+ """Value of S."""
+ return self._S
+ @S.setter
+ def S(self, S):
+ self._S = S
+
+ @property
+ def mus(self):
+ """Value of mus."""
+ return self._mus
+ @mus.setter
+ def mus(self, mus):
+ self._mus = np.array(mus, dtype = np.complex)
+
+ @property
+ def muBounds(self):
+ """Value of muBounds."""
+ return self._muBounds
+ @muBounds.setter
+ def muBounds(self, muBounds):
+ if len(muBounds) != 2:
+ raise Exception("2 limits must be specified.")
+ try:
+ muBounds = muBounds.tolist()
+ except:
+ muBounds = list(muBounds)
+ for j in range(2):
+ try:
+ len(muBounds[j])
+ except:
+ muBounds[j] = np.array([muBounds[j]])
+ if len(muBounds[0]) != len(muBounds[1]):
+ raise Exception("The bounds must have the same length.")
+ self._muBounds = muBounds
+
+ @property
+ def greedyTol(self):
+ """Value of greedyTol."""
+ return self._greedyTol
+ @greedyTol.setter
+ def greedyTol(self, greedyTol):
+ if greedyTol < 0:
+ raise ArithmeticError("greedyTol must be non-negative.")
+ if hasattr(self, "greedyTol"): greedyTolold = self.greedyTol
+ else: greedyTolold = -1
+ self._greedyTol = greedyTol
+ self._approxParameters["greedyTol"] = self.greedyTol
+ if greedyTolold != self.greedyTol:
+ self.resetSamples()
+
+ @property
+ def maxIter(self):
+ """Value of maxIter."""
+ return self._maxIter
+ @maxIter.setter
+ def maxIter(self, maxIter):
+ if maxIter <= 0: raise ArithmeticError("maxIter must be positive.")
+ if hasattr(self, "maxIter"): maxIterold = self.maxIter
+ else: maxIterold = -1
+ self._maxIter = maxIter
+ self._approxParameters["maxIter"] = self.maxIter
+ if maxIterold != self.maxIter:
+ self.resetSamples()
+
+ @property
+ def refinementRatio(self):
+ """Value of refinementRatio."""
+ return self._refinementRatio
+ @refinementRatio.setter
+ def refinementRatio(self, refinementRatio):
+ if refinementRatio <= 0. or refinementRatio > 1.:
+ raise ArithmeticError(("refinementRatio must be between 0 "
+ "(excluded) and 1."))
+ if hasattr(self, "refinementRatio"):
+ refinementRatioold = self.refinementRatio
+ else: refinementRatioold = -1
+ self._refinementRatio = refinementRatio
+ self._approxParameters["refinementRatio"] = self.refinementRatio
+ if refinementRatioold != self.refinementRatio:
+ self.resetSamples()
+
+ @property
+ def nTrainingPoints(self):
+ """Value of nTrainingPoints."""
+ return self._nTrainingPoints
+ @nTrainingPoints.setter
+ def nTrainingPoints(self, nTrainingPoints):
+ if nTrainingPoints <= 1:
+ raise ArithmeticError("nTrainingPoints must be greater than 1.")
+ if not np.isclose(nTrainingPoints, np.int(nTrainingPoints)):
+ raise ArithmeticError("nTrainingPoints must be an integer.")
+ nTrainingPoints = np.int(nTrainingPoints)
+ if hasattr(self, "nTrainingPoints"):
+ nTrainingPointsold = self.nTrainingPoints
+ else: nTrainingPointsold = -1
+ self._nTrainingPoints = nTrainingPoints
+ self._approxParameters["nTrainingPoints"] = self.nTrainingPoints
+ if nTrainingPointsold != self.nTrainingPoints:
+ self.resetSamples()
+
+ @property
+ def nTestPoints(self):
+ """Value of nTestPoints."""
+ return self._nTestPoints
+ @nTestPoints.setter
+ def nTestPoints(self, nTestPoints):
+ if nTestPoints <= 0:
+ raise ArithmeticError("nTestPoints must be positive.")
+ if not np.isclose(nTestPoints, np.int(nTestPoints)):
+ raise ArithmeticError("nTestPoints must be an integer.")
+ nTestPoints = np.int(nTestPoints)
+ if hasattr(self, "nTestPoints"):
+ nTestPointsold = self.nTestPoints
+ else: nTestPointsold = -1
+ self._nTestPoints = nTestPoints
+ self._approxParameters["nTestPoints"] = self.nTestPoints
+ if nTestPointsold != self.nTestPoints:
+ self.resetSamples()
+
+ @property
+ def trainingSetGenerator(self):
+ """Value of trainingSetGenerator."""
+ return self._trainingSetGenerator
+ @trainingSetGenerator.setter
+ def trainingSetGenerator(self, trainingSetGenerator):
+ if 'generatePoints' not in dir(trainingSetGenerator):
+ raise Exception("trainingSetGenerator type not recognized.")
+ if hasattr(self, '_trainingSetGenerator'):
+ trainingSetGeneratorOld = self.trainingSetGenerator
+ self._trainingSetGenerator = trainingSetGenerator
+ self._approxParameters["trainingSetGenerator"] = (
+ self.trainingSetGenerator)
+ if (not 'trainingSetGeneratorOld' in locals()
+ or trainingSetGeneratorOld != self.trainingSetGenerator):
+ self.resetSamples()
+
+ def resetSamples(self):
+ """Reset samples."""
+ super().resetSamples()
+ self._mus = []
+
+ def initMassMatrixInverse(self):
+ """Initialize mass matrix for error estimator."""
+ if not hasattr(self, "massInv"):
+ from scipy.sparse import csc_matrix, linalg as sla
+ from rrompy.hfengines.base import ProblemEngineBase as PEB
+ L2normer = PEB()
+ L2normer.V = self.HFEngine.V
+ L2normer.buildEnergyNormForm()
+ Mass = csc_matrix(L2normer.energyNormMatrix, dtype = np.complex)
+ self.massInv = sla.spilu(Mass)
+ del csc_matrix, sla, PEB
+
+ def errorEstimator(self, mus:List[np.complex],
+ homogeneized : bool = None) -> List[np.complex]:
+ """
+ Standard residual-based error estimator. Unreliable for unstable
+ problems (inf-sup constant is missing).
+ """
+ if homogeneized is None:
+ homogeneized = self.homogeneized
+ self.setupApprox()
+ nmus = len(mus)
+ err = [0.] * nmus
+ self.initMassMatrixInverse()
+ if self.HFEngine.nbs == 1:
+ RHS = self.getRHS(mus[0], homogeneized = homogeneized)
+ RHSNorm = self.HFEngine.norm(self.massInv.solve(RHS))
+ for j in range(nmus):
+ res = self.getRes(mus[j], homogeneized = homogeneized)
+ err[j] = (self.HFEngine.norm(self.massInv.solve(res))
+ / RHSNorm)
+ else:
+ for j in range(nmus):
+ res = self.getRes(mus[j], homogeneized = homogeneized)
+ RHS = self.getRHS(mus[j], homogeneized = homogeneized)
+ err[j] = (self.HFEngine.norm(self.massInv.solve(res))
+ / self.HFEngine.norm(self.massInv.solve(RHS)))
+ return err
+
+ def getMaxErrorEstimator(self) -> Tuple[float, int]:
+ """
+ Compute maximum of (and index of maximum of) error estimator over
+ training set.
+ """
+ errorEstTrain = self.errorEstimator(self.muTrain)
+ idxMaxEst = np.argmax(errorEstTrain)
+ maxEst = errorEstTrain[idxMaxEst]
+ return maxEst, idxMaxEst
+
+ def greedyNextSample(self, muidx:int, plotEst : bool = False):
+ """Compute next greedy snapshot of solution map."""
+ mu = self.muTrain[muidx]
+ if self.verbosity >= 2:
+ verbosityDepth("MAIN", ("Adding {}-th sample point at {} to test "
+ "set.").format(
+ self.samplingEngine.nsamples + 1, mu))
+ self.mus = np.append(self.mus, mu)
+ idxs = np.arange(len(self.muTrain))
+ mask = np.ones_like(idxs, dtype = bool)
+ mask[muidx] = False
+ idxs = idxs[mask]
+ self.muTrain = self.muTrain[idxs]
+ self.samplingEngine.nextSample(mu,
+ homogeneized = self.homogeneized)
+ errorEstTrain = self.errorEstimator(self.muTrain)
+ muidx = np.argmax(errorEstTrain)
+ maxErrorEst = errorEstTrain[muidx]
+ mu = self.muTrain[muidx]
+ if plotEst:
+ from matplotlib import pyplot as plt
+ plt.figure()
+ plt.semilogy(self.muTrain, errorEstTrain, 'k')
+ plt.semilogy(self.muTrain,
+ self.greedyTol * np.ones(len(self.muTrain)),
+ 'r--')
+ plt.semilogy(mu, maxErrorEst, 'xr')
+ plt.grid()
+ plt.show()
+ plt.close()
+ return errorEstTrain, muidx, maxErrorEst, mu
+
+ def greedy(self, plotEst : bool = False):
+ """Compute greedy snapshots of solution map."""
+ if self.samplingEngine.samples is None:
+ if self.verbosity >= 2:
+ verbosityDepth("INIT", "Starting computation of snapshots.")
+ self.resetSamples()
+ nTrain = self.nTrainingPoints
+ self.muTrain, _ = self.trainingSetGenerator.generatePoints(nTrain)
+ self.muTrain = np.array([x[0] for x in self.muTrain])
+ maxErrorEst, muidx, mu = np.inf, 0, self.muTrain[0]
+ while (self.samplingEngine.nsamples < self.maxIter
+ and (maxErrorEst > self.greedyTol
+ or self.samplingEngine.nsamples < self.nTestPoints)):
+ if (1. - self.refinementRatio) * nTrain > len(self.muTrain):
+ muTrainExtra, _ = self.trainingSetGenerator.refine(nTrain)
+ self.muTrain = np.sort(np.append(self.muTrain,
+ muTrainExtra))
+ nTrain += len(muTrainExtra)
+ if self.verbosity >= 5:
+ verbosityDepth("MAIN", ("Enriching training set by {} "
+ "elements.").format(
+ len(muTrainExtra)))
+ muTrainOld, maxErrorEstOld = self.muTrain, maxErrorEst
+ errorEstTrain, muidx, maxErrorEst, mu = self.greedyNextSample(
+ muidx, plotEst)
+ if self.verbosity >= 2:
+ verbosityDepth("MAIN", ("Uniform error estimate {:.4e}.")\
+ .format(maxErrorEst))
+ if maxErrorEstOld < maxErrorEst * self.TOL_INSTABILITY:
+ warn(("Instability in a posteriori estimator. Starting "
+ "preemptive greedy loop termination."))
+ self.muTrain = muTrainOld
+ self.mus = self.mus[:-1]
+ self.samplingEngine.popSample()
+ self.setupApprox()
+ maxErrorEst = maxErrorEstOld
+ break
+ if self.verbosity >= 2:
+ verbosityDepth("DEL", "Done computing snapshots.")
+
+ def checkComputedApprox(self) -> bool:
+ """
+ Check if setup of new approximant is not needed.
+
+ Returns:
+ True if new setup is not needed. False otherwise.
+ """
+ return (hasattr(self, "_S") and self.S == len(self.mus)
+ and super().checkComputedApprox())
+
diff --git a/rrompy/reduction_methods/taylor/approximant_taylor_pade.py b/rrompy/reduction_methods/taylor/approximant_taylor_pade.py
index 7a278e7..aad876e 100644
--- a/rrompy/reduction_methods/taylor/approximant_taylor_pade.py
+++ b/rrompy/reduction_methods/taylor/approximant_taylor_pade.py
@@ -1,588 +1,578 @@
# Copyright (C) 2018 by the RROMPy authors
#
# This file is part of RROMPy.
#
# RROMPy is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# RROMPy is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with RROMPy. If not, see .
#
from copy import copy
import numpy as np
+from rrompy.reduction_methods.base import checkRobustTolerance
from .generic_approximant_taylor import GenericApproximantTaylor
from rrompy.sampling.base.pod_engine import PODEngine
from rrompy.utilities.base.types import Np1D, Np2D, List, Tuple, DictAny
from rrompy.utilities.base.types import HFEng
from rrompy.utilities.base import purgeDict, verbosityDepth
from rrompy.utilities.warning_manager import warn
__all__ = ['ApproximantTaylorPade']
class ApproximantTaylorPade(GenericApproximantTaylor):
"""
ROM single-point fast Pade' approximant computation for parametric
problems.
Args:
HFEngine: HF problem solver.
mu0(optional): Default parameter. Defaults to 0.
approxParameters(optional): Dictionary containing values for main
parameters of approximant. Recognized keys are:
- 'POD': whether to compute POD of snapshots; defaults to True;
- 'rho': weight for computation of original Pade' approximant;
defaults to np.inf, i.e. fast approximant;
- 'M': degree of Pade' approximant numerator; defaults to 0;
- 'N': degree of Pade' approximant denominator; defaults to 0;
- 'E': total number of derivatives current approximant relies upon;
defaults to Emax;
- 'Emax': total number of derivatives of solution map to be
computed; defaults to E;
- 'robustTol': tolerance for robust Pade' denominator management;
defaults to 0;
- 'rescaleParam': whether to rescale parameter during denominator
computation; defaults to True;
- 'sampleType': label of sampling type; available values are:
- 'ARNOLDI': orthogonalization of solution derivatives through
Arnoldi algorithm;
- 'KRYLOV': standard computation of solution derivatives.
Defaults to 'KRYLOV'.
Defaults to empty dict.
Attributes:
HFEngine: HF problem solver.
mu0: Default parameter.
approxParameters: Dictionary containing values for main parameters of
approximant. Recognized keys are in parameterList.
parameterList: Recognized keys of approximant parameters:
- 'POD': whether to compute POD of snapshots;
- 'rho': weight for computation of original Pade' approximant;
- 'M': degree of Pade' approximant numerator;
- 'N': degree of Pade' approximant denominator;
- 'E': total number of derivatives current approximant relies upon;
- 'Emax': total number of derivatives of solution map to be
computed;
- 'robustTol': tolerance for robust Pade' denominator management;
- 'rescaleParam': whether to rescale parameter during denominator
computation;
- 'sampleType': label of sampling type.
POD: Whether to compute QR factorization of derivatives.
rho: Weight of approximant.
M: Numerator degree of approximant.
N: Denominator degree of approximant.
E: Number of solution derivatives over which current approximant is
based upon.
Emax: Total number of solution derivatives to be computed.
robustTol: Tolerance for robust Pade' denominator management.
sampleType: Label of sampling type.
initialHFData: HF problem initial data.
uHF: High fidelity solution with wavenumber lastSolvedHF as numpy
complex vector.
lastSolvedHF: Wavenumber corresponding to last computed high fidelity
solution.
G: Square Numpy 2D vector of size (N+1) corresponding to Pade'
denominator matrix (see paper).
Q: Numpy 1D vector containing complex coefficients of approximant
denominator.
P: Numpy 2D vector whose columns are FE dofs of coefficients of
approximant numerator.
uApp: Last evaluated approximant as numpy complex vector.
lastApproxParameters: List of parameters corresponding to last
computed approximant.
"""
def __init__(self, HFEngine:HFEng, mu0 : complex = 0,
- approxParameters : DictAny = {}, homogeneize : bool = False,
+ approxParameters : DictAny = {}, homogeneized : bool = False,
verbosity : int = 10):
self._preInit()
if not hasattr(self, "parameterList"):
self.parameterList = []
self.parameterList += ["M", "N", "robustTol", "rho", "rescaleParam"]
super().__init__(HFEngine = HFEngine, mu0 = mu0,
approxParameters = approxParameters,
- homogeneize = homogeneize,
+ homogeneized = homogeneized,
verbosity = verbosity)
self._postInit()
@property
def approxParameters(self):
"""Value of approximant parameters."""
return self._approxParameters
@approxParameters.setter
def approxParameters(self, approxParams):
approxParameters = purgeDict(approxParams, self.parameterList,
dictname = self.name() + ".approxParameters",
baselevel = 1)
approxParametersCopy = purgeDict(approxParameters,
["M", "N", "robustTol", "rho"],
True, True, baselevel = 1)
keyList = list(approxParameters.keys())
if "rho" in keyList:
self._rho = approxParameters["rho"]
elif hasattr(self, "rho"):
self._rho = self.rho
else:
self._rho = np.inf
GenericApproximantTaylor.approxParameters.fset(self,
approxParametersCopy)
self.rho = self._rho
if "robustTol" in keyList:
self.robustTol = approxParameters["robustTol"]
elif hasattr(self, "robustTol"):
self.robustTol = self.robustTol
else:
self.robustTol = 0
if "rescaleParam" in keyList:
self.rescaleParam = approxParameters["rescaleParam"]
elif hasattr(self, "rescaleParam"):
self.rescaleParam = self.rescaleParam
else:
self.rescaleParam = True
self._ignoreParWarnings = True
if "M" in keyList:
self.M = approxParameters["M"]
elif hasattr(self, "M"):
self.M = self.M
else:
self.M = 0
del self._ignoreParWarnings
if "N" in keyList:
self.N = approxParameters["N"]
elif hasattr(self, "N"):
self.N = self.N
else:
self.N = 0
@property
def rho(self):
"""Value of rho."""
return self._rho
@rho.setter
def rho(self, rho):
self._rho = np.abs(rho)
self._approxParameters["rho"] = self.rho
@property
def M(self):
"""Value of M. Its assignment may change Emax and E."""
return self._M
@M.setter
def M(self, M):
if M < 0: raise ArithmeticError("M must be non-negative.")
self._M = M
self._approxParameters["M"] = self.M
if not hasattr(self, "_ignoreParWarnings"):
self.checkMNEEmax()
@property
def N(self):
"""Value of N. Its assignment may change Emax."""
return self._N
@N.setter
def N(self, N):
if N < 0: raise ArithmeticError("N must be non-negative.")
self._N = N
self._approxParameters["N"] = self.N
self.checkMNEEmax()
def checkMNEEmax(self):
"""Check consistency of M, N, E, and Emax."""
M = self.M if hasattr(self, "_M") else 0
N = self.N if hasattr(self, "_N") else 0
E = self.E if hasattr(self, "_E") else M + N
Emax = self.Emax if hasattr(self, "_Emax") else M + N
if self.rho == np.inf:
if Emax < max(N, M):
warn(("Prescribed Emax is too small. Updating Emax to "
"max(M, N)."))
self.Emax = max(N, M)
if E < max(N, M):
warn("Prescribed E is too small. Updating E to max(M, N).")
self.E = max(N, M)
else:
if Emax < N + M:
warn("Prescribed Emax is too small. Updating Emax to M + N.")
self.Emax = self.N + M
if E < N + M:
warn("Prescribed E is too small. Updating E to M + N.")
self.E = self.N + M
@property
def robustTol(self):
"""Value of tolerance for robust Pade' denominator management."""
return self._robustTol
@robustTol.setter
def robustTol(self, robustTol):
if robustTol < 0.:
warn("Overriding prescribed negative robustness tolerance to 0.")
robustTol = 0.
self._robustTol = robustTol
self._approxParameters["robustTol"] = self.robustTol
@property
def rescaleParam(self):
"""Value of parameter rescaling during denominator computation."""
return self._rescaleParam
@rescaleParam.setter
def rescaleParam(self, rescaleParam):
if not isinstance(rescaleParam, (bool,)):
warn("Prescribed rescaleParam not recognized. Overriding to True.")
rescaleParam = True
self._rescaleParam = rescaleParam
self._approxParameters["rescaleParam"] = self.rescaleParam
@property
def E(self):
"""Value of E. Its assignment may change Emax."""
return self._E
@E.setter
def E(self, E):
if E < 0: raise ArithmeticError("E must be non-negative.")
self._E = E
self.checkMNEEmax()
self._approxParameters["E"] = self.E
if hasattr(self, "Emax") and self.Emax < self.E:
warn("Prescribed Emax is too small. Updating Emax to E.")
self.Emax = self.E
@property
def Emax(self):
"""Value of Emax. Its assignment may reset computed derivatives."""
return self._Emax
@Emax.setter
def Emax(self, Emax):
if Emax < 0: raise ArithmeticError("Emax must be non-negative.")
if hasattr(self, "Emax"): EmaxOld = self.Emax
else: EmaxOld = -1
if hasattr(self, "_N"): N = self.N
else: N = 0
if hasattr(self, "_M"): M = self.M
else: M = 0
if hasattr(self, "_E"): E = self.E
else: E = 0
if self.rho == np.inf:
if max(N, M, E) > Emax:
warn(("Prescribed Emax is too small. Updating Emax to "
"max(N, M, E)."))
Emax = max(N, M, E)
else:
if max(N + M, E) > Emax:
warn(("Prescribed Emax is too small. Updating Emax to "
"max(N + M, E)."))
Emax = max(N + M, E)
self._Emax = Emax
self._approxParameters["Emax"] = Emax
if EmaxOld >= self.Emax and self.samplingEngine.samples is not None:
self.samplingEngine.samples = self.samplingEngine.samples[:,
: self.Emax + 1]
if (self.sampleType == "ARNOLDI"
and self.samplingEngine.HArnoldi is not None):
self.samplingEngine.HArnoldi = self.samplingEngine.HArnoldi[
: self.Emax + 1,
: self.Emax + 1]
self.samplingEngine.RArnoldi = self.samplingEngine.RArnoldi[
: self.Emax + 1,
: self.Emax + 1]
def getTargetFunctionalOptimum(self):
"""Compute optimum of target LS functional."""
self.setupApprox()
- if self._funcOptimumWarn:
- warn("Target functional optimum numerically zero.")
- return np.abs(self._funcOptimum)
+ if self._funcOptimumIllCond:
+ if self.POD:
+ jOpt = np.linalg.norm(self.P[:, 0])
+ else:
+ jOpt = self.HFEngine.norm(self.samplingEngine.samples.dot(
+ self.P[:, 0]))
+ else:
+ jOpt = np.abs(self._funcOptimum)
+ return jOpt
def setupApprox(self):
"""
Compute Pade' approximant. SVD-based robust eigenvalue management.
"""
if not self.checkComputedApprox():
if self.verbosity >= 5:
verbosityDepth("INIT", "Setting up {}.". format(self.name()))
self.computeDerivatives()
if self.N > 0:
if self.verbosity >= 7:
verbosityDepth("INIT", ("Starting computation of "
"denominator."))
while self.N > 0:
if self.POD:
ev, eV = self.findeveVGQR()
self._funcOptimum = ev[0]
else:
ev, eV = self.findeveVGExplicit()
self._funcOptimum = ev[0] ** .5
if ev[0] <= 0 or np.isclose(np.abs(ev[0] / ev[-1]), 0.):
- self._funcOptimumWarn = True
- else:
- self._funcOptimumWarn = False
- ts = self.robustTol * np.linalg.norm(ev)
- Nnew = np.sum(np.abs(ev) >= ts)
- diff = self.N - Nnew
- if diff <= 0: break
- Enew = self.E - diff
- Mnew = min(self.M, Enew)
- if Mnew == self.M:
- strM = ""
+ self._funcOptimumIllCond = True
else:
- strM = ", M from {0} to {1},".format(self.M, Mnew)
- print(("Smallest {0} eigenvalues below tolerance.\n"
- "Reducing N from {1} to {2}{5} and E from {3} to "
- "{4}.").format(diff + 1, self.N, Nnew, self.E,
- Enew, strM))
- newParameters = {"N" : Nnew, "M" : Mnew, "E" : Enew}
+ self._funcOptimumIllCond = False
+ newParameters = checkRobustTolerance(ev, self.E,
+ self.robustTol)
+ if not newParameters:
+ break
self.approxParameters = newParameters
- if self.N <= 0:
- eV = np.ones((1, 1))
+ if self.N <= 0:
+ eV = np.ones((1, 1))
self.Q = np.poly1d(eV[:, 0])
if self.verbosity >= 7:
verbosityDepth("DEL", "Done computing denominator.")
else:
+ self._funcOptimumIllCond = True
self.Q = np.poly1d([1])
- if self.sampleType == "KRYLOV":
- self._funcOptimum = self.HFEngine.norm(
- self.samplingEngine.samples[:, self.E])
- else:
- self._funcOptimum = np.linalg.norm(
- self.samplingEngine.RArnoldi[:, self.E])
- self._funcOptimumWarn = False
if self.verbosity >= 7:
verbosityDepth("INIT", "Starting computation of numerator.")
self.P = np.zeros((self.Emax + 1, self.M + 1), dtype = np.complex)
for i in range(self.Q.order):
rng = np.arange(self.M + 1 - i)
self.P[rng, - 1 - rng - i] = self.Q[i]
if self.sampleType == "ARNOLDI":
self.P = self.samplingEngine.RArnoldi.dot(self.P)
if self.verbosity >= 7:
verbosityDepth("DEL", "Done computing numerator.")
self.lastApproxParameters = copy(self.approxParameters)
if hasattr(self, "lastSolvedApp"): del self.lastSolvedApp
if self.verbosity >= 5:
verbosityDepth("DEL", "Done setting up approximant.\n")
def rescaleParameter(self, R:Np2D, A : Np2D = None,
exponent : float = 1.) -> Np2D:
"""
Prepare parameter rescaling.
Args:
R: Matrix whose columns need rescaling.
A(optional): Matrix whose diagonal defines scaling factor. If None,
previous value of scaleFactor is used. Defaults to None.
exponent(optional): Exponent of scaling factor in matrix diagonal.
Defaults to 1.
Returns:
Rescaled matrix.
"""
if A is not None:
aDiag = np.diag(A)
scaleCoeffs = np.polyfit(np.arange(A.shape[1]),
np.log(aDiag), 1)
self.scaleFactor = np.exp(- scaleCoeffs[0] / exponent)
return np.multiply(R, np.power(self.scaleFactor,np.arange(R.shape[1])))
def buildG(self):
"""Assemble Pade' denominator matrix."""
self.computeDerivatives()
if self.rho == np.inf:
Nmin = self.E - self.N
else:
Nmin = self.M - self.N + 1
if self.sampleType == "KRYLOV":
DerE = self.samplingEngine.samples[:, Nmin : self.E + 1]
G = self.HFEngine.innerProduct(DerE, DerE)
if self.rescaleParam:
DerE = self.rescaleParameter(DerE, G, 2.)
G = self.HFEngine.innerProduct(DerE, DerE)
else:
RArnE = self.samplingEngine.RArnoldi[: self.E + 1,
Nmin : self.E + 1]
if self.rescaleParam:
RArnE = self.rescaleParameter(RArnE, RArnE[Nmin :, :])
G = RArnE.conj().T.dot(RArnE)
if self.rho == np.inf:
self.G = G
else:
Gbig = G
self.G = np.zeros((self.N + 1, self.N + 1), dtype = np.complex)
for k in range(self.E - self.M):
self.G += self.rho ** (2 * k) * Gbig[k : k + self.N + 1,
k : k + self.N + 1]
def findeveVGExplicit(self) -> Tuple[Np1D, Np2D]:
"""
Compute explicitly eigenvalues and eigenvectors of Pade' denominator
matrix.
"""
if self.verbosity >= 10:
verbosityDepth("INIT", "Building gramian matrix.")
self.buildG()
if self.verbosity >= 10:
verbosityDepth("DEL", "Done building gramian.")
if self.verbosity >= 7:
verbosityDepth("INIT",
"Solving eigenvalue problem for gramian matrix.")
ev, eV = np.linalg.eigh(self.G)
if self.rescaleParam:
eV = self.rescaleParameter(eV.T).T
if self.verbosity >= 7:
verbosityDepth("DEL", "Done solving eigenvalue problem.")
return ev, eV
def findeveVGQR(self) -> Tuple[Np1D, Np2D]:
"""
Compute eigenvalues and eigenvectors of Pade' denominator matrix
through SVD of R factor. See ``Householder triangularization of a
quasimatrix'', L.Trefethen, 2008 for QR algorithm.
Returns:
Eigenvalues in ascending order and corresponding eigenvector
matrix.
"""
self.computeDerivatives()
if self.rho == np.inf:
Nmin = self.E - self.N
else:
Nmin = self.M - self.N + 1
if self.sampleType == "KRYLOV":
A = copy(self.samplingEngine.samples[:, Nmin : self.E + 1])
self.PODEngine = PODEngine(self.HFEngine)
if self.verbosity >= 10:
verbosityDepth("INIT", "Orthogonalizing samples.", end = "")
R = self.PODEngine.QRHouseholder(A, only_R = True)
if self.verbosity >= 10:
verbosityDepth("DEL", " Done.", inline = True)
else:
R = self.samplingEngine.RArnoldi[: self.E + 1, Nmin : self.E + 1]
if self.rescaleParam:
R = self.rescaleParameter(R, R[Nmin :, :])
if self.rho == np.inf:
if self.verbosity >= 10:
verbosityDepth("INIT", ("Solving svd for square root of "
"gramian matrix."), end = "")
_, s, V = np.linalg.svd(R, full_matrices = False)
else:
if self.verbosity >= 10:
verbosityDepth("INIT", ("Building matrix stack for square "
"root of gramian."), end = "")
Rtower = np.zeros((R.shape[0] * (self.E - self.M), self.N + 1),
dtype = np.complex)
for k in range(self.E - self.M):
Rtower[k * R.shape[0] : (k + 1) * R.shape[0], :] = (
self.rho ** k
* R[:, self.M - self.N + 1 + k : self.M + 1 + k])
if self.verbosity >= 10:
verbosityDepth("DEL", " Done.", inline = True)
if self.verbosity >= 7:
verbosityDepth("INIT", ("Solving svd for square root of "
"gramian matrix."), end = "")
_, s, V = np.linalg.svd(Rtower, full_matrices = False)
eV = V.conj().T[:, ::-1]
if self.rescaleParam:
eV = self.rescaleParameter(eV.T).T
if self.verbosity >= 7:
verbosityDepth("DEL", " Done.", inline = True)
return s[::-1], eV
def radiusPade(self, mu:Np1D, mu0 : float = None) -> float:
"""
Compute translated radius to be plugged into Pade' approximant.
Args:
mu: Parameter(s) 1.
mu0: Parameter(s) 2. If None, set to self.mu0.
Returns:
Translated radius to be plugged into Pade' approximant.
"""
if mu0 is None: mu0 = self.mu0
return self.HFEngine.rescaling(mu) - self.HFEngine.rescaling(mu0)
def getPVal(self, mu:List[complex]):
"""
Evaluate Pade' numerator at arbitrary parameter.
Args:
mu: Target parameter.
"""
self.setupApprox()
if self.verbosity >= 10:
verbosityDepth("INIT",
"Evaluating numerator at mu = {}.".format(mu),
end = "")
try:
len(mu)
except:
mu = [mu]
powerlist = np.vander(self.radiusPade(mu), self.M + 1).T
p = self.P.dot(powerlist)
if len(mu) == 1:
p = p.flatten()
if self.verbosity >= 10:
verbosityDepth("DEL", " Done.", inline = True)
return p
def getQVal(self, mu:List[complex]):
"""
Evaluate Pade' denominator at arbitrary parameter.
Args:
mu: Target parameter.
"""
self.setupApprox()
if self.verbosity >= 10:
verbosityDepth("INIT",
"Evaluating denominator at mu = {}.".format(mu),
end = "")
q = self.Q(self.radiusPade(mu))
if self.verbosity >= 10:
verbosityDepth("DEL", " Done.", inline = True)
return q
def evalApproxReduced(self, mu:complex):
"""
Evaluate Pade' approximant at arbitrary parameter.
Args:
mu: Target parameter.
"""
self.setupApprox()
if (not hasattr(self, "lastSolvedApp")
or not np.isclose(self.lastSolvedApp, mu)):
if self.verbosity >= 5:
verbosityDepth("INIT",
"Evaluating approximant at mu = {}.".format(mu))
self.uAppReduced = self.getPVal(mu) / self.getQVal(mu)
self.lastSolvedApp = mu
if self.verbosity >= 5:
verbosityDepth("DEL", "Done evaluating approximant.")
def evalApprox(self, mu:complex):
"""
Evaluate approximant at arbitrary parameter.
Args:
mu: Target parameter.
"""
self.evalApproxReduced(mu)
self.uApp = self.samplingEngine.samples.dot(self.uAppReduced)
def getPoles(self) -> Np1D:
"""
Obtain approximant poles.
Returns:
Numpy complex vector of poles.
"""
self.setupApprox()
return self.HFEngine.rescalingInv(self.Q.r
+ self.HFEngine.rescaling(self.mu0))
diff --git a/rrompy/reduction_methods/taylor/approximant_taylor_rb.py b/rrompy/reduction_methods/taylor/approximant_taylor_rb.py
index aa209f4..5dba3d2 100644
--- a/rrompy/reduction_methods/taylor/approximant_taylor_rb.py
+++ b/rrompy/reduction_methods/taylor/approximant_taylor_rb.py
@@ -1,313 +1,313 @@
# Copyright (C) 2018 by the RROMPy authors
#
# This file is part of RROMPy.
#
# RROMPy is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# RROMPy is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with RROMPy. If not, see .
#
from copy import copy
import numpy as np
import scipy as sp
from .generic_approximant_taylor import GenericApproximantTaylor
from rrompy.sampling.base.pod_engine import PODEngine
from rrompy.utilities.base.types import Np1D, DictAny, HFEng
from rrompy.utilities.base import purgeDict, verbosityDepth
from rrompy.utilities.warning_manager import warn
__all__ = ['ApproximantTaylorRB']
class ApproximantTaylorRB(GenericApproximantTaylor):
"""
ROM single-point fast RB approximant computation for parametric problems
with polynomial dependence up to degree 2.
Args:
HFEngine: HF problem solver.
mu0(optional): Default parameter. Defaults to 0.
approxParameters(optional): Dictionary containing values for main
parameters of approximant. Recognized keys are:
- 'POD': whether to compute POD of snapshots; defaults to True;
- 'R': rank for Galerkin projection; defaults to E + 1;
- 'E': total number of derivatives current approximant relies upon;
defaults to Emax;
- 'Emax': total number of derivatives of solution map to be
computed; defaults to E;
- 'sampleType': label of sampling type; available values are:
- 'ARNOLDI': orthogonalization of solution derivatives through
Arnoldi algorithm;
- 'KRYLOV': standard computation of solution derivatives.
Defaults to 'KRYLOV'.
Defaults to empty dict.
Attributes:
HFEngine: HF problem solver.
mu0: Default parameter.
approxParameters: Dictionary containing values for main parameters of
approximant. Recognized keys are in parameterList.
parameterList: Recognized keys of approximant parameters:
- 'POD': whether to compute POD of snapshots;
- 'R': rank for Galerkin projection;
- 'E': total number of derivatives current approximant relies upon;
- 'Emax': total number of derivatives of solution map to be
computed;
- 'sampleType': label of sampling type.
POD: Whether to compute QR factorization of derivatives.
R: Rank for Galerkin projection.
E: Number of solution derivatives over which current approximant is
based upon.
Emax: Total number of solution derivatives to be computed.
sampleType: Label of sampling type, i.e. 'KRYLOV'.
uHF: High fidelity solution with wavenumber lastSolvedHF as numpy
complex vector.
lastSolvedHF: Wavenumber corresponding to last computed high fidelity
solution.
uApp: Last evaluated approximant as numpy complex vector.
lastApproxParameters: List of parameters corresponding to last
computed approximant.
projMat: Numpy matrix representing projection onto RB space.
projMat: Numpy matrix representing projection onto RB space.
As: List of sparse matrices (in CSC format) representing coefficients
of linear system matrix wrt mu.
bs: List of numpy vectors representing coefficients of linear system
RHS wrt mu.
ARBs: List of sparse matrices (in CSC format) representing RB
coefficients of linear system matrix wrt mu.
bRBs: List of numpy vectors representing RB coefficients of linear
system RHS wrt mu.
"""
def __init__(self, HFEngine:HFEng, mu0 : complex = 0,
- approxParameters : DictAny = {}, homogeneize : bool = False,
+ approxParameters : DictAny = {}, homogeneized : bool = False,
verbosity : int = 10):
self._preInit()
if not hasattr(self, "parameterList"):
self.parameterList = []
self.parameterList += ["R"]
super().__init__(HFEngine = HFEngine, mu0 = mu0,
approxParameters = approxParameters,
- homogeneize = homogeneize,
+ homogeneized = homogeneized,
verbosity = verbosity)
if self.verbosity >= 10:
verbosityDepth("INIT", "Computing affine blocks of system.")
self.As, self.thetaAs = self.HFEngine.affineBlocksA(self.mu0)
self.bs, self.thetabs = self.HFEngine.affineBlocksb(self.mu0,
- self.homogeneize)
+ self.homogeneized)
if self.verbosity >= 10:
verbosityDepth("DEL", "Done computing affine blocks.")
self._postInit()
def resetSamples(self):
"""Reset samples."""
super().resetSamples()
self.projMat = None
@property
def approxParameters(self):
"""
Value of approximant parameters. Its assignment may change M, N and S.
"""
return self._approxParameters
@approxParameters.setter
def approxParameters(self, approxParams):
approxParameters = purgeDict(approxParams, self.parameterList,
dictname = self.name() + ".approxParameters",
baselevel = 1)
approxParametersCopy = purgeDict(approxParameters, ["R"],
True, True, baselevel = 1)
GenericApproximantTaylor.approxParameters.fset(self,
approxParametersCopy)
keyList = list(approxParameters.keys())
if "R" in keyList:
self.R = approxParameters["R"]
else:
self.R = self.E + 1
@property
def POD(self):
"""Value of POD."""
return self._POD
@POD.setter
def POD(self, POD):
GenericApproximantTaylor.POD.fset(self, POD)
if (hasattr(self, "sampleType") and self.sampleType == "ARNOLDI"
and not self.POD):
warn(("Arnoldi sampling implicitly forces POD-type derivative "
"management."))
@property
def sampleType(self):
"""Value of sampleType."""
return self._sampleType
@sampleType.setter
def sampleType(self, sampleType):
GenericApproximantTaylor.sampleType.fset(self, sampleType)
if (hasattr(self, "POD") and not self.POD
and self.sampleType == "ARNOLDI"):
warn(("Arnoldi sampling implicitly forces POD-type derivative "
"management."))
@property
def R(self):
"""Value of R. Its assignment may change S."""
return self._R
@R.setter
def R(self, R):
if R < 0: raise ArithmeticError("R must be non-negative.")
self._R = R
self._approxParameters["R"] = self.R
if hasattr(self, "E") and self.E + 1 < self.R:
warn("Prescribed E is too small. Updating E to R - 1.")
self.E = self.R - 1
def setupApprox(self):
"""Setup RB system."""
if not self.checkComputedApprox():
if self.verbosity >= 5:
verbosityDepth("INIT", "Setting up {}.". format(self.name()))
self.computeDerivatives()
if self.verbosity >= 7:
verbosityDepth("INIT", "Computing projection matrix.",
end = "")
if self.POD and not self.sampleType == "ARNOLDI":
self.PODEngine = PODEngine(self.HFEngine)
self.projMatQ, self.projMatR = self.PODEngine.QRHouseholder(
self.samplingEngine.samples)
if self.POD:
if self.sampleType == "ARNOLDI":
self.projMatR = self.samplingEngine.RArnoldi
self.projMatQ = self.samplingEngine.samples
U, _, _ = np.linalg.svd(self.projMatR[: self.E + 1,
: self.E + 1])
self.projMat = self.projMatQ[:, : self.E + 1].dot(U[:,
: self.R])
else:
self.projMat = self.samplingEngine.samples[:, : self.R + 1]
if self.verbosity >= 7:
verbosityDepth("DEL", " Done.", inline = True)
self.lastApproxParameters = copy(self.approxParameters)
if hasattr(self, "lastSolvedApp"): del self.lastSolvedApp
self.assembleReducedSystem()
if self.verbosity >= 5:
verbosityDepth("DEL", "Done setting up approximant.\n")
def assembleReducedSystem(self):
"""Build affine blocks of RB linear system through projections."""
if not self.checkComputedApprox():
self.setupApprox()
if self.verbosity >= 10:
verbosityDepth("INIT", "Projecting affine terms of HF model.",
end = "")
projMatH = self.projMat.T.conj()
self.ARBs = [None] * len(self.As)
self.bRBs = [None] * len(self.bs)
if self.verbosity >= 10:
verbosityDepth("MAIN", ".", end = "", inline = True)
for j in range(len(self.As)):
self.ARBs[j] = projMatH.dot(self.As[j].dot(self.projMat))
if self.verbosity >= 10:
verbosityDepth("MAIN", ".", end = "", inline = True)
for j in range(len(self.bs)):
self.bRBs[j] = projMatH.dot(self.bs[j])
if self.verbosity >= 10:
verbosityDepth("DEL", "Done.", inline = True)
def solveReducedSystem(self, mu:complex) -> Np1D:
"""
Solve RB linear system.
Args:
mu: Target parameter.
Returns:
Solution of RB linear system.
"""
self.setupApprox()
if self.verbosity >= 10:
verbosityDepth("INIT",
"Assembling reduced model for mu = {}.".format(mu),
end = "")
ARBmu = self.thetaAs(mu, 0) * self.ARBs[0][:self.R,:self.R]
bRBmu = self.thetabs(mu, 0) * self.bRBs[0][:self.R]
if self.verbosity >= 10:
verbosityDepth("MAIN", ".", end = "", inline = True)
for j in range(1, len(self.ARBs)):
ARBmu += self.thetaAs(mu, j) * self.ARBs[j][:self.R, :self.R]
if self.verbosity >= 10:
verbosityDepth("MAIN", ".", end = "", inline = True)
for j in range(1, len(self.bRBs)):
bRBmu += self.thetabs(mu, j) * self.bRBs[j][:self.R]
if self.verbosity >= 10:
verbosityDepth("DEL", "Done.", inline = True)
if self.verbosity >= 5:
verbosityDepth("INIT",
"Solving reduced model for mu = {}.".format(mu),
end = "")
uRB = np.linalg.solve(ARBmu, bRBmu)
if self.verbosity >= 5:
verbosityDepth("DEL", " Done.", inline = True)
return uRB
def evalApproxReduced(self, mu:complex):
"""
Evaluate RB approximant at arbitrary wavenumber.
Args:
mu: Target parameter.
"""
self.setupApprox()
if (not hasattr(self, "lastSolvedApp")
or not np.isclose(self.lastSolvedApp, mu)):
if self.verbosity >= 5:
verbosityDepth("INIT",
"Computing RB solution at mu = {}.".format(mu))
self.uAppReduced = self.solveReducedSystem(mu)
self.lastSolvedApp = mu
if self.verbosity >= 5:
verbosityDepth("DEL", "Done computing RB solution.")
def evalApprox(self, mu:complex):
"""
Evaluate approximant at arbitrary parameter.
Args:
mu: Target parameter.
"""
self.evalApproxReduced(mu)
self.uApp = self.projMat[:, : self.R].dot(self.uAppReduced)
def getPoles(self) -> Np1D:
"""
Obtain approximant poles.
Returns:
Numpy complex vector of poles.
"""
warn(("Impossible to compute poles in general affine parameter "
"dependence. Results subject to interpretation/rescaling, or "
"possibly completely wrong."))
self.setupApprox()
if len(self.ARBs) < 2:
return
A = np.eye(self.ARBs[0].shape[0] * (len(self.ARBs) - 1),
dtype = np.complex)
B = np.zeros_like(A)
A[: self.ARBs[0].shape[0], : self.ARBs[0].shape[0]] = - self.ARBs[0]
for j in range(len(self.ARBs) - 1):
Aj = self.ARBs[j + 1]
B[: Aj.shape[0], j * Aj.shape[0] : (j + 1) * Aj.shape[0]] = Aj
II = np.arange(self.ARBs[0].shape[0],
self.ARBs[0].shape[0] * (len(self.ARBs) - 1))
B[II, II - self.ARBs[0].shape[0]] = 1.
return self.HFEngine.rescalingInv(sp.linalg.eigvals(A, B)
+ self.HFEngine.rescaling(self.mu0))
diff --git a/rrompy/reduction_methods/taylor/generic_approximant_taylor.py b/rrompy/reduction_methods/taylor/generic_approximant_taylor.py
index 6745854..25af630 100644
--- a/rrompy/reduction_methods/taylor/generic_approximant_taylor.py
+++ b/rrompy/reduction_methods/taylor/generic_approximant_taylor.py
@@ -1,243 +1,243 @@
# Copyright (C) 2018 by the RROMPy authors
#
# This file is part of RROMPy.
#
# RROMPy is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# RROMPy is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with RROMPy. If not, see .
#
import numpy as np
from rrompy.reduction_methods.base.generic_approximant import (
GenericApproximant)
from rrompy.utilities.base.types import DictAny, HFEng
from rrompy.utilities.base import purgeDict, verbosityDepth
from rrompy.utilities.warning_manager import warn
__all__ = ['GenericApproximantTaylor']
class GenericApproximantTaylor(GenericApproximant):
"""
ROM single-point approximant computation for parametric problems
(ABSTRACT).
Args:
HFEngine: HF problem solver.
mu0(optional): Default parameter. Defaults to 0.
approxParameters(optional): Dictionary containing values for main
parameters of approximant. Recognized keys are:
- 'POD': whether to compute POD of snapshots; defaults to True;
- 'E': total number of derivatives current approximant relies upon;
defaults to Emax;
- 'Emax': total number of derivatives of solution map to be
computed; defaults to E;
- 'sampleType': label of sampling type; available values are:
- 'ARNOLDI': orthogonalization of solution derivatives through
Arnoldi algorithm;
- 'KRYLOV': standard computation of solution derivatives.
Defaults to 'KRYLOV'.
Defaults to empty dict.
Attributes:
HFEngine: HF problem solver.
mu0: Default parameter.
approxParameters: Dictionary containing values for main parameters of
approximant. Recognized keys are in parameterList.
parameterList: Recognized keys of approximant parameters:
- 'POD': whether to compute POD of snapshots;
- 'E': total number of derivatives current approximant relies upon;
- 'Emax': total number of derivatives of solution map to be
computed;
- 'sampleType': label of sampling type.
POD: Whether to compute QR factorization of derivatives.
E: Number of solution derivatives over which current approximant is
based upon.
Emax: Total number of solution derivatives to be computed.
sampleType: Label of sampling type.
initialHFData: HF problem initial data.
samplingEngine: Sampling engine.
uHF: High fidelity solution with wavenumber lastSolvedHF as numpy
complex vector.
lastSolvedHF: Wavenumber corresponding to last computed high fidelity
solution.
uApp: Last evaluated approximant as numpy complex vector.
lastApproxParameters: List of parameters corresponding to last
computed approximant.
"""
def __init__(self, HFEngine:HFEng, mu0 : complex = 0,
- approxParameters : DictAny = {}, homogeneize : bool = False,
+ approxParameters : DictAny = {}, homogeneized : bool = False,
verbosity : int = 10):
self._preInit()
if not hasattr(self, "parameterList"):
self.parameterList = []
self.parameterList += ["E", "Emax", "sampleType"]
super().__init__(HFEngine = HFEngine, mu0 = mu0,
approxParameters = approxParameters,
- homogeneize = homogeneize,
+ homogeneized = homogeneized,
verbosity = verbosity)
self._postInit()
def setupSampling(self):
"""Setup sampling engine."""
if not hasattr(self, "sampleType"): return
if self.sampleType == "ARNOLDI":
from rrompy.sampling.scipy.sampling_engine_arnoldi import (
SamplingEngineArnoldi)
super().setupSampling(SamplingEngineArnoldi)
elif self.sampleType == "KRYLOV":
from rrompy.sampling.scipy.sampling_engine_krylov import (
SamplingEngineKrylov)
super().setupSampling(SamplingEngineKrylov)
else:
raise Exception("Sample type not recognized.")
@property
def approxParameters(self):
"""
Value of approximant parameters. Its assignment may change E and Emax.
"""
return self._approxParameters
@approxParameters.setter
def approxParameters(self, approxParams):
approxParameters = purgeDict(approxParams, self.parameterList,
dictname = self.name() + ".approxParameters",
baselevel = 1)
approxParametersCopy = purgeDict(approxParameters,
["E", "Emax", "sampleType"],
True, True, baselevel = 1)
GenericApproximant.approxParameters.fset(self, approxParametersCopy)
keyList = list(approxParameters.keys())
if "E" in keyList:
self._E = approxParameters["E"]
self._approxParameters["E"] = self.E
if "Emax" in keyList:
self.Emax = approxParameters["Emax"]
else:
if not hasattr(self, "Emax"):
self.Emax = self.E
else:
self.Emax = self.Emax
else:
if "Emax" in keyList:
self._E = approxParameters["Emax"]
self._approxParameters["E"] = self.E
self.Emax = self.E
else:
if not (hasattr(self, "Emax") and hasattr(self, "E")):
raise Exception("At least one of E and Emax must be set.")
if "sampleType" in keyList:
self.sampleType = approxParameters["sampleType"]
elif hasattr(self, "sampleType"):
self.sampleType = self.sampleType
else:
self.sampleType = "KRYLOV"
@property
def E(self):
"""Value of E. Its assignment may change Emax."""
return self._E
@E.setter
def E(self, E):
if E < 0: raise ArithmeticError("E must be non-negative.")
self._E = E
self._approxParameters["E"] = self.E
if hasattr(self, "Emax") and self.Emax < self.E:
warn("Prescribed E is too large. Updating Emax to E.")
self.Emax = self.E
@property
def Emax(self):
"""Value of Emax. Its assignment may reset computed derivatives."""
return self._Emax
@Emax.setter
def Emax(self, Emax):
if Emax < 0: raise ArithmeticError("Emax must be non-negative.")
if hasattr(self, "Emax"): EmaxOld = self.Emax
else: EmaxOld = -1
self._Emax = Emax
if hasattr(self, "E") and self.Emax < self.E:
warn("Prescribed Emax is too small. Updating Emax to E.")
self.Emax = self.E
else:
self._approxParameters["Emax"] = self.Emax
if (EmaxOld >= self.Emax
and self.samplingEngine.samples is not None):
self.samplingEngine.samples = self.samplingEngine.samples[:,
: self.Emax + 1]
if (self.sampleType == "ARNOLDI"
and self.samplingEngine.HArnoldi is not None):
self.samplingEngine.HArnoldi= self.samplingEngine.HArnoldi[
: self.Emax + 1,
: self.Emax + 1]
self.samplingEngine.RArnoldi= self.samplingEngine.RArnoldi[
: self.Emax + 1,
: self.Emax + 1]
else:
self.resetSamples()
@property
def sampleType(self):
"""Value of sampleType."""
return self._sampleType
@sampleType.setter
def sampleType(self, sampleType):
if hasattr(self, "sampleType"): sampleTypeOld = self.sampleType
else: sampleTypeOld = -1
try:
sampleType = sampleType.upper().strip().replace(" ","")
if sampleType not in ["ARNOLDI", "KRYLOV"]:
raise Exception("Sample type not recognized.")
self._sampleType = sampleType
except:
warn(("Prescribed sampleType not recognized. Overriding to "
"'KRYLOV'."))
self._sampleType = "KRYLOV"
self._approxParameters["sampleType"] = self.sampleType
if sampleTypeOld != self.sampleType:
self.resetSamples()
def computeDerivatives(self):
"""Compute derivatives of solution map starting from order 0."""
if self.samplingEngine.samples is None:
if self.verbosity >= 5:
verbosityDepth("INIT", "Starting computation of derivatives.")
self.samplingEngine.iterSample(self.mu0, self.Emax + 1,
- homogeneize = self.homogeneize)
+ homogeneized = self.homogeneized)
if self.verbosity >= 5:
verbosityDepth("DEL", "Done computing derivatives.")
def checkComputedApprox(self) -> bool:
"""
Check if setup of new approximant is not needed.
Returns:
True if new setup is not needed. False otherwise.
"""
return (self.samplingEngine.samples is not None
and super().checkComputedApprox())
def normApp(self, mu:complex, homogeneized : bool = False) -> float:
"""
Compute norm of approximant at arbitrary parameter.
Args:
mu: Target parameter.
homogeneized(optional): Whether to remove Dirichlet BC. Defaults to
False.
Returns:
Target norm of approximant.
"""
- if self.sampleType != "ARNOLDI" or self.homogeneize != homogeneized:
+ if self.sampleType != "ARNOLDI" or self.homogeneized != homogeneized:
return super().normApp(mu, homogeneized)
return np.linalg.norm(self.getAppReduced(mu, homogeneized))
diff --git a/rrompy/sampling/base/sampling_engine_base.py b/rrompy/sampling/base/sampling_engine_base.py
index 9a537f5..d245681 100644
--- a/rrompy/sampling/base/sampling_engine_base.py
+++ b/rrompy/sampling/base/sampling_engine_base.py
@@ -1,105 +1,122 @@
# Copyright (C) 2018 by the RROMPy authors
#
# This file is part of RROMPy.
#
# RROMPy is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# RROMPy is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with RROMPy. If not, see .
#
+import numpy as np
from rrompy.utilities.base.types import Np1D, HFEng, strLst
from rrompy.utilities.base import verbosityDepth
+from rrompy.utilities.warning_manager import warn
__all__ = ['SamplingEngineBase']
class SamplingEngineBase:
"""HERE"""
nameBase = 0
def __init__(self, HFEngine:HFEng, verbosity : int = 10):
self.verbosity = verbosity
if self.verbosity >= 10:
verbosityDepth("INIT",
"Initializing sampling engine of type {}.".format(
self.name()),
end = "")
self.HFEngine = HFEngine
if self.verbosity >= 10:
verbosityDepth("DEL", " Done.", inline = True)
def name(self) -> str:
return self.__class__.__name__
def __str__(self) -> str:
return self.name()
def __repr__(self) -> str:
return self.__str__() + " at " + hex(id(self))
def resetHistory(self):
self.samples = None
self.nsamples = 0
+ def popSample(self):
+ if hasattr(self, "nsamples") and self.nsamples > 1:
+ if self.samples.shape[1] > self.nsamples:
+ warn(("More than 'nsamples' memory allocated for samples. "
+ "Popping empty sample column."))
+ self.nsamples += 1
+ self.samples = self.samples[:, : -1]
+ self.nsamples -= 1
+ else:
+ self.resetHistory()
+
+ def preallocateSamples(self, u:Np1D, n:int):
+ self.samples = np.empty((u.size, n), dtype = u.dtype)
+ self.samples[:, 0] = u
+
@property
def HFEngine(self):
"""Value of HFEngine. Its assignment resets history."""
return self._HFEngine
@HFEngine.setter
def HFEngine(self, HFEngine):
self._HFEngine = HFEngine
self.resetHistory()
def solveLS(self, mu:complex, RHS : Np1D = None,
homogeneized : bool = False) -> Np1D:
"""
Solve linear system.
Args:
mu: Parameter value.
Returns:
Solution of system.
"""
if self.verbosity >= 5:
verbosityDepth("INIT", "Solving HF model for mu = {}.".format(mu))
u = self.HFEngine.solve(mu, RHS, homogeneized)
if self.verbosity >= 5:
verbosityDepth("DEL", "Done solving HF model.")
return u
def plotSamples(self, name : str = "u", save : str = None,
what : strLst = 'all', saveFormat : str = "eps",
saveDPI : int = 100, **figspecs):
"""
Do some nice plots of the samples.
Args:
name(optional): Name to be shown as title of the plots. Defaults to
'u'.
save(optional): Where to save plot(s). Defaults to None, i.e. no
saving.
what(optional): Which plots to do. If list, can contain 'ABS',
'PHASE', 'REAL', 'IMAG'. If str, same plus wildcard 'ALL'.
Defaults to 'ALL'.
saveFormat(optional): Format for saved plot(s). Defaults to "eps".
saveDPI(optional): DPI for saved plot(s). Defaults to 100.
figspecs(optional key args): Optional arguments for matplotlib
figure creation.
"""
for j in range(self.nsamples):
self.HFEngine.plot(self.samples[:, j],
name = "{}_{}".format(name, j + self.nameBase),
save = save, what = what,
saveFormat = saveFormat, saveDPI = saveDPI,
**figspecs)
diff --git a/rrompy/sampling/scipy/sampling_engine_arnoldi.py b/rrompy/sampling/scipy/sampling_engine_arnoldi.py
index d50de53..38f8133 100644
--- a/rrompy/sampling/scipy/sampling_engine_arnoldi.py
+++ b/rrompy/sampling/scipy/sampling_engine_arnoldi.py
@@ -1,129 +1,140 @@
# Copyright (C) 2018 by the RROMPy authors
#
# This file is part of RROMPy.
#
# RROMPy is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# RROMPy is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with RROMPy. If not, see .
#
from copy import copy
import numpy as np
from rrompy.sampling.base.pod_engine import PODEngine
from .sampling_engine_krylov import SamplingEngineKrylov
from rrompy.utilities.base.types import Np1D
from rrompy.utilities.base import verbosityDepth
__all__ = ['SamplingEngineArnoldi']
class SamplingEngineArnoldi(SamplingEngineKrylov):
"""HERE"""
def resetHistory(self):
super().resetHistory()
self.HArnoldi = None
self.RArnoldi = None
self.RHSs = None
self.samplesAug = None
+ def popSample(self):
+ if hasattr(self, "nsamples") and self.nsamples > 1:
+ self.HArnoldi = self.HArnoldi[: -1, : -1]
+ self.RArnoldi = self.RArnoldi[: -1, : -1]
+ if self.nsamples > 2:
+ self.RHSs = self.RHSs[:, : -1]
+ else:
+ self.RHSs = None
+ self.samplesAug = self.RHSs[self.HFEngine.V.dim() :, : -1]
+ super().popSample()
+
@property
def HFEngine(self):
"""Value of HFEngine. Its assignment resets history."""
return self._HFEngine
@HFEngine.setter
def HFEngine(self, HFEngine):
self._HFEngine = HFEngine
self.resetHistory()
self.PODEngine = PODEngine(self._HFEngine)
def preprocesssamples(self):
ns = self.nsamples
if ns <= 0: return
return self.samplesAug[:, ns - 1].reshape((-1,self.HFEngine.V.dim())).T
def preprocessb(self, mu:complex, overwrite : bool = False,
- homogeneize : bool = False):
+ homogeneized : bool = False):
ns = self.nsamples
- r = super().preprocessb(mu, overwrite, homogeneize)
+ r = super().preprocessb(mu, overwrite, homogeneized)
if ns == 0:
return r
elif ns == 1:
r = r / self.RArnoldi[0, 0]
else:
r = ((r - self.RHSs[:, :ns-1].dot(self.RArnoldi[:ns-1, ns-1]))
/ self.RArnoldi[ns-1, ns-1])
if overwrite:
self.RHSs[:, ns - 1] = r
else:
if ns == 1:
self.RHSs = r.reshape((- 1, 1))
else:
self.RHSs = np.hstack((self.RHSs, r[:, None]))
return r
def postprocessu(self, u:Np1D, overwrite : bool = False):
if self.verbosity >= 10:
verbosityDepth("INIT", "Starting orthogonalization.")
ns = self.nsamples
nsAug = (ns + 1) * self.HFEngine.V.dim()
if ns == 0:
u, h, _ = self.PODEngine.GS(u, np.empty((0, 0)))
r = h[0]
uAug = copy(u)
else:
uAug = np.concatenate((self.samplesAug[self.HFEngine.V.dim()
- nsAug :, ns - 1],
u), axis = None)
u, h, uAug = self.PODEngine.GS(u, self.samples[:, : ns], ns, uAug,
self.samplesAug[- nsAug :, : ns])
if overwrite:
self.HArnoldi[: ns + 1, ns] = h
if ns > 0:
r = self.HArnoldi[: ns + 1, 1 : ns + 1].dot(
self.RArnoldi[: ns, ns - 1])
self.RArnoldi[: ns + 1, ns] = r
self.samplesAug[- nsAug :, ns] = uAug
else:
if ns == 0:
self.HArnoldi = h.reshape((1, 1))
self.RArnoldi = r.reshape((1, 1))
self.samplesAug = uAug.reshape((-1, 1))
else:
self.HArnoldi=np.block([[ self.HArnoldi, h[:-1, None]],
[np.zeros((1, ns)), h[-1]]])
if ns > 0:
r = self.HArnoldi[: ns + 1, 1 : ns + 1].dot(
self.RArnoldi[: ns, ns - 1])
self.RArnoldi=np.block([[ self.RArnoldi, r[:-1, None]],
[np.zeros((1, ns)), r[-1]]])
self.samplesAug=np.vstack((np.zeros((self.HFEngine.V.dim(),
ns)),
self.samplesAug))
self.samplesAug = np.hstack((self.samplesAug, uAug[:, None]))
if self.verbosity >= 10:
verbosityDepth("DEL", "Done orthogonalizing.")
return u
def preallocateSamples(self, u:Np1D, n:int):
super().preallocateSamples(u, n)
h = self.HArnoldi
r = self.RArnoldi
saug = self.samplesAug
self.HArnoldi = np.zeros((n, n), dtype = u.dtype)
self.HArnoldi[0, 0] = h[0, 0]
self.RArnoldi = np.zeros((n, n), dtype = u.dtype)
self.RArnoldi[0, 0] = r[0, 0]
self.RHSs = np.empty((u.size, n - 1), dtype = u.dtype)
self.samplesAug = np.zeros((self.HFEngine.V.dim() * (n + 1), n),
dtype = u.dtype)
self.samplesAug[- self.HFEngine.V.dim() :, 0] = saug[:, 0]
diff --git a/rrompy/sampling/scipy/sampling_engine_krylov.py b/rrompy/sampling/scipy/sampling_engine_krylov.py
index 253d543..539eddf 100644
--- a/rrompy/sampling/scipy/sampling_engine_krylov.py
+++ b/rrompy/sampling/scipy/sampling_engine_krylov.py
@@ -1,87 +1,83 @@
# Copyright (C) 2018 by the RROMPy authors
#
# This file is part of RROMPy.
#
# RROMPy is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# RROMPy is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with RROMPy. If not, see .
#
import numpy as np
from rrompy.sampling.base.sampling_engine_base import SamplingEngineBase
from rrompy.utilities.base.types import Np1D, Np2D
from rrompy.utilities.base import verbosityDepth
__all__ = ['SamplingEngineKrylov']
class SamplingEngineKrylov(SamplingEngineBase):
"""HERE"""
def preprocesssamples(self):
if self.samples is None: return
return self.samples[:, : self.nsamples]
def preprocessb(self, mu:complex, overwrite : bool = False,
- homogeneize : bool = False):
- return self.HFEngine.b(mu, self.nsamples, homogeneized = homogeneize)
+ homogeneized : bool = False):
+ return self.HFEngine.b(mu, self.nsamples, homogeneized = homogeneized)
def postprocessu(self, u:Np1D, overwrite : bool = False):
return u
- def preallocateSamples(self, u:Np1D, n:int):
- self.samples = np.empty((u.size, n), dtype = u.dtype)
- self.samples[:, 0] = u
-
def nextSample(self, mu:complex, overwrite : bool = False,
- homogeneize : bool = False) -> Np1D:
+ homogeneized : bool = False) -> Np1D:
ns = self.nsamples
if self.verbosity >= 10:
verbosityDepth("INIT", ("Setting up computation of {}-th Taylor "
"coefficient.").format(ns))
samplesOld = self.preprocesssamples()
RHS = self.preprocessb(mu, overwrite = overwrite,
- homogeneize = homogeneize)
+ homogeneized = homogeneized)
for i in range(1, ns + 1):
RHS -= self.HFEngine.A(mu, i).dot(samplesOld[:, - i])
if self.verbosity >= 10:
verbosityDepth("DEL", "Done setting up for Taylor coefficient.")
u = self.postprocessu(self.solveLS(mu, RHS = RHS,
- homogeneized = homogeneize),
+ homogeneized = homogeneized),
overwrite = overwrite)
if overwrite:
self.samples[:, ns] = u
else:
if ns == 0:
self.samples = u[:, None]
else:
self.samples = np.hstack((self.samples, u[:, None]))
self.nsamples += 1
return u
def iterSample(self, mu:complex, n:int,
- homogeneize : bool = False) -> Np2D:
+ homogeneized : bool = False) -> Np2D:
if self.verbosity >= 5:
verbosityDepth("INIT", "Starting sampling iterations at mu = {}."\
.format(mu))
if n <= 0:
raise Exception(("Number of Krylov iterations must be positive."))
self.resetHistory()
- u = self.nextSample(mu, homogeneize = homogeneize)
+ u = self.nextSample(mu, homogeneized = homogeneized)
if n > 1:
self.preallocateSamples(u, n)
for _ in range(1, n):
self.nextSample(mu, overwrite = True,
- homogeneize = homogeneize)
+ homogeneized = homogeneized)
if self.verbosity >= 5:
verbosityDepth("DEL", "Finished sampling iterations.")
return self.samples
diff --git a/rrompy/sampling/scipy/sampling_engine_lagrange.py b/rrompy/sampling/scipy/sampling_engine_lagrange.py
index 8ff48f2..67d3f51 100644
--- a/rrompy/sampling/scipy/sampling_engine_lagrange.py
+++ b/rrompy/sampling/scipy/sampling_engine_lagrange.py
@@ -1,69 +1,65 @@
# Copyright (C) 2018 by the RROMPy authors
#
# This file is part of RROMPy.
#
# RROMPy is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# RROMPy is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with RROMPy. If not, see .
#
import numpy as np
from rrompy.sampling.base.sampling_engine_base import SamplingEngineBase
from rrompy.utilities.base.types import Np1D, Np2D
from rrompy.utilities.base import verbosityDepth
__all__ = ['SamplingEngineLagrange']
class SamplingEngineLagrange(SamplingEngineBase):
"""HERE"""
nameBase = 1
def postprocessu(self, u:Np1D, overwrite : bool = False):
return u
- def preallocateSamples(self, u:Np1D, n:int):
- self.samples = np.empty((u.size, n), dtype = u.dtype)
- self.samples[:, 0] = u
-
def nextSample(self, mu:complex, overwrite : bool = False,
- homogeneize : bool = False) -> Np1D:
+ homogeneized : bool = False) -> Np1D:
ns = self.nsamples
- u = self.postprocessu(self.solveLS(mu, homogeneized = homogeneize),
+ u = self.postprocessu(self.solveLS(mu, homogeneized = homogeneized),
overwrite = overwrite)
if overwrite:
self.samples[:, ns] = u
else:
if ns == 0:
self.samples = u[:, None]
else:
self.samples = np.hstack((self.samples, u[:, None]))
self.nsamples += 1
return u
- def iterSample(self, mus:Np1D, homogeneize : bool = False) -> Np2D:
+ def iterSample(self, mus:Np1D, homogeneized : bool = False) -> Np2D:
if self.verbosity >= 5:
verbosityDepth("INIT", "Starting sampling iterations.")
n = mus.size
if n <= 0:
raise Exception(("Number of samples must be positive."))
self.resetHistory()
- u = self.nextSample(mus[0], homogeneize = homogeneize)
+ u = self.nextSample(mus[0], homogeneized = homogeneized)
if n > 1:
self.preallocateSamples(u, n)
for j in range(1, n):
self.nextSample(mus[j], overwrite = True,
- homogeneize = homogeneize)
+ homogeneized = homogeneized)
if self.verbosity >= 5:
verbosityDepth("DEL", "Finished sampling iterations.")
return self.samples
diff --git a/rrompy/sampling/scipy/sampling_engine_lagrange_pod.py b/rrompy/sampling/scipy/sampling_engine_lagrange_pod.py
index b239295..273f17b 100644
--- a/rrompy/sampling/scipy/sampling_engine_lagrange_pod.py
+++ b/rrompy/sampling/scipy/sampling_engine_lagrange_pod.py
@@ -1,70 +1,75 @@
# Copyright (C) 2018 by the RROMPy authors
#
# This file is part of RROMPy.
#
# RROMPy is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# RROMPy is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with RROMPy. If not, see .
#
import numpy as np
from rrompy.sampling.base.pod_engine import PODEngine
from .sampling_engine_lagrange import SamplingEngineLagrange
from rrompy.utilities.base.types import Np1D
from rrompy.utilities.base import verbosityDepth
__all__ = ['SamplingEngineLagrangePOD']
class SamplingEngineLagrangePOD(SamplingEngineLagrange):
"""HERE"""
def resetHistory(self):
super().resetHistory()
self.RPOD = None
+ def popSample(self):
+ if hasattr(self, "nsamples") and self.nsamples > 1:
+ self.RPOD = self.RPOD[: -1, : -1]
+ super().popSample()
+
@property
def HFEngine(self):
"""Value of HFEngine. Its assignment resets history."""
return self._HFEngine
@HFEngine.setter
def HFEngine(self, HFEngine):
self._HFEngine = HFEngine
self.resetHistory()
self.PODEngine = PODEngine(self._HFEngine)
def postprocessu(self, u:Np1D, overwrite : bool = False):
if self.verbosity >= 10:
verbosityDepth("INIT", "Starting orthogonalization.")
ns = self.nsamples
if ns == 0:
u, r, _ = self.PODEngine.GS(u, np.empty((0, 0)))
r = r[0]
else:
u, r, _ = self.PODEngine.GS(u, self.samples[:, : ns], ns)
if overwrite:
self.RPOD[: ns + 1, ns] = r
else:
if ns == 0:
self.RPOD = r.reshape((1, 1))
else:
self.RPOD=np.block([[ self.RPOD, r[:-1, None]],
[np.zeros((1, ns)), r[-1]]])
if self.verbosity >= 10:
verbosityDepth("DEL", "Done orthogonalizing.")
return u
def preallocateSamples(self, u:Np1D, n:int):
super().preallocateSamples(u, n)
r = self.RPOD
self.RPOD = np.zeros((n, n), dtype = u.dtype)
self.RPOD[0, 0] = r[0, 0]
diff --git a/rrompy/utilities/parameter_sampling/fft_sampler.py b/rrompy/utilities/parameter_sampling/fft_sampler.py
index dfdfa58..57c6a88 100644
--- a/rrompy/utilities/parameter_sampling/fft_sampler.py
+++ b/rrompy/utilities/parameter_sampling/fft_sampler.py
@@ -1,61 +1,97 @@
# Copyright (C) 2018 by the RROMPy authors
#
# This file is part of RROMPy.
#
# RROMPy is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# RROMPy is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with RROMPy. If not, see .
#
import numpy as np
from rrompy.utilities.parameter_sampling.generic_sampler import GenericSampler
from rrompy.utilities.base.types import Np1D, List, Tuple
__all__ = ['FFTSampler']
class FFTSampler(GenericSampler):
"""Generator of FFT-type sample points on scaled roots of unity."""
def generatePoints(self, n:Np1D) -> Tuple[List[Np1D], Np1D]:
"""Array of sample points and array of weights."""
super().generatePoints(n)
d = len(self.lims[0])
try:
len(n)
except:
n = np.array([n])
if len(n) != d:
raise Exception(("Numbers of points must have same dimension as "
"limits."))
for j in range(d):
a, b = self.lims[0][j], self.lims[1][j]
if self.scaling is not None:
a, b = self.scaling[j](a), self.scaling[j](b)
c, r = (a + b) / 2., np.abs(a - b) / 2.
xj = c + r * np.exp(1.j *
np.linspace(0, 2 * np.pi, n[j] + 1)[:-1, None])
wj = r / n[j] * np.ones(n[j])
if self.scalingInv is not None:
xj = self.scalingInv[j](xj)
if j == 0:
x = xj
w = wj
xsize = n[0]
else:
x = np.concatenate((np.kron(np.ones(n[j])[:, None], x),
np.kron(xj, np.ones(xsize)[:, None])),
axis = 1)
w = np.multiply(np.kron(np.ones(n[j]), w),
np.kron(wj, np.ones(xsize)))
xsize = xsize * n[j]
return [y.flatten() for y in np.split(x, xsize)], w
+ def refine(self, n:int) -> Tuple[List[Np1D], Np1D]:
+ """
+ Apply refinement. If points are not nested, equivalent to
+ generatePoints([2 * x for x in n]).
+ """
+ super().generatePoints(n)
+ d = len(self.lims[0])
+ try:
+ len(n)
+ except:
+ n = np.array([n])
+ if len(n) != d:
+ raise Exception(("Numbers of points must have same dimension as "
+ "limits."))
+ for j in range(d):
+ a, b = self.lims[0][j], self.lims[1][j]
+ if self.scaling is not None:
+ a, b = self.scaling[j](a), self.scaling[j](b)
+ c, r = (a + b) / 2., np.abs(a - b) / 2.
+ xj = c + r * np.exp(1.j * (np.pi / n[j]
+ + np.linspace(0, 2 * np.pi, n[j] + 1)[:-1, None]))
+ wj = r / n[j] * np.ones(n[j])
+ if self.scalingInv is not None:
+ xj = self.scalingInv[j](xj)
+ if j == 0:
+ x = xj
+ w = wj
+ xsize = n[0]
+ else:
+ x = np.concatenate((np.kron(np.ones(n[j])[:, None], x),
+ np.kron(xj, np.ones(xsize)[:, None])),
+ axis = 1)
+ w = np.multiply(np.kron(np.ones(n[j]), w),
+ np.kron(wj, np.ones(xsize)))
+ xsize = xsize * n[j]
+ return [y.flatten() for y in np.split(x, xsize)], w
diff --git a/rrompy/utilities/parameter_sampling/generic_sampler.py b/rrompy/utilities/parameter_sampling/generic_sampler.py
index c3254f6..766d3cb 100644
--- a/rrompy/utilities/parameter_sampling/generic_sampler.py
+++ b/rrompy/utilities/parameter_sampling/generic_sampler.py
@@ -1,97 +1,107 @@
# Copyright (C) 2018 by the RROMPy authors
#
# This file is part of RROMPy.
#
# RROMPy is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# RROMPy is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with RROMPy. If not, see .
#
from abc import abstractmethod
import numpy as np
from rrompy.utilities.base.types import Np1D, Tuple, GenExpr, List
__all__ = ['GenericSampler']
class GenericSampler:
"""ABSTRACT. Generic generator of sample points."""
def __init__(self, lims:Tuple[Np1D, Np1D], scaling : List[callable] = None,
scalingInv : List[callable] = None):
self.lims = lims
self.scaling = scaling
self.scalingInv = scalingInv
def name(self) -> str:
return self.__class__.__name__
def __str__(self) -> str:
return "{}[{}{}]".format(self.name(),
np.array2string(self.lims[0], separator = '_'),
np.array2string(self.lims[1], separator = '_'))
def __repr__(self) -> str:
return self.__str__() + " at " + hex(id(self))
def __eq__(self, other) -> bool:
return self.__dict__ == other.__dict__
@property
def lims(self):
"""Value of lims."""
return self._lims
@lims.setter
def lims(self, lims):
if len(lims) != 2:
raise Exception("2 limits must be specified.")
try:
lims = lims.tolist()
except:
lims = list(lims)
for j in range(2):
try:
len(lims[j])
except:
lims[j] = np.array([lims[j]])
if len(lims[0]) != len(lims[1]):
raise Exception("The limits must have the same length.")
self._lims = lims
@property
def scaling(self):
"""Value of scaling."""
return self._scaling
@scaling.setter
def scaling(self, scaling):
if scaling is not None and not isinstance(scaling, list):
scaling = [scaling]
self._scaling = scaling
@property
def scalingInv(self):
"""Value of scalingInv."""
return self._scalingInv
@scalingInv.setter
def scalingInv(self, scalingInv):
if scalingInv is not None and not isinstance(scalingInv, list):
scalingInv = [scalingInv]
self._scalingInv = scalingInv
@abstractmethod
def generatePoints(self, n:GenExpr) -> Tuple[List[Np1D], Np1D]:
"""Array of points and array of weights."""
assert ((self.scaling is None
or len(self.scaling) == len(self.lims[0]))
and (self.scalingInv is None
or len(self.scalingInv) == len(self.lims[0])))
pass
+ def refine(self, n:int) -> Tuple[List[Np1D], Np1D]:
+ """
+ Apply refinement. If points are not nested, equivalent to
+ generatePoints([2 * x - 1 for x in n]).
+ """
+ try:
+ len(n)
+ except:
+ n = np.array([n])
+ return self.generatePoints([2 * nj - 1 for nj in n])
diff --git a/rrompy/utilities/parameter_sampling/manual_sampler.py b/rrompy/utilities/parameter_sampling/manual_sampler.py
index d06ea34..c069357 100644
--- a/rrompy/utilities/parameter_sampling/manual_sampler.py
+++ b/rrompy/utilities/parameter_sampling/manual_sampler.py
@@ -1,59 +1,81 @@
# Copyright (C) 2018 by the RROMPy authors
#
# This file is part of RROMPy.
#
# RROMPy is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# RROMPy is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with RROMPy. If not, see .
#
import numpy as np
from rrompy.utilities.parameter_sampling.generic_sampler import GenericSampler
from rrompy.utilities.base.types import Np1D, Tuple, List
from rrompy.utilities.warning_manager import warn
__all__ = ['ManualSampler']
class ManualSampler(GenericSampler):
"""Manual generator of sample points."""
def __init__(self, lims:Tuple[Np1D, Np1D], points:Np1D,
scaling : List[callable] = None,
scalingInv : List[callable] = None):
super().__init__(lims = lims, scaling = scaling,
scalingInv = scalingInv)
self.points = points
def __str__(self) -> str:
return "{}[{}]".format(self.name(), "_".join(map(str, self.points)))
def __repr__(self) -> str:
return self.__str__() + " at " + hex(id(self))
def generatePoints(self, n:int) -> Tuple[List[Np1D], Np1D]:
"""Array of quadrature points and array of weights."""
super().generatePoints(None)
size = 1. / n
for j in range(len(self.lims[0])):
a, b = self.lims[0][j], self.lims[1][j]
if self.scaling is not None:
a, b = self.scaling[j](a), self.scaling[j](b)
size *= np.abs(a - b)
if n > len(self.points):
warn(("Requested more points than given. Looping over first "
"points."))
pts = np.tile(self.points,
[np.int(np.ceil(len(self.points) / n))])[: n]
else:
pts = self.points[: n]
return pts, np.ones(n) * size
+ def refine(self, n:int) -> Tuple[List[Np1D], Np1D]:
+ """
+ Apply refinement. If points are not nested, equivalent to
+ generatePoints(2 * n - 1).
+ """
+ super().generatePoints(None)
+ size = 1. / (n - 1)
+ for j in range(len(self.lims[0])):
+ a, b = self.lims[0][j], self.lims[1][j]
+ if self.scaling is not None:
+ a, b = self.scaling[j](a), self.scaling[j](b)
+ size *= np.abs(a - b)
+ if 2 * n - 1 > len(self.points):
+ warn(("Requested more points than given. Looping over first "
+ "points."))
+ pts = np.tile(self.points,
+ [np.int(np.ceil(len(self.points) / (2 * n - 1)))]
+ )[n : 2 * n - 1]
+ else:
+ pts = self.points[n : 2 * n - 1]
+ return pts, np.ones(n - 1) * size
+
diff --git a/rrompy/utilities/parameter_sampling/quadrature_sampler.py b/rrompy/utilities/parameter_sampling/quadrature_sampler.py
index 3d46842..a46c586 100644
--- a/rrompy/utilities/parameter_sampling/quadrature_sampler.py
+++ b/rrompy/utilities/parameter_sampling/quadrature_sampler.py
@@ -1,93 +1,130 @@
# Copyright (C) 2018 by the RROMPy authors
#
# This file is part of RROMPy.
#
# RROMPy is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# RROMPy is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with RROMPy. If not, see .
#
import numpy as np
from rrompy.utilities.parameter_sampling.generic_sampler import GenericSampler
from rrompy.utilities.base.types import Np1D, Tuple, List
__all__ = ['QuadratureSampler']
class QuadratureSampler(GenericSampler):
"""Generator of quadrature sample points."""
allowedKinds = ["UNIFORM", "CHEBYSHEV", "GAUSSLEGENDRE"]
def __init__(self, lims:Tuple[Np1D, Np1D], kind : str = "UNIFORM",
scaling : List[callable] = None,
scalingInv : List[callable] = None):
super().__init__(lims = lims, scaling = scaling,
scalingInv = scalingInv)
self.kind = kind
def __str__(self) -> str:
return "{}_{}".format(super().__str__(), self.kind)
def __repr__(self) -> str:
return self.__str__() + " at " + hex(id(self))
@property
def kind(self):
"""Value of kind."""
return self._kind
@kind.setter
def kind(self, kind):
if kind.upper() not in self.allowedKinds:
raise Exception("Generator kind not recognized.")
self._kind = kind.upper()
def generatePoints(self, n:Np1D) -> Tuple[List[Np1D], Np1D]:
"""Array of quadrature points and array of weights."""
super().generatePoints(n)
d = len(self.lims[0])
try:
len(n)
except:
n = np.array([n])
if len(n) != d:
raise Exception(("Numbers of points must have same dimension as"
"limits."))
for j in range(d):
a, b = self.lims[0][j], self.lims[1][j]
if self.scaling is not None:
a, b = self.scaling[j](a), self.scaling[j](b)
if self.kind == "UNIFORM":
xj = np.linspace(a, b, n[j])[:, None]
- wj = np.abs(a - b) / n[j] * np.ones(n[j])
+ wj = np.abs(a - b) / (n[j] - 1) * np.ones(n[j])
elif self.kind == "CHEBYSHEV":
nodes, weights = np.polynomial.chebyshev.chebgauss(n[j])
xj = (a + b) / 2 + (a - b) / 2 * nodes[:, None]
wj = np.abs(a - b) / np.pi * weights
elif self.kind == "GAUSSLEGENDRE":
nodes, weights = np.polynomial.legendre.leggauss(n[j])
xj = (a + b) / 2 + (a - b) / 2 * nodes[:, None]
wj = np.abs(a - b) * weights
if self.scalingInv is not None:
xj = self.scalingInv[j](xj)
if j == 0:
x = xj
w = wj
xsize = n[0]
else:
x = np.concatenate((np.kron(np.ones(n[j])[:, None], x),
np.kron(xj, np.ones(xsize)[:, None])),
axis = 1)
w = np.multiply(np.kron(np.ones(n[j]), w),
np.kron(wj, np.ones(xsize)))
xsize = xsize * n[j]
return [y.flatten() for y in np.split(x, xsize)], w
+ def refine(self, n:int) -> Tuple[List[Np1D], Np1D]:
+ """
+ Apply refinement. If points are not nested, equivalent to
+ generatePoints([2 * x - 1 for x in n]).
+ """
+ if self.kind != "UNIFORM": return super().refine(n)
+ super().generatePoints(n)
+ d = len(self.lims[0])
+ try:
+ len(n)
+ except:
+ n = np.array([n])
+ if len(n) != d:
+ raise Exception(("Numbers of points must have same dimension as"
+ "limits."))
+ for j in range(d):
+ a, b = self.lims[0][j], self.lims[1][j]
+ if self.scaling is not None:
+ a, b = self.scaling[j](a), self.scaling[j](b)
+ xj = np.linspace(a + (b - a) / 2. / (n[j] - 1),
+ b + (a - b) / 2. / (n[j] - 1), n[j] - 1)[:, None]
+ wj = np.abs(a - b) / (n[j] - 2) * np.ones(n[j] - 1)
+ if self.scalingInv is not None:
+ xj = self.scalingInv[j](xj)
+ if j == 0:
+ x = xj
+ w = wj
+ xsize = n[0] - 1
+ else:
+ x = np.concatenate((np.kron(np.ones(n[j] - 1)[:, None], x),
+ np.kron(xj, np.ones(xsize)[:, None])),
+ axis = 1)
+ w = np.multiply(np.kron(np.ones(n[j] - 1), w),
+ np.kron(wj, np.ones(xsize)))
+ xsize = xsize * (n[j] - 1)
+ return [y.flatten() for y in np.split(x, xsize)], w
+
diff --git a/rrompy/utilities/parameter_sampling/random_sampler.py b/rrompy/utilities/parameter_sampling/random_sampler.py
index 94dcc57..a61763d 100644
--- a/rrompy/utilities/parameter_sampling/random_sampler.py
+++ b/rrompy/utilities/parameter_sampling/random_sampler.py
@@ -1,70 +1,80 @@
# Copyright (C) 2018 by the RROMPy authors
#
# This file is part of RROMPy.
#
# RROMPy is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# RROMPy is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with RROMPy. If not, see .
#
import numpy as np
from rrompy.utilities.base.sobol import sobolGenerate
from rrompy.utilities.parameter_sampling.generic_sampler import GenericSampler
from rrompy.utilities.base.types import Np1D, Tuple, List
__all__ = ['RandomSampler']
class RandomSampler(GenericSampler):
"""Generator of quadrature sample points."""
allowedKinds = ["UNIFORM", "SOBOL"]
def __init__(self, lims:Tuple[Np1D, Np1D], kind : str = "UNIFORM",
scaling : callable = None, scalingInv : callable = None):
super().__init__(lims = lims, scaling = scaling,
scalingInv = scalingInv)
self.kind = kind
def __str__(self) -> str:
return "{}_{}".format(super().__str__(), self.kind)
def __repr__(self) -> str:
return self.__str__() + " at " + hex(id(self))
@property
def kind(self):
"""Value of kind."""
return self._kind
@kind.setter
def kind(self, kind):
if kind.upper() not in self.allowedKinds:
raise Exception("Generator kind not recognized.")
self._kind = kind.upper()
def generatePoints(self, n:int, seed : int = 0) -> Tuple[List[Np1D], Np1D]:
"""Array of quadrature points and array of weights."""
super().generatePoints(n)
d = len(self.lims[0])
if self.kind == "UNIFORM":
np.random.seed(seed)
x = np.random.uniform(size = (n, d))
else:
x = sobolGenerate(d, n, seed)
for j in range(d):
a, b = self.lims[0][j], self.lims[1][j]
if self.scaling is not None:
a, b = self.scaling[j](a), self.scaling[j](b)
x[:, j] = a + (b - a) * x[:, j]
if self.scalingInv is not None:
x[:, j] = self.scalingInv[j](x[:, j])
return [y.flatten() for y in np.split(x, n)], np.ones(n)
+ def refine(self, n:int, seed : int = 420) -> Tuple[List[Np1D], Np1D]:
+ """
+ Apply refinement. If points are not nested, equivalent to
+ generatePoints([2 * x - 1 for x in n]).
+ """
+ try:
+ len(n)
+ except:
+ n = np.array([n])
+ return self.generatePoints([nj - 1 for nj in n], seed)
diff --git a/rrompy/utilities/parameter_sampling/warping_sampler.py b/rrompy/utilities/parameter_sampling/warping_sampler.py
index 84f0a9d..0f14302 100644
--- a/rrompy/utilities/parameter_sampling/warping_sampler.py
+++ b/rrompy/utilities/parameter_sampling/warping_sampler.py
@@ -1,114 +1,153 @@
# Copyright (C) 2018 by the RROMPy authors
#
# This file is part of RROMPy.
#
# RROMPy is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# RROMPy is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with RROMPy. If not, see .
#
import numpy as np
from rrompy.utilities.parameter_sampling.generic_sampler import GenericSampler
from rrompy.utilities.base.types import Np1D, Tuple, List
__all__ = ['WarpingSampler', 'WarpingFunction']
class WarpingFunction:
"""Wrapper for warping function."""
def __init__(self, other : "WarpingFunction" = None, **kwargs):
if other is not None:
self._call_ = other._call_
self._repr_ = other._repr_
else:
self._call_ = kwargs["call"]
self._repr_ = kwargs["repr"]
if not callable(self._call_):
self._call_ = lambda x: self._call_
if callable(self._repr_):
self._repr_ = self._repr_()
def __call__(self, x):
return self._call_(x)
def __str__(self):
return self._repr_
def __repr__(self) -> str:
return self.__str__() + " at " + hex(id(self))
class WarpingSampler(GenericSampler):
"""Generator of sample points from warping of uniform ones."""
def __init__(self, lims:Tuple[Np1D, Np1D], warping : List[callable],
scaling : List[callable] = None,
scalingInv : List[callable] = None):
super().__init__(lims = lims, scaling = scaling,
scalingInv = scalingInv)
self.warping = warping
def __str__(self) -> str:
return "{}_{}".format(super().__str__(), self.warping)
def __repr__(self) -> str:
return self.__str__() + " at " + hex(id(self))
@property
def warping(self):
"""Value of warping."""
return self._warping
@warping.setter
def warping(self, warping):
if not isinstance(warping, (list, tuple,)):
warping = [warping] * len(self.lims[0])
for j in range(len(warping)):
if not isinstance(warping[j], (WarpingFunction,)):
warping[j] = WarpingFunction(call = warping[j].__call__,
repr = warping[j].__repr__)
self._warping = warping
def generatePoints(self, n:Np1D) -> Tuple[List[Np1D], Np1D]:
"""Array of quadrature points and array of weights."""
super().generatePoints(n)
assert len(self.warping) == len(self.lims[0])
d = len(self.lims[0])
try:
len(n)
except:
n = np.array([n])
if len(n) != d:
raise Exception(("Numbers of points must have same dimension as"
"limits."))
for j in range(d):
a, b = self.lims[0][j], self.lims[1][j]
if self.scaling is not None:
a, b = self.scaling[j](a), self.scaling[j](b)
x0, sigma = (a + b) / 2, (a - b) / 2
xj0 = np.linspace(- 1., 1., n[j])[:, None]
xj = x0 + sigma * self.warping[j](xj0)
- wj = np.abs(a - b) / n[j] * np.ones(n[j])
+ wj = np.abs(a - b) / (n[j] - 1) * np.ones(n[j])
if self.scalingInv is not None:
xj = self.scalingInv[j](xj)
if j == 0:
x = xj
w = wj
xsize = n[0]
else:
x = np.concatenate((np.kron(np.ones(n[j])[:, None], x),
np.kron(xj, np.ones(xsize)[:, None])),
axis = 1)
w = np.multiply(np.kron(np.ones(n[j]), w),
np.kron(wj, np.ones(xsize)))
xsize = xsize * n[j]
return [y.flatten() for y in np.split(x, xsize)], w
+ def refine(self, n:int) -> Tuple[List[Np1D], Np1D]:
+ """
+ Apply refinement. If points are not nested, equivalent to
+ generatePoints([2 * x - 1 for x in n]).
+ """
+ super().generatePoints(n)
+ assert len(self.warping) == len(self.lims[0])
+ d = len(self.lims[0])
+ try:
+ len(n)
+ except:
+ n = np.array([n])
+ if len(n) != d:
+ raise Exception(("Numbers of points must have same dimension as"
+ "limits."))
+ for j in range(d):
+ a, b = self.lims[0][j], self.lims[1][j]
+ if self.scaling is not None:
+ a, b = self.scaling[j](a), self.scaling[j](b)
+ x0, sigma = (a + b) / 2, (a - b) / 2
+ xj0 = np.linspace(1. / (n[j] - 1) - 1., 1. - 1. / (n[j] - 1),
+ n[j] - 1)[:, None]
+ xj = x0 + sigma * self.warping[j](xj0)
+ wj = np.abs(a - b) / (n[j] - 2) * np.ones(n[j] - 1)
+ if self.scalingInv is not None:
+ xj = self.scalingInv[j](xj)
+ if j == 0:
+ x = xj
+ w = wj
+ xsize = n[0] - 1
+ else:
+ x = np.concatenate((np.kron(np.ones(n[j] - 1)[:, None], x),
+ np.kron(xj, np.ones(xsize)[:, None])),
+ axis = 1)
+ w = np.multiply(np.kron(np.ones(n[j] - 1), w),
+ np.kron(wj, np.ones(xsize)))
+ xsize = xsize * (n[j] - 1)
+ return [y.flatten() for y in np.split(x, xsize)], w
+