Page MenuHomec4science
Files F63633943

trained_model_rational_mls.py
No OneTemporary
Actions

File Metadata

Created
Tue, May 21, 12:19

trained_model_rational_mls.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 numpy as np
from .trained_model_rational import TrainedModelRational
from rrompy.utilities.base.types import Np1D, paramVal, paramList, sampList
from rrompy.utilities.base import verbosityManager as vbMng
from rrompy.utilities.poly_fitting.moving_least_squares import mlsweights
from rrompy.utilities.poly_fitting.polynomial import (
PolynomialInterpolator as PI)
from rrompy.utilities.numerical import customPInv
from rrompy.utilities.numerical.degree import degreeTotalToFull
from rrompy.parameter import checkParameterList
from rrompy.sampling import emptySampleList
from rrompy.utilities.exception_manager import RROMPyException, RROMPyWarning
__all__ = ['TrainedModelRationalMLS']
class TrainedModelRationalMLS(TrainedModelRational):
"""
ROM approximant evaluation for rational moving least squares approximant.
Attributes:
Data: dictionary with all that can be pickled.
"""
def reset(self):
super().reset()
self.lastSetupMu = None
def assembleReducedModel(self, mu:paramVal):
if not (hasattr(self.data, "lastSetupMu")
and self.data.lastSetupMu == mu):
vbMng(self, "INIT", "Assembling reduced model for mu = {}."\
.format(mu), 17)
vbMng(self, "INIT", "Starting computation of denominator.", 35)
muC = self.centerNormalize(mu)
muSC = self.centerNormalize(self.data.mus)
wQ = mlsweights(muC, muSC, self.data.radialBasisDen,
directionalWeights = self.data.radialWeightsDen,
nNearestNeighbor = self.data.nNearestNeighborDen)
if self.data.N > self.data.M:
PQVan = self.data.QVan
else:
PQVan = self.data.PVan
VQAdjW = PQVan.conj().T * wQ
VQAdjWVQ = VQAdjW.dot(PQVan)
interpPseudoInverse, info = customPInv(VQAdjWVQ, full = True,
rcond = self.data.interpRcond)
interpPseudoInverse = interpPseudoInverse.dot(VQAdjW).dot(
self.data.QBlocks)
if info[0] < interpPseudoInverse.shape[-1]:
q = np.zeros(interpPseudoInverse.shape[-1], dtype = np.complex)
q[0] = 1.
else:
halfGram = interpPseudoInverse[self.data.domQIdxs]
if self.data.POD:
Rstack = halfGram.reshape(-1, halfGram.shape[-1])
vbMng(self, "INIT",
"Solving svd for square root of gramian matrix.", 67)
try:
_, s, eV = np.linalg.svd(Rstack, full_matrices = False)
except np.linalg.LinAlgError as e:
raise RROMPyException(e)
condN = s[0] / s[-1]
q = eV[-1, :].T.conj()
vbMng(self, "MAIN",
("Solved svd problem of size {} x {} with condition "
"number {:.4e}.").format(*Rstack.shape, condN), 55)
vbMng(self, "DEL", "Done solving svd.", 67)
else:
RRstack = np.tensordot(self.trainedModel.gramian, halfGram,
1).reshape(-1, halfGram.shape[-1])
RLstack = halfGram.reshape(-1, halfGram.shape[-1])
gram = RLstack.T.conj().dot(RRstack)
vbMng(self, "INIT",
"Solving eigenvalue problem for gramian matrix.", 67)
try:
ev, eV = np.linalg.eigh(gram)
except np.linalg.LinAlgError as e:
raise RROMPyException(e)
condN = ev[-1] / ev[0]
q = eV[:, 0]
vbMng(self, "MAIN",
("Solved eigenvalue problem of size {} with "
"condition number {:.4e}.").format(gram.shape[0],
condN), 55)
vbMng(self, "DEL", "Done solving eigenvalue problem.", 67)
self.data.Q = PI()
self.data.Q.npar = self.npar
self.data.Q.polybasis = self.data.polybasis
if self.data.polydegreetype == "TOTAL":
self.data.Q.coeffs = degreeTotalToFull(
(self.data.N + 1,) * self.npar,
self.npar, q)
else:
self.data.Q.coeffs = q.reshape((self.data.N + 1,) * self.npar)
vbMng(self, "DEL", "Done computing denominator.", 35)
vbMng(self, "INIT", "Starting computation of numerator.", 35)
self.data.P = PI()
self.data.P.npar = self.npar
self.data.P.polybasis = self.data.polybasis
wP = mlsweights(muC, muSC, self.data.radialBasis,
directionalWeights = self.data.radialWeights,
nNearestNeighbor = self.data.nNearestNeighbor)
VAdjW = self.data.PVan.conj().T * wP
VAdjWV = VAdjW.dot(self.data.PVan)
interpPPseudoInverse = customPInv(VAdjWV, self.data.interpRcond)
Pcoeffs = np.tensordot(interpPPseudoInverse.dot(VAdjW),
self.data.QBlocks.dot(q), ([1], [1]))
if self.data.polydegreetype == "TOTAL":
self.data.P.coeffs = degreeTotalToFull(
(self.data.M + 1,) * self.npar
+ (self.data.QBlocks.shape[0],),
self.npar, Pcoeffs)
else:
self.data.P.coeffs = Pcoeffs.reshape(
(self.data.M + 1,) * self.npar
+ (self.data.QBlocks.shape[0],))
vbMng(self, "DEL", "Done computing numerator.", 35)
vbMng(self, "DEL", "Done assembling reduced model.", 17)
self.data.lastSetupMu = mu
def getApproxReduced(self, mu : paramList = []) -> sampList:
"""
Evaluate reduced representation of approximant at arbitrary parameter.
Args:
mu: Target parameter.
"""
mu = checkParameterList(mu, self.data.npar)[0]
if (not hasattr(self, "lastSolvedApproxReduced")
or self.lastSolvedApproxReduced != mu):
vbMng(self, "INIT",
"Evaluating approximant at mu = {}.".format(mu), 12)
self.uApproxReduced = emptySampleList()
for i in range(len(mu)):
self.assembleReducedModel(mu[i])
vbMng(self, "INIT",
"Solving reduced model for mu = {}.".format(mu[i]), 15)
Qv = self.getQVal(mu[i])
if Qv == 0.:
RROMPyWarning(("Adjusting approximation to avoid division "
"by numerically zero denominator."))
Qv = np.finfo(np.complex).eps / (1. + self.data.Q.deg[0])
uAppR = self.getPVal(mu[i]) / Qv
if i == 0:
self.uApproxReduced.reset((uAppR.shape[0], len(mu)),
dtype = uAppR.dtype)
self.uApproxReduced[i] = uAppR
vbMng(self, "DEL", "Done solving reduced model.", 15)
vbMng(self, "DEL", "Done evaluating approximant.", 12)
self.lastSolvedApproxReduced = mu
return self.uApproxReduced
def getPoles(self, *args, mu : paramVal = None, **kwargs) -> Np1D:
"""
Obtain approximant poles.
Returns:
Numpy complex vector of poles.
"""
if mu is None: mu = self.data.mu0
self.assembleReducedModel(mu)
return super().getPoles(*args, **kwargs)
def getResidues(self, *args, mu : paramVal = None, **kwargs) -> Np1D:
"""
Obtain approximant residues.
Returns:
Numpy matrix with residues as columns.
"""
if mu is None: mu = self.data.mu0
self.assembleReducedModel(mu)
return super().getResidues(*args, **kwargs)

Event Timeline

c4science · Help