diff --git a/rrompy/reduction_methods/trained_model/trained_model_data.py b/rrompy/reduction_methods/trained_model/trained_model_data.py index 766f236..00dd780 100644 --- a/rrompy/reduction_methods/trained_model/trained_model_data.py +++ b/rrompy/reduction_methods/trained_model/trained_model_data.py @@ -1,37 +1,36 @@ # 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 . # from copy import deepcopy as copy from rrompy.utilities.base.types import Np2D, List, paramVal from rrompy.utilities.exception_manager import RROMPyAssert __all__ = ['TrainedModelData'] class TrainedModelData: """ROM approximant evaluation data (must be pickle-able).""" - def __init__(self, name:str, mu0:paramVal, projMat:Np2D, + def __init__(self, mu0:paramVal, projMat:Np2D, scaleFactor : List[float] = [1.], rescalingExp : List[float] = [1.]): self.npar = len(rescalingExp) RROMPyAssert(mu0.shape[1], self.npar, "Number of parameters") - self.name = name self.mu0 = mu0 self.projMat = copy(projMat) self.scaleFactor = scaleFactor self.rescalingExp = rescalingExp diff --git a/rrompy/reduction_methods/trained_model/trained_model_pivoted_data.py b/rrompy/reduction_methods/trained_model/trained_model_pivoted_data.py index 79b8803..27b7a58 100644 --- a/rrompy/reduction_methods/trained_model/trained_model_pivoted_data.py +++ b/rrompy/reduction_methods/trained_model/trained_model_pivoted_data.py @@ -1,72 +1,72 @@ # 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 . # from .trained_model_data import TrainedModelData from rrompy.utilities.base.types import Np2D, List, ListAny, paramVal from rrompy.parameter import checkParameterList __all__ = ['TrainedModelPivotedData'] class TrainedModelPivotedData(TrainedModelData): """ROM approximant evaluation data (must be pickle-able).""" - def __init__(self, name:str, mu0:paramVal, projMat:Np2D, + def __init__(self, mu0:paramVal, projMat:Np2D, scaleFactor : ListAny = [1.], rescalingExp : List[float] = [1.], directionPivot : ListAny = [0]): - super().__init__(name, mu0, projMat, scaleFactor, rescalingExp) + super().__init__(mu0, projMat, scaleFactor, rescalingExp) self.directionPivot = directionPivot @property def directionMarginal(self): return tuple([x for x in range(self.npar) \ if x not in self.directionPivot]) @property def mu0Pivot(self): return checkParameterList(self.mu0(0, self.directionPivot), self.nparPivot)[0] @property def mu0Marginal(self): return checkParameterList(self.mu0(0, self.directionMarginal), self.nparMarginal)[0] @property def nparPivot(self): return len(self.directionPivot) @property def nparMarginal(self): return self.npar - self.nparPivot @property def rescalingExpPivot(self): return [self.rescalingExp[x] for x in self.directionPivot] @property def rescalingExpMarginal(self): return [self.rescalingExp[x] for x in self.directionMarginal] @property def scaleFactorPivot(self): return [self.scaleFactor[x] for x in self.directionPivot] @property def scaleFactorMarginal(self): return [self.scaleFactor[x] for x in self.directionMarginal] diff --git a/rrompy/reduction_methods/trained_model/trained_model_pivoted_rational.py b/rrompy/reduction_methods/trained_model/trained_model_pivoted_rational.py index 98f8117..0e22cbe 100644 --- a/rrompy/reduction_methods/trained_model/trained_model_pivoted_rational.py +++ b/rrompy/reduction_methods/trained_model/trained_model_pivoted_rational.py @@ -1,107 +1,107 @@ # 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 . # from .trained_model_pivoted_general import TrainedModelPivotedGeneral from .trained_model_rational import TrainedModelRational from rrompy.utilities.base.types import Np1D, List, paramList, sampList, HFEng from rrompy.utilities.base import verbosityManager as vbMng from rrompy.utilities.poly_fitting.heaviside import rational2heaviside from rrompy.utilities.exception_manager import RROMPyAssert from rrompy.parameter import checkParameterList from rrompy.sampling import sampleList __all__ = ['TrainedModelPivotedRational'] class TrainedModelPivotedRational(TrainedModelPivotedGeneral, TrainedModelRational): """ ROM approximant evaluation for pivoted Rational approximant (with pole matching). Attributes: Data: dictionary with all that can be pickled. """ def initializeFromRational(self, HFEngine:HFEng, matchingWeight : float = 1., POD : bool = True, is_state : bool = True): """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, basis, HFEngine, matchingWeight, POD, is_state) - self.data._temporary = False 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") - if not self.data._temporary: return super().getPVal(mu) + if not hasattr(self.data, "_temporary"): return super().getPVal(mu) mu = checkParameterList(mu, self.data.npar)[0] muP = checkParameterList(mu.data[:, self.data.directionPivot], self.data.nparPivot)[0] muM = checkParameterList(mu.data[:, self.data.directionMarginal], self.data.nparMarginal)[0] muP = checkParameterList(muP, self.data.nparPivot)[0] vbMng(self, "INIT", "Evaluating numerator at mu = {}.".format(muP), 17) muPC = self.centerNormalizePivot(muP) pP = [sampleList(P(muPC)) for P in self.data.Ps] vbMng(self, "DEL", "Done evaluating numerator.", 17) return self.interpolateMarginal(muM, pP) 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") - if not self.data._temporary: return super().getQVal(mu, der, scl) + if not hasattr(self.data, "_temporary"): + return super().getQVal(mu, der, scl) mu = checkParameterList(mu, self.data.npar)[0] muP = checkParameterList(mu.data[:, self.data.directionPivot], self.data.nparPivot)[0] muM = checkParameterList(mu.data[:, self.data.directionMarginal], self.data.nparMarginal)[0] if der is None: derP, derM = None, None else: derP = [der[x] for x in self.data.directionPivot] derM = [der[x] for x in self.data.directionMarginal] if scl is None: sclP, sclM = None, None else: sclP = [scl[x] for x in self.data.directionPivot] sclM = [scl[x] for x in self.data.directionMarginal] muP = checkParameterList(muP, self.data.nparPivot)[0] vbMng(self, "INIT", "Evaluating denominator at mu = {}.".format(muP), 17) muPC = self.centerNormalizePivot(muP) qP = [Q(muPC, derP, sclP).reshape(1, -1) for Q in self.data.Qs] vbMng(self, "DEL", "Done evaluating denominator.", 17) return self.interpolateMarginal(muM, qP, derM, sclM) diff --git a/rrompy/reduction_methods/trained_model/trained_model_rational_mls.py b/rrompy/reduction_methods/trained_model/trained_model_rational_mls.py index f012a55..180b232 100644 --- a/rrompy/reduction_methods/trained_model/trained_model_rational_mls.py +++ b/rrompy/reduction_methods/trained_model/trained_model_rational_mls.py @@ -1,174 +1,175 @@ # 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 . # 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, degreeTotalToFull from rrompy.parameter import checkParameterList from rrompy.sampling import emptySampleList __all__ = ['TrainedModelRationalMLS'] class TrainedModelRationalMLS(TrainedModelRational): """ ROM approximant evaluation for rational moving least squares approximant. Attributes: Data: dictionary with all that can be pickled. """ def assembleReducedModel(self, mu:paramVal): - if not hasattr(self, "lastSetupMu") or self.lastSetupMu != mu: + 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) _, s, eV = np.linalg.svd(Rstack, full_matrices = False) 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) ev, eV = np.linalg.eigh(gram) 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.lastSetupMu = mu + 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) uAppR = self.getPVal(mu[i]) / self.getQVal(mu[i]) if i == 0: #self.data.P.shape[-1], len(mu) self.uApproxReduced.reset((len(uAppR), 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) diff --git a/rrompy/reduction_methods/trained_model/trained_model_reduced_basis.py b/rrompy/reduction_methods/trained_model/trained_model_reduced_basis.py index 33e8364..975987e 100644 --- a/rrompy/reduction_methods/trained_model/trained_model_reduced_basis.py +++ b/rrompy/reduction_methods/trained_model/trained_model_reduced_basis.py @@ -1,115 +1,116 @@ # 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 . # import numpy as np from .trained_model import TrainedModel from rrompy.utilities.base.types import (Np1D, ListAny, paramVal, paramList, sampList) from rrompy.utilities.base import verbosityManager as vbMng, freepar as fp from rrompy.utilities.numerical import (eigvalsNonlinearDense, marginalizePolyList) from rrompy.utilities.expression import expressionEvaluator from rrompy.utilities.exception_manager import RROMPyException, RROMPyWarning from rrompy.parameter import checkParameter, checkParameterList from rrompy.sampling import emptySampleList __all__ = ['TrainedModelReducedBasis'] class TrainedModelReducedBasis(TrainedModel): """ ROM approximant evaluation for RB approximant. Attributes: Data: dictionary with all that can be pickled. """ def assembleReducedModel(self, mu:paramVal): mu = checkParameter(mu, self.data.npar) - if not hasattr(self, "lastSetupMu") or self.lastSetupMu != mu: + if not (hasattr(self.data, "lastSetupMu") + and self.data.lastSetupMu == mu): vbMng(self, "INIT", "Assembling reduced model for mu = {}."\ .format(mu), 17) muEff = mu ** self.data.rescalingExp self.data.ARBmu, self.data.bRBmu = 0., 0. for thA, ARB in zip(self.data.thAs, self.data.ARBs): self.data.ARBmu = (expressionEvaluator(thA[0], muEff) * ARB + self.data.ARBmu) for thb, bRB in zip(self.data.thbs, self.data.bRBs): self.data.bRBmu = (expressionEvaluator(thb[0], muEff) * bRB + self.data.bRBmu) vbMng(self, "DEL", "Done assembling reduced model.", 17) - self.lastSetupMu = mu + 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", "Computing RB solution 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) uAppR = np.linalg.solve(self.data.ARBmu, self.data.bRBmu) if i == 0: #self.data.ARBs[0].shape[-1], len(mu) self.uApproxReduced.reset((len(uAppR), len(mu)), dtype = uAppR.dtype) self.uApproxReduced[i] = uAppR vbMng(self, "DEL", "Done solving reduced model.", 15) vbMng(self, "DEL", "Done computing RB solution.", 12) self.lastSolvedApproxReduced = mu return self.uApproxReduced def getPoles(self, marginalVals : ListAny = [fp], jSupp : int = 1, **kwargs) -> Np1D: """ Obtain approximant poles. Returns: Numpy complex vector of poles. """ if not self.data.affinePoly: RROMPyWarning(("Unable to compute approximate poles due " "to parametric dependence (detected non-" "polynomial). Change HFEngine.affinePoly to True " "if necessary.")) return if not hasattr(marginalVals, "__len__"): marginalVals = [marginalVals] mVals = list(marginalVals) try: rDim = mVals.index(fp) if rDim < len(mVals) - 1 and fp in mVals[rDim + 1 :]: raise except: raise RROMPyException(("Exactly 1 'freepar' entry in " "marginalVals must be provided.")) ARBs = self.data.ARBs if self.data.npar > 1: mVals[rDim] = self.data.mu0(rDim) mVals = checkParameter(mVals).data.flatten() mVals[rDim] = fp ARBs = marginalizePolyList(ARBs, mVals, "auto") ev = eigvalsNonlinearDense(ARBs, jSupp = jSupp, **kwargs) return np.sort(np.power(ev, 1. / self.data.rescalingExp[rDim]))