rational_interpolant_pivoted.py
6643b38fd8aemaster

rational_interpolant_pivoted.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/>.
			#

			from copy import deepcopy as copy
			import numpy as np
			from rrompy.reduction_methods.base import checkRobustTolerance
			from .generic_pivoted_approximant import GenericPivotedApproximant
			from rrompy.reduction_methods.standard.rational_interpolant import (
			RationalInterpolant as RI)
			from rrompy.utilities.poly_fitting.polynomial import (
			polybases as ppb, polyfitname,
			polyvander as pvP, polyvanderTotal as pvTP,
			PolynomialInterpolator as PI)
			from rrompy.utilities.poly_fitting.radial_basis import (polybases as rbpb,
			RadialBasisInterpolator as RBI)
			from rrompy.utilities.poly_fitting.moving_least_squares import (
			polybases as mlspb,
			MovingLeastSquaresInterpolator as MLSI)
			from rrompy.reduction_methods.trained_model import (TrainedModelPivotedData,
			TrainedModelPivotedRational as tModel)
			from rrompy.utilities.base.types import (Np1D, Np2D, HFEng, DictAny, Tuple,
			List, ListAny, paramVal)
			from rrompy.utilities.base import verbosityManager as vbMng, freepar as fp
			from rrompy.utilities.numerical import (multifactorial, customPInv, dot,
			fullDegreeN, totalDegreeN,
			degreeTotalToFull, fullDegreeMaxMask,
			totalDegreeMaxMask,
			nextDerivativeIndices,
			hashDerivativeToIdx as hashD,
			hashIdxToDerivative as hashI)
			from rrompy.utilities.exception_manager import (RROMPyException, RROMPyAssert,
			RROMPyWarning)
			from rrompy.parameter import checkParameter

			__all__ = ['RationalInterpolantPivoted']

			class RationalInterpolantPivoted(GenericPivotedApproximant):
			"""
			ROM pivoted rational interpolant (with pole matching) computation for
			parametric problems.

			Args:
			HFEngine: HF problem solver.
			mu0(optional): Default parameter. Defaults to 0.
			directionPivot(optional): Pivot components. Defaults to [0].
			approxParameters(optional): Dictionary containing values for main
			parameters of approximant. Recognized keys are:
			- 'POD': whether to compute POD of snapshots; defaults to True;
			- 'matchingWeight': weight for pole matching optimization; defaults
			to 1;
			- 'cutOffTolerance': tolerance for ignoring parasitic poles;
			defaults to np.inf;
			- 'cutOffType': rule for tolerance computation for parasitic poles;
			defaults to 'MAGNITUDE';
			- 'S': total number of pivot samples current approximant relies
			upon;
			- 'samplerPivot': pivot sample point generator;
			- 'SMarginal': total number of marginal samples current approximant
			relies upon;
			- 'samplerMarginal': marginal sample point generator;
			- 'polybasisPivot': type of polynomial basis for pivot
			interpolation; defaults to 'MONOMIAL';
			- 'polybasisMarginal': type of polynomial basis for marginal
			interpolation; allowed values include 'MONOMIAL', 'CHEBYSHEV'
			and 'LEGENDRE'; defaults to 'MONOMIAL';
			- 'M': degree of rational interpolant numerator; defaults to 0;
			- 'N': degree of rational interpolant denominator; defaults to 0;
			- 'polydegreetype': type of polynomial degree; defaults to 'TOTAL';
			- 'MMarginal': degree of marginal interpolant; defaults to 0;
			- 'polydegreetypeMarginal': type of polynomial degree for marginal;
			defaults to 'TOTAL';
			- 'radialDirectionalWeightsPivot': radial basis weights for pivot
			numerator; defaults to 0, i.e. identity;
			- 'radialDirectionalWeightsMarginal': radial basis weights for
			marginal interpolant; defaults to 0, i.e. identity;
			- 'nNearestNeighborPivot': number of pivot nearest neighbors
			considered if polybasisPivot allows; defaults to -1;
			- 'nNearestNeighborMarginal': number of marginal nearest neighbors
			considered if polybasisMarginal allows; defaults to -1;
			- 'interpRcondPivot': tolerance for pivot interpolation; defaults
			to None;
			- 'interpRcondMarginal': tolerance for marginal interpolation;
			defaults to None;
			- 'robustTol': tolerance for robust rational denominator
			management; defaults to 0.
			Defaults to empty dict.
			verbosity(optional): Verbosity level. Defaults to 10.

			Attributes:
			HFEngine: HF problem solver.
			mu0: Default parameter.
			directionPivot: Pivot components.
			mus: Array of snapshot parameters.
			musPivot: Array of pivot snapshot parameters.
			musMarginal: Array of marginal snapshot parameters.
			approxParameters: Dictionary containing values for main parameters of
			approximant. Recognized keys are in parameterList.
			parameterListSoft: Recognized keys of soft approximant parameters:
			- 'POD': whether to compute POD of snapshots;
			- 'matchingWeight': weight for pole matching optimization;
			- 'cutOffTolerance': tolerance for ignoring parasitic poles;
			- 'cutOffType': rule for tolerance computation for parasitic poles;
			- 'polybasisPivot': type of polynomial basis for pivot
			interpolation;
			- 'polybasisMarginal': type of polynomial basis for marginal
			interpolation;
			- 'M': degree of rational interpolant numerator;
			- 'N': degree of rational interpolant denominator;
			- 'polydegreetype': type of polynomial degree;
			- 'MMarginal': degree of marginal interpolant;
			- 'polydegreetypeMarginal': type of polynomial degree for marginal;
			- 'radialDirectionalWeightsPivot': radial basis weights for pivot
			numerator;
			- 'radialDirectionalWeightsMarginal': radial basis weights for
			marginal interpolant;
			- 'nNearestNeighborPivot': number of pivot nearest neighbors
			considered if polybasisPivot allows;
			- 'nNearestNeighborMarginal': number of marginal nearest neighbors
			considered if polybasisMarginal allows;
			- 'interpRcondPivot': tolerance for pivot interpolation;
			- 'interpRcondMarginal': tolerance for marginal interpolation;
			- 'robustTol': tolerance for robust rational denominator
			management.
			parameterListCritical: Recognized keys of critical approximant
			parameters:
			- 'S': total number of pivot samples current approximant relies
			upon;
			- 'samplerPivot': pivot sample point generator;
			- 'SMarginal': total number of marginal samples current approximant
			relies upon;
			- 'samplerMarginal': marginal sample point generator.
			POD: Whether to compute POD of snapshots.
			matchingWeight: Weight for pole matching optimization.
			cutOffTolerance: Tolerance for ignoring parasitic poles.
			cutOffType: Rule for tolerance computation for parasitic poles.
			S: Total number of pivot samples current approximant relies upon.
			sampler: Pivot sample point generator.
			polybasisPivot: Type of polynomial basis for pivot interpolation.
			polybasisMarginal: Type of polynomial basis for marginal interpolation.
			M: Numerator degree of approximant.
			N: Denominator degree of approximant.
			polydegreetype: Type of polynomial degree.
			MMarginal: Degree of marginal interpolant.
			polydegreetypeMarginal: Type of polynomial degree for marginal.
			radialDirectionalWeightsPivot: Radial basis weights for pivot
			numerator.
			radialDirectionalWeightsMarginal: Radial basis weights for marginal
			interpolant.
			nNearestNeighborPivot: Number of pivot nearest neighbors considered if
			polybasisPivot allows.
			nNearestNeighborMarginal: Number of marginal nearest neighbors
			considered if polybasisMarginal allows.
			interpRcondPivot: Tolerance for pivot interpolation.
			interpRcondMarginal: Tolerance for marginal interpolation.
			robustTol: Tolerance for robust rational denominator management.
			muBoundsPivot: list of bounds for pivot parameter values.
			muBoundsMarginal: list of bounds for marginal parameter values.
			samplingEngine: Sampling engine.
			uHF: High fidelity solution(s) with parameter(s) lastSolvedHF as
			sampleList.
			lastSolvedHF: Parameter(s) corresponding to last computed high fidelity
			solution(s) as parameterList.
			uApproxReduced: Reduced approximate solution(s) with parameter(s)
			lastSolvedApprox as sampleList.
			lastSolvedApproxReduced: Parameter(s) corresponding to last computed
			reduced approximate solution(s) as parameterList.
			uApprox: Approximate solution(s) with parameter(s) lastSolvedApprox as
			sampleList.
			lastSolvedApprox: Parameter(s) corresponding to last computed
			approximate solution(s) as parameterList.
			Q: Numpy 1D vector containing complex coefficients of approximant
			denominator.
			P: Numpy 2D vector whose columns are FE dofs of coefficients of
			approximant numerator.
			"""

			def __init__(self, HFEngine:HFEng, mu0 : paramVal = None,
			directionPivot : ListAny = [0],
			approxParameters : DictAny = {}, verbosity : int = 10,
			timestamp : bool = True):
			self._preInit()
			self._addParametersToList(["polybasisPivot", "M", "N",
			"polydegreetype",
			"radialDirectionalWeightsPivot",
			"nNearestNeighborPivot",
			"interpRcondPivot", "robustTol"],
			["MONOMIAL", 0, 0, "TOTAL", 1, -1, -1, 0])
			super().__init__(HFEngine = HFEngine, mu0 = mu0,
			directionPivot = directionPivot,
			approxParameters = approxParameters,
			verbosity = verbosity, timestamp = timestamp)
			self._postInit()

			@property
			def polybasisPivot(self):
			"""Value of polybasisPivot."""
			return self._polybasisPivot
			@polybasisPivot.setter
			def polybasisPivot(self, polybasisPivot):
			try:
			polybasisPivot = polybasisPivot.upper().strip().replace(" ","")
			if polybasisPivot not in ppb + rbpb + mlspb:
			raise RROMPyException(
			"Prescribed pivot polybasis not recognized.")
			self._polybasisPivot = polybasisPivot
			except:
			RROMPyWarning(("Prescribed pivot polybasis not recognized. "
			"Overriding to 'MONOMIAL'."))
			self._polybasisPivot = "MONOMIAL"
			self._approxParameters["polybasisPivot"] = self.polybasisPivot

			@property
			def polybasisPivot0(self):
			if "_" in self.polybasisPivot:
			return self.polybasisPivot.split("_")[0]
			return self.polybasisPivot

			@property
			def radialDirectionalWeightsPivot(self):
			"""Value of radialDirectionalWeightsPivot."""
			return self._radialDirectionalWeightsPivot
			@radialDirectionalWeightsPivot.setter
			def radialDirectionalWeightsPivot(self, radialDirectionalWeightsPivot):
			self._radialDirectionalWeightsPivot = radialDirectionalWeightsPivot
			self._approxParameters["radialDirectionalWeightsPivot"] = (
			self.radialDirectionalWeightsPivot)

			@property
			def nNearestNeighborPivot(self):
			"""Value of nNearestNeighborPivot."""
			return self._nNearestNeighborPivot
			@nNearestNeighborPivot.setter
			def nNearestNeighborPivot(self, nNearestNeighborPivot):
			self._nNearestNeighborPivot = nNearestNeighborPivot
			self._approxParameters["nNearestNeighborPivot"] = (
			self.nNearestNeighborPivot)

			@property
			def interpRcondPivot(self):
			"""Value of interpRcondPivot."""
			return self._interpRcondPivot
			@interpRcondPivot.setter
			def interpRcondPivot(self, interpRcondPivot):
			self._interpRcondPivot = interpRcondPivot
			self._approxParameters["interpRcondPivot"] = self.interpRcondPivot

			@property
			def M(self):
			"""Value of M."""
			return self._M
			@M.setter
			def M(self, M):
			if M < 0: raise RROMPyException("M must be non-negative.")
			self._M = M
			self._approxParameters["M"] = self.M

			@property
			def N(self):
			"""Value of N."""
			return self._N
			@N.setter
			def N(self, N):
			if N < 0: raise RROMPyException("N must be non-negative.")
			self._N = N
			self._approxParameters["N"] = self.N

			@property
			def polydegreetype(self):
			"""Value of polydegreetype."""
			return self._polydegreetype
			@polydegreetype.setter
			def polydegreetype(self, polydegreetype):
			try:
			polydegreetype = polydegreetype.upper().strip().replace(" ","")
			if polydegreetype not in ["TOTAL", "FULL"]:
			raise RROMPyException(("Prescribed polydegreetype not "
			"recognized."))
			self._polydegreetype = polydegreetype
			except:
			RROMPyWarning(("Prescribed polydegreetype not recognized. "
			"Overriding to 'TOTAL'."))
			self._polydegreetype = "TOTAL"
			self._approxParameters["polydegreetype"] = self.polydegreetype

			@property
			def robustTol(self):
			"""Value of tolerance for robust rational denominator management."""
			return self._robustTol
			@robustTol.setter
			def robustTol(self, robustTol):
			if robustTol < 0.:
			RROMPyWarning(("Overriding prescribed negative robustness "
			"tolerance to 0."))
			robustTol = 0.
			self._robustTol = robustTol
			self._approxParameters["robustTol"] = self.robustTol

			def resetSamples(self):
			"""Reset samples."""
			super().resetSamples()
			self._musPUniqueCN = None
			self._derPIdxs = None
			self._reorderP = None

			def _setupPivotInterpolationIndices(self):
			"""Setup parameters for polyvander."""
			RROMPyAssert(self._mode,
			message = "Cannot setup interpolation indices.")
			if (self._musPUniqueCN is None
			or len(self._reorderP) != len(self.musPivot)):
			self._musPUniqueCN, musPIdxsTo, musPIdxs, musPCount = (
			self.trainedModel.centerNormalizePivot(self.musPivot).unique(
			return_index = True, return_inverse = True,
			return_counts = True))
			self._musPUnique = self.mus[musPIdxsTo]
			self._derPIdxs = [None] * len(self._musPUniqueCN)
			self._reorderP = np.empty(len(musPIdxs), dtype = int)
			filled = 0
			for j, cnt in enumerate(musPCount):
			self._derPIdxs[j] = nextDerivativeIndices([],
			self.nparPivot, cnt)
			jIdx = np.nonzero(musPIdxs == j)[0]
			self._reorderP[jIdx] = np.arange(filled, filled + cnt)
			filled += cnt

			def _setupDenominator(self):
			"""Compute rational denominator."""
			RROMPyAssert(self._mode, message = "Cannot setup denominator.")
			vbMng(self, "INIT", "Starting computation of denominator.", 7)
			NinvD = None
			N0 = copy(self.N)
			qs = []
			self.verbosity -= 10
			for j in range(len(self.musMarginal)):
			self._N = N0
			while self.N > 0:
			if NinvD != self.N:
			invD, fitinvP = self._computeInterpolantInverseBlocks()
			NinvD = self.N
			if self.POD:
			ev, eV = RI.findeveVGQR(self, self.samplingEngine.RPOD[j],
			invD)
			else:
			ev, eV = RI.findeveVGExplicit(self,
			self.samplingEngine.samples[j], invD)
			nevBad = checkRobustTolerance(ev, self.robustTol)
			if nevBad <= 1: break
			if self.catchInstability:
			raise RROMPyException(("Instability in denominator "
			"computation: eigenproblem is "
			"poorly conditioned."))
			RROMPyWarning(("Smallest {} eigenvalues below tolerance. "
			"Reducing N by 1.").format(nevBad))
			self.N = self.N - 1
			if self.N <= 0:
			self._N = 0
			eV = np.ones((1, 1))
			q = PI()
			q.npar = self.nparPivot
			q.polybasis = self.polybasisPivot0
			if self.polydegreetype == "TOTAL":
			q.coeffs = degreeTotalToFull(tuple([self.N + 1] * q.npar),
			q.npar, eV[:, 0])
			else:
			q.coeffs = eV[:, 0].reshape([self.N + 1] * q.npar)
			qs = qs + [copy(q)]
			self.verbosity += 10
			vbMng(self, "DEL", "Done computing denominator.", 7)
			return qs, fitinvP

			def _setupNumerator(self):
			"""Compute rational numerator."""
			RROMPyAssert(self._mode, message = "Cannot setup numerator.")
			vbMng(self, "INIT", "Starting computation of numerator.", 7)
			Qevaldiag = np.zeros((len(self.musPivot), len(self.musPivot)),
			dtype = np.complex)
			verb = self.trainedModel.verbosity
			self.trainedModel.verbosity = 0
			self._setupPivotInterpolationIndices()
			cfun = totalDegreeN if self.polydegreetype == "TOTAL" else fullDegreeN
			M = copy(self.M)
			while len(self.musPivot) < cfun(M, self.nparPivot): M -= 1
			if M < self.M:
			RROMPyWarning(("M too large compared to S. Reducing M by "
			"{}").format(self.M - M))
			self.M = M
			tensor_idx = 0
			ps = []
			for k, muM in enumerate(self.musMarginal):
			self._M = M
			idxGlob = 0
			for j, derIdxs in enumerate(self._derPIdxs):
			mujEff = [fp] * self.npar
			for jj, kk in enumerate(self.directionPivot):
			mujEff[kk] = self._musPUnique[j, jj]
			for jj, kk in enumerate(self.directionMarginal):
			mujEff[kk] = muM(0, jj)
			mujEff = checkParameter(mujEff, self.npar)
			nder = len(derIdxs)
			idxLoc = np.arange(len(self.musPivot))[
			(self._reorderP >= idxGlob)
			* (self._reorderP < idxGlob + nder)]
			idxGlob += nder
			Qval = [0] * nder
			for der in range(nder):
			derIdx = hashI(der, self.nparPivot)
			derIdxEff = [0] * self.npar
			sclEff = [0] * self.npar
			for jj, kk in enumerate(self.directionPivot):
			derIdxEff[kk] = derIdx[jj]
			sclEff[kk] = self.scaleFactorPivot[jj] ** -1.
			Qval[der] = (self.trainedModel.getQVal(mujEff, derIdxEff,
			scl = sclEff)
			/ multifactorial(derIdx))
			for derU, derUIdx in enumerate(derIdxs):
			for derQ, derQIdx in enumerate(derIdxs):
			diffIdx = [x - y for (x, y) in zip(derUIdx, derQIdx)]
			if all([x >= 0 for x in diffIdx]):
			diffj = hashD(diffIdx)
			Qevaldiag[idxLoc[derU], idxLoc[derQ]] = Qval[diffj]
			while self.M >= 0:
			if self.polybasisPivot in ppb:
			p = PI()
			wellCond, msg = p.setupByInterpolation(
			self._musPUniqueCN, Qevaldiag, self.M,
			self.polybasisPivot, self.verbosity >= 5,
			self.polydegreetype == "TOTAL",
			{"derIdxs": self._derPIdxs,
			"reorder": self._reorderP,
			"scl": np.power(self.scaleFactorPivot, -1.)},
			{"rcond": self.interpRcondPivot})
			elif self.polybasisPivot in rbpb:
			p = RBI()
			wellCond, msg = p.setupByInterpolation(
			self._musPUniqueCN, Qevaldiag, self.M,
			self.polybasisPivot,
			self.radialDirectionalWeightsPivot,
			self.verbosity >= 5,
			self.polydegreetype == "TOTAL",
			{"derIdxs": self._derPIdxs,
			"reorder": self._reorderP,
			"scl": np.power(self.scaleFactorPivot, -1.),
			"nNearestNeighbor" : self.nNearestNeighborPivot},
			{"rcond": self.interpRcondPivot})
			else:# if self.polybasisPivot in mlspb:
			p = MLSI()
			wellCond, msg = p.setupByInterpolation(
			self._musPUniqueCN, Qevaldiag, self.M,
			self.polybasisPivot,
			self.radialDirectionalWeightsPivot,
			self.verbosity >= 5,
			self.polydegreetype == "TOTAL",
			{"derIdxs": self._derPIdxs,
			"reorder": self._reorderP,
			"scl": np.power(self.scaleFactorPivot, -1.),
			"nNearestNeighbor" : self.nNearestNeighborPivot})
			vbMng(self, "MAIN", msg, 5)
			if wellCond: break
			if self.catchInstability:
			raise RROMPyException(("Instability in numerator "
			"computation: polyfit is "
			"poorly conditioned."))
			RROMPyWarning(("Polyfit is poorly conditioned. "
			"Reducing M by 1."))
			self.M = self.M - 1
			tensor_idx_new = tensor_idx + Qevaldiag.shape[1]
			if self.POD:
			p.postmultiplyTensorize(self.samplingEngine.RPODCoalesced.T[
			tensor_idx : tensor_idx_new, :])
			else:
			p.pad(tensor_idx, len(self.mus) - tensor_idx_new)
			tensor_idx = tensor_idx_new
			ps = ps + [copy(p)]
			self.trainedModel.verbosity = verb
			vbMng(self, "DEL", "Done computing numerator.", 7)
			return ps

			def setupApprox(self):
			"""
			Compute rational interpolant.
			SVD-based robust eigenvalue management.
			"""
			if self.checkComputedApprox():
			return
			RROMPyAssert(self._mode, message = "Cannot setup approximant.")
			vbMng(self, "INIT", "Setting up {}.". format(self.name()), 5)
			self.computeSnapshots()
			if self.trainedModel is None:
			self.trainedModel = tModel()
			self.trainedModel.verbosity = self.verbosity
			self.trainedModel.timestamp = self.timestamp
			data = TrainedModelPivotedData(self.trainedModel.name(), self.mu0,
			self.samplingEngine.samplesCoalesced,
			self.scaleFactor,
			self.HFEngine.rescalingExp,
			self.directionPivot)
			data.musPivot = copy(self.musPivot)
			data.musMarginal = copy(self.musMarginal)
			self.trainedModel.data = data
			else:
			self.trainedModel = self.trainedModel
			self.trainedModel.data.projMat = copy(
			self.samplingEngine.samplesCoalesced)
			self.trainedModel.data.marginalInterp = self._setupMarginalInterp()
			if self.N > 0:
			Qs = self._setupDenominator()[0]
			else:
			Q = PI()
			Q.npar = self.nparPivot
			Q.coeffs = np.ones(tuple([1] * Q.npar),
			dtype = self.musMarginal.dtype)
			Q.polybasis = self.polybasisPivot0
			Qs = [Q for _ in range(len(self.musMarginal))]
			self.trainedModel.data._temporary = True
			self.trainedModel.data.Qs = Qs
			self.trainedModel.data.Ps = self._setupNumerator()
			vbMng(self, "INIT", "Matching poles.", 10)
			self.trainedModel.initializeFromRational(self.HFEngine,
			self.matchingWeight, self.POD)
			vbMng(self, "DEL", "Done matching poles.", 10)
			if not np.isinf(self.cutOffTolerance):
			vbMng(self, "INIT", "Recompressing by cut-off.", 10)
			msg = self.trainedModel.recompressByCutOff([-1., 1.],
			self.cutOffTolerance,
			self.cutOffType)
			vbMng(self, "DEL", "Done recompressing." + msg, 10)
			self.trainedModel.data.approxParameters = copy(self.approxParameters)
			vbMng(self, "DEL", "Done setting up approximant.", 5)

			def _computeInterpolantInverseBlocks(self) -> Tuple[List[Np2D], Np2D]:
			"""
			Compute inverse factors for minimal interpolant target functional.
			"""
			RROMPyAssert(self._mode, message = "Cannot solve eigenvalue problem.")
			self._setupPivotInterpolationIndices()
			cfun = totalDegreeN if self.polydegreetype == "TOTAL" else fullDegreeN
			N = copy(self.N)
			while len(self.musPivot) < cfun(N, self.nparPivot): N -= 1
			if N < self.N:
			RROMPyWarning(("N too large compared to S. Reducing N by "
			"{}").format(self.N - N))
			self.N = N
			while self.N >= 0:
			if self.polydegreetype == "TOTAL":
			TE, _, argIdxs = pvTP(self._musPUniqueCN, self.N,
			self.polybasisPivot0, self._derPIdxs,
			self._reorderP,
			scl = np.power(self.scaleFactorPivot, -1.))
			TE = TE[:, argIdxs]
			idxsB = totalDegreeMaxMask(self.N, self.nparPivot)
			else: #if self.polydegreetype == "FULL":
			TE = pvP(self._musPUniqueCN, [self.N] * self.nparPivot,
			self.polybasisPivot0, self._derPIdxs, self._reorderP,
			scl = np.power(self.scaleFactorPivot, -1.))
			idxsB = fullDegreeMaxMask(self.N, self.nparPivot)
			fitOut = customPInv(TE, rcond = self.interpRcondPivot,
			full = True)
			vbMng(self, "MAIN",
			("Fitting {} samples with degree {} through {}... "
			"Conditioning of pseudoinverse system: {:.4e}.").format(
			TE.shape[0], self.N,
			polyfitname(self.polybasisPivot0),
			fitOut[1][1][0] / fitOut[1][1][-1]),
			5)
			if fitOut[1][0] == TE.shape[1]:
			fitinvP = fitOut[0][idxsB, :]
			break
			RROMPyWarning("Polyfit is poorly conditioned. Reducing N by 1.")
			self.N -= 1
			if self.N < 0:
			raise RROMPyException(("Instability in computation of "
			"denominator. Aborting."))
			TN, _, argIdxs = pvTP(self._musPUniqueCN, self.N, self.polybasisPivot0,
			self._derPIdxs, self._reorderP,
			scl = np.power(self.scaleFactorPivot, -1.))
			TN = TN[:, argIdxs]
			invD = [None] * (len(idxsB))
			for k in range(len(idxsB)):
			pseudoInv = np.diag(fitinvP[k, :])
			idxGlob = 0
			for j, derIdxs in enumerate(self._derPIdxs):
			nder = len(derIdxs)
			idxGlob += nder
			if nder > 1:
			idxLoc = np.arange(len(self.musPivot))[
			(self._reorderP >= idxGlob - nder)
			* (self._reorderP < idxGlob)]
			invLoc = fitinvP[k, idxLoc]
			pseudoInv[np.ix_(idxLoc, idxLoc)] = 0.
			for diffj, diffjIdx in enumerate(derIdxs):
			for derQ, derQIdx in enumerate(derIdxs):
			derUIdx = [x - y for (x, y) in
			zip(diffjIdx, derQIdx)]
			if all([x >= 0 for x in derUIdx]):
			derU = hashD(derUIdx)
			pseudoInv[idxLoc[derU], idxLoc[derQ]] = (
			invLoc[diffj])
			invD[k] = dot(pseudoInv, TN)
			return invD, fitinvP

			def getResidues(self, args, *kwargs) -> Np1D:
			"""
			Obtain approximant residues.

			Returns:
			Matrix with residues as columns.
			"""
			return self.trainedModel.getResidues(args, *kwargs)

RROMPy/rrompy/reduction_methods/pivoted/rational_interpolant_pivoted.py6643b38fd8aemaster

rational_interpolant_pivoted.pyView Options

RROMPy/rrompy/reduction_methods/pivoted/rational_interpolant_pivoted.py
6643b38fd8aemaster

rational_interpolant_pivoted.py
View Options