Page MenuHomec4science

sparse_grid_sampler.py
No OneTemporary

File Metadata

Created
Sun, May 12, 21:15

sparse_grid_sampler.py

# 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 copy import deepcopy as copy
from rrompy.parameter import checkParameterList
from rrompy.parameter.parameter_sampling.generic_sampler import GenericSampler
from rrompy.utilities.base.types import List, paramList
from rrompy.utilities.numerical import lowDiscrepancy
from rrompy.utilities.exception_manager import RROMPyException
_allowedSparseGridKinds = ["UNIFORM", "LOBATTO"]
__all__ = ['SparseGridSampler']
class SparseGridSampler(GenericSampler):
"""Generator of sparse grid sample points."""
def __init__(self, lims:paramList, kind : str = "UNIFORM",
scalingExp : List[float] = None):
super().__init__(lims = lims, scalingExp = scalingExp)
self._allowedKinds = _allowedSparseGridKinds
self.kind = kind
self.reset()
def __str__(self) -> str:
return "{}[{}_{}]_{}".format(self.name(), self.lims[0],
self.lims[1], self.kind)
@property
def npoints(self):
"""Number of points."""
return len(self.points)
@property
def kind(self):
"""Value of kind."""
return self._kind
@kind.setter
def kind(self, kind):
if kind.upper() not in self._allowedKinds:
raise RROMPyException("Generator kind not recognized.")
self._kind = kind.upper()
if self.kind == "UNIFORM":
self._phiD = lambda x: x
self._phiI = lambda y: y
elif self.kind == "LOBATTO":
self._phiD = lambda x: np.cos(.5 * np.pi * (x + 1.))
self._phiI = lambda y: 2. / np.pi * np.arccos(np.clip(y, -1., 1.,
y)) - 1.
def reset(self):
self.deltalims = .5 * np.abs(self.lims.data[1] ** self.scalingExp
- self.lims.data[0] ** self.scalingExp)
self.points = checkParameterList(self._phiD(.5 * (
self.lims[0] ** self.scalingExp
+ self.lims[1] ** self.scalingExp)
) ** (1. / self.scalingExp), self.npar)[0]
self.depth = np.zeros((1, self.npar))
self.deltadepth = np.zeros(1)
def refine(self, active : List[int] = None) -> List[int]:
if active is None: active = np.arange(self.npoints)
newIdxs = []
for act in active:
pnt, dpt = self.points[act], self.depth[act]
exhausted = False
while not exhausted:
ddp = self.deltadepth[act]
for jdelta in range(self.npar):
for signdelta in [-1., 1.]:
if np.isclose(dpt[jdelta], 1.):
gradj = (
self._phiI(pnt[jdelta] ** self.scalingExp[jdelta])
- self._phiI(self.points[0, jdelta] ** self.scalingExp[jdelta]))
if signdelta * gradj > 0:
continue
pntNew = copy(pnt)
pntNew[jdelta] = self._phiD(
self._phiI(pntNew[jdelta] ** self.scalingExp[jdelta])
+ .5 ** (dpt[jdelta] + ddp) * self.deltalims[jdelta]
* signdelta) ** (1. / self.scalingExp[jdelta])
dist = np.sum(np.abs(self.points.data
- pntNew.reshape(1, -1)), axis = 1)
samePt = np.where(np.isclose(dist, 0.))[0]
if len(samePt) > 0:
if samePt[0] in newIdxs: exhausted = True
continue
newIdxs += [self.npoints]
self.points.append(pntNew)
self.depth = np.append(self.depth, [dpt], 0)
self.depth[-1, jdelta] += 1 + ddp
self.deltadepth = np.append(self.deltadepth, [0])
exhausted = True
self.deltadepth[act] += 1
return newIdxs
def generatePoints(self, n:int, reorder : bool = True) -> paramList:
if self.npoints > n: self.reset()
idx = np.arange(self.npoints)
while self.npoints < n: idx = self.refine(idx)
x = self.points
if self.npoints > 1 and reorder:
fejerOrdering = lowDiscrepancy(self.npoints)
x = checkParameterList(x.data[fejerOrdering, :], self.npar)[0]
return x

Event Timeline