diff --git a/examples/5_anisotropic_square/anisotropic_square.py b/examples/5_anisotropic_square/anisotropic_square.py
index 136baa4..2dd8934 100644
--- a/examples/5_anisotropic_square/anisotropic_square.py
+++ b/examples/5_anisotropic_square/anisotropic_square.py
@@ -1,80 +1,80 @@
### example from Smetana, Zahm, Patera. Randomized residual-based error
### estimators for parametrized equations.
import numpy as np
import matplotlib.pyplot as plt
from itertools import product
from anisotropic_square_engine import (AnisotropicSquareEngine as engine,
AnisotropicSquareEnginePoles as plsEx)
from rrompy.reduction_methods import (
RationalInterpolantGreedyPivotedGreedyPoleMatch as RIGPG)
from rrompy.parameter.parameter_sampling import (QuadratureSampler as QS,
SparseGridSampler as SGS)
zs, Ls = [10., 50.], [.2, 1.2]
z0, L0, n = np.mean(zs), np.mean(Ls), 50
murange = [[zs[0], Ls[0]], [zs[-1], Ls[-1]]]
np.random.seed(4020)
mu = [zs[0] + np.random.rand() * (zs[-1] - zs[0]),
Ls[0] + np.random.rand() * (Ls[-1] - Ls[0])]
solver = engine(z0, L0, n)
fighandles = []
params = {"POD": True, "nTestPoints": 100, "greedyTol": 1e-4, "S": 3,
"polybasisMarginal": "MONOMIAL_WENDLAND",
"polybasis": "LEGENDRE", 'samplerPivot':QS(zs, "UNIFORM"),
'trainSetGenerator':QS(zs, "UNIFORM"),
'errorEstimatorKind':"LOOK_AHEAD_RES",
'errorEstimatorKindMarginal':"LOOK_AHEAD_RECOVER",
"SMarginal": 3, "paramsMarginal": {"MMarginal": 2,
"radialDirectionalWeightsMarginalAdapt": [1e9, 1e12]},
"greedyTolMarginal": 1e-2, "samplerMarginal":SGS(Ls),
"radialDirectionalWeightsMarginal": [4.], "matchingWeight": 1.}
for shared, tol in product([1., 0.], [1., 3.]):
print("Testing cutoff tolerance {} with shared ratio {}.".format(tol,
shared))
solver.cutOffPolesRMinRel = - 1. - tol
solver.cutOffPolesRMaxRel = 1. + tol
- params['sharedRatio'] = shared
+ params['matchingShared'] = shared
approx = RIGPG([0], solver, mu0 = [z0, L0], approxParameters = params,
verbosity = 5)
approx.setupApprox("ALL")
verb = approx.verbosity
approx.verbosity = 0
tspace = np.linspace(Ls[0], Ls[-1], 100)
for j, t in enumerate(tspace):
plsE = plsEx(t, 0., zs[-1])
pls = approx.getPoles([None, t])
pls[np.abs(np.imag(pls)) > 1e-5] = np.nan
if j == 0:
polesE = np.empty((len(tspace), len(plsE)))
poles = np.empty((len(tspace), len(pls)))
polesE[:] = np.nan
if len(plsE) > polesE.shape[1]:
nanR = np.empty((len(tspace), len(plsE) - polesE.shape[1]))
nanR[:] = np.nan
polesE = np.hstack((polesE, nanR))
polesE[j, : len(plsE)] = np.real(plsE)
poles[j] = np.real(pls)
approx.verbosity = verb
fighandles += [plt.figure(figsize = (17, 5))]
ax1 = fighandles[-1].add_subplot(1, 2, 1)
ax2 = fighandles[-1].add_subplot(1, 2, 2)
ax1.plot(poles, tspace)
ax1.set_ylim(Ls)
ax1.set_xlabel('mu_1')
ax1.set_ylabel('mu_2')
ax1.grid()
ax2.plot(polesE, tspace, 'k-.', linewidth = 1)
ax2.plot(poles, tspace)
for mm in approx.musMarginal:
ax2.plot(zs, [mm[0, 0]] * 2, 'k--', linewidth = 1)
ax2.set_xlim(zs)
ax2.set_ylim(Ls)
ax2.set_xlabel('mu_1')
ax2.set_ylabel('mu_2')
ax2.grid()
plt.show()
print("\n")
diff --git a/examples/9_active_remeshing/active_remeshing.py b/examples/9_active_remeshing/active_remeshing.py
index 6c033c4..f247787 100755
--- a/examples/9_active_remeshing/active_remeshing.py
+++ b/examples/9_active_remeshing/active_remeshing.py
@@ -1,118 +1,118 @@
import numpy as np
from pickle import load
from matplotlib import pyplot as plt
from active_remeshing_engine import ActiveRemeshingEngine
from rrompy.reduction_methods import (
RationalInterpolantGreedyPivotedNoMatch as RIGPNM,
RationalInterpolantGreedyPivotedGreedyPoleMatch as RIGPG)
from rrompy.parameter.parameter_sampling import (QuadratureSampler as QS,
SparseGridSampler as SGS)
zs, ts = [0., 100.], [0., .5]
z0, t0, n = np.mean(zs), np.mean(ts), 150
solver = ActiveRemeshingEngine(z0, t0, n)
mus = [[z0, ts[0]], [z0, ts[1]]]
for mu in mus:
u = solver.solve(mu, return_state = True)[0]
Y = solver.applyC(u, mu)[0][0]
_ = solver.plot(u, what = "REAL", name = "u(z={}, t={})".format(*mu),
is_state = True, figsize = (12, 4))
print("Y(z={}, t={}) = {} (solution norm squared)".format(*mu, Y))
fighandles = []
with open("./active_remeshing_hf_samples.pkl", "rb") as f:
zspace, tspace, Yex = load(f)
# zspace = np.linspace(zs[0], zs[-1], 198)
# tspace = np.linspace(ts[0], ts[-1], 50)
# Yex = np.log10([[solver.solve([z, t]) for t in tspace] for z in zspace])
# (from a ~2.5h simulation on one node of the EPFL Helvetios cluster)
for match in [0, 1]:
params = {"POD": True, "S": 5, "greedyTol": 1e-4, "nTestPoints": 500,
"polybasis": "LEGENDRE", "trainSetGenerator": QS(zs, "UNIFORM"),
"samplerPivot":QS(zs, "CHEBYSHEV"), "samplerMarginal":SGS(ts),
"errorEstimatorKind": "LOOK_AHEAD_OUTPUT"}
if match:
print("\nTesting output-based matching with weight 1.")
params["SMarginal"] = 3
params["maxIterMarginal"] = 25
params["greedyTolMarginal"] = 1e-2
params["matchingWeight"] = 1.
- params["sharedRatio"] = .75
+ params["matchingShared"] = .75
params["polybasisMarginal"] = "PIECEWISE_LINEAR_UNIFORM"
params["errorEstimatorKindMarginal"] = "LOOK_AHEAD_RECOVER"
algo = RIGPG
else:
print("\nTesting matching-free approach.")
params["SMarginal"] = 5
algo = RIGPNM
approx = algo([0], solver, mu0 = [z0, t0], approxParameters = params,
verbosity = 15)
approx.setupApprox("ALL" * match)
verb, verbTM = approx.verbosity, approx.trainedModel.verbosity
approx.verbosity, approx.trainedModel.verbosity = 0, 0
for j, t in enumerate(tspace):
out = approx.getApprox(np.pad(zspace.reshape(-1, 1), [(0, 0), (0, 1)],
"constant", constant_values = t))
pls = approx.getPoles([None, t])
pls[np.abs(np.imag(pls)) > 1e-5] = np.nan
if j == 0:
Ys = np.empty((len(zspace), len(tspace)))
poles = np.empty((len(tspace), len(pls)))
Ys[:, j] = np.log10(out.data)
if len(pls) > poles.shape[1]:
poles = np.pad(poles, [(0, 0), (0, len(pls) - poles.shape[1])],
"constant", constant_values = np.nan)
poles[j, : len(pls)] = np.real(pls)
approx.verbosity, approx.trainedModel.verbosity = verb, verbTM
Ymin, Ymax = min(np.min(Ys), np.min(Yex)), max(np.max(Ys), np.max(Yex))
fighandles += [plt.figure(figsize = (15, 5))]
ax1 = fighandles[-1].add_subplot(1, 2, 1)
ax2 = fighandles[-1].add_subplot(1, 2, 2)
if match:
ax1.plot(poles, tspace)
else:
ax1.plot(poles, tspace, "k.")
ax1.set_ylim(ts)
ax1.set_xlabel("z")
ax1.set_ylabel("t")
ax1.grid()
if match:
ax2.plot(poles, tspace)
else:
ax2.plot(poles, tspace, "k.")
for mm in approx.musMarginal:
ax2.plot(zs, [mm[0, 0]] * 2, "k--", linewidth = 1)
ax2.set_xlim(zs)
ax2.set_ylim(ts)
ax2.set_xlabel("z")
ax2.set_ylabel("t")
ax2.grid()
plt.show()
print("Approximate poles")
fighandles += [plt.figure(figsize = (15, 5))]
ax1 = fighandles[-1].add_subplot(1, 2, 1)
ax2 = fighandles[-1].add_subplot(1, 2, 2)
p = ax1.contourf(np.repeat(zspace.reshape(-1, 1), len(tspace), axis = 1),
np.repeat(tspace.reshape(1, -1), len(zspace), axis = 0),
Ys, vmin = Ymin, vmax = Ymax, cmap = "gray_r",
levels = np.linspace(Ymin, Ymax, 50))
plt.colorbar(p, ax = ax1)
ax1.set_xlabel("z")
ax1.set_ylabel("t")
ax1.grid()
p = ax2.contourf(np.repeat(zspace.reshape(-1, 1), len(tspace), axis = 1),
np.repeat(tspace.reshape(1, -1), len(zspace), axis = 0),
Yex, vmin = Ymin, vmax = Ymax, cmap = "gray_r",
levels = np.linspace(Ymin, Ymax, 50))
ax2.set_xlabel("z")
ax2.set_ylabel("t")
ax2.grid()
plt.colorbar(p, ax = ax2)
plt.show()
print("Approximate and exact output\n")
diff --git a/rrompy/hfengines/base/hfengine_base.py b/rrompy/hfengines/base/hfengine_base.py
index bce8ffa..ae452c1 100644
--- a/rrompy/hfengines/base/hfengine_base.py
+++ b/rrompy/hfengines/base/hfengine_base.py
@@ -1,402 +1,427 @@
# 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, List, DictAny, paramVal,
paramList, sampList)
from rrompy.utilities.numerical import solve as tsolve, dot
from rrompy.utilities.expression import expressionEvaluator
from rrompy.utilities.exception_manager import RROMPyException, 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
self.outputNormMatrix = 1.
def name(self) -> str:
return self.__class__.__name__
def __str__(self) -> str:
return self.name()
def __repr__(self) -> str:
return self.__str__() + " at " + hex(id(self))
def __dir_base__(self):
return [x for x in self.__dir__() if x[:2] != "__"]
def __deepcopy__(self, memo):
return softcopy(self)
@property
def npar(self):
"""Value of npar."""
return self._npar
@npar.setter
def npar(self, npar):
nparOld = self._npar if hasattr(self, "_npar") else -1
if npar != nparOld:
self.parameterMap = pMap(1., npar)
self._npar = npar
@property
def spacedim(self):
return 1
def checkParameter(self, mu:paramVal) -> paramVal:
muP = checkParameter(mu, self.npar)
if self.npar == 0: muP.reset((1, 0), muP.dtype)
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
def buildEnergyNormForm(self):
"""
Build sparse matrix (in CSR format) representative of scalar product.
"""
self.energyNormMatrix = 1.
def buildEnergyNormDualForm(self):
"""
Build sparse matrix (in CSR format) representative of dual scalar
product without duality.
"""
self.energyNormDualMatrix = 1.
def innerProduct(self, u:Np2D, v:Np2D, onlyDiag : bool = False,
dual : bool = False, is_state : bool = True) -> Np2D:
"""Scalar product."""
if is_state or self.isCEye:
if dual:
if not hasattr(self, "energyNormDualMatrix"):
self.buildEnergyNormDualForm()
energyMat = self.energyNormDualMatrix
else:
if not hasattr(self, "energyNormMatrix"):
self.buildEnergyNormForm()
energyMat = self.energyNormMatrix
else:
energyMat = self.outputNormMatrix
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
def setCQuadratic(self, quad : bool = 1):
if isinstance(quad, (str,)):
quad = quad.upper().strip().replace(" ","")
if quad == "SELFADJOINT":
quad = 2
elif quad in ["STANDARD", "GENERAL"]:
quad = 1
else:
raise RROMPyException(("String keyword for output symmetry "
"not recognized."))
if quad:
self._is_C_quadratic = int(quad)
elif hasattr(self, "_is_C_quadratic"):
del self._is_C_quadratic
@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."""
if not (hasattr(self, "_is_C_quadratic") and self._is_C_quadratic):
return dot(self.C(mu), u)
return self._applyCQuadratic(u, mu)
def _applyCQuadratic(self, u:sampList, mu:paramVal, v : sampList = None,
onlyDiag : bool = False):
"""Apply quadratic LHS of linear system."""
if not (hasattr(self, "_is_C_quadratic") and self._is_C_quadratic):
raise RROMPyException(("Cannot call quadratic output routine if "
"output is not quadratic."))
C = self.C(mu)
is_C_list = isinstance(C, list)
if ((is_C_list and np.any([Ci.ndim != 2 for Ci in C]))
or not (is_C_list or C.ndim in [2, 3])):
raise RROMPyException(("C array for quadratic output must have 2 "
"or 3 dimensions or be list of "
"2-dimensional arrays."))
symmetry = v is None # computing u^H * C * u
selfadjoint = self._is_C_quadratic == 2 and symmetry
if isinstance(u, sampleList): u = u.data
if isinstance(v, sampleList): v = v.data
while u.ndim < 2: u = np.expand_dims(u, -1)
if v is None: v = u
while v.ndim < 2: v = np.expand_dims(v, -1)
N, M = u.shape[1], v.shape[1]
if onlyDiag:
if N != M:
raise RROMPyException(("Cannot extract diagonal of "
"rectangular output."))
Ncol = (N,)
elif symmetry:
Ncol = (N * (N + 1) // 2,) if selfadjoint else (N ** 2,)
else:
Ncol = (N, M)
for j in range(N):
Rv = v[:, j] if onlyDiag else v[:, j :] if selfadjoint else v
if is_C_list:
cj = np.array([[dot(dot(Ci, u[:, j]), Rv.conj())] for Ci in C])
else:
cj = dot(dot(C, u[:, j]), Rv.conj())
if C.ndim == 2: cj = np.expand_dims(cj, 0)
if selfadjoint and (onlyDiag or N == 1): cj = np.real(cj)
if j == 0: res = np.empty((len(cj),) + Ncol, dtype = cj[0].dtype)
if not onlyDiag and symmetry: # map 2D idx to hierarchical 1D idx
if self._is_C_quadratic == 2:
res[:, j * (j + 1) // 2] = cj[:, 0]
iD = np.arange(j + 1, N) * np.arange(j + 4, N + 3) // 2 - j
res[:, iD] = 2 * cj[:, 1 :] # exploit symmetry
else:
res[:, j ** 2 : j * (j + 1)] = cj[:, : j]
iD = np.arange(j, N) * np.arange(j + 2, N + 2) - j
res[:, iD] = cj[:, j :]
else:
res[:, j] = cj
return res
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()
mu = self.checkParameterList(mu)
mu_loc, idx, sizes = listScatter(mu, return_sizes = True)
mu_loc = self.checkParameterList(mu_loc)
req, emptyCores = [], np.where(np.logical_not(sizes))[0]
if len(mu_loc) == 0:
uL, uT = recv(source = 0, tag = poolRank())
sol = np.empty((uL, 0), dtype = uT)
else:
applyCglob = (hasattr(self, "_is_C_quadratic")
and self._is_C_quadratic)
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 or applyCglob): u = self.applyC(u, mj)
if j == 0:
sol = np.empty((len(u), len(mu_loc)), dtype = u.dtype)
if masterCore():
for dest in emptyCores:
req += [isend((len(u), u.dtype), dest = dest,
tag = dest)]
sol[:, j] = u
for r in req: r.wait()
sol = matrixGatherv(sol, sizes)
if not return_state and applyCglob: sol = self.applyC(sol, mu)
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()
mu_loc, idx, sizes = listScatter(mu, return_sizes = True)
mu_loc = self.checkParameterList(mu_loc)
req, emptyCores = [], np.where(np.logical_not(sizes))[0]
if len(mu_loc) == 0:
uL, uT = recv(source = 0, tag = poolRank())
res = np.empty((uL, 0), dtype = uT)
else:
applyCglob = (hasattr(self, "_is_C_quadratic")
and self._is_C_quadratic)
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 and not applyCglob: r = self.applyC(r, mj)
if j == 0:
res = np.empty((len(r), len(mu_loc)), dtype = r.dtype)
if masterCore():
for dest in emptyCores:
req += [isend((len(r), r.dtype), dest = dest,
tag = dest)]
res[:, j] = r
for r in req: r.wait()
res = matrixGatherv(res, sizes)
if post_c and applyCglob: res = self.applyC(res, mu)
return sampleList(res)
cutOffPolesRMax,cutOffPolesRMin = np.inf, - np.inf
cutOffPolesRMaxRel, cutOffPolesRMinRel = np.inf, - np.inf
cutOffPolesIMax, cutOffPolesIMin = np.inf, - np.inf
cutOffPolesIMaxRel, cutOffPolesIMinRel = np.inf, - np.inf
cutOffResNormMin = -1
+ cutOffResAngleMin, cutOffResAngleMax = -1, np.pi + 1
+ @property
+ def _ignoreResidues(self):
+ return (self.cutOffResNormMin <= 0. and self.cutOffResAngleMin <= 0.
+ and self.cutOffResAngleMax >= np.pi)
+
def flagBadPolesResidues(self, poles:Np1D, residues : Np1D = None,
- relative : bool = False) -> Np1D:
+ relative : bool = False,
+ projMat : Np2D = None) -> Np1D:
"""
Flag (numerical) poles/residues which are impossible.
Args:
poles: poles to be judged.
- residues: residues to be judged.
+ residues: residues norms to be judged.
relative: whether relative values should be used for poles.
+ projMat: matrix for projection of residues.
"""
poles = np.array(poles).flatten()
flag = np.zeros(len(poles), dtype = bool)
- if residues is None:
- self._ignoreResidues = self.cutOffResNormMin <= 0.
- if relative:
- RMax, RMin = self.cutOffPolesRMaxRel, self.cutOffPolesRMinRel
- IMax, IMin = self.cutOffPolesIMaxRel, self.cutOffPolesIMinRel
- else:
- RMax, RMin = self.cutOffPolesRMax, self.cutOffPolesRMin
- IMax, IMin = self.cutOffPolesIMax, self.cutOffPolesIMin
- if not np.isinf(RMax):
- flag = np.logical_or(flag, np.real(poles) > RMax)
- if not np.isinf(RMin):
- flag = np.logical_or(flag, np.real(poles) < RMin)
- if not np.isinf(IMax):
- flag = np.logical_or(flag, np.imag(poles) > IMax)
- if not np.isinf(IMin):
- flag = np.logical_or(flag, np.imag(poles) < IMin)
+ if relative:
+ RMax, RMin = self.cutOffPolesRMaxRel, self.cutOffPolesRMinRel
+ IMax, IMin = self.cutOffPolesIMaxRel, self.cutOffPolesIMinRel
else:
- residues = np.array(residues).reshape(len(poles), -1)
+ RMax, RMin = self.cutOffPolesRMax, self.cutOffPolesRMin
+ IMax, IMin = self.cutOffPolesIMax, self.cutOffPolesIMin
+ if not np.isinf(RMax):
+ flag = np.logical_or(flag, np.real(poles) > RMax)
+ if not np.isinf(RMin):
+ flag = np.logical_or(flag, np.real(poles) < RMin)
+ if not np.isinf(IMax):
+ flag = np.logical_or(flag, np.imag(poles) > IMax)
+ if not np.isinf(IMin):
+ flag = np.logical_or(flag, np.imag(poles) < IMin)
+ if residues is not None and not self._ignoreResidues:
+ residues = np.array(residues).reshape(-1, len(poles))
+ if projMat is None:
+ resNorm = np.linalg.norm(residues, axis = 0)
+ else:
+ residues = projMat.dot(residues)
+ resNorm = self.norm(residues)
if self.cutOffResNormMin > 0.:
- if residues.shape[1] == self.spacedim:
- resEff = self.norm(residues.T)
+ resNormEff = resNorm / np.max(resNorm)
+ flag = np.logical_or(flag, resNormEff < self.cutOffResNormMin)
+ if self.cutOffResAngleMin > 0. or self.cutOffResAngleMax < np.pi:
+ if projMat is None:
+ angles = np.real(residues.T.conj().dot(residues))
else:
- resEff = np.linalg.norm(residues, axis = 1)
- resEff /= np.max(resEff)
- flag = np.logical_or(flag, resEff < self.cutOffResNormMin)
+ angles = np.real(self.innerProduct(residues, residues))
+ resNormEff = resNorm
+ resNormEff[np.isclose(resNormEff, 0.)] = 1.
+ angles = np.clip((angles / resNormEff).T / resNormEff, -1., 1.)
+ angles = np.arccos(angles)
+ badangles = np.logical_or(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] = 1
+ flag = np.logical_or(flag, idx)
return flag
diff --git a/rrompy/reduction_methods/pivoted/generic_pivoted_approximant.py b/rrompy/reduction_methods/pivoted/generic_pivoted_approximant.py
index afe3b99..41a88bb 100644
--- a/rrompy/reduction_methods/pivoted/generic_pivoted_approximant.py
+++ b/rrompy/reduction_methods/pivoted/generic_pivoted_approximant.py
@@ -1,783 +1,800 @@
# 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 os import mkdir, remove, rmdir
import numpy as np
from collections.abc import Iterable
from copy import deepcopy as copy
from rrompy.reduction_methods.base.generic_approximant import (
GenericApproximant)
+from .trained_model.convert_trained_model_pivoted import (
+ convertTrainedModelPivoted)
from rrompy.utilities.base.data_structures import purgeDict, getNewFilename
from rrompy.utilities.poly_fitting.polynomial import polybases as ppb
from rrompy.utilities.poly_fitting.radial_basis import polybases as rbpb
from rrompy.utilities.poly_fitting.piecewise_linear import sparsekinds as sk
from rrompy.utilities.base.types import Np2D, paramList, List, ListAny
from rrompy.utilities.base import verbosityManager as vbMng
from rrompy.utilities.numerical.degree import reduceDegreeN
from rrompy.utilities.exception_manager import RROMPyException, RROMPyWarning
from rrompy.parameter import checkParameterList
from rrompy.utilities.parallel import poolRank, bcast
__all__ = ['GenericPivotedApproximantNoMatch',
'GenericPivotedApproximantPoleMatch']
class GenericPivotedApproximantBase(GenericApproximant):
def __init__(self, directionPivot:ListAny, *args,
storeAllSamples : bool = False, **kwargs):
self._preInit()
if len(directionPivot) > 1:
raise RROMPyException(("Exactly 1 pivot parameter allowed in pole "
"matching."))
from rrompy.parameter.parameter_sampling import (EmptySampler as ES,
SparseGridSampler as SG)
self._addParametersToList(["radialDirectionalWeightsMarginal"], [1.],
["samplerPivot", "SMarginal",
"samplerMarginal"],
[ES(), 1, SG([[-1.], [1.]])],
toBeExcluded = ["sampler"])
self._directionPivot = directionPivot
self.storeAllSamples = storeAllSamples
if not hasattr(self, "_output_lvl"): self._output_lvl = []
self._output_lvl += [1 / 2]
super().__init__(*args, **kwargs)
self._postInit()
def setupSampling(self): super().setupSampling(False)
def initializeModelData(self, datadict):
if "directionPivot" in datadict.keys():
from .trained_model.trained_model_pivoted_data import (
TrainedModelPivotedData)
data = TrainedModelPivotedData(datadict["mu0"], datadict["mus"],
datadict.pop("projMat"),
datadict["scaleFactor"],
datadict.pop("parameterMap"),
datadict["directionPivot"])
if hasattr(self.HFEngine, "_is_C_quadratic") and not (
hasattr(self, "matchState") and self.matchState):
data._is_C_quadratic = self.HFEngine._is_C_quadratic
return (data, ["mu0", "scaleFactor", "directionPivot", "mus"])
else:
return super().initializeModelData(datadict)
@property
def npar(self):
"""Number of parameters."""
if hasattr(self, "_temporaryPivot"): return self.nparPivot
return super().npar
def checkParameterListPivot(self, mu:paramList,
check_if_single : bool = False) -> paramList:
return checkParameterList(mu, self.nparPivot, check_if_single)
def checkParameterListMarginal(self, mu:paramList,
check_if_single : bool = False) -> paramList:
return checkParameterList(mu, self.nparMarginal, check_if_single)
def mapParameterList(self, *args, **kwargs):
if hasattr(self, "_temporaryPivot"):
return self.mapParameterListPivot(*args, **kwargs)
return super().mapParameterList(*args, **kwargs)
def mapParameterListPivot(self, mu:paramList, direct : str = "F",
idx : List[int] = None):
if idx is None:
idx = self.directionPivot
else:
idx = [self.directionPivot[j] for j in idx]
return super().mapParameterList(mu, direct, idx)
def mapParameterListMarginal(self, mu:paramList, direct : str = "F",
idx : List[int] = None):
if idx is None:
idx = self.directionMarginal
else:
idx = [self.directionMarginal[j] for j in idx]
return super().mapParameterList(mu, direct, idx)
@property
def mu0(self):
"""Value of mu0."""
if hasattr(self, "_temporaryPivot"):
return self.checkParameterListPivot(self._mu0(self.directionPivot))
return self._mu0
@mu0.setter
def mu0(self, mu0):
GenericApproximant.mu0.fset(self, mu0)
@property
def mus(self):
"""Value of mus. Its assignment may reset snapshots."""
return self._mus
@mus.setter
def mus(self, mus):
mus = self.checkParameterList(mus)
musOld = copy(self.mus) if hasattr(self, '_mus') else None
if (musOld is None or len(mus) != len(musOld) or not mus == musOld):
self.resetSamples()
self._mus = mus
@property
def musMarginal(self):
"""Value of musMarginal. Its assignment may reset snapshots."""
return self._musMarginal
@musMarginal.setter
def musMarginal(self, musMarginal):
musMarginal = self.checkParameterListMarginal(musMarginal)
if hasattr(self, '_musMarginal'):
musMOld = copy(self.musMarginal)
else:
musMOld = None
if (musMOld is None or len(musMarginal) != len(musMOld)
or not musMarginal == musMOld):
self.resetSamples()
self._musMarginal = musMarginal
@property
def SMarginal(self):
"""Value of SMarginal."""
return self._SMarginal
@SMarginal.setter
def SMarginal(self, SMarginal):
if SMarginal <= 0:
raise RROMPyException("SMarginal must be positive.")
if hasattr(self, "_SMarginal") and self._SMarginal is not None:
Sold = self.SMarginal
else: Sold = -1
self._SMarginal = SMarginal
self._approxParameters["SMarginal"] = self.SMarginal
if Sold != self.SMarginal: self.resetSamples()
@property
def radialDirectionalWeightsMarginal(self):
"""Value of radialDirectionalWeightsMarginal."""
return self._radialDirectionalWeightsMarginal
@radialDirectionalWeightsMarginal.setter
def radialDirectionalWeightsMarginal(self, radialDirWeightsMarg):
if isinstance(radialDirWeightsMarg, Iterable):
radialDirWeightsMarg = list(radialDirWeightsMarg)
else:
radialDirWeightsMarg = [radialDirWeightsMarg]
self._radialDirectionalWeightsMarginal = radialDirWeightsMarg
self._approxParameters["radialDirectionalWeightsMarginal"] = (
self.radialDirectionalWeightsMarginal)
@property
def directionPivot(self):
"""Value of directionPivot. Its assignment may reset snapshots."""
return self._directionPivot
@directionPivot.setter
def directionPivot(self, directionPivot):
if hasattr(self, '_directionPivot'):
directionPivotOld = copy(self.directionPivot)
else:
directionPivotOld = None
if (directionPivotOld is None
or len(directionPivot) != len(directionPivotOld)
or not directionPivot == directionPivotOld):
self.resetSamples()
self._directionPivot = directionPivot
@property
def directionMarginal(self):
return [x for x in range(self.HFEngine.npar) \
if x not in self.directionPivot]
@property
def nparPivot(self):
return len(self.directionPivot)
@property
def nparMarginal(self):
return self.npar - self.nparPivot
@property
def muBounds(self):
"""Value of muBounds."""
return self.samplerPivot.lims
@property
def muBoundsMarginal(self):
"""Value of muBoundsMarginal."""
return self.samplerMarginal.lims
@property
def sampler(self):
"""Proxy of samplerPivot."""
return self._samplerPivot
@property
def samplerPivot(self):
"""Value of samplerPivot."""
return self._samplerPivot
@samplerPivot.setter
def samplerPivot(self, samplerPivot):
if 'generatePoints' not in dir(samplerPivot):
raise RROMPyException("Pivot sampler type not recognized.")
if hasattr(self, '_samplerPivot') and self._samplerPivot is not None:
samplerOld = self.samplerPivot
self._samplerPivot = samplerPivot
self._approxParameters["samplerPivot"] = self.samplerPivot
if not 'samplerOld' in locals() or samplerOld != self.samplerPivot:
self.resetSamples()
@property
def samplerMarginal(self):
"""Value of samplerMarginal."""
return self._samplerMarginal
@samplerMarginal.setter
def samplerMarginal(self, samplerMarginal):
if 'generatePoints' not in dir(samplerMarginal):
raise RROMPyException("Marginal sampler type not recognized.")
if (hasattr(self, '_samplerMarginal')
and self._samplerMarginal is not None):
samplerOld = self.samplerMarginal
self._samplerMarginal = samplerMarginal
self._approxParameters["samplerMarginal"] = self.samplerMarginal
if not 'samplerOld' in locals() or samplerOld != self.samplerMarginal:
self.resetSamples()
def computeScaleFactor(self):
"""Compute parameter rescaling factor."""
self.scaleFactorPivot = .5 * np.abs((
self.mapParameterListPivot(self.muBounds[0])
- self.mapParameterListPivot(self.muBounds[1]))[0])
self.scaleFactorMarginal = .5 * np.abs((
self.mapParameterListMarginal(self.muBoundsMarginal[0])
- self.mapParameterListMarginal(self.muBoundsMarginal[1]))[0])
self.scaleFactor = np.empty(self.npar)
self.scaleFactor[self.directionPivot] = self.scaleFactorPivot
self.scaleFactor[self.directionMarginal] = self.scaleFactorMarginal
def _setupTrainedModel(self, pMat:Np2D, pMatUpdate : bool = False,
pMatOld : Np2D = None, forceNew : bool = False):
if forceNew or self.trainedModel is None:
self.trainedModel = self.tModelType()
self.trainedModel.verbosity = self.verbosity
self.trainedModel.timestamp = self.timestamp
datadict = {"mu0": self.mu0, "mus": copy(self.mus),
"projMat": pMat, "scaleFactor": self.scaleFactor,
"parameterMap": self.HFEngine.parameterMap,
"directionPivot": self.directionPivot}
self.trainedModel.data = self.initializeModelData(datadict)[0]
else:
self.trainedModel = self.trainedModel
if pMatUpdate:
self.trainedModel.data.projMat = np.hstack(
(self.trainedModel.data.projMat, pMat))
else:
self.trainedModel.data.projMat = copy(pMat)
self.trainedModel.data.mus = copy(self.mus)
self.trainedModel.data.musMarginal = copy(self.musMarginal)
def normApprox(self, mu:paramList) -> float:
_PODOld, self._POD = self.POD, 0
result = super().normApprox(mu)
self._POD = _PODOld
return result
@property
def storedSamplesFilenames(self) -> List[str]:
if not hasattr(self, "_sampleBaseFilename"): return []
return [self._sampleBaseFilename
+ "{}_{}.pkl" .format(idx + 1, self.name())
for idx in range(len(self.musMarginal))]
def purgeStoredSamples(self):
if not hasattr(self, "_sampleBaseFilename"): return
for file in self.storedSamplesFilenames: remove(file)
rmdir(self._sampleBaseFilename[: -8])
def storeSamples(self, idx : int = None):
"""Store samples to file."""
if not hasattr(self, "_sampleBaseFilename"):
filenameBase = None
if poolRank() == 0:
foldername = getNewFilename(self.name(), "samples")
mkdir(foldername)
filenameBase = foldername + "/sample_"
self._sampleBaseFilename = bcast(filenameBase, force = True)
if idx is not None:
super().storeSamples(self._sampleBaseFilename + str(idx + 1),
False)
def loadTrainedModel(self, filename:str):
"""Load trained reduced model from file."""
super().loadTrainedModel(filename)
self._musMarginal = self.trainedModel.data.musMarginal
+ def setTrainedModel(self, model):
+ """Deepcopy approximation from trained model."""
+ super().setTrainedModel(model)
+ self.trainedModel = convertTrainedModelPivoted(self.trainedModel,
+ self.tModelType, self,
+ True)
+ self._finalizeMarginalization()
+ self.trainedModel.data.approxParameters = self.approxParameters
+
class GenericPivotedApproximantNoMatch(GenericPivotedApproximantBase):
"""
ROM pivoted approximant (without pole matching) computation for parametric
problems (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';
- '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;
- 'radialDirectionalWeightsMarginal': radial basis weights for
marginal interpolant; defaults to 1.
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;
- 'radialDirectionalWeightsMarginal': radial basis weights for
marginal interpolant.
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.
radialDirectionalWeightsMarginal: Radial basis weights for marginal
interpolant.
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.
"""
@property
def tModelType(self):
from .trained_model.trained_model_pivoted_rational_nomatch import (
TrainedModelPivotedRationalNoMatch)
return TrainedModelPivotedRationalNoMatch
def _finalizeMarginalization(self):
self.trainedModel.setupMarginalInterp(
[self.radialDirectionalWeightsMarginal])
self.trainedModel.data.approxParameters = copy(self.approxParameters)
def _preliminaryMarginalFinalization(self):
pass
class GenericPivotedApproximantPoleMatch(GenericPivotedApproximantBase):
"""
ROM pivoted approximant (with pole matching) computation for parametric
problems (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;
- - 'sharedRatio': required ratio of marginal points to share
+ - 'matchingShared': required ratio of marginal points to share
resonance; defaults to 1.;
- '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;
- '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_*';
- . 'interpRcondMarginal': tolerance for marginal interpolation;
+ . '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.
- 'radialDirectionalWeightsMarginal': radial basis weights for
marginal interpolant; defaults to 1.
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;
- - 'sharedRatio': required ratio of marginal points to share
+ - 'matchingShared': required ratio of marginal points to share
resonance;
- '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;
- . 'interpRcondMarginal': tolerance for marginal interpolation;
+ . 'interpTolMarginal': tolerance for marginal interpolation;
. 'radialDirectionalWeightsMarginalAdapt': bounds for adaptive
rescaling of marginal radial basis weights.
- 'radialDirectionalWeightsMarginal': radial basis weights for
marginal interpolant.
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.
- sharedRatio: Required ratio of marginal points to share resonance.
+ matchingShared: Required ratio of marginal points to share resonance.
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.
polybasisMarginal: Type of polynomial basis for marginal interpolation.
paramsMarginal: Dictionary of parameters for marginal interpolation.
radialDirectionalWeightsMarginal: Radial basis weights for marginal
interpolant.
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 __init__(self, *args, **kwargs):
self._preInit()
self._addParametersToList(["matchState", "matchingWeight",
- "sharedRatio", "polybasisMarginal",
+ "matchingShared", "polybasisMarginal",
"paramsMarginal"],
[False, 1., 1., "MONOMIAL", {}])
self.parameterMarginalList = ["MMarginal", "nNeighborsMarginal",
"polydegreetypeMarginal",
- "interpRcondMarginal",
+ "interpTolMarginal",
"radialDirectionalWeightsMarginalAdapt"]
super().__init__(*args, **kwargs)
self._postInit()
@property
def tModelType(self):
from .trained_model.trained_model_pivoted_rational_polematch import (
TrainedModelPivotedRationalPoleMatch)
return TrainedModelPivotedRationalPoleMatch
@property
def matchState(self):
"""Value of matchState."""
return self._matchState
@matchState.setter
def matchState(self, matchState):
self._matchState = matchState
self._approxParameters["matchState"] = self.matchState
@property
def matchingWeight(self):
"""Value of matchingWeight."""
return self._matchingWeight
@matchingWeight.setter
def matchingWeight(self, matchingWeight):
self._matchingWeight = matchingWeight
self._approxParameters["matchingWeight"] = self.matchingWeight
@property
- def sharedRatio(self):
- """Value of sharedRatio."""
- return self._sharedRatio
- @sharedRatio.setter
- def sharedRatio(self, sharedRatio):
- if sharedRatio > 1.:
+ def matchingShared(self):
+ """Value of matchingShared."""
+ return self._matchingShared
+ @matchingShared.setter
+ def matchingShared(self, matchingShared):
+ if matchingShared > 1.:
RROMPyWarning("Shared ratio too large. Clipping to 1.")
- sharedRatio = 1.
- elif sharedRatio < 0.:
+ matchingShared = 1.
+ elif matchingShared < 0.:
RROMPyWarning("Shared ratio too small. Clipping to 0.")
- sharedRatio = 0.
- self._sharedRatio = sharedRatio
- self._approxParameters["sharedRatio"] = self.sharedRatio
+ matchingShared = 0.
+ self._matchingShared = matchingShared
+ self._approxParameters["matchingShared"] = self.matchingShared
@property
def polybasisMarginal(self):
"""Value of polybasisMarginal."""
return self._polybasisMarginal
@polybasisMarginal.setter
def polybasisMarginal(self, polybasisMarginal):
try:
polybasisMarginal = polybasisMarginal.upper().strip().replace(" ",
"")
if polybasisMarginal not in ppb + rbpb + ["NEARESTNEIGHBOR"] + sk:
raise RROMPyException(
"Prescribed marginal polybasis not recognized.")
self._polybasisMarginal = polybasisMarginal
except:
RROMPyWarning(("Prescribed marginal polybasis not recognized. "
"Overriding to 'MONOMIAL'."))
self._polybasisMarginal = "MONOMIAL"
self._approxParameters["polybasisMarginal"] = self.polybasisMarginal
@property
def paramsMarginal(self):
"""Value of paramsMarginal."""
return self._paramsMarginal
@paramsMarginal.setter
def paramsMarginal(self, paramsMarginal):
paramsMarginal = purgeDict(paramsMarginal, self.parameterMarginalList,
dictname = self.name() + ".paramsMarginal",
baselevel = 1)
keyList = list(paramsMarginal.keys())
if not hasattr(self, "_paramsMarginal"): self._paramsMarginal = {}
if "MMarginal" in keyList:
MMarg = paramsMarginal["MMarginal"]
elif ("MMarginal" in self.paramsMarginal
and not hasattr(self, "_MMarginal_isauto")):
MMarg = self.paramsMarginal["MMarginal"]
else:
MMarg = "AUTO"
if isinstance(MMarg, str):
MMarg = MMarg.strip().replace(" ","")
if "-" not in MMarg: MMarg = MMarg + "-0"
self._MMarginal_isauto = True
self._MMarginal_shift = int(MMarg.split("-")[-1])
MMarg = 0
if MMarg < 0:
raise RROMPyException("MMarginal must be non-negative.")
self._paramsMarginal["MMarginal"] = MMarg
if "nNeighborsMarginal" in keyList:
self._paramsMarginal["nNeighborsMarginal"] = max(1,
paramsMarginal["nNeighborsMarginal"])
elif "nNeighborsMarginal" not in self.paramsMarginal:
self._paramsMarginal["nNeighborsMarginal"] = 1
if "polydegreetypeMarginal" in keyList:
try:
polydegtypeM = paramsMarginal["polydegreetypeMarginal"]\
.upper().strip().replace(" ","")
if polydegtypeM not in ["TOTAL", "FULL"]:
raise RROMPyException(("Prescribed polydegreetypeMarginal "
"not recognized."))
self._paramsMarginal["polydegreetypeMarginal"] = polydegtypeM
except:
RROMPyWarning(("Prescribed polydegreetypeMarginal not "
"recognized. Overriding to 'TOTAL'."))
self._paramsMarginal["polydegreetypeMarginal"] = "TOTAL"
elif "polydegreetypeMarginal" not in self.paramsMarginal:
self._paramsMarginal["polydegreetypeMarginal"] = "TOTAL"
- if "interpRcondMarginal" in keyList:
- self._paramsMarginal["interpRcondMarginal"] = (
- paramsMarginal["interpRcondMarginal"])
- elif "interpRcondMarginal" not in self.paramsMarginal:
- self._paramsMarginal["interpRcondMarginal"] = -1
+ if "interpTolMarginal" in keyList:
+ self._paramsMarginal["interpTolMarginal"] = (
+ paramsMarginal["interpTolMarginal"])
+ elif "interpTolMarginal" not in self.paramsMarginal:
+ self._paramsMarginal["interpTolMarginal"] = -1
if "radialDirectionalWeightsMarginalAdapt" in keyList:
self._paramsMarginal["radialDirectionalWeightsMarginalAdapt"] = (
paramsMarginal["radialDirectionalWeightsMarginalAdapt"])
elif "radialDirectionalWeightsMarginalAdapt" not in self.paramsMarginal:
self._paramsMarginal["radialDirectionalWeightsMarginalAdapt"] = [
-1., -1.]
self._approxParameters["paramsMarginal"] = self.paramsMarginal
def _setMMarginalAuto(self):
if (self.polybasisMarginal not in ppb + rbpb
or "MMarginal" not in self.paramsMarginal
or "polydegreetypeMarginal" not in self.paramsMarginal):
raise RROMPyException(("Cannot set MMarginal if "
"polybasisMarginal does not allow it."))
self.paramsMarginal["MMarginal"] = max(0, reduceDegreeN(
len(self.musMarginal), len(self.musMarginal),
self.nparMarginal,
self.paramsMarginal["polydegreetypeMarginal"])
- self._MMarginal_shift)
vbMng(self, "MAIN", ("Automatically setting MMarginal to {}.").format(
self.paramsMarginal["MMarginal"]), 25)
def purgeparamsMarginal(self):
self.paramsMarginal = {}
paramsMbadkeys = []
if self.polybasisMarginal in ppb + rbpb + sk:
paramsMbadkeys += ["nNeighborsMarginal"]
if self.polybasisMarginal not in rbpb:
paramsMbadkeys += ["radialDirectionalWeightsMarginalAdapt"]
if self.polybasisMarginal in ["NEARESTNEIGHBOR"] + sk:
paramsMbadkeys += ["MMarginal", "polydegreetypeMarginal",
- "interpRcondMarginal"]
+ "interpTolMarginal"]
if hasattr(self, "_MMarginal_isauto"): del self._MMarginal_isauto
if hasattr(self, "_MMarginal_shift"): del self._MMarginal_shift
for key in paramsMbadkeys:
if key in self._paramsMarginal: del self._paramsMarginal[key]
self._approxParameters["paramsMarginal"] = self.paramsMarginal
def _finalizeMarginalization(self):
vbMng(self, "INIT", "Checking shared ratio.", 10)
- msg = self.trainedModel.checkSharedRatio(self.sharedRatio)
+ msg = self.trainedModel.checkShared(self.matchingShared)
vbMng(self, "DEL", "Done checking." + msg, 10)
if self.polybasisMarginal in rbpb + ["NEARESTNEIGHBOR"]:
self.computeScaleFactor()
rDWMEff = np.array([w * f for w, f in zip(
self.radialDirectionalWeightsMarginal,
self.scaleFactorMarginal)])
if self.polybasisMarginal in ppb + rbpb + sk:
interpPars = [self.polybasisMarginal]
if self.polybasisMarginal in ppb + rbpb:
if self.polybasisMarginal in rbpb: interpPars += [rDWMEff]
interpPars += [self.verbosity >= 5,
self.paramsMarginal["polydegreetypeMarginal"] == "TOTAL"]
if self.polybasisMarginal in ppb:
interpPars += [{}]
else: # if self.polybasisMarginal in rbpb:
interpPars += [{"optimizeScalingBounds":self.paramsMarginal[
"radialDirectionalWeightsMarginalAdapt"]}]
interpPars += [
- {"rcond":self.paramsMarginal["interpRcondMarginal"]}]
+ {"rcond":self.paramsMarginal["interpTolMarginal"]}]
extraPar = hasattr(self, "_MMarginal_isauto")
else: # if self.polybasisMarginal in sk:
idxEff = [x for x in range(self.samplerMarginal.npoints)
if not hasattr(self.trainedModel, "_idxExcl")
or x not in self.trainedModel._idxExcl]
extraPar = self.samplerMarginal.depth[idxEff]
else: # if self.polybasisMarginal == "NEARESTNEIGHBOR":
interpPars = [self.paramsMarginal["nNeighborsMarginal"], rDWMEff]
extraPar = None
self.trainedModel.setupMarginalInterp(self, interpPars, extraPar)
self.trainedModel.data.approxParameters = copy(self.approxParameters)
def _preliminaryMarginalFinalization(self):
vbMng(self, "INIT", "Compressing and matching poles.", 10)
self.trainedModel.initializeFromRational(self.matchingWeight,
self.HFEngine,
self.matchState)
vbMng(self, "DEL", "Done compressing and matching poles.", 10)
def _postApplyC(self):
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."))
applyCglob = (hasattr(self.HFEngine, "_is_C_quadratic")
and self.HFEngine._is_C_quadratic)
vbMng(self, "INIT", "Extracting system output from state.", 35)
if applyCglob:
pMat, dirM = None, self.trainedModel.data.directionMarginal
for muM in self.trainedModel.data.musMarginal:
idx = np.where([np.allclose(mu(dirM)[0], muM[0])
for mu in self.trainedModel.data.mus])[0]
pMati = self.trainedModel.data.projMat[:, idx]
musi = self.trainedModel.data.mus[idx]
pMati = self.HFEngine.applyC(pMati, musi)
if pMat is None:
pMat = np.array(pMati)
else:
pMat = np.append(pMat, pMati, axis = 1)
else:
pMat = None
for j, mu in enumerate(musi):
pMij = self.trainedModel.data.projMat[:, j]
pMij = np.expand_dims(self.HFEngine.applyC(pMij, mu), -1)
if pMati is None:
pMat = np.array(pMij)
else:
pMat = np.append(pMat, pMij, axis = 1)
vbMng(self, "DEL", "Done extracting system output.", 35)
self.trainedModel.data.projMat = pMat
if hasattr(self.HFEngine, "_is_C_quadratic"):
self.trainedModel.data._is_C_quadratic = (
self.HFEngine._is_C_quadratic)
+ def setTrainedModel(self, model):
+ """Deepcopy approximation from trained model."""
+ super().setTrainedModel(model)
+ self._preliminaryMarginalFinalization()
+ self.trainedModel.data.approxParameters = self.approxParameters
+
@abstractmethod
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/greedy/generic_pivoted_greedy_approximant.py b/rrompy/reduction_methods/pivoted/greedy/generic_pivoted_greedy_approximant.py
index 74cdf02..7f00194 100644
--- a/rrompy/reduction_methods/pivoted/greedy/generic_pivoted_greedy_approximant.py
+++ b/rrompy/reduction_methods/pivoted/greedy/generic_pivoted_greedy_approximant.py
@@ -1,654 +1,654 @@
# 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,
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.point_matching import (pointMatching,
buildResiduesForDistance)
from rrompy.utilities.numerical.point_distances import chordalMetricAngleMatrix
from rrompy.utilities.exception_manager import (RROMPyException, RROMPyAssert,
RROMPyWarning)
from rrompy.parameter import emptyParameterList
from rrompy.utilities.parallel import (masterCore, indicesScatter,
arrayGatherv, isend)
__all__ = ['GenericPivotedGreedyApproximantPoleMatch']
class GenericPivotedGreedyApproximantBase(GenericPivotedApproximantBase):
_allowedEstimatorKindsMarginal = ["LOOK_AHEAD", "LOOK_AHEAD_RECOVER",
"NONE"]
def __init__(self, *args, **kwargs):
self._preInit()
self._addParametersToList(["matchingWeightError",
"errorEstimatorKindMarginal",
"greedyTolMarginal", "maxIterMarginal"],
[0., "NONE", 1e-1, 1e2])
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 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()
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()
def _getDistanceApp(self, polesEx:Np1D, resEx:Np2D, muTest:paramVal,
foci:Tuple[float, float], ground:float) -> float:
polesAp = self.trainedModel.interpolateMarginalPoles(muTest)[0]
if self.matchingWeightError != 0:
resAp = self.trainedModel.interpolateMarginalCoeffs(muTest)[0][
: len(polesAp), :]
if (hasattr(self.trainedModel.data, "_is_C_quadratic")
and self.trainedModel.data._is_C_quadratic):
self.trainedModel._setupQuadMapping()
projMapping = self.quad_mapping
projMappingReal = self.data._is_C_quadratic == 2
else:
projMapping, projMappingReal = None, False
resEx = buildResiduesForDistance(resEx,
self.trainedModel.data.projMat,
0, projMapping, projMappingReal)
resAp = buildResiduesForDistance(resAp,
self.trainedModel.data.projMat,
0, projMapping, projMappingReal)
else:
resAp = None
dist = chordalMetricAngleMatrix(polesEx, polesAp,
self.matchingWeightError, resEx, resAp,
self.HFEngine, False)
pmR, pmC = pointMatching(dist)
return np.mean(dist[pmR, pmC])
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
foci = self.samplerPivot.normalFoci()
ground = self.samplerPivot.groundPotential()
for j in idx:
muTest = self.trainedModel._musMExcl[j]
HITest = self.trainedModel._HIsExcl[j]
polesEx = HITest.poles
idxGood = np.logical_not(np.logical_or(np.isinf(polesEx),
np.isnan(polesEx)))
polesEx = polesEx[idxGood]
if self.matchingWeightError != 0:
resEx = HITest.coeffs[np.where(idxGood)[0]]
else:
resEx = None
if len(polesEx) == 0:
err += [0.]
continue
err += [self._getDistanceApp(polesEx, resEx, muTest,
foci, ground)]
self.verbosity += 25
self.trainedModel.verbosity += 25
return arrayGatherv(np.array(err), sizes)
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_max : 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 = (1. + self.greedyTolMarginal) * np.ones(nErr)
else:#if self.errorEstimatorKindMarginal[: 10] == "LOOK_AHEAD":
err = self.getErrorEstimatorMarginalLookAhead()
vbMng(self.trainedModel, "DEL", "Done evaluating error estimator.", 10)
if not return_max: return err
idxMaxEst = np.where(err > self.greedyTolMarginal)[0]
maxErr = err[idxMaxEst]
if self.errorEstimatorKindMarginal == "NONE": maxErr = None
return err, idxMaxEst, maxErr
def plotEstimatorMarginal(self, est:Np1D, idxMax:List[int],
estMax: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(idxMax) > 0 and estMax is not None:
maxrej = musre[idxMax, 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.))[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)
ax.semilogy(self.musMarginal.re(jpar),
(self.greedyTolMarginal,) * len(self.musMarginal),
'*m')
if len(idxMax) > 0 and estMax is not None:
ax.semilogy(maxrej, estMax, 'xr')
ax.set_xlim(*list(self.samplerMarginal.lims.re(jpar)))
ax.grid()
plt.tight_layout()
plt.show()
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, maxErrorEst = self.errorEstimatorMarginal(True)
if plotEst == "ALL":
self.plotEstimatorMarginal(errorEstTest, muidx, maxErrorEst)
return (errorEstTest, muidx, maxErrorEst,
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)
idx, sizes = indicesScatter(len(mus), return_sizes = True)
emptyCores = np.where(np.logical_not(sizes))[0]
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
padLeft = self.trainedModel.data.projMat.shape[1]
suppNew = np.append(0, np.cumsum(nsamples))
self._setupTrainedModel(pMat, padLeft > 0)
self.trainedModel.data.Qs += Qs
self.trainedModel.data.Ps += Ps
self.trainedModel.data.Psupp += list(padLeft + suppNew[: -1])
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]:
pMati = self.samplingEngine.projectionMatrix
musi = self.samplingEngine.mus
if not hasattr(self, "matchState") or 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)
if (hasattr(self.HFEngine, "_is_C_quadratic")
and self.HFEngine._is_C_quadratic):
pMati = self.HFEngine.applyC(pMati, musi)
else:
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), pMat.dtype, mus.dtype),
dest = dest, tag = dest)]
else:
mus = np.vstack((mus, musi.data))
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
RROMPyAssert(self._mode, message = "Cannot start greedy algorithm.")
vbMng(self, "INIT", "Starting computation of snapshots.", 3)
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, maxErrorEst, mu = \
self.greedyNextSampleMarginal(muidx, plotEst)
if maxErrorEst is None:
max2ErrorEst = 1. + self.greedyTolMarginal
else:
if len(maxErrorEst) > 0:
max2ErrorEst = np.max(maxErrorEst)
else:
max2ErrorEst = np.max(errorEstTest)
vbMng(self, "MAIN", ("Uniform testing error estimate "
"{:.4e}.").format(max2ErrorEst), 3)
if plotEst == "LAST":
self.plotEstimatorMarginal(errorEstTest, muidx, maxErrorEst)
vbMng(self, "DEL", ("Done computing snapshots (final snapshot count: "
"{}).").format(len(self.mus)), 3)
if (self.errorEstimatorKindMarginal == "LOOK_AHEAD_RECOVER"
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._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)
return 0
def checkComputedApproxPivoted(self) -> bool:
return (super().checkComputedApprox()
and len(self.musMarginal) == len(self.trainedModel.data.musMarginal))
class GenericPivotedGreedyApproximantPoleMatch(
GenericPivotedGreedyApproximantBase,
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;
- - 'sharedRatio': required ratio of marginal points to share
+ - 'matchingShared': required ratio of marginal points to share
resonance; defaults to 1.;
- 'matchingWeightError': weight for pole matching optimization in
error estimation; defaults to 0;
- '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_RECOVER',
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_*';
- . 'interpRcondMarginal': tolerance for marginal interpolation;
+ . '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.
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;
- - 'sharedRatio': required ratio of marginal points to share
+ - 'matchingShared': required ratio of marginal points to share
resonance;
- 'matchingWeightError': weight for pole matching optimization in
error estimation;
- '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;
- . 'interpRcondMarginal': tolerance for marginal interpolation;
+ . '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.
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.
- sharedRatio: Required ratio of marginal points to share resonance.
+ matchingShared: Required ratio of marginal points to share resonance.
matchingWeightError: Weight for pole matching optimization in error
estimation.
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.
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 _updateTrainedModelMarginalSamples(self, idx : ListAny = []):
self.trainedModel.updateEffectiveSamples(idx, self.matchingWeight,
self.HFEngine, False)
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
self._finalizeMarginalization()
if self.matchState: self._postApplyC()
return setupOK
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 d78a5d9..0dc1bc2 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,357 +1,355 @@
#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,
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 (
RationalInterpolantGreedyPivotedPoleMatch)
from rrompy.utilities.base.types import Np1D, Tuple, paramVal, 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__ = ['RationalInterpolantGreedyPivotedGreedyPoleMatch']
class RationalInterpolantGreedyPivotedGreedyBase(
GenericPivotedGreedyApproximantBase):
@property
def sampleBatchSize(self):
"""Value of sampleBatchSize."""
return 1
@property
def sampleBatchIdx(self):
"""Value of sampleBatchIdx."""
return self.S
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 _setSampleBatch(self, maxS:int):
return self.S
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.trainSetGenerator.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))
for k in range(self.S - 1):
muk = np.empty_like(self.mus[0])
muk[self.directionPivot] = musPivot[k]
muk[self.directionMarginal] = self.muMargLoc
self.mus[k] = muk
for k in range(len(muTestBasePivot)):
muk = np.empty_like(self.muTest[0])
muk[self.directionPivot] = muTestBasePivot[k]
muk[self.directionMarginal] = self.muMargLoc
self.muTest[k] = muk
muk = np.empty_like(self.mus[0])
muk[self.directionPivot] = musPivot[-1]
muk[self.directionMarginal] = self.muMargLoc
self.muTest[-1] = muk
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.M, self.N = ("AUTO",) * 2
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, pT, mT = recv(source = 0, tag = poolRank())
pMat = np.empty((pL, 0), dtype = pT)
musA = np.empty((0, self.mu0.shape[1]), dtype = mT)
else:
for i in idx:
self.muMargLoc = mus[i]
vbMng(self, "MAIN", "Building marginal model no. {} at "
"{}.".format(i + 1, self.muMargLoc), 25)
self.samplingEngine.resetHistory()
self.trainedModel = None
self.verbosity -= 5
self.samplingEngine.verbosity -= 10
RationalInterpolantGreedy.setupApprox(self, self._plotEstPivot)
self.verbosity += 5
self.samplingEngine.verbosity += 10
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)]
self._S = S0
del self.muMargLoc
for r in req: r.wait()
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 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;
- - 'sharedRatio': required ratio of marginal points to share
+ - 'matchingShared': required ratio of marginal points to share
resonance; defaults to 1.;
- 'matchingWeightError': weight for pole matching optimization in
error estimation; defaults to 0;
- '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' and 'LOOK_AHEAD_RECOVER';
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_*';
- . 'interpRcondMarginal': tolerance for marginal interpolation;
+ . '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;
- 'trainSetGenerator': training sample points generator; defaults
to sampler;
- '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;
- 'functionalSolve': strategy for minimization of denominator
functional; allowed values include 'NORM', 'DOMINANT', 'NODAL',
'BARYCENTRIC_NORM', and 'BARYCENTRIC[_AVERAGE]' (check pdf in
main folder for meaning); defaults to 'NORM';
- - 'interpRcond': tolerance for pivot interpolation; defaults to
- None;
- - 'robustTol': tolerance for robust rational denominator
- management; defaults to 0.
+ - 'interpTol': tolerance for pivot interpolation; defaults to None;
+ - '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;
- - 'sharedRatio': required ratio of marginal points to share
+ - 'matchingShared': required ratio of marginal points to share
resonance;
- 'matchingWeightError': weight for pole matching optimization in
error estimation;
- '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;
- . 'interpRcondMarginal': tolerance for marginal interpolation;
+ . '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;
- 'trainSetGenerator': 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;
- 'functionalSolve': strategy for minimization of denominator
functional;
- - 'interpRcond': tolerance for pivot interpolation;
- - 'robustTol': tolerance for robust rational denominator
- management.
+ - 'interpTol': tolerance for pivot interpolation;
+ - '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.
- sharedRatio: Required ratio of marginal points to share resonance.
+ matchingShared: Required ratio of marginal points to share resonance.
matchingWeightError: Weight for pole matching optimization in error
estimation.
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.
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.
trainSetGenerator: 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.
functionalSolve: Strategy for minimization of denominator functional.
- interpRcond: Tolerance for pivot interpolation.
- robustTol: Tolerance for robust rational denominator management.
+ interpTol: Tolerance for pivot interpolation.
+ 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/greedy/rational_interpolant_pivoted_greedy.py b/rrompy/reduction_methods/pivoted/greedy/rational_interpolant_pivoted_greedy.py
index e33f8b0..a6887a6 100644
--- a/rrompy/reduction_methods/pivoted/greedy/rational_interpolant_pivoted_greedy.py
+++ b/rrompy/reduction_methods/pivoted/greedy/rational_interpolant_pivoted_greedy.py
@@ -1,287 +1,285 @@
# 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
from numpy import empty, empty_like
from .generic_pivoted_greedy_approximant import (
GenericPivotedGreedyApproximantBase,
GenericPivotedGreedyApproximantPoleMatch)
from rrompy.reduction_methods.standard import RationalInterpolant
from rrompy.reduction_methods.pivoted import (
RationalInterpolantPivotedPoleMatch)
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__ = ['RationalInterpolantPivotedGreedyPoleMatch']
class RationalInterpolantPivotedGreedyBase(
GenericPivotedGreedyApproximantBase):
def computeSnapshots(self):
"""Compute snapshots of solution map."""
RROMPyAssert(self._mode,
message = "Cannot start snapshot computation.")
vbMng(self, "INIT", "Starting computation of snapshots.", 5)
self.samplingEngine.scaleFactor = self.scaleFactorDer
if not hasattr(self, "musPivot") or len(self.musPivot) != self.S:
self.musPivot = self.samplerPivot.generatePoints(self.S)
while len(self.musPivot) > self.S: self.musPivot.pop()
musLoc = emptyParameterList()
musLoc.reset((self.S, self.HFEngine.npar))
self.samplingEngine.resetHistory()
for k in range(self.S):
muk = empty_like(musLoc[0])
muk[self.directionPivot] = self.musPivot[k]
muk[self.directionMarginal] = self.muMargLoc
musLoc[k] = muk
self.samplingEngine.iterSample(musLoc)
vbMng(self, "DEL", "Done computing snapshots.", 5)
self._m_selfmus = copy(musLoc)
self._mus = self.musPivot
self._m_HFEparameterMap = copy(self.HFEngine.parameterMap)
self.HFEngine.parameterMap = {
"F": [self.HFEngine.parameterMap["F"][self.directionPivot[0]]],
"B": [self.HFEngine.parameterMap["B"][self.directionPivot[0]]]}
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)
idx, sizes, emptyCores = self._preSetupApproxPivoted(mus)
pMat, Ps, Qs, req, musA = None, [], [], [], None
if len(idx) == 0:
vbMng(self, "MAIN", "Idling.", 45)
if self.storeAllSamples: self.storeSamples()
pL, pT, mT = recv(source = 0, tag = poolRank())
pMat = empty((pL, 0), dtype = pT)
musA = empty((0, self.mu0.shape[1]), dtype = mT)
else:
for i in idx:
self.muMargLoc = mus[i]
vbMng(self, "MAIN", "Building marginal model no. {} at "
"{}.".format(i + 1, self.muMargLoc), 25)
self.samplingEngine.resetHistory()
self.trainedModel = None
self.verbosity -= 5
self.samplingEngine.verbosity -= 5
RationalInterpolant.setupApprox(self)
self.verbosity += 5
self.samplingEngine.verbosity += 5
self._mus = self._m_selfmus
self.HFEngine.parameterMap = self._m_HFEparameterMap
del self._m_selfmus, self._m_HFEparameterMap
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)]
del self.muMargLoc
for r in req: r.wait()
self._postSetupApproxPivoted(musA, pMat, Ps, Qs, sizes)
vbMng(self, "DEL", "Done setting up pivoted approximant.", 10)
return 0
class RationalInterpolantPivotedGreedyPoleMatch(
RationalInterpolantPivotedGreedyBase,
GenericPivotedGreedyApproximantPoleMatch,
RationalInterpolantPivotedPoleMatch):
"""
ROM 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;
- - 'sharedRatio': required ratio of marginal points to share
+ - 'matchingShared': required ratio of marginal points to share
resonance; defaults to 1.;
- 'matchingWeightError': weight for pole matching optimization in
error estimation; defaults to 0;
- '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' and 'LOOK_AHEAD_RECOVER';
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_*';
- . 'interpRcondMarginal': tolerance for marginal interpolation;
+ . '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;
- 'greedyTolMarginal': uniform error tolerance for marginal greedy
algorithm; defaults to 1e-1;
- 'maxIterMarginal': maximum number of marginal greedy steps;
defaults to 1e2;
- 'radialDirectionalWeights': radial basis weights for pivot
numerator; defaults to 1;
- 'radialDirectionalWeightsAdapt': bounds for adaptive rescaling of
radial basis weights; defaults to [-1, -1];
- 'radialDirectionalWeightsMarginal': radial basis weights for
marginal interpolant; defaults to 1;
- 'functionalSolve': strategy for minimization of denominator
functional; allowed values include 'NORM', 'DOMINANT', 'NODAL',
'BARYCENTRIC_NORM', and 'BARYCENTRIC[_AVERAGE]' (check pdf in
main folder for meaning); defaults to 'NORM';
- - 'interpRcond': tolerance for pivot interpolation; defaults to
- None;
- - 'robustTol': tolerance for robust rational denominator
- management; defaults to 0.
+ - 'interpTol': tolerance for pivot interpolation; defaults to None;
+ - '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;
- - 'sharedRatio': required ratio of marginal points to share
+ - 'matchingShared': required ratio of marginal points to share
resonance;
- 'matchingWeightError': weight for pole matching optimization in
error estimation;
- '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;
- . 'interpRcondMarginal': tolerance for marginal interpolation;
+ . '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;
- 'greedyTolMarginal': uniform error tolerance for marginal greedy
algorithm;
- 'maxIterMarginal': maximum number of marginal greedy steps;
- 'radialDirectionalWeights': radial basis weights for pivot
numerator;
- 'radialDirectionalWeightsAdapt': bounds for adaptive rescaling of
radial basis weights;
- 'radialDirectionalWeightsMarginal': radial basis weights for
marginal interpolant;
- 'functionalSolve': strategy for minimization of denominator
functional;
- - 'interpRcond': tolerance for pivot interpolation;
- - 'robustTol': tolerance for robust rational denominator
- management.
+ - 'interpTol': tolerance for pivot interpolation;
+ - '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.
- sharedRatio: Required ratio of marginal points to share resonance.
+ matchingShared: Required ratio of marginal points to share resonance.
matchingWeightError: Weight for pole matching optimization in error
estimation.
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.
polybasis: Type of polynomial basis for pivot interpolation.
polybasisMarginal: Type of polynomial basis for marginal interpolation.
paramsMarginal: Dictionary of parameters for marginal interpolation.
M: Degree of rational interpolant numerator.
N: Degree of rational interpolant denominator.
greedyTolMarginal: Uniform error tolerance for marginal greedy
algorithm.
maxIterMarginal: Maximum number of marginal greedy steps.
radialDirectionalWeights: Radial basis weights for pivot numerator.
radialDirectionalWeightsAdapt: Bounds for adaptive rescaling of radial
basis weights.
radialDirectionalWeightsMarginal: Radial basis weights for marginal
interpolant.
functionalSolve: Strategy for minimization of denominator functional.
- interpRcond: Tolerance for pivot interpolation.
- robustTol: Tolerance for robust rational denominator management.
+ interpTol: Tolerance for pivot interpolation.
+ 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 bd6b08f..64611cc 100644
--- a/rrompy/reduction_methods/pivoted/rational_interpolant_greedy_pivoted.py
+++ b/rrompy/reduction_methods/pivoted/rational_interpolant_greedy_pivoted.py
@@ -1,572 +1,569 @@
# 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,
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
from rrompy.utilities.base import verbosityManager as vbMng
from rrompy.utilities.poly_fitting.polynomial import polyvander as pv
from rrompy.utilities.exception_manager import RROMPyException, RROMPyAssert
from rrompy.parameter import emptyParameterList, parameterList
from rrompy.utilities.parallel import poolRank, indicesScatter, isend, recv
__all__ = ['RationalInterpolantGreedyPivotedNoMatch',
'RationalInterpolantGreedyPivotedPoleMatch']
class RationalInterpolantGreedyPivotedBase(GenericPivotedApproximantBase,
RationalInterpolantGreedy):
def __init__(self, *args, **kwargs):
self._preInit()
super().__init__(*args, **kwargs)
- self._ignoreResidues = self.nparPivot > 1
+ 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):
degEff = [0] * self.npar
degEff[self.directionPivot[0]] = deg
return pv(mus, degEff, *args)
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 = self._musUniqueCN.dtype)
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._S = self._setSampleBatch(self.S)
self.resetSamples()
self.samplingEngine.scaleFactor = self.scaleFactorDer
musPivot = self.trainSetGenerator.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])
self._mus = emptyParameterList()
self.mus.reset((self.S, self.npar + len(self.musMargLoc)))
muTestBase = emptyParameterList()
muTestBase.reset((len(muTestPivot), self.npar + len(self.musMargLoc)))
for k in range(self.S):
muk = np.empty_like(self.mus[0])
muk[self.directionPivot] = musPivot[k]
muk[self.directionMarginal] = self.musMargLoc
self.mus[k] = muk
for k in range(len(muTestPivot)):
muk = np.empty_like(muTestBase[0])
muk[self.directionPivot] = muTestPivot[k]
muk[self.directionMarginal] = self.musMargLoc
muTestBase[k] = muk
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 = parameterList(muTestBase)
self.muTest.append(muLast)
self.M, self.N = ("AUTO",) * 2
def setupApproxLocal(self) -> int:
"""Compute rational interpolant."""
self._marginalizeMiscellanea(True)
setupOK = super().setupApproxLocal()
self._marginalizeMiscellanea(False)
return setupOK
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._musMarginal = self.samplerMarginal.generatePoints(self.SMarginal)
while len(self.musMarginal) > self.SMarginal: self.musMarginal.pop()
S0 = copy(self.S)
idx, sizes = indicesScatter(len(self.musMarginal), return_sizes = True)
pMat, Ps, Qs, mus = None, [], [], None
req, emptyCores = [], np.where(np.logical_not(sizes))[0]
if len(idx) == 0:
vbMng(self, "MAIN", "Idling.", 25)
if self.storeAllSamples: self.storeSamples()
pL, pT, mT = recv(source = 0, tag = poolRank())
pMat = np.empty((pL, 0), dtype = pT)
mus = np.empty((0, self.mu0.shape[1]), dtype = mT)
else:
_scaleFactorOldPivot = copy(self.scaleFactor)
self.scaleFactor = self.scaleFactorPivot
self._temporaryPivot = 1
for i in idx:
self.musMargLoc = self.musMarginal[i]
vbMng(self, "MAIN",
"Building marginal model no. {} at {}.".format(i + 1,
self.musMargLoc), 5)
self.samplingEngine.resetHistory()
self.trainedModel = None
self.verbosity -= 5
self.samplingEngine.verbosity -= 10
RationalInterpolantGreedy.setupApprox(self, *args, **kwargs)
self.verbosity += 5
self.samplingEngine.verbosity += 10
if self.storeAllSamples: self.storeSamples(i)
musi = self.samplingEngine.mus
pMati = self.samplingEngine.projectionMatrix
if not hasattr(self, "matchState") or 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)
if (hasattr(self.HFEngine, "_is_C_quadratic")
and self.HFEngine._is_C_quadratic):
pMati = self.HFEngine.applyC(pMati, musi)
else:
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 i == 0:
for dest in emptyCores:
req += [isend((len(pMat), pMat.dtype, mus.dtype),
dest = dest, tag = dest)]
else:
mus = np.vstack((mus, musi.data))
pMat = np.hstack((pMat, pMati))
Ps += [copy(self.trainedModel.data.P)]
Qs += [copy(self.trainedModel.data.Q)]
self._S = S0
del self._temporaryPivot, self.musMargLoc
self.scaleFactor = _scaleFactorOldPivot
for r in req: r.wait()
pMat, Ps, Qs, mus, nsamples = gatherPivotedApproximant(pMat, Ps, Qs,
mus, sizes,
self.polybasis)
self._mus = self.checkParameterList(mus)
Psupp = np.append(0, np.cumsum(nsamples))
self._setupTrainedModel(pMat, forceNew = True)
self.trainedModel.data.Qs, self.trainedModel.data.Ps = Qs, Ps
self.trainedModel.data.Psupp = list(Psupp[: -1])
self._preliminaryMarginalFinalization()
self._finalizeMarginalization()
vbMng(self, "DEL", "Done setting up approximant.", 5)
return 0
class RationalInterpolantGreedyPivotedNoMatch(
RationalInterpolantGreedyPivotedBase,
GenericPivotedApproximantNoMatch):
"""
ROM pivoted rational interpolant (without 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';
- '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';
- '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;
- 'trainSetGenerator': training sample points generator; defaults
to sampler;
- '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', 'NODAL',
'BARYCENTRIC_NORM', and 'BARYCENTRIC[_AVERAGE]' (check pdf in
main folder for meaning); defaults to 'NORM';
- - 'interpRcond': tolerance for pivot interpolation; defaults to
+ - 'interpTol': tolerance for pivot interpolation; defaults to
None;
- - 'robustTol': tolerance for robust rational denominator
- management; 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;
- '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;
- 'trainSetGenerator': training sample points generator;
- 'errorEstimatorKind': kind of error estimator;
- 'radialDirectionalWeightsMarginal': radial basis weights for
marginal interpolant;
- 'functionalSolve': strategy for minimization of denominator
functional;
- - 'interpRcond': tolerance for pivot interpolation;
- - 'robustTol': tolerance for robust rational denominator
- management.
+ - 'interpTol': tolerance for pivot interpolation;
+ - '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.
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.
trainSetGenerator: training sample points generator.
errorEstimatorKind: kind of error estimator.
radialDirectionalWeightsMarginal: Radial basis weights for marginal
interpolant.
functionalSolve: Strategy for minimization of denominator functional.
- interpRcond: Tolerance for pivot interpolation.
- robustTol: Tolerance for robust rational denominator management.
+ interpTol: Tolerance for pivot interpolation.
+ 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 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;
- - 'sharedRatio': required ratio of marginal points to share
+ - 'matchingShared': required ratio of marginal points to share
resonance; defaults to 1.;
- '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_*';
- . 'interpRcondMarginal': tolerance for marginal interpolation;
+ . '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;
- 'trainSetGenerator': training sample points generator; defaults
to sampler;
- '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', 'NODAL',
'BARYCENTRIC_NORM', and 'BARYCENTRIC[_AVERAGE]' (check pdf in
main folder for meaning); defaults to 'NORM';
- - 'interpRcond': tolerance for pivot interpolation; defaults to
- None;
- - 'robustTol': tolerance for robust rational denominator
- management; defaults to 0.
+ - 'interpTol': tolerance for pivot interpolation; defaults to None;
+ - '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;
- - 'sharedRatio': required ratio of marginal points to share
+ - 'matchingShared': required ratio of marginal points to share
resonance;
- '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;
- . 'interpRcondMarginal': tolerance for marginal interpolation;
+ . '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;
- 'trainSetGenerator': training sample points generator;
- 'errorEstimatorKind': kind of error estimator;
- 'radialDirectionalWeightsMarginal': radial basis weights for
marginal interpolant;
- 'functionalSolve': strategy for minimization of denominator
functional;
- - 'interpRcond': tolerance for pivot interpolation;
- - 'robustTol': tolerance for robust rational denominator
- management.
+ - 'interpTol': tolerance for pivot interpolation;
+ - '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.
- sharedRatio: Required ratio of marginal points to share resonance.
+ matchingShared: Required ratio of marginal points to share resonance.
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.
trainSetGenerator: training sample points generator.
errorEstimatorKind: kind of error estimator.
radialDirectionalWeightsMarginal: Radial basis weights for marginal
interpolant.
functionalSolve: Strategy for minimization of denominator functional.
- interpRcond: Tolerance for pivot interpolation.
- robustTol: Tolerance for robust rational denominator management.
+ interpTol: Tolerance for pivot interpolation.
+ 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 95086fd..3142c71 100644
--- a/rrompy/reduction_methods/pivoted/rational_interpolant_pivoted.py
+++ b/rrompy/reduction_methods/pivoted/rational_interpolant_pivoted.py
@@ -1,487 +1,484 @@
# 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,
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.numerical.hash_derivative import nextDerivativeIndices
from rrompy.utilities.exception_manager import (RROMPyException, RROMPyAssert,
RROMPyWarning)
from rrompy.parameter import emptyParameterList
from rrompy.utilities.parallel import poolRank, indicesScatter, isend, recv
__all__ = ['RationalInterpolantPivotedNoMatch',
'RationalInterpolantPivotedPoleMatch']
class RationalInterpolantPivotedBase(GenericPivotedApproximantBase,
RationalInterpolant):
def __init__(self, *args, **kwargs):
self._preInit()
self._addParametersToList(toBeExcluded = ["polydegreetype"])
super().__init__(*args, **kwargs)
- self._ignoreResidues = self.nparPivot > 1
+ if self.nparPivot > 1: self.HFEngine._ignoreResidues = 1
self._postInit()
@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
@property
def polydegreetype(self):
"""Value of polydegreetype."""
return "TOTAL"
@polydegreetype.setter
def polydegreetype(self, polydegreetype):
RROMPyWarning(("polydegreetype is used just to simplify inheritance, "
"and its value cannot be changed from 'TOTAL'."))
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 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.musPivot = self.samplerPivot.generatePoints(self.S)
while len(self.musPivot) > self.S: self.musPivot.pop()
self._musMarginal = self.samplerMarginal.generatePoints(self.SMarginal)
while len(self.musMarginal) > self.SMarginal: self.musMarginal.pop()
self._mus = emptyParameterList()
self.mus.reset((self.S * self.SMarginal, self.HFEngine.npar))
for j, muMarg in enumerate(self.musMarginal):
for k in range(j * self.S, (j + 1) * self.S):
muk = np.empty_like(self.mus[0])
muk[self.directionPivot] = self.musPivot[k - j * self.S]
muk[self.directionMarginal] = muMarg
self.mus[k] = muk
N0 = copy(self.N)
self._setupTrainedModel(np.zeros((0, 0)), forceNew = True)
idx, sizes = indicesScatter(len(self.musMarginal), return_sizes = True)
pMat, Ps, Qs = None, [], []
req, emptyCores = [], np.where(np.logical_not(sizes))[0]
if len(idx) == 0:
vbMng(self, "MAIN", "Idling.", 30)
if self.storeAllSamples: self.storeSamples()
pL, pT = recv(source = 0, tag = poolRank())
pMat = np.empty((pL, 0), dtype = pT)
else:
_scaleFactorOldPivot = copy(self.scaleFactor)
self.scaleFactor = self.scaleFactorPivot
self._temporaryPivot = 1
for i in idx:
musi = self.mus[self.S * i : self.S * (i + 1)]
vbMng(self, "MAIN",
"Building marginal model no. {} at {}.".format(i + 1,
self.musMarginal[i]), 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 -= 10
self._setupRational(self._setupDenominator())
self.verbosity += 5
self.samplingEngine.verbosity += 10
if self.storeAllSamples: self.storeSamples(i)
pMati = self.samplingEngine.projectionMatrix
if not hasattr(self, "matchState") or 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)
if (hasattr(self.HFEngine, "_is_C_quadratic")
and self.HFEngine._is_C_quadratic):
pMati = self.HFEngine.applyC(pMati, musi)
else:
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:
pMat = copy(pMati)
if i == 0:
for dest in emptyCores:
req += [isend((len(pMat), pMat.dtype), dest = dest,
tag = dest)]
else:
pMat = np.hstack((pMat, pMati))
Ps += [copy(self.trainedModel.data.P)]
Qs += [copy(self.trainedModel.data.Q)]
del self.trainedModel.data.Q, self.trainedModel.data.P
self.N = N0
del self._temporaryPivot
self.scaleFactor = _scaleFactorOldPivot
for r in req: r.wait()
pMat, Ps, Qs, _, _ = gatherPivotedApproximant(pMat, Ps, Qs,
self.mus.data, sizes,
self.polybasis, False)
self._setupTrainedModel(pMat)
self.trainedModel.data.Qs, self.trainedModel.data.Ps = Qs, Ps
Psupp = np.arange(0, len(self.musMarginal) * self.S, self.S)
self.trainedModel.data.Psupp = list(Psupp)
self._preliminaryMarginalFinalization()
self._finalizeMarginalization()
vbMng(self, "DEL", "Done setting up approximant.", 5)
return 0
class RationalInterpolantPivotedNoMatch(RationalInterpolantPivotedBase,
GenericPivotedApproximantNoMatch):
"""
ROM pivoted rational interpolant (without 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';
- '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';
- '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;
- 'radialDirectionalWeights': radial basis weights for pivot
numerator; defaults to 1;
- 'radialDirectionalWeightsAdapt': bounds for adaptive rescaling of
radial basis weights; defaults to [-1, -1];
- 'radialDirectionalWeightsMarginal': radial basis weights for
marginal interpolant; defaults to 1;
- 'functionalSolve': strategy for minimization of denominator
functional; allowed values include 'NORM', 'DOMINANT', 'NODAL',
'BARYCENTRIC_NORM', and 'BARYCENTRIC[_AVERAGE]' (check pdf in
main folder for meaning); defaults to 'NORM';
- - 'interpRcond': tolerance for pivot interpolation; defaults to
+ - 'interpTol': tolerance for pivot interpolation; defaults to
None;
- - 'robustTol': tolerance for robust rational denominator
- management; 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;
- 'polybasis': type of polynomial basis for pivot
interpolation;
- 'M': degree of rational interpolant numerator;
- 'N': degree of rational interpolant denominator;
- 'radialDirectionalWeights': radial basis weights for pivot
numerator;
- 'radialDirectionalWeightsAdapt': bounds for adaptive rescaling of
radial basis weights;
- 'radialDirectionalWeightsMarginal': radial basis weights for
marginal interpolant;
- 'functionalSolve': strategy for minimization of denominator
functional;
- - 'interpRcond': tolerance for pivot interpolation;
- - 'robustTol': tolerance for robust rational denominator
- management.
+ - 'interpTol': tolerance for pivot interpolation;
+ - '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.
polybasis: Type of polynomial basis for pivot interpolation.
M: Numerator degree of approximant.
N: Denominator degree of approximant.
radialDirectionalWeights: Radial basis weights for pivot numerator.
radialDirectionalWeightsAdapt: Bounds for adaptive rescaling of radial
basis weights.
radialDirectionalWeightsMarginal: Radial basis weights for marginal
interpolant.
functionalSolve: Strategy for minimization of denominator functional.
- interpRcond: Tolerance for pivot interpolation.
- robustTol: Tolerance for robust rational denominator management.
+ interpTol: Tolerance for pivot interpolation.
+ 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 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;
- - 'sharedRatio': required ratio of marginal points to share
+ - 'matchingShared': required ratio of marginal points to share
resonance; defaults to 1.;
- '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_*';
- . 'interpRcondMarginal': tolerance for marginal interpolation;
+ . '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;
- 'radialDirectionalWeights': radial basis weights for pivot
numerator; defaults to 1;
- 'radialDirectionalWeightsAdapt': bounds for adaptive rescaling of
radial basis weights; defaults to [-1, -1];
- 'radialDirectionalWeightsMarginal': radial basis weights for
marginal interpolant; defaults to 1;
- 'functionalSolve': strategy for minimization of denominator
functional; allowed values include 'NORM', 'DOMINANT', 'NODAL',
'BARYCENTRIC_NORM', and 'BARYCENTRIC[_AVERAGE]' (check pdf in
main folder for meaning); defaults to 'NORM';
- - 'interpRcond': tolerance for pivot interpolation; defaults to
- None;
- - 'robustTol': tolerance for robust rational denominator
- management; defaults to 0.
+ - 'interpTol': tolerance for pivot interpolation; defaults to None;
+ - '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;
- - 'sharedRatio': required ratio of marginal points to share
+ - 'matchingShared': required ratio of marginal points to share
resonance;
- '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;
- . 'interpRcondMarginal': tolerance for marginal interpolation;
+ . '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;
- 'radialDirectionalWeights': radial basis weights for pivot
numerator;
- 'radialDirectionalWeightsAdapt': bounds for adaptive rescaling of
radial basis weights;
- 'radialDirectionalWeightsMarginal': radial basis weights for
marginal interpolant;
- 'functionalSolve': strategy for minimization of denominator
functional;
- - 'interpRcond': tolerance for pivot interpolation;
- - 'robustTol': tolerance for robust rational denominator
- management.
+ - 'interpTol': tolerance for pivot interpolation;
+ - '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.
- sharedRatio: Required ratio of marginal points to share resonance.
+ matchingShared: Required ratio of marginal points to share resonance.
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.
M: Numerator degree of approximant.
N: Denominator degree of approximant.
radialDirectionalWeights: Radial basis weights for pivot numerator.
radialDirectionalWeightsAdapt: Bounds for adaptive rescaling of radial
basis weights.
radialDirectionalWeightsMarginal: Radial basis weights for marginal
interpolant.
functionalSolve: Strategy for minimization of denominator functional.
- interpRcond: Tolerance for pivot interpolation.
- robustTol: Tolerance for robust rational denominator management.
+ interpTol: Tolerance for pivot interpolation.
+ 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/trained_model/convert_trained_model_pivoted.py b/rrompy/reduction_methods/pivoted/trained_model/convert_trained_model_pivoted.py
new file mode 100644
index 0000000..5a4a5eb
--- /dev/null
+++ b/rrompy/reduction_methods/pivoted/trained_model/convert_trained_model_pivoted.py
@@ -0,0 +1,68 @@
+# 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 .trained_model_pivoted_rational_nomatch import (
+ TrainedModelPivotedRationalNoMatch)
+from .trained_model_pivoted_rational_polematch import (
+ TrainedModelPivotedRationalPoleMatch)
+from rrompy.utilities.base import verbosityManager as vbMng
+from rrompy.utilities.exception_manager import RROMPyException, RROMPyWarning
+
+__all__ = ['convertTrainedModelPivoted']
+
+def convertTrainedModelPivoted(model, outType, verbObj = None,
+ muteWarnings : bool = False):
+ if isinstance(model, outType): return model
+ if ((isinstance(model, TrainedModelPivotedRationalNoMatch)
+ and outType == TrainedModelPivotedRationalPoleMatch)
+ or (isinstance(model, TrainedModelPivotedRationalPoleMatch)
+ and outType == TrainedModelPivotedRationalNoMatch)):
+ return convertTrainedModelPivotedMatchUnmatch(model, outType, verbObj,
+ muteWarnings)
+ raise RROMPyException(("Model type or conversion output type not "
+ "recognized."))
+
+def convertTrainedModelPivotedMatchUnmatch(model, outType, verbObj = None,
+ muteWarnings : bool = False):
+ if verbObj is not None:
+ sf, st = ["NoMatch", "PoleMatch"]
+ if outType == TrainedModelPivotedRationalPoleMatch:
+ msgw = "match poles, set up marginalInterp,"
+ else: #if outType == TrainedModelPivotedRationalNoMatch:
+ st, sf = sf, st
+ msgw = "set up marginalInterp"
+ vbMng(verbObj, "INIT",
+ "Starting model conversion from {} to {} model.".format(sf, st),
+ 10)
+ excludeDataKey = ["marginalInterp", "approxParameters"]
+ if outType == TrainedModelPivotedRationalPoleMatch:
+ modelC = TrainedModelPivotedRationalPoleMatch()
+ msgw = "match poles, set up marginalInterp,"
+ else: #if outType == TrainedModelPivotedRationalNoMatch:
+ modelC = TrainedModelPivotedRationalNoMatch()
+ msgw = "set up marginalInterp"
+ excludeDataKey += ["HIs", "suppEffPts", "suppEffIdx", "coeffsEff",
+ "polesEff"]
+ for key in model.__dict__.keys(): setattr(modelC, key, model.__dict__[key])
+ for key in excludeDataKey: delattr(modelC.data, key)
+ if verbObj is not None:
+ vbMng(verbObj, "DEL", "Finished model conversion.", 10)
+ if not muteWarnings:
+ RROMPyWarning(("Model conversion result not yet fuctional: must stil "
+ "{} and assign approxParameters.").format(msgw))
+ return modelC
diff --git a/rrompy/reduction_methods/pivoted/trained_model/trained_model_pivoted_rational_polematch.py b/rrompy/reduction_methods/pivoted/trained_model/trained_model_pivoted_rational_polematch.py
index 21ffc99..4a4f77a 100644
--- a/rrompy/reduction_methods/pivoted/trained_model/trained_model_pivoted_rational_polematch.py
+++ b/rrompy/reduction_methods/pivoted/trained_model/trained_model_pivoted_rational_polematch.py
@@ -1,503 +1,503 @@
# 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 warnings
import numpy as np
from scipy.special import factorial as fact
from scipy.sparse import csr_matrix, hstack, SparseEfficiencyWarning
from collections.abc import Iterable
from copy import deepcopy as copy
from itertools import combinations
from rrompy.reduction_methods.standard.trained_model.trained_model_rational \
import TrainedModelRational
from .trained_model_pivoted_rational_nomatch import (
TrainedModelPivotedRationalNoMatch)
from rrompy.utilities.base.types import (Np1D, Np2D, List, ListAny, paramVal,
paramList, sampList, HFEng)
from rrompy.utilities.base import verbosityManager as vbMng, freepar as fp
from rrompy.utilities.numerical import dot
from rrompy.utilities.numerical.point_matching import rationalFunctionMatching
from rrompy.utilities.numerical.degree import reduceDegreeN
from rrompy.utilities.poly_fitting.polynomial import (polybases as ppb,
PolynomialInterpolator as PI)
from rrompy.utilities.poly_fitting.radial_basis import (polybases as rbpb,
RadialBasisInterpolator as RBI)
from rrompy.utilities.poly_fitting.heaviside import (rational2heaviside,
heavisideUniformShape,
HeavisideInterpolator as HI)
from rrompy.utilities.poly_fitting.nearest_neighbor import (
NearestNeighborInterpolator as NNI)
from rrompy.utilities.poly_fitting.piecewise_linear import (sparsekinds,
PiecewiseLinearInterpolator as PLI)
from rrompy.utilities.exception_manager import RROMPyException, RROMPyAssert
from rrompy.sampling import sampleList, emptySampleList
__all__ = ['TrainedModelPivotedRationalPoleMatch']
class TrainedModelPivotedRationalPoleMatch(TrainedModelPivotedRationalNoMatch):
"""
ROM approximant evaluation for pivoted approximants based on interpolation
of rational approximants (with pole matching).
Attributes:
Data: dictionary with all that can be pickled.
"""
def compress(self, collapse : bool = False, tol : float = 0.,
returnRMat : bool = False, **compressMatrixkwargs):
Psupp = copy(self.data.Psupp)
RMat = super().compress(collapse, tol, True, **compressMatrixkwargs)
if RMat is None: return
for obj, suppj in zip(self.data.HIs, Psupp):
obj.postmultiplyTensorize(RMat.T[suppj : suppj + obj.shape[0]])
if hasattr(self, "_HIsExcl"):
for obj, suppj in zip(self._HIsExcl, Psupp):
obj.postmultiplyTensorize(RMat.T[suppj : suppj + obj.shape[0]])
if not hasattr(self, "_PsExcl"):
self._PsuppExcl = [0] * len(self._PsuppExcl)
if returnRMat: return RMat
def centerNormalizeMarginal(self, mu : paramList = [],
mu0 : paramVal = None) -> paramList:
"""
Compute normalized parameter to be plugged into approximant.
Args:
mu: Parameter(s) 1.
mu0: Parameter(s) 2. If None, set to self.data.mu0Marginal.
Returns:
Normalized parameter.
"""
mu = self.checkParameterListMarginal(mu)
if mu0 is None:
mu0 = self.checkParameterListMarginal(
self.data.mu0(0, self.data.directionMarginal))
return (self.mapParameterList(mu, idx = self.data.directionMarginal)
- self.mapParameterList(mu0, idx = self.data.directionMarginal)
) / [self.data.scaleFactor[x]
for x in self.data.directionMarginal]
def setupMarginalInterp(self, approx, interpPars:ListAny, extraPar = None):
vbMng(self, "INIT", "Starting computation of marginal interpolator.",
12)
musMCN = self.centerNormalizeMarginal(self.data.musMarginal)
nM, pbM = len(musMCN), approx.polybasisMarginal
if pbM in ppb + rbpb:
if extraPar: approx._setMMarginalAuto()
_MMarginalEff = approx.paramsMarginal["MMarginal"]
if pbM in ppb:
p = PI()
elif pbM in rbpb:
p = RBI()
else: # if pbM in sparsekinds + ["NEARESTNEIGHBOR"]:
if pbM == "NEARESTNEIGHBOR":
p = NNI()
else: # if pbM in sparsekinds:
pllims = [[-1.] * self.data.nparMarginal,
[1.] * self.data.nparMarginal]
p = PLI()
for ipts, pts in enumerate(self.data.suppEffPts):
if len(pts) == 0:
raise RROMPyException("Empty list of support points.")
musMCNEff, valsEff = musMCN[pts], np.eye(len(pts))
if pbM in ppb + rbpb:
if extraPar:
if ipts > 0:
verb = approx.verbosity
approx.verbosity = 0
_musM = approx.musMarginal
approx.musMarginal = musMCNEff
approx._setMMarginalAuto()
approx.musMarginal = _musM
approx.verbosity = verb
else:
approx.paramsMarginal["MMarginal"] = reduceDegreeN(
_MMarginalEff, len(musMCNEff), self.data.nparMarginal,
approx.paramsMarginal["polydegreetypeMarginal"])
MMEff = approx.paramsMarginal["MMarginal"]
while MMEff >= 0:
wellCond, msg = p.setupByInterpolation(musMCNEff, valsEff,
MMEff, *interpPars)
vbMng(self, "MAIN", msg, 30)
if wellCond: break
vbMng(self, "MAIN",
("Polyfit is poorly conditioned. Reducing "
"MMarginal by 1."), 35)
MMEff -= 1
if MMEff < 0:
raise RROMPyException(("Instability in computation of "
"interpolant. Aborting."))
if (pbM in rbpb and len(interpPars) > 4
and "optimizeScalingBounds" in interpPars[4].keys()):
interpPars[4]["optimizeScalingBounds"] = [-1., -1.]
elif pbM == "NEARESTNEIGHBOR":
if ipts > 0: interpPars[0] = 1
p.setupByInterpolation(musMCNEff, valsEff, *interpPars)
elif ipts == 0: # and pbM in sparsekinds:
p.setupByInterpolation(musMCNEff, valsEff, pllims,
extraPar[pts], *interpPars)
if ipts == 0:
self.data.marginalInterp = copy(p)
self.data.coeffsEff, self.data.polesEff = [], []
for hi, sup in zip(self.data.HIs, self.data.Psupp):
cEff = hi.coeffs
if (self.data._collapsed
or self.supportEnd == cEff.shape[1]):
cEff = copy(cEff)
else:
supC = self.supportEnd - sup - cEff.shape[1]
cEff = hstack((csr_matrix((len(cEff), sup)),
csr_matrix(cEff),
csr_matrix((len(cEff), supC))), "csr")
self.data.coeffsEff += [cEff]
self.data.polesEff += [copy(hi.poles)]
else:
ptsBad = [i for i in range(nM) if i not in pts]
idxBad = np.where(self.data.suppEffIdx == ipts)[0]
warnings.simplefilter('ignore', SparseEfficiencyWarning)
if pbM in sparsekinds:
for ij, j in enumerate(ptsBad):
nearest = pts[np.argmin(np.sum(np.abs(musMCNEff.data
- np.tile(musMCN[j], [len(pts), 1])
), axis = 1).flatten())]
self.data.coeffsEff[j][idxBad] = copy(
self.data.coeffsEff[nearest][idxBad])
self.data.polesEff[j][idxBad] = copy(
self.data.polesEff[nearest][idxBad])
else:
if (self.data._collapsed
or self.supportEnd == cEff.shape[1]):
cfBase = np.zeros((len(idxBad), cEff.shape[1]),
dtype = cEff.dtype)
else:
cfBase = csr_matrix((len(idxBad), self.supportEnd),
dtype = cEff.dtype)
valMuMBad = p(musMCN[ptsBad])
for ijb, jb in enumerate(ptsBad):
self.data.coeffsEff[jb][idxBad] = copy(cfBase)
self.data.polesEff[jb][idxBad] = 0.
for ij, j in enumerate(pts):
val = valMuMBad[ij][ijb]
if not np.isclose(val, 0.):
self.data.coeffsEff[jb][idxBad] += (val
* self.data.coeffsEff[j][idxBad])
self.data.polesEff[jb][idxBad] += (val
* self.data.polesEff[j][idxBad])
warnings.filters.pop(0)
if pbM in ppb + rbpb:
approx.paramsMarginal["MMarginal"] = _MMarginalEff
vbMng(self, "DEL", "Done computing marginal interpolator.", 12)
def updateEffectiveSamples(self, exclude:List[int], *args, **kwargs):
if hasattr(self, "_idxExcl"):
for j, excl in enumerate(self._idxExcl):
self.data.HIs.insert(excl, self._HIsExcl[j])
super().updateEffectiveSamples(exclude)
self._HIsExcl = []
for excl in self._idxExcl[::-1]:
self._HIsExcl = [self.data.HIs.pop(excl)] + self._HIsExcl
poles = [hi.poles for hi in self.data.HIs]
coeffs = [hi.coeffs for hi in self.data.HIs]
self.initializeFromLists(poles, coeffs, self.data.Psupp,
self.data.HIs[0].polybasis, *args, **kwargs)
def initializeFromRational(self, matchingWeight:float, HFEngine:HFEng,
is_state:bool):
"""Initialize Heaviside representation."""
RROMPyAssert(self.data.nparPivot, 1, "Number of pivot parameters")
poles, coeffs = [], []
for Q, P in zip(self.data.Qs, self.data.Ps):
cfs, pls, basis = rational2heaviside(P, Q)
poles += [pls]
coeffs += [cfs]
self.initializeFromLists(poles, coeffs, self.data.Psupp, basis,
matchingWeight, HFEngine, is_state)
def initializeFromLists(self, poles:ListAny, coeffs:ListAny, supps:ListAny,
basis:str, matchingWeight:float, HFEngine:HFEng,
is_state:bool):
"""Initialize Heaviside representation."""
Ns = [len(pls) for pls in poles]
poles, coeffs = heavisideUniformShape(poles, coeffs)
root = Ns.index(len(poles[0]))
if hasattr(self.data, "_is_C_quadratic") and self.data._is_C_quadratic:
csizemax = np.max([len(c) for c in coeffs])
#csizemax = np.max([c.shape[1] for c in coeffs])
if self.data._is_C_quadratic == 1:
csizemax = csizemax ** 2
else: # if self.data._is_C_quadratic == 2:
csizemax = csizemax * (csizemax + 1) // 2
TrainedModelRational._setupQuadMapping(self, csizemax)
projMapping = self.quad_mapping
projMappingReal = self.data._is_C_quadratic == 2
else:
projMapping, projMappingReal = None, False
poles, coeffs = rationalFunctionMatching(poles, coeffs,
self.data.musMarginal.data,
matchingWeight, supps,
self.data.projMat, HFEngine,
is_state, root, projMapping,
projMappingReal)
self.data.HIs = []
for pls, cfs in zip(poles, coeffs):
hsi = HI()
hsi.poles = pls
if len(cfs) == len(pls):
cfs = np.pad(cfs, ((0, 1), (0, 0)), "constant")
hsi.coeffs = cfs
hsi.npar = 1
hsi.polybasis = basis
self.data.HIs += [hsi]
self.data.suppEffPts = [np.arange(len(self.data.HIs))]
self.data.suppEffIdx = np.zeros(len(poles[0]), dtype = int)
- def checkSharedRatio(self, shared:float) -> str:
+ def checkShared(self, shared:float) -> str:
N = len(self.data.HIs[0].poles)
M = len(self.data.HIs)
goodLocPoles = np.array([np.logical_not(np.isinf(hi.poles)
) for hi in self.data.HIs])
self.data.suppEffPts = [np.arange(len(self.data.HIs))]
self.data.suppEffIdx = np.zeros(N, dtype = int)
if np.all(np.all(goodLocPoles)): return " No poles erased."
goodGlobPoles = np.sum(goodLocPoles, axis = 0)
goodEnoughPoles = goodGlobPoles >= max(1., 1. * shared * M)
keepPole = np.where(goodEnoughPoles)[0]
halfPole = np.where(np.logical_and(goodEnoughPoles,
goodGlobPoles < M))[0]
removePole = np.where(np.logical_not(goodEnoughPoles))[0]
if len(removePole) > 0:
keepCoeff = np.append(keepPole,
np.arange(N, len(self.data.HIs[0].coeffs)))
for hi in self.data.HIs:
for j in removePole:
if not np.isinf(hi.poles[j]):
hi.coeffs[N, :] -= hi.coeffs[j, :] / hi.poles[j]
hi.poles = hi.poles[keepPole]
hi.coeffs = hi.coeffs[keepCoeff, :]
for idxR in halfPole:
pts = np.where(goodLocPoles[:, idxR])[0]
idxEff = len(self.data.suppEffPts)
for idEff, prevPts in enumerate(self.data.suppEffPts):
if len(prevPts) == len(pts):
if np.allclose(prevPts, pts):
idxEff = idEff
break
if idxEff == len(self.data.suppEffPts):
self.data.suppEffPts += [pts]
self.data.suppEffIdx[idxR] = idxEff
self.data.suppEffIdx = self.data.suppEffIdx[keepPole]
return (" Hard-erased {} pole".format(len(removePole))
+ "s" * (len(removePole) != 1)
+ " and soft-erased {} pole".format(len(halfPole))
+ "s" * (len(halfPole) != 1) + ".")
def getApproxReduced(self, mu : paramList = []) -> sampList:
"""
Evaluate reduced representation of approximant at arbitrary parameter.
Args:
mu: Target parameter.
"""
RROMPyAssert(self.data.nparPivot, 1, "Number of pivot parameters")
mu = self.checkParameterList(mu)
if (not hasattr(self, "lastSolvedApproxReduced")
or self.lastSolvedApproxReduced != mu):
vbMng(self, "INIT",
"Evaluating approximant at mu = {}.".format(mu), 12)
muP = self.centerNormalizePivot(mu(self.data.directionPivot))
muM = mu(self.data.directionMarginal)
his = self.interpolateMarginalInterpolator(muM)
for i, (mP, hi) in enumerate(zip(muP, his)):
uAppR = hi(mP)[:, 0]
if i == 0:
uApproxR = np.empty((len(uAppR), len(mu)),
dtype = uAppR.dtype)
uApproxR[:, i] = uAppR
self.uApproxReduced = sampleList(uApproxR)
vbMng(self, "DEL", "Done evaluating approximant.", 12)
self.lastSolvedApproxReduced = mu
return self.uApproxReduced
def interpolateMarginalInterpolator(self, mu : paramList = []) -> ListAny:
"""Obtain interpolated approximant interpolator."""
mu = self.checkParameterListMarginal(mu)
vbMng(self, "INIT",
"Interpolating marginal models at mu = {}.".format(mu), 95)
his = []
muC = self.centerNormalizeMarginal(mu)
mIvals = self.data.marginalInterp(muC)
verb, self.verbosity = self.verbosity, 0
poless = self.interpolateMarginalPoles(mu, mIvals)
coeffss = self.interpolateMarginalCoeffs(mu, mIvals)
self.verbosity = verb
for j in range(len(mu)):
his += [HI()]
his[-1].poles = poless[j]
his[-1].coeffs = coeffss[j]
his[-1].npar = 1
his[-1].polybasis = self.data.HIs[0].polybasis
vbMng(self, "DEL", "Done interpolating marginal models.", 95)
return his
def interpolateMarginalPoles(self, mu : paramList = [],
mIvals : Np2D = None) -> ListAny:
"""Obtain interpolated approximant poles."""
mu = self.checkParameterListMarginal(mu)
vbMng(self, "INIT",
"Interpolating marginal poles at mu = {}.".format(mu), 95)
intMPoles = np.zeros((len(mu),) + self.data.polesEff[0].shape,
dtype = self.data.polesEff[0].dtype)
if mIvals is None:
muC = self.centerNormalizeMarginal(mu)
mIvals = self.data.marginalInterp(muC)
for pEff, mI in zip(self.data.polesEff, mIvals):
intMPoles += np.expand_dims(mI, - 1) * pEff
vbMng(self, "DEL", "Done interpolating marginal poles.", 95)
return intMPoles
def interpolateMarginalCoeffs(self, mu : paramList = [],
mIvals : Np2D = None) -> ListAny:
"""Obtain interpolated approximant coefficients."""
mu = self.checkParameterListMarginal(mu)
vbMng(self, "INIT",
"Interpolating marginal coefficients at mu = {}.".format(mu), 95)
intMCoeffs = np.zeros((len(mu),) + self.data.coeffsEff[0].shape,
dtype = self.data.coeffsEff[0].dtype)
if mIvals is None:
muC = self.centerNormalizeMarginal(mu)
mIvals = self.data.marginalInterp(muC)
for cEff, mI in zip(self.data.coeffsEff, mIvals):
for j, m in enumerate(mI): intMCoeffs[j] += m * cEff
vbMng(self, "DEL", "Done interpolating marginal coefficients.", 95)
return intMCoeffs
def getPVal(self, mu : paramList = []) -> sampList:
"""
Evaluate rational numerator at arbitrary parameter.
Args:
mu: Target parameter.
"""
RROMPyAssert(self.data.nparPivot, 1, "Number of pivot parameters")
mu = self.checkParameterList(mu)
p = emptySampleList()
muP = self.centerNormalizePivot(mu(self.data.directionPivot))
muM = mu(self.data.directionMarginal)
his = self.interpolateMarginalInterpolator(muM)
for i, (mP, hi) in enumerate(zip(muP, his)):
Pval = hi(mP) * np.prod(mP[0] - hi.poles)
if i == 0: p.reset((len(Pval), len(mu)), dtype = Pval.dtype)
p[i] = Pval
return p
def getQVal(self, mu:Np1D, der : List[int] = None,
scl : Np1D = None) -> Np1D:
"""
Evaluate rational denominator at arbitrary parameter.
Args:
mu: Target parameter.
der(optional): Derivatives to take before evaluation.
"""
RROMPyAssert(self.data.nparPivot, 1, "Number of pivot parameters")
mu = self.checkParameterList(mu)
muP = self.centerNormalizePivot(mu(self.data.directionPivot))
muM = mu(self.data.directionMarginal)
if der is None:
derP, derM = 0, [0]
else:
derP = der[self.data.directionPivot[0]]
derM = [der[x] for x in self.data.directionMarginal]
if np.any(np.array(derM) != 0):
raise RROMPyException(("Derivatives of Q with respect to marginal "
"parameters not allowed."))
sclP = 1 if scl is None else scl[self.data.directionPivot[0]]
derVal = np.zeros(len(mu), dtype = np.complex)
pls = self.interpolateMarginalPoles(muM)
for i, (mP, pl) in enumerate(zip(muP, pls)):
N = len(pl)
if derP == N: derVal[i] = 1.
elif derP >= 0 and derP < N:
plDist = mP[0] - pl
for terms in combinations(np.arange(N), N - derP):
derVal[i] += np.prod(plDist[list(terms)])
return sclP ** derP * fact(derP) * derVal
def getPoles(self, *args, **kwargs) -> Np1D:
"""
Obtain approximant poles.
Returns:
Numpy complex vector of poles.
"""
RROMPyAssert(self.data.nparPivot, 1, "Number of pivot parameters")
if len(args) + len(kwargs) > 1:
raise RROMPyException(("Wrong number of parameters passed. "
"Only 1 available."))
elif len(args) + len(kwargs) == 1:
if len(args) == 1:
mVals = args[0]
else:
mVals = kwargs["marginalVals"]
if not isinstance(mVals, Iterable): mVals = [mVals]
mVals = list(mVals)
else:
mVals = [fp]
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."))
if rDim != self.data.directionPivot[0]:
raise RROMPyException(("'freepar' entry in marginalVals must "
"coincide with pivot direction."))
mVals[rDim] = self.data.mu0(rDim)[0]
mMarg = [mVals[j] for j in range(len(mVals)) if j != rDim]
roots = (self.data.scaleFactor[rDim]
* self.interpolateMarginalPoles(mMarg)[0])
return self.mapParameterList(self.mapParameterList(self.data.mu0(rDim),
idx = [rDim])(0, 0)
+ roots, "B", [rDim])(0)
def getResidues(self, *args, **kwargs) -> Np2D:
"""
Obtain approximant residues.
Returns:
Numpy matrix with residues as columns.
"""
pls = self.getPoles(*args, **kwargs)
if len(args) == 1:
mVals = args[0]
elif len(args) == 0:
mVals = [None]
else:
mVals = kwargs["marginalVals"]
if not isinstance(mVals, Iterable): mVals = [mVals]
mVals = list(mVals)
rDim = mVals.index(fp)
mMarg = [mVals[j] for j in range(len(mVals)) if j != rDim]
res = self.interpolateMarginalCoeffs(mMarg)[0][: len(pls), :].T
if hasattr(self.data, "_is_C_quadratic") and self.data._is_C_quadratic:
self._setupQuadMapping()
res = res[self.quad_mapping[0]] * res[self.quad_mapping[1]].conj()
if not self.data._collapsed: res = dot(self.data.projMat, res).T
if (hasattr(self.data, "_is_C_quadratic")
and self.data._is_C_quadratic == 2):
res = np.real(res)
return pls, res
diff --git a/rrompy/reduction_methods/standard/greedy/rational_interpolant_greedy.py b/rrompy/reduction_methods/standard/greedy/rational_interpolant_greedy.py
index d8fa57d..794e59a 100644
--- a/rrompy/reduction_methods/standard/greedy/rational_interpolant_greedy.py
+++ b/rrompy/reduction_methods/standard/greedy/rational_interpolant_greedy.py
@@ -1,531 +1,528 @@
# 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 (polybases, polyfitname,
PolynomialInterpolator as PI,
polyvander)
from rrompy.utilities.numerical import dot
from rrompy.utilities.numerical.degree import totalDegreeN
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,
getVerbosityDepth, setVerbosityDepth)
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;
- '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;
- 'trainSetGenerator': training sample points generator; defaults
to sampler;
- 'polybasis': type of basis for interpolation; defaults to
'MONOMIAL';
- 'errorEstimatorKind': kind of error estimator; available values
include 'AFFINE', 'DISCREPANCY', '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', 'NODAL',
'BARYCENTRIC_NORM', and 'BARYCENTRIC[_AVERAGE]' (check pdf in
main folder for meaning); defaults to 'NORM';
- - 'interpRcond': tolerance for interpolation; defaults to None;
- - 'robustTol': tolerance for robust rational denominator
- management; defaults to 0.
+ - 'interpTol': tolerance for interpolation; defaults to None;
+ - '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;
- 'greedyTol': uniform error tolerance for greedy algorithm;
- 'collinearityTol': collinearity tolerance for greedy algorithm;
- 'maxIter': maximum number of greedy steps;
- 'nTestPoints': number of test points;
- 'trainSetGenerator': training sample points generator;
- 'errorEstimatorKind': kind of error estimator;
- 'functionalSolve': strategy for minimization of denominator
functional;
- - 'interpRcond': tolerance for interpolation;
- - 'robustTol': tolerance for robust rational denominator
- management.
+ - 'interpTol': tolerance for interpolation;
+ - '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.
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.
trainSetGenerator: training sample points generator.
- robustTol: tolerance for robust rational denominator management.
errorEstimatorKind: kind of error estimator.
functionalSolve: Strategy for minimization of denominator functional.
- interpRcond: tolerance for interpolation.
- robustTol: tolerance for robust rational denominator management.
+ interpTol: tolerance for interpolation.
+ 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 = ["M", "N", "polydegreetype",
"radialDirectionalWeights"])
super().__init__(*args, **kwargs)
self._postInit()
@property
def E(self):
"""Value of E."""
self._E = self.sampleBatchIdx - 1
return self._E
@E.setter
def E(self, E):
RROMPyWarning(("E is used just to simplify inheritance, and its value "
"cannot be changed from that of sampleBatchIdx - 1."))
def _setMAuto(self):
self.M = self.E
def _setNAuto(self):
self.N = self.E
@property
def polydegreetype(self):
"""Value of polydegreetype."""
return "TOTAL"
@polydegreetype.setter
def polydegreetype(self, polydegreetype):
RROMPyWarning(("polydegreetype is used just to simplify inheritance, "
"and its value cannot be changed from 'TOTAL'."))
@property
def polybasis(self):
"""Value of polybasis."""
return self._polybasis
@polybasis.setter
def polybasis(self, polybasis):
try:
polybasis = polybasis.upper().strip().replace(" ","")
if polybasis not in polybases:
raise RROMPyException("Sample type not recognized.")
self._polybasis = polybasis
except:
RROMPyWarning(("Prescribed polybasis not recognized. Overriding "
"to 'MONOMIAL'."))
self._polybasis = "MONOMIAL"
self._approxParameters["polybasis"] = self.polybasis
@property
def errorEstimatorKind(self):
"""Value of errorEstimatorKind."""
return self._errorEstimatorKind
@errorEstimatorKind.setter
def errorEstimatorKind(self, errorEstimatorKind):
errorEstimatorKind = errorEstimatorKind.upper()
if errorEstimatorKind not in self._allowedEstimatorKinds:
RROMPyWarning(("Error estimator kind not recognized. Overriding "
"to 'NONE'."))
errorEstimatorKind = "NONE"
self._errorEstimatorKind = errorEstimatorKind
self._approxParameters["errorEstimatorKind"] = self.errorEstimatorKind
def _polyvanderAuxiliary(self, mus, deg, *args):
return polyvander(mus, deg, *args)
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)
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)
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.E,
self.polybasis0, self._derIdxs,
self._reorder)
- interpPQ = customFit(vanTrain, radiusAbTrain,
- rcond = self.interpRcond)
+ interpPQ = customFit(vanTrain, radiusAbTrain, rcond = self.interpTol)
vanTest = self._polyvanderAuxiliary(muCTest, self.E, self.polybasis0)
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)
mu_muTestS = mus[idxMaxEst]
app_muTestSample = self.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:
applyCglob = (hasattr(self.HFEngine, "_is_C_quadratic")
and self.HFEngine._is_C_quadratic)
for j, mu in enumerate(mu_muTestS):
uEx = self.samplingEngine.nextSample(mu)
if what == "OUTPUT" and not applyCglob:
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 = app_muTS.dtype)
app_muTestS[j] = app_muTS
if j == 0:
solmu = emptySampleList()
solmu.reset((len(uEx), len(mu_muTestS)), dtype = uEx.dtype)
solmu[j] = uEx
if what == "OUTPUT":
if applyCglob:
solmu = sampleList(self.HFEngine.applyC(solmu, mu_muTestS))
app_muTestS = sampleList(self.HFEngine.applyC(
app_muTestSample, mu_muTestS))
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.E + 1, self.polybasis)
- fitOut = customFit(vanT, errT, full = True, rcond = self.interpRcond)
+ 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.E + 1,
polyfitname(self.polybasis),
fitOut[1][2][0] / fitOut[1][2][-1]), 15)
vanC = self._polyvanderAuxiliary(musC, self.E + 1, self.polybasis)
err = np.sum(np.abs(vanC.dot(fitOut[0])), axis = -1) / QTest
return err, idxMaxEst
def getErrorEstimatorNone(self, mus:Np1D) -> Np1D:
"""EIM-based residual estimator."""
err = np.max(self._EIMStep(mus, True), axis = 1)
err *= self.greedyTol / np.mean(err)
return err
def _EIMStep(self, mus:Np1D,
only_one : bool = False) -> Tuple[Np1D, Np1D, List[int]]:
"""Residual estimator based on look-ahead idea."""
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)
muCTest = self.trainedModel.centerNormalize(mus)
muCTrain = self.trainedModel.centerNormalize(self.mus)
self.trainedModel.verbosity = tMverb
vanTest = self._polyvanderAuxiliary(muCTest, self.E, self.polybasis)
vanTestNext = self._polyvanderAuxiliary(muCTest, self.E + 1,
self.polybasis)[:,
vanTest.shape[1] :]
idxsTest = np.arange(vanTestNext.shape[1])
basis = np.zeros((len(idxsTest), 0), dtype = float)
idxMaxEst = []
while len(idxsTest) > 0:
vanTrial = self._polyvanderAuxiliary(muCTrain, self.E,
self.polybasis)
vanTrialNext = self._polyvanderAuxiliary(muCTrain, self.E + 1,
self.polybasis)[:,
vanTrial.shape[1] :]
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
idxMaxErr = np.unravel_index(np.argmax(errTest), errTest.shape)
idxMaxEst += [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
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.empty(self.sampleBatchSize, dtype = int)
errCP = copy(err)
for j in range(self.sampleBatchSize):
k = np.argmax(errCP)
idxMaxEst[j] = k
if j + 1 < self.sampleBatchSize:
musZero = self.trainedModel.centerNormalize(mus, mus[k])
errCP *= np.linalg.norm(musZero.data, axis = 1)
return err, idxMaxEst, err[idxMaxEst]
def plotEstimator(self, *args, **kwargs):
super().plotEstimator(*args, **kwargs)
if self.errorEstimatorKind == "NONE":
vbMng(self, "MAIN",
("Warning! Error estimator has been arbitrarily normalized "
"before plotting."), 15)
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.")
self.sampleBatchIdx += 1
self.sampleBatchSize = totalDegreeN(self.npar - 1, self.sampleBatchIdx)
err, muidx, maxErr, muNext = super().greedyNextSample(*args, **kwargs)
if maxErr is not None and (np.any(np.isnan(maxErr))
or np.any(np.isinf(maxErr))):
self.sampleBatchIdx -= 1
self.sampleBatchSize = totalDegreeN(self.npar - 1,
self.sampleBatchIdx)
if (self.errorEstimatorKind == "NONE" and not np.isnan(maxErr)
and not np.isinf(maxErr)):
maxErr = None
return err, muidx, maxErr, muNext
def _setSampleBatch(self, maxS:int):
self.sampleBatchIdx, self.sampleBatchSize, S = -1, 0, 0
nextBatchSize = 1
while S + nextBatchSize <= maxS:
self.sampleBatchIdx += 1
self.sampleBatchSize = nextBatchSize
S += self.sampleBatchSize
nextBatchSize = totalDegreeN(self.npar - 1,
self.sampleBatchIdx + 1)
return S
def _preliminaryTraining(self):
"""Initialize starting snapshots of solution map."""
RROMPyAssert(self._mode, message = "Cannot start greedy algorithm.")
if self.samplingEngine.nsamples > 0: return
self._S = self._setSampleBatch(self.S)
super()._preliminaryTraining()
self.M, self.N = ("AUTO",) * 2
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)
pMatOld, pMat = None, self.samplingEngine.projectionMatrix
firstRun = self.trainedModel is None
applyCglob = (hasattr(self.HFEngine, "_is_C_quadratic")
and self.HFEngine._is_C_quadratic)
if not firstRun:
Sold = len(self.trainedModel.data.mus)
if applyCglob: pMatOld = pMat[:, : Sold]
pMat = pMat[:, Sold :]
self._setupTrainedModel(pMat, not firstRun, pMatOld)
self.catchInstability = 2
vbDepth = getVerbosityDepth()
unstable = 0
if self.E > 0:
try:
Q = self._setupDenominator()
except RROMPyException as RE:
if RE.critical: raise RE from None
setVerbosityDepth(vbDepth)
RROMPyWarning("Downgraded {}: {}".format(RE.__class__.__name__,
RE))
unstable = 1
else:
Q = PI()
Q.coeffs = np.ones((1,) * self.npar, dtype = np.complex)
Q.npar = self.npar
Q.polybasis = self.polybasis
if not unstable:
self.trainedModel.data.Q = copy(Q)
try:
P = copy(self._setupNumerator())
except RROMPyException as RE:
if RE.critical: raise RE from None
setVerbosityDepth(vbDepth)
RROMPyWarning("Downgraded {}: {}".format(RE.__class__.__name__,
RE))
unstable = 1
if not unstable:
self.trainedModel.data.P = copy(P)
self.trainedModel.data.approxParameters = copy(
self.approxParameters)
vbMng(self, "DEL", "Done setting up local approximant.", 5)
self.catchInstability = 0
self.verbosity += 10
return unstable
def setupApprox(self, plotEst : str = "NONE") -> int:
val = super().setupApprox(plotEst)
if val == 0:
if (self.errorEstimatorKind[:10] == "LOOK_AHEAD"
and len(self.mus) < self.samplingEngine.nsamples):
while len(self.mus) < self.samplingEngine.nsamples:
self.mus.append(self.samplingEngine.mus[len(self.mus)])
self.trainedModel = None
self._S = self._setSampleBatch(len(self.mus) + 1)
self.setupApproxLocal()
self._setupRational(self.trainedModel.data.Q,
self.trainedModel.data.P)
self.trainedModel.data.approxParameters = copy(
self.approxParameters)
return val
def loadTrainedModel(self, filename:str):
"""Load trained reduced model from file."""
super().loadTrainedModel(filename)
self._setSampleBatch(self.S + 1)
diff --git a/rrompy/reduction_methods/standard/rational_interpolant.py b/rrompy/reduction_methods/standard/rational_interpolant.py
index fe98b4f..ca32c08 100644
--- a/rrompy/reduction_methods/standard/rational_interpolant.py
+++ b/rrompy/reduction_methods/standard/rational_interpolant.py
@@ -1,836 +1,816 @@
# 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.linalg import eigvals
from collections.abc import Iterable
from .generic_standard_approximant import GenericStandardApproximant
from rrompy.utilities.poly_fitting.polynomial import (
polybases as ppb, polyfitname,
polyvander as pvP, polyTimes,
PolynomialInterpolator as PI,
PolynomialInterpolatorNodal as PIN)
from rrompy.utilities.poly_fitting.heaviside import rational2heaviside
from rrompy.utilities.poly_fitting.radial_basis import (polybases as rbpb,
RadialBasisInterpolator as RBI)
from rrompy.utilities.base.types import (Np1D, Np2D, Tuple, List, sampList,
interpEng)
from rrompy.utilities.base import verbosityManager as vbMng
from rrompy.utilities.numerical import (pseudoInverse, dot, potential,
distanceMatrix)
from rrompy.utilities.numerical.factorials import multifactorial
from rrompy.utilities.numerical.hash_derivative import (nextDerivativeIndices,
hashDerivativeToIdx as hashD,
hashIdxToDerivative as hashI)
from rrompy.utilities.numerical.degree import (reduceDegreeN,
degreeTotalToFull, fullDegreeMaxMask,
totalDegreeMaxMask)
from rrompy.solver import Np2DLikeGramian
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)
def vanderInvTable(vanInv:Np2D, idxs:List[int], reorder:List[int],
derIdxs:List[List[List[int]]]) -> Np2D:
"""Table of Vandermonde pseudo-inverse."""
S = len(reorder)
Ts = [None] * len(idxs)
for k in range(len(idxs)):
invLocs = [None] * len(derIdxs)
idxGlob = 0
for j, derIdx in enumerate(derIdxs):
nder = len(derIdx)
idxGlob += nder
idxLoc = np.arange(S)[np.logical_and(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) -> Np2D:
"""Table of derivative values for point confluence."""
S = len(reorder)
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)[np.logical_and(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]
T[idxLoc[i1], idxLoc[i2]] = val[i3]
return T
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;
- 'polybasis': type of polynomial basis for 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;
- 'polydegreetype': type of polynomial degree; defaults to 'TOTAL';
- 'radialDirectionalWeights': radial basis weights for interpolant
numerator; defaults to 1;
- 'radialDirectionalWeightsAdapt': bounds for adaptive rescaling of
radial basis weights; defaults to [-1, -1];
- 'functionalSolve': strategy for minimization of denominator
functional; allowed values include 'NORM', 'DOMINANT', 'NODAL',
'BARYCENTRIC_NORM', and 'BARYCENTRIC[_AVERAGE]' (check pdf in
main folder for explanation); defaults to 'NORM';
- - 'interpRcond': tolerance for interpolation; defaults to None;
- - 'robustTol': tolerance for robust rational denominator
- management; defaults to 0.
+ - 'interpTol': tolerance for interpolation; defaults to None;
+ - '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;
- 'polybasis': type of polynomial basis for interpolation;
- 'M': degree of rational interpolant numerator;
- 'N': degree of rational interpolant denominator;
- 'polydegreetype': type of polynomial degree;
- 'radialDirectionalWeights': radial basis weights for interpolant
numerator;
- 'radialDirectionalWeightsAdapt': bounds for adaptive rescaling of
radial basis weights;
- 'functionalSolve': strategy for minimization of denominator
functional;
- - 'interpRcond': tolerance for interpolation via numpy.polyfit;
- - 'robustTol': tolerance for robust rational denominator
- management.
+ - 'interpTol': tolerance for interpolation via numpy.polyfit;
+ - '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.
polybasis: type of polynomial basis for interpolation.
M: Numerator degree of approximant.
N: Denominator degree of approximant.
polydegreetype: Type of polynomial degree.
radialDirectionalWeights: Radial basis weights for interpolant
numerator.
radialDirectionalWeightsAdapt: Bounds for adaptive rescaling of radial
basis weights.
functionalSolve: Strategy for minimization of denominator functional.
- interpRcond: Tolerance for interpolation via numpy.polyfit.
- robustTol: Tolerance for robust rational denominator management.
+ interpTol: Tolerance for interpolation via numpy.polyfit.
+ 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", "NODAL",
"BARYCENTRIC_NORM", "BARYCENTRIC_AVERAGE"]
def __init__(self, *args, **kwargs):
self._preInit()
self._addParametersToList(["polybasis", "M", "N", "polydegreetype",
"radialDirectionalWeights",
"radialDirectionalWeightsAdapt",
- "functionalSolve", "interpRcond",
- "robustTol"],
+ "functionalSolve", "interpTol", "QTol"],
["MONOMIAL", "AUTO", "AUTO", "TOTAL", 1.,
[-1., -1.], "NORM", -1, 0.])
super().__init__(*args, **kwargs)
self.catchInstability = 0
self._postInit()
@property
def tModelType(self):
from .trained_model.trained_model_rational import TrainedModelRational
return TrainedModelRational
@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 ppb + rbpb:
raise RROMPyException("Prescribed polybasis not recognized.")
self._polybasis = polybasis
except:
RROMPyWarning(("Prescribed polybasis not recognized. Overriding "
"to 'MONOMIAL'."))
self._polybasis = "MONOMIAL"
self._approxParameters["polybasis"] = self.polybasis
@property
def polybasis0(self):
if "_" in self.polybasis:
return self.polybasis.split("_")[0]
return 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 += "_AVERAGE"
if functionalSolve not in self._allowedFunctionalSolveKinds:
raise RROMPyException(("Prescribed functionalSolve not "
"recognized."))
self._functionalSolve = functionalSolve
except:
RROMPyWarning(("Prescribed functionalSolve not recognized. "
"Overriding to 'NORM'."))
self._functionalSolve = "NORM"
self._approxParameters["functionalSolve"] = self.functionalSolve
@property
- def interpRcond(self):
- """Value of interpRcond."""
- return self._interpRcond
- @interpRcond.setter
- def interpRcond(self, interpRcond):
- self._interpRcond = interpRcond
- self._approxParameters["interpRcond"] = self.interpRcond
+ 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 radialDirectionalWeights(self):
"""Value of radialDirectionalWeights."""
return self._radialDirectionalWeights
@radialDirectionalWeights.setter
def radialDirectionalWeights(self, radialDirectionalWeights):
if isinstance(radialDirectionalWeights, Iterable):
radialDirectionalWeights = list(radialDirectionalWeights)
else:
radialDirectionalWeights = [radialDirectionalWeights]
self._radialDirectionalWeights = radialDirectionalWeights
self._approxParameters["radialDirectionalWeights"] = (
self.radialDirectionalWeights)
@property
def radialDirectionalWeightsAdapt(self):
"""Value of radialDirectionalWeightsAdapt."""
return self._radialDirectionalWeightsAdapt
@radialDirectionalWeightsAdapt.setter
def radialDirectionalWeightsAdapt(self, radialDirectionalWeightsAdapt):
self._radialDirectionalWeightsAdapt = radialDirectionalWeightsAdapt
self._approxParameters["radialDirectionalWeightsAdapt"] = (
self.radialDirectionalWeightsAdapt)
@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):
self.M = max(0, reduceDegreeN(self.S, self.S, self.npar,
self.polydegreetype) - 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):
self.N = max(0, reduceDegreeN(self.S, self.S, self.npar,
self.polydegreetype) - self._N_shift)
vbMng(self, "MAIN", "Automatically setting N to {}.".format(self.N),
25)
@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", "FULL"]:
raise RROMPyException(("Prescribed polydegreetype not "
"recognized."))
self._polydegreetype = polydegreetype
except:
RROMPyWarning(("Prescribed polydegreetype not recognized. "
"Overriding to 'TOTAL'."))
self._polydegreetype = "TOTAL"
self._approxParameters["polydegreetype"] = self.polydegreetype
@property
- def robustTol(self):
+ def QTol(self):
"""Value of tolerance for robust rational denominator management."""
- return self._robustTol
- @robustTol.setter
- def robustTol(self, robustTol):
- if robustTol < 0.:
+ return self._QTol
+ @QTol.setter
+ def QTol(self, QTol):
+ if QTol < 0.:
RROMPyWarning(("Overriding prescribed negative robustness "
"tolerance to 0."))
- robustTol = 0.
- self._robustTol = robustTol
- self._approxParameters["robustTol"] = self.robustTol
+ 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 hasattr(self, "_N_isauto"):
self._setNAuto()
else:
N = reduceDegreeN(self.N, self.S, self.npar, self.polydegreetype)
if N < self.N:
RROMPyWarning(("N too large compared to S. Reducing N by "
"{}").format(self.N - N))
self.N = N
while self.N > 0:
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.N + 1 < self.S):
RROMPyWarning(("Barycentric strategy cannot be applied with "
"Least Squares. Overriding to 'NORM'."))
self.functionalSolve = "NORM"
if self.functionalSolve[:11] == "BARYCENTRIC":
invD, TN = None, None
self._setupInterpolationIndices()
else:
invD, TN = self._computeInterpolantInverseBlocks()
if (self.functionalSolve in ["NODAL", "BARYCENTRIC_NORM",
"BARYCENTRIC_AVERAGE"]
- and len(self._musUnique) != len(self.mus)):
+ and len(self._musUnique) != self.S):
if self.functionalSolve[:11] == "BARYCENTRIC":
warnflag = "Barycentric"
invD, TN = self._computeInterpolantInverseBlocks()
else:
warnflag = "Iterative"
RROMPyWarning(("{} functional optimization cannot be applied "
"to repeated samples. Overriding to "
"'NORM'.").format(warnflag))
self.functionalSolve = "NORM"
idxSamplesEff = list(range(self.S))
if self.POD == 1:
ev, eV = self.findeveVGQR(
self.samplingEngine.Rscale[:, idxSamplesEff], invD, TN)
elif self.POD == 1/2:
ev, eV = self.findeveVGExplicit(
self.samplingEngine.samples_normal(idxSamplesEff),
invD, TN, self.samplingEngine.Rscale)
else:
ev, eV = self.findeveVGExplicit(
self.samplingEngine.samples(idxSamplesEff), invD, TN)
if self.functionalSolve == "NODAL": break
evR = ev / np.max(ev)
- ts = self.robustTol * np.linalg.norm(evR)
+ ts = self.QTol * np.linalg.norm(evR)
nevBad = len(ev) - np.sum(np.abs(evR) >= ts)
if nevBad <= (self.functionalSolve == "NORM"): break
if self.catchInstability:
raise RROMPyException(("Instability in denominator "
"computation: eigenproblem is poorly "
"conditioned."),
self.catchInstability == 1)
vbMng(self, "MAIN", ("Smallest {} eigenvalues below tolerance. "
"Reducing N by 1.").format(nevBad), 10)
self.N = self.N - 1
if self.N <= 0:
self.N = 0
eV = np.ones((1, 1))
if self.N > 0 and self.functionalSolve in ["NODAL", "BARYCENTRIC_NORM",
"BARYCENTRIC_AVERAGE"]:
q = PIN()
q.polybasis, q.nodes = self.polybasis0, eV
else:
q = PI()
q.npar = self.npar
q.polybasis = self.polybasis0
if self.polydegreetype == "TOTAL":
q.coeffs = degreeTotalToFull(tuple([self.N + 1] * self.npar),
self.npar, eV)
else:
q.coeffs = eV.reshape([self.N + 1] * self.npar)
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.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 = self.M
else:
M = reduceDegreeN(self.M, self.S, self.npar, self.polydegreetype)
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:
pParRest = [self.M, self.polybasis, self.verbosity >= 5,
self.polydegreetype == "TOTAL",
{"derIdxs": self._derIdxs, "reorder": self._reorder,
"scl": self.scaleFactorRel}]
if self.polybasis in ppb:
p = PI()
else:
self.computeScaleFactor()
rDWEff = np.array([w * f for w, f in zip(
self.radialDirectionalWeights,
self.scaleFactor)])
pParRest = pParRest[: 2] + [rDWEff] + pParRest[2 :]
pParRest[-1]["optimizeScalingBounds"] = (
self.radialDirectionalWeightsAdapt)
p = RBI()
if self.polybasis in ppb + rbpb:
- pParRest += [{"rcond": self.interpRcond}]
+ pParRest += [{"rcond": self.interpTol}]
wellCond, msg = p.setupByInterpolation(self._musUniqueCN,
Qevaldiag, *pParRest)
vbMng(self, "MAIN", msg, 5)
if wellCond: break
if self.catchInstability:
raise RROMPyException(("Instability in numerator computation: "
"polyfit is poorly conditioned."),
self.catchInstability == 1)
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."))
self.M = M
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()
self._setupTrainedModel(self.samplingEngine.projectionMatrix)
self._setupRational(self._setupDenominator())
self.trainedModel.data.approxParameters = copy(self.approxParameters)
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()
- if self.N > 0 and self.npar == 1:
- #check for bad poles
+ while self.N > 0 and self.npar == 1:
+ #check for bad poles (absolute)
pls = Q.roots()
- idxBad = self.HFEngine.flagBadPolesResidues(pls, relative = True)
- plsN = self.mapParameterList(self.mapParameterList(self.mu0)(0, 0)
- + self.scaleFactor * pls, "B")(0)
- idxBad = np.logical_or(self.HFEngine.flagBadPolesResidues(pls,
- relative = True),
- self.HFEngine.flagBadPolesResidues(plsN))
- if np.any(idxBad):
- vbMng(self, "MAIN",
- "Removing {} spurious poles out of {} due to poles."\
- .format(np.sum(idxBad), self.N), 10)
- if isinstance(Q, PIN):
- Q.nodes = Q.nodes[np.logical_not(idxBad)]
- else:
- Q = PI()
- Q.npar = self.npar
- Q.polybasis = self.polybasis0
- Q.coeffs = np.ones(1, dtype = np.complex)
- for pl in pls[np.logical_not(idxBad)]:
- 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()
- if (not hasattr(self.HFEngine, "_ignoreResidues")
- or not self.HFEngine._ignoreResidues):
- #check for bad residues
- cfs, pls, _ = rational2heaviside(P, Q)
- cfs = cfs[: self.N].T
+ plsA = self.mapParameterList(self.mapParameterList(self.mu0)(0, 0)
+ + self.scaleFactor * pls, "B")(0)
+ idxBadA = self.HFEngine.flagBadPolesResidues(plsA)
+ if self.HFEngine._ignoreResidues:
+ #check for bad poles (relative)
+ cfs, projMat = None, None
+ else:
+ #check for bad poles (relative) and residues
+ cfs = rational2heaviside(P, Q)[0][: self.N].T
if self.POD != 1:
- cfs = self.samplingEngine.projectionMatrix.dot(cfs)
+ projMat = self.samplingEngine.projectionMatrix
+ else:
+ projMat = None
foci = self.sampler.normalFoci()
ground = self.sampler.groundPotential()
potEff = potential(pls, foci) / ground
potEff[np.logical_or(potEff < 1., np.isinf(pls))] = 1.
cfs[:, np.isinf(pls)] = 0.
cfs /= potEff # rescale by potential
- idxBad = self.HFEngine.flagBadPolesResidues(pls, cfs)
- if np.any(idxBad):
- vbMng(self, "MAIN",
- ("Removing {} spurious poles out of {} due to "
- "residues.").format(np.sum(idxBad), self.N), 10)
- if isinstance(Q, PIN):
- Q.nodes = Q.nodes[np.logical_not(idxBad)]
- else:
- Q = PI()
- Q.npar = self.npar
- Q.polybasis = self.polybasis0
- Q.coeffs = np.ones(1, dtype = np.complex)
- for pl in pls[np.logical_not(idxBad)]:
- Q.coeffs = polyTimes(Q.coeffs, [- pl, 1.],
+ idxBadR = self.HFEngine.flagBadPolesResidues(pls, cfs, True,
+ projMat)
+ idxBad = np.logical_or(idxBadA, idxBadR)
+ if not np.any(idxBad): break
+ vbMng(self, "MAIN",
+ "Removing {} spurious poles out of {}.".format(
+ np.sum(idxBad), self.N), 10)
+ if isinstance(Q, PIN):
+ Q.nodes = Q.nodes[np.logical_not(idxBad)]
+ else:
+ Q = PI()
+ Q.npar = self.npar
+ Q.polybasis = self.polybasis0
+ Q.coeffs = np.ones(1, dtype = np.complex)
+ for pl in pls[np.logical_not(idxBad)]:
+ 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()
+ 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 minimal interpolant target functional.
"""
RROMPyAssert(self._mode, message = "Cannot solve eigenvalue problem.")
self._setupInterpolationIndices()
pvPPar = [self.polybasis0, self._derIdxs, self._reorder,
self.scaleFactorRel]
if hasattr(self, "_M_isauto"): self._setMAuto()
E = max(self.M, self.N)
fullE = E + 1 == len(self._reorder) == len(self._musUniqueCN)
if fullE:
mus = self._musUniqueCN[self._reorder]
dist = distanceMatrix(mus, magnitude = False)[..., 0]
dist[np.arange(E + 1), np.arange(E + 1)] = multifactorial([E])
fitinvE = np.prod(dist, axis = 1) ** -1
vbMng(self, "MAIN",
("Evaluating quasi-Lagrangian basis of degree {} at {} "
"sample points.").format(E, E + 1), 5)
invD = [np.diag(fitinvE)]
else:
while E >= 0:
if self.polydegreetype == "TOTAL":
Eeff = E
idxsB = totalDegreeMaxMask(E, self.npar)
else: #if self.polydegreetype == "FULL":
Eeff = [E] * self.npar
idxsB = fullDegreeMaxMask(E, self.npar)
TE = pvP(self._musUniqueCN, Eeff, *pvPPar)
- fitOut = pseudoInverse(TE, rcond = self.interpRcond,
- full = True)
+ fitOut = pseudoInverse(TE, rcond = self.interpTol, full = True)
vbMng(self, "MAIN",
("Fitting {} samples with degree {} through {}... "
"Conditioning of pseudoinverse system: {:.4e}.").format(
TE.shape[0], E,
polyfitname(self.polybasis0),
fitOut[1][1][0] / fitOut[1][1][-1]), 5)
if fitOut[1][0] == TE.shape[1]:
fitinv = fitOut[0][idxsB, :]
break
if self.catchInstability:
raise RROMPyException(("Instability in denominator "
"computation: polyfit is poorly "
"conditioned."),
self.catchInstability == 1)
EeqN = E == self.N
vbMng(self, "MAIN", ("Polyfit is poorly conditioned. Reducing "
"E {} by 1.").format("and N " * EeqN), 10)
if EeqN: self.N = self.N - 1
E -= 1
if self.N < 0:
raise RROMPyException(("Instability in computation of "
"denominator. Aborting."))
invD = vanderInvTable(fitinv, idxsB, self._reorder, self._derIdxs)
if self.N == E and not fullE:
TN = TE
else: #build TN from scratch
if self.polydegreetype == "TOTAL":
Neff = self.N
idxsB = totalDegreeMaxMask(self.N, self.npar)
else: #if self.polydegreetype == "FULL":
Neff = [self.N] * self.npar
idxsB = fullDegreeMaxMask(self.N, self.npar)
TN = pvP(self._musUniqueCN, Neff, *pvPPar)
return invD, TN
def findeveVGExplicit(self, sampleE:sampList, invD:List[Np2D], TN:Np2D,
Rscaling : Np1D = None) -> Tuple[Np1D, Np2D]:
"""
Compute explicitly eigenvalues and eigenvectors of rational denominator
matrix.
"""
RROMPyAssert(self._mode, message = "Cannot solve eigenvalue problem.")
vbMng(self, "INIT", "Building gramian matrix.", 10)
gramian = self.HFEngine.innerProduct(sampleE, sampleE, is_state = True)
if Rscaling is not None:
gramian = (gramian.T * Rscaling.conj()).T * Rscaling
if self.functionalSolve == "NODAL":
SEnd = invD[0].shape[1]
G0 = np.zeros((SEnd,) * 2, dtype = np.complex)
if self.functionalSolve[:11] == "BARYCENTRIC":
G = gramian
nEnd = len(gramian) - 1
else:
nEnd = TN.shape[1]
G = np.zeros((nEnd, nEnd), dtype = np.complex)
for k in range(len(invD)):
iDkN = dot(invD[k], TN)
G += dot(dot(gramian, iDkN).T, iDkN.conj()).T
if self.functionalSolve == "NODAL":
G0 += dot(dot(gramian, invD[k]).T, invD[k].conj()).T
vbMng(self, "DEL", "Done building gramian.", 10)
if self.functionalSolve in ["NORM", "BARYCENTRIC_NORM"]:
ev, eV = np.linalg.eigh(G)
eV = eV[:, 0]
if self.functionalSolve == "BARYCENTRIC_NORM":
eV = self.findeVBarycentric(eV)
problem = "eigenproblem"
else:
if self.functionalSolve == "BARYCENTRIC_AVERAGE":
- fitOut = pseudoInverse(G, rcond = self.interpRcond,
- full = True)
+ fitOut = pseudoInverse(G, rcond = self.interpTol, full = True)
eV = self.findeVBarycentric(np.sum(fitOut[0], axis = 1))
else:
- fitOut = pseudoInverse(G[:-1, :-1], rcond = self.interpRcond,
+ fitOut = pseudoInverse(G[:-1, :-1], rcond = self.interpTol,
full = True)
eV = np.append(fitOut[0].dot(G[:-1, -1]), -1.)
ev = fitOut[1][1][::-1]
problem = "linear problem"
vbMng(self, "MAIN",
("Solved {} of size {} with condition number "
"{:.4e}.").format(problem, nEnd, ev[-1] / ev[0]), 5)
if self.functionalSolve == "NODAL": eV = self.findeVNodal(eV, G0)
return ev, eV
def findeveVGQR(self, RPODE:Np2D, invD:List[Np2D],
TN:Np2D) -> Tuple[Np1D, Np2D]:
"""
Compute eigenvalues and eigenvectors of rational denominator matrix
through SVD of R factor.
"""
RROMPyAssert(self._mode, message = "Cannot solve eigenvalue problem.")
vbMng(self, "INIT", "Building half-gramian matrix stack.", 10)
if self.functionalSolve == "NODAL":
gramian = Np2DLikeGramian(None, dot(RPODE, invD[0]))
if self.functionalSolve[:11] == "BARYCENTRIC":
Rstack = RPODE
nEnd = RPODE.shape[1] - 1
else:
S, nEnd, eWidth = RPODE.shape[0], TN.shape[1], len(invD)
Rstack = np.zeros((S * eWidth, nEnd), dtype = np.complex)
for k in range(eWidth):
Rstack[k * S : (k + 1) * S, :] = dot(RPODE, dot(invD[k], TN))
vbMng(self, "DEL", "Done building half-gramian.", 10)
if self.functionalSolve in ["NORM", "BARYCENTRIC_NORM",
"BARYCENTRIC_AVERAGE"]:
_, ev, Vh = np.linalg.svd(Rstack, full_matrices = False)
if self.functionalSolve in ["NORM", "BARYCENTRIC_NORM"]:
eV = Vh[-1, :].conj()
ev = ev[::-1]
else: #if self.functionalSolve == "BARYCENTRIC_AVERAGE":
ev[np.logical_not(np.isclose(ev, 0.))] **= -2.
eV = Vh.T.conj().dot(ev * np.sum(Vh, axis = 1))
if self.functionalSolve[:11] == "BARYCENTRIC":
eV = self.findeVBarycentric(eV)
problem = "svd problem"
else:
- fitOut = pseudoInverse(Rstack[:, :-1], rcond = self.interpRcond,
+ fitOut = pseudoInverse(Rstack[:, :-1], rcond = self.interpTol,
full = True)
ev = fitOut[1][1][::-1]
eV = np.append(fitOut[0].dot(Rstack[:, -1]), -1.)
problem = "linear problem"
vbMng(self, "MAIN",
("Solved {} of size {} with condition number "
"{:.4e}.").format(problem, nEnd, ev[-1] / ev[0]), 5)
if self.functionalSolve == "NODAL": eV = self.findeVNodal(eV, gramian)
return ev, eV
def findeVBarycentric(self, baryWeights:Np1D) -> Np1D:
RROMPyAssert(self._mode,
message = "Cannot solve optimization problem.")
arrow = np.pad(np.diag(self._musUniqueCN[self._reorder].flatten()),
(1, 0), "constant", constant_values = 1.) + 0.j
arrow[0, 0] = 0.
arrow[0, 1:] = baryWeights
active = np.pad(np.eye(len(baryWeights)), (1, 0), "constant")
eV = eigvals(arrow, active)
return eV[np.logical_not(np.isinf(eV))]
def findeVNodal(self, q0coeffs:Np1D, gram:Np2D, maxiter : int = 25,
tolNewton : float = 1e-10) -> Np1D:
RROMPyAssert(self._mode,
message = "Cannot solve optimization problem.")
q = PI()
q.npar, q.polybasis, q.coeffs = self.npar, self.polybasis0, q0coeffs
nodes = q.roots()
N = len(nodes)
grad = np.zeros(N, dtype = np.complex)
hess = np.zeros((N, N), dtype = np.complex)
for niter in range(maxiter):
plDist = distanceMatrix(self._musUniqueCN[self._reorder], nodes,
magnitude = False)[:, :, 0]
q0, q = np.prod(plDist, axis = 1), []
for jS in range(N):
mask = np.arange(N) != jS
q += [np.prod(plDist[:, mask], axis = 1)]
grad[jS] = q[-1].conj().dot(gram.dot(q0))
for iS in range(jS + 1):
if iS == jS:
hij = 0.
sij = q[-1].conj().dot(gram.dot(q[-1]))
else:
mask = np.logical_and(np.arange(N) != iS,
np.arange(N) != jS)
qij = np.prod(plDist[:, mask], axis = 1)
hij = qij.conj().dot(gram.dot(q0))
sij = q[-1].conj().dot(gram.dot(q[iS]))
hess[jS, iS] = hij + sij
if iS < jS: hess[iS, jS] = hij + np.conj(sij)
dnodes = np.linalg.solve(hess, grad)
nodes += dnodes
tol = np.linalg.norm(dnodes) / np.linalg.norm(nodes)
if tol < tolNewton: break
if niter < maxiter:
vbMng(self, "MAIN",
("Newton's method for problem of size {} converged "
"(err = {:.4e}) in {} iteration{}.").format(
N + 1, tol, niter + 1, "s" * (niter > 0)), 5)
else:
RROMPyWarning(("Newton's method for problem of size {} did "
"not converge (err = {:.4e}) after {} "
"iterations.").format(N + 1, tol, niter + 1))
return nodes
def getResidues(self, *args, **kwargs) -> 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_rational.py b/rrompy/reduction_methods/standard/trained_model/trained_model_rational.py
index 543a05d..eb055c8 100644
--- a/rrompy/reduction_methods/standard/trained_model/trained_model_rational.py
+++ b/rrompy/reduction_methods/standard/trained_model/trained_model_rational.py
@@ -1,195 +1,197 @@
# 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.utilities.numerical import dot
from rrompy.utilities.numerical.compress_matrix import compressMatrix
from rrompy.utilities.base.types import (Np1D, Np2D, List, paramVal, paramList,
sampList)
from rrompy.utilities.base import verbosityManager as vbMng, freepar as fp
from rrompy.utilities.exception_manager import RROMPyException, RROMPyWarning
from rrompy.parameter import emptyParameterList
from rrompy.sampling import sampleList
__all__ = ['TrainedModelRational']
class TrainedModelRational(TrainedModel):
"""
ROM approximant evaluation for rational approximant.
Attributes:
Data: dictionary with all that can be pickled.
"""
def compress(self, collapse : bool = False, tol : float = 0., *args,
**kwargs):
if not collapse and tol <= 0.: return
if hasattr(self.data, "_is_C_quadratic") and self.data._is_C_quadratic:
raise RROMPyException(("Cannot compress model with quadratic "
"output."))
RMat = self.data.projMat
if not collapse:
if hasattr(self.data, "_compressTol"):
RROMPyWarning(("Recompressing already compressed model is "
"ineffective. Aborting."))
return
self.data.projMat, RMat, _ = compressMatrix(RMat, tol, *args,
**kwargs)
self.data.P.postmultiplyTensorize(RMat.T)
super().compress(collapse, tol)
def centerNormalize(self, mu : paramList = [],
mu0 : paramVal = None) -> paramList:
"""
Compute normalized parameter to be plugged into approximant.
Args:
mu: Parameter(s) 1.
mu0: Parameter(s) 2. If None, set to self.data.mu0.
Returns:
Normalized parameter.
"""
mu = self.checkParameterList(mu)
if mu0 is None: mu0 = self.data.mu0
return (self.mapParameterList(mu)
- self.mapParameterList(mu0)) / self.data.scaleFactor
def getPVal(self, mu : paramList = []) -> sampList:
"""
Evaluate rational numerator at arbitrary parameter.
Args:
mu: Target parameter.
"""
mu = self.checkParameterList(mu)
vbMng(self, "INIT", "Evaluating numerator at mu = {}.".format(mu), 17)
p = sampleList(self.data.P(self.centerNormalize(mu)))
vbMng(self, "DEL", "Done evaluating numerator.", 17)
return p
def getQVal(self, mu:Np1D, der : List[int] = None,
scl : Np1D = None) -> Np1D:
"""
Evaluate rational denominator at arbitrary parameter.
Args:
mu: Target parameter.
der(optional): Derivatives to take before evaluation.
"""
mu = self.checkParameterList(mu)
vbMng(self, "INIT", "Evaluating denominator at mu = {}.".format(mu),
17)
q = self.data.Q(self.centerNormalize(mu), der, scl)
vbMng(self, "DEL", "Done evaluating denominator.", 17)
return q
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",
"Evaluating approximant at mu = {}.".format(mu), 12)
QV = self.getQVal(mu)
QVzero = np.where(QV == 0.)[0]
if len(QVzero) > 0:
QV[QVzero] = np.finfo(np.complex).eps / (1.
+ self.data.Q.deg[0])
self.uApproxReduced = self.getPVal(mu) / QV
vbMng(self, "DEL", "Done evaluating approximant.", 12)
self.lastSolvedApproxReduced = mu
return self.uApproxReduced
def getPoles(self, *args, **kwargs) -> Np1D:
"""
Obtain approximant poles.
Returns:
Numpy complex vector of poles.
"""
if len(args) + len(kwargs) > 1:
raise RROMPyException(("Wrong number of parameters passed. "
"Only 1 available."))
elif len(args) + len(kwargs) == 1:
if len(args) == 1:
mVals = args[0]
else:
mVals = kwargs["marginalVals"]
if not isinstance(mVals, Iterable): mVals = [mVals]
mVals = list(mVals)
else:
mVals = [fp]
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."))
mVals[rDim] = self.data.mu0(rDim)
mVals = list(self.centerNormalize(mVals).data.flatten())
mVals[rDim] = fp
roots = self.data.scaleFactor[rDim] * self.data.Q.roots(mVals)
return self.mapParameterList(self.mapParameterList(self.data.mu0(rDim),
idx = [rDim])(0, 0)
+ roots, "B", [rDim])(0)
def getResidues(self, *args, **kwargs) -> Np2D:
"""
Obtain approximant residues.
Returns:
Numpy matrix with residues as columns.
"""
pls = self.getPoles(*args, **kwargs)
+ if len(pls) == 0:
+ return pls, np.empty((0, 0), dtype = self.data.P.coeffs.dtype)
if len(args) == 1:
mVals = args[0]
elif len(args) == 0:
mVals = [None]
else:
mVals = kwargs["marginalVals"]
if not isinstance(mVals, Iterable): mVals = [mVals]
mVals = list(mVals)
rDim = mVals.index(fp)
poles = emptyParameterList()
poles.reset((len(pls), self.data.npar), dtype = pls.dtype)
for k, pl in enumerate(pls):
mValsLoc = list(mVals)
mValsLoc[rDim] = pl
poles[k] = mValsLoc
QV = self.getQVal(poles, list(1 * (np.arange(self.data.npar) == rDim)))
QVzero = np.where(QV == 0.)[0]
if len(QVzero) > 0:
RROMPyWarning(("Adjusting residuals to avoid division by "
"numerically zero denominator."))
QV[QVzero] = np.finfo(np.complex).eps / (1. + self.data.Q.deg[0])
res = self.getPVal(poles).data / QV
if hasattr(self.data, "_is_C_quadratic") and self.data._is_C_quadratic:
self._setupQuadMapping()
res = res[self.quad_mapping[0]] * res[self.quad_mapping[1]].conj()
if not self.data._collapsed: res = dot(self.data.projMat, res).T
if (hasattr(self.data, "_is_C_quadratic")
and self.data._is_C_quadratic == 2):
res = np.real(res)
return pls, res
diff --git a/tests/3_reduction_methods_1D/rational_interpolant_1d.py b/tests/3_reduction_methods_1D/rational_interpolant_1d.py
index c747b0b..32711b5 100644
--- a/tests/3_reduction_methods_1D/rational_interpolant_1d.py
+++ b/tests/3_reduction_methods_1D/rational_interpolant_1d.py
@@ -1,67 +1,67 @@
# 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 matrix_fft import matrixFFT
from rrompy.reduction_methods import RationalInterpolant as RI
from rrompy.parameter.parameter_sampling import (QuadratureSampler as QS,
ManualSampler as MS)
from rrompy.parameter import checkParameterList
def test_monomials(capsys):
mu = 1.5
solver = matrixFFT()
- params = {"POD": False, "S": 10, "robustTol": 1e-6, "interpRcond": 1e-3,
+ params = {"POD": False, "S": 10, "QTol": 1e-6, "interpTol": 1e-3,
"polybasis": "MONOMIAL", "sampler": QS([1.5, 6.5], "UNIFORM")}
approx = RI(solver, 4., approxParameters = params, verbosity = 10)
approx.setupApprox()
out, err = capsys.readouterr()
assert "below tolerance. Reducing N " in out
assert "poorly conditioned. Reducing M " in out
assert len(err) == 0
assert np.isclose(approx.normErr(mu)[0], 1.4746e-05, atol = 1e-4)
def test_well_cond():
mu = 1.5
solver = matrixFFT()
- params = {"POD": True, "S": 10, "robustTol": 1e-14, "interpRcond": 1e-10,
+ params = {"POD": True, "S": 10, "QTol": 1e-14, "interpTol": 1e-10,
"polybasis": "CHEBYSHEV", "sampler": QS([1., 7.], "CHEBYSHEV")}
approx = RI(solver, 4., approxParameters = params, verbosity = 0)
approx.setupApprox()
poles = approx.getPoles()
for lambda_ in np.arange(1, 8):
assert np.isclose(np.min(np.abs(poles - lambda_)), 0., atol = 1e-4)
for mu in approx.mus:
assert np.isclose(approx.normErr(mu)[0], 0., atol = 1e-8)
def test_hermite():
mu = 1.5
solver = matrixFFT()
sampler0 = QS([1., 7.], "CHEBYSHEV")
points = checkParameterList(np.tile(sampler0.generatePoints(4)(0), 3))
params = {"POD": True, "S": 12, "polybasis": "CHEBYSHEV",
"sampler": MS([1., 7.], points = points)}
approx = RI(solver, 4., approxParameters = params, verbosity = 0)
approx.setupApprox()
poles = approx.getPoles()
for lambda_ in np.arange(1, 8):
assert np.isclose(np.min(np.abs(poles - lambda_)), 0., atol = 1e-4)
for mu in approx.mus:
assert np.isclose(approx.normErr(mu)[0], 0., atol = 1e-8)
diff --git a/tests/4_reduction_methods_multiD/pivoted_rational_2d.py b/tests/4_reduction_methods_multiD/pivoted_rational_2d.py
index 15c2c35..2a65e0a 100644
--- a/tests/4_reduction_methods_multiD/pivoted_rational_2d.py
+++ b/tests/4_reduction_methods_multiD/pivoted_rational_2d.py
@@ -1,112 +1,112 @@
# 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 matrix_random import matrixRandom
from rrompy.reduction_methods import (
RationalInterpolantPivotedPoleMatch as RIP,
RationalInterpolantGreedyPivotedPoleMatch as RIGP)
from rrompy.parameter.parameter_sampling import (QuadratureSampler as QS,
ManualSampler as MS)
def test_pivoted_uniform():
mu = [5.05, 7.1]
mu0 = [5., 7.]
solver = matrixRandom()
uh = solver.solve(mu)[0]
params = {"POD": True, "S": 5, "polybasis": "CHEBYSHEV",
"samplerPivot": QS([4.75, 5.25], "CHEBYSHEV"), "SMarginal": 5,
"polybasisMarginal": "MONOMIAL", "matchingWeight": 1.,
"samplerMarginal": QS([6.75, 7.25], "UNIFORM")}
approx = RIP([0], solver, mu0, approxParameters = params, verbosity = 0)
approx.setupApprox()
uhP1 = approx.getApprox(mu)[0]
errP = approx.getErr(mu)[0]
errNP = approx.normErr(mu)[0]
myerrP = uhP1 - uh
assert np.allclose(np.abs(errP - myerrP), 0., rtol = 1e-3)
assert np.isclose(solver.norm(errP), errNP, rtol = 1e-3)
resP = approx.getRes(mu)[0]
resNP = approx.normRes(mu)
assert np.isclose(solver.norm(resP), resNP, rtol = 1e-3)
assert np.allclose(np.abs(resP - (solver.b(mu) - solver.A(mu).dot(uhP1))),
0., rtol = 1e-3)
assert np.isclose(errNP / solver.norm(uh), 6.0631706e-04, rtol = 1)
def test_pivoted_manual_grid(capsys):
mu = [5.05, 7.1]
mu0 = [5., 7.]
solver = matrixRandom()
uh = solver.solve(mu)[0]
params = {"POD": False, "S": 5, "polybasis": "MONOMIAL",
"samplerPivot": MS([4.75, 5.25], np.array([5.]),
normalFoci = [0., 0.]), "SMarginal": 5,
"polybasisMarginal": "MONOMIAL", "matchingWeight": 1.,
"samplerMarginal": MS([6.75, 7.25], np.linspace(6.75, 7.25, 5)),
- "robustTol": 1e-6, "interpRcond": 1e-3}
+ "QTol": 1e-6, "interpTol": 1e-3}
approx = RIP([0], solver, mu0, approxParameters = params, verbosity = 0)
approx.setupApprox()
uhP1 = approx.getApprox(mu)[0]
errP = approx.getErr(mu)[0]
errNP = approx.normErr(mu)[0]
myerrP = uhP1 - uh
assert np.allclose(np.abs(errP - myerrP), 0., rtol = 1e-3)
assert np.isclose(solver.norm(errP), errNP, rtol = 1e-3)
resP = approx.getRes(mu)[0]
resNP = approx.normRes(mu)
assert np.isclose(solver.norm(resP), resNP, rtol = 1e-3)
assert np.allclose(np.abs(resP - (solver.b(mu) - solver.A(mu).dot(uhP1))),
0., rtol = 1e-3)
assert np.isclose(errNP / solver.norm(uh), .4763489, rtol = 1)
out, err = capsys.readouterr()
assert ("poorly conditioned" not in out)
assert len(err) == 0
def test_pivoted_greedy():
mu = [5.05, 7.1]
mu0 = [5., 7.]
solver = matrixRandom()
uh = solver.solve(mu)[0]
params = {"POD": True, "nTestPoints": 100, "greedyTol": 1e-4,
"collinearityTol": 1e8, "errorEstimatorKind": "DISCREPANCY",
"S": 5, "polybasis": "CHEBYSHEV",
"samplerPivot": QS([4.75, 5.25], "UNIFORM"),
"trainSetGenerator": QS([4.75, 5.25], "CHEBYSHEV"),
"SMarginal": 5, "polybasisMarginal": "MONOMIAL",
"samplerMarginal": QS([6.75, 7.25], "UNIFORM"),
"matchingWeight": 1.}
solver.cutOffPolesRMinRel, solver.cutOffPolesRMaxRel = -3., 3.
solver.cutOffPolesIMinRel, solver.cutOffPolesIMaxRel = -1.5, 1.5
approx = RIGP([0], solver, mu0, approxParameters = params, verbosity = 0)
approx.setupApprox()
uhP1 = approx.getApprox(mu)[0]
errP = approx.getErr(mu)[0]
errNP = approx.normErr(mu)[0]
myerrP = uhP1 - uh
assert np.allclose(np.abs(errP - myerrP), 0., rtol = 1e-3)
assert np.isclose(solver.norm(errP), errNP, rtol = 1e-3)
resP = approx.getRes(mu)[0]
resNP = approx.normRes(mu)
assert np.isclose(solver.norm(resP), resNP, rtol = 1e-3)
assert np.allclose(np.abs(resP - (solver.b(mu) - solver.A(mu).dot(uhP1))),
0., rtol = 1e-3)
assert np.isclose(errNP / solver.norm(uh), 1.181958e-02, rtol = 1)
diff --git a/tests/4_reduction_methods_multiD/rational_interpolant_2d.py b/tests/4_reduction_methods_multiD/rational_interpolant_2d.py
index 2f8242f..5e2cb73 100644
--- a/tests/4_reduction_methods_multiD/rational_interpolant_2d.py
+++ b/tests/4_reduction_methods_multiD/rational_interpolant_2d.py
@@ -1,77 +1,77 @@
# 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 matrix_random import matrixRandom
from rrompy.reduction_methods import RationalInterpolant as RI
from rrompy.parameter.parameter_sampling import (QuadratureSampler as QS,
ManualSampler as MS)
def test_monomials(capsys):
mu = [5.05, 7.1]
mu0 = [5., 7.]
solver = matrixRandom()
uh = solver.solve(mu)[0]
- params = {"POD": False, "S": 16, "robustTol": 1e-6,
- "interpRcond": 1e-3, "polybasis": "MONOMIAL",
+ params = {"POD": False, "S": 16, "QTol": 1e-6, "interpTol": 1e-3,
+ "polybasis": "MONOMIAL",
"sampler": QS([[4.9, 6.85], [5.1, 7.15]], "UNIFORM")}
approx = RI(solver, mu0, params, verbosity = 100)
approx.setupApprox()
uhP1 = approx.getApprox(mu)[0]
errP = approx.getErr(mu)[0]
errNP = approx.normErr(mu)[0]
myerrP = uhP1 - uh
assert np.allclose(np.abs(errP - myerrP), 0., rtol = 1e-3)
assert np.isclose(solver.norm(errP), errNP, rtol = 1e-3)
resP = approx.getRes(mu)[0]
resNP = approx.normRes(mu)
assert np.isclose(solver.norm(resP), resNP, rtol = 1e-3)
assert np.allclose(np.abs(resP - (solver.b(mu) - solver.A(mu).dot(uhP1))),
0., rtol = 1e-3)
assert np.isclose(errNP / solver.norm(uh), 5.2667e-05, rtol = 1)
out, err = capsys.readouterr()
assert ("poorly conditioned. Reducing E " in out)
assert len(err) == 0
def test_well_cond():
mu = [5.05, 7.1]
mu0 = [5., 7.]
solver = matrixRandom()
params = {"POD": True, "M": 3, "N": 3, "S": 16,
- "interpRcond": 1e-10, "polybasis": "CHEBYSHEV",
+ "interpTol": 1e-10, "polybasis": "CHEBYSHEV",
"sampler": QS([[4.9, 6.85], [5.1, 7.15]], "UNIFORM")}
approx = RI(solver, mu0, params, verbosity = 0)
approx.setupApprox()
assert np.isclose(approx.normErr(mu)[0] / approx.normHF(mu)[0],
5.98695e-05, rtol = 1e-1)
def test_hermite():
mu = [5.05, 7.1]
mu0 = [5., 7.]
solver = matrixRandom()
sampler0 = QS([[4.9, 6.85], [5.1, 7.15]], "UNIFORM")
params = {"POD": True, "M": 3, "N": 3, "S": 25, "polybasis": "CHEBYSHEV",
"sampler": MS([[4.9, 6.85], [5.1, 7.15]],
points = sampler0.generatePoints(9))}
approx = RI(solver, mu0, params, verbosity = 0)
approx.setupApprox()
assert np.isclose(approx.normErr(mu)[0] / approx.normHF(mu)[0],
5.50053e-05, rtol = 5e-1)