diff --git a/examples/3_sector_angle/sector_angle.py b/examples/3_sector_angle/sector_angle.py
index 837b0e9..fdbcedc 100644
--- a/examples/3_sector_angle/sector_angle.py
+++ b/examples/3_sector_angle/sector_angle.py
@@ -1,107 +1,107 @@
 import numpy as np
 import matplotlib.pyplot as plt
 from sector_angle_engine import SectorAngleEngine as engine
 from rrompy.reduction_methods import (NearestNeighbor as NN,
-                                      RationalInterpolantPivoted as RIP,
-                                      RationalInterpolantGreedyPivoted as RIGP)
+                             RationalInterpolantPivotedPoleMatch as RIP,
+                             RationalInterpolantGreedyPivotedPoleMatch as RIGP)
 from rrompy.parameter.parameter_sampling import (QuadratureSampler as QS,
                                                  EmptySampler as ES)
 
 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 = {'S':20, "paramsMarginal":{"MMarginal": 3}, 'SMarginal':11,
                   'POD':True, 'polybasis':"CHEBYSHEV",
                   'polybasisMarginal':"MONOMIAL_GAUSSIAN",
                   'radialDirectionalWeightsMarginal': 100.,
                   'matchingWeight':1., 'samplerPivot':QS(ks, "CHEBYSHEV", 2.),
                   'samplerMarginal':QS(ts, "UNIFORM")}
         algo = RIP
     if method == "RI_GREEDY":
         params = {'S':10, "paramsMarginal":{"MMarginal": 3}, 'SMarginal':11,
                   'POD':True, 'polybasis':"LEGENDRE",
                   'polybasisMarginal':"MONOMIAL_GAUSSIAN",
                   'radialDirectionalWeightsMarginal': 100.,
                   'matchingWeight':1., 'samplerPivot':QS(ks, "UNIFORM", 2.),
                   'greedyTol':1e-3, 'errorEstimatorKind':"LOOK_AHEAD_RES",
                   'trainSetGenerator':QS(ks, "CHEBYSHEV", 2.),
                   'samplerMarginal':QS(ts, "UNIFORM")}
         algo = RIGP
         
     approx = algo([0], solver, mu0 = [k0, t0], approxParameters = params,
                   verbosity = 10, storeAllSamples = True)
     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 ---")
     paramsNN = {'S':len(approx.mus), 'POD':True, 'sampler':ES()}
     approxNN = NN(solver, mu0 = [k0, t0], approxParameters = paramsNN,
                   verbosity = 0)
     approxNN.setSamples(approx.storedSamplesFilenames)
     approx.purgeStoredSamples()
     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/5_anisotropic_square/anisotropic_square.py b/examples/5_anisotropic_square/anisotropic_square.py
index 0ea7ff0..136baa4 100644
--- a/examples/5_anisotropic_square/anisotropic_square.py
+++ b/examples/5_anisotropic_square/anisotropic_square.py
@@ -1,80 +1,80 @@
 ### example from Smetana, Zahm, Patera. Randomized residual-based error
 ### estimators for parametrized equations.
 import numpy as np
 import matplotlib.pyplot as plt
 from itertools import product
 from anisotropic_square_engine import (AnisotropicSquareEngine as engine,
                                        AnisotropicSquareEnginePoles as plsEx)
 from rrompy.reduction_methods import (
-                               RationalInterpolantGreedyPivotedGreedy as RIGPG)
+                      RationalInterpolantGreedyPivotedGreedyPoleMatch as RIGPG)
 from rrompy.parameter.parameter_sampling import (QuadratureSampler as QS,
                                                  SparseGridSampler as SGS)
 
 zs, Ls = [10., 50.], [.2, 1.2]
 z0, L0, n = np.mean(zs), np.mean(Ls), 50
 murange = [[zs[0], Ls[0]], [zs[-1], Ls[-1]]]
 np.random.seed(4020)
 mu = [zs[0] + np.random.rand() * (zs[-1] - zs[0]),
       Ls[0] + np.random.rand() * (Ls[-1] - Ls[0])]
 solver = engine(z0, L0, n)
 
 fighandles = []
 params = {"POD": True, "nTestPoints": 100, "greedyTol": 1e-4, "S": 3,
           "polybasisMarginal": "MONOMIAL_WENDLAND",
           "polybasis": "LEGENDRE", 'samplerPivot':QS(zs, "UNIFORM"),
           'trainSetGenerator':QS(zs, "UNIFORM"),
           'errorEstimatorKind':"LOOK_AHEAD_RES",
           'errorEstimatorKindMarginal':"LOOK_AHEAD_RECOVER",
           "SMarginal": 3, "paramsMarginal": {"MMarginal": 2,
                          "radialDirectionalWeightsMarginalAdapt": [1e9, 1e12]},
           "greedyTolMarginal": 1e-2, "samplerMarginal":SGS(Ls),
           "radialDirectionalWeightsMarginal": [4.], "matchingWeight": 1.}
 
 for shared, tol in product([1., 0.], [1., 3.]):
     print("Testing cutoff tolerance {} with shared ratio {}.".format(tol,
                                                                      shared))
     solver.cutOffPolesRMinRel = - 1. - tol
     solver.cutOffPolesRMaxRel = 1. + tol
     params['sharedRatio'] = shared
 
     approx = RIGPG([0], solver, mu0 = [z0, L0], approxParameters = params,
                    verbosity = 5)
     approx.setupApprox("ALL")
     verb = approx.verbosity
     approx.verbosity = 0
     tspace = np.linspace(Ls[0], Ls[-1], 100)
     for j, t in enumerate(tspace):
         plsE = plsEx(t, 0., zs[-1])
         pls = approx.getPoles([None, t])
         pls[np.abs(np.imag(pls)) > 1e-5] = np.nan
         if j == 0:
             polesE = np.empty((len(tspace), len(plsE)))
             poles = np.empty((len(tspace), len(pls)))
             polesE[:] = np.nan
         if len(plsE) > polesE.shape[1]:
             nanR = np.empty((len(tspace), len(plsE) - polesE.shape[1]))
             nanR[:] = np.nan
             polesE = np.hstack((polesE, nanR))
         polesE[j, : len(plsE)] = np.real(plsE)
         poles[j] = np.real(pls)
     approx.verbosity = verb
     fighandles += [plt.figure(figsize = (17, 5))]
     ax1 = fighandles[-1].add_subplot(1, 2, 1)
     ax2 = fighandles[-1].add_subplot(1, 2, 2)
     ax1.plot(poles, tspace)
     ax1.set_ylim(Ls)
     ax1.set_xlabel('mu_1')
     ax1.set_ylabel('mu_2')
     ax1.grid()
     ax2.plot(polesE, tspace, 'k-.', linewidth = 1)
     ax2.plot(poles, tspace)
     for mm in approx.musMarginal:
         ax2.plot(zs, [mm[0, 0]] * 2, 'k--', linewidth = 1)
     ax2.set_xlim(zs)
     ax2.set_ylim(Ls)
     ax2.set_xlabel('mu_1')
     ax2.set_ylabel('mu_2')
     ax2.grid()
     plt.show()
 
     print("\n")
diff --git a/examples/8_damped_mass_chain/damped_mass_chain.py b/examples/8_damped_mass_chain/damped_mass_chain.py
index e73f135..cae7247 100644
--- a/examples/8_damped_mass_chain/damped_mass_chain.py
+++ b/examples/8_damped_mass_chain/damped_mass_chain.py
@@ -1,184 +1,184 @@
 ### example from Lohmann, Eid. Efficient Order Reduction of Parametric and
 ### Nonlinear Models by Superposition of Locally Reduced Models.
 from copy import deepcopy as copy
 import numpy as np
 import matplotlib.pyplot as plt
 from rrompy.reduction_methods import (NearestNeighbor as NN,
-                                      RationalInterpolant as RI,
-                                      RationalInterpolantGreedy as RIG,
-                                      RationalInterpolantPivoted as RIP,
-                                      RationalInterpolantGreedyPivoted as RIGP)
+                             RationalInterpolant as RI,
+                             RationalInterpolantGreedy as RIG,
+                             RationalInterpolantPivotedPoleMatch as RIP,
+                             RationalInterpolantGreedyPivotedPoleMatch as RIGP)
 from rrompy.parameter.parameter_sampling import (QuadratureSampler as QS,
                                                  EmptySampler as ES)
 from damped_mass_chain_engine import (bode as bode0, bodeLog, MassChainEngine,
      MassChainEngineLog, AugmentedMassChainEngine, AugmentedMassChainEngineLog)
 from rrompy.utilities.base.decorators import addWhiteNoise
 
 ##########################
 fullModelOrder = 1 #+ 1
 SMarginal = 1 #* 2
 state = 0 #+ 1
 noise_level = 0 #+ 1e-5
 LS = 0 #+ 6
 ##########################
 
 modelSign = "Surrogate modeling for frequency response of "
 if fullModelOrder == 1: modelSign += "augmented "
 modelSign += "damper-mass-spring model"
 if SMarginal > 1: modelSign += " with 1 design parameter"
 modelSign += ".\nOutput is "
 if state:
     modelSign += "vector of mass displacements. "
 else:
     modelSign += "displacement of last mass. "
 if LS:
     modelSign += "Least-squares: S - N - 1 = {}. ".format(LS)
 else:
     modelSign += "Interpolatory: S = N + 1. "
 modelSign += "Noise level: {}.\n".format(noise_level)
 print(modelSign)
 
 M = [np.array([1., 5., 25., 125.])]
 N = len(M[0])
 D = [np.zeros((N, N))]
 D[0][0, 1], D[0][1, 2], D[0][2, 3], D[0][3, 3] = .1, .4, 1.6, 0.
 D[0] = D[0] + D[0].T
 K = [np.zeros((N, N))]
 K[0][0, 1], K[0][1, 2], K[0][2, 3], K[0][0, 3], K[0][3, 3] = 9., 3., 1., 1., 2.
 K[0] = K[0] + K[0].T
 B = np.append(27., np.zeros(N - 1)).reshape(-1, 1)
 if SMarginal > 1:
     M += [np.zeros(N)]
     D += [np.zeros((N, N))]
     D[1][3, 3] = 1.
     D[1] = D[1] + D[1].T
     K += [np.zeros((N, N))]
     K[1][0, 3], K[1][3, 3] = 2., 2.
     K[1] = K[1] + K[1].T
 if state:
     C = np.eye(4)
 else:
     C = np.append(np.zeros(N - 1), 1.).reshape(1, -1)
 
 for logspace in range(2):
     print("Approximation in l{}space".format("og" * logspace
                                            + "in" * (not logspace)))
     if logspace:
         bode = bodeLog
         if fullModelOrder == 1:
             engine = AugmentedMassChainEngineLog
         else:
             engine = MassChainEngineLog
     else:
         bode = bode0
         if fullModelOrder == 1:
             engine = AugmentedMassChainEngine
         else:
             engine = MassChainEngine
     solver = addWhiteNoise(noise_level)(engine)(M, D, K, B, C)
 
     ss, mu = [1e-2, 1e1], []
     s0 = 10. ** np.mean(np.log10(ss))
     freq = np.logspace(np.log10(ss[0]), np.log10(ss[1]), 100)
     if logspace:
         ss, freq = [np.log10(ss[0]), np.log10(ss[1])], np.log10(freq)
         s0, parameterMap = np.log10(s0), 1.
     else:
         parameterMap = {"F": [("log10", "x")], "B": [(10., "**", "x")]}
     krange = [[ss[0]], [ss[-1]]]
     k0, srange = [s0], copy(krange)
     if SMarginal > 1:
         ms = [0., 1.]
         m0, mrange = np.mean(ms), [[ms[0]], [ms[-1]]]
         krange[0] += mrange[0]
         krange[1] += mrange[1]
         k0 += [m0]
         mu = [.5 * (ms[1] - ms[0]) / (SMarginal - 1)]
         if not logspace:
             parameterMap["F"] += [("x")]
             parameterMap["B"] += [("x")]
     
     for method in ["RI", "RI_GREEDY"]:
         print("Testing {} method".format(method))
         
         if method == "RI":
             params = {'S':15, 'POD':True, 'polybasis':"CHEBYSHEV"}
             if LS: params["N"] = params["S"] - 1 - LS
             if SMarginal > 1:
                 algo = RIP
             else:
                 params['sampler'] = QS(srange, "CHEBYSHEV", parameterMap)
                 algo = RI
         if method == "RI_GREEDY":
             params = {'S':5, 'POD':True, 'polybasis':"LEGENDRE",
                       'greedyTol':1e-2, 'errorEstimatorKind':"DISCREPANCY",
                       'trainSetGenerator':QS(srange, "CHEBYSHEV",
                                              parameterMap)}
             if SMarginal > 1:
                 algo = RIGP
             else:
                 params['sampler'] = QS(srange, "UNIFORM", parameterMap)
                 algo = RIG
         if SMarginal > 1:
             params["paramsMarginal"] = {"MMarginal": SMarginal - 1}
             params['SMarginal'] = SMarginal
             params['polybasisMarginal'] = "MONOMIAL"
             params['radialDirectionalWeightsMarginal'] = [2. / (ms[1] - ms[0])]
             params['matchingWeight'] = 1.
             params['samplerPivot'] = QS(srange, "UNIFORM", parameterMap)
             params['samplerMarginal'] = QS(mrange, "UNIFORM")
             approx = algo([0], solver, mu0 = k0, approxParameters = params,
                           verbosity = 5, storeAllSamples = True)
         else:
             approx = algo(solver, mu0 = k0, approxParameters = params,
                           verbosity = 5)
         if "GREEDY" in method:
             approx.setupApprox("LAST")
         else:
             approx.setupApprox()
         
         approxNN = NN(solver, mu0 = k0, verbosity = 5,
                       approxParameters = {'S':len(approx.mus),
                                           'POD':params['POD'], 'sampler':ES()})
         if SMarginal > 1:
             approxNN.setSamples(approx.storedSamplesFilenames)
             approx.purgeStoredSamples()
             for m in approx.musMarginal:
                 bode(freq, m[0],
                      [approx.getHF, approx.getApprox, approxNN.getApprox])
         else:
             approxNN.setSamples(approx.samplingEngine)
         bode(freq, mu, [approx.getHF, approx.getApprox, approxNN.getApprox])
         if SMarginal > 1:
             bode(freq, [1.5 * ms[1]],
                  [approx.getHF, approx.getApprox, approxNN.getApprox])
             bode(freq, [2. * ms[1]],
                  [approx.getHF, approx.getApprox, approxNN.getApprox])
             verb = approx.verbosity
             approx.verbosity = 0
             mspace = np.linspace(ms[0], ms[-1], 10)
             for j, t in enumerate(mspace):
                 pls = approx.getPoles([None, t])
                 if j == 0:
                     poles = np.empty((len(mspace), len(pls)),
                                      dtype = np.complex)
                 poles[j] = pls
             for j, t in enumerate(approx.musMarginal):
                 pls = approx.getPoles([None, t[0][0]])
                 if j == 0:
                     polesE = np.empty((SMarginal, len(pls)),
                                       dtype = np.complex)
                 polesE[j] = pls
             approx.verbosity = verb
             fig = plt.figure(figsize = (10, 6))
             ax = fig.add_subplot(1, 1, 1)
             ax.plot(np.real(poles), np.imag(poles), '--')
             ax.plot(np.real(polesE), np.imag(polesE), 'ko', markersize = 4)
             ax.set_xlabel('Real')
             ax.set_ylabel('Imag')
             ax.grid()
             plt.show()
         else:
             poles = approx.getPoles()
             print("Poles:\n{}".format(poles))
         print("\n")
diff --git a/examples/9_active_remeshing/active_remeshing.py b/examples/9_active_remeshing/active_remeshing.py
index 204b8c7..3071410 100755
--- a/examples/9_active_remeshing/active_remeshing.py
+++ b/examples/9_active_remeshing/active_remeshing.py
@@ -1,117 +1,119 @@
 import numpy as np
 from pickle import load
 from matplotlib import pyplot as plt
 from active_remeshing_engine import ActiveRemeshingEngine
 from rrompy.reduction_methods import (
-                      RationalInterpolantGreedyPivotedGreedyNoMatch as RIGPGNM,
-                      RationalInterpolantGreedyPivotedGreedy as RIGPG)
+                      RationalInterpolantGreedyPivotedNoMatch as RIGPNM,
+                      RationalInterpolantGreedyPivotedGreedyPoleMatch as RIGPG)
 from rrompy.parameter.parameter_sampling import (QuadratureSampler as QS,
                                                  SparseGridSampler as SGS)
 
 zs, ts = [0., 100.], [0., .5]
-z0, t0, n = np.mean(zs), np.mean(ts), 150
+z0, t0, n = np.mean(zs), np.mean(ts), 15
 solver = ActiveRemeshingEngine(z0, t0, n)
 
 mus = [[z0, ts[0]], [z0, ts[1]]]
 for mu in mus:
     u = solver.solve(mu, return_state = True)[0]
     Y = solver.applyC(u, mu)[0][0]
     _ = solver.plot(u, what = "REAL", name = "u(z={}, t={})".format(*mu),
                     is_state = True, figsize = (12, 4))
     print("Y(z={}, t={}) = {} (solution norm squared)".format(*mu, Y))
 
 fighandles = []
 with open("./active_remeshing_hf_samples.pkl", "rb") as f:
     zspace, tspace, Yex = load(f)
 # zspace = np.linspace(zs[0], zs[-1], 198)
 # tspace = np.linspace(ts[0], ts[-1], 50)
 # Yex = np.log10([[solver.solve([z, t]) for t in tspace] for z in zspace])
 # (from a ~2h45m simulation on one node of the EPFL Helvetios cluster)
 
-for match in [0, 1]:
-    params = {"POD": True, "S": 5, "SMarginal": 3, "greedyTol": 1e-4,
-              "nTestPoints": 500, "polybasis": "LEGENDRE",
-              "maxIterMarginal": 25, "trainSetGenerator": QS(zs, "UNIFORM"),
+for match in [1]:
+#for match in [0, 1]:
+    params = {"POD": True, "S": 5, "greedyTol": 1e-4, "nTestPoints": 500,
+              "polybasis": "LEGENDRE", "trainSetGenerator": QS(zs, "UNIFORM"),
               "samplerPivot":QS(zs, "CHEBYSHEV"), "samplerMarginal":SGS(ts),
-              "errorEstimatorKind": "LOOK_AHEAD_OUTPUT",
-              "errorEstimatorKindMarginal": "LOOK_AHEAD_RECOVER"}
+              "errorEstimatorKind": "LOOK_AHEAD_OUTPUT"}
     if match:
         print("\nTesting output-based matching with weight 1.")
+        params["SMarginal"] = 3
+        params["maxIterMarginal"] = 25
         params["greedyTolMarginal"] = 1e-2
         params["matchingWeight"] = 1.
         params["sharedRatio"] = .75
         params["polybasisMarginal"] = "PIECEWISE_LINEAR_UNIFORM"
+        params["errorEstimatorKindMarginal"] = "LOOK_AHEAD_RECOVER"
         algo = RIGPG
     else:
         print("\nTesting matching-free approach.")
-        params["greedyTolMarginal"] = 2e-2
-        algo = RIGPGNM
+        params["SMarginal"] = 5
+        algo = RIGPNM
 
     approx = algo([0], solver, mu0 = [z0, t0], approxParameters = params,
                   verbosity = 15)
-    approx.setupApprox("ALL")
+    approx.setupApprox("ALL" * match)
 
     verb, verbTM = approx.verbosity, approx.trainedModel.verbosity
     approx.verbosity, approx.trainedModel.verbosity = 0, 0
     for j, t in enumerate(tspace):
         out = approx.getApprox(np.pad(zspace.reshape(-1, 1), [(0, 0), (0, 1)],
                                       "constant", constant_values = t))
         pls = approx.getPoles([None, t])
         pls[np.abs(np.imag(pls)) > 1e-5] = np.nan
         if j == 0:
             Ys = np.empty((len(zspace), len(tspace)))
             poles = np.empty((len(tspace), len(pls)))
         Ys[:, j] = np.log10(out.data)
         if len(pls) > poles.shape[1]:
             poles = np.pad(poles, [(0, 0), (0, len(pls) - poles.shape[1])],
                            "constant", constant_values = np.nan)
         poles[j, : len(pls)] = np.real(pls)
     approx.verbosity, approx.trainedModel.verbosity = verb, verbTM
     Ymin, Ymax = min(np.min(Ys), np.min(Yex)), max(np.max(Ys), np.max(Yex))
     
     fighandles += [plt.figure(figsize = (15, 5))]
     ax1 = fighandles[-1].add_subplot(1, 2, 1)
     ax2 = fighandles[-1].add_subplot(1, 2, 2)
     if match:
         ax1.plot(poles, tspace)
     else:
         ax1.plot(poles, tspace, 'k.')
     ax1.set_ylim(ts)
     ax1.set_xlabel('z')
     ax1.set_ylabel('t')
     ax1.grid()
     if match:
         ax2.plot(poles, tspace)
     else:
         ax2.plot(poles, tspace, 'k.')
     for mm in approx.musMarginal:
         ax2.plot(zs, [mm[0, 0]] * 2, 'k--', linewidth = 1)
     ax2.set_xlim(zs)
     ax2.set_ylim(ts)
     ax2.set_xlabel('z')
     ax2.set_ylabel('t')
     ax2.grid()
     plt.show()
     print("Approximate poles")
 
     fighandles += [plt.figure(figsize = (15, 5))]
     ax1 = fighandles[-1].add_subplot(1, 2, 1)
     ax2 = fighandles[-1].add_subplot(1, 2, 2)
     p = ax1.contourf(np.repeat(zspace.reshape(-1, 1), len(tspace), axis = 1),
                      np.repeat(tspace.reshape(1, -1), len(zspace), axis = 0),
                      Ys, vmin = Ymin, vmax = Ymax, cmap = "gray_r",
                      levels = np.linspace(Ymin, Ymax, 50))
     plt.colorbar(p, ax = ax1)
     ax1.set_xlabel('z')
     ax1.set_ylabel('t')
     ax1.grid()
     p = ax2.contourf(np.repeat(zspace.reshape(-1, 1), len(tspace), axis = 1),
                     np.repeat(tspace.reshape(1, -1), len(zspace), axis = 0),
                     Yex, vmin = Ymin, vmax = Ymax, cmap = "gray_r",
                     levels = np.linspace(Ymin, Ymax, 50))
     ax2.set_xlabel('z')
     ax2.set_ylabel('t')
     ax2.grid()
     plt.colorbar(p, ax = ax2)
     plt.show()
     print("Approximate and exact output\n")
diff --git a/rational_interpolation_method.pdf b/rational_interpolation_method.pdf
index be1163d..a6088f0 100644
Binary files a/rational_interpolation_method.pdf and b/rational_interpolation_method.pdf differ
diff --git a/rational_interpolation_method.tex b/rational_interpolation_method.tex
index 441fa90..1182a8d 100644
--- a/rational_interpolation_method.tex
+++ b/rational_interpolation_method.tex
@@ -1,219 +1,225 @@
 \documentclass[10pt,a4paper]{article}
 \usepackage[left=1in,right=1in,top=1in,bottom=1in]{geometry}
 \usepackage{amsmath}
 \usepackage{amsfonts}
 \usepackage{amssymb}
 \usepackage{hyperref}
 \usepackage{xcolor}
 
 \setlength{\parindent}{0pt}
 \newcommand{\code}[1]{{\color{blue}\texttt{#1}}}
 \newcommand{\footpath}[1]{\footnote{\path{#1}}}
 \newcommand{\N}{\mathbb{N}}
 \newcommand{\R}{\mathbb{R}}
 \newcommand{\C}{\mathbb{C}}
 \newcommand{\I}{\mathcal{I}}
 \DeclareMathOperator*{\argmin}{arg\,min}
 \newcommand{\inner}[2]{\left\langle#1,#2\right\rangle_V}
 \newcommand{\norm}[1]{\left\|#1\right\|_V}
 
-\title{\bf The RROMPy rational interpolation method}
+\title{\bf The \code{RROMPy} rational interpolation method}
 \author{D. Pradovera, CSQI, EPF Lausanne -- \texttt{davide.pradovera@epfl.ch}}
 \date{}
+
+\hypersetup{ pdftitle={The RROMPy rational interpolation method}, pdfauthor={Davide Pradovera} }
+
 \begin{document}
 \maketitle
 
 \section*{Introduction}
 This document provides an explanation for the numerical method provided by the class \code{Rational Interpolant}\footpath{./rrompy/reduction_methods/standard/rational_interpolant.py} and daughters, e.g. \code{Rational Interpolant Greedy}\footpath{./rrompy/reduction_methods/standard/greedy/rational_interpolant_greedy.py}, as well as most of the pivoted approximants\footpath{./rrompy/reduction_methods/pivoted/{,greedy/}rational_interpolant_*.py}.
 
 We restrict the discussion to the single-parameter case, and most of the focus will be dedicated to the impact of the \code{functionalSolve} parameter, whose allowed values are
 \begin{itemize}
 \item \code{NORM} (default): see \ref{sec:norm}; allows for derivative information, i.e. repeated sample points.
 \item \code{DOMINANT}: see \ref{sec:dominant}; allows for derivative information, i.e. repeated sample points.
 \item \code{BARYCENTRIC\_NORM}: see \ref{sec:barycentricnorm}; does not allow for a Least Squares (LS) approach.
 \item \code{BARYCENTRIC\_AVERAGE}: see \ref{sec:barycentricaverage}; does not allow for a Least Squares (LS) approach.
 \item \code{NODAL}: see \ref{sec:nodal}; iterative method.
 \end{itemize}
 
 The main reference throughout the present document is \cite{Pradovera}.
 
 \section{Aim of approximation}
-We seek an approximation of $u:\C\to V$, with $(V,\inner{\cdot}{\cdot})$ a complex\footnote{The inner product is linear (resp. conjugate linear) in the first (resp. second) argument: $\inner{\alpha v}{\beta w}=\alpha\overline{\beta}\inner{v}{w}$.} Hilbert space (with endowed norm $\norm{\cdot}$), of the form $\widehat{p}/\widehat{q}$, where $\widehat{p}:\C\to V$ and $\widehat{q}:\C\to\C$. For a given denominator $\widehat{q}$, the numerator $\widehat{p}$ is found by interpolation (possibly, LS or based on radial basis functions) of $\widehat{q}u$. Hence, here we focus on the computation of the denominator $\widehat{q}$.
+We seek an approximation of $u:\C\to V$, with $(V,\inner{\cdot}{\cdot})$ a complex\footnote{The inner product is linear (resp. conjugate linear) in the first (resp. second) argument: $\inner{\alpha v}{\beta w}=\alpha\overline{\beta}\inner{v}{w}$.} Hilbert space (with induced norm $\norm{\cdot}$), of the form $\widehat{p}/\widehat{q}$, where $\widehat{p}:\C\to V$ and $\widehat{q}:\C\to\C$. For a given denominator $\widehat{q}$, the numerator $\widehat{p}$ is found by interpolation (possibly, LS or based on radial basis functions) of $\widehat{q}u$. Hence, here we focus on the computation of the denominator $\widehat{q}$.
 
 Other than the choice of target function $u$, the parameters which affect the computation of $\widehat{q}$ are:
 \begin{itemize}
 \item $\code{mus}\subset\C$ ($\{\mu_j\}_{j=1}^S$ below); for all \code{functionalSolve} values but \code{NORM} and \code{DOMINANT}, the $S$ points must be distinct.
 \item $\code{N}\in\N$ ($N$ below); for \code{BARYCENTRIC}, $N$ must equal $S-1$.
 \item $\code{polybasis0}\in\{\code{"CHEBYSHEV"}, \code{"LEGENDRE"}, \code{"MONOMIAL"}\}$; only for \code{NORM} and \code{DOMINANT}.
 \end{itemize}
 For simplicity, we will consider only the case of $S$ distinct sample points. One can deal with the case of confluent sample points by extending the standard (Lagrange) interpolation steps to Hermite-Lagrange ones.
 
 The main motivation behind the method involves the modified approximation problem
 \[u\approx\I^N\left(\Big(\big(\mu_j,\widehat{q}(\mu_j)u(\mu_j)\big)\Big)_{j=1}^S\right)\Big/\widehat{q},\]
 where $\widehat{q}:\C\to\C$ is a polynomial of degree $\leq N<S$, and $\I^N:(\C\times V)^S\to\mathbb{P}^N(\C;V)$ is a (LS) polynomial interpolation operator, which maps $S$ samples of a function (which lie in $V$) to a polynomial of degree $\leq N$ with coefficients in $V$.
 
 More precisely, let
 \[\mathbb{P}^N(\C;W)=\left\{\mu\mapsto\sum_{i=0}^N\alpha_i\mu^i\ :\ \alpha_0,\ldots,\alpha_N\in W\right\},\]
 with $W$ either $\C$ or $V$. We set
 \[\I^N\left(\big((\mu_j,\psi_j)\big)_{j=1}^S\right)\bigg|_\mu=\argmin_{p\in\mathbb{P}^N(\C;V)}\sum_{j=1}^S\norm{p(\mu_j)-\psi_j}^2.\]
-In RROMPy, we compute (LS-)interpolants by employing normal equations: given a basis $\{\phi_i\}_{i=0}^N$ of $\mathbb{P}^N(\C;\C)$, we expand
+In \code{RROMPy}, we compute (LS-)interpolants by employing normal equations: given a basis $\{\phi_i\}_{i=0}^N$ of $\mathbb{P}^N(\C;\C)$, we expand
 \[\I^N\left(\big((\mu_j,\psi_j)\big)_{j=1}^S\right)=\sum_{i=0}^Nc_i\phi_i\]
 and observe that, for optimality, the coefficients $\{c_i\}_{i=0}^N\subset V$ must satisfy
 \[\sum_{j=1}^S\sum_{l=0}^N\overline{\phi_i(\mu_j)}\phi_l(\mu_j)c_l=\sum_{j=1}^S\overline{\phi_i(\mu_j)}\psi_j\quad\forall i=0,\ldots,N,\]
 i.e., in matrix form\footnote{The superscript ${}^H$ denotes adjunction (conjugate transposition), i.e. $A^H=\overline{A}^\top$.},
 \[\underbrace{[c_i]_{i=0}^N}_{\in V^{N+1}}=\bigg(\underbrace{\Big[\overline{\phi_i(\mu_j)}\Big]_{i=0,j=1}^{N,S}}_{=:\Phi^H\in\C^{(N+1)\times S}}\underbrace{\Big[\phi_i(\mu_j)\Big]_{j=1,i=0}^{S,N}}_{=:\Phi\in\C^{S\times(N+1)}}\bigg)^{-1}\underbrace{\Big[\overline{\phi_i(\mu_j)}\Big]_{i=0,j=1}^{N,S}}_{=:\Phi^H\in\C^{(N+1)\times S}}\underbrace{[\psi_j]_{j=1}^S}_{\in V^{S}}.\]
 In practice the polynomial basis $\{\phi_i\}_{i=0}^N$ is determined by the value of $\code{polybasis0}$:
 \begin{itemize}
 \item If $\code{polybasis0}=\code{"CHEBYSHEV"}$, then $\phi_k(\mu)=\mu^k$ for $k\in\{0,1\}$ and $\phi_k(\mu)=2\mu\phi_{k-1}(\mu)-\phi_{k-2}(\mu)$ for $k\geq2$.
 \item If $\code{polybasis0}=\code{"LEGENDRE"}$, then $\phi_k(\mu)=\mu^k$ for $k\in\{0,1\}$ and $\phi_k(\mu)=(2-1/k)\mu\phi_{k-1}(\mu)-(1-1/k)\phi_{k-2}(\mu)$ for $k\geq2$.
 \item If $\code{polybasis0}=\code{"MONOMIAL"}$, then $\phi_k(\mu)=\mu^k$ for $k\geq0$.
 \end{itemize}
 
 The actual denominator $\widehat{q}$ is found as
 \begin{equation}\label{eq:glob1}
 \widehat{q}=\argmin_{\substack{q\in\mathbb{P}^N(\C;\C)\\(\star)}}\norm{\frac{\textrm{d}^N}{\textrm{d}\mu^N}\I^N\left(\Big(\big(\mu_j,q(\mu_j)u(\mu_j)\big)\Big)_{j=1}^S\right)}
 \end{equation}
 where $(\star)$ is a normalization condition (which changes depending on \code{functionalSolve}) to exclude the trivial minimizer $\widehat{q}=0$. 
 
 Broadly speaking, the five methods described differ in terms of the constraint $(\star)$, as well as of the degrees of freedom which are chosen to represent the denominator $q$.
 
 \section{Polynomial coefficients as degrees of freedom}
 
 If the polynomial basis $\{\phi_i\}_{i=0}^N$ is hierarchical (as the three ones above), then the $N$-th derivative of $\I^N$ coincides with the coefficient $c_N$, and we have
 \begin{equation}\label{eq:poly1}
 \widehat{q}=\argmin_{\substack{q\in\mathbb{P}^N(\C;\C)\\(\star)}}\norm{[\underbrace{0,\ldots,0}_{N},1]\big(\Phi^H\Phi\big)^{-1}\Phi^H[q(\mu_j)u(\mu_j)]_{j=1}^S}.
 \end{equation}
 
 Using the Kronecker delta ($\delta_{ij}=1$ if $i=j$ and $\delta_{ij}=0$ if $i\neq j$), the last term $[q(\mu_j)u(\mu_j)]_{j=1}^S\in V^S$ can be factored into
 \begin{equation}\label{eq:poly2}
 \Big[u(\mu_j)\delta_{jj'}\Big]_{j=1,j'=1}^{S,S}[q(\mu_j)]_{j=1}^S=\Big[u(\mu_j)\delta_{jj'}\Big]_{j=1,j'=1}^{S,S}\Phi[q_i]_{i=0}^N,
 \end{equation}
 where we have expanded the polynomial $q$ using the basis\footnote{In theory, nothing prevents us from using different bases for $\I^N$ and $q$, cf. \ref{sec:observations}.} $\{\phi_i\}_{i=0}^N$: $q(\mu)=\sum_{i=0}^Nq_i\phi_i(\mu)$, with coefficients $\{q_i\}_{i=0}^N\subset\C$.
 
 Combining \eqref{eq:poly1} and \eqref{eq:poly2}, it is useful to consider the $(N+1)\times(N+1)$ Hermitian matrix with entries ($0\leq i,i'\leq N$)
 \begin{equation}\label{eq:poly3}
 G_{ii'}=\inner{\sum_{j'=1}^S\left(\big(\Phi^H\Phi\big)^{-1}\Phi^H\right)_{Nj'}(\Phi)_{j'i'}u(\mu_{j'})}{\sum_{j=1}^S\left(\big(\Phi^H\Phi\big)^{-1}\Phi^H\right)_{Nj}(\Phi)_{ji}u(\mu_j)}.
 \end{equation}
 If $(\star)$ is quadratic (resp. linear) in $[q_i]_{i=0}^N$, then we can cast the computation of the denominator as a quadratically (resp. linearly) constrained quadratic program involving $G$.
 
 \subsection{Quadratic constraint}\label{sec:norm}
 We constrain $[\widehat{q}_i]_{i=0}^N$ to have unit (Euclidean) norm. The resulting optimization problem can be cast as a minimal (normalized) eigenvector problem for $G$ in \eqref{eq:poly3}. More explicitly,
 \[[\widehat{q}_i]_{i=0}^N=\argmin_{\substack{\mathbf{q}\in\C^{N+1}\\\|\mathbf{q}\|_2=1}}\mathbf{q}^HG\mathbf{q}.\]
 
 \subsection{Linear constraint}\label{sec:dominant}
 We constrain $\widehat{q}_N=1$, thus forcing $q$ to be monic, with degree exactly $N$. Given $G$ in \eqref{eq:poly3}, the resulting optimization problem can be solved rather easily as:
 \[[\widehat{q}_i]_{i=0}^N=\frac{G^{-1}\mathbf{e}_{N+1}}{\mathbf{e}_{N+1}^\top G^{-1}\mathbf{e}_{N+1}},\quad\text{with }\mathbf{e}_{N+1}=[0,\ldots,0,1]^\top\in\C^{N+1}.\]
 
 
-\section{Barycentric coefficients as degrees of freedom}
+\section{Barycentric coefficients as degrees of freedom}\label{sec:bary}
 Here we assume that the sample points are distinct, and such that $N=S-1$. We can choose for convenience a non-hierarchical basis, dependent on the sample points, for $q$ and $\I^N$, taking inspiration from barycentric interpolation:
 \begin{equation}\label{eq:bary1}
 \phi_i(\mu)=\prod_{\substack{j=1\\j\neq i+1}}^S(\mu-\mu_j).
 \end{equation}
 Since all elements of the basis are monic and of degree exactly $N$, the minimization problem can be cast as
 \begin{equation}\label{eq:bary2}
 \widehat{q}=\argmin_{\substack{q\in\mathbb{P}^N(\C;\C)\\(\star)}}\norm{[\underbrace{1,\ldots,1}_{N+1}]\big(\Phi^H\Phi\big)^{-1}\Phi^H\Big[u(\mu_j)\delta_{jj'}\Big]_{j=1,j'=1}^{S,S}\Phi[q_i]_{i=0}^N}.
 \end{equation}
 At the same time, it is easy to see from \eqref{eq:bary1} that the Vandermonde-like matrix $\Phi$ is diagonal, so that
 \[\big(\Phi^H\Phi\big)^{-1}\Phi^H\Big[u(\mu_j)\delta_{jj'}\Big]_{j=1,j'=1}^{S,S}\Phi=\big(\Phi^H\Phi\big)^{-1}\Phi^H\Phi\Big[u(\mu_j)\delta_{jj'}\Big]_{j=1,j'=1}^{S,S}=\Big[u(\mu_j)\delta_{jj'}\Big]_{j=1,j'=1}^{S,S},\]
 and
 \begin{equation}\label{eq:bary3}
 \widehat{q}=\argmin_{\substack{\mathbf{q}\in\C^{N+1}\\(\star)}}\norm{\sum_{i=0}^Nu(\mu_{i+1})q_i}.
 \end{equation}
 
 Considering \eqref{eq:bary3}, it is useful to define the $(N+1)\times(N+1)$ Hermitian (``snapshot Gramian'') matrix with entries ($0\leq i,i'\leq N$)
 \begin{equation}\label{eq:bary4}
 G_{ii'}=\inner{u(\mu_{i'+1})}{u(\mu_{i+1})}.
 \end{equation}
 So, once again, if $(\star)$ is quadratic (resp. linear) in $[q_i]_{i=0}^N$, then we can cast the computation of the denominator as a quadratically (resp. linearly) constrained quadratic program involving $G$.
 
-Before specifying the kind of normalization enforced, it is important to make a remark on numerical stability. The basis in \eqref{eq:bary1} is actually just a ($i$-dependent) factor away from being the Lagrangian one (for which $\phi_i(\mu_j)$ would equal $\delta_{(i+1)j}$ instead of $\delta_{(i+1)j}\prod_{k\neq i+1}(\mu_j-\mu_k)$ as is does in our case). As such, it is generally a bad idea to numerically evaluate $q$ starting from its expansion coefficients with respect to $\{\phi_i\}_{i=0}^N$. We get around this by exploiting the following trick, whose foundation is in \cite[Section 2.3.3]{Klein}: the roots of $\widehat{q}=\sum_{i=0}^N\widehat{q}_i\phi_i$ are the $N$ finite eigenvalues $\lambda$ of the generalized $(N+2)\times(N+2)$ eigenproblem
+Before specifying the kind of normalization enforced, it is important to make a remark on numerical stability. The basis in \eqref{eq:bary1} is actually just a ($i$-dependent) factor away from being the Lagrangian one (for which $\phi_i(\mu_j)$ would equal $\delta_{(i+1)j}$ instead of $\delta_{(i+1)j}\prod_{k\neq i+1}(\mu_j-\mu_k)$ as it does in our case). As such, it is generally a bad idea to numerically evaluate $q$ starting from its expansion coefficients with respect to $\{\phi_i\}_{i=0}^N$. We get around this by exploiting the following trick, whose foundation is in \cite[Section 2.3.3]{Klein}: the roots of $\widehat{q}=\sum_{i=0}^N\widehat{q}_i\phi_i$ are the $N$ finite eigenvalues $\lambda$ of the generalized $(N+2)\times(N+2)$ eigenproblem
 \begin{equation}
-\textrm{det}\left(\begin{bmatrix}
+\textrm{Det}\left(\begin{bmatrix}
 0 & \widehat{q}_0 & \cdots & \widehat{q}_N\\
 1 & \mu_1 &  & \\
 \vdots & & \ddots & \\
 1 & & & \mu_S
 \end{bmatrix}-\begin{bmatrix}
 0 & & & \\
  & 1 &  & \\
  & & \ddots & \\
  & & & 1
 \end{bmatrix}\lambda\right)=0.
 \end{equation}
-This computation is numerically more stable than most manipulations of a polynomial in the basis \eqref{eq:bary1}.
+This computation is numerically more stable than most manipulations of a polynomial in the basis \eqref{eq:bary1}. Unfortunately, the ensuing computation of $\widehat{p}$ invariably suffers from numerical instabilities as a consequence of this.
 
-Once the roots of $\widehat{q}$ have been computed, one can either convert it to nodal form \eqref{eq:pole1} or forgo using $\widehat{q}$ completely, in favor of a Heaviside like approximation involving the newly computed roots $\{\widehat{\lambda}_i\}_{i=1}^N$ as poles:
+Once the roots of $\widehat{q}$ have been computed, one can either convert it to nodal form \eqref{eq:pole1} or forgo using $\widehat{q}$ completely, in favor of a Heaviside-like approximation involving the newly computed roots $\{\widehat{\lambda}_i\}_{i=1}^N$ as poles:
 \[
 \frac{\widehat{p}(\mu)}{\widehat{q}(\mu)}\quad\rightsquigarrow\quad\sum_{i=1}^N\frac{\widehat{r}_i}{\mu-\widehat{\lambda}_i}+\widetilde{p}(\mu),
 \]
-with $\widetilde{p}$, e.g., a polynomial (of degree at most $S-N-1$) or a combination of radial basis functions.
+with $\widetilde{p}$, e.g., a polynomial (of degree at most $S-N-1$) or a combination of radial basis functions. See the final paragraph in \ref{sec:observations} for a slightly more detailed motivation of why the Heaviside form of the approximant might be more useful than the standard rational one $\widehat{p}/\widehat{q}$ in practice.
 
 \subsection{Quadratic constraint}\label{sec:barycentricnorm}
 We constrain $[\widehat{q}_i]_{i=0}^N$ to have unit (Euclidean) norm. The resulting optimization problem can be cast as a minimal (normalized) eigenvector problem for $G$ in \eqref{eq:bary4}. More explicitly,
 \[[\widehat{q}_i]_{i=0}^N=\argmin_{\substack{\mathbf{q}\in\C^{N+1}\\\|\mathbf{q}\|_2=1}}\mathbf{q}^HG\mathbf{q}.\]
 
 \subsection{Linear constraint}\label{sec:barycentricaverage}
 We constrain $\sum_{i=0}^N\widehat{q}_i=1$, so that the polynomial $\widehat{q}$ is monic. Given $G$ in \eqref{eq:bary4}, the resulting optimization problem can be solved rather easily as:
 \[[\widehat{q}_i]_{i=0}^N=\frac{G^{-1}\mathbf{1}_{N+1}}{\mathbf{1}_{N+1}^\top G^{-1}\mathbf{1}_{N+1}},\quad\text{with }\mathbf{1}_{N+1}=[1,\ldots,1]^\top\in\C^{N+1}.\]
 
 
 \section{Poles as degrees of freedom}\label{sec:nodal}
 Here we assume that the sample points are distinct. One can avoid expressing $q$ explicitly and work directly with its roots by employing a nodal form
 \begin{equation}\label{eq:pole1}
 q(\mu)=\prod_{i=1}^N(\mu-\lambda_i),
 \end{equation}
-with $\{\lambda_i\}_{i=1}^N\subset\C$. (It is important to note that we are constraining $q$ to have degree exactly $N$.)
+with $\{\lambda_i\}_{i=1}^N\subset\C$. (It is important to note that, as in \ref{sec:bary}, we are constraining $q$ to have degree exactly $N$.)
 
 The optimization problem that allows to find the desired poles can now be cast as
 \[\{\widehat{\lambda}_i\}_{i=1}^N=\argmin_{\{\lambda_i\}_{i=1}^N\subset\C}\norm{\frac{\textrm{d}^N}{\textrm{d}\mu^N}\I^N\left(\Big(\big(\mu_j,u(\mu_j)\prod_{i=1}^N(\mu_j-\lambda_i)\big)\Big)_{j=1}^S\right)},\]
 which, if $\I^N$ is expressed via a hierarchical basis, is equivalent to
-\begin{equation}\label{eq:pole2}
-\{\widehat{\lambda}_i\}_{i=1}^N=\argmin_{\{\lambda_i\}_{i=1}^N\subset\C}\norm{[\underbrace{0,\ldots,0}_{N},1]\big(\Phi^H\Phi\big)^{-1}\Phi^H\Big[u(\mu_j)\delta_{jj'}\Big]_{j=1,j'=1}^{S,S}\left[\prod_{i=1}^N(\mu_j-\lambda_i)\right]_{j=1}^S}^2
+\begin{align}
+\{\widehat{\lambda}_i\}_{i=1}^N=&\argmin_{\{\lambda_i\}_{i=1}^N\subset\C}\norm{[\underbrace{0,\ldots,0}_{N},1]\big(\Phi^H\Phi\big)^{-1}\Phi^H\Big[u(\mu_j)\delta_{jj'}\Big]_{j=1,j'=1}^{S,S}\left[\prod_{i=1}^N(\mu_j-\lambda_i)\right]_{j=1}^S}^2\label{eq:pole2}\\
+=&\argmin_{\{\lambda_i\}_{i=1}^N\subset\C}\norm{\mathbf{y}^H\left[q(\mu_j)\right]_{j=1}^S}^2\nonumber
+\end{align}
+(the outer square has been added for later convenience), with
+\begin{equation}\label{eq:poleadd}
+\mathbf{y}=\Big[\overline{u(\mu_j)}\delta_{jj'}\Big]_{j=1,j'=1}^{S,S}\Phi\big(\Phi^H\Phi\big)^{-1}[\underbrace{0,\ldots,0}_{N},1]^\top\in V^S.
 \end{equation}
-(the outer square has been added for later convenience). This problem is nonconvex and highly nonlinear with respect to $\{\lambda_i\}_{i=1}^N\subset\C$, and an iterative method is necessary for its solution.
+This problem is nonconvex and highly nonlinear with respect to $\{\lambda_i\}_{i=1}^N\subset\C$, and an iterative method is necessary for its solution.
 
 Two observations aid us here:
 \begin{itemize}
 \item The target functional has a relatively simple structure, allowing to compute exactly its Jacobian and Hessian with respect to the poles. Explicitly, for all $\{\lambda_i\}_{i=1}^N\subset\C$ and $k=1,\ldots,N$, the partial derivative with respect to $\lambda_k$ at $\{\lambda_i\}_{i=1}^N\subset\C$ equals
 \begin{equation}\label{eq:pole3}
--2\inner{\mathbf{y}^H\left[q(\mu_j)\right]_{j=1}^S}{\mathbf{y}^H\left[\frac{q(\mu_j)}{\mu_j-\lambda_k}\right]_{j=1}^S}
+-2\inner{\mathbf{y}^H\left[q(\mu_j)\right]_{j=1}^S}{\mathbf{y}^H\left[\frac{q(\mu_j)}{\mu_j-\lambda_k}\right]_{j=1}^S},
 \end{equation}
-where $q$ is given in \eqref{eq:pole1} and
-\[\mathbf{y}=\Big[\overline{u(\mu_j)}\delta_{jj'}\Big]_{j=1,j'=1}^{S,S}\Phi\big(\Phi^H\Phi\big)^{-1}[\underbrace{0,\ldots,0}_{N},1]^\top\in V^S\]
-does not change as the iterations proceed. Similarly, for all $\{\lambda_i\}_{i=1}^N\subset\C$ and $k,k'=1,\ldots,N$, the second order partial derivative with respect to $\lambda_k$ and $\lambda_{k'}$ at $\{\lambda_i\}_{i=1}^N\subset\C$ is
+where $\mathbf{y}$ from \eqref{eq:poleadd} does not change as the iterations proceed. Similarly, for all $\{\lambda_i\}_{i=1}^N\subset\C$ and $k,k'=1,\ldots,N$, the second order partial derivative with respect to $\lambda_k$ and $\lambda_{k'}$ at $\{\lambda_i\}_{i=1}^N\subset\C$ is
 \begin{multline}\label{eq:pole4}
 2\inner{\mathbf{y}^H\left[\frac{q(\mu_j)}{\mu_j-\lambda_{k'}}\right]_{j=1}^S}{\mathbf{y}^H\left[\frac{q(\mu_j)}{\mu_j-\lambda_k}\right]_{j=1}^S}+\\
 +2(1-\delta_{kk'})\inner{\mathbf{y}^H\left[q(\mu_j)\right]_{j=1}^S}{\mathbf{y}^H\left[\frac{q(\mu_j)}{(\mu_j-\lambda_k)(\mu_j-\lambda_{k'})}\right]_{j=1}^S}
 \end{multline}
 Hence, we can use Newton's method.
-\item Due to the iterative nature of the solution strategy, an initial guess of the optimal poles is necessary. Luckily, we can employ any of the other four methods to this aim. More precisely, we apply \code{DOMINANT} as a preliminary step, compute the roots of the resulting $\widehat{q}$, and use them as initial guess for the Newton iterations.
+\item Due to the iterative nature of the solution strategy, an initial guess of the optimal poles is necessary. Luckily, we can employ any of the other four methods to this aim. In \code{RROMPy}, we apply \code{DOMINANT} as a preliminary step, compute the roots of the resulting $\widehat{q}$, and use them as initial guess for the Newton iterations.
 \end{itemize}
 
-As stopping criterion for Newton's method, we use the (relative) norm of the increment:
-\[\sum_{i=1}^N\left|\lambda_i^{(\textrm{iter}+1)}-\lambda_i^{(\textrm{iter})}\right|^2\leq\underbrace{\textrm{tolerance}}_{\sim 10^{-10}}\sum_{i=1}^N\left|\lambda_i^{(\textrm{iter}+1)}\right|^2.\]
+As stopping criterion for Newton's method, we use the relative norm of the increment:
+\[\sum_{i=1}^N\left|\lambda_i^{(\textrm{iter}+1)}-\lambda_i^{(\textrm{iter})}\right|^2\leq\underbrace{\textrm{tolerance}}_{10^{-10}}\sum_{i=1}^N\left|\lambda_i^{(\textrm{iter}+1)}\right|^2.\]
 Since the initial guess is quite good, usually just a few (most often, 1) Newton iterations are necessary.
 
 
 \section{Minor observations}\label{sec:observations}
 \begin{itemize}
 \item For \code{BARYCENTRYC\_*}, a specific choice of polynomial basis for $\I^N$ was used to diagonalize the functional. Under the assumptions that the sample points are distinct and that $N=S-1$, one can employ the quasi-Lagrangian basis \eqref{eq:bary1} to expand $\I^N$ in the other approaches as well, thus simplifying significantly the structure of \eqref{eq:glob1}:
 \begin{equation*}
 \widehat{q}=\argmin_{\substack{q\in\mathbb{P}^N(\C;\C)\\(\star)}}\norm{\sum_{j=1}^Sq(\mu_j)u(\mu_j)\prod_{\substack{j'=1\\j'\neq j}}^S\frac1{\mu_j-\mu_{j'}}}.
 \end{equation*}
-Numerically, this has repercussions on the computation of the term $\big(\Phi^H\Phi\big)^{-1}\Phi^H$ in the Gramian \eqref{eq:poly3} and of its adjoint $\Phi\big(\Phi^H\Phi\big)^{-1}$ in the vector $\mathbf{y}$ in \eqref{eq:pole3}.
+Numerically, this has repercussions on the computation of the term $\big(\Phi^H\Phi\big)^{-1}\Phi^H$ in \eqref{eq:poly3} and of its adjoint $\Phi\big(\Phi^H\Phi\big)^{-1}$ in \eqref{eq:poleadd}.
 \item In general, \code{NORM} and \code{BARYCENTRIC\_NORM} can be expected to be more numerically stable than \code{DOMINANT} and \code{BARYCENTRIC\_AVERAGE}, respectively. This is due to the fact that the normalization is enforced in a more numerically robust fashion.
 \item If the snapshots are orthonormalized via \code{POD}\footpath{./rrompy/sampling/engines/sampling_engine_pod.py}, all the $V$-inner products (resp. norms) are recast as Euclidean inner products (resp. norms) involving the $R$ factor of the generalized ($V$-orthonormal) QR decomposition of the snapshots.
-\item If a univariate rational surrogate is built in the scope of multivariate pivoted approximation\footpath{./rrompy/reduction_methods/pivoted/*}, the rational approximant is converted into a Heaviside/nodal representation when different surrogates are combined. As such, the \code{BARYCENTRYC\_*} or \code{NODAL} approaches may be preferable to avoid extra computations, as well as additional round-off artifacts.
+\item If a univariate rational surrogate is built in the scope of multivariate pole-matching-based pivoted approximation\footpath{./rrompy/reduction_methods/pivoted/*}, the rational approximant is converted into a Heaviside/nodal representation when different surrogates are combined. As such, the \code{BARYCENTRYC\_*} or \code{NODAL} approaches may be preferable to avoid extra computations, as well as additional round-off artifacts.
 \end{itemize}
 
 \begin{thebibliography}{1}
 
 \bibitem{Pradovera}% 
  D.~Pradovera, Interpolatory rational model order reduction of parametric problems lacking uniform inf-sup stability, SIAM J. Numer. Anal. 58 (2020) 2265--2293. doi:10.1137/19M1269695.
 
 \bibitem{Klein}%
  G.~Klein, Applications of Linear Barycentric Rational Interpolation, PhD Thesis no. 1762, Universit\'e de Fribourg (2012).
 \end{thebibliography}
 
 \end{document}
\ No newline at end of file
diff --git a/rrompy/reduction_methods/__init__.py b/rrompy/reduction_methods/__init__.py
index 37b7ffa..cc4c46c 100644
--- a/rrompy/reduction_methods/__init__.py
+++ b/rrompy/reduction_methods/__init__.py
@@ -1,46 +1,42 @@
 # Copyright (C) 2018-2020 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 .standard import NearestNeighbor, RationalInterpolant, ReducedBasis
 from .standard.greedy import RationalInterpolantGreedy, ReducedBasisGreedy
 from .pivoted import (RationalInterpolantPivotedNoMatch,
-                      RationalInterpolantPivoted,
+                      RationalInterpolantPivotedPoleMatch,
                       RationalInterpolantGreedyPivotedNoMatch,
-                      RationalInterpolantGreedyPivoted)
-from .pivoted.greedy import (RationalInterpolantPivotedGreedyNoMatch,
-                             RationalInterpolantPivotedGreedy,
-                             RationalInterpolantGreedyPivotedGreedyNoMatch,
-                             RationalInterpolantGreedyPivotedGreedy)
+                      RationalInterpolantGreedyPivotedPoleMatch)
+from .pivoted.greedy import (RationalInterpolantPivotedGreedyPoleMatch,
+                             RationalInterpolantGreedyPivotedGreedyPoleMatch)
 
 __all__ = [
         'NearestNeighbor',
         'RationalInterpolant',
         'ReducedBasis',
         'RationalInterpolantGreedy',
         'ReducedBasisGreedy',
         'RationalInterpolantPivotedNoMatch',
-        'RationalInterpolantPivoted',
+        'RationalInterpolantPivotedPoleMatch',
         'RationalInterpolantGreedyPivotedNoMatch',
-        'RationalInterpolantGreedyPivoted',
-        'RationalInterpolantPivotedGreedyNoMatch',
-        'RationalInterpolantPivotedGreedy',
-        'RationalInterpolantGreedyPivotedGreedyNoMatch',
-        'RationalInterpolantGreedyPivotedGreedy'
+        'RationalInterpolantGreedyPivotedPoleMatch',
+        'RationalInterpolantPivotedGreedyPoleMatch',
+        'RationalInterpolantGreedyPivotedGreedyPoleMatch'
           ]
 
 
diff --git a/rrompy/reduction_methods/pivoted/__init__.py b/rrompy/reduction_methods/pivoted/__init__.py
index 80da5fb..ed82ffa 100644
--- a/rrompy/reduction_methods/pivoted/__init__.py
+++ b/rrompy/reduction_methods/pivoted/__init__.py
@@ -1,31 +1,31 @@
 # Copyright (C) 2018-2020 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 .rational_interpolant_pivoted import (RationalInterpolantPivotedNoMatch,
-                                           RationalInterpolantPivoted)
+                                           RationalInterpolantPivotedPoleMatch)
 from .rational_interpolant_greedy_pivoted import (RationalInterpolantGreedyPivotedNoMatch,
-                                                  RationalInterpolantGreedyPivoted)
+                                                  RationalInterpolantGreedyPivotedPoleMatch)
 
 __all__ = [
         'RationalInterpolantPivotedNoMatch',
-        'RationalInterpolantPivoted',
+        'RationalInterpolantPivotedPoleMatch',
         'RationalInterpolantGreedyPivotedNoMatch',
-        'RationalInterpolantGreedyPivoted'
+        'RationalInterpolantGreedyPivotedPoleMatch'
           ]
 
 
diff --git a/rrompy/reduction_methods/pivoted/generic_pivoted_approximant.py b/rrompy/reduction_methods/pivoted/generic_pivoted_approximant.py
index 219166b..afe3b99 100644
--- a/rrompy/reduction_methods/pivoted/generic_pivoted_approximant.py
+++ b/rrompy/reduction_methods/pivoted/generic_pivoted_approximant.py
@@ -1,784 +1,783 @@
 # Copyright (C) 2018-2020 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 os import mkdir, remove, rmdir
 import numpy as np
 from collections.abc import Iterable
 from copy import deepcopy as copy
 from rrompy.reduction_methods.base.generic_approximant import (
                                                             GenericApproximant)
 from rrompy.utilities.base.data_structures import purgeDict, getNewFilename
 from rrompy.utilities.poly_fitting.polynomial import polybases as ppb
 from rrompy.utilities.poly_fitting.radial_basis import polybases as rbpb
 from rrompy.utilities.poly_fitting.piecewise_linear import sparsekinds as sk
 from rrompy.utilities.base.types import Np2D, paramList, List, ListAny
 from rrompy.utilities.base import verbosityManager as vbMng
 from rrompy.utilities.numerical.degree import reduceDegreeN
 from rrompy.utilities.exception_manager import RROMPyException, RROMPyWarning
 from rrompy.parameter import checkParameterList
 from rrompy.utilities.parallel import poolRank, bcast
 
-__all__ = ['GenericPivotedApproximantNoMatch', 'GenericPivotedApproximant']
+__all__ = ['GenericPivotedApproximantNoMatch',
+           'GenericPivotedApproximantPoleMatch']
 
 class GenericPivotedApproximantBase(GenericApproximant):
     def __init__(self, directionPivot:ListAny, *args,
                  storeAllSamples : bool = False, **kwargs):
         self._preInit()
         if len(directionPivot) > 1:
             raise RROMPyException(("Exactly 1 pivot parameter allowed in pole "
                                    "matching."))
         from rrompy.parameter.parameter_sampling import (EmptySampler as ES,
                                                        SparseGridSampler as SG)
         self._addParametersToList(["radialDirectionalWeightsMarginal"], [1.],
                                   ["samplerPivot", "SMarginal",
                                    "samplerMarginal"],
                                   [ES(), 1, SG([[-1.], [1.]])],
                                   toBeExcluded = ["sampler"])
         self._directionPivot = directionPivot
         self.storeAllSamples = storeAllSamples
         if not hasattr(self, "_output_lvl"): self._output_lvl = []
         self._output_lvl += [1 / 2]
         super().__init__(*args, **kwargs)
         self._postInit()
 
     def setupSampling(self): super().setupSampling(False)
 
     def initializeModelData(self, datadict):
         if "directionPivot" in datadict.keys():
             from .trained_model.trained_model_pivoted_data import (
                                                        TrainedModelPivotedData)
             data = TrainedModelPivotedData(datadict["mu0"], datadict["mus"],
                                            datadict.pop("projMat"),
                                            datadict["scaleFactor"],
                                            datadict.pop("parameterMap"),
                                            datadict["directionPivot"])
             if hasattr(self.HFEngine, "_is_C_quadratic") and not (
                               hasattr(self, "matchState") and self.matchState):
                 data._is_C_quadratic = self.HFEngine._is_C_quadratic
             return (data, ["mu0", "scaleFactor", "directionPivot", "mus"])
         else:
             return super().initializeModelData(datadict)
 
     @property
     def npar(self):
         """Number of parameters."""
         if hasattr(self, "_temporaryPivot"): return self.nparPivot
         return super().npar
 
     def checkParameterListPivot(self, mu:paramList,
                                 check_if_single : bool = False) -> paramList:
         return checkParameterList(mu, self.nparPivot, check_if_single)
 
     def checkParameterListMarginal(self, mu:paramList,
                                   check_if_single : bool = False) -> paramList:
         return checkParameterList(mu, self.nparMarginal, check_if_single)
 
     def mapParameterList(self, *args, **kwargs):
         if hasattr(self, "_temporaryPivot"):
             return self.mapParameterListPivot(*args, **kwargs)
         return super().mapParameterList(*args, **kwargs)
 
     def mapParameterListPivot(self, mu:paramList, direct : str = "F",
                               idx : List[int] = None):
         if idx is None:
             idx = self.directionPivot
         else:
             idx = [self.directionPivot[j] for j in idx]
         return super().mapParameterList(mu, direct, idx)
 
     def mapParameterListMarginal(self, mu:paramList, direct : str = "F",
                                  idx : List[int] = None):
         if idx is None:
             idx = self.directionMarginal
         else:
             idx = [self.directionMarginal[j] for j in idx]
         return super().mapParameterList(mu, direct, idx)
 
     @property
     def mu0(self):
         """Value of mu0."""
         if hasattr(self, "_temporaryPivot"):
             return self.checkParameterListPivot(self._mu0(self.directionPivot))
         return self._mu0
     @mu0.setter
     def mu0(self, mu0):
         GenericApproximant.mu0.fset(self, mu0)
 
     @property
     def mus(self):
         """Value of mus. Its assignment may reset snapshots."""
         return self._mus
     @mus.setter
     def mus(self, mus):
         mus = self.checkParameterList(mus)
         musOld =  copy(self.mus) if hasattr(self, '_mus') else None
         if (musOld is None or len(mus) != len(musOld) or not mus == musOld):
             self.resetSamples()
             self._mus = mus
 
     @property
     def musMarginal(self):
         """Value of musMarginal. Its assignment may reset snapshots."""
         return self._musMarginal
     @musMarginal.setter
     def musMarginal(self, musMarginal):
         musMarginal = self.checkParameterListMarginal(musMarginal)
         if hasattr(self, '_musMarginal'):
             musMOld =  copy(self.musMarginal)
         else:
             musMOld = None
         if (musMOld is None or len(musMarginal) != len(musMOld)
                             or not musMarginal == musMOld):
             self.resetSamples()
             self._musMarginal = musMarginal
 
     @property
     def SMarginal(self):
         """Value of SMarginal."""
         return self._SMarginal
     @SMarginal.setter
     def SMarginal(self, SMarginal):
         if SMarginal <= 0:
             raise RROMPyException("SMarginal must be positive.")
         if hasattr(self, "_SMarginal") and self._SMarginal is not None:
             Sold = self.SMarginal
         else: Sold = -1
         self._SMarginal = SMarginal
         self._approxParameters["SMarginal"] = self.SMarginal
         if Sold != self.SMarginal: self.resetSamples()
 
     @property
     def radialDirectionalWeightsMarginal(self):
         """Value of radialDirectionalWeightsMarginal."""
         return self._radialDirectionalWeightsMarginal
     @radialDirectionalWeightsMarginal.setter
     def radialDirectionalWeightsMarginal(self, radialDirWeightsMarg):
         if isinstance(radialDirWeightsMarg, Iterable):
             radialDirWeightsMarg = list(radialDirWeightsMarg)
         else:
             radialDirWeightsMarg = [radialDirWeightsMarg]
         self._radialDirectionalWeightsMarginal = radialDirWeightsMarg
         self._approxParameters["radialDirectionalWeightsMarginal"] = (
                                          self.radialDirectionalWeightsMarginal)
 
     @property
     def directionPivot(self):
         """Value of directionPivot. Its assignment may reset snapshots."""
         return self._directionPivot
     @directionPivot.setter
     def directionPivot(self, directionPivot):
         if hasattr(self, '_directionPivot'):
             directionPivotOld = copy(self.directionPivot)
         else:
             directionPivotOld = None
         if (directionPivotOld is None
          or len(directionPivot) != len(directionPivotOld)
          or not directionPivot == directionPivotOld):
             self.resetSamples()
             self._directionPivot = directionPivot
 
     @property
     def directionMarginal(self):
         return [x for x in range(self.HFEngine.npar) \
                                                if x not in self.directionPivot]
 
     @property
     def nparPivot(self):
         return len(self.directionPivot)
 
     @property
     def nparMarginal(self):
         return self.npar - self.nparPivot
 
     @property
     def muBounds(self):
         """Value of muBounds."""
         return self.samplerPivot.lims
 
     @property
     def muBoundsMarginal(self):
         """Value of muBoundsMarginal."""
         return self.samplerMarginal.lims
 
     @property
     def sampler(self):
         """Proxy of samplerPivot."""
         return self._samplerPivot
 
     @property
     def samplerPivot(self):
         """Value of samplerPivot."""
         return self._samplerPivot
     @samplerPivot.setter
     def samplerPivot(self, samplerPivot):
         if 'generatePoints' not in dir(samplerPivot):
             raise RROMPyException("Pivot sampler type not recognized.")
         if hasattr(self, '_samplerPivot') and self._samplerPivot is not None:
             samplerOld = self.samplerPivot
         self._samplerPivot = samplerPivot
         self._approxParameters["samplerPivot"] = self.samplerPivot
         if not 'samplerOld' in locals() or samplerOld != self.samplerPivot:
             self.resetSamples()
     
     @property
     def samplerMarginal(self):
         """Value of samplerMarginal."""
         return self._samplerMarginal
     @samplerMarginal.setter
     def samplerMarginal(self, samplerMarginal):
         if 'generatePoints' not in dir(samplerMarginal):
             raise RROMPyException("Marginal sampler type not recognized.")
         if (hasattr(self, '_samplerMarginal')
         and self._samplerMarginal is not None):
             samplerOld = self.samplerMarginal
         self._samplerMarginal = samplerMarginal
         self._approxParameters["samplerMarginal"] = self.samplerMarginal
         if not 'samplerOld' in locals() or samplerOld != self.samplerMarginal:
             self.resetSamples()
     
     def computeScaleFactor(self):
         """Compute parameter rescaling factor."""
         self.scaleFactorPivot = .5 * np.abs((
                               self.mapParameterListPivot(self.muBounds[0])
                             - self.mapParameterListPivot(self.muBounds[1]))[0])
         self.scaleFactorMarginal = .5 * np.abs((
                    self.mapParameterListMarginal(self.muBoundsMarginal[0])
                  - self.mapParameterListMarginal(self.muBoundsMarginal[1]))[0])
         self.scaleFactor = np.empty(self.npar)
         self.scaleFactor[self.directionPivot] = self.scaleFactorPivot
         self.scaleFactor[self.directionMarginal] = self.scaleFactorMarginal
 
     def _setupTrainedModel(self, pMat:Np2D, pMatUpdate : bool = False,
                            pMatOld : Np2D = None, forceNew : bool = False):
         if forceNew or self.trainedModel is None:
             self.trainedModel = self.tModelType()
             self.trainedModel.verbosity = self.verbosity
             self.trainedModel.timestamp = self.timestamp
             datadict = {"mu0": self.mu0, "mus": copy(self.mus),
                         "projMat": pMat, "scaleFactor": self.scaleFactor,
                         "parameterMap": self.HFEngine.parameterMap,
                         "directionPivot": self.directionPivot}
             self.trainedModel.data = self.initializeModelData(datadict)[0]
         else:
             self.trainedModel = self.trainedModel
             if pMatUpdate:
                 self.trainedModel.data.projMat = np.hstack(
                                         (self.trainedModel.data.projMat, pMat))
             else:
                 self.trainedModel.data.projMat = copy(pMat)
             self.trainedModel.data.mus = copy(self.mus)
         self.trainedModel.data.musMarginal = copy(self.musMarginal)
 
     def normApprox(self, mu:paramList) -> float:
         _PODOld, self._POD = self.POD, 0
         result = super().normApprox(mu)
         self._POD = _PODOld
         return result
 
     @property
     def storedSamplesFilenames(self) -> List[str]:
         if not hasattr(self, "_sampleBaseFilename"): return []
         return [self._sampleBaseFilename
               + "{}_{}.pkl" .format(idx + 1, self.name())
                                        for idx in range(len(self.musMarginal))]
 
     def purgeStoredSamples(self):
         if not hasattr(self, "_sampleBaseFilename"): return
         for file in self.storedSamplesFilenames: remove(file)
         rmdir(self._sampleBaseFilename[: -8])
 
     def storeSamples(self, idx : int = None):
         """Store samples to file."""
         if not hasattr(self, "_sampleBaseFilename"):
             filenameBase = None
             if poolRank() == 0:
                 foldername = getNewFilename(self.name(), "samples")
                 mkdir(foldername)
                 filenameBase = foldername + "/sample_"
             self._sampleBaseFilename = bcast(filenameBase, force = True)
         if idx is not None:
             super().storeSamples(self._sampleBaseFilename + str(idx + 1),
                                  False)
 
     def loadTrainedModel(self, filename:str):
         """Load trained reduced model from file."""
         super().loadTrainedModel(filename)
         self._musMarginal = self.trainedModel.data.musMarginal
 
 class GenericPivotedApproximantNoMatch(GenericPivotedApproximantBase):
     """
     ROM pivoted approximant (without pole matching) computation for parametric
         problems (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': kind of snapshots orthogonalization; allowed values
                 include 0, 1/2, and 1; defaults to 1, i.e. POD;
             - 'scaleFactorDer': scaling factors for derivative computation;
                 defaults to 'AUTO';
             - '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;
             - 'radialDirectionalWeightsMarginal': radial basis weights for
                 marginal interpolant; defaults to 1.
             Defaults to empty dict.
         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': kind of snapshots orthogonalization;
             - 'scaleFactorDer': scaling factors for derivative computation;
             - '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.
         verbosity: Verbosity level.
         POD: Kind of snapshots orthogonalization.
         scaleFactorDer: Scaling factors for derivative computation.
         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.
         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.
     """
 
     @property
     def tModelType(self):
         from .trained_model.trained_model_pivoted_rational_nomatch import (
                                             TrainedModelPivotedRationalNoMatch)
         return TrainedModelPivotedRationalNoMatch
 
     def _finalizeMarginalization(self):
         self.trainedModel.setupMarginalInterp(
                                        [self.radialDirectionalWeightsMarginal])
         self.trainedModel.data.approxParameters = copy(self.approxParameters)
 
-    def _poleMatching(self):
-        vbMng(self, "INIT", "Compressing poles.", 10)
-        self.trainedModel.initializeFromRational()
-        vbMng(self, "DEL", "Done compressing poles.", 10)
+    def _preliminaryMarginalFinalization(self):
+        pass
 
-class GenericPivotedApproximant(GenericPivotedApproximantBase):
+class GenericPivotedApproximantPoleMatch(GenericPivotedApproximantBase):
     """
     ROM pivoted approximant (with pole matching) computation for parametric
         problems (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': kind of snapshots orthogonalization; allowed values
                 include 0, 1/2, and 1; defaults to 1, i.e. POD;
             - 'scaleFactorDer': scaling factors for derivative computation;
                 defaults to 'AUTO';
             - 'matchState': whether to match the system state rather than the
                 system output; defaults to False;
             - 'matchingWeight': weight for pole matching optimization; defaults
                 to 1;
             - 'sharedRatio': required ratio of marginal points to share
                 resonance; defaults to 1.;
             - '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;
             - '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.
             - 'radialDirectionalWeightsMarginal': radial basis weights for
                 marginal interpolant; defaults to 1.
             Defaults to empty dict.
         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': kind of snapshots orthogonalization;
             - 'scaleFactorDer': scaling factors for derivative computation;
             - 'matchState': whether to match the system state rather than the
                 system output;
             - 'matchingWeight': weight for pole matching optimization;
             - 'sharedRatio': required ratio of marginal points to share
                 resonance;
             - '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.
             - '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.
         verbosity: Verbosity level.
         POD: Kind of snapshots orthogonalization.
         scaleFactorDer: Scaling factors for derivative computation.
         matchState: Whether to match the system state rather than the system
             output.
         matchingWeight: Weight for pole matching optimization.
         sharedRatio: Required ratio of marginal points to share resonance.
         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.
         polybasisMarginal: Type of polynomial basis for marginal interpolation.
         paramsMarginal: Dictionary of parameters for marginal interpolation.
         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 __init__(self, *args, **kwargs):
         self._preInit()
         self._addParametersToList(["matchState", "matchingWeight",
                                    "sharedRatio", "polybasisMarginal",
                                    "paramsMarginal"],
                                   [False, 1., 1., "MONOMIAL", {}])
         self.parameterMarginalList = ["MMarginal", "nNeighborsMarginal",
                                       "polydegreetypeMarginal",
                                       "interpRcondMarginal",
                                       "radialDirectionalWeightsMarginalAdapt"]
         super().__init__(*args, **kwargs)
         self._postInit()
 
     @property
     def tModelType(self):
-        from .trained_model.trained_model_pivoted_rational import (
-                                                   TrainedModelPivotedRational)
-        return TrainedModelPivotedRational
+        from .trained_model.trained_model_pivoted_rational_polematch import (
+                                          TrainedModelPivotedRationalPoleMatch)
+        return TrainedModelPivotedRationalPoleMatch
 
     @property
     def matchState(self):
         """Value of matchState."""
         return self._matchState
     @matchState.setter
     def matchState(self, matchState):
         self._matchState = matchState
         self._approxParameters["matchState"] = self.matchState
 
     @property
     def matchingWeight(self):
         """Value of matchingWeight."""
         return self._matchingWeight
     @matchingWeight.setter
     def matchingWeight(self, matchingWeight):
         self._matchingWeight = matchingWeight
         self._approxParameters["matchingWeight"] = self.matchingWeight
 
     @property
     def sharedRatio(self):
         """Value of sharedRatio."""
         return self._sharedRatio
     @sharedRatio.setter
     def sharedRatio(self, sharedRatio):
         if sharedRatio > 1.:
             RROMPyWarning("Shared ratio too large. Clipping to 1.")
             sharedRatio = 1.
         elif sharedRatio < 0.:
             RROMPyWarning("Shared ratio too small. Clipping to 0.")
             sharedRatio = 0.
         self._sharedRatio = sharedRatio
         self._approxParameters["sharedRatio"] = self.sharedRatio
 
     @property
     def polybasisMarginal(self):
         """Value of polybasisMarginal."""
         return self._polybasisMarginal
     @polybasisMarginal.setter
     def polybasisMarginal(self, polybasisMarginal):
         try:
             polybasisMarginal = polybasisMarginal.upper().strip().replace(" ",
                                                                           "")
             if polybasisMarginal not in ppb + rbpb + ["NEARESTNEIGHBOR"] + sk:
                 raise RROMPyException(
                                "Prescribed marginal polybasis not recognized.")
             self._polybasisMarginal = polybasisMarginal
         except:
             RROMPyWarning(("Prescribed marginal polybasis not recognized. "
                            "Overriding to 'MONOMIAL'."))
             self._polybasisMarginal = "MONOMIAL"
         self._approxParameters["polybasisMarginal"] = self.polybasisMarginal
 
     @property
     def paramsMarginal(self):
         """Value of paramsMarginal."""
         return self._paramsMarginal
     @paramsMarginal.setter
     def paramsMarginal(self, paramsMarginal):
         paramsMarginal = purgeDict(paramsMarginal, self.parameterMarginalList,
                                    dictname = self.name() + ".paramsMarginal",
                                    baselevel = 1)
         keyList = list(paramsMarginal.keys())
         if not hasattr(self, "_paramsMarginal"): self._paramsMarginal = {}
         
         if "MMarginal" in keyList:
             MMarg = paramsMarginal["MMarginal"]
         elif ("MMarginal" in self.paramsMarginal
           and not hasattr(self, "_MMarginal_isauto")):
             MMarg = self.paramsMarginal["MMarginal"]
         else:
             MMarg = "AUTO"
         if isinstance(MMarg, str):
             MMarg = MMarg.strip().replace(" ","")
             if "-" not in MMarg: MMarg = MMarg + "-0"
             self._MMarginal_isauto = True
             self._MMarginal_shift = int(MMarg.split("-")[-1])
             MMarg = 0
         if MMarg < 0:
             raise RROMPyException("MMarginal must be non-negative.")
         self._paramsMarginal["MMarginal"] = MMarg
         
         if "nNeighborsMarginal" in keyList:
             self._paramsMarginal["nNeighborsMarginal"] = max(1,
                                           paramsMarginal["nNeighborsMarginal"])
         elif "nNeighborsMarginal" not in self.paramsMarginal:
             self._paramsMarginal["nNeighborsMarginal"] = 1
         
         if "polydegreetypeMarginal" in keyList:
             try:
                 polydegtypeM = paramsMarginal["polydegreetypeMarginal"]\
                                                .upper().strip().replace(" ","")
                 if polydegtypeM not in ["TOTAL", "FULL"]: 
                     raise RROMPyException(("Prescribed polydegreetypeMarginal "
                                            "not recognized."))
                 self._paramsMarginal["polydegreetypeMarginal"] = polydegtypeM
             except:
                 RROMPyWarning(("Prescribed polydegreetypeMarginal not "
                                "recognized. Overriding to 'TOTAL'."))
                 self._paramsMarginal["polydegreetypeMarginal"] = "TOTAL"
         elif "polydegreetypeMarginal" not in self.paramsMarginal:
             self._paramsMarginal["polydegreetypeMarginal"] = "TOTAL"
 
         if "interpRcondMarginal" in keyList:
             self._paramsMarginal["interpRcondMarginal"] = (
                                          paramsMarginal["interpRcondMarginal"])
         elif "interpRcondMarginal" not in self.paramsMarginal:
             self._paramsMarginal["interpRcondMarginal"] = -1
 
         if "radialDirectionalWeightsMarginalAdapt" in keyList:
             self._paramsMarginal["radialDirectionalWeightsMarginalAdapt"] = (
                        paramsMarginal["radialDirectionalWeightsMarginalAdapt"])
         elif "radialDirectionalWeightsMarginalAdapt" not in self.paramsMarginal:
             self._paramsMarginal["radialDirectionalWeightsMarginalAdapt"] = [
                                                                       -1., -1.]
         self._approxParameters["paramsMarginal"] = self.paramsMarginal
 
     def _setMMarginalAuto(self):
         if (self.polybasisMarginal not in ppb + rbpb
          or "MMarginal" not in self.paramsMarginal
          or "polydegreetypeMarginal" not in self.paramsMarginal):
             raise RROMPyException(("Cannot set MMarginal if "
                                    "polybasisMarginal does not allow it."))
         self.paramsMarginal["MMarginal"] = max(0, reduceDegreeN(
                                  len(self.musMarginal), len(self.musMarginal),
                                  self.nparMarginal,
                                  self.paramsMarginal["polydegreetypeMarginal"])
                                                 - self._MMarginal_shift)
         vbMng(self, "MAIN", ("Automatically setting MMarginal to {}.").format(
                                          self.paramsMarginal["MMarginal"]), 25)
 
     def purgeparamsMarginal(self):
         self.paramsMarginal = {}
         paramsMbadkeys = []
         if self.polybasisMarginal in ppb + rbpb + sk:
             paramsMbadkeys += ["nNeighborsMarginal"]
         if self.polybasisMarginal not in rbpb:
             paramsMbadkeys += ["radialDirectionalWeightsMarginalAdapt"]
         if self.polybasisMarginal in ["NEARESTNEIGHBOR"] + sk:
             paramsMbadkeys += ["MMarginal", "polydegreetypeMarginal",
                                "interpRcondMarginal"]
             if hasattr(self, "_MMarginal_isauto"): del self._MMarginal_isauto
             if hasattr(self, "_MMarginal_shift"): del self._MMarginal_shift
         for key in paramsMbadkeys:
             if key in self._paramsMarginal: del self._paramsMarginal[key]
         self._approxParameters["paramsMarginal"] = self.paramsMarginal
 
     def _finalizeMarginalization(self):
         vbMng(self, "INIT", "Checking shared ratio.", 10)
         msg = self.trainedModel.checkSharedRatio(self.sharedRatio)
         vbMng(self, "DEL", "Done checking." + msg, 10)
         if self.polybasisMarginal in rbpb + ["NEARESTNEIGHBOR"]:
             self.computeScaleFactor()
             rDWMEff = np.array([w * f for w, f in zip(
                                          self.radialDirectionalWeightsMarginal,
                                          self.scaleFactorMarginal)])
         if self.polybasisMarginal in ppb + rbpb + sk:
             interpPars = [self.polybasisMarginal]
             if self.polybasisMarginal in ppb + rbpb:
                 if self.polybasisMarginal in rbpb: interpPars += [rDWMEff]
                 interpPars += [self.verbosity >= 5,
                       self.paramsMarginal["polydegreetypeMarginal"] == "TOTAL"]
                 if self.polybasisMarginal in ppb:
                     interpPars += [{}]
                 else: # if self.polybasisMarginal in rbpb:
                     interpPars += [{"optimizeScalingBounds":self.paramsMarginal[
                                      "radialDirectionalWeightsMarginalAdapt"]}]
                 interpPars += [
                           {"rcond":self.paramsMarginal["interpRcondMarginal"]}]
                 extraPar = hasattr(self, "_MMarginal_isauto")
             else: # if self.polybasisMarginal in sk:
                 idxEff = [x for x in range(self.samplerMarginal.npoints)
                                   if not hasattr(self.trainedModel, "_idxExcl")
                                   or x not in self.trainedModel._idxExcl]
                 extraPar = self.samplerMarginal.depth[idxEff]
         else: # if self.polybasisMarginal == "NEARESTNEIGHBOR":
             interpPars = [self.paramsMarginal["nNeighborsMarginal"], rDWMEff]
             extraPar = None
         self.trainedModel.setupMarginalInterp(self, interpPars, extraPar)
         self.trainedModel.data.approxParameters = copy(self.approxParameters)
 
-    def _poleMatching(self):
+    def _preliminaryMarginalFinalization(self):
         vbMng(self, "INIT", "Compressing and matching poles.", 10)
         self.trainedModel.initializeFromRational(self.matchingWeight,
                                                  self.HFEngine,
                                                  self.matchState)
         vbMng(self, "DEL", "Done compressing and matching poles.", 10)
 
     def _postApplyC(self):
         if self.POD == 1 and not (
                         hasattr(self.HFEngine.C, "is_mu_independent")
                     and self.HFEngine.C.is_mu_independent in self._output_lvl):
             raise RROMPyException(("Cannot apply mu-dependent C to "
                                    "orthonormalized samples."))
         applyCglob = (hasattr(self.HFEngine, "_is_C_quadratic")
                   and self.HFEngine._is_C_quadratic)
         vbMng(self, "INIT", "Extracting system output from state.", 35)
         if applyCglob:
             pMat, dirM = None, self.trainedModel.data.directionMarginal
             for muM in self.trainedModel.data.musMarginal:
                 idx = np.where([np.allclose(mu(dirM)[0], muM[0])
                                       for mu in self.trainedModel.data.mus])[0]
                 pMati = self.trainedModel.data.projMat[:, idx]
                 musi = self.trainedModel.data.mus[idx]
                 pMati = self.HFEngine.applyC(pMati, musi)
                 if pMat is None:
                     pMat = np.array(pMati)
                 else:
                     pMat = np.append(pMat, pMati, axis = 1)
         else:
             pMat = None
             for j, mu in enumerate(musi):
                 pMij = self.trainedModel.data.projMat[:, j]
                 pMij = np.expand_dims(self.HFEngine.applyC(pMij, mu), -1)
                 if pMati is None:
                     pMat = np.array(pMij)
                 else:
                     pMat = np.append(pMat, pMij, axis = 1)
         vbMng(self, "DEL", "Done extracting system output.", 35)
         self.trainedModel.data.projMat = pMat
         if hasattr(self.HFEngine, "_is_C_quadratic"):
             self.trainedModel.data._is_C_quadratic = (
                                                  self.HFEngine._is_C_quadratic)
 
     @abstractmethod
     def setupApprox(self, *args, **kwargs) -> int:
         if self.checkComputedApprox(): return -1
         self.purgeparamsMarginal()
         setupOK = super().setupApprox(*args, **kwargs)
         if self.matchState: self._postApplyC()
         return setupOK
diff --git a/rrompy/reduction_methods/pivoted/greedy/__init__.py b/rrompy/reduction_methods/pivoted/greedy/__init__.py
index 65d9da4..80b0c7c 100644
--- a/rrompy/reduction_methods/pivoted/greedy/__init__.py
+++ b/rrompy/reduction_methods/pivoted/greedy/__init__.py
@@ -1,31 +1,27 @@
 # Copyright (C) 2018-2020 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 .rational_interpolant_pivoted_greedy import (RationalInterpolantPivotedGreedyNoMatch,
-                                                  RationalInterpolantPivotedGreedy)
-from .rational_interpolant_greedy_pivoted_greedy import (RationalInterpolantGreedyPivotedGreedyNoMatch,
-                                                         RationalInterpolantGreedyPivotedGreedy)
+from .rational_interpolant_pivoted_greedy import RationalInterpolantPivotedGreedyPoleMatch
+from .rational_interpolant_greedy_pivoted_greedy import RationalInterpolantGreedyPivotedGreedyPoleMatch
 
 __all__ = [
-        'RationalInterpolantPivotedGreedyNoMatch',
-        'RationalInterpolantPivotedGreedy',
-        'RationalInterpolantGreedyPivotedGreedyNoMatch',
-        'RationalInterpolantGreedyPivotedGreedy'
+        'RationalInterpolantPivotedGreedyPoleMatch',
+        'RationalInterpolantGreedyPivotedGreedyPoleMatch'
           ]
 
 
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 d837516..3790dbf 100644
--- a/rrompy/reduction_methods/pivoted/greedy/generic_pivoted_greedy_approximant.py
+++ b/rrompy/reduction_methods/pivoted/greedy/generic_pivoted_greedy_approximant.py
@@ -1,767 +1,654 @@
 # Copyright (C) 2018-2020 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 collections.abc import Iterable
 from matplotlib import pyplot as plt
 from rrompy.reduction_methods.pivoted.generic_pivoted_approximant import (
-                                              GenericPivotedApproximantBase,
-                                              GenericPivotedApproximantNoMatch,
-                                              GenericPivotedApproximant)
+                                            GenericPivotedApproximantBase,
+                                            GenericPivotedApproximantPoleMatch)
 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,
                                                       buildResiduesForDistance)
 from rrompy.utilities.numerical.point_distances import chordalMetricAngleMatrix
 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']
+__all__ = ['GenericPivotedGreedyApproximantPoleMatch']
 
 class GenericPivotedGreedyApproximantBase(GenericPivotedApproximantBase):
     _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 isinstance(scaleFactorDer, Iterable):
             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), :]
             if (hasattr(self.trainedModel.data, "_is_C_quadratic")
             and self.trainedModel.data._is_C_quadratic):
                 self.trainedModel._setupQuadMapping()
                 projMapping = self.quad_mapping
                 projMappingReal = self.data._is_C_quadratic == 2
             else:
                 projMapping, projMappingReal = None, False
             resEx = buildResiduesForDistance(resEx,
                                              self.trainedModel.data.projMat,
                                              0, projMapping, projMappingReal)
             resAp = buildResiduesForDistance(resAp,
                                              self.trainedModel.data.projMat,
                                              0, projMapping, projMappingReal)
         else:
             resAp = None
         dist = chordalMetricAngleMatrix(polesEx, polesAp,
                                         self.matchingWeightError, resEx, resAp,
                                         self.HFEngine, False)
         pmR, pmC = pointMatching(dist)
         return np.mean(dist[pmR, pmC])
     
     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 -= 25
             self.trainedModel.verbosity -= 25
             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 += 25
             self.trainedModel.verbosity += 25
         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 == "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()
         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 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)
                 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()
+        self._preliminaryMarginalFinalization()
         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)[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 -= 10
         self.samplingEngine.verbosity -= 10
         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 += 10
         self.samplingEngine.verbosity += 10
 
     def _localPivotedResult(self, pMat:Np2D, req:ListAny, emptyCores:ListAny,
                             mus:Np2D) -> Tuple[Np2D, ListAny, Np2D]:
         pMati = self.samplingEngine.projectionMatrix
         musi = self.samplingEngine.mus
         if not hasattr(self, "matchState") or not self.matchState:
             if self.POD == 1 and not (
                 hasattr(self.HFEngine.C, "is_mu_independent")
             and self.HFEngine.C.is_mu_independent in self._output_lvl):
                 raise RROMPyException(("Cannot apply mu-dependent C "
                                        "to orthonormalized samples."))
             vbMng(self, "INIT", "Extracting system output from state.", 35)
             if (hasattr(self.HFEngine, "_is_C_quadratic")
             and self.HFEngine._is_C_quadratic):
                 pMati = self.HFEngine.applyC(pMati, musi)
             else:
                 pMatiEff = None
                 for j, mu in enumerate(musi):
                     pMij = np.expand_dims(self.HFEngine.applyC(pMati[:, j],
                                                                mu), -1)
                     if pMatiEff is None:
                         pMatiEff = np.array(pMij)
                     else:
                         pMatiEff = np.append(pMatiEff, pMij, axis = 1)
                 pMati = pMatiEff
             vbMng(self, "DEL", "Done extracting system output.", 35)
         if pMat is None:
             mus = copy(musi.data)
             pMat = copy(pMati)
             if masterCore():
                 for dest in emptyCores:
                     req += [isend((len(pMat), pMat.dtype, mus.dtype),
                                   dest = dest, tag = dest)]
         else:
             mus = np.vstack((mus, musi.data))
             pMat = np.hstack((pMat, pMati))
         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 == "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(
+class GenericPivotedGreedyApproximantPoleMatch(
                                            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': kind of snapshots orthogonalization; allowed values
-                include 0, 1/2, and 1; defaults to 1, i.e. POD;
-            - '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 '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.
-        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': kind of snapshots orthogonalization;
-            - '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.
-        verbosity: Verbosity level.
-        POD: Kind of snapshots orthogonalization.
-        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):
+                                           GenericPivotedApproximantPoleMatch):
     """
     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': kind of snapshots orthogonalization; allowed values
                 include 0, 1/2, and 1; defaults to 1, i.e. POD;
             - 'scaleFactorDer': scaling factors for derivative computation;
                 defaults to 'AUTO';
             - 'matchState': whether to match the system state rather than the
                 system output; defaults to False;
             - 'matchingWeight': weight for pole matching optimization; defaults
                 to 1;
             - '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 '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.
         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': kind of snapshots orthogonalization;
             - 'scaleFactorDer': scaling factors for derivative computation;
             - 'matchState': whether to match the system state rather than the
                 system output;
             - 'matchingWeight': weight 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.
         verbosity: Verbosity level.
         POD: Kind of snapshots orthogonalization.
         scaleFactorDer: Scaling factors for derivative computation.
         matchState: Whether to match the system state rather than the system
             output.
         matchingWeight: Weight 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.HFEngine, False)
-        vbMng(self, "DEL", "Done compressing and matching poles.", 10)
-
     def _updateTrainedModelMarginalSamples(self, idx : ListAny = []):
         self.trainedModel.updateEffectiveSamples(idx, self.matchingWeight,
                                                  self.HFEngine, False)
 
     def setupApprox(self, *args, **kwargs) -> int:
         if self.checkComputedApprox(): return -1
         self.purgeparamsMarginal()
         _polybasisMarginal = self.polybasisMarginal
         self._polybasisMarginal = ("PIECEWISE_LINEAR_"
                                  + self.samplerMarginal.kind)
         setupOK = super().setupApprox(*args, **kwargs)
         self._polybasisMarginal = _polybasisMarginal
         self._finalizeMarginalization()
         if self.matchState: self._postApplyC()
         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 7dac0b1..d78a5d9 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,357 @@
 #Copyright (C) 2018-2020 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)
+                                      GenericPivotedGreedyApproximantBase,
+                                      GenericPivotedGreedyApproximantPoleMatch)
 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)
+                                     RationalInterpolantGreedyPivotedPoleMatch)
 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']
+__all__ = ['RationalInterpolantGreedyPivotedGreedyPoleMatch']
 
 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 -= 10
                 RationalInterpolantGreedy.setupApprox(self, self._plotEstPivot)
                 self.verbosity += 5
                 self.samplingEngine.verbosity += 10
                 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(
+class RationalInterpolantGreedyPivotedGreedyPoleMatch(
                                     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': kind of snapshots orthogonalization; allowed values
-                include 0, 1/2, and 1; defaults to 1, i.e. POD;
-            - '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 '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',
-                'BARYCENTRIC_NORM', and 'BARYCENTRIC[_AVERAGE]' (check pdf in
-                main folder for meaning); 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.
-        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': kind of snapshots orthogonalization;
-            - '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.
-        verbosity: Verbosity level.
-        POD: Kind of snapshots orthogonalization.
-        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):
+                                    GenericPivotedGreedyApproximantPoleMatch,
+                                    RationalInterpolantGreedyPivotedPoleMatch):
     """
     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': kind of snapshots orthogonalization; allowed values
                 include 0, 1/2, and 1; defaults to 1, i.e. POD;
             - 'scaleFactorDer': scaling factors for derivative computation;
                 defaults to 'AUTO';
             - 'matchState': whether to match the system state rather than the
                 system output; defaults to False;
             - 'matchingWeight': weight for pole matching optimization; defaults
                 to 1;
             - '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 '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',
                 'BARYCENTRIC_NORM', and 'BARYCENTRIC[_AVERAGE]' (check pdf in
                 main folder for meaning); 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.
         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': kind of snapshots orthogonalization;
             - 'scaleFactorDer': scaling factors for derivative computation;
             - 'matchState': whether to match the system state rather than the
                 system output;
             - 'matchingWeight': weight 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.
         verbosity: Verbosity level.
         POD: Kind of snapshots orthogonalization.
         scaleFactorDer: Scaling factors for derivative computation.
         matchState: Whether to match the system state rather than the system
             output.
         matchingWeight: Weight 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 91c5eda..e33f8b0 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,287 @@
 # Copyright (C) 2018-2020 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)
+                                      GenericPivotedGreedyApproximantBase,
+                                      GenericPivotedGreedyApproximantPoleMatch)
 from rrompy.reduction_methods.standard import RationalInterpolant
 from rrompy.reduction_methods.pivoted import (
-                                             RationalInterpolantPivotedNoMatch,
-                                             RationalInterpolantPivoted)
+                                           RationalInterpolantPivotedPoleMatch)
 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']
+__all__ = ['RationalInterpolantPivotedGreedyPoleMatch']
 
 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': kind of snapshots orthogonalization; allowed values
-                include 0, 1/2, and 1; defaults to 1, i.e. POD;
-            - '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 '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',
-                'BARYCENTRIC_NORM', and 'BARYCENTRIC[_AVERAGE]' (check pdf in
-                main folder for meaning); 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.
-        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': kind of snapshots orthogonalization;
-            - '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.
-        verbosity: Verbosity level.
-        POD: Kind of snapshots orthogonalization.
-        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):
+class RationalInterpolantPivotedGreedyPoleMatch(
+                                      RationalInterpolantPivotedGreedyBase,
+                                      GenericPivotedGreedyApproximantPoleMatch,
+                                      RationalInterpolantPivotedPoleMatch):
     """
     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': kind of snapshots orthogonalization; allowed values
                 include 0, 1/2, and 1; defaults to 1, i.e. POD;
             - 'scaleFactorDer': scaling factors for derivative computation;
                 defaults to 'AUTO';
             - 'matchState': whether to match the system state rather than the
                 system output; defaults to False;
             - 'matchingWeight': weight for pole matching optimization; defaults
                 to 1;
             - '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 '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',
                 'BARYCENTRIC_NORM', and 'BARYCENTRIC[_AVERAGE]' (check pdf in
                 main folder for meaning); 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.
         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': kind of snapshots orthogonalization;
             - 'scaleFactorDer': scaling factors for derivative computation;
             - 'matchState': whether to match the system state rather than the
                 system output;
             - 'matchingWeight': weight 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.
         verbosity: Verbosity level.
         POD: Kind of snapshots orthogonalization.
         scaleFactorDer: Scaling factors for derivative computation.
         matchState: Whether to match the system state rather than the system
             output.
         matchingWeight: Weight 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/rrompy/reduction_methods/pivoted/rational_interpolant_greedy_pivoted.py b/rrompy/reduction_methods/pivoted/rational_interpolant_greedy_pivoted.py
index 4267654..bd6b08f 100644
--- a/rrompy/reduction_methods/pivoted/rational_interpolant_greedy_pivoted.py
+++ b/rrompy/reduction_methods/pivoted/rational_interpolant_greedy_pivoted.py
@@ -1,571 +1,572 @@
 # Copyright (C) 2018-2020 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_approximant import (GenericPivotedApproximantBase,
                                           GenericPivotedApproximantNoMatch,
-                                          GenericPivotedApproximant)
+                                          GenericPivotedApproximantPoleMatch)
 from .gather_pivoted_approximant import gatherPivotedApproximant
 from rrompy.reduction_methods.standard.greedy.rational_interpolant_greedy \
                                                import RationalInterpolantGreedy
 from rrompy.reduction_methods.standard.greedy.generic_greedy_approximant \
                                                             import pruneSamples
 from rrompy.utilities.base.types import Np1D
 from rrompy.utilities.base import verbosityManager as vbMng
 from rrompy.utilities.poly_fitting.polynomial import polyvander as pv
 from rrompy.utilities.exception_manager import RROMPyException, RROMPyAssert
 from rrompy.parameter import emptyParameterList, parameterList
 from rrompy.utilities.parallel import poolRank, indicesScatter, isend, recv
 
 __all__ = ['RationalInterpolantGreedyPivotedNoMatch',
-           'RationalInterpolantGreedyPivoted']
+           'RationalInterpolantGreedyPivotedPoleMatch']
 
 class RationalInterpolantGreedyPivotedBase(GenericPivotedApproximantBase,
                                            RationalInterpolantGreedy):
     def __init__(self, *args, **kwargs):
         self._preInit()
         super().__init__(*args, **kwargs)
         self._ignoreResidues = self.nparPivot > 1
         self._postInit()
 
     @property
     def tModelType(self):
         if hasattr(self, "_temporaryPivot"):
             return RationalInterpolantGreedy.tModelType.fget(self)
         return super().tModelType
 
     def _polyvanderAuxiliary(self, mus, deg, *args):
         degEff = [0] * self.npar
         degEff[self.directionPivot[0]] = deg
         return pv(mus, degEff, *args)
 
     def _marginalizeMiscellanea(self, forward:bool):
         if forward:
             self._m_selfmus = copy(self.mus)
             self._m_HFEparameterMap = copy(self.HFEngine.parameterMap)
             self._mus = self.checkParameterListPivot(
                                                  self.mus(self.directionPivot))
             self.HFEngine.parameterMap = {
                 "F": [self.HFEngine.parameterMap["F"][self.directionPivot[0]]],
                 "B": [self.HFEngine.parameterMap["B"][self.directionPivot[0]]]}
         else:
             self._mus = self._m_selfmus
             self.HFEngine.parameterMap = self._m_HFEparameterMap
             del self._m_selfmus, self._m_HFEparameterMap
 
     def _marginalizeTrainedModel(self, forward:bool):
         if forward:
             del self._temporaryPivot
             self.trainedModel.data.mu0 = self.mu0
             self.trainedModel.data.scaleFactor = [1.] * self.npar
             self.trainedModel.data.scaleFactor[self.directionPivot[0]] = (
                                                            self.scaleFactor[0])
             self.trainedModel.data.parameterMap = self.HFEngine.parameterMap
             self._m_musUniqueCN = copy(self._musUniqueCN)
             musUniqueCNAux = np.zeros((self.S, self.npar),
                                       dtype = self._musUniqueCN.dtype)
             musUniqueCNAux[:, self.directionPivot[0]] = self._musUniqueCN(0)
             self._musUniqueCN = self.checkParameterList(musUniqueCNAux)
             self._m_derIdxs = copy(self._derIdxs)
             for j in range(len(self._derIdxs)):
                 for l in range(len(self._derIdxs[j])):
                     derjl = self._derIdxs[j][l][0]
                     self._derIdxs[j][l] = [0] * self.npar
                     self._derIdxs[j][l][self.directionPivot[0]] = derjl
             self.trainedModel.data.Q._dirPivot = self.directionPivot[0]
             self.trainedModel.data.P._dirPivot = self.directionPivot[0]
             # tell greedy error estimator that operator / RHS is pivot-affine
             if hasattr(self.HFEngine.A, "is_affine"):
                 self._A_is_affine = self.HFEngine.A.is_affine
             else:
                 self._A_is_affine = 0
             if hasattr(self.HFEngine.b, "is_affine"):
                 self._b_is_affine = self.HFEngine.b.is_affine
             else:
                 self._b_is_affine = 0
             if self._A_is_affine >= 1 / 2 and self._b_is_affine >= 1 / 2:
                 self._affine_lvl += [1 / 2]
         else:
             self._temporaryPivot = 1
             self.trainedModel.data.mu0 = self.checkParameterListPivot(
                                                  self.mu0(self.directionPivot))
             self.trainedModel.data.scaleFactor = self.scaleFactor
             self.trainedModel.data.parameterMap = {
                 "F": [self.HFEngine.parameterMap["F"][self.directionPivot[0]]],
                 "B": [self.HFEngine.parameterMap["B"][self.directionPivot[0]]]}
             self._musUniqueCN = copy(self._m_musUniqueCN)
             self._derIdxs = copy(self._m_derIdxs)
             del self._m_musUniqueCN, self._m_derIdxs
             del self.trainedModel.data.Q._dirPivot
             del self.trainedModel.data.P._dirPivot
             if self._A_is_affine >= 1 / 2 and self._b_is_affine >= 1 / 2:
                 self._affine_lvl.pop()
             del self._A_is_affine, self._b_is_affine
         self.trainedModel.data.npar = self.npar
 
     def errorEstimator(self, mus:Np1D, return_max : bool = False) -> Np1D:
         """Standard residual-based error estimator."""
         setupOK = self.setupApproxLocal()
         if setupOK > 0:
             err = np.empty(len(mus))
             err[:] = np.nan
             if not return_max: return err
             return err, [- setupOK], np.nan
         self._marginalizeTrainedModel(True)
         errRes = super().errorEstimator(mus, return_max)
         self._marginalizeTrainedModel(False)
         return errRes
 
     def _preliminaryTraining(self):
         """Initialize starting snapshots of solution map."""
         RROMPyAssert(self._mode, message = "Cannot start greedy algorithm.")
         self._S = self._setSampleBatch(self.S)
         self.resetSamples()
         self.samplingEngine.scaleFactor = self.scaleFactorDer
         musPivot = self.trainSetGenerator.generatePoints(self.S)
         while len(musPivot) > self.S: musPivot.pop()
         muTestPivot = self.samplerPivot.generatePoints(self.nTestPoints, False)
         idxPop = pruneSamples(self.mapParameterListPivot(muTestPivot),
                               self.mapParameterListPivot(musPivot),
                               1e-10 * self.scaleFactorPivot[0])
         self._mus = emptyParameterList()
         self.mus.reset((self.S, self.npar + len(self.musMargLoc)))
         muTestBase = emptyParameterList()
         muTestBase.reset((len(muTestPivot), self.npar + len(self.musMargLoc)))
         for k in range(self.S):
             muk = np.empty_like(self.mus[0])
             muk[self.directionPivot] = musPivot[k]
             muk[self.directionMarginal] = self.musMargLoc
             self.mus[k] = muk
         for k in range(len(muTestPivot)):
             muk = np.empty_like(muTestBase[0])
             muk[self.directionPivot] = muTestPivot[k]
             muk[self.directionMarginal] = self.musMargLoc
             muTestBase[k] = muk
         muTestBase.pop(idxPop)
         muLast = copy(self.mus[-1])
         self.mus.pop()
         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.muTest = parameterList(muTestBase)
         self.muTest.append(muLast)
         self.M, self.N = ("AUTO",) * 2
 
     def setupApproxLocal(self) -> int:
         """Compute rational interpolant."""
         self._marginalizeMiscellanea(True)
         setupOK = super().setupApproxLocal()
         self._marginalizeMiscellanea(False)
         return setupOK
 
     def setupApprox(self, *args, **kwargs) -> int:
         """Compute rational interpolant."""
         if self.checkComputedApprox(): return -1
         RROMPyAssert(self._mode, message = "Cannot setup approximant.")
         vbMng(self, "INIT", "Setting up {}.". format(self.name()), 5)
         self.computeScaleFactor()
         self._musMarginal = self.samplerMarginal.generatePoints(self.SMarginal)
         while len(self.musMarginal) > self.SMarginal: self.musMarginal.pop()
         S0 = copy(self.S)
         idx, sizes = indicesScatter(len(self.musMarginal), return_sizes = True)
         pMat, Ps, Qs, mus = None, [], [], None
         req, emptyCores = [], np.where(np.logical_not(sizes))[0]
         if len(idx) == 0:
             vbMng(self, "MAIN", "Idling.", 25)
             if self.storeAllSamples: self.storeSamples()
             pL, pT, mT = recv(source = 0, tag = poolRank())
             pMat = np.empty((pL, 0), dtype = pT)
             mus = np.empty((0, self.mu0.shape[1]), dtype = mT)
         else:
             _scaleFactorOldPivot = copy(self.scaleFactor)
             self.scaleFactor = self.scaleFactorPivot
             self._temporaryPivot = 1
             for i in idx:
                 self.musMargLoc = self.musMarginal[i]
                 vbMng(self, "MAIN",
                       "Building marginal model no. {} at {}.".format(i + 1,
                                                            self.musMargLoc), 5)
                 self.samplingEngine.resetHistory()
                 self.trainedModel = None
                 self.verbosity -= 5
                 self.samplingEngine.verbosity -= 10
                 RationalInterpolantGreedy.setupApprox(self, *args, **kwargs)
                 self.verbosity += 5
                 self.samplingEngine.verbosity += 10
                 if self.storeAllSamples: self.storeSamples(i)
                 musi = self.samplingEngine.mus
                 pMati = self.samplingEngine.projectionMatrix
                 if not hasattr(self, "matchState") or not self.matchState:
                     if self.POD == 1 and not (
                         hasattr(self.HFEngine.C, "is_mu_independent")
                     and self.HFEngine.C.is_mu_independent in self._output_lvl):
                         raise RROMPyException(("Cannot apply mu-dependent C "
                                                "to orthonormalized samples."))
                     vbMng(self, "INIT", "Extracting system output from state.",
                           35)
                     if (hasattr(self.HFEngine, "_is_C_quadratic")
                     and self.HFEngine._is_C_quadratic):
                         pMati = self.HFEngine.applyC(pMati, musi)
                     else:
                         pMatiEff = None
                         for j, mu in enumerate(musi):
                             pMij = np.expand_dims(self.HFEngine.applyC(
                                                           pMati[:, j], mu), -1)
                             if pMatiEff is None:
                                 pMatiEff = np.array(pMij)
                             else:
                                 pMatiEff = np.append(pMatiEff, pMij, axis = 1)
                         pMati = pMatiEff
                     vbMng(self, "DEL", "Done extracting system output.", 35)
                 if pMat is None:
                     mus = copy(musi.data)
                     pMat = copy(pMati)
                     if i == 0:
                         for dest in emptyCores:
                             req += [isend((len(pMat), pMat.dtype, mus.dtype),
                                           dest = dest, tag = dest)]
                 else:
                     mus = np.vstack((mus, musi.data))
                     pMat = np.hstack((pMat, pMati))
                 Ps += [copy(self.trainedModel.data.P)]
                 Qs += [copy(self.trainedModel.data.Q)]
                 self._S = S0
             del self._temporaryPivot, self.musMargLoc
             self.scaleFactor = _scaleFactorOldPivot
         for r in req: r.wait()
         pMat, Ps, Qs, mus, nsamples = gatherPivotedApproximant(pMat, Ps, Qs,
                                                                mus, sizes,
                                                                self.polybasis)
         self._mus = self.checkParameterList(mus)
         Psupp = np.append(0, np.cumsum(nsamples))
         self._setupTrainedModel(pMat, forceNew = True)
         self.trainedModel.data.Qs, self.trainedModel.data.Ps = Qs, Ps
         self.trainedModel.data.Psupp = list(Psupp[: -1])
-        self._poleMatching()
+        self._preliminaryMarginalFinalization()
         self._finalizeMarginalization()
         vbMng(self, "DEL", "Done setting up approximant.", 5)
         return 0
 
 class RationalInterpolantGreedyPivotedNoMatch(
                                           RationalInterpolantGreedyPivotedBase,
                                           GenericPivotedApproximantNoMatch):
     """
     ROM pivoted rational interpolant (without pole matching) computation for
         parametric problems.
 
     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': kind of snapshots orthogonalization; allowed values
                 include 0, 1/2, and 1; defaults to 1, i.e. POD;
             - 'scaleFactorDer': scaling factors for derivative computation;
                 defaults to 'AUTO';
             - '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;
             - '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';
             - 'radialDirectionalWeightsMarginal': radial basis weights for
                 marginal interpolant; defaults to 1;
             - 'functionalSolve': strategy for minimization of denominator
                 functional; allowed values include 'NORM', 'DOMINANT', 'NODAL',
                 'BARYCENTRIC_NORM', and 'BARYCENTRIC[_AVERAGE]' (check pdf in
                 main folder for meaning); 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.
         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': kind of snapshots orthogonalization;
             - 'scaleFactorDer': scaling factors for derivative computation;
             - '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;
             - '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.
         verbosity: Verbosity level.
         POD: Kind of snapshots orthogonalization.
         scaleFactorDer: Scaling factors for derivative computation.
         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.
         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.
         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.
         Q: Numpy 1D vector containing complex coefficients of approximant
             denominator.
         P: Numpy 2D vector whose columns are FE dofs of coefficients of
             approximant numerator.
     """
 
-class RationalInterpolantGreedyPivoted(RationalInterpolantGreedyPivotedBase,
-                                       GenericPivotedApproximant):
+class RationalInterpolantGreedyPivotedPoleMatch(
+                                          RationalInterpolantGreedyPivotedBase,
+                                          GenericPivotedApproximantPoleMatch):
     """
     ROM pivoted rational interpolant (with pole matching) computation for
         parametric problems.
 
     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': kind of snapshots orthogonalization; allowed values
                 include 0, 1/2, and 1; defaults to 1, i.e. POD;
             - 'scaleFactorDer': scaling factors for derivative computation;
                 defaults to 'AUTO';
             - 'matchState': whether to match the system state rather than the
                 system output; defaults to False;
             - 'matchingWeight': weight for pole matching optimization; defaults
                 to 1;
             - 'sharedRatio': required ratio of marginal points to share
                 resonance; defaults to 1.;
             - '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;
             - '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';
             - 'radialDirectionalWeightsMarginal': radial basis weights for
                 marginal interpolant; defaults to 1;
             - 'functionalSolve': strategy for minimization of denominator
                 functional; allowed values include 'NORM', 'DOMINANT', 'NODAL',
                 'BARYCENTRIC_NORM', and 'BARYCENTRIC[_AVERAGE]' (check pdf in
                 main folder for meaning); 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.
         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': kind of snapshots orthogonalization;
             - 'scaleFactorDer': scaling factors for derivative computation;
             - 'matchState': whether to match the system state rather than the
                 system output;
             - 'matchingWeight': weight for pole matching optimization;
             - 'sharedRatio': required ratio of marginal points to share
                 resonance;
             - '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;
             - '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.
         verbosity: Verbosity level.
         POD: Kind of snapshots orthogonalization.
         scaleFactorDer: Scaling factors for derivative computation.
         matchState: Whether to match the system state rather than the system
             output.
         matchingWeight: Weight for pole matching optimization.
         sharedRatio: Required ratio of marginal points to share resonance.
         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.
         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.
         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.
         Q: Numpy 1D vector containing complex coefficients of approximant
             denominator.
         P: Numpy 2D vector whose columns are FE dofs of coefficients of
             approximant numerator.
     """
 
     def setupApprox(self, *args, **kwargs) -> int:
         if self.checkComputedApprox(): return -1
         self.purgeparamsMarginal()
         setupOK = super().setupApprox(*args, **kwargs)
         if self.matchState: self._postApplyC()
         return setupOK
diff --git a/rrompy/reduction_methods/pivoted/rational_interpolant_pivoted.py b/rrompy/reduction_methods/pivoted/rational_interpolant_pivoted.py
index b997500..95086fd 100644
--- a/rrompy/reduction_methods/pivoted/rational_interpolant_pivoted.py
+++ b/rrompy/reduction_methods/pivoted/rational_interpolant_pivoted.py
@@ -1,486 +1,487 @@
 # Copyright (C) 2018-2020 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 collections.abc import Iterable
 from copy import deepcopy as copy
 from .generic_pivoted_approximant import (GenericPivotedApproximantBase,
                                           GenericPivotedApproximantNoMatch,
-                                          GenericPivotedApproximant)
+                                          GenericPivotedApproximantPoleMatch)
 from .gather_pivoted_approximant import gatherPivotedApproximant
 from rrompy.reduction_methods.standard.rational_interpolant import (
                                                            RationalInterpolant)
 from rrompy.utilities.base import verbosityManager as vbMng
 from rrompy.utilities.numerical.hash_derivative import nextDerivativeIndices
 from rrompy.utilities.exception_manager import (RROMPyException, RROMPyAssert,
                                                 RROMPyWarning)
 from rrompy.parameter import emptyParameterList
 from rrompy.utilities.parallel import poolRank, indicesScatter, isend, recv
 
-__all__ = ['RationalInterpolantPivotedNoMatch', 'RationalInterpolantPivoted']
+__all__ = ['RationalInterpolantPivotedNoMatch',
+           'RationalInterpolantPivotedPoleMatch']
 
 class RationalInterpolantPivotedBase(GenericPivotedApproximantBase,
                                      RationalInterpolant):
     def __init__(self, *args, **kwargs):
         self._preInit()
         self._addParametersToList(toBeExcluded = ["polydegreetype"])
         super().__init__(*args, **kwargs)
         self._ignoreResidues = self.nparPivot > 1
         self._postInit()
 
     @property
     def scaleFactorDer(self):
         """Value of scaleFactorDer."""
         if self._scaleFactorDer == "NONE": return 1.
         if self._scaleFactorDer == "AUTO": return self.scaleFactorPivot
         return self._scaleFactorDer
     @scaleFactorDer.setter
     def scaleFactorDer(self, scaleFactorDer):
         if isinstance(scaleFactorDer, (str,)):
             scaleFactorDer = scaleFactorDer.upper()
         elif isinstance(scaleFactorDer, Iterable):
             scaleFactorDer = list(scaleFactorDer)
         self._scaleFactorDer = scaleFactorDer
         self._approxParameters["scaleFactorDer"] = self._scaleFactorDer
 
     @property
     def polydegreetype(self):
         """Value of polydegreetype."""
         return "TOTAL"
     @polydegreetype.setter
     def polydegreetype(self, polydegreetype):
         RROMPyWarning(("polydegreetype is used just to simplify inheritance, "
                        "and its value cannot be changed from 'TOTAL'."))
 
     def _setupInterpolationIndices(self):
         """Setup parameters for polyvander."""
         RROMPyAssert(self._mode,
                      message = "Cannot setup interpolation indices.")
         if (self._musUniqueCN is None 
          or len(self._reorder) != len(self.musPivot)):
             try:
                 muPC = self.trainedModel.centerNormalizePivot(self.musPivot)
             except:
                 muPC = self.trainedModel.centerNormalize(self.musPivot)
             self._musUniqueCN, musIdxsTo, musIdxs, musCount = (muPC.unique(
                                     return_index = True, return_inverse = True,
                                     return_counts = True))
             self._musUnique = self.musPivot[musIdxsTo]
             self._derIdxs = [None] * len(self._musUniqueCN)
             self._reorder = np.empty(len(musIdxs), dtype = int)
             filled = 0
             for j, cnt in enumerate(musCount):
                 self._derIdxs[j] = nextDerivativeIndices([], self.nparPivot,
                                                          cnt)
                 jIdx = np.nonzero(musIdxs == j)[0]
                 self._reorder[jIdx] = np.arange(filled, filled + cnt)
                 filled += cnt
 
     def setupApprox(self) -> int:
         """Compute rational interpolant."""
         if self.checkComputedApprox(): return -1
         RROMPyAssert(self._mode, message = "Cannot setup approximant.")
         vbMng(self, "INIT", "Setting up {}.". format(self.name()), 5)
         self.computeScaleFactor()
         self.resetSamples()
         self.samplingEngine.scaleFactor = self.scaleFactorDer
         self.musPivot = self.samplerPivot.generatePoints(self.S)
         while len(self.musPivot) > self.S: self.musPivot.pop()
         self._musMarginal = self.samplerMarginal.generatePoints(self.SMarginal)
         while len(self.musMarginal) > self.SMarginal: self.musMarginal.pop()
         self._mus = emptyParameterList()
         self.mus.reset((self.S * self.SMarginal, self.HFEngine.npar))
         for j, muMarg in enumerate(self.musMarginal):
             for k in range(j * self.S, (j + 1) * self.S):
                 muk = np.empty_like(self.mus[0])
                 muk[self.directionPivot] = self.musPivot[k - j * self.S]
                 muk[self.directionMarginal] = muMarg
                 self.mus[k] = muk
         N0 = copy(self.N)
         self._setupTrainedModel(np.zeros((0, 0)), forceNew = True)
         idx, sizes = indicesScatter(len(self.musMarginal), return_sizes = True)
         pMat, Ps, Qs = None, [], []
         req, emptyCores = [], np.where(np.logical_not(sizes))[0]
         if len(idx) == 0:
             vbMng(self, "MAIN", "Idling.", 30)
             if self.storeAllSamples: self.storeSamples()
             pL, pT = recv(source = 0, tag = poolRank())
             pMat = np.empty((pL, 0), dtype = pT)
         else:
             _scaleFactorOldPivot = copy(self.scaleFactor)
             self.scaleFactor = self.scaleFactorPivot
             self._temporaryPivot = 1
             for i in idx:
                 musi = self.mus[self.S * i : self.S * (i + 1)]
                 vbMng(self, "MAIN",
                       "Building marginal model no. {} at {}.".format(i + 1,
                                                        self.musMarginal[i]), 5)
                 vbMng(self, "INIT", "Starting computation of snapshots.", 10)
                 self.samplingEngine.resetHistory()
                 self.samplingEngine.iterSample(musi)
                 vbMng(self, "DEL", "Done computing snapshots.", 10)
                 self.verbosity -= 5
                 self.samplingEngine.verbosity -= 10
                 self._setupRational(self._setupDenominator())
                 self.verbosity += 5
                 self.samplingEngine.verbosity += 10
                 if self.storeAllSamples: self.storeSamples(i)
                 pMati = self.samplingEngine.projectionMatrix
                 if not hasattr(self, "matchState") or not self.matchState:
                     if self.POD == 1 and not (
                         hasattr(self.HFEngine.C, "is_mu_independent")
                     and self.HFEngine.C.is_mu_independent in self._output_lvl):
                         raise RROMPyException(("Cannot apply mu-dependent C "
                                                "to orthonormalized samples."))
                     vbMng(self, "INIT", "Extracting system output from state.",
                           35)
                     if (hasattr(self.HFEngine, "_is_C_quadratic")
                     and self.HFEngine._is_C_quadratic):
                         pMati = self.HFEngine.applyC(pMati, musi)
                     else:
                         pMatiEff = None
                         for j, mu in enumerate(musi):
                             pMij = np.expand_dims(self.HFEngine.applyC(
                                                           pMati[:, j], mu), -1)
                             if pMatiEff is None:
                                 pMatiEff = np.array(pMij)
                             else:
                                 pMatiEff = np.append(pMatiEff, pMij, axis = 1)
                         pMati = pMatiEff
                     vbMng(self, "DEL", "Done extracting system output.", 35)
                 if pMat is None:
                     pMat = copy(pMati)
                     if i == 0:
                         for dest in emptyCores:
                             req += [isend((len(pMat), pMat.dtype), dest = dest,
                                           tag = dest)]
                 else:
                     pMat = np.hstack((pMat, pMati))
                 Ps += [copy(self.trainedModel.data.P)]
                 Qs += [copy(self.trainedModel.data.Q)]
                 del self.trainedModel.data.Q, self.trainedModel.data.P
                 self.N = N0
             del self._temporaryPivot
         self.scaleFactor = _scaleFactorOldPivot
         for r in req: r.wait()
         pMat, Ps, Qs, _, _ = gatherPivotedApproximant(pMat, Ps, Qs,
                                                       self.mus.data, sizes,
                                                       self.polybasis, False)
         self._setupTrainedModel(pMat)
         self.trainedModel.data.Qs, self.trainedModel.data.Ps = Qs, Ps
         Psupp = np.arange(0, len(self.musMarginal) * self.S, self.S)
         self.trainedModel.data.Psupp = list(Psupp)
-        self._poleMatching()
+        self._preliminaryMarginalFinalization()
         self._finalizeMarginalization()
         vbMng(self, "DEL", "Done setting up approximant.", 5)
         return 0
 
 class RationalInterpolantPivotedNoMatch(RationalInterpolantPivotedBase,
                                         GenericPivotedApproximantNoMatch):
     """
     ROM pivoted rational interpolant (without pole matching) computation for
         parametric problems.
 
     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': kind of snapshots orthogonalization; allowed values
                 include 0, 1/2, and 1; defaults to 1, i.e. POD;
             - 'scaleFactorDer': scaling factors for derivative computation;
                 defaults to 'AUTO';
             - '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;
             - '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;
             - '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',
                 'BARYCENTRIC_NORM', and 'BARYCENTRIC[_AVERAGE]' (check pdf in
                 main folder for meaning); 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.
         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': kind of snapshots orthogonalization;
             - 'scaleFactorDer': scaling factors for derivative computation;
             - 'polybasis': type of polynomial basis for pivot
                 interpolation;
             - 'M': degree of rational interpolant numerator;
             - 'N': degree of rational interpolant denominator;
             - '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.
         verbosity: Verbosity level.
         POD: Kind of snapshots orthogonalization.
         scaleFactorDer: Scaling factors for derivative computation.
         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.
         polybasis: Type of polynomial basis for pivot interpolation.
         M: Numerator degree of approximant.
         N: Denominator degree of approximant.
         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.
         Q: Numpy 1D vector containing complex coefficients of approximant
             denominator.
         P: Numpy 2D vector whose columns are FE dofs of coefficients of
             approximant numerator.
     """
     
-class RationalInterpolantPivoted(RationalInterpolantPivotedBase,
-                                 GenericPivotedApproximant):
+class RationalInterpolantPivotedPoleMatch(RationalInterpolantPivotedBase,
+                                          GenericPivotedApproximantPoleMatch):
     """
     ROM pivoted rational interpolant (with pole matching) computation for
         parametric problems.
 
     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': kind of snapshots orthogonalization; allowed values
                 include 0, 1/2, and 1; defaults to 1, i.e. POD;
             - 'scaleFactorDer': scaling factors for derivative computation;
                 defaults to 'AUTO';
             - 'matchState': whether to match the system state rather than the
                 system output; defaults to False;
             - 'matchingWeight': weight for pole matching optimization; defaults
                 to 1;
             - 'sharedRatio': required ratio of marginal points to share
                 resonance; defaults to 1.;
             - '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;
             - '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;
             - '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',
                 'BARYCENTRIC_NORM', and 'BARYCENTRIC[_AVERAGE]' (check pdf in
                 main folder for meaning); 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.
         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': kind of snapshots orthogonalization;
             - 'scaleFactorDer': scaling factors for derivative computation;
             - 'matchState': whether to match the system state rather than the
                 system output;
             - 'matchingWeight': weight for pole matching optimization;
             - 'sharedRatio': required ratio of marginal points to share
                 resonance;
             - '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;
             - '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.
         verbosity: Verbosity level.
         POD: Kind of snapshots orthogonalization.
         scaleFactorDer: Scaling factors for derivative computation.
         matchState: Whether to match the system state rather than the system
             output.
         matchingWeight: Weight for pole matching optimization.
         sharedRatio: Required ratio of marginal points to share resonance.
         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.
         polybasis: Type of polynomial basis for pivot interpolation.
         polybasisMarginal: Type of polynomial basis for marginal interpolation.
         paramsMarginal: Dictionary of parameters for marginal interpolation.
         M: Numerator degree of approximant.
         N: Denominator degree of approximant.
         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.
         Q: Numpy 1D vector containing complex coefficients of approximant
             denominator.
         P: Numpy 2D vector whose columns are FE dofs of coefficients of
             approximant numerator.
     """
 
     def setupApprox(self, *args, **kwargs) -> int:
         if self.checkComputedApprox(): return -1
         self.purgeparamsMarginal()
         setupOK = super().setupApprox(*args, **kwargs)
         if self.matchState: self._postApplyC()
         return setupOK
diff --git a/rrompy/reduction_methods/pivoted/trained_model/trained_model_pivoted_rational_nomatch.py b/rrompy/reduction_methods/pivoted/trained_model/trained_model_pivoted_rational_nomatch.py
index a2a5cc7..4eacc69 100644
--- a/rrompy/reduction_methods/pivoted/trained_model/trained_model_pivoted_rational_nomatch.py
+++ b/rrompy/reduction_methods/pivoted/trained_model/trained_model_pivoted_rational_nomatch.py
@@ -1,367 +1,284 @@
 # Copyright (C) 2018-2020 by the RROMPy authors
 #
 # This file is part of RROMPy.
 #
 # RROMPy is free software: you can redistribute it and/or modify
 # it under the terms of the GNU Lesser General Public License as published by
 # the Free Software Foundation, either version 3 of the License, or
 # (at your option) any later version.
 #
 # RROMPy is distributed in the hope that it will be useful,
 # but WITHOUT ANY WARRANTY; without even the implied warranty of
 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
 # GNU Lesser General Public License for more details.
 #
 # You should have received a copy of the GNU Lesser General Public License
 # along with RROMPy. If not, see <http://www.gnu.org/licenses/>.
 #
 
 import numpy as np
-from scipy.special import factorial as fact
 from collections.abc import Iterable
-from copy import deepcopy as copy
-from itertools import combinations
 from rrompy.reduction_methods.standard.trained_model.trained_model_rational \
                                                     import TrainedModelRational
 from rrompy.utilities.base.types import (Np1D, Np2D, List, ListAny, paramVal,
                                          paramList, sampList)
 from rrompy.utilities.base import verbosityManager as vbMng, freepar as fp
 from rrompy.utilities.numerical import dot
 from rrompy.utilities.numerical.compress_matrix import compressMatrix
-from rrompy.utilities.poly_fitting.heaviside import (rational2heaviside,
-                                                   HeavisideInterpolator as HI)
+from rrompy.utilities.poly_fitting.heaviside import rational2heaviside
 from rrompy.utilities.poly_fitting.nearest_neighbor import (
                                             NearestNeighborInterpolator as NNI)
 from rrompy.utilities.exception_manager import (RROMPyException, RROMPyAssert,
                                                 RROMPyWarning)
 from rrompy.parameter import checkParameterList
-from rrompy.sampling import sampleList, emptySampleList
+from rrompy.sampling import emptySampleList
 
 __all__ = ['TrainedModelPivotedRationalNoMatch']
 
 class TrainedModelPivotedRationalNoMatch(TrainedModelRational):
     """
     ROM approximant evaluation for pivoted approximants based on interpolation
         of rational approximants (without pole matching).
     
     Attributes:
         Data: dictionary with all that can be pickled.
     """
 
+    @property
+    def supportEnd(self):
+        lookAtExcl = hasattr(self, "_PsuppExcl") and len(self._PsuppExcl) > 0
+        idxIncl = np.argmax(self.data.Psupp)
+        if lookAtExcl: idxExcl = np.argmax(self._PsuppExcl)
+        if (not lookAtExcl
+         or self.data.Psupp[idxIncl] >= self._PsuppExcl[idxExcl]):
+            return self.data.Psupp[idxIncl] + self.data.Ps[idxIncl].shape[0]
+        return self._PsuppExcl[idxExcl] + self._PsExcl[idxExcl].shape[0]
+
     def checkParameterListPivot(self, mu:paramList,
                                 check_if_single : bool = False) -> paramList:
         return checkParameterList(mu, self.data.nparPivot, check_if_single)
 
     def checkParameterListMarginal(self, mu:paramList,
                                   check_if_single : bool = False) -> paramList:
         return checkParameterList(mu, self.data.nparMarginal, check_if_single)
 
-    def compress(self, collapse : bool = False, tol : float = 0., *args,
-                 **kwargs):
+    def compress(self, collapse : bool = False, tol : float = 0.,
+                 returnRMat : bool = False, **compressMatrixkwargs):
         if not collapse and tol <= 0.: return
         if hasattr(self.data, "_is_C_quadratic") and self.data._is_C_quadratic:
             raise RROMPyException(("Cannot compress model with quadratic "
                                    "output."))
         RMat = self.data.projMat
         if not collapse:
             if hasattr(self.data, "_compressTol"):
                 RROMPyWarning(("Recompressing already compressed model is "
                                "ineffective. Aborting."))
                 return
-            self.data.projMat, RMat, _ = compressMatrix(RMat, tol, *args,
-                                                        **kwargs)
-        for obj, suppj in zip(self.data.HIs, self.data.Psupp):
-            obj.postmultiplyTensorize(RMat.T[suppj : suppj + obj.shape[0]])
-        if hasattr(self, "_HIsExcl"):
-            for obj, suppj in zip(self._HIsExcl, self.data.Psupp):
-                obj.postmultiplyTensorize(RMat.T[suppj : suppj + obj.shape[0]])
+            self.data.projMat, RMat, _ = compressMatrix(RMat, tol,
+                                                        **compressMatrixkwargs)
         if hasattr(self.data, "Ps"):
             for obj, suppj in zip(self.data.Ps, self.data.Psupp):
                 obj.postmultiplyTensorize(RMat.T[suppj : suppj + obj.shape[0]])
         if hasattr(self, "_PsExcl"):
             for obj, suppj in zip(self._PsExcl, self._PsuppExcl):
                 obj.postmultiplyTensorize(RMat.T[suppj : suppj + obj.shape[0]])
-        if hasattr(self.data, "coeffsEff"):
-            for j in range(len(self.data.coeffsEff)):
-                self.data.coeffsEff[j] = dot(self.data.coeffsEff[j], RMat.T)
-        if hasattr(self, "_HIsExcl") or hasattr(self, "_PsExcl"):
             self._PsuppExcl = [0] * len(self._PsuppExcl)
         self.data.Psupp = [0] * len(self.data.Psupp)
         super(TrainedModelRational, self).compress(collapse, tol)
+        if returnRMat: return RMat
 
     def centerNormalizePivot(self, mu : paramList = [],
                              mu0 : paramVal = None) -> paramList:
         """
         Compute normalized parameter to be plugged into approximant.
 
         Args:
             mu: Parameter(s) 1.
             mu0: Parameter(s) 2. If None, set to self.data.mu0Pivot.
 
         Returns:
             Normalized parameter.
         """
         mu = self.checkParameterListPivot(mu)
         if mu0 is None:
             mu0 = self.checkParameterListPivot(
                                     self.data.mu0(0, self.data.directionPivot))
         return (self.mapParameterList(mu, idx = self.data.directionPivot)
               - self.mapParameterList(mu0, idx = self.data.directionPivot)
                ) / [self.data.scaleFactor[x] for x in self.data.directionPivot]
 
     def setupMarginalInterp(self, interpPars:ListAny):
         self.data.marginalInterp = NNI()
         self.data.marginalInterp.setupByInterpolation(self.data.musMarginal,
                                          np.arange(len(self.data.musMarginal)),
                                          1, *interpPars)
 
-    def updateEffectiveSamples(self, exclude:List[int], *args, **kwargs):
+    def updateEffectiveSamples(self, exclude:List[int]):
         if hasattr(self, "_idxExcl"):
             for j, excl in enumerate(self._idxExcl):
                 self.data.musMarginal.insert(self._musMExcl[j], excl)
-                self.data.HIs.insert(excl, self._HIsExcl[j])
                 self.data.Ps.insert(excl, self._PsExcl[j])
                 self.data.Qs.insert(excl, self._QsExcl[j])
                 self.data.Psupp.insert(excl, self._PsuppExcl[j])
         self._idxExcl, self._musMExcl = list(np.sort(exclude)), []
-        self._HIsExcl, self._PsExcl, self._QsExcl = [], [], []
-        self._PsuppExcl = []
+        self._PsExcl, self._QsExcl, self._PsuppExcl = [], [], []
         for excl in self._idxExcl[::-1]:
             self._musMExcl = [self.data.musMarginal[excl]] + self._musMExcl
             self.data.musMarginal.pop(excl)
-            self._HIsExcl = [self.data.HIs.pop(excl)] + self._HIsExcl
             self._PsExcl = [self.data.Ps.pop(excl)] + self._PsExcl
             self._QsExcl = [self.data.Qs.pop(excl)] + self._QsExcl
             self._PsuppExcl = [self.data.Psupp.pop(excl)] + self._PsuppExcl
-        poles = [hi.poles for hi in self.data.HIs]
-        coeffs = [hi.coeffs for hi in self.data.HIs]
-        self.initializeFromLists(poles, coeffs, self.data.Psupp,
-                                 self.data.HIs[0].polybasis, *args, **kwargs)
-
-    def initializeFromLists(self, poles:ListAny, coeffs:ListAny, supps:ListAny,
-                            basis:str, *args, **kwargs):
-        """Initialize Heaviside representation."""
-        self.data.HIs = []
-        for pls, cfs in zip(poles, coeffs):
-            hsi = HI()
-            hsi.poles = pls
-            if len(cfs) == len(pls):
-                cfs = np.pad(cfs, ((0, 1), (0, 0)), "constant")
-            hsi.coeffs = cfs
-            hsi.npar = 1
-            hsi.polybasis = basis
-            self.data.HIs += [hsi]
-        lookAtExcl = hasattr(self, "_PsuppExcl") and len(self._PsuppExcl) > 0
-        idxIncl = np.argmax(self.data.Psupp)
-        if lookAtExcl: idxExcl = np.argmax(self._PsuppExcl)
-        if (not lookAtExcl
-         or self.data.Psupp[idxIncl] >= self._PsuppExcl[idxExcl]):
-            sEnd = self.data.Psupp[idxIncl] + self.data.HIs[idxIncl].shape[0]
-        else:
-            sEnd = self._PsuppExcl[idxExcl] + self._HIsExcl[idxExcl].shape[0]
-        self.data.supportEnd = sEnd
-
-    def initializeFromRational(self, *args, **kwargs):
-        """Initialize Heaviside representation."""
-        RROMPyAssert(self.data.nparPivot, 1, "Number of pivot parameters")
-        poles, coeffs = [], []
-        for Q, P in zip(self.data.Qs, self.data.Ps):
-            cfs, pls, basis = rational2heaviside(P, Q)
-            poles += [pls]
-            coeffs += [cfs]
-        self.initializeFromLists(poles, coeffs, self.data.Psupp, basis, *args,
-                                 **kwargs)
-
-    def interpolateMarginalInterpolator(self, mu : paramList = []) -> ListAny:
-        """Obtain interpolated approximant interpolator."""
-        mu = self.checkParameterListMarginal(mu)
-        vbMng(self, "INIT", "Finding nearest neighbor to mu = {}.".format(mu),
-              95)
-        his = []
-        intM = np.array(self.data.marginalInterp(mu), dtype = int)
-        for j in range(len(mu)):
-            i = intM[j]
-            his += [HI()]
-            his[-1].poles = copy(self.data.HIs[i].poles)
-            his[-1].coeffs = copy(self.data.HIs[i].coeffs)
-            his[-1].npar = 1
-            his[-1].polybasis = self.data.HIs[0].polybasis
-            if not self.data._collapsed:
-                his[-1].pad(self.data.Psupp[i], self.data.supportEnd
-                                       - self.data.Psupp[i] - his[-1].shape[0])
-        vbMng(self, "DEL", "Done finding nearest neighbor.", 95)
-        return his
-
-    def interpolateMarginalPoles(self, mu : paramList = []) -> ListAny:
-        """Obtain interpolated approximant poles."""
-        interps = self.interpolateMarginalInterpolator(mu)
-        return [interp.poles for interp in interps]
 
-    def interpolateMarginalCoeffs(self, mu : paramList = []) -> ListAny:
-        """Obtain interpolated approximant poles."""
-        interps = self.interpolateMarginalInterpolator(mu)
-        return [interp.coeffs for interp in interps]
-
-    def getApproxReduced(self, mu : paramList = []) -> sampList:
-        """
-        Evaluate reduced representation of approximant at arbitrary parameter.
-
-        Args:
-            mu: Target parameter.
-        """
-        RROMPyAssert(self.data.nparPivot, 1, "Number of pivot parameters")
-        mu = self.checkParameterList(mu)
-        if (not hasattr(self, "lastSolvedApproxReduced")
-         or self.lastSolvedApproxReduced != mu):
-            vbMng(self, "INIT",
-                  "Evaluating approximant at mu = {}.".format(mu), 12)
-            muP = self.centerNormalizePivot(mu(self.data.directionPivot))
-            muM = mu(self.data.directionMarginal)
-            his = self.interpolateMarginalInterpolator(muM)
-            for i, (mP, hi) in enumerate(zip(muP, his)):
-                uAppR = hi(mP)[:, 0]
-                if i == 0:
-                    uApproxR = np.empty((len(uAppR), len(mu)),
-                                        dtype = uAppR.dtype)
-                uApproxR[:, i] = uAppR
-            self.uApproxReduced = sampleList(uApproxR)
-            vbMng(self, "DEL", "Done evaluating approximant.", 12)
-            self.lastSolvedApproxReduced = mu
-        return self.uApproxReduced
-    
     def _setupQuadMapping(self, N2 : List[int] = None):
         if not (hasattr(self.data, "_is_C_quadratic")
             and self.data._is_C_quadratic):
             RROMPyWarning(("Quadratic mapping is useless if output is not "
                            "quadratic. Aborting."))
             return
         if N2 is None:
-            N2 = np.array([hi.shape[0] for hi in self.data.HIs], dtype = int)
+            N2 = np.array([P.shape[0] for P in self.data.Ps], dtype = int)
             if self.data._is_C_quadratic == 1:
                 N2 = N2 ** 2
             else: # if self.data._is_C_quadratic == 2:
                 N2 = N2 * (N2 + 1) // 2
         if (hasattr(self, "quad_mapping")
         and self.quad_mapping.shape[1] == np.sum(N2)): return
         quad_mapping = np.empty((2, 0), dtype = int)
         for Nloc, suppj in zip(N2, self.data.Psupp):
             self.quad_mapping = np.empty((0, 0))
             super()._setupQuadMapping(Nloc)
             quad_mapping = np.append(quad_mapping, self.quad_mapping + suppj,
                                      axis = 1)
         self.quad_mapping = quad_mapping
 
     def getPVal(self, mu : paramList = []) -> sampList:
         """
         Evaluate rational numerator at arbitrary parameter.
 
         Args:
             mu: Target parameter.
         """
         RROMPyAssert(self.data.nparPivot, 1, "Number of pivot parameters")
         mu = self.checkParameterList(mu)
-        p = emptySampleList()
         muP = self.centerNormalizePivot(mu(self.data.directionPivot))
-        muM = mu(self.data.directionMarginal)
-        his = self.interpolateMarginalInterpolator(muM)
-        for i, (mP, hi) in enumerate(zip(muP, his)):
-            Pval = hi(mP) * np.prod(mP[0] - hi.poles)
-            if i == 0: p.reset((len(Pval), len(mu)), dtype = Pval.dtype)
-            p[i] = Pval
+        muM = self.checkParameterListMarginal(mu(self.data.directionMarginal))
+        intM = np.array(self.data.marginalInterp(muM), dtype = int)
+        p = emptySampleList()
+        vbMng(self, "INIT", "Evaluating numerator at mu = {}.".format(mu), 17)
+        for i, iM in enumerate(intM):
+            Pval, supp = self.data.Ps[iM](muP[i]), self.data.Psupp[iM]
+            if i == 0:
+                p.reset((self.supportEnd, len(mu)), dtype = Pval.dtype)
+                p.data[:] = 0.
+            p.data[supp : supp + len(Pval), i] = Pval[:, 0]
+        vbMng(self, "DEL", "Done evaluating numerator.", 17)
         return p
 
     def getQVal(self, mu:Np1D, der : List[int] = None,
                 scl : Np1D = None) -> Np1D:
         """
         Evaluate rational denominator at arbitrary parameter.
 
         Args:
             mu: Target parameter.
             der(optional): Derivatives to take before evaluation.
         """
         RROMPyAssert(self.data.nparPivot, 1, "Number of pivot parameters")
         mu = self.checkParameterList(mu)
         muP = self.centerNormalizePivot(mu(self.data.directionPivot))
-        muM = mu(self.data.directionMarginal)
+        muM = self.checkParameterListMarginal(mu(self.data.directionMarginal))
         if der is None:
             derP, derM = 0, [0]
         else:
             derP = der[self.data.directionPivot[0]]
             derM = [der[x] for x in self.data.directionMarginal]
         if np.any(np.array(derM) != 0):
             raise RROMPyException(("Derivatives of Q with respect to marginal "
                                    "parameters not allowed."))
         sclP = 1 if scl is None else scl[self.data.directionPivot[0]]
-        derVal = np.zeros(len(mu), dtype = np.complex)
-        pls = self.interpolateMarginalPoles(muM)
-        for i, (mP, pl) in enumerate(zip(muP, pls)):
-            N = len(pl)
-            if derP == N: derVal[i] = 1.
-            elif derP >= 0 and derP < N:
-                plDist = muP[0] - pl
-                for terms in combinations(np.arange(N), N - derP):
-                    derVal[i] += np.prod(plDist[list(terms)], axis = 1)
-        return sclP ** derP * fact(derP) * derVal
+        intM = np.array(self.data.marginalInterp(muM), dtype = int)
+        vbMng(self, "INIT", "Evaluating denominator at mu = {}.".format(mu),
+              17)
+        q = np.array([self.data.Qs[iM](mP, derP, sclP)[0]
+                                                 for iM, mP in zip(intM, muP)])
+        vbMng(self, "DEL", "Done evaluating denominator.", 17)
+        return q
 
     def getPoles(self, *args, **kwargs) -> Np1D:
         """
         Obtain approximant poles.
 
         Returns:
             Numpy complex vector of poles.
         """
         RROMPyAssert(self.data.nparPivot, 1, "Number of pivot parameters")
         if len(args) + len(kwargs) > 1:
             raise RROMPyException(("Wrong number of parameters passed. "
                                    "Only 1 available."))
         elif len(args) + len(kwargs) == 1:
             if len(args) == 1:
                 mVals = args[0]
             else:
                 mVals = kwargs["marginalVals"]
             if not isinstance(mVals, Iterable): mVals = [mVals]
             mVals = list(mVals)
         else:
             mVals = [fp]
         rDim = mVals.index(fp)
         if rDim < len(mVals) - 1 and fp in mVals[rDim + 1 :]:
             raise RROMPyException(("Exactly 1 'freepar' entry in "
                                    "marginalVals must be provided."))
         if rDim != self.data.directionPivot[0]:
             raise RROMPyException(("'freepar' entry in marginalVals must "
                                    "coincide with pivot direction."))
         mVals[rDim] = self.data.mu0(rDim)[0]
-        mMarg = [mVals[j] for j in range(len(mVals)) if j != rDim]
-        roots = (self.data.scaleFactor[rDim]
-               * self.interpolateMarginalPoles(mMarg)[0])
+        muM = self.checkParameterListMarginal([mVals[j]
+                                      for j in range(len(mVals)) if j != rDim])
+        iM = int(self.data.marginalInterp(muM))
+        roots = self.data.scaleFactor[rDim] * self.data.Qs[iM].roots()
         return self.mapParameterList(self.mapParameterList(self.data.mu0(rDim),
                                                            idx = [rDim])(0, 0)
                                    + roots, "B", [rDim])(0)
 
     def getResidues(self, *args, **kwargs) -> Np2D:
         """
         Obtain approximant residues.
 
         Returns:
             Numpy matrix with residues as columns.
         """
-        pls = self.getPoles(*args, **kwargs)
-        if len(args) == 1:
-            mVals = args[0]
-        elif len(args) == 0:
-            mVals = [None]
+        RROMPyAssert(self.data.nparPivot, 1, "Number of pivot parameters")
+        if len(args) + len(kwargs) > 1:
+            raise RROMPyException(("Wrong number of parameters passed. "
+                                   "Only 1 available."))
+        elif len(args) + len(kwargs) == 1:
+            if len(args) == 1:
+                mVals = args[0]
+            else:
+                mVals = kwargs["marginalVals"]
+            if not isinstance(mVals, Iterable): mVals = [mVals]
+            mVals = list(mVals)
         else:
-            mVals = kwargs["marginalVals"]
-        if not isinstance(mVals, Iterable): mVals = [mVals]
-        mVals = list(mVals)
+            mVals = [fp]
         rDim = mVals.index(fp)
-        mMarg = [mVals[j] for j in range(len(mVals)) if j != rDim]
-        res = self.interpolateMarginalCoeffs(mMarg)[0][: len(pls), :].T
+        if rDim < len(mVals) - 1 and fp in mVals[rDim + 1 :]:
+            raise RROMPyException(("Exactly 1 'freepar' entry in "
+                                   "marginalVals must be provided."))
+        if rDim != self.data.directionPivot[0]:
+            raise RROMPyException(("'freepar' entry in marginalVals must "
+                                   "coincide with pivot direction."))
+        mVals[rDim] = self.data.mu0(rDim)[0]
+        muM = self.checkParameterListMarginal([mVals[j]
+                                      for j in range(len(mVals)) if j != rDim])
+        iM = int(self.data.marginalInterp(muM))
+        res, pls, basis = rational2heaviside(self.data.Ps[iM],
+                                             self.data.Qs[iM])
+        res = res[: len(pls), :].T
         if hasattr(self.data, "_is_C_quadratic") and self.data._is_C_quadratic:
             self._setupQuadMapping()
             res = res[self.quad_mapping[0]] * res[self.quad_mapping[1]].conj()
         if not self.data._collapsed: res = dot(self.data.projMat, res).T
         if (hasattr(self.data, "_is_C_quadratic")
         and self.data._is_C_quadratic == 2):
             res = np.real(res)
         return pls, res
diff --git a/rrompy/reduction_methods/pivoted/trained_model/trained_model_pivoted_rational.py b/rrompy/reduction_methods/pivoted/trained_model/trained_model_pivoted_rational_polematch.py
similarity index 61%
rename from rrompy/reduction_methods/pivoted/trained_model/trained_model_pivoted_rational.py
rename to rrompy/reduction_methods/pivoted/trained_model/trained_model_pivoted_rational_polematch.py
index d0191f4..21ffc99 100644
--- a/rrompy/reduction_methods/pivoted/trained_model/trained_model_pivoted_rational.py
+++ b/rrompy/reduction_methods/pivoted/trained_model/trained_model_pivoted_rational_polematch.py
@@ -1,308 +1,503 @@
 # Copyright (C) 2018-2020 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 warnings
 import numpy as np
+from scipy.special import factorial as fact
 from scipy.sparse import csr_matrix, hstack, SparseEfficiencyWarning
+from collections.abc import Iterable
 from copy import deepcopy as copy
+from itertools import combinations
 from rrompy.reduction_methods.standard.trained_model.trained_model_rational \
                                                     import TrainedModelRational
 from .trained_model_pivoted_rational_nomatch import (
                                             TrainedModelPivotedRationalNoMatch)
-from rrompy.utilities.base.types import (Np2D, ListAny, paramVal, paramList,
-                                         HFEng)
-from rrompy.utilities.base import verbosityManager as vbMng
+from rrompy.utilities.base.types import (Np1D, Np2D, List, ListAny, paramVal,
+                                         paramList, sampList, HFEng)
+from rrompy.utilities.base import verbosityManager as vbMng, freepar as fp
+from rrompy.utilities.numerical import dot
 from rrompy.utilities.numerical.point_matching import rationalFunctionMatching
 from rrompy.utilities.numerical.degree import reduceDegreeN
 from rrompy.utilities.poly_fitting.polynomial import (polybases as ppb,
                                                   PolynomialInterpolator as PI)
 from rrompy.utilities.poly_fitting.radial_basis import (polybases as rbpb,
                                                 RadialBasisInterpolator as RBI)
-from rrompy.utilities.poly_fitting.heaviside import (heavisideUniformShape,
+from rrompy.utilities.poly_fitting.heaviside import (rational2heaviside,
+                                                   heavisideUniformShape,
                                                    HeavisideInterpolator as HI)
 from rrompy.utilities.poly_fitting.nearest_neighbor import (
                                             NearestNeighborInterpolator as NNI)
 from rrompy.utilities.poly_fitting.piecewise_linear import (sparsekinds,
                                             PiecewiseLinearInterpolator as PLI)
-from rrompy.utilities.exception_manager import RROMPyException
+from rrompy.utilities.exception_manager import RROMPyException, RROMPyAssert
+from rrompy.sampling import sampleList, emptySampleList
 
-__all__ = ['TrainedModelPivotedRational']
+__all__ = ['TrainedModelPivotedRationalPoleMatch']
 
-class TrainedModelPivotedRational(TrainedModelPivotedRationalNoMatch):
+class TrainedModelPivotedRationalPoleMatch(TrainedModelPivotedRationalNoMatch):
     """
     ROM approximant evaluation for pivoted approximants based on interpolation
         of rational approximants (with pole matching).
     
     Attributes:
         Data: dictionary with all that can be pickled.
     """
 
+    def compress(self, collapse : bool = False, tol : float = 0.,
+                 returnRMat : bool = False, **compressMatrixkwargs):
+        Psupp = copy(self.data.Psupp)
+        RMat = super().compress(collapse, tol, True, **compressMatrixkwargs)
+        if RMat is None: return
+        for obj, suppj in zip(self.data.HIs, Psupp):
+            obj.postmultiplyTensorize(RMat.T[suppj : suppj + obj.shape[0]])
+        if hasattr(self, "_HIsExcl"):
+            for obj, suppj in zip(self._HIsExcl, Psupp):
+                obj.postmultiplyTensorize(RMat.T[suppj : suppj + obj.shape[0]])
+            if not hasattr(self, "_PsExcl"):
+                self._PsuppExcl = [0] * len(self._PsuppExcl)
+        if returnRMat: return RMat
+    
     def centerNormalizeMarginal(self, mu : paramList = [],
                                 mu0 : paramVal = None) -> paramList:
         """
         Compute normalized parameter to be plugged into approximant.
 
         Args:
             mu: Parameter(s) 1.
             mu0: Parameter(s) 2. If None, set to self.data.mu0Marginal.
 
         Returns:
             Normalized parameter.
         """
         mu = self.checkParameterListMarginal(mu)
         if mu0 is None:
             mu0 = self.checkParameterListMarginal(
                                  self.data.mu0(0, self.data.directionMarginal))
         return (self.mapParameterList(mu, idx = self.data.directionMarginal)
               - self.mapParameterList(mu0, idx = self.data.directionMarginal)
                ) / [self.data.scaleFactor[x]
                                           for x in self.data.directionMarginal]
 
     def setupMarginalInterp(self, approx, interpPars:ListAny, extraPar = None):
         vbMng(self, "INIT", "Starting computation of marginal interpolator.",
               12)
         musMCN = self.centerNormalizeMarginal(self.data.musMarginal)
         nM, pbM = len(musMCN), approx.polybasisMarginal
         if pbM in ppb + rbpb:
             if extraPar: approx._setMMarginalAuto()
             _MMarginalEff = approx.paramsMarginal["MMarginal"]
         if pbM in ppb:
             p = PI()
         elif pbM in rbpb:
             p = RBI()
         else: # if pbM in sparsekinds + ["NEARESTNEIGHBOR"]:
             if pbM == "NEARESTNEIGHBOR":
                 p = NNI()
             else: # if pbM in sparsekinds:
                 pllims = [[-1.] * self.data.nparMarginal,
                           [1.] * self.data.nparMarginal]
                 p = PLI()
         for ipts, pts in enumerate(self.data.suppEffPts):
             if len(pts) == 0:
                 raise RROMPyException("Empty list of support points.")
             musMCNEff, valsEff = musMCN[pts], np.eye(len(pts))
             if pbM in ppb + rbpb:
                 if extraPar:
                     if ipts > 0:
                         verb = approx.verbosity
                         approx.verbosity = 0
                         _musM = approx.musMarginal
                         approx.musMarginal = musMCNEff
                         approx._setMMarginalAuto()
                         approx.musMarginal = _musM
                         approx.verbosity = verb
                 else:
                     approx.paramsMarginal["MMarginal"] = reduceDegreeN(
                          _MMarginalEff, len(musMCNEff), self.data.nparMarginal,
                          approx.paramsMarginal["polydegreetypeMarginal"])
                 MMEff = approx.paramsMarginal["MMarginal"]
                 while MMEff >= 0:
                     wellCond, msg = p.setupByInterpolation(musMCNEff, valsEff,
                                                            MMEff, *interpPars)
                     vbMng(self, "MAIN", msg, 30)
                     if wellCond: break
                     vbMng(self, "MAIN",
                           ("Polyfit is poorly conditioned. Reducing "
                            "MMarginal by 1."), 35)
                     MMEff -= 1
                 if MMEff < 0:
                     raise RROMPyException(("Instability in computation of "
                                            "interpolant. Aborting."))
                 if (pbM in rbpb and len(interpPars) > 4
                 and "optimizeScalingBounds" in interpPars[4].keys()):
                     interpPars[4]["optimizeScalingBounds"] = [-1., -1.]
             elif pbM == "NEARESTNEIGHBOR":
                 if ipts > 0: interpPars[0] = 1
                 p.setupByInterpolation(musMCNEff, valsEff, *interpPars)
             elif ipts == 0: # and pbM in sparsekinds:
                 p.setupByInterpolation(musMCNEff, valsEff, pllims,
                                        extraPar[pts], *interpPars)
             if ipts == 0:
                 self.data.marginalInterp = copy(p)
                 self.data.coeffsEff, self.data.polesEff = [], []
                 for hi, sup in zip(self.data.HIs, self.data.Psupp):
                     cEff = hi.coeffs
                     if (self.data._collapsed
-                     or self.data.supportEnd == cEff.shape[1]):
+                     or self.supportEnd == cEff.shape[1]):
                         cEff = copy(cEff)
                     else:
-                        supC = self.data.supportEnd - sup - cEff.shape[1]
+                        supC = self.supportEnd - sup - cEff.shape[1]
                         cEff = hstack((csr_matrix((len(cEff), sup)),
                                        csr_matrix(cEff),
                                        csr_matrix((len(cEff), supC))), "csr")
                     self.data.coeffsEff += [cEff]
                     self.data.polesEff += [copy(hi.poles)]
             else:
                 ptsBad = [i for i in range(nM) if i not in pts]
                 idxBad = np.where(self.data.suppEffIdx == ipts)[0]
                 warnings.simplefilter('ignore', SparseEfficiencyWarning)
                 if pbM in sparsekinds:
                     for ij, j in enumerate(ptsBad):
                         nearest = pts[np.argmin(np.sum(np.abs(musMCNEff.data
                                             - np.tile(musMCN[j], [len(pts), 1])
                                                 ), axis = 1).flatten())]
                         self.data.coeffsEff[j][idxBad] = copy(
                                           self.data.coeffsEff[nearest][idxBad])
                         self.data.polesEff[j][idxBad] = copy(
                                            self.data.polesEff[nearest][idxBad])
                 else:
                     if (self.data._collapsed
-                     or self.data.supportEnd == cEff.shape[1]):
+                     or self.supportEnd == cEff.shape[1]):
                         cfBase = np.zeros((len(idxBad), cEff.shape[1]),
                                           dtype = cEff.dtype)
                     else:
-                        cfBase = csr_matrix((len(idxBad),
-                                             self.data.supportEnd), 
+                        cfBase = csr_matrix((len(idxBad), self.supportEnd), 
                                             dtype = cEff.dtype)
                     valMuMBad = p(musMCN[ptsBad])
                     for ijb, jb in enumerate(ptsBad):
                         self.data.coeffsEff[jb][idxBad] = copy(cfBase)
                         self.data.polesEff[jb][idxBad] = 0.
                         for ij, j in enumerate(pts):
                             val = valMuMBad[ij][ijb]
                             if not np.isclose(val, 0.):
                                 self.data.coeffsEff[jb][idxBad] += (val
                                               * self.data.coeffsEff[j][idxBad])
                                 self.data.polesEff[jb][idxBad] += (val
                                                * self.data.polesEff[j][idxBad])
                 warnings.filters.pop(0)
         if pbM in ppb + rbpb:
             approx.paramsMarginal["MMarginal"] = _MMarginalEff
         vbMng(self, "DEL", "Done computing marginal interpolator.", 12)
 
+    def updateEffectiveSamples(self, exclude:List[int], *args, **kwargs):
+        if hasattr(self, "_idxExcl"):
+            for j, excl in enumerate(self._idxExcl):
+                self.data.HIs.insert(excl, self._HIsExcl[j])
+        super().updateEffectiveSamples(exclude)
+        self._HIsExcl = []
+        for excl in self._idxExcl[::-1]:
+            self._HIsExcl = [self.data.HIs.pop(excl)] + self._HIsExcl
+        poles = [hi.poles for hi in self.data.HIs]
+        coeffs = [hi.coeffs for hi in self.data.HIs]
+        self.initializeFromLists(poles, coeffs, self.data.Psupp,
+                                 self.data.HIs[0].polybasis, *args, **kwargs)
+
+    def initializeFromRational(self, matchingWeight:float, HFEngine:HFEng,
+                               is_state:bool):
+        """Initialize Heaviside representation."""
+        RROMPyAssert(self.data.nparPivot, 1, "Number of pivot parameters")
+        poles, coeffs = [], []
+        for Q, P in zip(self.data.Qs, self.data.Ps):
+            cfs, pls, basis = rational2heaviside(P, Q)
+            poles += [pls]
+            coeffs += [cfs]
+        self.initializeFromLists(poles, coeffs, self.data.Psupp, basis,
+                                 matchingWeight, HFEngine, is_state)
+
     def initializeFromLists(self, poles:ListAny, coeffs:ListAny, supps:ListAny,
                             basis:str, matchingWeight:float, HFEngine:HFEng,
                             is_state:bool):
         """Initialize Heaviside representation."""
         Ns = [len(pls) for pls in poles]
         poles, coeffs = heavisideUniformShape(poles, coeffs)
         root = Ns.index(len(poles[0]))
         if hasattr(self.data, "_is_C_quadratic") and self.data._is_C_quadratic:
             csizemax = np.max([len(c) for c in coeffs])
             #csizemax = np.max([c.shape[1] for c in coeffs])
             if self.data._is_C_quadratic == 1:
                 csizemax = csizemax ** 2
             else: # if self.data._is_C_quadratic == 2:
                 csizemax = csizemax * (csizemax + 1) // 2
             TrainedModelRational._setupQuadMapping(self, csizemax)
             projMapping = self.quad_mapping
             projMappingReal = self.data._is_C_quadratic == 2
         else:
             projMapping, projMappingReal = None, False
         poles, coeffs = rationalFunctionMatching(poles, coeffs,
                                                  self.data.musMarginal.data,
                                                  matchingWeight, supps,
                                                  self.data.projMat, HFEngine,
                                                  is_state, root, projMapping,
                                                  projMappingReal)
-        super().initializeFromLists(poles, coeffs, supps, basis)
+        self.data.HIs = []
+        for pls, cfs in zip(poles, coeffs):
+            hsi = HI()
+            hsi.poles = pls
+            if len(cfs) == len(pls):
+                cfs = np.pad(cfs, ((0, 1), (0, 0)), "constant")
+            hsi.coeffs = cfs
+            hsi.npar = 1
+            hsi.polybasis = basis
+            self.data.HIs += [hsi]
         self.data.suppEffPts = [np.arange(len(self.data.HIs))]
         self.data.suppEffIdx = np.zeros(len(poles[0]), dtype = int)
 
     def checkSharedRatio(self, shared:float) -> str:
         N = len(self.data.HIs[0].poles)
         M = len(self.data.HIs)
         goodLocPoles = np.array([np.logical_not(np.isinf(hi.poles)
                                                     ) for hi in self.data.HIs])
         self.data.suppEffPts = [np.arange(len(self.data.HIs))]
         self.data.suppEffIdx = np.zeros(N, dtype = int)
         if np.all(np.all(goodLocPoles)): return " No poles erased."
         goodGlobPoles = np.sum(goodLocPoles, axis = 0)
         goodEnoughPoles = goodGlobPoles >= max(1., 1. * shared * M)
         keepPole = np.where(goodEnoughPoles)[0]
         halfPole = np.where(np.logical_and(goodEnoughPoles,
                                            goodGlobPoles < M))[0]
         removePole = np.where(np.logical_not(goodEnoughPoles))[0]
         if len(removePole) > 0:
             keepCoeff = np.append(keepPole,
                                   np.arange(N, len(self.data.HIs[0].coeffs)))
             for hi in self.data.HIs:
                 for j in removePole:
                     if not np.isinf(hi.poles[j]):
                         hi.coeffs[N, :] -= hi.coeffs[j, :] / hi.poles[j]
                 hi.poles = hi.poles[keepPole]
                 hi.coeffs = hi.coeffs[keepCoeff, :]
         for idxR in halfPole:
             pts = np.where(goodLocPoles[:, idxR])[0]
             idxEff = len(self.data.suppEffPts)
             for idEff, prevPts in enumerate(self.data.suppEffPts):
                 if len(prevPts) == len(pts):
                     if np.allclose(prevPts, pts):
                         idxEff = idEff
                         break
             if idxEff == len(self.data.suppEffPts):
                 self.data.suppEffPts += [pts]
             self.data.suppEffIdx[idxR] = idxEff
         self.data.suppEffIdx = self.data.suppEffIdx[keepPole]
         return (" Hard-erased {} pole".format(len(removePole))
               + "s" * (len(removePole) != 1)
               + " and soft-erased {} pole".format(len(halfPole))
               + "s" * (len(halfPole) != 1) + ".")
 
+    def getApproxReduced(self, mu : paramList = []) -> sampList:
+        """
+        Evaluate reduced representation of approximant at arbitrary parameter.
+
+        Args:
+            mu: Target parameter.
+        """
+        RROMPyAssert(self.data.nparPivot, 1, "Number of pivot parameters")
+        mu = self.checkParameterList(mu)
+        if (not hasattr(self, "lastSolvedApproxReduced")
+         or self.lastSolvedApproxReduced != mu):
+            vbMng(self, "INIT",
+                  "Evaluating approximant at mu = {}.".format(mu), 12)
+            muP = self.centerNormalizePivot(mu(self.data.directionPivot))
+            muM = mu(self.data.directionMarginal)
+            his = self.interpolateMarginalInterpolator(muM)
+            for i, (mP, hi) in enumerate(zip(muP, his)):
+                uAppR = hi(mP)[:, 0]
+                if i == 0:
+                    uApproxR = np.empty((len(uAppR), len(mu)),
+                                        dtype = uAppR.dtype)
+                uApproxR[:, i] = uAppR
+            self.uApproxReduced = sampleList(uApproxR)
+            vbMng(self, "DEL", "Done evaluating approximant.", 12)
+            self.lastSolvedApproxReduced = mu
+        return self.uApproxReduced
+
     def interpolateMarginalInterpolator(self, mu : paramList = []) -> ListAny:
         """Obtain interpolated approximant interpolator."""
         mu = self.checkParameterListMarginal(mu)
         vbMng(self, "INIT",
               "Interpolating marginal models at mu = {}.".format(mu), 95)
         his = []
         muC = self.centerNormalizeMarginal(mu)
         mIvals = self.data.marginalInterp(muC)
         verb, self.verbosity = self.verbosity, 0
         poless = self.interpolateMarginalPoles(mu, mIvals)
         coeffss = self.interpolateMarginalCoeffs(mu, mIvals)
         self.verbosity = verb
         for j in range(len(mu)):
             his += [HI()]
             his[-1].poles = poless[j]
             his[-1].coeffs = coeffss[j]
             his[-1].npar = 1
             his[-1].polybasis = self.data.HIs[0].polybasis
         vbMng(self, "DEL", "Done interpolating marginal models.", 95)
         return his
 
     def interpolateMarginalPoles(self, mu : paramList = [],
                                  mIvals : Np2D = None) -> ListAny:
         """Obtain interpolated approximant poles."""
         mu = self.checkParameterListMarginal(mu)
         vbMng(self, "INIT",
               "Interpolating marginal poles at mu = {}.".format(mu), 95)
         intMPoles = np.zeros((len(mu),) + self.data.polesEff[0].shape,
                              dtype = self.data.polesEff[0].dtype)
         if mIvals is None:
             muC = self.centerNormalizeMarginal(mu)
             mIvals = self.data.marginalInterp(muC)
         for pEff, mI in zip(self.data.polesEff, mIvals):
             intMPoles += np.expand_dims(mI, - 1) * pEff
         vbMng(self, "DEL", "Done interpolating marginal poles.", 95)
         return intMPoles
 
     def interpolateMarginalCoeffs(self, mu : paramList = [],
                                   mIvals : Np2D = None) -> ListAny:
         """Obtain interpolated approximant coefficients."""
         mu = self.checkParameterListMarginal(mu)
         vbMng(self, "INIT",
               "Interpolating marginal coefficients at mu = {}.".format(mu), 95)
         intMCoeffs = np.zeros((len(mu),) + self.data.coeffsEff[0].shape,
                               dtype = self.data.coeffsEff[0].dtype)
         if mIvals is None:
             muC = self.centerNormalizeMarginal(mu)
             mIvals = self.data.marginalInterp(muC)
         for cEff, mI in zip(self.data.coeffsEff, mIvals):
             for j, m in enumerate(mI): intMCoeffs[j] += m * cEff
         vbMng(self, "DEL", "Done interpolating marginal coefficients.", 95)
         return intMCoeffs
+
+    def getPVal(self, mu : paramList = []) -> sampList:
+        """
+        Evaluate rational numerator at arbitrary parameter.
+
+        Args:
+            mu: Target parameter.
+        """
+        RROMPyAssert(self.data.nparPivot, 1, "Number of pivot parameters")
+        mu = self.checkParameterList(mu)
+        p = emptySampleList()
+        muP = self.centerNormalizePivot(mu(self.data.directionPivot))
+        muM = mu(self.data.directionMarginal)
+        his = self.interpolateMarginalInterpolator(muM)
+        for i, (mP, hi) in enumerate(zip(muP, his)):
+            Pval = hi(mP) * np.prod(mP[0] - hi.poles)
+            if i == 0: p.reset((len(Pval), len(mu)), dtype = Pval.dtype)
+            p[i] = Pval
+        return p
+
+    def getQVal(self, mu:Np1D, der : List[int] = None,
+                scl : Np1D = None) -> Np1D:
+        """
+        Evaluate rational denominator at arbitrary parameter.
+
+        Args:
+            mu: Target parameter.
+            der(optional): Derivatives to take before evaluation.
+        """
+        RROMPyAssert(self.data.nparPivot, 1, "Number of pivot parameters")
+        mu = self.checkParameterList(mu)
+        muP = self.centerNormalizePivot(mu(self.data.directionPivot))
+        muM = mu(self.data.directionMarginal)
+        if der is None:
+            derP, derM = 0, [0]
+        else:
+            derP = der[self.data.directionPivot[0]]
+            derM = [der[x] for x in self.data.directionMarginal]
+        if np.any(np.array(derM) != 0):
+            raise RROMPyException(("Derivatives of Q with respect to marginal "
+                                   "parameters not allowed."))
+        sclP = 1 if scl is None else scl[self.data.directionPivot[0]]
+        derVal = np.zeros(len(mu), dtype = np.complex)
+        pls = self.interpolateMarginalPoles(muM)
+        for i, (mP, pl) in enumerate(zip(muP, pls)):
+            N = len(pl)
+            if derP == N: derVal[i] = 1.
+            elif derP >= 0 and derP < N:
+                plDist = mP[0] - pl
+                for terms in combinations(np.arange(N), N - derP):
+                    derVal[i] += np.prod(plDist[list(terms)])
+        return sclP ** derP * fact(derP) * derVal
+
+    def getPoles(self, *args, **kwargs) -> Np1D:
+        """
+        Obtain approximant poles.
+
+        Returns:
+            Numpy complex vector of poles.
+        """
+        RROMPyAssert(self.data.nparPivot, 1, "Number of pivot parameters")
+        if len(args) + len(kwargs) > 1:
+            raise RROMPyException(("Wrong number of parameters passed. "
+                                   "Only 1 available."))
+        elif len(args) + len(kwargs) == 1:
+            if len(args) == 1:
+                mVals = args[0]
+            else:
+                mVals = kwargs["marginalVals"]
+            if not isinstance(mVals, Iterable): mVals = [mVals]
+            mVals = list(mVals)
+        else:
+            mVals = [fp]
+        rDim = mVals.index(fp)
+        if rDim < len(mVals) - 1 and fp in mVals[rDim + 1 :]:
+            raise RROMPyException(("Exactly 1 'freepar' entry in "
+                                   "marginalVals must be provided."))
+        if rDim != self.data.directionPivot[0]:
+            raise RROMPyException(("'freepar' entry in marginalVals must "
+                                   "coincide with pivot direction."))
+        mVals[rDim] = self.data.mu0(rDim)[0]
+        mMarg = [mVals[j] for j in range(len(mVals)) if j != rDim]
+        roots = (self.data.scaleFactor[rDim]
+               * self.interpolateMarginalPoles(mMarg)[0])
+        return self.mapParameterList(self.mapParameterList(self.data.mu0(rDim),
+                                                           idx = [rDim])(0, 0)
+                                   + roots, "B", [rDim])(0)
+
+    def getResidues(self, *args, **kwargs) -> Np2D:
+        """
+        Obtain approximant residues.
+
+        Returns:
+            Numpy matrix with residues as columns.
+        """
+        pls = self.getPoles(*args, **kwargs)
+        if len(args) == 1:
+            mVals = args[0]
+        elif len(args) == 0:
+            mVals = [None]
+        else:
+            mVals = kwargs["marginalVals"]
+        if not isinstance(mVals, Iterable): mVals = [mVals]
+        mVals = list(mVals)
+        rDim = mVals.index(fp)
+        mMarg = [mVals[j] for j in range(len(mVals)) if j != rDim]
+        res = self.interpolateMarginalCoeffs(mMarg)[0][: len(pls), :].T
+        if hasattr(self.data, "_is_C_quadratic") and self.data._is_C_quadratic:
+            self._setupQuadMapping()
+            res = res[self.quad_mapping[0]] * res[self.quad_mapping[1]].conj()
+        if not self.data._collapsed: res = dot(self.data.projMat, res).T
+        if (hasattr(self.data, "_is_C_quadratic")
+        and self.data._is_C_quadratic == 2):
+            res = np.real(res)
+        return pls, res
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 8550521..8c250e8 100644
--- a/tests/4_reduction_methods_multiD/greedy_pivoted_rational_2d.py
+++ b/tests/4_reduction_methods_multiD/greedy_pivoted_rational_2d.py
@@ -1,85 +1,85 @@
 # Copyright (C) 2018-2020 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)
+                      RationalInterpolantPivotedGreedyPoleMatch as RIPG,
+                      RationalInterpolantGreedyPivotedGreedyPoleMatch 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": "LOOK_AHEAD_RECOVER",
               "matchingWeight": 1., "samplerMarginal":SGS([6.75, 7.25])}
     approx = RIPG([0], solver, mu0, 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": "LOOK_AHEAD_RECOVER",
               "matchingWeight": 1., "samplerMarginal":SGS([6.75, 7.25])}
     approx = RIGPG([0], solver, mu0, 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/4_reduction_methods_multiD/pivoted_rational_2d.py b/tests/4_reduction_methods_multiD/pivoted_rational_2d.py
index 4552ba7..15c2c35 100644
--- a/tests/4_reduction_methods_multiD/pivoted_rational_2d.py
+++ b/tests/4_reduction_methods_multiD/pivoted_rational_2d.py
@@ -1,111 +1,112 @@
 # Copyright (C) 2018-2020 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 (RationalInterpolantPivoted as RIP,
-                                      RationalInterpolantGreedyPivoted as RIGP)
+from rrompy.reduction_methods import (
+                             RationalInterpolantPivotedPoleMatch as RIP,
+                             RationalInterpolantGreedyPivotedPoleMatch 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, "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, 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, "S": 5, "polybasis": "MONOMIAL",
               "samplerPivot": MS([4.75, 5.25], np.array([5.]),
                                  normalFoci = [0., 0.]), "SMarginal": 5,
               "polybasisMarginal": "MONOMIAL", "matchingWeight": 1.,
               "samplerMarginal": MS([6.75, 7.25], np.linspace(6.75, 7.25, 5)),
               "robustTol": 1e-6, "interpRcond": 1e-3}
     approx = RIP([0], solver, mu0, 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"),
               "SMarginal": 5, "polybasisMarginal": "MONOMIAL",
               "samplerMarginal": QS([6.75, 7.25], "UNIFORM"),
               "matchingWeight": 1.}
     solver.cutOffPolesRMinRel, solver.cutOffPolesRMaxRel = -3., 3.
     solver.cutOffPolesIMinRel, solver.cutOffPolesIMaxRel = -1.5, 1.5
     approx = RIGP([0], solver, mu0, 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)