diff --git a/examples/5_anisotropic_square/anisotropic_square.py b/examples/5_anisotropic_square/anisotropic_square.py index 136baa4..2dd8934 100644 --- a/examples/5_anisotropic_square/anisotropic_square.py +++ b/examples/5_anisotropic_square/anisotropic_square.py @@ -1,80 +1,80 @@ ### example from Smetana, Zahm, Patera. Randomized residual-based error ### estimators for parametrized equations. import numpy as np import matplotlib.pyplot as plt from itertools import product from anisotropic_square_engine import (AnisotropicSquareEngine as engine, AnisotropicSquareEnginePoles as plsEx) from rrompy.reduction_methods import ( RationalInterpolantGreedyPivotedGreedyPoleMatch as RIGPG) from rrompy.parameter.parameter_sampling import (QuadratureSampler as QS, SparseGridSampler as SGS) zs, Ls = [10., 50.], [.2, 1.2] z0, L0, n = np.mean(zs), np.mean(Ls), 50 murange = [[zs[0], Ls[0]], [zs[-1], Ls[-1]]] np.random.seed(4020) mu = [zs[0] + np.random.rand() * (zs[-1] - zs[0]), Ls[0] + np.random.rand() * (Ls[-1] - Ls[0])] solver = engine(z0, L0, n) fighandles = [] params = {"POD": True, "nTestPoints": 100, "greedyTol": 1e-4, "S": 3, "polybasisMarginal": "MONOMIAL_WENDLAND", "polybasis": "LEGENDRE", 'samplerPivot':QS(zs, "UNIFORM"), 'trainSetGenerator':QS(zs, "UNIFORM"), 'errorEstimatorKind':"LOOK_AHEAD_RES", 'errorEstimatorKindMarginal':"LOOK_AHEAD_RECOVER", "SMarginal": 3, "paramsMarginal": {"MMarginal": 2, "radialDirectionalWeightsMarginalAdapt": [1e9, 1e12]}, "greedyTolMarginal": 1e-2, "samplerMarginal":SGS(Ls), "radialDirectionalWeightsMarginal": [4.], "matchingWeight": 1.} for shared, tol in product([1., 0.], [1., 3.]): print("Testing cutoff tolerance {} with shared ratio {}.".format(tol, shared)) solver.cutOffPolesRMinRel = - 1. - tol solver.cutOffPolesRMaxRel = 1. + tol - params['sharedRatio'] = shared + params['matchingShared'] = shared approx = RIGPG([0], solver, mu0 = [z0, L0], approxParameters = params, verbosity = 5) approx.setupApprox("ALL") verb = approx.verbosity approx.verbosity = 0 tspace = np.linspace(Ls[0], Ls[-1], 100) for j, t in enumerate(tspace): plsE = plsEx(t, 0., zs[-1]) pls = approx.getPoles([None, t]) pls[np.abs(np.imag(pls)) > 1e-5] = np.nan if j == 0: polesE = np.empty((len(tspace), len(plsE))) poles = np.empty((len(tspace), len(pls))) polesE[:] = np.nan if len(plsE) > polesE.shape[1]: nanR = np.empty((len(tspace), len(plsE) - polesE.shape[1])) nanR[:] = np.nan polesE = np.hstack((polesE, nanR)) polesE[j, : len(plsE)] = np.real(plsE) poles[j] = np.real(pls) approx.verbosity = verb fighandles += [plt.figure(figsize = (17, 5))] ax1 = fighandles[-1].add_subplot(1, 2, 1) ax2 = fighandles[-1].add_subplot(1, 2, 2) ax1.plot(poles, tspace) ax1.set_ylim(Ls) ax1.set_xlabel('mu_1') ax1.set_ylabel('mu_2') ax1.grid() ax2.plot(polesE, tspace, 'k-.', linewidth = 1) ax2.plot(poles, tspace) for mm in approx.musMarginal: ax2.plot(zs, [mm[0, 0]] * 2, 'k--', linewidth = 1) ax2.set_xlim(zs) ax2.set_ylim(Ls) ax2.set_xlabel('mu_1') ax2.set_ylabel('mu_2') ax2.grid() plt.show() print("\n") diff --git a/examples/9_active_remeshing/active_remeshing.py b/examples/9_active_remeshing/active_remeshing.py index 6c033c4..f247787 100755 --- a/examples/9_active_remeshing/active_remeshing.py +++ b/examples/9_active_remeshing/active_remeshing.py @@ -1,118 +1,118 @@ import numpy as np from pickle import load from matplotlib import pyplot as plt from active_remeshing_engine import ActiveRemeshingEngine from rrompy.reduction_methods import ( RationalInterpolantGreedyPivotedNoMatch as RIGPNM, RationalInterpolantGreedyPivotedGreedyPoleMatch as RIGPG) from rrompy.parameter.parameter_sampling import (QuadratureSampler as QS, SparseGridSampler as SGS) zs, ts = [0., 100.], [0., .5] z0, t0, n = np.mean(zs), np.mean(ts), 150 solver = ActiveRemeshingEngine(z0, t0, n) mus = [[z0, ts[0]], [z0, ts[1]]] for mu in mus: u = solver.solve(mu, return_state = True)[0] Y = solver.applyC(u, mu)[0][0] _ = solver.plot(u, what = "REAL", name = "u(z={}, t={})".format(*mu), is_state = True, figsize = (12, 4)) print("Y(z={}, t={}) = {} (solution norm squared)".format(*mu, Y)) fighandles = [] with open("./active_remeshing_hf_samples.pkl", "rb") as f: zspace, tspace, Yex = load(f) # zspace = np.linspace(zs[0], zs[-1], 198) # tspace = np.linspace(ts[0], ts[-1], 50) # Yex = np.log10([[solver.solve([z, t]) for t in tspace] for z in zspace]) # (from a ~2.5h simulation on one node of the EPFL Helvetios cluster) for match in [0, 1]: params = {"POD": True, "S": 5, "greedyTol": 1e-4, "nTestPoints": 500, "polybasis": "LEGENDRE", "trainSetGenerator": QS(zs, "UNIFORM"), "samplerPivot":QS(zs, "CHEBYSHEV"), "samplerMarginal":SGS(ts), "errorEstimatorKind": "LOOK_AHEAD_OUTPUT"} if match: print("\nTesting output-based matching with weight 1.") params["SMarginal"] = 3 params["maxIterMarginal"] = 25 params["greedyTolMarginal"] = 1e-2 params["matchingWeight"] = 1. - params["sharedRatio"] = .75 + params["matchingShared"] = .75 params["polybasisMarginal"] = "PIECEWISE_LINEAR_UNIFORM" params["errorEstimatorKindMarginal"] = "LOOK_AHEAD_RECOVER" algo = RIGPG else: print("\nTesting matching-free approach.") params["SMarginal"] = 5 algo = RIGPNM approx = algo([0], solver, mu0 = [z0, t0], approxParameters = params, verbosity = 15) approx.setupApprox("ALL" * match) verb, verbTM = approx.verbosity, approx.trainedModel.verbosity approx.verbosity, approx.trainedModel.verbosity = 0, 0 for j, t in enumerate(tspace): out = approx.getApprox(np.pad(zspace.reshape(-1, 1), [(0, 0), (0, 1)], "constant", constant_values = t)) pls = approx.getPoles([None, t]) pls[np.abs(np.imag(pls)) > 1e-5] = np.nan if j == 0: Ys = np.empty((len(zspace), len(tspace))) poles = np.empty((len(tspace), len(pls))) Ys[:, j] = np.log10(out.data) if len(pls) > poles.shape[1]: poles = np.pad(poles, [(0, 0), (0, len(pls) - poles.shape[1])], "constant", constant_values = np.nan) poles[j, : len(pls)] = np.real(pls) approx.verbosity, approx.trainedModel.verbosity = verb, verbTM Ymin, Ymax = min(np.min(Ys), np.min(Yex)), max(np.max(Ys), np.max(Yex)) fighandles += [plt.figure(figsize = (15, 5))] ax1 = fighandles[-1].add_subplot(1, 2, 1) ax2 = fighandles[-1].add_subplot(1, 2, 2) if match: ax1.plot(poles, tspace) else: ax1.plot(poles, tspace, "k.") ax1.set_ylim(ts) ax1.set_xlabel("z") ax1.set_ylabel("t") ax1.grid() if match: ax2.plot(poles, tspace) else: ax2.plot(poles, tspace, "k.") for mm in approx.musMarginal: ax2.plot(zs, [mm[0, 0]] * 2, "k--", linewidth = 1) ax2.set_xlim(zs) ax2.set_ylim(ts) ax2.set_xlabel("z") ax2.set_ylabel("t") ax2.grid() plt.show() print("Approximate poles") fighandles += [plt.figure(figsize = (15, 5))] ax1 = fighandles[-1].add_subplot(1, 2, 1) ax2 = fighandles[-1].add_subplot(1, 2, 2) p = ax1.contourf(np.repeat(zspace.reshape(-1, 1), len(tspace), axis = 1), np.repeat(tspace.reshape(1, -1), len(zspace), axis = 0), Ys, vmin = Ymin, vmax = Ymax, cmap = "gray_r", levels = np.linspace(Ymin, Ymax, 50)) plt.colorbar(p, ax = ax1) ax1.set_xlabel("z") ax1.set_ylabel("t") ax1.grid() p = ax2.contourf(np.repeat(zspace.reshape(-1, 1), len(tspace), axis = 1), np.repeat(tspace.reshape(1, -1), len(zspace), axis = 0), Yex, vmin = Ymin, vmax = Ymax, cmap = "gray_r", levels = np.linspace(Ymin, Ymax, 50)) ax2.set_xlabel("z") ax2.set_ylabel("t") ax2.grid() plt.colorbar(p, ax = ax2) plt.show() print("Approximate and exact output\n") diff --git a/rrompy/hfengines/base/hfengine_base.py b/rrompy/hfengines/base/hfengine_base.py index bce8ffa..ae452c1 100644 --- a/rrompy/hfengines/base/hfengine_base.py +++ b/rrompy/hfengines/base/hfengine_base.py @@ -1,402 +1,427 @@ # Copyright (C) 2018-2020 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 import scipy.sparse as scsp from numbers import Number from collections.abc import Iterable from copy import copy as softcopy from rrompy.utilities.base.decorators import (nonaffine_construct, mu_independent) from rrompy.utilities.base.types import (Np1D, Np2D, List, DictAny, paramVal, paramList, sampList) from rrompy.utilities.numerical import solve as tsolve, dot from rrompy.utilities.expression import expressionEvaluator from rrompy.utilities.exception_manager import RROMPyException, RROMPyAssert from rrompy.sampling.sample_list import sampleList from rrompy.parameter import (checkParameter, checkParameterList, parameterList, parameterMap as pMap) from rrompy.solver.linear_solver import setupSolver from rrompy.utilities.parallel import (poolRank, masterCore, listScatter, matrixGatherv, isend, recv) __all__ = ['HFEngineBase'] class HFEngineBase: """Generic solver for parametric problems.""" def __init__(self, verbosity : int = 10, timestamp : bool = True): self.verbosity = verbosity self.timestamp = timestamp self.setSolver("SPSOLVE", {"use_umfpack" : False}) self.npar = 0 self._C = None self.outputNormMatrix = 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 __dir_base__(self): return [x for x in self.__dir__() if x[:2] != "__"] def __deepcopy__(self, memo): return softcopy(self) @property def npar(self): """Value of npar.""" return self._npar @npar.setter def npar(self, npar): nparOld = self._npar if hasattr(self, "_npar") else -1 if npar != nparOld: self.parameterMap = pMap(1., npar) self._npar = npar @property def spacedim(self): return 1 def checkParameter(self, mu:paramVal) -> paramVal: muP = checkParameter(mu, self.npar) if self.npar == 0: muP.reset((1, 0), muP.dtype) return muP def checkParameterList(self, mu:paramList, check_if_single : bool = False) -> paramList: muL = checkParameterList(mu, self.npar, check_if_single) return muL def mapParameterList(self, mu:paramList, direct : str = "F", idx : List[int] = None) -> paramList: if idx is None: idx = np.arange(self.npar) muMapped = checkParameterList(mu, len(idx)) for j, d in enumerate(idx): muMapped.data[:, j] = expressionEvaluator( self.parameterMap[direct][d], muMapped(j)).flatten() return muMapped def buildEnergyNormForm(self): """ Build sparse matrix (in CSR format) representative of scalar product. """ self.energyNormMatrix = 1. def buildEnergyNormDualForm(self): """ Build sparse matrix (in CSR format) representative of dual scalar product without duality. """ self.energyNormDualMatrix = 1. def innerProduct(self, u:Np2D, v:Np2D, onlyDiag : bool = False, dual : bool = False, is_state : bool = True) -> Np2D: """Scalar product.""" if is_state or self.isCEye: if dual: if not hasattr(self, "energyNormDualMatrix"): self.buildEnergyNormDualForm() energyMat = self.energyNormDualMatrix else: if not hasattr(self, "energyNormMatrix"): self.buildEnergyNormForm() energyMat = self.energyNormMatrix else: energyMat = self.outputNormMatrix if isinstance(u, (parameterList, sampleList)): u = u.data if isinstance(v, (parameterList, sampleList)): v = v.data if onlyDiag: return np.sum(dot(energyMat, u) * v.conj(), axis = 0) return dot(dot(energyMat, u).T, v.conj()).T def norm(self, u:Np2D, dual : bool = False, is_state : bool = True) -> Np1D: return np.abs(self.innerProduct(u, u, onlyDiag = True, dual = dual, is_state = is_state)) ** .5 def baselineA(self): """Return 0 of shape consistent with operator of linear system.""" if (hasattr(self, "As") and isinstance(self.As, Iterable) and self.As[0] is not None): d = self.As[0].shape[0] else: d = self.spacedim return scsp.csr_matrix((np.zeros(0), np.zeros(0), np.zeros(d + 1)), shape = (d, d), dtype = np.complex) def baselineb(self): """Return 0 of shape consistent with RHS of linear system.""" return np.zeros(self.spacedim, dtype = np.complex) @nonaffine_construct @abstractmethod def A(self, mu : paramVal = [], der : List[int] = 0) -> Np2D: """ Assemble terms of operator of linear system and return it (or its derivative) at a given parameter. """ return @nonaffine_construct @abstractmethod def b(self, mu : paramVal = [], der : List[int] = 0) -> Np1D: """ Assemble terms of RHS of linear system and return it (or its derivative) at a given parameter. """ return @mu_independent def C(self, mu:paramVal): """ Value of C. Should be overridden (with something like return self._C(mu) ) if a mu-dependent C is needed. """ if self._C is None: self._C = 1. return self._C def setCQuadratic(self, quad : bool = 1): if isinstance(quad, (str,)): quad = quad.upper().strip().replace(" ","") if quad == "SELFADJOINT": quad = 2 elif quad in ["STANDARD", "GENERAL"]: quad = 1 else: raise RROMPyException(("String keyword for output symmetry " "not recognized.")) if quad: self._is_C_quadratic = int(quad) elif hasattr(self, "_is_C_quadratic"): del self._is_C_quadratic @property def isCEye(self): """ Whether the action of C can be seen as a scalar multiplication. Should be overridden (with return True ) if a mu-dependent scalar C is used. """ return isinstance(self._C, Number) def applyC(self, u:sampList, mu:paramVal): """Apply LHS of linear system.""" if not (hasattr(self, "_is_C_quadratic") and self._is_C_quadratic): return dot(self.C(mu), u) return self._applyCQuadratic(u, mu) def _applyCQuadratic(self, u:sampList, mu:paramVal, v : sampList = None, onlyDiag : bool = False): """Apply quadratic LHS of linear system.""" if not (hasattr(self, "_is_C_quadratic") and self._is_C_quadratic): raise RROMPyException(("Cannot call quadratic output routine if " "output is not quadratic.")) C = self.C(mu) is_C_list = isinstance(C, list) if ((is_C_list and np.any([Ci.ndim != 2 for Ci in C])) or not (is_C_list or C.ndim in [2, 3])): raise RROMPyException(("C array for quadratic output must have 2 " "or 3 dimensions or be list of " "2-dimensional arrays.")) symmetry = v is None # computing u^H * C * u selfadjoint = self._is_C_quadratic == 2 and symmetry if isinstance(u, sampleList): u = u.data if isinstance(v, sampleList): v = v.data while u.ndim < 2: u = np.expand_dims(u, -1) if v is None: v = u while v.ndim < 2: v = np.expand_dims(v, -1) N, M = u.shape[1], v.shape[1] if onlyDiag: if N != M: raise RROMPyException(("Cannot extract diagonal of " "rectangular output.")) Ncol = (N,) elif symmetry: Ncol = (N * (N + 1) // 2,) if selfadjoint else (N ** 2,) else: Ncol = (N, M) for j in range(N): Rv = v[:, j] if onlyDiag else v[:, j :] if selfadjoint else v if is_C_list: cj = np.array([[dot(dot(Ci, u[:, j]), Rv.conj())] for Ci in C]) else: cj = dot(dot(C, u[:, j]), Rv.conj()) if C.ndim == 2: cj = np.expand_dims(cj, 0) if selfadjoint and (onlyDiag or N == 1): cj = np.real(cj) if j == 0: res = np.empty((len(cj),) + Ncol, dtype = cj[0].dtype) if not onlyDiag and symmetry: # map 2D idx to hierarchical 1D idx if self._is_C_quadratic == 2: res[:, j * (j + 1) // 2] = cj[:, 0] iD = np.arange(j + 1, N) * np.arange(j + 4, N + 3) // 2 - j res[:, iD] = 2 * cj[:, 1 :] # exploit symmetry else: res[:, j ** 2 : j * (j + 1)] = cj[:, : j] iD = np.arange(j, N) * np.arange(j + 2, N + 2) - j res[:, iD] = cj[:, j :] else: res[:, j] = cj return res def setSolver(self, solverType:str, solverArgs : DictAny = {}): """Choose solver type and parameters.""" self._solver, self._solverArgs = setupSolver(solverType, solverArgs) def solve(self, mu : paramList = [], RHS : sampList = None, return_state : bool = False) -> sampList: """ Find solution of linear system. Args: mu: parameter value. RHS: RHS of linear system. If None, defaults to that of parametric system. Defaults to None. return_state: whether to return state before multiplication by c. Defaults to False. """ from rrompy.sampling import sampleList, emptySampleList if mu == []: mu = self.mu0 mu = self.checkParameterList(mu) if len(mu) == 0: return emptySampleList() mu = self.checkParameterList(mu) mu_loc, idx, sizes = listScatter(mu, return_sizes = True) mu_loc = self.checkParameterList(mu_loc) req, emptyCores = [], np.where(np.logical_not(sizes))[0] if len(mu_loc) == 0: uL, uT = recv(source = 0, tag = poolRank()) sol = np.empty((uL, 0), dtype = uT) else: applyCglob = (hasattr(self, "_is_C_quadratic") and self._is_C_quadratic) if RHS is None: # build RHSs RHS = sampleList([self.b(m) for m in mu_loc]) else: RHS = sampleList(RHS) if len(RHS) > 1: RHS = sampleList([RHS[i] for i in idx]) mult = 0 if len(RHS) == 1 else 1 RROMPyAssert(mult * (len(mu_loc) - 1) + 1, len(RHS), "Sample size") for j, mj in enumerate(mu_loc): u = tsolve(self.A(mj), RHS[mult * j], self._solver, self._solverArgs) if not (return_state or applyCglob): u = self.applyC(u, mj) if j == 0: sol = np.empty((len(u), len(mu_loc)), dtype = u.dtype) if masterCore(): for dest in emptyCores: req += [isend((len(u), u.dtype), dest = dest, tag = dest)] sol[:, j] = u for r in req: r.wait() sol = matrixGatherv(sol, sizes) if not return_state and applyCglob: sol = self.applyC(sol, mu) return sampleList(sol) def residual(self, mu : paramList = [], u : sampList = None, post_c : bool = True) -> sampList: """ Find residual of linear system for given approximate solution. Args: mu: parameter value. u: numpy complex array with function dofs. If None, set to 0. post_c: whether to post-process using c. Defaults to True. """ from rrompy.sampling import sampleList, emptySampleList if mu == []: mu = self.mu0 mu = self.checkParameterList(mu) if len(mu) == 0: return emptySampleList() mu_loc, idx, sizes = listScatter(mu, return_sizes = True) mu_loc = self.checkParameterList(mu_loc) req, emptyCores = [], np.where(np.logical_not(sizes))[0] if len(mu_loc) == 0: uL, uT = recv(source = 0, tag = poolRank()) res = np.empty((uL, 0), dtype = uT) else: applyCglob = (hasattr(self, "_is_C_quadratic") and self._is_C_quadratic) v = sampleList(np.zeros((self.spacedim, len(mu_loc)))) if u is not None: u = sampleList(u) v = v + sampleList([u[i] for i in idx]) for j, (mj, vj) in enumerate(zip(mu_loc, v)): r = self.b(mj) - dot(self.A(mj), vj) if post_c and not applyCglob: r = self.applyC(r, mj) if j == 0: res = np.empty((len(r), len(mu_loc)), dtype = r.dtype) if masterCore(): for dest in emptyCores: req += [isend((len(r), r.dtype), dest = dest, tag = dest)] res[:, j] = r for r in req: r.wait() res = matrixGatherv(res, sizes) if post_c and applyCglob: res = self.applyC(res, mu) return sampleList(res) cutOffPolesRMax,cutOffPolesRMin = np.inf, - np.inf cutOffPolesRMaxRel, cutOffPolesRMinRel = np.inf, - np.inf cutOffPolesIMax, cutOffPolesIMin = np.inf, - np.inf cutOffPolesIMaxRel, cutOffPolesIMinRel = np.inf, - np.inf cutOffResNormMin = -1 + cutOffResAngleMin, cutOffResAngleMax = -1, np.pi + 1 + @property + def _ignoreResidues(self): + return (self.cutOffResNormMin <= 0. and self.cutOffResAngleMin <= 0. + and self.cutOffResAngleMax >= np.pi) + def flagBadPolesResidues(self, poles:Np1D, residues : Np1D = None, - relative : bool = False) -> Np1D: + relative : bool = False, + projMat : Np2D = None) -> Np1D: """ Flag (numerical) poles/residues which are impossible. Args: poles: poles to be judged. - residues: residues to be judged. + residues: residues norms to be judged. relative: whether relative values should be used for poles. + projMat: matrix for projection of residues. """ poles = np.array(poles).flatten() flag = np.zeros(len(poles), dtype = bool) - if residues is None: - self._ignoreResidues = self.cutOffResNormMin <= 0. - if relative: - RMax, RMin = self.cutOffPolesRMaxRel, self.cutOffPolesRMinRel - IMax, IMin = self.cutOffPolesIMaxRel, self.cutOffPolesIMinRel - else: - RMax, RMin = self.cutOffPolesRMax, self.cutOffPolesRMin - IMax, IMin = self.cutOffPolesIMax, self.cutOffPolesIMin - if not np.isinf(RMax): - flag = np.logical_or(flag, np.real(poles) > RMax) - if not np.isinf(RMin): - flag = np.logical_or(flag, np.real(poles) < RMin) - if not np.isinf(IMax): - flag = np.logical_or(flag, np.imag(poles) > IMax) - if not np.isinf(IMin): - flag = np.logical_or(flag, np.imag(poles) < IMin) + if relative: + RMax, RMin = self.cutOffPolesRMaxRel, self.cutOffPolesRMinRel + IMax, IMin = self.cutOffPolesIMaxRel, self.cutOffPolesIMinRel else: - residues = np.array(residues).reshape(len(poles), -1) + RMax, RMin = self.cutOffPolesRMax, self.cutOffPolesRMin + IMax, IMin = self.cutOffPolesIMax, self.cutOffPolesIMin + if not np.isinf(RMax): + flag = np.logical_or(flag, np.real(poles) > RMax) + if not np.isinf(RMin): + flag = np.logical_or(flag, np.real(poles) < RMin) + if not np.isinf(IMax): + flag = np.logical_or(flag, np.imag(poles) > IMax) + if not np.isinf(IMin): + flag = np.logical_or(flag, np.imag(poles) < IMin) + if residues is not None and not self._ignoreResidues: + residues = np.array(residues).reshape(-1, len(poles)) + if projMat is None: + resNorm = np.linalg.norm(residues, axis = 0) + else: + residues = projMat.dot(residues) + resNorm = self.norm(residues) if self.cutOffResNormMin > 0.: - if residues.shape[1] == self.spacedim: - resEff = self.norm(residues.T) + resNormEff = resNorm / np.max(resNorm) + flag = np.logical_or(flag, resNormEff < self.cutOffResNormMin) + if self.cutOffResAngleMin > 0. or self.cutOffResAngleMax < np.pi: + if projMat is None: + angles = np.real(residues.T.conj().dot(residues)) else: - resEff = np.linalg.norm(residues, axis = 1) - resEff /= np.max(resEff) - flag = np.logical_or(flag, resEff < self.cutOffResNormMin) + angles = np.real(self.innerProduct(residues, residues)) + resNormEff = resNorm + resNormEff[np.isclose(resNormEff, 0.)] = 1. + angles = np.clip((angles / resNormEff).T / resNormEff, -1., 1.) + angles = np.arccos(angles) + badangles = np.logical_or(angles < self.cutOffResAngleMin, + angles > self.cutOffResAngleMax) + badangles[np.arange(len(angles)), np.arange(len(angles))] = 0 + idx = np.zeros(len(angles), dtype = bool) + while np.sum(badangles) > 0: + idxn = np.argmax(np.sum(badangles, axis = 1)) + badangles[idxn], badangles[:, idxn] = 0, 0 + idx[idxn] = 1 + flag = np.logical_or(flag, idx) return flag diff --git a/rrompy/reduction_methods/pivoted/generic_pivoted_approximant.py b/rrompy/reduction_methods/pivoted/generic_pivoted_approximant.py index afe3b99..41a88bb 100644 --- a/rrompy/reduction_methods/pivoted/generic_pivoted_approximant.py +++ b/rrompy/reduction_methods/pivoted/generic_pivoted_approximant.py @@ -1,783 +1,800 @@ # Copyright (C) 2018-2020 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 os import mkdir, remove, rmdir import numpy as np from collections.abc import Iterable from copy import deepcopy as copy from rrompy.reduction_methods.base.generic_approximant import ( GenericApproximant) +from .trained_model.convert_trained_model_pivoted import ( + convertTrainedModelPivoted) from rrompy.utilities.base.data_structures import purgeDict, getNewFilename from rrompy.utilities.poly_fitting.polynomial import polybases as ppb from rrompy.utilities.poly_fitting.radial_basis import polybases as rbpb from rrompy.utilities.poly_fitting.piecewise_linear import sparsekinds as sk from rrompy.utilities.base.types import Np2D, paramList, List, ListAny from rrompy.utilities.base import verbosityManager as vbMng from rrompy.utilities.numerical.degree import reduceDegreeN from rrompy.utilities.exception_manager import RROMPyException, RROMPyWarning from rrompy.parameter import checkParameterList from rrompy.utilities.parallel import poolRank, bcast __all__ = ['GenericPivotedApproximantNoMatch', 'GenericPivotedApproximantPoleMatch'] class GenericPivotedApproximantBase(GenericApproximant): def __init__(self, directionPivot:ListAny, *args, storeAllSamples : bool = False, **kwargs): self._preInit() if len(directionPivot) > 1: raise RROMPyException(("Exactly 1 pivot parameter allowed in pole " "matching.")) from rrompy.parameter.parameter_sampling import (EmptySampler as ES, SparseGridSampler as SG) self._addParametersToList(["radialDirectionalWeightsMarginal"], [1.], ["samplerPivot", "SMarginal", "samplerMarginal"], [ES(), 1, SG([[-1.], [1.]])], toBeExcluded = ["sampler"]) self._directionPivot = directionPivot self.storeAllSamples = storeAllSamples if not hasattr(self, "_output_lvl"): self._output_lvl = [] self._output_lvl += [1 / 2] super().__init__(*args, **kwargs) self._postInit() def setupSampling(self): super().setupSampling(False) def initializeModelData(self, datadict): if "directionPivot" in datadict.keys(): from .trained_model.trained_model_pivoted_data import ( TrainedModelPivotedData) data = TrainedModelPivotedData(datadict["mu0"], datadict["mus"], datadict.pop("projMat"), datadict["scaleFactor"], datadict.pop("parameterMap"), datadict["directionPivot"]) if hasattr(self.HFEngine, "_is_C_quadratic") and not ( hasattr(self, "matchState") and self.matchState): data._is_C_quadratic = self.HFEngine._is_C_quadratic return (data, ["mu0", "scaleFactor", "directionPivot", "mus"]) else: return super().initializeModelData(datadict) @property def npar(self): """Number of parameters.""" if hasattr(self, "_temporaryPivot"): return self.nparPivot return super().npar def checkParameterListPivot(self, mu:paramList, check_if_single : bool = False) -> paramList: return checkParameterList(mu, self.nparPivot, check_if_single) def checkParameterListMarginal(self, mu:paramList, check_if_single : bool = False) -> paramList: return checkParameterList(mu, self.nparMarginal, check_if_single) def mapParameterList(self, *args, **kwargs): if hasattr(self, "_temporaryPivot"): return self.mapParameterListPivot(*args, **kwargs) return super().mapParameterList(*args, **kwargs) def mapParameterListPivot(self, mu:paramList, direct : str = "F", idx : List[int] = None): if idx is None: idx = self.directionPivot else: idx = [self.directionPivot[j] for j in idx] return super().mapParameterList(mu, direct, idx) def mapParameterListMarginal(self, mu:paramList, direct : str = "F", idx : List[int] = None): if idx is None: idx = self.directionMarginal else: idx = [self.directionMarginal[j] for j in idx] return super().mapParameterList(mu, direct, idx) @property def mu0(self): """Value of mu0.""" if hasattr(self, "_temporaryPivot"): return self.checkParameterListPivot(self._mu0(self.directionPivot)) return self._mu0 @mu0.setter def mu0(self, mu0): GenericApproximant.mu0.fset(self, mu0) @property def mus(self): """Value of mus. Its assignment may reset snapshots.""" return self._mus @mus.setter def mus(self, mus): mus = self.checkParameterList(mus) musOld = copy(self.mus) if hasattr(self, '_mus') else None if (musOld is None or len(mus) != len(musOld) or not mus == musOld): self.resetSamples() self._mus = mus @property def musMarginal(self): """Value of musMarginal. Its assignment may reset snapshots.""" return self._musMarginal @musMarginal.setter def musMarginal(self, musMarginal): musMarginal = self.checkParameterListMarginal(musMarginal) if hasattr(self, '_musMarginal'): musMOld = copy(self.musMarginal) else: musMOld = None if (musMOld is None or len(musMarginal) != len(musMOld) or not musMarginal == musMOld): self.resetSamples() self._musMarginal = musMarginal @property def SMarginal(self): """Value of SMarginal.""" return self._SMarginal @SMarginal.setter def SMarginal(self, SMarginal): if SMarginal <= 0: raise RROMPyException("SMarginal must be positive.") if hasattr(self, "_SMarginal") and self._SMarginal is not None: Sold = self.SMarginal else: Sold = -1 self._SMarginal = SMarginal self._approxParameters["SMarginal"] = self.SMarginal if Sold != self.SMarginal: self.resetSamples() @property def radialDirectionalWeightsMarginal(self): """Value of radialDirectionalWeightsMarginal.""" return self._radialDirectionalWeightsMarginal @radialDirectionalWeightsMarginal.setter def radialDirectionalWeightsMarginal(self, radialDirWeightsMarg): if isinstance(radialDirWeightsMarg, Iterable): radialDirWeightsMarg = list(radialDirWeightsMarg) else: radialDirWeightsMarg = [radialDirWeightsMarg] self._radialDirectionalWeightsMarginal = radialDirWeightsMarg self._approxParameters["radialDirectionalWeightsMarginal"] = ( self.radialDirectionalWeightsMarginal) @property def directionPivot(self): """Value of directionPivot. Its assignment may reset snapshots.""" return self._directionPivot @directionPivot.setter def directionPivot(self, directionPivot): if hasattr(self, '_directionPivot'): directionPivotOld = copy(self.directionPivot) else: directionPivotOld = None if (directionPivotOld is None or len(directionPivot) != len(directionPivotOld) or not directionPivot == directionPivotOld): self.resetSamples() self._directionPivot = directionPivot @property def directionMarginal(self): return [x for x in range(self.HFEngine.npar) \ if x not in self.directionPivot] @property def nparPivot(self): return len(self.directionPivot) @property def nparMarginal(self): return self.npar - self.nparPivot @property def muBounds(self): """Value of muBounds.""" return self.samplerPivot.lims @property def muBoundsMarginal(self): """Value of muBoundsMarginal.""" return self.samplerMarginal.lims @property def sampler(self): """Proxy of samplerPivot.""" return self._samplerPivot @property def samplerPivot(self): """Value of samplerPivot.""" return self._samplerPivot @samplerPivot.setter def samplerPivot(self, samplerPivot): if 'generatePoints' not in dir(samplerPivot): raise RROMPyException("Pivot sampler type not recognized.") if hasattr(self, '_samplerPivot') and self._samplerPivot is not None: samplerOld = self.samplerPivot self._samplerPivot = samplerPivot self._approxParameters["samplerPivot"] = self.samplerPivot if not 'samplerOld' in locals() or samplerOld != self.samplerPivot: self.resetSamples() @property def samplerMarginal(self): """Value of samplerMarginal.""" return self._samplerMarginal @samplerMarginal.setter def samplerMarginal(self, samplerMarginal): if 'generatePoints' not in dir(samplerMarginal): raise RROMPyException("Marginal sampler type not recognized.") if (hasattr(self, '_samplerMarginal') and self._samplerMarginal is not None): samplerOld = self.samplerMarginal self._samplerMarginal = samplerMarginal self._approxParameters["samplerMarginal"] = self.samplerMarginal if not 'samplerOld' in locals() or samplerOld != self.samplerMarginal: self.resetSamples() def computeScaleFactor(self): """Compute parameter rescaling factor.""" self.scaleFactorPivot = .5 * np.abs(( self.mapParameterListPivot(self.muBounds[0]) - self.mapParameterListPivot(self.muBounds[1]))[0]) self.scaleFactorMarginal = .5 * np.abs(( self.mapParameterListMarginal(self.muBoundsMarginal[0]) - self.mapParameterListMarginal(self.muBoundsMarginal[1]))[0]) self.scaleFactor = np.empty(self.npar) self.scaleFactor[self.directionPivot] = self.scaleFactorPivot self.scaleFactor[self.directionMarginal] = self.scaleFactorMarginal def _setupTrainedModel(self, pMat:Np2D, pMatUpdate : bool = False, pMatOld : Np2D = None, forceNew : bool = False): if forceNew or self.trainedModel is None: self.trainedModel = self.tModelType() self.trainedModel.verbosity = self.verbosity self.trainedModel.timestamp = self.timestamp datadict = {"mu0": self.mu0, "mus": copy(self.mus), "projMat": pMat, "scaleFactor": self.scaleFactor, "parameterMap": self.HFEngine.parameterMap, "directionPivot": self.directionPivot} self.trainedModel.data = self.initializeModelData(datadict)[0] else: self.trainedModel = self.trainedModel if pMatUpdate: self.trainedModel.data.projMat = np.hstack( (self.trainedModel.data.projMat, pMat)) else: self.trainedModel.data.projMat = copy(pMat) self.trainedModel.data.mus = copy(self.mus) self.trainedModel.data.musMarginal = copy(self.musMarginal) def normApprox(self, mu:paramList) -> float: _PODOld, self._POD = self.POD, 0 result = super().normApprox(mu) self._POD = _PODOld return result @property def storedSamplesFilenames(self) -> List[str]: if not hasattr(self, "_sampleBaseFilename"): return [] return [self._sampleBaseFilename + "{}_{}.pkl" .format(idx + 1, self.name()) for idx in range(len(self.musMarginal))] def purgeStoredSamples(self): if not hasattr(self, "_sampleBaseFilename"): return for file in self.storedSamplesFilenames: remove(file) rmdir(self._sampleBaseFilename[: -8]) def storeSamples(self, idx : int = None): """Store samples to file.""" if not hasattr(self, "_sampleBaseFilename"): filenameBase = None if poolRank() == 0: foldername = getNewFilename(self.name(), "samples") mkdir(foldername) filenameBase = foldername + "/sample_" self._sampleBaseFilename = bcast(filenameBase, force = True) if idx is not None: super().storeSamples(self._sampleBaseFilename + str(idx + 1), False) def loadTrainedModel(self, filename:str): """Load trained reduced model from file.""" super().loadTrainedModel(filename) self._musMarginal = self.trainedModel.data.musMarginal + def setTrainedModel(self, model): + """Deepcopy approximation from trained model.""" + super().setTrainedModel(model) + self.trainedModel = convertTrainedModelPivoted(self.trainedModel, + self.tModelType, self, + True) + self._finalizeMarginalization() + self.trainedModel.data.approxParameters = self.approxParameters + class GenericPivotedApproximantNoMatch(GenericPivotedApproximantBase): """ ROM pivoted approximant (without pole matching) computation for parametric problems (ABSTRACT). Args: HFEngine: HF problem solver. mu0(optional): Default parameter. Defaults to 0. directionPivot(optional): Pivot components. Defaults to [0]. approxParameters(optional): Dictionary containing values for main parameters of approximant. Recognized keys are: - 'POD': kind of snapshots orthogonalization; allowed values include 0, 1/2, and 1; defaults to 1, i.e. POD; - 'scaleFactorDer': scaling factors for derivative computation; defaults to 'AUTO'; - '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; - 'radialDirectionalWeightsMarginal': radial basis weights for marginal interpolant; defaults to 1. Defaults to empty dict. verbosity(optional): Verbosity level. Defaults to 10. Attributes: HFEngine: HF problem solver. mu0: Default parameter. directionPivot: Pivot components. mus: Array of snapshot parameters. musMarginal: Array of marginal snapshot parameters. approxParameters: Dictionary containing values for main parameters of approximant. Recognized keys are in parameterList. parameterListSoft: Recognized keys of soft approximant parameters: - 'POD': kind of snapshots orthogonalization; - 'scaleFactorDer': scaling factors for derivative computation; - 'radialDirectionalWeightsMarginal': radial basis weights for marginal interpolant. parameterListCritical: Recognized keys of critical approximant parameters: - 'S': total number of pivot samples current approximant relies upon; - 'samplerPivot': pivot sample point generator; - 'SMarginal': total number of marginal samples current approximant relies upon; - 'samplerMarginal': marginal sample point generator. verbosity: Verbosity level. POD: Kind of snapshots orthogonalization. scaleFactorDer: Scaling factors for derivative computation. 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. radialDirectionalWeightsMarginal: Radial basis weights for marginal interpolant. muBounds: list of bounds for pivot parameter values. muBoundsMarginal: list of bounds for marginal parameter values. samplingEngine: Sampling engine. uHF: High fidelity solution(s) with parameter(s) lastSolvedHF as sampleList. lastSolvedHF: Parameter(s) corresponding to last computed high fidelity solution(s) as parameterList. uApproxReduced: Reduced approximate solution(s) with parameter(s) lastSolvedApprox as sampleList. lastSolvedApproxReduced: Parameter(s) corresponding to last computed reduced approximate solution(s) as parameterList. uApprox: Approximate solution(s) with parameter(s) lastSolvedApprox as sampleList. lastSolvedApprox: Parameter(s) corresponding to last computed approximate solution(s) as parameterList. """ @property def tModelType(self): from .trained_model.trained_model_pivoted_rational_nomatch import ( TrainedModelPivotedRationalNoMatch) return TrainedModelPivotedRationalNoMatch def _finalizeMarginalization(self): self.trainedModel.setupMarginalInterp( [self.radialDirectionalWeightsMarginal]) self.trainedModel.data.approxParameters = copy(self.approxParameters) def _preliminaryMarginalFinalization(self): pass class GenericPivotedApproximantPoleMatch(GenericPivotedApproximantBase): """ ROM pivoted approximant (with pole matching) computation for parametric problems (ABSTRACT). Args: HFEngine: HF problem solver. mu0(optional): Default parameter. Defaults to 0. directionPivot(optional): Pivot components. Defaults to [0]. approxParameters(optional): Dictionary containing values for main parameters of approximant. Recognized keys are: - 'POD': kind of snapshots orthogonalization; allowed values include 0, 1/2, and 1; defaults to 1, i.e. POD; - 'scaleFactorDer': scaling factors for derivative computation; defaults to 'AUTO'; - 'matchState': whether to match the system state rather than the system output; defaults to False; - 'matchingWeight': weight for pole matching optimization; defaults to 1; - - 'sharedRatio': required ratio of marginal points to share + - 'matchingShared': required ratio of marginal points to share resonance; defaults to 1.; - 'S': total number of pivot samples current approximant relies upon; - 'samplerPivot': pivot sample point generator; - 'SMarginal': total number of marginal samples current approximant relies upon; - 'samplerMarginal': marginal sample point generator; - 'polybasisMarginal': type of polynomial basis for marginal interpolation; allowed values include 'MONOMIAL_*', 'CHEBYSHEV_*', 'LEGENDRE_*', 'NEARESTNEIGHBOR', and 'PIECEWISE_LINEAR_*'; defaults to 'MONOMIAL'; - 'paramsMarginal': dictionary of parameters for marginal interpolation; include: . 'MMarginal': degree of marginal interpolant; defaults to 'AUTO', i.e. maximum allowed; not for 'NEARESTNEIGHBOR' or 'PIECEWISE_LINEAR_*'; . 'nNeighborsMarginal': number of marginal nearest neighbors; defaults to 1; only for 'NEARESTNEIGHBOR'; . 'polydegreetypeMarginal': type of polynomial degree for marginal; defaults to 'TOTAL'; not for 'NEARESTNEIGHBOR' or 'PIECEWISE_LINEAR_*'; - . 'interpRcondMarginal': tolerance for marginal interpolation; + . 'interpTolMarginal': tolerance for marginal interpolation; defaults to None; not for 'NEARESTNEIGHBOR' or 'PIECEWISE_LINEAR_*'; . 'radialDirectionalWeightsMarginalAdapt': bounds for adaptive rescaling of marginal radial basis weights; only for radial basis. - 'radialDirectionalWeightsMarginal': radial basis weights for marginal interpolant; defaults to 1. Defaults to empty dict. verbosity(optional): Verbosity level. Defaults to 10. Attributes: HFEngine: HF problem solver. mu0: Default parameter. directionPivot: Pivot components. mus: Array of snapshot parameters. musMarginal: Array of marginal snapshot parameters. approxParameters: Dictionary containing values for main parameters of approximant. Recognized keys are in parameterList. parameterListSoft: Recognized keys of soft approximant parameters: - 'POD': kind of snapshots orthogonalization; - 'scaleFactorDer': scaling factors for derivative computation; - 'matchState': whether to match the system state rather than the system output; - 'matchingWeight': weight for pole matching optimization; - - 'sharedRatio': required ratio of marginal points to share + - 'matchingShared': required ratio of marginal points to share resonance; - 'polybasisMarginal': type of polynomial basis for marginal interpolation; - 'paramsMarginal': dictionary of parameters for marginal interpolation; include: . 'MMarginal': degree of marginal interpolant; . 'nNeighborsMarginal': number of marginal nearest neighbors; . 'polydegreetypeMarginal': type of polynomial degree for marginal; - . 'interpRcondMarginal': tolerance for marginal interpolation; + . 'interpTolMarginal': tolerance for marginal interpolation; . 'radialDirectionalWeightsMarginalAdapt': bounds for adaptive rescaling of marginal radial basis weights. - 'radialDirectionalWeightsMarginal': radial basis weights for marginal interpolant. parameterListCritical: Recognized keys of critical approximant parameters: - 'S': total number of pivot samples current approximant relies upon; - 'samplerPivot': pivot sample point generator; - 'SMarginal': total number of marginal samples current approximant relies upon; - 'samplerMarginal': marginal sample point generator. verbosity: Verbosity level. POD: Kind of snapshots orthogonalization. scaleFactorDer: Scaling factors for derivative computation. matchState: Whether to match the system state rather than the system output. matchingWeight: Weight for pole matching optimization. - sharedRatio: Required ratio of marginal points to share resonance. + matchingShared: Required ratio of marginal points to share resonance. S: Total number of pivot samples current approximant relies upon. samplerPivot: Pivot sample point generator. SMarginal: Total number of marginal samples current approximant relies upon. samplerMarginal: Marginal sample point generator. polybasisMarginal: Type of polynomial basis for marginal interpolation. paramsMarginal: Dictionary of parameters for marginal interpolation. radialDirectionalWeightsMarginal: Radial basis weights for marginal interpolant. muBounds: list of bounds for pivot parameter values. muBoundsMarginal: list of bounds for marginal parameter values. samplingEngine: Sampling engine. uHF: High fidelity solution(s) with parameter(s) lastSolvedHF as sampleList. lastSolvedHF: Parameter(s) corresponding to last computed high fidelity solution(s) as parameterList. uApproxReduced: Reduced approximate solution(s) with parameter(s) lastSolvedApprox as sampleList. lastSolvedApproxReduced: Parameter(s) corresponding to last computed reduced approximate solution(s) as parameterList. uApprox: Approximate solution(s) with parameter(s) lastSolvedApprox as sampleList. lastSolvedApprox: Parameter(s) corresponding to last computed approximate solution(s) as parameterList. """ def __init__(self, *args, **kwargs): self._preInit() self._addParametersToList(["matchState", "matchingWeight", - "sharedRatio", "polybasisMarginal", + "matchingShared", "polybasisMarginal", "paramsMarginal"], [False, 1., 1., "MONOMIAL", {}]) self.parameterMarginalList = ["MMarginal", "nNeighborsMarginal", "polydegreetypeMarginal", - "interpRcondMarginal", + "interpTolMarginal", "radialDirectionalWeightsMarginalAdapt"] super().__init__(*args, **kwargs) self._postInit() @property def tModelType(self): from .trained_model.trained_model_pivoted_rational_polematch import ( TrainedModelPivotedRationalPoleMatch) return TrainedModelPivotedRationalPoleMatch @property def matchState(self): """Value of matchState.""" return self._matchState @matchState.setter def matchState(self, matchState): self._matchState = matchState self._approxParameters["matchState"] = self.matchState @property def matchingWeight(self): """Value of matchingWeight.""" return self._matchingWeight @matchingWeight.setter def matchingWeight(self, matchingWeight): self._matchingWeight = matchingWeight self._approxParameters["matchingWeight"] = self.matchingWeight @property - def sharedRatio(self): - """Value of sharedRatio.""" - return self._sharedRatio - @sharedRatio.setter - def sharedRatio(self, sharedRatio): - if sharedRatio > 1.: + def matchingShared(self): + """Value of matchingShared.""" + return self._matchingShared + @matchingShared.setter + def matchingShared(self, matchingShared): + if matchingShared > 1.: RROMPyWarning("Shared ratio too large. Clipping to 1.") - sharedRatio = 1. - elif sharedRatio < 0.: + matchingShared = 1. + elif matchingShared < 0.: RROMPyWarning("Shared ratio too small. Clipping to 0.") - sharedRatio = 0. - self._sharedRatio = sharedRatio - self._approxParameters["sharedRatio"] = self.sharedRatio + matchingShared = 0. + self._matchingShared = matchingShared + self._approxParameters["matchingShared"] = self.matchingShared @property def polybasisMarginal(self): """Value of polybasisMarginal.""" return self._polybasisMarginal @polybasisMarginal.setter def polybasisMarginal(self, polybasisMarginal): try: polybasisMarginal = polybasisMarginal.upper().strip().replace(" ", "") if polybasisMarginal not in ppb + rbpb + ["NEARESTNEIGHBOR"] + sk: raise RROMPyException( "Prescribed marginal polybasis not recognized.") self._polybasisMarginal = polybasisMarginal except: RROMPyWarning(("Prescribed marginal polybasis not recognized. " "Overriding to 'MONOMIAL'.")) self._polybasisMarginal = "MONOMIAL" self._approxParameters["polybasisMarginal"] = self.polybasisMarginal @property def paramsMarginal(self): """Value of paramsMarginal.""" return self._paramsMarginal @paramsMarginal.setter def paramsMarginal(self, paramsMarginal): paramsMarginal = purgeDict(paramsMarginal, self.parameterMarginalList, dictname = self.name() + ".paramsMarginal", baselevel = 1) keyList = list(paramsMarginal.keys()) if not hasattr(self, "_paramsMarginal"): self._paramsMarginal = {} if "MMarginal" in keyList: MMarg = paramsMarginal["MMarginal"] elif ("MMarginal" in self.paramsMarginal and not hasattr(self, "_MMarginal_isauto")): MMarg = self.paramsMarginal["MMarginal"] else: MMarg = "AUTO" if isinstance(MMarg, str): MMarg = MMarg.strip().replace(" ","") if "-" not in MMarg: MMarg = MMarg + "-0" self._MMarginal_isauto = True self._MMarginal_shift = int(MMarg.split("-")[-1]) MMarg = 0 if MMarg < 0: raise RROMPyException("MMarginal must be non-negative.") self._paramsMarginal["MMarginal"] = MMarg if "nNeighborsMarginal" in keyList: self._paramsMarginal["nNeighborsMarginal"] = max(1, paramsMarginal["nNeighborsMarginal"]) elif "nNeighborsMarginal" not in self.paramsMarginal: self._paramsMarginal["nNeighborsMarginal"] = 1 if "polydegreetypeMarginal" in keyList: try: polydegtypeM = paramsMarginal["polydegreetypeMarginal"]\ .upper().strip().replace(" ","") if polydegtypeM not in ["TOTAL", "FULL"]: raise RROMPyException(("Prescribed polydegreetypeMarginal " "not recognized.")) self._paramsMarginal["polydegreetypeMarginal"] = polydegtypeM except: RROMPyWarning(("Prescribed polydegreetypeMarginal not " "recognized. Overriding to 'TOTAL'.")) self._paramsMarginal["polydegreetypeMarginal"] = "TOTAL" elif "polydegreetypeMarginal" not in self.paramsMarginal: self._paramsMarginal["polydegreetypeMarginal"] = "TOTAL" - if "interpRcondMarginal" in keyList: - self._paramsMarginal["interpRcondMarginal"] = ( - paramsMarginal["interpRcondMarginal"]) - elif "interpRcondMarginal" not in self.paramsMarginal: - self._paramsMarginal["interpRcondMarginal"] = -1 + if "interpTolMarginal" in keyList: + self._paramsMarginal["interpTolMarginal"] = ( + paramsMarginal["interpTolMarginal"]) + elif "interpTolMarginal" not in self.paramsMarginal: + self._paramsMarginal["interpTolMarginal"] = -1 if "radialDirectionalWeightsMarginalAdapt" in keyList: self._paramsMarginal["radialDirectionalWeightsMarginalAdapt"] = ( paramsMarginal["radialDirectionalWeightsMarginalAdapt"]) elif "radialDirectionalWeightsMarginalAdapt" not in self.paramsMarginal: self._paramsMarginal["radialDirectionalWeightsMarginalAdapt"] = [ -1., -1.] self._approxParameters["paramsMarginal"] = self.paramsMarginal def _setMMarginalAuto(self): if (self.polybasisMarginal not in ppb + rbpb or "MMarginal" not in self.paramsMarginal or "polydegreetypeMarginal" not in self.paramsMarginal): raise RROMPyException(("Cannot set MMarginal if " "polybasisMarginal does not allow it.")) self.paramsMarginal["MMarginal"] = max(0, reduceDegreeN( len(self.musMarginal), len(self.musMarginal), self.nparMarginal, self.paramsMarginal["polydegreetypeMarginal"]) - self._MMarginal_shift) vbMng(self, "MAIN", ("Automatically setting MMarginal to {}.").format( self.paramsMarginal["MMarginal"]), 25) def purgeparamsMarginal(self): self.paramsMarginal = {} paramsMbadkeys = [] if self.polybasisMarginal in ppb + rbpb + sk: paramsMbadkeys += ["nNeighborsMarginal"] if self.polybasisMarginal not in rbpb: paramsMbadkeys += ["radialDirectionalWeightsMarginalAdapt"] if self.polybasisMarginal in ["NEARESTNEIGHBOR"] + sk: paramsMbadkeys += ["MMarginal", "polydegreetypeMarginal", - "interpRcondMarginal"] + "interpTolMarginal"] if hasattr(self, "_MMarginal_isauto"): del self._MMarginal_isauto if hasattr(self, "_MMarginal_shift"): del self._MMarginal_shift for key in paramsMbadkeys: if key in self._paramsMarginal: del self._paramsMarginal[key] self._approxParameters["paramsMarginal"] = self.paramsMarginal def _finalizeMarginalization(self): vbMng(self, "INIT", "Checking shared ratio.", 10) - msg = self.trainedModel.checkSharedRatio(self.sharedRatio) + msg = self.trainedModel.checkShared(self.matchingShared) vbMng(self, "DEL", "Done checking." + msg, 10) if self.polybasisMarginal in rbpb + ["NEARESTNEIGHBOR"]: self.computeScaleFactor() rDWMEff = np.array([w * f for w, f in zip( self.radialDirectionalWeightsMarginal, self.scaleFactorMarginal)]) if self.polybasisMarginal in ppb + rbpb + sk: interpPars = [self.polybasisMarginal] if self.polybasisMarginal in ppb + rbpb: if self.polybasisMarginal in rbpb: interpPars += [rDWMEff] interpPars += [self.verbosity >= 5, self.paramsMarginal["polydegreetypeMarginal"] == "TOTAL"] if self.polybasisMarginal in ppb: interpPars += [{}] else: # if self.polybasisMarginal in rbpb: interpPars += [{"optimizeScalingBounds":self.paramsMarginal[ "radialDirectionalWeightsMarginalAdapt"]}] interpPars += [ - {"rcond":self.paramsMarginal["interpRcondMarginal"]}] + {"rcond":self.paramsMarginal["interpTolMarginal"]}] extraPar = hasattr(self, "_MMarginal_isauto") else: # if self.polybasisMarginal in sk: idxEff = [x for x in range(self.samplerMarginal.npoints) if not hasattr(self.trainedModel, "_idxExcl") or x not in self.trainedModel._idxExcl] extraPar = self.samplerMarginal.depth[idxEff] else: # if self.polybasisMarginal == "NEARESTNEIGHBOR": interpPars = [self.paramsMarginal["nNeighborsMarginal"], rDWMEff] extraPar = None self.trainedModel.setupMarginalInterp(self, interpPars, extraPar) self.trainedModel.data.approxParameters = copy(self.approxParameters) def _preliminaryMarginalFinalization(self): vbMng(self, "INIT", "Compressing and matching poles.", 10) self.trainedModel.initializeFromRational(self.matchingWeight, self.HFEngine, self.matchState) vbMng(self, "DEL", "Done compressing and matching poles.", 10) def _postApplyC(self): if self.POD == 1 and not ( hasattr(self.HFEngine.C, "is_mu_independent") and self.HFEngine.C.is_mu_independent in self._output_lvl): raise RROMPyException(("Cannot apply mu-dependent C to " "orthonormalized samples.")) applyCglob = (hasattr(self.HFEngine, "_is_C_quadratic") and self.HFEngine._is_C_quadratic) vbMng(self, "INIT", "Extracting system output from state.", 35) if applyCglob: pMat, dirM = None, self.trainedModel.data.directionMarginal for muM in self.trainedModel.data.musMarginal: idx = np.where([np.allclose(mu(dirM)[0], muM[0]) for mu in self.trainedModel.data.mus])[0] pMati = self.trainedModel.data.projMat[:, idx] musi = self.trainedModel.data.mus[idx] pMati = self.HFEngine.applyC(pMati, musi) if pMat is None: pMat = np.array(pMati) else: pMat = np.append(pMat, pMati, axis = 1) else: pMat = None for j, mu in enumerate(musi): pMij = self.trainedModel.data.projMat[:, j] pMij = np.expand_dims(self.HFEngine.applyC(pMij, mu), -1) if pMati is None: pMat = np.array(pMij) else: pMat = np.append(pMat, pMij, axis = 1) vbMng(self, "DEL", "Done extracting system output.", 35) self.trainedModel.data.projMat = pMat if hasattr(self.HFEngine, "_is_C_quadratic"): self.trainedModel.data._is_C_quadratic = ( self.HFEngine._is_C_quadratic) + def setTrainedModel(self, model): + """Deepcopy approximation from trained model.""" + super().setTrainedModel(model) + self._preliminaryMarginalFinalization() + self.trainedModel.data.approxParameters = self.approxParameters + @abstractmethod def setupApprox(self, *args, **kwargs) -> int: if self.checkComputedApprox(): return -1 self.purgeparamsMarginal() setupOK = super().setupApprox(*args, **kwargs) if self.matchState: self._postApplyC() return setupOK diff --git a/rrompy/reduction_methods/pivoted/greedy/generic_pivoted_greedy_approximant.py b/rrompy/reduction_methods/pivoted/greedy/generic_pivoted_greedy_approximant.py index 74cdf02..7f00194 100644 --- a/rrompy/reduction_methods/pivoted/greedy/generic_pivoted_greedy_approximant.py +++ b/rrompy/reduction_methods/pivoted/greedy/generic_pivoted_greedy_approximant.py @@ -1,654 +1,654 @@ # Copyright (C) 2018-2020 by the RROMPy authors # # This file is part of RROMPy. # # RROMPy is free software: you can redistribute it and/or modify # it under the terms of the GNU Lesser General Public License as published by # the Free Software Foundation, either version 3 of the License, or # (at your option) any later version. # # RROMPy is distributed in the hope that it will be useful, # but WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # GNU Lesser General Public License for more details. # # You should have received a copy of the GNU Lesser General Public License # along with RROMPy. If not, see . # from abc import abstractmethod from copy import deepcopy as copy import numpy as np from collections.abc import Iterable from matplotlib import pyplot as plt from rrompy.reduction_methods.pivoted.generic_pivoted_approximant import ( GenericPivotedApproximantBase, GenericPivotedApproximantPoleMatch) from rrompy.reduction_methods.pivoted.gather_pivoted_approximant import ( gatherPivotedApproximant) from rrompy.utilities.base.types import (Np1D, Np2D, Tuple, List, paramVal, paramList, ListAny) from rrompy.utilities.base import verbosityManager as vbMng from rrompy.utilities.numerical.point_matching import (pointMatching, buildResiduesForDistance) from rrompy.utilities.numerical.point_distances import chordalMetricAngleMatrix from rrompy.utilities.exception_manager import (RROMPyException, RROMPyAssert, RROMPyWarning) from rrompy.parameter import emptyParameterList from rrompy.utilities.parallel import (masterCore, indicesScatter, arrayGatherv, isend) __all__ = ['GenericPivotedGreedyApproximantPoleMatch'] class GenericPivotedGreedyApproximantBase(GenericPivotedApproximantBase): _allowedEstimatorKindsMarginal = ["LOOK_AHEAD", "LOOK_AHEAD_RECOVER", "NONE"] def __init__(self, *args, **kwargs): self._preInit() self._addParametersToList(["matchingWeightError", "errorEstimatorKindMarginal", "greedyTolMarginal", "maxIterMarginal"], [0., "NONE", 1e-1, 1e2]) super().__init__(*args, **kwargs) self._postInit() @property def scaleFactorDer(self): """Value of scaleFactorDer.""" if self._scaleFactorDer == "NONE": return 1. if self._scaleFactorDer == "AUTO": return self._scaleFactorOldPivot return self._scaleFactorDer @scaleFactorDer.setter def scaleFactorDer(self, scaleFactorDer): if isinstance(scaleFactorDer, (str,)): scaleFactorDer = scaleFactorDer.upper() elif isinstance(scaleFactorDer, Iterable): scaleFactorDer = list(scaleFactorDer) self._scaleFactorDer = scaleFactorDer self._approxParameters["scaleFactorDer"] = self._scaleFactorDer @property def samplerMarginal(self): """Value of samplerMarginal.""" return self._samplerMarginal @samplerMarginal.setter def samplerMarginal(self, samplerMarginal): if 'refine' not in dir(samplerMarginal): raise RROMPyException("Marginal sampler type not recognized.") GenericPivotedApproximantBase.samplerMarginal.fset(self, samplerMarginal) @property def errorEstimatorKindMarginal(self): """Value of errorEstimatorKindMarginal.""" return self._errorEstimatorKindMarginal @errorEstimatorKindMarginal.setter def errorEstimatorKindMarginal(self, errorEstimatorKindMarginal): errorEstimatorKindMarginal = errorEstimatorKindMarginal.upper() if errorEstimatorKindMarginal not in ( self._allowedEstimatorKindsMarginal): RROMPyWarning(("Marginal error estimator kind not recognized. " "Overriding to 'NONE'.")) errorEstimatorKindMarginal = "NONE" self._errorEstimatorKindMarginal = errorEstimatorKindMarginal self._approxParameters["errorEstimatorKindMarginal"] = ( self.errorEstimatorKindMarginal) @property def matchingWeightError(self): """Value of matchingWeightError.""" return self._matchingWeightError @matchingWeightError.setter def matchingWeightError(self, matchingWeightError): self._matchingWeightError = matchingWeightError self._approxParameters["matchingWeightError"] = ( self.matchingWeightError) @property def greedyTolMarginal(self): """Value of greedyTolMarginal.""" return self._greedyTolMarginal @greedyTolMarginal.setter def greedyTolMarginal(self, greedyTolMarginal): if greedyTolMarginal < 0: raise RROMPyException("greedyTolMarginal must be non-negative.") if (hasattr(self, "_greedyTolMarginal") and self.greedyTolMarginal is not None): greedyTolMarginalold = self.greedyTolMarginal else: greedyTolMarginalold = -1 self._greedyTolMarginal = greedyTolMarginal self._approxParameters["greedyTolMarginal"] = self.greedyTolMarginal if greedyTolMarginalold != self.greedyTolMarginal: self.resetSamples() @property def maxIterMarginal(self): """Value of maxIterMarginal.""" return self._maxIterMarginal @maxIterMarginal.setter def maxIterMarginal(self, maxIterMarginal): if maxIterMarginal <= 0: raise RROMPyException("maxIterMarginal must be positive.") if (hasattr(self, "_maxIterMarginal") and self.maxIterMarginal is not None): maxIterMarginalold = self.maxIterMarginal else: maxIterMarginalold = -1 self._maxIterMarginal = maxIterMarginal self._approxParameters["maxIterMarginal"] = self.maxIterMarginal if maxIterMarginalold != self.maxIterMarginal: self.resetSamples() def resetSamples(self): """Reset samples.""" super().resetSamples() if not hasattr(self, "_temporaryPivot"): self._mus = emptyParameterList() self._musMarginal = emptyParameterList() if hasattr(self, "samplerMarginal"): self.samplerMarginal.reset() if hasattr(self, "samplingEngine") and self.samplingEngine is not None: self.samplingEngine.resetHistory() def _getDistanceApp(self, polesEx:Np1D, resEx:Np2D, muTest:paramVal, foci:Tuple[float, float], ground:float) -> float: polesAp = self.trainedModel.interpolateMarginalPoles(muTest)[0] if self.matchingWeightError != 0: resAp = self.trainedModel.interpolateMarginalCoeffs(muTest)[0][ : len(polesAp), :] if (hasattr(self.trainedModel.data, "_is_C_quadratic") and self.trainedModel.data._is_C_quadratic): self.trainedModel._setupQuadMapping() projMapping = self.quad_mapping projMappingReal = self.data._is_C_quadratic == 2 else: projMapping, projMappingReal = None, False resEx = buildResiduesForDistance(resEx, self.trainedModel.data.projMat, 0, projMapping, projMappingReal) resAp = buildResiduesForDistance(resAp, self.trainedModel.data.projMat, 0, projMapping, projMappingReal) else: resAp = None dist = chordalMetricAngleMatrix(polesEx, polesAp, self.matchingWeightError, resEx, resAp, self.HFEngine, False) pmR, pmC = pointMatching(dist) return np.mean(dist[pmR, pmC]) def getErrorEstimatorMarginalLookAhead(self) -> Np1D: if not hasattr(self.trainedModel, "_musMExcl"): err = np.zeros(0) err[:] = np.inf self._musMarginalTestIdxs = np.zeros(0, dtype = int) return err self._musMarginalTestIdxs = np.array(self.trainedModel._idxExcl, dtype = int) idx, sizes = indicesScatter(len(self.trainedModel._musMExcl), return_sizes = True) err = [] if len(idx) > 0: self.verbosity -= 25 self.trainedModel.verbosity -= 25 foci = self.samplerPivot.normalFoci() ground = self.samplerPivot.groundPotential() for j in idx: muTest = self.trainedModel._musMExcl[j] HITest = self.trainedModel._HIsExcl[j] polesEx = HITest.poles idxGood = np.logical_not(np.logical_or(np.isinf(polesEx), np.isnan(polesEx))) polesEx = polesEx[idxGood] if self.matchingWeightError != 0: resEx = HITest.coeffs[np.where(idxGood)[0]] else: resEx = None if len(polesEx) == 0: err += [0.] continue err += [self._getDistanceApp(polesEx, resEx, muTest, foci, ground)] self.verbosity += 25 self.trainedModel.verbosity += 25 return arrayGatherv(np.array(err), sizes) def getErrorEstimatorMarginalNone(self) -> Np1D: nErr = len(self.trainedModel.data.musMarginal) self._musMarginalTestIdxs = np.arange(nErr) return (1. + self.greedyTolMarginal) * np.ones(nErr) def errorEstimatorMarginal(self, return_max : bool = False) -> Np1D: vbMng(self.trainedModel, "INIT", "Evaluating error estimator at mu = {}.".format( self.trainedModel.data.musMarginal), 10) if self.errorEstimatorKindMarginal == "NONE": nErr = len(self.trainedModel.data.musMarginal) self._musMarginalTestIdxs = np.arange(nErr) err = (1. + self.greedyTolMarginal) * np.ones(nErr) else:#if self.errorEstimatorKindMarginal[: 10] == "LOOK_AHEAD": err = self.getErrorEstimatorMarginalLookAhead() vbMng(self.trainedModel, "DEL", "Done evaluating error estimator.", 10) if not return_max: return err idxMaxEst = np.where(err > self.greedyTolMarginal)[0] maxErr = err[idxMaxEst] if self.errorEstimatorKindMarginal == "NONE": maxErr = None return err, idxMaxEst, maxErr def plotEstimatorMarginal(self, est:Np1D, idxMax:List[int], estMax:List[float]): if self.errorEstimatorKindMarginal == "NONE": return if (not (np.any(np.isnan(est)) or np.any(np.isinf(est))) and masterCore() and hasattr(self.trainedModel, "_musMExcl")): fig = plt.figure(figsize = plt.figaspect(1. / self.nparMarginal)) for jpar in range(self.nparMarginal): ax = fig.add_subplot(1, self.nparMarginal, 1 + jpar) musre = np.real(self.trainedModel._musMExcl) if len(idxMax) > 0 and estMax is not None: maxrej = musre[idxMax, jpar] errCP = copy(est) idx = np.delete(np.arange(self.nparMarginal), jpar) while len(musre) > 0: if self.nparMarginal == 1: currIdx = np.arange(len(musre)) else: currIdx = np.where(np.isclose(np.sum( np.abs(musre[:, idx] - musre[0, idx]), 1), 0.))[0] currIdxSorted = currIdx[np.argsort(musre[currIdx, jpar])] ax.semilogy(musre[currIdxSorted, jpar], errCP[currIdxSorted], 'k.-', linewidth = 1) musre = np.delete(musre, currIdx, 0) errCP = np.delete(errCP, currIdx) ax.semilogy(self.musMarginal.re(jpar), (self.greedyTolMarginal,) * len(self.musMarginal), '*m') if len(idxMax) > 0 and estMax is not None: ax.semilogy(maxrej, estMax, 'xr') ax.set_xlim(*list(self.samplerMarginal.lims.re(jpar))) ax.grid() plt.tight_layout() plt.show() def _addMarginalSample(self, mus:paramList): mus = self.checkParameterListMarginal(mus) if len(mus) == 0: return self._nmusOld, nmus = len(self.musMarginal), len(mus) if (hasattr(self, "trainedModel") and self.trainedModel is not None and hasattr(self.trainedModel, "_musMExcl")): self._nmusOld += len(self.trainedModel._musMExcl) vbMng(self, "MAIN", ("Adding marginal sample point{} no. {}{} at {} to training " "set.").format("s" * (nmus > 1), self._nmusOld + 1, "--{}".format(self._nmusOld + nmus) * (nmus > 1), mus), 3) self.musMarginal.append(mus) self.setupApproxPivoted(mus) self._preliminaryMarginalFinalization() del self._nmusOld if (self.errorEstimatorKindMarginal[: 10] == "LOOK_AHEAD" and not self.firstGreedyIterM): ubRange = len(self.trainedModel.data.musMarginal) if hasattr(self.trainedModel, "_idxExcl"): shRange = len(self.trainedModel._musMExcl) else: shRange = 0 testIdxs = list(range(ubRange + shRange - len(mus), ubRange + shRange)) for j in testIdxs[::-1]: self.musMarginal.pop(j - shRange) if hasattr(self.trainedModel, "_idxExcl"): testIdxs = self.trainedModel._idxExcl + testIdxs self._updateTrainedModelMarginalSamples(testIdxs) self._finalizeMarginalization() self._SMarginal = len(self.musMarginal) self._approxParameters["SMarginal"] = self.SMarginal self.trainedModel.data.approxParameters["SMarginal"] = self.SMarginal def greedyNextSampleMarginal(self, muidx:List[int], plotEst : str = "NONE") \ -> Tuple[Np1D, List[int], float, paramVal]: RROMPyAssert(self._mode, message = "Cannot add greedy sample.") muidx = self._musMarginalTestIdxs[muidx] if (self.errorEstimatorKindMarginal[: 10] == "LOOK_AHEAD" and not self.firstGreedyIterM): if not hasattr(self.trainedModel, "_idxExcl"): raise RROMPyException(("Sample index to be added not present " "in trained model.")) testIdxs = copy(self.trainedModel._idxExcl) skippedIdx = 0 for cj, j in enumerate(self.trainedModel._idxExcl): if j in muidx: testIdxs.pop(skippedIdx) self.musMarginal.insert(self.trainedModel._musMExcl[cj], j - skippedIdx) else: skippedIdx += 1 if len(self.trainedModel._idxExcl) < (len(muidx) + len(testIdxs)): raise RROMPyException(("Sample index to be added not present " "in trained model.")) self._updateTrainedModelMarginalSamples(testIdxs) self._SMarginal = len(self.musMarginal) self._approxParameters["SMarginal"] = self.SMarginal self.trainedModel.data.approxParameters["SMarginal"] = ( self.SMarginal) self.firstGreedyIterM = False idxAdded = self.samplerMarginal.refine(muidx)[0] self._addMarginalSample(self.samplerMarginal.points[idxAdded]) errorEstTest, muidx, maxErrorEst = self.errorEstimatorMarginal(True) if plotEst == "ALL": self.plotEstimatorMarginal(errorEstTest, muidx, maxErrorEst) return (errorEstTest, muidx, maxErrorEst, self.samplerMarginal.points[muidx]) def _preliminaryTrainingMarginal(self): """Initialize starting snapshots of solution map.""" RROMPyAssert(self._mode, message = "Cannot start greedy algorithm.") if np.sum(self.samplingEngine.nsamples) > 0: return self.resetSamples() self._addMarginalSample(self.samplerMarginal.generatePoints( self.SMarginal)) def _preSetupApproxPivoted(self, mus:paramList) \ -> Tuple[ListAny, ListAny, ListAny]: self.computeScaleFactor() if self.trainedModel is None: self._setupTrainedModel(np.zeros((0, 0))) self.trainedModel.data.Qs, self.trainedModel.data.Ps = [], [] self.trainedModel.data.Psupp = [] self._trainedModelOld = copy(self.trainedModel) self._scaleFactorOldPivot = copy(self.scaleFactor) self.scaleFactor = self.scaleFactorPivot self._temporaryPivot = 1 self._musLoc = copy(self.mus) idx, sizes = indicesScatter(len(mus), return_sizes = True) emptyCores = np.where(np.logical_not(sizes))[0] self.verbosity -= 10 self.samplingEngine.verbosity -= 10 return idx, sizes, emptyCores def _postSetupApproxPivoted(self, mus:Np2D, pMat:Np2D, Ps:ListAny, Qs:ListAny, sizes:ListAny): self.scaleFactor = self._scaleFactorOldPivot del self._scaleFactorOldPivot, self._temporaryPivot pMat, Ps, Qs, mus, nsamples = gatherPivotedApproximant(pMat, Ps, Qs, mus, sizes, self.polybasis) if len(self._musLoc) > 0: self._mus = self.checkParameterList(self._musLoc) self._mus.append(mus) else: self._mus = self.checkParameterList(mus) self.trainedModel = self._trainedModelOld del self._trainedModelOld padLeft = self.trainedModel.data.projMat.shape[1] suppNew = np.append(0, np.cumsum(nsamples)) self._setupTrainedModel(pMat, padLeft > 0) self.trainedModel.data.Qs += Qs self.trainedModel.data.Ps += Ps self.trainedModel.data.Psupp += list(padLeft + suppNew[: -1]) self.trainedModel.data.approxParameters = copy(self.approxParameters) self.verbosity += 10 self.samplingEngine.verbosity += 10 def _localPivotedResult(self, pMat:Np2D, req:ListAny, emptyCores:ListAny, mus:Np2D) -> Tuple[Np2D, ListAny, Np2D]: pMati = self.samplingEngine.projectionMatrix musi = self.samplingEngine.mus if not hasattr(self, "matchState") or not self.matchState: if self.POD == 1 and not ( hasattr(self.HFEngine.C, "is_mu_independent") and self.HFEngine.C.is_mu_independent in self._output_lvl): raise RROMPyException(("Cannot apply mu-dependent C " "to orthonormalized samples.")) vbMng(self, "INIT", "Extracting system output from state.", 35) if (hasattr(self.HFEngine, "_is_C_quadratic") and self.HFEngine._is_C_quadratic): pMati = self.HFEngine.applyC(pMati, musi) else: pMatiEff = None for j, mu in enumerate(musi): pMij = np.expand_dims(self.HFEngine.applyC(pMati[:, j], mu), -1) if pMatiEff is None: pMatiEff = np.array(pMij) else: pMatiEff = np.append(pMatiEff, pMij, axis = 1) pMati = pMatiEff vbMng(self, "DEL", "Done extracting system output.", 35) if pMat is None: mus = copy(musi.data) pMat = copy(pMati) if masterCore(): for dest in emptyCores: req += [isend((len(pMat), pMat.dtype, mus.dtype), dest = dest, tag = dest)] else: mus = np.vstack((mus, musi.data)) pMat = np.hstack((pMat, pMati)) return pMat, req, mus @abstractmethod def setupApproxPivoted(self, mus:paramList) -> int: if self.checkComputedApproxPivoted(): return -1 RROMPyAssert(self._mode, message = "Cannot setup approximant.") vbMng(self, "INIT", "Setting up pivoted approximant.", 10) self._preSetupApproxPivoted() data = [] pass self._postSetupApproxPivoted(mus, data) vbMng(self, "DEL", "Done setting up pivoted approximant.", 10) return 0 def setupApprox(self, plotEst : str = "NONE") -> int: """Compute greedy snapshots of solution map.""" if self.checkComputedApprox(): return -1 RROMPyAssert(self._mode, message = "Cannot start greedy algorithm.") vbMng(self, "INIT", "Starting computation of snapshots.", 3) max2ErrorEst, self.firstGreedyIterM = np.inf, True self._preliminaryTrainingMarginal() if self.errorEstimatorKindMarginal == "NONE": muidx = [] else:#if self.errorEstimatorKindMarginal[: 10] == "LOOK_AHEAD": muidx = np.arange(len(self.trainedModel.data.musMarginal)) self._musMarginalTestIdxs = np.array(muidx) while self.firstGreedyIterM or (max2ErrorEst > self.greedyTolMarginal and self.samplerMarginal.npoints < self.maxIterMarginal): errorEstTest, muidx, maxErrorEst, mu = \ self.greedyNextSampleMarginal(muidx, plotEst) if maxErrorEst is None: max2ErrorEst = 1. + self.greedyTolMarginal else: if len(maxErrorEst) > 0: max2ErrorEst = np.max(maxErrorEst) else: max2ErrorEst = np.max(errorEstTest) vbMng(self, "MAIN", ("Uniform testing error estimate " "{:.4e}.").format(max2ErrorEst), 3) if plotEst == "LAST": self.plotEstimatorMarginal(errorEstTest, muidx, maxErrorEst) vbMng(self, "DEL", ("Done computing snapshots (final snapshot count: " "{}).").format(len(self.mus)), 3) if (self.errorEstimatorKindMarginal == "LOOK_AHEAD_RECOVER" and hasattr(self.trainedModel, "_idxExcl") and len(self.trainedModel._idxExcl) > 0): vbMng(self, "INIT", "Recovering {} test models.".format( len(self.trainedModel._idxExcl)), 7) for j, mu in zip(self.trainedModel._idxExcl, self.trainedModel._musMExcl): self.musMarginal.insert(mu, j) self._updateTrainedModelMarginalSamples() self._finalizeMarginalization() self._SMarginal = len(self.musMarginal) self._approxParameters["SMarginal"] = self.SMarginal self.trainedModel.data.approxParameters["SMarginal"] = ( self.SMarginal) vbMng(self, "DEL", "Done recovering test models.", 7) return 0 def checkComputedApproxPivoted(self) -> bool: return (super().checkComputedApprox() and len(self.musMarginal) == len(self.trainedModel.data.musMarginal)) class GenericPivotedGreedyApproximantPoleMatch( GenericPivotedGreedyApproximantBase, GenericPivotedApproximantPoleMatch): """ ROM pivoted greedy interpolant computation for parametric problems (with pole matching) (ABSTRACT). Args: HFEngine: HF problem solver. mu0(optional): Default parameter. Defaults to 0. directionPivot(optional): Pivot components. Defaults to [0]. approxParameters(optional): Dictionary containing values for main parameters of approximant. Recognized keys are: - 'POD': kind of snapshots orthogonalization; allowed values include 0, 1/2, and 1; defaults to 1, i.e. POD; - 'scaleFactorDer': scaling factors for derivative computation; defaults to 'AUTO'; - 'matchState': whether to match the system state rather than the system output; defaults to False; - 'matchingWeight': weight for pole matching optimization; defaults to 1; - - 'sharedRatio': required ratio of marginal points to share + - 'matchingShared': required ratio of marginal points to share resonance; defaults to 1.; - 'matchingWeightError': weight for pole matching optimization in error estimation; defaults to 0; - 'S': total number of pivot samples current approximant relies upon; - 'samplerPivot': pivot sample point generator; - 'SMarginal': number of starting marginal samples; - 'samplerMarginal': marginal sample point generator via sparse grid; - 'errorEstimatorKindMarginal': kind of marginal error estimator; available values include 'LOOK_AHEAD', 'LOOK_AHEAD_RECOVER', and 'NONE'; defaults to 'NONE'; - 'polybasisMarginal': type of polynomial basis for marginal interpolation; allowed values include 'MONOMIAL_*', 'CHEBYSHEV_*', 'LEGENDRE_*', 'NEARESTNEIGHBOR', and 'PIECEWISE_LINEAR_*'; defaults to 'MONOMIAL'; - 'paramsMarginal': dictionary of parameters for marginal interpolation; include: . 'MMarginal': degree of marginal interpolant; defaults to 'AUTO', i.e. maximum allowed; not for 'NEARESTNEIGHBOR' or 'PIECEWISE_LINEAR_*'; . 'nNeighborsMarginal': number of marginal nearest neighbors; defaults to 1; only for 'NEARESTNEIGHBOR'; . 'polydegreetypeMarginal': type of polynomial degree for marginal; defaults to 'TOTAL'; not for 'NEARESTNEIGHBOR' or 'PIECEWISE_LINEAR_*'; - . 'interpRcondMarginal': tolerance for marginal interpolation; + . 'interpTolMarginal': tolerance for marginal interpolation; defaults to None; not for 'NEARESTNEIGHBOR' or 'PIECEWISE_LINEAR_*'; . 'radialDirectionalWeightsMarginalAdapt': bounds for adaptive rescaling of marginal radial basis weights; only for radial basis. - 'greedyTolMarginal': uniform error tolerance for marginal greedy algorithm; defaults to 1e-1; - 'maxIterMarginal': maximum number of marginal greedy steps; defaults to 1e2; - 'radialDirectionalWeightsMarginal': radial basis weights for marginal interpolant; defaults to 1. Defaults to empty dict. verbosity(optional): Verbosity level. Defaults to 10. Attributes: HFEngine: HF problem solver. mu0: Default parameter. directionPivot: Pivot components. mus: Array of snapshot parameters. musMarginal: Array of marginal snapshot parameters. approxParameters: Dictionary containing values for main parameters of approximant. Recognized keys are in parameterList. parameterListSoft: Recognized keys of soft approximant parameters: - 'POD': kind of snapshots orthogonalization; - 'scaleFactorDer': scaling factors for derivative computation; - 'matchState': whether to match the system state rather than the system output; - 'matchingWeight': weight for pole matching optimization; - - 'sharedRatio': required ratio of marginal points to share + - 'matchingShared': required ratio of marginal points to share resonance; - 'matchingWeightError': weight for pole matching optimization in error estimation; - 'errorEstimatorKindMarginal': kind of marginal error estimator; - 'polybasisMarginal': type of polynomial basis for marginal interpolation; - 'paramsMarginal': dictionary of parameters for marginal interpolation; include: . 'MMarginal': degree of marginal interpolant; . 'nNeighborsMarginal': number of marginal nearest neighbors; . 'polydegreetypeMarginal': type of polynomial degree for marginal; - . 'interpRcondMarginal': tolerance for marginal interpolation; + . 'interpTolMarginal': tolerance for marginal interpolation; . 'radialDirectionalWeightsMarginalAdapt': bounds for adaptive rescaling of marginal radial basis weights. - 'greedyTolMarginal': uniform error tolerance for marginal greedy algorithm; - 'maxIterMarginal': maximum number of marginal greedy steps; - 'radialDirectionalWeightsMarginal': radial basis weights for marginal interpolant. parameterListCritical: Recognized keys of critical approximant parameters: - 'S': total number of pivot samples current approximant relies upon; - 'samplerPivot': pivot sample point generator; - 'SMarginal': total number of marginal samples current approximant relies upon; - 'samplerMarginal': marginal sample point generator via sparse grid. verbosity: Verbosity level. POD: Kind of snapshots orthogonalization. scaleFactorDer: Scaling factors for derivative computation. matchState: Whether to match the system state rather than the system output. matchingWeight: Weight for pole matching optimization. - sharedRatio: Required ratio of marginal points to share resonance. + matchingShared: Required ratio of marginal points to share resonance. matchingWeightError: Weight for pole matching optimization in error estimation. S: Total number of pivot samples current approximant relies upon. samplerPivot: Pivot sample point generator. SMarginal: Total number of marginal samples current approximant relies upon. samplerMarginal: Marginal sample point generator via sparse grid. errorEstimatorKindMarginal: Kind of marginal error estimator. polybasisMarginal: Type of polynomial basis for marginal interpolation. paramsMarginal: Dictionary of parameters for marginal interpolation. greedyTolMarginal: Uniform error tolerance for marginal greedy algorithm. maxIterMarginal: Maximum number of marginal greedy steps. radialDirectionalWeightsMarginal: Radial basis weights for marginal interpolant. muBounds: list of bounds for pivot parameter values. muBoundsMarginal: list of bounds for marginal parameter values. samplingEngine: Sampling engine. uHF: High fidelity solution(s) with parameter(s) lastSolvedHF as sampleList. lastSolvedHF: Parameter(s) corresponding to last computed high fidelity solution(s) as parameterList. uApproxReduced: Reduced approximate solution(s) with parameter(s) lastSolvedApprox as sampleList. lastSolvedApproxReduced: Parameter(s) corresponding to last computed reduced approximate solution(s) as parameterList. uApprox: Approximate solution(s) with parameter(s) lastSolvedApprox as sampleList. lastSolvedApprox: Parameter(s) corresponding to last computed approximate solution(s) as parameterList. """ def _updateTrainedModelMarginalSamples(self, idx : ListAny = []): self.trainedModel.updateEffectiveSamples(idx, self.matchingWeight, self.HFEngine, False) def setupApprox(self, *args, **kwargs) -> int: if self.checkComputedApprox(): return -1 self.purgeparamsMarginal() _polybasisMarginal = self.polybasisMarginal self._polybasisMarginal = ("PIECEWISE_LINEAR_" + self.samplerMarginal.kind) setupOK = super().setupApprox(*args, **kwargs) self._polybasisMarginal = _polybasisMarginal self._finalizeMarginalization() if self.matchState: self._postApplyC() return setupOK diff --git a/rrompy/reduction_methods/pivoted/greedy/rational_interpolant_greedy_pivoted_greedy.py b/rrompy/reduction_methods/pivoted/greedy/rational_interpolant_greedy_pivoted_greedy.py index d78a5d9..0dc1bc2 100644 --- a/rrompy/reduction_methods/pivoted/greedy/rational_interpolant_greedy_pivoted_greedy.py +++ b/rrompy/reduction_methods/pivoted/greedy/rational_interpolant_greedy_pivoted_greedy.py @@ -1,357 +1,355 @@ #Copyright (C) 2018-2020 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_greedy_approximant import ( GenericPivotedGreedyApproximantBase, GenericPivotedGreedyApproximantPoleMatch) from rrompy.reduction_methods.standard.greedy import RationalInterpolantGreedy from rrompy.reduction_methods.standard.greedy.generic_greedy_approximant \ import pruneSamples from rrompy.reduction_methods.pivoted import ( RationalInterpolantGreedyPivotedPoleMatch) from rrompy.utilities.base.types import Np1D, Tuple, paramVal, paramList from rrompy.utilities.base import verbosityManager as vbMng from rrompy.utilities.exception_manager import RROMPyAssert from rrompy.parameter import emptyParameterList from rrompy.utilities.parallel import poolRank, recv __all__ = ['RationalInterpolantGreedyPivotedGreedyPoleMatch'] class RationalInterpolantGreedyPivotedGreedyBase( GenericPivotedGreedyApproximantBase): @property def sampleBatchSize(self): """Value of sampleBatchSize.""" return 1 @property def sampleBatchIdx(self): """Value of sampleBatchIdx.""" return self.S def greedyNextSample(self, muidx:int, plotEst : str = "NONE")\ -> Tuple[Np1D, int, float, paramVal]: """Compute next greedy snapshot of solution map.""" RROMPyAssert(self._mode, message = "Cannot add greedy sample.") mus = copy(self.muTest[muidx]) self.muTest.pop(muidx) for j, mu in enumerate(mus): vbMng(self, "MAIN", ("Adding sample point no. {} at {} to training " "set.").format(len(self.mus) + 1, mu), 3) self.mus.append(mu) self._S = len(self.mus) self._approxParameters["S"] = self.S if (self.samplingEngine.nsamples <= len(mus) - j - 1 or not np.allclose(mu, self.samplingEngine.mus[j - len(mus)])): self.samplingEngine.nextSample(mu) if self._isLastSampleCollinear(): vbMng(self, "MAIN", ("Collinearity above tolerance detected. Starting " "preemptive greedy loop termination."), 3) self._collinearityFlag = 1 errorEstTest = np.empty(len(self.muTest)) errorEstTest[:] = np.nan return errorEstTest, [-1], np.nan, np.nan errorEstTest, muidx, maxErrorEst = self.errorEstimator(self.muTest, True) if plotEst == "ALL": self.plotEstimator(errorEstTest, muidx, maxErrorEst) return errorEstTest, muidx, maxErrorEst, self.muTest[muidx] def _setSampleBatch(self, maxS:int): return self.S def _preliminaryTraining(self): """Initialize starting snapshots of solution map.""" RROMPyAssert(self._mode, message = "Cannot start greedy algorithm.") if self.samplingEngine.nsamples > 0: return self.resetSamples() self.samplingEngine.scaleFactor = self.scaleFactorDer musPivot = self.trainSetGenerator.generatePoints(self.S) while len(musPivot) > self.S: musPivot.pop() muTestBasePivot = self.samplerPivot.generatePoints(self.nTestPoints, False) idxPop = pruneSamples(self.mapParameterListPivot(muTestBasePivot), self.mapParameterListPivot(musPivot), 1e-10 * self.scaleFactorPivot[0]) muTestBasePivot.pop(idxPop) self._mus = emptyParameterList() self.mus.reset((self.S - 1, self.HFEngine.npar)) self.muTest = emptyParameterList() self.muTest.reset((len(muTestBasePivot) + 1, self.HFEngine.npar)) for k in range(self.S - 1): muk = np.empty_like(self.mus[0]) muk[self.directionPivot] = musPivot[k] muk[self.directionMarginal] = self.muMargLoc self.mus[k] = muk for k in range(len(muTestBasePivot)): muk = np.empty_like(self.muTest[0]) muk[self.directionPivot] = muTestBasePivot[k] muk[self.directionMarginal] = self.muMargLoc self.muTest[k] = muk muk = np.empty_like(self.mus[0]) muk[self.directionPivot] = musPivot[-1] muk[self.directionMarginal] = self.muMargLoc self.muTest[-1] = muk if len(self.mus) > 0: vbMng(self, "MAIN", ("Adding first {} sample point{} at {} to training " "set.").format(self.S - 1, "" + "s" * (self.S > 2), self.mus), 3) self.samplingEngine.iterSample(self.mus) self._S = len(self.mus) self._approxParameters["S"] = self.S self.M, self.N = ("AUTO",) * 2 def setupApproxPivoted(self, mus:paramList) -> int: if self.checkComputedApproxPivoted(): return -1 RROMPyAssert(self._mode, message = "Cannot setup approximant.") vbMng(self, "INIT", "Setting up pivoted approximant.", 10) if not hasattr(self, "_plotEstPivot"): self._plotEstPivot = "NONE" idx, sizes, emptyCores = self._preSetupApproxPivoted(mus) S0 = copy(self.S) pMat, Ps, Qs, req, musA = None, [], [], [], None if len(idx) == 0: vbMng(self, "MAIN", "Idling.", 45) if self.storeAllSamples: self.storeSamples() pL, pT, mT = recv(source = 0, tag = poolRank()) pMat = np.empty((pL, 0), dtype = pT) musA = np.empty((0, self.mu0.shape[1]), dtype = mT) else: for i in idx: self.muMargLoc = mus[i] vbMng(self, "MAIN", "Building marginal model no. {} at " "{}.".format(i + 1, self.muMargLoc), 25) self.samplingEngine.resetHistory() self.trainedModel = None self.verbosity -= 5 self.samplingEngine.verbosity -= 10 RationalInterpolantGreedy.setupApprox(self, self._plotEstPivot) self.verbosity += 5 self.samplingEngine.verbosity += 10 if self.storeAllSamples: self.storeSamples(i + self._nmusOld) pMat, req, musA = self._localPivotedResult(pMat, req, emptyCores, musA) Ps += [copy(self.trainedModel.data.P)] Qs += [copy(self.trainedModel.data.Q)] self._S = S0 del self.muMargLoc for r in req: r.wait() self._postSetupApproxPivoted(musA, pMat, Ps, Qs, sizes) vbMng(self, "DEL", "Done setting up pivoted approximant.", 10) return 0 def setupApprox(self, plotEst : str = "NONE") -> int: if self.checkComputedApprox(): return -1 if '_' not in plotEst: plotEst = plotEst + "_NONE" plotEstM, self._plotEstPivot = plotEst.split("_") val = super().setupApprox(plotEstM) return val class RationalInterpolantGreedyPivotedGreedyPoleMatch( RationalInterpolantGreedyPivotedGreedyBase, GenericPivotedGreedyApproximantPoleMatch, RationalInterpolantGreedyPivotedPoleMatch): """ ROM greedy pivoted greedy 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': kind of snapshots orthogonalization; allowed values include 0, 1/2, and 1; defaults to 1, i.e. POD; - 'scaleFactorDer': scaling factors for derivative computation; defaults to 'AUTO'; - 'matchState': whether to match the system state rather than the system output; defaults to False; - 'matchingWeight': weight for pole matching optimization; defaults to 1; - - 'sharedRatio': required ratio of marginal points to share + - 'matchingShared': required ratio of marginal points to share resonance; defaults to 1.; - 'matchingWeightError': weight for pole matching optimization in error estimation; defaults to 0; - 'S': total number of pivot samples current approximant relies upon; - 'samplerPivot': pivot sample point generator; - 'SMarginal': number of starting marginal samples; - 'samplerMarginal': marginal sample point generator via sparse grid; - 'errorEstimatorKindMarginal': kind of marginal error estimator; available values include 'LOOK_AHEAD' and 'LOOK_AHEAD_RECOVER'; defaults to 'NONE'; - 'polybasis': type of polynomial basis for pivot interpolation; defaults to 'MONOMIAL'; - 'polybasisMarginal': type of polynomial basis for marginal interpolation; allowed values include 'MONOMIAL_*', 'CHEBYSHEV_*', 'LEGENDRE_*', 'NEARESTNEIGHBOR', and 'PIECEWISE_LINEAR_*'; defaults to 'MONOMIAL'; - 'paramsMarginal': dictionary of parameters for marginal interpolation; include: . 'MMarginal': degree of marginal interpolant; defaults to 'AUTO', i.e. maximum allowed; not for 'NEARESTNEIGHBOR' or 'PIECEWISE_LINEAR_*'; . 'nNeighborsMarginal': number of marginal nearest neighbors; defaults to 1; only for 'NEARESTNEIGHBOR'; . 'polydegreetypeMarginal': type of polynomial degree for marginal; defaults to 'TOTAL'; not for 'NEARESTNEIGHBOR' or 'PIECEWISE_LINEAR_*'; - . 'interpRcondMarginal': tolerance for marginal interpolation; + . 'interpTolMarginal': tolerance for marginal interpolation; defaults to None; not for 'NEARESTNEIGHBOR' or 'PIECEWISE_LINEAR_*'; . 'radialDirectionalWeightsMarginalAdapt': bounds for adaptive rescaling of marginal radial basis weights; only for radial basis. - 'greedyTol': uniform error tolerance for greedy algorithm; defaults to 1e-2; - 'collinearityTol': collinearity tolerance for greedy algorithm; defaults to 0.; - 'maxIter': maximum number of greedy steps; defaults to 1e2; - 'nTestPoints': number of test points; defaults to 5e2; - 'trainSetGenerator': training sample points generator; defaults to sampler; - 'errorEstimatorKind': kind of error estimator; available values include 'AFFINE', 'DISCREPANCY', 'LOOK_AHEAD', 'LOOK_AHEAD_RES', and 'NONE'; defaults to 'NONE'; - 'greedyTolMarginal': uniform error tolerance for marginal greedy algorithm; defaults to 1e-1; - 'maxIterMarginal': maximum number of marginal greedy steps; defaults to 1e2; - 'radialDirectionalWeightsMarginal': radial basis weights for marginal interpolant; defaults to 1; - 'functionalSolve': strategy for minimization of denominator functional; allowed values include 'NORM', 'DOMINANT', 'NODAL', 'BARYCENTRIC_NORM', and 'BARYCENTRIC[_AVERAGE]' (check pdf in main folder for meaning); defaults to 'NORM'; - - 'interpRcond': tolerance for pivot interpolation; defaults to - None; - - 'robustTol': tolerance for robust rational denominator - management; defaults to 0. + - 'interpTol': tolerance for pivot interpolation; defaults to None; + - 'QTol': tolerance for robust rational denominator management; + defaults to 0. Defaults to empty dict. verbosity(optional): Verbosity level. Defaults to 10. Attributes: HFEngine: HF problem solver. mu0: Default parameter. directionPivot: Pivot components. mus: Array of snapshot parameters. musMarginal: Array of marginal snapshot parameters. approxParameters: Dictionary containing values for main parameters of approximant. Recognized keys are in parameterList. parameterListSoft: Recognized keys of soft approximant parameters: - 'POD': kind of snapshots orthogonalization; - 'scaleFactorDer': scaling factors for derivative computation; - 'matchState': whether to match the system state rather than the system output; - 'matchingWeight': weight for pole matching optimization; - - 'sharedRatio': required ratio of marginal points to share + - 'matchingShared': required ratio of marginal points to share resonance; - 'matchingWeightError': weight for pole matching optimization in error estimation; - 'errorEstimatorKindMarginal': kind of marginal error estimator; - 'polybasis': type of polynomial basis for pivot interpolation; - 'polybasisMarginal': type of polynomial basis for marginal interpolation; - 'paramsMarginal': dictionary of parameters for marginal interpolation; include: . 'MMarginal': degree of marginal interpolant; . 'nNeighborsMarginal': number of marginal nearest neighbors; . 'polydegreetypeMarginal': type of polynomial degree for marginal; - . 'interpRcondMarginal': tolerance for marginal interpolation; + . 'interpTolMarginal': tolerance for marginal interpolation; . 'radialDirectionalWeightsMarginalAdapt': bounds for adaptive rescaling of marginal radial basis weights. - 'greedyTol': uniform error tolerance for greedy algorithm; - 'collinearityTol': collinearity tolerance for greedy algorithm; - 'maxIter': maximum number of greedy steps; - 'nTestPoints': number of test points; - 'trainSetGenerator': training sample points generator; - 'errorEstimatorKind': kind of error estimator; - 'greedyTolMarginal': uniform error tolerance for marginal greedy algorithm; - 'maxIterMarginal': maximum number of marginal greedy steps; - 'radialDirectionalWeightsMarginal': radial basis weights for marginal interpolant; - 'functionalSolve': strategy for minimization of denominator functional; - - 'interpRcond': tolerance for pivot interpolation; - - 'robustTol': tolerance for robust rational denominator - management. + - 'interpTol': tolerance for pivot interpolation; + - 'QTol': 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 via sparse grid. verbosity: Verbosity level. POD: Kind of snapshots orthogonalization. scaleFactorDer: Scaling factors for derivative computation. matchState: Whether to match the system state rather than the system output. matchingWeight: Weight for pole matching optimization. - sharedRatio: Required ratio of marginal points to share resonance. + matchingShared: Required ratio of marginal points to share resonance. matchingWeightError: Weight for pole matching optimization in error estimation. S: Total number of pivot samples current approximant relies upon. samplerPivot: Pivot sample point generator. SMarginal: Total number of marginal samples current approximant relies upon. samplerMarginal: Marginal sample point generator via sparse grid. errorEstimatorKindMarginal: Kind of marginal error estimator. polybasis: Type of polynomial basis for pivot interpolation. polybasisMarginal: Type of polynomial basis for marginal interpolation. paramsMarginal: Dictionary of parameters for marginal interpolation. greedyTol: uniform error tolerance for greedy algorithm. collinearityTol: Collinearity tolerance for greedy algorithm. maxIter: maximum number of greedy steps. nTestPoints: number of starting training points. trainSetGenerator: training sample points generator. errorEstimatorKind: kind of error estimator. greedyTolMarginal: Uniform error tolerance for marginal greedy algorithm. maxIterMarginal: Maximum number of marginal greedy steps. radialDirectionalWeightsMarginal: Radial basis weights for marginal interpolant. functionalSolve: Strategy for minimization of denominator functional. - interpRcond: Tolerance for pivot interpolation. - robustTol: Tolerance for robust rational denominator management. + interpTol: Tolerance for pivot interpolation. + QTol: Tolerance for robust rational denominator management. muBounds: list of bounds for pivot parameter values. muBoundsMarginal: list of bounds for marginal parameter values. samplingEngine: Sampling engine. uHF: High fidelity solution(s) with parameter(s) lastSolvedHF as sampleList. lastSolvedHF: Parameter(s) corresponding to last computed high fidelity solution(s) as parameterList. uApproxReduced: Reduced approximate solution(s) with parameter(s) lastSolvedApprox as sampleList. lastSolvedApproxReduced: Parameter(s) corresponding to last computed reduced approximate solution(s) as parameterList. uApprox: Approximate solution(s) with parameter(s) lastSolvedApprox as sampleList. lastSolvedApprox: Parameter(s) corresponding to last computed approximate solution(s) as parameterList. """ diff --git a/rrompy/reduction_methods/pivoted/greedy/rational_interpolant_pivoted_greedy.py b/rrompy/reduction_methods/pivoted/greedy/rational_interpolant_pivoted_greedy.py index e33f8b0..a6887a6 100644 --- a/rrompy/reduction_methods/pivoted/greedy/rational_interpolant_pivoted_greedy.py +++ b/rrompy/reduction_methods/pivoted/greedy/rational_interpolant_pivoted_greedy.py @@ -1,287 +1,285 @@ # Copyright (C) 2018-2020 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 from numpy import empty, empty_like from .generic_pivoted_greedy_approximant import ( GenericPivotedGreedyApproximantBase, GenericPivotedGreedyApproximantPoleMatch) from rrompy.reduction_methods.standard import RationalInterpolant from rrompy.reduction_methods.pivoted import ( RationalInterpolantPivotedPoleMatch) from rrompy.utilities.base.types import paramList from rrompy.utilities.base import verbosityManager as vbMng from rrompy.utilities.exception_manager import RROMPyAssert from rrompy.parameter import emptyParameterList from rrompy.utilities.parallel import poolRank, recv __all__ = ['RationalInterpolantPivotedGreedyPoleMatch'] class RationalInterpolantPivotedGreedyBase( GenericPivotedGreedyApproximantBase): def computeSnapshots(self): """Compute snapshots of solution map.""" RROMPyAssert(self._mode, message = "Cannot start snapshot computation.") vbMng(self, "INIT", "Starting computation of snapshots.", 5) self.samplingEngine.scaleFactor = self.scaleFactorDer if not hasattr(self, "musPivot") or len(self.musPivot) != self.S: self.musPivot = self.samplerPivot.generatePoints(self.S) while len(self.musPivot) > self.S: self.musPivot.pop() musLoc = emptyParameterList() musLoc.reset((self.S, self.HFEngine.npar)) self.samplingEngine.resetHistory() for k in range(self.S): muk = empty_like(musLoc[0]) muk[self.directionPivot] = self.musPivot[k] muk[self.directionMarginal] = self.muMargLoc musLoc[k] = muk self.samplingEngine.iterSample(musLoc) vbMng(self, "DEL", "Done computing snapshots.", 5) self._m_selfmus = copy(musLoc) self._mus = self.musPivot self._m_HFEparameterMap = copy(self.HFEngine.parameterMap) self.HFEngine.parameterMap = { "F": [self.HFEngine.parameterMap["F"][self.directionPivot[0]]], "B": [self.HFEngine.parameterMap["B"][self.directionPivot[0]]]} def setupApproxPivoted(self, mus:paramList) -> int: if self.checkComputedApproxPivoted(): return -1 RROMPyAssert(self._mode, message = "Cannot setup approximant.") vbMng(self, "INIT", "Setting up pivoted approximant.", 10) idx, sizes, emptyCores = self._preSetupApproxPivoted(mus) pMat, Ps, Qs, req, musA = None, [], [], [], None if len(idx) == 0: vbMng(self, "MAIN", "Idling.", 45) if self.storeAllSamples: self.storeSamples() pL, pT, mT = recv(source = 0, tag = poolRank()) pMat = empty((pL, 0), dtype = pT) musA = empty((0, self.mu0.shape[1]), dtype = mT) else: for i in idx: self.muMargLoc = mus[i] vbMng(self, "MAIN", "Building marginal model no. {} at " "{}.".format(i + 1, self.muMargLoc), 25) self.samplingEngine.resetHistory() self.trainedModel = None self.verbosity -= 5 self.samplingEngine.verbosity -= 5 RationalInterpolant.setupApprox(self) self.verbosity += 5 self.samplingEngine.verbosity += 5 self._mus = self._m_selfmus self.HFEngine.parameterMap = self._m_HFEparameterMap del self._m_selfmus, self._m_HFEparameterMap if self.storeAllSamples: self.storeSamples(i + self._nmusOld) pMat, req, musA = self._localPivotedResult(pMat, req, emptyCores, musA) Ps += [copy(self.trainedModel.data.P)] Qs += [copy(self.trainedModel.data.Q)] del self.muMargLoc for r in req: r.wait() self._postSetupApproxPivoted(musA, pMat, Ps, Qs, sizes) vbMng(self, "DEL", "Done setting up pivoted approximant.", 10) return 0 class RationalInterpolantPivotedGreedyPoleMatch( RationalInterpolantPivotedGreedyBase, GenericPivotedGreedyApproximantPoleMatch, RationalInterpolantPivotedPoleMatch): """ ROM pivoted greedy 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': kind of snapshots orthogonalization; allowed values include 0, 1/2, and 1; defaults to 1, i.e. POD; - 'scaleFactorDer': scaling factors for derivative computation; defaults to 'AUTO'; - 'matchState': whether to match the system state rather than the system output; defaults to False; - 'matchingWeight': weight for pole matching optimization; defaults to 1; - - 'sharedRatio': required ratio of marginal points to share + - 'matchingShared': required ratio of marginal points to share resonance; defaults to 1.; - 'matchingWeightError': weight for pole matching optimization in error estimation; defaults to 0; - 'S': total number of pivot samples current approximant relies upon; - 'samplerPivot': pivot sample point generator; - 'SMarginal': number of starting marginal samples; - 'samplerMarginal': marginal sample point generator via sparse grid; - 'errorEstimatorKindMarginal': kind of marginal error estimator; available values include 'LOOK_AHEAD' and 'LOOK_AHEAD_RECOVER'; defaults to 'NONE'; - 'polybasis': type of polynomial basis for pivot interpolation; defaults to 'MONOMIAL'; - 'polybasisMarginal': type of polynomial basis for marginal interpolation; allowed values include 'MONOMIAL_*', 'CHEBYSHEV_*', 'LEGENDRE_*', 'NEARESTNEIGHBOR', and 'PIECEWISE_LINEAR_*'; defaults to 'MONOMIAL'; - 'paramsMarginal': dictionary of parameters for marginal interpolation; include: . 'MMarginal': degree of marginal interpolant; defaults to 'AUTO', i.e. maximum allowed; not for 'NEARESTNEIGHBOR' or 'PIECEWISE_LINEAR_*'; . 'nNeighborsMarginal': number of marginal nearest neighbors; defaults to 1; only for 'NEARESTNEIGHBOR'; . 'polydegreetypeMarginal': type of polynomial degree for marginal; defaults to 'TOTAL'; not for 'NEARESTNEIGHBOR' or 'PIECEWISE_LINEAR_*'; - . 'interpRcondMarginal': tolerance for marginal interpolation; + . 'interpTolMarginal': tolerance for marginal interpolation; defaults to None; not for 'NEARESTNEIGHBOR' or 'PIECEWISE_LINEAR_*'; . 'radialDirectionalWeightsMarginalAdapt': bounds for adaptive rescaling of marginal radial basis weights; only for radial basis. - 'M': degree of rational interpolant numerator; defaults to 'AUTO', i.e. maximum allowed; - 'N': degree of rational interpolant denominator; defaults to 'AUTO', i.e. maximum allowed; - 'greedyTolMarginal': uniform error tolerance for marginal greedy algorithm; defaults to 1e-1; - 'maxIterMarginal': maximum number of marginal greedy steps; defaults to 1e2; - 'radialDirectionalWeights': radial basis weights for pivot numerator; defaults to 1; - 'radialDirectionalWeightsAdapt': bounds for adaptive rescaling of radial basis weights; defaults to [-1, -1]; - 'radialDirectionalWeightsMarginal': radial basis weights for marginal interpolant; defaults to 1; - 'functionalSolve': strategy for minimization of denominator functional; allowed values include 'NORM', 'DOMINANT', 'NODAL', 'BARYCENTRIC_NORM', and 'BARYCENTRIC[_AVERAGE]' (check pdf in main folder for meaning); defaults to 'NORM'; - - 'interpRcond': tolerance for pivot interpolation; defaults to - None; - - 'robustTol': tolerance for robust rational denominator - management; defaults to 0. + - 'interpTol': tolerance for pivot interpolation; defaults to None; + - 'QTol': tolerance for robust rational denominator management; + defaults to 0. Defaults to empty dict. 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. approxParameters: Dictionary containing values for main parameters of approximant. Recognized keys are in parameterList. parameterListSoft: Recognized keys of soft approximant parameters: - 'POD': kind of snapshots orthogonalization; - 'scaleFactorDer': scaling factors for derivative computation; - 'matchState': whether to match the system state rather than the system output; - 'matchingWeight': weight for pole matching optimization; - - 'sharedRatio': required ratio of marginal points to share + - 'matchingShared': required ratio of marginal points to share resonance; - 'matchingWeightError': weight for pole matching optimization in error estimation; - 'errorEstimatorKindMarginal': kind of marginal error estimator; - 'polybasis': type of polynomial basis for pivot interpolation; - 'polybasisMarginal': type of polynomial basis for marginal interpolation; - 'paramsMarginal': dictionary of parameters for marginal interpolation; include: . 'MMarginal': degree of marginal interpolant; . 'nNeighborsMarginal': number of marginal nearest neighbors; . 'polydegreetypeMarginal': type of polynomial degree for marginal; - . 'interpRcondMarginal': tolerance for marginal interpolation; + . 'interpTolMarginal': tolerance for marginal interpolation; . 'radialDirectionalWeightsMarginalAdapt': bounds for adaptive rescaling of marginal radial basis weights. - 'M': degree of rational interpolant numerator; - 'N': degree of rational interpolant denominator; - 'greedyTolMarginal': uniform error tolerance for marginal greedy algorithm; - 'maxIterMarginal': maximum number of marginal greedy steps; - 'radialDirectionalWeights': radial basis weights for pivot numerator; - 'radialDirectionalWeightsAdapt': bounds for adaptive rescaling of radial basis weights; - 'radialDirectionalWeightsMarginal': radial basis weights for marginal interpolant; - 'functionalSolve': strategy for minimization of denominator functional; - - 'interpRcond': tolerance for pivot interpolation; - - 'robustTol': tolerance for robust rational denominator - management. + - 'interpTol': tolerance for pivot interpolation; + - 'QTol': 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 via sparse grid. verbosity: Verbosity level. POD: Kind of snapshots orthogonalization. scaleFactorDer: Scaling factors for derivative computation. matchState: Whether to match the system state rather than the system output. matchingWeight: Weight for pole matching optimization. - sharedRatio: Required ratio of marginal points to share resonance. + matchingShared: Required ratio of marginal points to share resonance. matchingWeightError: Weight for pole matching optimization in error estimation. S: Total number of pivot samples current approximant relies upon. samplerPivot: Pivot sample point generator. SMarginal: Total number of marginal samples current approximant relies upon. samplerMarginal: Marginal sample point generator via sparse grid. errorEstimatorKindMarginal: Kind of marginal error estimator. polybasis: Type of polynomial basis for pivot interpolation. polybasisMarginal: Type of polynomial basis for marginal interpolation. paramsMarginal: Dictionary of parameters for marginal interpolation. M: Degree of rational interpolant numerator. N: Degree of rational interpolant denominator. greedyTolMarginal: Uniform error tolerance for marginal greedy algorithm. maxIterMarginal: Maximum number of marginal greedy steps. radialDirectionalWeights: Radial basis weights for pivot numerator. radialDirectionalWeightsAdapt: Bounds for adaptive rescaling of radial basis weights. radialDirectionalWeightsMarginal: Radial basis weights for marginal interpolant. functionalSolve: Strategy for minimization of denominator functional. - interpRcond: Tolerance for pivot interpolation. - robustTol: Tolerance for robust rational denominator management. + interpTol: Tolerance for pivot interpolation. + QTol: Tolerance for robust rational denominator management. muBounds: list of bounds for pivot parameter values. muBoundsMarginal: list of bounds for marginal parameter values. samplingEngine: Sampling engine. uHF: High fidelity solution(s) with parameter(s) lastSolvedHF as sampleList. lastSolvedHF: Parameter(s) corresponding to last computed high fidelity solution(s) as parameterList. uApproxReduced: Reduced approximate solution(s) with parameter(s) lastSolvedApprox as sampleList. lastSolvedApproxReduced: Parameter(s) corresponding to last computed reduced approximate solution(s) as parameterList. uApprox: Approximate solution(s) with parameter(s) lastSolvedApprox as sampleList. lastSolvedApprox: Parameter(s) corresponding to last computed approximate solution(s) as parameterList. """ diff --git a/rrompy/reduction_methods/pivoted/rational_interpolant_greedy_pivoted.py b/rrompy/reduction_methods/pivoted/rational_interpolant_greedy_pivoted.py index bd6b08f..64611cc 100644 --- a/rrompy/reduction_methods/pivoted/rational_interpolant_greedy_pivoted.py +++ b/rrompy/reduction_methods/pivoted/rational_interpolant_greedy_pivoted.py @@ -1,572 +1,569 @@ # Copyright (C) 2018-2020 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 (GenericPivotedApproximantBase, GenericPivotedApproximantNoMatch, GenericPivotedApproximantPoleMatch) from .gather_pivoted_approximant import gatherPivotedApproximant from rrompy.reduction_methods.standard.greedy.rational_interpolant_greedy \ import RationalInterpolantGreedy from rrompy.reduction_methods.standard.greedy.generic_greedy_approximant \ import pruneSamples from rrompy.utilities.base.types import Np1D from rrompy.utilities.base import verbosityManager as vbMng from rrompy.utilities.poly_fitting.polynomial import polyvander as pv from rrompy.utilities.exception_manager import RROMPyException, RROMPyAssert from rrompy.parameter import emptyParameterList, parameterList from rrompy.utilities.parallel import poolRank, indicesScatter, isend, recv __all__ = ['RationalInterpolantGreedyPivotedNoMatch', 'RationalInterpolantGreedyPivotedPoleMatch'] class RationalInterpolantGreedyPivotedBase(GenericPivotedApproximantBase, RationalInterpolantGreedy): def __init__(self, *args, **kwargs): self._preInit() super().__init__(*args, **kwargs) - self._ignoreResidues = self.nparPivot > 1 + if self.nparPivot > 1: self.HFEngine._ignoreResidues = 1 self._postInit() @property def tModelType(self): if hasattr(self, "_temporaryPivot"): return RationalInterpolantGreedy.tModelType.fget(self) return super().tModelType def _polyvanderAuxiliary(self, mus, deg, *args): degEff = [0] * self.npar degEff[self.directionPivot[0]] = deg return pv(mus, degEff, *args) def _marginalizeMiscellanea(self, forward:bool): if forward: self._m_selfmus = copy(self.mus) self._m_HFEparameterMap = copy(self.HFEngine.parameterMap) self._mus = self.checkParameterListPivot( self.mus(self.directionPivot)) self.HFEngine.parameterMap = { "F": [self.HFEngine.parameterMap["F"][self.directionPivot[0]]], "B": [self.HFEngine.parameterMap["B"][self.directionPivot[0]]]} else: self._mus = self._m_selfmus self.HFEngine.parameterMap = self._m_HFEparameterMap del self._m_selfmus, self._m_HFEparameterMap def _marginalizeTrainedModel(self, forward:bool): if forward: del self._temporaryPivot self.trainedModel.data.mu0 = self.mu0 self.trainedModel.data.scaleFactor = [1.] * self.npar self.trainedModel.data.scaleFactor[self.directionPivot[0]] = ( self.scaleFactor[0]) self.trainedModel.data.parameterMap = self.HFEngine.parameterMap self._m_musUniqueCN = copy(self._musUniqueCN) musUniqueCNAux = np.zeros((self.S, self.npar), dtype = self._musUniqueCN.dtype) musUniqueCNAux[:, self.directionPivot[0]] = self._musUniqueCN(0) self._musUniqueCN = self.checkParameterList(musUniqueCNAux) self._m_derIdxs = copy(self._derIdxs) for j in range(len(self._derIdxs)): for l in range(len(self._derIdxs[j])): derjl = self._derIdxs[j][l][0] self._derIdxs[j][l] = [0] * self.npar self._derIdxs[j][l][self.directionPivot[0]] = derjl self.trainedModel.data.Q._dirPivot = self.directionPivot[0] self.trainedModel.data.P._dirPivot = self.directionPivot[0] # tell greedy error estimator that operator / RHS is pivot-affine if hasattr(self.HFEngine.A, "is_affine"): self._A_is_affine = self.HFEngine.A.is_affine else: self._A_is_affine = 0 if hasattr(self.HFEngine.b, "is_affine"): self._b_is_affine = self.HFEngine.b.is_affine else: self._b_is_affine = 0 if self._A_is_affine >= 1 / 2 and self._b_is_affine >= 1 / 2: self._affine_lvl += [1 / 2] else: self._temporaryPivot = 1 self.trainedModel.data.mu0 = self.checkParameterListPivot( self.mu0(self.directionPivot)) self.trainedModel.data.scaleFactor = self.scaleFactor self.trainedModel.data.parameterMap = { "F": [self.HFEngine.parameterMap["F"][self.directionPivot[0]]], "B": [self.HFEngine.parameterMap["B"][self.directionPivot[0]]]} self._musUniqueCN = copy(self._m_musUniqueCN) self._derIdxs = copy(self._m_derIdxs) del self._m_musUniqueCN, self._m_derIdxs del self.trainedModel.data.Q._dirPivot del self.trainedModel.data.P._dirPivot if self._A_is_affine >= 1 / 2 and self._b_is_affine >= 1 / 2: self._affine_lvl.pop() del self._A_is_affine, self._b_is_affine self.trainedModel.data.npar = self.npar def errorEstimator(self, mus:Np1D, return_max : bool = False) -> Np1D: """Standard residual-based error estimator.""" setupOK = self.setupApproxLocal() if setupOK > 0: err = np.empty(len(mus)) err[:] = np.nan if not return_max: return err return err, [- setupOK], np.nan self._marginalizeTrainedModel(True) errRes = super().errorEstimator(mus, return_max) self._marginalizeTrainedModel(False) return errRes def _preliminaryTraining(self): """Initialize starting snapshots of solution map.""" RROMPyAssert(self._mode, message = "Cannot start greedy algorithm.") self._S = self._setSampleBatch(self.S) self.resetSamples() self.samplingEngine.scaleFactor = self.scaleFactorDer musPivot = self.trainSetGenerator.generatePoints(self.S) while len(musPivot) > self.S: musPivot.pop() muTestPivot = self.samplerPivot.generatePoints(self.nTestPoints, False) idxPop = pruneSamples(self.mapParameterListPivot(muTestPivot), self.mapParameterListPivot(musPivot), 1e-10 * self.scaleFactorPivot[0]) self._mus = emptyParameterList() self.mus.reset((self.S, self.npar + len(self.musMargLoc))) muTestBase = emptyParameterList() muTestBase.reset((len(muTestPivot), self.npar + len(self.musMargLoc))) for k in range(self.S): muk = np.empty_like(self.mus[0]) muk[self.directionPivot] = musPivot[k] muk[self.directionMarginal] = self.musMargLoc self.mus[k] = muk for k in range(len(muTestPivot)): muk = np.empty_like(muTestBase[0]) muk[self.directionPivot] = muTestPivot[k] muk[self.directionMarginal] = self.musMargLoc muTestBase[k] = muk muTestBase.pop(idxPop) muLast = copy(self.mus[-1]) self.mus.pop() if len(self.mus) > 0: vbMng(self, "MAIN", ("Adding first {} sample point{} at {} to training " "set.").format(self.S - 1, "" + "s" * (self.S > 2), self.mus), 3) self.samplingEngine.iterSample(self.mus) self._S = len(self.mus) self._approxParameters["S"] = self.S self.muTest = parameterList(muTestBase) self.muTest.append(muLast) self.M, self.N = ("AUTO",) * 2 def setupApproxLocal(self) -> int: """Compute rational interpolant.""" self._marginalizeMiscellanea(True) setupOK = super().setupApproxLocal() self._marginalizeMiscellanea(False) return setupOK def setupApprox(self, *args, **kwargs) -> int: """Compute rational interpolant.""" if self.checkComputedApprox(): return -1 RROMPyAssert(self._mode, message = "Cannot setup approximant.") vbMng(self, "INIT", "Setting up {}.". format(self.name()), 5) self.computeScaleFactor() self._musMarginal = self.samplerMarginal.generatePoints(self.SMarginal) while len(self.musMarginal) > self.SMarginal: self.musMarginal.pop() S0 = copy(self.S) idx, sizes = indicesScatter(len(self.musMarginal), return_sizes = True) pMat, Ps, Qs, mus = None, [], [], None req, emptyCores = [], np.where(np.logical_not(sizes))[0] if len(idx) == 0: vbMng(self, "MAIN", "Idling.", 25) if self.storeAllSamples: self.storeSamples() pL, pT, mT = recv(source = 0, tag = poolRank()) pMat = np.empty((pL, 0), dtype = pT) mus = np.empty((0, self.mu0.shape[1]), dtype = mT) else: _scaleFactorOldPivot = copy(self.scaleFactor) self.scaleFactor = self.scaleFactorPivot self._temporaryPivot = 1 for i in idx: self.musMargLoc = self.musMarginal[i] vbMng(self, "MAIN", "Building marginal model no. {} at {}.".format(i + 1, self.musMargLoc), 5) self.samplingEngine.resetHistory() self.trainedModel = None self.verbosity -= 5 self.samplingEngine.verbosity -= 10 RationalInterpolantGreedy.setupApprox(self, *args, **kwargs) self.verbosity += 5 self.samplingEngine.verbosity += 10 if self.storeAllSamples: self.storeSamples(i) musi = self.samplingEngine.mus pMati = self.samplingEngine.projectionMatrix if not hasattr(self, "matchState") or not self.matchState: if self.POD == 1 and not ( hasattr(self.HFEngine.C, "is_mu_independent") and self.HFEngine.C.is_mu_independent in self._output_lvl): raise RROMPyException(("Cannot apply mu-dependent C " "to orthonormalized samples.")) vbMng(self, "INIT", "Extracting system output from state.", 35) if (hasattr(self.HFEngine, "_is_C_quadratic") and self.HFEngine._is_C_quadratic): pMati = self.HFEngine.applyC(pMati, musi) else: pMatiEff = None for j, mu in enumerate(musi): pMij = np.expand_dims(self.HFEngine.applyC( pMati[:, j], mu), -1) if pMatiEff is None: pMatiEff = np.array(pMij) else: pMatiEff = np.append(pMatiEff, pMij, axis = 1) pMati = pMatiEff vbMng(self, "DEL", "Done extracting system output.", 35) if pMat is None: mus = copy(musi.data) pMat = copy(pMati) if i == 0: for dest in emptyCores: req += [isend((len(pMat), pMat.dtype, mus.dtype), dest = dest, tag = dest)] else: mus = np.vstack((mus, musi.data)) pMat = np.hstack((pMat, pMati)) Ps += [copy(self.trainedModel.data.P)] Qs += [copy(self.trainedModel.data.Q)] self._S = S0 del self._temporaryPivot, self.musMargLoc self.scaleFactor = _scaleFactorOldPivot for r in req: r.wait() pMat, Ps, Qs, mus, nsamples = gatherPivotedApproximant(pMat, Ps, Qs, mus, sizes, self.polybasis) self._mus = self.checkParameterList(mus) Psupp = np.append(0, np.cumsum(nsamples)) self._setupTrainedModel(pMat, forceNew = True) self.trainedModel.data.Qs, self.trainedModel.data.Ps = Qs, Ps self.trainedModel.data.Psupp = list(Psupp[: -1]) self._preliminaryMarginalFinalization() self._finalizeMarginalization() vbMng(self, "DEL", "Done setting up approximant.", 5) return 0 class RationalInterpolantGreedyPivotedNoMatch( RationalInterpolantGreedyPivotedBase, GenericPivotedApproximantNoMatch): """ ROM pivoted rational interpolant (without pole matching) 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': kind of snapshots orthogonalization; allowed values include 0, 1/2, and 1; defaults to 1, i.e. POD; - 'scaleFactorDer': scaling factors for derivative computation; defaults to 'AUTO'; - '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; - 'polybasis': type of polynomial basis for pivot interpolation; defaults to 'MONOMIAL'; - 'greedyTol': uniform error tolerance for greedy algorithm; defaults to 1e-2; - 'collinearityTol': collinearity tolerance for greedy algorithm; defaults to 0.; - 'maxIter': maximum number of greedy steps; defaults to 1e2; - 'nTestPoints': number of test points; defaults to 5e2; - 'trainSetGenerator': training sample points generator; defaults to sampler; - 'errorEstimatorKind': kind of error estimator; available values include 'AFFINE', 'DISCREPANCY', 'LOOK_AHEAD', 'LOOK_AHEAD_RES', and 'NONE'; defaults to 'NONE'; - 'radialDirectionalWeightsMarginal': radial basis weights for marginal interpolant; defaults to 1; - 'functionalSolve': strategy for minimization of denominator functional; allowed values include 'NORM', 'DOMINANT', 'NODAL', 'BARYCENTRIC_NORM', and 'BARYCENTRIC[_AVERAGE]' (check pdf in main folder for meaning); defaults to 'NORM'; - - 'interpRcond': tolerance for pivot interpolation; defaults to + - 'interpTol': tolerance for pivot interpolation; defaults to None; - - 'robustTol': tolerance for robust rational denominator - management; defaults to 0. + - 'QTol': tolerance for robust rational denominator management; + defaults to 0. Defaults to empty dict. verbosity(optional): Verbosity level. Defaults to 10. Attributes: HFEngine: HF problem solver. mu0: Default parameter. directionPivot: Pivot components. mus: Array of snapshot parameters. musMarginal: Array of marginal snapshot parameters. approxParameters: Dictionary containing values for main parameters of approximant. Recognized keys are in parameterList. parameterListSoft: Recognized keys of soft approximant parameters: - 'POD': kind of snapshots orthogonalization; - 'scaleFactorDer': scaling factors for derivative computation; - 'polybasis': type of polynomial basis for pivot interpolation; - 'greedyTol': uniform error tolerance for greedy algorithm; - 'collinearityTol': collinearity tolerance for greedy algorithm; - 'maxIter': maximum number of greedy steps; - 'nTestPoints': number of test points; - 'trainSetGenerator': training sample points generator; - 'errorEstimatorKind': kind of error estimator; - 'radialDirectionalWeightsMarginal': radial basis weights for marginal interpolant; - 'functionalSolve': strategy for minimization of denominator functional; - - 'interpRcond': tolerance for pivot interpolation; - - 'robustTol': tolerance for robust rational denominator - management. + - 'interpTol': tolerance for pivot interpolation; + - 'QTol': 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. verbosity: Verbosity level. POD: Kind of snapshots orthogonalization. scaleFactorDer: Scaling factors for derivative computation. 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. polybasis: Type of polynomial basis for pivot interpolation. greedyTol: uniform error tolerance for greedy algorithm. collinearityTol: Collinearity tolerance for greedy algorithm. maxIter: maximum number of greedy steps. nTestPoints: number of starting training points. trainSetGenerator: training sample points generator. errorEstimatorKind: kind of error estimator. radialDirectionalWeightsMarginal: Radial basis weights for marginal interpolant. functionalSolve: Strategy for minimization of denominator functional. - interpRcond: Tolerance for pivot interpolation. - robustTol: Tolerance for robust rational denominator management. + interpTol: Tolerance for pivot interpolation. + QTol: Tolerance for robust rational denominator management. muBounds: list of bounds for pivot parameter values. muBoundsMarginal: list of bounds for marginal parameter values. samplingEngine: Sampling engine. uHF: High fidelity solution(s) with parameter(s) lastSolvedHF as sampleList. lastSolvedHF: Parameter(s) corresponding to last computed high fidelity solution(s) as parameterList. uApproxReduced: Reduced approximate solution(s) with parameter(s) lastSolvedApprox as sampleList. lastSolvedApproxReduced: Parameter(s) corresponding to last computed reduced approximate solution(s) as parameterList. uApprox: Approximate solution(s) with parameter(s) lastSolvedApprox as sampleList. lastSolvedApprox: Parameter(s) corresponding to last computed approximate solution(s) as parameterList. Q: Numpy 1D vector containing complex coefficients of approximant denominator. P: Numpy 2D vector whose columns are FE dofs of coefficients of approximant numerator. """ class RationalInterpolantGreedyPivotedPoleMatch( RationalInterpolantGreedyPivotedBase, GenericPivotedApproximantPoleMatch): """ ROM pivoted rational interpolant (with pole matching) 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': kind of snapshots orthogonalization; allowed values include 0, 1/2, and 1; defaults to 1, i.e. POD; - 'scaleFactorDer': scaling factors for derivative computation; defaults to 'AUTO'; - 'matchState': whether to match the system state rather than the system output; defaults to False; - 'matchingWeight': weight for pole matching optimization; defaults to 1; - - 'sharedRatio': required ratio of marginal points to share + - 'matchingShared': required ratio of marginal points to share resonance; defaults to 1.; - '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; - 'polybasis': type of polynomial basis for pivot interpolation; defaults to 'MONOMIAL'; - 'polybasisMarginal': type of polynomial basis for marginal interpolation; allowed values include 'MONOMIAL_*', 'CHEBYSHEV_*', 'LEGENDRE_*', 'NEARESTNEIGHBOR', and 'PIECEWISE_LINEAR_*'; defaults to 'MONOMIAL'; - 'paramsMarginal': dictionary of parameters for marginal interpolation; include: . 'MMarginal': degree of marginal interpolant; defaults to 'AUTO', i.e. maximum allowed; not for 'NEARESTNEIGHBOR' or 'PIECEWISE_LINEAR_*'; . 'nNeighborsMarginal': number of marginal nearest neighbors; defaults to 1; only for 'NEARESTNEIGHBOR'; . 'polydegreetypeMarginal': type of polynomial degree for marginal; defaults to 'TOTAL'; not for 'NEARESTNEIGHBOR' or 'PIECEWISE_LINEAR_*'; - . 'interpRcondMarginal': tolerance for marginal interpolation; + . 'interpTolMarginal': tolerance for marginal interpolation; defaults to None; not for 'NEARESTNEIGHBOR' or 'PIECEWISE_LINEAR_*'; . 'radialDirectionalWeightsMarginalAdapt': bounds for adaptive rescaling of marginal radial basis weights; only for radial basis. - 'greedyTol': uniform error tolerance for greedy algorithm; defaults to 1e-2; - 'collinearityTol': collinearity tolerance for greedy algorithm; defaults to 0.; - 'maxIter': maximum number of greedy steps; defaults to 1e2; - 'nTestPoints': number of test points; defaults to 5e2; - 'trainSetGenerator': training sample points generator; defaults to sampler; - 'errorEstimatorKind': kind of error estimator; available values include 'AFFINE', 'DISCREPANCY', 'LOOK_AHEAD', 'LOOK_AHEAD_RES', and 'NONE'; defaults to 'NONE'; - 'radialDirectionalWeightsMarginal': radial basis weights for marginal interpolant; defaults to 1; - 'functionalSolve': strategy for minimization of denominator functional; allowed values include 'NORM', 'DOMINANT', 'NODAL', 'BARYCENTRIC_NORM', and 'BARYCENTRIC[_AVERAGE]' (check pdf in main folder for meaning); defaults to 'NORM'; - - 'interpRcond': tolerance for pivot interpolation; defaults to - None; - - 'robustTol': tolerance for robust rational denominator - management; defaults to 0. + - 'interpTol': tolerance for pivot interpolation; defaults to None; + - 'QTol': tolerance for robust rational denominator management; + defaults to 0. Defaults to empty dict. verbosity(optional): Verbosity level. Defaults to 10. Attributes: HFEngine: HF problem solver. mu0: Default parameter. directionPivot: Pivot components. mus: Array of snapshot parameters. musMarginal: Array of marginal snapshot parameters. approxParameters: Dictionary containing values for main parameters of approximant. Recognized keys are in parameterList. parameterListSoft: Recognized keys of soft approximant parameters: - 'POD': kind of snapshots orthogonalization; - 'scaleFactorDer': scaling factors for derivative computation; - 'matchState': whether to match the system state rather than the system output; - 'matchingWeight': weight for pole matching optimization; - - 'sharedRatio': required ratio of marginal points to share + - 'matchingShared': required ratio of marginal points to share resonance; - 'polybasis': type of polynomial basis for pivot interpolation; - 'polybasisMarginal': type of polynomial basis for marginal interpolation; - 'paramsMarginal': dictionary of parameters for marginal interpolation; include: . 'MMarginal': degree of marginal interpolant; . 'nNeighborsMarginal': number of marginal nearest neighbors; . 'polydegreetypeMarginal': type of polynomial degree for marginal; - . 'interpRcondMarginal': tolerance for marginal interpolation; + . 'interpTolMarginal': tolerance for marginal interpolation; . 'radialDirectionalWeightsMarginalAdapt': bounds for adaptive rescaling of marginal radial basis weights. - 'greedyTol': uniform error tolerance for greedy algorithm; - 'collinearityTol': collinearity tolerance for greedy algorithm; - 'maxIter': maximum number of greedy steps; - 'nTestPoints': number of test points; - 'trainSetGenerator': training sample points generator; - 'errorEstimatorKind': kind of error estimator; - 'radialDirectionalWeightsMarginal': radial basis weights for marginal interpolant; - 'functionalSolve': strategy for minimization of denominator functional; - - 'interpRcond': tolerance for pivot interpolation; - - 'robustTol': tolerance for robust rational denominator - management. + - 'interpTol': tolerance for pivot interpolation; + - 'QTol': 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. verbosity: Verbosity level. POD: Kind of snapshots orthogonalization. scaleFactorDer: Scaling factors for derivative computation. matchState: Whether to match the system state rather than the system output. matchingWeight: Weight for pole matching optimization. - sharedRatio: Required ratio of marginal points to share resonance. + matchingShared: Required ratio of marginal points to share resonance. 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. polybasis: Type of polynomial basis for pivot interpolation. polybasisMarginal: Type of polynomial basis for marginal interpolation. paramsMarginal: Dictionary of parameters for marginal interpolation. greedyTol: uniform error tolerance for greedy algorithm. collinearityTol: Collinearity tolerance for greedy algorithm. maxIter: maximum number of greedy steps. nTestPoints: number of starting training points. trainSetGenerator: training sample points generator. errorEstimatorKind: kind of error estimator. radialDirectionalWeightsMarginal: Radial basis weights for marginal interpolant. functionalSolve: Strategy for minimization of denominator functional. - interpRcond: Tolerance for pivot interpolation. - robustTol: Tolerance for robust rational denominator management. + interpTol: Tolerance for pivot interpolation. + QTol: Tolerance for robust rational denominator management. muBounds: list of bounds for pivot parameter values. muBoundsMarginal: list of bounds for marginal parameter values. samplingEngine: Sampling engine. uHF: High fidelity solution(s) with parameter(s) lastSolvedHF as sampleList. lastSolvedHF: Parameter(s) corresponding to last computed high fidelity solution(s) as parameterList. uApproxReduced: Reduced approximate solution(s) with parameter(s) lastSolvedApprox as sampleList. lastSolvedApproxReduced: Parameter(s) corresponding to last computed reduced approximate solution(s) as parameterList. uApprox: Approximate solution(s) with parameter(s) lastSolvedApprox as sampleList. lastSolvedApprox: Parameter(s) corresponding to last computed approximate solution(s) as parameterList. 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, *args, **kwargs) -> int: if self.checkComputedApprox(): return -1 self.purgeparamsMarginal() setupOK = super().setupApprox(*args, **kwargs) if self.matchState: self._postApplyC() return setupOK diff --git a/rrompy/reduction_methods/pivoted/rational_interpolant_pivoted.py b/rrompy/reduction_methods/pivoted/rational_interpolant_pivoted.py index 95086fd..3142c71 100644 --- a/rrompy/reduction_methods/pivoted/rational_interpolant_pivoted.py +++ b/rrompy/reduction_methods/pivoted/rational_interpolant_pivoted.py @@ -1,487 +1,484 @@ # Copyright (C) 2018-2020 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 collections.abc import Iterable from copy import deepcopy as copy from .generic_pivoted_approximant import (GenericPivotedApproximantBase, GenericPivotedApproximantNoMatch, GenericPivotedApproximantPoleMatch) from .gather_pivoted_approximant import gatherPivotedApproximant from rrompy.reduction_methods.standard.rational_interpolant import ( RationalInterpolant) from rrompy.utilities.base import verbosityManager as vbMng from rrompy.utilities.numerical.hash_derivative import nextDerivativeIndices from rrompy.utilities.exception_manager import (RROMPyException, RROMPyAssert, RROMPyWarning) from rrompy.parameter import emptyParameterList from rrompy.utilities.parallel import poolRank, indicesScatter, isend, recv __all__ = ['RationalInterpolantPivotedNoMatch', 'RationalInterpolantPivotedPoleMatch'] class RationalInterpolantPivotedBase(GenericPivotedApproximantBase, RationalInterpolant): def __init__(self, *args, **kwargs): self._preInit() self._addParametersToList(toBeExcluded = ["polydegreetype"]) super().__init__(*args, **kwargs) - self._ignoreResidues = self.nparPivot > 1 + if self.nparPivot > 1: self.HFEngine._ignoreResidues = 1 self._postInit() @property def scaleFactorDer(self): """Value of scaleFactorDer.""" if self._scaleFactorDer == "NONE": return 1. if self._scaleFactorDer == "AUTO": return self.scaleFactorPivot return self._scaleFactorDer @scaleFactorDer.setter def scaleFactorDer(self, scaleFactorDer): if isinstance(scaleFactorDer, (str,)): scaleFactorDer = scaleFactorDer.upper() elif isinstance(scaleFactorDer, Iterable): scaleFactorDer = list(scaleFactorDer) self._scaleFactorDer = scaleFactorDer self._approxParameters["scaleFactorDer"] = self._scaleFactorDer @property def polydegreetype(self): """Value of polydegreetype.""" return "TOTAL" @polydegreetype.setter def polydegreetype(self, polydegreetype): RROMPyWarning(("polydegreetype is used just to simplify inheritance, " "and its value cannot be changed from 'TOTAL'.")) def _setupInterpolationIndices(self): """Setup parameters for polyvander.""" RROMPyAssert(self._mode, message = "Cannot setup interpolation indices.") if (self._musUniqueCN is None or len(self._reorder) != len(self.musPivot)): try: muPC = self.trainedModel.centerNormalizePivot(self.musPivot) except: muPC = self.trainedModel.centerNormalize(self.musPivot) self._musUniqueCN, musIdxsTo, musIdxs, musCount = (muPC.unique( return_index = True, return_inverse = True, return_counts = True)) self._musUnique = self.musPivot[musIdxsTo] self._derIdxs = [None] * len(self._musUniqueCN) self._reorder = np.empty(len(musIdxs), dtype = int) filled = 0 for j, cnt in enumerate(musCount): self._derIdxs[j] = nextDerivativeIndices([], self.nparPivot, cnt) jIdx = np.nonzero(musIdxs == j)[0] self._reorder[jIdx] = np.arange(filled, filled + cnt) filled += cnt def setupApprox(self) -> int: """Compute rational interpolant.""" if self.checkComputedApprox(): return -1 RROMPyAssert(self._mode, message = "Cannot setup approximant.") vbMng(self, "INIT", "Setting up {}.". format(self.name()), 5) self.computeScaleFactor() self.resetSamples() self.samplingEngine.scaleFactor = self.scaleFactorDer self.musPivot = self.samplerPivot.generatePoints(self.S) while len(self.musPivot) > self.S: self.musPivot.pop() self._musMarginal = self.samplerMarginal.generatePoints(self.SMarginal) while len(self.musMarginal) > self.SMarginal: self.musMarginal.pop() self._mus = emptyParameterList() self.mus.reset((self.S * self.SMarginal, self.HFEngine.npar)) for j, muMarg in enumerate(self.musMarginal): for k in range(j * self.S, (j + 1) * self.S): muk = np.empty_like(self.mus[0]) muk[self.directionPivot] = self.musPivot[k - j * self.S] muk[self.directionMarginal] = muMarg self.mus[k] = muk N0 = copy(self.N) self._setupTrainedModel(np.zeros((0, 0)), forceNew = True) idx, sizes = indicesScatter(len(self.musMarginal), return_sizes = True) pMat, Ps, Qs = None, [], [] req, emptyCores = [], np.where(np.logical_not(sizes))[0] if len(idx) == 0: vbMng(self, "MAIN", "Idling.", 30) if self.storeAllSamples: self.storeSamples() pL, pT = recv(source = 0, tag = poolRank()) pMat = np.empty((pL, 0), dtype = pT) else: _scaleFactorOldPivot = copy(self.scaleFactor) self.scaleFactor = self.scaleFactorPivot self._temporaryPivot = 1 for i in idx: musi = self.mus[self.S * i : self.S * (i + 1)] vbMng(self, "MAIN", "Building marginal model no. {} at {}.".format(i + 1, self.musMarginal[i]), 5) vbMng(self, "INIT", "Starting computation of snapshots.", 10) self.samplingEngine.resetHistory() self.samplingEngine.iterSample(musi) vbMng(self, "DEL", "Done computing snapshots.", 10) self.verbosity -= 5 self.samplingEngine.verbosity -= 10 self._setupRational(self._setupDenominator()) self.verbosity += 5 self.samplingEngine.verbosity += 10 if self.storeAllSamples: self.storeSamples(i) pMati = self.samplingEngine.projectionMatrix if not hasattr(self, "matchState") or not self.matchState: if self.POD == 1 and not ( hasattr(self.HFEngine.C, "is_mu_independent") and self.HFEngine.C.is_mu_independent in self._output_lvl): raise RROMPyException(("Cannot apply mu-dependent C " "to orthonormalized samples.")) vbMng(self, "INIT", "Extracting system output from state.", 35) if (hasattr(self.HFEngine, "_is_C_quadratic") and self.HFEngine._is_C_quadratic): pMati = self.HFEngine.applyC(pMati, musi) else: pMatiEff = None for j, mu in enumerate(musi): pMij = np.expand_dims(self.HFEngine.applyC( pMati[:, j], mu), -1) if pMatiEff is None: pMatiEff = np.array(pMij) else: pMatiEff = np.append(pMatiEff, pMij, axis = 1) pMati = pMatiEff vbMng(self, "DEL", "Done extracting system output.", 35) if pMat is None: pMat = copy(pMati) if i == 0: for dest in emptyCores: req += [isend((len(pMat), pMat.dtype), dest = dest, tag = dest)] else: pMat = np.hstack((pMat, pMati)) Ps += [copy(self.trainedModel.data.P)] Qs += [copy(self.trainedModel.data.Q)] del self.trainedModel.data.Q, self.trainedModel.data.P self.N = N0 del self._temporaryPivot self.scaleFactor = _scaleFactorOldPivot for r in req: r.wait() pMat, Ps, Qs, _, _ = gatherPivotedApproximant(pMat, Ps, Qs, self.mus.data, sizes, self.polybasis, False) self._setupTrainedModel(pMat) self.trainedModel.data.Qs, self.trainedModel.data.Ps = Qs, Ps Psupp = np.arange(0, len(self.musMarginal) * self.S, self.S) self.trainedModel.data.Psupp = list(Psupp) self._preliminaryMarginalFinalization() self._finalizeMarginalization() vbMng(self, "DEL", "Done setting up approximant.", 5) return 0 class RationalInterpolantPivotedNoMatch(RationalInterpolantPivotedBase, GenericPivotedApproximantNoMatch): """ ROM pivoted rational interpolant (without pole matching) 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': kind of snapshots orthogonalization; allowed values include 0, 1/2, and 1; defaults to 1, i.e. POD; - 'scaleFactorDer': scaling factors for derivative computation; defaults to 'AUTO'; - '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; - 'polybasis': type of polynomial basis for pivot interpolation; defaults to 'MONOMIAL'; - 'M': degree of rational interpolant numerator; defaults to 'AUTO', i.e. maximum allowed; - 'N': degree of rational interpolant denominator; defaults to 'AUTO', i.e. maximum allowed; - 'radialDirectionalWeights': radial basis weights for pivot numerator; defaults to 1; - 'radialDirectionalWeightsAdapt': bounds for adaptive rescaling of radial basis weights; defaults to [-1, -1]; - 'radialDirectionalWeightsMarginal': radial basis weights for marginal interpolant; defaults to 1; - 'functionalSolve': strategy for minimization of denominator functional; allowed values include 'NORM', 'DOMINANT', 'NODAL', 'BARYCENTRIC_NORM', and 'BARYCENTRIC[_AVERAGE]' (check pdf in main folder for meaning); defaults to 'NORM'; - - 'interpRcond': tolerance for pivot interpolation; defaults to + - 'interpTol': tolerance for pivot interpolation; defaults to None; - - 'robustTol': tolerance for robust rational denominator - management; defaults to 0. + - 'QTol': tolerance for robust rational denominator management; + defaults to 0. Defaults to empty dict. 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. approxParameters: Dictionary containing values for main parameters of approximant. Recognized keys are in parameterList. parameterListSoft: Recognized keys of soft approximant parameters: - 'POD': kind of snapshots orthogonalization; - 'scaleFactorDer': scaling factors for derivative computation; - 'polybasis': type of polynomial basis for pivot interpolation; - 'M': degree of rational interpolant numerator; - 'N': degree of rational interpolant denominator; - 'radialDirectionalWeights': radial basis weights for pivot numerator; - 'radialDirectionalWeightsAdapt': bounds for adaptive rescaling of radial basis weights; - 'radialDirectionalWeightsMarginal': radial basis weights for marginal interpolant; - 'functionalSolve': strategy for minimization of denominator functional; - - 'interpRcond': tolerance for pivot interpolation; - - 'robustTol': tolerance for robust rational denominator - management. + - 'interpTol': tolerance for pivot interpolation; + - 'QTol': 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. verbosity: Verbosity level. POD: Kind of snapshots orthogonalization. scaleFactorDer: Scaling factors for derivative computation. 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. polybasis: Type of polynomial basis for pivot interpolation. M: Numerator degree of approximant. N: Denominator degree of approximant. radialDirectionalWeights: Radial basis weights for pivot numerator. radialDirectionalWeightsAdapt: Bounds for adaptive rescaling of radial basis weights. radialDirectionalWeightsMarginal: Radial basis weights for marginal interpolant. functionalSolve: Strategy for minimization of denominator functional. - interpRcond: Tolerance for pivot interpolation. - robustTol: Tolerance for robust rational denominator management. + interpTol: Tolerance for pivot interpolation. + QTol: Tolerance for robust rational denominator management. muBounds: list of bounds for pivot parameter values. muBoundsMarginal: list of bounds for marginal parameter values. samplingEngine: Sampling engine. uHF: High fidelity solution(s) with parameter(s) lastSolvedHF as sampleList. lastSolvedHF: Parameter(s) corresponding to last computed high fidelity solution(s) as parameterList. uApproxReduced: Reduced approximate solution(s) with parameter(s) lastSolvedApprox as sampleList. lastSolvedApproxReduced: Parameter(s) corresponding to last computed reduced approximate solution(s) as parameterList. uApprox: Approximate solution(s) with parameter(s) lastSolvedApprox as sampleList. lastSolvedApprox: Parameter(s) corresponding to last computed approximate solution(s) as parameterList. Q: Numpy 1D vector containing complex coefficients of approximant denominator. P: Numpy 2D vector whose columns are FE dofs of coefficients of approximant numerator. """ class RationalInterpolantPivotedPoleMatch(RationalInterpolantPivotedBase, GenericPivotedApproximantPoleMatch): """ ROM pivoted rational interpolant (with pole matching) 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': kind of snapshots orthogonalization; allowed values include 0, 1/2, and 1; defaults to 1, i.e. POD; - 'scaleFactorDer': scaling factors for derivative computation; defaults to 'AUTO'; - 'matchState': whether to match the system state rather than the system output; defaults to False; - 'matchingWeight': weight for pole matching optimization; defaults to 1; - - 'sharedRatio': required ratio of marginal points to share + - 'matchingShared': required ratio of marginal points to share resonance; defaults to 1.; - '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; - 'polybasis': type of polynomial basis for pivot interpolation; defaults to 'MONOMIAL'; - 'polybasisMarginal': type of polynomial basis for marginal interpolation; allowed values include 'MONOMIAL_*', 'CHEBYSHEV_*', 'LEGENDRE_*', 'NEARESTNEIGHBOR', and 'PIECEWISE_LINEAR_*'; defaults to 'MONOMIAL'; - 'paramsMarginal': dictionary of parameters for marginal interpolation; include: . 'MMarginal': degree of marginal interpolant; defaults to 'AUTO', i.e. maximum allowed; not for 'NEARESTNEIGHBOR' or 'PIECEWISE_LINEAR_*'; . 'nNeighborsMarginal': number of marginal nearest neighbors; defaults to 1; only for 'NEARESTNEIGHBOR'; . 'polydegreetypeMarginal': type of polynomial degree for marginal; defaults to 'TOTAL'; not for 'NEARESTNEIGHBOR' or 'PIECEWISE_LINEAR_*'; - . 'interpRcondMarginal': tolerance for marginal interpolation; + . 'interpTolMarginal': tolerance for marginal interpolation; defaults to None; not for 'NEARESTNEIGHBOR' or 'PIECEWISE_LINEAR_*'; . 'radialDirectionalWeightsMarginalAdapt': bounds for adaptive rescaling of marginal radial basis weights; only for radial basis. - 'M': degree of rational interpolant numerator; defaults to 'AUTO', i.e. maximum allowed; - 'N': degree of rational interpolant denominator; defaults to 'AUTO', i.e. maximum allowed; - 'radialDirectionalWeights': radial basis weights for pivot numerator; defaults to 1; - 'radialDirectionalWeightsAdapt': bounds for adaptive rescaling of radial basis weights; defaults to [-1, -1]; - 'radialDirectionalWeightsMarginal': radial basis weights for marginal interpolant; defaults to 1; - 'functionalSolve': strategy for minimization of denominator functional; allowed values include 'NORM', 'DOMINANT', 'NODAL', 'BARYCENTRIC_NORM', and 'BARYCENTRIC[_AVERAGE]' (check pdf in main folder for meaning); defaults to 'NORM'; - - 'interpRcond': tolerance for pivot interpolation; defaults to - None; - - 'robustTol': tolerance for robust rational denominator - management; defaults to 0. + - 'interpTol': tolerance for pivot interpolation; defaults to None; + - 'QTol': tolerance for robust rational denominator management; + defaults to 0. Defaults to empty dict. 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. approxParameters: Dictionary containing values for main parameters of approximant. Recognized keys are in parameterList. parameterListSoft: Recognized keys of soft approximant parameters: - 'POD': kind of snapshots orthogonalization; - 'scaleFactorDer': scaling factors for derivative computation; - 'matchState': whether to match the system state rather than the system output; - 'matchingWeight': weight for pole matching optimization; - - 'sharedRatio': required ratio of marginal points to share + - 'matchingShared': required ratio of marginal points to share resonance; - 'polybasis': type of polynomial basis for pivot interpolation; - 'polybasisMarginal': type of polynomial basis for marginal interpolation; - 'paramsMarginal': dictionary of parameters for marginal interpolation; include: . 'MMarginal': degree of marginal interpolant; . 'nNeighborsMarginal': number of marginal nearest neighbors; . 'polydegreetypeMarginal': type of polynomial degree for marginal; - . 'interpRcondMarginal': tolerance for marginal interpolation; + . 'interpTolMarginal': tolerance for marginal interpolation; . 'radialDirectionalWeightsMarginalAdapt': bounds for adaptive rescaling of marginal radial basis weights. - 'M': degree of rational interpolant numerator; - 'N': degree of rational interpolant denominator; - 'radialDirectionalWeights': radial basis weights for pivot numerator; - 'radialDirectionalWeightsAdapt': bounds for adaptive rescaling of radial basis weights; - 'radialDirectionalWeightsMarginal': radial basis weights for marginal interpolant; - 'functionalSolve': strategy for minimization of denominator functional; - - 'interpRcond': tolerance for pivot interpolation; - - 'robustTol': tolerance for robust rational denominator - management. + - 'interpTol': tolerance for pivot interpolation; + - 'QTol': 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. verbosity: Verbosity level. POD: Kind of snapshots orthogonalization. scaleFactorDer: Scaling factors for derivative computation. matchState: Whether to match the system state rather than the system output. matchingWeight: Weight for pole matching optimization. - sharedRatio: Required ratio of marginal points to share resonance. + matchingShared: Required ratio of marginal points to share resonance. 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. polybasis: Type of polynomial basis for pivot interpolation. polybasisMarginal: Type of polynomial basis for marginal interpolation. paramsMarginal: Dictionary of parameters for marginal interpolation. M: Numerator degree of approximant. N: Denominator degree of approximant. radialDirectionalWeights: Radial basis weights for pivot numerator. radialDirectionalWeightsAdapt: Bounds for adaptive rescaling of radial basis weights. radialDirectionalWeightsMarginal: Radial basis weights for marginal interpolant. functionalSolve: Strategy for minimization of denominator functional. - interpRcond: Tolerance for pivot interpolation. - robustTol: Tolerance for robust rational denominator management. + interpTol: Tolerance for pivot interpolation. + QTol: Tolerance for robust rational denominator management. muBounds: list of bounds for pivot parameter values. muBoundsMarginal: list of bounds for marginal parameter values. samplingEngine: Sampling engine. uHF: High fidelity solution(s) with parameter(s) lastSolvedHF as sampleList. lastSolvedHF: Parameter(s) corresponding to last computed high fidelity solution(s) as parameterList. uApproxReduced: Reduced approximate solution(s) with parameter(s) lastSolvedApprox as sampleList. lastSolvedApproxReduced: Parameter(s) corresponding to last computed reduced approximate solution(s) as parameterList. uApprox: Approximate solution(s) with parameter(s) lastSolvedApprox as sampleList. lastSolvedApprox: Parameter(s) corresponding to last computed approximate solution(s) as parameterList. 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, *args, **kwargs) -> int: if self.checkComputedApprox(): return -1 self.purgeparamsMarginal() setupOK = super().setupApprox(*args, **kwargs) if self.matchState: self._postApplyC() return setupOK diff --git a/rrompy/reduction_methods/pivoted/trained_model/convert_trained_model_pivoted.py b/rrompy/reduction_methods/pivoted/trained_model/convert_trained_model_pivoted.py new file mode 100644 index 0000000..5a4a5eb --- /dev/null +++ b/rrompy/reduction_methods/pivoted/trained_model/convert_trained_model_pivoted.py @@ -0,0 +1,68 @@ +# Copyright (C) 2018-2020 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_pivoted_rational_nomatch import ( + TrainedModelPivotedRationalNoMatch) +from .trained_model_pivoted_rational_polematch import ( + TrainedModelPivotedRationalPoleMatch) +from rrompy.utilities.base import verbosityManager as vbMng +from rrompy.utilities.exception_manager import RROMPyException, RROMPyWarning + +__all__ = ['convertTrainedModelPivoted'] + +def convertTrainedModelPivoted(model, outType, verbObj = None, + muteWarnings : bool = False): + if isinstance(model, outType): return model + if ((isinstance(model, TrainedModelPivotedRationalNoMatch) + and outType == TrainedModelPivotedRationalPoleMatch) + or (isinstance(model, TrainedModelPivotedRationalPoleMatch) + and outType == TrainedModelPivotedRationalNoMatch)): + return convertTrainedModelPivotedMatchUnmatch(model, outType, verbObj, + muteWarnings) + raise RROMPyException(("Model type or conversion output type not " + "recognized.")) + +def convertTrainedModelPivotedMatchUnmatch(model, outType, verbObj = None, + muteWarnings : bool = False): + if verbObj is not None: + sf, st = ["NoMatch", "PoleMatch"] + if outType == TrainedModelPivotedRationalPoleMatch: + msgw = "match poles, set up marginalInterp," + else: #if outType == TrainedModelPivotedRationalNoMatch: + st, sf = sf, st + msgw = "set up marginalInterp" + vbMng(verbObj, "INIT", + "Starting model conversion from {} to {} model.".format(sf, st), + 10) + excludeDataKey = ["marginalInterp", "approxParameters"] + if outType == TrainedModelPivotedRationalPoleMatch: + modelC = TrainedModelPivotedRationalPoleMatch() + msgw = "match poles, set up marginalInterp," + else: #if outType == TrainedModelPivotedRationalNoMatch: + modelC = TrainedModelPivotedRationalNoMatch() + msgw = "set up marginalInterp" + excludeDataKey += ["HIs", "suppEffPts", "suppEffIdx", "coeffsEff", + "polesEff"] + for key in model.__dict__.keys(): setattr(modelC, key, model.__dict__[key]) + for key in excludeDataKey: delattr(modelC.data, key) + if verbObj is not None: + vbMng(verbObj, "DEL", "Finished model conversion.", 10) + if not muteWarnings: + RROMPyWarning(("Model conversion result not yet fuctional: must stil " + "{} and assign approxParameters.").format(msgw)) + return modelC diff --git a/rrompy/reduction_methods/pivoted/trained_model/trained_model_pivoted_rational_polematch.py b/rrompy/reduction_methods/pivoted/trained_model/trained_model_pivoted_rational_polematch.py index 21ffc99..4a4f77a 100644 --- a/rrompy/reduction_methods/pivoted/trained_model/trained_model_pivoted_rational_polematch.py +++ b/rrompy/reduction_methods/pivoted/trained_model/trained_model_pivoted_rational_polematch.py @@ -1,503 +1,503 @@ # Copyright (C) 2018-2020 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 warnings import numpy as np from scipy.special import factorial as fact from scipy.sparse import csr_matrix, hstack, SparseEfficiencyWarning from collections.abc import Iterable from copy import deepcopy as copy from itertools import combinations from rrompy.reduction_methods.standard.trained_model.trained_model_rational \ import TrainedModelRational from .trained_model_pivoted_rational_nomatch import ( TrainedModelPivotedRationalNoMatch) from rrompy.utilities.base.types import (Np1D, Np2D, List, ListAny, paramVal, paramList, sampList, HFEng) from rrompy.utilities.base import verbosityManager as vbMng, freepar as fp from rrompy.utilities.numerical import dot from rrompy.utilities.numerical.point_matching import rationalFunctionMatching from rrompy.utilities.numerical.degree import reduceDegreeN from rrompy.utilities.poly_fitting.polynomial import (polybases as ppb, PolynomialInterpolator as PI) from rrompy.utilities.poly_fitting.radial_basis import (polybases as rbpb, RadialBasisInterpolator as RBI) from rrompy.utilities.poly_fitting.heaviside import (rational2heaviside, heavisideUniformShape, HeavisideInterpolator as HI) from rrompy.utilities.poly_fitting.nearest_neighbor import ( NearestNeighborInterpolator as NNI) from rrompy.utilities.poly_fitting.piecewise_linear import (sparsekinds, PiecewiseLinearInterpolator as PLI) from rrompy.utilities.exception_manager import RROMPyException, RROMPyAssert from rrompy.sampling import sampleList, emptySampleList __all__ = ['TrainedModelPivotedRationalPoleMatch'] class TrainedModelPivotedRationalPoleMatch(TrainedModelPivotedRationalNoMatch): """ ROM approximant evaluation for pivoted approximants based on interpolation of rational approximants (with pole matching). Attributes: Data: dictionary with all that can be pickled. """ def compress(self, collapse : bool = False, tol : float = 0., returnRMat : bool = False, **compressMatrixkwargs): Psupp = copy(self.data.Psupp) RMat = super().compress(collapse, tol, True, **compressMatrixkwargs) if RMat is None: return for obj, suppj in zip(self.data.HIs, Psupp): obj.postmultiplyTensorize(RMat.T[suppj : suppj + obj.shape[0]]) if hasattr(self, "_HIsExcl"): for obj, suppj in zip(self._HIsExcl, Psupp): obj.postmultiplyTensorize(RMat.T[suppj : suppj + obj.shape[0]]) if not hasattr(self, "_PsExcl"): self._PsuppExcl = [0] * len(self._PsuppExcl) if returnRMat: return RMat 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.mu0Marginal. Returns: Normalized parameter. """ mu = self.checkParameterListMarginal(mu) if mu0 is None: mu0 = self.checkParameterListMarginal( self.data.mu0(0, self.data.directionMarginal)) return (self.mapParameterList(mu, idx = self.data.directionMarginal) - self.mapParameterList(mu0, idx = self.data.directionMarginal) ) / [self.data.scaleFactor[x] for x in self.data.directionMarginal] def setupMarginalInterp(self, approx, interpPars:ListAny, extraPar = None): vbMng(self, "INIT", "Starting computation of marginal interpolator.", 12) musMCN = self.centerNormalizeMarginal(self.data.musMarginal) nM, pbM = len(musMCN), approx.polybasisMarginal if pbM in ppb + rbpb: if extraPar: approx._setMMarginalAuto() _MMarginalEff = approx.paramsMarginal["MMarginal"] if pbM in ppb: p = PI() elif pbM in rbpb: p = RBI() else: # if pbM in sparsekinds + ["NEARESTNEIGHBOR"]: if pbM == "NEARESTNEIGHBOR": p = NNI() else: # if pbM in sparsekinds: pllims = [[-1.] * self.data.nparMarginal, [1.] * self.data.nparMarginal] p = PLI() for ipts, pts in enumerate(self.data.suppEffPts): if len(pts) == 0: raise RROMPyException("Empty list of support points.") musMCNEff, valsEff = musMCN[pts], np.eye(len(pts)) if pbM in ppb + rbpb: if extraPar: if ipts > 0: verb = approx.verbosity approx.verbosity = 0 _musM = approx.musMarginal approx.musMarginal = musMCNEff approx._setMMarginalAuto() approx.musMarginal = _musM approx.verbosity = verb else: approx.paramsMarginal["MMarginal"] = reduceDegreeN( _MMarginalEff, len(musMCNEff), self.data.nparMarginal, approx.paramsMarginal["polydegreetypeMarginal"]) MMEff = approx.paramsMarginal["MMarginal"] while MMEff >= 0: wellCond, msg = p.setupByInterpolation(musMCNEff, valsEff, MMEff, *interpPars) vbMng(self, "MAIN", msg, 30) if wellCond: break vbMng(self, "MAIN", ("Polyfit is poorly conditioned. Reducing " "MMarginal by 1."), 35) MMEff -= 1 if MMEff < 0: raise RROMPyException(("Instability in computation of " "interpolant. Aborting.")) if (pbM in rbpb and len(interpPars) > 4 and "optimizeScalingBounds" in interpPars[4].keys()): interpPars[4]["optimizeScalingBounds"] = [-1., -1.] elif pbM == "NEARESTNEIGHBOR": if ipts > 0: interpPars[0] = 1 p.setupByInterpolation(musMCNEff, valsEff, *interpPars) elif ipts == 0: # and pbM in sparsekinds: p.setupByInterpolation(musMCNEff, valsEff, pllims, extraPar[pts], *interpPars) if ipts == 0: self.data.marginalInterp = copy(p) self.data.coeffsEff, self.data.polesEff = [], [] for hi, sup in zip(self.data.HIs, self.data.Psupp): cEff = hi.coeffs if (self.data._collapsed or self.supportEnd == cEff.shape[1]): cEff = copy(cEff) else: supC = self.supportEnd - sup - cEff.shape[1] cEff = hstack((csr_matrix((len(cEff), sup)), csr_matrix(cEff), csr_matrix((len(cEff), supC))), "csr") self.data.coeffsEff += [cEff] self.data.polesEff += [copy(hi.poles)] else: ptsBad = [i for i in range(nM) if i not in pts] idxBad = np.where(self.data.suppEffIdx == ipts)[0] warnings.simplefilter('ignore', SparseEfficiencyWarning) if pbM in sparsekinds: for ij, j in enumerate(ptsBad): nearest = pts[np.argmin(np.sum(np.abs(musMCNEff.data - np.tile(musMCN[j], [len(pts), 1]) ), axis = 1).flatten())] self.data.coeffsEff[j][idxBad] = copy( self.data.coeffsEff[nearest][idxBad]) self.data.polesEff[j][idxBad] = copy( self.data.polesEff[nearest][idxBad]) else: if (self.data._collapsed or self.supportEnd == cEff.shape[1]): cfBase = np.zeros((len(idxBad), cEff.shape[1]), dtype = cEff.dtype) else: cfBase = csr_matrix((len(idxBad), self.supportEnd), dtype = cEff.dtype) valMuMBad = p(musMCN[ptsBad]) for ijb, jb in enumerate(ptsBad): self.data.coeffsEff[jb][idxBad] = copy(cfBase) self.data.polesEff[jb][idxBad] = 0. for ij, j in enumerate(pts): val = valMuMBad[ij][ijb] if not np.isclose(val, 0.): self.data.coeffsEff[jb][idxBad] += (val * self.data.coeffsEff[j][idxBad]) self.data.polesEff[jb][idxBad] += (val * self.data.polesEff[j][idxBad]) warnings.filters.pop(0) if pbM in ppb + rbpb: approx.paramsMarginal["MMarginal"] = _MMarginalEff vbMng(self, "DEL", "Done computing marginal interpolator.", 12) def updateEffectiveSamples(self, exclude:List[int], *args, **kwargs): if hasattr(self, "_idxExcl"): for j, excl in enumerate(self._idxExcl): self.data.HIs.insert(excl, self._HIsExcl[j]) super().updateEffectiveSamples(exclude) self._HIsExcl = [] for excl in self._idxExcl[::-1]: self._HIsExcl = [self.data.HIs.pop(excl)] + self._HIsExcl poles = [hi.poles for hi in self.data.HIs] coeffs = [hi.coeffs for hi in self.data.HIs] self.initializeFromLists(poles, coeffs, self.data.Psupp, self.data.HIs[0].polybasis, *args, **kwargs) def initializeFromRational(self, matchingWeight:float, HFEngine:HFEng, is_state:bool): """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, self.data.Psupp, basis, matchingWeight, HFEngine, is_state) def initializeFromLists(self, poles:ListAny, coeffs:ListAny, supps:ListAny, basis:str, matchingWeight:float, HFEngine:HFEng, is_state:bool): """Initialize Heaviside representation.""" Ns = [len(pls) for pls in poles] poles, coeffs = heavisideUniformShape(poles, coeffs) root = Ns.index(len(poles[0])) if hasattr(self.data, "_is_C_quadratic") and self.data._is_C_quadratic: csizemax = np.max([len(c) for c in coeffs]) #csizemax = np.max([c.shape[1] for c in coeffs]) if self.data._is_C_quadratic == 1: csizemax = csizemax ** 2 else: # if self.data._is_C_quadratic == 2: csizemax = csizemax * (csizemax + 1) // 2 TrainedModelRational._setupQuadMapping(self, csizemax) projMapping = self.quad_mapping projMappingReal = self.data._is_C_quadratic == 2 else: projMapping, projMappingReal = None, False poles, coeffs = rationalFunctionMatching(poles, coeffs, self.data.musMarginal.data, matchingWeight, supps, self.data.projMat, HFEngine, is_state, root, projMapping, projMappingReal) self.data.HIs = [] for pls, cfs in zip(poles, coeffs): hsi = HI() hsi.poles = pls if len(cfs) == len(pls): cfs = np.pad(cfs, ((0, 1), (0, 0)), "constant") hsi.coeffs = cfs hsi.npar = 1 hsi.polybasis = basis self.data.HIs += [hsi] self.data.suppEffPts = [np.arange(len(self.data.HIs))] self.data.suppEffIdx = np.zeros(len(poles[0]), dtype = int) - def checkSharedRatio(self, shared:float) -> str: + def checkShared(self, shared:float) -> str: N = len(self.data.HIs[0].poles) M = len(self.data.HIs) goodLocPoles = np.array([np.logical_not(np.isinf(hi.poles) ) for hi in self.data.HIs]) self.data.suppEffPts = [np.arange(len(self.data.HIs))] self.data.suppEffIdx = np.zeros(N, dtype = int) if np.all(np.all(goodLocPoles)): return " No poles erased." goodGlobPoles = np.sum(goodLocPoles, axis = 0) goodEnoughPoles = goodGlobPoles >= max(1., 1. * shared * M) keepPole = np.where(goodEnoughPoles)[0] halfPole = np.where(np.logical_and(goodEnoughPoles, goodGlobPoles < M))[0] removePole = np.where(np.logical_not(goodEnoughPoles))[0] if len(removePole) > 0: keepCoeff = np.append(keepPole, np.arange(N, len(self.data.HIs[0].coeffs))) for hi in self.data.HIs: for j in removePole: if not np.isinf(hi.poles[j]): hi.coeffs[N, :] -= hi.coeffs[j, :] / hi.poles[j] hi.poles = hi.poles[keepPole] hi.coeffs = hi.coeffs[keepCoeff, :] for idxR in halfPole: pts = np.where(goodLocPoles[:, idxR])[0] idxEff = len(self.data.suppEffPts) for idEff, prevPts in enumerate(self.data.suppEffPts): if len(prevPts) == len(pts): if np.allclose(prevPts, pts): idxEff = idEff break if idxEff == len(self.data.suppEffPts): self.data.suppEffPts += [pts] self.data.suppEffIdx[idxR] = idxEff self.data.suppEffIdx = self.data.suppEffIdx[keepPole] return (" Hard-erased {} pole".format(len(removePole)) + "s" * (len(removePole) != 1) + " and soft-erased {} pole".format(len(halfPole)) + "s" * (len(halfPole) != 1) + ".") 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 = self.checkParameterList(mu) if (not hasattr(self, "lastSolvedApproxReduced") or self.lastSolvedApproxReduced != mu): vbMng(self, "INIT", "Evaluating approximant at mu = {}.".format(mu), 12) muP = self.centerNormalizePivot(mu(self.data.directionPivot)) muM = mu(self.data.directionMarginal) his = self.interpolateMarginalInterpolator(muM) for i, (mP, hi) in enumerate(zip(muP, his)): uAppR = hi(mP)[:, 0] if i == 0: uApproxR = np.empty((len(uAppR), len(mu)), dtype = uAppR.dtype) uApproxR[:, i] = uAppR self.uApproxReduced = sampleList(uApproxR) vbMng(self, "DEL", "Done evaluating approximant.", 12) self.lastSolvedApproxReduced = mu return self.uApproxReduced def interpolateMarginalInterpolator(self, mu : paramList = []) -> ListAny: """Obtain interpolated approximant interpolator.""" mu = self.checkParameterListMarginal(mu) vbMng(self, "INIT", "Interpolating marginal models at mu = {}.".format(mu), 95) his = [] muC = self.centerNormalizeMarginal(mu) mIvals = self.data.marginalInterp(muC) verb, self.verbosity = self.verbosity, 0 poless = self.interpolateMarginalPoles(mu, mIvals) coeffss = self.interpolateMarginalCoeffs(mu, mIvals) self.verbosity = verb for j in range(len(mu)): his += [HI()] his[-1].poles = poless[j] his[-1].coeffs = coeffss[j] his[-1].npar = 1 his[-1].polybasis = self.data.HIs[0].polybasis vbMng(self, "DEL", "Done interpolating marginal models.", 95) return his def interpolateMarginalPoles(self, mu : paramList = [], mIvals : Np2D = None) -> ListAny: """Obtain interpolated approximant poles.""" mu = self.checkParameterListMarginal(mu) vbMng(self, "INIT", "Interpolating marginal poles at mu = {}.".format(mu), 95) intMPoles = np.zeros((len(mu),) + self.data.polesEff[0].shape, dtype = self.data.polesEff[0].dtype) if mIvals is None: muC = self.centerNormalizeMarginal(mu) mIvals = self.data.marginalInterp(muC) for pEff, mI in zip(self.data.polesEff, mIvals): intMPoles += np.expand_dims(mI, - 1) * pEff vbMng(self, "DEL", "Done interpolating marginal poles.", 95) return intMPoles def interpolateMarginalCoeffs(self, mu : paramList = [], mIvals : Np2D = None) -> ListAny: """Obtain interpolated approximant coefficients.""" mu = self.checkParameterListMarginal(mu) vbMng(self, "INIT", "Interpolating marginal coefficients at mu = {}.".format(mu), 95) intMCoeffs = np.zeros((len(mu),) + self.data.coeffsEff[0].shape, dtype = self.data.coeffsEff[0].dtype) if mIvals is None: muC = self.centerNormalizeMarginal(mu) mIvals = self.data.marginalInterp(muC) for cEff, mI in zip(self.data.coeffsEff, mIvals): for j, m in enumerate(mI): intMCoeffs[j] += m * cEff vbMng(self, "DEL", "Done interpolating marginal coefficients.", 95) return intMCoeffs def getPVal(self, mu : paramList = []) -> sampList: """ Evaluate rational numerator at arbitrary parameter. Args: mu: Target parameter. """ RROMPyAssert(self.data.nparPivot, 1, "Number of pivot parameters") mu = self.checkParameterList(mu) p = emptySampleList() muP = self.centerNormalizePivot(mu(self.data.directionPivot)) muM = mu(self.data.directionMarginal) his = self.interpolateMarginalInterpolator(muM) for i, (mP, hi) in enumerate(zip(muP, his)): Pval = hi(mP) * np.prod(mP[0] - hi.poles) if i == 0: p.reset((len(Pval), len(mu)), dtype = Pval.dtype) p[i] = Pval 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. """ RROMPyAssert(self.data.nparPivot, 1, "Number of pivot parameters") mu = self.checkParameterList(mu) muP = self.centerNormalizePivot(mu(self.data.directionPivot)) muM = mu(self.data.directionMarginal) 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 = np.complex) pls = self.interpolateMarginalPoles(muM) for i, (mP, pl) in enumerate(zip(muP, pls)): N = len(pl) if derP == N: derVal[i] = 1. elif derP >= 0 and derP < N: plDist = mP[0] - pl for terms in combinations(np.arange(N), N - derP): derVal[i] += np.prod(plDist[list(terms)]) return sclP ** derP * fact(derP) * derVal 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 isinstance(mVals, Iterable): mVals = [mVals] mVals = list(mVals) else: mVals = [fp] rDim = mVals.index(fp) if rDim < len(mVals) - 1 and fp in mVals[rDim + 1 :]: 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)[0] mMarg = [mVals[j] for j in range(len(mVals)) if j != rDim] roots = (self.data.scaleFactor[rDim] * self.interpolateMarginalPoles(mMarg)[0]) return self.mapParameterList(self.mapParameterList(self.data.mu0(rDim), idx = [rDim])(0, 0) + roots, "B", [rDim])(0) def getResidues(self, *args, **kwargs) -> Np2D: """ Obtain approximant residues. Returns: Numpy matrix with residues as columns. """ pls = self.getPoles(*args, **kwargs) if len(args) == 1: mVals = args[0] elif len(args) == 0: mVals = [None] else: mVals = kwargs["marginalVals"] if not isinstance(mVals, Iterable): mVals = [mVals] mVals = list(mVals) rDim = mVals.index(fp) mMarg = [mVals[j] for j in range(len(mVals)) if j != rDim] res = self.interpolateMarginalCoeffs(mMarg)[0][: len(pls), :].T if hasattr(self.data, "_is_C_quadratic") and self.data._is_C_quadratic: self._setupQuadMapping() res = res[self.quad_mapping[0]] * res[self.quad_mapping[1]].conj() if not self.data._collapsed: res = dot(self.data.projMat, res).T if (hasattr(self.data, "_is_C_quadratic") and self.data._is_C_quadratic == 2): res = np.real(res) return pls, res diff --git a/rrompy/reduction_methods/standard/greedy/rational_interpolant_greedy.py b/rrompy/reduction_methods/standard/greedy/rational_interpolant_greedy.py index d8fa57d..794e59a 100644 --- a/rrompy/reduction_methods/standard/greedy/rational_interpolant_greedy.py +++ b/rrompy/reduction_methods/standard/greedy/rational_interpolant_greedy.py @@ -1,531 +1,528 @@ # Copyright (C) 2018-2020 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.hfengines.base.linear_affine_engine import checkIfAffine from .generic_greedy_approximant import GenericGreedyApproximant from rrompy.utilities.poly_fitting.polynomial import (polybases, polyfitname, PolynomialInterpolator as PI, polyvander) from rrompy.utilities.numerical import dot from rrompy.utilities.numerical.degree import totalDegreeN from rrompy.utilities.expression import expressionEvaluator from rrompy.reduction_methods.standard import RationalInterpolant from rrompy.utilities.base.types import Np1D, Tuple, paramVal, List from rrompy.utilities.base.verbosity_depth import (verbosityManager as vbMng, getVerbosityDepth, setVerbosityDepth) from rrompy.utilities.poly_fitting import customFit from rrompy.utilities.exception_manager import (RROMPyWarning, RROMPyException, RROMPyAssert, RROMPy_FRAGILE) from rrompy.sampling import sampleList, emptySampleList __all__ = ['RationalInterpolantGreedy'] class RationalInterpolantGreedy(GenericGreedyApproximant, RationalInterpolant): """ ROM greedy rational interpolant 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': kind of snapshots orthogonalization; allowed values include 0, 1/2, and 1; defaults to 1, i.e. POD; - 'scaleFactorDer': scaling factors for derivative computation; defaults to 'AUTO'; - 'S': number of starting training points; - 'sampler': sample point generator; - 'greedyTol': uniform error tolerance for greedy algorithm; defaults to 1e-2; - 'collinearityTol': collinearity tolerance for greedy algorithm; defaults to 0.; - 'maxIter': maximum number of greedy steps; defaults to 1e2; - 'nTestPoints': number of test points; defaults to 5e2; - 'trainSetGenerator': training sample points generator; defaults to sampler; - 'polybasis': type of basis for interpolation; defaults to 'MONOMIAL'; - 'errorEstimatorKind': kind of error estimator; available values include 'AFFINE', 'DISCREPANCY', 'LOOK_AHEAD', 'LOOK_AHEAD_RES', 'LOOK_AHEAD_OUTPUT', and 'NONE'; defaults to 'NONE'; - 'functionalSolve': strategy for minimization of denominator functional; allowed values include 'NORM', 'DOMINANT', 'NODAL', 'BARYCENTRIC_NORM', and 'BARYCENTRIC[_AVERAGE]' (check pdf in main folder for meaning); defaults to 'NORM'; - - 'interpRcond': tolerance for interpolation; defaults to None; - - 'robustTol': tolerance for robust rational denominator - management; defaults to 0. + - 'interpTol': tolerance for interpolation; defaults to None; + - 'QTol': tolerance for robust rational denominator management; + defaults to 0. Defaults to empty dict. verbosity(optional): Verbosity level. Defaults to 10. Attributes: HFEngine: HF problem solver. mu0: Default parameter. mus: Array of snapshot parameters. approxParameters: Dictionary containing values for main parameters of approximant. Recognized keys are in parameterList. parameterListSoft: Recognized keys of soft approximant parameters: - 'POD': kind of snapshots orthogonalization; - 'scaleFactorDer': scaling factors for derivative computation; - 'greedyTol': uniform error tolerance for greedy algorithm; - 'collinearityTol': collinearity tolerance for greedy algorithm; - 'maxIter': maximum number of greedy steps; - 'nTestPoints': number of test points; - 'trainSetGenerator': training sample points generator; - 'errorEstimatorKind': kind of error estimator; - 'functionalSolve': strategy for minimization of denominator functional; - - 'interpRcond': tolerance for interpolation; - - 'robustTol': tolerance for robust rational denominator - management. + - 'interpTol': tolerance for interpolation; + - 'QTol': tolerance for robust rational denominator management. parameterListCritical: Recognized keys of critical approximant parameters: - 'S': total number of samples current approximant relies upon; - 'sampler': sample point generator. verbosity: Verbosity level. POD: Kind of snapshots orthogonalization. scaleFactorDer: Scaling factors for derivative computation. S: number of test points. sampler: Sample point generator. greedyTol: uniform error tolerance for greedy algorithm. collinearityTol: Collinearity tolerance for greedy algorithm. maxIter: maximum number of greedy steps. nTestPoints: number of starting training points. trainSetGenerator: training sample points generator. - robustTol: tolerance for robust rational denominator management. errorEstimatorKind: kind of error estimator. functionalSolve: Strategy for minimization of denominator functional. - interpRcond: tolerance for interpolation. - robustTol: tolerance for robust rational denominator management. + interpTol: tolerance for interpolation. + QTol: tolerance for robust rational denominator management. 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. """ _allowedEstimatorKinds = ["AFFINE", "DISCREPANCY", "LOOK_AHEAD", "LOOK_AHEAD_RES", "LOOK_AHEAD_OUTPUT", "NONE"] def __init__(self, *args, **kwargs): self._preInit() self._addParametersToList(["errorEstimatorKind"], ["DISCREPANCY"], toBeExcluded = ["M", "N", "polydegreetype", "radialDirectionalWeights"]) super().__init__(*args, **kwargs) self._postInit() @property def E(self): """Value of E.""" self._E = self.sampleBatchIdx - 1 return self._E @E.setter def E(self, E): RROMPyWarning(("E is used just to simplify inheritance, and its value " "cannot be changed from that of sampleBatchIdx - 1.")) def _setMAuto(self): self.M = self.E def _setNAuto(self): self.N = self.E @property def polydegreetype(self): """Value of polydegreetype.""" return "TOTAL" @polydegreetype.setter def polydegreetype(self, polydegreetype): RROMPyWarning(("polydegreetype is used just to simplify inheritance, " "and its value cannot be changed from 'TOTAL'.")) @property def polybasis(self): """Value of polybasis.""" return self._polybasis @polybasis.setter def polybasis(self, polybasis): try: polybasis = polybasis.upper().strip().replace(" ","") if polybasis not in polybases: raise RROMPyException("Sample type not recognized.") self._polybasis = polybasis except: RROMPyWarning(("Prescribed polybasis not recognized. Overriding " "to 'MONOMIAL'.")) self._polybasis = "MONOMIAL" self._approxParameters["polybasis"] = self.polybasis @property def errorEstimatorKind(self): """Value of errorEstimatorKind.""" return self._errorEstimatorKind @errorEstimatorKind.setter def errorEstimatorKind(self, errorEstimatorKind): errorEstimatorKind = errorEstimatorKind.upper() if errorEstimatorKind not in self._allowedEstimatorKinds: RROMPyWarning(("Error estimator kind not recognized. Overriding " "to 'NONE'.")) errorEstimatorKind = "NONE" self._errorEstimatorKind = errorEstimatorKind self._approxParameters["errorEstimatorKind"] = self.errorEstimatorKind def _polyvanderAuxiliary(self, mus, deg, *args): return polyvander(mus, deg, *args) def getErrorEstimatorDiscrepancy(self, mus:Np1D) -> Np1D: """Discrepancy-based residual estimator.""" checkIfAffine(self.HFEngine, "apply discrepancy-based error estimator", False, self._affine_lvl) mus = self.checkParameterList(mus) muCTest = self.trainedModel.centerNormalize(mus) tMverb, self.trainedModel.verbosity = self.trainedModel.verbosity, 0 QTest = self.trainedModel.getQVal(mus) QTzero = np.where(QTest == 0.)[0] if len(QTzero) > 0: RROMPyWarning(("Adjusting estimator to avoid division by " "numerically zero denominator.")) QTest[QTzero] = np.finfo(np.complex).eps / (1. + self.N) self.HFEngine.buildA() self.HFEngine.buildb() nAs, nbs = self.HFEngine.nAs, self.HFEngine.nbs muTrainEff = self.mapParameterList(self.mus) muTestEff = self.mapParameterList(mus) PTrain = self.trainedModel.getPVal(self.mus).data.T QTrain = self.trainedModel.getQVal(self.mus) QTzero = np.where(QTrain == 0.)[0] if len(QTzero) > 0: RROMPyWarning(("Adjusting estimator to avoid division by " "numerically zero denominator.")) QTrain[QTzero] = np.finfo(np.complex).eps / (1. + self.N) PTest = self.trainedModel.getPVal(mus).data self.trainedModel.verbosity = tMverb radiusAbTrain = np.empty((self.S, nAs * self.S + nbs), dtype = np.complex) radiusA = np.empty((self.S, nAs, len(mus)), dtype = np.complex) radiusb = np.empty((nbs, len(mus)), dtype = np.complex) for j, thA in enumerate(self.HFEngine.thAs): idxs = j * self.S + np.arange(self.S) radiusAbTrain[:, idxs] = expressionEvaluator(thA[0], muTrainEff, (self.S, 1)) * PTrain radiusA[:, j] = PTest * expressionEvaluator(thA[0], muTestEff, (len(mus),)) for j, thb in enumerate(self.HFEngine.thbs): idx = nAs * self.S + j radiusAbTrain[:, idx] = QTrain * expressionEvaluator(thb[0], muTrainEff, (self.S,)) radiusb[j] = QTest * expressionEvaluator(thb[0], muTestEff, (len(mus),)) QRHSNorm2 = self._affineResidualMatricesContraction(radiusb) vanTrain = self._polyvanderAuxiliary(self._musUniqueCN, self.E, self.polybasis0, self._derIdxs, self._reorder) - interpPQ = customFit(vanTrain, radiusAbTrain, - rcond = self.interpRcond) + interpPQ = customFit(vanTrain, radiusAbTrain, rcond = self.interpTol) vanTest = self._polyvanderAuxiliary(muCTest, self.E, self.polybasis0) DradiusAb = vanTest.dot(interpPQ) radiusA = (radiusA - DradiusAb[:, : - nbs].reshape(len(mus), -1, self.S).T) radiusb = radiusb - DradiusAb[:, - nbs :].T ff, Lf, LL = self._affineResidualMatricesContraction(radiusb, radiusA) err = np.abs((LL - 2. * np.real(Lf) + ff) / QRHSNorm2) ** .5 return err def getErrorEstimatorLookAhead(self, mus:Np1D, what : str = "") -> Tuple[Np1D, List[int]]: """Residual estimator based on look-ahead idea.""" errTest, QTest, idxMaxEst = self._EIMStep(mus) mu_muTestS = mus[idxMaxEst] app_muTestSample = self.getApproxReduced(mu_muTestS) if self._mode == RROMPy_FRAGILE: if what == "RES" and not self.HFEngine.isCEye: raise RROMPyException(("Cannot compute LOOK_AHEAD_RES " "estimator in fragile mode for " "non-scalar C.")) app_muTestSample = dot(self.trainedModel.data.projMat[:, : app_muTestSample.shape[0]], app_muTestSample) else: app_muTestSample = dot(self.samplingEngine.projectionMatrix, app_muTestSample) app_muTestSample = sampleList(app_muTestSample) if what == "RES": errmu = self.HFEngine.residual(mu_muTestS, app_muTestSample, post_c = False) solmu = self.HFEngine.residual(mu_muTestS, None, post_c = False) normSol = self.HFEngine.norm(solmu, dual = True) normErr = self.HFEngine.norm(errmu, dual = True) else: applyCglob = (hasattr(self.HFEngine, "_is_C_quadratic") and self.HFEngine._is_C_quadratic) for j, mu in enumerate(mu_muTestS): uEx = self.samplingEngine.nextSample(mu) if what == "OUTPUT" and not applyCglob: uEx = self.HFEngine.applyC(uEx, mu) app_muTS = self.HFEngine.applyC(app_muTestSample[j], mu) if j == 0: app_muTestS = emptySampleList() app_muTestS.reset((len(app_muTS), len(mu_muTestS)), dtype = app_muTS.dtype) app_muTestS[j] = app_muTS if j == 0: solmu = emptySampleList() solmu.reset((len(uEx), len(mu_muTestS)), dtype = uEx.dtype) solmu[j] = uEx if what == "OUTPUT": if applyCglob: solmu = sampleList(self.HFEngine.applyC(solmu, mu_muTestS)) app_muTestS = sampleList(self.HFEngine.applyC( app_muTestSample, mu_muTestS)) app_muTestSample = app_muTestS errmu = solmu - app_muTestSample normSol = self.HFEngine.norm(solmu, is_state = what != "OUTPUT") normErr = self.HFEngine.norm(errmu, is_state = what != "OUTPUT") errsamples = normErr / normSol musT = copy(self.mus) musT.append(mu_muTestS) musT = self.trainedModel.centerNormalize(musT) musC = self.trainedModel.centerNormalize(mus) errT = np.zeros((len(musT), len(mu_muTestS)), dtype = np.complex) errT[np.arange(len(self.mus), len(musT)), np.arange(len(mu_muTestS))] = errsamples * QTest[idxMaxEst] vanT = self._polyvanderAuxiliary(musT, self.E + 1, self.polybasis) - fitOut = customFit(vanT, errT, full = True, rcond = self.interpRcond) + fitOut = customFit(vanT, errT, full = True, rcond = self.interpTol) vbMng(self, "MAIN", ("Fitting {} samples with degree {} through {}... Conditioning " "of LS system: {:.4e}.").format(len(vanT), self.E + 1, polyfitname(self.polybasis), fitOut[1][2][0] / fitOut[1][2][-1]), 15) vanC = self._polyvanderAuxiliary(musC, self.E + 1, self.polybasis) err = np.sum(np.abs(vanC.dot(fitOut[0])), axis = -1) / QTest return err, idxMaxEst def getErrorEstimatorNone(self, mus:Np1D) -> Np1D: """EIM-based residual estimator.""" err = np.max(self._EIMStep(mus, True), axis = 1) err *= self.greedyTol / np.mean(err) return err def _EIMStep(self, mus:Np1D, only_one : bool = False) -> Tuple[Np1D, Np1D, List[int]]: """Residual estimator based on look-ahead idea.""" mus = self.checkParameterList(mus) tMverb, self.trainedModel.verbosity = self.trainedModel.verbosity, 0 QTest = self.trainedModel.getQVal(mus) QTzero = np.where(QTest == 0.)[0] if len(QTzero) > 0: RROMPyWarning(("Adjusting estimator to avoid division by " "numerically zero denominator.")) QTest[QTzero] = np.finfo(np.complex).eps / (1. + self.N) QTest = np.abs(QTest) muCTest = self.trainedModel.centerNormalize(mus) muCTrain = self.trainedModel.centerNormalize(self.mus) self.trainedModel.verbosity = tMverb vanTest = self._polyvanderAuxiliary(muCTest, self.E, self.polybasis) vanTestNext = self._polyvanderAuxiliary(muCTest, self.E + 1, self.polybasis)[:, vanTest.shape[1] :] idxsTest = np.arange(vanTestNext.shape[1]) basis = np.zeros((len(idxsTest), 0), dtype = float) idxMaxEst = [] while len(idxsTest) > 0: vanTrial = self._polyvanderAuxiliary(muCTrain, self.E, self.polybasis) vanTrialNext = self._polyvanderAuxiliary(muCTrain, self.E + 1, self.polybasis)[:, vanTrial.shape[1] :] vanTrial = np.hstack((vanTrial, vanTrialNext.dot(basis).reshape( len(vanTrialNext), basis.shape[1]))) valuesTrial = vanTrialNext[:, idxsTest] vanTestEff = np.hstack((vanTest, vanTestNext.dot(basis).reshape( len(vanTestNext), basis.shape[1]))) vanTestNextEff = vanTestNext[:, idxsTest] coeffTest = np.linalg.solve(vanTrial, valuesTrial) errTest = (np.abs(vanTestNextEff - vanTestEff.dot(coeffTest)) / np.expand_dims(QTest, 1)) if only_one: return errTest idxMaxErr = np.unravel_index(np.argmax(errTest), errTest.shape) idxMaxEst += [idxMaxErr[0]] muCTrain.append(muCTest[idxMaxErr[0]]) basis = np.pad(basis, [(0, 0), (0, 1)], "constant") basis[idxsTest[idxMaxErr[1]], -1] = 1. idxsTest = np.delete(idxsTest, idxMaxErr[1]) return errTest, QTest, idxMaxEst def errorEstimator(self, mus:Np1D, return_max : bool = False) -> Np1D: """Standard residual-based error estimator.""" setupOK = self.setupApproxLocal() if setupOK > 0: err = np.empty(len(mus)) err[:] = np.nan if not return_max: return err return err, [- setupOK], np.nan mus = self.checkParameterList(mus) vbMng(self.trainedModel, "INIT", "Evaluating error estimator at mu = {}.".format(mus), 10) if self.errorEstimatorKind == "AFFINE": err = self.getErrorEstimatorAffine(mus) else: self._setupInterpolationIndices() if self.errorEstimatorKind == "DISCREPANCY": err = self.getErrorEstimatorDiscrepancy(mus) elif self.errorEstimatorKind[: 10] == "LOOK_AHEAD": err, idxMaxEst = self.getErrorEstimatorLookAhead(mus, self.errorEstimatorKind[11 :]) else: #if self.errorEstimatorKind == "NONE": err = self.getErrorEstimatorNone(mus) vbMng(self.trainedModel, "DEL", "Done evaluating error estimator.", 10) if not return_max: return err if self.errorEstimatorKind[: 10] != "LOOK_AHEAD": idxMaxEst = np.empty(self.sampleBatchSize, dtype = int) errCP = copy(err) for j in range(self.sampleBatchSize): k = np.argmax(errCP) idxMaxEst[j] = k if j + 1 < self.sampleBatchSize: musZero = self.trainedModel.centerNormalize(mus, mus[k]) errCP *= np.linalg.norm(musZero.data, axis = 1) return err, idxMaxEst, err[idxMaxEst] def plotEstimator(self, *args, **kwargs): super().plotEstimator(*args, **kwargs) if self.errorEstimatorKind == "NONE": vbMng(self, "MAIN", ("Warning! Error estimator has been arbitrarily normalized " "before plotting."), 15) def greedyNextSample(self, *args, **kwargs) -> Tuple[Np1D, int, float, paramVal]: """Compute next greedy snapshot of solution map.""" RROMPyAssert(self._mode, message = "Cannot add greedy sample.") self.sampleBatchIdx += 1 self.sampleBatchSize = totalDegreeN(self.npar - 1, self.sampleBatchIdx) err, muidx, maxErr, muNext = super().greedyNextSample(*args, **kwargs) if maxErr is not None and (np.any(np.isnan(maxErr)) or np.any(np.isinf(maxErr))): self.sampleBatchIdx -= 1 self.sampleBatchSize = totalDegreeN(self.npar - 1, self.sampleBatchIdx) if (self.errorEstimatorKind == "NONE" and not np.isnan(maxErr) and not np.isinf(maxErr)): maxErr = None return err, muidx, maxErr, muNext def _setSampleBatch(self, maxS:int): self.sampleBatchIdx, self.sampleBatchSize, S = -1, 0, 0 nextBatchSize = 1 while S + nextBatchSize <= maxS: self.sampleBatchIdx += 1 self.sampleBatchSize = nextBatchSize S += self.sampleBatchSize nextBatchSize = totalDegreeN(self.npar - 1, self.sampleBatchIdx + 1) return S def _preliminaryTraining(self): """Initialize starting snapshots of solution map.""" RROMPyAssert(self._mode, message = "Cannot start greedy algorithm.") if self.samplingEngine.nsamples > 0: return self._S = self._setSampleBatch(self.S) super()._preliminaryTraining() self.M, self.N = ("AUTO",) * 2 def setupApproxLocal(self) -> int: """Compute rational interpolant.""" if self.checkComputedApprox(): return -1 RROMPyAssert(self._mode, message = "Cannot setup approximant.") self.verbosity -= 10 vbMng(self, "INIT", "Setting up local approximant.", 5) pMatOld, pMat = None, self.samplingEngine.projectionMatrix firstRun = self.trainedModel is None applyCglob = (hasattr(self.HFEngine, "_is_C_quadratic") and self.HFEngine._is_C_quadratic) if not firstRun: Sold = len(self.trainedModel.data.mus) if applyCglob: pMatOld = pMat[:, : Sold] pMat = pMat[:, Sold :] self._setupTrainedModel(pMat, not firstRun, pMatOld) self.catchInstability = 2 vbDepth = getVerbosityDepth() unstable = 0 if self.E > 0: try: Q = self._setupDenominator() except RROMPyException as RE: if RE.critical: raise RE from None setVerbosityDepth(vbDepth) RROMPyWarning("Downgraded {}: {}".format(RE.__class__.__name__, RE)) unstable = 1 else: Q = PI() Q.coeffs = np.ones((1,) * self.npar, dtype = np.complex) Q.npar = self.npar Q.polybasis = self.polybasis if not unstable: self.trainedModel.data.Q = copy(Q) try: P = copy(self._setupNumerator()) except RROMPyException as RE: if RE.critical: raise RE from None setVerbosityDepth(vbDepth) RROMPyWarning("Downgraded {}: {}".format(RE.__class__.__name__, RE)) unstable = 1 if not unstable: self.trainedModel.data.P = copy(P) self.trainedModel.data.approxParameters = copy( self.approxParameters) vbMng(self, "DEL", "Done setting up local approximant.", 5) self.catchInstability = 0 self.verbosity += 10 return unstable def setupApprox(self, plotEst : str = "NONE") -> int: val = super().setupApprox(plotEst) if val == 0: if (self.errorEstimatorKind[:10] == "LOOK_AHEAD" and len(self.mus) < self.samplingEngine.nsamples): while len(self.mus) < self.samplingEngine.nsamples: self.mus.append(self.samplingEngine.mus[len(self.mus)]) self.trainedModel = None self._S = self._setSampleBatch(len(self.mus) + 1) self.setupApproxLocal() self._setupRational(self.trainedModel.data.Q, self.trainedModel.data.P) self.trainedModel.data.approxParameters = copy( self.approxParameters) return val def loadTrainedModel(self, filename:str): """Load trained reduced model from file.""" super().loadTrainedModel(filename) self._setSampleBatch(self.S + 1) diff --git a/rrompy/reduction_methods/standard/rational_interpolant.py b/rrompy/reduction_methods/standard/rational_interpolant.py index fe98b4f..ca32c08 100644 --- a/rrompy/reduction_methods/standard/rational_interpolant.py +++ b/rrompy/reduction_methods/standard/rational_interpolant.py @@ -1,836 +1,816 @@ # Copyright (C) 2018-2020 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 scipy.linalg import eigvals from collections.abc import Iterable from .generic_standard_approximant import GenericStandardApproximant from rrompy.utilities.poly_fitting.polynomial import ( polybases as ppb, polyfitname, polyvander as pvP, polyTimes, PolynomialInterpolator as PI, PolynomialInterpolatorNodal as PIN) from rrompy.utilities.poly_fitting.heaviside import rational2heaviside from rrompy.utilities.poly_fitting.radial_basis import (polybases as rbpb, RadialBasisInterpolator as RBI) from rrompy.utilities.base.types import (Np1D, Np2D, Tuple, List, sampList, interpEng) from rrompy.utilities.base import verbosityManager as vbMng from rrompy.utilities.numerical import (pseudoInverse, dot, potential, distanceMatrix) from rrompy.utilities.numerical.factorials import multifactorial from rrompy.utilities.numerical.hash_derivative import (nextDerivativeIndices, hashDerivativeToIdx as hashD, hashIdxToDerivative as hashI) from rrompy.utilities.numerical.degree import (reduceDegreeN, degreeTotalToFull, fullDegreeMaxMask, totalDegreeMaxMask) from rrompy.solver import Np2DLikeGramian from rrompy.utilities.exception_manager import (RROMPyException, RROMPyAssert, RROMPyWarning) __all__ = ['RationalInterpolant'] def polyTimesTable(P:interpEng, mus:Np1D, reorder:List[int], derIdxs:List[List[List[int]]], scl : Np1D = None) -> Np2D: """Table of polynomial products.""" if not isinstance(P, PI): raise RROMPyException(("Polynomial to evaluate must be a polynomial " "interpolator.")) Pvals = [[0.] * len(derIdx) for derIdx in derIdxs] for j, derIdx in enumerate(derIdxs): nder = len(derIdx) for der in range(nder): derI = hashI(der, P.npar) Pvals[j][der] = P([mus[j]], derI, scl) / multifactorial(derI) return blockDiagDer(Pvals, reorder, derIdxs) def vanderInvTable(vanInv:Np2D, idxs:List[int], reorder:List[int], derIdxs:List[List[List[int]]]) -> Np2D: """Table of Vandermonde pseudo-inverse.""" S = len(reorder) Ts = [None] * len(idxs) for k in range(len(idxs)): invLocs = [None] * len(derIdxs) idxGlob = 0 for j, derIdx in enumerate(derIdxs): nder = len(derIdx) idxGlob += nder idxLoc = np.arange(S)[np.logical_and(reorder >= idxGlob - nder, reorder < idxGlob)] invLocs[j] = vanInv[k, idxLoc] Ts[k] = blockDiagDer(invLocs, reorder, derIdxs, [2, 1, 0]) return Ts def blockDiagDer(vals:List[Np1D], reorder:List[int], derIdxs:List[List[List[int]]], permute : List[int] = None) -> Np2D: """Table of derivative values for point confluence.""" S = len(reorder) T = np.zeros((S, S), dtype = np.complex) if permute is None: permute = [0, 1, 2] idxGlob = 0 for j, derIdx in enumerate(derIdxs): nder = len(derIdx) idxGlob += nder idxLoc = np.arange(S)[np.logical_and(reorder >= idxGlob - nder, reorder < idxGlob)] val = vals[j] for derI, derIdxI in enumerate(derIdx): for derJ, derIdxJ in enumerate(derIdx): diffIdx = [x - y for (x, y) in zip(derIdxI, derIdxJ)] if all([x >= 0 for x in diffIdx]): diffj = hashD(diffIdx) i1, i2, i3 = np.array([derI, derJ, diffj])[permute] T[idxLoc[i1], idxLoc[i2]] = val[i3] return T class RationalInterpolant(GenericStandardApproximant): """ ROM rational interpolant 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': kind of snapshots orthogonalization; allowed values include 0, 1/2, and 1; defaults to 1, i.e. POD; - 'scaleFactorDer': scaling factors for derivative computation; defaults to 'AUTO'; - 'S': total number of samples current approximant relies upon; - 'sampler': sample point generator; - 'polybasis': type of polynomial basis for interpolation; defaults to 'MONOMIAL'; - 'M': degree of rational interpolant numerator; defaults to 'AUTO', i.e. maximum allowed; - 'N': degree of rational interpolant denominator; defaults to 'AUTO', i.e. maximum allowed; - 'polydegreetype': type of polynomial degree; defaults to 'TOTAL'; - 'radialDirectionalWeights': radial basis weights for interpolant numerator; defaults to 1; - 'radialDirectionalWeightsAdapt': bounds for adaptive rescaling of radial basis weights; defaults to [-1, -1]; - 'functionalSolve': strategy for minimization of denominator functional; allowed values include 'NORM', 'DOMINANT', 'NODAL', 'BARYCENTRIC_NORM', and 'BARYCENTRIC[_AVERAGE]' (check pdf in main folder for explanation); defaults to 'NORM'; - - 'interpRcond': tolerance for interpolation; defaults to None; - - 'robustTol': tolerance for robust rational denominator - management; defaults to 0. + - 'interpTol': tolerance for interpolation; defaults to None; + - 'QTol': tolerance for robust rational denominator management; + defaults to 0. Defaults to empty dict. verbosity(optional): Verbosity level. Defaults to 10. Attributes: HFEngine: HF problem solver. mu0: Default parameter. mus: Array of snapshot parameters. approxParameters: Dictionary containing values for main parameters of approximant. Recognized keys are in parameterList. parameterListSoft: Recognized keys of soft approximant parameters: - 'POD': kind of snapshots orthogonalization; - 'scaleFactorDer': scaling factors for derivative computation; - 'polybasis': type of polynomial basis for interpolation; - 'M': degree of rational interpolant numerator; - 'N': degree of rational interpolant denominator; - 'polydegreetype': type of polynomial degree; - 'radialDirectionalWeights': radial basis weights for interpolant numerator; - 'radialDirectionalWeightsAdapt': bounds for adaptive rescaling of radial basis weights; - 'functionalSolve': strategy for minimization of denominator functional; - - 'interpRcond': tolerance for interpolation via numpy.polyfit; - - 'robustTol': tolerance for robust rational denominator - management. + - 'interpTol': tolerance for interpolation via numpy.polyfit; + - 'QTol': tolerance for robust rational denominator management. parameterListCritical: Recognized keys of critical approximant parameters: - 'S': total number of samples current approximant relies upon; - 'sampler': sample point generator. verbosity: Verbosity level. POD: Kind of snapshots orthogonalization. scaleFactorDer: Scaling factors for derivative computation. S: Number of solution snapshots over which current approximant is based upon. sampler: Sample point generator. polybasis: type of polynomial basis for interpolation. M: Numerator degree of approximant. N: Denominator degree of approximant. polydegreetype: Type of polynomial degree. radialDirectionalWeights: Radial basis weights for interpolant numerator. radialDirectionalWeightsAdapt: Bounds for adaptive rescaling of radial basis weights. functionalSolve: Strategy for minimization of denominator functional. - interpRcond: Tolerance for interpolation via numpy.polyfit. - robustTol: Tolerance for robust rational denominator management. + interpTol: Tolerance for interpolation via numpy.polyfit. + QTol: Tolerance for robust rational denominator management. 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. Q: Numpy 1D vector containing complex coefficients of approximant denominator. P: Numpy 2D vector whose columns are FE dofs of coefficients of approximant numerator. """ _allowedFunctionalSolveKinds = ["NORM", "DOMINANT", "NODAL", "BARYCENTRIC_NORM", "BARYCENTRIC_AVERAGE"] def __init__(self, *args, **kwargs): self._preInit() self._addParametersToList(["polybasis", "M", "N", "polydegreetype", "radialDirectionalWeights", "radialDirectionalWeightsAdapt", - "functionalSolve", "interpRcond", - "robustTol"], + "functionalSolve", "interpTol", "QTol"], ["MONOMIAL", "AUTO", "AUTO", "TOTAL", 1., [-1., -1.], "NORM", -1, 0.]) super().__init__(*args, **kwargs) self.catchInstability = 0 self._postInit() @property def tModelType(self): from .trained_model.trained_model_rational import TrainedModelRational return TrainedModelRational @property def polybasis(self): """Value of polybasis.""" return self._polybasis @polybasis.setter def polybasis(self, polybasis): try: polybasis = polybasis.upper().strip().replace(" ","") if polybasis not in ppb + rbpb: raise RROMPyException("Prescribed polybasis not recognized.") self._polybasis = polybasis except: RROMPyWarning(("Prescribed polybasis not recognized. Overriding " "to 'MONOMIAL'.")) self._polybasis = "MONOMIAL" self._approxParameters["polybasis"] = self.polybasis @property def polybasis0(self): if "_" in self.polybasis: return self.polybasis.split("_")[0] return self.polybasis @property def functionalSolve(self): """Value of functionalSolve.""" return self._functionalSolve @functionalSolve.setter def functionalSolve(self, functionalSolve): try: functionalSolve = functionalSolve.upper().strip().replace(" ","") if functionalSolve == "BARYCENTRIC": functionalSolve += "_AVERAGE" if functionalSolve not in self._allowedFunctionalSolveKinds: raise RROMPyException(("Prescribed functionalSolve not " "recognized.")) self._functionalSolve = functionalSolve except: RROMPyWarning(("Prescribed functionalSolve not recognized. " "Overriding to 'NORM'.")) self._functionalSolve = "NORM" self._approxParameters["functionalSolve"] = self.functionalSolve @property - def interpRcond(self): - """Value of interpRcond.""" - return self._interpRcond - @interpRcond.setter - def interpRcond(self, interpRcond): - self._interpRcond = interpRcond - self._approxParameters["interpRcond"] = self.interpRcond + def interpTol(self): + """Value of interpTol.""" + return self._interpTol + @interpTol.setter + def interpTol(self, interpTol): + self._interpTol = interpTol + self._approxParameters["interpTol"] = self.interpTol @property def radialDirectionalWeights(self): """Value of radialDirectionalWeights.""" return self._radialDirectionalWeights @radialDirectionalWeights.setter def radialDirectionalWeights(self, radialDirectionalWeights): if isinstance(radialDirectionalWeights, Iterable): radialDirectionalWeights = list(radialDirectionalWeights) else: radialDirectionalWeights = [radialDirectionalWeights] self._radialDirectionalWeights = radialDirectionalWeights self._approxParameters["radialDirectionalWeights"] = ( self.radialDirectionalWeights) @property def radialDirectionalWeightsAdapt(self): """Value of radialDirectionalWeightsAdapt.""" return self._radialDirectionalWeightsAdapt @radialDirectionalWeightsAdapt.setter def radialDirectionalWeightsAdapt(self, radialDirectionalWeightsAdapt): self._radialDirectionalWeightsAdapt = radialDirectionalWeightsAdapt self._approxParameters["radialDirectionalWeightsAdapt"] = ( self.radialDirectionalWeightsAdapt) @property def M(self): """Value of M.""" return self._M @M.setter def M(self, M): if isinstance(M, str): M = M.strip().replace(" ","") if "-" not in M: M = M + "-0" self._M_isauto, self._M_shift = True, int(M.split("-")[-1]) M = 0 if M < 0: raise RROMPyException("M must be non-negative.") self._M = M self._approxParameters["M"] = self.M def _setMAuto(self): self.M = max(0, reduceDegreeN(self.S, self.S, self.npar, self.polydegreetype) - self._M_shift) vbMng(self, "MAIN", "Automatically setting M to {}.".format(self.M), 25) @property def N(self): """Value of N.""" return self._N @N.setter def N(self, N): if isinstance(N, str): N = N.strip().replace(" ","") if "-" not in N: N = N + "-0" self._N_isauto, self._N_shift = True, int(N.split("-")[-1]) N = 0 if N < 0: raise RROMPyException("N must be non-negative.") self._N = N self._approxParameters["N"] = self.N def _setNAuto(self): self.N = max(0, reduceDegreeN(self.S, self.S, self.npar, self.polydegreetype) - self._N_shift) vbMng(self, "MAIN", "Automatically setting N to {}.".format(self.N), 25) @property def polydegreetype(self): """Value of polydegreetype.""" return self._polydegreetype @polydegreetype.setter def polydegreetype(self, polydegreetype): try: polydegreetype = polydegreetype.upper().strip().replace(" ","") if polydegreetype not in ["TOTAL", "FULL"]: raise RROMPyException(("Prescribed polydegreetype not " "recognized.")) self._polydegreetype = polydegreetype except: RROMPyWarning(("Prescribed polydegreetype not recognized. " "Overriding to 'TOTAL'.")) self._polydegreetype = "TOTAL" self._approxParameters["polydegreetype"] = self.polydegreetype @property - def robustTol(self): + def QTol(self): """Value of tolerance for robust rational denominator management.""" - return self._robustTol - @robustTol.setter - def robustTol(self, robustTol): - if robustTol < 0.: + return self._QTol + @QTol.setter + def QTol(self, QTol): + if QTol < 0.: RROMPyWarning(("Overriding prescribed negative robustness " "tolerance to 0.")) - robustTol = 0. - self._robustTol = robustTol - self._approxParameters["robustTol"] = self.robustTol + QTol = 0. + self._QTol = QTol + self._approxParameters["QTol"] = self.QTol def resetSamples(self): """Reset samples.""" super().resetSamples() self._musUniqueCN = None self._derIdxs = None self._reorder = None def _setupInterpolationIndices(self): """Setup parameters for polyvander.""" if self._musUniqueCN is None or len(self._reorder) != len(self.mus): self._musUniqueCN, musIdxsTo, musIdxs, musCount = ( self.trainedModel.centerNormalize(self.mus).unique( return_index = True, return_inverse = True, return_counts = True)) self._musUnique = self.mus[musIdxsTo] self._derIdxs = [None] * len(self._musUniqueCN) self._reorder = np.empty(len(musIdxs), dtype = int) filled = 0 for j, cnt in enumerate(musCount): self._derIdxs[j] = nextDerivativeIndices([], self.mus.shape[1], cnt) jIdx = np.nonzero(musIdxs == j)[0] self._reorder[jIdx] = np.arange(filled, filled + cnt) filled += cnt def _setupDenominator(self): """Compute rational denominator.""" RROMPyAssert(self._mode, message = "Cannot setup denominator.") vbMng(self, "INIT", "Starting computation of denominator.", 7) if hasattr(self, "_N_isauto"): self._setNAuto() else: N = reduceDegreeN(self.N, self.S, self.npar, self.polydegreetype) if N < self.N: RROMPyWarning(("N too large compared to S. Reducing N by " "{}").format(self.N - N)) self.N = N while self.N > 0: if self.functionalSolve != "NORM" and self.npar > 1: RROMPyWarning(("Strategy for functional optimization must be " "'NORM' for more than one parameter. " "Overriding to 'NORM'.")) self.functionalSolve = "NORM" if (self.functionalSolve[:11] == "BARYCENTRIC" and self.N + 1 < self.S): RROMPyWarning(("Barycentric strategy cannot be applied with " "Least Squares. Overriding to 'NORM'.")) self.functionalSolve = "NORM" if self.functionalSolve[:11] == "BARYCENTRIC": invD, TN = None, None self._setupInterpolationIndices() else: invD, TN = self._computeInterpolantInverseBlocks() if (self.functionalSolve in ["NODAL", "BARYCENTRIC_NORM", "BARYCENTRIC_AVERAGE"] - and len(self._musUnique) != len(self.mus)): + and len(self._musUnique) != self.S): if self.functionalSolve[:11] == "BARYCENTRIC": warnflag = "Barycentric" invD, TN = self._computeInterpolantInverseBlocks() else: warnflag = "Iterative" RROMPyWarning(("{} functional optimization cannot be applied " "to repeated samples. Overriding to " "'NORM'.").format(warnflag)) self.functionalSolve = "NORM" idxSamplesEff = list(range(self.S)) if self.POD == 1: ev, eV = self.findeveVGQR( self.samplingEngine.Rscale[:, idxSamplesEff], invD, TN) elif self.POD == 1/2: ev, eV = self.findeveVGExplicit( self.samplingEngine.samples_normal(idxSamplesEff), invD, TN, self.samplingEngine.Rscale) else: ev, eV = self.findeveVGExplicit( self.samplingEngine.samples(idxSamplesEff), invD, TN) if self.functionalSolve == "NODAL": break evR = ev / np.max(ev) - ts = self.robustTol * np.linalg.norm(evR) + ts = self.QTol * np.linalg.norm(evR) nevBad = len(ev) - np.sum(np.abs(evR) >= ts) if nevBad <= (self.functionalSolve == "NORM"): break if self.catchInstability: raise RROMPyException(("Instability in denominator " "computation: eigenproblem is poorly " "conditioned."), self.catchInstability == 1) vbMng(self, "MAIN", ("Smallest {} eigenvalues below tolerance. " "Reducing N by 1.").format(nevBad), 10) self.N = self.N - 1 if self.N <= 0: self.N = 0 eV = np.ones((1, 1)) if self.N > 0 and self.functionalSolve in ["NODAL", "BARYCENTRIC_NORM", "BARYCENTRIC_AVERAGE"]: q = PIN() q.polybasis, q.nodes = self.polybasis0, eV else: q = PI() q.npar = self.npar q.polybasis = self.polybasis0 if self.polydegreetype == "TOTAL": q.coeffs = degreeTotalToFull(tuple([self.N + 1] * self.npar), self.npar, eV) else: q.coeffs = eV.reshape([self.N + 1] * self.npar) vbMng(self, "DEL", "Done computing denominator.", 7) return q def _setupNumerator(self): """Compute rational numerator.""" RROMPyAssert(self._mode, message = "Cannot setup numerator.") vbMng(self, "INIT", "Starting computation of numerator.", 7) self._setupInterpolationIndices() Qevaldiag = polyTimesTable(self.trainedModel.data.Q, self._musUniqueCN, self._reorder, self._derIdxs, self.scaleFactorRel) if self.POD == 1: Qevaldiag = Qevaldiag.dot(self.samplingEngine.Rscale.T) elif self.POD == 1/2: Qevaldiag = Qevaldiag * self.samplingEngine.Rscale if hasattr(self, "_M_isauto"): self._setMAuto() M = self.M else: M = reduceDegreeN(self.M, self.S, self.npar, self.polydegreetype) if M < self.M: RROMPyWarning(("M too large compared to S. Reducing M by " "{}").format(self.M - M)) self.M = M while self.M >= 0: pParRest = [self.M, self.polybasis, self.verbosity >= 5, self.polydegreetype == "TOTAL", {"derIdxs": self._derIdxs, "reorder": self._reorder, "scl": self.scaleFactorRel}] if self.polybasis in ppb: p = PI() else: self.computeScaleFactor() rDWEff = np.array([w * f for w, f in zip( self.radialDirectionalWeights, self.scaleFactor)]) pParRest = pParRest[: 2] + [rDWEff] + pParRest[2 :] pParRest[-1]["optimizeScalingBounds"] = ( self.radialDirectionalWeightsAdapt) p = RBI() if self.polybasis in ppb + rbpb: - pParRest += [{"rcond": self.interpRcond}] + pParRest += [{"rcond": self.interpTol}] wellCond, msg = p.setupByInterpolation(self._musUniqueCN, Qevaldiag, *pParRest) vbMng(self, "MAIN", msg, 5) if wellCond: break if self.catchInstability: raise RROMPyException(("Instability in numerator computation: " "polyfit is poorly conditioned."), self.catchInstability == 1) vbMng(self, "MAIN", ("Polyfit is poorly conditioned. Reducing M " "by 1."), 10) self.M = self.M - 1 if self.M < 0: raise RROMPyException(("Instability in computation of numerator. " "Aborting.")) self.M = M vbMng(self, "DEL", "Done computing numerator.", 7) return p def setupApprox(self) -> int: """Compute rational interpolant.""" if self.checkComputedApprox(): return -1 RROMPyAssert(self._mode, message = "Cannot setup approximant.") vbMng(self, "INIT", "Setting up {}.". format(self.name()), 5) self.computeSnapshots() self._setupTrainedModel(self.samplingEngine.projectionMatrix) self._setupRational(self._setupDenominator()) self.trainedModel.data.approxParameters = copy(self.approxParameters) vbMng(self, "DEL", "Done setting up approximant.", 5) return 0 def _setupRational(self, Q:interpEng, P : interpEng = None): vbMng(self, "INIT", "Starting approximant finalization.", 5) self.trainedModel.data.Q = Q if P is None: P = self._setupNumerator() - if self.N > 0 and self.npar == 1: - #check for bad poles + while self.N > 0 and self.npar == 1: + #check for bad poles (absolute) pls = Q.roots() - idxBad = self.HFEngine.flagBadPolesResidues(pls, relative = True) - plsN = self.mapParameterList(self.mapParameterList(self.mu0)(0, 0) - + self.scaleFactor * pls, "B")(0) - idxBad = np.logical_or(self.HFEngine.flagBadPolesResidues(pls, - relative = True), - self.HFEngine.flagBadPolesResidues(plsN)) - if np.any(idxBad): - vbMng(self, "MAIN", - "Removing {} spurious poles out of {} due to poles."\ - .format(np.sum(idxBad), self.N), 10) - if isinstance(Q, PIN): - Q.nodes = Q.nodes[np.logical_not(idxBad)] - else: - Q = PI() - Q.npar = self.npar - Q.polybasis = self.polybasis0 - Q.coeffs = np.ones(1, dtype = np.complex) - for pl in pls[np.logical_not(idxBad)]: - Q.coeffs = polyTimes(Q.coeffs, [- pl, 1.], - Pbasis = Q.polybasis, Rbasis = Q.polybasis) - Q.coeffs /= np.linalg.norm(Q.coeffs) - self.trainedModel.data.Q = Q - self.N = Q.deg[0] - P = self._setupNumerator() - if (not hasattr(self.HFEngine, "_ignoreResidues") - or not self.HFEngine._ignoreResidues): - #check for bad residues - cfs, pls, _ = rational2heaviside(P, Q) - cfs = cfs[: self.N].T + plsA = self.mapParameterList(self.mapParameterList(self.mu0)(0, 0) + + self.scaleFactor * pls, "B")(0) + idxBadA = self.HFEngine.flagBadPolesResidues(plsA) + if self.HFEngine._ignoreResidues: + #check for bad poles (relative) + cfs, projMat = None, None + else: + #check for bad poles (relative) and residues + cfs = rational2heaviside(P, Q)[0][: self.N].T if self.POD != 1: - cfs = self.samplingEngine.projectionMatrix.dot(cfs) + projMat = self.samplingEngine.projectionMatrix + else: + projMat = None foci = self.sampler.normalFoci() ground = self.sampler.groundPotential() potEff = potential(pls, foci) / ground potEff[np.logical_or(potEff < 1., np.isinf(pls))] = 1. cfs[:, np.isinf(pls)] = 0. cfs /= potEff # rescale by potential - idxBad = self.HFEngine.flagBadPolesResidues(pls, cfs) - if np.any(idxBad): - vbMng(self, "MAIN", - ("Removing {} spurious poles out of {} due to " - "residues.").format(np.sum(idxBad), self.N), 10) - if isinstance(Q, PIN): - Q.nodes = Q.nodes[np.logical_not(idxBad)] - else: - Q = PI() - Q.npar = self.npar - Q.polybasis = self.polybasis0 - Q.coeffs = np.ones(1, dtype = np.complex) - for pl in pls[np.logical_not(idxBad)]: - Q.coeffs = polyTimes(Q.coeffs, [- pl, 1.], + idxBadR = self.HFEngine.flagBadPolesResidues(pls, cfs, True, + projMat) + idxBad = np.logical_or(idxBadA, idxBadR) + if not np.any(idxBad): break + vbMng(self, "MAIN", + "Removing {} spurious poles out of {}.".format( + np.sum(idxBad), self.N), 10) + if isinstance(Q, PIN): + Q.nodes = Q.nodes[np.logical_not(idxBad)] + else: + Q = PI() + Q.npar = self.npar + Q.polybasis = self.polybasis0 + Q.coeffs = np.ones(1, dtype = np.complex) + for pl in pls[np.logical_not(idxBad)]: + Q.coeffs = polyTimes(Q.coeffs, [- pl, 1.], Pbasis = Q.polybasis, Rbasis = Q.polybasis) - Q.coeffs /= np.linalg.norm(Q.coeffs) - self.trainedModel.data.Q = Q - self.N = Q.deg[0] - P = self._setupNumerator() + Q.coeffs /= np.linalg.norm(Q.coeffs) + self.trainedModel.data.Q = Q + self.N = Q.deg[0] + P = self._setupNumerator() self.trainedModel.data.P = P vbMng(self, "DEL", "Terminated approximant finalization.", 5) def _computeInterpolantInverseBlocks(self) -> Tuple[List[Np2D], Np2D]: """ Compute inverse factors for minimal interpolant target functional. """ RROMPyAssert(self._mode, message = "Cannot solve eigenvalue problem.") self._setupInterpolationIndices() pvPPar = [self.polybasis0, self._derIdxs, self._reorder, self.scaleFactorRel] if hasattr(self, "_M_isauto"): self._setMAuto() E = max(self.M, self.N) fullE = E + 1 == len(self._reorder) == len(self._musUniqueCN) if fullE: mus = self._musUniqueCN[self._reorder] dist = distanceMatrix(mus, magnitude = False)[..., 0] dist[np.arange(E + 1), np.arange(E + 1)] = multifactorial([E]) fitinvE = np.prod(dist, axis = 1) ** -1 vbMng(self, "MAIN", ("Evaluating quasi-Lagrangian basis of degree {} at {} " "sample points.").format(E, E + 1), 5) invD = [np.diag(fitinvE)] else: while E >= 0: if self.polydegreetype == "TOTAL": Eeff = E idxsB = totalDegreeMaxMask(E, self.npar) else: #if self.polydegreetype == "FULL": Eeff = [E] * self.npar idxsB = fullDegreeMaxMask(E, self.npar) TE = pvP(self._musUniqueCN, Eeff, *pvPPar) - fitOut = pseudoInverse(TE, rcond = self.interpRcond, - full = True) + fitOut = pseudoInverse(TE, rcond = self.interpTol, full = True) vbMng(self, "MAIN", ("Fitting {} samples with degree {} through {}... " "Conditioning of pseudoinverse system: {:.4e}.").format( TE.shape[0], E, polyfitname(self.polybasis0), fitOut[1][1][0] / fitOut[1][1][-1]), 5) if fitOut[1][0] == TE.shape[1]: fitinv = fitOut[0][idxsB, :] break if self.catchInstability: raise RROMPyException(("Instability in denominator " "computation: polyfit is poorly " "conditioned."), self.catchInstability == 1) EeqN = E == self.N vbMng(self, "MAIN", ("Polyfit is poorly conditioned. Reducing " "E {} by 1.").format("and N " * EeqN), 10) if EeqN: self.N = self.N - 1 E -= 1 if self.N < 0: raise RROMPyException(("Instability in computation of " "denominator. Aborting.")) invD = vanderInvTable(fitinv, idxsB, self._reorder, self._derIdxs) if self.N == E and not fullE: TN = TE else: #build TN from scratch if self.polydegreetype == "TOTAL": Neff = self.N idxsB = totalDegreeMaxMask(self.N, self.npar) else: #if self.polydegreetype == "FULL": Neff = [self.N] * self.npar idxsB = fullDegreeMaxMask(self.N, self.npar) TN = pvP(self._musUniqueCN, Neff, *pvPPar) return invD, TN def findeveVGExplicit(self, sampleE:sampList, invD:List[Np2D], TN:Np2D, Rscaling : Np1D = None) -> Tuple[Np1D, Np2D]: """ Compute explicitly eigenvalues and eigenvectors of rational denominator matrix. """ RROMPyAssert(self._mode, message = "Cannot solve eigenvalue problem.") vbMng(self, "INIT", "Building gramian matrix.", 10) gramian = self.HFEngine.innerProduct(sampleE, sampleE, is_state = True) if Rscaling is not None: gramian = (gramian.T * Rscaling.conj()).T * Rscaling if self.functionalSolve == "NODAL": SEnd = invD[0].shape[1] G0 = np.zeros((SEnd,) * 2, dtype = np.complex) if self.functionalSolve[:11] == "BARYCENTRIC": G = gramian nEnd = len(gramian) - 1 else: nEnd = TN.shape[1] G = np.zeros((nEnd, nEnd), dtype = np.complex) for k in range(len(invD)): iDkN = dot(invD[k], TN) G += dot(dot(gramian, iDkN).T, iDkN.conj()).T if self.functionalSolve == "NODAL": G0 += dot(dot(gramian, invD[k]).T, invD[k].conj()).T vbMng(self, "DEL", "Done building gramian.", 10) if self.functionalSolve in ["NORM", "BARYCENTRIC_NORM"]: ev, eV = np.linalg.eigh(G) eV = eV[:, 0] if self.functionalSolve == "BARYCENTRIC_NORM": eV = self.findeVBarycentric(eV) problem = "eigenproblem" else: if self.functionalSolve == "BARYCENTRIC_AVERAGE": - fitOut = pseudoInverse(G, rcond = self.interpRcond, - full = True) + fitOut = pseudoInverse(G, rcond = self.interpTol, full = True) eV = self.findeVBarycentric(np.sum(fitOut[0], axis = 1)) else: - fitOut = pseudoInverse(G[:-1, :-1], rcond = self.interpRcond, + fitOut = pseudoInverse(G[:-1, :-1], rcond = self.interpTol, full = True) eV = np.append(fitOut[0].dot(G[:-1, -1]), -1.) ev = fitOut[1][1][::-1] problem = "linear problem" vbMng(self, "MAIN", ("Solved {} of size {} with condition number " "{:.4e}.").format(problem, nEnd, ev[-1] / ev[0]), 5) if self.functionalSolve == "NODAL": eV = self.findeVNodal(eV, G0) return ev, eV def findeveVGQR(self, RPODE:Np2D, invD:List[Np2D], TN:Np2D) -> Tuple[Np1D, Np2D]: """ Compute eigenvalues and eigenvectors of rational denominator matrix through SVD of R factor. """ RROMPyAssert(self._mode, message = "Cannot solve eigenvalue problem.") vbMng(self, "INIT", "Building half-gramian matrix stack.", 10) if self.functionalSolve == "NODAL": gramian = Np2DLikeGramian(None, dot(RPODE, invD[0])) if self.functionalSolve[:11] == "BARYCENTRIC": Rstack = RPODE nEnd = RPODE.shape[1] - 1 else: S, nEnd, eWidth = RPODE.shape[0], TN.shape[1], len(invD) Rstack = np.zeros((S * eWidth, nEnd), dtype = np.complex) for k in range(eWidth): Rstack[k * S : (k + 1) * S, :] = dot(RPODE, dot(invD[k], TN)) vbMng(self, "DEL", "Done building half-gramian.", 10) if self.functionalSolve in ["NORM", "BARYCENTRIC_NORM", "BARYCENTRIC_AVERAGE"]: _, ev, Vh = np.linalg.svd(Rstack, full_matrices = False) if self.functionalSolve in ["NORM", "BARYCENTRIC_NORM"]: eV = Vh[-1, :].conj() ev = ev[::-1] else: #if self.functionalSolve == "BARYCENTRIC_AVERAGE": ev[np.logical_not(np.isclose(ev, 0.))] **= -2. eV = Vh.T.conj().dot(ev * np.sum(Vh, axis = 1)) if self.functionalSolve[:11] == "BARYCENTRIC": eV = self.findeVBarycentric(eV) problem = "svd problem" else: - fitOut = pseudoInverse(Rstack[:, :-1], rcond = self.interpRcond, + fitOut = pseudoInverse(Rstack[:, :-1], rcond = self.interpTol, full = True) ev = fitOut[1][1][::-1] eV = np.append(fitOut[0].dot(Rstack[:, -1]), -1.) problem = "linear problem" vbMng(self, "MAIN", ("Solved {} of size {} with condition number " "{:.4e}.").format(problem, nEnd, ev[-1] / ev[0]), 5) if self.functionalSolve == "NODAL": eV = self.findeVNodal(eV, gramian) return ev, eV def findeVBarycentric(self, baryWeights:Np1D) -> Np1D: RROMPyAssert(self._mode, message = "Cannot solve optimization problem.") arrow = np.pad(np.diag(self._musUniqueCN[self._reorder].flatten()), (1, 0), "constant", constant_values = 1.) + 0.j arrow[0, 0] = 0. arrow[0, 1:] = baryWeights active = np.pad(np.eye(len(baryWeights)), (1, 0), "constant") eV = eigvals(arrow, active) return eV[np.logical_not(np.isinf(eV))] def findeVNodal(self, q0coeffs:Np1D, gram:Np2D, maxiter : int = 25, tolNewton : float = 1e-10) -> Np1D: RROMPyAssert(self._mode, message = "Cannot solve optimization problem.") q = PI() q.npar, q.polybasis, q.coeffs = self.npar, self.polybasis0, q0coeffs nodes = q.roots() N = len(nodes) grad = np.zeros(N, dtype = np.complex) hess = np.zeros((N, N), dtype = np.complex) for niter in range(maxiter): plDist = distanceMatrix(self._musUniqueCN[self._reorder], nodes, magnitude = False)[:, :, 0] q0, q = np.prod(plDist, axis = 1), [] for jS in range(N): mask = np.arange(N) != jS q += [np.prod(plDist[:, mask], axis = 1)] grad[jS] = q[-1].conj().dot(gram.dot(q0)) for iS in range(jS + 1): if iS == jS: hij = 0. sij = q[-1].conj().dot(gram.dot(q[-1])) else: mask = np.logical_and(np.arange(N) != iS, np.arange(N) != jS) qij = np.prod(plDist[:, mask], axis = 1) hij = qij.conj().dot(gram.dot(q0)) sij = q[-1].conj().dot(gram.dot(q[iS])) hess[jS, iS] = hij + sij if iS < jS: hess[iS, jS] = hij + np.conj(sij) dnodes = np.linalg.solve(hess, grad) nodes += dnodes tol = np.linalg.norm(dnodes) / np.linalg.norm(nodes) if tol < tolNewton: break if niter < maxiter: vbMng(self, "MAIN", ("Newton's method for problem of size {} converged " "(err = {:.4e}) in {} iteration{}.").format( N + 1, tol, niter + 1, "s" * (niter > 0)), 5) else: RROMPyWarning(("Newton's method for problem of size {} did " "not converge (err = {:.4e}) after {} " "iterations.").format(N + 1, tol, niter + 1)) return nodes def getResidues(self, *args, **kwargs) -> Np2D: """ Obtain approximant residues. Returns: Matrix with residues as columns. """ return self.trainedModel.getResidues(*args, **kwargs) diff --git a/rrompy/reduction_methods/standard/trained_model/trained_model_rational.py b/rrompy/reduction_methods/standard/trained_model/trained_model_rational.py index 543a05d..eb055c8 100644 --- a/rrompy/reduction_methods/standard/trained_model/trained_model_rational.py +++ b/rrompy/reduction_methods/standard/trained_model/trained_model_rational.py @@ -1,195 +1,197 @@ # Copyright (C) 2018-2020 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 collections.abc import Iterable from rrompy.reduction_methods.base.trained_model.trained_model import ( TrainedModel) from rrompy.utilities.numerical import dot from rrompy.utilities.numerical.compress_matrix import compressMatrix from rrompy.utilities.base.types import (Np1D, Np2D, List, paramVal, paramList, sampList) from rrompy.utilities.base import verbosityManager as vbMng, freepar as fp from rrompy.utilities.exception_manager import RROMPyException, RROMPyWarning from rrompy.parameter import emptyParameterList from rrompy.sampling import sampleList __all__ = ['TrainedModelRational'] class TrainedModelRational(TrainedModel): """ ROM approximant evaluation for rational approximant. Attributes: Data: dictionary with all that can be pickled. """ def compress(self, collapse : bool = False, tol : float = 0., *args, **kwargs): if not collapse and tol <= 0.: return if hasattr(self.data, "_is_C_quadratic") and self.data._is_C_quadratic: raise RROMPyException(("Cannot compress model with quadratic " "output.")) RMat = self.data.projMat if not collapse: if hasattr(self.data, "_compressTol"): RROMPyWarning(("Recompressing already compressed model is " "ineffective. Aborting.")) return self.data.projMat, RMat, _ = compressMatrix(RMat, tol, *args, **kwargs) self.data.P.postmultiplyTensorize(RMat.T) super().compress(collapse, tol) def centerNormalize(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 = self.checkParameterList(mu) if mu0 is None: mu0 = self.data.mu0 return (self.mapParameterList(mu) - self.mapParameterList(mu0)) / self.data.scaleFactor def getPVal(self, mu : paramList = []) -> sampList: """ Evaluate rational numerator at arbitrary parameter. Args: mu: Target parameter. """ mu = self.checkParameterList(mu) vbMng(self, "INIT", "Evaluating numerator at mu = {}.".format(mu), 17) p = sampleList(self.data.P(self.centerNormalize(mu))) vbMng(self, "DEL", "Done evaluating numerator.", 17) 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. """ mu = self.checkParameterList(mu) vbMng(self, "INIT", "Evaluating denominator at mu = {}.".format(mu), 17) q = self.data.Q(self.centerNormalize(mu), der, scl) vbMng(self, "DEL", "Done evaluating denominator.", 17) return q def getApproxReduced(self, mu : paramList = []) -> sampList: """ Evaluate reduced representation of approximant at arbitrary parameter. Args: mu: Target parameter. """ mu = self.checkParameterList(mu) if (not hasattr(self, "lastSolvedApproxReduced") or self.lastSolvedApproxReduced != mu): vbMng(self, "INIT", "Evaluating approximant at mu = {}.".format(mu), 12) QV = self.getQVal(mu) QVzero = np.where(QV == 0.)[0] if len(QVzero) > 0: QV[QVzero] = np.finfo(np.complex).eps / (1. + self.data.Q.deg[0]) self.uApproxReduced = self.getPVal(mu) / QV 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. """ 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 isinstance(mVals, Iterable): mVals = [mVals] mVals = list(mVals) else: mVals = [fp] rDim = mVals.index(fp) if rDim < len(mVals) - 1 and fp in mVals[rDim + 1 :]: raise RROMPyException(("Exactly 1 'freepar' entry in " "marginalVals must be provided.")) mVals[rDim] = self.data.mu0(rDim) mVals = list(self.centerNormalize(mVals).data.flatten()) mVals[rDim] = fp roots = self.data.scaleFactor[rDim] * self.data.Q.roots(mVals) return self.mapParameterList(self.mapParameterList(self.data.mu0(rDim), idx = [rDim])(0, 0) + roots, "B", [rDim])(0) def getResidues(self, *args, **kwargs) -> Np2D: """ Obtain approximant residues. Returns: Numpy matrix with residues as columns. """ pls = self.getPoles(*args, **kwargs) + if len(pls) == 0: + return pls, np.empty((0, 0), dtype = self.data.P.coeffs.dtype) if len(args) == 1: mVals = args[0] elif len(args) == 0: mVals = [None] else: mVals = kwargs["marginalVals"] if not isinstance(mVals, Iterable): mVals = [mVals] mVals = list(mVals) rDim = mVals.index(fp) poles = emptyParameterList() poles.reset((len(pls), self.data.npar), dtype = pls.dtype) for k, pl in enumerate(pls): mValsLoc = list(mVals) mValsLoc[rDim] = pl poles[k] = mValsLoc QV = self.getQVal(poles, list(1 * (np.arange(self.data.npar) == rDim))) QVzero = np.where(QV == 0.)[0] if len(QVzero) > 0: RROMPyWarning(("Adjusting residuals to avoid division by " "numerically zero denominator.")) QV[QVzero] = np.finfo(np.complex).eps / (1. + self.data.Q.deg[0]) res = self.getPVal(poles).data / QV if hasattr(self.data, "_is_C_quadratic") and self.data._is_C_quadratic: self._setupQuadMapping() res = res[self.quad_mapping[0]] * res[self.quad_mapping[1]].conj() if not self.data._collapsed: res = dot(self.data.projMat, res).T if (hasattr(self.data, "_is_C_quadratic") and self.data._is_C_quadratic == 2): res = np.real(res) return pls, res diff --git a/tests/3_reduction_methods_1D/rational_interpolant_1d.py b/tests/3_reduction_methods_1D/rational_interpolant_1d.py index c747b0b..32711b5 100644 --- a/tests/3_reduction_methods_1D/rational_interpolant_1d.py +++ b/tests/3_reduction_methods_1D/rational_interpolant_1d.py @@ -1,67 +1,67 @@ # Copyright (C) 2018-2020 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 import RationalInterpolant as RI from rrompy.parameter.parameter_sampling import (QuadratureSampler as QS, ManualSampler as MS) from rrompy.parameter import checkParameterList def test_monomials(capsys): mu = 1.5 solver = matrixFFT() - params = {"POD": False, "S": 10, "robustTol": 1e-6, "interpRcond": 1e-3, + params = {"POD": False, "S": 10, "QTol": 1e-6, "interpTol": 1e-3, "polybasis": "MONOMIAL", "sampler": QS([1.5, 6.5], "UNIFORM")} approx = RI(solver, 4., approxParameters = params, verbosity = 10) approx.setupApprox() out, err = capsys.readouterr() assert "below tolerance. Reducing N " in out assert "poorly conditioned. Reducing M " in out assert len(err) == 0 assert np.isclose(approx.normErr(mu)[0], 1.4746e-05, atol = 1e-4) def test_well_cond(): mu = 1.5 solver = matrixFFT() - params = {"POD": True, "S": 10, "robustTol": 1e-14, "interpRcond": 1e-10, + params = {"POD": True, "S": 10, "QTol": 1e-14, "interpTol": 1e-10, "polybasis": "CHEBYSHEV", "sampler": QS([1., 7.], "CHEBYSHEV")} approx = RI(solver, 4., approxParameters = params, verbosity = 0) approx.setupApprox() poles = approx.getPoles() for lambda_ in np.arange(1, 8): assert np.isclose(np.min(np.abs(poles - lambda_)), 0., atol = 1e-4) for mu in approx.mus: assert np.isclose(approx.normErr(mu)[0], 0., atol = 1e-8) 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, "polybasis": "CHEBYSHEV", "sampler": MS([1., 7.], points = points)} approx = RI(solver, 4., approxParameters = params, verbosity = 0) approx.setupApprox() poles = approx.getPoles() for lambda_ in np.arange(1, 8): assert np.isclose(np.min(np.abs(poles - lambda_)), 0., atol = 1e-4) for mu in approx.mus: assert np.isclose(approx.normErr(mu)[0], 0., atol = 1e-8) diff --git a/tests/4_reduction_methods_multiD/pivoted_rational_2d.py b/tests/4_reduction_methods_multiD/pivoted_rational_2d.py index 15c2c35..2a65e0a 100644 --- a/tests/4_reduction_methods_multiD/pivoted_rational_2d.py +++ b/tests/4_reduction_methods_multiD/pivoted_rational_2d.py @@ -1,112 +1,112 @@ # Copyright (C) 2018-2020 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 import ( RationalInterpolantPivotedPoleMatch as RIP, RationalInterpolantGreedyPivotedPoleMatch as RIGP) from rrompy.parameter.parameter_sampling import (QuadratureSampler as QS, ManualSampler as MS) def test_pivoted_uniform(): mu = [5.05, 7.1] mu0 = [5., 7.] solver = matrixRandom() uh = solver.solve(mu)[0] params = {"POD": True, "S": 5, "polybasis": "CHEBYSHEV", "samplerPivot": QS([4.75, 5.25], "CHEBYSHEV"), "SMarginal": 5, "polybasisMarginal": "MONOMIAL", "matchingWeight": 1., "samplerMarginal": QS([6.75, 7.25], "UNIFORM")} approx = RIP([0], solver, mu0, approxParameters = params, verbosity = 0) approx.setupApprox() uhP1 = approx.getApprox(mu)[0] errP = approx.getErr(mu)[0] errNP = approx.normErr(mu)[0] myerrP = uhP1 - uh assert np.allclose(np.abs(errP - myerrP), 0., rtol = 1e-3) assert np.isclose(solver.norm(errP), errNP, rtol = 1e-3) resP = approx.getRes(mu)[0] resNP = approx.normRes(mu) assert np.isclose(solver.norm(resP), resNP, rtol = 1e-3) assert np.allclose(np.abs(resP - (solver.b(mu) - solver.A(mu).dot(uhP1))), 0., rtol = 1e-3) assert np.isclose(errNP / solver.norm(uh), 6.0631706e-04, rtol = 1) def test_pivoted_manual_grid(capsys): mu = [5.05, 7.1] mu0 = [5., 7.] solver = matrixRandom() uh = solver.solve(mu)[0] params = {"POD": False, "S": 5, "polybasis": "MONOMIAL", "samplerPivot": MS([4.75, 5.25], np.array([5.]), normalFoci = [0., 0.]), "SMarginal": 5, "polybasisMarginal": "MONOMIAL", "matchingWeight": 1., "samplerMarginal": MS([6.75, 7.25], np.linspace(6.75, 7.25, 5)), - "robustTol": 1e-6, "interpRcond": 1e-3} + "QTol": 1e-6, "interpTol": 1e-3} approx = RIP([0], solver, mu0, approxParameters = params, verbosity = 0) approx.setupApprox() uhP1 = approx.getApprox(mu)[0] errP = approx.getErr(mu)[0] errNP = approx.normErr(mu)[0] myerrP = uhP1 - uh assert np.allclose(np.abs(errP - myerrP), 0., rtol = 1e-3) assert np.isclose(solver.norm(errP), errNP, rtol = 1e-3) resP = approx.getRes(mu)[0] resNP = approx.normRes(mu) assert np.isclose(solver.norm(resP), resNP, rtol = 1e-3) assert np.allclose(np.abs(resP - (solver.b(mu) - solver.A(mu).dot(uhP1))), 0., rtol = 1e-3) assert np.isclose(errNP / solver.norm(uh), .4763489, rtol = 1) out, err = capsys.readouterr() assert ("poorly conditioned" not in out) assert len(err) == 0 def test_pivoted_greedy(): mu = [5.05, 7.1] mu0 = [5., 7.] solver = matrixRandom() uh = solver.solve(mu)[0] params = {"POD": True, "nTestPoints": 100, "greedyTol": 1e-4, "collinearityTol": 1e8, "errorEstimatorKind": "DISCREPANCY", "S": 5, "polybasis": "CHEBYSHEV", "samplerPivot": QS([4.75, 5.25], "UNIFORM"), "trainSetGenerator": QS([4.75, 5.25], "CHEBYSHEV"), "SMarginal": 5, "polybasisMarginal": "MONOMIAL", "samplerMarginal": QS([6.75, 7.25], "UNIFORM"), "matchingWeight": 1.} solver.cutOffPolesRMinRel, solver.cutOffPolesRMaxRel = -3., 3. solver.cutOffPolesIMinRel, solver.cutOffPolesIMaxRel = -1.5, 1.5 approx = RIGP([0], solver, mu0, approxParameters = params, verbosity = 0) approx.setupApprox() uhP1 = approx.getApprox(mu)[0] errP = approx.getErr(mu)[0] errNP = approx.normErr(mu)[0] myerrP = uhP1 - uh assert np.allclose(np.abs(errP - myerrP), 0., rtol = 1e-3) assert np.isclose(solver.norm(errP), errNP, rtol = 1e-3) resP = approx.getRes(mu)[0] resNP = approx.normRes(mu) assert np.isclose(solver.norm(resP), resNP, rtol = 1e-3) assert np.allclose(np.abs(resP - (solver.b(mu) - solver.A(mu).dot(uhP1))), 0., rtol = 1e-3) assert np.isclose(errNP / solver.norm(uh), 1.181958e-02, rtol = 1) diff --git a/tests/4_reduction_methods_multiD/rational_interpolant_2d.py b/tests/4_reduction_methods_multiD/rational_interpolant_2d.py index 2f8242f..5e2cb73 100644 --- a/tests/4_reduction_methods_multiD/rational_interpolant_2d.py +++ b/tests/4_reduction_methods_multiD/rational_interpolant_2d.py @@ -1,77 +1,77 @@ # Copyright (C) 2018-2020 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 import RationalInterpolant as RI from rrompy.parameter.parameter_sampling import (QuadratureSampler as QS, ManualSampler as MS) def test_monomials(capsys): mu = [5.05, 7.1] mu0 = [5., 7.] solver = matrixRandom() uh = solver.solve(mu)[0] - params = {"POD": False, "S": 16, "robustTol": 1e-6, - "interpRcond": 1e-3, "polybasis": "MONOMIAL", + params = {"POD": False, "S": 16, "QTol": 1e-6, "interpTol": 1e-3, + "polybasis": "MONOMIAL", "sampler": QS([[4.9, 6.85], [5.1, 7.15]], "UNIFORM")} approx = RI(solver, mu0, params, verbosity = 100) approx.setupApprox() uhP1 = approx.getApprox(mu)[0] errP = approx.getErr(mu)[0] errNP = approx.normErr(mu)[0] myerrP = uhP1 - uh assert np.allclose(np.abs(errP - myerrP), 0., rtol = 1e-3) assert np.isclose(solver.norm(errP), errNP, rtol = 1e-3) resP = approx.getRes(mu)[0] resNP = approx.normRes(mu) assert np.isclose(solver.norm(resP), resNP, rtol = 1e-3) assert np.allclose(np.abs(resP - (solver.b(mu) - solver.A(mu).dot(uhP1))), 0., rtol = 1e-3) assert np.isclose(errNP / solver.norm(uh), 5.2667e-05, rtol = 1) out, err = capsys.readouterr() assert ("poorly conditioned. Reducing E " in out) assert len(err) == 0 def test_well_cond(): mu = [5.05, 7.1] mu0 = [5., 7.] solver = matrixRandom() params = {"POD": True, "M": 3, "N": 3, "S": 16, - "interpRcond": 1e-10, "polybasis": "CHEBYSHEV", + "interpTol": 1e-10, "polybasis": "CHEBYSHEV", "sampler": QS([[4.9, 6.85], [5.1, 7.15]], "UNIFORM")} approx = RI(solver, mu0, params, verbosity = 0) approx.setupApprox() assert np.isclose(approx.normErr(mu)[0] / approx.normHF(mu)[0], 5.98695e-05, rtol = 1e-1) def test_hermite(): mu = [5.05, 7.1] mu0 = [5., 7.] solver = matrixRandom() sampler0 = QS([[4.9, 6.85], [5.1, 7.15]], "UNIFORM") params = {"POD": True, "M": 3, "N": 3, "S": 25, "polybasis": "CHEBYSHEV", "sampler": MS([[4.9, 6.85], [5.1, 7.15]], points = sampler0.generatePoints(9))} approx = RI(solver, mu0, params, verbosity = 0) approx.setupApprox() assert np.isclose(approx.normErr(mu)[0] / approx.normHF(mu)[0], 5.50053e-05, rtol = 5e-1)