Page MenuHomec4science
Files F61594270

trained_model_pivoted_rational_match.py
No OneTemporary
Actions

File Metadata

Created
Tue, May 7, 16:58

trained_model_pivoted_rational_match.py
View Options

# 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 <http://www.gnu.org/licenses/>.
#
import numpy as np
from copy import deepcopy as copy
from .trained_model_pivoted_rational_nomatch import (
TrainedModelPivotedRationalNoMatch)
from rrompy.utilities.base.types import (Tuple, 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 polynomialMatching
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
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, RROMPyAssert
from rrompy.sampling import sampleList
__all__ = ['TrainedModelPivotedRationalMatch']
class TrainedModelPivotedRationalMatch(TrainedModelPivotedRationalNoMatch):
"""
ROM approximant evaluation for pivoted approximants based on interpolation
of rational approximants (with some 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 = self.checkParameterListMarginal(mu)
if mu0 is None:
mu0 = self.data.mu0(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()
if pbM in ppb + rbpb:
if not extraPar:
approx.paramsMarginal["MMarginal"] = reduceDegreeN(
_MMarginalEff, len(musMCN), self.data.nparMarginal,
approx.paramsMarginal["polydegreetypeMarginal"])
MMEff = approx.paramsMarginal["MMarginal"]
while MMEff >= 0:
wellCond, msg = p.setupByInterpolation(musMCN, np.eye(nM),
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":
p.setupByInterpolation(musMCN, np.eye(nM), *interpPars)
else: # if pbM in sparsekinds:
p.setupByInterpolation(musMCN, np.eye(nM), pllims, extraPar,
*interpPars)
self.data.marginalInterp = p
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):
super().updateEffectiveSamples(exclude)
self.initializeFromRational(*args, **kwargs)
def initializeFromRational(self, matchingWeight:float, matchingKind:str,
HFEngine:HFEng, is_state:bool):
"""Initialize rational representation."""
RROMPyAssert(self.data.nparPivot, 1, "Number of pivot parameters")
Qs = [Q.coeffs for Q in self.data.Qs]
Ps = [P.coeffs for P in self.data.Ps]
N = len(Qs)
degQ = np.max([Q.shape[0] for Q in Qs])
degP = np.max([P.shape[0] for P in Ps])
for j in range(N):
if Qs[j].shape[0] < degQ:
Qs[j] = np.pad(Qs[j], (0, degQ - Qs[j].shape[0]), "constant")
if Ps[j].shape[0] < degP:
Ps[j] = np.pad(Ps[j], [(0, degP - Ps[j].shape[0]), (0, 0)],
"constant")
Qs, Ps = polynomialMatching(Qs, Ps, self.data.musMarginal.data,
matchingWeight, self.data.Psupp,
self.data.projMat, HFEngine, is_state, 0,
matchingKind)
for j in range(N):
if not isinstance(self.data.Qs[j], PI):
q = PI()
q.npar = self.data.Qs[j].npar
q.polybasis = self.data.Qs[j].polybasis
self.data.Qs[j] = q
if not isinstance(self.data.Ps[j], PI):
p = PI()
p.npar = self.data.Ps[j].npar
p.polybasis = self.data.Ps[j].polybasis
self.data.Ps[j] = p
self.data.Qs[j].coeffs, self.data.Ps[j].coeffs = Qs[j], Ps[j]
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 = mu(self.data.directionPivot)
muMC = self.centerNormalizeMarginal(
mu(self.data.directionMarginal))
mIvals = self.data.marginalInterp(muMC)
QV = self.getQVal(muP, mIvals = mIvals)
QVzero = np.where(QV == 0.)[0]
if len(QVzero) > 0:
QV[QVzero] = np.finfo(np.complex).eps / (1.
+ self.data.Qs[0].deg[0])
self.uApproxReduced = self.getPVal(muP, mIvals = mIvals) / QV
vbMng(self, "DEL", "Done evaluating approximant.", 12)
self.lastSolvedApproxReduced = mu
return self.uApproxReduced
def interpolateMarginalP(self, mu : paramList = [],
mIvals : Np2D = None) -> ListAny:
"""Obtain interpolated approximant numerator."""
mu = self.checkParameterListMarginal(mu)
vbMng(self, "INIT", "Interpolating marginal P at mu = {}.".format(mu),
95)
if self.data._collapsed:
outShape = self.data.Ps[0].coeffs.shape
else:
outShape = (self.data.Ps[0].coeffs.shape[0],
self.data.projMat.shape[1])
intMP = np.zeros((len(mu),) + outShape,
dtype = self.data.Ps[0].coeffs.dtype)
if mIvals is None:
muC = self.centerNormalizeMarginal(mu)
mIvals = self.data.marginalInterp(muC)
for P, Psupp, mI in zip(self.data.Ps, self.data.Psupp, mIvals):
iL, iR = Psupp, Psupp + P.shape[0]
for j, m in enumerate(mI): intMP[j, :, iL : iR] += m * P.coeffs
vbMng(self, "DEL", "Done interpolating marginal P.", 95)
return intMP
def interpolateMarginalQ(self, mu : paramList = [],
mIvals : Np2D = None) -> ListAny:
"""Obtain interpolated approximant denominator."""
mu = self.checkParameterListMarginal(mu)
vbMng(self, "INIT", "Interpolating marginal Q at mu = {}.".format(mu),
95)
intMQ = np.zeros((len(mu),) + self.data.Qs[0].coeffs.shape,
dtype = self.data.Qs[0].coeffs.dtype)
if mIvals is None:
muC = self.centerNormalizeMarginal(mu)
mIvals = self.data.marginalInterp(muC)
for Q, mI in zip(self.data.Qs, mIvals):
for j, m in enumerate(mI): intMQ[j] += m * Q.coeffs
vbMng(self, "DEL", "Done interpolating marginal Q.", 95)
return intMQ
def getPVal(self, mu : paramList = [],
mIvals : List[int] = None) -> sampList:
"""
Evaluate rational numerator at arbitrary parameter.
Args:
mu: Target parameter.
"""
RROMPyAssert(self.data.nparPivot, 1, "Number of pivot parameters")
if mIvals is None:
mu = self.checkParameterList(mu)
muP = self.centerNormalizePivot(mu(self.data.directionPivot))
muM = self.centerNormalizeMarginal(mu(self.data.directionMarginal))
mIvals = self.data.marginalInterp(muM)
else:
mu = self.checkParameterListPivot(mu)
muP = self.centerNormalizePivot(mu)
if self.data._collapsed:
outShape = self.data.Ps[0].shape
else:
outShape = (self.data.projMat.shape[1],)
p = np.zeros(outShape + (len(mu),),
dtype = self.data.Ps[0].coeffs.dtype)
for P, Psupp, mI in zip(self.data.Ps, self.data.Psupp, mIvals):
iL, iR = Psupp, Psupp + P.shape[0]
for j, m in enumerate(mI):
p[iL : iR, j] += m * P(muP[j])[:, 0]
return sampleList(p)
def getQVal(self, mu : paramList = [], der : List[int] = None,
scl : Np1D = None, mIvals : List[int] = 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")
if mIvals is None:
mu = self.checkParameterList(mu)
muP = self.centerNormalizePivot(mu(self.data.directionPivot))
muM = self.centerNormalizeMarginal(mu(self.data.directionMarginal))
mIvals = self.data.marginalInterp(muM)
else:
mu = self.checkParameterListPivot(mu)
muP = self.centerNormalizePivot(mu)
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]]
q = np.zeros(len(mu), dtype = np.complex)
for Q, mI in zip(self.data.Qs, mIvals):
for j, m in enumerate(mI): q[j] = q[j] + m * Q(muP[j], derP, sclP)
return q
def getPoles(self, marginalVals : ListAny = [fp]) -> paramList:
"""
Obtain approximant poles.
Returns:
Numpy complex vector of poles.
"""
RROMPyAssert(self.data.nparPivot, 1, "Number of pivot parameters")
mVals = list(marginalVals)
rDim = mVals.index(fp)
if rDim < len(mVals) - 1 and fp in mVals[rDim + 1 :]:
raise RROMPyException(("Exactly 1 'freepar' entry in "
"marginalVals must be provided."))
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]
Q = copy(self.data.Qs[0])
Q.coeffs = self.interpolateMarginalQ(mMarg)[0]
roots = self.data.scaleFactor[rDim] * Q.roots()
return self.mapParameterList(self.mapParameterList(self.data.mu0(rDim),
idx = [rDim])(0, 0)
+ roots, "B", [rDim])(0)
def getResidues(self, marginalVals : ListAny = [fp]) -> Tuple[paramList,
Np2D]:
"""
Obtain approximant residues.
Returns:
Numpy matrix with residues as columns.
"""
RROMPyAssert(self.data.nparPivot, 1, "Number of pivot parameters")
mVals = list(marginalVals)
rDim = mVals.index(fp)
if rDim < len(mVals) - 1 and fp in mVals[rDim + 1 :]:
raise RROMPyException(("Exactly 1 'freepar' entry in "
"marginalVals must be provided."))
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]
P, Q = copy(self.data.Ps[0]), copy(self.data.Qs[0])
P.coeffs = self.interpolateMarginalP(mMarg)[0]
Q.coeffs = self.interpolateMarginalQ(mMarg)[0]
res, pls, _ = rational2heaviside(P, Q)
if len(pls) == 0:
return pls, np.empty((0, 0), dtype = self.data.Ps[0].coeffs.dtype)
res = res[: len(pls), :].T
pls = self.mapParameterList(self.mapParameterList(self.data.mu0(rDim),
idx = [rDim])(0, 0)
+ self.data.scaleFactor[rDim] * pls, "B",
[rDim])(0)
if not self.data._collapsed: res = dot(self.data.projMat, res).T
return pls, res

Event Timeline

c4science · Help