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 f130eec..6e9fda0 100644
--- a/rrompy/parameter/parameter_sampling/sparse_grid/sparse_grid_sampler.py
+++ b/rrompy/parameter/parameter_sampling/sparse_grid/sparse_grid_sampler.py
@@ -1,112 +1,111 @@
 # 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 [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 = 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/reduction_methods/base/reduced_basis_utils.py b/rrompy/reduction_methods/base/reduced_basis_utils.py
index 2d1b8cd..cac7a3d 100644
--- a/rrompy/reduction_methods/base/reduced_basis_utils.py
+++ b/rrompy/reduction_methods/base/reduced_basis_utils.py
@@ -1,69 +1,66 @@
 # 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, Np2D, Tuple, List, ListAny,
-                                         sampList)
+from rrompy.utilities.base.types import Np1D, Np2D, Tuple, List, sampList
 from rrompy.utilities.numerical import dot
-from rrompy.utilities.numerical.hash_derivative import (
-                                                  hashIdxToDerivative as hashI)
 from rrompy.utilities.exception_manager import RROMPyAssert
 from rrompy.sampling import sampleList
 
 __all__ = ['projectAffineDecomposition']
 
 def projectAffineDecomposition(As:List[Np2D], bs:List[Np1D], pMat:sampList,
                                ARBsOld : List[Np2D] = None,
                                bRBsOld : List[Np1D] = None,
                                pMatOld : sampList = None)\
                                               -> Tuple[List[Np2D], List[Np1D]]:
     """Project affine decomposition of linear system onto basis."""
     RROMPyAssert((ARBsOld is None, bRBsOld is None),
                  (pMatOld is None, pMatOld is None),
                  "Old affine projected terms")
     if isinstance(pMat, (sampleList,)): pMat = pMat.data
     pMatC = pMat.conj()
     ARBs = [None] * len(As)
     bRBs = [None] * len(bs)
     if pMatOld is None:
         for j in range(max(len(As), len(bs))):
             if j < len(As):
                 ARBs[j] = dot(dot(As[j], pMat).T, pMatC).T
             if j < len(bs):
                 bRBs[j] = dot(pMatC.T, bs[j])
     else:
         RROMPyAssert((len(ARBsOld), len(bRBsOld)), (len(As), len(bs)),
                      "Old affine projected terms")
         if isinstance(pMatOld, (sampleList,)): pMatOld = pMatOld.data
         pMatOldC = pMatOld.conj()
         Sold = pMatOld.shape[1]
         Snew = pMat.shape[1]
         for j in range(max(len(As), len(bs))):
             if j < len(As):
                 ARBs[j] = np.empty((Sold + Snew, Sold + Snew),
                                    dtype = ARBsOld[j].dtype)
                 ARBs[j][: Sold, : Sold] = ARBsOld[j]
                 ARBs[j][: Sold, Sold :] = dot(dot(As[j], pMat).T, pMatOldC).T
                 ARBs[j][Sold :, : Sold] = dot(dot(As[j], pMatOld).T, pMatC).T
                 ARBs[j][Sold :, Sold :] = dot(dot(As[j], pMat).T, pMatC).T
             if j < len(bs):
                 bRBs[j] = np.empty((Sold + Snew), dtype = bRBsOld[j].dtype)
                 bRBs[j][: Sold] = bRBsOld[j]
                 bRBs[j][Sold :] = dot(pMatC.T, bs[j])
     return ARBs, bRBs
 
diff --git a/rrompy/utilities/numerical/compress_matrix.py b/rrompy/utilities/numerical/compress_matrix.py
index a68a285..76fe175 100644
--- a/rrompy/utilities/numerical/compress_matrix.py
+++ b/rrompy/utilities/numerical/compress_matrix.py
@@ -1,39 +1,38 @@
 # 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.numerical.tensor_la import dot
 from rrompy.utilities.base.types import Np2D, Tuple, HFEng
-from rrompy.utilities.exception_manager import RROMPyWarning, RROMPyException
 
 __all__ = ["compressMatrix"]
 
 def compressMatrix(A:Np2D, tol : float = 0., HFEngine : HFEng = None,
                    is_state : bool = True) -> Tuple[Np2D, Np2D, float]:
     if HFEngine is None or not is_state:
         U, s, _ = np.linalg.svd(A.T.conj().dot(A))
     else:
         U, s, _ = np.linalg.svd(HFEngine.innerProduct(A, A,
                                                       is_state = is_state))
     remove = np.where(s < tol * s[0])[0]
     ncut = len(s) if len(remove) == 0 else remove[0]
     sums = np.sum(s)
     s = s[: ncut] ** .5
     R = (U[:, : ncut].conj() * s).T
     U = dot(A, U[:, : ncut] * s ** -1.)
     return U, R, 1. - np.linalg.norm(s) / sums
diff --git a/rrompy/utilities/poly_fitting/heaviside/heaviside_to_from_affine.py b/rrompy/utilities/poly_fitting/heaviside/heaviside_to_from_affine.py
index 0045e23..fd1c49c 100644
--- a/rrompy/utilities/poly_fitting/heaviside/heaviside_to_from_affine.py
+++ b/rrompy/utilities/poly_fitting/heaviside/heaviside_to_from_affine.py
@@ -1,96 +1,95 @@
 # 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 scipy.special import binom
 import scipy.sparse as sp
 from rrompy.utilities.base.types import (Np1D, Np2D, List, ListAny, Tuple,
                                          paramVal)
 from rrompy.utilities.numerical import dot, solve
 from rrompy.utilities.numerical.nonlinear_eigenproblem import eigNonlinearDense
 from rrompy.utilities.exception_manager import RROMPyException
 from rrompy.parameter import checkParameter
 
 __all__ = ['heaviside2affine', 'affine2heaviside']
 
 def heaviside2affine(c:Np2D, poles:Np1D, mu : paramVal = [],
                      basis : str = "MONOMIAL_HEAVISIDE",
                      sparse : bool = False) \
                                         -> Tuple[Np2D, List[Np2D], List[Np1D]]:
     mu = checkParameter(mu, 1)(0, 0)
     n, d = len(poles), len(c) - len(poles)
     basisN = basis.split("_")[0]
     if basisN not in ["MONOMIAL", "CHEBYSHEV", "LEGENDRE"]:
         raise RROMPyException("Polynomial basis not recognized.")
     if sparse:
         A0 = sp.spdiags([np.concatenate((- mu - poles, np.ones(d)))],
                          [0], n + d, n + d)
         A1 = sp.spdiags([np.concatenate((np.ones(n), np.zeros(d)))],
                          [0], n + d, n + d)
     else:
         A0 = np.diag(np.concatenate((mu - poles, np.ones(d))))
         A1 = np.diag(np.concatenate((np.ones(n), np.zeros(d))))
     As = [A0, A1]
     bs = np.zeros((d, n + d), dtype = poles.dtype)
     bs[0, :] = 1.
     if d > 0:
         bs[0, n + 1 :] = 0.
         if d > 1:
             bs[1, n + 1] = 1.
             for j in range(2, d):
                 if basisN == "MONOMIAL":
                     bs[j, n + j] = 1.
                 else:
                     alpha = - 1. if basisN == "CHEBYSHEV" else 1. / j - 1.
                     bs[:, n + j] = alpha * bs[:, n + j - 2]
                     bs[1 :, n + j] += (1. - alpha) * bs[: -1, n + j - 1]
     bs = list(bs)
     return c.reshape(c.shape[0], -1).T, As, bs
 
 def affine2heaviside(As:ListAny, bs:ListAny,
                      jSupp : int = 1) -> Tuple[Np2D, Np1D, str]:
     if jSupp != 1 and not (isinstance(jSupp, (int,))
                        and jSupp.upper() == "COMPANION"):
         raise RROMPyException(("Affine to heaviside conversion does not allow "
                                "nonlinear eigenproblem support outside first "
                                "block row."))
     N = len(As)
     M = len(bs)
     n = As[0].shape[0]
     if N == 1:
         poles = np.empty(0, dtype = np.complex)
         Q = np.eye(n)
     else:
         basis = "MONOMIAL_HEAVISIDE"
         poles, P, Q = eigNonlinearDense(As, jSupp = jSupp,
                                         return_inverse = True)
         P = P[- n :, :]
         Q = Q[:, : n]
     bEffs = np.array([dot(Q, solve(As[-1], b, np.linalg.solve,
                                    {})) for b in bs])
     if N == 1:
         c = bEffs
     else:
         c = np.zeros((len(poles) + M - 1, As[0].shape[1]), dtype = np.complex)
         for l, pl in enumerate(poles):
             for i in range(M):
                 c[l, :] = pl ** i * bEffs[i, l] * P[:, l]
         for l in range(M - 1):
             for i in range(l + 1, M):
                 c[len(poles) + l, :] = dot(P, poles ** (i- 1 - l) * bEffs[i, :])
     return c, poles, basis
diff --git a/rrompy/utilities/poly_fitting/heaviside/vander.py b/rrompy/utilities/poly_fitting/heaviside/vander.py
index 3d0a2f7..d762104 100644
--- a/rrompy/utilities/poly_fitting/heaviside/vander.py
+++ b/rrompy/utilities/poly_fitting/heaviside/vander.py
@@ -1,87 +1,86 @@
 # 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.poly_fitting.polynomial.vander import (polyvander as pvP,
                                                        polyvanderTotal as pvTP)
-from rrompy.utilities.base.types import Np1D, Np2D, Tuple, List, paramList
+from rrompy.utilities.base.types import Np1D, Np2D, List, paramList
 from rrompy.parameter import checkParameterList
 from rrompy.utilities.exception_manager import RROMPyException
 
 __all__ = ['heavisidevander', 'polyvander', 'polyvanderTotal']
 
 def heavisidevander(x:paramList, poles:Np1D,
                     reorder : List[int] = None) -> Np2D:
     """Compute Heaviside-Vandermonde matrix."""
     x = checkParameterList(x, 1)[0]
     x_un, idx_un = x.unique(return_inverse = True)
     nx = len(x)
     if len(x_un) < nx:
         raise RROMPyException("Sample points must be distinct.")
     del x_un
     x = x.data.flatten()
     if reorder is not None: x = x[reorder]
     poles = poles.flatten()
     Van = np.empty((len(x), len(poles)), dtype = np.complex)
     for j in range(len(x)):
         Van[j, :] = 1. / (x[j] - poles)
     return Van
     
 def polyvander(x:paramList, poles:Np1D, degs : List[int] = None,
                basis : str = "MONOMIAL_HEAVISIDE",
                derIdxs : List[List[List[int]]] = None,
                reorder : List[int] = None, scl : Np1D = None) -> Np2D:
     """
     Compute full Hermite-Vandermonde matrix with specified derivative
         directions.
     """
     if derIdxs is not None and np.sum(np.sum(derIdxs)) > 0:
         raise RROMPyException(("Cannot take derivatives of heaviside "
                                "function."))
     basisp = basis.split("_")[0]
     VanH = heavisidevander(x, poles, reorder = reorder)
     if degs is None or np.sum(degs) < 0:
         VanP = np.empty((len(x), 0))
     else:
         VanP = pvP(x, degs, basisp, derIdxs = derIdxs, reorder = reorder,
                    scl = scl)
     return np.block([[VanH, VanP]])
 
 def polyvanderTotal(x:paramList, poles:Np1D, deg : int = None,
                     basis : str = "MONOMIAL_HEAVISIDE",
                     derIdxs : List[List[List[int]]] = None,
                     reorder : List[int] = None, scl : Np1D = None) -> Np2D:
     """
     Compute full total degree Hermite-Vandermonde matrix with specified
         derivative directions.
     """
     if derIdxs is not None and np.sum(np.sum(derIdxs)) > 0:
         raise RROMPyException(("Cannot take derivatives of radial basis "
                                "function."))
     basisp = basis.split("_")[0]
     VanR = heavisidevander(x, poles, reorder = reorder)
     if deg is None or deg < 0:
         VanP = np.empty((len(x), 0))
-        derIdxs, ordIdxs = np.zeros(0, dtype = int), np.zeros(0, dtype = int)
+        derIdxs, _ = np.zeros(0, dtype = int), np.zeros(0, dtype = int)
     else:
         VanP = pvTP(x, deg, basisp, derIdxs = derIdxs,
                     reorder = reorder, scl = scl)
-    ordIdxsEff = np.concatenate((np.arange(len(VanR)), ordIdxs + len(VanR)))
     return np.block([[VanR, VanP],
                      [VanP.T.conj(), np.zeros(tuple([VanP.shape[1]] * 2))]])