Homec4science
Diffusion RationalROMPy (develop) Browse

RROMPy/rrompy/reduction_methods/pivoted/trained_model/trained_model_pivoted_rational.py
22f2e4dda93ddevelop

Last Change
Raw File

trained_model_pivoted_rational.py
View Options

# Copyright (C) 2018 by the RROMPy authors
#
# This file is part of RROMPy.
#
# RROMPy is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# RROMPy is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with RROMPy. If not, see <http://www.gnu.org/licenses/>.
#
import warnings
import numpy as np
from scipy.sparse import csr_matrix, hstack, SparseEfficiencyWarning
from copy import deepcopy as copy
from .trained_model_pivoted_rational_nomatch import (
TrainedModelPivotedRationalNoMatch)
from rrompy.utilities.base.types import (Np1D, Np2D, List, ListAny, paramVal,
paramList, HFEng)
from rrompy.utilities.base import verbosityManager as vbMng
from rrompy.utilities.numerical.point_matching import (potential,
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 (heavisideUniformShape,
HeavisideInterpolator as HI)
from rrompy.utilities.poly_fitting.nearest_neighbor import (
NearestNeighborInterpolator as NNI)
from rrompy.utilities.poly_fitting.piecewise_linear import (
PiecewiseLinearInterpolator as PLI)
from rrompy.utilities.exception_manager import RROMPyException, RROMPyWarning
from rrompy.parameter import checkParameterList
__all__ = ['TrainedModelPivotedRational']
class TrainedModelPivotedRational(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 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 = checkParameterList(mu, self.data.nparMarginal)[0]
if mu0 is None: mu0 = self.data.mu0Marginal
rad = ((mu ** self.data.rescalingExpMarginal
- mu0 ** self.data.rescalingExpMarginal)
/ self.data.scaleFactorMarginal)
return rad
def setupMarginalInterp(self, approx, interpPars:ListAny,
MMAuto:bool, rDWM : Np1D = None,
extraPar = True):
vbMng(self, "INIT", "Starting computation of marginal interpolator.",
12)
musMCN = self.centerNormalizeMarginal(self.data.musMarginal)
nM = len(musMCN)
pbM = approx.polybasisMarginal
for ipts, pts in enumerate(self.data.suppEffPts):
if pbM in ppb:
p = PI()
elif pbM in rbpb:
p = RBI()
else: # #if pbM in sparsekinds + ["NEARESTNEIGHBOR"]:
if ipts > 0 or pbM == "NEARESTNEIGHBOR":
p = NNI()
else: #if ipts == 0 or pbM in sparsekinds:
pllims = [[-1.] * self.data.nparMarginal,
[1.] * self.data.nparMarginal]
p = PLI()
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 MMAuto:
if ipts > 0:
verb = approx.verbosity
approx.verbosity = 0
_musM = approx.musMarginal
approx.musMarginal = musMCNEff
approx._setMMarginalAuto()
if ipts > 0:
approx.musMarginal = _musM
approx.verbosity = verb
if ipts == 0:
_MMarginalEff = approx.paramsMarginal["MMarginal"]
if not MMAuto:
approx.paramsMarginal["MMarginal"] = _MMarginalEff
MM = reduceDegreeN(approx.paramsMarginal["MMarginal"],
len(musMCNEff), self.data.nparMarginal,
approx.paramsMarginal["polydegreetypeMarginal"])
if MM < approx.paramsMarginal["MMarginal"]:
if ipts == 0 and not extraPar:
RROMPyWarning(
("MMarginal too large compared to SMarginal. Reducing "
"MMarginal by {}").format(
approx.paramsMarginal["MMarginal"] - MM))
approx.paramsMarginal["MMarginal"] = MM
MMEff = approx.paramsMarginal["MMarginal"]
while MMEff >= 0:
pParRest = copy(interpPars)
if pbM in rbpb:
pParRest = [rDWM] + pParRest
wellCond, msg = p.setupByInterpolation(musMCNEff,
valsEff, MMEff,
pbM, *pParRest)
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."))
elif ipts > 0 or pbM == "NEARESTNEIGHBOR":
nNeigh = interpPars[0] if ipts == 0 else 1
p.setupByInterpolation(musMCNEff, valsEff, nNeigh, rDWM)
else: #if ipts == 0 and pbM in sparsekinds:
wellCond, msg = p.setupByInterpolation(musMCNEff, valsEff,
pllims, extraPar[pts],
pbM, *interpPars)
vbMng(self, "MAIN", msg, 30)
if not wellCond:
vbMng(self, "MAIN",
"Warning: polyfit is poorly conditioned.", 35)
if ipts == 0:
self.data.marginalInterp = copy(p)
self.data.polesEff = [copy(hi.poles) for hi in self.data.HIs]
self.data.coeffsEff = []
for hi, sup in zip(self.data.HIs, self.data.Psupp):
cEff = hi.coeffs
if (self.data._collapsed
or self.data.projMat.shape[1] == cEff.shape[1]):
cEff = copy(cEff)
else:
cEff = hstack((csr_matrix((len(cEff), sup)),
csr_matrix(cEff),
csr_matrix((len(cEff),
self.data.projMat.shape[1]
- sup - cEff.shape[1]))),
"csr")
self.data.coeffsEff += [cEff]
else:
ptsBad = [i for i in range(nM) if i not in pts]
musMCNBad = musMCN[ptsBad]
idxBad = np.where(self.data.suppEffIdx == ipts)[0]
valMuMBad = p(musMCNBad)
pls = np.empty((len(ptsBad), len(idxBad)),
dtype = self.data.polesEff[0].dtype)
warnings.simplefilter('ignore', SparseEfficiencyWarning)
if (self.data._collapsed
or self.data.projMat.shape[1] == cEff.shape[1]):
cfs = [0.] * len(ptsBad)
else:
cfs = [None] * len(ptsBad)
for jj in range(len(ptsBad)):
cfs[jj] = csr_matrix(
(len(idxBad), self.data.projMat.shape[1]),
dtype = self.data.coeffsEff[0].dtype)
for ij, j in enumerate(pts):
pls += (np.expand_dims(valMuMBad[ij], -1)
* self.data.polesEff[j][idxBad])
for jj, val in enumerate(valMuMBad[ij]):
cfs[jj] = (cfs[jj]
+ val * self.data.coeffsEff[j][idxBad])
for ij, j in enumerate(ptsBad):
self.data.polesEff[j][idxBad] = pls[ij]
self.data.coeffsEff[j][idxBad] = cfs[ij]
if pbM in ppb + rbpb:
approx.paramsMarginal["MMarginal"] = _MMarginalEff
vbMng(self, "DEL", "Done computing marginal interpolator.", 12)
def initializeFromLists(self, poles:ListAny, coeffs:ListAny, supps:ListAny,
basis:str, matchingWeight:float, matchingMode:str,
HFEngine:HFEng, is_state:bool):
"""Initialize Heaviside representation."""
poles, coeffs = rationalFunctionMatching(
*heavisideUniformShape(poles, coeffs),
self.data.musMarginal.data, matchingWeight,
matchingMode, supps, self.data.projMat,
HFEngine, is_state)
super().initializeFromLists(poles, coeffs, supps, basis)
self.data.suppEffPts = [np.arange(len(self.data.HIs))]
self.data.suppEffIdx = np.zeros(len(poles[0]), dtype = int)
def recompressByCutOff(self, tol:float, shared:float,
foci:List[np.complex], ground:float) -> str:
N = len(self.data.HIs[0].poles)
M = len(self.data.HIs)
mu0 = np.mean(foci)
goodLocPoles = np.array([potential(hi.poles, foci) - ground
<= tol * ground 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.append([N],
np.arange(N + 1, len(self.data.HIs[0].coeffs))))
for hi in self.data.HIs:
polyCorrection = np.zeros_like(hi.coeffs[0, :])
for j in removePole:
if not np.isinf(hi.poles[j]):
polyCorrection += hi.coeffs[j, :] / (mu0 - hi.poles[j])
if len(hi.coeffs) == N:
hi.coeffs = np.vstack((hi.coeffs, polyCorrection))
else:
hi.coeffs[N, :] += polyCorrection
hi.poles = hi.poles[keepPole]
hi.coeffs = hi.coeffs[keepCoeff, :]
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 interpolateMarginalInterpolator(self, mu : paramList = []) -> ListAny:
"""Obtain interpolated approximant interpolator."""
mu = checkParameterList(mu, self.data.nparMarginal)[0]
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 = checkParameterList(mu, self.data.nparMarginal)[0]
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 = checkParameterList(mu, self.data.nparMarginal)[0]
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
c4science · Help