diff --git a/rrompy/hfengines/base/hfengine_base.py b/rrompy/hfengines/base/hfengine_base.py
index 79076f0..1fbbeab 100644
--- a/rrompy/hfengines/base/hfengine_base.py
+++ b/rrompy/hfengines/base/hfengine_base.py
@@ -1,397 +1,397 @@
# Copyright (C) 2018-2020 by the RROMPy authors
#
# This file is part of RROMPy.
#
# RROMPy is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# RROMPy is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with RROMPy. If not, see .
#
from abc import abstractmethod
import numpy as np
import scipy.sparse as scsp
from numbers import Number
from collections.abc import Iterable
from copy import copy as softcopy
from rrompy.utilities.base.decorators import (nonaffine_construct,
mu_independent)
from rrompy.utilities.base.types import (Np1D, Np2D, Tuple, List, DictAny,
paramVal, paramList, sampList)
from rrompy.utilities.numerical import solve as tsolve, dot, potential
from rrompy.utilities.expression import expressionEvaluator
from rrompy.utilities.exception_manager import RROMPyAssert
from rrompy.sampling.sample_list import sampleList
from rrompy.parameter import (checkParameter, checkParameterList,
parameterList, parameterMap as pMap)
from rrompy.solver.linear_solver import setupSolver
from rrompy.utilities.parallel import (poolRank, masterCore, listScatter,
matrixGatherv, isend, recv)
__all__ = ['HFEngineBase']
class HFEngineBase:
"""Generic solver for parametric problems."""
def __init__(self, verbosity : int = 10, timestamp : bool = True):
self.verbosity = verbosity
self.timestamp = timestamp
self.setSolver("SPSOLVE", {"use_umfpack" : False})
self.npar = 0
self._C = None
def name(self) -> str:
return self.__class__.__name__
def __str__(self) -> str:
return self.name()
def __repr__(self) -> str:
return self.__str__() + " at " + hex(id(self))
def __dir_base__(self):
return [x for x in self.__dir__() if x[:2] != "__"]
def __deepcopy__(self, memo):
return softcopy(self)
@property
def npar(self):
"""Value of npar."""
return self._npar
@npar.setter
def npar(self, npar):
nparOld = self._npar if hasattr(self, "_npar") else -1
if npar != nparOld:
self.parameterMap = pMap(1., npar)
self._npar = npar
@property
def spacedim(self):
return 1
@property
def outputdim(self):
if self.isCEye: return self.spacedim
return self._C.shape[0]
def checkParameter(self, mu:paramVal) -> paramVal:
muP = checkParameter(mu, self.npar)
if self.npar == 0: muP.reset((1, 0), np.complex)
return muP
def checkParameterList(self, mu:paramList,
check_if_single : bool = False) -> paramList:
muL = checkParameterList(mu, self.npar, check_if_single)
return muL
def mapParameterList(self, mu:paramList, direct : str = "F",
idx : List[int] = None) -> paramList:
if idx is None: idx = np.arange(self.npar)
muMapped = checkParameterList(mu, len(idx))
for j, d in enumerate(idx):
muMapped.data[:, j] = expressionEvaluator(
self.parameterMap[direct][d],
muMapped(j)).flatten()
return muMapped
@property
def energyNormMatrix(self):
if not hasattr(self, "_energyNormMatrix"):
self.buildEnergyNormForm()
return self._energyNormMatrix
def buildEnergyNormForm(self):
"""
Build sparse matrix (in CSR format) representative of scalar product.
"""
self._energyNormMatrix = 1.
@property
def energyNormDualMatrix(self):
if not hasattr(self, "_energyNormDualMatrix"):
self.buildEnergyNormDualForm()
return self._energyNormDualMatrix
def buildEnergyNormDualForm(self):
"""
Build sparse matrix (in CSR format) representative of dual scalar
product without duality.
"""
self._energyNormDualMatrix = 1.
@property
def energyNormOutputMatrix(self):
if not hasattr(self, "_energyNormOutputMatrix"):
self.buildEnergyNormOutput()
return self._energyNormOutputMatrix
def buildEnergyNormOutput(self):
"""
Build sparse matrix (in CSR format) representative of scalar product
over output space.
"""
self._energyNormOutputMatrix = 1.
def innerProduct(self, u:Np2D, v:Np2D, onlyDiag : bool = False,
dual : bool = False, is_state : bool = True) -> Np2D:
"""Scalar product."""
if is_state or self.isCEye:
if dual:
energyMat = self.energyNormDualMatrix
else:
energyMat = self.energyNormMatrix
else:
energyMat = self.energyNormOutputMatrix
if isinstance(u, (parameterList, sampleList)): u = u.data
if isinstance(v, (parameterList, sampleList)): v = v.data
if onlyDiag:
return np.sum(dot(energyMat, u) * v.conj(), axis = 0)
return dot(dot(energyMat, u).T, v.conj()).T
def norm(self, u:Np2D, dual : bool = False,
is_state : bool = True) -> Np1D:
return np.abs(self.innerProduct(u, u, onlyDiag = True, dual = dual,
is_state = is_state)) ** .5
def baselineA(self):
"""Return 0 of shape consistent with operator of linear system."""
if (hasattr(self, "As") and isinstance(self.As, Iterable)
and self.As[0] is not None):
d = self.As[0].shape[0]
else:
d = self.spacedim
return scsp.csr_matrix((np.zeros(0), np.zeros(0), np.zeros(d + 1)),
shape = (d, d), dtype = np.complex)
def baselineb(self):
"""Return 0 of shape consistent with RHS of linear system."""
return np.zeros(self.spacedim, dtype = np.complex)
@nonaffine_construct
@abstractmethod
def A(self, mu : paramVal = [], der : List[int] = 0) -> Np2D:
"""
Assemble terms of operator of linear system and return it (or its
derivative) at a given parameter.
"""
return
@nonaffine_construct
@abstractmethod
def b(self, mu : paramVal = [], der : List[int] = 0) -> Np1D:
"""
Assemble terms of RHS of linear system and return it (or its
derivative) at a given parameter.
"""
return
@mu_independent
def C(self, mu:paramVal):
"""
Value of C. Should be overridden (with something like
return self._C(mu)
) if a mu-dependent C is needed.
"""
if self._C is None: self._C = 1.
return self._C
@property
def isCEye(self):
"""
Whether the action of C can be seen as a scalar multiplication. Should
be overridden (with
return True
) if a mu-dependent scalar C is used.
"""
return isinstance(self._C, Number)
def applyC(self, u:sampList, mu:paramVal):
"""Apply LHS of linear system."""
return dot(self.C(mu), u)
def setSolver(self, solverType:str, solverArgs : DictAny = {}):
"""Choose solver type and parameters."""
self._solver, self._solverArgs = setupSolver(solverType, solverArgs)
def solve(self, mu : paramList = [], RHS : sampList = None,
return_state : bool = False) -> sampList:
"""
Find solution of linear system.
Args:
mu: parameter value.
RHS: RHS of linear system. If None, defaults to that of parametric
system. Defaults to None.
return_state: whether to return state before multiplication by c.
Defaults to False.
"""
from rrompy.sampling import sampleList, emptySampleList
if mu == []: mu = self.mu0
mu = self.checkParameterList(mu)
if len(mu) == 0: return emptySampleList()
+ req, is_master = [], masterCore()
mu_loc, idx, sizes = listScatter(mu, return_sizes = True)
mu_loc = self.checkParameterList(mu_loc)
- req, emptyCores = [], np.where(sizes == 0)[0]
+ emptyCores = np.where(np.array(sizes) == 0)[0] if is_master else []
if len(mu_loc) == 0:
uL = recv(source = 0, tag = poolRank())
sol = np.empty((uL, 0), dtype = np.complex)
else:
if RHS is None: # build RHSs
RHS = sampleList([self.b(m) for m in mu_loc])
else:
RHS = sampleList(RHS)
if len(RHS) > 1: RHS = sampleList([RHS[i] for i in idx])
mult = 0 if len(RHS) == 1 else 1
RROMPyAssert(mult * (len(mu_loc) - 1) + 1, len(RHS), "Sample size")
for j, mj in enumerate(mu_loc):
u = tsolve(self.A(mj), RHS[mult * j], self._solver,
self._solverArgs)
if not return_state: u = self.applyC(u, mj)
if j == 0:
sol = np.empty((len(u), len(mu_loc)), dtype = np.complex)
- if masterCore():
- for dest in emptyCores:
- req += [isend(len(u), dest = dest, tag = dest)]
+ for dest in emptyCores:
+ req += [isend(len(u), dest = dest, tag = dest)]
sol[:, j] = u
for r in req: r.wait()
sol = matrixGatherv(sol, sizes)
return sampleList(sol)
def residual(self, mu : paramList = [], u : sampList = None,
post_c : bool = True) -> sampList:
"""
Find residual of linear system for given approximate solution.
Args:
mu: parameter value.
u: numpy complex array with function dofs. If None, set to 0.
post_c: whether to post-process using c. Defaults to True.
"""
from rrompy.sampling import sampleList, emptySampleList
if mu == []: mu = self.mu0
mu = self.checkParameterList(mu)
if len(mu) == 0: return emptySampleList()
+ req, is_master = [], masterCore()
mu_loc, idx, sizes = listScatter(mu, return_sizes = True)
mu_loc = self.checkParameterList(mu_loc)
- req, emptyCores = [], np.where(sizes == 0)[0]
+ emptyCores = np.where(np.array(sizes) == 0)[0] if is_master else []
if len(mu_loc) == 0:
uL = recv(source = 0, tag = poolRank())
res = np.empty((uL, 0), dtype = np.complex)
else:
v = sampleList(np.zeros((self.spacedim, len(mu_loc))))
if u is not None:
u = sampleList(u)
v = v + sampleList([u[i] for i in idx])
for j, (mj, vj) in enumerate(zip(mu_loc, v)):
r = self.b(mj) - dot(self.A(mj), vj)
if post_c: r = self.applyC(r, mj)
if j == 0:
res = np.empty((len(r), len(mu_loc)), dtype = np.complex)
- if masterCore():
- for dest in emptyCores:
- req += [isend(len(r), dest = dest, tag = dest)]
+ for dest in emptyCores:
+ req += [isend(len(r), dest = dest, tag = dest)]
res[:, j] = r
for r in req: r.wait()
res = matrixGatherv(res, sizes)
return sampleList(res)
cutOffPolesRMax,cutOffPolesRMin = np.inf, - np.inf
cutOffPolesIMax, cutOffPolesIMin = np.inf, - np.inf
def flagBadPolesResiduesAbsolute(self, poles:Np1D, residues : Np1D = None,
projMat : Np2D = None) -> Np1D:
"""
Flag (numerical) poles/residues which are impossible.
Args:
poles: poles to be judged.
residues: residues norms to be judged.
projMat: matrix for projection of residues.
"""
poles = np.array(poles).flatten()
flag = np.zeros(len(poles), dtype = bool)
RMax, RMin = self.cutOffPolesRMax, self.cutOffPolesRMin
IMax, IMin = self.cutOffPolesIMax, self.cutOffPolesIMin
if not np.isinf(RMax): flag = flag + (np.real(poles) > RMax)
if not np.isinf(RMin): flag = flag + (np.real(poles) < RMin)
if not np.isinf(IMax): flag = flag + (np.imag(poles) > IMax)
if not np.isinf(IMin): flag = flag + (np.imag(poles) < IMin)
return flag
cutOffPolesPotentialMax = np.inf
cutOffPolesRMaxRel, cutOffPolesRMinRel = np.inf, - np.inf
cutOffPolesIMaxRel, cutOffPolesIMinRel = np.inf, - np.inf
cutOffResNormMin = -1
cutOffResAngleMin, cutOffResAngleMax = -1, np.pi + 1
def flagBadPolesResiduesRelative(self, poles:Np1D, residues : Np1D = None,
projMat : Np2D = None,
foci : Tuple[float, float] = [-1., 1.]) \
-> Np1D:
"""
Flag (numerical) poles/residues which are impossible.
Args:
poles: poles to be judged.
residues: residues norms to be judged.
projMat: matrix for projection of residues.
foci: foci for potential evaluation.
"""
poles = np.array(poles).flatten()
flag = np.zeros(len(poles), dtype = bool)
potMax = self.cutOffPolesPotentialMax
RMax, RMin = self.cutOffPolesRMaxRel, self.cutOffPolesRMinRel
IMax, IMin = self.cutOffPolesIMaxRel, self.cutOffPolesIMinRel
if not np.isinf(potMax) or (residues is not None
and not self._ignoreResidues):
plsInf = np.isinf(poles)
pot = potential(poles, foci)
if not np.isinf(potMax): flag = flag + (pot > potMax)
if not np.isinf(RMax): flag = flag + (np.real(poles) > RMax)
if not np.isinf(RMin): flag = flag + (np.real(poles) < RMin)
if not np.isinf(IMax): flag = flag + (np.imag(poles) > IMax)
if not np.isinf(IMin): flag = flag + (np.imag(poles) < IMin)
if residues is not None and not self._ignoreResidues:
residues = np.array(residues).reshape(-1, len(poles))
resGood = np.where(flag + plsInf == False)[0]
if len(resGood) > 0:
residues = residues[:, resGood] / pot[resGood]
if projMat is None:
resNorm = np.linalg.norm(residues, axis = 0)
else:
residues = projMat.dot(residues)
resNorm = self.norm(residues)
if self.cutOffResNormMin > 0.:
flag[resGood[resNorm < self.cutOffResNormMin
* np.max(resNorm)]] = 1
resGood = np.where(flag + plsInf == False)[0]
if len(resGood) > 0 and (self.cutOffResAngleMin > 0.
or self.cutOffResAngleMax < np.pi):
if projMat is None:
angles = np.real(residues.T.conj().dot(residues))
else:
angles = np.real(self.innerProduct(residues, residues))
resNormEff = resNorm
resNormEff[np.isclose(resNormEff, 0., atol = 1e-15)] = 1.
angles = np.clip((angles / resNormEff).T / resNormEff, -1., 1.)
angles = np.arccos(angles)
badangles = ((angles < self.cutOffResAngleMin)
+ (angles > self.cutOffResAngleMax))
badangles[np.arange(len(angles)), np.arange(len(angles))] = 0
idx = np.zeros(len(angles), dtype = bool)
while np.sum(badangles) > 0:
idxn = np.argmax(np.sum(badangles, axis = 1))
badangles[idxn], badangles[:, idxn] = 0, 0
idx[idxn] = True
flag[resGood[idx]] = 1
return flag > 0
@property
def _ignoreResidues(self):
return (self.cutOffResNormMin <= 0. and self.cutOffResAngleMin <= 0.
and self.cutOffResAngleMax >= np.pi)
diff --git a/rrompy/parameter/parameter_sampling/segment/quadrature_sampler.py b/rrompy/parameter/parameter_sampling/segment/quadrature_sampler.py
index d177c67..d7ea0ab 100644
--- a/rrompy/parameter/parameter_sampling/segment/quadrature_sampler.py
+++ b/rrompy/parameter/parameter_sampling/segment/quadrature_sampler.py
@@ -1,48 +1,48 @@
# Copyright (C) 2018-2020 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.generic_quadrature_sampler import (
GenericQuadratureSampler)
from rrompy.utilities.base.types import paramList
from rrompy.utilities.numerical import (lowDiscrepancy, kroneckerer,
quadraturePointsGenerate)
__all__ = ['QuadratureSampler']
class QuadratureSampler(GenericQuadratureSampler):
"""Generator of quadrature sample points."""
def generatePoints(self, n:int, reorder : bool = True) -> paramList:
"""Array of sample points."""
n1d = int(np.ceil(n ** (1. / self.npar)))
nleft, nright = 1, n1d ** self.npar
xmat = np.empty((nright, self.npar), dtype = np.complex)
limsE = self.mapParameterList(self.lims)
for d in range(self.npar):
nright //= n1d
a, b = limsE(d)
- c, r = (a + b) / 2., (a - b) / 2.
+ c, r = .5 * (a + b), .5 * (b - a)
xd = c + r * quadraturePointsGenerate(n1d, self.kind)
xmat[:, d] = kroneckerer(xd, nleft, nright)
nleft *= n1d
nright = n1d ** self.npar
if nright > 1 and reorder:
fejerOrdering = [nright - 1] + lowDiscrepancy(nright - 1)
xmat = xmat[fejerOrdering, :]
return self.mapParameterList(xmat, "B")
diff --git a/rrompy/parameter/parameter_sampling/shape/fft_sampler.py b/rrompy/parameter/parameter_sampling/shape/fft_sampler.py
index ab3b869..49f5a72 100644
--- a/rrompy/parameter/parameter_sampling/shape/fft_sampler.py
+++ b/rrompy/parameter/parameter_sampling/shape/fft_sampler.py
@@ -1,48 +1,48 @@
# Copyright (C) 2018-2020 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_shape_sampler import GenericShapeSampler
from rrompy.utilities.base.types import paramList
from rrompy.utilities.numerical import lowDiscrepancy, kroneckerer
__all__ = ['FFTSampler']
class FFTSampler(GenericShapeSampler):
"""Generator of FFT-type sample points on scaled roots of unity."""
def generatePoints(self, n:int, reorder : bool = True) -> paramList:
"""Array of sample points."""
n1d = int(np.ceil(n ** (1. / self.npar)))
nleft, nright = 1, n1d ** self.npar
xmat = np.empty((nright, self.npar), dtype = np.complex)
limsE = self.mapParameterList(self.lims)
for d in range(self.npar):
nright //= n1d
a, b = limsE(d)
- c, r = (a + b) / 2., (a - b) / 2.
+ c, r = .5 * (a + b), .5 * (b - a)
unitpts = np.exp(1.j * np.linspace(0, 2 * np.pi, n1d + 1)[:-1])
unitpts = (np.real(unitpts)
+ 1.j * self.axisRatios[d] * np.imag(unitpts))
xd = c + r * unitpts
if n1d > 1 and reorder:
fejerOrdering = [n1d - 1] + lowDiscrepancy(n1d - 1)
xd = xd[fejerOrdering]
xmat[:, d] = kroneckerer(xd, nleft, nright)
nleft *= n1d
return self.mapParameterList(xmat, "B")
diff --git a/rrompy/parameter/parameter_sampling/shape/quadrature_box_sampler.py b/rrompy/parameter/parameter_sampling/shape/quadrature_box_sampler.py
index c1ef160..b9f3836 100644
--- a/rrompy/parameter/parameter_sampling/shape/quadrature_box_sampler.py
+++ b/rrompy/parameter/parameter_sampling/shape/quadrature_box_sampler.py
@@ -1,55 +1,55 @@
# Copyright (C) 2018-2020 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_shape_quadrature_sampler import GenericShapeQuadratureSampler
from rrompy.utilities.base.types import paramList
from rrompy.utilities.numerical import (lowDiscrepancy, kroneckerer,
quadraturePointsGenerate)
__all__ = ['QuadratureBoxSampler']
class QuadratureBoxSampler(GenericShapeQuadratureSampler):
"""Generator of quadrature sample points on boxes."""
def generatePoints(self, n:int, reorder : bool = True) -> paramList:
"""Array of sample points."""
n1d = int(np.ceil(n ** (1. / self.npar)))
xds, nds = [None] * self.npar, [None] * self.npar
limsE = self.mapParameterList(self.lims)
for d in range(self.npar):
ax = self.axisRatios[d]
a, b = limsE(d)
- c, r = (a + b) / 2., (a - b) / 2.
+ c, r = .5 * (a + b), .5 * (b - a)
n1dx = int(np.ceil((n1d / ax) ** .5))
n1dy = int(np.ceil(n1d / n1dx))
Xdx, Xdy = np.meshgrid(quadraturePointsGenerate(n1dx, self.kind),
quadraturePointsGenerate(n1dy, self.kind))
Z = Xdx.flatten() + 1.j * ax * Xdy.flatten()
xds[d] = c + r * Z
nds[d] = len(Z)
nleft, nright = 1, np.prod(nds)
xmat = np.empty((nright, self.npar), dtype = np.complex)
for d in range(self.npar):
nright //= nds[d]
xmat[:, d] = kroneckerer(xds[d], nleft, nright)
nleft *= nds[d]
if len(xmat) > 1 and reorder:
fejerOrdering = [len(xmat) - 1] + lowDiscrepancy(len(xmat) - 1)
xmat = xmat[fejerOrdering, :]
return self.mapParameterList(xmat, "B")
diff --git a/rrompy/parameter/parameter_sampling/shape/quadrature_circle_sampler.py b/rrompy/parameter/parameter_sampling/shape/quadrature_circle_sampler.py
index 235925a..13ebaa9 100644
--- a/rrompy/parameter/parameter_sampling/shape/quadrature_circle_sampler.py
+++ b/rrompy/parameter/parameter_sampling/shape/quadrature_circle_sampler.py
@@ -1,61 +1,61 @@
# Copyright (C) 2018-2020 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_shape_quadrature_sampler import GenericShapeQuadratureSampler
from rrompy.utilities.base.types import paramList
from rrompy.utilities.numerical import (lowDiscrepancy, kroneckerer,
potential, quadraturePointsGenerate)
__all__ = ['QuadratureCircleSampler']
class QuadratureCircleSampler(GenericShapeQuadratureSampler):
"""Generator of quadrature sample points on ellipses."""
def generatePoints(self, n:int, reorder : bool = True) -> paramList:
"""Array of sample points."""
n1d = int(np.ceil(n ** (1. / self.npar)))
xds, nds = [None] * self.npar, [None] * self.npar
limsE = self.mapParameterList(self.lims)
for d in range(self.npar):
ax = self.axisRatios[d]
a, b = limsE(d)
- c, r = (a + b) / 2., (a - b) / 2.
+ c, r = .5 * (a + b), .5 * (b - a)
n1dx = int(np.ceil((n1d * 4. / np.pi / ax) ** .5))
n1dy = int(np.ceil(n1d * 4. / np.pi / n1dx))
even, Z = True, []
while len(Z) < n1d:
Xdx, Xdy = np.meshgrid(
quadraturePointsGenerate(n1dx, self.kind),
quadraturePointsGenerate(n1dy, self.kind))
Z = Xdx.flatten() + 1.j * ax * Xdy.flatten()
Z = Z[potential(Z, self.normalFoci(d)) <= 1.]
if even: n1dx += 1
else: n1dy += 1
xds[d] = c + r * Z
nds[d] = len(Z)
nleft, nright = 1, np.prod(nds)
xmat = np.empty((nright, self.npar), dtype = np.complex)
for d in range(self.npar):
nright //= nds[d]
xmat[:, d] = kroneckerer(xds[d], nleft, nright)
nleft *= nds[d]
if len(xmat) > 1 and reorder:
fejerOrdering = [len(xmat) - 1] + lowDiscrepancy(len(xmat) - 1)
xmat = xmat[fejerOrdering, :]
return self.mapParameterList(xmat, "B")
diff --git a/rrompy/reduction_methods/pivoted/greedy/generic_pivoted_greedy_approximant.py b/rrompy/reduction_methods/pivoted/greedy/generic_pivoted_greedy_approximant.py
index a405400..2449b6d 100644
--- a/rrompy/reduction_methods/pivoted/greedy/generic_pivoted_greedy_approximant.py
+++ b/rrompy/reduction_methods/pivoted/greedy/generic_pivoted_greedy_approximant.py
@@ -1,942 +1,942 @@
# Copyright (C) 2018-2020 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 copy import deepcopy as copy
import numpy as np
from collections.abc import Iterable
from matplotlib import pyplot as plt
from rrompy.reduction_methods.pivoted.generic_pivoted_approximant import (
GenericPivotedApproximantBase,
GenericPivotedApproximantPolyMatch,
GenericPivotedApproximantPoleMatch)
from rrompy.reduction_methods.pivoted.gather_pivoted_approximant import (
gatherPivotedApproximant)
from rrompy.utilities.base.types import (Np1D, Np2D, Tuple, List, paramVal,
paramList, ListAny)
from rrompy.utilities.base import verbosityManager as vbMng
from rrompy.utilities.numerical import dot
from rrompy.utilities.numerical.point_matching import pointMatching
from rrompy.utilities.numerical.point_distances import doubleDistanceMatrix
from rrompy.utilities.exception_manager import (RROMPyException, RROMPyAssert,
RROMPyWarning)
from rrompy.parameter import emptyParameterList
from rrompy.utilities.parallel import (masterCore, indicesScatter,
arrayGatherv, isend)
__all__ = ['GenericPivotedGreedyApproximantPolyMatch',
'GenericPivotedGreedyApproximantPoleMatch']
class GenericPivotedGreedyApproximantBase(GenericPivotedApproximantBase):
_allowedEstimatorKindsMarginal = ["LOOK_AHEAD", "LOOK_AHEAD_DOERFLER",
"NONE"]
def __init__(self, *args, **kwargs):
self._preInit()
self._addParametersToList(["matchingWeightError",
"matchingErrorRelative",
"errorEstimatorKindMarginal",
"greedyTolMarginal", "maxIterMarginal",
"autoCollapse"],
[0., False, "NONE", 1e-1, 1e2, False])
super().__init__(*args, **kwargs)
self._postInit()
@property
def scaleFactorDer(self):
"""Value of scaleFactorDer."""
if self._scaleFactorDer == "NONE": return 1.
if self._scaleFactorDer == "AUTO": return self._scaleFactorOldPivot
return self._scaleFactorDer
@scaleFactorDer.setter
def scaleFactorDer(self, scaleFactorDer):
if isinstance(scaleFactorDer, (str,)):
scaleFactorDer = scaleFactorDer.upper()
elif isinstance(scaleFactorDer, Iterable):
scaleFactorDer = list(scaleFactorDer)
self._scaleFactorDer = scaleFactorDer
self._approxParameters["scaleFactorDer"] = self._scaleFactorDer
@property
def samplerMarginal(self):
"""Value of samplerMarginal."""
return self._samplerMarginal
@samplerMarginal.setter
def samplerMarginal(self, samplerMarginal):
if 'refine' not in dir(samplerMarginal):
raise RROMPyException("Marginal sampler type not recognized.")
GenericPivotedApproximantBase.samplerMarginal.fset(self,
samplerMarginal)
@property
def errorEstimatorKindMarginal(self):
"""Value of errorEstimatorKindMarginal."""
return self._errorEstimatorKindMarginal
@errorEstimatorKindMarginal.setter
def errorEstimatorKindMarginal(self, errorEstimatorKindMarginal):
errorEstimatorKindMarginal = errorEstimatorKindMarginal.upper()
if errorEstimatorKindMarginal not in (
self._allowedEstimatorKindsMarginal):
RROMPyWarning(("Marginal error estimator kind not recognized. "
"Overriding to 'NONE'."))
errorEstimatorKindMarginal = "NONE"
self._errorEstimatorKindMarginal = errorEstimatorKindMarginal
self._approxParameters["errorEstimatorKindMarginal"] = (
self.errorEstimatorKindMarginal)
@property
def matchingWeightError(self):
"""Value of matchingWeightError."""
return self._matchingWeightError
@matchingWeightError.setter
def matchingWeightError(self, matchingWeightError):
self._matchingWeightError = matchingWeightError
self._approxParameters["matchingWeightError"] = (
self.matchingWeightError)
@property
def matchingErrorRelative(self):
"""Value of matchingErrorRelative."""
return self._matchingErrorRelative
@matchingErrorRelative.setter
def matchingErrorRelative(self, matchingErrorRelative):
self._matchingErrorRelative = matchingErrorRelative
self._approxParameters["matchingErrorRelative"] = (
self.matchingErrorRelative)
@property
def greedyTolMarginal(self):
"""Value of greedyTolMarginal."""
return self._greedyTolMarginal
@greedyTolMarginal.setter
def greedyTolMarginal(self, greedyTolMarginal):
if greedyTolMarginal < 0:
raise RROMPyException("greedyTolMarginal must be non-negative.")
if (hasattr(self, "_greedyTolMarginal")
and self.greedyTolMarginal is not None):
greedyTolMarginalold = self.greedyTolMarginal
else:
greedyTolMarginalold = -1
self._greedyTolMarginal = greedyTolMarginal
self._approxParameters["greedyTolMarginal"] = self.greedyTolMarginal
if greedyTolMarginalold != self.greedyTolMarginal:
self.resetSamples()
@property
def maxIterMarginal(self):
"""Value of maxIterMarginal."""
return self._maxIterMarginal
@maxIterMarginal.setter
def maxIterMarginal(self, maxIterMarginal):
if maxIterMarginal <= 0:
raise RROMPyException("maxIterMarginal must be positive.")
if (hasattr(self, "_maxIterMarginal")
and self.maxIterMarginal is not None):
maxIterMarginalold = self.maxIterMarginal
else:
maxIterMarginalold = -1
self._maxIterMarginal = maxIterMarginal
self._approxParameters["maxIterMarginal"] = self.maxIterMarginal
if maxIterMarginalold != self.maxIterMarginal:
self.resetSamples()
@property
def autoCollapse(self):
"""Value of autoCollapse."""
return self._autoCollapse
@autoCollapse.setter
def autoCollapse(self, autoCollapse):
self._autoCollapse = autoCollapse
self._approxParameters["autoCollapse"] = self.autoCollapse
def resetSamples(self):
"""Reset samples."""
super().resetSamples()
if not hasattr(self, "_temporaryPivot"):
self._mus = emptyParameterList()
self._musMarginal = emptyParameterList()
if hasattr(self, "samplerMarginal"): self.samplerMarginal.reset()
if hasattr(self, "samplingEngine") and self.samplingEngine is not None:
self.samplingEngine.resetHistory()
@abstractmethod
def getErrorEstimatorMarginalLookAhead(self) -> Np1D:
pass
def getErrorEstimatorMarginalNone(self) -> Np1D:
nErr = len(self.trainedModel.data.musMarginal)
self._musMarginalTestIdxs = np.arange(nErr)
return (1. + self.greedyTolMarginal) * np.ones(nErr)
def errorEstimatorMarginal(self, return_hi : bool = False) -> Np1D:
vbMng(self.trainedModel, "INIT",
"Evaluating error estimator at mu = {}.".format(
self.trainedModel.data.musMarginal), 10)
if self.errorEstimatorKindMarginal == "NONE":
nErr = len(self.trainedModel.data.musMarginal)
self._musMarginalTestIdxs = np.arange(nErr)
err = [None] * nErr
else:#if self.errorEstimatorKindMarginal[: 10] == "LOOK_AHEAD":
err = self.getErrorEstimatorMarginalLookAhead()
vbMng(self.trainedModel, "DEL", "Done evaluating error estimator.", 10)
if not return_hi: return err
if self.errorEstimatorKindMarginal == "NONE":
idxHiEst, hiErr = np.arange(len(err)), None
else:
if self.errorEstimatorKindMarginal == "LOOK_AHEAD_DOERFLER":
idxErr = np.argsort(err)
errC = np.cumsum(err[idxErr])
nRel = np.where(errC > self.greedyTolMarginal * errC[-1])[0]
if len(nRel) > 0:
idxHiEst = np.sort(idxErr[nRel[0] :])
else:
idxHiEst = []
else:#if self.errorEstimatorKindMarginal == "LOOK_AHEAD":
idxHiEst = np.where(err > self.greedyTolMarginal)[0]
hiErr = err[idxHiEst]
return err, idxHiEst, hiErr
def plotEstimatorMarginal(self, est:Np1D, idxHi:List[int],
estHi:List[float]):
if self.errorEstimatorKindMarginal == "NONE": return
if (not (np.any(np.isnan(est)) or np.any(np.isinf(est)))
and masterCore() and hasattr(self.trainedModel, "_musMExcl")):
fig = plt.figure(figsize = plt.figaspect(1. / self.nparMarginal))
for jpar in range(self.nparMarginal):
ax = fig.add_subplot(1, self.nparMarginal, 1 + jpar)
musre = np.real(self.trainedModel._musMExcl)
if len(idxHi) > 0 and estHi is not None:
hirej = musre[idxHi, jpar]
errCP = copy(est)
idx = np.delete(np.arange(self.nparMarginal), jpar)
while len(musre) > 0:
if self.nparMarginal == 1:
currIdx = np.arange(len(musre))
else:
currIdx = np.where(np.isclose(np.sum(
np.abs(musre[:, idx] - musre[0, idx]), 1),
0., atol = 1e-15))[0]
currIdxSorted = currIdx[np.argsort(musre[currIdx, jpar])]
ax.semilogy(musre[currIdxSorted, jpar],
errCP[currIdxSorted], 'k.-', linewidth = 1)
musre = np.delete(musre, currIdx, 0)
errCP = np.delete(errCP, currIdx)
if self.errorEstimatorKindMarginal == "LOOK_AHEAD_DOERFLER":
height = self.greedyTolMarginal * np.mean(est)
else: #if self.errorEstimatorKindMarginal == "LOOK_AHEAD":
height = self.greedyTolMarginal
ax.semilogy(self.musMarginal.re(jpar),
(height,) * len(self.musMarginal), '*m')
if len(idxHi) > 0 and estHi is not None:
ax.semilogy(hirej, estHi, 'xr')
ax.set_xlim(*list(self.samplerMarginal.lims.re(jpar)))
ax.grid()
plt.tight_layout()
plt.show()
@abstractmethod
def _updateTrainedModelMarginalSamples(self, idx : ListAny = []):
pass
def _addMarginalSample(self, mus:paramList):
mus = self.checkParameterListMarginal(mus)
if len(mus) == 0: return
self._nmusOld, nmus = len(self.musMarginal), len(mus)
if (hasattr(self, "trainedModel") and self.trainedModel is not None
and hasattr(self.trainedModel, "_musMExcl")):
self._nmusOld += len(self.trainedModel._musMExcl)
vbMng(self, "MAIN",
("Adding marginal sample point{} no. {}{} at {} to training "
"set.").format("s" * (nmus > 1), self._nmusOld + 1,
"--{}".format(self._nmusOld + nmus) * (nmus > 1),
mus), 3)
self.musMarginal.append(mus)
self.setupApproxPivoted(mus)
self._preliminaryMarginalFinalization()
del self._nmusOld
if (self.errorEstimatorKindMarginal[: 10] == "LOOK_AHEAD"
and not self.firstGreedyIterM):
ubRange = len(self.trainedModel.data.musMarginal)
if hasattr(self.trainedModel, "_idxExcl"):
shRange = len(self.trainedModel._musMExcl)
else:
shRange = 0
testIdxs = list(range(ubRange + shRange - len(mus),
ubRange + shRange))
for j in testIdxs[::-1]:
self.musMarginal.pop(j - shRange)
if hasattr(self.trainedModel, "_idxExcl"):
testIdxs = self.trainedModel._idxExcl + testIdxs
self._updateTrainedModelMarginalSamples(testIdxs)
self._finalizeMarginalization()
self._SMarginal = len(self.musMarginal)
self._approxParameters["SMarginal"] = self.SMarginal
self.trainedModel.data.approxParameters["SMarginal"] = self.SMarginal
def greedyNextSampleMarginal(self, muidx:List[int],
plotEst : str = "NONE") \
-> Tuple[Np1D, List[int], float, paramVal]:
RROMPyAssert(self._mode, message = "Cannot add greedy sample.")
muidx = self._musMarginalTestIdxs[muidx]
if (self.errorEstimatorKindMarginal[: 10] == "LOOK_AHEAD"
and not self.firstGreedyIterM):
if not hasattr(self.trainedModel, "_idxExcl"):
raise RROMPyException(("Sample index to be added not present "
"in trained model."))
testIdxs = copy(self.trainedModel._idxExcl)
skippedIdx = 0
for cj, j in enumerate(self.trainedModel._idxExcl):
if j in muidx:
testIdxs.pop(skippedIdx)
self.musMarginal.insert(self.trainedModel._musMExcl[cj],
j - skippedIdx)
else:
skippedIdx += 1
if len(self.trainedModel._idxExcl) < (len(muidx)
+ len(testIdxs)):
raise RROMPyException(("Sample index to be added not present "
"in trained model."))
self._updateTrainedModelMarginalSamples(testIdxs)
self._SMarginal = len(self.musMarginal)
self._approxParameters["SMarginal"] = self.SMarginal
self.trainedModel.data.approxParameters["SMarginal"] = (
self.SMarginal)
self.firstGreedyIterM = False
idxAdded = self.samplerMarginal.refine(muidx)[0]
self._addMarginalSample(self.samplerMarginal.points[idxAdded])
errorEstTest, muidx, hiErrorEst = self.errorEstimatorMarginal(True)
if plotEst == "ALL":
self.plotEstimatorMarginal(errorEstTest, muidx, hiErrorEst)
return (errorEstTest, muidx, hiErrorEst,
self.samplerMarginal.points[muidx])
def _preliminaryTrainingMarginal(self):
"""Initialize starting snapshots of solution map."""
RROMPyAssert(self._mode, message = "Cannot start greedy algorithm.")
if np.sum(self.samplingEngine.nsamples) > 0: return
self.resetSamples()
self._addMarginalSample(self.samplerMarginal.generatePoints(
self.SMarginal))
def _preSetupApproxPivoted(self, mus:paramList) \
-> Tuple[ListAny, ListAny, ListAny]:
self.computeScaleFactor()
if self.trainedModel is None:
self._setupTrainedModel(np.zeros((0, 0)))
self.trainedModel.data.Qs, self.trainedModel.data.Ps = [], []
self.trainedModel.data.Psupp = []
self._trainedModelOld = copy(self.trainedModel)
self._scaleFactorOldPivot = copy(self.scaleFactor)
self.scaleFactor = self.scaleFactorPivot
self._temporaryPivot = 1
self._musLoc = copy(self.mus)
+ is_master = masterCore()
idx, sizes = indicesScatter(len(mus), return_sizes = True)
- emptyCores = np.where(sizes == 0)[0]
+ emptyCores = np.where(np.array(sizes) == 0)[0] if is_master else []
self.verbosity -= 10
self.samplingEngine.verbosity -= 10
return idx, sizes, emptyCores
def _postSetupApproxPivoted(self, mus:Np2D, pMat:Np2D, Ps:ListAny,
Qs:ListAny, sizes:ListAny):
self.scaleFactor = self._scaleFactorOldPivot
del self._scaleFactorOldPivot, self._temporaryPivot
pMat, Ps, Qs, mus, nsamples = gatherPivotedApproximant(pMat, Ps, Qs,
mus, sizes,
self.polybasis)
if len(self._musLoc) > 0:
self._mus = self.checkParameterList(self._musLoc)
self._mus.append(mus)
else:
self._mus = self.checkParameterList(mus)
self.trainedModel = self._trainedModelOld
del self._trainedModelOld
if not self.matchState and self.autoCollapse:
pMat, padLeft, suppNew = 1., 0, [0] * len(nsamples)
else:
padLeft = self.trainedModel.data.projMat.shape[1]
suppNew = list(padLeft + np.append(0, np.cumsum(nsamples[: -1])))
self._setupTrainedModel(pMat, padLeft > 0,
collapsed = (not self.matchState
and self.autoCollapse))
self.trainedModel.data.Qs += Qs
self.trainedModel.data.Ps += Ps
self.trainedModel.data.Psupp += suppNew
self.trainedModel.data.approxParameters = copy(self.approxParameters)
self.verbosity += 10
self.samplingEngine.verbosity += 10
def _localPivotedResult(self, pMat:Np2D, req:ListAny, emptyCores:ListAny,
mus:Np2D) -> Tuple[Np2D, ListAny, Np2D]:
musi = self.samplingEngine.mus
if self.rationalMode[-6:] == "OUTPUT":
pMati = np.empty((1, 0), dtype = np.complex)
else:
pMati = self.samplingEngine.projectionMatrix
if not self.matchState:
if self.POD == 1 and not (
hasattr(self.HFEngine.C, "is_mu_independent")
and self.HFEngine.C.is_mu_independent in self._output_lvl):
raise RROMPyException(("Cannot apply mu-dependent C "
"to orthonormalized samples."))
vbMng(self, "INIT", "Extracting system output from state.", 35)
pMatiEff = None
for j, mu in enumerate(musi):
pMij = np.expand_dims(self.HFEngine.applyC(pMati[:, j],
mu), -1)
if pMatiEff is None:
pMatiEff = np.array(pMij)
else:
pMatiEff = np.append(pMatiEff, pMij, axis = 1)
pMati = pMatiEff
vbMng(self, "DEL", "Done extracting system output.", 35)
if pMat is None:
mus = copy(musi.data)
pMat = copy(pMati)
- if masterCore():
- for dest in emptyCores:
- req += [isend(len(pMat), dest = dest, tag = dest)]
+ for dest in emptyCores:
+ req += [isend(len(pMat), dest = dest, tag = dest)]
else:
mus = np.vstack((mus, musi.data))
if not self.matchState and self.autoCollapse:
pMat = pMati
else:
pMat = np.hstack((pMat, pMati))
return pMat, req, mus
@abstractmethod
def setupApproxPivoted(self, mus:paramList) -> int:
if self.checkComputedApproxPivoted(): return -1
RROMPyAssert(self._mode, message = "Cannot setup approximant.")
vbMng(self, "INIT", "Setting up pivoted approximant.", 10)
self._preSetupApproxPivoted()
data = []
pass
self._postSetupApproxPivoted(mus, data)
vbMng(self, "DEL", "Done setting up pivoted approximant.", 10)
return 0
def setupApprox(self, plotEst : str = "NONE") -> int:
"""Compute greedy snapshots of solution map."""
if self.checkComputedApprox(): return -1
if self.rationalMode[-6:] == "OUTPUT": self.autoCollapse = True
RROMPyAssert(self._mode, message = "Cannot start greedy algorithm.")
vbMng(self, "INIT", "Setting up {}.". format(self.name()), 5)
vbMng(self, "INIT", "Starting computation of snapshots.", 5)
max2ErrorEst, self.firstGreedyIterM = np.inf, True
self._preliminaryTrainingMarginal()
if self.errorEstimatorKindMarginal == "NONE":
muidx = []
else:#if self.errorEstimatorKindMarginal[: 10] == "LOOK_AHEAD":
muidx = np.arange(len(self.trainedModel.data.musMarginal))
self._musMarginalTestIdxs = np.array(muidx)
while self.firstGreedyIterM or (max2ErrorEst > self.greedyTolMarginal
and self.samplerMarginal.npoints < self.maxIterMarginal):
errorEstTest, muidx, hiErrorEst, mu = \
self.greedyNextSampleMarginal(muidx, plotEst)
if hiErrorEst is not None:
if len(hiErrorEst) > 0:
max2ErrorEst = np.max(hiErrorEst)
else:
max2ErrorEst = np.max(errorEstTest)
vbMng(self, "MAIN", ("Uniform testing error estimate "
"{:.4e}.").format(max2ErrorEst), 5)
if self.errorEstimatorKindMarginal in ["LOOK_AHEAD_DOERFLER",
"NONE"]:
max2ErrorEst = np.inf
if plotEst == "LAST":
self.plotEstimatorMarginal(errorEstTest, muidx, hiErrorEst)
vbMng(self, "DEL", ("Done computing snapshots (final snapshot count: "
"{}).").format(len(self.mus)), 5)
if (self.errorEstimatorKindMarginal != "NONE"
and hasattr(self.trainedModel, "_idxExcl")
and len(self.trainedModel._idxExcl) > 0):
vbMng(self, "INIT", "Recovering {} test models.".format(
len(self.trainedModel._idxExcl)), 7)
for j, mu in zip(self.trainedModel._idxExcl,
self.trainedModel._musMExcl):
self.musMarginal.insert(mu, j)
self._preliminaryMarginalFinalization()
self._updateTrainedModelMarginalSamples()
self._finalizeMarginalization()
self._SMarginal = len(self.musMarginal)
self._approxParameters["SMarginal"] = self.SMarginal
self.trainedModel.data.approxParameters["SMarginal"] = (
self.SMarginal)
vbMng(self, "DEL", "Done recovering test models.", 7)
vbMng(self, "DEL", "Done setting up approximant.", 5)
return 0
def checkComputedApproxPivoted(self) -> bool:
return (super().checkComputedApprox()
and len(self.musMarginal) == len(self.trainedModel.data.musMarginal))
class GenericPivotedGreedyApproximantPolyMatch(
GenericPivotedGreedyApproximantBase,
GenericPivotedApproximantPolyMatch):
"""
ROM pivoted greedy interpolant computation for parametric problems (with
polynomial matching) (ABSTRACT).
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': kind of snapshots orthogonalization; allowed values
include 0, 1/2, and 1; defaults to 1, i.e. POD;
- 'scaleFactorDer': scaling factors for derivative computation;
defaults to 'AUTO';
- 'matchState': whether to match the system state rather than the
system output; defaults to False;
- 'matchingWeight': weight for pole matching optimization; defaults
to 1;
- 'matchingKind': kind of matching; allowed values include 'ROTATE'
and 'PROJECT'; defaults to 'ROTATE';
- 'matchingWeightError': weight for matching in error estimation;
defaults to 0;
- 'matchingErrorRelative': whether error estimation is relative;
defaults to False, i.e. absolute error;
- 'S': total number of pivot samples current approximant relies
upon;
- 'samplerPivot': pivot sample point generator;
- 'SMarginal': number of starting marginal samples;
- 'samplerMarginal': marginal sample point generator via sparse
grid;
- 'errorEstimatorKindMarginal': kind of marginal error estimator;
available values include 'LOOK_AHEAD', 'LOOK_AHEAD_DOERFLER',
and 'NONE'; defaults to 'NONE';
- 'polybasisMarginal': type of polynomial basis for marginal
interpolation; allowed values include 'MONOMIAL_*',
'CHEBYSHEV_*', 'LEGENDRE_*', 'NEARESTNEIGHBOR', and
'PIECEWISE_LINEAR_*'; defaults to 'MONOMIAL';
- 'paramsMarginal': dictionary of parameters for marginal
interpolation; include:
. 'MMarginal': degree of marginal interpolant; defaults to
'AUTO', i.e. maximum allowed; not for 'NEARESTNEIGHBOR' or
'PIECEWISE_LINEAR_*';
. 'nNeighborsMarginal': number of marginal nearest neighbors;
defaults to 1; only for 'NEARESTNEIGHBOR';
. 'polydegreetypeMarginal': type of polynomial degree for
marginal; defaults to 'TOTAL'; not for 'NEARESTNEIGHBOR' or
'PIECEWISE_LINEAR_*';
. 'polyTruncateTolMarginal': tolerance for truncation of
marginal interpolator; defaults to 0;
. 'interpTolMarginal': tolerance for marginal interpolation;
defaults to None; not for 'NEARESTNEIGHBOR' or
'PIECEWISE_LINEAR_*';
. 'radialDirectionalWeightsMarginalAdapt': bounds for adaptive
rescaling of marginal radial basis weights; only for
radial basis.
- 'greedyTolMarginal': uniform error tolerance for marginal greedy
algorithm; defaults to 1e-1;
- 'maxIterMarginal': maximum number of marginal greedy steps;
defaults to 1e2;
- 'radialDirectionalWeightsMarginal': radial basis weights for
marginal interpolant; defaults to 1;
- 'autoCollapse': whether to collapse trained reduced model as soon
as it is built; defaults to False.
Defaults to empty dict.
verbosity(optional): Verbosity level. Defaults to 10.
Attributes:
HFEngine: HF problem solver.
mu0: Default parameter.
directionPivot: Pivot components.
mus: Array of snapshot parameters.
musMarginal: Array of marginal snapshot parameters.
approxParameters: Dictionary containing values for main parameters of
approximant. Recognized keys are in parameterList.
parameterListSoft: Recognized keys of soft approximant parameters:
- 'POD': kind of snapshots orthogonalization;
- 'scaleFactorDer': scaling factors for derivative computation;
- 'matchState': whether to match the system state rather than the
system output;
- 'matchingWeight': weight for matching;
- 'matchingKind': kind of matching;
- 'matchingWeightError': weight for matching in error estimation;
- 'matchingErrorRelative': whether error estimation is relative;
- 'errorEstimatorKindMarginal': kind of marginal error estimator;
- 'polybasisMarginal': type of polynomial basis for marginal
interpolation;
- 'paramsMarginal': dictionary of parameters for marginal
interpolation; include:
. 'MMarginal': degree of marginal interpolant;
. 'nNeighborsMarginal': number of marginal nearest neighbors;
. 'polydegreetypeMarginal': type of polynomial degree for
marginal;
. 'polyTruncateTolMarginal': tolerance for truncation of
marginal interpolator;
. 'interpTolMarginal': tolerance for marginal interpolation;
. 'radialDirectionalWeightsMarginalAdapt': bounds for adaptive
rescaling of marginal radial basis weights.
- 'greedyTolMarginal': uniform error tolerance for marginal greedy
algorithm;
- 'maxIterMarginal': maximum number of marginal greedy steps;
- 'radialDirectionalWeightsMarginal': radial basis weights for
marginal interpolant;
- 'autoCollapse': whether to collapse trained reduced model as soon
as it is built.
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 via sparse
grid.
verbosity: Verbosity level.
POD: Kind of snapshots orthogonalization.
scaleFactorDer: Scaling factors for derivative computation.
matchState: Whether to match the system state rather than the system
output.
matchingWeight: Weight for pole matching optimization.
matchingKind: Kind of matching.
matchingWeightError: Weight for pole matching optimization in error
estimation.
matchingErrorRelative: Whether error estimation is relative.
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 via sparse grid.
errorEstimatorKindMarginal: Kind of marginal error estimator.
polybasisMarginal: Type of polynomial basis for marginal interpolation.
paramsMarginal: Dictionary of parameters for marginal interpolation.
greedyTolMarginal: Uniform error tolerance for marginal greedy
algorithm.
maxIterMarginal: Maximum number of marginal greedy steps.
radialDirectionalWeightsMarginal: Radial basis weights for marginal
interpolant.
autoCollapse: Whether to collapse trained reduced model as soon as it
is built.
muBounds: 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.
"""
def getErrorEstimatorMarginalLookAhead(self) -> Np1D:
if not hasattr(self.trainedModel, "_musMExcl"):
err = np.zeros(0)
err[:] = np.inf
self._musMarginalTestIdxs = np.zeros(0, dtype = int)
return err
self._musMarginalTestIdxs = np.array(self.trainedModel._idxExcl,
dtype = int)
idx, sizes = indicesScatter(len(self.trainedModel._musMExcl),
return_sizes = True)
err = []
if len(idx) > 0:
self.verbosity -= 25
self.trainedModel.verbosity -= 25
for j in idx:
muTest = self.trainedModel._musMExcl[j]
QEx = self.trainedModel._QsExcl[j].coeffs
QAp = self.trainedModel.interpolateMarginalQ(muTest)[0]
no2Ex = np.sum(np.abs(QEx) ** 2.)
no2Ap = np.sum(np.abs(QAp) ** 2.)
inner = np.sum([QEx[j] * QAp[j].conj()
for j in range(min(len(QEx), len(QAp)))])
if self.matchingWeightError != 0:
PEx = self.trainedModel._PsExcl[j].coeffs.T
PAp = self.trainedModel.interpolateMarginalP(muTest)[0].T
if not self.trainedModel.data._collapsed:
PsuppTest = self.trainedModel._PsuppExcl[j]
PEx = dot(self.trainedModel.data.projMat[:,
PsuppTest : PsuppTest + PEx.shape[0]],
PEx)
PAp = dot(self.trainedModel.data.projMat, PAp)
no2PEx = np.sum(self.HFEngine.norm(PEx,
is_state = self.matchState) ** 2.)
no2PAp = np.sum(self.HFEngine.norm(PAp,
is_state = self.matchState) ** 2.)
innerP = np.sum([self.HFEngine.innerProduct(
PEx[:, j], PAp[:, j], is_state = self.matchState)
for j in range(min(PEx.shape[1], PAp.shape[1]))])
no2Ex = no2Ex + self.matchingWeightError * no2PEx
no2Ap = no2Ap + self.matchingWeightError * no2PAp
inner = inner + self.matchingWeightError * innerP
#component of QEx orthogonal to QAp
dist2 = no2Ex - np.abs(inner) ** 2. / no2Ap
if self.matchingErrorRelative:
dist2 /= no2Ex
else:
dist2 /= 1. + self.matchingWeightError
err += [dist2 ** .5]
self.verbosity += 25
self.trainedModel.verbosity += 25
return arrayGatherv(np.array(err), sizes)
def _updateTrainedModelMarginalSamples(self, idx : ListAny = []):
self.trainedModel.updateEffectiveSamples(idx, self.matchingWeight,
self.matchingKind,
self.HFEngine,
self.matchState)
def setupApprox(self, *args, **kwargs) -> int:
if self.checkComputedApprox(): return -1
self.purgeparamsMarginal()
_polybasisMarginal = self.polybasisMarginal
self._polybasisMarginal = ("PIECEWISE_LINEAR_"
+ self.samplerMarginal.kind)
setupOK = super().setupApprox(*args, **kwargs)
self._polybasisMarginal = _polybasisMarginal
if self.matchState: self._postApplyC()
return setupOK
class GenericPivotedGreedyApproximantPoleMatch(
GenericPivotedGreedyApproximantPolyMatch,
GenericPivotedApproximantPoleMatch):
"""
ROM pivoted greedy interpolant computation for parametric problems (with
pole matching) (ABSTRACT).
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': kind of snapshots orthogonalization; allowed values
include 0, 1/2, and 1; defaults to 1, i.e. POD;
- 'scaleFactorDer': scaling factors for derivative computation;
defaults to 'AUTO';
- 'matchState': whether to match the system state rather than the
system output; defaults to False;
- 'matchingWeight': weight for pole matching optimization; defaults
to 1;
- 'matchingShared': required ratio of marginal points to share
resonance; defaults to 1.;
- 'badPoleCorrection': strategy for correction of bad poles;
available values include 'ERASE', 'RATIONAL', and 'POLYNOMIAL';
defaults to 'ERASE';
- 'matchingWeightError': weight for pole matching optimization in
error estimation; defaults to 0;
- 'matchingErrorRelative': whether error estimation is relative;
defaults to False, i.e. absolute error;
- 'S': total number of pivot samples current approximant relies
upon;
- 'samplerPivot': pivot sample point generator;
- 'SMarginal': number of starting marginal samples;
- 'samplerMarginal': marginal sample point generator via sparse
grid;
- 'errorEstimatorKindMarginal': kind of marginal error estimator;
available values include 'LOOK_AHEAD', 'LOOK_AHEAD_DOERFLER',
and 'NONE'; defaults to 'NONE';
- 'polybasisMarginal': type of polynomial basis for marginal
interpolation; allowed values include 'MONOMIAL_*',
'CHEBYSHEV_*', 'LEGENDRE_*', 'NEARESTNEIGHBOR', and
'PIECEWISE_LINEAR_*'; defaults to 'MONOMIAL';
- 'paramsMarginal': dictionary of parameters for marginal
interpolation; include:
. 'MMarginal': degree of marginal interpolant; defaults to
'AUTO', i.e. maximum allowed; not for 'NEARESTNEIGHBOR' or
'PIECEWISE_LINEAR_*';
. 'nNeighborsMarginal': number of marginal nearest neighbors;
defaults to 1; only for 'NEARESTNEIGHBOR';
. 'polydegreetypeMarginal': type of polynomial degree for
marginal; defaults to 'TOTAL'; not for 'NEARESTNEIGHBOR' or
'PIECEWISE_LINEAR_*';
. 'polyTruncateTolMarginal': tolerance for truncation of
marginal interpolator; defaults to 0;
. 'interpTolMarginal': tolerance for marginal interpolation;
defaults to None; not for 'NEARESTNEIGHBOR' or
'PIECEWISE_LINEAR_*';
. 'radialDirectionalWeightsMarginalAdapt': bounds for adaptive
rescaling of marginal radial basis weights; only for
radial basis.
- 'greedyTolMarginal': uniform error tolerance for marginal greedy
algorithm; defaults to 1e-1;
- 'maxIterMarginal': maximum number of marginal greedy steps;
defaults to 1e2;
- 'radialDirectionalWeightsMarginal': radial basis weights for
marginal interpolant; defaults to 1;
- 'autoCollapse': whether to collapse trained reduced model as soon
as it is built; defaults to False.
Defaults to empty dict.
verbosity(optional): Verbosity level. Defaults to 10.
Attributes:
HFEngine: HF problem solver.
mu0: Default parameter.
directionPivot: Pivot components.
mus: Array of snapshot parameters.
musMarginal: Array of marginal snapshot parameters.
approxParameters: Dictionary containing values for main parameters of
approximant. Recognized keys are in parameterList.
parameterListSoft: Recognized keys of soft approximant parameters:
- 'POD': kind of snapshots orthogonalization;
- 'scaleFactorDer': scaling factors for derivative computation;
- 'matchState': whether to match the system state rather than the
system output;
- 'matchingWeight': weight for pole matching optimization;
- 'matchingShared': required ratio of marginal points to share
resonance;
- 'badPoleCorrection': strategy for correction of bad poles;
- 'matchingWeightError': weight for pole matching optimization in
error estimation;
- 'matchingErrorRelative': whether error estimation is relative;
- 'errorEstimatorKindMarginal': kind of marginal error estimator;
- 'polybasisMarginal': type of polynomial basis for marginal
interpolation;
- 'paramsMarginal': dictionary of parameters for marginal
interpolation; include:
. 'MMarginal': degree of marginal interpolant;
. 'nNeighborsMarginal': number of marginal nearest neighbors;
. 'polydegreetypeMarginal': type of polynomial degree for
marginal;
. 'polyTruncateTolMarginal': tolerance for truncation of
marginal interpolator;
. 'interpTolMarginal': tolerance for marginal interpolation;
. 'radialDirectionalWeightsMarginalAdapt': bounds for adaptive
rescaling of marginal radial basis weights.
- 'greedyTolMarginal': uniform error tolerance for marginal greedy
algorithm;
- 'maxIterMarginal': maximum number of marginal greedy steps;
- 'radialDirectionalWeightsMarginal': radial basis weights for
marginal interpolant;
- 'autoCollapse': whether to collapse trained reduced model as soon
as it is built.
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 via sparse
grid.
verbosity: Verbosity level.
POD: Kind of snapshots orthogonalization.
scaleFactorDer: Scaling factors for derivative computation.
matchState: Whether to match the system state rather than the system
output.
matchingWeight: Weight for pole matching optimization.
matchingShared: Required ratio of marginal points to share resonance.
badPoleCorrection: Strategy for correction of bad poles.
matchingWeightError: Weight for pole matching optimization in error
estimation.
matchingErrorRelative: Whether error estimation is relative.
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 via sparse grid.
errorEstimatorKindMarginal: Kind of marginal error estimator.
polybasisMarginal: Type of polynomial basis for marginal interpolation.
paramsMarginal: Dictionary of parameters for marginal interpolation.
greedyTolMarginal: Uniform error tolerance for marginal greedy
algorithm.
maxIterMarginal: Maximum number of marginal greedy steps.
radialDirectionalWeightsMarginal: Radial basis weights for marginal
interpolant.
autoCollapse: Whether to collapse trained reduced model as soon as it
is built.
muBounds: 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.
"""
def getErrorEstimatorMarginalLookAhead(self) -> Np1D:
if not hasattr(self.trainedModel, "_musMExcl"):
err = np.zeros(0)
err[:] = np.inf
self._musMarginalTestIdxs = np.zeros(0, dtype = int)
return err
self._musMarginalTestIdxs = np.array(self.trainedModel._idxExcl,
dtype = int)
idx, sizes = indicesScatter(len(self.trainedModel._musMExcl),
return_sizes = True)
err = []
if len(idx) > 0:
self.verbosity -= 25
self.trainedModel.verbosity -= 25
for j in idx:
muTest = self.trainedModel._musMExcl[j]
HITest = self.trainedModel._HIsExcl[j]
polesEx = HITest.poles
idxGood = np.isinf(polesEx) + np.isnan(polesEx) == False
polesEx = polesEx[idxGood]
if len(polesEx) == 0:
err += [0.]
continue
polesAp = self.trainedModel.interpolateMarginalPoles(muTest)[0]
if self.matchingWeightError != 0:
resEx = HITest.coeffs[np.where(idxGood)[0]].T
resAp = self.trainedModel.interpolateMarginalCoeffs(
muTest)[0][: len(polesAp), :].T
if not self.trainedModel.data._collapsed:
PsuppTest = self.trainedModel._PsuppExcl[j]
resEx = dot(self.trainedModel.data.projMat[:,
PsuppTest : PsuppTest + resEx.shape[0]],
resEx)
resAp = dot(self.trainedModel.data.projMat, resAp)
else:
resEx, resAp = None, None
#match Ap to Ex
distMat = doubleDistanceMatrix(polesEx, polesAp,
self.matchingWeightError,
resEx, resAp, self.HFEngine,
self.matchState)
pmR, pmC = pointMatching(distMat)
dist = np.linalg.norm(distMat[pmR, pmC].flatten())
if self.matchingErrorRelative:
if self.matchingWeightError != 0:
resEx0 = resEx[:, pmR]
res0 = np.zeros_like(resEx[:, [0]])
else:
resEx0, res0 = None, None
dist0 = doubleDistanceMatrix(polesEx[pmR], [0.],
self.matchingWeightError,
resEx0, res0, self.HFEngine,
self.matchState).flatten()
dist /= np.linalg.norm(dist0)
err += [dist]
self.verbosity += 25
self.trainedModel.verbosity += 25
return arrayGatherv(np.array(err), sizes)
def _updateTrainedModelMarginalSamples(self, idx : ListAny = []):
self.trainedModel.updateEffectiveSamples(idx, self.matchingWeight,
self.HFEngine,
self.matchState)
diff --git a/rrompy/reduction_methods/pivoted/greedy/rational_interpolant_greedy_pivoted_greedy.py b/rrompy/reduction_methods/pivoted/greedy/rational_interpolant_greedy_pivoted_greedy.py
index 7198768..5a188a8 100644
--- a/rrompy/reduction_methods/pivoted/greedy/rational_interpolant_greedy_pivoted_greedy.py
+++ b/rrompy/reduction_methods/pivoted/greedy/rational_interpolant_greedy_pivoted_greedy.py
@@ -1,582 +1,551 @@
#Copyright (C) 2018-2020 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_greedy_approximant import (
GenericPivotedGreedyApproximantBase,
GenericPivotedGreedyApproximantPolyMatch,
GenericPivotedGreedyApproximantPoleMatch)
from rrompy.reduction_methods.standard.greedy import RationalInterpolantGreedy
from rrompy.reduction_methods.standard.greedy.generic_greedy_approximant \
import pruneSamples
from rrompy.reduction_methods.pivoted import (
RationalInterpolantGreedyPivotedPolyMatch,
RationalInterpolantGreedyPivotedPoleMatch)
-from rrompy.utilities.base.types import Np1D, Tuple, paramVal, paramList
+from rrompy.utilities.base.types import paramList
from rrompy.utilities.base import verbosityManager as vbMng
from rrompy.utilities.exception_manager import RROMPyAssert
from rrompy.parameter import emptyParameterList
from rrompy.utilities.parallel import poolRank, recv
__all__ = ['RationalInterpolantGreedyPivotedGreedyPolyMatch',
'RationalInterpolantGreedyPivotedGreedyPoleMatch']
class RationalInterpolantGreedyPivotedGreedyBase(
GenericPivotedGreedyApproximantBase):
- def greedyNextSample(self, muidx:int, plotEst : str = "NONE")\
- -> Tuple[Np1D, int, float, paramVal]:
- """Compute next greedy snapshot of solution map."""
- RROMPyAssert(self._mode, message = "Cannot add greedy sample.")
- mus = copy(self.muTest[muidx])
- self.muTest.pop(muidx)
- for j, mu in enumerate(mus):
- vbMng(self, "MAIN",
- ("Adding sample point no. {} at {} to training "
- "set.").format(len(self.mus) + 1, mu), 3)
- self.mus.append(mu)
- self._S = len(self.mus)
- self._approxParameters["S"] = self.S
- if (self.samplingEngine.nsamples <= len(mus) - j - 1
- or not np.allclose(mu, self.samplingEngine.mus[j - len(mus)])):
- self.samplingEngine.nextSample(mu)
- if self._isLastSampleCollinear():
- vbMng(self, "MAIN",
- ("Collinearity above tolerance detected. Starting "
- "preemptive greedy loop termination."), 3)
- self._collinearityFlag = 1
- errorEstTest = np.empty(len(self.muTest))
- errorEstTest[:] = np.nan
- return errorEstTest, [-1], np.nan, np.nan
- errorEstTest, muidx, maxErrorEst = self.errorEstimator(self.muTest,
- True)
- if plotEst == "ALL":
- self.plotEstimator(errorEstTest, muidx, maxErrorEst)
- return errorEstTest, muidx, maxErrorEst, self.muTest[muidx]
-
def _preliminaryTraining(self):
"""Initialize starting snapshots of solution map."""
RROMPyAssert(self._mode, message = "Cannot start greedy algorithm.")
if self.samplingEngine.nsamples > 0: return
self.resetSamples()
self.samplingEngine.scaleFactor = self.scaleFactorDer
musPivot = self.samplerTrainSet.generatePoints(self.S)
while len(musPivot) > self.S: musPivot.pop()
muTestBasePivot = self.samplerPivot.generatePoints(self.nTestPoints,
False)
idxPop = pruneSamples(self.mapParameterListPivot(muTestBasePivot),
self.mapParameterListPivot(musPivot),
1e-10 * self.scaleFactorPivot[0])
muTestBasePivot.pop(idxPop)
self._mus = emptyParameterList()
self.mus.reset((self.S - 1, self.HFEngine.npar))
self.muTest = emptyParameterList()
self.muTest.reset((len(muTestBasePivot) + 1, self.HFEngine.npar))
self.mus.data[:, self.directionPivot] = musPivot[: -1]
self.mus.data[:, self.directionMarginal] = np.repeat(self.muMargLoc,
self.S - 1, axis = 0)
self.muTest.data[: -1, self.directionPivot] = muTestBasePivot.data
self.muTest.data[-1, self.directionPivot] = musPivot[-1]
self.muTest.data[:, self.directionMarginal] = np.repeat(self.muMargLoc,
len(muTestBasePivot) + 1,
axis = 0)
- if len(self.mus) > 0:
- vbMng(self, "MAIN",
- ("Adding first {} sample point{} at {} to training "
- "set.").format(self.S - 1, "" + "s" * (self.S > 2),
- self.mus), 3)
- self.samplingEngine.iterSample(self.mus)
+ vbMng(self, "MAIN",
+ ("Adding first {} sample point{} at {} to training "
+ "set.").format(self.S - 1, "" + "s" * (self.S > 2), self.mus),
+ 3)
+ self.samplingEngine.iterSample(self.mus)
self._S = len(self.mus)
self._approxParameters["S"] = self.S
def setupApproxPivoted(self, mus:paramList) -> int:
if self.checkComputedApproxPivoted(): return -1
RROMPyAssert(self._mode, message = "Cannot setup approximant.")
vbMng(self, "INIT", "Setting up pivoted approximant.", 10)
if not hasattr(self, "_plotEstPivot"): self._plotEstPivot = "NONE"
idx, sizes, emptyCores = self._preSetupApproxPivoted(mus)
S0 = copy(self.S)
pMat, Ps, Qs, req, musA = None, [], [], [], None
if len(idx) == 0:
vbMng(self, "MAIN", "Idling.", 45)
if self.storeAllSamples: self.storeSamples()
pL = recv(source = 0, tag = poolRank())
pMat = np.empty((pL, 0), dtype = np.complex)
musA = np.empty((0, self.mu0.shape[1]), dtype = np.complex)
else:
for i in idx:
self.muMargLoc = mus[[i]]
vbMng(self, "MAIN", "Building marginal model no. {} at "
"{}.".format(i + 1, self.muMargLoc[0]), 25)
self.samplingEngine.resetHistory()
self.trainedModel = None
self.verbosity -= 5
self.samplingEngine.verbosity -= 5
RationalInterpolantGreedy.setupApprox(self, self._plotEstPivot)
self.verbosity += 5
self.samplingEngine.verbosity += 5
if self.storeAllSamples: self.storeSamples(i + self._nmusOld)
pMat, req, musA = self._localPivotedResult(pMat, req,
emptyCores, musA)
Ps += [copy(self.trainedModel.data.P)]
Qs += [copy(self.trainedModel.data.Q)]
if not self.matchState and self.autoCollapse:
Ps[-1].postmultiplyTensorize(pMat.T)
self._S = S0
del self.muMargLoc
for r in req: r.wait()
if not self.matchState and self.autoCollapse: pMat = pMat[:, : 0]
self._postSetupApproxPivoted(musA, pMat, Ps, Qs, sizes)
vbMng(self, "DEL", "Done setting up pivoted approximant.", 10)
return 0
def setupApprox(self, plotEst : str = "NONE") -> int:
if self.checkComputedApprox(): return -1
if '_' not in plotEst: plotEst = plotEst + "_NONE"
plotEstM, self._plotEstPivot = plotEst.split("_")
val = super().setupApprox(plotEstM)
return val
class RationalInterpolantGreedyPivotedGreedyPolyMatch(
RationalInterpolantGreedyPivotedGreedyBase,
GenericPivotedGreedyApproximantPolyMatch,
RationalInterpolantGreedyPivotedPolyMatch):
"""
ROM greedy pivoted greedy rational interpolant computation for parametric
problems (with polynomial 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': kind of snapshots orthogonalization; allowed values
include 0, 1/2, and 1; defaults to 1, i.e. POD;
- 'scaleFactorDer': scaling factors for derivative computation;
defaults to 'AUTO';
- 'matchState': whether to match the system state rather than the
system output; defaults to False;
- 'matchingWeight': weight for matching; defaults to 1;
- 'matchingKind': kind of matching; allowed values include 'ROTATE'
and 'PROJECT'; defaults to 'ROTATE';
- 'matchingWeightError': weight for matching in error estimation;
defaults to 0;
- 'matchingErrorRelative': whether error estimation is relative;
defaults to False, i.e. absolute error;
- 'S': total number of pivot samples current approximant relies
upon;
- 'N': degree of rational interpolant denominator; defaults to
'AUTO', i.e. maximum allowed;
- 'samplerPivot': pivot sample point generator;
- 'SMarginal': number of starting marginal samples;
- 'samplerMarginal': marginal sample point generator via sparse
grid;
- 'errorEstimatorKindMarginal': kind of marginal error estimator;
available values include 'LOOK_AHEAD', 'LOOK_AHEAD_DOERFLER',
and 'NONE'; defaults to 'NONE';
- 'polybasis': type of polynomial basis for pivot interpolation;
defaults to 'MONOMIAL';
- 'polybasisMarginal': type of polynomial basis for marginal
interpolation; allowed values include 'MONOMIAL_*',
'CHEBYSHEV_*', 'LEGENDRE_*', 'NEARESTNEIGHBOR', and
'PIECEWISE_LINEAR_*'; defaults to 'MONOMIAL';
- 'paramsMarginal': dictionary of parameters for marginal
interpolation; include:
. 'MMarginal': degree of marginal interpolant; defaults to
'AUTO', i.e. maximum allowed; not for 'NEARESTNEIGHBOR' or
'PIECEWISE_LINEAR_*';
. 'nNeighborsMarginal': number of marginal nearest neighbors;
defaults to 1; only for 'NEARESTNEIGHBOR';
. 'polydegreetypeMarginal': type of polynomial degree for
marginal; defaults to 'TOTAL'; not for 'NEARESTNEIGHBOR' or
'PIECEWISE_LINEAR_*';
. 'polyTruncateTolMarginal': tolerance for truncation of
marginal interpolator; defaults to 0;
. 'interpTolMarginal': tolerance for marginal interpolation;
defaults to None; not for 'NEARESTNEIGHBOR' or
'PIECEWISE_LINEAR_*';
. 'radialDirectionalWeightsMarginalAdapt': bounds for adaptive
rescaling of marginal radial basis weights; only for
radial basis.
- 'greedyTol': uniform error tolerance for greedy algorithm;
defaults to 1e-2;
- 'collinearityTol': collinearity tolerance for greedy algorithm;
defaults to 0.;
- 'maxIter': maximum number of greedy steps; defaults to 1e2;
- 'nTestPoints': number of test points; defaults to 5e2;
- 'samplerTrainSet': training sample points generator; defaults to
samplerPivot;
- 'errorEstimatorKind': kind of error estimator; available values
include 'AFFINE', 'DISCREPANCY', 'LOOK_AHEAD',
'LOOK_AHEAD_RES', and 'NONE'; defaults to 'NONE';
- 'greedyTolMarginal': uniform error tolerance for marginal greedy
algorithm; defaults to 1e-1;
- 'maxIterMarginal': maximum number of marginal greedy steps;
defaults to 1e2;
- 'radialDirectionalWeightsMarginal': radial basis weights for
marginal interpolant; defaults to 1;
- 'autoCollapse': whether to collapse trained reduced model as soon
as it is built; defaults to False;
- 'functionalSolve': strategy for minimization of denominator
functional; allowed values include 'NORM', 'DOMINANT',
'BARYCENTRIC[_NORM]', and 'BARYCENTRIC_AVERAGE' (check pdf in
main folder for explanation); defaults to 'NORM';
- 'interpTol': tolerance for pivot interpolation; defaults to None;
- 'forceQReal': force denominator to have real coefficients;
defaults to False;
- 'polyTruncateTol': tolerance for truncation of rational terms;
defaults to 0;
- 'QTol': tolerance for robust rational denominator management;
defaults to 0.
Defaults to empty dict.
verbosity(optional): Verbosity level. Defaults to 10.
Attributes:
HFEngine: HF problem solver.
mu0: Default parameter.
directionPivot: Pivot components.
mus: Array of snapshot parameters.
musMarginal: Array of marginal snapshot parameters.
approxParameters: Dictionary containing values for main parameters of
approximant. Recognized keys are in parameterList.
parameterListSoft: Recognized keys of soft approximant parameters:
- 'POD': kind of snapshots orthogonalization;
- 'scaleFactorDer': scaling factors for derivative computation;
- 'matchState': whether to match the system state rather than the
system output;
- 'matchingWeight': weight for matching;
- 'matchingKind': kind of matching;
- 'matchingWeightError': weight for matching in error estimation;
- 'matchingErrorRelative': whether error estimation is relative;
- 'N': degree of rational interpolant denominator;
- 'errorEstimatorKindMarginal': kind of marginal error estimator;
- 'polybasis': type of polynomial basis for pivot interpolation;
- 'polybasisMarginal': type of polynomial basis for marginal
interpolation;
- 'paramsMarginal': dictionary of parameters for marginal
interpolation; include:
. 'MMarginal': degree of marginal interpolant;
. 'nNeighborsMarginal': number of marginal nearest neighbors;
. 'polydegreetypeMarginal': type of polynomial degree for
marginal;
. 'polyTruncateTolMarginal': tolerance for truncation of
marginal interpolator;
. 'interpTolMarginal': tolerance for marginal interpolation;
. 'radialDirectionalWeightsMarginalAdapt': bounds for adaptive
rescaling of marginal radial basis weights.
- 'greedyTol': uniform error tolerance for greedy algorithm;
- 'collinearityTol': collinearity tolerance for greedy algorithm;
- 'maxIter': maximum number of greedy steps;
- 'nTestPoints': number of test points;
- 'samplerTrainSet': training sample points generator;
- 'errorEstimatorKind': kind of error estimator;
- 'greedyTolMarginal': uniform error tolerance for marginal greedy
algorithm;
- 'maxIterMarginal': maximum number of marginal greedy steps;
- 'radialDirectionalWeightsMarginal': radial basis weights for
marginal interpolant;
- 'autoCollapse': whether to collapse trained reduced model as soon
as it is built;
- 'functionalSolve': strategy for minimization of denominator
functional;
- 'interpTol': tolerance for pivot interpolation;
- 'forceQReal': force denominator to have real coefficients;
- 'polyTruncateTol': tolerance for truncation of rational terms;
- 'QTol': 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 via sparse
grid.
verbosity: Verbosity level.
POD: Kind of snapshots orthogonalization.
scaleFactorDer: Scaling factors for derivative computation.
matchState: Whether to match the system state rather than the system
output.
matchingWeight: Weight for matching.
matchingKind: Kind of matching.
matchingWeightError: Weight for matching in error estimation.
matchingErrorRelative: Whether error estimation is relative.
S: Total number of pivot samples current approximant relies upon.
N: Denominator degree of approximant.
samplerPivot: Pivot sample point generator.
SMarginal: Total number of marginal samples current approximant relies
upon.
samplerMarginal: Marginal sample point generator via sparse grid.
errorEstimatorKindMarginal: Kind of marginal error estimator.
polybasis: Type of polynomial basis for pivot interpolation.
polybasisMarginal: Type of polynomial basis for marginal interpolation.
paramsMarginal: Dictionary of parameters for marginal interpolation.
greedyTol: uniform error tolerance for greedy algorithm.
collinearityTol: Collinearity tolerance for greedy algorithm.
maxIter: maximum number of greedy steps.
nTestPoints: number of starting training points.
samplerTrainSet: training sample points generator.
errorEstimatorKind: kind of error estimator.
greedyTolMarginal: Uniform error tolerance for marginal greedy
algorithm.
maxIterMarginal: Maximum number of marginal greedy steps.
radialDirectionalWeightsMarginal: Radial basis weights for marginal
interpolant.
autoCollapse: Whether to collapse trained reduced model as soon as it
is built.
functionalSolve: Strategy for minimization of denominator functional.
interpTol: Tolerance for pivot interpolation.
forceQReal: Force denominator to have real coefficients.
polyTruncateTol: Tolerance for truncation of rational terms.
QTol: Tolerance for robust rational denominator management.
muBounds: 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.
"""
class RationalInterpolantGreedyPivotedGreedyPoleMatch(
RationalInterpolantGreedyPivotedGreedyBase,
GenericPivotedGreedyApproximantPoleMatch,
RationalInterpolantGreedyPivotedPoleMatch):
"""
ROM greedy pivoted greedy 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': kind of snapshots orthogonalization; allowed values
include 0, 1/2, and 1; defaults to 1, i.e. POD;
- 'scaleFactorDer': scaling factors for derivative computation;
defaults to 'AUTO';
- 'matchState': whether to match the system state rather than the
system output; defaults to False;
- 'matchingWeight': weight for pole matching optimization; defaults
to 1;
- 'matchingShared': required ratio of marginal points to share
resonance; defaults to 1.;
- 'badPoleCorrection': strategy for correction of bad poles;
available values include 'ERASE', 'RATIONAL', and 'POLYNOMIAL';
defaults to 'ERASE';
- 'matchingWeightError': weight for pole matching optimization in
error estimation; defaults to 0;
- 'matchingErrorRelative': whether error estimation is relative;
defaults to False, i.e. absolute error;
- 'S': total number of pivot samples current approximant relies
upon;
- 'N': degree of rational interpolant denominator; defaults to
'AUTO', i.e. maximum allowed;
- 'samplerPivot': pivot sample point generator;
- 'SMarginal': number of starting marginal samples;
- 'samplerMarginal': marginal sample point generator via sparse
grid;
- 'errorEstimatorKindMarginal': kind of marginal error estimator;
available values include 'LOOK_AHEAD', 'LOOK_AHEAD_DOERFLER',
and 'NONE'; defaults to 'NONE';
- 'polybasis': type of polynomial basis for pivot interpolation;
defaults to 'MONOMIAL';
- 'polybasisMarginal': type of polynomial basis for marginal
interpolation; allowed values include 'MONOMIAL_*',
'CHEBYSHEV_*', 'LEGENDRE_*', 'NEARESTNEIGHBOR', and
'PIECEWISE_LINEAR_*'; defaults to 'MONOMIAL';
- 'paramsMarginal': dictionary of parameters for marginal
interpolation; include:
. 'MMarginal': degree of marginal interpolant; defaults to
'AUTO', i.e. maximum allowed; not for 'NEARESTNEIGHBOR' or
'PIECEWISE_LINEAR_*';
. 'nNeighborsMarginal': number of marginal nearest neighbors;
defaults to 1; only for 'NEARESTNEIGHBOR';
. 'polydegreetypeMarginal': type of polynomial degree for
marginal; defaults to 'TOTAL'; not for 'NEARESTNEIGHBOR' or
'PIECEWISE_LINEAR_*';
. 'polyTruncateTolMarginal': tolerance for truncation of
marginal interpolator; defaults to 0;
. 'interpTolMarginal': tolerance for marginal interpolation;
defaults to None; not for 'NEARESTNEIGHBOR' or
'PIECEWISE_LINEAR_*';
. 'radialDirectionalWeightsMarginalAdapt': bounds for adaptive
rescaling of marginal radial basis weights; only for
radial basis.
- 'greedyTol': uniform error tolerance for greedy algorithm;
defaults to 1e-2;
- 'collinearityTol': collinearity tolerance for greedy algorithm;
defaults to 0.;
- 'maxIter': maximum number of greedy steps; defaults to 1e2;
- 'nTestPoints': number of test points; defaults to 5e2;
- 'samplerTrainSet': training sample points generator; defaults to
samplerPivot;
- 'errorEstimatorKind': kind of error estimator; available values
include 'AFFINE', 'DISCREPANCY', 'LOOK_AHEAD',
'LOOK_AHEAD_RES', and 'NONE'; defaults to 'NONE';
- 'greedyTolMarginal': uniform error tolerance for marginal greedy
algorithm; defaults to 1e-1;
- 'maxIterMarginal': maximum number of marginal greedy steps;
defaults to 1e2;
- 'radialDirectionalWeightsMarginal': radial basis weights for
marginal interpolant; defaults to 1;
- 'autoCollapse': whether to collapse trained reduced model as soon
as it is built; defaults to False;
- 'functionalSolve': strategy for minimization of denominator
functional; allowed values include 'NORM', 'DOMINANT',
'BARYCENTRIC[_NORM]', and 'BARYCENTRIC_AVERAGE' (check pdf in
main folder for explanation); defaults to 'NORM';
- 'interpTol': tolerance for pivot interpolation; defaults to None;
- 'forceQReal': force denominator to have real coefficients;
defaults to False;
- 'polyTruncateTol': tolerance for truncation of rational terms;
defaults to 0;
- 'QTol': tolerance for robust rational denominator management;
defaults to 0.
Defaults to empty dict.
verbosity(optional): Verbosity level. Defaults to 10.
Attributes:
HFEngine: HF problem solver.
mu0: Default parameter.
directionPivot: Pivot components.
mus: Array of snapshot parameters.
musMarginal: Array of marginal snapshot parameters.
approxParameters: Dictionary containing values for main parameters of
approximant. Recognized keys are in parameterList.
parameterListSoft: Recognized keys of soft approximant parameters:
- 'POD': kind of snapshots orthogonalization;
- 'scaleFactorDer': scaling factors for derivative computation;
- 'matchState': whether to match the system state rather than the
system output;
- 'matchingWeight': weight for pole matching optimization;
- 'matchingShared': required ratio of marginal points to share
resonance;
- 'badPoleCorrection': strategy for correction of bad poles;
- 'matchingWeightError': weight for pole matching optimization in
error estimation;
- 'matchingErrorRelative': whether error estimation is relative;
- 'N': degree of rational interpolant denominator;
- 'errorEstimatorKindMarginal': kind of marginal error estimator;
- 'polybasis': type of polynomial basis for pivot interpolation;
- 'polybasisMarginal': type of polynomial basis for marginal
interpolation;
- 'paramsMarginal': dictionary of parameters for marginal
interpolation; include:
. 'MMarginal': degree of marginal interpolant;
. 'nNeighborsMarginal': number of marginal nearest neighbors;
. 'polydegreetypeMarginal': type of polynomial degree for
marginal;
. 'polyTruncateTolMarginal': tolerance for truncation of
marginal interpolator;
. 'interpTolMarginal': tolerance for marginal interpolation;
. 'radialDirectionalWeightsMarginalAdapt': bounds for adaptive
rescaling of marginal radial basis weights.
- 'greedyTol': uniform error tolerance for greedy algorithm;
- 'collinearityTol': collinearity tolerance for greedy algorithm;
- 'maxIter': maximum number of greedy steps;
- 'nTestPoints': number of test points;
- 'samplerTrainSet': training sample points generator;
- 'errorEstimatorKind': kind of error estimator;
- 'greedyTolMarginal': uniform error tolerance for marginal greedy
algorithm;
- 'maxIterMarginal': maximum number of marginal greedy steps;
- 'radialDirectionalWeightsMarginal': radial basis weights for
marginal interpolant;
- 'autoCollapse': whether to collapse trained reduced model as soon
as it is built;
- 'functionalSolve': strategy for minimization of denominator
functional;
- 'interpTol': tolerance for pivot interpolation;
- 'forceQReal': force denominator to have real coefficients;
- 'polyTruncateTol': tolerance for truncation of rational terms;
- 'QTol': 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 via sparse
grid.
verbosity: Verbosity level.
POD: Kind of snapshots orthogonalization.
scaleFactorDer: Scaling factors for derivative computation.
matchState: Whether to match the system state rather than the system
output.
matchingWeight: Weight for pole matching optimization.
matchingShared: Required ratio of marginal points to share resonance.
badPoleCorrection: Strategy for correction of bad poles.
matchingWeightError: Weight for pole matching optimization in error
estimation.
matchingErrorRelative: Whether error estimation is relative.
S: Total number of pivot samples current approximant relies upon.
N: Denominator degree of approximant.
samplerPivot: Pivot sample point generator.
SMarginal: Total number of marginal samples current approximant relies
upon.
samplerMarginal: Marginal sample point generator via sparse grid.
errorEstimatorKindMarginal: Kind of marginal error estimator.
polybasis: Type of polynomial basis for pivot interpolation.
polybasisMarginal: Type of polynomial basis for marginal interpolation.
paramsMarginal: Dictionary of parameters for marginal interpolation.
greedyTol: uniform error tolerance for greedy algorithm.
collinearityTol: Collinearity tolerance for greedy algorithm.
maxIter: maximum number of greedy steps.
nTestPoints: number of starting training points.
samplerTrainSet: training sample points generator.
errorEstimatorKind: kind of error estimator.
greedyTolMarginal: Uniform error tolerance for marginal greedy
algorithm.
maxIterMarginal: Maximum number of marginal greedy steps.
radialDirectionalWeightsMarginal: Radial basis weights for marginal
interpolant.
autoCollapse: Whether to collapse trained reduced model as soon as it
is built.
functionalSolve: Strategy for minimization of denominator functional.
interpTol: Tolerance for pivot interpolation.
forceQReal: Force denominator to have real coefficients.
polyTruncateTol: Tolerance for truncation of rational terms.
QTol: Tolerance for robust rational denominator management.
muBounds: 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.
"""
diff --git a/rrompy/reduction_methods/pivoted/rational_interpolant_greedy_pivoted.py b/rrompy/reduction_methods/pivoted/rational_interpolant_greedy_pivoted.py
index d1b30a9..c8b8f7e 100644
--- a/rrompy/reduction_methods/pivoted/rational_interpolant_greedy_pivoted.py
+++ b/rrompy/reduction_methods/pivoted/rational_interpolant_greedy_pivoted.py
@@ -1,824 +1,830 @@
# Copyright (C) 2018-2020 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 (GenericPivotedApproximantBase,
GenericPivotedApproximantNoMatch,
GenericPivotedApproximantPolyMatch,
GenericPivotedApproximantPoleMatch)
from .gather_pivoted_approximant import gatherPivotedApproximant
from rrompy.reduction_methods.standard.greedy.rational_interpolant_greedy \
import RationalInterpolantGreedy
from rrompy.reduction_methods.standard.greedy.generic_greedy_approximant \
import pruneSamples
from rrompy.utilities.base.types import Np1D, paramList
from rrompy.utilities.base import verbosityManager as vbMng
from rrompy.utilities.poly_fitting.polynomial import polyvander as pv
from rrompy.utilities.poly_fitting.piecewise_linear import sparsekinds as sk
from rrompy.utilities.exception_manager import (RROMPyException, RROMPyAssert,
RROMPyWarning)
from rrompy.parameter import emptyParameterList
from rrompy.utilities.parallel import (masterCore, poolRank, indicesScatter,
isend, recv)
__all__ = ['RationalInterpolantGreedyPivotedNoMatch',
'RationalInterpolantGreedyPivotedPolyMatch',
'RationalInterpolantGreedyPivotedPoleMatch']
class RationalInterpolantGreedyPivotedBase(GenericPivotedApproximantBase,
RationalInterpolantGreedy):
def __init__(self, *args, **kwargs):
self._preInit()
+ self._addParametersToList(toBeExcluded = ["MAuxiliary", "NAuxiliary"])
super().__init__(*args, **kwargs)
if self.nparPivot > 1: self.HFEngine._ignoreResidues = 1
self._postInit()
@property
def tModelType(self):
if hasattr(self, "_temporaryPivot"):
return RationalInterpolantGreedy.tModelType.fget(self)
return super().tModelType
- def _polyvanderAuxiliary(self, mus, deg, *args):
+ @property
+ def MAuxiliary(self): return 0
+
+ @property
+ def NAuxiliary(self): return 0
+
+ def _polyvanderAuxiliary(self, mus, deg, *args, **kwargs):
degEff = [0] * self.npar
- degEff[self.directionPivot[0]] = deg
- return pv(mus, degEff, *args)
+ degEff[self.directionPivot[0]] = deg[0]
+ return pv(mus, degEff, *args, **kwargs)
def _marginalizeMiscellanea(self, forward:bool):
if forward:
self._m_selfmus = copy(self.mus)
self._m_HFEparameterMap = copy(self.HFEngine.parameterMap)
self._mus = self.checkParameterListPivot(
self.mus(self.directionPivot))
self.HFEngine.parameterMap = {
"F": [self.HFEngine.parameterMap["F"][self.directionPivot[0]]],
"B": [self.HFEngine.parameterMap["B"][self.directionPivot[0]]]}
else:
self._mus = self._m_selfmus
self.HFEngine.parameterMap = self._m_HFEparameterMap
del self._m_selfmus, self._m_HFEparameterMap
def _marginalizeTrainedModel(self, forward:bool):
if forward:
del self._temporaryPivot
self.trainedModel.data.mu0 = self.mu0
self.trainedModel.data.scaleFactor = [1.] * self.npar
self.trainedModel.data.scaleFactor[self.directionPivot[0]] = (
self.scaleFactor[0])
self.trainedModel.data.parameterMap = self.HFEngine.parameterMap
self._m_musUniqueCN = copy(self._musUniqueCN)
musUniqueCNAux = np.zeros((self.S, self.npar), dtype = np.complex)
musUniqueCNAux[:, self.directionPivot[0]] = self._musUniqueCN(0)
self._musUniqueCN = self.checkParameterList(musUniqueCNAux)
self._m_derIdxs = copy(self._derIdxs)
for j in range(len(self._derIdxs)):
for l in range(len(self._derIdxs[j])):
derjl = self._derIdxs[j][l][0]
self._derIdxs[j][l] = [0] * self.npar
self._derIdxs[j][l][self.directionPivot[0]] = derjl
self.trainedModel.data.Q._dirPivot = self.directionPivot[0]
self.trainedModel.data.P._dirPivot = self.directionPivot[0]
# tell greedy error estimator that operator / RHS is pivot-affine
if hasattr(self.HFEngine.A, "is_affine"):
self._A_is_affine = self.HFEngine.A.is_affine
else:
self._A_is_affine = 0
if hasattr(self.HFEngine.b, "is_affine"):
self._b_is_affine = self.HFEngine.b.is_affine
else:
self._b_is_affine = 0
if self._A_is_affine >= 1 / 2 and self._b_is_affine >= 1 / 2:
self._affine_lvl += [1 / 2]
else:
self._temporaryPivot = 1
self.trainedModel.data.mu0 = self.checkParameterListPivot(
self.mu0(self.directionPivot))
self.trainedModel.data.scaleFactor = self.scaleFactor
self.trainedModel.data.parameterMap = {
"F": [self.HFEngine.parameterMap["F"][self.directionPivot[0]]],
"B": [self.HFEngine.parameterMap["B"][self.directionPivot[0]]]}
self._musUniqueCN = copy(self._m_musUniqueCN)
self._derIdxs = copy(self._m_derIdxs)
del self._m_musUniqueCN, self._m_derIdxs
del self.trainedModel.data.Q._dirPivot
del self.trainedModel.data.P._dirPivot
if self._A_is_affine >= 1 / 2 and self._b_is_affine >= 1 / 2:
self._affine_lvl.pop()
del self._A_is_affine, self._b_is_affine
self.trainedModel.data.npar = self.npar
def errorEstimator(self, mus:Np1D, return_max : bool = False) -> Np1D:
"""Standard residual-based error estimator."""
setupOK = self.setupApproxLocal()
if setupOK > 0:
err = np.empty(len(mus))
err[:] = np.nan
if not return_max: return err
return err, [- setupOK], np.nan
self._marginalizeTrainedModel(True)
errRes = super().errorEstimator(mus, return_max)
self._marginalizeTrainedModel(False)
return errRes
def _preliminaryTraining(self):
"""Initialize starting snapshots of solution map."""
RROMPyAssert(self._mode, message = "Cannot start greedy algorithm.")
self.resetSamples()
self.samplingEngine.scaleFactor = self.scaleFactorDer
musPivot = self.samplerTrainSet.generatePoints(self.S)
while len(musPivot) > self.S: musPivot.pop()
muTestPivot = self.samplerPivot.generatePoints(self.nTestPoints, False)
idxPop = pruneSamples(self.mapParameterListPivot(muTestPivot),
self.mapParameterListPivot(musPivot),
1e-10 * self.scaleFactorPivot[0])
muTestPivot.pop(idxPop)
self._mus = emptyParameterList()
self.mus.reset((self.S - 1, self.HFEngine.npar))
self.muTest = emptyParameterList()
self.muTest.reset((len(muTestPivot) + 1, self.HFEngine.npar))
self.mus.data[:, self.directionPivot] = musPivot[: -1]
self.mus.data[:, self.directionMarginal] = np.repeat(self.muMargLoc,
self.S - 1, axis = 0)
self.muTest.data[: -1, self.directionPivot] = muTestPivot.data
self.muTest.data[-1, self.directionPivot] = musPivot[-1]
self.muTest.data[:, self.directionMarginal] = np.repeat(self.muMargLoc,
len(muTestPivot) + 1,
axis = 0)
- if len(self.mus) > 0:
- vbMng(self, "MAIN",
- ("Adding first {} sample point{} at {} to training "
- "set.").format(self.S - 1, "" + "s" * (self.S > 2),
- self.mus), 3)
- self.samplingEngine.iterSample(self.mus)
+ vbMng(self, "MAIN",
+ ("Adding first {} sample point{} at {} to training "
+ "set.").format(self.S - 1, "" + "s" * (self.S > 2), self.mus),
+ 3)
+ self.samplingEngine.iterSample(self.mus)
self._S = len(self.mus)
self._approxParameters["S"] = self.S
def setupApproxLocal(self) -> int:
"""Compute rational interpolant."""
self._marginalizeMiscellanea(True)
setupOK = super().setupApproxLocal()
self._marginalizeMiscellanea(False)
return setupOK
def addMarginalSamplePoints(self, musMarginal:paramList, *args,
**kwargs) -> int:
"""Add marginal sample points to reduced model."""
RROMPyAssert(self._mode, message = "Cannot add sample points.")
musMarginal = self.checkParameterListMarginal(musMarginal)
vbMng(self, "INIT",
"Adding marginal sample point{} at {}.".format(
"s" * (len(musMarginal) > 1), musMarginal), 5)
if (self.SMarginal > 0 and hasattr(self, "polybasisMarginal")
and self.polybasisMarginal in sk):
RROMPyWarning(("Manually adding new samples with piecewise linear "
"marginal interpolation is dangerous. Sample depth "
"in samplerMarginal must be managed correctly."))
_musOld = self.mus
self._musMarginal.append(musMarginal)
S0 = copy(self.S)
+ req, is_master = [], masterCore()
idx, sizes = indicesScatter(len(musMarginal), return_sizes = True)
_trainedModelOld = copy(self.trainedModel)
_collapsed = (_trainedModelOld is not None
and _trainedModelOld.data._collapsed)
pMat, Ps, Qs, mus = None, [], [], None
- req, emptyCores = [], np.where(sizes == 0)[0]
+ emptyCores = np.where(np.array(sizes) == 0)[0] if is_master else []
if len(idx) == 0:
vbMng(self, "MAIN", "Idling.", 25)
if self.storeAllSamples: self.storeSamples()
pL = recv(source = 0, tag = poolRank())
pMat = np.empty((pL, 0), dtype = np.complex)
mus = np.empty((0, self.mu0.shape[1]), dtype = np.complex)
else:
_scaleFactorOldPivot = copy(self.scaleFactor)
self.scaleFactor = self.scaleFactorPivot
self._temporaryPivot = 1
for i in idx:
self.muMargLoc = self.musMarginal[[i + self.SMarginal]]
vbMng(self, "MAIN",
"Building marginal model no. {} at {}.".format(
i + self.SMarginal + 1,
self.musMarginal[i + self.SMarginal]), 5)
self.samplingEngine.resetHistory()
self.trainedModel = None
self.verbosity -= 5
self.samplingEngine.verbosity -= 5
RationalInterpolantGreedy.setupApprox(self, *args, **kwargs)
self.verbosity += 5
self.samplingEngine.verbosity += 5
musi = self.samplingEngine.mus
pMati = self.samplingEngine.projectionMatrix
if self.storeAllSamples: self.storeSamples(i + self.SMarginal)
if not self.matchState:
if self.POD == 1 and not (
hasattr(self.HFEngine.C, "is_mu_independent")
and self.HFEngine.C.is_mu_independent in self._output_lvl):
raise RROMPyException(("Cannot apply mu-dependent C "
"to orthonormalized samples."))
vbMng(self, "INIT", "Extracting system output from state.",
35)
pMatiEff = None
for j, mu in enumerate(musi):
pMij = np.expand_dims(self.HFEngine.applyC(
pMati[:, j], mu), -1)
if pMatiEff is None:
pMatiEff = np.array(pMij)
else:
pMatiEff = np.append(pMatiEff, pMij, axis = 1)
pMati = pMatiEff
vbMng(self, "DEL", "Done extracting system output.", 35)
if i == 0:
mus = copy(musi.data)
- if masterCore():
- for dest in emptyCores:
- req += [isend(len(pMati), dest = dest, tag = dest)]
+ for dest in emptyCores:
+ req += [isend(len(pMati), dest = dest, tag = dest)]
else:
mus = np.vstack((mus, musi.data))
if not _collapsed:
if pMat is None:
pMat = copy(pMati)
else:
pMat = np.hstack((pMat, pMati))
Ps += [copy(self.trainedModel.data.P)]
Qs += [copy(self.trainedModel.data.Q)]
if _collapsed: # collapse by hand
Ps[-1].postmultiplyTensorize(pMati.T)
self._S = S0
del self._temporaryPivot, self.muMargLoc
self.scaleFactor = _scaleFactorOldPivot
for r in req: r.wait()
if _collapsed: pMat = pMati[:, : 0]
pMat, Ps, Qs, mus, nsamples = gatherPivotedApproximant(pMat, Ps, Qs,
mus, sizes, self.polybasis)
self._mus = _musOld
self.mus.append(mus)
Psupp = np.append(0, np.cumsum(nsamples[: -1]).astype(int))
if _trainedModelOld is None:
self._setupTrainedModel(pMat, forceNew = True)
self.trainedModel.data.Qs, self.trainedModel.data.Ps = [], []
self.trainedModel.data.Psupp = []
else:
self._trainedModel = _trainedModelOld
if _collapsed:
self._setupTrainedModel(1.)
Psupp = [0] * len(musMarginal)
else:
Psupp = Psupp + self.trainedModel.data.projMat.shape[1]
self._setupTrainedModel(pMat, 1)
self._SMarginal += len(musMarginal)
self.trainedModel.data.Qs += Qs
self.trainedModel.data.Ps += Ps
self.trainedModel.data.Psupp += list(Psupp)
self._preliminaryMarginalFinalization()
self._finalizeMarginalization()
vbMng(self, "DEL", "Done adding marginal sample points.", 5)
return 0
def setupApprox(self, *args, **kwargs) -> int:
"""Compute rational interpolant."""
if self.checkComputedApprox(): return -1
RROMPyAssert(self._mode, message = "Cannot setup approximant.")
vbMng(self, "INIT", "Setting up {}.". format(self.name()), 5)
self.computeScaleFactor()
self._mus = emptyParameterList()
self._musMarginal = emptyParameterList()
musMarginal = self.samplerMarginal.generatePoints(self.SMarginal)
while len(musMarginal) > self.SMarginal: musMarginal.pop()
self._SMarginal = 0
val = self.addMarginalSamplePoints(musMarginal, *args, **kwargs)
vbMng(self, "DEL", "Done setting up approximant.", 5)
return val
class RationalInterpolantGreedyPivotedNoMatch(
RationalInterpolantGreedyPivotedBase,
GenericPivotedApproximantNoMatch):
"""
ROM pivoted rational interpolant (without matching) 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': kind of snapshots orthogonalization; allowed values
include 0, 1/2, and 1; defaults to 1, i.e. POD;
- 'scaleFactorDer': scaling factors for derivative computation;
defaults to 'AUTO';
- 'S': total number of pivot samples current approximant relies
upon;
- 'N': degree of rational interpolant denominator; defaults to
'AUTO', i.e. maximum allowed;
- 'samplerPivot': pivot sample point generator;
- 'SMarginal': total number of marginal samples current approximant
relies upon;
- 'samplerMarginal': marginal sample point generator;
- 'polybasis': type of polynomial basis for pivot
interpolation; defaults to 'MONOMIAL';
- 'greedyTol': uniform error tolerance for greedy algorithm;
defaults to 1e-2;
- 'collinearityTol': collinearity tolerance for greedy algorithm;
defaults to 0.;
- 'maxIter': maximum number of greedy steps; defaults to 1e2;
- 'nTestPoints': number of test points; defaults to 5e2;
- 'samplerTrainSet': training sample points generator; defaults to
samplerPivot;
- 'errorEstimatorKind': kind of error estimator; available values
include 'AFFINE', 'DISCREPANCY', 'LOOK_AHEAD',
'LOOK_AHEAD_RES', and 'NONE'; defaults to 'NONE';
- 'radialDirectionalWeightsMarginal': radial basis weights for
marginal interpolant; defaults to 1;
- 'functionalSolve': strategy for minimization of denominator
functional; allowed values include 'NORM', 'DOMINANT',
'BARYCENTRIC[_NORM]', and 'BARYCENTRIC_AVERAGE' (check pdf in
main folder for explanation); defaults to 'NORM';
- 'interpTol': tolerance for pivot interpolation; defaults to
None;
- 'forceQReal': force denominator to have real coefficients;
defaults to False;
- 'polyTruncateTol': tolerance for truncation of rational terms;
defaults to 0;
- 'QTol': tolerance for robust rational denominator management;
defaults to 0.
Defaults to empty dict.
verbosity(optional): Verbosity level. Defaults to 10.
Attributes:
HFEngine: HF problem solver.
mu0: Default parameter.
directionPivot: Pivot components.
mus: Array of snapshot parameters.
musMarginal: Array of marginal snapshot parameters.
approxParameters: Dictionary containing values for main parameters of
approximant. Recognized keys are in parameterList.
parameterListSoft: Recognized keys of soft approximant parameters:
- 'POD': kind of snapshots orthogonalization;
- 'scaleFactorDer': scaling factors for derivative computation;
- 'N': degree of rational interpolant denominator;
- 'polybasis': type of polynomial basis for pivot interpolation;
- 'greedyTol': uniform error tolerance for greedy algorithm;
- 'collinearityTol': collinearity tolerance for greedy algorithm;
- 'maxIter': maximum number of greedy steps;
- 'nTestPoints': number of test points;
- 'samplerTrainSet': training sample points generator;
- 'errorEstimatorKind': kind of error estimator;
- 'radialDirectionalWeightsMarginal': radial basis weights for
marginal interpolant;
- 'functionalSolve': strategy for minimization of denominator
functional;
- 'interpTol': tolerance for pivot interpolation;
- 'forceQReal': force denominator to have real coefficients;
- 'polyTruncateTol': tolerance for truncation of rational terms;
- 'QTol': 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.
verbosity: Verbosity level.
POD: Kind of snapshots orthogonalization.
scaleFactorDer: Scaling factors for derivative computation.
S: Total number of pivot samples current approximant relies upon.
N: Denominator degree of approximant.
samplerPivot: Pivot sample point generator.
SMarginal: Total number of marginal samples current approximant relies
upon.
samplerMarginal: Marginal sample point generator.
polybasis: Type of polynomial basis for pivot interpolation.
greedyTol: uniform error tolerance for greedy algorithm.
collinearityTol: Collinearity tolerance for greedy algorithm.
maxIter: maximum number of greedy steps.
nTestPoints: number of starting training points.
samplerTrainSet: training sample points generator.
errorEstimatorKind: kind of error estimator.
radialDirectionalWeightsMarginal: Radial basis weights for marginal
interpolant.
functionalSolve: Strategy for minimization of denominator functional.
interpTol: Tolerance for pivot interpolation.
forceQReal: Force denominator to have real coefficients.
polyTruncateTol: Tolerance for truncation of rational terms.
QTol: Tolerance for robust rational denominator management.
muBounds: 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.
"""
class RationalInterpolantGreedyPivotedPolyMatch(
RationalInterpolantGreedyPivotedBase,
GenericPivotedApproximantPolyMatch):
"""
ROM pivoted rational interpolant (with polynomial matching) 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': kind of snapshots orthogonalization; allowed values
include 0, 1/2, and 1; defaults to 1, i.e. POD;
- 'scaleFactorDer': scaling factors for derivative computation;
defaults to 'AUTO';
- 'matchState': whether to match the system state rather than the
system output; defaults to False;
- 'matchingWeight': weight for matching; defaults to 1;
- 'matchingKind': kind of matching; allowed values include 'ROTATE'
and 'PROJECT'; defaults to 'ROTATE';
- 'S': total number of pivot samples current approximant relies
upon;
- 'N': degree of rational interpolant denominator; defaults to
'AUTO', i.e. maximum allowed;
- 'samplerPivot': pivot sample point generator;
- 'SMarginal': total number of marginal samples current approximant
relies upon;
- 'samplerMarginal': marginal sample point generator;
- 'polybasis': type of polynomial basis for pivot
interpolation; defaults to 'MONOMIAL';
- 'polybasisMarginal': type of polynomial basis for marginal
interpolation; allowed values include 'MONOMIAL_*',
'CHEBYSHEV_*', 'LEGENDRE_*', 'NEARESTNEIGHBOR', and
'PIECEWISE_LINEAR_*'; defaults to 'MONOMIAL';
- 'paramsMarginal': dictionary of parameters for marginal
interpolation; include:
. 'MMarginal': degree of marginal interpolant; defaults to
'AUTO', i.e. maximum allowed; not for 'NEARESTNEIGHBOR' or
'PIECEWISE_LINEAR_*';
. 'nNeighborsMarginal': number of marginal nearest neighbors;
defaults to 1; only for 'NEARESTNEIGHBOR';
. 'polydegreetypeMarginal': type of polynomial degree for
marginal; defaults to 'TOTAL'; not for 'NEARESTNEIGHBOR' or
'PIECEWISE_LINEAR_*';
. 'polyTruncateTolMarginal': tolerance for truncation of
marginal interpolator; defaults to 0;
. 'interpTolMarginal': tolerance for marginal interpolation;
defaults to None; not for 'NEARESTNEIGHBOR' or
'PIECEWISE_LINEAR_*';
. 'radialDirectionalWeightsMarginalAdapt': bounds for adaptive
rescaling of marginal radial basis weights; only for
radial basis.
- 'greedyTol': uniform error tolerance for greedy algorithm;
defaults to 1e-2;
- 'collinearityTol': collinearity tolerance for greedy algorithm;
defaults to 0.;
- 'maxIter': maximum number of greedy steps; defaults to 1e2;
- 'nTestPoints': number of test points; defaults to 5e2;
- 'samplerTrainSet': training sample points generator; defaults to
samplerPivot;
- 'errorEstimatorKind': kind of error estimator; available values
include 'AFFINE', 'DISCREPANCY', 'LOOK_AHEAD',
'LOOK_AHEAD_RES', and 'NONE'; defaults to 'NONE';
- 'radialDirectionalWeightsMarginal': radial basis weights for
marginal interpolant; defaults to 1;
- 'functionalSolve': strategy for minimization of denominator
functional; allowed values include 'NORM', 'DOMINANT',
'BARYCENTRIC[_NORM]', and 'BARYCENTRIC_AVERAGE' (check pdf in
main folder for explanation); defaults to 'NORM';
- 'interpTol': tolerance for pivot interpolation; defaults to None;
- 'forceQReal': force denominator to have real coefficients;
defaults to False;
- 'polyTruncateTol': tolerance for truncation of rational terms;
defaults to 0;
- 'QTol': tolerance for robust rational denominator management;
defaults to 0.
Defaults to empty dict.
verbosity(optional): Verbosity level. Defaults to 10.
Attributes:
HFEngine: HF problem solver.
mu0: Default parameter.
directionPivot: Pivot components.
mus: Array of snapshot parameters.
musMarginal: Array of marginal snapshot parameters.
approxParameters: Dictionary containing values for main parameters of
approximant. Recognized keys are in parameterList.
parameterListSoft: Recognized keys of soft approximant parameters:
- 'POD': kind of snapshots orthogonalization;
- 'scaleFactorDer': scaling factors for derivative computation;
- 'matchState': whether to match the system state rather than the
system output;
- 'matchingWeight': weight for matching;
- 'matchingKind': kind of matching;
- 'N': degree of rational interpolant denominator;
- 'polybasis': type of polynomial basis for pivot interpolation;
- 'polybasisMarginal': type of polynomial basis for marginal
interpolation;
- 'paramsMarginal': dictionary of parameters for marginal
interpolation; include:
. 'MMarginal': degree of marginal interpolant;
. 'nNeighborsMarginal': number of marginal nearest neighbors;
. 'polydegreetypeMarginal': type of polynomial degree for
marginal;
. 'polyTruncateTolMarginal': tolerance for truncation of
marginal interpolator;
. 'interpTolMarginal': tolerance for marginal interpolation;
. 'radialDirectionalWeightsMarginalAdapt': bounds for adaptive
rescaling of marginal radial basis weights.
- 'greedyTol': uniform error tolerance for greedy algorithm;
- 'collinearityTol': collinearity tolerance for greedy algorithm;
- 'maxIter': maximum number of greedy steps;
- 'nTestPoints': number of test points;
- 'samplerTrainSet': training sample points generator;
- 'errorEstimatorKind': kind of error estimator;
- 'radialDirectionalWeightsMarginal': radial basis weights for
marginal interpolant;
- 'functionalSolve': strategy for minimization of denominator
functional;
- 'interpTol': tolerance for pivot interpolation;
- 'forceQReal': force denominator to have real coefficients;
- 'polyTruncateTol': tolerance for truncation of rational terms;
- 'QTol': 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.
verbosity: Verbosity level.
POD: Kind of snapshots orthogonalization.
scaleFactorDer: Scaling factors for derivative computation.
matchState: Whether to match the system state rather than the system
output.
matchingWeight: Weight for matching.
matchingKind: Kind of matching.
S: Total number of pivot samples current approximant relies upon.
N: Denominator degree of approximant.
samplerPivot: Pivot sample point generator.
SMarginal: Total number of marginal samples current approximant relies
upon.
samplerMarginal: Marginal sample point generator.
polybasis: Type of polynomial basis for pivot interpolation.
polybasisMarginal: Type of polynomial basis for marginal interpolation.
paramsMarginal: Dictionary of parameters for marginal interpolation.
greedyTol: uniform error tolerance for greedy algorithm.
collinearityTol: Collinearity tolerance for greedy algorithm.
maxIter: maximum number of greedy steps.
nTestPoints: number of starting training points.
samplerTrainSet: training sample points generator.
errorEstimatorKind: kind of error estimator.
radialDirectionalWeightsMarginal: Radial basis weights for marginal
interpolant.
functionalSolve: Strategy for minimization of denominator functional.
interpTol: Tolerance for pivot interpolation.
forceQReal: Force denominator to have real coefficients.
polyTruncateTol: Tolerance for truncation of rational terms.
QTol: Tolerance for robust rational denominator management.
muBounds: 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, *args, **kwargs) -> int:
if self.checkComputedApprox(): return -1
self.purgeparamsMarginal()
setupOK = super().setupApprox(*args, **kwargs)
if self.matchState: self._postApplyC()
return setupOK
class RationalInterpolantGreedyPivotedPoleMatch(
RationalInterpolantGreedyPivotedBase,
GenericPivotedApproximantPoleMatch):
"""
ROM pivoted rational interpolant (with pole matching) 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': kind of snapshots orthogonalization; allowed values
include 0, 1/2, and 1; defaults to 1, i.e. POD;
- 'scaleFactorDer': scaling factors for derivative computation;
defaults to 'AUTO';
- 'matchState': whether to match the system state rather than the
system output; defaults to False;
- 'matchingWeight': weight for pole matching optimization; defaults
to 1;
- 'matchingShared': required ratio of marginal points to share
resonance; defaults to 1.;
- 'badPoleCorrection': strategy for correction of bad poles;
available values include 'ERASE', 'RATIONAL', and 'POLYNOMIAL';
defaults to 'ERASE';
- '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;
- 'polybasis': type of polynomial basis for pivot
interpolation; defaults to 'MONOMIAL';
- 'polybasisMarginal': type of polynomial basis for marginal
interpolation; allowed values include 'MONOMIAL_*',
'CHEBYSHEV_*', 'LEGENDRE_*', 'NEARESTNEIGHBOR', and
'PIECEWISE_LINEAR_*'; defaults to 'MONOMIAL';
- 'paramsMarginal': dictionary of parameters for marginal
interpolation; include:
. 'MMarginal': degree of marginal interpolant; defaults to
'AUTO', i.e. maximum allowed; not for 'NEARESTNEIGHBOR' or
'PIECEWISE_LINEAR_*';
. 'nNeighborsMarginal': number of marginal nearest neighbors;
defaults to 1; only for 'NEARESTNEIGHBOR';
. 'polydegreetypeMarginal': type of polynomial degree for
marginal; defaults to 'TOTAL'; not for 'NEARESTNEIGHBOR' or
'PIECEWISE_LINEAR_*';
. 'polyTruncateTolMarginal': tolerance for truncation of
marginal interpolator; defaults to 0;
. 'interpTolMarginal': tolerance for marginal interpolation;
defaults to None; not for 'NEARESTNEIGHBOR' or
'PIECEWISE_LINEAR_*';
. 'radialDirectionalWeightsMarginalAdapt': bounds for adaptive
rescaling of marginal radial basis weights; only for
radial basis.
- 'greedyTol': uniform error tolerance for greedy algorithm;
defaults to 1e-2;
- 'collinearityTol': collinearity tolerance for greedy algorithm;
defaults to 0.;
- 'maxIter': maximum number of greedy steps; defaults to 1e2;
- 'nTestPoints': number of test points; defaults to 5e2;
- 'samplerTrainSet': training sample points generator; defaults to
samplerPivot;
- 'errorEstimatorKind': kind of error estimator; available values
include 'AFFINE', 'DISCREPANCY', 'LOOK_AHEAD',
'LOOK_AHEAD_RES', and 'NONE'; defaults to 'NONE';
- 'radialDirectionalWeightsMarginal': radial basis weights for
marginal interpolant; defaults to 1;
- 'functionalSolve': strategy for minimization of denominator
functional; allowed values include 'NORM', 'DOMINANT',
'BARYCENTRIC[_NORM]', and 'BARYCENTRIC_AVERAGE' (check pdf in
main folder for explanation); defaults to 'NORM';
- 'interpTol': tolerance for pivot interpolation; defaults to None;
- 'forceQReal': force denominator to have real coefficients;
defaults to False;
- 'polyTruncateTol': tolerance for truncation of rational terms;
defaults to 0;
- 'QTol': tolerance for robust rational denominator management;
defaults to 0.
Defaults to empty dict.
verbosity(optional): Verbosity level. Defaults to 10.
Attributes:
HFEngine: HF problem solver.
mu0: Default parameter.
directionPivot: Pivot components.
mus: Array of snapshot parameters.
musMarginal: Array of marginal snapshot parameters.
approxParameters: Dictionary containing values for main parameters of
approximant. Recognized keys are in parameterList.
parameterListSoft: Recognized keys of soft approximant parameters:
- 'POD': kind of snapshots orthogonalization;
- 'scaleFactorDer': scaling factors for derivative computation;
- 'matchState': whether to match the system state rather than the
system output;
- 'matchingWeight': weight for pole matching optimization;
- 'matchingShared': required ratio of marginal points to share
resonance;
- 'badPoleCorrection': strategy for correction of bad poles;
- 'polybasis': type of polynomial basis for pivot
interpolation;
- 'polybasisMarginal': type of polynomial basis for marginal
interpolation;
- 'paramsMarginal': dictionary of parameters for marginal
interpolation; include:
. 'MMarginal': degree of marginal interpolant;
. 'nNeighborsMarginal': number of marginal nearest neighbors;
. 'polydegreetypeMarginal': type of polynomial degree for
marginal;
. 'polyTruncateTolMarginal': tolerance for truncation of
marginal interpolator;
. 'interpTolMarginal': tolerance for marginal interpolation;
. 'radialDirectionalWeightsMarginalAdapt': bounds for adaptive
rescaling of marginal radial basis weights.
- 'greedyTol': uniform error tolerance for greedy algorithm;
- 'collinearityTol': collinearity tolerance for greedy algorithm;
- 'maxIter': maximum number of greedy steps;
- 'nTestPoints': number of test points;
- 'samplerTrainSet': training sample points generator;
- 'errorEstimatorKind': kind of error estimator;
- 'radialDirectionalWeightsMarginal': radial basis weights for
marginal interpolant;
- 'functionalSolve': strategy for minimization of denominator
functional;
- 'interpTol': tolerance for pivot interpolation;
- 'forceQReal': force denominator to have real coefficients;
- 'polyTruncateTol': tolerance for truncation of rational terms;
- 'QTol': 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.
verbosity: Verbosity level.
POD: Kind of snapshots orthogonalization.
scaleFactorDer: Scaling factors for derivative computation.
matchState: Whether to match the system state rather than the system
output.
matchingWeight: Weight for pole matching optimization.
matchingShared: Required ratio of marginal points to share resonance.
badPoleCorrection: Strategy for correction of bad 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.
polybasis: Type of polynomial basis for pivot interpolation.
polybasisMarginal: Type of polynomial basis for marginal interpolation.
paramsMarginal: Dictionary of parameters for marginal interpolation.
greedyTol: uniform error tolerance for greedy algorithm.
collinearityTol: Collinearity tolerance for greedy algorithm.
maxIter: maximum number of greedy steps.
nTestPoints: number of starting training points.
samplerTrainSet: training sample points generator.
errorEstimatorKind: kind of error estimator.
radialDirectionalWeightsMarginal: Radial basis weights for marginal
interpolant.
functionalSolve: Strategy for minimization of denominator functional.
interpTol: Tolerance for pivot interpolation.
forceQReal: Force denominator to have real coefficients.
polyTruncateTol: Tolerance for truncation of rational terms.
QTol: Tolerance for robust rational denominator management.
muBounds: 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, *args, **kwargs) -> int:
if self.checkComputedApprox(): return -1
self.purgeparamsMarginal()
setupOK = super().setupApprox(*args, **kwargs)
if self.matchState: self._postApplyC()
return setupOK
diff --git a/rrompy/reduction_methods/pivoted/rational_interpolant_pivoted.py b/rrompy/reduction_methods/pivoted/rational_interpolant_pivoted.py
index 613f190..42131a8 100755
--- a/rrompy/reduction_methods/pivoted/rational_interpolant_pivoted.py
+++ b/rrompy/reduction_methods/pivoted/rational_interpolant_pivoted.py
@@ -1,727 +1,724 @@
# Copyright (C) 2018-2020 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 collections.abc import Iterable
from copy import deepcopy as copy
from .generic_pivoted_approximant import (GenericPivotedApproximantBase,
GenericPivotedApproximantNoMatch,
GenericPivotedApproximantPolyMatch,
GenericPivotedApproximantPoleMatch)
from .gather_pivoted_approximant import gatherPivotedApproximant
from rrompy.reduction_methods.standard.rational_interpolant import (
RationalInterpolant)
from rrompy.utilities.base import verbosityManager as vbMng
from rrompy.utilities.base.types import paramList
from rrompy.utilities.numerical.hash_derivative import nextDerivativeIndices
from rrompy.utilities.poly_fitting.piecewise_linear import sparsekinds as sk
from rrompy.utilities.exception_manager import (RROMPyException, RROMPyAssert,
RROMPyWarning)
from rrompy.parameter import emptyParameterList
from rrompy.utilities.parallel import (masterCore, poolRank, indicesScatter,
isend, recv)
__all__ = ['RationalInterpolantPivotedNoMatch',
'RationalInterpolantPivotedPolyMatch',
'RationalInterpolantPivotedPoleMatch']
class RationalInterpolantPivotedBase(GenericPivotedApproximantBase,
RationalInterpolant):
def __init__(self, *args, **kwargs):
self._preInit()
self._addParametersToList(toBeExcluded = ["polydegreetype",
- "indexPrimary", "MAuxiliary",
- "NAuxiliary"])
+ "MAuxiliary", "NAuxiliary"])
super().__init__(*args, **kwargs)
if self.nparPivot > 1: self.HFEngine._ignoreResidues = 1
self._postInit()
- @property
- def indexPrimary(self): return 0
-
@property
def MAuxiliary(self): return 0
@property
def NAuxiliary(self): return 0
@property
- def polydegreetype(self): return "TENSOR"
+ def polydegreetype(self): return "TENSOR_TOTAL"
@property
def scaleFactorDer(self):
"""Value of scaleFactorDer."""
if self._scaleFactorDer == "NONE": return 1.
if self._scaleFactorDer == "AUTO": return self.scaleFactorPivot
return self._scaleFactorDer
@scaleFactorDer.setter
def scaleFactorDer(self, scaleFactorDer):
if isinstance(scaleFactorDer, (str,)):
scaleFactorDer = scaleFactorDer.upper()
elif isinstance(scaleFactorDer, Iterable):
scaleFactorDer = list(scaleFactorDer)
self._scaleFactorDer = scaleFactorDer
self._approxParameters["scaleFactorDer"] = self._scaleFactorDer
def _setupInterpolationIndices(self):
"""Setup parameters for polyvander."""
RROMPyAssert(self._mode,
message = "Cannot setup interpolation indices.")
if (self._musUniqueCN is None
or len(self._reorder) != len(self.musPivot)):
try:
muPC = self.trainedModel.centerNormalizePivot(self.musPivot)
except:
muPC = self.trainedModel.centerNormalize(self.musPivot)
self._musUniqueCN, musIdxsTo, musIdxs, musCount = (muPC.unique(
return_index = True, return_inverse = True,
return_counts = True))
self._musUnique = self.musPivot[musIdxsTo]
self._derIdxs = [None] * len(self._musUniqueCN)
self._reorder = np.empty(len(musIdxs), dtype = int)
filled = 0
for j, cnt in enumerate(musCount):
self._derIdxs[j] = nextDerivativeIndices([], self.nparPivot,
cnt)
jIdx = np.nonzero(musIdxs == j)[0]
self._reorder[jIdx] = np.arange(filled, filled + cnt)
filled += cnt
def addMarginalSamplePoints(self, musMarginal:paramList) -> int:
"""Add marginal sample points to reduced model."""
RROMPyAssert(self._mode, message = "Cannot add sample points.")
musMarginal = self.checkParameterListMarginal(musMarginal)
vbMng(self, "INIT",
"Adding marginal sample point{} at {}.".format(
"s" * (len(musMarginal) > 1), musMarginal), 5)
if (self.SMarginal > 0 and hasattr(self, "polybasisMarginal")
and self.polybasisMarginal in sk):
RROMPyWarning(("Manually adding new samples with piecewise linear "
"marginal interpolation is dangerous. Sample depth "
"in samplerMarginal must be managed correctly."))
mus = np.empty((self.S * len(musMarginal), self.HFEngine.npar),
dtype = np.complex)
mus[:, self.directionPivot] = np.tile(self.musPivot.data,
(len(musMarginal), 1))
mus[:, self.directionMarginal] = np.repeat(musMarginal.data, self.S,
axis = 0)
self._mus.append(mus)
self._musMarginal.append(musMarginal)
N0 = copy(self.N)
+ req, is_master = [], masterCore()
idx, sizes = indicesScatter(len(musMarginal), return_sizes = True)
_collapsed = self.trainedModel.data._collapsed
pMat, Ps, Qs = None, [], []
- req, emptyCores = [], np.where(sizes == 0)[0]
+ emptyCores = np.where(np.array(sizes) == 0)[0] if is_master else []
if len(idx) == 0:
vbMng(self, "MAIN", "Idling.", 30)
if self.storeAllSamples: self.storeSamples()
pL = recv(source = 0, tag = poolRank())
pMat = np.empty((pL, 0), dtype = np.complex)
else:
_scaleFactorOldPivot = copy(self.scaleFactor)
self.scaleFactor = self.scaleFactorPivot
self._temporaryPivot = 1
for i in idx:
musi = self.mus[self.S * (i + self.SMarginal)
: self.S * (i + self.SMarginal + 1)]
vbMng(self, "MAIN",
"Building marginal model no. {} at {}.".format(
i + self.SMarginal + 1,
self.musMarginal[i + self.SMarginal]), 5)
vbMng(self, "INIT", "Starting computation of snapshots.", 10)
self.samplingEngine.resetHistory()
self.samplingEngine.iterSample(musi)
vbMng(self, "DEL", "Done computing snapshots.", 10)
self.verbosity -= 5
self.samplingEngine.verbosity -= 5
if self.rationalMode[-6:] == "OUTPUT":
vbMng(self, "INIT", "Extracting system output from state.",
35)
self.samplingEngine.samples_output = self.HFEngine.applyC(
self.samplingEngine.samples, self.mus)
_POD, self._POD = self.POD, 1
vbMng(self, "DEL", "Done extracting system output.", 35)
self._setupRational(self._setupDenominator())
if self.rationalMode[-6:] == "OUTPUT":
self._POD = _POD
pMati = np.empty((1, 0), dtype = np.complex)
else:
pMati = self.samplingEngine.projectionMatrix
self.verbosity += 5
self.samplingEngine.verbosity += 5
if self.storeAllSamples: self.storeSamples(i + self.SMarginal)
if self.rationalMode[-6:] != "OUTPUT" and not self.matchState:
if self.POD == 1 and not (
hasattr(self.HFEngine.C, "is_mu_independent")
and self.HFEngine.C.is_mu_independent in self._output_lvl):
raise RROMPyException(("Cannot apply mu-dependent C "
"to orthonormalized samples."))
vbMng(self, "INIT", "Extracting system output from state.",
35)
pMatiEff = None
for j, mu in enumerate(musi):
pMij = np.expand_dims(self.HFEngine.applyC(
pMati[:, j], mu), -1)
if pMatiEff is None:
pMatiEff = np.array(pMij)
else:
pMatiEff = np.append(pMatiEff, pMij, axis = 1)
pMati = pMatiEff
vbMng(self, "DEL", "Done extracting system output.", 35)
- if i == 0 and masterCore():
+ if i == 0:
for dest in emptyCores:
req += [isend(len(pMati), dest = dest, tag = dest)]
if not _collapsed:
if pMat is None:
pMat = copy(pMati)
else:
pMat = np.hstack((pMat, pMati))
Ps += [copy(self.trainedModel.data.P)]
Qs += [copy(self.trainedModel.data.Q)]
if _collapsed and self.rationalMode[-6:] != "OUTPUT":
Ps[-1].postmultiplyTensorize(pMati.T) # collapse by hand
del self.trainedModel.data.Q, self.trainedModel.data.P
self.N = N0
del self._temporaryPivot
self.scaleFactor = _scaleFactorOldPivot
if _collapsed: pMat = pMati[:, : 0]
for r in req: r.wait()
pMat, Ps, Qs, _, _ = gatherPivotedApproximant(pMat, Ps, Qs,
self.mus.data, sizes,
self.polybasis, False)
if _collapsed:
self._setupTrainedModel(1.)
Psupp = [0] * len(musMarginal)
else:
self._setupTrainedModel(pMat,
len(self.trainedModel.data.projMat) > 0)
Psupp = (self.SMarginal + np.arange(0, len(musMarginal))) * self.S
self._SMarginal += len(musMarginal)
self.trainedModel.data.Qs += Qs
self.trainedModel.data.Ps += Ps
self.trainedModel.data.Psupp += list(Psupp)
self._preliminaryMarginalFinalization()
self._finalizeMarginalization()
vbMng(self, "DEL", "Done adding marginal sample points.", 5)
return 0
def setupApprox(self) -> int:
"""Compute rational interpolant."""
if self.checkComputedApprox(): return -1
RROMPyAssert(self._mode, message = "Cannot setup approximant.")
vbMng(self, "INIT", "Setting up {}.". format(self.name()), 5)
self.computeScaleFactor()
self.resetSamples()
self.samplingEngine.scaleFactor = self.scaleFactorDer
self._mus = emptyParameterList()
self._musMarginal = emptyParameterList()
self.musPivot = self.samplerPivot.generatePoints(self.S)
while len(self.musPivot) > self.S: self.musPivot.pop()
musMarginal = self.samplerMarginal.generatePoints(self.SMarginal)
while len(musMarginal) > self.SMarginal: musMarginal.pop()
self._setupTrainedModel(np.zeros((0, 0)), forceNew = True,
collapsed = self.rationalMode[-6:] == "OUTPUT")
self.trainedModel.data.Qs, self.trainedModel.data.Ps = [], []
self.trainedModel.data.Psupp = []
self._SMarginal = 0
val = self.addMarginalSamplePoints(musMarginal)
vbMng(self, "DEL", "Done setting up approximant.", 5)
return val
class RationalInterpolantPivotedNoMatch(RationalInterpolantPivotedBase,
GenericPivotedApproximantNoMatch):
"""
ROM pivoted rational interpolant (without matching) 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': kind of snapshots orthogonalization; allowed values
include 0, 1/2, and 1; defaults to 1, i.e. POD;
- 'scaleFactorDer': scaling factors for derivative computation;
defaults to 'AUTO';
- '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;
- 'rationalMode': mode of rational approximation; allowed values
include 'MINIMAL[_STATE]', 'MINIMAL_OUTPUT',
'STANDARD[_STATE]', and 'STANDARD_OUTPUT'; defaults to
'MINIMAL';
- 'polybasis': type of polynomial basis for pivot
interpolation; defaults to 'MONOMIAL';
- 'M': degree of rational interpolant numerator; defaults to
'AUTO', i.e. maximum allowed;
- 'N': degree of rational interpolant denominator; defaults to
'AUTO', i.e. maximum allowed;
- 'radialDirectionalWeightsMarginal': radial basis weights for
marginal interpolant; defaults to 1;
- 'functionalSolve': strategy for minimization of denominator
functional; allowed values include 'NORM', 'DOMINANT',
'BARYCENTRIC[_NORM]', and 'BARYCENTRIC_AVERAGE' (check pdf in
main folder for explanation); defaults to 'NORM';
- 'interpTol': tolerance for pivot interpolation; defaults to
None;
- 'forceQReal': force denominator to have real coefficients;
defaults to False;
- 'polyTruncateTol': tolerance for truncation of rational terms;
defaults to 0;
- 'QTol': tolerance for robust rational denominator management;
defaults to 0.
Defaults to empty dict.
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.
approxParameters: Dictionary containing values for main parameters of
approximant. Recognized keys are in parameterList.
parameterListSoft: Recognized keys of soft approximant parameters:
- 'POD': kind of snapshots orthogonalization;
- 'scaleFactorDer': scaling factors for derivative computation;
- 'rationalMode': mode of rational approximation;
- 'polybasis': type of polynomial basis for pivot
interpolation;
- 'M': degree of rational interpolant numerator;
- 'N': degree of rational interpolant denominator;
- 'radialDirectionalWeightsMarginal': radial basis weights for
marginal interpolant;
- 'functionalSolve': strategy for minimization of denominator
functional;
- 'interpTol': tolerance for pivot interpolation;
- 'forceQReal': force denominator to have real coefficients;
- 'polyTruncateTol': tolerance for truncation of rational terms;
- 'QTol': 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.
verbosity: Verbosity level.
POD: Kind of snapshots orthogonalization.
scaleFactorDer: Scaling factors for derivative computation.
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.
rationalMode: Mode of rational approximation.
polybasis: Type of polynomial basis for pivot interpolation.
M: Numerator degree of approximant.
N: Denominator degree of approximant.
radialDirectionalWeightsMarginal: Radial basis weights for marginal
interpolant.
functionalSolve: Strategy for minimization of denominator functional.
interpTol: Tolerance for pivot interpolation.
forceQReal: Force denominator to have real coefficients.
polyTruncateTol: Tolerance for truncation of rational terms.
QTol: Tolerance for robust rational denominator management.
muBounds: 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.
"""
class RationalInterpolantPivotedPolyMatch(RationalInterpolantPivotedBase,
GenericPivotedApproximantPolyMatch):
"""
ROM pivoted rational interpolant (with polynomial matching) 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': kind of snapshots orthogonalization; allowed values
include 0, 1/2, and 1; defaults to 1, i.e. POD;
- 'scaleFactorDer': scaling factors for derivative computation;
defaults to 'AUTO';
- 'matchState': whether to match the system state rather than the
system output; defaults to False;
- 'matchingWeight': weight for matching; defaults to 1;
- 'matchingKind': kind of matching; allowed values include 'ROTATE'
and 'PROJECT'; defaults to 'ROTATE';
- '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;
- 'rationalMode': mode of rational approximation; allowed values
include 'MINIMAL[_STATE]', 'MINIMAL_OUTPUT',
'STANDARD[_STATE]', and 'STANDARD_OUTPUT'; defaults to
'MINIMAL';
- 'polybasis': type of polynomial basis for pivot
interpolation; defaults to 'MONOMIAL';
- 'polybasisMarginal': type of polynomial basis for marginal
interpolation; allowed values include 'MONOMIAL_*',
'CHEBYSHEV_*', 'LEGENDRE_*', 'NEARESTNEIGHBOR', and
'PIECEWISE_LINEAR_*'; defaults to 'MONOMIAL';
- 'paramsMarginal': dictionary of parameters for marginal
interpolation; include:
. 'MMarginal': degree of marginal interpolant; defaults to
'AUTO', i.e. maximum allowed; not for 'NEARESTNEIGHBOR' or
'PIECEWISE_LINEAR_*';
. 'nNeighborsMarginal': number of marginal nearest neighbors;
defaults to 1; only for 'NEARESTNEIGHBOR';
. 'polydegreetypeMarginal': type of polynomial degree for
marginal; defaults to 'TOTAL'; not for 'NEARESTNEIGHBOR' or
'PIECEWISE_LINEAR_*';
. 'polyTruncateTolMarginal': tolerance for truncation of
marginal interpolator; defaults to 0;
. 'interpTolMarginal': tolerance for marginal interpolation;
defaults to None; not for 'NEARESTNEIGHBOR' or
'PIECEWISE_LINEAR_*';
. 'radialDirectionalWeightsMarginalAdapt': bounds for adaptive
rescaling of marginal radial basis weights; only for
radial basis.
- 'M': degree of rational interpolant numerator; defaults to
'AUTO', i.e. maximum allowed;
- 'N': degree of rational interpolant denominator; defaults to
'AUTO', i.e. maximum allowed;
- 'radialDirectionalWeightsMarginal': radial basis weights for
marginal interpolant; defaults to 1;
- 'functionalSolve': strategy for minimization of denominator
functional; allowed values include 'NORM', 'DOMINANT',
'BARYCENTRIC[_NORM]', and 'BARYCENTRIC_AVERAGE' (check pdf in
main folder for explanation); defaults to 'NORM';
- 'interpTol': tolerance for pivot interpolation; defaults to None;
- 'forceQReal': force denominator to have real coefficients;
defaults to False;
- 'polyTruncateTol': tolerance for truncation of rational terms;
defaults to 0;
- 'QTol': tolerance for robust rational denominator management;
defaults to 0.
Defaults to empty dict.
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.
approxParameters: Dictionary containing values for main parameters of
approximant. Recognized keys are in parameterList.
parameterListSoft: Recognized keys of soft approximant parameters:
- 'POD': kind of snapshots orthogonalization;
- 'scaleFactorDer': scaling factors for derivative computation;
- 'matchState': whether to match the system state rather than the
system output;
- 'matchingWeight': weight for matching;
- 'matchingKind': kind of matching;
- 'rationalMode': mode of rational approximation;
- 'polybasis': type of polynomial basis for pivot interpolation;
- 'polybasisMarginal': type of polynomial basis for marginal
interpolation;
- 'paramsMarginal': dictionary of parameters for marginal
interpolation; include:
. 'MMarginal': degree of marginal interpolant;
. 'nNeighborsMarginal': number of marginal nearest neighbors;
. 'polydegreetypeMarginal': type of polynomial degree for
marginal;
. 'polyTruncateTolMarginal': tolerance for truncation of
marginal interpolator;
. 'interpTolMarginal': tolerance for marginal interpolation;
. 'radialDirectionalWeightsMarginalAdapt': bounds for adaptive
rescaling of marginal radial basis weights.
- 'M': degree of rational interpolant numerator;
- 'N': degree of rational interpolant denominator;
- 'radialDirectionalWeightsMarginal': radial basis weights for
marginal interpolant;
- 'functionalSolve': strategy for minimization of denominator
functional;
- 'interpTol': tolerance for pivot interpolation;
- 'forceQReal': force denominator to have real coefficients;
- 'polyTruncateTol': tolerance for truncation of rational terms;
- 'QTol': 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.
verbosity: Verbosity level.
POD: Kind of snapshots orthogonalization.
scaleFactorDer: Scaling factors for derivative computation.
matchState: Whether to match the system state rather than the system
output.
matchingWeight: Weight for matching.
matchingKind: Kind of matching.
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.
rationalMode: Mode of rational approximation.
polybasis: Type of polynomial basis for pivot interpolation.
polybasisMarginal: Type of polynomial basis for marginal interpolation.
paramsMarginal: Dictionary of parameters for marginal interpolation.
M: Numerator degree of approximant.
N: Denominator degree of approximant.
radialDirectionalWeightsMarginal: Radial basis weights for marginal
interpolant.
functionalSolve: Strategy for minimization of denominator functional.
interpTol: Tolerance for pivot interpolation.
forceQReal: Force denominator to have real coefficients.
polyTruncateTol: Tolerance for truncation of rational terms.
QTol: Tolerance for robust rational denominator management.
muBounds: 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, *args, **kwargs) -> int:
if self.checkComputedApprox(): return -1
if self.rationalMode[-6:] == "OUTPUT": self.matchState = False
self.purgeparamsMarginal()
setupOK = super().setupApprox(*args, **kwargs)
if self.matchState: self._postApplyC()
return setupOK
class RationalInterpolantPivotedPoleMatch(RationalInterpolantPivotedBase,
GenericPivotedApproximantPoleMatch):
"""
ROM pivoted rational interpolant (with pole matching) 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': kind of snapshots orthogonalization; allowed values
include 0, 1/2, and 1; defaults to 1, i.e. POD;
- 'scaleFactorDer': scaling factors for derivative computation;
defaults to 'AUTO';
- 'matchState': whether to match the system state rather than the
system output; defaults to False;
- 'matchingWeight': weight for pole matching optimization; defaults
to 1;
- 'matchingShared': required ratio of marginal points to share
resonance; defaults to 1.;
- 'badPoleCorrection': strategy for correction of bad poles;
available values include 'ERASE', 'RATIONAL', and 'POLYNOMIAL';
defaults to 'ERASE';
- '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;
- 'rationalMode': mode of rational approximation; allowed values
include 'MINIMAL[_STATE]', 'MINIMAL_OUTPUT',
'STANDARD[_STATE]', and 'STANDARD_OUTPUT'; defaults to
'MINIMAL';
- 'polybasis': type of polynomial basis for pivot
interpolation; defaults to 'MONOMIAL';
- 'polybasisMarginal': type of polynomial basis for marginal
interpolation; allowed values include 'MONOMIAL_*',
'CHEBYSHEV_*', 'LEGENDRE_*', 'NEARESTNEIGHBOR', and
'PIECEWISE_LINEAR_*'; defaults to 'MONOMIAL';
- 'paramsMarginal': dictionary of parameters for marginal
interpolation; include:
. 'MMarginal': degree of marginal interpolant; defaults to
'AUTO', i.e. maximum allowed; not for 'NEARESTNEIGHBOR' or
'PIECEWISE_LINEAR_*';
. 'nNeighborsMarginal': number of marginal nearest neighbors;
defaults to 1; only for 'NEARESTNEIGHBOR';
. 'polydegreetypeMarginal': type of polynomial degree for
marginal; defaults to 'TOTAL'; not for 'NEARESTNEIGHBOR' or
'PIECEWISE_LINEAR_*';
. 'polyTruncateTolMarginal': tolerance for truncation of
marginal interpolator; defaults to 0;
. 'interpTolMarginal': tolerance for marginal interpolation;
defaults to None; not for 'NEARESTNEIGHBOR' or
'PIECEWISE_LINEAR_*';
. 'radialDirectionalWeightsMarginalAdapt': bounds for adaptive
rescaling of marginal radial basis weights; only for
radial basis.
- 'M': degree of rational interpolant numerator; defaults to
'AUTO', i.e. maximum allowed;
- 'N': degree of rational interpolant denominator; defaults to
'AUTO', i.e. maximum allowed;
- 'radialDirectionalWeightsMarginal': radial basis weights for
marginal interpolant; defaults to 1;
- 'functionalSolve': strategy for minimization of denominator
functional; allowed values include 'NORM', 'DOMINANT',
'BARYCENTRIC[_NORM]', and 'BARYCENTRIC_AVERAGE' (check pdf in
main folder for explanation); defaults to 'NORM';
- 'interpTol': tolerance for pivot interpolation; defaults to None;
- 'forceQReal': force denominator to have real coefficients;
defaults to False;
- 'polyTruncateTol': tolerance for truncation of rational terms;
defaults to 0;
- 'QTol': tolerance for robust rational denominator management;
defaults to 0.
Defaults to empty dict.
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.
approxParameters: Dictionary containing values for main parameters of
approximant. Recognized keys are in parameterList.
parameterListSoft: Recognized keys of soft approximant parameters:
- 'POD': kind of snapshots orthogonalization;
- 'scaleFactorDer': scaling factors for derivative computation;
- 'matchState': whether to match the system state rather than the
system output;
- 'matchingWeight': weight for pole matching optimization;
- 'matchingShared': required ratio of marginal points to share
resonance;
- 'badPoleCorrection': strategy for correction of bad poles;
- 'rationalMode': mode of rational approximation;
- 'polybasis': type of polynomial basis for pivot
interpolation;
- 'polybasisMarginal': type of polynomial basis for marginal
interpolation;
- 'paramsMarginal': dictionary of parameters for marginal
interpolation; include:
. 'MMarginal': degree of marginal interpolant;
. 'nNeighborsMarginal': number of marginal nearest neighbors;
. 'polydegreetypeMarginal': type of polynomial degree for
marginal;
. 'polyTruncateTolMarginal': tolerance for truncation of
marginal interpolator;
. 'interpTolMarginal': tolerance for marginal interpolation;
. 'radialDirectionalWeightsMarginalAdapt': bounds for adaptive
rescaling of marginal radial basis weights.
- 'M': degree of rational interpolant numerator;
- 'N': degree of rational interpolant denominator;
- 'functionalSolve': strategy for minimization of denominator
functional;
- 'interpTol': tolerance for pivot interpolation;
- 'forceQReal': force denominator to have real coefficients;
- 'polyTruncateTol': tolerance for truncation of rational terms;
- 'QTol': 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.
verbosity: Verbosity level.
POD: Kind of snapshots orthogonalization.
scaleFactorDer: Scaling factors for derivative computation.
matchState: Whether to match the system state rather than the system
output.
matchingWeight: Weight for pole matching optimization.
matchingShared: Required ratio of marginal points to share resonance.
badPoleCorrection: Strategy for correction of bad 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.
rationalMode: Mode of rational approximation.
polybasis: Type of polynomial basis for pivot interpolation.
polybasisMarginal: Type of polynomial basis for marginal interpolation.
paramsMarginal: Dictionary of parameters for marginal interpolation.
M: Numerator degree of approximant.
N: Denominator degree of approximant.
radialDirectionalWeightsMarginal: Radial basis weights for marginal
interpolant.
functionalSolve: Strategy for minimization of denominator functional.
interpTol: Tolerance for pivot interpolation.
forceQReal: Force denominator to have real coefficients.
polyTruncateTol: Tolerance for truncation of rational terms.
QTol: Tolerance for robust rational denominator management.
muBounds: 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, *args, **kwargs) -> int:
if self.checkComputedApprox(): return -1
self.purgeparamsMarginal()
setupOK = super().setupApprox(*args, **kwargs)
if self.matchState: self._postApplyC()
return setupOK
diff --git a/rrompy/reduction_methods/standard/greedy/generic_greedy_approximant.py b/rrompy/reduction_methods/standard/greedy/generic_greedy_approximant.py
index 40496a9..78a6ab8 100644
--- a/rrompy/reduction_methods/standard/greedy/generic_greedy_approximant.py
+++ b/rrompy/reduction_methods/standard/greedy/generic_greedy_approximant.py
@@ -1,633 +1,639 @@
# Copyright (C) 2018-2020 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 copy import deepcopy as copy
import numpy as np
from matplotlib import pyplot as plt
from rrompy.hfengines.base.linear_affine_engine import checkIfAffine
from rrompy.reduction_methods.standard.generic_standard_approximant import (
GenericStandardApproximant)
from rrompy.utilities.base.types import (Np1D, Np2D, Tuple, List, paramVal,
paramList, sampList)
from rrompy.utilities.base import verbosityManager as vbMng
from rrompy.utilities.numerical import dot
from rrompy.utilities.expression import expressionEvaluator
from rrompy.utilities.exception_manager import (RROMPyException, RROMPyAssert,
RROMPyWarning)
from rrompy.sampling.sample_list import sampleList
from rrompy.parameter import emptyParameterList, parameterList
from rrompy.utilities.parallel import masterCore
__all__ = ['GenericGreedyApproximant']
def localL2Distance(mus:Np2D, badmus:Np2D) -> Np2D:
return np.linalg.norm(np.tile(mus[..., np.newaxis], [1, 1, len(badmus)])
- badmus[..., np.newaxis].T, axis = 1)
def pruneSamples(mus:paramList, badmus:paramList,
tol : float = 1e-8) -> Np1D:
"""Remove from mus all the elements which are too close to badmus."""
if isinstance(mus, (parameterList, sampleList)): mus = mus.data
if isinstance(badmus, (parameterList, sampleList)): badmus = badmus.data
if len(badmus) == 0: return np.arange(len(mus))
proximity = np.min(localL2Distance(mus, badmus), axis = 1)
return np.where(proximity <= tol)[0]
class GenericGreedyApproximant(GenericStandardApproximant):
"""
ROM greedy interpolant computation for parametric problems
(ABSTRACT).
Args:
HFEngine: HF problem solver.
mu0(optional): Default parameter. Defaults to 0.
approxParameters(optional): Dictionary containing values for main
parameters of approximant. Recognized keys are:
- 'POD': kind of snapshots orthogonalization; allowed values
include 0, 1/2, and 1; defaults to 1, i.e. POD;
- 'scaleFactorDer': scaling factors for derivative computation;
defaults to 'AUTO';
- 'S': number of starting training points;
- 'sampler': sample point generator;
- 'greedyTol': uniform error tolerance for greedy algorithm;
defaults to 1e-2;
- 'collinearityTol': collinearity tolerance for greedy algorithm;
defaults to 0.;
- 'maxIter': maximum number of greedy steps; defaults to 1e2;
- 'nTestPoints': number of test points; defaults to 5e2;
- 'samplerTrainSet': training sample points generator; defaults to
sampler.
Defaults to empty dict.
verbosity(optional): Verbosity level. Defaults to 10.
Attributes:
HFEngine: HF problem solver.
mu0: Default parameter.
mus: Array of snapshot parameters.
approxParameters: Dictionary containing values for main parameters of
approximant. Recognized keys are in parameterList.
parameterListSoft: Recognized keys of soft approximant parameters:
- 'POD': kind of snapshots orthogonalization;
- 'scaleFactorDer': scaling factors for derivative computation;
- 'greedyTol': uniform error tolerance for greedy algorithm;
- 'collinearityTol': collinearity tolerance for greedy algorithm;
- 'maxIter': maximum number of greedy steps;
- 'nTestPoints': number of test points;
- 'samplerTrainSet': training sample points generator.
parameterListCritical: Recognized keys of critical approximant
parameters:
- 'S': total number of samples current approximant relies upon;
- 'sampler': sample point generator.
verbosity: Verbosity level.
POD: Kind of snapshots orthogonalization.
scaleFactorDer: Scaling factors for derivative computation.
S: number of test points.
sampler: Sample point generator.
greedyTol: Uniform error tolerance for greedy algorithm.
collinearityTol: Collinearity tolerance for greedy algorithm.
maxIter: maximum number of greedy steps.
nTestPoints: number of starting training points.
samplerTrainSet: training sample points generator.
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.
"""
def __init__(self, *args, **kwargs):
self._preInit()
if not hasattr(self, "_affine_lvl"): self._affine_lvl = []
self._affine_lvl += [1]
self._addParametersToList(["greedyTol", "collinearityTol", "maxIter",
"nTestPoints", "samplerTrainSet"],
[1e-2, 0., 1e2, 5e2, "AUTO"])
super().__init__(*args, **kwargs)
self._postInit()
@property
def greedyTol(self):
"""Value of greedyTol."""
return self._greedyTol
@greedyTol.setter
def greedyTol(self, greedyTol):
if greedyTol < 0:
raise RROMPyException("greedyTol must be non-negative.")
if hasattr(self, "_greedyTol") and self.greedyTol is not None:
greedyTolold = self.greedyTol
else:
greedyTolold = -1
self._greedyTol = greedyTol
self._approxParameters["greedyTol"] = self.greedyTol
if greedyTolold != self.greedyTol:
self.resetSamples()
@property
def collinearityTol(self):
"""Value of collinearityTol."""
return self._collinearityTol
@collinearityTol.setter
def collinearityTol(self, collinearityTol):
if collinearityTol < 0:
raise RROMPyException("collinearityTol must be non-negative.")
if (hasattr(self, "_collinearityTol")
and self.collinearityTol is not None):
collinearityTolold = self.collinearityTol
else:
collinearityTolold = -1
self._collinearityTol = collinearityTol
self._approxParameters["collinearityTol"] = self.collinearityTol
if collinearityTolold != self.collinearityTol:
self.resetSamples()
@property
def maxIter(self):
"""Value of maxIter."""
return self._maxIter
@maxIter.setter
def maxIter(self, maxIter):
if maxIter <= 0: raise RROMPyException("maxIter must be positive.")
if hasattr(self, "_maxIter") and self.maxIter is not None:
maxIterold = self.maxIter
else:
maxIterold = -1
self._maxIter = maxIter
self._approxParameters["maxIter"] = self.maxIter
if maxIterold != self.maxIter:
self.resetSamples()
@property
def nTestPoints(self):
"""Value of nTestPoints."""
return self._nTestPoints
@nTestPoints.setter
def nTestPoints(self, nTestPoints):
if nTestPoints <= 0:
raise RROMPyException("nTestPoints must be positive.")
if not np.isclose(nTestPoints, np.int(nTestPoints)):
raise RROMPyException("nTestPoints must be an integer.")
nTestPoints = np.int(nTestPoints)
if hasattr(self, "_nTestPoints") and self.nTestPoints is not None:
nTestPointsold = self.nTestPoints
else:
nTestPointsold = -1
self._nTestPoints = nTestPoints
self._approxParameters["nTestPoints"] = self.nTestPoints
if nTestPointsold != self.nTestPoints:
self.resetSamples()
@property
def samplerTrainSet(self):
"""Value of samplerTrainSet."""
return self._samplerTrainSet
@samplerTrainSet.setter
def samplerTrainSet(self, samplerTrainSet):
if (isinstance(samplerTrainSet, (str,))
and samplerTrainSet.upper() == "AUTO"):
samplerTrainSet = self.sampler
if 'generatePoints' not in dir(samplerTrainSet):
raise RROMPyException("samplerTrainSet type not recognized.")
if (hasattr(self, '_samplerTrainSet')
and self.samplerTrainSet not in [None, "AUTO"]):
samplerTrainSetOld = self.samplerTrainSet
self._samplerTrainSet = copy(samplerTrainSet)
self._approxParameters["samplerTrainSet"] = self.samplerTrainSet
if (not 'samplerTrainSetOld' in locals()
or samplerTrainSetOld != self.samplerTrainSet):
self.resetSamples()
def resetSamples(self):
"""Reset samples."""
super().resetSamples()
self._mus = emptyParameterList()
def _affineResidualMatricesContraction(self, rb:Np2D, rA : Np2D = None) \
-> Tuple[Np1D, Np1D, Np1D]:
self.assembleReducedResidualBlocks(full = rA is not None)
# 'ij,jk,ik->k', resbb, radiusb, radiusb.conj()
ff = np.sum(self.trainedModel.data.resbb.dot(rb) * rb.conj(), axis = 0)
if rA is None: return ff
# 'ijk,jkl,il->l', resAb, radiusA, radiusb.conj()
Lf = np.sum(np.tensordot(self.trainedModel.data.resAb, rA, 2)
* rb.conj(), axis = 0)
# 'ijkl,klt,ijt->t', resAA, radiusA, radiusA.conj()
LL = np.sum(np.tensordot(self.trainedModel.data.resAA, rA, 2)
* rA.conj(), axis = (0, 1))
return ff, Lf, LL
def getErrorEstimatorAffine(self, mus:Np1D) -> Np1D:
"""Standard residual estimator."""
checkIfAffine(self.HFEngine, "apply affinity-based error estimator",
False, self._affine_lvl)
self.HFEngine.buildA()
self.HFEngine.buildb()
mus = self.checkParameterList(mus)
tMverb, self.trainedModel.verbosity = self.trainedModel.verbosity, 0
uApproxRs = self.getApproxReduced(mus).data
self.trainedModel.verbosity = tMverb
muTestEff = self.mapParameterList(mus)
radiusA = np.empty((len(self.HFEngine.thAs), len(mus)),
dtype = np.complex)
radiusb = np.empty((len(self.HFEngine.thbs), len(mus)),
dtype = np.complex)
for j, thA in enumerate(self.HFEngine.thAs):
radiusA[j] = expressionEvaluator(thA[0], muTestEff)
for j, thb in enumerate(self.HFEngine.thbs):
radiusb[j] = expressionEvaluator(thb[0], muTestEff)
radiusA = np.expand_dims(uApproxRs, 1) * radiusA
ff, Lf, LL = self._affineResidualMatricesContraction(radiusb, radiusA)
err = np.abs((LL - 2. * np.real(Lf) + ff) / ff) ** .5
return err
def errorEstimator(self, mus:Np1D, return_max : bool = False) -> Np1D:
setupOK = self.setupApproxLocal()
if setupOK > 0:
err = np.empty(len(mus))
err[:] = np.nan
if not return_max: return err
return err, [- setupOK], np.nan
mus = self.checkParameterList(mus)
vbMng(self.trainedModel, "INIT",
"Evaluating error estimator at mu = {}.".format(mus), 10)
err = self.getErrorEstimatorAffine(mus)
vbMng(self.trainedModel, "DEL", "Done evaluating error estimator.", 10)
if not return_max: return err
idxMaxEst = [np.argmax(err)]
return err, idxMaxEst, err[idxMaxEst]
def _isLastSampleCollinear(self) -> bool:
"""Check collinearity of last sample."""
if self.collinearityTol <= 0.: return False
if self.POD == 1:
reff = self.samplingEngine.Rscale[:, -1]
else:
RROMPyWarning(("Repeated orthogonalization of the samples for "
"collinearity check. Consider setting POD to "
"True."))
if not hasattr(self, "_PODEngine"):
from rrompy.sampling import PODEngine
self._PODEngine = PODEngine(self.HFEngine)
reff = self._PODEngine.generalizedQR(self.samplingEngine.samples,
only_R = True,
is_state = True)[:, -1]
cLevel = np.abs(reff[-1]) / np.linalg.norm(reff)
cLevel = np.inf if np.isclose(cLevel, 0., atol = 1e-15) else 1 / cLevel
vbMng(self, "MAIN", "Collinearity indicator {:.4e}.".format(cLevel), 3)
return cLevel > self.collinearityTol
- def plotEstimator(self, est:Np1D, idxMax:List[int], estMax:List[float]):
+ def plotEstimator(self, est:Np1D, idxMax:List[int], estMax:Np1D):
if (not (np.any(np.isnan(est)) or np.any(np.isinf(est)))
and masterCore()):
fig = plt.figure(figsize = plt.figaspect(1. / self.npar))
for jpar in range(self.npar):
ax = fig.add_subplot(1, self.npar, 1 + jpar)
musre = np.array(self.muTest.re.data)
errCP = copy(est)
idx = np.delete(np.arange(self.npar), jpar)
while len(musre) > 0:
if self.npar == 1:
currIdx = np.arange(len(musre))
else:
currIdx = np.where(np.isclose(np.sum(
np.abs(musre[:, idx] - musre[0, idx]), 1),
0., atol = 1e-15))[0]
ax.semilogy(musre[currIdx, jpar], errCP[currIdx], 'k',
linewidth = 1)
musre = np.delete(musre, currIdx, 0)
errCP = np.delete(errCP, currIdx)
ax.semilogy([self.muBounds.re(0, jpar),
self.muBounds.re(-1, jpar)],
[self.greedyTol] * 2, 'r--')
ax.semilogy(self.mus.re(jpar),
2. * self.greedyTol * np.ones(len(self.mus)), '*m')
if len(idxMax) > 0 and estMax is not None:
ax.semilogy(self.muTest.re(idxMax, jpar), estMax, 'xr')
ax.set_xlim(*list(self.sampler.lims.re(jpar)))
ax.grid()
plt.tight_layout()
plt.show()
def greedyNextSample(self, muidx:int, plotEst : str = "NONE")\
-> Tuple[Np1D, int, float, paramVal]:
"""Compute next greedy snapshot of solution map."""
RROMPyAssert(self._mode, message = "Cannot add greedy sample.")
- mus = copy(self.muTest[muidx])
+ mus = parameterList(self.muTest[muidx])
self.muTest.pop(muidx)
+ if len(mus) == 1:
+ msg = ("Adding sample point no. {} at {} to training "
+ "set.").format(len(self.mus) + 1, mus[0])
+ else:
+ msg = ("Adding sample points no. {}--{} at {} to training "
+ "set.").format(len(self.mus) + 1, len(self.mus) + len(mus),
+ mus)
+ vbMng(self, "MAIN", msg, 3)
for j, mu in enumerate(mus):
- vbMng(self, "MAIN",
- ("Adding sample point no. {} at {} to training "
- "set.").format(len(self.mus) + 1, mu), 3)
self.mus.append(mu)
self._S = len(self.mus)
self._approxParameters["S"] = self.S
if (self.samplingEngine.nsamples <= len(mus) - j - 1
- or not np.allclose(mu, self.samplingEngine.mus[j - len(mus)])):
+ or not np.allclose(mu.data,
+ self.samplingEngine.mus[j - len(mus)])):
self.samplingEngine.nextSample(mu)
if self._isLastSampleCollinear():
vbMng(self, "MAIN",
("Collinearity above tolerance detected. Starting "
"preemptive greedy loop termination."), 3)
self._collinearityFlag = 1
errorEstTest = np.empty(len(self.muTest))
errorEstTest[:] = np.nan
return errorEstTest, [-1], np.nan, np.nan
errorEstTest, muidx, maxErrorEst = self.errorEstimator(self.muTest,
True)
if plotEst == "ALL":
self.plotEstimator(errorEstTest, muidx, maxErrorEst)
return errorEstTest, muidx, maxErrorEst, self.muTest[muidx]
def _preliminaryTraining(self):
"""Initialize starting snapshots of solution map."""
RROMPyAssert(self._mode, message = "Cannot start greedy algorithm.")
if self.samplingEngine.nsamples > 0: return
self.resetSamples()
self.computeScaleFactor()
self.samplingEngine.scaleFactor = self.scaleFactorDer
self.mus = self.samplerTrainSet.generatePoints(self.S)
while len(self.mus) > self.S: self.mus.pop()
muTestBase = self.sampler.generatePoints(self.nTestPoints, False)
idxPop = pruneSamples(self.mapParameterList(muTestBase),
self.mapParameterList(self.mus),
1e-10 * self.scaleFactor[0])
muTestBase.pop(idxPop)
muLast = copy(self.mus[-1])
self.mus.pop()
if len(self.mus) > 0:
vbMng(self, "MAIN",
("Adding first {} sample point{} at {} to training "
"set.").format(self.S - 1, "" + "s" * (self.S > 2),
self.mus), 3)
self.samplingEngine.iterSample(self.mus)
self._S = len(self.mus)
self._approxParameters["S"] = self.S
self.muTest = emptyParameterList()
self.muTest.reset((len(muTestBase) + 1, self.mus.shape[1]))
self.muTest.data[: -1] = muTestBase.data
self.muTest.data[-1] = muLast.data
@abstractmethod
def setupApproxLocal(self) -> int:
if self.checkComputedApprox(): return -1
RROMPyAssert(self._mode, message = "Cannot setup approximant.")
vbMng(self, "INIT", "Setting up local approximant.", 5)
pass
vbMng(self, "DEL", "Done setting up local approximant.", 5)
return 0
def addSamplePoints(self, mus:paramList):
"""Add sample points to reduced model."""
raise RROMPyException("Cannot add samples to greedy reduced model.")
def setupApprox(self, plotEst : str = "NONE") -> int:
"""Compute greedy snapshots of solution map."""
if self.checkComputedApprox(): return -1
RROMPyAssert(self._mode, message = "Cannot start greedy algorithm.")
vbMng(self, "INIT", "Setting up {}.". format(self.name()), 5)
vbMng(self, "INIT", "Starting computation of snapshots.", 5)
self._collinearityFlag = 0
self._preliminaryTraining()
muidx, self.firstGreedyIter = [len(self.muTest) - 1], True
errorEstTest, maxErrorEst = [np.inf], np.inf
max2ErrorEst, trainedModelOld = np.inf, None
while self.firstGreedyIter or (len(self.muTest) > 0
and (maxErrorEst is None or max2ErrorEst > self.greedyTol)
and self.samplingEngine.nsamples < self.maxIter):
muTestOld, errorEstTestOld = self.muTest, errorEstTest
muidxOld, maxErrorEstOld = muidx, maxErrorEst
errorEstTest, muidx, maxErrorEst, mu = self.greedyNextSample(muidx,
plotEst)
if maxErrorEst is not None and (np.any(np.isnan(maxErrorEst))
or np.any(np.isinf(maxErrorEst))):
if self._collinearityFlag == 0 and not self.firstGreedyIter:
RROMPyWarning(("Instability in a posteriori "
"estimator. Starting preemptive greedy "
"loop termination."))
self.muTest, errorEstTest = muTestOld, errorEstTestOld
if self.firstGreedyIter and muidx[0] < 0:
self.trainedModel = None
if self._collinearityFlag:
raise RROMPyException(("Starting sample points too "
"collinear. Aborting greedy "
"iterations."))
raise RROMPyException(("Instability in approximant "
"computation. Aborting greedy "
"iterations."))
self._S = trainedModelOld.data.approxParameters["S"]
self._approxParameters["S"] = self.S
while self.samplingEngine.nsamples > self.S:
self.samplingEngine.popSample()
while len(self.mus) > self.S: self.mus.pop(-1)
muidx, maxErrorEst = muidxOld, maxErrorEstOld
break
if maxErrorEst is not None:
max2ErrorEst = np.max(maxErrorEst)
vbMng(self, "MAIN", ("Uniform testing error estimate "
"{:.4e}.").format(max2ErrorEst), 5)
if self.firstGreedyIter:
trainedModelOld = copy(self.trainedModel)
else:
trainedModelOld.data = copy(self.trainedModel.data)
self.firstGreedyIter = False
vbMng(self, "DEL", ("Done computing snapshots (final snapshot count: "
"{}).").format(self.samplingEngine.nsamples), 5)
if (maxErrorEst is None or np.any(np.isnan(maxErrorEst))
or np.any(np.isinf(maxErrorEst))):
while self.samplingEngine.nsamples > self.S:
self.samplingEngine.popSample()
while len(self.mus) > self.S: self.mus.pop(-1)
else:
self._S = self.samplingEngine.nsamples
while len(self.mus) < self.S:
self.mus.append(self.samplingEngine.mus[len(self.mus)])
self.trainedModel = None
self.setupApproxLocal()
if plotEst == "LAST":
self.plotEstimator(errorEstTest, muidx, maxErrorEst)
vbMng(self, "DEL", "Done setting up approximant.", 5)
return 0
def assembleReducedResidualGramian(self, pMat:sampList):
"""
Build residual gramian of reduced linear system through projections.
"""
if (not hasattr(self.trainedModel.data, "gramian")
or self.trainedModel.data.gramian is None):
gramian = self.HFEngine.innerProduct(pMat, pMat, dual = True)
else:
Sold = self.trainedModel.data.gramian.shape[0]
S = len(self.mus)
if Sold > S:
gramian = self.trainedModel.data.gramian[: S, : S]
else:
idxOld = list(range(Sold))
idxNew = list(range(Sold, S))
gramian = np.empty((S, S), dtype = np.complex)
gramian[: Sold, : Sold] = self.trainedModel.data.gramian
gramian[: Sold, Sold :] = self.HFEngine.innerProduct(
pMat(idxNew), pMat(idxOld), dual = True)
gramian[Sold :, : Sold] = gramian[: Sold, Sold :].T.conj()
gramian[Sold :, Sold :] = self.HFEngine.innerProduct(
pMat(idxNew), pMat(idxNew), dual = True)
self.trainedModel.data.gramian = gramian
def assembleReducedResidualBlocksbb(self, bs:List[Np1D]):
"""
Build blocks (of type bb) of reduced linear system through projections.
"""
nbs = len(bs)
if (not hasattr(self.trainedModel.data, "resbb")
or self.trainedModel.data.resbb is None):
resbb = np.empty((nbs, nbs), dtype = np.complex)
for i in range(nbs):
Mbi = bs[i]
resbb[i, i] = self.HFEngine.innerProduct(Mbi, Mbi, dual = True)
for j in range(i):
Mbj = bs[j]
resbb[i, j] = self.HFEngine.innerProduct(Mbj, Mbi,
dual = True)
for i in range(nbs):
for j in range(i + 1, nbs):
resbb[i, j] = resbb[j, i].conj()
self.trainedModel.data.resbb = resbb
def assembleReducedResidualBlocksAb(self, As:List[Np2D], bs:List[Np1D],
pMat:sampList):
"""
Build blocks (of type Ab) of reduced linear system through projections.
"""
nAs = len(As)
nbs = len(bs)
S = len(self.mus)
if (not hasattr(self.trainedModel.data, "resAb")
or self.trainedModel.data.resAb is None):
if isinstance(pMat, (parameterList, sampleList)): pMat = pMat.data
resAb = np.empty((nbs, S, nAs), dtype = np.complex)
for j in range(nAs):
MAj = dot(As[j], pMat)
for i in range(nbs):
Mbi = bs[i]
resAb[i, :, j] = self.HFEngine.innerProduct(MAj, Mbi,
dual = True)
else:
Sold = self.trainedModel.data.resAb.shape[1]
if Sold == S: return
if Sold > S:
resAb = self.trainedModel.data.resAb[:, : S, :]
else:
if isinstance(pMat, (parameterList, sampleList)):
pMat = pMat.data
resAb = np.empty((nbs, S, nAs), dtype = np.complex)
resAb[:, : Sold, :] = self.trainedModel.data.resAb
for j in range(nAs):
MAj = dot(As[j], pMat[:, Sold :])
for i in range(nbs):
Mbi = bs[i]
resAb[i, Sold :, j] = self.HFEngine.innerProduct(
MAj, Mbi, dual = True)
self.trainedModel.data.resAb = resAb
def assembleReducedResidualBlocksAA(self, As:List[Np2D], pMat:sampList):
"""
Build blocks (of type AA) of reduced linear system through projections.
"""
nAs = len(As)
S = len(self.mus)
if (not hasattr(self.trainedModel.data, "resAA")
or self.trainedModel.data.resAA is None):
if isinstance(pMat, (parameterList, sampleList)): pMat = pMat.data
resAA = np.empty((S, nAs, S, nAs), dtype = np.complex)
for i in range(nAs):
MAi = dot(As[i], pMat)
resAA[:, i, :, i] = self.HFEngine.innerProduct(MAi, MAi,
dual = True)
for j in range(i):
MAj = dot(As[j], pMat)
resAA[:, i, :, j] = self.HFEngine.innerProduct(MAj, MAi,
dual = True)
for i in range(nAs):
for j in range(i + 1, nAs):
resAA[:, i, :, j] = resAA[:, j, :, i].T.conj()
else:
Sold = self.trainedModel.data.resAA.shape[0]
if Sold == S: return
if Sold > S:
resAA = self.trainedModel.data.resAA[: S, :, : S, :]
else:
if isinstance(pMat, (parameterList, sampleList)):
pMat = pMat.data
resAA = np.empty((S, nAs, S, nAs), dtype = np.complex)
resAA[: Sold, :, : Sold, :] = self.trainedModel.data.resAA
for i in range(nAs):
MAi = dot(As[i], pMat)
resAA[: Sold, i, Sold :, i] = self.HFEngine.innerProduct(
MAi[:, Sold :], MAi[:, : Sold], dual = True)
resAA[Sold :, i, : Sold, i] = resAA[: Sold, i,
Sold :, i].T.conj()
resAA[Sold :, i, Sold :, i] = self.HFEngine.innerProduct(
MAi[:, Sold :], MAi[:, Sold :], dual = True)
for j in range(i):
MAj = dot(As[j], pMat)
resAA[: Sold, i, Sold :, j] = (
self.HFEngine.innerProduct(MAj[:, Sold :],
MAi[:, : Sold],
dual = True))
resAA[Sold :, i, : Sold, j] = (
self.HFEngine.innerProduct(MAj[:, : Sold],
MAi[:, Sold :],
dual = True))
resAA[Sold :, i, Sold :, j] = (
self.HFEngine.innerProduct(MAj[:, Sold :],
MAi[:, Sold :],
dual = True))
for i in range(nAs):
for j in range(i + 1, nAs):
resAA[: Sold, i, Sold :, j] = (
resAA[Sold :, j, : Sold, i].T.conj())
resAA[Sold :, i, : Sold, j] = (
resAA[: Sold, j, Sold :, i].T.conj())
resAA[Sold :, i, Sold :, j] = (
resAA[Sold :, j, Sold :, i].T.conj())
self.trainedModel.data.resAA = resAA
def assembleReducedResidualBlocks(self, full : bool = False):
"""Build affine blocks of affine decomposition of residual."""
if full:
checkIfAffine(self.HFEngine, "assemble reduced residual blocks",
False, self._affine_lvl)
else:
checkIfAffine(self.HFEngine, "assemble reduced RHS blocks", True,
self._affine_lvl)
self.HFEngine.buildb()
self.assembleReducedResidualBlocksbb(self.HFEngine.bs)
if full:
pMat = self.samplingEngine.projectionMatrix
self.HFEngine.buildA()
self.assembleReducedResidualBlocksAb(self.HFEngine.As,
self.HFEngine.bs, pMat)
self.assembleReducedResidualBlocksAA(self.HFEngine.As, pMat)
diff --git a/rrompy/reduction_methods/standard/greedy/rational_interpolant_greedy.py b/rrompy/reduction_methods/standard/greedy/rational_interpolant_greedy.py
index ce73886..fdf9cb6 100755
--- a/rrompy/reduction_methods/standard/greedy/rational_interpolant_greedy.py
+++ b/rrompy/reduction_methods/standard/greedy/rational_interpolant_greedy.py
@@ -1,480 +1,440 @@
# Copyright (C) 2018-2020 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.hfengines.base.linear_affine_engine import checkIfAffine
from .generic_greedy_approximant import GenericGreedyApproximant
-from rrompy.utilities.poly_fitting.polynomial import (polyfitname, polyvander,
- PolynomialInterpolator as PI)
+from rrompy.utilities.poly_fitting.polynomial import polyvander
from rrompy.utilities.numerical import dot
+from rrompy.utilities.numerical.degree import degreeToDOFs
from rrompy.utilities.expression import expressionEvaluator
from rrompy.reduction_methods.standard import RationalInterpolant
from rrompy.utilities.base.types import Np1D, Tuple, paramVal, List
from rrompy.utilities.base.verbosity_depth import verbosityManager as vbMng
from rrompy.utilities.poly_fitting import customFit
from rrompy.utilities.exception_manager import (RROMPyWarning, RROMPyException,
RROMPyAssert, RROMPy_FRAGILE)
from rrompy.sampling import sampleList, emptySampleList
__all__ = ['RationalInterpolantGreedy']
class RationalInterpolantGreedy(GenericGreedyApproximant, RationalInterpolant):
"""
ROM greedy rational interpolant computation for parametric problems.
Args:
HFEngine: HF problem solver.
mu0(optional): Default parameter. Defaults to 0.
approxParameters(optional): Dictionary containing values for main
parameters of approximant. Recognized keys are:
- 'POD': kind of snapshots orthogonalization; allowed values
include 0, 1/2, and 1; defaults to 1, i.e. POD;
- 'scaleFactorDer': scaling factors for derivative computation;
defaults to 'AUTO';
- 'S': number of starting training points;
+ - 'polydegreetype': type of polynomial degree; allowed values are
+ '*' or '*_*', with * either 'TENSOR' or 'TOTAL'; defaults to
+ 'TOTAL';
- 'N': degree of rational interpolant denominator; defaults to
'AUTO', i.e. maximum allowed;
+ - 'MAuxiliary': degree of rational interpolant numerator with
+ respect to non-pivot parameters; defaults to 0;
+ - 'NAuxiliary': degree of rational interpolant denominator with
+ respect to non-pivot parameters; defaults to 0;
- 'sampler': sample point generator;
- 'greedyTol': uniform error tolerance for greedy algorithm;
defaults to 1e-2;
- 'collinearityTol': collinearity tolerance for greedy algorithm;
defaults to 0.;
- 'maxIter': maximum number of greedy steps; defaults to 1e2;
- 'nTestPoints': number of test points; defaults to 5e2;
- 'samplerTrainSet': training sample points generator; defaults to
sampler;
- 'polybasis': type of basis for interpolation; defaults to
'MONOMIAL';
- 'errorEstimatorKind': kind of error estimator; available values
include 'AFFINE', 'DISCREPANCY', 'LOOK_AHEAD',
'LOOK_AHEAD_RES', 'LOOK_AHEAD_OUTPUT', and 'NONE'; defaults to
'NONE';
- 'functionalSolve': strategy for minimization of denominator
functional; allowed values include 'NORM', 'DOMINANT',
'BARYCENTRIC[_NORM]', and 'BARYCENTRIC_AVERAGE' (check pdf in
main folder for explanation); defaults to 'NORM';
- 'interpTol': tolerance for interpolation; defaults to None;
- 'forceQReal': force denominator to have real coefficients;
defaults to False;
- 'polyTruncateTol': tolerance for truncation of rational terms;
defaults to 0;
- 'QTol': tolerance for robust rational denominator management;
defaults to 0.
Defaults to empty dict.
verbosity(optional): Verbosity level. Defaults to 10.
Attributes:
HFEngine: HF problem solver.
mu0: Default parameter.
mus: Array of snapshot parameters.
approxParameters: Dictionary containing values for main parameters of
approximant. Recognized keys are in parameterList.
parameterListSoft: Recognized keys of soft approximant parameters:
- 'POD': kind of snapshots orthogonalization;
- 'scaleFactorDer': scaling factors for derivative computation;
+ - 'polydegreetype': type of polynomial degree;
- 'N': degree of rational interpolant denominator;
+ - 'MAuxiliary': degree of rational interpolant numerator with
+ respect to non-pivot parameters; defaults to 0;
+ - 'NAuxiliary': degree of rational interpolant denominator with
+ respect to non-pivot parameters; defaults to 0;
- 'greedyTol': uniform error tolerance for greedy algorithm;
- 'collinearityTol': collinearity tolerance for greedy algorithm;
- 'maxIter': maximum number of greedy steps;
- 'nTestPoints': number of test points;
- 'samplerTrainSet': training sample points generator;
- 'errorEstimatorKind': kind of error estimator;
- 'functionalSolve': strategy for minimization of denominator
functional;
- 'interpTol': tolerance for interpolation;
- 'forceQReal': force denominator to have real coefficients;
- 'polyTruncateTol': tolerance for truncation of rational terms;
- 'QTol': tolerance for robust rational denominator management.
parameterListCritical: Recognized keys of critical approximant
parameters:
- 'S': total number of samples current approximant relies upon;
- 'sampler': sample point generator.
verbosity: Verbosity level.
POD: Kind of snapshots orthogonalization.
scaleFactorDer: Scaling factors for derivative computation.
S: number of test points.
+ polydegreetype: Type of polynomial degree.
N: Denominator degree of approximant.
+ MAuxiliary: Degree of rational interpolant numerator with respect to
+ non-pivot parameters.
+ NAuxiliary: Degree of rational interpolant denominator with respect to
+ non-pivot parameters.
sampler: Sample point generator.
greedyTol: uniform error tolerance for greedy algorithm.
collinearityTol: Collinearity tolerance for greedy algorithm.
maxIter: maximum number of greedy steps.
nTestPoints: number of starting training points.
samplerTrainSet: training sample points generator.
errorEstimatorKind: kind of error estimator.
functionalSolve: Strategy for minimization of denominator functional.
interpTol: tolerance for interpolation.
forceQReal: Force denominator to have real coefficients.
polyTruncateTol: Tolerance for truncation of rational terms.
QTol: tolerance for robust rational denominator management.
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.
"""
_allowedEstimatorKinds = ["AFFINE", "DISCREPANCY", "LOOK_AHEAD",
"LOOK_AHEAD_RES", "LOOK_AHEAD_OUTPUT", "NONE"]
def __init__(self, *args, **kwargs):
self._preInit()
self._addParametersToList(["errorEstimatorKind"], ["DISCREPANCY"],
- toBeExcluded = ["rationalMode", "M",
- "polydegreetype",
- "indexPrimary", "MAuxiliary",
- "NAuxiliary"])
+ toBeExcluded = ["rationalMode", "M"])
super().__init__(*args, **kwargs)
self.M = "AUTO"
self._postInit()
@property
def rationalMode(self): return "MINIMAL_STATE"
- def _setMAuto(self):
- self.M = self.S - 1
-
- @property
- def N(self):
- """Value of N."""
- return self._N
- _N_fix = None
- @N.setter
- def N(self, N):
- self._N_isauto, self._N_shift = True, 0
- if isinstance(N, str):
- N = N.strip().replace(" ","")
- if "-" in N: self._N_shift = int(N.split("-")[-1])
- self._N_fix, N = -1, 0
- elif self._N_fix is None:
- self._N_fix = N
- if N < 0: raise RROMPyException("N must be non-negative.")
- self._N = N
- self._approxParameters["N"] = self.N
-
- def _setNAuto(self):
- N = max(0, self.S - self._N_shift - 1)
- if self._N_fix >= 0: N = min(N, self._N_fix)
- self.N = N
-
- @property
- def indexPrimary(self): return 0
-
- @property
- def MAuxiliary(self): return 0
-
- @property
- def NAuxiliary(self): return 0
-
- @property
- def polydegreetype(self): return "TENSOR"
-
@property
def errorEstimatorKind(self):
"""Value of errorEstimatorKind."""
return self._errorEstimatorKind
@errorEstimatorKind.setter
def errorEstimatorKind(self, errorEstimatorKind):
errorEstimatorKind = errorEstimatorKind.upper()
if errorEstimatorKind not in self._allowedEstimatorKinds:
RROMPyWarning(("Error estimator kind not recognized. Overriding "
"to 'NONE'."))
errorEstimatorKind = "NONE"
self._errorEstimatorKind = errorEstimatorKind
self._approxParameters["errorEstimatorKind"] = self.errorEstimatorKind
- def _polyvanderAuxiliary(self, mus, deg, *args):
- return polyvander(mus, deg, *args)
+ def _polyvanderAuxiliary(self, mus, deg, *args, **kwargs):
+ return polyvander(mus, deg, *args, **kwargs)
def getErrorEstimatorDiscrepancy(self, mus:Np1D) -> Np1D:
"""Discrepancy-based residual estimator."""
checkIfAffine(self.HFEngine, "apply discrepancy-based error estimator",
False, self._affine_lvl)
mus = self.checkParameterList(mus)
muCTest = self.trainedModel.centerNormalize(mus)
tMverb, self.trainedModel.verbosity = self.trainedModel.verbosity, 0
QTest = self.trainedModel.getQVal(mus)
QTzero = np.where(QTest == 0.)[0]
if len(QTzero) > 0:
RROMPyWarning(("Adjusting estimator to avoid division by "
"numerically zero denominator."))
- QTest[QTzero] = np.finfo(np.complex).eps / (1. + self.N)
+ QTest[QTzero] = (np.finfo(np.complex).eps
+ / len(self.trainedModel.data.Q))
self.HFEngine.buildA()
self.HFEngine.buildb()
nAs, nbs = self.HFEngine.nAs, self.HFEngine.nbs
muTrainEff = self.mapParameterList(self.mus)
muTestEff = self.mapParameterList(mus)
PTrain = self.trainedModel.getPVal(self.mus).data.T
QTrain = self.trainedModel.getQVal(self.mus)
QTzero = np.where(QTrain == 0.)[0]
if len(QTzero) > 0:
RROMPyWarning(("Adjusting estimator to avoid division by "
"numerically zero denominator."))
- QTrain[QTzero] = np.finfo(np.complex).eps / (1. + self.N)
+ QTrain[QTzero] = (np.finfo(np.complex).eps
+ / len(self.trainedModel.data.Q))
PTest = self.trainedModel.getPVal(mus).data
self.trainedModel.verbosity = tMverb
radiusAbTrain = np.empty((self.S, nAs * self.S + nbs),
dtype = np.complex)
radiusA = np.empty((self.S, nAs, len(mus)), dtype = np.complex)
radiusb = np.empty((nbs, len(mus)), dtype = np.complex)
for j, thA in enumerate(self.HFEngine.thAs):
idxs = j * self.S + np.arange(self.S)
radiusAbTrain[:, idxs] = expressionEvaluator(thA[0], muTrainEff,
(self.S, 1)) * PTrain
radiusA[:, j] = PTest * expressionEvaluator(thA[0], muTestEff,
(len(mus),))
for j, thb in enumerate(self.HFEngine.thbs):
idx = nAs * self.S + j
radiusAbTrain[:, idx] = QTrain * expressionEvaluator(thb[0],
muTrainEff, (self.S,))
radiusb[j] = QTest * expressionEvaluator(thb[0], muTestEff,
(len(mus),))
QRHSNorm2 = self._affineResidualMatricesContraction(radiusb)
- vanTrain = self._polyvanderAuxiliary(self._musUniqueCN, self.S - 1,
+ deg = [self.M] + [self.MAuxiliary] * (self.npar - 1)
+ vanTrain = self._polyvanderAuxiliary(self._musUniqueCN, deg,
self.polybasis, self._derIdxs,
self._reorder)
interpPQ = customFit(vanTrain, radiusAbTrain, rcond = self.interpTol)
- vanTest = self._polyvanderAuxiliary(muCTest, self.S - 1,
- self.polybasis)
+ vanTest = self._polyvanderAuxiliary(muCTest, deg, self.polybasis)
DradiusAb = vanTest.dot(interpPQ)
radiusA = (radiusA
- DradiusAb[:, : - nbs].reshape(len(mus), -1, self.S).T)
radiusb = radiusb - DradiusAb[:, - nbs :].T
ff, Lf, LL = self._affineResidualMatricesContraction(radiusb, radiusA)
err = np.abs((LL - 2. * np.real(Lf) + ff) / QRHSNorm2) ** .5
return err
def getErrorEstimatorLookAhead(self, mus:Np1D,
what : str = "") -> Tuple[Np1D, List[int]]:
"""Residual estimator based on look-ahead idea."""
- errTest, QTest, idxMaxEst = self._EIMStep(mus)
+ err, idxMaxEst = self._EIMStep(mus)
mu_muTestS = mus[idxMaxEst]
app_muTestSample = self.trainedModel.getApproxReduced(mu_muTestS)
if self._mode == RROMPy_FRAGILE:
if what == "RES" and not self.HFEngine.isCEye:
raise RROMPyException(("Cannot compute LOOK_AHEAD_RES "
"estimator in fragile mode for "
"non-scalar C."))
app_muTestSample = dot(self.trainedModel.data.projMat[:,
: app_muTestSample.shape[0]],
app_muTestSample)
else:
app_muTestSample = dot(self.samplingEngine.projectionMatrix,
app_muTestSample)
app_muTestSample = sampleList(app_muTestSample)
if what == "RES":
errmu = self.HFEngine.residual(mu_muTestS, app_muTestSample,
post_c = False)
solmu = self.HFEngine.residual(mu_muTestS, None, post_c = False)
normSol = self.HFEngine.norm(solmu, dual = True)
normErr = self.HFEngine.norm(errmu, dual = True)
else:
for j, mu in enumerate(mu_muTestS):
uEx = self.samplingEngine.nextSample(mu)
if what == "OUTPUT":
uEx = self.HFEngine.applyC(uEx, mu)
app_muTS = self.HFEngine.applyC(app_muTestSample[j], mu)
if j == 0:
app_muTestS = emptySampleList()
app_muTestS.reset((len(app_muTS), len(mu_muTestS)),
dtype = np.complex)
app_muTestS[j] = app_muTS
if j == 0:
solmu = emptySampleList()
solmu.reset((len(uEx), len(mu_muTestS)),
dtype = np.complex)
solmu[j] = uEx
if what == "OUTPUT": app_muTestSample = app_muTestS
errmu = solmu - app_muTestSample
normSol = self.HFEngine.norm(solmu, is_state = what != "OUTPUT")
normErr = self.HFEngine.norm(errmu, is_state = what != "OUTPUT")
- errsamples = normErr / normSol
- musT = copy(self.mus)
- musT.append(mu_muTestS)
- musT = self.trainedModel.centerNormalize(musT)
- musC = self.trainedModel.centerNormalize(mus)
- errT = np.zeros((len(musT), len(mu_muTestS)), dtype = np.complex)
- errT[np.arange(len(self.mus), len(musT)),
- np.arange(len(mu_muTestS))] = errsamples * QTest[idxMaxEst]
- vanT = self._polyvanderAuxiliary(musT, self.S, self.polybasis)
- fitOut = customFit(vanT, errT, full = True, rcond = self.interpTol)
- vbMng(self, "MAIN",
- ("Fitting {} samples with degree {} through {}... Conditioning "
- "of LS system: {:.4e}.").format(len(vanT), self.S,
- polyfitname(self.polybasis),
- fitOut[1][2][0] / fitOut[1][2][-1]), 15)
- vanC = self._polyvanderAuxiliary(musC, self.S, self.polybasis)
- err = np.sum(np.abs(vanC.dot(fitOut[0])), axis = -1) / QTest
- return err, idxMaxEst
+ errTestMean = np.mean(err[idxMaxEst])
+ errMean = np.mean(normErr / normSol)
+ return err * errMean / errTestMean, idxMaxEst
def getErrorEstimatorNone(self, mus:Np1D) -> Np1D:
"""EIM-based residual estimator."""
- err = np.max(self._EIMStep(mus, True), axis = 1)
+ err = self._EIMStep(mus, True)
err *= self.greedyTol / np.mean(err)
return err
def _EIMStep(self, mus:Np1D,
- only_one : bool = False) -> Tuple[Np1D, Np1D, List[int]]:
+ only_one : bool = False) -> Tuple[Np1D, List[int]]:
"""EIM step to find next magic point."""
mus = self.checkParameterList(mus)
tMverb, self.trainedModel.verbosity = self.trainedModel.verbosity, 0
QTest = self.trainedModel.getQVal(mus)
QTzero = np.where(QTest == 0.)[0]
if len(QTzero) > 0:
RROMPyWarning(("Adjusting estimator to avoid division by "
"numerically zero denominator."))
- QTest[QTzero] = np.finfo(np.complex).eps / (1. + self.N)
- QTest = np.abs(QTest)
+ QTest[QTzero] = (np.finfo(np.complex).eps
+ / len(self.trainedModel.data.Q))
muCTest = self.trainedModel.centerNormalize(mus)
muCTrain = self.trainedModel.centerNormalize(self.mus)
self.trainedModel.verbosity = tMverb
- vanTest = self._polyvanderAuxiliary(muCTest, self.S - 1,
- self.polybasis)
- vanTestNext = self._polyvanderAuxiliary(muCTest, self.S,
- self.polybasis)[:,
- vanTest.shape[1] :]
+ deg = [self.M] + [self.MAuxiliary] * (self.npar - 1)
+ S0 = degreeToDOFs(deg, self.polydegreetypeM)
+ if S0 != len(muCTrain):
+ raise RROMPyException(("Cannot evaluate error estimator for LS "
+ "interpolation."))
+ deg[0] = self.M + 1
+ vanTest = self._polyvanderAuxiliary(muCTest, deg, self.polybasis,
+ degreetype = self.polydegreetypeM)
+ vanTest, vanTestNext = vanTest[:, : S0], vanTest[:, S0 :]
idxsTest = np.arange(vanTestNext.shape[1])
basis = np.zeros((len(idxsTest), 0), dtype = float)
- idxMaxEst = []
+ idxMaxEst = np.empty(len(idxsTest), dtype = int)
while len(idxsTest) > 0:
- vanTrial = self._polyvanderAuxiliary(muCTrain, self.S - 1,
- self.polybasis)
- vanTrialNext = self._polyvanderAuxiliary(muCTrain, self.S,
- self.polybasis)[:,
- vanTrial.shape[1] :]
+ vanTrial = self._polyvanderAuxiliary(muCTrain, deg, self.polybasis,
+ degreetype = self.polydegreetypeM)
+ vanTrial, vanTrialNext = vanTrial[:, : S0], vanTrial[:, S0 :]
vanTrial = np.hstack((vanTrial, vanTrialNext.dot(basis).reshape(
len(vanTrialNext), basis.shape[1])))
valuesTrial = vanTrialNext[:, idxsTest]
vanTestEff = np.hstack((vanTest, vanTestNext.dot(basis).reshape(
len(vanTestNext), basis.shape[1])))
vanTestNextEff = vanTestNext[:, idxsTest]
coeffTest = np.linalg.solve(vanTrial, valuesTrial)
- errTest = (np.abs(vanTestNextEff - vanTestEff.dot(coeffTest))
- / np.expand_dims(QTest, 1))
- if only_one: return errTest
+ errTest = np.abs((vanTestNextEff - vanTestEff.dot(coeffTest))
+ / np.expand_dims(QTest, 1))
+ if only_one:
+ return errTest
+ elif basis.shape[1] == 0:
+ errTest0 = errTest[:, 0]
idxMaxErr = np.unravel_index(np.argmax(errTest), errTest.shape)
- idxMaxEst += [idxMaxErr[0]]
+ idxMaxEst[idxsTest[idxMaxErr[1]]] = idxMaxErr[0]
muCTrain.append(muCTest[idxMaxErr[0]])
basis = np.pad(basis, [(0, 0), (0, 1)], "constant")
basis[idxsTest[idxMaxErr[1]], -1] = 1.
idxsTest = np.delete(idxsTest, idxMaxErr[1])
- return errTest, QTest, idxMaxEst
+ return errTest0, idxMaxEst
def errorEstimator(self, mus:Np1D, return_max : bool = False) -> Np1D:
"""Standard residual-based error estimator."""
setupOK = self.setupApproxLocal()
if setupOK > 0:
err = np.empty(len(mus))
err[:] = np.nan
if not return_max: return err
return err, [- setupOK], np.nan
mus = self.checkParameterList(mus)
vbMng(self.trainedModel, "INIT",
"Evaluating error estimator at mu = {}.".format(mus), 10)
if self.errorEstimatorKind == "AFFINE":
err = self.getErrorEstimatorAffine(mus)
else:
self._setupInterpolationIndices()
if self.errorEstimatorKind == "DISCREPANCY":
err = self.getErrorEstimatorDiscrepancy(mus)
elif self.errorEstimatorKind[: 10] == "LOOK_AHEAD":
err, idxMaxEst = self.getErrorEstimatorLookAhead(mus,
self.errorEstimatorKind[11 :])
else: #if self.errorEstimatorKind == "NONE":
err = self.getErrorEstimatorNone(mus)
vbMng(self.trainedModel, "DEL", "Done evaluating error estimator.", 10)
if not return_max: return err
if self.errorEstimatorKind[: 10] != "LOOK_AHEAD":
idxMaxEst = [np.argmax(err)]
return err, idxMaxEst, err[idxMaxEst]
_warnPlottingNormalization = 1
def plotEstimator(self, *args, **kwargs):
super().plotEstimator(*args, **kwargs)
if (self.errorEstimatorKind == "NONE"
and self._warnPlottingNormalization):
RROMPyWarning(("Error estimator arbitrarily normalized before "
"plotting."))
self._warnPlottingNormalization = 0
def greedyNextSample(self, *args,
**kwargs) -> Tuple[Np1D, int, float, paramVal]:
"""Compute next greedy snapshot of solution map."""
RROMPyAssert(self._mode, message = "Cannot add greedy sample.")
err, muidx, maxErr, muNext = super().greedyNextSample(*args, **kwargs)
if (self.errorEstimatorKind == "NONE" and not np.isnan(maxErr)
and not np.isinf(maxErr)):
maxErr = None
return err, muidx, maxErr, muNext
def _preliminaryTraining(self):
"""Initialize starting snapshots of solution map."""
RROMPyAssert(self._mode, message = "Cannot start greedy algorithm.")
if self.samplingEngine.nsamples > 0: return
- if self.npar > 1:
- raise RROMPyException(("Cannot apply minimal rational interpolant "
- "in multivariate case."))
super()._preliminaryTraining()
def setupApproxLocal(self) -> int:
"""Compute rational interpolant."""
if self.checkComputedApprox(): return -1
RROMPyAssert(self._mode, message = "Cannot setup approximant.")
self.verbosity -= 10
vbMng(self, "INIT", "Setting up local approximant.", 5)
pMat = self.samplingEngine.projectionMatrix
firstRun = self.trainedModel is None
if not firstRun: pMat = pMat[:, len(self.trainedModel.data.mus) :]
self._setupTrainedModel(pMat, not firstRun)
+ Q = self._setupDenominator()
unstable = 0
- if self.S > 1:
- Q = self._setupDenominator()
- else:
- Q = PI()
- Q.coeffs = np.ones(1, dtype = np.complex)
- Q.npar, Q.polybasis = 1, self.polybasis
- if not unstable:
- self._setupRational(Q)
- if self.M < self.S - 1:
- RROMPyWarning(("Instability in numerator computation. "
- "Aborting."))
- unstable = 1
+ if not firstRun: _M = self.M
+ self._setupRational(Q)
+ if not firstRun and self.M < _M:
+ RROMPyWarning(("Instability in numerator computation. "
+ "Aborting."))
+ unstable = 1
if not unstable:
self.trainedModel.data.approxParameters = copy(
self.approxParameters)
vbMng(self, "DEL", "Done setting up local approximant.", 5)
self.verbosity += 10
return unstable
def setupApprox(self, plotEst : str = "NONE") -> int:
val = super().setupApprox(plotEst)
if val == 0:
self._setupRational(self.trainedModel.data.Q,
self.trainedModel.data.P)
self.trainedModel.data.approxParameters = copy(
self.approxParameters)
return val
diff --git a/rrompy/reduction_methods/standard/rational_interpolant.py b/rrompy/reduction_methods/standard/rational_interpolant.py
index b492c9f..b1ef6b9 100755
--- a/rrompy/reduction_methods/standard/rational_interpolant.py
+++ b/rrompy/reduction_methods/standard/rational_interpolant.py
@@ -1,1026 +1,1016 @@
# Copyright (C) 2018-2020 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 scipy.sparse import coo_matrix
from scipy.linalg import eig, eigvals, khatri_rao
from .generic_standard_approximant import GenericStandardApproximant
from rrompy.utilities.poly_fitting.polynomial import (
polybases, polyfitname,
polyvander as pvP, polyTimes,
PolynomialInterpolator as PI,
PolynomialInterpolatorNodal as PIN)
from rrompy.utilities.poly_fitting.heaviside import rational2heaviside
from rrompy.utilities.base.types import (Np1D, Np2D, Tuple, List, paramList,
interpEng)
from rrompy.utilities.base import verbosityManager as vbMng
from rrompy.utilities.numerical import pseudoInverse, dot, baseDistanceMatrix
from rrompy.utilities.numerical.factorials import multifactorial
from rrompy.utilities.numerical.degree import (mapTotalToTensor, degreeToDOFs,
reduceDegreeMN, reduceDegreeN)
from rrompy.utilities.numerical.hash_derivative import (nextDerivativeIndices,
hashDerivativeToIdx as hashD,
hashIdxToDerivative as hashI)
from rrompy.utilities.exception_manager import (RROMPyException, RROMPyAssert,
RROMPyWarning)
__all__ = ['RationalInterpolant']
def polyTimesTable(P:interpEng, mus:Np1D, reorder:List[int],
derIdxs:List[List[List[int]]], scl : Np1D = None) -> Np2D:
"""Table of polynomial products."""
if not isinstance(P, PI):
raise RROMPyException(("Polynomial to evaluate must be a polynomial "
"interpolator."))
Pvals = [[0.] * len(derIdx) for derIdx in derIdxs]
for j, derIdx in enumerate(derIdxs):
nder = len(derIdx)
for der in range(nder):
derI = hashI(der, P.npar)
Pvals[j][der] = P([mus[j]], derI, scl) / multifactorial(derI)
return blockDiagDer(Pvals, reorder, derIdxs, format = "dense")
def vanderInvTable(vanInv:Np2D, reorder:List[int],
derIdxs:List[List[List[int]]]) -> Np2D:
"""Table of Vandermonde pseudo-inverse."""
S = len(reorder)
Ts = [None] * len(vanInv)
for k in range(len(vanInv)):
invLocs = [None] * len(derIdxs)
idxGlob = 0
for j, derIdx in enumerate(derIdxs):
nder = len(derIdx)
idxGlob += nder
idxLoc = np.arange(S)[(reorder >= idxGlob - nder)
* (reorder < idxGlob)]
invLocs[j] = vanInv[k, idxLoc]
Ts[k] = blockDiagDer(invLocs, reorder, derIdxs, [2, 1, 0])
return Ts
def blockDiagDer(vals:List[Np1D], reorder:List[int],
derIdxs:List[List[List[int]]],
permute : List[int] = None, format : str = "sparse") -> Np2D:
"""Table of derivative values for point confluence."""
S = len(reorder)
if format == "sparse":
data, idxI, idxJ = [], [], []
else: #format == "dense":
T = np.zeros((S, S), dtype = np.complex)
if permute is None: permute = [0, 1, 2]
idxGlob = 0
for j, derIdx in enumerate(derIdxs):
nder = len(derIdx)
idxGlob += nder
idxLoc = np.arange(S)[(reorder >= idxGlob - nder)
* (reorder < idxGlob)]
val = vals[j]
for derI, derIdxI in enumerate(derIdx):
for derJ, derIdxJ in enumerate(derIdx):
diffIdx = [x - y for (x, y) in zip(derIdxI, derIdxJ)]
if all([x >= 0 for x in diffIdx]):
diffj = hashD(diffIdx)
i1, i2, i3 = np.array([derI, derJ, diffj])[permute]
if format == "sparse":
data += [val[i3]]
idxI += [idxLoc[i1]]
idxJ += [idxLoc[i2]]
else:
T[idxLoc[i1], idxLoc[i2]] = val[i3]
if format == "sparse":
T = coo_matrix((data, (idxI, idxJ)), shape = (S, S)).tocsr()
return T
-def fullDegreeMaxMask(degs:List[int], index:int) -> List[int]:
- deg, npar = degs[index], len(degs)
- degS = degs[0] if index or npar < 2 else degs[1]
- nLeft = (1 + degS) ** (npar - index - 1)
- nRight = (1 + degS) ** (index)
- baseidx = deg * nLeft + np.arange(nLeft)
- return (np.tile(baseidx.reshape(-1, 1), (1, nRight))
- + np.arange(0, (deg + 1) * nLeft * nRight, (deg + 1) * nLeft)
- ).T.flatten()
-
def barycentricRoots(coeffs:Np1D, supp:Np1D,
return_vectors : bool = False) -> Np1D:
m = len(supp)
if len(coeffs) == m: coeffs = np.append(0., coeffs)
arrow = np.diag(np.append(0., supp) + 0.j)
arrow[0], arrow[1 :, 0] = coeffs, 1.
active = np.diag(np.append(0., np.ones(m)))
if return_vectors:
roots, vectors = eig(arrow, active)
else:
roots = eigvals(arrow, active)
badRoots = np.isinf(roots) + np.isnan(roots) == False
if return_vectors: return roots[badRoots], vectors[:, badRoots]
return roots[badRoots]
class RationalInterpolant(GenericStandardApproximant):
"""
ROM rational interpolant computation for parametric problems.
Args:
HFEngine: HF problem solver.
mu0(optional): Default parameter. Defaults to 0.
approxParameters(optional): Dictionary containing values for main
parameters of approximant. Recognized keys are:
- 'POD': kind of snapshots orthogonalization; allowed values
include 0, 1/2, and 1; defaults to 1, i.e. POD;
- 'scaleFactorDer': scaling factors for derivative computation;
defaults to 'AUTO';
- 'S': total number of samples current approximant relies upon;
- 'sampler': sample point generator;
- 'rationalMode': mode of rational approximation; allowed values
include 'MINIMAL[_STATE]', 'MINIMAL_OUTPUT',
'STANDARD[_STATE]', and 'STANDARD_OUTPUT'; defaults to
'MINIMAL';
- 'polybasis': type of polynomial basis for interpolation; defaults
to 'MONOMIAL';
- - 'polydegreetype': type of polynomial degree; allowed values
- include 'TENSOR' and 'TOTAL'; defaults to 'TOTAL';
+ - 'polydegreetype': type of polynomial degree; allowed values are
+ '*' or '*_*', with * either 'TENSOR' or 'TOTAL'; defaults to
+ 'TOTAL';
- 'M': degree of rational interpolant numerator; defaults to
'AUTO', i.e. maximum allowed;
- 'N': degree of rational interpolant denominator; defaults to
'AUTO', i.e. maximum allowed;
- - 'indexPrimary': index of pivot parameter; defaults to 0;
- 'MAuxiliary': degree of rational interpolant numerator with
respect to non-pivot parameters; defaults to 0;
- 'NAuxiliary': degree of rational interpolant denominator with
respect to non-pivot parameters; defaults to 0;
- 'functionalSolve': strategy for minimization of denominator
functional; allowed values include 'NORM', 'DOMINANT',
'BARYCENTRIC[_NORM]', and 'BARYCENTRIC_AVERAGE' (check pdf in
main folder for explanation); defaults to 'NORM';
- 'interpTol': tolerance for interpolation; defaults to None;
- 'forceQReal': force denominator to have real coefficients;
defaults to False;
- 'polyTruncateTol': tolerance for truncation of rational terms;
defaults to 0;
- 'QTol': tolerance for robust rational denominator management;
defaults to 0.
Defaults to empty dict.
verbosity(optional): Verbosity level. Defaults to 10.
Attributes:
HFEngine: HF problem solver.
mu0: Default parameter.
mus: Array of snapshot parameters.
approxParameters: Dictionary containing values for main parameters of
approximant. Recognized keys are in parameterList.
parameterListSoft: Recognized keys of soft approximant parameters:
- 'POD': kind of snapshots orthogonalization;
- 'scaleFactorDer': scaling factors for derivative computation;
- 'rationalMode': mode of rational approximation;
- 'polybasis': type of polynomial basis for interpolation;
- 'polydegreetype': type of polynomial degree;
- 'M': degree of rational interpolant numerator;
- 'N': degree of rational interpolant denominator;
- - 'indexPrimary': index of pivot parameter; defaults to 0;
- 'MAuxiliary': degree of rational interpolant numerator with
respect to non-pivot parameters; defaults to 0;
- 'NAuxiliary': degree of rational interpolant denominator with
respect to non-pivot parameters; defaults to 0;
- 'functionalSolve': strategy for minimization of denominator
functional;
- 'interpTol': tolerance for interpolation via numpy.polyfit;
- 'forceQReal': force denominator to have real coefficients;
- 'polyTruncateTol': tolerance for truncation of rational terms;
- 'QTol': tolerance for robust rational denominator management.
parameterListCritical: Recognized keys of critical approximant
parameters:
- 'S': total number of samples current approximant relies upon;
- 'sampler': sample point generator.
verbosity: Verbosity level.
POD: Kind of snapshots orthogonalization.
scaleFactorDer: Scaling factors for derivative computation.
S: Number of solution snapshots over which current approximant is
based upon.
sampler: Sample point generator.
rationalMode: Mode of rational approximation.
polybasis: Type of polynomial basis for interpolation.
polydegreetype: Type of polynomial degree.
M: Numerator degree of approximant.
N: Denominator degree of approximant.
- indexPrimary: Index of pivot parameter.
MAuxiliary: Degree of rational interpolant numerator with respect to
non-pivot parameters.
NAuxiliary: Degree of rational interpolant denominator with respect to
non-pivot parameters.
functionalSolve: Strategy for minimization of denominator functional.
interpTol: Tolerance for interpolation via numpy.polyfit.
forceQReal: Force denominator to have real coefficients.
polyTruncateTol: Tolerance for truncation of rational terms.
QTol: Tolerance for robust rational denominator management.
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.
Q: Numpy 1D vector containing complex coefficients of approximant
denominator.
P: Numpy 2D vector whose columns are FE dofs of coefficients of
approximant numerator.
"""
_allowedFunctionalSolveKinds = ["NORM", "DOMINANT", "BARYCENTRIC_NORM",
"BARYCENTRIC_AVERAGE"]
def __init__(self, *args, **kwargs):
self._preInit()
self._addParametersToList(["rationalMode", "polybasis", "M", "N",
- "polydegreetype", "indexPrimary",
- "MAuxiliary", "NAuxiliary",
- "functionalSolve", "interpTol",
- "forceQReal", "polyTruncateTol", "QTol"],
+ "polydegreetype", "MAuxiliary",
+ "NAuxiliary", "functionalSolve",
+ "interpTol", "forceQReal",
+ "polyTruncateTol", "QTol"],
["MINIMAL", "MONOMIAL", "AUTO", "AUTO",
- "TOTAL", 0, 0, 0, "NORM", -1, 0, 0., 0.])
+ "TOTAL", 0, 0, "NORM", -1, 0, 0., 0.])
super().__init__(*args, **kwargs)
self._postInit()
@property
def tModelType(self):
from .trained_model.trained_model_rational import TrainedModelRational
return TrainedModelRational
@property
def rationalMode(self):
"""Value of rationalMode."""
return self._rationalMode
@rationalMode.setter
def rationalMode(self, rationalMode):
try:
rationalMode = rationalMode.upper().strip().replace(" ","")
if rationalMode in ["MINIMAL", "STANDARD"]:
rationalMode += "_STATE"
if rationalMode not in ["MINIMAL_STATE", "MINIMAL_OUTPUT",
"STANDARD_STATE", "STANDARD_OUTPUT"]:
raise RROMPyException("Prescribed mode not recognized.")
except Exception as E:
RROMPyWarning(str(E) + " Overriding to 'MINIMAL'.")
rationalMode = "MINIMAL_STATE"
self._rationalMode = rationalMode
self._approxParameters["rationalMode"] = self.rationalMode
@property
def polybasis(self):
"""Value of polybasis."""
return self._polybasis
@polybasis.setter
def polybasis(self, polybasis):
try:
polybasis = polybasis.upper().strip().replace(" ","")
if polybasis not in polybases:
raise RROMPyException("Prescribed polybasis not recognized.")
except Exception as E:
RROMPyWarning(str(E) + " Overriding to 'MONOMIAL'.")
polybasis = "MONOMIAL"
self._polybasis = polybasis
self._approxParameters["polybasis"] = self.polybasis
@property
def functionalSolve(self):
"""Value of functionalSolve."""
return self._functionalSolve
@functionalSolve.setter
def functionalSolve(self, functionalSolve):
try:
functionalSolve = functionalSolve.upper().strip().replace(" ","")
if functionalSolve == "BARYCENTRIC": functionalSolve += "_NORM"
if functionalSolve not in self._allowedFunctionalSolveKinds:
raise RROMPyException(("Prescribed functionalSolve not "
"recognized."))
except Exception as E:
RROMPyWarning(str(E) + " Overriding to 'NORM'.")
functionalSolve = "NORM"
self._functionalSolve = functionalSolve
self._approxParameters["functionalSolve"] = self.functionalSolve
@property
def interpTol(self):
"""Value of interpTol."""
return self._interpTol
@interpTol.setter
def interpTol(self, interpTol):
self._interpTol = interpTol
self._approxParameters["interpTol"] = self.interpTol
@property
def M(self):
"""Value of M."""
return self._M
@M.setter
def M(self, M):
if isinstance(M, str):
M = M.strip().replace(" ","")
if "-" not in M: M = M + "-0"
self._M_isauto, self._M_shift = True, int(M.split("-")[-1])
M = 0
if M < 0: raise RROMPyException("M must be non-negative.")
self._M = M
self._approxParameters["M"] = self.M
def _setMAuto(self):
if (self.rationalMode != "STANDARD_OUTPUT"
or self.functionalSolve[:11] == "BARYCENTRIC"):
SEff = self.S
else:
dim = self.HFEngine.outputdim
if hasattr(self, "_N_isauto") and self._N_isauto:
M, N = reduceDegreeMN(self.S - 2, self.S - 1, self.S,
self.npar, self.polydegreetype,
self.MAuxiliary, self.NAuxiliary, dim, 0)
self.M = max(0, M - self._M_shift)
self.N = max(0, N - self._N_shift)
vbMng(self, "MAIN",
"Automatically setting M to {} and N to {}.".format(
self.M, self.N), 25)
return
degN = degreeToDOFs([self.N] + [self.NAuxiliary] * (self.npar - 1),
- self.polydegreetype) - 1
+ self.polydegreetypeM) - 1
SEff = int(np.floor(self.S - degN / dim))
- M = reduceDegreeN(SEff - 1, SEff, self.npar, self.polydegreetype,
+ M = reduceDegreeN(SEff - 1, SEff, self.npar, self.polydegreetypeM,
self.MAuxiliary)
self.M = max(0, M - self._M_shift)
vbMng(self, "MAIN", "Automatically setting M to {}.".format(self.M),
25)
@property
def N(self):
"""Value of N."""
return self._N
@N.setter
def N(self, N):
if isinstance(N, str):
N = N.strip().replace(" ","")
if "-" not in N: N = N + "-0"
self._N_isauto, self._N_shift = True, int(N.split("-")[-1])
N = 0
if N < 0: raise RROMPyException("N must be non-negative.")
self._N = N
self._approxParameters["N"] = self.N
def _setNAuto(self):
if self.rationalMode[-6:] == "OUTPUT":
dim = self.HFEngine.outputdim
else:
dim = self.HFEngine.spacedim
if (self.rationalMode[:7] == "MINIMAL"
or self.functionalSolve[:11] == "BARYCENTRIC"):
if self.rationalMode[:7] == "MINIMAL":
SEff = min(self.S, dim + 1)
else: #if self.rationalMode[:7] == "STANDARD":
SEff = int(np.floor(dim * self.S / (dim + 1.))) + 1
else:
if hasattr(self, "_M_isauto") and self._M_isauto:
if self.rationalMode == "STANDARD_OUTPUT":
return self._setMAuto()
M = reduceDegreeN(self.S - 2, self.S - 1, self.npar,
- self.polydegreetype, self.MAuxiliary)
+ self.polydegreetypeM, self.MAuxiliary)
else:
M = reduceDegreeN(self.M, self.S - 1, self.npar,
- self.polydegreetype, self.MAuxiliary)
+ self.polydegreetypeM, self.MAuxiliary)
degM = degreeToDOFs([M] + [self.MAuxiliary] * (self.npar - 1),
- self.polydegreetype)
+ self.polydegreetypeN)
SEff = min(self.S, dim * (self.S - degM) + 1)
- N = reduceDegreeN(SEff - 1, SEff, self.npar, self.polydegreetype,
+ N = reduceDegreeN(SEff - 1, SEff, self.npar, self.polydegreetypeN,
self.NAuxiliary)
self.N = max(0, N - self._N_shift)
vbMng(self, "MAIN", "Automatically setting N to {}.".format(self.N),
25)
- @property
- def indexPrimary(self):
- """Value of indexPrimary."""
- return self._indexPrimary
- @indexPrimary.setter
- def indexPrimary(self, indexPrimary):
- if indexPrimary < 0 or indexPrimary > self.npar:
- raise RROMPyException("indexPrimary must be between 0 and npar.")
- self._indexPrimary = indexPrimary
- self._approxParameters["indexPrimary"] = self.indexPrimary
@property
def MAuxiliary(self):
"""Value of MAuxiliary."""
return self._MAuxiliary
@MAuxiliary.setter
def MAuxiliary(self, MAuxiliary):
if MAuxiliary < 0:
raise RROMPyException("MAuxiliary must be non-negative.")
self._MAuxiliary = MAuxiliary
self._approxParameters["MAuxiliary"] = self.MAuxiliary
@property
def NAuxiliary(self):
"""Value of NAuxiliary."""
return self._NAuxiliary
@NAuxiliary.setter
def NAuxiliary(self, NAuxiliary):
if NAuxiliary < 0:
raise RROMPyException("NAuxiliary must be non-negative.")
self._NAuxiliary = NAuxiliary
self._approxParameters["NAuxiliary"] = self.NAuxiliary
@property
def polydegreetype(self):
"""Value of polydegreetype."""
return self._polydegreetype
@polydegreetype.setter
def polydegreetype(self, polydegreetype):
try:
polydegreetype = polydegreetype.upper().strip().replace(" ","")
- if polydegreetype not in ["TOTAL", "TENSOR"]:
+ if "_" not in polydegreetype:
+ polydegreetype = "{0}_{0}".format(polydegreetype)
+ polydegreetypes = polydegreetype.split("_")
+ if not np.all([pdt in ["TOTAL", "TENSOR"]
+ for pdt in polydegreetypes]):
raise RROMPyException(("Prescribed polydegreetype not "
"recognized."))
except Exception as E:
- RROMPyWarning(str(E) + " Overriding to 'TOTAL'.")
- polydegreetype = "TOTAL"
+ RROMPyWarning(str(E) + " Overriding to 'TENSOR_TOTAL'.")
+ polydegreetype = "TENSOR_TOTAL"
self._polydegreetype = polydegreetype
self._approxParameters["polydegreetype"] = self.polydegreetype
+ @property
+ def polydegreetypeM(self): return self.polydegreetype.split("_")[0]
+
+ @property
+ def polydegreetypeN(self): return self.polydegreetype.split("_")[1]
+
@property
def forceQReal(self):
"""Value of forceQReal."""
return self._forceQReal
@forceQReal.setter
def forceQReal(self, forceQReal):
self._forceQReal = forceQReal
self._approxParameters["forceQReal"] = self.forceQReal
@property
def polyTruncateTol(self):
"""Value of tolerance for truncation of rational terms."""
return self._polyTruncateTol
@polyTruncateTol.setter
def polyTruncateTol(self, polyTruncateTol):
if polyTruncateTol < 0.:
RROMPyWarning(("Overriding prescribed negative numerator "
"tolerance to 0."))
polyTruncateTol = 0.
self._polyTruncateTol = polyTruncateTol
self._approxParameters["polyTruncateTol"] = self.polyTruncateTol
@property
def QTol(self):
"""Value of tolerance for robust rational denominator management."""
return self._QTol
@QTol.setter
def QTol(self, QTol):
if QTol < 0.:
RROMPyWarning(("Overriding prescribed negative robustness "
"tolerance to 0."))
QTol = 0.
self._QTol = QTol
self._approxParameters["QTol"] = self.QTol
def resetSamples(self):
"""Reset samples."""
super().resetSamples()
self._musUniqueCN = None
self._derIdxs = None
self._reorder = None
def _setupInterpolationIndices(self):
"""Setup parameters for polyvander."""
if self._musUniqueCN is None or len(self._reorder) != len(self.mus):
self._musUniqueCN, musIdxsTo, musIdxs, musCount = (
self.trainedModel.centerNormalize(self.mus).unique(
return_index = True, return_inverse = True,
return_counts = True))
self._musUnique = self.mus[musIdxsTo]
self._derIdxs = [None] * len(self._musUniqueCN)
self._reorder = np.empty(len(musIdxs), dtype = int)
filled = 0
for j, cnt in enumerate(musCount):
self._derIdxs[j] = nextDerivativeIndices([], self.mus.shape[1],
cnt)
jIdx = np.nonzero(musIdxs == j)[0]
self._reorder[jIdx] = np.arange(filled, filled + cnt)
filled += cnt
def _setupDenominator(self):
"""Compute rational denominator."""
RROMPyAssert(self._mode, message = "Cannot setup denominator.")
vbMng(self, "INIT", "Starting computation of denominator.", 7)
if self.functionalSolve != "NORM" and self.npar > 1:
RROMPyWarning(("Strategy for functional optimization must be "
"'NORM' for more than one parameter. Overriding to "
"'NORM'."))
self.functionalSolve = "NORM"
if (self.functionalSolve[:11] == "BARYCENTRIC"
and self.rationalMode[:7] == "MINIMAL"
and not hasattr(self, "_N_isauto") and self.N + 1 < self.S):
RROMPyWarning(("Minimal barycentric strategy cannot be applied "
"with Least Squares. Overriding to 'NORM'."))
self.functionalSolve = "NORM"
if self.functionalSolve[:11] == "BARYCENTRIC":
self._setupInterpolationIndices()
if len(self._musUnique) != self.S:
RROMPyWarning(("Barycentric functional optimization cannot be "
"applied to repeated samples. Overriding to "
"'NORM'."))
self.functionalSolve = "NORM"
if hasattr(self, "_N_isauto"): self._setNAuto()
if (self.rationalMode[:7] != "MINIMAL" and hasattr(self, "_M_isauto")
and self.functionalSolve[:11] != "BARYCENTRIC"):
self._setMAuto()
if self.rationalMode[-6:] == "OUTPUT":
dim = self.HFEngine.outputdim
else:
dim = self.HFEngine.spacedim
if (self.rationalMode[:7] == "MINIMAL"
or self.functionalSolve[:11] == "BARYCENTRIC"):
if self.rationalMode[:7] == "MINIMAL":
SEff = min(self.S, dim + 1)
else: #if self.rationalMode[:7] == "STANDARD":
SEff = int(np.floor(dim * self.S / (dim + 1.))) + 1
- N = reduceDegreeN(self.N, SEff, self.npar, self.polydegreetype,
+ N = reduceDegreeN(self.N, SEff, self.npar, self.polydegreetypeN,
self.NAuxiliary)
if N < self.N:
RROMPyWarning(("N too large compared to S. Reducing N by "
"{}.").format(self.N - N))
self.N = N
else:
N = reduceDegreeMN(self.M, self.N, self.S, self.npar,
self.polydegreetype, self.MAuxiliary,
self.NAuxiliary, dim)[1]
if N < self.N:
RROMPyWarning(("N too large compared to S. Reducing N by "
"{}.").format(self.N - N))
self.N = N
invD, TN = None, None
if self.functionalSolve[:11] != "BARYCENTRIC":
invD, TN = self._computeInterpolantInverseBlocks()
elif self.rationalMode[:7] != "MINIMAL":
musCN = self._musUniqueCN[self._reorder]
TN = np.empty((0, self.S + 1), dtype = np.complex)
Rscaling = None
if self.rationalMode[-6:] == "OUTPUT":
sampleE = self.samplingEngine.samples_output
elif self.POD == 1:
sampleE = self.samplingEngine.Rscale
elif self.POD == 1/2:
sampleE = self.samplingEngine.samples_normal
Rscaling = self.samplingEngine.Rscale
else:
sampleE = self.samplingEngine.samples
while self.N > 0:
if self.functionalSolve[:11] != "BARYCENTRIC":
NEff = degreeToDOFs([self.N]
+ [self.NAuxiliary] * (self.npar - 1),
- self.polydegreetype)
+ self.polydegreetypeN)
TN = TN[:, : NEff]
elif self.rationalMode[:7] != "MINIMAL":
NOld = TN.shape[1] - 1
dTN = (musCN[self.N : NOld]
- musCN[: self.N].reshape(1, -1)) ** -1
dTN = np.pad(dTN, [(0, 0), (1, 0)], "constant",
constant_values = 1)
TN = np.vstack((TN[:, : self.N + 1], dTN))
if self.rationalMode[:7] == "MINIMAL":
ev, eV = self.findeveVGMinimal(sampleE, invD, TN, Rscaling)
else: #if self.rationalMode[: 8] == "STANDARD":
ev, eV = self.findeveVGStandard(sampleE, invD, TN, Rscaling)
if (self.rationalMode[:7] == "MINIMAL"
and self.functionalSolve[:11] == "BARYCENTRIC"):
break
nevBad = np.sum(np.abs(ev / ev[-1]) < self.QTol)
if not nevBad: break
dN = int(np.ceil(nevBad / (1 + self.NAuxiliary)
** (self.npar - 1)))
vbMng(self, "MAIN",
("Smallest {} eigenvalue{} below tolerance. Reducing N by "
"{}.").format(nevBad, "s" * (nevBad > 1), dN), 10)
self.N = self.N - dN
if self.rationalMode[-5:] == "STATE" and self.POD != 1: del self._gram
if self.N <= 0:
self.N = 0
- deg = [self.NAuxiliary + 1] * self.npar
- deg[self.indexPrimary] = 1
- eV = np.zeros(deg, dtype = np.complex)
+ eV = np.zeros((1,) + (self.NAuxiliary + 1,) * (self.npar - 1),
+ dtype = np.complex)
eV[(0,) * self.npar] = 1.
if self.N > 0 and self.functionalSolve[:11] == "BARYCENTRIC":
q = PIN()
q.polybasis, q.nodes = self.polybasis, eV
else:
q = PI()
q.npar, q.polybasis = self.npar, self.polybasis
- deg = [self.NAuxiliary + 1] * self.npar
- deg[self.indexPrimary] = self.N + 1
- if self.polydegreetype == "TENSOR":
+ deg = [self.N + 1] + [self.NAuxiliary + 1] * (self.npar - 1)
+ if self.polydegreetypeN == "TENSOR":
q.coeffs = eV.reshape(deg)
- else: #if self.polydegreetype == "TOTAL":
+ else: #if self.polydegreetypeN == "TOTAL":
q.coeffs = mapTotalToTensor(deg, self.npar, eV)
q.truncate(self.polyTruncateTol)
vbMng(self, "DEL", "Done computing denominator.", 7)
return q
def _setupNumerator(self):
"""Compute rational numerator."""
RROMPyAssert(self._mode, message = "Cannot setup numerator.")
vbMng(self, "INIT", "Starting computation of numerator.", 7)
self._setupInterpolationIndices()
Qevaldiag = polyTimesTable(self.trainedModel.data.Q, self._musUniqueCN,
self._reorder, self._derIdxs,
self.scaleFactorRel)
if self.rationalMode[-6:] == "OUTPUT":
Qevaldiag = Qevaldiag.dot(self.samplingEngine.samples_output.T)
else:
if self.POD == 1:
Qevaldiag = Qevaldiag.dot(self.samplingEngine.Rscale.T)
elif self.POD == 1/2:
Qevaldiag = Qevaldiag * self.samplingEngine.Rscale
if hasattr(self, "_M_isauto"): self._setMAuto()
- M = reduceDegreeN(self.M, self.S, self.npar, self.polydegreetype,
+ M = reduceDegreeN(self.M, self.S, self.npar, self.polydegreetypeM,
self.MAuxiliary)
if M < self.M:
RROMPyWarning(("M too large compared to S. Reducing M by "
"{}").format(self.M - M))
self.M = M
while self.M >= 0:
- deg = [self.MAuxiliary] * self.npar
- deg[self.indexPrimary] = self.M
+ deg = [self.M] + [self.MAuxiliary] * (self.npar - 1)
pParRest = [deg, self.polybasis, self.verbosity >= 5,
- self.polydegreetype, self.polyTruncateTol,
+ self.polydegreetypeM, self.polyTruncateTol,
{"derIdxs": self._derIdxs, "reorder": self._reorder,
"scl": self.scaleFactorRel,
- "degreetype": self.polydegreetype},
+ "degreetype": self.polydegreetypeM},
{"rcond": self.interpTol}]
p = PI()
wellCond, msg = p.setupByInterpolation(self._musUniqueCN,
Qevaldiag, *pParRest)
vbMng(self, "MAIN", msg, 5)
if wellCond: break
vbMng(self, "MAIN", ("Polyfit is poorly conditioned. Reducing M "
"by 1."), 10)
self.M = self.M - 1
if self.M < 0:
raise RROMPyException(("Instability in computation of numerator. "
"Aborting."))
vbMng(self, "DEL", "Done computing numerator.", 7)
return p
def setupApprox(self) -> int:
"""Compute rational interpolant."""
if self.checkComputedApprox(): return -1
RROMPyAssert(self._mode, message = "Cannot setup approximant.")
vbMng(self, "INIT", "Setting up {}.".format(self.name()), 5)
self.computeSnapshots()
if self.rationalMode[-6:] == "OUTPUT":
vbMng(self, "INIT", "Extracting system output from state.", 35)
self.samplingEngine.samples_output = self.HFEngine.applyC(
self.samplingEngine.samples, self.mus)
_POD, self._POD = self.POD, 1
vbMng(self, "DEL", "Done extracting system output.", 35)
self._setupTrainedModel(self.samplingEngine.projectionMatrix,
collapsed = (self.rationalMode[-6:] == "OUTPUT"))
self._setupRational(self._setupDenominator())
self.trainedModel.data.approxParameters = copy(self.approxParameters)
if self.rationalMode[-6:] == "OUTPUT": self._POD = _POD
vbMng(self, "DEL", "Done setting up approximant.", 5)
return 0
def _setupRational(self, Q:interpEng, P : interpEng = None):
vbMng(self, "INIT", "Starting approximant finalization.", 5)
self.trainedModel.data.Q = Q
if P is None: P = self._setupNumerator()
while self.N > 0 and self.npar == 1:
if self.HFEngine._ignoreResidues:
pls = Q.roots()
cfs, projMat = None, None
else:
cfs, pls, _ = rational2heaviside(P, Q)
cfs = cfs[: len(pls)].T
if self.POD != 1:
projMat = self.samplingEngine.projectionMatrix
else:
projMat = None
foci = self.sampler.normalFoci()
plsA = self.mapParameterList(self.mapParameterList(self.mu0)(0, 0)
+ self.scaleFactor * pls, "B")(0)
idxBad = self.HFEngine.flagBadPolesResiduesAbsolute(plsA, cfs,
projMat)
if not self.HFEngine._ignoreResidues: cfs[:, idxBad] = 0.
idxBad += self.HFEngine.flagBadPolesResiduesRelative(pls, cfs,
projMat, foci)
idxBad = idxBad > 0
if not np.any(idxBad): break
vbMng(self, "MAIN",
"Removing {} spurious pole{} out of {}.".format(
np.sum(idxBad), "s" * (np.sum(idxBad) > 1), self.N), 10)
if isinstance(Q, PIN):
Q.nodes = Q.nodes[idxBad == False]
else:
Q = PI()
Q.npar, Q.polybasis = 1, self.polybasis
Q.coeffs = np.ones(1, dtype = np.complex)
for pl in pls[idxBad == False]:
Q.coeffs = polyTimes(Q.coeffs, [- pl, 1.],
Pbasis = Q.polybasis, Rbasis = Q.polybasis)
Q.coeffs /= np.linalg.norm(Q.coeffs)
self.trainedModel.data.Q = Q
self.N = Q.deg[0]
P = self._setupNumerator()
self.trainedModel.data.P = P
vbMng(self, "DEL", "Terminated approximant finalization.", 5)
def _computeInterpolantInverseBlocks(self) -> Tuple[List[Np2D], Np2D]:
"""
Compute inverse factors for rational interpolant target functional.
"""
RROMPyAssert(self._mode, message = "Cannot solve eigenvalue problem.")
self._setupInterpolationIndices()
pvPPar = [self.polybasis, self._derIdxs, self._reorder,
- self.scaleFactorRel, self.polydegreetype]
+ self.scaleFactorRel, self.polydegreetypeN]
lagrange = self.npar == 1 and self.S == len(self._musUniqueCN)
if self.rationalMode[:7] == "MINIMAL" and not lagrange:
E = reduceDegreeN(self.S - 1, self.S, self.npar,
- self.polydegreetype, self.NAuxiliary)
+ self.polydegreetypeN, self.NAuxiliary)
TN = None
deg = [self.NAuxiliary] * self.npar
while self.N >= 0:
if TN is None:
- deg[self.indexPrimary] = self.N
+ deg[0] = self.N
TE = TN = pvP(self._musUniqueCN, deg, *pvPPar)
if self.rationalMode[:7] == "MINIMAL" and not lagrange:
- deg[self.indexPrimary] = E
+ deg[0] = E
if self.N == E:
TE = TN
else:
TE = pvP(self._musUniqueCN, deg, *pvPPar)
fitOut = pseudoInverse(TE, rcond = self.interpTol,
full = True)
degE = deg[0] if self.npar == 1 else deg
vbMng(self, "MAIN",
("Fitting {} samples with degree {} through {}... "
"Conditioning of pseudoinverse system:"
"{:.4e}.").format(TE.shape[0], degE,
polyfitname(self.polybasis),
fitOut[1][1][0] / fitOut[1][1][-1]), 5)
if fitOut[1][0] == TE.shape[1]:
if self.rationalMode[:7] == "MINIMAL" and not lagrange:
- if self.polydegreetype == "TENSOR":
- fitinv = fitOut[0][fullDegreeMaxMask(deg,
- self.indexPrimary)]
- else: #if self.polydegreetype == "TOTAL":
- deg0 = [E - 1] + [self.NAuxiliary] * (self.npar - 1)
- fitinv = fitOut[0][degreeToDOFs(deg0, "TOTAL") :]
+ deg0 = [E - 1] + [self.NAuxiliary] * (self.npar - 1)
+ S0 = degreeToDOFs(deg0, self.polydegreetypeN)
+ fitinv = fitOut[0][S0 :]
break
if self.rationalMode[:7] == "MINIMAL" and not lagrange:
E -= 1
if E >= self.N: continue
vbMng(self, "MAIN", ("Polyfit is poorly conditioned. Reducing "
"N by 1."), 10)
self.N, TN = self.N - 1, None
if self.N < 0:
raise RROMPyException(("Instability in computation of "
"denominator. Aborting."))
if self.rationalMode[:7] == "MINIMAL":
if lagrange:
mus = self._musUniqueCN[self._reorder]
dist = baseDistanceMatrix(mus, magnitude = False)[..., 0]
dist[np.arange(self.S), np.arange(self.S)] = 1.
fitinvE = np.prod(dist, axis = 1) ** -1
vbMng(self, "MAIN",
("Evaluating quasi-Lagrangian basis of degree {} at {} "
"sample points.").format(self.S - 1, self.S), 5)
invD = [np.diag(fitinvE)]
else:
invD = vanderInvTable(fitinv, self._reorder, self._derIdxs)
else: #if self.rationalMode[: 8] == "STANDARD":
if self.rationalMode[-6:] == "OUTPUT":
dim = self.HFEngine.outputdim
else:
dim = self.HFEngine.spacedim
degN = degreeToDOFs([self.N] + [self.NAuxiliary] * (self.npar - 1),
- self.polydegreetype) - 1
+ self.polydegreetypeN) - 1
SEff = int(np.floor(self.S - degN / dim))
MEff = reduceDegreeN(SEff - 1, SEff, self.npar,
- self.polydegreetype, self.MAuxiliary)
+ self.polydegreetypeM, self.MAuxiliary)
deg = [self.MAuxiliary] * self.npar
while MEff >= 0:
- deg[self.indexPrimary] = MEff
+ deg[0] = MEff
if MEff == self.N and self.MAuxiliary == self.NAuxiliary:
TM = TN
else:
TM = pvP(self._musUniqueCN, deg, *pvPPar)
Q, R = np.linalg.qr(TM, 'complete')
fitOut = pseudoInverse(R[: R.shape[1]], rcond = self.interpTol,
full = True)
degM = deg[0] if self.npar == 1 else deg
vbMng(self, "MAIN",
("Fitting {} samples with degree {} through {}... "
"Conditioning of pseudoinverse system:"
"{:.4e}.").format(TM.shape[0], degM,
polyfitname(self.polybasis),
fitOut[1][1][0] / fitOut[1][1][-1]),
15)
if fitOut[1][0] == R.shape[1]: break
MEff -= 1
Qo = Q[:, TM.shape[1] :].T.conj()
vbMng(self, "MAIN",
("Loss of orthogonality for adjoint vandermonde: "
"{:.4e}.").format(np.linalg.norm(Qo.dot(TM))), 5)
invD = [Qo]
return invD, TN
def findeveVGMinimal(self, sampleE:Np2D, invD:List[Np2D], TN:Np2D,
Rscaling : Np1D = None) -> Tuple[Np1D, Np2D]:
"""
Compute eigenvalues and eigenvectors of rational denominator matrix, or
of its right chol factor if POD.
"""
RROMPyAssert(self._mode, message = "Cannot solve spectral problem.")
if self.POD == 1 or self.rationalMode == "MINIMAL_OUTPUT":
if self.functionalSolve[:11] == "BARYCENTRIC":
Rstack = sampleE
else:
vbMng(self, "INIT", "Building generalized half-gramian.", 10)
S, eWidth = sampleE.shape[0], len(invD)
Rstack = np.zeros((S * eWidth, TN.shape[1]),
dtype = np.complex)
for k in range(eWidth):
Rstack[k * S : (k + 1) * S, :] = dot(sampleE,
dot(invD[k], TN))
vbMng(self, "DEL", "Done building half-gramian.", 10)
if self.forceQReal:
Rstack = np.vstack((np.real(Rstack), np.imag(Rstack)))
_, s, Vh = np.linalg.svd(Rstack,
full_matrices = Rstack.shape[0] < Rstack.shape[1])
evG, eVG = s[::-1], Vh[::-1].T.conj()
if len(evG) < eVG.shape[1]:
evG = np.append(np.zeros(eVG.shape[1] - len(evG)), evG)
evExp, probKind = -2., "svd "
else:
if not hasattr(self, "_gram"):
vbMng(self, "INIT", "Building gramian matrix.", 10)
self._gram = self.HFEngine.innerProduct(sampleE, sampleE,
is_state = True)
if Rscaling is not None:
self._gram = (self._gram.T * Rscaling.conj()).T * Rscaling
vbMng(self, "DEL", "Done building gramian.", 10)
if self.functionalSolve[:11] == "BARYCENTRIC":
G = self._gram
else:
vbMng(self, "INIT", "Building generalized gramian.", 10)
G = np.zeros((TN.shape[1],) * 2, dtype = np.complex)
for k in range(len(invD)):
iDkN = dot(invD[k], TN)
G += dot(dot(self._gram, iDkN).T, iDkN.conj()).T
vbMng(self, "DEL", "Done building gramian.", 10)
if self.forceQReal: G = np.real(G)
evG, eVG = np.linalg.eigh(G)
evExp, probKind = -1., "eigen"
if (self.functionalSolve in ["NORM", "BARYCENTRIC_NORM"]
or np.sum(np.abs(evG) < np.finfo(float).eps * np.abs(evG[-1])
* len(evG)) == 1):
eV = eVG[:, 0]
elif self.functionalSolve == "BARYCENTRIC_AVERAGE":
eV = eVG.dot(evG ** evExp * np.sum(eVG, axis = 0).conj())
else:
eV = eVG.dot(evG ** evExp * eVG[0].conj())
if self.functionalSolve[:11] == "BARYCENTRIC": _evG = np.array(evG)
evG[1 :] -= evG[0]
eV /= np.linalg.norm(eV)
+ if self.npar == 1:
+ degE = self.N
+ else:
+ degE = [self.N] + [self.NAuxiliary] * (self.npar - 1)
cond = evG[-1] / evG[1] if evG[0] > 0 and evG[1] > 0 else np.inf
vbMng(self, "MAIN",
- ("Solved {}problem of size {} with condition number "
- "{:.4e}.").format(probKind, len(evG) - 1, cond), 5)
+ ("Solved {}problem of size {} (degree {}) with condition number "
+ "{:.4e}.").format(probKind, len(evG) - 1, degE, cond), 5)
if self.functionalSolve[:11] == "BARYCENTRIC":
S, mus = len(eV), self._musUniqueCN[self._reorder].flatten()
poles, qTm1 = barycentricRoots(eV, mus, True)
self.N = len(poles)
if self.QTol > 0:
# compare optimal score with N poles to those obtained by
# removing one of the poles
qTm1 = qTm1[1 :].conj() ** -1.
dists = mus.reshape(-1, 1) - mus
dists[np.arange(S), np.arange(S)] = 1.
dists = np.prod(dists, axis = 1).conj() ** -1.
qComp = np.empty((self.N + 1, S), dtype = np.complex)
qComp[0] = dists * np.prod(qTm1, axis = 1)
for j in range(self.N):
qTmj = np.prod(qTm1[:, np.arange(self.N) != j], axis = 1)
qComp[j + 1] = dists * qTmj
Lqs = qComp.dot(eVG)
scores = np.real(np.sum(Lqs * _evG ** -evExp * Lqs.conj(),
axis = 1))
evBad = scores[1 :] < self.QTol * scores[0]
nevBad = np.sum(evBad)
if nevBad:
vbMng(self, "MAIN",
("Suboptimal pole{} detected. Reducing N by "
"{}.").format("s" * (nevBad > 1), nevBad), 10)
self.N = self.N - nevBad
poles = poles[evBad == False]
eV = poles
return evG[1 :], eV
def findeveVGStandard(self, sampleE:Np2D, invD:List[Np2D], TN:Np2D,
Rscaling : Np1D = None) -> Tuple[Np1D, Np2D]:
"""
Compute eigenvalues and eigenvectors of rational denominator matrix, or
of its right chol factor if POD.
"""
RROMPyAssert(self._mode, message = "Cannot solve spectral problem.")
if self.POD == 1 or self.rationalMode == "STANDARD_OUTPUT":
vbMng(self, "INIT", "Building generalized half-gramian.", 10)
if self.functionalSolve[:11] == "BARYCENTRIC":
Rold = np.pad(sampleE[: self.N, : self.N], [(0, 0), (1, 0)],
"constant") # adjust for bias
Rnew = sampleE[: self.N, self.N :]
S, eWidth = Rnew.shape
Rstack = np.zeros((S * eWidth, self.N + 1), dtype = np.complex)
for k in range(eWidth):
Gloc = (Rold.T - Rnew[:, k]).T
Rstack[k * S : (k + 1) * S, :] = Gloc * TN[k]
else:
Rstack = khatri_rao(sampleE, invD[0]).dot(TN)
vbMng(self, "DEL", "Done building half-gramian.", 10)
if self.forceQReal:
Rstack = np.vstack((np.real(Rstack), np.imag(Rstack)))
_, s, Vh = np.linalg.svd(Rstack,
full_matrices = Rstack.shape[0] < Rstack.shape[1])
evG, eVG = s[::-1], Vh[::-1].T.conj()
if len(evG) < eVG.shape[1]:
evG = np.append(np.zeros(eVG.shape[1] - len(evG)), evG)
evExp, probKind = -2., "svd "
else:
if not hasattr(self, "_gram"):
vbMng(self, "INIT", "Building gramian matrix.", 10)
self._gram = self.HFEngine.innerProduct(sampleE, sampleE,
is_state = True)
if Rscaling is not None:
self._gram = (self._gram.T * Rscaling.conj()).T * Rscaling
vbMng(self, "DEL", "Done building gramian.", 10)
vbMng(self, "INIT", "Building generalized gramian.", 10)
if self.functionalSolve[:11] == "BARYCENTRIC":
Idiag = np.arange(self.N, self._gram.shape[0])
GNN = self._gram[Idiag, Idiag]
GON = np.pad(self._gram[: self.N, self.N :], [(1, 0), (0, 0)],
"constant") # adjust for bias
GNO = GON.T.conj()
GOO = np.pad(self._gram[: self.N, : self.N], (1, 0),
"constant") # adjust for bias
G = ((TN.T.conj() * GNN).dot(TN) #L^H @ diag(GNN) @ L
- (GON * TN.T.conj()).dot(TN) #(L^H * GON) @ L
- TN.T.conj().dot(GNO * TN) #L^H @ (GNO * L)
+ GOO * (TN.T.conj().dot(TN))) #GOO * (L^H @ L)
else:
G = self._gram * invD[0].T.conj().dot(invD[0])
G = TN.T.conj().dot(G.dot(TN))
vbMng(self, "DEL", "Done building gramian.", 10)
if self.forceQReal: G = np.real(G)
evG, eVG = np.linalg.eigh(G)
evExp, probKind = -1., "eigen"
if (self.functionalSolve in ["NORM", "BARYCENTRIC_NORM"]
or np.sum(np.abs(evG) < np.finfo(float).eps * np.abs(evG[-1])
* len(evG)) == 1):
eV = eVG[:, 0]
else:
eV = eVG.dot(evG ** evExp * eVG[0].conj())
evG[1 :] -= evG[0]
eV /= np.linalg.norm(eV)
+ if self.npar == 1:
+ degE = self.N
+ else:
+ degE = [self.N] + [self.NAuxiliary] * (self.npar - 1)
cond = evG[-1] / evG[1] if evG[0] > 0 and evG[1] > 0 else np.inf
vbMng(self, "MAIN",
- ("Solved {}problem of size {} with condition number "
- "{:.4e}.").format(probKind, len(evG) - 1, cond), 5)
+ ("Solved {}problem of size {} (degree {}) with condition number "
+ "{:.4e}.").format(probKind, len(evG) - 1, degE, cond), 5)
if self.functionalSolve[:11] == "BARYCENTRIC":
mus = self._musUniqueCN[self._reorder][: self.N, 0]
eV = barycentricRoots(eV, mus)
self.N = len(eV)
return evG[1 :], eV
def getResidues(self, *args, **kwargs) -> Tuple[paramList, Np2D]:
"""
Obtain approximant residues.
Returns:
Matrix with residues as columns.
"""
return self.trainedModel.getResidues(*args, **kwargs)
diff --git a/rrompy/reduction_methods/standard/trained_model/trained_model_reduced_basis.py b/rrompy/reduction_methods/standard/trained_model/trained_model_reduced_basis.py
index 905cd72..52761f4 100644
--- a/rrompy/reduction_methods/standard/trained_model/trained_model_reduced_basis.py
+++ b/rrompy/reduction_methods/standard/trained_model/trained_model_reduced_basis.py
@@ -1,150 +1,149 @@
# Copyright (C) 2018-2020 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 collections.abc import Iterable
from rrompy.reduction_methods.base.trained_model.trained_model import (
TrainedModel)
from rrompy.reduction_methods.standard.reduced_basis_utils import (
projectAffineDecomposition)
from rrompy.utilities.base.types import ListAny, paramVal, paramList, sampList
from rrompy.utilities.base import verbosityManager as vbMng, freepar as fp
from rrompy.utilities.numerical.compress_matrix import compressMatrix
from rrompy.utilities.numerical.marginalize_poly_list import (
marginalizePolyList)
from rrompy.utilities.numerical.nonlinear_eigenproblem import (
eigvalsNonlinearDense)
from rrompy.utilities.expression import expressionEvaluator
from rrompy.utilities.exception_manager import (RROMPyException, RROMPyAssert,
RROMPyWarning)
from rrompy.parameter import checkParameter
from rrompy.sampling import sampleList
from rrompy.utilities.parallel import (poolRank, masterCore, listScatter,
matrixGatherv, isend, recv)
__all__ = ['TrainedModelReducedBasis']
class TrainedModelReducedBasis(TrainedModel):
"""
ROM approximant evaluation for RB approximant.
Attributes:
Data: dictionary with all that can be pickled.
"""
def reset(self):
super().reset()
if hasattr(self, "data") and hasattr(self.data, "lastSetupMu"):
self.data.lastSetupMu = None
def compress(self, collapse : bool = False, tol : float = 0., *args,
**kwargs):
if collapse:
raise RROMPyException("Cannot collapse implicit surrogates.")
if tol <= 0.: return
super().compress(collapse, tol)
self.data.projMat, RMat, _ = compressMatrix(self.data.projMat, tol,
*args, **kwargs)
self.data.ARBs, self.data.bRBs = projectAffineDecomposition(
self.data.ARBs, self.data.bRBs, RMat)
def assembleReducedModel(self, mu:paramVal):
mu = checkParameter(mu, self.data.npar)
if not (hasattr(self.data, "lastSetupMu")
and self.data.lastSetupMu == mu):
vbMng(self, "INIT", "Assembling reduced model for mu = {}."\
.format(mu), 17)
muEff = self.mapParameterList(mu)
self.data.ARBmu, self.data.bRBmu = 0., 0.
for thA, ARB in zip(self.data.thAs, self.data.ARBs):
self.data.ARBmu = (expressionEvaluator(thA[0], muEff) * ARB
+ self.data.ARBmu)
for thb, bRB in zip(self.data.thbs, self.data.bRBs):
self.data.bRBmu = (expressionEvaluator(thb[0], muEff) * bRB
+ self.data.bRBmu)
vbMng(self, "DEL", "Done assembling reduced model.", 17)
self.data.lastSetupMu = mu
def getApproxReduced(self, mu : paramList = []) -> sampList:
"""
Evaluate reduced representation of approximant at arbitrary parameter.
Args:
mu: Target parameter.
"""
mu = self.checkParameterList(mu)
if (not hasattr(self, "lastSolvedApproxReduced")
or self.lastSolvedApproxReduced != mu):
vbMng(self, "INIT",
"Computing RB solution at mu = {}.".format(mu), 12)
+ req, is_master = [], masterCore()
mu, _, sizes = listScatter(mu, return_sizes = True)
mu = self.checkParameterList(mu)
- req, emptyCores = [], np.where(sizes == 0)[0]
+ emptyCores = np.where(np.array(sizes) == 0)[0] if is_master else []
if len(mu) == 0:
vbMng(self, "MAIN", "Idling.", 37)
uL = recv(source = 0, tag = poolRank())
uApproxR = np.empty((uL, 0), dtype = np.complex)
else:
for j, mj in enumerate(mu):
self.assembleReducedModel(mj)
uAppR = np.linalg.solve(self.data.ARBmu, self.data.bRBmu)
if j == 0:
uApproxR = np.empty((len(uAppR), len(mu)),
dtype = np.complex)
- if masterCore():
- for dest in emptyCores:
- req += [isend(len(uAppR), dest = dest,
- tag = dest)]
+ for dest in emptyCores:
+ req += [isend(len(uAppR), dest = dest, tag = dest)]
uApproxR[:, j] = uAppR
for r in req: r.wait()
uApproxR = matrixGatherv(uApproxR, sizes)
self.uApproxReduced = sampleList(uApproxR)
vbMng(self, "DEL", "Done computing RB solution.", 12)
self.lastSolvedApproxReduced = mu
return self.uApproxReduced
def getPoles(self, marginalVals : ListAny = [fp], jSupp : int = 1,
**kwargs) -> paramList:
"""
Obtain approximant poles.
Returns:
Numpy complex vector of poles.
"""
RROMPyAssert(self.data.npar, len(marginalVals), "Number of parameters")
if not self.data.affinePoly:
RROMPyWarning(("Unable to compute approximate poles due "
"to parametric dependence (detected non-"
"polynomial). Change HFEngine.affinePoly to True "
"if necessary."))
return
if not isinstance(marginalVals, Iterable):
marginalVals = [marginalVals]
mVals = list(marginalVals)
rDim = mVals.index(fp)
if rDim < len(mVals) - 1 and fp in mVals[rDim + 1 :]:
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(0, rDim)
mVals = checkParameter(mVals, return_data = True).flatten()
mVals[rDim] = fp
ARBs = marginalizePolyList(ARBs, mVals, "auto")
ev = eigvalsNonlinearDense(ARBs, jSupp = jSupp, **kwargs)
return self.mapParameterList(ev, "B", [rDim])(0)
diff --git a/rrompy/utilities/numerical/degree.py b/rrompy/utilities/numerical/degree.py
index 4ca24b8..7059a83 100644
--- a/rrompy/utilities/numerical/degree.py
+++ b/rrompy/utilities/numerical/degree.py
@@ -1,116 +1,117 @@
# Copyright (C) 2018-2020 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 binom
from rrompy.utilities.base.types import Np2D, Tuple, List
from .kroneckerer import kroneckerer
from rrompy.utilities.exception_manager import RROMPyException
__all__ = ['degreeToDOFs', 'reduceDegreeMN', 'reduceDegreeN', 'degreeSet',
'mapTotalToTensor']
-def findPrimaryDegree(degs:List[int]) -> int:
+def checkPrimaryDegree(degs:List[int]):
idx = np.where(np.array(degs) != degs[0])[0]
- if len(idx) == 0 or len(idx) == len(degs) - 1: return 0
- if len(idx) == 1: return idx[0]
- raise RROMPyException("Must have at most one different degree.")
+ if len(idx) not in [0, len(degs) - 1]:
+ raise RROMPyException("Primary degree must be first entry.")
def degreeToDOFs(degs:List[int], degtype:str) -> int:
npar = len(degs)
if npar < 1: return 0
if npar < 2: return 1 + degs[0] * npar
- index = findPrimaryDegree(degs)
- deg = degs[index]
- degS = degs[0] if index else degs[1]
+ checkPrimaryDegree(degs)
+ deg, degS = degs[0], degs[1]
degtype = degtype.upper().strip().replace(" ","")
if degtype == "TENSOR": return (1 + deg) * (1 + degS) ** (npar - 1)
if degtype == "TOTAL":
tot = 0
for j in range(min(npar, deg // max(1, degS + 1) + 1)):
tot += ((1 - 2 * (j % 2)) * int(binom(npar - 1, j))
* int(binom(deg - (degS + 1) * j + npar, npar)))
return tot
raise RROMPyException("Degree type not recognized.")
def reduceDegreeMN(degM:int, degN:int, S:int, npar:int, degtype:str,
degMS : int = None, degNS : int = None,
weight : float = 1., Mfirst : bool = 1) -> int:
+ if "_" not in degtype: degtype = "{0}_{0}".format(degtype)
+ degtypeM, degtypeN = degtype.split("_")
if degMS is None: degMS = degM
if degNS is None: degNS = degN
degsM, degsN = [degMS] * (npar - 1), [degNS] * (npar - 1)
- dM = degreeToDOFs([degM] + degsM, degtype) - 1
- dN = degreeToDOFs([degN] + degsN, degtype) - 1
+ dM = degreeToDOFs([degM] + degsM, degtypeM) - 1
+ dN = degreeToDOFs([degN] + degsN, degtypeN) - 1
while degM > 0 or degN > 0:
if Mfirst:
if S >= dM + dN / weight + 1: break
if degM > 0:
degM -= 1
- dM = degreeToDOFs([degM] + degsM, degtype) - 1
+ dM = degreeToDOFs([degM] + degsM, degtypeM) - 1
Mfirst = 1
if S >= dM + dN / weight + 1: break
if degN > 0:
degN -= 1
- dN = degreeToDOFs([degN] + degsN, degtype) - 1
+ dN = degreeToDOFs([degN] + degsN, degtypeN) - 1
return degM, degN
def reduceDegreeN(deg:int, S:int, npar:int, degtype:str,
degS : int = None) -> int:
if degS is None: degS = deg
degs = [degS] * (npar - 1)
while deg > 0 and S < degreeToDOFs([deg] + degs, degtype): deg -= 1
return deg
def degreeSet(degs:List[int], degtype:str,
return_mask : bool = False) -> List[List[int]]:
npar = len(degs)
+ checkPrimaryDegree(degs)
degtype = degtype.upper().strip().replace(" ","")
if degtype == "TENSOR":
idxs = np.empty((degreeToDOFs(degs, degtype), npar), dtype = int)
for j in range(npar):
nL = degreeToDOFs(degs[: j], "TENSOR") if j else 1
nR = degreeToDOFs(degs[j + 1 :], "TENSOR") if npar - 1 - j else 1
idxs[:, j] = kroneckerer(np.arange(degs[j] + 1), nL, nR)
if return_mask: return idxs, np.ones(len(idxs), dtype = bool)
return idxs
if degtype == "TOTAL":
idxs, mask = degreeSet(degs, "TENSOR", True)
if npar > 1:
- mask = np.sum(idxs, axis = 1) <= degs[findPrimaryDegree(degs)]
+ mask = np.sum(idxs, axis = 1) <= degs[0]
idxs = idxs[mask]
if return_mask: return idxs, mask
return idxs
raise RROMPyException("Degree type not recognized.")
def mapTotalToTensor(shapeTensor:Tuple[int], npar:int, coeffs:Np2D) -> Np2D:
if npar < 2: return coeffs.reshape(shapeTensor)
- index = findPrimaryDegree(shapeTensor[: npar])
+ checkPrimaryDegree(shapeTensor[: npar])
from .hash_derivative import hashIdxToDerivative as hashI
- degS = shapeTensor[0] - 1 if index else shapeTensor[1] - 1
- mask = np.arange(npar) != index
+ degS = shapeTensor[1] - 1
+ mask = np.arange(npar) != 0
full = np.zeros(shapeTensor, dtype = coeffs.dtype)
ij, j = 0, 0
while j < len(coeffs):
idx = np.array(hashI(ij, npar))
- if np.sum(idx) >= shapeTensor[index]:
+ if np.sum(idx) >= shapeTensor[0]:
raise RROMPyException("Too many coefficients for prescribed size.")
if np.max(idx[mask]) <= degS:
full[tuple(idx)] = coeffs[j]
j += 1
ij += 1
return full
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 799ead6..b600a40 100644
--- a/rrompy/utilities/poly_fitting/heaviside/heaviside_to_from_rational.py
+++ b/rrompy/utilities/poly_fitting/heaviside/heaviside_to_from_rational.py
@@ -1,75 +1,74 @@
# Copyright (C) 2018-2020 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, PolynomialInterpolatorNodal)
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, Tuple, interpEng
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,
truncateTol : float = 0) -> Tuple[interpEng, interpEng]:
basisN = basis.split("_")[0]
if basisD is None: basisD = basisN
num = PolynomialInterpolator()
den = PolynomialInterpolatorNodal()
den.polybasis, den.nodes = basisD, poles
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, :]]
xAux = np.linspace(murange[0], murange[1], len(c) * 2)
valsAux = den(xAux) * polyval(xAux, c, poles, basis)
num.setupByInterpolation(xAux, valsAux, len(c) - 1, basisN,
truncateTol = truncateTol)
return num, den
def rational2heaviside(num:interpEng, den:interpEng,
murange : Np1D = np.array([-1., 1.]), scl : Np1D = None,
badPoleTol : float = -1) -> 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"
poles = den.roots()
if badPoleTol > 0.: poles = poles[np.abs(den(poles)) < badPoleTol]
c = np.zeros((max(len(poles), len(num.coeffs)),) + num.coeffs.shape[1 :],
dtype = np.complex)
if len(poles) > 0: c[: len(poles)] = (num(poles) / den(poles, 1, scl)).T
if len(c) > len(poles):
xAux = np.linspace(murange[0], murange[1], len(c) * 2)
valsAux = num(xAux) / den(xAux)
- if len(poles) > 0:
- valsAux -= polyval(xAux, c, poles, basis)
+ if len(poles) > 0: valsAux = valsAux - polyval(xAux, c, poles, basis)
VanAux = polyvander(xAux, [len(c) - len(poles) - 1], num.polybasis)
c[len(poles) :] = customFit(VanAux, valsAux.T)
return c, poles, basis