Page MenuHomec4science
Files F61532470

quadrature_sampler.py
No OneTemporary
Actions

File Metadata

Created
Tue, May 7, 06:26

quadrature_sampler.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 rrompy.utilities.parameter_sampling.generic_sampler import GenericSampler
from rrompy.utilities.base.types import Np1D, Tuple
from rrompy.utilities.exception_manager import RROMPyException
from rrompy.utilities.base import lowDiscrepancy
__all__ = ['QuadratureSampler']
class QuadratureSampler(GenericSampler):
"""Generator of quadrature sample points."""
allowedKinds = ["UNIFORM", "CHEBYSHEV", "GAUSSLEGENDRE", "CLENSHAWCURTIS"]
def __init__(self, lims:Np1D, kind : str = "UNIFORM",
scaling : callable = None, scalingInv : callable = None):
super().__init__(lims = lims, scaling = scaling,
scalingInv = scalingInv)
self.kind = kind
def __str__(self) -> str:
return "{}_{}".format(super().__str__(), self.kind)
def __repr__(self) -> str:
return self.__str__() + " at " + hex(id(self))
@property
def kind(self):
"""Value of kind."""
return self._kind
@kind.setter
def kind(self, kind):
if kind.upper() not in self.allowedKinds:
raise RROMPyException("Generator kind not recognized.")
self._kind = kind.upper()
def generatePoints(self, n:int) -> Tuple[Np1D, Np1D]:
"""Array of quadrature points and array of weights."""
a, b = self.lims[0], self.lims[1]
if self.scaling is not None:
a, b = self.scaling(a), self.scaling(b)
if self.kind == "UNIFORM":
x = np.linspace(a, b, n)
w = np.abs(a - b) / n * np.ones(n)
elif self.kind == "CHEBYSHEV":
nodes, weights = np.polynomial.chebyshev.chebgauss(n)
x = (a + b) / 2 + (a - b) / 2 * nodes
w = np.abs(a - b) / np.pi * weights[:]
elif self.kind == "GAUSSLEGENDRE":
nodes, weights = np.polynomial.legendre.leggauss(n)
x = (a + b) / 2 + (a - b) / 2 * nodes[::-1]
w = np.abs(a - b) * weights[::-1]
elif self.kind == "CLENSHAWCURTIS":
thetas = np.pi / (n - 1) * np.arange(n)
nodes = np.cos(thetas)
weights = np.ones(n)
if n == 1:
weights[0] = 2.
else:
for j in range((n - 1) // 2):
bw = 1. + 1. * (2 * (j + 1) != n - 1)
weights -= (bw * np.cos(2. * (j + 1) * thetas)
/ (4. * j * (j + 2) + 3))
weights /= (n - 1)
weights[1 : -1] *= 2.
x = (a + b) / 2 + (a - b) / 2 * nodes
w = np.abs(a - b) / 2 * weights
if len(x) > 1:
fejerOrdering = [len(x) - 1] + lowDiscrepancy(len(x) - 1)
x = x[fejerOrdering]
w = w[fejerOrdering]
if self.scalingInv is not None:
x = self.scalingInv(x)
return x, w

Event Timeline

c4science · Help