diff --git a/rrompy/parameter/parameter_sampling/generic_quadrature_sampler.py b/rrompy/parameter/parameter_sampling/generic_quadrature_sampler.py index 60a77c5..e9f0f96 100644 --- a/rrompy/parameter/parameter_sampling/generic_quadrature_sampler.py +++ b/rrompy/parameter/parameter_sampling/generic_quadrature_sampler.py @@ -1,52 +1,49 @@ # 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 . # from .generic_sampler import GenericSampler from rrompy.utilities.base.types import DictAny, paramList from rrompy.utilities.exception_manager import RROMPyException _allowedQuadratureKinds = ["UNIFORM", "CHEBYSHEV", "EXTENDEDCHEBYSHEV", "GAUSSLEGENDRE", "EXTENDEDGAUSSLEGENDRE"] __all__ = ['GenericQuadratureSampler'] class GenericQuadratureSampler(GenericSampler): """Generator of quadrature sample points.""" def __init__(self, lims:paramList, kind : str = "UNIFORM", parameterMap : DictAny = 1.): super().__init__(lims = lims, parameterMap = parameterMap) self._allowedKinds = _allowedQuadratureKinds self.kind = kind def __str__(self) -> str: return "{}_{}".format(super().__str__(), self.kind) - def __repr__(self) -> str: - return self.__str__() + " at " + hex(id(self)) - @property def kind(self): """Value of kind.""" return self._kind @kind.setter def kind(self, kind): kind = kind.upper().strip().replace(" ","") if kind not in self._allowedKinds: raise RROMPyException("Generator kind not recognized.") self._kind = kind diff --git a/rrompy/parameter/parameter_sampling/generic_random_sampler.py b/rrompy/parameter/parameter_sampling/generic_random_sampler.py index e10b832..8341c46 100644 --- a/rrompy/parameter/parameter_sampling/generic_random_sampler.py +++ b/rrompy/parameter/parameter_sampling/generic_random_sampler.py @@ -1,51 +1,78 @@ # 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 . # +from abc import abstractmethod +import numpy as np from .generic_sampler import GenericSampler -from rrompy.utilities.base.types import DictAny, paramList +from rrompy.utilities.base.types import List, DictAny, paramList from rrompy.utilities.exception_manager import RROMPyException +from rrompy.parameter.parameter_list import emptyParameterList _allowedRandomKinds = ["UNIFORM", "HALTON", "SOBOL"] __all__ = ['GenericRandomSampler'] class GenericRandomSampler(GenericSampler): """Generator of random sample points.""" def __init__(self, lims:paramList, kind : str = "UNIFORM", - parameterMap : DictAny = 1., seed : int = 42): + parameterMap : DictAny = 1., refinementFactor : float = 1., + seed : int = 42): super().__init__(lims = lims, parameterMap = parameterMap) self._allowedKinds = _allowedRandomKinds self.kind = kind + self.refinementFactor = refinementFactor self.seed = seed - + self.reset() + def __str__(self) -> str: return "{}_{}".format(super().__str__(), self.kind) - def __repr__(self) -> str: - return self.__str__() + " at " + hex(id(self)) + @property + def npoints(self): + """Number of points.""" + return len(self.points) @property def kind(self): """Value of kind.""" return self._kind @kind.setter def kind(self, kind): if kind.upper() not in self._allowedKinds: raise RROMPyException("Generator kind not recognized.") self._kind = kind.upper() + + def reset(self): + self.points = emptyParameterList() + + def generatePoints(self, n:int, reorder : bool = True) -> paramList: + """Array of quadrature points.""" + if self.kind == "UNIFORM": + np.random.seed(self.seed) + else: + self.seedLoc = self.seed + self.reset() + rF, self.refinementFactor = self.refinementFactor, 1. + _ = self.refine([None] * n) + self.refinementFactor = rF + return self.points + + @abstractmethod + def refine(self, active : List[int] = None) -> List[int]: + pass diff --git a/rrompy/parameter/parameter_sampling/manual_sampler.py b/rrompy/parameter/parameter_sampling/manual_sampler.py index 8a1afdd..a00c0f5 100644 --- a/rrompy/parameter/parameter_sampling/manual_sampler.py +++ b/rrompy/parameter/parameter_sampling/manual_sampler.py @@ -1,64 +1,61 @@ # 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 . # import numpy as np from copy import deepcopy as copy from .generic_sampler import GenericSampler from rrompy.utilities.base.types import List, DictAny, paramList from rrompy.parameter import checkParameterList __all__ = ['ManualSampler'] class ManualSampler(GenericSampler): """Manual generator of sample points.""" def __init__(self, lims:paramList, points:paramList, parameterMap : DictAny = 1., normalFoci : List[np.complex] = [-1., 1.]): super().__init__(lims = lims, parameterMap = parameterMap) self.points = points self._normalFoci = normalFoci def normalFoci(self, d : int = 0): return self._normalFoci @property def points(self): """Value of points.""" return self._points @points.setter def points(self, points): points = checkParameterList(points, self.npar) self._points = points def __str__(self) -> str: return "{}[{}]".format(self.name(), "_".join(map(str, self.points))) - def __repr__(self) -> str: - return self.__str__() + " at " + hex(id(self)) - def generatePoints(self, n:int, reorder : bool = True) -> paramList: """Array of sample points.""" if n > len(self.points): pts = copy(self.points) for j in range(int(np.ceil(n / len(self.points)))): pts.append(self.points) else: pts = self.points x = checkParameterList(pts[list(range(n))], self.npar) return x diff --git a/rrompy/parameter/parameter_sampling/segment/random_sampler.py b/rrompy/parameter/parameter_sampling/segment/random_sampler.py index 3e0c8fb..c1515f5 100644 --- a/rrompy/parameter/parameter_sampling/segment/random_sampler.py +++ b/rrompy/parameter/parameter_sampling/segment/random_sampler.py @@ -1,43 +1,55 @@ # 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 . # +from numbers import Number import numpy as np from rrompy.parameter.parameter_sampling.generic_random_sampler import ( GenericRandomSampler) -from rrompy.utilities.numerical import haltonGenerate, sobolGenerate -from rrompy.utilities.base.types import paramList +from rrompy.utilities.numerical.halton import haltonGenerate +from rrompy.utilities.numerical.sobol import sobolGenerate +from rrompy.utilities.base.types import List __all__ = ['RandomSampler'] class RandomSampler(GenericRandomSampler): """Generator of (quasi-)random sample points.""" - def generatePoints(self, n:int, reorder : bool = True) -> paramList: - """Array of quadrature points.""" + def refine(self, active : List[int] = None) -> List[int]: + if active is None: + n = self.npoints + elif isinstance(active, (Number,)): + n = active + else: + n = len(active) + n = int(n * self.refinementFactor) if self.kind == "UNIFORM": - np.random.seed(self.seed) xmat = np.random.uniform(size = (n, self.npar)) elif self.kind == "HALTON": - xmat = haltonGenerate(self.npar, n, self.seed) + xmat, self.seedLoc = haltonGenerate(self.npar, n, self.seedLoc, + return_seed = True) else: - xmat = sobolGenerate(self.npar, n, self.seed) + xmat, self.seedLoc = sobolGenerate(self.npar, n, self.seedLoc, + return_seed = True) limsE = self.mapParameterList(self.lims) for d in range(self.npar): a, b = limsE(d) xmat[:, d] = a + (b - a) * xmat[:, d] - return self.mapParameterList(xmat, "B") + pts = self.mapParameterList(xmat, "B") + idx = np.arange(n, dtype = int) + len(self.points) + for pj in pts: self.points.append(pj) + return list(idx) diff --git a/rrompy/parameter/parameter_sampling/shape/generic_shape_quadrature_sampler.py b/rrompy/parameter/parameter_sampling/shape/generic_shape_quadrature_sampler.py index 6e9ab16..65a4701 100644 --- a/rrompy/parameter/parameter_sampling/shape/generic_shape_quadrature_sampler.py +++ b/rrompy/parameter/parameter_sampling/shape/generic_shape_quadrature_sampler.py @@ -1,52 +1,33 @@ # 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 . # from .generic_shape_sampler import GenericShapeSampler from rrompy.parameter.parameter_sampling.generic_quadrature_sampler import ( - _allowedQuadratureKinds) + GenericQuadratureSampler) from rrompy.utilities.base.types import List, DictAny, paramList -from rrompy.utilities.exception_manager import RROMPyException __all__ = ['GenericShapeQuadratureSampler'] -class GenericShapeQuadratureSampler(GenericShapeSampler): +class GenericShapeQuadratureSampler(GenericQuadratureSampler, + GenericShapeSampler): """Generator of quadrature sample points on shapes.""" def __init__(self, lims:paramList, kind : str = "UNIFORM", axisRatios : List[float] = None, parameterMap : DictAny = 1.): - super().__init__(lims = lims, axisRatios = axisRatios, - parameterMap = parameterMap) - self._allowedKinds = _allowedQuadratureKinds - self.kind = kind - - def __str__(self) -> str: - return "{}_{}".format(super().__str__(), self.kind) - - def __repr__(self) -> str: - return self.__str__() + " at " + hex(id(self)) - - @property - def kind(self): - """Value of kind.""" - return self._kind - @kind.setter - def kind(self, kind): - kind = kind.upper().strip().replace(" ","") - if kind not in self._allowedKinds: - raise RROMPyException("Generator kind not recognized.") - self._kind = kind + super().__init__(lims = lims, kind = kind, parameterMap = parameterMap) + self.axisRatios = axisRatios diff --git a/rrompy/parameter/parameter_sampling/shape/generic_shape_random_sampler.py b/rrompy/parameter/parameter_sampling/shape/generic_shape_random_sampler.py index 536f01b..75c4fa3 100644 --- a/rrompy/parameter/parameter_sampling/shape/generic_shape_random_sampler.py +++ b/rrompy/parameter/parameter_sampling/shape/generic_shape_random_sampler.py @@ -1,53 +1,35 @@ # 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 . # from .generic_shape_sampler import GenericShapeSampler from rrompy.parameter.parameter_sampling.generic_random_sampler import ( - _allowedRandomKinds) + GenericRandomSampler) from rrompy.utilities.base.types import List, DictAny, paramList -from rrompy.utilities.exception_manager import RROMPyException __all__ = ['GenericShapeRandomSampler'] -class GenericShapeRandomSampler(GenericShapeSampler): +class GenericShapeRandomSampler(GenericRandomSampler, GenericShapeSampler): """Generator of random sample points on shapes.""" def __init__(self, lims:paramList, kind : str = "UNIFORM", axisRatios : List[float] = None, parameterMap : DictAny = 1., - seed : int = 42): - super().__init__(lims = lims, axisRatios = axisRatios, - parameterMap = parameterMap) - self._allowedKinds = _allowedRandomKinds - self.kind = kind - self.seed = seed - - def __str__(self) -> str: - return "{}_{}".format(super().__str__(), self.kind) - - def __repr__(self) -> str: - return self.__str__() + " at " + hex(id(self)) - - @property - def kind(self): - """Value of kind.""" - return self._kind - @kind.setter - def kind(self, kind): - if kind.upper() not in self._allowedKinds: - raise RROMPyException("Generator kind not recognized.") - self._kind = kind.upper() + refinementFactor : float = 1., seed : int = 42): + super().__init__(lims = lims, kind = kind, parameterMap = parameterMap, + refinementFactor = refinementFactor) + self.axisRatios = axisRatios + \ No newline at end of file diff --git a/rrompy/parameter/parameter_sampling/shape/random_box_sampler.py b/rrompy/parameter/parameter_sampling/shape/random_box_sampler.py index cc48aea..f1e872f 100644 --- a/rrompy/parameter/parameter_sampling/shape/random_box_sampler.py +++ b/rrompy/parameter/parameter_sampling/shape/random_box_sampler.py @@ -1,56 +1,69 @@ # 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 . # +from numbers import Number import numpy as np from .generic_shape_random_sampler import GenericShapeRandomSampler from rrompy.utilities.numerical import haltonGenerate, sobolGenerate -from rrompy.utilities.base.types import paramList +from rrompy.utilities.base.types import List __all__ = ['RandomBoxSampler'] class RandomBoxSampler(GenericShapeRandomSampler): """Generator of (quasi-)random sample points on boxes.""" - def generatePoints(self, n:int, reorder : bool = True) -> paramList: - """Array of quadrature points.""" + def refine(self, active : List[int] = None) -> List[int]: + if active is None: + n = self.npoints + elif isinstance(active, (Number,)): + n = int(active) + else: + n = len(active) + n = int(n * self.refinementFactor) nEff = int(np.ceil(n * np.prod( [max(x, 1. / x) for x in self.axisRatios]))) xmat2 = [] while len(xmat2) < n: if self.kind == "UNIFORM": - np.random.seed(self.seed) xmat2 = np.random.uniform(size = (nEff, 2 * self.npar)) elif self.kind == "HALTON": - xmat2 = haltonGenerate(2 * self.npar, nEff, self.seed) + xmat2, self.seedLoc = haltonGenerate(2 * self.npar, nEff, + self.seedLoc, + return_seed = True) else: - xmat2 = sobolGenerate(2 * self.npar, nEff, self.seed) + xmat2, self.seedLoc = sobolGenerate(2 * self.npar, nEff, + self.seed, + return_seed = True) for d in range(self.npar): ax = self.axisRatios[d] if ax <= 1.: xmat2 = xmat2[xmat2[:, 2 * d + 1] <= ax] else: xmat2 = xmat2[xmat2[:, 2 * d] <= 1. / ax] xmat2[:, 2 * d : 2 * d + 2] *= ax nEff += 1 xmat = np.empty((n, self.npar), dtype = np.complex) limsE = self.mapParameterList(self.lims) for d in range(self.npar): a, b = limsE(d) xmat[:, d] = a + (b - a) * (xmat2[: n, 2 * d] + 1.j * self.axisRatios[d] * (xmat2[: n, 2 * d + 1] - .5)) - return self.mapParameterList(xmat, "B") + pts = self.mapParameterList(xmat, "B") + idx = np.arange(n, dtype = int) + len(self.points) + for pj in pts: self.points.append(pj) + return list(idx) diff --git a/rrompy/parameter/parameter_sampling/shape/random_circle_sampler.py b/rrompy/parameter/parameter_sampling/shape/random_circle_sampler.py index 1c237b2..5a1de1f 100644 --- a/rrompy/parameter/parameter_sampling/shape/random_circle_sampler.py +++ b/rrompy/parameter/parameter_sampling/shape/random_circle_sampler.py @@ -1,59 +1,72 @@ # 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 . # +from numbers import Number import numpy as np from .generic_shape_random_sampler import GenericShapeRandomSampler from rrompy.utilities.numerical import (haltonGenerate, sobolGenerate, potential) -from rrompy.utilities.base.types import paramList +from rrompy.utilities.base.types import List __all__ = ['RandomCircleSampler'] class RandomCircleSampler(GenericShapeRandomSampler): """Generator of (quasi-)random sample points on ellipses.""" - def generatePoints(self, n:int, reorder : bool = True) -> paramList: - """Array of quadrature points.""" + def refine(self, active : List[int] = None) -> List[int]: + if active is None: + n = self.npoints + elif isinstance(active, (Number,)): + n = active + else: + n = len(active) + n = int(n * self.refinementFactor) nEff = int(np.ceil(n * (4. / np.pi) ** self.npar * np.prod( [max(x, 1. / x) for x in self.axisRatios]))) xmat2 = [] while len(xmat2) < n: if self.kind == "UNIFORM": - np.random.seed(self.seed) xmat2 = np.random.uniform(size = (nEff, 2 * self.npar)) elif self.kind == "HALTON": - xmat2 = haltonGenerate(2 * self.npar, nEff, self.seed) + xmat2, self.seedLoc = haltonGenerate(2 * self.npar, nEff, + self.seedLoc, + return_seed = True) else: - xmat2 = sobolGenerate(2 * self.npar, nEff, self.seed) + xmat2, self.seedLoc = sobolGenerate(2 * self.npar, nEff, + self.seed, + return_seed = True) xmat2 = xmat2 * 2. - 1. for d in range(self.npar): ax = self.axisRatios[d] if ax > 1.: xmat2[:, 2 * d : 2 * d + 2] *= ax Z = xmat2[:, 2 * d] + 1.j * ax * xmat2[:, 2 * d + 1] ptscore = potential(Z, self.normalFoci(d)) xmat2 = xmat2[ptscore <= self.groundPotential(d)] nEff += 1 xmat = np.empty((n, self.npar), dtype = np.complex) limsE = self.mapParameterList(self.lims) for d in range(self.npar): ax = self.axisRatios[d] a, b = limsE(d) c, r = (a + b) / 2., (a - b) / 2. xmat[:, d] = c + r * (xmat2[: n, 2 * d] + 1.j * ax * xmat2[: n, 2 * d + 1]) - return self.mapParameterList(xmat, "B") + pts = self.mapParameterList(xmat, "B") + idx = np.arange(n, dtype = int) + len(self.points) + for pj in pts: self.points.append(pj) + return list(idx) diff --git a/rrompy/reduction_methods/pivoted/greedy/generic_pivoted_greedy_approximant.py b/rrompy/reduction_methods/pivoted/greedy/generic_pivoted_greedy_approximant.py index 42b1ed0..9b571f7 100644 --- a/rrompy/reduction_methods/pivoted/greedy/generic_pivoted_greedy_approximant.py +++ b/rrompy/reduction_methods/pivoted/greedy/generic_pivoted_greedy_approximant.py @@ -1,836 +1,837 @@ # 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 . # from abc import abstractmethod from copy import deepcopy as copy import numpy as np from matplotlib import pyplot as plt from rrompy.reduction_methods.pivoted.generic_pivoted_approximant import ( GenericPivotedApproximantBase, GenericPivotedApproximantNoMatch, GenericPivotedApproximant) from rrompy.reduction_methods.pivoted.gather_pivoted_approximant import ( gatherPivotedApproximant) from rrompy.utilities.base.types import (Np1D, Np2D, Tuple, List, paramVal, paramList, ListAny) from rrompy.utilities.base import verbosityManager as vbMng from rrompy.utilities.numerical.point_matching import (pointMatching, chordalMetricAdjusted, potential) from rrompy.utilities.exception_manager import (RROMPyException, RROMPyAssert, RROMPyWarning) from rrompy.parameter import emptyParameterList from rrompy.utilities.parallel import (masterCore, indicesScatter, arrayGatherv, isend) __all__ = ['GenericPivotedGreedyApproximantNoMatch', 'GenericPivotedGreedyApproximant'] class GenericPivotedGreedyApproximantBase(GenericPivotedApproximantBase): _allowedEstimatorKindsMarginal = ["LEAVE_ONE_OUT", "LOOK_AHEAD", "LOOK_AHEAD_RECOVER", "NONE"] def __init__(self, *args, **kwargs): self._preInit() self._addParametersToList(["matchingWeightError", "cutOffToleranceError", "errorEstimatorKindMarginal", "greedyTolMarginal", "maxIterMarginal"], [0., "AUTO", "NONE", 1e-1, 1e2]) super().__init__(*args, **kwargs) self._postInit() @property def scaleFactorDer(self): """Value of scaleFactorDer.""" if self._scaleFactorDer == "NONE": return 1. if self._scaleFactorDer == "AUTO": return self._scaleFactorOldPivot return self._scaleFactorDer @scaleFactorDer.setter def scaleFactorDer(self, scaleFactorDer): if isinstance(scaleFactorDer, (str,)): scaleFactorDer = scaleFactorDer.upper() elif hasattr(scaleFactorDer, "__len__"): scaleFactorDer = list(scaleFactorDer) self._scaleFactorDer = scaleFactorDer self._approxParameters["scaleFactorDer"] = self._scaleFactorDer @property def samplerMarginal(self): """Value of samplerMarginal.""" return self._samplerMarginal @samplerMarginal.setter def samplerMarginal(self, samplerMarginal): if 'refine' not in dir(samplerMarginal): raise RROMPyException("Marginal sampler type not recognized.") GenericPivotedApproximantBase.samplerMarginal.fset(self, samplerMarginal) @property def errorEstimatorKindMarginal(self): """Value of errorEstimatorKindMarginal.""" return self._errorEstimatorKindMarginal @errorEstimatorKindMarginal.setter def errorEstimatorKindMarginal(self, errorEstimatorKindMarginal): errorEstimatorKindMarginal = errorEstimatorKindMarginal.upper() if errorEstimatorKindMarginal not in ( self._allowedEstimatorKindsMarginal): RROMPyWarning(("Marginal error estimator kind not recognized. " "Overriding to 'NONE'.")) errorEstimatorKindMarginal = "NONE" self._errorEstimatorKindMarginal = errorEstimatorKindMarginal self._approxParameters["errorEstimatorKindMarginal"] = ( self.errorEstimatorKindMarginal) @property def matchingWeightError(self): """Value of matchingWeightError.""" return self._matchingWeightError @matchingWeightError.setter def matchingWeightError(self, matchingWeightError): self._matchingWeightError = matchingWeightError self._approxParameters["matchingWeightError"] = ( self.matchingWeightError) @property def cutOffToleranceError(self): """Value of cutOffToleranceError.""" return self._cutOffToleranceError @cutOffToleranceError.setter def cutOffToleranceError(self, cutOffToleranceError): if isinstance(cutOffToleranceError, (str,)): cutOffToleranceError = cutOffToleranceError.upper()\ .strip().replace(" ","") if cutOffToleranceError != "AUTO": RROMPyWarning(("String value of cutOffToleranceError not " "recognized. Overriding to 'AUTO'.")) cutOffToleranceError == "AUTO" self._cutOffToleranceError = cutOffToleranceError self._approxParameters["cutOffToleranceError"] = ( self.cutOffToleranceError) @property def greedyTolMarginal(self): """Value of greedyTolMarginal.""" return self._greedyTolMarginal @greedyTolMarginal.setter def greedyTolMarginal(self, greedyTolMarginal): if greedyTolMarginal < 0: raise RROMPyException("greedyTolMarginal must be non-negative.") if (hasattr(self, "_greedyTolMarginal") and self.greedyTolMarginal is not None): greedyTolMarginalold = self.greedyTolMarginal else: greedyTolMarginalold = -1 self._greedyTolMarginal = greedyTolMarginal self._approxParameters["greedyTolMarginal"] = self.greedyTolMarginal if greedyTolMarginalold != self.greedyTolMarginal: self.resetSamples() @property def maxIterMarginal(self): """Value of maxIterMarginal.""" return self._maxIterMarginal @maxIterMarginal.setter def maxIterMarginal(self, maxIterMarginal): if maxIterMarginal <= 0: raise RROMPyException("maxIterMarginal must be positive.") if (hasattr(self, "_maxIterMarginal") and self.maxIterMarginal is not None): maxIterMarginalold = self.maxIterMarginal else: maxIterMarginalold = -1 self._maxIterMarginal = maxIterMarginal self._approxParameters["maxIterMarginal"] = self.maxIterMarginal if maxIterMarginalold != self.maxIterMarginal: self.resetSamples() def resetSamples(self): """Reset samples.""" super().resetSamples() if not hasattr(self, "_temporaryPivot"): self._mus = emptyParameterList() self._musMarginal = emptyParameterList() if hasattr(self, "samplerMarginal"): self.samplerMarginal.reset() if hasattr(self, "samplingEngine") and self.samplingEngine is not None: self.samplingEngine.resetHistory() def _getPolesResExact(self, HITest, foci:Tuple[float, float], ground:float) -> Tuple[Np1D, Np2D]: if self.cutOffToleranceError == "AUTO": cutOffTolErr = self.cutOffTolerance else: cutOffTolErr = self.cutOffToleranceError polesEx = copy(HITest.poles) idxExEff = np.where(potential(polesEx, foci) - ground <= cutOffTolErr * ground)[0] if self.matchingWeightError != 0: resEx = HITest.coeffs[idxExEff] else: resEx = None return polesEx[idxExEff], resEx def _getDistanceApp(self, polesEx:Np1D, resEx:Np2D, muTest:paramVal, foci:Tuple[float, float], ground:float) -> float: if self.cutOffToleranceError == "AUTO": cutOffTolErr = self.cutOffTolerance else: cutOffTolErr = self.cutOffToleranceError polesAp = self.trainedModel.interpolateMarginalPoles(muTest)[0] idxApEff = np.where(potential(polesAp, foci) - ground <= cutOffTolErr * ground)[0] polesAp = polesAp[idxApEff] if self.matchingWeightError != 0: resAp = self.trainedModel.interpolateMarginalCoeffs(muTest)[0][ idxApEff, :] resEx = self.trainedModel.data.projMat[:, : resEx.shape[1]].dot(resEx.T) resAp = self.trainedModel.data.projMat[:, : resAp.shape[1]].dot(resAp.T) else: resAp = None dist = chordalMetricAdjusted(polesEx, polesAp, self.matchingWeightError, resEx, resAp, self.HFEngine, False) pmR, pmC = pointMatching(dist) return np.mean(dist[pmR, pmC]) def getErrorEstimatorMarginalLeaveOneOut(self) -> Np1D: nTest = len(self.trainedModel.data.musMarginal) self._musMarginalTestIdxs = np.arange(nTest) if nTest <= 1: err = np.empty(nTest) err[:] = np.inf return err idx, sizes = indicesScatter(nTest, return_sizes = True) err = [] if len(idx) > 0: _tMdataFull = copy(self.trainedModel.data) _musMExcl = None self.verbosity -= 35 self.trainedModel.verbosity -= 35 foci = self.samplerPivot.normalFoci() ground = self.samplerPivot.groundPotential() for i, j in enumerate(idx): jEff = j - (i > 0) muTest = self.trainedModel.data.musMarginal[jEff] polesEx, resEx = self._getPolesResExact( self.trainedModel.data.HIs[jEff], foci, ground) if i > 0: self.musMarginal.insert(_musMExcl, j - 1) _musMExcl = self.musMarginal[j] self.musMarginal.pop(j) if len(polesEx) == 0: err += [0.] continue self._updateTrainedModelMarginalSamples([j]) self._finalizeMarginalization() err += [self._getDistanceApp(polesEx, resEx, muTest, foci, ground)] self._updateTrainedModelMarginalSamples() self.musMarginal.insert(_musMExcl, idx[-1]) self.verbosity += 35 self.trainedModel.verbosity += 35 self.trainedModel.data = _tMdataFull return arrayGatherv(np.array(err), sizes) def getErrorEstimatorMarginalLookAhead(self) -> Np1D: if not hasattr(self.trainedModel, "_musMExcl"): err = np.zeros(0) err[:] = np.inf self._musMarginalTestIdxs = np.zeros(0, dtype = int) return err self._musMarginalTestIdxs = np.array(self.trainedModel._idxExcl, dtype = int) idx, sizes = indicesScatter(len(self.trainedModel._musMExcl), return_sizes = True) err = [] if len(idx) > 0: self.verbosity -= 35 self.trainedModel.verbosity -= 35 foci = self.samplerPivot.normalFoci() ground = self.samplerPivot.groundPotential() for j in idx: muTest = self.trainedModel._musMExcl[j] HITest = self.trainedModel._HIsExcl[j] polesEx, resEx = self._getPolesResExact(HITest, foci, ground) if len(polesEx) == 0: err += [0.] continue err += [self._getDistanceApp(polesEx, resEx, muTest, foci, ground)] self.verbosity += 35 self.trainedModel.verbosity += 35 return arrayGatherv(np.array(err), sizes) def getErrorEstimatorMarginalNone(self) -> Np1D: nErr = len(self.trainedModel.data.musMarginal) self._musMarginalTestIdxs = np.arange(nErr) return (1. + self.greedyTolMarginal) * np.ones(nErr) def errorEstimatorMarginal(self, return_max : bool = False) -> Np1D: vbMng(self.trainedModel, "INIT", "Evaluating error estimator at mu = {}.".format( self.trainedModel.data.musMarginal), 10) if self.errorEstimatorKindMarginal == "LEAVE_ONE_OUT": err = self.getErrorEstimatorMarginalLeaveOneOut() elif self.errorEstimatorKindMarginal[: 10] == "LOOK_AHEAD": err = self.getErrorEstimatorMarginalLookAhead() else:#if self.errorEstimatorKindMarginal == "NONE": err = self.getErrorEstimatorMarginalNone() vbMng(self.trainedModel, "DEL", "Done evaluating error estimator", 10) if not return_max: return err idxMaxEst = np.where(err > self.greedyTolMarginal)[0] maxErr = err[idxMaxEst] if self.errorEstimatorKindMarginal == "NONE": maxErr = None return err, idxMaxEst, maxErr def plotEstimatorMarginal(self, est:Np1D, idxMax:List[int], estMax:List[float]): if self.errorEstimatorKindMarginal == "NONE": return if (not (np.any(np.isnan(est)) or np.any(np.isinf(est))) and masterCore()): fig = plt.figure(figsize = plt.figaspect(1. / self.nparMarginal)) for jpar in range(self.nparMarginal): ax = fig.add_subplot(1, self.nparMarginal, 1 + jpar) if self.errorEstimatorKindMarginal == "LEAVE_ONE_OUT": musre = copy(self.trainedModel.data.musMarginal.re.data) else:#if self.errorEstimatorKindMarginal[: 10] == "LOOK_AHEAD": if not hasattr(self.trainedModel, "_musMExcl"): return musre = np.real(self.trainedModel._musMExcl) if len(idxMax) > 0 and estMax is not None: maxrej = musre[idxMax, jpar] errCP = copy(est) idx = np.delete(np.arange(self.nparMarginal), jpar) while len(musre) > 0: if self.nparMarginal == 1: currIdx = np.arange(len(musre)) else: currIdx = np.where(np.isclose(np.sum( np.abs(musre[:, idx] - musre[0, idx]), 1), 0.))[0] currIdxSorted = currIdx[np.argsort(musre[currIdx, jpar])] ax.semilogy(musre[currIdxSorted, jpar], errCP[currIdxSorted], 'k.-', linewidth = 1) musre = np.delete(musre, currIdx, 0) errCP = np.delete(errCP, currIdx) ax.semilogy(self.musMarginal.re(jpar), (self.greedyTolMarginal,) * len(self.musMarginal), '*m') if len(idxMax) > 0 and estMax is not None: ax.semilogy(maxrej, estMax, 'xr') + ax.set_xlim(*list(self.samplerMarginal.lims.re(jpar))) ax.grid() plt.tight_layout() plt.show() def _addMarginalSample(self, mus:paramList): mus = self.checkParameterListMarginal(mus) if len(mus) == 0: return self._nmusOld, nmus = len(self.musMarginal), len(mus) if (hasattr(self, "trainedModel") and self.trainedModel is not None and hasattr(self.trainedModel, "_musMExcl")): self._nmusOld += len(self.trainedModel._musMExcl) vbMng(self, "MAIN", ("Adding marginal sample point{} no. {}{} at {} to training " "set.").format("s" * (nmus > 1), self._nmusOld + 1, "--{}".format(self._nmusOld + nmus) * (nmus > 1), mus), 3) self.musMarginal.append(mus) self.setupApproxPivoted(mus) self._poleMatching() del self._nmusOld if (self.errorEstimatorKindMarginal[: 10] == "LOOK_AHEAD" and not self.firstGreedyIterM): ubRange = len(self.trainedModel.data.musMarginal) if hasattr(self.trainedModel, "_idxExcl"): shRange = len(self.trainedModel._musMExcl) else: shRange = 0 testIdxs = list(range(ubRange + shRange - len(mus), ubRange + shRange)) for j in testIdxs[::-1]: self.musMarginal.pop(j - shRange) if hasattr(self.trainedModel, "_idxExcl"): testIdxs = self.trainedModel._idxExcl + testIdxs self._updateTrainedModelMarginalSamples(testIdxs) self._finalizeMarginalization() self._SMarginal = len(self.musMarginal) self._approxParameters["SMarginal"] = self.SMarginal self.trainedModel.data.approxParameters["SMarginal"] = self.SMarginal def greedyNextSampleMarginal(self, muidx:List[int], plotEst : str = "NONE") \ -> Tuple[Np1D, List[int], float, paramVal]: RROMPyAssert(self._mode, message = "Cannot add greedy sample.") muidx = self._musMarginalTestIdxs[muidx] if (self.errorEstimatorKindMarginal[: 10] == "LOOK_AHEAD" and not self.firstGreedyIterM): if not hasattr(self.trainedModel, "_idxExcl"): raise RROMPyException(("Sample index to be added not present " "in trained model.")) testIdxs = copy(self.trainedModel._idxExcl) skippedIdx = 0 for cj, j in enumerate(self.trainedModel._idxExcl): if j in muidx: testIdxs.pop(skippedIdx) self.musMarginal.insert(self.trainedModel._musMExcl[cj], j - skippedIdx) else: skippedIdx += 1 if len(self.trainedModel._idxExcl) < (len(muidx) + len(testIdxs)): raise RROMPyException(("Sample index to be added not present " "in trained model.")) self._updateTrainedModelMarginalSamples(testIdxs) self._SMarginal = len(self.musMarginal) self._approxParameters["SMarginal"] = self.SMarginal self.trainedModel.data.approxParameters["SMarginal"] = ( self.SMarginal) self.firstGreedyIterM = False idxAdded = self.samplerMarginal.refine(muidx) self._addMarginalSample(self.samplerMarginal.points[idxAdded]) errorEstTest, muidx, maxErrorEst = self.errorEstimatorMarginal(True) if plotEst == "ALL": self.plotEstimatorMarginal(errorEstTest, muidx, maxErrorEst) return (errorEstTest, muidx, maxErrorEst, self.samplerMarginal.points[muidx]) def _preliminaryTrainingMarginal(self): """Initialize starting snapshots of solution map.""" RROMPyAssert(self._mode, message = "Cannot start greedy algorithm.") if np.sum(self.samplingEngine.nsamples) > 0: return self.resetSamples() self._addMarginalSample(self.samplerMarginal.generatePoints( self.SMarginal)) def _preSetupApproxPivoted(self, mus:paramList) \ -> Tuple[ListAny, ListAny, ListAny]: self.computeScaleFactor() if self.trainedModel is None: self._setupTrainedModel(np.zeros((0, 0))) self.trainedModel.data.Qs, self.trainedModel.data.Ps = [], [] self.trainedModel.data.Psupp = [] self._trainedModelOld = copy(self.trainedModel) self._scaleFactorOldPivot = copy(self.scaleFactor) self.scaleFactor = self.scaleFactorPivot self._temporaryPivot = 1 self._musLoc = copy(self.mus) idx, sizes = indicesScatter(len(mus), return_sizes = True) emptyCores = np.where(np.logical_not(sizes))[0] self.verbosity -= 15 return idx, sizes, emptyCores def _postSetupApproxPivoted(self, mus:Np2D, pMat:Np2D, Ps:ListAny, Qs:ListAny, sizes:ListAny): self.scaleFactor = self._scaleFactorOldPivot del self._scaleFactorOldPivot, self._temporaryPivot pMat, Ps, Qs, mus, nsamples = gatherPivotedApproximant(pMat, Ps, Qs, mus, sizes, self.polybasis) if len(self._musLoc) > 0: self._mus = self.checkParameterList(self._musLoc) self._mus.append(mus) else: self._mus = self.checkParameterList(mus) self.trainedModel = self._trainedModelOld del self._trainedModelOld padLeft = self.trainedModel.data.projMat.shape[1] suppNew = np.append(0, np.cumsum(nsamples)) self._setupTrainedModel(pMat, padLeft > 0) self.trainedModel.data.Qs += Qs self.trainedModel.data.Ps += Ps self.trainedModel.data.Psupp += list(padLeft + suppNew[: -1]) self.trainedModel.data.approxParameters = copy(self.approxParameters) self.verbosity += 15 def _localPivotedResult(self, pMat:Np2D, req:ListAny, emptyCores:ListAny, mus:Np2D) -> Tuple[Np2D, ListAny, Np2D]: if pMat is None: mus = copy(self.samplingEngine.mus.data) pMat = copy(self.samplingEngine.projectionMatrix) if masterCore(): for dest in emptyCores: req += [isend((len(pMat), pMat.dtype, mus.dtype), dest = dest, tag = dest)] else: mus = np.vstack((mus, self.samplingEngine.mus.data)) pMat = np.hstack((pMat, self.samplingEngine.projectionMatrix)) return pMat, req, mus @abstractmethod def setupApproxPivoted(self, mus:paramList) -> int: if self.checkComputedApproxPivoted(): return -1 RROMPyAssert(self._mode, message = "Cannot setup approximant.") vbMng(self, "INIT", "Setting up pivoted approximant.", 10) self._preSetupApproxPivoted() data = [] pass self._postSetupApproxPivoted(mus, data) vbMng(self, "DEL", "Done setting up pivoted approximant.", 10) return 0 def setupApprox(self, plotEst : str = "NONE") -> int: """Compute greedy snapshots of solution map.""" if self.checkComputedApprox(): return -1 RROMPyAssert(self._mode, message = "Cannot start greedy algorithm.") vbMng(self, "INIT", "Starting computation of snapshots.", 3) max2ErrorEst, self.firstGreedyIterM = np.inf, True self._preliminaryTrainingMarginal() if self.errorEstimatorKindMarginal[: 10] == "LOOK_AHEAD": muidx = np.arange(len(self.trainedModel.data.musMarginal)) else:#if self.errorEstimatorKindMarginal in ["LEAVE_ONE_OUT", "NONE"]: muidx = [] self._musMarginalTestIdxs = np.array(muidx) while self.firstGreedyIterM or (max2ErrorEst > self.greedyTolMarginal and self.samplerMarginal.npoints < self.maxIterMarginal): errorEstTest, muidx, maxErrorEst, mu = \ self.greedyNextSampleMarginal(muidx, plotEst) if maxErrorEst is None: max2ErrorEst = 1. + self.greedyTolMarginal else: if len(maxErrorEst) > 0: max2ErrorEst = np.max(maxErrorEst) vbMng(self, "MAIN", ("Uniform testing error estimate " "{:.4e}.").format(max2ErrorEst), 3) else: max2ErrorEst = 0. if plotEst == "LAST": self.plotEstimatorMarginal(errorEstTest, muidx, maxErrorEst) vbMng(self, "DEL", ("Done computing snapshots (final snapshot count: " "{}).").format(len(self.mus)), 3) if (self.errorEstimatorKindMarginal == "LOOK_AHEAD_RECOVER" and hasattr(self.trainedModel, "_idxExcl") and len(self.trainedModel._idxExcl) > 0): vbMng(self, "INIT", "Recovering {} test models.".format( len(self.trainedModel._idxExcl)), 7) for j, mu in zip(self.trainedModel._idxExcl, self.trainedModel._musMExcl): self.musMarginal.insert(mu, j) self._updateTrainedModelMarginalSamples() self._finalizeMarginalization() self._SMarginal = len(self.musMarginal) self._approxParameters["SMarginal"] = self.SMarginal self.trainedModel.data.approxParameters["SMarginal"] = ( self.SMarginal) vbMng(self, "DEL", "Done recovering test models.", 7) return 0 def checkComputedApproxPivoted(self) -> bool: return (super().checkComputedApprox() and len(self.musMarginal) == len(self.trainedModel.data.musMarginal)) class GenericPivotedGreedyApproximantNoMatch( GenericPivotedGreedyApproximantBase, GenericPivotedApproximantNoMatch): """ ROM pivoted greedy interpolant computation for parametric problems (without pole matching) (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; - 'scaleFactorDer': scaling factors for derivative computation; defaults to 'AUTO'; - 'cutOffTolerance': tolerance for ignoring parasitic poles; defaults to np.inf; - 'matchingWeightError': weight for pole matching optimization in error estimation; defaults to 0; - 'cutOffToleranceError': tolerance for ignoring parasitic poles in error estimation; defaults to 'AUTO', i.e. cutOffTolerance; - 'S': total number of pivot samples current approximant relies upon; - 'samplerPivot': pivot sample point generator; - 'SMarginal': number of starting marginal samples; - 'samplerMarginal': marginal sample point generator via sparse grid; - 'errorEstimatorKindMarginal': kind of marginal error estimator; available values include 'LEAVE_ONE_OUT', 'LOOK_AHEAD', 'LOOK_AHEAD_RECOVER', and 'NONE'; defaults to 'NONE'; - 'greedyTolMarginal': uniform error tolerance for marginal greedy algorithm; defaults to 1e-1; - 'maxIterMarginal': maximum number of marginal greedy steps; defaults to 1e2; - 'radialDirectionalWeightsMarginal': radial basis weights for marginal interpolant; defaults to 1. 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. directionPivot: Pivot components. mus: Array of snapshot parameters. musMarginal: Array of marginal 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; - 'scaleFactorDer': scaling factors for derivative computation; - 'cutOffTolerance': tolerance for ignoring parasitic poles; - 'matchingWeightError': weight for pole matching optimization in error estimation; - 'cutOffToleranceError': tolerance for ignoring parasitic poles in error estimation; - 'errorEstimatorKindMarginal': kind of marginal error estimator; - 'greedyTolMarginal': uniform error tolerance for marginal greedy algorithm; - 'maxIterMarginal': maximum number of marginal greedy steps; - 'radialDirectionalWeightsMarginal': radial basis weights for marginal interpolant. 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 via sparse grid. approx_state: Whether to approximate state. verbosity: Verbosity level. POD: Whether to compute POD of snapshots. scaleFactorDer: Scaling factors for derivative computation. cutOffTolerance: Tolerance for ignoring parasitic poles. matchingWeightError: Weight for pole matching optimization in error estimation. cutOffToleranceError: Tolerance for ignoring parasitic poles in error estimation. 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 via sparse grid. errorEstimatorKindMarginal: Kind of marginal error estimator. greedyTolMarginal: Uniform error tolerance for marginal greedy algorithm. maxIterMarginal: Maximum number of marginal greedy steps. radialDirectionalWeightsMarginal: Radial basis weights for marginal interpolant. muBounds: 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 _poleMatching(self): vbMng(self, "INIT", "Compressing poles.", 10) self.trainedModel.initializeFromRational() vbMng(self, "DEL", "Done compressing poles.", 10) def _updateTrainedModelMarginalSamples(self, idx : ListAny = []): self.trainedModel.updateEffectiveSamples(idx) class GenericPivotedGreedyApproximant(GenericPivotedGreedyApproximantBase, GenericPivotedApproximant): """ ROM pivoted greedy interpolant computation for parametric problems (with pole matching) (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; - 'scaleFactorDer': scaling factors for derivative computation; defaults to 'AUTO'; - 'matchingWeight': weight for pole matching optimization; defaults to 1; - 'matchingMode': mode for pole matching optimization; allowed values include 'NONE' and 'SHIFT'; defaults to 'NONE'; - 'cutOffTolerance': tolerance for ignoring parasitic poles; defaults to np.inf; - 'cutOffSharedRatio': required ratio of marginal points to share resonance in cut off strategy; defaults to 1.; - 'matchingWeightError': weight for pole matching optimization in error estimation; defaults to 0; - 'cutOffToleranceError': tolerance for ignoring parasitic poles in error estimation; defaults to 'AUTO', i.e. cutOffTolerance; - 'S': total number of pivot samples current approximant relies upon; - 'samplerPivot': pivot sample point generator; - 'SMarginal': number of starting marginal samples; - 'samplerMarginal': marginal sample point generator via sparse grid; - 'errorEstimatorKindMarginal': kind of marginal error estimator; available values include 'LEAVE_ONE_OUT', 'LOOK_AHEAD', 'LOOK_AHEAD_RECOVER', and 'NONE'; defaults to 'NONE'; - 'polybasisMarginal': type of polynomial basis for marginal interpolation; allowed values include 'MONOMIAL_*', 'CHEBYSHEV_*', 'LEGENDRE_*', 'NEARESTNEIGHBOR', and 'PIECEWISE_LINEAR_*'; defaults to 'MONOMIAL'; - 'paramsMarginal': dictionary of parameters for marginal interpolation; include: . 'MMarginal': degree of marginal interpolant; defaults to 'AUTO', i.e. maximum allowed; not for 'NEARESTNEIGHBOR' or 'PIECEWISE_LINEAR_*'; . 'nNeighborsMarginal': number of marginal nearest neighbors; defaults to 1; only for 'NEARESTNEIGHBOR'; . 'polydegreetypeMarginal': type of polynomial degree for marginal; defaults to 'TOTAL'; not for 'NEARESTNEIGHBOR' or 'PIECEWISE_LINEAR_*'; . 'interpRcondMarginal': tolerance for marginal interpolation; defaults to None; not for 'NEARESTNEIGHBOR'; . 'radialDirectionalWeightsMarginalAdapt': bounds for adaptive rescaling of marginal radial basis weights; only for radial basis. - 'greedyTolMarginal': uniform error tolerance for marginal greedy algorithm; defaults to 1e-1; - 'maxIterMarginal': maximum number of marginal greedy steps; defaults to 1e2; - 'radialDirectionalWeightsMarginal': radial basis weights for marginal interpolant; defaults to 1. 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. directionPivot: Pivot components. mus: Array of snapshot parameters. musMarginal: Array of marginal 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; - 'scaleFactorDer': scaling factors for derivative computation; - 'matchingWeight': weight for pole matching optimization; - 'matchingMode': mode for pole matching optimization; - 'cutOffTolerance': tolerance for ignoring parasitic poles; - 'cutOffSharedRatio': required ratio of marginal points to share resonance in cut off strategy; - 'matchingWeightError': weight for pole matching optimization in error estimation; - 'cutOffToleranceError': tolerance for ignoring parasitic poles in error estimation; - 'errorEstimatorKindMarginal': kind of marginal error estimator; - 'polybasisMarginal': type of polynomial basis for marginal interpolation; - 'paramsMarginal': dictionary of parameters for marginal interpolation; include: . 'MMarginal': degree of marginal interpolant; . 'nNeighborsMarginal': number of marginal nearest neighbors; . 'polydegreetypeMarginal': type of polynomial degree for marginal; . 'interpRcondMarginal': tolerance for marginal interpolation; . 'radialDirectionalWeightsMarginalAdapt': bounds for adaptive rescaling of marginal radial basis weights. - 'greedyTolMarginal': uniform error tolerance for marginal greedy algorithm; - 'maxIterMarginal': maximum number of marginal greedy steps; - 'radialDirectionalWeightsMarginal': radial basis weights for marginal interpolant. 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 via sparse grid. approx_state: Whether to approximate state. verbosity: Verbosity level. POD: Whether to compute POD of snapshots. scaleFactorDer: Scaling factors for derivative computation. matchingWeight: Weight for pole matching optimization. matchingMode: Mode for pole matching optimization. cutOffTolerance: Tolerance for ignoring parasitic poles. cutOffSharedRatio: Required ratio of marginal points to share resonance in cut off strategy. matchingWeightError: Weight for pole matching optimization in error estimation. cutOffToleranceError: Tolerance for ignoring parasitic poles in error estimation. 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 via sparse grid. errorEstimatorKindMarginal: Kind of marginal error estimator. polybasisMarginal: Type of polynomial basis for marginal interpolation. paramsMarginal: Dictionary of parameters for marginal interpolation. greedyTolMarginal: Uniform error tolerance for marginal greedy algorithm. maxIterMarginal: Maximum number of marginal greedy steps. radialDirectionalWeightsMarginal: Radial basis weights for marginal interpolant. muBounds: 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 _poleMatching(self): vbMng(self, "INIT", "Compressing and matching poles.", 10) self.trainedModel.initializeFromRational(self.matchingWeight, self.matchingMode, self.HFEngine, False) vbMng(self, "DEL", "Done compressing and matching poles.", 10) def _updateTrainedModelMarginalSamples(self, idx : ListAny = []): self.trainedModel.updateEffectiveSamples(idx, self.matchingWeight, self.matchingMode, self.HFEngine, False) def getErrorEstimatorMarginalLeaveOneOut(self) -> Np1D: if self.polybasisMarginal != "NEARESTNEIGHBOR": if not hasattr(self, "_MMarginal_isauto"): if not hasattr(self, "_MMarginalOriginal"): self._MMarginalOriginal = self.paramsMarginal["MMarginal"] self.paramsMarginal["MMarginal"] = self._MMarginalOriginal self._reduceDegreeNNoWarn = 1 err = super().getErrorEstimatorMarginalLeaveOneOut() if self.polybasisMarginal != "NEARESTNEIGHBOR": del self._reduceDegreeNNoWarn return err def setupApprox(self, *args, **kwargs) -> int: if self.checkComputedApprox(): return -1 self.purgeparamsMarginal() return super().setupApprox(*args, **kwargs) diff --git a/rrompy/reduction_methods/standard/greedy/generic_greedy_approximant.py b/rrompy/reduction_methods/standard/greedy/generic_greedy_approximant.py index aee4fd9..c37383b 100644 --- a/rrompy/reduction_methods/standard/greedy/generic_greedy_approximant.py +++ b/rrompy/reduction_methods/standard/greedy/generic_greedy_approximant.py @@ -1,650 +1,651 @@ # 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 . # from abc import abstractmethod from copy import deepcopy as copy import numpy as np from matplotlib import pyplot as plt from rrompy.hfengines.base.linear_affine_engine import checkIfAffine from rrompy.reduction_methods.standard.generic_standard_approximant import ( GenericStandardApproximant) from rrompy.utilities.base.types import (Np1D, Np2D, 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.sampling.sample_list import sampleList from rrompy.parameter import emptyParameterList, parameterList from rrompy.utilities.parallel import masterCore __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 isinstance(mus, (parameterList, sampleList)): mus = mus.data if isinstance(badmus, (parameterList, sampleList)): badmus = badmus.data if len(badmus) == 0: return np.arange(len(mus)) proximity = np.min(localL2Distance(mus, badmus), axis = 1) return np.where(proximity <= tol)[0] 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; - 'scaleFactorDer': scaling factors for derivative computation; defaults to 'AUTO'; - '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.; - 'maxIter': maximum number of greedy steps; defaults to 1e2; - '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; - 'scaleFactorDer': scaling factors for derivative computation; - 'greedyTol': uniform error tolerance for greedy algorithm; - 'collinearityTol': collinearity tolerance for greedy algorithm; - 'maxIter': maximum number of greedy steps; - '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. scaleFactorDer: Scaling factors for derivative computation. S: number of test points. sampler: Sample point generator. greedyTol: Uniform error tolerance for greedy algorithm. collinearityTol: Collinearity tolerance for greedy algorithm. maxIter: maximum number of greedy steps. 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. """ def __init__(self, *args, **kwargs): self._preInit() self._addParametersToList(["greedyTol", "collinearityTol", "maxIter", "nTestPoints"], [1e-2, 0., 1e2, 5e2], ["trainSetGenerator"], ["AUTO"]) super().__init__(*args, **kwargs) 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 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 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 self.approx_state: if not hasattr(self.HFEngine, "energyNormDualMatrix"): self.HFEngine.buildEnergyNormDualForm() estimatorEnergyMatrix = self.HFEngine.energyNormDualMatrix else: estimatorEnergyMatrix = self.HFEngine.outputNormMatrix else: if hasattr(normEngn, "buildEnergyNormDualForm"): if not hasattr(normEngn, "energyNormDualMatrix"): normEngn.buildEnergyNormDualForm() estimatorEnergyMatrix = normEngn.energyNormDualMatrix else: estimatorEnergyMatrix = normEngn self.estimatorNormEngine = normEngine(estimatorEnergyMatrix) def _affineResidualMatricesContraction(self, rb:Np2D, rA : Np2D = None) \ -> Tuple[Np1D, Np1D, Np1D]: self.assembleReducedResidualBlocks(full = rA is not None) # '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 getErrorEstimatorAffine(self, mus:Np1D) -> Np1D: """Standard residual estimator.""" checkIfAffine(self.HFEngine, "apply affinity-based error estimator") self.HFEngine.buildA() self.HFEngine.buildb() mus = self.checkParameterList(mus) tMverb, self.trainedModel.verbosity = self.trainedModel.verbosity, 0 uApproxRs = self.getApproxReduced(mus).data self.trainedModel.verbosity = tMverb muTestEff = self.HFEngine.mapParameterList(mus) 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, 1) * radiusA ff, Lf, LL = self._affineResidualMatricesContraction(radiusb, radiusA) err = np.abs((LL - 2. * np.real(Lf) + ff) / ff) ** .5 return err def errorEstimator(self, mus:Np1D, return_max : bool = False) -> Np1D: setupOK = self.setupApproxLocal() if setupOK > 0: err = np.empty(len(mus)) err[:] = np.nan if not return_max: return err return err, [- setupOK], np.nan mus = self.checkParameterList(mus) vbMng(self.trainedModel, "INIT", "Evaluating error estimator at mu = {}.".format(mus), 10) err = self.getErrorEstimatorAffine(mus) vbMng(self.trainedModel, "DEL", "Done evaluating error estimator", 10) if not return_max: return err idxMaxEst = [np.argmax(err)] return err, idxMaxEst, err[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 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) cLevel = np.inf if np.isclose(cLevel, 0.) else cLevel ** -1. vbMng(self, "MAIN", "Collinearity indicator {:.4e}.".format(cLevel), 3) return cLevel > self.collinearityTol def plotEstimator(self, est:Np1D, idxMax:List[int], estMax:List[float]): if (not (np.any(np.isnan(est)) or np.any(np.isinf(est))) and masterCore()): fig = plt.figure(figsize = plt.figaspect(1. / self.npar)) for jpar in range(self.npar): ax = fig.add_subplot(1, self.npar, 1 + jpar) musre = np.array(self.muTest.re.data) errCP = copy(est) idx = np.delete(np.arange(self.npar), jpar) while len(musre) > 0: if self.npar == 1: currIdx = np.arange(len(musre)) else: currIdx = np.where(np.isclose(np.sum( np.abs(musre[:, idx] - musre[0, idx]), 1), 0.))[0] ax.semilogy(musre[currIdx, jpar], errCP[currIdx], 'k', linewidth = 1) musre = np.delete(musre, currIdx, 0) errCP = np.delete(errCP, currIdx) ax.semilogy([self.muBounds.re(0, jpar), self.muBounds.re(-1, jpar)], [self.greedyTol] * 2, 'r--') ax.semilogy(self.mus.re(jpar), 2. * self.greedyTol * np.ones(len(self.mus)), '*m') if len(idxMax) > 0 and estMax is not None: ax.semilogy(self.muTest.re(idxMax, jpar), estMax, 'xr') + ax.set_xlim(*list(self.sampler.lims.re(jpar))) ax.grid() plt.tight_layout() plt.show() def greedyNextSample(self, muidx:int, plotEst : str = "NONE")\ -> 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 j, mu in enumerate(mus): vbMng(self, "MAIN", ("Adding sample point no. {} at {} to training " "set.").format(len(self.mus) + 1, mu), 3) self.mus.append(mu) self._S = len(self.mus) self._approxParameters["S"] = self.S if (self.samplingEngine.nsamples <= len(mus) - j - 1 or not np.allclose(mu, self.samplingEngine.mus[j - len(mus)])): self.samplingEngine.nextSample(mu) if self._isLastSampleCollinear(): vbMng(self, "MAIN", ("Collinearity above tolerance detected. Starting " "preemptive greedy loop termination."), 3) self._collinearityFlag = 1 errorEstTest = np.empty(len(self.muTest)) errorEstTest[:] = np.nan return errorEstTest, [-1], np.nan, np.nan errorEstTest, muidx, maxErrorEst = self.errorEstimator(self.muTest, True) if plotEst == "ALL": self.plotEstimator(errorEstTest, muidx, maxErrorEst) 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.resetSamples() self.computeScaleFactor() self.samplingEngine.scaleFactor = self.scaleFactorDer self.mus = self.trainSetGenerator.generatePoints(self.S) while len(self.mus) > self.S: self.mus.pop() muTestBase = self.sampler.generatePoints(self.nTestPoints, False) idxPop = pruneSamples(self.HFEngine.mapParameterList(muTestBase), self.HFEngine.mapParameterList(self.mus), 1e-10 * self.scaleFactor[0]) muTestBase.pop(idxPop) 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), 3) 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 @abstractmethod def setupApproxLocal(self) -> int: if self.checkComputedApprox(): return -1 RROMPyAssert(self._mode, message = "Cannot setup approximant.") vbMng(self, "INIT", "Setting up local approximant.", 5) pass vbMng(self, "DEL", "Done setting up local approximant.", 5) return 0 def setupApprox(self, plotEst : str = "NONE") -> int: """Compute greedy snapshots of solution map.""" if self.checkComputedApprox(): return -1 RROMPyAssert(self._mode, message = "Cannot start greedy algorithm.") vbMng(self, "INIT", "Starting computation of snapshots.", 3) self._collinearityFlag = 0 self._preliminaryTraining() muidx, self.firstGreedyIter = [len(self.muTest) - 1], True errorEstTest, maxErrorEst = [np.inf], np.inf max2ErrorEst, trainedModelOld = np.inf, None while self.firstGreedyIter or (len(self.muTest) > 0 and (maxErrorEst is None or max2ErrorEst > self.greedyTol) and self.samplingEngine.nsamples < self.maxIter): muTestOld, errorEstTestOld = self.muTest, errorEstTest muidxOld, maxErrorEstOld = muidx, maxErrorEst errorEstTest, muidx, maxErrorEst, mu = self.greedyNextSample( muidx, plotEst) if maxErrorEst is not None and (np.any(np.isnan(maxErrorEst)) or np.any(np.isinf(maxErrorEst))): if self._collinearityFlag == 0 and not self.firstGreedyIter: RROMPyWarning(("Instability in a posteriori " "estimator. Starting preemptive greedy " "loop termination.")) self.muTest, errorEstTest = muTestOld, errorEstTestOld if self.firstGreedyIter: self.mus.pop(-1) self.samplingEngine.popSample() if muidx[0] < 0: self.trainedModel = None raise RROMPyException(("Instability in approximant " "computation. Aborting greedy " "iterations.")) else: self._approxParameters = ( trainedModelOld.data.approxParameters) self._S = trainedModelOld.data.approxParameters["S"] self._approxParameters["S"] = self.S self.trainedModel.data = copy(trainedModelOld.data) muidx, maxErrorEst = muidxOld, maxErrorEstOld break if maxErrorEst is not None: max2ErrorEst = np.max(maxErrorEst) vbMng(self, "MAIN", ("Uniform testing error estimate " "{:.4e}.").format(max2ErrorEst), 3) if self.firstGreedyIter: trainedModelOld = copy(self.trainedModel) else: trainedModelOld.data = copy(self.trainedModel.data) self.firstGreedyIter = False if (maxErrorEst is None or max2ErrorEst <= self.greedyTol or np.any(np.isnan(maxErrorEst)) or np.any(np.isinf(maxErrorEst))): while self.samplingEngine.nsamples > self.S: self.samplingEngine.popSample() while len(self.mus) > self.S: self.mus.pop(-1) else: self._S = self.samplingEngine.nsamples self._approxParameters["S"] = self.S while len(self.mus) < self.S: self.mus.append(self.samplingEngine.mus[len(self.mus)]) self.setupApproxLocal() if plotEst == "LAST": self.plotEstimator(errorEstTest, muidx, maxErrorEst) vbMng(self, "DEL", ("Done computing snapshots (final snapshot count: " "{}).").format(self.samplingEngine.nsamples), 3) return 0 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 isinstance(pMat, (parameterList, sampleList)): 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 isinstance(pMat, (parameterList, sampleList)): 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 isinstance(pMat, (parameterList, sampleList)): 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 isinstance(pMat, (parameterList, sampleList)): 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.""" if full: checkIfAffine(self.HFEngine, "assemble reduced residual blocks") else: checkIfAffine(self.HFEngine, "assemble reduced RHS blocks", True) self.HFEngine.buildb() self.assembleReducedResidualBlocksbb(self.HFEngine.bs) if full: pMat = self.samplingEngine.projectionMatrix self.HFEngine.buildA() self.assembleReducedResidualBlocksAb(self.HFEngine.As, self.HFEngine.bs, pMat) self.assembleReducedResidualBlocksAA(self.HFEngine.As, pMat)