generic_approximant_lagrange_greedy.py
No OneTemporary
Actions

Subscribers

None

File Metadata

Created: Fri, May 3, 18:00

generic_approximant_lagrange_greedy.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/>.
	#

	import numpy as np
	from rrompy.reduction_methods.base.generic_approximant import (
	GenericApproximant)
	from rrompy.reduction_methods.lagrange.generic_approximant_lagrange import (
	GenericApproximantLagrange)
	from rrompy.utilities.base.types import DictAny, HFEng, Tuple, List
	from rrompy.utilities.base import purgeDict, verbosityDepth
	from rrompy.utilities.warning_manager import warn

	__all__ = ['GenericApproximantLagrangeGreedy']

	class GenericApproximantLagrangeGreedy(GenericApproximantLagrange):
	"""
	ROM greedy Lagrange 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]];
	- '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 training points to be exhausted
	before training set refinement; defaults to 0.2;
	- 'nTrainingPoints': number of training points; defaults to
	maxIter / refinementRatio;
	- 'nTestPoints': number of starting test points; defaults to 1;
	- 'trainingSetGenerator': training sample points generator;
	defaults to uniform sampler within muBounds;
	- 'robustTol': tolerance for robust Pade' denominator management;
	defaults to 0.
	Defaults to empty dict.
	homogeneized: Whether to homogeneize Dirichlet BCs. Defaults to False.
	verbosity(optional): Verbosity level. Defaults to 10.

	Attributes:
	HFEngine: HF problem solver.
	mu0: Default parameter.
	mus: Array of snapshot parameters.
	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;
	- '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;
	- 'nTrainingPoints': number of training points;
	- 'nTestPoints': number of starting test points;
	- 'trainingSetGenerator': training sample points generator.
	- 'robustTol': tolerance for robust Pade' denominator management.
	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.
	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.
	nTrainingPoints: number of training points.
	nTestPoints: number of starting test points.
	trainingSetGenerator: training sample points generator.
	robustTol: tolerance for robust Pade' denominator management.
	samplingEngine: Sampling engine.
	uHF: High fidelity solution with wavenumber lastSolvedHF as numpy
	complex vector.
	lastSolvedHF: Wavenumber corresponding to last computed high fidelity
	solution.
	uApp: Last evaluated approximant as numpy complex vector.
	lastApproxParameters: List of parameters corresponding to last
	computed approximant.
	"""

	TOL_INSTABILITY = 1e-6

	def __init__(self, HFEngine:HFEng, mu0 : complex = 0.,
	approxParameters : DictAny = {}, homogeneized : bool = False,
	verbosity : int = 10):
	self._preInit()
	self._addParametersToList(["muBounds","greedyTol", "interactive",
	"maxIter", "refinementRatio",
	"nTrainingPoints", "nTestPoints",
	"trainingSetGenerator"])
	super(GenericApproximantLagrange, self).__init__(
	HFEngine = HFEngine, mu0 = mu0,
	approxParameters = approxParameters,
	homogeneized = homogeneized,
	verbosity = verbosity)
	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,
	["muBounds","greedyTol",
	"interactive", "maxIter",
	"refinementRatio", "nTrainingPoints",
	"nTestPoints",
	"trainingSetGenerator"],
	True, True, baselevel = 1)
	GenericApproximant.approxParameters.fset(self, approxParametersCopy)
	keyList = list(approxParameters.keys())
	if "muBounds" in keyList:
	self.muBounds = approxParameters["muBounds"]
	elif hasattr(self, "muBounds"):
	self.muBounds = self.muBounds
	else:
	self.muBounds = [[0.], [1.]]
	if "greedyTol" in keyList:
	self.greedyTol = approxParameters["greedyTol"]
	elif hasattr(self, "greedyTol"):
	self.greedyTol = self.greedyTol
	else:
	self.greedyTol = 1e-2
	if "interactive" in keyList:
	self.interactive = approxParameters["interactive"]
	elif hasattr(self, "interactive"):
	self.interactive = self.interactive
	else:
	self.interactive = False
	if "maxIter" in keyList:
	self.maxIter = approxParameters["maxIter"]
	elif hasattr(self, "maxIter"):
	self.maxIter = self.maxIter
	else:
	self.maxIter = 1e2
	if "refinementRatio" in keyList:
	self.refinementRatio = approxParameters["refinementRatio"]
	elif hasattr(self, "refinementRatio"):
	self.refinementRatio = self.refinementRatio
	else:
	self.refinementRatio = 0.2
	if "nTrainingPoints" in keyList:
	self.nTrainingPoints = approxParameters["nTrainingPoints"]
	elif hasattr(self, "nTrainingPoints"):
	self.nTrainingPoints = self.nTrainingPoints
	else:
	self.nTrainingPoints = np.int(np.ceil(self.maxIter
	/ self.refinementRatio))
	if "nTestPoints" in keyList:
	self.nTestPoints = approxParameters["nTestPoints"]
	elif hasattr(self, "nTestPoints"):
	self.nTestPoints = self.nTestPoints
	else:
	self.nTestPoints = 1
	if "trainingSetGenerator" in keyList:
	self.trainingSetGenerator = (
	approxParameters["trainingSetGenerator"])
	elif hasattr(self, "trainingSetGenerator"):
	self.trainingSetGenerator = self.trainingSetGenerator
	else:
	from rrompy.utilities.parameter_sampling import QuadratureSampler
	self.trainingSetGenerator = QuadratureSampler(self.muBounds,
	"UNIFORM")
	del QuadratureSampler

	@property
	def S(self):
	"""Value of S."""
	return self._S
	@S.setter
	def S(self, S):
	self._S = S

	@property
	def mus(self):
	"""Value of mus."""
	return self._mus
	@mus.setter
	def mus(self, mus):
	self._mus = np.array(mus, dtype = np.complex)

	@property
	def muBounds(self):
	"""Value of muBounds."""
	return self._muBounds
	@muBounds.setter
	def muBounds(self, muBounds):
	if len(muBounds) != 2:
	raise Exception("2 limits must be specified.")
	try:
	muBounds = muBounds.tolist()
	except:
	muBounds = list(muBounds)
	for j in range(2):
	try:
	len(muBounds[j])
	except:
	muBounds[j] = np.array([muBounds[j]])
	if len(muBounds[0]) != len(muBounds[1]):
	raise Exception("The bounds must have the same length.")
	self._muBounds = muBounds

	@property
	def greedyTol(self):
	"""Value of greedyTol."""
	return self._greedyTol
	@greedyTol.setter
	def greedyTol(self, greedyTol):
	if greedyTol < 0:
	raise ArithmeticError("greedyTol must be non-negative.")
	if hasattr(self, "greedyTol"): 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 ArithmeticError("maxIter must be positive.")
	if hasattr(self, "maxIter"): 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 ArithmeticError(("refinementRatio must be between 0 "
	"(excluded) and 1."))
	if hasattr(self, "refinementRatio"):
	refinementRatioold = self.refinementRatio
	else: refinementRatioold = -1
	self._refinementRatio = refinementRatio
	self._approxParameters["refinementRatio"] = self.refinementRatio
	if refinementRatioold != self.refinementRatio:
	self.resetSamples()

	@property
	def nTrainingPoints(self):
	"""Value of nTrainingPoints."""
	return self._nTrainingPoints
	@nTrainingPoints.setter
	def nTrainingPoints(self, nTrainingPoints):
	if nTrainingPoints <= 1:
	raise ArithmeticError("nTrainingPoints must be greater than 1.")
	if not np.isclose(nTrainingPoints, np.int(nTrainingPoints)):
	raise ArithmeticError("nTrainingPoints must be an integer.")
	nTrainingPoints = np.int(nTrainingPoints)
	if hasattr(self, "nTrainingPoints"):
	nTrainingPointsold = self.nTrainingPoints
	else: nTrainingPointsold = -1
	self._nTrainingPoints = nTrainingPoints
	self._approxParameters["nTrainingPoints"] = self.nTrainingPoints
	if nTrainingPointsold != self.nTrainingPoints:
	self.resetSamples()

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

	@property
	def trainingSetGenerator(self):
	"""Value of trainingSetGenerator."""
	return self._trainingSetGenerator
	@trainingSetGenerator.setter
	def trainingSetGenerator(self, trainingSetGenerator):
	if 'generatePoints' not in dir(trainingSetGenerator):
	raise Exception("trainingSetGenerator type not recognized.")
	if hasattr(self, '_trainingSetGenerator'):
	trainingSetGeneratorOld = self.trainingSetGenerator
	self._trainingSetGenerator = trainingSetGenerator
	self._approxParameters["trainingSetGenerator"] = (
	self.trainingSetGenerator)
	if (not 'trainingSetGeneratorOld' in locals()
	or trainingSetGeneratorOld != self.trainingSetGenerator):
	self.resetSamples()

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

	def initEstNormer(self):
	"""Initialize estimator norm engine."""
	if not hasattr(self, "estNormer"):
	from rrompy.hfengines.base import ProblemEngineBase as PEB
	self.estNormer = PEB() # L2 norm
	self.estNormer.V = self.HFEngine.V
	self.estNormer.buildEnergyNormForm()

	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.estNormer.norm(RHS)
	for j in range(nmus):
	res = self.getRes(mus[j], homogeneized = self.homogeneized)
	err[j] = self.estNormer.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.estNormer.norm(res) / self.estNormer.norm(RHS)
	return np.abs(err)

	def getMaxErrorEstimator(self) -> Tuple[float, int]:
	"""
	Compute maximum of (and index of maximum of) error estimator over
	training set.
	"""
	errorEstTrain = self.errorEstimator(self.muTrain)
	idxMaxEst = np.argmax(errorEstTrain)
	maxEst = errorEstTrain[idxMaxEst]
	return maxEst, idxMaxEst

	def greedyNextSample(self, muidx:int, plotEst : bool = False):
	"""Compute next greedy snapshot of solution map."""
	mu = self.muTrain[muidx]
	if self.verbosity >= 2:
	verbosityDepth("MAIN", ("Adding {}-th sample point at {} to "
	"training set.").format(
	self.samplingEngine.nsamples + 1, mu))
	self.mus = np.append(self.mus, mu)
	idxs = np.arange(len(self.muTrain))
	mask = np.ones_like(idxs, dtype = bool)
	mask[muidx] = False
	idxs = idxs[mask]
	self.muTrain = self.muTrain[idxs]
	self.samplingEngine.nextSample(mu,
	homogeneized = self.homogeneized)
	errorEstTrain = self.errorEstimator(self.muTrain)
	muidx = np.argmax(errorEstTrain)
	maxErrorEst = errorEstTrain[muidx]
	mu = self.muTrain[muidx]
	if plotEst and not np.all(np.isinf(errorEstTrain)):
	from matplotlib import pyplot as plt
	plt.figure()
	plt.semilogy(np.real(self.muTrain), errorEstTrain, 'k')
	plt.semilogy(np.real(self.muTrain),
	self.greedyTol * np.ones(len(self.muTrain)), 'r--')
	plt.semilogy(np.real(self.mus),
	2. * self.greedyTol * np.ones(len(self.mus)), '*m')
	plt.semilogy(np.real(mu), maxErrorEst, 'xr')
	plt.grid()
	plt.show()
	plt.close()
	return errorEstTrain, muidx, maxErrorEst, mu

	def greedy(self, plotEst : bool = False):
	"""Compute greedy snapshots of solution map."""
	if self.samplingEngine.samples is None:
	if self.verbosity >= 2:
	verbosityDepth("INIT", "Starting computation of snapshots.")
	self.resetSamples()
	self.mus, _ = self.trainingSetGenerator.generatePoints(
	self.nTestPoints)
	self.mus = np.array([x[0] for x in self.mus], dtype = np.complex)
	nTrain = self.nTrainingPoints
	muTrainBase, _ = self.trainingSetGenerator.generatePoints(nTrain)
	self.muTrain = np.empty(len(muTrainBase) + 1, dtype = np.complex)
	j = 0
	for mu in muTrainBase:
	if not np.any(np.isclose(self.mus, mu)):
	self.muTrain[j] = mu[0]
	j += 1
	self.muTrain[j] = self.mus[-1]
	self.muTrain = self.muTrain[: j + 1]
	self.mus = self.mus[:-1]
	for j in range(len(self.mus)):
	if self.verbosity >= 2:
	verbosityDepth("MAIN", ("Adding {}-th sample point at {} "
	"to training set.").format(
	self.samplingEngine.nsamples + 1,
	self.mus[j]))
	self.samplingEngine.nextSample(self.mus[j],
	homogeneized = self.homogeneized)
	errorEstTrain, muidx, maxErrorEst, mu = self.greedyNextSample(-1,
	plotEst)
	if self.verbosity >= 2:
	verbosityDepth("MAIN", ("Uniform error estimate {:.4e}.")\
	.format(maxErrorEst))
	while (self.samplingEngine.nsamples < self.maxIter
	and maxErrorEst > self.greedyTol):
	if (1. - self.refinementRatio) * nTrain > len(self.muTrain):
	muTrainExtra, _ = self.trainingSetGenerator.refine(nTrain)
	self.muTrain = np.sort(np.append(self.muTrain,
	muTrainExtra))
	nTrain += len(muTrainExtra)
	if self.verbosity >= 5:
	verbosityDepth("MAIN", ("Enriching training set by {} "
	"elements.").format(
	len(muTrainExtra)))
	muTrainOld, maxErrorEstOld = self.muTrain, maxErrorEst
	errorEstTrain, muidx, maxErrorEst, mu = self.greedyNextSample(
	muidx, plotEst)
	if self.verbosity >= 2:
	verbosityDepth("MAIN", ("Uniform error estimate {:.4e}.")\
	.format(maxErrorEst))
	if (np.isnan(maxErrorEst) or np.isinf(maxErrorEst)
	or maxErrorEstOld < maxErrorEst * self.TOL_INSTABILITY):
	warn(("Instability in a posteriori estimator. Starting "
	"preemptive greedy loop termination."))
	maxErrorEst = maxErrorEstOld
	self.muTrain = muTrainOld
	self.mus = self.mus[:-1]
	self.samplingEngine.popSample()
	self.setupApprox()
	break
	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))
	increasemaxIter = input()
	if increasemaxIter.upper() == "Y":
	verbosityDepth("MAIN", "Doubling value of maxIter...")
	self.maxIter *= 2
	if (self.interactive and maxErrorEst <= self.greedyTol):
	verbosityDepth("MAIN", ("Required tolerance {} achieved. "
	"Want to decrease greedyTol and "
	"continue? Y/N").format(
	self.greedyTol))
	increasemaxIter = input()
	if increasemaxIter.upper() == "Y":
	verbosityDepth("MAIN", "Halving value of greedyTol...")
	self.greedyTol *= .5
	if self.verbosity >= 2:
	verbosityDepth("DEL", ("Done computing snapshots (final "
	"snapshot count: {}).").format(
	self.samplingEngine.nsamples))

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

	Returns:
	True if new setup is not needed. False otherwise.
	"""
	return (hasattr(self, "_S") and self.S == len(self.mus)
	and super().checkComputedApprox())

	def computeScaleFactor(self):
	"""Compute parameter rescaling factor."""
	self.scaleFactor = .5 * np.abs(self.HFEngine.rescaling(
	self.trainingSetGenerator.lims[0][0])
	- self.HFEngine.rescaling(
	self.trainingSetGenerator.lims[1][0]))

generic_approximant_lagrange_greedy.pyNo OneTemporaryActions

File Metadata

generic_approximant_lagrange_greedy.pyView Options

Event Timeline

generic_approximant_lagrange_greedy.py
No OneTemporary
Actions

generic_approximant_lagrange_greedy.py
View Options