generic_greedy_approximant.py
9dc9eb5a8081master

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 = copy(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

			def addSamplePoints(self, mus:paramList):
			"""Add sample points to reduced model."""
			raise RROMPyException("Cannot add samples to greedy reduced model.")

			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)
			else:
			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)

RROMPy/rrompy/reduction_methods/standard/greedy/generic_greedy_approximant.py9dc9eb5a8081master

generic_greedy_approximant.pyView Options

RROMPy/rrompy/reduction_methods/standard/greedy/generic_greedy_approximant.py
9dc9eb5a8081master

generic_greedy_approximant.py
View Options