diff --git a/rational_interpolation_method.pdf b/rational_interpolation_method.pdf
index a6088f0..9223c50 100644
Binary files a/rational_interpolation_method.pdf and b/rational_interpolation_method.pdf differ
diff --git a/rational_interpolation_method.tex b/rational_interpolation_method.tex
index 1182a8d..ce72be3 100644
--- a/rational_interpolation_method.tex
+++ b/rational_interpolation_method.tex
@@ -1,225 +1,190 @@
 \documentclass[10pt,a4paper]{article}
 \usepackage[left=1in,right=1in,top=1in,bottom=1in]{geometry}
 \usepackage{amsmath}
 \usepackage{amsfonts}
 \usepackage{amssymb}
 \usepackage{hyperref}
 \usepackage{xcolor}
 
 \setlength{\parindent}{0pt}
 \newcommand{\code}[1]{{\color{blue}\texttt{#1}}}
 \newcommand{\footpath}[1]{\footnote{\path{#1}}}
 \newcommand{\N}{\mathbb{N}}
 \newcommand{\R}{\mathbb{R}}
 \newcommand{\C}{\mathbb{C}}
 \newcommand{\I}{\mathcal{I}}
 \DeclareMathOperator*{\argmin}{arg\,min}
 \newcommand{\inner}[2]{\left\langle#1,#2\right\rangle_V}
 \newcommand{\norm}[1]{\left\|#1\right\|_V}
 
 \title{\bf The \code{RROMPy} rational interpolation method}
 \author{D. Pradovera, CSQI, EPF Lausanne -- \texttt{davide.pradovera@epfl.ch}}
 \date{}
 
 \hypersetup{ pdftitle={The RROMPy rational interpolation method}, pdfauthor={Davide Pradovera} }
 
 \begin{document}
 \maketitle
 
 \section*{Introduction}
 This document provides an explanation for the numerical method provided by the class \code{Rational Interpolant}\footpath{./rrompy/reduction_methods/standard/rational_interpolant.py} and daughters, e.g. \code{Rational Interpolant Greedy}\footpath{./rrompy/reduction_methods/standard/greedy/rational_interpolant_greedy.py}, as well as most of the pivoted approximants\footpath{./rrompy/reduction_methods/pivoted/{,greedy/}rational_interpolant_*.py}.
 
 We restrict the discussion to the single-parameter case, and most of the focus will be dedicated to the impact of the \code{functionalSolve} parameter, whose allowed values are
 \begin{itemize}
 \item \code{NORM} (default): see \ref{sec:norm}; allows for derivative information, i.e. repeated sample points.
 \item \code{DOMINANT}: see \ref{sec:dominant}; allows for derivative information, i.e. repeated sample points.
 \item \code{BARYCENTRIC\_NORM}: see \ref{sec:barycentricnorm}; does not allow for a Least Squares (LS) approach.
 \item \code{BARYCENTRIC\_AVERAGE}: see \ref{sec:barycentricaverage}; does not allow for a Least Squares (LS) approach.
-\item \code{NODAL}: see \ref{sec:nodal}; iterative method.
 \end{itemize}
 
 The main reference throughout the present document is \cite{Pradovera}.
 
 \section{Aim of approximation}
 We seek an approximation of $u:\C\to V$, with $(V,\inner{\cdot}{\cdot})$ a complex\footnote{The inner product is linear (resp. conjugate linear) in the first (resp. second) argument: $\inner{\alpha v}{\beta w}=\alpha\overline{\beta}\inner{v}{w}$.} Hilbert space (with induced norm $\norm{\cdot}$), of the form $\widehat{p}/\widehat{q}$, where $\widehat{p}:\C\to V$ and $\widehat{q}:\C\to\C$. For a given denominator $\widehat{q}$, the numerator $\widehat{p}$ is found by interpolation (possibly, LS or based on radial basis functions) of $\widehat{q}u$. Hence, here we focus on the computation of the denominator $\widehat{q}$.
 
 Other than the choice of target function $u$, the parameters which affect the computation of $\widehat{q}$ are:
 \begin{itemize}
 \item $\code{mus}\subset\C$ ($\{\mu_j\}_{j=1}^S$ below); for all \code{functionalSolve} values but \code{NORM} and \code{DOMINANT}, the $S$ points must be distinct.
-\item $\code{N}\in\N$ ($N$ below); for \code{BARYCENTRIC}, $N$ must equal $S-1$.
+\item $\code{N}\in\N$ ($N$ below); for \code{BARYCENTRIC\_*}, $N$ must equal $S-1$.
 \item $\code{polybasis0}\in\{\code{"CHEBYSHEV"}, \code{"LEGENDRE"}, \code{"MONOMIAL"}\}$; only for \code{NORM} and \code{DOMINANT}.
 \end{itemize}
 For simplicity, we will consider only the case of $S$ distinct sample points. One can deal with the case of confluent sample points by extending the standard (Lagrange) interpolation steps to Hermite-Lagrange ones.
 
 The main motivation behind the method involves the modified approximation problem
 \[u\approx\I^N\left(\Big(\big(\mu_j,\widehat{q}(\mu_j)u(\mu_j)\big)\Big)_{j=1}^S\right)\Big/\widehat{q},\]
 where $\widehat{q}:\C\to\C$ is a polynomial of degree $\leq N<S$, and $\I^N:(\C\times V)^S\to\mathbb{P}^N(\C;V)$ is a (LS) polynomial interpolation operator, which maps $S$ samples of a function (which lie in $V$) to a polynomial of degree $\leq N$ with coefficients in $V$.
 
 More precisely, let
 \[\mathbb{P}^N(\C;W)=\left\{\mu\mapsto\sum_{i=0}^N\alpha_i\mu^i\ :\ \alpha_0,\ldots,\alpha_N\in W\right\},\]
 with $W$ either $\C$ or $V$. We set
 \[\I^N\left(\big((\mu_j,\psi_j)\big)_{j=1}^S\right)\bigg|_\mu=\argmin_{p\in\mathbb{P}^N(\C;V)}\sum_{j=1}^S\norm{p(\mu_j)-\psi_j}^2.\]
 In \code{RROMPy}, we compute (LS-)interpolants by employing normal equations: given a basis $\{\phi_i\}_{i=0}^N$ of $\mathbb{P}^N(\C;\C)$, we expand
 \[\I^N\left(\big((\mu_j,\psi_j)\big)_{j=1}^S\right)=\sum_{i=0}^Nc_i\phi_i\]
 and observe that, for optimality, the coefficients $\{c_i\}_{i=0}^N\subset V$ must satisfy
 \[\sum_{j=1}^S\sum_{l=0}^N\overline{\phi_i(\mu_j)}\phi_l(\mu_j)c_l=\sum_{j=1}^S\overline{\phi_i(\mu_j)}\psi_j\quad\forall i=0,\ldots,N,\]
 i.e., in matrix form\footnote{The superscript ${}^H$ denotes adjunction (conjugate transposition), i.e. $A^H=\overline{A}^\top$.},
 \[\underbrace{[c_i]_{i=0}^N}_{\in V^{N+1}}=\bigg(\underbrace{\Big[\overline{\phi_i(\mu_j)}\Big]_{i=0,j=1}^{N,S}}_{=:\Phi^H\in\C^{(N+1)\times S}}\underbrace{\Big[\phi_i(\mu_j)\Big]_{j=1,i=0}^{S,N}}_{=:\Phi\in\C^{S\times(N+1)}}\bigg)^{-1}\underbrace{\Big[\overline{\phi_i(\mu_j)}\Big]_{i=0,j=1}^{N,S}}_{=:\Phi^H\in\C^{(N+1)\times S}}\underbrace{[\psi_j]_{j=1}^S}_{\in V^{S}}.\]
 In practice the polynomial basis $\{\phi_i\}_{i=0}^N$ is determined by the value of $\code{polybasis0}$:
 \begin{itemize}
 \item If $\code{polybasis0}=\code{"CHEBYSHEV"}$, then $\phi_k(\mu)=\mu^k$ for $k\in\{0,1\}$ and $\phi_k(\mu)=2\mu\phi_{k-1}(\mu)-\phi_{k-2}(\mu)$ for $k\geq2$.
 \item If $\code{polybasis0}=\code{"LEGENDRE"}$, then $\phi_k(\mu)=\mu^k$ for $k\in\{0,1\}$ and $\phi_k(\mu)=(2-1/k)\mu\phi_{k-1}(\mu)-(1-1/k)\phi_{k-2}(\mu)$ for $k\geq2$.
 \item If $\code{polybasis0}=\code{"MONOMIAL"}$, then $\phi_k(\mu)=\mu^k$ for $k\geq0$.
 \end{itemize}
 
 The actual denominator $\widehat{q}$ is found as
 \begin{equation}\label{eq:glob1}
 \widehat{q}=\argmin_{\substack{q\in\mathbb{P}^N(\C;\C)\\(\star)}}\norm{\frac{\textrm{d}^N}{\textrm{d}\mu^N}\I^N\left(\Big(\big(\mu_j,q(\mu_j)u(\mu_j)\big)\Big)_{j=1}^S\right)}
 \end{equation}
 where $(\star)$ is a normalization condition (which changes depending on \code{functionalSolve}) to exclude the trivial minimizer $\widehat{q}=0$. 
 
-Broadly speaking, the five methods described differ in terms of the constraint $(\star)$, as well as of the degrees of freedom which are chosen to represent the denominator $q$.
+Broadly speaking, the methods described differ in terms of the constraint $(\star)$, as well as of the degrees of freedom which are chosen to represent the denominator $q$.
 
 \section{Polynomial coefficients as degrees of freedom}
 
 If the polynomial basis $\{\phi_i\}_{i=0}^N$ is hierarchical (as the three ones above), then the $N$-th derivative of $\I^N$ coincides with the coefficient $c_N$, and we have
 \begin{equation}\label{eq:poly1}
 \widehat{q}=\argmin_{\substack{q\in\mathbb{P}^N(\C;\C)\\(\star)}}\norm{[\underbrace{0,\ldots,0}_{N},1]\big(\Phi^H\Phi\big)^{-1}\Phi^H[q(\mu_j)u(\mu_j)]_{j=1}^S}.
 \end{equation}
 
 Using the Kronecker delta ($\delta_{ij}=1$ if $i=j$ and $\delta_{ij}=0$ if $i\neq j$), the last term $[q(\mu_j)u(\mu_j)]_{j=1}^S\in V^S$ can be factored into
 \begin{equation}\label{eq:poly2}
 \Big[u(\mu_j)\delta_{jj'}\Big]_{j=1,j'=1}^{S,S}[q(\mu_j)]_{j=1}^S=\Big[u(\mu_j)\delta_{jj'}\Big]_{j=1,j'=1}^{S,S}\Phi[q_i]_{i=0}^N,
 \end{equation}
 where we have expanded the polynomial $q$ using the basis\footnote{In theory, nothing prevents us from using different bases for $\I^N$ and $q$, cf. \ref{sec:observations}.} $\{\phi_i\}_{i=0}^N$: $q(\mu)=\sum_{i=0}^Nq_i\phi_i(\mu)$, with coefficients $\{q_i\}_{i=0}^N\subset\C$.
 
 Combining \eqref{eq:poly1} and \eqref{eq:poly2}, it is useful to consider the $(N+1)\times(N+1)$ Hermitian matrix with entries ($0\leq i,i'\leq N$)
 \begin{equation}\label{eq:poly3}
 G_{ii'}=\inner{\sum_{j'=1}^S\left(\big(\Phi^H\Phi\big)^{-1}\Phi^H\right)_{Nj'}(\Phi)_{j'i'}u(\mu_{j'})}{\sum_{j=1}^S\left(\big(\Phi^H\Phi\big)^{-1}\Phi^H\right)_{Nj}(\Phi)_{ji}u(\mu_j)}.
 \end{equation}
 If $(\star)$ is quadratic (resp. linear) in $[q_i]_{i=0}^N$, then we can cast the computation of the denominator as a quadratically (resp. linearly) constrained quadratic program involving $G$.
 
 \subsection{Quadratic constraint}\label{sec:norm}
 We constrain $[\widehat{q}_i]_{i=0}^N$ to have unit (Euclidean) norm. The resulting optimization problem can be cast as a minimal (normalized) eigenvector problem for $G$ in \eqref{eq:poly3}. More explicitly,
 \[[\widehat{q}_i]_{i=0}^N=\argmin_{\substack{\mathbf{q}\in\C^{N+1}\\\|\mathbf{q}\|_2=1}}\mathbf{q}^HG\mathbf{q}.\]
 
 \subsection{Linear constraint}\label{sec:dominant}
 We constrain $\widehat{q}_N=1$, thus forcing $q$ to be monic, with degree exactly $N$. Given $G$ in \eqref{eq:poly3}, the resulting optimization problem can be solved rather easily as:
 \[[\widehat{q}_i]_{i=0}^N=\frac{G^{-1}\mathbf{e}_{N+1}}{\mathbf{e}_{N+1}^\top G^{-1}\mathbf{e}_{N+1}},\quad\text{with }\mathbf{e}_{N+1}=[0,\ldots,0,1]^\top\in\C^{N+1}.\]
 
 
 \section{Barycentric coefficients as degrees of freedom}\label{sec:bary}
 Here we assume that the sample points are distinct, and such that $N=S-1$. We can choose for convenience a non-hierarchical basis, dependent on the sample points, for $q$ and $\I^N$, taking inspiration from barycentric interpolation:
 \begin{equation}\label{eq:bary1}
 \phi_i(\mu)=\prod_{\substack{j=1\\j\neq i+1}}^S(\mu-\mu_j).
 \end{equation}
 Since all elements of the basis are monic and of degree exactly $N$, the minimization problem can be cast as
 \begin{equation}\label{eq:bary2}
-\widehat{q}=\argmin_{\substack{q\in\mathbb{P}^N(\C;\C)\\(\star)}}\norm{[\underbrace{1,\ldots,1}_{N+1}]\big(\Phi^H\Phi\big)^{-1}\Phi^H\Big[u(\mu_j)\delta_{jj'}\Big]_{j=1,j'=1}^{S,S}\Phi[q_i]_{i=0}^N}.
+\widehat{\mathbf{q}}=\argmin_{\substack{\mathbf{q}\in\C^{N+1}\\(\star)}}\norm{[\underbrace{1,\ldots,1}_{N+1}]\big(\Phi^H\Phi\big)^{-1}\Phi^H\Big[u(\mu_j)\delta_{jj'}\Big]_{j=1,j'=1}^{S,S}\Phi\mathbf{q}\,}.
 \end{equation}
 At the same time, it is easy to see from \eqref{eq:bary1} that the Vandermonde-like matrix $\Phi$ is diagonal, so that
 \[\big(\Phi^H\Phi\big)^{-1}\Phi^H\Big[u(\mu_j)\delta_{jj'}\Big]_{j=1,j'=1}^{S,S}\Phi=\big(\Phi^H\Phi\big)^{-1}\Phi^H\Phi\Big[u(\mu_j)\delta_{jj'}\Big]_{j=1,j'=1}^{S,S}=\Big[u(\mu_j)\delta_{jj'}\Big]_{j=1,j'=1}^{S,S},\]
 and
 \begin{equation}\label{eq:bary3}
-\widehat{q}=\argmin_{\substack{\mathbf{q}\in\C^{N+1}\\(\star)}}\norm{\sum_{i=0}^Nu(\mu_{i+1})q_i}.
+\widehat{\mathbf{q}}=\argmin_{\substack{\mathbf{q}\in\C^{N+1}\\(\star)}}\norm{\sum_{i=0}^Nu(\mu_{i+1})q_i}.
 \end{equation}
 
 Considering \eqref{eq:bary3}, it is useful to define the $(N+1)\times(N+1)$ Hermitian (``snapshot Gramian'') matrix with entries ($0\leq i,i'\leq N$)
 \begin{equation}\label{eq:bary4}
 G_{ii'}=\inner{u(\mu_{i'+1})}{u(\mu_{i+1})}.
 \end{equation}
 So, once again, if $(\star)$ is quadratic (resp. linear) in $[q_i]_{i=0}^N$, then we can cast the computation of the denominator as a quadratically (resp. linearly) constrained quadratic program involving $G$.
 
-Before specifying the kind of normalization enforced, it is important to make a remark on numerical stability. The basis in \eqref{eq:bary1} is actually just a ($i$-dependent) factor away from being the Lagrangian one (for which $\phi_i(\mu_j)$ would equal $\delta_{(i+1)j}$ instead of $\delta_{(i+1)j}\prod_{k\neq i+1}(\mu_j-\mu_k)$ as it does in our case). As such, it is generally a bad idea to numerically evaluate $q$ starting from its expansion coefficients with respect to $\{\phi_i\}_{i=0}^N$. We get around this by exploiting the following trick, whose foundation is in \cite[Section 2.3.3]{Klein}: the roots of $\widehat{q}=\sum_{i=0}^N\widehat{q}_i\phi_i$ are the $N$ finite eigenvalues $\lambda$ of the generalized $(N+2)\times(N+2)$ eigenproblem
+Before specifying the kind of normalization enforced, it is important to make a remark on numerical stability. The basis in \eqref{eq:bary1} is actually just a ($i$-dependent) factor away from being the Lagrangian one (for which $\phi_i(\mu_j)$ would equal $\delta_{(i+1)j}$ instead of
+\[\delta_{(i+1)j}\prod_{k\neq i+1}(\mu_j-\mu_k),\]
+as it does in our case). As such, it is generally a bad idea to numerically evaluate $q$ starting from its expansion coefficients with respect to $\{\phi_i\}_{i=0}^N$. We get around this by exploiting the following trick, whose foundation is in \cite[Section 2.3.3]{Klein}: the roots of $\widehat{q}=\sum_{i=0}^N\widehat{q}_i\phi_i$ are the $N$ finite eigenvalues $\lambda$ of the generalized $(N+2)\times(N+2)$ eigenproblem
 \begin{equation}
 \textrm{Det}\left(\begin{bmatrix}
 0 & \widehat{q}_0 & \cdots & \widehat{q}_N\\
 1 & \mu_1 &  & \\
 \vdots & & \ddots & \\
 1 & & & \mu_S
 \end{bmatrix}-\begin{bmatrix}
 0 & & & \\
  & 1 &  & \\
  & & \ddots & \\
  & & & 1
 \end{bmatrix}\lambda\right)=0.
 \end{equation}
-This computation is numerically more stable than most manipulations of a polynomial in the basis \eqref{eq:bary1}. Unfortunately, the ensuing computation of $\widehat{p}$ invariably suffers from numerical instabilities as a consequence of this.
+This computation is numerically more stable than most other manipulations of a polynomial in the basis \eqref{eq:bary1}.
 
-Once the roots of $\widehat{q}$ have been computed, one can either convert it to nodal form \eqref{eq:pole1} or forgo using $\widehat{q}$ completely, in favor of a Heaviside-like approximation involving the newly computed roots $\{\widehat{\lambda}_i\}_{i=1}^N$ as poles:
+Once the roots of $\widehat{q}$ have been computed, one can either convert it to nodal form
+\begin{equation}\label{eq:pole1}
+\widehat{q}(\mu)\propto\prod_{i=1}^N(\mu-\widehat{\lambda}_i),
+\end{equation}
+or forgo using $\widehat{q}$ completely, in favor of a Heaviside-like approximation involving the newly computed roots $\{\widehat{\lambda}_i\}_{i=1}^N$ as poles:
 \[
 \frac{\widehat{p}(\mu)}{\widehat{q}(\mu)}\quad\rightsquigarrow\quad\sum_{i=1}^N\frac{\widehat{r}_i}{\mu-\widehat{\lambda}_i}+\widetilde{p}(\mu),
 \]
 with $\widetilde{p}$, e.g., a polynomial (of degree at most $S-N-1$) or a combination of radial basis functions. See the final paragraph in \ref{sec:observations} for a slightly more detailed motivation of why the Heaviside form of the approximant might be more useful than the standard rational one $\widehat{p}/\widehat{q}$ in practice.
 
 \subsection{Quadratic constraint}\label{sec:barycentricnorm}
 We constrain $[\widehat{q}_i]_{i=0}^N$ to have unit (Euclidean) norm. The resulting optimization problem can be cast as a minimal (normalized) eigenvector problem for $G$ in \eqref{eq:bary4}. More explicitly,
 \[[\widehat{q}_i]_{i=0}^N=\argmin_{\substack{\mathbf{q}\in\C^{N+1}\\\|\mathbf{q}\|_2=1}}\mathbf{q}^HG\mathbf{q}.\]
 
 \subsection{Linear constraint}\label{sec:barycentricaverage}
 We constrain $\sum_{i=0}^N\widehat{q}_i=1$, so that the polynomial $\widehat{q}$ is monic. Given $G$ in \eqref{eq:bary4}, the resulting optimization problem can be solved rather easily as:
 \[[\widehat{q}_i]_{i=0}^N=\frac{G^{-1}\mathbf{1}_{N+1}}{\mathbf{1}_{N+1}^\top G^{-1}\mathbf{1}_{N+1}},\quad\text{with }\mathbf{1}_{N+1}=[1,\ldots,1]^\top\in\C^{N+1}.\]
 
 
-\section{Poles as degrees of freedom}\label{sec:nodal}
-Here we assume that the sample points are distinct. One can avoid expressing $q$ explicitly and work directly with its roots by employing a nodal form
-\begin{equation}\label{eq:pole1}
-q(\mu)=\prod_{i=1}^N(\mu-\lambda_i),
-\end{equation}
-with $\{\lambda_i\}_{i=1}^N\subset\C$. (It is important to note that, as in \ref{sec:bary}, we are constraining $q$ to have degree exactly $N$.)
-
-The optimization problem that allows to find the desired poles can now be cast as
-\[\{\widehat{\lambda}_i\}_{i=1}^N=\argmin_{\{\lambda_i\}_{i=1}^N\subset\C}\norm{\frac{\textrm{d}^N}{\textrm{d}\mu^N}\I^N\left(\Big(\big(\mu_j,u(\mu_j)\prod_{i=1}^N(\mu_j-\lambda_i)\big)\Big)_{j=1}^S\right)},\]
-which, if $\I^N$ is expressed via a hierarchical basis, is equivalent to
-\begin{align}
-\{\widehat{\lambda}_i\}_{i=1}^N=&\argmin_{\{\lambda_i\}_{i=1}^N\subset\C}\norm{[\underbrace{0,\ldots,0}_{N},1]\big(\Phi^H\Phi\big)^{-1}\Phi^H\Big[u(\mu_j)\delta_{jj'}\Big]_{j=1,j'=1}^{S,S}\left[\prod_{i=1}^N(\mu_j-\lambda_i)\right]_{j=1}^S}^2\label{eq:pole2}\\
-=&\argmin_{\{\lambda_i\}_{i=1}^N\subset\C}\norm{\mathbf{y}^H\left[q(\mu_j)\right]_{j=1}^S}^2\nonumber
-\end{align}
-(the outer square has been added for later convenience), with
-\begin{equation}\label{eq:poleadd}
-\mathbf{y}=\Big[\overline{u(\mu_j)}\delta_{jj'}\Big]_{j=1,j'=1}^{S,S}\Phi\big(\Phi^H\Phi\big)^{-1}[\underbrace{0,\ldots,0}_{N},1]^\top\in V^S.
-\end{equation}
-This problem is nonconvex and highly nonlinear with respect to $\{\lambda_i\}_{i=1}^N\subset\C$, and an iterative method is necessary for its solution.
-
-Two observations aid us here:
-\begin{itemize}
-\item The target functional has a relatively simple structure, allowing to compute exactly its Jacobian and Hessian with respect to the poles. Explicitly, for all $\{\lambda_i\}_{i=1}^N\subset\C$ and $k=1,\ldots,N$, the partial derivative with respect to $\lambda_k$ at $\{\lambda_i\}_{i=1}^N\subset\C$ equals
-\begin{equation}\label{eq:pole3}
--2\inner{\mathbf{y}^H\left[q(\mu_j)\right]_{j=1}^S}{\mathbf{y}^H\left[\frac{q(\mu_j)}{\mu_j-\lambda_k}\right]_{j=1}^S},
-\end{equation}
-where $\mathbf{y}$ from \eqref{eq:poleadd} does not change as the iterations proceed. Similarly, for all $\{\lambda_i\}_{i=1}^N\subset\C$ and $k,k'=1,\ldots,N$, the second order partial derivative with respect to $\lambda_k$ and $\lambda_{k'}$ at $\{\lambda_i\}_{i=1}^N\subset\C$ is
-\begin{multline}\label{eq:pole4}
-2\inner{\mathbf{y}^H\left[\frac{q(\mu_j)}{\mu_j-\lambda_{k'}}\right]_{j=1}^S}{\mathbf{y}^H\left[\frac{q(\mu_j)}{\mu_j-\lambda_k}\right]_{j=1}^S}+\\
-+2(1-\delta_{kk'})\inner{\mathbf{y}^H\left[q(\mu_j)\right]_{j=1}^S}{\mathbf{y}^H\left[\frac{q(\mu_j)}{(\mu_j-\lambda_k)(\mu_j-\lambda_{k'})}\right]_{j=1}^S}
-\end{multline}
-Hence, we can use Newton's method.
-\item Due to the iterative nature of the solution strategy, an initial guess of the optimal poles is necessary. Luckily, we can employ any of the other four methods to this aim. In \code{RROMPy}, we apply \code{DOMINANT} as a preliminary step, compute the roots of the resulting $\widehat{q}$, and use them as initial guess for the Newton iterations.
-\end{itemize}
-
-As stopping criterion for Newton's method, we use the relative norm of the increment:
-\[\sum_{i=1}^N\left|\lambda_i^{(\textrm{iter}+1)}-\lambda_i^{(\textrm{iter})}\right|^2\leq\underbrace{\textrm{tolerance}}_{10^{-10}}\sum_{i=1}^N\left|\lambda_i^{(\textrm{iter}+1)}\right|^2.\]
-Since the initial guess is quite good, usually just a few (most often, 1) Newton iterations are necessary.
-
-
 \section{Minor observations}\label{sec:observations}
 \begin{itemize}
 \item For \code{BARYCENTRYC\_*}, a specific choice of polynomial basis for $\I^N$ was used to diagonalize the functional. Under the assumptions that the sample points are distinct and that $N=S-1$, one can employ the quasi-Lagrangian basis \eqref{eq:bary1} to expand $\I^N$ in the other approaches as well, thus simplifying significantly the structure of \eqref{eq:glob1}:
 \begin{equation*}
 \widehat{q}=\argmin_{\substack{q\in\mathbb{P}^N(\C;\C)\\(\star)}}\norm{\sum_{j=1}^Sq(\mu_j)u(\mu_j)\prod_{\substack{j'=1\\j'\neq j}}^S\frac1{\mu_j-\mu_{j'}}}.
 \end{equation*}
-Numerically, this has repercussions on the computation of the term $\big(\Phi^H\Phi\big)^{-1}\Phi^H$ in \eqref{eq:poly3} and of its adjoint $\Phi\big(\Phi^H\Phi\big)^{-1}$ in \eqref{eq:poleadd}.
+Numerically, this has repercussions on the computation of the term $\big(\Phi^H\Phi\big)^{-1}\Phi^H$ in \eqref{eq:poly3}.
 \item In general, \code{NORM} and \code{BARYCENTRIC\_NORM} can be expected to be more numerically stable than \code{DOMINANT} and \code{BARYCENTRIC\_AVERAGE}, respectively. This is due to the fact that the normalization is enforced in a more numerically robust fashion.
 \item If the snapshots are orthonormalized via \code{POD}\footpath{./rrompy/sampling/engines/sampling_engine_pod.py}, all the $V$-inner products (resp. norms) are recast as Euclidean inner products (resp. norms) involving the $R$ factor of the generalized ($V$-orthonormal) QR decomposition of the snapshots.
-\item If a univariate rational surrogate is built in the scope of multivariate pole-matching-based pivoted approximation\footpath{./rrompy/reduction_methods/pivoted/*}, the rational approximant is converted into a Heaviside/nodal representation when different surrogates are combined. As such, the \code{BARYCENTRYC\_*} or \code{NODAL} approaches may be preferable to avoid extra computations, as well as additional round-off artifacts.
+\item If a univariate rational surrogate is built in the scope of multivariate pole-matching-based pivoted approximation\footpath{./rrompy/reduction_methods/pivoted/*}, the rational approximant is converted into a Heaviside/nodal representation when different surrogates are combined. As such, the \code{BARYCENTRIC\_*} approach may be preferable to avoid extra computations, as well as additional round-off artifacts.
 \end{itemize}
 
 \begin{thebibliography}{1}
 
 \bibitem{Pradovera}% 
  D.~Pradovera, Interpolatory rational model order reduction of parametric problems lacking uniform inf-sup stability, SIAM J. Numer. Anal. 58 (2020) 2265--2293. doi:10.1137/19M1269695.
 
 \bibitem{Klein}%
  G.~Klein, Applications of Linear Barycentric Rational Interpolation, PhD Thesis no. 1762, Universit\'e de Fribourg (2012).
 \end{thebibliography}
 
-\end{document}
\ No newline at end of file
+\end{document}
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 c3035a9..1659aa5 100644
--- a/rrompy/reduction_methods/pivoted/greedy/generic_pivoted_greedy_approximant.py
+++ b/rrompy/reduction_methods/pivoted/greedy/generic_pivoted_greedy_approximant.py
@@ -1,650 +1,651 @@
 # Copyright (C) 2018-2020 by the RROMPy authors
 #
 # This file is part of RROMPy.
 #
 # RROMPy is free software: you can redistribute it and/or modify
 # it under the terms of the GNU Lesser General Public License as published by
 # the Free Software Foundation, either version 3 of the License, or
 # (at your option) any later version.
 #
 # RROMPy is distributed in the hope that it will be useful,
 # but WITHOUT ANY WARRANTY; without even the implied warranty of
 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
 # GNU Lesser General Public License for more details.
 #
 # You should have received a copy of the GNU Lesser General Public License
 # along with RROMPy. If not, see <http://www.gnu.org/licenses/>.
 #
 
 from abc import abstractmethod
 from copy import deepcopy as copy
 import numpy as np
 from collections.abc import Iterable
 from matplotlib import pyplot as plt
 from rrompy.reduction_methods.pivoted.generic_pivoted_approximant import (
                                             GenericPivotedApproximantBase,
                                             GenericPivotedApproximantPoleMatch)
 from rrompy.reduction_methods.pivoted.gather_pivoted_approximant import (
                                                       gatherPivotedApproximant)
 from rrompy.utilities.base.types import (Np1D, Np2D, Tuple, List, paramVal,
                                          paramList, ListAny)
 from rrompy.utilities.base import verbosityManager as vbMng
 from rrompy.utilities.numerical import dot
 from rrompy.utilities.numerical.point_matching import pointMatching
 from rrompy.utilities.numerical.point_distances import doubleDistanceMatrix
 from rrompy.utilities.exception_manager import (RROMPyException, RROMPyAssert,
                                                 RROMPyWarning)
 from rrompy.parameter import emptyParameterList
 from rrompy.utilities.parallel import (masterCore, indicesScatter,
                                        arrayGatherv, isend)
 
 __all__ = ['GenericPivotedGreedyApproximantPoleMatch']
 
 class GenericPivotedGreedyApproximantBase(GenericPivotedApproximantBase):
     _allowedEstimatorKindsMarginal = ["LOOK_AHEAD", "LOOK_AHEAD_RECOVER",
                                       "NONE"]
 
     def __init__(self, *args, **kwargs):
         self._preInit()
         self._addParametersToList(["matchingWeightError",
                                    "errorEstimatorKindMarginal",
                                    "greedyTolMarginal", "maxIterMarginal"],
                                   [0., "NONE", 1e-1, 1e2])
         super().__init__(*args, **kwargs)
         self._postInit()
 
     @property
     def scaleFactorDer(self):
         """Value of scaleFactorDer."""
         if self._scaleFactorDer == "NONE": return 1.
         if self._scaleFactorDer == "AUTO": return self._scaleFactorOldPivot
         return self._scaleFactorDer
     @scaleFactorDer.setter
     def scaleFactorDer(self, scaleFactorDer):
         if isinstance(scaleFactorDer, (str,)):
             scaleFactorDer = scaleFactorDer.upper()
         elif isinstance(scaleFactorDer, Iterable):
             scaleFactorDer = list(scaleFactorDer)
         self._scaleFactorDer = scaleFactorDer
         self._approxParameters["scaleFactorDer"] = self._scaleFactorDer
 
     @property
     def samplerMarginal(self):
         """Value of samplerMarginal."""
         return self._samplerMarginal
     @samplerMarginal.setter
     def samplerMarginal(self, samplerMarginal):
         if 'refine' not in dir(samplerMarginal):
             raise RROMPyException("Marginal sampler type not recognized.")
         GenericPivotedApproximantBase.samplerMarginal.fset(self,
                                                            samplerMarginal)
 
     @property
     def errorEstimatorKindMarginal(self):
         """Value of errorEstimatorKindMarginal."""
         return self._errorEstimatorKindMarginal
     @errorEstimatorKindMarginal.setter
     def errorEstimatorKindMarginal(self, errorEstimatorKindMarginal):
         errorEstimatorKindMarginal = errorEstimatorKindMarginal.upper()
         if errorEstimatorKindMarginal not in (
                                           self._allowedEstimatorKindsMarginal):
             RROMPyWarning(("Marginal error estimator kind not recognized. "
                            "Overriding to 'NONE'."))
             errorEstimatorKindMarginal = "NONE"
         self._errorEstimatorKindMarginal = errorEstimatorKindMarginal
         self._approxParameters["errorEstimatorKindMarginal"] = (
                                                self.errorEstimatorKindMarginal)
 
     @property
     def matchingWeightError(self):
         """Value of matchingWeightError."""
         return self._matchingWeightError
     @matchingWeightError.setter
     def matchingWeightError(self, matchingWeightError):
         self._matchingWeightError = matchingWeightError
         self._approxParameters["matchingWeightError"] = (
                                                       self.matchingWeightError)
 
     @property
     def greedyTolMarginal(self):
         """Value of greedyTolMarginal."""
         return self._greedyTolMarginal
     @greedyTolMarginal.setter
     def greedyTolMarginal(self, greedyTolMarginal):
         if greedyTolMarginal < 0:
             raise RROMPyException("greedyTolMarginal must be non-negative.")
         if (hasattr(self, "_greedyTolMarginal")
         and self.greedyTolMarginal is not None):
             greedyTolMarginalold = self.greedyTolMarginal
         else:
             greedyTolMarginalold = -1
         self._greedyTolMarginal = greedyTolMarginal
         self._approxParameters["greedyTolMarginal"] = self.greedyTolMarginal
         if greedyTolMarginalold != self.greedyTolMarginal:
             self.resetSamples()
 
     @property
     def maxIterMarginal(self):
         """Value of maxIterMarginal."""
         return self._maxIterMarginal
     @maxIterMarginal.setter
     def maxIterMarginal(self, maxIterMarginal):
         if maxIterMarginal <= 0:
             raise RROMPyException("maxIterMarginal must be positive.")
         if (hasattr(self, "_maxIterMarginal")
         and self.maxIterMarginal is not None):
             maxIterMarginalold = self.maxIterMarginal
         else:
             maxIterMarginalold = -1
         self._maxIterMarginal = maxIterMarginal
         self._approxParameters["maxIterMarginal"] = self.maxIterMarginal
         if maxIterMarginalold != self.maxIterMarginal:
             self.resetSamples()
 
     def resetSamples(self):
         """Reset samples."""
         super().resetSamples()
         if not hasattr(self, "_temporaryPivot"):
             self._mus = emptyParameterList()
             self._musMarginal = emptyParameterList()
             if hasattr(self, "samplerMarginal"): self.samplerMarginal.reset()
         if hasattr(self, "samplingEngine") and self.samplingEngine is not None:
             self.samplingEngine.resetHistory()
 
     def _getDistanceApp(self, polesEx:Np1D, resEx:Np2D,
                         muTest:paramVal) -> float:
         polesAp = self.trainedModel.interpolateMarginalPoles(muTest)[0]
         if self.matchingWeightError != 0:
             resAp = self.trainedModel.interpolateMarginalCoeffs(muTest)[0][
                                                              : len(polesAp), :]
             resEx = dot(self.trainedModel.data.projMat, resEx)
             resAp = dot(self.trainedModel.data.projMat, resAp)
         else:
             resAp = None
         dist = doubleDistanceMatrix(polesEx, polesAp, self.matchingWeightError,
                                     resEx, resAp, self.HFEngine, False,
                                     self.trainedModel.data.chordalRadius)
         pmR, pmC = pointMatching(dist)
         return np.mean(dist[pmR, pmC])
     
     def getErrorEstimatorMarginalLookAhead(self) -> Np1D:
         if not hasattr(self.trainedModel, "_musMExcl"):
             err = np.zeros(0)
             err[:] = np.inf
             self._musMarginalTestIdxs = np.zeros(0, dtype = int)
             return err
         self._musMarginalTestIdxs = np.array(self.trainedModel._idxExcl,
                                              dtype = int)
         idx, sizes = indicesScatter(len(self.trainedModel._musMExcl),
                                     return_sizes = True)
         err = []
         if len(idx) > 0:
             self.verbosity -= 25
             self.trainedModel.verbosity -= 25
             for j in idx:
                 muTest = self.trainedModel._musMExcl[j]
                 HITest = self.trainedModel._HIsExcl[j]
                 polesEx = HITest.poles
                 idxGood = np.isinf(polesEx) + np.isnan(polesEx) == False
                 polesEx = polesEx[idxGood]
                 if self.matchingWeightError != 0:
                     resEx = HITest.coeffs[np.where(idxGood)[0]]
                 else:
                     resEx = None
                 if len(polesEx) == 0:
                     err += [0.]
                     continue
                 err += [self._getDistanceApp(polesEx, resEx, muTest)]
             self.verbosity += 25
             self.trainedModel.verbosity += 25
         return arrayGatherv(np.array(err), sizes)
 
     def getErrorEstimatorMarginalNone(self) -> Np1D:
         nErr = len(self.trainedModel.data.musMarginal)
         self._musMarginalTestIdxs = np.arange(nErr)
         return (1. + self.greedyTolMarginal) * np.ones(nErr)
 
     def errorEstimatorMarginal(self, return_max : bool = False) -> Np1D:
         vbMng(self.trainedModel, "INIT",
               "Evaluating error estimator at mu = {}.".format(
                                        self.trainedModel.data.musMarginal), 10)
         if self.errorEstimatorKindMarginal == "NONE":
             nErr = len(self.trainedModel.data.musMarginal)
             self._musMarginalTestIdxs = np.arange(nErr)
             err = (1. + self.greedyTolMarginal) * np.ones(nErr)
         else:#if self.errorEstimatorKindMarginal[: 10] == "LOOK_AHEAD":
             err = self.getErrorEstimatorMarginalLookAhead()
         vbMng(self.trainedModel, "DEL", "Done evaluating error estimator.", 10)
         if not return_max: return err
         idxMaxEst = np.where(err > self.greedyTolMarginal)[0]
         maxErr = err[idxMaxEst]
         if self.errorEstimatorKindMarginal == "NONE": maxErr = None
         return err, idxMaxEst, maxErr
 
     def plotEstimatorMarginal(self, est:Np1D, idxMax:List[int],
                               estMax:List[float]):
         if self.errorEstimatorKindMarginal == "NONE": return
         if (not (np.any(np.isnan(est)) or np.any(np.isinf(est)))
         and masterCore() and hasattr(self.trainedModel, "_musMExcl")):
             fig = plt.figure(figsize = plt.figaspect(1. / self.nparMarginal))
             for jpar in range(self.nparMarginal):
                 ax = fig.add_subplot(1, self.nparMarginal, 1 + jpar)
                 musre = np.real(self.trainedModel._musMExcl)
                 if len(idxMax) > 0 and estMax is not None:
                     maxrej = musre[idxMax, jpar]
                 errCP = copy(est)
                 idx = np.delete(np.arange(self.nparMarginal), jpar)
                 while len(musre) > 0:
                     if self.nparMarginal == 1:
                         currIdx = np.arange(len(musre))
                     else:
                         currIdx = np.where(np.isclose(np.sum(
                                      np.abs(musre[:, idx] - musre[0, idx]), 1),
                                                       0., atol = 1e-15))[0]
                     currIdxSorted = currIdx[np.argsort(musre[currIdx, jpar])]
                     ax.semilogy(musre[currIdxSorted, jpar],
                                 errCP[currIdxSorted], 'k.-', linewidth = 1)
                     musre = np.delete(musre, currIdx, 0)
                     errCP = np.delete(errCP, currIdx)
                 ax.semilogy(self.musMarginal.re(jpar),
                             (self.greedyTolMarginal,) * len(self.musMarginal),
                             '*m')
                 if len(idxMax) > 0 and estMax is not None:
                     ax.semilogy(maxrej, estMax, 'xr')
                 ax.set_xlim(*list(self.samplerMarginal.lims.re(jpar)))
                 ax.grid()
             plt.tight_layout()
             plt.show()
 
     def _addMarginalSample(self, mus:paramList):
         mus = self.checkParameterListMarginal(mus)
         if len(mus) == 0: return
         self._nmusOld, nmus = len(self.musMarginal), len(mus)
         if (hasattr(self, "trainedModel") and self.trainedModel is not None
         and hasattr(self.trainedModel, "_musMExcl")):
             self._nmusOld += len(self.trainedModel._musMExcl)
         vbMng(self, "MAIN",
               ("Adding marginal sample point{} no. {}{} at {} to training "
                "set.").format("s" * (nmus > 1), self._nmusOld + 1,
                               "--{}".format(self._nmusOld + nmus) * (nmus > 1),
                               mus), 3)
         self.musMarginal.append(mus)
         self.setupApproxPivoted(mus)
         self._preliminaryMarginalFinalization()
         del self._nmusOld
         if (self.errorEstimatorKindMarginal[: 10] == "LOOK_AHEAD"
         and not self.firstGreedyIterM):
             ubRange = len(self.trainedModel.data.musMarginal)
             if hasattr(self.trainedModel, "_idxExcl"):
                 shRange = len(self.trainedModel._musMExcl)
             else:
                 shRange = 0
             testIdxs = list(range(ubRange + shRange - len(mus),
                                   ubRange + shRange))
             for j in testIdxs[::-1]:
                 self.musMarginal.pop(j - shRange)
             if hasattr(self.trainedModel, "_idxExcl"):
                 testIdxs = self.trainedModel._idxExcl + testIdxs
             self._updateTrainedModelMarginalSamples(testIdxs)
         self._finalizeMarginalization()
         self._SMarginal = len(self.musMarginal)
         self._approxParameters["SMarginal"] = self.SMarginal
         self.trainedModel.data.approxParameters["SMarginal"] = self.SMarginal
 
     def greedyNextSampleMarginal(self, muidx:List[int],
                                  plotEst : str = "NONE") \
                                     -> Tuple[Np1D, List[int], float, paramVal]:
         RROMPyAssert(self._mode, message = "Cannot add greedy sample.")
         muidx = self._musMarginalTestIdxs[muidx]
         if (self.errorEstimatorKindMarginal[: 10] == "LOOK_AHEAD"
         and not self.firstGreedyIterM):
             if not hasattr(self.trainedModel, "_idxExcl"):
                 raise RROMPyException(("Sample index to be added not present "
                                        "in trained model."))
             testIdxs = copy(self.trainedModel._idxExcl)
             skippedIdx = 0
             for cj, j in enumerate(self.trainedModel._idxExcl):
                 if j in muidx:
                     testIdxs.pop(skippedIdx)
                     self.musMarginal.insert(self.trainedModel._musMExcl[cj],
                                             j - skippedIdx)
                 else:
                     skippedIdx += 1
             if len(self.trainedModel._idxExcl) < (len(muidx)
                                                 + len(testIdxs)):
                 raise RROMPyException(("Sample index to be added not present "
                                        "in trained model."))
             self._updateTrainedModelMarginalSamples(testIdxs)
             self._SMarginal = len(self.musMarginal)
             self._approxParameters["SMarginal"] = self.SMarginal
             self.trainedModel.data.approxParameters["SMarginal"] = (
                                                                 self.SMarginal)
         self.firstGreedyIterM = False
         idxAdded = self.samplerMarginal.refine(muidx)[0]
         self._addMarginalSample(self.samplerMarginal.points[idxAdded])
         errorEstTest, muidx, maxErrorEst = self.errorEstimatorMarginal(True)
         if plotEst == "ALL":
             self.plotEstimatorMarginal(errorEstTest, muidx, maxErrorEst)
         return (errorEstTest, muidx, maxErrorEst,
                 self.samplerMarginal.points[muidx])
                                               
     def _preliminaryTrainingMarginal(self):
         """Initialize starting snapshots of solution map."""
         RROMPyAssert(self._mode, message = "Cannot start greedy algorithm.")
         if np.sum(self.samplingEngine.nsamples) > 0: return
         self.resetSamples()
         self._addMarginalSample(self.samplerMarginal.generatePoints(
                                                                self.SMarginal))
 
     def _preSetupApproxPivoted(self, mus:paramList) \
                                            -> Tuple[ListAny, ListAny, ListAny]:
         self.computeScaleFactor()
         if self.trainedModel is None:
             self._setupTrainedModel(np.zeros((0, 0)))
             self.trainedModel.data.Qs, self.trainedModel.data.Ps = [], []
             self.trainedModel.data.Psupp = []
         self._trainedModelOld = copy(self.trainedModel)
         self._scaleFactorOldPivot = copy(self.scaleFactor)
         self.scaleFactor = self.scaleFactorPivot
         self._temporaryPivot = 1
         self._musLoc = copy(self.mus)
         idx, sizes = indicesScatter(len(mus), return_sizes = True)
         emptyCores = np.where(sizes == 0)[0]
         self.verbosity -= 10
         self.samplingEngine.verbosity -= 10
         return idx, sizes, emptyCores
 
     def _postSetupApproxPivoted(self, mus:Np2D, pMat:Np2D, Ps:ListAny,
                                 Qs:ListAny, sizes:ListAny):
         self.scaleFactor = self._scaleFactorOldPivot
         del self._scaleFactorOldPivot, self._temporaryPivot
         pMat, Ps, Qs, mus, nsamples = gatherPivotedApproximant(pMat, Ps, Qs,
                                                                mus, sizes,
                                                                self.polybasis)
         if len(self._musLoc) > 0:
             self._mus = self.checkParameterList(self._musLoc)
             self._mus.append(mus)
         else:
             self._mus = self.checkParameterList(mus)
         self.trainedModel = self._trainedModelOld
         del self._trainedModelOld
         padLeft = self.trainedModel.data.projMat.shape[1]
         suppNew = np.append(0, np.cumsum(nsamples))
         self._setupTrainedModel(pMat, padLeft > 0)
         self.trainedModel.data.Qs += Qs
         self.trainedModel.data.Ps += Ps
         self.trainedModel.data.Psupp += list(padLeft + suppNew[: -1])
         self.trainedModel.data.approxParameters = copy(self.approxParameters)
         self.verbosity += 10
         self.samplingEngine.verbosity += 10
 
     def _localPivotedResult(self, pMat:Np2D, req:ListAny, emptyCores:ListAny,
                             mus:Np2D) -> Tuple[Np2D, ListAny, Np2D]:
         pMati = self.samplingEngine.projectionMatrix
         musi = self.samplingEngine.mus
         if not hasattr(self, "matchState") or not self.matchState:
             if self.POD == 1 and not (
                 hasattr(self.HFEngine.C, "is_mu_independent")
             and self.HFEngine.C.is_mu_independent in self._output_lvl):
                 raise RROMPyException(("Cannot apply mu-dependent C "
                                        "to orthonormalized samples."))
             vbMng(self, "INIT", "Extracting system output from state.", 35)
             pMatiEff = None
             for j, mu in enumerate(musi):
                 pMij = np.expand_dims(self.HFEngine.applyC(pMati[:, j], mu),
                                       -1)
                 if pMatiEff is None:
                     pMatiEff = np.array(pMij)
                 else:
                     pMatiEff = np.append(pMatiEff, pMij, axis = 1)
             pMati = pMatiEff
             vbMng(self, "DEL", "Done extracting system output.", 35)
         if pMat is None:
             mus = copy(musi.data)
             pMat = copy(pMati)
             if masterCore():
                 for dest in emptyCores:
                     req += [isend((len(pMat), pMat.dtype, mus.dtype),
                                   dest = dest, tag = dest)]
         else:
             mus = np.vstack((mus, musi.data))
             pMat = np.hstack((pMat, pMati))
         return pMat, req, mus
     
     @abstractmethod
     def setupApproxPivoted(self, mus:paramList) -> int:
         if self.checkComputedApproxPivoted(): return -1
         RROMPyAssert(self._mode, message = "Cannot setup approximant.")
         vbMng(self, "INIT", "Setting up pivoted approximant.", 10)
         self._preSetupApproxPivoted()
         data = []
         pass
         self._postSetupApproxPivoted(mus, data)
         vbMng(self, "DEL", "Done setting up pivoted approximant.", 10)
         return 0
 
     def setupApprox(self, plotEst : str = "NONE") -> int:
         """Compute greedy snapshots of solution map."""
         if self.checkComputedApprox(): return -1
         RROMPyAssert(self._mode, message = "Cannot start greedy algorithm.")
         vbMng(self, "INIT", "Setting up {}.". format(self.name()), 5)
         vbMng(self, "INIT", "Starting computation of snapshots.", 5)
         max2ErrorEst, self.firstGreedyIterM = np.inf, True
         self._preliminaryTrainingMarginal()
         if self.errorEstimatorKindMarginal == "NONE":
             muidx = []
         else:#if self.errorEstimatorKindMarginal[: 10] == "LOOK_AHEAD":
             muidx = np.arange(len(self.trainedModel.data.musMarginal))
         self._musMarginalTestIdxs = np.array(muidx)
         while self.firstGreedyIterM or (max2ErrorEst > self.greedyTolMarginal
                       and self.samplerMarginal.npoints < self.maxIterMarginal):
             errorEstTest, muidx, maxErrorEst, mu = \
                                   self.greedyNextSampleMarginal(muidx, plotEst)
             if maxErrorEst is None:
                 max2ErrorEst = 1. + self.greedyTolMarginal
             else:
                 if len(maxErrorEst) > 0:
                     max2ErrorEst = np.max(maxErrorEst)
                 else:
                     max2ErrorEst = np.max(errorEstTest)
                 vbMng(self, "MAIN", ("Uniform testing error estimate "
                                      "{:.4e}.").format(max2ErrorEst), 5)
         if plotEst == "LAST":
             self.plotEstimatorMarginal(errorEstTest, muidx, maxErrorEst)
         vbMng(self, "DEL", ("Done computing snapshots (final snapshot count: "
                             "{}).").format(len(self.mus)), 5)
         if (self.errorEstimatorKindMarginal == "LOOK_AHEAD_RECOVER"
         and hasattr(self.trainedModel, "_idxExcl")
         and len(self.trainedModel._idxExcl) > 0):
             vbMng(self, "INIT", "Recovering {} test models.".format(
                                            len(self.trainedModel._idxExcl)), 7)
             for j, mu in zip(self.trainedModel._idxExcl,
                              self.trainedModel._musMExcl):
                 self.musMarginal.insert(mu, j)
             self._preliminaryMarginalFinalization()
             self._updateTrainedModelMarginalSamples()
             self._finalizeMarginalization()
             self._SMarginal = len(self.musMarginal)
             self._approxParameters["SMarginal"] = self.SMarginal
             self.trainedModel.data.approxParameters["SMarginal"] = (
                                                                 self.SMarginal)
             vbMng(self, "DEL", "Done recovering test models.", 7)
         vbMng(self, "DEL", "Done setting up approximant.", 5)
         return 0
 
     def checkComputedApproxPivoted(self) -> bool:
         return (super().checkComputedApprox()
           and len(self.musMarginal) == len(self.trainedModel.data.musMarginal))
 
 class GenericPivotedGreedyApproximantPoleMatch(
                                            GenericPivotedGreedyApproximantBase,
                                            GenericPivotedApproximantPoleMatch):
     """
     ROM pivoted greedy interpolant computation for parametric problems (with
         pole matching) (ABSTRACT).
 
     Args:
         HFEngine: HF problem solver.
         mu0(optional): Default parameter. Defaults to 0.
         directionPivot(optional): Pivot components. Defaults to [0].
         approxParameters(optional): Dictionary containing values for main
             parameters of approximant. Recognized keys are:
             - 'POD': kind of snapshots orthogonalization; allowed values
                 include 0, 1/2, and 1; defaults to 1, i.e. POD;
             - 'scaleFactorDer': scaling factors for derivative computation;
                 defaults to 'AUTO';
             - 'matchState': whether to match the system state rather than the
                 system output; defaults to False;
             - 'matchingWeight': weight for pole matching optimization; defaults
                 to 1;
             - 'matchingChordalRadius': radius to be used in chordal metric for
                 poles and residues; if <= 0, Euclidean metric is used; if
                 'AUTO', automatically selected; defaults to -1;
             - 'matchingShared': required ratio of marginal points to share
                 resonance; defaults to 1.;
-            - 'badPoleCorrectionTol': tolerance for correction of bad poles;
-                defaults to 1., i.e. all corrections allowed;
+            - 'badPoleCorrection': strategy for correction of bad poles;
+                available values include 'ERASE', 'RATIONAL', and 'POLYNOMIAL';
+                defaults to 'ERASE';
             - 'matchingWeightError': weight for pole matching optimization in
                 error estimation; defaults to 0;
             - 'S': total number of pivot samples current approximant relies
                 upon;
             - 'samplerPivot': pivot sample point generator;
             - 'SMarginal': number of starting marginal samples;
             - 'samplerMarginal': marginal sample point generator via sparse
                 grid;
             - 'errorEstimatorKindMarginal': kind of marginal error estimator;
                 available values include 'LOOK_AHEAD', 'LOOK_AHEAD_RECOVER',
                 and 'NONE'; defaults to 'NONE';
             - 'polybasisMarginal': type of polynomial basis for marginal
                 interpolation; allowed values include 'MONOMIAL_*',
                 'CHEBYSHEV_*', 'LEGENDRE_*', 'NEARESTNEIGHBOR', and 
                 'PIECEWISE_LINEAR_*'; defaults to 'MONOMIAL';
             - 'paramsMarginal': dictionary of parameters for marginal
                 interpolation; include:
                 . 'MMarginal': degree of marginal interpolant; defaults to
                     'AUTO', i.e. maximum allowed; not for 'NEARESTNEIGHBOR' or
                     'PIECEWISE_LINEAR_*';
                 . 'nNeighborsMarginal': number of marginal nearest neighbors;
                      defaults to 1; only for 'NEARESTNEIGHBOR';
                 . 'polydegreetypeMarginal': type of polynomial degree for
                     marginal; defaults to 'TOTAL'; not for 'NEARESTNEIGHBOR' or
                     'PIECEWISE_LINEAR_*';
                 . 'interpTolMarginal': tolerance for marginal interpolation;
                     defaults to None; not for 'NEARESTNEIGHBOR' or
                     'PIECEWISE_LINEAR_*';
                 . 'radialDirectionalWeightsMarginalAdapt': bounds for adaptive
                     rescaling of marginal radial basis weights; only for
                     radial basis.
             - 'greedyTolMarginal': uniform error tolerance for marginal greedy
                algorithm; defaults to 1e-1;
             - 'maxIterMarginal': maximum number of marginal greedy steps;
                defaults to 1e2;
             - 'radialDirectionalWeightsMarginal': radial basis weights for
                 marginal interpolant; defaults to 1.
             Defaults to empty dict.
         verbosity(optional): Verbosity level. Defaults to 10.
 
     Attributes:
         HFEngine: HF problem solver.
         mu0: Default parameter.
         directionPivot: Pivot components.
         mus: Array of snapshot parameters.
         musMarginal: Array of marginal snapshot parameters.
         approxParameters: Dictionary containing values for main parameters of
             approximant. Recognized keys are in parameterList.
         parameterListSoft: Recognized keys of soft approximant parameters:
             - 'POD': kind of snapshots orthogonalization;
             - 'scaleFactorDer': scaling factors for derivative computation;
             - 'matchState': whether to match the system state rather than the
                 system output;
             - 'matchingWeight': weight for pole matching optimization;
             - 'matchingChordalRadius': radius to be used in chordal metric for
                 poles and residues;
             - 'matchingShared': required ratio of marginal points to share
                 resonance;
             - 'badPoleCorrection': strategy for correction of bad poles;
             - 'matchingWeightError': weight for pole matching optimization in
                 error estimation;
             - 'errorEstimatorKindMarginal': kind of marginal error estimator;
             - 'polybasisMarginal': type of polynomial basis for marginal
                 interpolation;
             - 'paramsMarginal': dictionary of parameters for marginal
                 interpolation; include:
                 . 'MMarginal': degree of marginal interpolant;
                 . 'nNeighborsMarginal': number of marginal nearest neighbors;
                 . 'polydegreetypeMarginal': type of polynomial degree for
                     marginal;
                 . 'interpTolMarginal': tolerance for marginal interpolation;
                 . 'radialDirectionalWeightsMarginalAdapt': bounds for adaptive
                     rescaling of marginal radial basis weights.
             - 'greedyTolMarginal': uniform error tolerance for marginal greedy
                 algorithm;
             - 'maxIterMarginal': maximum number of marginal greedy steps;
             - 'radialDirectionalWeightsMarginal': radial basis weights for
                 marginal interpolant.
         parameterListCritical: Recognized keys of critical approximant
             parameters:
             - 'S': total number of pivot samples current approximant relies
                 upon;
             - 'samplerPivot': pivot sample point generator;
             - 'SMarginal': total number of marginal samples current approximant
                 relies upon;
             - 'samplerMarginal': marginal sample point generator via sparse
                 grid.
         verbosity: Verbosity level.
         POD: Kind of snapshots orthogonalization.
         scaleFactorDer: Scaling factors for derivative computation.
         matchState: Whether to match the system state rather than the system
             output.
         matchingWeight: Weight for pole matching optimization.
         matchingChordalRadius: Radius to be used in chordal metric for poles
             and residues.
         matchingShared: Required ratio of marginal points to share resonance.
         badPoleCorrection: Strategy for correction of bad poles.
         matchingWeightError: Weight for pole matching optimization in error
             estimation.
         S: Total number of pivot samples current approximant relies upon.
         samplerPivot: Pivot sample point generator.
         SMarginal: Total number of marginal samples current approximant relies
             upon.
         samplerMarginal: Marginal sample point generator via sparse grid.
         errorEstimatorKindMarginal: Kind of marginal error estimator.
         polybasisMarginal: Type of polynomial basis for marginal interpolation.
         paramsMarginal: Dictionary of parameters for marginal interpolation.
         greedyTolMarginal: Uniform error tolerance for marginal greedy
             algorithm.
         maxIterMarginal: Maximum number of marginal greedy steps.
         radialDirectionalWeightsMarginal: Radial basis weights for marginal
             interpolant.
         muBounds: list of bounds for pivot parameter values.
         muBoundsMarginal: list of bounds for marginal parameter values.
         samplingEngine: Sampling engine.
         uHF: High fidelity solution(s) with parameter(s) lastSolvedHF as
             sampleList.
         lastSolvedHF: Parameter(s) corresponding to last computed high fidelity
             solution(s) as parameterList.
         uApproxReduced: Reduced approximate solution(s) with parameter(s)
             lastSolvedApprox as sampleList.
         lastSolvedApproxReduced: Parameter(s) corresponding to last computed
             reduced approximate solution(s) as parameterList.
         uApprox: Approximate solution(s) with parameter(s) lastSolvedApprox as
             sampleList.
         lastSolvedApprox: Parameter(s) corresponding to last computed
             approximate solution(s) as parameterList.
     """
 
     def _updateTrainedModelMarginalSamples(self, idx : ListAny = []):
         self.trainedModel.updateEffectiveSamples(idx, self.matchingWeight,
                                                  self.HFEngine, False,
                                                  self.matchingChordalRadius)
 
     def setupApprox(self, *args, **kwargs) -> int:
         if self.checkComputedApprox(): return -1
         self.purgeparamsMarginal()
         _polybasisMarginal = self.polybasisMarginal
         self._polybasisMarginal = ("PIECEWISE_LINEAR_"
                                  + self.samplerMarginal.kind)
         setupOK = super().setupApprox(*args, **kwargs)
         self._polybasisMarginal = _polybasisMarginal
         if self.matchState: self._postApplyC()
         return setupOK
diff --git a/rrompy/reduction_methods/pivoted/greedy/rational_interpolant_greedy_pivoted_greedy.py b/rrompy/reduction_methods/pivoted/greedy/rational_interpolant_greedy_pivoted_greedy.py
index 2070725..66760b6 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,361 +1,361 @@
 #Copyright (C) 2018-2020 by the RROMPy authors
 #
 # This file is part of RROMPy.
 #
 # RROMPy is free software: you can redistribute it and/or modify
 # it under the terms of the GNU Lesser General Public License as published by
 # the Free Software Foundation, either version 3 of the License, or
 # (at your option) any later version.
 #
 # RROMPy is distributed in the hope that it will be useful,
 # but WITHOUT ANY WARRANTY; without even the implied warranty of
 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
 # GNU Lesser General Public License for more details.
 #
 # You should have received a copy of the GNU Lesser General Public License
 # along with RROMPy. If not, see <http://www.gnu.org/licenses/>.
 #
 
 from copy import deepcopy as copy
 import numpy as np
 from .generic_pivoted_greedy_approximant import (
                                       GenericPivotedGreedyApproximantBase,
                                       GenericPivotedGreedyApproximantPoleMatch)
 from rrompy.reduction_methods.standard.greedy import RationalInterpolantGreedy
 from rrompy.reduction_methods.standard.greedy.generic_greedy_approximant \
                                                             import pruneSamples
 from rrompy.reduction_methods.pivoted import (
                                      RationalInterpolantGreedyPivotedPoleMatch)
 from rrompy.utilities.base.types import Np1D, Tuple, paramVal, paramList
 from rrompy.utilities.base import verbosityManager as vbMng
 from rrompy.utilities.exception_manager import RROMPyAssert
 from rrompy.parameter import emptyParameterList
 from rrompy.utilities.parallel import poolRank, recv
 
 __all__ = ['RationalInterpolantGreedyPivotedGreedyPoleMatch']
 
 class RationalInterpolantGreedyPivotedGreedyBase(
                                           GenericPivotedGreedyApproximantBase):
     @property
     def sampleBatchSize(self):
         """Value of sampleBatchSize."""
         return 1
 
     @property
     def sampleBatchIdx(self):
         """Value of sampleBatchIdx."""
         return self.S
 
     def greedyNextSample(self, muidx:int, plotEst : str = "NONE")\
                                           -> Tuple[Np1D, int, float, paramVal]:
         """Compute next greedy snapshot of solution map."""
         RROMPyAssert(self._mode, message = "Cannot add greedy sample.")
         mus = copy(self.muTest[muidx])
         self.muTest.pop(muidx)
         for j, mu in enumerate(mus):
             vbMng(self, "MAIN",
                   ("Adding sample point no. {} at {} to training "
                    "set.").format(len(self.mus) + 1, mu), 3)
             self.mus.append(mu)
             self._S = len(self.mus)
             self._approxParameters["S"] = self.S
             if (self.samplingEngine.nsamples <= len(mus) - j - 1
              or not np.allclose(mu, self.samplingEngine.mus[j - len(mus)])):
                 self.samplingEngine.nextSample(mu)
             if self._isLastSampleCollinear():
                 vbMng(self, "MAIN",
                       ("Collinearity above tolerance detected. Starting "
                        "preemptive greedy loop termination."), 3)
                 self._collinearityFlag = 1
                 errorEstTest = np.empty(len(self.muTest))
                 errorEstTest[:] = np.nan
                 return errorEstTest, [-1], np.nan, np.nan
         errorEstTest, muidx, maxErrorEst = self.errorEstimator(self.muTest,
                                                                True)
         if plotEst == "ALL":
             self.plotEstimator(errorEstTest, muidx, maxErrorEst)
         return errorEstTest, muidx, maxErrorEst, self.muTest[muidx]
 
     def _setSampleBatch(self, maxS:int):
         return self.S
 
     def _preliminaryTraining(self):
         """Initialize starting snapshots of solution map."""
         RROMPyAssert(self._mode, message = "Cannot start greedy algorithm.")
         if self.samplingEngine.nsamples > 0: return
         self.resetSamples()
         self.samplingEngine.scaleFactor = self.scaleFactorDer
         musPivot = self.samplerTrainSet.generatePoints(self.S)
         while len(musPivot) > self.S: musPivot.pop()
         muTestBasePivot = self.samplerPivot.generatePoints(self.nTestPoints,
                                                            False)
         idxPop = pruneSamples(self.mapParameterListPivot(muTestBasePivot),
                               self.mapParameterListPivot(musPivot),
                               1e-10 * self.scaleFactorPivot[0])
         muTestBasePivot.pop(idxPop)
         self._mus = emptyParameterList()
         self.mus.reset((self.S - 1, self.HFEngine.npar))
         self.muTest = emptyParameterList()
         self.muTest.reset((len(muTestBasePivot) + 1, self.HFEngine.npar))
         self.mus.data[:, self.directionPivot] = musPivot[: -1]
         self.mus.data[:, self.directionMarginal] = np.repeat(self.muMargLoc,
                                                           self.S - 1, axis = 0)
         self.muTest.data[: -1, self.directionPivot] = muTestBasePivot.data
         self.muTest.data[-1, self.directionPivot] = musPivot[-1]
         self.muTest.data[:, self.directionMarginal] = np.repeat(self.muMargLoc,
                                                       len(muTestBasePivot) + 1,
                                                       axis = 0)
         if len(self.mus) > 0:
             vbMng(self, "MAIN", 
                   ("Adding first {} sample point{} at {} to training "
                    "set.").format(self.S - 1, "" + "s" * (self.S > 2),
                                   self.mus), 3)
             self.samplingEngine.iterSample(self.mus)
         self._S = len(self.mus)
         self._approxParameters["S"] = self.S
         self.M, self.N = ("AUTO",) * 2
 
     def setupApproxPivoted(self, mus:paramList) -> int:
         if self.checkComputedApproxPivoted(): return -1
         RROMPyAssert(self._mode, message = "Cannot setup approximant.")
         vbMng(self, "INIT", "Setting up pivoted approximant.", 10)
         if not hasattr(self, "_plotEstPivot"): self._plotEstPivot = "NONE"
         idx, sizes, emptyCores = self._preSetupApproxPivoted(mus)
         S0 = copy(self.S)
         pMat, Ps, Qs, req, musA = None, [], [], [], None
         if len(idx) == 0:
             vbMng(self, "MAIN", "Idling.", 45)
             if self.storeAllSamples: self.storeSamples()
             pL, pT, mT = recv(source = 0, tag = poolRank())
             pMat = np.empty((pL, 0), dtype = pT)
             musA = np.empty((0, self.mu0.shape[1]), dtype = mT)
         else:
             for i in idx:
                 self.muMargLoc = mus[[i]]
                 vbMng(self, "MAIN", "Building marginal model no. {} at "
                                     "{}.".format(i + 1, self.muMargLoc), 25)
                 self.samplingEngine.resetHistory()
                 self.trainedModel = None
                 self.verbosity -= 5
                 self.samplingEngine.verbosity -= 10
                 RationalInterpolantGreedy.setupApprox(self, self._plotEstPivot)
                 self.verbosity += 5
                 self.samplingEngine.verbosity += 10
                 if self.storeAllSamples: self.storeSamples(i + self._nmusOld)
                 pMat, req, musA = self._localPivotedResult(pMat, req,
                                                            emptyCores, musA)
                 Ps += [copy(self.trainedModel.data.P)]
                 Qs += [copy(self.trainedModel.data.Q)]
                 self._S = S0
             del self.muMargLoc
         for r in req: r.wait()
         self._postSetupApproxPivoted(musA, pMat, Ps, Qs, sizes)
         vbMng(self, "DEL", "Done setting up pivoted approximant.", 10)
         return 0
 
     def setupApprox(self, plotEst : str = "NONE") -> int:
         if self.checkComputedApprox(): return -1
         if '_' not in plotEst: plotEst = plotEst + "_NONE"
         plotEstM, self._plotEstPivot = plotEst.split("_")
         val = super().setupApprox(plotEstM)
         return val
 
 class RationalInterpolantGreedyPivotedGreedyPoleMatch(
                                     RationalInterpolantGreedyPivotedGreedyBase,
                                     GenericPivotedGreedyApproximantPoleMatch,
                                     RationalInterpolantGreedyPivotedPoleMatch):
     """
     ROM greedy pivoted greedy rational interpolant computation for parametric
         problems (with pole matching).
 
     Args:
         HFEngine: HF problem solver.
         mu0(optional): Default parameter. Defaults to 0.
         directionPivot(optional): Pivot components. Defaults to [0].
         approxParameters(optional): Dictionary containing values for main
             parameters of approximant. Recognized keys are:
             - 'POD': kind of snapshots orthogonalization; allowed values
                 include 0, 1/2, and 1; defaults to 1, i.e. POD;
             - 'scaleFactorDer': scaling factors for derivative computation;
                 defaults to 'AUTO';
             - 'matchState': whether to match the system state rather than the
                 system output; defaults to False;
             - 'matchingWeight': weight for pole matching optimization; defaults
                 to 1;
             - 'matchingChordalRadius': radius to be used in chordal metric for
                 poles and residues; if <= 0, Euclidean metric is used; if
                 'AUTO', automatically selected; defaults to -1;
             - 'matchingShared': required ratio of marginal points to share
                 resonance; defaults to 1.;
             - 'badPoleCorrection': strategy for correction of bad poles;
                 available values include 'ERASE', 'RATIONAL', and 'POLYNOMIAL';
                 defaults to 'ERASE';
             - 'matchingWeightError': weight for pole matching optimization in
                 error estimation; defaults to 0;
             - 'S': total number of pivot samples current approximant relies
                 upon;
             - 'samplerPivot': pivot sample point generator;
             - 'SMarginal': number of starting marginal samples;
             - 'samplerMarginal': marginal sample point generator via sparse
                 grid;
             - 'errorEstimatorKindMarginal': kind of marginal error estimator;
                 available values include 'LOOK_AHEAD' and 'LOOK_AHEAD_RECOVER';
                 defaults to 'NONE';
             - 'polybasis': type of polynomial basis for pivot interpolation;
                 defaults to 'MONOMIAL';
             - 'polybasisMarginal': type of polynomial basis for marginal
                 interpolation; allowed values include 'MONOMIAL_*',
                 'CHEBYSHEV_*', 'LEGENDRE_*', 'NEARESTNEIGHBOR', and 
                 'PIECEWISE_LINEAR_*'; defaults to 'MONOMIAL';
             - 'paramsMarginal': dictionary of parameters for marginal
                 interpolation; include:
                 . 'MMarginal': degree of marginal interpolant; defaults to
                     'AUTO', i.e. maximum allowed; not for 'NEARESTNEIGHBOR' or
                     'PIECEWISE_LINEAR_*';
                 . 'nNeighborsMarginal': number of marginal nearest neighbors;
                      defaults to 1; only for 'NEARESTNEIGHBOR';
                 . 'polydegreetypeMarginal': type of polynomial degree for
                     marginal; defaults to 'TOTAL'; not for 'NEARESTNEIGHBOR' or
                     'PIECEWISE_LINEAR_*';
                 . 'interpTolMarginal': tolerance for marginal interpolation;
                     defaults to None; not for 'NEARESTNEIGHBOR' or
                     'PIECEWISE_LINEAR_*';
                 . 'radialDirectionalWeightsMarginalAdapt': bounds for adaptive
                     rescaling of marginal radial basis weights; only for
                     radial basis.
             - 'greedyTol': uniform error tolerance for greedy algorithm;
                 defaults to 1e-2;
             - 'collinearityTol': collinearity tolerance for greedy algorithm;
                 defaults to 0.;
             - 'maxIter': maximum number of greedy steps; defaults to 1e2;
             - 'nTestPoints': number of test points; defaults to 5e2;
             - 'samplerTrainSet': training sample points generator; defaults to
                 samplerPivot;
             - 'errorEstimatorKind': kind of error estimator; available values
                 include 'AFFINE', 'DISCREPANCY', 'LOOK_AHEAD',
                 'LOOK_AHEAD_RES', and 'NONE'; defaults to 'NONE';
             - 'greedyTolMarginal': uniform error tolerance for marginal greedy
                algorithm; defaults to 1e-1;
             - 'maxIterMarginal': maximum number of marginal greedy steps;
                defaults to 1e2;
             - 'radialDirectionalWeightsMarginal': radial basis weights for
                 marginal interpolant; defaults to 1;
             - 'functionalSolve': strategy for minimization of denominator
-                functional; allowed values include 'NORM', 'DOMINANT', 'NODAL',
-                'BARYCENTRIC_NORM', and 'BARYCENTRIC[_AVERAGE]' (check pdf in
-                main folder for meaning); defaults to 'NORM';
+                functional; allowed values include 'NORM', 'DOMINANT',
+                'BARYCENTRIC[_NORM]', and 'BARYCENTRIC_AVERAGE' (check pdf in
+                main folder for explanation); defaults to 'NORM';
             - 'interpTol': tolerance for pivot interpolation; defaults to None;
             - 'QTol': tolerance for robust rational denominator management;
                 defaults to 0.
             Defaults to empty dict.
         verbosity(optional): Verbosity level. Defaults to 10.
             
     Attributes:
         HFEngine: HF problem solver.
         mu0: Default parameter.
         directionPivot: Pivot components.
         mus: Array of snapshot parameters.
         musMarginal: Array of marginal snapshot parameters.
         approxParameters: Dictionary containing values for main parameters of
             approximant. Recognized keys are in parameterList.
         parameterListSoft: Recognized keys of soft approximant parameters:
             - 'POD': kind of snapshots orthogonalization;
             - 'scaleFactorDer': scaling factors for derivative computation;
             - 'matchState': whether to match the system state rather than the
                 system output;
             - 'matchingWeight': weight for pole matching optimization;
             - 'matchingChordalRadius': radius to be used in chordal metric for
                 poles and residues;
             - 'matchingShared': required ratio of marginal points to share
                 resonance;
             - 'badPoleCorrection': strategy for correction of bad poles;
             - 'matchingWeightError': weight for pole matching optimization in
                 error estimation;
             - 'errorEstimatorKindMarginal': kind of marginal error estimator;
             - 'polybasis': type of polynomial basis for pivot interpolation;
             - 'polybasisMarginal': type of polynomial basis for marginal
                 interpolation;
             - 'paramsMarginal': dictionary of parameters for marginal
                 interpolation; include:
                 . 'MMarginal': degree of marginal interpolant;
                 . 'nNeighborsMarginal': number of marginal nearest neighbors;
                 . 'polydegreetypeMarginal': type of polynomial degree for
                     marginal;
                 . 'interpTolMarginal': tolerance for marginal interpolation;
                 . 'radialDirectionalWeightsMarginalAdapt': bounds for adaptive
                     rescaling of marginal radial basis weights.
             - 'greedyTol': uniform error tolerance for greedy algorithm;
             - 'collinearityTol': collinearity tolerance for greedy algorithm;
             - 'maxIter': maximum number of greedy steps;
             - 'nTestPoints': number of test points;
             - 'samplerTrainSet': training sample points generator;
             - 'errorEstimatorKind': kind of error estimator;
             - 'greedyTolMarginal': uniform error tolerance for marginal greedy
                 algorithm;
             - 'maxIterMarginal': maximum number of marginal greedy steps;
             - 'radialDirectionalWeightsMarginal': radial basis weights for
                 marginal interpolant;
             - 'functionalSolve': strategy for minimization of denominator
                 functional;
             - 'interpTol': tolerance for pivot interpolation;
             - 'QTol': tolerance for robust rational denominator management.
         parameterListCritical: Recognized keys of critical approximant
             parameters:
             - 'S': total number of pivot samples current approximant relies
                 upon;
             - 'samplerPivot': pivot sample point generator;
             - 'SMarginal': total number of marginal samples current approximant
                 relies upon;
             - 'samplerMarginal': marginal sample point generator via sparse
                 grid.
         verbosity: Verbosity level.
         POD: Kind of snapshots orthogonalization.
         scaleFactorDer: Scaling factors for derivative computation.
         matchState: Whether to match the system state rather than the system
             output.
         matchingWeight: Weight for pole matching optimization.
         matchingChordalRadius: Radius to be used in chordal metric for poles
             and residues.
         matchingShared: Required ratio of marginal points to share resonance.
         badPoleCorrection: Strategy for correction of bad poles.
         matchingWeightError: Weight for pole matching optimization in error
             estimation.
         S: Total number of pivot samples current approximant relies upon.
         samplerPivot: Pivot sample point generator.
         SMarginal: Total number of marginal samples current approximant relies
             upon.
         samplerMarginal: Marginal sample point generator via sparse grid.
         errorEstimatorKindMarginal: Kind of marginal error estimator.
         polybasis: Type of polynomial basis for pivot interpolation.
         polybasisMarginal: Type of polynomial basis for marginal interpolation.
         paramsMarginal: Dictionary of parameters for marginal interpolation.
         greedyTol: uniform error tolerance for greedy algorithm.
         collinearityTol: Collinearity tolerance for greedy algorithm.
         maxIter: maximum number of greedy steps.
         nTestPoints: number of starting training points.
         samplerTrainSet: training sample points generator.
         errorEstimatorKind: kind of error estimator.
         greedyTolMarginal: Uniform error tolerance for marginal greedy
             algorithm.
         maxIterMarginal: Maximum number of marginal greedy steps.
         radialDirectionalWeightsMarginal: Radial basis weights for marginal
             interpolant.
         functionalSolve: Strategy for minimization of denominator functional.
         interpTol: Tolerance for pivot interpolation.
         QTol: Tolerance for robust rational denominator management.
         muBounds: list of bounds for pivot parameter values.
         muBoundsMarginal: list of bounds for marginal parameter values.
         samplingEngine: Sampling engine.
         uHF: High fidelity solution(s) with parameter(s) lastSolvedHF as
             sampleList.
         lastSolvedHF: Parameter(s) corresponding to last computed high fidelity
             solution(s) as parameterList.
         uApproxReduced: Reduced approximate solution(s) with parameter(s)
             lastSolvedApprox as sampleList.
         lastSolvedApproxReduced: Parameter(s) corresponding to last computed
             reduced approximate solution(s) as parameterList.
         uApprox: Approximate solution(s) with parameter(s) lastSolvedApprox as
             sampleList.
         lastSolvedApprox: Parameter(s) corresponding to last computed
             approximate solution(s) as parameterList.
     """
diff --git a/rrompy/reduction_methods/pivoted/greedy/rational_interpolant_pivoted_greedy.py b/rrompy/reduction_methods/pivoted/greedy/rational_interpolant_pivoted_greedy.py
index 3330a75..1aad812 100644
--- a/rrompy/reduction_methods/pivoted/greedy/rational_interpolant_pivoted_greedy.py
+++ b/rrompy/reduction_methods/pivoted/greedy/rational_interpolant_pivoted_greedy.py
@@ -1,295 +1,295 @@
 # Copyright (C) 2018-2020 by the RROMPy authors
 #
 # This file is part of RROMPy.
 #
 # RROMPy is free software: you can redistribute it and/or modify
 # it under the terms of the GNU Lesser General Public License as published by
 # the Free Software Foundation, either version 3 of the License, or
 # (at your option) any later version.
 #
 # RROMPy is distributed in the hope that it will be useful,
 # but WITHOUT ANY WARRANTY; without even the implied warranty of
 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
 # GNU Lesser General Public License for more details.
 #
 # You should have received a copy of the GNU Lesser General Public License
 # along with RROMPy. If not, see <http://www.gnu.org/licenses/>.
 #
 
 from copy import deepcopy as copy
 import numpy as np
 from .generic_pivoted_greedy_approximant import (
                                       GenericPivotedGreedyApproximantBase,
                                       GenericPivotedGreedyApproximantPoleMatch)
 from rrompy.reduction_methods.standard import RationalInterpolant
 from rrompy.reduction_methods.pivoted import (
                                            RationalInterpolantPivotedPoleMatch)
 from rrompy.utilities.base.types import paramList
 from rrompy.utilities.base import verbosityManager as vbMng
 from rrompy.utilities.exception_manager import RROMPyAssert
 from rrompy.parameter import emptyParameterList
 from rrompy.utilities.parallel import poolRank, recv
 
 __all__ = ['RationalInterpolantPivotedGreedyPoleMatch']
 
 class RationalInterpolantPivotedGreedyBase(
                                           GenericPivotedGreedyApproximantBase):
     
     def computeSnapshots(self):
         """Compute snapshots of solution map."""
         RROMPyAssert(self._mode,
                      message = "Cannot start snapshot computation.")
         vbMng(self, "INIT", "Starting computation of snapshots.", 5)
         self.samplingEngine.scaleFactor = self.scaleFactorDer
         if not hasattr(self, "musPivot") or len(self.musPivot) != self.S:
             self.musPivot = self.samplerPivot.generatePoints(self.S)
             while len(self.musPivot) > self.S: self.musPivot.pop()
         musLoc = emptyParameterList()
         musLoc.reset((self.S, self.HFEngine.npar))
         self.samplingEngine.resetHistory()
         musLoc.data[:, self.directionPivot] = self.musPivot.data
         musLoc.data[:, self.directionMarginal] = np.repeat(self.muMargLoc,
                                                            self.S, axis = 0)        
         self.samplingEngine.iterSample(musLoc)
         vbMng(self, "DEL", "Done computing snapshots.", 5)
         self._m_selfmus = copy(musLoc)
         self._mus = self.musPivot
         self._m_HFEparameterMap = copy(self.HFEngine.parameterMap)
         self.HFEngine.parameterMap = {
                 "F": [self.HFEngine.parameterMap["F"][self.directionPivot[0]]],
                 "B": [self.HFEngine.parameterMap["B"][self.directionPivot[0]]]}
 
     def setupApproxPivoted(self, mus:paramList) -> int:
         if self.checkComputedApproxPivoted(): return -1
         RROMPyAssert(self._mode, message = "Cannot setup approximant.")
         vbMng(self, "INIT", "Setting up pivoted approximant.", 10)
         idx, sizes, emptyCores = self._preSetupApproxPivoted(mus)
         pMat, Ps, Qs, req, musA = None, [], [], [], None
         if len(idx) == 0:
             vbMng(self, "MAIN", "Idling.", 45)
             if self.storeAllSamples: self.storeSamples()
             pL, pT, mT = recv(source = 0, tag = poolRank())
             pMat = np.empty((pL, 0), dtype = pT)
             musA = np.empty((0, self.mu0.shape[1]), dtype = mT)
         else:
             for i in idx:
                 self.muMargLoc = mus[[i]]
                 vbMng(self, "MAIN", "Building marginal model no. {} at "
                                     "{}.".format(i + 1, self.muMargLoc), 25)
                 self.samplingEngine.resetHistory()
                 self.trainedModel = None
                 self.verbosity -= 5
                 self.samplingEngine.verbosity -= 5
                 RationalInterpolant.setupApprox(self)
                 self.verbosity += 5
                 self.samplingEngine.verbosity += 5
                 self._mus = self._m_selfmus
                 self.HFEngine.parameterMap = self._m_HFEparameterMap
                 del self._m_selfmus, self._m_HFEparameterMap
                 if self.storeAllSamples: self.storeSamples(i + self._nmusOld)
                 pMat, req, musA = self._localPivotedResult(pMat, req,
                                                            emptyCores, musA)
                 Ps += [copy(self.trainedModel.data.P)]
                 Qs += [copy(self.trainedModel.data.Q)]
             del self.muMargLoc
         for r in req: r.wait()
         self._postSetupApproxPivoted(musA, pMat, Ps, Qs, sizes)
         vbMng(self, "DEL", "Done setting up pivoted approximant.", 10)
         return 0
 
 class RationalInterpolantPivotedGreedyPoleMatch(
                                       RationalInterpolantPivotedGreedyBase,
                                       GenericPivotedGreedyApproximantPoleMatch,
                                       RationalInterpolantPivotedPoleMatch):
     """
     ROM pivoted greedy rational interpolant computation for parametric
         problems (with pole matching).
 
     Args:
         HFEngine: HF problem solver.
         mu0(optional): Default parameter. Defaults to 0.
         directionPivot(optional): Pivot components. Defaults to [0].
         approxParameters(optional): Dictionary containing values for main
             parameters of approximant. Recognized keys are:
             - 'POD': kind of snapshots orthogonalization; allowed values
                 include 0, 1/2, and 1; defaults to 1, i.e. POD;
             - 'scaleFactorDer': scaling factors for derivative computation;
                 defaults to 'AUTO';
             - 'matchState': whether to match the system state rather than the
                 system output; defaults to False;
             - 'matchingWeight': weight for pole matching optimization; defaults
                 to 1;
             - 'matchingChordalRadius': radius to be used in chordal metric for
                 poles and residues; if <= 0, Euclidean metric is used; if
                 'AUTO', automatically selected; defaults to -1;
             - 'matchingShared': required ratio of marginal points to share
                 resonance; defaults to 1.;
             - 'badPoleCorrection': strategy for correction of bad poles;
                 available values include 'ERASE', 'RATIONAL', and 'POLYNOMIAL';
                 defaults to 'ERASE';
             - 'matchingWeightError': weight for pole matching optimization in
                 error estimation; defaults to 0;
             - 'S': total number of pivot samples current approximant relies
                 upon;
             - 'samplerPivot': pivot sample point generator;
             - 'SMarginal': number of starting marginal samples;
             - 'samplerMarginal': marginal sample point generator via sparse
                 grid;
             - 'errorEstimatorKindMarginal': kind of marginal error estimator;
                 available values include 'LOOK_AHEAD' and 'LOOK_AHEAD_RECOVER';
                 defaults to 'NONE';
             - 'polybasis': type of polynomial basis for pivot interpolation;
                 defaults to 'MONOMIAL';
             - 'polybasisMarginal': type of polynomial basis for marginal
                 interpolation; allowed values include 'MONOMIAL_*',
                 'CHEBYSHEV_*', 'LEGENDRE_*', 'NEARESTNEIGHBOR', and 
                 'PIECEWISE_LINEAR_*'; defaults to 'MONOMIAL';
             - 'paramsMarginal': dictionary of parameters for marginal
                 interpolation; include:
                 . 'MMarginal': degree of marginal interpolant; defaults to
                     'AUTO', i.e. maximum allowed; not for 'NEARESTNEIGHBOR' or
                     'PIECEWISE_LINEAR_*';
                 . 'nNeighborsMarginal': number of marginal nearest neighbors;
                      defaults to 1; only for 'NEARESTNEIGHBOR';
                 . 'polydegreetypeMarginal': type of polynomial degree for
                     marginal; defaults to 'TOTAL'; not for 'NEARESTNEIGHBOR' or
                     'PIECEWISE_LINEAR_*';
                 . 'interpTolMarginal': tolerance for marginal interpolation;
                     defaults to None; not for 'NEARESTNEIGHBOR' or
                     'PIECEWISE_LINEAR_*';
                 . 'radialDirectionalWeightsMarginalAdapt': bounds for adaptive
                     rescaling of marginal radial basis weights; only for
                     radial basis.
             - 'M': degree of rational interpolant numerator; defaults to
                 'AUTO', i.e. maximum allowed;
             - 'N': degree of rational interpolant denominator; defaults to
                 'AUTO', i.e. maximum allowed;
             - 'greedyTolMarginal': uniform error tolerance for marginal greedy
                algorithm; defaults to 1e-1;
             - 'maxIterMarginal': maximum number of marginal greedy steps;
                defaults to 1e2;
             - 'radialDirectionalWeights': radial basis weights for pivot
                 numerator; defaults to 1;
             - 'radialDirectionalWeightsAdapt': bounds for adaptive rescaling of
                 radial basis weights; defaults to [-1, -1];
             - 'radialDirectionalWeightsMarginal': radial basis weights for
                 marginal interpolant; defaults to 1;
             - 'functionalSolve': strategy for minimization of denominator
-                functional; allowed values include 'NORM', 'DOMINANT', 'NODAL',
-                'BARYCENTRIC_NORM', and 'BARYCENTRIC[_AVERAGE]' (check pdf in
-                main folder for meaning); defaults to 'NORM';
+                functional; allowed values include 'NORM', 'DOMINANT',
+                'BARYCENTRIC[_NORM]', and 'BARYCENTRIC_AVERAGE' (check pdf in
+                main folder for explanation); defaults to 'NORM';
             - 'interpTol': tolerance for pivot interpolation; defaults to None;
             - 'QTol': tolerance for robust rational denominator management;
                 defaults to 0.
             Defaults to empty dict.
         verbosity(optional): Verbosity level. Defaults to 10.
             
     Attributes:
         HFEngine: HF problem solver.
         mu0: Default parameter.
         directionPivot: Pivot components.
         mus: Array of snapshot parameters.
         musPivot: Array of pivot snapshot parameters.
         musMarginal: Array of marginal snapshot parameters.
         approxParameters: Dictionary containing values for main parameters of
             approximant. Recognized keys are in parameterList.
         parameterListSoft: Recognized keys of soft approximant parameters:
             - 'POD': kind of snapshots orthogonalization;
             - 'scaleFactorDer': scaling factors for derivative computation;
             - 'matchState': whether to match the system state rather than the
                 system output;
             - 'matchingWeight': weight for pole matching optimization;
             - 'matchingChordalRadius': radius to be used in chordal metric for
                 poles and residues;
             - 'matchingShared': required ratio of marginal points to share
                 resonance;
             - 'badPoleCorrection': strategy for correction of bad poles;
             - 'matchingWeightError': weight for pole matching optimization in
                 error estimation;
             - 'errorEstimatorKindMarginal': kind of marginal error estimator;
             - 'polybasis': type of polynomial basis for pivot interpolation;
             - 'polybasisMarginal': type of polynomial basis for marginal
                 interpolation;
             - 'paramsMarginal': dictionary of parameters for marginal
                 interpolation; include:
                 . 'MMarginal': degree of marginal interpolant;
                 . 'nNeighborsMarginal': number of marginal nearest neighbors;
                 . 'polydegreetypeMarginal': type of polynomial degree for
                     marginal;
                 . 'interpTolMarginal': tolerance for marginal interpolation;
                 . 'radialDirectionalWeightsMarginalAdapt': bounds for adaptive
                     rescaling of marginal radial basis weights.
             - 'M': degree of rational interpolant numerator;
             - 'N': degree of rational interpolant denominator;
             - 'greedyTolMarginal': uniform error tolerance for marginal greedy
                 algorithm;
             - 'maxIterMarginal': maximum number of marginal greedy steps;
             - 'radialDirectionalWeights': radial basis weights for pivot
                 numerator;
             - 'radialDirectionalWeightsAdapt': bounds for adaptive rescaling of
                 radial basis weights;
             - 'radialDirectionalWeightsMarginal': radial basis weights for
                 marginal interpolant;
             - 'functionalSolve': strategy for minimization of denominator
                 functional;
             - 'interpTol': tolerance for pivot interpolation;
             - 'QTol': tolerance for robust rational denominator management.
         parameterListCritical: Recognized keys of critical approximant
             parameters:
             - 'S': total number of pivot samples current approximant relies
                 upon;
             - 'samplerPivot': pivot sample point generator;
             - 'SMarginal': total number of marginal samples current approximant
                 relies upon;
             - 'samplerMarginal': marginal sample point generator via sparse
                 grid.
         verbosity: Verbosity level.
         POD: Kind of snapshots orthogonalization.
         scaleFactorDer: Scaling factors for derivative computation.
         matchState: Whether to match the system state rather than the system
             output.
         matchingWeight: Weight for pole matching optimization.
         matchingChordalRadius: Radius to be used in chordal metric for poles
             and residues.
         matchingShared: Required ratio of marginal points to share resonance.
         badPoleCorrection: Strategy for correction of bad poles.
         matchingWeightError: Weight for pole matching optimization in error
             estimation.
         S: Total number of pivot samples current approximant relies upon.
         samplerPivot: Pivot sample point generator.
         SMarginal: Total number of marginal samples current approximant relies
             upon.
         samplerMarginal: Marginal sample point generator via sparse grid.
         errorEstimatorKindMarginal: Kind of marginal error estimator.
         polybasis: Type of polynomial basis for pivot interpolation.
         polybasisMarginal: Type of polynomial basis for marginal interpolation.
         paramsMarginal: Dictionary of parameters for marginal interpolation.
         M: Degree of rational interpolant numerator.
         N: Degree of rational interpolant denominator.
         greedyTolMarginal: Uniform error tolerance for marginal greedy
             algorithm.
         maxIterMarginal: Maximum number of marginal greedy steps.
         radialDirectionalWeights: Radial basis weights for pivot numerator.
         radialDirectionalWeightsAdapt: Bounds for adaptive rescaling of radial
             basis weights.
         radialDirectionalWeightsMarginal: Radial basis weights for marginal
             interpolant.
         functionalSolve: Strategy for minimization of denominator functional.
         interpTol: Tolerance for pivot interpolation.
         QTol: Tolerance for robust rational denominator management.
         muBounds: list of bounds for pivot parameter values.
         muBoundsMarginal: list of bounds for marginal parameter values.
         samplingEngine: Sampling engine.
         uHF: High fidelity solution(s) with parameter(s) lastSolvedHF as
             sampleList.
         lastSolvedHF: Parameter(s) corresponding to last computed high fidelity
             solution(s) as parameterList.
         uApproxReduced: Reduced approximate solution(s) with parameter(s)
             lastSolvedApprox as sampleList.
         lastSolvedApproxReduced: Parameter(s) corresponding to last computed
             reduced approximate solution(s) as parameterList.
         uApprox: Approximate solution(s) with parameter(s) lastSolvedApprox as
             sampleList.
         lastSolvedApprox: Parameter(s) corresponding to last computed
             approximate solution(s) as parameterList.
     """
diff --git a/rrompy/reduction_methods/pivoted/rational_interpolant_greedy_pivoted.py b/rrompy/reduction_methods/pivoted/rational_interpolant_greedy_pivoted.py
index 626ef1d..2e7cb2b 100644
--- a/rrompy/reduction_methods/pivoted/rational_interpolant_greedy_pivoted.py
+++ b/rrompy/reduction_methods/pivoted/rational_interpolant_greedy_pivoted.py
@@ -1,571 +1,571 @@
 # Copyright (C) 2018-2020 by the RROMPy authors
 #
 # This file is part of RROMPy.
 #
 # RROMPy is free software: you can redistribute it and/or modify
 # it under the terms of the GNU Lesser General Public License as published by
 # the Free Software Foundation, either version 3 of the License, or
 # (at your option) any later version.
 #
 # RROMPy is distributed in the hope that it will be useful,
 # but WITHOUT ANY WARRANTY; without even the implied warranty of
 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
 # GNU Lesser General Public License for more details.
 #
 # You should have received a copy of the GNU Lesser General Public License
 # along with RROMPy. If not, see <http://www.gnu.org/licenses/>.
 #
 
 from copy import deepcopy as copy
 import numpy as np
 from .generic_pivoted_approximant import (GenericPivotedApproximantBase,
                                           GenericPivotedApproximantNoMatch,
                                           GenericPivotedApproximantPoleMatch)
 from .gather_pivoted_approximant import gatherPivotedApproximant
 from rrompy.reduction_methods.standard.greedy.rational_interpolant_greedy \
                                                import RationalInterpolantGreedy
 from rrompy.reduction_methods.standard.greedy.generic_greedy_approximant \
                                                             import pruneSamples
 from rrompy.utilities.base.types import Np1D
 from rrompy.utilities.base import verbosityManager as vbMng
 from rrompy.utilities.poly_fitting.polynomial import polyvander as pv
 from rrompy.utilities.exception_manager import RROMPyException, RROMPyAssert
 from rrompy.parameter import emptyParameterList
 from rrompy.utilities.parallel import poolRank, indicesScatter, isend, recv
 
 __all__ = ['RationalInterpolantGreedyPivotedNoMatch',
            'RationalInterpolantGreedyPivotedPoleMatch']
 
 class RationalInterpolantGreedyPivotedBase(GenericPivotedApproximantBase,
                                            RationalInterpolantGreedy):
     def __init__(self, *args, **kwargs):
         self._preInit()
         super().__init__(*args, **kwargs)
         if self.nparPivot > 1: self.HFEngine._ignoreResidues = 1
         self._postInit()
 
     @property
     def tModelType(self):
         if hasattr(self, "_temporaryPivot"):
             return RationalInterpolantGreedy.tModelType.fget(self)
         return super().tModelType
 
     def _polyvanderAuxiliary(self, mus, deg, *args):
         degEff = [0] * self.npar
         degEff[self.directionPivot[0]] = deg
         return pv(mus, degEff, *args)
 
     def _marginalizeMiscellanea(self, forward:bool):
         if forward:
             self._m_selfmus = copy(self.mus)
             self._m_HFEparameterMap = copy(self.HFEngine.parameterMap)
             self._mus = self.checkParameterListPivot(
                                                  self.mus(self.directionPivot))
             self.HFEngine.parameterMap = {
                 "F": [self.HFEngine.parameterMap["F"][self.directionPivot[0]]],
                 "B": [self.HFEngine.parameterMap["B"][self.directionPivot[0]]]}
         else:
             self._mus = self._m_selfmus
             self.HFEngine.parameterMap = self._m_HFEparameterMap
             del self._m_selfmus, self._m_HFEparameterMap
 
     def _marginalizeTrainedModel(self, forward:bool):
         if forward:
             del self._temporaryPivot
             self.trainedModel.data.mu0 = self.mu0
             self.trainedModel.data.scaleFactor = [1.] * self.npar
             self.trainedModel.data.scaleFactor[self.directionPivot[0]] = (
                                                            self.scaleFactor[0])
             self.trainedModel.data.parameterMap = self.HFEngine.parameterMap
             self._m_musUniqueCN = copy(self._musUniqueCN)
             musUniqueCNAux = np.zeros((self.S, self.npar),
                                       dtype = self._musUniqueCN.dtype)
             musUniqueCNAux[:, self.directionPivot[0]] = self._musUniqueCN(0)
             self._musUniqueCN = self.checkParameterList(musUniqueCNAux)
             self._m_derIdxs = copy(self._derIdxs)
             for j in range(len(self._derIdxs)):
                 for l in range(len(self._derIdxs[j])):
                     derjl = self._derIdxs[j][l][0]
                     self._derIdxs[j][l] = [0] * self.npar
                     self._derIdxs[j][l][self.directionPivot[0]] = derjl
             self.trainedModel.data.Q._dirPivot = self.directionPivot[0]
             self.trainedModel.data.P._dirPivot = self.directionPivot[0]
             # tell greedy error estimator that operator / RHS is pivot-affine
             if hasattr(self.HFEngine.A, "is_affine"):
                 self._A_is_affine = self.HFEngine.A.is_affine
             else:
                 self._A_is_affine = 0
             if hasattr(self.HFEngine.b, "is_affine"):
                 self._b_is_affine = self.HFEngine.b.is_affine
             else:
                 self._b_is_affine = 0
             if self._A_is_affine >= 1 / 2 and self._b_is_affine >= 1 / 2:
                 self._affine_lvl += [1 / 2]
         else:
             self._temporaryPivot = 1
             self.trainedModel.data.mu0 = self.checkParameterListPivot(
                                                  self.mu0(self.directionPivot))
             self.trainedModel.data.scaleFactor = self.scaleFactor
             self.trainedModel.data.parameterMap = {
                 "F": [self.HFEngine.parameterMap["F"][self.directionPivot[0]]],
                 "B": [self.HFEngine.parameterMap["B"][self.directionPivot[0]]]}
             self._musUniqueCN = copy(self._m_musUniqueCN)
             self._derIdxs = copy(self._m_derIdxs)
             del self._m_musUniqueCN, self._m_derIdxs
             del self.trainedModel.data.Q._dirPivot
             del self.trainedModel.data.P._dirPivot
             if self._A_is_affine >= 1 / 2 and self._b_is_affine >= 1 / 2:
                 self._affine_lvl.pop()
             del self._A_is_affine, self._b_is_affine
         self.trainedModel.data.npar = self.npar
 
     def errorEstimator(self, mus:Np1D, return_max : bool = False) -> Np1D:
         """Standard residual-based error estimator."""
         setupOK = self.setupApproxLocal()
         if setupOK > 0:
             err = np.empty(len(mus))
             err[:] = np.nan
             if not return_max: return err
             return err, [- setupOK], np.nan
         self._marginalizeTrainedModel(True)
         errRes = super().errorEstimator(mus, return_max)
         self._marginalizeTrainedModel(False)
         return errRes
 
     def _preliminaryTraining(self):
         """Initialize starting snapshots of solution map."""
         RROMPyAssert(self._mode, message = "Cannot start greedy algorithm.")
         self._S = self._setSampleBatch(self.S)
         self.resetSamples()
         self.samplingEngine.scaleFactor = self.scaleFactorDer
         musPivot = self.samplerTrainSet.generatePoints(self.S)
         while len(musPivot) > self.S: musPivot.pop()
         muTestPivot = self.samplerPivot.generatePoints(self.nTestPoints, False)
         idxPop = pruneSamples(self.mapParameterListPivot(muTestPivot),
                               self.mapParameterListPivot(musPivot),
                               1e-10 * self.scaleFactorPivot[0])
         muTestPivot.pop(idxPop)
         self._mus = emptyParameterList()
         self.mus.reset((self.S - 1, self.HFEngine.npar))
         self.muTest = emptyParameterList()
         self.muTest.reset((len(muTestPivot) + 1, self.HFEngine.npar))
         self.mus.data[:, self.directionPivot] = musPivot[: -1]
         self.mus.data[:, self.directionMarginal] = np.repeat(self.muMargLoc,
                                                           self.S - 1, axis = 0)
         self.muTest.data[: -1, self.directionPivot] = muTestPivot.data
         self.muTest.data[-1, self.directionPivot] = musPivot[-1]
         self.muTest.data[:, self.directionMarginal] = np.repeat(self.muMargLoc,
                                                           len(muTestPivot) + 1,
                                                           axis = 0)
         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 setupApproxLocal(self) -> int:
         """Compute rational interpolant."""
         self._marginalizeMiscellanea(True)
         setupOK = super().setupApproxLocal()
         self._marginalizeMiscellanea(False)
         return setupOK
 
     def setupApprox(self, *args, **kwargs) -> int:
         """Compute rational interpolant."""
         if self.checkComputedApprox(): return -1
         RROMPyAssert(self._mode, message = "Cannot setup approximant.")
         vbMng(self, "INIT", "Setting up {}.". format(self.name()), 5)
         self.computeScaleFactor()
         self._musMarginal = self.samplerMarginal.generatePoints(self.SMarginal)
         while len(self.musMarginal) > self.SMarginal: self.musMarginal.pop()
         S0 = copy(self.S)
         idx, sizes = indicesScatter(len(self.musMarginal), return_sizes = True)
         pMat, Ps, Qs, mus = None, [], [], None
         req, emptyCores = [], np.where(sizes == 0)[0]
         if len(idx) == 0:
             vbMng(self, "MAIN", "Idling.", 25)
             if self.storeAllSamples: self.storeSamples()
             pL, pT, mT = recv(source = 0, tag = poolRank())
             pMat = np.empty((pL, 0), dtype = pT)
             mus = np.empty((0, self.mu0.shape[1]), dtype = mT)
         else:
             _scaleFactorOldPivot = copy(self.scaleFactor)
             self.scaleFactor = self.scaleFactorPivot
             self._temporaryPivot = 1
             for i in idx:
                 self.muMargLoc = self.musMarginal[[i]]
                 vbMng(self, "MAIN",
                       "Building marginal model no. {} at {}.".format(i + 1,
                                                             self.muMargLoc), 5)
                 self.samplingEngine.resetHistory()
                 self.trainedModel = None
                 self.verbosity -= 5
                 self.samplingEngine.verbosity -= 10
                 RationalInterpolantGreedy.setupApprox(self, *args, **kwargs)
                 self.verbosity += 5
                 self.samplingEngine.verbosity += 10
                 if self.storeAllSamples: self.storeSamples(i)
                 musi = self.samplingEngine.mus
                 pMati = self.samplingEngine.projectionMatrix
                 if not hasattr(self, "matchState") or not self.matchState:
                     if self.POD == 1 and not (
                         hasattr(self.HFEngine.C, "is_mu_independent")
                     and self.HFEngine.C.is_mu_independent in self._output_lvl):
                         raise RROMPyException(("Cannot apply mu-dependent C "
                                                "to orthonormalized samples."))
                     vbMng(self, "INIT", "Extracting system output from state.",
                           35)
                     pMatiEff = None
                     for j, mu in enumerate(musi):
                         pMij = np.expand_dims(self.HFEngine.applyC(
                                                           pMati[:, j], mu), -1)
                         if pMatiEff is None:
                             pMatiEff = np.array(pMij)
                         else:
                             pMatiEff = np.append(pMatiEff, pMij, axis = 1)
                     pMati = pMatiEff
                     vbMng(self, "DEL", "Done extracting system output.", 35)
                 if pMat is None:
                     mus = copy(musi.data)
                     pMat = copy(pMati)
                     if i == 0:
                         for dest in emptyCores:
                             req += [isend((len(pMat), pMat.dtype, mus.dtype),
                                           dest = dest, tag = dest)]
                 else:
                     mus = np.vstack((mus, musi.data))
                     pMat = np.hstack((pMat, pMati))
                 Ps += [copy(self.trainedModel.data.P)]
                 Qs += [copy(self.trainedModel.data.Q)]
                 self._S = S0
             del self._temporaryPivot, self.muMargLoc
             self.scaleFactor = _scaleFactorOldPivot
         for r in req: r.wait()
         pMat, Ps, Qs, mus, nsamples = gatherPivotedApproximant(pMat, Ps, Qs,
                                                                mus, sizes,
                                                                self.polybasis)
         self._mus = self.checkParameterList(mus)
         Psupp = np.append(0, np.cumsum(nsamples))
         self._setupTrainedModel(pMat, forceNew = True)
         self.trainedModel.data.Qs, self.trainedModel.data.Ps = Qs, Ps
         self.trainedModel.data.Psupp = list(Psupp[: -1])
         self._preliminaryMarginalFinalization()
         self._finalizeMarginalization()
         vbMng(self, "DEL", "Done setting up approximant.", 5)
         return 0
 
 class RationalInterpolantGreedyPivotedNoMatch(
                                           RationalInterpolantGreedyPivotedBase,
                                           GenericPivotedApproximantNoMatch):
     """
     ROM pivoted rational interpolant (without pole matching) computation for
         parametric problems.
 
     Args:
         HFEngine: HF problem solver.
         mu0(optional): Default parameter. Defaults to 0.
         directionPivot(optional): Pivot components. Defaults to [0].
         approxParameters(optional): Dictionary containing values for main
             parameters of approximant. Recognized keys are:
             - 'POD': kind of snapshots orthogonalization; allowed values
                 include 0, 1/2, and 1; defaults to 1, i.e. POD;
             - 'scaleFactorDer': scaling factors for derivative computation;
                 defaults to 'AUTO';
             - 'S': total number of pivot samples current approximant relies
                 upon;
             - 'samplerPivot': pivot sample point generator;
             - 'SMarginal': total number of marginal samples current approximant
                 relies upon;
             - 'samplerMarginal': marginal sample point generator;
             - 'polybasis': type of polynomial basis for pivot
                 interpolation; defaults to 'MONOMIAL';
             - 'greedyTol': uniform error tolerance for greedy algorithm;
                 defaults to 1e-2;
             - 'collinearityTol': collinearity tolerance for greedy algorithm;
                 defaults to 0.;
             - 'maxIter': maximum number of greedy steps; defaults to 1e2;
             - 'nTestPoints': number of test points; defaults to 5e2;
             - 'samplerTrainSet': training sample points generator; defaults to
                 samplerPivot;
             - 'errorEstimatorKind': kind of error estimator; available values
                 include 'AFFINE', 'DISCREPANCY', 'LOOK_AHEAD',
                 'LOOK_AHEAD_RES', and 'NONE'; defaults to 'NONE';
             - 'radialDirectionalWeightsMarginal': radial basis weights for
                 marginal interpolant; defaults to 1;
             - 'functionalSolve': strategy for minimization of denominator
-                functional; allowed values include 'NORM', 'DOMINANT', 'NODAL',
-                'BARYCENTRIC_NORM', and 'BARYCENTRIC[_AVERAGE]' (check pdf in
-                main folder for meaning); defaults to 'NORM';
+                functional; allowed values include 'NORM', 'DOMINANT',
+                'BARYCENTRIC[_NORM]', and 'BARYCENTRIC_AVERAGE' (check pdf in
+                main folder for explanation); defaults to 'NORM';
             - 'interpTol': tolerance for pivot interpolation; defaults to
                 None;
             - 'QTol': tolerance for robust rational denominator management;
                 defaults to 0.
             Defaults to empty dict.
         verbosity(optional): Verbosity level. Defaults to 10.
             
     Attributes:
         HFEngine: HF problem solver.
         mu0: Default parameter.
         directionPivot: Pivot components.
         mus: Array of snapshot parameters.
         musMarginal: Array of marginal snapshot parameters.
         approxParameters: Dictionary containing values for main parameters of
             approximant. Recognized keys are in parameterList.
         parameterListSoft: Recognized keys of soft approximant parameters:
             - 'POD': kind of snapshots orthogonalization;
             - 'scaleFactorDer': scaling factors for derivative computation;
             - 'polybasis': type of polynomial basis for pivot
                 interpolation;
             - 'greedyTol': uniform error tolerance for greedy algorithm;
             - 'collinearityTol': collinearity tolerance for greedy algorithm;
             - 'maxIter': maximum number of greedy steps;
             - 'nTestPoints': number of test points;
             - 'samplerTrainSet': training sample points generator;
             - 'errorEstimatorKind': kind of error estimator;
             - 'radialDirectionalWeightsMarginal': radial basis weights for
                 marginal interpolant;
             - 'functionalSolve': strategy for minimization of denominator
                 functional;
             - 'interpTol': tolerance for pivot interpolation;
             - 'QTol': tolerance for robust rational denominator management.
         parameterListCritical: Recognized keys of critical approximant
             parameters:
             - 'S': total number of pivot samples current approximant relies
                 upon;
             - 'samplerPivot': pivot sample point generator;
             - 'SMarginal': total number of marginal samples current approximant
                 relies upon;
             - 'samplerMarginal': marginal sample point generator.
         verbosity: Verbosity level.
         POD: Kind of snapshots orthogonalization.
         scaleFactorDer: Scaling factors for derivative computation.
         S: Total number of pivot samples current approximant relies upon.
         samplerPivot: Pivot sample point generator.
         SMarginal: Total number of marginal samples current approximant relies
             upon.
         samplerMarginal: Marginal sample point generator.
         polybasis: Type of polynomial basis for pivot interpolation.
         greedyTol: uniform error tolerance for greedy algorithm.
         collinearityTol: Collinearity tolerance for greedy algorithm.
         maxIter: maximum number of greedy steps.
         nTestPoints: number of starting training points.
         samplerTrainSet: training sample points generator.
         errorEstimatorKind: kind of error estimator.
         radialDirectionalWeightsMarginal: Radial basis weights for marginal
             interpolant.
         functionalSolve: Strategy for minimization of denominator functional.
         interpTol: Tolerance for pivot interpolation.
         QTol: Tolerance for robust rational denominator management.
         muBounds: list of bounds for pivot parameter values.
         muBoundsMarginal: list of bounds for marginal parameter values.
         samplingEngine: Sampling engine.
         uHF: High fidelity solution(s) with parameter(s) lastSolvedHF as
             sampleList.
         lastSolvedHF: Parameter(s) corresponding to last computed high fidelity
             solution(s) as parameterList.
         uApproxReduced: Reduced approximate solution(s) with parameter(s)
             lastSolvedApprox as sampleList.
         lastSolvedApproxReduced: Parameter(s) corresponding to last computed
             reduced approximate solution(s) as parameterList.
         uApprox: Approximate solution(s) with parameter(s) lastSolvedApprox as
             sampleList.
         lastSolvedApprox: Parameter(s) corresponding to last computed
             approximate solution(s) as parameterList.
         Q: Numpy 1D vector containing complex coefficients of approximant
             denominator.
         P: Numpy 2D vector whose columns are FE dofs of coefficients of
             approximant numerator.
     """
 
 class RationalInterpolantGreedyPivotedPoleMatch(
                                           RationalInterpolantGreedyPivotedBase,
                                           GenericPivotedApproximantPoleMatch):
     """
     ROM pivoted rational interpolant (with pole matching) computation for
         parametric problems.
 
     Args:
         HFEngine: HF problem solver.
         mu0(optional): Default parameter. Defaults to 0.
         directionPivot(optional): Pivot components. Defaults to [0].
         approxParameters(optional): Dictionary containing values for main
             parameters of approximant. Recognized keys are:
             - 'POD': kind of snapshots orthogonalization; allowed values
                 include 0, 1/2, and 1; defaults to 1, i.e. POD;
             - 'scaleFactorDer': scaling factors for derivative computation;
                 defaults to 'AUTO';
             - 'matchState': whether to match the system state rather than the
                 system output; defaults to False;
             - 'matchingWeight': weight for pole matching optimization; defaults
                 to 1;
             - 'matchingChordalRadius': radius to be used in chordal metric for
                 poles and residues; if <= 0, Euclidean metric is used; if
                 'AUTO', automatically selected; defaults to -1;
             - 'matchingShared': required ratio of marginal points to share
                 resonance; defaults to 1.;
             - 'badPoleCorrection': strategy for correction of bad poles;
                 available values include 'ERASE', 'RATIONAL', and 'POLYNOMIAL';
                 defaults to 'ERASE';
             - 'S': total number of pivot samples current approximant relies
                 upon;
             - 'samplerPivot': pivot sample point generator;
             - 'SMarginal': total number of marginal samples current approximant
                 relies upon;
             - 'samplerMarginal': marginal sample point generator;
             - 'polybasis': type of polynomial basis for pivot
                 interpolation; defaults to 'MONOMIAL';
             - 'polybasisMarginal': type of polynomial basis for marginal
                 interpolation; allowed values include 'MONOMIAL_*',
                 'CHEBYSHEV_*', 'LEGENDRE_*', 'NEARESTNEIGHBOR', and 
                 'PIECEWISE_LINEAR_*'; defaults to 'MONOMIAL';
             - 'paramsMarginal': dictionary of parameters for marginal
                 interpolation; include:
                 . 'MMarginal': degree of marginal interpolant; defaults to
                     'AUTO', i.e. maximum allowed; not for 'NEARESTNEIGHBOR' or
                     'PIECEWISE_LINEAR_*';
                 . 'nNeighborsMarginal': number of marginal nearest neighbors;
                      defaults to 1; only for 'NEARESTNEIGHBOR';
                 . 'polydegreetypeMarginal': type of polynomial degree for
                     marginal; defaults to 'TOTAL'; not for 'NEARESTNEIGHBOR' or
                     'PIECEWISE_LINEAR_*';
                 . 'interpTolMarginal': tolerance for marginal interpolation;
                     defaults to None; not for 'NEARESTNEIGHBOR' or
                     'PIECEWISE_LINEAR_*';
                 . 'radialDirectionalWeightsMarginalAdapt': bounds for adaptive
                     rescaling of marginal radial basis weights; only for
                     radial basis.
             - 'greedyTol': uniform error tolerance for greedy algorithm;
                 defaults to 1e-2;
             - 'collinearityTol': collinearity tolerance for greedy algorithm;
                 defaults to 0.;
             - 'maxIter': maximum number of greedy steps; defaults to 1e2;
             - 'nTestPoints': number of test points; defaults to 5e2;
             - 'samplerTrainSet': training sample points generator; defaults to
                 samplerPivot;
             - 'errorEstimatorKind': kind of error estimator; available values
                 include 'AFFINE', 'DISCREPANCY', 'LOOK_AHEAD',
                 'LOOK_AHEAD_RES', and 'NONE'; defaults to 'NONE';
             - 'radialDirectionalWeightsMarginal': radial basis weights for
                 marginal interpolant; defaults to 1;
             - 'functionalSolve': strategy for minimization of denominator
-                functional; allowed values include 'NORM', 'DOMINANT', 'NODAL',
-                'BARYCENTRIC_NORM', and 'BARYCENTRIC[_AVERAGE]' (check pdf in
-                main folder for meaning); defaults to 'NORM';
+                functional; allowed values include 'NORM', 'DOMINANT',
+                'BARYCENTRIC[_NORM]', and 'BARYCENTRIC_AVERAGE' (check pdf in
+                main folder for explanation); defaults to 'NORM';
             - 'interpTol': tolerance for pivot interpolation; defaults to None;
             - 'QTol': tolerance for robust rational denominator management;
                 defaults to 0.
             Defaults to empty dict.
         verbosity(optional): Verbosity level. Defaults to 10.
             
     Attributes:
         HFEngine: HF problem solver.
         mu0: Default parameter.
         directionPivot: Pivot components.
         mus: Array of snapshot parameters.
         musMarginal: Array of marginal snapshot parameters.
         approxParameters: Dictionary containing values for main parameters of
             approximant. Recognized keys are in parameterList.
         parameterListSoft: Recognized keys of soft approximant parameters:
             - 'POD': kind of snapshots orthogonalization;
             - 'scaleFactorDer': scaling factors for derivative computation;
             - 'matchState': whether to match the system state rather than the
                 system output;
             - 'matchingWeight': weight for pole matching optimization;
             - 'matchingChordalRadius': radius to be used in chordal metric for
                 poles and residues;
             - 'matchingShared': required ratio of marginal points to share
                 resonance;
             - 'badPoleCorrection': strategy for correction of bad poles;
             - 'polybasis': type of polynomial basis for pivot
                 interpolation;
             - 'polybasisMarginal': type of polynomial basis for marginal
                 interpolation;
             - 'paramsMarginal': dictionary of parameters for marginal
                 interpolation; include:
                 . 'MMarginal': degree of marginal interpolant;
                 . 'nNeighborsMarginal': number of marginal nearest neighbors;
                 . 'polydegreetypeMarginal': type of polynomial degree for
                     marginal;
                 . 'interpTolMarginal': tolerance for marginal interpolation;
                 . 'radialDirectionalWeightsMarginalAdapt': bounds for adaptive
                     rescaling of marginal radial basis weights.
             - 'greedyTol': uniform error tolerance for greedy algorithm;
             - 'collinearityTol': collinearity tolerance for greedy algorithm;
             - 'maxIter': maximum number of greedy steps;
             - 'nTestPoints': number of test points;
             - 'samplerTrainSet': training sample points generator;
             - 'errorEstimatorKind': kind of error estimator;
             - 'radialDirectionalWeightsMarginal': radial basis weights for
                 marginal interpolant;
             - 'functionalSolve': strategy for minimization of denominator
                 functional;
             - 'interpTol': tolerance for pivot interpolation;
             - 'QTol': tolerance for robust rational denominator management.
         parameterListCritical: Recognized keys of critical approximant
             parameters:
             - 'S': total number of pivot samples current approximant relies
                 upon;
             - 'samplerPivot': pivot sample point generator;
             - 'SMarginal': total number of marginal samples current approximant
                 relies upon;
             - 'samplerMarginal': marginal sample point generator.
         verbosity: Verbosity level.
         POD: Kind of snapshots orthogonalization.
         scaleFactorDer: Scaling factors for derivative computation.
         matchState: Whether to match the system state rather than the system
             output.
         matchingWeight: Weight for pole matching optimization.
         matchingChordalRadius: Radius to be used in chordal metric for poles
             and residues.
         matchingShared: Required ratio of marginal points to share resonance.
         badPoleCorrection: Strategy for correction of bad poles.
         S: Total number of pivot samples current approximant relies upon.
         samplerPivot: Pivot sample point generator.
         SMarginal: Total number of marginal samples current approximant relies
             upon.
         samplerMarginal: Marginal sample point generator.
         polybasis: Type of polynomial basis for pivot interpolation.
         polybasisMarginal: Type of polynomial basis for marginal interpolation.
         paramsMarginal: Dictionary of parameters for marginal interpolation.
         greedyTol: uniform error tolerance for greedy algorithm.
         collinearityTol: Collinearity tolerance for greedy algorithm.
         maxIter: maximum number of greedy steps.
         nTestPoints: number of starting training points.
         samplerTrainSet: training sample points generator.
         errorEstimatorKind: kind of error estimator.
         radialDirectionalWeightsMarginal: Radial basis weights for marginal
             interpolant.
         functionalSolve: Strategy for minimization of denominator functional.
         interpTol: Tolerance for pivot interpolation.
         QTol: Tolerance for robust rational denominator management.
         muBounds: list of bounds for pivot parameter values.
         muBoundsMarginal: list of bounds for marginal parameter values.
         samplingEngine: Sampling engine.
         uHF: High fidelity solution(s) with parameter(s) lastSolvedHF as
             sampleList.
         lastSolvedHF: Parameter(s) corresponding to last computed high fidelity
             solution(s) as parameterList.
         uApproxReduced: Reduced approximate solution(s) with parameter(s)
             lastSolvedApprox as sampleList.
         lastSolvedApproxReduced: Parameter(s) corresponding to last computed
             reduced approximate solution(s) as parameterList.
         uApprox: Approximate solution(s) with parameter(s) lastSolvedApprox as
             sampleList.
         lastSolvedApprox: Parameter(s) corresponding to last computed
             approximate solution(s) as parameterList.
         Q: Numpy 1D vector containing complex coefficients of approximant
             denominator.
         P: Numpy 2D vector whose columns are FE dofs of coefficients of
             approximant numerator.
     """
 
     def setupApprox(self, *args, **kwargs) -> int:
         if self.checkComputedApprox(): return -1
         self.purgeparamsMarginal()
         setupOK = super().setupApprox(*args, **kwargs)
         if self.matchState: self._postApplyC()
         return setupOK
diff --git a/rrompy/reduction_methods/pivoted/rational_interpolant_pivoted.py b/rrompy/reduction_methods/pivoted/rational_interpolant_pivoted.py
index 816e2ca..bb07351 100644
--- a/rrompy/reduction_methods/pivoted/rational_interpolant_pivoted.py
+++ b/rrompy/reduction_methods/pivoted/rational_interpolant_pivoted.py
@@ -1,491 +1,491 @@
 # Copyright (C) 2018-2020 by the RROMPy authors
 #
 # This file is part of RROMPy.
 #
 # RROMPy is free software: you can redistribute it and/or modify
 # it under the terms of the GNU Lesser General Public License as published by
 # the Free Software Foundation, either version 3 of the License, or
 # (at your option) any later version.
 #
 # RROMPy is distributed in the hope that it will be useful,
 # but WITHOUT ANY WARRANTY; without even the implied warranty of
 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
 # GNU Lesser General Public License for more details.
 #
 # You should have received a copy of the GNU Lesser General Public License
 # along with RROMPy. If not, see <http://www.gnu.org/licenses/>.
 #
 
 import numpy as np
 from collections.abc import Iterable
 from copy import deepcopy as copy
 from .generic_pivoted_approximant import (GenericPivotedApproximantBase,
                                           GenericPivotedApproximantNoMatch,
                                           GenericPivotedApproximantPoleMatch)
 from .gather_pivoted_approximant import gatherPivotedApproximant
 from rrompy.reduction_methods.standard.rational_interpolant import (
                                                            RationalInterpolant)
 from rrompy.utilities.base import verbosityManager as vbMng
 from rrompy.utilities.numerical.hash_derivative import nextDerivativeIndices
 from rrompy.utilities.exception_manager import (RROMPyException, RROMPyAssert,
                                                 RROMPyWarning)
 from rrompy.parameter import emptyParameterList
 from rrompy.utilities.parallel import poolRank, indicesScatter, isend, recv
 
 __all__ = ['RationalInterpolantPivotedNoMatch',
            'RationalInterpolantPivotedPoleMatch']
 
 class RationalInterpolantPivotedBase(GenericPivotedApproximantBase,
                                      RationalInterpolant):
     def __init__(self, *args, **kwargs):
         self._preInit()
         self._addParametersToList(toBeExcluded = ["polydegreetype"])
         super().__init__(*args, **kwargs)
         if self.nparPivot > 1: self.HFEngine._ignoreResidues = 1
         self._postInit()
 
     @property
     def scaleFactorDer(self):
         """Value of scaleFactorDer."""
         if self._scaleFactorDer == "NONE": return 1.
         if self._scaleFactorDer == "AUTO": return self.scaleFactorPivot
         return self._scaleFactorDer
     @scaleFactorDer.setter
     def scaleFactorDer(self, scaleFactorDer):
         if isinstance(scaleFactorDer, (str,)):
             scaleFactorDer = scaleFactorDer.upper()
         elif isinstance(scaleFactorDer, Iterable):
             scaleFactorDer = list(scaleFactorDer)
         self._scaleFactorDer = scaleFactorDer
         self._approxParameters["scaleFactorDer"] = self._scaleFactorDer
 
     @property
     def polydegreetype(self):
         """Value of polydegreetype."""
         return "TOTAL"
     @polydegreetype.setter
     def polydegreetype(self, polydegreetype):
         RROMPyWarning(("polydegreetype is used just to simplify inheritance, "
                        "and its value cannot be changed from 'TOTAL'."))
 
     def _setupInterpolationIndices(self):
         """Setup parameters for polyvander."""
         RROMPyAssert(self._mode,
                      message = "Cannot setup interpolation indices.")
         if (self._musUniqueCN is None 
          or len(self._reorder) != len(self.musPivot)):
             try:
                 muPC = self.trainedModel.centerNormalizePivot(self.musPivot)
             except:
                 muPC = self.trainedModel.centerNormalize(self.musPivot)
             self._musUniqueCN, musIdxsTo, musIdxs, musCount = (muPC.unique(
                                     return_index = True, return_inverse = True,
                                     return_counts = True))
             self._musUnique = self.musPivot[musIdxsTo]
             self._derIdxs = [None] * len(self._musUniqueCN)
             self._reorder = np.empty(len(musIdxs), dtype = int)
             filled = 0
             for j, cnt in enumerate(musCount):
                 self._derIdxs[j] = nextDerivativeIndices([], self.nparPivot,
                                                          cnt)
                 jIdx = np.nonzero(musIdxs == j)[0]
                 self._reorder[jIdx] = np.arange(filled, filled + cnt)
                 filled += cnt
 
     def setupApprox(self) -> int:
         """Compute rational interpolant."""
         if self.checkComputedApprox(): return -1
         RROMPyAssert(self._mode, message = "Cannot setup approximant.")
         vbMng(self, "INIT", "Setting up {}.". format(self.name()), 5)
         self.computeScaleFactor()
         self.resetSamples()
         self.samplingEngine.scaleFactor = self.scaleFactorDer
         self.musPivot = self.samplerPivot.generatePoints(self.S)
         while len(self.musPivot) > self.S: self.musPivot.pop()
         self._musMarginal = self.samplerMarginal.generatePoints(self.SMarginal)
         while len(self.musMarginal) > self.SMarginal: self.musMarginal.pop()
         self._mus = emptyParameterList()
         self.mus.reset((self.S * self.SMarginal, self.HFEngine.npar))
         self.mus.data[:, self.directionPivot] = np.tile(self.musPivot.data,
                                                         (self.SMarginal, 1))
         self.mus.data[:, self.directionMarginal] = np.repeat(
                                                          self.musMarginal.data,
                                                          self.S, axis = 0)
         N0 = copy(self.N)
         self._setupTrainedModel(np.zeros((0, 0)), forceNew = True)
         idx, sizes = indicesScatter(len(self.musMarginal), return_sizes = True)
         pMat, Ps, Qs = None, [], []
         req, emptyCores = [], np.where(sizes == 0)[0]
         if len(idx) == 0:
             vbMng(self, "MAIN", "Idling.", 30)
             if self.storeAllSamples: self.storeSamples()
             pL, pT = recv(source = 0, tag = poolRank())
             pMat = np.empty((pL, 0), dtype = pT)
         else:
             _scaleFactorOldPivot = copy(self.scaleFactor)
             self.scaleFactor = self.scaleFactorPivot
             self._temporaryPivot = 1
             for i in idx:
                 musi = self.mus[self.S * i : self.S * (i + 1)]
                 vbMng(self, "MAIN",
                       "Building marginal model no. {} at {}.".format(i + 1,
                                                        self.musMarginal[i]), 5)
                 vbMng(self, "INIT", "Starting computation of snapshots.", 10)
                 self.samplingEngine.resetHistory()
                 self.samplingEngine.iterSample(musi)
                 vbMng(self, "DEL", "Done computing snapshots.", 10)
                 self.verbosity -= 5
                 self.samplingEngine.verbosity -= 10
                 self._setupRational(self._setupDenominator())
                 self.verbosity += 5
                 self.samplingEngine.verbosity += 10
                 if self.storeAllSamples: self.storeSamples(i)
                 pMati = self.samplingEngine.projectionMatrix
                 if not hasattr(self, "matchState") or not self.matchState:
                     if self.POD == 1 and not (
                         hasattr(self.HFEngine.C, "is_mu_independent")
                     and self.HFEngine.C.is_mu_independent in self._output_lvl):
                         raise RROMPyException(("Cannot apply mu-dependent C "
                                                "to orthonormalized samples."))
                     vbMng(self, "INIT", "Extracting system output from state.",
                           35)
                     pMatiEff = None
                     for j, mu in enumerate(musi):
                         pMij = np.expand_dims(self.HFEngine.applyC(
                                                           pMati[:, j], mu), -1)
                         if pMatiEff is None:
                             pMatiEff = np.array(pMij)
                         else:
                             pMatiEff = np.append(pMatiEff, pMij, axis = 1)
                     pMati = pMatiEff
                     vbMng(self, "DEL", "Done extracting system output.", 35)
                 if pMat is None:
                     pMat = copy(pMati)
                     if i == 0:
                         for dest in emptyCores:
                             req += [isend((len(pMat), pMat.dtype), dest = dest,
                                           tag = dest)]
                 else:
                     pMat = np.hstack((pMat, pMati))
                 Ps += [copy(self.trainedModel.data.P)]
                 Qs += [copy(self.trainedModel.data.Q)]
                 del self.trainedModel.data.Q, self.trainedModel.data.P
                 self.N = N0
             del self._temporaryPivot
         self.scaleFactor = _scaleFactorOldPivot
         for r in req: r.wait()
         pMat, Ps, Qs, _, _ = gatherPivotedApproximant(pMat, Ps, Qs,
                                                       self.mus.data, sizes,
                                                       self.polybasis, False)
         self._setupTrainedModel(pMat)
         self.trainedModel.data.Qs, self.trainedModel.data.Ps = Qs, Ps
         Psupp = np.arange(0, len(self.musMarginal) * self.S, self.S)
         self.trainedModel.data.Psupp = list(Psupp)
         self._preliminaryMarginalFinalization()
         self._finalizeMarginalization()
         vbMng(self, "DEL", "Done setting up approximant.", 5)
         return 0
 
 class RationalInterpolantPivotedNoMatch(RationalInterpolantPivotedBase,
                                         GenericPivotedApproximantNoMatch):
     """
     ROM pivoted rational interpolant (without pole matching) computation for
         parametric problems.
 
     Args:
         HFEngine: HF problem solver.
         mu0(optional): Default parameter. Defaults to 0.
         directionPivot(optional): Pivot components. Defaults to [0].
         approxParameters(optional): Dictionary containing values for main
             parameters of approximant. Recognized keys are:
             - 'POD': kind of snapshots orthogonalization; allowed values
                 include 0, 1/2, and 1; defaults to 1, i.e. POD;
             - 'scaleFactorDer': scaling factors for derivative computation;
                 defaults to 'AUTO';
             - 'S': total number of pivot samples current approximant relies
                 upon;
             - 'samplerPivot': pivot sample point generator;
             - 'SMarginal': total number of marginal samples current approximant
                 relies upon;
             - 'samplerMarginal': marginal sample point generator;
             - 'polybasis': type of polynomial basis for pivot
                 interpolation; defaults to 'MONOMIAL';
             - 'M': degree of rational interpolant numerator; defaults to
                 'AUTO', i.e. maximum allowed;
             - 'N': degree of rational interpolant denominator; defaults to
                 'AUTO', i.e. maximum allowed;
             - 'radialDirectionalWeights': radial basis weights for pivot
                 numerator; defaults to 1;
             - 'radialDirectionalWeightsAdapt': bounds for adaptive rescaling of
                 radial basis weights; defaults to [-1, -1];
             - 'radialDirectionalWeightsMarginal': radial basis weights for
                 marginal interpolant; defaults to 1;
             - 'functionalSolve': strategy for minimization of denominator
-                functional; allowed values include 'NORM', 'DOMINANT', 'NODAL',
-                'BARYCENTRIC_NORM', and 'BARYCENTRIC[_AVERAGE]' (check pdf in
-                main folder for meaning); defaults to 'NORM';
+                functional; allowed values include 'NORM', 'DOMINANT',
+                'BARYCENTRIC[_NORM]', and 'BARYCENTRIC_AVERAGE' (check pdf in
+                main folder for explanation); defaults to 'NORM';
             - 'interpTol': tolerance for pivot interpolation; defaults to
                 None;
             - 'QTol': tolerance for robust rational denominator management;
                 defaults to 0.
             Defaults to empty dict.
         verbosity(optional): Verbosity level. Defaults to 10.
             
     Attributes:
         HFEngine: HF problem solver.
         mu0: Default parameter.
         directionPivot: Pivot components.
         mus: Array of snapshot parameters.
         musPivot: Array of pivot snapshot parameters.
         musMarginal: Array of marginal snapshot parameters.
         approxParameters: Dictionary containing values for main parameters of
             approximant. Recognized keys are in parameterList.
         parameterListSoft: Recognized keys of soft approximant parameters:
             - 'POD': kind of snapshots orthogonalization;
             - 'scaleFactorDer': scaling factors for derivative computation;
             - 'polybasis': type of polynomial basis for pivot
                 interpolation;
             - 'M': degree of rational interpolant numerator;
             - 'N': degree of rational interpolant denominator;
             - 'radialDirectionalWeights': radial basis weights for pivot
                 numerator;
             - 'radialDirectionalWeightsAdapt': bounds for adaptive rescaling of
                 radial basis weights;
             - 'radialDirectionalWeightsMarginal': radial basis weights for
                 marginal interpolant;
             - 'functionalSolve': strategy for minimization of denominator
                 functional;
             - 'interpTol': tolerance for pivot interpolation;
             - 'QTol': tolerance for robust rational denominator management.
         parameterListCritical: Recognized keys of critical approximant
             parameters:
             - 'S': total number of pivot samples current approximant relies
                 upon;
             - 'samplerPivot': pivot sample point generator;
             - 'SMarginal': total number of marginal samples current approximant
                 relies upon;
             - 'samplerMarginal': marginal sample point generator.
         verbosity: Verbosity level.
         POD: Kind of snapshots orthogonalization.
         scaleFactorDer: Scaling factors for derivative computation.
         S: Total number of pivot samples current approximant relies upon.
         samplerPivot: Pivot sample point generator.
         SMarginal: Total number of marginal samples current approximant relies
             upon.
         samplerMarginal: Marginal sample point generator.
         polybasis: Type of polynomial basis for pivot interpolation.
         M: Numerator degree of approximant.
         N: Denominator degree of approximant.
         radialDirectionalWeights: Radial basis weights for pivot numerator.
         radialDirectionalWeightsAdapt: Bounds for adaptive rescaling of radial
             basis weights.
         radialDirectionalWeightsMarginal: Radial basis weights for marginal
             interpolant.
         functionalSolve: Strategy for minimization of denominator functional.
         interpTol: Tolerance for pivot interpolation.
         QTol: Tolerance for robust rational denominator management.
         muBounds: list of bounds for pivot parameter values.
         muBoundsMarginal: list of bounds for marginal parameter values.
         samplingEngine: Sampling engine.
         uHF: High fidelity solution(s) with parameter(s) lastSolvedHF as
             sampleList.
         lastSolvedHF: Parameter(s) corresponding to last computed high fidelity
             solution(s) as parameterList.
         uApproxReduced: Reduced approximate solution(s) with parameter(s)
             lastSolvedApprox as sampleList.
         lastSolvedApproxReduced: Parameter(s) corresponding to last computed
             reduced approximate solution(s) as parameterList.
         uApprox: Approximate solution(s) with parameter(s) lastSolvedApprox as
             sampleList.
         lastSolvedApprox: Parameter(s) corresponding to last computed
             approximate solution(s) as parameterList.
         Q: Numpy 1D vector containing complex coefficients of approximant
             denominator.
         P: Numpy 2D vector whose columns are FE dofs of coefficients of
             approximant numerator.
     """
     
 class RationalInterpolantPivotedPoleMatch(RationalInterpolantPivotedBase,
                                           GenericPivotedApproximantPoleMatch):
     """
     ROM pivoted rational interpolant (with pole matching) computation for
         parametric problems.
 
     Args:
         HFEngine: HF problem solver.
         mu0(optional): Default parameter. Defaults to 0.
         directionPivot(optional): Pivot components. Defaults to [0].
         approxParameters(optional): Dictionary containing values for main
             parameters of approximant. Recognized keys are:
             - 'POD': kind of snapshots orthogonalization; allowed values
                 include 0, 1/2, and 1; defaults to 1, i.e. POD;
             - 'scaleFactorDer': scaling factors for derivative computation;
                 defaults to 'AUTO';
             - 'matchState': whether to match the system state rather than the
                 system output; defaults to False;
             - 'matchingWeight': weight for pole matching optimization; defaults
                 to 1;
             - 'matchingChordalRadius': radius to be used in chordal metric for
                 poles and residues; if <= 0, Euclidean metric is used; if
                 'AUTO', automatically selected; defaults to -1;
             - 'matchingShared': required ratio of marginal points to share
                 resonance; defaults to 1.;
             - 'badPoleCorrection': strategy for correction of bad poles;
                 available values include 'ERASE', 'RATIONAL', and 'POLYNOMIAL';
                 defaults to 'ERASE';
             - 'S': total number of pivot samples current approximant relies
                 upon;
             - 'samplerPivot': pivot sample point generator;
             - 'SMarginal': total number of marginal samples current approximant
                 relies upon;
             - 'samplerMarginal': marginal sample point generator;
             - 'polybasis': type of polynomial basis for pivot
                 interpolation; defaults to 'MONOMIAL';
             - 'polybasisMarginal': type of polynomial basis for marginal
                 interpolation; allowed values include 'MONOMIAL_*',
                 'CHEBYSHEV_*', 'LEGENDRE_*', 'NEARESTNEIGHBOR', and 
                 'PIECEWISE_LINEAR_*'; defaults to 'MONOMIAL';
             - 'paramsMarginal': dictionary of parameters for marginal
                 interpolation; include:
                 . 'MMarginal': degree of marginal interpolant; defaults to
                     'AUTO', i.e. maximum allowed; not for 'NEARESTNEIGHBOR' or
                     'PIECEWISE_LINEAR_*';
                 . 'nNeighborsMarginal': number of marginal nearest neighbors;
                      defaults to 1; only for 'NEARESTNEIGHBOR';
                 . 'polydegreetypeMarginal': type of polynomial degree for
                     marginal; defaults to 'TOTAL'; not for 'NEARESTNEIGHBOR' or
                     'PIECEWISE_LINEAR_*';
                 . 'interpTolMarginal': tolerance for marginal interpolation;
                     defaults to None; not for 'NEARESTNEIGHBOR' or
                     'PIECEWISE_LINEAR_*';
                 . 'radialDirectionalWeightsMarginalAdapt': bounds for adaptive
                     rescaling of marginal radial basis weights; only for
                     radial basis.
             - 'M': degree of rational interpolant numerator; defaults to
                 'AUTO', i.e. maximum allowed;
             - 'N': degree of rational interpolant denominator; defaults to
                 'AUTO', i.e. maximum allowed;
             - 'radialDirectionalWeights': radial basis weights for pivot
                 numerator; defaults to 1;
             - 'radialDirectionalWeightsAdapt': bounds for adaptive rescaling of
                 radial basis weights; defaults to [-1, -1];
             - 'radialDirectionalWeightsMarginal': radial basis weights for
                 marginal interpolant; defaults to 1;
             - 'functionalSolve': strategy for minimization of denominator
-                functional; allowed values include 'NORM', 'DOMINANT', 'NODAL',
-                'BARYCENTRIC_NORM', and 'BARYCENTRIC[_AVERAGE]' (check pdf in
-                main folder for meaning); defaults to 'NORM';
+                functional; allowed values include 'NORM', 'DOMINANT',
+                'BARYCENTRIC[_NORM]', and 'BARYCENTRIC_AVERAGE' (check pdf in
+                main folder for explanation); defaults to 'NORM';
             - 'interpTol': tolerance for pivot interpolation; defaults to None;
             - 'QTol': tolerance for robust rational denominator management;
                 defaults to 0.
             Defaults to empty dict.
         verbosity(optional): Verbosity level. Defaults to 10.
             
     Attributes:
         HFEngine: HF problem solver.
         mu0: Default parameter.
         directionPivot: Pivot components.
         mus: Array of snapshot parameters.
         musPivot: Array of pivot snapshot parameters.
         musMarginal: Array of marginal snapshot parameters.
         approxParameters: Dictionary containing values for main parameters of
             approximant. Recognized keys are in parameterList.
         parameterListSoft: Recognized keys of soft approximant parameters:
             - 'POD': kind of snapshots orthogonalization;
             - 'scaleFactorDer': scaling factors for derivative computation;
             - 'matchState': whether to match the system state rather than the
                 system output;
             - 'matchingWeight': weight for pole matching optimization;
             - 'matchingChordalRadius': radius to be used in chordal metric for
                 poles and residues;
             - 'matchingShared': required ratio of marginal points to share
                 resonance;
             - 'badPoleCorrection': strategy for correction of bad poles;
             - 'polybasis': type of polynomial basis for pivot
                 interpolation;
             - 'polybasisMarginal': type of polynomial basis for marginal
                 interpolation;
             - 'paramsMarginal': dictionary of parameters for marginal
                 interpolation; include:
                 . 'MMarginal': degree of marginal interpolant;
                 . 'nNeighborsMarginal': number of marginal nearest neighbors;
                 . 'polydegreetypeMarginal': type of polynomial degree for
                     marginal;
                 . 'interpTolMarginal': tolerance for marginal interpolation;
                 . 'radialDirectionalWeightsMarginalAdapt': bounds for adaptive
                     rescaling of marginal radial basis weights.
             - 'M': degree of rational interpolant numerator;
             - 'N': degree of rational interpolant denominator;
             - 'radialDirectionalWeights': radial basis weights for pivot
                 numerator;
             - 'radialDirectionalWeightsAdapt': bounds for adaptive rescaling of
                 radial basis weights;
             - 'radialDirectionalWeightsMarginal': radial basis weights for
                 marginal interpolant;
             - 'functionalSolve': strategy for minimization of denominator
                 functional;
             - 'interpTol': tolerance for pivot interpolation;
             - 'QTol': tolerance for robust rational denominator management.
         parameterListCritical: Recognized keys of critical approximant
             parameters:
             - 'S': total number of pivot samples current approximant relies
                 upon;
             - 'samplerPivot': pivot sample point generator;
             - 'SMarginal': total number of marginal samples current approximant
                 relies upon;
             - 'samplerMarginal': marginal sample point generator.
         verbosity: Verbosity level.
         POD: Kind of snapshots orthogonalization.
         scaleFactorDer: Scaling factors for derivative computation.
         matchState: Whether to match the system state rather than the system
             output.
         matchingWeight: Weight for pole matching optimization.
         matchingChordalRadius: Radius to be used in chordal metric for poles
             and residues.
         matchingShared: Required ratio of marginal points to share resonance.
         badPoleCorrection: Strategy for correction of bad poles.
         S: Total number of pivot samples current approximant relies upon.
         samplerPivot: Pivot sample point generator.
         SMarginal: Total number of marginal samples current approximant relies
             upon.
         samplerMarginal: Marginal sample point generator.
         polybasis: Type of polynomial basis for pivot interpolation.
         polybasisMarginal: Type of polynomial basis for marginal interpolation.
         paramsMarginal: Dictionary of parameters for marginal interpolation.
         M: Numerator degree of approximant.
         N: Denominator degree of approximant.
         radialDirectionalWeights: Radial basis weights for pivot numerator.
         radialDirectionalWeightsAdapt: Bounds for adaptive rescaling of radial
             basis weights.
         radialDirectionalWeightsMarginal: Radial basis weights for marginal
             interpolant.
         functionalSolve: Strategy for minimization of denominator functional.
         interpTol: Tolerance for pivot interpolation.
         QTol: Tolerance for robust rational denominator management.
         muBounds: list of bounds for pivot parameter values.
         muBoundsMarginal: list of bounds for marginal parameter values.
         samplingEngine: Sampling engine.
         uHF: High fidelity solution(s) with parameter(s) lastSolvedHF as
             sampleList.
         lastSolvedHF: Parameter(s) corresponding to last computed high fidelity
             solution(s) as parameterList.
         uApproxReduced: Reduced approximate solution(s) with parameter(s)
             lastSolvedApprox as sampleList.
         lastSolvedApproxReduced: Parameter(s) corresponding to last computed
             reduced approximate solution(s) as parameterList.
         uApprox: Approximate solution(s) with parameter(s) lastSolvedApprox as
             sampleList.
         lastSolvedApprox: Parameter(s) corresponding to last computed
             approximate solution(s) as parameterList.
         Q: Numpy 1D vector containing complex coefficients of approximant
             denominator.
         P: Numpy 2D vector whose columns are FE dofs of coefficients of
             approximant numerator.
     """
 
     def setupApprox(self, *args, **kwargs) -> int:
         if self.checkComputedApprox(): return -1
         self.purgeparamsMarginal()
         setupOK = super().setupApprox(*args, **kwargs)
         if self.matchState: self._postApplyC()
         return setupOK
diff --git a/rrompy/reduction_methods/standard/greedy/rational_interpolant_greedy.py b/rrompy/reduction_methods/standard/greedy/rational_interpolant_greedy.py
index 366b989..4d7356c 100644
--- a/rrompy/reduction_methods/standard/greedy/rational_interpolant_greedy.py
+++ b/rrompy/reduction_methods/standard/greedy/rational_interpolant_greedy.py
@@ -1,500 +1,500 @@
 # Copyright (C) 2018-2020 by the RROMPy authors
 #
 # This file is part of RROMPy.
 #
 # RROMPy is free software: you can redistribute it and/or modify
 # it under the terms of the GNU Lesser General Public License as published by
 # the Free Software Foundation, either version 3 of the License, or
 # (at your option) any later version.
 #
 # RROMPy is distributed in the hope that it will be useful,
 # but WITHOUT ANY WARRANTY; without even the implied warranty of
 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
 # GNU Lesser General Public License for more details.
 #
 # You should have received a copy of the GNU Lesser General Public License
 # along with RROMPy. If not, see <http://www.gnu.org/licenses/>.
 #
 
 from copy import deepcopy as copy
 import numpy as np
 from rrompy.hfengines.base.linear_affine_engine import checkIfAffine
 from .generic_greedy_approximant import GenericGreedyApproximant
 from rrompy.utilities.poly_fitting.polynomial import (polybases, polyfitname,
                                                   PolynomialInterpolator as PI,
                                                   polyvander)
 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.verbosity_depth import verbosityManager as vbMng
 from rrompy.utilities.poly_fitting import customFit
 from rrompy.utilities.exception_manager import (RROMPyWarning, RROMPyException,
                                                 RROMPyAssert, RROMPy_FRAGILE)
 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': kind of snapshots orthogonalization; allowed values
                 include 0, 1/2, and 1; defaults to 1, i.e. POD;
             - 'scaleFactorDer': scaling factors for derivative computation;
                 defaults to 'AUTO';
             - 'S': 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;
             - 'samplerTrainSet': 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';
             - 'functionalSolve': strategy for minimization of denominator
-                functional; allowed values include 'NORM', 'DOMINANT', 'NODAL',
-                'BARYCENTRIC_NORM', and 'BARYCENTRIC[_AVERAGE]' (check pdf in
-                main folder for meaning); defaults to 'NORM';
+                functional; allowed values include 'NORM', 'DOMINANT',
+                'BARYCENTRIC[_NORM]', and 'BARYCENTRIC_AVERAGE' (check pdf in
+                main folder for explanation); defaults to 'NORM';
             - 'interpTol': tolerance for interpolation; defaults to None;
             - 'QTol': tolerance for robust rational denominator management;
                 defaults to 0.
             Defaults to empty dict.
         verbosity(optional): Verbosity level. Defaults to 10.
             
     Attributes:
         HFEngine: HF problem solver.
         mu0: Default parameter.
         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': kind of snapshots orthogonalization;
             - '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;
             - 'samplerTrainSet': training sample points generator;
             - 'errorEstimatorKind': kind of error estimator;
             - 'functionalSolve': strategy for minimization of denominator
                 functional;
             - 'interpTol': tolerance for interpolation;
             - 'QTol': 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.
         verbosity: Verbosity level.
         POD: Kind of snapshots orthogonalization.
         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.
         samplerTrainSet: training sample points generator.
         errorEstimatorKind: kind of error estimator.
         functionalSolve: Strategy for minimization of denominator functional.
         interpTol: tolerance for interpolation.
         QTol: 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.
     """
     
     _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"])
         super().__init__(*args, **kwargs)
         self._postInit()
 
     @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
 
     def _polyvanderAuxiliary(self, mus, deg, *args):
         return polyvander(mus, deg, *args)
 
     def getErrorEstimatorDiscrepancy(self, mus:Np1D) -> Np1D:
         """Discrepancy-based residual estimator."""
         checkIfAffine(self.HFEngine, "apply discrepancy-based error estimator",
                       False, self._affine_lvl)
         mus = self.checkParameterList(mus)
         muCTest = self.trainedModel.centerNormalize(mus)
         tMverb, 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.mapParameterList(self.mus)
         muTestEff = self.mapParameterList(mus)
         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
         self.trainedModel.verbosity = tMverb
         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.interpTol)
         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
         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)
         mu_muTestS = mus[idxMaxEst]
         app_muTestSample = self.getApproxReduced(mu_muTestS)
         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)
         else:
             app_muTestSample = dot(self.samplingEngine.projectionMatrix,
                                    app_muTestSample)
         app_muTestSample = sampleList(app_muTestSample)
         if what == "RES":
             errmu = self.HFEngine.residual(mu_muTestS, app_muTestSample,
                                            post_c = False)
             solmu = self.HFEngine.residual(mu_muTestS, None, post_c = False)
             normSol = self.HFEngine.norm(solmu, dual = True)
             normErr = self.HFEngine.norm(errmu, dual = True)
         else:
             for j, mu in enumerate(mu_muTestS):
                 uEx = self.samplingEngine.nextSample(mu)
                 if what == "OUTPUT":
                     uEx = self.HFEngine.applyC(uEx, mu)
                     app_muTS = self.HFEngine.applyC(app_muTestSample[j], mu)
                     if j == 0:
                         app_muTestS = emptySampleList()
                         app_muTestS.reset((len(app_muTS), len(mu_muTestS)),
                                           dtype = app_muTS.dtype)
                     app_muTestS[j] = app_muTS
                 if j == 0:
                     solmu = emptySampleList()
                     solmu.reset((len(uEx), len(mu_muTestS)), dtype = uEx.dtype)
                 solmu[j] = uEx
             if what == "OUTPUT": app_muTestSample = app_muTestS
             errmu = solmu - app_muTestSample
             normSol = self.HFEngine.norm(solmu, is_state = what != "OUTPUT")
             normErr = self.HFEngine.norm(errmu, is_state = what != "OUTPUT")
         errsamples = normErr / normSol
         musT = copy(self.mus)
         musT.append(mu_muTestS)
         musT = self.trainedModel.centerNormalize(musT)
         musC = self.trainedModel.centerNormalize(mus)
         errT = np.zeros((len(musT), len(mu_muTestS)), dtype = np.complex)
         errT[np.arange(len(self.mus), len(musT)),
              np.arange(len(mu_muTestS))] = errsamples * QTest[idxMaxEst]
         vanT = self._polyvanderAuxiliary(musT, self.E + 1, self.polybasis)
         fitOut = customFit(vanT, errT, full = True, rcond = self.interpTol)
         vbMng(self, "MAIN",
               ("Fitting {} samples with degree {} through {}... Conditioning "
                "of LS system: {:.4e}.").format(len(vanT), self.E + 1,
                                        polyfitname(self.polybasis),
                                        fitOut[1][2][0] / fitOut[1][2][-1]), 15)
         vanC = self._polyvanderAuxiliary(musC, self.E + 1, self.polybasis)
         err = np.sum(np.abs(vanC.dot(fitOut[0])), axis = -1) / QTest
         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 = self.checkParameterList(mus)
         tMverb, 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)
         self.trainedModel.verbosity = tMverb
         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]
             coeffTest = np.linalg.solve(vanTrial, valuesTrial)
             errTest = (np.abs(vanTestNextEff - vanTestEff.dot(coeffTest))
                      / np.expand_dims(QTest, 1))
             if only_one: 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])
         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 = self.checkParameterList(mus)
         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 _setSampleBatch(self, maxS:int):
         self.sampleBatchIdx, self.sampleBatchSize, S = -1, 0, 0
         nextBatchSize = 1
         while S + nextBatchSize <= maxS:
             self.sampleBatchIdx += 1
             self.sampleBatchSize = nextBatchSize
             S += self.sampleBatchSize
             nextBatchSize = totalDegreeN(self.npar - 1,
                                          self.sampleBatchIdx + 1)
         return S
 
     def _preliminaryTraining(self):
         """Initialize starting snapshots of solution map."""
         RROMPyAssert(self._mode, message = "Cannot start greedy algorithm.")
         if self.samplingEngine.nsamples > 0: return
         self._S = self._setSampleBatch(self.S)
         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.projectionMatrix
         firstRun = self.trainedModel is None
         if not firstRun: pMat = pMat[:, len(self.trainedModel.data.mus) :]
         self._setupTrainedModel(pMat, not firstRun)
         unstable = 0
         if self.E > 0:
             Q = self._setupDenominator()
         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._setupRational(Q)
             if self.M < self.E:
                 RROMPyWarning(("Instability in numerator computation. "
                                "Aborting."))
                 unstable = 1
         if not unstable:
             self.trainedModel.data.approxParameters = copy(
                                                          self.approxParameters)
         vbMng(self, "DEL", "Done setting up local approximant.", 5)
         self.verbosity += 10
         return unstable
         
     def setupApprox(self, plotEst : str = "NONE") -> int:
         val = super().setupApprox(plotEst)
         if val == 0:
             if (self.errorEstimatorKind[:10] == "LOOK_AHEAD"
             and len(self.mus) < self.samplingEngine.nsamples):
                 while len(self.mus) < self.samplingEngine.nsamples:
                     self.mus.append(self.samplingEngine.mus[len(self.mus)])
                 self.trainedModel = None
                 self._S = self._setSampleBatch(len(self.mus) + 1)
                 self.setupApproxLocal()
             self._setupRational(self.trainedModel.data.Q,
                                 self.trainedModel.data.P)
             self.trainedModel.data.approxParameters = copy(
                                                          self.approxParameters)
         return val
         
     def loadTrainedModel(self, filename:str):
         """Load trained reduced model from file."""
         super().loadTrainedModel(filename)
         self._setSampleBatch(self.S + 1)
diff --git a/rrompy/reduction_methods/standard/rational_interpolant.py b/rrompy/reduction_methods/standard/rational_interpolant.py
index 4c0470c..a72e39d 100644
--- a/rrompy/reduction_methods/standard/rational_interpolant.py
+++ b/rrompy/reduction_methods/standard/rational_interpolant.py
@@ -1,798 +1,707 @@
 # Copyright (C) 2018-2020 by the RROMPy authors
 #
 # This file is part of RROMPy.
 #
 # RROMPy is free software: you can redistribute it and/or modify
 # it under the terms of the GNU Lesser General Public License as published by
 # the Free Software Foundation, either version 3 of the License, or
 # (at your option) any later version.
 #
 # RROMPy is distributed in the hope that it will be useful,
 # but WITHOUT ANY WARRANTY; without even the implied warranty of
 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
 # GNU Lesser General Public License for more details.
 #
 # You should have received a copy of the GNU Lesser General Public License
 # along with RROMPy. If not, see <http://www.gnu.org/licenses/>.
 #
 
 from copy import deepcopy as copy
 import numpy as np
-from scipy.linalg import eigvals
+from scipy.linalg import eig
 from collections.abc import Iterable
 from .generic_standard_approximant import GenericStandardApproximant
 from rrompy.utilities.poly_fitting.polynomial import (
                                             polybases as ppb, polyfitname,
                                             polyvander as pvP, polyTimes,
                                             PolynomialInterpolator as PI,
                                             PolynomialInterpolatorNodal as PIN)
 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.base.types import (Np1D, Np2D, Tuple, List, sampList,
-                                         paramList, interpEng)
+from rrompy.utilities.base.types import (Np1D, Np2D, Tuple, List, paramList,
+                                         interpEng)
 from rrompy.utilities.base import verbosityManager as vbMng
 from rrompy.utilities.numerical import pseudoInverse, dot, baseDistanceMatrix
 from rrompy.utilities.numerical.factorials import multifactorial
 from rrompy.utilities.numerical.hash_derivative import (nextDerivativeIndices,
                                                   hashDerivativeToIdx as hashD,
                                                   hashIdxToDerivative as hashI)
 from rrompy.utilities.numerical.degree import (reduceDegreeN,
-                                          degreeTotalToFull, fullDegreeMaxMask,
-                                          totalDegreeMaxMask)
-from rrompy.solver import Np2DLikeGramian
+                                               degreeTotalToFull,
+                                               fullDegreeMaxMask,
+                                               totalDegreeMaxMask)
 from rrompy.utilities.exception_manager import (RROMPyException, RROMPyAssert,
                                                 RROMPyWarning)
 
 __all__ = ['RationalInterpolant']
 
 def polyTimesTable(P:interpEng, mus:Np1D, reorder:List[int],
                    derIdxs:List[List[List[int]]], scl : Np1D = None) -> Np2D:
     """Table of polynomial products."""
     if not isinstance(P, PI):
         raise RROMPyException(("Polynomial to evaluate must be a polynomial "
                                "interpolator."))
     Pvals = [[0.] * len(derIdx) for derIdx in derIdxs]
     for j, derIdx in enumerate(derIdxs):
         nder = len(derIdx)
         for der in range(nder):
             derI = hashI(der, P.npar)
             Pvals[j][der] = P([mus[j]], derI, scl) / multifactorial(derI)
     return blockDiagDer(Pvals, reorder, derIdxs)
 
 def vanderInvTable(vanInv:Np2D, idxs:List[int], reorder:List[int],
                    derIdxs:List[List[List[int]]]) -> Np2D:
     """Table of Vandermonde pseudo-inverse."""
     S = len(reorder)
     Ts = [None] * len(idxs)
     for k in range(len(idxs)):
         invLocs = [None] * len(derIdxs)
         idxGlob = 0
         for j, derIdx in enumerate(derIdxs):
             nder = len(derIdx)
             idxGlob += nder
             idxLoc = np.arange(S)[(reorder >= idxGlob - nder)
                                 * (reorder < idxGlob)]
             invLocs[j] = vanInv[k, idxLoc]
         Ts[k] = blockDiagDer(invLocs, reorder, derIdxs, [2, 1, 0])
     return Ts
 
 def blockDiagDer(vals:List[Np1D], reorder:List[int],
                  derIdxs:List[List[List[int]]],
                  permute : List[int] = None) -> Np2D:
     """Table of derivative values for point confluence."""
     S = len(reorder)
     T = np.zeros((S, S), dtype = np.complex)
     if permute is None: permute = [0, 1, 2]
     idxGlob = 0
     for j, derIdx in enumerate(derIdxs):
         nder = len(derIdx)
         idxGlob += nder
         idxLoc = np.arange(S)[(reorder >= idxGlob - nder)
                             * (reorder < idxGlob)]
         val = vals[j]
         for derI, derIdxI in enumerate(derIdx):
             for derJ, derIdxJ in enumerate(derIdx):
                 diffIdx = [x - y for (x, y) in zip(derIdxI, derIdxJ)]
                 if all([x >= 0 for x in diffIdx]):
                     diffj = hashD(diffIdx)
                     i1, i2, i3 = np.array([derI, derJ, diffj])[permute]
                     T[idxLoc[i1], idxLoc[i2]] = val[i3]
     return T
 
 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': kind of snapshots orthogonalization; allowed values
                 include 0, 1/2, and 1; defaults to 1, i.e. POD;
             - 'scaleFactorDer': scaling factors for derivative computation;
                 defaults to 'AUTO';
             - 'S': total number of 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;
             - 'radialDirectionalWeightsAdapt': bounds for adaptive rescaling of
                 radial basis weights; defaults to [-1, -1];
             - 'functionalSolve': strategy for minimization of denominator
-                functional; allowed values include 'NORM', 'DOMINANT', 'NODAL',
-                'BARYCENTRIC_NORM', and 'BARYCENTRIC[_AVERAGE]' (check pdf in
+                functional; allowed values include 'NORM', 'DOMINANT',
+                'BARYCENTRIC[_NORM]', and 'BARYCENTRIC_AVERAGE' (check pdf in
                 main folder for explanation); defaults to 'NORM';
             - 'interpTol': tolerance for interpolation; defaults to None;
             - 'QTol': tolerance for robust rational denominator management;
                 defaults to 0.
             Defaults to empty dict.
         verbosity(optional): Verbosity level. Defaults to 10.
             
     Attributes:
         HFEngine: HF problem solver.
         mu0: Default parameter.
         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': kind of snapshots orthogonalization;
             - '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;
             - 'radialDirectionalWeightsAdapt': bounds for adaptive rescaling of
                 radial basis weights;
             - 'functionalSolve': strategy for minimization of denominator
                 functional;
             - 'interpTol': tolerance for interpolation via numpy.polyfit;
             - 'QTol': 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.
         verbosity: Verbosity level.
         POD: Kind of snapshots orthogonalization.
         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.
         radialDirectionalWeightsAdapt: Bounds for adaptive rescaling of radial
             basis weights.
         functionalSolve: Strategy for minimization of denominator functional.
         interpTol: Tolerance for interpolation via numpy.polyfit.
         QTol: 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.
     """
 
-    _allowedFunctionalSolveKinds = ["NORM", "DOMINANT", "NODAL",
-                                    "BARYCENTRIC_NORM", "BARYCENTRIC_AVERAGE"]
+    _allowedFunctionalSolveKinds = ["NORM", "DOMINANT", "BARYCENTRIC_NORM",
+                                    "BARYCENTRIC_AVERAGE"]
     
     def __init__(self, *args, **kwargs):
         self._preInit()
         self._addParametersToList(["polybasis", "M", "N", "polydegreetype",
                                    "radialDirectionalWeights",
                                    "radialDirectionalWeightsAdapt",
                                    "functionalSolve", "interpTol", "QTol"],
                                   ["MONOMIAL", "AUTO", "AUTO", "TOTAL", 1.,
                                    [-1., -1.], "NORM", -1, 0.])
         super().__init__(*args, **kwargs)
         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: 
                 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 functionalSolve(self):
         """Value of functionalSolve."""
         return self._functionalSolve
     @functionalSolve.setter
     def functionalSolve(self, functionalSolve):
         try:
             functionalSolve = functionalSolve.upper().strip().replace(" ","")
-            if functionalSolve == "BARYCENTRIC": functionalSolve += "_AVERAGE"
+            if functionalSolve == "BARYCENTRIC": functionalSolve += "_NORM"
             if functionalSolve not in self._allowedFunctionalSolveKinds:
                 raise RROMPyException(("Prescribed functionalSolve not "
                                        "recognized."))
             self._functionalSolve = functionalSolve
         except:
             RROMPyWarning(("Prescribed functionalSolve not recognized. "
                            "Overriding to 'NORM'."))
             self._functionalSolve = "NORM"
         self._approxParameters["functionalSolve"] = self.functionalSolve
 
     @property
     def interpTol(self):
         """Value of interpTol."""
         return self._interpTol
     @interpTol.setter
     def interpTol(self, interpTol):
         self._interpTol = interpTol
         self._approxParameters["interpTol"] = self.interpTol
 
     @property
     def radialDirectionalWeights(self):
         """Value of radialDirectionalWeights."""
         return self._radialDirectionalWeights
     @radialDirectionalWeights.setter
     def radialDirectionalWeights(self, radialDirectionalWeights):
         if isinstance(radialDirectionalWeights, Iterable):
             radialDirectionalWeights = list(radialDirectionalWeights)
         else:
             radialDirectionalWeights = [radialDirectionalWeights]
         self._radialDirectionalWeights = radialDirectionalWeights
         self._approxParameters["radialDirectionalWeights"] = (
                                                  self.radialDirectionalWeights)
 
     @property
     def radialDirectionalWeightsAdapt(self):
         """Value of radialDirectionalWeightsAdapt."""
         return self._radialDirectionalWeightsAdapt
     @radialDirectionalWeightsAdapt.setter
     def radialDirectionalWeightsAdapt(self, radialDirectionalWeightsAdapt):
         self._radialDirectionalWeightsAdapt = radialDirectionalWeightsAdapt
         self._approxParameters["radialDirectionalWeightsAdapt"] = (
                                             self.radialDirectionalWeightsAdapt)
 
     @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 QTol(self):
         """Value of tolerance for robust rational denominator management."""
         return self._QTol
     @QTol.setter
     def QTol(self, QTol):
         if QTol < 0.:
             RROMPyWarning(("Overriding prescribed negative robustness "
                            "tolerance to 0."))
             QTol = 0.
         self._QTol = QTol
         self._approxParameters["QTol"] = self.QTol
         
     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:
             if self.functionalSolve != "NORM" and self.npar > 1:
                 RROMPyWarning(("Strategy for functional optimization must be "
                                "'NORM' for more than one parameter. "
                                "Overriding to 'NORM'."))
                 self.functionalSolve = "NORM"
             if (self.functionalSolve[:11] == "BARYCENTRIC"
             and self.N + 1 < self.S):
                 RROMPyWarning(("Barycentric strategy cannot be applied with "
                                "Least Squares. Overriding to 'NORM'."))
                 self.functionalSolve = "NORM"
             if self.functionalSolve[:11] == "BARYCENTRIC":
                 invD, TN = None, None
                 self._setupInterpolationIndices()
-            else:
+                if len(self._musUnique) != self.S:
+                    RROMPyWarning(("Barycentric functional optimization "
+                                   "cannot be applied to repeated samples. "
+                                   "Overriding to 'NORM'."))
+                    self.functionalSolve = "NORM"
+            if self.functionalSolve[:11] != "BARYCENTRIC":
                 invD, TN = self._computeInterpolantInverseBlocks()
-            if (self.functionalSolve in ["NODAL", "BARYCENTRIC_NORM",
-                                         "BARYCENTRIC_AVERAGE"]
-            and len(self._musUnique) != self.S):
-                if self.functionalSolve[:11] == "BARYCENTRIC":
-                    warnflag = "Barycentric"
-                    invD, TN = self._computeInterpolantInverseBlocks()
-                else:
-                    warnflag = "Iterative"
-                RROMPyWarning(("{} functional optimization cannot be applied "
-                               "to repeated samples. Overriding to "
-                               "'NORM'.").format(warnflag))
-                self.functionalSolve = "NORM"
-            idxSamplesEff = list(range(self.S))
             if self.POD == 1:
-                ev, eV = self.findeveVGQR(
-                        self.samplingEngine.Rscale[:, idxSamplesEff], invD, TN)
+                sampleE = self.samplingEngine.Rscale
+                Rscaling = None
+            elif self.POD == 1/2:
+                sampleE = self.samplingEngine.samples_normal
+                Rscaling = self.samplingEngine.Rscale
             else:
-                if self.POD == 1/2:
-                    sampleE = self.samplingEngine.samples_normal(idxSamplesEff)
-                    Rscaling = self.samplingEngine.Rscale
-                else:
-                    sampleE = self.samplingEngine.samples(idxSamplesEff)
-                    Rscaling = None
-                ev, eV = self.findeveVGExplicit(sampleE, invD, TN, Rscaling)
-            if self.functionalSolve == "NODAL": break
-            evR = ev / np.max(ev)
-            ts = self.QTol * np.linalg.norm(evR)
-            nevBad = len(ev) - np.sum(np.abs(evR) >= ts)
-            if nevBad <= (self.functionalSolve == "NORM"): break
-            if self.catchInstability:
-                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
+                sampleE = self.samplingEngine.samples
+                Rscaling = None
+            ev, eV = self.findeveVG(sampleE, invD, TN, Rscaling)
+            if self.functionalSolve[:11] == "BARYCENTRIC":
+                evBad = np.abs(ev) < self.QTol * np.linalg.norm(ev)
+            else:
+                evBad = np.abs(ev) < self.QTol * np.abs(ev[-1])
+            nevBad = np.sum(evBad)
+            if not nevBad: break
+            if self.npar == 1:
+                dN = nevBad
+            else: #if self.npar > 1 and self.functionalSolve == "NORM":
+                dN = self.N - reduceDegreeN(self.N, len(eV) - nevBad,
+                                            self.npar, self.polydegreetype)
+            vbMng(self, "MAIN",
+                  ("Smallest {} eigenvalue{} below tolerance. Reducing N by "
+                   "{}.").format(nevBad, "s" * (nevBad > 1), dN), 10)
+            self.N = self.N - dN
+            if self.functionalSolve[:11] == "BARYCENTRIC":
+                eV = eV[evBad == False]
+                break
+        if hasattr(self, "_gram"): del self._gram
         if self.N <= 0:
-            self.N = 0
-            eV = np.ones((1, 1))
-        if self.N > 0 and self.functionalSolve in ["NODAL", "BARYCENTRIC_NORM",
-                                                   "BARYCENTRIC_AVERAGE"]:
-            eV = eV[np.isinf(eV) + np.isnan(eV) == False]
+            self.N, eV = 0, np.ones((1,) * self.npar, dtype = np.complex)
+        if self.N > 0 and self.functionalSolve[:11] == "BARYCENTRIC":
             q = PIN()
             q.polybasis, q.nodes = self.polybasis0, eV
         else:
             q = PI()
             q.npar, q.polybasis = self.npar, self.polybasis0
             if self.polydegreetype == "TOTAL":
                 q.coeffs = degreeTotalToFull(tuple([self.N + 1] * self.npar),
                                              self.npar, eV)
             else:
                 q.coeffs = eV.reshape([self.N + 1] * self.npar)
         vbMng(self, "DEL", "Done computing denominator.", 7)
         return q
 
     def _setupNumerator(self):
         """Compute rational numerator."""
         RROMPyAssert(self._mode, message = "Cannot setup numerator.")
         vbMng(self, "INIT", "Starting computation of numerator.", 7)
         self._setupInterpolationIndices()
         Qevaldiag = polyTimesTable(self.trainedModel.data.Q, self._musUniqueCN,
                                    self._reorder, self._derIdxs,
                                    self.scaleFactorRel)
         if self.POD == 1:
             Qevaldiag = Qevaldiag.dot(self.samplingEngine.Rscale.T)
         elif self.POD == 1/2:
             Qevaldiag = Qevaldiag * self.samplingEngine.Rscale
         if hasattr(self, "_M_isauto"):
             self._setMAuto()
             M = self.M
         else:
             M = reduceDegreeN(self.M, self.S, self.npar, self.polydegreetype)
             if M < self.M:
                 RROMPyWarning(("M too large compared to S. Reducing M by "
                                "{}").format(self.M - M))
                 self.M = M
         while self.M >= 0:
             pParRest = [self.M, self.polybasis, self.verbosity >= 5,
                         self.polydegreetype == "TOTAL",
                         {"derIdxs": self._derIdxs, "reorder": self._reorder,
                          "scl": self.scaleFactorRel}]
             if self.polybasis in ppb:
                 p = PI()
             else:
                 self.computeScaleFactor()
                 rDWEff = np.array([w * f for w, f in zip(
                                                  self.radialDirectionalWeights,
                                                  self.scaleFactor)])
                 pParRest = pParRest[: 2] + [rDWEff] + pParRest[2 :]
                 pParRest[-1]["optimizeScalingBounds"] = (
                                             self.radialDirectionalWeightsAdapt)
                 p = RBI()
             if self.polybasis in ppb + rbpb:
                 pParRest += [{"rcond": self.interpTol}]
             wellCond, msg = p.setupByInterpolation(self._musUniqueCN,
                                                    Qevaldiag, *pParRest)
             vbMng(self, "MAIN", msg, 5)
             if wellCond: break
             vbMng(self, "MAIN", ("Polyfit is poorly conditioned. Reducing M "
                                  "by 1."), 10)
             self.M = self.M - 1
         if self.M < 0:
             raise RROMPyException(("Instability in computation of numerator. "
                                    "Aborting."))
         self.M = M
         vbMng(self, "DEL", "Done computing numerator.", 7)
         return p
 
     def setupApprox(self) -> int:
         """Compute rational interpolant."""
         if self.checkComputedApprox(): return -1
         RROMPyAssert(self._mode, message = "Cannot setup approximant.")
         vbMng(self, "INIT", "Setting up {}.". format(self.name()), 5)
         self.computeSnapshots()
         self._setupTrainedModel(self.samplingEngine.projectionMatrix)
         self._setupRational(self._setupDenominator())
         self.trainedModel.data.approxParameters = copy(self.approxParameters)
         vbMng(self, "DEL", "Done setting up approximant.", 5)
         return 0
 
     def _setupRational(self, Q:interpEng, P : interpEng = None):
         vbMng(self, "INIT", "Starting approximant finalization.", 5)
         self.trainedModel.data.Q = Q
         if P is None: P = self._setupNumerator()
         while self.N > 0 and self.npar == 1:
             if self.HFEngine._ignoreResidues:
                 pls = Q.roots()
                 cfs, projMat = None, None
             else:
                 cfs, pls, _ = rational2heaviside(P, Q)
                 cfs = cfs[: self.N].T
                 if self.POD != 1:
                     projMat = self.samplingEngine.projectionMatrix
                 else:
                     projMat = None
             foci = self.sampler.normalFoci()
             plsA = self.mapParameterList(self.mapParameterList(self.mu0)(0, 0)
                                        + self.scaleFactor * pls, "B")(0)
             idxBad = self.HFEngine.flagBadPolesResiduesAbsolute(plsA, cfs,
                                                                 projMat)
             if not self.HFEngine._ignoreResidues: cfs[:, idxBad] = 0.
             idxBad += self.HFEngine.flagBadPolesResiduesRelative(pls, cfs,
                                                                  projMat, foci)
             idxBad = idxBad > 0
             if not np.any(idxBad): break
             vbMng(self, "MAIN",
                   "Removing {} spurious pole{} out of {}.".format(
                        np.sum(idxBad), "s" * (np.sum(idxBad) > 1), self.N), 10)
             if isinstance(Q, PIN):
                 Q.nodes = Q.nodes[idxBad == False]
             else:
                 Q = PI()
                 Q.npar = self.npar
                 Q.polybasis = self.polybasis0
                 Q.coeffs = np.ones(1, dtype = np.complex)
                 for pl in pls[idxBad == False]:
                     Q.coeffs = polyTimes(Q.coeffs, [- pl, 1.],
                                     Pbasis = Q.polybasis, Rbasis = Q.polybasis)
                 Q.coeffs /= np.linalg.norm(Q.coeffs)
             self.trainedModel.data.Q = Q
             self.N = Q.deg[0]
             P = self._setupNumerator()
         self.trainedModel.data.P = P
         vbMng(self, "DEL", "Terminated approximant finalization.", 5)
 
     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()
         pvPPar = [self.polybasis0, self._derIdxs, self._reorder,
                   self.scaleFactorRel]
-        if hasattr(self, "_M_isauto"): self._setMAuto()
-        E = max(self.M, self.N)
-        fullE = E + 1 == len(self._reorder) == len(self._musUniqueCN)
-        if fullE:
+        full = self.N + 1 == self.S == len(self._musUniqueCN)
+        if full:
             mus = self._musUniqueCN[self._reorder]
             dist = baseDistanceMatrix(mus, magnitude = False)[..., 0]
-            dist[np.arange(E + 1), np.arange(E + 1)] = multifactorial([E])
+            dist[np.arange(self.N + 1),
+                 np.arange(self.N + 1)] = multifactorial([self.N])
             fitinvE = np.prod(dist, axis = 1) ** -1
             vbMng(self, "MAIN",
                   ("Evaluating quasi-Lagrangian basis of degree {} at {} "
-                   "sample points.").format(E, E + 1), 5)
+                   "sample points.").format(self.N, self.N + 1), 5)
             invD = [np.diag(fitinvE)]
+            TN = pvP(self._musUniqueCN, self.N, *pvPPar)
         else:
-            while E >= 0:
+            while self.N >= 0:
                 if self.polydegreetype == "TOTAL":
-                    Eeff = E
-                    idxsB = totalDegreeMaxMask(E, self.npar)
+                    Neff = self.N
+                    idxsB = totalDegreeMaxMask(self.N, self.npar)
                 else: #if self.polydegreetype == "FULL":
-                    Eeff = [E] * self.npar
-                    idxsB = fullDegreeMaxMask(E, self.npar)
-                TE = pvP(self._musUniqueCN, Eeff, *pvPPar)
-                fitOut = pseudoInverse(TE, rcond = self.interpTol, full = True)
+                    Neff = [self.N] * self.npar
+                    idxsB = fullDegreeMaxMask(self.N, self.npar)
+                TN = pvP(self._musUniqueCN, Neff, *pvPPar)
+                fitOut = pseudoInverse(TN, rcond = self.interpTol, full = True)
                 vbMng(self, "MAIN",
                       ("Fitting {} samples with degree {} through {}... "
                        "Conditioning of pseudoinverse system: {:.4e}.").format(
-                                        TE.shape[0], E,
+                                        TN.shape[0], self.N,
                                         polyfitname(self.polybasis0),
                                         fitOut[1][1][0] / fitOut[1][1][-1]), 5)
-                if fitOut[1][0] == TE.shape[1]:
+                if fitOut[1][0] == TN.shape[1]:
                     fitinv = fitOut[0][idxsB, :]
                     break
-                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
+                vbMng(self, "MAIN", 
+                      "Polyfit is poorly conditioned. Reducing N by 1.", 10)
+                self.N = self.N - 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 and not fullE:
-            TN = TE
-        else: #build TN from scratch
-            if self.polydegreetype == "TOTAL":
-                Neff = self.N
-            else: #if self.polydegreetype == "FULL":
-                Neff = [self.N] * self.npar
-            TN = pvP(self._musUniqueCN, Neff, *pvPPar)
         return invD, TN
 
-    def findeveVGExplicit(self, sampleE:sampList, invD:List[Np2D], TN:Np2D,
-                          Rscaling : Np1D = None) -> Tuple[Np1D, Np2D]:
+    def findeveVG(self, sampleE:Np2D, invD:List[Np2D], TN:Np2D,
+                  Rscaling : Np1D = None) -> Tuple[Np1D, Np2D]:
         """
-        Compute explicitly eigenvalues and eigenvectors of rational denominator
-            matrix.
+        Compute eigenvalues and eigenvectors of rational denominator matrix, or
+            of its right chol factor if POD.
         """
-        RROMPyAssert(self._mode, message = "Cannot solve eigenvalue problem.")
-        vbMng(self, "INIT", "Building gramian matrix.", 10)
-        gramian = self.HFEngine.innerProduct(sampleE, sampleE, is_state = True)
-        if Rscaling is not None:
-            gramian = (gramian.T * Rscaling.conj()).T * Rscaling
-        if self.functionalSolve == "NODAL":
-            SEnd = invD[0].shape[1]
-            G0 = np.zeros((SEnd,) * 2, dtype = np.complex)
-        if self.functionalSolve[:11] == "BARYCENTRIC":
-            G = gramian
-            nEnd = len(gramian) - 1
+        RROMPyAssert(self._mode, message = "Cannot solve spectral problem.")
+        if self.POD == 1:
+            if self.functionalSolve[:11] == "BARYCENTRIC":
+                Rstack = sampleE
+            else:
+                vbMng(self, "INIT", "Building generalized half-gramian.",
+                      10)
+                S, eWidth = sampleE.shape[0], len(invD)
+                Rstack = np.zeros((S * eWidth, TN.shape[1]),
+                                  dtype = np.complex)
+                for k in range(eWidth):
+                    Rstack[k * S : (k + 1) * S, :] = dot(sampleE, dot(invD[k],
+                                                                      TN))
+                vbMng(self, "DEL", "Done building half-gramian.", 10)
+            _, s, Vh = np.linalg.svd(Rstack, full_matrices = False)
+            ev, eV = s[::-1], Vh[::-1].T.conj()
+            evExp, probKind = -2., "svd "
         else:
-            nEnd = TN.shape[1]
-            G = np.zeros((nEnd, nEnd), dtype = np.complex)
-            for k in range(len(invD)):
-                iDkN = dot(invD[k], TN)
-                G += dot(dot(gramian, iDkN).T, iDkN.conj()).T
-                if self.functionalSolve == "NODAL":
-                    G0 += dot(dot(gramian, invD[k]).T, invD[k].conj()).T
-        vbMng(self, "DEL", "Done building gramian.", 10)
-        if self.functionalSolve in ["NORM", "BARYCENTRIC_NORM"]:
+            if not hasattr(self, "_gram"):
+                vbMng(self, "INIT", "Building gramian matrix.", 10)
+                self._gram = self.HFEngine.innerProduct(sampleE, sampleE,
+                                                        is_state = True)
+                if Rscaling is not None:
+                    self._gram = (self._gram.T * Rscaling.conj()).T * Rscaling
+                vbMng(self, "DEL", "Done building gramian.", 10)
+            if self.functionalSolve[:11] == "BARYCENTRIC":
+                G = self._gram
+            else:
+                vbMng(self, "INIT", "Building generalized gramian.", 10)
+                G = np.zeros((TN.shape[1],) * 2, dtype = np.complex)
+                for k in range(len(invD)):
+                    iDkN = dot(invD[k], TN)
+                    G += dot(dot(self._gram, iDkN).T, iDkN.conj()).T
+                vbMng(self, "DEL", "Done building gramian.", 10)
             ev, eV = np.linalg.eigh(G)
+            evExp, probKind = -1., "eigen"
+        if (self.functionalSolve in ["NORM", "BARYCENTRIC_NORM"]
+         or np.sum(np.abs(ev) < np.finfo(float).eps * np.abs(ev[-1])
+                                                    * len(ev)) == 1):
             eV = eV[:, 0]
-            if self.functionalSolve == "BARYCENTRIC_NORM":
-                eV = self.findeVBarycentric(eV)
-            problem = "eigenproblem"
+        elif self.functionalSolve == "BARYCENTRIC_AVERAGE":
+            eV = eV.dot(ev ** evExp * np.sum(eV, axis = 0).conj())
         else:
-            if self.functionalSolve == "BARYCENTRIC_AVERAGE":
-                fitOut = pseudoInverse(G, rcond = self.interpTol, full = True)
-                eV = self.findeVBarycentric(np.sum(fitOut[0], axis = 1))
-            else:
-                fitOut = pseudoInverse(G[:-1, :-1], rcond = self.interpTol,
-                                       full = True)
-                eV = np.append(fitOut[0].dot(G[:-1, -1]), -1.)
-            ev = fitOut[1][1][::-1]
-            problem = "linear problem"
+            eV = eV.dot(ev ** evExp * eV[0].conj())
+        ev = ev[1 :]
         vbMng(self, "MAIN",
-              ("Solved {} of size {} with condition number "
-               "{:.4e}.").format(problem, nEnd, ev[-1] / ev[0]), 5)
-        if self.functionalSolve == "NODAL": eV = self.findeVNodal(eV, G0)
-        return ev, eV
-
-    def findeveVGQR(self, RPODE:Np2D, invD:List[Np2D],
-                    TN: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.")
-        vbMng(self, "INIT", "Building half-gramian matrix stack.", 10)
-        if self.functionalSolve == "NODAL":
-            gramian = Np2DLikeGramian(None, dot(RPODE, invD[0]))
+              ("Solved {}problem of size {} with condition number "
+               "{:.4e}.").format(probKind, len(ev), ev[-1] / ev[0]), 5)
         if self.functionalSolve[:11] == "BARYCENTRIC":
-            Rstack = RPODE
-            nEnd = RPODE.shape[1] - 1
-        else:
-            S, nEnd, eWidth = RPODE.shape[0], TN.shape[1], len(invD)
-            Rstack = np.zeros((S * eWidth, nEnd), dtype = np.complex)
-            for k in range(eWidth):
-                Rstack[k * S : (k + 1) * S, :] = dot(RPODE, dot(invD[k], TN))
-        vbMng(self, "DEL", "Done building half-gramian.", 10)
-        if self.functionalSolve in ["NORM", "BARYCENTRIC_NORM",
-                                    "BARYCENTRIC_AVERAGE"]:
-            _, ev, Vh = np.linalg.svd(Rstack, full_matrices = False)
-            if self.functionalSolve in ["NORM", "BARYCENTRIC_NORM"]:
-                eV = Vh[-1, :].conj()
-                ev = ev[::-1]
-            else: #if self.functionalSolve == "BARYCENTRIC_AVERAGE":
-                ev[ev > 0.] **= -2.
-                eV = Vh.T.conj().dot(ev * np.sum(Vh, axis = 1))
-            if self.functionalSolve[:11] == "BARYCENTRIC":
-                eV = self.findeVBarycentric(eV)
-            problem = "svd problem"
-        else:
-            fitOut = pseudoInverse(Rstack[:, :-1], rcond = self.interpTol,
-                                   full = True)
-            ev = fitOut[1][1][::-1]
-            eV = np.append(fitOut[0].dot(Rstack[:, -1]), -1.)
-            problem = "linear problem"
-        vbMng(self, "MAIN",
-              ("Solved {} of size {} with condition number "
-               "{:.4e}.").format(problem, nEnd, ev[-1] / ev[0]), 5)
-        if self.functionalSolve == "NODAL": eV = self.findeVNodal(eV, gramian)
+            N = len(eV)
+            arrow = np.zeros((N + 1,) * 2, dtype = np.complex)
+            arrow[1 :, 0] = 1.
+            arrow[0, 1 :] = eV
+            arrow[np.arange(1, N + 1),
+                  np.arange(1, N + 1)] = self._musUniqueCN[self._reorder[: N]
+                                                                    ].flatten()
+            active = np.eye(N + 1)
+            active[0, 0] = 0.
+            Aev, AeV = eig(arrow, active)
+            AevGood = np.isinf(Aev) + np.isnan(Aev) == False
+            eV, AeV = Aev[AevGood], AeV[:, AevGood]
+            ev = np.sum(np.abs(arrow.dot(AeV) - eV * active.dot(AeV)) ** 2.,
+                        axis = 0) ** .5
         return ev, eV
-
-    def findeVBarycentric(self, baryWeights:Np1D) -> Np1D:
-        RROMPyAssert(self._mode,
-                     message = "Cannot solve optimization problem.")
-        arrow = np.pad(np.diag(self._musUniqueCN[self._reorder].flatten()),
-                       (1, 0), "constant", constant_values = 1.) + 0.j
-        arrow[0, 0] = 0.
-        arrow[0, 1:] = baryWeights
-        active = np.pad(np.eye(len(baryWeights)), (1, 0), "constant")
-        eV = eigvals(arrow, active)
-        return eV[np.isinf(eV) + np.isnan(eV) == False]
         
-    def findeVNodal(self, q0coeffs:Np1D, gram:Np2D, maxiter : int = 25,
-                    tolNewton : float = 1e-10) -> Np1D:
-        RROMPyAssert(self._mode,
-                     message = "Cannot solve optimization problem.")
-        q = PI()
-        q.npar, q.polybasis, q.coeffs = self.npar, self.polybasis0, q0coeffs
-        nodes = q.roots()
-        N = len(nodes)
-        grad = np.zeros(N, dtype = np.complex)
-        hess = np.zeros((N, N), dtype = np.complex)
-        for niter in range(maxiter):
-            plDist = baseDistanceMatrix(self._musUniqueCN[self._reorder],
-                                        nodes, magnitude = False)[:, :, 0]
-            q0, q = np.prod(plDist, axis = 1), []
-            for jS in range(N):
-                mask = np.arange(N) != jS
-                q += [np.prod(plDist[:, mask], axis = 1)]
-                grad[jS] = q[-1].conj().dot(gram.dot(q0))
-                for iS in range(jS + 1):
-                    if iS == jS:
-                        hij = 0.
-                        sij = q[-1].conj().dot(gram.dot(q[-1]))
-                    else:
-                        mask = (np.arange(N) != iS) * (np.arange(N) != jS)
-                        qij = np.prod(plDist[:, mask], axis = 1)
-                        hij = qij.conj().dot(gram.dot(q0))
-                        sij = q[-1].conj().dot(gram.dot(q[iS]))
-                    hess[jS, iS] = hij + sij
-                    if iS < jS: hess[iS, jS] = hij + np.conj(sij)
-            dnodes = np.linalg.solve(hess, grad)
-            nodes += dnodes
-            tol = np.linalg.norm(dnodes) / np.linalg.norm(nodes)
-            if tol < tolNewton: break
-        if niter < maxiter:
-            vbMng(self, "MAIN",
-                  ("Newton's method for problem of size {} converged "
-                   "(err = {:.4e}) in {} iteration{}.").format(
-                                  N + 1, tol, niter + 1, "s" * (niter > 0)), 5)
-        else:
-            RROMPyWarning(("Newton's method for problem of size {} did "
-                           "not converge (err = {:.4e}) after {} "
-                           "iterations.").format(N + 1, tol, niter + 1))
-        return nodes
-
     def getResidues(self, *args, **kwargs) -> Tuple[paramList, Np2D]:
         """
         Obtain approximant residues.
 
         Returns:
             Matrix with residues as columns.
         """
         return self.trainedModel.getResidues(*args, **kwargs)
diff --git a/tests/3_reduction_methods_1D/rational_interpolant_1d.py b/tests/3_reduction_methods_1D/rational_interpolant_1d.py
index 1899737..52779d6 100644
--- a/tests/3_reduction_methods_1D/rational_interpolant_1d.py
+++ b/tests/3_reduction_methods_1D/rational_interpolant_1d.py
@@ -1,67 +1,67 @@
 # Copyright (C) 2018-2020 by the RROMPy authors
 #
 # This file is part of RROMPy.
 #
 # RROMPy is free software: you can redistribute it and/or modify
 # it under the terms of the GNU Lesser General Public License as published by
 # the Free Software Foundation, either version 3 of the License, or
 # (at your option) any later version.
 #
 # RROMPy is distributed in the hope that it will be useful,
 # but WITHOUT ANY WARRANTY; without even the implied warranty of
 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
 # GNU Lesser General Public License for more details.
 #
 # You should have received a copy of the GNU Lesser General Public License
 # along with RROMPy. If not, see <http://www.gnu.org/licenses/>.
 #
 
 import numpy as np
 from matrix_fft import matrixFFT
 from rrompy.reduction_methods import RationalInterpolant as RI
 from rrompy.parameter.parameter_sampling import (QuadratureSampler as QS,
                                                  ManualSampler as MS)
 from rrompy.parameter import checkParameterList
 
 def test_monomials(capsys):
     mu = 1.5
     solver = matrixFFT()
     params = {"POD": False, "S": 10, "QTol": 1e-6, "interpTol": 1e-3,
               "polybasis": "MONOMIAL", "sampler": QS([1.5, 6.5], "UNIFORM")}
     approx = RI(solver, 4., approxParameters = params, verbosity = 10)
     approx.setupApprox()
 
     out, err = capsys.readouterr()
     assert "below tolerance. Reducing N" in out
     assert "poorly conditioned. Reducing M " in out
     assert len(err) == 0
-    assert np.isclose(approx.normErr(mu)[0], 1.4746e-05, atol = 1e-4)
+    assert np.isclose(approx.normErr(mu)[0], 3.454e-2, atol = 1e-2)
 
 def test_well_cond():
     mu = 1.5
     solver = matrixFFT()
     params = {"POD": True, "S": 10, "QTol": 1e-14, "interpTol": 1e-10,
               "polybasis": "CHEBYSHEV", "sampler": QS([1., 7.], "CHEBYSHEV")}
     approx = RI(solver, 4., approxParameters = params, verbosity = 0)
     approx.setupApprox()
     poles = approx.getPoles()
     for lambda_ in np.arange(1, 8):
         assert np.isclose(np.min(np.abs(poles - lambda_)), 0., atol = 1e-4)
     for mu in approx.mus:
         assert np.isclose(approx.normErr(mu)[0], 0., atol = 1e-8)
 
 def test_hermite():
     mu = 1.5
     solver = matrixFFT()
     sampler0 = QS([1., 7.], "CHEBYSHEV")
     points = checkParameterList(np.tile(sampler0.generatePoints(4)(0), 3))
     params = {"POD": True, "S": 12, "polybasis": "CHEBYSHEV",
               "sampler": MS([1., 7.], points = points)}
     approx = RI(solver, 4., approxParameters = params, verbosity = 0)
     approx.setupApprox()
     poles = approx.getPoles()
     for lambda_ in np.arange(1, 8):
         assert np.isclose(np.min(np.abs(poles - lambda_)), 0., atol = 1e-4)
     for mu in approx.mus:
         assert np.isclose(approx.normErr(mu)[0], 0., atol = 1e-8)
 
diff --git a/tests/4_reduction_methods_multiD/rational_interpolant_2d.py b/tests/4_reduction_methods_multiD/rational_interpolant_2d.py
index 5e2cb73..3148314 100644
--- a/tests/4_reduction_methods_multiD/rational_interpolant_2d.py
+++ b/tests/4_reduction_methods_multiD/rational_interpolant_2d.py
@@ -1,77 +1,77 @@
 # Copyright (C) 2018-2020 by the RROMPy authors
 #
 # This file is part of RROMPy.
 #
 # RROMPy is free software: you can redistribute it and/or modify
 # it under the terms of the GNU Lesser General Public License as published by
 # the Free Software Foundation, either version 3 of the License, or
 # (at your option) any later version.
 #
 # RROMPy is distributed in the hope that it will be useful,
 # but WITHOUT ANY WARRANTY; without even the implied warranty of
 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
 # GNU Lesser General Public License for more details.
 #
 # You should have received a copy of the GNU Lesser General Public License
 # along with RROMPy. If not, see <http://www.gnu.org/licenses/>.
 #
 
 import numpy as np
 from matrix_random import matrixRandom
 from rrompy.reduction_methods import RationalInterpolant as RI
 from rrompy.parameter.parameter_sampling import (QuadratureSampler as QS,
                                                  ManualSampler as MS)
 
 def test_monomials(capsys):
     mu = [5.05, 7.1]
     mu0 = [5., 7.]
     solver = matrixRandom()
     uh = solver.solve(mu)[0]
     params = {"POD": False, "S": 16, "QTol": 1e-6, "interpTol": 1e-3,
               "polybasis": "MONOMIAL",
               "sampler": QS([[4.9, 6.85], [5.1, 7.15]], "UNIFORM")}
     approx = RI(solver, mu0, params, verbosity = 100)
     approx.setupApprox()
 
     uhP1 = approx.getApprox(mu)[0]
     errP = approx.getErr(mu)[0]
     errNP = approx.normErr(mu)[0]
     myerrP = uhP1 - uh
     assert np.allclose(np.abs(errP - myerrP), 0., rtol = 1e-3)
     assert np.isclose(solver.norm(errP), errNP, rtol = 1e-3)
     resP = approx.getRes(mu)[0]
     resNP = approx.normRes(mu)
     assert np.isclose(solver.norm(resP), resNP, rtol = 1e-3)
     assert np.allclose(np.abs(resP - (solver.b(mu) - solver.A(mu).dot(uhP1))),
                        0., rtol = 1e-3)
     assert np.isclose(errNP / solver.norm(uh), 5.2667e-05, rtol = 1)
 
     out, err = capsys.readouterr()
-    assert ("poorly conditioned. Reducing E " in out)
+    assert ("poorly conditioned. Reducing N " in out)
     assert len(err) == 0
     
 def test_well_cond():
     mu = [5.05, 7.1]
     mu0 = [5., 7.]
     solver = matrixRandom()
     params = {"POD": True, "M": 3, "N": 3, "S": 16,
               "interpTol": 1e-10, "polybasis": "CHEBYSHEV",
               "sampler": QS([[4.9, 6.85], [5.1, 7.15]], "UNIFORM")}
     approx = RI(solver, mu0, params, verbosity = 0)
     approx.setupApprox()
     assert np.isclose(approx.normErr(mu)[0] / approx.normHF(mu)[0],
                       5.98695e-05, rtol = 1e-1)
 
 def test_hermite():
     mu = [5.05, 7.1]
     mu0 = [5., 7.]
     solver = matrixRandom()
     sampler0 = QS([[4.9, 6.85], [5.1, 7.15]], "UNIFORM")
     params = {"POD": True, "M": 3, "N": 3, "S": 25, "polybasis": "CHEBYSHEV",
               "sampler": MS([[4.9, 6.85], [5.1, 7.15]],
                             points = sampler0.generatePoints(9))}
     approx = RI(solver, mu0, params, verbosity = 0)
     approx.setupApprox()
     assert np.isclose(approx.normErr(mu)[0] / approx.normHF(mu)[0],
                       5.50053e-05, rtol = 5e-1)