\chapter{Evaluation}

In this chapter we will examine the techniques used to ensure that the implementation of the linear programming algorithm \ref{final-algo} that decides consistency is correct. We begin by focusing on the methodology of generating consistent/inconsistent input for both function and derivative information, which will lead to an effective method of performing boundary testing \cite{boundarytesting}. Also, in the particular case when the domain is two-dimensional, we will use the visualisation of the 3D piecewise linear surfaces to further validate our implementation subject to varying some parameters. We then conclude with the challenges and limitations associated with this project.

\section{Input Generation}

Recall that at the end of Section \ref{Input Format and Data Generation} from the previous chapter we established that it is rather impractical to manually input data in the form of closed and compact intervals for function and derivation information, respectively. Even in the case of a two-dimensional domain, the task of writing values into a file is certainly cumbersome, let alone in higher dimensional spaces.
\\

The problem becomes even more serious when we would like to distinguish between consistent and inconsistent input for the purposes of verification. This is because we want to guarantee that the implementation of the LP algorithm is correctly identifying consistent input as being consistent and, conversely, that inconsistent input is classified as being inconsistent. Thus, the absence of any false alarms when testing on a large number of randomly-generated inputs will provide a strong justification towards the correctness of our implementation.
\\

As a result, the task of generating input is of paramount importance to a successful evaluation of this project. We start by looking at how we can automatically produce consistent input, based either on randomly chosen values or randomly-generated polynomials. We will then see how we can create inconsistent input by making a minimal and non-trivial modification to any consistent pair of function and derivative approximation. 

Throughout this analysis, we will make extensive use of the notations developed in Section \ref{Input Format and Data Generation} from the previous chapter. If needed, the reader is strongly encouraged to refer back to that section, as the following discussion will be rather notationally-heavy.

\subsection{Generation of consistent input} 

The problem of generating consistent function and derivative information reduces to reverse engineering the constraints of the linear programming algorithms \ref{LP-min} and \ref{LP-max} from Chapter 3. For this reason, we start by randomly generating heights $h_{k_1k_2\ldots k_n}$ at each of the of the grid points $G_{k_1k_2\ldots k_n}$, where $0\le k_i\le p_i - 1$ and $1\le i\le n$. We implemented two distinct ways for generating such arbitrary heights:
		
	\begin{itemize}
		\item Either by drawing $p_1 \times p_2\times\ldots \times p_n$ values from a uniform distribution over a fixed interval;
		\item Or, more interestingly, by constructing a random polynomial in $n$ variables from which we calculate the values at each of the $p_1 \times p_2\times\ldots \times p_n$ grid points $G_{k_1k_2\ldots k_n}$.
	\end{itemize}
	
\textbf{Consistent Function Information.} After having the heights in place for all the $p_1 \times p_2\times\ldots \times p_n$ grid points, we firstly need to ensure that we construct large enough intervals such that all the heights within a sub-hyper-rectangle satisfy:

	\[c_{k_1k_2\ldots k_n}^- \le h_{k_1 + b_1, \ldots, k_n + b_n} \le  c_{k_1k_2\ldots k_n}^+,\\[0.3cm]\]

\noindent where $0\le k_i\le k_i - 2$ and $b_{i}\in\{0, 1\}$, for $1\le i\le n$. It is now easy to see that we can guarantee consistent function information by choosing the tightest closed and compact intervals 
	
	\[\Big[c_{k_1k_2\ldots k_n}^-, c_{k_1k_2\ldots k_n}^+\Big],\\[0.3cm]\] 

\noindent where $0\le k_1\le p_i - 2$ and $1\le i\le n$. This is trivial because for each of the $(p_1 - 1) \times (p_2 - 1)\times\ldots \times (p_n - 1)$ sub-hyper-rectangles inside the unit-square $U$, the corresponding heights will be contained within the most rigid $(n+1)$-dimensional box. 
\\

However, for added flexibility, we can introduce a parameter $\varepsilon$ to allow for more loose bounding boxes at each of the sub-hyper-rectangles within the grid, so we will define intervals of the form:
	
	\[\Big[c_{k_1k_2\ldots k_n}^- - \varepsilon, c_{k_1k_2\ldots k_n}^+ + \varepsilon\Big],\\[0.3cm]\]
	
	\noindent for all $0\le k_i\le p_i - 2$, $1\le i\le n$.\\
	
	A simple implementation of the ideas discussed thus far is given by the pseudocode procedure \ref{function-algo}. This function returns an $n$-dimensional array of intervals representing consistent function information within each sub-hyper-rectangle of $U$, given random heights $h_{\text{rand}}$, grid information $grid$ and tolerance parameter $\varepsilon$.
	
\vspace{0.3cm}
	
\begin{algorithm}
\caption{Generating Consistent Function Information}
\label{function-algo}
\begin{algorithmic}[1]
\Procedure {GenerateFunctionInfo}{$h_{\text{rand}}$, $grid$, $\varepsilon$}
	\State $f \leftarrow \emptyset$
	\ForAll {$index_{\text{curr}} \in \textsc{GetGridIndices}(grid)$}
		\State $interval \leftarrow (0, 0)$
		\If {$index_{\text{curr}} \notin \textsc{BorderIndexSet}(grid)$}
			\State $(h_{\min}, h_{\max}) \leftarrow (\infty, -\infty)$
			\ForAll {$index_{\text{next}} \in \textsc{NextGridIndices}(index_{\text{curr}})$}
				\State $h \leftarrow h_{\text{rand}}[index_{\text{next}}]$
				\State $(h_{\min}, h_{\max}) \leftarrow (\textsc{Min}(h_{\min},\, h), \textsc{Max}(h_{\max},\, h))$
			\EndFor
		\State $interval \leftarrow (h_{\min} - \varepsilon, h_{\max} + \varepsilon)$
		\EndIf
		\State $f[index_{\text{curr}}] \leftarrow interval$
	\EndFor
	\State \textbf{return} $f$
\EndProcedure
\end{algorithmic}
\end{algorithm}
	
\vspace{0.3cm}	

\textbf{Consistent Derivative Information.} Finally, we also need to deal with the construction of the $n$-dimensional array that will provide the derivative information. Given that we derived the minimal bounding boxes that guarantee consistent function information, we could very easily compute the minimal and maximal slopes along $n$ intersecting hyper-rectangular faces which are pairwise orthogonal. However, this approach may result in some loss of information since we specifically started from random heights at the grid points and most of these heights will be strictly included in the minimal bounding boxes computed previously.
\\

As a result, we will take into account all the randomly generated heights in our strategy for deriving the tightest closed and compact intervals along each of the $n$ axis and for each of the $(p_1 - 1) \times (p_2 - 1)\times\ldots \times (p_n - 1)$ sub-hyper-rectangles in the unit-grid. For convenience, we will explain this method for the case when $n = 2$, as it is similar for higher dimensions.
\\

Recall Figure \ref{derivative-info} in which we had a single sub-rectangle defined by the lower-left corner at $(p_i, q_j)\in\mathbb{R}^2$ and $4$ heights $h_{st}$, $s\in\{i, i+1\}$, $t\in\{t, t+1\}$. As demonstrated by \ref{deriv-x} in Section \ref{nd-cons-section}, the gradients along the $x$ axis are given by the following slopes:
\begin{align}
\label{ratios}
\frac {h_{i+1,j} - h_{ij}}{p_{i+1} - p_i} && \frac {h_{i+1,j+1} - h_{i,j+1}}{p_{i+1} - p_i}.
\end{align}

\vspace{0.3cm}

\begin{algorithm}[H]
\caption{Generating Consistent Derivative Information}
\label{derivative-algo}
\begin{algorithmic}[1]
\Procedure {GenerateDerivativeInfo}{$h_{\text{rand}}$, $grid$, $\varepsilon$}
	\State $d \gets \emptyset$
	\ForAll {$index_{\text{curr}} \in \textsc{GetGridIndices}(grid)$}
		\LineComment {$n$-dimensional tuple of pairs}
		\State $intervals \gets ((0, 0),\ldots (0, 0))$
		\If {$index_{\text{curr}} \notin \textsc{BorderIndexSet}(grid)$}
			\For {$k\gets 1, n$}
				\LineComment{Calculate the edge-length between $index_{\text{curr}}$ and the \phantom{a b c d e f g....\,\,} next index along the $k^{\rm{th}}$ axis for the current sub-hyper- \phantom{.........................\,}rectangle in the grid}
				\State $edge_{\text{len}} \gets \textsc{GetEdgeLength}(grid,\, index_{\text{curr}},\, k)$
				\State $(s_{\min}, s_{\max}) \leftarrow (\infty, -\infty)$
				\State $I_{\text{edges}} \gets \textsc{IndicesForParallelEdges}(grid,\, index_{\text{curr}},\, k)$
				\ForAll {$(index_{\text{left}}, index_{\text{right}}) \in I_{\text{edges}}$}
					\State $(h_{\text{left}}, h_{\text{right}}) \gets (h_{\text{rand}}[index_{\text{left}}], h_{\text{rand}}[index_{\text{right}}])$
					\State $s \gets (h_{\text{right}} - h_{\text{left}}) / edge_{\text{len}}$
					\State $(s_{\min}, s_{\max}) \leftarrow (\textsc{Min}(s_{\min},\, s), \textsc{Max}(s_{\max},\, s))$
				\EndFor 
				\State $intervals[i - 1] \gets (s_{\min} - \varepsilon, s_{\max} + \varepsilon)$
			\EndFor
		\EndIf
		\State $d[index_{\text{curr}}] \gets intervals$
	\EndFor
	\State \textbf{return} $d$
\EndProcedure
\end{algorithmic}
\end{algorithm}

Notice that both ratios have in common the distance given by $p_{i+1} - p_i$, which represents the edge-length of the sub-rectangle in the $x$ direction. Also, it is obvious that in this simple two-dimensional case we would just construct a consistent interval from the values of these 2 ratios.
\\

With these details in mind, we can explain the most challenging parts in the pseudocode procedure \ref{derivative-algo}. We use the same input parameters as before and we output an $n$-dimensional array of tuples, with each such tuple consisting of $n$ pairs. Now, on line $9$ we effectively determine the common edge-length for the current sub-hyper-rectangle in the direction of the $k^{\rm{th}}$ derivative. Secondly, on line $11$ we collect all $2^{n-1}$ pairs of grid indices for which the corresponding heights will need to be subtracted, in the style shown by (\ref{ratios}) for the 2D case. In the particular instance that we just analyzed, the \textsc{IndicesForParallelEdges} procedure would return the set with $2$ elements:

\[\left\{\Big((i, j), (i+1, j)\Big), \Big((i, j+1), (i+1, j+1)\Big)\right\}\\[0.3cm]\]
\noindent since $index_{\text{curr}} = (i, j)$ and $k = 1$ (i.e. the partial derivative along the $x$ axis). Lastly, one can specify the tolerance parameter $\varepsilon$ to accommodate wider intervals for each partial derivative, in the same manner as we did when extending the minimal bounding boxes for function information.     

\subsection{Generation of inconsistent input} 

Having seen a rather involved procedure for constructing consistent input, one may think that generating inconsistent function or derivative information would be much easier. This is indeed the case, and we have two ways for tackling this problem:

\begin{itemize}
	\item We can either break the function information by making two adjacent bounding boxes to have empty intersection -- thus, no continuous witness will exist in such a case;
	\item Or, we can make a minimal change to break the derivative information as follows: we consider to have generated the most constrained intervals for both function and derivative information, as shown in Algorithms \ref{function-algo} and \ref{derivative-algo} when $\varepsilon = 0$ is supplied to both procedures. This means that we have the most rigid $(n+1)$-dimensional hyper-rectangles and also the tightest intervals across all partial derivatives. Now we can pick a random sub-hyper-rectangle in the unit-square $U$ at grid index $G_{r_1r_2\ldots r_n}$, as well as a random axis $r_a$, for which we have the following interval:
	\begin{equation}
		\label{b-interval}
		\Big[b^{r_a^-}_{r_1r_2\ldots r_n}, b^{r_a^+}_{r_1r_2\ldots r_n}\Big]
	\end{equation}
	where $0\le r_i\le p_i - 2$, $0\le i\le n$ and $1\le r\le n$. Now, since the corresponding minimal bounding box is given in terms of the closed and compact interval:
	\[\Big[c^{-}_{r_1r_2\ldots r_n}, c^{+}_{r_1r_2\ldots r_n}\Big],\]
	it means that the greatest consistent interval for the $r_a^{\rm{th}}$ partial derivative is $S^{r_a}\equiv \Big[-s^{r_a}_{\max}, s^{r_a}_{\max}\Big]$, where
	\[s^{r_a}_{\max} \equiv \frac {c^{+}_{r_1r_2\ldots r_n} - c^{-}_{r_1r_2\ldots r_n}}{G^{r_a}_{r_1\ldots (r_a + 1)\ldots r_n} - G^{r_a}_{r_1\ldots r_a\ldots r_n}}\]
	represents the largest achievable slope within the corresponding sub-hyper-rectangle along this randomly chosen derivative. Notice that we use $G_{k_1k_2\ldots k_n}^{k}$ to denote the $k^{\rm{th}}$ coordinate of the grid point $G_{k_1k_2\ldots k_n}$, where $0\le k_i\le p_i - 1$, $1\le i\le n$, $1\le k\le n$. As a matter of fact, the denominator of $s_{\max}^{r_a}$ would be the output of the virtual procedure \textsc{GetEdgeLength} within \textsc{GenerateDerivativeInfo}, whose pseudocode was given as part of Algorithm \ref{derivative-algo}.\\[0.4cm]
	At last, coming back to the single modification that is needed in order to obtain inconsistent derivative information, we need to see that the interval from (\ref{b-interval}) will always be included within $S^{r_a}$, i.e.
	\[\Big[b^{r_a^-}_{r_1r_2\ldots r_n}, b^{r_a^+}_{r_1r_2\ldots r_n}\Big] \subseteq \Big[-s^{r_a}_{\max}, s^{r_a}_{\max}\Big].\]
	Hence, by choosing any interval that is strictly outside of $S^{r_a}$ it will be sufficient to guarantee inconsistent input. This is because the slope along the $r_a^{\rm{th}}$ partial derivative within the minimal $(n+1)$-dimensional bounding box under consideration will not have the necessary freedom to oscillate beyond the limits of $S^{r_a}$.
	
\end{itemize}

For the purposes of verifying our implementation described in Algorithm \ref{final-algo}, we decided to choose for second alternative when generating inconsistent data. This is easily motivated because it enables us to perform effective boundary testing \cite{boundarytesting}, as we are able to create inconsistent input by making the minimal non-trivial change to a consistent pair of function and derivative information, respectively, upon altering the latter.

\section{Results}

Having developed robust methods for automatically generating both consistent and inconsistent input, it was now possible to assess the correctness of our implementation of the Algorithm \ref{final-algo}. 
\\

We checked the implementation against \textbf{hundreds} of randomly generated input files for dimensions $n\in\{2, 3, 4, 5, 6\}$ and considering at most $12$ points along each axis. Furthermore, we allowed the parameter $\varepsilon\ge 0$ to disturb both the function and derivative information whenever generating consistent input -- which was either coming from random polynomials or from arbitrary values drawn from the interval $[10, 30]$. We chose $\varepsilon\in \{0.0, 0.1, 0.5, 1.0, 10.0, 50.0\}$ in order to cover a wider spectrum of test cases, although it is only when $\varepsilon \to 0$ that the generated test cases are more relevant.
\\

The results are indeed very encouraging. After a couple of weeks spent to address various issues in the implementation outlined in Algorithm \ref{final-algo}, we managed to successfully classify each type of input when running a substantial suite of tests for all dimensions considered under analysis ($\approx 1000$ tests in total). As such, in the absence of any counterexamples, we have the confidence that we delivered a solid implementation of the linear programming algorithm that decides consistency in the rectangular case. 
\\

We do not, however, exclude the possibility of any bugs being missed throughout this process: despite having generated as many various tests as possible, there might still be rather unusual situations which may very well uncover bugs in our framework, e.g. precision errors due to floating-point arithmetic, bugs in the CVXOPT solver etc.

\subsection{Examples of 3D Piecewise Linear Surfaces}

As a final method of evaluation, let us now visually inspect and comment on some of the piecewise linear maps that we have briefly seen as part of the Front-End Section in the previous chapter. We will start by looking at the simplest example in the two-dimensional domain, that is, when the grid is given by the corners of the unit-square $U = [0, 1]^2$. We will then gradually increase the granularity of the grid, together with varying the $\varepsilon$ parameter. In all the plots that follow, the randomly-generated surfaces from which consistent data has been constructed are coloured with \textcolor{qqzzqq}{green}. \textcolor{blue}{Blue} and \textcolor{red}{red} are reserved for the minimal and maximal consistent surfaces, respectively.

Figures \ref{2-2-e0} and \ref{2-2-e1} show two examples of consistent input when the grid coincides with the unit-square. The randomly generated surface in the middle is in fact the same in both cases, but the constraints for function and derivative approximation are altered by means of $\epsilon\in\{0, 1\}$. 
\\

Note that when $\varepsilon=0$, the least and greatest piecewise linear witnesses will just touch the original surface at its global minimum and maximum, respectively. This is indeed expected as these extreme values determine the most rigid rectangular parallelepiped in which any witness will lie. In addition, we should also observe that the minimal and maximal surfaces are slanted towards the interior of this tightest bounding box, because the intervals for derivative constraints did not contain $0$ for sure -- as otherwise these witnesses would have been parallel to the $xy$-plane.
\\

\begin{figure}[H]
    \centerline{\includegraphics[width=0.95\linewidth]{2-2-e0}}
    \caption{$2\times 2$ grid and $\varepsilon = 0$.}
    \label{2-2-e0}
\end{figure}

Now, when $\varepsilon = 1$ in Figure \ref{2-2-e1}, we get to see that the minimal and maximal surfaces are separated from the initial one by exactly this value of $\varepsilon$ at the extreme points of the original witness. Moreover, the global bounding surfaces are flatter this time compared to the ones in Figure \ref{2-2-e0} where $\varepsilon = 0$. This is because an $\varepsilon > 0$ will result in wider intervals for the derivative constraints, thus allowing the gradients to have more freedom.

\begin{figure}[H]
    \centerline{\includegraphics[width=0.95\linewidth]{2-2-e1}}
    \caption{$2\times 2$ grid and $\varepsilon = 1$}
     \label{2-2-e1}
\end{figure}

Next, Figures \ref{3-3-e0}, \ref{3-3-e1} and \ref{3-3-e5} show the case of a $3\times 3$ grid with equally spaced points along each axis. In all of these scenarios we end up with much richer surfaces given the more granular partitioning of the unit-square.
\\

\begin{figure}[H]
    \centerline{\includegraphics[width=0.95\linewidth]{3-3-e0}}
    \caption{$3\times 3$ grid and $\varepsilon = 0$}
    \label{3-3-e0}
\end{figure}

\begin{figure}[H]
    \centerline{\includegraphics[width=0.95\linewidth]{3-3-e1}}
    \caption{$3\times 3$ grid and $\varepsilon = 1$}
    \label{3-3-e1}
\end{figure}

As in the case of a $2\times 2$ grid, the witnesses are more apart from each other as the $\varepsilon$ increases. Similarly, the maximal and minimal surfaces are flatter, due to more permissive constraints on the derivative information.
\\

\begin{figure}[H]
    \centerline{\includegraphics[width=0.95\linewidth]{3-3-e5}}
    \caption{$3\times 3$ grid and $\varepsilon = 5$}
    \label{3-3-e5}
\end{figure}

Apart from randomly generating values within an interval, we also have the ability to generate consistent input from random polynomials. In the two-dimensional case,  these polynomials have the following general form:

\[P(x_1, x_2) \equiv \sum_{i=1}^{10} \alpha_i x_1^{e_i^1} x_2^{e_i^2},\\[0.3cm]\]

\noindent where the coefficients $\alpha_i\in [-40, 40]$ and exponents $e_i^1, e_i^2\in [0, 10]$ are randomly chosen, for all $1\le i\le 10$. We opted for a rather small number of terms when generating polynomials, because otherwise evaluating $P$ for each point in the unit grid would require a great deal of unnecessary arithmetic operations.
\\

As we can see in Figures \ref{6-5-e0-poly}, \ref{11-11-e0-poly} and \ref{11-11-e30-poly}, the resulting surfaces feature a much smoother geometry than before. Note that $\varepsilon = 0$ for the first two plots. Since there is no additional offset introduced into the constraints, the original \textcolor{qqzzqq}{green} surface is now ``sandwiched'' between the least and greatest consistent witnesses, for all the sub-rectangles in the unit square $[0, 1]^2$. Finally, in Figure \ref{11-11-e30-poly}, a larger offset $\varepsilon = 30$ produces separate minimal and maximal surfaces which is what we anticipated.
\\

\begin{figure}[H]
    \centerline{\includegraphics[width=0.95\linewidth]{6-5-e0-poly}}
    \caption{$6\times 5$ grid and $\varepsilon = 0$}
    \label{6-5-e0-poly}
\end{figure}
 

 \begin{figure}[H]
    \centerline{\includegraphics[width=0.95\linewidth]{11-11-e0-poly}}
    \caption{$11\times 11$ grid and $\varepsilon = 0$}
	\label{11-11-e0-poly}
\end{figure}

In addition, we also implemented random generation of grid points, as they do not necessarily need to be equally spaced along the axis of the grid. This is depicted in Figures \ref{11-11-e0-poly} and \ref{11-11-e30-poly}.
\\

\begin{figure}[H]
    \centerline{\includegraphics[width=0.95\linewidth]{11-11-e30-poly}}
    \caption{$11\times 11$ grid and $\varepsilon = 30$}
    \label{11-11-e30-poly}
\end{figure}

\begin{figure}[H]
    \centerline{\includegraphics[width=0.95\linewidth]{11-11-e10}}
    \caption{$11\times 11$ grid and $\varepsilon = 10$}
    \label{11-11-e10}
\end{figure}

For the remaining plots, only arbitrary values are used to create the input. Figures \ref{11-11-e10}, \ref{11-11-e50} show an $11\times 11$ grid in which $\varepsilon\in \{10, 50\}$. As expected, the surfaces are much closer together for the smaller value of the $\varepsilon$.  
\\

\begin{figure}[H]
    \centerline{\includegraphics[width=0.95\linewidth]{11-11-e50}}
    \caption{$11\times 11$ grid and $\varepsilon = 50$}
    \label{11-11-e50}
\end{figure}

\begin{figure}[H]
    \centerline{\includegraphics[width=0.95\linewidth]{20-20-e80}}
    \caption{$20\times 20$ grid and $\varepsilon = 80$}
    \label{20-20-e80}
\end{figure}

Upon increasing the grain size of the grid, the piecewise linear objects present a much richer geometrical structure, as we can examine in Figures \ref{20-20-e80} and \ref{40-40-e80}, respectively.
\\

\begin{figure}[H]
    \centerline{\includegraphics[width=0.95\linewidth]{40-40-e80}}
    \caption{$40\times 40$ grid and $\varepsilon = 80$}
    \label{40-40-e80}
\end{figure}

\section{Challenges \& Limitations}

In this section, we will outline some of the limitations that arise in the implementation of our framework that decides consistency in the rectangular case. The major bottleneck concerns the size of the linear programming problem which involves no less than $p_1\times p_2\times\ldots \times p_n$ number of decision variables, where $p_i$ is the number of points along the $i^{\rm{th}}$ axis, $1\le i\le n$.
\\

To show why this presents an enormous problem, we carried out a very simple performance analysis to evaluate the time spent in each of the core components of our implementation. We run $10$ separate tests for each dimension ranging from $n=2$ to $n=9$, where we chose the hyper-rectangle $[0, 1]^n$ as the grid information -- hence each LP problem involves $2^n$ decision variables. The averaged results across all of these experiments are shown in Figure \ref{performance} below: 
\\

\begin{figure}[H]
    \centerline{\includegraphics[width=0.95\linewidth]{performance}}
    \caption{Performance analysis for dimensions $n=\overline{2, 9}$, where the grid coincides with the $n$-dimensional unit square. Measurements are averaged across $10$ separate runs using randomly generated input.}
    \label{performance}
\end{figure}
 
\vspace{1cm}

Even in this simplest scenario where we have the least number of points along each axis, i.e. $p_i = 2$, for $1\le i \le n$, we can notice that the total time spent within each of the main components (input generation, construction of matrix form $\bm{Gx}\le \bm{b}$ and both LP algorithms) nearly doubles with every increment of the domain dimension $n$. This is clearly expected, as we need to build an exponentially bigger input which will also reflect in the size of the resulting linear programming problem.
\\

However, as it obvious from the previous plot, the time spent in the components prior to invoking the CVXOPT solver dominate more than half of the total running time. Thus, there is a scope for improvement in the areas of input generation and construction of final matrix form, which would involve careful parallelisation of the current codebase. The main challenge here would be to decompose the overlapping structure of the function and derivative approximation into disjoint components, as most of the heights use information not only from their designated sub-hyper-rectangle, but also from the adjacent ones. 
