approximant_lagrange_pade.py
No OneTemporary
Actions

Subscribers

None

File Metadata

Created: Thu, May 2, 12:22

approximant_lagrange_pade.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 copy
	import numpy as np
	from scipy.special import factorial as fact
	from rrompy.reduction_methods.base import checkRobustTolerance
	from .generic_approximant_lagrange import GenericApproximantLagrange
	from rrompy.utilities.poly_fitting import (polybases, polyvander, polyfitname,
	customFit)
	from rrompy.reduction_methods.trained_model import TrainedModelPade as tModel
	from rrompy.reduction_methods.trained_model import TrainedModelData
	from rrompy.utilities.base.types import Np1D, Np2D, HFEng, DictAny, Tuple
	from rrompy.utilities.base import verbosityDepth, purgeDict
	from rrompy.utilities.exception_manager import (RROMPyException, modeAssert,
	RROMPyWarning)

	__all__ = ['ApproximantLagrangePade']

	class ApproximantLagrangePade(GenericApproximantLagrange):
	"""
	ROM Lagrange Pade' interpolant computation for parametric problems.

	Args:
	HFEngine: HF problem solver.
	mu0(optional): Default parameter. 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;
	- 'S': total number of samples current approximant relies upon;
	defaults to 2;
	- 'sampler': sample point generator; defaults to uniform sampler on
	[0, 1];
	- 'polybasis': type of polynomial basis for interpolation; allowed
	values include 'MONOMIAL', 'CHEBYSHEV' and 'LEGENDRE'; defaults
	to 'MONOMIAL';
	- 'E': coefficient of interpolant to be minimized; defaults to
	min(S, M + 1);
	- 'M': degree of Pade' interpolant numerator; defaults to 0;
	- 'N': degree of Pade' interpolant denominator; defaults to 0;
	- 'interpRcond': tolerance for interpolation via numpy.polyfit;
	defaults to None;
	- 'robustTol': tolerance for robust Pade' denominator management;
	defaults to 0.
	Defaults to empty dict.
	homogeneized: Whether to homogeneize Dirichlet BCs. Defaults to False.
	verbosity(optional): Verbosity level. Defaults to 10.

	Attributes:
	HFEngine: HF problem solver.
	mu0: Default parameter.
	mus: Array of snapshot parameters.
	ws: Array of snapshot weigths.
	homogeneized: Whether to homogeneize Dirichlet BCs.
	approxParameters: Dictionary containing values for main parameters of
	approximant. Recognized keys are in parameterList.
	parameterList: Recognized keys of approximant parameters:
	- 'POD': whether to compute POD of snapshots;
	- 'S': total number of samples current approximant relies upon;
	- 'sampler': sample point generator;
	- 'polybasis': type of polynomial basis for interpolation;
	- 'E': coefficient of interpolant to be minimized;
	- 'M': degree of Pade' interpolant numerator;
	- 'N': degree of Pade' interpolant denominator;
	- 'interpRcond': tolerance for interpolation via numpy.polyfit;
	- 'robustTol': tolerance for robust Pade' denominator management.
	extraApproxParameters: List of approxParameters keys in addition to
	mother class's.
	S: Number of solution snapshots over which current approximant is
	based upon.
	sampler: Sample point generator.
	M: Numerator degree of approximant.
	N: Denominator degree of approximant.
	POD: Whether to compute POD of snapshots.
	interpRcond: Tolerance for interpolation via numpy.polyfit.
	robustTol: Tolerance for robust Pade' denominator management.
	samplingEngine: Sampling engine.
	uHF: High fidelity solution with wavenumber lastSolvedHF as numpy
	complex vector.
	lastSolvedHF: Wavenumber corresponding to last computed high fidelity
	solution.
	Q: Numpy 1D vector containing complex coefficients of approximant
	denominator.
	P: Numpy 2D vector whose columns are FE dofs of coefficients of
	approximant numerator.
	uApp: Last evaluated approximant as numpy complex vector.
	"""

	def __init__(self, HFEngine:HFEng, mu0 : complex = 0.,
	approxParameters : DictAny = {}, homogeneized : bool = False,
	verbosity : int = 10, timestamp : bool = True):
	self.__preInit()
	self.__addParametersToList(["polybasis", "E", "M", "N",
	"interpRcond", "robustTol"])
	super().__init__(HFEngine = HFEngine, mu0 = mu0,
	approxParameters = approxParameters,
	homogeneized = homogeneized,
	verbosity = verbosity, timestamp = timestamp)
	self.__postInit()

	@property
	def approxParameters(self):
	"""
	Value of approximant parameters.
	"""
	return self._approxParameters
	@approxParameters.setter
	def approxParameters(self, approxParams):
	approxParameters = purgeDict(approxParams, self.parameterList,
	dictname = self.name() + ".approxParameters",
	baselevel = 1)
	approxParametersCopy = purgeDict(approxParameters, ["polybasis",
	"E", "M", "N",
	"interpRcond",
	"robustTol"],
	True, True, baselevel = 1)
	if hasattr(self, "_M") and self.M is not None:
	Mold = self.M
	self._M = 0
	if hasattr(self, "_N") and self.N is not None:
	Nold = self.N
	self._N = 0
	if hasattr(self, "_E") and self.E is not None:
	self._E = 0
	GenericApproximantLagrange.approxParameters.fset(self,
	approxParametersCopy)
	keyList = list(approxParameters.keys())
	if "polybasis" in keyList:
	self.polybasis = approxParameters["polybasis"]
	elif not hasattr(self, "_polybasis") or self._polybasis is None:
	self.polybasis = "MONOMIAL"
	if "interpRcond" in keyList:
	self.interpRcond = approxParameters["interpRcond"]
	elif not hasattr(self, "interpRcond") or self.interpRcond is None:
	self.interpRcond = None
	if "robustTol" in keyList:
	self.robustTol = approxParameters["robustTol"]
	elif not hasattr(self, "_robustTol") or self._robustTol is None:
	self.robustTol = 0
	if "M" in keyList:
	self.M = approxParameters["M"]
	elif hasattr(self, "_M") and self.M is not None:
	self.M = Mold
	else:
	self.M = 0
	if "N" in keyList:
	self.N = approxParameters["N"]
	elif hasattr(self, "_N") and self.N is not None:
	self.N = Nold
	else:
	self.N = 0
	if "E" in keyList:
	self.E = approxParameters["E"]
	else:
	self.E = min(self.S - 1, self.M + 1)

	@property
	def polybasis(self):
	"""Value of polybasis."""
	return self._polybasis
	@polybasis.setter
	def polybasis(self, polybasis):
	try:
	polybasis = polybasis.upper().strip().replace(" ","")
	if polybasis not in polybases:
	raise RROMPyException("Prescribed polybasis not recognized.")
	self._polybasis = polybasis
	except:
	RROMPyWarning(("Prescribed polybasis not recognized. Overriding "
	"to 'MONOMIAL'."))
	self._sampleType = "MONOMIAL"
	self._approxParameters["polybasis"] = self.polybasis

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

	@property
	def M(self):
	"""Value of M. Its assignment may change S."""
	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
	if hasattr(self, "_S") and self.S < self.M + 1:
	RROMPyWarning("Prescribed S is too small. Updating S to M + 1.")
	self.S = self.M + 1

	@property
	def N(self):
	"""Value of N. Its assignment may change S."""
	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
	if hasattr(self, "_S") and self.S < self.N + 1:
	RROMPyWarning("Prescribed S is too small. Updating S to N + 1.")
	self.S = self.N + 1

	@property
	def E(self):
	"""Value of E. Its assignment may change S."""
	return self._E
	@E.setter
	def E(self, E):
	if E < 0: raise RROMPyException("E must be non-negative.")
	self._E = E
	self._approxParameters["E"] = self.E
	if hasattr(self, "_S") and self.S < self.E + 1:
	RROMPyWarning("Prescribed S is too small. Updating S to E + 1.")
	self.S = self.E + 1

	@property
	def robustTol(self):
	"""Value of tolerance for robust Pade' 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

	@property
	def S(self):
	"""Value of S."""
	return self._S
	@S.setter
	def S(self, S):
	if S <= 0: raise RROMPyException("S must be positive.")
	if hasattr(self, "_S"): Sold = self.S
	else: Sold = -1
	vals, label = [0] * 3, {0:"M", 1:"N", 2:"E"}
	if hasattr(self, "_M") and self._M is not None: vals[0] = self.M
	if hasattr(self, "_N") and self._N is not None: vals[1] = self.N
	if hasattr(self, "_E") and self._E is not None: vals[2] = self.E
	idxmax = np.argmax(vals)
	if vals[idxmax] + 1 > S:
	RROMPyWarning(("Prescribed S is too small. Updating S to {} + "
	"1.").format(label[idxmax]))
	self.S = vals[idxmax] + 1
	else:
	self._S = S
	self._approxParameters["S"] = self.S
	if Sold != self.S:
	self.resetSamples()

	def __setupDenominator(self):
	"""Compute Pade' denominator."""
	if self.verbosity >= 7:
	verbosityDepth("INIT", "Starting computation of denominator.",
	timestamp = self.timestamp)
	while self.N > 0:
	TE = polyvander[self.polybasis](self.radiusPade(self.mus), self.E,
	scl = 1. / self.scaleFactor)
	TE = (TE.T * self.ws).T
	RHS = np.zeros(self.E + 1)
	RHS[-1] = 1.
	fitOut = customFit(TE.T, RHS, full = True,
	rcond = self.interpRcond)
	if self.verbosity >= 5:
	condfit = fitOut[1][2][0] / fitOut[1][2][-1]
	verbosityDepth("MAIN", ("Fitting {} samples with degree {} "
	"through {}... Conditioning of LS "
	"system: {:.4e}.").format(
	self.S, self.E,
	polyfitname[self.polybasis],
	condfit),
	timestamp = self.timestamp)
	if fitOut[1][1] < self.E + 1:
	Enew = fitOut[1][1] - 1
	Nnew = min(self.N, Enew)
	Mnew = min(self.M, Enew)
	if Nnew == self.N:
	strN = ""
	else:
	strN = "N from {} to {} and ".format(self.N, Nnew)
	if Mnew == self.M:
	strM = ""
	else:
	strM = "M from {} to {} and ".format(self.M, Mnew)
	RROMPyWarning(("Polyfit is poorly conditioned.\nReducing {}{}"
	"E from {} to {}.").format(strN, strM,
	self.E, Enew))
	newParams = {"N" : Nnew, "M" : Mnew, "E" : Enew}
	self.approxParameters = newParams
	continue
	mus_un, idx_un, cnt_un = np.unique(self.mus, return_inverse = True,
	return_counts = True)
	TE = polyvander[self.polybasis](self.radiusPade(self.mus), self.N,
	scl = 1. / self.scaleFactor)
	TE = (TE.T * self.ws).T
	if len(mus_un) == len(self.mus):
	Ghalf = (TE.T * fitOut[0]).T
	else:
	pseudoInv = np.zeros((len(self.mus), len(self.mus)),
	dtype = np.complex)
	for j in range(len(mus_un)):
	pseudoInv_loc = np.zeros((cnt_un[j], cnt_un[j]),
	dtype = np.complex)
	mask = np.arange(len(self.mus))[idx_un == j]
	for der in range(cnt_un[j]):
	fitderj = fitOut[0][mask[der]]
	pseudoInv_loc = (pseudoInv_loc + fitderj
	* np.diag(np.ones(1 + der),
	k = der - cnt_un[j] + 1))
	I = np.ix_(mask, mask)
	pseudoInv[I] = np.flipud(pseudoInv_loc)
	Ghalf = pseudoInv.dot(TE)
	if self.POD:
	self.Ghalf = self.samplingEngine.RPOD.dot(Ghalf)
	ev, eV = self.findeveVGQR()
	else:
	self.Ghalf = self.samplingEngine.samples.dot(Ghalf)
	ev, eV = self.findeveVGExplicit()
	newParams = checkRobustTolerance(ev, self.E, self.robustTol)
	if not newParams:
	break
	self.approxParameters = newParams
	if self.N <= 0:
	self._N = 0
	eV = np.ones((1, 1))
	if self.verbosity >= 7:
	verbosityDepth("DEL", "Done computing denominator.",
	timestamp = self.timestamp)
	return eV[:, 0]

	def __setupNumerator(self):
	"""Compute Pade' numerator."""
	if self.verbosity >= 7:
	verbosityDepth("INIT", "Starting computation of numerator.",
	timestamp = self.timestamp)
	Qevaldiag = np.diag(self.trainedModel.getQVal(self.mus))
	verb = self.trainedModel.verbosity
	self.trainedModel.verbosity = 0
	mus_un, idx_un, cnt_un = np.unique(self.mus, return_inverse = True,
	return_counts = True)
	for j in range(len(mus_un)):
	if cnt_un[j] > 1:
	Q_loc = np.zeros((cnt_un[j], cnt_un[j]), dtype = np.complex)
	for der in range(1, cnt_un[j]):
	Qderj = (self.trainedModel.getQVal(mus_un[j], der,
	scl = 1. / self.scaleFactor)
	/ fact(der))
	Q_loc = Q_loc + Qderj * np.diag(np.ones(cnt_un[j] - der),
	k = - der)
	I = np.ix_(idx_un == j, idx_un == j)
	Qevaldiag[I] = Qevaldiag[I] + Q_loc
	self.trainedModel.verbosity = verb
	while self.M >= 0:
	fitVander = polyvander[self.polybasis](self.radiusPade(self.mus),
	self.M,
	scl = 1. / self.scaleFactor)
	fitOut = customFit(fitVander, Qevaldiag, w = self.ws, full = True,
	rcond = self.interpRcond)
	if self.verbosity >= 5:
	condfit = fitOut[1][2][0] / fitOut[1][2][-1]
	verbosityDepth("MAIN", ("Fitting {} samples with degree {} "
	"through {}... Conditioning of LS "
	"system: {:.4e}.").format(
	self.S, self.M,
	polyfitname[self.polybasis],
	condfit),
	timestamp = self.timestamp)
	if fitOut[1][1] == self.M + 1:
	P = fitOut[0].T
	break
	RROMPyWarning(("Polyfit is poorly conditioned. Reducing M from {} "
	"to {}. Exact snapshot interpolation not "
	"guaranteed.").format(self.M, fitOut[1][1] - 1))
	self.M = fitOut[1][1] - 1
	if self.M <= 0:
	raise RROMPyException(("Instability in computation of numerator. "
	"Aborting."))
	if self.verbosity >= 7:
	verbosityDepth("DEL", "Done computing numerator.",
	timestamp = self.timestamp)
	return np.atleast_2d(P)

	def setupApprox(self):
	"""
	Compute Pade' interpolant.
	SVD-based robust eigenvalue management.
	"""
	if self.checkComputedApprox():
	return
	modeAssert(self._mode, message = "Cannot setup approximant.")
	if self.verbosity >= 5:
	verbosityDepth("INIT", "Setting up {}.". format(self.name()),
	timestamp = self.timestamp)
	self.computeScaleFactor()
	self.computeSnapshots()
	if self.trainedModel is None:
	self.trainedModel = tModel()
	self.trainedModel.verbosity = self.verbosity
	self.trainedModel.timestamp = self.timestamp
	data = TrainedModelData(self.trainedModel.name(), self.mu0,
	np.copy(self.samplingEngine.samples),
	self.HFEngine.rescalingExp)
	data.polytype = self.polybasis
	data.scaleFactor = self.scaleFactor
	data.mus = np.copy(self.mus)
	self.trainedModel.data = data
	else:
	self.trainedModel.data.projMat = np.copy(
	self.samplingEngine.samples)
	if self.N > 0:
	Q = self.__setupDenominator()
	else:
	Q = np.ones(1, dtype = np.complex)
	self.trainedModel.data.Q = np.copy(Q)
	P = self.__setupNumerator()
	if self.POD:
	P = self.samplingEngine.RPOD.dot(P)
	self.trainedModel.data.P = np.copy(P)
	self.trainedModel.data.approxParameters = copy(self.approxParameters)
	if self.verbosity >= 5:
	verbosityDepth("DEL", "Done setting up approximant.",
	timestamp = self.timestamp)

	def findeveVGExplicit(self, verbOutput : int = 5) -> Tuple[Np1D, Np2D]:
	"""
	Compute explicitly eigenvalues and eigenvectors of Pade' denominator
	matrix.
	"""
	modeAssert(self._mode, message = "Cannot solve eigenvalue problem.")
	if self.verbosity >= 10:
	verbosityDepth("INIT", "Building gramian matrix.",
	timestamp = self.timestamp)
	self.G = self.HFEngine.innerProduct(self.Ghalf, self.Ghalf)
	if self.verbosity >= 10:
	verbosityDepth("DEL", "Done building gramian.",
	timestamp = self.timestamp)
	if self.verbosity >= 7:
	verbosityDepth("INIT", ("Solving eigenvalue problem for gramian "
	"matrix."), timestamp = self.timestamp)
	ev, eV = np.linalg.eigh(self.G)
	if self.verbosity >= verbOutput:
	try: condev = ev[-1] / ev[0]
	except: condev = np.inf
	verbosityDepth("MAIN", ("Solved eigenvalue problem of size {} "
	"with condition number {:.4e}.").format(
	self.N + 1, condev),
	timestamp = self.timestamp)
	if self.verbosity >= 7:
	verbosityDepth("DEL", "Done solving eigenvalue problem.",
	timestamp = self.timestamp)
	return ev, eV

	def findeveVGQR(self, verbOutput : int = 5) -> Tuple[Np1D, Np2D]:
	"""
	Compute eigenvalues and eigenvectors of Pade' denominator matrix
	through SVD of R factor.
	"""
	if self.verbosity >= 7:
	verbosityDepth("INIT", ("Solving svd for square root of gramian "
	"matrix."), timestamp = self.timestamp)
	_, s, eV = np.linalg.svd(self.Ghalf, full_matrices = False)
	ev = s[::-1]
	eV = eV[::-1, :].T.conj()
	if self.verbosity >= verbOutput:
	try: condev = s[0] / s[-1]
	except: condev = np.inf
	verbosityDepth("MAIN", ("Solved svd problem of size {} x {} with "
	"condition number {:.4e}.").format(
	self.S, self.N + 1, condev),
	timestamp = self.timestamp)
	if self.verbosity >= 7:
	verbosityDepth("DEL", "Done solving eigenvalue problem.",
	timestamp = self.timestamp)
	return ev, eV

	def radiusPade(self, mu:Np1D, mu0 : float = None) -> float:
	"""
	Compute translated radius to be plugged into Pade' approximant.

	Args:
	mu: Parameter(s) 1.
	mu0: Parameter(s) 2. If None, set to self.mu0.

	Returns:
	Translated radius to be plugged into Pade' approximant.
	"""
	return self.trainedModel.radiusPade(mu, mu0)

	def getResidues(self) -> Np1D:
	"""
	Obtain approximant residues.

	Returns:
	Matrix with residues as columns.
	"""
	return self.trainedModel.getResidues()

approximant_lagrange_pade.pyNo OneTemporaryActions

File Metadata

approximant_lagrange_pade.pyView Options

Event Timeline

approximant_lagrange_pade.py
No OneTemporary
Actions

approximant_lagrange_pade.py
View Options