Page MenuHomec4science
Files F61376236

quadrature_circle_sampler.py
No OneTemporary
Actions

File Metadata

Created
Mon, May 6, 07:05

quadrature_circle_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 .generic_shape_quadrature_sampler import GenericShapeQuadratureSampler
from rrompy.utilities.base.types import paramList
from rrompy.utilities.numerical import (lowDiscrepancy, kroneckerer,
potential, quadraturePointsGenerate)
__all__ = ['QuadratureCircleSampler']
class QuadratureCircleSampler(GenericShapeQuadratureSampler):
"""Generator of quadrature sample points on ellipses."""
def generatePoints(self, n:int, reorder : bool = True) -> paramList:
"""Array of sample points."""
n1d = int(np.ceil(n ** (1. / self.npar)))
xds, nds = [None] * self.npar, [None] * self.npar
limsE = self.mapParameterList(self.lims)
for d in range(self.npar):
ax = self.axisRatios[d]
scorebase = self.groundPotential(d)
a, b = limsE(d)
c, r = (a + b) / 2., (a - b) / 2.
n1dx = int(np.ceil((n1d * 4. / np.pi / ax) ** .5))
n1dy = int(np.ceil(n1d * 4. / np.pi / n1dx))
even, Z = True, []
while len(Z) < n1d:
Xdx, Xdy = np.meshgrid(
quadraturePointsGenerate(n1dx, self.kind),
quadraturePointsGenerate(n1dy, self.kind))
Z = Xdx.flatten() + 1.j * ax * Xdy.flatten()
ptscore = potential(Z, self.normalFoci(d))
Z = Z[ptscore <= scorebase]
if even: n1dx += 1
else: n1dy += 1
xds[d] = c + r * Z
nds[d] = len(Z)
nleft, nright = 1, np.prod(nds)
xmat = np.empty((nright, self.npar), dtype = np.complex)
for d in range(self.npar):
nright //= nds[d]
xmat[:, d] = kroneckerer(xds[d], nleft, nright)
nleft *= nds[d]
if len(xmat) > 1 and reorder:
fejerOrdering = [len(xmat) - 1] + lowDiscrepancy(len(xmat) - 1)
xmat = xmat[fejerOrdering, :]
return self.mapParameterList(xmat, "B")

Event Timeline

c4science · Help