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 8fb1fc4..63754a6 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,318 +1,318 @@ # 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 (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 (sparsekinds, PiecewiseLinearInterpolator as PLI) from rrompy.utilities.exception_manager import RROMPyException __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 = self.checkParameterListMarginal(mu) if mu0 is None: mu0 = self.checkParameterListMarginal( self.data.mu0(0, 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, nP = len(musMCN), len(self.data.HIs[0].poles) pbM = 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() for ipts, pts in enumerate(self.data.suppEffPts): 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 extraPar: if ipts > 0: verb = approx.verbosity approx.verbosity = 0 _musM = approx.musMarginal approx.musMarginal = musMCNEff approx._setMMarginalAuto() approx.musMarginal = _musM approx.verbosity = verb else: approx.paramsMarginal["MMarginal"] = reduceDegreeN( _MMarginalEff, len(musMCNEff), self.data.nparMarginal, approx.paramsMarginal["polydegreetypeMarginal"]) MMEff = approx.paramsMarginal["MMarginal"] while MMEff >= 0: wellCond, msg = p.setupByInterpolation(musMCNEff, valsEff, 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": if ipts > 0: interpPars[0] = 1 p.setupByInterpolation(musMCNEff, valsEff, *interpPars) elif ipts == 0: # and pbM in sparsekinds: wellCond, msg = p.setupByInterpolation(musMCNEff, valsEff, pllims, extraPar[pts], *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.coeffsEff, self.data.polesEff = [], [] 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: supC = self.data.projMat.shape[1] - sup - cEff.shape[1] cEff = hstack((csr_matrix((len(cEff), sup)), csr_matrix(cEff), csr_matrix((len(cEff), supC))), "csr") self.data.coeffsEff += [cEff] self.data.polesEff += [copy(hi.poles)] else: ptsBad = [i for i in range(nM) if i not in pts] idxBad = np.where(self.data.suppEffIdx == ipts)[0] if pbM in sparsekinds: for ij, j in enumerate(ptsBad): for jb in idxBad: poleBad = self.data.polesEff[j][jb] if not np.isinf(poleBad): self.data.coeffsEff[j][nP] -= ( self.data.coeffsEff[j][jb] / poleBad) bdsL = np.append([0], self.data.coeffsEff[j].indptr[idxBad + 1]) bdsR = np.append(self.data.coeffsEff[j].indptr[idxBad], self.data.coeffsEff[j].indptr[-1]) self.data.coeffsEff[j].data = np.concatenate([ self.data.coeffsEff[j].data[l : r] for l, r in zip(bdsL, bdsR)]) self.data.coeffsEff[j].indices = np.concatenate([ self.data.coeffsEff[j].indices[l : r] for l, r in zip(bdsL, bdsR)]) for ijb, jb in enumerate(idxBad): self.data.coeffsEff[j].indptr[jb + 1 :] -= ( bdsL[ijb + 1] - bdsR[ijb]) self.data.polesEff[j][idxBad] = np.inf else: warnings.simplefilter('ignore', SparseEfficiencyWarning) if (self.data._collapsed or self.data.projMat.shape[1] == cEff.shape[1]): cfBase = np.zeros((len(idxBad), cEff.shape[1]), dtype = cEff.dtype) else: cfBase = csr_matrix((len(idxBad), self.data.projMat.shape[1]), dtype = cEff.dtype) valMuMBad = p(musMCN[ptsBad]) for ijb, jb in enumerate(ptsBad): self.data.coeffsEff[jb][idxBad] = copy(cfBase) self.data.polesEff[jb][idxBad] = 0. for ij, j in enumerate(pts): val = valMuMBad[ij][ijb] if not np.isclose(val, 0.): self.data.coeffsEff[jb][idxBad] += (val * self.data.coeffsEff[j][idxBad]) self.data.polesEff[jb][idxBad] += (val * self.data.polesEff[j][idxBad]) warnings.filters.pop(0) 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) goodLocPoles = np.array([np.logical_and( np.logical_not(np.isinf(hi.poles)), 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.arange(N, len(self.data.HIs[0].coeffs))) for hi in self.data.HIs: for j in removePole: if not np.isinf(hi.poles[j]): hi.coeffs[N, :] -= hi.coeffs[j, :] / hi.poles[j] 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 = self.checkParameterListMarginal(mu) 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 = self.checkParameterListMarginal(mu) 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): mIBad, pEffBad = np.logical_not(np.isclose(mI, 0.)), np.isinf(pEff) pEffGood = np.logical_not(pEffBad) if np.sum(pEffGood) > 0: intMPoles[:, pEffGood] += ( np.expand_dims(mI, - 1) * pEff[pEffGood]) if np.sum(mIBad) * np.sum(pEffBad) > 0: - intMPoles[mIBad, pEffBad] = np.inf + intMPoles[np.ix_(mIBad, pEffBad)] = np.inf vbMng(self, "DEL", "Done interpolating marginal poles.", 95) return intMPoles def interpolateMarginalCoeffs(self, mu : paramList = [], mIvals : Np2D = None) -> ListAny: """Obtain interpolated approximant coefficients.""" mu = self.checkParameterListMarginal(mu) 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