diff --git a/rrompy/reduction_methods/distributed_greedy/generic_distributed_greedy_approximant.py b/rrompy/reduction_methods/distributed_greedy/generic_distributed_greedy_approximant.py index 389621b..0d938bd 100644 --- a/rrompy/reduction_methods/distributed_greedy/generic_distributed_greedy_approximant.py +++ b/rrompy/reduction_methods/distributed_greedy/generic_distributed_greedy_approximant.py @@ -1,657 +1,651 @@ # Copyright (C) 2018 by the RROMPy authors # # This file is part of RROMPy. # # RROMPy is free software: you can redistribute it and/or modify # it under the terms of the GNU Lesser General Public License as published by # the Free Software Foundation, either version 3 of the License, or # (at your option) any later version. # # RROMPy is distributed in the hope that it will be useful, # but WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # GNU Lesser General Public License for more details. # # You should have received a copy of the GNU Lesser General Public License # along with RROMPy. If not, see . # from copy import deepcopy as copy import numpy as np from rrompy.reduction_methods.distributed.generic_distributed_approximant \ import GenericDistributedApproximant from rrompy.utilities.base.types import Np1D, Np2D, DictAny, HFEng, Tuple, List from rrompy.utilities.base import purgeDict, verbosityDepth from rrompy.utilities.exception_manager import (RROMPyException, modeAssert, RROMPyWarning) __all__ = ['GenericDistributedGreedyApproximant'] class normEngine: def __init__(self, energyNormMatrix:Np2D): self.energyNormMatrix = copy(energyNormMatrix) def innerProduct(self, u:Np2D, v:Np2D, onlyDiag : bool = False) -> Np2D: if onlyDiag: return np.sum(self.energyNormMatrix.dot(u) * v.conj(), axis = 0) return v.T.conj().dot(self.energyNormMatrix.dot(u)) def norm(self, u:Np2D) -> Np1D: return np.abs(self.innerProduct(u, u, onlyDiag = True)) ** .5 +def pruneSamples(mus:Np1D, badmus:Np1D, tol : float = 1e-8) -> Np1D: + """Remove from mus all the elements which are too close to badmus.""" + proximity = np.min(np.abs(mus.reshape(-1, 1) + - np.tile(badmus.reshape(1, -1), [len(mus), 1])), + axis = 1).flatten() + proxMask = proximity > tol + return mus[proxMask] + + class GenericDistributedGreedyApproximant(GenericDistributedApproximant): """ ROM greedy interpolant computation for parametric problems (ABSTRACT). Args: HFEngine: HF problem solver. mu0(optional): Default parameter. Defaults to 0. approxParameters(optional): Dictionary containing values for main parameters of approximant. Recognized keys are: - 'POD': whether to compute POD of snapshots; defaults to True; - 'muBounds': list of bounds for parameter values; defaults to [0, 1]; - 'S': number of starting training points; defaults to 2; - 'sampler': sample point generator; defaults to uniform sampler on muBounds; - 'greedyTol': uniform error tolerance for greedy algorithm; defaults to 1e-2; - 'interactive': whether to interactively terminate greedy algorithm; defaults to False; - 'maxIter': maximum number of greedy steps; defaults to 1e2; - 'refinementRatio': ratio of test points to be exhausted before test set refinement; defaults to 0.2; - 'nTestPoints': number of test points; defaults to maxIter / refinementRatio; - 'trainSetGenerator': training sample points generator; defaults to Chebyshev sampler within muBounds. Defaults to empty dict. homogeneized(optional): Whether to homogeneize Dirichlet BCs. Defaults to False. verbosity(optional): Verbosity level. Defaults to 10. Attributes: HFEngine: HF problem solver. mu0: Default parameter. mus: Array of snapshot parameters. homogeneized: Whether to homogeneize Dirichlet BCs. approxParameters: Dictionary containing values for main parameters of approximant. Recognized keys are in parameterList. parameterList: Recognized keys of approximant parameters: - 'POD': whether to compute POD of snapshots; - 'muBounds': list of bounds for parameter values; - 'S': number of starting training points; - 'sampler': sample point generator; - 'greedyTol': uniform error tolerance for greedy algorithm; - 'interactive': whether to interactively terminate greedy algorithm; - 'maxIter': maximum number of greedy steps; - 'refinementRatio': ratio of test points to be exhausted before test set refinement; - 'nTestPoints': number of test points; - 'trainSetGenerator': training sample points generator. extraApproxParameters: List of approxParameters keys in addition to mother class's. POD: whether to compute POD of snapshots. muBounds: list of bounds for parameter values. S: number of test points. sampler: Sample point generator. greedyTol: uniform error tolerance for greedy algorithm. maxIter: maximum number of greedy steps. refinementRatio: ratio of training points to be exhausted before training set refinement. nTestPoints: number of starting training points. trainSetGenerator: training sample points generator. robustTol: tolerance for robust Pade' denominator management. samplingEngine: Sampling engine. estimatorNormEngine: Engine for estimator norm computation. uHF: High fidelity solution with wavenumber lastSolvedHF as numpy complex vector. lastSolvedHF: Wavenumber corresponding to last computed high fidelity solution. uApp: Last evaluated approximant as numpy complex vector. """ TOL_INSTABILITY = 1e-6 def __init__(self, HFEngine:HFEng, mu0 : complex = 0., approxParameters : DictAny = {}, homogeneized : bool = False, verbosity : int = 10, timestamp : bool = True): self._preInit() self._addParametersToList(["greedyTol", "interactive", "maxIter", "refinementRatio", "nTestPoints", "trainSetGenerator"]) super().__init__(HFEngine = HFEngine, mu0 = mu0, approxParameters = approxParameters, homogeneized = homogeneized, verbosity = verbosity, timestamp = timestamp) self._postInit() @property def approxParameters(self): """Value of approximant parameters. Its assignment may change S.""" return self._approxParameters @approxParameters.setter def approxParameters(self, approxParams): approxParameters = purgeDict(approxParams, self.parameterList, dictname = self.name() + ".approxParameters", baselevel = 1) approxParametersCopy = purgeDict(approxParameters, ["greedyTol", "interactive", "maxIter", "refinementRatio", "nTestPoints", "trainSetGenerator"], True, True, baselevel = 1) GenericDistributedApproximant.approxParameters.fset(self, approxParametersCopy) keyList = list(approxParameters.keys()) if "greedyTol" in keyList: self.greedyTol = approxParameters["greedyTol"] elif not hasattr(self, "_greedyTol") or self.greedyTol is None: self.greedyTol = 1e-2 if "interactive" in keyList: self.interactive = approxParameters["interactive"] elif not hasattr(self, "interactive") or self.interactive is None: self.interactive = False if "maxIter" in keyList: self.maxIter = approxParameters["maxIter"] elif not hasattr(self, "_maxIter") or self.maxIter is None: self.maxIter = 1e2 if "refinementRatio" in keyList: self.refinementRatio = approxParameters["refinementRatio"] elif (not hasattr(self, "_refinementRatio") or self.refinementRatio is None): self.refinementRatio = 0.2 if "nTestPoints" in keyList: self.nTestPoints = approxParameters["nTestPoints"] elif (not hasattr(self, "_nTestPoints") or self.nTestPoints is None): self.nTestPoints = np.int(np.ceil(self.maxIter / self.refinementRatio)) if "trainSetGenerator" in keyList: self.trainSetGenerator = approxParameters["trainSetGenerator"] elif (not hasattr(self, "_trainSetGenerator") or self.trainSetGenerator is None): from rrompy.utilities.parameter_sampling import QuadratureSampler self.trainSetGenerator = QuadratureSampler(self.muBounds, "CHEBYSHEV") del QuadratureSampler @property def mus(self): """Value of mus.""" return self._mus @mus.setter def mus(self, mus): self._mus = np.array(mus, dtype = np.complex) @property def greedyTol(self): """Value of greedyTol.""" return self._greedyTol @greedyTol.setter def greedyTol(self, greedyTol): if greedyTol < 0: raise RROMPyException("greedyTol must be non-negative.") if hasattr(self, "_greedyTol") and self.greedyTol is not None: greedyTolold = self.greedyTol else: greedyTolold = -1 self._greedyTol = greedyTol self._approxParameters["greedyTol"] = self.greedyTol if greedyTolold != self.greedyTol: self.resetSamples() @property def maxIter(self): """Value of maxIter.""" return self._maxIter @maxIter.setter def maxIter(self, maxIter): if maxIter <= 0: raise RROMPyException("maxIter must be positive.") if hasattr(self, "_maxIter") and self.maxIter is not None: maxIterold = self.maxIter else: maxIterold = -1 self._maxIter = maxIter self._approxParameters["maxIter"] = self.maxIter if maxIterold != self.maxIter: self.resetSamples() @property def refinementRatio(self): """Value of refinementRatio.""" return self._refinementRatio @refinementRatio.setter def refinementRatio(self, refinementRatio): if refinementRatio <= 0. or refinementRatio > 1.: raise RROMPyException(("refinementRatio must be between 0 " "(excluded) and 1.")) if (hasattr(self, "_refinementRatio") and self.refinementRatio is not None): refinementRatioold = self.refinementRatio else: refinementRatioold = -1 self._refinementRatio = refinementRatio self._approxParameters["refinementRatio"] = self.refinementRatio if refinementRatioold != self.refinementRatio: self.resetSamples() @property def nTestPoints(self): """Value of nTestPoints.""" return self._nTestPoints @nTestPoints.setter def nTestPoints(self, nTestPoints): if nTestPoints <= 0: raise RROMPyException("nTestPoints must be positive.") if not np.isclose(nTestPoints, np.int(nTestPoints)): raise RROMPyException("nTestPoints must be an integer.") nTestPoints = np.int(nTestPoints) if hasattr(self, "_nTestPoints") and self.nTestPoints is not None: nTestPointsold = self.nTestPoints else: nTestPointsold = -1 self._nTestPoints = nTestPoints self._approxParameters["nTestPoints"] = self.nTestPoints if nTestPointsold != self.nTestPoints: self.resetSamples() @property def trainSetGenerator(self): """Value of trainSetGenerator.""" return self._trainSetGenerator @trainSetGenerator.setter def trainSetGenerator(self, trainSetGenerator): if 'generatePoints' not in dir(trainSetGenerator): raise RROMPyException("trainSetGenerator type not recognized.") if (hasattr(self, '_trainSetGenerator') and self.trainSetGenerator is not None): trainSetGeneratorOld = self.trainSetGenerator self._trainSetGenerator = trainSetGenerator self._approxParameters["trainSetGenerator"] = self.trainSetGenerator if (not 'trainSetGeneratorOld' in locals() or trainSetGeneratorOld != self.trainSetGenerator): self.resetSamples() - + def resetSamples(self): """Reset samples.""" super().resetSamples() self._mus = [] def initEstimatorNormEngine(self, normEng = None): """Initialize estimator norm engine.""" if (normEng is not None or not hasattr(self, "estimatorNormEngine") or self.estimatorNormEngine is None): if normEng is None: if not hasattr(self.HFEngine, "energyNormMatrix"): self.HFEngine.buildEnergyNormForm() estimatorEnergyMatrix = self.HFEngine.energyNormMatrix else: if hasattr(normEng, "buildEnergyNormForm"): if not hasattr(normEng, "energyNormMatrix"): normEng.buildEnergyNormForm() estimatorEnergyMatrix = normEng.energyNormMatrix else: estimatorEnergyMatrix = normEng self.estimatorNormEngine = normEngine(estimatorEnergyMatrix) def errorEstimator(self, mus:List[np.complex]) -> List[np.complex]: """ Standard residual-based error estimator with explicit residual computation. """ self.setupApprox() nmus = len(mus) err = np.empty(nmus) if self.HFEngine.nbs == 1: RHS = self.getRHS(mus[0], homogeneized = self.homogeneized) RHSNorm = self.estimatorNormEngine.norm(RHS) for j in range(nmus): res = self.getRes(mus[j], homogeneized = self.homogeneized) err[j] = self.estimatorNormEngine.norm(res) / RHSNorm else: for j in range(nmus): res = self.getRes(mus[j], homogeneized = self.homogeneized) RHS = self.getRHS(mus[j], homogeneized = self.homogeneized) err[j] = (self.estimatorNormEngine.norm(res) / self.estimatorNormEngine.norm(RHS)) return np.abs(err) def getMaxErrorEstimator(self, mus:List[np.complex], plot : bool = False) -> Tuple[Np1D, int, float]: """ Compute maximum of (and index of maximum of) error estimator over given parameters. """ errorEstTest = self.errorEstimator(mus) idxMaxEst = np.argmax(errorEstTest) maxEst = errorEstTest[idxMaxEst] if plot and not np.all(np.isinf(errorEstTest)): from matplotlib import pyplot as plt plt.figure() plt.semilogy(np.real(mus), errorEstTest, 'k') plt.semilogy(np.real(mus[[0, -1]]), [self.greedyTol] * 2, 'r--') plt.semilogy(np.real(self.mus), 2. * self.greedyTol * np.ones(len(self.mus)), '*m') plt.semilogy(np.real(mus[idxMaxEst]), maxEst, 'xr') plt.grid() plt.show() plt.close() return errorEstTest, idxMaxEst, maxEst def greedyNextSample(self, muidx:int, plotEst : bool = False)\ -> Tuple[Np1D, int, float, complex]: """Compute next greedy snapshot of solution map.""" modeAssert(self._mode, message = "Cannot add greedy sample.") mu = self.muTest[muidx] + self.muTest = np.delete(self.muTest, muidx) if self.verbosity >= 2: verbosityDepth("MAIN", ("Adding {}-th sample point at {} to " "training set.").format( self.samplingEngine.nsamples + 1, mu), timestamp = self.timestamp) self.mus = np.append(self.mus, mu) - idxs = np.arange(len(self.muTest)) - mask = np.ones_like(idxs, dtype = bool) - mask[muidx] = False - idxs = idxs[mask] - self.muTest = self.muTest[idxs] self.samplingEngine.nextSample(mu, homogeneized = self.homogeneized) errorEstTest, muidx, maxErrorEst = self.getMaxErrorEstimator( self.muTest, plotEst) return errorEstTest, muidx, maxErrorEst, self.muTest[muidx] def _preliminaryTraining(self): """Initialize starting snapshots of solution map.""" modeAssert(self._mode, message = "Cannot start greedy algorithm.") + self.computeScaleFactor() if self.samplingEngine.samples is not None: return if self.verbosity >= 2: verbosityDepth("INIT", "Starting computation of snapshots.", timestamp = self.timestamp) self.resetSamples() self.mus, _ = self.trainSetGenerator.generatePoints(self.S) muTestBase, _ = self.sampler.generatePoints(self.nTestPoints) - proxVal = np.min(np.abs(muTestBase.reshape(-1, 1) - - np.tile(self.mus.reshape(1, -1), - [self.nTestPoints, 1])), - axis = 1) - proxMask = ~(proxVal < 1e-12 * np.abs(muTestBase[0] - muTestBase[-1])) - self.muTest = np.empty(np.sum(proxMask) + 1, dtype = np.complex) - self.muTest[:-1] = np.sort(muTestBase[proxMask]).flatten() - self.muTest[-1] = self.mus[-1] - self.mus = self.mus[:-1] - for j in range(len(self.mus)): + for j in range(len(self.mus) - 1): if self.verbosity >= 2: verbosityDepth("MAIN", ("Adding {}-th sample point at {} to training " "set.").format(self.samplingEngine.nsamples+ 1, self.mus[j]), timestamp = self.timestamp) self.samplingEngine.nextSample(self.mus[j], homogeneized = self.homogeneized) + muTestBase = np.sort(pruneSamples(muTestBase, self.mus[: -1], + 1e-10 * self.scaleFactor)) + self.muTest = np.empty(len(muTestBase) + 1, dtype = np.complex) + self.muTest[: -1] = muTestBase + self.muTest[-1] = self.mus[-1] + self.mus = self.mus[: -1] def _enrichTestSet(self, nTest:int): - """Double number of elements of test set.""" + """Add extra elements to test set.""" muTestExtra, _ = self.sampler.generatePoints(2 * nTest) - muGiven = np.append(self.mus, self.muTest).reshape(1, -1) - proxVal = np.min(np.abs(muTestExtra.reshape(-1, 1) - - np.tile(muGiven, [2 * nTest, 1])), - axis = 1) - proxMask = ~(proxVal < 1e-12 * np.abs(muTestExtra[0]-muTestExtra[-1])) - muTestNew = np.empty(len(self.muTest) + np.sum(proxMask), + muTestExtra = pruneSamples(muTestExtra, + np.append(self.mus, self.muTest), + 1e-10 * self.scaleFactor) + muTestNew = np.empty(len(self.muTest) + len(muTestExtra), dtype = np.complex) muTestNew[: len(self.muTest)] = self.muTest - muTestNew[len(self.muTest) :] = muTestExtra[proxMask] + muTestNew[len(self.muTest) :] = muTestExtra self.muTest = np.sort(muTestNew) if self.verbosity >= 5: verbosityDepth("MAIN", "Enriching test set by {} elements.".format( - np.sum(proxMask)), + len(muTestExtra)), timestamp = self.timestamp) def greedy(self, plotEst : bool = False): """Compute greedy snapshots of solution map.""" modeAssert(self._mode, message = "Cannot start greedy algorithm.") if self.samplingEngine.samples is not None: return self._preliminaryTraining() nTest = self.nTestPoints errorEstTest, muidx, maxErrorEst, mu = self.greedyNextSample(-1, plotEst) if self.verbosity >= 2: verbosityDepth("MAIN", ("Uniform testing error estimate " "{:.4e}.").format(maxErrorEst), timestamp = self.timestamp) trainedModelOld = copy(self.trainedModel) while (self.samplingEngine.nsamples < self.maxIter and maxErrorEst > self.greedyTol): if (1. - self.refinementRatio) * nTest > len(self.muTest): self._enrichTestSet(nTest) nTest = len(self.muTest) muTestOld, maxErrorEstOld = self.muTest, maxErrorEst errorEstTest, muidx, maxErrorEst, mu = self.greedyNextSample( muidx, plotEst) if self.verbosity >= 2: verbosityDepth("MAIN", ("Uniform testing error estimate " "{:.4e}.").format(maxErrorEst), timestamp = self.timestamp) if (np.isnan(maxErrorEst) or np.isinf(maxErrorEst) or maxErrorEstOld < maxErrorEst * self.TOL_INSTABILITY): RROMPyWarning(("Instability in a posteriori estimator. " "Starting preemptive greedy loop termination.")) maxErrorEst = maxErrorEstOld self.muTest = muTestOld self.mus = self.mus[:-1] self.samplingEngine.popSample() self.trainedModel.data = copy(trainedModelOld.data) break trainedModelOld.data = copy(self.trainedModel.data) if (self.interactive and maxErrorEst <= self.greedyTol): verbosityDepth("MAIN", ("Required tolerance {} achieved. Want " "to decrease greedyTol and continue? " "Y/N").format(self.greedyTol), timestamp = self.timestamp, end = "") increasemaxIter = input() if increasemaxIter.upper() == "Y": verbosityDepth("MAIN", "Reducing value of greedyTol...", timestamp = self.timestamp) while maxErrorEst <= self._greedyTol: self._greedyTol *= .5 if (self.interactive and self.samplingEngine.nsamples >= self.maxIter): verbosityDepth("MAIN", ("Maximum number of iterations {} " "reached. Want to increase maxIter " "and continue? Y/N").format( self.maxIter), timestamp = self.timestamp, end = "") increasemaxIter = input() if increasemaxIter.upper() == "Y": verbosityDepth("MAIN", "Doubling value of maxIter...", timestamp = self.timestamp) self._maxIter *= 2 if self.verbosity >= 2: verbosityDepth("DEL", ("Done computing snapshots (final snapshot " "count: {}).").format( self.samplingEngine.nsamples), timestamp = self.timestamp) def checkComputedApprox(self) -> bool: """ Check if setup of new approximant is not needed. Returns: True if new setup is not needed. False otherwise. """ return (super().checkComputedApprox() and len(self.mus) == self.trainedModel.data.projMat.shape[1]) - def computeScaleFactor(self): - """Compute parameter rescaling factor.""" - modeAssert(self._mode, message = "Cannot compute rescaling factor.") - self.scaleFactor= .5 * np.abs( - np.power(self.muBounds[0], self.HFEngine.rescalingExp) - - np.power(self.muBounds[1], self.HFEngine.rescalingExp)) - def assembleReducedResidualGramian(self, pMat:Np2D): """ Build residual gramian of reduced linear system through projections. """ self.initEstimatorNormEngine() if (not hasattr(self.trainedModel.data, "gramian") or self.trainedModel.data.gramian is None): gramian = self.estimatorNormEngine.innerProduct(pMat, pMat) else: Sold = self.trainedModel.data.gramian.shape[0] S = len(self.mus) if Sold > S: gramian = self.trainedModel.data.gramian[: S, : S] else: gramian = np.empty((S, S), dtype = np.complex) gramian[: Sold, : Sold] = self.trainedModel.data.gramian gramian[: Sold, Sold :] = ( self.estimatorNormEngine.innerProduct(pMat[:, Sold :], pMat[:, : Sold])) gramian[Sold :, : Sold] = gramian[: Sold, Sold :].T.conj() gramian[Sold :, Sold :] = ( self.estimatorNormEngine.innerProduct(pMat[:, Sold :], pMat[:, Sold :])) self.trainedModel.data.gramian = gramian def assembleReducedResidualBlocksbb(self, bs:List[Np1D], pMat:Np2D, scaling : float = 1.): """ Build blocks (of type bb) of reduced linear system through projections. """ self.initEstimatorNormEngine() nbs = len(bs) if (not hasattr(self.trainedModel.data, "resbb") or self.trainedModel.data.resbb is None): resbb = np.empty((nbs, nbs), dtype = np.complex) for i in range(nbs): Mbi = scaling ** i * bs[i] resbb[i, i] = self.estimatorNormEngine.innerProduct(Mbi, Mbi) for j in range(i): Mbj = scaling ** j * bs[j] resbb[i, j] = self.estimatorNormEngine.innerProduct(Mbj, Mbi) for i in range(nbs): for j in range(i + 1, nbs): resbb[i, j] = resbb[j, i].conj() self.trainedModel.data.resbb = resbb def assembleReducedResidualBlocksAb(self, As:List[Np2D], bs:List[Np1D], pMat:Np2D, scaling : float = 1.): """ Build blocks (of type Ab) of reduced linear system through projections. """ self.initEstimatorNormEngine() nAs = len(As) nbs = len(bs) S = len(self.mus) if (not hasattr(self.trainedModel.data, "resAb") or self.trainedModel.data.resAb is None): resAb = np.empty((nbs, S, nAs), dtype = np.complex) for j in range(nAs): MAj = scaling ** (j + 1) * As[j].dot(pMat) for i in range(nbs): Mbi = scaling ** (i + 1) * bs[i] resAb[i, :, j] = self.estimatorNormEngine.innerProduct(MAj, Mbi) else: Sold = self.trainedModel.data.resAb.shape[1] if Sold == S: return if Sold > S: resAb = self.trainedModel.data.resAb[:, : S, :] else: resAb = np.empty((nbs, S, nAs), dtype = np.complex) resAb[:, : Sold, :] = self.trainedModel.data.resAb for j in range(nAs): MAj = scaling ** (j + 1) * As[j].dot(pMat[:, Sold :]) for i in range(nbs): Mbi = scaling ** (i + 1) * bs[i] resAb[i, Sold :, j] = ( self.estimatorNormEngine.innerProduct(MAj, Mbi)) self.trainedModel.data.resAb = resAb def assembleReducedResidualBlocksAA(self, As:List[Np2D], pMat:Np2D, scaling : float = 1., basic : bool = False): """ Build blocks (of type AA) of reduced linear system through projections. """ self.initEstimatorNormEngine() nAs = len(As) S = len(self.mus) if (not hasattr(self.trainedModel.data, "resAA") or self.trainedModel.data.resAA is None): if basic: MAEnd = scaling ** nAs * As[-1].dot(pMat) resAA = self.estimatorNormEngine.innerProduct(MAEnd, MAEnd) else: resAA = np.empty((S, nAs, S, nAs), dtype = np.complex) for i in range(nAs): MAi = scaling ** (i + 1) * As[i].dot(pMat) resAA[:, i, :, i] = ( self.estimatorNormEngine.innerProduct(MAi, MAi)) for j in range(i): MAj = scaling ** (j + 1) * As[j].dot(pMat) resAA[:, i, :, j] = ( self.estimatorNormEngine.innerProduct(MAj, MAi)) for i in range(nAs): for j in range(i + 1, nAs): resAA[:, i, :, j] = resAA[:, j, :, i].T.conj() else: Sold = self.trainedModel.data.resAA.shape[0] if Sold == S: return if Sold > S: if basic: resAA = self.trainedModel.data.resAA[: S, : S] else: resAA = self.trainedModel.data.resAA[: S, :, : S, :] else: if basic: resAA = np.empty((S, S), dtype = np.complex) resAA[: Sold, : Sold] = self.trainedModel.data.resAA MAi = scaling ** nAs * As[-1].dot(pMat) resAA[: Sold, Sold :] = ( self.estimatorNormEngine.innerProduct(MAi[:, Sold :], MAi[:, : Sold])) resAA[Sold :, : Sold] = resAA[: Sold, Sold :].T.conj() resAA[Sold :, Sold :] = ( self.estimatorNormEngine.innerProduct(MAi[:, Sold :], MAi[:, Sold :])) else: resAA = np.empty((S, nAs, S, nAs), dtype = np.complex) resAA[: Sold, :, : Sold, :] = self.trainedModel.data.resAA for i in range(nAs): MAi = scaling ** (i + 1) * As[i].dot(pMat) resAA[: Sold, i, Sold :, i] = ( self.estimatorNormEngine.innerProduct(MAi[:, Sold :], MAi[:, : Sold])) resAA[Sold :, i, : Sold, i] = resAA[: Sold, i, Sold :, i].T.conj() resAA[Sold :, i, Sold :, i] = ( self.estimatorNormEngine.innerProduct(MAi[:, Sold :], MAi[:, Sold :])) for j in range(i): MAj = scaling ** (j + 1) * As[j].dot(pMat) resAA[: Sold, i, Sold :, j] = ( self.estimatorNormEngine.innerProduct(MAj[:, Sold :], MAi[:, : Sold])) resAA[Sold :, i, : Sold, j] = ( self.estimatorNormEngine.innerProduct(MAj[:, : Sold], MAi[:, Sold :])) resAA[Sold :, i, Sold :, j] = ( self.estimatorNormEngine.innerProduct(MAj[:, Sold :], MAi[:, Sold :])) for i in range(nAs): for j in range(i + 1, nAs): resAA[: Sold, i, Sold :, j] = ( resAA[Sold :, j, : Sold, i].T.conj()) resAA[Sold :, i, : Sold, j] = ( resAA[: Sold, j, Sold :, i].T.conj()) resAA[Sold :, i, Sold :, j] = ( resAA[Sold :, j, Sold :, i].T.conj()) self.trainedModel.data.resAA = resAA diff --git a/tests/test_2_hfengines/helmholtz_external.py b/tests/test_2_hfengines/helmholtz_external.py index aa1e4b5..e297ebc 100644 --- a/tests/test_2_hfengines/helmholtz_external.py +++ b/tests/test_2_hfengines/helmholtz_external.py @@ -1,45 +1,45 @@ # Copyright (C) 2018 by the RROMPy authors # # This file is part of RROMPy. # # RROMPy is free software: you can redistribute it and/or modify # it under the terms of the GNU Lesser General Public License as published by # the Free Software Foundation, either version 3 of the License, or # (at your option) any later version. # # RROMPy is distributed in the hope that it will be useful, # but WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # GNU Lesser General Public License for more details. # # You should have received a copy of the GNU Lesser General Public License # along with RROMPy. If not, see . # import numpy as np from rrompy.hfengines.linear_problem import ( HelmholtzCavityScatteringProblemEngine, HelmholtzBoxScatteringProblemEngine) def test_helmholtz_square_scattering(): solver = HelmholtzCavityScatteringProblemEngine(kappa = 4, gamma = 2., - n = 50, verbosity = 0) + n = 20, verbosity = 0) mu = 5 uh = solver.solve(mu) - assert np.isclose(solver.norm(uh), 20.6980450234954, rtol = 1e-5) + assert np.isclose(solver.norm(uh), 20.719752682674923, rtol = 1e-5) assert np.isclose(solver.norm(solver.residual(uh, mu)), 4.25056407e-13, rtol = 1e-1) def test_helmholtz_box_scattering(): solver = HelmholtzBoxScatteringProblemEngine(R = 2, kappa = 10., - theta = np.pi * 30 / 180, n = 50, verbosity = 0) + theta = np.pi * 30 / 180, n = 20, verbosity = 0) mu = 15 uh = solver.solve(mu) solver.plotmesh(show = False, figsize = (7, 7)) - assert np.isclose(solver.norm(uh), 64.05173319241996, rtol = 1e-5) + assert np.isclose(solver.norm(uh), 63.98946657389119, rtol = 1e-5) assert np.isclose(solver.norm(solver.residual(uh, mu)), 9.62989935e-13, rtol = 1e-1) from matplotlib import pyplot as plt plt.close('all')