generic_pivoted_approximant.py
16a7048180acmaster

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

			import numpy as np
			from copy import deepcopy as copy
			from rrompy.reduction_methods.base.generic_approximant import (
			GenericApproximant)
			from rrompy.hfengines.base import MarginalProxyEngine
			from rrompy.utilities.poly_fitting.polynomial import (polybases,
			nextDerivativeIndices,
			PolynomialInterpolator as PI)
			from rrompy.utilities.poly_fitting.radial_basis import (
			RadialBasisInterpolator as RBI)
			from rrompy.utilities.poly_fitting.radial_basis import rbbases
			from rrompy.sampling.pivoted import (SamplingEnginePivoted,
			SamplingEnginePivotedPOD)
			from rrompy.utilities.base.types import (ListAny, DictAny, HFEng, paramVal,
			paramList)
			from rrompy.utilities.base import verbosityManager as vbMng, freepar as fp
			from rrompy.utilities.exception_manager import (RROMPyException, RROMPyAssert,
			RROMPyWarning)
			from rrompy.parameter import emptyParameterList

			__all__ = ['GenericPivotedApproximant']

			class GenericPivotedApproximant(GenericApproximant):
			"""
			ROM pivoted approximant computation for parametric problems (ABSTRACT).

			Args:
			HFEngine: HF problem solver.
			mu0(optional): Default parameter. Defaults to 0.
			directionPivot(optional): Pivot components. 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': total number of pivot samples current approximant relies
			upon;
			- 'samplerPivot': pivot sample point generator;
			- 'SMarginal': total number of marginal samples current approximant
			relies upon;
			- 'samplerMarginal': marginal sample point generator;
			- 'polybasisMarginal': type of polynomial basis for marginal
			interpolation; allowed values include 'MONOMIAL', 'CHEBYSHEV'
			and 'LEGENDRE'; defaults to 'MONOMIAL';
			- 'MMarginal': degree of marginal interpolant; defaults to 0;
			- 'radialBasisMarginal': radial basis family for marginal
			interpolant; defaults to 0, i.e. no radial basis;
			- 'radialBasisWeightsMarginal': radial basis weights for marginal
			interpolant; defaults to 0, i.e. identity;
			- 'interpRcondMarginal': tolerance for marginal interpolation;
			defaults to None.
			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.
			directionPivot: Pivot components.
			mus: Array of snapshot parameters.
			musPivot: Array of pivot snapshot parameters.
			musMarginal: Array of marginal snapshot parameters.
			homogeneized: Whether to homogeneize Dirichlet BCs.
			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;
			- 'polybasisMarginal': type of polynomial basis for marginal
			interpolation;
			- 'MMarginal': degree of marginal interpolant;
			- 'radialBasisMarginal': radial basis family for marginal
			interpolant;
			- 'radialBasisWeightsMarginal': radial basis weights for marginal
			interpolant;
			- 'interpRcondMarginal': tolerance for marginal interpolation.
			parameterListCritical: Recognized keys of critical approximant
			parameters:
			- 'S': total number of pivot samples current approximant relies
			upon;
			- 'samplerPivot': pivot sample point generator;
			- 'SMarginal': total number of marginal samples current approximant
			relies upon;
			- 'samplerMarginal': marginal sample point generator.
			POD: Whether to compute POD of snapshots.
			S: Total number of pivot samples current approximant relies upon.
			samplerPivot: Pivot sample point generator.
			SMarginal: Total number of marginal samples current approximant relies
			upon.
			samplerMarginal: Marginal sample point generator.
			polybasisMarginal: Type of polynomial basis for marginal interpolation.
			MMarginal: Degree of marginal interpolant.
			radialBasisMarginal: Radial basis family for marginal interpolant.
			radialBasisWeightsMarginal: Radial basis weights for marginal
			interpolant.
			interpRcondMarginal: Tolerance for marginal interpolation.
			muBoundsPivot: list of bounds for pivot parameter values.
			muBoundsMarginal: list of bounds for marginal 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, HFEngine:HFEng, mu0 : paramVal = None,
			directionPivot : ListAny = [0],
			approxParameters : DictAny = {}, homogeneized : bool = False,
			verbosity : int = 10, timestamp : bool = True):
			self._preInit()
			from rrompy.parameter.parameter_sampling import QuadratureSampler as QS
			QSBase = QS([[0], [1]], "UNIFORM")
			self._addParametersToList(
			["polybasisMarginal", "MMarginal", "radialBasisMarginal",
			"radialBasisWeightsMarginal", "interpRcondMarginal"],
			["MONOMIAL", 0, 0, 1, -1],
			["samplerPivot", "SMarginal", "samplerMarginal"],
			[QSBase, [1], QSBase])
			del QS
			self._directionPivot = directionPivot
			super().__init__(HFEngine = HFEngine, mu0 = mu0,
			approxParameters = approxParameters,
			homogeneized = homogeneized,
			verbosity = verbosity, timestamp = timestamp)
			self._postInit()

			def setupSampling(self):
			"""Setup sampling engine."""
			RROMPyAssert(self._mode, message = "Cannot setup sampling engine.")
			if not hasattr(self, "_POD") or self._POD is None: return
			if self.POD:
			SamplingEngine = SamplingEnginePivotedPOD
			else:
			SamplingEngine = SamplingEnginePivoted
			self.samplingEngine = SamplingEngine(self.HFEngine,
			self.directionPivot,
			verbosity = self.verbosity,
			allowRepeatedSamples = True)

			@property
			def SMarginal(self):
			"""Value of SMarginal."""
			return self._SMarginal
			@SMarginal.setter
			def SMarginal(self, SMarginal):
			if not hasattr(SMarginal, "__len__"): SMarginal = [SMarginal]
			if any([s <= 0 for s in SMarginal]):
			raise RROMPyException("SMarginal must be positive.")
			if hasattr(self, "_SMarginal") and self._SMarginal is not None:
			Sold = tuple(self.SMarginal)
			else: Sold = -1
			self._SMarginal = SMarginal
			self._approxParameters["SMarginal"] = self.SMarginal
			if Sold != tuple(self.SMarginal):
			self.resetSamples()

			@property
			def polybasisMarginal(self):
			"""Value of polybasisMarginal."""
			return self._polybasisMarginal
			@polybasisMarginal.setter
			def polybasisMarginal(self, polybasisMarginal):
			try:
			polybasisMarginal = polybasisMarginal.upper().strip().replace(" ",
			"")
			if polybasisMarginal not in polybases:
			raise RROMPyException(
			"Prescribed marginal polybasis not recognized.")
			self._polybasisMarginal = polybasisMarginal
			except:
			RROMPyWarning(("Prescribed marginal polybasis not recognized. "
			"Overriding to 'MONOMIAL'."))
			self._polybasisMarginal = "MONOMIAL"
			self._approxParameters["polybasisMarginal"] = self.polybasisMarginal

			@property
			def MMarginal(self):
			"""Value of MMarginal."""
			return self._MMarginal
			@MMarginal.setter
			def MMarginal(self, MMarginal):
			if MMarginal < 0:
			raise RROMPyException("MMarginal must be non-negative.")
			self._MMarginal = MMarginal
			self._approxParameters["MMarginal"] = self.MMarginal

			@property
			def radialBasisMarginal(self):
			"""Value of radialBasisMarginal."""
			return self._radialBasisMarginal
			@radialBasisMarginal.setter
			def radialBasisMarginal(self, radialBasisMarginal):
			try:
			if radialBasisMarginal != 0:
			radialBasisMarginal = radialBasisMarginal.upper().strip()\
			.replace(" ","")
			if radialBasisMarginal not in rbbases:
			raise RROMPyException(("Prescribed marginal radialBasis "
			"not recognized."))
			self._radialBasisMarginal = radialBasisMarginal
			except:
			RROMPyWarning(("Prescribed marginal radialBasis not recognized. "
			"Overriding to 0."))
			self._radialBasisMarginal = 0
			self._approxParameters["radialBasisMarginal"] = (
			self.radialBasisMarginal)

			@property
			def radialBasisWeightsMarginal(self):
			"""Value of radialBasisWeightsMarginal."""
			return self._radialBasisWeightsMarginal
			@radialBasisWeightsMarginal.setter
			def radialBasisWeightsMarginal(self, radialBasisWeightsMarginal):
			self._radialBasisWeightsMarginal = radialBasisWeightsMarginal
			self._approxParameters["radialBasisWeightsMarginal"] = (
			self.radialBasisWeightsMarginal)

			@property
			def polybasisMarginalP(self):
			if self.radialBasisMarginal == 0:
			return self._polybasisMarginal
			return self._polybasisMarginal + "_" + self.radialBasisMarginal

			@property
			def interpRcondMarginal(self):
			"""Value of interpRcondMarginal."""
			return self._interpRcondMarginal
			@interpRcondMarginal.setter
			def interpRcondMarginal(self, interpRcondMarginal):
			self._interpRcondMarginal = interpRcondMarginal
			self._approxParameters["interpRcondMarginal"] = (
			self.interpRcondMarginal)

			@property
			def directionPivot(self):
			"""Value of directionPivot. Its assignment may reset snapshots."""
			return self._directionPivot
			@directionPivot.setter
			def directionPivot(self, directionPivot):
			if hasattr(self, '_directionPivot'):
			directionPivotOld = copy(self.directionPivot)
			else:
			directionPivotOld = None
			if (directionPivotOld is None
			or len(directionPivot) != len(directionPivotOld)
			or not directionPivot == directionPivotOld):
			self.resetSamples()
			self._directionPivot = directionPivot

			@property
			def directionMarginal(self):
			return [x for x in range(self.HFEngine.npar) \
			if x not in self.directionPivot]

			@property
			def rescalingExpPivot(self):
			return [self.HFEngine.rescalingExp[x] for x in self.directionPivot]

			@property
			def rescalingExpMarginal(self):
			return [self.HFEngine.rescalingExp[x] for x in self.directionMarginal]

			@property
			def muBoundsPivot(self):
			"""Value of muBoundsPivot."""
			return self.samplerPivot.lims

			@property
			def muBoundsMarginal(self):
			"""Value of muBoundsMarginal."""
			return self.samplerMarginal.lims

			@property
			def samplerPivot(self):
			"""Value of samplerPivot."""
			return self._samplerPivot
			@samplerPivot.setter
			def samplerPivot(self, samplerPivot):
			if 'generatePoints' not in dir(samplerPivot):
			raise RROMPyException("Pivot sampler type not recognized.")
			if hasattr(self, '_samplerPivot') and self._samplerPivot is not None:
			samplerOld = self.samplerPivot
			self._samplerPivot = samplerPivot
			self._approxParameters["samplerPivot"] = self.samplerPivot.__str__()
			if not 'samplerOld' in locals() or samplerOld != self.samplerPivot:
			self.resetSamples()

			@property
			def samplerMarginal(self):
			"""Value of samplerMarginal."""
			return self._samplerMarginal
			@samplerMarginal.setter
			def samplerMarginal(self, samplerMarginal):
			if 'generatePoints' not in dir(samplerMarginal):
			raise RROMPyException("Marginal sampler type not recognized.")
			if (hasattr(self, '_samplerMarginal')
			and self._samplerMarginal is not None):
			samplerOld = self.samplerMarginal
			self._samplerMarginal = samplerMarginal
			self._approxParameters["samplerMarginal"] = (
			self.samplerMarginal.__str__())
			if not 'samplerOld' in locals() or samplerOld != self.samplerMarginal:
			self.resetSamples()

			def resetSamples(self):
			"""Reset samples."""
			super().resetSamples()
			self._musMUniqueCN = None
			self._derMIdxs = None
			self._reorderM = None

			def setSamples(self, samplingEngine):
			"""Copy samplingEngine and samples."""
			super().setSamples(samplingEngine)
			self.mus = copy(self.samplingEngine[0].mus)
			for sEj in self.samplingEngine[1:]:
			self.mus.append(sEj.mus)

			def computeSnapshots(self):
			"""Compute snapshots of solution map."""
			RROMPyAssert(self._mode,
			message = "Cannot start snapshot computation.")
			Sp = np.prod(self.S)
			Seff = Sp * np.prod(self.SMarginal)
			if self.samplingEngine.nsamplesTot != Seff:
			self.resetSamples()
			vbMng(self, "INIT", "Starting computation of snapshots.", 5)
			self.musPivot = self.samplerPivot.generatePoints(self.S)
			self.musMarginal = self.samplerMarginal.generatePoints(
			self.SMarginal)
			self.mus = emptyParameterList()
			self.mus.reset((Seff, self.HFEngine.npar))
			self.samplingEngine.resetHistory(Seff // Sp)
			for j, muMarg in enumerate(self.musMarginal):
			for k in range(j * Sp, (j + 1) * Sp):
			self.mus.data[k, self.directionPivot] = (
			self.musPivot[k - j*Sp].data)
			self.mus.data[k, self.directionMarginal] = muMarg.data
			muMEff = [fp] * self.npar
			for k, x in enumerate(self.directionMarginal):
			muMEff[x] = muMarg(0, k)
			self.samplingEngine.HFEngineMarginalized = MarginalProxyEngine(
			self.HFEngine, list(muMEff))
			self.samplingEngine.iterSample(self.musPivot, j,
			homogeneized = self.homogeneized)
			if self.POD:
			self.samplingEngine.coalesceSamples(self.interpRcondMarginal)
			else:
			self.samplingEngine.coalesceSamples()
			vbMng(self, "DEL", "Done computing snapshots.", 5)

			def _setupMarginalInterpolationIndices(self):
			"""Setup parameters for polyvander."""
			RROMPyAssert(self._mode,
			message = "Cannot setup interpolation indices.")
			if (self._musMUniqueCN is None
			or len(self._reorderM) != len(self.musMarginal)):
			self._musMUniqueCN, musMIdxsTo, musMIdxs, musMCount = (
			self.centerNormalizeMarginal(self.musMarginal).unique(
			return_index = True,
			return_inverse = True,
			return_counts = True))
			self._musMUnique = self.musMarginal[musMIdxsTo]
			self._derMIdxs = [None] * len(self._musMUniqueCN)
			self._reorderM = np.empty(len(musMIdxs), dtype = int)
			filled = 0
			for j, cnt in enumerate(musMCount):
			self._derMIdxs[j] = nextDerivativeIndices([],
			self.musMarginal.shape[1], cnt)
			jIdx = np.nonzero(musMIdxs == j)[0]
			self._reorderM[jIdx] = np.arange(filled, filled + cnt)
			filled += cnt

			def _setupMarginalInterp(self):
			"""Compute marginal interpolator."""
			RROMPyAssert(self._mode, message = "Cannot setup numerator.")
			vbMng(self, "INIT", "Starting computation of marginal interpolator.",
			7)
			self._setupMarginalInterpolationIndices()
			MMarginal0 = copy(self.MMarginal)
			mI = []
			for j in range(len(self.musMarginal)):
			canonicalj = 1. * (np.arange(len(self.musMarginal)) == j)
			self._MMarginal = MMarginal0
			while self.MMarginal >= 0:
			if self.radialBasisMarginal == 0:
			p = PI()
			wellCond, msg = p.setupByInterpolation(
			self._musMUniqueCN, canonicalj, self.MMarginal,
			self.polybasisMarginal, self.verbosity >= 5, True,
			{"derIdxs": self._derMIdxs,
			"reorder": self._reorderM,
			"scl": np.power(self.scaleFactorMarginal, -1.)},
			{"rcond": self.interpRcondMarginal})
			else:
			p = RBI()
			wellCond, msg = p.setupByInterpolation(
			self._musMUniqueCN, canonicalj, self.MMarginal,
			self.polybasisMarginalP,
			self.radialBasisWeightsMarginal,
			self.verbosity >= 5, True,
			{"derIdxs": self._derMIdxs,
			"reorder": self._reorderM,
			"scl": np.power(self.scaleFactorMarginal, -1.)},
			{"rcond": self.interpRcondMarginal})
			vbMng(self, "MAIN", msg, 5)
			if wellCond: break
			RROMPyWarning(("Polyfit is poorly conditioned. Reducing "
			"MMarginal by 1."))
			self.MMarginal -= 1
			if self.MMarginal < 0:
			raise RROMPyException(("Instability in computation of "
			"marginal interpolator. Aborting."))
			mI = mI + [copy(p)]
			vbMng(self, "DEL", "Done computing marginal interpolator.", 7)
			return mI

			def normApprox(self, mu:paramList, homogeneized : bool = False) -> float:
			"""
			Compute norm of approximant at arbitrary parameter.

			Args:
			mu: Target parameter.
			homogeneized(optional): Whether to remove Dirichlet BC. Defaults to
			False.

			Returns:
			Target norm of approximant.
			"""
			if not self.POD or self.homogeneized != homogeneized:
			return super().normApprox(mu, homogeneized)
			return np.linalg.norm(self.getApproxReduced(mu).data, axis = 0)

			def computeScaleFactor(self):
			"""Compute parameter rescaling factor."""
			RROMPyAssert(self._mode, message = "Cannot compute rescaling factor.")
			self.scaleFactorPivot = .5 * np.abs(
			self.muBoundsPivot[0] ** self.rescalingExpPivot
			- self.muBoundsPivot[1] ** self.rescalingExpPivot)
			self.scaleFactorMarginal = .5 * np.abs(
			self.muBoundsMarginal[0] ** self.rescalingExpMarginal
			- self.muBoundsMarginal[1] ** self.rescalingExpMarginal)
			self.scaleFactor = np.empty(self.npar)
			self.scaleFactor[self.directionPivot] = self.scaleFactorPivot
			self.scaleFactor[self.directionMarginal] = self.scaleFactorMarginal

			def centerNormalizeMarginal(self, mu : paramList = [],
			mu0 : paramVal = None) -> paramList:
			"""
			Compute normalized parameter to be plugged into approximant.

			Args:
			mu: Parameter(s) 1.
			mu0: Parameter(s) 2. If None, set to self.mu0.

			Returns:
			Normalized parameter.
			"""
			return self.trainedModel.centerNormalizeMarginal(mu, mu0)

RROMPy/rrompy/reduction_methods/pivoting/generic_pivoted_approximant.py16a7048180acmaster

generic_pivoted_approximant.pyView Options

RROMPy/rrompy/reduction_methods/pivoting/generic_pivoted_approximant.py
16a7048180acmaster

generic_pivoted_approximant.py
View Options