Homec4science
Diffusion RationalROMPy (master) Browse

RROMPy/rrompy/hfengines/base/marginal_proxy_engine.py
59aa187eb896master

Last Change
Raw File

marginal_proxy_engine.py
View Options

# 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 <http://www.gnu.org/licenses/>.
#
import numpy as np
import scipy.sparse as scsp
from copy import deepcopy as copy, copy as softcopy
from rrompy.utilities.base.types import (Np1D, Np2D, ScOp, List, paramVal,
paramList, HFEng, sampList)
from rrompy.utilities.base import freepar as fp
from rrompy.utilities.numerical import multibinom, marginalizePolyList
from rrompy.utilities.poly_fitting.polynomial import (
hashDerivativeToIdx as hashD, hashIdxToDerivative as hashI)
from rrompy.utilities.exception_manager import RROMPyAssert
from rrompy.parameter import checkParameter, checkParameterList
from rrompy.sampling import sampleList, emptySampleList
from rrompy.solver.scipy import tensorizeLS, detensorizeLS
__all__ = ['MarginalProxyEngine']
class MarginalProxyEngine:
"""
Marginalized should prescribe fixed value for the marginalized parameters
and leave freepar/None elsewhere.
"""
def __init__(self, HFEngine:HFEng, marginalized:Np1D):
self.baseHF = HFEngine
self.marg = marginalized
self._setupAs()
self._setupbs()
@property
def freeLocations(self):
return tuple([x for x in range(self.nparBase) if self.marg[x] == fp])
@property
def fixedLocations(self):
return tuple([x for x in range(self.nparBase) if self.marg[x] != fp])
@property
def muFixed(self):
muF = checkParameter([m for m in self.marg if m != fp])
if self.nparBase - self.npar > 0: muF = muF[0]
return muF
@property
def rescalingExp(self):
return [self.baseHF.rescalingExp[x] for x in self.freeLocations]
@property
def rescalingExpFixed(self):
return [self.baseHF.rescalingExp[x] for x in self.fixedLocations]
@property
def V(self):
return self.baseHF.V
def name(self) -> str:
return "{}-proxy for {}".format(self.freeLocations, self.baseHF.name())
def __str__(self) -> str:
return self.name()
def __repr__(self) -> str:
return self.__str__() + " at " + hex(id(self))
def __dir_base__(self):
return [x for x in self.__dir__() if x[:2] != "__"]
def __deepcopy__(self, memo):
return softcopy(self)
@property
def npar(self):
"""Value of npar."""
return len(self.freeLocations)
@property
def nparBase(self):
"""Value of npar."""
return self.baseHF.npar
@property
def nAs(self):
"""Value of nAs."""
return self._nAs
@property
def nbs(self):
"""Value of nbs."""
return self._nbs
@property
def nbsH(self) -> int:
return max(self.nbs, self.nAs)
def spacedim(self):
return self.As[0].shape[1]
def checkParameter(self, mu:paramList):
return checkParameter(mu, self.npar)
def checkParameterList(self, mu:paramList):
return checkParameterList(mu, self.npar)
def buildEnergyNormForm(self):
self.baseHF.buildEnergyNormForm()
self.energyNormMatrix = self.baseHF.energyNormMatrix
def buildEnergyNormDualForm(self):
self.baseHF.buildEnergyNormDualForm()
self.energyNormDualMatrix = self.baseHF.energyNormDualMatrix
def buildDualityPairingForm(self):
self.baseHF.buildDualityPairingForm()
self.dualityMatrix = self.baseHF.dualityMatrix
def buildEnergyNormPartialDualForm(self):
self.baseHF.buildEnergyNormPartialDualForm()
self.energyNormPartialDualMatrix = (
self.baseHF.energyNormPartialDualMatrix)
def innerProduct(self, u:Np2D, v:Np2D, onlyDiag : bool = False,
dual : bool = False, duality : bool = True) -> Np2D:
return self.baseHF.innerProduct(u, v, onlyDiag, dual, duality)
def norm(self, u:Np2D, dual : bool = False, duality : bool = True) -> Np1D:
return self.baseHF.norm(u, dual, duality)
def checkAInBounds(self, derI : int = 0):
"""Check if derivative index is oob for operator of linear system."""
if derI < 0 or derI >= self.nAs:
d = self.spacedim()
return scsp.csr_matrix((np.zeros(0), np.zeros(0), np.zeros(d + 1)),
shape = (d, d), dtype = np.complex)
def checkbInBounds(self, derI : int = 0, homogeneized : bool = False):
"""Check if derivative index is oob for RHS of linear system."""
nbs = self.nbsH if homogeneized else self.nbs
if derI < 0 or derI >= nbs:
return np.zeros(self.spacedim(), dtype = np.complex)
def _assembleA(self, mu : paramVal = [], der : List[int] = 0,
derI : int = None) -> ScOp:
"""Assemble (derivative of) operator of linear system."""
mu = self.checkParameter(mu)
if self.npar > 0: mu = mu[0]
if not hasattr(der, "__len__"): der = [der] * self.npar
if derI is None: derI = hashD(der)
Anull = self.checkAInBounds(derI)
if Anull is not None: return Anull
rExp = self.rescalingExp
A = copy(self.As[derI])
for j in range(derI + 1, self.nAs):
derIdx = hashI(j, self.npar)
diffIdx = [x - y for (x, y) in zip(derIdx, der)]
if np.all([x >= 0 for x in diffIdx]):
A = A + (np.prod((mu ** rExp) ** diffIdx)
* multibinom(derIdx, diffIdx) * self.As[j])
return A
def A(self, mu : paramVal = [], der : List[int] = 0) -> ScOp:
"""
Assemble terms of operator of linear system and return it (or its
derivative) at a given parameter.
"""
mu = self.checkParameter(mu)
if not hasattr(der, "__len__"): der = [der] * self.npar
derI = hashD(der)
return self._assembleA(mu, der, derI)
def _setupAs(self):
mVals = [fp] * self.nparBase
muFixedEff = self.muFixed ** self.rescalingExpFixed
for i, mFE in zip(self.fixedLocations, muFixedEff):
mVals[i] = mFE
self.As = marginalizePolyList(self.baseHF.As, mVals, "auto")
self._nAs = len(self.As)
def affineLinearSystemA(self, mu : paramVal = []) -> List[Np2D]:
"""
Assemble affine blocks of operator of linear system (just linear
blocks).
"""
As = [None] * self.nAs
for j in range(self.nAs):
As[j] = self.A(mu, hashI(j, self.npar))
return As
def _assembleb(self, mu : paramVal = [], der : List[int] = 0,
derI : int = None, homogeneized : bool = False) -> ScOp:
"""Assemble (derivative of) (homogeneized) RHS of linear system."""
mu = self.checkParameter(mu)
if self.npar > 0: mu = mu[0]
if not hasattr(der, "__len__"): der = [der] * self.npar
if derI is None: derI = hashD(der)
bnull = self.checkbInBounds(derI, homogeneized)
if bnull is not None: return bnull
bs = self.bsH if homogeneized else self.bs
rExp = self.rescalingExp
b = copy(bs[derI])
for j in range(derI + 1, len(bs)):
derIdx = hashI(j, self.npar)
diffIdx = [x - y for (x, y) in zip(derIdx, der)]
if np.all([x >= 0 for x in diffIdx]):
b = b + (np.prod((mu ** rExp) ** diffIdx)
* multibinom(derIdx, diffIdx) * bs[j])
return b
def b(self, mu : paramVal = [], der : List[int] = 0,
homogeneized : bool = False) -> Np1D:
"""
Assemble terms of (homogeneized) RHS of linear system and return it (or
its derivative) at a given parameter.
"""
mu = self.checkParameter(mu)
if not hasattr(der, "__len__"): der = [der] * self.npar
derI = hashD(der)
return self._assembleb(mu, der, derI, homogeneized)
def _setupbs(self):
mVals = [fp] * self.nparBase
muFixedEff = self.muFixed ** self.rescalingExpFixed
for i, mFE in zip(self.fixedLocations, muFixedEff):
mVals[i] = mFE
self.bs = marginalizePolyList(self.baseHF.bs, mVals, "auto")
self._nbs = len(self.bs)
def affineLinearSystemb(self, mu : paramVal = [],
homogeneized : bool = False) -> List[Np1D]:
"""
Assemble affine blocks of RHS of linear system (just linear blocks).
"""
nbs = self.nbsH if homogeneized else self.nbs
bs = [None] * nbs
for j in range(nbs):
bs[j] = self.b(mu, hashI(j, self.npar), homogeneized)
return bs
def solve(self, mu : paramList = [], RHS : sampList = None,
homogeneized : bool = False) -> sampList:
"""
Find solution of linear system.
Args:
mu: parameter value.
RHS: RHS of linear system. If None, defaults to that of parametric
system. Defaults to None.
"""
mu = self.checkParameterList(mu)[0]
if mu.shape[0] == 0: mu.reset((1, self.npar), mu.dtype)
sol = emptySampleList()
if len(mu) > 0:
if RHS is None:
RHS = [self.b(m, homogeneized = homogeneized) for m in mu]
RHS = sampleList(RHS)
mult = 0 if len(RHS) == 1 else 1
RROMPyAssert(mult * (len(mu) - 1) + 1, len(RHS), "Sample size")
u = self.baseHF._solver(self.A(mu[0]), RHS[0],
self.baseHF._solverArgs)
sol.reset((len(u), len(mu)), dtype = u.dtype)
sol[0] = u
for j in range(1, len(mu), self.baseHF._solveBatchSize):
if self.baseHF._solveBatchSize != 1:
iRange = list(range(j, min(j + self.baseHF._solveBatchSize,
len(mu))))
As = [self.A(mu[i]) for i in iRange]
bs = [RHS[mult * i] for i in iRange]
A, b = tensorizeLS(As, bs)
else:
A, b = self.A(mu[j]), RHS[mult * j]
solStack = self.baseHF._solver(A, b, self.baseHF._solverArgs)
if self.baseHF._solveBatchSize != 1:
sol[iRange] = detensorizeLS(solStack, len(iRange))
else:
sol[j] = solStack
return sol
def residual(self, u:sampList, mu : paramList = [],
homogeneized : bool = False,
duality : bool = True) -> sampList:
"""
Find residual of linear system for given approximate solution.
Args:
u: numpy complex array with function dofs. If None, set to 0.
mu: parameter value.
"""
mu = self.checkParameterList(mu)[0]
if mu.shape[0] == 0: mu.reset((1, self.npar), mu.dtype)
if u is not None:
u = sampleList(u)
mult = 0 if len(u) == 1 else 1
RROMPyAssert(mult * (len(mu) - 1) + 1, len(u), "Sample size")
res = emptySampleList()
if duality and not hasattr(self, "dualityMatrix"):
self.buildDualityPairingForm()
for j in range(len(mu)):
b = self.b(mu[j], homogeneized = homogeneized)
if u is None:
r = b
else:
r = b - self.A(mu[j]).dot(u[mult * j])
if j == 0:
res.reset((len(r), len(mu)), dtype = r.dtype)
if duality:
r = self.dualityMatrix.dot(r)
res[j] = r
return res
def _rayleighQuotient(self, *args, **kwargs) -> float:
return self.baseHF._rayleighQuotient(*args, **kwargs)
def stabilityFactor(self, u:sampList, mu : paramList = [],
nIterP : int = 10, nIterR : int = 10) -> sampList:
"""
Find stability factor of matrix of linear system using iterative
inverse power iteration- and Rayleigh quotient-based procedure.
Args:
u: numpy complex arrays with function dofs.
mu: parameter values.
nIterP: number of iterations of power method.
nIterR: number of iterations of Rayleigh quotient method.
"""
mu = self.checkParameterList(mu)[0]
if mu.shape[0] == 0: mu.reset((1, self.npar), mu.dtype)
u = sampleList(u)
mult = 0 if len(u) == 1 else 1
RROMPyAssert(mult * (len(mu) - 1) + 1, len(u), "Sample size")
stabFact = np.empty(len(mu), dtype = float)
if not hasattr(self, "energyNormMatrix"):
self.buildEnergyNormForm()
for j in range(len(mu)):
stabFact[j] = self._rayleighQuotient(self.A(mu[j]), u[mult * j],
self.energyNormMatrix,
0., nIterP, nIterR)
return stabFact
def plot(self, *args, **kwargs):
"""
Do some nice plots of the complex-valued function with given dofs.
"""
self.baseHF.plot(*args, **kwargs)
def plotmesh(self, *args, **kwargs):
"""
Do a nice plot of the mesh.
"""
self.baseHF.plotmesh(*args, **kwargs)
def outParaview(self, *args, **kwargs):
"""
Output complex-valued function with given dofs to ParaView file.
"""
self.baseHF.outParaview(*args, **kwargs)
def outParaviewTimeDomain(self, *args, **kwargs):
"""
Output complex-valued function with given dofs to ParaView file,
converted to time domain.
"""
self.baseHF.outParaviewTimeDomain(*args, **kwargs)
c4science · Help