generic_distributed_greedy_approximant.py
No OneTemporary
Actions

Subscribers

None

File Metadata

Created: Fri, May 3, 08:45

generic_distributed_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.distributed.generic_distributed_approximant \
	import GenericDistributedApproximant
	from rrompy.utilities.base.types import (Np1D, Np2D, DictAny, HFEng, Tuple,
	List, normEng, paramList, sampList)
	from rrompy.utilities.base import purgeDict, verbosityDepth
	from rrompy.solver import normEngine
	from rrompy.utilities.exception_manager import (RROMPyException, RROMPyAssert,
	RROMPyWarning)
	from rrompy.parameter import checkParameterList, emptyParameterList

	__all__ = ['GenericDistributedGreedyApproximant']


	def pruneSamples(mus:paramList, badmus:paramList,
	tol : float = 1e-8) -> paramList:
	"""Remove from mus all the elements which are too close to badmus."""
	if len(badmus) == 0: return mus
	musNp = np.array(mus(0))
	badmus = np.array(badmus(0))
	proximity = np.min(np.abs(musNp.reshape(-1, 1)
	- np.tile(badmus.reshape(1, -1), [len(mus), 1])),
	axis = 1).flatten()
	idxPop = np.arange(len(mus))[proximity <= tol]
	for i, j in enumerate(idxPop):
	mus.pop(j - i)
	return mus

	class GenericDistributedGreedyApproximant(GenericDistributedApproximant):
	"""
	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;
	- 'muBounds': list of bounds for parameter values; defaults to
	[0, 1];
	- 'S': number of starting training points; defaults to 2;
	- 'sampler': sample point generator; defaults to uniform sampler on
	muBounds;
	- 'greedyTol': uniform error tolerance for greedy algorithm;
	defaults to 1e-2;
	- '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 maxIter /
	refinementRatio;
	- 'trainSetGenerator': training sample points generator; defaults
	to Chebyshev sampler within muBounds.
	Defaults to empty dict.
	homogeneized(optional): 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.
	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;
	- 'muBounds': list of bounds for parameter values;
	- 'S': number of starting training points;
	- 'sampler': sample point generator;
	- 'greedyTol': uniform error 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.
	extraApproxParameters: List of approxParameters keys in addition to
	mother class's.
	POD: whether to compute POD of snapshots.
	muBounds: list of bounds for parameter values.
	S: number of test points.
	sampler: Sample point generator.
	greedyTol: uniform error tolerance for 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.
	robustTol: tolerance for robust Pade' denominator management.
	samplingEngine: Sampling engine.
	estimatorNormEngine: Engine for estimator norm computation.
	uHF: High fidelity solution with wavenumber lastSolvedHF as numpy
	complex vector.
	lastSolvedHF: Wavenumber corresponding to last computed high fidelity
	solution.
	uApp: Last evaluated approximant as numpy complex vector.
	"""

	TOL_INSTABILITY = 1e-6

	def __init__(self, HFEngine:HFEng, mu0 : complex = 0.,
	approxParameters : DictAny = {}, homogeneized : bool = False,
	verbosity : int = 10, timestamp : bool = True):
	self._preInit()
	self._addParametersToList(["greedyTol", "interactive", "maxIter",
	"refinementRatio", "nTestPoints",
	"trainSetGenerator"])
	super().__init__(HFEngine = HFEngine, mu0 = mu0,
	approxParameters = approxParameters,
	homogeneized = homogeneized, verbosity = verbosity,
	timestamp = timestamp)
	self._postInit()

	@property
	def approxParameters(self):
	"""Value of approximant parameters. Its assignment may change S."""
	return self._approxParameters
	@approxParameters.setter
	def approxParameters(self, approxParams):
	approxParameters = purgeDict(approxParams, self.parameterList,
	dictname = self.name() + ".approxParameters",
	baselevel = 1)
	approxParametersCopy = purgeDict(approxParameters,
	["greedyTol", "interactive",
	"maxIter", "refinementRatio",
	"nTestPoints", "trainSetGenerator"],
	True, True, baselevel = 1)
	GenericDistributedApproximant.approxParameters.fset(self,
	approxParametersCopy)
	keyList = list(approxParameters.keys())
	if "greedyTol" in keyList:
	self.greedyTol = approxParameters["greedyTol"]
	elif not hasattr(self, "_greedyTol") or self.greedyTol is None:
	self.greedyTol = 1e-2
	if "interactive" in keyList:
	self.interactive = approxParameters["interactive"]
	elif not hasattr(self, "interactive") or self.interactive is None:
	self.interactive = False
	if "maxIter" in keyList:
	self.maxIter = approxParameters["maxIter"]
	elif not hasattr(self, "_maxIter") or self.maxIter is None:
	self.maxIter = 1e2
	if "refinementRatio" in keyList:
	self.refinementRatio = approxParameters["refinementRatio"]
	elif (not hasattr(self, "_refinementRatio")
	or self.refinementRatio is None):
	self.refinementRatio = 0.2
	if "nTestPoints" in keyList:
	self.nTestPoints = approxParameters["nTestPoints"]
	elif (not hasattr(self, "_nTestPoints")
	or self.nTestPoints is None):
	self.nTestPoints = np.int(np.ceil(self.maxIter
	/ self.refinementRatio))
	if "trainSetGenerator" in keyList:
	self.trainSetGenerator = approxParameters["trainSetGenerator"]
	elif (not hasattr(self, "_trainSetGenerator")
	or self.trainSetGenerator is None):
	from rrompy.parameter.parameter_sampling import QuadratureSampler
	self.trainSetGenerator = QuadratureSampler(self.muBounds,
	"CHEBYSHEV")
	del QuadratureSampler

	@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 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 'generatePoints' not in dir(trainSetGenerator):
	raise RROMPyException("trainSetGenerator type not recognized.")
	if (hasattr(self, '_trainSetGenerator')
	and self.trainSetGenerator is not None):
	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, "energyNormMatrix"):
	self.HFEngine.buildEnergyNormForm()
	estimatorEnergyMatrix = self.HFEngine.energyNormMatrix
	else:
	if hasattr(normEngn, "buildEnergyNormForm"):
	if not hasattr(normEngn, "energyNormMatrix"):
	normEngn.buildEnergyNormForm()
	estimatorEnergyMatrix = normEngn.energyNormMatrix
	else:
	estimatorEnergyMatrix = normEngn
	self.estimatorNormEngine = normEngine(estimatorEnergyMatrix)

	def errorEstimator(self, mus:List[np.complex]) -> List[np.complex]:
	"""
	Standard residual-based error estimator with explicit residual
	computation.
	"""
	self.setupApprox()
	nmus = len(mus)
	err = np.empty(nmus)
	if self.HFEngine.nbs == 1:
	RHS = self.getRHS(mus[0], homogeneized = self.homogeneized)
	RHSNorm = self.estimatorNormEngine.norm(RHS)
	for j in range(nmus):
	res = self.getRes(mus[j], homogeneized = self.homogeneized)
	err[j] = self.estimatorNormEngine.norm(res) / RHSNorm
	else:
	for j in range(nmus):
	res = self.getRes(mus[j], homogeneized = self.homogeneized)
	RHS = self.getRHS(mus[j], homogeneized = self.homogeneized)
	err[j] = (self.estimatorNormEngine.norm(res)
	/ self.estimatorNormEngine.norm(RHS))
	return np.abs(err)

	def getMaxErrorEstimator(self, mus:paramList,
	plot : bool = False) -> Tuple[Np1D, int, float]:
	"""
	Compute maximum of (and index of maximum of) error estimator over given
	parameters.
	"""
	errorEstTest = self.errorEstimator(mus(0))
	idxMaxEst = np.argmax(errorEstTest)
	maxEst = errorEstTest[idxMaxEst]
	if plot and not np.all(np.isinf(errorEstTest)):
	musre = mus.re(0)
	from matplotlib import pyplot as plt
	plt.figure()
	plt.semilogy(musre, errorEstTest, 'k')
	plt.semilogy([musre[0], musre[-1]], [self.greedyTol] * 2, 'r--')
	plt.semilogy(self.mus(0),
	2. * self.greedyTol * np.ones(len(self.mus)), '*m')
	plt.semilogy(musre[idxMaxEst], maxEst, 'xr')
	plt.grid()
	plt.show()
	plt.close()
	return errorEstTest, idxMaxEst, maxEst

	def greedyNextSample(self, muidx:int, plotEst : bool = False)\
	-> Tuple[Np1D, int, float, complex]:
	"""Compute next greedy snapshot of solution map."""
	RROMPyAssert(self._mode, message = "Cannot add greedy sample.")
	mu = copy(self.muTest[muidx])
	self.muTest.pop(muidx)
	if self.verbosity >= 2:
	verbosityDepth("MAIN", ("Adding {}-th sample point at {} to "
	"training set.").format(
	self.samplingEngine.nsamples + 1, mu),
	timestamp = self.timestamp)
	self.mus.append(mu)
	self.samplingEngine.nextSample(mu, homogeneized = self.homogeneized)
	errorEstTest, muidx, maxErrorEst = self.getMaxErrorEstimator(
	self.muTest, plotEst)
	return errorEstTest, muidx, maxErrorEst, self.muTest[muidx]

	def _preliminaryTraining(self):
	"""Initialize starting snapshots of solution map."""
	RROMPyAssert(self._mode, message = "Cannot start greedy algorithm.")
	self.computeScaleFactor()
	if self.samplingEngine.nsamples > 0:
	return
	if self.verbosity >= 2:
	verbosityDepth("INIT", "Starting computation of snapshots.",
	timestamp = self.timestamp)
	self.resetSamples()
	self.mus, _ = self.trainSetGenerator.generatePoints(self.S)
	muLast = copy(self.mus[-1])
	self.mus.pop()
	muTestBase, _ = self.sampler.generatePoints(self.nTestPoints)
	if len(self.mus) > 1:
	if self.verbosity >= 2:
	verbosityDepth("MAIN",
	("Adding first {} samples point at {} to "
	"training set.").format(self.S - 1, self.mus),
	timestamp = self.timestamp)
	self.samplingEngine.iterSample(self.mus,
	homogeneized = self.homogeneized)
	muTestBase = np.sort(pruneSamples(muTestBase, self.mus,
	1e-10 * self.scaleFactor))
	self.muTest = emptyParameterList()
	self.muTest.reset(len(muTestBase) + 1)
	self.muTest[: -1] = muTestBase
	self.muTest[-1] = muLast

	def _enrichTestSet(self, nTest:int):
	"""Add extra elements to test set."""

	muTestExtra, _ = self.sampler.generatePoints(2 * nTest)
	muTotal = copy(self.mus)
	muTotal.append(self.muTest)
	muTestExtra = pruneSamples(muTestExtra, muTotal,
	1e-10 * self.scaleFactor)
	muTestNew = np.empty(len(self.muTest) + len(muTestExtra),
	dtype = np.complex)
	muTestNew[: len(self.muTest)] = self.muTest(0)
	muTestNew[len(self.muTest) :] = muTestExtra(0)
	self.muTest = checkParameterList(np.sort(muTestNew))
	if self.verbosity >= 5:
	verbosityDepth("MAIN", "Enriching test set by {} elements.".format(
	len(muTestExtra)),
	timestamp = self.timestamp)

	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
	self._preliminaryTraining()
	nTest = self.nTestPoints
	errorEstTest, muidx, maxErrorEst, mu = self.greedyNextSample(-1,
	plotEst)
	if self.verbosity >= 2:
	verbosityDepth("MAIN", ("Uniform testing error estimate "
	"{:.4e}.").format(maxErrorEst),
	timestamp = self.timestamp)
	trainedModelOld = copy(self.trainedModel)
	while (self.samplingEngine.nsamples < self.maxIter
	and maxErrorEst > self.greedyTol):
	if (1. - self.refinementRatio) * nTest > len(self.muTest):
	self._enrichTestSet(nTest)
	nTest = len(self.muTest)
	muTestOld, maxErrorEstOld = self.muTest, maxErrorEst
	errorEstTest, muidx, maxErrorEst, mu = self.greedyNextSample(
	muidx, plotEst)
	if self.verbosity >= 2:
	verbosityDepth("MAIN", ("Uniform testing error estimate "
	"{:.4e}.").format(maxErrorEst),
	timestamp = self.timestamp)
	if (np.isnan(maxErrorEst) or np.isinf(maxErrorEst)
	or maxErrorEstOld < maxErrorEst * self.TOL_INSTABILITY):
	RROMPyWarning(("Instability in a posteriori estimator. "
	"Starting preemptive greedy loop termination."))
	maxErrorEst = maxErrorEstOld
	self.muTest = muTestOld
	self.mus = self.mus[:-1]
	self.samplingEngine.popSample()
	self.trainedModel.data = copy(trainedModelOld.data)
	break
	trainedModelOld.data = copy(self.trainedModel.data)
	if (self.interactive and maxErrorEst <= self.greedyTol):
	verbosityDepth("MAIN", ("Required tolerance {} achieved. Want "
	"to decrease greedyTol and continue? "
	"Y/N").format(self.greedyTol),
	timestamp = self.timestamp, end = "")
	increasemaxIter = input()
	if increasemaxIter.upper() == "Y":
	verbosityDepth("MAIN", "Reducing value of greedyTol...",
	timestamp = self.timestamp)
	while maxErrorEst <= self._greedyTol:
	self._greedyTol *= .5
	if (self.interactive
	and self.samplingEngine.nsamples >= self.maxIter):
	verbosityDepth("MAIN", ("Maximum number of iterations {} "
	"reached. Want to increase maxIter "
	"and continue? Y/N").format(
	self.maxIter),
	timestamp = self.timestamp, end = "")
	increasemaxIter = input()
	if increasemaxIter.upper() == "Y":
	verbosityDepth("MAIN", "Doubling value of maxIter...",
	timestamp = self.timestamp)
	self._maxIter *= 2
	if self.verbosity >= 2:
	verbosityDepth("DEL", ("Done computing snapshots (final snapshot "
	"count: {}).").format(
	self.samplingEngine.nsamples),
	timestamp = self.timestamp)

	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],
	scaling : float = 1.):
	"""
	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 = scaling ** i * bs[i]
	resbb[i, i] = self.estimatorNormEngine.innerProduct(Mbi, Mbi)
	for j in range(i):
	Mbj = scaling ** j * 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, scaling : float = 1.):
	"""
	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 = scaling ** (j + 1) * As[j].dot(pMat)
	for i in range(nbs):
	Mbi = scaling ** (i + 1) * 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 = scaling ** (j + 1) * As[j].dot(pMat[:, Sold :])
	for i in range(nbs):
	Mbi = scaling ** (i + 1) * bs[i]
	resAb[i, Sold :, j] = (
	self.estimatorNormEngine.innerProduct(MAj, Mbi))
	self.trainedModel.data.resAb = resAb

	def assembleReducedResidualBlocksAA(self, As:List[Np2D], pMat:sampList,
	scaling : float = 1.,
	basic : bool = False):
	"""
	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
	if basic:
	MAEnd = scaling ** nAs * As[-1].dot(pMat)
	resAA = self.estimatorNormEngine.innerProduct(MAEnd, MAEnd)
	else:
	resAA = np.empty((S, nAs, S, nAs), dtype = np.complex)
	for i in range(nAs):
	MAi = scaling ** (i + 1) * As[i].dot(pMat)
	resAA[:, i, :, i] = (
	self.estimatorNormEngine.innerProduct(MAi, MAi))
	for j in range(i):
	MAj = scaling ** (j + 1) * As[j].dot(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:
	if basic:
	resAA = self.trainedModel.data.resAA[: S, : S]
	else:
	resAA = self.trainedModel.data.resAA[: S, :, : S, :]
	else:
	if not isinstance(pMat, (np.ndarray,)): pMat = pMat.data
	if basic:
	resAA = np.empty((S, S), dtype = np.complex)
	resAA[: Sold, : Sold] = self.trainedModel.data.resAA
	MAi = scaling ** nAs * As[-1].dot(pMat)
	resAA[: Sold, Sold :] = (
	self.estimatorNormEngine.innerProduct(MAi[:, Sold :],
	MAi[:, : Sold]))
	resAA[Sold :, : Sold] = resAA[: Sold, Sold :].T.conj()
	resAA[Sold :, Sold :] = (
	self.estimatorNormEngine.innerProduct(MAi[:, Sold :],
	MAi[:, Sold :]))
	else:
	resAA = np.empty((S, nAs, S, nAs), dtype = np.complex)
	resAA[: Sold, :, : Sold, :] = self.trainedModel.data.resAA
	for i in range(nAs):
	MAi = scaling ** (i + 1) * As[i].dot(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 = scaling ** (j + 1) * As[j].dot(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

generic_distributed_greedy_approximant.pyNo OneTemporaryActions

File Metadata

generic_distributed_greedy_approximant.pyView Options

Event Timeline

generic_distributed_greedy_approximant.py
No OneTemporary
Actions

generic_distributed_greedy_approximant.py
View Options