diff --git a/rrompy/parameter/parameter_sampling/__init__.py b/rrompy/parameter/parameter_sampling/__init__.py index 2ed053b..b44b6da 100644 --- a/rrompy/parameter/parameter_sampling/__init__.py +++ b/rrompy/parameter/parameter_sampling/__init__.py @@ -1,40 +1,41 @@ # 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 . # from .empty_sampler import EmptySampler from .manual_sampler import ManualSampler from .segment import QuadratureSampler, QuadratureSamplerTotal, RandomSampler from .shape import (FFTSampler, QuadratureBoxSampler, QuadratureCircleSampler, RandomBoxSampler, RandomCircleSampler) -from .sparse_grid import SparseGridSampler +from .sparse_grid import SparseGridSampler, SparseGridSamplerTwoWay __all__ = [ 'EmptySampler', 'ManualSampler', 'QuadratureSampler', 'QuadratureSamplerTotal', 'RandomSampler', 'FFTSampler', 'QuadratureBoxSampler', 'QuadratureCircleSampler', 'RandomBoxSampler', 'RandomCircleSampler', - 'SparseGridSampler' + 'SparseGridSampler', + 'SparseGridSamplerTwoWay' ] diff --git a/rrompy/parameter/parameter_sampling/sparse_grid/__init__.py b/rrompy/parameter/parameter_sampling/sparse_grid/__init__.py index 64373b4..f69816e 100644 --- a/rrompy/parameter/parameter_sampling/sparse_grid/__init__.py +++ b/rrompy/parameter/parameter_sampling/sparse_grid/__init__.py @@ -1,25 +1,27 @@ # 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 . # from .sparse_grid_sampler import SparseGridSampler +from .sparse_grid_sampler_two_way import SparseGridSamplerTwoWay __all__ = [ - 'SparseGridSampler' + 'SparseGridSampler', + 'SparseGridSamplerTwoWay' ] diff --git a/rrompy/parameter/parameter_sampling/sparse_grid/sparse_grid_sampler.py b/rrompy/parameter/parameter_sampling/sparse_grid/sparse_grid_sampler.py index 5293190..419ba12 100644 --- a/rrompy/parameter/parameter_sampling/sparse_grid/sparse_grid_sampler.py +++ b/rrompy/parameter/parameter_sampling/sparse_grid/sparse_grid_sampler.py @@ -1,110 +1,105 @@ # 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 . # from copy import deepcopy as copy from itertools import product import numpy as np from rrompy.parameter.parameter_sampling.generic_sampler import GenericSampler from rrompy.utilities.poly_fitting.piecewise_linear import (sparsekinds, sparseMap) -from rrompy.utilities.base.types import Tuple, List, DictAny, paramList +from rrompy.utilities.base.types import Tuple, List, Np1D, DictAny, paramList from rrompy.utilities.exception_manager import RROMPyException __all__ = ['SparseGridSampler'] class SparseGridSampler(GenericSampler): """Generator of sparse grid sample points.""" def __init__(self, lims:paramList, kind : str = "UNIFORM", parameterMap : DictAny = 1.): super().__init__(lims = lims, parameterMap = parameterMap) 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 [sk.split("_")[2] + extra for sk, extra in - product(sparsekinds, ["", "-NOBOUNDARY"])]: + product(sparsekinds, ["", "-HAAR"])]: raise RROMPyException("Generator kind not recognized.") self._kind = kind.upper() - self._noBoundary = "NOBOUNDARY" in self._kind + self._noBoundary = "HAAR" in self._kind def reset(self): limsE = self.mapParameterList(self.lims) centerEff = .5 * (limsE[0] + limsE[1]) self.points = self.mapParameterList(centerEff, "B") - self.depth = np.zeros((1, self.npar), dtype = int) + self.depth = np.array([[self._noBoundary] * self.npar], dtype = int) def refine(self, active : List[int] = None) -> Tuple[List[int], List[int]]: if active is None: active = np.arange(self.npoints) active = np.array(active) if np.any(active < 0) or np.any(active >= self.npoints): raise RROMPyException(("Active indices must be between 0 " "(included) and npoints (excluded).")) - limsX = self.mapParameterList(self.lims) newIdxs, oldIdxs = [], [] for act in active: point, dpt = self.points[act], self.depth[act] for jdelta, signdelta in product(range(self.npar), [-1., 1.]): - Pointj = sparseMap(self.mapParameterList(point[jdelta], - idx = [jdelta])(0, 0), - limsX(jdelta), self.kind, False) - if not self._noBoundary and dpt[jdelta] == 1: - Centerj = sparseMap( - self.mapParameterList(self.points[0][jdelta], - idx = [jdelta])(0, 0), - limsX(jdelta), self.kind, False) - gradj = Pointj - Centerj - if signdelta * gradj > 0: - continue - pointj = copy(point) - Pointj = Pointj + .5 ** (dpt[jdelta] - + self._noBoundary) * signdelta - pointj[jdelta] = self.mapParameterList( - sparseMap(Pointj, limsX(jdelta), self.kind), - "B", [jdelta])(0, 0) - dist = np.sum(np.abs(self.points.data - - pointj.reshape(1, -1)), axis = 1) - samePt = np.where(np.isclose(dist, 0.))[0] - if len(samePt) > 0: - oldIdxs += [samePt[0]] - continue - newIdxs += [self.npoints] - self.points.append(pointj) - self.depth = np.append(self.depth, [dpt], 0) - self.depth[-1, jdelta] += 1 + idx = self.addForwardPoint(point, dpt, jdelta, signdelta) + if idx is not None: + if idx > 0: newIdxs += [idx] + else: oldIdxs += [- idx] return newIdxs, oldIdxs + + def addForwardPoint(self, basepoint:Np1D, basedepth:Np1D, index:int, + sign:float) -> int: + if basedepth[index] < self._noBoundary: return None + limd = self.mapParameterList(self.lims(index), idx = [index])(0) + xd0 = sparseMap(self.mapParameterList(basepoint[index], + idx = [index])(0, 0), + limd, self.kind, False) + .5 ** basedepth[index] * sign + if np.abs(xd0) >= 1. + 1e-15 * (1 - 2 * self._noBoundary): return None + pt = copy(basepoint) + pt[index] = self.mapParameterList(sparseMap(xd0, limd, self.kind), + "B", [index])(0, 0) + dist = np.sum(np.abs(self.points.data - pt.reshape(1, -1)), axis = 1) + samePt = np.where(np.isclose(dist, 0.))[0] + if len(samePt) > 0: return - samePt[0] + self.points.append(pt) + self.depth = np.append(self.depth, [basedepth], 0) + self.depth[-1, index] += 1 + return self.npoints - 1 def generatePoints(self, n:int, reorder = None) -> paramList: if self.npoints > n: self.reset() idx = np.arange(self.npoints) while self.npoints < n: idx = self.refine(idx)[0] return self.points diff --git a/rrompy/parameter/parameter_sampling/sparse_grid/sparse_grid_sampler_two_way.py b/rrompy/parameter/parameter_sampling/sparse_grid/sparse_grid_sampler_two_way.py new file mode 100644 index 0000000..12a379b --- /dev/null +++ b/rrompy/parameter/parameter_sampling/sparse_grid/sparse_grid_sampler_two_way.py @@ -0,0 +1,47 @@ +# 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 . +# + +from itertools import product +import numpy as np +from .sparse_grid_sampler import SparseGridSampler +from rrompy.utilities.base.types import Tuple, List +from rrompy.utilities.exception_manager import RROMPyException + +__all__ = ['SparseGridSamplerTwoWay'] + +class SparseGridSamplerTwoWay(SparseGridSampler): + """Generator of sparse grid sample points with two-way refinement.""" + + def refine(self, active : List[int] = None) -> Tuple[List[int], List[int]]: + if active is None: active = np.arange(self.npoints) + active = np.array(active) + if np.any(active < 0) or np.any(active >= self.npoints): + raise RROMPyException(("Active indices must be between 0 " + "(included) and npoints (excluded).")) + newIdxs, oldIdxs = [], [] + for act in active: + point, dpt = self.points[act], self.depth[act] + for jdelta, signdelta, backwards in product(range(self.npar), + [-1., 1.], [0, 1]): + if backwards: dpt[jdelta] -= 1 + idx = self.addForwardPoint(point, dpt, jdelta, signdelta) + if idx is not None: + if idx > 0: newIdxs += [idx] + else: oldIdxs += [- idx] + if backwards: dpt[jdelta] += 1 + return newIdxs, oldIdxs diff --git a/rrompy/utilities/poly_fitting/piecewise_linear/kernel.py b/rrompy/utilities/poly_fitting/piecewise_linear/kernel.py index 1b995a1..60c5f69 100644 --- a/rrompy/utilities/poly_fitting/piecewise_linear/kernel.py +++ b/rrompy/utilities/poly_fitting/piecewise_linear/kernel.py @@ -1,65 +1,63 @@ # 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 . # import numpy as np from .base import centerNormalize, sparseMap from rrompy.utilities.base.types import Np1D, Np2D, paramList from rrompy.parameter import checkParameterList __all__ = ['hatFunction', 'val', 'vander'] def hatFunctionRef(x:Np1D, supp:float, depth:int, kind:str) -> Np1D: - if depth < 0: y = np.zeros_like(x) - elif depth == 0: y = np.ones_like(x) + noBdr = "HAAR" in kind.upper() + if depth < noBdr: return np.zeros_like(x) + if depth == noBdr: return np.ones_like(x) + suppEff = sparseMap(supp, [-1., 1.], kind, False) + suppLREff = suppEff + .5 ** (depth - 1) * np.array([-1., 1.]) + widthL, widthR = sparseMap(suppLREff, [-1., 1.], kind) - supp + xC = np.array(x - supp) + if np.isclose(widthL, 0.) or supp + (.95 + .1 * noBdr) * widthL < - 1.: + isleft, isright = 0, 1 + elif np.isclose(widthR, 0.) or supp + (.95 + .1 * noBdr) * widthR > 1.: + isleft, isright = 1, 0 else: - noBdr = "NOBOUNDARY" in kind.upper() - suppEff = sparseMap(supp, [-1., 1.], kind, False) - exp = depth if noBdr else depth - 1 - suppLREff = suppEff + .5 ** exp * np.array([-1., 1.]) - widthL, widthR = sparseMap(suppLREff, [-1., 1.], kind) - supp - xC = np.array(x - supp) - if np.isclose(widthL, 0.) or supp + (.95 + .1 * noBdr) * widthL < - 1.: - isleft, isright = 0, 1 - elif np.isclose(widthR, 0.) or supp + (.95 + .1 * noBdr) * widthR > 1.: - isleft, isright = 1, 0 - else: - isleft, isright = xC < 0., xC >= 0. - y = 1. - xC / (widthL * isleft + widthR * isright) + isleft, isright = xC < 0., xC >= 0. + y = 1. - xC / (widthL * isleft + widthR * isright) return np.clip(y, 0., 1., y) def hatFunction(x:paramList, supportPoints:paramList, depths:Np2D, kind:str, lims:paramList) -> Np2D: x = checkParameterList(x) supportPoints = checkParameterList(supportPoints, x.shape[1]) lims = checkParameterList(lims, x.shape[1]) res = np.ones((len(supportPoints), len(x))) for d in range(x.shape[1]): x0 = centerNormalize(x(d), lims(d), False) for j in range(len(supportPoints)): supp = centerNormalize(supportPoints(j, d), lims(d), False) res[j] *= hatFunctionRef(x0, supp, depths[j, d], kind) return res.T def vander(supportPoints:paramList, depths:Np2D, kind:str, lims:paramList) -> Np2D: return hatFunction(supportPoints, supportPoints, depths, kind, lims) def val(x:paramList, c:Np2D, supportPoints:paramList, depths:Np2D, kind:str, lims:paramList) -> Np2D: van = hatFunction(x, supportPoints, depths, kind, lims) return np.tensordot(c, van, (0, -1))