diff --git a/rrompy/hfengines/base/matrix_engine_base.py b/rrompy/hfengines/base/matrix_engine_base.py
index 4076a7d..b5cac5d 100644
--- a/rrompy/hfengines/base/matrix_engine_base.py
+++ b/rrompy/hfengines/base/matrix_engine_base.py
@@ -1,514 +1,519 @@
# Copyright (C) 2018 by the RROMPy authors
#
# This file is part of RROMPy.
#
# RROMPy is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# RROMPy is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with RROMPy. If not, see .
#
from abc import abstractmethod
import numpy as np
import scipy.sparse as scsp
from numbers import Number
from matplotlib import pyplot as plt
from copy import deepcopy as copy, copy as softcopy
from rrompy.utilities.base.types import (Np1D, Np2D, ScOp, strLst, TupleAny,
List, ListAny, DictAny, paramVal,
paramList, sampList)
from rrompy.utilities.base import purgeList, getNewFilename
from rrompy.utilities.expression import (expressionEvaluator, createMonomial,
createMonomialList)
from rrompy.utilities.numerical import (hashDerivativeToIdx as hashD,
solve as tsolve, dot, customPInv)
from rrompy.utilities.exception_manager import RROMPyException, RROMPyAssert
from rrompy.parameter import checkParameter, checkParameterList
from rrompy.sampling import sampleList, emptySampleList
from rrompy.solver import setupSolver
__all__ = ['MatrixEngineBase']
class MatrixEngineBase:
"""
Generic solver for parametric matrix problems.
Attributes:
verbosity: Verbosity level.
As: Scipy sparse array representation (in CSC format) of As.
bs: Numpy array representation of bs.
cs: Numpy array representation of cs.
energyNormMatrix: Scipy sparse matrix representing inner product.
energyNormDualMatrix: Scipy sparse matrix representing dual inner
product.
energyNormPartialDualMatrix: Scipy sparse matrix representing dual
inner product without duality.
"""
def __init__(self, verbosity : int = 10, timestamp : bool = True):
self.verbosity = verbosity
self.timestamp = timestamp
self._affinePoly = True
self.nAs, self.nbs = 1, 1
self._C = None
self.setSolver("SPSOLVE", {"use_umfpack" : False})
self.npar = 0
self.outputNormMatrix = 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 __dir_base__(self):
return [x for x in self.__dir__() if x[:2] != "__"]
def __deepcopy__(self, memo):
return softcopy(self)
@property
def npar(self):
"""Value of npar."""
return self._npar
@npar.setter
def npar(self, npar):
nparOld = self._npar if hasattr(self, "_npar") else -1
if npar != nparOld:
self.rescalingExp = [1.] * npar
self._npar = npar
@property
def nAs(self):
"""Value of nAs."""
return self._nAs
@nAs.setter
def nAs(self, nAs):
self._nAs = nAs
self.resetAs()
@property
def nbs(self):
"""Value of nbs."""
return self._nbs
@nbs.setter
def nbs(self, nbs):
self._nbs = nbs
self.resetbs()
@property
def C(self):
"""Value of C."""
if self._C is None: self.buildC()
return self._C
@property
def isCEye(self):
return isinstance(self.C, Number)
@property
def affinePoly(self):
return self._affinePoly
@property
def spacedim(self):
return self.bs[0].shape[0]
def checkParameter(self, mu:paramVal):
return checkParameter(mu, self.npar)
def checkParameterList(self, mu:paramList):
return checkParameterList(mu, self.npar)
def buildEnergyNormForm(self):
"""
Build sparse matrix (in CSR format) representative of scalar product.
"""
self.energyNormMatrix = 1.
def buildEnergyNormDualForm(self):
"""
Build sparse matrix (in CSR format) representative of dual scalar
product.
"""
self.energyNormDualMatrix = 1.
def buildEnergyNormPartialDualForm(self):
"""
Build sparse matrix (in CSR format) representative of dual scalar
product without duality.
"""
self.energyNormPartialDualMatrix = 1.
def innerProduct(self, u:Np2D, v:Np2D, onlyDiag : bool = False,
dual : bool = False, is_state : bool = True) -> Np2D:
"""Scalar product."""
if is_state or self.isCEye:
if dual:
if not hasattr(self, "energyNormPartialDualMatrix"):
self.buildEnergyNormPartialDualForm()
energyMat = self.energyNormPartialDualMatrix
else:
if not hasattr(self, "energyNormMatrix"):
self.buildEnergyNormForm()
energyMat = self.energyNormMatrix
else:
energyMat = self.outputNormMatrix
if not isinstance(u, (np.ndarray,)): u = u.data
if not isinstance(v, (np.ndarray,)): v = v.data
if onlyDiag:
return np.sum(dot(energyMat, u) * v.conj(), axis = 0)
return dot(dot(energyMat, u).T, v.conj()).T
def norm(self, u:Np2D, dual : bool = False,
is_state : bool = True) -> Np1D:
return np.abs(self.innerProduct(u, u, onlyDiag = True, dual = dual,
is_state = is_state)) ** .5
def checkAInBounds(self, derI : int = 0):
"""Check if derivative index is oob for operator of linear system."""
if derI < 0:
d = self.spacedim
return scsp.csr_matrix((np.zeros(0), np.zeros(0), np.zeros(d + 1)),
shape = (d, d), dtype = np.complex)
def checkbInBounds(self, derI : int = 0):
"""Check if derivative index is oob for RHS of linear system."""
if derI < 0:
return np.zeros(self.spacedim, dtype = np.complex)
def resetAs(self):
"""Reset (derivatives of) operator of linear system."""
self.setAs([None] * self.nAs)
self.setthAs([None] * self.nAs)
def resetbs(self):
"""Reset (derivatives of) RHS of linear system."""
self.setbs([None] * self.nbs)
self.setthbs([None] * self.nbs)
def getMonomialSingleWeight(self, deg:List[int]):
return createMonomial(deg, True)
def getMonomialWeights(self, n:int):
return createMonomialList(n, self.npar, True)
def setAs(self, As:List[Np2D]):
"""Assign terms of operator of linear system."""
if len(As) != self.nAs:
raise RROMPyException(("Expected number {} of terms of As not "
"matching given list length {}.").format(self.nAs,
len(As)))
self.As = [copy(A) for A in As]
def setthAs(self, thAs:List[List[TupleAny]]):
"""Assign terms of operator of linear system."""
if len(thAs) != self.nAs:
raise RROMPyException(("Expected number {} of terms of thAs not "
"matching given list length {}.").format(self.nAs,
len(thAs)))
self.thAs = copy(thAs)
def setbs(self, bs:List[Np1D]):
"""Assign terms of RHS of linear system."""
if len(bs) != self.nbs:
raise RROMPyException(("Expected number {} of terms of bs not "
"matching given list length {}.").format(self.nbs,
len(bs)))
self.bs = [copy(b) for b in bs]
def setthbs(self, thbs:List[List[TupleAny]]):
"""Assign terms of RHS of linear system."""
if len(thbs) != self.nbs:
raise RROMPyException(("Expected number {} of terms of thbs not "
"matching given list length {}.").format(self.nbs,
len(thbs)))
self.thbs = copy(thbs)
def _assembleObject(self, mu:paramVal, objs:ListAny, th:ListAny,
derI:int) -> ScOp:
"""Assemble (derivative of) object from list of derivatives."""
mu = self.checkParameter(mu)
rExp = self.rescalingExp
muE = mu ** rExp
obj = None
for j in range(len(objs)):
if len(th[j]) <= derI and th[j][-1] is not None:
raise RROMPyException(("Cannot assemble operator. Non enough "
"derivatives of theta provided."))
if len(th[j]) > derI and th[j][derI] is not None:
expr = expressionEvaluator(th[j][derI], muE)
if hasattr(expr, "__len__"):
if len(expr) > 1:
raise RROMPyException(("Size mismatch in value of "
"theta function. Only scalars "
"allowed."))
expr = expr[0]
if obj is None:
obj = expr * objs[j]
else:
obj = obj + expr * objs[j]
return obj
@abstractmethod
def buildA(self):
"""Build terms of operator of linear system."""
if self.thAs[0] is None: self.thAs = self.getMonomialWeights(self.nAs)
if self.As[0] is None:
self.As[0] = scsp.eye(self.spacedim, format = "csr")
for j in range(self.nAs):
if self.As[j] is None: self.As[j] = self.checkAInBounds(-1)
def A(self, mu : paramVal = [], der : List[int] = 0) -> ScOp:
"""
Assemble terms of operator of linear system and return it (or its
derivative) at a given parameter.
"""
derI = hashD(der) if hasattr(der, "__len__") else der
Anull = self.checkAInBounds(derI)
if Anull is not None: return Anull
self.buildA()
assembledA = self._assembleObject(mu, self.As, self.thAs, derI)
if assembledA is None: return self.checkAInBounds(-1)
return assembledA
@abstractmethod
def buildb(self):
"""Build terms of RHS of linear system."""
if self.thbs[0] is None: self.thbs = self.getMonomialWeights(self.nbs)
for j in range(self.nbs):
if self.bs[j] is None: self.bs[j] = self.checkbInBounds(-1)
def b(self, mu : paramVal = [], der : List[int] = 0) -> Np1D:
"""
Assemble terms of RHS of linear system and return it (or its
derivative) at a given parameter.
"""
derI = hashD(der) if hasattr(der, "__len__") else der
bnull = self.checkbInBounds(derI)
if bnull is not None: return bnull
self.buildb()
assembledb = self._assembleObject(mu, self.bs, self.thbs, derI)
if assembledb is None: return self.checkbInBounds(-1)
return assembledb
def buildC(self):
"""Build terms of LHS of linear system."""
if self._C is None: self._C = 1.
def applyC(self, u:sampList):
"""Apply LHS of linear system."""
return dot(self.C, u)
def applyCpInv(self, u:sampList):
"""Apply pseudoinverse of LHS of linear system."""
return dot(customPInv(self.C), u)
+ _isStateShiftZero = True
+ def stateShift(self, mu : paramVal = []) -> Np1D:
+ return np.zeros((self.spacedim, len(mu)))
+
+ _isOutputShiftZero = True
+ def outputShift(self, mu : paramVal = []) -> Np1D:
+ return self.applyC(self.stateShift(mu))
+
def setSolver(self, solverType:str, solverArgs : DictAny = {}):
"""Choose solver type and parameters."""
self._solver, self._solverArgs = setupSolver(solverType, solverArgs)
def solve(self, mu : paramList = [], RHS : sampList = None,
force_state : bool = False, verbose : bool = False) -> sampList:
"""
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.
force_state: whether to return state before multiplication by c.
Defaults to False.
verbose: whether to notify for each solution computed. Defaults to
False.
"""
if mu == []: mu = self.mu0
mu = self.checkParameterList(mu)[0]
if self.npar == 0: mu.reset((1, self.npar), mu.dtype)
if len(mu) == 0: return emptySampleList()
if RHS is None: RHS = [self.b(m) for m in mu]
RHS = sampleList(RHS)
mult = 0 if len(RHS) == 1 else 1
RROMPyAssert(mult * (len(mu) - 1) + 1, len(RHS), "Sample size")
for j in range(len(mu)):
u = tsolve(self.A(mu[j]), RHS[mult * j], self._solver,
self._solverArgs)
if force_state:
if j == 0:
sol = emptySampleList()
sol.reset((len(u), len(mu)), dtype = u.dtype)
sol[j] = u
else:
if j == 0: sol = np.empty((len(u), len(mu)), dtype = u.dtype)
sol[:, j] = u
if verbose:
print("." * (j % 5 != 4) + "*" * (j % 5 == 4), end = "")
- if not force_state: sol = sampleList(self.applyC(sol))
+ if not force_state:
+ sol = sampleList(self.applyC(sol) - self.outputShift(mu)) #FIXME
+ else:
+ sol = sampleList(sol - self.stateShift(mu)) #FIXME
if verbose: print()
return sol
def residual(self, mu : paramList = [], u : sampList = None,
- is_state : bool = False, post_c : bool = True) -> sampList:
+ post_c : bool = True) -> sampList:
"""
Find residual of linear system for given approximate solution.
Args:
mu: parameter value.
u: numpy complex array with function dofs. If None, set to 0.
- is_state: whether given u is value before multiplication by c.
- Defaults to False.
post_c: whether to post-process using c. Defaults to True.
"""
if mu == []: mu = self.mu0
mu = self.checkParameterList(mu)[0]
if self.npar == 0: mu.reset((1, self.npar), mu.dtype)
- if u is not None:
- u = sampleList(u)
- mult = 0 if len(u) == 1 else 1
- RROMPyAssert(mult * (len(mu) - 1) + 1, len(u), "Sample size")
if len(mu) == 0: return emptySampleList()
- if u is not None:
- v = u if is_state else self.applyCpInv(u.data.T)
+ v = sampleList(self.stateShift(mu))
+ if u is not None: v = v + sampleList(u) #FIXME
for j in range(len(mu)):
- r = self.b(mu[j])
- if u is not None:
- r = r - dot(self.A(mu[j]), v[mult * j])
+ r = self.b(mu[j]) - dot(self.A(mu[j]), v[j])
if post_c:
if j == 0: res = np.empty((len(r), len(mu)), dtype = r.dtype)
res[:, j] = r
else:
if j == 0:
res = emptySampleList()
res.reset((len(r), len(mu)), dtype = r.dtype)
res[j] = r
if post_c: res = sampleList(self.applyC(res))
return res
def _rayleighQuotient(self, A:Np2D, v0:Np1D, M:Np2D, sigma : float = 0.,
nIterP : int = 10, nIterR : int = 10) -> float:
nIterP = min(nIterP, len(v0) // 2)
nIterR = min(nIterR, (len(v0) + 1) // 2)
v0 /= dot(dot(M, v0).T, v0.conj()) ** .5
for j in range(nIterP):
v0 = tsolve(A - sigma * M, dot(M, v0), self._solver,
self._solverArgs)
v0 /= dot(dot(M, v0).T, v0.conj()) ** .5
l0 = dot(A.dot(v0).T, v0.conj())
for j in range(nIterR):
v0 = tsolve(A - l0 * M, dot(M, v0), self._solver, self._solverArgs)
v0 /= dot(dot(M, v0).T, v0.conj()) ** .5
l0 = dot(A.dot(v0).T, v0.conj())
if np.isnan(l0): l0 = np.finfo(float).eps
return np.abs(l0)
def stabilityFactor(self, mu : paramList = [], u : sampList = None,
nIterP : int = 10, nIterR : int = 10) -> sampList:
"""
Find stability factor of matrix of linear system using iterative
inverse power iteration- and Rayleigh quotient-based procedure.
Args:
mu: parameter values.
u: numpy complex arrays with function dofs.
nIterP: number of iterations of power method.
nIterR: number of iterations of Rayleigh quotient method.
"""
if mu == []: mu = self.mu0
mu = self.checkParameterList(mu)[0]
if self.npar == 0: mu.reset((1, self.npar), mu.dtype)
u = sampleList(u)
- mult = 0 if len(u) == 1 else 1
- RROMPyAssert(mult * (len(mu) - 1) + 1, len(u), "Sample size")
+ solShift = self.stateShift(mu)
+ if len(u) == len(mu):
+ u = u + solShift #FIXME
+ else:
+ u = sampleList(solShift) + np.tile(u.data, (1, len(mu))) #FIXME
stabFact = np.empty(len(mu), dtype = float)
if not hasattr(self, "energyNormMatrix"):
self.buildEnergyNormForm()
- v = sampleList(self.applyCInv(u.data))
for j in range(len(mu)):
- stabFact[j] = self._rayleighQuotient(self.A(mu[j]), v[mult * j],
+ stabFact[j] = self._rayleighQuotient(self.A(mu[j]), u[j],
self.energyNormMatrix,
0., nIterP, nIterR)
return stabFact
def plot(self, u:Np1D, warping : List[callable] = None, name : str = "u",
save : str = None, what : strLst = 'all',
saveFormat : str = "eps", saveDPI : int = 100, show : bool = True,
pyplotArgs : dict = {}, **figspecs) -> str:
"""
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.
show(optional): Whether to show figure. Defaults to True.
pyplotArgs(optional): Optional arguments for pyplot.
figspecs(optional key args): Optional arguments for matplotlib
figure creation.
Returns:
Output filename.
"""
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
idxs = np.arange(self.spacedim)
if warping is not None:
idxs = warping[0](np.arange(self.spacedim))
plt.figure(**figspecs)
plt.jet()
if 'ABS' in what:
subplotcode = subplotcode + 1
plt.subplot(subplotcode)
plt.plot(idxs, np.abs(u).flatten(), **pyplotArgs)
plt.title("|{0}|".format(name))
if 'PHASE' in what:
subplotcode = subplotcode + 1
plt.subplot(subplotcode)
plt.plot(idxs, np.angle(u).flatten(), **pyplotArgs)
plt.title("phase({0})".format(name))
if 'REAL' in what:
subplotcode = subplotcode + 1
plt.subplot(subplotcode)
plt.plot(idxs, np.real(u).flatten(), **pyplotArgs)
plt.title("Re({0})".format(name))
if 'IMAG' in what:
subplotcode = subplotcode + 1
plt.subplot(subplotcode)
plt.plot(idxs, np.imag(u).flatten(), **pyplotArgs)
plt.title("Im({0})".format(name))
if save is not None:
save = save.strip()
fileOut = getNewFilename("{}_fig_".format(save), saveFormat)
plt.savefig(fileOut, format = saveFormat, dpi = saveDPI)
else: fileOut = None
if show:
plt.show()
plt.close()
return fileOut
diff --git a/rrompy/reduction_methods/base/generic_approximant.py b/rrompy/reduction_methods/base/generic_approximant.py
index 6da93ce..5bca9bc 100644
--- a/rrompy/reduction_methods/base/generic_approximant.py
+++ b/rrompy/reduction_methods/base/generic_approximant.py
@@ -1,903 +1,906 @@
# Copyright (C) 2018 by the RROMPy authors
#
# This file is part of RROMPy.
#
# RROMPy is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# RROMPy is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with RROMPy. If not, see .
#
from abc import abstractmethod
import numpy as np
from itertools import product as iterprod
from copy import deepcopy as copy
from os import remove as osrm
from rrompy.sampling.standard import (SamplingEngineStandard,
SamplingEngineStandardPOD)
from rrompy.utilities.base.types import (Np1D, DictAny, HFEng, List, Tuple,
ListAny, strLst, paramVal, paramList,
sampList)
from rrompy.utilities.base import (purgeDict, verbosityManager as vbMng,
getNewFilename)
from rrompy.utilities.exception_manager import (RROMPyException, RROMPyAssert,
RROMPy_READY, RROMPy_FRAGILE)
from rrompy.utilities.base import pickleDump, pickleLoad
from rrompy.parameter import (emptyParameterList, checkParameter,
checkParameterList)
from rrompy.sampling import sampleList, emptySampleList
__all__ = ['GenericApproximant']
def addNormFieldToClass(self, fieldName):
def objFunc(self, mu:paramList, *args, **kwargs) -> Np1D:
uV = getattr(self.__class__, "get" + fieldName)(self, mu)
kwargs["is_state"] = False
val = self.HFEngine.norm(uV, *args, **kwargs)
return val
setattr(self.__class__, "norm" + fieldName, objFunc)
def addNormDualFieldToClass(self, fieldName):
def objFunc(self, mu:paramList, *args, **kwargs) -> Np1D:
uV = getattr(self.__class__, "get" + fieldName)(self, mu)
kwargs["is_state"] = True
if "dual" not in kwargs.keys(): kwargs["dual"] = True
val = self.HFEngine.norm(uV, *args, **kwargs)
return val
setattr(self.__class__, "norm" + fieldName, objFunc)
def addPlotFieldToClass(self, fieldName):
def objFunc(self, mu:paramList, *args, **kwargs):
uV = getattr(self.__class__, "get" + fieldName)(self, mu)
kwargsCopy = copy(kwargs)
filesOut = []
for j, u in enumerate(uV):
if "name" in kwargs.keys():
kwargsCopy["name"] = kwargs["name"] + str(j)
filesOut += [self.HFEngine.plot(u, *args, **kwargsCopy)]
if filesOut[0] is None: return None
return filesOut
setattr(self.__class__, "plot" + fieldName, objFunc)
def addPlotDualFieldToClass(self, fieldName):
def objFunc(self, mu:paramList, *args, **kwargs):
uV = getattr(self.__class__, "get" + fieldName)(self, mu)
kwargsCopy = copy(kwargs)
filesOut = []
for j, u in enumerate(uV):
if "name" in kwargs.keys():
kwargsCopy["name"] = kwargs["name"] + str(j)
filesOut += [self.HFEngine.plot(u, *args, **kwargsCopy)]
if filesOut[0] is None: return None
return filesOut
setattr(self.__class__, "plot" + fieldName, objFunc)
def addOutParaviewFieldToClass(self, fieldName):
def objFunc(self, mu:paramVal, *args, **kwargs):
if not hasattr(self.HFEngine, "outParaview"):
raise RROMPyException(("High fidelity engine cannot output to "
"Paraview."))
uV = getattr(self.__class__, "get" + fieldName)(self, mu)
kwargsCopy = copy(kwargs)
filesOut = []
for j, u in enumerate(uV):
if "name" in kwargs.keys():
kwargsCopy["name"] = kwargs["name"] + str(j)
filesOut += [self.HFEngine.outParaview(u, *args, **kwargsCopy)]
if filesOut[0] is None: return None
return filesOut
setattr(self.__class__, "outParaview" + fieldName, objFunc)
def addOutParaviewTimeDomainFieldToClass(self, fieldName):
def objFunc(self, mu:paramVal, *args, **kwargs):
if not hasattr(self.HFEngine, "outParaviewTimeDomain"):
raise RROMPyException(("High fidelity engine cannot output to "
"Paraview."))
uV = getattr(self.__class__, "get" + fieldName)(self, mu)
omega = args.pop(0) if len(args) > 0 else np.real(mu)
kwargsCopy = copy(kwargs)
filesOut = []
for j, u in enumerate(uV):
if "name" in kwargs.keys():
kwargsCopy["name"] = kwargs["name"] + str(j)
filesOut += [self.HFEngine.outParaviewTimeDomain(u, omega, *args,
**kwargsCopy)]
if filesOut[0] is None: return None
return filesOut
setattr(self.__class__, "outParaviewTimeDomain" + fieldName, objFunc)
def getTrainedModelClass(name):
from importlib import import_module as im
try:
return getattr(im("rrompy.reduction_methods.trained_model"), name)
except:
raise RROMPyException("Trained model name not recognized.")
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;
- 'S': total number of samples current approximant relies upon.
Defaults to empty dict.
force_state(optional): Whether to approximate state. Defaults to False.
verbosity(optional): Verbosity level. Defaults to 10.
Attributes:
HFEngine: HF problem solver.
trainedModel: Trained model evaluator.
mu0: Default parameter.
approxParameters: Dictionary containing values for main parameters of
approximant. Recognized keys are in parameterList{Soft,Critical}.
parameterListSoft: Recognized keys of soft approximant parameters:
- 'POD': whether to compute POD of snapshots.
parameterListCritical: Recognized keys of critical approximant
parameters:
- 'S': total number of samples current approximant relies upon.
force_state: Whether to approximate state.
verbosity: Verbosity level.
POD: Whether to compute POD of snapshots.
S: Number of solution snapshots over which current approximant is
based upon.
samplingEngine: Sampling engine.
uHF: High fidelity solution(s) with parameter(s) lastSolvedHF as
sampleList.
lastSolvedHF: Parameter(s) corresponding to last computed high fidelity
solution(s) as parameterList.
uApproxReduced: Reduced approximate solution(s) with parameter(s)
lastSolvedApprox as sampleList.
lastSolvedApproxReduced: Parameter(s) corresponding to last computed
reduced approximate solution(s) as parameterList.
uApprox: Approximate solution(s) with parameter(s) lastSolvedApprox as
sampleList.
lastSolvedApprox: Parameter(s) corresponding to last computed
approximate solution(s) as parameterList.
"""
__all__ += [ftype + dtype for ftype, dtype in iterprod(
["norm", "plot", "outParaview", "outParaviewTimeDomain"],
["HF", "RHS", "Approx", "Res", "Err"])]
def __init__(self, HFEngine:HFEng, mu0 : paramVal = None,
approxParameters : DictAny = {}, force_state : bool = False,
verbosity : int = 10, timestamp : bool = True):
self._preInit()
self._mode = RROMPy_READY
self.force_state = force_state
self.verbosity = verbosity
self.timestamp = timestamp
vbMng(self, "INIT",
"Initializing engine of type {}.".format(self.name()), 10)
self._HFEngine = HFEngine
self.trainedModel = None
self.lastSolvedHF = emptyParameterList()
self.uHF = emptySampleList()
self._addParametersToList(["POD"], [True], ["S"], [1])
if mu0 is None:
if hasattr(self.HFEngine, "mu0"):
self.mu0 = checkParameter(self.HFEngine.mu0)
else:
raise RROMPyException(("Center of approximation cannot be "
"inferred from HF engine. Parameter "
"required"))
else:
self.mu0 = checkParameter(mu0, self.HFEngine.npar)
self.resetSamples()
self.approxParameters = approxParameters
self._postInit()
### add norm{HF,Err} methods
"""
Compute norm of * at arbitrary parameter.
Args:
mu: Target parameter.
Returns:
Target norm of *.
"""
for objName in ["HF", "Err"]:
addNormFieldToClass(self, objName)
### add norm{RHS,Res} methods
"""
Compute norm of * at arbitrary parameter.
Args:
mu: Target parameter.
Returns:
Target norm of *.
"""
for objName in ["RHS", "Res"]:
addNormDualFieldToClass(self, objName)
### add plot{HF,Approx,Err} methods
"""
Do some nice plots of * 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.
show(optional): Whether to show figure. Defaults to True.
figspecs(optional key args): Optional arguments for matplotlib
figure creation.
"""
for objName in ["HF", "Approx", "Err"]:
addPlotFieldToClass(self, objName)
### add plot{RHS,Res} methods
"""
Do some nice plots of * 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.
show(optional): Whether to show figure. Defaults to True.
figspecs(optional key args): Optional arguments for matplotlib
figure creation.
"""
for objName in ["RHS", "Res"]:
addPlotDualFieldToClass(self, objName)
### add outParaview{HF,RHS,Approx,Res,Err} methods
"""
Output * to ParaView file.
Args:
mu: Target parameter.
name(optional): Base name to be used for data output.
filename(optional): Name of output file.
time(optional): Timestamp.
what(optional): Which plots to do. If list, can contain 'MESH',
'ABS', 'PHASE', 'REAL', 'IMAG'. If str, same plus wildcard
'ALL'. Defaults to 'ALL'.
forceNewFile(optional): Whether to create new output file.
filePW(optional): Fenics File entity (for time series).
"""
for objName in ["HF", "RHS", "Approx", "Res", "Err"]:
addOutParaviewFieldToClass(self, objName)
### add outParaviewTimeDomain{HF,RHS,Approx,Res,Err} methods
"""
Output * to ParaView file, converted to time domain.
Args:
mu: Target parameter.
omega(optional): frequency.
timeFinal(optional): final time of simulation.
periodResolution(optional): number of time steps per period.
name(optional): Base name to be used for data output.
filename(optional): Name of output file.
forceNewFile(optional): Whether to create new output file.
"""
for objName in ["HF", "RHS", "Approx", "Res", "Err"]:
addOutParaviewTimeDomainFieldToClass(self, objName)
def _preInit(self):
if not hasattr(self, "depth"): self.depth = 0
else: self.depth += 1
@property
def tModelType(self):
raise RROMPyException("No trainedModel type assigned.")
def initializeModelData(self, datadict):
from rrompy.reduction_methods.trained_model import TrainedModelData
return (TrainedModelData(datadict["mu0"], datadict.pop("projMat"),
datadict["scaleFactor"],
datadict.pop("rescalingExp")),
["mu0", "scaleFactor", "mus"])
@property
def parameterList(self):
"""Value of parameterListSoft + parameterListCritical."""
return self.parameterListSoft + self.parameterListCritical
def _addParametersToList(self, whatSoft:strLst, defaultSoft:ListAny,
whatCritical : strLst = [],
defaultCritical : ListAny = [],
toBeExcluded : strLst = []):
if not hasattr(self, "parameterToBeExcluded"):
self.parameterToBeExcluded = []
self.parameterToBeExcluded = toBeExcluded + self.parameterToBeExcluded
if not hasattr(self, "parameterListSoft"):
self.parameterListSoft = []
if not hasattr(self, "parameterDefaultSoft"):
self.parameterDefaultSoft = {}
if not hasattr(self, "parameterListCritical"):
self.parameterListCritical = []
if not hasattr(self, "parameterDefaultCritical"):
self.parameterDefaultCritical = {}
for j, what in enumerate(whatSoft):
if what not in self.parameterToBeExcluded:
self.parameterListSoft = [what] + self.parameterListSoft
self.parameterDefaultSoft[what] = defaultSoft[j]
for j, what in enumerate(whatCritical):
if what not in self.parameterToBeExcluded:
self.parameterListCritical = ([what]
+ self.parameterListCritical)
self.parameterDefaultCritical[what] = defaultCritical[j]
def _postInit(self):
if self.depth == 0:
vbMng(self, "DEL", "Done initializing.", 10)
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):
"""Setup sampling engine."""
RROMPyAssert(self._mode, message = "Cannot setup sampling engine.")
if not hasattr(self, "_POD") or self._POD is None: return
if self.POD:
SamplingEngine = SamplingEngineStandardPOD
else:
SamplingEngine = SamplingEngineStandard
self.samplingEngine = SamplingEngine(self.HFEngine,
force_state = self.force_state,
verbosity = self.verbosity)
@property
def HFEngine(self):
"""Value of HFEngine."""
return self._HFEngine
@HFEngine.setter
def HFEngine(self, HFEngine):
raise RROMPyException("Cannot change HFEngine.")
@property
def mu0(self):
"""Value of mu0."""
return self._mu0
@mu0.setter
def mu0(self, mu0):
mu0 = checkParameter(mu0)
if not hasattr(self, "_mu0") or mu0 != self.mu0:
self.resetSamples()
self._mu0 = mu0
@property
def npar(self):
"""Number of parameters."""
return self.mu0.shape[1]
@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())
for key in self.parameterListCritical:
if key in keyList:
setattr(self, "_" + key, self.parameterDefaultCritical[key])
for key in self.parameterListSoft:
if key in keyList:
setattr(self, "_" + key, self.parameterDefaultSoft[key])
fragile = False
for key in self.parameterListCritical:
if key in keyList:
val = approxParameters[key]
else:
val = getattr(self, "_" + key, None)
if val is None:
val = self.parameterDefaultCritical[key]
getattr(self.__class__, key, None).fset(self, val)
fragile = fragile or val is None
for key in self.parameterListSoft:
if key in keyList:
val = approxParameters[key]
else:
val = getattr(self, "_" + key, None)
if val is None:
val = self.parameterDefaultSoft[key]
getattr(self.__class__, key, None).fset(self, val)
if fragile:
self._mode = RROMPy_FRAGILE
@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.samplingEngine = None
self.resetSamples()
@property
def force_state(self):
"""Value of force_state."""
return self._force_state
@force_state.setter
def force_state(self, force_state):
if hasattr(self, "_force_state"): force_stateold = self.force_state
else: force_stateold = -1
self._force_state = force_state
if force_stateold != self.force_state:
self.samplingEngine = None
self.resetSamples()
@property
def S(self):
"""Value of S."""
return self._S
@S.setter
def S(self, S):
if S <= 0: raise RROMPyException("S must be positive.")
if hasattr(self, "_S") and self._S is not None: Sold = self.S
else: Sold = -1
self._S = S
self._approxParameters["S"] = self.S
if Sold != self.S: self.resetSamples()
@property
def trainedModel(self):
"""Value of trainedModel."""
return self._trainedModel
@trainedModel.setter
def trainedModel(self, trainedModel):
self._trainedModel = trainedModel
if self._trainedModel is not None:
self._trainedModel.lastSolvedApproxReduced = emptyParameterList()
self._trainedModel.lastSolvedApprox = emptyParameterList()
self.lastSolvedApproxReduced = emptyParameterList()
self.lastSolvedApprox = emptyParameterList()
self.uApproxReduced = emptySampleList()
self.uApprox = emptySampleList()
def resetSamples(self):
if hasattr(self, "samplingEngine") and self.samplingEngine is not None:
self.samplingEngine.resetHistory()
else:
self.setupSampling()
self._mode = RROMPy_READY
def plotSamples(self, warping : List[callable] = None, name : str = "u",
save : str = None, what : strLst = 'all',
saveFormat : str = "eps", saveDPI : int = 100,
show : bool = True, plotArgs : dict = {},
**figspecs) -> List[str]:
"""
Do some nice plots of the samples.
Args:
warping(optional): Domain warping functions.
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.
show(optional): Whether to show figure. Defaults to True.
plotArgs(optional): Optional arguments for fen/pyplot.
figspecs(optional key args): Optional arguments for matplotlib
figure creation.
Returns:
Output filenames.
"""
RROMPyAssert(self._mode, message = "Cannot plot samples.")
return self.samplingEngine.plotSamples(warping, name, save, what,
saveFormat, saveDPI, show,
plotArgs, **figspecs)
def outParaviewSamples(self, name : str = "u", filename : str = "out",
times : Np1D = None, what : strLst = 'all',
forceNewFile : bool = True, folders : bool = False,
filePW = None) -> List[str]:
"""
Output samples to ParaView file.
Args:
name(optional): Base name to be used for data output.
filename(optional): Name of output file.
times(optional): Timestamps.
what(optional): Which plots to do. If list, can contain 'MESH',
'ABS', 'PHASE', 'REAL', 'IMAG'. If str, same plus wildcard
'ALL'. Defaults to 'ALL'.
forceNewFile(optional): Whether to create new output file.
folders(optional): Whether to split output in folders.
filePW(optional): Fenics File entity (for time series).
Returns:
Output filenames.
"""
RROMPyAssert(self._mode, message = "Cannot output samples.")
return self.samplingEngine.outParaviewSamples(name, folders, filename,
times, what,
forceNewFile, filePW)
def outParaviewTimeDomainSamples(self, omegas : Np1D = None,
timeFinal : Np1D = None,
periodResolution : int = 20,
name : str = "u",
filename : str = "out",
forceNewFile : bool = True,
folders : bool = False) -> List[str]:
"""
Output samples to ParaView file, converted to time domain.
Args:
omegas(optional): frequencies.
timeFinal(optional): final time of simulation.
periodResolution(optional): number of time steps per period.
name(optional): Base name to be used for data output.
filename(optional): Name of output file.
forceNewFile(optional): Whether to create new output file.
folders(optional): Whether to split output in folders.
Returns:
Output filenames.
"""
RROMPyAssert(self._mode, message = "Cannot output samples.")
return self.samplingEngine.outParaviewTimeDomainSamples(
omegas, timeFinal,
periodResolution,
name, folders, filename,
forceNewFile)
def setSamples(self, samplingEngine):
"""Copy samplingEngine and samples."""
vbMng(self, "INIT", "Transfering samples.", 10)
self.samplingEngine = copy(samplingEngine)
vbMng(self, "DEL", "Done transfering samples.", 10)
def setTrainedModel(self, model):
"""Deepcopy approximation from trained model."""
if hasattr(model, "storeTrainedModel"):
verb = model.verbosity
model.verbosity = 0
fileOut = model.storeTrainedModel()
model.verbosity = verb
else:
try:
fileOut = getNewFilename("trained_model", "pkl")
pickleDump(model.data.__dict__, fileOut)
except:
raise RROMPyException(("Failed to store model data. Parameter "
"model must have either "
"storeTrainedModel or "
"data.__dict__ properties."))
self.loadTrainedModel(fileOut)
osrm(fileOut)
@abstractmethod
def setupApprox(self):
"""
Setup approximant. (ABSTRACT)
Any specialization should include something like
if self.checkComputedApprox():
return
RROMPyAssert(self._mode, message = "Cannot setup approximant.")
...
self.trainedModel = ...
self.trainedModel.data = ...
self.trainedModel.data.approxParameters = 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 self._mode == RROMPy_FRAGILE or (self.trainedModel is not None
and self.trainedModel.data.approxParameters == self.approxParameters)
def _pruneBeforeEval(self, mu:paramList, field:str, append:bool,
prune:bool) -> Tuple[paramList, Np1D, Np1D, bool]:
mu = checkParameterList(mu, self.npar)[0]
idx = np.empty(len(mu), dtype = np.int)
if prune:
jExtra = np.zeros(len(mu), dtype = bool)
muExtra = emptyParameterList()
lastSolvedMus = getattr(self, "lastSolved" + field)
if (len(mu) > 0 and len(mu) == len(lastSolvedMus)
and mu == lastSolvedMus):
idx = np.arange(len(mu), dtype = np.int)
return muExtra, jExtra, idx, True
muKeep = copy(muExtra)
for j in range(len(mu)):
jPos = lastSolvedMus.find(mu[j])
if jPos is not None:
idx[j] = jPos
muKeep.append(mu[j])
else:
jExtra[j] = True
muExtra.append(mu[j])
if len(muKeep) > 0 and not append:
lastSolvedu = getattr(self, "u" + field)
idx[~jExtra] = getattr(self.__class__, "set" + field)(self,
muKeep, lastSolvedu[idx[~jExtra]], append)
append = True
else:
jExtra = np.ones(len(mu), dtype = bool)
muExtra = mu
return muExtra, jExtra, idx, append
def _setObject(self, mu:paramList, field:str, object:sampList,
append:bool) -> List[int]:
newMus = checkParameterList(mu, self.npar)[0]
newObj = sampleList(object)
if append:
getattr(self, "lastSolved" + field).append(newMus)
getattr(self, "u" + field).append(newObj)
Ltot = len(getattr(self, "u" + field))
return list(range(Ltot - len(newObj), Ltot))
setattr(self, "lastSolved" + field, copy(newMus))
setattr(self, "u" + field, copy(newObj))
return list(range(len(getattr(self, "u" + field))))
def setHF(self, muHF:paramList, uHF:sampleList,
append : bool = False) -> List[int]:
"""Assign high fidelity solution."""
return self._setObject(muHF, "HF", uHF, append)
def evalHF(self, mu:paramList, append : bool = False,
prune : bool = True) -> List[int]:
"""
Find high fidelity solution with original parameters and arbitrary
parameter.
Args:
mu: Target parameter.
append(optional): Whether to append new HF solutions to old ones.
prune(optional): Whether to remove duplicates of already appearing
HF solutions.
"""
muExtra, jExtra, idx, append = self._pruneBeforeEval(mu, "HF", append,
prune)
if len(muExtra) > 0:
vbMng(self, "INIT", "Solving HF model for mu = {}.".format(mu),
15)
newuHFs = self.HFEngine.solve(muExtra)
vbMng(self, "DEL", "Done solving HF model.", 15)
idx[jExtra] = self.setHF(muExtra, newuHFs, append)
return list(idx)
def setApproxReduced(self, muApproxR:paramList, uApproxR:sampleList,
append : bool = False) -> List[int]:
"""Assign high fidelity solution."""
return self._setObject(muApproxR, "ApproxReduced", uApproxR, append)
def evalApproxReduced(self, mu:paramList, append : bool = False,
prune : bool = True) -> List[int]:
"""
Evaluate reduced representation of approximant at arbitrary parameter.
Args:
mu: Target parameter.
append(optional): Whether to append new HF solutions to old ones.
prune(optional): Whether to remove duplicates of already appearing
HF solutions.
"""
self.setupApprox()
muExtra, jExtra, idx, append = self._pruneBeforeEval(mu,
"ApproxReduced",
append, prune)
if len(muExtra) > 0:
newuApproxs = self.trainedModel.getApproxReduced(muExtra)
idx[jExtra] = self.setApproxReduced(muExtra, newuApproxs, append)
return list(idx)
def setApprox(self, muApprox:paramList, uApprox:sampleList,
append : bool = False) -> List[int]:
"""Assign high fidelity solution."""
return self._setObject(muApprox, "Approx", uApprox, append)
def evalApprox(self, mu:paramList, append : bool = False,
prune : bool = True) -> List[int]:
"""
Evaluate approximant at arbitrary parameter.
Args:
mu: Target parameter.
append(optional): Whether to append new HF solutions to old ones.
prune(optional): Whether to remove duplicates of already appearing
HF solutions.
"""
self.setupApprox()
muExtra, jExtra, idx, append = self._pruneBeforeEval(mu, "Approx",
append, prune)
if len(muExtra) > 0:
newuApproxs = self.trainedModel.getApprox(muExtra)
idx[jExtra] = self.setApprox(muExtra, newuApproxs, append)
return list(idx)
def getHF(self, mu:paramList, append : bool = False,
prune : bool = True) -> sampList:
"""
Get HF solution at arbitrary parameter.
Args:
mu: Target parameter.
Returns:
HFsolution.
"""
mu = checkParameterList(mu, self.npar)[0]
idx = self.evalHF(mu, append = append, prune = prune)
return self.uHF(idx)
def getRHS(self, mu:paramList) -> sampList:
"""
Get linear system RHS at arbitrary parameter.
Args:
mu: Target parameter.
Returns:
Linear system RHS.
"""
return self.HFEngine.residual(mu, None)
def getApproxReduced(self, mu:paramList, append : bool = False,
prune : bool = True) -> sampList:
"""
Get approximant at arbitrary parameter.
Args:
mu: Target parameter.
Returns:
Reduced approximant.
"""
mu = checkParameterList(mu, self.npar)[0]
idx = self.evalApproxReduced(mu, append = append, prune = prune)
return self.uApproxReduced(idx)
def getApprox(self, mu:paramList, append : bool = False,
prune : bool = True) -> sampList:
"""
Get approximant at arbitrary parameter.
Args:
mu: Target parameter.
Returns:
Approximant.
"""
mu = checkParameterList(mu, self.npar)[0]
idx = self.evalApprox(mu, append = append, prune = prune)
return self.uApprox(idx)
def getRes(self, mu:paramList) -> sampList:
"""
Get residual at arbitrary parameter.
Args:
mu: Target parameter.
Returns:
Approximant residual.
"""
- return self.HFEngine.residual(mu, self.getApprox(mu))
+ if not (self.force_state or self.HFEngine.isCEye):
+ raise RROMPyException(("Residual of solution with non-scalar C "
+ "not computable."))
+ return self.HFEngine.residual(mu, self.getApprox(mu) / self.HFEngine.C)
def getErr(self, mu:paramList, append : bool = False,
prune : bool = True) -> sampList:
"""
Get error at arbitrary parameter.
Args:
mu: Target parameter.
Returns:
Approximant error.
"""
return (self.getApprox(mu, append = append, prune =prune)
- self.getHF(mu, append = append, prune = prune))
def normApprox(self, mu:paramList) -> float:
"""
Compute norm of approximant at arbitrary parameter.
Args:
mu: Target parameter.
Returns:
Target norm of approximant.
"""
if not (self.POD and self.HFEngine.isCEye):
return self.HFEngine.norm(self.getApprox(mu), is_state = False)
return np.linalg.norm(self.C * self.getApproxReduced(mu).data,
axis = 0)
def getPoles(self, *args, **kwargs) -> Np1D:
"""
Obtain approximant poles.
Returns:
Numpy complex vector of poles.
"""
self.setupApprox()
vbMng(self, "INIT", "Computing poles of model.", 20)
poles = self.trainedModel.getPoles(*args, **kwargs)
vbMng(self, "DEL", "Done computing poles.", 20)
return poles
def storeTrainedModel(self, filenameBase : str = "trained_model",
forceNewFile : bool = True) -> str:
"""Store trained reduced model to file."""
self.setupApprox()
vbMng(self, "INIT", "Storing trained model to file.", 20)
if forceNewFile:
filename = getNewFilename(filenameBase, "pkl")
else:
filename = "{}.pkl".format(filenameBase)
pickleDump(self.trainedModel.data.__dict__, filename)
vbMng(self, "DEL", "Done storing trained model.", 20)
return filename
def loadTrainedModel(self, filename:str):
"""Load trained reduced model from file."""
vbMng(self, "INIT", "Loading pre-trained model from file.", 20)
datadict = pickleLoad(filename)
self.mu0 = datadict["mu0"]
self.scaleFactor = datadict["scaleFactor"]
self.mus = datadict["mus"]
trainedModel = self.tModelType()
trainedModel.verbosity = self.verbosity
trainedModel.timestamp = self.timestamp
data, selfkeys = self.initializeModelData(datadict)
for key in selfkeys: setattr(self, key, datadict.pop(key))
approxParameters = datadict.pop("approxParameters")
data.approxParameters = copy(approxParameters)
for apkey in data.approxParameters.keys():
self._approxParameters[apkey] = approxParameters.pop(apkey)
for key in datadict: setattr(data, key, datadict[key])
trainedModel.data = data
self.trainedModel = trainedModel
self._mode = RROMPy_FRAGILE
vbMng(self, "DEL", "Done loading pre-trained model.", 20)
diff --git a/rrompy/reduction_methods/greedy/rational_interpolant_greedy.py b/rrompy/reduction_methods/greedy/rational_interpolant_greedy.py
index d0c0270..92de7cb 100644
--- a/rrompy/reduction_methods/greedy/rational_interpolant_greedy.py
+++ b/rrompy/reduction_methods/greedy/rational_interpolant_greedy.py
@@ -1,497 +1,497 @@
# Copyright (C) 2018 by the RROMPy authors
#
# This file is part of RROMPy.
#
# RROMPy is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# RROMPy is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with RROMPy. If not, see .
#
from copy import deepcopy as copy
import numpy as np
from .generic_greedy_approximant import (GenericGreedyApproximant,
localL2Distance as lL2D)
from rrompy.utilities.poly_fitting.polynomial import (polybases, polydomcoeff,
PolynomialInterpolator as PI,
polyvanderTotal as pvT)
from rrompy.utilities.numerical import totalDegreeN, dot
from rrompy.utilities.expression import expressionEvaluator
from rrompy.reduction_methods.standard import RationalInterpolant
from rrompy.utilities.base.types import (Np1D, Tuple, DictAny, HFEng, paramVal,
paramList)
from rrompy.utilities.base import verbosityManager as vbMng
from rrompy.utilities.poly_fitting import customFit
from rrompy.utilities.exception_manager import (RROMPyWarning, RROMPyException,
RROMPyAssert)
from rrompy.parameter import checkParameterList
__all__ = ['RationalInterpolantGreedy']
class RationalInterpolantGreedy(GenericGreedyApproximant, RationalInterpolant):
"""
ROM greedy rational 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': number of starting training points;
- 'sampler': sample point generator;
- 'greedyTol': uniform error tolerance for greedy algorithm;
defaults to 1e-2;
- 'collinearityTol': collinearity tolerance for greedy algorithm;
defaults to 0.;
- 'interactive': whether to interactively terminate greedy
algorithm; defaults to False;
- 'maxIter': maximum number of greedy steps; defaults to 1e2;
- 'refinementRatio': ratio of test points to be exhausted before
test set refinement; defaults to 0.2;
- 'nTestPoints': number of test points; defaults to 5e2;
- 'trainSetGenerator': training sample points generator; defaults
to sampler;
- 'polybasis': type of basis for interpolation; defaults to
'MONOMIAL';
- 'errorEstimatorKind': kind of error estimator; available values
include 'AFFINE', 'DISCREPANCY', 'INTERPOLATORY',
'EIM_INTERPOLATORY', and 'EIM_DIAGONAL'; defaults to 'AFFINE';
- 'interpRcond': tolerance for interpolation; defaults to None;
- 'robustTol': tolerance for robust rational denominator
management; defaults to 0.
Defaults to empty dict.
force_state(optional): Whether to approximate state. Defaults and must
be True.
verbosity(optional): Verbosity level. Defaults to 10.
Attributes:
HFEngine: HF problem solver.
mu0: Default parameter.
mus: Array of snapshot parameters.
approxParameters: Dictionary containing values for main parameters of
approximant. Recognized keys are in parameterList.
parameterListSoft: Recognized keys of soft approximant parameters:
- 'POD': whether to compute POD of snapshots.
- 'greedyTol': uniform error tolerance for greedy algorithm;
- 'collinearityTol': collinearity tolerance for greedy algorithm;
- 'interactive': whether to interactively terminate greedy
algorithm;
- 'maxIter': maximum number of greedy steps;
- 'refinementRatio': ratio of test points to be exhausted before
test set refinement;
- 'nTestPoints': number of test points;
- 'trainSetGenerator': training sample points generator;
- 'errorEstimatorKind': kind of error estimator;
- 'interpRcond': tolerance for interpolation;
- 'robustTol': tolerance for robust rational denominator
management.
parameterListCritical: Recognized keys of critical approximant
parameters:
- 'S': total number of samples current approximant relies upon;
- 'sampler': sample point generator.
force_state: Whether to approximate state.
verbosity: Verbosity level.
POD: whether to compute POD of snapshots.
S: number of test points.
sampler: Sample point generator.
greedyTol: uniform error tolerance for greedy algorithm.
collinearityTol: Collinearity tolerance for greedy algorithm.
interactive: whether to interactively terminate greedy algorithm.
maxIter: maximum number of greedy steps.
refinementRatio: ratio of training points to be exhausted before
training set refinement.
nTestPoints: number of starting training points.
trainSetGenerator: training sample points generator.
robustTol: tolerance for robust rational denominator management.
errorEstimatorKind: kind of error estimator.
interpRcond: tolerance for interpolation.
robustTol: tolerance for robust rational denominator management.
muBounds: list of bounds for parameter values.
samplingEngine: Sampling engine.
estimatorNormEngine: Engine for estimator norm computation.
uHF: High fidelity solution(s) with parameter(s) lastSolvedHF as
sampleList.
lastSolvedHF: Parameter(s) corresponding to last computed high fidelity
solution(s) as parameterList.
uApproxReduced: Reduced approximate solution(s) with parameter(s)
lastSolvedApprox as sampleList.
lastSolvedApproxReduced: Parameter(s) corresponding to last computed
reduced approximate solution(s) as parameterList.
uApprox: Approximate solution(s) with parameter(s) lastSolvedApprox as
sampleList.
lastSolvedApprox: Parameter(s) corresponding to last computed
approximate solution(s) as parameterList.
"""
_allowedEstimatorKinds = ["AFFINE", "DISCREPANCY", "INTERPOLATORY",
"EIM_INTERPOLATORY", "EIM_DIAGONAL"]
def __init__(self, HFEngine:HFEng, mu0 : paramVal = None,
approxParameters : DictAny = {}, force_state : bool = True,
verbosity : int = 10, timestamp : bool = True):
if not force_state: RROMPyWarning("Overriding force_state to True.")
self._preInit()
self._addParametersToList(["errorEstimatorKind"], ["AFFINE"],
toBeExcluded = ["M", "N", "polydegreetype",
"radialDirectionalWeights",
"nNearestNeighbor"])
super().__init__(HFEngine = HFEngine, mu0 = mu0,
approxParameters = approxParameters,
force_state = True, verbosity = verbosity,
timestamp = timestamp)
self._postInit()
@property
def E(self):
"""Value of E."""
self._E = self.sampleBatchIdx - 1
return self._E
@E.setter
def E(self, E):
RROMPyWarning(("E is used just to simplify inheritance, and its value "
"cannot be changed from that of sampleBatchIdx - 1."))
@property
def polydegreetype(self):
"""Value of polydegreetype."""
return "TOTAL"
@polydegreetype.setter
def polydegreetype(self, polydegreetype):
RROMPyWarning(("polydegreetype is used just to simplify inheritance, "
"and its value cannot be changed from 'TOTAL'."))
@property
def polybasis(self):
"""Value of polybasis."""
return self._polybasis
@polybasis.setter
def polybasis(self, polybasis):
try:
polybasis = polybasis.upper().strip().replace(" ","")
if polybasis not in polybases:
raise RROMPyException("Sample type not recognized.")
self._polybasis = polybasis
except:
RROMPyWarning(("Prescribed polybasis not recognized. Overriding "
"to 'MONOMIAL'."))
self._polybasis = "MONOMIAL"
self._approxParameters["polybasis"] = self.polybasis
@property
def errorEstimatorKind(self):
"""Value of errorEstimatorKind."""
return self._errorEstimatorKind
@errorEstimatorKind.setter
def errorEstimatorKind(self, errorEstimatorKind):
errorEstimatorKind = errorEstimatorKind.upper()
if errorEstimatorKind not in self._allowedEstimatorKinds:
RROMPyWarning(("Error estimator kind not recognized. Overriding "
"to 'AFFINE'."))
errorEstimatorKind = "AFFINE"
self._errorEstimatorKind = errorEstimatorKind
self._approxParameters["errorEstimatorKind"] = self.errorEstimatorKind
def errorEstimator(self, mus:Np1D) -> Np1D:
"""Standard residual-based error estimator."""
if self.errorEstimatorKind == "AFFINE":
return super().errorEstimator(mus)
setupOK = self.setupApprox()
if not setupOK:
err = np.empty(len(mus))
err[:] = np.nan
return err
if self.errorEstimatorKind == "DIAGONAL":
return self.errorEstimatorEIM(mus)
mus = checkParameterList(mus, self.npar)[0]
muCTest = self.trainedModel.centerNormalize(mus)
vbMng(self.trainedModel, "INIT",
"Evaluating error estimator at mu = {}.".format(mus), 10)
verb = self.trainedModel.verbosity
self.trainedModel.verbosity = 0
QTest = self.trainedModel.getQVal(mus)
if self.errorEstimatorKind == "DISCREPANCY":
nAs, nbs = len(self.HFEngine.thAs), len(self.HFEngine.thbs)
muTrainEff = self.mus ** self.HFEngine.rescalingExp
muTestEff = mus ** self.HFEngine.rescalingExp
PTrain = self.trainedModel.getPVal(self.mus).data.T
QTrain = self.trainedModel.getQVal(self.mus)
PTest = self.trainedModel.getPVal(mus).data
radiusAbTrain = np.empty((self.S, nAs * self.S + nbs),
dtype = np.complex)
radiusA = np.empty((self.S, nAs, len(mus)), dtype = np.complex)
radiusb = np.empty((nbs, len(mus)), dtype = np.complex)
for j, thA in enumerate(self.HFEngine.thAs):
idxs = j * self.S + np.arange(self.S)
radiusAbTrain[:, idxs] = expressionEvaluator(thA[0],
muTrainEff, (self.S, 1)) * PTrain
radiusA[:, j] = PTest * expressionEvaluator(thA[0],
muTestEff, (len(mus),))
for j, thb in enumerate(self.HFEngine.thbs):
idx = nAs * self.S + j
radiusAbTrain[:, idx] = QTrain * expressionEvaluator(thb[0],
muTrainEff, (self.S,))
radiusb[j] = QTest * expressionEvaluator(thb[0], muTestEff,
(len(mus),))
QRHSNorm2 = self._affineResidualMatricesContraction(radiusb)
vanTrain, _, vanTrainIdxs = pvT(self._musUniqueCN, self.N,
self.polybasis0, self._derIdxs,
self._reorder)
interpPQ = customFit(vanTrain[:, vanTrainIdxs], radiusAbTrain,
rcond = self.interpRcond)
vanTest, _, vanTestIdxs = pvT(muCTest, self.N, self.polybasis0)
DradiusAb = vanTest[:, vanTestIdxs].dot(interpPQ)
radiusA = (radiusA
- DradiusAb[:, : - nbs].reshape(len(mus), -1, self.S).T)
radiusb = radiusb - DradiusAb[:, - nbs :].T
ff, Lf, LL = self._affineResidualMatricesContraction(radiusb,
radiusA)
err = np.abs((LL - 2. * np.real(Lf) + ff) / QRHSNorm2) ** .5
else: #if self.errorEstimatorKind == "INTERPOLATORY":
muCTrain = self.trainedModel.centerNormalize(self.mus)
samplingRatio = np.prod(lL2D(muCTest.data, muCTrain.data),
axis = 1) / np.abs(QTest)
self.initEstimatorNormEngine()
QTest = np.abs(QTest)
sampRCP = copy(samplingRatio)
idx_muTestSample = np.empty(self.sampleBatchSize,
dtype = int)
for j in range(self.sampleBatchSize):
k = np.argmax(sampRCP)
idx_muTestSample[j] = k
if j + 1 < self.sampleBatchSize:
musZero = self.trainedModel.centerNormalize(mus, mus[k])
sampRCP *= np.linalg.norm(musZero.data, axis = 1)
mu_muTestSample = mus[idx_muTestSample]
app_muTestSample = self.getApproxReduced(mu_muTestSample)
app_muTestSample = dot(self.samplingEngine.samples,
app_muTestSample)
resmu = self.HFEngine.residual(mu_muTestSample, app_muTestSample,
- is_state = True, post_c = False)
+ post_c = False)
RHSmu = self.HFEngine.residual(mu_muTestSample, None,
post_c = False)
ressamples = (self.estimatorNormEngine.norm(resmu)
/ self.estimatorNormEngine.norm(RHSmu))
musT = copy(self.mus)
musT.append(mu_muTestSample)
musT = self.trainedModel.centerNormalize(musT)
musC = self.trainedModel.centerNormalize(mus)
resT = np.zeros(len(musT), dtype = np.complex)
err = np.zeros(len(mus))
for l in range(len(mu_muTestSample)):
resT[len(self.mus) + l] = (ressamples[l]
* QTest[idx_muTestSample[l]])
p = PI()
wellCond, msg = p.setupByInterpolation(musT, resT, self.M + 1,
self.polybasis, self.verbosity >= 15,
True, {}, {"rcond": self.interpRcond})
err += np.abs(p(musC))
resT[len(self.mus) + l] = 0.
err /= QTest
vbMng(self, "MAIN", msg, 15)
self.trainedModel.verbosity = verb
vbMng(self.trainedModel, "DEL", "Done evaluating error estimator", 10)
return err
def errorEstimatorEIM(self, mus:Np1D,
return_max_idxs : bool = False) -> Np1D:
"""EIM-like interpolation error estimator."""
setupOK = self.setupApprox()
if not setupOK:
err = np.empty(len(mus))
err[:] = np.nan
return err
mus = checkParameterList(mus, self.npar)[0]
vbMng(self.trainedModel, "INIT",
"Evaluating error estimator at mu = {}.".format(mus), 10)
verb = self.trainedModel.verbosity
self.trainedModel.verbosity = 0
QTest = self.trainedModel.getQVal(mus)
muCTest = self.trainedModel.centerNormalize(mus)
muCTrain = self.trainedModel.centerNormalize(self.mus)
vanderTest, _, vanderTestR = pvT(muCTest, self.E, self.polybasis)
vanderTest = vanderTest[:, vanderTestR]
vanderTestNext, _, vanderTestNextR = pvT(muCTest, self.E + 1,
self.polybasis)
vanderTestNext = vanderTestNext[:, vanderTestNextR[
vanderTest.shape[1] :]]
idxsTest = np.arange(vanderTestNext.shape[1])
basis = np.zeros((len(idxsTest), 0), dtype = float)
idxMaxEst = []
err = None
while len(idxsTest) > 0:
vanderTrial, _, vanderTrialR = pvT(muCTrain, self.E,
self.polybasis)
vanderTrial = vanderTrial[:, vanderTrialR]
vanderTrialNext, _, vanderTrialNextR = pvT(muCTrain, self.E + 1,
self.polybasis)
vanderTrialNext = vanderTrialNext[:, vanderTrialNextR[
vanderTrial.shape[1] :]]
vanderTrial = np.hstack((vanderTrial,
vanderTrialNext.dot(basis).reshape(
len(vanderTrialNext),
basis.shape[1])))
valuesTrial = vanderTrialNext[:, idxsTest]
vanderTestEff = np.hstack((vanderTest,
vanderTestNext.dot(basis).reshape(
len(vanderTestNext),
basis.shape[1])))
vanderTestNextEff = vanderTestNext[:, idxsTest]
coeffTest = np.linalg.solve(vanderTrial, valuesTrial)
errTest = np.abs((vanderTestNextEff - vanderTestEff.dot(coeffTest))
/ np.expand_dims(QTest, 1))
idxMaxErr = np.unravel_index(np.argmax(errTest), errTest.shape)
idxMaxEst += [idxMaxErr[0]]
if err is None: err = np.max(errTest, axis = 1)
if not return_max_idxs: break
muCTrain.append(muCTest[idxMaxErr[0]])
basis = np.pad(basis, [(0, 0), (0, 1)], "constant")
basis[idxsTest[idxMaxErr[1]], -1] = 1.
idxsTest = np.delete(idxsTest, idxMaxErr[1])
if self.errorEstimatorKind == "EIM_DIAGONAL":
self.assembleReducedResidualBlocks(full = False)
muTestEff = mus ** self.HFEngine.rescalingExp
radiusb = np.empty((len(self.HFEngine.thbs), len(mus)),
dtype = np.complex)
for j, thb in enumerate(self.HFEngine.thbs):
radiusb[j] = expressionEvaluator(thb[0], muTestEff)
bresb = self._affineResidualMatricesContraction(radiusb)
self.assembleReducedResidualGramian(self.trainedModel.data.projMat)
pDom = (polydomcoeff(self.E, self.polybasis)
* self.trainedModel.data.P[(-1,) + (0,) * (self.npar - 1)])
LL = pDom.conj().dot(self.trainedModel.data.gramian.dot(pDom))
if not hasattr(self, "Anorm2Approx"):
if self.HFEngine.nAs > 1:
Ader = self.HFEngine.A(self.mu0,
[1] + [0] * (self.npar - 1))
try:
Adiag = self.scaleFactor[0] * Ader.diagonal()
except:
Adiag = self.scaleFactor[0] * np.diagonal(Ader)
self.Anorm2Approx = np.mean(np.abs(Adiag) ** 2.)
if (np.isclose(self.Anorm2Approx, 0.)
or self.HFEngine.nAs <= 1):
self.Anorm2Approx = 1
jOpt = np.abs(self.Anorm2Approx * LL / bresb) ** .5
err = jOpt * err
else: #if self.errorEstimatorKind == "EIM_INTERPOLATORY":
self.initEstimatorNormEngine()
mu_muTestSample = mus[idxMaxEst[0]]
app_muTestSample = self.getApproxReduced(mu_muTestSample)
app_muTestSample = dot(self.samplingEngine.samples,
app_muTestSample)
resmu = self.HFEngine.residual(mu_muTestSample, app_muTestSample,
- is_state = True, post_c = False)
+ post_c = False)
RHSmu = self.HFEngine.residual(mu_muTestSample, None,
post_c = False)
jOpt = np.abs(self.estimatorNormEngine.norm(resmu)[0]
/ err[idxMaxEst[0]]
/ self.estimatorNormEngine.norm(RHSmu)[0])
err = jOpt * err
self.trainedModel.verbosity = verb
vbMng(self.trainedModel, "DEL", "Done evaluating error estimator", 10)
if return_max_idxs: return err, idxMaxEst
return err
def getMaxErrorEstimator(self, mus:paramList) -> Tuple[Np1D, int, float]:
"""
Compute maximum of (and index of maximum of) error estimator over given
parameters.
"""
if self.errorEstimatorKind[: 4] == "EIM_":
errorEstTest, idxMaxEst = self.errorEstimatorEIM(mus, True)
else:
errorEstTest = self.errorEstimator(mus)
idxMaxEst = np.empty(self.sampleBatchSize, dtype = int)
errCP = copy(errorEstTest)
for j in range(self.sampleBatchSize):
k = np.argmax(errCP)
idxMaxEst[j] = k
if j + 1 < self.sampleBatchSize:
musZero = self.trainedModel.centerNormalize(mus, mus[k])
errCP *= np.linalg.norm(musZero.data, axis = 1)
maxEst = errorEstTest[idxMaxEst]
return errorEstTest, idxMaxEst, maxEst
def greedyNextSample(self, muidx:int, plotEst : bool = False)\
-> Tuple[Np1D, int, float, paramVal]:
"""Compute next greedy snapshot of solution map."""
RROMPyAssert(self._mode, message = "Cannot add greedy sample.")
self.sampleBatchIdx += 1
self.sampleBatchSize = totalDegreeN(self.npar - 1, self.sampleBatchIdx)
return super().greedyNextSample(muidx, plotEst)
def _preliminaryTraining(self):
"""Initialize starting snapshots of solution map."""
RROMPyAssert(self._mode, message = "Cannot start greedy algorithm.")
if self.samplingEngine.nsamples > 0:
return
S = self.S
self.sampleBatchIdx, self.sampleBatchSize, self._S = -1, 0, 0
nextBatchSize = 1
while self._S + nextBatchSize <= S:
self.sampleBatchIdx += 1
self.sampleBatchSize = nextBatchSize
self._S += self.sampleBatchSize
nextBatchSize = totalDegreeN(self.npar - 1,
self.sampleBatchIdx + 1)
super()._preliminaryTraining()
def setupApprox(self, plotEst : bool = False):
"""
Compute rational interpolant.
SVD-based robust eigenvalue management.
"""
if self.checkComputedApprox():
return True
RROMPyAssert(self._mode, message = "Cannot setup approximant.")
vbMng(self, "INIT", "Setting up {}.".format(self.name()), 5)
self.greedy(plotEst)
self._S = len(self.mus)
self._N, self._M = [self.E] * 2
pMat = self.samplingEngine.samples.data
pMatEff = dot(self.HFEngine.C, pMat) if self.force_state else pMat
if self.trainedModel is None:
self.trainedModel = self.tModelType()
self.trainedModel.verbosity = self.verbosity
self.trainedModel.timestamp = self.timestamp
datadict = {"mu0": self.mu0, "projMat": pMatEff,
"scaleFactor": self.scaleFactor,
"rescalingExp": self.HFEngine.rescalingExp}
self.trainedModel.data = self.initializeModelData(datadict)[0]
else:
self.trainedModel = self.trainedModel
self.trainedModel.data.projMat = copy(pMatEff)
self.trainedModel.data.mus = copy(self.mus)
self.trainedModel.data.mus = copy(self.mus)
self.catchInstability = True
if self.N > 0:
try:
Q = self._setupDenominator()[0]
except RROMPyException as RE:
RROMPyWarning(RE)
vbMng(self, "DEL", "Done setting up approximant.", 5)
return False
else:
Q = PI()
Q.coeffs = np.ones(1, dtype = np.complex)
Q.npar = 1
Q.polybasis = self.polybasis
self.trainedModel.data.Q = copy(Q)
try:
self.trainedModel.data.P = copy(self._setupNumerator())
except RROMPyException as RE:
RROMPyWarning(RE)
vbMng(self, "DEL", "Done setting up approximant.", 5)
return False
self.trainedModel.data.approxParameters = copy(self.approxParameters)
vbMng(self, "DEL", "Done setting up approximant.", 5)
return True
diff --git a/rrompy/sampling/pivoted/sampling_engine_pivoted.py b/rrompy/sampling/pivoted/sampling_engine_pivoted.py
index e751811..bc2012c 100644
--- a/rrompy/sampling/pivoted/sampling_engine_pivoted.py
+++ b/rrompy/sampling/pivoted/sampling_engine_pivoted.py
@@ -1,128 +1,131 @@
# Copyright (C) 2018 by the RROMPy authors
#
# This file is part of RROMPy.
#
# RROMPy is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# RROMPy is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with RROMPy. If not, see .
#
from copy import deepcopy as copy
import numpy as np
from rrompy.sampling.base.sampling_engine_base_pivoted import (
SamplingEngineBasePivoted)
from rrompy.hfengines.base import MarginalProxyEngine
from rrompy.utilities.base.types import Np1D, paramVal, paramList, sampList
from rrompy.utilities.base import verbosityManager as vbMng, freepar as fp
from rrompy.utilities.exception_manager import RROMPyException
from rrompy.utilities.numerical import nextDerivativeIndices, dot
from rrompy.parameter import checkParameter, checkParameterList
from rrompy.sampling import sampleList
__all__ = ['SamplingEnginePivoted']
class SamplingEnginePivoted(SamplingEngineBasePivoted):
"""HERE"""
def preprocesssamples(self, idxs:Np1D, j:int) -> sampList:
if self.samples[j] is None or len(self.samples[j]) == 0: return
return self.samples[j](idxs)
def setsample(self, u:sampList, j:int, overwrite : bool = False) -> Np1D:
if overwrite:
self.samples[j][self.nsamples[j]] = u
else:
if self.nsamples[j] == 0:
self.samples[j] = sampleList(u)
else:
self.samples[j].append(u)
def postprocessu(self, u:sampList, j:int, overwrite : bool = False):
self.setsample(u, j, overwrite)
def postprocessuBulk(self, j:int):
pass
def _getSampleConcurrence(self, mu:paramVal, j:int,
previous:Np1D) -> sampList:
if not (self.force_state or self.HFEngine.isCEye):
raise RROMPyException(("Derivatives of solution with non-scalar "
"C not computable."))
+ if not self.HFEngine._isStateShiftZero:
+ raise RROMPyException(("Derivatives of solution with non-zero "
+ "solution shift not computable."))
if len(previous) >= len(self._derIdxs[j]):
self._derIdxs[j] += nextDerivativeIndices(
self._derIdxs[j], self.nPivot,
len(previous) + 1 - len(self._derIdxs[j]))
derIdx = self._derIdxs[j][len(previous)]
mu = checkParameter(mu, self.nPivot)
samplesOld = self.preprocesssamples(previous, j)
RHS = self.HFEngineMarginalized.b(mu, derIdx)
for j, derP in enumerate(self._derIdxs[j][: len(previous)]):
diffP = [x - y for (x, y) in zip(derIdx, derP)]
if np.all([x >= 0 for x in diffP]):
RHS -= dot(self.HFEngineMarginalized.A(mu, diffP),
samplesOld[j])
return self.solveLS(mu, RHS = RHS)
def nextSample(self, mu:paramVal, j:int, overwrite : bool = False,
postprocess : bool = True) -> Np1D:
mu = checkParameter(mu, self.nPivot)
muidxs = self.mus[j].findall(mu[0])
if len(muidxs) > 0:
u = self._getSampleConcurrence(mu, j, np.sort(muidxs))
else:
u = self.solveLS(mu)
if postprocess:
self.postprocessu(u, j, overwrite = overwrite)
else:
self.setsample(u, j, overwrite)
if overwrite:
self.mus[j][self.nsamples[j]] = mu[0]
else:
self.mus[j].append(mu)
self.nsamples[j] += 1
return self.samples[j][self.nsamples[j] - 1]
def iterSample(self, mus:paramList, musM:paramList) -> sampList:
mus = checkParameterList(mus, self.nPivot)[0]
musM = checkParameterList(musM, self.nMarginal)[0]
vbMng(self, "INIT", "Starting sampling iterations.", 5)
n = len(mus)
m = len(musM)
if n <= 0:
raise RROMPyException("Number of samples must be positive.")
if m <= 0:
raise RROMPyException(("Number of marginal samples must be "
"positive."))
repeatedSamples = len(mus.unique()) != n
for j in range(m):
muMEff = [fp] * self.HFEngine.npar
for k, x in enumerate(self.directionMarginal):
muMEff[x] = musM(j, k)
self.HFEngineMarginalized = MarginalProxyEngine(self.HFEngine,
list(muMEff))
if repeatedSamples:
for k in range(n):
vbMng(self, "MAIN",
"Computing sample {} / {} for marginal {} / {}."\
.format(k + 1, n, j, m), 10)
self.nextSample(mus[k], j, overwrite = (k > 0),
postprocess = False)
if n > 1 and k == 0:
self.preallocateSamples(self.samples[j][0], mus[0],
n, j)
else:
self.samples[j] = self.postprocessuBulk(self.solveLS(mus), j)
self.mus[j] = copy(mus)
self.nsamples[j] = n
self.postprocessuBulk(j)
vbMng(self, "DEL", "Finished sampling iterations.", 5)
return self.samples[j]
diff --git a/rrompy/sampling/standard/sampling_engine_standard.py b/rrompy/sampling/standard/sampling_engine_standard.py
index 6dc03ae..d6c0c83 100644
--- a/rrompy/sampling/standard/sampling_engine_standard.py
+++ b/rrompy/sampling/standard/sampling_engine_standard.py
@@ -1,112 +1,115 @@
# Copyright (C) 2018 by the RROMPy authors
#
# This file is part of RROMPy.
#
# RROMPy is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# RROMPy is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with RROMPy. If not, see .
#
from copy import deepcopy as copy
import numpy as np
from rrompy.sampling.base.sampling_engine_base import SamplingEngineBase
from rrompy.utilities.base.types import Np1D, paramVal, paramList, sampList
from rrompy.utilities.base import verbosityManager as vbMng
from rrompy.utilities.exception_manager import RROMPyException
from rrompy.utilities.numerical import nextDerivativeIndices, dot
from rrompy.parameter import checkParameter, checkParameterList
from rrompy.sampling import sampleList
__all__ = ['SamplingEngineStandard']
class SamplingEngineStandard(SamplingEngineBase):
"""HERE"""
def preprocesssamples(self, idxs:Np1D) -> sampList:
if self.samples is None or len(self.samples) == 0: return
return self.samples(idxs)
def setsample(self, u:sampList, overwrite : bool = False):
if overwrite:
self.samples[self.nsamples] = u
else:
if self.nsamples == 0:
self.samples = sampleList(u)
else:
self.samples.append(u)
def postprocessu(self, u:sampList, overwrite : bool = False):
self.setsample(u, overwrite)
def postprocessuBulk(self):
pass
def _getSampleConcurrence(self, mu:paramVal, previous:Np1D) -> sampList:
if not (self.force_state or self.HFEngine.isCEye):
raise RROMPyException(("Derivatives of solution with non-scalar "
"C not computable."))
+ if not self.HFEngine._isStateShiftZero:
+ raise RROMPyException(("Derivatives of solution with non-zero "
+ "solution shift not computable."))
if len(previous) >= len(self._derIdxs):
self._derIdxs += nextDerivativeIndices(self._derIdxs,
self.HFEngine.npar,
len(previous) + 1 - len(self._derIdxs))
derIdx = self._derIdxs[len(previous)]
mu = checkParameter(mu, self.HFEngine.npar)
samplesOld = self.preprocesssamples(previous)
RHS = self.HFEngine.b(mu, derIdx)
for j, derP in enumerate(self._derIdxs[: len(previous)]):
diffP = [x - y for (x, y) in zip(derIdx, derP)]
if np.all([x >= 0 for x in diffP]):
RHS -= dot(self.HFEngine.A(mu, diffP), samplesOld[j])
return self.solveLS(mu, RHS = RHS)
def nextSample(self, mu : paramVal = [], overwrite : bool = False,
postprocess : bool = True) -> Np1D:
mu = checkParameter(mu, self.HFEngine.npar)
muidxs = self.mus.findall(mu[0])
if len(muidxs) > 0:
u = self._getSampleConcurrence(mu, np.sort(muidxs))
else:
u = self.solveLS(mu)
if postprocess:
self.postprocessu(u, overwrite = overwrite)
else:
self.setsample(u, overwrite)
if overwrite:
self.mus[self.nsamples] = mu[0]
else:
self.mus.append(mu)
self.nsamples += 1
return self.samples[self.nsamples - 1]
def iterSample(self, mus:paramList) -> sampList:
mus = checkParameterList(mus, self.HFEngine.npar)[0]
vbMng(self, "INIT", "Starting sampling iterations.", 5)
n = len(mus)
if n <= 0:
raise RROMPyException(("Number of samples must be positive."))
self.resetHistory()
if len(mus.unique()) != n:
for j in range(n):
vbMng(self, "MAIN",
"Computing sample {} / {}.".format(j + 1, n), 7)
self.nextSample(mus[j], overwrite = (j > 0),
postprocess = False)
if n > 1 and j == 0:
self.preallocateSamples(self.samples[0], mus[0], n)
else:
self.setsample(self.solveLS(mus), overwrite = False)
self.mus = copy(mus)
self.nsamples = n
self.postprocessuBulk()
vbMng(self, "DEL", "Finished sampling iterations.", 5)
return self.samples