\chapter{Implementation}

In this chapter we will present a detailed overview of the implementation and design decisions behind the linear test that decides consistency for the rectangular case in $n$-dimensions, as explained by Theorem \ref{decidability}. In the case when the test indicates consistency for a pair of functions $(f, g)$, representing function and derivative information, respectively, we make use of the linear programming algorithms \ref{LP-min} and \ref{LP-max} to determine the minimal and maximal bounding surfaces. We also implemented a simple GUI that allows the user of this program to visualise the 3D piecewise-linear surfaces whenever the algorithm reports consistency for the two-dimensional case.

\section{Tools}

The entire codebase of the project has been written in the Python programming language. As mentioned in Section \ref{cvxopt-section}, this choice was motivated by the existence of a robust and developer-friendly linear programming framework called CVXOPT, which is implemented in Python. The language features both functional and object oriented programming styles that are expressed using a lightweight and intuitive syntax. 
\\

In addition to being a cross-platform language, Python is a very popular choice among programmers due to the plethora of readily-available packages which can be easily integrated in projects of any scale and complexity. Python also enjoys tremendous support from the academical community, which has boosted the development of a great deal of libraries. Apart from CVXOPT, we used a couple of other dependencies including:

\begin{itemize}
	\item \lit*{numpy} \cite{numpy}: scientific package which offers a variety of convenient features such as powerful $n$-dimensional arrays with flexibility for arbitrary data types, efficient numerical computation and random number capabilities.
	
	\item \lit*{sympy} \cite{sympy}: written entirely in Python, this module enables effective \textit{symbolic computation} for mathematical expressions. Moreover, it has excellent support for polynomials with an arbitrary number of variables e.g. computing the value of a polynomial at a specific point, or the value of a derivative at any given point.
	
	\item \lit*{matplotlib} \cite{matplotlib}: a library that is rather destined for 2D plots, but also features decent support for basic 3D graphs via the \lit*{mplot3d} toolkit. This package allowed very fast development of the front-end component of our framework, which will be detailed later in this chapter.
\end{itemize}

\section{Back-end}

We will now explain the two main components which make up the core implementation of the framework that decides consistency when the derivative approximation is constrained by hyper-rectangles. Despite having dealt with the 2D case in the previous chapter, the implementation below is described as much as possible for the general setting where the domain is $n$-dimensional. For the sake of clarity, we may occasionally illustrate with examples in the two-dimensional case whenever the general instance becomes notationally heavy or harder to reason about. 

\subsection{Input Format and Data Generation} \label{Input Format and Data Generation}

In order to formulate the constraints for the linear programming problems given by Algorithms \ref{LP-min} and \ref{LP-max}, we first need to specify three essential elements to our program: grid information, function information and derivative information respectively. For simplicity, these are all provided in a single file that will be loaded when the final program is run. We impose the following format for each element:

\begin{enumerate}
	\item \textbf{Grid information}: this is given as $n$ separate lines, corresponding to the $n$ axis of the domain. Each line is represented as a sequence of strictly increasing rational numbers which specify the divisions within $[0, 1]$ (including the endpoints) along the $i^{\rm{th}}$ axis, $1\le i\le n$. Thus, each line is assumed to contain at least $2$ values. \\[0.4cm]
	For presentation purposes, let us denote by $p_i\ge 2$ the number of points along the $i^{\rm{th}}$ axis so that any point of the grid can be written as $G_{k_1k_2\ldots k_n}$, where $0\le k_i\le p_i - 1$, for every $1\le i\le n$ (e.g. $G_{00\ldots 0}$ will correspond to the origin $O$ of $\mathbb{R}^n$, while $G_{p_1-1, \ldots, p_n-1}$ will be mapped to $(1, 1, \ldots, 1)\in\mathbb{R}^n$; more specifically, the $i^{\rm{th}}$ coordinate of such a point is given by the $k_i^{\rm{th}}$ rational value located on the $i^{\rm{th}}$ line above).
	
	\item \textbf{Function information}: this is specified as an $n$-dimensional array of intervals, whose size $p_1 \times p_2\times\ldots \times p_n$ is determined by the grid information provided earlier. An entry at an array index $(k_1, k_2, \ldots k_n)$, where $0\le k_i\le p_i-2$ is given as a pair of floating point numbers:
	
		 \[\Big(c_{k_1k_2\ldots k_n}^-, c_{k_1k_2\ldots k_n}^+\Big)\]
		  
		 with $c_{k_1k_2\ldots k_n}^-\le c_{k_1k_2\ldots k_n}^+$, such that it represents the closed and compact interval constraint inside the sub-hyper-rectangle with $2^n$ vertices $G_{k_1 + b_1, \ldots, k_n + b_n}$, where $b_{i}\in\{0, 1\}$, for $1\le i\le n$.\\[0.4cm]
		 It is now crucial to observe that a point $G_{k_1k_2\ldots k_n}$ for which at least some $k_i = p_i - 1$, where $1\le i\le n$, then the function information is already provided within one or more sub-hyper-rectangles that lie on the border of the unit-square $U = [0,1]^n$. We will refer to such points as being \textit{border} grid points (not to be confused with the interpretation given in Section \ref{nd-cons-section} from Chapter 2, which is fundamentally different). For simplicity, we will provide a default pair of $(0, 0)$ at each of these $p_1 \times p_2\times\ldots \times p_n - (p_1 - 1)\times (p_2 - 1)\times\ldots \times (p_n - 1)$ border grid-points so that the $n$-dimensional array is fully initialised. \\[0.4cm]
		Finally, in order to store the $n$-dimensional array within the input file, we use a convenient 2D matrix layout with $p_1\times p_2\times\ldots \times p_{n-1}$ rows and $p_n$ columns.
	
	\item \textbf{Derivative information}: this is also constructed as an $n$-dimensional array of size $p_1 \times p_2\times\ldots \times p_n$ similar to the one presented above for the function information. The only distinction here is that each entry will be given as a tuple of $n$ pairs, where each pair will represent a closed and compact interval constraint for each partial derivative within the corresponding sub-hyper-rectangle. Specifically, an entry belonging to the derivative information array at index $(k_1, k_2, \ldots k_n)$ features the following form:
	\vspace{0.1cm}
		\[\bigg(
			\Big(b^{1-}_{k_1k_2\ldots k_n}, b^{1+}_{k_1k_2\ldots k_n}\Big), 
			\Big(b^{2-}_{k_1k_2\ldots k_n}, b^{2+}_{k_1k_2\ldots k_n}\Big), \ldots,
			\Big(b^{n-}_{k_1k_2\ldots k_n}, b^{n+}_{k_1k_2\ldots k_n}\Big)
		  \bigg),\]		  
		where $0\le k_i\le p_i - 2$ and $b^{i-}_{k_1k_2\ldots k_n}\le b^{i+}_{k_1k_2\ldots k_n}$ are given rational values for each $1\le i\le n$. Similarly, for the indices corresponding to border grid points in the unit-square, we need to initialise the $n$-dimensional array with a default value -- e.g. we can consider each element of the $n$-dimensional tuple to be the interval $(0, 0)$.
\end{enumerate}

\begin{lstlisting}[frame=lines, language=C, captionpos=b, float, caption=Sample input file, label=info]
# Grid information:
0.0 0.6 1.0 
0.0 1.0

# Function information:
# Array shape: (3, 2)
(2.89, 29.64) (0.0, 0.0)
(2.89, 29.73) (0.0, 0.0)
( 0.0,   0.0) (0.0, 0.0)

# Derivative information:
# Array shape: (3, 2)
((-21.16,  5.97), (-13.55, 10.76)) ((0.0, 0.0), (0.0, 0.0))
(( -5.55, 18.32), (-15.08,  6.44)) ((0.0, 0.0), (0.0, 0.0))
((   0.0,   0.0), (   0.0,   0.0)) ((0.0, 0.0), (0.0, 0.0))
\end{lstlisting}

A very minimalistic example is provided in Listing \ref{info} where we considered three points along the first axis, and only the endpoints $0$ and $1$ along the second axis. Observe that for the border points previously defined as having at least one coordinate equal to $1$, the default values of \lit*{(0.0, 0.0)} and \lit*{((0.0, 0.0), (0.0, 0.0))} are used for function and derivative information, respectively.\\

Although the format developed above for specifying grid, function and derivation information can be scaled to an arbitrary dimension $n\ge 2$, it would become extremely cumbersome if we were to input such data manually into a file. Also, apart from very simple cases like the one shown above, we would not have any guarantees as to whether the input is consistent or not. This will immediately pose a significant problem, especially when dealing with a fine-grained structure of the $n$-dimensional grid contained inside the unit-square $U = [0, 1]^n$.
\\

In order to address these issues, we had to come up with a robust method of generating data automatically in the above format. Also, for the purposes of evaluation, we had to ensure that we can actually generate both consistent and inconsistent input so that the correctness of our implementation can then be thoroughly validated. We will present the input generation mechanisms in very great detail in the next chapter dedicated to Evaluation.

\subsection{Linear Programming Algorithms}

Given that we devised a scalable method for representing all the necessary input required for the linear programming algorithms \ref{LP-min} and \ref{LP-max}, we can now dive into their implementation. 
\\

It is first important to notice that the Theorem \ref{decidability} allows us to test for whether a pair of functions $(f, g)$ is consistent or not by constructing an LP problem which will only check for the feasibility of a solution, rather than for optimising a linear objective function. However, either algorithms \ref{LP-min} or \ref{LP-max} can verify the existence of a solution given the provided constraint set while at the same time trying to derive the least and greatest surfaces that are witnesses to consistency. 
\\

Therefore, it makes complete sense to use one or the other optimisation problems given by \ref{LP-min} and \ref{LP-max} for assessing the consistency of a given pair $(f, g)$. Without loss of generality, we will first ask the LP solver to deal with the minimisation problem \ref{LP-min}. In case consistent input has been detected -- and implicitly the minimal witness constructed -- we will then proceed to solving the converse maximisation problem. According to Theorem \ref{min-max}, the latter will be nothing else than solving the minimisation problem again but flipping the sign of the objective function upon multiplication with $-1$, as the constraint sets are otherwise identical for both problems. 
\\

On the other hand, given our reliance on the CVXOPT library which we covered in Section \ref{cvxopt-section}, we had to ensure that the Algorithm \ref{LP-min} (or the corresponding result in higher dimensions) is mapped to the following matrix form \ref{LP-cvxopt} expected by the LP solver (note that we stripped off the equality constraints $\bm{Ax} = \bm{b}$ which are not relevant in this case):

\begin{equation}
\label{cvxopt-simple}
\begin{array}{ll}
\text{Minimise}   & \bm{c}^\top \bm{x} \\
\text{subject to} & \bm{Gx} \le \bm{h}.
\end{array}\\[0.2cm]
\end{equation}

Firstly, it should be clear than all the $p_1 \times p_2\times\ldots \times p_n$ heights in the unit-square $U$ will be the decision variables corresponding to the vector $\bm{x}$ in \ref{cvxopt-simple} above (not to be confused with $\bm{h}$, which is the column vector for the inequality constraints that will be discussed further). As we are trying to minimise the sum of all these heights, it immediately follows that all the $p_1 \times p_2\times\ldots \times p_n$ components of $\bm{c}$ will have to be equal to $1$.
\\

Before moving on to the inequality constraints of the LP problem, we need to make a crucial assumption that will be used throughout the construction of the coefficient matrix $\bm{G}$ and the column vector of right hand sides $\bm{h}$. As pointed out earlier, the vector $\bm{x}$ of decision variables will include the heights $h_{k_1k_2\ldots k_n}$ at each grid point $G_{k_1k_2\ldots k_n}$ in the unit-square $U$, where $0\le k_i\le p_i - 1$, $1\le i\le n$. Since the product between the coefficient matrix $\bm{G}$ and the column vector $\bm{x}$ will involve all the dot products between each row of the matrix and the vector of decision variables, we need to have a precise ordering of the entries $h_{k_1k_2\ldots k_n}$ within $\bm{x}$. This is equivalent to imposing a specific ordering on the indices $k_1k_2\ldots k_n$, where, as before, $0\le k_i\le p_i - 1$ and $1\le i\le n$. The most convenient arrangement is to consider the lexicographical ordering for these indices, that is:
\[(k_1, k_2, \ldots, k_n)\le_{l} (k_1^{\prime}, k_2^{\prime}, \ldots, k_n^{\prime})\] 
holds if and only if:
\[\exists\ m > 0.\ \forall\ i < m.\ \left(k_i = k_i^{\prime} \wedge k_m < k_m^{\prime}\right).\]

With this assumption in mind, we are now ready to tackle the inequality constraints from Algorithm \ref{LP-min}, which need to be written in the matrix form $\bm{Gx} \le \bm{h}$. We should now observe that there are only two types of inequality constraints that can arise from the LP Algorithm \ref{LP-min} (written below for the arbitrary case $n\ge 2$):
\begin{enumerate}
	\item From the function information, the heights corresponding to each grid-point in the unit-square will be constrained by lower and upper bounds of the form:
	\begin{equation}
	\label{type1}
	C_{k_1k_2\ldots k_n}^- \le h_{k_1k_2\ldots k_n} \le C_{k_1k_2\ldots k_n}^+,
	\end{equation}
	
	where $0\le k_i\le p_i - 1$, $1\le i\le n$.
	\item From the derivative information, the differences between adjacent heights at a sub-hyper-rectangle will satisfy inequalities of the form:
	\begin{equation}
	\label{type2}
	B_{k_1k_2\ldots k_n}^{k-} \le h_{k_1^{\prime}k_2^{\prime}\ldots k_n^{\prime}} - h_{k_1k_2\ldots k_n} \le B_{k_1k_2\ldots k_n}^{k+},
	\end{equation}
	where $k_1k_2\ldots k_n$ and $k_1^{\prime}k_2^{\prime}\ldots k_n^{\prime}$ are appropriately chosen indices along the $k^{\rm{th}}$ axis within the sub-hyper-rectangle with vertices $G_{k_1 + b_1, \ldots, k_n + b_n}$, $b_{i}\in\{0, 1\}$, for $0\le k_i\le p_i - 2$, $1\le i\le n$.
\end{enumerate}

However, it is important to observe that the representation \ref{cvxopt-simple} enforced by CVXOPT requires that any linear combination of decision variables to be bounded above by a corresponding upper bound from the column vector $\bm{h}$. Thus, we will need to flip the lower bounds from the above constraints and write them as follows:
\begin{equation}
\label{lower-bounds}
\begin{aligned}
	- h_{k_1k_2\ldots k_n} \le - C_{k_1k_2\ldots k_n}^- \quad && \quad 
	h_{k_1k_2\ldots k_n} - h_{k_1^{\prime}k_2^{\prime}\ldots k_n^{\prime}} \le - B_{k_1k_2\ldots k_n}^{k-}.
\end{aligned}
\end{equation}
\noindent As a result, it suffices to create the coefficient matrix $\bm{G}$ for the right hand sides of (\ref{type1}) and (\ref{type2}), because then $\bm{G}^{\prime} = -\bm{G}$ will be the corresponding coefficient matrix for the lower bounds written in the form \ref{lower-bounds}. We can actually think of the matrix $\bm{G}$ as being split into sub-blocks $\bm{G}_1$ and $\bm{G}_2$, which represent the coefficient matrices for the upper bounds of (\ref{type1}) and the upper bounds of (\ref{type2}), respectively. Corresponding to these $\bm{G}_1$ and $\bm{G}_2$ we will have some column vectors of right hand sides $\bm{h}_1^{+}$ and $\bm{h}_2^{+}$, respectively, whose entries will be nothing else than the values $C_{k_1k_2\ldots k_n}^+$ and $B^{k+}_{k_1k_2\ldots k_n}$, respectively. Similarly, for the inverted lower bounds given by the inequalities (\ref{lower-bounds}), we will have coefficient sub-matrices $-\bm{G}_1$ and $-\bm{G}_2$ with corresponding upper bounds given by $-C_{k_1k_2\ldots k_n}^{-}$ and $-B^{k-}_{k_1k_2\ldots k_n}$, that will be gathered in some column vectors $\bm{h}_1^{-}$ and $\bm{h}_2^{-}$, respectively. 
\\

In summary, by putting together all the constraints from (\ref{type1}) and (\ref{type2}) in the form required by (\ref{cvxopt-simple}), we get that:
\begin{equation}
\label{final-all}
\begin{gathered}
\begin{aligned}
 \bm{G}_1\bm{x}  \le \bm{h}_1^+ \\
-\bm{G}_1\bm{x}  \le \bm{h}_1^- \\
 \bm{G}_2\bm{x}  \le \bm{h}_2^+ \\
-\bm{G}_2\bm{x}  \le \bm{h}_2^-	
\end{aligned}
\end{gathered}
\end{equation}

Note that the construction of the column vectors $\bm{h}_k^-$ and $\bm{h}_k^+$, $k\in\{1, 2\}$ is equivalent to finding the the best bounds $C_{k_1k_2\ldots k_n}^{-}$, $C_{k_1k_2\ldots k_n}^{+}$, $B_{k_1k_2\ldots k_n}^{k-}$, $B_{k_1k_2\ldots k_n}^{k+}$, where $1\le k\le n$, which one can express in terms of the given $c_{k_1k_2\ldots k_n}^{-}$, $c_{k_1k_2\ldots k_n}^{+}$, $b_{k_1k_2\ldots k_n}^{k-}$, $b_{k_1k_2\ldots k_n}^{k+}$ as well as the grid information -- the latter being needed only when computing the gradients as a difference of appropriately chosen consecutive heights. We have described this procedure for the two-dimensional case in Chapter 3, Section \ref{nd-cons-section}. This is not more technically or conceptually involved in the $n$-dimensional setting, as we follow the same procedure described in the particular case: we consider all the sub-hyper-rectangles in the unit-grid and iteratively optimise the bounds for each height and for each particular difference of heights corresponding to all the $n$ partial derivatives.
\\

Thus, it remains to address the construction of the coefficient matrices $\bm{G}_1$ and $\bm{G}_2$, respectively. We will also briefly discuss them in relation to the column vectors of right hand sides $\bm{h}_1^+$ as well as $\bm{h}_2^-$ (note that the case for the inverted lower bounds is similar). Let us take each of them in turn:
\\

\textbf{Construction of $\bm{G}_1$.} According to the upper bound of (\ref{type1}), $\bm{G}_1$ will simply be the \textit{identity matrix} of size $p_1 \times p_2\times\ldots \times p_n$. This is easy to see because relation (\ref{type1}) has $p_1 \times p_2\times\ldots \times p_n$ independent inequalities, each corresponding to a distinct height defined at grid point $G_{k_1k_2\ldots k_n}$. It is also obvious that, since we imposed the lexicographical ordering within the vector $\bm{x}$ of heights, then the entries of type $C_{k_1k_2\ldots k_n}^+$ in the column vector of right hand sides $\bm{h}_1^+$ will follow the same lexicographical ordering for the subscripts $k_1k_2\ldots k_n$.\\

\textbf{Construction of $\bm{G}_2$.} Finally, for the construction of $\bm{G}_2$ we need the following important observation. From the right hand side of (\ref{type2}), we see that it is convenient to consider the inequalities for each $i^{\rm{th}}$ partial derivative in turn. Thus, matrix $\bm{G}_2$ can striped horizontally in $n$ sub-matrices $\bm{G}^i_2$, each corresponding to the $i^{\rm{th}}$ partial derivative, where $1\le i\le n$.\\[0.4cm]
		Let us now focus on such a sub-matrix $\bm{G}_2^d$, for some fixed $d$, $1\le d\le n$. Note that the upper bounds for the inequalities (\ref{type2}) require all the differences between heights at adjacent grid points for \textit{all} the sub-hyper-rectangles within the unit-square $[0, 1]^n$. Since we examine the constraints along the $d^{\rm{th}}$ partial derivative, we can observe that there are exactly:
		\[N_d\equiv (p_d - 1)\cdot\prod_{i \ne d} p_i\]
		such inequalities, given that those differences are only undefined when we have reached a \textit{border} grid point $G_{k_1k_2\ldots k_n}$ for which $k_d = p_d - 1$ (since we cannot step outside the unit-grid to obtain an additional difference of two heights). \\[0.4cm]
		To be more specific, for this  $d^{\rm{th}}$ partial derivative, we can rewrite the upper bounds of (\ref{type2}) as follows:
		\begin{equation}
			\label{all-diff-d}
			\begin{split}
				 h_{k_1, \ldots, k_{d-1}, k_d + 1, k_{d+1}, \ldots, k_n} - h_{k_1, \ldots, k_{d-1}, k_d, k_{d+1}, \ldots, k_n} \le B^{d+}_{k_1k_2\ldots k_n},
			\end{split}
		\end{equation}
		
		where $0\le k_i\le p_i - 1$ for $1\le i\le n$, $i\ne d$ and $0\le k_d\le p_d -2$. \\[0.4cm]
		Finally, we can construct the upper coefficient matrix $\bm{G}_2^+$ since we now have the mathematical expression of all the $N_d$ inequalities given by (\ref{all-diff-d}). Hence, for each of the $N_d$ rows of the sub-matrix in question, we need to determine the column indexes, running from $0$ to $p_1 \times p_2\times\ldots \times p_n$, at which the values $-1$ and $1$ will be placed, so that taking the dot product with $\bm{x}$ will yield all the inequalities (\ref{all-diff-d}). Recall that $\bm{x}$ contains all the decision variables whose indices are ordered lexicographically. Thus, in order to determine the appropriate column indices previously mentioned, we simply need to find the sequencing/row number of each height $h_{k_1k_2\ldots k_n}$ within the column vector $\bm{x}$. But this is now trivial, as given any index $k_1k_2\ldots k_n$, we can compute its $0$-indexed position within the vector $\bm{x}$ as follows:
		\begin{equation}
			\label{pos-index}	
			k_1 \cdot \prod_{i\ge 1} p_i + k_2 \cdot \prod_{i\ge 2} p_i + \ldots + k_{n-1}\cdot p_n + k_n.
		\end{equation}
		Needless to say, the vector of right hand sides $\bm{h}_2^+$ will be similarly row striped in $d$ sub-vectors and the position within such a stripe at which $B^{d+}_{k_1k_2\ldots k_n}$ will be inserted can be computed using the same (\ref{pos-index}) derived earlier.
\\

To better illustrate this rather convoluted discussion, we will show in full the sparse matrix $\bm{G}_2$, together with the the column vector $\bm{h}_2^+$ for a simple example in the two-dimensional case. We will consider, for the sake of presentation, that we have $p_1 = 3$ and $p_2 = 4$ points along the $x$ and $y$ axis, respectively. Thus, there are $3 \times 4 = 12$ decision variable in total and the system $\bm{G}_2\bm{x}\le \bm{h}_2^+$ is given by:

\[
\left[
\arraycolsep=4pt
\begin{array}{rrrr|rrrr|rrrr}
	-1 &  0 &  0 &  0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
	 0 & -1 &  0 &  0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \\
	 0 &  0 & -1 &  0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 \\
	 0 &  0 &  0 & -1 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 \\ \hline
 
     0 &  0 &  0 &  0 & -1 &  0 &  0 &  0 & 1 & 0 & 0 & 0 \\
	 0 &  0 &  0 &  0 &  0 & -1 &  0 &  0 & 0 & 1 & 0 & 0 \\
	 0 &  0 &  0 &  0 &  0 &  0 & -1 &  0 & 0 & 0 & 1 & 0 \\
	 0 &  0 &  0 &  0 &  0 &  0 &  0 & -1 & 0 & 0 & 0 & 1 \\ \hline \hline

	 
	-1 &  1 &  0 &  0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
	 0 & -1 &  1 &  0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 
	 0 &  0 & -1 &  1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ \hline
	 
	 0 &  0 &  0 &  0 & -1 &  1 &  0 & 0 & 0 & 0 & 0 & 0 \\
	 0 &  0 &  0 &  0 & 0  & -1 &  1 & 0 & 0 & 0 & 0 & 0 \\ 
	 0 &  0 &  0 &  1 & 0  &  0 & -1 & 1 & 0 & 0 & 0 & 0 \\ \hline
	 
	 0 &  0 &  0 &  0 & 0 & 0 & 0 & 0 & -1 &  1 &  0 & 0 \\
	 0 &  0 &  0 &  0 & 0 & 0 & 0 & 0 &  0 & -1 &  1 & 0 \\ 
	 0 &  0 &  0 &  1 & 0 & 0 & 0 & 0 &  0 &  0 & -1 & 1 \\
 
\end{array}\right] \left[
\def\arraystretch{1.45}
\begin{array}{r}
	h_{00} \\
	h_{01} \\
	h_{02} \\
	h_{03} \\ \hline
	h_{10} \\
	h_{11} \\
	h_{12} \\
	h_{13} \\ \hline
	h_{20} \\
	h_{21} \\
	h_{22} \\
	h_{23}
\end{array}\right] \le \left[
\def\arraystretch{1.45}
\begin{array}{r}
	B^+_{00} \\
	B^+_{01} \\
	B^+_{02} \\
	B^+_{03} \\ \hline
	B^+_{10} \\
	B^+_{11} \\
	B^+_{12} \\
	B^+_{13} \\ \hline
	B^+_{20} \\
	B^+_{21} \\
	B^+_{22} \\
	B^+_{23}
\end{array}
\right]
\]

\vspace{0.5cm}

\textbf{Putting together the algorithm.} After the thorough analysis developed so far, we now possess all the necessary information to write a procedure that, given $n$-dimensional function and derivative constraints for each of the sub-hyper-rectangles in the unit square $U$, determines whether consistency holds or not. 
\\

The pseudocode Algorithm \ref{final-algo} defines a procedure \textsc{SolveLPAlgorithm} which takes as input interval-valued function information $f$ and rectangular derivative-information $g$ within an $n$-dimensional $grid$. The program builds the constraints in the matrix form \ref{final-all} and then solves the minimisation LP problem to decide consistency. If found, the procedure also solves the converse optimisation problem by inverting the objective function. In such a case, the output will be a pair of $n$-dimensional arrays representing the minimal and maximal bounding surfaces which are witnesses to consistency. Otherwise, an empty pair is returned signalling that the input given by the pair $(f, g)$ is inconsistent.

\begin{algorithm}
\caption{The Final Linear Programming Algorithm}
\label{final-algo}
\begin{algorithmic}[1]
\Procedure {SolveLPAlgorithm}{$f$, $g$, $grid$}
	
	\State $result \gets (\emptyset, \emptyset)$
	\\
	\LineComment {\textsl{Coefficient matrices for function and derivative constraints.}}
	\State $\bm{G}_1 \gets \textsc{BuildFunctionCoefMatrix}(f,\, grid)$
	\State $\bm{G}_2 \gets \textsc{BuildDerivativeCoefMatrix}(g,\, grid)$
	\\
	
	\LineComment {\textsl{Column vectors of optimised lower and upper bound constraints. \phantom{c b c} For function information, this works by taking the $\max$ of all the \phantom{c b c} lower bounds as well as the $\min$ of all the upper bounds for each \phantom{c b c} height that is defined in adjacent sub-hyper-rectangles within the \phantom{c b c} grid. A similar reasoning applies for optimising the overlapping \phantom{c b c} differences of heights using the derivative constraints.}}
	\State $\mathrlap{(C^-_{\text{vec}}, C^+_{\text{vec}})}\phantom{(B^-_{\text{vec}}, B^+_{\text{vec}})} \gets \textsc{BuildFunctionConstraints}(f,\, grid)$
	\State $(B^-_{\text{vec}}, B^+_{\text{vec}}) \gets \textsc{BuildDerivativeConstraints}(g,\, grid)$
	\\
	
	\LineComment {\textsl{Right hand sides corresponding to \ref{final-all}.}}
	\State $(\bm{h}_1^+, \bm{h}_1^-) \gets (C^+_{\text{vec}}, -C^-_{\text{vec}})$
	\State $(\bm{h}_2^+, \bm{h}_2^-) \gets (B^+_{\text{vec}}, -B^-_{\text{vec}})$
	\\
	
	\LineComment {\textsl{Solve the LP minimisation problem by combining constraints \ref{final-all} \phantom{c b c} into single matrix form $\bm{Gx}\le \bm{h}$.}}
	\State $\mathrlap{\bm{c}}\phantom{\bm{G}} \gets \textsc{Vector}(1,\, grid_{\text{noPoints}})$
	\State $\bm{G} \gets \textsc{JoinMatrices}(\bm{G}_1,\, -\bm{G}_1,\, \bm{G}_2,\, -\bm{G}_2)$
	\State $\mathrlap{\bm{h}}\phantom{\bm{G}} \gets \textsc{JoinVectors}(\bm{h}_1^+,\, \bm{h}_1^-,\, \bm{h}_2^+,\, \bm{h}_2^-)$
	\State $\bm{x}_{\min} \gets \textsc{SolveLP}(\bm{c},\, \bm{G},\, \bm{h})$	
	\\
	
	\LineComment {\textsl{If consistency is found, then solve the LP maximisation problem.}}
	\If {$\bm{x}_{\min} \ne \emptyset$}
		\State $\bm{x}_{\max} \gets \textsc{SolveLP}(-\bm{c},\, \bm{G},\, \bm{h})$
		\State $result \gets (\bm{x}_{\min}, \bm{x}_{\max})$
	\EndIf
	\\
	
	\State \textbf{return} $result$
\EndProcedure
\end{algorithmic}
\end{algorithm}

\section{Front-end}

Lastly, we will cover some details regarding the graphical user interface that we implemented alongside the LP algorithms described above. Given that for an $n$-dimensional domain the resulting piecewise linear surfaces are in dimension $(n+1)$, we could only provide 3D visualisation of the geometric objects for a consistent pair defined in the two-dimensional unit square.
\\

In order to achieve this, we leveraged the \lit*{mplot3d} toolkit from the \lit*{matplotlib}\footnote{\url{http://matplotlib.org/mpl_toolkits/mplot3d/}} package, which provides rendering of tri-surface plots via the following convenient API: \footnote{\url{http://matplotlib.org/mpl_toolkits/mplot3d/tutorial.html\#mpl_toolkits.mplot3d.Axes3D.plot_trisurf}}

\begin{centerverbatim}
	
	  Axes3D.plot_tri_surf(X, Y, Z, color, cmap, norm, vmin, vmax, 
	                       shade, triangles),
\end{centerverbatim}

\noindent where \lit*{X}, \lit*{Y} and \lit*{Z} are the data values represented as flat, one-dimensional arrays. All the next parameters offer enough flexibility for customising the look and feel of the resulting plot. We will dedicate a special attention to the last argument, i.e. \lit*{triangles}, which can be used to specify a particular arrangement of the triangular patches that are finally rendered in the 3D scene.

\subsection{Triangulation of Sub-rectangles}

Recall that in Section \ref{nd-cons-section} from Chapter 3, when discussing the construction of the piecewise linear surfaces, we chose to divide each 2D sub-rectangle of the grid along the diagonal joining the lower-right and upper-left corners of any such sub-rectangle. This allows us to linearly interpolate the $z$-values at the vertices of all of these triangles in a consistent manner.
\\

However, without the last parameter, the \lit*{Axes3D.plot_tri_surf} function will employ Delaunay triangulation \cite{delaunay} which performs the partitioning of sub-rectangles in such a way so that it maximises the minimal angle of the resulting 3D piecewise linear triangles. As a consequence, this technique allows for more aesthetically pleasing surfaces because it avoids skinny or close to degenerate triangles. Despite the fact that this triangulation also yields a consistent piecewise linear witness, we nonetheless wanted to follow our own construction that we previously described. To do this, one can override the default Delaunay triangulation by providing a list of vertex-labels upon which the triangulation of any rectangular quadrilateral will be uniquely determined. We illustrate how this labelling works using a simple example. Consider a $3\times 4$ grid within the square $[0, 1]\times [0, 1]$ as follows:

\begin{figure}[H]
\centering
\definecolor{qqqqff}{rgb}{0.,0.,1.}
\definecolor{uuuuuu}{rgb}{0.26666666666666666,0.26666666666666666,0.26666666666666666}
\definecolor{xdxdff}{rgb}{0.49019607843137253,0.49019607843137253,1.}
\begin{tikzpicture}[line cap=round,line join=round,>=triangle 45,x=10.0cm,y=10.0cm,scale=0.78]
\clip(-0.2, -0.1) rectangle (1.2, 1.1);

\draw (0.,1.)-- (0.,0.);
\draw (0.,0.)-- (1.,0.);
\draw (1.,0.)-- (1.,1.);
\draw (1.,1.)-- (0.,1.);

\draw [dash pattern=on 5pt off 5pt] (0.,1.)-- (0.2924724568271015,0.636437351334242);
\draw [dash pattern=on 5pt off 5pt] (0.2924724568271015,1.)-- (1.,0.636437351334242);
\draw [dash pattern=on 5pt off 5pt] (0.,0.636437351334242)-- (0.2924724568271015,0.21324696655940126);
\draw [dash pattern=on 5pt off 5pt] (0.2924724568271015,0.636437351334242)-- (1.,0.21324696655940126);
\draw [dash pattern=on 5pt off 5pt] (0.2924724568271015,0.21324696655940126)-- (1.,0.);
\draw [dash pattern=on 5pt off 5pt] (0.,0.21324696655940126)-- (0.2924724568271015,0.);

\draw (0.2924724568271015,1.)-- (0.2924724568271015,0.);
\draw (0.,0.636437351334242)-- (1.,0.636437351334242);
\draw (1.,0.21324696655940126)-- (0.,0.21324696655940126);
\begin{scriptsize}

% Corner points.
\draw[fill=red]  (0, 0) circle (2.5pt);
\draw[color=red] (0, 0) node [anchor=north east] {$(p_0, q_0)$};
\draw[color=red] (0, 0) node [anchor=south west] {$\boxed{0}$};

\draw[fill=red]  (1, 0) circle (2.5pt);
\draw[color=red] (1, 0) node [anchor=north east] {$(p_2, q_0)$};
\draw[color=red] (1, 0) node [anchor=south west] {$\boxed{8}$};

\draw[fill=red]  (0, 1) circle (2.5pt);
\draw[color=red] (0, 1) node [anchor=north east] {$(p_0, q_3)$};
\draw[color=red] (0, 1) node [anchor=south west] {$\boxed{3}$};

\draw[fill=red]  (1, 1) circle (2.5pt);
\draw[color=red] (1, 1) node [anchor=north east] {$(p_2, q_3)$};
\draw[color=red] (1, 1) node [anchor=south west] {$\boxed{11}$};

% Border points.
\draw [fill=qqzzqq] (0, 0.21324696655940126) circle (2.5pt);
\draw[color=qqzzqq] (0, 0.21324696655940126) node [anchor=north east] {$(p_0, q_1)$};
\draw[color=qqzzqq] (0, 0.21324696655940126) node [anchor=south west] {$\boxed{1}$};

\draw [fill=qqzzqq] (0, 0.636437351334242) circle (2.5pt);
\draw[color=qqzzqq] (0, 0.636437351334242) node [anchor=north east] {$(p_0, q_2)$};
\draw[color=qqzzqq] (0, 0.636437351334242) node [anchor=south west] {$\boxed{2}$};

\draw [fill=qqzzqq] (0.2924724568271015,0.) circle (2.5pt);
\draw[color=qqzzqq] (0.2924724568271015,0.) node [anchor=north east] {$(p_1, q_0)$};
\draw[color=qqzzqq] (0.2924724568271015,0.) node [anchor=south west] {$\boxed{4}$};

\draw [fill=qqzzqq] (0.2924724568271015, 1) circle (2.5pt);
\draw[color=qqzzqq] (0.2924724568271015, 1) node [anchor=north east] {$(p_1, q_3)$};
\draw[color=qqzzqq] (0.2924724568271015, 1) node [anchor=south west] {$\boxed{7}$};

\draw [fill=qqzzqq] (1, 0.636437351334242) circle (2.5pt);
\draw[color=qqzzqq] (1, 0.636437351334242) node [anchor=north east] {$(p_2, q_2)$};
\draw[color=qqzzqq] (1, 0.636437351334242) node [anchor=south west] {$\boxed{10}$};

\draw [fill=qqzzqq] (1, 0.21324696655940126) circle (2.5pt);
\draw[color=qqzzqq] (1, 0.21324696655940126) node [anchor=north east] {$(p_2, q_1)$};
\draw[color=qqzzqq] (1, 0.21324696655940126) node [anchor=south west] {$\boxed{9}$};

% Interior points.
\draw [fill=blue] (0.2924724568271015, 0.636437351334242) circle (2.5pt);
\draw[color=blue] (0.2924724568271015, 0.636437351334242) node [anchor=north east] {$(p_1, q_2)$};
\draw[color=blue] (0.2924724568271015, 0.636437351334242) node [anchor=south west] {$\boxed{6}$};

\draw [fill=blue] (0.2924724568271015, 0.21324696655940126) circle (2.5pt);
\draw[color=blue] (0.2924724568271015, 0.21324696655940126) node [anchor=north east] {$(p_1, q_1)$};
\draw[color=blue] (0.2924724568271015, 0.21324696655940126) node [anchor=south west] {$\boxed{5}$};

\end{scriptsize}
\end{tikzpicture}
\caption{Labelling for a $3\times 4$ grid.}
\label{labels}
\end{figure}

We will now assume that the $z$-values for the heights $h_{ij}$ at each grid-point $(p_i, q_j)$, $i = \overline{0, 2}$, $j=\overline{0, 3}$, have been arranged as a one-dimensional array such that the collection of indices $ij$ are ordered lexicographically. If this is the case, then each two-dimensional grid point $(p_i, q_j)$ will be assigned a label corresponding to the index at which the height $h_{ij}$ is located in the flat array of values. In the above Figure \ref{labels}, the labels are shown in the box next to each grid-point.
\\

Finally, the triangulation expected by \lit*{Axes3D.plot_tri_surf} will simply require a two-dimensional list of all the triangle labels defined using the above procedure. In this example, the \lit*{triangles} parameter would therefore be set to:

\begin{centerverbatim}
	[[0, 4, 1], [4, 5, 1], [1, 5, 2], [5, 6, 2], [2, 6, 3], 
	 [6, 7, 3], [4, 8, 5], [8, 9, 5], [5, 9, 6], [9, 10, 6], 
	 [6, 10, 7], [10, 11, 7]].
\end{centerverbatim}

\vspace{0.3cm}

\begin{figure}[H]
	\centerline{\includegraphics[width=1.05\linewidth]{piecewise-linear}}
	\caption{Triangulation for a consistent piecewise linear witness}
\label{piecewise-linear}
\end{figure}

\vspace{1cm}

Figure \ref{piecewise-linear} shows an example of this triangulation for a piecewise linear surface within a similar $3\times 4$ grid. Note that, for presentation purposes, we chose equally spaced points along the $[0, 1]$ segments. In addition, the 3D scene in the running application can be zoomed, panned or even rotated so that one can analyse the piecewise linear surfaces at a greater level of detail.

\subsection{Minimal and Maximal Surfaces}

We can now easily integrate the visualisation of the least and greatest witnesses to consistency in our framework. Assuming consistent input has been supplied, we retrieve both the minimal and maximal heights via Algorithm \ref{final-algo} and we then make two separate calls to \lit*{Axes3D.plot_tri_surf} by using the same \lit*{x} and \lit*{y} parameters, but different \lit*{z} arguments. 
\\

Figure \ref{min-max-surfaces} shows an example of the least and greatest bounding surfaces that are witnesses to consistency within a $3\times 4$ grid as before:

\begin{figure}[H]
	\centerline{\includegraphics[width=1.05\linewidth]{min-max-surfaces}}
	\caption{Example of minimal and maximal piecewise linear surfaces}
\label{min-max-surfaces}
\end{figure}





 
