diff --git a/examples/2d/pod/plot_inf_set.py b/examples/2d/pod/plot_inf_set.py
index a57dd29..5f5c891 100644
--- a/examples/2d/pod/plot_inf_set.py
+++ b/examples/2d/pod/plot_inf_set.py
@@ -1,323 +1,335 @@
 import warnings
 import numpy as np
 from matplotlib import pyplot as plt
 
 def plotInfSet1FromData(mus, Z, T, R, E, beta, murange, approx, mu0,
                         exp = 2., normalizeDen = False):
     if hasattr(approx, "mus"):
         mu2x = approx.mus(0) ** exp
     else:
         mu2x = mu0[0] ** exp
     murangeExp = [[murange[0][0] ** exp], [murange[1][0] ** exp]]
     mu1 = np.real(np.power(mus, exp))
-    eta = R / beta / E
     ZTmin, ZTmax = min(np.min(Z), np.min(T)), max(np.max(Z), np.max(T))
     Rmin, Rmax = np.min(R), np.max(R)
     Emin, Emax = np.min(E), np.max(E)
-    betamin, betamax = np.min(beta), np.max(beta)
-    etamin, etamax = np.min(eta), np.max(eta)
+    if not np.isnan(beta[0]):
+        eta = R / beta / E
+        betamin, betamax = np.min(beta), np.max(beta)
+        etamin, etamax = np.min(eta), np.max(eta)
     plt.figure(figsize = (15, 7))
     plt.jet()
     plt.semilogy(mu1, Z)
     plt.semilogy(mu1, T, '--')
     for l_ in approx.trainedModel.getPoles():
         plt.plot([np.real(l_ ** exp)] * 2, [ZTmin, ZTmax], 'b:')
     plt.plot(mu2x, [ZTmin] * len(mu2x), 'kx')
     plt.plot([murangeExp[0][0]] * 2, [ZTmin, ZTmax], 'm:')
     plt.plot([murangeExp[1][0]] * 2, [ZTmin, ZTmax], 'm:')
     plt.xlim(mu1[0], mu1[-1])
     plt.title("|u(mu)|, |u_app(mu)|")
     plt.grid()
     plt.show()
 
     plt.figure(figsize = (15, 7))
     plt.jet()
     plt.semilogy(mu1, R)
     for l_ in approx.trainedModel.getPoles():
         plt.plot([np.real(l_ ** exp)] * 2, [Rmin, Rmax], 'b:')
     plt.plot(mu2x, [Rmax] * len(mu2x), 'kx')
     plt.plot([murangeExp[0][0]] * 2, [Rmin, Rmax], 'm:')
     plt.plot([murangeExp[1][0]] * 2, [Rmin, Rmax], 'm:')
     plt.xlim(mu1[0], mu1[-1])
     if normalizeDen:
         plt.title("|Q(mu)res(mu)|")
     else:
         plt.title("|res(mu)|")
     plt.grid()
     plt.show()
     
     plt.figure(figsize = (15, 7))
     plt.jet()
     plt.semilogy(mu1, E)
     for l_ in approx.trainedModel.getPoles():
         plt.plot([np.real(l_ ** exp)] * 2, [Emin, Emax], 'b:')
     plt.plot(mu2x, [Emax] * len(mu2x), 'kx')
     plt.plot([murangeExp[0][0]] * 2, [Emin, Emax], 'm:')
     plt.plot([murangeExp[1][0]] * 2, [Emin, Emax], 'm:')
     plt.xlim(mu1[0], mu1[-1])
     if normalizeDen:
         plt.title("|Q(mu)err(mu)|")
     else:
         plt.title("|err(mu)|")
     plt.grid()
     plt.show()
     
-    plt.figure(figsize = (15, 7))
-    plt.jet()
-    plt.plot(mu1, beta)
-    for l_ in approx.trainedModel.getPoles():
-        plt.plot([np.real(l_ ** exp)] * 2, [betamin, betamax], 'b:')
-    plt.plot(mu2x, [betamax] * len(mu2x), 'kx')
-    plt.plot([murangeExp[0][0]] * 2, [betamin, betamax], 'm:')
-    plt.plot([murangeExp[1][0]] * 2, [betamin, betamax], 'm:')
-    plt.xlim(mu1[0], mu1[-1])
-    plt.title("beta(mu)")
-    plt.grid()
-    plt.show()
+    if not np.isnan(beta[0]):
+        plt.figure(figsize = (15, 7))
+        plt.jet()
+        plt.plot(mu1, beta)
+        for l_ in approx.trainedModel.getPoles():
+            plt.plot([np.real(l_ ** exp)] * 2, [betamin, betamax], 'b:')
+        plt.plot(mu2x, [betamax] * len(mu2x), 'kx')
+        plt.plot([murangeExp[0][0]] * 2, [betamin, betamax], 'm:')
+        plt.plot([murangeExp[1][0]] * 2, [betamin, betamax], 'm:')
+        plt.xlim(mu1[0], mu1[-1])
+        plt.title("beta(mu)")
+        plt.grid()
+        plt.show()
     
-    plt.figure(figsize = (15, 7))
-    plt.jet()
-    plt.semilogy(mu1, R / beta)
-    plt.semilogy(mu1, E, '--')
-    for l_ in approx.trainedModel.getPoles():
-        plt.plot([np.real(l_ ** exp)] * 2, [Emin, Emax], 'b:')
-    plt.plot(mu2x, [Emax] * len(mu2x), 'kx')
-    plt.plot([murangeExp[0][0]] * 2, [Emin, Emax], 'm:')
-    plt.plot([murangeExp[1][0]] * 2, [Emin, Emax], 'm:')
-    plt.xlim(mu1[0], mu1[-1])
-    if normalizeDen:
-        plt.title("|Q(mu)res(mu)/beta(mu)|")
-    else:
-        plt.title("|res(mu)/beta(mu)|")
-    plt.grid()
-    plt.show()
-    
-    plt.figure(figsize = (15, 7))
-    plt.jet()
-    plt.semilogy(mu1, eta)
-    for l_ in approx.trainedModel.getPoles():
-        plt.plot([np.real(l_ ** exp)] * 2, [etamin, etamax], 'b:')
-    plt.plot(mu2x, [etamax] * len(mu2x), 'kx')
-    plt.plot([murangeExp[0][0]] * 2, [etamin, etamax], 'm:')
-    plt.plot([murangeExp[1][0]] * 2, [etamin, etamax], 'm:')
-    plt.xlim(mu1[0], mu1[-1])
-    plt.title("eta(mu)")
-    plt.grid()
-    plt.show()
+        plt.figure(figsize = (15, 7))
+        plt.jet()
+        plt.semilogy(mu1, R / beta)
+        plt.semilogy(mu1, E, '--')
+        for l_ in approx.trainedModel.getPoles():
+            plt.plot([np.real(l_ ** exp)] * 2, [Emin, Emax], 'b:')
+        plt.plot(mu2x, [Emax] * len(mu2x), 'kx')
+        plt.plot([murangeExp[0][0]] * 2, [Emin, Emax], 'm:')
+        plt.plot([murangeExp[1][0]] * 2, [Emin, Emax], 'm:')
+        plt.xlim(mu1[0], mu1[-1])
+        if normalizeDen:
+            plt.title("|Q(mu)res(mu)/beta(mu)|")
+        else:
+            plt.title("|res(mu)/beta(mu)|")
+        plt.grid()
+        plt.show()
+        
+        plt.figure(figsize = (15, 7))
+        plt.jet()
+        plt.semilogy(mu1, eta)
+        for l_ in approx.trainedModel.getPoles():
+            plt.plot([np.real(l_ ** exp)] * 2, [etamin, etamax], 'b:')
+        plt.plot(mu2x, [etamax] * len(mu2x), 'kx')
+        plt.plot([murangeExp[0][0]] * 2, [etamin, etamax], 'm:')
+        plt.plot([murangeExp[1][0]] * 2, [etamin, etamax], 'm:')
+        plt.xlim(mu1[0], mu1[-1])
+        plt.title("eta(mu)")
+        plt.grid()
+        plt.show()
     
 def plotInfSet1(murange, murangeEff, approx, mu0, nSamples = 20, exp = 2.,
-                relative = True, normalizeDen = False):
+                relative = True, normalizeDen = False, nobeta = False):
     mu1 = np.linspace(murangeEff[0][0] ** exp, murangeEff[1][0] ** exp,
                       nSamples)
     mus = np.power(mu1, 1. / exp)
     Z = approx.normHF(mus)
     T = approx.normApprox(mus)
     R = approx.normRes(mus)
     E = approx.normErr(mus)
     if relative:
         F = approx.normRHS(mus)
         R /= F
         E /= Z
     if normalizeDen:
         Qvals = np.abs(approx.trainedModel.getQVal(mus))
         R *= Qvals
         E *= Qvals
-    beta = approx.HFEngine.stabilityFactor(approx.getHF(mus), mus)
+    if nobeta:
+        beta = np.empty(len(mus))
+        beta[:] = np.nan
+    else:
+        beta = approx.HFEngine.stabilityFactor(approx.getHF(mus), mus)
     plotInfSet1FromData(mus, Z, T, R, E, beta, murange, approx, mu0,
                         exp, normalizeDen)
     return mus, Z, T, R, E, beta
     
 def plotInfSet2FromData(mus, Ze, Te, Re, Ee, beta, murange, approx, mu0,
                         exps = [2., 2.], clip = -1, normalizeDen = False):
     if hasattr(approx, "mus"):
         mu2x, mu2y = approx.mus(0) ** exps[0], approx.mus(1) ** exps[1]
     else:
         mu2x, mu2y = mu0[0] ** exps[0], mu0[1] ** exps[1]
     murangeExp = [[murange[0][0] ** exps[0], murange[0][1] ** exps[1]],
                   [murange[1][0] ** exps[0], murange[1][1] ** exps[1]]]
     mu1s = np.unique([m[0] for m in mus])
     mu2s = np.unique([m[1] for m in mus])
     mu1 = np.power(mu1s, exps[0])
     mu2 = np.power(mu2s, exps[1])
     Mu1, Mu2 = np.meshgrid(np.real(mu1), np.real(mu2))
 
     Z = np.log10(Ze)
     T = np.log10(Te)
     R = np.log10(Re)
     E = np.log10(Ee)
     ZTmin, ZTmax = min(np.min(Z), np.min(T)), max(np.max(Z), np.max(T))
     Rmin, Rmax = np.min(R), np.max(R)
     Emin, Emax = np.min(E), np.max(E)
-    betamin, betamax = np.min(beta), np.max(beta)
+    if not np.isnan(beta[0, 0]):
+        betamin, betamax = np.min(beta), np.max(beta)
     if clip > 0:
         ZTmax -= clip * (ZTmax - ZTmin)
         cmap = plt.cm.bone
     else:
         cmap = plt.cm.jet
     warnings.simplefilter("ignore", category = UserWarning)
     plt.figure(figsize = (15, 7))
     plt.jet()
     p = plt.contourf(Mu1, Mu2, Z, cmap = cmap,
                      levels = np.linspace(ZTmin, ZTmax, 50))
     if clip > 0:
         plt.contour(Mu1, Mu2, Z, [ZTmin])
     plt.plot(mu2x, mu2y, 'kx')
     plt.plot([murangeExp[0][0]] * 2,
              [murangeExp[0][1], murangeExp[1][1]], 'm:')
     plt.plot([murangeExp[0][0], murangeExp[1][0]],
              [murangeExp[1][1]] * 2, 'm:')
     plt.plot([murangeExp[1][0]] * 2,
              [murangeExp[1][1], murangeExp[0][1]], 'm:')
     plt.plot([murangeExp[1][0], murangeExp[0][0]],
              [murangeExp[0][1]] * 2, 'm:')
     plt.colorbar(p)
     plt.title("log10|u(mu)|")
     plt.grid()
     plt.show()
 
     plt.figure(figsize = (15, 7))
     plt.jet()
     p = plt.contourf(Mu1, Mu2, T, cmap = cmap,
                      levels = np.linspace(ZTmin, ZTmax, 50))
     if clip > 0:
         plt.contour(Mu1, Mu2, T, [ZTmin])
     plt.plot(mu2x, mu2y, 'kx')
     plt.plot([murangeExp[0][0]] * 2,
              [murangeExp[0][1], murangeExp[1][1]], 'm:')
     plt.plot([murangeExp[0][0], murangeExp[1][0]],
              [murangeExp[1][1]] * 2, 'm:')
     plt.plot([murangeExp[1][0]] * 2,
              [murangeExp[1][1], murangeExp[0][1]], 'm:')
     plt.plot([murangeExp[1][0], murangeExp[0][0]],
              [murangeExp[0][1]] * 2, 'm:')
     plt.title("log10|u_app(mu)|")
     plt.colorbar(p)
     plt.grid()
     plt.show()
     
     plt.figure(figsize = (15, 7))
     plt.jet()
     p = plt.contourf(Mu1, Mu2, R, levels = np.linspace(Rmin, Rmax, 50))
     plt.plot(mu2x, mu2y, 'kx')
     plt.plot([murangeExp[0][0]] * 2,
              [murangeExp[0][1], murangeExp[1][1]], 'm:')
     plt.plot([murangeExp[0][0], murangeExp[1][0]],
              [murangeExp[1][1]] * 2, 'm:')
     plt.plot([murangeExp[1][0]] * 2,
              [murangeExp[1][1], murangeExp[0][1]], 'm:')
     plt.plot([murangeExp[1][0], murangeExp[0][0]],
              [murangeExp[0][1]] * 2, 'm:')
     if normalizeDen:
         plt.title("log10|Q(mu)res(mu)|")
     else:
         plt.title("log10|res(mu)|")
     plt.colorbar(p)
     plt.grid()
     plt.show()
     
     plt.figure(figsize = (15, 7))
     plt.jet()
     p = plt.contourf(Mu1, Mu2, E, levels = np.linspace(Emin, Emax, 50))
     plt.plot(mu2x, mu2y, 'kx')
     plt.plot([murangeExp[0][0]] * 2,
              [murangeExp[0][1], murangeExp[1][1]], 'm:')
     plt.plot([murangeExp[0][0], murangeExp[1][0]],
              [murangeExp[1][1]] * 2, 'm:')
     plt.plot([murangeExp[1][0]] * 2,
              [murangeExp[1][1], murangeExp[0][1]], 'm:')
     plt.plot([murangeExp[1][0], murangeExp[0][0]],
              [murangeExp[0][1]] * 2, 'm:')
     if normalizeDen:
         plt.title("log10|Q(mu)err(mu)|")
     else:
         plt.title("log10|err(mu)|")
     plt.colorbar(p)
     plt.grid()
     plt.show()
     
-    plt.figure(figsize = (15, 7))
-    plt.jet()
-    p = plt.contourf(Mu1, Mu2, beta,
-                     levels = np.linspace(betamin, betamax, 50))
-    plt.plot(mu2x, mu2y, 'kx')
-    plt.plot([murangeExp[0][0]] * 2,
-             [murangeExp[0][1], murangeExp[1][1]], 'm:')
-    plt.plot([murangeExp[0][0], murangeExp[1][0]],
-             [murangeExp[1][1]] * 2, 'm:')
-    plt.plot([murangeExp[1][0]] * 2,
-             [murangeExp[1][1], murangeExp[0][1]], 'm:')
-    plt.plot([murangeExp[1][0], murangeExp[0][0]],
-             [murangeExp[0][1]] * 2, 'm:')
-    plt.title("beta(mu)")
-    plt.colorbar(p)
-    plt.grid()
-    plt.show()
-    
-    plt.figure(figsize = (15, 7))
-    plt.jet()
-    p = plt.contourf(Mu1, Mu2, R - np.log10(beta),
-                     levels = np.linspace(Emin, Emax, 50))
-    plt.plot(mu2x, mu2y, 'kx')
-    plt.plot([murangeExp[0][0]] * 2,
-             [murangeExp[0][1], murangeExp[1][1]], 'm:')
-    plt.plot([murangeExp[0][0], murangeExp[1][0]],
-             [murangeExp[1][1]] * 2, 'm:')
-    plt.plot([murangeExp[1][0]] * 2,
-             [murangeExp[1][1], murangeExp[0][1]], 'm:')
-    plt.plot([murangeExp[1][0], murangeExp[0][0]],
-             [murangeExp[0][1]] * 2, 'm:')
-    if normalizeDen:
-        plt.title("log10|Q(mu)res(mu)/beta(mu)|")
-    else:
-        plt.title("log10|res(mu)/beta(mu)|")
-    plt.colorbar(p)
-    plt.grid()
-    plt.show()
-    
-    plt.figure(figsize = (15, 7))
-    plt.jet()
-    p = plt.contourf(Mu1, Mu2, R - np.log10(beta) - E, 50)
-    plt.plot(mu2x, mu2y, 'kx')
-    plt.plot([murangeExp[0][0]] * 2,
-             [murangeExp[0][1], murangeExp[1][1]], 'm:')
-    plt.plot([murangeExp[0][0], murangeExp[1][0]],
-             [murangeExp[1][1]] * 2, 'm:')
-    plt.plot([murangeExp[1][0]] * 2,
-             [murangeExp[1][1], murangeExp[0][1]], 'm:')
-    plt.plot([murangeExp[1][0], murangeExp[0][0]],
-             [murangeExp[0][1]] * 2, 'm:')
-    plt.title("log10|eta(mu)|")
-    plt.colorbar(p)
-    plt.grid()
-    plt.show()
+    if not np.isnan(beta[0, 0]):
+        plt.figure(figsize = (15, 7))
+        plt.jet()
+        p = plt.contourf(Mu1, Mu2, beta,
+                         levels = np.linspace(betamin, betamax, 50))
+        plt.plot(mu2x, mu2y, 'kx')
+        plt.plot([murangeExp[0][0]] * 2,
+                 [murangeExp[0][1], murangeExp[1][1]], 'm:')
+        plt.plot([murangeExp[0][0], murangeExp[1][0]],
+                 [murangeExp[1][1]] * 2, 'm:')
+        plt.plot([murangeExp[1][0]] * 2,
+                 [murangeExp[1][1], murangeExp[0][1]], 'm:')
+        plt.plot([murangeExp[1][0], murangeExp[0][0]],
+                 [murangeExp[0][1]] * 2, 'm:')
+        plt.title("beta(mu)")
+        plt.colorbar(p)
+        plt.grid()
+        plt.show()
+        
+        plt.figure(figsize = (15, 7))
+        plt.jet()
+        p = plt.contourf(Mu1, Mu2, R - np.log10(beta),
+                         levels = np.linspace(Emin, Emax, 50))
+        plt.plot(mu2x, mu2y, 'kx')
+        plt.plot([murangeExp[0][0]] * 2,
+                 [murangeExp[0][1], murangeExp[1][1]], 'm:')
+        plt.plot([murangeExp[0][0], murangeExp[1][0]],
+                 [murangeExp[1][1]] * 2, 'm:')
+        plt.plot([murangeExp[1][0]] * 2,
+                 [murangeExp[1][1], murangeExp[0][1]], 'm:')
+        plt.plot([murangeExp[1][0], murangeExp[0][0]],
+                 [murangeExp[0][1]] * 2, 'm:')
+        if normalizeDen:
+            plt.title("log10|Q(mu)res(mu)/beta(mu)|")
+        else:
+            plt.title("log10|res(mu)/beta(mu)|")
+        plt.colorbar(p)
+        plt.grid()
+        plt.show()
+        
+        plt.figure(figsize = (15, 7))
+        plt.jet()
+        p = plt.contourf(Mu1, Mu2, R - np.log10(beta) - E, 50)
+        plt.plot(mu2x, mu2y, 'kx')
+        plt.plot([murangeExp[0][0]] * 2,
+                 [murangeExp[0][1], murangeExp[1][1]], 'm:')
+        plt.plot([murangeExp[0][0], murangeExp[1][0]],
+                 [murangeExp[1][1]] * 2, 'm:')
+        plt.plot([murangeExp[1][0]] * 2,
+                 [murangeExp[1][1], murangeExp[0][1]], 'm:')
+        plt.plot([murangeExp[1][0], murangeExp[0][0]],
+                 [murangeExp[0][1]] * 2, 'm:')
+        plt.title("log10|eta(mu)|")
+        plt.colorbar(p)
+        plt.grid()
+        plt.show()
     
 def plotInfSet2(murange, murangeEff, approx, mu0, nSamples = 200,
                 exps = [2., 2.], clip = -1, relative = True,
-                normalizeDen = False):
+                normalizeDen = False, nobeta = False):
     mu1 = np.linspace(murangeEff[0][0] ** exps[0], murangeEff[1][0] ** exps[0],
                       nSamples)
     mu2 = np.linspace(murangeEff[0][1] ** exps[1], murangeEff[1][1] ** exps[1],
                       nSamples)
     mu1s = np.power(mu1, 1. / exps[0])
     mu2s = np.power(mu2, 1. / exps[1])
     mus = [(m1, m2) for m2 in mu2s for m1 in mu1s]
 
     Ze = approx.normHF(mus).reshape((nSamples, nSamples))
     Te = approx.normApprox(mus).reshape((nSamples, nSamples))
     Re = approx.normRes(mus).reshape((nSamples, nSamples))
     Ee = approx.normErr(mus).reshape((nSamples, nSamples))
     if relative:
         Fe = approx.normRHS(mus).reshape((nSamples, nSamples))
         Re /= Fe
         Ee /= Ze
     if normalizeDen:
         Qvals = np.abs(approx.trainedModel.getQVal(mus).reshape(
                                                          (nSamples, nSamples)))
         Re *= Qvals
         Ee *= Qvals
-    betae = approx.HFEngine.stabilityFactor(approx.getHF(mus), mus).reshape(
-                                                          (nSamples, nSamples))
+    if nobeta:
+        betae = np.empty((nSamples, nSamples))
+        betae[:, :] = np.nan
+    else:
+        betae = approx.HFEngine.stabilityFactor(approx.getHF(mus), mus)\
+                                                 .reshape((nSamples, nSamples))
     plotInfSet2FromData(mus, Ze, Te, Re, Ee, betae, murange,
                         approx, mu0, exps, clip, normalizeDen)
     return mus, Ze, Te, Re, Ee, betae
 
diff --git a/examples/3d/pod/plot_inf_set_3.py b/examples/3d/pod/plot_inf_set_3.py
index 880b562..4ca96db 100644
--- a/examples/3d/pod/plot_inf_set_3.py
+++ b/examples/3d/pod/plot_inf_set_3.py
@@ -1,347 +1,359 @@
 import warnings
 from copy import deepcopy as copy
 import numpy as np
 from matplotlib import pyplot as plt
 
 def plotInfSet31FromData(mNone, mus, Z, T, R, E, beta, murange, approx, mu0,
                          exp = 2., normalizeDen = False):
     if hasattr(approx, "mus"):
         mu2x = approx.mus(mNone[0]) ** exp
     else:
         mu2x = mu0[0] ** exp
     murangeExp = [[murange[0][mNone[0]] ** exp],
                   [murange[1][mNone[0]] ** exp]]
     mu1 = np.real(np.power([m[mNone[0]] for m in mus], exp))
-    eta = R / beta / E
     ZTmin, ZTmax = min(np.min(Z), np.min(T)), max(np.max(Z), np.max(T))
     Rmin, Rmax = np.min(R), np.max(R)
     Emin, Emax = np.min(E), np.max(E)
-    betamin, betamax = np.min(beta), np.max(beta)
-    etamin, etamax = np.min(eta), np.max(eta)
+    if not np.isnan(beta[0]):
+        eta = R / beta / E
+        betamin, betamax = np.min(beta), np.max(beta)
+        etamin, etamax = np.min(eta), np.max(eta)
     plt.figure(figsize = (15, 7))
     plt.jet()
     plt.semilogy(mu1, Z)
     plt.semilogy(mu1, T, '--')
     mTMP = mus[0]
     mTMP[mNone[0]] = None
     for l_ in approx.trainedModel.getPoles(mTMP):
         plt.plot([np.real(l_ ** exp)] * 2, [ZTmin, ZTmax], 'b:')
     plt.plot(mu2x, [ZTmin] * len(mu2x), 'kx')
     plt.plot([murangeExp[0][0]] * 2, [ZTmin, ZTmax], 'm:')
     plt.plot([murangeExp[1][0]] * 2, [ZTmin, ZTmax], 'm:')
     plt.xlim(mu1[0], mu1[-1])
     plt.title("|u(mu)|, |u_app(mu)|")
     plt.grid()
     plt.show()
 
     plt.figure(figsize = (15, 7))
     plt.jet()
     plt.semilogy(mu1, R)
     for l_ in approx.trainedModel.getPoles(mTMP):
         plt.plot([np.real(l_ ** exp)] * 2, [Rmin, Rmax], 'b:')
     plt.plot(mu2x, [Rmax] * len(mu2x), 'kx')
     plt.plot([murangeExp[0][0]] * 2, [Rmin, Rmax], 'm:')
     plt.plot([murangeExp[1][0]] * 2, [Rmin, Rmax], 'm:')
     plt.xlim(mu1[0], mu1[-1])
     if normalizeDen:
         plt.title("|Q(mu)res(mu)|")
     else:
         plt.title("|res(mu)|")
     plt.grid()
     plt.show()
     
     plt.figure(figsize = (15, 7))
     plt.jet()
     plt.semilogy(mu1, E)
     for l_ in approx.trainedModel.getPoles(mTMP):
         plt.plot([np.real(l_ ** exp)] * 2, [Emin, Emax], 'b:')
     plt.plot(mu2x, [Emax] * len(mu2x), 'kx')
     plt.plot([murangeExp[0][0]] * 2, [Emin, Emax], 'm:')
     plt.plot([murangeExp[1][0]] * 2, [Emin, Emax], 'm:')
     plt.xlim(mu1[0], mu1[-1])
     if normalizeDen:
         plt.title("|Q(mu)err(mu)|")
     else:
         plt.title("|err(mu)|")
     plt.grid()
     plt.show()
     
-    plt.figure(figsize = (15, 7))
-    plt.jet()
-    plt.plot(mu1, beta)
-    for l_ in approx.trainedModel.getPoles(mTMP):
-        plt.plot([np.real(l_ ** exp)] * 2, [betamin, betamax], 'b:')
-    plt.plot(mu2x, [betamax] * len(mu2x), 'kx')
-    plt.plot([murangeExp[0][0]] * 2, [betamin, betamax], 'm:')
-    plt.plot([murangeExp[1][0]] * 2, [betamin, betamax], 'm:')
-    plt.xlim(mu1[0], mu1[-1])
-    plt.title("beta(mu)")
-    plt.grid()
-    plt.show()
-    
-    plt.figure(figsize = (15, 7))
-    plt.jet()
-    plt.semilogy(mu1, R / beta)
-    plt.semilogy(mu1, E, '--')
-    for l_ in approx.trainedModel.getPoles(mTMP):
-        plt.plot([np.real(l_ ** exp)] * 2, [Emin, Emax], 'b:')
-    plt.plot(mu2x, [Emax] * len(mu2x), 'kx')
-    plt.plot([murangeExp[0][0]] * 2, [Emin, Emax], 'm:')
-    plt.plot([murangeExp[1][0]] * 2, [Emin, Emax], 'm:')
-    plt.xlim(mu1[0], mu1[-1])
-    if normalizeDen:
-        plt.title("|Q(mu)res(mu)/beta(mu)|")
-    else:
-        plt.title("|res(mu)/beta(mu)|")
-    plt.grid()
-    plt.show()
-    
-    plt.figure(figsize = (15, 7))
-    plt.jet()
-    plt.semilogy(mu1, eta)
-    for l_ in approx.trainedModel.getPoles(mTMP):
-        plt.plot([np.real(l_ ** exp)] * 2, [etamin, etamax], 'b:')
-    plt.plot(mu2x, [etamax] * len(mu2x), 'kx')
-    plt.plot([murangeExp[0][0]] * 2, [etamin, etamax], 'm:')
-    plt.plot([murangeExp[1][0]] * 2, [etamin, etamax], 'm:')
-    plt.xlim(mu1[0], mu1[-1])
-    plt.title("eta(mu)")
-    plt.grid()
-    plt.show()
+    if not np.isnan(beta[0]):
+        plt.figure(figsize = (15, 7))
+        plt.jet()
+        plt.plot(mu1, beta)
+        for l_ in approx.trainedModel.getPoles(mTMP):
+            plt.plot([np.real(l_ ** exp)] * 2, [betamin, betamax], 'b:')
+        plt.plot(mu2x, [betamax] * len(mu2x), 'kx')
+        plt.plot([murangeExp[0][0]] * 2, [betamin, betamax], 'm:')
+        plt.plot([murangeExp[1][0]] * 2, [betamin, betamax], 'm:')
+        plt.xlim(mu1[0], mu1[-1])
+        plt.title("beta(mu)")
+        plt.grid()
+        plt.show()
+        
+        plt.figure(figsize = (15, 7))
+        plt.jet()
+        plt.semilogy(mu1, R / beta)
+        plt.semilogy(mu1, E, '--')
+        for l_ in approx.trainedModel.getPoles(mTMP):
+            plt.plot([np.real(l_ ** exp)] * 2, [Emin, Emax], 'b:')
+        plt.plot(mu2x, [Emax] * len(mu2x), 'kx')
+        plt.plot([murangeExp[0][0]] * 2, [Emin, Emax], 'm:')
+        plt.plot([murangeExp[1][0]] * 2, [Emin, Emax], 'm:')
+        plt.xlim(mu1[0], mu1[-1])
+        if normalizeDen:
+            plt.title("|Q(mu)res(mu)/beta(mu)|")
+        else:
+            plt.title("|res(mu)/beta(mu)|")
+        plt.grid()
+        plt.show()
+        
+        plt.figure(figsize = (15, 7))
+        plt.jet()
+        plt.semilogy(mu1, eta)
+        for l_ in approx.trainedModel.getPoles(mTMP):
+            plt.plot([np.real(l_ ** exp)] * 2, [etamin, etamax], 'b:')
+        plt.plot(mu2x, [etamax] * len(mu2x), 'kx')
+        plt.plot([murangeExp[0][0]] * 2, [etamin, etamax], 'm:')
+        plt.plot([murangeExp[1][0]] * 2, [etamin, etamax], 'm:')
+        plt.xlim(mu1[0], mu1[-1])
+        plt.title("eta(mu)")
+        plt.grid()
+        plt.show()
 
 def plotInfSet32FromData(mNone, mus, Ze, Te, Re, Ee, beta, murange, approx,
                          mu0, exps = [2., 1.], clip = -1,
                          normalizeDen = False):
     if hasattr(approx, "mus"):
         mu2x = approx.mus(mNone[0]) ** exps[0]
         mu2y = approx.mus(mNone[1]) ** exps[1]
     else:
         mu2x, mu2y = mu0[mNone[0]] ** exps[0], mu0[mNone[1]] ** exps[1]
     murangeExp = [[murange[0][mNone[0]] ** exps[0],
                    murange[0][mNone[1]] ** exps[1]],
                   [murange[1][mNone[0]] ** exps[0],
                    murange[1][mNone[1]] ** exps[1]]]
     mu1s = np.unique([m[mNone[0]] for m in mus])
     mu2s = np.unique([m[mNone[1]] for m in mus])
     mu1 = np.power(mu1s, exps[0])
     mu2 = np.power(mu2s, exps[1])
     Mu1, Mu2 = np.meshgrid(np.real(mu1), np.real(mu2))
 
     Z = np.log10(Ze)
     T = np.log10(Te)
     R = np.log10(Re)
     E = np.log10(Ee)
     ZTmin, ZTmax = min(np.min(Z), np.min(T)), max(np.max(Z), np.max(T))
     Rmin, Rmax = np.min(R), np.max(R)
     Emin, Emax = np.min(E), np.max(E)
-    betamin, betamax = np.min(beta), np.max(beta)
+    if not np.isnan(beta[0, 0]):
+        betamin, betamax = np.min(beta), np.max(beta)
     if clip > 0:
         ZTmax -= clip * (ZTmax - ZTmin)
         cmap = plt.cm.bone
     else:
         cmap = plt.cm.jet
     warnings.simplefilter("ignore", category = UserWarning)
     plt.figure(figsize = (15, 7))
     plt.jet()
     p = plt.contourf(Mu1, Mu2, Z, cmap = cmap,
                      levels = np.linspace(ZTmin, ZTmax, 50))
     if clip > 0:
         plt.contour(Mu1, Mu2, Z, [ZTmin])
     plt.plot(mu2x, mu2y, 'kx')
     plt.plot([murangeExp[0][0]] * 2,
              [murangeExp[0][1], murangeExp[1][1]], 'm:')
     plt.plot([murangeExp[0][0], murangeExp[1][0]],
              [murangeExp[1][1]] * 2, 'm:')
     plt.plot([murangeExp[1][0]] * 2,
              [murangeExp[1][1], murangeExp[0][1]], 'm:')
     plt.plot([murangeExp[1][0], murangeExp[0][0]],
              [murangeExp[0][1]] * 2, 'm:')
     plt.colorbar(p)
     plt.title("log10|u(mu)|")
     plt.grid()
     plt.show()
 
     plt.figure(figsize = (15, 7))
     plt.jet()
     p = plt.contourf(Mu1, Mu2, T, cmap = cmap,
                      levels = np.linspace(ZTmin, ZTmax, 50))
     if clip > 0:
         plt.contour(Mu1, Mu2, T, [ZTmin])
     plt.plot(mu2x, mu2y, 'kx')
     plt.plot([murangeExp[0][0]] * 2,
              [murangeExp[0][1], murangeExp[1][1]], 'm:')
     plt.plot([murangeExp[0][0], murangeExp[1][0]],
              [murangeExp[1][1]] * 2, 'm:')
     plt.plot([murangeExp[1][0]] * 2,
              [murangeExp[1][1], murangeExp[0][1]], 'm:')
     plt.plot([murangeExp[1][0], murangeExp[0][0]],
              [murangeExp[0][1]] * 2, 'm:')
     plt.title("log10|u_app(mu)|")
     plt.colorbar(p)
     plt.grid()
     plt.show()
     
     plt.figure(figsize = (15, 7))
     plt.jet()
     p = plt.contourf(Mu1, Mu2, R, levels = np.linspace(Rmin, Rmax, 50))
     plt.plot(mu2x, mu2y, 'kx')
     plt.plot([murangeExp[0][0]] * 2,
              [murangeExp[0][1], murangeExp[1][1]], 'm:')
     plt.plot([murangeExp[0][0], murangeExp[1][0]],
              [murangeExp[1][1]] * 2, 'm:')
     plt.plot([murangeExp[1][0]] * 2,
              [murangeExp[1][1], murangeExp[0][1]], 'm:')
     plt.plot([murangeExp[1][0], murangeExp[0][0]],
              [murangeExp[0][1]] * 2, 'm:')
     if normalizeDen:
         plt.title("log10|Q(mu)res(mu)|")
     else:
         plt.title("log10|res(mu)|")
     plt.colorbar(p)
     plt.grid()
     plt.show()
     
     plt.figure(figsize = (15, 7))
     plt.jet()
     p = plt.contourf(Mu1, Mu2, E, levels = np.linspace(Emin, Emax, 50))
     plt.plot(mu2x, mu2y, 'kx')
     plt.plot([murangeExp[0][0]] * 2,
              [murangeExp[0][1], murangeExp[1][1]], 'm:')
     plt.plot([murangeExp[0][0], murangeExp[1][0]],
              [murangeExp[1][1]] * 2, 'm:')
     plt.plot([murangeExp[1][0]] * 2,
              [murangeExp[1][1], murangeExp[0][1]], 'm:')
     plt.plot([murangeExp[1][0], murangeExp[0][0]],
              [murangeExp[0][1]] * 2, 'm:')
     if normalizeDen:
         plt.title("log10|Q(mu)err(mu)|")
     else:
         plt.title("log10|err(mu)|")
     plt.colorbar(p)
     plt.grid()
     plt.show()
     
-    plt.figure(figsize = (15, 7))
-    plt.jet()
-    p = plt.contourf(Mu1, Mu2, beta,
-                     levels = np.linspace(betamin, betamax, 50))
-    plt.plot(mu2x, mu2y, 'kx')
-    plt.plot([murangeExp[0][0]] * 2,
-             [murangeExp[0][1], murangeExp[1][1]], 'm:')
-    plt.plot([murangeExp[0][0], murangeExp[1][0]],
-             [murangeExp[1][1]] * 2, 'm:')
-    plt.plot([murangeExp[1][0]] * 2,
-             [murangeExp[1][1], murangeExp[0][1]], 'm:')
-    plt.plot([murangeExp[1][0], murangeExp[0][0]],
-             [murangeExp[0][1]] * 2, 'm:')
-    plt.title("beta(mu)")
-    plt.colorbar(p)
-    plt.grid()
-    plt.show()
-    
-    plt.figure(figsize = (15, 7))
-    plt.jet()
-    p = plt.contourf(Mu1, Mu2, R - np.log10(beta),
-                     levels = np.linspace(Emin, Emax, 50))
-    plt.plot(mu2x, mu2y, 'kx')
-    plt.plot([murangeExp[0][0]] * 2,
-             [murangeExp[0][1], murangeExp[1][1]], 'm:')
-    plt.plot([murangeExp[0][0], murangeExp[1][0]],
-             [murangeExp[1][1]] * 2, 'm:')
-    plt.plot([murangeExp[1][0]] * 2,
-             [murangeExp[1][1], murangeExp[0][1]], 'm:')
-    plt.plot([murangeExp[1][0], murangeExp[0][0]],
-             [murangeExp[0][1]] * 2, 'm:')
-    if normalizeDen:
-        plt.title("log10|Q(mu)res(mu)/beta(mu)|")
-    else:
-        plt.title("log10|res(mu)/beta(mu)|")
-    plt.colorbar(p)
-    plt.grid()
-    plt.show()
-    
-    plt.figure(figsize = (15, 7))
-    plt.jet()
-    p = plt.contourf(Mu1, Mu2, R - np.log10(beta) - E, 50)
-    plt.plot(mu2x, mu2y, 'kx')
-    plt.plot([murangeExp[0][0]] * 2,
-             [murangeExp[0][1], murangeExp[1][1]], 'm:')
-    plt.plot([murangeExp[0][0], murangeExp[1][0]],
-             [murangeExp[1][1]] * 2, 'm:')
-    plt.plot([murangeExp[1][0]] * 2,
-             [murangeExp[1][1], murangeExp[0][1]], 'm:')
-    plt.plot([murangeExp[1][0], murangeExp[0][0]],
-             [murangeExp[0][1]] * 2, 'm:')
-    plt.title("log10|eta(mu)|")
-    plt.colorbar(p)
-    plt.grid()
-    plt.show()
+    if not np.isnan(beta[0, 0]):
+        plt.figure(figsize = (15, 7))
+        plt.jet()
+        p = plt.contourf(Mu1, Mu2, beta,
+                         levels = np.linspace(betamin, betamax, 50))
+        plt.plot(mu2x, mu2y, 'kx')
+        plt.plot([murangeExp[0][0]] * 2,
+                 [murangeExp[0][1], murangeExp[1][1]], 'm:')
+        plt.plot([murangeExp[0][0], murangeExp[1][0]],
+                 [murangeExp[1][1]] * 2, 'm:')
+        plt.plot([murangeExp[1][0]] * 2,
+                 [murangeExp[1][1], murangeExp[0][1]], 'm:')
+        plt.plot([murangeExp[1][0], murangeExp[0][0]],
+                 [murangeExp[0][1]] * 2, 'm:')
+        plt.title("beta(mu)")
+        plt.colorbar(p)
+        plt.grid()
+        plt.show()
+        
+        plt.figure(figsize = (15, 7))
+        plt.jet()
+        p = plt.contourf(Mu1, Mu2, R - np.log10(beta),
+                         levels = np.linspace(Emin, Emax, 50))
+        plt.plot(mu2x, mu2y, 'kx')
+        plt.plot([murangeExp[0][0]] * 2,
+                 [murangeExp[0][1], murangeExp[1][1]], 'm:')
+        plt.plot([murangeExp[0][0], murangeExp[1][0]],
+                 [murangeExp[1][1]] * 2, 'm:')
+        plt.plot([murangeExp[1][0]] * 2,
+                 [murangeExp[1][1], murangeExp[0][1]], 'm:')
+        plt.plot([murangeExp[1][0], murangeExp[0][0]],
+                 [murangeExp[0][1]] * 2, 'm:')
+        if normalizeDen:
+            plt.title("log10|Q(mu)res(mu)/beta(mu)|")
+        else:
+            plt.title("log10|res(mu)/beta(mu)|")
+        plt.colorbar(p)
+        plt.grid()
+        plt.show()
+        
+        plt.figure(figsize = (15, 7))
+        plt.jet()
+        p = plt.contourf(Mu1, Mu2, R - np.log10(beta) - E, 50)
+        plt.plot(mu2x, mu2y, 'kx')
+        plt.plot([murangeExp[0][0]] * 2,
+                 [murangeExp[0][1], murangeExp[1][1]], 'm:')
+        plt.plot([murangeExp[0][0], murangeExp[1][0]],
+                 [murangeExp[1][1]] * 2, 'm:')
+        plt.plot([murangeExp[1][0]] * 2,
+                 [murangeExp[1][1], murangeExp[0][1]], 'm:')
+        plt.plot([murangeExp[1][0], murangeExp[0][0]],
+                 [murangeExp[0][1]] * 2, 'm:')
+        plt.title("log10|eta(mu)|")
+        plt.colorbar(p)
+        plt.grid()
+        plt.show()
     
 def plotInfSet3(murange, murangeEff, approx, mu0, nSamples = 200,
                 marginal = [None, None, .05], exps = [2., 1., 1.], clip = -1,
-                relative = True, normalizeDen = False):
+                relative = True, normalizeDen = False, nobeta = False):
     mNone = [i for i, m in enumerate(marginal) if m is None]
     nNone = len(mNone)
     if nNone not in [1, 2]:
         raise
     if nNone == 1:
         exp = exps[mNone[0]]
         mu1 = np.linspace(murangeEff[0][mNone[0]] ** exp,
                           murangeEff[1][mNone[0]] ** exp, nSamples)
         mu1s = np.power(mu1, 1. / exp)
         mus = []
         mBase = copy(marginal)
         for m1 in mu1s:
             mBase[mNone[0]] = m1
             mus += [copy(mBase)]
         Z = approx.normHF(mus)
         T = approx.normApprox(mus)
         R = approx.normRes(mus)
         E = approx.normErr(mus)
         if relative:
             F = approx.normRHS(mus)
             R /= F
             E /= Z
         if normalizeDen:
             Qvals = np.abs(approx.trainedModel.getQVal(mus))
             R *= Qvals
             E *= Qvals
-        beta = approx.HFEngine.stabilityFactor(approx.getHF(mus), mus)
+        if nobeta:
+            beta = np.empty(len(mus))
+            beta[:] = np.nan
+        else:
+            beta = approx.HFEngine.stabilityFactor(approx.getHF(mus), mus)
         plotInfSet31FromData(mNone, mus, Z, T, R, E, beta, murange, approx,
                              mu0, exp, normalizeDen)
         return mus, Z, T, R, E, beta
     if nNone == 2:
         exps = [exps[m] for m in mNone]
         mu1 = np.linspace(murangeEff[0][mNone[0]] ** exps[0],
                           murangeEff[1][mNone[0]] ** exps[0], nSamples)
         mu2 = np.linspace(murangeEff[0][mNone[1]] ** exps[1],
                           murangeEff[1][mNone[1]] ** exps[1], nSamples)
         mu1s = np.power(mu1, 1. / exps[0])
         mu2s = np.power(mu2, 1. / exps[1])
         Mu1, Mu2 = np.meshgrid(np.real(mu1), np.real(mu2))
         mus = []
         mBase = copy(marginal)
         for m2 in mu2s:
             mBase[mNone[1]] = m2
             for m1 in mu1s:
                 mBase[mNone[0]] = m1
                 mus += [copy(mBase)]
         Ze = approx.normHF(mus).reshape((nSamples, nSamples))
         Te = approx.normApprox(mus).reshape((nSamples, nSamples))
         Re = approx.normRes(mus).reshape((nSamples, nSamples))
         Ee = approx.normErr(mus).reshape((nSamples, nSamples))
         if relative:
             Fe = approx.normRHS(mus).reshape((nSamples, nSamples))
             Re /= Fe
             Ee /= Ze
         if normalizeDen:
             Qvals = np.abs(approx.trainedModel.getQVal(mus).reshape(
                                                          (nSamples, nSamples)))
             Re *= Qvals
             Ee *= Qvals
-        betae = approx.HFEngine.stabilityFactor(approx.getHF(mus), mus)\
+        if nobeta:
+            betae = np.empty((nSamples, nSamples))
+            betae[:, :] = np.nan
+        else:
+            betae = approx.HFEngine.stabilityFactor(approx.getHF(mus), mus)\
                                                  .reshape((nSamples, nSamples))
         plotInfSet32FromData(mNone, mus, Ze, Te, Re, Ee, betae, murange,
                              approx, mu0, exps, clip, normalizeDen)
         return mus, Ze, Te, Re, Ee, betae
 
diff --git a/rrompy/reduction_methods/distributed/rb_distributed.py b/rrompy/reduction_methods/distributed/rb_distributed.py
index 478569a..aa64217 100644
--- a/rrompy/reduction_methods/distributed/rb_distributed.py
+++ b/rrompy/reduction_methods/distributed/rb_distributed.py
@@ -1,212 +1,210 @@
 # Copyright (C) 2018 by the RROMPy authors
 #
 # This file is part of RROMPy.
 #
 # RROMPy is free software: you can redistribute it and/or modify
 # it under the terms of the GNU Lesser General Public License as published by
 # the Free Software Foundation, either version 3 of the License, or
 # (at your option) any later version.
 #
 # RROMPy is distributed in the hope that it will be useful,
 # but WITHOUT ANY WARRANTY; without even the implied warranty of
 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
 # GNU Lesser General Public License for more details.
 #
 # You should have received a copy of the GNU Lesser General Public License
 # along with RROMPy. If not, see <http://www.gnu.org/licenses/>.
 #
 
 from copy import deepcopy as copy
 import numpy as np
 from .generic_distributed_approximant import GenericDistributedApproximant
 from rrompy.reduction_methods.trained_model import TrainedModelRB as tModel
 from rrompy.reduction_methods.trained_model import TrainedModelData
 from rrompy.reduction_methods.base.rb_utils import projectAffineDecomposition
 from rrompy.utilities.base.types import (Np1D, Np2D, List, Tuple, DictAny,
                                          HFEng, paramVal, sampList)
 from rrompy.utilities.base import verbosityManager as vbMng
 from rrompy.utilities.exception_manager import (RROMPyWarning, RROMPyException,
                                                 RROMPyAssert)
 
 __all__ = ['RBDistributed']
 
 class RBDistributed(GenericDistributedApproximant):
     """
     ROM RB approximant computation for parametric problems.
 
     Args:
         HFEngine: HF problem solver.
         mu0(optional): Default parameter. Defaults to 0.
         approxParameters(optional): Dictionary containing values for main
             parameters of approximant. Recognized keys are:
             - 'POD': whether to compute POD of snapshots; defaults to True;
             - 'S': total number of samples current approximant relies upon;
             - 'sampler': sample point generator;
             - 'R': rank for Galerkin projection; defaults to prod(S);
             - 'PODTolerance': tolerance for snapshots POD; defaults to -1.
             Defaults to empty dict.
         homogeneized(optional): Whether to homogeneize Dirichlet BCs. Defaults
             to False.
         verbosity(optional): Verbosity level. Defaults to 10.
             
     Attributes:
         HFEngine: HF problem solver.
         mu0: Default parameter.
         mus: Array of snapshot parameters.
         homogeneized: Whether to homogeneize Dirichlet BCs.
         approxRadius: Dummy radius of approximant (i.e. distance from mu0 to
             farthest sample point).
         approxParameters: Dictionary containing values for main parameters of
             approximant. Recognized keys are in parameterList.
         parameterListSoft: Recognized keys of soft approximant parameters:
             - 'POD': whether to compute POD of snapshots.
             - 'R': rank for Galerkin projection;
             - 'PODTolerance': tolerance for snapshots POD.
         parameterListCritical: Recognized keys of critical approximant
             parameters:
             - 'S': total number of samples current approximant relies upon;
             - 'sampler': sample point generator;
         POD: Whether to compute POD of snapshots.
         S: Number of solution snapshots over which current approximant is
             based upon.
         sampler: Sample point generator.
         R: Rank for Galerkin projection.
         muBounds: list of bounds for 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.
         As: List of sparse matrices (in CSC format) representing coefficients
             of linear system matrix wrt theta(mu).
         bs: List of numpy vectors representing coefficients of linear system
            RHS wrt theta(mu).
         thetaAs: List of callables representing coefficients of linear system
             matrix wrt mu.
         thetabs: List of callables representing coefficients of linear system
             RHS wrt mu.
         ARBs: List of sparse matrices (in CSC format) representing coefficients
             of compressed linear system matrix wrt theta(mu).
         bRBs: List of numpy vectors representing coefficients of compressed 
             linear system RHS wrt theta(mu).
     """
     
     def __init__(self, HFEngine:HFEng, mu0 : paramVal = None,
                  approxParameters : DictAny = {}, homogeneized : bool = False,
                  verbosity : int = 10, timestamp : bool = True):
         self._preInit()
         self._addParametersToList(["R", "PODTolerance"], [1, -1])
         super().__init__(HFEngine = HFEngine, mu0 = mu0,
                          approxParameters = approxParameters,
                          homogeneized = homogeneized,
                          verbosity = verbosity, timestamp = timestamp)
         vbMng(self, "INIT", "Computing affine blocks of system.", 10)
         self.As = self.HFEngine.affineLinearSystemA(self.mu0)
         self.bs = self.HFEngine.affineLinearSystemb(self.mu0,
                                                     self.homogeneized)
         vbMng(self, "DEL", "Done computing affine blocks.", 10)
         self._postInit()
 
     @property
     def S(self):
         """Value of S."""
         return self._S
     @S.setter
     def S(self, S):
         GenericDistributedApproximant.S.fset(self, S)
         if not hasattr(self, "_R"): self._R = np.prod(self.S)
 
     @property
     def R(self):
         """Value of R. Its assignment may change S."""
         return self._R
     @R.setter
     def R(self, R):
         if R < 0: raise RROMPyException("R must be non-negative.")
         self._R = R
         self._approxParameters["R"] = self.R
         if hasattr(self, "_S") and np.prod(self.S) < self.R:
             RROMPyWarning("Prescribed S is too small. Decreasing R.")
             self.R = np.prod(self.S)
 
     @property
     def PODTolerance(self):
         """Value of PODTolerance."""
         return self._PODTolerance
     @PODTolerance.setter
     def PODTolerance(self, PODTolerance):
         self._PODTolerance = PODTolerance
         self._approxParameters["PODTolerance"] = self.PODTolerance
 
     def setupApprox(self):
         """Compute RB projection matrix."""
         if self.checkComputedApprox():
             return
         RROMPyAssert(self._mode, message = "Cannot setup approximant.")
         vbMng(self, "INIT", "Setting up {}.". format(self.name()), 5)
         self.computeSnapshots()
         vbMng(self, "INIT", "Computing projection matrix.", 7)
         if self.POD:
             U, s, _ = np.linalg.svd(self.samplingEngine.RPOD)
             s = s ** 2.
         else:
             Gramian = self.HFEngine.innerProduct(self.samplingEngine.samples,
                                                  self.samplingEngine.samples)
             U, s, _ = np.linalg.svd(Gramian)
         nsamples = self.samplingEngine.nsamples
         snorm = np.cumsum(s[::-1]) / np.sum(s)
         nPODTrunc = min(nsamples - np.argmax(snorm > self.PODTolerance),
                         self.R)
         pMat = self.samplingEngine.samples.dot(U[:, : nPODTrunc])
-        from rrompy.sampling import sampleList
-        pMat = sampleList(self.samplingEngine.samples.data[:, : nPODTrunc])
         vbMng(self, "MAIN",
               ("Assembling {}x{} projection matrix from {} "
                "samples.").format(*(pMat.shape), nsamples), 5)
         if self.trainedModel is None:
             self.trainedModel = tModel()
             self.trainedModel.verbosity = self.verbosity
             self.trainedModel.timestamp = self.timestamp
             data = TrainedModelData(self.trainedModel.name(), self.mu0,
                                     pMat, self.HFEngine.rescalingExp)
             data.thetaAs = self.HFEngine.affineWeightsA(self.mu0)
             data.thetabs = self.HFEngine.affineWeightsb(self.mu0,
                                                         self.homogeneized)
             data.mus = copy(self.mus)
             self.trainedModel.data = data
         else:
             self.trainedModel = self.trainedModel
             self.trainedModel.data.projMat = copy(pMat)
         ARBs, bRBs = self.assembleReducedSystem(pMat)
         self.trainedModel.data.ARBs = ARBs
         self.trainedModel.data.bRBs = bRBs
 
         vbMng(self, "DEL", "Done computing projection matrix.", 7)
         self.trainedModel.data.approxParameters = copy(self.approxParameters)
         vbMng(self, "DEL", "Done setting up approximant.", 5)
 
     def assembleReducedSystem(self, pMat : sampList = None,
                               pMatOld : sampList = None)\
                                               -> Tuple[List[Np2D], List[Np1D]]:
         """Build affine blocks of RB linear system through projections."""
         if pMat is None:
             self.setupApprox()
             ARBs = self.trainedModel.data.ARBs
             bRBs = self.trainedModel.data.bRBs
         else:
             vbMng(self, "INIT", "Projecting affine terms of HF model.", 10)
             ARBsOld = None if pMatOld is None else self.trainedModel.data.ARBs
             bRBsOld = None if pMatOld is None else self.trainedModel.data.bRBs
             ARBs, bRBs = projectAffineDecomposition(self.As, self.bs, pMat,
                                                     ARBsOld, bRBsOld, pMatOld)
             vbMng(self, "DEL", "Done projecting affine terms.", 10)
         return ARBs, bRBs