Page MenuHomec4science
Files F61004621

quadrature_sampler_total.py
No OneTemporary
Actions

File Metadata

Created
Fri, May 3, 22:34

quadrature_sampler_total.py
View Options

# Copyright (C) 2018-2020 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, factorial as fact
from .quadrature_sampler import QuadratureSampler
from rrompy.utilities.base.types import paramList
from rrompy.utilities.numerical import lowDiscrepancy
from rrompy.parameter import checkParameterList
__all__ = ['QuadratureSamplerTotal']
class QuadratureSamplerTotal(QuadratureSampler):
"""
Generator of quadrature sample points for total degree polynomial
computations.
"""
def generatePoints(self, n:int, reorder : bool = True) -> paramList:
"""Array of sample points."""
d = self.npar
n1d = int((fact(d) * n) ** (1. / d))
while binom(n1d + d - 1, d) > n: n1d -= 1
x = super().generatePoints(n1d ** d, reorder = False)
nTot = n1d ** d
indicesBase = np.zeros(nTot, dtype = int)
if reorder:
idxBase = ([y + 1 for y in lowDiscrepancy(n1d - 1, inverse = True)]
+ [0])
else:
idxBase = list(range(n1d))
linearIdxs = np.array(idxBase)
nleft, nright = 1, nTot
for j in range(d):
nright //= n1d
kronIn = np.repeat(linearIdxs, nright)
indicesBase += np.tile(kronIn, nleft)
nleft *= n1d
keepIdxs = np.zeros(nTot, dtype = bool)
keepIdxs[indicesBase < n1d] = True
xmat = x[keepIdxs]
if reorder:
nx = len(xmat)
fejerOrdering = [nx - 1] + lowDiscrepancy(nx - 1, inverse = True)
xmat = xmat[fejerOrdering, :]
return checkParameterList(xmat, d)

Event Timeline

c4science · Help