diff --git a/rrompy/parameter/parameter_sampling/generic_random_sampler.py b/rrompy/parameter/parameter_sampling/generic_random_sampler.py
index 8341c46..c1cdbe1 100644
--- a/rrompy/parameter/parameter_sampling/generic_random_sampler.py
+++ b/rrompy/parameter/parameter_sampling/generic_random_sampler.py
@@ -1,78 +1,78 @@
 # 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/>.
 #
 
 from abc import abstractmethod
 import numpy as np
 from .generic_sampler import GenericSampler
-from rrompy.utilities.base.types import List, DictAny, paramList
+from rrompy.utilities.base.types import Tuple, List, DictAny, paramList
 from rrompy.utilities.exception_manager import RROMPyException
 from rrompy.parameter.parameter_list import emptyParameterList
 
 _allowedRandomKinds = ["UNIFORM", "HALTON", "SOBOL"]
 
 __all__ = ['GenericRandomSampler']
 
 class GenericRandomSampler(GenericSampler):
     """Generator of random sample points."""
 
     def __init__(self, lims:paramList, kind : str = "UNIFORM",
                  parameterMap : DictAny = 1., refinementFactor : float = 1.,
                  seed : int = 42):
         super().__init__(lims = lims, parameterMap = parameterMap)
         self._allowedKinds = _allowedRandomKinds
         self.kind = kind
         self.refinementFactor = refinementFactor
         self.seed = seed
         self.reset()
 
     def __str__(self) -> str:
         return "{}_{}".format(super().__str__(), 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()
 
     def reset(self):
         self.points = emptyParameterList()
 
     def generatePoints(self, n:int, reorder : bool = True) -> paramList:
         """Array of quadrature points."""
         if self.kind == "UNIFORM":
             np.random.seed(self.seed)
         else:
             self.seedLoc = self.seed
         self.reset()
         rF, self.refinementFactor = self.refinementFactor, 1.
         _ = self.refine([None] * n)
         self.refinementFactor = rF
         return self.points
 
     @abstractmethod
-    def refine(self, active : List[int] = None) -> List[int]:
+    def refine(self, active : List[int] = None) -> Tuple[List[int], List[int]]:
         pass
diff --git a/rrompy/parameter/parameter_sampling/segment/random_sampler.py b/rrompy/parameter/parameter_sampling/segment/random_sampler.py
index c1515f5..a79cc26 100644
--- a/rrompy/parameter/parameter_sampling/segment/random_sampler.py
+++ b/rrompy/parameter/parameter_sampling/segment/random_sampler.py
@@ -1,55 +1,55 @@
 # 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/>.
 #
 
 from numbers import Number
 import numpy as np
 from rrompy.parameter.parameter_sampling.generic_random_sampler import (
                                                           GenericRandomSampler)
 from rrompy.utilities.numerical.halton import haltonGenerate
 from rrompy.utilities.numerical.sobol import sobolGenerate
-from rrompy.utilities.base.types import List
+from rrompy.utilities.base.types import Tuple, List
 
 __all__ = ['RandomSampler']
 
 class RandomSampler(GenericRandomSampler):
     """Generator of (quasi-)random sample points."""
 
-    def refine(self, active : List[int] = None) -> List[int]:
+    def refine(self, active : List[int] = None) -> Tuple[List[int], List[int]]:
         if active is None:
             n = self.npoints
         elif isinstance(active, (Number,)):
             n = active
         else:
             n = len(active)
         n = int(n * self.refinementFactor)
         if self.kind == "UNIFORM":
             xmat = np.random.uniform(size = (n, self.npar))
         elif self.kind == "HALTON":
             xmat, self.seedLoc = haltonGenerate(self.npar, n, self.seedLoc,
                                                 return_seed = True)
         else:
             xmat, self.seedLoc = sobolGenerate(self.npar, n, self.seedLoc,
                                                return_seed = True)
         limsE = self.mapParameterList(self.lims)
         for d in range(self.npar):
             a, b = limsE(d)
             xmat[:, d] = a + (b - a) * xmat[:, d]
         pts = self.mapParameterList(xmat, "B")
         idx = np.arange(n, dtype = int) + len(self.points)
         for pj in pts: self.points.append(pj)
-        return list(idx)
+        return list(idx), []
diff --git a/rrompy/parameter/parameter_sampling/shape/random_box_sampler.py b/rrompy/parameter/parameter_sampling/shape/random_box_sampler.py
index f1e872f..4142563 100644
--- a/rrompy/parameter/parameter_sampling/shape/random_box_sampler.py
+++ b/rrompy/parameter/parameter_sampling/shape/random_box_sampler.py
@@ -1,69 +1,69 @@
 # 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/>.
 #
 
 from numbers import Number
 import numpy as np
 from .generic_shape_random_sampler import GenericShapeRandomSampler
 from rrompy.utilities.numerical import haltonGenerate, sobolGenerate
-from rrompy.utilities.base.types import List
+from rrompy.utilities.base.types import Tuple, List
 
 __all__ = ['RandomBoxSampler']
 
 class RandomBoxSampler(GenericShapeRandomSampler):
     """Generator of (quasi-)random sample points on boxes."""
 
-    def refine(self, active : List[int] = None) -> List[int]:
+    def refine(self, active : List[int] = None) -> Tuple[List[int], List[int]]:
         if active is None:
             n = self.npoints
         elif isinstance(active, (Number,)):
             n = int(active)
         else:
             n = len(active)
         n = int(n * self.refinementFactor)
         nEff = int(np.ceil(n * np.prod(
                                    [max(x, 1. / x) for x in self.axisRatios])))
         xmat2 = []
         while len(xmat2) < n:
             if self.kind == "UNIFORM":
                 xmat2 = np.random.uniform(size = (nEff, 2 * self.npar))
             elif self.kind == "HALTON":
                 xmat2, self.seedLoc = haltonGenerate(2 * self.npar, nEff,
                                                      self.seedLoc,
                                                      return_seed = True)
             else:
                 xmat2, self.seedLoc = sobolGenerate(2 * self.npar, nEff,
                                                     self.seed,
                                                     return_seed = True)
             for d in range(self.npar):
                 ax = self.axisRatios[d]
                 if ax <= 1.:
                     xmat2 = xmat2[xmat2[:, 2 * d + 1] <= ax]
                 else:
                     xmat2 = xmat2[xmat2[:, 2 * d] <= 1. / ax]
                     xmat2[:, 2 * d : 2 * d + 2] *= ax
             nEff += 1
         xmat = np.empty((n, self.npar), dtype = np.complex)
         limsE = self.mapParameterList(self.lims)
         for d in range(self.npar):
             a, b = limsE(d)
             xmat[:, d] = a + (b - a) * (xmat2[: n, 2 * d]
                      + 1.j * self.axisRatios[d] * (xmat2[: n, 2 * d + 1] - .5))
         pts = self.mapParameterList(xmat, "B")
         idx = np.arange(n, dtype = int) + len(self.points)
         for pj in pts: self.points.append(pj)
-        return list(idx)
+        return list(idx), []
diff --git a/rrompy/parameter/parameter_sampling/shape/random_circle_sampler.py b/rrompy/parameter/parameter_sampling/shape/random_circle_sampler.py
index 5a1de1f..e484597 100644
--- a/rrompy/parameter/parameter_sampling/shape/random_circle_sampler.py
+++ b/rrompy/parameter/parameter_sampling/shape/random_circle_sampler.py
@@ -1,72 +1,72 @@
 # 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/>.
 #
 
 from numbers import Number
 import numpy as np
 from .generic_shape_random_sampler import GenericShapeRandomSampler
 from rrompy.utilities.numerical import (haltonGenerate, sobolGenerate,
                                         potential)
-from rrompy.utilities.base.types import List
+from rrompy.utilities.base.types import Tuple, List
 
 __all__ = ['RandomCircleSampler']
 
 class RandomCircleSampler(GenericShapeRandomSampler):
     """Generator of (quasi-)random sample points on ellipses."""
 
-    def refine(self, active : List[int] = None) -> List[int]:
+    def refine(self, active : List[int] = None) -> Tuple[List[int], List[int]]:
         if active is None:
             n = self.npoints
         elif isinstance(active, (Number,)):
             n = active
         else:
             n = len(active)
         n = int(n * self.refinementFactor)
         nEff = int(np.ceil(n * (4. / np.pi) ** self.npar * np.prod(
                                    [max(x, 1. / x) for x in self.axisRatios])))
         xmat2 = []
         while len(xmat2) < n:
             if self.kind == "UNIFORM":
                 xmat2 = np.random.uniform(size = (nEff, 2 * self.npar))
             elif self.kind == "HALTON":
                 xmat2, self.seedLoc = haltonGenerate(2 * self.npar, nEff,
                                                      self.seedLoc,
                                                      return_seed = True)
             else:
                 xmat2, self.seedLoc = sobolGenerate(2 * self.npar, nEff,
                                                     self.seed,
                                                     return_seed = True)
             xmat2 = xmat2 * 2. - 1.
             for d in range(self.npar):
                 ax = self.axisRatios[d]
                 if ax > 1.: xmat2[:, 2 * d : 2 * d + 2] *= ax
                 Z = xmat2[:, 2 * d] + 1.j * ax * xmat2[:, 2 * d + 1]
                 ptscore = potential(Z, self.normalFoci(d))
                 xmat2 = xmat2[ptscore <= self.groundPotential(d)]
             nEff += 1
         xmat = np.empty((n, self.npar), dtype = np.complex)
         limsE = self.mapParameterList(self.lims)
         for d in range(self.npar):
             ax = self.axisRatios[d]
             a, b = limsE(d)
             c, r = (a + b) / 2., (a - b) / 2.
             xmat[:, d] = c + r * (xmat2[: n, 2 * d]
                                 + 1.j * ax * xmat2[: n, 2 * d + 1])
         pts = self.mapParameterList(xmat, "B")
         idx = np.arange(n, dtype = int) + len(self.points)
         for pj in pts: self.points.append(pj)
-        return list(idx)
+        return list(idx), []
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 bdce203..5293190 100644
--- a/rrompy/parameter/parameter_sampling/sparse_grid/sparse_grid_sampler.py
+++ b/rrompy/parameter/parameter_sampling/sparse_grid/sparse_grid_sampler.py
@@ -1,119 +1,110 @@
 # 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/>.
 #
 
 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 List, DictAny, paramList
+from rrompy.utilities.base.types import Tuple, List, 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"])]:
             raise RROMPyException("Generator kind not recognized.")
         self._kind = kind.upper()
         self._noBoundary = "NOBOUNDARY" 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.deltadepth = np.zeros(1, dtype = int)
 
-    def refine(self, active : List[int] = None) -> List[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 = []
+        newIdxs, oldIdxs = [], []
         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(
-                                   self.mapParameterList(point[jdelta],
+            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(
+                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] + ddp
-                                               + self._noBoundary) * signdelta
-                        pointj[jdelta] = self.mapParameterList(
+                    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:
-                            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
+                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
+        return newIdxs, oldIdxs
 
     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)
+        while self.npoints < n: idx = self.refine(idx)[0]
         return self.points
diff --git a/rrompy/reduction_methods/pivoted/greedy/generic_pivoted_greedy_approximant.py b/rrompy/reduction_methods/pivoted/greedy/generic_pivoted_greedy_approximant.py
index 170b5e1..e398688 100644
--- a/rrompy/reduction_methods/pivoted/greedy/generic_pivoted_greedy_approximant.py
+++ b/rrompy/reduction_methods/pivoted/greedy/generic_pivoted_greedy_approximant.py
@@ -1,790 +1,735 @@
 # 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/>.
 #
 
 from abc import abstractmethod
 from copy import deepcopy as copy
 import numpy as np
 from matplotlib import pyplot as plt
 from rrompy.reduction_methods.pivoted.generic_pivoted_approximant import (
                                               GenericPivotedApproximantBase,
                                               GenericPivotedApproximantNoMatch,
                                               GenericPivotedApproximant)
 from rrompy.reduction_methods.pivoted.gather_pivoted_approximant import (
                                                       gatherPivotedApproximant)
 from rrompy.utilities.base.types import (Np1D, Np2D, Tuple, List, paramVal,
                                          paramList, ListAny)
 from rrompy.utilities.base import verbosityManager as vbMng
 from rrompy.utilities.numerical.point_matching import (pointMatching,
                                                        chordalMetricAdjusted)
 from rrompy.utilities.exception_manager import (RROMPyException, RROMPyAssert,
                                                 RROMPyWarning)
 from rrompy.parameter import emptyParameterList
 from rrompy.utilities.parallel import (masterCore, indicesScatter,
                                        arrayGatherv, isend)
 
 __all__ = ['GenericPivotedGreedyApproximantNoMatch',
            'GenericPivotedGreedyApproximant']
 
 class GenericPivotedGreedyApproximantBase(GenericPivotedApproximantBase):
-    _allowedEstimatorKindsMarginal = ["LEAVE_ONE_OUT", "LOOK_AHEAD",
-                                      "LOOK_AHEAD_RECOVER", "NONE"]
+    _allowedEstimatorKindsMarginal = ["LOOK_AHEAD", "LOOK_AHEAD_RECOVER",
+                                      "NONE"]
 
     def __init__(self, *args, **kwargs):
         self._preInit()
         self._addParametersToList(["matchingWeightError",
                                    "errorEstimatorKindMarginal",
                                    "greedyTolMarginal", "maxIterMarginal"],
                                   [0., "NONE", 1e-1, 1e2])
         super().__init__(*args, **kwargs)
         self._postInit()
 
     @property
     def scaleFactorDer(self):
         """Value of scaleFactorDer."""
         if self._scaleFactorDer == "NONE": return 1.
         if self._scaleFactorDer == "AUTO": return self._scaleFactorOldPivot
         return self._scaleFactorDer
     @scaleFactorDer.setter
     def scaleFactorDer(self, scaleFactorDer):
         if isinstance(scaleFactorDer, (str,)):
             scaleFactorDer = scaleFactorDer.upper()
         elif hasattr(scaleFactorDer, "__len__"):
             scaleFactorDer = list(scaleFactorDer)
         self._scaleFactorDer = scaleFactorDer
         self._approxParameters["scaleFactorDer"] = self._scaleFactorDer
 
     @property
     def samplerMarginal(self):
         """Value of samplerMarginal."""
         return self._samplerMarginal
     @samplerMarginal.setter
     def samplerMarginal(self, samplerMarginal):
         if 'refine' not in dir(samplerMarginal):
             raise RROMPyException("Marginal sampler type not recognized.")
         GenericPivotedApproximantBase.samplerMarginal.fset(self,
                                                            samplerMarginal)
 
     @property
     def errorEstimatorKindMarginal(self):
         """Value of errorEstimatorKindMarginal."""
         return self._errorEstimatorKindMarginal
     @errorEstimatorKindMarginal.setter
     def errorEstimatorKindMarginal(self, errorEstimatorKindMarginal):
         errorEstimatorKindMarginal = errorEstimatorKindMarginal.upper()
         if errorEstimatorKindMarginal not in (
                                           self._allowedEstimatorKindsMarginal):
             RROMPyWarning(("Marginal error estimator kind not recognized. "
                            "Overriding to 'NONE'."))
             errorEstimatorKindMarginal = "NONE"
         self._errorEstimatorKindMarginal = errorEstimatorKindMarginal
         self._approxParameters["errorEstimatorKindMarginal"] = (
                                                self.errorEstimatorKindMarginal)
 
     @property
     def matchingWeightError(self):
         """Value of matchingWeightError."""
         return self._matchingWeightError
     @matchingWeightError.setter
     def matchingWeightError(self, matchingWeightError):
         self._matchingWeightError = matchingWeightError
         self._approxParameters["matchingWeightError"] = (
                                                       self.matchingWeightError)
 
     @property
     def greedyTolMarginal(self):
         """Value of greedyTolMarginal."""
         return self._greedyTolMarginal
     @greedyTolMarginal.setter
     def greedyTolMarginal(self, greedyTolMarginal):
         if greedyTolMarginal < 0:
             raise RROMPyException("greedyTolMarginal must be non-negative.")
         if (hasattr(self, "_greedyTolMarginal")
         and self.greedyTolMarginal is not None):
             greedyTolMarginalold = self.greedyTolMarginal
         else:
             greedyTolMarginalold = -1
         self._greedyTolMarginal = greedyTolMarginal
         self._approxParameters["greedyTolMarginal"] = self.greedyTolMarginal
         if greedyTolMarginalold != self.greedyTolMarginal:
             self.resetSamples()
 
     @property
     def maxIterMarginal(self):
         """Value of maxIterMarginal."""
         return self._maxIterMarginal
     @maxIterMarginal.setter
     def maxIterMarginal(self, maxIterMarginal):
         if maxIterMarginal <= 0:
             raise RROMPyException("maxIterMarginal must be positive.")
         if (hasattr(self, "_maxIterMarginal")
         and self.maxIterMarginal is not None):
             maxIterMarginalold = self.maxIterMarginal
         else:
             maxIterMarginalold = -1
         self._maxIterMarginal = maxIterMarginal
         self._approxParameters["maxIterMarginal"] = self.maxIterMarginal
         if maxIterMarginalold != self.maxIterMarginal:
             self.resetSamples()
 
     def resetSamples(self):
         """Reset samples."""
         super().resetSamples()
         if not hasattr(self, "_temporaryPivot"):
             self._mus = emptyParameterList()
             self._musMarginal = emptyParameterList()
             if hasattr(self, "samplerMarginal"): self.samplerMarginal.reset()
         if hasattr(self, "samplingEngine") and self.samplingEngine is not None:
             self.samplingEngine.resetHistory()
 
     def _getDistanceApp(self, polesEx:Np1D, resEx:Np2D, muTest:paramVal,
                         foci:Tuple[float, float], ground:float) -> float:
         polesAp = self.trainedModel.interpolateMarginalPoles(muTest)[0]
         if self.matchingWeightError != 0:
             resAp = self.trainedModel.interpolateMarginalCoeffs(muTest)[0][
                                                              : len(polesAp), :]
             resEx = self.trainedModel.data.projMat[:,
                                                  : resEx.shape[1]].dot(resEx.T)
             resAp = self.trainedModel.data.projMat[:,
                                                  : resAp.shape[1]].dot(resAp.T)
         else:
             resAp = None
         dist = chordalMetricAdjusted(polesEx, polesAp,
                                      self.matchingWeightError, resEx, resAp,
                                      self.HFEngine, False)
         pmR, pmC = pointMatching(dist)
         return np.mean(dist[pmR, pmC])
     
-    def getErrorEstimatorMarginalLeaveOneOut(self) -> Np1D:
-        nTest = len(self.trainedModel.data.musMarginal)
-        self._musMarginalTestIdxs = np.arange(nTest)
-        if nTest <= 1:
-            err = np.empty(nTest)
-            err[:] = np.inf
-            return err
-        idx, sizes = indicesScatter(nTest, return_sizes = True)
-        err = []
-        if len(idx) > 0:
-            _tMdataFull = copy(self.trainedModel.data)
-            _musMExcl = None
-            self.verbosity -= 35
-            self.trainedModel.verbosity -= 35
-            foci = self.samplerPivot.normalFoci()
-            ground = self.samplerPivot.groundPotential()
-            for i, j in enumerate(idx):
-                jEff = j - (i > 0)
-                muTest = self.trainedModel.data.musMarginal[jEff]
-                polesEx = copy(self.trainedModel.data.HIs[jEff].poles)
-                idxGood = np.logical_not(np.logical_or(np.isinf(polesEx),
-                                                       np.isnan(polesEx)))
-                polesEx = polesEx[idxGood]
-                if self.matchingWeightError != 0:
-                    resEx = self.trainedModel.data.HIs[jEff].coeffs[
-                                                          np.where(idxGood)[0]]
-                else:
-                    resEx = None
-                if i > 0: self.musMarginal.insert(_musMExcl, j - 1)
-                _musMExcl = self.musMarginal[j]
-                self.musMarginal.pop(j)
-                if len(polesEx) == 0:
-                    err += [0.]
-                    continue
-                self._updateTrainedModelMarginalSamples([j])
-                self._finalizeMarginalization()
-                err += [self._getDistanceApp(polesEx, resEx, muTest,
-                                             foci, ground)]
-            self._updateTrainedModelMarginalSamples()
-            self.musMarginal.insert(_musMExcl, idx[-1])
-            self.verbosity += 35
-            self.trainedModel.verbosity += 35
-            self.trainedModel.data = _tMdataFull
-        return arrayGatherv(np.array(err), sizes)
-
     def getErrorEstimatorMarginalLookAhead(self) -> Np1D:
         if not hasattr(self.trainedModel, "_musMExcl"):
             err = np.zeros(0)
             err[:] = np.inf
             self._musMarginalTestIdxs = np.zeros(0, dtype = int)
             return err
         self._musMarginalTestIdxs = np.array(self.trainedModel._idxExcl,
                                              dtype = int)
         idx, sizes = indicesScatter(len(self.trainedModel._musMExcl),
                                     return_sizes = True)
         err = []
         if len(idx) > 0:
             self.verbosity -= 35
             self.trainedModel.verbosity -= 35
             foci = self.samplerPivot.normalFoci()
             ground = self.samplerPivot.groundPotential()
             for j in idx:
                 muTest = self.trainedModel._musMExcl[j]
                 HITest = self.trainedModel._HIsExcl[j]
                 polesEx = HITest.poles
                 idxGood = np.logical_not(np.logical_or(np.isinf(polesEx),
                                                        np.isnan(polesEx)))
                 polesEx = polesEx[idxGood]
                 if self.matchingWeightError != 0:
                     resEx = HITest.coeffs[np.where(idxGood)[0]]
                 else:
                     resEx = None
                 if len(polesEx) == 0:
                     err += [0.]
                     continue
                 err += [self._getDistanceApp(polesEx, resEx, muTest,
                                              foci, ground)]
             self.verbosity += 35
             self.trainedModel.verbosity += 35
         return arrayGatherv(np.array(err), sizes)
 
     def getErrorEstimatorMarginalNone(self) -> Np1D:
         nErr = len(self.trainedModel.data.musMarginal)
         self._musMarginalTestIdxs = np.arange(nErr)
         return (1. + self.greedyTolMarginal) * np.ones(nErr)
 
     def errorEstimatorMarginal(self, return_max : bool = False) -> Np1D:
         vbMng(self.trainedModel, "INIT",
               "Evaluating error estimator at mu = {}.".format(
                                        self.trainedModel.data.musMarginal), 10)
-        if self.errorEstimatorKindMarginal == "LEAVE_ONE_OUT":
-            err = self.getErrorEstimatorMarginalLeaveOneOut()
-        elif self.errorEstimatorKindMarginal[: 10] == "LOOK_AHEAD":
+        if self.errorEstimatorKindMarginal == "NONE":
+            nErr = len(self.trainedModel.data.musMarginal)
+            self._musMarginalTestIdxs = np.arange(nErr)
+            err = (1. + self.greedyTolMarginal) * np.ones(nErr)
+        else:#if self.errorEstimatorKindMarginal[: 10] == "LOOK_AHEAD":
             err = self.getErrorEstimatorMarginalLookAhead()
-        else:#if self.errorEstimatorKindMarginal == "NONE":
-            err = self.getErrorEstimatorMarginalNone()
         vbMng(self.trainedModel, "DEL", "Done evaluating error estimator", 10)
         if not return_max: return err
         idxMaxEst = np.where(err > self.greedyTolMarginal)[0]
         maxErr = err[idxMaxEst]
         if self.errorEstimatorKindMarginal == "NONE": maxErr = None
         return err, idxMaxEst, maxErr
 
     def plotEstimatorMarginal(self, est:Np1D, idxMax:List[int],
                               estMax:List[float]):
         if self.errorEstimatorKindMarginal == "NONE": return
         if (not (np.any(np.isnan(est)) or np.any(np.isinf(est)))
-        and masterCore()):
+        and masterCore() and hasattr(self.trainedModel, "_musMExcl")):
             fig = plt.figure(figsize = plt.figaspect(1. / self.nparMarginal))
             for jpar in range(self.nparMarginal):
                 ax = fig.add_subplot(1, self.nparMarginal, 1 + jpar)
-                if self.errorEstimatorKindMarginal == "LEAVE_ONE_OUT":
-                    musre = copy(self.trainedModel.data.musMarginal.re.data)
-                else:#if self.errorEstimatorKindMarginal[: 10] == "LOOK_AHEAD":
-                    if not hasattr(self.trainedModel, "_musMExcl"): return
-                    musre = np.real(self.trainedModel._musMExcl)
+                musre = np.real(self.trainedModel._musMExcl)
                 if len(idxMax) > 0 and estMax is not None:
                     maxrej = musre[idxMax, jpar]
                 errCP = copy(est)
                 idx = np.delete(np.arange(self.nparMarginal), jpar)
                 while len(musre) > 0:
                     if self.nparMarginal == 1:
                         currIdx = np.arange(len(musre))
                     else:
                         currIdx = np.where(np.isclose(np.sum(
                              np.abs(musre[:, idx] - musre[0, idx]), 1), 0.))[0]
                     currIdxSorted = currIdx[np.argsort(musre[currIdx, jpar])]
                     ax.semilogy(musre[currIdxSorted, jpar],
                                 errCP[currIdxSorted], 'k.-', linewidth = 1)
                     musre = np.delete(musre, currIdx, 0)
                     errCP = np.delete(errCP, currIdx)
                 ax.semilogy(self.musMarginal.re(jpar),
                             (self.greedyTolMarginal,) * len(self.musMarginal),
                             '*m')
                 if len(idxMax) > 0 and estMax is not None:
                     ax.semilogy(maxrej, estMax, 'xr')
                 ax.set_xlim(*list(self.samplerMarginal.lims.re(jpar)))
                 ax.grid()
             plt.tight_layout()
             plt.show()
 
     def _addMarginalSample(self, mus:paramList):
         mus = self.checkParameterListMarginal(mus)
         if len(mus) == 0: return
         self._nmusOld, nmus = len(self.musMarginal), len(mus)
         if (hasattr(self, "trainedModel") and self.trainedModel is not None
         and hasattr(self.trainedModel, "_musMExcl")):
             self._nmusOld += len(self.trainedModel._musMExcl)
         vbMng(self, "MAIN",
               ("Adding marginal sample point{} no. {}{} at {} to training "
                "set.").format("s" * (nmus > 1), self._nmusOld + 1,
                               "--{}".format(self._nmusOld + nmus) * (nmus > 1),
                               mus), 3)
         self.musMarginal.append(mus)
         self.setupApproxPivoted(mus)
         self._poleMatching()
         del self._nmusOld
         if (self.errorEstimatorKindMarginal[: 10] == "LOOK_AHEAD"
         and not self.firstGreedyIterM):
             ubRange = len(self.trainedModel.data.musMarginal)
             if hasattr(self.trainedModel, "_idxExcl"):
                 shRange = len(self.trainedModel._musMExcl)
             else:
                 shRange = 0
             testIdxs = list(range(ubRange + shRange - len(mus),
                                   ubRange + shRange))
             for j in testIdxs[::-1]:
                 self.musMarginal.pop(j - shRange)
             if hasattr(self.trainedModel, "_idxExcl"):
                 testIdxs = self.trainedModel._idxExcl + testIdxs
             self._updateTrainedModelMarginalSamples(testIdxs)
         self._finalizeMarginalization()
         self._SMarginal = len(self.musMarginal)
         self._approxParameters["SMarginal"] = self.SMarginal
         self.trainedModel.data.approxParameters["SMarginal"] = self.SMarginal
 
     def greedyNextSampleMarginal(self, muidx:List[int],
                                  plotEst : str = "NONE") \
                                     -> Tuple[Np1D, List[int], float, paramVal]:
         RROMPyAssert(self._mode, message = "Cannot add greedy sample.")
         muidx = self._musMarginalTestIdxs[muidx]
         if (self.errorEstimatorKindMarginal[: 10] == "LOOK_AHEAD"
         and not self.firstGreedyIterM):
             if not hasattr(self.trainedModel, "_idxExcl"):
                 raise RROMPyException(("Sample index to be added not present "
                                        "in trained model."))
             testIdxs = copy(self.trainedModel._idxExcl)
             skippedIdx = 0
             for cj, j in enumerate(self.trainedModel._idxExcl):
                 if j in muidx:
                     testIdxs.pop(skippedIdx)
                     self.musMarginal.insert(self.trainedModel._musMExcl[cj],
                                             j - skippedIdx)
                 else:
                     skippedIdx += 1
             if len(self.trainedModel._idxExcl) < (len(muidx)
                                                 + len(testIdxs)):
                 raise RROMPyException(("Sample index to be added not present "
                                        "in trained model."))
             self._updateTrainedModelMarginalSamples(testIdxs)
             self._SMarginal = len(self.musMarginal)
             self._approxParameters["SMarginal"] = self.SMarginal
             self.trainedModel.data.approxParameters["SMarginal"] = (
                                                                 self.SMarginal)
         self.firstGreedyIterM = False
-        idxAdded = self.samplerMarginal.refine(muidx)
+        idxAdded = self.samplerMarginal.refine(muidx)[0]
         self._addMarginalSample(self.samplerMarginal.points[idxAdded])
         errorEstTest, muidx, maxErrorEst = self.errorEstimatorMarginal(True)
         if plotEst == "ALL":
             self.plotEstimatorMarginal(errorEstTest, muidx, maxErrorEst)
         return (errorEstTest, muidx, maxErrorEst,
                 self.samplerMarginal.points[muidx])
                                               
     def _preliminaryTrainingMarginal(self):
         """Initialize starting snapshots of solution map."""
         RROMPyAssert(self._mode, message = "Cannot start greedy algorithm.")
         if np.sum(self.samplingEngine.nsamples) > 0: return
         self.resetSamples()
         self._addMarginalSample(self.samplerMarginal.generatePoints(
                                                                self.SMarginal))
 
     def _preSetupApproxPivoted(self, mus:paramList) \
                                            -> Tuple[ListAny, ListAny, ListAny]:
         self.computeScaleFactor()
         if self.trainedModel is None:
             self._setupTrainedModel(np.zeros((0, 0)))
             self.trainedModel.data.Qs, self.trainedModel.data.Ps = [], []
             self.trainedModel.data.Psupp = []
         self._trainedModelOld = copy(self.trainedModel)
         self._scaleFactorOldPivot = copy(self.scaleFactor)
         self.scaleFactor = self.scaleFactorPivot
         self._temporaryPivot = 1
         self._musLoc = copy(self.mus)
         idx, sizes = indicesScatter(len(mus), return_sizes = True)
         emptyCores = np.where(np.logical_not(sizes))[0]
         self.verbosity -= 15
         return idx, sizes, emptyCores
 
     def _postSetupApproxPivoted(self, mus:Np2D, pMat:Np2D, Ps:ListAny,
                                 Qs:ListAny, sizes:ListAny):
         self.scaleFactor = self._scaleFactorOldPivot
         del self._scaleFactorOldPivot, self._temporaryPivot
         pMat, Ps, Qs, mus, nsamples = gatherPivotedApproximant(pMat, Ps, Qs,
                                                                mus, sizes,
                                                                self.polybasis)
         if len(self._musLoc) > 0:
             self._mus = self.checkParameterList(self._musLoc)
             self._mus.append(mus)
         else:
             self._mus = self.checkParameterList(mus)
         self.trainedModel = self._trainedModelOld
         del self._trainedModelOld
         padLeft = self.trainedModel.data.projMat.shape[1]
         suppNew = np.append(0, np.cumsum(nsamples))
         self._setupTrainedModel(pMat, padLeft > 0)
         self.trainedModel.data.Qs += Qs
         self.trainedModel.data.Ps += Ps
         self.trainedModel.data.Psupp += list(padLeft + suppNew[: -1])
         self.trainedModel.data.approxParameters = copy(self.approxParameters)
         self.verbosity += 15
 
     def _localPivotedResult(self, pMat:Np2D, req:ListAny, emptyCores:ListAny,
                             mus:Np2D) -> Tuple[Np2D, ListAny, Np2D]:
         if pMat is None:
             mus = copy(self.samplingEngine.mus.data)
             pMat = copy(self.samplingEngine.projectionMatrix)
             if masterCore():
                 for dest in emptyCores:
                     req += [isend((len(pMat), pMat.dtype, mus.dtype),
                                   dest = dest, tag = dest)]
         else:
             mus = np.vstack((mus, self.samplingEngine.mus.data))
             pMat = np.hstack((pMat,
                               self.samplingEngine.projectionMatrix))
         return pMat, req, mus
     
     @abstractmethod
     def setupApproxPivoted(self, mus:paramList) -> int:
         if self.checkComputedApproxPivoted(): return -1
         RROMPyAssert(self._mode, message = "Cannot setup approximant.")
         vbMng(self, "INIT", "Setting up pivoted approximant.", 10)
         self._preSetupApproxPivoted()
         data = []
         pass
         self._postSetupApproxPivoted(mus, data)
         vbMng(self, "DEL", "Done setting up pivoted approximant.", 10)
         return 0
 
     def setupApprox(self, plotEst : str = "NONE") -> int:
         """Compute greedy snapshots of solution map."""
         if self.checkComputedApprox(): return -1
         RROMPyAssert(self._mode, message = "Cannot start greedy algorithm.")
         vbMng(self, "INIT", "Starting computation of snapshots.", 3)
         max2ErrorEst, self.firstGreedyIterM = np.inf, True
         self._preliminaryTrainingMarginal()
-        if self.errorEstimatorKindMarginal[: 10] == "LOOK_AHEAD":
-            muidx = np.arange(len(self.trainedModel.data.musMarginal))
-        else:#if self.errorEstimatorKindMarginal in ["LEAVE_ONE_OUT", "NONE"]:
+        if self.errorEstimatorKindMarginal == "NONE":
             muidx = []
+        else:#if self.errorEstimatorKindMarginal[: 10] == "LOOK_AHEAD":
+            muidx = np.arange(len(self.trainedModel.data.musMarginal))
         self._musMarginalTestIdxs = np.array(muidx)
         while self.firstGreedyIterM or (max2ErrorEst > self.greedyTolMarginal
                       and self.samplerMarginal.npoints < self.maxIterMarginal):
             errorEstTest, muidx, maxErrorEst, mu = \
                                   self.greedyNextSampleMarginal(muidx, plotEst)
             if maxErrorEst is None:
                 max2ErrorEst = 1. + self.greedyTolMarginal
             else:
                 if len(maxErrorEst) > 0:
                     max2ErrorEst = np.max(maxErrorEst)
                 else:
                     max2ErrorEst = np.max(errorEstTest)
                 vbMng(self, "MAIN", ("Uniform testing error estimate "
                                      "{:.4e}.").format(max2ErrorEst), 3)
         if plotEst == "LAST":
             self.plotEstimatorMarginal(errorEstTest, muidx, maxErrorEst)
         vbMng(self, "DEL", ("Done computing snapshots (final snapshot count: "
                             "{}).").format(len(self.mus)), 3)
         if (self.errorEstimatorKindMarginal == "LOOK_AHEAD_RECOVER"
         and hasattr(self.trainedModel, "_idxExcl")
         and len(self.trainedModel._idxExcl) > 0):
             vbMng(self, "INIT", "Recovering {} test models.".format(
                                            len(self.trainedModel._idxExcl)), 7)
             for j, mu in zip(self.trainedModel._idxExcl,
                              self.trainedModel._musMExcl):
                 self.musMarginal.insert(mu, j)
             self._updateTrainedModelMarginalSamples()
             self._finalizeMarginalization()
             self._SMarginal = len(self.musMarginal)
             self._approxParameters["SMarginal"] = self.SMarginal
             self.trainedModel.data.approxParameters["SMarginal"] = (
                                                                 self.SMarginal)
             vbMng(self, "DEL", "Done recovering test models.", 7)
         return 0
 
     def checkComputedApproxPivoted(self) -> bool:
         return (super().checkComputedApprox()
           and len(self.musMarginal) == len(self.trainedModel.data.musMarginal))
 
 class GenericPivotedGreedyApproximantNoMatch(
                                            GenericPivotedGreedyApproximantBase,
                                            GenericPivotedApproximantNoMatch):
     """
     ROM pivoted greedy interpolant computation for parametric problems (without
         pole matching) (ABSTRACT).
 
     Args:
         HFEngine: HF problem solver.
         mu0(optional): Default parameter. Defaults to 0.
         directionPivot(optional): Pivot components. Defaults to [0].
         approxParameters(optional): Dictionary containing values for main
             parameters of approximant. Recognized keys are:
             - 'POD': whether to compute POD of snapshots; defaults to True;
             - 'scaleFactorDer': scaling factors for derivative computation;
                 defaults to 'AUTO';
             - 'matchingWeightError': weight for pole matching optimization in
                 error estimation; defaults to 0;
             - 'S': total number of pivot samples current approximant relies
                 upon;
             - 'samplerPivot': pivot sample point generator;
             - 'SMarginal': number of starting marginal samples;
             - 'samplerMarginal': marginal sample point generator via sparse
                 grid;
             - 'errorEstimatorKindMarginal': kind of marginal error estimator;
-                available values include 'LEAVE_ONE_OUT', 'LOOK_AHEAD',
-                'LOOK_AHEAD_RECOVER', and 'NONE'; defaults to 'NONE';
+                available values include 'LOOK_AHEAD', 'LOOK_AHEAD_RECOVER',
+                and 'NONE'; defaults to 'NONE';
             - 'greedyTolMarginal': uniform error tolerance for marginal greedy
                algorithm; defaults to 1e-1;
             - 'maxIterMarginal': maximum number of marginal greedy steps;
                defaults to 1e2;
             - 'radialDirectionalWeightsMarginal': radial basis weights for
                 marginal interpolant; defaults to 1.
             Defaults to empty dict.
         approx_state(optional): Whether to approximate state. Defaults to
             False.
         verbosity(optional): Verbosity level. Defaults to 10.
 
     Attributes:
         HFEngine: HF problem solver.
         mu0: Default parameter.
         directionPivot: Pivot components.
         mus: Array of snapshot parameters.
         musMarginal: Array of marginal snapshot parameters.
         approxParameters: Dictionary containing values for main parameters of
             approximant. Recognized keys are in parameterList.
         parameterListSoft: Recognized keys of soft approximant parameters:
             - 'POD': whether to compute POD of snapshots;
             - 'scaleFactorDer': scaling factors for derivative computation;
             - 'matchingWeightError': weight for pole matching optimization in
                 error estimation;
             - 'errorEstimatorKindMarginal': kind of marginal error estimator;
             - 'greedyTolMarginal': uniform error tolerance for marginal greedy
                 algorithm;
             - 'maxIterMarginal': maximum number of marginal greedy steps;
             - 'radialDirectionalWeightsMarginal': radial basis weights for
                 marginal interpolant.
         parameterListCritical: Recognized keys of critical approximant
             parameters:
             - 'S': total number of pivot samples current approximant relies
                 upon;
             - 'samplerPivot': pivot sample point generator;
             - 'SMarginal': total number of marginal samples current approximant
                 relies upon;
             - 'samplerMarginal': marginal sample point generator via sparse
                 grid.
         approx_state: Whether to approximate state.
         verbosity: Verbosity level.
         POD: Whether to compute POD of snapshots.
         scaleFactorDer: Scaling factors for derivative computation.
         matchingWeightError: Weight for pole matching optimization in error
             estimation.
         S: Total number of pivot samples current approximant relies upon.
         samplerPivot: Pivot sample point generator.
         SMarginal: Total number of marginal samples current approximant relies
             upon.
         samplerMarginal: Marginal sample point generator via sparse grid.
         errorEstimatorKindMarginal: Kind of marginal error estimator.
         greedyTolMarginal: Uniform error tolerance for marginal greedy
             algorithm.
         maxIterMarginal: Maximum number of marginal greedy steps.
         radialDirectionalWeightsMarginal: Radial basis weights for marginal
             interpolant.
         muBounds: list of bounds for pivot parameter values.
         muBoundsMarginal: list of bounds for marginal parameter values.
         samplingEngine: Sampling engine.
         uHF: High fidelity solution(s) with parameter(s) lastSolvedHF as
             sampleList.
         lastSolvedHF: Parameter(s) corresponding to last computed high fidelity
             solution(s) as parameterList.
         uApproxReduced: Reduced approximate solution(s) with parameter(s)
             lastSolvedApprox as sampleList.
         lastSolvedApproxReduced: Parameter(s) corresponding to last computed
             reduced approximate solution(s) as parameterList.
         uApprox: Approximate solution(s) with parameter(s) lastSolvedApprox as
             sampleList.
         lastSolvedApprox: Parameter(s) corresponding to last computed
             approximate solution(s) as parameterList.
     """
 
     def _poleMatching(self):
         vbMng(self, "INIT", "Compressing poles.", 10)
         self.trainedModel.initializeFromRational()
         vbMng(self, "DEL", "Done compressing poles.", 10)
 
     def _updateTrainedModelMarginalSamples(self, idx : ListAny = []):
         self.trainedModel.updateEffectiveSamples(idx)
         
 class GenericPivotedGreedyApproximant(GenericPivotedGreedyApproximantBase,
                                       GenericPivotedApproximant):
     """
     ROM pivoted greedy interpolant computation for parametric problems (with
         pole matching) (ABSTRACT).
 
     Args:
         HFEngine: HF problem solver.
         mu0(optional): Default parameter. Defaults to 0.
         directionPivot(optional): Pivot components. Defaults to [0].
         approxParameters(optional): Dictionary containing values for main
             parameters of approximant. Recognized keys are:
             - 'POD': whether to compute POD of snapshots; defaults to True;
             - 'scaleFactorDer': scaling factors for derivative computation;
                 defaults to 'AUTO';
             - 'matchingWeight': weight for pole matching optimization; defaults
                 to 1;
             - 'matchingMode': mode for pole matching optimization; allowed
                 values include 'NONE' and 'SHIFT'; defaults to 'NONE';
             - 'sharedRatio': required ratio of marginal points to share
                 resonance; defaults to 1.;
             - 'matchingWeightError': weight for pole matching optimization in
                 error estimation; defaults to 0;
             - 'S': total number of pivot samples current approximant relies
                 upon;
             - 'samplerPivot': pivot sample point generator;
             - 'SMarginal': number of starting marginal samples;
             - 'samplerMarginal': marginal sample point generator via sparse
                 grid;
             - 'errorEstimatorKindMarginal': kind of marginal error estimator;
-                available values include 'LEAVE_ONE_OUT', 'LOOK_AHEAD',
-                'LOOK_AHEAD_RECOVER', and 'NONE'; defaults to 'NONE';
+                available values include 'LOOK_AHEAD', 'LOOK_AHEAD_RECOVER',
+                and 'NONE'; defaults to 'NONE';
             - 'polybasisMarginal': type of polynomial basis for marginal
                 interpolation; allowed values include 'MONOMIAL_*',
                 'CHEBYSHEV_*', 'LEGENDRE_*', 'NEARESTNEIGHBOR', and 
                 'PIECEWISE_LINEAR_*'; defaults to 'MONOMIAL';
             - 'paramsMarginal': dictionary of parameters for marginal
                 interpolation; include:
                 . 'MMarginal': degree of marginal interpolant; defaults to
                     'AUTO', i.e. maximum allowed; not for 'NEARESTNEIGHBOR' or
                     'PIECEWISE_LINEAR_*';
                 . 'nNeighborsMarginal': number of marginal nearest neighbors;
                      defaults to 1; only for 'NEARESTNEIGHBOR';
                 . 'polydegreetypeMarginal': type of polynomial degree for
                     marginal; defaults to 'TOTAL'; not for 'NEARESTNEIGHBOR' or
                     'PIECEWISE_LINEAR_*';
                 . 'interpRcondMarginal': tolerance for marginal interpolation;
                     defaults to None; not for 'NEARESTNEIGHBOR' or
                     'PIECEWISE_LINEAR_*';
                 . 'radialDirectionalWeightsMarginalAdapt': bounds for adaptive
                     rescaling of marginal radial basis weights; only for
                     radial basis.
             - 'greedyTolMarginal': uniform error tolerance for marginal greedy
                algorithm; defaults to 1e-1;
             - 'maxIterMarginal': maximum number of marginal greedy steps;
                defaults to 1e2;
             - 'radialDirectionalWeightsMarginal': radial basis weights for
                 marginal interpolant; defaults to 1.
             Defaults to empty dict.
         approx_state(optional): Whether to approximate state. Defaults to
             False.
         verbosity(optional): Verbosity level. Defaults to 10.
 
     Attributes:
         HFEngine: HF problem solver.
         mu0: Default parameter.
         directionPivot: Pivot components.
         mus: Array of snapshot parameters.
         musMarginal: Array of marginal snapshot parameters.
         approxParameters: Dictionary containing values for main parameters of
             approximant. Recognized keys are in parameterList.
         parameterListSoft: Recognized keys of soft approximant parameters:
             - 'POD': whether to compute POD of snapshots;
             - 'scaleFactorDer': scaling factors for derivative computation;
             - 'matchingWeight': weight for pole matching optimization;
             - 'matchingMode': mode for pole matching optimization;
             - 'sharedRatio': required ratio of marginal points to share
                 resonance;
             - 'matchingWeightError': weight for pole matching optimization in
                 error estimation;
             - 'errorEstimatorKindMarginal': kind of marginal error estimator;
             - 'polybasisMarginal': type of polynomial basis for marginal
                 interpolation;
             - 'paramsMarginal': dictionary of parameters for marginal
                 interpolation; include:
                 . 'MMarginal': degree of marginal interpolant;
                 . 'nNeighborsMarginal': number of marginal nearest neighbors;
                 . 'polydegreetypeMarginal': type of polynomial degree for
                     marginal;
                 . 'interpRcondMarginal': tolerance for marginal interpolation;
                 . 'radialDirectionalWeightsMarginalAdapt': bounds for adaptive
                     rescaling of marginal radial basis weights.
             - 'greedyTolMarginal': uniform error tolerance for marginal greedy
                 algorithm;
             - 'maxIterMarginal': maximum number of marginal greedy steps;
             - 'radialDirectionalWeightsMarginal': radial basis weights for
                 marginal interpolant.
         parameterListCritical: Recognized keys of critical approximant
             parameters:
             - 'S': total number of pivot samples current approximant relies
                 upon;
             - 'samplerPivot': pivot sample point generator;
             - 'SMarginal': total number of marginal samples current approximant
                 relies upon;
             - 'samplerMarginal': marginal sample point generator via sparse
                 grid.
         approx_state: Whether to approximate state.
         verbosity: Verbosity level.
         POD: Whether to compute POD of snapshots.
         scaleFactorDer: Scaling factors for derivative computation.
         matchingWeight: Weight for pole matching optimization.
         matchingMode: Mode for pole matching optimization.
         sharedRatio: Required ratio of marginal points to share resonance.
         matchingWeightError: Weight for pole matching optimization in error
             estimation.
         S: Total number of pivot samples current approximant relies upon.
         samplerPivot: Pivot sample point generator.
         SMarginal: Total number of marginal samples current approximant relies
             upon.
         samplerMarginal: Marginal sample point generator via sparse grid.
         errorEstimatorKindMarginal: Kind of marginal error estimator.
         polybasisMarginal: Type of polynomial basis for marginal interpolation.
         paramsMarginal: Dictionary of parameters for marginal interpolation.
         greedyTolMarginal: Uniform error tolerance for marginal greedy
             algorithm.
         maxIterMarginal: Maximum number of marginal greedy steps.
         radialDirectionalWeightsMarginal: Radial basis weights for marginal
             interpolant.
         muBounds: list of bounds for pivot parameter values.
         muBoundsMarginal: list of bounds for marginal parameter values.
         samplingEngine: Sampling engine.
         uHF: High fidelity solution(s) with parameter(s) lastSolvedHF as
             sampleList.
         lastSolvedHF: Parameter(s) corresponding to last computed high fidelity
             solution(s) as parameterList.
         uApproxReduced: Reduced approximate solution(s) with parameter(s)
             lastSolvedApprox as sampleList.
         lastSolvedApproxReduced: Parameter(s) corresponding to last computed
             reduced approximate solution(s) as parameterList.
         uApprox: Approximate solution(s) with parameter(s) lastSolvedApprox as
             sampleList.
         lastSolvedApprox: Parameter(s) corresponding to last computed
             approximate solution(s) as parameterList.
     """
 
     def _poleMatching(self):
         vbMng(self, "INIT", "Compressing and matching poles.", 10)
         self.trainedModel.initializeFromRational(self.matchingWeight,
                                                  self.matchingMode,
                                                  self.HFEngine, False)
         vbMng(self, "DEL", "Done compressing and matching poles.", 10)
 
     def _updateTrainedModelMarginalSamples(self, idx : ListAny = []):
         self.trainedModel.updateEffectiveSamples(idx, self.matchingWeight,
                                                  self.matchingMode,
                                                  self.HFEngine, False)
 
-    def getErrorEstimatorMarginalLeaveOneOut(self) -> Np1D:
-        if self.polybasisMarginal != "NEARESTNEIGHBOR":
-            if not hasattr(self, "_MMarginal_isauto"):
-                if not hasattr(self, "_MMarginalOriginal"):
-                    self._MMarginalOriginal = self.paramsMarginal["MMarginal"]
-                self.paramsMarginal["MMarginal"] = self._MMarginalOriginal
-            self._reduceDegreeNNoWarn = 1
-        err = super().getErrorEstimatorMarginalLeaveOneOut()
-        if self.polybasisMarginal != "NEARESTNEIGHBOR":
-            del self._reduceDegreeNNoWarn
-        return err
-
     def setupApprox(self, *args, **kwargs) -> int:
         if self.checkComputedApprox(): return -1
         self.purgeparamsMarginal()
-        return super().setupApprox(*args, **kwargs)
+        _polybasisMarginal = self.polybasisMarginal
+        self._polybasisMarginal = ("PIECEWISE_LINEAR_"
+                                 + self.samplerMarginal.kind)
+        setupOK = super().setupApprox(*args, **kwargs)
+        self._polybasisMarginal = _polybasisMarginal
+        self._finalizeMarginalization()
+        return setupOK
diff --git a/rrompy/reduction_methods/pivoted/greedy/rational_interpolant_greedy_pivoted_greedy.py b/rrompy/reduction_methods/pivoted/greedy/rational_interpolant_greedy_pivoted_greedy.py
index c1949e0..1e04f2a 100644
--- a/rrompy/reduction_methods/pivoted/greedy/rational_interpolant_greedy_pivoted_greedy.py
+++ b/rrompy/reduction_methods/pivoted/greedy/rational_interpolant_greedy_pivoted_greedy.py
@@ -1,502 +1,502 @@
 #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/>.
 #
 
 from copy import deepcopy as copy
 import numpy as np
 from .generic_pivoted_greedy_approximant import (
                                         GenericPivotedGreedyApproximantBase,
                                         GenericPivotedGreedyApproximantNoMatch,
                                         GenericPivotedGreedyApproximant)
 from rrompy.reduction_methods.standard.greedy import RationalInterpolantGreedy
 from rrompy.reduction_methods.standard.greedy.generic_greedy_approximant \
                                                             import pruneSamples
 from rrompy.reduction_methods.pivoted import (
                                        RationalInterpolantGreedyPivotedNoMatch,
                                        RationalInterpolantGreedyPivoted)
 from rrompy.utilities.base.types import Np1D, Tuple, paramVal, paramList
 from rrompy.utilities.base import verbosityManager as vbMng
 from rrompy.utilities.exception_manager import RROMPyAssert
 from rrompy.parameter import emptyParameterList
 from rrompy.utilities.parallel import poolRank, recv
 
 __all__ = ['RationalInterpolantGreedyPivotedGreedyNoMatch',
            'RationalInterpolantGreedyPivotedGreedy']
 
 class RationalInterpolantGreedyPivotedGreedyBase(
                                           GenericPivotedGreedyApproximantBase):
     @property
     def sampleBatchSize(self):
         """Value of sampleBatchSize."""
         return 1
 
     @property
     def sampleBatchIdx(self):
         """Value of sampleBatchIdx."""
         return self.S
 
     def greedyNextSample(self, muidx:int, plotEst : str = "NONE")\
                                           -> Tuple[Np1D, int, float, paramVal]:
         """Compute next greedy snapshot of solution map."""
         RROMPyAssert(self._mode, message = "Cannot add greedy sample.")
         mus = copy(self.muTest[muidx])
         self.muTest.pop(muidx)
         for j, mu in enumerate(mus):
             vbMng(self, "MAIN",
                   ("Adding sample point no. {} at {} to training "
                    "set.").format(len(self.mus) + 1, mu), 3)
             self.mus.append(mu)
             self._S = len(self.mus)
             self._approxParameters["S"] = self.S
             if (self.samplingEngine.nsamples <= len(mus) - j - 1
              or not np.allclose(mu, self.samplingEngine.mus[j - len(mus)])):
                 self.samplingEngine.nextSample(mu)
             if self._isLastSampleCollinear():
                 vbMng(self, "MAIN",
                       ("Collinearity above tolerance detected. Starting "
                        "preemptive greedy loop termination."), 3)
                 self._collinearityFlag = 1
                 errorEstTest = np.empty(len(self.muTest))
                 errorEstTest[:] = np.nan
                 return errorEstTest, [-1], np.nan, np.nan
         errorEstTest, muidx, maxErrorEst = self.errorEstimator(self.muTest,
                                                                True)
         if plotEst == "ALL":
             self.plotEstimator(errorEstTest, muidx, maxErrorEst)
         return errorEstTest, muidx, maxErrorEst, self.muTest[muidx]
 
     def _setSampleBatch(self, maxS:int):
         return self.S
 
     def _preliminaryTraining(self):
         """Initialize starting snapshots of solution map."""
         RROMPyAssert(self._mode, message = "Cannot start greedy algorithm.")
         if self.samplingEngine.nsamples > 0: return
         self.resetSamples()
         self.samplingEngine.scaleFactor = self.scaleFactorDer
         musPivot = self.trainSetGenerator.generatePoints(self.S)
         while len(musPivot) > self.S: musPivot.pop()
         muTestBasePivot = self.samplerPivot.generatePoints(self.nTestPoints,
                                                            False)
         idxPop = pruneSamples(self.mapParameterListPivot(muTestBasePivot),
                               self.mapParameterListPivot(musPivot),
                               1e-10 * self.scaleFactorPivot[0])
         muTestBasePivot.pop(idxPop)
         self._mus = emptyParameterList()
         self.mus.reset((self.S - 1, self.HFEngine.npar))
         self.muTest = emptyParameterList()
         self.muTest.reset((len(muTestBasePivot) + 1, self.HFEngine.npar))
         for k in range(self.S - 1):
             muk = np.empty_like(self.mus[0])
             muk[self.directionPivot] = musPivot[k]
             muk[self.directionMarginal] = self.muMargLoc
             self.mus[k] = muk
         for k in range(len(muTestBasePivot)):
             muk = np.empty_like(self.muTest[0])
             muk[self.directionPivot] = muTestBasePivot[k]
             muk[self.directionMarginal] = self.muMargLoc
             self.muTest[k] = muk
         muk = np.empty_like(self.mus[0])
         muk[self.directionPivot] = musPivot[-1]
         muk[self.directionMarginal] = self.muMargLoc
         self.muTest[-1] = muk
         if len(self.mus) > 0:
             vbMng(self, "MAIN", 
                   ("Adding first {} sample point{} at {} to training "
                    "set.").format(self.S - 1, "" + "s" * (self.S > 2),
                                   self.mus), 3)
             self.samplingEngine.iterSample(self.mus)
         self._S = len(self.mus)
         self._approxParameters["S"] = self.S
         self.M, self.N = ("AUTO",) * 2
 
     def setupApproxPivoted(self, mus:paramList) -> int:
         if self.checkComputedApproxPivoted(): return -1
         RROMPyAssert(self._mode, message = "Cannot setup approximant.")
         vbMng(self, "INIT", "Setting up pivoted approximant.", 10)
         if not hasattr(self, "_plotEstPivot"): self._plotEstPivot = "NONE"
         idx, sizes, emptyCores = self._preSetupApproxPivoted(mus)
         S0 = copy(self.S)
         pMat, Ps, Qs, req, musA = None, [], [], [], None
         if len(idx) == 0:
             vbMng(self, "MAIN", "Idling.", 45)
             if self.storeAllSamples: self.storeSamples()
             pL, pT, mT = recv(source = 0, tag = poolRank())
             pMat = np.empty((pL, 0), dtype = pT)
             musA = np.empty((0, self.mu0.shape[1]), dtype = mT)
         else:
             for i in idx:
                 self.muMargLoc = mus[i]
                 vbMng(self, "MAIN", "Building marginal model no. {} at "
                                     "{}.".format(i + 1, self.muMargLoc), 25)
                 self.samplingEngine.resetHistory()
                 self.trainedModel = None
                 self.verbosity -= 5
                 self.samplingEngine.verbosity -= 5
                 RationalInterpolantGreedy.setupApprox(self, self._plotEstPivot)
                 self.verbosity += 5
                 self.samplingEngine.verbosity += 5
                 if self.storeAllSamples: self.storeSamples(i + self._nmusOld)
                 pMat, req, musA = self._localPivotedResult(pMat, req,
                                                            emptyCores, musA)
                 Ps += [copy(self.trainedModel.data.P)]
                 Qs += [copy(self.trainedModel.data.Q)]
                 self._S = S0
             del self.muMargLoc
         for r in req: r.wait()
         self._postSetupApproxPivoted(musA, pMat, Ps, Qs, sizes)
         vbMng(self, "DEL", "Done setting up pivoted approximant.", 10)
         return 0
 
     def setupApprox(self, plotEst : str = "NONE") -> int:
         if self.checkComputedApprox(): return -1
         if '_' not in plotEst: plotEst = plotEst + "_NONE"
         plotEstM, self._plotEstPivot = plotEst.split("_")
         val = super().setupApprox(plotEstM)
         return val
 
 class RationalInterpolantGreedyPivotedGreedyNoMatch(
                                     RationalInterpolantGreedyPivotedGreedyBase,
                                     GenericPivotedGreedyApproximantNoMatch,
                                     RationalInterpolantGreedyPivotedNoMatch):
     """
     ROM greedy pivoted greedy rational interpolant computation for parametric
         problems (without pole matching).
 
     Args:
         HFEngine: HF problem solver.
         mu0(optional): Default parameter. Defaults to 0.
         directionPivot(optional): Pivot components. Defaults to [0].
         approxParameters(optional): Dictionary containing values for main
             parameters of approximant. Recognized keys are:
             - 'POD': whether to compute POD of snapshots; defaults to True;
             - 'scaleFactorDer': scaling factors for derivative computation;
                 defaults to 'AUTO';
             - 'matchingWeightError': weight for pole matching optimization in
                 error estimation; defaults to 0;
             - 'S': total number of pivot samples current approximant relies
                 upon;
             - 'samplerPivot': pivot sample point generator;
             - 'SMarginal': number of starting marginal samples;
             - 'samplerMarginal': marginal sample point generator via sparse
                 grid;
             - 'errorEstimatorKindMarginal': kind of marginal error estimator;
-                available values include 'LEAVE_ONE_OUT', 'LOOK_AHEAD', and
-                'LOOK_AHEAD_RECOVER'; defaults to 'LEAVE_ONE_OUT';
+                available values include 'LOOK_AHEAD' and 'LOOK_AHEAD_RECOVER';
+                defaults to 'NONE';
             - 'polybasis': type of polynomial basis for pivot interpolation;
                 defaults to 'MONOMIAL';
             - 'greedyTol': uniform error tolerance for greedy algorithm;
                 defaults to 1e-2;
             - 'collinearityTol': collinearity tolerance for greedy algorithm;
                 defaults to 0.;
             - 'maxIter': maximum number of greedy steps; defaults to 1e2;
             - 'nTestPoints': number of test points; defaults to 5e2;
             - 'trainSetGenerator': training sample points generator; defaults
                 to sampler;
             - 'errorEstimatorKind': kind of error estimator; available values
                 include 'AFFINE', 'DISCREPANCY', 'LOOK_AHEAD',
                 'LOOK_AHEAD_RES', and 'NONE'; defaults to 'NONE';
             - 'greedyTolMarginal': uniform error tolerance for marginal greedy
                algorithm; defaults to 1e-1;
             - 'maxIterMarginal': maximum number of marginal greedy steps;
                defaults to 1e2;
             - 'radialDirectionalWeightsMarginal': radial basis weights for
                 marginal interpolant; defaults to 1;
             - 'functionalSolve': strategy for minimization of denominator
                 functional; allowed values include 'NORM', 'DOMINANT', 'NODAL',
                 'LOEWNER', and 'BARYCENTRIC'; defaults to 'NORM';
             - 'interpRcond': tolerance for pivot interpolation; defaults to
                 None;
             - 'robustTol': tolerance for robust rational denominator
                 management; defaults to 0.
             Defaults to empty dict.
         approx_state(optional): Whether to approximate state. Defaults to
             False.
         verbosity(optional): Verbosity level. Defaults to 10.
             
     Attributes:
         HFEngine: HF problem solver.
         mu0: Default parameter.
         directionPivot: Pivot components.
         mus: Array of snapshot parameters.
         musMarginal: Array of marginal snapshot parameters.
         approxParameters: Dictionary containing values for main parameters of
             approximant. Recognized keys are in parameterList.
         parameterListSoft: Recognized keys of soft approximant parameters:
             - 'POD': whether to compute POD of snapshots;
             - 'scaleFactorDer': scaling factors for derivative computation;
             - 'matchingWeightError': weight for pole matching optimization in
                 error estimation;
             - 'errorEstimatorKindMarginal': kind of marginal error estimator;
             - 'polybasis': type of polynomial basis for pivot interpolation;
             - 'greedyTol': uniform error tolerance for greedy algorithm;
             - 'collinearityTol': collinearity tolerance for greedy algorithm;
             - 'maxIter': maximum number of greedy steps;
             - 'nTestPoints': number of test points;
             - 'trainSetGenerator': training sample points generator;
             - 'errorEstimatorKind': kind of error estimator;
             - 'greedyTolMarginal': uniform error tolerance for marginal greedy
                 algorithm;
             - 'maxIterMarginal': maximum number of marginal greedy steps;
             - 'radialDirectionalWeightsMarginal': radial basis weights for
                 marginal interpolant;
             - 'functionalSolve': strategy for minimization of denominator
                 functional;
             - 'interpRcond': tolerance for pivot interpolation;
             - 'robustTol': tolerance for robust rational denominator
                 management.
         parameterListCritical: Recognized keys of critical approximant
             parameters:
             - 'S': total number of pivot samples current approximant relies
                 upon;
             - 'samplerPivot': pivot sample point generator;
             - 'SMarginal': total number of marginal samples current approximant
                 relies upon;
             - 'samplerMarginal': marginal sample point generator via sparse
                 grid.
         approx_state: Whether to approximate state.
         verbosity: Verbosity level.
         POD: Whether to compute POD of snapshots.
         scaleFactorDer: Scaling factors for derivative computation.
         matchingWeightError: Weight for pole matching optimization in error
             estimation.
         S: Total number of pivot samples current approximant relies upon.
         samplerPivot: Pivot sample point generator.
         SMarginal: Total number of marginal samples current approximant relies
             upon.
         samplerMarginal: Marginal sample point generator via sparse grid.
         errorEstimatorKindMarginal: Kind of marginal error estimator.
         polybasis: Type of polynomial basis for pivot interpolation.
         greedyTol: uniform error tolerance for greedy algorithm.
         collinearityTol: Collinearity tolerance for greedy algorithm.
         maxIter: maximum number of greedy steps.
         nTestPoints: number of starting training points.
         trainSetGenerator: training sample points generator.
         errorEstimatorKind: kind of error estimator.
         greedyTolMarginal: Uniform error tolerance for marginal greedy
             algorithm.
         maxIterMarginal: Maximum number of marginal greedy steps.
         radialDirectionalWeightsMarginal: Radial basis weights for marginal
             interpolant.
         functionalSolve: Strategy for minimization of denominator functional.
         interpRcond: Tolerance for pivot interpolation.
         robustTol: Tolerance for robust rational denominator management.
         muBounds: list of bounds for pivot parameter values.
         muBoundsMarginal: list of bounds for marginal parameter values.
         samplingEngine: Sampling engine.
         uHF: High fidelity solution(s) with parameter(s) lastSolvedHF as
             sampleList.
         lastSolvedHF: Parameter(s) corresponding to last computed high fidelity
             solution(s) as parameterList.
         uApproxReduced: Reduced approximate solution(s) with parameter(s)
             lastSolvedApprox as sampleList.
         lastSolvedApproxReduced: Parameter(s) corresponding to last computed
             reduced approximate solution(s) as parameterList.
         uApprox: Approximate solution(s) with parameter(s) lastSolvedApprox as
             sampleList.
         lastSolvedApprox: Parameter(s) corresponding to last computed
             approximate solution(s) as parameterList.
     """
 
 class RationalInterpolantGreedyPivotedGreedy(
                                     RationalInterpolantGreedyPivotedGreedyBase,
                                     GenericPivotedGreedyApproximant,
                                     RationalInterpolantGreedyPivoted):
     """
     ROM greedy pivoted greedy rational interpolant computation for parametric
         problems (with pole matching).
 
     Args:
         HFEngine: HF problem solver.
         mu0(optional): Default parameter. Defaults to 0.
         directionPivot(optional): Pivot components. Defaults to [0].
         approxParameters(optional): Dictionary containing values for main
             parameters of approximant. Recognized keys are:
             - 'POD': whether to compute POD of snapshots; defaults to True;
             - 'scaleFactorDer': scaling factors for derivative computation;
                 defaults to 'AUTO';
             - 'matchingWeight': weight for pole matching optimization; defaults
                 to 1;
             - 'matchingMode': mode for pole matching optimization; allowed
                 values include 'NONE' and 'SHIFT'; defaults to 'NONE';
             - 'sharedRatio': required ratio of marginal points to share
                 resonance; defaults to 1.;
             - 'matchingWeightError': weight for pole matching optimization in
                 error estimation; defaults to 0;
             - 'S': total number of pivot samples current approximant relies
                 upon;
             - 'samplerPivot': pivot sample point generator;
             - 'SMarginal': number of starting marginal samples;
             - 'samplerMarginal': marginal sample point generator via sparse
                 grid;
             - 'errorEstimatorKindMarginal': kind of marginal error estimator;
-                available values include 'LEAVE_ONE_OUT', 'LOOK_AHEAD', and
-                'LOOK_AHEAD_RECOVER'; defaults to 'LEAVE_ONE_OUT';
+                available values include 'LOOK_AHEAD' and 'LOOK_AHEAD_RECOVER';
+                defaults to 'NONE';
             - 'polybasis': type of polynomial basis for pivot interpolation;
                 defaults to 'MONOMIAL';
             - 'polybasisMarginal': type of polynomial basis for marginal
                 interpolation; allowed values include 'MONOMIAL_*',
                 'CHEBYSHEV_*', 'LEGENDRE_*', 'NEARESTNEIGHBOR', and 
                 'PIECEWISE_LINEAR_*'; defaults to 'MONOMIAL';
             - 'paramsMarginal': dictionary of parameters for marginal
                 interpolation; include:
                 . 'MMarginal': degree of marginal interpolant; defaults to
                     'AUTO', i.e. maximum allowed; not for 'NEARESTNEIGHBOR' or
                     'PIECEWISE_LINEAR_*';
                 . 'nNeighborsMarginal': number of marginal nearest neighbors;
                      defaults to 1; only for 'NEARESTNEIGHBOR';
                 . 'polydegreetypeMarginal': type of polynomial degree for
                     marginal; defaults to 'TOTAL'; not for 'NEARESTNEIGHBOR' or
                     'PIECEWISE_LINEAR_*';
                 . 'interpRcondMarginal': tolerance for marginal interpolation;
                     defaults to None; not for 'NEARESTNEIGHBOR' or
                     'PIECEWISE_LINEAR_*';
                 . 'radialDirectionalWeightsMarginalAdapt': bounds for adaptive
                     rescaling of marginal radial basis weights; only for
                     radial basis.
             - 'greedyTol': uniform error tolerance for greedy algorithm;
                 defaults to 1e-2;
             - 'collinearityTol': collinearity tolerance for greedy algorithm;
                 defaults to 0.;
             - 'maxIter': maximum number of greedy steps; defaults to 1e2;
             - 'nTestPoints': number of test points; defaults to 5e2;
             - 'trainSetGenerator': training sample points generator; defaults
                 to sampler;
             - 'errorEstimatorKind': kind of error estimator; available values
                 include 'AFFINE', 'DISCREPANCY', 'LOOK_AHEAD',
                 'LOOK_AHEAD_RES', and 'NONE'; defaults to 'NONE';
             - 'greedyTolMarginal': uniform error tolerance for marginal greedy
                algorithm; defaults to 1e-1;
             - 'maxIterMarginal': maximum number of marginal greedy steps;
                defaults to 1e2;
             - 'radialDirectionalWeightsMarginal': radial basis weights for
                 marginal interpolant; defaults to 1;
             - 'functionalSolve': strategy for minimization of denominator
                 functional; allowed values include 'NORM', 'DOMINANT', 'NODAL',
                 'LOEWNER', and 'BARYCENTRIC'; defaults to 'NORM';
             - 'interpRcond': tolerance for pivot interpolation; defaults to
                 None;
             - 'robustTol': tolerance for robust rational denominator
                 management; defaults to 0.
             Defaults to empty dict.
         approx_state(optional): Whether to approximate state. Defaults to
             False.
         verbosity(optional): Verbosity level. Defaults to 10.
             
     Attributes:
         HFEngine: HF problem solver.
         mu0: Default parameter.
         directionPivot: Pivot components.
         mus: Array of snapshot parameters.
         musMarginal: Array of marginal snapshot parameters.
         approxParameters: Dictionary containing values for main parameters of
             approximant. Recognized keys are in parameterList.
         parameterListSoft: Recognized keys of soft approximant parameters:
             - 'POD': whether to compute POD of snapshots;
             - 'scaleFactorDer': scaling factors for derivative computation;
             - 'matchingWeight': weight for pole matching optimization;
             - 'matchingMode': mode for pole matching optimization;
             - 'sharedRatio': required ratio of marginal points to share
                 resonance;
             - 'matchingWeightError': weight for pole matching optimization in
                 error estimation;
             - 'errorEstimatorKindMarginal': kind of marginal error estimator;
             - 'polybasis': type of polynomial basis for pivot interpolation;
             - 'polybasisMarginal': type of polynomial basis for marginal
                 interpolation;
             - 'paramsMarginal': dictionary of parameters for marginal
                 interpolation; include:
                 . 'MMarginal': degree of marginal interpolant;
                 . 'nNeighborsMarginal': number of marginal nearest neighbors;
                 . 'polydegreetypeMarginal': type of polynomial degree for
                     marginal;
                 . 'interpRcondMarginal': tolerance for marginal interpolation;
                 . 'radialDirectionalWeightsMarginalAdapt': bounds for adaptive
                     rescaling of marginal radial basis weights.
             - 'greedyTol': uniform error tolerance for greedy algorithm;
             - 'collinearityTol': collinearity tolerance for greedy algorithm;
             - 'maxIter': maximum number of greedy steps;
             - 'nTestPoints': number of test points;
             - 'trainSetGenerator': training sample points generator;
             - 'errorEstimatorKind': kind of error estimator;
             - 'greedyTolMarginal': uniform error tolerance for marginal greedy
                 algorithm;
             - 'maxIterMarginal': maximum number of marginal greedy steps;
             - 'radialDirectionalWeightsMarginal': radial basis weights for
                 marginal interpolant;
             - 'functionalSolve': strategy for minimization of denominator
                 functional;
             - 'interpRcond': tolerance for pivot interpolation;
             - 'robustTol': tolerance for robust rational denominator
                 management.
         parameterListCritical: Recognized keys of critical approximant
             parameters:
             - 'S': total number of pivot samples current approximant relies
                 upon;
             - 'samplerPivot': pivot sample point generator;
             - 'SMarginal': total number of marginal samples current approximant
                 relies upon;
             - 'samplerMarginal': marginal sample point generator via sparse
                 grid.
         approx_state: Whether to approximate state.
         verbosity: Verbosity level.
         POD: Whether to compute POD of snapshots.
         scaleFactorDer: Scaling factors for derivative computation.
         matchingWeight: Weight for pole matching optimization.
         matchingMode: Mode for pole matching optimization.
         sharedRatio: Required ratio of marginal points to share resonance.
         matchingWeightError: Weight for pole matching optimization in error
             estimation.
         S: Total number of pivot samples current approximant relies upon.
         samplerPivot: Pivot sample point generator.
         SMarginal: Total number of marginal samples current approximant relies
             upon.
         samplerMarginal: Marginal sample point generator via sparse grid.
         errorEstimatorKindMarginal: Kind of marginal error estimator.
         polybasis: Type of polynomial basis for pivot interpolation.
         polybasisMarginal: Type of polynomial basis for marginal interpolation.
         paramsMarginal: Dictionary of parameters for marginal interpolation.
         greedyTol: uniform error tolerance for greedy algorithm.
         collinearityTol: Collinearity tolerance for greedy algorithm.
         maxIter: maximum number of greedy steps.
         nTestPoints: number of starting training points.
         trainSetGenerator: training sample points generator.
         errorEstimatorKind: kind of error estimator.
         greedyTolMarginal: Uniform error tolerance for marginal greedy
             algorithm.
         maxIterMarginal: Maximum number of marginal greedy steps.
         radialDirectionalWeightsMarginal: Radial basis weights for marginal
             interpolant.
         functionalSolve: Strategy for minimization of denominator functional.
         interpRcond: Tolerance for pivot interpolation.
         robustTol: Tolerance for robust rational denominator management.
         muBounds: list of bounds for pivot parameter values.
         muBoundsMarginal: list of bounds for marginal parameter values.
         samplingEngine: Sampling engine.
         uHF: High fidelity solution(s) with parameter(s) lastSolvedHF as
             sampleList.
         lastSolvedHF: Parameter(s) corresponding to last computed high fidelity
             solution(s) as parameterList.
         uApproxReduced: Reduced approximate solution(s) with parameter(s)
             lastSolvedApprox as sampleList.
         lastSolvedApproxReduced: Parameter(s) corresponding to last computed
             reduced approximate solution(s) as parameterList.
         uApprox: Approximate solution(s) with parameter(s) lastSolvedApprox as
             sampleList.
         lastSolvedApprox: Parameter(s) corresponding to last computed
             approximate solution(s) as parameterList.
     """
diff --git a/rrompy/reduction_methods/pivoted/greedy/rational_interpolant_pivoted_greedy.py b/rrompy/reduction_methods/pivoted/greedy/rational_interpolant_pivoted_greedy.py
index 4e0c59d..58b3408 100644
--- a/rrompy/reduction_methods/pivoted/greedy/rational_interpolant_pivoted_greedy.py
+++ b/rrompy/reduction_methods/pivoted/greedy/rational_interpolant_pivoted_greedy.py
@@ -1,428 +1,428 @@
 # 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/>.
 #
 
 from copy import deepcopy as copy
 from numpy import empty, empty_like
 from .generic_pivoted_greedy_approximant import (
                                         GenericPivotedGreedyApproximantBase,
                                         GenericPivotedGreedyApproximantNoMatch,
                                         GenericPivotedGreedyApproximant)
 from rrompy.reduction_methods.standard import RationalInterpolant
 from rrompy.reduction_methods.pivoted import (
                                              RationalInterpolantPivotedNoMatch,
                                              RationalInterpolantPivoted)
 from rrompy.utilities.base.types import paramList
 from rrompy.utilities.base import verbosityManager as vbMng
 from rrompy.utilities.exception_manager import RROMPyAssert
 from rrompy.parameter import emptyParameterList
 from rrompy.utilities.parallel import poolRank, recv
 
 __all__ = ['RationalInterpolantPivotedGreedyNoMatch',
            'RationalInterpolantPivotedGreedy']
 
 class RationalInterpolantPivotedGreedyBase(
                                           GenericPivotedGreedyApproximantBase):
     
     def computeSnapshots(self):
         """Compute snapshots of solution map."""
         RROMPyAssert(self._mode,
                      message = "Cannot start snapshot computation.")
         vbMng(self, "INIT", "Starting computation of snapshots.", 5)
         self.samplingEngine.scaleFactor = self.scaleFactorDer
         if not hasattr(self, "musPivot") or len(self.musPivot) != self.S:
             self.musPivot = self.samplerPivot.generatePoints(self.S)
             while len(self.musPivot) > self.S: self.musPivot.pop()
         musLoc = emptyParameterList()
         musLoc.reset((self.S, self.HFEngine.npar))
         self.samplingEngine.resetHistory()
         for k in range(self.S):
             muk = empty_like(musLoc[0])
             muk[self.directionPivot] = self.musPivot[k]
             muk[self.directionMarginal] = self.muMargLoc
             musLoc[k] = muk
         self.samplingEngine.iterSample(musLoc)
         vbMng(self, "DEL", "Done computing snapshots.", 5)
         self._m_selfmus = copy(musLoc)
         self._mus = self.musPivot
         self._m_HFEparameterMap = copy(self.HFEngine.parameterMap)
         self.HFEngine.parameterMap = {
                 "F": [self.HFEngine.parameterMap["F"][self.directionPivot[0]]],
                 "B": [self.HFEngine.parameterMap["B"][self.directionPivot[0]]]}
 
     def setupApproxPivoted(self, mus:paramList) -> int:
         if self.checkComputedApproxPivoted(): return -1
         RROMPyAssert(self._mode, message = "Cannot setup approximant.")
         vbMng(self, "INIT", "Setting up pivoted approximant.", 10)
         idx, sizes, emptyCores = self._preSetupApproxPivoted(mus)
         pMat, Ps, Qs, req, musA = None, [], [], [], None
         if len(idx) == 0:
             vbMng(self, "MAIN", "Idling.", 45)
             if self.storeAllSamples: self.storeSamples()
             pL, pT, mT = recv(source = 0, tag = poolRank())
             pMat = empty((pL, 0), dtype = pT)
             musA = empty((0, self.mu0.shape[1]), dtype = mT)
         else:
             for i in idx:
                 self.muMargLoc = mus[i]
                 vbMng(self, "MAIN", "Building marginal model no. {} at "
                                     "{}.".format(i + 1, self.muMargLoc), 25)
                 self.samplingEngine.resetHistory()
                 self.trainedModel = None
                 self.verbosity -= 5
                 self.samplingEngine.verbosity -= 5
                 RationalInterpolant.setupApprox(self)
                 self.verbosity += 5
                 self.samplingEngine.verbosity += 5
                 self._mus = self._m_selfmus
                 self.HFEngine.parameterMap = self._m_HFEparameterMap
                 del self._m_selfmus, self._m_HFEparameterMap
                 if self.storeAllSamples: self.storeSamples(i + self._nmusOld)
                 pMat, req, musA = self._localPivotedResult(pMat, req,
                                                            emptyCores, musA)
                 Ps += [copy(self.trainedModel.data.P)]
                 Qs += [copy(self.trainedModel.data.Q)]
             del self.muMargLoc
         for r in req: r.wait()
         self._postSetupApproxPivoted(musA, pMat, Ps, Qs, sizes)
         vbMng(self, "DEL", "Done setting up pivoted approximant.", 10)
         return 0
 
 class RationalInterpolantPivotedGreedyNoMatch(
                                         RationalInterpolantPivotedGreedyBase,
                                         GenericPivotedGreedyApproximantNoMatch,
                                         RationalInterpolantPivotedNoMatch):
     """
     ROM pivoted greedy rational interpolant computation for parametric
         problems (without pole matching).
 
     Args:
         HFEngine: HF problem solver.
         mu0(optional): Default parameter. Defaults to 0.
         directionPivot(optional): Pivot components. Defaults to [0].
         approxParameters(optional): Dictionary containing values for main
             parameters of approximant. Recognized keys are:
             - 'POD': whether to compute POD of snapshots; defaults to True;
             - 'scaleFactorDer': scaling factors for derivative computation;
                 defaults to 'AUTO';
             - 'matchingWeightError': weight for pole matching optimization in
                 error estimation; defaults to 0;
             - 'S': total number of pivot samples current approximant relies
                 upon;
             - 'samplerPivot': pivot sample point generator;
             - 'SMarginal': number of starting marginal samples;
             - 'samplerMarginal': marginal sample point generator via sparse
                 grid;
             - 'errorEstimatorKindMarginal': kind of marginal error estimator;
-                available values include 'LEAVE_ONE_OUT', 'LOOK_AHEAD', and
-                'LOOK_AHEAD_RECOVER'; defaults to 'LEAVE_ONE_OUT';
+                available values include 'LOOK_AHEAD' and 'LOOK_AHEAD_RECOVER';
+                defaults to 'NONE';
             - 'polybasis': type of polynomial basis for pivot interpolation;
                 defaults to 'MONOMIAL';
             - 'M': degree of rational interpolant numerator; defaults to
                 'AUTO', i.e. maximum allowed;
             - 'N': degree of rational interpolant denominator; defaults to
                 'AUTO', i.e. maximum allowed;
             - 'greedyTolMarginal': uniform error tolerance for marginal greedy
                algorithm; defaults to 1e-1;
             - 'maxIterMarginal': maximum number of marginal greedy steps;
                defaults to 1e2;
             - 'radialDirectionalWeights': radial basis weights for pivot
                 numerator; defaults to 1;
             - 'radialDirectionalWeightsAdapt': bounds for adaptive rescaling of
                 radial basis weights; defaults to [-1, -1];
             - 'radialDirectionalWeightsMarginal': radial basis weights for
                 marginal interpolant; defaults to 1;
             - 'functionalSolve': strategy for minimization of denominator
                 functional; allowed values include 'NORM', 'DOMINANT', 'NODAL',
                 'LOEWNER', and 'BARYCENTRIC'; defaults to 'NORM';
             - 'interpRcond': tolerance for pivot interpolation; defaults to
                 None;
             - 'robustTol': tolerance for robust rational denominator
                 management; defaults to 0.
             Defaults to empty dict.
         approx_state(optional): Whether to approximate state. Defaults to
             False.
         verbosity(optional): Verbosity level. Defaults to 10.
             
     Attributes:
         HFEngine: HF problem solver.
         mu0: Default parameter.
         directionPivot: Pivot components.
         mus: Array of snapshot parameters.
         musPivot: Array of pivot snapshot parameters.
         musMarginal: Array of marginal snapshot parameters.
         approxParameters: Dictionary containing values for main parameters of
             approximant. Recognized keys are in parameterList.
         parameterListSoft: Recognized keys of soft approximant parameters:
             - 'POD': whether to compute POD of snapshots;
             - 'scaleFactorDer': scaling factors for derivative computation;
             - 'matchingWeightError': weight for pole matching optimization in
                 error estimation;
             - 'errorEstimatorKindMarginal': kind of marginal error estimator;
             - 'polybasis': type of polynomial basis for pivot interpolation;
             - 'M': degree of rational interpolant numerator;
             - 'N': degree of rational interpolant denominator;
             - 'greedyTolMarginal': uniform error tolerance for marginal greedy
                 algorithm;
             - 'maxIterMarginal': maximum number of marginal greedy steps;
             - 'radialDirectionalWeights': radial basis weights for pivot
                 numerator;
             - 'radialDirectionalWeightsAdapt': bounds for adaptive rescaling of
                 radial basis weights;
             - 'radialDirectionalWeightsMarginal': radial basis weights for
                 marginal interpolant;
             - 'functionalSolve': strategy for minimization of denominator
                 functional;
             - 'interpRcond': tolerance for pivot interpolation;
             - 'robustTol': tolerance for robust rational denominator
                 management.
         parameterListCritical: Recognized keys of critical approximant
             parameters:
             - 'S': total number of pivot samples current approximant relies
                 upon;
             - 'samplerPivot': pivot sample point generator;
             - 'SMarginal': total number of marginal samples current approximant
                 relies upon;
             - 'samplerMarginal': marginal sample point generator via sparse
                 grid.
         approx_state: Whether to approximate state.
         verbosity: Verbosity level.
         POD: Whether to compute POD of snapshots.
         scaleFactorDer: Scaling factors for derivative computation.
         matchingWeightError: Weight for pole matching optimization in error
             estimation.
         S: Total number of pivot samples current approximant relies upon.
         samplerPivot: Pivot sample point generator.
         SMarginal: Total number of marginal samples current approximant relies
             upon.
         samplerMarginal: Marginal sample point generator via sparse grid.
         errorEstimatorKindMarginal: Kind of marginal error estimator.
         polybasis: Type of polynomial basis for pivot interpolation.
         M: Degree of rational interpolant numerator.
         N: Degree of rational interpolant denominator.
         greedyTolMarginal: Uniform error tolerance for marginal greedy
             algorithm.
         maxIterMarginal: Maximum number of marginal greedy steps.
         radialDirectionalWeights: Radial basis weights for pivot numerator.
         radialDirectionalWeightsAdapt: Bounds for adaptive rescaling of radial
             basis weights.
         radialDirectionalWeightsMarginal: Radial basis weights for marginal
             interpolant.
         functionalSolve: Strategy for minimization of denominator functional.
         interpRcond: Tolerance for pivot interpolation.
         robustTol: Tolerance for robust rational denominator management.
         muBounds: list of bounds for pivot parameter values.
         muBoundsMarginal: list of bounds for marginal parameter values.
         samplingEngine: Sampling engine.
         uHF: High fidelity solution(s) with parameter(s) lastSolvedHF as
             sampleList.
         lastSolvedHF: Parameter(s) corresponding to last computed high fidelity
             solution(s) as parameterList.
         uApproxReduced: Reduced approximate solution(s) with parameter(s)
             lastSolvedApprox as sampleList.
         lastSolvedApproxReduced: Parameter(s) corresponding to last computed
             reduced approximate solution(s) as parameterList.
         uApprox: Approximate solution(s) with parameter(s) lastSolvedApprox as
             sampleList.
         lastSolvedApprox: Parameter(s) corresponding to last computed
             approximate solution(s) as parameterList.
     """
 
 class RationalInterpolantPivotedGreedy(RationalInterpolantPivotedGreedyBase,
                                        GenericPivotedGreedyApproximant,
                                        RationalInterpolantPivoted):
     """
     ROM pivoted greedy rational interpolant computation for parametric
         problems (with pole matching).
 
     Args:
         HFEngine: HF problem solver.
         mu0(optional): Default parameter. Defaults to 0.
         directionPivot(optional): Pivot components. Defaults to [0].
         approxParameters(optional): Dictionary containing values for main
             parameters of approximant. Recognized keys are:
             - 'POD': whether to compute POD of snapshots; defaults to True;
             - 'scaleFactorDer': scaling factors for derivative computation;
                 defaults to 'AUTO';
             - 'matchingWeight': weight for pole matching optimization; defaults
                 to 1;
             - 'matchingMode': mode for pole matching optimization; allowed
                 values include 'NONE' and 'SHIFT'; defaults to 'NONE';
             - 'sharedRatio': required ratio of marginal points to share
                 resonance; defaults to 1.;
             - 'matchingWeightError': weight for pole matching optimization in
                 error estimation; defaults to 0;
             - 'S': total number of pivot samples current approximant relies
                 upon;
             - 'samplerPivot': pivot sample point generator;
             - 'SMarginal': number of starting marginal samples;
             - 'samplerMarginal': marginal sample point generator via sparse
                 grid;
             - 'errorEstimatorKindMarginal': kind of marginal error estimator;
-                available values include 'LEAVE_ONE_OUT', 'LOOK_AHEAD', and
-                'LOOK_AHEAD_RECOVER'; defaults to 'LEAVE_ONE_OUT';
+                available values include 'LOOK_AHEAD' and 'LOOK_AHEAD_RECOVER';
+                defaults to 'NONE';
             - 'polybasis': type of polynomial basis for pivot interpolation;
                 defaults to 'MONOMIAL';
             - 'polybasisMarginal': type of polynomial basis for marginal
                 interpolation; allowed values include 'MONOMIAL_*',
                 'CHEBYSHEV_*', 'LEGENDRE_*', 'NEARESTNEIGHBOR', and 
                 'PIECEWISE_LINEAR_*'; defaults to 'MONOMIAL';
             - 'paramsMarginal': dictionary of parameters for marginal
                 interpolation; include:
                 . 'MMarginal': degree of marginal interpolant; defaults to
                     'AUTO', i.e. maximum allowed; not for 'NEARESTNEIGHBOR' or
                     'PIECEWISE_LINEAR_*';
                 . 'nNeighborsMarginal': number of marginal nearest neighbors;
                      defaults to 1; only for 'NEARESTNEIGHBOR';
                 . 'polydegreetypeMarginal': type of polynomial degree for
                     marginal; defaults to 'TOTAL'; not for 'NEARESTNEIGHBOR' or
                     'PIECEWISE_LINEAR_*';
                 . 'interpRcondMarginal': tolerance for marginal interpolation;
                     defaults to None; not for 'NEARESTNEIGHBOR' or
                     'PIECEWISE_LINEAR_*';
                 . 'radialDirectionalWeightsMarginalAdapt': bounds for adaptive
                     rescaling of marginal radial basis weights; only for
                     radial basis.
             - 'M': degree of rational interpolant numerator; defaults to
                 'AUTO', i.e. maximum allowed;
             - 'N': degree of rational interpolant denominator; defaults to
                 'AUTO', i.e. maximum allowed;
             - 'greedyTolMarginal': uniform error tolerance for marginal greedy
                algorithm; defaults to 1e-1;
             - 'maxIterMarginal': maximum number of marginal greedy steps;
                defaults to 1e2;
             - 'radialDirectionalWeights': radial basis weights for pivot
                 numerator; defaults to 1;
             - 'radialDirectionalWeightsAdapt': bounds for adaptive rescaling of
                 radial basis weights; defaults to [-1, -1];
             - 'radialDirectionalWeightsMarginal': radial basis weights for
                 marginal interpolant; defaults to 1;
             - 'functionalSolve': strategy for minimization of denominator
                 functional; allowed values include 'NORM', 'DOMINANT', 'NODAL',
                 'LOEWNER', and 'BARYCENTRIC'; defaults to 'NORM';
             - 'interpRcond': tolerance for pivot interpolation; defaults to
                 None;
             - 'robustTol': tolerance for robust rational denominator
                 management; defaults to 0.
             Defaults to empty dict.
         approx_state(optional): Whether to approximate state. Defaults to
             False.
         verbosity(optional): Verbosity level. Defaults to 10.
             
     Attributes:
         HFEngine: HF problem solver.
         mu0: Default parameter.
         directionPivot: Pivot components.
         mus: Array of snapshot parameters.
         musPivot: Array of pivot snapshot parameters.
         musMarginal: Array of marginal snapshot parameters.
         approxParameters: Dictionary containing values for main parameters of
             approximant. Recognized keys are in parameterList.
         parameterListSoft: Recognized keys of soft approximant parameters:
             - 'POD': whether to compute POD of snapshots;
             - 'scaleFactorDer': scaling factors for derivative computation;
             - 'matchingWeight': weight for pole matching optimization;
             - 'matchingMode': mode for pole matching optimization;
             - 'sharedRatio': required ratio of marginal points to share
                 resonance;
             - 'matchingWeightError': weight for pole matching optimization in
                 error estimation;
             - 'errorEstimatorKindMarginal': kind of marginal error estimator;
             - 'polybasis': type of polynomial basis for pivot interpolation;
             - 'polybasisMarginal': type of polynomial basis for marginal
                 interpolation;
             - 'paramsMarginal': dictionary of parameters for marginal
                 interpolation; include:
                 . 'MMarginal': degree of marginal interpolant;
                 . 'nNeighborsMarginal': number of marginal nearest neighbors;
                 . 'polydegreetypeMarginal': type of polynomial degree for
                     marginal;
                 . 'interpRcondMarginal': tolerance for marginal interpolation;
                 . 'radialDirectionalWeightsMarginalAdapt': bounds for adaptive
                     rescaling of marginal radial basis weights.
             - 'M': degree of rational interpolant numerator;
             - 'N': degree of rational interpolant denominator;
             - 'greedyTolMarginal': uniform error tolerance for marginal greedy
                 algorithm;
             - 'maxIterMarginal': maximum number of marginal greedy steps;
             - 'radialDirectionalWeights': radial basis weights for pivot
                 numerator;
             - 'radialDirectionalWeightsAdapt': bounds for adaptive rescaling of
                 radial basis weights;
             - 'radialDirectionalWeightsMarginal': radial basis weights for
                 marginal interpolant;
             - 'functionalSolve': strategy for minimization of denominator
                 functional;
             - 'interpRcond': tolerance for pivot interpolation;
             - 'robustTol': tolerance for robust rational denominator
                 management.
         parameterListCritical: Recognized keys of critical approximant
             parameters:
             - 'S': total number of pivot samples current approximant relies
                 upon;
             - 'samplerPivot': pivot sample point generator;
             - 'SMarginal': total number of marginal samples current approximant
                 relies upon;
             - 'samplerMarginal': marginal sample point generator via sparse
                 grid.
         approx_state: Whether to approximate state.
         verbosity: Verbosity level.
         POD: Whether to compute POD of snapshots.
         scaleFactorDer: Scaling factors for derivative computation.
         matchingWeight: Weight for pole matching optimization.
         matchingMode: Mode for pole matching optimization.
         sharedRatio: Required ratio of marginal points to share resonance.
         matchingWeightError: Weight for pole matching optimization in error
             estimation.
         S: Total number of pivot samples current approximant relies upon.
         samplerPivot: Pivot sample point generator.
         SMarginal: Total number of marginal samples current approximant relies
             upon.
         samplerMarginal: Marginal sample point generator via sparse grid.
         errorEstimatorKindMarginal: Kind of marginal error estimator.
         polybasis: Type of polynomial basis for pivot interpolation.
         polybasisMarginal: Type of polynomial basis for marginal interpolation.
         paramsMarginal: Dictionary of parameters for marginal interpolation.
         M: Degree of rational interpolant numerator.
         N: Degree of rational interpolant denominator.
         greedyTolMarginal: Uniform error tolerance for marginal greedy
             algorithm.
         maxIterMarginal: Maximum number of marginal greedy steps.
         radialDirectionalWeights: Radial basis weights for pivot numerator.
         radialDirectionalWeightsAdapt: Bounds for adaptive rescaling of radial
             basis weights.
         radialDirectionalWeightsMarginal: Radial basis weights for marginal
             interpolant.
         functionalSolve: Strategy for minimization of denominator functional.
         interpRcond: Tolerance for pivot interpolation.
         robustTol: Tolerance for robust rational denominator management.
         muBounds: list of bounds for pivot parameter values.
         muBoundsMarginal: list of bounds for marginal parameter values.
         samplingEngine: Sampling engine.
         uHF: High fidelity solution(s) with parameter(s) lastSolvedHF as
             sampleList.
         lastSolvedHF: Parameter(s) corresponding to last computed high fidelity
             solution(s) as parameterList.
         uApproxReduced: Reduced approximate solution(s) with parameter(s)
             lastSolvedApprox as sampleList.
         lastSolvedApproxReduced: Parameter(s) corresponding to last computed
             reduced approximate solution(s) as parameterList.
         uApprox: Approximate solution(s) with parameter(s) lastSolvedApprox as
             sampleList.
         lastSolvedApprox: Parameter(s) corresponding to last computed
             approximate solution(s) as parameterList.
     """
diff --git a/tests/4_reduction_methods_multiD/greedy_pivoted_rational_2d.py b/tests/4_reduction_methods_multiD/greedy_pivoted_rational_2d.py
index 25e239f..448a5d4 100644
--- a/tests/4_reduction_methods_multiD/greedy_pivoted_rational_2d.py
+++ b/tests/4_reduction_methods_multiD/greedy_pivoted_rational_2d.py
@@ -1,87 +1,87 @@
 # 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 matrix_random import matrixRandom
 from rrompy.reduction_methods import (
                                RationalInterpolantPivotedGreedy as RIPG,
                                RationalInterpolantGreedyPivotedGreedy as RIGPG)
 from rrompy.parameter.parameter_sampling import (QuadratureSampler as QS,
                                                  SparseGridSampler as SGS)
 
 def test_pivoted_greedy():
     mu = [5.05, 7.1]
     mu0 = [5., 7.]
     solver = matrixRandom()
     uh = solver.solve(mu)[0]
     params = {"POD": True, "S": 5, "polybasis": "CHEBYSHEV",
               "samplerPivot": QS([4.75, 5.25], "CHEBYSHEV"),
               "SMarginal": 3, "greedyTolMarginal": 1e-2,
               "radialDirectionalWeightsMarginal": 2.,
               "polybasisMarginal": "MONOMIAL_GAUSSIAN",
               "paramsMarginal":{"MMarginal": 1},
-              "errorEstimatorKindMarginal": "LEAVE_ONE_OUT",
+              "errorEstimatorKindMarginal": "LOOK_AHEAD_RECOVER",
               "matchingWeight": 1., "samplerMarginal":SGS([6.75, 7.25])}
     approx = RIPG([0], solver, mu0, approx_state = True,
                   approxParameters = params, verbosity = 0)
     approx.setupApprox()
 
     uhP1 = approx.getApprox(mu)[0]
     errP = approx.getErr(mu)[0]
     errNP = approx.normErr(mu)[0]
     myerrP = uhP1 - uh
     assert np.allclose(np.abs(errP - myerrP), 0., rtol = 1e-3)
     assert np.isclose(solver.norm(errP), errNP, rtol = 1e-3)
     resP = approx.getRes(mu)[0]
     resNP = approx.normRes(mu)
     assert np.isclose(solver.norm(resP), resNP, rtol = 1e-3)
     assert np.allclose(np.abs(resP - (solver.b(mu) - solver.A(mu).dot(uhP1))),
                        0., rtol = 1e-3)
     assert np.isclose(errNP / solver.norm(uh), 6.0631706e-04, rtol = 1)
 
 def test_greedy_pivoted_greedy():
     mu = [5.05, 7.1]
     mu0 = [5., 7.]
     solver = matrixRandom()
     uh = solver.solve(mu)[0]
     params = {"POD": True, "nTestPoints": 100, "greedyTol": 1e-3, "S": 2,
               "polybasis": "CHEBYSHEV", 
               "samplerPivot": QS([4.75, 5.25], "CHEBYSHEV"),
               "trainSetGenerator": QS([4.75, 5.25], "CHEBYSHEV"),
               "SMarginal": 3, "greedyTolMarginal": 1e-2,
               "radialDirectionalWeightsMarginal": 2.,
               "polybasisMarginal": "MONOMIAL_GAUSSIAN",
               "paramsMarginal":{"MMarginal": 1},
-              "errorEstimatorKindMarginal": "LEAVE_ONE_OUT",
+              "errorEstimatorKindMarginal": "LOOK_AHEAD_RECOVER",
               "matchingWeight": 1., "samplerMarginal":SGS([6.75, 7.25])}
     approx = RIGPG([0], solver, mu0, approx_state = True,
                    approxParameters = params, verbosity = 0)
     approx.setupApprox()
 
     uhP1 = approx.getApprox(mu)[0]
     errP = approx.getErr(mu)[0]
     errNP = approx.normErr(mu)[0]
     myerrP = uhP1 - uh
     assert np.allclose(np.abs(errP - myerrP), 0., rtol = 1e-3)
     assert np.isclose(solver.norm(errP), errNP, rtol = 1e-3)
     resP = approx.getRes(mu)[0]
     resNP = approx.normRes(mu)
     assert np.isclose(solver.norm(resP), resNP, rtol = 1e-3)
     assert np.allclose(np.abs(resP - (solver.b(mu) - solver.A(mu).dot(uhP1))),
                        0., rtol = 1e-3)
     assert np.isclose(errNP / solver.norm(uh), 6.0631706e-04, rtol = 1)