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 6e9fda0..c50aaf7 100644 --- a/rrompy/parameter/parameter_sampling/sparse_grid/sparse_grid_sampler.py +++ b/rrompy/parameter/parameter_sampling/sparse_grid/sparse_grid_sampler.py @@ -1,111 +1,115 @@ # 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 copy import deepcopy as copy +from itertools import product +import numpy as np from rrompy.parameter import checkParameterList 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 List, 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", scalingExp : List[float] = None): super().__init__(lims = lims, scalingExp = scalingExp) 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] for sk in sparsekinds]: + if kind.upper() not in [sk.split("_")[2] + extra for sk, extra in + product(sparsekinds, ["", "-NOBOUNDARY"])]: raise RROMPyException("Generator kind not recognized.") self._kind = kind.upper() + self._noBoundary = "NOBOUNDARY" in self._kind def reset(self): centerEff = .5 * (self.lims[0] ** self.scalingExp + self.lims[1] ** self.scalingExp) self.points = checkParameterList(centerEff ** (1. / self.scalingExp), self.npar)[0] self.depth = np.zeros((1, self.npar), dtype = int) self.deltadepth = np.zeros(1, dtype = int) def refine(self, active : List[int] = None) -> List[int]: if active is None: active = np.arange(self.npoints) limsX = [self.lims(j) ** self.scalingExp[j] for j in range(self.npar)] newIdxs = [] for act in active: point, 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.]: Pointj = sparseMap( point[jdelta] ** self.scalingExp[jdelta], limsX[jdelta], self.kind, False) - if dpt[jdelta] == 1: + if not self._noBoundary and dpt[jdelta] == 1: Centerj = sparseMap( self.points(0, jdelta) ** self.scalingExp[jdelta], limsX[jdelta], self.kind, False) gradj = Pointj - Centerj if signdelta * gradj > 0: continue pointj = copy(point) - Pointj = Pointj + .5 ** (dpt[jdelta] + ddp) * signdelta + Pointj = Pointj + .5 ** (dpt[jdelta] + ddp + + self._noBoundary) * signdelta pointj[jdelta] = sparseMap(Pointj, limsX[jdelta], self.kind) ** (1. / self.scalingExp[jdelta]) 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: if samePt[0] in newIdxs: exhausted = True continue newIdxs += [self.npoints] self.points.append(pointj) 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 = None) -> paramList: if self.npoints > n: self.reset() idx = np.arange(self.npoints) while self.npoints < n: idx = self.refine(idx) return self.points diff --git a/rrompy/utilities/poly_fitting/piecewise_linear/base.py b/rrompy/utilities/poly_fitting/piecewise_linear/base.py index 1471e7c..eb9dd52 100644 --- a/rrompy/utilities/poly_fitting/piecewise_linear/base.py +++ b/rrompy/utilities/poly_fitting/piecewise_linear/base.py @@ -1,47 +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 . # import numpy as np from rrompy.utilities.base.types import Np1D, Tuple from rrompy.utilities.exception_manager import RROMPyException __all__ = ['sparsekinds', 'sparseMap'] sparsekinds = ["PIECEWISE_LINEAR_" + k for k in ["UNIFORM", "CLENSHAWCURTIS"]] def centerNormalize(x:Np1D, lims:Tuple[np.complex, np.complex], forward : bool = True) -> Np1D: """forward: X([-1, 1]) -> X(lims)""" center, width = .5 * (lims[0] + lims[-1]), .5 * (lims[-1] - lims[0]) if forward: return width * x + center return np.real((x - center) / width) def sparseMap(x:Np1D, lims:Tuple[np.complex, np.complex], kind:str, forward : bool = True) -> Np1D: """forward: U([-1, 1]) -> lims""" - kind = kind.upper().strip().replace(" ", "").split("_")[-1] + kind = kind.upper().strip().replace(" ", "").split("_")[-1].split("-")[0] if kind == "UNIFORM": return centerNormalize(x, lims, forward) elif kind == "CLENSHAWCURTIS": if forward: x0 = np.cos(.5 * np.pi * (1. - x)) return centerNormalize(x0, lims, forward) x0 = centerNormalize(x, lims, forward) return 1. - 2. / np.pi * np.arccos(np.clip(x0, -1., 1.)) else: raise RROMPyException("Sparse map kind not recognized.") diff --git a/rrompy/utilities/poly_fitting/piecewise_linear/kernel.py b/rrompy/utilities/poly_fitting/piecewise_linear/kernel.py index ce50937..8bc6a42 100644 --- a/rrompy/utilities/poly_fitting/piecewise_linear/kernel.py +++ b/rrompy/utilities/poly_fitting/piecewise_linear/kernel.py @@ -1,66 +1,65 @@ # 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.utilities.exception_manager import RROMPyException 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) - elif depth == 1: - if not np.isclose(np.abs(supp), 1.): - raise RROMPyException(("Wrong depth. Depth 1 points must be on " - "boundary.")) - y = np.array(np.sign(supp) * x) else: suppEff = sparseMap(supp, [-1., 1.], kind, False) - suppLREff = suppEff + .5 ** (depth - 1.) * np.array([-1., 1.]) + exp = depth if "NOBOUNDARY" in kind.upper() else depth - 1 + suppLREff = suppEff + .5 ** exp * np.array([-1., 1.]) widthL, widthR = sparseMap(suppLREff, [-1., 1.], kind) - supp - if supp + widthL < -1. + .25*widthL or supp + widthR > 1. + .25*widthR: - raise RROMPyException(("Wrong depth or support point. Support " - "point too close to boundary.")) - x = np.array(x) - y = 1. - (x - supp) / (widthL * (x < supp) + widthR * (x >= supp)) + xC = np.array(x - supp) + if np.isclose(widthL, 0.) or supp + widthL < - 1. - .1 * widthL: + isleft, isright = 0, 1 + elif np.isclose(widthR, 0.) or supp + widthR > 1. - .1 * widthR: + isleft, isright = 1, 0 + else: + 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)[0] supportPoints = checkParameterList(supportPoints, x.shape[1])[0] lims = checkParameterList(lims, x.shape[1])[0] 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))