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 9f62f2c..f130eec 100644
--- a/rrompy/parameter/parameter_sampling/sparse_grid/sparse_grid_sampler.py
+++ b/rrompy/parameter/parameter_sampling/sparse_grid/sparse_grid_sampler.py
@@ -1,116 +1,112 @@
 # 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.poly_fitting.piecewise_linear import (sparsekinds,
                                                             sparseMap)
 from rrompy.utilities.base.types import List, paramList
 from rrompy.utilities.numerical import lowDiscrepancy
 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 sparsekinds:
+        if kind.upper() not in [sk.split("_")[2] for sk in sparsekinds]:
             raise RROMPyException("Generator kind not recognized.")
         self._kind = kind.upper()
 
     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:
                             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[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 : bool = True) -> paramList:
+    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)
-        x = self.points
-        if self.npoints > 1 and reorder:
-            fejerOrdering = lowDiscrepancy(self.npoints)
-            x = checkParameterList(x.data[fejerOrdering, :], self.npar)[0]
-        return x
+        return self.points
diff --git a/rrompy/utilities/poly_fitting/piecewise_linear/base.py b/rrompy/utilities/poly_fitting/piecewise_linear/base.py
index fc55c42..112234a 100644
--- a/rrompy/utilities/poly_fitting/piecewise_linear/base.py
+++ b/rrompy/utilities/poly_fitting/piecewise_linear/base.py
@@ -1,41 +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 <http://www.gnu.org/licenses/>.
 #
 
 import numpy as np
 from rrompy.utilities.base.types import Np1D, Tuple
 from rrompy.utilities.exception_manager import RROMPyException
 
 __all__ = ['sparsekinds', 'sparseMap']
 
-sparsekinds = ["UNIFORM", "LOBATTO"]
+sparsekinds = ["PIECEWISE_LINEAR_" + k for k in ["UNIFORM", "LOBATTO"]]
 
 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(" ", "")
+    kind = kind.upper().strip().replace(" ", "").split("_")[-1]
     if kind == "UNIFORM":
         if forward: return .5 * ((lims[-1] - lims[0]) * x + lims[0] + lims[-1])
         return np.real((2. * x - lims[0] - lims[-1]) / (lims[-1] - lims[0]))
     elif kind == "LOBATTO":
         if forward:
             x0 = np.cos(.5 * np.pi * (1. - x))
             return sparseMap(x0, lims, "UNIFORM", True)
         x0 = sparseMap(x, lims, "UNIFORM", False)
         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 eb643d8..1942660 100644
--- a/rrompy/utilities/poly_fitting/piecewise_linear/kernel.py
+++ b/rrompy/utilities/poly_fitting/piecewise_linear/kernel.py
@@ -1,63 +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 <http://www.gnu.org/licenses/>.
 #
 
 import numpy as np
 from .base import 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:
         width = .5 ** (depth - 1.)
         if np.abs(supp) + width > 1. + .25 * width:
             raise RROMPyException(("Wrong depth or support point. Support "
                                    "point too close to boundary."))
         y = np.array(1. - np.abs(x - supp) / width)
     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 = sparseMap(x(d), lims(d), kind, False)
         for j in range(len(supportPoints)):
             supp = sparseMap(supportPoints(j, d), lims(d), kind, 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(van, c, 1)
+    return np.tensordot(c, van, (0, -1))
diff --git a/rrompy/utilities/poly_fitting/piecewise_linear/piecewise_linear_interpolator.py b/rrompy/utilities/poly_fitting/piecewise_linear/piecewise_linear_interpolator.py
index 25a220c..164b19f 100644
--- a/rrompy/utilities/poly_fitting/piecewise_linear/piecewise_linear_interpolator.py
+++ b/rrompy/utilities/poly_fitting/piecewise_linear/piecewise_linear_interpolator.py
@@ -1,101 +1,101 @@
 # 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.utilities.base.types import (List, ListAny, DictAny, Np1D, Np2D,
                                          paramList)
 from rrompy.utilities.poly_fitting.interpolator import GenericInterpolator
 from rrompy.utilities.poly_fitting.custom_fit import customFit
 from .kernel import vander, val
 from rrompy.utilities.numerical import dot
 from rrompy.utilities.exception_manager import RROMPyException, RROMPyAssert
 from rrompy.parameter import checkParameterList
 
 __all__ = ['PiecewiseLinearInterpolator']
 
 class PiecewiseLinearInterpolator(GenericInterpolator):
     def __init__(self, other = None):
         if other is None: return
         self.support = other.support
         self.lims = other.lims
         self.coeffs = other.coeffs
         self.depths = other.depths
         self.npar = other.npar
         self.kind = other.kind
 
     @property
     def shape(self):
         sh = self.coeffs.shape[1 :] if self.coeffs.ndim > 1 else 1
         return sh
 
     def __call__(self, mu:paramList, der : List[int] = None,
                  scl : Np1D = None):
         if der is not None and np.sum(der) > 0:
             raise RROMPyException(("Cannot take derivatives of piecewise "
                                    "linear function."))           
         return val(mu, self.coeffs, self.support, self.depths, self.kind,
                    self.lims)
         
     def __copy__(self):
-        return RadialBasisInterpolator(self)
+        return PiecewiseLinearInterpolator(self)
 
     def __deepcopy__(self, memo):
-        other = RadialBasisInterpolator()
+        other = PiecewiseLinearInterpolator()
         (other.support, other.lims, other.coeffs, other.depths, other.npar,
          other.kind) = copy((self.support, self.lims, self.coeffs, self.depths,
                              self.npar, self.kind), memo)
         return other
 
     def postmultiplyTensorize(self, A:Np2D):
         RROMPyAssert(A.shape[0], self.shape[-1], "Shape of output")
         self.coeffs = dot(self.coeffs, A)
 
     def pad(self, nleft : List[int] = None, nright : List[int] = None):
         if nleft is None: nleft = [0] * len(self.shape)
         if nright is None: nright = [0] * len(self.shape)
         if not hasattr(nleft, "__len__"): nleft = [nleft]
         if not hasattr(nright, "__len__"): nright = [nright]
         RROMPyAssert(len(self.shape), len(nleft), "Shape of output")
         RROMPyAssert(len(self.shape), len(nright), "Shape of output")
         padwidth = [(0, 0)] + [(l, r) for l, r in zip(nleft, nright)]
         self.coeffs = np.pad(self.coeffs, padwidth, "constant",
                              constant_values = (0., 0.))
         padwidth = [(0, 0)] * (self.npar - 1) + padwidth
 
     def setupByInterpolation(self, support:paramList, values:ListAny,
                              lims:paramList, depths:Np2D,
-                             kind : str = "UNIFORM", verbose : bool = True,
-                             fitCoeffs : DictAny = {}):
+                             kind : str = "PIECEWISE_LINEAR_UNIFORM",
+                             verbose : bool = True, fitCoeffs : DictAny = {}):
         support = checkParameterList(support)[0]
         self.support = copy(support)
         self.npar = support.shape[1]
         lims = checkParameterList(lims, self.npar)[0]
         self.lims = copy(lims)
         self.depths = copy(depths)
         self.kind = kind
         values = values.reshape(values.shape[0], -1)
         van = vander(support, depths, kind, lims)
         fitOut = customFit(van, values, full = True, **fitCoeffs)
         if verbose:
             msg = ("Fitting {} samples through piecewise linear "
                    "interpolator... Conditioning of system: {:.4e}.").format(
                               len(support), fitOut[1][2][0] / fitOut[1][2][-1])
         else: msg = None
         self.coeffs = fitOut[0]
         return fitOut[1][1] == van.shape[1], msg