diff --git a/examples/3d/pod/fracture3_pod.py b/examples/3d/pod/fracture3_pod.py
index d1eb22f..b1b9a60 100644
--- a/examples/3d/pod/fracture3_pod.py
+++ b/examples/3d/pod/fracture3_pod.py
@@ -1,195 +1,195 @@
 import numpy as np
 from mpl_toolkits.mplot3d import Axes3D
 from matplotlib import pyplot as plt
 from rrompy.hfengines.linear_problem.tridimensional import \
                                                  MembraneFractureEngine3 as MFE
 from rrompy.reduction_methods.distributed import RationalInterpolant as RI
 from rrompy.reduction_methods.distributed import RBDistributed as RBD
 from rrompy.parameter.parameter_sampling import (QuadratureSampler as QS,
                                                  QuadratureSamplerTotal as QST,
                                                  ManualSampler as MS,
                                                  RandomSampler as RS)
 
 verb = 50
 size = 4
 show_sample = False
 show_norm = False
 clip = -1
 #clip = .4
 #clip = .6
 homogeneize = False
 #homogeneize = True
 
 Delta = 0
-MN = 10
+MN = 5
 R = (MN + 3) * (MN + 2) * (MN + 1) // 6
 S = [int(np.ceil(R ** (1. / 3.)))] * 3
 PODTol = 1e-8
 
 samples = "centered"
 samples = "distributed"
 algo = "rational"
 #algo = "RB"
 sampling = "quadrature"
 sampling = "quadrature_total"
 #sampling = "random"
 
 if samples == "distributed":
     radial = 0
-#    radial = "gaussian"
+    radial = "gaussian"
 #    radial = "thinplate"
 #    radial = "multiquadric"
     rW0 = 5.
-    radialWeight = [rW0] * 2
+    radialWeight = [rW0] * 3
 
 assert Delta <= 0
 
 if size == 1:
     mu0 = [47.5 ** .5, .4, .05]
     mutar = [50 ** .5, .45, .07]
     murange = [[40 ** .5, .3, .025], [55 ** .5, .5, .075]]
 if size == 2:
     mu0 = [50 ** .5, .3, .02]
     mutar = [55 ** .5, .4, .03]
     murange = [[40 ** .5, .1, .005], [60 ** .5, .5, .035]]
 if size == 3:
     mu0 = [45 ** .5, .5, .05]
     mutar = [47 ** .5, .4, .045]
     murange = [[40 ** .5, .3, .04], [50 ** .5, .7, .06]]
 if size == 4:
     mu0 = [45 ** .5, .4, .05]
     mutar = [47 ** .5, .45, .045]
     murange = [[40 ** .5, .3, .04], [50 ** .5, .5, .06]]
 
 aEff = 1.#25
 bEff = 1. - aEff
 murangeEff = [[(aEff*murange[0][0]**2. + bEff*murange[1][0]**2.) ** .5,
                 aEff*murange[0][1] + bEff*murange[1][1],
                 aEff*murange[0][2] + bEff*murange[1][2]],
               [(aEff*murange[1][0]**2. + bEff*murange[0][0]**2.) ** .5,
                 aEff*murange[1][1] + bEff*murange[0][1],
                 aEff*murange[1][2] + bEff*murange[0][2]]]
 
 H = 1.
 L = .75
 delta = .05
 n = 50
 
 solver = MFE(mu0 = mu0, H = H, L = L, delta = delta, n = n, verbosity = verb)
 
 rescaling = [lambda x: np.power(x, 2.), lambda x: x, lambda x: x]
 rescalingInv = [lambda x: np.power(x, .5), lambda x: x, lambda x: x]
 if algo == "rational":
     params = {'N':MN, 'M':MN + Delta, 'S':S, 'POD':True}
     if samples == "distributed":
         params['polybasis'] = "CHEBYSHEV"
 #            params['polybasis'] = "LEGENDRE"
 #            params['polybasis'] = "MONOMIAL"
         params['E'] = MN
         params['radialBasis'] = radial
         params['radialBasisWeights'] = radialWeight
         method = RI
     elif samples == "centered":
         params['polybasis'] = "MONOMIAL"
         params['S'] = R
         method = RI
 else: #if algo == "RB":
     params = {'R':(MN + 3 + Delta) * (MN + 2 + Delta) * (MN + 1 + Delta) // 6,
               'S':S, 'POD':True, 'PODTolerance':PODTol}
     if samples == "distributed":
         method = RBD
     elif samples == "centered":
         params['S'] = R
         method = RBD
 
 if samples == "distributed":
     if sampling == "quadrature":
         params['sampler'] = QS(murange, "CHEBYSHEV", scaling = rescaling,
                                scalingInv = rescalingInv)
 #        params['sampler'] = QS(murange, "GAUSSLEGENDRE", scaling = rescaling,
 #                               scalingInv = rescalingInv)
 #        params['sampler'] = QS(murange, "UNIFORM", scaling = rescaling,
 #                               scalingInv = rescalingInv)
         params['S'] = [max(j, MN + 1) for j in params['S']]
     elif sampling == "quadrature_total":
         params['sampler'] = QST(murange, "CHEBYSHEV", scaling = rescaling,
                                 scalingInv = rescalingInv)
         params['S'] = R
     else: # if sampling == "random":
         params['sampler'] = RS(murange, "HALTON", scaling = rescaling,
                                scalingInv = rescalingInv)
         params['S'] = R
         
 elif samples == "centered":
     params['sampler'] = MS(murange, points = [mu0], scaling = rescaling,
                            scalingInv = rescalingInv)
 
 approx = method(solver, mu0 = mu0, approxParameters = params,
                 verbosity = verb, homogeneized = homogeneize)
 if samples == "distributed": approx.samplingEngine.allowRepeatedSamples = False
 
 approx.setupApprox()
 if show_sample:
     from fracture3_warping import fracture3_warping
     warps = fracture3_warping(solver.V.mesh(), L, mutar[1], delta, mutar[2])
     approx.plotApprox(mutar, warps, name = 'u_app',
                       homogeneized = False, what = "REAL")
     approx.plotHF(mutar, warps, name = 'u_HF',
                   homogeneized = False, what = "REAL")
     approx.plotErr(mutar, warps, name = 'err',
                    homogeneized = False, what = "REAL")
 #    approx.plotRes(mutar, warps, name = 'res',
 #                   homogeneized = False, what = "REAL")
     appErr = approx.normErr(mutar, homogeneized = homogeneize)
     solNorm = approx.normHF(mutar, homogeneized = homogeneize)
     resNorm = approx.normRes(mutar, homogeneized = homogeneize)
     RHSNorm = approx.normRHS(mutar, homogeneized = homogeneize)
     print(('SolNorm:\t{}\nErr:\t{}\nErrRel:\t{}').format(solNorm, appErr,
                                                np.divide(appErr, solNorm)))
     print(('RHSNorm:\t{}\nRes:\t{}\nResRel:\t{}').format(RHSNorm, resNorm,
                                               np.divide(resNorm, RHSNorm)))
 
 fig = plt.figure(figsize = (8, 6))
 ax = Axes3D(fig)
 ax.scatter(approx.trainedModel.data.mus(0) ** 2.,
            approx.trainedModel.data.mus(1),
            approx.trainedModel.data.mus(2), '.')
 plt.show()
 plt.close()
 
 approx.verbosity = 0
 approx.trainedModel.verbosity = 0
 if algo == "rational" and approx.N > 0:
     from plot_zero_set_3 import plotZeroSet3
 #    muZeroVals, Qvals = plotZeroSet3(murange, murangeEff, approx, mu0, 200,
 #                                     [None, mu0[1], mu0[2]], [2., 1., 1.],
 #                                     clip = clip)
 #    muZeroVals, Qvals = plotZeroSet3(murange, murangeEff, approx, mu0, 200,
 #                                     [None, None, mu0[2]], [2., 1., 1.],
 #                                     clip = clip)
 #    muZeroVals, Qvals = plotZeroSet3(murange, murangeEff, approx, mu0, 200,
 #                                     [None, mu0[1], None], [2., 1., 1.],
 #                                     clip = clip)
     muZeroScatter = plotZeroSet3(murange, murangeEff, approx, mu0, 50,
                                  [None, None, None], [2., 1., 1.], clip = clip)
 
 if show_norm:
     solver._solveBatchSize = 25
     from plot_inf_set_3 import plotInfSet3
     muInfVals, normEx, normApp, normRes, normErr, beta = plotInfSet3(
                                        murange, murangeEff, approx, mu0, 25,
                                        [None, mu0[1], mu0[2]], [2., 1., 1.],
                                        clip = clip, relative = False,
                                        normalizeDen = True)
     muInfVals, normEx, normApp, normRes, normErr, beta = plotInfSet3(
                                        murange, murangeEff, approx, mu0, 25,
                                        [None, None, mu0[2]], [2., 1., 1.],
                                        clip = clip, relative = False,
                                        normalizeDen = True)
     muInfVals, normEx, normApp, normRes, normErr, beta = plotInfSet3(
                                        murange, murangeEff, approx, mu0, 25,
                                        [None, mu0[1], None], [2., 1., 1.],
                                        clip = clip, relative = False,
                                        normalizeDen = True)
 
 print(approx.getPoles([None, .5, .05]) ** 2.)
diff --git a/rrompy/reduction_methods/distributed/rational_interpolant.py b/rrompy/reduction_methods/distributed/rational_interpolant.py
index 3845d8e..d8ce41e 100644
--- a/rrompy/reduction_methods/distributed/rational_interpolant.py
+++ b/rrompy/reduction_methods/distributed/rational_interpolant.py
@@ -1,571 +1,548 @@
 # Copyright (C) 2018 by the RROMPy authors
 #
 # This file is part of RROMPy.
 #
 # RROMPy is free software: you can redistribute it and/or modify
 # it under the terms of the GNU Lesser General Public License as published by
 # the Free Software Foundation, either version 3 of the License, or
 # (at your option) any later version.
 #
 # RROMPy is distributed in the hope that it will be useful,
 # but WITHOUT ANY WARRANTY; without even the implied warranty of
 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
 # GNU Lesser General Public License for more details.
 #
 # You should have received a copy of the GNU Lesser General Public License
 # along with RROMPy. If not, see <http://www.gnu.org/licenses/>.
 #
 
 from copy import deepcopy as copy
 import numpy as np
 from rrompy.reduction_methods.base import checkRobustTolerance
 from .generic_distributed_approximant import GenericDistributedApproximant
-from rrompy.utilities.poly_fitting import customFit, customPInv
+from rrompy.utilities.poly_fitting import customPInv
 from rrompy.utilities.poly_fitting.polynomial import (polybases, polyfitname,
                                                 nextDerivativeIndices,
                                                 hashDerivativeToIdx as hashD,
                                                 hashIdxToDerivative as hashI,
                                                 homogeneizedpolyvander as hpvP,
-                                                homogeneizedToFull)
+                                                homogeneizedToFull,
+                                                PolynomialInterpolator as PI)
 from rrompy.utilities.poly_fitting.radial_basis import (rbbases,
-                                                radialFunction,
-                                                polyfitname as polyfitnameR,
-                                                homogeneizedpolyvander as hpvR)
+                                                RadialBasisInterpolator as RBI)
 from rrompy.reduction_methods.trained_model import TrainedModelPade as tModel
 from rrompy.reduction_methods.trained_model import TrainedModelData
 from rrompy.utilities.base.types import (Np1D, Np2D, HFEng, DictAny, Tuple,
                                          List, paramVal, paramList, sampList)
 from rrompy.utilities.base import verbosityManager as vbMng, multifactorial
 from rrompy.utilities.exception_manager import (RROMPyException, RROMPyAssert,
                                                 RROMPyWarning)
 
 __all__ = ['RationalInterpolant']
 
 class RationalInterpolant(GenericDistributedApproximant):
     """
     ROM rational interpolant 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;
             - 'polybasis': type of polynomial basis for interpolation; allowed
                 values include 'MONOMIAL', 'CHEBYSHEV' and 'LEGENDRE'; defaults
                 to 'MONOMIAL';
             - 'E': number of derivatives used to compute interpolant; defaults
                 to 0;
             - 'M': degree of Pade' interpolant numerator; defaults to 0;
             - 'N': degree of Pade' interpolant denominator; defaults to 0;
             - 'radialBasis': radial basis family for interpolant numerator;
                 defaults to 0, i.e. no radial basis;
             - 'radialBasisWeights': radial basis weights for interpolant
                 numerator; defaults to 0, i.e. identity;
             - 'interpRcond': tolerance for interpolation; defaults to None;
             - 'robustTol': tolerance for robust Pade' denominator management; 
                 defaults to 0.
             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.
         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;
             - 'polybasis': type of polynomial basis for interpolation;
             - 'E': number of derivatives used to compute interpolant;
             - 'M': degree of Pade' interpolant numerator;
             - 'N': degree of Pade' interpolant denominator;
             - 'radialBasis': radial basis family for interpolant numerator;
             - 'radialBasisWeights': radial basis weights for interpolant
                 numerator;
             - 'interpRcond': tolerance for interpolation via numpy.polyfit;
             - 'robustTol': tolerance for robust Pade' denominator management.
         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.
         polybasis: type of polynomial basis for interpolation.
         M: Numerator degree of approximant.
         N: Denominator degree of approximant.
         radialBasis: Radial basis family for interpolant numerator.
         radialBasisWeights: Radial basis weights for interpolant numerator.
         interpRcond: Tolerance for interpolation via numpy.polyfit.
         robustTol: Tolerance for robust Pade' denominator management.
         E: Complete derivative depth level of S.
         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.
         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 __init__(self, HFEngine:HFEng, mu0 : paramVal = None,
                  approxParameters : DictAny = {}, homogeneized : bool = False,
                  verbosity : int = 10, timestamp : bool = True):
         self._preInit()
         self._addParametersToList(["polybasis", "E", "M", "N", "radialBasis",
                                    "radialBasisWeights", "interpRcond",
                                    "robustTol"],
                                   ["MONOMIAL", -1, 0, 0, 0, 1, -1, 0])
         super().__init__(HFEngine = HFEngine, mu0 = mu0,
                          approxParameters = approxParameters,
                          homogeneized = homogeneized,
                          verbosity = verbosity, timestamp = timestamp)
         self._postInit()
 
     @property
     def polybasis(self):
         """Value of polybasis."""
         return self._polybasis
     @polybasis.setter
     def polybasis(self, polybasis):
         try:
             polybasis = polybasis.upper().strip().replace(" ","")
             if polybasis not in polybases: 
                 raise RROMPyException("Prescribed polybasis not recognized.")
             self._polybasis = polybasis
         except:
             RROMPyWarning(("Prescribed polybasis not recognized. Overriding "
                            "to 'MONOMIAL'."))
             self._polybasis = "MONOMIAL"
         self._approxParameters["polybasis"] = self.polybasis
 
     @property
     def radialBasis(self):
         """Value of radialBasis."""
         return self._radialBasis
     @radialBasis.setter
     def radialBasis(self, radialBasis):
         try:
             if radialBasis != 0:
                 radialBasis = radialBasis.upper().strip().replace(" ","")
                 if radialBasis not in rbbases: 
                     raise RROMPyException(("Prescribed radialBasis not "
                                            "recognized."))
             self._radialBasis = radialBasis
         except:
             RROMPyWarning(("Prescribed radialBasis not recognized. Overriding "
                            "to 0."))
             self._radialBasis = 0
         self._approxParameters["radialBasis"] = self.radialBasis
 
     @property
     def polybasisP(self):
         if self.radialBasis == 0:
             return self._polybasis
         return self._polybasis + "_" + self.radialBasis
 
     @property
     def interpRcond(self):
         """Value of interpRcond."""
         return self._interpRcond
     @interpRcond.setter
     def interpRcond(self, interpRcond):
         self._interpRcond = interpRcond
         self._approxParameters["interpRcond"] = self.interpRcond
 
     @property
     def radialBasisWeights(self):
         """Value of radialBasisWeights."""
         return self._radialBasisWeights
     @radialBasisWeights.setter
     def radialBasisWeights(self, radialBasisWeights):
         self._radialBasisWeights = radialBasisWeights
         self._approxParameters["radialBasisWeights"] = self.radialBasisWeights
 
     @property
     def M(self):
         """Value of M. Its assignment may change S."""
         return self._M
     @M.setter
     def M(self, M):
         if M < 0: raise RROMPyException("M must be non-negative.")
         self._M = M
         self._approxParameters["M"] = self.M
         if hasattr(self, "_E") and self.E >= 0 and self.E < self.M:
             RROMPyWarning("Prescribed S is too small. Decreasing M.")
             self.M = self.E
 
     @property
     def N(self):
         """Value of N. Its assignment may change S."""
         return self._N
     @N.setter
     def N(self, N):
         if N < 0: raise RROMPyException("N must be non-negative.")
         self._N = N
         self._approxParameters["N"] = self.N
         if hasattr(self, "_E") and self.E >= 0 and self.E < self.N:
             RROMPyWarning("Prescribed S is too small. Decreasing N.")
             self.N = self.E
 
     @property
     def E(self):
         """Value of E."""
         return self._E
     @E.setter
     def E(self, E):
         if E < 0:
             if not hasattr(self, "_S"):
                 raise RROMPyException(("Value of E must be positive if S is "
                                        "not yet assigned."))
             E = np.sum(hashI(np.prod(self.S), self.npar)) - 1
         self._E = E
         self._approxParameters["E"] = self.E
         if (hasattr(self, "_S")
         and self.E >= np.sum(hashI(np.prod(self.S), self.npar))):
             RROMPyWarning("Prescribed S is too small. Decreasing E.")
             self.E = -1
         if hasattr(self, "_M"): self.M = self.M
         if hasattr(self, "_N"): self.N = self.N
 
     @property
     def robustTol(self):
         """Value of tolerance for robust Pade' denominator management."""
         return self._robustTol
     @robustTol.setter
     def robustTol(self, robustTol):
         if robustTol < 0.:
             RROMPyWarning(("Overriding prescribed negative robustness "
                            "tolerance to 0."))
             robustTol = 0.
         self._robustTol = robustTol
         self._approxParameters["robustTol"] = self.robustTol
 
     @property
     def S(self):
         """Value of S."""
         return self._S
     @S.setter
     def S(self, S):
         GenericDistributedApproximant.S.fset(self, S)
         if hasattr(self, "_M"): self.M = self.M
         if hasattr(self, "_N"): self.N = self.N
         if hasattr(self, "_E"): self.E = self.E
 
     def resetSamples(self):
         """Reset samples."""
         super().resetSamples()
         self._musUniqueCN = None
         self._derIdxs = None
         self._reorder = None
 
     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.mus):
             self._musUniqueCN, musIdxsTo, musIdxs, musCount = (
                    self.centerNormalize(self.mus).unique(return_index = True, 
                                                          return_inverse = True,
                                                          return_counts = True))
             self._musUnique = self.mus[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.mus.shape[1],
                                                          cnt)
                 jIdx = np.nonzero(musIdxs == j)[0]
                 self._reorder[jIdx] = np.arange(filled, filled + cnt)
                 filled += cnt
 
     def _setupDenominator(self):
         """Compute Pade' denominator."""
         RROMPyAssert(self._mode, message = "Cannot setup denominator.")
         vbMng(self, "INIT", "Starting computation of denominator.", 7)
         while self.N > 0:
             invD = self._computeInterpolantInverseBlocks()
             if self.POD:
                 ev, eV = self.findeveVGQR(self.samplingEngine.RPOD, invD)
             else:
                 ev, eV = self.findeveVGExplicit(self.samplingEngine.samples,
                                                 invD)
             newParams = checkRobustTolerance(ev, self.N, self.robustTol)
             if not newParams:
                 break
             self.approxParameters = newParams
         if self.N <= 0:
             self._N = 0
             eV = np.ones((1, 1))
-        q = homogeneizedToFull(tuple([self.N + 1] * self.npar), self.npar, 
-                               eV[:, 0])
+        q = PI()
+        q.npar = self.npar
+        q.polybasis = self.polybasis
+        q.coeffs = homogeneizedToFull(tuple([self.N + 1] * self.npar),
+                                      self.npar, eV[:, 0])
         vbMng(self, "DEL", "Done computing denominator.", 7)
         return q
 
     def _setupNumerator(self):
         """Compute Pade' numerator."""
         RROMPyAssert(self._mode, message = "Cannot setup numerator.")
         vbMng(self, "INIT", "Starting computation of numerator.", 7)
         Qevaldiag = np.zeros((len(self.mus), len(self.mus)),
                              dtype = np.complex)
         verb = self.trainedModel.verbosity
         self.trainedModel.verbosity = 0
 
         self._setupInterpolationIndices()
         idxGlob = 0
         for j, derIdxs in enumerate(self._derIdxs):
             nder = len(derIdxs)
             idxLoc = np.arange(len(self.mus))[(self._reorder >= idxGlob)
                                             * (self._reorder < idxGlob + nder)]
             idxGlob += nder
             Qval = [0] * nder
             for der in range(nder):
                 derIdx = hashI(der, self.npar)
                 Qval[der] = (self.trainedModel.getQVal(
                                          self._musUnique[j], derIdx,
                                          scl = np.power(self.scaleFactor, -1.))
                            / multifactorial(derIdx))
             for derU, derUIdx in enumerate(derIdxs):
                 for derQ, derQIdx in enumerate(derIdxs):
                     diffIdx = [x - y for (x, y) in zip(derUIdx, derQIdx)]
                     if all([x >= 0 for x in diffIdx]):
                         diffj = hashD(diffIdx)
                         Qevaldiag[idxLoc[derU], idxLoc[derQ]] = Qval[diffj]
         if self.POD:
             Qevaldiag = Qevaldiag.dot(self.samplingEngine.RPOD.T)
         self.trainedModel.verbosity = verb
         while self.M >= 0:
             if self.radialBasis == 0:
-                fitVander, _, argIdxs = hpvP(self._musUniqueCN, self.M,
-                                         self.polybasisP, self._derIdxs,
-                                         self._reorder,
-                                         scl = np.power(self.scaleFactor, -1.))
-                fitnameEff = polyfitname(self.polybasisP)
-                nsamplesPrint = "{}".format(len(fitVander))
+                p = PI()
+                wellCond, msg = p.setupByInterpolation(
+                           self._musUniqueCN, Qevaldiag, self.M, 
+                           self.polybasisP, self.verbosity >= 5, True,
+                           {"derIdxs": self._derIdxs, "reorder": self._reorder,
+                            "scl": np.power(self.scaleFactor, -1.)},
+                           {"rcond": self.interpRcond})
             else:
-                fitVander, _, argIdxs = hpvR(self._musUniqueCN, self.M,
-                                        self.polybasisP, self._derIdxs,
-                                        self._reorder, self.radialBasisWeights,
-                                        scl = np.power(self.scaleFactor, -1.))
-                fitnameEff = polyfitnameR(self.polybasisP)
-                nConstraints = len(fitVander) - len(Qevaldiag)
-                if nConstraints > 0:
-                    Qevaldiag = np.pad(Qevaldiag, ((0, nConstraints), (0, 0)),
-                                       "constant")
-                elif nConstraints < 0:
-                    Qevaldiag = Qevaldiag[: len(fitVander)]
-                fitVander = fitVander[argIdxs]
-                nsamplesPrint = "{}+{}".format(len(self.mus),
-                                               len(fitVander) - len(self.mus))
-            fitVander = fitVander[:, argIdxs]
-            fitOut = customFit(fitVander, Qevaldiag, full = True,
-                               rcond = self.interpRcond)
-            vbMng(self, "MAIN",
-                  ("Fitting {} samples with degree {} through {}... "
-                   "Conditioning of LS system: {:.4e}.").format(
-                                           nsamplesPrint, self.M, fitnameEff,
-                                           fitOut[1][2][0] / fitOut[1][2][-1]),
-                  5)
-            if fitOut[1][1] == fitVander.shape[1]:
-                P = fitOut[0]
-                break
+                p = RBI()
+                wellCond, msg = p.setupByInterpolation(
+                         self._musUniqueCN, Qevaldiag, self.M, self.polybasisP,
+                         self.radialBasisWeights, self.verbosity >= 5, True,
+                         {"derIdxs": self._derIdxs, "reorder": self._reorder,
+                          "scl": np.power(self.scaleFactor, -1.)},
+                         {"rcond": self.interpRcond})
+            vbMng(self, "MAIN", msg, 5)
+            if wellCond: break
             RROMPyWarning("Polyfit is poorly conditioned. Reducing M by 1.")
             self.M -= 1
         if self.M < 0:
             raise RROMPyException(("Instability in computation of numerator. "
                                    "Aborting."))
-        if self.radialBasis == 0:
-            p = homogeneizedToFull(tuple([self.M + 1] * self.npar)
-                                 + (P.shape[1],), self.npar, P)
-        else:
-            pGlob = homogeneizedToFull(
-                               tuple([self.M + 1] * self.npar) + (P.shape[1],),
-                               self.npar, P[len(self.mus) :])
-            p = radialFunction(self._musUniqueCN[self._reorder],
-                               self.radialBasisWeights, 
-                               P[: len(self.mus)], pGlob)
         vbMng(self, "DEL", "Done computing numerator.", 7)
         return p
 
     def setupApprox(self):
         """
         Compute Pade' interpolant.
         SVD-based robust eigenvalue management.
         """
         if self.checkComputedApprox():
             return
         RROMPyAssert(self._mode, message = "Cannot setup approximant.")
         vbMng(self, "INIT", "Setting up {}.". format(self.name()), 5)
         self.computeScaleFactor()
         self.computeSnapshots()
         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,
                                     self.samplingEngine.samples,
                                     self.HFEngine.rescalingExp)
-            data.polytype = self.polybasis
-            data.polytypeP = self.polybasisP
             data.scaleFactor = self.scaleFactor
             data.mus = copy(self.mus)
             self.trainedModel.data = data
         else:
             self.trainedModel = self.trainedModel
             self.trainedModel.data.projMat = copy(self.samplingEngine.samples)
         if self.N > 0:
             Q = self._setupDenominator()
         else:
-            Q = np.ones(tuple([1] * self.npar), dtype = np.complex)
-        self.trainedModel.data.Q = copy(Q)
-        self.trainedModel.data.P = copy(self._setupNumerator())
+            Q = PI()
+            Q.coeffs = np.ones(tuple([1] * self.npar), dtype = np.complex)
+            Q.npar = self.npar
+            Q.polybasis = self.polybasis
+        self.trainedModel.data.Q = Q
+        self.trainedModel.data.P = self._setupNumerator()
         self.trainedModel.data.approxParameters = copy(self.approxParameters)
         vbMng(self, "DEL", "Done setting up approximant.", 5)
 
     def _computeInterpolantInverseBlocks(self) -> List[Np2D]:
         """
         Compute inverse factors for minimal interpolant target functional.
         """
         RROMPyAssert(self._mode, message = "Cannot solve eigenvalue problem.")
         self._setupInterpolationIndices()
         while self.E >= 0:
             eWidth = (hashD([self.E + 1] + [0] * (self.npar - 1))
                     - hashD([self.E] + [0] * (self.npar - 1)))
             TE, _, argIdxs = hpvP(self._musUniqueCN, self.E, self.polybasis,
                                   self._derIdxs, self._reorder, 
                                   scl = np.power(self.scaleFactor, -1.))
             fitOut = customPInv(TE[:, argIdxs], rcond = self.interpRcond,
                                 full = True)
             vbMng(self, "MAIN",
                   ("Fitting {} samples with degree {} through {}... "
                    "Conditioning of pseudoinverse system: {:.4e}.").format(
                                            TE.shape[0], self.E,
                                            polyfitname(self.polybasis),
                                            fitOut[1][1][0] / fitOut[1][1][-1]),
                   5)
             if fitOut[1][0] == len(argIdxs):
                 self._fitinv = fitOut[0][- eWidth : , :]
                 break
             RROMPyWarning("Polyfit is poorly conditioned. Reducing E by 1.")
             self.E -= 1
         if self.E < 0:
             raise RROMPyException(("Instability in computation of "
                                    "denominator. Aborting."))
         TN, _, argIdxs = hpvP(self._musUniqueCN, self.N, self.polybasis,
                               self._derIdxs, self._reorder, 
                               scl = np.power(self.scaleFactor, -1.))
         TN = TN[:, argIdxs]
         invD = [None] * (eWidth)
         for k in range(eWidth):
             pseudoInv = np.diag(self._fitinv[k, :])
             idxGlob = 0
             for j, derIdxs in enumerate(self._derIdxs):
                 nder = len(derIdxs)
                 idxGlob += nder
                 if nder > 1:
                     idxLoc = np.arange(len(self.mus))[
                                               (self._reorder >= idxGlob - nder)
                                             * (self._reorder < idxGlob)]
                     invLoc = self._fitinv[k, idxLoc]
                     pseudoInv[np.ix_(idxLoc, idxLoc)] = 0.
                     for diffj, diffjIdx in enumerate(derIdxs):
                         for derQ, derQIdx in enumerate(derIdxs):
                             derUIdx = [x - y for (x, y) in 
                                                     zip(diffjIdx, derQIdx)]
                             if all([x >= 0 for x in derUIdx]):
                                 derU = hashD(derUIdx)
                                 pseudoInv[idxLoc[derU], idxLoc[derQ]] = (
                                                              invLoc[diffj])
             invD[k] = pseudoInv.dot(TN)
         return invD
 
     def findeveVGExplicit(self, sampleE:sampList,
                           invD:List[Np2D]) -> Tuple[Np1D, Np2D]:
         """
         Compute explicitly eigenvalues and eigenvectors of Pade' denominator
             matrix.
         """
         RROMPyAssert(self._mode, message = "Cannot solve eigenvalue problem.")
         nEnd = invD[0].shape[1]
         eWidth = len(invD)
         vbMng(self, "INIT", "Building gramian matrix.", 10)
         gramian = self.HFEngine.innerProduct(sampleE, sampleE)
         G = np.zeros((nEnd, nEnd), dtype = np.complex)
         for k in range(eWidth):
              G += invD[k].T.conj().dot(gramian.dot(invD[k]))
         vbMng(self, "DEL", "Done building gramian.", 10)
         vbMng(self, "INIT", "Solving eigenvalue problem for gramian matrix.",
               7)
         ev, eV = np.linalg.eigh(G)
         vbMng(self, "MAIN",
               ("Solved eigenvalue problem of size {} with condition number "
                "{:.4e}.").format(nEnd, ev[-1] / ev[0]), 5)
         vbMng(self, "DEL", "Done solving eigenvalue problem.", 7)
         return ev, eV
 
     def findeveVGQR(self, RPODE:Np2D, invD:List[Np2D]) -> Tuple[Np1D, Np2D]:
         """
         Compute eigenvalues and eigenvectors of Pade' denominator matrix
             through SVD of R factor.
         """
         RROMPyAssert(self._mode, message = "Cannot solve eigenvalue problem.")
         nEnd = invD[0].shape[1]
         S = RPODE.shape[0]
         eWidth = len(invD)
         vbMng(self, "INIT", "Building half-gramian matrix stack.", 10)
         Rstack = np.zeros((S * eWidth, nEnd), dtype = np.complex)
         for k in range(eWidth):
             Rstack[k * S : (k + 1) * S, :] = RPODE.dot(invD[k])
         vbMng(self, "DEL", "Done building half-gramian.", 10)
         vbMng(self, "INIT", "Solving svd for square root of gramian matrix.",
               7)
         _, s, eV = np.linalg.svd(Rstack, full_matrices = False)
         ev = s[::-1]
         eV = eV[::-1, :].T.conj()
         vbMng(self, "MAIN", 
               ("Solved svd problem of size {} x {} with condition number "
                "{:.4e}.").format(*Rstack.shape, s[0] / s[-1]), 5)
         vbMng(self, "DEL", "Done solving svd.", 7)
         return ev, eV
 
     def centerNormalize(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.mu0.
 
         Returns:
             Normalized parameter.
         """
         return self.trainedModel.centerNormalize(mu, mu0)
 
     def getResidues(self) -> Np1D:
         """
         Obtain approximant residues.
 
         Returns:
             Matrix with residues as columns.
         """
         return self.trainedModel.getResidues()
 
diff --git a/rrompy/reduction_methods/distributed_greedy/rational_interpolant_greedy.py b/rrompy/reduction_methods/distributed_greedy/rational_interpolant_greedy.py
index b777bb0..0b2a6f2 100644
--- a/rrompy/reduction_methods/distributed_greedy/rational_interpolant_greedy.py
+++ b/rrompy/reduction_methods/distributed_greedy/rational_interpolant_greedy.py
@@ -1,414 +1,410 @@
 # 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_greedy_approximant import \
                                             GenericDistributedGreedyApproximant
-from rrompy.utilities.poly_fitting.polynomial import polybases, polydomcoeff
+from rrompy.utilities.poly_fitting.polynomial import polybases
 from rrompy.reduction_methods.distributed import RationalInterpolant
 from rrompy.reduction_methods.trained_model import TrainedModelPade as tModel
 from rrompy.reduction_methods.trained_model import TrainedModelData
 from rrompy.utilities.base.types import Np1D, Np2D, DictAny, HFEng, paramVal
 from rrompy.utilities.base import verbosityManager as vbMng
 from rrompy.utilities.exception_manager import RROMPyWarning
 from rrompy.utilities.exception_manager import RROMPyException, RROMPyAssert
 
 __all__ = ['RationalInterpolantGreedy']
 
 class RationalInterpolantGreedy(GenericDistributedGreedyApproximant,
                                 RationalInterpolant):
     """
     ROM greedy rational interpolant 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': number of starting training points;
             - 'sampler': sample point generator;
             - 'radialBasis': radial basis family for interpolant numerator;
                 defaults to 0, i.e. no radial basis;
             - 'radialBasisWeights': radial basis weights for interpolant
                 numerator; defaults to 0, i.e. identity;
             - 'greedyTol': uniform error tolerance for greedy algorithm;
                 defaults to 1e-2;
             - 'interactive': whether to interactively terminate greedy
                 algorithm; defaults to False;
             - 'maxIter': maximum number of greedy steps; defaults to 1e2;
             - 'refinementRatio': ratio of training points to be exhausted
                 before training set refinement; defaults to 0.2;
             - 'nTestPoints': number of test points; defaults to 5e2;
             - 'trainSetGenerator': training sample points generator; defaults
                 to Chebyshev sampler within muBounds;
             - 'polybasis': type of basis for interpolation; allowed values
                 include 'MONOMIAL', 'CHEBYSHEV' and 'LEGENDRE'; defaults to
                 'MONOMIAL';
             - 'Delta': difference between M and N in rational approximant;
                 defaults to 0;
             - 'errorEstimatorKind': kind of error estimator; available values
                 include 'EXACT', 'BASIC', and 'BARE'; defaults to 'EXACT';
             - 'interpRcond': tolerance for interpolation; defaults to None;
             - 'robustTol': tolerance for robust Pade' denominator management;
                 defaults to 0.
             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.
         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.
             - 'radialBasis': radial basis family for interpolant numerator;
                 defaults to 0, i.e. no radial basis;
             - 'radialBasisWeights': radial basis weights for interpolant
                 numerator; defaults to 0, i.e. identity;
             - 'greedyTol': uniform error tolerance for greedy algorithm;
             - 'interactive': whether to interactively terminate greedy
                 algorithm;
             - 'maxIter': maximum number of greedy steps;
             - 'refinementRatio': ratio of test points to be exhausted before
                 test set refinement;
             - 'nTestPoints': number of test points;
             - 'trainSetGenerator': training sample points generator;
             - 'Delta': difference between M and N in rational approximant;
             - 'errorEstimatorKind': kind of error estimator;
             - 'interpRcond': tolerance for interpolation;
             - 'robustTol': tolerance for robust Pade' denominator management.
         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 test points.
         sampler: Sample point generator.
         radialBasis: Radial basis family for interpolant numerator.
         radialBasisWeights: Radial basis weights for interpolant numerator.
         greedyTol: uniform error tolerance for greedy algorithm.
         interactive: whether to interactively terminate greedy algorithm.
         maxIter: maximum number of greedy steps.
         refinementRatio: ratio of training points to be exhausted before
             training set refinement.
         nTestPoints: number of starting training points.
         trainSetGenerator: training sample points generator.
         robustTol: tolerance for robust Pade' denominator management.
         Delta: difference between M and N in rational approximant.
         errorEstimatorKind: kind of error estimator.
         interpRcond: tolerance for interpolation.
         robustTol: tolerance for robust Pade' denominator management.
         muBounds: list of bounds for parameter values.
         samplingEngine: Sampling engine.
         estimatorNormEngine: Engine for estimator norm computation.
         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.
     """
     
     _allowedEstimatorKinds = ["EXACT", "BASIC", "BARE"]
     
     def __init__(self, HFEngine:HFEng, mu0 : paramVal = None,
                  approxParameters : DictAny = {}, homogeneized : bool = False,
                  verbosity : int = 10, timestamp : bool = True):
         self._preInit()
         self._addParametersToList(["Delta", "polybasis", "errorEstimatorKind",
                                    "interpRcond", "robustTol"], 
                                   [0, "MONOMIAL", "EXACT", -1, 0],
                                   toBeExcluded = ["E"])
         super().__init__(HFEngine = HFEngine, mu0 = mu0,
                          approxParameters = approxParameters,
                          homogeneized = homogeneized,
                          verbosity = verbosity, timestamp = timestamp)
         vbMng(self, "INIT", "Computing Taylor blocks of system.", 7)
         nAs = self.HFEngine.nAs - 1
         nbs = max(self.HFEngine.nbs, (nAs + 1) * self.homogeneized)
         self.As = [self.HFEngine.A(self.mu0, j + 1) for j in range(nAs)]
         self.bs = [self.HFEngine.b(self.mu0, j, self.homogeneized)
                                                        for j in range(nbs)]
         vbMng(self, "DEL", "Done computing Taylor blocks.", 7)
         self._postInit()
 
     @property
     def polybasis(self):
         """Value of polybasis."""
         return self._polybasis
     @polybasis.setter
     def polybasis(self, polybasis):
         try:
             polybasis = polybasis.upper().strip().replace(" ","")
             if polybasis not in polybases: 
                 raise RROMPyException("Sample type not recognized.")
             self._polybasis = polybasis
         except:
             RROMPyWarning(("Prescribed polybasis not recognized. Overriding "
                            "to 'MONOMIAL'."))
             self._polybasis = "MONOMIAL"
         self._approxParameters["polybasis"] = self.polybasis
 
     @property
     def Delta(self):
         """Value of Delta."""
         return self._Delta
     @Delta.setter
     def Delta(self, Delta):
         if not np.isclose(Delta, np.floor(Delta)):
             raise RROMPyException("Delta must be an integer.")
         if Delta < 0:
             RROMPyWarning(("Error estimator unreliable for Delta < 0. "
                            "Overloading of errorEstimator is suggested."))
         else:
             Deltamin = (max(self.HFEngine.nbs,
                             self.HFEngine.nAs * self.homogeneized)
                       - 1 - 1 * (self.HFEngine.nAs > 1))
             if Delta < Deltamin:
                 RROMPyWarning(("Method may be unreliable for selected Delta. "
                                "Suggested minimal value of Delta: {}.").format(
                                                                      Deltamin))
         self._Delta = Delta
         self._approxParameters["Delta"] = self.Delta
 
     @property
     def errorEstimatorKind(self):
         """Value of errorEstimatorKind."""
         return self._errorEstimatorKind
     @errorEstimatorKind.setter
     def errorEstimatorKind(self, errorEstimatorKind):
         errorEstimatorKind = errorEstimatorKind.upper()
         if errorEstimatorKind not in self._allowedEstimatorKinds:
             RROMPyWarning(("Error estimator kind not recognized. Overriding "
                            "to 'EXACT'."))
             errorEstimatorKind = "EXACT"
         self._errorEstimatorKind = errorEstimatorKind
         self._approxParameters["errorEstimatorKind"] = self.errorEstimatorKind
 
     @property
     def nTestPoints(self):
         """Value of nTestPoints."""
         return self._nTestPoints
     @nTestPoints.setter
     def nTestPoints(self, nTestPoints):
         if nTestPoints <= np.abs(self.Delta):
             RROMPyWarning(("nTestPoints must be at least abs(Delta) + 1. "
                            "Increasing value to abs(Delta) + 1."))
             nTestPoints = np.abs(self.Delta) + 1
         if not np.isclose(nTestPoints, np.int(nTestPoints)):
             raise RROMPyException("nTestPoints must be an integer.")
         nTestPoints = np.int(nTestPoints)
         if hasattr(self, "_nTestPoints") and self.nTestPoints is not None:
             nTestPointsold = self.nTestPoints
         else: nTestPointsold = -1
         self._nTestPoints = nTestPoints
         self._approxParameters["nTestPoints"] = self.nTestPoints
         if nTestPointsold != self.nTestPoints:
             self.resetSamples()
 
     def _errorSamplingRatio(self, mus:Np1D, muRTest:Np1D,
                             muRTrain:Np1D) -> Np1D:
         """Scalar ratio in explicit error estimator."""
         self.setupApprox()
         testTile = np.tile(np.reshape(muRTest, (-1, 1)), [1, len(muRTrain)])
         nodalVals = np.prod(testTile - np.reshape(muRTrain, (1, -1)), axis = 1)
         denVals = self.trainedModel.getQVal(mus)
         return np.abs(nodalVals / denVals)
 
     def _RHSNorms(self, radiusb0:Np2D) -> Np1D:
         """High fidelity system RHS norms."""
         self.assembleReducedResidualBlocks(full = False)
         # 'ij,jk,ik->k', resbb, radiusb0, radiusb0.conj()
         b0resb0 = np.sum(self.trainedModel.data.resbb.dot(radiusb0)
                        * radiusb0.conj(), axis = 0)
         RHSnorms = np.power(np.abs(b0resb0), .5)
         return RHSnorms
 
     def _errorEstimatorBare(self) -> Np1D:
         """Bare residual-based error estimator."""
         self.setupApprox()
         self.assembleReducedResidualGramian(self.trainedModel.data.projMat)
-        pDom = self.trainedModel.data.P[-1, :]
+        pDom = self.trainedModel.data.P.domCoeff
         LL = pDom.conj().dot(self.trainedModel.data.gramian.dot(pDom))
         Adiag = self.As[0].diagonal()
         LL = ((self.scaleFactor[0] * np.linalg.norm(Adiag)) ** 2. * LL
             / np.size(Adiag))
-        scalingDom = polydomcoeff(len(self.mus) - 1, self.polybasis)
-        return scalingDom * np.power(np.abs(LL), .5)
+        return np.power(np.abs(LL), .5)
 
     def _errorEstimatorBasic(self, muTest:paramVal, ratioTest:complex) -> Np1D:
         """Basic residual-based error estimator."""
         resmu = self.HFEngine.residual(self.trainedModel.getApprox(muTest),
                                        muTest, self.homogeneized,
                                        duality = False)[0]
         return np.abs(self.estimatorNormEngine.norm(resmu) / ratioTest)
 
     def _errorEstimatorExact(self, muRTrain:Np1D, vanderBase:Np2D) -> Np1D:
         """Exact residual-based error estimator."""
         self.setupApprox()
         self.assembleReducedResidualBlocks(full = True)
         nAs = self.HFEngine.nAs - 1
         nbs = max(self.HFEngine.nbs - 1, nAs * self.homogeneized)
-        delta = len(self.mus) - len(self.trainedModel.data.Q)
+        delta = len(self.mus) - self.trainedModel.data.Q.deg[0]
         nbsEff = max(0, nbs - delta)
         momentQ = np.zeros(nbsEff, dtype = np.complex)
         momentQu = np.zeros((len(self.mus), nAs), dtype = np.complex)
         radiusbTen = np.zeros((nbsEff, nbsEff, vanderBase.shape[1]),
                               dtype = np.complex)
         radiusATen = np.zeros((nAs, nAs, vanderBase.shape[1]),
                               dtype = np.complex)
         if nbsEff > 0:
-            momentQ[0] = self.trainedModel.data.Q[-1]
+            momentQ[0] = self.trainedModel.data.Q.domCoeff
             radiusbTen[0, :, :] = vanderBase[: nbsEff, :]
-        momentQu[:, 0] = self.trainedModel.data.P[-1, :]
+        momentQu[:, 0] = self.trainedModel.data.P.domCoeff
         radiusATen[0, :, :] = vanderBase[: nAs, :]
         Qvals = self.trainedModel.getQVal(self.mus)
         for k in range(1, max(nAs, nbs * (nbsEff > 0))):
             Qvals = Qvals * muRTrain
             if k > delta and k < nbs:
                 momentQ[k - delta] = self._fitinv.dot(Qvals)
                 radiusbTen[k - delta, k :, :] = (
                                          radiusbTen[0, : delta - k, :])
             if k < nAs:
                 momentQu[:, k] = Qvals * self._fitinv
                 radiusATen[k, k :, :] = radiusATen[0, : - k, :]
         if self.POD and nAs > 1:
             momentQu[:, 1 :] = self.samplingEngine.RPOD.dot(
                                                       momentQu[:, 1 :])
         radiusA = np.tensordot(momentQu, radiusATen, 1)
         if nbsEff > 0:
             radiusb = np.tensordot(momentQ, radiusbTen, 1)
             # 'ij,jk,ik->k', resbb, radiusb, radiusb.conj()
             ff = np.sum(self.trainedModel.data.resbb[delta + 1 :, delta + 1 :]\
                         .dot(radiusb) * radiusb.conj(), axis = 0)
             # 'ijk,jkl,il->l', resAb, radiusA, radiusb.conj()
             Lf = np.sum(np.tensordot(
                        self.trainedModel.data.resAb[delta :, :, :], radiusA, 2)
                       * radiusb.conj(), axis = 0)
         else:
             ff, Lf = 0., 0.
         # 'ijkl,klt,ijt->t', resAA, radiusA, radiusA.conj()
         LL = np.sum(np.tensordot(self.trainedModel.data.resAA, radiusA, 2)
                   * radiusA.conj(), axis = (0, 1))
-        scalingDom = polydomcoeff(len(self.mus) - 1, self.polybasis)
-        return scalingDom * np.power(np.abs(ff - 2. * np.real(Lf) + LL), .5)
+        return np.power(np.abs(ff - 2. * np.real(Lf) + LL), .5)
 
     def errorEstimator(self, mus:Np1D) -> Np1D:
         """Standard residual-based error estimator."""
         self.setupApprox()
-        if (np.any(np.isnan(self.trainedModel.data.P[-1, :]))
-         or np.any(np.isinf(self.trainedModel.data.P[-1, :]))):
+        if (np.any(np.isnan(self.trainedModel.data.P.domCoeff))
+         or np.any(np.isinf(self.trainedModel.data.P.domCoeff))):
             err = np.empty(len(mus))
             err[:] = np.inf
             return err
         nAs = self.HFEngine.nAs - 1
         nbs = max(self.HFEngine.nbs - 1, nAs * self.homogeneized)
         muRTest = self.centerNormalize(mus)(0)
         muRTrain = self.centerNormalize(self.mus)(0)
         samplingRatio = self._errorSamplingRatio(mus, muRTest, muRTrain)
         vanderBase = np.polynomial.polynomial.polyvander(muRTest,
                                                          max(nAs, nbs)).T
         RHSnorms = self._RHSNorms(vanderBase[: nbs + 1, :])
         if self.errorEstimatorKind == "BARE":
             jOpt = self._errorEstimatorBare()
         elif self.errorEstimatorKind == "BASIC":
             idx_muTestSample = np.argmax(samplingRatio)
             jOpt = self._errorEstimatorBasic(mus[idx_muTestSample],
                                              samplingRatio[idx_muTestSample])
         else: #if self.errorEstimatorKind == "EXACT":
             jOpt = self._errorEstimatorExact(muRTrain, vanderBase[: -1, :])
         return jOpt * samplingRatio / RHSnorms
 
     def setupApprox(self, plotEst : bool = False):
         """
         Compute Pade' interpolant.
         SVD-based robust eigenvalue management.
         """
         if self.checkComputedApprox():
             return
         RROMPyAssert(self._mode, message = "Cannot setup approximant.")
         vbMng(self, "INIT", "Setting up {}.". format(self.name()), 5)
         self.computeScaleFactor()
         self.greedy(plotEst)
         self._S = len(self.mus)
         self._N, self._M, self._E = [self._S - 1] * 3
         if self.Delta < 0:
             self._M += self.Delta
         else:
             self._N -= self.Delta
         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,
                                     self.samplingEngine.samples,
                                     self.HFEngine.rescalingExp)
-            data.polytype = self.polybasis
-            data.polytypeP = self.polybasisP
             data.scaleFactor = self.scaleFactor
             data.mus = copy(self.mus)
             self.trainedModel.data = data
         else:
             self.trainedModel = self.trainedModel
             self.trainedModel.data.projMat = copy(self.samplingEngine.samples)
             self.trainedModel.data.mus = copy(self.mus)
         if min(self.M, self.N) < 0:
             vbMng(self, "MAIN", "Minimal sample size not achieved.", 5)
             Q = np.empty(tuple([max(self.N, 0) + 1] * self.npar),
                          dtype = np.complex)
             P = np.empty(tuple([max(self.M, 0) + 1] * self.npar)
                        + (len(self.mus),), dtype = np.complex)
             Q[:] = np.nan
             P[:] = np.nan
             self.trainedModel.data.Q = copy(Q)
             self.trainedModel.data.P = copy(P)
             self.trainedModel.data.approxParameters = copy(
                                                          self.approxParameters)
             vbMng(self, "DEL", "Aborting computation of approximant.", 5)
             return
         if self.N > 0:
             Q = self._setupDenominator()
         else:
             Q = np.ones(tuple([1] * self.npar), dtype = np.complex)
         self.trainedModel.data.Q = copy(Q)
         self.trainedModel.data.P = copy(self._setupNumerator())
         self.trainedModel.data.approxParameters = copy(self.approxParameters)
         vbMng(self, "DEL", "Done setting up approximant.", 5)
 
     def assembleReducedResidualBlocks(self, full : bool = False):
         """Build affine blocks of reduced linear system through projections."""
         scaling = self.trainedModel.data.scaleFactor[0]
         self.assembleReducedResidualBlocksbb(self.bs, scaling)
         if full:
             pMat = self.trainedModel.data.projMat
             self.assembleReducedResidualBlocksAb(self.As, self.bs[1 :],
                                                  pMat, scaling)
             self.assembleReducedResidualBlocksAA(self.As, pMat, scaling)
 
diff --git a/rrompy/reduction_methods/trained_model/trained_model_pade.py b/rrompy/reduction_methods/trained_model/trained_model_pade.py
index 6e4b392..cbda031 100644
--- a/rrompy/reduction_methods/trained_model/trained_model_pade.py
+++ b/rrompy/reduction_methods/trained_model/trained_model_pade.py
@@ -1,159 +1,156 @@
 # Copyright (C) 2018 by the RROMPy authors
 #
 # This file is part of RROMPy.
 #
 # RROMPy is free software: you can redistribute it and/or modify
 # it under the terms of the GNU Lesser General Public License as published by
 # the Free Software Foundation, either version 3 of the License, or
 # (at your option) any later version.
 #
 # RROMPy is distributed in the hope that it will be useful,
 # but WITHOUT ANY WARRANTY; without even the implied warranty of
 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
 # GNU Lesser General Public License for more details.
 #
 # You should have received a copy of the GNU Lesser General Public License
 # along with RROMPy. If not, see <http://www.gnu.org/licenses/>.
 #
 
 import numpy as np
 from . import TrainedModel
 from rrompy.utilities.base.types import (Np1D, List, paramVal, paramList,
                                          sampList)
 from rrompy.utilities.base import verbosityManager as vbMng
 from rrompy.utilities.poly_fitting.polynomial import polyval, polyroots
 from rrompy.utilities.poly_fitting.radial_basis import polyval as polyvalR
-from rrompy.utilities.exception_manager import RROMPyAssert, RROMPyException
+from rrompy.utilities.exception_manager import RROMPyException
 from rrompy.parameter import checkParameter, checkParameterList
 from rrompy.sampling import sampleList
 
 __all__ = ['TrainedModelPade']
 
 class TrainedModelPade(TrainedModel):
     """
     ROM approximant evaluation for Pade' approximant.
     
     Attributes:
         Data: dictionary with all that can be pickled.
     """
 
     def centerNormalize(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.mu0.
 
         Returns:
             Normalized parameter.
         """
         mu, _ = checkParameterList(mu, self.data.npar)
         if mu0 is None: mu0 = self.data.mu0
         rad = ((mu ** self.data.rescalingExp - mu0 ** self.data.rescalingExp)
              / self.data.scaleFactor)
         return rad
 
-    def getPVal(self, mu : paramList = [], der : List[int] = None) -> sampList:
+    def getPVal(self, mu : paramList = [], der : List[int] = None,
+                scl : Np1D = None) -> sampList:
         """
         Evaluate Pade' numerator at arbitrary parameter.
 
         Args:
             mu: Target parameter.
             der(optional): Derivatives to take before evaluation.
         """
         mu, _ = checkParameterList(mu, self.data.npar)
         vbMng(self, "INIT", "Evaluating numerator at mu = {}.".format(mu), 17)
         muCenter = self.centerNormalize(mu)
-        if "_" in self.data.polytypeP:
-            p = sampleList(polyvalR(muCenter, self.data.P,
-                                    self.data.polytypeP, der))
-        else:
-            p = sampleList(polyval(muCenter, self.data.P,
-                                   self.data.polytypeP, der))
+        p = sampleList(self.data.P(muCenter, der, scl))
         vbMng(self, "DEL", "Done evaluating numerator.", 17)
         return p
 
     def getQVal(self, mu:Np1D, der : List[int] = None,
                 scl : Np1D = None) -> Np1D:
         """
         Evaluate Pade' denominator at arbitrary parameter.
 
         Args:
             mu: Target parameter.
             der(optional): Derivatives to take before evaluation.
         """
         mu, _ = checkParameterList(mu, self.data.npar)
         vbMng(self, "INIT", "Evaluating denominator at mu = {}.".format(mu),
               17)
         muCenter = self.centerNormalize(mu)
-        q = polyval(muCenter, self.data.Q, self.data.polytype, der, scl)
+        q = self.data.Q(muCenter, der, scl)
         vbMng(self, "DEL", "Done evaluating denominator.", 17)
         return q
 
     def getApproxReduced(self, mu : paramList = []) -> sampList:
         """
         Evaluate reduced representation of approximant at arbitrary parameter.
 
         Args:
             mu: Target parameter.
         """
         mu, _ = checkParameterList(mu, self.data.npar)
         if (not hasattr(self, "lastSolvedApproxReduced")
          or self.lastSolvedApproxReduced != mu):
             vbMng(self, "INIT",
                   "Evaluating approximant at mu = {}.".format(mu), 12)
             self.uApproxReduced = self.getPVal(mu) / self.getQVal(mu)
             vbMng(self, "DEL", "Done evaluating approximant.", 12)
             self.lastSolvedApproxReduced = mu
         return self.uApproxReduced
     
     def getPoles(self, *args, **kwargs) -> Np1D:
         """
         Obtain approximant poles.
 
         Returns:
             Numpy complex vector of poles.
         """
         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 hasattr(mVals, "__len__"): mVals = [mVals]
             mVals = list(mVals)
         else:
             mVals = [None]
         try:
             rDim = mVals.index(None)
             if rDim < len(mVals) - 1 and None in mVals[rDim + 1 :]:
                 raise
         except:
             raise RROMPyException(("Exactly 1 'None' entry in "
                                    "marginalVals must be provided."))
         mVals[rDim] = self.data.mu0(rDim)
         mVals = self.centerNormalize(checkParameter(mVals, len(mVals)))
         mVals = list(mVals.data.flatten())
         mVals[rDim] = None
         return np.power(self.data.mu0(rDim) ** self.data.rescalingExp[rDim]
                       + self.data.scaleFactor[rDim]
-                      * polyroots(self.data.Q, self.data.polytype, mVals),
+                      * polyroots(self.data.Q.coeffs,
+                                  self.data.Q.polybasis, mVals),
                         1. / self.data.rescalingExp[rDim])
 
     def getResidues(self) -> Np1D:
         """
         Obtain approximant residues.
 
         Returns:
             Numpy matrix with residues as columns.
         """
         pls = self.getPoles()
         poles, _ = checkParameterList(pls, 1)
         res = (self.data.projMat.dot(self.getPVal(poles).data)
              / self.getQVal(poles, 1))
         return pls, res
 
diff --git a/rrompy/utilities/base/types.py b/rrompy/utilities/base/types.py
index b84e189..06f9835 100644
--- a/rrompy/utilities/base/types.py
+++ b/rrompy/utilities/base/types.py
@@ -1,61 +1,60 @@
 # 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 typing import TypeVar, List, Tuple, Dict, Any
 
 __all__ = ['TupleAny','ListAny','DictAny','ScOp','Np1D','Np2D','Np1DLst',
            'N2FSExpr','FenExpr','FenFunc','FenFuncSpace','FenBC','HFEng',
            'ROMEng','sampleEng','normEng','paramVal','paramList', 'sampList',
            'GenExpr','strLst', 'BfSExpr']
 
 # ANY
 TupleAny = Tuple[Any]
 ListAny = List[Any]
 DictAny = Dict[Any, Any]
 
 # SCIPY
 ScOp = TypeVar("Scipy sparse matrix for space operator")
 
 # NUMPY
 Np1D = TypeVar("NumPy 1D array")
 Np2D = TypeVar("NumPy 2D array-like")
 Np1DLst = TypeVar("NumPy 1D array or list of NumPy 1D array")
 N2FSExpr = TypeVar("NumPy 2D array, float or str")
 
 # FENICS
 FenExpr = TypeVar("FEniCS expression")
 FenFunc = TypeVar("FEniCS function")
 FenFuncSpace = TypeVar("FEniCS function space")
 FenBC = TypeVar("FEniCS boundary condition")
 
 # ENGINES
 HFEng = TypeVar("High fidelity engine")
 ROMEng = TypeVar("ROM engine")
 sampleEng = TypeVar("Sampling engine")
 normEng = TypeVar("Norm engine")
 
 # CUSTOM TYPES
 paramVal = TypeVar("Parameter value tuple")
 paramList = TypeVar("Parameter value tuple list")
 sampList = TypeVar("Sample list")
-radialFun = TypeVar("Radial basis function")
 
 # OTHERS
 GenExpr = TypeVar("Generic expression")
 strLst = TypeVar("str or list of str")
 BfSExpr = TypeVar("Boolean function or string")
diff --git a/rrompy/utilities/poly_fitting/__init__.py b/rrompy/utilities/poly_fitting/__init__.py
index c9d4457..0e297b8 100644
--- a/rrompy/utilities/poly_fitting/__init__.py
+++ b/rrompy/utilities/poly_fitting/__init__.py
@@ -1,27 +1,29 @@
 # 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 .custom_fit import customFit
 from .custom_pinv import customPInv
+from .interpolator import GenericInterpolator
 
 __all__ = [
         'customFit',
-        'customPInv'
+        'customPInv',
+        'GenericInterpolator'
            ]
 
 
diff --git a/rrompy/utilities/poly_fitting/radial_basis/der.py b/rrompy/utilities/poly_fitting/interpolator.py
similarity index 61%
rename from rrompy/utilities/poly_fitting/radial_basis/der.py
rename to rrompy/utilities/poly_fitting/interpolator.py
index 89f2aed..fd08cd3 100644
--- a/rrompy/utilities/poly_fitting/radial_basis/der.py
+++ b/rrompy/utilities/poly_fitting/interpolator.py
@@ -1,30 +1,41 @@
 # Copyright (C) 2018 by the RROMPy authors
 #
 # This file is part of RROMPy.
 #
 # RROMPy is free software: you can redistribute it and/or modify
 # it under the terms of the GNU Lesser General Public License as published by
 # the Free Software Foundation, either version 3 of the License, or
 # (at your option) any later version.
 #
 # RROMPy is distributed in the hope that it will be useful,
 # but WITHOUT ANY WARRANTY; without even the implied warranty of
 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
 # GNU Lesser General Public License for more details.
 #
 # You should have received a copy of the GNU Lesser General Public License
 # along with RROMPy. If not, see <http://www.gnu.org/licenses/>.
 #
 
-import numpy as np
-from rrompy.utilities.base.types import Np1D, List, radialFun
-from rrompy.utilities.exception_manager import RROMPyException
+from abc import abstractmethod
+from rrompy.utilities.base.types import List, paramList
 
-__all__ = ['polyder']
+__all__ = ['GenericInterpolator']
 
-def polyder(c:radialFun, basis:str, m : List[int] = None,
-            scl : Np1D = None) -> radialFun:
-    if m is not None and np.sum(m) > 0:
-        raise RROMPyException(("Cannot take derivatives of radial basis "
-                               "function."))
-    return c
+class GenericInterpolator:
+    """HERE"""
+
+    @abstractmethod
+    def __init__(self, other = None):
+        pass
+
+    @abstractmethod
+    def __call__(self, mu:paramList, der : List[int] = None):
+        pass
+
+    @abstractmethod
+    def __copy__(self):
+        pass
+
+    @abstractmethod
+    def __deepcopy__(self, memo):
+        pass
diff --git a/rrompy/utilities/poly_fitting/polynomial/__init__.py b/rrompy/utilities/poly_fitting/polynomial/__init__.py
index 7f8faec..747c8b5 100644
--- a/rrompy/utilities/poly_fitting/polynomial/__init__.py
+++ b/rrompy/utilities/poly_fitting/polynomial/__init__.py
@@ -1,47 +1,49 @@
 # 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 .base import (polybases, polyfitname, polydomcoeff)
 from .der import polyder
 from .val import polyval
 from .marginalize import polymarginalize
 from .vander import polyvander
 from .roots import polyroots
 from .derivative import nextDerivativeIndices
 from .hash_derivative import hashDerivativeToIdx, hashIdxToDerivative
 from .homogeneization import (homogeneizationMask, homogeneizedpolyvander,
                               homogeneizedToFull)
+from .polynomial_interpolator import PolynomialInterpolator
 
 __all__ = [
         'polybases',
         'polyfitname',
         'polydomcoeff',
         'polyder',
         'polyval',
         'polymarginalize',
         'polyvander',
         'polyroots',
         'nextDerivativeIndices',
         'hashDerivativeToIdx',
         'hashIdxToDerivative',
         'homogeneizationMask',
         'homogeneizedpolyvander',
-        'homogeneizedToFull'
+        'homogeneizedToFull',
+        'PolynomialInterpolator'
            ]
 
 
diff --git a/rrompy/utilities/poly_fitting/polynomial/polynomial_interpolator.py b/rrompy/utilities/poly_fitting/polynomial/polynomial_interpolator.py
new file mode 100644
index 0000000..6a864be
--- /dev/null
+++ b/rrompy/utilities/poly_fitting/polynomial/polynomial_interpolator.py
@@ -0,0 +1,124 @@
+# Copyright (C) 2018 by the RROMPy authors
+#
+# This file is part of RROMPy.
+#
+# RROMPy is free software: you can redistribute it and/or modify
+# it under the terms of the GNU Lesser General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# RROMPy is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public License
+# along with RROMPy. If not, see <http://www.gnu.org/licenses/>.
+#
+
+from copy import deepcopy as copy
+from rrompy.utilities.base.types import List, ListAny, DictAny, Np1D, paramList
+from rrompy.utilities.poly_fitting.interpolator import GenericInterpolator
+from rrompy.utilities.poly_fitting.custom_fit import customFit
+from rrompy.utilities.poly_fitting.polynomial.base import (polyfitname,
+                                                           polydomcoeff)
+from rrompy.utilities.poly_fitting.polynomial.val import polyval
+from rrompy.utilities.poly_fitting.polynomial.roots import polyroots
+from rrompy.utilities.poly_fitting.polynomial.homogeneization import (
+                                                 homogeneizedpolyvander as hpv,
+                                                 homogeneizedToFull)
+from rrompy.utilities.poly_fitting.polynomial.vander import polyvander as pv
+from rrompy.utilities.exception_manager import RROMPyAssert, RROMPyException
+from rrompy.parameter import checkParameterList
+
+__all__ = ['PolynomialInterpolator']
+
+class PolynomialInterpolator(GenericInterpolator):
+    """HERE"""
+
+    def __init__(self, other = None):
+        if other is None: return
+        self.coeffs = other.coeffs
+        self.npar = other.npar
+        self.polybasis = other.polybasis
+
+    @property
+    def shape(self):
+        if self.coeffs.ndim > self.npar:
+            sh = self.coeffs.shape[self.npar :]
+        else: sh = 1
+        return sh
+
+    @property
+    def deg(self):
+        return [x - 1 for x in self.coeffs.shape[: self.npar]]
+
+    def __getitem__(self, key):
+        return self.coeffs[key]
+
+    def __call__(self, mu:paramList, der : List[int] = None,
+                 scl : Np1D = None):
+        return polyval(mu, self.coeffs, self.polybasis, der, scl)
+
+    def __copy__(self):
+        return PolynomialInterpolator(self)
+
+    def __deepcopy__(self, memo):
+        other = PolynomialInterpolator()
+        other.coeffs, other.npar, other.polybasis = copy(
+                                (self.coeffs, self.npar, self.polybasis), memo)
+        return other
+
+    @property
+    def domCoeff(self):
+        RROMPyAssert(self.npar, 1, "Number of parameters")
+        return polydomcoeff(self.deg, self.polybasis) * self[-1]
+
+    def setupByInterpolation(self, support:paramList, values:ListAny,
+                             deg:int, polybasis : str = "MONOMIAL",
+                             verbose : bool = True, homogeneized : bool = True,
+                             vanderCoeffs : DictAny = {},
+                             fitCoeffs : DictAny = {}):
+        support, _ = checkParameterList(support)
+        self.npar = support.shape[1]
+        self.polybasis = polybasis
+        if homogeneized:
+            vander, _, reorder = hpv(support, deg = deg, basis = polybasis,
+                                     **vanderCoeffs)
+            vander = vander[:, reorder]
+        else:
+            if not hasattr(deg, "__len__"): deg = [deg] * self.npar
+            vander = pv(support, deg = deg, basis = polybasis,
+                        **vanderCoeffs)
+        outDim = values.shape[1:]
+        values = values.reshape(values.shape[0], -1)
+        fitOut = customFit(vander, values, full = True, **fitCoeffs)
+        P = fitOut[0]
+        if verbose:
+            msg = ("Fitting {} samples with degree {} through {}... "
+                   "Conditioning of LS system: {:.4e}.").format(
+                                            len(vander), deg,
+                                            polyfitname(self.polybasis),
+                                            fitOut[1][2][0] / fitOut[1][2][-1])
+        else: msg = None
+        if homogeneized:
+            self.coeffs = homogeneizedToFull(
+                                         tuple([deg + 1] * self.npar) + outDim,
+                                         self.npar, P)
+        else:
+            self.coeffs = P.reshape(tuple([d + 1 for d in deg]) + outDim)
+        return fitOut[1][1] == vander.shape[1], msg
+
+    def roots(self, marginalVals : ListAny = [None]):
+        RROMPyAssert(self.shape, 1, "Shape of output")
+        RROMPyAssert(len(ListAny), self.npar, "Number of parameters")
+        try:
+            rDim = marginalVals.index(None)
+            if (rDim < len(marginalVals) - 1
+            and None in marginalVals[rDim + 1 :]):
+                raise
+        except:
+            raise RROMPyException(("Exactly 1 'None' entry in "
+                                   "marginalVals must be provided."))
+        return polyroots(self.coeffs, self.polybasis, marginalVals)
+
diff --git a/rrompy/utilities/poly_fitting/radial_basis/__init__.py b/rrompy/utilities/poly_fitting/radial_basis/__init__.py
index e4d799d..108132b 100644
--- a/rrompy/utilities/poly_fitting/radial_basis/__init__.py
+++ b/rrompy/utilities/poly_fitting/radial_basis/__init__.py
@@ -1,42 +1,41 @@
 # 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 .kernel import radialGaussian, thinPlateSpline, multiQuadric
-from .base import rbbases, polybases, polyfitname, polydomcoeff, radialFunction
-from .der import polyder
+from .base import rbbases, polybases, polyfitname, polydomcoeff
 from .val import polyval
 from .vander import rbvander, polyvander
 from .homogeneization import homogeneizedpolyvander
+from .radial_basis_interpolator import RadialBasisInterpolator
 
 __all__ = [
         'radialGaussian',
         'thinPlateSpline',
         'multiQuadric',
         'rbbases',
         'polybases',
         'polyfitname',
         'polydomcoeff',
-        'radialFunction',
-        'polyder',
         'polyval',
         'rbvander',
         'polyvander',
-        'homogeneizedpolyvander'
+        'homogeneizedpolyvander',
+        'RadialBasisInterpolator'
            ]
 
 
diff --git a/rrompy/utilities/poly_fitting/radial_basis/base.py b/rrompy/utilities/poly_fitting/radial_basis/base.py
index d5ace7d..130c979 100644
--- a/rrompy/utilities/poly_fitting/radial_basis/base.py
+++ b/rrompy/utilities/poly_fitting/radial_basis/base.py
@@ -1,55 +1,46 @@
 # 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 itertools import product
 from rrompy.utilities.base.types import Np1D, Np2D, paramList
 from rrompy.utilities.exception_manager import RROMPyException
 from rrompy.utilities.poly_fitting.polynomial.base import polydomcoeff as \
                                                                   polydomcoeffB
 
-__all__ = ['rbbases', 'polybases', 'polyfitname', 'polydomcoeff',
-           'radialFunction']
+__all__ = ['rbbases', 'polybases', 'polyfitname', 'polydomcoeff']
 
 rbbases = ["GAUSSIAN", "THINPLATE", "MULTIQUADRIC"]
 
 polybases = [x + "_" + y for x, y in product(
                                ["CHEBYSHEV", "LEGENDRE", "MONOMIAL"], rbbases)]
 
 def polyfitname(basis:str) -> str:
     fitpnames = {"CHEBYSHEV" : "chebfit", "LEGENDRE" : "legfit",
                  "MONOMIAL" : "polyfit"}
     fitrnames = {"GAUSSIAN" : "gaussian", "THINPLATE" : "thinplate",
                  "MULTIQUADRIC" : "multiquadric"}
     basisp, basisr = basis.split("_")
     try:
         return fitpnames[basisp] + "_" + fitrnames[basisr]
     except:
         raise RROMPyException("Polynomial-radial basis combination not "
                               "recognized.")
 
 def polydomcoeff(n:int, basis:str) -> float:
     return polydomcoeffB(n, basis.split("_")[0])
 
-class radialFunction:
-    def __init__(self, supportPoints : paramList = None,
-                 directionalWeights : Np1D = None, localCoeffs : Np1D = None,
-                 globalCoeffs : Np2D = None):
-        self.supportPoints = supportPoints
-        self.directionalWeights = directionalWeights
-        self.localCoeffs = localCoeffs
-        self.globalCoeffs = globalCoeffs
diff --git a/rrompy/utilities/poly_fitting/radial_basis/radial_basis_interpolator.py b/rrompy/utilities/poly_fitting/radial_basis/radial_basis_interpolator.py
new file mode 100644
index 0000000..98f9b62
--- /dev/null
+++ b/rrompy/utilities/poly_fitting/radial_basis/radial_basis_interpolator.py
@@ -0,0 +1,123 @@
+# Copyright (C) 2018 by the RROMPy authors
+#
+# This file is part of RROMPy.
+#
+# RROMPy is free software: you can redistribute it and/or modify
+# it under the terms of the GNU Lesser General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# RROMPy is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public License
+# along with RROMPy. If not, see <http://www.gnu.org/licenses/>.
+#
+
+import numpy as np
+from copy import deepcopy as copy
+from rrompy.utilities.base.types import List, ListAny, DictAny, Np1D, paramList
+from rrompy.utilities.poly_fitting.interpolator import GenericInterpolator
+from rrompy.utilities.poly_fitting.custom_fit import customFit
+from rrompy.utilities.poly_fitting.radial_basis.base import polyfitname
+from rrompy.utilities.poly_fitting.radial_basis.val import polyval
+from rrompy.utilities.poly_fitting.radial_basis.homogeneization import (
+                                                 homogeneizedpolyvander as hpv)
+from rrompy.utilities.poly_fitting.radial_basis.vander import polyvander as pv
+from rrompy.utilities.poly_fitting.polynomial.homogeneization import (
+                                                            homogeneizedToFull)
+from rrompy.utilities.exception_manager import RROMPyException
+from rrompy.parameter import checkParameterList
+
+__all__ = ['RadialBasisInterpolator']
+
+class RadialBasisInterpolator(GenericInterpolator):
+    """HERE"""
+
+    def __init__(self, other = None):
+        if other is None: return
+        self.support = other.support
+        self.coeffsGlobal = other.coeffsGlobal
+        self.coeffsLocal = other.coeffsLocal
+        self.directionalWeights = other.directionalWeights
+        self.npar = other.npar
+        self.polybasis = other.polybasis
+
+    @property
+    def shape(self):
+        if self.coeffsLocal.ndim > 1:
+            sh = self.coeffsLocal.shape[1 :]
+        else: sh = 1
+        return sh
+
+    @property
+    def deg(self):
+        return [x - 1 for x in self.coeffsGlobal.shape[: self.npar]]
+
+    def __call__(self, mu:paramList, der : List[int] = None,
+                 scl : Np1D = None):
+        if der is not None and np.sum(der) > 0:
+            raise RROMPyException(("Cannot take derivatives of radial basis "
+                                   "function."))           
+        return polyval(mu, self.coeffsGlobal, self.coeffsLocal, self.support,
+                       self.directionalWeights, self.polybasis)
+        
+    def __copy__(self):
+        return RadialBasisInterpolator(self)
+
+    def __deepcopy__(self, memo):
+        other = RadialBasisInterpolator()
+        (other.support, other.coeffsGlobal, other.coeffsLocal,
+         other.directionalWeights, other.npar, other.polybasis) = copy(
+                    (self.support, self.coeffsGlobal, self.coeffsLocal,
+                     self.directionalWeights, self.npar, self.polybasis), memo)
+        return other
+
+    def setupByInterpolation(self, support:paramList, values:ListAny,
+                             deg:int, polybasis : str = "MONOMIAL_GAUSSIAN",
+                             directionalWeights : Np1D = None,
+                             verbose : bool = True, homogeneized : bool = True,
+                             vanderCoeffs : DictAny = {},
+                             fitCoeffs : DictAny = {}):
+        support, _ = checkParameterList(support)
+        self.support = copy(support)
+        self.npar = support.shape[1]
+        if directionalWeights is None:
+            directionalWeights = np.ones(self.npar)
+        self.directionalWeights = directionalWeights
+        self.polybasis = polybasis
+        if homogeneized:
+            vander, _, reorder = hpv(support, deg = deg, basis = polybasis,
+                                  directionalWeights = self.directionalWeights,
+                                  **vanderCoeffs)
+            vander = vander[reorder]
+            vander = vander[:, reorder]
+        else:
+            if not hasattr(deg, "__len__"): deg = [deg] * self.npar
+            vander = pv(support, deg = deg, basis = polybasis,
+                        **vanderCoeffs)
+        outDim = values.shape[1:]
+        values = values.reshape(values.shape[0], -1)
+        values = np.pad(values, ((0, len(vander) - len(values)), (0, 0)),
+                        "constant")
+        fitOut = customFit(vander, values, full = True, **fitCoeffs)
+        P = fitOut[0][len(support) :]
+        if verbose:
+            msg = ("Fitting {}+{} samples with degree {} through {}... "
+                   "Conditioning of LS system: {:.4e}.").format(
+                                            len(support),
+                                            len(vander) - len(support), deg,
+                                            polyfitname(self.polybasis),
+                                            fitOut[1][2][0] / fitOut[1][2][-1])
+        else: msg = None
+        self.coeffsLocal = fitOut[0][: len(support)]
+        if homogeneized:
+            self.coeffsGlobal = homogeneizedToFull(
+                                         tuple([deg + 1] * self.npar) + outDim,
+                                         self.npar, P)
+        else:
+            self.coeffsGlobal = P.reshape(tuple([d + 1 for d in deg]) + outDim)
+        return fitOut[1][1] == vander.shape[1], msg
+
diff --git a/rrompy/utilities/poly_fitting/radial_basis/val.py b/rrompy/utilities/poly_fitting/radial_basis/val.py
index 51f4c5b..81d07f9 100644
--- a/rrompy/utilities/poly_fitting/radial_basis/val.py
+++ b/rrompy/utilities/poly_fitting/radial_basis/val.py
@@ -1,54 +1,52 @@
 # Copyright (C) 2018 by the RROMPy authors
 #
 # This file is part of RROMPy.
 #
 # RROMPy is free software: you can redistribute it and/or modify
 # it under the terms of the GNU Lesser General Public License as published by
 # the Free Software Foundation, either version 3 of the License, or
 # (at your option) any later version.
 #
 # RROMPy is distributed in the hope that it will be useful,
 # but WITHOUT ANY WARRANTY; without even the implied warranty of
 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
 # GNU Lesser General Public License for more details.
 #
 # You should have received a copy of the GNU Lesser General Public License
 # along with RROMPy. If not, see <http://www.gnu.org/licenses/>.
 #
 
 from copy import deepcopy as copy
 import numpy as np
-from rrompy.utilities.base.types import Np1D, Np2D, List, paramList, radialFun
+from rrompy.utilities.base.types import Np1D, Np2D, paramList
 from rrompy.parameter import checkParameterList
 from rrompy.utilities.exception_manager import RROMPyException
 from rrompy.utilities.poly_fitting.polynomial import polyval as pvP
-from .der import polyder
 from .kernel import radialGaussian, thinPlateSpline, multiQuadric
 
 __all__ = ['polyval']
 
-def polyval(x:paramList, c:radialFun, basis:str, m : List[int] = None,
-            scl : Np1D = None) -> Np2D:
-    cFun = polyder(c, basis, m = m, scl = scl)
+def polyval(x:paramList, cG:Np2D, cL:Np2D, supportPoints:paramList,
+            directionalWeights:Np1D, basis:str) -> Np2D:
     x, _ = checkParameterList(x)
     basisp, basisr = basis.split("_")
-    c = pvP(x, cFun.globalCoeffs, basisp, m, scl)
+    c = pvP(x, cG, basisp)
     try:
         radialvalbase = {"GAUSSIAN" : radialGaussian,
                          "THINPLATE" : thinPlateSpline,
                          "MULTIQUADRIC" : multiQuadric
                       }[basisr.upper()]
     except:
         raise RROMPyException("Radial basis not recognized.")
     csh = copy(c.shape)
     if len(csh) == 1: c = c.reshape(1, -1)
     for j in range(len(x)):
-        muDiff = (cFun.supportPoints.data - x[j]) * cFun.directionalWeights
+        muDiff = (supportPoints.data - x[j]) * directionalWeights
         r2j = np.sum(np.abs(muDiff) ** 2., axis = 1).reshape(1, -1)
-        val = radialvalbase(r2j).dot(cFun.localCoeffs)
+        val = radialvalbase(r2j).dot(cL)
         try:
             c[..., j] += val
         except:
             c[..., j] += val.flatten()
     if len(csh) == 1: c = c.flatten()
     return c
diff --git a/tests/utilities/radial_fitting.py b/tests/utilities/radial_fitting.py
index 766e0a9..d10448b 100644
--- a/tests/utilities/radial_fitting.py
+++ b/tests/utilities/radial_fitting.py
@@ -1,161 +1,167 @@
 # Copyright (C) 2018 by the RROMPy authors
 #
 # This file is part of RROMPy.
 #
 # RROMPy is free software: you can redistribute it and/or modify
 # it under the terms of the GNU Lesser General Public License as published by
 # the Free Software Foundation, either version 3 of the License, or
 # (at your option) any later version.
 #
 # RROMPy is distributed in the hope that it will be useful,
 # but WITHOUT ANY WARRANTY; without even the implied warranty of
 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
 # GNU Lesser General Public License for more details.
 #
 # You should have received a copy of the GNU Lesser General Public License
 # along with RROMPy. If not, see <http://www.gnu.org/licenses/>.
 #
 
 import numpy as np
 from rrompy.utilities.poly_fitting import customFit
 from rrompy.utilities.poly_fitting.radial_basis import (radialGaussian,
                                                         thinPlateSpline,
                                                         multiQuadric,
                                                         polybases, polyfitname,
-                                                        polydomcoeff, 
-                                                        radialFunction,
+                                                        polydomcoeff,
                                                         polyval, polyvander,
                                                         homogeneizedpolyvander)
 from rrompy.utilities.poly_fitting.polynomial import homogeneizedToFull
 from rrompy.parameter import checkParameterList
 
 def test_monomial_gaussian():
     polyrbname = "MONOMIAL_GAUSSIAN"
     assert polyrbname in polybases
     fitname = polyfitname(polyrbname)
     domcoeff = polydomcoeff(5, polyrbname)
     assert fitname == "polyfit_gaussian"
     assert np.isclose(domcoeff, 1., rtol = 1e-5)
 
     directionalWeights = np.array([5.])
     xSupp = checkParameterList(np.arange(-1, 3), 1)[0]
     cRBCoeffs = np.array([-1., 3., -3., 1., 1., 2., -.5])
-    cRB = radialFunction(directionalWeights = directionalWeights,
-                         supportPoints = xSupp, localCoeffs = cRBCoeffs[: 4],
-                         globalCoeffs = cRBCoeffs[4 :])
+    
+    globalCoeffs = cRBCoeffs[4 :]
+    localCoeffs = cRBCoeffs[: 4]
+
     ySupp = 1 + 2. * xSupp.data - .5 * xSupp.data ** 2.
     xx = np.linspace(-2., 3., 100)
-    yy = polyval(checkParameterList(xx, 1)[0], cRB, polyrbname)
+    yy = polyval(checkParameterList(xx, 1)[0], globalCoeffs, localCoeffs,
+                 xSupp, directionalWeights, polyrbname)
     yyman = 1. + 2. * xx - .5 * xx ** 2.
     for j, xc in enumerate(np.arange(-1, 3)):
         r2j = (5. * (xx - xc)) ** 2.
         rbj = radialGaussian(r2j)
         assert np.allclose(rbj, np.exp(-.5 * r2j))
-        yyman += cRB.localCoeffs[j] * rbj
-        ySupp += cRB.localCoeffs[j] * radialGaussian((directionalWeights[0]
-                                                    * (xSupp.data - xc)) ** 2.)
+        yyman += localCoeffs[j] * rbj
+        ySupp += localCoeffs[j] * radialGaussian((directionalWeights[0]
+                                               * (xSupp.data - xc)) ** 2.)
     assert np.allclose(yy, yyman, atol = 1e-5)
     
     VanT = polyvander(xSupp, [2], polyrbname,
                       directionalWeights = directionalWeights)
     ySupp = np.pad(ySupp.flatten(), (0, len(VanT) - len(xSupp)), "constant")
     out = customFit(VanT, ySupp)
     assert np.allclose(out, cRBCoeffs, atol = 1e-8)
 
 def test_legendre_thinplate():
     polyrbname = "LEGENDRE_THINPLATE"
     assert polyrbname in polybases
     fitname = polyfitname(polyrbname)
     domcoeff = polydomcoeff(5, polyrbname)
     assert fitname == "legfit_thinplate"
     assert np.isclose(domcoeff, 63. / 8, rtol = 1e-5)
 
     directionalWeights = np.array([.5])
     xSupp = checkParameterList(np.arange(-1, 3), 1)[0]
     cRBCoeffs = np.array([-1., 3., -3., 1., 1., 2., -.5])
-    cRB = radialFunction(directionalWeights = directionalWeights,
-                         supportPoints = xSupp, localCoeffs = cRBCoeffs[: 4],
-                         globalCoeffs = cRBCoeffs[4 :])
+    
+    localCoeffs = cRBCoeffs[: 4]
+    globalCoeffs = cRBCoeffs[4 :]
+    
     ySupp = 1 + 2. * xSupp.data - .5 * (.5 * (3. * xSupp.data ** 2. - 1.))
     xx = np.linspace(-2., 3., 100)
-    yy = polyval(checkParameterList(xx, 1)[0], cRB, polyrbname)
+    yy = polyval(checkParameterList(xx, 1)[0], globalCoeffs, localCoeffs,
+                 xSupp, directionalWeights, polyrbname)
     yyman = 1. + 2. * xx - .5 * (.5 * (3. * xx ** 2. - 1.))
     for j, xc in enumerate(np.arange(-1, 3)):
         r2j = (directionalWeights[0] * (xx - xc)) ** 2.
         rbj = thinPlateSpline(r2j)
         assert np.allclose(rbj, .5 * r2j * np.log(np.finfo(float).eps + r2j))
-        yyman += cRB.localCoeffs[j] * rbj
-        ySupp += cRB.localCoeffs[j] * thinPlateSpline((directionalWeights[0]
+        yyman += localCoeffs[j] * rbj
+        ySupp += localCoeffs[j] * thinPlateSpline((directionalWeights[0]
                                                     * (xSupp.data - xc)) ** 2.)
     assert np.allclose(yy, yyman, atol = 1e-5)
 
     VanT = polyvander(xSupp, [2], polyrbname,
                       directionalWeights = directionalWeights)
     ySupp = np.pad(ySupp.flatten(), (0, len(VanT) - len(xSupp)), "constant")
     out = customFit(VanT, ySupp)
     assert np.allclose(out, cRBCoeffs, atol = 1e-8)
 
 def test_chebyshev_multiquadric():
     polyrbname = "CHEBYSHEV_MULTIQUADRIC"
     assert polyrbname in polybases
     fitname = polyfitname(polyrbname)
     domcoeff = polydomcoeff(5, polyrbname)
     assert fitname == "chebfit_multiquadric"
     assert np.isclose(domcoeff, 16, rtol = 1e-5)
 
     directionalWeights = np.array([1.])
     xSupp = checkParameterList(np.arange(-1, 3), 1)[0]
     cRBCoeffs = np.array([-1., 3., -3., 1., 1., 2., -.5])
-    cRB = radialFunction(directionalWeights = directionalWeights,
-                         supportPoints = xSupp, localCoeffs = cRBCoeffs[: 4],
-                         globalCoeffs = cRBCoeffs[4 :])
+    
+    localCoeffs = cRBCoeffs[: 4]
+    globalCoeffs = cRBCoeffs[4 :]
+    
     ySupp = 1 + 2. * xSupp.data - .5 * (2. * xSupp.data ** 2. - 1.)
     xx = np.linspace(-2., 3., 100)
-    yy = polyval(checkParameterList(xx, 1)[0], cRB, polyrbname)
+    yy = polyval(checkParameterList(xx, 1)[0], globalCoeffs, localCoeffs,
+                 xSupp, directionalWeights, polyrbname)
     yyman = 1. + 2. * xx - .5 * (2. * xx ** 2. - 1.)
     for j, xc in enumerate(np.arange(-1, 3)):
         r2j = (directionalWeights[0] * (xx - xc)) ** 2.
         rbj = multiQuadric(r2j)
         assert np.allclose(rbj, np.power(r2j + 1, -.5))
-        yyman += cRB.localCoeffs[j] * rbj
-        ySupp += cRB.localCoeffs[j] * multiQuadric((directionalWeights[0]
-                                                  * (xSupp.data - xc)) ** 2.)
+        yyman += localCoeffs[j] * rbj
+        ySupp += localCoeffs[j] * multiQuadric((directionalWeights[0]
+                                             * (xSupp.data - xc)) ** 2.)
     assert np.allclose(yy, yyman, atol = 1e-5)
     
     VanT = polyvander(xSupp, [2], polyrbname,
                       directionalWeights = directionalWeights)
     ySupp = np.pad(ySupp.flatten(), (0, len(VanT) - len(xSupp)), "constant")
     out = customFit(VanT, ySupp)
     assert np.allclose(out, cRBCoeffs, atol = 1e-8)
 
 def test_total_degree_2d():
     values = lambda x, y: (x - 3.) ** 2. * y - (x + 1.) * y ** 2.
 
     polyrbname = "CHEBYSHEV_GAUSSIAN"
 
     xs, ys = np.meshgrid(np.linspace(0., 4., 5), np.linspace(0., 4., 4))
     xySupp = np.concatenate((xs.reshape(-1, 1), ys.reshape(-1, 1)), axis = 1)
     zs = values(xs, ys)
     zSupp = zs.flatten()
 
     deg = 3
     directionalWeights = [2., 1.]
     VanT, _, reidxs = homogeneizedpolyvander(xySupp, deg, polyrbname,
                                        directionalWeights = directionalWeights)
     VanT = VanT[reidxs]
     VanT = VanT[:, reidxs]
     cFit = np.linalg.solve(VanT, np.pad(zSupp, (0, len(VanT) - len(zSupp)),
                                         'constant'))
     globCoeff = homogeneizedToFull([deg + 1] * 2, 2, cFit[len(zSupp) :])
-    cRB = radialFunction(directionalWeights = directionalWeights,
-                         supportPoints = xySupp,
-                         localCoeffs = cFit[: len(zSupp)],
-                         globalCoeffs = globCoeff)
+
+    localCoeffs = cFit[: len(zSupp)]
+    globalCoeffs = globCoeff
+
     xx, yy = np.meshgrid(np.linspace(0., 4., 100), np.linspace(0., 4., 100))
     xxyy = np.concatenate((xx.reshape(-1, 1), yy.reshape(-1, 1)), axis = 1)
-    zz = polyval(xxyy, cRB, polyrbname).reshape(xx.shape)
+    zz = polyval(xxyy, globalCoeffs, localCoeffs, xySupp, directionalWeights,
+                 polyrbname).reshape(xx.shape)
     zzex = values(xx, yy)
 
     error = np.abs(zz - zzex)
     print(np.max(error))
     assert np.max(error) < 1e-10