Page MenuHomec4science
Files F61039899

quadrature_sampler_total.py
No OneTemporary
Actions

File Metadata

Created
Sat, May 4, 03:53

quadrature_sampler_total.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 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."""
if hasattr(n, "__len__"): n = n[0]
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)
idxBase = [x + 1 for x in lowDiscrepancy(n1d - 1, inverse = True)]
linearIdxs = np.array(idxBase + [0])
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.data[keepIdxs, :]
if reorder:
fejerTot = np.array([nTot - 1] + list(lowDiscrepancy(nTot - 1)))
xmat = xmat[np.argsort(np.argsort(fejerTot[keepIdxs])), :]
x = checkParameterList(xmat, d)[0]
return x

Event Timeline

c4science · Help