Page MenuHomec4science

quadrature_sampler.py
No OneTemporary

File Metadata

Created
Thu, May 9, 21:44

quadrature_sampler.py

# 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 .generic_sampler import GenericSampler
from rrompy.utilities.base.types import List, paramList
from rrompy.utilities.exception_manager import RROMPyException
from rrompy.utilities.numerical import lowDiscrepancy, kroneckerer
from rrompy.parameter import checkParameterList
__all__ = ['QuadratureSampler']
class QuadratureSampler(GenericSampler):
"""Generator of quadrature sample points."""
allowedKinds = ["UNIFORM", "CHEBYSHEV", "GAUSSLEGENDRE"]
def __init__(self, lims:paramList, kind : str = "UNIFORM",
scalingExp : List[float] = None):
super().__init__(lims = lims, scalingExp = scalingExp)
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:List[int], reorder : bool = True) -> paramList:
"""Array of sample points."""
if not hasattr(n, "__len__"): n = [n]
super().generatePoints(n)
nleft, nright = 1, np.prod(n)
xmat = np.empty((nright, self.npar), dtype = self.lims.dtype)
for d in range(self.npar):
nright //= n[d]
a = self.lims(0, d) ** self.scalingExp[d]
b = self.lims(1, d) ** self.scalingExp[d]
c, r = (a + b) / 2., (a - b) / 2.
if self.kind == "UNIFORM":
xd = np.linspace(a, b, n[d])
elif self.kind == "CHEBYSHEV":
nodes, _ = np.polynomial.chebyshev.chebgauss(n[d])
xd = c + r * nodes
elif self.kind == "GAUSSLEGENDRE":
nodes, _ = np.polynomial.legendre.leggauss(n[d])
xd = c + r * nodes[::-1]
xd **= 1. / self.scalingExp[d]
if n[d] > 1 and reorder:
fejerOrdering = [n[d] - 1] + lowDiscrepancy(n[d] - 1)
xd = xd[fejerOrdering]
xmat[:, d] = kroneckerer(xd, nleft, nright)
nleft *= n[d]
x = checkParameterList(xmat, self.npar)[0]
return x

Event Timeline