Diffusion RationalROMPy (develop) Browse

generic_greedy_approximant.py
6091a93ff6f4develop

generic_greedy_approximant.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.standard.generic_standard_approximant \
			import GenericStandardApproximant
			from rrompy.utilities.base.types import (Np1D, Np2D, DictAny, HFEng, Tuple,
			List, normEng, 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.solver import normEngine
			from rrompy.utilities.exception_manager import (RROMPyException, RROMPyAssert,
			RROMPyWarning)
			from rrompy.parameter import checkParameterList, emptyParameterList

			__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 len(badmus) == 0: return mus
			proximity = np.min(localL2Distance(mus.data, badmus.data), axis = 1)
			return np.arange(len(mus))[proximity <= tol]

			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': whether to compute POD of snapshots; defaults to True;
			- '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.;
			- 'interactive': whether to interactively terminate greedy
			algorithm; defaults to False;
			- 'maxIter': maximum number of greedy steps; defaults to 1e2;
			- 'refinementRatio': ratio of test points to be exhausted before
			test set refinement; defaults to 0.2;
			- 'nTestPoints': number of test points; defaults to 5e2;
			- 'trainSetGenerator': training sample points generator; defaults
			to sampler.
			Defaults to empty dict.
			approx_state(optional): Whether to approximate state. Defaults to
			False.
			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': whether to compute POD of snapshots.
			- 'greedyTol': uniform error tolerance for greedy algorithm;
			- 'collinearityTol': collinearity tolerance for greedy algorithm;
			- 'interactive': whether to interactively terminate greedy
			algorithm;
			- 'maxIter': maximum number of greedy steps;
			- 'refinementRatio': ratio of test points to be exhausted before
			test set refinement;
			- 'nTestPoints': number of test points;
			- 'trainSetGenerator': training sample points generator.
			parameterListCritical: Recognized keys of critical approximant
			parameters:
			- 'S': total number of samples current approximant relies upon;
			- 'sampler': sample point generator.
			approx_state: Whether to approximate state.
			verbosity: Verbosity level.
			POD: whether to compute POD of snapshots.
			S: number of test points.
			sampler: Sample point generator.
			greedyTol: Uniform error tolerance for greedy algorithm.
			collinearityTol: Collinearity tolerance for greedy algorithm.
			interactive: whether to interactively terminate greedy algorithm.
			maxIter: maximum number of greedy steps.
			refinementRatio: ratio of training points to be exhausted before
			training set refinement.
			nTestPoints: number of starting training points.
			trainSetGenerator: training sample points generator.
			muBounds: list of bounds for parameter values.
			samplingEngine: Sampling engine.
			estimatorNormEngine: Engine for estimator norm computation.
			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.
			"""

			TOL_INSTABILITY = 1e-6

			def __init__(self, HFEngine:HFEng, mu0 : paramVal = None,
			approxParameters : DictAny = {}, approx_state : bool = False,
			verbosity : int = 10, timestamp : bool = True):
			self._preInit()
			self._addParametersToList(["greedyTol", "collinearityTol",
			"interactive", "maxIter", "refinementRatio",
			"nTestPoints"],
			[1e-2, 0., False, 1e2, .2, 5e2],
			["trainSetGenerator"], ["AUTO"])
			super().__init__(HFEngine = HFEngine, mu0 = mu0,
			approxParameters = approxParameters,
			approx_state = approx_state, verbosity = verbosity,
			timestamp = timestamp)
			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 interactive(self):
			"""Value of interactive."""
			return self._interactive
			@interactive.setter
			def interactive(self, interactive):
			self._interactive = interactive

			@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 refinementRatio(self):
			"""Value of refinementRatio."""
			return self._refinementRatio
			@refinementRatio.setter
			def refinementRatio(self, refinementRatio):
			if refinementRatio <= 0. or refinementRatio > 1.:
			raise RROMPyException(("refinementRatio must be between 0 "
			"(excluded) and 1."))
			if (hasattr(self, "_refinementRatio")
			and self.refinementRatio is not None):
			refinementRatioold = self.refinementRatio
			else:
			refinementRatioold = -1
			self._refinementRatio = refinementRatio
			self._approxParameters["refinementRatio"] = self.refinementRatio
			if refinementRatioold != self.refinementRatio:
			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 trainSetGenerator(self):
			"""Value of trainSetGenerator."""
			return self._trainSetGenerator
			@trainSetGenerator.setter
			def trainSetGenerator(self, trainSetGenerator):
			if (isinstance(trainSetGenerator, (str,))
			and trainSetGenerator.upper() == "AUTO"):
			trainSetGenerator = self.sampler
			if 'generatePoints' not in dir(trainSetGenerator):
			raise RROMPyException("trainSetGenerator type not recognized.")
			if (hasattr(self, '_trainSetGenerator')
			and self.trainSetGenerator not in [None, "AUTO"]):
			trainSetGeneratorOld = self.trainSetGenerator
			self._trainSetGenerator = trainSetGenerator
			self._approxParameters["trainSetGenerator"] = self.trainSetGenerator
			if (not 'trainSetGeneratorOld' in locals()
			or trainSetGeneratorOld != self.trainSetGenerator):
			self.resetSamples()

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

			def initEstimatorNormEngine(self, normEngn : normEng = None):
			"""Initialize estimator norm engine."""
			if (normEngn is not None or not hasattr(self, "estimatorNormEngine")
			or self.estimatorNormEngine is None):
			if normEngn is None:
			if not hasattr(self.HFEngine, "energyNormPartialDualMatrix"):
			self.HFEngine.buildEnergyNormPartialDualForm()
			estimatorEnergyMatrix = (
			self.HFEngine.energyNormPartialDualMatrix)
			else:
			if hasattr(normEngn, "buildEnergyNormPartialDualForm"):
			if not hasattr(normEngn, "energyNormPartialDualMatrix"):
			normEngn.buildEnergyNormPartialDualForm()
			estimatorEnergyMatrix = (
			normEngn.energyNormPartialDualMatrix)
			else:
			estimatorEnergyMatrix = normEngn
			self.estimatorNormEngine = normEngine(estimatorEnergyMatrix)

			def _affineResidualMatricesContraction(self, rb:Np2D, rA : Np2D = None) \
			-> Tuple[Np1D, Np1D, Np1D]:
			self.assembleReducedResidualBlocks(full = True)
			# '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 errorEstimator(self, mus:Np1D) -> Np1D:
			"""Standard residual-based error estimator."""
			self.setupApprox()
			self.HFEngine.buildA()
			self.HFEngine.buildb()
			mus = checkParameterList(mus, self.npar)[0]
			vbMng(self.trainedModel, "INIT",
			"Evaluating error estimator at mu = {}.".format(mus), 10)
			verb = self.trainedModel.verbosity
			self.trainedModel.verbosity = 0
			uApproxRs = self.getApproxReduced(mus)
			muTestEff = mus ** self.HFEngine.rescalingExp
			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.data, 1) * radiusA
			ff, Lf, LL = self._affineResidualMatricesContraction(radiusb, radiusA)
			err = np.abs((LL - 2. * np.real(Lf) + ff) / ff) ** .5
			self.trainedModel.verbosity = verb
			vbMng(self.trainedModel, "DEL", "Done evaluating error estimator", 10)
			return err

			def getMaxErrorEstimator(self, mus:paramList) -> Tuple[Np1D, int, float]:
			"""
			Compute maximum of (and index of maximum of) error estimator over given
			parameters.
			"""
			errorEstTest = self.errorEstimator(mus)
			idxMaxEst = [np.argmax(errorEstTest)]
			return errorEstTest, idxMaxEst, errorEstTest[idxMaxEst]

			def _isLastSampleCollinear(self) -> bool:
			"""Check collinearity of last sample."""
			if self.collinearityTol <= 0.: return False
			if self.POD:
			reff = self.samplingEngine.RPOD[:, -1]
			else:
			RROMPyWarning(("Repeated orthogonalization of the samples for "
			"collinearity check. Consider setting POD to "
			"True."))
			if not hasattr(self, "_PODEngine"):
			from rrompy.sampling.base.pod_engine 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)
			vbMng(self, "MAIN", "Collinearity indicator {:.4e}.".format(cLevel), 5)
			return cLevel < self.collinearityTol

			def greedyNextSample(self, muidx:int, plotEst : bool = False)\
			-> 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 mu in mus:
			vbMng(self, "MAIN",
			("Adding sample point no. {} at {} to training "
			"set.").format(self.samplingEngine.nsamples + 1, mu), 2)
			self.mus.append(mu)
			self._S = len(self.mus)
			self._approxParameters["S"] = self.S
			self.samplingEngine.nextSample(mu)
			if self._isLastSampleCollinear():
			RROMPyWarning("Collinearity above tolerance detected.")
			errorEstTest = np.empty(len(self.muTest))
			errorEstTest[:] = np.nan
			return errorEstTest, [-1], np.nan, np.nan
			errorEstTest, muidx, maxErrorEst = self.getMaxErrorEstimator(
			self.muTest)
			if (plotEst and not np.any(np.isnan(errorEstTest))
			and not np.any(np.isinf(errorEstTest))):
			musre = copy(self.muTest.re.data)
			from matplotlib import pyplot as plt
			plt.figure()
			errCP = copy(errorEstTest)
			while len(musre) > 0:
			if self.npar == 1:
			currIdx = np.arange(len(musre))
			else:
			currIdx = np.where(np.isclose(np.sum(
			np.abs(musre[:, 1 :] - musre[0, 1 :]), 1), 0.))[0]
			plt.semilogy(musre[currIdx, 0], errCP[currIdx], 'k',
			linewidth = 1)
			musre = np.delete(musre, currIdx, 0)
			errCP = np.delete(errCP, currIdx)
			plt.semilogy([self.muTest.re(0, 0), self.muTest.re(-1, 0)],
			[self.greedyTol] * 2, 'r--')
			plt.semilogy(self.mus.re(0),
			2. * self.greedyTol * np.ones(len(self.mus)), '*m')
			plt.semilogy(self.muTest.re(muidx, 0), maxErrorEst, 'xr')
			plt.grid()
			plt.show()
			plt.close()
			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.computeScaleFactor()
			self.resetSamples()
			self.mus = self.trainSetGenerator.generatePoints(self.S)[
			list(range(self.S))]
			muTestBase = self.sampler.generatePoints(self.nTestPoints)
			idxPop = pruneSamples(muTestBase ** self.HFEngine.rescalingExp,
			self.mus ** self.HFEngine.rescalingExp,
			1e-10 * self.scaleFactor[0])
			muTestBase.pop(idxPop)
			muTestBase = muTestBase.sort()
			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), 2)
			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[: -1] = muTestBase.data
			self.muTest[-1] = muLast.data

			def _enrichTestSet(self, nTest:int):
			"""Add extra elements to test set."""
			RROMPyAssert(self._mode, message = "Cannot enrich test set.")
			muTestExtra = self.sampler.generatePoints(2 * nTest)
			muTotal = copy(self.mus)
			muTotal.append(self.muTest)
			idxPop = pruneSamples(muTestExtra ** self.HFEngine.rescalingExp,
			muTotal ** self.HFEngine.rescalingExp,
			1e-10 * self.scaleFactor[0])
			muTestExtra.pop(idxPop)
			muTestNew = np.empty((len(self.muTest) + len(muTestExtra),
			self.muTest.shape[1]), dtype = np.complex)
			muTestNew[: len(self.muTest)] = self.muTest.data
			muTestNew[len(self.muTest) :] = muTestExtra.data
			self.muTest = checkParameterList(muTestNew, self.npar)[0].sort()
			vbMng(self, "MAIN",
			"Enriching test set by {} elements.".format(len(muTestExtra)), 5)

			def greedy(self, plotEst : bool = False):
			"""Compute greedy snapshots of solution map."""
			RROMPyAssert(self._mode, message = "Cannot start greedy algorithm.")
			if self.samplingEngine.nsamples > 0:
			return
			vbMng(self, "INIT", "Starting computation of snapshots.", 2)
			self._preliminaryTraining()
			nTest = self.nTestPoints
			muT0 = copy(self.muTest[-1])
			errorEstTest, muidx, maxErrorEst, mu = self.greedyNextSample(
			[len(self.muTest) - 1], plotEst)
			if np.any(np.isnan(maxErrorEst)):
			RROMPyWarning(("Instability in a posteriori estimator. "
			"Starting preemptive greedy loop termination."))
			self.muTest.append(muT0)
			self.mus.pop(-1)
			self.samplingEngine.popSample()
			self.setupApprox()
			else:
			vbMng(self, "MAIN", ("Uniform testing error estimate "
			"{:.4e}.").format(np.max(maxErrorEst)), 2)
			trainedModelOld = copy(self.trainedModel)
			while (self.samplingEngine.nsamples < self.maxIter
			and np.max(maxErrorEst) > self.greedyTol):
			if (1. - self.refinementRatio) * nTest > len(self.muTest):
			self._enrichTestSet(nTest)
			nTest = len(self.muTest)
			muTestOld, maxErrorEstOld = self.muTest, np.max(maxErrorEst)
			errorEstTest, muidx, maxErrorEst, mu = self.greedyNextSample(
			muidx, plotEst)
			vbMng(self, "MAIN", ("Uniform testing error estimate "
			"{:.4e}.").format(np.max(maxErrorEst)), 2)
			if (np.any(np.isnan(maxErrorEst))
			or np.any(np.isinf(maxErrorEst))
			or maxErrorEstOld < (np.max(maxErrorEst)
			* self.TOL_INSTABILITY)):
			RROMPyWarning(("Instability in a posteriori estimator. "
			"Starting preemptive greedy loop "
			"termination."))
			self.muTest = muTestOld
			self._approxParameters = (
			trainedModelOld.data.approxParameters)
			self._S -= 1
			self.mus.pop(-1)
			self.samplingEngine.popSample()
			self.trainedModel.data = copy(trainedModelOld.data)
			break
			trainedModelOld.data = copy(self.trainedModel.data)
			if (self.interactive and np.max(maxErrorEst) <= self.greedyTol):
			vbMng(self, "MAIN",
			("Required tolerance {} achieved. Want to decrease "
			"greedyTol and continue? "
			"Y/N").format(self.greedyTol), 0, end = "")
			increasemaxIter = input()
			if increasemaxIter.upper() == "Y":
			vbMng(self, "MAIN", "Reducing value of greedyTol...",
			0)
			while np.max(maxErrorEst) <= self._greedyTol:
			self._greedyTol *= .5
			if (self.interactive
			and self.samplingEngine.nsamples >= self.maxIter):
			vbMng(self, "MAIN",
			("Maximum number of iterations {} reached. Want to "
			"increase maxIter and continue? "
			"Y/N").format(self.maxIter), 0, end = "")
			increasemaxIter = input()
			if increasemaxIter.upper() == "Y":
			vbMng(self, "MAIN", "Doubling value of maxIter...", 0)
			self._maxIter *= 2
			vbMng(self, "DEL",
			("Done computing snapshots (final snapshot count: "
			"{}).").format(self.samplingEngine.nsamples), 2)

			def checkComputedApprox(self) -> bool:
			"""
			Check if setup of new approximant is not needed.

			Returns:
			True if new setup is not needed. False otherwise.
			"""
			return (super().checkComputedApprox()
			and len(self.mus) == self.trainedModel.data.projMat.shape[1])

			def assembleReducedResidualGramian(self, pMat:sampList):
			"""
			Build residual gramian of reduced linear system through projections.
			"""
			self.initEstimatorNormEngine()
			if (not hasattr(self.trainedModel.data, "gramian")
			or self.trainedModel.data.gramian is None):
			gramian = self.estimatorNormEngine.innerProduct(pMat, pMat)
			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.estimatorNormEngine.innerProduct(pMat(idxNew),
			pMat(idxOld)))
			gramian[Sold :, : Sold] = gramian[: Sold, Sold :].T.conj()
			gramian[Sold :, Sold :] = (
			self.estimatorNormEngine.innerProduct(pMat(idxNew),
			pMat(idxNew)))
			self.trainedModel.data.gramian = gramian

			def assembleReducedResidualBlocksbb(self, bs:List[Np1D]):
			"""
			Build blocks (of type bb) of reduced linear system through projections.
			"""
			self.initEstimatorNormEngine()
			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.estimatorNormEngine.innerProduct(Mbi, Mbi)
			for j in range(i):
			Mbj = bs[j]
			resbb[i, j] = self.estimatorNormEngine.innerProduct(Mbj,
			Mbi)
			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.
			"""
			self.initEstimatorNormEngine()
			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 not isinstance(pMat, (np.ndarray,)): 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.estimatorNormEngine.innerProduct(MAj,
			Mbi)
			else:
			Sold = self.trainedModel.data.resAb.shape[1]
			if Sold == S: return
			if Sold > S:
			resAb = self.trainedModel.data.resAb[:, : S, :]
			else:
			if not isinstance(pMat, (np.ndarray,)): 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.estimatorNormEngine.innerProduct(MAj, Mbi))
			self.trainedModel.data.resAb = resAb

			def assembleReducedResidualBlocksAA(self, As:List[Np2D], pMat:sampList):
			"""
			Build blocks (of type AA) of reduced linear system through projections.
			"""
			self.initEstimatorNormEngine()
			nAs = len(As)
			S = len(self.mus)
			if (not hasattr(self.trainedModel.data, "resAA")
			or self.trainedModel.data.resAA is None):
			if not isinstance(pMat, (np.ndarray,)): 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.estimatorNormEngine.innerProduct(MAi, MAi))
			for j in range(i):
			MAj = dot(As[j], pMat)
			resAA[:, i, :, j] = (
			self.estimatorNormEngine.innerProduct(MAj, MAi))
			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 not isinstance(pMat, (np.ndarray,)): 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.estimatorNormEngine.innerProduct(MAi[:, Sold :],
			MAi[:, : Sold]))
			resAA[Sold :, i, : Sold, i] = resAA[: Sold, i,
			Sold :, i].T.conj()
			resAA[Sold :, i, Sold :, i] = (
			self.estimatorNormEngine.innerProduct(MAi[:, Sold :],
			MAi[:, Sold :]))
			for j in range(i):
			MAj = dot(As[j], pMat)
			resAA[: Sold, i, Sold :, j] = (
			self.estimatorNormEngine.innerProduct(MAj[:, Sold :],
			MAi[:, : Sold]))
			resAA[Sold :, i, : Sold, j] = (
			self.estimatorNormEngine.innerProduct(MAj[:, : Sold],
			MAi[:, Sold :]))
			resAA[Sold :, i, Sold :, j] = (
			self.estimatorNormEngine.innerProduct(MAj[:, Sold :],
			MAi[:, Sold :]))
			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."""
			self.HFEngine.buildb()
			self.assembleReducedResidualBlocksbb(self.HFEngine.bs)
			if full:
			pMat = self.samplingEngine.samples
			self.HFEngine.buildA()
			self.assembleReducedResidualBlocksAb(self.HFEngine.As,
			self.HFEngine.bs, pMat)
			self.assembleReducedResidualBlocksAA(self.HFEngine.As, pMat)

RROMPy/rrompy/reduction_methods/greedy/generic_greedy_approximant.py6091a93ff6f4develop

generic_greedy_approximant.pyView Options

RROMPy/rrompy/reduction_methods/greedy/generic_greedy_approximant.py
6091a93ff6f4develop

generic_greedy_approximant.py
View Options