diff --git a/rrompy/parameter/parameter_sampling/fft_sampler.py b/rrompy/parameter/parameter_sampling/fft_sampler.py
index 5ea7e20..4199886 100644
--- a/rrompy/parameter/parameter_sampling/fft_sampler.py
+++ b/rrompy/parameter/parameter_sampling/fft_sampler.py
@@ -1,50 +1,48 @@
# Copyright (C) 2018 by the RROMPy authors
#
# This file is part of RROMPy.
#
# RROMPy is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# RROMPy is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with RROMPy. If not, see .
#
import numpy as np
from .generic_sampler import GenericSampler
from rrompy.utilities.base.types import List, paramList
from rrompy.utilities.numerical import lowDiscrepancy, kroneckerer
from rrompy.parameter import checkParameterList
__all__ = ['FFTSampler']
class FFTSampler(GenericSampler):
"""Generator of FFT-type sample points on scaled roots of unity."""
def generatePoints(self, n:List[int], reorder : bool = True) -> paramList:
"""Array of sample points."""
if not hasattr(n, "__len__"): n = [n]
super().generatePoints(n)
nleft, nright = 1, np.prod(n)
xmat = np.empty((nright, self.npar), dtype = np.complex)
for d in range(self.npar):
nright //= n[d]
- a, b = self.lims(0, d), self.lims(1, d)
- if self.scaling is not None:
- a, b = self.scaling[d](a), self.scaling[d](b)
+ a = self.lims(0, d) ** self.scalingExp[d]
+ b = self.lims(1, d) ** self.scalingExp[d]
c, r = (a + b) / 2., np.abs(a - b) / 2.
xd = c + r * np.exp(1.j * np.linspace(0, 2 * np.pi, n[d] + 1)[:-1])
- if self.scalingInv is not None:
- xd = self.scalingInv[d](xd)
+ xd **= 1. / self.scalingExp[d]
xmat[:, d] = kroneckerer(xd, nleft, nright)
nleft *= n[d]
if reorder:
xmat = xmat[lowDiscrepancy(np.prod(n)), :]
x = checkParameterList(xmat, self.npar)[0]
return x
diff --git a/rrompy/parameter/parameter_sampling/generic_sampler.py b/rrompy/parameter/parameter_sampling/generic_sampler.py
index 9282e18..30f631a 100644
--- a/rrompy/parameter/parameter_sampling/generic_sampler.py
+++ b/rrompy/parameter/parameter_sampling/generic_sampler.py
@@ -1,98 +1,79 @@
# 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
from rrompy.utilities.base.types import List, paramList
from rrompy.utilities.exception_manager import RROMPyException, RROMPyAssert
from rrompy.parameter import checkParameterList
__all__ = ['GenericSampler']
class GenericSampler:
"""ABSTRACT. Generic generator of sample points."""
- def __init__(self, lims:paramList, scaling : List[callable] = None,
- scalingInv : List[callable] = None):
+ def __init__(self, lims:paramList, scalingExp : List[float] = None):
self.lims = lims
- self.scaling = scaling
- self.scalingInv = scalingInv
+ self.scalingExp = scalingExp
def name(self) -> str:
return self.__class__.__name__
def __str__(self) -> str:
return "{}[{}_{}]".format(self.name(), self.lims[0], self.lims[1])
def __repr__(self) -> str:
return self.__str__() + " at " + hex(id(self))
def __eq__(self, other) -> bool:
return self.__dict__ == other.__dict__
@property
def npar(self):
"""Number of parameters."""
return self._lims.shape[1]
@property
def lims(self):
"""Value of lims."""
return self._lims
@lims.setter
def lims(self, lims):
lims = checkParameterList(lims)[0]
if len(lims) != 2:
raise RROMPyException("2 limits must be specified.")
self._lims = lims
@property
- def scaling(self):
- """Value of scaling."""
- return self._scaling
- @scaling.setter
- def scaling(self, scaling):
- if scaling is not None:
- if not hasattr(scaling, "__len__"): scaling = [scaling]
- RROMPyAssert(self.npar, len(scaling), "Number of scaling terms")
- if not all([callable(s) for s in scaling]):
- raise RROMPyException(("Each value of scaling must be a "
- "callable."))
- self._scaling = scaling
-
- @property
- def scalingInv(self):
- """Value of scalingInv."""
- return self._scalingInv
- @scalingInv.setter
- def scalingInv(self, scalingInv):
- if scalingInv is not None:
- if not hasattr(scalingInv, "__len__"): scalingInv = [scalingInv]
- RROMPyAssert(self.npar, len(scalingInv),
- "Number of scalingInv terms")
- if not all([callable(sInv) for sInv in scalingInv]):
- raise RROMPyException(("Each value of scalingInv must be a "
- "callable."))
- self._scalingInv = scalingInv
+ def scalingExp(self):
+ """Value of scalingExp."""
+ return self._scalingExp
+ @scalingExp.setter
+ def scalingExp(self, scalingExp):
+ if scalingExp is None:
+ scalingExp = [1.] * self.npar
+ if not hasattr(scalingExp, "__len__"): scalingExp = [scalingExp]
+ RROMPyAssert(self.npar, len(scalingExp), "Number of scaling terms")
+ self._scalingExp = scalingExp
@abstractmethod
def generatePoints(self, n:List[int]) -> paramList:
"""Array of points."""
if not hasattr(n, "__len__"): n = [n]
RROMPyAssert(self.npar, len(n), "Point number")
pass
diff --git a/rrompy/parameter/parameter_sampling/manual_sampler.py b/rrompy/parameter/parameter_sampling/manual_sampler.py
index 4f4be17..42aa93d 100644
--- a/rrompy/parameter/parameter_sampling/manual_sampler.py
+++ b/rrompy/parameter/parameter_sampling/manual_sampler.py
@@ -1,63 +1,61 @@
# Copyright (C) 2018 by the RROMPy authors
#
# This file is part of RROMPy.
#
# RROMPy is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# RROMPy is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with RROMPy. If not, see .
#
import numpy as np
from copy import deepcopy as copy
from .generic_sampler import GenericSampler
from rrompy.utilities.base.types import List, paramList
from rrompy.parameter import checkParameterList
__all__ = ['ManualSampler']
class ManualSampler(GenericSampler):
"""Manual generator of sample points."""
def __init__(self, lims:paramList, points:paramList,
- scaling : List[callable] = None,
- scalingInv : List[callable] = None):
- super().__init__(lims = lims, scaling = scaling,
- scalingInv = scalingInv)
+ scalingExp : List[float] = None):
+ super().__init__(lims = lims, scalingExp = scalingExp)
self.points = points
@property
def points(self):
"""Value of points."""
return self._points
@points.setter
def points(self, points):
points = checkParameterList(points, self.npar)[0]
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) -> paramList:
"""Array of sample points."""
if hasattr(n, "__len__"): n = n[0]
if n > len(self.points):
pts = copy(self.points)
for j in range(int(np.ceil(n / len(self.points)))):
pts.append(self.points)
else:
pts = self.points
x = checkParameterList(pts[list(range(n))], self.npar)[0]
return x
diff --git a/rrompy/parameter/parameter_sampling/quadrature_sampler.py b/rrompy/parameter/parameter_sampling/quadrature_sampler.py
index 5fb7ca9..a08adf8 100644
--- a/rrompy/parameter/parameter_sampling/quadrature_sampler.py
+++ b/rrompy/parameter/parameter_sampling/quadrature_sampler.py
@@ -1,87 +1,83 @@
# Copyright (C) 2018 by the RROMPy authors
#
# This file is part of RROMPy.
#
# RROMPy is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# RROMPy is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with RROMPy. If not, see .
#
import numpy as np
from .generic_sampler import GenericSampler
from rrompy.utilities.base.types import List, paramList
from rrompy.utilities.exception_manager import RROMPyException
from rrompy.utilities.numerical import lowDiscrepancy, kroneckerer
from rrompy.parameter import checkParameterList
__all__ = ['QuadratureSampler']
class QuadratureSampler(GenericSampler):
"""Generator of quadrature sample points."""
allowedKinds = ["UNIFORM", "CHEBYSHEV", "GAUSSLEGENDRE"]
def __init__(self, lims:paramList, kind : str = "UNIFORM",
- scaling : List[callable] = None,
- scalingInv : List[callable] = None):
- super().__init__(lims = lims, scaling = scaling,
- scalingInv = scalingInv)
+ scalingExp : List[float] = None):
+ super().__init__(lims = lims, scalingExp = scalingExp)
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 RROMPyException("Generator kind not recognized.")
self._kind = kind.upper()
def generatePoints(self, n:List[int], reorder : bool = True) -> paramList:
"""Array of sample points."""
if not hasattr(n, "__len__"): n = [n]
super().generatePoints(n)
nleft, nright = 1, np.prod(n)
xmat = np.empty((nright, self.npar), dtype = self.lims.dtype)
for d in range(self.npar):
nright //= n[d]
- a, b = self.lims(0, d), self.lims(1, d)
- if self.scaling is not None:
- a, b = self.scaling[d](a), self.scaling[d](b)
+ a = self.lims(0, d) ** self.scalingExp[d]
+ b = self.lims(1, d) ** self.scalingExp[d]
c, r = (a + b) / 2., (a - b) / 2.
dAbs = 2. * np.abs(r)
if self.kind == "UNIFORM":
xd = np.linspace(a, b, n[d])
elif self.kind == "CHEBYSHEV":
nodes, _ = np.polynomial.chebyshev.chebgauss(n[d])
xd = c + r * nodes
elif self.kind == "GAUSSLEGENDRE":
nodes, _ = np.polynomial.legendre.leggauss(n[d])
xd = c + r * nodes[::-1]
- if self.scalingInv is not None:
- xd = self.scalingInv[d](xd)
+ xd **= 1. / self.scalingExp[d]
xmat[:, d] = kroneckerer(xd, nleft, nright)
nleft *= n[d]
nright = np.prod(n)
if nright > 1 and reorder:
fejerOrdering = [nright - 1] + lowDiscrepancy(nright - 1)
xmat = xmat[fejerOrdering, :]
x = checkParameterList(xmat, self.npar)[0]
return x
diff --git a/rrompy/parameter/parameter_sampling/random_sampler.py b/rrompy/parameter/parameter_sampling/random_sampler.py
index 460d4e8..68625d2 100644
--- a/rrompy/parameter/parameter_sampling/random_sampler.py
+++ b/rrompy/parameter/parameter_sampling/random_sampler.py
@@ -1,76 +1,72 @@
# Copyright (C) 2018 by the RROMPy authors
#
# This file is part of RROMPy.
#
# RROMPy is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# RROMPy is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with RROMPy. If not, see .
#
import numpy as np
from .generic_sampler import GenericSampler
from rrompy.utilities.numerical import haltonGenerate, sobolGenerate
from rrompy.utilities.base.types import List, paramList
from rrompy.utilities.exception_manager import RROMPyException
from rrompy.parameter import checkParameterList
__all__ = ['RandomSampler']
class RandomSampler(GenericSampler):
"""Generator of quadrature sample points."""
allowedKinds = ["UNIFORM", "HALTON", "SOBOL"]
def __init__(self, lims:paramList, kind : str = "UNIFORM",
- scaling : List[callable] = None,
- scalingInv : List[callable] = None, seed : int = 42):
- super().__init__(lims = lims, scaling = scaling,
- scalingInv = scalingInv)
+ scalingExp : List[float] = None, seed : int = 42):
+ super().__init__(lims = lims, scalingExp = scalingExp)
self.kind = kind
self.seed = seed
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 RROMPyException("Generator kind not recognized.")
self._kind = kind.upper()
def generatePoints(self, n:int) -> paramList:
"""Array of quadrature points."""
if hasattr(n, "__len__"): n = n[0]
if self.kind == "UNIFORM":
np.random.seed(self.seed)
xmat = np.random.uniform(size = (n, self.npar))
elif self.kind == "HALTON":
xmat = haltonGenerate(self.npar, n, self.seed)
else:
xmat = sobolGenerate(self.npar, n, self.seed)
for d in range(self.npar):
- a, b = self.lims(0, d), self.lims(1, d)
- if self.scaling is not None:
- a, b = self.scaling[d](a), self.scaling[d](b)
+ a = self.lims(0, d) ** self.scalingExp[d]
+ b = self.lims(1, d) ** self.scalingExp[d]
xmat[:, d] = a + (b - a) * xmat[:, d]
- if self.scalingInv is not None:
- xmat[:, d] = self.scalingInv[d](xmat[:, d])
+ xmat[:, d] **= 1. / self.scalingExp[d]
x = checkParameterList(xmat, self.npar)[0]
return x
diff --git a/rrompy/reduction_methods/base/__init__.py b/rrompy/reduction_methods/base/__init__.py
index 98d784c..f6aa877 100644
--- a/rrompy/reduction_methods/base/__init__.py
+++ b/rrompy/reduction_methods/base/__init__.py
@@ -1,29 +1,29 @@
# Copyright (C) 2018 by the RROMPy authors
#
# This file is part of RROMPy.
#
# RROMPy is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# RROMPy is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with RROMPy. If not, see .
#
from .generic_approximant import GenericApproximant
-from .pade_utils import checkRobustTolerance
-from .rb_utils import projectAffineDecomposition
+from .rational_interpolant_utils import checkRobustTolerance
+from .reduced_basis_utils import projectAffineDecomposition
__all__ = [
'GenericApproximant',
'checkRobustTolerance',
'projectAffineDecomposition'
]
diff --git a/rrompy/reduction_methods/base/generic_approximant.py b/rrompy/reduction_methods/base/generic_approximant.py
index 53ca38f..294f1ed 100644
--- a/rrompy/reduction_methods/base/generic_approximant.py
+++ b/rrompy/reduction_methods/base/generic_approximant.py
@@ -1,868 +1,868 @@
# 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, homogeneized : bool = False) -> Np1D:
uV = getattr(self.__class__, "get" + fieldName)(self, mu, homogeneized)
val = self.HFEngine.norm(uV)
return val
setattr(self.__class__, "norm" + fieldName, objFunc)
def addPlotFieldToClass(self, fieldName):
def objFunc(self, mu:paramList, *args, homogeneized : bool = False,
**kwargs):
uV = getattr(self.__class__, "get" + fieldName)(self, mu, homogeneized)
kwargsCopy = copy(kwargs)
for j, u in enumerate(uV):
if "name" in kwargs.keys():
kwargsCopy["name"] = kwargs["name"] + str(j)
self.HFEngine.plot(u, *args, **kwargs)
setattr(self.__class__, "plot" + fieldName, objFunc)
def addOutParaviewFieldToClass(self, fieldName):
def objFunc(self, mu:paramVal, *args, homogeneized : bool = False,
**kwargs):
if not hasattr(self.HFEngine, "outParaview"):
raise RROMPyException(("High fidelity engine cannot output to "
"Paraview."))
uV = getattr(self.__class__, "get" + fieldName)(self, mu, homogeneized)
kwargsCopy = copy(kwargs)
for j, u in enumerate(uV):
if "name" in kwargs.keys():
kwargsCopy["name"] = kwargs["name"] + str(j)
self.HFEngine.outParaview(u, *args, **kwargsCopy)
setattr(self.__class__, "outParaview" + fieldName, objFunc)
def addOutParaviewTimeDomainFieldToClass(self, fieldName):
def objFunc(self, mu:paramVal, *args,
homogeneized : bool = False, **kwargs):
if not hasattr(self.HFEngine, "outParaviewTimeDomain"):
raise RROMPyException(("High fidelity engine cannot output to "
"Paraview."))
uV = getattr(self.__class__, "get" + fieldName)(self, mu, homogeneized)
omega = args.pop(0) if len(args) > 0 else np.real(mu)
kwargsCopy = copy(kwargs)
for j, u in enumerate(uV):
if "name" in kwargs.keys():
kwargsCopy["name"] = kwargs["name"] + str(j)
self.HFEngine.outParaviewTimeDomain(u, omega, *args,
**kwargsCopy)
setattr(self.__class__, "outParaviewTimeDomain" + fieldName, objFunc)
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.
homogeneized(optional): Whether to homogeneize Dirichlet BCs. Defaults
to False.
verbosity(optional): Verbosity level. Defaults to 10.
Attributes:
HFEngine: HF problem solver.
trainedModel: Trained model evaluator.
mu0: Default parameter.
homogeneized: Whether to homogeneize Dirichlet BCs.
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.
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 = {}, homogeneized : bool = False,
verbosity : int = 10, timestamp : bool = True):
self._preInit()
self._mode = RROMPy_READY
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.homogeneized = homogeneized
self.approxParameters = approxParameters
self._postInit()
### add norm{HF,RHS,Approx,Res,Err} methods
"""
Compute norm of * at arbitrary parameter.
Args:
mu: Target parameter.
homogeneized(optional): Whether to remove Dirichlet BC. Defaults to
False.
Returns:
Target norm of *.
"""
for objName in ["HF", "RHS", "Err"]:
addNormFieldToClass(self, objName)
def objFunc(self, mu:paramList, homogeneized : bool = False) -> Np1D:
# uV = getattr(self.__class__, "getRes")(self, mu, homogeneized,
# duality = False)
uV = self.getRes(mu, homogeneized, duality = False)
val = self.HFEngine.norm(uV, dual = True, duality = False)
return val
setattr(self.__class__, "normRes", objFunc)
if not hasattr(self, "normApprox"):
addNormFieldToClass(self, "Approx")
### add plot{HF,RHS,Approx,Res,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.
homogeneized(optional): Whether to remove Dirichlet BC. Defaults to
False.
figspecs(optional key args): Optional arguments for matplotlib
figure creation.
"""
for objName in ["HF", "RHS", "Approx", "Res", "Err"]:
addPlotFieldToClass(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).
homogeneized(optional): Whether to remove Dirichlet BC. Defaults to
False.
"""
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.
homogeneized(optional): Whether to remove Dirichlet BC. Defaults to
False.
"""
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 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
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.parameterDefaultSoft[what] = defaultSoft[j]
for j, what in enumerate(whatCritical):
if what not in self.parameterToBeExcluded:
self.parameterListCritical += [what]
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,
verbosity = self.verbosity,
allowRepeatedSamples = True)
@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 S(self):
"""Value of S."""
return self._S
@S.setter
def S(self, S):
if not hasattr(S, "__len__"): S = [S]
if any([s <= 0 for s in S]):
raise RROMPyException("S must be positive.")
if hasattr(self, "_S") and self._S is not None: Sold = tuple(self.S)
else: Sold = -1
self._S = S
self._approxParameters["S"] = self.S
if Sold != tuple(self.S):
self.resetSamples()
@property
def homogeneized(self):
"""Value of homogeneized."""
return self._homogeneized
@homogeneized.setter
def homogeneized(self, homogeneized):
if not hasattr(self, "_homogeneized"):
self._homogeneized = None
if homogeneized != self.homogeneized:
self._homogeneized = homogeneized
self.resetSamples()
@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):
"""
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.
"""
RROMPyAssert(self._mode, message = "Cannot plot samples.")
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):
"""
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).
"""
RROMPyAssert(self._mode, message = "Cannot output samples.")
self.samplingEngine.outParaviewSamples(name = name,
filename = filename,
times = times, what = what,
forceNewFile = forceNewFile,
folders = folders,
filePW = 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):
"""
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.
"""
RROMPyAssert(self._mode, message = "Cannot output samples.")
self.samplingEngine.outParaviewTimeDomainSamples(omegas = omegas,
timeFinal = timeFinal,
periodResolution = periodResolution,
name = name, filename = filename,
forceNewFile = forceNewFile,
folders = folders)
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,
homogeneized = self.homogeneized)
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, homogeneized : bool = False,
append : bool = False, prune : bool = True) -> sampList:
"""
Get HF solution at arbitrary parameter.
Args:
mu: Target parameter.
homogeneized(optional): Whether to remove Dirichlet BC. Defaults to
False.
Returns:
HFsolution.
"""
mu = checkParameterList(mu, self.npar)[0]
idx = self.evalHF(mu, append = append, prune = prune)
uHFs = self.uHF(idx)
if self.homogeneized and not homogeneized:
for j, m in enumerate(mu):
uHFs[j] += self.HFEngine.liftDirichletData(m)
if not self.homogeneized and homogeneized:
for j, m in enumerate(mu):
uHFs[j] -= self.HFEngine.liftDirichletData(m)
return uHFs
def getRHS(self, mu:paramList, homogeneized : bool = False,
duality : bool = True) -> sampList:
"""
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,
duality = duality)
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)
uApproxRs = self.uApproxReduced(idx)
return uApproxRs
def getApprox(self, mu:paramList, homogeneized : bool = False,
append : bool = False, prune : bool = True) -> sampList:
"""
Get approximant at arbitrary parameter.
Args:
mu: Target parameter.
homogeneized(optional): Whether to remove Dirichlet BC. Defaults to
False.
Returns:
Approximant.
"""
mu = checkParameterList(mu, self.npar)[0]
idx = self.evalApprox(mu, append = append, prune = prune)
uApproxs = self.uApprox(idx)
if self.homogeneized and not homogeneized:
for j, m in enumerate(mu):
uApproxs[j] += self.HFEngine.liftDirichletData(m)
if not self.homogeneized and homogeneized:
for j, m in enumerate(mu):
uApproxs[j] -= self.HFEngine.liftDirichletData(m)
return uApproxs
def getRes(self, mu:paramList, homogeneized : bool = False,
duality : bool = True) -> sampList:
"""
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.getApprox(mu, homogeneized), mu,
homogeneized = homogeneized,
duality = duality)
def getErr(self, mu:paramList, homogeneized : bool = False,
append : bool = False, prune : bool = True) -> sampList:
"""
Get error at arbitrary parameter.
Args:
mu: Target parameter.
homogeneized(optional): Whether to remove Dirichlet BC. Defaults to
False.
Returns:
Approximant error.
"""
return (self.getApprox(mu, homogeneized, append = append, prune =prune)
- self.getHF(mu, homogeneized, append = append, prune = prune))
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)
name = datadict.pop("name")
if name == "TrainedModelRational":
from rrompy.reduction_methods.trained_model import \
TrainedModelRational as tModel
- elif name == "TrainedModelRB":
+ elif name == "TrainedModelReducedBasis":
from rrompy.reduction_methods.trained_model import \
- TrainedModelRB as tModel
+ TrainedModelReducedBasis as tModel
else:
raise RROMPyException(("Trained model name not recognized. "
"Loading failed."))
self.mu0 = datadict.pop("mu0")
from rrompy.reduction_methods.trained_model import TrainedModelData
trainedModel = tModel()
trainedModel.verbosity = self.verbosity
trainedModel.timestamp = self.timestamp
self.scaleFactor = datadict.pop("scaleFactor")
data = TrainedModelData(name, self.mu0, datadict.pop("projMat"),
self.scaleFactor, datadict.pop("rescalingExp"))
if "mus" in datadict:
data.mus = datadict.pop("mus")
approxParameters = datadict.pop("approxParameters")
data.approxParameters = copy(approxParameters)
if "sampler" in approxParameters:
self._approxParameters["sampler"] = approxParameters.pop("sampler")
self.approxParameters = copy(approxParameters)
if "mus" in data.__dict__:
self.mus = copy(data.mus)
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/base/pade_utils.py b/rrompy/reduction_methods/base/rational_interpolant_utils.py
similarity index 100%
rename from rrompy/reduction_methods/base/pade_utils.py
rename to rrompy/reduction_methods/base/rational_interpolant_utils.py
diff --git a/rrompy/reduction_methods/base/rb_utils.py b/rrompy/reduction_methods/base/reduced_basis_utils.py
similarity index 100%
rename from rrompy/reduction_methods/base/rb_utils.py
rename to rrompy/reduction_methods/base/reduced_basis_utils.py
diff --git a/rrompy/reduction_methods/distributed/__init__.py b/rrompy/reduction_methods/distributed/__init__.py
index 1d7c1fc..b136774 100644
--- a/rrompy/reduction_methods/distributed/__init__.py
+++ b/rrompy/reduction_methods/distributed/__init__.py
@@ -1,29 +1,29 @@
# Copyright (C) 2018 by the RROMPy authors
#
# This file is part of RROMPy.
#
# RROMPy is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# RROMPy is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with RROMPy. If not, see .
#
from .generic_distributed_approximant import GenericDistributedApproximant
from .rational_interpolant import RationalInterpolant
-from .rb_distributed import RBDistributed
+from .reduced_basis import ReducedBasis
__all__ = [
'GenericDistributedApproximant',
'RationalInterpolant',
- 'RBDistributed'
+ 'ReducedBasis'
]
diff --git a/rrompy/reduction_methods/distributed/rb_distributed.py b/rrompy/reduction_methods/distributed/reduced_basis.py
similarity index 96%
rename from rrompy/reduction_methods/distributed/rb_distributed.py
rename to rrompy/reduction_methods/distributed/reduced_basis.py
index 3f142ed..d63bbf4 100644
--- a/rrompy/reduction_methods/distributed/rb_distributed.py
+++ b/rrompy/reduction_methods/distributed/reduced_basis.py
@@ -1,207 +1,209 @@
# 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_distributed_approximant import GenericDistributedApproximant
-from rrompy.reduction_methods.trained_model import TrainedModelRB as tModel
+from rrompy.reduction_methods.trained_model import \
+ TrainedModelReducedBasis as tModel
from rrompy.reduction_methods.trained_model import TrainedModelData
-from rrompy.reduction_methods.base.rb_utils import projectAffineDecomposition
+from rrompy.reduction_methods.base.reduced_basis_utils import \
+ projectAffineDecomposition
from rrompy.utilities.base.types import (Np1D, Np2D, List, ListAny, Tuple,
DictAny, HFEng, paramVal, sampList)
from rrompy.utilities.base import verbosityManager as vbMng
from rrompy.utilities.exception_manager import (RROMPyWarning, RROMPyException,
RROMPyAssert)
-__all__ = ['RBDistributed']
+__all__ = ['ReducedBasis']
-class RBDistributed(GenericDistributedApproximant):
+class ReducedBasis(GenericDistributedApproximant):
"""
ROM RB approximant computation for parametric problems.
Args:
HFEngine: HF problem solver.
mu0(optional): Default parameter. Defaults to 0.
approxParameters(optional): Dictionary containing values for main
parameters of approximant. Recognized keys are:
- 'POD': whether to compute POD of snapshots; defaults to True;
- 'S': total number of samples current approximant relies upon;
- 'sampler': sample point generator;
- 'R': rank for Galerkin projection; defaults to prod(S);
- 'PODTolerance': tolerance for snapshots POD; defaults to -1.
Defaults to empty dict.
homogeneized(optional): Whether to homogeneize Dirichlet BCs. Defaults
to False.
verbosity(optional): Verbosity level. Defaults to 10.
Attributes:
HFEngine: HF problem solver.
mu0: Default parameter.
mus: Array of snapshot parameters.
homogeneized: Whether to homogeneize Dirichlet BCs.
approxRadius: Dummy radius of approximant (i.e. distance from mu0 to
farthest sample point).
approxParameters: Dictionary containing values for main parameters of
approximant. Recognized keys are in parameterList.
parameterListSoft: Recognized keys of soft approximant parameters:
- 'POD': whether to compute POD of snapshots.
- 'R': rank for Galerkin projection;
- 'PODTolerance': tolerance for snapshots POD.
parameterListCritical: Recognized keys of critical approximant
parameters:
- 'S': total number of samples current approximant relies upon;
- 'sampler': sample point generator;
POD: Whether to compute POD of snapshots.
S: Number of solution snapshots over which current approximant is
based upon.
sampler: Sample point generator.
R: Rank for Galerkin projection.
muBounds: list of bounds for parameter values.
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.
As: List of sparse matrices (in CSC format) representing coefficients
of linear system matrix.
bs: List of numpy vectors representing coefficients of linear system
RHS.
ARBs: List of sparse matrices (in CSC format) representing coefficients
of compressed linear system matrix.
bRBs: List of numpy vectors representing coefficients of compressed
linear system RHS.
"""
def __init__(self, HFEngine:HFEng, mu0 : paramVal = None,
approxParameters : DictAny = {}, homogeneized : bool = False,
verbosity : int = 10, timestamp : bool = True):
self._preInit()
self._addParametersToList(["R", "PODTolerance"], [1, -1])
super().__init__(HFEngine = HFEngine, mu0 = mu0,
approxParameters = approxParameters,
homogeneized = homogeneized,
verbosity = verbosity, timestamp = timestamp)
vbMng(self, "INIT", "Computing affine blocks of system.", 10)
vbMng(self, "DEL", "Done computing affine blocks.", 10)
self._postInit()
@property
def S(self):
"""Value of S."""
return self._S
@S.setter
def S(self, S):
GenericDistributedApproximant.S.fset(self, S)
if not hasattr(self, "_R"): self._R = np.prod(self.S)
@property
def R(self):
"""Value of R. Its assignment may change S."""
return self._R
@R.setter
def R(self, R):
if R < 0: raise RROMPyException("R must be non-negative.")
self._R = R
self._approxParameters["R"] = self.R
if hasattr(self, "_S") and np.prod(self.S) < self.R:
RROMPyWarning("Prescribed S is too small. Decreasing R.")
self.R = np.prod(self.S)
@property
def PODTolerance(self):
"""Value of PODTolerance."""
return self._PODTolerance
@PODTolerance.setter
def PODTolerance(self, PODTolerance):
self._PODTolerance = PODTolerance
self._approxParameters["PODTolerance"] = self.PODTolerance
def setupApprox(self):
"""Compute RB projection matrix."""
if self.checkComputedApprox():
return
RROMPyAssert(self._mode, message = "Cannot setup approximant.")
vbMng(self, "INIT", "Setting up {}.". format(self.name()), 5)
self.computeScaleFactor()
self.computeSnapshots()
vbMng(self, "INIT", "Computing projection matrix.", 7)
if self.POD:
U, s, _ = np.linalg.svd(self.samplingEngine.RPOD)
s = s ** 2.
else:
Gramian = self.HFEngine.innerProduct(self.samplingEngine.samples,
self.samplingEngine.samples)
U, s, _ = np.linalg.svd(Gramian)
nsamples = self.samplingEngine.nsamples
snorm = np.cumsum(s[::-1]) / np.sum(s)
nPODTrunc = min(nsamples - np.argmax(snorm > self.PODTolerance),
self.R)
pMat = self.samplingEngine.samples.dot(U[:, : nPODTrunc])
vbMng(self, "MAIN",
("Assembling {}x{} projection matrix from {} "
"samples.").format(*(pMat.shape), nsamples), 5)
if self.trainedModel is None:
self.trainedModel = tModel()
self.trainedModel.verbosity = self.verbosity
self.trainedModel.timestamp = self.timestamp
data = TrainedModelData(self.trainedModel.name(), self.mu0,
pMat, self.scaleFactor,
self.HFEngine.rescalingExp)
data.mus = copy(self.mus)
self.trainedModel.data = data
else:
self.trainedModel = self.trainedModel
self.trainedModel.data.projMat = copy(pMat)
vbMng(self, "INIT", "Computing affine blocks of system.", 10)
self.As = self.HFEngine.affineLinearSystemA(self.mu0, self.scaleFactor)
self.bs = self.HFEngine.affineLinearSystemb(self.mu0, self.scaleFactor,
self.homogeneized)
vbMng(self, "DEL", "Done computing affine blocks.", 10)
ARBs, bRBs = self.assembleReducedSystem(pMat)
self.trainedModel.data.ARBs = ARBs
self.trainedModel.data.bRBs = bRBs
vbMng(self, "DEL", "Done computing projection matrix.", 7)
self.trainedModel.data.approxParameters = copy(self.approxParameters)
vbMng(self, "DEL", "Done setting up approximant.", 5)
def assembleReducedSystem(self, pMat : sampList = None,
pMatOld : sampList = None)\
-> Tuple[List[Np2D], List[Np1D]]:
"""Build affine blocks of RB linear system through projections."""
if pMat is None:
self.setupApprox()
ARBs = self.trainedModel.data.ARBs
bRBs = self.trainedModel.data.bRBs
else:
vbMng(self, "INIT", "Projecting affine terms of HF model.", 10)
ARBsOld = None if pMatOld is None else self.trainedModel.data.ARBs
bRBsOld = None if pMatOld is None else self.trainedModel.data.bRBs
ARBs, bRBs = projectAffineDecomposition(self.As, self.bs, pMat,
ARBsOld, bRBsOld, pMatOld)
vbMng(self, "DEL", "Done projecting affine terms.", 10)
return ARBs, bRBs
diff --git a/rrompy/reduction_methods/distributed_greedy/__init__.py b/rrompy/reduction_methods/distributed_greedy/__init__.py
index ba0008b..ccdf783 100644
--- a/rrompy/reduction_methods/distributed_greedy/__init__.py
+++ b/rrompy/reduction_methods/distributed_greedy/__init__.py
@@ -1,30 +1,30 @@
# Copyright (C) 2018 by the RROMPy authors
#
# This file is part of RROMPy.
#
# RROMPy is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# RROMPy is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with RROMPy. If not, see .
#
from .generic_distributed_greedy_approximant import \
GenericDistributedGreedyApproximant
from .rational_interpolant_greedy import RationalInterpolantGreedy
-from .rb_distributed_greedy import RBDistributedGreedy
+from .reduced_basis_greedy import ReducedBasisGreedy
__all__ = [
'GenericDistributedGreedyApproximant',
'RationalInterpolantGreedy',
- 'RBDistributedGreedy'
+ 'ReducedBasisGreedy'
]
diff --git a/rrompy/reduction_methods/distributed_greedy/rb_distributed_greedy.py b/rrompy/reduction_methods/distributed_greedy/reduced_basis_greedy.py
similarity index 97%
rename from rrompy/reduction_methods/distributed_greedy/rb_distributed_greedy.py
rename to rrompy/reduction_methods/distributed_greedy/reduced_basis_greedy.py
index 1bb837c..9fc0247 100644
--- a/rrompy/reduction_methods/distributed_greedy/rb_distributed_greedy.py
+++ b/rrompy/reduction_methods/distributed_greedy/reduced_basis_greedy.py
@@ -1,238 +1,239 @@
# Copyright (C) 2018 by the RROMPy authors
#
# This file is part of RROMPy.
#
# RROMPy is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# RROMPy is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with RROMPy. If not, see .
#
import numpy as np
from copy import deepcopy as copy
from .generic_distributed_greedy_approximant import \
GenericDistributedGreedyApproximant
-from rrompy.reduction_methods.distributed import RBDistributed
-from rrompy.reduction_methods.trained_model import TrainedModelRB as tModel
+from rrompy.reduction_methods.distributed import ReducedBasis
+from rrompy.reduction_methods.trained_model import \
+ TrainedModelReducedBasis as tModel
from rrompy.reduction_methods.trained_model import TrainedModelData
from rrompy.utilities.base.types import Np1D, DictAny, HFEng, paramVal
from rrompy.utilities.base import verbosityManager as vbMng
from rrompy.utilities.poly_fitting.polynomial import (
hashIdxToDerivative as hashI)
from rrompy.utilities.exception_manager import RROMPyWarning, RROMPyAssert
from rrompy.parameter import checkParameterList
-__all__ = ['RBDistributedGreedy']
+__all__ = ['ReducedBasisGreedy']
-class RBDistributedGreedy(GenericDistributedGreedyApproximant, RBDistributed):
+class ReducedBasisGreedy(GenericDistributedGreedyApproximant, ReducedBasis):
"""
ROM greedy RB approximant computation for parametric problems.
Args:
HFEngine: HF problem solver.
mu0(optional): Default parameter. Defaults to 0.
approxParameters(optional): Dictionary containing values for main
parameters of approximant. Recognized keys are:
- 'POD': whether to compute POD of snapshots; defaults to True;
- 'S': number of starting training points;
- 'sampler': sample point generator;
- 'greedyTol': uniform error tolerance for greedy algorithm;
defaults to 1e-2;
- '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 Chebyshev sampler within muBounds.
Defaults to empty dict.
homogeneized(optional): Whether to homogeneize Dirichlet BCs. Defaults
to False.
verbosity(optional): Verbosity level. Defaults to 10.
Attributes:
HFEngine: HF problem solver.
mu0: Default parameter.
mus: Array of snapshot parameters.
homogeneized: Whether to homogeneize Dirichlet BCs.
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;
- '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.
parameterListCritical: Recognized keys of critical approximant
parameters:
- 'S': total number of samples current approximant relies upon;
- 'sampler': sample point generator.
POD: whether to compute POD of snapshots.
S: number of test points.
sampler: Sample point generator.
greedyTol: uniform error 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.
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.
As: List of sparse matrices (in CSC format) representing coefficients
of linear system matrix.
bs: List of numpy vectors representing coefficients of linear system
RHS.
ARBs: List of sparse matrices (in CSC format) representing coefficients
of compressed linear system matrix.
bRBs: List of numpy vectors representing coefficients of compressed
linear system RHS.
"""
def __init__(self, HFEngine:HFEng, mu0 : paramVal = None,
approxParameters : DictAny = {}, homogeneized : bool = False,
verbosity : int = 10, timestamp : bool = True):
self._preInit()
super().__init__(HFEngine = HFEngine, mu0 = mu0,
approxParameters = approxParameters,
homogeneized = homogeneized,
verbosity = verbosity, timestamp = timestamp)
self._postInit()
@property
def R(self):
"""Value of R."""
self._R = np.prod(self._S)
return self._R
@R.setter
def R(self, R):
RROMPyWarning(("R is used just to simplify inheritance, and its value "
"cannot be changed from that of prod(S)."))
@property
def PODTolerance(self):
"""Value of PODTolerance."""
self._PODTolerance = -1
return self._PODTolerance
@PODTolerance.setter
def PODTolerance(self, PODTolerance):
RROMPyWarning(("PODTolerance is used just to simplify inheritance, "
"and its value cannot be changed from -1."))
def errorEstimator(self, mus:Np1D) -> Np1D:
"""
Standard residual-based error estimator. Unreliable for unstable
problems (inf-sup constant is missing).
"""
self.setupApprox()
self.assembleReducedResidualBlocks(full = True)
nmus = len(mus)
nAs = self.trainedModel.data.resAA.shape[1]
nbs = self.trainedModel.data.resbb.shape[0]
radiusA = np.empty((len(self.mus), nAs, nmus), dtype = np.complex)
radiusb = np.empty((nbs, nmus), dtype = np.complex)
verb = self.trainedModel.verbosity
self.trainedModel.verbosity = 0
if verb >= 5:
mustr = mus
if nmus > 2:
mustr = "[{} ..({}).. {}]".format(mus[0], nmus - 2, mus[-1])
vbMng(self, "INIT",
"Computing RB solution at mu = {}.".format(mustr), 0)
parmus = checkParameterList(mus, self.npar)[0]
uApproxRs = self.getApproxReduced(parmus)
mu0Eff = self.mu0(0, 0) ** self.HFEngine.rescalingExp[0]
for j, muPL in enumerate(parmus):
uApproxR = uApproxRs[j]
for i in range(max(nAs, nbs)):
thetai = ((muPL[0] ** self.HFEngine.rescalingExp[0]
- mu0Eff) / self.scaleFactor[0]) ** i
if i < nAs:
radiusA[:, i, j] = thetai * uApproxR
if i < nbs:
radiusb[i, j] = thetai
vbMng(self, "DEL", "Done computing RB solution.", 5)
self.trainedModel.verbosity = verb
# 'ij,jk,ik->k', resbb, radiusb, radiusb.conj()
ff = np.sum(self.trainedModel.data.resbb.dot(radiusb) * radiusb.conj(),
axis = 0)
# 'ijk,jkl,il->l', resAb, radiusA, radiusb.conj()
Lf = np.sum(np.tensordot(self.trainedModel.data.resAb, radiusA, 2)
* radiusb.conj(), axis = 0)
# 'ijkl,klt,ijt->t', resAA, radiusA, radiusA.conj()
LL = np.sum(np.tensordot(self.trainedModel.data.resAA, radiusA, 2)
* radiusA.conj(), axis = (0, 1))
return np.abs((LL - 2. * np.real(Lf) + ff) / ff) ** .5
def setupApprox(self, plotEst : bool = False):
"""Compute RB projection matrix."""
if self.checkComputedApprox():
return
RROMPyAssert(self._mode, message = "Cannot setup approximant.")
vbMng(self, "INIT", "Setting up {}.". format(self.name()), 5)
self.computeScaleFactor()
self.greedy(plotEst)
vbMng(self, "INIT", "Computing projection matrix.", 7)
pMat = self.samplingEngine.samples
vbMng(self, "INIT", "Computing affine blocks of system.", 10)
self.As = self.HFEngine.affineLinearSystemA(self.mu0, self.scaleFactor)
self.bs = self.HFEngine.affineLinearSystemb(self.mu0, self.scaleFactor,
self.homogeneized)
vbMng(self, "DEL", "Done computing affine blocks.", 10)
if self.trainedModel is None:
self.trainedModel = tModel()
self.trainedModel.verbosity = self.verbosity
self.trainedModel.timestamp = self.timestamp
data = TrainedModelData(self.trainedModel.name(), self.mu0,
pMat, self.scaleFactor,
self.HFEngine.rescalingExp)
ARBs, bRBs = self.assembleReducedSystem(pMat)
self.trainedModel.data = data
else:
self.trainedModel = self.trainedModel
pMatOld = self.trainedModel.data.projMat
Sold = pMatOld.shape[1]
idxNew = list(range(Sold, pMat.shape[1]))
ARBs, bRBs = self.assembleReducedSystem(pMat(idxNew), pMatOld)
self.trainedModel.data.projMat = copy(pMat)
self.trainedModel.data.mus = copy(self.mus)
self.trainedModel.data.ARBs = ARBs
self.trainedModel.data.bRBs = bRBs
vbMng(self, "DEL", "Done computing projection matrix.", 7)
self.trainedModel.data.approxParameters = copy(self.approxParameters)
vbMng(self, "DEL", "Done setting up approximant.", 5)
def assembleReducedResidualBlocks(self, full : bool = False):
"""Build affine blocks of RB linear system through projections."""
self.assembleReducedResidualBlocksbb(self.bs)
if full:
pMat = self.trainedModel.data.projMat
self.assembleReducedResidualBlocksAb(self.As, self.bs, pMat)
self.assembleReducedResidualBlocksAA(self.As, pMat)
diff --git a/rrompy/reduction_methods/pivoting/__init__.py b/rrompy/reduction_methods/pivoting/__init__.py
index 9830125..6237733 100644
--- a/rrompy/reduction_methods/pivoting/__init__.py
+++ b/rrompy/reduction_methods/pivoting/__init__.py
@@ -1,29 +1,29 @@
# Copyright (C) 2018 by the RROMPy authors
#
# This file is part of RROMPy.
#
# RROMPy is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# RROMPy is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with RROMPy. If not, see .
#
from .generic_pivoted_approximant import GenericPivotedApproximant
from .rational_interpolant_pivoted import RationalInterpolantPivoted
-from .rb_distributed_pivoted import RBDistributedPivoted
+from .reduced_basis_pivoted import ReducedBasisPivoted
__all__ = [
'GenericPivotedApproximant',
'RationalInterpolantPivoted',
- 'RBDistributedPivoted'
+ 'ReducedBasisPivoted'
]
diff --git a/rrompy/reduction_methods/pivoting/rb_distributed_pivoted.py b/rrompy/reduction_methods/pivoting/reduced_basis_pivoted.py
similarity index 97%
rename from rrompy/reduction_methods/pivoting/rb_distributed_pivoted.py
rename to rrompy/reduction_methods/pivoting/reduced_basis_pivoted.py
index 4fabd41..1200f83 100644
--- a/rrompy/reduction_methods/pivoting/rb_distributed_pivoted.py
+++ b/rrompy/reduction_methods/pivoting/reduced_basis_pivoted.py
@@ -1,284 +1,285 @@
# 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_pivoted_approximant import GenericPivotedApproximant
from rrompy.hfengines.base import MarginalProxyEngine
from rrompy.reduction_methods.trained_model import (TrainedModelPivotedData,
- TrainedModelPivotedRB as tModel)
-from rrompy.reduction_methods.base.rb_utils import projectAffineDecomposition
+ TrainedModelPivotedReducedBasis as tModel)
+from rrompy.reduction_methods.base.reduced_basis_utils import \
+ projectAffineDecomposition
from rrompy.utilities.base.types import (Np1D, Np2D, List, ListAny, Tuple,
DictAny, HFEng, paramVal, sampList)
from rrompy.utilities.base import verbosityManager as vbMng, freepar as fp
from rrompy.utilities.numerical import customPInv
from rrompy.utilities.exception_manager import (RROMPyWarning, RROMPyException,
RROMPyAssert)
-__all__ = ['RBDistributedPivoted']
+__all__ = ['ReducedBasisPivoted']
-class RBDistributedPivoted(GenericPivotedApproximant):
+class ReducedBasisPivoted(GenericPivotedApproximant):
"""
ROM pivoted RB approximant computation for parametric problems.
Args:
HFEngine: HF problem solver.
mu0(optional): Default parameter. Defaults to 0.
directionPivot(optional): Pivot components. 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 pivot samples current approximant relies
upon;
- 'samplerPivot': pivot sample point generator;
- 'SMarginal': total number of marginal samples current approximant
relies upon;
- 'samplerMarginal': marginal sample point generator;
- 'R': rank for pivot Galerkin projection; defaults to prod(S);
- 'PODTolerance': tolerance for pivot snapshots POD; defaults to
-1;
- 'polybasisMarginal': type of polynomial basis for marginal
interpolation; allowed values include 'MONOMIAL', 'CHEBYSHEV'
and 'LEGENDRE'; defaults to 'MONOMIAL';
- 'MMarginal': degree of marginal interpolant; defaults to 0;
- 'radialBasisMarginal': radial basis family for marginal
interpolant; defaults to 0, i.e. no radial basis;
- 'radialBasisWeightsMarginal': radial basis weights for marginal
interpolant; defaults to 0, i.e. identity;
- 'interpRcondMarginal': tolerance for marginal interpolation;
defaults to None.
Defaults to empty dict.
homogeneized(optional): Whether to homogeneize Dirichlet BCs. Defaults
to False.
verbosity(optional): Verbosity level. Defaults to 10.
Attributes:
HFEngine: HF problem solver.
mu0: Default parameter.
directionPivot: Pivot components.
mus: Array of snapshot parameters.
homogeneized: Whether to homogeneize Dirichlet BCs.
approxRadius: Dummy radius of approximant (i.e. distance from mu0 to
farthest sample point).
approxParameters: Dictionary containing values for main parameters of
approximant. Recognized keys are in parameterList.
parameterListSoft: Recognized keys of soft approximant parameters:
- 'POD': whether to compute POD of snapshots;
- 'R': rank for Galerkin projection;
- 'PODTolerance': tolerance for snapshots POD;
- 'polybasisMarginal': type of polynomial basis for marginal
interpolation;
- 'MMarginal': degree of marginal interpolant;
- 'radialBasisMarginal': radial basis family for marginal
interpolant;
- 'radialBasisWeightsMarginal': radial basis weights for marginal
interpolant;
- 'interpRcondMarginal': tolerance for marginal interpolation.
parameterListCritical: Recognized keys of critical approximant
parameters:
- 'S': total number of pivot samples current approximant relies
upon;
- 'samplerPivot': pivot sample point generator;
- 'SMarginal': total number of marginal samples current approximant
relies upon;
- 'samplerMarginal': marginal sample point generator.
POD: Whether to compute POD of snapshots.
S: Total number of pivot samples current approximant relies upon.
samplerPivot: Pivot sample point generator.
SMarginal: Total number of marginal samples current approximant relies
upon.
samplerMarginal: Marginal sample point generator.
R: Rank for Galerkin projection.
PODTolerance: Tolerance for pivot snapshots POD.
polybasisMarginal: Type of polynomial basis for marginal interpolation.
MMarginal: Degree of marginal interpolant.
radialBasisMarginal: Radial basis family for marginal interpolant.
radialBasisWeightsMarginal: Radial basis weights for marginal
interpolant.
interpRcondMarginal: Tolerance for marginal interpolation.
muBoundsPivot: list of bounds for pivot parameter values.
muBoundsMarginal: list of bounds for marginal parameter values.
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.
As: List of sparse matrices (in CSC format) representing coefficients
of linear system matrix.
bs: List of numpy vectors representing coefficients of linear system
RHS.
ARBs: List of sparse matrices (in CSC format) representing coefficients
of compressed linear system matrix.
bRBs: List of numpy vectors representing coefficients of compressed
linear system RHS.
"""
def __init__(self, HFEngine:HFEng, mu0 : paramVal = None,
directionPivot : ListAny = [0],
approxParameters : DictAny = {}, homogeneized : bool = False,
verbosity : int = 10, timestamp : bool = True):
self._preInit()
self._addParametersToList(["R", "PODTolerance"], [1, -1])
super().__init__(HFEngine = HFEngine, mu0 = mu0,
directionPivot = directionPivot,
approxParameters = approxParameters,
homogeneized = homogeneized,
verbosity = verbosity, timestamp = timestamp)
self._postInit()
@property
def S(self):
"""Value of S."""
return self._S
@S.setter
def S(self, S):
GenericPivotedApproximant.S.fset(self, S)
if not hasattr(self, "_R"):
self._R = np.prod(self.S) * np.prod(self.SMarginal)
@property
def R(self):
"""Value of R. Its assignment may change S."""
return self._R
@R.setter
def R(self, R):
if R < 0: raise RROMPyException("R must be non-negative.")
self._R = R
self._approxParameters["R"] = self.R
if (hasattr(self, "_S") and hasattr(self, "_SMarginal")
and np.prod(self.S) * np.prod(self.SMarginal) < self.R):
RROMPyWarning(("Prescribed S and SMarginal are too small. "
"Decreasing R."))
self.R = np.prod(self.S) * np.prod(self.SMarginal)
@property
def PODTolerance(self):
"""Value of PODTolerance."""
return self._PODTolerance
@PODTolerance.setter
def PODTolerance(self, PODTolerance):
self._PODTolerance = PODTolerance
self._approxParameters["PODTolerance"] = self.PODTolerance
def _setupProjectionMatrix(self):
"""Compute projection matrix."""
RROMPyAssert(self._mode, message = "Cannot setup numerator.")
vbMng(self, "INIT", "Starting computation of projection matrix.", 7)
if self.POD:
U, s, _ = np.linalg.svd(self.samplingEngine.RPODCoalesced)
s = s ** 2.
else:
Gramian = self.HFEngine.innerProduct(
self.samplingEngine.samplesCoalesced,
self.samplingEngine.samplesCoalesced)
U, s, _ = np.linalg.svd(Gramian)
nsamples = self.samplingEngine.samplesCoalesced.shape[1]
snorm = np.cumsum(s[::-1]) / np.sum(s)
nPODTrunc = min(nsamples - np.argmax(snorm > self.PODTolerance),
self.R)
pMat = self.samplingEngine.samplesCoalesced.dot(U[:, : nPODTrunc])
vbMng(self, "MAIN",
"Assembling {}x{} projection matrix from {} samples.".format(
*(pMat.shape), nsamples), 5)
vbMng(self, "DEL", "Done computing projection matrix.", 7)
return pMat
def _setupAffineBlocks(self):
"""Compute list of marginalized affine blocks of system."""
hasAs = hasattr(self, "AsList") and self.AsList is not None
hasbs = hasattr(self, "bsList") and self.bsList is not None
if hasAs and hasbs: return
vbMng(self, "INIT", "Computing affine blocks of system.", 10)
mu0Eff = self.mu0.data[0, self.directionPivot]
if not hasAs: self.AsList = [None] * len(self.musMarginal)
if not hasbs: self.bsList = [None] * len(self.musMarginal)
for k, muM in enumerate(self.musMarginal):
muEff = [fp] * self.npar
for jj, kk in enumerate(self.directionMarginal):
muEff[kk] = muM(0, jj)
MEnginek = MarginalProxyEngine(self.HFEngine, muEff)
if not hasAs:
self.AsList[k] = MEnginek.affineLinearSystemA(mu0Eff,
self.scaleFactorPivot)
if not hasbs:
self.bsList[k] = MEnginek.affineLinearSystemb(mu0Eff,
self.scaleFactorPivot,
self.homogeneized)
vbMng(self, "DEL", "Done computing affine blocks.", 10)
def setupApprox(self):
"""Compute RB projection matrix."""
if self.checkComputedApprox():
return
RROMPyAssert(self._mode, message = "Cannot setup approximant.")
vbMng(self, "INIT", "Setting up {}.".format(self.name()), 5)
self.computeScaleFactor()
self.computeSnapshots()
self._setupAffineBlocks()
pMat = self._setupProjectionMatrix()
ARBsList, bRBsList, pList = self.assembleReducedSystemMarginal(pMat)
if self.trainedModel is None:
self.trainedModel = tModel()
self.trainedModel.verbosity = self.verbosity
self.trainedModel.timestamp = self.timestamp
data = TrainedModelPivotedData(self.trainedModel.name(), self.mu0,
pList, self.scaleFactor,
self.HFEngine.rescalingExp,
self.directionPivot)
data.musPivot = copy(self.musPivot)
data.musMarginal = copy(self.musMarginal)
self.trainedModel.data = data
else:
self.trainedModel = self.trainedModel
self.trainedModel.data.projMat = copy(pList)
self.trainedModel.data.marginalInterp = self._setupMarginalInterp()
self.trainedModel.data.ARBsList = ARBsList
self.trainedModel.data.bRBsList = bRBsList
self.trainedModel.data.approxParameters = copy(self.approxParameters)
vbMng(self, "DEL", "Done setting up approximant.", 5)
def assembleReducedSystemMarginal(self, pMat : sampList = None)\
-> Tuple[List[List[Np2D]],
List[List[Np1D]], List[Np2D]]:
"""Build affine blocks of RB linear system through projections."""
if pMat is None:
self.setupApprox()
ARBsList = self.trainedModel.data.ARBsList
bRBsList = self.trainedModel.data.bRBsList
else:
vbMng(self, "INIT", "Projecting affine terms of HF model.", 10)
ARBsList = [None] * len(self.musMarginal)
bRBsList = [None] * len(self.musMarginal)
projList = [None] * len(self.musMarginal)
for k, (As, bs, sample) in enumerate(zip(self.AsList, self.bsList,
self.samplingEngine.samples)):
compLocal = self.HFEngine.innerProduct(sample, pMat)
projList[k] = sample.dot(customPInv(compLocal))
ARBsList[k], bRBsList[k] = projectAffineDecomposition(As, bs,
projList[k])
vbMng(self, "DEL", "Done projecting affine terms.", 10)
return ARBsList, bRBsList, projList
diff --git a/rrompy/reduction_methods/pole_matching/__init__.py b/rrompy/reduction_methods/pole_matching/__init__.py
index ae1e7db..0370cf1 100644
--- a/rrompy/reduction_methods/pole_matching/__init__.py
+++ b/rrompy/reduction_methods/pole_matching/__init__.py
@@ -1,29 +1,29 @@
# Copyright (C) 2018 by the RROMPy authors
#
# This file is part of RROMPy.
#
# RROMPy is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# RROMPy is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with RROMPy. If not, see .
#
from .generic_pole_matching_approximant import GenericPoleMatchingApproximant
from .rational_interpolant_pole_matching import RationalInterpolantPoleMatching
-from .rb_distributed_pole_matching import RBDistributedPoleMatching
+from .reduced_basis_pole_matching import ReducedBasisPoleMatching
__all__ = [
'GenericPoleMatchingApproximant',
'RationalInterpolantPoleMatching',
- 'RBDistributedPoleMatching'
+ 'ReducedBasisPoleMatching'
]
diff --git a/rrompy/reduction_methods/pole_matching/rational_interpolant_pole_matching.py b/rrompy/reduction_methods/pole_matching/rational_interpolant_pole_matching.py
index bc8eeec..736d68b 100644
--- a/rrompy/reduction_methods/pole_matching/rational_interpolant_pole_matching.py
+++ b/rrompy/reduction_methods/pole_matching/rational_interpolant_pole_matching.py
@@ -1,222 +1,222 @@
# 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.reduction_methods.pivoting.rational_interpolant_pivoted import \
RationalInterpolantPivoted
from .generic_pole_matching_approximant import GenericPoleMatchingApproximant
from rrompy.utilities.poly_fitting.polynomial import (
PolynomialInterpolator as PI)
from rrompy.reduction_methods.trained_model import (TrainedModelPivotedData,
- TrainedModelPivotedRationalPoleMatching as tModel)
+ TrainedModelPoleMatchingRational as tModel)
from rrompy.utilities.base import verbosityManager as vbMng
from rrompy.utilities.exception_manager import RROMPyAssert
__all__ = ['RationalInterpolantPoleMatching']
class RationalInterpolantPoleMatching(GenericPoleMatchingApproximant,
RationalInterpolantPivoted):
"""
ROM pivoted rational interpolant computation for parametric problems with
pole matching.
Args:
HFEngine: HF problem solver.
mu0(optional): Default parameter. Defaults to 0.
directionPivot(optional): Pivot components. 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;
- 'matchingWeight': weight for pole matching optimization; defaults
to 1;
- 'cutOffTolerance': tolerance for ignoring parasitic poles;
defaults to np.inf;
- 'cutOffType': rule for tolerance computation for parasitic poles;
defaults to 'MAGNITUDE';
- 'S': total number of pivot samples current approximant relies
upon;
- 'samplerPivot': pivot sample point generator;
- 'SMarginal': total number of marginal samples current approximant
relies upon;
- 'samplerMarginal': marginal sample point generator;
- 'polybasisPivot': type of polynomial basis for pivot
interpolation; allowed values include 'MONOMIAL', 'CHEBYSHEV'
and 'LEGENDRE'; defaults to 'MONOMIAL';
- 'polybasisMarginal': type of polynomial basis for marginal
interpolation; allowed values include 'MONOMIAL', 'CHEBYSHEV'
and 'LEGENDRE'; defaults to 'MONOMIAL';
- 'E': number of derivatives used to compute interpolant; defaults
to 0;
- 'M': degree of rational interpolant numerator; defaults to 0;
- 'N': degree of rational interpolant denominator; defaults to 0;
- 'MMarginal': degree of marginal interpolant; defaults to 0;
- 'radialBasisPivot': radial basis family for pivot numerator;
defaults to 0, i.e. no radial basis;
- 'radialBasisMarginal': radial basis family for marginal
interpolant; defaults to 0, i.e. no radial basis;
- 'radialBasisWeightsPivot': radial basis weights for pivot
numerator; defaults to 0, i.e. identity;
- 'radialBasisWeightsMarginal': radial basis weights for marginal
interpolant; defaults to 0, i.e. identity;
- 'interpRcondPivot': tolerance for pivot interpolation; defaults
to None;
- 'interpRcondMarginal': tolerance for marginal interpolation;
defaults to None;
- 'robustTol': tolerance for robust rational denominator
management; defaults to 0.
Defaults to empty dict.
homogeneized(optional): Whether to homogeneize Dirichlet BCs. Defaults
to False.
verbosity(optional): Verbosity level. Defaults to 10.
Attributes:
HFEngine: HF problem solver.
mu0: Default parameter.
directionPivot: Pivot components.
mus: Array of snapshot parameters.
musPivot: Array of pivot snapshot parameters.
musMarginal: Array of marginal snapshot parameters.
homogeneized: Whether to homogeneize Dirichlet BCs.
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;
- 'matchingWeight': weight for pole matching optimization;
- 'cutOffTolerance': tolerance for ignoring parasitic poles;
- 'cutOffType': rule for tolerance computation for parasitic poles;
- 'polybasisPivot': type of polynomial basis for pivot
interpolation;
- 'polybasisMarginal': type of polynomial basis for marginal
interpolation;
- 'E': number of derivatives used to compute interpolant;
- 'M': degree of rational interpolant numerator;
- 'N': degree of rational interpolant denominator;
- 'MMarginal': degree of marginal interpolant;
- 'radialBasisPivot': radial basis family for pivot numerator;
- 'radialBasisMarginal': radial basis family for marginal
interpolant;
- 'radialBasisWeightsPivot': radial basis weights for pivot
numerator;
- 'radialBasisWeightsMarginal': radial basis weights for marginal
interpolant;
- 'interpRcondPivot': tolerance for pivot interpolation;
- 'interpRcondMarginal': tolerance for marginal interpolation;
- 'robustTol': tolerance for robust rational denominator
management.
parameterListCritical: Recognized keys of critical approximant
parameters:
- 'S': total number of pivot samples current approximant relies
upon;
- 'samplerPivot': pivot sample point generator;
- 'SMarginal': total number of marginal samples current approximant
relies upon;
- 'samplerMarginal': marginal sample point generator.
POD: Whether to compute POD of snapshots.
matchingWeight: Weight for pole matching optimization.
cutOffTolerance: Tolerance for ignoring parasitic poles.
cutOffType: Rule for tolerance computation for parasitic poles.
S: Total number of pivot samples current approximant relies upon.
sampler: Pivot sample point generator.
polybasisPivot: Type of polynomial basis for pivot interpolation.
polybasisMarginal: Type of polynomial basis for marginal interpolation.
M: Numerator degree of approximant.
N: Denominator degree of approximant.
MMarginal: Degree of marginal interpolant.
radialBasisPivot: Radial basis family for pivot numerator.
radialBasisMarginal: Radial basis family for marginal interpolant.
radialBasisWeightsPivot: Radial basis weights for pivot numerator.
radialBasisWeightsMarginal: Radial basis weights for marginal
interpolant.
interpRcondPivot: Tolerance for pivot interpolation.
interpRcondMarginal: Tolerance for marginal interpolation.
robustTol: Tolerance for robust rational denominator management.
E: Complete derivative depth level of S.
muBoundsPivot: list of bounds for pivot parameter values.
muBoundsMarginal: list of bounds for marginal parameter values.
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.
Q: Numpy 1D vector containing complex coefficients of approximant
denominator.
P: Numpy 2D vector whose columns are FE dofs of coefficients of
approximant numerator.
"""
def setupApprox(self):
"""
Compute rational interpolant.
SVD-based robust eigenvalue management.
"""
if self.checkComputedApprox():
return
RROMPyAssert(self._mode, message = "Cannot setup approximant.")
vbMng(self, "INIT", "Setting up {}.". format(self.name()), 5)
self.computeScaleFactor()
self.computeSnapshots()
if self.trainedModel is None:
self.trainedModel = tModel()
self.trainedModel.verbosity = self.verbosity
self.trainedModel.timestamp = self.timestamp
data = TrainedModelPivotedData(self.trainedModel.name(), self.mu0,
self.samplingEngine.samplesCoalesced,
self.scaleFactor,
self.HFEngine.rescalingExp,
self.directionPivot)
data.musPivot = copy(self.musPivot)
data.musMarginal = copy(self.musMarginal)
self.trainedModel.data = data
else:
self.trainedModel = self.trainedModel
self.trainedModel.data.projMat = copy(
self.samplingEngine.samplesCoalesced)
self.trainedModel.data.marginalInterp = self._setupMarginalInterp()
if self.N > 0:
Qs = self._setupDenominator()
else:
Q = PI()
Q.npar = self.musPivot.shape[1]
Q.coeffs = np.ones(tuple([1] * Q.npar), dtype = np.complex)
Q.polybasis = self.polybasisPivot
Qs = [Q for _ in range(len(self.musMarginal))]
self.trainedModel.data._temporary = True
self.trainedModel.data.Qs = Qs
self.trainedModel.data.Ps = self._setupNumerator()
vbMng(self, "INIT", "Matching poles.", 10)
self.trainedModel.initializeFromRational(self.HFEngine,
self.matchingWeight, self.POD)
vbMng(self, "DEL", "Done matching poles.", 10)
if not np.isinf(self.cutOffTolerance):
vbMng(self, "INIT", "Recompressing by cut-off.", 10)
msg = self.trainedModel.recompressByCutOff([-1., 1.],
self.cutOffTolerance,
self.cutOffType)
vbMng(self, "DEL", "Done recompressing." + msg, 10)
self.trainedModel.data.approxParameters = copy(self.approxParameters)
vbMng(self, "DEL", "Done setting up approximant.", 5)
diff --git a/rrompy/reduction_methods/pole_matching/rb_distributed_pole_matching.py b/rrompy/reduction_methods/pole_matching/reduced_basis_pole_matching.py
similarity index 96%
rename from rrompy/reduction_methods/pole_matching/rb_distributed_pole_matching.py
rename to rrompy/reduction_methods/pole_matching/reduced_basis_pole_matching.py
index d211dd3..d3c7aa4 100644
--- a/rrompy/reduction_methods/pole_matching/rb_distributed_pole_matching.py
+++ b/rrompy/reduction_methods/pole_matching/reduced_basis_pole_matching.py
@@ -1,190 +1,190 @@
# 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 numpy import isinf
from copy import deepcopy as copy
-from rrompy.reduction_methods.pivoting.rb_distributed_pivoted import \
- RBDistributedPivoted
+from rrompy.reduction_methods.pivoting.reduced_basis_pivoted import \
+ ReducedBasisPivoted
from .generic_pole_matching_approximant import GenericPoleMatchingApproximant
from rrompy.reduction_methods.trained_model import (TrainedModelPivotedData,
- TrainedModelPivotedRBPoleMatching as tModel)
+ TrainedModelPoleMatchingReducedBasis as tModel)
from rrompy.utilities.base import verbosityManager as vbMng
from rrompy.utilities.exception_manager import RROMPyAssert
-__all__ = ['RBDistributedPoleMatching']
+__all__ = ['ReducedBasisPoleMatching']
-class RBDistributedPoleMatching(GenericPoleMatchingApproximant,
- RBDistributedPivoted):
+class ReducedBasisPoleMatching(GenericPoleMatchingApproximant,
+ ReducedBasisPivoted):
"""
ROM pivoted RB approximant computation for parametric problems with pole
matching.
Args:
HFEngine: HF problem solver.
mu0(optional): Default parameter. Defaults to 0.
directionPivot(optional): Pivot components. 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;
- 'matchingWeight': weight for pole matching optimization; defaults
to 1;
- 'cutOffTolerance': tolerance for ignoring parasitic poles;
defaults to np.inf;
- 'cutOffType': rule for tolerance computation for parasitic poles;
defaults to 'MAGNITUDE';
- 'S': total number of pivot samples current approximant relies
upon;
- 'samplerPivot': pivot sample point generator;
- 'SMarginal': total number of marginal samples current approximant
relies upon;
- 'samplerMarginal': marginal sample point generator;
- 'R': rank for pivot Galerkin projection; defaults to prod(S);
- 'PODTolerance': tolerance for pivot snapshots POD; defaults to
-1;
- 'polybasisMarginal': type of polynomial basis for marginal
interpolation; allowed values include 'MONOMIAL', 'CHEBYSHEV'
and 'LEGENDRE'; defaults to 'MONOMIAL';
- 'MMarginal': degree of marginal interpolant; defaults to 0;
- 'radialBasisMarginal': radial basis family for marginal
interpolant; defaults to 0, i.e. no radial basis;
- 'radialBasisWeightsMarginal': radial basis weights for marginal
interpolant; defaults to 0, i.e. identity;
- 'interpRcondMarginal': tolerance for marginal interpolation;
defaults to None.
Defaults to empty dict.
homogeneized(optional): Whether to homogeneize Dirichlet BCs. Defaults
to False.
verbosity(optional): Verbosity level. Defaults to 10.
Attributes:
HFEngine: HF problem solver.
mu0: Default parameter.
directionPivot: Pivot components.
mus: Array of snapshot parameters.
musPivot: Array of pivot snapshot parameters.
musMarginal: Array of marginal snapshot parameters.
homogeneized: Whether to homogeneize Dirichlet BCs.
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;
- 'matchingWeight': weight for pole matching optimization;
- 'cutOffTolerance': tolerance for ignoring parasitic poles;
- 'cutOffType': rule for tolerance computation for parasitic poles;
- 'R': rank for Galerkin projection;
- 'PODTolerance': tolerance for snapshots POD;
- 'polybasisMarginal': type of polynomial basis for marginal
interpolation;
- 'MMarginal': degree of marginal interpolant;
- 'radialBasisMarginal': radial basis family for marginal
interpolant;
- 'radialBasisWeightsMarginal': radial basis weights for marginal
interpolant;
- 'interpRcondMarginal': tolerance for marginal interpolation.
parameterListCritical: Recognized keys of critical approximant
parameters:
- 'S': total number of pivot samples current approximant relies
upon;
- 'samplerPivot': pivot sample point generator;
- 'SMarginal': total number of marginal samples current approximant
relies upon;
- 'samplerMarginal': marginal sample point generator.
POD: Whether to compute POD of snapshots.
matchingWeight: Weight for pole matching optimization.
cutOffTolerance: Tolerance for ignoring parasitic poles.
cutOffType: Rule for tolerance computation for parasitic poles.
S: Total number of pivot samples current approximant relies upon.
samplerPivot: Pivot sample point generator.
SMarginal: Total number of marginal samples current approximant relies
upon.
samplerMarginal: Marginal sample point generator.
R: Rank for Galerkin projection.
PODTolerance: Tolerance for pivot snapshots POD.
polybasisMarginal: Type of polynomial basis for marginal interpolation.
MMarginal: Degree of marginal interpolant.
radialBasisMarginal: Radial basis family for marginal interpolant.
radialBasisWeightsMarginal: Radial basis weights for marginal
interpolant.
interpRcondMarginal: Tolerance for marginal interpolation.
muBoundsPivot: list of bounds for pivot parameter values.
muBoundsMarginal: list of bounds for marginal parameter values.
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.
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.
As: List of sparse matrices (in CSC format) representing coefficients
of linear system matrix.
bs: List of numpy vectors representing coefficients of linear system
RHS.
ARBs: List of sparse matrices (in CSC format) representing coefficients
of compressed linear system matrix.
bRBs: List of numpy vectors representing coefficients of compressed
linear system RHS.
"""
def setupApprox(self):
"""Compute RB projection matrix."""
if self.checkComputedApprox():
return
RROMPyAssert(self._mode, message = "Cannot setup approximant.")
vbMng(self, "INIT", "Setting up {}.".format(self.name()), 5)
self.computeScaleFactor()
self.computeSnapshots()
self._setupAffineBlocks()
pMat = self._setupProjectionMatrix()
ARBsList, bRBsList, pList = self.assembleReducedSystemMarginal(pMat)
if self.trainedModel is None:
self.trainedModel = tModel()
self.trainedModel.verbosity = self.verbosity
self.trainedModel.timestamp = self.timestamp
data = TrainedModelPivotedData(self.trainedModel.name(), self.mu0,
pList, self.scaleFactor,
self.HFEngine.rescalingExp,
self.directionPivot)
data.musPivot = copy(self.musPivot)
data.musMarginal = copy(self.musMarginal)
self.trainedModel.data = data
else:
self.trainedModel = self.trainedModel
self.trainedModel.data.projMat = copy(pList)
self.trainedModel.data.marginalInterp = self._setupMarginalInterp()
self.trainedModel.data.ARBsList = ARBsList
self.trainedModel.data.bRBsList = bRBsList
vbMng(self, "INIT", "Matching poles.", 10)
self.trainedModel.initializeFromAffine(self.HFEngine,
self.matchingWeight, self.POD)
vbMng(self, "DEL", "Done matching poles.", 10)
if not isinf(self.cutOffTolerance):
vbMng(self, "INIT", "Recompressing by cut-off.", 10)
msg = self.trainedModel.recompressByCutOff([- 1., 1.],
self.cutOffTolerance,
self.cutOffType)
vbMng(self, "DEL", "Done recompressing." + msg, 10)
self.trainedModel.data.approxParameters = copy(self.approxParameters)
vbMng(self, "DEL", "Done setting up approximant.", 5)
diff --git a/rrompy/reduction_methods/trained_model/__init__.py b/rrompy/reduction_methods/trained_model/__init__.py
index e16c89a..a80b481 100644
--- a/rrompy/reduction_methods/trained_model/__init__.py
+++ b/rrompy/reduction_methods/trained_model/__init__.py
@@ -1,45 +1,47 @@
# 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 .trained_model import TrainedModel
from .trained_model_data import TrainedModelData
from .trained_model_rational import TrainedModelRational
-from .trained_model_rb import TrainedModelRB
+from .trained_model_reduced_basis import TrainedModelReducedBasis
from .trained_model_pivoted_data import TrainedModelPivotedData
from .trained_model_pivoted_rational import TrainedModelPivotedRational
-from .trained_model_pivoted_rb import TrainedModelPivotedRB
-from .trained_model_pivoted_general_pole_matching import \
- TrainedModelPivotedGeneralPoleMatching
-from .trained_model_pivoted_rational_pole_matching import \
- TrainedModelPivotedRationalPoleMatching
-from .trained_model_pivoted_rb_pole_matching import \
- TrainedModelPivotedRBPoleMatching
+from .trained_model_pivoted_reduced_basis import \
+ TrainedModelPivotedReducedBasis
+from .trained_model_pole_matching_general import \
+ TrainedModelPoleMatchingGeneral
+from .trained_model_pole_matching_rational import \
+ TrainedModelPoleMatchingRational
+from .trained_model_pole_matching_reduced_basis import \
+ TrainedModelPoleMatchingReducedBasis
__all__ = [
'TrainedModel',
'TrainedModelData',
'TrainedModelRational',
- 'TrainedModelRB',
+ 'TrainedModelReducedBasis',
'TrainedModelPivotedData',
'TrainedModelPivotedRational',
- 'TrainedModelPivotedRB',
- 'TrainedModelPivotedRationalPoleMatching',
- 'TrainedModelPivotedRBPoleMatching'
+ 'TrainedModelPivotedReducedBasis',
+ 'TrainedModelPoleMatchingGeneral',
+ 'TrainedModelPoleMatchingRational',
+ 'TrainedModelPoleMatchingReducedBasis'
]
diff --git a/rrompy/reduction_methods/trained_model/trained_model_pivoted_rb.py b/rrompy/reduction_methods/trained_model/trained_model_pivoted_reduced_basis.py
similarity index 97%
rename from rrompy/reduction_methods/trained_model/trained_model_pivoted_rb.py
rename to rrompy/reduction_methods/trained_model/trained_model_pivoted_reduced_basis.py
index 33c669d..490a1b0 100644
--- a/rrompy/reduction_methods/trained_model/trained_model_pivoted_rb.py
+++ b/rrompy/reduction_methods/trained_model/trained_model_pivoted_reduced_basis.py
@@ -1,202 +1,202 @@
# Copyright (C) 2018 by the RROMPy authors
#
# This file is part of RROMPy.
#
# RROMPy is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# RROMPy is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with RROMPy. If not, see .
#
import numpy as np
-from .trained_model_rb import TrainedModelRB
+from .trained_model_reduced_basis import TrainedModelReducedBasis
from rrompy.utilities.base.types import (Np1D, ListAny, paramVal, paramList,
sampList)
from rrompy.utilities.base import verbosityManager as vbMng, freepar as fp
from rrompy.utilities.poly_fitting.polynomial import (
hashIdxToDerivative as hashI)
from rrompy.utilities.numerical import (eigvalsNonlinearDense,
marginalizePolyList)
from rrompy.utilities.exception_manager import RROMPyException
from rrompy.parameter import checkParameter, checkParameterList
from rrompy.sampling import emptySampleList, sampleList
-__all__ = ['TrainedModelPivotedRB']
+__all__ = ['TrainedModelPivotedReducedBasis']
-class TrainedModelPivotedRB(TrainedModelRB):
+class TrainedModelPivotedReducedBasis(TrainedModelReducedBasis):
"""
ROM approximant evaluation for pivoted RB approximant.
Attributes:
Data: dictionary with all that can be pickled.
"""
def centerNormalizePivot(self, mu : paramList = [],
mu0 : paramVal = None) -> paramList:
"""
Compute normalized parameter to be plugged into approximant.
Args:
mu: Parameter(s) 1.
mu0: Parameter(s) 2. If None, set to self.data.mu0.
Returns:
Normalized parameter.
"""
mu = checkParameterList(mu, self.data.nparPivot)[0]
if mu0 is None: mu0 = self.data.mu0Pivot
rad = ((mu ** self.data.rescalingExpPivot
- mu0 ** self.data.rescalingExpPivot)
/ self.data.scaleFactorPivot)
return rad
def centerNormalizeMarginal(self, mu : paramList = [],
mu0 : paramVal = None) -> paramList:
"""
Compute normalized parameter to be plugged into approximant.
Args:
mu: Parameter(s) 1.
mu0: Parameter(s) 2. If None, set to self.data.mu0.
Returns:
Normalized parameter.
"""
mu = checkParameterList(mu, self.data.nparMarginal)[0]
if mu0 is None: mu0 = self.data.mu0Marginal
rad = ((mu ** self.data.rescalingExpMarginal
- mu0 ** self.data.rescalingExpMarginal)
/ self.data.scaleFactorMarginal)
return rad
def interpolateMarginal(self, mu : paramVal = [],
samples : ListAny = []) -> ListAny:
"""
Evaluate marginal interpolator at arbitrary marginal parameter.
Args:
mu: Target parameter.
samples: Objects to interpolate.
"""
mu = checkParameter(mu, self.data.nparMarginal)
sList = isinstance(samples[0], sampleList)
sEff = [None] * len(samples)
for j in range(len(samples)):
if sList:
sEff[j] = samples[j].data
else:
sEff[j] = samples[j]
vbMng(self, "INIT", "Interpolating marginal at mu = {}.".format(mu),
95)
muC = self.centerNormalizeMarginal(mu)
p = np.zeros_like(sEff[0], dtype = sEff[0].dtype)
for mIj, spj in zip(self.data.marginalInterp, sEff):
p += mIj(muC) * spj
vbMng(self, "DEL", "Done interpolating marginal.", 95)
# if sList: p = sampleList(p)
return p
def getApproxReduced(self, mu : paramList = []) -> sampList:
"""
Evaluate reduced representation of approximant at arbitrary parameter.
Args:
mu: Target parameter.
"""
mu = checkParameterList(mu, self.data.npar)[0]
if (not hasattr(self, "lastSolvedApproxReduced")
or self.lastSolvedApproxReduced != mu):
vbMng(self, "INIT",
"Computing RB solution at mu = {}.".format(mu), 12)
self.uApproxReduced = emptySampleList()
self.uApproxReduced.reset((self.data.ARBsList[0][0].shape[0],
len(mu)), self.data.projMat[0].dtype)
for i, muPL in enumerate(mu):
muP = checkParameter([muPL(0, x) \
for x in self.data.directionPivot])
muM = [muPL(0, x) for x in self.data.directionMarginal]
vbMng(self, "INIT",
"Assembling reduced model for mu = {}.".format(muPL), 87)
ARBs = self.interpolateMarginal(muM, self.data.ARBsList)
bRBs = self.interpolateMarginal(muM, self.data.bRBsList)
ARBmu = ARBs[0]
bRBmu = bRBs[0]
for j in range(1, max(len(ARBs), len(bRBs))):
theta = np.prod(self.centerNormalizePivot(muP)
** hashI(j, self.data.nparPivot))
if j < len(ARBs):
ARBmu += theta * ARBs[j]
if j < len(bRBs):
bRBmu += theta * bRBs[j]
vbMng(self, "DEL", "Done assembling reduced model.", 87)
vbMng(self, "INIT",
"Solving reduced model for mu = {}.".format(muPL), 75)
self.uApproxReduced[i] = np.linalg.solve(ARBmu, bRBmu)
vbMng(self, "DEL", "Done solving reduced model.", 75)
vbMng(self, "DEL", "Done computing RB solution.", 12)
self.lastSolvedApproxReduced = mu
return self.uApproxReduced
def getApprox(self, mu : paramList = []) -> sampList:
"""
Evaluate approximant at arbitrary parameter.
Args:
mu: Target parameter.
"""
mu = checkParameterList(mu, self.data.npar)[0]
if (not hasattr(self, "lastSolvedApprox")
or self.lastSolvedApprox != mu):
uApproxR = self.getApproxReduced(mu)
self.uApprox = emptySampleList()
self.uApprox.reset((self.data.projMat[0].shape[0], len(mu)),
self.data.projMat[0].dtype)
for i in range(len(mu)):
muM = [mu(i, x) for x in self.data.directionMarginal]
projLoc = self.interpolateMarginal(muM, self.data.projMat)[0]
self.uApprox[i] = projLoc.dot(uApproxR[i])
self.lastSolvedApprox = mu
return self.uApprox
def getPoles(self, marginalVals : ListAny = [fp], jSupp : int = 1,
**kwargs) -> Np1D:
"""
Obtain approximant poles.
Returns:
Numpy complex vector of poles.
"""
if not hasattr(marginalVals, "__len__"): marginalVals = [marginalVals]
marginalVals = list(marginalVals)
try:
rDim = marginalVals.index(fp)
if rDim < len(marginalVals) - 1 and fp in marginalVals[rDim + 1 :]:
raise
except:
raise RROMPyException(("Exactly 1 'freepar' entry in "
"marginalVals must be provided."))
if rDim in self.data.directionMarginal:
raise RROMPyException(("'freepar' entry in marginalVals must be "
"in pivot dimension."))
mValsP = [marginalVals[x] for x in self.data.directionPivot]
rDimP = mValsP.index(fp)
mValsM = [marginalVals[x] for x in self.data.directionMarginal]
ARBs = self.interpolateMarginal(mValsM, self.data.ARBsList)
if self.data.nparPivot > 1:
mValsP[rDimP] = self.data.mu0Pivot(rDimP)
mValsP = self.centerNormalizePivot(mValsP)
mValsP = mValsP.data.flatten()
mValsP[rDimP] = fp
ARBs = marginalizePolyList(ARBs, mValsP, "auto")
ev = eigvalsNonlinearDense(ARBs, jSupp = jSupp, **kwargs)
return np.power(self.data.mu0(rDim) ** self.data.rescalingExp[rDim]
+ ev * self.scaleFactor[rDim],
1. / self.data.rescalingExp[rDim])
diff --git a/rrompy/reduction_methods/trained_model/trained_model_pivoted_general_pole_matching.py b/rrompy/reduction_methods/trained_model/trained_model_pole_matching_general.py
similarity index 98%
rename from rrompy/reduction_methods/trained_model/trained_model_pivoted_general_pole_matching.py
rename to rrompy/reduction_methods/trained_model/trained_model_pole_matching_general.py
index fc1720c..3f088c5 100644
--- a/rrompy/reduction_methods/trained_model/trained_model_pivoted_general_pole_matching.py
+++ b/rrompy/reduction_methods/trained_model/trained_model_pole_matching_general.py
@@ -1,244 +1,244 @@
# Copyright (C) 2018 by the RROMPy authors
#
# This file is part of RROMPy.
#
# RROMPy is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# RROMPy is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with RROMPy. If not, see .
#
import numpy as np
from .trained_model import TrainedModel
from rrompy.utilities.base.types import (Np1D, Tuple, ListAny, paramVal,
paramList, sampList, HFEng)
from rrompy.utilities.base import verbosityManager as vbMng, freepar as fp
from rrompy.utilities.numerical import pointMatching
from rrompy.utilities.poly_fitting.heaviside import HeavisideInterpolator as HI
from rrompy.utilities.exception_manager import RROMPyException, RROMPyAssert
from rrompy.parameter import checkParameter, checkParameterList
from rrompy.sampling import emptySampleList
-__all__ = ['TrainedModelPivotedGeneralPoleMatching']
+__all__ = ['TrainedModelPoleMatchingGeneral']
-class TrainedModelPivotedGeneralPoleMatching(TrainedModel):
+class TrainedModelPoleMatchingGeneral(TrainedModel):
"""
ROM approximant evaluation for pivoted approximants with pole matching.
Attributes:
Data: dictionary with all that can be pickled.
"""
def initializeFromLists(self, poles:ListAny, coeffs:ListAny, basis:str,
HFEngine:HFEng, matchingWeight : float = 1.,
POD : bool = True):
"""Initialize Heaviside representation."""
musM = self.data.musMarginal
margAbsDist = np.sum(np.abs(np.tile(musM.data.reshape(len(musM),1,-1),
[1, len(musM), 1])
- np.tile(musM.data.reshape(1,len(musM),-1),
[len(musM), 1, 1])), axis = 2)
N = len(poles[0])
explored = [0]
unexplored = list(range(1, len(musM)))
for _ in range(1, len(musM)):
minIdx = np.argmin(np.concatenate([margAbsDist[ex, unexplored] \
for ex in explored]))
minIex = explored[minIdx // len(unexplored)]
minIunex = unexplored[minIdx % len(unexplored)]
dist = np.abs(np.tile(poles[minIex].reshape(-1, 1), N)
- poles[minIunex].reshape(1, -1))
if matchingWeight != 0:
resex = coeffs[minIex][: N]
resunex = coeffs[minIunex][: N]
if POD:
distR = resex.dot(resunex.T.conj())
distR = (distR.T / np.linalg.norm(resex, axis = 1)).T
distR = distR / np.linalg.norm(resunex, axis = 1)
else:
resex = self.data.projMat.dot(resex.T)
resunex = self.data.projMat.dot(resunex.T)
distR = HFEngine.innerProduct(resex, resunex).T
distR = (distR.T / HFEngine.norm(resex)).T
distR = distR / HFEngine.norm(resunex)
distR = np.abs(distR)
distR[distR > 1.] = 1.
dist += 2. / np.pi * matchingWeight * np.arccos(distR)
reordering = pointMatching(dist, verbObj = self)
poles[minIunex] = poles[minIunex][reordering]
coeffs[minIunex][: N] = coeffs[minIunex][reordering]
explored += [minIunex]
unexplored.remove(minIunex)
# explored = np.append(explored, minIunex)
# unexplored = unexplored[unexplored != minIunex]
HIs = []
for pls, cfs in zip(poles, coeffs):
hsi = HI()
hsi.poles = pls
hsi.coeffs = cfs
hsi.npar = 1
hsi.polybasis = basis
HIs += [hsi]
self.data.HIs = HIs
def recompressByCutOff(self, murange : Tuple[float, float] = [- 1., 1.],
tol : float = np.inf, rtype : str = "MAGNITUDE"):
if np.isinf(tol): return " No poles erased."
N = len(self.data.HIs[0].poles)
mu0 = np.mean(murange)
musig = murange[0] - mu0
if np.isclose(musig, 0.):
radius = lambda x: np.abs(x - mu0)
else:
if rtype == "MAGNITUDE":
murdir = (murange[0] - mu0) / np.abs(musig)
def radius(x):
scalprod = (x - mu0) * murdir.conj() / np.abs(musig)
rescalepar = np.abs(np.real(scalprod))
rescaleort = np.abs(np.imag(scalprod))
return ((rescalepar - 1.) ** 2. * (rescalepar > 1.)
+ rescaleort ** 2.) ** .5
else:#if rtype == "POTENTIAL":
def radius(x):
rescale = (x - mu0) / musig
return np.max(np.abs(rescale * np.array([-1., 1.])
+ (rescale ** 2. - 1) ** .5)) - 1.
keepPole, removePole = [], []
for j in range(N):
for hi in self.data.HIs:
if radius(hi.poles[j]) <= tol:
keepPole += [j]
break
if len(keepPole) == 0 or keepPole[-1] != j: removePole += [j]
if len(keepPole) == N: return " No poles erased."
keepCoeff = keepPole + [N]
keepCoeff = keepCoeff + list(range(N + 1,len(self.data.HIs[0].coeffs)))
for hi in self.data.HIs:
polyCorrection = np.zeros_like(hi.coeffs[0, :])
for j in removePole:
polyCorrection += hi.coeffs[j, :] / (mu0 - hi.poles[j])
if len(hi.coeffs) == N:
hi.coeffs = np.vstack((hi.coeffs, polyCorrection))
else:
hi.coeffs[N, :] += polyCorrection
hi.poles = hi.poles[keepPole]
hi.coeffs = hi.coeffs[keepCoeff, :]
return (" Erased {} pole".format(len(removePole))
+ "s" * (len(removePole) > 1) + ".")
def interpolateMarginalInterpolator(self, mu : paramVal = []) -> Np1D:
"""Obtain interpolated approximant interpolator."""
mu = checkParameter(mu, self.data.nparMarginal)[0]
hsi = HI()
hsi.poles = self.interpolateMarginalPoles(mu)
hsi.coeffs = self.interpolateMarginalCoeffs(mu)
hsi.npar = 1
hsi.polybasis = self.data.HIs[0].polybasis
return hsi
def interpolateMarginalPoles(self, mu : paramList = []) -> Np1D:
"""Obtain interpolated approximant poles."""
mu = checkParameterList(mu, self.data.nparMarginal)[0]
return self.interpolateMarginal(mu, [hi.poles for hi in self.data.HIs])
def interpolateMarginalCoeffs(self, mu : paramList = []) -> Np1D:
"""Obtain interpolated approximant coefficients."""
mu = checkParameterList(mu, self.data.nparMarginal)[0]
cs = self.interpolateMarginal(mu, [hi.coeffs for hi in self.data.HIs])
if isinstance(cs, (list, tuple,)): cs = np.array(cs)
return cs.reshape(self.data.HIs[0].coeffs.shape)
def getApproxReduced(self, mu : paramList = []) -> sampList:
"""
Evaluate reduced representation of approximant at arbitrary parameter.
Args:
mu: Target parameter.
"""
RROMPyAssert(self.data.nparPivot, 1, "Number of pivot parameters")
mu = checkParameterList(mu, self.data.npar)[0]
if (not hasattr(self, "lastSolvedApproxReduced")
or self.lastSolvedApproxReduced != mu):
vbMng(self, "INIT",
"Evaluating approximant at mu = {}.".format(mu), 12)
self.uApproxReduced = emptySampleList()
self.uApproxReduced.reset((self.data.HIs[0].coeffs.shape[1],
len(mu)), self.data.projMat[0].dtype)
for i, muPL in enumerate(mu):
muL = self.centerNormalizePivot([muPL(0, x) \
for x in self.data.directionPivot])
muM = [muPL(0, x) for x in self.data.directionMarginal]
vbMng(self, "INIT",
"Assembling reduced model for mu = {}.".format(muPL), 87)
hsL = self.interpolateMarginalInterpolator(muM)
vbMng(self, "DEL", "Done assembling reduced model.", 87)
self.uApproxReduced[i] = hsL(muL)
vbMng(self, "DEL", "Done evaluating approximant.", 12)
self.lastSolvedApproxReduced = mu
return self.uApproxReduced
def getPoles(self, *args, **kwargs) -> Np1D:
"""
Obtain approximant poles.
Returns:
Numpy complex vector of poles.
"""
RROMPyAssert(self.data.nparPivot, 1, "Number of pivot parameters")
if len(args) + len(kwargs) > 1:
raise RROMPyException(("Wrong number of parameters passed. "
"Only 1 available."))
elif len(args) + len(kwargs) == 1:
if len(args) == 1:
mVals = args[0]
else:
mVals = kwargs["marginalVals"]
if not hasattr(mVals, "__len__"): mVals = [mVals]
mVals = list(mVals)
else:
mVals = [fp]
try:
rDim = mVals.index(fp)
if rDim < len(mVals) - 1 and fp in mVals[rDim + 1 :]:
raise
except:
raise RROMPyException(("Exactly 1 'freepar' entry in "
"marginalVals must be provided."))
if rDim != self.data.directionPivot[0]:
raise RROMPyException(("'freepar' entry in marginalVals must "
"coincide with pivot direction."))
mVals[rDim] = self.data.mu0(rDim)
mMarg = [mVals[j] for j in range(len(mVals)) if j != rDim]
roots = np.array(self.interpolateMarginalPoles(mMarg))
return np.power(self.data.mu0(rDim) ** self.data.rescalingExp[rDim]
+ self.data.scaleFactor[rDim] * roots,
1. / self.data.rescalingExp[rDim])
def getResidues(self, *args, **kwargs) -> Np1D:
"""
Obtain approximant residues.
Returns:
Numpy matrix with residues as columns.
"""
pls = self.getPoles(*args, **kwargs)
if len(args) == 1:
mVals = args[0]
else:
mVals = kwargs["marginalVals"]
if not hasattr(mVals, "__len__"): mVals = [mVals]
mVals = list(mVals)
rDim = mVals.index(fp)
mMarg = [mVals[j] for j in range(len(mVals)) if j != rDim]
residues = self.interpolateMarginalCoeffs(mMarg)[: len(pls)]
res = self.data.projMat.dot(residues)
return pls, res
diff --git a/rrompy/reduction_methods/trained_model/trained_model_pivoted_rational_pole_matching.py b/rrompy/reduction_methods/trained_model/trained_model_pole_matching_rational.py
similarity index 92%
rename from rrompy/reduction_methods/trained_model/trained_model_pivoted_rational_pole_matching.py
rename to rrompy/reduction_methods/trained_model/trained_model_pole_matching_rational.py
index a7e9fba..c24171f 100644
--- a/rrompy/reduction_methods/trained_model/trained_model_pivoted_rational_pole_matching.py
+++ b/rrompy/reduction_methods/trained_model/trained_model_pole_matching_rational.py
@@ -1,110 +1,109 @@
# Copyright (C) 2018 by the RROMPy authors
#
# This file is part of RROMPy.
#
# RROMPy is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# RROMPy is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with RROMPy. If not, see .
#
import numpy as np
from scipy.special import factorial as fact
from itertools import combinations
-from .trained_model_pivoted_general_pole_matching import \
- TrainedModelPivotedGeneralPoleMatching
+from .trained_model_pole_matching_general import \
+ TrainedModelPoleMatchingGeneral
from .trained_model_pivoted_rational import TrainedModelPivotedRational
from rrompy.utilities.base.types import Np1D, List, paramList, sampList, HFEng
from rrompy.utilities.poly_fitting.heaviside import rational2heaviside
from rrompy.utilities.exception_manager import RROMPyException, RROMPyAssert
from rrompy.parameter import checkParameterList
from rrompy.sampling import emptySampleList
-__all__ = ['TrainedModelPivotedRationalPoleMatching']
+__all__ = ['TrainedModelPoleMatchingRational']
-class TrainedModelPivotedRationalPoleMatching(
- TrainedModelPivotedGeneralPoleMatching,
- TrainedModelPivotedRational):
+class TrainedModelPoleMatchingRational(TrainedModelPoleMatchingGeneral,
+ TrainedModelPivotedRational):
"""
ROM approximant evaluation for pivoted Rational approximant with pole
matching.
Attributes:
Data: dictionary with all that can be pickled.
"""
def initializeFromRational(self, HFEngine:HFEng,
matchingWeight : float = 1., POD : bool = True):
"""Initialize Heaviside representation."""
RROMPyAssert(self.data.nparPivot, 1, "Number of pivot parameters")
poles, coeffs = [], []
for Q, P in zip(self.data.Qs, self.data.Ps):
cfs, pls, basis = rational2heaviside(P, Q)
poles += [pls]
coeffs += [cfs]
self.initializeFromLists(poles, coeffs, basis, HFEngine,
matchingWeight, POD)
self.data._temporary = False
def getPVal(self, mu : paramList = []) -> sampList:
"""
Evaluate rational numerator at arbitrary parameter.
Args:
mu: Target parameter.
"""
if self.data._temporary: return super().getPVal(mu)
RROMPyAssert(self.data.nparPivot, 1, "Number of pivot parameters")
mu = checkParameterList(mu, self.data.npar)[0]
p = emptySampleList()
p.reset((len(self.data.HIs[0].coeffs.shape[1]), len(mu)))
for i, muPL in enumerate(mu):
muL = self.centerNormalizePivot([muPL(0, x) \
for x in self.data.directionPivot])
muM = [muPL(0, x) for x in self.data.directionMarginal]
hsL = self.interpolateMarginalInterpolator(muM)
p[i] = hsL(muL) * np.prod(muL(0, 0) - hsL.poles)
return p
def getQVal(self, mu:Np1D, der : List[int] = None,
scl : Np1D = None) -> Np1D:
"""
Evaluate rational denominator at arbitrary parameter.
Args:
mu: Target parameter.
der(optional): Derivatives to take before evaluation.
"""
if self.data._temporary: return super().getQVal(mu, der, scl)
RROMPyAssert(self.data.nparPivot, 1, "Number of pivot parameters")
mu = checkParameterList(mu, self.data.npar)[0]
muP = self.centerNormalizePivot(checkParameterList(
mu.data[:, self.data.directionPivot], 1)[0])
muM = checkParameterList(mu.data[:, self.data.directionMarginal],
self.data.npar - 1)[0]
if der is None:
derP, derM = 0, [0]
else:
derP = der[self.data.directionPivot[0]]
derM = [der[x] for x in self.data.directionMarginal]
if np.any(np.array(derM) != 0):
raise RROMPyException(("Derivatives of Q with respect to marginal "
"parameters not allowed."))
sclP = 1 if scl is None else scl[self.data.directionPivot[0]]
derVal = np.zeros(len(mu), dtype = self.data.HIs[0].poles.dtype)
N = len(self.data.HIs[0].poles)
if derP == N: derVal[:] = 1.
elif derP >= 0 and derP < N:
pls = self.interpolateMarginalPoles(muM).reshape(-1, len(mu)).T
plsDist = muP.data.reshape(-1, 1) - pls
for terms in combinations(np.arange(N), N - derP):
derVal += np.prod(plsDist[:, list(terms)], axis = 1)
return sclP ** derP * fact(derP) * derVal
diff --git a/rrompy/reduction_methods/trained_model/trained_model_pivoted_rb_pole_matching.py b/rrompy/reduction_methods/trained_model/trained_model_pole_matching_reduced_basis.py
similarity index 81%
rename from rrompy/reduction_methods/trained_model/trained_model_pivoted_rb_pole_matching.py
rename to rrompy/reduction_methods/trained_model/trained_model_pole_matching_reduced_basis.py
index f2cf22d..23802b1 100644
--- a/rrompy/reduction_methods/trained_model/trained_model_pivoted_rb_pole_matching.py
+++ b/rrompy/reduction_methods/trained_model/trained_model_pole_matching_reduced_basis.py
@@ -1,55 +1,56 @@
# Copyright (C) 2018 by the RROMPy authors
#
# This file is part of RROMPy.
#
# RROMPy is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# RROMPy is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with RROMPy. If not, see .
#
import numpy as np
-from .trained_model_pivoted_general_pole_matching import \
- TrainedModelPivotedGeneralPoleMatching
-from .trained_model_pivoted_rb import TrainedModelPivotedRB
+from .trained_model_pole_matching_general import \
+ TrainedModelPoleMatchingGeneral
+from .trained_model_pivoted_reduced_basis import \
+ TrainedModelPivotedReducedBasis
from rrompy.utilities.base.types import HFEng
from rrompy.utilities.poly_fitting.heaviside import affine2heaviside
from rrompy.utilities.exception_manager import RROMPyAssert
-__all__ = ['TrainedModelPivotedRBPoleMatching']
+__all__ = ['TrainedModelPoleMatchingReducedBasis']
-class TrainedModelPivotedRBPoleMatching(TrainedModelPivotedGeneralPoleMatching,
- TrainedModelPivotedRB):
+class TrainedModelPoleMatchingReducedBasis(TrainedModelPoleMatchingGeneral,
+ TrainedModelPivotedReducedBasis):
"""
ROM approximant evaluation for pivoted RB approximant with pole matching.
Attributes:
Data: dictionary with all that can be pickled.
"""
def initializeFromAffine(self, HFEngine:HFEng, matchingWeight : float = 1.,
POD : bool = True, jSupp : int = 1):
"""Initialize Heaviside representation."""
RROMPyAssert(self.data.nparPivot, 1, "Number of pivot parameters")
poles, coeffs = [], []
nbadpls = 0
for As, bs in zip(self.data.ARBsList, self.data.bRBsList):
cfs, pls, basis = affine2heaviside(As, bs, jSupp)
poles += [pls]
coeffs += [cfs]
nbadpls = max(nbadpls, np.sum(np.isinf(pls)))
if nbadpls > 0:
for j in range(len(poles)):
plsgood = np.argsort(np.abs(pls))[: - nbadpls]
poles[j] = poles[j][plsgood]
coeffs[j] = coeffs[j][plsgood, :]
self.initializeFromLists(poles, coeffs, basis, HFEngine,
matchingWeight, POD)
diff --git a/rrompy/reduction_methods/trained_model/trained_model_rb.py b/rrompy/reduction_methods/trained_model/trained_model_reduced_basis.py
similarity index 98%
rename from rrompy/reduction_methods/trained_model/trained_model_rb.py
rename to rrompy/reduction_methods/trained_model/trained_model_reduced_basis.py
index bd5ad7f..6754e66 100644
--- a/rrompy/reduction_methods/trained_model/trained_model_rb.py
+++ b/rrompy/reduction_methods/trained_model/trained_model_reduced_basis.py
@@ -1,108 +1,108 @@
# Copyright (C) 2018 by the RROMPy authors
#
# This file is part of RROMPy.
#
# RROMPy is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# RROMPy is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with RROMPy. If not, see .
#
import numpy as np
from .trained_model import TrainedModel
from rrompy.utilities.base.types import Np1D, ListAny, paramList, sampList
from rrompy.utilities.base import verbosityManager as vbMng, freepar as fp
from rrompy.utilities.poly_fitting.polynomial import (
hashIdxToDerivative as hashI)
from rrompy.utilities.numerical import (eigvalsNonlinearDense,
marginalizePolyList)
from rrompy.utilities.exception_manager import RROMPyException
from rrompy.parameter import checkParameter, checkParameterList
from rrompy.sampling import emptySampleList
-__all__ = ['TrainedModelRB']
+__all__ = ['TrainedModelReducedBasis']
-class TrainedModelRB(TrainedModel):
+class TrainedModelReducedBasis(TrainedModel):
"""
ROM approximant evaluation for RB approximant.
Attributes:
Data: dictionary with all that can be pickled.
"""
def getApproxReduced(self, mu : paramList = []) -> sampList:
"""
Evaluate reduced representation of approximant at arbitrary parameter.
Args:
mu: Target parameter.
"""
mu = checkParameterList(mu, self.data.npar)[0]
if (not hasattr(self, "lastSolvedApproxReduced")
or self.lastSolvedApproxReduced != mu):
vbMng(self, "INIT",
"Computing RB solution at mu = {}.".format(mu), 12)
ARBs, bRBs = self.data.ARBs, self.data.bRBs
self.uApproxReduced = emptySampleList()
self.uApproxReduced.reset((ARBs[0].shape[0], len(mu)),
self.data.projMat.dtype)
mu0Eff = (self.data.mu0 ** self.data.rescalingExp).data
for i in range(len(mu)):
vbMng(self, "INIT", "Assembling reduced model for mu = {}."\
.format(mu[i]), 17)
ARBmu = 1. * ARBs[0]
bRBmu = 1. * bRBs[0]
for j in range(1, max(len(ARBs), len(bRBs))):
derIdx = hashI(j, self.data.npar)
theta = np.prod(((mu[i] ** self.data.rescalingExp - mu0Eff)
/ self.data.scaleFactor) ** derIdx)
if j < len(ARBs):
ARBmu += theta * ARBs[j]
if j < len(bRBs):
bRBmu += theta * bRBs[j]
vbMng(self, "DEL", "Done assembling reduced model.", 17)
vbMng(self, "INIT",
"Solving reduced model for mu = {}.".format(mu[i]), 15)
self.uApproxReduced[i] = np.linalg.solve(ARBmu, bRBmu)
vbMng(self, "DEL", "Done solving reduced model.", 15)
vbMng(self, "DEL", "Done computing RB solution.", 12)
self.lastSolvedApproxReduced = mu
return self.uApproxReduced
def getPoles(self, marginalVals : ListAny = [fp], jSupp : int = 1,
**kwargs) -> Np1D:
"""
Obtain approximant poles.
Returns:
Numpy complex vector of poles.
"""
if not hasattr(marginalVals, "__len__"): marginalVals = [marginalVals]
mVals = list(marginalVals)
try:
rDim = mVals.index(fp)
if rDim < len(mVals) - 1 and fp in mVals[rDim + 1 :]:
raise
except:
raise RROMPyException(("Exactly 1 'freepar' entry in "
"marginalVals must be provided."))
ARBs = self.data.ARBs
if self.data.npar > 1:
mVals[rDim] = self.data.mu0(rDim)
mVals = (checkParameter(mVals) ** self.data.rescalingExp
- self.data.mu0 ** self.data.rescalingExp)
mVals = mVals.data.flatten()
mVals[rDim] = fp
ARBs = marginalizePolyList(ARBs, mVals, "auto")
ev = eigvalsNonlinearDense(ARBs, jSupp = jSupp, **kwargs)
return np.power(self.data.mu0(rDim) ** self.data.rescalingExp[rDim]
+ ev * self.data.scaleFactor[rDim],
1. / self.data.rescalingExp[rDim])
diff --git a/rrompy/utilities/poly_fitting/heaviside/heaviside_interpolator.py b/rrompy/utilities/poly_fitting/heaviside/heaviside_interpolator.py
index 79334ef..f0aca96 100644
--- a/rrompy/utilities/poly_fitting/heaviside/heaviside_interpolator.py
+++ b/rrompy/utilities/poly_fitting/heaviside/heaviside_interpolator.py
@@ -1,78 +1,77 @@
# Copyright (C) 2018 by the RROMPy authors
#
# This file is part of RROMPy.
#
# RROMPy is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# RROMPy is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with RROMPy. If not, see .
#
import numpy as np
from copy import deepcopy as copy
from rrompy.utilities.base.types import (List, ListAny, Np1D, paramList,
interpEng)
from rrompy.utilities.base import freepar as fp
from rrompy.utilities.poly_fitting.polynomial.polynomial_interpolator import (
PolynomialInterpolator)
from rrompy.utilities.poly_fitting.polynomial.roots import polyroots
from rrompy.utilities.poly_fitting.heaviside.val import polyval
from rrompy.utilities.poly_fitting.heaviside.heaviside_to_from_rational import\
heaviside2rational, rational2heaviside
from rrompy.utilities.exception_manager import RROMPyAssert
__all__ = ['HeavisideInterpolator']
class HeavisideInterpolator(PolynomialInterpolator):
"""HERE"""
def __init__(self, other = None):
if other is None: return
self.poles = other.poles
super().__init__(other)
def __call__(self, mu:paramList, der : List[int] = None,
scl : Np1D = None):
return polyval(mu, self.coeffs, self.poles, self.polybasis)
def __copy__(self):
return HeavisideInterpolator(self)
def __deepcopy__(self, memo):
other = HeavisideInterpolator()
other.poles, other.coeffs, other.npar, other.polybasis = copy(
(self.poles, self.coeffs, self.npar, self.polybasis), memo)
return other
def pad(self, nleft : List[int] = None, nright : List[int] = None):
if nleft is None: nleft = [0] * len(self.shape)
if nright is None: nright = [0] * len(self.shape)
if not hasattr(nleft, "__len__"): nleft = [nleft]
if not hasattr(nright, "__len__"): nright = [nright]
RROMPyAssert(nleft[0], 0, "Padding in free direction")
super().pad(nleft, nright)
def setupFromRational(self, num:interpEng, den:interpEng,
murange : Np1D = np.array([-1., 1.]),
- scl : Np1D = None, scaling : List[callable] = None,
- scalingInv : List[callable] = None):
+ scl : Np1D = None, scalingExp : List[float] = None):
self.coeffs, self.poles, self.polybasis = rational2heaviside(num, den,
- murange, scl, scaling, scalingInv)
+ murange, scl,
+ scalingExp)
def roots(self, marginalVals : ListAny = [fp], murange : Np1D = None,
- scaling : List[callable] = None,
- scalingInv : List[callable] = None):
+ scalingExp : List[float] = None):
RROMPyAssert(self.shape, (1,), "Shape of output")
RROMPyAssert(marginalVals, [fp], "Marginal values")
basisN = self.polybasis.split("_")[0]
coeffsN = heaviside2rational(self.coeffs, self.poles, murange, basisN,
- scaling = scaling, scalingInv = scalingInv)[0]
+ scalingExp = scalingExp)[0]
return polyroots(coeffsN, basisN)
diff --git a/rrompy/utilities/poly_fitting/heaviside/heaviside_to_from_rational.py b/rrompy/utilities/poly_fitting/heaviside/heaviside_to_from_rational.py
index 697ac43..89f7e3d 100644
--- a/rrompy/utilities/poly_fitting/heaviside/heaviside_to_from_rational.py
+++ b/rrompy/utilities/poly_fitting/heaviside/heaviside_to_from_rational.py
@@ -1,99 +1,97 @@
# Copyright (C) 2018 by the RROMPy authors
#
# This file is part of RROMPy.
#
# RROMPy is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# RROMPy is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with RROMPy. If not, see .
#
import numpy as np
from rrompy.utilities.poly_fitting.polynomial.polynomial_interpolator import \
PolynomialInterpolator
from rrompy.utilities.poly_fitting.polynomial.vander import polyvander
from rrompy.utilities.poly_fitting.custom_fit import customFit
from rrompy.utilities.base.types import Np1D, Np2D, List, Tuple, interpEng
from rrompy.parameter.parameter_sampling import (RandomSampler as RS,
QuadratureSampler as QS)
from .val import polyval
from rrompy.utilities.exception_manager import RROMPyException, RROMPyAssert
__all__ = ['heaviside2rational', 'rational2heaviside']
def heaviside2rational(c:Np2D, poles:Np1D, murange : Np1D = None,
basis : str = "MONOMIAL_HEAVISIDE", basisD : str = None,
- scaling : List[callable] = None,
- scalingInv : List[callable] = None) \
+ scalingExp : List[float] = None) \
-> Tuple[interpEng, interpEng]:
num = PolynomialInterpolator()
den = PolynomialInterpolator()
basisN = basis.split("_")[0]
if basisD is None: basisD = basisN
if murange is None:
multiplier = [1., 1.j]
avgs = [np.mean(np.real(poles)), np.mean(np.imag(poles))]
ranges = np.array([[np.min(np.real(poles)), np.max(np.real(poles))],
[np.min(np.imag(poles)), np.max(np.imag(poles))]])
domIdx = np.argmax(ranges[:, 1] - ranges[:, 0])
murange = [multiplier[domIdx] * x
+ multiplier[1 - domIdx] * avgs[1 - domIdx]
for x in ranges[domIdx, :]]
extraPt = None
while extraPt is None or np.any(np.isclose(extraPt, poles)):
extraPt = murange[0] + (murange[1] - murange[0]) * np.random.rand(1)
denAuxPts = np.concatenate((poles, extraPt))
denAuxVals = np.concatenate((np.zeros(len(poles)), [1.]))
den.setupByInterpolation(denAuxPts, denAuxVals, len(poles), basisD)
den.coeffs /= np.linalg.norm(den.coeffs)
if basis == "CHEBYSHEV":
- sampler = QS(murange, "CHEBYSHEV", scaling, scalingInv)
+ sampler = QS(murange, "CHEBYSHEV", scalingExp)
elif basis == "LEGENDRE":
- sampler = QS(murange, "GAUSSLEGENDRE", scaling, scalingInv)
+ sampler = QS(murange, "GAUSSLEGENDRE", scalingExp)
else:
- sampler = RS(murange, "HALTON", scaling, scalingInv)
+ sampler = RS(murange, "HALTON", scalingExp)
xAux = sampler.generatePoints(len(c))
valsAux = den(xAux) * polyval(xAux, c, poles, basis)
num.setupByInterpolation(xAux, valsAux, len(c) - 1, basisN)
return num, den
def rational2heaviside(num:interpEng, den:interpEng,
murange : Np1D = np.array([-1., 1.]), scl : Np1D = None,
- scaling : List[callable] = None,
- scalingInv : List[callable] = None) \
+ scalingExp : List[float] = None) \
-> Tuple[Np2D, Np1D, str]:
if (not isinstance(num, PolynomialInterpolator)
or not isinstance(den, PolynomialInterpolator)):
raise RROMPyException(("Rational numerator and denominator must be "
"polynomial interpolators."))
RROMPyAssert(num.npar, 1, "Number of parameters")
RROMPyAssert(den.npar, 1, "Number of parameters")
basis = num.polybasis + "_HEAVISIDE"
c = np.zeros_like(num.coeffs)
poles = den.roots()
Psp = num(poles)
Qspder = den(poles, 1, scl)
c[: len(poles)] = (Psp / Qspder).T
if len(c) > len(poles):
from rrompy.parameter.parameter_sampling import (RandomSampler as RS,
QuadratureSampler as QS)
if num.polybasis == "CHEBYSHEV":
- sampler = QS(murange, "CHEBYSHEV", scaling, scalingInv)
+ sampler = QS(murange, "CHEBYSHEV", scalingExp)
elif num.polybasis == "LEGENDRE":
- sampler = QS(murange, "GAUSSLEGENDRE", scaling, scalingInv)
+ sampler = QS(murange, "GAUSSLEGENDRE", scalingExp)
else:
- sampler = RS(murange, "HALTON", scaling, scalingInv)
+ sampler = RS(murange, "HALTON", scalingExp)
xAux = sampler.generatePoints(len(c))
valsAux = (num(xAux) / den(xAux)
- polyval(xAux, c, poles, basis)).T
VanAux = polyvander(xAux, [len(c) - len(poles) - 1], num.polybasis)
c[len(poles) :] = customFit(VanAux, valsAux)
return c, poles, basis
diff --git a/setup.py b/setup.py
index 32b1d45..9a05576 100644
--- a/setup.py
+++ b/setup.py
@@ -1,52 +1,52 @@
-# Copyright (C) 2015-2018 by the RBniCS authors
+# Copyright (C) 2018 by the RROMPy authors
#
-# This file is part of RBniCS.
+# This file is part of RROMPy.
#
-# RBniCS is free software: you can redistribute it and/or modify
+# 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.
#
-# RBniCS is distributed in the hope that it will be useful,
+# 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 RBniCS. If not, see .
+# along with RROMPy. If not, see .
#
import os
from setuptools import find_packages, setup
rrompy_directory = os.path.abspath(os.path.dirname(os.path.realpath(__file__)))
#rrompy_directory = os.path.join(rrompy_directory, 'rrompy')
setup(name="RROMPy",
description="Rational reduced order modelling in Python",
long_description="Rational reduced order modelling in Python",
author="Davide Pradovera",
author_email="davide.pradovera@epfl.ch",
version="1.6",
license="GNU Library or Lesser General Public License (LGPL)",
classifiers=[
"Development Status :: 1 - Planning"
"Intended Audience :: Developers",
"Intended Audience :: Science/Research",
"Programming Language :: Python :: 3",
"Programming Language :: Python :: 3.4",
"Programming Language :: Python :: 3.5",
"Programming Language :: Python :: 3.6",
"License :: OSI Approved :: GNU Library or Lesser General Public License (LGPL)",
"Topic :: Scientific/Engineering :: Mathematics",
"Topic :: Software Development :: Libraries :: Python Modules",
],
packages=find_packages(rrompy_directory),
setup_requires=[
"pytest-runner"
],
tests_require=[
"pytest"
],
zip_safe=False
)
diff --git a/tests/reduction_methods_1D/rb_distributed_1d.py b/tests/reduction_methods_1D/reduced_basis_distributed_1d.py
similarity index 89%
rename from tests/reduction_methods_1D/rb_distributed_1d.py
rename to tests/reduction_methods_1D/reduced_basis_distributed_1d.py
index 3b1f288..6905e65 100644
--- a/tests/reduction_methods_1D/rb_distributed_1d.py
+++ b/tests/reduction_methods_1D/reduced_basis_distributed_1d.py
@@ -1,57 +1,57 @@
# Copyright (C) 2018 by the RROMPy authors
#
# This file is part of RROMPy.
#
# RROMPy is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# RROMPy is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with RROMPy. If not, see .
#
import numpy as np
from matrix_fft import matrixFFT
-from rrompy.reduction_methods.distributed import RBDistributed as RBD
+from rrompy.reduction_methods.distributed import ReducedBasis as RB
from rrompy.parameter.parameter_sampling import (QuadratureSampler as QS,
ManualSampler as MS)
from rrompy.parameter import checkParameterList
def test_LS():
solver = matrixFFT()
params = {"POD": True, "R": 5, "S": 10,
"sampler": QS([1., 7.], "CHEBYSHEV")}
- approx = RBD(solver, 4., params, verbosity = 0)
+ approx = RB(solver, 4., params, verbosity = 0)
approx.setupApprox()
for mu in approx.mus:
assert not np.isclose(approx.normErr(mu)[0], 0., atol = 1e-7)
approx.POD = False
approx.setupApprox()
for mu in approx.mus[approx.R :]:
assert not np.isclose(approx.normErr(mu)[0], 0., atol = 1e-3)
def test_interp():
solver = matrixFFT()
params = {"POD": False, "S": 10, "sampler": QS([1., 7.], "CHEBYSHEV")}
- approx = RBD(solver, 4., params, verbosity = 0)
+ approx = RB(solver, 4., params, verbosity = 0)
approx.setupApprox()
for mu in approx.mus:
assert np.isclose(approx.normErr(mu)[0], 0., atol = 1e-7)
def test_hermite():
mu = 1.5
solver = matrixFFT()
sampler0 = QS([1., 7.], "CHEBYSHEV")
points, _ = checkParameterList(np.tile(sampler0.generatePoints(4)(0), 3))
params = {"POD": True, "S": 12, "sampler": MS([1., 7.], points = points)}
- approx = RBD(solver, 4., params, verbosity = 0)
+ approx = RB(solver, 4., params, verbosity = 0)
approx.setupApprox()
for mu in approx.mus:
assert np.isclose(approx.normErr(mu)[0], 0., atol = 1e-8)
diff --git a/tests/reduction_methods_1D/rb_distributed_greedy_1d.py b/tests/reduction_methods_1D/reduced_basis_distributed_greedy_1d.py
similarity index 91%
rename from tests/reduction_methods_1D/rb_distributed_greedy_1d.py
rename to tests/reduction_methods_1D/reduced_basis_distributed_greedy_1d.py
index 512a049..e0a2462 100644
--- a/tests/reduction_methods_1D/rb_distributed_greedy_1d.py
+++ b/tests/reduction_methods_1D/reduced_basis_distributed_greedy_1d.py
@@ -1,60 +1,60 @@
# Copyright (C) 2018 by the RROMPy authors
#
# This file is part of RROMPy.
#
# RROMPy is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# RROMPy is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with RROMPy. If not, see .
#
import numpy as np
from matrix_fft import matrixFFT
-from rrompy.reduction_methods.distributed_greedy import RBDistributedGreedy \
- as RBDG
+from rrompy.reduction_methods.distributed_greedy import ReducedBasisGreedy \
+ as RBG
from rrompy.parameter.parameter_sampling import QuadratureSampler as QS
def test_lax_tolerance(capsys):
mu = 2.25
solver = matrixFFT()
params = {"POD": True, "sampler": QS([1.5, 6.5], "UNIFORM"),
"S": 4, "greedyTol": 1e-2,
"trainSetGenerator": QS([1.5, 6.5], "CHEBYSHEV")}
- approx = RBDG(solver, 4., params, verbosity = 10)
+ approx = RBG(solver, 4., params, verbosity = 10)
approx.greedy()
out, err = capsys.readouterr()
assert "Done computing snapshots (final snapshot count: 10)." in out
assert len(err) == 0
assert len(approx.mus) == 10
_, _, maxEst = approx.getMaxErrorEstimator(approx.muTest)
assert maxEst < 1e-2
assert np.isclose(approx.normErr(mu)[0], 1.5056e-05, rtol = 1e-1)
def test_samples_at_poles():
solver = matrixFFT()
params = {"POD": True, "sampler": QS([1.5, 6.5], "UNIFORM"),
"S": 4, "nTestPoints": 100, "greedyTol": 1e-5,
"trainSetGenerator": QS([1.5, 6.5], "CHEBYSHEV")}
- approx = RBDG(solver, 4., params, verbosity = 0)
+ approx = RBG(solver, 4., params, verbosity = 0)
approx.greedy()
for mu in approx.mus:
assert np.isclose(approx.normErr(mu)[0] / (1e-15+approx.normHF(mu)[0]),
0., atol = 1e-4)
poles = approx.getPoles()
for lambda_ in range(2, 7):
assert np.isclose(np.min(np.abs(poles - lambda_)), 0., atol = 1e-3)
assert np.isclose(np.min(np.abs(np.array(approx.mus(0)) - lambda_)),
0., atol = 1e-1)
diff --git a/tests/reduction_methods_multiD/rb_distributed_2d.py b/tests/reduction_methods_multiD/reduced_basis_distributed_2d.py
similarity index 90%
rename from tests/reduction_methods_multiD/rb_distributed_2d.py
rename to tests/reduction_methods_multiD/reduced_basis_distributed_2d.py
index 8182e2d..2008338 100644
--- a/tests/reduction_methods_multiD/rb_distributed_2d.py
+++ b/tests/reduction_methods_multiD/reduced_basis_distributed_2d.py
@@ -1,62 +1,62 @@
# Copyright (C) 2018 by the RROMPy authors
#
# This file is part of RROMPy.
#
# RROMPy is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# RROMPy is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with RROMPy. If not, see .
#
import numpy as np
from matrix_random import matrixRandom
-from rrompy.reduction_methods.distributed import RBDistributed as RBD
+from rrompy.reduction_methods.distributed import ReducedBasis as RB
from rrompy.parameter.parameter_sampling import (QuadratureSampler as QS,
ManualSampler as MS)
from rrompy.parameter import checkParameterList
def test_LS():
mu0 = [2, 3]
solver = matrixRandom()
params = {"POD": True, "R": 5, "S": [3, 3],
"sampler": QS([[0., 4.], [1., 5.]], "CHEBYSHEV")}
- approx = RBD(solver, mu0, params, verbosity = 0)
+ approx = RB(solver, mu0, params, verbosity = 0)
approx.setupApprox()
for mu in approx.mus:
assert not np.isclose(approx.normErr(mu)[0], 0., atol = 1e-7)
approx.POD = False
approx.setupApprox()
for mu in approx.mus[approx.R :]:
assert not np.isclose(approx.normErr(mu)[0], 0., atol = 1e-3)
def test_interp():
mu0 = [2, 3]
solver = matrixRandom()
params = {"POD": False, "S": [3, 3],
"sampler": QS([[0., 4.], [1., 5.]], "CHEBYSHEV")}
- approx = RBD(solver, mu0, params, verbosity = 0)
+ approx = RB(solver, mu0, params, verbosity = 0)
approx.setupApprox()
for mu in approx.mus:
assert np.isclose(approx.normErr(mu)[0], 0., atol = 1e-7)
def test_hermite():
mu0 = [2, 3]
solver = matrixRandom()
sampler0 = QS([[0., 4.], [1., 5.]], "CHEBYSHEV")
points, _ = checkParameterList(np.tile(
sampler0.generatePoints([2, 2]).data, [3, 1]))
params = {"POD": True, "S": [12],
"sampler": MS([[0., 4.], [1., 5.]], points = points)}
- approx = RBD(solver, mu0, params, verbosity = 0)
+ approx = RB(solver, mu0, params, verbosity = 0)
approx.setupApprox()
for mu in approx.mus:
assert np.isclose(approx.normErr(mu)[0], 0., atol = 1e-8)
diff --git a/tests/utilities/parameter_sampling.py b/tests/utilities/parameter_sampling.py
index b94792e..b1ea088 100644
--- a/tests/utilities/parameter_sampling.py
+++ b/tests/utilities/parameter_sampling.py
@@ -1,60 +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 .
#
import numpy as np
from rrompy.parameter.parameter_sampling import (ManualSampler,
QuadratureSampler, RandomSampler, FFTSampler)
from rrompy.parameter import checkParameter
def test_manual():
sampler = ManualSampler(lims = [0., 3.], points = np.linspace(0, 3, 101),
- scaling = lambda x: np.power(x, 2.),
- scalingInv = lambda x: np.power(x, .5))
+ scalingExp = 2.)
assert sampler.name() == "ManualSampler"
x = sampler.generatePoints(10)
assert np.allclose(x(0), np.linspace(0, 3, 101)[:10], rtol = 1e-5)
def test_quadrature():
sampler = QuadratureSampler(lims = [0., 3.], kind = "CHEBYSHEV")
x = sampler.generatePoints(9)
assert np.isclose(x(0)[2], 1.5, rtol = 1e-5)
def test_random():
sampler = RandomSampler(lims = [0., 3.], kind = "SOBOL", seed = 13432)
x = sampler.generatePoints(100)
assert np.isclose(x(0)[47], 0.55609130859375, rtol = 1e-5)
def test_fft():
sampler = FFTSampler(lims = [-1., 1.])
x = sampler.generatePoints(100)
assert np.allclose(np.power(x(0), 100), 1., rtol = 1e-5)
def test_2D():
sampler = QuadratureSampler(lims = [(0., 0.), (3., 1.)],
kind = "GAUSSLEGENDRE")
x = sampler.generatePoints((9, 5))
assert sum(np.isclose(x(0), 1.5)) == 5
assert sum(np.isclose(x(1), .5)) == 9
def test_4D():
sampler = RandomSampler(lims = [tuple([0.] * 4), tuple([1.] * 4)],
kind = "UNIFORM", seed = 1234)
x = sampler.generatePoints(10)
assert x.shape[1] == 4
assert checkParameter([x[0]]) == checkParameter([(0.191519450378892,
0.622108771039832, 0.437727739007115, 0.785358583713769)])