diff --git a/examples/8_damped_mass_chain/damped_mass_chain.py b/examples/8_damped_mass_chain/damped_mass_chain.py index 9b4f09c..a2cda76 100644 --- a/examples/8_damped_mass_chain/damped_mass_chain.py +++ b/examples/8_damped_mass_chain/damped_mass_chain.py @@ -1,157 +1,180 @@ ### example from Lohmann, Eid. Efficient Order Reduction of Parametric and ### Nonlinear Models by Superposition of Locally Reduced Models. from copy import deepcopy as copy import numpy as np import matplotlib.pyplot as plt from rrompy.reduction_methods import (NearestNeighbor as NN, RationalInterpolant as RI, RationalInterpolantGreedy as RIG, RationalInterpolantPivoted as RIP, RationalInterpolantGreedyPivoted as RIGP) from rrompy.parameter.parameter_sampling import (QuadratureSampler as QS, EmptySampler as ES) +from damped_mass_chain_engine import (bode as bode0, bodeLog, MassChainEngine, + MassChainEngineLog, AugmentedMassChainEngine, AugmentedMassChainEngineLog) +from rrompy.utilities.base.decorators import addWhiteNoise + +########################## +fullModelOrder = 1 #+ 1 +SMarginal = 1 #* 2 +state = 0 #+ 1 +noise_level = 0 #+ 1e-5 +########################## + +modelSign = "Surrogate modeling for frequency response of " +if fullModelOrder == 1: modelSign += "augmented " +modelSign += "damper-mass-spring model" +if SMarginal > 1: modelSign += " with 1 design parameter" +modelSign += ".\nOutput is " +if state: + modelSign += "vector of mass displacements. " +else: + modelSign += "displacement of last mass. " +modelSign += "Noise level: {}.\n".format(noise_level) +print(modelSign) -N, order, SMarginal = 4, 1, 2 M = [np.array([1., 5., 25., 125.])] +N = len(M[0]) D = [np.zeros((N, N))] D[0][0, 1], D[0][1, 2], D[0][2, 3], D[0][3, 3] = .1, .4, 1.6, 0. D[0] = D[0] + D[0].T K = [np.zeros((N, N))] K[0][0, 1], K[0][1, 2], K[0][2, 3], K[0][0, 3], K[0][3, 3] = 9., 3., 1., 1., 2. K[0] = K[0] + K[0].T -B, C = np.zeros((N, 1)), np.zeros((1, N)) -B[0, 0], C[0, 3] = 27., 1. +B = np.append(27., np.zeros(N - 1)).reshape(-1, 1) if SMarginal > 1: M += [np.zeros(N)] D += [np.zeros((N, N))] D[1][3, 3] = 1. D[1] = D[1] + D[1].T K += [np.zeros((N, N))] K[1][0, 3], K[1][3, 3] = 2., 2. K[1] = K[1] + K[1].T +if state: + C = np.eye(4) +else: + C = np.append(np.zeros(N - 1), 1.).reshape(1, -1) for logspace in range(2): print("Approximation in l{}space".format("og" * logspace + "in" * (not logspace))) if logspace: - from damped_mass_chain_engine import bodeLog as bode - if order == 1: - from damped_mass_chain_engine import ( - AugmentedMassChainEngineLog as engine) + bode = bodeLog + if fullModelOrder == 1: + engine = AugmentedMassChainEngineLog else: - from damped_mass_chain_engine import MassChainEngineLog as engine + engine = MassChainEngineLog else: - from damped_mass_chain_engine import bode - if order == 1: - from damped_mass_chain_engine import ( - AugmentedMassChainEngine as engine) + bode = bode0 + if fullModelOrder == 1: + engine = AugmentedMassChainEngine else: - from damped_mass_chain_engine import MassChainEngine as engine - solver = engine(M, D, K, B, C) + engine = MassChainEngine + solver = addWhiteNoise(noise_level)(engine)(M, D, K, B, C) ss, mu = [1e-2, 1e1], [] s0 = 10. ** np.mean(np.log10(ss)) freq = np.logspace(np.log10(ss[0]), np.log10(ss[1]), 100) if logspace: ss, freq = [np.log10(ss[0]), np.log10(ss[1])], np.log10(freq) s0, parameterMap = np.log10(s0), 1. else: parameterMap = {"F": [("log10", "x")], "B": [(10., "**", "x")]} krange = [[ss[0]], [ss[-1]]] k0, srange = [s0], copy(krange) if SMarginal > 1: ms = [0., 1.] m0, mrange = np.mean(ms), [[ms[0]], [ms[-1]]] krange[0] += mrange[0] krange[1] += mrange[1] k0 += [m0] mu = [.5 * (ms[1] - ms[0]) / (SMarginal - 1)] if not logspace: parameterMap["F"] += [("x")] parameterMap["B"] += [("x")] for method in ["RI", "RI_GREEDY"]: print("Testing {} method".format(method)) if method == "RI": - params = {'S':8, 'POD':True, 'polybasis':"CHEBYSHEV"} + params = {'S':15, 'N': 8, 'POD':True, 'polybasis':"CHEBYSHEV"} if SMarginal > 1: algo = RIP else: params['sampler'] = QS(srange, "CHEBYSHEV", parameterMap) algo = RI if method == "RI_GREEDY": params = {'S':5, 'POD':True, 'polybasis':"LEGENDRE", 'greedyTol':1e-2, 'errorEstimatorKind':"DISCREPANCY", 'trainSetGenerator':QS(srange, "CHEBYSHEV", parameterMap)} if SMarginal > 1: algo = RIGP else: params['sampler'] = QS(srange, "UNIFORM", parameterMap) algo = RIG if SMarginal > 1: params["paramsMarginal"] = {"MMarginal": SMarginal - 1} params['SMarginal'] = SMarginal params['polybasisMarginal'] = "MONOMIAL" params['radialDirectionalWeightsMarginal'] = [2. / (ms[1] - ms[0])] params['matchingWeight'] = 1. #params['cutOffTolerance'] = 2. params['samplerPivot'] = QS(srange, "UNIFORM", parameterMap) params['samplerMarginal'] = QS(mrange, "UNIFORM") approx = algo([0], solver, mu0 = k0, approx_state = True, approxParameters = params, verbosity = 5, storeAllSamples = True) else: approx = algo(solver, mu0 = k0, approx_state = True, approxParameters = params, verbosity = 5) if "GREEDY" in method: approx.setupApprox("LAST") else: approx.setupApprox() approxNN = NN(solver, mu0 = k0, approx_state = True, verbosity = 5, - approxParameters = {'S':len(approx.mus), 'POD':True, - 'sampler':ES()}) + approxParameters = {'S':len(approx.mus), + 'POD':params['POD'], 'sampler':ES()}) if SMarginal > 1: approxNN.setSamples(approx.storedSamplesFilenames) approx.purgeStoredSamples() for m in approx.musMarginal: bode(freq, m[0], [approx.getHF, approx.getApprox, approxNN.getApprox]) else: approxNN.setSamples(approx.samplingEngine) bode(freq, mu, [approx.getHF, approx.getApprox, approxNN.getApprox]) if SMarginal > 1: bode(freq, [1.5 * ms[1]], [approx.getHF, approx.getApprox, approxNN.getApprox]) bode(freq, [2. * ms[1]], [approx.getHF, approx.getApprox, approxNN.getApprox]) verb = approx.verbosity approx.verbosity = 0 mspace = np.linspace(ms[0], ms[-1], 10) for j, t in enumerate(mspace): pls = approx.getPoles([None, t]) if j == 0: poles = np.empty((len(mspace), len(pls)), dtype = np.complex) poles[j] = pls for j, t in enumerate(approx.musMarginal): pls = approx.getPoles([None, t[0][0]]) if j == 0: polesE = np.empty((SMarginal, len(pls)), dtype = np.complex) polesE[j] = pls approx.verbosity = verb fig = plt.figure(figsize = (10, 6)) ax = fig.add_subplot(1, 1, 1) ax.plot(np.real(poles), np.imag(poles), '--') ax.plot(np.real(polesE), np.imag(polesE), 'ko', markersize = 4) ax.set_xlabel('Real') ax.set_ylabel('Imag') ax.grid() plt.show() else: poles = approx.getPoles() print("Poles:\n{}".format(poles)) print("\n") diff --git a/rrompy/sampling/engines/pod_engine.py b/rrompy/sampling/engines/pod_engine.py index c9ebe8f..1f146bd 100644 --- a/rrompy/sampling/engines/pod_engine.py +++ b/rrompy/sampling/engines/pod_engine.py @@ -1,144 +1,137 @@ # Copyright (C) 2018 by the RROMPy authors # # This file is part of RROMPy. # # RROMPy is free software: you can redistribute it and/or modify # it under the terms of the GNU Lesser General Public License as published by # the Free Software Foundation, either version 3 of the License, or # (at your option) any later version. # # RROMPy is distributed in the hope that it will be useful, # but WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # GNU Lesser General Public License for more details. # # You should have received a copy of the GNU Lesser General Public License # along with RROMPy. If not, see . # import numpy as np from copy import deepcopy as copy from warnings import catch_warnings from rrompy.utilities.base.types import Np1D, Np2D, Tuple, HFEng, sampList from rrompy.utilities.numerical import dot from rrompy.sampling import sampleList __all__ = ['PODEngine'] class PODEngine: """ POD engine for general matrix orthogonalization. """ def __init__(self, HFEngine:HFEng): self.HFEngine = HFEngine 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 GS(self, a:Np1D, Q:sampList, n : int = -1, is_state : bool = True) -> Tuple[Np1D, Np1D, bool]: """ Compute 1 Gram-Schmidt step with given projector. Args: a: vector to be projected; Q: orthogonal projection matrix; n: number of columns of Q to be considered; is_state: whether state-norm should be used. Returns: Resulting normalized vector, coefficients of a wrt the updated basis, whether computation is ill-conditioned. """ if n == -1: n = Q.shape[1] r = np.zeros((n + 1,), dtype = Q.dtype) if n > 0: Q = Q[: n] for j in range(2): # twice is enough! nu = self.HFEngine.innerProduct(a, Q, is_state = is_state) a = a - dot(Q, nu) r[:-1] = r[:-1] + nu.flatten() r[-1] = self.HFEngine.norm(a, is_state = is_state) ill_cond = False with catch_warnings(record = True) as w: snr = np.abs(r[-1]) / np.linalg.norm(r) if len(w) > 0 or np.isclose(snr, 0.): ill_cond = True r[-1] = 1. a = a / r[-1] return a, r, ill_cond def generalizedQR(self, A:sampList, Q0 : sampList = None, only_R : bool = False, genTrials : int = 10, is_state : bool = True) -> Tuple[sampList, Np2D]: """ Compute generalized QR decomposition of a matrix through Householder - method. + method; see Trefethen, IMA J.N.A., 2010. Args: A: matrix to be decomposed; Q0(optional): initial orthogonal guess for Q; defaults to random; only_R(optional): whether to skip reconstruction of Q; defaults to False. genTrials(optional): number of trials of generation of linearly independent vector; defaults to 10. is_state(optional): whether state-norm should be used; defaults to True. Returns: Resulting (orthogonal and )upper-triangular factor(s). """ Nh, N = A.shape B = copy(A) - V = copy(A) + V = sampleList(np.zeros(A.shape, dtype = A.dtype)) R = np.zeros((N, N), dtype = A.dtype) - if Q0 is None: - Q = sampleList(np.zeros(A.shape, dtype = A.dtype) - + np.random.randn(*(A.shape))) - else: - Q = copy(Q0) + Q = copy(V) if Q0 is None else copy(Q0) for k in range(N): a = B[k] R[k, k] = self.HFEngine.norm(a, is_state = is_state) if Q0 is None and k < Nh: for _ in range(genTrials): Q[k], _, illC = self.GS(np.random.randn(Nh), Q, k, is_state) if not illC: break else: - illC = False - if illC or k >= Nh: - Q[k] = np.zeros(Nh, dtype = Q.dtype) - alpha = 0. + illC = k >= Nh + if illC: + if Q0 is not None or k < Nh: Q.data[:, k] = 0. else: alpha = self.HFEngine.innerProduct(a, Q[k], is_state = is_state) if np.isclose(np.abs(alpha), 0.): s = 1. - else: s = alpha / np.abs(alpha) - Q[k] = s * Q[k] - V[k], _, _ = self.GS(R[k, k] * Q[k] - a, Q, k, is_state) + else: s = - alpha / np.abs(alpha) + Q.data[:, k] *= s + V[k], _, _ = self.GS(R[k, k] * Q[k] - a, Q, k, is_state) J = np.arange(k + 1, N) vtB = self.HFEngine.innerProduct(B[J], V[k], is_state = is_state) B.data[:, J] -= 2 * np.outer(V[k], vtB) - if illC: - R[k, J] = 0. - else: + if not illC: R[k, J] = self.HFEngine.innerProduct(B[J], Q[k], is_state = is_state) B.data[:, J] -= np.outer(Q[k], R[k, J]) if only_R: return R - for k in range(N - 1, -1, -1): - J = list(range(k, N)) + for k in range(min(Nh, N) - 1, -1, -1): + J = np.arange(k, min(Nh, N)) vtQ = self.HFEngine.innerProduct(Q[J], V[k], is_state = is_state) Q.data[:, J] -= 2 * np.outer(V[k], vtQ) return Q, R