Homec4science
Diffusion RationalROMPy (develop) Browse

RROMPy/rrompy/parameter/parameter_sampling/shape/random_circle_sampler.py
d6451a36626adevelop

Last Change
Raw File

random_circle_sampler.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/>.
#
from numbers import Number
import numpy as np
from .generic_shape_random_sampler import GenericShapeRandomSampler
from rrompy.utilities.numerical import (haltonGenerate, sobolGenerate,
potential)
from rrompy.utilities.base.types import Tuple, List
__all__ = ['RandomCircleSampler']
class RandomCircleSampler(GenericShapeRandomSampler):
"""Generator of (quasi-)random sample points on ellipses."""
def refine(self, active : List[int] = None) -> Tuple[List[int], List[int]]:
if active is None:
n = self.npoints
elif isinstance(active, (Number,)):
n = active
else:
n = len(active)
n = int(n * self.refinementFactor)
nEff = int(np.ceil(n * (4. / np.pi) ** self.npar * np.prod(
[max(x, 1. / x) for x in self.axisRatios])))
xmat2 = []
while len(xmat2) < n:
if self.kind == "UNIFORM":
xmat2 = np.random.uniform(size = (nEff, 2 * self.npar))
elif self.kind == "HALTON":
xmat2, self.seedLoc = haltonGenerate(2 * self.npar, nEff,
self.seedLoc,
return_seed = True)
else:
xmat2, self.seedLoc = sobolGenerate(2 * self.npar, nEff,
self.seed,
return_seed = True)
xmat2 = xmat2 * 2. - 1.
for d in range(self.npar):
ax = self.axisRatios[d]
if ax > 1.: xmat2[:, 2 * d : 2 * d + 2] *= ax
Z = xmat2[:, 2 * d] + 1.j * ax * xmat2[:, 2 * d + 1]
xmat2 = xmat2[potential(Z, self.normalFoci(d)) <= 1.]
nEff += 1
xmat = np.empty((n, self.npar), dtype = np.complex)
limsE = self.mapParameterList(self.lims)
for d in range(self.npar):
ax = self.axisRatios[d]
a, b = limsE(d)
c, r = (a + b) / 2., (a - b) / 2.
xmat[:, d] = c + r * (xmat2[: n, 2 * d]
+ 1.j * ax * xmat2[: n, 2 * d + 1])
pts = self.mapParameterList(xmat, "B")
idx = np.arange(n, dtype = int) + len(self.points)
for pj in pts: self.points.append(pj)
return list(idx), []
c4science · Help