Diffusion RationalROMPy (feature_greedy_coarsening) Browse

generic_approximant_lagrange_greedy.py
b68d96516204feature_greedy_coarseni…

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;
			- '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.

			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.
			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", "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", "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 "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
			# class L2InverseNormer(PEB):
			# def innerProduct(self, u:Np1D, v:Np1D) -> float:
			# if not hasattr(self, "energyNormMatrix"):
			# self.buildEnergyNormForm()
			# return v.conj().T.dot(self.energyNormMatrix.solve(u))

			# def buildEnergyNormForm(self):
			# super().buildEnergyNormForm()
			# from scipy.sparse import csc_matrix, linalg as sla
			# Mass = csc_matrix(self.energyNormMatrix,
			# dtype = np.complex)
			# self.energyNormMatrix = sla.spilu(Mass)
			# self.estNormer = L2InverseNormer() # inverse L2 norm
			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 test "
			"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 test 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
			or self.samplingEngine.nsamples < self.nTestPoints)):
			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 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.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 computeRescaleParameter(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]))

RROMPy/rrompy/reduction_methods/lagrange_greedy/generic_approximant_lagrange_greedy.pyb68d96516204feature_greedy_coarseni…

generic_approximant_lagrange_greedy.pyView Options

RROMPy/rrompy/reduction_methods/lagrange_greedy/generic_approximant_lagrange_greedy.py
b68d96516204feature_greedy_coarseni…

generic_approximant_lagrange_greedy.py
View Options