diff --git a/examples/3d/base/membrane_fracture_engine3.py b/examples/3d/base/membrane_fracture_engine3.py
index f07eb34..29db01c 100644
--- a/examples/3d/base/membrane_fracture_engine3.py
+++ b/examples/3d/base/membrane_fracture_engine3.py
@@ -1,262 +1,262 @@
 # 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
 import fenics as fen
 import ufl
 from rrompy.utilities.base.types import paramVal
 from rrompy.solver.fenics import fenZERO, fenONE
 from rrompy.hfengines.linear_problem.helmholtz_problem_engine import (
                                                         HelmholtzProblemEngine)
 from rrompy.utilities.base import verbosityManager as vbMng
 from rrompy.utilities.numerical import (nextDerivativeIndices,
                                         hashDerivativeToIdx as hashD)
 from rrompy.solver.fenics import fenics2Sparse
 
 def gen_mesh_fracture(W, delta, H, n):
     n = 2 * (n // 2) + 1
     xL = np.linspace(-.5 * W, - .5 * delta, n // 2 + 1)
     xR = - xL[::-1]
     nEff = int(np.ceil(delta / (xL[1] - xL[0])))
     xC = np.linspace(-.5 * delta, .5 * delta, nEff + 1)[1:-1]
     yA = np.linspace(-.5 * H, .5 * H, n)
     yL = np.linspace(-.5 * H, 0., n // 2 + 1)
     editor = fen.MeshEditor()
     mesh = fen.Mesh()
     editor.open(mesh, "triangle", 2, 2)
     editor.init_vertices((2 * n + nEff - 1) * (n // 2 + 1))
     editor.init_cells(2 * (2 * n - 2 + nEff) * (n // 2))
     for j, y in enumerate(yA):
         for i, x in enumerate(xL):
             editor.add_vertex(i + j * (n // 2 + 1), np.array([x, y]))
     for j in range(n - 1):
         for i in range(n // 2):
             editor.add_cell(2 * (i + j * (n // 2)),
                       np.array([i + j * (n // 2 + 1), i + 1 + j * (n // 2 + 1),
                                 i + (j + 1) * (n // 2 + 1)], dtype=np.uintp))
             editor.add_cell(2 * (i + j * (n // 2)) + 1,
             np.array([i + 1 + j * (n // 2 + 1), i + 1 + (j + 1) * (n // 2 + 1),
                       i + (j + 1) * (n // 2 + 1)], dtype=np.uintp))
     vBase1, cBase1 = n * (n // 2 + 1), 2 * (n - 1) * (n // 2)
     for j, y in enumerate(yA):
         for i, x in enumerate(xR):
             editor.add_vertex(vBase1 + i + j * (n // 2 + 1), np.array([x, y]))
     for j in range(n - 1):
         for i in range(n // 2):
             editor.add_cell(cBase1 + 2 * (i + j * (n // 2)),
     np.array([vBase1 + i + j * (n // 2 + 1), vBase1 + i + 1 + j * (n // 2 + 1),
               vBase1 + i + (j + 1) * (n // 2 + 1)], dtype=np.uintp))
             editor.add_cell(cBase1 + 2 * (i + j * (n // 2)) + 1,
                np.array([vBase1 + i + 1 + j * (n // 2 + 1),
                          vBase1 + i + 1 + (j + 1) * (n // 2 + 1),
                          vBase1 + i + (j + 1) * (n // 2 + 1)], dtype=np.uintp))
     vBase2, cBase2 = 2 * n * (n // 2 + 1), 4 * (n - 1) * (n // 2)
     for j, y in enumerate(yL):
         for i, x in enumerate(xC):
             editor.add_vertex(vBase2 + i + j * (nEff - 1), np.array([x, y]))
     for j in range(n // 2):
         for i in range(nEff - 2):
             editor.add_cell(cBase2 + 2 * (i + j * (nEff - 2)),
         np.array([vBase2 + i + j * (nEff - 1), vBase2 + i + 1 + j * (nEff - 1),
                   vBase2 + i + (j + 1) * (nEff - 1)], dtype=np.uintp))
             editor.add_cell(cBase2 + 2 * (i + j * (nEff - 2)) + 1,
                  np.array([vBase2 + i + 1 + j * (nEff - 1),
                            vBase2 + i + 1 + (j + 1) * (nEff - 1),
                            vBase2 + i + (j + 1) * (nEff - 1)], dtype=np.uintp))
     if nEff == 1:
         for j in range(n // 2):
             editor.add_cell(cBase2 + 2 * j,
                np.array([(j + 1) * (n // 2 + 1) - 1, vBase1 + j * (n // 2 + 1),
                          (j + 2) * (n // 2 + 1) - 1], dtype=np.uintp))
             editor.add_cell(cBase2 + 2 * j + 1,
           np.array([vBase1 + j * (n // 2 + 1), vBase1 + (j + 1) * (n // 2 + 1),
                     (j + 2) * (n // 2 + 1) - 1], dtype=np.uintp))
     else:
         cBase3 = 2 * (2 * n + nEff - 4) * (n // 2)
         for j in range(n // 2):
             editor.add_cell(cBase3 + 2 * j,
                  np.array([(j + 1) * (n // 2 + 1) - 1, vBase2 + j * (nEff - 1),
                            (j + 2) * (n // 2 + 1) - 1], dtype=np.uintp))
             editor.add_cell(cBase3 + 2 * j + 1,
               np.array([vBase2 + j * (nEff - 1), vBase2 + (j + 1) * (nEff - 1),
                         (j + 2) * (n // 2 + 1) - 1], dtype=np.uintp))
         cBase4 = 2 * (2 * n + nEff - 3) * (n // 2)
         for j in range(n // 2):
             editor.add_cell(cBase4 + 2 * j,
         np.array([vBase2 + (j + 1) * (nEff - 1) - 1, vBase1 + j * (n // 2 + 1),
                   vBase2 + (j + 2) * (nEff - 1) - 1], dtype=np.uintp))
             editor.add_cell(cBase4 + 2 * j + 1,
           np.array([vBase1 + j * (n // 2 + 1), vBase1 + (j + 1) * (n // 2 + 1),
                     vBase2 + (j + 2) * (nEff - 1) - 1], dtype=np.uintp))
     editor.close()
     return mesh
 
 class MembraneFractureEngine3(HelmholtzProblemEngine):
 
     def __init__(self, mu0 : paramVal = [20. ** .5, .6, .08], H : float = 1.,
                  L : float = .75, delta : float = .05, n : int = 50,
                  degree_threshold : int = np.inf, verbosity : int = 10,
                  timestamp : bool = True):
         super().__init__(mu0 = mu0, degree_threshold = degree_threshold,
                          verbosity = verbosity, timestamp = timestamp)
         self._affinePoly = False
         self.nAs = 6
         self.npar = 3
         self._H = H
         self._delta = delta
         self._L = L
         self._W = 2 * L + delta
         self.rescalingExp = [2., 1., 1.]
         mesh = gen_mesh_fracture(self._W, delta, H, n)
         self.V = fen.FunctionSpace(mesh, "P", 1)
         self.NeumannBoundary = lambda x, on_b: (on_b and x[1] >= - H / 4.
                                and x[0] >= - .5 * delta and x[0] <= .5 * delta)
         self.DirichletBoundary = "REST"
 
         x, y = fen.SpatialCoordinate(mesh)[:]
         self._aboveIndicator = ufl.conditional(ufl.gt(y, 0.), fenONE, fenZERO)
         self._belowCenterIndicator = ufl.conditional(
                                                ufl.And(ufl.ge(x, - .5 * delta),
                                                        ufl.le(x, .5 * delta)),
                                                fenONE, fenZERO)
         self._belowSidesIndicator = (fenONE - self._aboveIndicator
                                             - self._belowCenterIndicator)
         self.DirichletDatum = [fen.exp(- 10. * (H / 2. + y) / H
                             - .5 * ((x + .4 * L + .5 * delta) / (.1 * L)) ** 2.
                                       ) * self._belowSidesIndicator, fenZERO]
 
     def buildA(self):
         """Build terms of operator of linear system."""
         if self.As[1] is None or self.As[2] is None or self.As[5] is None:
             thy2 = [(1.,), ('x', '()', 1, '*', 4., '-', 2. * self._H),
                     ('x', '()', 1, '**', 2., '*', 6., '-',
                      ('x', '()', 1, '*', 6. * self._H), '+', self._H ** 2.),
                     ('x', '()', 1, '**', 2., '*', 4., '-',
                      ('x', '()', 1, '*', 3. * self._H), '+', self._H ** 2.,
                      '*', ('x', '()', 1), '*', 2.),
                     ('x', '()', 1, '*', -1., '+', self._H, '*', ('x', '()', 1),
                      '**', 2.)]
         if self.As[3] is None or self.As[4] is None or self.As[5] is None:
             thz2 = [(1.,), ('x', '()', 2, '*', 4., '-', 2. * self._W),
                     ('x', '()', 2, '**', 2., '*', 6., '-',
                      ('x', '()', 2, '*', 6. * self._W), '+', self._W ** 2.),
                     ('x', '()', 2, '**', 2., '*', 4., '-',
                      ('x', '()', 2, '*', 3. * self._W), '+', self._W ** 2.,
                      '*', ('x', '()', 2), '*', 2.),
                     ('x', '()', 2, '*', -1., '+', self._W, '*', ('x', '()', 2),
                      '**', 2.)]
         if self.As[0] is None:
             vbMng(self, "INIT", "Assembling operator term A0.", 20)
             DirichletBC0 = fen.DirichletBC(self.V, fenZERO,
                                            self.DirichletBoundary)
             a0Re = fenZERO * fen.dot(self.u, self.v) * fen.dx
             self.As[0] = fenics2Sparse(a0Re, {}, DirichletBC0, 1)
             self.thAs[0] = self.getMonomialSingleWeight([0])
             vbMng(self, "DEL", "Done assembling operator term.", 20)
         if self.As[1] is None:
             vbMng(self, "INIT", "Assembling operator term A1.", 20)
             DirichletBC0 = fen.DirichletBC(self.V, fenZERO,
                                            self.DirichletBoundary)
             ajRe = ((self._aboveIndicator + self._belowSidesIndicator)
                   * 16. * self._L ** 2. * self.u.dx(0) * self.v.dx(0) * fen.dx)
             self.As[1] = fenics2Sparse(ajRe, {}, DirichletBC0, 0)
             thz0 = [(1.,), ('x', '()', 2, '*', 2.), ('x', '()', 2, '**', 2.)]
             self.thAs[1] = []
             idxs = nextDerivativeIndices([], 3, hashD([0, 4, 2]) + 1)
             for idx in idxs:
                 dy, dz = 4 - idx[1], 2 - idx[2]
                 if idx[0] == 0 and dy >= 0 and dz >= 0:
-                    self.thAs[1] += [thy2[dy] + ('*') + (thz0[dz],)]
+                    self.thAs[1] += [thy2[dy] + ('*',) + (thz0[dz],)]
                 else:
                     self.thAs[1] += [(0.,)]
             self.thAs[1] += [None]
             vbMng(self, "DEL", "Done assembling operator term.", 20)
         if self.As[2] is None:
             vbMng(self, "INIT", "Assembling operator term A2.", 20)
             DirichletBC0 = fen.DirichletBC(self.V, fenZERO,
                                            self.DirichletBoundary)
             ajRe = (4. * self._delta ** 2. * self._belowCenterIndicator
                   * self.u.dx(0) * self.v.dx(0) * fen.dx)
             self.As[2] = fenics2Sparse(ajRe, {}, DirichletBC0, 0)
             thz1 = [(1.,), ('x', '()', 2, '*', 2., '-', 2. * self._W),
                     ('x', '()', 2, '*', -1., '+', self._W, '**', 2.)]
             self.thAs[2] = []
             idxs = nextDerivativeIndices([], 3, hashD([0, 4, 2]) + 1)
             for idx in idxs:
                 dy, dz = 4 - idx[1], 2 - idx[2]
                 if idx[0] == 0 and dy >= 0 and dz >= 0:
-                    self.thAs[2] += [thy2[dy] + ('*') + (thz1[dz],)]
+                    self.thAs[2] += [thy2[dy] + ('*',) + (thz1[dz],)]
                 else:
                     self.thAs[2] += [(0.,)]
             self.thAs[2] += [None]
             vbMng(self, "DEL", "Done assembling operator term.", 20)
         if self.As[3] is None:
             vbMng(self, "INIT", "Assembling operator term A3.", 20)
             DirichletBC0 = fen.DirichletBC(self.V, fenZERO,
                                            self.DirichletBoundary)
             ajRe = (self._H ** 2. * self._aboveIndicator
                   * self.u.dx(1) * self.v.dx(1) * fen.dx)
             self.As[3] = fenics2Sparse(ajRe, {}, DirichletBC0, 0)
             thy1 = [(1.,), ('x', '()', 1, '*', 2., '-', 2. * self._H),
                     (self._H, '-', ('x', '()', 1), '**', 2.)]
             self.thAs[3] = []
             idxs = nextDerivativeIndices([], 3, hashD([0, 2, 4]) + 1)
             for idx in idxs:
                 dy, dz = 2 - idx[1], 4 - idx[2]
                 if idx[0] == 0 and dy >= 0 and dz >= 0:
-                    self.thAs[3] += [thy1[dy] + ('*') + (thz2[dz],)]
+                    self.thAs[3] += [thy1[dy] + ('*',) + (thz2[dz],)]
                 else:
                     self.thAs[3] += [(0.,)]
             self.thAs[3] += [None]
             vbMng(self, "DEL", "Done assembling operator term.", 20)
         if self.As[4] is None:
             vbMng(self, "INIT", "Assembling operator term A4.", 20)
             DirichletBC0 = fen.DirichletBC(self.V, fenZERO,
                                            self.DirichletBoundary)
             ajRe = ((self._belowSidesIndicator + self._belowCenterIndicator)
                   * self._H ** 2. * self.u.dx(1) * self.v.dx(1) * fen.dx)
             self.As[4] = fenics2Sparse(ajRe, {}, DirichletBC0, 0)
             thy0 = [(1.,), ('x', '()', 1, '*', 2.), ('x', '()', 1, '**', 2.)]
             self.thAs[4] = []
             idxs = nextDerivativeIndices([], 3, hashD([0, 2, 4]) + 1)
             for idx in idxs:
                 dy, dz = 2 - idx[1], 4 - idx[2]
                 if idx[0] == 0 and dy >= 0 and dz >= 0:
-                    self.thAs[4] += [thy0[dy] + ('*') + (thz2[dz],)]
+                    self.thAs[4] += [thy0[dy] + ('*',) + (thz2[dz],)]
                 else:
                     self.thAs[4] += [(0.,)]
             self.thAs[4] += [None]
             vbMng(self, "DEL", "Done assembling operator term.", 20)
         if self.As[5] is None:
             vbMng(self, "INIT", "Assembling operator term A5.", 20)
             DirichletBC0 = fen.DirichletBC(self.V, fenZERO,
                                            self.DirichletBoundary)
             ajRe = - 4. * self.u * self.v * fen.dx
             self.As[5] = fenics2Sparse(ajRe, {}, DirichletBC0, 0)
             thx = [(1.,), ('x', '()', 0)]
             self.thAs[5] = []
             idxs = nextDerivativeIndices([], 3, hashD([1, 4, 4]) + 1)
             for idx in idxs:
                 dx, dy, dz = 1 - idx[0], 4 - idx[1], 4 - idx[2]
                 if dx >= 0 and dy >= 0 and dz >= 0:
-                    self.thAs[5] += [thx[dx] + ('*') + (thy2[dy],)
-                                             + ('*') + (thz2[dz],)]
+                    self.thAs[5] += [thx[dx] + ('*',) + (thy2[dy],)
+                                             + ('*',) + (thz2[dz],)]
                 else:
                     self.thAs[5] += [(0.,)]
             self.thAs[5] += [None]
             vbMng(self, "DEL", "Done assembling operator term.", 20)
diff --git a/rrompy/hfengines/base/marginal_proxy_engine.py b/rrompy/hfengines/base/marginal_proxy_engine.py
index 2837328..4e4b762 100644
--- a/rrompy/hfengines/base/marginal_proxy_engine.py
+++ b/rrompy/hfengines/base/marginal_proxy_engine.py
@@ -1,138 +1,143 @@
 # 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 inspect
 import numpy as np
 from copy import copy as softcopy
 from rrompy.utilities.base.types import Np1D, paramVal, paramList, HFEng
 from rrompy.utilities.base import freepar as fp
 from rrompy.utilities.exception_manager import RROMPyException
 from rrompy.parameter import checkParameter, checkParameterList
 
 __all__ = ['MarginalProxyEngine']
 
 class MarginalProxyEngine:
     """
     Marginalized should prescribe fixed value for the marginalized parameters
         and leave freepar/None elsewhere.
     """
     
     _allowedMuDependencies = ["A", "b", "checkParameter", "checkParameterList",
                               "_assembleObject", "solve", "residual",
-                              "stabilityFactor"]
+                              "stabilityFactor", "stateShift", "outputShift"]
     
     def __init__(self, HFEngine:HFEng, marginalized:Np1D):
         self.baseHF = HFEngine
         self.marg = marginalized
         for name in HFEngine.__dir_base__():
             att = getattr(HFEngine, name)
             if inspect.ismethod(att):
                 attargs = inspect.getfullargspec(att).args
                 if "mu" not in attargs:
                     setattr(self.__class__, name, getattr(HFEngine, name))
                 else:
                     if name not in self._allowedMuDependencies:
                         raise RROMPyException(("Function {} depends on mu "
                                                "and was not accounted for. "
                                                "Must override.").format(name))
     
     @property
     def affinePoly(self):
         return self.nparFixed == 0 and self.baseHF.affinePoly
 
     @property
     def freeLocations(self):
-        return tuple([x for x in range(self.baseHF.npar) \
-                                                        if self.marg[x] == fp])
+        return [x for x in range(self.baseHF.npar) if self.marg[x] == fp]
 
     @property
     def fixedLocations(self):
-        return tuple([x for x in range(self.baseHF.npar) \
-                                                        if self.marg[x] != fp])
+        return [x for x in range(self.baseHF.npar) if self.marg[x] != fp]
 
     @property
     def _freeLocationsInsert(self):
         return np.cumsum([m == fp for m in self.marg])[self.fixedLocations]
 
     @property
     def muFixed(self):
         muF = checkParameter([m for m in self.marg if m != fp])
         if self.baseHF.npar - self.nparFree > 0: muF = muF[0]
         return muF
 
     @property
     def nparFree(self):
         """Value of nparFree."""
         return len(self.freeLocations)
 
     @property
     def nparFixed(self):
         """Value of nparFixed."""
         return len(self.fixedLocations)
 
     def name(self) -> str:
         return "{}-proxy for {}".format(self.freeLocations, self.baseHF.name())
         
     def __str__(self) -> str:
         return self.name()
 
     def __repr__(self) -> str:
         return self.__str__() + " at " + hex(id(self))
 
     def __dir_base__(self):
         return [x for x in self.__dir__() if x[:2] != "__"]
 
     def __deepcopy__(self, memo):
         return softcopy(self)
 
     def completeMu(self, mu:paramVal) -> paramVal:
         mu = checkParameter(mu, self.nparFree)
         return np.insert(mu.data, self._freeLocationsInsert, self.muFixed,
                          axis = 1)
         
     def completeMuList(self, mu:paramList) -> paramList:
         mu = checkParameterList(mu, self.nparFree)[0]
         return np.insert(mu.data, self._freeLocationsInsert, self.muFixed,
                          axis = 1)
 
     def A(self, mu : paramVal = [], *args, **kwargs):
         return self.baseHF.A(self.completeMu(mu), *args, **kwargs)
 
     def b(self, mu : paramVal = [], *args, **kwargs):
         return self.baseHF.b(self.completeMu(mu), *args, **kwargs)
 
     def checkParameter(self, mu : paramVal = [], *args, **kwargs):
         return self.baseHF.checkParameter(self.completeMu(mu), *args, **kwargs)
 
     def checkParameterList(self, mu : paramList = [], *args, **kwargs):
         return self.baseHF.checkParameterList(self.completeMuList(mu),
                                               *args, **kwargs)
 
     def _assembleObject(self, mu : paramVal = [], *args, **kwargs):
         return self.baseHF._assembleObject(self.completeMu(mu),
                                            *args, **kwargs)
 
+    def stateShift(self, mu : paramVal = [], *args, **kwargs):
+        return self.baseHF.stateShift(self.completeMuList(mu), *args, **kwargs)
+
+    def outputShift(self, mu : paramVal = [], *args, **kwargs):
+        return self.baseHF.outputShift(self.completeMuList(mu),
+                                       *args, **kwargs)
+
     def solve(self, mu : paramList = [], *args, **kwargs):
         return self.baseHF.solve(self.completeMuList(mu), *args, **kwargs)
 
     def residual(self, mu : paramList = [], *args, **kwargs):
         return self.baseHF.residual(self.completeMuList(mu), *args, **kwargs)
 
     def stabilityFactor(self, mu : paramList = [], *args, **kwargs):
         return self.baseHF.stabilityFactor(self.completeMuList(mu),
                                            *args, **kwargs)
diff --git a/rrompy/sampling/pivoted/sampling_engine_pivoted.py b/rrompy/sampling/pivoted/sampling_engine_pivoted.py
index 343fa3c..58d539d 100644
--- a/rrompy/sampling/pivoted/sampling_engine_pivoted.py
+++ b/rrompy/sampling/pivoted/sampling_engine_pivoted.py
@@ -1,131 +1,131 @@
 # 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 rrompy.sampling.base.sampling_engine_base_pivoted import (
                                                      SamplingEngineBasePivoted)
 from rrompy.hfengines.base import MarginalProxyEngine
 from rrompy.utilities.base.types import Np1D, paramVal, paramList, sampList
 from rrompy.utilities.base import verbosityManager as vbMng, freepar as fp
 from rrompy.utilities.exception_manager import RROMPyException
 from rrompy.utilities.numerical import nextDerivativeIndices, dot
 from rrompy.parameter import checkParameter, checkParameterList
 from rrompy.sampling import sampleList
 
 __all__ = ['SamplingEnginePivoted']
 
 class SamplingEnginePivoted(SamplingEngineBasePivoted):
     """HERE"""
 
     def preprocesssamples(self, idxs:Np1D, j:int) -> sampList:
         if self.samples[j] is None or len(self.samples[j]) == 0: return
         return self.samples[j](idxs)
 
     def setsample(self, u:sampList, j:int, overwrite : bool = False) -> Np1D:
         if overwrite:
             self.samples[j][self.nsamples[j]] = u
         else:
             if self.nsamples[j] == 0:
                 self.samples[j] = sampleList(u)
             else:
                 self.samples[j].append(u)
 
     def postprocessu(self, u:sampList, j:int, overwrite : bool = False):
         self.setsample(u, j, overwrite)
 
     def postprocessuBulk(self, j:int):
         pass
 
     def _getSampleConcurrence(self, mu:paramVal, j:int,
                               previous:Np1D) -> sampList:
         if not (self.sample_state or self.HFEngine.isCEye):
             raise RROMPyException(("Derivatives of solution with non-scalar "
                                    "C not computable."))
         if not self.HFEngine._isStateShiftZero:
             raise RROMPyException(("Derivatives of solution with non-zero "
                                    "solution shift not computable."))
         if len(previous) >= len(self._derIdxs[j]):
             self._derIdxs[j] += nextDerivativeIndices(
                                      self._derIdxs[j], self.nPivot,
                                      len(previous) + 1 - len(self._derIdxs[j]))
         derIdx = self._derIdxs[j][len(previous)]
         mu = checkParameter(mu, self.nPivot)
         samplesOld = self.preprocesssamples(previous, j)
         RHS = self.HFEngineMarginalized.b(mu, derIdx)
         for j, derP in enumerate(self._derIdxs[j][: len(previous)]):
             diffP = [x - y for (x, y) in zip(derIdx, derP)]
             if np.all([x >= 0 for x in diffP]):
                 RHS -= dot(self.HFEngineMarginalized.A(mu, diffP),
                            samplesOld[j])
         return self.solveLS(mu, RHS = RHS)
 
     def nextSample(self, mu:paramVal, j:int, overwrite : bool = False,
                    postprocess : bool = True) -> Np1D:
         mu = checkParameter(mu, self.nPivot)
         muidxs = self.mus[j].findall(mu[0])
         if len(muidxs) > 0:
             u = self._getSampleConcurrence(mu, j, np.sort(muidxs))
         else:
             u = self.solveLS(mu)
         if postprocess:
             self.postprocessu(u, j, overwrite = overwrite)
         else:
             self.setsample(u, j, overwrite)
         if overwrite:
             self.mus[j][self.nsamples[j]] = mu[0]
         else:
             self.mus[j].append(mu)
         self.nsamples[j] += 1
         return self.samples[j][self.nsamples[j] - 1]
 
     def iterSample(self, mus:paramList, musM:paramList) -> sampList:
         mus = checkParameterList(mus, self.nPivot)[0]
         musM = checkParameterList(musM, self.nMarginal)[0]
         vbMng(self, "INIT", "Starting sampling iterations.", 5)
         n = len(mus)
         m = len(musM)
         if n <= 0:
             raise RROMPyException("Number of samples must be positive.")
         if m <= 0:
             raise RROMPyException(("Number of marginal samples must be "
                                    "positive."))
         repeatedSamples = len(mus.unique()) != n
         for j in range(m):
             muMEff = [fp] * self.HFEngine.npar
             for k, x in enumerate(self.directionMarginal):
                 muMEff[x] = musM(j, k)
             self.HFEngineMarginalized = MarginalProxyEngine(self.HFEngine,
                                                             list(muMEff))
             if repeatedSamples:
                 for k in range(n):
                     vbMng(self, "MAIN",
                           "Computing sample {} / {} for marginal {} / {}."\
                                                    .format(k + 1, n, j, m), 10)
                     self.nextSample(mus[k], j, overwrite = (k > 0),
                                     postprocess = False)
                     if n > 1 and k == 0:
                         self.preallocateSamples(self.samples[j][0], mus[0],
                                                 n, j)
             else:
-                self.samples[j] = self.postprocessuBulk(self.solveLS(mus), j)
+                self.setsample(self.solveLS(mus), j, overwrite = False)
                 self.mus[j] = copy(mus)
                 self.nsamples[j] = n
             self.postprocessuBulk(j)
         vbMng(self, "DEL", "Finished sampling iterations.", 5)
         return self.samples[j]