diff --git a/examples/1_symmetric_disk/symmetric_disk.py b/examples/1_symmetric_disk/symmetric_disk.py
index 8e4b003..285e395 100644
--- a/examples/1_symmetric_disk/symmetric_disk.py
+++ b/examples/1_symmetric_disk/symmetric_disk.py
@@ -1,87 +1,87 @@
 import numpy as np
 from symmetric_disk_engine import SymmetricDiskEngine as engine
 from rrompy.sampling.standard import SamplingEngineStandard as SES
 from rrompy.reduction_methods.standard import (NearestNeighbor as NN,
                                                RationalInterpolant as RI,
                                                ReducedBasis as RB)
 from rrompy.reduction_methods.standard.greedy import (
                                               RationalInterpolantGreedy as RIG,
                                               ReducedBasisGreedy as RBG)
 from rrompy.parameter.parameter_sampling import QuadratureSampler as QS
 
 ks = [10., 20.]
 k0, n = np.mean(np.power(ks, 2.)) ** .5, 150
 solver = engine(k0, n)
 k = 12.
 
 for method in ["RI", "RB", "RI_GREEDY", "RB_GREEDY"]:
     print("Testing {} method".format(method))
     
     if method == "RI":
-        params = {'N':39, 'M':39, 'S':40, 'POD':True, 'polybasis':"CHEBYSHEV",
+        params = {'S':40, 'POD':True, 'polybasis':"CHEBYSHEV",
                   'sampler':QS(ks, "CHEBYSHEV", scalingExp = 2.)}
         algo = RI
     if method == "RB":
-        params = {'R':40, 'S':40, 'POD':True,
+        params = {'S':40, 'POD':True,
                   'sampler':QS(ks, "CHEBYSHEV", scalingExp = 2.)}
         algo = RB
     if method == "RI_GREEDY":
         params = {'S':10, 'POD':True, 'polybasis':"LEGENDRE", 'greedyTol':1e-2,
                   'sampler':QS(ks, "UNIFORM", scalingExp = 2.),
                   'errorEstimatorKind':"DISCREPANCY"}
         algo = RIG
     if method == "RB_GREEDY":
         params = {'S':10, 'POD':True, 'greedyTol':1e-2,
                   'sampler':QS(ks, "UNIFORM", scalingExp = 2.)}
         algo = RBG
         
     approx = algo(solver, mu0 = k0, approx_state = True,
                   approxParameters = params, verbosity = 20)
     if len(method) == 2:
         approx.setupApprox()
     else:
         approx.setupApprox("LAST")
     
     print("--- Approximant ---")
     approx.plotApprox(k, name = 'u_app')
     approx.plotHF(k, name = 'u_HF')
     approx.plotErr(k, name = 'err_app')
     approx.plotRes(k, name = 'res_app')
     normErr = approx.normErr(k)[0]
     normSol = approx.normHF(k)[0]
     normRes = approx.normRes(k)[0]
     normRHS = approx.normRHS(k)[0]
     print("SolNorm:\t{:.5e}\nErr_app: \t{:.5e}\nErrRel_app:\t{:.5e}".format(
                                           normSol, normErr, normErr / normSol))
     print("RHSNorm:\t{:.5e}\nRes_app: \t{:.5e}\nResRel_app:\t{:.5e}".format(
                                           normRHS, normRes, normRes / normRHS))
     
     print("--- Closest snapshot ---")
     eng = SES(solver, verbosity = 0)
     approxNN = NN(solver, mu0 = k0, approx_state = True, verbosity = 0,
                   approxParameters = {'S':approx.S, 'POD':True})
     approxNN.setSamples(approx.samplingEngine)
     approxNN.plotApprox(k, name = 'u_close')
     approxNN.plotHF(k, name = 'u_HF')
     approxNN.plotErr(k, name = 'err_close')
     approxNN.plotRes(k, name = 'res_close')
     normErr = approxNN.normErr(k)[0]
     normSol = approxNN.normHF(k)[0]
     normRes = approxNN.normRes(k)[0]
     normRHS = approxNN.normRHS(k)[0]
     print("SolNorm:\t{:.5e}\nErr_close:\t{:.5e}\nErrRel_close:\t{:.5e}".format(
                                           normSol, normErr, normErr / normSol))
     print("RHSNorm:\t{:.5e}\nRes_close:\t{:.5e}\nResRel_close:\t{:.5e}".format(
                                           normRHS, normRes, normRes / normRHS))
 
     if method[:2] == "RI":
         poles, residues = approx.getResidues()
     if method[:2] == "RB":
         poles = approx.getPoles()
     print("Poles:\n{}".format(poles))
     if method[:2] == "RI":
         for pol, res in zip(poles, residues.T):
             solver.plot(res)
             print("pole = {:.5e}".format(pol))
 
     print("\n")
diff --git a/examples/2_double_slit/double_slit.py b/examples/2_double_slit/double_slit.py
index 283c717..bf0c84b 100644
--- a/examples/2_double_slit/double_slit.py
+++ b/examples/2_double_slit/double_slit.py
@@ -1,84 +1,84 @@
 import numpy as np
 from double_slit_engine import DoubleSlitEngine as engine
 from rrompy.sampling.standard import SamplingEngineStandard as SES
 from rrompy.reduction_methods.standard import (NearestNeighbor as NN,
                                                RationalInterpolant as RI)
 from rrompy.reduction_methods.standard.greedy import (
                                               RationalInterpolantGreedy as RIG)
 from rrompy.parameter.parameter_sampling import QuadratureSampler as QS
 from rrompy.solver.fenics import interp_project
 
 ks = [10., 15.]
 k0, n = np.mean(ks), 150
 solver = engine(k0, n)
 k = 11.
 
 for method in ["RI", "RI_GREEDY"]:
     print("Testing {} method".format(method))
 
     if method == "RI":
-        params = {'N':19, 'M':19, 'S':20, 'POD':True, 'polybasis':"CHEBYSHEV",
+        params = {'S':20, 'POD':True, 'polybasis':"CHEBYSHEV",
                   'sampler':QS(ks, "CHEBYSHEV")}
         algo = RI
     if method == "RI_GREEDY":
         params = {'S':10, 'POD':True, 'polybasis':"LEGENDRE", 'greedyTol':1e-2,
                   'sampler':QS(ks, "UNIFORM"),
                   'errorEstimatorKind':"LOOK_AHEAD",
                   'trainSetGenerator':QS(ks, "CHEBYSHEV")}
         algo = RIG
         
     approx = algo(solver, mu0 = k0, approx_state = True,
                   approxParameters = params, verbosity = 20)
     if len(method) == 2:
         approx.setupApprox()
     else:
         approx.setupApprox("LAST")
     
     print("--- Approximant ---")
     approx.plotApprox(k, name = 'u_app')
     approx.plotHF(k, name = 'u_HF')
     approx.plotErr(k, name = 'err_app')
     approx.plotRes(k, name = 'res_app')
     normErr = approx.normErr(k)[0]
     normSol = approx.normHF(k)[0]
     normRes = approx.normRes(k)[0]
     normRHS = approx.normRHS(k)[0]
     print("SolNorm:\t{:.5e}\nErr_app: \t{:.5e}\nErrRel_app:\t{:.5e}".format(
                                           normSol, normErr, normErr / normSol))
     print("RHSNorm:\t{:.5e}\nRes_app: \t{:.5e}\nResRel_app:\t{:.5e}".format(
                                           normRHS, normRes, normRes / normRHS))
     
     print("--- Closest snapshot ---")
     eng = SES(solver, verbosity = 0)
     approxNN = NN(solver, mu0 = k0, approx_state = True, verbosity = 0,
                   approxParameters = {'S':approx.S, 'POD':True})
     approxNN.setSamples(approx.samplingEngine)
     approxNN.plotApprox(k, name = 'u_close')
     approxNN.plotHF(k, name = 'u_HF')
     approxNN.plotErr(k, name = 'err_close')
     approxNN.plotRes(k, name = 'res_close')
     normErr = approxNN.normErr(k)[0]
     normSol = approxNN.normHF(k)[0]
     normRes = approxNN.normRes(k)[0]
     normRHS = approxNN.normRHS(k)[0]
     print("SolNorm:\t{:.5e}\nErr_close:\t{:.5e}\nErrRel_close:\t{:.5e}".format(
                                           normSol, normErr, normErr / normSol))
     print("RHSNorm:\t{:.5e}\nRes_close:\t{:.5e}\nResRel_close:\t{:.5e}".format(
                                           normRHS, normRes, normRes / normRHS))
 
     uIncR, uIncI = solver.getDirichletValues(k)
     uIncR = interp_project(uIncR, solver.V)
     uIncI = interp_project(uIncI, solver.V)
     uInc = np.array(uIncR.vector()) + 1.j * np.array(uIncI.vector())
     uEx = approx.getHF(k)[0] - uInc
     uApp = approx.getApprox(k)[0] - uInc
     solver.plot(uEx, name = 'uex_tot')
     solver.plot(uApp, name = 'u_app_tot')
     
     poles, residues = approx.getResidues()
     print("Poles:\n{}".format(poles))
     for pol, res in zip(poles, residues.T):
         solver.plot(res)
         print("pole = {:.5e}".format(pol))
 
     print("\n")
diff --git a/examples/3_sector_angle/sector_angle.py b/examples/3_sector_angle/sector_angle.py
index eb73ecf..da0776b 100644
--- a/examples/3_sector_angle/sector_angle.py
+++ b/examples/3_sector_angle/sector_angle.py
@@ -1,114 +1,114 @@
 import numpy as np
 import matplotlib.pyplot as plt
 from sector_angle_engine import SectorAngleEngine as engine
 from rrompy.sampling.standard import SamplingEngineStandard as SES
 from rrompy.reduction_methods.standard import NearestNeighbor as NN
 from rrompy.reduction_methods.pivoted import (
                                       RationalInterpolantPivoted as RIP,
                                       RationalInterpolantGreedyPivoted as RIGP)
 from rrompy.parameter.parameter_sampling import QuadratureSampler as QS
 
 ks, ts = [10., 15.], [.4, .6]
 k0, t0, n = np.mean(np.power(ks, 2.)) ** .5, np.mean(ts), 50
 solver = engine(k0, t0, n)
 
 murange = [[ks[0], ts[0]], [ks[-1], ts[-1]]]
 mu = [12., .535]
 
 fighandles = []
 for method in ["RI", "RI_GREEDY"]:
     print("Testing {} method".format(method))
 
     if method == "RI":
-        params = {'N':19, 'M':19, 'S':20, 'MMarginal':3, 'SMarginal':11,
+        params = {'S':20, 'MMarginal':3, 'SMarginal':11,
                   'POD':True, 'polybasis':"CHEBYSHEV",
                   'polybasisMarginal':"MONOMIAL_GAUSSIAN",
                   'radialDirectionalWeightsMarginal': 10.,
                   'matchingWeight':1., 'cutOffTolerance': 2.,
                   'samplerPivot':QS(ks, "CHEBYSHEV", 2.),
                   'samplerMarginal':QS(ts, "UNIFORM")}
         algo = RIP
     if method == "RI_GREEDY":
         params = {'S':10, 'MMarginal':3, 'SMarginal':11, 'POD':True,
                   'polybasis':"LEGENDRE",
                   'polybasisMarginal':"MONOMIAL_GAUSSIAN",
                   'radialDirectionalWeightsMarginal': 10.,
                   'matchingWeight':1., 'cutOffTolerance': 2.,
                   'samplerPivot':QS(ks, "UNIFORM", 2.),
                   'greedyTol':1e-3, 'errorEstimatorKind':"INTERPOLATORY",
                   'trainSetGenerator':QS(ks, "CHEBYSHEV", 2.),
                   'samplerMarginal':QS(ts, "UNIFORM")}
         algo = RIGP
         
     approx = algo([0], solver, mu0 = [k0, t0], approx_state = True,
                   approxParameters = params, verbosity = 10)
     if len(method) == 2:
         approx.setupApprox()
     else:
         approx.setupApprox("LAST")
     
     print("--- Approximant ---")
     approx.plotApprox(mu, name = 'u_app')
     approx.plotHF(mu, name = 'u_HF')
     approx.plotErr(mu, name = 'err_app')
     approx.plotRes(mu, name = 'res_app')
     normErr = approx.normErr(mu)[0]
     normSol = approx.normHF(mu)[0]
     normRes = approx.normRes(mu)[0]
     normRHS = approx.normRHS(mu)[0]
     print("SolNorm:\t{:.5e}\nErr_app: \t{:.5e}\nErrRel_app:\t{:.5e}".format(
                                           normSol, normErr, normErr / normSol))
     print("RHSNorm:\t{:.5e}\nRes_app: \t{:.5e}\nResRel_app:\t{:.5e}".format(
                                           normRHS, normRes, normRes / normRHS))
     
     print("--- Closest snapshot ---")
     eng = SES(solver, verbosity = 0)
     eng.nsamples = approx.samplingEngine.nsamplesCoalesced
     eng.mus = approx.samplingEngine.musCoalesced
     eng.samples = approx.samplingEngine.samples_fullCoalesced
     paramsNN = {'S':eng.nsamples, 'POD':False,
                 'sampler':QS([[ks[0], ts[0]], [ks[-1], ts[-1]]], "UNIFORM")}
     approxNN = NN(solver, mu0 = [k0, t0], approx_state = True,
                   approxParameters = paramsNN, verbosity = 0)
     approxNN.setSamples(eng)
     approxNN.plotApprox(mu, name = 'u_close')
     approxNN.plotHF(mu, name = 'u_HF')
     approxNN.plotErr(mu, name = 'err_close')
     approxNN.plotRes(mu, name = 'res_close')
     normErr = approxNN.normErr(mu)[0]
     normSol = approxNN.normHF(mu)[0]
     normRes = approxNN.normRes(mu)[0]
     normRHS = approxNN.normRHS(mu)[0]
     print("SolNorm:\t{:.5e}\nErr_close:\t{:.5e}\nErrRel_close:\t{:.5e}".format(
                                           normSol, normErr, normErr / normSol))
     print("RHSNorm:\t{:.5e}\nRes_close:\t{:.5e}\nResRel_close:\t{:.5e}".format(
                                           normRHS, normRes, normRes / normRHS))
 
     verb = approx.verbosity
     approx.verbosity = 0
     tspace = np.linspace(ts[0], ts[-1], 100)
     for j, t in enumerate(tspace):
         pls = approx.getPoles([None, t])
         pls[np.abs(np.imag(pls ** 2.)) > 1e-5] = np.nan
         if j == 0: poles = np.empty((len(tspace), len(pls)))
         poles[j] = np.real(pls)
     approx.verbosity = verb
     fighandles += [plt.figure(figsize = (12, 5))]
     ax1 = fighandles[-1].add_subplot(1, 2, 1)
     ax2 = fighandles[-1].add_subplot(1, 2, 2)
     ax1.plot(poles, tspace)
     ax1.set_ylim(ts)
     ax1.set_xlabel('mu_1')
     ax1.set_ylabel('mu_2')
     ax1.grid()
     ax2.plot(poles, tspace)
     for mm in approx.musMarginal:
         ax2.plot(ks, [mm[0, 0]] * 2, 'k--', linewidth = 1)
     ax2.set_xlim(ks)
     ax2.set_ylim(ts)
     ax2.set_xlabel('mu_1')
     ax2.set_ylabel('mu_2')
     ax2.grid()
     plt.show()
     
     print("\n")
diff --git a/examples/4_funnel_output/funnel_output.py b/examples/4_funnel_output/funnel_output.py
index e618ce6..bfff8df 100644
--- a/examples/4_funnel_output/funnel_output.py
+++ b/examples/4_funnel_output/funnel_output.py
@@ -1,59 +1,59 @@
 import numpy as np
 from funnel_output_engine import FunnelOutputEngine as engine
 from rrompy.sampling.standard import SamplingEngineStandard as SES
 from rrompy.reduction_methods.standard import (NearestNeighbor as NN,
                                                RationalInterpolant as RI)
 from rrompy.reduction_methods.standard.greedy import (
                                               RationalInterpolantGreedy as RIG)
 from rrompy.parameter.parameter_sampling import QuadratureSampler as QS
 
 ks = [5., 10.]
 k0, n = np.mean(ks), 50
 solver = engine(k0, n)
 k = 6.5
 
 for method in ["RI", "RI_OUTPUT", "RI_GREEDY", "RI_GREEDY_OUTPUT"]:
     print("Testing {} method".format(method))
     
     if "GREEDY" not in method:
-        params = {'N':19, 'M':19, 'S':20, 'POD':True, 'polybasis':"CHEBYSHEV",
+        params = {'S':20, 'POD':True, 'polybasis':"CHEBYSHEV",
                   'sampler':QS(ks, "CHEBYSHEV")}
         algo = RI
     if "GREEDY" in method:
         params = {'S':2, 'POD':True, 'polybasis':"LEGENDRE", 'greedyTol':1e-1,
                   'maxIter':25, 'sampler':QS(ks, "UNIFORM"),
                   'errorEstimatorKind':"LOOK_AHEAD_OUTPUT"}
         algo = RIG
         
     approx = algo(solver, mu0 = k0, approx_state = method[-7 :] != "_OUTPUT",
                   approxParameters = params, verbosity = 5)
     if "GREEDY" not in method:
         approx.setupApprox()
     else:
         approx.setupApprox("LAST")
     
     print("--- Approximant ---")
     approx.plotApprox(k, name = 'u_app')
     approx.plotHF(k, name = 'u_HF')
     approx.plotErr(k, name = 'err_app')
     normErr = approx.normErr(k)[0]
     normSol = approx.normHF(k)[0]
     print("SolNorm:\t{:.5e}\nErr_app: \t{:.5e}\nErrRel_app:\t{:.5e}".format(
                                           normSol, normErr, normErr / normSol))
     
     print("--- Closest snapshot ---")
     eng = SES(solver, verbosity = 0)
     approxNN = NN(solver, mu0 = k0, approx_state = method[-7 :] != "_OUTPUT",
                   approxParameters = {'S':approx.S, 'POD':True}, verbosity = 0)
     approxNN.setSamples(approx.samplingEngine)
     approxNN.plotApprox(k, name = 'u_close')
     approxNN.plotHF(k, name = 'u_HF')
     approxNN.plotErr(k, name = 'err_close')
     normErr = approxNN.normErr(k)[0]
     normSol = approxNN.normHF(k)[0]
     print("SolNorm:\t{:.5e}\nErr_close:\t{:.5e}\nErrRel_close:\t{:.5e}".format(
                                           normSol, normErr, normErr / normSol))
 
     print("Poles:\n{}".format(approx.getPoles()))
 
     print("\n")
diff --git a/tests/reduction_methods_1D/rational_interpolant_1d.py b/tests/reduction_methods_1D/rational_interpolant_1d.py
index 9fc1cf9..3e1b673 100644
--- a/tests/reduction_methods_1D/rational_interpolant_1d.py
+++ b/tests/reduction_methods_1D/rational_interpolant_1d.py
@@ -1,71 +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/>.
 #
 
 import numpy as np
 from matrix_fft import matrixFFT
 from rrompy.reduction_methods.standard import RationalInterpolant as RI
 from rrompy.parameter.parameter_sampling import (QuadratureSampler as QS,
                                                  ManualSampler as MS)
 from rrompy.parameter import checkParameterList
 
 def test_monomials(capsys):
     mu = 1.5
     solver = matrixFFT()
-    params = {"POD": False, "M": 9, "N": 9, "S": 10, "robustTol": 1e-6,
-              "interpRcond": 1e-3, "polybasis": "MONOMIAL",
-              "sampler": QS([1.5, 6.5], "UNIFORM")}
+    params = {"POD": False, "S": 10, "robustTol": 1e-6, "interpRcond": 1e-3,
+              "polybasis": "MONOMIAL", "sampler": QS([1.5, 6.5], "UNIFORM")}
     approx = RI(solver, 4., approx_state = True, approxParameters = params,
                 verbosity = 10)
     approx.setupApprox()
 
     out, err = capsys.readouterr()
     assert "poorly conditioned. Reducing N " in out
     assert len(err) == 0
     assert np.isclose(approx.normErr(mu)[0], 1.4746e-05, atol = 1e-4)
 
 def test_well_cond():
     mu = 1.5
     solver = matrixFFT()
-    params = {"POD": True, "M": 9, "N": 9, "S": 10, "robustTol": 1e-14,
-              "interpRcond": 1e-10, "polybasis": "CHEBYSHEV",
-              "sampler": QS([1., 7.], "CHEBYSHEV")}
+    params = {"POD": True, "S": 10, "robustTol": 1e-14, "interpRcond": 1e-10,
+              "polybasis": "CHEBYSHEV", "sampler": QS([1., 7.], "CHEBYSHEV")}
     approx = RI(solver, 4., approx_state = True, approxParameters = params,
                 verbosity = 0)
     approx.setupApprox()
     poles = approx.getPoles()
     for lambda_ in np.arange(1, 8):
         assert np.isclose(np.min(np.abs(poles - lambda_)), 0., atol = 1e-4)
     for mu in approx.mus:
         assert np.isclose(approx.normErr(mu)[0], 0., atol = 1e-8)
 
 def test_hermite():
     mu = 1.5
     solver = matrixFFT()
     sampler0 = QS([1., 7.], "CHEBYSHEV")
     points, _ = checkParameterList(np.tile(sampler0.generatePoints(4)(0), 3))
-    params = {"POD": True, "M": 11, "N": 11, "S": 12, "polybasis": "CHEBYSHEV",
+    params = {"POD": True, "S": 12, "polybasis": "CHEBYSHEV",
               "sampler": MS([1., 7.], points = points)}
     approx = RI(solver, 4., approx_state = True, approxParameters = params,
                 verbosity = 0)
     approx.setupApprox()
     poles = approx.getPoles()
     for lambda_ in np.arange(1, 8):
         assert np.isclose(np.min(np.abs(poles - lambda_)), 0., atol = 1e-4)
     for mu in approx.mus:
         assert np.isclose(approx.normErr(mu)[0], 0., atol = 1e-8)
 
diff --git a/tests/reduction_methods_multiD/greedy_pivoted_rational_2d.py b/tests/reduction_methods_multiD/greedy_pivoted_rational_2d.py
index dc9509e..6d07613 100644
--- a/tests/reduction_methods_multiD/greedy_pivoted_rational_2d.py
+++ b/tests/reduction_methods_multiD/greedy_pivoted_rational_2d.py
@@ -1,83 +1,83 @@
 # 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.pivoted.greedy import (
                                RationalInterpolantPivotedGreedy as RIPG,
                                RationalInterpolantGreedyPivotedGreedy as RIGPG)
 from rrompy.parameter.parameter_sampling import QuadratureSampler as QS
 from rrompy.parameter import localSparseGrid as LSG
 
 def test_greedy_pivoted():
     mu = [5.05, 7.1]
     mu0 = [5., 7.]
     solver = matrixRandom()
     uh = solver.solve(mu)[0]
-    params = {"POD": True, "M": 4, "N": 4, "S": 5, "polybasis": "CHEBYSHEV",
+    params = {"POD": True, "S": 5, "polybasis": "CHEBYSHEV",
               "samplerPivot": QS([4.75, 5.25], "CHEBYSHEV"),
               "MMarginal": 1, "SMarginal": 3, "greedyTolMarginal": 1e-2,
               "radialDirectionalWeightsMarginal": 2.,
               "polybasisMarginal": "MONOMIAL_GAUSSIAN",
               "matchingWeight": 1., "samplerMarginalGrid":LSG([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"),
               "MMarginal": 1, "SMarginal": 3, "greedyTolMarginal": 1e-2,
               "radialDirectionalWeightsMarginal": 2.,
               "polybasisMarginal": "MONOMIAL_GAUSSIAN",
               "matchingWeight": 1., "samplerMarginalGrid":LSG([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)
diff --git a/tests/reduction_methods_multiD/pivoted_rational_2d.py b/tests/reduction_methods_multiD/pivoted_rational_2d.py
index 9e24654..804a2d8 100644
--- a/tests/reduction_methods_multiD/pivoted_rational_2d.py
+++ b/tests/reduction_methods_multiD/pivoted_rational_2d.py
@@ -1,114 +1,112 @@
 # Copyright (C) 2018 by the RROMPy authors
 #
 # This file is part of RROMPy.
 #
 # RROMPy is free software: you can redistribute it and/or modify
 # it under the terms of the GNU Lesser General Public License as published by
 # the Free Software Foundation, either version 3 of the License, or
 # (at your option) any later version.
 #
 # RROMPy is distributed in the hope that it will be useful,
 # but WITHOUT ANY WARRANTY; without even the implied warranty of
 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
 # GNU Lesser General Public License for more details.
 #
 # You should have received a copy of the GNU Lesser General Public License
 # along with RROMPy. If not, see <http://www.gnu.org/licenses/>.
 #
 
 import numpy as np
 from matrix_random import matrixRandom
 from rrompy.reduction_methods.pivoted import (
                                       RationalInterpolantPivoted as RIP,
                                       RationalInterpolantGreedyPivoted as RIGP)
 from rrompy.parameter.parameter_sampling import (QuadratureSampler as QS,
                                                  ManualSampler as MS)
 
 def test_pivoted_uniform():
     mu = [5.05, 7.1]
     mu0 = [5., 7.]
     solver = matrixRandom()
     uh = solver.solve(mu)[0]
-    params = {"POD": True, "M": 4, "N": 4, "S": 5, "polybasis": "CHEBYSHEV",
-              "samplerPivot": QS([4.75, 5.25], "CHEBYSHEV"),
-              "MMarginal": 4, "SMarginal": 5, "polybasisMarginal": "MONOMIAL",
-              "samplerMarginal": QS([6.75, 7.25], "UNIFORM"),
-              "matchingWeight": 1.}
+    params = {"POD": True, "S": 5, "polybasis": "CHEBYSHEV",
+              "samplerPivot": QS([4.75, 5.25], "CHEBYSHEV"), "SMarginal": 5,
+              "polybasisMarginal": "MONOMIAL", "matchingWeight": 1.,
+              "samplerMarginal": QS([6.75, 7.25], "UNIFORM")}
     approx = RIP([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_pivoted_manual_grid(capsys):
     mu = [5.05, 7.1]
     mu0 = [5., 7.]
     solver = matrixRandom()
     uh = solver.solve(mu)[0]
-    params = {"POD": False, "M": 4, "N": 4, "S": 5, "polybasis": "MONOMIAL",
-              "samplerPivot": MS([4.75, 5.25], np.array([5.])),
-              "MMarginal": 4, "SMarginal": 5, "polybasisMarginal": "MONOMIAL",
+    params = {"POD": False, "S": 5, "polybasis": "MONOMIAL",
+              "samplerPivot": MS([4.75, 5.25], np.array([5.])), "SMarginal": 5,
+              "polybasisMarginal": "MONOMIAL", "matchingWeight": 1.,
               "samplerMarginal": MS([6.75, 7.25], np.linspace(6.75, 7.25, 5)),
-              "matchingWeight": 1., "robustTol": 1e-6, "interpRcond": 1e-3,
-              "cutOffTolerance": 1.}
+              "robustTol": 1e-6, "interpRcond": 1e-3, "cutOffTolerance": 1.}
     approx = RIP([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), .4763489, rtol = 1)
 
     out, err = capsys.readouterr()
     assert ("poorly conditioned" not in out)
     assert len(err) == 0
 
 def test_pivoted_greedy():
     mu = [5.05, 7.1]
     mu0 = [5., 7.]
     solver = matrixRandom()
     uh = solver.solve(mu)[0]
     params = {"POD": True, "nTestPoints": 100, "greedyTol": 1e-4,
               "collinearityTol": 1e8, "errorEstimatorKind": "DISCREPANCY", 
               "S": 5, "polybasis": "CHEBYSHEV",
               "samplerPivot": QS([4.75, 5.25], "UNIFORM"),
               "trainSetGenerator": QS([4.75, 5.25], "CHEBYSHEV"),
-              "MMarginal": 4, "SMarginal": 5, "polybasisMarginal": "MONOMIAL",
+              "SMarginal": 5, "polybasisMarginal": "MONOMIAL",
               "samplerMarginal": QS([6.75, 7.25], "UNIFORM"),
               "matchingWeight": 1., "cutOffTolerance": 1.5}
     approx = RIGP([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), 1.181958e-02, rtol = 1)
diff --git a/tests/reduction_methods_multiD/rational_interpolant_2d.py b/tests/reduction_methods_multiD/rational_interpolant_2d.py
index 067880c..4ef22fb 100644
--- a/tests/reduction_methods_multiD/rational_interpolant_2d.py
+++ b/tests/reduction_methods_multiD/rational_interpolant_2d.py
@@ -1,77 +1,77 @@
 # 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.standard import RationalInterpolant as RI
 from rrompy.parameter.parameter_sampling import (QuadratureSampler as QS,
                                                  ManualSampler as MS)
 
 def test_monomials(capsys):
     mu = [5.05, 7.1]
     mu0 = [5., 7.]
     solver = matrixRandom()
     uh = solver.solve(mu)[0]
-    params = {"POD": False, "M": 4, "N": 4, "S": 16, "robustTol": 1e-6,
+    params = {"POD": False, "S": 16, "robustTol": 1e-6,
               "interpRcond": 1e-3, "polybasis": "MONOMIAL",
               "sampler": QS([[4.9, 6.85], [5.1, 7.15]], "UNIFORM")}
-    approx = RI(solver, mu0, params, verbosity = 0)
+    approx = RI(solver, mu0, params, verbosity = 100)
     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), 5.2667e-05, rtol = 1)
 
     out, err = capsys.readouterr()
     assert ("poorly conditioned. Reducing N " in out)
     assert len(err) == 0
     
 def test_well_cond():
     mu = [5.05, 7.1]
     mu0 = [5., 7.]
     solver = matrixRandom()
     params = {"POD": True, "M": 3, "N": 3, "S": 16,
               "interpRcond": 1e-10, "polybasis": "CHEBYSHEV",
               "sampler": QS([[4.9, 6.85], [5.1, 7.15]], "UNIFORM")}
     approx = RI(solver, mu0, params, verbosity = 0)
     approx.setupApprox()
     assert np.isclose(approx.normErr(mu)[0] / approx.normHF(mu)[0],
                       5.98695e-05, rtol = 1e-1)
 
 def test_hermite():
     mu = [5.05, 7.1]
     mu0 = [5., 7.]
     solver = matrixRandom()
     sampler0 = QS([[4.9, 6.85], [5.1, 7.15]], "UNIFORM")
     params = {"POD": True, "M": 3, "N": 3, "S": 25, "polybasis": "CHEBYSHEV",
               "sampler": MS([[4.9, 6.85], [5.1, 7.15]],
                             points = sampler0.generatePoints(9))}
     approx = RI(solver, mu0, params, verbosity = 0)
     approx.setupApprox()
     assert np.isclose(approx.normErr(mu)[0] / approx.normHF(mu)[0],
                       0.000224, rtol = 5e-1)