Homec4science
Diffusion RationalROMPy (develop) Browse

RROMPy/rrompy/utilities/poly_fitting/heaviside/heaviside_to_from_affine.py
3cacafea472fdevelop

Last Change
Raw File

heaviside_to_from_affine.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
from scipy.special import binom
import scipy.sparse as sp
from rrompy.utilities.base.types import (Np1D, Np2D, List, ListAny, Tuple,
paramVal)
from rrompy.utilities.numerical import dot, solve
from rrompy.utilities.numerical.nonlinear_eigenproblem import eigNonlinearDense
from rrompy.utilities.exception_manager import RROMPyException
from rrompy.parameter import checkParameter
__all__ = ['heaviside2affine', 'affine2heaviside']
def heaviside2affine(c:Np2D, poles:Np1D, mu : paramVal = [],
basis : str = "MONOMIAL_HEAVISIDE",
sparse : bool = False) \
-> Tuple[Np2D, List[Np2D], List[Np1D]]:
mu = checkParameter(mu, 1)(0, 0)
n, d = len(poles), len(c) - len(poles)
basisN = basis.split("_")[0]
if basisN not in ["MONOMIAL", "CHEBYSHEV", "LEGENDRE"]:
raise RROMPyException("Polynomial basis not recognized.")
if sparse:
A0 = sp.spdiags([np.concatenate((- mu - poles, np.ones(d)))],
[0], n + d, n + d)
A1 = sp.spdiags([np.concatenate((np.ones(n), np.zeros(d)))],
[0], n + d, n + d)
else:
A0 = np.diag(np.concatenate((mu - poles, np.ones(d))))
A1 = np.diag(np.concatenate((np.ones(n), np.zeros(d))))
As = [A0, A1]
bs = np.zeros((d, n + d), dtype = poles.dtype)
bs[0, :] = 1.
if d > 0:
bs[0, n + 1 :] = 0.
if d > 1:
bs[1, n + 1] = 1.
for j in range(2, d):
if basisN == "MONOMIAL":
bs[j, n + j] = 1.
else:
alpha = - 1. if basisN == "CHEBYSHEV" else 1. / j - 1.
bs[:, n + j] = alpha * bs[:, n + j - 2]
bs[1 :, n + j] += (1. - alpha) * bs[: -1, n + j - 1]
bs = list(bs)
return c.reshape(c.shape[0], -1).T, As, bs
def affine2heaviside(As:ListAny, bs:ListAny,
jSupp : int = 1) -> Tuple[Np2D, Np1D, str]:
if jSupp != 1 and not (isinstance(jSupp, (int,))
and jSupp.upper() == "COMPANION"):
raise RROMPyException(("Affine to heaviside conversion does not allow "
"nonlinear eigenproblem support outside first "
"block row."))
N = len(As)
M = len(bs)
n = As[0].shape[0]
if N == 1:
poles = np.empty(0, dtype = np.complex)
Q = np.eye(n)
else:
basis = "MONOMIAL_HEAVISIDE"
poles, P, Q = eigNonlinearDense(As, jSupp = jSupp,
return_inverse = True)
P = P[- n :, :]
Q = Q[:, : n]
bEffs = np.array([dot(Q, solve(As[-1], b, np.linalg.solve,
{})) for b in bs])
if N == 1:
c = bEffs
else:
c = np.zeros((len(poles) + M - 1, As[0].shape[1]), dtype = np.complex)
for l, pl in enumerate(poles):
for i in range(M):
c[l, :] = pl ** i * bEffs[i, l] * P[:, l]
for l in range(M - 1):
for i in range(l + 1, M):
c[len(poles) + l, :] = dot(P, poles ** (i- 1 - l) * bEffs[i, :])
return c, poles, basis
c4science · Help