\chapter{Consistency of function and derivative constraints}

In this section we focus on the main idea of the proposed project, which involves the concept of \textit{consistency} between function and derivative information for a given real-valued map defined on an $n$-dimensional domain. More specifically, the basic idea is to determine whether a map $h : \mathbb{R}^n \to \mathbb{R}$ can be constructed given restrictions on the function value and derivative information, respectively. The nature of the derivative constraints can pose a major challenge towards the linear-programming algorithms that will be developed further, depending on whether we consider hyper-rectangles or convex-polyhedra in an $n$-dimensional setup. 
\\

We begin by stating the main notations and terminology on interval-valued maps, Lipschitz functions and their derivatives (see \cite{ref2}). We then proceed to defining the consistency relationship in a domain-theoretical setting and towards the end of the chapter we present the algorithms that can decide consistency for maps defined on $\mathbb{R}^n$, $n\ge 1$. In what follows, we consider $n\ge 1$, unless otherwise stated.

\section{Notations and terminology}

We denote by $\mathbb{R}$ the set of real numbers and by $\mathbf{I}\mathbb{R} = \{[a, b]\ |\ a\le b\in\mathbb{R}\}\cup\{\mathbb{R}\}$ the interval domain, i.e. the set of compact, nonempty intervals, equipped with a least element $\bot = \mathbb{R}$, ordered by reverse inclusion. It has a canonical basis consisting of all compact intervals with rational end points augmented with $\bot$. We write a non-bottom element $v\in\mathbf{I}\mathbb{R}$ as $v = [v^-, v^+]$ and we identify any real number $x\in\mathbb{R}$ with the singleton $\{x\}\subset\mathbb{R}$.
\\

Also, we denote by $\mathbf{I}\mathbb{R}^n$ the product domain consisting of all non-empty compact hyper-rectangles with faces parallel to the standard coordinate planes ordered with reverse inclusion and augmented with the whole space $\mathbb{R}^n$ as the bottom element. It has a canonical basis consisting of all its rational (compact) hyper-rectangles and the bottom element. We denote the continuous Scott domain of the nonempty, compact and convex subsets of $\mathbb{R}^n$, taken together with $\mathbb{R}^n$ as the bottom element and ordered by reverse inclusion, by $\mathbf{C}\mathbb{R}^n$. We will use a canonical basis of $\mathbf{C}\mathbb{R}^n$, consisting of rational convex compact polyhedra together with the set $\mathbb{R}^n$ as the bottom element.
\\

For an open subset $U\subset\mathbb{R}^n$, let $C^0(U)$ be the function space of all continuous functions of type $U\to\mathbb{R}$. We will also use domains of function spaces of the form $(U\to D)$ where $D$ is a countably based continuous dcpo, which is either $\mathbf{I}\mathbb{R}$, $\mathbf{I}\mathbb{R}^n$ or $\mathbf{C}\mathbb{R}^n$ in our case. For the sake of convenience, denote $D^0(U) = U\to\mathbf{I}\mathbb{R}$. A function $f\in D^0(U)$ is given by a pair of respectively lower and upper semi-continuous functions $f^-, f^+:U\to\mathbb{R}$ with $f(x) = [f^-(x), f^+(x)]$ when $f(x)\ne\bot$ for all $x\in U$. Recall that given an open subset $a\subset U$ and an element $b\in D$, the single step function $b\chi_a : X\to D$ is deﬁned as $(b\chi_a)(x) = b$ if $x\in a$ and $\bot$ otherwise. Single-step functions are continuous with respect to the Scott topology. Any finite set of single-step functions that are bounded in the function space $U\to D$ has a least upper bound, called a step function; the set of step functions provides a basis for the continuous Scott domain $U\to D$. This basis in turn gives a countable and canonical basis of rational step functions for $U\to D$, where $D = \mathbf{I}\mathbb{R}$, $\mathbf{I}\mathbb{R}^n$ or $\mathbf{C}\mathbb{R}^n$, generated by single-step functions of the form $b\chi_a$ where $a$ is a rational open hyper-rectangle with faces parallel to the coordinate hyper-planes of $\mathbb{R}^n$ and $b$ is a rational interval for $D = \mathbf{I}\mathbb{R}$, a rational hyper-rectangle for $D = \mathbf{I}\mathbb{R}^n$ and a rational compact convex polyhedron in $\mathbb{R}^n$ for $D = \mathbf{C}\mathbb{R}^n$. Finally, the set of elements above an element $c$ in a domain will be denoted by $\uparrow c$.
\\

Furthermore, we will make use of two operations from interval arithmetic which extend the conventional addition and multiplication of numbers by point-wise application to sets of points. Recall that $\norm{\bm{x}} = \sqrt{\sum_{i=1}^nx_i^2}$ is the standard Euclidean norm of $\bm{x} = [x_1, x_2, \ldots x_n]^\top\in\mathbb{R}^n$. Then the Euclidean norm is extended point-wise to $b\in\mathbf{C}\mathbb{R}^n$ by $\norm{\bm{b}} = \max\{\norm{\bm{x}}\ :\ x\in \bm{b}\}$. We will also consider the extension $- \cdot - : \mathbf{C}\mathbb{R}^n \times \mathbb{R}^n \to\ \mathbf{I}\mathbb{R}$ of the scalar product which is defined point-wise $\bm{b}\cdot \bm{x} = \{\bm{y} \cdot \bm{x}\ :\ \bm{y} \in \bm{b}\}$.
\\

Let us now recall the several well-known definitions from calculus regarding directional derivatives and Lipschitz functions.

\begin{definition} (Partial Derivative)
	If $f : \mathbb{R}^n\to\mathbb{R}$, then its partial derivative with respect to dimension $i$ is defined as:
	
	\[\frac {\partial f(\bm{x})}{\partial x_i} = \lim_{h\to 0} \frac {f(x_1, \ldots, x_i + h \ldots x_n) - f(x_1, \ldots, x_i, \ldots, x_n)}h\\[0.2cm]\]
\end{definition}

\begin{definition} (Gradient)
	If $f : \mathbb{R}^n\to\mathbb{R}$, then the vector of partial derivatives is called the \textit{gradient} and is defined by:
	
	\[\nabla f(\bm{x}) = 
		\begin{bmatrix}
 			\frac {\partial f(\bm{x})}{\partial x_1} \\\\
 			\frac {\partial f(\bm{x})}{\partial x_2} \\\\
 			\vdots \\\\
 			\frac {\partial f(\bm{x})}{\partial x_n}
 		\end{bmatrix}.\\[0.2cm]
\]
\end{definition}

\begin{definition} (Directional Derivative)
	Let $f:\mathbb{R}^n\to\mathbb{R}$ be a real valued function, and let $\bm{d}\in\mathbb{R}^n\backslash\{\bm{0}\}$. The directional derivative of $f$ in the direction $\bm{d}$ is defined as:
	
	\[\frac {\partial f}{\partial \bm{d}}(\bm{x}) = \lim_{h\to 0} \frac {f(\bm{x} + h\bm{d}) - f(\bm{x})}h.\\[0.2cm]\]
		
\end{definition}

\begin{definition} (Lipschitz Function) \cite{clarke} A function $f : S \to\mathbb{R}^m$ is called \textit{Lipschitz continuous} if there is a constant $L$ such that:

\[\norm{f(\bm{x}) - f(\bm{y})} \le L\cdot\norm{\bm{x} - \bm{y}},\ \forall\ \bm{x}, \bm{y}\in S.\\[0.2cm]\] 
	
\end{definition}

\begin{definition} 
	Let $U\subset \mathbb{R}^n$ be an open subset of a $\mathbb{R}^n$ and let $f : U\to\mathbb{R}$ be Lipschitz near $\bm{x}\in U$ and $\bm{d}\in \mathbb{R}^n$. The \textit{generalised directional derivative} \cite{clarke} of $f$ at $\bm{x}$ in the direction of $\bm{d}$ is given by:
	
	\[f^{\circ}(\bm{x};\bm{d}) = \lim_{\bm{y}\to\bm{x}}\sup_{h\downarrow 0} \frac{f(\bm{y} +h\bm{d}) - f(\bm{y})}h\]
\end{definition}

\begin{definition} (Clarke Gradient)
	The \textit{generalised gradient} of $f$ at $\bm{x}$, denoted by $\partial f(\bm{x})$ is the subset of $\mathbb{R}^n$ given by:
	
	\[\{A\in X^\star: f^{\circ}(\bm{x}; \bm{d})\ge A(\bm{d})\text{ for all }\bm{v}\in \mathbb{R}^n\}.\\[0.2cm]\]
\end{definition}

\noindent In addition, $\partial f(\bm{x})$ is a non-empty, convex and compact subset of $\mathbb{R}^n$ and for $\bm{d}\in \mathbb{R}^n$ we have:

\[f^{\circ}(\bm{x};\bm{d}) = \max\{A(\bm{d}) : A\in\partial f(\bm{x})\}\\[0.2cm]\]

\hfill

The following definitions, which are given in \cite{ref2}, provide the necessary background on the $\mathcal{L}$-derivative and the domain of Lipschitz maps, specialised to finite dimensions. We begin by introducing the following extension to Lipschitz functions:

\begin{definition} (Set-valued Lipschitz functions) The continuous function $f:U\to\mathbb{R}$ has a \textit{non-empty, convex and compact set-valued Lipschitz constant} $\bm{b}\in\mathbf{C}\mathbb{R}^n$ in an open subset $a\subset U$ if for all $\bm{x}, \bm{y}\in a$ we have: 

\[\bm{b}\cdot (\bm{x}-\bm{y}) \sqsubseteq f(\bm{x}) - f(\bm{y}).\\[0.2cm]\] 

\noindent The single step tie $\delta (a, \bm{b})\subseteq C^0(U)$ of $a$ with $\bm{b}$ is the collection of all partial functions $f$ on $U$ with $a\subset \text{dom}(f)\subset U$ in $C^0(U)$ which have $\bm{b}$ as non-empty convex compact set-values Lipschitz constant in $a$.
\end{definition}

For the rest of this subsection, we will assume that $n\ge 2$. It is possible, as shown in \cite{ref12} to give an equivalent definition for the Clarke gradient which is expressed only in terms of elementary set notions. In what follows, the equivalent generalised derivative for Lipschitz continuous maps will be referred to as $\mathcal{L}$-derivative. Also, from this points onwards we assume $n\ge 2$ and, in a slight abuse of notation, we will write $C^0$ instead of $C^0(U)$. 

\begin{definition}
	A \textit{step tie} of $C^0$ is any finite intersection $\bigcap_{i\in I} \delta(a_i, b_i)\subset C^0$, where $I$ is a finite indexing set. A \textit{tie} of $C^0$ is any intersection $\Delta = \bigcap_{i\in I} \delta(a_i, b_i)\subset C^0$, for an arbitrary index set $I$. The \textit{domain} of a non-empty tie $\Delta$ is defined as $\text{dom}(\Delta) = \bigcup_{i\in I} \{a_i\ |\ b_i\ne\bot\}$.
\end{definition}

\begin{definition}
	The \textit{primitive map} $\int : (U \to \mathbf{C}\mathbb{R}^n) \to T^1(U)$ is defined by $\int (g) = \bigcap_{i\in I} \delta(a_i, b_i)$, where $g = \sup_{i\in I} b_i\chi_{a_i}$. We usually write $\int (f)$ as $\int f$ and call it the set of \textit{primitives} of $f$.
\end{definition}

\begin{definition}
	A map $g : U\to\mathbf{C}\mathbb{R}^n$ is said to be integrable if $\int g\ne \emptyset$.
\end{definition}
 
Given a continuous function $f : U \to\mathbb{R}$, the relation $f ∈ \delta(a, \bm{b})$ gives finitary information about the local differential properties of $f$. The collection of all such information defines the $\mathcal{L}$-derivative of $f$, as follows:

\begin{definition}
	The \textit{derivative} of a continuous function $f : U\to\mathbb{R}$ is defined as:
	
	\[\mathcal{L}f = \bigsqcup_{f\in\delta (a, \bm{b})} b_{\chi_a} : U\to\mathbf{C}\mathbb{R}^n \\[0.2cm]\]
\end{definition}
 
\begin{theorem}
	In the case of arbitrary (possibly infinite) dimension, it holds that (see \cite{ref12}):
	\begin{enumerate}
		\item $\mathcal{L}f$ is well-defined and Scott continuous.
		\item If $f\in C^1$ then $\mathcal{L}f = f^\prime$ .
		\item $f\in \delta(a, \bm{b})$ if and only if $b\chi_a \sqsubseteq \mathcal{L}f$.
	\end{enumerate}	
\end{theorem}

The following interesting corollary holds, which can be seen as an alternative of the Fundamental Theorem of Calculus for domain-theoretic functions. Moreover, we can state the duality between Clarke gradient and the previously defined $\mathcal{L}$-derivative (see \cite{ref12}):

\begin{corollary}
	$f \in \int g \iff g \sqsubseteq \mathcal{L}f$.
\end{corollary}

\begin{theorem}
	In finite dimensional Euclidean spaces, the $\mathcal{L}$-derivative coincides with the Clarke gradient.
\end{theorem}

\section{The property of consistency}

In this section, we will give a formal definition to the concept of consistency which we briefly mentioned at the beginning of this chapter and is the central piece of the present work. We will also present the algorithms that can decide whether consistency holds given constraints on the function and derivative information, respectively. 
\\

Informally, a pair of functions $(f, g)$, representing function and derivative information respectively, is called \textit{consistent} if there exists a third function $w$ whose function value is approximated by the first component of the pair and whose derivative information is approximated by the second component.
\\

The function information $f$ is given by a finite set of \textit{step functions} represented as $\{a_i, b_i\}_{i\in I}$, where $a_i\subseteq\mathbb{R}^n$ is a rational hyper-rectangle and $b_i\subseteq \mathbb{R}$ is a compact interval such that $b_i$ and $b_j$ have non empty intersection whenever this is the case for the interiors of $a_i$ and $a_j$. Equivalently, we can write that $f\in D^0(U)$, i.e. $f$ is interval-valued within each hyper-rectangle defined in its domain. On the other hand, the derivative information $g$ is given for each of the $n$ partial derivatives in the form $\{a_i, b_i\}_{i\in I}$, where $a_i$ are as before, but we allow $b_i$ to be, in the most general setting, rational compact polyhedra. Simply put, $g : U \to\mathbf{C}\mathbb{R}^n$, i.e. the range of is a convex polyhedra for each $n$-dimensional hyper-rectangle inside the domain.
\\

As we will see, in the particular case when the derivative constraints lie within hyper-rectangles with faces parallel to the coordinate planes, we can decide whether a \textit{consistency witness} $w$ exists or not, via a linear programming algorithm. In addition, when the algorithm indicates that $(f, g)$ is consistent, we can also construct the least and greatest consistent maps -- in the sense that the heights of the minimal and maximal witnesses are minimised and maximised, respectively:

\[w_\min \le w\le w_\max.\\[0.2cm]\]

\noindent As further explained, the construction of such witnesses will give rise to \textit{piecewise linear} surfaces, because we approximate a witness by linear interpolating between the heights at the corners of any $n$-dimensional hyper-rectangle defined in the domain. Besides, the resulting $(n+1)$-dimensional surfaces (hyper-planes) will be kinked when transitioning between adjacent hyper-rectangles in the domain.
\\

It should be noted that choosing rectangular derivative information may result in some loss of information, but it will vastly simplify the LP algorithms and provides a practical framework for implementation. The more challenging problem of convex derivative information will also be analysed in a simple 2-dimensional setting, for which a linear test can also be developed.
\\

Let us now state the consistency relation more formally, as presented in \cite{ref2}, using the terminologies developed so far. Later, we will look at the algorithms that decide consistency, firstly in the rectangular setting for the 1D and $n$D cases, $n\ge 2$ (as shown in \cite{ref1} and \cite{ref2}, respectively).

\begin{definition}
	The consistency relation $\mathsf{Cons} \subset D^0(U)\times (U\to\mathbf{C}\mathbb{R}^n)$ is defined by:
	
	\[
		(f, g)\in\mathsf{Cons} \text{ if } \uparrow f\cap \textstyle\int g \ne \emptyset.\\[0.2cm]
	\]
\end{definition}

Also the least and greatest consistency witnesses can be defined through the following:

\begin{proposition} Let $O$ be a connected component $\text{dom}(g)$ and let $R(U)$ be the set of partial maps of $U$ into the extended real line $\mathbb{R}\cup \{-\infty, \infty\}$. Consider the two dcpos $(R(U), \le)$ and $(R(U), \ge)$ having pointwise ordering inherited from the extended real line. Then the maps $s : D^0(O) \times (U\to \mathbf{C}\mathbb{R}^n)\to (R(U), \le)$ and $t : D^0(O) \times (U\to \mathbf{C}\mathbb{R}^n)\to (R(U), \ge)$ defined by:

	\[
	\begin{array}{c}
		\displaystyle s(f, g) = \inf \left\{h:\text{dom}(g) \to\mathbb{R}\ \left|\ h\in\int g, h\ge f^-\right.\right\} \\\\
		\displaystyle t(f, g) = \sup \left\{h:\text{dom}(g) \to\mathbb{R}\ \left|\ h\in\int g, h\le f^+\right.\right\}
	\end{array}\\[0.2cm]
	\]
	
\noindent represent the least primitive map of $g$ that is greater than the lower part of $f$ and the greatest primitive map of $g$ that is less than the upper part of $f$, respectively. We also get from here that the following 3 conditions are equivalent:
	\begin{enumerate}
		\item $(f, g)\in \mathsf{Cons}$;
		\item $s(f, g) \le t(f, g)$;
		\item There exists a locally Lipschitz function $h:\text{dom}(g)\to\mathbb{R}$ with $g\sqsubseteq\mathcal{L}h$ and $f\sqsubseteq h$ on $\text{dom}(g)$.  
	\end{enumerate}

	Furthermore, the maps $s$ and $t$ are Scott continuous and the relation $\mathsf{Cons
}$ is Scott closed (see \cite{ref2}).

	\end{proposition}

\begin{corollary}
	Let $(f, g)\in \mathsf{Cons}$. Then in each connected component $O$ of the domain of definition of $g$ which intersects the domain of definition of $f$ there exist two locally Lipschitz functions $s(f, g) : O\to\mathbb{R}$ and $t(f, g) : O\to\mathbb{R}$ such that $s(f, g), t(f, g)\in\ \uparrow f \cap \int g$ and for each $w\in\ \uparrow f\cap \int g$, it holds that:
	
	\[s(f, g)(x) \le w(x) \le t(f, g)(x),\]
	
	\noindent for all $x\in O$.
\end{corollary}

\subsection{Consistency for the one-dimensional case}

We will now explain the framework for consistency in the case when the domain is unidimensional. For simplicity, we can assume, without loss of generality, that the domain $U = [0, 1]$. Otherwise, one can easily convert a closed interval $[a, b]$ into $[0, 1]$ upon a simple rescaling of the axis. Notice that the range of the derivative map $g$ will be $\mathbf{I}\mathbb{R}$, as we are only dealing with the derivative concerning a single variable. 
\\

Let us consider a partition $0 = y_0 < y_1 < \ldots < y_n = 1$ of the domain $[0, 1]$ in which interval-valued functions $f, g:[0, 1] \to \mathbf{I}\mathbb{R}$, representing approximations for function and derivative information, respectively, are given. Then the following algorithm (which appears in Section 3 from \cite{ref1}) with linear complexity in the number of partitions induced by $(f, g)$ can be developed as a test for consistency and also for determining the least map $s(f, g)\equiv w_\min$ which is also a witness to consistency. A analogous algorithm computes the greatest witness, i.e. $t(f, g)\equiv w_\max$.

\begin{algo} 
	The function updating algorithm consists of an initialisation step and two other main steps (see Figure \ref{1dAlgo}). The initialisation process determines the common partition points $\{y_0, \ldots y_n\}$ of $(f, g)$. On each interval $(y_{k-1}, y_k)$, the functions $g^-$ and $g^+$ are constant, with $g^-\vert_{(y_{k-1}, y_k)} = \lambda t.e_k^-$ and $g^+\vert_{(y_{k-1}, y_k)} = \lambda t.e_k^+$, where $e_k^-, e_k^+\in\mathbb{R}$. Furthermore, on each interval $(y_{k-1}, y_k)$, the map $f^-$ has a constant slope, $a_k$ say, i.e. $f^-\vert_{(y_{k-1}, y_k)} = f_k^-$, with $f_k^-(x) = a_kx + b_k$.
\end{algo}

\begin{figure}[H]
	\centerline{\includegraphics[scale=0.5]{pic8}}
	\caption{The function updating algorithm} 
	\label{1dAlgo}
\end{figure}

\vspace{2cm}

\noindent \textbf{Input:} $f, g :[0, 1]\to\mathbf{I}\mathbb{R}$, where $f$ is a linear step function and $g$ is a step function.
\\

\noindent \textbf{Output:} Continuous function $s(f, g) : [0, 1]\to\mathbf{I}\mathbb{R}$ which represents the least function consistent with the information from $f$ and $g$.
\\

\noindent \textbf{Initialisation:} 

\noindent $\{y_0,\ldots y_n\}$ \lit*{:  induced-partition-of} $(f, g)$ 
\\

\noindent \textbf{Part 1:}

\noindent $u(y_0) := f^-(y_0^+)$

\noindent \lit*{for } $k = 1\ldots n$ and $\forall x\in [y_{k-1}, y_k)$

\noindent \hspace{0.5cm} $u(x) := \max\{f^-(x), u(y_{k-1}) + (x - y_{k-1}) e_k^-\}$

\noindent \hspace{0.5cm} $u(y_k) := \max\{\overline{\lim} f^-(y_k), u(y_{k-1}) + (y_k - y_{k-1}) e_k^-\}$

\noindent \textbf{Part 2:}

\noindent $s(y_n) := u(y_n)$

\noindent \lit*{for } $k = n\ldots 1$ and $\forall x\in [y_{k-1}, y_k)$

\noindent \hspace{0.5cm} $s(f, g)(x) := \max\{u(x), s(y_k) + (x - y_{k-1}) e_k^+\}$

\begin{figure}[H]
\centerline{\includegraphics[scale=0.36]{pic7}}
\caption{Consistency for the 1D case. The least witness $s$ and the greatest witness to consistency $t$ are shown (they are both made up of piecewise-linear segments). Note that any other witness can vary within the shaded region.}
\label{1d}
\end{figure}

\subsection{Consistency for the $n$-dimensional case, $n\ge 2$} \label{nd-cons-section}

After getting a flavour of the algorithm in the one-dimensional setup, we can now proceed to the central piece of the present work by explaining the general framework that decides consistency in higher dimensional spaces. As with the previous case, we will assume, for convenience, that the derivative information is rectangular, that is, contained within hyper-rectangles (as we are referring to higher dimensional spaces). 
\\

We will begin our analysis with the case $n = 2$, as it provides a simpler way to reason about the geometric structure of the problem. The algorithm for this case is described in Section 3 from \cite{ref2}, but we will follow a much simpler geometric approach to derive the linear programming algorithm (without using theory of cones). Afterwards, it will become straightforward to generalise the algorithm for any arbitrary $n\ge 2$. 
\\

Let us therefore consider consider the two-dimensional case. We will assume, as in the uni-dimensional setting, that each variable is constrained to lie within the closed interval $[0, 1]$ and as such we can consider $U = [0, 1]^2$ as the domain of definition. Now, both the function approximation $f$ and derivative approximation $g$ will be defined for rational rectangles inside $U$, so we can therefore enforce a grid $(p_0, p_1, \ldots, p_{k-1})\times (q_0, q_1, \ldots, q_{l-1})$ within the unit square, where all points $p_i$ and $q_j$ lie on the $x$ and $y$ axis, respectively, $i=\overline{0, k - 1}$, $j=\overline{0, l - 1}$. By this construction, we have: $p_0 = q_0 = 0$ and $p_{k-1} = q_{l-1} = 1$. Then, for every sub-rectangle $R_{ij} = (p_i, p_{i+1})\times (q_j, q_{j+1})$ formed by adjacent grid points, the functions $f : U\to\mathbf{I}\mathbb{R}$ and $g : U\to\mathbf{I}\mathbb{R}^2$ are given as follows:
\begin{align}
\label{func-info} 
\restr{f}{R_{ij}} & = \left[c_{ij}^-, c_{ij}^+\right] = c_{ij} \in \mathbf{I}\mathbb{R} \\ \nonumber \\
\label{derivative-info}
	\restr{g}{R_{ij}} & = b_{ij}^1 \times b_{ij}^2 = b_{ij} \in \mathbf{I}\mathbb{R}^2,
\end{align}
$ $\\[-0.8cm]

\noindent for all $i=\overline{0, k - 2}$, $j=\overline{0, l - 2}$. Note that each $b_{ij}^k$, $k\in\{1, 2\}$ is also an interval, as we restrict the value of each partial derivative to lie within a closed and compact interval, so: $b_{ij}^k = \left[b_{ij}^{k-}, b_{ij}^{k+}\right]$ for every $k$.
\\

Therefore, the problem of checking consistency for the pair of step functions $(f, g)$ in the two-dimensional case reduces to determining the existence of heights $h_{ij}\equiv h((p_i, q_j))$ at all the grid points $(p_i, q_j)$, $i=\overline{0, k - 1}$, $j=\overline{0, l - 1}$, subject to the given constraints (\ref{func-info}) and (\ref{derivative-info}).

\definecolor{qqzzqq}{rgb}{0.,0.6,0.}

Let us first tackle the constraints for function information. By (\ref{func-info}) we see that all $4$ heights at the corners of any of the $k\times l$ sub-rectangles $R_{ij}$ must lie within $c_{ij}\in\mathbf{I}\mathbb{R}$ and hence:

\begin{equation}
\label{h-st}
c_{ij}^- \le h_{st} \le c_{ij}^+,\text{ for } s\in \{i, i+1\}, t\in\{j, j+1\}
\end{equation}
$ $\\[-0.7cm]

\noindent where $i=\overline{0, k-2}$, $j=\overline{0, l-2}$. An observant reader would now notice that most of the inequalities (\ref{h-st}) have the potential to overlap, due to the fact that a height $h_{ij}$ can de defined for at most $4$ sub-rectangles at a time (e.g. $h_{11}$ is constrained by $c_{00}$, $c_{10}$, $c_{01}$ as well as $c_{11}$). Therefore, we need to account for $3$ types of grid points $(p_i, q_j)$: \textcolor{red}{corner}, \textcolor{qqzzqq}{border} and \textcolor{blue}{interior} (see Figure \ref{grid-info} below), so that, upon intersecting all possible constraints, we can derive the tightest bounds for each height $h_{ij}$.
\\

\begin{figure}[H]
\centering
\definecolor{xdxdff}{rgb}{0.49019607843137253,0.49019607843137253,1.}
\definecolor{ffqqqq}{rgb}{1.,0.,0.}
\definecolor{qqqqff}{rgb}{0.,0.,1.}
\definecolor{cqcqcq}{rgb}{0.7529411764705882,0.7529411764705882,0.7529411764705882}
\begin{tikzpicture}[line cap=round,line join=round,>=triangle 45,x=1.0cm,y=1.0cm, scale=0.65]
\draw [color=cqcqcq,dash pattern=on 5pt off 5pt, xstep=2.0cm,ystep=2.0cm] (-1.1,-1) grid (20.1,18.1);

\draw[->,color=black] (0.,0.) -- (20.1,0.);
\draw[color=black] (19.446453960984204,0.11390969621757836) node [anchor=south west] {$x$};
\draw[->,color=black] (0.,0.) -- (0.,18.1);
\draw[color=black] (0.14238712027197378,17.278051347485054) node [anchor=south west] {$y$};

\clip(-1.1,-1) rectangle (20.1,18.1);
\draw (4.,2.)-- (0.,2.);
\draw (0.,4.)-- (4.,4.);
\draw (2.,4.)-- (2.,0.);
\draw (4.,4.)-- (4.,0.);
\draw (0.,14.)-- (2.,14.);
\draw (2.,14.)-- (2.,16.);
\draw (16.,2.)-- (16.,0.);
\draw (18.,16.)-- (18.,12.);
\draw (16.,14.)-- (18.,14.);
\draw (18.,4.)-- (18.,0.);
\draw (16.,2.)-- (18.,2.);
\draw (18.,14.)-- (18.,2.);
\draw (8.,12.)-- (12.,12.);
\draw (12.,8.)-- (8.,8.);
\draw (8.,12.)-- (8.,8.);
\draw (12.,12.)-- (12.,8.);
\draw (10.,12.)-- (10.,8.);
\draw (8.,10.)-- (12.,10.);
\draw (8.,2.)-- (12.,2.);
\draw (12.,2.)-- (12.,0.);
\draw (8.,2.)-- (8.,0.);
\draw (0.,16.)-- (18.,16.);
\draw (16.,16.)-- (16.,14.);
\draw (0.,16.)-- (0.,0.);
\draw (0.,0.)-- (18.,0.);
\draw (10.,2.)-- (10.,0.);
\begin{scriptsize}

\draw [fill=ffqqqq] (0, 0) circle (4.0pt);
\draw[color=ffqqqq] (0, 0.4) node [anchor=west] {$h_{00}$};
\draw[color=ffqqqq] (0, -0.4) node [anchor=west] {$p_0$};
\draw[color=ffqqqq] (0, 0.4) node [anchor=east] {$q_0$};

\draw [fill=qqzzqq] (2, 0) circle (4.0pt);
\draw[color=qqzzqq] (2, 0.4) node [anchor=west] {$h_{10}$};
\draw[color=qqzzqq] (2, -0.4) node [anchor=west] {$p_1$};

\draw [fill=qqzzqq] (4, 0) circle (4.0pt);
\draw[color=qqzzqq] (4, 0.4) node [anchor=west] {$h_{20}$};
\draw[color=qqzzqq] (4, -0.4) node [anchor=west] {$p_2$};

\draw [fill=qqzzqq] (0, 2) circle (4.0pt);
\draw[color=qqzzqq] (0, 2.4) node [anchor=west] {$h_{01}$};
\draw[color=qqzzqq] (0, 2.4) node [anchor=east] {$q_1$};

\draw [fill=qqqqff] (2, 2) circle (4.0pt);
\draw[color=qqqqff] (2, 2.4) node [anchor=west] {$h_{11}$};

\draw [fill=qqqqff] (4, 2) circle (4.0pt);
\draw[color=qqqqff] (4, 2.4) node [anchor=west] {$h_{21}$};

\draw [fill=qqzzqq] (0, 4) circle (4.0pt);
\draw[color=qqzzqq] (0, 4.4) node [anchor=west] {$h_{02}$};
\draw[color=qqzzqq] (0, 4.4) node [anchor=east] {$q_2$};

\draw [fill=qqqqff] (2, 4) circle (4.0pt);
\draw[color=qqqqff] (2, 4.4) node [anchor=west] {$h_{12}$};

\draw [fill=qqqqff] (4, 4) circle (4.0pt);
\draw[color=qqqqff] (4, 4.4) node [anchor=west] {$h_{22}$};

\draw [fill=ffqqqq] (0, 16) circle (4.0pt);
\draw[color=ffqqqq] (0, 16.4) node [anchor=west] {$h_{0,l-1}$};
\draw[color=ffqqqq] (0, 16.4) node [anchor=east] {$q_{l-1}$};

\draw [fill=qqzzqq] (0, 14) circle (4.0pt);
\draw[color=qqzzqq] (0, 14.4) node [anchor=west] {$h_{0,l-2}$};
\draw[color=qqzzqq] (0, 14.4) node [anchor=east] {$q_{l-2}$};

\draw [fill=qqqqff] (2, 14) circle (4.0pt);
\draw[color=qqqqff] (2, 14.4) node [anchor=west] {$h_{1,l-2}$};

\draw [fill=qqzzqq] (2, 16) circle (4.0pt);
\draw[color=qqzzqq] (2, 16.4) node [anchor=west] {$h_{1,l-1}$};

\draw [fill=qqqqff] (10, 10) circle (4.0pt);
\draw[color=qqqqff] (10, 10.4) node [anchor=west] {$h_{ij}$};

\draw [fill=ffqqqq] (18, 0) circle (4.0pt);
\draw[color=ffqqqq] (18, 0.4) node [anchor=west] {$h_{k-1,0}$};
\draw[color=ffqqqq] (18, -0.4) node [anchor=west] {$p_{k-1}$};

\draw [fill=ffqqqq] (18, 16) circle (4.0pt);
\draw[color=ffqqqq] (18, 16.4) node [anchor=west] {$h_{k-1,l-1}$};

\draw [fill=qqqqff] (16 , 14) circle (4.0pt);
\draw[color=qqqqff] (16, 14.4) node [anchor=west] {$h_{k-2, l-2}$};

\draw [fill=qqqqff] (16, 2) circle (4.0pt);
\draw[color=qqqqff] (16, 2.4) node [anchor=west] {$h_{k-2,1}$};

\draw [fill=qqzzqq] (16, 0) circle (4.0pt);
\draw[color=qqzzqq] (16, 0.4) node [anchor=west] {$h_{k-2, 0}$};
\draw[color=qqzzqq] (16, -0.4) node [anchor=west] {$p_{k-2}$};

\draw [fill=qqzzqq] (18, 14) circle (4.0pt);
\draw[color=qqzzqq] (18, 14.4) node [anchor=west] {$h_{k-1, l-2}$};

\draw [fill=qqzzqq] (18, 2) circle (4.0pt);
\draw[color=qqzzqq] (18, 2.4) node [anchor=west] {$h_{k-1, 1}$};

\draw [fill=qqqqff] (8, 12) circle (4.0pt);
\draw[color=qqqqff] (8, 12.4) node [anchor=west] {$h_{i-1, j+1}$};

\draw [fill=qqqqff] (12, 12) circle (4.0pt);
\draw[color=qqqqff] (12, 12.4) node [anchor=west] {$h_{i+1, j+1}$};

\draw [fill=qqqqff] (12, 8) circle (4.0pt);
\draw[color=qqqqff] (12, 8.4) node [anchor=west] {$h_{i+1, j-1}$};

\draw [fill=qqqqff] (8, 8) circle (4.0pt);
\draw[color=qqqqff] (8, 8.4) node [anchor=west] {$h_{i-1, j-1}$};

\draw [fill=qqqqff] (10, 12) circle (4.0pt);
\draw[color=qqqqff] (10, 12.4) node [anchor=west] {$h_{i, j+1}$};

\draw [fill=qqqqff] (10, 8.) circle (4.0pt);
\draw[color=qqqqff] (10, 8.4) node [anchor=west] {$h_{i, j-1}$};

\draw [fill=qqqqff] (8, 10) circle (4.0pt);
\draw[color=qqqqff] (8, 10.4) node [anchor=west] {$h_{i-1, j}$};

\draw [fill=qqqqff] (12, 10) circle (4.0pt);
\draw[color=qqqqff] (12, 10.4) node [anchor=west] {$h_{i+1, j}$};

\draw [fill=qqqqff] (8, 2) circle (4.0pt);
\draw[color=qqqqff] (8, 2.4) node [anchor=west] {$h_{i-1, 1}$};

\draw [fill=qqqqff] (12, 2) circle (4.0pt);
\draw[color=qqqqff] (12, 2.4) node [anchor=west] {$h_{i+1, 1}$};

\draw [fill=qqzzqq] (12, 0.) circle (4.0pt);
\draw[color=qqzzqq] (12, 0.4) node [anchor=west] {$h_{i+1,0}$};
\draw[color=qqzzqq] (12, -0.4) node [anchor=west] {$p_{i+1}$};

\draw [fill=qqzzqq] (8, 0) circle (4.0pt);
\draw[color=qqzzqq] (8, 0.4) node [anchor=west] {$h_{i-1,0}$};
\draw[color=qqzzqq] (8, -0.4) node [anchor=west] {$p_{i-1}$};

\draw [fill=qqqqff] (10, 2) circle (4.0pt);
\draw[color=qqqqff] (10, 2.4) node [anchor=west] {$h_{i1}$};

\draw [fill=qqzzqq] (16, 16) circle (4.0pt);
\draw[color=qqzzqq] (16, 16.4) node [anchor=west]  {$h_{k-2, l-1}$};

\draw [fill=qqzzqq] (10, 0) circle (4.0pt);
\draw[color=qqzzqq] (10, 0.4) node [anchor=west]  {$h_{i0}$};
\draw[color=qqzzqq] (10, -0.4) node [anchor=west] {$p_i$};

% c_ij numbers.
\draw[color=black] (1, 1)  node {$c_{00}$};
\draw[color=black] (3, 1)  node {$c_{10}$};
\draw[color=black] (1, 3)  node {$c_{01}$};
\draw[color=black] (3, 3)  node {$c_{11}$};
\draw[color=black] (1, 15) node {$c_{0,l-2}$};
\draw[color=black] (9, 1)  node {$c_{i-1,0}$};
\draw[color=black] (11, 1) node {$c_{i0}$};
\draw[color=black] (11, 9) node {$c_{i,j-1}$};
\draw[color=black] (9, 9)  node {$c_{i-1,j-1}$};
\draw[color=black] (9, 11) node {$c_{i-1,j}$};
\draw[color=black] (11, 11) node {$c_{ij}$};
\draw[color=black] (17, 15) node {$c_{k-2,l-2}$};
\draw[color=black] (17, 1) node {$c_{k-2,0}$};
\end{scriptsize}
\end{tikzpicture}
\caption{Function information for each sub-rectangle of the grid.}
\label{grid-info}
\end{figure}

\noindent We now consider each case in turn:

\begin{enumerate}
	\item If $(p_i, q_j)$ is a {\color{red}{corner}} point of the grid, then $h_{ij}$ will only use the constraint within its own sub-rectangle, namely:
	\begin{equation}
	\label{corner-ineq}
	\begin{gathered}
		\begin{array}{rcccl}
			c_{00}^-       & \le & h_{00}       & \le & c_{00}^+       \\\\
			c_{k-1, 0}^-   & \le & h_{k-1, 0}   & \le & c_{k-1, 0}^+   \\\\
			c_{0, l-1}^-   & \le & h_{0, l-1}   & \le & c_{0, l-1}^+   \\\\
			c_{k-1, l-1}^- & \le & h_{k-1, l-1} & \le & c_{k-1, l-1}^+.
		\end{array}
	\end{gathered}
	\end{equation}
	
	\item Now, if $(p_i, q_j)$ is a \textcolor{qqzzqq}{border} point of the grid, then, apart from the constraint defined for its sub-rectangle, $h_{ij}$ will also require the constraint from the adjacent sub-rectangle. For example, consider the grid point $(p_i, q_0)$ along the $x$ axis. Then the constraint for $h_{i0}$ will be derived by intersecting the constraints for the sub-rectangles $R_{i-1, 0}$ and $R_{i,0}$. We can proceed similarly for all the other border points along the sides of the unit-square to get the following inequalities:
	\begin{equation}
	\label{border-ineq}
	\begin{gathered}
		\begin{array}{rcccl}
			\max\left\{c_{i-1, 0}^-, c_{i0}^-\right\}       & \le & h_{i0}     & \le & \min\left\{c_{i-1, 0}^+, c_{i0}^+\right\}       \\\\
			\max\left\{c_{i-1, l-1}^-, c_{i, l-1}^-\right\} & \le & h_{i, l-1} & \le & \min\left\{c_{i-1, l-1}^+, c_{i, l-1}^+\right\} \\\\
			\max\left\{c_{0, j-1}^-, c_{0j}^-\right\}       & \le & h_{0j}     & \le & \min\left\{c_{0, j-1}^+, c_{0j}^+\right\}       \\\\
			\max\left\{c_{k-1, j-1}^-, c_{k-1, j}^-\right\} & \le & h_{k-1, j} & \le & \min\left\{c_{k-1, j-1}^+, c_{k-1, j}^+\right\},
		\end{array}
	\end{gathered}
	\end{equation}

	\noindent where $i=\overline{1, k-2}$ and $j=\overline{1, l-2}$.
	
	\item Finally, if $(p_i, q_j)$ is an \textcolor{blue}{interior} point of the grid, then in order to derive the best bounds for $h_{ij}$ we will need to account for all four adjacent sub-rectangles with common vertex $(p_i, q_j)$. As such, we can optimise the constraints in this case by taking the intersection of the intervals $c_{i-1, j-1}$, $c_{i, j-1}$, $c_{i-1, j}$ and $c_{ij}$, respectively, so we can thus write:
	\[	
		\max\left\{c_{i-1, j-1}^-, c_{i, j-1}^-, c_{i-1, j}^-, c_{ij}^-\right\} \le h_{ij} \le \min\left\{c_{i-1, j-1}^+, c_{i, j-1}^+, c_{i-1, j}^+, c_{ij}^+\right\}, \\
	\]
	\begin{equation}
	\label{interior-ineq}
	\begin{array}{l}
		\hspace*{-6.38cm} \text{for all } i=\overline{1, k-2} \text{ and } j=\overline{1, l-2}.
	\end{array}
	\end{equation}
	
\end{enumerate}

In the analysis developed so far, we made the tacit assumption that all adjacent intervals $c_{ij}$ have non-empty intersection. This is an absolute requirement for consistency to hold, because otherwise the rectangular parallelepipeds defined by adjacent sub-rectangles in the grid would not intersect and will naturally give rise to discontinuities. If that is the case, then the function information will be classified as being inconsistent. In equivalent terms, inconsistent function information will translate to at least one of the inequalities from above being false (in such a scenario the greatest lower bound for some $h_{ij}$ will be strictly greater than its lowest upper bound, contradiction). 
\\

Having broken down the function information constraints into inequalities of the type: $\verb|lower_bound|\le h_{ij}\le \verb|upper_bound|$, for all suitable indices $i$ and $j$, let us now focus our attention on the derivative constraints given by function $g$. 
\\

Consider an arbitrary sub-rectangle with lower left corner starting at $(p_i, q_j)$ in the square $U$, where $1\le i\le k-2$, $1\le j\le l-2$. The corresponding heights at each of the $4$ corners of the sub-rectangle $R_{ij}$ are, respectively, $h_{ij}$, $h_{i+1,j}$, $h_{i,j+1}$, $h_{i+1,j+1}$, all bounded within the rectangular box with height $c_{ij}^+ - c_{ij}^-$ (see the diagram from Figure \ref{piecewise-linear}). Since a consistent witness $h$ needs to pass through all the $4$ heights $h_{st}$, where $s\in\{i, i+1\}$, $t\in\{j, j+1\}$, it means that the resulting surface contained inside the rectangular parallelepiped will be piecewise-linear and can be obtained by interpolating the heights as follows:

\begin{itemize}
	\item In the triangle with vertices $(p_i, q_j)$, $(p_i, q_{j+1})$ and $(p_{i+1}, q_j)$, the map $h$ linearly interpolates between the values $h_{ij}$, $h_{i,j+1}$ and $h_{i+1,j}$ at these vertices respectively.
	\item In the triangle with vertices $(p_{i+1}, q_{j+1})$, $(p_i, q_{j+1})$ and $(p_{i+1}, q_j)$, the map $h$ linearly interpolates between the values $h_{i+1, j+1}$, $h_{i,j+1}$ and $h_{i+1,j}$ at these vertices respectively.
\end{itemize}

\noindent This construction is depicted in \ref{piecewise-linear}. Alternatively, the piecewise linear surface can be constructed by interpolating the values given at the vertices of the other two triangles, that is, one with vertices $(p_i, q_j)$, $(p_i, q_{j+1})$ and $(p_{i+1}, q_{j+1})$, and the other with vertices $(p_i , q_j)$, $(p_{i+1}, q_j)$ and $(p_{i+1} , q_{j+1})$ -- from a geometric perspective, this means that the slanted line of intersection between the two triangular surfaces would join $h_{ij}$ and $h_{i+1, j+1}$ as opposed to $h_{i, j+1}$ and $h_{i+1, j}$.

\begin{figure}[H]
\centering
\vspace{-1.8cm}
\definecolor{ffqqtt}{rgb}{1.,0.,0.2}
\definecolor{zzttqq}{rgb}{0.6,0.2,0.}
\definecolor{xdxdff}{rgb}{0.49019607843137253,0.49019607843137253,1.}
\definecolor{qqqqff}{rgb}{0.,0.,1.}
\definecolor{uuuuuu}{rgb}{0.26666666666666666,0.26666666666666666,0.26666666666666666}
\begin{tikzpicture}[line cap=round,line join=round,>=triangle 45,x=1.0cm,y=1.0cm, scale=1]
\clip(-0.6,-1) rectangle (14 ,13);
\fill[color=zzttqq,fill=zzttqq,pattern=north east lines,pattern color=zzttqq] (6.184386950551921,5.656551936950526) -- (7.407051402737721,8.556621621284727) -- (9.964386950551921,6.633559061719129) -- cycle;
\fill[color=zzttqq,fill=zzttqq,pattern=north east lines,pattern color=zzttqq] (7.407051402737721,8.556621621284727) -- (11.187051402737723,7.492774735097489) -- (9.964386950551921,6.633559061719129) -- cycle;
\fill[color=ffqqtt,fill=qqzzqq, fill opacity=0.5] (6.184386950551921,5.656551936950526) -- (9.964386950551921,5.656551936950523) -- (9.964386950551921,6.633559061719129) -- cycle;

\draw [dash pattern=on 2pt off 2pt] (7.407051402737721,4.237873618697574)-- (11.18705140273772,4.237873618697574);
\draw [dash pattern=on 2pt off 2pt] (11.18705140273772,4.237873618697574)-- (9.964386950551921,2.44313680814961);
\draw [dash pattern=on 2pt off 2pt] (9.964386950551921,2.44313680814961)-- (6.184386950551921,2.44313680814961);
\draw [dash pattern=on 2pt off 2pt] (7.407051402737721,4.237873618697574)-- (6.184386950551921,2.44313680814961);

\draw [dotted] (1.6643869505519218,2.44313680814961)-- (6.184386950551921,2.44313680814961);
\draw [dotted] (7.407051402737721,4.237873618697574)-- (2.8870514027377223,4.237873618697574);
\draw [dotted] (6.184386950551921,2.44313680814961)-- (4.52,0.);
\draw [dotted] (9.964386950551921,2.44313680814961)-- (8.3,0.);

\draw [color=qqqqff] (7.407051402737721,9.38)-- (11.18705140273772,9.38);
\draw [color=qqqqff] (6.184386950551922,7.585263189452039)-- (9.964386950551921,7.585263189452038);
\draw [color=qqqqff] (6.184386950551922,7.585263189452039)-- (7.407051402737721,9.38);
\draw [color=qqqqff] (9.964386950551921,7.585263189452038)-- (11.18705140273772,9.38);
\draw [color=qqqqff] (6.184386950551921,4.937671735849466)-- (7.407051402737721,6.732408546397429);
\draw [dotted,color=qqqqff] (7.407051402737721,6.732408546397429)-- (11.187051402737723,6.732408546397429);
\draw [color=qqqqff] (11.187051402737723,6.732408546397429)-- (9.96438695055192,4.937671735849462);
\draw [color=qqqqff] (9.96438695055192,4.937671735849462)-- (6.184386950551921,4.937671735849466);
\draw [color=qqqqff] (6.184386950551922,7.585263189452039)-- (6.184386950551921,4.937671735849466);
\draw [dotted,color=qqqqff] (7.407051402737721,9.38)-- (7.407051402737721,6.732408546397429);
\draw [color=qqqqff] (11.18705140273772,9.38)-- (11.187051402737723,6.732408546397429);
\draw [color=qqqqff] (9.964386950551921,7.585263189452038)-- (9.96438695055192,4.937671735849462);
\draw [dotted] (6.184386950551921,4.937671735849466)-- (6.184386950551921,2.44313680814961);
\draw [dotted] (7.407051402737721,6.732408546397429)-- (7.407051402737721,4.237873618697574);
\draw [dotted] (9.96438695055192,4.937671735849462)-- (9.964386950551921,2.44313680814961);
\draw [dotted] (11.187051402737723,6.732408546397429)-- (11.18705140273772,4.237873618697574);
\draw (6.184386950551921,5.656551936950526)-- (9.964386950551921,6.633559061719129);
\draw (9.964386950551921,6.633559061719129)-- (11.187051402737723,7.492774735097489);
\draw (11.187051402737723,7.492774735097489)-- (7.407051402737721,8.556621621284727);
\draw (6.184386950551921,5.656551936950526)-- (7.407051402737721,8.556621621284727);
\draw (7.407051402737721,8.556621621284727)-- (9.964386950551921,6.633559061719129);
\draw [color=zzttqq] (6.184386950551921,5.656551936950526)-- (7.407051402737721,8.556621621284727);
\draw [color=zzttqq] (7.407051402737721,8.556621621284727)-- (9.964386950551921,6.633559061719129);
\draw [color=zzttqq] (9.964386950551921,6.633559061719129)-- (6.184386950551921,5.656551936950526);
\draw [color=zzttqq] (7.407051402737721,8.556621621284727)-- (11.187051402737723,7.492774735097489);
\draw [color=zzttqq] (11.187051402737723,7.492774735097489)-- (9.964386950551921,6.633559061719129);
\draw [color=zzttqq] (9.964386950551921,6.633559061719129)-- (7.407051402737721,8.556621621284727);
\draw [dotted] (2.8870514027377223,9.38)-- (2.8870514027377223,4.237873618697574);
\draw [dotted] (1.664386950551922,7.58526318945204)-- (1.6643869505519218,2.44313680814961);
\draw [dotted] (7.407051402737721,9.38)-- (2.8870514027377223,9.38);
\draw [dotted] (6.184386950551922,7.585263189452039)-- (1.664386950551922,7.58526318945204);
\draw [dotted] (2.8870514027377223,9.38)-- (0.,5.142126381302433);
\draw [dotted] (2.8870514027377228,6.73240854639743)-- (0.,2.494534927699856);
\draw [dotted] (6.184386950551921,4.937671735849466)-- (1.6643869505519218,4.937671735849466);
\draw [dotted] (7.407051402737721,6.732408546397429)-- (2.8870514027377228,6.73240854639743);
\draw (6.184386950551921,5.656551936950526)-- (9.964386950551921,5.656551936950523);
\draw [color=ffqqtt] (6.184386950551921,5.656551936950526)-- (9.964386950551921,5.656551936950523);
\draw [color=ffqqtt] (9.964386950551921,5.656551936950523)-- (9.964386950551921,6.633559061719129);
\draw [color=ffqqtt] (9.964386950551921,6.633559061719129)-- (6.184386950551921,5.656551936950526);

% Axis.
%\draw [->] (0, 0) -- (7.678880847842544, 11.271751703255111);
\draw [->] (0, 0) -- (7.222766136228866, 10.602225521069895);
\draw [->] (0, 0) -- (0, 10.6);
\draw [->] (0, 0) -- (12.5, 0);

\begin{scriptsize}

% Origin.
\draw [fill=black] (0, 0) circle (1.5pt);
\draw[color=black] (0, 0) node [anchor=north east] {$\mbox{\normalsize\( O \)}$};

% p_i, p_(i+1), q_j, q_(j+1).
\draw [fill=black] (1.6643869505519218,2.44313680814961) circle (1.5pt);
\draw[color=black] (1.6643869505519218,2.44313680814961) node [anchor=north west] {$\mbox{\normalsize\( q_j \)}$};

\draw [fill=black] (2.8870514027377223,4.237873618697574) circle (1.5pt);
\draw[color=black] (2.8870514027377223,4.237873618697574) node [anchor=north west] {$\mbox{\normalsize\( q_{j+1} \)}$};

\draw [fill=black] (4.52, 0) circle (1.5pt);
\draw[color=black] (4.52, 0) node [anchor=north east] {$\mbox{\normalsize\( p_i \)}$};

\draw [fill=black] (8.3, 0) circle (1.5pt);
\draw[color=black] (8.3, 0) node [anchor=north east] {$\mbox{\normalsize\( p_{i+1} \)}$};

% Sub-rectangle.
\draw [fill=uuuuuu] (7.407051402737721,4.237873618697574) circle (1.5pt);
\draw[color=uuuuuu] (7.407051402737721,4.237873618697574) node [anchor=north east] {$\mbox{\normalsize\( (p_i, q_{j+1}) \)}$};

\draw [fill=uuuuuu] (11.18705140273772,4.237873618697574) circle (1.5pt);
\draw[color=uuuuuu] (11.18705140273772,4.237873618697574) node [anchor=north east] {$\mbox{\normalsize\( (p_{i+1}, q_{j+1}) \)}$};

\draw [fill=uuuuuu] (9.964386950551921, 2.44313680814961) circle (1.5pt);
\draw[color=uuuuuu] (9.964386950551921, 2.44313680814961) node [anchor=north east] {$\mbox{\normalsize\( (p_{i+1}, q_j) \)}$};

\draw [fill=uuuuuu] (6.184386950551921,2.44313680814961) circle (1.5pt);
\draw[color=uuuuuu] (6.184386950551921,2.44313680814961) node [anchor=north east] {$\mbox{\normalsize\( (p_i, q_j) \)}$};

% Box upper rectangle.
\draw [fill=black] (7.407051402737721,9.38) circle (1.5pt);
\draw [fill=black] (11.18705140273772,9.38) circle (1.5pt);
\draw [fill=black] (9.964386950551921,7.585263189452038) circle (1.5pt);
\draw [fill=black] (6.184386950551922,7.585263189452039) circle (1.5pt);

% Box lower rectangle.
\draw [fill=black] (7.407051402737721,6.732408546397429) circle (1.5pt);
\draw [fill=black] (11.187051402737723,6.732408546397429) circle (1.5pt);
\draw [fill=black] (6.184386950551921,4.937671735849466) circle (1.5pt);
\draw [fill=black] (9.96438695055192,4.937671735849462) circle (1.5pt);

% Heights.
\draw [fill=blue] (6.184386950551921,5.656551936950526) circle (2.5pt);
\draw[color=blue] (6.184386950551921,5.656551936950526) node [anchor=north east] {$\mbox{\normalsize\( h_{ij} \)}$};

\draw [fill=blue] (9.964386950551921,6.633559061719129) circle (2.5pt);
\draw[color=blue] (9.964386950551921,6.633559061719129) node [anchor=north west] {$\mbox{\normalsize\( h_{i+1, j} \)}$};

\draw [fill=blue] (11.187051402737723,7.492774735097489) circle (2.5pt);
\draw[color=blue] (11.187051402737723,7.492774735097489) node [anchor=south west] {$\mbox{\normalsize\( h_{i+1, j+1} \)}$};

\draw [fill=blue] (7.407051402737721,8.556621621284727) circle (2.5pt);
\draw[color=blue] (7.407051402737721,8.556621621284727) node [anchor=south west] {$\mbox{\normalsize\( h_{i, j+1} \)}$};

% Projection of left face.
\draw [fill=black] (2.8870514027377223,9.38) circle (1.5pt);
\draw [fill=black] (1.664386950551922,7.58526318945204) circle (1.5pt);
\draw [fill=black] (2.8870514027377228,6.73240854639743) circle (1.5pt);
\draw [fill=black] (1.6643869505519218,4.937671735849466) circle (1.5pt);

% C_ij- and c_ij+.
\draw [fill=uuuuuu] (0, 5.142126381302433) circle (1.5pt);
\draw[color=uuuuuu] (0, 5.142126381302433) node [anchor=east] {$\mbox{\normalsize\( c_{ij}^+ \)}$};

\draw [fill=uuuuuu] (0, 2.494534927699856) circle (1.5pt);
\draw[color=uuuuuu] (0, 2.494534927699856) node [anchor=east] {$\mbox{\normalsize\( c_{ij}^- \)}$};

% Right angle of red triangle.
\draw [fill=uuuuuu] (9.964386950551921,5.656551936950523) circle (1.5pt);

% Axis labels.
\draw[color=black] (12.8, 0) node {$\mbox{\normalsize\( x \)}$};
\draw[color=black] (7.4, 10.85) node {$\mbox{\normalsize\( y \)}$};
\draw[color=black] (0, 10.9) node {$\mbox{\normalsize\( z \)}$};

\end{scriptsize}
\end{tikzpicture}
\caption{Piecewise-linear surfaces obtained by interpolating along the heights of a consistent witness.}
\label{piecewise-linear}
\end{figure}

Hence, the piecewise linear map constructed above (in either scenario) satisfies the derivative constraints if and only if the slopes of each of the $4$ lines obtained by joining adjacent heights lie within the appropriate part of $b_{ij} = b_{ij}^1\times b_{ij}^2 \in \mathbf{I}\mathbb{R}^2$. For example, the slope of the line passing through the points $((p_i, q_j), h_{ij})$ and $((p_{i+1}, q_j), h_{i+1, j})$ must be included in the interval $b_{ij}^1$ corresponding to the partial derivative with respect to $x$ (see the \textcolor{qqzzqq}{green} triangle from Figure \ref{piecewise-linear}). Similarly, the slope of the line joining $((p_i, q_{j+1}), h_{i,j+1})$ and $((p_{i+1}, q_{j+1}), h_{i+1, j+1})$ must belong to $b_{ij}^1$, so we can thus write:
\begin{align}
\label{deriv-x}
	b_{ij}^{1-} \le \frac {h_{i+1,j} - h_{ij}}{p_{i+1} - p_i} \le b_{ij}^{1+} && b_{ij}^{1-} \le \frac {h_{i+1,j+1} - h_{i,j+1}}{p_{i+1} - p_i} \le b_{ij}^{1+},
\end{align}
$ $\\[-0.7cm]

\noindent where $i=\overline{0, k-2}$ and $j=\overline{0, l-1}$ (but in the second case we require that $j=\overline{0, l-2}$).

In a similar fashion, the constraints for the partial derivative with respect to $y$ translate to:
\begin{align}
\label{deriv-y}
	b_{ij}^{2-} \le \frac {h_{i,j+1} - h_{ij}}{q_{i+1} - q_i} \le b_{ij}^{2+} && b_{ij}^{2-} \le \frac {h_{i+1,j+1} - h_{i+1,j}}{q_{i+1} - q_i} \le b_{ij}^{2+},
\end{align}
$ $\\[-0.7cm]

\noindent where $i=\overline{0, k-1}$ and $j=\overline{0, l-2}$ (but $i=\overline{0, k-2}$ in the second case).
\\

We can already notice that the inequalities given by (\ref{deriv-x}) and (\ref{deriv-y}) feature the same kind of overlapping structure as with the inequalities (\ref{h-st}) defined within each sub-rectangle $R_{ij}$. As a result, there is a potential for optimising these constraints in the same manner as we did earlier for the function approximation (see Figure \ref{grid-info-deriv} below).
\\

\begin{figure}[H]
\centering
\definecolor{xdxdff}{rgb}{0.49019607843137253,0.49019607843137253,1.}
\definecolor{ffqqqq}{rgb}{1.,0.,0.}
\definecolor{qqqqff}{rgb}{0.,0.,1.}
\definecolor{cqcqcq}{rgb}{0.7529411764705882,0.7529411764705882,0.7529411764705882}
\begin{tikzpicture}[line cap=round,line join=round,>=triangle 45,x=1.0cm,y=1.0cm, scale=0.65]
\draw [color=cqcqcq,dash pattern=on 5pt off 5pt, xstep=2.0cm,ystep=2.0cm] (-1.1,-1) grid (20.1,18.1);

\draw[->,color=black] (0.,0.) -- (20.1,0.);
\draw[color=black] (19.446453960984204,0.11390969621757836) node [anchor=south west] {$x$};
\draw[->,color=black] (0.,0.) -- (0.,18.1);
\draw[color=black] (0.14238712027197378,17.278051347485054) node [anchor=south west] {$y$};

\clip(-1.1,-1) rectangle (20.1, 18.1);


\draw[color=blue] (4, 2) -- (0, 2);
\draw[color=blue] (0.,4.)-- (4.,4.);
\draw[color=blue] (2.,4.)-- (2.,0.);
\draw[color=blue] (4.,4.)-- (4.,0.);
\draw[color=blue] (0.,14.)-- (2.,14.);
\draw[color=blue] (2.,14.)-- (2.,16.);
\draw[color=blue] (16.,2.)-- (16.,0.);
\draw[color=blue] (18.,16.)-- (18.,12.);
\draw[color=blue] (16.,14.)-- (18.,14.);
\draw[color=blue] (18.,4.)-- (18.,0.);
\draw[color=blue] (16.,2.)-- (18.,2.);
\draw[color=blue] (18.,14.)-- (18.,2.);
\draw[color=blue] (8.,12.)-- (12.,12.);
\draw[color=blue] (12.,8.)-- (8.,8.);
\draw[color=blue] (8.,12.)-- (8.,8.);
\draw[color=blue] (12.,12.)-- (12.,8.);
\draw[color=blue] (10.,12.)-- (10.,8.);
\draw[color=blue] (8.,10.)-- (12.,10.);
\draw[color=blue] (8.,2.)-- (12.,2.);
\draw[color=blue] (12.,2.)-- (12.,0.);
\draw[color=blue] (8.,2.)-- (8.,0.);
\draw[color=blue] (0.,16.)-- (18.,16.);
\draw[color=blue] (16.,16.)-- (16.,14.);
\draw[color=blue] (0.,16.)-- (0.,0.);
\draw[color=blue] (0.,0.)-- (18.,0.);
\draw[color=blue] (10.,2.)-- (10.,0.);

% Green edges.
\draw[color=qqzzqq] (0, 0) -- (0, 16);
\draw[color=qqzzqq] (0, 0) -- (18, 0);
\draw[color=qqzzqq] (18, 0) -- (18, 16);
\draw[color=qqzzqq] (0, 16) -- (18, 16);
\begin{scriptsize}

\draw[fill=qqzzqq]  (0, 0) circle (4.0pt);
\draw[color=qqzzqq] (0, 0.4) node [anchor=west] {$h_{00}$};
\draw[color=qqzzqq] (0, -0.4) node [anchor=west] {$p_0$};
\draw[color=qqzzqq] (0, 0.4) node [anchor=east] {$q_0$};

\draw [fill=qqzzqq] (2, 0) circle (4.0pt);
\draw[color=qqzzqq] (2, 0.4) node [anchor=west] {$h_{10}$};
\draw[color=qqzzqq] (2, -0.4) node [anchor=west] {$p_1$};

\draw [fill=qqzzqq] (4, 0) circle (4.0pt);
\draw[color=qqzzqq] (4, 0.4) node [anchor=west] {$h_{20}$};
\draw[color=qqzzqq] (4, -0.4) node [anchor=west] {$p_2$};

\draw [fill=qqzzqq] (0, 2) circle (4.0pt);
\draw[color=qqzzqq] (0, 2.4) node [anchor=west] {$h_{01}$};
\draw[color=qqzzqq] (0, 2.4) node [anchor=east] {$q_1$};

\draw [fill=qqqqff] (2, 2) circle (4.0pt);
\draw[color=qqqqff] (2, 2.4) node [anchor=west] {$h_{11}$};

\draw [fill=qqqqff] (4, 2) circle (4.0pt);
\draw[color=qqqqff] (4, 2.4) node [anchor=west] {$h_{21}$};

\draw [fill=qqzzqq] (0, 4) circle (4.0pt);
\draw[color=qqzzqq] (0, 4.4) node [anchor=west] {$h_{02}$};
\draw[color=qqzzqq] (0, 4.4) node [anchor=east] {$q_2$};

\draw [fill=qqqqff] (2, 4) circle (4.0pt);
\draw[color=qqqqff] (2, 4.4) node [anchor=west] {$h_{12}$};

\draw [fill=qqqqff] (4, 4) circle (4.0pt);
\draw[color=qqqqff] (4, 4.4) node [anchor=west] {$h_{22}$};

\draw [fill=qqzzqq] (0, 16) circle (4.0pt);
\draw[color=qqzzqq] (0, 16.4) node [anchor=west] {$h_{0,l-1}$};
\draw[color=qqzzqq] (0, 16.4) node [anchor=east] {$q_{l-1}$};

\draw [fill=qqzzqq] (0, 14) circle (4.0pt);
\draw[color=qqzzqq] (0, 14.4) node [anchor=west] {$h_{0,l-2}$};
\draw[color=qqzzqq] (0, 14.4) node [anchor=east] {$q_{l-2}$};

\draw [fill=qqqqff] (2, 14) circle (4.0pt);
\draw[color=qqqqff] (2, 14.4) node [anchor=west] {$h_{1,l-2}$};

\draw [fill=qqzzqq] (2, 16) circle (4.0pt);
\draw[color=qqzzqq] (2, 16.4) node [anchor=west] {$h_{1,l-1}$};

\draw [fill=qqqqff] (10, 10) circle (4.0pt);
\draw[color=qqqqff] (10, 10.4) node [anchor=west] {$h_{ij}$};

\draw [fill=qqzzqq] (18, 0) circle (4.0pt);
\draw[color=qqzzqq] (18, 0.4) node [anchor=west] {$h_{k-1,0}$};
\draw[color=qqzzqq] (18, -0.4) node [anchor=west] {$p_{k-1}$};

\draw [fill=qqzzqq] (18, 16) circle (4.0pt);
\draw[color=qqzzqq] (18, 16.4) node [anchor=west] {$h_{k-1,l-1}$};

\draw [fill=qqqqff] (16 , 14) circle (4.0pt);
\draw[color=qqqqff] (16, 14.4) node [anchor=west] {$h_{k-2, l-2}$};

\draw [fill=qqqqff] (16, 2) circle (4.0pt);
\draw[color=qqqqff] (16, 2.4) node [anchor=west] {$h_{k-2,1}$};

\draw [fill=qqzzqq] (16, 0) circle (4.0pt);
\draw[color=qqzzqq] (16, 0.4) node [anchor=west] {$h_{k-2, 0}$};
\draw[color=qqzzqq] (16, -0.4) node [anchor=west] {$p_{k-2}$};

\draw [fill=qqzzqq] (18, 14) circle (4.0pt);
\draw[color=qqzzqq] (18, 14.4) node [anchor=west] {$h_{k-1, l-2}$};

\draw [fill=qqzzqq] (18, 2) circle (4.0pt);
\draw[color=qqzzqq] (18, 2.4) node [anchor=west] {$h_{k-1, 1}$};

\draw [fill=qqqqff] (8, 12) circle (4.0pt);
\draw[color=qqqqff] (8, 12.4) node [anchor=west] {$h_{i-1, j+1}$};

\draw [fill=qqqqff] (12, 12) circle (4.0pt);
\draw[color=qqqqff] (12, 12.4) node [anchor=west] {$h_{i+1, j+1}$};

\draw [fill=qqqqff] (12, 8) circle (4.0pt);
\draw[color=qqqqff] (12, 8.4) node [anchor=west] {$h_{i+1, j-1}$};

\draw [fill=qqqqff] (8, 8) circle (4.0pt);
\draw[color=qqqqff] (8, 8.4) node [anchor=west] {$h_{i-1, j-1}$};

\draw [fill=qqqqff] (10, 12) circle (4.0pt);
\draw[color=qqqqff] (10, 12.4) node [anchor=west] {$h_{i, j+1}$};

\draw [fill=qqqqff] (10, 8.) circle (4.0pt);
\draw[color=qqqqff] (10, 8.4) node [anchor=west] {$h_{i, j-1}$};

\draw [fill=qqqqff] (8, 10) circle (4.0pt);
\draw[color=qqqqff] (8, 10.4) node [anchor=west] {$h_{i-1, j}$};

\draw [fill=qqqqff] (12, 10) circle (4.0pt);
\draw[color=qqqqff] (12, 10.4) node [anchor=west] {$h_{i+1, j}$};

\draw [fill=qqqqff] (8, 2) circle (4.0pt);
\draw[color=qqqqff] (8, 2.4) node [anchor=west] {$h_{i-1, 1}$};

\draw [fill=qqqqff] (12, 2) circle (4.0pt);
\draw[color=qqqqff] (12, 2.4) node [anchor=west] {$h_{i+1, 1}$};

\draw [fill=qqzzqq] (12, 0.) circle (4.0pt);
\draw[color=qqzzqq] (12, 0.4) node [anchor=west] {$h_{i+1,0}$};
\draw[color=qqzzqq] (12, -0.4) node [anchor=west] {$p_{i+1}$};

\draw [fill=qqzzqq] (8, 0) circle (4.0pt);
\draw[color=qqzzqq] (8, 0.4) node [anchor=west] {$h_{i-1,0}$};
\draw[color=qqzzqq] (8, -0.4) node [anchor=west] {$p_{i-1}$};

\draw [fill=qqqqff] (10, 2) circle (4.0pt);
\draw[color=qqqqff] (10, 2.4) node [anchor=west] {$h_{i1}$};

\draw [fill=qqzzqq] (16, 16) circle (4.0pt);
\draw[color=qqzzqq] (16, 16.4) node [anchor=west]  {$h_{k-2, l-1}$};

\draw [fill=qqzzqq] (10, 0) circle (4.0pt);
\draw[color=qqzzqq] (10, 0.4) node [anchor=west]  {$h_{i0}$};
\draw[color=qqzzqq] (10, -0.4) node [anchor=west] {$p_i$};

% b_ij numbers.
\draw[color=black] (1, 1)  node {$b_{00}$};
\draw[color=black] (3, 1)  node {$b_{10}$};
\draw[color=black] (1, 3)  node {$b_{01}$};
\draw[color=black] (3, 3)  node {$b_{11}$};
\draw[color=black] (1, 15) node {$b_{0,l-2}$};
\draw[color=black] (9, 1)  node {$b_{i-1,0}$};
\draw[color=black] (11, 1) node {$b_{i0}$};
\draw[color=black] (11, 9) node {$b_{i,j-1}$};
\draw[color=black] (9, 9)  node {$b_{i-1,j-1}$};
\draw[color=black] (9, 11) node {$b_{i-1,j}$};
\draw[color=black] (11, 11) node {$b_{ij}$};
\draw[color=black] (17, 15) node {$b_{k-2,l-2}$};
\draw[color=black] (17, 1) node {$b_{k-2,0}$};
\end{scriptsize}
\end{tikzpicture}
\caption{Derivative information for each sub-rectangle of the grid.}
\label{grid-info-deriv}
\end{figure}

The difference here is that we only need to account for \textcolor{blue}{interior} sides of sub-rectangles, because the constraint coming from an adjacent sub-rectangle may potentially contribute towards improving the bounds of a slope of type $(h_{i+1,j}-h_{ij}) / (p_{i+1}-p_i)$. On the other hand, the slopes defined along the \textcolor{qqzzqq}{border} edges of the grid will only use the constraint within their corresponding sub-rectangle.
\\

By writing the inequalities in full and taking great care of indices, we get in each case that:

\begin{enumerate}
	\item The slopes along the \textcolor{qqzzqq}{border} sides of sub-rectangles are simply:
		\begin{equation}
		\label{border-deriv-ineq}
		\begin{gathered}
			\begin{array}{rcccl}
				b_{i0}^{1-}     & \le & \displaystyle \frac {h_{i+1, 0}   - h_{i0}    } {p_{i+1} - p_i} & \le & b_{i0}^{1+}     \\\\
				b_{i,l-2}^{1-}  & \le & \displaystyle \frac {h_{i+1, l-1} - h_{i, l-1}} {p_{i+1} - p_i} & \le & b_{i, l-2}^{1+} \\\\
		    	b_{0j}^{2-}     & \le & \displaystyle \frac {h_{0, j+1}   - h_{0j}    } {q_{j+1} - q_j} & \le & b_{0j}^{2+}     \\\\
				b_{k-2, j}^{2-} & \le & \displaystyle \frac {h_{k-1, j+1} - h_{k-1, j}} {q_{j+1} - q_j} & \le & b_{k-2, j}^{2+},
			\end{array}
		\end{gathered}
		\end{equation}
		
		\noindent where $i = \overline{0, k-2}$ and $j = \overline{0, l-2}$.
		
	\item The slopes along the \textcolor{blue}{interior} sides of sub-rectangles need also to take into account that the same slope occurs in an adjacent sub-rectangle and hence we get the following inequalities:
		\begin{equation}
		\label{interior-deriv-ineq}
		\begin{gathered}
			\begin{array}{rcccl}
				\max\left\{b_{i, j-1}^{1-}, b_{ij}^{1-}\right\} & \le & \displaystyle \frac {h_{i+1, j} - h_{ij}} {p_{i+1} - p_i} & \le & \min\left\{b_{i, j-1}^{1+}, b_{ij}^{1+}\right\} \\\\
				\max\left\{b_{i-1, j}^{2-}, b_{ij}^{2-}\right\} & \le & \displaystyle \frac {h_{i, j+1} - h_{ij}} {q_{j+1} - q_j} & \le & \min\left\{b_{i-1, j}^{2+}, b_{ij}^{2+}\right\},
			\end{array}
		\end{gathered}
		\end{equation}
		
		\noindent where $i = \overline{0, k-2}$ and $j = \overline{1, l-2}$ for the first chain of inequalities, while for the second we impose $i = \overline{1, k-2}$ and $j = \overline{0, l-2}$.
\end{enumerate}

From the discussion developed so far, we observe that all the inequalities given by (\ref{corner-ineq}), (\ref{border-ineq}), (\ref{interior-ineq}), (\ref{border-deriv-ineq}) and (\ref{interior-deriv-ineq}) represent, in fact, necessary and sufficient conditions for deciding whether a pair of functions $(f, g)$ is consistent or not. We can therefore introduce the following:

\begin{theorem}
	\label{decidability}
	(Decidability of consistency) A pair $(f, g)\in (U \to \mathbf{I}\mathbb{R}) \times (U \to \mathbf{I}\mathbb{R}^2)$, representing function and derivative information respectively, is consistent if and only if we can find values $h_{ij}\in\mathbb{R}$ at the grid points $(p_i, q_j)$ such that the inequalities {\rm{(\ref{corner-ineq})}}, {\rm{(\ref{border-ineq})}}, {\rm{(\ref{interior-ineq})}}, {\rm{(\ref{border-deriv-ineq})}} and {\rm{(\ref{interior-deriv-ineq})}} are simultaneously satisfied.
\end{theorem}

Since $f$ and $g$ are given in terms of rational numbers, the question of consistency boils down to solving a finite set of inequalities with rational coefficients for the $k\times l$ unknowns $h_{ij}$, which is decidable. In fact, it represents the set of constraints for a linear programming problem in which only a feasible solution is required to be found in order to guarantee the existence of a consistent witness.
\\

Finally, we will show that the construction of the minimal and maximal surfaces which are also witnesses to consistency can be derived upon minimising or maximising $\sum_{\substack{0\, \le\, i\, \le\, k-1 \\ 0\, \le\, j\, \le\, l-1}} h_{ij}$, respectively. Let us therefore present:
\\

\begin{algo} \label{minimal-witness}
	(Minimal consistent witness) For a consistent pair $(f, g)\in (U \to \mathbf{I}\mathbb{R}) \times (U \to \mathbf{I}\mathbb{R}^2)$, representing function and derivative information respectively, consider the following linear programming problem:
	\begin{equation}
	\label{LP-min}
	\begin{array}{ll}
		\rm{Minimise}   & \displaystyle\sum_{\substack{0\, \le\, i\, \le\, k-1 \\ 0\, \le\, j\, \le\, l-1}} h_{ij} \\\\
		\rm{subject\ to} & \rm{Constraints\ (\ref{corner-ineq}),\  (\ref{border-ineq}),\ (\ref{interior-ineq}),\ (\ref{border-deriv-ineq})\ and\ (\ref{interior-deriv-ineq}),}
	\end{array}	
	\end{equation} 
	
	\noindent which minimises the sum of heights at each of the given grid points $(p_i, q_j)$ in the unit square. Then the linear programming problem {\rm{(\ref{LP-min})}} has an optimal solution given by:
	\[w_{\min} = \left\{h_{ij}^\star\mid 0\le i\le k-1, 0\le j\le l-1\right\}\\[0.2cm]\]
	\noindent and furthermore $w_{\min}\le w$, for any other consistent witness $w$.
\end{algo}

\begin{proof} Let $W$ be the set of all witnesses consistent with the pair $(f, g)$. By Theorem \ref{decidability} we clearly have that $W\ne \emptyset$. According to Corollary 6.4 from \cite{ref2}, the minimal and maximal surfaces are guaranteed to exist for a consistent pair $(f, g)$. Thus, $w^\star = \inf W$ exists and is a witness. In particular, the heights $h_{ij}^\star$ of $w^\star$ are optimal and will therefore minimise the objective function. Now, if $w_{\min}\le w^\star$ then there would be some $h_{ij} < h_{ij}^\star$, thus contradicting the optimality of $h_{ij}^\star$ and also the optimality of the objective function. Consequently, $w^\star = w_{\min}$, as desired.
\end{proof}

A similar recipe applies for deriving the algorithm that will determine the greatest consistent witness $w_{\max}$, as we only need to convert (\ref{LP-min}) to become a maximisation problem:

\begin{algo} \label{maximal-witness}
	(Maximal consistent witness) For a consistent pair $(f, g)\in (U \to \mathbf{I}\mathbb{R}) \times (U \to \mathbf{I}\mathbb{R}^2)$, representing function and derivative information respectively, consider the following linear programming problem:
	\begin{equation}
	\label{LP-max}
	\begin{array}{ll}
		\rm{Maximise}    & \displaystyle\sum_{\substack{0\, \le\, i\, \le\, k-1 \\ 0\, \le\, j\, \le\, l-1}} h_{ij} \\\\
		\rm{subject\ to} & \rm{Constraints\ (\ref{corner-ineq}),\  (\ref{border-ineq}),\ (\ref{interior-ineq}),\ (\ref{border-deriv-ineq})\ and\ (\ref{interior-deriv-ineq}),}
	\end{array}	
	\end{equation} 
	
	\noindent which maximises the sum of heights at each of the given grid points $(p_i, q_j)$ in the unit square. Then the linear programming problem {\rm{(\ref{LP-min})}} has an optimal solution given by:
	\[w_{\max} = \left\{h_{ij}^\star\mid 0\le i\le k-1, 0\le j\le l-1\right\}\\[0.2cm]\]
	\noindent and furthermore $w\le w_{\max}$, for any other consistent witness $w$.\\
\end{algo}

The above algorithms (\ref{LP-min}) and (\ref{LP-max}), together with the Theorem of decidability of consistency \ref{decidability} can now be easily extended to an $n$-dimensional setting. We consider the unit cube $U\subset \mathbb{R}^n$ and a number of points along each of the $n$ edges intersecting at the origin of the cube so that $f\in\mathbf{I}\mathbb{R}$ and $g\in\mathbf{I}\mathbb{R}^n$ are constant in each of the resulting sub-hyper-rectangles. We then derive similar constraints to {\rm{(\ref{corner-ineq})}}, {\rm{(\ref{border-ineq})}} and {\rm{(\ref{interior-ineq})}} from the function approximation by intersecting all the overlapping intervals at grid-points contained in adjacent sub-hyper-rectangles. Similar to the inequalities {\rm{(\ref{border-deriv-ineq})}} and {\rm{(\ref{interior-deriv-ineq})}}, we need to account for each of the $n$ partial derivatives when constructing the derivative constraints -- the slope between heights at adjacent grid-points must lie within the intersection of all the intervals at sub-hyper-rectangles that contain the edges defined by neighbouring grid-points.
\\

Finally, the algorithms for minimising and maximising the least and greatest consistent witnesses are identical to the ones developed for the 2D case: we require that sum of heights for all grid-points to be either minimised or maximised, respectively, subject to the constraints derived from the pair .

\subsection{Consistency for a triangle with convex \\
derivative information in the two-dimensional case}

Let us now examine a more general framework for deciding consistency, which is presented in \cite{abbas}. Suppose we have a non-empty convex and compact polygon $B$ of the plane as the derivative information. Assume three distinct points $v_i = (x_i, y_i)$ with $i \in \{ 1, 2, 3 \}$ and $T_{123}$, or simply $T$, is the closed region defined by these three points. Suppose that we have $h_i \in \mathbb{R}$ for $i = \{1, 2, 3\}$. We aim to establish if there is a Lipschitz witness $z = w(x, y)$ with $w : T \rightarrow \mathbb{R}$, that goes through the three points $(v_i, h_i)$, for $i \in \{1, 2, 3\}$, whose derivative everywhere is contained in $B$. Let $z = P (x, y) = \alpha x + \beta y + \gamma$ with $P : T \rightarrow \mathbb{R}$ be the plane that goes through three points $(v_i, h_i)$, for $i \in\{1, 2, 3\}$ with $\nabla P = (\alpha, \beta) = b$. Then, $w$ satisfies our requirements if and only if $w - P$ goes through $(v_i , 0)$, for $i = \{1, 2, 3\}$, with its derivative consistent with $B - b$. Thus, we can equivalently consider the latter problem of finding a Lipschitz maps that goes through $(x_i, y_i , 0)$, for $i \in\{ 1, 2, 3 \}$, with its derivative contained everywhere in $B^{\prime} \equiv B - b$.
\\

Consider $e_{ij} \equiv v_j - v_i$ , with $ij$ in the cyclic order $1, 2, 3$ and let $e_{ij}^{\perp}$ be unit vector orthogonal to $e_{ij}$ in the direction into the triangle $T$. Now, by the mean value theorem (MVT) for Lipschitz maps, for a witness to exist, it is necessary that:
\[0 \in B^{\prime} \cdot e_{ij}\]
for all distinct pairs $i, j \in\{1, 2, 3\}$. Assume therefore that these conditions, called the MVT conditions, hold. Thus, there exist $b_{ij} \in B^{\prime}$ for cyclic ordered pairs $ij$ such that $b_{ij} \cdot e_{ij} = 0$, i.e., $b_{ij} = ke_{ij}^{\perp}$ for some $k \in \mathbb{R}$. Let $[b_{ij}^{-}, b_{ij}^{+}]$ be the interval along the direction $e_{ij}^{\perp}$ that is contained in $B$. Let $P_{ij}^{-}$ be the plane with $\nabla P_{ij}^{-} = b_{ij}^{-}$ that contains $e_{ij}$ for each cyclic order $ij$.

\begin{proposition} \cite{abbas} Suppose the MVT conditions hold. Then, there is Lipschitz map $w^\star : T \rightarrow \mathbb{R}$ that goes
through the three points $(x_i, y_i, h_i)$, with $i = 1, 2, 3$ and whose derivative is contained in $B$, such that for every other witness $w$ with these properties we have: 
\[\min w \leq w^\star \leq \max w.\]
\end{proposition} 

\begin{proof} If $b \in B$ (i.e. $0 \in B^{\prime}$) then the plane $w^\star = P$ satisfies our conditions; see the thick dashed triangle in the figure. Suppose, therefore, that $b \notin B$. By considering witnesses of the form $w - P$ we can equivalently consider the reduced problem with $h_i = 0$ and $B^{\prime} = B - b$. Then $0 \notin B^{\prime}$. In particular, the convex hull of the three segments $[b_{ij}^-, b_{ij}^+]$ along $e_{ij}^{\perp}$ does not contain $0$. This means that there exist two pairs $i_1 i_2$ and $i_3i_1$ such that precisely one of the following two conditions hold:

\begin{itemize}
  \item Both $b_{i_1i_2}^{-}$ and $b_{i_3i_1}^{-}$ have positive components along $e_{i_1i_2}^{\perp}$ and $e_{i_3i_1}^{\perp}$ respectively while $b_{i_2 i_3}^{+}$ has negative component on $e_{i_2i_3}^{+}$.
  \item Both $b_{i_1 i_2}^{-}$ and $b_{i_3 i_1}^{-}$ have negative components along $e_{i_1 i_2}^{\perp}$ and $e_{i_3 i_1}^{\perp}$ respectively while $b_{i_2 i_3}^{+}$ has positive component on $e_{i_2i_3}^{+}$.
\end{itemize}

\vspace{-0.3cm}

\begin{figure}[H]
    \centerline{\includegraphics[width=0.95\linewidth]{triangle}}
    \caption{Consistency for triangle}
    \label{triangle-convex}
\end{figure}

We consider the first case depicted in Figure \ref{triangle-convex} as the second is similar. Therefore, assume that there exist two pairs $i_1i_2$ and $i_3i_1$ such that $b_{i_1i_2}^{-}$ and $b_{i_3i_1}^{-}$ have positive components along $e_{i_1i_2}^{\perp}$ and $e_{i_3i_1}^{\perp}$ respectively while $b_{i_2i_3}^{+}$ has negative component on $e_{i_2i_3}^{\perp}$.
\\

Let $c_0 \equiv b_{i_3i_1}^-$, $c_1, \ldots, c_{k-1}, c_k \equiv b_{i_1i_2}^-$ be the vertices of $B^{\prime}$ from $b_{i_3 i_1}^{-}$ to $b_{i_1i_2}^{-}$ on the same side of the origin with respect to the line $l$ that goes through $b_{i_3i_1}^-$ and $b_{i_1 i_2}^-$. For any $c \in \mathbb{R}^2$ let $z = P_c(x, y) = c_1 x + c_2 y + \gamma_{c}$ be the plane through $(v_i, h_i)$ with $\nabla P_c = (c_1, c_2) = c$. Therefore, using our previous notation, we have $P_{c_0} = P_{i_3 i_1}$ and $P_{c_k} = P_{i_1i_2}$. Let $w^\star \equiv \min \left\{P_{c_j}^{-} : 0 \le j \le k \right\}$.
\\

Now let $w$ be any Lipschitz map through the three points $v_i = (x_i, y_i, 0)$ for $i \in\{ 1, 2, 3 \}$ and consistent with $B^{\prime}$. Then, $w$ is differentiable almost everywhere, and is equal to the integral of its derivative. Consider any path $p : [0, 1] \rightarrow T$, given by $p(t) = v_{i_1} + t\cdot (r - v_{i_1})$ from vertex $v_{i_1}$ to a point $r$ on the opposite edge $E$ of $T$. Then $\inf B^{\prime} \cdot (r - v_{i_1}) = c_j \cdot (r - v_{i_1})$ for some $j$ with $0 \le j \leq k$ and $\restr{w^\star}{[v_{i_1}, r]} = P_{c_j}$, where $[v_{i_1} , r]$ is the line segment from $v_{i_1}$ to $r$ in $T$. Thus we immediately obtain that:
\[
\begin{array}{lll}
	
w(r) - w(v_{i_1}) & = & \displaystyle \int_{0}^{1} w^{\prime} (v_{i_1} + t (r - v_{i_1})) \cdot (r - v_{i_1})\, \text{d}t  \\\\
                 & \ge & \displaystyle \int_0^1 c_j \cdot (r - v_{i_1})\, \text{d}t \\\\
                 & = & \displaystyle \int_{0}^{1} w^\star (v_{i_1} + t (r - v_{i_1})) \cdot (r - v_{i1})\, \text{d}t \\\\
                 & = & w^\star (r) - w^\star (v_{i_1})

\end{array}
\]

\noindent On the other hand, $w(v_{i_1}) = w^\star (v_{i_1}) = h_{i_1} = 0$. But $\max w^\star = \max_{r\in E} w^\star(r)$ and hence:

\[\max w \ge \max_{r\in E} w(r) \ge \max_{r\in E} w^\star (r) = \max w^\star\]

\noindent Consequently, $\max w \ge \max w^\star$ and since $\min w \le \min \left\{h_i \mid i \in \{1, 2, 3\}\right\} = \min w^\star$, the result follows.
\end{proof}

Now recall the original problem with three points $(x_i, y_i, h_i)$, for $i \in \{1, 2, 3\}$, and derivative information $B$. Suppose lower and upper limits $c^- \le c^+$ are given. Let $P : T \rightarrow \mathbb{R}$ be the plane that goes through three points $(x_i, y_i, h_i)$, for $i \in \{1, 2, 3\}$ with $\nabla P = b$ and put $B^{\prime} = B - b$. Consider $w^\star$ constructed above.

\begin{theorem} \cite{abbas} There is a witness to consistency if and only if the following two conditions hold:
\begin{itemize}
  \item \textit{For all distinct pairs $i, j \in\{1, 2, 3\}$, we have: $0 \in B^{\prime} \cdot e_{ij}$.}
  \item $c^- \leq w^\star + P \le c^+$.
\end{itemize}
\end{theorem}

This shows that consistency is semi-decidable, that is, for any given $h = (h_1, h_2, h_3) \in \mathbb{R}^3$ we can decide if there is a witness for consistency with heights $h_i$ at vertex $v_i$ for $i \in\{1, 2, 3\}$. Returning to the original problem we can now state the following:

\begin{theorem} \cite{abbas} Let rational points $v_i$ for $i \in\{1, 2, 3\}$, forming a triangle in the plane, a rational convex polygon $B$ and rational numbers $c^{-} \leq c^{+}$ be given. Then it is decidable that there exist heights $h_i$ for which there exists a Lipschitz witness going through $(v_i, h_i)$ consistent with $B$ and the bound $c^{-}$ and $c^{+}$ in the closed region bounded by the triangle, where $i\in\{1, 2, 3\}$.

\end{theorem}

\begin{proof} We first check if there exists $h \in \mathbb{R}^3$ such that the plane $z = P (x, y) = \alpha x + \beta y + \gamma$ going through the three points $(v_i, h_i)$ for $i\in\{1, 2, 3\}$ has gradient $(\alpha, \beta) \in B$. If this condition holds then $P$ is clearly a witness and we are done. Otherwise, we know from the construction presented in this section that there is a witness $w$ if and only if $w^\star + P$ is a witness, where $w^\star$ is the piecewise linear surface constructed above using $B^{\prime} = B - \nabla P$ and $z = P (x, y)$ is the plane passing through the three points $(v_i, h_i)$, with $h_i = w(v_i)$ where $i\in\{1, 2, 3\}$. Let us now show that for each consecutive pair of vertices $c_j$ and $c_{j+1}$ in the construction of $w^\star$ above, the \textit{slanted} line of intersection of the two planes $P_{c_j}$ and $P_{c_{j+1}}$ -- equivalently the line of intersection of the two planes $P_{j}^{\prime} \equiv P_{c_j} + P$ and $P_{j+1}^{\prime} \equiv P_{c_{j+1}} + P$ -- is perpendicular to the line segment $c_j c_{j+1}$. Recall that all planes $P_{c_j}$ pass through the point $(v, 0)$ with $v \equiv v_{i_1}$. This allows us to write that:

\[P_j (v) = P_{j + 1} (v) = 0.\\[0.3cm]\]

\noindent Now, for some $u = (x, y) \in T$ we can successively write that:

\[
\begin{array}{c}
	P_j^{\prime}(u) = P^{\prime}_{j + 1} (u) \iff P_{c_j} (u) = P_{c_{j + 1}}(u) \\\\
	\iff P_{c_j} (u) - P_{c_j} (v) = P_{c_{j + 1}} (u) - P_{j + 1} (v) \\\\
	\iff c_j \cdot (u - v) = c_{j + 1} \cdot (u - v) \\\\
	\iff (u - v) \cdot (c_j - c_{j + 1}) = 0,
\end{array}
\]

\noindent as claimed. It follows that the piecewise linear witness $w$ will linearly interpolate in each of the triangles with a vertex at $v_{i_1}$ and sides with common vertex $v_{i_1}$ perpendicular to the faces $c_j c_{j+1}$ for $j = \overline{0, k}$. This gives a simple triangulation of $T$. Let the vertices of the all the triangles in the triangulation be denoted by $u_0, u_1, \ldots, u_N$ which includes the vertices $v_1, v_2$ and $v_3$ and suppose $t_j \equiv w(u_j)$ for $j = \overline{1, N}$.
\end{proof}

