Page MenuHomec4science
Files F60157359

fft_sampler.py
No OneTemporary
Actions

File Metadata

Created
Sat, Apr 27, 22:46

fft_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, List, Tuple
from rrompy.utilities.exception_manager import RROMPyException
from rrompy.utilities.base import lowDiscrepancy
__all__ = ['FFTSampler']
class FFTSampler(GenericSampler):
"""Generator of FFT-type sample points on scaled roots of unity."""
def generatePoints(self, n:Np1D) -> Tuple[List[Np1D], Np1D]:
"""Array of sample points and array of weights."""
super().generatePoints(n)
d = len(self.lims[0])
try:
len(n)
except:
n = np.array([n])
if len(n) != d:
raise RROMPyException(("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.scaling is not None:
a, b = self.scaling[j](a), self.scaling[j](b)
c, r = (a + b) / 2., np.abs(a - b) / 2.
xj = c + r * np.exp(1.j *
np.linspace(0, 2 * np.pi, n[j] + 1)[:-1, None])
wj = r / n[j] * np.ones(n[j])
fejerOrdering = lowDiscrepancy(len(xj))
xj = xj[fejerOrdering]
wj = wj[fejerOrdering]
if self.scalingInv is not None:
xj = self.scalingInv[j](xj)
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
def refine(self, n:int) -> Tuple[List[Np1D], Np1D]:
"""
Apply refinement. If points are not nested, equivalent to
generatePoints([2 * x for x in n]).
"""
super().generatePoints(n)
d = len(self.lims[0])
try:
len(n)
except:
n = np.array([n])
if len(n) != d:
raise RROMPyException(("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.scaling is not None:
a, b = self.scaling[j](a), self.scaling[j](b)
c, r = (a + b) / 2., np.abs(a - b) / 2.
xj = c + r * np.exp(1.j * (np.pi / n[j]
+ np.linspace(0, 2 * np.pi, n[j] + 1)[:-1, None]))
wj = r / n[j] * np.ones(n[j])
fejerOrdering = lowDiscrepancy(len(xj))
xj = xj[fejerOrdering]
wj = wj[fejerOrdering]
if self.scalingInv is not None:
xj = self.scalingInv[j](xj)
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

Event Timeline

c4science · Help