Page MenuHomec4science
Files F60628375

vander.py
No OneTemporary
Actions

File Metadata

Created
Wed, May 1, 14:02

vander.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 .kernel import kernels
from rrompy.utilities.poly_fitting.polynomial.vander import polyvander as pvP
from rrompy.utilities.base.types import Np1D, Np2D, Tuple, List, paramList
from rrompy.parameter import checkParameterList
from rrompy.utilities.exception_manager import RROMPyException, RROMPyAssert
__all__ = ['rbvander', 'polyvander']
def rbvander(x:paramList, basis:str, reorder : List[int] = None,
directionalWeights : Np1D = None) -> Np2D:
"""Compute radial-basis-Vandermonde matrix."""
x = checkParameterList(x)
x_un = x.unique()
nx = len(x)
if len(x_un) < nx:
raise RROMPyException("Sample points must be distinct.")
del x_un
x = x.data
if directionalWeights is None:
directionalWeights = np.ones(x.shape[1])
elif not hasattr(directionalWeights, "__len__"):
directionalWeights = directionalWeights * np.ones(x.shape[1])
RROMPyAssert(len(directionalWeights), x.shape[1],
"Number of directional weights")
try:
radialker = kernels[basis.upper()]
except:
raise RROMPyException("Radial basis not recognized.")
Van = np.zeros((nx, nx))
if reorder is not None: x = x[reorder]
xDiff2V, xDiff2I, xDiff2J = np.zeros(0), [], []
for i in range(nx - 1):
xiD2Loc = np.sum(np.abs((x[i + 1 :] - x[i])
* directionalWeights) ** 2., axis = 1)
xiD2Good = np.where(xiD2Loc <= radialker.threshold)[0]
xDiff2V = np.append(xDiff2V, xiD2Loc[xiD2Good])
xDiff2I += [i] * len(xiD2Good)
xDiff2J += list(i + 1 + xiD2Good)
kernelV = radialker(xDiff2V, apply_threshold = False)
Van = np.eye(nx)
Van[xDiff2I, xDiff2J] = kernelV
Van[xDiff2J, xDiff2I] = kernelV
return Van
def polyvander(x:paramList, degs:List[int], basis:str,
derIdxs : List[List[List[int]]] = None,
reorder : List[int] = None, directionalWeights : Np1D = None,
scl : Np1D = None, forceTotalDegree : bool = False,
optimizeScalingBounds : List[float] = [-1., -1.],
maxOptimizeIter : int = 10) -> Tuple[Np2D, Np1D]:
"""
Compute full Hermite-Vandermonde matrix with specified derivative
directions.
"""
if derIdxs is not None and np.sum(np.sum(derIdxs)) > 0:
raise RROMPyException(("Cannot take derivatives of radial basis "
"function."))
basisp, basisr = basis.split("_")
VanP = pvP(x, degs, basisp, derIdxs = derIdxs, reorder = reorder,
scl = scl, forceTotalDegree = forceTotalDegree)
VanZ = np.zeros([VanP.shape[1]] * 2)
optDir, optFact = 0, 100.
for it in range(maxOptimizeIter):
VanR = rbvander(x, basisr, reorder = reorder,
directionalWeights = directionalWeights)
Van = np.block([[VanR, VanP], [VanP.T.conj(), VanZ]])
if optimizeScalingBounds[0] < 0. or it == maxOptimizeIter - 1: break
ncond = np.linalg.cond(Van)
if ncond < optimizeScalingBounds[0]:
if optDir != -1: optFact **= .5
optDir, directionalWeights = -1, directionalWeights / optFact
elif ncond > optimizeScalingBounds[1]:
if optDir != 1: optFact **= .5
optDir, directionalWeights = 1, directionalWeights * optFact
else: break
return Van, directionalWeights

Event Timeline

c4science · Help