diff --git a/rrompy/reduction_methods/base/generic_approximant.py b/rrompy/reduction_methods/base/generic_approximant.py
index e0607fe..bef162b 100644
--- a/rrompy/reduction_methods/base/generic_approximant.py
+++ b/rrompy/reduction_methods/base/generic_approximant.py
@@ -1,849 +1,884 @@
 # 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 abc import abstractmethod
 import numpy as np
 from itertools import product as iterprod
 from copy import deepcopy as copy
 from os import remove as osrm
 from rrompy.sampling import SamplingEngineStandard, SamplingEngineStandardPOD
 from rrompy.utilities.base.types import (Np1D, DictAny, HFEng, List, Tuple,
                                          ListAny, strLst, paramVal, paramList,
                                          sampList)
 from rrompy.utilities.base.data_structures import purgeDict, getNewFilename
 from rrompy.utilities.base import verbosityManager as vbMng
 from rrompy.utilities.exception_manager import (RROMPyException, RROMPyAssert,
                                                 RROMPy_READY, RROMPy_FRAGILE)
 from rrompy.utilities.base.pickle_utilities import pickleDump, pickleLoad
 from rrompy.parameter import (emptyParameterList, checkParameter,
                               checkParameterList)
 from rrompy.sampling import sampleList, emptySampleList
 
 __all__ = ['GenericApproximant']
 
 def addNormFieldToClass(self, fieldName):
     def objFunc(self, mu:paramList, *args, **kwargs) -> Np1D:
         uV = getattr(self.__class__, "get" + fieldName)(self, mu)
         kwargs["is_state"] = False
         val = self.HFEngine.norm(uV, *args, **kwargs)
         return val
     setattr(self.__class__, "norm" + fieldName, objFunc)
 
 def addNormDualFieldToClass(self, fieldName):
     def objFunc(self, mu:paramList, *args, **kwargs) -> Np1D:
         uV = getattr(self.__class__, "get" + fieldName)(self, mu)
         kwargs["is_state"] = True
         if "dual" not in kwargs.keys(): kwargs["dual"] = True
         val = self.HFEngine.norm(uV, *args, **kwargs)
         return val
     setattr(self.__class__, "norm" + fieldName, objFunc)
 
 def addPlotFieldToClass(self, fieldName):
     def objFunc(self, mu:paramList, *args, **kwargs):
         uV = getattr(self.__class__, "get" + fieldName)(self, mu)
         if "name" in kwargs.keys(): nameBase = copy(kwargs["name"])
         filesOut = []
         for j, u in enumerate(uV):
             if "name" in kwargs.keys(): kwargs["name"] = nameBase + str(j)
             filesOut += [self.HFEngine.plot(u, *args, **kwargs)]
         if "name" in kwargs.keys(): kwargs["name"] = nameBase
         if filesOut[0] is None: return None
         return filesOut
     setattr(self.__class__, "plot" + fieldName, objFunc)
 
 def addPlotDualFieldToClass(self, fieldName):
     def objFunc(self, mu:paramList, *args, **kwargs):
         uV = getattr(self.__class__, "get" + fieldName)(self, mu)
         if "name" in kwargs.keys(): nameBase = copy(kwargs["name"])
         filesOut = []
         for j, u in enumerate(uV):
             if "name" in kwargs.keys(): kwargs["name"] = nameBase + str(j)
             filesOut += [self.HFEngine.plot(u, *args, **kwargs)]
         if "name" in kwargs.keys(): kwargs["name"] = nameBase
         if filesOut[0] is None: return None
         return filesOut
     setattr(self.__class__, "plot" + fieldName, objFunc)
 
 def addOutParaviewFieldToClass(self, fieldName):
     def objFunc(self, mu:paramVal, *args, **kwargs):
         if not hasattr(self.HFEngine, "outParaview"):
             raise RROMPyException(("High fidelity engine cannot output to "
                                    "Paraview."))
         uV = getattr(self.__class__, "get" + fieldName)(self, mu)
         if "name" in kwargs.keys(): nameBase = copy(kwargs["name"])
         filesOut = []
         for j, u in enumerate(uV):
             if "name" in kwargs.keys(): kwargs["name"] = nameBase + str(j)
             filesOut += [self.HFEngine.outParaview(u, *args, **kwargs)]
         if "name" in kwargs.keys(): kwargs["name"] = nameBase
         if filesOut[0] is None: return None
         return filesOut
     setattr(self.__class__, "outParaview" + fieldName, objFunc)
 
 def addOutParaviewTimeDomainFieldToClass(self, fieldName):
     def objFunc(self, mu:paramVal, *args, **kwargs):
         if not hasattr(self.HFEngine, "outParaviewTimeDomain"):
             raise RROMPyException(("High fidelity engine cannot output to "
                                    "Paraview."))
         uV = getattr(self.__class__, "get" + fieldName)(self, mu)
         omega = args.pop(0) if len(args) > 0 else np.real(mu)
         if "name" in kwargs.keys(): nameBase = copy(kwargs["name"])
         filesOut = []
         for j, u in enumerate(uV):
             if "name" in kwargs.keys(): kwargs["name"] = nameBase + str(j)
             filesOut += [self.HFEngine.outParaviewTimeDomain(u, omega, *args,
                                                              **kwargs)]
         if "name" in kwargs.keys(): kwargs["name"] = nameBase
         if filesOut[0] is None: return None
         return filesOut
     setattr(self.__class__, "outParaviewTimeDomain" + fieldName, objFunc)
 
 class GenericApproximant:
     """
     ABSTRACT
     ROM approximant computation for parametric problems.
     
     Args:
         HFEngine: HF problem solver.
         mu0(optional): Default parameter. Defaults to 0.
         approxParameters(optional): Dictionary containing values for main
             parameters of approximant. Recognized keys are:
             - 'POD': whether to compute POD of snapshots; defaults to True;
+            - 'scaleFactorDer': scaling factors for derivative computation;
+                defaults to 'AUTO';
             - 'S': total number of samples current approximant relies upon.
             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.
         trainedModel: Trained model evaluator.
         mu0: Default parameter.
         approxParameters: Dictionary containing values for main parameters of
             approximant. Recognized keys are in parameterList{Soft,Critical}.
         parameterListSoft: Recognized keys of soft approximant parameters:
-            - 'POD': whether to compute POD of snapshots.
+            - 'POD': whether to compute POD of snapshots;
+            - 'scaleFactorDer': scaling factors for derivative computation.
         parameterListCritical: Recognized keys of critical approximant
             parameters:
             - 'S': total number of samples current approximant relies upon.
         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.
         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.
     """
 
     __all__ += [ftype + dtype for ftype, dtype in iterprod(
                       ["norm", "plot", "outParaview", "outParaviewTimeDomain"],
                       ["HF", "RHS", "Approx", "Res", "Err"])]
 
     def __init__(self, HFEngine:HFEng, mu0 : paramVal = None,
                  approxParameters : DictAny = {}, approx_state : bool = False,
                  verbosity : int = 10, timestamp : bool = True):
         self._preInit()
         self._mode = RROMPy_READY
         self.approx_state = approx_state
         self.verbosity = verbosity
         self.timestamp = timestamp
         vbMng(self, "INIT",
               "Initializing engine of type {}.".format(self.name()), 10)
         self._HFEngine = HFEngine
         self.trainedModel = None
         self.lastSolvedHF = emptyParameterList()
         self.uHF = emptySampleList()
-        self._addParametersToList(["POD"], [True], ["S"], [1])
+        self._addParametersToList(["POD", "scaleFactorDer"], [True, "AUTO"],
+                                  ["S"], [1])
         if mu0 is None:
             if hasattr(self.HFEngine, "mu0"):
                 self.mu0 = checkParameter(self.HFEngine.mu0)
             else:
                 raise RROMPyException(("Center of approximation cannot be "
                                        "inferred from HF engine. Parameter "
                                        "required"))
         else:
             self.mu0 = checkParameter(mu0, self.HFEngine.npar)
         self.resetSamples()
         self.approxParameters = approxParameters
         self._postInit()
 
         ### add norm{HF,Err} methods
         """
         Compute norm of * at arbitrary parameter.
         Args:
             mu: Target parameter.
         Returns:
             Target norm of *.
         """
         for objName in ["HF", "Err"]:
             addNormFieldToClass(self, objName)
 
         ### add norm{RHS,Res} methods
         """
         Compute norm of * at arbitrary parameter.
         Args:
             mu: Target parameter.
         Returns:
             Target norm of *.
         """
         for objName in ["RHS", "Res"]:
             addNormDualFieldToClass(self, objName)
 
         ### add plot{HF,Approx,Err} methods
         """
         Do some nice plots of * at arbitrary parameter.
 
         Args:
             mu: Target parameter.
             name(optional): Name to be shown as title of the plots. Defaults to
                 'u'.
             what(optional): Which plots to do. If list, can contain 'ABS',
                 'PHASE', 'REAL', 'IMAG'. If str, same plus wildcard 'ALL'.
                 Defaults to 'ALL'.
             save(optional): Where to save plot(s). Defaults to None, i.e. no
                 saving.
             saveFormat(optional): Format for saved plot(s). Defaults to "eps".
             saveDPI(optional): DPI for saved plot(s). Defaults to 100.
             show(optional): Whether to show figure. Defaults to True.
             figspecs(optional key args): Optional arguments for matplotlib
                 figure creation.
         """
         for objName in ["HF", "Approx", "Err"]:
             addPlotFieldToClass(self, objName)
 
         ### add plot{RHS,Res} methods
         """
         Do some nice plots of * at arbitrary parameter.
 
         Args:
             mu: Target parameter.
             name(optional): Name to be shown as title of the plots. Defaults to
                 'u'.
             what(optional): Which plots to do. If list, can contain 'ABS',
                 'PHASE', 'REAL', 'IMAG'. If str, same plus wildcard 'ALL'.
                 Defaults to 'ALL'.
             save(optional): Where to save plot(s). Defaults to None, i.e. no
                 saving.
             saveFormat(optional): Format for saved plot(s). Defaults to "eps".
             saveDPI(optional): DPI for saved plot(s). Defaults to 100.
             show(optional): Whether to show figure. Defaults to True.
             figspecs(optional key args): Optional arguments for matplotlib
                 figure creation.
         """
         for objName in ["RHS", "Res"]:
             addPlotDualFieldToClass(self, objName)
 
         ### add outParaview{HF,RHS,Approx,Res,Err} methods
         """
         Output * to ParaView file.
 
         Args:
             mu: Target parameter.
             name(optional): Base name to be used for data output.
             filename(optional): Name of output file.
             time(optional): Timestamp.
             what(optional): Which plots to do. If list, can contain 'MESH',
                 'ABS', 'PHASE', 'REAL', 'IMAG'. If str, same plus wildcard
                 'ALL'. Defaults to 'ALL'.
             forceNewFile(optional): Whether to create new output file.
             filePW(optional): Fenics File entity (for time series).
         """
         for objName in ["HF", "RHS", "Approx", "Res", "Err"]:
             addOutParaviewFieldToClass(self, objName)
 
         ### add outParaviewTimeDomain{HF,RHS,Approx,Res,Err} methods
         """
         Output * to ParaView file, converted to time domain.
 
         Args:
             mu: Target parameter.
             omega(optional): frequency.
             timeFinal(optional): final time of simulation.
             periodResolution(optional): number of time steps per period.
             name(optional): Base name to be used for data output.
             filename(optional): Name of output file.
             forceNewFile(optional): Whether to create new output file.
         """
         for objName in ["HF", "RHS", "Approx", "Res", "Err"]:
             addOutParaviewTimeDomainFieldToClass(self, objName)
 
     def _preInit(self):
         if not hasattr(self, "depth"): self.depth = 0
         else: self.depth += 1
 
     @property
     def tModelType(self):
         raise RROMPyException("No trainedModel type assigned.")
 
     def initializeModelData(self, datadict):
         from .trained_model.trained_model_data import TrainedModelData
         return (TrainedModelData(datadict["mu0"], datadict.pop("projMat"),
                                  datadict["scaleFactor"],
                                  datadict.pop("rescalingExp")),
                 ["mu0", "scaleFactor", "mus"])
 
     @property
     def parameterList(self):
         """Value of parameterListSoft + parameterListCritical."""
         return self.parameterListSoft + self.parameterListCritical
 
     def _addParametersToList(self, whatSoft : strLst = [],
                              defaultSoft : ListAny = [],
                              whatCritical : strLst = [],
                              defaultCritical : ListAny = [],
                              toBeExcluded : strLst = []):
         if not hasattr(self, "parameterToBeExcluded"):
             self.parameterToBeExcluded = []
         self.parameterToBeExcluded = toBeExcluded + self.parameterToBeExcluded
         if not hasattr(self, "parameterListSoft"):
             self.parameterListSoft = []
         if not hasattr(self, "parameterDefaultSoft"):
             self.parameterDefaultSoft = {}
         if not hasattr(self, "parameterListCritical"):
             self.parameterListCritical = []
         if not hasattr(self, "parameterDefaultCritical"):
             self.parameterDefaultCritical = {}
         for j, what in enumerate(whatSoft):
             if what not in self.parameterToBeExcluded:
                 self.parameterListSoft = [what] + self.parameterListSoft
                 self.parameterDefaultSoft[what] = defaultSoft[j]
         for j, what in enumerate(whatCritical):
             if what not in self.parameterToBeExcluded:
                 self.parameterListCritical = ([what]
                                             + self.parameterListCritical)
                 self.parameterDefaultCritical[what] = defaultCritical[j]
 
     def _postInit(self):
         if self.depth == 0:
             vbMng(self, "DEL", "Done initializing.", 10)
             del self.depth
         else: self.depth -= 1
 
     def name(self) -> str:
         return self.__class__.__name__
         
     def __str__(self) -> str:
         return self.name()
 
     def __repr__(self) -> str:
         return self.__str__() + " at " + hex(id(self))
 
     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:
             SamplingEngine = SamplingEngineStandardPOD
         else:
             SamplingEngine = SamplingEngineStandard
         self.samplingEngine = SamplingEngine(self.HFEngine,
                                              sample_state = self.approx_state,
                                              verbosity = self.verbosity)
         self.resetSamples()
         
     @property
     def HFEngine(self):
         """Value of HFEngine."""
         return self._HFEngine
     @HFEngine.setter
     def HFEngine(self, HFEngine):
         raise RROMPyException("Cannot change HFEngine.")
 
     @property
     def mu0(self):
         """Value of mu0."""
         return self._mu0
     @mu0.setter
     def mu0(self, mu0):
         mu0 = checkParameter(mu0)
         if not hasattr(self, "_mu0") or mu0 != self.mu0:
             self.resetSamples()
         self._mu0 = mu0
 
     @property
     def npar(self):
         """Number of parameters."""
         return self.mu0.shape[1]
 
     @property
     def approxParameters(self):
         """Value of approximant parameters."""
         return self._approxParameters
     @approxParameters.setter
     def approxParameters(self, approxParams):
         if not hasattr(self, "approxParameters"):
             self._approxParameters = {}
         approxParameters = purgeDict(approxParams, self.parameterList,
                                   dictname = self.name() + ".approxParameters",
                                   baselevel = 1)
         keyList = list(approxParameters.keys())
         for key in self.parameterListCritical:
             if key in keyList:
                 setattr(self, "_" + key, self.parameterDefaultCritical[key])
         for key in self.parameterListSoft:
             if key in keyList:
                 setattr(self, "_" + key, self.parameterDefaultSoft[key])
         fragile = False
         for key in self.parameterListCritical:
             if key in keyList:
                 val = approxParameters[key]
             else:
                 val = getattr(self, "_" + key, None)
                 if val is None:
                     val = self.parameterDefaultCritical[key]
             getattr(self.__class__, key, None).fset(self, val)
             fragile = fragile or val is None
         for key in self.parameterListSoft:
             if key in keyList:
                 val = approxParameters[key]
             else:
                 val = getattr(self, "_" + key, None)
                 if val is None:
                     val = self.parameterDefaultSoft[key]
             getattr(self.__class__, key, None).fset(self, val)
         if fragile:
             self._mode = RROMPy_FRAGILE
 
     @property
     def POD(self):
         """Value of POD."""
         return self._POD
     @POD.setter
     def POD(self, POD):
         if hasattr(self, "_POD"): PODold = self.POD
         else: PODold = -1
         self._POD = POD
         self._approxParameters["POD"] = self.POD
         if PODold != self.POD:
             self.samplingEngine = None
             self.resetSamples()
 
+    @property
+    def scaleFactorDer(self):
+        """Value of scaleFactorDer."""
+        if self._scaleFactorDer == "NONE": return 1.
+        if self._scaleFactorDer == "AUTO": return self.scaleFactor
+        return self._scaleFactorDer
+    @scaleFactorDer.setter
+    def scaleFactorDer(self, scaleFactorDer):
+        if hasattr(self, "_scaleFactorDer"):
+            scaleFactorDerold = self.scaleFactorDer
+        else: scaleFactorDerold = -1
+        if isinstance(scaleFactorDer, (str,)):
+            scaleFactorDer = scaleFactorDer.upper()
+        self._scaleFactorDer = scaleFactorDer
+        self._approxParameters["scaleFactorDer"] = self._scaleFactorDer
+        if scaleFactorDerold != self._scaleFactorDer:
+            self.samplingEngine = None
+            self.resetSamples()
+
+    @property
+    def scaleFactorRel(self):
+        """Value of scaleFactorDer / scaleFactor."""
+        if self._scaleFactorDer == "AUTO": return None
+        try:
+            return np.divide(self.scaleFactorDer, self.scaleFactor)
+        except:
+            raise RROMPyException(("Error in computation of relative scaling "
+                                   "factor. Make sure that scaleFactor is "
+                                   "properly initialized."))
+
     @property
     def approx_state(self):
         """Value of approx_state."""
         return self._approx_state
     @approx_state.setter
     def approx_state(self, approx_state):
         if hasattr(self, "_approx_state"): approx_stateold = self.approx_state
         else: approx_stateold = -1
         self._approx_state = approx_state
         if approx_stateold != self.approx_state:
             self.samplingEngine = None
             self.resetSamples()
 
     @property
     def S(self):
         """Value of S."""
         return self._S
     @S.setter
     def S(self, S):
         if S <= 0: raise RROMPyException("S must be positive.")
         if hasattr(self, "_S") and self._S is not None: Sold = self.S
         else: Sold = -1
         self._S = S
         self._approxParameters["S"] = self.S
         if Sold != self.S: self.resetSamples()
 
     @property
     def trainedModel(self):
         """Value of trainedModel."""
         return self._trainedModel
     @trainedModel.setter
     def trainedModel(self, trainedModel):
         self._trainedModel = trainedModel
         if self._trainedModel is not None:
             self._trainedModel.reset()
         self.lastSolvedApproxReduced = emptyParameterList()
         self.lastSolvedApprox = emptyParameterList()
         self.uApproxReduced = emptySampleList()
         self.uApprox = emptySampleList()
 
     def resetSamples(self):
         if hasattr(self, "samplingEngine") and self.samplingEngine is not None:
             self.samplingEngine.resetHistory()
         else:
             self.setupSampling()
         self._mode = RROMPy_READY
 
     def plotSamples(self, *args, **kwargs) -> List[str]:
         """
         Do some nice plots of the samples.
 
         Returns:
             Output filenames.
         """
         RROMPyAssert(self._mode, message = "Cannot plot samples.")
         return self.samplingEngine.plotSamples(*args, **kwargs)
         
     def outParaviewSamples(self, *args, **kwargs) -> List[str]:
         """
         Output samples to ParaView file.
 
         Returns:
             Output filenames.
         """
         RROMPyAssert(self._mode, message = "Cannot output samples.")
         return self.samplingEngine.outParaviewSamples(*args, **kwargs)
         
     def outParaviewTimeDomainSamples(self, *args, **kwargs) -> List[str]:
         """
         Output samples to ParaView file, converted to time domain.
 
         Returns:
             Output filenames.
         """
         RROMPyAssert(self._mode, message = "Cannot output samples.")
         return self.samplingEngine.outParaviewTimeDomainSamples(*args,
                                                                 **kwargs)
     
     def setSamples(self, samplingEngine):
         """Copy samplingEngine and samples."""
         vbMng(self, "INIT", "Transfering samples.", 10)
         self.samplingEngine = copy(samplingEngine)
         vbMng(self, "DEL", "Done transfering samples.", 10)
 
     def setTrainedModel(self, model):
         """Deepcopy approximation from trained model."""
         if hasattr(model, "storeTrainedModel"):
             verb = model.verbosity
             model.verbosity = 0
             fileOut = model.storeTrainedModel()
             model.verbosity = verb
         else:
             try:
                 fileOut = getNewFilename("trained_model", "pkl")
                 pickleDump(model.data.__dict__, fileOut)
             except:
                 raise RROMPyException(("Failed to store model data. Parameter "
                                        "model must have either "
                                        "storeTrainedModel or "
                                        "data.__dict__ properties."))
         self.loadTrainedModel(fileOut)
         osrm(fileOut)
 
     @abstractmethod
     def setupApprox(self) -> int:
         """
         Setup approximant. (ABSTRACT)
 
         Any specialization should include something like
             self.trainedModel = ...
             self.trainedModel.data = ...
             self.trainedModel.data.approxParameters = copy(
                                                          self.approxParameters)
             
         Returns > 0 if error was encountered, < 0 if no computation was
             necessary.
         """
         if self.checkComputedApprox(): return -1
         RROMPyAssert(self._mode, message = "Cannot setup approximant.")
         vbMng(self, "INIT", "Setting up {}.". format(self.name()), 5)
         pass
         vbMng(self, "DEL", "Done setting up approximant.", 5)
         return 0
 
     def checkComputedApprox(self) -> bool:
         """
         Check if setup of new approximant is not needed.
 
         Returns:
             True if new setup is not needed. False otherwise.
         """
         return self._mode == RROMPy_FRAGILE or (self.trainedModel is not None
           and self.trainedModel.data.approxParameters == self.approxParameters)
 
     def _pruneBeforeEval(self, mu:paramList, field:str, append:bool,
                          prune:bool) -> Tuple[paramList, Np1D, Np1D, bool]:
         mu = checkParameterList(mu, self.npar)[0]
         idx = np.empty(len(mu), dtype = np.int)
         if prune:
             jExtra = np.zeros(len(mu), dtype = bool)
             muExtra = emptyParameterList()
             lastSolvedMus = getattr(self, "lastSolved" + field)
             if (len(mu) > 0 and len(mu) == len(lastSolvedMus)
             and mu == lastSolvedMus):
                 idx = np.arange(len(mu), dtype = np.int)
                 return muExtra, jExtra, idx, True
             muKeep = copy(muExtra)
             for j in range(len(mu)):
                 jPos = lastSolvedMus.find(mu[j])
                 if jPos is not None:
                     idx[j] = jPos
                     muKeep.append(mu[j])
                 else:
                     jExtra[j] = True
                     muExtra.append(mu[j])
             if len(muKeep) > 0 and not append:
                 lastSolvedu = getattr(self, "u" + field)
                 idx[~jExtra] = getattr(self.__class__, "set" + field)(self,
                                      muKeep, lastSolvedu[idx[~jExtra]], append)
                 append = True
         else:
             jExtra = np.ones(len(mu), dtype = bool)
             muExtra = mu
         return muExtra, jExtra, idx, append
 
     def _setObject(self, mu:paramList, field:str, object:sampList,
                    append:bool) -> List[int]:
         newMus = checkParameterList(mu, self.npar)[0]
         newObj = sampleList(object)
         if append:
             getattr(self, "lastSolved" + field).append(newMus)
             getattr(self, "u" + field).append(newObj)
             Ltot = len(getattr(self, "u" + field))
             return list(range(Ltot - len(newObj), Ltot))
         setattr(self, "lastSolved" + field,  copy(newMus))
         setattr(self, "u" + field, copy(newObj))
         return list(range(len(getattr(self, "u" + field))))
 
     def setHF(self, muHF:paramList, uHF:sampleList,
               append : bool = False) -> List[int]:
         """Assign high fidelity solution."""
         return self._setObject(muHF, "HF", uHF, append)
 
     def evalHF(self, mu:paramList, append : bool = False,
                prune : bool = True) -> List[int]:
         """
         Find high fidelity solution with original parameters and arbitrary
             parameter.
 
         Args:
             mu: Target parameter.
             append(optional): Whether to append new HF solutions to old ones.
             prune(optional): Whether to remove duplicates of already appearing
                 HF solutions.
         """
         muExtra, jExtra, idx, append = self._pruneBeforeEval(mu, "HF", append,
                                                              prune)
         if len(muExtra) > 0:
             vbMng(self, "INIT", "Solving HF model for mu = {}.".format(mu),
                   15)
             newuHFs = self.HFEngine.solve(muExtra)
             vbMng(self, "DEL", "Done solving HF model.", 15)
             idx[jExtra] = self.setHF(muExtra, newuHFs, append)
         return list(idx)
 
     def setApproxReduced(self, muApproxR:paramList, uApproxR:sampleList,
                          append : bool = False) -> List[int]:
         """Assign high fidelity solution."""
         return self._setObject(muApproxR, "ApproxReduced", uApproxR, append)
 
     def evalApproxReduced(self, mu:paramList, append : bool = False,
                           prune : bool = True) -> List[int]:
         """
         Evaluate reduced representation of approximant at arbitrary parameter.
 
         Args:
             mu: Target parameter.
             append(optional): Whether to append new HF solutions to old ones.
             prune(optional): Whether to remove duplicates of already appearing
                 HF solutions.
         """
         self.setupApprox()
         muExtra, jExtra, idx, append = self._pruneBeforeEval(mu,
                                                              "ApproxReduced",
                                                              append, prune)
         if len(muExtra) > 0:
             newuApproxs = self.trainedModel.getApproxReduced(muExtra)
             idx[jExtra] = self.setApproxReduced(muExtra, newuApproxs, append)
         return list(idx)
 
     def setApprox(self, muApprox:paramList, uApprox:sampleList, 
                   append : bool = False) -> List[int]:
         """Assign high fidelity solution."""
         return self._setObject(muApprox, "Approx", uApprox, append)
 
     def evalApprox(self, mu:paramList, append : bool = False,
                    prune : bool = True) -> List[int]:
         """
         Evaluate approximant at arbitrary parameter.
 
         Args:
             mu: Target parameter.
             append(optional): Whether to append new HF solutions to old ones.
             prune(optional): Whether to remove duplicates of already appearing
                 HF solutions.
         """
         self.setupApprox()
         muExtra, jExtra, idx, append = self._pruneBeforeEval(mu, "Approx",
                                                              append, prune)
         if len(muExtra) > 0:
             newuApproxs = self.trainedModel.getApprox(muExtra)
             idx[jExtra] = self.setApprox(muExtra, newuApproxs, append)
         return list(idx)
 
     def getHF(self, mu:paramList, append : bool = False,
               prune : bool = True) -> sampList:
         """
         Get HF solution at arbitrary parameter.
 
         Args:
             mu: Target parameter.
 
         Returns:
             HFsolution.
         """
         mu = checkParameterList(mu, self.npar)[0]
         idx = self.evalHF(mu, append = append, prune = prune)
         return self.uHF(idx)
 
     def getRHS(self, mu:paramList) -> sampList:
         """
         Get linear system RHS at arbitrary parameter.
 
         Args:
             mu: Target parameter.
 
         Returns:
             Linear system RHS.
         """
         return self.HFEngine.residual(mu, None)
 
     def getApproxReduced(self, mu:paramList, append : bool = False,
                          prune : bool = True) -> sampList:
         """
         Get approximant at arbitrary parameter.
 
         Args:
             mu: Target parameter.
 
         Returns:
             Reduced approximant.
         """
         mu = checkParameterList(mu, self.npar)[0]
         idx = self.evalApproxReduced(mu, append = append, prune = prune)
         return self.uApproxReduced(idx)
 
     def getApprox(self, mu:paramList, append : bool = False,
                   prune : bool = True) -> sampList:
         """
         Get approximant at arbitrary parameter.
 
         Args:
             mu: Target parameter.
 
         Returns:
             Approximant.
         """
         mu = checkParameterList(mu, self.npar)[0]
         idx = self.evalApprox(mu, append = append, prune = prune)
         return self.uApprox(idx)
 
     def getRes(self, mu:paramList) -> sampList:
         """
         Get residual at arbitrary parameter.
 
         Args:
             mu: Target parameter.
 
         Returns:
             Approximant residual.
         """
         if not self.HFEngine.isCEye:
             raise RROMPyException(("Residual of solution with non-scalar C "
                                    "not computable."))
         return self.HFEngine.residual(mu, self.getApprox(mu) / self.HFEngine.C)
 
     def getErr(self, mu:paramList, append : bool = False,
                prune : bool = True) -> sampList:
         """
         Get error at arbitrary parameter.
 
         Args:
             mu: Target parameter.
 
         Returns:
             Approximant error.
         """
         return (self.getApprox(mu, append = append, prune =prune)
               - self.getHF(mu, append = append, prune = prune))
 
     def normApprox(self, mu:paramList) -> float:
         """
         Compute norm of approximant at arbitrary parameter.
 
         Args:
             mu: Target parameter.
 
         Returns:
             Target norm of approximant.
         """
         if not (self.POD and self.HFEngine.isCEye):
             return self.HFEngine.norm(self.getApprox(mu), is_state = False)
         return np.linalg.norm(self.HFEngine.C * self.getApproxReduced(mu).data,
                               axis = 0)
 
     def getPoles(self, *args, **kwargs) -> Np1D:
         """
         Obtain approximant poles.
 
         Returns:
             Numpy complex vector of poles.
         """
         self.setupApprox()
         vbMng(self, "INIT", "Computing poles of model.", 20)
         poles = self.trainedModel.getPoles(*args, **kwargs)
         vbMng(self, "DEL", "Done computing poles.", 20)
         return poles
 
     def storeTrainedModel(self, filenameBase : str = "trained_model",
                           forceNewFile : bool = True) -> str:
         """Store trained reduced model to file."""
         self.setupApprox()
         vbMng(self, "INIT", "Storing trained model to file.", 20)
         if forceNewFile:
             filename = getNewFilename(filenameBase, "pkl")
         else:
             filename = "{}.pkl".format(filenameBase)
         pickleDump(self.trainedModel.data.__dict__, filename)
         vbMng(self, "DEL", "Done storing trained model.", 20)
         return filename
 
     def loadTrainedModel(self, filename:str):
         """Load trained reduced model from file."""
         vbMng(self, "INIT", "Loading pre-trained model from file.", 20)
         datadict = pickleLoad(filename)
         self.mu0 = datadict["mu0"]
         self.scaleFactor = datadict["scaleFactor"]
         self.mus = datadict["mus"]
         trainedModel = self.tModelType()
         trainedModel.verbosity = self.verbosity
         trainedModel.timestamp = self.timestamp
         data, selfkeys = self.initializeModelData(datadict)
         for key in selfkeys: setattr(self, key, datadict.pop(key))
         approxParameters = datadict.pop("approxParameters")
         data.approxParameters = copy(approxParameters)
         for apkey in data.approxParameters.keys():
             self._approxParameters[apkey] = approxParameters.pop(apkey)
             setattr(self, "_" + apkey, self._approxParameters[apkey])
         for key in datadict: setattr(data, key, datadict[key])
         trainedModel.data = data
         self.trainedModel = trainedModel
         self._mode = RROMPy_FRAGILE
         vbMng(self, "DEL", "Done loading pre-trained model.", 20)
diff --git a/rrompy/reduction_methods/pivoted/generic_pivoted_approximant.py b/rrompy/reduction_methods/pivoted/generic_pivoted_approximant.py
index b086c2f..1eec8b0 100644
--- a/rrompy/reduction_methods/pivoted/generic_pivoted_approximant.py
+++ b/rrompy/reduction_methods/pivoted/generic_pivoted_approximant.py
@@ -1,486 +1,490 @@
 # 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.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.sampling import (SamplingEnginePivoted, SamplingEnginePivotedPOD,
                              SamplingEnginePivotedPODGlobal)
 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__ = ['GenericPivotedApproximant', 'PODGlobal']
 
 PODGlobal = 2
 
 class GenericPivotedApproximant(GenericApproximant):
     """
     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;
             - 'cutOffKind': kind of cut off strategy; available values
                 include 'SOFT' and 'HARD'; defaults to 'HARD';
             - '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;
             - 'nNearestNeighborMarginal': number of marginal nearest neighbors
                 considered if polybasisMarginal allows; 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;
             - 'cutOffKind': kind of 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;
             - 'nNearestNeighborMarginal': number of marginal nearest neighbors
                 considered if polybasisMarginal allows;
             - '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.
         cutOffKind: Kind of 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.
         nNearestNeighborMarginal: Number of marginal nearest neighbors
             considered if polybasisMarginal allows.
         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, 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(["matchingWeight", "cutOffTolerance",
                                    "cutOffKind", "polybasisMarginal",
                                    "MMarginal", "polydegreetypeMarginal",
                                    "radialDirectionalWeightsMarginal",
                                    "nNearestNeighborMarginal",
                                    "interpRcondMarginal"],
                                   [1, np.inf, "HARD", "MONOMIAL", "AUTO",
                                    "TOTAL", [1], -1, -1], ["samplerPivot",
                                    "SMarginal", "samplerMarginal"],
                                   [QSBase, [1], QSBase])
         del QS
         self._directionPivot = directionPivot
         super().__init__(*args, **kwargs)
         self._postInit()
 
     @property
     def tModelType(self):
         from .trained_model.trained_model_pivoted import TrainedModelPivoted
         return TrainedModelPivoted
 
     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 matchingWeight(self):
         """Value of matchingWeight."""
         return self._matchingWeight
     @matchingWeight.setter
     def matchingWeight(self, matchingWeight):
         self._matchingWeight = matchingWeight
         self._approxParameters["matchingWeight"] = self.matchingWeight
 
     @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 cutOffKind(self):
         """Value of cutOffKind."""
         return self._cutOffKind
     @cutOffKind.setter
     def cutOffKind(self, cutOffKind):
         cutOffKind = cutOffKind.upper()
         if cutOffKind not in ["SOFT", "HARD"]:
             RROMPyWarning(("Cut off kind not recognized. Overriding to "
                            "'HARD'."))
             cutOffKind = "HARD"
         self._cutOffKind = cutOffKind
         self._approxParameters["cutOffKind"] = self.cutOffKind
 
     @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 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: 
                 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)
 
     @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 nNearestNeighborMarginal(self):
         """Value of nNearestNeighborMarginal."""
         return self._nNearestNeighborMarginal
     @nNearestNeighborMarginal.setter
     def nNearestNeighborMarginal(self, nNearestNeighborMarginal):
         self._nNearestNeighborMarginal = nNearestNeighborMarginal
         self._approxParameters["nNearestNeighborMarginal"] = (
                                                  self.nNearestNeighborMarginal)
 
     @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.__str__()
         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.__str__())
         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 _finalizeMarginalization(self):
         vbMng(self, "INIT", "Recompressing by cut off.", 10)
         msg = self.trainedModel.recompressByCutOff(
                                          self.cutOffTolerance, self.cutOffKind,
                                          self.samplerPivot.normalFoci(),
                                          self.samplerPivot.groundPotential())
         vbMng(self, "DEL", "Done recompressing." + msg, 10)
         interpPars = [self.verbosity >= 5,
                       self.polydegreetypeMarginal == "TOTAL", {}]
         if self.polybasisMarginal not in ppb:
             interpPars[-1]["nNearestNeighbor"] = self.nNearestNeighborMarginal
         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)
 
     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
     
diff --git a/rrompy/reduction_methods/pivoted/greedy/generic_pivoted_greedy_approximant.py b/rrompy/reduction_methods/pivoted/greedy/generic_pivoted_greedy_approximant.py
index 51aa91e..fdbab32 100644
--- a/rrompy/reduction_methods/pivoted/greedy/generic_pivoted_greedy_approximant.py
+++ b/rrompy/reduction_methods/pivoted/greedy/generic_pivoted_greedy_approximant.py
@@ -1,461 +1,484 @@
 # 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 abc import abstractmethod
 from copy import deepcopy as copy
 import numpy as np
 from matplotlib import pyplot as plt
 from rrompy.reduction_methods.pivoted.generic_pivoted_approximant import (
                                           GenericPivotedApproximant, PODGlobal)
 from rrompy.utilities.base.types import Np1D, Tuple, List, paramVal, paramList
 from rrompy.utilities.base import verbosityManager as vbMng
 from rrompy.utilities.numerical.point_matching import (pointMatching,
                                               chordalMetricAdjusted, potential)
 from rrompy.utilities.exception_manager import (RROMPyException, RROMPyAssert,
                                                 RROMPyWarning)
 from rrompy.parameter import checkParameterList, emptyParameterList
 
 __all__ = ['GenericPivotedGreedyApproximant']
 
 class GenericPivotedGreedyApproximant(GenericPivotedApproximant):
     """
     ROM pivoted greedy interpolant 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;
             - 'cutOffKind': kind of cut off strategy; available values
                 include 'SOFT' and 'HARD'; defaults to 'HARD';
             - 'matchingWeightError': weight for pole matching optimization in
                 error estimation; defaults to 0;
             - 'cutOffToleranceError': tolerance for ignoring parasitic poles
                 in error estimation; defaults to 'AUTO', i.e. cutOffTolerance;
             - 'S': total number of pivot samples current approximant relies
                 upon;
             - 'samplerPivot': pivot sample point generator;
             - 'SMarginal': number of starting marginal samples;
             - 'samplerMarginalGrid': marginal sample point generator via sparse
                 grid;
             - '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;
             - 'greedyTolMarginal': uniform error tolerance for marginal greedy
                algorithm; defaults to 1e-1;
             - 'maxIterMarginal': maximum number of marginal greedy steps;
                defaults to 1e2;
             - 'polydegreetypeMarginal': type of polynomial degree for marginal;
                 defaults to 'TOTAL';
             - 'radialDirectionalWeightsMarginal': radial basis weights for
                 marginal interpolant; defaults to 1;
             - 'nNearestNeighborMarginal': number of marginal nearest neighbors
                 considered if polybasisMarginal allows; 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;
             - 'cutOffKind': kind of cut off strategy;
             - 'matchingWeightError': weight for pole matching optimization in
                 error estimation;
             - 'cutOffToleranceError': tolerance for ignoring parasitic poles
                 in error estimation;
             - 'polybasisMarginal': type of polynomial basis for marginal
                 interpolation;
             - 'MMarginal': degree of marginal interpolant;
             - 'greedyTolMarginal': uniform error tolerance for marginal greedy
                 algorithm;
             - 'maxIterMarginal': maximum number of marginal greedy steps;
             - 'polydegreetypeMarginal': type of polynomial degree for marginal;
             - 'radialDirectionalWeightsMarginal': radial basis weights for
                 marginal interpolant;
             - 'nNearestNeighborMarginal': number of marginal nearest neighbors
                 considered if polybasisMarginal allows;
             - '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;
             - 'samplerMarginalGrid': marginal sample point generator via sparse
                 grid.
         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.
         cutOffKind: Kind of cut off strategy.
         matchingWeightError: Weight for pole matching optimization in error
             estimation.
         cutOffToleranceError: Tolerance for ignoring parasitic poles in error
             estimation.
         S: Total number of pivot samples current approximant relies upon.
         samplerPivot: Pivot sample point generator.
         SMarginal: Total number of marginal samples current approximant relies
             upon.
         samplerMarginalGrid: Marginal sample point generator via sparse grid.
         polybasisMarginal: Type of polynomial basis for marginal interpolation.
         MMarginal: Degree of marginal interpolant.
         greedyTolMarginal: Uniform error tolerance for marginal greedy
             algorithm.
         maxIterMarginal: Maximum number of marginal greedy steps.
         polydegreetypeMarginal: Type of polynomial degree for marginal.
         radialDirectionalWeightsMarginal: Radial basis weights for marginal
             interpolant.
         nNearestNeighborMarginal: Number of marginal nearest neighbors
             considered if polybasisMarginal allows.
         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()
         from rrompy.parameter import localSparseGrid as SG
         SGBase = SG([[0.], [1.]], "UNIFORM")
         self._addParametersToList(["matchingWeightError",
                                    "cutOffToleranceError", "greedyTolMarginal",
                                    "maxIterMarginal"], [0., "AUTO", 1e-1, 1e2],
                                   ["samplerMarginalGrid"], [SGBase],
                                   toBeExcluded = ["samplerMarginal"])
         super().__init__(*args, **kwargs)
         self._postInit()
 
     @property
     def muBoundsMarginal(self):
         """Value of muBoundsMarginal."""
         return self.samplerMarginalGrid.lims
 
+    @property
+    def scaleFactorDer(self):
+        """Value of scaleFactorDer."""
+        if self._scaleFactorDer == "NONE": return 1.
+        if self._scaleFactorDer == "AUTO": return self._scaleFactorOldPivot
+        return self._scaleFactorDer
+    @scaleFactorDer.setter
+    def scaleFactorDer(self, scaleFactorDer):
+        if hasattr(self, "_scaleFactorDer"):
+            scaleFactorDerold = self.scaleFactorDer
+        else: scaleFactorDerold = -1
+        if isinstance(scaleFactorDer, (str,)):
+            scaleFactorDer = scaleFactorDer.upper()
+        self._scaleFactorDer = scaleFactorDer
+        self._approxParameters["scaleFactorDer"] = self._scaleFactorDer
+        if scaleFactorDerold != self._scaleFactorDer:
+            self.samplingEngine = None
+            self.resetSamples()
+
     @property
     def samplerMarginalGrid(self):
         """Value of samplerMarginalGrid."""
         return self._samplerMarginalGrid
     @samplerMarginalGrid.setter
     def samplerMarginalGrid(self, samplerMarginalGrid):
         if 'refine' not in dir(samplerMarginalGrid):
             raise RROMPyException("Marginal sampler type not recognized.")
         if (hasattr(self, '_samplerMarginalGrid')
         and self._samplerMarginalGrid is not None):
             samplerOld = self.samplerMarginalGrid
         self._samplerMarginalGrid = samplerMarginalGrid
         self._approxParameters["samplerMarginalGrid"] = (
                                             self.samplerMarginalGrid.__str__())
         if (not 'samplerOld' in locals()
          or samplerOld != self.samplerMarginalGrid):
             self.resetSamples()
     
     @property
     def matchingWeightError(self):
         """Value of matchingWeightError."""
         return self._matchingWeightError
     @matchingWeightError.setter
     def matchingWeightError(self, matchingWeightError):
         self._matchingWeightError = matchingWeightError
         self._approxParameters["matchingWeightError"] = (
                                                       self.matchingWeightError)
 
     @property
     def cutOffToleranceError(self):
         """Value of cutOffToleranceError."""
         return self._cutOffToleranceError
     @cutOffToleranceError.setter
     def cutOffToleranceError(self, cutOffToleranceError):
         if isinstance(cutOffToleranceError, (str,)):
             cutOffToleranceError = cutOffToleranceError.upper()\
                                                        .strip().replace(" ","")
             if cutOffToleranceError != "AUTO":
                 RROMPyWarning(("String value of cutOffToleranceError not "
                                "recognized. Overriding to 'AUTO'."))
                 cutOffToleranceError == "AUTO"
         self._cutOffToleranceError = cutOffToleranceError
         self._approxParameters["cutOffToleranceError"] = (
                                                      self.cutOffToleranceError)
 
     @property
     def greedyTolMarginal(self):
         """Value of greedyTolMarginal."""
         return self._greedyTolMarginal
     @greedyTolMarginal.setter
     def greedyTolMarginal(self, greedyTolMarginal):
         if greedyTolMarginal < 0:
             raise RROMPyException("greedyTolMarginal must be non-negative.")
         if (hasattr(self, "_greedyTolMarginal")
         and self.greedyTolMarginal is not None):
             greedyTolMarginalold = self.greedyTolMarginal
         else:
             greedyTolMarginalold = -1
         self._greedyTolMarginal = greedyTolMarginal
         self._approxParameters["greedyTolMarginal"] = self.greedyTolMarginal
         if greedyTolMarginalold != self.greedyTolMarginal:
             self.resetSamples()
 
     @property
     def maxIterMarginal(self):
         """Value of maxIterMarginal."""
         return self._maxIterMarginal
     @maxIterMarginal.setter
     def maxIterMarginal(self, maxIterMarginal):
         if maxIterMarginal <= 0:
             raise RROMPyException("maxIterMarginal must be positive.")
         if (hasattr(self, "_maxIterMarginal")
         and self.maxIterMarginal is not None):
             maxIterMarginalold = self.maxIterMarginal
         else:
             maxIterMarginalold = -1
         self._maxIterMarginal = maxIterMarginal
         self._approxParameters["maxIterMarginal"] = self.maxIterMarginal
         if maxIterMarginalold != self.maxIterMarginal:
             self.resetSamples()
 
     def resetSamples(self):
         """Reset samples."""
         super().resetSamples()
         if not hasattr(self, "_temporaryPivot"):
             self._mus = emptyParameterList()
             self.musMarginal = emptyParameterList()
             if hasattr(self, "samplerMarginalGrid"):
                 self.samplerMarginalGrid.reset()
         if hasattr(self, "samplingEngine") and self.samplingEngine is not None:
             self.samplingEngine.resetHistory()
 
     def errorEstimatorMarginal(self, return_max : bool = False) -> Np1D:
         vbMng(self, "INIT", "Matching poles.", 10)
         self.trainedModel.initializeFromRational(
                                       self.HFEngine, self.matchingWeight,
                                       self.POD == PODGlobal, self.approx_state)
         vbMng(self, "DEL", "Done matching poles.", 10)
         self._finalizeMarginalization()
         _tMdataFull = copy(self.trainedModel.data)
         vbMng(self.trainedModel, "INIT",
               "Evaluating error estimator at mu = {}.".format(
                                        self.trainedModel.data.musMarginal), 10)
         err = np.zeros(len(self.trainedModel.data.musMarginal))
         if len(err) <= 1: err[:] = np.inf
         else:
             if self.cutOffToleranceError == "AUTO":
                 cutOffTolErr = self.cutOffTolerance
             else:
                 cutOffTolErr = self.cutOffToleranceError
             if not hasattr(self, "_MMarginal_isauto"):
                 if not hasattr(self, "_MMarginalOriginal"):
                     self._MMarginalOriginal = self.MMarginal
                 self.MMarginal = self._MMarginalOriginal
             _musMExcl = None
             self.verbosity -= 35
             self.trainedModel.verbosity -= 35
             foci = self.samplerPivot.normalFoci()
             ground = self.samplerPivot.groundPotential()
             for j in range(len(err)):
                 jEff = j - (j > 0)
                 muTest = self.trainedModel.data.musMarginal[jEff]
                 polesEx = self.trainedModel.data.HIs[jEff].poles
                 idxExEff = np.where(potential(polesEx, foci) - ground
                                                    <= cutOffTolErr * ground)[0]
                 polesEx = polesEx[idxExEff]
                 if self.matchingWeightError != 0:
                     resEx = self.trainedModel.data.HIs[jEff].coeffs[idxExEff]
                 else:
                     resEx = None
                 if j > 0: self.musMarginal.insert(_musMExcl, j - 1)
                 _musMExcl = self.musMarginal[j]
                 self.musMarginal.pop(j)
                 if len(polesEx) == 0: continue
                 self.trainedModel.updateEffectiveSamples(
                                       self.HFEngine, [j], self.matchingWeight,
                                       self.POD == PODGlobal, self.approx_state)
                 self._reduceDegreeNNoWarn = 1
                 self._finalizeMarginalization()
                 polesAp = self.trainedModel.interpolateMarginalPoles(muTest)[
                                                                         ..., 0]
                 idxApEff = np.where(potential(polesAp, foci) - ground
                                                    <= cutOffTolErr * ground)[0]
                 polesAp = polesAp[idxApEff]
                 if self.matchingWeightError != 0:
                     resAp = self.trainedModel.interpolateMarginalCoeffs(
                                                         muTest)[idxApEff, :, 0]
                     if self.POD != PODGlobal:
                         resEx = self.trainedModel.data.projMat[:,
                                                  : resEx.shape[1]].dot(resEx.T)
                         resAp = self.trainedModel.data.projMat[:,
                                                  : resAp.shape[1]].dot(resAp.T)
                 else:
                     resAp = None
                 dist = chordalMetricAdjusted(
                                 polesEx, polesAp, self.matchingWeightError,
                                 resEx, resAp, self.HFEngine, self.approx_state)
                 pmR, pmC = pointMatching(dist)
                 err[j] = np.mean(dist[pmR, pmC])
             self.trainedModel.updateEffectiveSamples(self.HFEngine, None,
                  self.matchingWeight, self.POD == PODGlobal, self.approx_state)
             if not hasattr(self, "_MMarginal_isauto"):
                 self.MMarginal = self._MMarginalOriginal
             self.musMarginal.append(_musMExcl)
             self.verbosity += 35
             self.trainedModel.verbosity += 35
         self.trainedModel.data = _tMdataFull
         del self._reduceDegreeNNoWarn
         vbMng(self.trainedModel, "DEL", "Done evaluating error estimator", 10)
         if not return_max: return err
         idxMaxEst = np.where(err > self.greedyTolMarginal)[0]
         return err, idxMaxEst, err[idxMaxEst]
 
     def plotEstimatorMarginal(self, est:Np1D, idxMax:List[int],
                               estMax:List[float]):
         if not (np.any(np.isnan(est)) or np.any(np.isinf(est))):
             fig = plt.figure(figsize = plt.figaspect(1. / self.nparMarginal))
             for jpar in range(self.nparMarginal):
                 ax = fig.add_subplot(1, self.nparMarginal, 1 + jpar)
                 musre = copy(self.trainedModel.data.musMarginal.re.data)
                 errCP = copy(est)
                 idx = np.delete(np.arange(self.nparMarginal), jpar)
                 while len(musre) > 0:
                     if self.nparMarginal == 1:
                         currIdx = np.arange(len(musre))
                     else:
                         currIdx = np.where(np.isclose(np.sum(
                              np.abs(musre[:, idx] - musre[0, idx]), 1), 0.))[0]
                     currIdxSorted = currIdx[np.argsort(musre[currIdx, jpar])]
                     ax.semilogy(musre[currIdxSorted, jpar],
                                 errCP[currIdxSorted], 'k.-', linewidth = 1)
                     musre = np.delete(musre, currIdx, 0)
                     errCP = np.delete(errCP, currIdx)
                 ax.semilogy(self.musMarginal.re(jpar),
                             (self.greedyTolMarginal,) * len(self.musMarginal),
                             '*m')
                 if len(idxMax) > 0 and estMax is not None:
                     ax.semilogy(self.trainedModel.data.musMarginal.re(
                                                    idxMax, jpar), estMax, 'xr')
                 ax.grid()
             plt.tight_layout()
             plt.show()
 
     def _addMarginalSample(self, mus:paramList):
         mus = checkParameterList(mus, self.nparMarginal)[0]
         if len(mus) == 0: return
         nmus = len(mus)
         vbMng(self, "MAIN",
               ("Adding marginal sample point{} no. {}{} at {} to training "
                "set.").format("s" * (nmus > 1), len(self.musMarginal) + 1,
                       "--{}".format(len(self.musMarginal) + nmus) * (nmus > 1),
                       mus), 3)
         self.musMarginal.append(mus)
         self.setupApproxPivoted(mus)
         self._SMarginal = len(self.musMarginal)
         self._approxParameters["SMarginal"] = self.SMarginal
 
     def greedyNextSampleMarginal(self, muidx:int, plotEst : str = "NONE") \
                                           -> Tuple[Np1D, int, float, paramVal]:
         RROMPyAssert(self._mode, message = "Cannot add greedy sample.")
         idxAdded = self.samplerMarginalGrid.refine(muidx)
         self._addMarginalSample(self.samplerMarginalGrid.points[idxAdded])
         errorEstTest, muidx, maxErrorEst = self.errorEstimatorMarginal(True)
         if plotEst == "ALL":
             self.plotEstimatorMarginal(errorEstTest, muidx, maxErrorEst)
         return (errorEstTest, muidx, maxErrorEst,
                 self.samplerMarginalGrid.points[muidx])
                                               
     def _preliminaryTrainingMarginal(self):
         """Initialize starting snapshots of solution map."""
         RROMPyAssert(self._mode, message = "Cannot start greedy algorithm.")
         if np.sum(self.samplingEngine.nsamples) > 0: return
         self.resetSamples()
         idx = [0]
         while self.samplerMarginalGrid.npoints < self.SMarginal:
             idx = self.samplerMarginalGrid.refine(idx)
         self._addMarginalSample(self.samplerMarginalGrid.points)
 
     @abstractmethod
     def setupApproxPivoted(self, mu:paramVal) -> int:
         if self.checkComputedApproxPivoted(): return -1
         RROMPyAssert(self._mode, message = "Cannot setup approximant.")
         vbMng(self, "INIT", "Setting up pivoted approximant.", 10)
         pass
         vbMng(self, "DEL", "Done setting up pivoted approximant.", 10)
         return 0
 
     def setupApprox(self, plotEst : str = "NONE") -> int:
         """Compute greedy snapshots of solution map."""
         if self.checkComputedApprox(): return -1
         RROMPyAssert(self._mode, message = "Cannot start greedy algorithm.")
         vbMng(self, "INIT", "Starting computation of snapshots.", 3)
         self._preliminaryTrainingMarginal()
         muidx, max2ErrorEst, firstGreedyIter = [], np.inf, True
         while firstGreedyIter or (max2ErrorEst > self.greedyTolMarginal
                   and self.samplerMarginalGrid.npoints < self.maxIterMarginal):
             errorEstTest, muidx, maxErrorEst, mu = \
                                   self.greedyNextSampleMarginal(muidx, plotEst)
             if len(maxErrorEst) > 0:
                 max2ErrorEst = np.max(maxErrorEst)
                 vbMng(self, "MAIN", ("Uniform testing error estimate "
                                      "{:.4e}.").format(max2ErrorEst), 3)
             else:
                 max2ErrorEst = 0.
             firstGreedyIter = False
         if plotEst == "LAST":
             self.plotEstimatorMarginal(errorEstTest, muidx, maxErrorEst)
         vbMng(self, "DEL", 
               ("Done computing snapshots (final snapshot count: "
                "{}).").format(np.sum(self.samplingEngine.nsamples)), 3)
         return 0
 
     def checkComputedApprox(self) -> bool:
         return (super().checkComputedApprox()
             and len(self.mus) == len(self.trainedModel.data.mus))
 
     def checkComputedApproxPivoted(self) -> bool:
         return (super().checkComputedApprox()
           and len(self.musMarginal) == len(self.trainedModel.data.musMarginal))
diff --git a/rrompy/reduction_methods/pivoted/greedy/rational_interpolant_greedy_pivoted_greedy.py b/rrompy/reduction_methods/pivoted/greedy/rational_interpolant_greedy_pivoted_greedy.py
index 07bf27b..246232d 100644
--- a/rrompy/reduction_methods/pivoted/greedy/rational_interpolant_greedy_pivoted_greedy.py
+++ b/rrompy/reduction_methods/pivoted/greedy/rational_interpolant_greedy_pivoted_greedy.py
@@ -1,375 +1,379 @@
 #Copyright (C) 2018 by the RROMPy authors
 #
 # This file is part of RROMPy.
 #
 # RROMPy is free software: you can redistribute it and/or modify
 # it under the terms of the GNU Lesser General Public License as published by
 # the Free Software Foundation, either version 3 of the License, or
 # (at your option) any later version.
 #
 # RROMPy is distributed in the hope that it will be useful,
 # but WITHOUT ANY WARRANTY; without even the implied warranty of
 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
 # GNU Lesser General Public License for more details.
 #
 # You should have received a copy of the GNU Lesser General Public License
 # along with RROMPy. If not, see <http://www.gnu.org/licenses/>.
 #
 
 from copy import deepcopy as copy
 import numpy as np
 from .generic_pivoted_greedy_approximant import GenericPivotedGreedyApproximant
 from rrompy.utilities.numerical import dot
 from rrompy.reduction_methods.standard.greedy import RationalInterpolantGreedy
 from rrompy.reduction_methods.standard.greedy.generic_greedy_approximant \
                                                             import pruneSamples
 from rrompy.reduction_methods.pivoted import  RationalInterpolantGreedyPivoted
 from rrompy.reduction_methods.pivoted.generic_pivoted_approximant import (
                                                                      PODGlobal)
 from rrompy.utilities.base.types import Np1D, Tuple, paramVal, paramList
 from rrompy.utilities.base import verbosityManager as vbMng
 from rrompy.utilities.exception_manager import RROMPyAssert
 from rrompy.parameter import emptyParameterList
 
 __all__ = ['RationalInterpolantGreedyPivotedGreedy']
 
 class RationalInterpolantGreedyPivotedGreedy(GenericPivotedGreedyApproximant,
                                              RationalInterpolantGreedyPivoted):
     """
     ROM greedy pivoted greedy rational interpolant computation for parametric
         problems.
 
     Args:
         HFEngine: HF problem solver.
         mu0(optional): Default parameter. Defaults to 0.
         directionPivot(optional): Pivot components. Defaults to [0].
         approxParameters(optional): Dictionary containing values for main
             parameters of approximant. Recognized keys are:
             - 'POD': 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;
             - 'cutOffKind': kind of cut off strategy; available values
                 include 'SOFT' and 'HARD'; defaults to 'HARD';
             - 'matchingWeightError': weight for pole matching optimization in
                 error estimation; defaults to 0;
             - 'cutOffToleranceError': tolerance for ignoring parasitic poles
                 in error estimation; defaults to 'AUTO', i.e. cutOffTolerance;
             - 'S': total number of pivot samples current approximant relies
                 upon;
             - 'samplerPivot': pivot sample point generator;
             - 'SMarginal': number of starting marginal samples;
             - 'samplerMarginalGrid': marginal sample point generator via sparse
                 grid;
             - 'polybasis': type of polynomial basis for pivot interpolation;
                 defaults to 'MONOMIAL';
             - 'polybasisMarginal': type of polynomial basis for marginal
                 interpolation; allowed values include 'MONOMIAL', 'CHEBYSHEV'
                 and 'LEGENDRE'; defaults to 'MONOMIAL';
             - 'greedyTol': uniform error tolerance for greedy algorithm;
                 defaults to 1e-2;
             - 'collinearityTol': collinearity tolerance for greedy algorithm;
                 defaults to 0.;
             - 'maxIter': maximum number of greedy steps; defaults to 1e2;
             - 'nTestPoints': number of test points; defaults to 5e2;
             - 'trainSetGenerator': training sample points generator; defaults
                 to sampler;
             - 'errorEstimatorKind': kind of error estimator; available values
                 include 'AFFINE', 'DISCREPANCY', 'LOOK_AHEAD',
                 'LOOK_AHEAD_RES', and 'NONE'; defaults to 'NONE';
             - 'MMarginal': degree of marginal interpolant; defaults to 'AUTO',
                 i.e. maximum allowed;
             - 'greedyTolMarginal': uniform error tolerance for marginal greedy
                algorithm; defaults to 1e-1;
             - 'maxIterMarginal': maximum number of marginal greedy steps;
                defaults to 1e2;
             - 'polydegreetypeMarginal': type of polynomial degree for marginal;
                 defaults to 'TOTAL';
             - 'radialDirectionalWeightsMarginal': radial basis weights for
                 marginal interpolant; defaults to 1;
             - 'nNearestNeighborMarginal': number of marginal nearest neighbors
                 considered if polybasisMarginal allows; defaults to -1;
             - 'interpRcond': tolerance for pivot interpolation; defaults to
                 None;
             - 'interpRcondMarginal': tolerance for marginal 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.
         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;
             - 'cutOffKind': kind of cut off strategy;
             - 'matchingWeightError': weight for pole matching optimization in
                 error estimation;
             - 'cutOffToleranceError': tolerance for ignoring parasitic poles
                 in error estimation;
             - 'polybasis': type of polynomial basis for pivot interpolation;
             - 'polybasisMarginal': type of polynomial basis for marginal
                 interpolation;
             - 'greedyTol': uniform error tolerance for greedy algorithm;
             - 'collinearityTol': collinearity tolerance for greedy algorithm;
             - 'maxIter': maximum number of greedy steps;
             - 'nTestPoints': number of test points;
             - 'trainSetGenerator': training sample points generator;
             - 'errorEstimatorKind': kind of error estimator;
             - 'MMarginal': degree of marginal interpolant;
             - 'greedyTolMarginal': uniform error tolerance for marginal greedy
                 algorithm;
             - 'maxIterMarginal': maximum number of marginal greedy steps;
             - 'polydegreetypeMarginal': type of polynomial degree for marginal;
             - 'radialDirectionalWeightsMarginal': radial basis weights for
                 marginal interpolant;
             - 'nNearestNeighborMarginal': number of marginal nearest neighbors
                 considered if polybasisMarginal allows;
             - 'interpRcond': tolerance for pivot interpolation;
             - 'interpRcondMarginal': tolerance for marginal interpolation;
             - 'robustTol': tolerance for robust rational denominator
                 management.
         parameterListCritical: Recognized keys of critical approximant
             parameters:
             - 'S': total number of pivot samples current approximant relies
                 upon;
             - 'samplerPivot': pivot sample point generator;
             - 'SMarginal': total number of marginal samples current approximant
                 relies upon;
             - 'samplerMarginalGrid': marginal sample point generator via sparse
                 grid.
         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.
         cutOffKind: Kind of cut off strategy.
         matchingWeightError: Weight for pole matching optimization in error
             estimation.
         cutOffToleranceError: Tolerance for ignoring parasitic poles in error
             estimation.
         S: Total number of pivot samples current approximant relies upon.
         samplerPivot: Pivot sample point generator.
         SMarginal: Total number of marginal samples current approximant relies
             upon.
         samplerMarginalGrid: Marginal sample point generator via sparse grid.
         polybasis: Type of polynomial basis for pivot interpolation.
         polybasisMarginal: Type of polynomial basis for marginal interpolation.
         greedyTol: uniform error tolerance for greedy algorithm.
         collinearityTol: Collinearity tolerance for greedy algorithm.
         maxIter: maximum number of greedy steps.
         nTestPoints: number of starting training points.
         trainSetGenerator: training sample points generator.
         errorEstimatorKind: kind of error estimator.
         MMarginal: Degree of marginal interpolant.
         greedyTolMarginal: Uniform error tolerance for marginal greedy
             algorithm.
         maxIterMarginal: Maximum number of marginal greedy steps.
         polydegreetypeMarginal: Type of polynomial degree for marginal.
         radialDirectionalWeightsMarginal: Radial basis weights for marginal
             interpolant.
         nNearestNeighborMarginal: Number of marginal nearest neighbors
             considered if polybasisMarginal allows.
         interpRcond: Tolerance for pivot interpolation.
         interpRcondMarginal: Tolerance for marginal interpolation.
         robustTol: Tolerance for robust rational denominator management.
         muBounds: list of bounds for pivot parameter values.
         muBoundsMarginal: list of bounds for marginal parameter values.
         samplingEngine: Sampling engine.
         uHF: High fidelity solution(s) with parameter(s) lastSolvedHF as
             sampleList.
         lastSolvedHF: Parameter(s) corresponding to last computed high fidelity
             solution(s) as parameterList.
         uApproxReduced: Reduced approximate solution(s) with parameter(s)
             lastSolvedApprox as sampleList.
         lastSolvedApproxReduced: Parameter(s) corresponding to last computed
             reduced approximate solution(s) as parameterList.
         uApprox: Approximate solution(s) with parameter(s) lastSolvedApprox as
             sampleList.
         lastSolvedApprox: Parameter(s) corresponding to last computed
             approximate solution(s) as parameterList.
     """
     
     @property
     def sampleBatchSize(self):
         """Value of sampleBatchSize."""
         return 1
 
     @property
     def sampleBatchIdx(self):
         """Value of sampleBatchIdx."""
         return self.S
 
     def _finalizeSnapshots(self):
         self.samplingEngine = self._samplingEngineOld
         for muM, sEN in zip(self.musMargLoc, self.samplingEngs):
             self.samplingEngine.samples += [sEN.samples]
             self.samplingEngine.nsamples += [sEN.nsamples]
             self.samplingEngine.mus += [sEN.mus]
             self.samplingEngine.musMarginal.append(muM)
             self.samplingEngine._derIdxs += [[(0,) * self.npar] 
                                                   for _ in range(sEN.nsamples)]
             if self.POD:
                 self.samplingEngine.RPOD += [sEN.RPOD]
                 self.samplingEngine.samples_full += [copy(sEN.samples_full)]
         if self.POD == PODGlobal:
             self.samplingEngine.coalesceSamples(self.interpRcondMarginal)
         else:
             self.samplingEngine.coalesceSamples()
 
     def greedyNextSample(self, muidx:int, plotEst : str = "NONE")\
                                           -> Tuple[Np1D, int, float, paramVal]:
         """Compute next greedy snapshot of solution map."""
         RROMPyAssert(self._mode, message = "Cannot add greedy sample.")
         mus = copy(self.muTest[muidx])
         self.muTest.pop(muidx)
         for j, mu in enumerate(mus):
             vbMng(self, "MAIN",
                   ("Adding sample point no. {} at {} to training "
                    "set.").format(len(self.mus) + 1, mu), 3)
             self.mus.append(mu)
             self._S = len(self.mus)
             self._approxParameters["S"] = self.S
             if (self.samplingEngine.nsamples <= len(mus) - j - 1
              or not np.allclose(mu,
                                 self.samplingEngine.mus.data[j - len(mus)])):
                 self.samplingEngine.nextSample(mu)
             if self._isLastSampleCollinear():
                 vbMng(self, "MAIN",
                       ("Collinearity above tolerance detected. Starting "
                        "preemptive greedy loop termination."), 3)
                 self._collinearityFlag = 1
                 errorEstTest = np.empty(len(self.muTest))
                 errorEstTest[:] = np.nan
                 return errorEstTest, [-1], np.nan, np.nan
         errorEstTest, muidx, maxErrorEst = self.errorEstimator(self.muTest,
                                                                True)
         if plotEst == "ALL":
             self.plotEstimator(errorEstTest, muidx, maxErrorEst)
         return errorEstTest, muidx, maxErrorEst, self.muTest[muidx]
 
     def _preliminaryTraining(self):
         """Initialize starting snapshots of solution map."""
         RROMPyAssert(self._mode, message = "Cannot start greedy algorithm.")
         if self.samplingEngine.nsamples > 0: return
         self.resetSamples()
-        self.samplingEngine.scaleFactor = self._scaleFactorOldPivot
+        self.samplingEngine.scaleFactor = self.scaleFactorDer
         musPivot = self.trainSetGenerator.generatePoints(self.S)
         while len(musPivot) > self.S: musPivot.pop()
         muTestBasePivot = self.samplerPivot.generatePoints(self.nTestPoints,
                                                            False)
         idxPop = pruneSamples(
          muTestBasePivot ** self.HFEngine.rescalingExp[self.directionPivot[0]],
          musPivot ** self.HFEngine.rescalingExp[self.directionPivot[0]],
          1e-10 * self.scaleFactor[0])
         muTestBasePivot.pop(idxPop)
         self.mus = emptyParameterList()
         self.mus.reset((self.S - 1, self.HFEngine.npar))
         self.muTest = emptyParameterList()
         self.muTest.reset((len(muTestBasePivot) + 1, self.HFEngine.npar))
         for k in range(self.S - 1):
             self.mus.data[k, self.directionPivot] = musPivot[k].data
             self.mus.data[k, self.directionMarginal] = self.musMargLoc[-1].data
         for k in range(len(muTestBasePivot)):
             self.muTest.data[k, self.directionPivot] = muTestBasePivot[k].data
             self.muTest.data[k, self.directionMarginal] = (
                                                       self.musMargLoc[-1].data)
         self.muTest.data[-1, self.directionPivot] = musPivot[-1].data
         self.muTest.data[-1, self.directionMarginal] = self.musMargLoc[-1].data
         if len(self.mus) > 0:
             vbMng(self, "MAIN", 
                   ("Adding first {} sample point{} at {} to training "
                    "set.").format(self.S - 1, "" + "s" * (self.S > 2),
                                   self.mus), 3)
             self.samplingEngine.iterSample(self.mus)
         self._S = len(self.mus)
         self._approxParameters["S"] = self.S
         self.M, self.N = ("AUTO",) * 2
 
     def setupApproxPivoted(self, mus:paramList) -> int:
         if self.checkComputedApproxPivoted(): return -1
         RROMPyAssert(self._mode, message = "Cannot setup approximant.")
         vbMng(self, "INIT", "Setting up pivoted approximant.", 10)
         if not hasattr(self, "_plotEstPivot"): self._plotEstPivot = "NONE"
         self.computeScaleFactor()
         if self.trainedModel is None:
             self.trainedModel = self.tModelType()
             self.trainedModel.verbosity = self.verbosity
             self.trainedModel.timestamp = self.timestamp
             datadict = {"mu0": self.mu0, "projMat": np.zeros((0, 0)),
                         "scaleFactor": self.scaleFactor,
                         "rescalingExp": self.HFEngine.rescalingExp,
                         "directionPivot": self.directionPivot}
             self.trainedModel.data = self.initializeModelData(datadict)[0]
             self.trainedModel.data.Qs, self.trainedModel.data.Ps = [], []
         _trainedModelOld = copy(self.trainedModel)
         self._scaleFactorOldPivot = copy(self.scaleFactor)
         self.scaleFactor = self.scaleFactorPivot
         self._temporaryPivot = 1
         self._samplingEngineOld = copy(self.samplingEngine)
         self.musMargLoc, self.samplingEngs = [], [None] * len(mus)
         Qs, Ps = [None] * len(mus), [None] * len(mus)
         self.verbosity -= 15
         S0 = copy(self.S)
         for j, mu in enumerate(mus):
             RationalInterpolantGreedy.setupSampling(self)
             self.trainedModel = None
             self.musMargLoc += [mu]
             RationalInterpolantGreedy.setupApprox(self, self._plotEstPivot)
             self.samplingEngs[j] = copy(self.samplingEngine)
             Qs[j] = copy(self.trainedModel.data.Q)
             Ps[j] = copy(self.trainedModel.data.P)
             self._S = S0
         self.scaleFactor = self._scaleFactorOldPivot
         del self._scaleFactorOldPivot, self._temporaryPivot
         self._finalizeSnapshots()
         del self._samplingEngineOld, self.musMargLoc, self.samplingEngs
         self._mus = self.samplingEngine.musCoalesced
         self.trainedModel = _trainedModelOld
         self.trainedModel.data.mus = copy(self.mus)
         self.trainedModel.data.musMarginal = copy(self.musMarginal)
         padRight = (self.samplingEngine.nsamplesTot
                   - self.trainedModel.data.projMat.shape[1])
         nmusOld = len(self.trainedModel.data.Ps)
         for j in range(nmusOld):
             nsj = self.samplingEngine.nsamples[j]
             self.trainedModel.data.Ps[j].pad(0, padRight)
             self.trainedModel.data.HIs[j].pad(0, padRight)
         padLeft = self.trainedModel.data.projMat.shape[1]
         for j in range(len(mus)):
             nsj = self.samplingEngine.nsamples[nmusOld + j]
             if self.POD == PODGlobal:
                 rRightj = self.samplingEngine.RPODCPart[:,
                                                        padLeft : padLeft + nsj]
                 Ps[j].postmultiplyTensorize(rRightj.T)
             else:
                 padRight -= nsj
                 Ps[j].pad(padLeft, padRight)
             padLeft += nsj
         pMat = self.samplingEngine.samplesCoalesced.data
         pMatEff = dot(self.HFEngine.C, pMat) if self.approx_state else pMat
         self.trainedModel.data.projMat = pMatEff
         self.trainedModel.data.Qs += Qs
         self.trainedModel.data.Ps += Ps
         self.trainedModel.data.approxParameters = copy(self.approxParameters)
         self.verbosity += 15
         vbMng(self, "DEL", "Done setting up pivoted approximant.", 10)
         return 0
 
     def setupApprox(self, plotEst : str = "NONE") -> int:
         if self.checkComputedApprox(): return -1
         if '_' not in plotEst: plotEst = plotEst + "_NONE"
         plotEstM, self._plotEstPivot = plotEst.split("_")
         val = super().setupApprox(plotEstM)
         return val
diff --git a/rrompy/reduction_methods/pivoted/greedy/rational_interpolant_pivoted_greedy.py b/rrompy/reduction_methods/pivoted/greedy/rational_interpolant_pivoted_greedy.py
index bdd7e1a..3a42400 100644
--- a/rrompy/reduction_methods/pivoted/greedy/rational_interpolant_pivoted_greedy.py
+++ b/rrompy/reduction_methods/pivoted/greedy/rational_interpolant_pivoted_greedy.py
@@ -1,321 +1,325 @@
 # Copyright (C) 2018 by the RROMPy authors
 #
 # This file is part of RROMPy.
 #
 # RROMPy is free software: you can redistribute it and/or modify
 # it under the terms of the GNU Lesser General Public License as published by
 # the Free Software Foundation, either version 3 of the License, or
 # (at your option) any later version.
 #
 # RROMPy is distributed in the hope that it will be useful,
 # but WITHOUT ANY WARRANTY; without even the implied warranty of
 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
 # GNU Lesser General Public License for more details.
 #
 # You should have received a copy of the GNU Lesser General Public License
 # along with RROMPy. If not, see <http://www.gnu.org/licenses/>.
 #
 
 from copy import deepcopy as copy
 import numpy as np
 from .generic_pivoted_greedy_approximant import GenericPivotedGreedyApproximant
 from rrompy.utilities.numerical import dot
 from rrompy.reduction_methods.standard import RationalInterpolant
 from rrompy.reduction_methods.pivoted import RationalInterpolantPivoted
 from rrompy.reduction_methods.pivoted.generic_pivoted_approximant import (
                                                                      PODGlobal)
 from rrompy.utilities.base.types import paramList
 from rrompy.utilities.base import verbosityManager as vbMng
 from rrompy.utilities.exception_manager import RROMPyAssert
 from rrompy.parameter import checkParameterList, emptyParameterList
 
 __all__ = ['RationalInterpolantPivotedGreedy']
 
 class RationalInterpolantPivotedGreedy(GenericPivotedGreedyApproximant,
                                        RationalInterpolantPivoted):
     """
     ROM pivoted greedy rational interpolant computation for parametric
         problems.
 
     Args:
         HFEngine: HF problem solver.
         mu0(optional): Default parameter. Defaults to 0.
         directionPivot(optional): Pivot components. Defaults to [0].
         approxParameters(optional): Dictionary containing values for main
             parameters of approximant. Recognized keys are:
             - 'POD': 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;
             - 'cutOffKind': kind of cut off strategy; available values
                 include 'SOFT' and 'HARD'; defaults to 'HARD';
             - 'matchingWeightError': weight for pole matching optimization in
                 error estimation; defaults to 0;
             - 'cutOffToleranceError': tolerance for ignoring parasitic poles
                 in error estimation; defaults to 'AUTO', i.e. cutOffTolerance;
             - 'S': total number of pivot samples current approximant relies
                 upon;
             - 'samplerPivot': pivot sample point generator;
             - 'SMarginal': number of starting marginal samples;
             - 'samplerMarginalGrid': marginal sample point generator via sparse
                 grid;
             - 'polybasis': type of polynomial basis for pivot interpolation;
                 defaults to 'MONOMIAL';
             - 'polybasisMarginal': type of polynomial basis for marginal
                 interpolation; allowed values include 'MONOMIAL', 'CHEBYSHEV'
                 and 'LEGENDRE'; 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;
             - 'MMarginal': degree of marginal interpolant; defaults to 'AUTO',
                 i.e. maximum allowed;
             - 'greedyTolMarginal': uniform error tolerance for marginal greedy
                algorithm; defaults to 1e-1;
             - 'maxIterMarginal': maximum number of marginal greedy steps;
                defaults to 1e2;
             - 'polydegreetypeMarginal': type of polynomial degree for marginal;
                 defaults to 'TOTAL';
             - 'radialDirectionalWeights': radial basis weights for pivot
                 numerator; defaults to 1;
             - 'radialDirectionalWeightsMarginal': radial basis weights for
                 marginal interpolant; defaults to 1;
             - 'nNearestNeighbor': number of pivot nearest neighbors considered
                 if polybasis allows; defaults to -1;
             - 'nNearestNeighborMarginal': number of marginal nearest neighbors
                 considered if polybasisMarginal allows; defaults to -1;
             - 'interpRcond': tolerance for pivot interpolation; defaults to
                 None;
             - 'interpRcondMarginal': tolerance for marginal 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.
         directionPivot: Pivot components.
         mus: Array of snapshot parameters.
         musPivot: Array of pivot snapshot parameters.
         musMarginal: Array of marginal snapshot parameters.
         approxParameters: Dictionary containing values for main parameters of
             approximant. Recognized keys are in parameterList.
         parameterListSoft: Recognized keys of soft approximant parameters:
             - 'POD': whether to compute POD of snapshots;
+            - 'scaleFactorDer': scaling factors for derivative computation;
             - 'matchingWeight': weight for pole matching optimization;
             - 'cutOffTolerance': tolerance for ignoring parasitic poles;
             - 'cutOffKind': kind of cut off strategy;
             - 'matchingWeightError': weight for pole matching optimization in
                 error estimation;
             - 'cutOffToleranceError': tolerance for ignoring parasitic poles
                 in error estimation;
             - 'polybasis': type of polynomial basis for pivot interpolation;
             - 'polybasisMarginal': type of polynomial basis for marginal
                 interpolation;
             - 'M': degree of rational interpolant numerator;
             - 'N': degree of rational interpolant denominator;
             - 'MMarginal': degree of marginal interpolant;
             - 'greedyTolMarginal': uniform error tolerance for marginal greedy
                 algorithm;
             - 'maxIterMarginal': maximum number of marginal greedy steps;
             - 'polydegreetypeMarginal': type of polynomial degree for marginal;
             - 'radialDirectionalWeights': radial basis weights for pivot
                 numerator;
             - 'radialDirectionalWeightsMarginal': radial basis weights for
                 marginal interpolant;
             - 'nNearestNeighbor': number of pivot nearest neighbors considered
                 if polybasis allows;
             - 'nNearestNeighborMarginal': number of marginal nearest neighbors
                 considered if polybasisMarginal allows;
             - 'interpRcond': tolerance for pivot interpolation;
             - 'interpRcondMarginal': tolerance for marginal interpolation;
             - '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 pivot samples current approximant relies
                 upon;
             - 'samplerPivot': pivot sample point generator;
             - 'SMarginal': total number of marginal samples current approximant
                 relies upon;
             - 'samplerMarginalGrid': marginal sample point generator via sparse
                 grid.
         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.
         cutOffKind: Kind of cut off strategy.
         matchingWeightError: Weight for pole matching optimization in error
             estimation.
         cutOffToleranceError: Tolerance for ignoring parasitic poles in error
             estimation.
         S: Total number of pivot samples current approximant relies upon.
         samplerPivot: Pivot sample point generator.
         SMarginal: Total number of marginal samples current approximant relies
             upon.
         samplerMarginalGrid: Marginal sample point generator via sparse grid.
         polybasis: Type of polynomial basis for pivot interpolation.
         polybasisMarginal: Type of polynomial basis for marginal interpolation.
         M: Degree of rational interpolant numerator.
         N: Degree of rational interpolant denominator.
         MMarginal: Degree of marginal interpolant.
         greedyTolMarginal: Uniform error tolerance for marginal greedy
             algorithm.
         maxIterMarginal: Maximum number of marginal greedy steps.
         polydegreetypeMarginal: Type of polynomial degree for marginal.
         radialDirectionalWeights: Radial basis weights for pivot numerator.
         radialDirectionalWeightsMarginal: Radial basis weights for marginal
             interpolant.
         nNearestNeighbor: Number of pivot nearest neighbors considered if
             polybasis allows.
         nNearestNeighborMarginal: Number of marginal nearest neighbors
             considered if polybasisMarginal allows.
         interpRcond: Tolerance for pivot interpolation.
         interpRcondMarginal: Tolerance for marginal interpolation.
         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 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 _finalizeSnapshots(self):
         self.samplingEngine = self._samplingEngineOld
         for muM, sEN in zip(self.musMargLoc, self.samplingEngs):
             self.samplingEngine.samples += [sEN.samples]
             self.samplingEngine.nsamples += [sEN.nsamples]
             self.samplingEngine.mus += [sEN.mus]
             self.samplingEngine.musMarginal.append(muM)
             self.samplingEngine._derIdxs += [[(0,) * self.npar] 
                                                   for _ in range(sEN.nsamples)]
             if self.POD:
                 self.samplingEngine.RPOD += [sEN.RPOD]
                 self.samplingEngine.samples_full += [copy(sEN.samples_full)]
         if self.POD == PODGlobal:
             self.samplingEngine.coalesceSamples(self.interpRcondMarginal)
         else:
             self.samplingEngine.coalesceSamples()
 
     def computeSnapshots(self):
         """Compute snapshots of solution map."""
         RROMPyAssert(self._mode,
                      message = "Cannot start snapshot computation.")
         vbMng(self, "INIT", "Starting computation of snapshots.", 5)
-        self.samplingEngine.scaleFactor = self._scaleFactorOldPivot
+        self.samplingEngine.scaleFactor = self.scaleFactorDer
         self.musPivot = self.samplerPivot.generatePoints(self.S)
         while len(self.musPivot) > self.S: self.musPivot.pop()
         self.mus = emptyParameterList()
         self.mus.reset((self.S, self.HFEngine.npar))
         self.samplingEngine.resetHistory()
         for k in range(self.S):
             self.mus.data[k, self.directionPivot] = self.musPivot[k].data
             self.mus.data[k, self.directionMarginal] = self.musMargLoc[-1].data
         self.samplingEngine.iterSample(self.mus)
         vbMng(self, "DEL", "Done computing snapshots.", 5)
         self._m_selfmus = copy(self.mus)
         self._mus = self.musPivot
         self._m_mu0 = copy(self.mu0)
         self._m_HFErescalingExp = copy(self.HFEngine.rescalingExp)
         self._mu0 = checkParameterList(self.mu0(self.directionPivot), 1)[0]
         self.HFEngine.rescalingExp = [self.HFEngine.rescalingExp[
                                                        self.directionPivot[0]]]
 
     def setupApproxPivoted(self, mus:paramList) -> int:
         if self.checkComputedApproxPivoted(): return -1
         RROMPyAssert(self._mode, message = "Cannot setup approximant.")
         vbMng(self, "INIT", "Setting up pivoted approximant.", 10)
         self.computeScaleFactor()
         if self.trainedModel is None:
             self.trainedModel = self.tModelType()
             self.trainedModel.verbosity = self.verbosity
             self.trainedModel.timestamp = self.timestamp
             datadict = {"mu0": self.mu0, "projMat": np.zeros((0, 0)),
                         "scaleFactor": self.scaleFactor,
                         "rescalingExp": self.HFEngine.rescalingExp,
                         "directionPivot": self.directionPivot}
             self.trainedModel.data = self.initializeModelData(datadict)[0]
             self.trainedModel.data.Qs, self.trainedModel.data.Ps = [], []
         _trainedModelOld = copy(self.trainedModel)
         self._scaleFactorOldPivot = copy(self.scaleFactor)
         self.scaleFactor = self.scaleFactorPivot
         self._temporaryPivot = 1
         self._samplingEngineOld = copy(self.samplingEngine)
         self.musMargLoc, self.samplingEngs = [], [None] * len(mus)
         Qs, Ps = [None] * len(mus), [None] * len(mus)
         self.verbosity -= 15
         for j, mu in enumerate(mus):
             RationalInterpolant.setupSampling(self)
             self.trainedModel = None
             self.musMargLoc += [mu]
             RationalInterpolant.setupApprox(self)
             self._mu0 = self._m_mu0
             self._mus = self._m_selfmus
             self.HFEngine.rescalingExp = self._m_HFErescalingExp
             del self._m_mu0, self._m_selfmus, self._m_HFErescalingExp
             self.samplingEngs[j] = copy(self.samplingEngine)
             Qs[j] = copy(self.trainedModel.data.Q)
             Ps[j] = copy(self.trainedModel.data.P)
         self.scaleFactor = self._scaleFactorOldPivot
         del self._scaleFactorOldPivot, self._temporaryPivot
         self._finalizeSnapshots()
         del self._samplingEngineOld, self.musMargLoc, self.samplingEngs
         self._mus = self.samplingEngine.musCoalesced
         self.trainedModel = _trainedModelOld
         self.trainedModel.data.mus = copy(self.mus)
         self.trainedModel.data.musMarginal = copy(self.musMarginal)
         padRight = (self.samplingEngine.nsamplesTot
                   - self.trainedModel.data.projMat.shape[1])
         nmusOld = len(self.trainedModel.data.Ps)
         for j in range(nmusOld):
             nsj = self.samplingEngine.nsamples[j]
             self.trainedModel.data.Ps[j].pad(0, padRight)
             self.trainedModel.data.HIs[j].pad(0, padRight)
         padLeft = self.trainedModel.data.projMat.shape[1]
         for j in range(len(mus)):
             nsj = self.samplingEngine.nsamples[nmusOld + j]
             if self.POD == PODGlobal:
                 rRightj = self.samplingEngine.RPODCPart[:,
                                                        padLeft : padLeft + nsj]
                 Ps[j].postmultiplyTensorize(rRightj.T)
             else:
                 padRight -= nsj
                 Ps[j].pad(padLeft, padRight)
             padLeft += nsj
         pMat = self.samplingEngine.samplesCoalesced.data
         pMatEff = dot(self.HFEngine.C, pMat) if self.approx_state else pMat
         self.trainedModel.data.projMat = pMatEff
         self.trainedModel.data.Qs += Qs
         self.trainedModel.data.Ps += Ps
         self.verbosity += 15
         vbMng(self, "DEL", "Done setting up pivoted approximant.", 10)
         return 0
diff --git a/rrompy/reduction_methods/pivoted/rational_interpolant_greedy_pivoted.py b/rrompy/reduction_methods/pivoted/rational_interpolant_greedy_pivoted.py
index a1e293b..cd73f00 100644
--- a/rrompy/reduction_methods/pivoted/rational_interpolant_greedy_pivoted.py
+++ b/rrompy/reduction_methods/pivoted/rational_interpolant_greedy_pivoted.py
@@ -1,438 +1,442 @@
 # Copyright (C) 2018 by the RROMPy authors
 #
 # This file is part of RROMPy.
 #
 # RROMPy is free software: you can redistribute it and/or modify
 # it under the terms of the GNU Lesser General Public License as published by
 # the Free Software Foundation, either version 3 of the License, or
 # (at your option) any later version.
 #
 # RROMPy is distributed in the hope that it will be useful,
 # but WITHOUT ANY WARRANTY; without even the implied warranty of
 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
 # GNU Lesser General Public License for more details.
 #
 # You should have received a copy of the GNU Lesser General Public License
 # along with RROMPy. If not, see <http://www.gnu.org/licenses/>.
 #
 
 from copy import deepcopy as copy
 import numpy as np
 from .generic_pivoted_approximant import GenericPivotedApproximant, PODGlobal
 from rrompy.reduction_methods.standard.greedy.rational_interpolant_greedy \
                                                import RationalInterpolantGreedy
 from rrompy.reduction_methods.standard.greedy.generic_greedy_approximant \
                                                             import pruneSamples
 from rrompy.utilities.base.types import Np1D
 from rrompy.utilities.base import verbosityManager as vbMng
 from rrompy.utilities.numerical import dot
 from rrompy.utilities.numerical.degree import totalDegreeN
 from rrompy.utilities.poly_fitting.polynomial import polyvander as pv
 from rrompy.utilities.exception_manager import RROMPyAssert
 from rrompy.parameter import emptyParameterList, checkParameterList
 
 __all__ = ['RationalInterpolantGreedyPivoted']
 
 class RationalInterpolantGreedyPivoted(GenericPivotedApproximant,
                                        RationalInterpolantGreedy):
     """
     ROM pivoted rational interpolant (with pole matching) computation for
         parametric problems.
 
     Args:
         HFEngine: HF problem solver.
         mu0(optional): Default parameter. Defaults to 0.
         directionPivot(optional): Pivot components. Defaults to [0].
         approxParameters(optional): Dictionary containing values for main
             parameters of approximant. Recognized keys are:
             - 'POD': 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;
             - 'cutOffKind': kind of cut off strategy; available values
                 include 'SOFT' and 'HARD'; defaults to 'HARD';
             - 'S': total number of pivot samples current approximant relies
                 upon;
             - 'samplerPivot': pivot sample point generator;
             - 'SMarginal': total number of marginal samples current approximant
                 relies upon;
             - 'samplerMarginal': marginal sample point generator;
             - 'polybasis': type of polynomial basis for pivot
                 interpolation; defaults to 'MONOMIAL';
             - 'polybasisMarginal': type of polynomial basis for marginal
                 interpolation; allowed values include 'MONOMIAL', 'CHEBYSHEV'
                 and 'LEGENDRE'; defaults to 'MONOMIAL';
             - 'greedyTol': uniform error tolerance for greedy algorithm;
                 defaults to 1e-2;
             - 'collinearityTol': collinearity tolerance for greedy algorithm;
                 defaults to 0.;
             - 'maxIter': maximum number of greedy steps; defaults to 1e2;
             - 'nTestPoints': number of test points; defaults to 5e2;
             - 'trainSetGenerator': training sample points generator; defaults
                 to sampler;
             - 'errorEstimatorKind': kind of error estimator; available values
                 include 'AFFINE', 'DISCREPANCY', 'LOOK_AHEAD',
                 'LOOK_AHEAD_RES', and 'NONE'; defaults to 'NONE';
             - '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;
             - 'nNearestNeighborMarginal': number of marginal nearest neighbors
                 considered if polybasisMarginal allows; defaults to -1;
             - 'interpRcond': tolerance for pivot interpolation; defaults to
                 None;
             - 'interpRcondMarginal': tolerance for marginal 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.
         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;
             - 'cutOffKind': kind of cut off strategy;
             - 'polybasis': type of polynomial basis for pivot
                 interpolation;
             - 'polybasisMarginal': type of polynomial basis for marginal
                 interpolation;
             - 'greedyTol': uniform error tolerance for greedy algorithm;
             - 'collinearityTol': collinearity tolerance for greedy algorithm;
             - 'maxIter': maximum number of greedy steps;
             - 'nTestPoints': number of test points;
             - 'trainSetGenerator': training sample points generator;
             - 'errorEstimatorKind': kind of error estimator;
             - 'MMarginal': degree of marginal interpolant;
             - 'polydegreetypeMarginal': type of polynomial degree for marginal;
             - 'radialDirectionalWeightsMarginal': radial basis weights for
                 marginal interpolant;
             - 'nNearestNeighbor': number of pivot nearest neighbors considered
                 if polybasis allows;
             - 'nNearestNeighborMarginal': number of marginal nearest neighbors
                 considered if polybasisMarginal allows;
             - 'interpRcond': tolerance for pivot interpolation;
             - 'interpRcondMarginal': tolerance for marginal interpolation;
             - '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 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.
         cutOffKind: Kind of 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.
         polybasis: Type of polynomial basis for pivot interpolation.
         polybasisMarginal: Type of polynomial basis for marginal interpolation.
         greedyTol: uniform error tolerance for greedy algorithm.
         collinearityTol: Collinearity tolerance for greedy algorithm.
         maxIter: maximum number of greedy steps.
         nTestPoints: number of starting training points.
         trainSetGenerator: training sample points generator.
         errorEstimatorKind: kind of error estimator.
         MMarginal: Degree of marginal interpolant.
         polydegreetypeMarginal: Type of polynomial degree for marginal.
         radialDirectionalWeightsMarginal: Radial basis weights for marginal
             interpolant.
         nNearestNeighborMarginal: Number of marginal nearest neighbors
             considered if polybasisMarginal allows.
         interpRcond: Tolerance for pivot interpolation.
         interpRcondMarginal: Tolerance for marginal interpolation.
         robustTol: Tolerance for robust rational denominator management.
         muBounds: list of bounds for pivot parameter values.
         muBoundsMarginal: list of bounds for marginal parameter values.
         samplingEngine: Sampling engine.
         uHF: High fidelity solution(s) with parameter(s) lastSolvedHF as
             sampleList.
         lastSolvedHF: Parameter(s) corresponding to last computed high fidelity
             solution(s) as parameterList.
         uApproxReduced: Reduced approximate solution(s) with parameter(s)
             lastSolvedApprox as sampleList.
         lastSolvedApproxReduced: Parameter(s) corresponding to last computed
             reduced approximate solution(s) as parameterList.
         uApprox: Approximate solution(s) with parameter(s) lastSolvedApprox as
             sampleList.
         lastSolvedApprox: Parameter(s) corresponding to last computed
             approximate solution(s) as parameterList.
         Q: Numpy 1D vector containing complex coefficients of approximant
             denominator.
         P: Numpy 2D vector whose columns are FE dofs of coefficients of
             approximant numerator.
     """
     
     def __init__(self, *args, **kwargs):
         self._preInit()
         self._addParametersToList(toBeExcluded = ["sampler"])
         super().__init__(*args, **kwargs)
         self._postInit()
 
     @property
     def tModelType(self):
         if hasattr(self, "_temporaryPivot"):
             return RationalInterpolantGreedy.tModelType.fget(self)
         return super().tModelType
 
     @property
     def polybasis0(self):
         if "_" in self.polybasis:
             return self.polybasis.split("_")[0]
         return self.polybasis
 
     def _polyvanderAuxiliary(self, mus, deg, *args):
         degEff = [0] * self.npar
         degEff[self.directionPivot[0]] = deg
         return pv(mus, degEff, *args)
 
     def _marginalizeMiscellanea(self, forward:bool):
         if forward:
             self._m_mu0 = copy(self.mu0)
             self._m_selfmus = copy(self.mus)
             self._m_HFErescalingExp = copy(self.HFEngine.rescalingExp)
             self._mu0 = checkParameterList(self.mu0(self.directionPivot), 1)[0]
             self._mus = checkParameterList(self.mus(self.directionPivot), 1)[0]
             self.HFEngine.rescalingExp = [self.HFEngine.rescalingExp[
                                                        self.directionPivot[0]]]
         else:
             self._mu0 = self._m_mu0
             self._mus = self._m_selfmus
             self.HFEngine.rescalingExp = self._m_HFErescalingExp
             del self._m_mu0, self._m_selfmus, self._m_HFErescalingExp
 
     def _marginalizeTrainedModel(self, forward:bool):
         if forward:
             del self._temporaryPivot
             self.trainedModel.data.mu0 = self.mu0
             self.trainedModel.data.scaleFactor = [1.] * self.npar
             self.trainedModel.data.scaleFactor[self.directionPivot[0]] = (
                                                            self.scaleFactor[0])
             self.trainedModel.data.rescalingExp = self.HFEngine.rescalingExp
             Qc = np.zeros((len(self.trainedModel.data.Q.coeffs),) * self.npar,
                           dtype = self.trainedModel.data.Q.coeffs.dtype)
             Pc = np.zeros((len(self.trainedModel.data.P.coeffs),) * self.npar
                         + (self.trainedModel.data.P.coeffs.shape[1],),
                           dtype = self.trainedModel.data.P.coeffs.dtype)
             for j in range(len(self.trainedModel.data.Q.coeffs)):
                 Qc[(0,) * self.directionPivot[0] + (j,)
                  + (0,) * (self.npar - self.directionPivot[0] - 1)] = (
                                             self.trainedModel.data.Q.coeffs[j])
             for j in range(len(self.trainedModel.data.P.coeffs)):
                 for k in range(self.trainedModel.data.P.coeffs.shape[1]):
                     Pc[(0,) * self.directionPivot[0] + (j,)
                      + (0,) * (self.npar - self.directionPivot[0] - 1)
                      + (k,)] = self.trainedModel.data.P.coeffs[j, k]
             self.trainedModel.data.Q.coeffs = Qc
             self.trainedModel.data.P.coeffs = Pc
 
             self._m_musUniqueCN = copy(self._musUniqueCN)
             musUniqueCNAux = np.zeros((self.S, self.npar),
                                       dtype = self._musUniqueCN.dtype)
             musUniqueCNAux[:, self.directionPivot[0]] = self._musUniqueCN(0)
             self._musUniqueCN = checkParameterList(musUniqueCNAux,
                                                    self.npar)[0]
             self._m_derIdxs = copy(self._derIdxs)
             for j in range(len(self._derIdxs)):
                 for l in range(len(self._derIdxs[j])):
                     derjl = self._derIdxs[j][l][0]
                     self._derIdxs[j][l] = [0] * self.npar
                     self._derIdxs[j][l][self.directionPivot[0]] = derjl
         else:
             self._temporaryPivot = 1
             self.trainedModel.data.mu0 = checkParameterList(
                                            self.mu0(self.directionPivot), 1)[0]
             self.trainedModel.data.scaleFactor = self.scaleFactor
             self.trainedModel.data.rescalingExp = self.HFEngine.rescalingExp[
                                                         self.directionPivot[0]]
             self.trainedModel.data.Q.coeffs = self.trainedModel.data.Q.coeffs[
                                                (0,) * self.directionPivot[0]
                                              + (slice(None),)
                                              + (0,) * (self.HFEngine.npar - 1
                                                      - self.directionPivot[0])]
             self.trainedModel.data.P.coeffs = self.trainedModel.data.P.coeffs[
                                                (0,) * self.directionPivot[0]
                                              + (slice(None),)
                                              + (0,) * (self.HFEngine.npar - 1
                                                      - self.directionPivot[0])]
 
             self._musUniqueCN = copy(self._m_musUniqueCN)
             self._derIdxs = copy(self._m_derIdxs)
             del self._m_musUniqueCN, self._m_derIdxs
         self.trainedModel.data.npar = self.npar
         self.trainedModel.data.Q.npar = self.npar
         self.trainedModel.data.P.npar = self.npar
 
     def errorEstimator(self, mus:Np1D, return_max : bool = False) -> Np1D:
         """Standard residual-based error estimator."""
         self._marginalizeMiscellanea(True)
         setupOK = self.setupApproxLocal()
         self._marginalizeMiscellanea(False)
         if setupOK > 0:
             err = np.empty(len(mus))
             err[:] = np.nan
             if not return_max: return err
             return err, [- setupOK], np.nan
         self._marginalizeTrainedModel(True)
         errRes = super().errorEstimator(mus, return_max)
         self._marginalizeTrainedModel(False)
         return errRes
 
     def _preliminaryTraining(self):
         """Initialize starting snapshots of solution map."""
         RROMPyAssert(self._mode, message = "Cannot start greedy algorithm.")
         S = self.S
         self.sampleBatchIdx, self.sampleBatchSize, self._S = -1, 0, 0
         nextBatchSize = 1
         while self._S + nextBatchSize <= S:
             self.sampleBatchIdx += 1
             self.sampleBatchSize = nextBatchSize
             self._S += self.sampleBatchSize
             nextBatchSize = totalDegreeN(self.npar - 1,
                                          self.sampleBatchIdx + 1)
         self.resetSamples()
-        self.samplingEngine.scaleFactor = self.scaleFactor
+        self.samplingEngine.scaleFactor = self.scaleFactorDer
         musPivot = self.trainSetGenerator.generatePoints(self.S)
         while len(musPivot) > self.S: musPivot.pop()
         muTestPivot = self.samplerPivot.generatePoints(self.nTestPoints, False)
         idxPop = pruneSamples(muTestPivot ** self.HFEngine.rescalingExp[
                                                        self.directionPivot[0]],
                               musPivot ** self.HFEngine.rescalingExp[
                                                        self.directionPivot[0]],
                               1e-10 * self.scaleFactor[0])
         self.mus = emptyParameterList()
         self.mus.reset((self.S, self.npar + len(self.musMargLoc)))
         muTestBase = emptyParameterList()
         muTestBase.reset((len(muTestPivot), self.npar + len(self.musMargLoc)))
         for k in range(self.S):
             self.mus.data[k, self.directionPivot] = musPivot[k].data
             self.mus.data[k, self.directionMarginal] = self.musMargLoc.data
         for k in range(len(muTestPivot)):
             muTestBase.data[k, self.directionPivot] = muTestPivot[k].data
             muTestBase.data[k, self.directionMarginal] = self.musMargLoc.data
         muTestBase.pop(idxPop)
         muLast = copy(self.mus[-1])
         self.mus.pop()
         if len(self.mus) > 0:
             vbMng(self, "MAIN", 
                   ("Adding first {} sample point{} at {} to training "
                    "set.").format(self.S - 1, "" + "s" * (self.S > 2),
                                   self.mus), 3)
             self.samplingEngine.iterSample(self.mus)
         self._S = len(self.mus)
         self._approxParameters["S"] = self.S
         self.muTest = emptyParameterList()
         self.muTest.reset((len(muTestBase) + 1, self.mus.shape[1]))
         self.muTest.data[: -1] = muTestBase.data
         self.muTest.data[-1] = muLast.data
         self.M, self.N = ("AUTO",) * 2
 
     def _finalizeSnapshots(self):
         self.setupSampling()
         self.samplingEngine.resetHistory(len(self.musMarginal))
         for j in range(len(self.musMarginal)):
             self.samplingEngine.setsample(self.samplingEngs[j].samples,
                                           j, False)
             self.samplingEngine.mus[j] = copy(self.samplingEngs[j].mus)
             self.samplingEngine.musMarginal[j] = copy(self.musMarginal[j])
             self.samplingEngine.nsamples[j] = self.samplingEngs[j].nsamples
             if self.POD:
                 self.samplingEngine.RPOD[j] = self.samplingEngs[j].RPOD
                 self.samplingEngine.samples_full[j].data = (
                                         self.samplingEngs[j].samples_full.data)
         if self.POD == PODGlobal:
             self.samplingEngine.coalesceSamples(self.interpRcondMarginal)
         else:
             self.samplingEngine.coalesceSamples()
 
     def setupApprox(self, *args, **kwargs) -> int:
         """Compute rational interpolant."""
         if self.checkComputedApprox(): return -1
         RROMPyAssert(self._mode, message = "Cannot setup approximant.")
         vbMng(self, "INIT", "Setting up {}.". format(self.name()), 5)
         self.musMarginal = self.samplerMarginal.generatePoints(self.SMarginal)
         while len(self.musMarginal) > self.SMarginal: self.musMarginal.pop()
         S0 = copy(self.S)
         Qs, Ps = [None] * len(self.musMarginal), [None] * len(self.musMarginal)
         self.samplingEngs = [None] * len(self.musMarginal)
         self.computeScaleFactor()
         self._scaleFactorOldPivot = copy(self.scaleFactor)
         self.scaleFactor = self.scaleFactorPivot
         self._temporaryPivot = 1
         for j in range(len(self.musMarginal)):
             self._S = S0
             self.musMargLoc = self.musMarginal[j]
             RationalInterpolantGreedy.setupSampling(self)
             self.trainedModel = None
             self.verbosity -= 5
             self.samplingEngine.verbosity -= 5
             super().setupApprox(*args, **kwargs)
             self.verbosity += 5
             self.samplingEngine.verbosity += 5
             self.samplingEngs[j] = copy(self.samplingEngine)
             Qs[j] = copy(self.trainedModel.data.Q)
             Ps[j] = copy(self.trainedModel.data.P)
         self.scaleFactor = self._scaleFactorOldPivot
         del self._scaleFactorOldPivot, self._temporaryPivot
         self._finalizeSnapshots()
         del self.musMargLoc, self.samplingEngs
         self._mus = self.samplingEngine.musCoalesced
         padLeft = 0
         if self.POD != PODGlobal: padRight = self.samplingEngine.nsamplesTot
         for j in range(len(self.musMarginal)):
             nsj = self.samplingEngine.nsamples[j]
             if self.POD == PODGlobal:
                 rRightj = self.samplingEngine.RPODCPart[:,
                                                        padLeft : padLeft + nsj]
                 Ps[j].postmultiplyTensorize(rRightj.T)
             else:
                 padRight -= nsj
                 Ps[j].pad(padLeft, padRight)
             padLeft += nsj
         pMat = self.samplingEngine.samplesCoalesced.data
         pMatEff = dot(self.HFEngine.C, pMat) if self.approx_state else pMat
         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,
                     "directionPivot": self.directionPivot}
         self.trainedModel.data = self.initializeModelData(datadict)[0]
         self.trainedModel.data.mus = copy(self.mus)
         self.trainedModel.data.musMarginal = copy(self.musMarginal)
         self.trainedModel.data.Qs, self.trainedModel.data.Ps = Qs, Ps
         vbMng(self, "INIT", "Matching poles.", 10)
         self.trainedModel.initializeFromRational(
                                       self.HFEngine, self.matchingWeight,
                                       self.POD == PODGlobal, self.approx_state)
         vbMng(self, "DEL", "Done matching poles.", 10)
         self._finalizeMarginalization()
         vbMng(self, "DEL", "Done setting up approximant.", 5)
         return 0
diff --git a/rrompy/reduction_methods/pivoted/rational_interpolant_pivoted.py b/rrompy/reduction_methods/pivoted/rational_interpolant_pivoted.py
index 6352391..da09218 100644
--- a/rrompy/reduction_methods/pivoted/rational_interpolant_pivoted.py
+++ b/rrompy/reduction_methods/pivoted/rational_interpolant_pivoted.py
@@ -1,363 +1,386 @@
 # Copyright (C) 2018 by the RROMPy authors
 #
 # This file is part of RROMPy.
 #
 # RROMPy is free software: you can redistribute it and/or modify
 # it under the terms of the GNU Lesser General Public License as published by
 # the Free Software Foundation, either version 3 of the License, or
 # (at your option) any later version.
 #
 # RROMPy is distributed in the hope that it will be useful,
 # but WITHOUT ANY WARRANTY; without even the implied warranty of
 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
 # GNU Lesser General Public License for more details.
 #
 # You should have received a copy of the GNU Lesser General Public License
 # along with RROMPy. If not, see <http://www.gnu.org/licenses/>.
 #
 
 from copy import deepcopy as copy
 import numpy as np
 from .generic_pivoted_approximant import GenericPivotedApproximant, PODGlobal
 from rrompy.reduction_methods.standard.rational_interpolant import (
                                                            RationalInterpolant)
 from rrompy.utilities.base import verbosityManager as vbMng
 from rrompy.utilities.numerical import dot
 from rrompy.utilities.numerical.hash_derivative import nextDerivativeIndices
 from rrompy.utilities.exception_manager import RROMPyAssert, RROMPyWarning
 from rrompy.parameter import emptyParameterList
 
 __all__ = ['RationalInterpolantPivoted']
 
 class RationalInterpolantPivoted(GenericPivotedApproximant,
                                  RationalInterpolant):
     """
     ROM pivoted rational interpolant (with pole matching) computation for
         parametric problems.
 
     Args:
         HFEngine: HF problem solver.
         mu0(optional): Default parameter. Defaults to 0.
         directionPivot(optional): Pivot components. Defaults to [0].
         approxParameters(optional): Dictionary containing values for main
             parameters of approximant. Recognized keys are:
             - 'POD': 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;
             - 'cutOffKind': kind of cut off strategy; available values
                 include 'SOFT' and 'HARD'; defaults to 'HARD';
             - 'S': total number of pivot samples current approximant relies
                 upon;
             - 'samplerPivot': pivot sample point generator;
             - 'SMarginal': total number of marginal samples current approximant
                 relies upon;
             - 'samplerMarginal': marginal sample point generator;
             - 'polybasis': type of polynomial basis for pivot
                 interpolation; defaults to 'MONOMIAL';
             - 'polybasisMarginal': type of polynomial basis for marginal
                 interpolation; allowed values include 'MONOMIAL', 'CHEBYSHEV'
                 and 'LEGENDRE'; 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;
             - 'MMarginal': degree of marginal interpolant; defaults to 'AUTO',
                 i.e. maximum allowed;
             - 'polydegreetypeMarginal': type of polynomial degree for marginal;
                 defaults to 'TOTAL';
             - 'radialDirectionalWeights': radial basis weights for pivot
                 numerator; defaults to 1;
             - 'radialDirectionalWeightsMarginal': radial basis weights for
                 marginal interpolant; defaults to 1;
             - 'nNearestNeighbor': number of pivot nearest neighbors considered
                 if polybasis allows; defaults to -1;
             - 'nNearestNeighborMarginal': number of marginal nearest neighbors
                 considered if polybasisMarginal allows; defaults to -1;
             - 'interpRcond': tolerance for pivot interpolation; defaults to
                 None;
             - 'interpRcondMarginal': tolerance for marginal 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.
         directionPivot: Pivot components.
         mus: Array of snapshot parameters.
         musPivot: Array of pivot snapshot parameters.
         musMarginal: Array of marginal snapshot parameters.
         approxParameters: Dictionary containing values for main parameters of
             approximant. Recognized keys are in parameterList.
         parameterListSoft: Recognized keys of soft approximant parameters:
             - 'POD': whether to compute POD of snapshots;
+            - 'scaleFactorDer': scaling factors for derivative computation;
             - 'matchingWeight': weight for pole matching optimization;
             - 'cutOffTolerance': tolerance for ignoring parasitic poles;
             - 'cutOffKind': kind of cut off strategy;
             - 'polybasis': type of polynomial basis for pivot
                 interpolation;
             - 'polybasisMarginal': type of polynomial basis for marginal
                 interpolation;
             - 'M': degree of rational interpolant numerator;
             - 'N': degree of rational interpolant denominator;
             - 'MMarginal': degree of marginal interpolant;
             - 'polydegreetypeMarginal': type of polynomial degree for marginal;
             - 'radialDirectionalWeights': radial basis weights for pivot
                 numerator;
             - 'radialDirectionalWeightsMarginal': radial basis weights for
                 marginal interpolant;
             - 'nNearestNeighbor': number of pivot nearest neighbors considered
                 if polybasis allows;
             - 'nNearestNeighborMarginal': number of marginal nearest neighbors
                 considered if polybasisMarginal allows;
             - 'interpRcond': tolerance for pivot interpolation;
             - 'interpRcondMarginal': tolerance for marginal interpolation;
             - '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 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.
         cutOffKind: Kind of 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.
         polybasis: Type of polynomial basis for pivot interpolation.
         polybasisMarginal: Type of polynomial basis for marginal interpolation.
         M: Numerator degree of approximant.
         N: Denominator degree of approximant.
         MMarginal: Degree of marginal interpolant.
         polydegreetypeMarginal: Type of polynomial degree for marginal.
         radialDirectionalWeights: Radial basis weights for pivot numerator.
         radialDirectionalWeightsMarginal: Radial basis weights for marginal
             interpolant.
         nNearestNeighbor: Number of pivot nearest neighbors considered if
             polybasis allows.
         nNearestNeighborMarginal: Number of marginal nearest neighbors
             considered if polybasisMarginal allows.
         interpRcond: Tolerance for pivot interpolation.
         interpRcondMarginal: Tolerance for marginal interpolation.
         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 pivot parameter values.
         muBoundsMarginal: list of bounds for marginal parameter values.
         samplingEngine: Sampling engine.
         uHF: High fidelity solution(s) with parameter(s) lastSolvedHF as
             sampleList.
         lastSolvedHF: Parameter(s) corresponding to last computed high fidelity
             solution(s) as parameterList.
         uApproxReduced: Reduced approximate solution(s) with parameter(s)
             lastSolvedApprox as sampleList.
         lastSolvedApproxReduced: Parameter(s) corresponding to last computed
             reduced approximate solution(s) as parameterList.
         uApprox: Approximate solution(s) with parameter(s) lastSolvedApprox as
             sampleList.
         lastSolvedApprox: Parameter(s) corresponding to last computed
             approximate solution(s) as parameterList.
         Q: Numpy 1D vector containing complex coefficients of approximant
             denominator.
         P: Numpy 2D vector whose columns are FE dofs of coefficients of
             approximant numerator.
     """
     
     def __init__(self, *args, **kwargs):
         self._preInit()
         self._addParametersToList(toBeExcluded = ["polydegreetype", "sampler"])
         super().__init__(*args, **kwargs)
         self._postInit()
 
+    @property
+    def scaleFactorDer(self):
+        """Value of scaleFactorDer."""
+        if self._scaleFactorDer == "NONE": return 1.
+        if self._scaleFactorDer == "AUTO": return self.scaleFactorPivot
+        return self._scaleFactorDer
+    @scaleFactorDer.setter
+    def scaleFactorDer(self, scaleFactorDer):
+        if hasattr(self, "_scaleFactorDer"):
+            scaleFactorDerold = self.scaleFactorDer
+        else: scaleFactorDerold = -1
+        if isinstance(scaleFactorDer, (str,)):
+            scaleFactorDer = scaleFactorDer.upper()
+        self._scaleFactorDer = scaleFactorDer
+        self._approxParameters["scaleFactorDer"] = self._scaleFactorDer
+        if scaleFactorDerold != self._scaleFactorDer:
+            self.samplingEngine = None
+            self.resetSamples()
+
     @property
     def polydegreetype(self):
         """Value of polydegreetype."""
         return "TOTAL"
     @polydegreetype.setter
     def polydegreetype(self, polydegreetype):
         RROMPyWarning(("polydegreetype is used just to simplify inheritance, "
                        "and its value cannot be changed from 'TOTAL'."))
         
     @property
     def polybasis0(self):
         if "_" in self.polybasis:
             return self.polybasis.split("_")[0]
         return self.polybasis
 
     @property
     def correctorTol(self):
         """Value of correctorTol."""
         return self._correctorTol
     @correctorTol.setter
     def correctorTol(self, correctorTol):
         if correctorTol < 0. or (correctorTol > 0. and self.nparPivot > 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.nparPivot > 1):
             RROMPyWarning(("Overriding prescribed max number of corrector "
                            "iterations to 1."))
             correctorMaxIter = 1
         self._correctorMaxIter = correctorMaxIter
         self._approxParameters["correctorMaxIter"] = self.correctorMaxIter
         
     def _setupInterpolationIndices(self):
         """Setup parameters for polyvander."""
         RROMPyAssert(self._mode,
                      message = "Cannot setup interpolation indices.")
         if (self._musUniqueCN is None 
          or len(self._reorder) != len(self.musPivot)):
             try:
                 muPC = self.trainedModel.centerNormalizePivot(self.musPivot)
             except:
                 muPC = self.trainedModel.centerNormalize(self.musPivot)
             self._musUniqueCN, musIdxsTo, musIdxs, musCount = (muPC.unique(
                                     return_index = True, return_inverse = True,
                                     return_counts = True))
             self._musUnique = self.musPivot[musIdxsTo]
             self._derIdxs = [None] * len(self._musUniqueCN)
             self._reorder = np.empty(len(musIdxs), dtype = int)
             filled = 0
             for j, cnt in enumerate(musCount):
                 self._derIdxs[j] = nextDerivativeIndices([], self.nparPivot,
                                                          cnt)
                 jIdx = np.nonzero(musIdxs == j)[0]
                 self._reorder[jIdx] = np.arange(filled, filled + cnt)
                 filled += cnt
 
     def computeSnapshots(self):
         """Compute snapshots of solution map."""
         RROMPyAssert(self._mode,
                      message = "Cannot start snapshot computation.")
         self.computeScaleFactor()
         if self.samplingEngine.nsamplesTot != self.S * self.SMarginal:
             self.resetSamples()
-            self.samplingEngine.scaleFactor = self.scaleFactorPivot
+            self.samplingEngine.scaleFactor = self.scaleFactorDer
             vbMng(self, "INIT", "Starting computation of snapshots.", 5)
             self.musPivot = self.samplerPivot.generatePoints(self.S)
             while len(self.musPivot) > self.S: self.musPivot.pop()
             self.musMarginal = self.samplerMarginal.generatePoints(
                                                                 self.SMarginal)
             while len(self.musMarginal) > self.SMarginal:
                 self.musMarginal.pop()
             self.mus = emptyParameterList()
             self.mus.reset((self.S * self.SMarginal, self.HFEngine.npar))
             self.samplingEngine.resetHistory(self.SMarginal)
             for j, muMarg in enumerate(self.musMarginal):
                 for k in range(j * self.S, (j + 1) * self.S):
                     self.mus.data[k, self.directionPivot] = (
                                             self.musPivot[k - j * self.S].data)
                     self.mus.data[k, self.directionMarginal] = muMarg.data
             self.samplingEngine.iterSample(self.musPivot, self.musMarginal)
             self._finalizeSnapshots()
             vbMng(self, "DEL", "Done computing snapshots.", 5)
 
     def _finalizeSnapshots(self):
         if self.POD == PODGlobal:
             self.samplingEngine.coalesceSamples(self.interpRcondMarginal)
         else:
             self.samplingEngine.coalesceSamples()
 
     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.samplesCoalesced.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,
                         "directionPivot": self.directionPivot}
             self.trainedModel.data = self.initializeModelData(datadict)[0]
         else:
             self.trainedModel = self.trainedModel
             self.trainedModel.data.projMat = copy(pMatEff)
         N0 = copy(self.N)
         Qs, Ps = [None] * len(self.musMarginal), [None] * len(self.musMarginal)
         self._temporaryPivot = 1
         padLeft = 0
         if self.POD:
             self._RPODOldPivot = copy(self.samplingEngine.RPODCoalesced)
         else:
             self._samplesOldPivot = copy(self.samplingEngine.samples)
             padRight = self.samplingEngine.nsamplesTot
         self._scaleFactorOldPivot = copy(self.scaleFactor)
         self.scaleFactor = self.scaleFactorPivot
         for j in range(len(self.musMarginal)):
             self.N = N0
             if self.POD:
                 self.samplingEngine.RPOD = (
                              self._RPODOldPivot[:, padLeft : padLeft + self.S])
             else:
                 self.samplingEngine.samples = self._samplesOldPivot[j]
                 padRight -= self.S
             self.verbosity -= 5
             self._iterCorrector()
             self.verbosity += 5
             Qs[j] = copy(self.trainedModel.data.Q)
             Ps[j] = copy(self.trainedModel.data.P)
             del self.trainedModel.data.Q, self.trainedModel.data.P
             if not self.POD: Ps[j].pad(padLeft, padRight)
             padLeft += self.S
         if self.POD:
             self.samplingEngine.RPODCoalesced = copy(self._RPODOldPivot)
             del self._RPODOldPivot
         else:
             self.samplingEngine.samples = copy(self._samplesOldPivot)
             del self._samplesOldPivot
         self.scaleFactor = self._scaleFactorOldPivot
         del self._temporaryPivot, self._scaleFactorOldPivot
         self.trainedModel.data.mus = copy(self.mus)
         self.trainedModel.data.musPivot = copy(self.musPivot)
         self.trainedModel.data.musMarginal = copy(self.musMarginal)
         self.trainedModel.data.Qs, self.trainedModel.data.Ps = Qs, Ps
         vbMng(self, "INIT", "Matching poles.", 10)
         self.trainedModel.initializeFromRational(
                                       self.HFEngine, self.matchingWeight,
                                       self.POD == PODGlobal, self.approx_state)
         vbMng(self, "DEL", "Done matching poles.", 10)
         self._finalizeMarginalization()
         vbMng(self, "DEL", "Done setting up approximant.", 5)
         return 0
diff --git a/rrompy/reduction_methods/standard/generic_standard_approximant.py b/rrompy/reduction_methods/standard/generic_standard_approximant.py
index b52d38e..09fd64a 100644
--- a/rrompy/reduction_methods/standard/generic_standard_approximant.py
+++ b/rrompy/reduction_methods/standard/generic_standard_approximant.py
@@ -1,144 +1,148 @@
 # 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.utilities.base import verbosityManager as vbMng
 from rrompy.utilities.exception_manager import RROMPyException, RROMPyAssert
 from rrompy.parameter import checkParameterList
 
 __all__ = ['GenericStandardApproximant']
 
 class GenericStandardApproximant(GenericApproximant):
     """
     ROM interpolant computation for parametric problems (ABSTRACT).
 
     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.
             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.
+            - 'POD': whether to compute POD of snapshots;
+            - 'scaleFactorDer': scaling factors for derivative computation.
         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.
         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.
     """
 
     def __init__(self, *args, **kwargs):
         self._preInit()
         from rrompy.parameter.parameter_sampling import QuadratureSampler as QS
         self._addParametersToList([], [], ["sampler"],
                                   [QS([[0], [1]], "UNIFORM")])
         del QS
         super().__init__(*args, **kwargs)
         self._postInit()
 
     @property
     def mus(self):
         """Value of mus. Its assignment may reset snapshots."""
         return self._mus
     @mus.setter
     def mus(self, mus):
         mus = checkParameterList(mus, self.npar)[0]
         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 muBounds(self):
         """Value of muBounds."""
         return self.sampler.lims
 
     @property
     def sampler(self):
         """Value of sampler."""
         return self._sampler
     @sampler.setter
     def sampler(self, sampler):
         if 'generatePoints' not in dir(sampler):
             raise RROMPyException("Sampler type not recognized.")
         if hasattr(self, '_sampler') and self._sampler is not None:
             samplerOld = self.sampler
         self._sampler = sampler
         self._approxParameters["sampler"] = self.sampler.__str__()
         if not 'samplerOld' in locals() or samplerOld != self.sampler:
             self.resetSamples()
     
     def setSamples(self, samplingEngine):
         """Copy samplingEngine and samples."""
         self.mus = copy(samplingEngine.mus)
         super().setSamples(samplingEngine)
 
     def computeSnapshots(self):
         """Compute snapshots of solution map."""
         RROMPyAssert(self._mode,
                      message = "Cannot start snapshot computation.")
         self.computeScaleFactor()
-        self.samplingEngine.scaleFactor = self.scaleFactor
+        self.samplingEngine.scaleFactor = self.scaleFactorDer
         if self.samplingEngine.nsamples != self.S:
             vbMng(self, "INIT", "Starting computation of snapshots.", 5)
             self.mus = self.sampler.generatePoints(self.S)
             while len(self.mus) > self.S: self.mus.pop()
             self.samplingEngine.iterSample(self.mus)
             vbMng(self, "DEL", "Done computing snapshots.", 5)
 
     def computeScaleFactor(self):
         """Compute parameter rescaling factor."""
         RROMPyAssert(self._mode, message = "Cannot compute rescaling factor.")
         self.scaleFactor = .5 * np.abs(
                                 self.muBounds[0] ** self.HFEngine.rescalingExp
                               - self.muBounds[1] ** self.HFEngine.rescalingExp)
diff --git a/rrompy/reduction_methods/standard/greedy/generic_greedy_approximant.py b/rrompy/reduction_methods/standard/greedy/generic_greedy_approximant.py
index 2c14df1..0630cbe 100644
--- a/rrompy/reduction_methods/standard/greedy/generic_greedy_approximant.py
+++ b/rrompy/reduction_methods/standard/greedy/generic_greedy_approximant.py
@@ -1,651 +1,655 @@
 # 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 abc import abstractmethod
 from copy import deepcopy as copy
 import numpy as np
 from matplotlib import pyplot as plt
 from rrompy.hfengines.base.linear_affine_engine import checkIfAffine
 from rrompy.reduction_methods.standard.generic_standard_approximant import (
                                                     GenericStandardApproximant)
 from rrompy.utilities.base.types import (Np1D, Np2D, Tuple, List, normEng,
                                          paramVal, paramList, sampList)
 from rrompy.utilities.base import verbosityManager as vbMng
 from rrompy.utilities.numerical import dot
 from rrompy.utilities.expression import expressionEvaluator
 from rrompy.solver import normEngine
 from rrompy.utilities.exception_manager import (RROMPyException, RROMPyAssert,
                                                 RROMPyWarning)
 from rrompy.parameter import checkParameterList, emptyParameterList
 
 __all__ = ['GenericGreedyApproximant']
 
 def localL2Distance(mus:Np2D, badmus:Np2D) -> Np2D:
     return np.linalg.norm(np.tile(mus[..., np.newaxis], [1, 1, len(badmus)])
                         - badmus[..., np.newaxis].T, axis = 1)
 
 def pruneSamples(mus:paramList, badmus:paramList,
                  tol : float = 1e-8) -> Np1D:
     """Remove from mus all the elements which are too close to badmus."""
     if len(badmus) == 0: return mus
     proximity = np.min(localL2Distance(mus.data, badmus.data), axis = 1)
     return np.arange(len(mus))[proximity <= tol]
 
 class GenericGreedyApproximant(GenericStandardApproximant):
     """
     ROM greedy interpolant computation for parametric problems
         (ABSTRACT).
 
     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': number of starting training points;
             - 'sampler': sample point generator;
             - 'greedyTol': uniform error tolerance for greedy algorithm;
                 defaults to 1e-2;
             - 'collinearityTol': collinearity tolerance for greedy algorithm;
                 defaults to 0.;
             - 'maxIter': maximum number of greedy steps; defaults to 1e2;
             - 'nTestPoints': number of test points; defaults to 5e2;
             - 'trainSetGenerator': training sample points generator; defaults
                 to sampler.
             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.
+            - 'POD': whether to compute POD of snapshots;
+            - 'scaleFactorDer': scaling factors for derivative computation;
             - 'greedyTol': uniform error tolerance for greedy algorithm;
             - 'collinearityTol': collinearity tolerance for greedy algorithm;
             - 'maxIter': maximum number of greedy steps;
             - 'nTestPoints': number of test points;
             - 'trainSetGenerator': training sample points generator.
         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 test points.
         sampler: Sample point generator.
         greedyTol: Uniform error tolerance for greedy algorithm.
         collinearityTol: Collinearity tolerance for greedy algorithm.
         maxIter: maximum number of greedy steps.
         nTestPoints: number of starting training points.
         trainSetGenerator: training sample points generator.
         muBounds: list of bounds for parameter values.
         samplingEngine: Sampling engine.
         estimatorNormEngine: Engine for estimator norm computation.
         uHF: High fidelity solution(s) with parameter(s) lastSolvedHF as
             sampleList.
         lastSolvedHF: Parameter(s) corresponding to last computed high fidelity
             solution(s) as parameterList.
         uApproxReduced: Reduced approximate solution(s) with parameter(s)
             lastSolvedApprox as sampleList.
         lastSolvedApproxReduced: Parameter(s) corresponding to last computed
             reduced approximate solution(s) as parameterList.
         uApprox: Approximate solution(s) with parameter(s) lastSolvedApprox as
             sampleList.
         lastSolvedApprox: Parameter(s) corresponding to last computed
             approximate solution(s) as parameterList.
     """
     
     def __init__(self, *args, **kwargs):
         self._preInit()
         self._addParametersToList(["greedyTol", "collinearityTol", "maxIter",
                                    "nTestPoints"], [1e-2, 0., 1e2, 5e2],
                                    ["trainSetGenerator"], ["AUTO"])
         super().__init__(*args, **kwargs)
         self._postInit()
 
     @property
     def greedyTol(self):
         """Value of greedyTol."""
         return self._greedyTol
     @greedyTol.setter
     def greedyTol(self, greedyTol):
         if greedyTol < 0:
             raise RROMPyException("greedyTol must be non-negative.")
         if hasattr(self, "_greedyTol") and self.greedyTol is not None:
             greedyTolold = self.greedyTol
         else:
             greedyTolold = -1
         self._greedyTol = greedyTol
         self._approxParameters["greedyTol"] = self.greedyTol
         if greedyTolold != self.greedyTol:
             self.resetSamples()
 
     @property
     def collinearityTol(self):
         """Value of collinearityTol."""
         return self._collinearityTol
     @collinearityTol.setter
     def collinearityTol(self, collinearityTol):
         if collinearityTol < 0:
             raise RROMPyException("collinearityTol must be non-negative.")
         if (hasattr(self, "_collinearityTol")
         and self.collinearityTol is not None):
             collinearityTolold = self.collinearityTol
         else:
             collinearityTolold = -1
         self._collinearityTol = collinearityTol
         self._approxParameters["collinearityTol"] = self.collinearityTol
         if collinearityTolold != self.collinearityTol:
             self.resetSamples()
 
     @property
     def maxIter(self):
         """Value of maxIter."""
         return self._maxIter
     @maxIter.setter
     def maxIter(self, maxIter):
         if maxIter <= 0: raise RROMPyException("maxIter must be positive.")
         if hasattr(self, "_maxIter") and self.maxIter is not None:
             maxIterold = self.maxIter
         else:
             maxIterold = -1
         self._maxIter = maxIter
         self._approxParameters["maxIter"] = self.maxIter
         if maxIterold != self.maxIter:
             self.resetSamples()
 
     @property
     def nTestPoints(self):
         """Value of nTestPoints."""
         return self._nTestPoints
     @nTestPoints.setter
     def nTestPoints(self, nTestPoints):
         if nTestPoints <= 0:
             raise RROMPyException("nTestPoints must be positive.")
         if not np.isclose(nTestPoints, np.int(nTestPoints)):
             raise RROMPyException("nTestPoints must be an integer.")
         nTestPoints = np.int(nTestPoints)
         if hasattr(self, "_nTestPoints") and self.nTestPoints is not None:
             nTestPointsold = self.nTestPoints
         else:
             nTestPointsold = -1
         self._nTestPoints = nTestPoints
         self._approxParameters["nTestPoints"] = self.nTestPoints
         if nTestPointsold != self.nTestPoints:
             self.resetSamples()
 
     @property
     def trainSetGenerator(self):
         """Value of trainSetGenerator."""
         return self._trainSetGenerator
     @trainSetGenerator.setter
     def trainSetGenerator(self, trainSetGenerator):
         if (isinstance(trainSetGenerator, (str,))
         and trainSetGenerator.upper() == "AUTO"):
             trainSetGenerator = self.sampler
         if 'generatePoints' not in dir(trainSetGenerator):
             raise RROMPyException("trainSetGenerator type not recognized.")
         if (hasattr(self, '_trainSetGenerator')
         and self.trainSetGenerator not in [None, "AUTO"]):
             trainSetGeneratorOld = self.trainSetGenerator
         self._trainSetGenerator = trainSetGenerator
         self._approxParameters["trainSetGenerator"] = self.trainSetGenerator
         if (not 'trainSetGeneratorOld' in locals()
          or trainSetGeneratorOld != self.trainSetGenerator):
             self.resetSamples()
 
     def resetSamples(self):
         """Reset samples."""
         super().resetSamples()
         self._mus = emptyParameterList()
 
     def initEstimatorNormEngine(self, normEngn : normEng = None):
         """Initialize estimator norm engine."""
         if (normEngn is not None or not hasattr(self, "estimatorNormEngine")
                                  or self.estimatorNormEngine is None):
             if normEngn is None:
                 if self.approx_state:
                     if not hasattr(self.HFEngine, "energyNormDualMatrix"):
                         self.HFEngine.buildEnergyNormDualForm()
                     estimatorEnergyMatrix =  self.HFEngine.energyNormDualMatrix
                 else:
                     estimatorEnergyMatrix = self.HFEngine.outputNormMatrix
             else:
                 if hasattr(normEngn, "buildEnergyNormDualForm"):
                     if not hasattr(normEngn, "energyNormDualMatrix"):
                         normEngn.buildEnergyNormDualForm()
                     estimatorEnergyMatrix = normEngn.energyNormDualMatrix
                 else:
                     estimatorEnergyMatrix = normEngn
             self.estimatorNormEngine = normEngine(estimatorEnergyMatrix)
 
     def _affineResidualMatricesContraction(self, rb:Np2D, rA : Np2D = None) \
                                                     -> Tuple[Np1D, Np1D, Np1D]:
         self.assembleReducedResidualBlocks(full = rA is not None)
         # 'ij,jk,ik->k', resbb, radiusb, radiusb.conj()
         ff = np.sum(self.trainedModel.data.resbb.dot(rb) * rb.conj(), axis = 0)
         if rA is None: return ff
         # 'ijk,jkl,il->l', resAb, radiusA, radiusb.conj()
         Lf = np.sum(np.tensordot(self.trainedModel.data.resAb, rA, 2)
                   * rb.conj(), axis = 0)
         # 'ijkl,klt,ijt->t', resAA, radiusA, radiusA.conj()
         LL = np.sum(np.tensordot(self.trainedModel.data.resAA, rA, 2)
                   * rA.conj(), axis = (0, 1))
         return ff, Lf, LL
 
     def getErrorEstimatorAffine(self, mus:Np1D) -> Np1D:
         """Standard residual estimator."""
         checkIfAffine(self.HFEngine, "apply affinity-based error estimator")
         self.HFEngine.buildA()
         self.HFEngine.buildb()
         mus = checkParameterList(mus, self.npar)[0]
         verb = self.trainedModel.verbosity
         self.trainedModel.verbosity = 0
         uApproxRs = self.getApproxReduced(mus)
         muTestEff = mus ** self.HFEngine.rescalingExp
         radiusA = np.empty((len(self.HFEngine.thAs), len(mus)),
                            dtype = np.complex)
         radiusb = np.empty((len(self.HFEngine.thbs), len(mus)),
                            dtype = np.complex)
         for j, thA in enumerate(self.HFEngine.thAs):
             radiusA[j] = expressionEvaluator(thA[0], muTestEff)
         for j, thb in enumerate(self.HFEngine.thbs):
             radiusb[j] = expressionEvaluator(thb[0], muTestEff)
         radiusA = np.expand_dims(uApproxRs.data, 1) * radiusA
         ff, Lf, LL = self._affineResidualMatricesContraction(radiusb, radiusA)
         err = np.abs((LL - 2. * np.real(Lf) + ff) / ff) ** .5
         self.trainedModel.verbosity = verb
         return err
 
     def errorEstimator(self, mus:Np1D, return_max : bool = False) -> Np1D:
         setupOK = self.setupApproxLocal()
         if setupOK > 0:
             err = np.empty(len(mus))
             err[:] = np.nan
             if not return_max: return err
             return err, [- setupOK], np.nan
         mus = checkParameterList(mus, self.npar)[0]
         vbMng(self.trainedModel, "INIT",
               "Evaluating error estimator at mu = {}.".format(mus), 10)
         err = self.getErrorEstimatorAffine(mus)
         vbMng(self.trainedModel, "DEL", "Done evaluating error estimator", 10)
         if not return_max: return err
         idxMaxEst = [np.argmax(err)]
         return err, idxMaxEst, err[idxMaxEst]
 
     def _isLastSampleCollinear(self) -> bool:
         """Check collinearity of last sample."""
         if self.collinearityTol <= 0.: return False
         if self.POD:
             reff = self.samplingEngine.RPOD[:, -1]
         else:
             RROMPyWarning(("Repeated orthogonalization of the samples for "
                            "collinearity check. Consider setting POD to "
                            "True."))
             if not hasattr(self, "_PODEngine"):
                 from rrompy.sampling.base.pod_engine import PODEngine
                 self._PODEngine = PODEngine(self.HFEngine)
             reff = self._PODEngine.generalizedQR(self.samplingEngine.samples,
                                                  only_R = True,
                                                  is_state = True)[:, -1]
         cLevel = np.abs(reff[-1]) / np.linalg.norm(reff)
         cLevel = np.inf if np.isclose(cLevel, 0.) else cLevel ** -1.
         vbMng(self, "MAIN", "Collinearity indicator {:.4e}.".format(cLevel), 3)
         return cLevel > self.collinearityTol
 
     def plotEstimator(self, est:Np1D, idxMax:List[int], estMax:List[float]):
         if not (np.any(np.isnan(est)) or np.any(np.isinf(est))):
             fig = plt.figure(figsize = plt.figaspect(1. / self.npar))
             for jpar in range(self.npar):
                 ax = fig.add_subplot(1, self.npar, 1 + jpar)
                 musre = copy(self.muTest.re.data)
                 errCP = copy(est)
                 idx = np.delete(np.arange(self.npar), jpar)
                 while len(musre) > 0:
                     if self.npar == 1:
                         currIdx = np.arange(len(musre))
                     else:
                         currIdx = np.where(np.isclose(np.sum(
                              np.abs(musre[:, idx] - musre[0, idx]), 1), 0.))[0]
                     ax.semilogy(musre[currIdx, jpar], errCP[currIdx], 'k',
                                 linewidth = 1)
                     musre = np.delete(musre, currIdx, 0)
                     errCP = np.delete(errCP, currIdx)
                 ax.semilogy([self.muBounds.re(0, jpar),
                              self.muBounds.re(-1, jpar)],
                             [self.greedyTol] * 2, 'r--')
                 ax.semilogy(self.mus.re(jpar),
                             2. * self.greedyTol * np.ones(len(self.mus)), '*m')
                 if len(idxMax) > 0 and estMax is not None:
                     ax.semilogy(self.muTest.re(idxMax, jpar), estMax, 'xr')
                 ax.grid()
             plt.tight_layout()
             plt.show()
     
     def greedyNextSample(self, muidx:int, plotEst : str = "NONE")\
                                           -> Tuple[Np1D, int, float, paramVal]:
         """Compute next greedy snapshot of solution map."""
         RROMPyAssert(self._mode, message = "Cannot add greedy sample.")
         mus = copy(self.muTest[muidx])
         self.muTest.pop(muidx)
         for j, mu in enumerate(mus):
             vbMng(self, "MAIN",
                   ("Adding sample point no. {} at {} to training "
                    "set.").format(len(self.mus) + 1, mu), 3)
             self.mus.append(mu)
             self._S = len(self.mus)
             self._approxParameters["S"] = self.S
             if (self.samplingEngine.nsamples <= len(mus) - j - 1
              or not np.allclose(mu,
                                 self.samplingEngine.mus.data[j - len(mus)])):
                 self.samplingEngine.nextSample(mu)
             if self._isLastSampleCollinear():
                 vbMng(self, "MAIN",
                       ("Collinearity above tolerance detected. Starting "
                        "preemptive greedy loop termination."), 3)
                 self._collinearityFlag = 1
                 errorEstTest = np.empty(len(self.muTest))
                 errorEstTest[:] = np.nan
                 return errorEstTest, [-1], np.nan, np.nan
         errorEstTest, muidx, maxErrorEst = self.errorEstimator(self.muTest,
                                                                True)
         if plotEst == "ALL":
             self.plotEstimator(errorEstTest, muidx, maxErrorEst)
         return errorEstTest, muidx, maxErrorEst, self.muTest[muidx]
 
     def _preliminaryTraining(self):
         """Initialize starting snapshots of solution map."""
         RROMPyAssert(self._mode, message = "Cannot start greedy algorithm.")
         self.computeScaleFactor()
         if self.samplingEngine.nsamples > 0: return
         self.resetSamples()
-        self.samplingEngine.scaleFactor = self.scaleFactor
+        self.samplingEngine.scaleFactor = self.scaleFactorDer
         self.mus = self.trainSetGenerator.generatePoints(self.S)
         while len(self.mus) > self.S: self.mus.pop()
         muTestBase = self.sampler.generatePoints(self.nTestPoints, False)
         idxPop = pruneSamples(muTestBase ** self.HFEngine.rescalingExp,
                               self.mus ** self.HFEngine.rescalingExp,
                               1e-10 * self.scaleFactor[0])
         muTestBase.pop(idxPop)
         muLast = copy(self.mus[-1])
         self.mus.pop()
         if len(self.mus) > 0:
             vbMng(self, "MAIN", 
                   ("Adding first {} sample point{} at {} to training "
                    "set.").format(self.S - 1, "" + "s" * (self.S > 2),
                                   self.mus), 3)
             self.samplingEngine.iterSample(self.mus)
         self._S = len(self.mus)
         self._approxParameters["S"] = self.S
         self.muTest = emptyParameterList()
         self.muTest.reset((len(muTestBase) + 1, self.mus.shape[1]))
         self.muTest[: -1] = muTestBase.data
         self.muTest[-1] = muLast.data
 
     @abstractmethod
     def setupApproxLocal(self) -> int:
         if self.checkComputedApprox(): return -1
         RROMPyAssert(self._mode, message = "Cannot setup approximant.")
         vbMng(self, "INIT", "Setting up local approximant.", 5)
         pass
         vbMng(self, "DEL", "Done setting up local approximant.", 5)
         return 0
 
     def setupApprox(self, plotEst : str = "NONE") -> int:
         """Compute greedy snapshots of solution map."""
         if self.checkComputedApprox(): return -1
         RROMPyAssert(self._mode, message = "Cannot start greedy algorithm.")
         vbMng(self, "INIT", "Starting computation of snapshots.", 3)
         self._collinearityFlag = 0
         self._preliminaryTraining()
         muidx, firstGreedyIter = [len(self.muTest) - 1], True
         errorEstTest, maxErrorEst = [np.inf], np.inf
         max2ErrorEst, trainedModelOld = np.inf, None
         while firstGreedyIter or (len(self.muTest) > 0
                      and (maxErrorEst is None or max2ErrorEst > self.greedyTol)
                      and self.samplingEngine.nsamples < self.maxIter):
             muTestOld, errorEstTestOld = self.muTest, errorEstTest
             muidxOld, maxErrorEstOld = muidx, maxErrorEst
             errorEstTest, muidx, maxErrorEst, mu = self.greedyNextSample(
                                                             muidx, plotEst)
             if maxErrorEst is not None and (np.any(np.isnan(maxErrorEst))
                                          or np.any(np.isinf(maxErrorEst))):
                 if self._collinearityFlag == 0 and not firstGreedyIter:
                     RROMPyWarning(("Instability in a posteriori "
                                    "estimator. Starting preemptive greedy "
                                    "loop termination."))
                 self.muTest, errorEstTest = muTestOld, errorEstTestOld
                 if firstGreedyIter:
                     self.mus.pop(-1)
                     self.samplingEngine.popSample()
                     if muidx[0] < 0:
                         self.trainedModel = None
                         raise RROMPyException(("Instability in approximant "
                                                "computation. Aborting greedy "
                                                "iterations."))
                 else:
                     self._approxParameters = (
                                          trainedModelOld.data.approxParameters)
                     self._S = trainedModelOld.data.approxParameters["S"]
                     self._approxParameters["S"] = self.S
                     self.trainedModel.data = copy(trainedModelOld.data)
                 muidx, maxErrorEst = muidxOld, maxErrorEstOld
                 break
             if maxErrorEst is not None:
                 max2ErrorEst = np.max(maxErrorEst)
                 vbMng(self, "MAIN", ("Uniform testing error estimate "
                                      "{:.4e}.").format(max2ErrorEst), 3)
             if firstGreedyIter:
                 trainedModelOld = copy(self.trainedModel)
             else:
                 trainedModelOld.data = copy(self.trainedModel.data)
             firstGreedyIter = False
         if (maxErrorEst is None or max2ErrorEst <= self.greedyTol
          or np.any(np.isnan(maxErrorEst)) or np.any(np.isinf(maxErrorEst))):
             while self.samplingEngine.nsamples > self.S:
                 self.samplingEngine.popSample()
             while len(self.mus) > self.S: self.mus.pop(-1)
         else:
             self._S = self.samplingEngine.nsamples
             self._approxParameters["S"] = self.S
             while len(self.mus) < self.S:
                 self.mus.append(self.samplingEngine.mus[len(self.mus)])
             self.setupApproxLocal()
         if plotEst == "LAST":
             self.plotEstimator(errorEstTest, muidx, maxErrorEst)
         vbMng(self, "DEL", 
               ("Done computing snapshots (final snapshot count: "
                "{}).").format(self.samplingEngine.nsamples), 3)
         return 0
 
     def checkComputedApprox(self) -> bool:
         """
         Check if setup of new approximant is not needed.
 
         Returns:
             True if new setup is not needed. False otherwise.
         """
         return (super().checkComputedApprox()
             and len(self.mus) == self.trainedModel.data.projMat.shape[1])
 
     def assembleReducedResidualGramian(self, pMat:sampList):
         """
         Build residual gramian of reduced linear system through projections.
         """
         self.initEstimatorNormEngine()
         if (not hasattr(self.trainedModel.data, "gramian")
          or self.trainedModel.data.gramian is None):
             gramian = self.estimatorNormEngine.innerProduct(pMat, pMat)
         else:
             Sold = self.trainedModel.data.gramian.shape[0]
             S = len(self.mus)
             if Sold > S:
                 gramian = self.trainedModel.data.gramian[: S, : S]
             else:
                 idxOld = list(range(Sold))
                 idxNew = list(range(Sold, S))
                 gramian = np.empty((S, S), dtype = np.complex)
                 gramian[: Sold, : Sold] = self.trainedModel.data.gramian
                 gramian[: Sold, Sold :] = (
                            self.estimatorNormEngine.innerProduct(pMat(idxNew),
                                                                  pMat(idxOld)))
                 gramian[Sold :, : Sold] = gramian[: Sold, Sold :].T.conj()
                 gramian[Sold :, Sold :] = (
                            self.estimatorNormEngine.innerProduct(pMat(idxNew),
                                                                  pMat(idxNew)))
         self.trainedModel.data.gramian = gramian
 
     def assembleReducedResidualBlocksbb(self, bs:List[Np1D]):
         """
         Build blocks (of type bb) of reduced linear system through projections.
         """
         self.initEstimatorNormEngine()
         nbs = len(bs)
         if (not hasattr(self.trainedModel.data, "resbb")
          or self.trainedModel.data.resbb is None):
             resbb = np.empty((nbs, nbs), dtype = np.complex)
             for i in range(nbs):
                 Mbi = bs[i]
                 resbb[i, i] = self.estimatorNormEngine.innerProduct(Mbi, Mbi)
                 for j in range(i):
                     Mbj = bs[j]
                     resbb[i, j] = self.estimatorNormEngine.innerProduct(Mbj,
                                                                         Mbi)
             for i in range(nbs):
                 for j in range(i + 1, nbs):
                     resbb[i, j] = resbb[j, i].conj()
             self.trainedModel.data.resbb = resbb
 
     def assembleReducedResidualBlocksAb(self, As:List[Np2D], bs:List[Np1D],
                                         pMat:sampList):
         """
         Build blocks (of type Ab) of reduced linear system through projections.
         """
         self.initEstimatorNormEngine()
         nAs = len(As)
         nbs = len(bs)
         S = len(self.mus)
         if (not hasattr(self.trainedModel.data, "resAb")
          or self.trainedModel.data.resAb is None):
             if not isinstance(pMat, (np.ndarray,)): pMat = pMat.data
             resAb = np.empty((nbs, S, nAs), dtype = np.complex)
             for j in range(nAs):
                 MAj = dot(As[j], pMat)
                 for i in range(nbs):
                     Mbi = bs[i]
                     resAb[i, :, j] = self.estimatorNormEngine.innerProduct(MAj,
                                                                            Mbi)
         else:
             Sold = self.trainedModel.data.resAb.shape[1]
             if Sold == S: return
             if Sold > S:
                 resAb = self.trainedModel.data.resAb[:, : S, :]
             else:
                 if not isinstance(pMat, (np.ndarray,)): pMat = pMat.data
                 resAb = np.empty((nbs, S, nAs), dtype = np.complex)
                 resAb[:, : Sold, :] = self.trainedModel.data.resAb
                 for j in range(nAs):
                     MAj = dot(As[j], pMat[:, Sold :])
                     for i in range(nbs):
                         Mbi = bs[i]
                         resAb[i, Sold :, j] = (
                                self.estimatorNormEngine.innerProduct(MAj, Mbi))
         self.trainedModel.data.resAb = resAb
 
     def assembleReducedResidualBlocksAA(self, As:List[Np2D], pMat:sampList):
         """
         Build blocks (of type AA) of reduced linear system through projections.
         """
         self.initEstimatorNormEngine()
         nAs = len(As)
         S = len(self.mus)
         if (not hasattr(self.trainedModel.data, "resAA")
          or self.trainedModel.data.resAA is None):
             if not isinstance(pMat, (np.ndarray,)): pMat = pMat.data
             resAA = np.empty((S, nAs, S, nAs), dtype = np.complex)
             for i in range(nAs):
                 MAi = dot(As[i], pMat)
                 resAA[:, i, :, i] = (
                                self.estimatorNormEngine.innerProduct(MAi, MAi))
                 for j in range(i):
                     MAj = dot(As[j], pMat)
                     resAA[:, i, :, j] = (
                                self.estimatorNormEngine.innerProduct(MAj, MAi))
             for i in range(nAs):
                 for j in range(i + 1, nAs):
                     resAA[:, i, :, j] = resAA[:, j, :, i].T.conj()
         else:
             Sold = self.trainedModel.data.resAA.shape[0]
             if Sold == S: return
             if Sold > S:
                 resAA = self.trainedModel.data.resAA[: S, :, : S, :]
             else:
                 if not isinstance(pMat, (np.ndarray,)): pMat = pMat.data
                 resAA = np.empty((S, nAs, S, nAs), dtype = np.complex)
                 resAA[: Sold, :, : Sold, :] = self.trainedModel.data.resAA
                 for i in range(nAs):
                     MAi = dot(As[i], pMat)
                     resAA[: Sold, i, Sold :, i] = (
                          self.estimatorNormEngine.innerProduct(MAi[:, Sold :],
                                                                MAi[:, : Sold]))
                     resAA[Sold :, i, : Sold, i] = resAA[: Sold, i,
                                                         Sold :, i].T.conj()
                     resAA[Sold :, i, Sold :, i] = (
                          self.estimatorNormEngine.innerProduct(MAi[:, Sold :],
                                                                MAi[:, Sold :]))
                     for j in range(i):
                         MAj = dot(As[j], pMat)
                         resAA[: Sold, i, Sold :, j] = (
                          self.estimatorNormEngine.innerProduct(MAj[:, Sold :],
                                                                MAi[:, : Sold]))
                         resAA[Sold :, i, : Sold, j] = (
                          self.estimatorNormEngine.innerProduct(MAj[:, : Sold],
                                                                MAi[:, Sold :]))
                         resAA[Sold :, i, Sold :, j] = (
                          self.estimatorNormEngine.innerProduct(MAj[:, Sold :],
                                                                MAi[:, Sold :]))
                 for i in range(nAs):
                     for j in range(i + 1, nAs):
                         resAA[: Sold, i, Sold :, j] = (
                                           resAA[Sold :, j, : Sold, i].T.conj())
                         resAA[Sold :, i, : Sold, j] = (
                                           resAA[: Sold, j, Sold :, i].T.conj())
                         resAA[Sold :, i, Sold :, j] = (
                                           resAA[Sold :, j, Sold :, i].T.conj())
         self.trainedModel.data.resAA = resAA
 
     def assembleReducedResidualBlocks(self, full : bool = False):
         """Build affine blocks of affine decomposition of residual."""
         if full:
             checkIfAffine(self.HFEngine, "assemble reduced residual blocks")
         else:
             checkIfAffine(self.HFEngine, "assemble reduced RHS blocks", True)
         self.HFEngine.buildb()
         self.assembleReducedResidualBlocksbb(self.HFEngine.bs)
         if full:
             pMat = self.samplingEngine.samples
             self.HFEngine.buildA()
             self.assembleReducedResidualBlocksAb(self.HFEngine.As,
                                                  self.HFEngine.bs, pMat)
             self.assembleReducedResidualBlocksAA(self.HFEngine.As, pMat)
diff --git a/rrompy/reduction_methods/standard/greedy/rational_interpolant_greedy.py b/rrompy/reduction_methods/standard/greedy/rational_interpolant_greedy.py
index 84b2b99..068b5e1 100644
--- a/rrompy/reduction_methods/standard/greedy/rational_interpolant_greedy.py
+++ b/rrompy/reduction_methods/standard/greedy/rational_interpolant_greedy.py
@@ -1,548 +1,552 @@
 # 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.hfengines.base.linear_affine_engine import checkIfAffine
 from .generic_greedy_approximant import GenericGreedyApproximant
 from rrompy.utilities.poly_fitting.polynomial import (polybases,
                                                   PolynomialInterpolator as PI,
                                                   polyvanderTotal as pvT)
 from rrompy.utilities.numerical import dot
 from rrompy.utilities.numerical.degree import totalDegreeN
 from rrompy.utilities.expression import expressionEvaluator
 from rrompy.reduction_methods.standard import RationalInterpolant
 from rrompy.utilities.base.types import Np1D, Tuple, paramVal, List
 from rrompy.utilities.base import verbosityManager as vbMng
 from rrompy.utilities.poly_fitting import customFit
 from rrompy.utilities.exception_manager import (RROMPyWarning, RROMPyException,
                                                 RROMPyAssert, RROMPy_FRAGILE)
 from rrompy.parameter import checkParameterList
 from rrompy.sampling import sampleList, emptySampleList
 
 __all__ = ['RationalInterpolantGreedy']
 
 class RationalInterpolantGreedy(GenericGreedyApproximant, RationalInterpolant):
     """
     ROM greedy rational interpolant computation for parametric problems.
 
     Args:
         HFEngine: HF problem solver.
         mu0(optional): Default parameter. Defaults to 0.
         approxParameters(optional): Dictionary containing values for main
             parameters of approximant. Recognized keys are:
             - 'POD': whether to compute POD of snapshots; defaults to True;
+            - 'scaleFactorDer': scaling factors for derivative computation;
+                defaults to 'AUTO';
             - 'S': number of starting training points;
             - 'sampler': sample point generator;
             - 'greedyTol': uniform error tolerance for greedy algorithm;
                 defaults to 1e-2;
             - 'collinearityTol': collinearity tolerance for greedy algorithm;
                 defaults to 0.;
             - 'maxIter': maximum number of greedy steps; defaults to 1e2;
             - 'nTestPoints': number of test points; defaults to 5e2;
             - 'trainSetGenerator': training sample points generator; defaults
                 to sampler;
             - 'polybasis': type of basis for interpolation; defaults to
                 'MONOMIAL';
             - 'errorEstimatorKind': kind of error estimator; available values
                 include 'AFFINE', 'DISCREPANCY', 'LOOK_AHEAD',
                 'LOOK_AHEAD_RES', 'LOOK_AHEAD_OUTPUT', and 'NONE'; defaults to
                 'NONE';
             - '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 and must
             be True.
         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.
+            - 'POD': whether to compute POD of snapshots;
+            - 'scaleFactorDer': scaling factors for derivative computation;
             - 'greedyTol': uniform error tolerance for greedy algorithm;
             - 'collinearityTol': collinearity tolerance for greedy algorithm;
             - 'maxIter': maximum number of greedy steps;
             - 'nTestPoints': number of test points;
             - 'trainSetGenerator': training sample points generator;
             - 'errorEstimatorKind': kind of error estimator;
             - 'interpRcond': tolerance for interpolation;
             - '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 test points.
         sampler: Sample point generator.
         greedyTol: uniform error tolerance for greedy algorithm.
         collinearityTol: Collinearity tolerance for greedy algorithm.
         maxIter: maximum number of greedy steps.
         nTestPoints: number of starting training points.
         trainSetGenerator: training sample points generator.
         robustTol: tolerance for robust rational denominator management.
         errorEstimatorKind: kind of error estimator.
         interpRcond: tolerance for interpolation.
         robustTol: tolerance for robust rational denominator management.
         muBounds: list of bounds for parameter values.
         samplingEngine: Sampling engine.
         estimatorNormEngine: Engine for estimator norm computation.
         uHF: High fidelity solution(s) with parameter(s) lastSolvedHF as
             sampleList.
         lastSolvedHF: Parameter(s) corresponding to last computed high fidelity
             solution(s) as parameterList.
         uApproxReduced: Reduced approximate solution(s) with parameter(s)
             lastSolvedApprox as sampleList.
         lastSolvedApproxReduced: Parameter(s) corresponding to last computed
             reduced approximate solution(s) as parameterList.
         uApprox: Approximate solution(s) with parameter(s) lastSolvedApprox as
             sampleList.
         lastSolvedApprox: Parameter(s) corresponding to last computed
             approximate solution(s) as parameterList.
     """
     
     _allowedEstimatorKinds = ["AFFINE", "DISCREPANCY", "LOOK_AHEAD",
                               "LOOK_AHEAD_RES", "LOOK_AHEAD_OUTPUT", "NONE"]
     
     def __init__(self, *args, **kwargs):
         self._preInit()
         self._addParametersToList(["errorEstimatorKind"], ["DISCREPANCY"],
                                   toBeExcluded = ["M", "N", "polydegreetype", 
                                                   "radialDirectionalWeights",
                                                   "nNearestNeighbor"])
         super().__init__(*args, **kwargs)
         if not self.approx_state and self.errorEstimatorKind not in [
                                     "LOOK_AHEAD", "LOOK_AHEAD_OUTPUT", "NONE"]:
             raise RROMPyException(("Must compute greedy approximation of "
                                    "state, unless error estimator allows "
                                    "otherwise."))
         self._postInit()
 
     @property
     def approx_state(self):
         """Value of approx_state."""
         return self._approx_state
     @approx_state.setter
     def approx_state(self, approx_state):
         RationalInterpolant.approx_state.fset(self, approx_state)
         if (not self.approx_state and hasattr(self, "_errorEstimatorKind")
                                   and self.errorEstimatorKind not in [
                                    "LOOK_AHEAD", "LOOK_AHEAD_OUTPUT", "NONE"]):
             raise RROMPyException(("Must compute greedy approximation of "
                                    "state, unless error estimator allows "
                                    "otherwise."))
 
     @property
     def E(self):
         """Value of E."""
         self._E = self.sampleBatchIdx - 1
         return self._E
     @E.setter
     def E(self, E):
         RROMPyWarning(("E is used just to simplify inheritance, and its value "
                        "cannot be changed from that of sampleBatchIdx - 1."))
         
     def _setMAuto(self):
         self.M = self.E
 
     def _setNAuto(self):
         self.N = self.E
 
     @property
     def polydegreetype(self):
         """Value of polydegreetype."""
         return "TOTAL"
     @polydegreetype.setter
     def polydegreetype(self, polydegreetype):
         RROMPyWarning(("polydegreetype is used just to simplify inheritance, "
                        "and its value cannot be changed from 'TOTAL'."))
         
     @property
     def polybasis(self):
         """Value of polybasis."""
         return self._polybasis
     @polybasis.setter
     def polybasis(self, polybasis):
         try:
             polybasis = polybasis.upper().strip().replace(" ","")
             if polybasis not in polybases: 
                 raise RROMPyException("Sample type not recognized.")
             self._polybasis = polybasis
         except:
             RROMPyWarning(("Prescribed polybasis not recognized. Overriding "
                            "to 'MONOMIAL'."))
             self._polybasis = "MONOMIAL"
         self._approxParameters["polybasis"] = self.polybasis
 
     @property
     def errorEstimatorKind(self):
         """Value of errorEstimatorKind."""
         return self._errorEstimatorKind
     @errorEstimatorKind.setter
     def errorEstimatorKind(self, errorEstimatorKind):
         errorEstimatorKind = errorEstimatorKind.upper()
         if errorEstimatorKind not in self._allowedEstimatorKinds:
             RROMPyWarning(("Error estimator kind not recognized. Overriding "
                            "to 'NONE'."))
             errorEstimatorKind = "NONE"
         self._errorEstimatorKind = errorEstimatorKind
         self._approxParameters["errorEstimatorKind"] = self.errorEstimatorKind
         if (self.errorEstimatorKind not in [
                                      "LOOK_AHEAD", "LOOK_AHEAD_OUTPUT", "NONE"]
         and hasattr(self, "_approx_state") and not self.approx_state):
             raise RROMPyException(("Must compute greedy approximation of "
                                    "state, unless error estimator allows "
                                    "otherwise."))
 
     def _polyvanderAuxiliary(self, mus, deg, *args):
         return pvT(mus, deg, *args)
 
     def getErrorEstimatorDiscrepancy(self, mus:Np1D) -> Np1D:
         """Discrepancy-based residual estimator."""
         checkIfAffine(self.HFEngine, "apply discrepancy-based error estimator")
         mus = checkParameterList(mus, self.npar)[0]
         muCTest = self.trainedModel.centerNormalize(mus)
         verb = self.trainedModel.verbosity
         self.trainedModel.verbosity = 0
         QTest = self.trainedModel.getQVal(mus)
         QTzero = np.where(QTest == 0.)[0]
         if len(QTzero) > 0:
             RROMPyWarning(("Adjusting estimator to avoid division by "
                            "numerically zero denominator."))
             QTest[QTzero] = np.finfo(np.complex).eps / (1. + self.N)
         self.HFEngine.buildA()
         self.HFEngine.buildb()
         nAs, nbs = self.HFEngine.nAs, self.HFEngine.nbs
         muTrainEff = self.mus ** self.HFEngine.rescalingExp
         muTestEff = mus ** self.HFEngine.rescalingExp
         PTrain = self.trainedModel.getPVal(self.mus).data.T
         QTrain = self.trainedModel.getQVal(self.mus)
         QTzero = np.where(QTrain == 0.)[0]
         if len(QTzero) > 0:
             RROMPyWarning(("Adjusting estimator to avoid division by "
                            "numerically zero denominator."))
             QTrain[QTzero] = np.finfo(np.complex).eps / (1. + self.N)
         PTest = self.trainedModel.getPVal(mus).data
         radiusAbTrain = np.empty((self.S, nAs * self.S + nbs),
                                  dtype = np.complex)
         radiusA = np.empty((self.S, nAs, len(mus)), dtype = np.complex)
         radiusb = np.empty((nbs, len(mus)), dtype = np.complex)
         for j, thA in enumerate(self.HFEngine.thAs):
             idxs = j * self.S + np.arange(self.S)
             radiusAbTrain[:, idxs] = expressionEvaluator(thA[0], muTrainEff,
                                                          (self.S, 1)) * PTrain
             radiusA[:, j] = PTest * expressionEvaluator(thA[0], muTestEff,
                                                         (len(mus),))
         for j, thb in enumerate(self.HFEngine.thbs):
             idx = nAs * self.S + j
             radiusAbTrain[:, idx] = QTrain * expressionEvaluator(thb[0],
                                                          muTrainEff, (self.S,))
             radiusb[j] = QTest * expressionEvaluator(thb[0], muTestEff,
                                                      (len(mus),))
         QRHSNorm2 = self._affineResidualMatricesContraction(radiusb)
         vanTrain = self._polyvanderAuxiliary(self._musUniqueCN, self.E,
                                              self.polybasis0, self._derIdxs,
                                              self._reorder)
         interpPQ = customFit(vanTrain, radiusAbTrain,
                              rcond = self.interpRcond)
         vanTest = self._polyvanderAuxiliary(muCTest, self.E,
                                             self.polybasis0)
         DradiusAb = vanTest.dot(interpPQ)
         radiusA = (radiusA
                  - DradiusAb[:, : - nbs].reshape(len(mus), -1, self.S).T)
         radiusb = radiusb - DradiusAb[:, - nbs :].T
         ff, Lf, LL = self._affineResidualMatricesContraction(radiusb, radiusA)
         err = np.abs((LL - 2. * np.real(Lf) + ff) / QRHSNorm2) ** .5
         self.trainedModel.verbosity = verb
         return err
     
     def getErrorEstimatorLookAhead(self, mus:Np1D,
                                    what : str = "") -> Tuple[Np1D, List[int]]:
         """Residual estimator based on look-ahead idea."""
         errTest, QTest, idxMaxEst = self._EIMStep(mus)
         _approx_state_old = self.approx_state
         if what == "OUTPUT" and _approx_state_old: self._approx_state = False
         self.initEstimatorNormEngine()
         self._approx_state = _approx_state_old
         mu_muTestSample = mus[idxMaxEst]
         app_muTestSample = self.getApproxReduced(mu_muTestSample)
         if self._mode == RROMPy_FRAGILE:
             if what == "RES" and not self.HFEngine.isCEye:
                 raise RROMPyException(("Cannot compute LOOK_AHEAD_RES "
                                        "estimator in fragile mode for "
                                        "non-scalar C."))
             app_muTestSample = dot(self.trainedModel.data.projMat[:,
                                                   : app_muTestSample.shape[0]],
                                    app_muTestSample.data)
         else:
             app_muTestSample = dot(self.samplingEngine.samples,
                                    app_muTestSample)
         if what == "RES":
             errmu = self.HFEngine.residual(mu_muTestSample, app_muTestSample,
                                post_c = False)
             solmu = self.HFEngine.residual(mu_muTestSample, None,
                                            post_c = False)
         else:
             for j, mu in enumerate(mu_muTestSample):
                 uEx = self.samplingEngine.nextSample(mu)
                 if hasattr(self.samplingEngine, "samples_full"):
                     uEx = self.samplingEngine.samples_full[-1]
                 if j == 0:
                     solmu = emptySampleList()
                     solmu.reset((len(uEx), len(mu_muTestSample)),
                                 dtype = uEx.dtype)
                 solmu[j] = uEx
             if what == "OUTPUT" and self.approx_state:
                 solmu = sampleList(self.HFEngine.applyC(solmu.data))
                 app_muTestSample = sampleList(self.HFEngine.applyC(
                                                         app_muTestSample.data))
             errmu = solmu - app_muTestSample
         errsamples = (self.estimatorNormEngine.norm(errmu)
                     / self.estimatorNormEngine.norm(solmu))
         musT = copy(self.mus)
         musT.append(mu_muTestSample)
         musT = self.trainedModel.centerNormalize(musT)
         musC = self.trainedModel.centerNormalize(mus)
         errT = np.zeros(len(musT), dtype = np.complex)
         err = np.zeros(len(mus))
         for l in range(len(mu_muTestSample)):
             errT[len(self.mus) + l] = errsamples[l] * QTest[idxMaxEst[l]]
             p = PI()
             wellCond, msg = p.setupByInterpolation(musT, errT, self.E + 1,
                                          self.polybasis, self.verbosity >= 15,
                                          True, {}, {"rcond": self.interpRcond})
             err += np.abs(p(musC))
             errT[len(self.mus) + l] = 0.
         err /= QTest
         vbMng(self, "MAIN", msg, 15)
         return err, idxMaxEst
     
     def getErrorEstimatorNone(self, mus:Np1D) -> Np1D:
         """EIM-based residual estimator."""
         err = np.max(self._EIMStep(mus, True), axis = 1)
         err *= self.greedyTol / np.mean(err)
         return err
 
     def _EIMStep(self, mus:Np1D,
                  only_one : bool = False) -> Tuple[Np1D, Np1D, List[int]]:
         """Residual estimator based on look-ahead idea."""
         mus = checkParameterList(mus, self.npar)[0]
         verb = self.trainedModel.verbosity
         self.trainedModel.verbosity = 0
         QTest = self.trainedModel.getQVal(mus)
         QTzero = np.where(QTest == 0.)[0]
         if len(QTzero) > 0:
             RROMPyWarning(("Adjusting estimator to avoid division by "
                            "numerically zero denominator."))
             QTest[QTzero] = np.finfo(np.complex).eps / (1. + self.N)
         QTest = np.abs(QTest)
         muCTest = self.trainedModel.centerNormalize(mus)
         muCTrain = self.trainedModel.centerNormalize(self.mus)
         vanTest = self._polyvanderAuxiliary(muCTest, self.E, self.polybasis)
         vanTestNext = self._polyvanderAuxiliary(muCTest, self.E + 1,
                                                 self.polybasis)[:,
                                                             vanTest.shape[1] :]
         idxsTest = np.arange(vanTestNext.shape[1])
         basis = np.zeros((len(idxsTest), 0), dtype = float)
         idxMaxEst = []
         while len(idxsTest) > 0:
             vanTrial = self._polyvanderAuxiliary(muCTrain, self.E,
                                                  self.polybasis)
             vanTrialNext = self._polyvanderAuxiliary(muCTrain, self.E + 1,
                                                      self.polybasis)[:,
                                                            vanTrial.shape[1] :]
             vanTrial = np.hstack((vanTrial, vanTrialNext.dot(basis).reshape(
                                            len(vanTrialNext), basis.shape[1])))
             valuesTrial = vanTrialNext[:, idxsTest]
             vanTestEff = np.hstack((vanTest, vanTestNext.dot(basis).reshape(
                                             len(vanTestNext), basis.shape[1])))
             vanTestNextEff = vanTestNext[:, idxsTest]
             try:
                 coeffTest = np.linalg.solve(vanTrial, valuesTrial)
             except np.linalg.LinAlgError as e:
                 raise RROMPyException(e)
             errTest = (np.abs(vanTestNextEff - vanTestEff.dot(coeffTest))
                      / np.expand_dims(QTest, 1))
             if only_one:
                 self.trainedModel.verbosity = verb
                 return errTest
             idxMaxErr = np.unravel_index(np.argmax(errTest), errTest.shape)
             idxMaxEst += [idxMaxErr[0]]
             muCTrain.append(muCTest[idxMaxErr[0]])
             basis = np.pad(basis, [(0, 0), (0, 1)], "constant")
             basis[idxsTest[idxMaxErr[1]], -1] = 1.
             idxsTest = np.delete(idxsTest, idxMaxErr[1])
         self.trainedModel.verbosity = verb
         return errTest, QTest, idxMaxEst
     
     def errorEstimator(self, mus:Np1D, return_max : bool = False) -> Np1D:
         """Standard residual-based error estimator."""
         setupOK = self.setupApproxLocal()
         if setupOK > 0:
             err = np.empty(len(mus))
             err[:] = np.nan
             if not return_max: return err
             return err, [- setupOK], np.nan
         mus = checkParameterList(mus, self.npar)[0]
         vbMng(self.trainedModel, "INIT",
               "Evaluating error estimator at mu = {}.".format(mus), 10)
         if self.errorEstimatorKind == "AFFINE":
             err = self.getErrorEstimatorAffine(mus)
         else:
             self._setupInterpolationIndices()
             if self.errorEstimatorKind == "DISCREPANCY":
                 err = self.getErrorEstimatorDiscrepancy(mus)
             elif self.errorEstimatorKind[: 10] == "LOOK_AHEAD":
                 err, idxMaxEst = self.getErrorEstimatorLookAhead(mus,
                                                  self.errorEstimatorKind[11 :])
             else: #if self.errorEstimatorKind == "NONE":
                 err = self.getErrorEstimatorNone(mus)
         vbMng(self.trainedModel, "DEL", "Done evaluating error estimator", 10)
         if not return_max: return err
         if self.errorEstimatorKind[: 10] != "LOOK_AHEAD":
             idxMaxEst = np.empty(self.sampleBatchSize, dtype = int)
             errCP = copy(err)
             for j in range(self.sampleBatchSize):
                 k = np.argmax(errCP)
                 idxMaxEst[j] = k
                 if j + 1 < self.sampleBatchSize:
                     musZero = self.trainedModel.centerNormalize(mus, mus[k])
                     errCP *= np.linalg.norm(musZero.data, axis = 1)
         return err, idxMaxEst, err[idxMaxEst]
 
     def plotEstimator(self, *args, **kwargs):
         super().plotEstimator(*args, **kwargs)
         if self.errorEstimatorKind == "NONE":
             vbMng(self, "MAIN",
                   ("Warning! Error estimator has been arbitrarily normalized "
                    "before plotting."), 15)
 
     def greedyNextSample(self, *args,
                          **kwargs) -> Tuple[Np1D, int, float, paramVal]:
         """Compute next greedy snapshot of solution map."""
         RROMPyAssert(self._mode, message = "Cannot add greedy sample.")
         self.sampleBatchIdx += 1
         self.sampleBatchSize = totalDegreeN(self.npar - 1, self.sampleBatchIdx)
         err, muidx, maxErr, muNext =  super().greedyNextSample(*args, **kwargs)
         if maxErr is not None and (np.any(np.isnan(maxErr))
                                 or np.any(np.isinf(maxErr))):
             self.sampleBatchIdx -= 1
             self.sampleBatchSize = totalDegreeN(self.npar - 1,
                                                 self.sampleBatchIdx)
         if (self.errorEstimatorKind == "NONE" and not np.isnan(maxErr)
                                               and not np.isinf(maxErr)):
             maxErr = None
         return err, muidx, maxErr, muNext
 
     def _preliminaryTraining(self):
         """Initialize starting snapshots of solution map."""
         RROMPyAssert(self._mode, message = "Cannot start greedy algorithm.")
         if self.samplingEngine.nsamples > 0:
             return
         S = self.S
         self.sampleBatchIdx, self.sampleBatchSize, self._S = -1, 0, 0
         nextBatchSize = 1
         while self._S + nextBatchSize <= S:
             self.sampleBatchIdx += 1
             self.sampleBatchSize = nextBatchSize
             self._S += self.sampleBatchSize
             nextBatchSize = totalDegreeN(self.npar - 1,
                                          self.sampleBatchIdx + 1)
         super()._preliminaryTraining()
         self.M, self.N = ("AUTO",) * 2
 
     def setupApproxLocal(self) -> int:
         """Compute rational interpolant."""
         if self.checkComputedApprox(): return -1
         RROMPyAssert(self._mode, message = "Cannot setup approximant.")
         self.verbosity -= 10
         vbMng(self, "INIT", "Setting up local approximant.", 5)
         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.trainedModel.data.mus = copy(self.mus)
         self.catchInstability = 2
         unstable = False
         if self.E > 0:
             try:
                 Q = self._setupDenominator()[0]
             except RROMPyException as RE:
                 RROMPyWarning("Downgraded {}: {}".format(RE.__class__.__name__,
                                                          RE))
                 vbMng(self, "DEL", "", 7)
                 unstable = True
         else:
             Q = PI()
             Q.coeffs = np.ones((1,) * self.npar, dtype = np.complex)
             Q.npar = self.npar
             Q.polybasis = self.polybasis
         if not unstable:
             self.trainedModel.data.Q = copy(Q)
             try:
                 P = copy(self._setupNumerator())
             except RROMPyException as RE:
                 RROMPyWarning("Downgraded {}: {}".format(RE.__class__.__name__,
                                                          RE))
                 vbMng(self, "DEL", "", 7)
                 unstable = True
         if not unstable:
             self.trainedModel.data.P = copy(P)
             self.trainedModel.data.approxParameters = copy(
                                                          self.approxParameters)
         vbMng(self, "DEL", "Done setting up local approximant.", 5)
         self.catchInstability = 0
         self.verbosity += 10
         return 1 * unstable
         
     def setupApprox(self, plotEst : str = "NONE") -> int:
         val = super().setupApprox(plotEst)
         if val == 0:
             self._iterCorrector()
             self.trainedModel.data.approxParameters = copy(
                                                          self.approxParameters)
         return val
         
     def loadTrainedModel(self, filename:str):
         """Load trained reduced model from file."""
         super().loadTrainedModel(filename)
         self.sampleBatchIdx, self.sampleBatchSize, _S = -1, 0, 0
         nextBatchSize = 1
         while _S + nextBatchSize <= self.S + 1:
             self.sampleBatchIdx += 1
             self.sampleBatchSize = nextBatchSize
             _S += self.sampleBatchSize
             nextBatchSize = totalDegreeN(self.npar - 1,
                                          self.sampleBatchIdx + 1)
diff --git a/rrompy/reduction_methods/standard/greedy/reduced_basis_greedy.py b/rrompy/reduction_methods/standard/greedy/reduced_basis_greedy.py
index c239357..70b532d 100644
--- a/rrompy/reduction_methods/standard/greedy/reduced_basis_greedy.py
+++ b/rrompy/reduction_methods/standard/greedy/reduced_basis_greedy.py
@@ -1,161 +1,165 @@
 # Copyright (C) 2018 by the RROMPy authors
 #
 # This file is part of RROMPy.
 #
 # RROMPy is free software: you can redistribute it and/or modify
 # it under the terms of the GNU Lesser General Public License as published by
 # the Free Software Foundation, either version 3 of the License, or
 # (at your option) any later version.
 #
 # RROMPy is distributed in the hope that it will be useful,
 # but WITHOUT ANY WARRANTY; without even the implied warranty of
 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
 # GNU Lesser General Public License for more details.
 #
 # You should have received a copy of the GNU Lesser General Public License
 # along with RROMPy. If not, see <http://www.gnu.org/licenses/>.
 #
 
 from copy import deepcopy as copy
 from .generic_greedy_approximant import GenericGreedyApproximant
 from rrompy.reduction_methods.standard import ReducedBasis
 from rrompy.utilities.base import verbosityManager as vbMng
 from rrompy.utilities.numerical import dot
 from rrompy.utilities.exception_manager import RROMPyWarning, RROMPyAssert
 
 __all__ = ['ReducedBasisGreedy']
 
 class ReducedBasisGreedy(GenericGreedyApproximant, ReducedBasis):
     """
     ROM greedy RB approximant computation for parametric problems.
 
     Args:
         HFEngine: HF problem solver.
         mu0(optional): Default parameter. Defaults to 0.
         approxParameters(optional): Dictionary containing values for main
             parameters of approximant. Recognized keys are:
             - 'POD': whether to compute POD of snapshots; defaults to True;
+            - 'scaleFactorDer': scaling factors for derivative computation;
+                defaults to 'AUTO';
             - 'S': number of starting training points;
             - 'sampler': sample point generator;
             - 'greedyTol': uniform error tolerance for greedy algorithm;
                 defaults to 1e-2;
             - 'collinearityTol': collinearity tolerance for greedy algorithm;
                 defaults to 0.;
             - 'maxIter': maximum number of greedy steps; defaults to 1e2;
             - 'nTestPoints': number of test points; defaults to 5e2;
             - 'trainSetGenerator': training sample points generator; defaults
                 to sampler.
             Defaults to empty dict.
         approx_state(optional): Whether to approximate state. Defaults and must
             be True.
         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.
+            - 'POD': whether to compute POD of snapshots;
+            - 'scaleFactorDer': scaling factors for derivative computation;
             - 'greedyTol': uniform error tolerance for greedy algorithm;
             - 'collinearityTol': collinearity tolerance for greedy algorithm;
             - 'maxIter': maximum number of greedy steps;
             - 'nTestPoints': number of test points;
             - 'trainSetGenerator': training sample points generator.
         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 test points.
         sampler: Sample point generator.
         greedyTol: uniform error tolerance for greedy algorithm.
         collinearityTol: Collinearity tolerance for greedy algorithm.
         maxIter: maximum number of greedy steps.
         nTestPoints: number of starting training points.
         trainSetGenerator: training sample points generator.
         muBounds: list of bounds for parameter values.
         samplingEngine: Sampling engine.
         estimatorNormEngine: Engine for estimator norm computation.
         uHF: High fidelity solution(s) with parameter(s) lastSolvedHF as
             sampleList.
         lastSolvedHF: Parameter(s) corresponding to last computed high fidelity
             solution(s) as parameterList.
         uApproxReduced: Reduced approximate solution(s) with parameter(s)
             lastSolvedApprox as sampleList.
         lastSolvedApproxReduced: Parameter(s) corresponding to last computed
             reduced approximate solution(s) as parameterList.
         uApprox: Approximate solution(s) with parameter(s) lastSolvedApprox as
             sampleList.
         lastSolvedApprox: Parameter(s) corresponding to last computed
             approximate solution(s) as parameterList.
         As: List of sparse matrices (in CSC format) representing coefficients
             of linear system matrix.
         bs: List of numpy vectors representing coefficients of linear system
            RHS.
         ARBs: List of sparse matrices (in CSC format) representing coefficients
             of compressed linear system matrix.
         bRBs: List of numpy vectors representing coefficients of compressed 
             linear system RHS.
     """
     
     def __init__(self, *args, **kwargs):
         self._preInit()
         self._addParametersToList(toBeExcluded = ["R", "PODTolerance"])
         super().__init__(*args, **kwargs)
         self.R = "AUTO"
         self._postInit()
 
     def _setRAuto(self):
         self.R = self.S
 
     @property
     def PODTolerance(self):
         """Value of PODTolerance."""
         self._PODTolerance = -1
         return self._PODTolerance
     @PODTolerance.setter
     def PODTolerance(self, PODTolerance):
         RROMPyWarning(("PODTolerance is used just to simplify inheritance, "
                        "and its value cannot be changed from -1."))
         
     def setupApproxLocal(self) -> int:
         """Compute RB projection matrix."""
         if self.checkComputedApprox(): return -1
         RROMPyAssert(self._mode, message = "Cannot setup approximant.")
         self.verbosity -= 10
         vbMng(self, "INIT", "Setting up local approximant.", 5)
         vbMng(self, "INIT", "Computing projection matrix.", 7)
         pMat = self.samplingEngine.samples.data
         pMatEff = dot(self.HFEngine.C, 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]
             ARBs, bRBs = self.assembleReducedSystem(pMat)
             data.affinePoly = self.HFEngine.affinePoly
             self.HFEngine.buildA()
             self.HFEngine.buildb()
             data.thAs, data.thbs = self.HFEngine.thAs, self.HFEngine.thbs
             self.trainedModel.data = data
         else:
             self.trainedModel = self.trainedModel
             Sold = self.trainedModel.data.projMat.shape[1]
             ARBs, bRBs = self.assembleReducedSystem(pMat[:, Sold :],
                                                     pMat[:, : Sold])
             self.trainedModel.data.projMat = copy(pMatEff)
         self.trainedModel.data.mus = copy(self.mus)
         self.trainedModel.data.ARBs = ARBs
         self.trainedModel.data.bRBs = bRBs
         vbMng(self, "DEL", "Done computing projection matrix.", 7)
         self.trainedModel.data.approxParameters = copy(self.approxParameters)
         vbMng(self, "DEL", "Done setting up local approximant.", 5)
         self.verbosity += 10
         return 0
diff --git a/rrompy/reduction_methods/standard/nearest_neighbor.py b/rrompy/reduction_methods/standard/nearest_neighbor.py
index 4028a93..08ff30f 100644
--- a/rrompy/reduction_methods/standard/nearest_neighbor.py
+++ b/rrompy/reduction_methods/standard/nearest_neighbor.py
@@ -1,110 +1,114 @@
 # Copyright (C) 2018 by the RROMPy authors
 #
 # This file is part of RROMPy.
 #
 # RROMPy is free software: you can redistribute it and/or modify
 # it under the terms of the GNU Lesser General Public License as published by
 # the Free Software Foundation, either version 3 of the License, or
 # (at your option) any later version.
 #
 # RROMPy is distributed in the hope that it will be useful,
 # but WITHOUT ANY WARRANTY; without even the implied warranty of
 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
 # GNU Lesser General Public License for more details.
 #
 # You should have received a copy of the GNU Lesser General Public License
 # along with RROMPy. If not, see <http://www.gnu.org/licenses/>.
 #
 
 from copy import deepcopy as copy
 from .generic_standard_approximant import GenericStandardApproximant
 from rrompy.utilities.base import verbosityManager as vbMng
 from rrompy.utilities.numerical import dot
 from rrompy.utilities.exception_manager import RROMPyAssert
 
 __all__ = ['NearestNeighbor']
 
 class NearestNeighbor(GenericStandardApproximant):
     """
     ROM nearest neighbor approximant computation for parametric problems.
 
     Args:
         HFEngine: HF problem solver.
         mu0(optional): Default parameter. Defaults to 0.
         approxParameters(optional): Dictionary containing values for main
             parameters of approximant. Recognized keys are:
             - 'POD': whether to compute POD of snapshots; defaults to True;
+            - 'scaleFactorDer': scaling factors for derivative computation;
+                defaults to 'AUTO';
             - 'S': total number of samples current approximant relies upon;
             - 'sampler': sample point generator.
             Defaults to empty dict.
         approx_state(optional): Whether to approximate state. Defaults and must
             be True.
         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.
+            - 'POD': whether to compute POD of snapshots;
+            - 'scaleFactorDer': scaling factors for derivative computation.
         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.
         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.
     """
     
     @property
     def tModelType(self):
         from .trained_model.trained_model_nearest_neighbor import (
                                                    TrainedModelNearestNeighbor)
         return TrainedModelNearestNeighbor
 
     def setupApprox(self) -> int:
         """Compute RB projection matrix."""
         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)
         if self.POD:
             self.trainedModel.data.matR = copy(self.samplingEngine.RPOD)
         else:
             self.trainedModel.data.matR = 1.
         self.trainedModel.data.mus = copy(self.mus)
         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_interpolant.py b/rrompy/reduction_methods/standard/rational_interpolant.py
index ffb57f8..fdb5864 100644
--- a/rrompy/reduction_methods/standard/rational_interpolant.py
+++ b/rrompy/reduction_methods/standard/rational_interpolant.py
@@ -1,699 +1,707 @@
 # 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;
             - 'nNearestNeighbor': mumber of nearest neighbors considered in
                 numerator if polybasis allows; 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;
             - 'nNearestNeighbor': mumber of nearest neighbors considered in
                 numerator if polybasis allows;
             - '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.
         nNearestNeighbor: Number of nearest neighbors considered in numerator
             if polybasis allows.
         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",
                                    "nNearestNeighbor", "interpRcond",
                                    "robustTol", "correctorForce",
                                    "correctorTol", "correctorMaxIter"],
                                   ["MONOMIAL", "AUTO", "AUTO", "TOTAL", [1],
                                    -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: 
                 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 nNearestNeighbor(self):
         """Value of nNearestNeighbor."""
         return self._nNearestNeighbor
     @nNearestNeighbor.setter
     def nNearestNeighbor(self, nNearestNeighbor):
         self._nNearestNeighbor = nNearestNeighbor
         self._approxParameters["nNearestNeighbor"] = self.nNearestNeighbor
 
     @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) #SCALE
+                                   self._reorder, self._derIdxs,
+                                   self.scaleFactorRel)
         if self.POD:
             Qevaldiag = Qevaldiag.dot(self.samplingEngine.RPOD.T)
         if hasattr(self, "radialDirectionalWeights"):
             rDW = copy(self.radialDirectionalWeights)
         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 and (not hasattr(self, "radialDirectionalWeights")
             or self.radialDirectionalWeights[0] <= rDW[0] * 2 ** 6)):
             pParRest = [self.verbosity >= 5, self.polydegreetype == "TOTAL",
-                        {"derIdxs": self._derIdxs, "reorder": self._reorder}]
+                        {"derIdxs": self._derIdxs, "reorder": self._reorder,
+                         "scl": self.scaleFactorRel}]
             if self.polybasis in ppb:
                 p = PI()
             else:
                 pParRest = [self.radialDirectionalWeights] + pParRest
                 pParRest[-1]["nNearestNeighbor"] = self.nNearestNeighbor
                 p = RBI() if self.polybasis in rbpb else MLSI()
             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)
             if self.polybasis in ppb:
                 vbMng(self, "MAIN", ("Polyfit is poorly conditioned. Reducing "
                                      "M by 1."), 10)
                 self.M = self.M - 1
             else:
                 vbMng(self, "MAIN", ("Polyfit is poorly conditioned. "
                                      "Multiplying radialDirectionalWeights by "
                                      "2."), 10)
                 for j in range(self.npar):
                     self._radialDirectionalWeights[j] *= 2.
         if self.M < 0 or (hasattr(self, "radialDirectionalWeights")
                       and self.radialDirectionalWeights[0] > rDW[0] * 2 ** 6):
             raise RROMPyException(("Instability in computation of numerator. "
                                    "Aborting."))
         if self.polybasis in ppb:
             self.M = M
         else:
             self.radialDirectionalWeights = rDW
         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).T #SCALE
+                              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]
+        TEGenPar = [self.polybasis0, self._derIdxs, self._reorder,
+                    self.scaleFactorRel]
         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
index 3cd263a..19d4ffe 100644
--- a/rrompy/reduction_methods/standard/rational_moving_least_squares.py
+++ b/rrompy/reduction_methods/standard/rational_moving_least_squares.py
@@ -1,327 +1,337 @@
 # 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;
             - 'nNearestNeighbor': number of nearest neighbors considered in 
                 numerator if radialBasis allows; defaults to -1;
             - 'radialBasisDen': denominator radial basis type; defaults to
                 'GAUSSIAN';
             - 'radialDirectionalWeightsDen': radial basis weights for
                 interpolant denominator; defaults to 1;
             - 'nNearestNeighborDen': number of nearest neighbors considered in
                 denominator if radialBasisDen allows; 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;
             - 'nNearestNeighbor': number of nearest neighbors considered in
                 numerator if radialBasis allows;
             - 'radialBasisDen': denominator radial basis type;
             - 'radialDirectionalWeightsDen': radial basis weights for
                 interpolant denominator;
             - 'nNearestNeighborDen': number of nearest neighbors considered in
                 denominator if radialBasisDen allows;
             - '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.
         nNearestNeighbor: Number of nearest neighbors considered in numerator
             if radialBasis allows.
         radialBasisDen: Denominator radial basis type.
         radialDirectionalWeightsDen: Radial basis weights for interpolant
             denominator.
         nNearestNeighborDen: Number of nearest neighbors considered in
             denominator if radialBasisDen allows.
         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",
                                    "nNearestNeighborDen"],
                                   ["GAUSSIAN", "GAUSSIAN", 1, -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)
 
     @property
     def nNearestNeighborDen(self):
         """Value of nNearestNeighborDen."""
         return self._nNearestNeighborDen
     @nNearestNeighborDen.setter
     def nNearestNeighborDen(self, nNearestNeighborDen):
         self._nNearestNeighborDen = nNearestNeighborDen
         self._approxParameters["nNearestNeighborDen"] = (
                                                       self.nNearestNeighborDen)
 
     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]
-        if self.polydegreetype == "FULL": pPar[1] = [self.N] * self.npar
-        TN = pvP(*pPar) #SCALE
+                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]
-        if self.polydegreetype == "FULL": pPar[1] = [self.M] * self.npar
-        TM = pvP(*pPar) #SCALE
+                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.nNearestNeighbor = self.nNearestNeighbor
             data.radialBasisDen = self.radialBasisDen
             data.radialWeightsDen = self.radialDirectionalWeightsDen
             data.nNearestNeighborDen = self.nNearestNeighborDen
             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 391b8eb..b1d6b8d 100644
--- a/rrompy/reduction_methods/standard/rational_pade.py
+++ b/rrompy/reduction_methods/standard/rational_pade.py
@@ -1,325 +1,331 @@
 # 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 .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;
             - 'nNearestNeighbor': mumber of nearest neighbors considered in
                 numerator if polybasis allows; 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;
             - 'nNearestNeighbor': mumber of nearest neighbors considered in
                 numerator if polybasis allows;
             - '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.
         nNearestNeighbor: Number of nearest neighbors considered in numerator
             if polybasis allows.
         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) #SCALE
+                                   self._reorder, self._derIdxs,
+                                   self.scaleFactorRel)
         if self.POD:
             Qevaldiag = Qevaldiag.dot(self.samplingEngine.RPOD.T)
         if hasattr(self, "radialDirectionalWeights"):
             rDW = copy(self.radialDirectionalWeights)
         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 and (not hasattr(self, "radialDirectionalWeights")
             or self.radialDirectionalWeights[0] <= rDW[0] * 2 ** 6)):
             Seff = cfun(self.M, self.npar)
             pParRest = [self.verbosity >= 5, self.polydegreetype == "TOTAL",
                         {"derIdxs": [self._derIdxs[0][: Seff]],
-                         "reorder": self._reorder[: Seff]}]
+                         "reorder": self._reorder[: Seff],
+                         "scl": self.scaleFactorRel}]
             if self.polybasis in ppb:
                 p = PI()
             else:
                 pParRest = [self.radialDirectionalWeights] + pParRest
                 pParRest[-1]["nNearestNeighbor"] = self.nNearestNeighbor
                 p = RBI() if self.polybasis in rbpb else MLSI()
             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)
             if self.polybasis in ppb:
                 vbMng(self, "MAIN", ("Polyfit is poorly conditioned. Reducing "
                                      "M by 1."), 10)
                 self.M = self.M - 1
             else:
                 vbMng(self, "MAIN", ("Polyfit is poorly conditioned. "
                                      "Multiplying radialDirectionalWeights "
                                      "by 2."), 10)
                 for j in range(self.npar):
                     self._radialDirectionalWeights[j] *= 2.
         if self.M < 0 or (hasattr(self, "radialDirectionalWeights")
                       and self.radialDirectionalWeights[0] > rDW[0] * 2 ** 6):
             raise RROMPyException(("Instability in computation of numerator. "
                                    "Aborting."))
         if self.polybasis in ppb:
             self.M = M
         else:
             self.radialDirectionalWeights = rDW
         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]] #SCALE
+                        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/reduced_basis.py b/rrompy/reduction_methods/standard/reduced_basis.py
index d6abce7..5ef4b96 100644
--- a/rrompy/reduction_methods/standard/reduced_basis.py
+++ b/rrompy/reduction_methods/standard/reduced_basis.py
@@ -1,214 +1,218 @@
 # Copyright (C) 2018 by the RROMPy authors
 #
 # This file is part of RROMPy.
 #
 # RROMPy is free software: you can redistribute it and/or modify
 # it under the terms of the GNU Lesser General Public License as published by
 # the Free Software Foundation, either version 3 of the License, or
 # (at your option) any later version.
 #
 # RROMPy is distributed in the hope that it will be useful,
 # but WITHOUT ANY WARRANTY; without even the implied warranty of
 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
 # GNU Lesser General Public License for more details.
 #
 # You should have received a copy of the GNU Lesser General Public License
 # along with RROMPy. If not, see <http://www.gnu.org/licenses/>.
 #
 
 from copy import deepcopy as copy
 import numpy as np
 from .generic_standard_approximant import GenericStandardApproximant
 from rrompy.hfengines.base.linear_affine_engine import checkIfAffine
 from rrompy.reduction_methods.base.reduced_basis_utils import \
                                                      projectAffineDecomposition
 from rrompy.utilities.base.types import Np1D, Np2D, List, Tuple, sampList
 from rrompy.utilities.base import verbosityManager as vbMng
 from rrompy.utilities.numerical import dot
 from rrompy.utilities.exception_manager import (RROMPyWarning, RROMPyException,
                                                 RROMPyAssert)
 
 __all__ = ['ReducedBasis']
 
 class ReducedBasis(GenericStandardApproximant):
     """
     ROM RB approximant computation for parametric problems.
 
     Args:
         HFEngine: HF problem solver.
         mu0(optional): Default parameter. Defaults to 0.
         approxParameters(optional): Dictionary containing values for main
             parameters of approximant. Recognized keys are:
             - 'POD': whether to compute POD of snapshots; defaults to True;
+            - 'scaleFactorDer': scaling factors for derivative computation;
+                defaults to 'AUTO';
             - 'S': total number of samples current approximant relies upon;
             - 'sampler': sample point generator;
             - 'R': rank for Galerkin projection; defaults to 'AUTO', i.e.
                 maximum allowed;
             - 'PODTolerance': tolerance for snapshots POD; defaults to -1.
             Defaults to empty dict.
         approx_state(optional): Whether to approximate state. Defaults and must
             be True.
         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;
             - 'R': rank for Galerkin projection;
             - 'PODTolerance': tolerance for snapshots POD.
         parameterListCritical: Recognized keys of critical approximant
             parameters:
             - 'S': total number of samples current approximant relies upon;
             - 'sampler': sample point generator.
         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.
         R: Rank for Galerkin projection.
         muBounds: list of bounds for parameter values.
         samplingEngine: Sampling engine.
         uHF: High fidelity solution(s) with parameter(s) lastSolvedHF as
             sampleList.
         lastSolvedHF: Parameter(s) corresponding to last computed high fidelity
             solution(s) as parameterList.
         uApproxReduced: Reduced approximate solution(s) with parameter(s)
             lastSolvedApprox as sampleList.
         lastSolvedApproxReduced: Parameter(s) corresponding to last computed
             reduced approximate solution(s) as parameterList.
         uApprox: Approximate solution(s) with parameter(s) lastSolvedApprox as
             sampleList.
         lastSolvedApprox: Parameter(s) corresponding to last computed
             approximate solution(s) as parameterList.
     """
     
     def __init__(self, *args, **kwargs):
         self._preInit()
         self._addParametersToList(["R", "PODTolerance"], ["AUTO", -1])
         super().__init__(*args, **kwargs)
         checkIfAffine(self.HFEngine, "apply RB method")
         if not self.approx_state:
             raise RROMPyException("Must compute RB approximation of state.")
         self._postInit()
 
     @property
     def tModelType(self):
         from .trained_model.trained_model_reduced_basis import (
                                                       TrainedModelReducedBasis)
         return TrainedModelReducedBasis
 
     @property
     def R(self):
         """Value of R. Its assignment may change S."""
         return self._R
     @R.setter
     def R(self, R):
         if isinstance(R, str):
             R = R.strip().replace(" ","")
             if "-" not in R: R = R + "-0"
             self._R_isauto, self._R_shift = True, int(R.split("-")[-1])
             R = 0
         if R < 0: raise RROMPyException("R must be non-negative.")
         self._R = R
         self._approxParameters["R"] = self.R
 
     def _setRAuto(self):
         self.R = max(0, self.S - self._R_shift)
         vbMng(self, "MAIN", "Automatically setting R to {}.".format(self.R),
               25)
 
     @property
     def PODTolerance(self):
         """Value of PODTolerance."""
         return self._PODTolerance
     @PODTolerance.setter
     def PODTolerance(self, PODTolerance):
         self._PODTolerance = PODTolerance
         self._approxParameters["PODTolerance"] = self.PODTolerance
 
     def _setupProjectionMatrix(self):
         """Compute projection matrix."""
         RROMPyAssert(self._mode, message = "Cannot setup numerator.")
         vbMng(self, "INIT", "Starting computation of projection matrix.", 7)
         if hasattr(self, "_R_isauto"):
             self._setRAuto()
         else:
             if self.S < self.R:
                 RROMPyWarning(("R too large compared to S. Reducing R by "
                                "{}").format(self.R - self.S))
                 self.S = self.S
         try:
             if self.POD:
                 U, s, _ = np.linalg.svd(self.samplingEngine.RPOD)
                 s = s ** 2.
             else:
                 Gramian = self.HFEngine.innerProduct(
                                   self.samplingEngine.samples,
                                   self.samplingEngine.samples, is_state = True)
                 U, s, _ = np.linalg.svd(Gramian)
         except np.linalg.LinAlgError as e:
             raise RROMPyException(e)
         snorm = np.cumsum(s[::-1]) / np.sum(s)
         nPODTrunc = min(self.S - np.argmax(snorm > self.PODTolerance),
                         self.R)
         pMat = dot(self.samplingEngine.samples, U[:, : nPODTrunc])
         vbMng(self, "MAIN",
               ("Assembling {}x{} projection matrix from {} "
                "samples.").format(*(pMat.shape), self.S), 5)
         vbMng(self, "DEL", "Done computing projection matrix.", 7)
         return pMat
     
     def setupApprox(self) -> int:
         """Compute RB projection matrix."""
         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._setupProjectionMatrix().data
         pMatEff = dot(self.HFEngine.C, 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.affinePoly = self.HFEngine.affinePoly
             data.thAs, data.thbs = self.HFEngine.thAs, self.HFEngine.thbs
             self.trainedModel.data = data
         else:
             self.trainedModel = self.trainedModel
             self.trainedModel.data.projMat = copy(pMatEff)
         self.trainedModel.data.mus = copy(self.mus)
         ARBs, bRBs = self.assembleReducedSystem(pMat)
         self.trainedModel.data.ARBs = ARBs
         self.trainedModel.data.bRBs = bRBs
         self.trainedModel.data.approxParameters = copy(self.approxParameters)
         vbMng(self, "DEL", "Done setting up approximant.", 5)
         return 0
 
     def assembleReducedSystem(self, pMat : sampList = None,
                               pMatOld : sampList = None)\
                                               -> Tuple[List[Np2D], List[Np1D]]:
         """Build affine blocks of RB linear system through projections."""
         if pMat is None:
             self.setupApprox()
             ARBs = self.trainedModel.data.ARBs
             bRBs = self.trainedModel.data.bRBs
         else:
             self.HFEngine.buildA()
             self.HFEngine.buildb()
             vbMng(self, "INIT", "Projecting affine terms of HF model.", 10)
             ARBsOld = None if pMatOld is None else self.trainedModel.data.ARBs
             bRBsOld = None if pMatOld is None else self.trainedModel.data.bRBs
             ARBs, bRBs = projectAffineDecomposition(self.HFEngine.As,
                                                     self.HFEngine.bs, pMat,
                                                     ARBsOld, bRBsOld, pMatOld)
             vbMng(self, "DEL", "Done projecting affine terms.", 10)
         return ARBs, bRBs