diff --git a/examples/scipy/LagrangeSweep.py b/examples/scipy/LagrangeSweep.py
index 7194251..d1d09c9 100644
--- a/examples/scipy/LagrangeSweep.py
+++ b/examples/scipy/LagrangeSweep.py
@@ -1,85 +1,102 @@
from copy import copy
import numpy as np
-from rrompy.hfengines.scipy import HelmholtzSquareBubbleProblemEngine as HSBPE
-from rrompy.hfengines.scipy import HelmholtzBoxScatteringProblemEngine as HBSPE
+#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
-testNo = 1
npoints = 100
homog = True
homog = False
+LSratio = 2. / 3.
+sampling = "Uniform"
+#sampling = "Cheby"
+#sampling = "SinC"
-if testNo == 1:
- ks = np.power([24 + 5.j, 34 + 5.j], .5)
- solver = HSBPE(kappa = 32 ** .5, theta = np.pi / 3, n = 15)
- solver.omega = np.real(np.mean(ks))
- mutars = np.linspace(21**.5, 42**.5, npoints)
- filenamebase = '../data/output/HelmholtzBubbleSHIFTLSweep'
- homog = False
-elif testNo == 2:
- ks = [10 + 0.j, 14 + 0.j]
- solver = HBSPE(R = 5, kappa = 3, theta = - np.pi * 75 / 180, n = 10)
- solver.omega = np.real(np.mean(ks))
- mutars = np.linspace(9, 15, npoints)
- homogMSG = ""
- if not homog: homogMSG = "Non"
- filenamebase = '../data/output/HelmholtzBoxLSweep' + homogMSG + "Homog"
+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
k0 = np.mean(ks)
-shift = 6
-nsets = 5
-stride = 3
+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)
-paramsPade = {'S':Smax, 'POD':True,
- 'sampler':QS(ks, "CHEBYSHEV", rescaling, rescalingInv)}
+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': stride * i + shift + 1,
- 'M': stride * i + shift + 1,
+ 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': stride * i + shift + 2,
+ paramsSetsRB[i] = {'R': np.int(LSratio * (stride * i + shift + 1)),
'S': stride * i + shift + 2}
paramsSetsPoly[i] = {'N': 0,
- 'M': stride * i + shift + 1,
+ 'M': np.int(LSratio * (stride * i + shift + 1)),
'S': stride * i + shift + 2}
appPade = Pade(solver, mu0 = k0, approxParameters = paramsPade,
verbosity = verb, homogeneize = homog)
appRB = RB(solver, mu0 = k0, approxParameters = paramsRB,
verbosity = verb, homogeneize = homog)
appPoly = Pade(solver, mu0 = k0, approxParameters = paramsPoly,
verbosity = verb, homogeneize = 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',
- saveFormat = "png", labels = ["Pade'", "RB", "Poly"])
+ 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',
- saveFormat = "png", labels = ["Pade'", "RB", "Poly"])
+ 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',
- saveFormat = "png", labels = ["Pade'", "RB", "Poly"])
-
+ ylims = {'top' : 1e1}, saveFormat = "png",
+ labels = ["Pade'", "RB", "Poly"],
+# figsize = (5, 3.75))
+ figsize = (10, 7.5))
diff --git a/examples/scipy/PadeLagrange.py b/examples/scipy/PadeLagrange.py
index 6fdcec3..e3f330f 100644
--- a/examples/scipy/PadeLagrange.py
+++ b/examples/scipy/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
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)
approx.setupApprox()
- approx.plotSamples()
+# 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)
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/scipy/PadeTaylor.py
index 7bae70c..4fe5a30 100644
--- a/examples/scipy/PadeTaylor.py
+++ b/examples/scipy/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 = 4
+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)
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)
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/rrompy/hfengines/base/boundary_conditions.py b/rrompy/hfengines/base/boundary_conditions.py
index 48bebb9..f92ae87 100644
--- a/rrompy/hfengines/base/boundary_conditions.py
+++ b/rrompy/hfengines/base/boundary_conditions.py
@@ -1,110 +1,113 @@
# 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 <http://www.gnu.org/licenses/>.
#
from rrompy.utilities.base.types import GenExpr
from rrompy.utilities.base.fenics import bdrTrue, bdrFalse
__all__ = ['BoundaryConditions']
class BoundaryConditions:
"""Boundary conditions manager."""
allowedKinds = ["Dirichlet", "Neumann", "Robin"]
def __init__(self, kind : str = None):
if kind is None: return
kind = kind[0].upper() + kind[1:].lower()
if kind in self.allowedKinds:
getattr(self.__class__, kind + "Boundary", None).fset(self, "ALL")
else:
raise Exception("BC kind not recognized.")
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 _generalManagement(self, kind:str, value:GenExpr):
if isinstance(value, (str,)):
value = value.upper()
if value.upper() == "ALL":
self._complementaryManagementAll(kind)
elif value.upper() == "REST":
self._complementaryManagementRest(kind)
else:
raise Exception("Wildcard not recognized.")
elif callable(value):
self._standardManagement(kind, value)
else:
raise Exception(kind + "Boundary type not recognized.")
def _complementaryManagementAll(self, kind:str):
if kind not in self.allowedKinds:
raise Exception("BC kind not recognized.")
for k in self.allowedKinds:
if k != kind:
getattr(self.__class__, k + "Boundary").fset(self, bdrFalse)
super().__setattr__("_" + kind + "Boundary", bdrTrue)
if hasattr(self, "_" + kind + "Rest"):
super().__delattr__("_" + kind + "Rest")
def _complementaryManagementRest(self, kind:str):
if kind not in self.allowedKinds:
raise Exception("BC kind not recognized.")
otherBCs = []
for k in self.allowedKinds:
if k != kind:
if hasattr(self, "_" + k + "Rest"):
raise Exception("Only 1 'REST' wildcard can be specified.")
otherBCs += [getattr(self, k + "Boundary")]
def restCall(x, on_boundary):
return (on_boundary
and not any([bc(x, on_boundary) for bc in otherBCs]))
super().__setattr__("_" + kind + "Boundary", restCall)
super().__setattr__("_" + kind + "Rest", 1)
def _standardManagement(self, kind:str, bc:callable):
super().__setattr__("_" + kind + "Boundary", bc)
if hasattr(self, "_" + kind + "Rest"):
super().__delattr__("_" + kind + "Rest")
@property
def DirichletBoundary(self):
"""Function handle to DirichletBoundary."""
return self._DirichletBoundary
@DirichletBoundary.setter
def DirichletBoundary(self, DirichletBoundary):
self._generalManagement("Dirichlet", DirichletBoundary)
@property
def NeumannBoundary(self):
"""Function handle to NeumannBoundary."""
return self._NeumannBoundary
@NeumannBoundary.setter
def NeumannBoundary(self, NeumannBoundary):
self._generalManagement("Neumann", NeumannBoundary)
@property
def RobinBoundary(self):
"""Function handle to RobinBoundary."""
return self._RobinBoundary
@RobinBoundary.setter
def RobinBoundary(self, RobinBoundary):
self._generalManagement("Robin", RobinBoundary)
diff --git a/rrompy/hfengines/base/problem_engine_base.py b/rrompy/hfengines/base/problem_engine_base.py
index fd909ab..3767bcd 100644
--- a/rrompy/hfengines/base/problem_engine_base.py
+++ b/rrompy/hfengines/base/problem_engine_base.py
@@ -1,333 +1,336 @@
# 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 <http://www.gnu.org/licenses/>.
#
from abc import abstractmethod
import fenics as fen
import numpy as np
from scipy.sparse import csr_matrix
import scipy.sparse as scsp
import scipy.sparse.linalg as scspla
from matplotlib import pyplot as plt
from rrompy.utilities.base.types import (Np1D, ScOp, strLst, FenFunc,
Tuple, List)
from rrompy.utilities.base import purgeList, getNewFilename, verbosityDepth
from .boundary_conditions import BoundaryConditions
__all__ = ['ProblemEngineBase']
class ProblemEngineBase:
"""
Generic solver for parametric problems.
"""
nAs, nbs = 1, 1
functional = lambda self, u: 0.
def __init__(self, degree_threshold : int = np.inf, verbosity : int = 10):
self.BCManager = BoundaryConditions("Dirichlet")
self.V = fen.FunctionSpace(fen.UnitSquareMesh(10, 10), "P", 1)
self.verbosity = verbosity
self.resetAs()
self.resetbs()
self.bsmu = np.nan
self.liftDirichletDatamu = np.nan
self.mu0BC = np.nan
self.degree_threshold = degree_threshold
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 __dir_base__(self):
return [x for x in self.__dir__() if x[:2] != "__"]
def innerProduct(self, u:Np1D, v:Np1D) -> float:
"""Hilbert space scalar product."""
if not hasattr(self, "energyNormMatrix"):
self.buildEnergyNormForm()
return v.conj().T.dot(self.energyNormMatrix.dot(u))
def buildEnergyNormForm(self):
"""
Build sparse matrix (in CSR format) representative of scalar product.
"""
if self.verbosity >= 20:
verbosityDepth("INIT", "Assembling energy matrix.", end = "")
normMat = fen.assemble(fen.dot(self.u, self.v) * fen.dx)
normMatFen = fen.as_backend_type(normMat).mat()
self.energyNormMatrix = csr_matrix(normMatFen.getValuesCSR()[::-1],
shape = normMat.size)
if self.verbosity >= 20:
verbosityDepth("DEL", " Done.", inline = True)
def norm(self, u:Np1D) -> float:
return np.abs(self.innerProduct(u, u)) ** .5
def rescaling(self, x:Np1D) -> Np1D:
"""Rescaling in parameter dependence."""
return x
def rescalingInv(self, x:Np1D) -> Np1D:
"""Inverse rescaling in parameter dependence."""
return x
def checkAInBounds(self, der : int = 0):
"""Check if derivative index is oob for operator of linear system."""
if der < 0 or der >= self.nAs:
d = self.V.dim()
return scsp.csr_matrix((np.zeros(0), np.zeros(0), np.zeros(d + 1)),
shape = (d, d), dtype = np.complex)
def checkbInBounds(self, der : int = 0, homogeneized : bool = False):
"""Check if derivative index is oob for RHS of linear system."""
if der < 0 or der >= max(self.nbs, self.nAs * homogeneized):
return np.zeros(self.V.dim(), dtype = np.complex)
def setDirichletDatum(self, mu:complex):
"""Set Dirichlet datum if parametric."""
if hasattr(self, "liftedDirichletDatum"):
self.liftDirichletDatamu = mu
def liftDirichletData(self, mu:complex) -> Np1D:
"""Lift Dirichlet datum."""
self.setDirichletDatum(mu)
if not np.isclose(self.liftDirichletDatamu, mu):
try:
liftRe = fen.interpolate(self.DirichletDatum[0], self.V)
except:
liftRe = fen.project(self.DirichletDatum[0], self.V)
try:
liftIm = fen.interpolate(self.DirichletDatum[1], self.V)
except:
liftIm = fen.project(self.DirichletDatum[1], self.V)
self.liftedDirichletDatum = (np.array(liftRe.vector())
+ 1.j * np.array(liftIm.vector()))
return self.liftedDirichletDatum
def resetAs(self):
"""Reset (derivatives of) operator of linear system."""
self.As = [None] * self.nAs
def resetbs(self):
"""Reset (derivatives of) RHS of linear system."""
self.bs = {True: [None] * max(self.nbs, self.nAs),
False: [None] * self.nbs}
def reduceQuadratureDegree(self, fun:FenFunc, name:str):
"""Check whether to reduce compiler parameters to degree threshold."""
if not np.isinf(self.degree_threshold):
from ufl.algorithms.estimate_degrees import (
estimate_total_polynomial_degree as ETPD)
try:
deg = ETPD(fun)
except:
return False
if deg > self.degree_threshold:
if self.verbosity >= 15:
verbosityDepth("MAIN", ("Reducing quadrature degree from "
"{} to {} for {}.").format(
deg,
self.degree_threshold,
name))
return True
return False
def iterReduceQuadratureDegree(self, funsNames:List[Tuple[FenFunc, str]]):
"""
Iterate reduceQuadratureDegree over list and define reduce compiler
parameters.
"""
if funsNames is not None:
for fun, name in funsNames:
if self.reduceQuadratureDegree(fun, name):
return {"quadrature_degree" : self.degree_threshold}
return {}
@abstractmethod
def A(self, mu:complex, der : int = 0) -> ScOp:
"""Assemble (derivative of) operator of linear system."""
Anull = self.checkAInBounds(der)
if Anull is not None: return Anull
if self.As[der] is None:
self.As[der] = 0.
return self.As[der]
@abstractmethod
def b(self, mu:complex, der : int = 0,
homogeneized : bool = False) -> Np1D:
"""Assemble (derivative of) RHS of linear system."""
bnull = self.checkbInBounds(der, homogeneized)
if bnull is not None: return bnull
if self.bs[homogeneized][der] is None:
self.bs[homogeneized][der] = 0.
return self.bs[homogeneized][der]
def affineBlocksA(self, mu : complex = 0.) -> Tuple[List[Np1D], callable]:
"""Assemble affine blocks of operator of linear system."""
def lambdas(x, j):
if j == 0: return np.ones(np.size(x))
if j in range(1, self.nAs): return np.power(self.rescaling(x)
- self.rescaling(mu), j)
raise Exception("Wrong j value.")
As = [None] * self.nAs
for j in range(self.nAs):
As[j] = self.A(mu, j)
return As, lambdas
def affineBlocksb(self, mu : complex = 0., homogeneized : bool = False)\
-> Tuple[List[Np1D], callable]:
"""Assemble affine blocks of RHS of linear system."""
def lambdas(x, j):
if j == 0: return np.ones(np.size(x))
if j in range(1, self.nbsEff): return np.power(self.rescaling(x)
- self.rescaling(mu),
j)
raise Exception("Wrong j value.")
if homogeneized:
self.nbsEff = max(self.nAs, self.nbs)
else:
self.nbsEff = self.nbs
bs = [None] * self.nbsEff
for j in range(self.nbsEff):
bs[j] = self.b(mu, j, homogeneized)
return bs, lambdas
def solve(self, mu:complex, RHS : Np1D = None,
homogeneized : bool = False) -> Np1D:
"""
Find solution of linear system.
Args:
mu: parameter value.
RHS: RHS of linear system. If None, defaults to that of parametric
system. Defaults to None.
"""
A = self.A(mu)
if RHS is None: RHS = self.b(mu, 0, homogeneized)
return scspla.spsolve(A, RHS)
def residual(self, u:Np1D, mu:complex,
homogeneized : bool = False) -> Np1D:
"""
Find residual of linear system for given approximate solution.
Args:
u: numpy complex array with function dofs. If None, set to 0.
mu: parameter value.
"""
A = self.A(mu)
RHS = self.b(mu, 0, homogeneized)
if u is None: return RHS
return RHS - A.dot(u)
def plot(self, u:Np1D, name : str = "u", save : str = None,
what : strLst = 'all', saveFormat : str = "eps",
saveDPI : int = 100, **figspecs):
"""
Do some nice plots of the complex-valued function with given dofs.
Args:
u: numpy complex array with function dofs.
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.
"""
if isinstance(what, (str,)):
if what.upper() == 'ALL':
what = ['ABS', 'PHASE', 'REAL', 'IMAG']
else:
what = [what]
what = purgeList(what, ['ABS', 'PHASE', 'REAL', 'IMAG'],
listname = self.name() + ".what", baselevel = 1)
if len(what) == 0: return
if 'figsize' not in figspecs.keys():
figspecs['figsize'] = (13. * len(what) / 4, 3)
subplotcode = 100 + len(what) * 10
plt.figure(**figspecs)
plt.jet()
if 'ABS' in what:
uAb = fen.Function(self.V)
uAb.vector()[:] = np.array(np.abs(u), dtype = float)
subplotcode = subplotcode + 1
plt.subplot(subplotcode)
p = fen.plot(uAb, title = "|{0}|".format(name))
plt.colorbar(p)
if 'PHASE' in what:
uPh = fen.Function(self.V)
uPh.vector()[:] = np.array(np.angle(u), dtype = float)
subplotcode = subplotcode + 1
plt.subplot(subplotcode)
p = fen.plot(uPh, title = "phase({0})".format(name))
plt.colorbar(p)
if 'REAL' in what:
uRe = fen.Function(self.V)
uRe.vector()[:] = np.array(np.real(u), dtype = float)
subplotcode = subplotcode + 1
plt.subplot(subplotcode)
p = fen.plot(uRe, title = "Re({0})".format(name))
plt.colorbar(p)
if 'IMAG' in what:
uIm = fen.Function(self.V)
uIm.vector()[:] = np.array(np.imag(u), dtype = float)
subplotcode = subplotcode + 1
plt.subplot(subplotcode)
p = fen.plot(uIm, title = "Im({0})".format(name))
plt.colorbar(p)
if save is not None:
save = save.strip()
plt.savefig(getNewFilename("{}_fig_".format(save), saveFormat),
format = saveFormat, dpi = saveDPI)
plt.show()
plt.close()
def plotmesh(self, name : str = "Mesh", save : str = None,
saveFormat : str = "eps", saveDPI : int = 100, **figspecs):
"""
Do a nice plot of the mesh.
Args:
u: numpy complex array with function dofs.
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.
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.
"""
plt.figure(**figspecs)
fen.plot(self.V.mesh())
if save is not None:
save = save.strip()
plt.savefig(getNewFilename("{}_msh_".format(save), saveFormat),
format = saveFormat, dpi = saveDPI)
plt.show()
plt.close()
diff --git a/rrompy/reduction_methods/base/generic_approximant.py b/rrompy/reduction_methods/base/generic_approximant.py
index 924764e..1dbdbb1 100644
--- a/rrompy/reduction_methods/base/generic_approximant.py
+++ b/rrompy/reduction_methods/base/generic_approximant.py
@@ -1,499 +1,502 @@
# 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 <http://www.gnu.org/licenses/>.
#
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,
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.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
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)
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 evalApprox(self, mu:complex):
"""
Evaluate approximant at arbitrary parameter. (ABSTRACT)
Any specialization should include something like
self.setupApprox()
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:
return self.uHF + self.HFEngine.liftDirichletData(mu)
if not self.homogeneize 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 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:
return self.uApp + self.HFEngine.liftDirichletData(mu)
if not self.homogeneize 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 (translated) 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/lagrange/approximant_lagrange_pade.py b/rrompy/reduction_methods/lagrange/approximant_lagrange_pade.py
index 438501f..870bf68 100644
--- a/rrompy/reduction_methods/lagrange/approximant_lagrange_pade.py
+++ b/rrompy/reduction_methods/lagrange/approximant_lagrange_pade.py
@@ -1,358 +1,420 @@
# 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 <http://www.gnu.org/licenses/>.
#
from copy import copy
import numpy as np
from .generic_approximant_lagrange import GenericApproximantLagrange
-from rrompy.utilities.base.types import Np1D, DictAny, HFEng
+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.
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;
- '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.
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,
verbosity : int = 10):
self._preInit()
if not hasattr(self, "parameterList"):
self.parameterList = []
- self.parameterList += ["M", "N", "robustTol"]
+ self.parameterList += ["E", "M", "N", "robustTol"]
super().__init__(HFEngine = HFEngine, mu0 = mu0,
approxParameters = approxParameters,
homogeneize = homogeneize,
verbosity = verbosity)
self._postInit()
@property
def approxParameters(self):
"""
- Value of approximant parameters. Its assignment may change M, N and S.
+ Value of approximant parameters. Its assignment may change E, 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, ["M", "N"],
+ approxParametersCopy = purgeDict(approxParameters, ["E", "M", "N"],
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 "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 = self.M
+ self.M = Mold
else:
self.M = 0
if "N" in keyList:
self.N = approxParameters["N"]
elif hasattr(self, "N"):
- self.N = 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] * 2, {0:"M", 1:"N"}
+ 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)
+
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.verbosity >= 7:
- verbosityDepth("INIT", "Starting computation of denominator.")
- while True:
- self.computeSnapshots()
- TN = np.vander(self.radiusPade(self.mus), N = self.N + 1,
- increasing = True)
+ 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.T
+ else:
+ data = self.samplingEngine.samples.T
+ RHSFull = np.empty((self.S, data.shape[1] * (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)
+ G = G[0, :].reshape((self.N + 1, data.shape[1])).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
+ else:
+ self._funcOptimumWarn = 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}
+ self.approxParameters = newParameters
+ 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.T
else:
data = self.samplingEngine.samples.T
- RHSFull = np.empty((self.S, data.shape[1] * (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.N, w = self.ws)[0, :].reshape(
- (self.N + 1, data.shape[1])).T
+ g = np.polyfit(self.radiusPade(self.mus), data,
+ deg = self.E, w = self.ws)[0, :].flatten()
if self.POD:
- if self.verbosity >= 7:
- verbosityDepth("INIT", ("Solving svd for square root "
- "of gramian matrix."),
- end = "")
- _, ev, V = np.linalg.svd(G, full_matrices = False)
- ev = ev[::-1]
- eV = V[::-1, :].conj().T
- 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)
- 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
- Snew = self.S - diff
- Mnew = min(self.M, Snew - 1)
- if Mnew == self.M:
- strM = ""
+ self._funcOptimum = np.linalg.norm(g)
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 S from {3} to {4}.")\
- .format(diff + 1, self.N, Nnew, self.S, Snew, strM))
- newParameters = {"N" : Nnew, "M" : Mnew, "S" : Snew}
- self.approxParameters = newParameters
- self.Q = eV[:, 0]
+ self._funcOptimum = self.HFEngine.norm(g)
+ self._funcOptimumWarn = False
+ self.Q = np.poly1d([1])
if self.verbosity >= 7:
- verbosityDepth("DEL", "Done computing denominator.")
verbosityDepth("INIT", "Starting computation of numerator.")
- TNQ = TN.dot(self.Q).reshape((self.S, 1))
- if self.POD:
- data = self.samplingEngine.samples.dot(
- self.samplingEngine.RPOD)
- else:
- data = self.samplingEngine.samples
- RHS = np.multiply(data.T, TNQ)
+ Qeval = self.Q(self.radiusPade(self.mus)).reshape((self.S, 1))
+ RHS = np.multiply(data, Qeval)
self.P = np.polyfit(self.radiusPade(self.mus), RHS, deg = self.M,
- w = self.ws)[::-1, :].T
+ w = self.ws).T
+ if self.POD:
+ self.P = self.samplingEngine.samples.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:complex):
+ 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 = "")
- powerlist = np.power(self.radiusPade(mu), np.arange(self.M + 1))
+ 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:complex):
+ 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 = "")
- powerlist = np.power(self.radiusPade(mu), np.arange(self.N + 1))
- q = self.Q.dot(powerlist)
+ q = self.Q(self.radiusPade(mu))
if self.verbosity >= 10:
verbosityDepth("DEL", " Done.", inline = True)
return q
def evalApprox(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))
try:
self.uApp = self.getPVal(mu) / self.getQVal(mu)
except:
self.uApp = np.empty(self.P.shape[0])
self.uApp[:] = np.nan
self.lastSolvedApp = mu
if self.verbosity >= 5:
verbosityDepth("DEL", "Done evaluating approximant.")
def getPoles(self) -> Np1D:
"""
Obtain approximant poles.
Returns:
Numpy complex vector of poles.
"""
self.setupApprox()
- return self.HFEngine.rescalingInv(np.roots(self.Q[::-1])
+ return self.HFEngine.rescalingInv(self.Q.r
+ self.HFEngine.rescaling(self.mu0))
diff --git a/rrompy/reduction_methods/taylor/approximant_taylor_pade.py b/rrompy/reduction_methods/taylor/approximant_taylor_pade.py
index 9095d2a..b454091 100644
--- a/rrompy/reduction_methods/taylor/approximant_taylor_pade.py
+++ b/rrompy/reduction_methods/taylor/approximant_taylor_pade.py
@@ -1,563 +1,580 @@
# 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 <http://www.gnu.org/licenses/>.
#
from copy import copy
import numpy as np
from .generic_approximant_taylor import GenericApproximantTaylor
from rrompy.sampling.base.pod_engine import PODEngine
-from rrompy.utilities.base.types import Np1D, Np2D, Tuple, DictAny
+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,
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,
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)
+
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 True:
- if self.N == 0:
- if self.verbosity >= 7:
- verbosityDepth("DEL",
- "Done computing denominator.")
- self.Q = np.ones(1)
- break
+ 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 = ""
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.approxParameters = newParameters
- self.Q = eV[::-1, 0]
+ 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.Q = np.ones(1)
+ 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.")
- QToeplitz = np.zeros((self.Emax + 1, self.M + 1),
+ QHankel = np.zeros((self.Emax + 1, self.M + 1),
dtype = np.complex)
- for i in range(len(self.Q)):
+ for i in range(self.Q.order):
rng = np.arange(self.M + 1 - i)
- QToeplitz[rng, rng + i] = self.Q[i]
+ QHankel[rng, - 1 - rng - i] = self.Q[i]
if self.sampleType == "ARNOLDI":
- QToeplitz = self.samplingEngine.RArnoldi.dot(QToeplitz)
- self.P = self.samplingEngine.samples.dot(QToeplitz)
+ QHankel = self.samplingEngine.RArnoldi.dot(QHankel)
+ self.P = self.samplingEngine.samples.dot(QHankel)
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:complex):
+ 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 = "")
- powerlist = np.power(self.radiusPade(mu), np.arange(self.M + 1))
+ 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:complex):
+ 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 = "")
- powerlist = np.power(self.radiusPade(mu), np.arange(self.N + 1))
- q = self.Q.dot(powerlist)
+ q = self.Q(self.radiusPade(mu))
if self.verbosity >= 10:
verbosityDepth("DEL", " Done.", inline = True)
return q
def evalApprox(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))
- try:
- self.uApp = self.getPVal(mu) / self.getQVal(mu)
- except:
- self.uApp = np.empty(self.P.shape[0])
- self.uApp[:] = np.nan
+ self.uApp = self.getPVal(mu) / self.getQVal(mu)
self.lastSolvedApp = mu
if self.verbosity >= 5:
verbosityDepth("DEL", "Done evaluating approximant.")
def getPoles(self) -> Np1D:
"""
Obtain approximant poles.
Returns:
Numpy complex vector of poles.
"""
self.setupApprox()
- return self.HFEngine.rescalingInv(np.roots(self.Q[::-1])
+ return self.HFEngine.rescalingInv(self.Q.r
+ self.HFEngine.rescaling(self.mu0))
diff --git a/rrompy/sampling/base/pod_engine.py b/rrompy/sampling/base/pod_engine.py
index a376ea7..da00bf1 100644
--- a/rrompy/sampling/base/pod_engine.py
+++ b/rrompy/sampling/base/pod_engine.py
@@ -1,147 +1,150 @@
# 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 <http://www.gnu.org/licenses/>.
#
import numpy as np
from copy import copy
from rrompy.utilities.base.types import Np1D, Np2D, Tuple, HFEng
__all__ = ['PODEngine']
class PODEngine:
"""
POD engine for general matrix orthogonalization.
"""
def __init__(self, HFEngine:HFEng):
self.HFEngine = HFEngine
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 norm(self, a:Np1D) -> float:
"""Compute norm of a Hilbert space object."""
pass
def GS(self, a:Np1D, Q:Np2D, n : int = None,
aA:Np1D = None, QA:Np2D = None) -> Tuple[Np1D, Np1D, Np1D]:
"""
Compute 1 Gram-Schmidt step with given projector.
Args:
a: vector to be projected;
Q: orthogonal projection matrix;
n: number of columns of Q to be considered;
aA: augmented components of vector to be projected;
QA: augmented components of projection matrix.
Returns:
Resulting normalized vector, coefficients of a wrt the updated
basis.
"""
if n is None:
n = Q.shape[1]
if aA is None != QA is None:
raise Exception(("Either both or none of augmented components "
"must be provided."))
r = np.zeros((n + 1,), dtype = a.dtype)
if n > 0:
Q = Q[:, : n]
for j in range(2): # twice is enough!
nu = self.HFEngine.innerProduct(a, Q)
a = a - Q.dot(nu)
if aA is not None:
aA = aA - QA.dot(nu)
r[:-1] = r[:-1] + nu
r[-1] = self.HFEngine.norm(a)
if np.isclose(np.abs(r[-1]), 0.):
r[-1] = 1.
a = a / r[-1]
if aA is not None:
aA = aA / r[-1]
return a, r, aA
def QRGramSchmidt(self, A:Np2D,
only_R : bool = False) -> Tuple[Np1D, Np1D]:
"""
Compute QR decomposition of a matrix through Gram-Schmidt method.
Args:
A: matrix to be decomposed;
only_R(optional): whether to skip reconstruction of Q; defaults to
False.
Returns:
Resulting orthogonal and upper-triangular factors.
"""
N = A.shape[1]
Q = np.zeros_like(A, dtype = A.dtype)
R = np.zeros((N, N), dtype = A.dtype)
for k in range(N):
Q[:, k], R[: k + 1, k], _ = self.GS(A[:, k], Q, k)
if only_R:
return R
return Q, R
def QRHouseholder(self, A:Np2D, Q0 : Np2D = None,
only_R : bool = False) -> Tuple[Np1D, Np1D]:
"""
Compute QR decomposition of a matrix through Householder method.
Args:
A: matrix to be decomposed;
Q0(optional): initial orthogonal guess for Q; defaults to random;
only_R(optional): whether to skip reconstruction of Q; defaults to
False.
Returns:
Resulting (orthogonal and )upper-triangular factor(s).
"""
B = copy(A)
N = B.shape[1]
V = np.zeros_like(B, dtype = B.dtype)
R = np.zeros((N, N), dtype = B.dtype)
if Q0 is None:
Q = np.zeros_like(B, dtype = B.dtype) + np.random.randn(*(B.shape))
else:
Q = copy(Q0)
for k in range(N):
if Q0 is None:
Q[:, k], _, _ = self.GS(Q[:, k], Q, k)
a = B[:, k]
R[k, k] = self.HFEngine.norm(a)
alpha = self.HFEngine.innerProduct(a, Q[:, k])
if np.isclose(np.abs(alpha), 0.): s = 1.
else: s = - alpha / np.abs(alpha)
Q[:, k] = s * Q[:, k]
V[:, k], _, _ = self.GS(R[k, k] * Q[:, k] - a, Q, k)
J = np.arange(k + 1, N)
vtB = self.HFEngine.innerProduct(B[:, J], V[:, k])
B[:, J] = B[:, J] - 2 * np.outer(V[:, k], vtB)
R[k, J] = self.HFEngine.innerProduct(B[:, J], Q[:, k])
B[:, J] = B[:, J] - np.outer(Q[:, k], R[k, J])
if only_R:
return R
for k in range(N - 1, -1, -1):
J = np.arange(k, N)
vtQ = self.HFEngine.innerProduct(Q[:, J], V[:, k])
Q[:, J] = Q[:, J] - 2 * np.outer(V[:, k], vtQ)
return Q, R
diff --git a/rrompy/sampling/base/sampling_engine_base.py b/rrompy/sampling/base/sampling_engine_base.py
index b5c6066..9a537f5 100644
--- a/rrompy/sampling/base/sampling_engine_base.py
+++ b/rrompy/sampling/base/sampling_engine_base.py
@@ -1,102 +1,105 @@
# 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 <http://www.gnu.org/licenses/>.
#
from rrompy.utilities.base.types import Np1D, HFEng, strLst
from rrompy.utilities.base import verbosityDepth
__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
@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_lagrange.py b/rrompy/sampling/scipy/sampling_engine_lagrange.py
index 81bf598..8ff48f2 100644
--- a/rrompy/sampling/scipy/sampling_engine_lagrange.py
+++ b/rrompy/sampling/scipy/sampling_engine_lagrange.py
@@ -1,69 +1,69 @@
# 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 <http://www.gnu.org/licenses/>.
#
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:
ns = self.nsamples
u = self.postprocessu(self.solveLS(mu, homogeneized = homogeneize),
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:Np1D, homogeneize : bool = False) -> Np2D:
+ def iterSample(self, mus:Np1D, homogeneize : bool = False) -> Np2D:
if self.verbosity >= 5:
verbosityDepth("INIT", "Starting sampling iterations.")
- n = mu.size
+ n = mus.size
if n <= 0:
raise Exception(("Number of samples must be positive."))
self.resetHistory()
- u = self.nextSample(mu[0], homogeneize = homogeneize)
+ u = self.nextSample(mus[0], homogeneize = homogeneize)
if n > 1:
self.preallocateSamples(u, n)
for j in range(1, n):
- self.nextSample(mu[j], overwrite = True,
+ self.nextSample(mus[j], overwrite = True,
homogeneize = homogeneize)
if self.verbosity >= 5:
verbosityDepth("DEL", "Finished sampling iterations.")
return self.samples
diff --git a/rrompy/utilities/parameter_sampling/__init__.py b/rrompy/utilities/parameter_sampling/__init__.py
index fe293d3..9577ccb 100644
--- a/rrompy/utilities/parameter_sampling/__init__.py
+++ b/rrompy/utilities/parameter_sampling/__init__.py
@@ -1,33 +1,36 @@
# 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 <http://www.gnu.org/licenses/>.
#
from .generic_sampler import GenericSampler
-from .quadrature_sampler import QuadratureSampler
-from .random_sampler import RandomSampler
from .fft_sampler import FFTSampler
from .manual_sampler import ManualSampler
+from .quadrature_sampler import QuadratureSampler
+from .random_sampler import RandomSampler
+from .warping_sampler import WarpingFunction, WarpingSampler
__all__ = [
'GenericSampler',
+ 'FFTSampler',
+ 'ManualSampler',
'QuadratureSampler',
'RandomSampler',
- 'FFTSampler',
- 'ManualSampler'
+ 'WarpingFunction',
+ 'WarpingSampler'
]
diff --git a/rrompy/utilities/parameter_sampling/generic_sampler.py b/rrompy/utilities/parameter_sampling/generic_sampler.py
index 4746b77..c3254f6 100644
--- a/rrompy/utilities/parameter_sampling/generic_sampler.py
+++ b/rrompy/utilities/parameter_sampling/generic_sampler.py
@@ -1,94 +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 <http://www.gnu.org/licenses/>.
#
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
diff --git a/rrompy/utilities/parameter_sampling/manual_sampler.py b/rrompy/utilities/parameter_sampling/manual_sampler.py
index 431d28a..d06ea34 100644
--- a/rrompy/utilities/parameter_sampling/manual_sampler.py
+++ b/rrompy/utilities/parameter_sampling/manual_sampler.py
@@ -1,48 +1,59 @@
# 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 <http://www.gnu.org/licenses/>.
#
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)
- return self.points[: n], np.ones(n) / size
+ 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
diff --git a/rrompy/utilities/parameter_sampling/quadrature_sampler.py b/rrompy/utilities/parameter_sampling/quadrature_sampler.py
index 13135f7..3d46842 100644
--- a/rrompy/utilities/parameter_sampling/quadrature_sampler.py
+++ b/rrompy/utilities/parameter_sampling/quadrature_sampler.py
@@ -1,90 +1,93 @@
# 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 <http://www.gnu.org/licenses/>.
#
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])
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
diff --git a/rrompy/utilities/parameter_sampling/random_sampler.py b/rrompy/utilities/parameter_sampling/random_sampler.py
index d47d8e3..94dcc57 100644
--- a/rrompy/utilities/parameter_sampling/random_sampler.py
+++ b/rrompy/utilities/parameter_sampling/random_sampler.py
@@ -1,67 +1,70 @@
# 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 <http://www.gnu.org/licenses/>.
#
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)
diff --git a/rrompy/utilities/parameter_sampling/quadrature_sampler.py b/rrompy/utilities/parameter_sampling/warping_sampler.py
similarity index 55%
copy from rrompy/utilities/parameter_sampling/quadrature_sampler.py
copy to rrompy/utilities/parameter_sampling/warping_sampler.py
index 13135f7..84f0a9d 100644
--- a/rrompy/utilities/parameter_sampling/quadrature_sampler.py
+++ b/rrompy/utilities/parameter_sampling/warping_sampler.py
@@ -1,90 +1,114 @@
# 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 <http://www.gnu.org/licenses/>.
#
import numpy as np
from rrompy.utilities.parameter_sampling.generic_sampler import GenericSampler
from rrompy.utilities.base.types import Np1D, Tuple, List
-__all__ = ['QuadratureSampler']
+__all__ = ['WarpingSampler', 'WarpingFunction']
-class QuadratureSampler(GenericSampler):
- """Generator of quadrature sample points."""
+class WarpingFunction:
+ """Wrapper for warping function."""
- allowedKinds = ["UNIFORM", "CHEBYSHEV", "GAUSSLEGENDRE"]
-
- def __init__(self, lims:Tuple[Np1D, Np1D], kind : str = "UNIFORM",
+ 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.kind = kind
+ self.warping = warping
def __str__(self) -> str:
- return "{}_{}".format(super().__str__(), self.kind)
+ return "{}_{}".format(super().__str__(), self.warping)
+
+ 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 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)
- if self.kind == "UNIFORM":
- xj = np.linspace(a, b, n[j])[:, None]
- wj = np.abs(a - b) / n[j] * 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
+ 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])
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_sweeper/parameter_sweeper.py b/rrompy/utilities/parameter_sweeper/parameter_sweeper.py
index a4d9167..5692624 100644
--- a/rrompy/utilities/parameter_sweeper/parameter_sweeper.py
+++ b/rrompy/utilities/parameter_sweeper/parameter_sweeper.py
@@ -1,489 +1,499 @@
# 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 <http://www.gnu.org/licenses/>.
#
from copy import copy
import itertools
import csv
import numpy as np
from matplotlib import pyplot as plt
from rrompy.utilities.base.types import Np1D, DictAny, List, ROMEng
from rrompy.utilities.base import purgeList, getNewFilename, verbosityDepth
from rrompy.utilities.warning_manager import warn
__all__ = ['ParameterSweeper']
def C2R2csv(x):
x = np.ravel(x)
y = np.concatenate((np.real(x), np.imag(x)))
z = np.ravel(np.reshape(y, [2, np.size(x)]).T)
return np.array2string(z, separator = '_', suppress_small = False,
max_line_width = np.inf, sign = '+',
formatter = {'all' : lambda x : "{:.15E}".format(x)}
)[1 : -1]
class ParameterSweeper:
"""
ROM approximant parameter sweeper.
Args:
ROMEngine(optional): Generic approximant class. Defaults to None.
mutars(optional): Array of parameter values to sweep. Defaults to empty
array.
params(optional): List of parameter settings (each as a dict) to
explore. Defaults to single empty set.
mostExpensive(optional): String containing label of most expensive
step, to be executed fewer times. Allowed options are 'HF' and
'Approx'. Defaults to 'HF'.
Attributes:
ROMEngine: Generic approximant class.
mutars: Array of parameter values to sweep.
params: List of parameter settings (each as a dict) to explore.
mostExpensive: String containing label of most expensive step, to be
executed fewer times.
"""
allowedOutputsStandard = ["normHF", "normApp", "normRes", "normResRel",
"normErr", "normErrRel"]
allowedOutputs = allowedOutputsStandard + ["HFFunc", "AppFunc",
"ErrFunc", "ErrFuncRel"]
allowedOutputsFull = allowedOutputs + ["poles"]
def __init__(self, ROMEngine : ROMEng = None, mutars : Np1D = np.array([]),
params : List[DictAny] = [{}], mostExpensive : str = "HF"):
self.ROMEngine = ROMEngine
self.mutars = mutars
self.params = params
self.mostExpensive = mostExpensive
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))
+
@property
def mostExpensive(self):
"""Value of mostExpensive."""
return self._mostExpensive
@mostExpensive.setter
def mostExpensive(self, mostExpensive:str):
mostExpensive = mostExpensive.upper()
if mostExpensive not in ["HF", "APPROX"]:
warn(("Value of mostExpensive not recognized. Overriding to "
"'APPROX'."))
mostExpensive = "APPROX"
self._mostExpensive = mostExpensive
def checkValues(self) -> bool:
"""Check if sweep can be performed."""
if self.ROMEngine is None:
raise Exception("ROMEngine is missing. Aborting.")
if len(self.mutars) == 0:
raise Exception("Empty target parameter vector. Aborting.")
if len(self.params) == 0:
raise Exception("Empty method parameters vector. Aborting.")
def sweep(self, filename : str = "out.dat", outputs : List[str] = [],
verbose : int = 10):
self.checkValues()
try:
if outputs.upper() == "ALL":
outputs = self.allowedOutputsFull
except:
if len(outputs) == 0:
outputs = self.allowedOutputsStandard
outputs = purgeList(outputs, self.allowedOutputsFull,
listname = self.name() + ".outputs",
baselevel = 1)
poles = ("poles" in outputs)
if len(outputs) == 0:
raise Exception("Empty outputs. Aborting.")
outParList = self.ROMEngine.parameterList
Nparams = len(self.params)
if poles: polesCheckList = []
allowedParams = self.ROMEngine.parameterList
dotPos = filename.rfind('.')
if dotPos in [-1, len(filename) - 1]:
filename = getNewFilename(filename[:dotPos])
else:
filename = getNewFilename(filename[:dotPos], filename[dotPos + 1:])
append_write = "w"
initial_row = (outParList + ["muRe", "muIm"]
+ [x for x in self.allowedOutputs if x in outputs]
+ ["type"] + ["poles"] * poles)
with open(filename, append_write, buffering = 1) as fout:
writer = csv.writer(fout, delimiter=",")
writer.writerow(initial_row)
if self.mostExpensive == "HF":
outerSet = self.mutars
innerSet = self.params
elif self.mostExpensive == "APPROX":
outerSet = self.params
innerSet = self.mutars
for outerIdx, outerPar in enumerate(outerSet):
if self.mostExpensive == "HF":
i, mutar = outerIdx, outerPar
elif self.mostExpensive == "APPROX":
j, par = outerIdx, outerPar
self.ROMEngine.approxParameters = {k: par[k] for k in\
par.keys() & allowedParams}
self.ROMEngine.setupApprox()
for innerIdx, innerPar in enumerate(innerSet):
if self.mostExpensive == "APPROX":
i, mutar = innerIdx, innerPar
elif self.mostExpensive == "HF":
j, par = innerIdx, innerPar
self.ROMEngine.approxParameters = {k: par[k] for k in\
par.keys() & allowedParams}
self.ROMEngine.setupApprox()
if verbose >= 5:
verbosityDepth("INIT", "Set {}/{}\tmu_{} = {:.10f}"\
.format(j + 1, Nparams, i, mutar))
outData = []
if "normHF" in outputs:
valNorm = self.ROMEngine.normHF(mutar)
outData = outData + [valNorm]
if "normApp" in outputs:
val = self.ROMEngine.normApp(mutar)
outData = outData + [val]
if "normRes" in outputs:
valNRes = self.ROMEngine.normRes(mutar)
outData = outData + [valNRes]
if "normResRel" in outputs:
if "normRes" not in outputs:
valNRes = self.ROMEngine.normRes(mutar)
val = self.ROMEngine.normRHS(mutar)
outData = outData + [valNRes / val]
if "normErr" in outputs:
valNErr = self.ROMEngine.normErr(mutar)
outData = outData + [valNErr]
if "normErrRel" in outputs:
if "normHF" not in outputs:
valNorm = self.ROMEngine.normHF(mutar)
if "normErr" not in outputs:
valNErr = self.ROMEngine.normErr(mutar)
outData = outData + [valNErr / valNorm]
if "HFFunc" in outputs:
valFunc = self.ROMEngine.HFEngine.functional(
self.ROMEngine.getHF(mutar))
outData = outData + [valFunc]
if "AppFunc" in outputs:
valFApp = self.ROMEngine.HFEngine.functional(
self.ROMEngine.getApp(mutar))
outData = outData + [valFApp]
if "ErrFunc" in outputs:
if "HFFunc" not in outputs:
valFunc = self.ROMEngine.HFEngine.functional(
self.ROMEngine.getHF(mutar))
if "AppFunc" not in outputs:
valFApp = self.ROMEngine.HFEngine.functional(
self.ROMEngine.getApp(mutar))
valFErr = np.abs(valFApp - valFunc)
outData = outData + [valFErr]
if "ErrFuncRel" in outputs:
if not ("HFFunc" in outputs or "ErrFunc" in outputs):
valFunc = self.ROMEngine.HFEngine.functional(
self.ROMEngine.getHF(mutar))
if not ("AppFunc" in outputs or "ErrFunc" in outputs):
valFApp = self.ROMEngine.HFEngine.functional(
self.ROMEngine.getApp(mutar))
val = np.nan
if not np.isclose(valFunc, 0.):
val = valFApp / valFunc
outData = outData + [val]
writeData = []
for parn in outParList:
writeData = (writeData
+ [self.ROMEngine.approxParameters[parn]])
writeData = (writeData + [mutar.real, mutar.imag]
+ outData + [self.ROMEngine.name()])
if poles:
if j not in polesCheckList:
polesCheckList += [j]
writeData = writeData + [C2R2csv(
self.ROMEngine.getPoles())]
else:
writeData = writeData + [""]
writer.writerow(str(x) for x in writeData)
if verbose >= 5:
verbosityDepth("DEL", "", end = "", inline = "")
if verbose >= 5:
if self.mostExpensive == "APPROX":
out = "Set {}/{}\tdone.\n".format(j + 1, Nparams)
elif self.mostExpensive == "HF":
out = "Point mu_{} = {:.10f}\tdone.\n".format(i, mutar)
verbosityDepth("INIT", out)
verbosityDepth("DEL", "", end = "", inline = "")
self.filename = filename
return self.filename
def read(self, filename:str, restrictions : DictAny = {},
outputs : List[str] = []) -> DictAny:
"""
Execute a query on a custom format CSV.
Args:
filename: CSV filename.
restrictions(optional): Parameter configurations to output.
Defaults to empty dictionary, i.e. output all.
outputs(optional): Values to output. Defaults to empty list, i.e.
no output.
Returns:
Dictionary of desired results, with a key for each entry of
outputs, and a numpy 1D array as corresponding value.
"""
with open(filename, 'r') as f:
reader = csv.reader(f, delimiter=',')
header = next(reader)
restrIndices, outputIndices, outputData = {}, {}, {}
for key in restrictions.keys():
try:
restrIndices[key] = header.index(key)
if not isinstance(restrictions[key], list):
restrictions[key] = [restrictions[key]]
restrictions[key] = copy(restrictions[key])
except:
warn("Ignoring key {} from restrictions.".format(key))
for key in outputs:
try:
outputIndices[key] = header.index(key)
outputData[key] = np.array([])
except:
warn("Ignoring key {} from outputs.".format(key))
for row in reader:
restrTrue = True
for key in restrictions.keys():
if row[restrIndices[key]] == restrictions[key]:
continue
try:
if np.any(np.isclose(float(row[restrIndices[key]]),
[float(x) for x in restrictions[key]])):
continue
except: pass
restrTrue = False
if restrTrue:
for key in outputIndices.keys():
try:
val = row[outputIndices[key]]
val = float(val)
finally:
outputData[key] = np.append(outputData[key], val)
return outputData
def plot(self, filename:str, xs:List[str], ys:List[str], zs:List[str],
onePlot : bool = False, save : str = None,
saveFormat : str = "eps", saveDPI : int = 100, **figspecs):
"""
Perform plots from data in filename.
Args:
filename: CSV filename.
xs: Values to put on x axes.
ys: Values to put on y axes.
zs: Meta-values for constraints.
onePlot(optional): Whether to create a single figure per x.
Defaults to False.
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.
"""
if save is not None:
save = save.strip()
zsVals = self.read(filename, outputs = zs)
zs = list(zsVals.keys())
zss = None
for key in zs:
vals = np.unique(zsVals[key])
if zss is None:
zss = copy(vals)
else:
zss = list(itertools.product(zss, vals))
lzs = len(zs)
for z in zss:
if lzs <= 1:
constr = {zs[0] : z}
else:
constr = {zs[j] : z[j] for j in range(len(zs))}
data = self.read(filename, restrictions = constr, outputs = xs+ys)
if onePlot:
for x in xs:
xVals = data[x]
p = plt.figure(**figspecs)
logScale = False
for y in ys:
yVals = data[y]
label = '{} vs {} for {}'.format(y, x, constr)
if np.min(yVals) <= - np.finfo(float).eps:
plt.plot(xVals, yVals, label = label)
else:
plt.plot(xVals, yVals, label = label)
if np.log10(np.max(yVals) / np.min(yVals)) > 1.:
logScale = True
if logScale:
ax = p.get_axes()[0]
ax.set_yscale('log')
plt.legend()
plt.grid()
if save is not None:
prefix = "{}_{}_vs_{}_{}".format(save, ys, x, constr)
plt.savefig(getNewFilename(prefix, saveFormat),
format = saveFormat, dpi = saveDPI)
plt.show()
plt.close()
else:
for x, y in itertools.product(xs, ys):
xVals, yVals = data[x], data[y]
label = '{} vs {} for {}'.format(y, x, constr)
p = plt.figure(**figspecs)
if np.min(yVals) <= - np.finfo(float).eps:
plt.plot(xVals, yVals, label = label)
else:
plt.plot(xVals, yVals, label = label)
if np.log10(np.max(yVals) / np.min(yVals)) > 1.:
ax = p.get_axes()[0]
ax.set_yscale('log')
plt.legend()
plt.grid()
if save is not None:
prefix = "{}_{}_vs_{}_{}".format(save, y, x, constr)
plt.savefig(getNewFilename(prefix, saveFormat),
format = saveFormat, dpi = saveDPI)
plt.show()
plt.close()
def plotCompare(self, filenames:List[str], xs:List[str], ys:List[str],
zs:List[str], onePlot : bool = False, save : str = None,
- saveFormat : str = "eps", saveDPI : int = 100,
- labels : List[str] = None, **figspecs):
+ ylims : dict = None, saveFormat : str = "eps",
+ saveDPI : int = 100, labels : List[str] = None,
+ **figspecs):
"""
Perform plots from data in filename1 and filename2.
Args:
filenames: CSV filenames.
xs: Values to put on x axes.
ys: Values to put on y axes.
zs: Meta-values for constraints.
onePlot(optional): Whether to create a single figure per x.
Defaults to False.
save(optional): Where to save plot(s). Defaults to None, i.e. no
saving.
+ clip(optional): Custom y axis limits. If None, automatic values are
+ kept. Defaults to None.
saveFormat(optional): Format for saved plot(s). Defaults to "eps".
saveDPI(optional): DPI for saved plot(s). Defaults to 100.
labels: Label for each dataset.
figspecs(optional key args): Optional arguments for matplotlib
figure creation.
"""
nfiles = len(filenames)
if save is not None:
save = save.strip()
if labels is None:
labels = ["{}".format(j + 1) for j in range(nfiles)]
zsVals = self.read(filenames[0], outputs = zs)
zs = list(zsVals.keys())
zss = None
for key in zs:
vals = np.unique(zsVals[key])
if zss is None:
zss = copy(vals)
else:
zss = list(itertools.product(zss, vals))
lzs = len(zs)
for z in zss:
if lzs <= 1:
constr = {zs[0] : z}
else:
constr = {zs[j] : z[j] for j in range(len(zs))}
data = [None] * nfiles
for j in range(nfiles):
data[j] = self.read(filenames[j], restrictions = constr,
outputs = xs + ys)
if onePlot:
for x in xs:
xVals = [None] * nfiles
for j in range(nfiles):
try:
xVals[j] = data[j][x]
except:
pass
p = plt.figure(**figspecs)
logScale = False
for y in ys:
for j in range(nfiles):
try:
yVals = data[j][y]
except:
pass
l = '{} vs {} for {}, {}'.format(y, x, constr,
labels[j])
if np.min(yVals) <= - np.finfo(float).eps:
plt.plot(xVals[j], yVals, label = l)
else:
plt.plot(xVals[j], yVals, label = l)
if np.log10(np.max(yVals)/np.min(yVals)) > 1.:
logScale = True
if logScale:
ax = p.get_axes()[0]
ax.set_yscale('log')
+ if ylims is not None:
+ plt.ylim(**ylims)
plt.legend()
plt.grid()
if save is not None:
prefix = "{}_{}_vs_{}_{}".format(save, ys, x, constr)
plt.savefig(getNewFilename(prefix, saveFormat),
format = saveFormat, dpi = saveDPI)
plt.show()
plt.close()
else:
for x, y in itertools.product(xs, ys):
p = plt.figure(**figspecs)
logScale = False
for j in range(nfiles):
xVals, yVals = data[j][x], data[j][y]
l = '{} vs {} for {}, {}'.format(y, x, constr,
labels[j])
if np.min(yVals) <= - np.finfo(float).eps:
plt.plot(xVals, yVals, label = l)
else:
plt.plot(xVals, yVals, label = l)
if np.log10(np.max(yVals)/np.min(yVals)) > 1.:
logScale = True
if logScale:
ax = p.get_axes()[0]
ax.set_yscale('log')
+ if ylims is not None:
+ plt.ylim(**ylims)
plt.legend()
plt.grid()
if save is not None:
prefix = "{}_{}_vs_{}_{}".format(save, y, x, constr)
plt.savefig(getNewFilename(prefix, saveFormat),
format = saveFormat, dpi = saveDPI)
plt.show()
plt.close()