diff --git a/README.md b/README.md index 2ae0ba7..4998fbb 100644 --- a/README.md +++ b/README.md @@ -1,31 +1,36 @@ RROMPy ============ Module for the solution and rational model order reduction of parametric PDE-based problem. Coded in Python 3.6. ## Prerequisites -**RBniCS** requires -* **numpy**. -* **scipy**. +**RROMPy** requires +* **numpy**; +* **scipy**; +* **fenics**. + +### Fenics +Most of the PDE engines already provided rely on [FEniCS](http://fenicsproject.org/). If you do not have FEniCS installed, you may want to create an [Anaconda3/Miniconda3](http://anaconda.org/) environment using the provided [file](conda-fenics.yml) by running the command +``` +conda env create --file conda-fenics.yml +``` +This will create an environment where Fenics can be used. In order to use FEniCS, the environment must be activated through +``` +source activate fenicsenv +``` ## Installing Clone the repository ``` git clone https://c4science.ch/source/RROMPy.git ``` and install the package by typing ``` python3 setup.py install ``` -### Fenics -Most of the PDE engines already provided rely on [FEniCS](https://fenicsproject.org/). If you do not have FEniCS installed, you may want to create an [Anaconda3/Miniconda3](http://anaconda.org/) environment using the provided [file](./conda-fenics.yml) by running the command -``` -conda env create --file conda-fenics.yml -``` -This will create an environment where Fenics can be used. In order to use FEniCS, the environment must be activated through -``` -source activate fenicsenv -``` ## License -This project is licensed under the GNU GENERAL PUBLIC LICENSE license - see the [LICENSE](./LICENSE) file for details. +This project is licensed under the GNU GENERAL PUBLIC LICENSE license - see the [LICENSE](LICENSE) file for details. + +## Acknowledgments +Part of the funding that made this module possible has been provided by the Swiss National Science Foundation through the FNS Research Project No. 200021. diff --git a/examples/scipy/solver.py b/examples/base/solver.py similarity index 100% rename from examples/scipy/solver.py rename to examples/base/solver.py diff --git a/examples/greedy/squareBubbleHomog.py b/examples/greedy/squareBubbleHomog.py new file mode 100644 index 0000000..53abf14 --- /dev/null +++ b/examples/greedy/squareBubbleHomog.py @@ -0,0 +1,95 @@ +import numpy as np +from rrompy.hfengines.scipy import HelmholtzSquareBubbleProblemEngine as HSBPE +from rrompy.reduction_methods.lagrange_greedy import ApproximantLagrangePadeGreedy as Pade +from rrompy.reduction_methods.lagrange_greedy import ApproximantLagrangeRBGreedy as RB +from rrompy.utilities.base import squareResonances + +verb = 2 +algo = "Pade" +#algo = "RB" + +k0s = np.power(np.linspace(95, 127, 25), .5) +k0 = np.mean(np.power(k0s, 2.)) ** .5 +kl, kr = min(k0s), max(k0s) + +polesexact = np.unique(np.power(squareResonances(kl**2., kr**2., False), .5)) + +rescaling = lambda x: np.power(x, 2.) +rescalingInv = lambda x: np.power(x, .5) +params = {'muBounds':[kl, kr], 'nTrainingPoints':500, + 'greedyTol':1e-2, 'nTestPoints':5} + +solver = HSBPE(kappa = 12 ** .5, theta = np.pi / 3, n = 30, + verbosity = verb) +solver.omega = np.real(k0) +if algo == "Pade": + approx = Pade(solver, mu0 = k0, approxParameters = params, + verbosity = verb) +else: + approx = RB(solver, mu0 = k0, approxParameters = params, verbosity = verb) + +#from time import clock +#start_time = clock() +approx.greedy(True) +#print("Time: ", clock() - start_time) + +approx.samplingEngine.verbosity = 0 +approx.verbosity = 0 +from matplotlib import pyplot as plt +MassInv = approx.massInv +normApp = np.zeros_like(k0s) +norm = np.zeros_like(k0s) +res = np.zeros_like(k0s) +err = np.zeros_like(k0s) +for j in range(len(k0s)): + normApp[j] = approx.normApp(k0s[j]) + norm[j] = approx.normHF(k0s[j]) + res[j] = (solver.norm(MassInv.solve(approx.getRes(k0s[j]))) + / solver.norm(MassInv.solve(approx.getRHS(k0s[j])))) + err[j] = approx.normErr(k0s[j]) / approx.normHF(k0s[j]) +resApp = approx.errorEstimator(k0s) + +plt.figure() +plt.semilogy(k0s, norm) +plt.semilogy(k0s, normApp, '--') +plt.semilogy(polesexact, + 2.*np.max(norm)*np.ones_like(polesexact, dtype = float), 'm.') +plt.semilogy(np.real(approx.mus), + 4.*np.max(norm)*np.ones_like(approx.mus, dtype = float), 'rx') +plt.xlim([kl, kr]) +plt.grid() +plt.show() +plt.close() + +plt.figure() +plt.semilogy(k0s, res) +plt.semilogy(k0s, resApp, '--') +plt.semilogy(polesexact, + 2.*np.max(resApp)*np.ones_like(polesexact, dtype = float), 'm.') +plt.semilogy(np.real(approx.mus), + 4.*np.max(resApp)*np.ones_like(approx.mus, dtype = float), 'rx') +plt.xlim([kl, kr]) +plt.grid() +plt.show() +plt.close() + +plt.figure() +plt.semilogy(k0s, err) +plt.semilogy(polesexact, + 2.*np.max(err)*np.ones_like(polesexact, dtype = float), 'm.') +plt.xlim([kl, kr]) +plt.grid() +plt.show() +plt.close() + +polesApp = approx.getPoles() +plt.figure() +plt.plot(np.real(polesApp), np.imag(polesApp), 'kx') +plt.plot(np.real(polesexact), np.imag(polesexact), 'm.') +plt.xlim(kl, kr) +plt.ylim(- (kr - kl) * .5, + (kr - kl) * .5) +#plt.axis('scaled') +plt.grid() +plt.show() +plt.close() + diff --git a/examples/greedy/squareTransmissionNonHomog.py b/examples/greedy/squareTransmissionNonHomog.py new file mode 100644 index 0000000..51a00f0 --- /dev/null +++ b/examples/greedy/squareTransmissionNonHomog.py @@ -0,0 +1,77 @@ +import numpy as np +from rrompy.hfengines.scipy import HelmholtzSquareTransmissionProblemEngine as HSTPE +from rrompy.reduction_methods.lagrange_greedy import ApproximantLagrangePadeGreedy as Pade + +verb = 7 +homog = True +#homog = False + +k0s = np.power(np.linspace(4, 15, 200), .5) +k0 = np.mean(np.power(k0s, 2.)) ** .5 +kl, kr = min(k0s), max(k0s) + +rescaling = lambda x: np.power(x, 2.) +rescalingInv = lambda x: np.power(x, .5) +params = {'muBounds':[kl, kr], 'nTrainingPoints':500, + 'greedyTol':1e-3, 'nTestPoints':10} + +solver = HSTPE(nT = 1, nB = 2, theta = np.pi * 20 / 180, kappa = 4., + n = 20, verbosity = verb) +solver.omega = np.real(k0) +approx = Pade(solver, mu0 = k0, approxParameters = params, + verbosity = verb, homogeneized = homog) +approx.greedy(True) + +approx.samplingEngine.verbosity = 0 +approx.verbosity = 0 +from matplotlib import pyplot as plt +MassInv = approx.massInv +normApp = np.zeros_like(k0s) +norm = np.zeros_like(k0s) +res = np.zeros_like(k0s) +err = np.zeros_like(k0s) +for j in range(len(k0s)): + normApp[j] = approx.normApp(k0s[j]) + norm[j] = approx.normHF(k0s[j]) + res[j] = (solver.norm(MassInv.solve(approx.getRes(k0s[j], + homogeneized = homog))) + / solver.norm(MassInv.solve(approx.getRHS(k0s[j], + homogeneized = homog)))) + err[j] = (approx.normErr(k0s[j], homogeneized = homog) + / approx.normHF(k0s[j], homogeneized = homog)) +resApp = approx.errorEstimator(k0s, homogeneized = homog) + +polesApp = approx.getPoles() +polesApp = polesApp[np.abs(np.imag(polesApp)) < 1e-3] +plt.figure() +plt.semilogy(k0s, norm) +plt.semilogy(k0s, normApp, '--') +plt.semilogy(np.real(approx.mus), + 4.*np.max(norm)*np.ones_like(approx.mus, dtype = float), 'rx') +plt.semilogy(np.real(polesApp), + 2.*np.max(norm)*np.ones_like(polesApp, dtype = float), 'k.') +plt.xlim([kl, kr]) +plt.grid() +plt.show() +plt.close() + +plt.figure() +plt.semilogy(k0s, res) +plt.semilogy(k0s, resApp, '--') +plt.semilogy(np.real(approx.mus), + 4.*np.max(resApp)*np.ones_like(approx.mus, dtype = float), 'rx') +plt.semilogy(np.real(polesApp), + 2.*np.max(resApp)*np.ones_like(polesApp, dtype = float), 'k.') +plt.xlim([kl, kr]) +plt.grid() +plt.show() +plt.close() + +plt.figure() +plt.semilogy(k0s, err) +plt.semilogy(np.real(polesApp), + 2.*np.max(err)*np.ones_like(polesApp, dtype = float), 'k.') +plt.xlim([kl, kr]) +plt.grid() +plt.show() +plt.close() diff --git a/examples/pod/LagrangePoles.py b/examples/pod/LagrangePoles.py new file mode 100644 index 0000000..141a5f6 --- /dev/null +++ b/examples/pod/LagrangePoles.py @@ -0,0 +1,48 @@ +from matplotlib import pyplot as plt +import numpy as np +from rrompy.hfengines.scipy import HelmholtzSquareBubbleProblemEngine as HSBPE +from rrompy.reduction_methods.lagrange import ApproximantLagrangePade as Pade +from rrompy.utilities.parameter_sampling import QuadratureSampler as QS +from rrompy.utilities.base import squareResonances + +verb = 0 +sampling = "Uniform" + +ks = [1, 46 ** .5] +solver = HSBPE(kappa = 12 ** .5, theta = np.pi / 3, n = 20, + verbosity = verb) +k0 = np.mean(np.power(ks, 2.)) ** .5 +k0 = 3.46104724 +solver.omega = np.real(k0) + +rescaling = lambda x: np.power(x, 2.) +rescalingInv = lambda x: np.power(x, .5) +if sampling == "Uniform": + sampler = QS(ks, "UNIFORM", rescaling, rescalingInv) + +nsets = 15 +paramsPade = {'S':2, 'POD':True, 'sampler':sampler} +approx = Pade(solver, mu0 = k0, approxParameters = paramsPade, + verbosity = verb) + +poles = [None] * (nsets - 1) +polesexact = np.unique(np.power(squareResonances(ks[0]**2., ks[1]**2., False), + .5)) +for i in range(1, nsets): + print("N = {}".format(4 * i)) + approx.approxParameters = {'N': 4 * i, 'M': 4 * i, 'S': 4 * i + 1} + approx.setupApprox() + poles[i - 1] = approx.getPoles() + + +for i in range(1, nsets): + plt.figure() + plt.plot(np.real(poles[i - 1]), np.imag(poles[i - 1]), 'kx') + plt.plot(polesexact, np.zeros_like(polesexact), 'm.') + plt.plot(k0, 0, 'r*') + plt.xlim(ks) + plt.ylim((ks[0] - ks[1]) / 2., (ks[1] - ks[0]) / 2.) + plt.title("N = {}, Neff = {}".format(4 * i, len(poles[i - 1]))) + plt.grid() + plt.show() + plt.close() \ No newline at end of file diff --git a/examples/scipy/LagrangeSweep.py b/examples/pod/LagrangeSweep.py similarity index 91% rename from examples/scipy/LagrangeSweep.py rename to examples/pod/LagrangeSweep.py index d1d09c9..45be186 100644 --- a/examples/scipy/LagrangeSweep.py +++ b/examples/pod/LagrangeSweep.py @@ -1,102 +1,102 @@ from copy import copy import numpy as np #from rrompy.hfengines.scipy import HelmholtzBoxScatteringProblemEngine as HBSPE from rrompy.hfengines.scipy import HelmholtzCavityScatteringProblemEngine as HBSPE from rrompy.reduction_methods.lagrange import ApproximantLagrangePade as Pade from rrompy.reduction_methods.lagrange import ApproximantLagrangeRB as RB from rrompy.utilities.parameter_sweeper import ParameterSweeper as Sweeper from rrompy.utilities.parameter_sampling import QuadratureSampler as QS from rrompy.utilities.parameter_sampling import WarpingFunction as WF from rrompy.utilities.parameter_sampling import WarpingSampler as WS verb = 0 npoints = 100 homog = True homog = False LSratio = 2. / 3. sampling = "Uniform" #sampling = "Cheby" #sampling = "SinC" assert LSratio <= 1. + np.finfo(float).eps ks = [10 + 0.j, 14 + 0.j] solver = HBSPE(kappa = 3, n = 15) solver.omega = np.real(np.mean(ks)) mutars = np.linspace(9, 15, npoints) homogMSG = "Homog" if not homog: homogMSG = "Non" + homogMSG -#filenamebase = '../data/output/HelmholtzBoxLSweep' -filenamebase = '../data/plots/LagrangeScatteringSquare1p5' -filenamebase = filenamebase + sampling + "/HelmholtzBoxLSweep" + homogMSG +filenamebase = '../data/output/ScatteringSquareLSweep' +#filenamebase = '../data/plots/LagrangeScatteringSquare1p5' +#filenamebase = filenamebase + sampling + "/HelmholtzBoxLSweep" + homogMSG k0 = np.mean(ks) shift = 7 shift = np.int(8 / LSratio - 1) nsets = 3 stride = np.int(8 / LSratio) Smax = stride * (nsets - 1) + shift + 2 rescaling = lambda x: np.power(x, 2.) rescalingInv = lambda x: np.power(x, .5) if sampling == "Uniform": sampler = QS(ks, "UNIFORM", rescaling, rescalingInv) if sampling == "Cheby": sampler = QS(ks, "CHEBYSHEV", rescaling, rescalingInv) if sampling == "SinC": warping = WF(call = lambda x: x - 2. * (1. - LSratio) / np.pi * np.sin(np.pi * x), repr = "x-{}*sin(pi*x)".format(2. * (1. - LSratio) / np.pi)) sampler = WS(ks, warping, rescaling, rescalingInv) paramsPade = {'S':Smax, 'POD':True, 'sampler':sampler} paramsRB = copy(paramsPade) paramsPoly = copy(paramsPade) paramsSetsPade = [None] * nsets paramsSetsRB = [None] * nsets paramsSetsPoly = [None] * nsets for i in range(nsets): paramsSetsPade[i] = {'N': np.int(LSratio * (stride * i + shift + 1)), 'M': np.int(LSratio * (stride * i + shift + 1)), 'S': stride * i + shift + 2} paramsSetsRB[i] = {'R': np.int(LSratio * (stride * i + shift + 1)), 'S': stride * i + shift + 2} paramsSetsPoly[i] = {'N': 0, 'M': np.int(LSratio * (stride * i + shift + 1)), 'S': stride * i + shift + 2} appPade = Pade(solver, mu0 = k0, approxParameters = paramsPade, - verbosity = verb, homogeneize = homog) + verbosity = verb, homogeneized = homog) appRB = RB(solver, mu0 = k0, approxParameters = paramsRB, - verbosity = verb, homogeneize = homog) + verbosity = verb, homogeneized = homog) appPoly = Pade(solver, mu0 = k0, approxParameters = paramsPoly, - verbosity = verb, homogeneize = homog) + verbosity = verb, homogeneized = homog) sweeper = Sweeper(mutars = mutars, mostExpensive = 'Approx') sweeper.ROMEngine = appPade sweeper.params = paramsSetsPade filenamePade = sweeper.sweep(filenamebase + 'Pade.dat') sweeper.ROMEngine = appRB sweeper.params = paramsSetsRB filenameRB = sweeper.sweep(filenamebase + 'RB.dat') sweeper.ROMEngine = appPoly sweeper.params = paramsSetsPoly filenamePoly = sweeper.sweep(filenamebase + 'Poly.dat') sweeper.plotCompare([filenamePade, filenameRB, filenamePoly], ['muRe'], ['normHF', 'normApp'], ['S'], onePlot = True, save = filenamebase + 'Norm', ylims = {'top' : 1e1}, saveFormat = "png", labels = ["Pade'", "RB", "Poly"], # figsize = (5, 3.75)) figsize = (10, 7.5)) sweeper.plotCompare([filenamePade, filenameRB, filenamePoly], ['muRe'], ['normResRel'], ['S'], save = filenamebase + 'Res', ylims = {'top' : 1e1}, saveFormat = "png", labels = ["Pade'", "RB", "Poly"], # figsize = (5, 3.75)) figsize = (10, 7.5)) sweeper.plotCompare([filenamePade, filenameRB, filenamePoly], ['muRe'], ['normErrRel'], ['S'], save = filenamebase + 'Err', ylims = {'top' : 1e1}, saveFormat = "png", labels = ["Pade'", "RB", "Poly"], # figsize = (5, 3.75)) figsize = (10, 7.5)) diff --git a/examples/scipy/LagrangeSweepUnstable.py b/examples/pod/LagrangeSweepUnstable.py similarity index 100% rename from examples/scipy/LagrangeSweepUnstable.py rename to examples/pod/LagrangeSweepUnstable.py diff --git a/examples/scipy/PadeLagrange.py b/examples/pod/PadeLagrange.py similarity index 96% rename from examples/scipy/PadeLagrange.py rename to examples/pod/PadeLagrange.py index e3f330f..300165c 100644 --- a/examples/scipy/PadeLagrange.py +++ b/examples/pod/PadeLagrange.py @@ -1,107 +1,107 @@ import numpy as np from rrompy.hfengines.scipy import HelmholtzSquareBubbleProblemEngine as HSBPE from rrompy.hfengines.scipy import HelmholtzSquareTransmissionProblemEngine as HSTPE from rrompy.hfengines.scipy import HelmholtzBoxScatteringProblemEngine as HBSPE from rrompy.reduction_methods.lagrange import ApproximantLagrangePade as Pade from rrompy.utilities.parameter_sampling import QuadratureSampler as QS -testNo = 2 +testNo = 1 verb = 0 homog = True homog = False if testNo == 1: k0s = np.power([10 + 0.j, 14 + 0.j], .5) k0 = np.mean(k0s) ktar = (11 + .5j) ** .5 rescaling = lambda x: np.power(x, 2.) rescalingInv = lambda x: np.power(x, .5) params = {'N':4, 'M':3, 'S':5, 'POD':True, 'sampler':QS(k0s, "CHEBYSHEV", rescaling, rescalingInv)} solver = HSBPE(kappa = 12 ** .5, theta = np.pi / 3, n = 40, verbosity = verb) solver.omega = np.real(k0) approx = Pade(solver, mu0 = k0, approxParameters = params, verbosity = verb) approx.setupApprox() # approx.plotSamples() approx.plotApp(ktar, name = 'u_Pade''') approx.plotHF(ktar, name = 'u_HF') approx.plotErr(ktar, name = 'err') approx.plotRes(ktar, name = 'res') appErr, solNorm = approx.normErr(ktar), approx.normHF(ktar) resNorm, RHSNorm = approx.normRes(ktar), approx.normRHS(ktar) print(('SolNorm:\t{}\nErr:\t{}\nErrRel:\t{}').format(solNorm, appErr, np.divide(appErr, solNorm))) print(('RHSNorm:\t{}\nRes:\t{}\nResRel:\t{}').format(RHSNorm, resNorm, np.divide(resNorm, RHSNorm))) print('\nPoles Pade'':') print(approx.getPoles()) ############ elif testNo == 2: k0s = [3.85 + 0.j, 4.15 + 0.j] k0 = np.mean(k0s) ktar = 4 + .15j rescaling = lambda x: np.power(x, 2.) rescalingInv = lambda x: np.power(x, .5) params = {'N':8, 'M':9, 'S':10, 'POD':True, 'sampler':QS(k0s, "CHEBYSHEV", rescaling, rescalingInv)} solver = HSTPE(nT = 2, nB = 1, theta = np.pi * 45/180, kappa = 4., n = 50, verbosity = verb) solver.omega = np.real(k0) approx = Pade(solver, mu0 = k0, approxParameters = params, - verbosity = verb, homogeneize = homog) + verbosity = verb, homogeneized = homog) approx.setupApprox() # approx.plotSamples() approx.plotApp(ktar, name = 'u_Pade''') approx.plotHF(ktar, name = 'u_HF') approx.plotErr(ktar, name = 'err') approx.plotRes(ktar, name = 'res') appErr, solNorm = approx.normErr(ktar), approx.normHF(ktar) resNorm, RHSNorm = approx.normRes(ktar), approx.normRHS(ktar) print(('SolNorm:\t{}\nErr:\t{}\nErrRel:\t{}').format(solNorm, appErr, np.divide(appErr, solNorm))) print(('RHSNorm:\t{}\nRes:\t{}\nResRel:\t{}').format(RHSNorm, resNorm, np.divide(resNorm, RHSNorm))) print('\nPoles Pade'':') print(approx.getPoles()) ############ elif testNo == 3: k0s = [2, 5] k0 = np.mean(k0s) ktar = 4.5 - 0.j params = {'N':10, 'M':11, 'S':15, 'POD':True, 'sampler':QS(k0s, "CHEBYSHEV")} solver = HBSPE(R = 7, kappa = 3, theta = - np.pi * 75 / 180, n = 40, verbosity = verb) solver.omega = np.real(k0) approx = Pade(solver, mu0 = k0, approxParameters = params, - verbosity = verb, homogeneize = homog) + verbosity = verb, homogeneized = homog) approx.setupApprox() approx.plotSamples() approx.plotApp(ktar, name = 'u_Pade''') approx.plotHF(ktar, name = 'u_HF') approx.plotErr(ktar, name = 'err') approx.plotRes(ktar, name = 'res') appErr, solNorm = approx.normErr(ktar), approx.normHF(ktar) resNorm, RHSNorm = approx.normRes(ktar), approx.normRHS(ktar) print(('SolNorm:\t{}\nErr:\t{}\nErrRel:\t{}').format(solNorm, appErr, np.divide(appErr, solNorm))) print(('RHSNorm:\t{}\nRes:\t{}\nResRel:\t{}').format(RHSNorm, resNorm, np.divide(resNorm, RHSNorm))) print('\nPoles Pade'':') print(approx.getPoles()) diff --git a/examples/scipy/PadeTaylor.py b/examples/pod/PadeTaylor.py similarity index 96% rename from examples/scipy/PadeTaylor.py rename to examples/pod/PadeTaylor.py index 4fe5a30..5461db0 100644 --- a/examples/scipy/PadeTaylor.py +++ b/examples/pod/PadeTaylor.py @@ -1,97 +1,97 @@ import numpy as np from rrompy.hfengines.scipy import HelmholtzSquareBubbleProblemEngine as HSBPE from rrompy.hfengines.scipy import ( HelmholtzSquareTransmissionProblemEngine as HSTPE) from rrompy.hfengines.scipy import HelmholtzBoxScatteringProblemEngine as HBSPE from rrompy.reduction_methods.taylor import ApproximantTaylorPade as Pade testNo = 2 verb = 0 homog = True #homog = False if testNo == 1: params = {'N':4, 'M':3, 'E':4, 'sampleType':'Arnoldi', 'POD':True} k0 = 12 ** .5 ktar = 10.5 ** .5 solver = HSBPE(kappa = 12 ** .5, theta = np.pi / 3, n = 40, verbosity = verb) solver.omega = np.real(k0) approx = Pade(solver, mu0 = k0, approxParameters = params, verbosity = verb) approx.setupApprox() # approx.plotSamples() approx.plotApp(ktar, name = 'u_Pade''') approx.plotHF(ktar, name = 'u_HF') approx.plotErr(ktar, name = 'err') approx.plotRes(ktar, name = 'res') appErr, solNorm = approx.normErr(ktar), approx.normHF(ktar) resNorm, RHSNorm = approx.normRes(ktar), approx.normRHS(ktar) print(('SolNorm:\t{}\nErr:\t{}\nErrRel:\t{}').format(solNorm, appErr, np.divide(appErr, solNorm))) print(('RHSNorm:\t{}\nRes:\t{}\nResRel:\t{}').format(RHSNorm, resNorm, np.divide(resNorm, RHSNorm))) print('\nPoles Pade'':') print(approx.getPoles()) ############ elif testNo == 2: params = {'N':6, 'M':7, 'E':7, 'sampleType':'Arnoldi', 'POD':True} k0 = 16 ** .5 ktar = 15 ** .5 solver = HSTPE(nT = 2, nB = 1, theta = np.pi * 45/180, kappa = 4., n = 50, verbosity = verb) solver.omega = np.real(k0) approx = Pade(solver, mu0 = k0, approxParameters = params, - verbosity = verb, homogeneize = homog) + verbosity = verb, homogeneized = homog) approx.setupApprox() # approx.plotSamples() approx.plotApp(ktar, name = 'u_Pade''') approx.plotHF(ktar, name = 'u_HF') approx.plotErr(ktar, name = 'err') approx.plotRes(ktar, name = 'res') appErr, solNorm = approx.normErr(ktar), approx.normHF(ktar) resNorm, RHSNorm = approx.normRes(ktar), approx.normRHS(ktar) print(('SolNorm:\t{}\nErr:\t{}\nErrRel:\t{}').format(solNorm, appErr, np.divide(appErr, solNorm))) print(('RHSNorm:\t{}\nRes:\t{}\nResRel:\t{}').format(RHSNorm, resNorm, np.divide(resNorm, RHSNorm))) print('\nPoles Pade'':') print(approx.getPoles()) ############ elif testNo in [3, 4]: if testNo == 3: params = {'N':7, 'M':8, 'E':8, 'sampleType':'Krylov', 'POD':True} else: params = {'N':7, 'M':8, 'E':8, 'sampleType':'Arnoldi', 'POD':True} k0 = 3 ktar = 4.+0.j solver = HBSPE(R = 5, kappa = 3, theta = - np.pi * 75 / 180, n = 30, verbosity = verb) solver.omega = np.real(k0) approx = Pade(solver, mu0 = k0, approxParameters = params, - verbosity = verb, homogeneize = homog) + verbosity = verb, homogeneized = homog) approx.setupApprox() approx.plotSamples() approx.plotApp(ktar, name = 'u_Pade''') approx.plotHF(ktar, name = 'u_HF') approx.plotErr(ktar, name = 'err') approx.plotRes(ktar, name = 'res') appErr, solNorm = approx.normErr(ktar), approx.normHF(ktar) resNorm, RHSNorm = approx.normRes(ktar), approx.normRHS(ktar) print(('SolNorm:\t{}\nErr:\t{}\nErrRel:\t{}').format(solNorm, appErr, np.divide(appErr, solNorm))) print(('RHSNorm:\t{}\nRes:\t{}\nResRel:\t{}').format(RHSNorm, resNorm, np.divide(resNorm, RHSNorm))) print('\nPoles Pade'':') print(approx.getPoles()) diff --git a/examples/scipy/RBLagrange.py b/examples/pod/RBLagrange.py similarity index 97% rename from examples/scipy/RBLagrange.py rename to examples/pod/RBLagrange.py index bbfe2bc..ac5474f 100644 --- a/examples/scipy/RBLagrange.py +++ b/examples/pod/RBLagrange.py @@ -1,99 +1,99 @@ import numpy as np from rrompy.hfengines.scipy import HelmholtzSquareBubbleProblemEngine as HSBPE from rrompy.hfengines.scipy import HelmholtzSquareTransmissionProblemEngine as HSTPE from rrompy.hfengines.scipy import HelmholtzBoxScatteringProblemEngine as HBSPE from rrompy.reduction_methods.lagrange import ApproximantLagrangeRB as RB from rrompy.utilities.parameter_sampling import QuadratureSampler as QS testNo = 3 verb = 0 homog = True homog = False if testNo == 1: k0s = np.power([10 + 0.j, 14 + 0.j], .5) k0 = np.mean(k0s) ktar = (11 + .5j) ** .5 rescaling = lambda x: np.power(x, 2.) rescalingInv = lambda x: np.power(x, .5) params = {'S':5, 'R':4, 'POD':True, 'sampler':QS(k0s, "CHEBYSHEV", rescaling, rescalingInv)} solver = HSBPE(kappa = 12 ** .5, theta = np.pi / 3, n = 40, verbosity = verb) solver.omega = np.real(k0) approx = RB(solver, mu0 = k0, approxParameters = params, verbosity = verb) approx.setupApprox() # approx.plotSamples() approx.plotApp(ktar, name = 'u_RB') approx.plotHF(ktar, name = 'u_HF') approx.plotErr(ktar, name = 'err') approx.plotRes(ktar, name = 'res') appErr, solNorm = approx.normErr(ktar), approx.normHF(ktar) resNorm, RHSNorm = approx.normRes(ktar), approx.normRHS(ktar) print(('SolNorm:\t{}\nErr:\t{}\nErrRel:\t{}').format(solNorm, appErr, np.divide(appErr, solNorm))) print(('RHSNorm:\t{}\nRes:\t{}\nResRel:\t{}').format(RHSNorm, resNorm, np.divide(resNorm, RHSNorm))) ############ elif testNo == 2: k0s = [3.85 + 0.j, 4.15 + 0.j] k0 = np.mean(k0s) ktar = 4 + .15j rescaling = lambda x: np.power(x, 2.) rescalingInv = lambda x: np.power(x, .5) params = {'S':10, 'R':9, 'POD':True, 'sampler':QS(k0s, "CHEBYSHEV", rescaling, rescalingInv)} solver = HSTPE(nT = 2, nB = 1, theta = np.pi * 45/180, kappa = 4., n = 50, verbosity = verb) solver.omega = np.real(k0) approx = RB(solver, mu0 = k0, approxParameters = params, - verbosity = verb, homogeneize = homog) + verbosity = verb, homogeneized = homog) approx.setupApprox() # approx.plotSamples() approx.plotApp(ktar, name = 'u_RB') approx.plotHF(ktar, name = 'u_HF') approx.plotErr(ktar, name = 'err') approx.plotRes(ktar, name = 'res') appErr, solNorm = approx.normErr(ktar), approx.normHF(ktar) resNorm, RHSNorm = approx.normRes(ktar), approx.normRHS(ktar) print(('SolNorm:\t{}\nErr:\t{}\nErrRel:\t{}').format(solNorm, appErr, np.divide(appErr, solNorm))) print(('RHSNorm:\t{}\nRes:\t{}\nResRel:\t{}').format(RHSNorm, resNorm, np.divide(resNorm, RHSNorm))) ############ elif testNo == 3: k0s = [2, 5] k0 = np.mean(k0s) ktar = 4.5 - 0.j params = {'S':15, 'R':10, 'POD':True, 'sampler':QS(k0s, "CHEBYSHEV")} solver = HBSPE(R = 7, kappa = 3, theta = - np.pi * 75 / 180, n = 40, verbosity = verb) solver.omega = np.real(k0) approx = RB(solver, mu0 = k0, approxParameters = params, - verbosity = verb, homogeneize = homog) + verbosity = verb, homogeneized = homog) approx.setupApprox() # approx.plotSamples() approx.plotApp(ktar, name = 'u_RB') approx.plotHF(ktar, name = 'u_HF') approx.plotErr(ktar, name = 'err') approx.plotRes(ktar, name = 'res') appErr, solNorm = approx.normErr(ktar), approx.normHF(ktar) resNorm, RHSNorm = approx.normRes(ktar), approx.normRHS(ktar) print(('SolNorm:\t{}\nErr:\t{}\nErrRel:\t{}').format(solNorm, appErr, np.divide(appErr, solNorm))) print(('RHSNorm:\t{}\nRes:\t{}\nResRel:\t{}').format(RHSNorm, resNorm, np.divide(resNorm, RHSNorm))) diff --git a/examples/scipy/RBTaylor.py b/examples/pod/RBTaylor.py similarity index 96% rename from examples/scipy/RBTaylor.py rename to examples/pod/RBTaylor.py index b95241c..00b729a 100644 --- a/examples/scipy/RBTaylor.py +++ b/examples/pod/RBTaylor.py @@ -1,90 +1,90 @@ import numpy as np from rrompy.hfengines.scipy import HelmholtzSquareBubbleProblemEngine as HSBPE from rrompy.hfengines.scipy import HelmholtzSquareTransmissionProblemEngine as HSTPE from rrompy.hfengines.scipy import HelmholtzBoxScatteringProblemEngine as HBSPE from rrompy.reduction_methods.taylor import ApproximantTaylorRB as RB testNo = 4 verb = 0 homog = True #homog = False if testNo == 1: params = {'E':4, 'R':4, 'sampleType':'Arnoldi', 'POD':True} k0 = 12 ** .5 ktar = 10.5 ** .5 solver = HSBPE(kappa = 12 ** .5, theta = np.pi / 3, n = 40, verbosity = verb) solver.omega = np.real(k0) approx = RB(solver, mu0 = k0, approxParameters = params, verbosity = verb) approx.setupApprox() # approx.plotSamples() approx.plotApp(ktar, name = 'u_RB') approx.plotHF(ktar, name = 'u_HF') approx.plotErr(ktar, name = 'err') approx.plotRes(ktar, name = 'res') appErr, solNorm = approx.normErr(ktar), approx.normHF(ktar) resNorm, RHSNorm = approx.normRes(ktar), approx.normRHS(ktar) print(('SolNorm:\t{}\nErr:\t{}\nErrRel:\t{}').format(solNorm, appErr, np.divide(appErr, solNorm))) print(('RHSNorm:\t{}\nRes:\t{}\nResRel:\t{}').format(RHSNorm, resNorm, np.divide(resNorm, RHSNorm))) ############ elif testNo == 2: params = {'E':7, 'R':7, 'sampleType':'Arnoldi', 'POD':True} k0 = 16**.5 ktar = 15**.5 solver = HSTPE(nT = 2, nB = 1, theta = np.pi * 45/180, kappa = 3., n = 50, verbosity = verb) solver.omega = np.real(k0) approx = RB(solver, mu0 = k0, approxParameters = params, - verbosity = verb, homogeneize = homog) + verbosity = verb, homogeneized = homog) approx.setupApprox() # approx.plotSamples() approx.plotApp(ktar, name = 'u_RB') approx.plotHF(ktar, name = 'u_HF') approx.plotErr(ktar, name = 'err') approx.plotRes(ktar, name = 'res') appErr, solNorm = approx.normErr(ktar), approx.normHF(ktar) resNorm, RHSNorm = approx.normRes(ktar), approx.normRHS(ktar) print(('SolNorm:\t{}\nErr:\t{}\nErrRel:\t{}').format(solNorm, appErr, np.divide(appErr, solNorm))) print(('RHSNorm:\t{}\nRes:\t{}\nResRel:\t{}').format(RHSNorm, resNorm, np.divide(resNorm, RHSNorm))) ############ elif testNo in [3, 4]: if testNo == 3: params = {'E':8, 'sampleType':'Krylov', 'POD':True} else: params = {'E':8, 'sampleType':'Arnoldi', 'POD':True} k0 = 3 ktar = 4.25+.5j solver = HBSPE(R = 5, kappa = 3, theta = - np.pi * 75 / 180, n = 30, verbosity = verb) solver.omega = np.real(k0) approx = RB(solver, mu0 = k0, approxParameters = params, - verbosity = verb, homogeneize = homog) + verbosity = verb, homogeneized = homog) approx.setupApprox() # approx.plotSamples() approx.plotApp(ktar, name = 'u_RB') approx.plotHF(ktar, name = 'u_HF') approx.plotErr(ktar, name = 'err') approx.plotRes(ktar, name = 'res') appErr, solNorm = approx.normErr(ktar), approx.normHF(ktar) resNorm, RHSNorm = approx.normRes(ktar), approx.normRHS(ktar) print(('SolNorm:\t{}\nErr:\t{}\nErrRel:\t{}').format(solNorm, appErr, np.divide(appErr, solNorm))) print(('RHSNorm:\t{}\nRes:\t{}\nResRel:\t{}').format(RHSNorm, resNorm, np.divide(resNorm, RHSNorm))) diff --git a/examples/scipy/TaylorPoles.py b/examples/pod/TaylorPoles.py similarity index 100% rename from examples/scipy/TaylorPoles.py rename to examples/pod/TaylorPoles.py diff --git a/examples/scipy/TaylorSweep.py b/examples/pod/TaylorSweep.py similarity index 94% rename from examples/scipy/TaylorSweep.py rename to examples/pod/TaylorSweep.py index 9e9c8a3..3cf37e6 100644 --- a/examples/scipy/TaylorSweep.py +++ b/examples/pod/TaylorSweep.py @@ -1,77 +1,77 @@ import numpy as np from rrompy.hfengines.scipy import HelmholtzSquareBubbleProblemEngine as HSBPE from rrompy.hfengines.scipy import HelmholtzBoxScatteringProblemEngine as HBSPE from rrompy.reduction_methods.taylor import ApproximantTaylorPade as Pade from rrompy.reduction_methods.taylor import ApproximantTaylorRB as RB from rrompy.utilities.parameter_sweeper import ParameterSweeper as Sweeper testNo = 1 verb = 0 homog = True homog = False k0 = (12 + 0.j) ** .5 npoints = 100 shift = 5 nsets = 2 stride = 2 Emax = stride * (nsets - 1) + shift + 1 if testNo == 1: solver = HSBPE(kappa = 12 ** .5, theta = np.pi / 3, n = 15, verbosity = verb) solver.omega = np.real(k0) params = {'Emax':Emax, 'sampleType':'ARNOLDI', 'POD':True} ktars = np.linspace(7**.5, 16**.5, npoints) filenamebase = '../data/output/HelmholtzBubbleTSweep' homog = False elif testNo == 2: solver = HBSPE(R = 5, kappa = 3, theta = - np.pi * 75 / 180, n = 10, verbosity = verb) solver.omega = np.real(k0) params = {'Emax':Emax, 'sampleType':'Arnoldi', 'POD':True} ktars = np.linspace(7**.5, 17**.5, npoints) homogMSG = "" if not homog: homogMSG = "Non" filenamebase = '../data/output/HelmholtzBoxTSweep' + homogMSG + "Homog" paramsSetsPade = [None] * nsets paramsSetsRB = [None] * nsets paramsSetsPoly = [None] * nsets for i in range(nsets): paramsSetsPade[i] = {'N':stride*i+shift+1, 'M':stride*i+shift+1, 'E':stride*i+shift+1} paramsSetsRB[i] = {'E':stride*i+shift+1,'R':stride*i+shift+2} paramsSetsPoly[i] = {'N':0, 'M':stride*i+shift+1, 'E':stride*i+shift+1} appPade = Pade(solver, mu0 = k0, approxParameters = params, - verbosity = verb, homogeneize = homog) + verbosity = verb, homogeneized = homog) appRB = RB(solver, mu0 = k0, approxParameters = params, - verbosity = verb, homogeneize = homog) + verbosity = verb, homogeneized = homog) appPoly = Pade(solver, mu0 = k0, approxParameters = params, - verbosity = verb, homogeneize = homog) + verbosity = verb, homogeneized = homog) sweeper = Sweeper(mutars = ktars, mostExpensive = 'Approx') sweeper.ROMEngine = appPade sweeper.params = paramsSetsPade filenamePade = sweeper.sweep(filenamebase + 'Pade.dat') sweeper.ROMEngine = appRB sweeper.params = paramsSetsRB filenameRB = sweeper.sweep(filenamebase + 'RB.dat') sweeper.ROMEngine = appPoly sweeper.params = paramsSetsPoly filenamePoly = sweeper.sweep(filenamebase + 'Poly.dat') sweeper.plotCompare([filenamePade, filenameRB, filenamePoly], ['muRe'], ['normHF', 'normApp'], ['E'], onePlot = True, save = filenamebase + 'Norm', saveFormat = "png", labels = ["Pade'", "RB", "Poly"]) sweeper.plotCompare([filenamePade, filenameRB, filenamePoly], ['muRe'], ['normResRel'], ['E'], save = filenamebase + 'Res', saveFormat = "png", labels = ["Pade'", "RB", "Poly"]) sweeper.plotCompare([filenamePade, filenameRB, filenamePoly], ['muRe'], ['normErrRel'], ['E'], save = filenamebase + 'Err', saveFormat = "png", labels = ["Pade'", "RB", "Poly"]) diff --git a/examples/scipy/airfoil.py b/examples/pod/airfoil.py similarity index 94% rename from examples/scipy/airfoil.py rename to examples/pod/airfoil.py index 498aa50..2121fff 100644 --- a/examples/scipy/airfoil.py +++ b/examples/pod/airfoil.py @@ -1,212 +1,212 @@ from copy import deepcopy as copy import numpy as np import fenics as fen import ufl from rrompy.hfengines.scipy import HelmholtzBoxScatteringProblemEngine as HSP from rrompy.reduction_methods.taylor import ApproximantTaylorPade as TP from rrompy.reduction_methods.taylor import ApproximantTaylorRB as TRB from rrompy.reduction_methods.lagrange import ApproximantLagrangePade as LP from rrompy.reduction_methods.lagrange import ApproximantLagrangeRB as LRB from rrompy.utilities.parameter_sweeper import ParameterSweeper as Sweeper from rrompy.utilities.parameter_sampling import QuadratureSampler as QS from rrompy.utilities.base.fenics import fenONE from operator import itemgetter def subdict(d, ks): return dict(zip(ks, itemgetter(*ks)(d))) verb = 0 #################### homog = True #homog = False #################### test = "solve" test = "Taylor" test = "Lagrange" test = "TaylorSweep" test = "LagrangeSweep" plotSamples = True k0 = 10 kLeft, kRight = 8 + 0.j, 12 + 0.j ktar = 11 ktars = np.linspace(8, 12, 21) + 0.j PI = np.pi R = 2 def Dboundary(x, on_boundary): return on_boundary and (x[0]**2+x[1]**2)**.5 < .95 * R kappa = 10 theta = PI * - 45 / 180. mu = 1.1 epsilon = .1 mesh = fen.Mesh('../data/mesh/airfoil.xml') c, s = np.cos(theta), np.sin(theta) x, y = fen.SpatialCoordinate(mesh)[:] u0R = - fen.cos(kappa * (c * x + s * y)) u0I = - fen.sin(kappa * (c * x + s * y)) checkReal = x**2-x+y**2 rhop5 = ((x**2+y**2)/((x-1)**2+y**2))**.25 phiroot1 = fen.atan(-y/(x**2-x+y**2)) / 2 phiroot2 = fen.atan(-y/(x**2-x+y**2)) / 2 - PI * ufl.sign(-y/(x**2-x+y**2)) / 2 kappam1 = (((rhop5*fen.cos(phiroot1)+.5)**2.+(rhop5*fen.sin(phiroot1))**2.)/ ((rhop5*fen.cos(phiroot1)-1.)**2.+(rhop5*fen.sin(phiroot1))**2.) )**.5 - mu kappam2 = (((rhop5*fen.cos(phiroot2)+.5)**2.+(rhop5*fen.sin(phiroot2))**2.)/ ((rhop5*fen.cos(phiroot2)-1.)**2.+(rhop5*fen.sin(phiroot2))**2.) )**.5 - mu Heps1 = .9 * .5 * (1 + kappam1/epsilon + fen.sin(PI*kappam1/epsilon) / PI) + .1 Heps2 = .9 * .5 * (1 + kappam2/epsilon + fen.sin(PI*kappam2/epsilon) / PI) + .1 cTT = ufl.conditional(ufl.le(kappam1, epsilon), Heps1, fenONE) c_F = fen.Constant(.1) cFT = ufl.conditional(ufl.le(kappam2, epsilon), Heps2, fenONE) c_F = fen.Constant(.1) cT = ufl.conditional(ufl.ge(kappam1, - epsilon), cTT, c_F) cF = ufl.conditional(ufl.ge(kappam2, - epsilon), cFT, c_F) a = ufl.conditional(ufl.ge(checkReal, 0.), cT, cF) ### solver = HSP(R, np.abs(k0), theta, n = 1, verbosity = verb, degree_threshold = 8) solver.V = fen.FunctionSpace(mesh, "P", 3) solver.diffusivity = a solver.DirichletBoundary = Dboundary solver.RobinBoundary = "REST" solver.DirichletDatum = [u0R, u0I] ### if test == "solve": uinc = solver.liftDirichletData(k0) if homog: uhtot = solver.solve(k0, homogeneized = homog) uh = uhtot + uinc else: uh = solver.solve(k0, homogeneized = homog) uhtot = uh - uinc print(solver.norm(uh)) print(solver.norm(uhtot)) solver.plot(fen.project(a, solver.V).vector(), what = 'Real', name = 'a') solver.plot(uinc, what = 'Real', name = 'u_inc') solver.plot(uh, what = 'ABS') solver.plot(uhtot, what = 'ABS', name = 'u + u_inc') elif test in ["Taylor", "Lagrange"]: if test == "Taylor": params = {'N':8, 'M':8, 'R':8, 'E':8, 'sampleType':'Arnoldi', 'POD':True} parPade = subdict(params, ['N', 'M', 'E', 'sampleType', 'POD']) parRB = subdict(params, ['R', 'E', 'sampleType', 'POD']) approxPade = TP(solver, mu0 = k0, approxParameters = parPade, - verbosity = verb, homogeneize = homog) + verbosity = verb, homogeneized = homog) approxRB = TRB(solver, mu0 = k0, approxParameters = parRB, - verbosity = verb, homogeneize = homog) + verbosity = verb, homogeneized = homog) else: params = {'N':8, 'M':8, 'R':9, 'S':9, 'POD':True, 'sampler':QS([kLeft, kRight], "CHEBYSHEV")} parPade = subdict(params, ['N', 'M', 'S', 'POD', 'sampler']) parRB = subdict(params, ['R', 'S', 'POD', 'sampler']) approxPade = LP(solver, mu0 = np.mean([kLeft, kRight]), approxParameters = parPade, verbosity = verb, - homogeneize = homog) + homogeneized = homog) approxRB = LRB(solver, mu0 = np.mean([kLeft, kRight]), approxParameters = parRB, verbosity = verb, - homogeneize = homog) + homogeneized = homog) approxPade.setupApprox() approxRB.setupApprox() if plotSamples: approxPade.plotSamples() approxPade.plotHF(ktar, name = 'u_HF') approxPade.plotApp(ktar, name = 'u_Pade''') approxPade.plotErr(ktar, name = 'err_Pade''') approxPade.plotRes(ktar, name = 'res_Pade''') approxRB.plotApp(ktar, name = 'u_RB') approxRB.plotErr(ktar, name = 'err_RB') approxRB.plotRes(ktar, name = 'res_RB') HFNorm, RHSNorm = approxPade.normHF(ktar), approxPade.normRHS(ktar) PadeRes, PadeErr = approxPade.normRes(ktar), approxPade.normErr(ktar) RBRes, RBErr = approxRB.normRes(ktar), approxRB.normErr(ktar) print('HFNorm:\t{}\nRHSNorm:\t{}'.format(HFNorm, RHSNorm)) print('PadeRes:\t{}\nPadeErr:\t{}'.format(PadeRes, PadeErr)) print('RBRes:\t{}\nRBErr:\t{}'.format(RBRes, RBErr)) print('\nPoles Pade'':') print(approxPade.getPoles()) elif test in ["TaylorSweep", "LagrangeSweep"]: if test == "TaylorSweep": shift = 10 nsets = 2 stride = 3 Emax = stride * (nsets - 1) + shift + 1 params = {'Emax':Emax, 'sampleType':'Arnoldi', 'POD':True} paramsSetsPade = [None] * nsets paramsSetsRB = [None] * nsets for i in range(nsets): paramsSetsPade[i] = {'N':stride*i+shift + 1, 'M':stride*i+shift+1, 'E':stride*i+shift + 1} paramsSetsRB[i] = {'E':stride*i+shift + 1,'R':stride*i+shift + 2} approxPade = TP(solver, mu0 = kappa,approxParameters = params, - verbosity = verb, homogeneize = homog) + verbosity = verb, homogeneized = homog) approxRB = TRB(solver, mu0 = kappa, approxParameters = params, - verbosity = verb, homogeneize = homog) + verbosity = verb, homogeneized = homog) else: shift = 10 nsets = 2 stride = 3 Smax = stride * (nsets - 1) + shift + 2 paramsPade = {'S':Smax, 'POD':True, 'sampler':QS([kLeft, kRight], "CHEBYSHEV")} paramsRB = copy(paramsPade) paramsSetsPade = [None] * nsets paramsSetsRB = [None] * nsets for i in range(nsets): paramsSetsPade[i] = {'N': stride*i+shift+1, 'M': stride*i+shift+1, 'S': stride*i+shift+2} paramsSetsRB[i] = {'R': stride*i+shift+2, 'S': stride*i+shift+2} approxPade = LP(solver, mu0 = np.mean([kLeft, kRight]), approxParameters = paramsPade, verbosity = verb, - homogeneize = homog) + homogeneized = homog) approxRB = LRB(solver, mu0 = np.mean([kLeft, kRight]), approxParameters = paramsRB, verbosity = verb, - homogeneize = homog) + homogeneized = homog) sweeper = Sweeper(mutars = ktars, mostExpensive = 'Approx') homogMSG = "" if not homog: homogMSG = "Non" filenamebase = '../data/output/airfoil' + test[:-5] + homogMSG + "Homog" sweeper.ROMEngine = approxPade sweeper.params = paramsSetsPade filenamePade = sweeper.sweep(filenamebase +'Pade.dat') sweeper.ROMEngine = approxRB sweeper.params = paramsSetsRB filenameRB = sweeper.sweep(filenamebase +'RB.dat') if test == "TaylorSweep": constr = ['E'] else: constr = ['S'] sweeper.plotCompare([filenamePade, filenameRB], ['muRe'], ['normHF', 'normApp'], constr, onePlot = True, save = filenamebase + 'Norm', saveFormat = "png", labels = ["Pade'", "RB"]) sweeper.plotCompare([filenamePade, filenameRB], ['muRe'], ['normResRel'], constr, save = filenamebase + 'Res', saveFormat = "png", labels = ["Pade'", "RB"]) sweeper.plotCompare([filenamePade, filenameRB], ['muRe'], ['normErrRel'], constr, save = filenamebase + 'Err', saveFormat = "png", labels = ["Pade'", "RB"]) diff --git a/examples/scipy/laplaceGaussianTaylor.py b/examples/pod/laplaceGaussianTaylor.py similarity index 100% rename from examples/scipy/laplaceGaussianTaylor.py rename to examples/pod/laplaceGaussianTaylor.py diff --git a/examples/scipy/membraneTaylor.py b/examples/pod/membraneTaylor.py similarity index 100% rename from examples/scipy/membraneTaylor.py rename to examples/pod/membraneTaylor.py diff --git a/examples/scipy/parametricDomain.py b/examples/pod/parametricDomain.py similarity index 100% rename from examples/scipy/parametricDomain.py rename to examples/pod/parametricDomain.py diff --git a/examples/scipy/scatteringSquare.py b/examples/pod/scatteringSquare.py similarity index 100% rename from examples/scipy/scatteringSquare.py rename to examples/pod/scatteringSquare.py diff --git a/rrompy/reduction_methods/base/__init__.py b/rrompy/reduction_methods/base/__init__.py index fdc4ada..993f983 100644 --- a/rrompy/reduction_methods/base/__init__.py +++ b/rrompy/reduction_methods/base/__init__.py @@ -1,25 +1,28 @@ # Copyright (C) 2018 by the RROMPy authors # # This file is part of RROMPy. # # RROMPy is free software: you can redistribute it and/or modify # it under the terms of the GNU Lesser General Public License as published by # the Free Software Foundation, either version 3 of the License, or # (at your option) any later version. # # RROMPy is distributed in the hope that it will be useful, # but WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # GNU Lesser General Public License for more details. # # You should have received a copy of the GNU Lesser General Public License # along with RROMPy. If not, see . # from .generic_approximant import GenericApproximant +from .pade_utils import checkRobustTolerance, checkPolyfitRank __all__ = [ - 'GenericApproximant' + 'GenericApproximant', + 'checkRobustTolerance', + 'checkPolyfitRank' ] diff --git a/rrompy/reduction_methods/base/generic_approximant.py b/rrompy/reduction_methods/base/generic_approximant.py index 804e5f6..9b25cae 100644 --- a/rrompy/reduction_methods/base/generic_approximant.py +++ b/rrompy/reduction_methods/base/generic_approximant.py @@ -1,530 +1,530 @@ # Copyright (C) 2018 by the RROMPy authors # # This file is part of RROMPy. # # RROMPy is free software: you can redistribute it and/or modify # it under the terms of the GNU Lesser General Public License as published by # the Free Software Foundation, either version 3 of the License, or # (at your option) any later version. # # RROMPy is distributed in the hope that it will be useful, # but WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # GNU Lesser General Public License for more details. # # You should have received a copy of the GNU Lesser General Public License # along with RROMPy. If not, see . # from abc import abstractmethod import numpy as np from copy import copy from rrompy.sampling.base.sampling_engine_base import SamplingEngineBase from rrompy.utilities.base.types import Np1D, DictAny, HFEng, sampleEng, strLst from rrompy.utilities.base import purgeDict, verbosityDepth __all__ = ['GenericApproximant'] class GenericApproximant: """ ABSTRACT ROM approximant computation for parametric problems. Args: HFEngine: HF problem solver. mu0(optional): Default parameter. Defaults to 0. approxParameters(optional): Dictionary containing values for main parameters of approximant. Recognized keys are: - 'POD': whether to compute POD of snapshots; defaults to True. Defaults to empty dict. verbosity(optional): Verbosity level. Defaults to 10. Attributes: HFEngine: HF problem solver. mu0: Default parameter. w: Weight for norm computation. 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. verbosity: Verbosity level. POD: Whether to compute POD of snapshots. samplingEngine: Sampling engine. 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. lastApproxParameters: List of parameters corresponding to last computed approximant. """ def __init__(self, HFEngine:HFEng, mu0 : complex = 0, - approxParameters : DictAny = {}, homogeneize : bool = False, + approxParameters : DictAny = {}, homogeneized : bool = False, verbosity : int = 10): self._preInit() self.verbosity = verbosity if self.verbosity >= 10: verbosityDepth("INIT", ("Initializing approximant engine of " "type {}.").format(self.name())) self.HFEngine = HFEngine self._HFEngine0 = copy(HFEngine) if not hasattr(self, "parameterList"): self.parameterList = [] self.parameterList += ["POD"] self.mu0 = mu0 - self.homogeneize = homogeneize + self.homogeneized = homogeneized self.approxParameters = approxParameters self._postInit() def _preInit(self): if not hasattr(self, "depth"): self.depth = 0 else: self.depth += 1 def _postInit(self): if self.depth == 0: if self.verbosity >= 10: verbosityDepth("DEL", "Done initializing.\n") del self.depth else: self.depth -= 1 def name(self) -> str: return self.__class__.__name__ def __str__(self) -> str: return self.name() def __repr__(self) -> str: return self.__str__() + " at " + hex(id(self)) def setupSampling(self, SamplingEngine : sampleEng = SamplingEngineBase): """Setup sampling engine.""" self.samplingEngine = SamplingEngine(self.HFEngine, verbosity = self.verbosity) @property def approxParameters(self): """Value of approximant parameters.""" return self._approxParameters @approxParameters.setter def approxParameters(self, approxParams): if not hasattr(self, "approxParameters"): self._approxParameters = {} approxParameters = purgeDict(approxParams, self.parameterList, dictname = self.name() + ".approxParameters", baselevel = 1) keyList = list(approxParameters.keys()) if "POD" in keyList: self.POD = approxParameters["POD"] elif hasattr(self, "POD"): self.POD = self.POD else: self.POD = True @property def POD(self): """Value of POD.""" return self._POD @POD.setter def POD(self, POD): if hasattr(self, "POD"): PODold = self.POD else: PODold = -1 self._POD = POD self._approxParameters["POD"] = self.POD if PODold != self.POD: self.setupSampling() self.resetSamples() @property - def homogeneize(self): - """Value of homogeneize.""" - return self._homogeneize - @homogeneize.setter - def homogeneize(self, homogeneize): - if not hasattr(self, "_homogeneize"): - self._homogeneize = None - if homogeneize != self.homogeneize: - self._homogeneize = homogeneize + def homogeneized(self): + """Value of homogeneized.""" + return self._homogeneized + @homogeneized.setter + def homogeneized(self, homogeneized): + if not hasattr(self, "_homogeneized"): + self._homogeneized = None + if homogeneized != self.homogeneized: + self._homogeneized = homogeneized self.resetSamples() def solveHF(self, mu : complex = None): """ Find high fidelity solution with original parameters and arbitrary parameter. Args: mu: Target parameter. """ if mu is None: mu = self.mu0 if (not hasattr(self, "lastSolvedHF") or not np.isclose(self.lastSolvedHF, mu)): self.uHF = self.samplingEngine.solveLS(mu, - homogeneized = self.homogeneize) + homogeneized = self.homogeneized) self.lastSolvedHF = mu def resetSamples(self): """Reset samples.""" if hasattr(self, "samplingEngine"): self.samplingEngine.resetHistory() else: self.setupSampling() def plotSamples(self, name : str = "u", save : str = None, what : strLst = 'all', saveFormat : str = "eps", saveDPI : int = 100, **figspecs): """ Do some nice plots of the samples. Args: name(optional): Name to be shown as title of the plots. Defaults to 'u'. what(optional): Which plots to do. If list, can contain 'ABS', 'PHASE', 'REAL', 'IMAG'. If str, same plus wildcard 'ALL'. Defaults to 'ALL'. save(optional): Where to save plot(s). Defaults to None, i.e. no saving. saveFormat(optional): Format for saved plot(s). Defaults to "eps". saveDPI(optional): DPI for saved plot(s). Defaults to 100. figspecs(optional key args): Optional arguments for matplotlib figure creation. """ self.samplingEngine.plotSamples(name = name, save = save, what = what, saveFormat = saveFormat, saveDPI = saveDPI, **figspecs) @abstractmethod def setupApprox(self): """ Setup approximant. (ABSTRACT) Any specialization should include something like self.computeDerivatives() if not self.checkComputedApprox(): ... self.lastApproxParameters = copy(self.approxParameters) """ pass def checkComputedApprox(self) -> bool: """ Check if setup of new approximant is not needed. Returns: True if new setup is not needed. False otherwise. """ return (hasattr(self, "lastApproxParameters") and self.approxParameters == self.lastApproxParameters) @abstractmethod def evalApproxReduced(self, mu:complex): """ Evaluate reduced representation of approximant at arbitrary parameter. (ABSTRACT) Any specialization should include something like self.setupApprox() self.uAppReduced = ... Args: mu: Target parameter. """ pass @abstractmethod def evalApprox(self, mu:complex): """ Evaluate approximant at arbitrary parameter. (ABSTRACT) Any specialization should include something like self.evalApproxReduced(mu) self.uApp = ... Args: mu: Target parameter. """ pass def getHF(self, mu:complex, homogeneized : bool = False) -> Np1D: """ Get HF solution at arbitrary parameter. Args: mu: Target parameter. homogeneized(optional): Whether to remove Dirichlet BC. Defaults to False. Returns: HFsolution. """ self.solveHF(mu) - if self.homogeneize and not homogeneized: + if self.homogeneized and not homogeneized: return self.uHF + self.HFEngine.liftDirichletData(mu) - if not self.homogeneize and homogeneized: + if not self.homogeneized and homogeneized: return self.uHF - self.HFEngine.liftDirichletData(mu) return self.uHF def getRHS(self, mu:complex, homogeneized : bool = False) -> Np1D: """ Get linear system RHS at arbitrary parameter. Args: mu: Target parameter. homogeneized(optional): Whether to remove Dirichlet BC. Defaults to False. Returns: Linear system RHS. """ return self.HFEngine.residual(None, mu, homogeneized = homogeneized) def getAppReduced(self, mu:complex) -> Np1D: """ Get approximant at arbitrary parameter. Args: mu: Target parameter. Returns: Reduced approximant. """ self.evalApproxReduced(mu) return self.uAppReduced def getApp(self, mu:complex, homogeneized : bool = False) -> Np1D: """ Get approximant at arbitrary parameter. Args: mu: Target parameter. homogeneized(optional): Whether to remove Dirichlet BC. Defaults to False. Returns: Approximant. """ self.evalApprox(mu) - if self.homogeneize and not homogeneized: + if self.homogeneized and not homogeneized: return self.uApp + self.HFEngine.liftDirichletData(mu) - if not self.homogeneize and homogeneized: + if not self.homogeneized and homogeneized: return self.uApp - self.HFEngine.liftDirichletData(mu) return self.uApp def getRes(self, mu:complex, homogeneized : bool = False) -> Np1D: """ Get residual at arbitrary parameter. Args: mu: Target parameter. homogeneized(optional): Whether to remove Dirichlet BC. Defaults to False. Returns: Approximant residual. """ return self.HFEngine.residual(self.getApp(mu, homogeneized), mu, homogeneized = homogeneized) def getErr(self, mu:complex, homogeneized : bool = False) -> Np1D: """ Get error at arbitrary parameter. Args: mu: Target parameter. homogeneized(optional): Whether to remove Dirichlet BC. Defaults to False. Returns: Approximant error. """ return self.getApp(mu, homogeneized) - self.getHF(mu, homogeneized) @abstractmethod def getPoles(self) -> Np1D: """ Obtain approximant poles. Returns: Numpy complex vector of poles. """ pass def normHF(self, mu:complex, homogeneized : bool = False) -> float: """ Compute norm of HF solution at arbitrary parameter. Args: mu: Target parameter. homogeneized(optional): Whether to remove Dirichlet BC. Defaults to False. Returns: Target norm of HFsolution. """ return self.HFEngine.norm(self.getHF(mu, homogeneized)) def normRHS(self, mu:complex, homogeneized : bool = False) -> Np1D: """ Compute norm of linear system RHS at arbitrary parameter. Args: mu: Target parameter. homogeneized(optional): Whether to remove Dirichlet BC. Defaults to False. Returns: Norm of linear system RHS. """ return self.HFEngine.norm(self.getRHS(mu, homogeneized)) def normApp(self, mu:complex, homogeneized : bool = False) -> float: """ Compute norm of approximant at arbitrary parameter. Args: mu: Target parameter. homogeneized(optional): Whether to remove Dirichlet BC. Defaults to False. Returns: Target norm of approximant. """ return self.HFEngine.norm(self.getApp(mu, homogeneized)) def normRes(self, mu:complex, homogeneized : bool = False) -> float: """ Compute norm of approximant residual at arbitrary parameter. Args: mu: Target parameter. homogeneized(optional): Whether to remove Dirichlet BC. Defaults to False. Returns: Target norm of (A(mu)app(mu) - f(mu)). """ return self.HFEngine.norm(self.getRes(mu, homogeneized)) def normErr(self, mu:complex, homogeneized : bool = False) -> float: """ Compute norm of approximant error at arbitrary parameter. Args: mu: Target parameter. homogeneized(optional): Whether to remove Dirichlet BC. Defaults to False. Returns: Target norm of (approximant - HFsolution). """ return self.HFEngine.norm(self.getErr(mu, homogeneized)) def plotHF(self, mu:complex, name : str = "uHF", save : str = None, what : strLst = 'all', saveFormat : str = "eps", saveDPI : int = 100, homogeneized : bool = False, **figspecs): """ Do some nice plots of the HF solution at arbitrary parameter. Args: mu: Target parameter. name(optional): Name to be shown as title of the plots. Defaults to 'u'. what(optional): Which plots to do. If list, can contain 'ABS', 'PHASE', 'REAL', 'IMAG'. If str, same plus wildcard 'ALL'. Defaults to 'ALL'. save(optional): Where to save plot(s). Defaults to None, i.e. no saving. saveFormat(optional): Format for saved plot(s). Defaults to "eps". saveDPI(optional): DPI for saved plot(s). Defaults to 100. homogeneized(optional): Whether to remove Dirichlet BC. Defaults to False. figspecs(optional key args): Optional arguments for matplotlib figure creation. """ uHF = self.getHF(mu, homogeneized) self.HFEngine.plot(uHF, name = name, save = save, what = what, saveFormat = saveFormat, saveDPI = saveDPI, **figspecs) def plotApp(self, mu:complex, name : str = "uApp", save : str = None, what : strLst = 'all', saveFormat : str = "eps", saveDPI : int = 100, homogeneized : bool = False, **figspecs): """ Do some nice plots of approximant at arbitrary parameter. Args: mu: Target parameter. name(optional): Name to be shown as title of the plots. Defaults to 'u'. what(optional): Which plots to do. If list, can contain 'ABS', 'PHASE', 'REAL', 'IMAG'. If str, same plus wildcard 'ALL'. Defaults to 'ALL'. save(optional): Where to save plot(s). Defaults to None, i.e. no saving. saveFormat(optional): Format for saved plot(s). Defaults to "eps". saveDPI(optional): DPI for saved plot(s). Defaults to 100. homogeneized(optional): Whether to remove Dirichlet BC. Defaults to False. figspecs(optional key args): Optional arguments for matplotlib figure creation. """ uApp = self.getApp(mu, homogeneized) self.HFEngine.plot(uApp, name = name, save = save, what = what, saveFormat = saveFormat, saveDPI = saveDPI, **figspecs) def plotRes(self, mu:complex, name : str = "res", save : str = None, what : strLst = 'all', saveFormat : str = "eps", saveDPI : int = 100, homogeneized : bool = False, **figspecs): """ Do some nice plots of approximation residual at arbitrary parameter. Args: mu: Target parameter. name(optional): Name to be shown as title of the plots. Defaults to 'u'. what(optional): Which plots to do. If list, can contain 'ABS', 'PHASE', 'REAL', 'IMAG'. If str, same plus wildcard 'ALL'. Defaults to 'ALL'. save(optional): Where to save plot(s). Defaults to None, i.e. no saving. saveFormat(optional): Format for saved plot(s). Defaults to "eps". saveDPI(optional): DPI for saved plot(s). Defaults to 100. homogeneized(optional): Whether to remove Dirichlet BC. Defaults to False. figspecs(optional key args): Optional arguments for matplotlib figure creation. """ uRes = self.getRes(mu, homogeneized) self.HFEngine.plot(uRes, name = name, save = save, what = what, saveFormat = saveFormat, saveDPI = saveDPI, **figspecs) def plotErr(self, mu:complex, name : str = "err", save : str = None, what : strLst = 'all', saveFormat : str = "eps", saveDPI : int = 100, homogeneized : bool = False, **figspecs): """ Do some nice plots of approximation error at arbitrary parameter. Args: mu: Target parameter. name(optional): Name to be shown as title of the plots. Defaults to 'u'. what(optional): Which plots to do. If list, can contain 'ABS', 'PHASE', 'REAL', 'IMAG'. If str, same plus wildcard 'ALL'. Defaults to 'ALL'. save(optional): Where to save plot(s). Defaults to None, i.e. no saving. saveFormat(optional): Format for saved plot(s). Defaults to "eps". saveDPI(optional): DPI for saved plot(s). Defaults to 100. homogeneized(optional): Whether to remove Dirichlet BC. Defaults to False. figspecs(optional key args): Optional arguments for matplotlib figure creation. """ uErr = self.getErr(mu, homogeneized) self.HFEngine.plot(uErr, name = name, save = save, what = what, saveFormat = saveFormat, saveDPI = saveDPI, **figspecs) diff --git a/rrompy/reduction_methods/base/pade_utils.py b/rrompy/reduction_methods/base/pade_utils.py new file mode 100644 index 0000000..87744e9 --- /dev/null +++ b/rrompy/reduction_methods/base/pade_utils.py @@ -0,0 +1,46 @@ +# Copyright (C) 2018 by the RROMPy authors +# +# This file is part of RROMPy. +# +# RROMPy is free software: you can redistribute it and/or modify +# it under the terms of the GNU Lesser General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# RROMPy is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU Lesser General Public License for more details. +# +# You should have received a copy of the GNU Lesser General Public License +# along with RROMPy. If not, see . +# + +import numpy as np +from rrompy.utilities.base.types import Np1D, Np2D, Tuple +from rrompy.utilities.warning_manager import warn + +machineEps = np.finfo(float).eps + +def checkRobustTolerance(ev:Np1D, E:int, tol:float) -> dict: + N = len(ev) - 1 + ts = tol * np.linalg.norm(ev) + diff = N - np.sum(np.abs(ev) >= ts) + if diff <= 0: return {} + Enew = E - diff + Nnew = min(N, Enew) + if Nnew == N: + strN = "" + else: + strN = "N from {} to {} and ".format(N, Nnew) + warn(("Smallest {} eigenvalues below tolerance.\nReducing {}E from {} to " + "{}.").format(diff + 1, strN, E, Enew)) + newParameters = {"N" : Nnew, "E" : Enew} + return newParameters + +def checkPolyfitRank(xs:Np1D, ys:Np2D, E:int, rcond:float, + w : Np1D = None) -> Tuple[bool, Np2D]: + G = np.polyfit(xs, ys, deg = E, w = w, rcond = rcond, full = True) + if G[2] < E + 1: return False, G[2] - 1 + return True, G[0] + diff --git a/rrompy/reduction_methods/lagrange/approximant_lagrange_pade.py b/rrompy/reduction_methods/lagrange/approximant_lagrange_pade.py index 132a1d5..f77b827 100644 --- a/rrompy/reduction_methods/lagrange/approximant_lagrange_pade.py +++ b/rrompy/reduction_methods/lagrange/approximant_lagrange_pade.py @@ -1,425 +1,455 @@ # 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 copy import numpy as np +from rrompy.reduction_methods.base import (checkRobustTolerance, + checkPolyfitRank) from .generic_approximant_lagrange import GenericApproximantLagrange from rrompy.utilities.base.types import Np1D, DictAny, List, HFEng from rrompy.utilities.base import purgeDict, verbosityDepth from rrompy.utilities.warning_manager import warn __all__ = ['ApproximantLagrangePade'] class ApproximantLagrangePade(GenericApproximantLagrange): """ ROM Lagrange Pade' 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': whether to compute POD of snapshots; defaults to True; - 'S': total number of samples current approximant relies upon; defaults to 2; - 'sampler': sample point generator; defaults to uniform sampler on [0, 1]; - 'E': coefficient of interpolant to be minimized; defaults to min(S, M + 1); - 'M': degree of Pade' interpolant numerator; defaults to 0; - - 'N': degree of Pade' interpolant denominator; defaults to 0. + - 'N': degree of Pade' interpolant denominator; defaults to 0; + - 'interpRcond': tolerance for interpolation via numpy.polyfit; + defaults to None; + - 'robustTol': tolerance for robust Pade' denominator management; + defaults to 0. Defaults to empty dict. Attributes: HFEngine: HF problem solver. mu0: Default parameter. mus: Array of snapshot parameters. ws: Array of snapshot weigths. 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; - 'S': total number of samples current approximant relies upon; - 'sampler': sample point generator; - 'E': coefficient of interpolant to be minimized; - 'M': degree of Pade' interpolant numerator; - 'N': degree of Pade' interpolant denominator; + - 'interpRcond': tolerance for interpolation via numpy.polyfit; - 'robustTol': tolerance for robust Pade' denominator management. extraApproxParameters: List of approxParameters keys in addition to mother class's. S: Number of solution snapshots over which current approximant is based upon. sampler: Sample point generator. M: Numerator degree of approximant. N: Denominator degree of approximant. POD: Whether to compute POD of snapshots. + interpRcond: Tolerance for interpolation via numpy.polyfit. robustTol: Tolerance for robust Pade' denominator management. samplingEngine: Sampling engine. uHF: High fidelity solution with wavenumber lastSolvedHF as numpy complex vector. lastSolvedHF: Wavenumber corresponding to last computed high fidelity solution. Q: Numpy 1D vector containing complex coefficients of approximant denominator. P: Numpy 2D vector whose columns are FE dofs of coefficients of approximant numerator. uApp: Last evaluated approximant as numpy complex vector. lastApproxParameters: List of parameters corresponding to last computed approximant. """ def __init__(self, HFEngine:HFEng, mu0 : complex = 0., - approxParameters : DictAny = {}, homogeneize : bool = False, + approxParameters : DictAny = {}, homogeneized : bool = False, verbosity : int = 10): self._preInit() if not hasattr(self, "parameterList"): self.parameterList = [] - self.parameterList += ["E", "M", "N", "robustTol"] + self.parameterList += ["E", "M", "N", "interpRcond", "robustTol"] super().__init__(HFEngine = HFEngine, mu0 = mu0, approxParameters = approxParameters, - homogeneize = homogeneize, + homogeneized = homogeneized, verbosity = verbosity) self._postInit() @property def approxParameters(self): """ - Value of approximant parameters. Its assignment may change E, M, N and - S. + Value of approximant parameters. Its assignment may change E, M, N, + robustTol and S. """ return self._approxParameters @approxParameters.setter def approxParameters(self, approxParams): approxParameters = purgeDict(approxParams, self.parameterList, dictname = self.name() + ".approxParameters", baselevel = 1) - approxParametersCopy = purgeDict(approxParameters, ["E", "M", "N"], + approxParametersCopy = purgeDict(approxParameters, ["E", "M", "N", + "interpRcond", + "robustTol"], True, True, baselevel = 1) if hasattr(self, "M"): Mold = self.M self._M = 0 if hasattr(self, "N"): Nold = self.N self._N = 0 if hasattr(self, "E"): self._E = 0 GenericApproximantLagrange.approxParameters.fset(self, approxParametersCopy) keyList = list(approxParameters.keys()) + if "interpRcond" in keyList: + self.interpRcond = approxParameters["interpRcond"] + elif hasattr(self, "interpRcond"): + self.interpRcond = self.interpRcond + else: + self.interpRcond = None if "robustTol" in keyList: self.robustTol = approxParameters["robustTol"] elif hasattr(self, "robustTol"): self.robustTol = self.robustTol else: self.robustTol = 0 if "M" in keyList: self.M = approxParameters["M"] elif hasattr(self, "M"): self.M = Mold else: self.M = 0 if "N" in keyList: self.N = approxParameters["N"] elif hasattr(self, "N"): self.N = Nold else: self.N = 0 if "E" in keyList: self.E = approxParameters["E"] else: self.E = min(self.S - 1, self.M + 1) @property def M(self): """Value of M. Its assignment may change S.""" return self._M @M.setter def M(self, M): if M < 0: raise ArithmeticError("M must be non-negative.") self._M = M self._approxParameters["M"] = self.M if hasattr(self, "S") and self.S < self.M + 1: warn("Prescribed S is too small. Updating S to M + 1.") self.S = self.M + 1 @property def N(self): """Value of N. Its assignment may change S.""" return self._N @N.setter def N(self, N): if N < 0: raise ArithmeticError("N must be non-negative.") self._N = N self._approxParameters["N"] = self.N if hasattr(self, "S") and self.S < self.N + 1: warn("Prescribed S is too small. Updating S to N + 1.") self.S = self.N + 1 @property def E(self): """Value of E. Its assignment may change S.""" return self._E @E.setter def E(self, E): if E < 0: raise ArithmeticError("E must be non-negative.") self._E = E self._approxParameters["E"] = self.E if hasattr(self, "S") and self.S < self.E + 1: warn("Prescribed S is too small. Updating S to E + 1.") self.S = self.E + 1 @property def robustTol(self): """Value of tolerance for robust Pade' denominator management.""" return self._robustTol @robustTol.setter def robustTol(self, robustTol): if robustTol < 0.: warn("Overriding prescribed negative robustness tolerance to 0.") robustTol = 0. self._robustTol = robustTol self._approxParameters["robustTol"] = self.robustTol @property def S(self): """Value of S.""" return self._S @S.setter def S(self, S): if S <= 0: raise ArithmeticError("S must be positive.") if hasattr(self, "S"): Sold = self.S else: Sold = -1 vals, label = [0] * 3, {0:"M", 1:"N", 2:"E"} if hasattr(self, "M"): vals[0] = self.M if hasattr(self, "N"): vals[1] = self.N if hasattr(self, "E"): vals[2] = self.E idxmax = np.argmax(vals) if vals[idxmax] + 1 > S: warn("Prescribed S is too small. Updating S to {} + 1."\ .format(label[idxmax])) self.S = vals[idxmax] + 1 else: self._S = S self._approxParameters["S"] = self.S if Sold != self.S: self.resetSamples() def getTargetFunctionalOptimum(self): """Compute optimum of target LS functional.""" self.setupApprox() - if self._funcOptimumWarn: - warn("Target functional optimum numerically zero.") - return np.abs(self._funcOptimum) + if self._funcOptimumIllCond: + if self.POD: + jOpt = np.linalg.norm(self.P[:, 0]) + else: + jOpt = self.HFEngine.norm(self.samplingEngine.samples.dot( + self.P[:, 0])) + else: + jOpt = np.abs(self._funcOptimum) + return jOpt def setupApprox(self): """ Compute Pade' interpolant. SVD-based robust eigenvalue management. """ if not self.checkComputedApprox(): if self.verbosity >= 5: verbosityDepth("INIT", "Setting up {}.". format(self.name())) self.computeSnapshots() if self.N > 0: if self.verbosity >= 7: verbosityDepth("INIT", ("Starting computation of " "denominator.")) while self.N > 0: TN = np.vander(self.radiusPade(self.mus), N = self.N + 1) if self.POD: data = self.samplingEngine.RPOD else: data = self.samplingEngine.samples RHSFull = np.empty((self.S, data.shape[0] * (self.N + 1)), dtype = np.complex) for j in range(self.S): RHSFull[j, :] = np.kron(TN[j, :], data[:, j]) - G = np.polyfit(self.radiusPade(self.mus), RHSFull, - deg = self.E, w = self.ws) + fullR, G = checkPolyfitRank(self.radiusPade(self.mus), + RHSFull, self.E, w = self.ws, + rcond = self.interpRcond) + if not fullR: + Enew = G + Nnew = min(self.N, Enew) + Mnew = min(self.M, Enew) + if Nnew == self.N: + strN = "" + else: + strN = "N from {} to {} and ".format(self.N, Nnew) + if Mnew == self.M: + strM = "" + else: + strM = "M from {} to {} and ".format(self.M, Mnew) + warn(("Polyfit is poorly conditioned.\nReducing {}{}E " + "from {} to {}.").format(strN, strM, + self.E, Enew)) + newParameters = {"N" : Nnew, "M" : Mnew, "E" : Enew} + self.approxParameters = newParameters + continue G = G[0, :].reshape((self.N + 1, data.shape[0])).T if self.POD: if self.verbosity >= 7: verbosityDepth("INIT", ("Solving svd for square " "root of gramian matrix."), end = "") _, ev, eV = np.linalg.svd(G, full_matrices = False) ev = ev[::-1] eV = eV[::-1, :].conj().T self._funcOptimum = ev[0] else: if self.verbosity >= 10: verbosityDepth("INIT", "Building gramian matrix.", end = "") G2 = self.HFEngine.innerProduct(G, G) if self.verbosity >= 10: verbosityDepth("DEL", "Done building gramian.", inline = True) if self.verbosity >= 7: verbosityDepth("INIT", ("Solving eigenvalue " "problem for gramian " "matrix."), end = "") ev, eV = np.linalg.eigh(G2) self._funcOptimum = ev[0] ** .5 if ev[0] <= 0 or np.isclose(np.abs(ev[0] / ev[-1]), 0.): - self._funcOptimumWarn = True + self._funcOptimumIllCond = True else: - self._funcOptimumWarn = False + self._funcOptimumIllCond = False if self.verbosity >= 7: verbosityDepth("DEL", " Done.", inline = True) - ts = self.robustTol * np.linalg.norm(ev) - Nnew = np.sum(np.abs(ev) >= ts) - diff = self.N - Nnew - if diff <= 0: break - Enew = self.E - diff - Mnew = min(self.M, Enew - 1) - if Mnew == self.M: - strM = "" - else: - strM = ", M from {0} to {1},".format(self.M, Mnew) - print(("Smallest {0} eigenvalues below tolerance.\n" - "Reducing N from {1} to {2}{5} and E from {3} to " - "{4}.").format(diff + 1, self.N, Nnew, - self.E, Enew, strM)) - newParameters = {"N" : Nnew, "M" : Mnew, "E" : Enew} + newParameters = checkRobustTolerance(ev, self.E, + self.robustTol) + if not newParameters: + break self.approxParameters = newParameters - if self.N <= 0: - eV = np.ones((1, 1)) + if self.N <= 0: + eV = np.ones((1, 1)) self.Q = np.poly1d(eV[:, 0]) if self.verbosity >= 7: verbosityDepth("DEL", "Done computing denominator.") else: - if self.POD: - data = self.samplingEngine.RPOD - else: - data = self.samplingEngine.samples - g = np.polyfit(self.radiusPade(self.mus), data.T, - deg = self.E, w = self.ws)[0, :].flatten() - if self.POD: - self._funcOptimum = np.linalg.norm(g) - else: - self._funcOptimum = self.HFEngine.norm(g) - self._funcOptimumWarn = False + self._funcOptimumIllCond = True self.Q = np.poly1d([1]) if self.verbosity >= 7: verbosityDepth("INIT", "Starting computation of numerator.") Qevaldiag = np.diag(self.Q(self.radiusPade(self.mus))) - P = np.polyfit(self.radiusPade(self.mus), Qevaldiag, deg = self.M, - w = self.ws).T + + while self.M >= 0: + fullR, P = checkPolyfitRank(self.radiusPade(self.mus), + Qevaldiag, self.M, w = self.ws, + rcond = self.interpRcond) + if fullR: + P = P.T + break + warn(("Polyfit is poorly conditioned. Reducing M from {} to " + "{}. Exact snapshot interpolation not guaranteed.")\ + .format(self.M, P)) + self.M = P self.P = np.atleast_2d(P) if self.POD: self.P = self.samplingEngine.RPOD.dot(self.P) if self.verbosity >= 7: verbosityDepth("DEL", "Done computing numerator.") self.lastApproxParameters = copy(self.approxParameters) if hasattr(self, "lastSolvedApp"): del self.lastSolvedApp if self.verbosity >= 5: verbosityDepth("DEL", "Done setting up approximant.\n") def radiusPade(self, mu:Np1D, mu0 : float = None) -> float: """ Compute translated radius to be plugged into Pade' approximant. Args: mu: Parameter(s) 1. mu0: Parameter(s) 2. If None, set to self.mu0. Returns: Translated radius to be plugged into Pade' approximant. """ if mu0 is None: mu0 = self.mu0 return self.HFEngine.rescaling(mu) - self.HFEngine.rescaling(mu0) def getPVal(self, mu:List[complex]): """ Evaluate Pade' numerator at arbitrary parameter. Args: mu: Target parameter. """ self.setupApprox() if self.verbosity >= 10: verbosityDepth("INIT", "Evaluating numerator at mu = {}.".format(mu), end = "") try: len(mu) except: mu = [mu] - powerlist = np.vander(self.radiusPade(mu), self.M + 1).T + powerlist = np.vander(self.radiusPade(mu), self.P.shape[1]).T p = self.P.dot(powerlist) if len(mu) == 1: p = p.flatten() if self.verbosity >= 10: verbosityDepth("DEL", " Done.", inline = True) return p def getQVal(self, mu:List[complex]): """ Evaluate Pade' denominator at arbitrary parameter. Args: mu: Target parameter. """ self.setupApprox() if self.verbosity >= 10: verbosityDepth("INIT", "Evaluating denominator at mu = {}.".format(mu), end = "") q = self.Q(self.radiusPade(mu)) if self.verbosity >= 10: verbosityDepth("DEL", " Done.", inline = True) return q def evalApproxReduced(self, mu:complex): """ Evaluate Pade' approximant at arbitrary parameter. Args: mu: Target parameter. """ self.setupApprox() if (not hasattr(self, "lastSolvedApp") or not np.isclose(self.lastSolvedApp, mu)): if self.verbosity >= 5: verbosityDepth("INIT", "Evaluating approximant at mu = {}.".format(mu)) self.uAppReduced = self.getPVal(mu) / self.getQVal(mu) self.lastSolvedApp = mu if self.verbosity >= 5: verbosityDepth("DEL", "Done evaluating approximant.") def evalApprox(self, mu:complex): """ Evaluate approximant at arbitrary parameter. Args: mu: Target parameter. """ self.evalApproxReduced(mu) self.uApp = self.samplingEngine.samples.dot(self.uAppReduced) def getPoles(self) -> Np1D: """ Obtain approximant poles. Returns: Numpy complex vector of poles. """ self.setupApprox() return self.HFEngine.rescalingInv(self.Q.r + self.HFEngine.rescaling(self.mu0)) diff --git a/rrompy/reduction_methods/lagrange/approximant_lagrange_rb.py b/rrompy/reduction_methods/lagrange/approximant_lagrange_rb.py index 8d0751f..d9a725b 100644 --- a/rrompy/reduction_methods/lagrange/approximant_lagrange_rb.py +++ b/rrompy/reduction_methods/lagrange/approximant_lagrange_rb.py @@ -1,292 +1,291 @@ # 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 copy import numpy as np import scipy as sp from .generic_approximant_lagrange import GenericApproximantLagrange from rrompy.utilities.base.types import Np1D, DictAny, HFEng from rrompy.utilities.base import purgeDict, verbosityDepth from rrompy.utilities.warning_manager import warn __all__ = ['ApproximantLagrangeRB'] class ApproximantLagrangeRB(GenericApproximantLagrange): """ ROM RB approximant computation for parametric problems. Args: HFEngine: HF problem solver. mu0(optional): Default parameter. Defaults to 0. approxParameters(optional): Dictionary containing values for main parameters of approximant. Recognized keys are: - 'POD': whether to compute POD of snapshots; defaults to True; - 'S': total number of samples current approximant relies upon; defaults to 2; - 'sampler': sample point generator; defaults to uniform sampler on [0, 1]; - 'R': rank for Galerkin projection; defaults to S. Defaults to empty dict. Attributes: HFEngine: HF problem solver. mu0: Default parameter. mus: Array of snapshot parameters. ws: Array of snapshot weigths (unused). approxRadius: Dummy radius of approximant (i.e. distance from mu0 to farthest sample point). approxParameters: Dictionary containing values for main parameters of approximant. Recognized keys are in parameterList. parameterList: Recognized keys of approximant parameters: - 'POD': whether to compute POD of snapshots; - 'S': total number of samples current approximant relies upon; - 'sampler': sample point generator; - 'R': rank for Galerkin projection. extraApproxParameters: List of approxParameters keys in addition to mother class's. S: Number of solution snapshots over which current approximant is based upon. sampler: Sample point generator. R: Rank for Galerkin projection. POD: Whether to compute POD of snapshots. samplingEngine: Sampling engine. projMat: Projection matrix for RB system assembly. 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. lastApproxParameters: List of parameters corresponding to last computed approximant. As: List of sparse matrices (in CSC format) representing coefficients of linear system matrix wrt theta(mu). bs: List of numpy vectors representing coefficients of linear system RHS wrt theta(mu). thetaAs: List of callables representing coefficients of linear system matrix wrt mu. thetabs: List of callables representing coefficients of linear system RHS wrt mu. ARBs: List of sparse matrices (in CSC format) representing coefficients of compressed linear system matrix wrt theta(mu). bRBs: List of numpy vectors representing coefficients of compressed linear system RHS wrt theta(mu). """ def __init__(self, HFEngine:HFEng, mu0 : complex = 0., - approxParameters : DictAny = {}, homogeneize : bool = False, + approxParameters : DictAny = {}, homogeneized : bool = False, verbosity : int = 10): self._preInit() if not hasattr(self, "parameterList"): self.parameterList = [] self.parameterList += ["R"] super().__init__(HFEngine = HFEngine, mu0 = mu0, approxParameters = approxParameters, - homogeneize = homogeneize, + homogeneized = homogeneized, verbosity = verbosity) if self.verbosity >= 10: verbosityDepth("INIT", "Computing affine blocks of system.") self.As, self.thetaAs = self.HFEngine.affineBlocksA(self.mu0) self.bs, self.thetabs = self.HFEngine.affineBlocksb(self.mu0, - self.homogeneize) + self.homogeneized) if self.verbosity >= 10: verbosityDepth("DEL", "Done computing affine blocks.") self._postInit() def resetSamples(self): """Reset samples.""" super().resetSamples() self.projMat = None @property def approxParameters(self): """ Value of approximant parameters. Its assignment may change M, N and S. """ return self._approxParameters @approxParameters.setter def approxParameters(self, approxParams): approxParameters = purgeDict(approxParams, self.parameterList, dictname = self.name() + ".approxParameters", baselevel = 1) approxParametersCopy = purgeDict(approxParameters, ["R"], True, True, baselevel = 1) GenericApproximantLagrange.approxParameters.fset(self, approxParametersCopy) keyList = list(approxParameters.keys()) if "R" in keyList: self.R = approxParameters["R"] elif hasattr(self, "R"): self.R = self.R else: self.R = self.S @property def R(self): """Value of R. Its assignment may change S.""" return self._R @R.setter def R(self, R): if R < 0: raise ArithmeticError("R must be non-negative.") self._R = R self._approxParameters["R"] = self.R if hasattr(self, "S") and self.S < self.R: warn("Prescribed S is too small. Updating S to R.") self.S = self.R def checkComputedApprox(self) -> bool: """ Check if setup of new approximant is not needed. Returns: True if new setup is not needed. False otherwise. """ return (self.projMat is not None and super().checkComputedApprox()) def setupApprox(self): """Compute RB projection matrix.""" if not self.checkComputedApprox(): if self.verbosity >= 5: verbosityDepth("INIT", "Setting up {}.". format(self.name())) self.computeSnapshots() if self.verbosity >= 7: verbosityDepth("INIT", "Computing projection matrix.", end = "") if self.POD: U, _, _ = np.linalg.svd(self.samplingEngine.RPOD, full_matrices = False) self.projMat = self.samplingEngine.samples.dot(U[:, : self.R]) else: self.projMat = self.samplingEngine.samples[:, : self.R] if self.verbosity >= 7: verbosityDepth("DEL", " Done.", inline = True) self.lastApproxParameters = copy(self.approxParameters) if hasattr(self, "lastSolvedApp"): del self.lastSolvedApp - self.assembleReducedSystem() if self.verbosity >= 5: verbosityDepth("DEL", "Done setting up approximant.\n") def assembleReducedSystem(self): """Build affine blocks of RB linear system through projections.""" self.setupApprox() if self.verbosity >= 10: verbosityDepth("INIT", "Projecting affine terms of HF model.", end = "") - projMatH = self.projMat.T.conjugate() + projMatH = self.projMat.T.conj() self.ARBs = [None] * len(self.As) self.bRBs = [None] * len(self.bs) if self.verbosity >= 10: verbosityDepth("MAIN", ".", end = "", inline = True) for j in range(len(self.As)): self.ARBs[j] = projMatH.dot(self.As[j].dot(self.projMat)) if self.verbosity >= 10: verbosityDepth("MAIN", ".", end = "", inline = True) for j in range(len(self.bs)): self.bRBs[j] = projMatH.dot(self.bs[j]) if self.verbosity >= 10: verbosityDepth("DEL", "Done.", inline = True) def solveReducedSystem(self, mu:complex) -> Np1D: """ Solve RB linear system. Args: mu: Target parameter. Returns: Solution of RB linear system. """ - self.setupApprox() + self.assembleReducedSystem() if self.verbosity >= 10: verbosityDepth("INIT", "Assembling reduced model for mu = {}.".format(mu), end = "") ARBmu = self.thetaAs(mu, 0) * self.ARBs[0][:self.R,:self.R] bRBmu = self.thetabs(mu, 0) * self.bRBs[0][:self.R] if self.verbosity >= 10: verbosityDepth("MAIN", ".", end = "", inline = True) for j in range(1, len(self.ARBs)): ARBmu += self.thetaAs(mu, j) * self.ARBs[j][:self.R, :self.R] if self.verbosity >= 10: verbosityDepth("MAIN", ".", end = "", inline = True) for j in range(1, len(self.bRBs)): bRBmu += self.thetabs(mu, j) * self.bRBs[j][:self.R] if self.verbosity >= 10: verbosityDepth("DEL", "Done.", inline = True) if self.verbosity >= 5: verbosityDepth("INIT", "Solving reduced model for mu = {}.".format(mu), end = "") uRB = np.linalg.solve(ARBmu, bRBmu) if self.verbosity >= 5: verbosityDepth("DEL", " Done.", inline = True) return uRB def evalApproxReduced(self, mu:complex): """ Evaluate RB approximant at arbitrary wavenumber. Args: mu: Target parameter. """ - self.setupApprox() + self.assembleReducedSystem() if (not hasattr(self, "lastSolvedApp") or not np.isclose(self.lastSolvedApp, mu)): if self.verbosity >= 5: verbosityDepth("INIT", "Computing RB solution at mu = {}.".format(mu)) self.uAppReduced = self.solveReducedSystem(mu) self.lastSolvedApp = mu if self.verbosity >= 5: verbosityDepth("DEL", "Done computing RB solution.") def evalApprox(self, mu:complex): """ Evaluate approximant at arbitrary parameter. Args: mu: Target parameter. """ self.evalApproxReduced(mu) self.uApp = self.projMat[:, :self.R].dot(self.uAppReduced) def getPoles(self) -> Np1D: """ Obtain approximant poles. Returns: Numpy complex vector of poles. """ warn(("Impossible to compute poles in general affine parameter " "dependence. Results subject to interpretation/rescaling, or " "possibly completely wrong.")) - self.setupApprox() + self.assembleReducedSystem() if len(self.ARBs) < 2: return A = np.eye(self.ARBs[0].shape[0] * (len(self.ARBs) - 1), dtype = np.complex) B = np.zeros_like(A) A[: self.ARBs[0].shape[0], : self.ARBs[0].shape[0]] = - self.ARBs[0] for j in range(len(self.ARBs) - 1): Aj = self.ARBs[j + 1] B[: Aj.shape[0], j * Aj.shape[0] : (j + 1) * Aj.shape[0]] = Aj II = np.arange(self.ARBs[0].shape[0], self.ARBs[0].shape[0] * (len(self.ARBs) - 1)) B[II, II - self.ARBs[0].shape[0]] = 1. return self.HFEngine.rescalingInv(sp.linalg.eigvals(A, B) + self.HFEngine.rescaling(self.mu0)) diff --git a/rrompy/reduction_methods/lagrange/generic_approximant_lagrange.py b/rrompy/reduction_methods/lagrange/generic_approximant_lagrange.py index 8c3c94b..00756f3 100644 --- a/rrompy/reduction_methods/lagrange/generic_approximant_lagrange.py +++ b/rrompy/reduction_methods/lagrange/generic_approximant_lagrange.py @@ -1,217 +1,202 @@ # Copyright (C) 2018 by the RROMPy authors # # This file is part of RROMPy. # # RROMPy is free software: you can redistribute it and/or modify # it under the terms of the GNU Lesser General Public License as published by # the Free Software Foundation, either version 3 of the License, or # (at your option) any later version. # # RROMPy is distributed in the hope that it will be useful, # but WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # GNU Lesser General Public License for more details. # # You should have received a copy of the GNU Lesser General Public License # along with RROMPy. If not, see . # import numpy as np from rrompy.reduction_methods.base.generic_approximant import ( GenericApproximant) from rrompy.utilities.base.types import DictAny, HFEng from rrompy.utilities.base import purgeDict, verbosityDepth __all__ = ['GenericApproximantLagrange'] class GenericApproximantLagrange(GenericApproximant): """ ROM Lagrange 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; - 'S': total number of samples current approximant relies upon; - 'sampler': sample point generator; defaults to uniform sampler on [0, 1]. Defaults to empty dict. Attributes: HFEngine: HF problem solver. mu0: Default parameter. mus: Array of snapshot parameters. ws: Array of snapshot weigths. 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; - 'S': total number of snapshots current approximant relies upon; - 'sampler': sample point generator. extraApproxParameters: List of approxParameters keys in addition to mother class's. S: Number of solution snapshots over which current approximant is based upon. sampler: Sample point generator. POD: Whether to compute POD of snapshots. samplingEngine: Sampling engine. 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. lastApproxParameters: List of parameters corresponding to last computed approximant. """ def __init__(self, HFEngine:HFEng, mu0 : complex = 0., - approxParameters : DictAny = {}, homogeneize : bool = False, + approxParameters : DictAny = {}, homogeneized : bool = False, verbosity : int = 10): self._preInit() if not hasattr(self, "parameterList"): self.parameterList = [] self.parameterList += ["S", "sampler"] super().__init__(HFEngine = HFEngine, mu0 = mu0, approxParameters = approxParameters, - homogeneize = homogeneize, + homogeneized = homogeneized, verbosity = verbosity) self._postInit() def setupSampling(self): """Setup sampling engine.""" if not hasattr(self, "POD"): return if self.POD: from rrompy.sampling.scipy.sampling_engine_lagrange_pod import ( SamplingEngineLagrangePOD) super().setupSampling(SamplingEngineLagrangePOD) else: from rrompy.sampling.scipy.sampling_engine_lagrange import ( SamplingEngineLagrange) super().setupSampling(SamplingEngineLagrange) @property def mus(self): """Value of mus. Its assignment may reset snapshots.""" return self._mus @mus.setter def mus(self, mus): musOld = self.mus if hasattr(self, 'mus') else None self._mus = np.array(mus) _, musCounts = np.unique(self._mus, return_counts = True) if len(np.where(musCounts > 1)[0]) > 0: raise Exception("Repeated sample points not allowed.") if (musOld is None or len(self.mus) != len(musOld) or not np.allclose(self.mus, musOld, 1e-14)): self.resetSamples() self.autoNode = None @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, ["S", "sampler"], True, True, baselevel = 1) GenericApproximant.approxParameters.fset(self, approxParametersCopy) keyList = list(approxParameters.keys()) if "S" in keyList: self.S = approxParameters["S"] elif hasattr(self, "S"): self.S = self.S else: self.S = 2 if "sampler" in keyList: self.sampler = approxParameters["sampler"] elif not hasattr(self, "S"): from rrompy.utilities.parameter_sampling import QuadratureSampler self.sampler = QuadratureSampler([0., 1.], "UNIFORM") del QuadratureSampler - - @property - def POD(self): - """Value of POD.""" - return self._POD - @POD.setter - def POD(self, POD): - if hasattr(self, "POD"): PODold = self.POD - else: PODold = -1 - self._POD = POD - self._approxParameters["POD"] = self.POD - if PODold != self.POD: - self.setupSampling() - self.resetSamples() @property def S(self): """Value of S.""" return self._S @S.setter def S(self, S): if S <= 0: raise ArithmeticError("S must be positive.") if hasattr(self, "S"): Sold = self.S else: Sold = -1 self._S = S self._approxParameters["S"] = self.S if Sold != self.S: self.resetSamples() @property def sampler(self): """Value of sampler.""" return self._sampler @sampler.setter def sampler(self, sampler): if 'generatePoints' not in dir(sampler): raise Exception("Sampler type not recognized.") if hasattr(self, '_sampler'): samplerOld = self.sampler self._sampler = sampler self._approxParameters["sampler"] = self.sampler if not 'samplerOld' in locals() or samplerOld != self.sampler: self.resetSamples() def computeSnapshots(self): - """ - Compute snapshots of solution map. - """ + """Compute snapshots of solution map.""" if self.samplingEngine.samples is None: if self.verbosity >= 5: verbosityDepth("INIT", "Starting computation of snapshots.") self.mus, self.ws = self.sampler.generatePoints(self.S) self.mus = np.array([x[0] for x in self.mus]) self.samplingEngine.iterSample(self.mus, - homogeneize = self.homogeneize) + homogeneized = self.homogeneized) if self.verbosity >= 5: verbosityDepth("DEL", "Done computing snapshots.") def checkComputedApprox(self) -> bool: """ Check if setup of new approximant is not needed. Returns: True if new setup is not needed. False otherwise. """ return (self.samplingEngine.samples is not None and super().checkComputedApprox()) def normApp(self, mu:complex, homogeneized : bool = False) -> float: """ Compute norm of approximant at arbitrary parameter. Args: mu: Target parameter. homogeneized(optional): Whether to remove Dirichlet BC. Defaults to False. Returns: Target norm of approximant. """ - if not self.POD or self.homogeneize != homogeneized: + if not self.POD or self.homogeneized != homogeneized: return super().normApp(mu, homogeneized) return np.linalg.norm(self.getAppReduced(mu)) + diff --git a/rrompy/reduction_methods/hermite/__init__.py b/rrompy/reduction_methods/lagrange_greedy/__init__.py similarity index 62% rename from rrompy/reduction_methods/hermite/__init__.py rename to rrompy/reduction_methods/lagrange_greedy/__init__.py index ed60590..f926ff6 100644 --- a/rrompy/reduction_methods/hermite/__init__.py +++ b/rrompy/reduction_methods/lagrange_greedy/__init__.py @@ -1,18 +1,30 @@ # Copyright (C) 2018 by the RROMPy authors # # This file is part of RROMPy. # # RROMPy is free software: you can redistribute it and/or modify # it under the terms of the GNU Lesser General Public License as published by # the Free Software Foundation, either version 3 of the License, or # (at your option) any later version. # # RROMPy is distributed in the hope that it will be useful, # but WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # GNU Lesser General Public License for more details. # # You should have received a copy of the GNU Lesser General Public License # along with RROMPy. If not, see . # +from .generic_approximant_lagrange_greedy import ( + GenericApproximantLagrangeGreedy) +from .approximant_lagrange_greedy_pade import ApproximantLagrangePadeGreedy +from .approximant_lagrange_greedy_rb import ApproximantLagrangeRBGreedy + +__all__ = [ + 'GenericApproximantLagrangeGreedy', + 'ApproximantLagrangePadeGreedy', + 'ApproximantLagrangeRBGreedy' + ] + + diff --git a/rrompy/reduction_methods/lagrange_greedy/approximant_lagrange_greedy_pade.py b/rrompy/reduction_methods/lagrange_greedy/approximant_lagrange_greedy_pade.py new file mode 100644 index 0000000..7b5cdbd --- /dev/null +++ b/rrompy/reduction_methods/lagrange_greedy/approximant_lagrange_greedy_pade.py @@ -0,0 +1,294 @@ +# 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 copy +import numpy as np +from rrompy.reduction_methods.base import (checkRobustTolerance, + checkPolyfitRank) +from .generic_approximant_lagrange_greedy import ( + GenericApproximantLagrangeGreedy) +from rrompy.reduction_methods.lagrange import ApproximantLagrangePade +from rrompy.utilities.base.types import DictAny, List, HFEng +from rrompy.utilities.base import purgeDict, verbosityDepth +from rrompy.utilities.warning_manager import warn + +__all__ = ['ApproximantLagrangePadeGreedy'] + +class ApproximantLagrangePadeGreedy(GenericApproximantLagrangeGreedy, + ApproximantLagrangePade): + """ + ROM greedy Lagrange Pade' 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': whether to compute POD of snapshots; defaults to True; + - 'muBounds': list of bounds for parameter values; defaults to + [[0], [1]]; + - 'greedyTol': uniform error tolerance for greedy algorithm; + defaults to 1e-2; + - 'maxIter': maximum number of greedy steps; defaults to 1e2; + - 'refinementRatio': ratio of training points to be exhausted + before training set refinement; defaults to 0.2; + - 'nTrainingPoints': number of training points; defaults to + maxIter / refinementRatio; + - 'nTestPoints': number of starting test points; defaults to 1; + - 'trainingSetGenerator': training sample points generator; + defaults to uniform sampler within muBounds; + - 'interpRcond': tolerance for interpolation via numpy.polyfit; + defaults to None; + - 'robustTol': tolerance for robust Pade' denominator management; + defaults to 0. + Defaults to empty dict. + + 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. + parameterList: Recognized keys of approximant parameters: + - 'POD': whether to compute POD of snapshots; + - 'muBounds': list of bounds for parameter values; + - '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; + - 'nTrainingPoints': number of training points; + - 'nTestPoints': number of starting test points; + - 'trainingSetGenerator': training sample points generator; + - 'interpRcond': tolerance for interpolation via numpy.polyfit; + - 'robustTol': tolerance for robust Pade' denominator management. + 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. + 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. + nTrainingPoints: number of training points. + nTestPoints: number of starting test points. + trainingSetGenerator: training sample points generator. + interpRcond: Tolerance for interpolation via numpy.polyfit. + robustTol: Tolerance for robust Pade' denominator management. + samplingEngine: Sampling engine. + uHF: High fidelity solution with wavenumber lastSolvedHF as numpy + complex vector. + lastSolvedHF: Wavenumber corresponding to last computed high fidelity + solution. + Q: Numpy 1D vector containing complex coefficients of approximant + denominator. + P: Numpy 2D vector whose columns are FE dofs of coefficients of + approximant numerator. + uApp: Last evaluated approximant as numpy complex vector. + lastApproxParameters: List of parameters corresponding to last + computed approximant. + """ + + def __init__(self, HFEngine:HFEng, mu0 : complex = 0., + approxParameters : DictAny = {}, homogeneized : bool = False, + verbosity : int = 10): + self._preInit() + if not hasattr(self, "parameterList"): + self.parameterList = [] + self.parameterList += ["interpRcond", "robustTol"] + super().__init__(HFEngine = HFEngine, mu0 = mu0, + approxParameters = approxParameters, + homogeneized = homogeneized, + verbosity = verbosity) + self._postInit() + + @property + def approxParameters(self): + """ + Value of approximant parameters. Its assignment may change robustTol. + """ + return self._approxParameters + @approxParameters.setter + def approxParameters(self, approxParams): + approxParameters = purgeDict(approxParams, self.parameterList, + dictname = self.name() + ".approxParameters", + baselevel = 1) + approxParametersCopy = purgeDict(approxParameters, ["interpRcond", + "robustTol"], + True, True, baselevel = 1) + GenericApproximantLagrangeGreedy.approxParameters.fset(self, + approxParametersCopy) + keyList = list(approxParameters.keys()) + if "interpRcond" in keyList: + self.interpRcond = approxParameters["interpRcond"] + elif hasattr(self, "interpRcond"): + self.interpRcond = self.interpRcond + else: + self.interpRcond = None + if "robustTol" in keyList: + self.robustTol = approxParameters["robustTol"] + elif hasattr(self, "robustTol"): + self.robustTol = self.robustTol + else: + self.robustTol = 0 + + def errorEstimator(self, mus:List[np.complex], + homogeneized : bool = None) -> List[np.complex]: + """ + Standard residual-based error estimator. Unreliable for unstable + problems. + """ + if homogeneized is None: + homogeneized = self.homogeneized + self.setupApprox() + self.initMassMatrixInverse() + nmus = len(mus) + jOpt = self.getTargetFunctionalOptimum() + if np.isnan(jOpt): + err = np.empty_like(mus) + err[:] = np.inf + return err + musTile = np.tile(self.HFEngine.rescaling(mus).reshape(-1, 1), + [1, len(self.mus)]) + smusCol = self.HFEngine.rescaling(self.mus).reshape(1, -1) + num = np.abs(np.prod(musTile - smusCol, axis = 1)) + den = np.abs(self.getQVal(mus)) + if self.HFEngine.nbs == 1: + RHS = self.getRHS(mus[0], homogeneized = homogeneized) + err = (jOpt * num / den + / self.HFEngine.norm(self.massInv.solve(RHS))) + else: + err = [0.] * nmus + for j in range(nmus): + RHS = self.getRHS(mus[j], homogeneized = homogeneized) + err[j] = (jOpt * num[j] / den[j] + / self.HFEngine.norm(self.massInv.solve(RHS))) + return err + + def setupApprox(self): + """ + Compute Pade' interpolant. + SVD-based robust eigenvalue management. + """ + if not self.checkComputedApprox(): + if self.verbosity >= 5: + verbosityDepth("INIT", "Setting up {}.". format(self.name())) + self.S = len(self.mus) + M = self.S - 1 + N = self.S - 1 + self.greedy() + if N > 0: + if self.verbosity >= 7: + verbosityDepth("INIT", ("Starting computation of " + "denominator.")) + while N > 0: + TN = np.vander(self.radiusPade(self.mus), N = N + 1) + if self.POD: + data = self.samplingEngine.RPOD + else: + data = self.samplingEngine.samples + RHSFull = np.empty((self.S, data.shape[0] * (N + 1)), + dtype = np.complex) + for j in range(self.S): + RHSFull[j, :] = np.kron(data[:, j], TN[j, :]) + fullR, G = checkPolyfitRank(self.radiusPade(self.mus), + RHSFull, M, + rcond = self.interpRcond) + if not fullR: + warn(("Polyfit is poorly conditioned. Starting " + "preemptive termination of computation of " + "approximant.")) + self._funcOptimum = np.nan + self._funcOptimumIllCond = False + self.Q = np.poly1d([1.]) + self.P = np.diag([np.nan] * len(self.mus)) + self.lastApproxParameters = copy(self.approxParameters) + if hasattr(self, "lastSolvedApp"): + del self.lastSolvedApp + if self.verbosity >= 7: + verbosityDepth("DEL", ("Aborting computation of " + "denominator.")) + if self.verbosity >= 5: + verbosityDepth("DEL", ("Aborting computation of " + "approximant.\n")) + return + G = G[0, :].reshape((data.shape[0], N + 1)) + if self.POD: + if self.verbosity >= 7: + verbosityDepth("INIT", ("Solving svd for square " + "root of gramian matrix."), + end = "") + _, ev, eV = np.linalg.svd(G, full_matrices = False) + ev = ev[::-1] + eV = eV[::-1, :].conj().T + self._funcOptimum = ev[0] + else: + if self.verbosity >= 10: + verbosityDepth("INIT", "Building gramian matrix.", + end = "") + G2 = self.HFEngine.innerProduct(G, G) + if self.verbosity >= 10: + verbosityDepth("DEL", "Done building gramian.", + inline = True) + if self.verbosity >= 7: + verbosityDepth("INIT", ("Solving eigenvalue " + "problem for gramian " + "matrix."), end = "") + ev, eV = np.linalg.eigh(G2) + self._funcOptimum = ev[0] ** .5 + if ev[0] <= 0 or np.isclose(np.abs(ev[0] / ev[-1]), 0.): + self._funcOptimumIllCond = True + else: + self._funcOptimumIllCond = False + if self.verbosity >= 7: + verbosityDepth("DEL", " Done.", inline = True) + newParameters = checkRobustTolerance(ev, M, self.robustTol) + if not newParameters: + break + N = newParameters["N"] + M = newParameters["E"] + if N <= 0: + eV = np.ones((1, 1)) + self.Q = np.poly1d(eV[:, 0]) + if self.verbosity >= 7: + verbosityDepth("DEL", "Done computing denominator.") + else: + self._funcOptimumIllCond = True + self.Q = np.poly1d([1]) + if self.verbosity >= 7: + verbosityDepth("INIT", "Starting computation of numerator.") + Qevaldiag = np.diag(self.Q(self.radiusPade(self.mus))) + while M >= 0: + fullR, P = checkPolyfitRank(self.radiusPade(self.mus), + Qevaldiag, M, + rcond = self.interpRcond) + if fullR: + P = P.T + break + warn(("Polyfit is poorly conditioned. Reducing M from {} to " + "{}. Exact snapshot interpolation not guaranteed.")\ + .format(M, P)) + M = P + self.P = np.atleast_2d(P) + if self.POD: + self.P = self.samplingEngine.RPOD.dot(self.P) + if self.verbosity >= 7: + verbosityDepth("DEL", "Done computing numerator.") + self.lastApproxParameters = copy(self.approxParameters) + if hasattr(self, "lastSolvedApp"): del self.lastSolvedApp + if self.verbosity >= 5: + verbosityDepth("DEL", "Done setting up approximant.\n") diff --git a/rrompy/reduction_methods/lagrange_greedy/approximant_lagrange_greedy_rb.py b/rrompy/reduction_methods/lagrange_greedy/approximant_lagrange_greedy_rb.py new file mode 100644 index 0000000..7781de6 --- /dev/null +++ b/rrompy/reduction_methods/lagrange_greedy/approximant_lagrange_greedy_rb.py @@ -0,0 +1,319 @@ +# Copyright (C) 2018 by the RROMPy authors +# +# This file is part of RROMPy. +# +# RROMPy is free software: you can redistribute it and/or modify +# it under the terms of the GNU Lesser General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# RROMPy is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU Lesser General Public License for more details. +# +# You should have received a copy of the GNU Lesser General Public License +# along with RROMPy. If not, see . +# + +import numpy as np +from copy import copy +from .generic_approximant_lagrange_greedy import ( + GenericApproximantLagrangeGreedy) +from rrompy.reduction_methods.lagrange import ApproximantLagrangeRB +from rrompy.utilities.base.types import DictAny, HFEng, List +from rrompy.utilities.base import verbosityDepth +from rrompy.utilities.warning_manager import warn + +__all__ = ['ApproximantLagrangeRBGreedy'] + +class ApproximantLagrangeRBGreedy(GenericApproximantLagrangeGreedy, + ApproximantLagrangeRB): + """ + ROM greedy RB approximant computation for parametric problems. + + Args: + HFEngine: HF problem solver. + mu0(optional): Default parameter. Defaults to 0. + approxParameters(optional): Dictionary containing values for main + parameters of approximant. Recognized keys are: + - 'POD': whether to compute POD of snapshots; defaults to True; + - 'muBounds': list of bounds for parameter values; defaults to + [[0], [1]]; + - 'greedyTol': uniform error tolerance for greedy algorithm; + defaults to 1e-2; + - 'maxIter': maximum number of greedy steps; defaults to 1e2; + - 'refinementRatio': ratio of training points to be exhausted + before training set refinement; defaults to 0.2; + - 'nTrainingPoints': number of training points; defaults to + maxIter / refinementRatio; + - 'nTestPoints': number of starting test points; defaults to 1; + - 'trainingSetGenerator': training sample points generator; + defaults to uniform sampler within muBounds. + Defaults to empty dict. + + 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. + parameterList: Recognized keys of approximant parameters: + - 'POD': whether to compute POD of snapshots; + - 'muBounds': list of bounds for parameter values; + - '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; + - 'nTrainingPoints': number of training points; + - 'nTestPoints': number of starting test points; + - 'trainingSetGenerator': training sample points generator; + - 'robustTol': tolerance for robust Pade' denominator management. + 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. + 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. + nTrainingPoints: number of training points. + nTestPoints: number of starting test points. + trainingSetGenerator: training sample points generator. + samplingEngine: Sampling engine. + projMat: Projection matrix for RB system assembly. + 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. + lastApproxParameters: List of parameters corresponding to last + computed approximant. + As: List of sparse matrices (in CSC format) representing coefficients + of linear system matrix wrt theta(mu). + bs: List of numpy vectors representing coefficients of linear system + RHS wrt theta(mu). + thetaAs: List of callables representing coefficients of linear system + matrix wrt mu. + thetabs: List of callables representing coefficients of linear system + RHS wrt mu. + ARBs: List of sparse matrices (in CSC format) representing coefficients + of compressed linear system matrix wrt theta(mu). + bRBs: List of numpy vectors representing coefficients of compressed + linear system RHS wrt theta(mu). + """ + + def __init__(self, HFEngine:HFEng, mu0 : complex = 0., + approxParameters : DictAny = {}, homogeneized : bool = False, + verbosity : int = 10): + self._preInit() + super().__init__(HFEngine = HFEngine, mu0 = mu0, + approxParameters = approxParameters, + homogeneized = homogeneized, + verbosity = verbosity) + if self.verbosity >= 10: + verbosityDepth("INIT", "Computing affine blocks of system.") + self.As, self.thetaAs = self.HFEngine.affineBlocksA(self.mu0) + self.bs, self.thetabs = self.HFEngine.affineBlocksb(self.mu0, + self.homogeneized) + if self.verbosity >= 10: + verbosityDepth("DEL", "Done computing affine blocks.") + self._postInit() + + @property + def S(self): + """Value of S.""" + return self._S + @S.setter + def S(self, S): + self._S = S + + @property + def R(self): + """Value of R.""" + return self._S + @R.setter + def R(self, R): + warn(("R is used just to simplify inheritance, and its value cannot " + "be changed from that of S.")) + + def resetSamples(self): + """Reset samples.""" + super().resetSamples() + self.projMat = None + self.resbb = None + self.resAb = None + self.resAA = None + + def checkComputedApprox(self) -> bool: + """ + Check if setup of new approximant is not needed. + + Returns: + True if new setup is not needed. False otherwise. + """ + return (self.projMat is not None and super().checkComputedApprox()) + + def errorEstimator(self, mus:List[np.complex], + homogeneized : bool = None) -> List[np.complex]: + """ + Standard residual-based error estimator. Unreliable for unstable + problems (inf-sup constant is missing). + """ + if homogeneized is None: + homogeneized = self.homogeneized + nmus = len(mus) + err = [0.] * nmus + self.assembleReducedSystem() + self.assembleReducedResidualBlocks(homogeneized) + assert homogeneized == self.homogeneized + nAs = len(self.As) + nbs = len(self.bs) + for j in range(nmus): + mu = mus[j] + uAppReduced = self.getAppReduced(mu) + prodbb = 0. + prodAb = 0. + prodAA = 0. + for i1 in range(nbs): + rhobi1 = self.thetabs(mu, i1) + for i2 in range(nbs): + rhobi2 = self.thetabs(mu, i2).conj() + prodbb += (rhobi1 * rhobi2 + * self.resbb[homogeneized][i1][i2]) + for i1 in range(nAs): + rhoAi1 = self.thetaAs(mu, i1) + for i2 in range(nbs): + rhobi2 = self.thetabs(mu, i2).conj() + prodAb += (rhoAi1 * rhobi2 + * self.resAb[homogeneized][i1][i2]) + for i1 in range(nAs): + rhoAi1 = self.thetaAs(mu, i1) + for i2 in range(nAs): + rhoAi2 = self.thetaAs(mu, i2).conj() + prodAA += (rhoAi1 * rhoAi2 + * self.resAA[homogeneized][i1][i2]) + err[j] = np.abs(((uAppReduced.T.conj().dot(prodAA.dot(uAppReduced)) + - 2. * np.real(prodAb.dot(uAppReduced))) / prodbb + + 1.)[0]) ** .5 + return err + + def setupApprox(self): + """Compute RB projection matrix.""" + if not self.checkComputedApprox(): + if self.verbosity >= 5: + verbosityDepth("INIT", "Setting up {}.". format(self.name())) + self.greedy() + self.S = len(self.mus) + if self.verbosity >= 7: + verbosityDepth("INIT", "Computing projection matrix.", + end = "") + self.projMat = self.samplingEngine.samples + if self.verbosity >= 7: + verbosityDepth("DEL", " Done.", inline = True) + self.lastApproxParameters = copy(self.approxParameters) + if hasattr(self, "lastSolvedApp"): del self.lastSolvedApp + if self.verbosity >= 5: + verbosityDepth("DEL", "Done setting up approximant.\n") + + def assembleReducedResidualBlocks(self, homogeneized:bool): + """Build affine blocks of RB linear system through projections.""" + self.initMassMatrixInverse() + self.assembleReducedSystem() + computeResbb = (self.resbb is None + or homogeneized not in self.resbb.keys()) + computeResAb = (self.resAb is None + or homogeneized not in self.resAb.keys() + or self.resAb[homogeneized][0][0].shape[0] != self.S) + computeResAA = (self.resAA is None + or homogeneized not in self.resAA.keys() + or self.resAA[homogeneized][0][0].shape[0] != self.S) + if computeResbb or computeResAb or computeResAA: + if self.verbosity >= 7: + verbosityDepth("INIT", "Projecting affine terms of residual.", + end = "") + nAs = len(self.As) + nbs = len(self.bs) + if computeResbb: + if self.resbb is None: self.resbb = {} + self.resbb[homogeneized] = [] + for i in range(nbs): + resbbi = [] + Mbi = self.massInv.solve(self.bs[i]) + for j in range(i): + Mbj = self.massInv.solve(self.bs[j]) + resbbi += [self.HFEngine.innerProduct(Mbi, Mbj)] + resbbi += [self.HFEngine.innerProduct(Mbi, Mbi)] + self.resbb[homogeneized] += [resbbi] + for i in range(nbs): + for j in range(i + 1, nbs): + self.resbb[homogeneized][i] += [ + self.resbb[homogeneized][j][i].conj()] + if computeResAb: + if self.resAb is None: self.resAb = {} + if homogeneized not in self.resAb.keys(): + self.resAb[homogeneized] = [] + for i in range(nAs): + resAbi = [] + MAi = self.massInv.solve(self.As[i].dot(self.projMat)) + for j in range(nbs): + Mbj = self.massInv.solve(self.bs[j]) + resAbi += [self.HFEngine.innerProduct(MAi, Mbj)] + self.resAb[homogeneized] += [resAbi] + else: + Sold = self.resAb[homogeneized][0][0].shape[0] + for i in range(nAs): + MAi = self.massInv.solve(self.As[i].dot( + self.projMat[:, Sold :])) + for j in range(nbs): + Mbj = self.massInv.solve(self.bs[j]) + self.resAb[homogeneized][i][j] = np.append( + self.resAb[homogeneized][i][j], + self.HFEngine.innerProduct(MAi, Mbj)) + if computeResAA: + if self.resAA is None: self.resAA = {} + if homogeneized not in self.resAA.keys(): + self.resAA[homogeneized] = [] + for i in range(nAs): + resAAi = [] + MAi = self.massInv.solve(self.As[i].dot(self.projMat)) + for j in range(i): + MAj = self.massInv.solve(self.As[j].dot( + self.projMat)) + resAAi += [self.HFEngine.innerProduct(MAi, MAj)] + resAAi += [self.HFEngine.innerProduct(MAi, MAi)] + self.resAA[homogeneized] += [resAAi] + for i in range(nAs): + for j in range(i + 1, nAs): + self.resAA[homogeneized][i] += [ + self.resAA[homogeneized][j][i].conj()] + else: + Sold = self.resAA[homogeneized][0][0].shape[0] + for i in range(nAs): + MAi = self.massInv.solve(self.As[i].dot(self.projMat)) + for j in range(i): + MAj = self.massInv.solve(self.As[j].dot( + self.projMat)) + resAAcross1 = self.HFEngine.innerProduct( + MAi[:, Sold :], MAj[:, : Sold]) + resAAcross2 = self.HFEngine.innerProduct( + MAi[:, : Sold], MAj[:, Sold :]) + resAAdiag = self.HFEngine.innerProduct( + MAi[:, Sold :], MAj[:, Sold :]) + self.resAA[homogeneized][i][j] = np.block( + [[self.resAA[homogeneized][i][j], resAAcross1], + [ resAAcross2, resAAdiag]]) + resAAcross = self.HFEngine.innerProduct(MAi[:, Sold :], + MAi[:, : Sold]) + resAAdiag = self.HFEngine.innerProduct(MAi[:, Sold :], + MAi[:, Sold :]) + self.resAA[homogeneized][i][i] = np.block( + [[self.resAA[homogeneized][i][i], resAAcross], + [ resAAcross.T.conj(), resAAdiag]]) + for i in range(nAs): + for j in range(i + 1, nAs): + self.resAA[homogeneized][i][j] = ( + self.resAA[homogeneized][j][i].conj()) + if self.verbosity >= 7: + verbosityDepth("DEL", ("Done setting up affine " + "decomposition of residual.")) diff --git a/rrompy/reduction_methods/lagrange_greedy/generic_approximant_lagrange_greedy.py b/rrompy/reduction_methods/lagrange_greedy/generic_approximant_lagrange_greedy.py new file mode 100644 index 0000000..7b4d4ae --- /dev/null +++ b/rrompy/reduction_methods/lagrange_greedy/generic_approximant_lagrange_greedy.py @@ -0,0 +1,452 @@ +# 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.reduction_methods.base.generic_approximant import ( + GenericApproximant) +from rrompy.reduction_methods.lagrange.generic_approximant_lagrange import ( + GenericApproximantLagrange) +from rrompy.utilities.base.types import DictAny, HFEng, Tuple, List +from rrompy.utilities.base import purgeDict, verbosityDepth +from rrompy.utilities.warning_manager import warn + +__all__ = ['GenericApproximantLagrangeGreedy'] + +class GenericApproximantLagrangeGreedy(GenericApproximantLagrange): + """ + ROM greedy Lagrange 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]]; + - 'greedyTol': uniform error tolerance for greedy algorithm; + defaults to 1e-2; + - 'maxIter': maximum number of greedy steps; defaults to 1e2; + - 'refinementRatio': ratio of training points to be exhausted + before training set refinement; defaults to 0.2; + - 'nTrainingPoints': number of training points; defaults to + maxIter / refinementRatio; + - 'nTestPoints': number of starting test points; defaults to 1; + - 'trainingSetGenerator': training sample points generator; + defaults to uniform sampler within muBounds; + - 'robustTol': tolerance for robust Pade' denominator management; + defaults to 0. + Defaults to empty dict. + + 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. + parameterList: Recognized keys of approximant parameters: + - 'POD': whether to compute POD of snapshots; + - 'muBounds': list of bounds for parameter values; + - '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; + - 'nTrainingPoints': number of training points; + - 'nTestPoints': number of starting test points; + - 'trainingSetGenerator': training sample points generator. + - 'robustTol': tolerance for robust Pade' denominator management. + 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. + 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. + nTrainingPoints: number of training points. + nTestPoints: number of starting test points. + trainingSetGenerator: training sample points generator. + robustTol: tolerance for robust Pade' denominator management. + samplingEngine: Sampling engine. + 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. + lastApproxParameters: List of parameters corresponding to last + computed approximant. + """ + + TOL_INSTABILITY = 1e-6 + + def __init__(self, HFEngine:HFEng, mu0 : complex = 0., + approxParameters : DictAny = {}, homogeneized : bool = False, + verbosity : int = 10): + self._preInit() + if not hasattr(self, "parameterList"): + self.parameterList = [] + self.parameterList += ["muBounds","greedyTol", "maxIter", + "refinementRatio", "nTrainingPoints", + "nTestPoints", "trainingSetGenerator"] + super(GenericApproximantLagrange, self).__init__( + HFEngine = HFEngine, mu0 = mu0, + approxParameters = approxParameters, + homogeneized = homogeneized, + verbosity = verbosity) + 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, + ["muBounds","greedyTol", "maxIter", + "refinementRatio", "nTrainingPoints", + "nTestPoints", + "trainingSetGenerator"], + True, True, baselevel = 1) + GenericApproximant.approxParameters.fset(self, approxParametersCopy) + keyList = list(approxParameters.keys()) + if "muBounds" in keyList: + self.muBounds = approxParameters["muBounds"] + elif hasattr(self, "muBounds"): + self.muBounds = self.muBounds + else: + self.muBounds = [[0.], [1.]] + if "greedyTol" in keyList: + self.greedyTol = approxParameters["greedyTol"] + elif hasattr(self, "greedyTol"): + self.greedyTol = self.greedyTol + else: + self.greedyTol = 1e-2 + if "maxIter" in keyList: + self.maxIter = approxParameters["maxIter"] + elif hasattr(self, "maxIter"): + self.maxIter = self.maxIter + else: + self.maxIter = 1e2 + if "refinementRatio" in keyList: + self.refinementRatio = approxParameters["refinementRatio"] + elif hasattr(self, "refinementRatio"): + self.refinementRatio = self.refinementRatio + else: + self.refinementRatio = 0.2 + if "nTrainingPoints" in keyList: + self.nTrainingPoints = approxParameters["nTrainingPoints"] + elif hasattr(self, "nTrainingPoints"): + self.nTrainingPoints = self.nTrainingPoints + else: + self.nTrainingPoints = np.int(np.ceil(self.maxIter + / self.refinementRatio)) + if "nTestPoints" in keyList: + self.nTestPoints = approxParameters["nTestPoints"] + elif hasattr(self, "nTestPoints"): + self.nTestPoints = self.nTestPoints + else: + self.nTestPoints = 1 + if "trainingSetGenerator" in keyList: + self.trainingSetGenerator = ( + approxParameters["trainingSetGenerator"]) + elif hasattr(self, "trainingSetGenerator"): + self.trainingSetGenerator = self.trainingSetGenerator + else: + from rrompy.utilities.parameter_sampling import QuadratureSampler + self.trainingSetGenerator = QuadratureSampler(self.muBounds, + "UNIFORM") + del QuadratureSampler + + @property + def S(self): + """Value of S.""" + return self._S + @S.setter + def S(self, S): + self._S = S + + @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 muBounds(self): + """Value of muBounds.""" + return self._muBounds + @muBounds.setter + def muBounds(self, muBounds): + if len(muBounds) != 2: + raise Exception("2 limits must be specified.") + try: + muBounds = muBounds.tolist() + except: + muBounds = list(muBounds) + for j in range(2): + try: + len(muBounds[j]) + except: + muBounds[j] = np.array([muBounds[j]]) + if len(muBounds[0]) != len(muBounds[1]): + raise Exception("The bounds must have the same length.") + self._muBounds = muBounds + + @property + def greedyTol(self): + """Value of greedyTol.""" + return self._greedyTol + @greedyTol.setter + def greedyTol(self, greedyTol): + if greedyTol < 0: + raise ArithmeticError("greedyTol must be non-negative.") + if hasattr(self, "greedyTol"): 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 ArithmeticError("maxIter must be positive.") + if hasattr(self, "maxIter"): 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 ArithmeticError(("refinementRatio must be between 0 " + "(excluded) and 1.")) + if hasattr(self, "refinementRatio"): + refinementRatioold = self.refinementRatio + else: refinementRatioold = -1 + self._refinementRatio = refinementRatio + self._approxParameters["refinementRatio"] = self.refinementRatio + if refinementRatioold != self.refinementRatio: + self.resetSamples() + + @property + def nTrainingPoints(self): + """Value of nTrainingPoints.""" + return self._nTrainingPoints + @nTrainingPoints.setter + def nTrainingPoints(self, nTrainingPoints): + if nTrainingPoints <= 1: + raise ArithmeticError("nTrainingPoints must be greater than 1.") + if not np.isclose(nTrainingPoints, np.int(nTrainingPoints)): + raise ArithmeticError("nTrainingPoints must be an integer.") + nTrainingPoints = np.int(nTrainingPoints) + if hasattr(self, "nTrainingPoints"): + nTrainingPointsold = self.nTrainingPoints + else: nTrainingPointsold = -1 + self._nTrainingPoints = nTrainingPoints + self._approxParameters["nTrainingPoints"] = self.nTrainingPoints + if nTrainingPointsold != self.nTrainingPoints: + self.resetSamples() + + @property + def nTestPoints(self): + """Value of nTestPoints.""" + return self._nTestPoints + @nTestPoints.setter + def nTestPoints(self, nTestPoints): + if nTestPoints <= 0: + raise ArithmeticError("nTestPoints must be positive.") + if not np.isclose(nTestPoints, np.int(nTestPoints)): + raise ArithmeticError("nTestPoints must be an integer.") + nTestPoints = np.int(nTestPoints) + if hasattr(self, "nTestPoints"): + nTestPointsold = self.nTestPoints + else: nTestPointsold = -1 + self._nTestPoints = nTestPoints + self._approxParameters["nTestPoints"] = self.nTestPoints + if nTestPointsold != self.nTestPoints: + self.resetSamples() + + @property + def trainingSetGenerator(self): + """Value of trainingSetGenerator.""" + return self._trainingSetGenerator + @trainingSetGenerator.setter + def trainingSetGenerator(self, trainingSetGenerator): + if 'generatePoints' not in dir(trainingSetGenerator): + raise Exception("trainingSetGenerator type not recognized.") + if hasattr(self, '_trainingSetGenerator'): + trainingSetGeneratorOld = self.trainingSetGenerator + self._trainingSetGenerator = trainingSetGenerator + self._approxParameters["trainingSetGenerator"] = ( + self.trainingSetGenerator) + if (not 'trainingSetGeneratorOld' in locals() + or trainingSetGeneratorOld != self.trainingSetGenerator): + self.resetSamples() + + def resetSamples(self): + """Reset samples.""" + super().resetSamples() + self._mus = [] + + def initMassMatrixInverse(self): + """Initialize mass matrix for error estimator.""" + if not hasattr(self, "massInv"): + from scipy.sparse import csc_matrix, linalg as sla + from rrompy.hfengines.base import ProblemEngineBase as PEB + L2normer = PEB() + L2normer.V = self.HFEngine.V + L2normer.buildEnergyNormForm() + Mass = csc_matrix(L2normer.energyNormMatrix, dtype = np.complex) + self.massInv = sla.spilu(Mass) + del csc_matrix, sla, PEB + + def errorEstimator(self, mus:List[np.complex], + homogeneized : bool = None) -> List[np.complex]: + """ + Standard residual-based error estimator. Unreliable for unstable + problems (inf-sup constant is missing). + """ + if homogeneized is None: + homogeneized = self.homogeneized + self.setupApprox() + nmus = len(mus) + err = [0.] * nmus + self.initMassMatrixInverse() + if self.HFEngine.nbs == 1: + RHS = self.getRHS(mus[0], homogeneized = homogeneized) + RHSNorm = self.HFEngine.norm(self.massInv.solve(RHS)) + for j in range(nmus): + res = self.getRes(mus[j], homogeneized = homogeneized) + err[j] = (self.HFEngine.norm(self.massInv.solve(res)) + / RHSNorm) + else: + for j in range(nmus): + res = self.getRes(mus[j], homogeneized = homogeneized) + RHS = self.getRHS(mus[j], homogeneized = homogeneized) + err[j] = (self.HFEngine.norm(self.massInv.solve(res)) + / self.HFEngine.norm(self.massInv.solve(RHS))) + return err + + def getMaxErrorEstimator(self) -> Tuple[float, int]: + """ + Compute maximum of (and index of maximum of) error estimator over + training set. + """ + errorEstTrain = self.errorEstimator(self.muTrain) + idxMaxEst = np.argmax(errorEstTrain) + maxEst = errorEstTrain[idxMaxEst] + return maxEst, idxMaxEst + + def greedyNextSample(self, muidx:int, plotEst : bool = False): + """Compute next greedy snapshot of solution map.""" + mu = self.muTrain[muidx] + if self.verbosity >= 2: + verbosityDepth("MAIN", ("Adding {}-th sample point at {} to test " + "set.").format( + self.samplingEngine.nsamples + 1, mu)) + self.mus = np.append(self.mus, mu) + idxs = np.arange(len(self.muTrain)) + mask = np.ones_like(idxs, dtype = bool) + mask[muidx] = False + idxs = idxs[mask] + self.muTrain = self.muTrain[idxs] + self.samplingEngine.nextSample(mu, + homogeneized = self.homogeneized) + errorEstTrain = self.errorEstimator(self.muTrain) + muidx = np.argmax(errorEstTrain) + maxErrorEst = errorEstTrain[muidx] + mu = self.muTrain[muidx] + if plotEst: + from matplotlib import pyplot as plt + plt.figure() + plt.semilogy(self.muTrain, errorEstTrain, 'k') + plt.semilogy(self.muTrain, + self.greedyTol * np.ones(len(self.muTrain)), + 'r--') + plt.semilogy(mu, maxErrorEst, 'xr') + plt.grid() + plt.show() + plt.close() + return errorEstTrain, muidx, maxErrorEst, mu + + def greedy(self, plotEst : bool = False): + """Compute greedy snapshots of solution map.""" + if self.samplingEngine.samples is None: + if self.verbosity >= 2: + verbosityDepth("INIT", "Starting computation of snapshots.") + self.resetSamples() + nTrain = self.nTrainingPoints + self.muTrain, _ = self.trainingSetGenerator.generatePoints(nTrain) + self.muTrain = np.array([x[0] for x in self.muTrain]) + maxErrorEst, muidx, mu = np.inf, 0, self.muTrain[0] + while (self.samplingEngine.nsamples < self.maxIter + and (maxErrorEst > self.greedyTol + or self.samplingEngine.nsamples < self.nTestPoints)): + if (1. - self.refinementRatio) * nTrain > len(self.muTrain): + muTrainExtra, _ = self.trainingSetGenerator.refine(nTrain) + self.muTrain = np.sort(np.append(self.muTrain, + muTrainExtra)) + nTrain += len(muTrainExtra) + if self.verbosity >= 5: + verbosityDepth("MAIN", ("Enriching training set by {} " + "elements.").format( + len(muTrainExtra))) + muTrainOld, maxErrorEstOld = self.muTrain, maxErrorEst + errorEstTrain, muidx, maxErrorEst, mu = self.greedyNextSample( + muidx, plotEst) + if self.verbosity >= 2: + verbosityDepth("MAIN", ("Uniform error estimate {:.4e}.")\ + .format(maxErrorEst)) + if maxErrorEstOld < maxErrorEst * self.TOL_INSTABILITY: + warn(("Instability in a posteriori estimator. Starting " + "preemptive greedy loop termination.")) + self.muTrain = muTrainOld + self.mus = self.mus[:-1] + self.samplingEngine.popSample() + self.setupApprox() + maxErrorEst = maxErrorEstOld + break + if self.verbosity >= 2: + verbosityDepth("DEL", "Done computing snapshots.") + + def checkComputedApprox(self) -> bool: + """ + Check if setup of new approximant is not needed. + + Returns: + True if new setup is not needed. False otherwise. + """ + return (hasattr(self, "_S") and self.S == len(self.mus) + and super().checkComputedApprox()) + diff --git a/rrompy/reduction_methods/taylor/approximant_taylor_pade.py b/rrompy/reduction_methods/taylor/approximant_taylor_pade.py index 7a278e7..aad876e 100644 --- a/rrompy/reduction_methods/taylor/approximant_taylor_pade.py +++ b/rrompy/reduction_methods/taylor/approximant_taylor_pade.py @@ -1,588 +1,578 @@ # 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 copy import numpy as np +from rrompy.reduction_methods.base import checkRobustTolerance from .generic_approximant_taylor import GenericApproximantTaylor from rrompy.sampling.base.pod_engine import PODEngine from rrompy.utilities.base.types import Np1D, Np2D, List, Tuple, DictAny from rrompy.utilities.base.types import HFEng from rrompy.utilities.base import purgeDict, verbosityDepth from rrompy.utilities.warning_manager import warn __all__ = ['ApproximantTaylorPade'] class ApproximantTaylorPade(GenericApproximantTaylor): """ ROM single-point fast Pade' approximant computation for parametric problems. Args: HFEngine: HF problem solver. mu0(optional): Default parameter. Defaults to 0. approxParameters(optional): Dictionary containing values for main parameters of approximant. Recognized keys are: - 'POD': whether to compute POD of snapshots; defaults to True; - 'rho': weight for computation of original Pade' approximant; defaults to np.inf, i.e. fast approximant; - 'M': degree of Pade' approximant numerator; defaults to 0; - 'N': degree of Pade' approximant denominator; defaults to 0; - 'E': total number of derivatives current approximant relies upon; defaults to Emax; - 'Emax': total number of derivatives of solution map to be computed; defaults to E; - 'robustTol': tolerance for robust Pade' denominator management; defaults to 0; - 'rescaleParam': whether to rescale parameter during denominator computation; defaults to True; - 'sampleType': label of sampling type; available values are: - 'ARNOLDI': orthogonalization of solution derivatives through Arnoldi algorithm; - 'KRYLOV': standard computation of solution derivatives. Defaults to 'KRYLOV'. Defaults to empty dict. Attributes: HFEngine: HF problem solver. mu0: Default parameter. 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; - 'rho': weight for computation of original Pade' approximant; - 'M': degree of Pade' approximant numerator; - 'N': degree of Pade' approximant denominator; - 'E': total number of derivatives current approximant relies upon; - 'Emax': total number of derivatives of solution map to be computed; - 'robustTol': tolerance for robust Pade' denominator management; - 'rescaleParam': whether to rescale parameter during denominator computation; - 'sampleType': label of sampling type. POD: Whether to compute QR factorization of derivatives. rho: Weight of approximant. M: Numerator degree of approximant. N: Denominator degree of approximant. E: Number of solution derivatives over which current approximant is based upon. Emax: Total number of solution derivatives to be computed. robustTol: Tolerance for robust Pade' denominator management. sampleType: Label of sampling type. initialHFData: HF problem initial data. uHF: High fidelity solution with wavenumber lastSolvedHF as numpy complex vector. lastSolvedHF: Wavenumber corresponding to last computed high fidelity solution. G: Square Numpy 2D vector of size (N+1) corresponding to Pade' denominator matrix (see paper). Q: Numpy 1D vector containing complex coefficients of approximant denominator. P: Numpy 2D vector whose columns are FE dofs of coefficients of approximant numerator. uApp: Last evaluated approximant as numpy complex vector. lastApproxParameters: List of parameters corresponding to last computed approximant. """ def __init__(self, HFEngine:HFEng, mu0 : complex = 0, - approxParameters : DictAny = {}, homogeneize : bool = False, + approxParameters : DictAny = {}, homogeneized : bool = False, verbosity : int = 10): self._preInit() if not hasattr(self, "parameterList"): self.parameterList = [] self.parameterList += ["M", "N", "robustTol", "rho", "rescaleParam"] super().__init__(HFEngine = HFEngine, mu0 = mu0, approxParameters = approxParameters, - homogeneize = homogeneize, + homogeneized = homogeneized, verbosity = verbosity) self._postInit() @property def approxParameters(self): """Value of approximant parameters.""" return self._approxParameters @approxParameters.setter def approxParameters(self, approxParams): approxParameters = purgeDict(approxParams, self.parameterList, dictname = self.name() + ".approxParameters", baselevel = 1) approxParametersCopy = purgeDict(approxParameters, ["M", "N", "robustTol", "rho"], True, True, baselevel = 1) keyList = list(approxParameters.keys()) if "rho" in keyList: self._rho = approxParameters["rho"] elif hasattr(self, "rho"): self._rho = self.rho else: self._rho = np.inf GenericApproximantTaylor.approxParameters.fset(self, approxParametersCopy) self.rho = self._rho if "robustTol" in keyList: self.robustTol = approxParameters["robustTol"] elif hasattr(self, "robustTol"): self.robustTol = self.robustTol else: self.robustTol = 0 if "rescaleParam" in keyList: self.rescaleParam = approxParameters["rescaleParam"] elif hasattr(self, "rescaleParam"): self.rescaleParam = self.rescaleParam else: self.rescaleParam = True self._ignoreParWarnings = True if "M" in keyList: self.M = approxParameters["M"] elif hasattr(self, "M"): self.M = self.M else: self.M = 0 del self._ignoreParWarnings if "N" in keyList: self.N = approxParameters["N"] elif hasattr(self, "N"): self.N = self.N else: self.N = 0 @property def rho(self): """Value of rho.""" return self._rho @rho.setter def rho(self, rho): self._rho = np.abs(rho) self._approxParameters["rho"] = self.rho @property def M(self): """Value of M. Its assignment may change Emax and E.""" return self._M @M.setter def M(self, M): if M < 0: raise ArithmeticError("M must be non-negative.") self._M = M self._approxParameters["M"] = self.M if not hasattr(self, "_ignoreParWarnings"): self.checkMNEEmax() @property def N(self): """Value of N. Its assignment may change Emax.""" return self._N @N.setter def N(self, N): if N < 0: raise ArithmeticError("N must be non-negative.") self._N = N self._approxParameters["N"] = self.N self.checkMNEEmax() def checkMNEEmax(self): """Check consistency of M, N, E, and Emax.""" M = self.M if hasattr(self, "_M") else 0 N = self.N if hasattr(self, "_N") else 0 E = self.E if hasattr(self, "_E") else M + N Emax = self.Emax if hasattr(self, "_Emax") else M + N if self.rho == np.inf: if Emax < max(N, M): warn(("Prescribed Emax is too small. Updating Emax to " "max(M, N).")) self.Emax = max(N, M) if E < max(N, M): warn("Prescribed E is too small. Updating E to max(M, N).") self.E = max(N, M) else: if Emax < N + M: warn("Prescribed Emax is too small. Updating Emax to M + N.") self.Emax = self.N + M if E < N + M: warn("Prescribed E is too small. Updating E to M + N.") self.E = self.N + M @property def robustTol(self): """Value of tolerance for robust Pade' denominator management.""" return self._robustTol @robustTol.setter def robustTol(self, robustTol): if robustTol < 0.: warn("Overriding prescribed negative robustness tolerance to 0.") robustTol = 0. self._robustTol = robustTol self._approxParameters["robustTol"] = self.robustTol @property def rescaleParam(self): """Value of parameter rescaling during denominator computation.""" return self._rescaleParam @rescaleParam.setter def rescaleParam(self, rescaleParam): if not isinstance(rescaleParam, (bool,)): warn("Prescribed rescaleParam not recognized. Overriding to True.") rescaleParam = True self._rescaleParam = rescaleParam self._approxParameters["rescaleParam"] = self.rescaleParam @property def E(self): """Value of E. Its assignment may change Emax.""" return self._E @E.setter def E(self, E): if E < 0: raise ArithmeticError("E must be non-negative.") self._E = E self.checkMNEEmax() self._approxParameters["E"] = self.E if hasattr(self, "Emax") and self.Emax < self.E: warn("Prescribed Emax is too small. Updating Emax to E.") self.Emax = self.E @property def Emax(self): """Value of Emax. Its assignment may reset computed derivatives.""" return self._Emax @Emax.setter def Emax(self, Emax): if Emax < 0: raise ArithmeticError("Emax must be non-negative.") if hasattr(self, "Emax"): EmaxOld = self.Emax else: EmaxOld = -1 if hasattr(self, "_N"): N = self.N else: N = 0 if hasattr(self, "_M"): M = self.M else: M = 0 if hasattr(self, "_E"): E = self.E else: E = 0 if self.rho == np.inf: if max(N, M, E) > Emax: warn(("Prescribed Emax is too small. Updating Emax to " "max(N, M, E).")) Emax = max(N, M, E) else: if max(N + M, E) > Emax: warn(("Prescribed Emax is too small. Updating Emax to " "max(N + M, E).")) Emax = max(N + M, E) self._Emax = Emax self._approxParameters["Emax"] = Emax if EmaxOld >= self.Emax and self.samplingEngine.samples is not None: self.samplingEngine.samples = self.samplingEngine.samples[:, : self.Emax + 1] if (self.sampleType == "ARNOLDI" and self.samplingEngine.HArnoldi is not None): self.samplingEngine.HArnoldi = self.samplingEngine.HArnoldi[ : self.Emax + 1, : self.Emax + 1] self.samplingEngine.RArnoldi = self.samplingEngine.RArnoldi[ : self.Emax + 1, : self.Emax + 1] def getTargetFunctionalOptimum(self): """Compute optimum of target LS functional.""" self.setupApprox() - if self._funcOptimumWarn: - warn("Target functional optimum numerically zero.") - return np.abs(self._funcOptimum) + if self._funcOptimumIllCond: + if self.POD: + jOpt = np.linalg.norm(self.P[:, 0]) + else: + jOpt = self.HFEngine.norm(self.samplingEngine.samples.dot( + self.P[:, 0])) + else: + jOpt = np.abs(self._funcOptimum) + return jOpt def setupApprox(self): """ Compute Pade' approximant. SVD-based robust eigenvalue management. """ if not self.checkComputedApprox(): if self.verbosity >= 5: verbosityDepth("INIT", "Setting up {}.". format(self.name())) self.computeDerivatives() if self.N > 0: if self.verbosity >= 7: verbosityDepth("INIT", ("Starting computation of " "denominator.")) while self.N > 0: if self.POD: ev, eV = self.findeveVGQR() self._funcOptimum = ev[0] else: ev, eV = self.findeveVGExplicit() self._funcOptimum = ev[0] ** .5 if ev[0] <= 0 or np.isclose(np.abs(ev[0] / ev[-1]), 0.): - self._funcOptimumWarn = True - else: - self._funcOptimumWarn = False - ts = self.robustTol * np.linalg.norm(ev) - Nnew = np.sum(np.abs(ev) >= ts) - diff = self.N - Nnew - if diff <= 0: break - Enew = self.E - diff - Mnew = min(self.M, Enew) - if Mnew == self.M: - strM = "" + self._funcOptimumIllCond = True else: - strM = ", M from {0} to {1},".format(self.M, Mnew) - print(("Smallest {0} eigenvalues below tolerance.\n" - "Reducing N from {1} to {2}{5} and E from {3} to " - "{4}.").format(diff + 1, self.N, Nnew, self.E, - Enew, strM)) - newParameters = {"N" : Nnew, "M" : Mnew, "E" : Enew} + self._funcOptimumIllCond = False + newParameters = checkRobustTolerance(ev, self.E, + self.robustTol) + if not newParameters: + break self.approxParameters = newParameters - if self.N <= 0: - eV = np.ones((1, 1)) + if self.N <= 0: + eV = np.ones((1, 1)) self.Q = np.poly1d(eV[:, 0]) if self.verbosity >= 7: verbosityDepth("DEL", "Done computing denominator.") else: + self._funcOptimumIllCond = True self.Q = np.poly1d([1]) - if self.sampleType == "KRYLOV": - self._funcOptimum = self.HFEngine.norm( - self.samplingEngine.samples[:, self.E]) - else: - self._funcOptimum = np.linalg.norm( - self.samplingEngine.RArnoldi[:, self.E]) - self._funcOptimumWarn = False if self.verbosity >= 7: verbosityDepth("INIT", "Starting computation of numerator.") self.P = np.zeros((self.Emax + 1, self.M + 1), dtype = np.complex) for i in range(self.Q.order): rng = np.arange(self.M + 1 - i) self.P[rng, - 1 - rng - i] = self.Q[i] if self.sampleType == "ARNOLDI": self.P = self.samplingEngine.RArnoldi.dot(self.P) if self.verbosity >= 7: verbosityDepth("DEL", "Done computing numerator.") self.lastApproxParameters = copy(self.approxParameters) if hasattr(self, "lastSolvedApp"): del self.lastSolvedApp if self.verbosity >= 5: verbosityDepth("DEL", "Done setting up approximant.\n") def rescaleParameter(self, R:Np2D, A : Np2D = None, exponent : float = 1.) -> Np2D: """ Prepare parameter rescaling. Args: R: Matrix whose columns need rescaling. A(optional): Matrix whose diagonal defines scaling factor. If None, previous value of scaleFactor is used. Defaults to None. exponent(optional): Exponent of scaling factor in matrix diagonal. Defaults to 1. Returns: Rescaled matrix. """ if A is not None: aDiag = np.diag(A) scaleCoeffs = np.polyfit(np.arange(A.shape[1]), np.log(aDiag), 1) self.scaleFactor = np.exp(- scaleCoeffs[0] / exponent) return np.multiply(R, np.power(self.scaleFactor,np.arange(R.shape[1]))) def buildG(self): """Assemble Pade' denominator matrix.""" self.computeDerivatives() if self.rho == np.inf: Nmin = self.E - self.N else: Nmin = self.M - self.N + 1 if self.sampleType == "KRYLOV": DerE = self.samplingEngine.samples[:, Nmin : self.E + 1] G = self.HFEngine.innerProduct(DerE, DerE) if self.rescaleParam: DerE = self.rescaleParameter(DerE, G, 2.) G = self.HFEngine.innerProduct(DerE, DerE) else: RArnE = self.samplingEngine.RArnoldi[: self.E + 1, Nmin : self.E + 1] if self.rescaleParam: RArnE = self.rescaleParameter(RArnE, RArnE[Nmin :, :]) G = RArnE.conj().T.dot(RArnE) if self.rho == np.inf: self.G = G else: Gbig = G self.G = np.zeros((self.N + 1, self.N + 1), dtype = np.complex) for k in range(self.E - self.M): self.G += self.rho ** (2 * k) * Gbig[k : k + self.N + 1, k : k + self.N + 1] def findeveVGExplicit(self) -> Tuple[Np1D, Np2D]: """ Compute explicitly eigenvalues and eigenvectors of Pade' denominator matrix. """ if self.verbosity >= 10: verbosityDepth("INIT", "Building gramian matrix.") self.buildG() if self.verbosity >= 10: verbosityDepth("DEL", "Done building gramian.") if self.verbosity >= 7: verbosityDepth("INIT", "Solving eigenvalue problem for gramian matrix.") ev, eV = np.linalg.eigh(self.G) if self.rescaleParam: eV = self.rescaleParameter(eV.T).T if self.verbosity >= 7: verbosityDepth("DEL", "Done solving eigenvalue problem.") return ev, eV def findeveVGQR(self) -> Tuple[Np1D, Np2D]: """ Compute eigenvalues and eigenvectors of Pade' denominator matrix through SVD of R factor. See ``Householder triangularization of a quasimatrix'', L.Trefethen, 2008 for QR algorithm. Returns: Eigenvalues in ascending order and corresponding eigenvector matrix. """ self.computeDerivatives() if self.rho == np.inf: Nmin = self.E - self.N else: Nmin = self.M - self.N + 1 if self.sampleType == "KRYLOV": A = copy(self.samplingEngine.samples[:, Nmin : self.E + 1]) self.PODEngine = PODEngine(self.HFEngine) if self.verbosity >= 10: verbosityDepth("INIT", "Orthogonalizing samples.", end = "") R = self.PODEngine.QRHouseholder(A, only_R = True) if self.verbosity >= 10: verbosityDepth("DEL", " Done.", inline = True) else: R = self.samplingEngine.RArnoldi[: self.E + 1, Nmin : self.E + 1] if self.rescaleParam: R = self.rescaleParameter(R, R[Nmin :, :]) if self.rho == np.inf: if self.verbosity >= 10: verbosityDepth("INIT", ("Solving svd for square root of " "gramian matrix."), end = "") _, s, V = np.linalg.svd(R, full_matrices = False) else: if self.verbosity >= 10: verbosityDepth("INIT", ("Building matrix stack for square " "root of gramian."), end = "") Rtower = np.zeros((R.shape[0] * (self.E - self.M), self.N + 1), dtype = np.complex) for k in range(self.E - self.M): Rtower[k * R.shape[0] : (k + 1) * R.shape[0], :] = ( self.rho ** k * R[:, self.M - self.N + 1 + k : self.M + 1 + k]) if self.verbosity >= 10: verbosityDepth("DEL", " Done.", inline = True) if self.verbosity >= 7: verbosityDepth("INIT", ("Solving svd for square root of " "gramian matrix."), end = "") _, s, V = np.linalg.svd(Rtower, full_matrices = False) eV = V.conj().T[:, ::-1] if self.rescaleParam: eV = self.rescaleParameter(eV.T).T if self.verbosity >= 7: verbosityDepth("DEL", " Done.", inline = True) return s[::-1], eV def radiusPade(self, mu:Np1D, mu0 : float = None) -> float: """ Compute translated radius to be plugged into Pade' approximant. Args: mu: Parameter(s) 1. mu0: Parameter(s) 2. If None, set to self.mu0. Returns: Translated radius to be plugged into Pade' approximant. """ if mu0 is None: mu0 = self.mu0 return self.HFEngine.rescaling(mu) - self.HFEngine.rescaling(mu0) def getPVal(self, mu:List[complex]): """ Evaluate Pade' numerator at arbitrary parameter. Args: mu: Target parameter. """ self.setupApprox() if self.verbosity >= 10: verbosityDepth("INIT", "Evaluating numerator at mu = {}.".format(mu), end = "") try: len(mu) except: mu = [mu] powerlist = np.vander(self.radiusPade(mu), self.M + 1).T p = self.P.dot(powerlist) if len(mu) == 1: p = p.flatten() if self.verbosity >= 10: verbosityDepth("DEL", " Done.", inline = True) return p def getQVal(self, mu:List[complex]): """ Evaluate Pade' denominator at arbitrary parameter. Args: mu: Target parameter. """ self.setupApprox() if self.verbosity >= 10: verbosityDepth("INIT", "Evaluating denominator at mu = {}.".format(mu), end = "") q = self.Q(self.radiusPade(mu)) if self.verbosity >= 10: verbosityDepth("DEL", " Done.", inline = True) return q def evalApproxReduced(self, mu:complex): """ Evaluate Pade' approximant at arbitrary parameter. Args: mu: Target parameter. """ self.setupApprox() if (not hasattr(self, "lastSolvedApp") or not np.isclose(self.lastSolvedApp, mu)): if self.verbosity >= 5: verbosityDepth("INIT", "Evaluating approximant at mu = {}.".format(mu)) self.uAppReduced = self.getPVal(mu) / self.getQVal(mu) self.lastSolvedApp = mu if self.verbosity >= 5: verbosityDepth("DEL", "Done evaluating approximant.") def evalApprox(self, mu:complex): """ Evaluate approximant at arbitrary parameter. Args: mu: Target parameter. """ self.evalApproxReduced(mu) self.uApp = self.samplingEngine.samples.dot(self.uAppReduced) def getPoles(self) -> Np1D: """ Obtain approximant poles. Returns: Numpy complex vector of poles. """ self.setupApprox() return self.HFEngine.rescalingInv(self.Q.r + self.HFEngine.rescaling(self.mu0)) diff --git a/rrompy/reduction_methods/taylor/approximant_taylor_rb.py b/rrompy/reduction_methods/taylor/approximant_taylor_rb.py index aa209f4..5dba3d2 100644 --- a/rrompy/reduction_methods/taylor/approximant_taylor_rb.py +++ b/rrompy/reduction_methods/taylor/approximant_taylor_rb.py @@ -1,313 +1,313 @@ # 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 copy import numpy as np import scipy as sp from .generic_approximant_taylor import GenericApproximantTaylor from rrompy.sampling.base.pod_engine import PODEngine from rrompy.utilities.base.types import Np1D, DictAny, HFEng from rrompy.utilities.base import purgeDict, verbosityDepth from rrompy.utilities.warning_manager import warn __all__ = ['ApproximantTaylorRB'] class ApproximantTaylorRB(GenericApproximantTaylor): """ ROM single-point fast RB approximant computation for parametric problems with polynomial dependence up to degree 2. 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; - 'R': rank for Galerkin projection; defaults to E + 1; - 'E': total number of derivatives current approximant relies upon; defaults to Emax; - 'Emax': total number of derivatives of solution map to be computed; defaults to E; - 'sampleType': label of sampling type; available values are: - 'ARNOLDI': orthogonalization of solution derivatives through Arnoldi algorithm; - 'KRYLOV': standard computation of solution derivatives. Defaults to 'KRYLOV'. Defaults to empty dict. Attributes: HFEngine: HF problem solver. mu0: Default parameter. 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; - 'R': rank for Galerkin projection; - 'E': total number of derivatives current approximant relies upon; - 'Emax': total number of derivatives of solution map to be computed; - 'sampleType': label of sampling type. POD: Whether to compute QR factorization of derivatives. R: Rank for Galerkin projection. E: Number of solution derivatives over which current approximant is based upon. Emax: Total number of solution derivatives to be computed. sampleType: Label of sampling type, i.e. 'KRYLOV'. 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. lastApproxParameters: List of parameters corresponding to last computed approximant. projMat: Numpy matrix representing projection onto RB space. projMat: Numpy matrix representing projection onto RB space. As: List of sparse matrices (in CSC format) representing coefficients of linear system matrix wrt mu. bs: List of numpy vectors representing coefficients of linear system RHS wrt mu. ARBs: List of sparse matrices (in CSC format) representing RB coefficients of linear system matrix wrt mu. bRBs: List of numpy vectors representing RB coefficients of linear system RHS wrt mu. """ def __init__(self, HFEngine:HFEng, mu0 : complex = 0, - approxParameters : DictAny = {}, homogeneize : bool = False, + approxParameters : DictAny = {}, homogeneized : bool = False, verbosity : int = 10): self._preInit() if not hasattr(self, "parameterList"): self.parameterList = [] self.parameterList += ["R"] super().__init__(HFEngine = HFEngine, mu0 = mu0, approxParameters = approxParameters, - homogeneize = homogeneize, + homogeneized = homogeneized, verbosity = verbosity) if self.verbosity >= 10: verbosityDepth("INIT", "Computing affine blocks of system.") self.As, self.thetaAs = self.HFEngine.affineBlocksA(self.mu0) self.bs, self.thetabs = self.HFEngine.affineBlocksb(self.mu0, - self.homogeneize) + self.homogeneized) if self.verbosity >= 10: verbosityDepth("DEL", "Done computing affine blocks.") self._postInit() def resetSamples(self): """Reset samples.""" super().resetSamples() self.projMat = None @property def approxParameters(self): """ Value of approximant parameters. Its assignment may change M, N and S. """ return self._approxParameters @approxParameters.setter def approxParameters(self, approxParams): approxParameters = purgeDict(approxParams, self.parameterList, dictname = self.name() + ".approxParameters", baselevel = 1) approxParametersCopy = purgeDict(approxParameters, ["R"], True, True, baselevel = 1) GenericApproximantTaylor.approxParameters.fset(self, approxParametersCopy) keyList = list(approxParameters.keys()) if "R" in keyList: self.R = approxParameters["R"] else: self.R = self.E + 1 @property def POD(self): """Value of POD.""" return self._POD @POD.setter def POD(self, POD): GenericApproximantTaylor.POD.fset(self, POD) if (hasattr(self, "sampleType") and self.sampleType == "ARNOLDI" and not self.POD): warn(("Arnoldi sampling implicitly forces POD-type derivative " "management.")) @property def sampleType(self): """Value of sampleType.""" return self._sampleType @sampleType.setter def sampleType(self, sampleType): GenericApproximantTaylor.sampleType.fset(self, sampleType) if (hasattr(self, "POD") and not self.POD and self.sampleType == "ARNOLDI"): warn(("Arnoldi sampling implicitly forces POD-type derivative " "management.")) @property def R(self): """Value of R. Its assignment may change S.""" return self._R @R.setter def R(self, R): if R < 0: raise ArithmeticError("R must be non-negative.") self._R = R self._approxParameters["R"] = self.R if hasattr(self, "E") and self.E + 1 < self.R: warn("Prescribed E is too small. Updating E to R - 1.") self.E = self.R - 1 def setupApprox(self): """Setup RB system.""" if not self.checkComputedApprox(): if self.verbosity >= 5: verbosityDepth("INIT", "Setting up {}.". format(self.name())) self.computeDerivatives() if self.verbosity >= 7: verbosityDepth("INIT", "Computing projection matrix.", end = "") if self.POD and not self.sampleType == "ARNOLDI": self.PODEngine = PODEngine(self.HFEngine) self.projMatQ, self.projMatR = self.PODEngine.QRHouseholder( self.samplingEngine.samples) if self.POD: if self.sampleType == "ARNOLDI": self.projMatR = self.samplingEngine.RArnoldi self.projMatQ = self.samplingEngine.samples U, _, _ = np.linalg.svd(self.projMatR[: self.E + 1, : self.E + 1]) self.projMat = self.projMatQ[:, : self.E + 1].dot(U[:, : self.R]) else: self.projMat = self.samplingEngine.samples[:, : self.R + 1] if self.verbosity >= 7: verbosityDepth("DEL", " Done.", inline = True) self.lastApproxParameters = copy(self.approxParameters) if hasattr(self, "lastSolvedApp"): del self.lastSolvedApp self.assembleReducedSystem() if self.verbosity >= 5: verbosityDepth("DEL", "Done setting up approximant.\n") def assembleReducedSystem(self): """Build affine blocks of RB linear system through projections.""" if not self.checkComputedApprox(): self.setupApprox() if self.verbosity >= 10: verbosityDepth("INIT", "Projecting affine terms of HF model.", end = "") projMatH = self.projMat.T.conj() self.ARBs = [None] * len(self.As) self.bRBs = [None] * len(self.bs) if self.verbosity >= 10: verbosityDepth("MAIN", ".", end = "", inline = True) for j in range(len(self.As)): self.ARBs[j] = projMatH.dot(self.As[j].dot(self.projMat)) if self.verbosity >= 10: verbosityDepth("MAIN", ".", end = "", inline = True) for j in range(len(self.bs)): self.bRBs[j] = projMatH.dot(self.bs[j]) if self.verbosity >= 10: verbosityDepth("DEL", "Done.", inline = True) def solveReducedSystem(self, mu:complex) -> Np1D: """ Solve RB linear system. Args: mu: Target parameter. Returns: Solution of RB linear system. """ self.setupApprox() if self.verbosity >= 10: verbosityDepth("INIT", "Assembling reduced model for mu = {}.".format(mu), end = "") ARBmu = self.thetaAs(mu, 0) * self.ARBs[0][:self.R,:self.R] bRBmu = self.thetabs(mu, 0) * self.bRBs[0][:self.R] if self.verbosity >= 10: verbosityDepth("MAIN", ".", end = "", inline = True) for j in range(1, len(self.ARBs)): ARBmu += self.thetaAs(mu, j) * self.ARBs[j][:self.R, :self.R] if self.verbosity >= 10: verbosityDepth("MAIN", ".", end = "", inline = True) for j in range(1, len(self.bRBs)): bRBmu += self.thetabs(mu, j) * self.bRBs[j][:self.R] if self.verbosity >= 10: verbosityDepth("DEL", "Done.", inline = True) if self.verbosity >= 5: verbosityDepth("INIT", "Solving reduced model for mu = {}.".format(mu), end = "") uRB = np.linalg.solve(ARBmu, bRBmu) if self.verbosity >= 5: verbosityDepth("DEL", " Done.", inline = True) return uRB def evalApproxReduced(self, mu:complex): """ Evaluate RB approximant at arbitrary wavenumber. Args: mu: Target parameter. """ self.setupApprox() if (not hasattr(self, "lastSolvedApp") or not np.isclose(self.lastSolvedApp, mu)): if self.verbosity >= 5: verbosityDepth("INIT", "Computing RB solution at mu = {}.".format(mu)) self.uAppReduced = self.solveReducedSystem(mu) self.lastSolvedApp = mu if self.verbosity >= 5: verbosityDepth("DEL", "Done computing RB solution.") def evalApprox(self, mu:complex): """ Evaluate approximant at arbitrary parameter. Args: mu: Target parameter. """ self.evalApproxReduced(mu) self.uApp = self.projMat[:, : self.R].dot(self.uAppReduced) def getPoles(self) -> Np1D: """ Obtain approximant poles. Returns: Numpy complex vector of poles. """ warn(("Impossible to compute poles in general affine parameter " "dependence. Results subject to interpretation/rescaling, or " "possibly completely wrong.")) self.setupApprox() if len(self.ARBs) < 2: return A = np.eye(self.ARBs[0].shape[0] * (len(self.ARBs) - 1), dtype = np.complex) B = np.zeros_like(A) A[: self.ARBs[0].shape[0], : self.ARBs[0].shape[0]] = - self.ARBs[0] for j in range(len(self.ARBs) - 1): Aj = self.ARBs[j + 1] B[: Aj.shape[0], j * Aj.shape[0] : (j + 1) * Aj.shape[0]] = Aj II = np.arange(self.ARBs[0].shape[0], self.ARBs[0].shape[0] * (len(self.ARBs) - 1)) B[II, II - self.ARBs[0].shape[0]] = 1. return self.HFEngine.rescalingInv(sp.linalg.eigvals(A, B) + self.HFEngine.rescaling(self.mu0)) diff --git a/rrompy/reduction_methods/taylor/generic_approximant_taylor.py b/rrompy/reduction_methods/taylor/generic_approximant_taylor.py index 6745854..25af630 100644 --- a/rrompy/reduction_methods/taylor/generic_approximant_taylor.py +++ b/rrompy/reduction_methods/taylor/generic_approximant_taylor.py @@ -1,243 +1,243 @@ # 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.reduction_methods.base.generic_approximant import ( GenericApproximant) from rrompy.utilities.base.types import DictAny, HFEng from rrompy.utilities.base import purgeDict, verbosityDepth from rrompy.utilities.warning_manager import warn __all__ = ['GenericApproximantTaylor'] class GenericApproximantTaylor(GenericApproximant): """ ROM single-point approximant 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; - 'E': total number of derivatives current approximant relies upon; defaults to Emax; - 'Emax': total number of derivatives of solution map to be computed; defaults to E; - 'sampleType': label of sampling type; available values are: - 'ARNOLDI': orthogonalization of solution derivatives through Arnoldi algorithm; - 'KRYLOV': standard computation of solution derivatives. Defaults to 'KRYLOV'. Defaults to empty dict. Attributes: HFEngine: HF problem solver. mu0: Default parameter. 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; - 'E': total number of derivatives current approximant relies upon; - 'Emax': total number of derivatives of solution map to be computed; - 'sampleType': label of sampling type. POD: Whether to compute QR factorization of derivatives. E: Number of solution derivatives over which current approximant is based upon. Emax: Total number of solution derivatives to be computed. sampleType: Label of sampling type. initialHFData: HF problem initial data. samplingEngine: Sampling engine. 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. lastApproxParameters: List of parameters corresponding to last computed approximant. """ def __init__(self, HFEngine:HFEng, mu0 : complex = 0, - approxParameters : DictAny = {}, homogeneize : bool = False, + approxParameters : DictAny = {}, homogeneized : bool = False, verbosity : int = 10): self._preInit() if not hasattr(self, "parameterList"): self.parameterList = [] self.parameterList += ["E", "Emax", "sampleType"] super().__init__(HFEngine = HFEngine, mu0 = mu0, approxParameters = approxParameters, - homogeneize = homogeneize, + homogeneized = homogeneized, verbosity = verbosity) self._postInit() def setupSampling(self): """Setup sampling engine.""" if not hasattr(self, "sampleType"): return if self.sampleType == "ARNOLDI": from rrompy.sampling.scipy.sampling_engine_arnoldi import ( SamplingEngineArnoldi) super().setupSampling(SamplingEngineArnoldi) elif self.sampleType == "KRYLOV": from rrompy.sampling.scipy.sampling_engine_krylov import ( SamplingEngineKrylov) super().setupSampling(SamplingEngineKrylov) else: raise Exception("Sample type not recognized.") @property def approxParameters(self): """ Value of approximant parameters. Its assignment may change E and Emax. """ return self._approxParameters @approxParameters.setter def approxParameters(self, approxParams): approxParameters = purgeDict(approxParams, self.parameterList, dictname = self.name() + ".approxParameters", baselevel = 1) approxParametersCopy = purgeDict(approxParameters, ["E", "Emax", "sampleType"], True, True, baselevel = 1) GenericApproximant.approxParameters.fset(self, approxParametersCopy) keyList = list(approxParameters.keys()) if "E" in keyList: self._E = approxParameters["E"] self._approxParameters["E"] = self.E if "Emax" in keyList: self.Emax = approxParameters["Emax"] else: if not hasattr(self, "Emax"): self.Emax = self.E else: self.Emax = self.Emax else: if "Emax" in keyList: self._E = approxParameters["Emax"] self._approxParameters["E"] = self.E self.Emax = self.E else: if not (hasattr(self, "Emax") and hasattr(self, "E")): raise Exception("At least one of E and Emax must be set.") if "sampleType" in keyList: self.sampleType = approxParameters["sampleType"] elif hasattr(self, "sampleType"): self.sampleType = self.sampleType else: self.sampleType = "KRYLOV" @property def E(self): """Value of E. Its assignment may change Emax.""" return self._E @E.setter def E(self, E): if E < 0: raise ArithmeticError("E must be non-negative.") self._E = E self._approxParameters["E"] = self.E if hasattr(self, "Emax") and self.Emax < self.E: warn("Prescribed E is too large. Updating Emax to E.") self.Emax = self.E @property def Emax(self): """Value of Emax. Its assignment may reset computed derivatives.""" return self._Emax @Emax.setter def Emax(self, Emax): if Emax < 0: raise ArithmeticError("Emax must be non-negative.") if hasattr(self, "Emax"): EmaxOld = self.Emax else: EmaxOld = -1 self._Emax = Emax if hasattr(self, "E") and self.Emax < self.E: warn("Prescribed Emax is too small. Updating Emax to E.") self.Emax = self.E else: self._approxParameters["Emax"] = self.Emax if (EmaxOld >= self.Emax and self.samplingEngine.samples is not None): self.samplingEngine.samples = self.samplingEngine.samples[:, : self.Emax + 1] if (self.sampleType == "ARNOLDI" and self.samplingEngine.HArnoldi is not None): self.samplingEngine.HArnoldi= self.samplingEngine.HArnoldi[ : self.Emax + 1, : self.Emax + 1] self.samplingEngine.RArnoldi= self.samplingEngine.RArnoldi[ : self.Emax + 1, : self.Emax + 1] else: self.resetSamples() @property def sampleType(self): """Value of sampleType.""" return self._sampleType @sampleType.setter def sampleType(self, sampleType): if hasattr(self, "sampleType"): sampleTypeOld = self.sampleType else: sampleTypeOld = -1 try: sampleType = sampleType.upper().strip().replace(" ","") if sampleType not in ["ARNOLDI", "KRYLOV"]: raise Exception("Sample type not recognized.") self._sampleType = sampleType except: warn(("Prescribed sampleType not recognized. Overriding to " "'KRYLOV'.")) self._sampleType = "KRYLOV" self._approxParameters["sampleType"] = self.sampleType if sampleTypeOld != self.sampleType: self.resetSamples() def computeDerivatives(self): """Compute derivatives of solution map starting from order 0.""" if self.samplingEngine.samples is None: if self.verbosity >= 5: verbosityDepth("INIT", "Starting computation of derivatives.") self.samplingEngine.iterSample(self.mu0, self.Emax + 1, - homogeneize = self.homogeneize) + homogeneized = self.homogeneized) if self.verbosity >= 5: verbosityDepth("DEL", "Done computing derivatives.") def checkComputedApprox(self) -> bool: """ Check if setup of new approximant is not needed. Returns: True if new setup is not needed. False otherwise. """ return (self.samplingEngine.samples is not None and super().checkComputedApprox()) def normApp(self, mu:complex, homogeneized : bool = False) -> float: """ Compute norm of approximant at arbitrary parameter. Args: mu: Target parameter. homogeneized(optional): Whether to remove Dirichlet BC. Defaults to False. Returns: Target norm of approximant. """ - if self.sampleType != "ARNOLDI" or self.homogeneize != homogeneized: + if self.sampleType != "ARNOLDI" or self.homogeneized != homogeneized: return super().normApp(mu, homogeneized) return np.linalg.norm(self.getAppReduced(mu, homogeneized)) diff --git a/rrompy/sampling/base/sampling_engine_base.py b/rrompy/sampling/base/sampling_engine_base.py index 9a537f5..d245681 100644 --- a/rrompy/sampling/base/sampling_engine_base.py +++ b/rrompy/sampling/base/sampling_engine_base.py @@ -1,105 +1,122 @@ # Copyright (C) 2018 by the RROMPy authors # # This file is part of RROMPy. # # RROMPy is free software: you can redistribute it and/or modify # it under the terms of the GNU Lesser General Public License as published by # the Free Software Foundation, either version 3 of the License, or # (at your option) any later version. # # RROMPy is distributed in the hope that it will be useful, # but WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # GNU Lesser General Public License for more details. # # You should have received a copy of the GNU Lesser General Public License # along with RROMPy. If not, see . # +import numpy as np from rrompy.utilities.base.types import Np1D, HFEng, strLst from rrompy.utilities.base import verbosityDepth +from rrompy.utilities.warning_manager import warn __all__ = ['SamplingEngineBase'] class SamplingEngineBase: """HERE""" nameBase = 0 def __init__(self, HFEngine:HFEng, verbosity : int = 10): self.verbosity = verbosity if self.verbosity >= 10: verbosityDepth("INIT", "Initializing sampling engine of type {}.".format( self.name()), end = "") self.HFEngine = HFEngine if self.verbosity >= 10: verbosityDepth("DEL", " Done.", inline = True) 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 resetHistory(self): self.samples = None self.nsamples = 0 + def popSample(self): + if hasattr(self, "nsamples") and self.nsamples > 1: + if self.samples.shape[1] > self.nsamples: + warn(("More than 'nsamples' memory allocated for samples. " + "Popping empty sample column.")) + self.nsamples += 1 + self.samples = self.samples[:, : -1] + self.nsamples -= 1 + else: + self.resetHistory() + + def preallocateSamples(self, u:Np1D, n:int): + self.samples = np.empty((u.size, n), dtype = u.dtype) + self.samples[:, 0] = u + @property def HFEngine(self): """Value of HFEngine. Its assignment resets history.""" return self._HFEngine @HFEngine.setter def HFEngine(self, HFEngine): self._HFEngine = HFEngine self.resetHistory() def solveLS(self, mu:complex, RHS : Np1D = None, homogeneized : bool = False) -> Np1D: """ Solve linear system. Args: mu: Parameter value. Returns: Solution of system. """ if self.verbosity >= 5: verbosityDepth("INIT", "Solving HF model for mu = {}.".format(mu)) u = self.HFEngine.solve(mu, RHS, homogeneized) if self.verbosity >= 5: verbosityDepth("DEL", "Done solving HF model.") return u def plotSamples(self, name : str = "u", save : str = None, what : strLst = 'all', saveFormat : str = "eps", saveDPI : int = 100, **figspecs): """ Do some nice plots of the samples. Args: name(optional): Name to be shown as title of the plots. Defaults to 'u'. save(optional): Where to save plot(s). Defaults to None, i.e. no saving. what(optional): Which plots to do. If list, can contain 'ABS', 'PHASE', 'REAL', 'IMAG'. If str, same plus wildcard 'ALL'. Defaults to 'ALL'. saveFormat(optional): Format for saved plot(s). Defaults to "eps". saveDPI(optional): DPI for saved plot(s). Defaults to 100. figspecs(optional key args): Optional arguments for matplotlib figure creation. """ for j in range(self.nsamples): self.HFEngine.plot(self.samples[:, j], name = "{}_{}".format(name, j + self.nameBase), save = save, what = what, saveFormat = saveFormat, saveDPI = saveDPI, **figspecs) diff --git a/rrompy/sampling/scipy/sampling_engine_arnoldi.py b/rrompy/sampling/scipy/sampling_engine_arnoldi.py index d50de53..38f8133 100644 --- a/rrompy/sampling/scipy/sampling_engine_arnoldi.py +++ b/rrompy/sampling/scipy/sampling_engine_arnoldi.py @@ -1,129 +1,140 @@ # 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 copy import numpy as np from rrompy.sampling.base.pod_engine import PODEngine from .sampling_engine_krylov import SamplingEngineKrylov from rrompy.utilities.base.types import Np1D from rrompy.utilities.base import verbosityDepth __all__ = ['SamplingEngineArnoldi'] class SamplingEngineArnoldi(SamplingEngineKrylov): """HERE""" def resetHistory(self): super().resetHistory() self.HArnoldi = None self.RArnoldi = None self.RHSs = None self.samplesAug = None + def popSample(self): + if hasattr(self, "nsamples") and self.nsamples > 1: + self.HArnoldi = self.HArnoldi[: -1, : -1] + self.RArnoldi = self.RArnoldi[: -1, : -1] + if self.nsamples > 2: + self.RHSs = self.RHSs[:, : -1] + else: + self.RHSs = None + self.samplesAug = self.RHSs[self.HFEngine.V.dim() :, : -1] + super().popSample() + @property def HFEngine(self): """Value of HFEngine. Its assignment resets history.""" return self._HFEngine @HFEngine.setter def HFEngine(self, HFEngine): self._HFEngine = HFEngine self.resetHistory() self.PODEngine = PODEngine(self._HFEngine) def preprocesssamples(self): ns = self.nsamples if ns <= 0: return return self.samplesAug[:, ns - 1].reshape((-1,self.HFEngine.V.dim())).T def preprocessb(self, mu:complex, overwrite : bool = False, - homogeneize : bool = False): + homogeneized : bool = False): ns = self.nsamples - r = super().preprocessb(mu, overwrite, homogeneize) + r = super().preprocessb(mu, overwrite, homogeneized) if ns == 0: return r elif ns == 1: r = r / self.RArnoldi[0, 0] else: r = ((r - self.RHSs[:, :ns-1].dot(self.RArnoldi[:ns-1, ns-1])) / self.RArnoldi[ns-1, ns-1]) if overwrite: self.RHSs[:, ns - 1] = r else: if ns == 1: self.RHSs = r.reshape((- 1, 1)) else: self.RHSs = np.hstack((self.RHSs, r[:, None])) return r def postprocessu(self, u:Np1D, overwrite : bool = False): if self.verbosity >= 10: verbosityDepth("INIT", "Starting orthogonalization.") ns = self.nsamples nsAug = (ns + 1) * self.HFEngine.V.dim() if ns == 0: u, h, _ = self.PODEngine.GS(u, np.empty((0, 0))) r = h[0] uAug = copy(u) else: uAug = np.concatenate((self.samplesAug[self.HFEngine.V.dim() - nsAug :, ns - 1], u), axis = None) u, h, uAug = self.PODEngine.GS(u, self.samples[:, : ns], ns, uAug, self.samplesAug[- nsAug :, : ns]) if overwrite: self.HArnoldi[: ns + 1, ns] = h if ns > 0: r = self.HArnoldi[: ns + 1, 1 : ns + 1].dot( self.RArnoldi[: ns, ns - 1]) self.RArnoldi[: ns + 1, ns] = r self.samplesAug[- nsAug :, ns] = uAug else: if ns == 0: self.HArnoldi = h.reshape((1, 1)) self.RArnoldi = r.reshape((1, 1)) self.samplesAug = uAug.reshape((-1, 1)) else: self.HArnoldi=np.block([[ self.HArnoldi, h[:-1, None]], [np.zeros((1, ns)), h[-1]]]) if ns > 0: r = self.HArnoldi[: ns + 1, 1 : ns + 1].dot( self.RArnoldi[: ns, ns - 1]) self.RArnoldi=np.block([[ self.RArnoldi, r[:-1, None]], [np.zeros((1, ns)), r[-1]]]) self.samplesAug=np.vstack((np.zeros((self.HFEngine.V.dim(), ns)), self.samplesAug)) self.samplesAug = np.hstack((self.samplesAug, uAug[:, None])) if self.verbosity >= 10: verbosityDepth("DEL", "Done orthogonalizing.") return u def preallocateSamples(self, u:Np1D, n:int): super().preallocateSamples(u, n) h = self.HArnoldi r = self.RArnoldi saug = self.samplesAug self.HArnoldi = np.zeros((n, n), dtype = u.dtype) self.HArnoldi[0, 0] = h[0, 0] self.RArnoldi = np.zeros((n, n), dtype = u.dtype) self.RArnoldi[0, 0] = r[0, 0] self.RHSs = np.empty((u.size, n - 1), dtype = u.dtype) self.samplesAug = np.zeros((self.HFEngine.V.dim() * (n + 1), n), dtype = u.dtype) self.samplesAug[- self.HFEngine.V.dim() :, 0] = saug[:, 0] diff --git a/rrompy/sampling/scipy/sampling_engine_krylov.py b/rrompy/sampling/scipy/sampling_engine_krylov.py index 253d543..539eddf 100644 --- a/rrompy/sampling/scipy/sampling_engine_krylov.py +++ b/rrompy/sampling/scipy/sampling_engine_krylov.py @@ -1,87 +1,83 @@ # Copyright (C) 2018 by the RROMPy authors # # This file is part of RROMPy. # # RROMPy is free software: you can redistribute it and/or modify # it under the terms of the GNU Lesser General Public License as published by # the Free Software Foundation, either version 3 of the License, or # (at your option) any later version. # # RROMPy is distributed in the hope that it will be useful, # but WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # GNU Lesser General Public License for more details. # # You should have received a copy of the GNU Lesser General Public License # along with RROMPy. If not, see . # import numpy as np from rrompy.sampling.base.sampling_engine_base import SamplingEngineBase from rrompy.utilities.base.types import Np1D, Np2D from rrompy.utilities.base import verbosityDepth __all__ = ['SamplingEngineKrylov'] class SamplingEngineKrylov(SamplingEngineBase): """HERE""" def preprocesssamples(self): if self.samples is None: return return self.samples[:, : self.nsamples] def preprocessb(self, mu:complex, overwrite : bool = False, - homogeneize : bool = False): - return self.HFEngine.b(mu, self.nsamples, homogeneized = homogeneize) + homogeneized : bool = False): + return self.HFEngine.b(mu, self.nsamples, homogeneized = homogeneized) def postprocessu(self, u:Np1D, overwrite : bool = False): return u - def preallocateSamples(self, u:Np1D, n:int): - self.samples = np.empty((u.size, n), dtype = u.dtype) - self.samples[:, 0] = u - def nextSample(self, mu:complex, overwrite : bool = False, - homogeneize : bool = False) -> Np1D: + homogeneized : bool = False) -> Np1D: ns = self.nsamples if self.verbosity >= 10: verbosityDepth("INIT", ("Setting up computation of {}-th Taylor " "coefficient.").format(ns)) samplesOld = self.preprocesssamples() RHS = self.preprocessb(mu, overwrite = overwrite, - homogeneize = homogeneize) + homogeneized = homogeneized) for i in range(1, ns + 1): RHS -= self.HFEngine.A(mu, i).dot(samplesOld[:, - i]) if self.verbosity >= 10: verbosityDepth("DEL", "Done setting up for Taylor coefficient.") u = self.postprocessu(self.solveLS(mu, RHS = RHS, - homogeneized = homogeneize), + homogeneized = homogeneized), overwrite = overwrite) if overwrite: self.samples[:, ns] = u else: if ns == 0: self.samples = u[:, None] else: self.samples = np.hstack((self.samples, u[:, None])) self.nsamples += 1 return u def iterSample(self, mu:complex, n:int, - homogeneize : bool = False) -> Np2D: + homogeneized : bool = False) -> Np2D: if self.verbosity >= 5: verbosityDepth("INIT", "Starting sampling iterations at mu = {}."\ .format(mu)) if n <= 0: raise Exception(("Number of Krylov iterations must be positive.")) self.resetHistory() - u = self.nextSample(mu, homogeneize = homogeneize) + u = self.nextSample(mu, homogeneized = homogeneized) if n > 1: self.preallocateSamples(u, n) for _ in range(1, n): self.nextSample(mu, overwrite = True, - homogeneize = homogeneize) + homogeneized = homogeneized) if self.verbosity >= 5: verbosityDepth("DEL", "Finished sampling iterations.") return self.samples diff --git a/rrompy/sampling/scipy/sampling_engine_lagrange.py b/rrompy/sampling/scipy/sampling_engine_lagrange.py index 8ff48f2..67d3f51 100644 --- a/rrompy/sampling/scipy/sampling_engine_lagrange.py +++ b/rrompy/sampling/scipy/sampling_engine_lagrange.py @@ -1,69 +1,65 @@ # 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.sampling.base.sampling_engine_base import SamplingEngineBase from rrompy.utilities.base.types import Np1D, Np2D from rrompy.utilities.base import verbosityDepth __all__ = ['SamplingEngineLagrange'] class SamplingEngineLagrange(SamplingEngineBase): """HERE""" nameBase = 1 def postprocessu(self, u:Np1D, overwrite : bool = False): return u - def preallocateSamples(self, u:Np1D, n:int): - self.samples = np.empty((u.size, n), dtype = u.dtype) - self.samples[:, 0] = u - def nextSample(self, mu:complex, overwrite : bool = False, - homogeneize : bool = False) -> Np1D: + homogeneized : bool = False) -> Np1D: ns = self.nsamples - u = self.postprocessu(self.solveLS(mu, homogeneized = homogeneize), + u = self.postprocessu(self.solveLS(mu, homogeneized = homogeneized), overwrite = overwrite) if overwrite: self.samples[:, ns] = u else: if ns == 0: self.samples = u[:, None] else: self.samples = np.hstack((self.samples, u[:, None])) self.nsamples += 1 return u - def iterSample(self, mus:Np1D, homogeneize : bool = False) -> Np2D: + def iterSample(self, mus:Np1D, homogeneized : bool = False) -> Np2D: if self.verbosity >= 5: verbosityDepth("INIT", "Starting sampling iterations.") n = mus.size if n <= 0: raise Exception(("Number of samples must be positive.")) self.resetHistory() - u = self.nextSample(mus[0], homogeneize = homogeneize) + u = self.nextSample(mus[0], homogeneized = homogeneized) if n > 1: self.preallocateSamples(u, n) for j in range(1, n): self.nextSample(mus[j], overwrite = True, - homogeneize = homogeneize) + homogeneized = homogeneized) if self.verbosity >= 5: verbosityDepth("DEL", "Finished sampling iterations.") return self.samples diff --git a/rrompy/sampling/scipy/sampling_engine_lagrange_pod.py b/rrompy/sampling/scipy/sampling_engine_lagrange_pod.py index b239295..273f17b 100644 --- a/rrompy/sampling/scipy/sampling_engine_lagrange_pod.py +++ b/rrompy/sampling/scipy/sampling_engine_lagrange_pod.py @@ -1,70 +1,75 @@ # 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.sampling.base.pod_engine import PODEngine from .sampling_engine_lagrange import SamplingEngineLagrange from rrompy.utilities.base.types import Np1D from rrompy.utilities.base import verbosityDepth __all__ = ['SamplingEngineLagrangePOD'] class SamplingEngineLagrangePOD(SamplingEngineLagrange): """HERE""" def resetHistory(self): super().resetHistory() self.RPOD = None + def popSample(self): + if hasattr(self, "nsamples") and self.nsamples > 1: + self.RPOD = self.RPOD[: -1, : -1] + super().popSample() + @property def HFEngine(self): """Value of HFEngine. Its assignment resets history.""" return self._HFEngine @HFEngine.setter def HFEngine(self, HFEngine): self._HFEngine = HFEngine self.resetHistory() self.PODEngine = PODEngine(self._HFEngine) def postprocessu(self, u:Np1D, overwrite : bool = False): if self.verbosity >= 10: verbosityDepth("INIT", "Starting orthogonalization.") ns = self.nsamples if ns == 0: u, r, _ = self.PODEngine.GS(u, np.empty((0, 0))) r = r[0] else: u, r, _ = self.PODEngine.GS(u, self.samples[:, : ns], ns) if overwrite: self.RPOD[: ns + 1, ns] = r else: if ns == 0: self.RPOD = r.reshape((1, 1)) else: self.RPOD=np.block([[ self.RPOD, r[:-1, None]], [np.zeros((1, ns)), r[-1]]]) if self.verbosity >= 10: verbosityDepth("DEL", "Done orthogonalizing.") return u def preallocateSamples(self, u:Np1D, n:int): super().preallocateSamples(u, n) r = self.RPOD self.RPOD = np.zeros((n, n), dtype = u.dtype) self.RPOD[0, 0] = r[0, 0] diff --git a/rrompy/utilities/parameter_sampling/fft_sampler.py b/rrompy/utilities/parameter_sampling/fft_sampler.py index dfdfa58..57c6a88 100644 --- a/rrompy/utilities/parameter_sampling/fft_sampler.py +++ b/rrompy/utilities/parameter_sampling/fft_sampler.py @@ -1,61 +1,97 @@ # Copyright (C) 2018 by the RROMPy authors # # This file is part of RROMPy. # # RROMPy is free software: you can redistribute it and/or modify # it under the terms of the GNU Lesser General Public License as published by # the Free Software Foundation, either version 3 of the License, or # (at your option) any later version. # # RROMPy is distributed in the hope that it will be useful, # but WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # GNU Lesser General Public License for more details. # # You should have received a copy of the GNU Lesser General Public License # along with RROMPy. If not, see . # import numpy as np from rrompy.utilities.parameter_sampling.generic_sampler import GenericSampler from rrompy.utilities.base.types import Np1D, List, Tuple __all__ = ['FFTSampler'] class FFTSampler(GenericSampler): """Generator of FFT-type sample points on scaled roots of unity.""" def generatePoints(self, n:Np1D) -> Tuple[List[Np1D], Np1D]: """Array of sample points and array of weights.""" super().generatePoints(n) d = len(self.lims[0]) try: len(n) except: n = np.array([n]) if len(n) != d: raise Exception(("Numbers of points must have same dimension as " "limits.")) for j in range(d): a, b = self.lims[0][j], self.lims[1][j] if self.scaling is not None: a, b = self.scaling[j](a), self.scaling[j](b) c, r = (a + b) / 2., np.abs(a - b) / 2. xj = c + r * np.exp(1.j * np.linspace(0, 2 * np.pi, n[j] + 1)[:-1, None]) wj = r / n[j] * np.ones(n[j]) if self.scalingInv is not None: xj = self.scalingInv[j](xj) if j == 0: x = xj w = wj xsize = n[0] else: x = np.concatenate((np.kron(np.ones(n[j])[:, None], x), np.kron(xj, np.ones(xsize)[:, None])), axis = 1) w = np.multiply(np.kron(np.ones(n[j]), w), np.kron(wj, np.ones(xsize))) xsize = xsize * n[j] return [y.flatten() for y in np.split(x, xsize)], w + def refine(self, n:int) -> Tuple[List[Np1D], Np1D]: + """ + Apply refinement. If points are not nested, equivalent to + generatePoints([2 * x for x in n]). + """ + super().generatePoints(n) + d = len(self.lims[0]) + try: + len(n) + except: + n = np.array([n]) + if len(n) != d: + raise Exception(("Numbers of points must have same dimension as " + "limits.")) + for j in range(d): + a, b = self.lims[0][j], self.lims[1][j] + if self.scaling is not None: + a, b = self.scaling[j](a), self.scaling[j](b) + c, r = (a + b) / 2., np.abs(a - b) / 2. + xj = c + r * np.exp(1.j * (np.pi / n[j] + + np.linspace(0, 2 * np.pi, n[j] + 1)[:-1, None])) + wj = r / n[j] * np.ones(n[j]) + if self.scalingInv is not None: + xj = self.scalingInv[j](xj) + if j == 0: + x = xj + w = wj + xsize = n[0] + else: + x = np.concatenate((np.kron(np.ones(n[j])[:, None], x), + np.kron(xj, np.ones(xsize)[:, None])), + axis = 1) + w = np.multiply(np.kron(np.ones(n[j]), w), + np.kron(wj, np.ones(xsize))) + xsize = xsize * n[j] + return [y.flatten() for y in np.split(x, xsize)], w diff --git a/rrompy/utilities/parameter_sampling/generic_sampler.py b/rrompy/utilities/parameter_sampling/generic_sampler.py index c3254f6..766d3cb 100644 --- a/rrompy/utilities/parameter_sampling/generic_sampler.py +++ b/rrompy/utilities/parameter_sampling/generic_sampler.py @@ -1,97 +1,107 @@ # Copyright (C) 2018 by the RROMPy authors # # This file is part of RROMPy. # # RROMPy is free software: you can redistribute it and/or modify # it under the terms of the GNU Lesser General Public License as published by # the Free Software Foundation, either version 3 of the License, or # (at your option) any later version. # # RROMPy is distributed in the hope that it will be useful, # but WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # GNU Lesser General Public License for more details. # # You should have received a copy of the GNU Lesser General Public License # along with RROMPy. If not, see . # from abc import abstractmethod import numpy as np from rrompy.utilities.base.types import Np1D, Tuple, GenExpr, List __all__ = ['GenericSampler'] class GenericSampler: """ABSTRACT. Generic generator of sample points.""" def __init__(self, lims:Tuple[Np1D, Np1D], scaling : List[callable] = None, scalingInv : List[callable] = None): self.lims = lims self.scaling = scaling self.scalingInv = scalingInv def name(self) -> str: return self.__class__.__name__ def __str__(self) -> str: return "{}[{}{}]".format(self.name(), np.array2string(self.lims[0], separator = '_'), np.array2string(self.lims[1], separator = '_')) def __repr__(self) -> str: return self.__str__() + " at " + hex(id(self)) def __eq__(self, other) -> bool: return self.__dict__ == other.__dict__ @property def lims(self): """Value of lims.""" return self._lims @lims.setter def lims(self, lims): if len(lims) != 2: raise Exception("2 limits must be specified.") try: lims = lims.tolist() except: lims = list(lims) for j in range(2): try: len(lims[j]) except: lims[j] = np.array([lims[j]]) if len(lims[0]) != len(lims[1]): raise Exception("The limits must have the same length.") self._lims = lims @property def scaling(self): """Value of scaling.""" return self._scaling @scaling.setter def scaling(self, scaling): if scaling is not None and not isinstance(scaling, list): scaling = [scaling] self._scaling = scaling @property def scalingInv(self): """Value of scalingInv.""" return self._scalingInv @scalingInv.setter def scalingInv(self, scalingInv): if scalingInv is not None and not isinstance(scalingInv, list): scalingInv = [scalingInv] self._scalingInv = scalingInv @abstractmethod def generatePoints(self, n:GenExpr) -> Tuple[List[Np1D], Np1D]: """Array of points and array of weights.""" assert ((self.scaling is None or len(self.scaling) == len(self.lims[0])) and (self.scalingInv is None or len(self.scalingInv) == len(self.lims[0]))) pass + def refine(self, n:int) -> Tuple[List[Np1D], Np1D]: + """ + Apply refinement. If points are not nested, equivalent to + generatePoints([2 * x - 1 for x in n]). + """ + try: + len(n) + except: + n = np.array([n]) + return self.generatePoints([2 * nj - 1 for nj in n]) diff --git a/rrompy/utilities/parameter_sampling/manual_sampler.py b/rrompy/utilities/parameter_sampling/manual_sampler.py index d06ea34..c069357 100644 --- a/rrompy/utilities/parameter_sampling/manual_sampler.py +++ b/rrompy/utilities/parameter_sampling/manual_sampler.py @@ -1,59 +1,81 @@ # Copyright (C) 2018 by the RROMPy authors # # This file is part of RROMPy. # # RROMPy is free software: you can redistribute it and/or modify # it under the terms of the GNU Lesser General Public License as published by # the Free Software Foundation, either version 3 of the License, or # (at your option) any later version. # # RROMPy is distributed in the hope that it will be useful, # but WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # GNU Lesser General Public License for more details. # # You should have received a copy of the GNU Lesser General Public License # along with RROMPy. If not, see . # import numpy as np from rrompy.utilities.parameter_sampling.generic_sampler import GenericSampler from rrompy.utilities.base.types import Np1D, Tuple, List from rrompy.utilities.warning_manager import warn __all__ = ['ManualSampler'] class ManualSampler(GenericSampler): """Manual generator of sample points.""" def __init__(self, lims:Tuple[Np1D, Np1D], points:Np1D, scaling : List[callable] = None, scalingInv : List[callable] = None): super().__init__(lims = lims, scaling = scaling, scalingInv = scalingInv) self.points = points def __str__(self) -> str: return "{}[{}]".format(self.name(), "_".join(map(str, self.points))) def __repr__(self) -> str: return self.__str__() + " at " + hex(id(self)) def generatePoints(self, n:int) -> Tuple[List[Np1D], Np1D]: """Array of quadrature points and array of weights.""" super().generatePoints(None) size = 1. / n for j in range(len(self.lims[0])): a, b = self.lims[0][j], self.lims[1][j] if self.scaling is not None: a, b = self.scaling[j](a), self.scaling[j](b) size *= np.abs(a - b) if n > len(self.points): warn(("Requested more points than given. Looping over first " "points.")) pts = np.tile(self.points, [np.int(np.ceil(len(self.points) / n))])[: n] else: pts = self.points[: n] return pts, np.ones(n) * size + def refine(self, n:int) -> Tuple[List[Np1D], Np1D]: + """ + Apply refinement. If points are not nested, equivalent to + generatePoints(2 * n - 1). + """ + super().generatePoints(None) + size = 1. / (n - 1) + for j in range(len(self.lims[0])): + a, b = self.lims[0][j], self.lims[1][j] + if self.scaling is not None: + a, b = self.scaling[j](a), self.scaling[j](b) + size *= np.abs(a - b) + if 2 * n - 1 > len(self.points): + warn(("Requested more points than given. Looping over first " + "points.")) + pts = np.tile(self.points, + [np.int(np.ceil(len(self.points) / (2 * n - 1)))] + )[n : 2 * n - 1] + else: + pts = self.points[n : 2 * n - 1] + return pts, np.ones(n - 1) * size + diff --git a/rrompy/utilities/parameter_sampling/quadrature_sampler.py b/rrompy/utilities/parameter_sampling/quadrature_sampler.py index 3d46842..a46c586 100644 --- a/rrompy/utilities/parameter_sampling/quadrature_sampler.py +++ b/rrompy/utilities/parameter_sampling/quadrature_sampler.py @@ -1,93 +1,130 @@ # Copyright (C) 2018 by the RROMPy authors # # This file is part of RROMPy. # # RROMPy is free software: you can redistribute it and/or modify # it under the terms of the GNU Lesser General Public License as published by # the Free Software Foundation, either version 3 of the License, or # (at your option) any later version. # # RROMPy is distributed in the hope that it will be useful, # but WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # GNU Lesser General Public License for more details. # # You should have received a copy of the GNU Lesser General Public License # along with RROMPy. If not, see . # import numpy as np from rrompy.utilities.parameter_sampling.generic_sampler import GenericSampler from rrompy.utilities.base.types import Np1D, Tuple, List __all__ = ['QuadratureSampler'] class QuadratureSampler(GenericSampler): """Generator of quadrature sample points.""" allowedKinds = ["UNIFORM", "CHEBYSHEV", "GAUSSLEGENDRE"] def __init__(self, lims:Tuple[Np1D, Np1D], kind : str = "UNIFORM", scaling : List[callable] = None, scalingInv : List[callable] = None): super().__init__(lims = lims, scaling = scaling, scalingInv = scalingInv) self.kind = kind def __str__(self) -> str: return "{}_{}".format(super().__str__(), self.kind) def __repr__(self) -> str: return self.__str__() + " at " + hex(id(self)) @property def kind(self): """Value of kind.""" return self._kind @kind.setter def kind(self, kind): if kind.upper() not in self.allowedKinds: raise Exception("Generator kind not recognized.") self._kind = kind.upper() def generatePoints(self, n:Np1D) -> Tuple[List[Np1D], Np1D]: """Array of quadrature points and array of weights.""" super().generatePoints(n) d = len(self.lims[0]) try: len(n) except: n = np.array([n]) if len(n) != d: raise Exception(("Numbers of points must have same dimension as" "limits.")) for j in range(d): a, b = self.lims[0][j], self.lims[1][j] if self.scaling is not None: a, b = self.scaling[j](a), self.scaling[j](b) if self.kind == "UNIFORM": xj = np.linspace(a, b, n[j])[:, None] - wj = np.abs(a - b) / n[j] * np.ones(n[j]) + wj = np.abs(a - b) / (n[j] - 1) * np.ones(n[j]) elif self.kind == "CHEBYSHEV": nodes, weights = np.polynomial.chebyshev.chebgauss(n[j]) xj = (a + b) / 2 + (a - b) / 2 * nodes[:, None] wj = np.abs(a - b) / np.pi * weights elif self.kind == "GAUSSLEGENDRE": nodes, weights = np.polynomial.legendre.leggauss(n[j]) xj = (a + b) / 2 + (a - b) / 2 * nodes[:, None] wj = np.abs(a - b) * weights if self.scalingInv is not None: xj = self.scalingInv[j](xj) if j == 0: x = xj w = wj xsize = n[0] else: x = np.concatenate((np.kron(np.ones(n[j])[:, None], x), np.kron(xj, np.ones(xsize)[:, None])), axis = 1) w = np.multiply(np.kron(np.ones(n[j]), w), np.kron(wj, np.ones(xsize))) xsize = xsize * n[j] return [y.flatten() for y in np.split(x, xsize)], w + def refine(self, n:int) -> Tuple[List[Np1D], Np1D]: + """ + Apply refinement. If points are not nested, equivalent to + generatePoints([2 * x - 1 for x in n]). + """ + if self.kind != "UNIFORM": return super().refine(n) + super().generatePoints(n) + d = len(self.lims[0]) + try: + len(n) + except: + n = np.array([n]) + if len(n) != d: + raise Exception(("Numbers of points must have same dimension as" + "limits.")) + for j in range(d): + a, b = self.lims[0][j], self.lims[1][j] + if self.scaling is not None: + a, b = self.scaling[j](a), self.scaling[j](b) + xj = np.linspace(a + (b - a) / 2. / (n[j] - 1), + b + (a - b) / 2. / (n[j] - 1), n[j] - 1)[:, None] + wj = np.abs(a - b) / (n[j] - 2) * np.ones(n[j] - 1) + if self.scalingInv is not None: + xj = self.scalingInv[j](xj) + if j == 0: + x = xj + w = wj + xsize = n[0] - 1 + else: + x = np.concatenate((np.kron(np.ones(n[j] - 1)[:, None], x), + np.kron(xj, np.ones(xsize)[:, None])), + axis = 1) + w = np.multiply(np.kron(np.ones(n[j] - 1), w), + np.kron(wj, np.ones(xsize))) + xsize = xsize * (n[j] - 1) + return [y.flatten() for y in np.split(x, xsize)], w + diff --git a/rrompy/utilities/parameter_sampling/random_sampler.py b/rrompy/utilities/parameter_sampling/random_sampler.py index 94dcc57..a61763d 100644 --- a/rrompy/utilities/parameter_sampling/random_sampler.py +++ b/rrompy/utilities/parameter_sampling/random_sampler.py @@ -1,70 +1,80 @@ # Copyright (C) 2018 by the RROMPy authors # # This file is part of RROMPy. # # RROMPy is free software: you can redistribute it and/or modify # it under the terms of the GNU Lesser General Public License as published by # the Free Software Foundation, either version 3 of the License, or # (at your option) any later version. # # RROMPy is distributed in the hope that it will be useful, # but WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # GNU Lesser General Public License for more details. # # You should have received a copy of the GNU Lesser General Public License # along with RROMPy. If not, see . # import numpy as np from rrompy.utilities.base.sobol import sobolGenerate from rrompy.utilities.parameter_sampling.generic_sampler import GenericSampler from rrompy.utilities.base.types import Np1D, Tuple, List __all__ = ['RandomSampler'] class RandomSampler(GenericSampler): """Generator of quadrature sample points.""" allowedKinds = ["UNIFORM", "SOBOL"] def __init__(self, lims:Tuple[Np1D, Np1D], kind : str = "UNIFORM", scaling : callable = None, scalingInv : callable = None): super().__init__(lims = lims, scaling = scaling, scalingInv = scalingInv) self.kind = kind def __str__(self) -> str: return "{}_{}".format(super().__str__(), self.kind) def __repr__(self) -> str: return self.__str__() + " at " + hex(id(self)) @property def kind(self): """Value of kind.""" return self._kind @kind.setter def kind(self, kind): if kind.upper() not in self.allowedKinds: raise Exception("Generator kind not recognized.") self._kind = kind.upper() def generatePoints(self, n:int, seed : int = 0) -> Tuple[List[Np1D], Np1D]: """Array of quadrature points and array of weights.""" super().generatePoints(n) d = len(self.lims[0]) if self.kind == "UNIFORM": np.random.seed(seed) x = np.random.uniform(size = (n, d)) else: x = sobolGenerate(d, n, seed) for j in range(d): a, b = self.lims[0][j], self.lims[1][j] if self.scaling is not None: a, b = self.scaling[j](a), self.scaling[j](b) x[:, j] = a + (b - a) * x[:, j] if self.scalingInv is not None: x[:, j] = self.scalingInv[j](x[:, j]) return [y.flatten() for y in np.split(x, n)], np.ones(n) + def refine(self, n:int, seed : int = 420) -> Tuple[List[Np1D], Np1D]: + """ + Apply refinement. If points are not nested, equivalent to + generatePoints([2 * x - 1 for x in n]). + """ + try: + len(n) + except: + n = np.array([n]) + return self.generatePoints([nj - 1 for nj in n], seed) diff --git a/rrompy/utilities/parameter_sampling/warping_sampler.py b/rrompy/utilities/parameter_sampling/warping_sampler.py index 84f0a9d..0f14302 100644 --- a/rrompy/utilities/parameter_sampling/warping_sampler.py +++ b/rrompy/utilities/parameter_sampling/warping_sampler.py @@ -1,114 +1,153 @@ # Copyright (C) 2018 by the RROMPy authors # # This file is part of RROMPy. # # RROMPy is free software: you can redistribute it and/or modify # it under the terms of the GNU Lesser General Public License as published by # the Free Software Foundation, either version 3 of the License, or # (at your option) any later version. # # RROMPy is distributed in the hope that it will be useful, # but WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # GNU Lesser General Public License for more details. # # You should have received a copy of the GNU Lesser General Public License # along with RROMPy. If not, see . # import numpy as np from rrompy.utilities.parameter_sampling.generic_sampler import GenericSampler from rrompy.utilities.base.types import Np1D, Tuple, List __all__ = ['WarpingSampler', 'WarpingFunction'] class WarpingFunction: """Wrapper for warping function.""" def __init__(self, other : "WarpingFunction" = None, **kwargs): if other is not None: self._call_ = other._call_ self._repr_ = other._repr_ else: self._call_ = kwargs["call"] self._repr_ = kwargs["repr"] if not callable(self._call_): self._call_ = lambda x: self._call_ if callable(self._repr_): self._repr_ = self._repr_() def __call__(self, x): return self._call_(x) def __str__(self): return self._repr_ def __repr__(self) -> str: return self.__str__() + " at " + hex(id(self)) class WarpingSampler(GenericSampler): """Generator of sample points from warping of uniform ones.""" def __init__(self, lims:Tuple[Np1D, Np1D], warping : List[callable], scaling : List[callable] = None, scalingInv : List[callable] = None): super().__init__(lims = lims, scaling = scaling, scalingInv = scalingInv) self.warping = warping def __str__(self) -> str: return "{}_{}".format(super().__str__(), self.warping) def __repr__(self) -> str: return self.__str__() + " at " + hex(id(self)) @property def warping(self): """Value of warping.""" return self._warping @warping.setter def warping(self, warping): if not isinstance(warping, (list, tuple,)): warping = [warping] * len(self.lims[0]) for j in range(len(warping)): if not isinstance(warping[j], (WarpingFunction,)): warping[j] = WarpingFunction(call = warping[j].__call__, repr = warping[j].__repr__) self._warping = warping def generatePoints(self, n:Np1D) -> Tuple[List[Np1D], Np1D]: """Array of quadrature points and array of weights.""" super().generatePoints(n) assert len(self.warping) == len(self.lims[0]) d = len(self.lims[0]) try: len(n) except: n = np.array([n]) if len(n) != d: raise Exception(("Numbers of points must have same dimension as" "limits.")) for j in range(d): a, b = self.lims[0][j], self.lims[1][j] if self.scaling is not None: a, b = self.scaling[j](a), self.scaling[j](b) x0, sigma = (a + b) / 2, (a - b) / 2 xj0 = np.linspace(- 1., 1., n[j])[:, None] xj = x0 + sigma * self.warping[j](xj0) - wj = np.abs(a - b) / n[j] * np.ones(n[j]) + wj = np.abs(a - b) / (n[j] - 1) * np.ones(n[j]) if self.scalingInv is not None: xj = self.scalingInv[j](xj) if j == 0: x = xj w = wj xsize = n[0] else: x = np.concatenate((np.kron(np.ones(n[j])[:, None], x), np.kron(xj, np.ones(xsize)[:, None])), axis = 1) w = np.multiply(np.kron(np.ones(n[j]), w), np.kron(wj, np.ones(xsize))) xsize = xsize * n[j] return [y.flatten() for y in np.split(x, xsize)], w + def refine(self, n:int) -> Tuple[List[Np1D], Np1D]: + """ + Apply refinement. If points are not nested, equivalent to + generatePoints([2 * x - 1 for x in n]). + """ + super().generatePoints(n) + assert len(self.warping) == len(self.lims[0]) + d = len(self.lims[0]) + try: + len(n) + except: + n = np.array([n]) + if len(n) != d: + raise Exception(("Numbers of points must have same dimension as" + "limits.")) + for j in range(d): + a, b = self.lims[0][j], self.lims[1][j] + if self.scaling is not None: + a, b = self.scaling[j](a), self.scaling[j](b) + x0, sigma = (a + b) / 2, (a - b) / 2 + xj0 = np.linspace(1. / (n[j] - 1) - 1., 1. - 1. / (n[j] - 1), + n[j] - 1)[:, None] + xj = x0 + sigma * self.warping[j](xj0) + wj = np.abs(a - b) / (n[j] - 2) * np.ones(n[j] - 1) + if self.scalingInv is not None: + xj = self.scalingInv[j](xj) + if j == 0: + x = xj + w = wj + xsize = n[0] - 1 + else: + x = np.concatenate((np.kron(np.ones(n[j] - 1)[:, None], x), + np.kron(xj, np.ones(xsize)[:, None])), + axis = 1) + w = np.multiply(np.kron(np.ones(n[j] - 1), w), + np.kron(wj, np.ones(xsize))) + xsize = xsize * (n[j] - 1) + return [y.flatten() for y in np.split(x, xsize)], w +