diff --git a/rrompy/reduction_methods/pivoted/trained_model/trained_model_pivoted_rational.py b/rrompy/reduction_methods/pivoted/trained_model/trained_model_pivoted_rational.py index e6086ed..4c503cb 100644 --- a/rrompy/reduction_methods/pivoted/trained_model/trained_model_pivoted_rational.py +++ b/rrompy/reduction_methods/pivoted/trained_model/trained_model_pivoted_rational.py @@ -1,335 +1,341 @@ # 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 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 (pointMatching, chordalMetricAdjusted, potential) 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, rational2heaviside, 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, RROMPyAssert, 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 sup is None: + 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) - cfs = [None] * len(ptsBad) warnings.simplefilter('ignore', SparseEfficiencyWarning) - for jj in range(len(ptsBad)): - cfs[jj] = csr_matrix( + 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] += val * self.data.coeffsEff[j][idxBad] + 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, HFEngine:HFEng, matchingWeight : float = 1., is_state : bool = True): """Initialize Heaviside representation.""" musM = self.data.musMarginal margAbsDist = np.sum(np.abs(np.repeat(musM.data, len(musM), 0) - np.tile(musM.data, [len(musM), 1]) ), axis = 1).reshape(len(musM), len(musM)) explored = [0] unexplored = list(range(1, len(musM))) poles, coeffs = heavisideUniformShape(poles, coeffs) N = len(poles[0]) for _ in range(1, len(musM)): minIdx = np.argmin(np.concatenate([margAbsDist[ex, unexplored] \ for ex in explored])) minIex = explored[minIdx // len(unexplored)] minIunex = unexplored[minIdx % len(unexplored)] resex = coeffs[minIex][: N].T resunex = coeffs[minIunex][: N].T if matchingWeight != 0 and not self.data._collapsed: suppex, suppunex = supps[minIex], supps[minIunex] resex = self.data.projMat[:, suppex : suppex + len(resex)].dot(resex) resunex = self.data.projMat[:, suppunex : suppunex + len(resunex)].dot(resunex) dist = chordalMetricAdjusted(poles[minIex], poles[minIunex], matchingWeight, resex, resunex, HFEngine, is_state) reordering = pointMatching(dist)[1] poles[minIunex] = poles[minIunex][reordering] coeffs[minIunex][: N] = coeffs[minIunex][reordering] explored += [minIunex] unexplored.remove(minIunex) super().initializeFromLists(poles, coeffs, supps, basis) self.data.suppEffPts = [np.arange(len(self.data.HIs))] self.data.suppEffIdx = np.zeros(N, dtype = int) def initializeFromRational(self, HFEngine:HFEng, matchingWeight : float = 1., 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, self.data.Psupp, basis, HFEngine, matchingWeight, is_state) 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] muC = self.centerNormalizeMarginal(mu) mIvals = self.data.marginalInterp(muC) poless = self.interpolateMarginalPoles(mu, mIvals) coeffss = self.interpolateMarginalCoeffs(mu, mIvals) his = [None] * len(mu) for j in range(len(mu)): his[j] = HI() his[j].poles = poless[j] his[j].coeffs = coeffss[j] his[j].npar = 1 his[j].polybasis = self.data.HIs[0].polybasis 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: mIvals = self.data.marginalInterp(self.centerNormalizeMarginal(mu)) 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: mIvals = self.data.marginalInterp(self.centerNormalizeMarginal(mu)) 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 diff --git a/rrompy/reduction_methods/pivoted/trained_model/trained_model_pivoted_rational_nomatch.py b/rrompy/reduction_methods/pivoted/trained_model/trained_model_pivoted_rational_nomatch.py index 62f6053..4013279 100644 --- a/rrompy/reduction_methods/pivoted/trained_model/trained_model_pivoted_rational_nomatch.py +++ b/rrompy/reduction_methods/pivoted/trained_model/trained_model_pivoted_rational_nomatch.py @@ -1,341 +1,342 @@ # 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 copy import deepcopy as copy from scipy.special import factorial as fact from itertools import combinations from rrompy.reduction_methods.standard.trained_model.trained_model_rational \ import TrainedModelRational from rrompy.utilities.base.types import (Np1D, Np2D, List, ListAny, paramVal, paramList, sampList) from rrompy.utilities.base import verbosityManager as vbMng, freepar as fp from rrompy.utilities.numerical import dot from rrompy.utilities.numerical.compress_matrix import compressMatrix from rrompy.utilities.numerical.point_matching import potential from rrompy.utilities.poly_fitting.heaviside import (rational2heaviside, HeavisideInterpolator as HI) from rrompy.utilities.poly_fitting.nearest_neighbor import ( NearestNeighborInterpolator as NNI) from rrompy.utilities.exception_manager import (RROMPyException, RROMPyAssert, RROMPyWarning) from rrompy.parameter import checkParameterList from rrompy.sampling import emptySampleList __all__ = ['TrainedModelPivotedRationalNoMatch'] class TrainedModelPivotedRationalNoMatch(TrainedModelRational): """ ROM approximant evaluation for pivoted approximants based on interpolation of rational approximants (without pole matching). 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 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) for obj, suppj in zip(self.data.HIs, self.data.Psupp): obj.postmultiplyTensorize(RMat.T[suppj : suppj + obj.shape[0]]) if hasattr(self, "_HIsExcl"): for obj, suppj in zip(self._HIsExcl, self.data.Psupp): obj.postmultiplyTensorize(RMat.T[suppj : suppj + obj.shape[0]]) if hasattr(self.data, "Ps"): for obj, suppj in zip(self.data.Ps, self.data.Psupp): obj.postmultiplyTensorize(RMat.T[suppj : suppj + obj.shape[0]]) if hasattr(self, "_PsExcl"): for obj, suppj in zip(self._PsExcl, self._PsuppExcl): obj.postmultiplyTensorize(RMat.T[suppj : suppj + obj.shape[0]]) if hasattr(self.data, "coeffsEff"): for j in range(len(self.data.coeffsEff)): self.data.coeffsEff[j] = dot(self.data.coeffsEff[j], RMat.T) if hasattr(self, "_HIsExcl") or hasattr(self, "_PsExcl"): self._PsuppExcl = [0] * len(self._PsuppExcl) self.data.Psupp = [0] * len(self.data.Psupp) super(TrainedModelRational, self).compress(collapse, tol) def centerNormalizePivot(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.mu0Pivot. Returns: Normalized parameter. """ mu = checkParameterList(mu, self.data.nparPivot)[0] if mu0 is None: mu0 = self.data.mu0Pivot rad = ((mu ** self.data.rescalingExpPivot - mu0 ** self.data.rescalingExpPivot) / self.data.scaleFactorPivot) return rad def setupMarginalInterp(self, rDWM : Np1D = None): self.data.NN = NNI() self.data.NN.setupByInterpolation(self.data.musMarginal, np.arange(len(self.data.musMarginal)), 1, rDWM) def updateEffectiveSamples(self, exclude:List[int], *args, **kwargs): if hasattr(self, "_idxExcl"): for j, excl in enumerate(self._idxExcl): self.data.musMarginal.insert(self._musMExcl[j], excl) self.data.HIs.insert(excl, self._HIsExcl[j]) self.data.Ps.insert(excl, self._PsExcl[j]) self.data.Qs.insert(excl, self._QsExcl[j]) self.data.Psupp.insert(excl, self._PsuppExcl[j]) self._idxExcl, self._musMExcl = list(np.sort(exclude)), [] self._HIsExcl, self._PsExcl, self._QsExcl = [], [], [] self._PsuppExcl = [] for excl in self._idxExcl[::-1]: self._musMExcl = [self.data.musMarginal[excl]] + self._musMExcl self.data.musMarginal.pop(excl) self._HIsExcl = [self.data.HIs.pop(excl)] + self._HIsExcl self._PsExcl = [self.data.Ps.pop(excl)] + self._PsExcl self._QsExcl = [self.data.Qs.pop(excl)] + self._QsExcl self._PsuppExcl = [self.data.Psupp.pop(excl)] + self._PsuppExcl 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 initializeFromLists(self, poles:ListAny, coeffs:ListAny, supps:ListAny, basis:str): """Initialize Heaviside representation.""" self.data.HIs = [] for pls, cfs in zip(poles, coeffs): hsi = HI() hsi.poles = pls hsi.coeffs = cfs hsi.npar = 1 hsi.polybasis = basis self.data.HIs += [hsi] def initializeFromRational(self): """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) def recompressByCutOff(self, tol:float, foci:List[np.complex], ground:float) -> str: mu0 = np.mean(foci) gLocPoles = [potential(hi.poles, foci) - ground <= tol * ground for hi in self.data.HIs] nRemPole = np.sum([np.sum(np.logical_not(gLPi)) for gLPi in gLocPoles]) if nRemPole == 0: return " No poles erased." for hi, gLocPolesi in zip(self.data.HIs, gLocPoles): N = len(hi.poles) polyCorrection = np.zeros_like(hi.coeffs[0, :]) for j, goodj in enumerate(gLocPolesi): if not goodj and 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[gLocPolesi] gLocCoeffi = np.append(gLocPolesi, np.ones(len(hi.coeffs) - N, dtype = bool)) hi.coeffs = hi.coeffs[gLocCoeffi, :] return " Erased {} pole{}.".format(nRemPole, "s" * (nRemPole != 1)) def interpolateMarginalInterpolator(self, mu : paramList = []) -> ListAny: """Obtain interpolated approximant interpolator.""" mu = checkParameterList(mu, self.data.nparMarginal)[0] vbMng(self, "INIT", "Finding nearest neighbor to mu = {}.".format(mu), 95) intM = np.array(self.data.NN(mu), dtype = int) vbMng(self, "DEL", "Done finding nearest neighbor.", 95) his = [None] * len(mu) for j in range(len(mu)): i = intM[j] his[j] = HI() his[j].poles = copy(self.data.HIs[i].poles) his[j].coeffs = copy(self.data.HIs[i].coeffs) his[j].npar = 1 his[j].polybasis = self.data.HIs[0].polybasis - his[j].pad(self.data.Psupp[i], - self.data.projMat.shape[1] - self.data.Psupp[i] - - his[j].shape[0]) + if not self.data._collapsed: + his[j].pad(self.data.Psupp[i], + self.data.projMat.shape[1] - self.data.Psupp[i] + - his[j].shape[0]) return his def interpolateMarginalPoles(self, mu : paramList = []) -> ListAny: """Obtain interpolated approximant poles.""" interps = self.interpolateMarginalInterpolator(mu) return [interp.poles for interp in interps] def interpolateMarginalCoeffs(self, mu : paramList = []) -> ListAny: """Obtain interpolated approximant poles.""" interps = self.interpolateMarginalInterpolator(mu) return [interp.coeffs for interp in interps] 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 = 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() muP = self.centerNormalizePivot(mu.data[:, self.data.directionPivot]) muM = mu.data[:, self.data.directionMarginal] his = self.interpolateMarginalInterpolator(muM) for i, (mP, hi) in enumerate(zip(muP, his)): uAppR = hi(mP) if i == 0: self.uApproxReduced.reset((len(uAppR), len(mu)), dtype = uAppR.dtype) self.uApproxReduced[i] = uAppR vbMng(self, "DEL", "Done evaluating approximant.", 12) self.lastSolvedApproxReduced = mu return self.uApproxReduced 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 = checkParameterList(mu, self.data.npar)[0] p = emptySampleList() muP = self.centerNormalizePivot(mu.data[:, self.data.directionPivot]) muM = mu.data[:, 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 = checkParameterList(mu, self.data.npar)[0] muP = self.centerNormalizePivot(mu.data[:, self.data.directionPivot]) muM = mu.data[:, 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 = muP[0] - pl for terms in combinations(np.arange(N), N - derP): derVal[i] += np.prod(plDist[list(terms)], axis = 1) 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 hasattr(mVals, "__len__"): mVals = [mVals] mVals = list(mVals) else: mVals = [fp] 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.")) if rDim != self.data.directionPivot[0]: raise RROMPyException(("'freepar' entry in marginalVals must " "coincide with pivot direction.")) mVals[rDim] = self.data.mu0(rDim) mMarg = [mVals[j] for j in range(len(mVals)) if j != rDim] roots = self.interpolateMarginalPoles(mMarg)[0] return np.power(self.data.mu0(rDim) ** self.data.rescalingExp[rDim] + self.data.scaleFactor[rDim] * roots, 1. / self.data.rescalingExp[rDim]) 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 hasattr(mVals, "__len__"): 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), :] if not self.data._collapsed: res = self.data.projMat.dot(res.T).T return pls, res