diff --git a/rrompy/parameter/parameter_sampling/fft_sampler.py b/rrompy/parameter/parameter_sampling/fft_sampler.py index 5ea7e20..4199886 100644 --- a/rrompy/parameter/parameter_sampling/fft_sampler.py +++ b/rrompy/parameter/parameter_sampling/fft_sampler.py @@ -1,50 +1,48 @@ # 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 .generic_sampler import GenericSampler from rrompy.utilities.base.types import List, paramList from rrompy.utilities.numerical import lowDiscrepancy, kroneckerer from rrompy.parameter import checkParameterList __all__ = ['FFTSampler'] class FFTSampler(GenericSampler): """Generator of FFT-type sample points on scaled roots of unity.""" def generatePoints(self, n:List[int], reorder : bool = True) -> paramList: """Array of sample points.""" if not hasattr(n, "__len__"): n = [n] super().generatePoints(n) nleft, nright = 1, np.prod(n) xmat = np.empty((nright, self.npar), dtype = np.complex) for d in range(self.npar): nright //= n[d] - a, b = self.lims(0, d), self.lims(1, d) - if self.scaling is not None: - a, b = self.scaling[d](a), self.scaling[d](b) + a = self.lims(0, d) ** self.scalingExp[d] + b = self.lims(1, d) ** self.scalingExp[d] c, r = (a + b) / 2., np.abs(a - b) / 2. xd = c + r * np.exp(1.j * np.linspace(0, 2 * np.pi, n[d] + 1)[:-1]) - if self.scalingInv is not None: - xd = self.scalingInv[d](xd) + xd **= 1. / self.scalingExp[d] xmat[:, d] = kroneckerer(xd, nleft, nright) nleft *= n[d] if reorder: xmat = xmat[lowDiscrepancy(np.prod(n)), :] x = checkParameterList(xmat, self.npar)[0] return x diff --git a/rrompy/parameter/parameter_sampling/generic_sampler.py b/rrompy/parameter/parameter_sampling/generic_sampler.py index 9282e18..30f631a 100644 --- a/rrompy/parameter/parameter_sampling/generic_sampler.py +++ b/rrompy/parameter/parameter_sampling/generic_sampler.py @@ -1,98 +1,79 @@ # 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 rrompy.utilities.base.types import List, paramList from rrompy.utilities.exception_manager import RROMPyException, RROMPyAssert from rrompy.parameter import checkParameterList __all__ = ['GenericSampler'] class GenericSampler: """ABSTRACT. Generic generator of sample points.""" - def __init__(self, lims:paramList, scaling : List[callable] = None, - scalingInv : List[callable] = None): + def __init__(self, lims:paramList, scalingExp : List[float] = None): self.lims = lims - self.scaling = scaling - self.scalingInv = scalingInv + self.scalingExp = scalingExp def name(self) -> str: return self.__class__.__name__ def __str__(self) -> str: return "{}[{}_{}]".format(self.name(), self.lims[0], self.lims[1]) def __repr__(self) -> str: return self.__str__() + " at " + hex(id(self)) def __eq__(self, other) -> bool: return self.__dict__ == other.__dict__ @property def npar(self): """Number of parameters.""" return self._lims.shape[1] @property def lims(self): """Value of lims.""" return self._lims @lims.setter def lims(self, lims): lims = checkParameterList(lims)[0] if len(lims) != 2: raise RROMPyException("2 limits must be specified.") self._lims = lims @property - def scaling(self): - """Value of scaling.""" - return self._scaling - @scaling.setter - def scaling(self, scaling): - if scaling is not None: - if not hasattr(scaling, "__len__"): scaling = [scaling] - RROMPyAssert(self.npar, len(scaling), "Number of scaling terms") - if not all([callable(s) for s in scaling]): - raise RROMPyException(("Each value of scaling must be a " - "callable.")) - self._scaling = scaling - - @property - def scalingInv(self): - """Value of scalingInv.""" - return self._scalingInv - @scalingInv.setter - def scalingInv(self, scalingInv): - if scalingInv is not None: - if not hasattr(scalingInv, "__len__"): scalingInv = [scalingInv] - RROMPyAssert(self.npar, len(scalingInv), - "Number of scalingInv terms") - if not all([callable(sInv) for sInv in scalingInv]): - raise RROMPyException(("Each value of scalingInv must be a " - "callable.")) - self._scalingInv = scalingInv + def scalingExp(self): + """Value of scalingExp.""" + return self._scalingExp + @scalingExp.setter + def scalingExp(self, scalingExp): + if scalingExp is None: + scalingExp = [1.] * self.npar + if not hasattr(scalingExp, "__len__"): scalingExp = [scalingExp] + RROMPyAssert(self.npar, len(scalingExp), "Number of scaling terms") + self._scalingExp = scalingExp @abstractmethod def generatePoints(self, n:List[int]) -> paramList: """Array of points.""" if not hasattr(n, "__len__"): n = [n] RROMPyAssert(self.npar, len(n), "Point number") pass diff --git a/rrompy/parameter/parameter_sampling/manual_sampler.py b/rrompy/parameter/parameter_sampling/manual_sampler.py index 4f4be17..42aa93d 100644 --- a/rrompy/parameter/parameter_sampling/manual_sampler.py +++ b/rrompy/parameter/parameter_sampling/manual_sampler.py @@ -1,63 +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, paramList from rrompy.parameter import checkParameterList __all__ = ['ManualSampler'] class ManualSampler(GenericSampler): """Manual generator of sample points.""" def __init__(self, lims:paramList, points:paramList, - scaling : List[callable] = None, - scalingInv : List[callable] = None): - super().__init__(lims = lims, scaling = scaling, - scalingInv = scalingInv) + scalingExp : List[float] = None): + super().__init__(lims = lims, scalingExp = scalingExp) self.points = points @property def points(self): """Value of points.""" return self._points @points.setter def points(self, points): points = checkParameterList(points, self.npar)[0] 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) -> paramList: """Array of sample points.""" if hasattr(n, "__len__"): n = n[0] 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)[0] return x diff --git a/rrompy/parameter/parameter_sampling/quadrature_sampler.py b/rrompy/parameter/parameter_sampling/quadrature_sampler.py index 5fb7ca9..a08adf8 100644 --- a/rrompy/parameter/parameter_sampling/quadrature_sampler.py +++ b/rrompy/parameter/parameter_sampling/quadrature_sampler.py @@ -1,87 +1,83 @@ # 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 .generic_sampler import GenericSampler from rrompy.utilities.base.types import List, paramList from rrompy.utilities.exception_manager import RROMPyException from rrompy.utilities.numerical import lowDiscrepancy, kroneckerer from rrompy.parameter import checkParameterList __all__ = ['QuadratureSampler'] class QuadratureSampler(GenericSampler): """Generator of quadrature sample points.""" allowedKinds = ["UNIFORM", "CHEBYSHEV", "GAUSSLEGENDRE"] def __init__(self, lims:paramList, kind : str = "UNIFORM", - scaling : List[callable] = None, - scalingInv : List[callable] = None): - super().__init__(lims = lims, scaling = scaling, - scalingInv = scalingInv) + scalingExp : List[float] = None): + super().__init__(lims = lims, scalingExp = scalingExp) 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): if kind.upper() not in self.allowedKinds: raise RROMPyException("Generator kind not recognized.") self._kind = kind.upper() def generatePoints(self, n:List[int], reorder : bool = True) -> paramList: """Array of sample points.""" if not hasattr(n, "__len__"): n = [n] super().generatePoints(n) nleft, nright = 1, np.prod(n) xmat = np.empty((nright, self.npar), dtype = self.lims.dtype) for d in range(self.npar): nright //= n[d] - a, b = self.lims(0, d), self.lims(1, d) - if self.scaling is not None: - a, b = self.scaling[d](a), self.scaling[d](b) + a = self.lims(0, d) ** self.scalingExp[d] + b = self.lims(1, d) ** self.scalingExp[d] c, r = (a + b) / 2., (a - b) / 2. dAbs = 2. * np.abs(r) if self.kind == "UNIFORM": xd = np.linspace(a, b, n[d]) elif self.kind == "CHEBYSHEV": nodes, _ = np.polynomial.chebyshev.chebgauss(n[d]) xd = c + r * nodes elif self.kind == "GAUSSLEGENDRE": nodes, _ = np.polynomial.legendre.leggauss(n[d]) xd = c + r * nodes[::-1] - if self.scalingInv is not None: - xd = self.scalingInv[d](xd) + xd **= 1. / self.scalingExp[d] xmat[:, d] = kroneckerer(xd, nleft, nright) nleft *= n[d] nright = np.prod(n) if nright > 1 and reorder: fejerOrdering = [nright - 1] + lowDiscrepancy(nright - 1) xmat = xmat[fejerOrdering, :] x = checkParameterList(xmat, self.npar)[0] return x diff --git a/rrompy/parameter/parameter_sampling/random_sampler.py b/rrompy/parameter/parameter_sampling/random_sampler.py index 460d4e8..68625d2 100644 --- a/rrompy/parameter/parameter_sampling/random_sampler.py +++ b/rrompy/parameter/parameter_sampling/random_sampler.py @@ -1,76 +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 . # import numpy as np from .generic_sampler import GenericSampler from rrompy.utilities.numerical import haltonGenerate, sobolGenerate from rrompy.utilities.base.types import List, paramList from rrompy.utilities.exception_manager import RROMPyException from rrompy.parameter import checkParameterList __all__ = ['RandomSampler'] class RandomSampler(GenericSampler): """Generator of quadrature sample points.""" allowedKinds = ["UNIFORM", "HALTON", "SOBOL"] def __init__(self, lims:paramList, kind : str = "UNIFORM", - scaling : List[callable] = None, - scalingInv : List[callable] = None, seed : int = 42): - super().__init__(lims = lims, scaling = scaling, - scalingInv = scalingInv) + scalingExp : List[float] = None, seed : int = 42): + super().__init__(lims = lims, scalingExp = scalingExp) 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() def generatePoints(self, n:int) -> paramList: """Array of quadrature points.""" if hasattr(n, "__len__"): n = n[0] 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) else: xmat = sobolGenerate(self.npar, n, self.seed) for d in range(self.npar): - a, b = self.lims(0, d), self.lims(1, d) - if self.scaling is not None: - a, b = self.scaling[d](a), self.scaling[d](b) + a = self.lims(0, d) ** self.scalingExp[d] + b = self.lims(1, d) ** self.scalingExp[d] xmat[:, d] = a + (b - a) * xmat[:, d] - if self.scalingInv is not None: - xmat[:, d] = self.scalingInv[d](xmat[:, d]) + xmat[:, d] **= 1. / self.scalingExp[d] x = checkParameterList(xmat, self.npar)[0] return x diff --git a/rrompy/reduction_methods/base/__init__.py b/rrompy/reduction_methods/base/__init__.py index 98d784c..f6aa877 100644 --- a/rrompy/reduction_methods/base/__init__.py +++ b/rrompy/reduction_methods/base/__init__.py @@ -1,29 +1,29 @@ # 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_approximant import GenericApproximant -from .pade_utils import checkRobustTolerance -from .rb_utils import projectAffineDecomposition +from .rational_interpolant_utils import checkRobustTolerance +from .reduced_basis_utils import projectAffineDecomposition __all__ = [ 'GenericApproximant', 'checkRobustTolerance', 'projectAffineDecomposition' ] diff --git a/rrompy/reduction_methods/base/generic_approximant.py b/rrompy/reduction_methods/base/generic_approximant.py index 53ca38f..294f1ed 100644 --- a/rrompy/reduction_methods/base/generic_approximant.py +++ b/rrompy/reduction_methods/base/generic_approximant.py @@ -1,868 +1,868 @@ # 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 itertools import product as iterprod from copy import deepcopy as copy from os import remove as osrm from rrompy.sampling.standard import (SamplingEngineStandard, SamplingEngineStandardPOD) from rrompy.utilities.base.types import (Np1D, DictAny, HFEng, List, Tuple, ListAny, strLst, paramVal, paramList, sampList) from rrompy.utilities.base import (purgeDict, verbosityManager as vbMng, getNewFilename) from rrompy.utilities.exception_manager import (RROMPyException, RROMPyAssert, RROMPy_READY, RROMPy_FRAGILE) from rrompy.utilities.base import pickleDump, pickleLoad from rrompy.parameter import (emptyParameterList, checkParameter, checkParameterList) from rrompy.sampling import sampleList, emptySampleList __all__ = ['GenericApproximant'] def addNormFieldToClass(self, fieldName): def objFunc(self, mu:paramList, homogeneized : bool = False) -> Np1D: uV = getattr(self.__class__, "get" + fieldName)(self, mu, homogeneized) val = self.HFEngine.norm(uV) return val setattr(self.__class__, "norm" + fieldName, objFunc) def addPlotFieldToClass(self, fieldName): def objFunc(self, mu:paramList, *args, homogeneized : bool = False, **kwargs): uV = getattr(self.__class__, "get" + fieldName)(self, mu, homogeneized) kwargsCopy = copy(kwargs) for j, u in enumerate(uV): if "name" in kwargs.keys(): kwargsCopy["name"] = kwargs["name"] + str(j) self.HFEngine.plot(u, *args, **kwargs) setattr(self.__class__, "plot" + fieldName, objFunc) def addOutParaviewFieldToClass(self, fieldName): def objFunc(self, mu:paramVal, *args, homogeneized : bool = False, **kwargs): if not hasattr(self.HFEngine, "outParaview"): raise RROMPyException(("High fidelity engine cannot output to " "Paraview.")) uV = getattr(self.__class__, "get" + fieldName)(self, mu, homogeneized) kwargsCopy = copy(kwargs) for j, u in enumerate(uV): if "name" in kwargs.keys(): kwargsCopy["name"] = kwargs["name"] + str(j) self.HFEngine.outParaview(u, *args, **kwargsCopy) setattr(self.__class__, "outParaview" + fieldName, objFunc) def addOutParaviewTimeDomainFieldToClass(self, fieldName): def objFunc(self, mu:paramVal, *args, homogeneized : bool = False, **kwargs): if not hasattr(self.HFEngine, "outParaviewTimeDomain"): raise RROMPyException(("High fidelity engine cannot output to " "Paraview.")) uV = getattr(self.__class__, "get" + fieldName)(self, mu, homogeneized) omega = args.pop(0) if len(args) > 0 else np.real(mu) kwargsCopy = copy(kwargs) for j, u in enumerate(uV): if "name" in kwargs.keys(): kwargsCopy["name"] = kwargs["name"] + str(j) self.HFEngine.outParaviewTimeDomain(u, omega, *args, **kwargsCopy) setattr(self.__class__, "outParaviewTimeDomain" + fieldName, objFunc) class GenericApproximant: """ ABSTRACT ROM approximant computation for parametric problems. 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; - 'S': total number of samples current approximant relies upon. 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. trainedModel: Trained model evaluator. mu0: Default parameter. homogeneized: Whether to homogeneize Dirichlet BCs. approxParameters: Dictionary containing values for main parameters of approximant. Recognized keys are in parameterList{Soft,Critical}. parameterListSoft: Recognized keys of soft approximant parameters: - 'POD': whether to compute POD of snapshots. parameterListCritical: Recognized keys of critical approximant parameters: - 'S': total number of samples current approximant relies upon. verbosity: Verbosity level. POD: Whether to compute POD of snapshots. S: Number of solution snapshots over which current approximant is based upon. 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. """ __all__ += [ftype + dtype for ftype, dtype in iterprod( ["norm", "plot", "outParaview", "outParaviewTimeDomain"], ["HF", "RHS", "Approx", "Res", "Err"])] def __init__(self, HFEngine:HFEng, mu0 : paramVal = None, approxParameters : DictAny = {}, homogeneized : bool = False, verbosity : int = 10, timestamp : bool = True): self._preInit() self._mode = RROMPy_READY self.verbosity = verbosity self.timestamp = timestamp vbMng(self, "INIT", "Initializing engine of type {}.".format(self.name()), 10) self._HFEngine = HFEngine self.trainedModel = None self.lastSolvedHF = emptyParameterList() self.uHF = emptySampleList() self._addParametersToList(["POD"], [True], ["S"], [[1]]) if mu0 is None: if hasattr(self.HFEngine, "mu0"): self.mu0 = checkParameter(self.HFEngine.mu0) else: raise RROMPyException(("Center of approximation cannot be " "inferred from HF engine. Parameter " "required")) else: self.mu0 = checkParameter(mu0, self.HFEngine.npar) self.resetSamples() self.homogeneized = homogeneized self.approxParameters = approxParameters self._postInit() ### add norm{HF,RHS,Approx,Res,Err} methods """ Compute norm of * at arbitrary parameter. Args: mu: Target parameter. homogeneized(optional): Whether to remove Dirichlet BC. Defaults to False. Returns: Target norm of *. """ for objName in ["HF", "RHS", "Err"]: addNormFieldToClass(self, objName) def objFunc(self, mu:paramList, homogeneized : bool = False) -> Np1D: # uV = getattr(self.__class__, "getRes")(self, mu, homogeneized, # duality = False) uV = self.getRes(mu, homogeneized, duality = False) val = self.HFEngine.norm(uV, dual = True, duality = False) return val setattr(self.__class__, "normRes", objFunc) if not hasattr(self, "normApprox"): addNormFieldToClass(self, "Approx") ### add plot{HF,RHS,Approx,Res,Err} methods """ Do some nice plots of * at arbitrary parameter. Args: mu: Target parameter. name(optional): Name to be shown as title of the plots. Defaults to 'u'. what(optional): Which plots to do. If list, can contain 'ABS', 'PHASE', 'REAL', 'IMAG'. If str, same plus wildcard 'ALL'. Defaults to 'ALL'. save(optional): Where to save plot(s). Defaults to None, i.e. no saving. saveFormat(optional): Format for saved plot(s). Defaults to "eps". saveDPI(optional): DPI for saved plot(s). Defaults to 100. show(optional): Whether to show figure. Defaults to True. homogeneized(optional): Whether to remove Dirichlet BC. Defaults to False. figspecs(optional key args): Optional arguments for matplotlib figure creation. """ for objName in ["HF", "RHS", "Approx", "Res", "Err"]: addPlotFieldToClass(self, objName) ### add outParaview{HF,RHS,Approx,Res,Err} methods """ Output * to ParaView file. Args: mu: Target parameter. name(optional): Base name to be used for data output. filename(optional): Name of output file. time(optional): Timestamp. what(optional): Which plots to do. If list, can contain 'MESH', 'ABS', 'PHASE', 'REAL', 'IMAG'. If str, same plus wildcard 'ALL'. Defaults to 'ALL'. forceNewFile(optional): Whether to create new output file. filePW(optional): Fenics File entity (for time series). homogeneized(optional): Whether to remove Dirichlet BC. Defaults to False. """ for objName in ["HF", "RHS", "Approx", "Res", "Err"]: addOutParaviewFieldToClass(self, objName) ### add outParaviewTimeDomain{HF,RHS,Approx,Res,Err} methods """ Output * to ParaView file, converted to time domain. Args: mu: Target parameter. omega(optional): frequency. timeFinal(optional): final time of simulation. periodResolution(optional): number of time steps per period. name(optional): Base name to be used for data output. filename(optional): Name of output file. forceNewFile(optional): Whether to create new output file. homogeneized(optional): Whether to remove Dirichlet BC. Defaults to False. """ for objName in ["HF", "RHS", "Approx", "Res", "Err"]: addOutParaviewTimeDomainFieldToClass(self, objName) def _preInit(self): if not hasattr(self, "depth"): self.depth = 0 else: self.depth += 1 @property def parameterList(self): """Value of parameterListSoft + parameterListCritical.""" return self.parameterListSoft + self.parameterListCritical def _addParametersToList(self, whatSoft:strLst, defaultSoft:ListAny, whatCritical : strLst = [], defaultCritical : ListAny = [], toBeExcluded : strLst = []): if not hasattr(self, "parameterToBeExcluded"): self.parameterToBeExcluded = [] self.parameterToBeExcluded += toBeExcluded if not hasattr(self, "parameterListSoft"): self.parameterListSoft = [] if not hasattr(self, "parameterDefaultSoft"): self.parameterDefaultSoft = {} if not hasattr(self, "parameterListCritical"): self.parameterListCritical = [] if not hasattr(self, "parameterDefaultCritical"): self.parameterDefaultCritical = {} for j, what in enumerate(whatSoft): if what not in self.parameterToBeExcluded: self.parameterListSoft += [what] self.parameterDefaultSoft[what] = defaultSoft[j] for j, what in enumerate(whatCritical): if what not in self.parameterToBeExcluded: self.parameterListCritical += [what] self.parameterDefaultCritical[what] = defaultCritical[j] def _postInit(self): if self.depth == 0: vbMng(self, "DEL", "Done initializing.", 10) del self.depth else: self.depth -= 1 def name(self) -> str: return self.__class__.__name__ def __str__(self) -> str: return self.name() def __repr__(self) -> str: return self.__str__() + " at " + hex(id(self)) 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 = SamplingEngineStandardPOD else: SamplingEngine = SamplingEngineStandard self.samplingEngine = SamplingEngine(self.HFEngine, verbosity = self.verbosity, allowRepeatedSamples = True) @property def HFEngine(self): """Value of HFEngine.""" return self._HFEngine @HFEngine.setter def HFEngine(self, HFEngine): raise RROMPyException("Cannot change HFEngine.") @property def mu0(self): """Value of mu0.""" return self._mu0 @mu0.setter def mu0(self, mu0): mu0 = checkParameter(mu0) if not hasattr(self, "_mu0") or mu0 != self.mu0: self.resetSamples() self._mu0 = mu0 @property def npar(self): """Number of parameters.""" return self.mu0.shape[1] @property def approxParameters(self): """Value of approximant parameters.""" return self._approxParameters @approxParameters.setter def approxParameters(self, approxParams): if not hasattr(self, "approxParameters"): self._approxParameters = {} approxParameters = purgeDict(approxParams, self.parameterList, dictname = self.name() + ".approxParameters", baselevel = 1) keyList = list(approxParameters.keys()) for key in self.parameterListCritical: if key in keyList: setattr(self, "_" + key, self.parameterDefaultCritical[key]) for key in self.parameterListSoft: if key in keyList: setattr(self, "_" + key, self.parameterDefaultSoft[key]) fragile = False for key in self.parameterListCritical: if key in keyList: val = approxParameters[key] else: val = getattr(self, "_" + key, None) if val is None: val = self.parameterDefaultCritical[key] getattr(self.__class__, key, None).fset(self, val) fragile = fragile or val is None for key in self.parameterListSoft: if key in keyList: val = approxParameters[key] else: val = getattr(self, "_" + key, None) if val is None: val = self.parameterDefaultSoft[key] getattr(self.__class__, key, None).fset(self, val) if fragile: self._mode = RROMPy_FRAGILE @property def POD(self): """Value of POD.""" return self._POD @POD.setter def POD(self, POD): if hasattr(self, "_POD"): PODold = self.POD else: PODold = -1 self._POD = POD self._approxParameters["POD"] = self.POD if PODold != self.POD: self.samplingEngine = None self.resetSamples() @property def S(self): """Value of S.""" return self._S @S.setter def S(self, S): if not hasattr(S, "__len__"): S = [S] if any([s <= 0 for s in S]): raise RROMPyException("S must be positive.") if hasattr(self, "_S") and self._S is not None: Sold = tuple(self.S) else: Sold = -1 self._S = S self._approxParameters["S"] = self.S if Sold != tuple(self.S): self.resetSamples() @property def homogeneized(self): """Value of homogeneized.""" return self._homogeneized @homogeneized.setter def homogeneized(self, homogeneized): if not hasattr(self, "_homogeneized"): self._homogeneized = None if homogeneized != self.homogeneized: self._homogeneized = homogeneized self.resetSamples() @property def trainedModel(self): """Value of trainedModel.""" return self._trainedModel @trainedModel.setter def trainedModel(self, trainedModel): self._trainedModel = trainedModel if self._trainedModel is not None: self._trainedModel.lastSolvedApproxReduced = emptyParameterList() self._trainedModel.lastSolvedApprox = emptyParameterList() self.lastSolvedApproxReduced = emptyParameterList() self.lastSolvedApprox = emptyParameterList() self.uApproxReduced = emptySampleList() self.uApprox = emptySampleList() def resetSamples(self): if hasattr(self, "samplingEngine") and self.samplingEngine is not None: self.samplingEngine.resetHistory() else: self.setupSampling() self._mode = RROMPy_READY def plotSamples(self, warping : List[callable] = None, name : str = "u", save : str = None, what : strLst = 'all', saveFormat : str = "eps", saveDPI : int = 100, show : bool = True, plotArgs : dict = {}, **figspecs): """ Do some nice plots of the samples. Args: warping(optional): Domain warping functions. name(optional): Name to be shown as title of the plots. Defaults to 'u'. what(optional): Which plots to do. If list, can contain 'ABS', 'PHASE', 'REAL', 'IMAG'. If str, same plus wildcard 'ALL'. Defaults to 'ALL'. save(optional): Where to save plot(s). Defaults to None, i.e. no saving. saveFormat(optional): Format for saved plot(s). Defaults to "eps". saveDPI(optional): DPI for saved plot(s). Defaults to 100. show(optional): Whether to show figure. Defaults to True. plotArgs(optional): Optional arguments for fen/pyplot. figspecs(optional key args): Optional arguments for matplotlib figure creation. """ RROMPyAssert(self._mode, message = "Cannot plot samples.") self.samplingEngine.plotSamples(warping, name, save, what, saveFormat, saveDPI, show, plotArgs, **figspecs) def outParaviewSamples(self, name : str = "u", filename : str = "out", times : Np1D = None, what : strLst = 'all', forceNewFile : bool = True, folders : bool = False, filePW = None): """ Output samples to ParaView file. Args: name(optional): Base name to be used for data output. filename(optional): Name of output file. times(optional): Timestamps. what(optional): Which plots to do. If list, can contain 'MESH', 'ABS', 'PHASE', 'REAL', 'IMAG'. If str, same plus wildcard 'ALL'. Defaults to 'ALL'. forceNewFile(optional): Whether to create new output file. folders(optional): Whether to split output in folders. filePW(optional): Fenics File entity (for time series). """ RROMPyAssert(self._mode, message = "Cannot output samples.") self.samplingEngine.outParaviewSamples(name = name, filename = filename, times = times, what = what, forceNewFile = forceNewFile, folders = folders, filePW = filePW) def outParaviewTimeDomainSamples(self, omegas : Np1D = None, timeFinal : Np1D = None, periodResolution : int = 20, name : str = "u", filename : str = "out", forceNewFile : bool = True, folders : bool = False): """ Output samples to ParaView file, converted to time domain. Args: omegas(optional): frequencies. timeFinal(optional): final time of simulation. periodResolution(optional): number of time steps per period. name(optional): Base name to be used for data output. filename(optional): Name of output file. forceNewFile(optional): Whether to create new output file. folders(optional): Whether to split output in folders. """ RROMPyAssert(self._mode, message = "Cannot output samples.") self.samplingEngine.outParaviewTimeDomainSamples(omegas = omegas, timeFinal = timeFinal, periodResolution = periodResolution, name = name, filename = filename, forceNewFile = forceNewFile, folders = folders) def setSamples(self, samplingEngine): """Copy samplingEngine and samples.""" vbMng(self, "INIT", "Transfering samples.", 10) self.samplingEngine = copy(samplingEngine) vbMng(self, "DEL", "Done transfering samples.", 10) def setTrainedModel(self, model): """Deepcopy approximation from trained model.""" if hasattr(model, "storeTrainedModel"): verb = model.verbosity model.verbosity = 0 fileOut = model.storeTrainedModel() model.verbosity = verb else: try: fileOut = getNewFilename("trained_model", "pkl") pickleDump(model.data.__dict__, fileOut) except: raise RROMPyException(("Failed to store model data. Parameter " "model must have either " "storeTrainedModel or " "data.__dict__ properties.")) self.loadTrainedModel(fileOut) osrm(fileOut) @abstractmethod def setupApprox(self): """ Setup approximant. (ABSTRACT) Any specialization should include something like if self.checkComputedApprox(): return RROMPyAssert(self._mode, message = "Cannot setup approximant.") ... self.trainedModel = ... self.trainedModel.data = ... self.trainedModel.data.approxParameters = copy( self.approxParameters) """ pass def checkComputedApprox(self) -> bool: """ Check if setup of new approximant is not needed. Returns: True if new setup is not needed. False otherwise. """ return self._mode == RROMPy_FRAGILE or (self.trainedModel is not None and self.trainedModel.data.approxParameters == self.approxParameters) def _pruneBeforeEval(self, mu:paramList, field:str, append:bool, prune:bool) -> Tuple[paramList, Np1D, Np1D, bool]: mu = checkParameterList(mu, self.npar)[0] idx = np.empty(len(mu), dtype = np.int) if prune: jExtra = np.zeros(len(mu), dtype = bool) muExtra = emptyParameterList() lastSolvedMus = getattr(self, "lastSolved" + field) if (len(mu) > 0 and len(mu) == len(lastSolvedMus) and mu == lastSolvedMus): idx = np.arange(len(mu), dtype = np.int) return muExtra, jExtra, idx, True muKeep = copy(muExtra) for j in range(len(mu)): jPos = lastSolvedMus.find(mu[j]) if jPos is not None: idx[j] = jPos muKeep.append(mu[j]) else: jExtra[j] = True muExtra.append(mu[j]) if len(muKeep) > 0 and not append: lastSolvedu = getattr(self, "u" + field) idx[~jExtra] = getattr(self.__class__, "set" + field)(self, muKeep, lastSolvedu[idx[~jExtra]], append) append = True else: jExtra = np.ones(len(mu), dtype = bool) muExtra = mu return muExtra, jExtra, idx, append def _setObject(self, mu:paramList, field:str, object:sampList, append:bool) -> List[int]: newMus = checkParameterList(mu, self.npar)[0] newObj = sampleList(object) if append: getattr(self, "lastSolved" + field).append(newMus) getattr(self, "u" + field).append(newObj) Ltot = len(getattr(self, "u" + field)) return list(range(Ltot - len(newObj), Ltot)) setattr(self, "lastSolved" + field, copy(newMus)) setattr(self, "u" + field, copy(newObj)) return list(range(len(getattr(self, "u" + field)))) def setHF(self, muHF:paramList, uHF:sampleList, append : bool = False) -> List[int]: """Assign high fidelity solution.""" return self._setObject(muHF, "HF", uHF, append) def evalHF(self, mu:paramList, append : bool = False, prune : bool = True) -> List[int]: """ Find high fidelity solution with original parameters and arbitrary parameter. Args: mu: Target parameter. append(optional): Whether to append new HF solutions to old ones. prune(optional): Whether to remove duplicates of already appearing HF solutions. """ muExtra, jExtra, idx, append = self._pruneBeforeEval(mu, "HF", append, prune) if len(muExtra) > 0: vbMng(self, "INIT", "Solving HF model for mu = {}.".format(mu), 15) newuHFs = self.HFEngine.solve(muExtra, homogeneized = self.homogeneized) vbMng(self, "DEL", "Done solving HF model.", 15) idx[jExtra] = self.setHF(muExtra, newuHFs, append) return list(idx) def setApproxReduced(self, muApproxR:paramList, uApproxR:sampleList, append : bool = False) -> List[int]: """Assign high fidelity solution.""" return self._setObject(muApproxR, "ApproxReduced", uApproxR, append) def evalApproxReduced(self, mu:paramList, append : bool = False, prune : bool = True) -> List[int]: """ Evaluate reduced representation of approximant at arbitrary parameter. Args: mu: Target parameter. append(optional): Whether to append new HF solutions to old ones. prune(optional): Whether to remove duplicates of already appearing HF solutions. """ self.setupApprox() muExtra, jExtra, idx, append = self._pruneBeforeEval(mu, "ApproxReduced", append, prune) if len(muExtra) > 0: newuApproxs = self.trainedModel.getApproxReduced(muExtra) idx[jExtra] = self.setApproxReduced(muExtra, newuApproxs, append) return list(idx) def setApprox(self, muApprox:paramList, uApprox:sampleList, append : bool = False) -> List[int]: """Assign high fidelity solution.""" return self._setObject(muApprox, "Approx", uApprox, append) def evalApprox(self, mu:paramList, append : bool = False, prune : bool = True) -> List[int]: """ Evaluate approximant at arbitrary parameter. Args: mu: Target parameter. append(optional): Whether to append new HF solutions to old ones. prune(optional): Whether to remove duplicates of already appearing HF solutions. """ self.setupApprox() muExtra, jExtra, idx, append = self._pruneBeforeEval(mu, "Approx", append, prune) if len(muExtra) > 0: newuApproxs = self.trainedModel.getApprox(muExtra) idx[jExtra] = self.setApprox(muExtra, newuApproxs, append) return list(idx) def getHF(self, mu:paramList, homogeneized : bool = False, append : bool = False, prune : bool = True) -> sampList: """ Get HF solution at arbitrary parameter. Args: mu: Target parameter. homogeneized(optional): Whether to remove Dirichlet BC. Defaults to False. Returns: HFsolution. """ mu = checkParameterList(mu, self.npar)[0] idx = self.evalHF(mu, append = append, prune = prune) uHFs = self.uHF(idx) if self.homogeneized and not homogeneized: for j, m in enumerate(mu): uHFs[j] += self.HFEngine.liftDirichletData(m) if not self.homogeneized and homogeneized: for j, m in enumerate(mu): uHFs[j] -= self.HFEngine.liftDirichletData(m) return uHFs def getRHS(self, mu:paramList, homogeneized : bool = False, duality : bool = True) -> sampList: """ Get linear system RHS at arbitrary parameter. Args: mu: Target parameter. homogeneized(optional): Whether to remove Dirichlet BC. Defaults to False. Returns: Linear system RHS. """ return self.HFEngine.residual(None, mu, homogeneized = homogeneized, duality = duality) def getApproxReduced(self, mu:paramList, append : bool = False, prune : bool = True) -> sampList: """ Get approximant at arbitrary parameter. Args: mu: Target parameter. Returns: Reduced approximant. """ mu = checkParameterList(mu, self.npar)[0] idx = self.evalApproxReduced(mu, append = append, prune = prune) uApproxRs = self.uApproxReduced(idx) return uApproxRs def getApprox(self, mu:paramList, homogeneized : bool = False, append : bool = False, prune : bool = True) -> sampList: """ Get approximant at arbitrary parameter. Args: mu: Target parameter. homogeneized(optional): Whether to remove Dirichlet BC. Defaults to False. Returns: Approximant. """ mu = checkParameterList(mu, self.npar)[0] idx = self.evalApprox(mu, append = append, prune = prune) uApproxs = self.uApprox(idx) if self.homogeneized and not homogeneized: for j, m in enumerate(mu): uApproxs[j] += self.HFEngine.liftDirichletData(m) if not self.homogeneized and homogeneized: for j, m in enumerate(mu): uApproxs[j] -= self.HFEngine.liftDirichletData(m) return uApproxs def getRes(self, mu:paramList, homogeneized : bool = False, duality : bool = True) -> sampList: """ Get residual at arbitrary parameter. Args: mu: Target parameter. homogeneized(optional): Whether to remove Dirichlet BC. Defaults to False. Returns: Approximant residual. """ return self.HFEngine.residual(self.getApprox(mu, homogeneized), mu, homogeneized = homogeneized, duality = duality) def getErr(self, mu:paramList, homogeneized : bool = False, append : bool = False, prune : bool = True) -> sampList: """ Get error at arbitrary parameter. Args: mu: Target parameter. homogeneized(optional): Whether to remove Dirichlet BC. Defaults to False. Returns: Approximant error. """ return (self.getApprox(mu, homogeneized, append = append, prune =prune) - self.getHF(mu, homogeneized, append = append, prune = prune)) def getPoles(self, *args, **kwargs) -> Np1D: """ Obtain approximant poles. Returns: Numpy complex vector of poles. """ self.setupApprox() vbMng(self, "INIT", "Computing poles of model.", 20) poles = self.trainedModel.getPoles(*args, **kwargs) vbMng(self, "DEL", "Done computing poles.", 20) return poles def storeTrainedModel(self, filenameBase : str = "trained_model", forceNewFile : bool = True) -> str: """Store trained reduced model to file.""" self.setupApprox() vbMng(self, "INIT", "Storing trained model to file.", 20) if forceNewFile: filename = getNewFilename(filenameBase, "pkl") else: filename = "{}.pkl".format(filenameBase) pickleDump(self.trainedModel.data.__dict__, filename) vbMng(self, "DEL", "Done storing trained model.", 20) return filename def loadTrainedModel(self, filename:str): """Load trained reduced model from file.""" vbMng(self, "INIT", "Loading pre-trained model from file.", 20) datadict = pickleLoad(filename) name = datadict.pop("name") if name == "TrainedModelRational": from rrompy.reduction_methods.trained_model import \ TrainedModelRational as tModel - elif name == "TrainedModelRB": + elif name == "TrainedModelReducedBasis": from rrompy.reduction_methods.trained_model import \ - TrainedModelRB as tModel + TrainedModelReducedBasis as tModel else: raise RROMPyException(("Trained model name not recognized. " "Loading failed.")) self.mu0 = datadict.pop("mu0") from rrompy.reduction_methods.trained_model import TrainedModelData trainedModel = tModel() trainedModel.verbosity = self.verbosity trainedModel.timestamp = self.timestamp self.scaleFactor = datadict.pop("scaleFactor") data = TrainedModelData(name, self.mu0, datadict.pop("projMat"), self.scaleFactor, datadict.pop("rescalingExp")) if "mus" in datadict: data.mus = datadict.pop("mus") approxParameters = datadict.pop("approxParameters") data.approxParameters = copy(approxParameters) if "sampler" in approxParameters: self._approxParameters["sampler"] = approxParameters.pop("sampler") self.approxParameters = copy(approxParameters) if "mus" in data.__dict__: self.mus = copy(data.mus) for key in datadict: setattr(data, key, datadict[key]) trainedModel.data = data self.trainedModel = trainedModel self._mode = RROMPy_FRAGILE vbMng(self, "DEL", "Done loading pre-trained model.", 20) diff --git a/rrompy/reduction_methods/base/pade_utils.py b/rrompy/reduction_methods/base/rational_interpolant_utils.py similarity index 100% rename from rrompy/reduction_methods/base/pade_utils.py rename to rrompy/reduction_methods/base/rational_interpolant_utils.py diff --git a/rrompy/reduction_methods/base/rb_utils.py b/rrompy/reduction_methods/base/reduced_basis_utils.py similarity index 100% rename from rrompy/reduction_methods/base/rb_utils.py rename to rrompy/reduction_methods/base/reduced_basis_utils.py diff --git a/rrompy/reduction_methods/distributed/__init__.py b/rrompy/reduction_methods/distributed/__init__.py index 1d7c1fc..b136774 100644 --- a/rrompy/reduction_methods/distributed/__init__.py +++ b/rrompy/reduction_methods/distributed/__init__.py @@ -1,29 +1,29 @@ # 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_distributed_approximant import GenericDistributedApproximant from .rational_interpolant import RationalInterpolant -from .rb_distributed import RBDistributed +from .reduced_basis import ReducedBasis __all__ = [ 'GenericDistributedApproximant', 'RationalInterpolant', - 'RBDistributed' + 'ReducedBasis' ] diff --git a/rrompy/reduction_methods/distributed/rb_distributed.py b/rrompy/reduction_methods/distributed/reduced_basis.py similarity index 96% rename from rrompy/reduction_methods/distributed/rb_distributed.py rename to rrompy/reduction_methods/distributed/reduced_basis.py index 3f142ed..d63bbf4 100644 --- a/rrompy/reduction_methods/distributed/rb_distributed.py +++ b/rrompy/reduction_methods/distributed/reduced_basis.py @@ -1,207 +1,209 @@ # 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 copy import deepcopy as copy import numpy as np from .generic_distributed_approximant import GenericDistributedApproximant -from rrompy.reduction_methods.trained_model import TrainedModelRB as tModel +from rrompy.reduction_methods.trained_model import \ + TrainedModelReducedBasis as tModel from rrompy.reduction_methods.trained_model import TrainedModelData -from rrompy.reduction_methods.base.rb_utils import projectAffineDecomposition +from rrompy.reduction_methods.base.reduced_basis_utils import \ + projectAffineDecomposition from rrompy.utilities.base.types import (Np1D, Np2D, List, ListAny, Tuple, DictAny, HFEng, paramVal, sampList) from rrompy.utilities.base import verbosityManager as vbMng from rrompy.utilities.exception_manager import (RROMPyWarning, RROMPyException, RROMPyAssert) -__all__ = ['RBDistributed'] +__all__ = ['ReducedBasis'] -class RBDistributed(GenericDistributedApproximant): +class ReducedBasis(GenericDistributedApproximant): """ ROM RB approximant computation for parametric problems. 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; - 'S': total number of samples current approximant relies upon; - 'sampler': sample point generator; - 'R': rank for Galerkin projection; defaults to prod(S); - 'PODTolerance': tolerance for snapshots POD; defaults to -1. Defaults to empty dict. homogeneized(optional): Whether to homogeneize Dirichlet BCs. Defaults to False. verbosity(optional): Verbosity level. Defaults to 10. Attributes: HFEngine: HF problem solver. mu0: Default parameter. mus: Array of snapshot parameters. homogeneized: Whether to homogeneize Dirichlet BCs. approxRadius: Dummy radius of approximant (i.e. distance from mu0 to farthest sample point). 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. - 'R': rank for Galerkin projection; - 'PODTolerance': tolerance for snapshots POD. parameterListCritical: Recognized keys of critical approximant parameters: - 'S': total number of samples current approximant relies upon; - 'sampler': sample point generator; POD: Whether to compute POD of snapshots. S: Number of solution snapshots over which current approximant is based upon. sampler: Sample point generator. R: Rank for Galerkin projection. muBounds: list of bounds for 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. As: List of sparse matrices (in CSC format) representing coefficients of linear system matrix. bs: List of numpy vectors representing coefficients of linear system RHS. ARBs: List of sparse matrices (in CSC format) representing coefficients of compressed linear system matrix. bRBs: List of numpy vectors representing coefficients of compressed linear system RHS. """ def __init__(self, HFEngine:HFEng, mu0 : paramVal = None, approxParameters : DictAny = {}, homogeneized : bool = False, verbosity : int = 10, timestamp : bool = True): self._preInit() self._addParametersToList(["R", "PODTolerance"], [1, -1]) super().__init__(HFEngine = HFEngine, mu0 = mu0, approxParameters = approxParameters, homogeneized = homogeneized, verbosity = verbosity, timestamp = timestamp) vbMng(self, "INIT", "Computing affine blocks of system.", 10) vbMng(self, "DEL", "Done computing affine blocks.", 10) self._postInit() @property def S(self): """Value of S.""" return self._S @S.setter def S(self, S): GenericDistributedApproximant.S.fset(self, S) if not hasattr(self, "_R"): self._R = np.prod(self.S) @property def R(self): """Value of R. Its assignment may change S.""" return self._R @R.setter def R(self, R): if R < 0: raise RROMPyException("R must be non-negative.") self._R = R self._approxParameters["R"] = self.R if hasattr(self, "_S") and np.prod(self.S) < self.R: RROMPyWarning("Prescribed S is too small. Decreasing R.") self.R = np.prod(self.S) @property def PODTolerance(self): """Value of PODTolerance.""" return self._PODTolerance @PODTolerance.setter def PODTolerance(self, PODTolerance): self._PODTolerance = PODTolerance self._approxParameters["PODTolerance"] = self.PODTolerance def setupApprox(self): """Compute RB projection matrix.""" if self.checkComputedApprox(): return RROMPyAssert(self._mode, message = "Cannot setup approximant.") vbMng(self, "INIT", "Setting up {}.". format(self.name()), 5) self.computeScaleFactor() self.computeSnapshots() vbMng(self, "INIT", "Computing projection matrix.", 7) if self.POD: U, s, _ = np.linalg.svd(self.samplingEngine.RPOD) s = s ** 2. else: Gramian = self.HFEngine.innerProduct(self.samplingEngine.samples, self.samplingEngine.samples) U, s, _ = np.linalg.svd(Gramian) nsamples = self.samplingEngine.nsamples snorm = np.cumsum(s[::-1]) / np.sum(s) nPODTrunc = min(nsamples - np.argmax(snorm > self.PODTolerance), self.R) pMat = self.samplingEngine.samples.dot(U[:, : nPODTrunc]) vbMng(self, "MAIN", ("Assembling {}x{} projection matrix from {} " "samples.").format(*(pMat.shape), nsamples), 5) if self.trainedModel is None: self.trainedModel = tModel() self.trainedModel.verbosity = self.verbosity self.trainedModel.timestamp = self.timestamp data = TrainedModelData(self.trainedModel.name(), self.mu0, pMat, self.scaleFactor, self.HFEngine.rescalingExp) data.mus = copy(self.mus) self.trainedModel.data = data else: self.trainedModel = self.trainedModel self.trainedModel.data.projMat = copy(pMat) vbMng(self, "INIT", "Computing affine blocks of system.", 10) self.As = self.HFEngine.affineLinearSystemA(self.mu0, self.scaleFactor) self.bs = self.HFEngine.affineLinearSystemb(self.mu0, self.scaleFactor, self.homogeneized) vbMng(self, "DEL", "Done computing affine blocks.", 10) ARBs, bRBs = self.assembleReducedSystem(pMat) self.trainedModel.data.ARBs = ARBs self.trainedModel.data.bRBs = bRBs vbMng(self, "DEL", "Done computing projection matrix.", 7) self.trainedModel.data.approxParameters = copy(self.approxParameters) vbMng(self, "DEL", "Done setting up approximant.", 5) def assembleReducedSystem(self, pMat : sampList = None, pMatOld : sampList = None)\ -> Tuple[List[Np2D], List[Np1D]]: """Build affine blocks of RB linear system through projections.""" if pMat is None: self.setupApprox() ARBs = self.trainedModel.data.ARBs bRBs = self.trainedModel.data.bRBs else: vbMng(self, "INIT", "Projecting affine terms of HF model.", 10) ARBsOld = None if pMatOld is None else self.trainedModel.data.ARBs bRBsOld = None if pMatOld is None else self.trainedModel.data.bRBs ARBs, bRBs = projectAffineDecomposition(self.As, self.bs, pMat, ARBsOld, bRBsOld, pMatOld) vbMng(self, "DEL", "Done projecting affine terms.", 10) return ARBs, bRBs diff --git a/rrompy/reduction_methods/distributed_greedy/__init__.py b/rrompy/reduction_methods/distributed_greedy/__init__.py index ba0008b..ccdf783 100644 --- a/rrompy/reduction_methods/distributed_greedy/__init__.py +++ b/rrompy/reduction_methods/distributed_greedy/__init__.py @@ -1,30 +1,30 @@ # 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_distributed_greedy_approximant import \ GenericDistributedGreedyApproximant from .rational_interpolant_greedy import RationalInterpolantGreedy -from .rb_distributed_greedy import RBDistributedGreedy +from .reduced_basis_greedy import ReducedBasisGreedy __all__ = [ 'GenericDistributedGreedyApproximant', 'RationalInterpolantGreedy', - 'RBDistributedGreedy' + 'ReducedBasisGreedy' ] diff --git a/rrompy/reduction_methods/distributed_greedy/rb_distributed_greedy.py b/rrompy/reduction_methods/distributed_greedy/reduced_basis_greedy.py similarity index 97% rename from rrompy/reduction_methods/distributed_greedy/rb_distributed_greedy.py rename to rrompy/reduction_methods/distributed_greedy/reduced_basis_greedy.py index 1bb837c..9fc0247 100644 --- a/rrompy/reduction_methods/distributed_greedy/rb_distributed_greedy.py +++ b/rrompy/reduction_methods/distributed_greedy/reduced_basis_greedy.py @@ -1,238 +1,239 @@ # 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_distributed_greedy_approximant import \ GenericDistributedGreedyApproximant -from rrompy.reduction_methods.distributed import RBDistributed -from rrompy.reduction_methods.trained_model import TrainedModelRB as tModel +from rrompy.reduction_methods.distributed import ReducedBasis +from rrompy.reduction_methods.trained_model import \ + TrainedModelReducedBasis as tModel from rrompy.reduction_methods.trained_model import TrainedModelData from rrompy.utilities.base.types import Np1D, DictAny, HFEng, paramVal from rrompy.utilities.base import verbosityManager as vbMng from rrompy.utilities.poly_fitting.polynomial import ( hashIdxToDerivative as hashI) from rrompy.utilities.exception_manager import RROMPyWarning, RROMPyAssert from rrompy.parameter import checkParameterList -__all__ = ['RBDistributedGreedy'] +__all__ = ['ReducedBasisGreedy'] -class RBDistributedGreedy(GenericDistributedGreedyApproximant, RBDistributed): +class ReducedBasisGreedy(GenericDistributedGreedyApproximant, ReducedBasis): """ ROM greedy RB approximant computation for parametric problems. 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; - 'S': number of starting training points; - 'sampler': sample point generator; - 'greedyTol': uniform error tolerance for greedy algorithm; defaults to 1e-2; - 'interactive': whether to interactively terminate greedy algorithm; defaults to False; - 'maxIter': maximum number of greedy steps; defaults to 1e2; - 'refinementRatio': ratio of test points to be exhausted before test set refinement; defaults to 0.2; - 'nTestPoints': number of test points; defaults to 5e2; - 'trainSetGenerator': training sample points generator; defaults to Chebyshev sampler within muBounds. Defaults to empty dict. homogeneized(optional): Whether to homogeneize Dirichlet BCs. Defaults to False. verbosity(optional): Verbosity level. Defaults to 10. Attributes: HFEngine: HF problem solver. mu0: Default parameter. mus: Array of snapshot parameters. homogeneized: Whether to homogeneize Dirichlet BCs. approxParameters: Dictionary containing values for main parameters of approximant. Recognized keys are in parameterList. parameterListSoft: Recognized keys of soft approximant parameters: - 'POD': whether to compute POD of snapshots. - 'greedyTol': uniform error tolerance for greedy algorithm; - 'interactive': whether to interactively terminate greedy algorithm; - 'maxIter': maximum number of greedy steps; - 'refinementRatio': ratio of test points to be exhausted before test set refinement; - 'nTestPoints': number of test points; - 'trainSetGenerator': training sample points generator. parameterListCritical: Recognized keys of critical approximant parameters: - 'S': total number of samples current approximant relies upon; - 'sampler': sample point generator. POD: whether to compute POD of snapshots. S: number of test points. sampler: Sample point generator. greedyTol: uniform error tolerance for greedy algorithm. interactive: whether to interactively terminate greedy algorithm. maxIter: maximum number of greedy steps. refinementRatio: ratio of training points to be exhausted before training set refinement. 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. As: List of sparse matrices (in CSC format) representing coefficients of linear system matrix. bs: List of numpy vectors representing coefficients of linear system RHS. ARBs: List of sparse matrices (in CSC format) representing coefficients of compressed linear system matrix. bRBs: List of numpy vectors representing coefficients of compressed linear system RHS. """ def __init__(self, HFEngine:HFEng, mu0 : paramVal = None, approxParameters : DictAny = {}, homogeneized : bool = False, verbosity : int = 10, timestamp : bool = True): self._preInit() super().__init__(HFEngine = HFEngine, mu0 = mu0, approxParameters = approxParameters, homogeneized = homogeneized, verbosity = verbosity, timestamp = timestamp) self._postInit() @property def R(self): """Value of R.""" self._R = np.prod(self._S) return self._R @R.setter def R(self, R): RROMPyWarning(("R is used just to simplify inheritance, and its value " "cannot be changed from that of prod(S).")) @property def PODTolerance(self): """Value of PODTolerance.""" self._PODTolerance = -1 return self._PODTolerance @PODTolerance.setter def PODTolerance(self, PODTolerance): RROMPyWarning(("PODTolerance is used just to simplify inheritance, " "and its value cannot be changed from -1.")) def errorEstimator(self, mus:Np1D) -> Np1D: """ Standard residual-based error estimator. Unreliable for unstable problems (inf-sup constant is missing). """ self.setupApprox() self.assembleReducedResidualBlocks(full = True) nmus = len(mus) nAs = self.trainedModel.data.resAA.shape[1] nbs = self.trainedModel.data.resbb.shape[0] radiusA = np.empty((len(self.mus), nAs, nmus), dtype = np.complex) radiusb = np.empty((nbs, nmus), dtype = np.complex) verb = self.trainedModel.verbosity self.trainedModel.verbosity = 0 if verb >= 5: mustr = mus if nmus > 2: mustr = "[{} ..({}).. {}]".format(mus[0], nmus - 2, mus[-1]) vbMng(self, "INIT", "Computing RB solution at mu = {}.".format(mustr), 0) parmus = checkParameterList(mus, self.npar)[0] uApproxRs = self.getApproxReduced(parmus) mu0Eff = self.mu0(0, 0) ** self.HFEngine.rescalingExp[0] for j, muPL in enumerate(parmus): uApproxR = uApproxRs[j] for i in range(max(nAs, nbs)): thetai = ((muPL[0] ** self.HFEngine.rescalingExp[0] - mu0Eff) / self.scaleFactor[0]) ** i if i < nAs: radiusA[:, i, j] = thetai * uApproxR if i < nbs: radiusb[i, j] = thetai vbMng(self, "DEL", "Done computing RB solution.", 5) self.trainedModel.verbosity = verb # 'ij,jk,ik->k', resbb, radiusb, radiusb.conj() ff = np.sum(self.trainedModel.data.resbb.dot(radiusb) * radiusb.conj(), axis = 0) # 'ijk,jkl,il->l', resAb, radiusA, radiusb.conj() Lf = np.sum(np.tensordot(self.trainedModel.data.resAb, radiusA, 2) * radiusb.conj(), axis = 0) # 'ijkl,klt,ijt->t', resAA, radiusA, radiusA.conj() LL = np.sum(np.tensordot(self.trainedModel.data.resAA, radiusA, 2) * radiusA.conj(), axis = (0, 1)) return np.abs((LL - 2. * np.real(Lf) + ff) / ff) ** .5 def setupApprox(self, plotEst : bool = False): """Compute RB projection matrix.""" if self.checkComputedApprox(): return RROMPyAssert(self._mode, message = "Cannot setup approximant.") vbMng(self, "INIT", "Setting up {}.". format(self.name()), 5) self.computeScaleFactor() self.greedy(plotEst) vbMng(self, "INIT", "Computing projection matrix.", 7) pMat = self.samplingEngine.samples vbMng(self, "INIT", "Computing affine blocks of system.", 10) self.As = self.HFEngine.affineLinearSystemA(self.mu0, self.scaleFactor) self.bs = self.HFEngine.affineLinearSystemb(self.mu0, self.scaleFactor, self.homogeneized) vbMng(self, "DEL", "Done computing affine blocks.", 10) if self.trainedModel is None: self.trainedModel = tModel() self.trainedModel.verbosity = self.verbosity self.trainedModel.timestamp = self.timestamp data = TrainedModelData(self.trainedModel.name(), self.mu0, pMat, self.scaleFactor, self.HFEngine.rescalingExp) ARBs, bRBs = self.assembleReducedSystem(pMat) self.trainedModel.data = data else: self.trainedModel = self.trainedModel pMatOld = self.trainedModel.data.projMat Sold = pMatOld.shape[1] idxNew = list(range(Sold, pMat.shape[1])) ARBs, bRBs = self.assembleReducedSystem(pMat(idxNew), pMatOld) self.trainedModel.data.projMat = copy(pMat) self.trainedModel.data.mus = copy(self.mus) self.trainedModel.data.ARBs = ARBs self.trainedModel.data.bRBs = bRBs vbMng(self, "DEL", "Done computing projection matrix.", 7) self.trainedModel.data.approxParameters = copy(self.approxParameters) vbMng(self, "DEL", "Done setting up approximant.", 5) def assembleReducedResidualBlocks(self, full : bool = False): """Build affine blocks of RB linear system through projections.""" self.assembleReducedResidualBlocksbb(self.bs) if full: pMat = self.trainedModel.data.projMat self.assembleReducedResidualBlocksAb(self.As, self.bs, pMat) self.assembleReducedResidualBlocksAA(self.As, pMat) diff --git a/rrompy/reduction_methods/pivoting/__init__.py b/rrompy/reduction_methods/pivoting/__init__.py index 9830125..6237733 100644 --- a/rrompy/reduction_methods/pivoting/__init__.py +++ b/rrompy/reduction_methods/pivoting/__init__.py @@ -1,29 +1,29 @@ # 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_pivoted_approximant import GenericPivotedApproximant from .rational_interpolant_pivoted import RationalInterpolantPivoted -from .rb_distributed_pivoted import RBDistributedPivoted +from .reduced_basis_pivoted import ReducedBasisPivoted __all__ = [ 'GenericPivotedApproximant', 'RationalInterpolantPivoted', - 'RBDistributedPivoted' + 'ReducedBasisPivoted' ] diff --git a/rrompy/reduction_methods/pivoting/rb_distributed_pivoted.py b/rrompy/reduction_methods/pivoting/reduced_basis_pivoted.py similarity index 97% rename from rrompy/reduction_methods/pivoting/rb_distributed_pivoted.py rename to rrompy/reduction_methods/pivoting/reduced_basis_pivoted.py index 4fabd41..1200f83 100644 --- a/rrompy/reduction_methods/pivoting/rb_distributed_pivoted.py +++ b/rrompy/reduction_methods/pivoting/reduced_basis_pivoted.py @@ -1,284 +1,285 @@ # 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 copy import deepcopy as copy import numpy as np from .generic_pivoted_approximant import GenericPivotedApproximant from rrompy.hfengines.base import MarginalProxyEngine from rrompy.reduction_methods.trained_model import (TrainedModelPivotedData, - TrainedModelPivotedRB as tModel) -from rrompy.reduction_methods.base.rb_utils import projectAffineDecomposition + TrainedModelPivotedReducedBasis as tModel) +from rrompy.reduction_methods.base.reduced_basis_utils import \ + projectAffineDecomposition from rrompy.utilities.base.types import (Np1D, Np2D, List, ListAny, Tuple, DictAny, HFEng, paramVal, sampList) from rrompy.utilities.base import verbosityManager as vbMng, freepar as fp from rrompy.utilities.numerical import customPInv from rrompy.utilities.exception_manager import (RROMPyWarning, RROMPyException, RROMPyAssert) -__all__ = ['RBDistributedPivoted'] +__all__ = ['ReducedBasisPivoted'] -class RBDistributedPivoted(GenericPivotedApproximant): +class ReducedBasisPivoted(GenericPivotedApproximant): """ ROM pivoted RB approximant computation for parametric problems. 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; - 'R': rank for pivot Galerkin projection; defaults to prod(S); - 'PODTolerance': tolerance for pivot snapshots POD; defaults to -1; - '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. homogeneized: Whether to homogeneize Dirichlet BCs. approxRadius: Dummy radius of approximant (i.e. distance from mu0 to farthest sample point). 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; - 'R': rank for Galerkin projection; - 'PODTolerance': tolerance for snapshots POD; - '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. R: Rank for Galerkin projection. PODTolerance: Tolerance for pivot snapshots POD. 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. As: List of sparse matrices (in CSC format) representing coefficients of linear system matrix. bs: List of numpy vectors representing coefficients of linear system RHS. ARBs: List of sparse matrices (in CSC format) representing coefficients of compressed linear system matrix. bRBs: List of numpy vectors representing coefficients of compressed linear system RHS. """ def __init__(self, HFEngine:HFEng, mu0 : paramVal = None, directionPivot : ListAny = [0], approxParameters : DictAny = {}, homogeneized : bool = False, verbosity : int = 10, timestamp : bool = True): self._preInit() self._addParametersToList(["R", "PODTolerance"], [1, -1]) super().__init__(HFEngine = HFEngine, mu0 = mu0, directionPivot = directionPivot, approxParameters = approxParameters, homogeneized = homogeneized, verbosity = verbosity, timestamp = timestamp) self._postInit() @property def S(self): """Value of S.""" return self._S @S.setter def S(self, S): GenericPivotedApproximant.S.fset(self, S) if not hasattr(self, "_R"): self._R = np.prod(self.S) * np.prod(self.SMarginal) @property def R(self): """Value of R. Its assignment may change S.""" return self._R @R.setter def R(self, R): if R < 0: raise RROMPyException("R must be non-negative.") self._R = R self._approxParameters["R"] = self.R if (hasattr(self, "_S") and hasattr(self, "_SMarginal") and np.prod(self.S) * np.prod(self.SMarginal) < self.R): RROMPyWarning(("Prescribed S and SMarginal are too small. " "Decreasing R.")) self.R = np.prod(self.S) * np.prod(self.SMarginal) @property def PODTolerance(self): """Value of PODTolerance.""" return self._PODTolerance @PODTolerance.setter def PODTolerance(self, PODTolerance): self._PODTolerance = PODTolerance self._approxParameters["PODTolerance"] = self.PODTolerance def _setupProjectionMatrix(self): """Compute projection matrix.""" RROMPyAssert(self._mode, message = "Cannot setup numerator.") vbMng(self, "INIT", "Starting computation of projection matrix.", 7) if self.POD: U, s, _ = np.linalg.svd(self.samplingEngine.RPODCoalesced) s = s ** 2. else: Gramian = self.HFEngine.innerProduct( self.samplingEngine.samplesCoalesced, self.samplingEngine.samplesCoalesced) U, s, _ = np.linalg.svd(Gramian) nsamples = self.samplingEngine.samplesCoalesced.shape[1] snorm = np.cumsum(s[::-1]) / np.sum(s) nPODTrunc = min(nsamples - np.argmax(snorm > self.PODTolerance), self.R) pMat = self.samplingEngine.samplesCoalesced.dot(U[:, : nPODTrunc]) vbMng(self, "MAIN", "Assembling {}x{} projection matrix from {} samples.".format( *(pMat.shape), nsamples), 5) vbMng(self, "DEL", "Done computing projection matrix.", 7) return pMat def _setupAffineBlocks(self): """Compute list of marginalized affine blocks of system.""" hasAs = hasattr(self, "AsList") and self.AsList is not None hasbs = hasattr(self, "bsList") and self.bsList is not None if hasAs and hasbs: return vbMng(self, "INIT", "Computing affine blocks of system.", 10) mu0Eff = self.mu0.data[0, self.directionPivot] if not hasAs: self.AsList = [None] * len(self.musMarginal) if not hasbs: self.bsList = [None] * len(self.musMarginal) for k, muM in enumerate(self.musMarginal): muEff = [fp] * self.npar for jj, kk in enumerate(self.directionMarginal): muEff[kk] = muM(0, jj) MEnginek = MarginalProxyEngine(self.HFEngine, muEff) if not hasAs: self.AsList[k] = MEnginek.affineLinearSystemA(mu0Eff, self.scaleFactorPivot) if not hasbs: self.bsList[k] = MEnginek.affineLinearSystemb(mu0Eff, self.scaleFactorPivot, self.homogeneized) vbMng(self, "DEL", "Done computing affine blocks.", 10) def setupApprox(self): """Compute RB projection matrix.""" if self.checkComputedApprox(): return RROMPyAssert(self._mode, message = "Cannot setup approximant.") vbMng(self, "INIT", "Setting up {}.".format(self.name()), 5) self.computeScaleFactor() self.computeSnapshots() self._setupAffineBlocks() pMat = self._setupProjectionMatrix() ARBsList, bRBsList, pList = self.assembleReducedSystemMarginal(pMat) if self.trainedModel is None: self.trainedModel = tModel() self.trainedModel.verbosity = self.verbosity self.trainedModel.timestamp = self.timestamp data = TrainedModelPivotedData(self.trainedModel.name(), self.mu0, pList, self.scaleFactor, self.HFEngine.rescalingExp, self.directionPivot) data.musPivot = copy(self.musPivot) data.musMarginal = copy(self.musMarginal) self.trainedModel.data = data else: self.trainedModel = self.trainedModel self.trainedModel.data.projMat = copy(pList) self.trainedModel.data.marginalInterp = self._setupMarginalInterp() self.trainedModel.data.ARBsList = ARBsList self.trainedModel.data.bRBsList = bRBsList self.trainedModel.data.approxParameters = copy(self.approxParameters) vbMng(self, "DEL", "Done setting up approximant.", 5) def assembleReducedSystemMarginal(self, pMat : sampList = None)\ -> Tuple[List[List[Np2D]], List[List[Np1D]], List[Np2D]]: """Build affine blocks of RB linear system through projections.""" if pMat is None: self.setupApprox() ARBsList = self.trainedModel.data.ARBsList bRBsList = self.trainedModel.data.bRBsList else: vbMng(self, "INIT", "Projecting affine terms of HF model.", 10) ARBsList = [None] * len(self.musMarginal) bRBsList = [None] * len(self.musMarginal) projList = [None] * len(self.musMarginal) for k, (As, bs, sample) in enumerate(zip(self.AsList, self.bsList, self.samplingEngine.samples)): compLocal = self.HFEngine.innerProduct(sample, pMat) projList[k] = sample.dot(customPInv(compLocal)) ARBsList[k], bRBsList[k] = projectAffineDecomposition(As, bs, projList[k]) vbMng(self, "DEL", "Done projecting affine terms.", 10) return ARBsList, bRBsList, projList diff --git a/rrompy/reduction_methods/pole_matching/__init__.py b/rrompy/reduction_methods/pole_matching/__init__.py index ae1e7db..0370cf1 100644 --- a/rrompy/reduction_methods/pole_matching/__init__.py +++ b/rrompy/reduction_methods/pole_matching/__init__.py @@ -1,29 +1,29 @@ # 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_pole_matching_approximant import GenericPoleMatchingApproximant from .rational_interpolant_pole_matching import RationalInterpolantPoleMatching -from .rb_distributed_pole_matching import RBDistributedPoleMatching +from .reduced_basis_pole_matching import ReducedBasisPoleMatching __all__ = [ 'GenericPoleMatchingApproximant', 'RationalInterpolantPoleMatching', - 'RBDistributedPoleMatching' + 'ReducedBasisPoleMatching' ] diff --git a/rrompy/reduction_methods/pole_matching/rational_interpolant_pole_matching.py b/rrompy/reduction_methods/pole_matching/rational_interpolant_pole_matching.py index bc8eeec..736d68b 100644 --- a/rrompy/reduction_methods/pole_matching/rational_interpolant_pole_matching.py +++ b/rrompy/reduction_methods/pole_matching/rational_interpolant_pole_matching.py @@ -1,222 +1,222 @@ # 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 copy import deepcopy as copy import numpy as np from rrompy.reduction_methods.pivoting.rational_interpolant_pivoted import \ RationalInterpolantPivoted from .generic_pole_matching_approximant import GenericPoleMatchingApproximant from rrompy.utilities.poly_fitting.polynomial import ( PolynomialInterpolator as PI) from rrompy.reduction_methods.trained_model import (TrainedModelPivotedData, - TrainedModelPivotedRationalPoleMatching as tModel) + TrainedModelPoleMatchingRational as tModel) from rrompy.utilities.base import verbosityManager as vbMng from rrompy.utilities.exception_manager import RROMPyAssert __all__ = ['RationalInterpolantPoleMatching'] class RationalInterpolantPoleMatching(GenericPoleMatchingApproximant, RationalInterpolantPivoted): """ ROM pivoted rational interpolant computation for parametric problems with pole matching. 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; - 'matchingWeight': weight for pole matching optimization; defaults to 1; - 'cutOffTolerance': tolerance for ignoring parasitic poles; defaults to np.inf; - 'cutOffType': rule for tolerance computation for parasitic poles; defaults to 'MAGNITUDE'; - '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; - 'polybasisPivot': type of polynomial basis for pivot interpolation; allowed values include 'MONOMIAL', 'CHEBYSHEV' and 'LEGENDRE'; defaults to 'MONOMIAL'; - 'polybasisMarginal': type of polynomial basis for marginal interpolation; allowed values include 'MONOMIAL', 'CHEBYSHEV' and 'LEGENDRE'; defaults to 'MONOMIAL'; - 'E': number of derivatives used to compute interpolant; defaults to 0; - 'M': degree of rational interpolant numerator; defaults to 0; - 'N': degree of rational interpolant denominator; defaults to 0; - 'MMarginal': degree of marginal interpolant; defaults to 0; - 'radialBasisPivot': radial basis family for pivot numerator; defaults to 0, i.e. no radial basis; - 'radialBasisMarginal': radial basis family for marginal interpolant; defaults to 0, i.e. no radial basis; - 'radialBasisWeightsPivot': radial basis weights for pivot numerator; defaults to 0, i.e. identity; - 'radialBasisWeightsMarginal': radial basis weights for marginal interpolant; defaults to 0, i.e. identity; - 'interpRcondPivot': tolerance for pivot interpolation; defaults to None; - 'interpRcondMarginal': tolerance for marginal interpolation; defaults to None; - 'robustTol': tolerance for robust rational denominator management; defaults to 0. 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; - 'matchingWeight': weight for pole matching optimization; - 'cutOffTolerance': tolerance for ignoring parasitic poles; - 'cutOffType': rule for tolerance computation for parasitic poles; - 'polybasisPivot': type of polynomial basis for pivot interpolation; - 'polybasisMarginal': type of polynomial basis for marginal interpolation; - 'E': number of derivatives used to compute interpolant; - 'M': degree of rational interpolant numerator; - 'N': degree of rational interpolant denominator; - 'MMarginal': degree of marginal interpolant; - 'radialBasisPivot': radial basis family for pivot numerator; - 'radialBasisMarginal': radial basis family for marginal interpolant; - 'radialBasisWeightsPivot': radial basis weights for pivot numerator; - 'radialBasisWeightsMarginal': radial basis weights for marginal interpolant; - 'interpRcondPivot': tolerance for pivot interpolation; - 'interpRcondMarginal': tolerance for marginal interpolation; - 'robustTol': tolerance for robust rational denominator management. 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. matchingWeight: Weight for pole matching optimization. cutOffTolerance: Tolerance for ignoring parasitic poles. cutOffType: Rule for tolerance computation for parasitic poles. S: Total number of pivot samples current approximant relies upon. sampler: Pivot sample point generator. polybasisPivot: Type of polynomial basis for pivot interpolation. polybasisMarginal: Type of polynomial basis for marginal interpolation. M: Numerator degree of approximant. N: Denominator degree of approximant. MMarginal: Degree of marginal interpolant. radialBasisPivot: Radial basis family for pivot numerator. radialBasisMarginal: Radial basis family for marginal interpolant. radialBasisWeightsPivot: Radial basis weights for pivot numerator. radialBasisWeightsMarginal: Radial basis weights for marginal interpolant. interpRcondPivot: Tolerance for pivot interpolation. interpRcondMarginal: Tolerance for marginal interpolation. robustTol: Tolerance for robust rational denominator management. E: Complete derivative depth level of S. 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. Q: Numpy 1D vector containing complex coefficients of approximant denominator. P: Numpy 2D vector whose columns are FE dofs of coefficients of approximant numerator. """ def setupApprox(self): """ Compute rational interpolant. SVD-based robust eigenvalue management. """ if self.checkComputedApprox(): return RROMPyAssert(self._mode, message = "Cannot setup approximant.") vbMng(self, "INIT", "Setting up {}.". format(self.name()), 5) self.computeScaleFactor() self.computeSnapshots() if self.trainedModel is None: self.trainedModel = tModel() self.trainedModel.verbosity = self.verbosity self.trainedModel.timestamp = self.timestamp data = TrainedModelPivotedData(self.trainedModel.name(), self.mu0, self.samplingEngine.samplesCoalesced, self.scaleFactor, self.HFEngine.rescalingExp, self.directionPivot) data.musPivot = copy(self.musPivot) data.musMarginal = copy(self.musMarginal) self.trainedModel.data = data else: self.trainedModel = self.trainedModel self.trainedModel.data.projMat = copy( self.samplingEngine.samplesCoalesced) self.trainedModel.data.marginalInterp = self._setupMarginalInterp() if self.N > 0: Qs = self._setupDenominator() else: Q = PI() Q.npar = self.musPivot.shape[1] Q.coeffs = np.ones(tuple([1] * Q.npar), dtype = np.complex) Q.polybasis = self.polybasisPivot Qs = [Q for _ in range(len(self.musMarginal))] self.trainedModel.data._temporary = True self.trainedModel.data.Qs = Qs self.trainedModel.data.Ps = self._setupNumerator() vbMng(self, "INIT", "Matching poles.", 10) self.trainedModel.initializeFromRational(self.HFEngine, self.matchingWeight, self.POD) vbMng(self, "DEL", "Done matching poles.", 10) if not np.isinf(self.cutOffTolerance): vbMng(self, "INIT", "Recompressing by cut-off.", 10) msg = self.trainedModel.recompressByCutOff([-1., 1.], self.cutOffTolerance, self.cutOffType) vbMng(self, "DEL", "Done recompressing." + msg, 10) self.trainedModel.data.approxParameters = copy(self.approxParameters) vbMng(self, "DEL", "Done setting up approximant.", 5) diff --git a/rrompy/reduction_methods/pole_matching/rb_distributed_pole_matching.py b/rrompy/reduction_methods/pole_matching/reduced_basis_pole_matching.py similarity index 96% rename from rrompy/reduction_methods/pole_matching/rb_distributed_pole_matching.py rename to rrompy/reduction_methods/pole_matching/reduced_basis_pole_matching.py index d211dd3..d3c7aa4 100644 --- a/rrompy/reduction_methods/pole_matching/rb_distributed_pole_matching.py +++ b/rrompy/reduction_methods/pole_matching/reduced_basis_pole_matching.py @@ -1,190 +1,190 @@ # 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 numpy import isinf from copy import deepcopy as copy -from rrompy.reduction_methods.pivoting.rb_distributed_pivoted import \ - RBDistributedPivoted +from rrompy.reduction_methods.pivoting.reduced_basis_pivoted import \ + ReducedBasisPivoted from .generic_pole_matching_approximant import GenericPoleMatchingApproximant from rrompy.reduction_methods.trained_model import (TrainedModelPivotedData, - TrainedModelPivotedRBPoleMatching as tModel) + TrainedModelPoleMatchingReducedBasis as tModel) from rrompy.utilities.base import verbosityManager as vbMng from rrompy.utilities.exception_manager import RROMPyAssert -__all__ = ['RBDistributedPoleMatching'] +__all__ = ['ReducedBasisPoleMatching'] -class RBDistributedPoleMatching(GenericPoleMatchingApproximant, - RBDistributedPivoted): +class ReducedBasisPoleMatching(GenericPoleMatchingApproximant, + ReducedBasisPivoted): """ ROM pivoted RB approximant computation for parametric problems with pole matching. 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; - 'matchingWeight': weight for pole matching optimization; defaults to 1; - 'cutOffTolerance': tolerance for ignoring parasitic poles; defaults to np.inf; - 'cutOffType': rule for tolerance computation for parasitic poles; defaults to 'MAGNITUDE'; - '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; - 'R': rank for pivot Galerkin projection; defaults to prod(S); - 'PODTolerance': tolerance for pivot snapshots POD; defaults to -1; - '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; - 'matchingWeight': weight for pole matching optimization; - 'cutOffTolerance': tolerance for ignoring parasitic poles; - 'cutOffType': rule for tolerance computation for parasitic poles; - 'R': rank for Galerkin projection; - 'PODTolerance': tolerance for snapshots POD; - '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. matchingWeight: Weight for pole matching optimization. cutOffTolerance: Tolerance for ignoring parasitic poles. cutOffType: Rule for tolerance computation for parasitic poles. 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. R: Rank for Galerkin projection. PODTolerance: Tolerance for pivot snapshots POD. 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. 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. As: List of sparse matrices (in CSC format) representing coefficients of linear system matrix. bs: List of numpy vectors representing coefficients of linear system RHS. ARBs: List of sparse matrices (in CSC format) representing coefficients of compressed linear system matrix. bRBs: List of numpy vectors representing coefficients of compressed linear system RHS. """ def setupApprox(self): """Compute RB projection matrix.""" if self.checkComputedApprox(): return RROMPyAssert(self._mode, message = "Cannot setup approximant.") vbMng(self, "INIT", "Setting up {}.".format(self.name()), 5) self.computeScaleFactor() self.computeSnapshots() self._setupAffineBlocks() pMat = self._setupProjectionMatrix() ARBsList, bRBsList, pList = self.assembleReducedSystemMarginal(pMat) if self.trainedModel is None: self.trainedModel = tModel() self.trainedModel.verbosity = self.verbosity self.trainedModel.timestamp = self.timestamp data = TrainedModelPivotedData(self.trainedModel.name(), self.mu0, pList, self.scaleFactor, self.HFEngine.rescalingExp, self.directionPivot) data.musPivot = copy(self.musPivot) data.musMarginal = copy(self.musMarginal) self.trainedModel.data = data else: self.trainedModel = self.trainedModel self.trainedModel.data.projMat = copy(pList) self.trainedModel.data.marginalInterp = self._setupMarginalInterp() self.trainedModel.data.ARBsList = ARBsList self.trainedModel.data.bRBsList = bRBsList vbMng(self, "INIT", "Matching poles.", 10) self.trainedModel.initializeFromAffine(self.HFEngine, self.matchingWeight, self.POD) vbMng(self, "DEL", "Done matching poles.", 10) if not isinf(self.cutOffTolerance): vbMng(self, "INIT", "Recompressing by cut-off.", 10) msg = self.trainedModel.recompressByCutOff([- 1., 1.], self.cutOffTolerance, self.cutOffType) vbMng(self, "DEL", "Done recompressing." + msg, 10) self.trainedModel.data.approxParameters = copy(self.approxParameters) vbMng(self, "DEL", "Done setting up approximant.", 5) diff --git a/rrompy/reduction_methods/trained_model/__init__.py b/rrompy/reduction_methods/trained_model/__init__.py index e16c89a..a80b481 100644 --- a/rrompy/reduction_methods/trained_model/__init__.py +++ b/rrompy/reduction_methods/trained_model/__init__.py @@ -1,45 +1,47 @@ # 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 .trained_model import TrainedModel from .trained_model_data import TrainedModelData from .trained_model_rational import TrainedModelRational -from .trained_model_rb import TrainedModelRB +from .trained_model_reduced_basis import TrainedModelReducedBasis from .trained_model_pivoted_data import TrainedModelPivotedData from .trained_model_pivoted_rational import TrainedModelPivotedRational -from .trained_model_pivoted_rb import TrainedModelPivotedRB -from .trained_model_pivoted_general_pole_matching import \ - TrainedModelPivotedGeneralPoleMatching -from .trained_model_pivoted_rational_pole_matching import \ - TrainedModelPivotedRationalPoleMatching -from .trained_model_pivoted_rb_pole_matching import \ - TrainedModelPivotedRBPoleMatching +from .trained_model_pivoted_reduced_basis import \ + TrainedModelPivotedReducedBasis +from .trained_model_pole_matching_general import \ + TrainedModelPoleMatchingGeneral +from .trained_model_pole_matching_rational import \ + TrainedModelPoleMatchingRational +from .trained_model_pole_matching_reduced_basis import \ + TrainedModelPoleMatchingReducedBasis __all__ = [ 'TrainedModel', 'TrainedModelData', 'TrainedModelRational', - 'TrainedModelRB', + 'TrainedModelReducedBasis', 'TrainedModelPivotedData', 'TrainedModelPivotedRational', - 'TrainedModelPivotedRB', - 'TrainedModelPivotedRationalPoleMatching', - 'TrainedModelPivotedRBPoleMatching' + 'TrainedModelPivotedReducedBasis', + 'TrainedModelPoleMatchingGeneral', + 'TrainedModelPoleMatchingRational', + 'TrainedModelPoleMatchingReducedBasis' ] diff --git a/rrompy/reduction_methods/trained_model/trained_model_pivoted_rb.py b/rrompy/reduction_methods/trained_model/trained_model_pivoted_reduced_basis.py similarity index 97% rename from rrompy/reduction_methods/trained_model/trained_model_pivoted_rb.py rename to rrompy/reduction_methods/trained_model/trained_model_pivoted_reduced_basis.py index 33c669d..490a1b0 100644 --- a/rrompy/reduction_methods/trained_model/trained_model_pivoted_rb.py +++ b/rrompy/reduction_methods/trained_model/trained_model_pivoted_reduced_basis.py @@ -1,202 +1,202 @@ # 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 .trained_model_rb import TrainedModelRB +from .trained_model_reduced_basis import TrainedModelReducedBasis from rrompy.utilities.base.types import (Np1D, ListAny, paramVal, paramList, sampList) from rrompy.utilities.base import verbosityManager as vbMng, freepar as fp from rrompy.utilities.poly_fitting.polynomial import ( hashIdxToDerivative as hashI) from rrompy.utilities.numerical import (eigvalsNonlinearDense, marginalizePolyList) from rrompy.utilities.exception_manager import RROMPyException from rrompy.parameter import checkParameter, checkParameterList from rrompy.sampling import emptySampleList, sampleList -__all__ = ['TrainedModelPivotedRB'] +__all__ = ['TrainedModelPivotedReducedBasis'] -class TrainedModelPivotedRB(TrainedModelRB): +class TrainedModelPivotedReducedBasis(TrainedModelReducedBasis): """ ROM approximant evaluation for pivoted RB approximant. Attributes: Data: dictionary with all that can be pickled. """ def centerNormalizePivot(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.data.mu0. Returns: Normalized parameter. """ mu = checkParameterList(mu, self.data.nparPivot)[0] if mu0 is None: mu0 = self.data.mu0Pivot rad = ((mu ** self.data.rescalingExpPivot - mu0 ** self.data.rescalingExpPivot) / self.data.scaleFactorPivot) return rad 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.data.mu0. Returns: Normalized parameter. """ mu = checkParameterList(mu, self.data.nparMarginal)[0] if mu0 is None: mu0 = self.data.mu0Marginal rad = ((mu ** self.data.rescalingExpMarginal - mu0 ** self.data.rescalingExpMarginal) / self.data.scaleFactorMarginal) return rad def interpolateMarginal(self, mu : paramVal = [], samples : ListAny = []) -> ListAny: """ Evaluate marginal interpolator at arbitrary marginal parameter. Args: mu: Target parameter. samples: Objects to interpolate. """ mu = checkParameter(mu, self.data.nparMarginal) sList = isinstance(samples[0], sampleList) sEff = [None] * len(samples) for j in range(len(samples)): if sList: sEff[j] = samples[j].data else: sEff[j] = samples[j] vbMng(self, "INIT", "Interpolating marginal at mu = {}.".format(mu), 95) muC = self.centerNormalizeMarginal(mu) p = np.zeros_like(sEff[0], dtype = sEff[0].dtype) for mIj, spj in zip(self.data.marginalInterp, sEff): p += mIj(muC) * spj vbMng(self, "DEL", "Done interpolating marginal.", 95) # if sList: p = sampleList(p) return p def getApproxReduced(self, mu : paramList = []) -> sampList: """ Evaluate reduced representation of approximant at arbitrary parameter. Args: mu: Target parameter. """ mu = checkParameterList(mu, self.data.npar)[0] if (not hasattr(self, "lastSolvedApproxReduced") or self.lastSolvedApproxReduced != mu): vbMng(self, "INIT", "Computing RB solution at mu = {}.".format(mu), 12) self.uApproxReduced = emptySampleList() self.uApproxReduced.reset((self.data.ARBsList[0][0].shape[0], len(mu)), self.data.projMat[0].dtype) for i, muPL in enumerate(mu): muP = checkParameter([muPL(0, x) \ for x in self.data.directionPivot]) muM = [muPL(0, x) for x in self.data.directionMarginal] vbMng(self, "INIT", "Assembling reduced model for mu = {}.".format(muPL), 87) ARBs = self.interpolateMarginal(muM, self.data.ARBsList) bRBs = self.interpolateMarginal(muM, self.data.bRBsList) ARBmu = ARBs[0] bRBmu = bRBs[0] for j in range(1, max(len(ARBs), len(bRBs))): theta = np.prod(self.centerNormalizePivot(muP) ** hashI(j, self.data.nparPivot)) if j < len(ARBs): ARBmu += theta * ARBs[j] if j < len(bRBs): bRBmu += theta * bRBs[j] vbMng(self, "DEL", "Done assembling reduced model.", 87) vbMng(self, "INIT", "Solving reduced model for mu = {}.".format(muPL), 75) self.uApproxReduced[i] = np.linalg.solve(ARBmu, bRBmu) vbMng(self, "DEL", "Done solving reduced model.", 75) vbMng(self, "DEL", "Done computing RB solution.", 12) self.lastSolvedApproxReduced = mu return self.uApproxReduced def getApprox(self, mu : paramList = []) -> sampList: """ Evaluate approximant at arbitrary parameter. Args: mu: Target parameter. """ mu = checkParameterList(mu, self.data.npar)[0] if (not hasattr(self, "lastSolvedApprox") or self.lastSolvedApprox != mu): uApproxR = self.getApproxReduced(mu) self.uApprox = emptySampleList() self.uApprox.reset((self.data.projMat[0].shape[0], len(mu)), self.data.projMat[0].dtype) for i in range(len(mu)): muM = [mu(i, x) for x in self.data.directionMarginal] projLoc = self.interpolateMarginal(muM, self.data.projMat)[0] self.uApprox[i] = projLoc.dot(uApproxR[i]) self.lastSolvedApprox = mu return self.uApprox def getPoles(self, marginalVals : ListAny = [fp], jSupp : int = 1, **kwargs) -> Np1D: """ Obtain approximant poles. Returns: Numpy complex vector of poles. """ if not hasattr(marginalVals, "__len__"): marginalVals = [marginalVals] marginalVals = list(marginalVals) try: rDim = marginalVals.index(fp) if rDim < len(marginalVals) - 1 and fp in marginalVals[rDim + 1 :]: raise except: raise RROMPyException(("Exactly 1 'freepar' entry in " "marginalVals must be provided.")) if rDim in self.data.directionMarginal: raise RROMPyException(("'freepar' entry in marginalVals must be " "in pivot dimension.")) mValsP = [marginalVals[x] for x in self.data.directionPivot] rDimP = mValsP.index(fp) mValsM = [marginalVals[x] for x in self.data.directionMarginal] ARBs = self.interpolateMarginal(mValsM, self.data.ARBsList) if self.data.nparPivot > 1: mValsP[rDimP] = self.data.mu0Pivot(rDimP) mValsP = self.centerNormalizePivot(mValsP) mValsP = mValsP.data.flatten() mValsP[rDimP] = fp ARBs = marginalizePolyList(ARBs, mValsP, "auto") ev = eigvalsNonlinearDense(ARBs, jSupp = jSupp, **kwargs) return np.power(self.data.mu0(rDim) ** self.data.rescalingExp[rDim] + ev * self.scaleFactor[rDim], 1. / self.data.rescalingExp[rDim]) diff --git a/rrompy/reduction_methods/trained_model/trained_model_pivoted_general_pole_matching.py b/rrompy/reduction_methods/trained_model/trained_model_pole_matching_general.py similarity index 98% rename from rrompy/reduction_methods/trained_model/trained_model_pivoted_general_pole_matching.py rename to rrompy/reduction_methods/trained_model/trained_model_pole_matching_general.py index fc1720c..3f088c5 100644 --- a/rrompy/reduction_methods/trained_model/trained_model_pivoted_general_pole_matching.py +++ b/rrompy/reduction_methods/trained_model/trained_model_pole_matching_general.py @@ -1,244 +1,244 @@ # 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 .trained_model import TrainedModel from rrompy.utilities.base.types import (Np1D, Tuple, ListAny, paramVal, paramList, sampList, HFEng) from rrompy.utilities.base import verbosityManager as vbMng, freepar as fp from rrompy.utilities.numerical import pointMatching from rrompy.utilities.poly_fitting.heaviside import HeavisideInterpolator as HI from rrompy.utilities.exception_manager import RROMPyException, RROMPyAssert from rrompy.parameter import checkParameter, checkParameterList from rrompy.sampling import emptySampleList -__all__ = ['TrainedModelPivotedGeneralPoleMatching'] +__all__ = ['TrainedModelPoleMatchingGeneral'] -class TrainedModelPivotedGeneralPoleMatching(TrainedModel): +class TrainedModelPoleMatchingGeneral(TrainedModel): """ ROM approximant evaluation for pivoted approximants with pole matching. Attributes: Data: dictionary with all that can be pickled. """ def initializeFromLists(self, poles:ListAny, coeffs:ListAny, basis:str, HFEngine:HFEng, matchingWeight : float = 1., POD : bool = True): """Initialize Heaviside representation.""" musM = self.data.musMarginal margAbsDist = np.sum(np.abs(np.tile(musM.data.reshape(len(musM),1,-1), [1, len(musM), 1]) - np.tile(musM.data.reshape(1,len(musM),-1), [len(musM), 1, 1])), axis = 2) N = len(poles[0]) explored = [0] unexplored = list(range(1, len(musM))) for _ in range(1, len(musM)): minIdx = np.argmin(np.concatenate([margAbsDist[ex, unexplored] \ for ex in explored])) minIex = explored[minIdx // len(unexplored)] minIunex = unexplored[minIdx % len(unexplored)] dist = np.abs(np.tile(poles[minIex].reshape(-1, 1), N) - poles[minIunex].reshape(1, -1)) if matchingWeight != 0: resex = coeffs[minIex][: N] resunex = coeffs[minIunex][: N] if POD: distR = resex.dot(resunex.T.conj()) distR = (distR.T / np.linalg.norm(resex, axis = 1)).T distR = distR / np.linalg.norm(resunex, axis = 1) else: resex = self.data.projMat.dot(resex.T) resunex = self.data.projMat.dot(resunex.T) distR = HFEngine.innerProduct(resex, resunex).T distR = (distR.T / HFEngine.norm(resex)).T distR = distR / HFEngine.norm(resunex) distR = np.abs(distR) distR[distR > 1.] = 1. dist += 2. / np.pi * matchingWeight * np.arccos(distR) reordering = pointMatching(dist, verbObj = self) poles[minIunex] = poles[minIunex][reordering] coeffs[minIunex][: N] = coeffs[minIunex][reordering] explored += [minIunex] unexplored.remove(minIunex) # explored = np.append(explored, minIunex) # unexplored = unexplored[unexplored != minIunex] HIs = [] for pls, cfs in zip(poles, coeffs): hsi = HI() hsi.poles = pls hsi.coeffs = cfs hsi.npar = 1 hsi.polybasis = basis HIs += [hsi] self.data.HIs = HIs def recompressByCutOff(self, murange : Tuple[float, float] = [- 1., 1.], tol : float = np.inf, rtype : str = "MAGNITUDE"): if np.isinf(tol): return " No poles erased." N = len(self.data.HIs[0].poles) mu0 = np.mean(murange) musig = murange[0] - mu0 if np.isclose(musig, 0.): radius = lambda x: np.abs(x - mu0) else: if rtype == "MAGNITUDE": murdir = (murange[0] - mu0) / np.abs(musig) def radius(x): scalprod = (x - mu0) * murdir.conj() / np.abs(musig) rescalepar = np.abs(np.real(scalprod)) rescaleort = np.abs(np.imag(scalprod)) return ((rescalepar - 1.) ** 2. * (rescalepar > 1.) + rescaleort ** 2.) ** .5 else:#if rtype == "POTENTIAL": def radius(x): rescale = (x - mu0) / musig return np.max(np.abs(rescale * np.array([-1., 1.]) + (rescale ** 2. - 1) ** .5)) - 1. keepPole, removePole = [], [] for j in range(N): for hi in self.data.HIs: if radius(hi.poles[j]) <= tol: keepPole += [j] break if len(keepPole) == 0 or keepPole[-1] != j: removePole += [j] if len(keepPole) == N: return " No poles erased." keepCoeff = keepPole + [N] keepCoeff = keepCoeff + list(range(N + 1,len(self.data.HIs[0].coeffs))) for hi in self.data.HIs: polyCorrection = np.zeros_like(hi.coeffs[0, :]) for j in removePole: polyCorrection += hi.coeffs[j, :] / (mu0 - hi.poles[j]) if len(hi.coeffs) == N: hi.coeffs = np.vstack((hi.coeffs, polyCorrection)) else: hi.coeffs[N, :] += polyCorrection hi.poles = hi.poles[keepPole] hi.coeffs = hi.coeffs[keepCoeff, :] return (" Erased {} pole".format(len(removePole)) + "s" * (len(removePole) > 1) + ".") def interpolateMarginalInterpolator(self, mu : paramVal = []) -> Np1D: """Obtain interpolated approximant interpolator.""" mu = checkParameter(mu, self.data.nparMarginal)[0] hsi = HI() hsi.poles = self.interpolateMarginalPoles(mu) hsi.coeffs = self.interpolateMarginalCoeffs(mu) hsi.npar = 1 hsi.polybasis = self.data.HIs[0].polybasis return hsi def interpolateMarginalPoles(self, mu : paramList = []) -> Np1D: """Obtain interpolated approximant poles.""" mu = checkParameterList(mu, self.data.nparMarginal)[0] return self.interpolateMarginal(mu, [hi.poles for hi in self.data.HIs]) def interpolateMarginalCoeffs(self, mu : paramList = []) -> Np1D: """Obtain interpolated approximant coefficients.""" mu = checkParameterList(mu, self.data.nparMarginal)[0] cs = self.interpolateMarginal(mu, [hi.coeffs for hi in self.data.HIs]) if isinstance(cs, (list, tuple,)): cs = np.array(cs) return cs.reshape(self.data.HIs[0].coeffs.shape) def getApproxReduced(self, mu : paramList = []) -> sampList: """ Evaluate reduced representation of approximant at arbitrary parameter. Args: mu: Target parameter. """ RROMPyAssert(self.data.nparPivot, 1, "Number of pivot parameters") mu = checkParameterList(mu, self.data.npar)[0] if (not hasattr(self, "lastSolvedApproxReduced") or self.lastSolvedApproxReduced != mu): vbMng(self, "INIT", "Evaluating approximant at mu = {}.".format(mu), 12) self.uApproxReduced = emptySampleList() self.uApproxReduced.reset((self.data.HIs[0].coeffs.shape[1], len(mu)), self.data.projMat[0].dtype) for i, muPL in enumerate(mu): muL = self.centerNormalizePivot([muPL(0, x) \ for x in self.data.directionPivot]) muM = [muPL(0, x) for x in self.data.directionMarginal] vbMng(self, "INIT", "Assembling reduced model for mu = {}.".format(muPL), 87) hsL = self.interpolateMarginalInterpolator(muM) vbMng(self, "DEL", "Done assembling reduced model.", 87) self.uApproxReduced[i] = hsL(muL) vbMng(self, "DEL", "Done evaluating approximant.", 12) self.lastSolvedApproxReduced = mu return self.uApproxReduced def getPoles(self, *args, **kwargs) -> Np1D: """ Obtain approximant poles. Returns: Numpy complex vector of poles. """ RROMPyAssert(self.data.nparPivot, 1, "Number of pivot parameters") if len(args) + len(kwargs) > 1: raise RROMPyException(("Wrong number of parameters passed. " "Only 1 available.")) elif len(args) + len(kwargs) == 1: if len(args) == 1: mVals = args[0] else: mVals = kwargs["marginalVals"] if not hasattr(mVals, "__len__"): mVals = [mVals] mVals = list(mVals) else: mVals = [fp] try: rDim = mVals.index(fp) if rDim < len(mVals) - 1 and fp in mVals[rDim + 1 :]: raise except: raise RROMPyException(("Exactly 1 'freepar' entry in " "marginalVals must be provided.")) if rDim != self.data.directionPivot[0]: raise RROMPyException(("'freepar' entry in marginalVals must " "coincide with pivot direction.")) mVals[rDim] = self.data.mu0(rDim) mMarg = [mVals[j] for j in range(len(mVals)) if j != rDim] roots = np.array(self.interpolateMarginalPoles(mMarg)) return np.power(self.data.mu0(rDim) ** self.data.rescalingExp[rDim] + self.data.scaleFactor[rDim] * roots, 1. / self.data.rescalingExp[rDim]) def getResidues(self, *args, **kwargs) -> Np1D: """ Obtain approximant residues. Returns: Numpy matrix with residues as columns. """ pls = self.getPoles(*args, **kwargs) if len(args) == 1: mVals = args[0] else: mVals = kwargs["marginalVals"] if not hasattr(mVals, "__len__"): mVals = [mVals] mVals = list(mVals) rDim = mVals.index(fp) mMarg = [mVals[j] for j in range(len(mVals)) if j != rDim] residues = self.interpolateMarginalCoeffs(mMarg)[: len(pls)] res = self.data.projMat.dot(residues) return pls, res diff --git a/rrompy/reduction_methods/trained_model/trained_model_pivoted_rational_pole_matching.py b/rrompy/reduction_methods/trained_model/trained_model_pole_matching_rational.py similarity index 92% rename from rrompy/reduction_methods/trained_model/trained_model_pivoted_rational_pole_matching.py rename to rrompy/reduction_methods/trained_model/trained_model_pole_matching_rational.py index a7e9fba..c24171f 100644 --- a/rrompy/reduction_methods/trained_model/trained_model_pivoted_rational_pole_matching.py +++ b/rrompy/reduction_methods/trained_model/trained_model_pole_matching_rational.py @@ -1,110 +1,109 @@ # 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 scipy.special import factorial as fact from itertools import combinations -from .trained_model_pivoted_general_pole_matching import \ - TrainedModelPivotedGeneralPoleMatching +from .trained_model_pole_matching_general import \ + TrainedModelPoleMatchingGeneral from .trained_model_pivoted_rational import TrainedModelPivotedRational from rrompy.utilities.base.types import Np1D, List, paramList, sampList, HFEng from rrompy.utilities.poly_fitting.heaviside import rational2heaviside from rrompy.utilities.exception_manager import RROMPyException, RROMPyAssert from rrompy.parameter import checkParameterList from rrompy.sampling import emptySampleList -__all__ = ['TrainedModelPivotedRationalPoleMatching'] +__all__ = ['TrainedModelPoleMatchingRational'] -class TrainedModelPivotedRationalPoleMatching( - TrainedModelPivotedGeneralPoleMatching, - TrainedModelPivotedRational): +class TrainedModelPoleMatchingRational(TrainedModelPoleMatchingGeneral, + TrainedModelPivotedRational): """ ROM approximant evaluation for pivoted Rational approximant with pole matching. Attributes: Data: dictionary with all that can be pickled. """ def initializeFromRational(self, HFEngine:HFEng, matchingWeight : float = 1., POD : bool = True): """Initialize Heaviside representation.""" RROMPyAssert(self.data.nparPivot, 1, "Number of pivot parameters") poles, coeffs = [], [] for Q, P in zip(self.data.Qs, self.data.Ps): cfs, pls, basis = rational2heaviside(P, Q) poles += [pls] coeffs += [cfs] self.initializeFromLists(poles, coeffs, basis, HFEngine, matchingWeight, POD) self.data._temporary = False def getPVal(self, mu : paramList = []) -> sampList: """ Evaluate rational numerator at arbitrary parameter. Args: mu: Target parameter. """ if self.data._temporary: return super().getPVal(mu) RROMPyAssert(self.data.nparPivot, 1, "Number of pivot parameters") mu = checkParameterList(mu, self.data.npar)[0] p = emptySampleList() p.reset((len(self.data.HIs[0].coeffs.shape[1]), len(mu))) for i, muPL in enumerate(mu): muL = self.centerNormalizePivot([muPL(0, x) \ for x in self.data.directionPivot]) muM = [muPL(0, x) for x in self.data.directionMarginal] hsL = self.interpolateMarginalInterpolator(muM) p[i] = hsL(muL) * np.prod(muL(0, 0) - hsL.poles) return p def getQVal(self, mu:Np1D, der : List[int] = None, scl : Np1D = None) -> Np1D: """ Evaluate rational denominator at arbitrary parameter. Args: mu: Target parameter. der(optional): Derivatives to take before evaluation. """ if self.data._temporary: return super().getQVal(mu, der, scl) RROMPyAssert(self.data.nparPivot, 1, "Number of pivot parameters") mu = checkParameterList(mu, self.data.npar)[0] muP = self.centerNormalizePivot(checkParameterList( mu.data[:, self.data.directionPivot], 1)[0]) muM = checkParameterList(mu.data[:, self.data.directionMarginal], self.data.npar - 1)[0] if der is None: derP, derM = 0, [0] else: derP = der[self.data.directionPivot[0]] derM = [der[x] for x in self.data.directionMarginal] if np.any(np.array(derM) != 0): raise RROMPyException(("Derivatives of Q with respect to marginal " "parameters not allowed.")) sclP = 1 if scl is None else scl[self.data.directionPivot[0]] derVal = np.zeros(len(mu), dtype = self.data.HIs[0].poles.dtype) N = len(self.data.HIs[0].poles) if derP == N: derVal[:] = 1. elif derP >= 0 and derP < N: pls = self.interpolateMarginalPoles(muM).reshape(-1, len(mu)).T plsDist = muP.data.reshape(-1, 1) - pls for terms in combinations(np.arange(N), N - derP): derVal += np.prod(plsDist[:, list(terms)], axis = 1) return sclP ** derP * fact(derP) * derVal diff --git a/rrompy/reduction_methods/trained_model/trained_model_pivoted_rb_pole_matching.py b/rrompy/reduction_methods/trained_model/trained_model_pole_matching_reduced_basis.py similarity index 81% rename from rrompy/reduction_methods/trained_model/trained_model_pivoted_rb_pole_matching.py rename to rrompy/reduction_methods/trained_model/trained_model_pole_matching_reduced_basis.py index f2cf22d..23802b1 100644 --- a/rrompy/reduction_methods/trained_model/trained_model_pivoted_rb_pole_matching.py +++ b/rrompy/reduction_methods/trained_model/trained_model_pole_matching_reduced_basis.py @@ -1,55 +1,56 @@ # 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 .trained_model_pivoted_general_pole_matching import \ - TrainedModelPivotedGeneralPoleMatching -from .trained_model_pivoted_rb import TrainedModelPivotedRB +from .trained_model_pole_matching_general import \ + TrainedModelPoleMatchingGeneral +from .trained_model_pivoted_reduced_basis import \ + TrainedModelPivotedReducedBasis from rrompy.utilities.base.types import HFEng from rrompy.utilities.poly_fitting.heaviside import affine2heaviside from rrompy.utilities.exception_manager import RROMPyAssert -__all__ = ['TrainedModelPivotedRBPoleMatching'] +__all__ = ['TrainedModelPoleMatchingReducedBasis'] -class TrainedModelPivotedRBPoleMatching(TrainedModelPivotedGeneralPoleMatching, - TrainedModelPivotedRB): +class TrainedModelPoleMatchingReducedBasis(TrainedModelPoleMatchingGeneral, + TrainedModelPivotedReducedBasis): """ ROM approximant evaluation for pivoted RB approximant with pole matching. Attributes: Data: dictionary with all that can be pickled. """ def initializeFromAffine(self, HFEngine:HFEng, matchingWeight : float = 1., POD : bool = True, jSupp : int = 1): """Initialize Heaviside representation.""" RROMPyAssert(self.data.nparPivot, 1, "Number of pivot parameters") poles, coeffs = [], [] nbadpls = 0 for As, bs in zip(self.data.ARBsList, self.data.bRBsList): cfs, pls, basis = affine2heaviside(As, bs, jSupp) poles += [pls] coeffs += [cfs] nbadpls = max(nbadpls, np.sum(np.isinf(pls))) if nbadpls > 0: for j in range(len(poles)): plsgood = np.argsort(np.abs(pls))[: - nbadpls] poles[j] = poles[j][plsgood] coeffs[j] = coeffs[j][plsgood, :] self.initializeFromLists(poles, coeffs, basis, HFEngine, matchingWeight, POD) diff --git a/rrompy/reduction_methods/trained_model/trained_model_rb.py b/rrompy/reduction_methods/trained_model/trained_model_reduced_basis.py similarity index 98% rename from rrompy/reduction_methods/trained_model/trained_model_rb.py rename to rrompy/reduction_methods/trained_model/trained_model_reduced_basis.py index bd5ad7f..6754e66 100644 --- a/rrompy/reduction_methods/trained_model/trained_model_rb.py +++ b/rrompy/reduction_methods/trained_model/trained_model_reduced_basis.py @@ -1,108 +1,108 @@ # 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 .trained_model import TrainedModel from rrompy.utilities.base.types import Np1D, ListAny, paramList, sampList from rrompy.utilities.base import verbosityManager as vbMng, freepar as fp from rrompy.utilities.poly_fitting.polynomial import ( hashIdxToDerivative as hashI) from rrompy.utilities.numerical import (eigvalsNonlinearDense, marginalizePolyList) from rrompy.utilities.exception_manager import RROMPyException from rrompy.parameter import checkParameter, checkParameterList from rrompy.sampling import emptySampleList -__all__ = ['TrainedModelRB'] +__all__ = ['TrainedModelReducedBasis'] -class TrainedModelRB(TrainedModel): +class TrainedModelReducedBasis(TrainedModel): """ ROM approximant evaluation for RB approximant. Attributes: Data: dictionary with all that can be pickled. """ def getApproxReduced(self, mu : paramList = []) -> sampList: """ Evaluate reduced representation of approximant at arbitrary parameter. Args: mu: Target parameter. """ mu = checkParameterList(mu, self.data.npar)[0] if (not hasattr(self, "lastSolvedApproxReduced") or self.lastSolvedApproxReduced != mu): vbMng(self, "INIT", "Computing RB solution at mu = {}.".format(mu), 12) ARBs, bRBs = self.data.ARBs, self.data.bRBs self.uApproxReduced = emptySampleList() self.uApproxReduced.reset((ARBs[0].shape[0], len(mu)), self.data.projMat.dtype) mu0Eff = (self.data.mu0 ** self.data.rescalingExp).data for i in range(len(mu)): vbMng(self, "INIT", "Assembling reduced model for mu = {}."\ .format(mu[i]), 17) ARBmu = 1. * ARBs[0] bRBmu = 1. * bRBs[0] for j in range(1, max(len(ARBs), len(bRBs))): derIdx = hashI(j, self.data.npar) theta = np.prod(((mu[i] ** self.data.rescalingExp - mu0Eff) / self.data.scaleFactor) ** derIdx) if j < len(ARBs): ARBmu += theta * ARBs[j] if j < len(bRBs): bRBmu += theta * bRBs[j] vbMng(self, "DEL", "Done assembling reduced model.", 17) vbMng(self, "INIT", "Solving reduced model for mu = {}.".format(mu[i]), 15) self.uApproxReduced[i] = np.linalg.solve(ARBmu, bRBmu) vbMng(self, "DEL", "Done solving reduced model.", 15) vbMng(self, "DEL", "Done computing RB solution.", 12) self.lastSolvedApproxReduced = mu return self.uApproxReduced def getPoles(self, marginalVals : ListAny = [fp], jSupp : int = 1, **kwargs) -> Np1D: """ Obtain approximant poles. Returns: Numpy complex vector of poles. """ if not hasattr(marginalVals, "__len__"): marginalVals = [marginalVals] mVals = list(marginalVals) try: rDim = mVals.index(fp) if rDim < len(mVals) - 1 and fp in mVals[rDim + 1 :]: raise except: raise RROMPyException(("Exactly 1 'freepar' entry in " "marginalVals must be provided.")) ARBs = self.data.ARBs if self.data.npar > 1: mVals[rDim] = self.data.mu0(rDim) mVals = (checkParameter(mVals) ** self.data.rescalingExp - self.data.mu0 ** self.data.rescalingExp) mVals = mVals.data.flatten() mVals[rDim] = fp ARBs = marginalizePolyList(ARBs, mVals, "auto") ev = eigvalsNonlinearDense(ARBs, jSupp = jSupp, **kwargs) return np.power(self.data.mu0(rDim) ** self.data.rescalingExp[rDim] + ev * self.data.scaleFactor[rDim], 1. / self.data.rescalingExp[rDim]) diff --git a/rrompy/utilities/poly_fitting/heaviside/heaviside_interpolator.py b/rrompy/utilities/poly_fitting/heaviside/heaviside_interpolator.py index 79334ef..f0aca96 100644 --- a/rrompy/utilities/poly_fitting/heaviside/heaviside_interpolator.py +++ b/rrompy/utilities/poly_fitting/heaviside/heaviside_interpolator.py @@ -1,78 +1,77 @@ # 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 rrompy.utilities.base.types import (List, ListAny, Np1D, paramList, interpEng) from rrompy.utilities.base import freepar as fp from rrompy.utilities.poly_fitting.polynomial.polynomial_interpolator import ( PolynomialInterpolator) from rrompy.utilities.poly_fitting.polynomial.roots import polyroots from rrompy.utilities.poly_fitting.heaviside.val import polyval from rrompy.utilities.poly_fitting.heaviside.heaviside_to_from_rational import\ heaviside2rational, rational2heaviside from rrompy.utilities.exception_manager import RROMPyAssert __all__ = ['HeavisideInterpolator'] class HeavisideInterpolator(PolynomialInterpolator): """HERE""" def __init__(self, other = None): if other is None: return self.poles = other.poles super().__init__(other) def __call__(self, mu:paramList, der : List[int] = None, scl : Np1D = None): return polyval(mu, self.coeffs, self.poles, self.polybasis) def __copy__(self): return HeavisideInterpolator(self) def __deepcopy__(self, memo): other = HeavisideInterpolator() other.poles, other.coeffs, other.npar, other.polybasis = copy( (self.poles, self.coeffs, self.npar, self.polybasis), memo) return other def pad(self, nleft : List[int] = None, nright : List[int] = None): if nleft is None: nleft = [0] * len(self.shape) if nright is None: nright = [0] * len(self.shape) if not hasattr(nleft, "__len__"): nleft = [nleft] if not hasattr(nright, "__len__"): nright = [nright] RROMPyAssert(nleft[0], 0, "Padding in free direction") super().pad(nleft, nright) def setupFromRational(self, num:interpEng, den:interpEng, murange : Np1D = np.array([-1., 1.]), - scl : Np1D = None, scaling : List[callable] = None, - scalingInv : List[callable] = None): + scl : Np1D = None, scalingExp : List[float] = None): self.coeffs, self.poles, self.polybasis = rational2heaviside(num, den, - murange, scl, scaling, scalingInv) + murange, scl, + scalingExp) def roots(self, marginalVals : ListAny = [fp], murange : Np1D = None, - scaling : List[callable] = None, - scalingInv : List[callable] = None): + scalingExp : List[float] = None): RROMPyAssert(self.shape, (1,), "Shape of output") RROMPyAssert(marginalVals, [fp], "Marginal values") basisN = self.polybasis.split("_")[0] coeffsN = heaviside2rational(self.coeffs, self.poles, murange, basisN, - scaling = scaling, scalingInv = scalingInv)[0] + scalingExp = scalingExp)[0] return polyroots(coeffsN, basisN) diff --git a/rrompy/utilities/poly_fitting/heaviside/heaviside_to_from_rational.py b/rrompy/utilities/poly_fitting/heaviside/heaviside_to_from_rational.py index 697ac43..89f7e3d 100644 --- a/rrompy/utilities/poly_fitting/heaviside/heaviside_to_from_rational.py +++ b/rrompy/utilities/poly_fitting/heaviside/heaviside_to_from_rational.py @@ -1,99 +1,97 @@ # 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 rrompy.utilities.poly_fitting.polynomial.polynomial_interpolator import \ PolynomialInterpolator from rrompy.utilities.poly_fitting.polynomial.vander import polyvander from rrompy.utilities.poly_fitting.custom_fit import customFit from rrompy.utilities.base.types import Np1D, Np2D, List, Tuple, interpEng from rrompy.parameter.parameter_sampling import (RandomSampler as RS, QuadratureSampler as QS) from .val import polyval from rrompy.utilities.exception_manager import RROMPyException, RROMPyAssert __all__ = ['heaviside2rational', 'rational2heaviside'] def heaviside2rational(c:Np2D, poles:Np1D, murange : Np1D = None, basis : str = "MONOMIAL_HEAVISIDE", basisD : str = None, - scaling : List[callable] = None, - scalingInv : List[callable] = None) \ + scalingExp : List[float] = None) \ -> Tuple[interpEng, interpEng]: num = PolynomialInterpolator() den = PolynomialInterpolator() basisN = basis.split("_")[0] if basisD is None: basisD = basisN if murange is None: multiplier = [1., 1.j] avgs = [np.mean(np.real(poles)), np.mean(np.imag(poles))] ranges = np.array([[np.min(np.real(poles)), np.max(np.real(poles))], [np.min(np.imag(poles)), np.max(np.imag(poles))]]) domIdx = np.argmax(ranges[:, 1] - ranges[:, 0]) murange = [multiplier[domIdx] * x + multiplier[1 - domIdx] * avgs[1 - domIdx] for x in ranges[domIdx, :]] extraPt = None while extraPt is None or np.any(np.isclose(extraPt, poles)): extraPt = murange[0] + (murange[1] - murange[0]) * np.random.rand(1) denAuxPts = np.concatenate((poles, extraPt)) denAuxVals = np.concatenate((np.zeros(len(poles)), [1.])) den.setupByInterpolation(denAuxPts, denAuxVals, len(poles), basisD) den.coeffs /= np.linalg.norm(den.coeffs) if basis == "CHEBYSHEV": - sampler = QS(murange, "CHEBYSHEV", scaling, scalingInv) + sampler = QS(murange, "CHEBYSHEV", scalingExp) elif basis == "LEGENDRE": - sampler = QS(murange, "GAUSSLEGENDRE", scaling, scalingInv) + sampler = QS(murange, "GAUSSLEGENDRE", scalingExp) else: - sampler = RS(murange, "HALTON", scaling, scalingInv) + sampler = RS(murange, "HALTON", scalingExp) xAux = sampler.generatePoints(len(c)) valsAux = den(xAux) * polyval(xAux, c, poles, basis) num.setupByInterpolation(xAux, valsAux, len(c) - 1, basisN) return num, den def rational2heaviside(num:interpEng, den:interpEng, murange : Np1D = np.array([-1., 1.]), scl : Np1D = None, - scaling : List[callable] = None, - scalingInv : List[callable] = None) \ + scalingExp : List[float] = None) \ -> Tuple[Np2D, Np1D, str]: if (not isinstance(num, PolynomialInterpolator) or not isinstance(den, PolynomialInterpolator)): raise RROMPyException(("Rational numerator and denominator must be " "polynomial interpolators.")) RROMPyAssert(num.npar, 1, "Number of parameters") RROMPyAssert(den.npar, 1, "Number of parameters") basis = num.polybasis + "_HEAVISIDE" c = np.zeros_like(num.coeffs) poles = den.roots() Psp = num(poles) Qspder = den(poles, 1, scl) c[: len(poles)] = (Psp / Qspder).T if len(c) > len(poles): from rrompy.parameter.parameter_sampling import (RandomSampler as RS, QuadratureSampler as QS) if num.polybasis == "CHEBYSHEV": - sampler = QS(murange, "CHEBYSHEV", scaling, scalingInv) + sampler = QS(murange, "CHEBYSHEV", scalingExp) elif num.polybasis == "LEGENDRE": - sampler = QS(murange, "GAUSSLEGENDRE", scaling, scalingInv) + sampler = QS(murange, "GAUSSLEGENDRE", scalingExp) else: - sampler = RS(murange, "HALTON", scaling, scalingInv) + sampler = RS(murange, "HALTON", scalingExp) xAux = sampler.generatePoints(len(c)) valsAux = (num(xAux) / den(xAux) - polyval(xAux, c, poles, basis)).T VanAux = polyvander(xAux, [len(c) - len(poles) - 1], num.polybasis) c[len(poles) :] = customFit(VanAux, valsAux) return c, poles, basis diff --git a/setup.py b/setup.py index 32b1d45..9a05576 100644 --- a/setup.py +++ b/setup.py @@ -1,52 +1,52 @@ -# Copyright (C) 2015-2018 by the RBniCS authors +# Copyright (C) 2018 by the RROMPy authors # -# This file is part of RBniCS. +# This file is part of RROMPy. # -# RBniCS is free software: you can redistribute it and/or modify +# 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. # -# RBniCS is distributed in the hope that it will be useful, +# 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 RBniCS. If not, see . +# along with RROMPy. If not, see . # import os from setuptools import find_packages, setup rrompy_directory = os.path.abspath(os.path.dirname(os.path.realpath(__file__))) #rrompy_directory = os.path.join(rrompy_directory, 'rrompy') setup(name="RROMPy", description="Rational reduced order modelling in Python", long_description="Rational reduced order modelling in Python", author="Davide Pradovera", author_email="davide.pradovera@epfl.ch", version="1.6", license="GNU Library or Lesser General Public License (LGPL)", classifiers=[ "Development Status :: 1 - Planning" "Intended Audience :: Developers", "Intended Audience :: Science/Research", "Programming Language :: Python :: 3", "Programming Language :: Python :: 3.4", "Programming Language :: Python :: 3.5", "Programming Language :: Python :: 3.6", "License :: OSI Approved :: GNU Library or Lesser General Public License (LGPL)", "Topic :: Scientific/Engineering :: Mathematics", "Topic :: Software Development :: Libraries :: Python Modules", ], packages=find_packages(rrompy_directory), setup_requires=[ "pytest-runner" ], tests_require=[ "pytest" ], zip_safe=False ) diff --git a/tests/reduction_methods_1D/rb_distributed_1d.py b/tests/reduction_methods_1D/reduced_basis_distributed_1d.py similarity index 89% rename from tests/reduction_methods_1D/rb_distributed_1d.py rename to tests/reduction_methods_1D/reduced_basis_distributed_1d.py index 3b1f288..6905e65 100644 --- a/tests/reduction_methods_1D/rb_distributed_1d.py +++ b/tests/reduction_methods_1D/reduced_basis_distributed_1d.py @@ -1,57 +1,57 @@ # 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 matrix_fft import matrixFFT -from rrompy.reduction_methods.distributed import RBDistributed as RBD +from rrompy.reduction_methods.distributed import ReducedBasis as RB from rrompy.parameter.parameter_sampling import (QuadratureSampler as QS, ManualSampler as MS) from rrompy.parameter import checkParameterList def test_LS(): solver = matrixFFT() params = {"POD": True, "R": 5, "S": 10, "sampler": QS([1., 7.], "CHEBYSHEV")} - approx = RBD(solver, 4., params, verbosity = 0) + approx = RB(solver, 4., params, verbosity = 0) approx.setupApprox() for mu in approx.mus: assert not np.isclose(approx.normErr(mu)[0], 0., atol = 1e-7) approx.POD = False approx.setupApprox() for mu in approx.mus[approx.R :]: assert not np.isclose(approx.normErr(mu)[0], 0., atol = 1e-3) def test_interp(): solver = matrixFFT() params = {"POD": False, "S": 10, "sampler": QS([1., 7.], "CHEBYSHEV")} - approx = RBD(solver, 4., params, verbosity = 0) + approx = RB(solver, 4., params, verbosity = 0) approx.setupApprox() for mu in approx.mus: assert np.isclose(approx.normErr(mu)[0], 0., atol = 1e-7) def test_hermite(): mu = 1.5 solver = matrixFFT() sampler0 = QS([1., 7.], "CHEBYSHEV") points, _ = checkParameterList(np.tile(sampler0.generatePoints(4)(0), 3)) params = {"POD": True, "S": 12, "sampler": MS([1., 7.], points = points)} - approx = RBD(solver, 4., params, verbosity = 0) + approx = RB(solver, 4., params, verbosity = 0) approx.setupApprox() for mu in approx.mus: assert np.isclose(approx.normErr(mu)[0], 0., atol = 1e-8) diff --git a/tests/reduction_methods_1D/rb_distributed_greedy_1d.py b/tests/reduction_methods_1D/reduced_basis_distributed_greedy_1d.py similarity index 91% rename from tests/reduction_methods_1D/rb_distributed_greedy_1d.py rename to tests/reduction_methods_1D/reduced_basis_distributed_greedy_1d.py index 512a049..e0a2462 100644 --- a/tests/reduction_methods_1D/rb_distributed_greedy_1d.py +++ b/tests/reduction_methods_1D/reduced_basis_distributed_greedy_1d.py @@ -1,60 +1,60 @@ # 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 matrix_fft import matrixFFT -from rrompy.reduction_methods.distributed_greedy import RBDistributedGreedy \ - as RBDG +from rrompy.reduction_methods.distributed_greedy import ReducedBasisGreedy \ + as RBG from rrompy.parameter.parameter_sampling import QuadratureSampler as QS def test_lax_tolerance(capsys): mu = 2.25 solver = matrixFFT() params = {"POD": True, "sampler": QS([1.5, 6.5], "UNIFORM"), "S": 4, "greedyTol": 1e-2, "trainSetGenerator": QS([1.5, 6.5], "CHEBYSHEV")} - approx = RBDG(solver, 4., params, verbosity = 10) + approx = RBG(solver, 4., params, verbosity = 10) approx.greedy() out, err = capsys.readouterr() assert "Done computing snapshots (final snapshot count: 10)." in out assert len(err) == 0 assert len(approx.mus) == 10 _, _, maxEst = approx.getMaxErrorEstimator(approx.muTest) assert maxEst < 1e-2 assert np.isclose(approx.normErr(mu)[0], 1.5056e-05, rtol = 1e-1) def test_samples_at_poles(): solver = matrixFFT() params = {"POD": True, "sampler": QS([1.5, 6.5], "UNIFORM"), "S": 4, "nTestPoints": 100, "greedyTol": 1e-5, "trainSetGenerator": QS([1.5, 6.5], "CHEBYSHEV")} - approx = RBDG(solver, 4., params, verbosity = 0) + approx = RBG(solver, 4., params, verbosity = 0) approx.greedy() for mu in approx.mus: assert np.isclose(approx.normErr(mu)[0] / (1e-15+approx.normHF(mu)[0]), 0., atol = 1e-4) poles = approx.getPoles() for lambda_ in range(2, 7): assert np.isclose(np.min(np.abs(poles - lambda_)), 0., atol = 1e-3) assert np.isclose(np.min(np.abs(np.array(approx.mus(0)) - lambda_)), 0., atol = 1e-1) diff --git a/tests/reduction_methods_multiD/rb_distributed_2d.py b/tests/reduction_methods_multiD/reduced_basis_distributed_2d.py similarity index 90% rename from tests/reduction_methods_multiD/rb_distributed_2d.py rename to tests/reduction_methods_multiD/reduced_basis_distributed_2d.py index 8182e2d..2008338 100644 --- a/tests/reduction_methods_multiD/rb_distributed_2d.py +++ b/tests/reduction_methods_multiD/reduced_basis_distributed_2d.py @@ -1,62 +1,62 @@ # 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 matrix_random import matrixRandom -from rrompy.reduction_methods.distributed import RBDistributed as RBD +from rrompy.reduction_methods.distributed import ReducedBasis as RB from rrompy.parameter.parameter_sampling import (QuadratureSampler as QS, ManualSampler as MS) from rrompy.parameter import checkParameterList def test_LS(): mu0 = [2, 3] solver = matrixRandom() params = {"POD": True, "R": 5, "S": [3, 3], "sampler": QS([[0., 4.], [1., 5.]], "CHEBYSHEV")} - approx = RBD(solver, mu0, params, verbosity = 0) + approx = RB(solver, mu0, params, verbosity = 0) approx.setupApprox() for mu in approx.mus: assert not np.isclose(approx.normErr(mu)[0], 0., atol = 1e-7) approx.POD = False approx.setupApprox() for mu in approx.mus[approx.R :]: assert not np.isclose(approx.normErr(mu)[0], 0., atol = 1e-3) def test_interp(): mu0 = [2, 3] solver = matrixRandom() params = {"POD": False, "S": [3, 3], "sampler": QS([[0., 4.], [1., 5.]], "CHEBYSHEV")} - approx = RBD(solver, mu0, params, verbosity = 0) + approx = RB(solver, mu0, params, verbosity = 0) approx.setupApprox() for mu in approx.mus: assert np.isclose(approx.normErr(mu)[0], 0., atol = 1e-7) def test_hermite(): mu0 = [2, 3] solver = matrixRandom() sampler0 = QS([[0., 4.], [1., 5.]], "CHEBYSHEV") points, _ = checkParameterList(np.tile( sampler0.generatePoints([2, 2]).data, [3, 1])) params = {"POD": True, "S": [12], "sampler": MS([[0., 4.], [1., 5.]], points = points)} - approx = RBD(solver, mu0, params, verbosity = 0) + approx = RB(solver, mu0, params, verbosity = 0) approx.setupApprox() for mu in approx.mus: assert np.isclose(approx.normErr(mu)[0], 0., atol = 1e-8) diff --git a/tests/utilities/parameter_sampling.py b/tests/utilities/parameter_sampling.py index b94792e..b1ea088 100644 --- a/tests/utilities/parameter_sampling.py +++ b/tests/utilities/parameter_sampling.py @@ -1,60 +1,59 @@ # 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 rrompy.parameter.parameter_sampling import (ManualSampler, QuadratureSampler, RandomSampler, FFTSampler) from rrompy.parameter import checkParameter def test_manual(): sampler = ManualSampler(lims = [0., 3.], points = np.linspace(0, 3, 101), - scaling = lambda x: np.power(x, 2.), - scalingInv = lambda x: np.power(x, .5)) + scalingExp = 2.) assert sampler.name() == "ManualSampler" x = sampler.generatePoints(10) assert np.allclose(x(0), np.linspace(0, 3, 101)[:10], rtol = 1e-5) def test_quadrature(): sampler = QuadratureSampler(lims = [0., 3.], kind = "CHEBYSHEV") x = sampler.generatePoints(9) assert np.isclose(x(0)[2], 1.5, rtol = 1e-5) def test_random(): sampler = RandomSampler(lims = [0., 3.], kind = "SOBOL", seed = 13432) x = sampler.generatePoints(100) assert np.isclose(x(0)[47], 0.55609130859375, rtol = 1e-5) def test_fft(): sampler = FFTSampler(lims = [-1., 1.]) x = sampler.generatePoints(100) assert np.allclose(np.power(x(0), 100), 1., rtol = 1e-5) def test_2D(): sampler = QuadratureSampler(lims = [(0., 0.), (3., 1.)], kind = "GAUSSLEGENDRE") x = sampler.generatePoints((9, 5)) assert sum(np.isclose(x(0), 1.5)) == 5 assert sum(np.isclose(x(1), .5)) == 9 def test_4D(): sampler = RandomSampler(lims = [tuple([0.] * 4), tuple([1.] * 4)], kind = "UNIFORM", seed = 1234) x = sampler.generatePoints(10) assert x.shape[1] == 4 assert checkParameter([x[0]]) == checkParameter([(0.191519450378892, 0.622108771039832, 0.437727739007115, 0.785358583713769)])