Homec4science
Diffusion RationalROMPy (master) Browse

RROMPy/rrompy/sampling/quadrature_sampler.py
c6af74908c2cmaster

Last Change
Raw File

quadrature_sampler.py
View Options

#!/usr/bin/python
# 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 rrompy.sampling.generic_sampler import GenericSampler
from rrompy.utilities.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):
super().__init__(self, lims = lims)
self.kind = kind
@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) -> List[Np1D]:
"""Array of quadrature points and array of weights."""
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.kind == "UNIFORM":
xj = np.linspace(a, b, n[j])[:, None]
wj = np.abs(a - b) / n[j] * 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 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
c4science · Help