diff --git a/rrompy/reduction_methods/__init__.py b/rrompy/reduction_methods/__init__.py
index f91242d..9e3d83f 100644
--- a/rrompy/reduction_methods/__init__.py
+++ b/rrompy/reduction_methods/__init__.py
@@ -1,49 +1,48 @@
 # 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 .standard import (NearestNeighbor, RationalInterpolant,
-                       RationalMovingLeastSquares, RationalPade, ReducedBasis)
+from .standard import (NearestNeighbor, RationalInterpolant, RationalPade,
+                       ReducedBasis)
 from .standard.greedy import RationalInterpolantGreedy, ReducedBasisGreedy
 from .pivoted import (RationalInterpolantPivotedNoMatch,
                       RationalInterpolantPivoted,
                       RationalInterpolantGreedyPivotedNoMatch,
                       RationalInterpolantGreedyPivoted)
 from .pivoted.greedy import (RationalInterpolantPivotedGreedyNoMatch,
                              RationalInterpolantPivotedGreedy,
                              RationalInterpolantGreedyPivotedGreedyNoMatch,
                              RationalInterpolantGreedyPivotedGreedy)
 
 __all__ = [
         'NearestNeighbor',
         'RationalInterpolant',
-        'RationalMovingLeastSquares',
         'RationalPade',
         'ReducedBasis',
         'RationalInterpolantGreedy',
         'ReducedBasisGreedy',
         'RationalInterpolantPivotedNoMatch',
         'RationalInterpolantPivoted',
         'RationalInterpolantGreedyPivotedNoMatch',
         'RationalInterpolantGreedyPivoted',
         'RationalInterpolantPivotedGreedyNoMatch',
         'RationalInterpolantPivotedGreedy',
         'RationalInterpolantGreedyPivotedGreedyNoMatch',
         'RationalInterpolantGreedyPivotedGreedy'
           ]
 
 
diff --git a/rrompy/reduction_methods/pivoted/generic_pivoted_approximant.py b/rrompy/reduction_methods/pivoted/generic_pivoted_approximant.py
index 3a539d8..c6aa142 100644
--- a/rrompy/reduction_methods/pivoted/generic_pivoted_approximant.py
+++ b/rrompy/reduction_methods/pivoted/generic_pivoted_approximant.py
@@ -1,579 +1,577 @@
 # 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.reduction_methods.base.generic_approximant import (
                                                             GenericApproximant)
 from rrompy.sampling import (SamplingEnginePivoted, SamplingEnginePivotedPOD,
                              SamplingEnginePivotedPODGlobal)
 from rrompy.utilities.poly_fitting.polynomial import polybases as ppb
 from rrompy.utilities.poly_fitting.radial_basis import polybases as rbpb
-from rrompy.utilities.poly_fitting.moving_least_squares import (
-                                                            polybases as mlspb)
 from rrompy.utilities.base.types import paramList, ListAny
 from rrompy.utilities.base import verbosityManager as vbMng
 from rrompy.utilities.numerical.degree import reduceDegreeN
 from rrompy.utilities.exception_manager import (RROMPyException, RROMPyAssert,
                                                 RROMPyWarning)
 
 __all__ = ['GenericPivotedApproximantNoMatch', 'GenericPivotedApproximant',
            'PODGlobal']
 
 PODGlobal = 2
 
 class GenericPivotedApproximantBase(GenericApproximant):
     def __init__(self, directionPivot:ListAny, *args, **kwargs):
         self._preInit()
         if len(directionPivot) > 1:
             raise RROMPyException(("Exactly 1 pivot parameter allowed in pole "
                                    "matching."))
         from rrompy.parameter.parameter_sampling import QuadratureSampler as QS
         QSBase = QS([[0], [1]], "UNIFORM")
         self._addParametersToList(["cutOffTolerance",
                                    "radialDirectionalWeightsMarginal",
                                    "interpRcondMarginal"], [np.inf, [1.], -1],
                                   ["samplerPivot", "SMarginal",
                                    "samplerMarginal"], [QSBase, [1], QSBase])
         del QS
         self._directionPivot = directionPivot
         super().__init__(*args, **kwargs)
         self._postInit()
 
     def setupSampling(self):
         """Setup sampling engine."""
         RROMPyAssert(self._mode, message = "Cannot setup sampling engine.")
         if not hasattr(self, "_POD") or self._POD is None: return
         if self.POD:
             if self.POD == PODGlobal:
                 SamplingEngine = SamplingEnginePivotedPODGlobal
             else:
                 SamplingEngine = SamplingEnginePivotedPOD
         else:
             SamplingEngine = SamplingEnginePivoted
         self.samplingEngine = SamplingEngine(self.HFEngine,
                                              self.directionPivot,
                                              sample_state = self.approx_state,
                                              verbosity = self.verbosity)
 
     def initializeModelData(self, datadict):
         if "directionPivot" in datadict.keys():
             from .trained_model.trained_model_pivoted_data import (
                                                        TrainedModelPivotedData)
             return (TrainedModelPivotedData(datadict["mu0"],
                                             datadict.pop("projMat"),
                                             datadict["scaleFactor"],
                                             datadict.pop("rescalingExp"),
                                             datadict["directionPivot"]),
                     ["mu0", "scaleFactor", "directionPivot", "mus"])
         else:
             return super().initializeModelData(datadict)
 
     @property
     def npar(self):
         """Number of parameters."""
         if hasattr(self, "_temporaryPivot"): return self.nparPivot
         return super().npar
 
     @property
     def mus(self):
         """Value of mus. Its assignment may reset snapshots."""
         return self._mus
     @mus.setter
     def mus(self, mus):
         musOld =  copy(self.mus) if hasattr(self, '_mus') else None
         if (musOld is None or len(mus) != len(musOld) or not mus == musOld):
             self.resetSamples()
             self._mus = mus
 
     @property
     def cutOffTolerance(self):
         """Value of cutOffTolerance."""
         return self._cutOffTolerance
     @cutOffTolerance.setter
     def cutOffTolerance(self, cutOffTolerance):
         self._cutOffTolerance = cutOffTolerance
         self._approxParameters["cutOffTolerance"] = self.cutOffTolerance
 
     @property
     def SMarginal(self):
         """Value of SMarginal."""
         return self._SMarginal
     @SMarginal.setter
     def SMarginal(self, SMarginal):
         if SMarginal <= 0:
             raise RROMPyException("SMarginal must be positive.")
         if hasattr(self, "_SMarginal") and self._SMarginal is not None:
             Sold = self.SMarginal
         else: Sold = -1
         self._SMarginal = SMarginal
         self._approxParameters["SMarginal"] = self.SMarginal
         if Sold != self.SMarginal: self.resetSamples()
 
     @property
     def radialDirectionalWeightsMarginal(self):
         """Value of radialDirectionalWeightsMarginal."""
         return self._radialDirectionalWeightsMarginal
     @radialDirectionalWeightsMarginal.setter
     def radialDirectionalWeightsMarginal(self, radialDirWeightsMarginal):
         if not hasattr(radialDirWeightsMarginal, "__len__"):
             radialDirWeightsMarginal = [radialDirWeightsMarginal]
         self._radialDirectionalWeightsMarginal = radialDirWeightsMarginal
         self._approxParameters["radialDirectionalWeightsMarginal"] = (
                                          self.radialDirectionalWeightsMarginal)
 
     @property
     def interpRcondMarginal(self):
         """Value of interpRcondMarginal."""
         return self._interpRcondMarginal
     @interpRcondMarginal.setter
     def interpRcondMarginal(self, interpRcondMarginal):
         self._interpRcondMarginal = interpRcondMarginal
         self._approxParameters["interpRcondMarginal"] = (
                                                       self.interpRcondMarginal)
 
     @property
     def directionPivot(self):
         """Value of directionPivot. Its assignment may reset snapshots."""
         return self._directionPivot
     @directionPivot.setter
     def directionPivot(self, directionPivot):
         if hasattr(self, '_directionPivot'):
             directionPivotOld = copy(self.directionPivot)
         else:
             directionPivotOld = None
         if (directionPivotOld is None
          or len(directionPivot) != len(directionPivotOld)
          or not directionPivot == directionPivotOld):
             self.resetSamples()
             self._directionPivot = directionPivot
 
     @property
     def directionMarginal(self):
         return [x for x in range(self.HFEngine.npar) \
                                                if x not in self.directionPivot]
 
     @property
     def nparPivot(self):
         return len(self.directionPivot)
 
     @property
     def nparMarginal(self):
         return self.npar - self.nparPivot
 
     @property
     def rescalingExpPivot(self):
         return [self.HFEngine.rescalingExp[x] for x in self.directionPivot]
 
     @property
     def rescalingExpMarginal(self):
         return [self.HFEngine.rescalingExp[x] for x in self.directionMarginal]
 
     @property
     def muBounds(self):
         """Value of muBounds."""
         return self.samplerPivot.lims
 
     @property
     def muBoundsMarginal(self):
         """Value of muBoundsMarginal."""
         return self.samplerMarginal.lims
 
     @property
     def samplerPivot(self):
         """Value of samplerPivot."""
         return self._samplerPivot
     @samplerPivot.setter
     def samplerPivot(self, samplerPivot):
         if 'generatePoints' not in dir(samplerPivot):
             raise RROMPyException("Pivot sampler type not recognized.")
         if hasattr(self, '_samplerPivot') and self._samplerPivot is not None:
             samplerOld = self.samplerPivot
         self._samplerPivot = samplerPivot
         self._approxParameters["samplerPivot"] = self.samplerPivot
         if not 'samplerOld' in locals() or samplerOld != self.samplerPivot:
             self.resetSamples()
     
     @property
     def samplerMarginal(self):
         """Value of samplerMarginal."""
         return self._samplerMarginal
     @samplerMarginal.setter
     def samplerMarginal(self, samplerMarginal):
         if 'generatePoints' not in dir(samplerMarginal):
             raise RROMPyException("Marginal sampler type not recognized.")
         if (hasattr(self, '_samplerMarginal')
         and self._samplerMarginal is not None):
             samplerOld = self.samplerMarginal
         self._samplerMarginal = samplerMarginal
         self._approxParameters["samplerMarginal"] = self.samplerMarginal
         if not 'samplerOld' in locals() or samplerOld != self.samplerMarginal:
             self.resetSamples()
     
     def setSamples(self, samplingEngine):
         """Copy samplingEngine and samples."""
         self.mus = copy(samplingEngine.mus[0])
         for sEj in samplingEngine.mus[1:]:
             self.mus.append(sEj)
         super().setSamples(samplingEngine)
 
     def computeScaleFactor(self):
         """Compute parameter rescaling factor."""
         RROMPyAssert(self._mode, message = "Cannot compute rescaling factor.")
         self.scaleFactorPivot = .5 * np.abs(
                                     self.muBounds[0] ** self.rescalingExpPivot
                                   - self.muBounds[1] ** self.rescalingExpPivot)
         self.scaleFactorMarginal = .5 * np.abs(
                          self.muBoundsMarginal[0] ** self.rescalingExpMarginal
                        - self.muBoundsMarginal[1] ** self.rescalingExpMarginal)
         self.scaleFactor = np.empty(self.npar)
         self.scaleFactor[self.directionPivot] = self.scaleFactorPivot
         self.scaleFactor[self.directionMarginal] = self.scaleFactorMarginal
 
     def normApprox(self, mu:paramList) -> float:
         _PODOld = self.POD
         self._POD = self.POD == PODGlobal
         result = super().normApprox(mu)
         self._POD = _PODOld
         return result
 
 class GenericPivotedApproximantNoMatch(GenericPivotedApproximantBase):
     """
     ROM pivoted approximant (without pole matching) computation for parametric
         problems (ABSTRACT).
 
     Args:
         HFEngine: HF problem solver.
         mu0(optional): Default parameter. Defaults to 0.
         directionPivot(optional): Pivot components. Defaults to [0].
         approxParameters(optional): Dictionary containing values for main
             parameters of approximant. Recognized keys are:
             - 'POD': whether to compute POD of snapshots; defaults to True;
             - 'scaleFactorDer': scaling factors for derivative computation;
                 defaults to 'AUTO';
             - 'cutOffTolerance': tolerance for ignoring parasitic poles;
                 defaults to np.inf;
             - 'S': total number of pivot samples current approximant relies
                 upon;
             - 'samplerPivot': pivot sample point generator;
             - 'SMarginal': total number of marginal samples current approximant
                 relies upon;
             - 'samplerMarginal': marginal sample point generator;
             - 'radialDirectionalWeightsMarginal': radial basis weights for
                 marginal interpolant; defaults to 1;
             - 'interpRcondMarginal': tolerance for marginal interpolation;
                 defaults to None.
             Defaults to empty dict.
         approx_state(optional): Whether to approximate state. Defaults to
             False.
         verbosity(optional): Verbosity level. Defaults to 10.
 
     Attributes:
         HFEngine: HF problem solver.
         mu0: Default parameter.
         directionPivot: Pivot components.
         mus: Array of snapshot parameters.
         musMarginal: Array of marginal snapshot parameters.
         approxParameters: Dictionary containing values for main parameters of
             approximant. Recognized keys are in parameterList.
         parameterListSoft: Recognized keys of soft approximant parameters:
             - 'POD': whether to compute POD of snapshots;
             - 'scaleFactorDer': scaling factors for derivative computation;
             - 'cutOffTolerance': tolerance for ignoring parasitic poles;
             - 'radialDirectionalWeightsMarginal': radial basis weights for
                 marginal interpolant;
             - 'interpRcondMarginal': tolerance for marginal interpolation.
         parameterListCritical: Recognized keys of critical approximant
             parameters:
             - 'S': total number of pivot samples current approximant relies
                 upon;
             - 'samplerPivot': pivot sample point generator;
             - 'SMarginal': total number of marginal samples current approximant
                 relies upon;
             - 'samplerMarginal': marginal sample point generator.
         approx_state: Whether to approximate state.
         verbosity: Verbosity level.
         POD: Whether to compute POD of snapshots.
         scaleFactorDer: Scaling factors for derivative computation.
         cutOffTolerance: Tolerance for ignoring parasitic poles.
         S: Total number of pivot samples current approximant relies upon.
         samplerPivot: Pivot sample point generator.
         SMarginal: Total number of marginal samples current approximant relies
             upon.
         samplerMarginal: Marginal sample point generator.
         radialDirectionalWeightsMarginal: Radial basis weights for marginal
             interpolant.
         interpRcondMarginal: Tolerance for marginal interpolation.
         muBounds: list of bounds for pivot parameter values.
         muBoundsMarginal: list of bounds for marginal parameter values.
         samplingEngine: Sampling engine.
         uHF: High fidelity solution(s) with parameter(s) lastSolvedHF as
             sampleList.
         lastSolvedHF: Parameter(s) corresponding to last computed high fidelity
             solution(s) as parameterList.
         uApproxReduced: Reduced approximate solution(s) with parameter(s)
             lastSolvedApprox as sampleList.
         lastSolvedApproxReduced: Parameter(s) corresponding to last computed
             reduced approximate solution(s) as parameterList.
         uApprox: Approximate solution(s) with parameter(s) lastSolvedApprox as
             sampleList.
         lastSolvedApprox: Parameter(s) corresponding to last computed
             approximate solution(s) as parameterList.
     """
 
     @property
     def tModelType(self):
         from .trained_model.trained_model_pivoted_rational_nomatch import (
                                             TrainedModelPivotedRationalNoMatch)
         return TrainedModelPivotedRationalNoMatch
     def _finalizeMarginalization(self):
         vbMng(self, "INIT", "Recompressing by cut off.", 10)
         msg = self.trainedModel.recompressByCutOff(self.cutOffTolerance,
                                            self.samplerPivot.normalFoci(),
                                            self.samplerPivot.groundPotential())
         vbMng(self, "DEL", "Done recompressing." + msg, 10)
         self.trainedModel.setupMarginalInterp(self,
                                          self.radialDirectionalWeightsMarginal)
         self.trainedModel.data.approxParameters = copy(self.approxParameters)
 
 class GenericPivotedApproximant(GenericPivotedApproximantBase):
     """
     ROM pivoted approximant (with pole matching) computation for parametric
         problems (ABSTRACT).
 
     Args:
         HFEngine: HF problem solver.
         mu0(optional): Default parameter. Defaults to 0.
         directionPivot(optional): Pivot components. Defaults to [0].
         approxParameters(optional): Dictionary containing values for main
             parameters of approximant. Recognized keys are:
             - 'POD': whether to compute POD of snapshots; defaults to True;
             - 'scaleFactorDer': scaling factors for derivative computation;
                 defaults to 'AUTO';
             - 'matchingWeight': weight for pole matching optimization; defaults
                 to 1;
             - 'cutOffTolerance': tolerance for ignoring parasitic poles;
                 defaults to np.inf;
             - 'cutOffSharedRatio': required ratio of marginal points to share
                 resonance in cut off strategy; defaults to 1.;
             - 'S': total number of pivot samples current approximant relies
                 upon;
             - 'samplerPivot': pivot sample point generator;
             - 'SMarginal': total number of marginal samples current approximant
                 relies upon;
             - 'samplerMarginal': marginal sample point generator;
             - 'polybasisMarginal': type of polynomial basis for marginal
                 interpolation; allowed values include 'MONOMIAL', 'CHEBYSHEV'
                 and 'LEGENDRE'; defaults to 'MONOMIAL';
             - 'MMarginal': degree of marginal interpolant; defaults to 'AUTO',
                 i.e. maximum allowed;
             - 'polydegreetypeMarginal': type of polynomial degree for marginal;
                 defaults to 'TOTAL';
             - 'radialDirectionalWeightsMarginal': radial basis weights for
                 marginal interpolant; defaults to 1;
             - 'interpRcondMarginal': tolerance for marginal interpolation;
                 defaults to None.
             Defaults to empty dict.
         approx_state(optional): Whether to approximate state. Defaults to
             False.
         verbosity(optional): Verbosity level. Defaults to 10.
 
     Attributes:
         HFEngine: HF problem solver.
         mu0: Default parameter.
         directionPivot: Pivot components.
         mus: Array of snapshot parameters.
         musMarginal: Array of marginal snapshot parameters.
         approxParameters: Dictionary containing values for main parameters of
             approximant. Recognized keys are in parameterList.
         parameterListSoft: Recognized keys of soft approximant parameters:
             - 'POD': whether to compute POD of snapshots;
             - 'scaleFactorDer': scaling factors for derivative computation;
             - 'matchingWeight': weight for pole matching optimization;
             - 'cutOffTolerance': tolerance for ignoring parasitic poles;
             - 'cutOffSharedRatio': required ratio of marginal points to share
                 resonance in cut off strategy;
             - 'polybasisMarginal': type of polynomial basis for marginal
                 interpolation;
             - 'MMarginal': degree of marginal interpolant;
             - 'polydegreetypeMarginal': type of polynomial degree for marginal;
             - 'radialDirectionalWeightsMarginal': radial basis weights for
                 marginal interpolant;
             - 'interpRcondMarginal': tolerance for marginal interpolation.
         parameterListCritical: Recognized keys of critical approximant
             parameters:
             - 'S': total number of pivot samples current approximant relies
                 upon;
             - 'samplerPivot': pivot sample point generator;
             - 'SMarginal': total number of marginal samples current approximant
                 relies upon;
             - 'samplerMarginal': marginal sample point generator.
         approx_state: Whether to approximate state.
         verbosity: Verbosity level.
         POD: Whether to compute POD of snapshots.
         scaleFactorDer: Scaling factors for derivative computation.
         matchingWeight: Weight for pole matching optimization.
         cutOffTolerance: Tolerance for ignoring parasitic poles.
         cutOffSharedRatio: Required ratio of marginal points to share resonance
             in cut off strategy.
         S: Total number of pivot samples current approximant relies upon.
         samplerPivot: Pivot sample point generator.
         SMarginal: Total number of marginal samples current approximant relies
             upon.
         samplerMarginal: Marginal sample point generator.
         polybasisMarginal: Type of polynomial basis for marginal interpolation.
         MMarginal: Degree of marginal interpolant.
         polydegreetypeMarginal: Type of polynomial degree for marginal.
         radialDirectionalWeightsMarginal: Radial basis weights for marginal
             interpolant.
         interpRcondMarginal: Tolerance for marginal interpolation.
         muBounds: list of bounds for pivot parameter values.
         muBoundsMarginal: list of bounds for marginal parameter values.
         samplingEngine: Sampling engine.
         uHF: High fidelity solution(s) with parameter(s) lastSolvedHF as
             sampleList.
         lastSolvedHF: Parameter(s) corresponding to last computed high fidelity
             solution(s) as parameterList.
         uApproxReduced: Reduced approximate solution(s) with parameter(s)
             lastSolvedApprox as sampleList.
         lastSolvedApproxReduced: Parameter(s) corresponding to last computed
             reduced approximate solution(s) as parameterList.
         uApprox: Approximate solution(s) with parameter(s) lastSolvedApprox as
             sampleList.
         lastSolvedApprox: Parameter(s) corresponding to last computed
             approximate solution(s) as parameterList.
     """
 
     def __init__(self, *args, **kwargs):
         self._preInit()
         self._addParametersToList(["matchingWeight", "cutOffSharedRatio",
                                    "polybasisMarginal", "MMarginal",
                                    "polydegreetypeMarginal"],
                                   [1., 1., "MONOMIAL", "AUTO", "TOTAL"])
         super().__init__(*args, **kwargs)
         self._postInit()
 
     @property
     def tModelType(self):
         from .trained_model.trained_model_pivoted_rational import (
                                                    TrainedModelPivotedRational)
         return TrainedModelPivotedRational
 
     @property
     def matchingWeight(self):
         """Value of matchingWeight."""
         return self._matchingWeight
     @matchingWeight.setter
     def matchingWeight(self, matchingWeight):
         self._matchingWeight = matchingWeight
         self._approxParameters["matchingWeight"] = self.matchingWeight
 
     @property
     def cutOffSharedRatio(self):
         """Value of cutOffSharedRatio."""
         return self._cutOffSharedRatio
     @cutOffSharedRatio.setter
     def cutOffSharedRatio(self, cutOffSharedRatio):
         if cutOffSharedRatio > 1.:
             RROMPyWarning("Cut off shared ratio too large. Clipping to 1.")
             cutOffSharedRatio = 1.
         elif cutOffSharedRatio < 0.:
             RROMPyWarning("Cut off shared ratio too small. Clipping to 0.")
             cutOffSharedRatio = 0.
         self._cutOffSharedRatio = cutOffSharedRatio
         self._approxParameters["cutOffSharedRatio"] = self.cutOffSharedRatio
 
     @property
     def polybasisMarginal(self):
         """Value of polybasisMarginal."""
         return self._polybasisMarginal
     @polybasisMarginal.setter
     def polybasisMarginal(self, polybasisMarginal):
         try:
             polybasisMarginal = polybasisMarginal.upper().strip().replace(" ",
                                                                           "")
-            if polybasisMarginal not in ppb + rbpb + mlspb: 
+            if polybasisMarginal not in ppb + rbpb: 
                 raise RROMPyException(
                                "Prescribed marginal polybasis not recognized.")
             self._polybasisMarginal = polybasisMarginal
         except:
             RROMPyWarning(("Prescribed marginal polybasis not recognized. "
                            "Overriding to 'MONOMIAL'."))
             self._polybasisMarginal = "MONOMIAL"
         self._approxParameters["polybasisMarginal"] = self.polybasisMarginal
 
     @property
     def MMarginal(self):
         """Value of MMarginal."""
         return self._MMarginal
     @MMarginal.setter
     def MMarginal(self, MMarginal):
         if isinstance(MMarginal, str):
             MMarginal = MMarginal.strip().replace(" ","")
             if "-" not in MMarginal: MMarginal = MMarginal + "-0"
             self._MMarginal_isauto = True
             self._MMarginal_shift = int(MMarginal.split("-")[-1])
             MMarginal = 0
         if MMarginal < 0:
             raise RROMPyException("MMarginal must be non-negative.")
         self._MMarginal = MMarginal
         self._approxParameters["MMarginal"] = self.MMarginal
 
     def _setMMarginalAuto(self):
         self.MMarginal = max(0, reduceDegreeN(
                                  len(self.musMarginal), len(self.musMarginal),
                                  self.nparMarginal, self.polydegreetypeMarginal
                                              ) - self._MMarginal_shift)
         vbMng(self, "MAIN", ("Automatically setting MMarginal to "
                              "{}.").format(self.MMarginal), 25)
 
     @property
     def polydegreetypeMarginal(self):
         """Value of polydegreetypeMarginal."""
         return self._polydegreetypeMarginal
     @polydegreetypeMarginal.setter
     def polydegreetypeMarginal(self, polydegreetypeM):
         try:
             polydegreetypeM = polydegreetypeM.upper().strip().replace(" ","")
             if polydegreetypeM not in ["TOTAL", "FULL"]: 
                 raise RROMPyException(("Prescribed polydegreetypeMarginal not "
                                        "recognized."))
             self._polydegreetypeMarginal = polydegreetypeM
         except:
             RROMPyWarning(("Prescribed polydegreetypeMarginal not recognized. "
                            "Overriding to 'TOTAL'."))
             self._polydegreetypeMarginal = "TOTAL"
         self._approxParameters["polydegreetypeMarginal"] = (
                                                    self.polydegreetypeMarginal)
 
     def _finalizeMarginalization(self):
         vbMng(self, "INIT", "Recompressing by cut off.", 10)
         msg = self.trainedModel.recompressByCutOff(self.cutOffTolerance,
                                            self.cutOffSharedRatio,
                                            self.samplerPivot.normalFoci(),
                                            self.samplerPivot.groundPotential())
         vbMng(self, "DEL", "Done recompressing." + msg, 10)
         interpPars = [self.verbosity >= 5,
                       self.polydegreetypeMarginal == "TOTAL", {}]
         if self.polybasisMarginal in ppb + rbpb:
             interpPars += [{"rcond": self.interpRcondMarginal}]
         self.trainedModel.setupMarginalInterp(self, interpPars,
                                          hasattr(self, "_MMarginal_isauto"),
                                          self.radialDirectionalWeightsMarginal,
                                          hasattr(self, "_reduceDegreeNNoWarn"))
         self.trainedModel.data.approxParameters = copy(self.approxParameters)
diff --git a/rrompy/reduction_methods/pivoted/trained_model/trained_model_pivoted_rational.py b/rrompy/reduction_methods/pivoted/trained_model/trained_model_pivoted_rational.py
index ded2be9..7f877ae 100644
--- a/rrompy/reduction_methods/pivoted/trained_model/trained_model_pivoted_rational.py
+++ b/rrompy/reduction_methods/pivoted/trained_model/trained_model_pivoted_rational.py
@@ -1,282 +1,280 @@
 # 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 .trained_model_pivoted_rational_nomatch import (
                                             TrainedModelPivotedRationalNoMatch)
 from rrompy.utilities.base.types import (Np1D, Np2D, List, ListAny, paramVal,
                                          paramList, HFEng)
 from rrompy.utilities.base import verbosityManager as vbMng
 from rrompy.utilities.numerical.point_matching import (pointMatching,
                                               chordalMetricAdjusted, potential)
 from rrompy.utilities.numerical.degree import reduceDegreeN
 from rrompy.utilities.poly_fitting.polynomial import (polybases as ppb,
                                                   PolynomialInterpolator as PI)
-from rrompy.utilities.poly_fitting.radial_basis import (polybases as rbpb,
+from rrompy.utilities.poly_fitting.radial_basis import (
                                                 RadialBasisInterpolator as RBI)
-from rrompy.utilities.poly_fitting.moving_least_squares import (
-                                        MovingLeastSquaresInterpolator as MLSI)
 from rrompy.utilities.poly_fitting.heaviside import (heavisideUniformShape,
                                                    rational2heaviside,
                                                    HeavisideInterpolator as HI)
 from rrompy.utilities.exception_manager import (RROMPyException, RROMPyAssert,
                                                 RROMPyWarning)
 from rrompy.parameter import checkParameterList
 
 __all__ = ['TrainedModelPivotedRational']
 
 class TrainedModelPivotedRational(TrainedModelPivotedRationalNoMatch):
     """
     ROM approximant evaluation for pivoted approximants based on interpolation
         of rational approximants (with pole matching).
     
     Attributes:
         Data: dictionary with all that can be pickled.
     """
 
     def centerNormalizeMarginal(self, mu : paramList = [],
                                 mu0 : paramVal = None) -> paramList:
         """
         Compute normalized parameter to be plugged into approximant.
 
         Args:
             mu: Parameter(s) 1.
             mu0: Parameter(s) 2. If None, set to self.data.mu0Marginal.
 
         Returns:
             Normalized parameter.
         """
         mu = checkParameterList(mu, self.data.nparMarginal)[0]
         if mu0 is None: mu0 = self.data.mu0Marginal
         rad = ((mu ** self.data.rescalingExpMarginal
               - mu0 ** self.data.rescalingExpMarginal)
               / self.data.scaleFactorMarginal)
         return rad
 
     def setupMarginalInterp(self, approx, interpPars:ListAny,
                             MMAuto:bool, rDWM : Np1D = None,
                             noWarnReduceAuto : bool = True):
         vbMng(self, "INIT", "Starting computation of marginal interpolator.",
               12)
         musMCN = self.centerNormalizeMarginal(self.data.musMarginal)
         pbM = approx.polybasisMarginal
         self.data.marginalInterp = []
         for ipts, pts in enumerate(self.data.suppEffPts):
             if len(pts) > 0:
                 musMCNEff = musMCN[pts]
                 if MMAuto:
                     if ipts > 0:
                         verb = approx.verbosity
                         approx.verbosity = 0
                         _musM = approx.musMarginal
                         approx.musMarginal = musMCNEff
                     approx._setMMarginalAuto()
                     if ipts > 0:
                         approx.musMarginal = _musM
                         approx.verbosity = verb
                 if ipts == 0: _MMarginalEffective = approx.MMarginal
                 if not MMAuto:
                     approx.MMarginal = _MMarginalEffective
                     MM = reduceDegreeN(approx.MMarginal, len(musMCNEff),
                                        self.data.nparMarginal, 
                                        approx.polydegreetypeMarginal)
                     if MM < approx.MMarginal:
                         if ipts == 0 and not noWarnReduceAuto:
                             RROMPyWarning(("MMarginal too large compared to "
                                            "SMarginal. Reducing MMarginal by "
                                            "{}").format(approx.MMarginal - MM))
                         approx.MMarginal = MM
                 MMEff = approx.MMarginal
                 while MMEff >= 0:
                     pParRest = copy(interpPars)
                     if pbM in ppb:
                         p = PI()
                     else:
                         pParRest = [rDWM] + pParRest
-                        p = RBI() if pbM in rbpb else MLSI()
+                        p = RBI()
                     wellCond, msg = p.setupByInterpolation(musMCNEff,
                                                         np.eye(len(musMCNEff)),
                                                         MMEff, pbM, *pParRest)
                     vbMng(self, "MAIN", msg, 30)
                     if wellCond: break
                     vbMng(self, "MAIN",
                           ("Polyfit is poorly conditioned. Reducing "
                            "MMarginal by 1."), 35)
                     MMEff -= 1
                 if MMEff < 0:
                     raise RROMPyException(("Instability in computation of "
                                            "interpolant. Aborting."))
             self.data.marginalInterp += [copy(p)]
         approx.MMarginal = _MMarginalEffective
         vbMng(self, "DEL", "Done computing marginal interpolator.", 12)
 
     def initializeFromLists(self, poles:ListAny, coeffs:ListAny, basis:str,
                             HFEngine:HFEng, matchingWeight : float = 1.,
                             POD : bool = True, is_state : bool = True):
         """Initialize Heaviside representation."""
         musM = self.data.musMarginal
         margAbsDist = np.sum(np.abs(np.repeat(musM.data, len(musM), 0)
                                   - np.tile(musM.data, [len(musM), 1])
                                    ), axis = 1).reshape(len(musM), len(musM))
         explored = [0]
         unexplored = list(range(1, len(musM)))
         poles, coeffs = heavisideUniformShape(poles, coeffs)
         N = len(poles[0])
         for _ in range(1, len(musM)):
             minIdx = np.argmin(np.concatenate([margAbsDist[ex, unexplored] \
                                                           for ex in explored]))
             minIex = explored[minIdx // len(unexplored)]
             minIunex = unexplored[minIdx % len(unexplored)]
             resex = coeffs[minIex][: N]
             resunex = coeffs[minIunex][: N]
             if matchingWeight != 0 and not POD:
                 resex = self.data.projMat[:, : resex.shape[1]].dot(resex.T)
                 resunex = self.data.projMat[:, : resunex.shape[1]].dot(
                                                                      resunex.T)
             dist = chordalMetricAdjusted(poles[minIex], poles[minIunex],
                                          matchingWeight, resex, resunex,
                                          HFEngine, is_state)
             reordering = pointMatching(dist)[1]
             poles[minIunex] = poles[minIunex][reordering]
             coeffs[minIunex][: N] = coeffs[minIunex][reordering]
             explored += [minIunex]
             unexplored.remove(minIunex)
         super().initializeFromLists(poles, coeffs, basis)
         self.data.suppEffPts = [np.arange(len(self.data.HIs))]
         self.data.suppEffIdx = np.zeros(N, dtype = int)
 
     def initializeFromRational(self, HFEngine:HFEng,
                                matchingWeight : float = 1., POD : bool = True,
                                is_state : bool = True):
         """Initialize Heaviside representation."""
         RROMPyAssert(self.data.nparPivot, 1, "Number of pivot parameters")
         poles, coeffs = [], []
         for Q, P in zip(self.data.Qs, self.data.Ps):
             cfs, pls, basis = rational2heaviside(P, Q)
             poles += [pls]
             coeffs += [cfs]
         self.initializeFromLists(poles, coeffs, basis, HFEngine,
                                  matchingWeight, POD, is_state)
 
     def recompressByCutOff(self, tol:float, shared:float,
                            foci:List[np.complex], ground:float) -> str:
         N = len(self.data.HIs[0].poles)
         M = len(self.data.HIs)
         mu0 = np.mean(foci)
         goodLocPoles = np.array([potential(hi.poles, foci) - ground
                                       <= tol * ground for hi in self.data.HIs])
         self.data.suppEffPts = [np.arange(len(self.data.HIs))]
         self.data.suppEffIdx = np.zeros(N, dtype = int)
         if np.all(np.all(goodLocPoles)): return " No poles erased."
         goodGlobPoles = np.sum(goodLocPoles, axis = 0)
         goodEnoughPoles = goodGlobPoles >= max(1., 1. * shared * M)
         keepPole = np.where(goodEnoughPoles)[0]
         halfPole = np.where(np.logical_and(goodEnoughPoles,
                                            goodGlobPoles < M))[0]
         removePole = np.where(np.logical_not(goodEnoughPoles))[0]
         if len(removePole) > 0:
             keepCoeff = np.append(keepPole, np.append([N],
                                np.arange(N + 1, len(self.data.HIs[0].coeffs))))
             for hi in self.data.HIs:
                 polyCorrection = np.zeros_like(hi.coeffs[0, :])
                 for j in removePole:
                     if not np.isinf(hi.poles[j]):
                         polyCorrection += hi.coeffs[j, :] / (mu0 - hi.poles[j])
                 if len(hi.coeffs) == N:
                     hi.coeffs = np.vstack((hi.coeffs, polyCorrection))
                 else:
                     hi.coeffs[N, :] += polyCorrection
                 hi.poles = hi.poles[keepPole]
                 hi.coeffs = hi.coeffs[keepCoeff, :]
         for idxR in halfPole:
             pts = np.where(goodLocPoles[:, idxR])[0]
             idxEff = len(self.data.suppEffPts)
             for idEff, prevPts in enumerate(self.data.suppEffPts):
                 if len(prevPts) == len(pts):
                     if np.allclose(prevPts, pts):
                         idxEff = idEff
                         break
             if idxEff == len(self.data.suppEffPts):
                 self.data.suppEffPts += [pts]
             self.data.suppEffIdx[idxR] = idxEff
         self.data.suppEffIdx = self.data.suppEffIdx[keepPole]
         return (" Hard-erased {} pole".format(len(removePole))
               + "s" * (len(removePole) != 1)
               + " and soft-erased {} pole".format(len(halfPole))
               + "s" * (len(halfPole) != 1) + ".")
 
     def _interpolateMarginal(self, mu:paramList, objs:ListAny,
                              mIvals : List[Np2D] = None) -> Np2D:
         res = np.zeros(objs[0].shape + (len(mu),), dtype = objs[0].dtype)
         if mIvals is None: muC = self.centerNormalizeMarginal(mu)
         for suppIdx, suppEffPt in enumerate(self.data.suppEffPts):
             i = np.where(self.data.suppEffIdx == suppIdx)[0]
             if suppIdx == 0:
                 i = np.append(i, np.arange(len(self.data.suppEffIdx),
                                            len(res)))
             if len(i) > 0:
                 if mIvals is None:
                     mIval = self.data.marginalInterp[suppIdx](muC)
                 else:
                     mIval = mIvals[suppIdx]
                 for ij, j in enumerate(suppEffPt):
                     res[i] += np.expand_dims(objs[j][i], - 1) * mIval[ij]
         return res
 
     def interpolateMarginalInterpolator(self, mu : paramList = []):
         """Obtain interpolated approximant interpolator."""
         mu = checkParameterList(mu, self.data.nparMarginal)[0]
         muC = self.centerNormalizeMarginal(mu)
         mIvals = [None] * len(self.data.suppEffPts)
         for suppIdx in range(len(self.data.suppEffPts)):
             mIvals[suppIdx] = self.data.marginalInterp[suppIdx](muC)
         poless = self.interpolateMarginalPoles(mu, mIvals)
         coeffss = self.interpolateMarginalCoeffs(mu, mIvals)
         his = [None] * len(mu)
         for j in range(len(mu)):
             his[j] = HI()
             his[j].poles = poless[..., j]
             his[j].coeffs = coeffss[..., j]
             his[j].npar = 1
             his[j].polybasis = self.data.HIs[0].polybasis
         return his
 
     def interpolateMarginalPoles(self, mu : paramList = [],
                                  mIvals : Np2D = None) -> Np1D:
         """Obtain interpolated approximant poles."""
         mu = checkParameterList(mu, self.data.nparMarginal)[0]
         vbMng(self, "INIT",
               "Interpolating marginal poles at mu = {}.".format(mu), 95)
         intMPoles = self._interpolateMarginal(mu,
                                             [hi.poles for hi in self.data.HIs],
                                             mIvals)
         vbMng(self, "DEL", "Done interpolating marginal poles.", 95)
         return intMPoles
 
     def interpolateMarginalCoeffs(self, mu : paramList = [],
                                   mIvals : Np2D = None) -> Np1D:
         """Obtain interpolated approximant coefficients."""
         mu = checkParameterList(mu, self.data.nparMarginal)[0]
         vbMng(self, "INIT",
               "Interpolating marginal coefficients at mu = {}.".format(mu), 95)
         intMCoeffs = self._interpolateMarginal(mu,
                                            [hi.coeffs for hi in self.data.HIs],
                                            mIvals)
         vbMng(self, "DEL", "Done interpolating marginal coefficients.", 95)
         return intMCoeffs
diff --git a/rrompy/reduction_methods/standard/__init__.py b/rrompy/reduction_methods/standard/__init__.py
index 5155dab..eb87cda 100644
--- a/rrompy/reduction_methods/standard/__init__.py
+++ b/rrompy/reduction_methods/standard/__init__.py
@@ -1,33 +1,31 @@
 # 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 .nearest_neighbor import NearestNeighbor
 from .rational_interpolant import RationalInterpolant
-from .rational_moving_least_squares import RationalMovingLeastSquares
 from .rational_pade import RationalPade
 from .reduced_basis import ReducedBasis
 
 __all__ = [
         'NearestNeighbor',
         'RationalInterpolant',
-        'RationalMovingLeastSquares',
         'RationalPade',
         'ReducedBasis'
           ]
 
 
diff --git a/rrompy/reduction_methods/standard/rational_interpolant.py b/rrompy/reduction_methods/standard/rational_interpolant.py
index 9969d38..c7d73da 100644
--- a/rrompy/reduction_methods/standard/rational_interpolant.py
+++ b/rrompy/reduction_methods/standard/rational_interpolant.py
@@ -1,677 +1,674 @@
 # 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_standard_approximant import GenericStandardApproximant
 from rrompy.utilities.poly_fitting.polynomial import (
                                     polybases as ppb, polyfitname,
                                     polyvander as pvP, polyvanderTotal as pvTP,
                                     polyTimes, polyTimesTable, vanderInvTable,
                                     PolynomialInterpolator as PI)
 from rrompy.utilities.poly_fitting.heaviside import rational2heaviside
 from rrompy.utilities.poly_fitting.radial_basis import (polybases as rbpb,
                                                 RadialBasisInterpolator as RBI)
-from rrompy.utilities.poly_fitting.moving_least_squares import (
-                                        polybases as mlspb,
-                                        MovingLeastSquaresInterpolator as MLSI)
 from rrompy.utilities.base.types import Np1D, Np2D, Tuple, List, sampList
 from rrompy.utilities.base import verbosityManager as vbMng
 from rrompy.utilities.numerical import customPInv, dot, potential
 from rrompy.utilities.numerical.hash_derivative import nextDerivativeIndices
 from rrompy.utilities.numerical.degree import (reduceDegreeN,
                                           degreeTotalToFull, fullDegreeMaxMask,
                                           totalDegreeMaxMask)
 from rrompy.utilities.exception_manager import (RROMPyException, RROMPyAssert,
                                                 RROMPyWarning)
 
 __all__ = ['RationalInterpolant']
 
 class RationalInterpolant(GenericStandardApproximant):
     """
     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;
             - 'scaleFactorDer': scaling factors for derivative computation;
                 defaults to 'AUTO';
             - 'S': total number of samples current approximant relies upon;
             - 'sampler': sample point generator;
             - 'polybasis': type of polynomial basis for interpolation; defaults
                 to 'MONOMIAL';
             - 'M': degree of rational interpolant numerator; defaults to
                 'AUTO', i.e. maximum allowed;
             - 'N': degree of rational interpolant denominator; defaults to
                 'AUTO', i.e. maximum allowed;
             - 'polydegreetype': type of polynomial degree; defaults to 'TOTAL';
             - 'radialDirectionalWeights': radial basis weights for interpolant
                 numerator; defaults to 1;
             - 'interpRcond': tolerance for interpolation; defaults to None;
             - 'robustTol': tolerance for robust rational denominator
                 management; defaults to 0;
             - 'correctorForce': whether corrector should forcefully delete bad
                 poles; defaults to False;
             - 'correctorTol': tolerance for corrector step; defaults to 0.,
                 i.e. no bad poles;
             - 'correctorMaxIter': maximum number of corrector iterations;
                 defaults to 1.
             Defaults to empty dict.
         approx_state(optional): Whether to approximate state. Defaults to
             False.
         verbosity(optional): Verbosity level. Defaults to 10.
             
     Attributes:
         HFEngine: HF problem solver.
         mu0: Default parameter.
         mus: Array of snapshot parameters.
         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;
             - 'scaleFactorDer': scaling factors for derivative computation;
             - 'polybasis': type of polynomial basis for interpolation;
             - 'M': degree of rational interpolant numerator;
             - 'N': degree of rational interpolant denominator;
             - 'polydegreetype': type of polynomial degree;
             - 'radialDirectionalWeights': radial basis weights for interpolant
                 numerator;
             - 'interpRcond': tolerance for interpolation via numpy.polyfit;
             - 'robustTol': tolerance for robust rational denominator
                 management;
             - 'correctorForce': whether corrector should forcefully delete bad
                 poles;
             - 'correctorTol': tolerance for corrector step;
             - 'correctorMaxIter': maximum number of corrector iterations.
         parameterListCritical: Recognized keys of critical approximant
             parameters:
             - 'S': total number of samples current approximant relies upon;
             - 'sampler': sample point generator.
         approx_state: Whether to approximate state.
         verbosity: Verbosity level.
         POD: Whether to compute POD of snapshots.
         scaleFactorDer: Scaling factors for derivative computation.
         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.
         polydegreetype: Type of polynomial degree.
         radialDirectionalWeights: Radial basis weights for interpolant
             numerator.
         interpRcond: Tolerance for interpolation via numpy.polyfit.
         robustTol: Tolerance for robust rational denominator management.
         correctorForce: Whether corrector should forcefully delete bad poles.
         correctorTol: Tolerance for corrector step.
         correctorMaxIter: Maximum number of corrector iterations.
         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, *args, **kwargs):
         self._preInit()
         self._addParametersToList(["polybasis", "M", "N", "polydegreetype",
                                    "radialDirectionalWeights", "interpRcond",
                                    "robustTol", "correctorForce",
                                    "correctorTol", "correctorMaxIter"],
                                   ["MONOMIAL", "AUTO", "AUTO", "TOTAL", [1.],
                                    -1, 0., False, 0., 1])
         super().__init__(*args, **kwargs)
         self.catchInstability = 0
         self._postInit()
 
     @property
     def tModelType(self):
         from .trained_model.trained_model_rational import TrainedModelRational
         return TrainedModelRational
 
     @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 ppb + rbpb + mlspb: 
+            if polybasis not in ppb + rbpb: 
                 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 polybasis0(self):
         if "_" in self.polybasis:
             return self.polybasis.split("_")[0]
         return self.polybasis
 
     @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 radialDirectionalWeights(self):
         """Value of radialDirectionalWeights."""
         return self._radialDirectionalWeights
     @radialDirectionalWeights.setter
     def radialDirectionalWeights(self, radialDirectionalWeights):
         if not hasattr(radialDirectionalWeights, "__len__"):
             radialDirectionalWeights = [radialDirectionalWeights]
         self._radialDirectionalWeights = radialDirectionalWeights
         self._approxParameters["radialDirectionalWeights"] = (
                                                  self.radialDirectionalWeights)
 
     @property
     def M(self):
         """Value of M."""
         return self._M
     @M.setter
     def M(self, M):
         if isinstance(M, str):
             M = M.strip().replace(" ","")
             if "-" not in M: M = M + "-0"
             self._M_isauto, self._M_shift = True, int(M.split("-")[-1])
             M = 0
         if M < 0: raise RROMPyException("M must be non-negative.")
         self._M = M
         self._approxParameters["M"] = self.M
 
     def _setMAuto(self):
         self.M = max(0, reduceDegreeN(self.S, self.S, self.npar,
                                       self.polydegreetype) - self._M_shift)
         vbMng(self, "MAIN", "Automatically setting M to {}.".format(self.M),
               25)
 
     @property
     def N(self):
         """Value of N."""
         return self._N
     @N.setter
     def N(self, N):
         if isinstance(N, str):
             N = N.strip().replace(" ","")
             if "-" not in N: N = N + "-0"
             self._N_isauto, self._N_shift = True, int(N.split("-")[-1])
             N = 0
         if N < 0: raise RROMPyException("N must be non-negative.")
         self._N = N
         self._approxParameters["N"] = self.N
 
     def _setNAuto(self):
         self.N = max(0, reduceDegreeN(self.S, self.S, self.npar,
                                       self.polydegreetype) - self._N_shift)
         vbMng(self, "MAIN", "Automatically setting N to {}.".format(self.N),
               25)
 
     @property
     def polydegreetype(self):
         """Value of polydegreetype."""
         return self._polydegreetype
     @polydegreetype.setter
     def polydegreetype(self, polydegreetype):
         try:
             polydegreetype = polydegreetype.upper().strip().replace(" ","")
             if polydegreetype not in ["TOTAL", "FULL"]: 
                 raise RROMPyException(("Prescribed polydegreetype not "
                                        "recognized."))
             self._polydegreetype = polydegreetype
         except:
             RROMPyWarning(("Prescribed polydegreetype not recognized. "
                            "Overriding to 'TOTAL'."))
             self._polydegreetype = "TOTAL"
         self._approxParameters["polydegreetype"] = self.polydegreetype
 
     @property
     def robustTol(self):
         """Value of tolerance for robust rational 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 correctorForce(self):
         """Value of correctorForce."""
         return self._correctorForce
     @correctorForce.setter
     def correctorForce(self, correctorForce):
         self._correctorForce = correctorForce
         self._approxParameters["correctorForce"] = self.correctorForce
         
     @property
     def correctorTol(self):
         """Value of correctorTol."""
         return self._correctorTol
     @correctorTol.setter
     def correctorTol(self, correctorTol):
         if correctorTol < 0. or (correctorTol > 0. and self.npar > 1):
             RROMPyWarning(("Overriding prescribed corrector tolerance "
                            "to 0."))
             correctorTol = 0.
         self._correctorTol = correctorTol
         self._approxParameters["correctorTol"] = self.correctorTol
 
     @property
     def correctorMaxIter(self):
         """Value of correctorMaxIter."""
         return self._correctorMaxIter
     @correctorMaxIter.setter
     def correctorMaxIter(self, correctorMaxIter):
         if correctorMaxIter < 1 or (correctorMaxIter > 1 and self.npar > 1):
             RROMPyWarning(("Overriding prescribed max number of corrector "
                            "iterations to 1."))
             correctorMaxIter = 1
         self._correctorMaxIter = correctorMaxIter
         self._approxParameters["correctorMaxIter"] = self.correctorMaxIter
         
     def resetSamples(self):
         """Reset samples."""
         super().resetSamples()
         self._musUniqueCN = None
         self._derIdxs = None
         self._reorder = None
 
     def _setupInterpolationIndices(self):
         """Setup parameters for polyvander."""
         if self._musUniqueCN is None or len(self._reorder) != len(self.mus):
             self._musUniqueCN, musIdxsTo, musIdxs, musCount = (
                             self.trainedModel.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 rational denominator."""
         RROMPyAssert(self._mode, message = "Cannot setup denominator.")
         vbMng(self, "INIT", "Starting computation of denominator.", 7)
         if hasattr(self, "_N_isauto"):
             self._setNAuto()
         else:
             N = reduceDegreeN(self.N, self.S, self.npar, self.polydegreetype)
             if N < self.N:
                 RROMPyWarning(("N too large compared to S. Reducing N by "
                                "{}").format(self.N - N))
                 self.N = N
         while self.N > 0:
             invD, fitinv = self._computeInterpolantInverseBlocks()
             idxSamplesEff = list(range(self.S))
             if self.POD:
                 ev, eV = self.findeveVGQR(
                               self.samplingEngine.RPOD[:, idxSamplesEff], invD)
             else:
                 ev, eV = self.findeveVGExplicit(
                               self.samplingEngine.samples(idxSamplesEff), invD)
             nevBad = checkRobustTolerance(ev, self.robustTol)
             if nevBad <= 1: break
             if self.catchInstability > 0:
                 raise RROMPyException(("Instability in denominator "
                                        "computation: eigenproblem is poorly "
                                        "conditioned."),
                                       self.catchInstability == 1)
             vbMng(self, "MAIN", ("Smallest {} eigenvalues below tolerance. "
                                  "Reducing N by 1.").format(nevBad), 10)
             self.N = self.N - 1
         if self.N <= 0:
             self.N = 0
             eV = np.ones((1, 1))
         q = PI()
         q.npar = self.npar
         q.polybasis = self.polybasis0
         if self.polydegreetype == "TOTAL":
             q.coeffs = degreeTotalToFull(tuple([self.N + 1] * self.npar),
                                          self.npar, eV[:, 0])
         else:
             q.coeffs = eV[:, 0].reshape([self.N + 1] * self.npar)
         vbMng(self, "DEL", "Done computing denominator.", 7)
         return q, fitinv
 
     def _setupNumerator(self):
         """Compute rational numerator."""
         RROMPyAssert(self._mode, message = "Cannot setup numerator.")
         vbMng(self, "INIT", "Starting computation of numerator.", 7)
         self._setupInterpolationIndices()
         Qevaldiag = polyTimesTable(self.trainedModel.data.Q, self._musUniqueCN,
                                    self._reorder, self._derIdxs,
                                    self.scaleFactorRel)
         if self.POD:
             Qevaldiag = Qevaldiag.dot(self.samplingEngine.RPOD.T)
         if hasattr(self, "_M_isauto"):
             self._setMAuto()
             M = self.M
         else:
             M = reduceDegreeN(self.M, self.S, self.npar, self.polydegreetype)
             if M < self.M:
                 RROMPyWarning(("M too large compared to S. Reducing M by "
                                "{}").format(self.M - M))
                 self.M = M
         while self.M >= 0:
             pParRest = [self.verbosity >= 5, self.polydegreetype == "TOTAL",
                         {"derIdxs": self._derIdxs, "reorder": self._reorder,
                          "scl": self.scaleFactorRel}]
             if self.polybasis in ppb:
                 p = PI()
             else:
                 pParRest = [self.radialDirectionalWeights] + pParRest
-                p = RBI() if self.polybasis in rbpb else MLSI()
+                p = RBI()
             if self.polybasis in ppb + rbpb:
                 pParRest += [{"rcond": self.interpRcond}]
             wellCond, msg = p.setupByInterpolation(self._musUniqueCN,
                                                    Qevaldiag, self.M, 
                                                    self.polybasis, *pParRest)
             vbMng(self, "MAIN", msg, 5)
             if wellCond: break
             if self.catchInstability > 0:
                 raise RROMPyException(("Instability in numerator computation: "
                                        "polyfit is poorly conditioned."),
                                       self.catchInstability == 1)
             vbMng(self, "MAIN", ("Polyfit is poorly conditioned. Reducing M "
                                  "by 1."), 10)
             self.M = self.M - 1
         if self.M < 0:
             raise RROMPyException(("Instability in computation of numerator. "
                                    "Aborting."))
         self.M = M
         vbMng(self, "DEL", "Done computing numerator.", 7)
         return p
 
     def setupApprox(self) -> int:
         """Compute rational interpolant."""
         if self.checkComputedApprox(): return -1
         RROMPyAssert(self._mode, message = "Cannot setup approximant.")
         vbMng(self, "INIT", "Setting up {}.". format(self.name()), 5)
         self.computeSnapshots()
         pMat = self.samplingEngine.samples.data
         pMatEff = dot(self.HFEngine.C, pMat) if self.approx_state else pMat
         if self.trainedModel is None:
             self.trainedModel = self.tModelType()
             self.trainedModel.verbosity = self.verbosity
             self.trainedModel.timestamp = self.timestamp
             datadict = {"mu0": self.mu0, "projMat": pMatEff,
                         "scaleFactor": self.scaleFactor,
                         "rescalingExp": self.HFEngine.rescalingExp}
             self.trainedModel.data = self.initializeModelData(datadict)[0]
         else:
             self.trainedModel = self.trainedModel
             self.trainedModel.data.projMat = copy(pMatEff)
         self.trainedModel.data.mus = copy(self.mus)
         self._iterCorrector()
         self.trainedModel.data.approxParameters = copy(self.approxParameters)
         vbMng(self, "DEL", "Done setting up approximant.", 5)
         return 0
 
     def _iterCorrector(self):
         if self.correctorTol > 0. and (self.correctorMaxIter > 1
                                     or self.correctorForce):
             vbMng(self, "INIT", "Starting correction iterations.", 5)
             self._Qhat = PI()
             self._Qhat.npar = self.npar
             self._Qhat.polybasis = "MONOMIAL"
             self._Qhat.coeffs = np.ones(1, dtype = np.complex)
             if self.POD:
                 self._RPODOld = copy(self.samplingEngine.RPOD)
             else:
                 self._samplesOld = copy(self.samplingEngine.samples)
         else: vbMng(self, "INIT", "Starting approximant finalization.", 5)
         for j in range(self.correctorMaxIter):
             if self.N > 0 or (hasattr(self, "_N_isauto")
                           and self.S > self.npar):
                 Q = self._setupDenominator()[0]
             else:
                 Q = PI()
                 Q.coeffs = np.ones((1,) * self.npar, dtype = np.complex)
                 Q.npar = self.npar
                 Q.polybasis = self.polybasis
                 self.N = 0
             if j == 0: _N0 = self.N
             self.trainedModel.data.Q = Q
             self.trainedModel.data.P = self._setupNumerator()
             self._applyCorrector(j)
             if self.N <= 0: break
         self.N = _N0
         if self.correctorTol <= 0. or (self.correctorMaxIter <= 1
                                    and not self.correctorForce):
             vbMng(self, "DEL", "Terminated approximant finalization.", 5)
             return
         if self.POD:
             self.samplingEngine.RPOD = self._RPODOld
             del self._RPODOld
         else:
             self.samplingEngine.samples = self._samplesOld
             del self._samplesOld
         if self.correctorForce:
             QOld, QOldBasis = [1.], "MONOMIAL"
         else:
             QOld, QOldBasis = Q.coeffs, self.polybasis
         Q = polyTimes(self._Qhat.coeffs, QOld, Pbasis = self._Qhat.polybasis,
                       Qbasis = QOldBasis, Rbasis = self.polybasis)
         del self._Qhat
         gamma = np.linalg.norm(Q)
         self.trainedModel.data.Q.coeffs = np.pad(Q, (0, self.N - len(Q) + 1),
                                                  "constant") / gamma
         if self.correctorForce:
             self.trainedModel.data.P = self._setupNumerator()
         else:
             self.trainedModel.data.P.coeffs /= gamma
         vbMng(self, "DEL", "Terminated correction iterations.", 5)
 
     def _applyCorrector(self, j:int):
         if self.correctorTol <= 0. or (j >= self.correctorMaxIter - 1
                                    and not self.correctorForce):
             self.N = 0
             return
         cfs, pls, _ = rational2heaviside(self.trainedModel.data.P,
                                          self.trainedModel.data.Q)
         cfs = cfs[: self.N]
         if self.POD:
             resEff = np.linalg.norm(cfs, axis = 1)
         else:    
             resEff = self.HFEngine.norm(self.samplingEngine.samples.dot(cfs.T),
                                         is_state = self.approx_state)
         goodPole = np.logical_not(np.isinf(pls))
         potentialGood = (potential(pls[goodPole], self.sampler.normalFoci())
                        / self.sampler.groundPotential())
         potentialGood[potentialGood < 1.] = 1.
         resEff[goodPole] /= potentialGood
         resEff /= np.max(resEff)
         idxKeep = np.where(resEff >= self.correctorTol)[0]
         vbMng(self, "MAIN",
               ("Correction iteration no. {}: {} out of {} residuals satisfy "
                "tolerance.").format(j + 1, len(idxKeep), self.N), 10)
         self.N -= len(idxKeep)
         if self.N <= 0 and not self.correctorForce: return
         for i in idxKeep:
             self._Qhat.coeffs = polyTimes(self._Qhat.coeffs, [- pls[i], 1.],
                                           Pbasis = self._Qhat.polybasis,
                                           Rbasis = self._Qhat.polybasis)
         self._Qhat.coeffs /= np.linalg.norm(self._Qhat.coeffs)
         if self.N <= 0: return
         vbMng(self, "MAIN",
               ("Removing poles at (normalized positions): "
                "{}.").format(pls[resEff < self.correctorTol]), 65)
         That = polyTimesTable(self._Qhat, self._musUniqueCN,
                               self._reorder, self._derIdxs,
                               self.scaleFactorRel).T
         if self.POD:
             self.samplingEngine.RPOD = self._RPODOld.dot(That)
         else:
             self.samplingEngine.samples = self._samplesOld.dot(That)
 
     def _computeInterpolantInverseBlocks(self) -> Tuple[List[Np2D], Np2D]:
         """
         Compute inverse factors for minimal interpolant target functional.
         """
         RROMPyAssert(self._mode, message = "Cannot solve eigenvalue problem.")
         self._setupInterpolationIndices()
         TEGen = pvTP if self.polydegreetype == "TOTAL" else pvP
         TEGenPar = [self.polybasis0, self._derIdxs, self._reorder,
                     self.scaleFactorRel]
         if hasattr(self, "_M_isauto"): self._setMAuto()
         E = max(self.M, self.N)
         while E >= 0:
             if self.polydegreetype == "TOTAL":
                 Eeff = E
                 idxsB = totalDegreeMaxMask(E, self.npar)
             else: #if self.polydegreetype == "FULL":
                 Eeff = [E] * self.npar
                 idxsB = fullDegreeMaxMask(E, self.npar)
             TE = TEGen(self._musUniqueCN, Eeff, *TEGenPar)
             fitOut = customPInv(TE, rcond = self.interpRcond, full = True)
             vbMng(self, "MAIN",
                   ("Fitting {} samples with degree {} through {}... "
                    "Conditioning of pseudoinverse system: {:.4e}.").format(
                                            TE.shape[0], E,
                                            polyfitname(self.polybasis0),
                                            fitOut[1][1][0] / fitOut[1][1][-1]),
                   5)
             if fitOut[1][0] == TE.shape[1]:
                 fitinv = fitOut[0][idxsB, :]
                 break
             if self.catchInstability > 0:
                 raise RROMPyException(("Instability in denominator "
                                        "computation: polyfit is poorly "
                                        "conditioned."),
                                       self.catchInstability == 1)
             EeqN = E == self.N
             vbMng(self, "MAIN", ("Polyfit is poorly conditioned. Reducing E {}"
                                  "by 1.").format("and N " * EeqN), 10)
             if EeqN: self.N = self.N - 1
             E -= 1
         if self.N < 0:
             raise RROMPyException(("Instability in computation of "
                                    "denominator. Aborting."))
         invD = vanderInvTable(fitinv, idxsB, self._reorder, self._derIdxs)
         if self.N == E:
             TN = TE
         else:
             if self.polydegreetype == "TOTAL":
                 Neff = self.N
                 idxsB = totalDegreeMaxMask(self.N, self.npar)
             else: #if self.polydegreetype == "FULL":
                 Neff = [self.N] * self.npar
                 idxsB = fullDegreeMaxMask(self.N, self.npar)
             TN = TEGen(self._musUniqueCN, Neff, *TEGenPar)
         for k in range(len(invD)): invD[k] = dot(invD[k], TN)
         return invD, fitinv
 
     def findeveVGExplicit(self, sampleE:sampList,
                           invD:List[Np2D]) -> Tuple[Np1D, Np2D]:
         """
         Compute explicitly eigenvalues and eigenvectors of rational 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,
                                              is_state = self.approx_state)
         G = np.zeros((nEnd, nEnd), dtype = np.complex)
         for k in range(eWidth):
              G += dot(dot(gramian, invD[k]).T, invD[k].conj()).T
         vbMng(self, "DEL", "Done building gramian.", 10)
         vbMng(self, "INIT", "Solving eigenvalue problem for gramian matrix.",
               7)
         try:
             ev, eV = np.linalg.eigh(G)
         except np.linalg.LinAlgError as e:
             raise RROMPyException(e)
         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 rational 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, :] = dot(RPODE, invD[k])
         vbMng(self, "DEL", "Done building half-gramian.", 10)
         vbMng(self, "INIT", "Solving svd for square root of gramian matrix.",
               7)
         try:
             _, s, eV = np.linalg.svd(Rstack, full_matrices = False)
         except np.linalg.LinAlgError as e:
             raise RROMPyException(e)
         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 getResidues(self, *args, **kwargs) -> Np1D:
         """
         Obtain approximant residues.
 
         Returns:
             Matrix with residues as columns.
         """
         return self.trainedModel.getResidues(*args, **kwargs)
diff --git a/rrompy/reduction_methods/standard/rational_moving_least_squares.py b/rrompy/reduction_methods/standard/rational_moving_least_squares.py
deleted file mode 100644
index fecdbec..0000000
--- a/rrompy/reduction_methods/standard/rational_moving_least_squares.py
+++ /dev/null
@@ -1,312 +0,0 @@
-# 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 .rational_interpolant import RationalInterpolant
-from rrompy.utilities.poly_fitting.polynomial import (polybases as ppb,
-                                                      polyvander as pvP,
-                                                      polyvanderTotal as pvTP)
-from rrompy.utilities.base.types import Np2D
-from rrompy.utilities.base import verbosityManager as vbMng
-from rrompy.utilities.numerical import dot
-from rrompy.utilities.numerical.degree import (fullDegreeMaxMask,
-                                               totalDegreeMaxMask)
-from rrompy.utilities.exception_manager import (RROMPyException, RROMPyAssert,
-                                                RROMPyWarning)
-
-__all__ = ['RationalMovingLeastSquares']
-
-class RationalMovingLeastSquares(RationalInterpolant):
-    """
-    ROM rational moving LS 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;
-            - 'scaleFactorDer': scaling factors for derivative computation;
-                defaults to 'AUTO';
-            - 'S': total number of samples current approximant relies upon;
-            - 'sampler': sample point generator;
-            - 'polybasis': type of polynomial basis for interpolation; defaults
-                to 'MONOMIAL';
-            - 'M': degree of rational interpolant numerator; defaults to
-                'AUTO', i.e. maximum allowed;
-            - 'N': degree of rational interpolant denominator; defaults to
-                'AUTO', i.e. maximum allowed;
-            - 'polydegreetype': type of polynomial degree; defaults to 'TOTAL';
-            - 'radialBasis': numerator radial basis type; defaults to
-                'GAUSSIAN';
-            - 'radialDirectionalWeights': radial basis weights for interpolant
-                numerator; defaults to 1;
-            - 'radialBasisDen': denominator radial basis type; defaults to
-                'GAUSSIAN';
-            - 'radialDirectionalWeightsDen': radial basis weights for
-                interpolant denominator; defaults to 1;
-            - 'interpRcond': tolerance for interpolation; defaults to None;
-            - 'robustTol': tolerance for robust rational denominator
-                management; defaults to 0.
-            Defaults to empty dict.
-        approx_state(optional): Whether to approximate state. Defaults to
-            False.
-        verbosity(optional): Verbosity level. Defaults to 10.
-            
-    Attributes:
-        HFEngine: HF problem solver.
-        mu0: Default parameter.
-        mus: Array of snapshot parameters.
-        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;
-            - 'scaleFactorDer': scaling factors for derivative computation;
-            - 'polybasis': type of polynomial basis for interpolation;
-            - 'M': degree of rational interpolant numerator;
-            - 'N': degree of rational interpolant denominator;
-            - 'polydegreetype': type of polynomial degree;
-            - 'radialBasis': numerator radial basis type;
-            - 'radialDirectionalWeights': radial basis weights for interpolant
-                numerator;
-            - 'radialBasisDen': denominator radial basis type;
-            - 'radialDirectionalWeightsDen': radial basis weights for
-                interpolant denominator;
-            - 'interpRcond': tolerance for interpolation via numpy.polyfit;
-            - 'robustTol': tolerance for robust rational denominator
-                management.
-        parameterListCritical: Recognized keys of critical approximant
-            parameters:
-            - 'S': total number of samples current approximant relies upon;
-            - 'sampler': sample point generator.
-        approx_state: Whether to approximate state.
-        verbosity: Verbosity level.
-        POD: Whether to compute POD of snapshots.
-        scaleFactorDer: Scaling factors for derivative computation.
-        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.
-        polydegreetype: Type of polynomial degree.
-        radialBasis: Numerator radial basis type.
-        radialDirectionalWeights: Radial basis weights for interpolant
-            numerator.
-        radialBasisDen: Denominator radial basis type.
-        radialDirectionalWeightsDen: Radial basis weights for interpolant
-            denominator.
-        interpRcond: Tolerance for interpolation via numpy.polyfit.
-        robustTol: Tolerance for robust rational denominator management.
-        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, *args, **kwargs):
-        self._preInit()
-        self._addParametersToList(["radialBasis", "radialBasisDen",
-                                   "radialDirectionalWeightsDen"],
-                                  ["GAUSSIAN", "GAUSSIAN", 1.],
-                                  toBeExcluded = ["correctorForce",
-                                                  "correctorTol",
-                                                  "correctorMaxIter"])
-        super().__init__(*args, **kwargs)
-        self.catchInstability = 0
-        self._postInit()
-
-    @property
-    def correctorForce(self):
-        """Value of correctorForce."""
-        return False
-    @correctorForce.setter
-    def correctorForce(self, correctorForce):
-        RROMPyWarning(("correctorForce is used just to simplify inheritance, "
-                       "and its value cannot be changed from False."))
-        
-    @property
-    def correctorTol(self):
-        """Value of correctorTol."""
-        return 0.
-    @correctorTol.setter
-    def correctorTol(self, correctorTol):
-        RROMPyWarning(("correctorTol is used just to simplify inheritance, "
-                       "and its value cannot be changed from 0."))
-        
-    @property
-    def correctorMaxIter(self):
-        """Value of correctorMaxIter."""
-        return 1
-    @correctorMaxIter.setter
-    def correctorMaxIter(self, correctorMaxIter):
-        RROMPyWarning(("correctorMaxIter is used just to simplify "
-                       "inheritance, and its value cannot be changed from 1."))
-        
-    @property
-    def tModelType(self):
-        from .trained_model.trained_model_rational_mls import (
-                                                       TrainedModelRationalMLS)
-        return TrainedModelRationalMLS
-
-    @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 ppb:
-                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):
-        self._radialBasis = radialBasis
-        self._approxParameters["radialBasis"] = self.radialBasis
-
-    @property
-    def radialBasisDen(self):
-        """Value of radialBasisDen."""
-        return self._radialBasisDen
-    @radialBasisDen.setter
-    def radialBasisDen(self, radialBasisDen):
-        self._radialBasisDen = radialBasisDen
-        self._approxParameters["radialBasisDen"] = self.radialBasisDen
-
-    @property
-    def radialDirectionalWeightsDen(self):
-        """Value of radialDirectionalWeightsDen."""
-        return self._radialDirectionalWeightsDen
-    @radialDirectionalWeightsDen.setter
-    def radialDirectionalWeightsDen(self, radialDirectionalWeightsDen):
-        self._radialDirectionalWeightsDen = radialDirectionalWeightsDen
-        self._approxParameters["radialDirectionalWeightsDen"] = (
-                                              self.radialDirectionalWeightsDen)
-
-    def _setupDenominator(self) -> Np2D:
-        """Compute rational denominator."""
-        RROMPyAssert(self._mode, message = "Cannot setup denominator.")
-        vbMng(self, "INIT",
-              "Starting computation of denominator-related blocks.", 7)
-        self._setupInterpolationIndices()
-        pPar = [self._musUniqueCN, self.N, self.polybasis0, self._derIdxs,
-                self._reorder, self.scaleFactorRel]
-        if self.polydegreetype == "TOTAL":
-            TN = pvTP(*pPar)
-        else: #if self.polydegreetype == "FULL":
-            pPar[1] = [pPar[1]] * self.npar
-            TN = pvP(*pPar)
-        TNTen = np.zeros((self.S, self.S, TN.shape[1]), dtype = TN.dtype)
-        TNTen[np.arange(self.S), np.arange(self.S)] = TN
-        if self.POD: TNTen = dot(self.samplingEngine.RPOD, TNTen)
-        vbMng(self, "DEL", "Done computing denominator-related blocks.", 7)
-        return TN, TNTen
-
-    def _setupNumerator(self) -> Np2D:
-        """Compute rational numerator."""
-        RROMPyAssert(self._mode, message = "Cannot setup numerator.")
-        vbMng(self, "INIT",
-              "Starting computation of denominator-related blocks.", 7)
-        self._setupInterpolationIndices()
-        pPar = [self._musUniqueCN, self.M, self.polybasis0, self._derIdxs,
-                self._reorder, self.scaleFactorRel]
-        if self.polydegreetype == "TOTAL":
-            TM = pvTP(*pPar)
-        else: #if self.polydegreetype == "FULL":
-            pPar[1] = [pPar[1]] * self.npar
-            TM = pvP(*pPar)
-        vbMng(self, "DEL", "Done computing denominator-related blocks.", 7)
-        return TM
-
-    def setupApprox(self) -> int:
-        """Compute rational interpolant."""
-        if self.checkComputedApprox(): return -1
-        RROMPyAssert(self._mode, message = "Cannot setup approximant.")
-        vbMng(self, "INIT", "Setting up {}.". format(self.name()), 5)
-        self.computeSnapshots()
-        pMat = self.samplingEngine.samples.data
-        pMatEff = dot(self.HFEngine.C, pMat) if self.approx_state else pMat
-        if self.trainedModel is None:
-            self.trainedModel = self.tModelType()
-            self.trainedModel.verbosity = self.verbosity
-            self.trainedModel.timestamp = self.timestamp
-            datadict = {"mu0": self.mu0, "projMat": pMatEff,
-                        "scaleFactor": self.scaleFactor,
-                        "rescalingExp": self.HFEngine.rescalingExp}
-            data = self.initializeModelData(datadict)[0]
-            data.POD = self.POD
-            data.polybasis = self.polybasis
-            data.polydegreetype = self.polydegreetype
-            data.radialBasis = self.radialBasis
-            data.radialWeights = self.radialDirectionalWeights
-            data.radialBasisDen = self.radialBasisDen
-            data.radialWeightsDen = self.radialDirectionalWeightsDen
-            data.interpRcond = self.interpRcond
-            self.trainedModel.data = data
-        else:
-            self.trainedModel = self.trainedModel
-            self.trainedModel.data.projMat = copy(pMatEff)
-        if not self.POD:
-            self.trainedModel.data.gramian = self.HFEngine.innerProduct(
-                                                  self.samplingEngine.samples,
-                                                  self.samplingEngine.samples,
-                                                  is_state = self.approx_state)
-        self.trainedModel.data.mus = copy(self.mus)
-        self.trainedModel.data.M = self.M
-        self.trainedModel.data.N = self.N
-        QVan, self.trainedModel.data.QBlocks = self._setupDenominator()
-        self.trainedModel.data.PVan = self._setupNumerator()
-        if self.polydegreetype == "TOTAL":
-            degreeMaxMask = totalDegreeMaxMask
-        else: #if self.polydegreetype == "FULL":
-            degreeMaxMask = fullDegreeMaxMask
-        if self.N > self.M:
-            self.trainedModel.data.QVan = QVan
-            self.trainedModel.data.domQIdxs = degreeMaxMask(self.N, self.npar)
-        else:
-            self.trainedModel.data.domQIdxs = degreeMaxMask(self.M, self.npar)
-        self.trainedModel.data.approxParameters = copy(self.approxParameters)
-        vbMng(self, "DEL", "Done setting up approximant.", 5)
-        return 0
-
diff --git a/rrompy/reduction_methods/standard/rational_pade.py b/rrompy/reduction_methods/standard/rational_pade.py
index d75461a..6b48d90 100644
--- a/rrompy/reduction_methods/standard/rational_pade.py
+++ b/rrompy/reduction_methods/standard/rational_pade.py
@@ -1,309 +1,307 @@
 # 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.reduction_methods.base import checkRobustTolerance
 from .rational_interpolant import RationalInterpolant
 from rrompy.utilities.poly_fitting.polynomial import (
                                     polybases as ppb, polyfitname,
                                     polyvander as pvP, polyvanderTotal as pvTP,
                                     polyTimesTable, vanderInvTable,
                                     PolynomialInterpolator as PI)
 from rrompy.utilities.poly_fitting.radial_basis import (polybases as rbpb,
                                                 RadialBasisInterpolator as RBI)
-from rrompy.utilities.poly_fitting.moving_least_squares import (
-                                        MovingLeastSquaresInterpolator as MLSI)
 from rrompy.utilities.base.types import Np2D, Tuple, List
 from rrompy.utilities.base import verbosityManager as vbMng
 from rrompy.utilities.numerical import customPInv, dot
 from rrompy.utilities.numerical.degree import (fullDegreeN, totalDegreeN,
                                          reduceDegreeN, degreeTotalToFull,
                                          fullDegreeMaxMask, totalDegreeMaxMask)
 from rrompy.utilities.exception_manager import (RROMPyException, RROMPyAssert,
                                                 RROMPyWarning)
 
 __all__ = ['RationalPade']
 
 class RationalPade(RationalInterpolant):
     """
     ROM rational Pade' 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;
             - 'scaleFactorDer': scaling factors for derivative computation;
                 defaults to 'AUTO';
             - 'S': total number of samples current approximant relies upon;
             - 'sampler': sample point generator;
             - 'polybasis': type of polynomial basis for interpolation; defaults
                 to 'MONOMIAL';
             - 'M': degree of rational interpolant numerator; defaults to
                 'AUTO', i.e. maximum allowed;
             - 'N': degree of rational interpolant denominator; defaults to
                 'AUTO', i.e. maximum allowed;
             - 'polydegreetype': type of polynomial degree; defaults to 'TOTAL';
             - 'radialDirectionalWeights': radial basis weights for interpolant
                 numerator; defaults to 1;
             - 'interpRcond': tolerance for interpolation; defaults to None;
             - 'robustTol': tolerance for robust rational denominator
                 management; defaults to 0;
             - 'correctorForce': whether corrector should forcefully delete bad
                 poles; defaults to False;
             - 'correctorTol': tolerance for corrector step; defaults to 0.,
                 i.e. no bad poles;
             - 'correctorMaxIter': maximum number of corrector iterations;
                 defaults to 1.
             Defaults to empty dict.
         approx_state(optional): Whether to approximate state. Defaults to
             False.
         verbosity(optional): Verbosity level. Defaults to 10.
 
     Attributes:
         HFEngine: HF problem solver.
         mu0: Default parameter.
         mus: Array of snapshot parameters.
         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;
             - 'scaleFactorDer': scaling factors for derivative computation;
             - 'polybasis': type of polynomial basis for interpolation;
             - 'M': degree of rational interpolant numerator;
             - 'N': degree of rational interpolant denominator;
             - 'polydegreetype': type of polynomial degree;
             - 'radialDirectionalWeights': radial basis weights for interpolant
                 numerator;
             - 'interpRcond': tolerance for interpolation via numpy.polyfit;
             - 'robustTol': tolerance for robust rational denominator
                 management;
             - 'correctorForce': whether corrector should forcefully delete bad
                 poles;
             - 'correctorTol': tolerance for corrector step;
             - 'correctorMaxIter': maximum number of corrector iterations.
         parameterListCritical: Recognized keys of critical approximant
             parameters:
             - 'S': total number of samples current approximant relies upon;
             - 'sampler': sample point generator.
         approx_state: Whether to approximate state.
         verbosity: Verbosity level.
         POD: Whether to compute POD of snapshots.
         scaleFactorDer: Scaling factors for derivative computation.
         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.
         polydegreetype: Type of polynomial degree.
         radialDirectionalWeights: Radial basis weights for interpolant
             numerator.
         interpRcond: Tolerance for interpolation via numpy.polyfit.
         robustTol: Tolerance for robust rational denominator management.
         correctorForce: Whether corrector should forcefully delete bad poles.
         correctorTol: Tolerance for corrector step.
         correctorMaxIter: Maximum number of corrector iterations.
         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 _setupInterpolationIndices(self):
         """Setup parameters for polyvander."""
         super()._setupInterpolationIndices()
         if len(self._musUniqueCN) > 1:
             raise RROMPyException(("Cannot apply centered-like method with "
                                    "more than one distinct sample point."))
 
     def _setupDenominator(self):
         """Compute rational denominator."""
         RROMPyAssert(self._mode, message = "Cannot setup denominator.")
         vbMng(self, "INIT", "Starting computation of denominator.", 7)
         cfun = totalDegreeN if self.polydegreetype == "TOTAL" else fullDegreeN
         if hasattr(self, "_N_isauto"):
             self._setNAuto()
         else:
             N = reduceDegreeN(self.N, self.S, self.npar, self.polydegreetype)
             if N < self.N:
                 RROMPyWarning(("N too large compared to S. Reducing N by "
                                "{}").format(self.N - N))
                 self.N = N
         while self.N > 0:
             invD, fitinv = self._computeInterpolantInverseBlocks()
             Seff = cfun(self.N, self.npar)
             idxSamplesEff = list(range(self.S - Seff, self.S))
             if self.POD:
                 ev, eV = self.findeveVGQR(
                               self.samplingEngine.RPOD[:, idxSamplesEff], invD)
             else:
                 ev, eV = self.findeveVGExplicit(
                               self.samplingEngine.samples(idxSamplesEff), invD)
             nevBad = checkRobustTolerance(ev, self.robustTol)
             if nevBad <= 1: break
             if self.catchInstability > 0:
                 raise RROMPyException(("Instability in denominator "
                                        "computation: eigenproblem is poorly "
                                        "conditioned."),
                                       self.catchInstability == 1)
             RROMPyWarning(("Smallest {} eigenvalues below tolerance. Reducing "
                            "N by 1.").format(nevBad))
             self.N = self.N - 1
         if self.N <= 0:
             self.N = 0
             eV = np.ones((1, 1))
         q = PI()
         q.npar = self.npar
         q.polybasis = self.polybasis0
         if self.polydegreetype == "TOTAL":
             q.coeffs = degreeTotalToFull(tuple([self.N + 1] * self.npar),
                                          self.npar, eV[:, 0])
         else:
             q.coeffs = eV[:, 0].reshape([self.N + 1] * self.npar)
         vbMng(self, "DEL", "Done computing denominator.", 7)
         return q, fitinv
 
     def _setupNumerator(self):
         """Compute rational numerator."""
         RROMPyAssert(self._mode, message = "Cannot setup numerator.")
         vbMng(self, "INIT", "Starting computation of numerator.", 7)
         self._setupInterpolationIndices()
         Qevaldiag = polyTimesTable(self.trainedModel.data.Q, self._musUniqueCN,
                                    self._reorder, self._derIdxs,
                                    self.scaleFactorRel)
         if self.POD:
             Qevaldiag = Qevaldiag.dot(self.samplingEngine.RPOD.T)
         cfun = totalDegreeN if self.polydegreetype == "TOTAL" else fullDegreeN
         if hasattr(self, "_M_isauto"):
             self._setMAuto()
             M = self.M
         else:
             M = reduceDegreeN(self.M, self.S, self.npar, self.polydegreetype)
             if M < self.M:
                 RROMPyWarning(("M too large compared to S. Reducing M by "
                                "{}").format(self.M - M))
                 self.M = M
         while self.M >= 0:
             Seff = cfun(self.M, self.npar)
             pParRest = [self.verbosity >= 5, self.polydegreetype == "TOTAL",
                         {"derIdxs": [self._derIdxs[0][: Seff]],
                          "reorder": self._reorder[: Seff],
                          "scl": self.scaleFactorRel}]
             if self.polybasis in ppb:
                 p = PI()
             else:
                 pParRest = [self.radialDirectionalWeights] + pParRest
-                p = RBI() if self.polybasis in rbpb else MLSI()
+                p = RBI()
             if self.polybasis in ppb + rbpb:
                 pParRest += [{"rcond": self.interpRcond}]
             wellCond, msg = p.setupByInterpolation(self._musUniqueCN,
                                                    Qevaldiag[: Seff, : Seff],
                                                    self.M, self.polybasis,
                                                    *pParRest)
             vbMng(self, "MAIN", msg, 5)
             if wellCond: break
             if self.catchInstability > 0:
                 raise RROMPyException(("Instability in numerator computation: "
                                        "polyfit is poorly conditioned."),
                                       self.catchInstability == 1)
             vbMng(self, "MAIN", ("Polyfit is poorly conditioned. Reducing M "
                                  "by 1."), 10)
             self.M = self.M - 1
         if self.M < 0:
             raise RROMPyException(("Instability in computation of numerator. "
                                    "Aborting."))
         self.M = M
         vbMng(self, "DEL", "Done computing numerator.", 7)
         return p
 
     def _computeInterpolantInverseBlocks(self) -> Tuple[List[Np2D], Np2D]:
         """
         Compute inverse factors for minimal interpolant target functional.
         """
         RROMPyAssert(self._mode, message = "Cannot solve eigenvalue problem.")
         self._setupInterpolationIndices()
         if self.polydegreetype == "TOTAL":
             cfun, TEGen = totalDegreeN, pvTP
         else:
             cfun, TEGen = fullDegreeN, pvP
         E = max(self.M, self.N)
         while E >= 0:
             Seff = cfun(E, self.npar)
             TEGenPar = [self.polybasis0, [self._derIdxs[0][: Seff]],
                         self._reorder[: Seff], self.scaleFactorRel]
             if self.polydegreetype == "TOTAL":
                 Eeff = E
                 idxsB = totalDegreeMaxMask(E, self.npar)
             else: #if self.polydegreetype == "FULL":
                 Eeff = [E] * self.npar
                 idxsB = fullDegreeMaxMask(E, self.npar)
             TE = TEGen(self._musUniqueCN, Eeff, *TEGenPar)
             fitOut = customPInv(TE, rcond = self.interpRcond, full = True)
             vbMng(self, "MAIN",
                   ("Fitting {} samples with degree {} through {}... "
                    "Conditioning of pseudoinverse system: {:.4e}.").format(
                                            TE.shape[0], E,
                                            polyfitname(self.polybasis0),
                                            fitOut[1][1][0] / fitOut[1][1][-1]),
                   5)
             if fitOut[1][0] == TE.shape[1]:
                 fitinv = fitOut[0][idxsB, :]
                 break
             if self.catchInstability > 0:
                 raise RROMPyException(("Instability in denominator "
                                        "computation: polyfit is poorly "
                                        "conditioned."),
                                       self.catchInstability == 1)
             EeqN = E == self.N
             vbMng(self, "MAIN", ("Polyfit is poorly conditioned. Reducing E {}"
                                  "by 1.").format("and N " * EeqN), 10)
             if EeqN: self.N = self.N - 1
             E -= 1
         if self.N < 0:
             raise RROMPyException(("Instability in computation of "
                                    "denominator. Aborting."))
         invD = vanderInvTable(fitinv, idxsB, self._reorder[: Seff],
                               [self._derIdxs[0][: Seff]])
         if self.N == E:
             TN = TE
         else:
             if self.polydegreetype == "TOTAL":
                 Neff = self.N
                 idxsB = totalDegreeMaxMask(self.N, self.npar)
             else: #if self.polydegreetype == "FULL":
                 Neff = [self.N] * self.npar
                 idxsB = fullDegreeMaxMask(self.N, self.npar)
             TN = TEGen(self._musUniqueCN, Neff, *TEGenPar)
         for k in range(len(invD)): invD[k] = dot(invD[k], TN)
         return invD, fitinv
diff --git a/rrompy/reduction_methods/standard/trained_model/trained_model_rational_mls.py b/rrompy/reduction_methods/standard/trained_model/trained_model_rational_mls.py
deleted file mode 100644
index d2fcaa3..0000000
--- a/rrompy/reduction_methods/standard/trained_model/trained_model_rational_mls.py
+++ /dev/null
@@ -1,189 +0,0 @@
-# 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 .trained_model_rational import TrainedModelRational
-from rrompy.utilities.base.types import Np1D, paramVal, paramList, sampList
-from rrompy.utilities.base import verbosityManager as vbMng
-from rrompy.utilities.poly_fitting.moving_least_squares import mlsweights
-from rrompy.utilities.poly_fitting.polynomial import (
-                                                  PolynomialInterpolator as PI)
-from rrompy.utilities.numerical import customPInv
-from rrompy.utilities.numerical.degree import degreeTotalToFull
-from rrompy.parameter import checkParameterList
-from rrompy.sampling import emptySampleList
-from rrompy.utilities.exception_manager import RROMPyException, RROMPyWarning
-
-__all__ = ['TrainedModelRationalMLS']
-
-class TrainedModelRationalMLS(TrainedModelRational):
-    """
-    ROM approximant evaluation for rational moving least squares approximant.
-    
-    Attributes:
-        Data: dictionary with all that can be pickled.
-    """
-
-    def reset(self):
-        super().reset()
-        self.lastSetupMu = None
-
-    def assembleReducedModel(self, mu:paramVal):
-        if not (hasattr(self.data, "lastSetupMu")
-            and self.data.lastSetupMu == mu):
-            vbMng(self, "INIT", "Assembling reduced model for mu = {}."\
-                                                               .format(mu), 17)
-            vbMng(self, "INIT", "Starting computation of denominator.", 35)
-            muC = self.centerNormalize(mu)
-            muSC = self.centerNormalize(self.data.mus)
-            wQ = mlsweights(muC, muSC, self.data.radialBasisDen,
-                            directionalWeights = self.data.radialWeightsDen)
-            if self.data.N > self.data.M:
-                PQVan = self.data.QVan
-            else:
-                PQVan = self.data.PVan
-            VQAdjW = PQVan.conj().T * wQ
-            VQAdjWVQ = VQAdjW.dot(PQVan)
-            interpPseudoInverse, info = customPInv(VQAdjWVQ, full = True,
-                                                 rcond = self.data.interpRcond)
-            interpPseudoInverse = interpPseudoInverse.dot(VQAdjW).dot(
-                                                             self.data.QBlocks)
-            if info[0] < interpPseudoInverse.shape[-1]:
-                q = np.zeros(interpPseudoInverse.shape[-1], dtype = np.complex)
-                q[0] = 1.
-            else:
-                halfGram = interpPseudoInverse[self.data.domQIdxs]
-                if self.data.POD:
-                    Rstack = halfGram.reshape(-1, halfGram.shape[-1])
-                    vbMng(self, "INIT",
-                          "Solving svd for square root of gramian matrix.", 67)
-                    try:
-                        _, s, eV = np.linalg.svd(Rstack, full_matrices = False)
-                    except np.linalg.LinAlgError as e:
-                        raise RROMPyException(e)
-                    condN = s[0] / s[-1]
-                    q = eV[-1, :].T.conj()
-                    vbMng(self, "MAIN", 
-                          ("Solved svd problem of size {} x {} with condition "
-                           "number {:.4e}.").format(*Rstack.shape, condN), 55)
-                    vbMng(self, "DEL", "Done solving svd.", 67)
-                else:
-                    RRstack = np.tensordot(self.trainedModel.gramian, halfGram,
-                                           1).reshape(-1, halfGram.shape[-1])
-                    RLstack = halfGram.reshape(-1, halfGram.shape[-1])
-                    gram = RLstack.T.conj().dot(RRstack)
-                    vbMng(self, "INIT",
-                          "Solving eigenvalue problem for gramian matrix.", 67)
-                    try:
-                        ev, eV = np.linalg.eigh(gram)
-                    except np.linalg.LinAlgError as e:
-                        raise RROMPyException(e)
-                    condN = ev[-1] / ev[0]
-                    q = eV[:, 0]
-                    vbMng(self, "MAIN",
-                          ("Solved eigenvalue problem of size {} with "
-                           "condition number {:.4e}.").format(gram.shape[0],
-                                                              condN), 55)
-                    vbMng(self, "DEL", "Done solving eigenvalue problem.", 67)
-            self.data.Q = PI()
-            self.data.Q.npar = self.npar
-            self.data.Q.polybasis = self.data.polybasis
-            if self.data.polydegreetype == "TOTAL":
-                self.data.Q.coeffs = degreeTotalToFull(
-                                                (self.data.N + 1,) * self.npar,
-                                                self.npar, q)
-            else:
-                self.data.Q.coeffs = q.reshape((self.data.N + 1,) * self.npar)
-            vbMng(self, "DEL", "Done computing denominator.", 35)
-            vbMng(self, "INIT", "Starting computation of numerator.", 35)
-            self.data.P = PI()
-            self.data.P.npar = self.npar
-            self.data.P.polybasis = self.data.polybasis
-            wP = mlsweights(muC, muSC, self.data.radialBasis,
-                            directionalWeights = self.data.radialWeights)
-            VAdjW = self.data.PVan.conj().T * wP
-            VAdjWV = VAdjW.dot(self.data.PVan)
-            interpPPseudoInverse = customPInv(VAdjWV, self.data.interpRcond)
-            Pcoeffs = np.tensordot(interpPPseudoInverse.dot(VAdjW),
-                                   self.data.QBlocks.dot(q), ([1], [1]))
-            if self.data.polydegreetype == "TOTAL":
-                self.data.P.coeffs = degreeTotalToFull(
-                                                 (self.data.M + 1,) * self.npar
-                                               + (self.data.QBlocks.shape[0],),
-                                                 self.npar, Pcoeffs)
-            else:
-                self.data.P.coeffs = Pcoeffs.reshape(
-                                                 (self.data.M + 1,) * self.npar
-                                               + (self.data.QBlocks.shape[0],))
-            vbMng(self, "DEL", "Done computing numerator.", 35)
-            vbMng(self, "DEL", "Done assembling reduced model.", 17)
-            self.data.lastSetupMu = mu
-
-    def getApproxReduced(self, mu : paramList = []) -> sampList:
-        """
-        Evaluate reduced representation of approximant at arbitrary parameter.
-
-        Args:
-            mu: Target parameter.
-        """
-        mu = checkParameterList(mu, self.data.npar)[0]
-        if (not hasattr(self, "lastSolvedApproxReduced")
-         or self.lastSolvedApproxReduced != mu):
-            vbMng(self, "INIT",
-                  "Evaluating approximant at mu = {}.".format(mu), 12)
-            self.uApproxReduced = emptySampleList()
-            for i in range(len(mu)):
-                self.assembleReducedModel(mu[i])
-                vbMng(self, "INIT",
-                      "Solving reduced model for mu = {}.".format(mu[i]), 15)
-                Qv = self.getQVal(mu[i])
-                if Qv == 0.:
-                    RROMPyWarning(("Adjusting approximation to avoid division "
-                                   "by numerically zero denominator."))
-                    Qv = np.finfo(np.complex).eps / (1. + self.data.Q.deg[0])
-                uAppR = self.getPVal(mu[i]) / Qv
-                if i == 0:
-                    self.uApproxReduced.reset((uAppR.shape[0], len(mu)),
-                                              dtype = uAppR.dtype)
-                self.uApproxReduced[i] = uAppR
-                vbMng(self, "DEL", "Done solving reduced model.", 15)
-            vbMng(self, "DEL", "Done evaluating approximant.", 12)
-            self.lastSolvedApproxReduced = mu
-        return self.uApproxReduced
-    
-    def getPoles(self, *args, mu : paramVal = None, **kwargs) -> Np1D:
-        """
-        Obtain approximant poles.
-
-        Returns:
-            Numpy complex vector of poles.
-        """
-        if mu is None: mu = self.data.mu0
-        self.assembleReducedModel(mu)
-        return super().getPoles(*args, **kwargs)
-
-    def getResidues(self, *args, mu : paramVal = None, **kwargs) -> Np1D:
-        """
-        Obtain approximant residues.
-
-        Returns:
-            Numpy matrix with residues as columns.
-        """
-        if mu is None: mu = self.data.mu0
-        self.assembleReducedModel(mu)
-        return super().getResidues(*args, **kwargs)
diff --git a/rrompy/utilities/poly_fitting/moving_least_squares/__init__.py b/rrompy/utilities/poly_fitting/moving_least_squares/__init__.py
deleted file mode 100644
index c126e42..0000000
--- a/rrompy/utilities/poly_fitting/moving_least_squares/__init__.py
+++ /dev/null
@@ -1,40 +0,0 @@
-# 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,
-                     localWendland)
-from .base import mlsbases, polybases, polyfitname, polydomcoeff
-from .vander import mlsweights, polyvander, polyvanderTotal
-from .moving_least_squares_interpolator import MovingLeastSquaresInterpolator
-
-__all__ = [
-        'radialGaussian',
-        'thinPlateSpline',
-        'multiQuadric',
-        'localWendland',
-        'mlsbases',
-        'polybases',
-        'polyfitname',
-        'polydomcoeff',
-        'mlsweights',
-        'polyvander',
-        'polyvanderTotal',
-        'MovingLeastSquaresInterpolator'
-           ]
-
-
diff --git a/rrompy/utilities/poly_fitting/moving_least_squares/base.py b/rrompy/utilities/poly_fitting/moving_least_squares/base.py
deleted file mode 100644
index 155b780..0000000
--- a/rrompy/utilities/poly_fitting/moving_least_squares/base.py
+++ /dev/null
@@ -1,27 +0,0 @@
-# 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 rrompy.utilities.poly_fitting.radial_basis.base import (
-                         rbbases, polybases as rbpb, polyfitname, polydomcoeff)
-
-__all__ = ['mlsbases', 'polybases', 'polyfitname', 'polydomcoeff']
-
-mlsbases = [rbb + "_MLS" for rbb in rbbases]
-
-polybases = [rbp + "_MLS" for rbp in rbpb]
-
diff --git a/rrompy/utilities/poly_fitting/moving_least_squares/kernel.py b/rrompy/utilities/poly_fitting/moving_least_squares/kernel.py
deleted file mode 100644
index 207a8f1..0000000
--- a/rrompy/utilities/poly_fitting/moving_least_squares/kernel.py
+++ /dev/null
@@ -1,24 +0,0 @@
-# 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 rrompy.utilities.poly_fitting.radial_basis.kernel import (radialGaussian,
-                                  thinPlateSpline, multiQuadric, localWendland)
-
-__all__ = ['radialGaussian', 'thinPlateSpline', 'multiQuadric',
-           'localWendland']
-
diff --git a/rrompy/utilities/poly_fitting/moving_least_squares/moving_least_squares_interpolator.py b/rrompy/utilities/poly_fitting/moving_least_squares/moving_least_squares_interpolator.py
deleted file mode 100644
index 37440d2..0000000
--- a/rrompy/utilities/poly_fitting/moving_least_squares/moving_least_squares_interpolator.py
+++ /dev/null
@@ -1,139 +0,0 @@
-# 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, Np2D,
-                                         paramList)
-from rrompy.utilities.numerical import customPInv, dot
-from .vander import mlsweights
-from rrompy.utilities.poly_fitting.interpolator import GenericInterpolator
-from rrompy.utilities.poly_fitting.polynomial.vander import (polyvander as pv,
-                                                        polyvanderTotal as pvT)
-from rrompy.utilities.exception_manager import RROMPyException, RROMPyAssert
-from rrompy.parameter import checkParameterList
-
-__all__ = ['MovingLeastSquaresInterpolator']
-
-class MovingLeastSquaresInterpolator(GenericInterpolator):
-    def __init__(self, other = None):
-        if other is None: return
-        self.support = other.support
-        self.localProjector = other.localProjector
-        self.localVanders = other.localVanders
-        self.suppValues = other.suppValues
-        self.directionalWeights = other.directionalWeights
-        self.degree = other.degree
-        self.npar = other.npar
-        self.radialbasis = other.radialbasis
-        self.polybasis = other.polybasis
-        self.evalParams = other.evalParams
-        self.totalDegree = other.totalDegree
-
-    @property
-    def shape(self):
-        sh = self.suppValues.shape[1 :] if self.suppValues.ndim > 1 else 1
-        return sh
-
-    @property
-    def deg(self):
-        return self.degree
-
-    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 moving least "
-                                   "squares function."))
-        mu = checkParameterList(mu, self.npar)[0]
-        sh = self.shape
-        if sh == 1: sh = tuple([])
-        values = np.empty((len(mu),) + sh, dtype = np.complex)
-        for i, m in enumerate(mu):
-            weights = mlsweights(m, self.support, self.radialbasis,
-                                 directionalWeights = self.directionalWeights)
-            weights /= np.linalg.norm(weights)
-            vanderLS = np.sum(self.localVanders * weights, axis = 2)
-            RHSLS = dot(self.localProjector * weights, self.suppValues)
-            if self.totalDegree:
-                vanderEval = pvT(m, self.deg[0], self.polybasis,
-                                 **self.evalParams)
-            else:
-                vanderEval = pv(m, self.deg, self.polybasis, **self.evalParams)
-            vanderEval = vanderEval.flatten()
-            values[i] = dot(vanderEval, dot(customPInv(vanderLS), RHSLS))
-        return values
-        
-    def __copy__(self):
-        return MovingLeastSquaresInterpolator(self)
-
-    def __deepcopy__(self, memo):
-        other = MovingLeastSquaresInterpolator()
-        (other.support, other.localProjector, other.localVanders,
-         other.suppValues, other.directionalWeights, other.degree, other.npar,
-         other.radialbasis, other.polybasis, other.evalParams,
-         other.totalDegree) = copy(
-                        (self.support, self.localProjector, self.localVanders,
-                         self.suppValues, self.directionalWeights, self.degree,
-                         self.npar, self.radialbasis, self.polybasis,
-                         self.evalParams, self.totalDegree), memo)
-        return other
-
-    def postmultiplyTensorize(self, A:Np2D):
-        RROMPyAssert(A.shape[0], self.shape[-1], "Shape of output")
-        self.suppValues = np.tensordot(self.suppValues, A, 1)
-
-    def pad(self, nleft : List[int] = None, nright : List[int] = None):
-        if nleft is None: nleft = [0] * len(self.shape)
-        if nright is None: nright = [0] * len(self.shape)
-        if not hasattr(nleft, "__len__"): nleft = [nleft]
-        if not hasattr(nright, "__len__"): nright = [nright]
-        RROMPyAssert(len(self.shape), len(nleft), "Shape of output")
-        RROMPyAssert(len(self.shape), len(nright), "Shape of output")
-        padwidth = [(0, 0)] + [(l, r) for l, r in zip(nleft, nright)]
-        self.suppValues = np.pad(self.suppValues, padwidth, "constant",
-                                 constant_values = (0., 0.))
-
-    def setupByInterpolation(self, support:paramList, values:ListAny,
-                             deg:int, polybasis : str = "MONOMIAL_GAUSSIAN",
-                             directionalWeights : Np1D = None,
-                             totalDegree : bool = True,
-                             vanderCoeffs : DictAny = {}):
-        support = checkParameterList(support)[0]
-        self.support = copy(support)
-        if "reorder" in vanderCoeffs.keys():
-            self.support = self.support[vanderCoeffs["reorder"]]
-        self.npar = support.shape[1]
-        if directionalWeights is None:
-            directionalWeights = np.ones(self.npar)
-        self.directionalWeights = directionalWeights
-        self.polybasis, self.radialbasis, _ = polybasis.split("_")
-        self.totalDegree = totalDegree
-        self.evalParams = vanderCoeffs
-        if totalDegree:
-            vander = pvT(support, deg, self.polybasis, **vanderCoeffs)
-            if not hasattr(deg, "__len__"): deg = [deg] * self.npar
-        else:
-            if not hasattr(deg, "__len__"): deg = [deg] * self.npar
-            vander = pv(support, deg, self.polybasis, **vanderCoeffs)
-        self.degree = deg
-        self.localProjector = vander.T.conj()
-        self.localVanders = np.array([np.outer(van, van.conj()) \
-                                                            for van in vander])
-        self.localVanders = np.swapaxes(self.localVanders, 0, 2)
-        self.suppValues = np.array(values)
-
diff --git a/rrompy/utilities/poly_fitting/moving_least_squares/vander.py b/rrompy/utilities/poly_fitting/moving_least_squares/vander.py
deleted file mode 100644
index f378f57..0000000
--- a/rrompy/utilities/poly_fitting/moving_least_squares/vander.py
+++ /dev/null
@@ -1,86 +0,0 @@
-# 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 .kernel import (radialGaussian, thinPlateSpline, multiQuadric,
-                     localWendland)
-from rrompy.utilities.poly_fitting.polynomial.vander import (polyvander as pvP,
-                                                       polyvanderTotal as pvTP)
-from rrompy.utilities.base.types import (Np1D, Np2D, Tuple, List, paramVal,
-                                         paramList)
-from rrompy.parameter import checkParameter, checkParameterList
-from rrompy.utilities.exception_manager import RROMPyException, RROMPyAssert
-
-__all__ = ['mlsweights', 'polyvander', 'polyvanderTotal']
-
-def mlsweights(x:paramVal, xSupp:paramList, basis:str,
-               reorder : List[int] = None,
-               directionalWeights : Np1D = None) -> Np2D:
-    """Compute moving least squares weight vector."""
-    x = checkParameter(x)
-    xSupp = checkParameterList(xSupp)[0]
-    x = x.data
-    xSupp = xSupp.data
-    if directionalWeights is None:
-        directionalWeights = np.ones(x.shape[1])
-    elif not hasattr(directionalWeights, "__len__"):
-        directionalWeights = directionalWeights * np.ones(x.shape[1])
-    RROMPyAssert(len(directionalWeights), x.shape[1],
-                 "Number of directional weights")
-    try:
-        radialkernel = {"GAUSSIAN" : radialGaussian,
-                        "THINPLATE" : thinPlateSpline,
-                        "MULTIQUADRIC" : multiQuadric,
-                        "WENDLAND" : localWendland}[basis.upper()]
-    except:
-        raise RROMPyException("Radial basis not recognized.")
-    if reorder is not None: xSupp = xSupp[reorder]
-    muDiff = (xSupp.data - x) * directionalWeights
-    r2 = np.sum(np.abs(muDiff) ** 2., axis = 1).reshape(1, -1)
-    return radialkernel(r2)[0]
-    
-def polyvander(x:paramVal, xSupp:paramList, degs:List[int], basis:str,
-               derIdxs : List[List[List[int]]] = None,
-               reorder : List[int] = None, directionalWeights : Np1D = None,
-               scl : Np1D = None) -> Tuple[Np2D, Np2D]:
-    """
-    Compute full Hermite-Vandermonde matrix with specified derivative
-        directions.
-    """
-    basisp, basisr, _ = basis.split("_")
-    Weights = mlsweights(x, xSupp, basisr, reorder, directionalWeights)
-    VanP = pvP(xSupp, degs, basisp, derIdxs = derIdxs, reorder = reorder,
-               scl = scl)
-    RHP = VanP.T.conj() * Weights
-    return RHP.dot(VanP), RHP
-
-def polyvanderTotal(x:paramList, xSupp:paramList, deg:int, basis:str,
-                    derIdxs : List[List[List[int]]] = None,
-                    reorder : List[int] = None,
-                    directionalWeights : Np1D = None,
-                    scl : Np1D = None) -> Tuple[Np2D, Np2D]:
-    """
-    Compute full total degree Hermite-Vandermonde matrix with specified
-        derivative directions.
-    """
-    basisp, basisr, _ = basis.split("_")
-    Weights = mlsweights(x, xSupp, basisr, reorder, directionalWeights)
-    VanP = pvTP(x, deg, basisp, derIdxs = derIdxs,
-                reorder = reorder, scl = scl)
-    RHP = VanP.T.conj() * Weights
-    return RHP.dot(VanP), RHP