generic_greedy_approximant.py
No OneTemporary
Actions

Subscribers

None

File Metadata

Created: Sat, Apr 27, 14:52

generic_greedy_approximant.py
View Options

	# Copyright (C) 2018-2020 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 abc import abstractmethod
	from copy import deepcopy as copy
	import numpy as np
	from matplotlib import pyplot as plt
	from rrompy.hfengines.base.linear_affine_engine import checkIfAffine
	from rrompy.reduction_methods.standard.generic_standard_approximant import (
	GenericStandardApproximant)
	from rrompy.utilities.base.types import (Np1D, Np2D, Tuple, List, paramVal,
	paramList, sampList)
	from rrompy.utilities.base import verbosityManager as vbMng
	from rrompy.utilities.numerical import dot
	from rrompy.utilities.expression import expressionEvaluator
	from rrompy.utilities.exception_manager import (RROMPyException, RROMPyAssert,
	RROMPyWarning)
	from rrompy.sampling.sample_list import sampleList
	from rrompy.parameter import emptyParameterList, parameterList
	from rrompy.utilities.parallel import masterCore

	__all__ = ['GenericGreedyApproximant']

	def localL2Distance(mus:Np2D, badmus:Np2D) -> Np2D:
	return np.linalg.norm(np.tile(mus[..., np.newaxis], [1, 1, len(badmus)])
	- badmus[..., np.newaxis].T, axis = 1)

	def pruneSamples(mus:paramList, badmus:paramList,
	tol : float = 1e-8) -> Np1D:
	"""Remove from mus all the elements which are too close to badmus."""
	if isinstance(mus, (parameterList, sampleList)): mus = mus.data
	if isinstance(badmus, (parameterList, sampleList)): badmus = badmus.data
	if len(badmus) == 0: return np.arange(len(mus))
	proximity = np.min(localL2Distance(mus, badmus), axis = 1)
	return np.where(proximity <= tol)[0]

	class GenericGreedyApproximant(GenericStandardApproximant):
	"""
	ROM greedy interpolant computation for parametric problems
	(ABSTRACT).

	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': kind of snapshots orthogonalization; allowed values
	include 0, 1/2, and 1; defaults to 1, i.e. POD;
	- 'scaleFactorDer': scaling factors for derivative computation;
	defaults to 'AUTO';
	- 'S': number of starting training points;
	- 'sampler': sample point generator;
	- 'greedyTol': uniform error tolerance for greedy algorithm;
	defaults to 1e-2;
	- 'collinearityTol': collinearity tolerance for greedy algorithm;
	defaults to 0.;
	- 'maxIter': maximum number of greedy steps; defaults to 1e2;
	- 'nTestPoints': number of test points; defaults to 5e2;
	- 'samplerTrainSet': training sample points generator; defaults to
	sampler.
	Defaults to empty dict.
	verbosity(optional): Verbosity level. Defaults to 10.

	Attributes:
	HFEngine: HF problem solver.
	mu0: Default parameter.
	mus: Array of snapshot parameters.
	approxParameters: Dictionary containing values for main parameters of
	approximant. Recognized keys are in parameterList.
	parameterListSoft: Recognized keys of soft approximant parameters:
	- 'POD': kind of snapshots orthogonalization;
	- 'scaleFactorDer': scaling factors for derivative computation;
	- 'greedyTol': uniform error tolerance for greedy algorithm;
	- 'collinearityTol': collinearity tolerance for greedy algorithm;
	- 'maxIter': maximum number of greedy steps;
	- 'nTestPoints': number of test points;
	- 'samplerTrainSet': training sample points generator.
	parameterListCritical: Recognized keys of critical approximant
	parameters:
	- 'S': total number of samples current approximant relies upon;
	- 'sampler': sample point generator.
	verbosity: Verbosity level.
	POD: Kind of snapshots orthogonalization.
	scaleFactorDer: Scaling factors for derivative computation.
	S: number of test points.
	sampler: Sample point generator.
	greedyTol: Uniform error tolerance for greedy algorithm.
	collinearityTol: Collinearity tolerance for greedy algorithm.
	maxIter: maximum number of greedy steps.
	nTestPoints: number of starting training points.
	samplerTrainSet: training sample points generator.
	muBounds: list of bounds for 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.
	"""

	def __init__(self, args, *kwargs):
	self._preInit()
	if not hasattr(self, "_affine_lvl"): self._affine_lvl = []
	self._affine_lvl += [1]
	self._addParametersToList(["greedyTol", "collinearityTol", "maxIter",
	"nTestPoints", "samplerTrainSet"],
	[1e-2, 0., 1e2, 5e2, "AUTO"])
	super().__init__(args, *kwargs)
	self._postInit()

	@property
	def greedyTol(self):
	"""Value of greedyTol."""
	return self._greedyTol
	@greedyTol.setter
	def greedyTol(self, greedyTol):
	if greedyTol < 0:
	raise RROMPyException("greedyTol must be non-negative.")
	if hasattr(self, "_greedyTol") and self.greedyTol is not None:
	greedyTolold = self.greedyTol
	else:
	greedyTolold = -1
	self._greedyTol = greedyTol
	self._approxParameters["greedyTol"] = self.greedyTol
	if greedyTolold != self.greedyTol:
	self.resetSamples()

	@property
	def collinearityTol(self):
	"""Value of collinearityTol."""
	return self._collinearityTol
	@collinearityTol.setter
	def collinearityTol(self, collinearityTol):
	if collinearityTol < 0:
	raise RROMPyException("collinearityTol must be non-negative.")
	if (hasattr(self, "_collinearityTol")
	and self.collinearityTol is not None):
	collinearityTolold = self.collinearityTol
	else:
	collinearityTolold = -1
	self._collinearityTol = collinearityTol
	self._approxParameters["collinearityTol"] = self.collinearityTol
	if collinearityTolold != self.collinearityTol:
	self.resetSamples()

	@property
	def maxIter(self):
	"""Value of maxIter."""
	return self._maxIter
	@maxIter.setter
	def maxIter(self, maxIter):
	if maxIter <= 0: raise RROMPyException("maxIter must be positive.")
	if hasattr(self, "_maxIter") and self.maxIter is not None:
	maxIterold = self.maxIter
	else:
	maxIterold = -1
	self._maxIter = maxIter
	self._approxParameters["maxIter"] = self.maxIter
	if maxIterold != self.maxIter:
	self.resetSamples()

	@property
	def nTestPoints(self):
	"""Value of nTestPoints."""
	return self._nTestPoints
	@nTestPoints.setter
	def nTestPoints(self, nTestPoints):
	if nTestPoints <= 0:
	raise RROMPyException("nTestPoints must be positive.")
	if not np.isclose(nTestPoints, np.int(nTestPoints)):
	raise RROMPyException("nTestPoints must be an integer.")
	nTestPoints = np.int(nTestPoints)
	if hasattr(self, "_nTestPoints") and self.nTestPoints is not None:
	nTestPointsold = self.nTestPoints
	else:
	nTestPointsold = -1
	self._nTestPoints = nTestPoints
	self._approxParameters["nTestPoints"] = self.nTestPoints
	if nTestPointsold != self.nTestPoints:
	self.resetSamples()

	@property
	def samplerTrainSet(self):
	"""Value of samplerTrainSet."""
	return self._samplerTrainSet
	@samplerTrainSet.setter
	def samplerTrainSet(self, samplerTrainSet):
	if (isinstance(samplerTrainSet, (str,))
	and samplerTrainSet.upper() == "AUTO"):
	samplerTrainSet = self.sampler
	if 'generatePoints' not in dir(samplerTrainSet):
	raise RROMPyException("samplerTrainSet type not recognized.")
	if (hasattr(self, '_samplerTrainSet')
	and self.samplerTrainSet not in [None, "AUTO"]):
	samplerTrainSetOld = self.samplerTrainSet
	self._samplerTrainSet = samplerTrainSet
	self._approxParameters["samplerTrainSet"] = self.samplerTrainSet
	if (not 'samplerTrainSetOld' in locals()
	or samplerTrainSetOld != self.samplerTrainSet):
	self.resetSamples()

	def resetSamples(self):
	"""Reset samples."""
	super().resetSamples()
	self._mus = emptyParameterList()

	def _affineResidualMatricesContraction(self, rb:Np2D, rA : Np2D = None) \
	-> Tuple[Np1D, Np1D, Np1D]:
	self.assembleReducedResidualBlocks(full = rA is not None)
	# 'ij,jk,ik->k', resbb, radiusb, radiusb.conj()
	ff = np.sum(self.trainedModel.data.resbb.dot(rb) * rb.conj(), axis = 0)
	if rA is None: return ff
	# 'ijk,jkl,il->l', resAb, radiusA, radiusb.conj()
	Lf = np.sum(np.tensordot(self.trainedModel.data.resAb, rA, 2)
	* rb.conj(), axis = 0)
	# 'ijkl,klt,ijt->t', resAA, radiusA, radiusA.conj()
	LL = np.sum(np.tensordot(self.trainedModel.data.resAA, rA, 2)
	* rA.conj(), axis = (0, 1))
	return ff, Lf, LL

	def getErrorEstimatorAffine(self, mus:Np1D) -> Np1D:
	"""Standard residual estimator."""
	checkIfAffine(self.HFEngine, "apply affinity-based error estimator",
	False, self._affine_lvl)
	self.HFEngine.buildA()
	self.HFEngine.buildb()
	mus = self.checkParameterList(mus)
	tMverb, self.trainedModel.verbosity = self.trainedModel.verbosity, 0
	uApproxRs = self.getApproxReduced(mus).data
	self.trainedModel.verbosity = tMverb
	muTestEff = self.mapParameterList(mus)
	radiusA = np.empty((len(self.HFEngine.thAs), len(mus)),
	dtype = np.complex)
	radiusb = np.empty((len(self.HFEngine.thbs), len(mus)),
	dtype = np.complex)
	for j, thA in enumerate(self.HFEngine.thAs):
	radiusA[j] = expressionEvaluator(thA[0], muTestEff)
	for j, thb in enumerate(self.HFEngine.thbs):
	radiusb[j] = expressionEvaluator(thb[0], muTestEff)
	radiusA = np.expand_dims(uApproxRs, 1) * radiusA
	ff, Lf, LL = self._affineResidualMatricesContraction(radiusb, radiusA)
	err = np.abs((LL - 2. * np.real(Lf) + ff) / ff) ** .5
	return err

	def errorEstimator(self, mus:Np1D, return_max : bool = False) -> Np1D:
	setupOK = self.setupApproxLocal()
	if setupOK > 0:
	err = np.empty(len(mus))
	err[:] = np.nan
	if not return_max: return err
	return err, [- setupOK], np.nan
	mus = self.checkParameterList(mus)
	vbMng(self.trainedModel, "INIT",
	"Evaluating error estimator at mu = {}.".format(mus), 10)
	err = self.getErrorEstimatorAffine(mus)
	vbMng(self.trainedModel, "DEL", "Done evaluating error estimator.", 10)
	if not return_max: return err
	idxMaxEst = [np.argmax(err)]
	return err, idxMaxEst, err[idxMaxEst]

	def _isLastSampleCollinear(self) -> bool:
	"""Check collinearity of last sample."""
	if self.collinearityTol <= 0.: return False
	if self.POD == 1:
	reff = self.samplingEngine.Rscale[:, -1]
	else:
	RROMPyWarning(("Repeated orthogonalization of the samples for "
	"collinearity check. Consider setting POD to "
	"True."))
	if not hasattr(self, "_PODEngine"):
	from rrompy.sampling import PODEngine
	self._PODEngine = PODEngine(self.HFEngine)
	reff = self._PODEngine.generalizedQR(self.samplingEngine.samples,
	only_R = True,
	is_state = True)[:, -1]
	cLevel = np.abs(reff[-1]) / np.linalg.norm(reff)
	cLevel = np.inf if np.isclose(cLevel, 0., atol = 1e-15) else 1 / cLevel
	vbMng(self, "MAIN", "Collinearity indicator {:.4e}.".format(cLevel), 3)
	return cLevel > self.collinearityTol

	def plotEstimator(self, est:Np1D, idxMax:List[int], estMax:List[float]):
	if (not (np.any(np.isnan(est)) or np.any(np.isinf(est)))
	and masterCore()):
	fig = plt.figure(figsize = plt.figaspect(1. / self.npar))
	for jpar in range(self.npar):
	ax = fig.add_subplot(1, self.npar, 1 + jpar)
	musre = np.array(self.muTest.re.data)
	errCP = copy(est)
	idx = np.delete(np.arange(self.npar), jpar)
	while len(musre) > 0:
	if self.npar == 1:
	currIdx = np.arange(len(musre))
	else:
	currIdx = np.where(np.isclose(np.sum(
	np.abs(musre[:, idx] - musre[0, idx]), 1),
	0., atol = 1e-15))[0]
	ax.semilogy(musre[currIdx, jpar], errCP[currIdx], 'k',
	linewidth = 1)
	musre = np.delete(musre, currIdx, 0)
	errCP = np.delete(errCP, currIdx)
	ax.semilogy([self.muBounds.re(0, jpar),
	self.muBounds.re(-1, jpar)],
	[self.greedyTol] * 2, 'r--')
	ax.semilogy(self.mus.re(jpar),
	2. * self.greedyTol * np.ones(len(self.mus)), '*m')
	if len(idxMax) > 0 and estMax is not None:
	ax.semilogy(self.muTest.re(idxMax, jpar), estMax, 'xr')
	ax.set_xlim(*list(self.sampler.lims.re(jpar)))
	ax.grid()
	plt.tight_layout()
	plt.show()

	def greedyNextSample(self, muidx:int, plotEst : str = "NONE")\
	-> Tuple[Np1D, int, float, paramVal]:
	"""Compute next greedy snapshot of solution map."""
	RROMPyAssert(self._mode, message = "Cannot add greedy sample.")
	mus = copy(self.muTest[muidx])
	self.muTest.pop(muidx)
	for j, mu in enumerate(mus):
	vbMng(self, "MAIN",
	("Adding sample point no. {} at {} to training "
	"set.").format(len(self.mus) + 1, mu), 3)
	self.mus.append(mu)
	self._S = len(self.mus)
	self._approxParameters["S"] = self.S
	if (self.samplingEngine.nsamples <= len(mus) - j - 1
	or not np.allclose(mu, self.samplingEngine.mus[j - len(mus)])):
	self.samplingEngine.nextSample(mu)
	if self._isLastSampleCollinear():
	vbMng(self, "MAIN",
	("Collinearity above tolerance detected. Starting "
	"preemptive greedy loop termination."), 3)
	self._collinearityFlag = 1
	errorEstTest = np.empty(len(self.muTest))
	errorEstTest[:] = np.nan
	return errorEstTest, [-1], np.nan, np.nan
	errorEstTest, muidx, maxErrorEst = self.errorEstimator(self.muTest,
	True)
	if plotEst == "ALL":
	self.plotEstimator(errorEstTest, muidx, maxErrorEst)
	return errorEstTest, muidx, maxErrorEst, self.muTest[muidx]

	def _preliminaryTraining(self):
	"""Initialize starting snapshots of solution map."""
	RROMPyAssert(self._mode, message = "Cannot start greedy algorithm.")
	if self.samplingEngine.nsamples > 0: return
	self.resetSamples()
	self.computeScaleFactor()
	self.samplingEngine.scaleFactor = self.scaleFactorDer
	self.mus = self.samplerTrainSet.generatePoints(self.S)
	while len(self.mus) > self.S: self.mus.pop()
	muTestBase = self.sampler.generatePoints(self.nTestPoints, False)
	idxPop = pruneSamples(self.mapParameterList(muTestBase),
	self.mapParameterList(self.mus),
	1e-10 * self.scaleFactor[0])
	muTestBase.pop(idxPop)
	muLast = copy(self.mus[-1])
	self.mus.pop()
	if len(self.mus) > 0:
	vbMng(self, "MAIN",
	("Adding first {} sample point{} at {} to training "
	"set.").format(self.S - 1, "" + "s" * (self.S > 2),
	self.mus), 3)
	self.samplingEngine.iterSample(self.mus)
	self._S = len(self.mus)
	self._approxParameters["S"] = self.S
	self.muTest = emptyParameterList()
	self.muTest.reset((len(muTestBase) + 1, self.mus.shape[1]))
	self.muTest.data[: -1] = muTestBase.data
	self.muTest.data[-1] = muLast.data

	@abstractmethod
	def setupApproxLocal(self) -> int:
	if self.checkComputedApprox(): return -1
	RROMPyAssert(self._mode, message = "Cannot setup approximant.")
	vbMng(self, "INIT", "Setting up local approximant.", 5)
	pass
	vbMng(self, "DEL", "Done setting up local approximant.", 5)
	return 0

	_postGreedyRecover = 1
	def setupApprox(self, plotEst : str = "NONE") -> int:
	"""Compute greedy snapshots of solution map."""
	if self.checkComputedApprox(): return -1
	RROMPyAssert(self._mode, message = "Cannot start greedy algorithm.")
	vbMng(self, "INIT", "Setting up {}.". format(self.name()), 5)
	vbMng(self, "INIT", "Starting computation of snapshots.", 5)
	self._collinearityFlag = 0
	self._preliminaryTraining()
	muidx, self.firstGreedyIter = [len(self.muTest) - 1], True
	errorEstTest, maxErrorEst = [np.inf], np.inf
	max2ErrorEst, trainedModelOld = np.inf, None
	while self.firstGreedyIter or (len(self.muTest) > 0
	and (maxErrorEst is None or max2ErrorEst > self.greedyTol)
	and self.samplingEngine.nsamples < self.maxIter):
	muTestOld, errorEstTestOld = self.muTest, errorEstTest
	muidxOld, maxErrorEstOld = muidx, maxErrorEst
	errorEstTest, muidx, maxErrorEst, mu = self.greedyNextSample(muidx,
	plotEst)
	if maxErrorEst is not None and (np.any(np.isnan(maxErrorEst))
	or np.any(np.isinf(maxErrorEst))):
	if self._collinearityFlag == 0 and not self.firstGreedyIter:
	RROMPyWarning(("Instability in a posteriori "
	"estimator. Starting preemptive greedy "
	"loop termination."))
	self.muTest, errorEstTest = muTestOld, errorEstTestOld
	if self.firstGreedyIter and muidx[0] < 0:
	self.trainedModel = None
	if self._collinearityFlag:
	raise RROMPyException(("Starting sample points too "
	"collinear. Aborting greedy "
	"iterations."))
	raise RROMPyException(("Instability in approximant "
	"computation. Aborting greedy "
	"iterations."))
	self._S = trainedModelOld.data.approxParameters["S"]
	self._approxParameters["S"] = self.S
	while self.samplingEngine.nsamples > self.S:
	self.samplingEngine.popSample()
	while len(self.mus) > self.S: self.mus.pop(-1)
	muidx, maxErrorEst = muidxOld, maxErrorEstOld
	break
	if maxErrorEst is not None:
	max2ErrorEst = np.max(maxErrorEst)
	vbMng(self, "MAIN", ("Uniform testing error estimate "
	"{:.4e}.").format(max2ErrorEst), 5)
	if self.firstGreedyIter:
	trainedModelOld = copy(self.trainedModel)
	else:
	trainedModelOld.data = copy(self.trainedModel.data)
	self.firstGreedyIter = False
	vbMng(self, "DEL", ("Done computing snapshots (final snapshot count: "
	"{}).").format(self.samplingEngine.nsamples), 5)
	if (maxErrorEst is None or np.any(np.isnan(maxErrorEst))
	or np.any(np.isinf(maxErrorEst))):
	while self.samplingEngine.nsamples > self.S:
	self.samplingEngine.popSample()
	while len(self.mus) > self.S: self.mus.pop(-1)
	elif self._postGreedyRecover:
	self._S = self.samplingEngine.nsamples
	while len(self.mus) < self.S:
	self.mus.append(self.samplingEngine.mus[len(self.mus)])
	self.trainedModel = None
	self.setupApproxLocal()
	if plotEst == "LAST":
	self.plotEstimator(errorEstTest, muidx, maxErrorEst)
	vbMng(self, "DEL", "Done setting up approximant.", 5)
	return 0

	def assembleReducedResidualGramian(self, pMat:sampList):
	"""
	Build residual gramian of reduced linear system through projections.
	"""
	if (not hasattr(self.trainedModel.data, "gramian")
	or self.trainedModel.data.gramian is None):
	gramian = self.HFEngine.innerProduct(pMat, pMat, dual = True)
	else:
	Sold = self.trainedModel.data.gramian.shape[0]
	S = len(self.mus)
	if Sold > S:
	gramian = self.trainedModel.data.gramian[: S, : S]
	else:
	idxOld = list(range(Sold))
	idxNew = list(range(Sold, S))
	gramian = np.empty((S, S), dtype = np.complex)
	gramian[: Sold, : Sold] = self.trainedModel.data.gramian
	gramian[: Sold, Sold :] = self.HFEngine.innerProduct(
	pMat(idxNew), pMat(idxOld), dual = True)
	gramian[Sold :, : Sold] = gramian[: Sold, Sold :].T.conj()
	gramian[Sold :, Sold :] = self.HFEngine.innerProduct(
	pMat(idxNew), pMat(idxNew), dual = True)
	self.trainedModel.data.gramian = gramian

	def assembleReducedResidualBlocksbb(self, bs:List[Np1D]):
	"""
	Build blocks (of type bb) of reduced linear system through projections.
	"""
	nbs = len(bs)
	if (not hasattr(self.trainedModel.data, "resbb")
	or self.trainedModel.data.resbb is None):
	resbb = np.empty((nbs, nbs), dtype = np.complex)
	for i in range(nbs):
	Mbi = bs[i]
	resbb[i, i] = self.HFEngine.innerProduct(Mbi, Mbi, dual = True)
	for j in range(i):
	Mbj = bs[j]
	resbb[i, j] = self.HFEngine.innerProduct(Mbj, Mbi,
	dual = True)
	for i in range(nbs):
	for j in range(i + 1, nbs):
	resbb[i, j] = resbb[j, i].conj()
	self.trainedModel.data.resbb = resbb

	def assembleReducedResidualBlocksAb(self, As:List[Np2D], bs:List[Np1D],
	pMat:sampList):
	"""
	Build blocks (of type Ab) of reduced linear system through projections.
	"""
	nAs = len(As)
	nbs = len(bs)
	S = len(self.mus)
	if (not hasattr(self.trainedModel.data, "resAb")
	or self.trainedModel.data.resAb is None):
	if isinstance(pMat, (parameterList, sampleList)): pMat = pMat.data
	resAb = np.empty((nbs, S, nAs), dtype = np.complex)
	for j in range(nAs):
	MAj = dot(As[j], pMat)
	for i in range(nbs):
	Mbi = bs[i]
	resAb[i, :, j] = self.HFEngine.innerProduct(MAj, Mbi,
	dual = True)
	else:
	Sold = self.trainedModel.data.resAb.shape[1]
	if Sold == S: return
	if Sold > S:
	resAb = self.trainedModel.data.resAb[:, : S, :]
	else:
	if isinstance(pMat, (parameterList, sampleList)):
	pMat = pMat.data
	resAb = np.empty((nbs, S, nAs), dtype = np.complex)
	resAb[:, : Sold, :] = self.trainedModel.data.resAb
	for j in range(nAs):
	MAj = dot(As[j], pMat[:, Sold :])
	for i in range(nbs):
	Mbi = bs[i]
	resAb[i, Sold :, j] = self.HFEngine.innerProduct(
	MAj, Mbi, dual = True)
	self.trainedModel.data.resAb = resAb

	def assembleReducedResidualBlocksAA(self, As:List[Np2D], pMat:sampList):
	"""
	Build blocks (of type AA) of reduced linear system through projections.
	"""
	nAs = len(As)
	S = len(self.mus)
	if (not hasattr(self.trainedModel.data, "resAA")
	or self.trainedModel.data.resAA is None):
	if isinstance(pMat, (parameterList, sampleList)): pMat = pMat.data
	resAA = np.empty((S, nAs, S, nAs), dtype = np.complex)
	for i in range(nAs):
	MAi = dot(As[i], pMat)
	resAA[:, i, :, i] = self.HFEngine.innerProduct(MAi, MAi,
	dual = True)
	for j in range(i):
	MAj = dot(As[j], pMat)
	resAA[:, i, :, j] = self.HFEngine.innerProduct(MAj, MAi,
	dual = True)
	for i in range(nAs):
	for j in range(i + 1, nAs):
	resAA[:, i, :, j] = resAA[:, j, :, i].T.conj()
	else:
	Sold = self.trainedModel.data.resAA.shape[0]
	if Sold == S: return
	if Sold > S:
	resAA = self.trainedModel.data.resAA[: S, :, : S, :]
	else:
	if isinstance(pMat, (parameterList, sampleList)):
	pMat = pMat.data
	resAA = np.empty((S, nAs, S, nAs), dtype = np.complex)
	resAA[: Sold, :, : Sold, :] = self.trainedModel.data.resAA
	for i in range(nAs):
	MAi = dot(As[i], pMat)
	resAA[: Sold, i, Sold :, i] = self.HFEngine.innerProduct(
	MAi[:, Sold :], MAi[:, : Sold], dual = True)
	resAA[Sold :, i, : Sold, i] = resAA[: Sold, i,
	Sold :, i].T.conj()
	resAA[Sold :, i, Sold :, i] = self.HFEngine.innerProduct(
	MAi[:, Sold :], MAi[:, Sold :], dual = True)
	for j in range(i):
	MAj = dot(As[j], pMat)
	resAA[: Sold, i, Sold :, j] = (
	self.HFEngine.innerProduct(MAj[:, Sold :],
	MAi[:, : Sold],
	dual = True))
	resAA[Sold :, i, : Sold, j] = (
	self.HFEngine.innerProduct(MAj[:, : Sold],
	MAi[:, Sold :],
	dual = True))
	resAA[Sold :, i, Sold :, j] = (
	self.HFEngine.innerProduct(MAj[:, Sold :],
	MAi[:, Sold :],
	dual = True))
	for i in range(nAs):
	for j in range(i + 1, nAs):
	resAA[: Sold, i, Sold :, j] = (
	resAA[Sold :, j, : Sold, i].T.conj())
	resAA[Sold :, i, : Sold, j] = (
	resAA[: Sold, j, Sold :, i].T.conj())
	resAA[Sold :, i, Sold :, j] = (
	resAA[Sold :, j, Sold :, i].T.conj())
	self.trainedModel.data.resAA = resAA

	def assembleReducedResidualBlocks(self, full : bool = False):
	"""Build affine blocks of affine decomposition of residual."""
	if full:
	checkIfAffine(self.HFEngine, "assemble reduced residual blocks",
	False, self._affine_lvl)
	else:
	checkIfAffine(self.HFEngine, "assemble reduced RHS blocks", True,
	self._affine_lvl)
	self.HFEngine.buildb()
	self.assembleReducedResidualBlocksbb(self.HFEngine.bs)
	if full:
	pMat = self.samplingEngine.projectionMatrix
	self.HFEngine.buildA()
	self.assembleReducedResidualBlocksAb(self.HFEngine.As,
	self.HFEngine.bs, pMat)
	self.assembleReducedResidualBlocksAA(self.HFEngine.As, pMat)

generic_greedy_approximant.pyNo OneTemporaryActions

File Metadata

generic_greedy_approximant.pyView Options

Event Timeline

generic_greedy_approximant.py
No OneTemporary
Actions

generic_greedy_approximant.py
View Options