android / toolchain / cloog / 98972d5434ffcb4d11d2c81a46600e9a1cda9110 / . / cloog-0.17.0 / isl / doc / implementation.tex

\section{Sets and Relations} | |

\begin{definition}[Polyhedral Set] | |

A {\em polyhedral set}\index{polyhedral set} $S$ is a finite union of basic sets | |

$S = \bigcup_i S_i$, each of which can be represented using affine | |

constraints | |

$$ | |

S_i : \Z^n \to 2^{\Z^d} : \vec s \mapsto | |

S_i(\vec s) = | |

\{\, \vec x \in \Z^d \mid \exists \vec z \in \Z^e : | |

A \vec x + B \vec s + D \vec z + \vec c \geq \vec 0 \,\} | |

, | |

$$ | |

with $A \in \Z^{m \times d}$, | |

$B \in \Z^{m \times n}$, | |

$D \in \Z^{m \times e}$ | |

and $\vec c \in \Z^m$. | |

\end{definition} | |

\begin{definition}[Parameter Domain of a Set] | |

Let $S \in \Z^n \to 2^{\Z^d}$ be a set. | |

The {\em parameter domain} of $S$ is the set | |

$$\pdom S \coloneqq \{\, \vec s \in \Z^n \mid S(\vec s) \ne \emptyset \,\}.$$ | |

\end{definition} | |

\begin{definition}[Polyhedral Relation] | |

A {\em polyhedral relation}\index{polyhedral relation} | |

$R$ is a finite union of basic relations | |

$R = \bigcup_i R_i$ of type | |

$\Z^n \to 2^{\Z^{d_1+d_2}}$, | |

each of which can be represented using affine | |

constraints | |

$$ | |

R_i = \vec s \mapsto | |

R_i(\vec s) = | |

\{\, \vec x_1 \to \vec x_2 \in \Z^{d_1} \times \Z^{d_2} | |

\mid \exists \vec z \in \Z^e : | |

A_1 \vec x_1 + A_2 \vec x_2 + B \vec s + D \vec z + \vec c \geq \vec 0 \,\} | |

, | |

$$ | |

with $A_i \in \Z^{m \times d_i}$, | |

$B \in \Z^{m \times n}$, | |

$D \in \Z^{m \times e}$ | |

and $\vec c \in \Z^m$. | |

\end{definition} | |

\begin{definition}[Parameter Domain of a Relation] | |

Let $R \in \Z^n \to 2^{\Z^{d+d}}$ be a relation. | |

The {\em parameter domain} of $R$ is the set | |

$$\pdom R \coloneqq \{\, \vec s \in \Z^n \mid R(\vec s) \ne \emptyset \,\}.$$ | |

\end{definition} | |

\begin{definition}[Domain of a Relation] | |

Let $R \in \Z^n \to 2^{\Z^{d+d}}$ be a relation. | |

The {\em domain} of $R$ is the polyhedral set | |

$$\domain R \coloneqq \vec s \mapsto | |

\{\, \vec x_1 \in \Z^{d_1} \mid \exists \vec x_2 \in \Z^{d_2} : | |

(\vec x_1, \vec x_2) \in R(\vec s) \,\} | |

. | |

$$ | |

\end{definition} | |

\begin{definition}[Range of a Relation] | |

Let $R \in \Z^n \to 2^{\Z^{d+d}}$ be a relation. | |

The {\em range} of $R$ is the polyhedral set | |

$$ | |

\range R \coloneqq \vec s \mapsto | |

\{\, \vec x_2 \in \Z^{d_2} \mid \exists \vec x_1 \in \Z^{d_1} : | |

(\vec x_1, \vec x_2) \in R(\vec s) \,\} | |

. | |

$$ | |

\end{definition} | |

\begin{definition}[Composition of Relations] | |

Let $R \in \Z^n \to 2^{\Z^{d_1+d_2}}$ and | |

$S \in \Z^n \to 2^{\Z^{d_2+d_3}}$ be two relations, | |

then the composition of | |

$R$ and $S$ is defined as | |

$$ | |

S \circ R \coloneqq | |

\vec s \mapsto | |

\{\, \vec x_1 \to \vec x_3 \in \Z^{d_1} \times \Z^{d_3} | |

\mid \exists \vec x_2 \in \Z^{d_2} : | |

\vec x_1 \to \vec x_2 \in R(\vec s) \wedge | |

\vec x_2 \to \vec x_3 \in S(\vec s) | |

\,\} | |

. | |

$$ | |

\end{definition} | |

\begin{definition}[Difference Set of a Relation] | |

Let $R \in \Z^n \to 2^{\Z^{d+d}}$ be a relation. | |

The difference set ($\Delta \, R$) of $R$ is the set | |

of differences between image elements and the corresponding | |

domain elements, | |

$$ | |

\diff R \coloneqq | |

\vec s \mapsto | |

\{\, \vec \delta \in \Z^{d} \mid \exists \vec x \to \vec y \in R : | |

\vec \delta = \vec y - \vec x | |

\,\} | |

$$ | |

\end{definition} | |

\section{Simple Hull}\label{s:simple hull} | |

It is sometimes useful to have a single | |

basic set or basic relation that contains a given set or relation. | |

For rational sets, the obvious choice would be to compute the | |

(rational) convex hull. For integer sets, the obvious choice | |

would be the integer hull. | |

However, {\tt isl} currently does not support an integer hull operation | |

and even if it did, it would be fairly expensive to compute. | |

The convex hull operation is supported, but it is also fairly | |

expensive to compute given only an implicit representation. | |

Usually, it is not required to compute the exact integer hull, | |

and an overapproximation of this hull is sufficient. | |

The ``simple hull'' of a set is such an overapproximation | |

and it is defined as the (inclusion-wise) smallest basic set | |

that is described by constraints that are translates of | |

the constraints in the input set. | |

This means that the simple hull is relatively cheap to compute | |

and that the number of constraints in the simple hull is no | |

larger than the number of constraints in the input. | |

\begin{definition}[Simple Hull of a Set] | |

The {\em simple hull} of a set | |

$S = \bigcup_{1 \le i \le v} S_i$, with | |

$$ | |

S : \Z^n \to 2^{\Z^d} : \vec s \mapsto | |

S(\vec s) = | |

\left\{\, \vec x \in \Z^d \mid \exists \vec z \in \Z^e : | |

\bigvee_{1 \le i \le v} | |

A_i \vec x + B_i \vec s + D_i \vec z + \vec c_i \geq \vec 0 \,\right\} | |

$$ | |

is the set | |

$$ | |

H : \Z^n \to 2^{\Z^d} : \vec s \mapsto | |

S(\vec s) = | |

\left\{\, \vec x \in \Z^d \mid \exists \vec z \in \Z^e : | |

\bigwedge_{1 \le i \le v} | |

A_i \vec x + B_i \vec s + D_i \vec z + \vec c_i + \vec K_i \geq \vec 0 | |

\,\right\} | |

, | |

$$ | |

with $\vec K_i$ the (component-wise) smallest non-negative integer vectors | |

such that $S \subseteq H$. | |

\end{definition} | |

The $\vec K_i$ can be obtained by solving a number of | |

LP problems, one for each element of each $\vec K_i$. | |

If any LP problem is unbounded, then the corresponding constraint | |

is dropped. | |

\section{Parametric Integer Programming} | |

\subsection{Introduction}\label{s:intro} | |

Parametric integer programming \shortcite{Feautrier88parametric} | |

is used to solve many problems within the context of the polyhedral model. | |

Here, we are mainly interested in dependence analysis \shortcite{Fea91} | |

and in computing a unique representation for existentially quantified | |

variables. The latter operation has been used for counting elements | |

in sets involving such variables | |

\shortcite{BouletRe98,Verdoolaege2005experiences} and lies at the core | |

of the internal representation of {\tt isl}. | |

Parametric integer programming was first implemented in \texttt{PipLib}. | |

An alternative method for parametric integer programming | |

was later implemented in {\tt barvinok} \cite{barvinok-0.22}. | |

This method is not based on Feautrier's algorithm, but on rational | |

generating functions \cite{Woods2003short} and was inspired by the | |

``digging'' technique of \shortciteN{DeLoera2004Three} for solving | |

non-parametric integer programming problems. | |

In the following sections, we briefly recall the dual simplex | |

method combined with Gomory cuts and describe some extensions | |

and optimizations. The main algorithm is applied to a matrix | |

data structure known as a tableau. In case of parametric problems, | |

there are two tableaus, one for the main problem and one for | |

the constraints on the parameters, known as the context tableau. | |

The handling of the context tableau is described in \autoref{s:context}. | |

\subsection{The Dual Simplex Method} | |

Tableaus can be represented in several slightly different ways. | |

In {\tt isl}, the dual simplex method uses the same representation | |

as that used by its incremental LP solver based on the \emph{primal} | |

simplex method. The implementation of this LP solver is based | |

on that of {\tt Simplify} \shortcite{Detlefs2005simplify}, which, in turn, | |

was derived from the work of \shortciteN{Nelson1980phd}. | |

In the original \shortcite{Nelson1980phd}, the tableau was implemented | |

as a sparse matrix, but neither {\tt Simplify} nor the current | |

implementation of {\tt isl} does so. | |

Given some affine constraints on the variables, | |

$A \vec x + \vec b \ge \vec 0$, the tableau represents the relationship | |

between the variables $\vec x$ and non-negative variables | |

$\vec y = A \vec x + \vec b$ corresponding to the constraints. | |

The initial tableau contains $\begin{pmatrix} | |

\vec b & A | |

\end{pmatrix}$ and expresses the constraints $\vec y$ in the rows in terms | |

of the variables $\vec x$ in the columns. The main operation defined | |

on a tableau exchanges a column and a row variable and is called a pivot. | |

During this process, some coefficients may become rational. | |

As in the \texttt{PipLib} implementation, | |

{\tt isl} maintains a shared denominator per row. | |

The sample value of a tableau is one where each column variable is assigned | |

zero and each row variable is assigned the constant term of the row. | |

This sample value represents a valid solution if each constraint variable | |

is assigned a non-negative value, i.e., if the constant terms of | |

rows corresponding to constraints are all non-negative. | |

The dual simplex method starts from an initial sample value that | |

may be invalid, but that is known to be (lexicographically) no | |

greater than any solution, and gradually increments this sample value | |

through pivoting until a valid solution is obtained. | |

In particular, each pivot exchanges a row variable | |

$r = -n + \sum_i a_i \, c_i$ with negative | |

sample value $-n$ with a column variable $c_j$ | |

such that $a_j > 0$. Since $c_j = (n + r - \sum_{i\ne j} a_i \, c_i)/a_j$, | |

the new row variable will have a positive sample value $n$. | |

If no such column can be found, then the problem is infeasible. | |

By always choosing the column that leads to the (lexicographically) | |

smallest increment in the variables $\vec x$, | |

the first solution found is guaranteed to be the (lexicographically) | |

minimal solution \cite{Feautrier88parametric}. | |

In order to be able to determine the smallest increment, the tableau | |

is (implicitly) extended with extra rows defining the original | |

variables in terms of the column variables. | |

If we assume that all variables are non-negative, then we know | |

that the zero vector is no greater than the minimal solution and | |

then the initial extended tableau looks as follows. | |

$$ | |

\begin{tikzpicture} | |

\matrix (m) [matrix of math nodes] | |

{ | |

& {} & 1 & \vec c \\ | |

\vec x && |(top)| \vec 0 & I \\ | |

\vec r && \vec b & |(bottom)|A \\ | |

}; | |

\begin{pgfonlayer}{background} | |

\node (core) [inner sep=0pt,fill=black!20,right delimiter=),left delimiter=(,fit=(top)(bottom)] {}; | |

\end{pgfonlayer} | |

\end{tikzpicture} | |

$$ | |

Each column in this extended tableau is lexicographically positive | |

and will remain so because of the column choice explained above. | |

It is then clear that the value of $\vec x$ will increase in each step. | |

Note that there is no need to store the extra rows explicitly. | |

If a given $x_i$ is a column variable, then the corresponding row | |

is the unit vector $e_i$. If, on the other hand, it is a row variable, | |

then the row already appears somewhere else in the tableau. | |

In case of parametric problems, the sign of the constant term | |

may depend on the parameters. Each time the constant term of a constraint row | |

changes, we therefore need to check whether the new term can attain | |

negative and/or positive values over the current set of possible | |

parameter values, i.e., the context. | |

If all these terms can only attain non-negative values, the current | |

state of the tableau represents a solution. If one of the terms | |

can only attain non-positive values and is not identically zero, | |

the corresponding row can be pivoted. | |

Otherwise, we pick one of the terms that can attain both positive | |

and negative values and split the context into a part where | |

it only attains non-negative values and a part where it only attains | |

negative values. | |

\subsection{Gomory Cuts} | |

The solution found by the dual simplex method may have | |

non-integral coordinates. If so, some rational solutions | |

(including the current sample value), can be cut off by | |

applying a (parametric) Gomory cut. | |

Let $r = b(\vec p) + \sp {\vec a} {\vec c}$ be the row | |

corresponding to the first non-integral coordinate of $\vec x$, | |

with $b(\vec p)$ the constant term, an affine expression in the | |

parameters $\vec p$, i.e., $b(\vec p) = \sp {\vec f} {\vec p} + g$. | |

Note that only row variables can attain | |

non-integral values as the sample value of the column variables is zero. | |

Consider the expression | |

$b(\vec p) - \ceil{b(\vec p)} + \sp {\fract{\vec a}} {\vec c}$, | |

with $\ceil\cdot$ the ceiling function and $\fract\cdot$ the | |

fractional part. This expression is negative at the sample value | |

since $\vec c = \vec 0$ and $r = b(\vec p)$ is fractional, i.e., | |

$\ceil{b(\vec p)} > b(\vec p)$. On the other hand, for each integral | |

value of $r$ and $\vec c \ge 0$, the expression is non-negative | |

because $b(\vec p) - \ceil{b(\vec p)} > -1$. | |

Imposing this expression to be non-negative therefore does not | |

invalidate any integral solutions, while it does cut away the current | |

fractional sample value. To be able to formulate this constraint, | |

a new variable $q = \floor{-b(\vec p)} = - \ceil{b(\vec p)}$ is added | |

to the context. This integral variable is uniquely defined by the constraints | |

$0 \le -d \, b(\vec p) - d \, q \le d - 1$, with $d$ the common | |

denominator of $\vec f$ and $g$. In practice, the variable | |

$q' = \floor{\sp {\fract{-f}} {\vec p} + \fract{-g}}$ is used instead | |

and the coefficients of the new constraint are adjusted accordingly. | |

The sign of the constant term of this new constraint need not be determined | |

as it is non-positive by construction. | |

When several of these extra context variables are added, it is important | |

to avoid adding duplicates. | |

Recent versions of {\tt PipLib} also check for such duplicates. | |

\subsection{Negative Unknowns and Maximization} | |

There are two places in the above algorithm where the unknowns $\vec x$ | |

are assumed to be non-negative: the initial tableau starts from | |

sample value $\vec x = \vec 0$ and $\vec c$ is assumed to be non-negative | |

during the construction of Gomory cuts. | |

To deal with negative unknowns, \shortciteN[Appendix A.2]{Fea91} | |

proposed to use a ``big parameter'', say $M$, that is taken to be | |

an arbitrarily large positive number. Instead of looking for the | |

lexicographically minimal value of $\vec x$, we search instead | |

for the lexicographically minimal value of $\vec x' = \vec M + \vec x$. | |

The sample value $\vec x' = \vec 0$ of the initial tableau then | |

corresponds to $\vec x = -\vec M$, which is clearly not greater than | |

any potential solution. The sign of the constant term of a row | |

is determined lexicographically, with the coefficient of $M$ considered | |

first. That is, if the coefficient of $M$ is not zero, then its sign | |

is the sign of the entire term. Otherwise, the sign is determined | |

by the remaining affine expression in the parameters. | |

If the original problem has a bounded optimum, then the final sample | |

value will be of the form $\vec M + \vec v$ and the optimal value | |

of the original problem is then $\vec v$. | |

Maximization problems can be handled in a similar way by computing | |

the minimum of $\vec M - \vec x$. | |

When the optimum is unbounded, the optimal value computed for | |

the original problem will involve the big parameter. | |

In the original implementation of {\tt PipLib}, the big parameter could | |

even appear in some of the extra variables $\vec q$ created during | |

the application of a Gomory cut. The final result could then contain | |

implicit conditions on the big parameter through conditions on such | |

$\vec q$ variables. This problem was resolved in later versions | |

of {\tt PipLib} by taking $M$ to be divisible by any positive number. | |

The big parameter can then never appear in any $\vec q$ because | |

$\fract {\alpha M } = 0$. It should be noted, though, that an unbounded | |

problem usually (but not always) | |

indicates an incorrect formulation of the problem. | |

The original version of {\tt PipLib} required the user to ``manually'' | |

add a big parameter, perform the reformulation and interpret the result | |

\shortcite{Feautrier02}. Recent versions allow the user to simply | |

specify that the unknowns may be negative or that the maximum should | |

be computed and then these transformations are performed internally. | |

Although there are some application, e.g., | |

that of \shortciteN{Feautrier92multi}, | |

where it is useful to have explicit control over the big parameter, | |

negative unknowns and maximization are by far the most common applications | |

of the big parameter and we believe that the user should not be bothered | |

with such implementation issues. | |

The current version of {\tt isl} therefore does not | |

provide any interface for specifying big parameters. Instead, the user | |

can specify whether a maximum needs to be computed and no assumptions | |

are made on the sign of the unknowns. Instead, the sign of the unknowns | |

is checked internally and a big parameter is automatically introduced when | |

needed. For compatibility with {\tt PipLib}, the {\tt isl\_pip} tool | |

does explicitly add non-negativity constraints on the unknowns unless | |

the \verb+Urs_unknowns+ option is specified. | |

Currently, there is also no way in {\tt isl} of expressing a big | |

parameter in the output. Even though | |

{\tt isl} makes the same divisibility assumption on the big parameter | |

as recent versions of {\tt PipLib}, it will therefore eventually | |

produce an error if the problem turns out to be unbounded. | |

\subsection{Preprocessing} | |

In this section, we describe some transformations that are | |

or can be applied in advance to reduce the running time | |

of the actual dual simplex method with Gomory cuts. | |

\subsubsection{Feasibility Check and Detection of Equalities} | |

Experience with the original {\tt PipLib} has shown that Gomory cuts | |

do not perform very well on problems that are (non-obviously) empty, | |

i.e., problems with rational solutions, but no integer solutions. | |

In {\tt isl}, we therefore first perform a feasibility check on | |

the original problem considered as a non-parametric problem | |

over the combined space of unknowns and parameters. | |

In fact, we do not simply check the feasibility, but we also | |

check for implicit equalities among the integer points by computing | |

the integer affine hull. The algorithm used is the same as that | |

described in \autoref{s:GBR} below. | |

Computing the affine hull is fairly expensive, but it can | |

bring huge benefits if any equalities can be found or if the problem | |

turns out to be empty. | |

\subsubsection{Constraint Simplification} | |

If the coefficients of the unknown and parameters in a constraint | |

have a common factor, then this factor should be removed, possibly | |

rounding down the constant term. For example, the constraint | |

$2 x - 5 \ge 0$ should be simplified to $x - 3 \ge 0$. | |

{\tt isl} performs such simplifications on all sets and relations. | |

Recent versions of {\tt PipLib} also perform this simplification | |

on the input. | |

\subsubsection{Exploiting Equalities}\label{s:equalities} | |

If there are any (explicit) equalities in the input description, | |

{\tt PipLib} converts each into a pair of inequalities. | |

It is also possible to write $r$ equalities as $r+1$ inequalities | |

\shortcite{Feautrier02}, but it is even better to \emph{exploit} the | |

equalities to reduce the dimensionality of the problem. | |

Given an equality involving at least one unknown, we pivot | |

the row corresponding to the equality with the column corresponding | |

to the last unknown with non-zero coefficient. The new column variable | |

can then be removed completely because it is identically zero, | |

thereby reducing the dimensionality of the problem by one. | |

The last unknown is chosen to ensure that the columns of the initial | |

tableau remain lexicographically positive. In particular, if | |

the equality is of the form $b + \sum_{i \le j} a_i \, x_i = 0$ with | |

$a_j \ne 0$, then the (implicit) top rows of the initial tableau | |

are changed as follows | |

$$ | |

\begin{tikzpicture} | |

\matrix [matrix of math nodes] | |

{ | |

& {} & |(top)| 0 & I_1 & |(j)| & \\ | |

j && 0 & & 1 & \\ | |

&& 0 & & & |(bottom)|I_2 \\ | |

}; | |

\node[overlay,above=2mm of j,anchor=south]{j}; | |

\begin{pgfonlayer}{background} | |

\node (m) [inner sep=0pt,fill=black!20,right delimiter=),left delimiter=(,fit=(top)(bottom)] {}; | |

\end{pgfonlayer} | |

\begin{scope}[xshift=4cm] | |

\matrix [matrix of math nodes] | |

{ | |

& {} & |(top)| 0 & I_1 & \\ | |

j && |(left)| -b/a_j & -a_i/a_j & \\ | |

&& 0 & & |(bottom)|I_2 \\ | |

}; | |

\begin{pgfonlayer}{background} | |

\node (m2) [inner sep=0pt,fill=black!20,right delimiter=),left delimiter=(,fit=(top)(bottom)(left)] {}; | |

\end{pgfonlayer} | |

\end{scope} | |

\draw [shorten >=7mm,-to,thick,decorate, | |

decoration={snake,amplitude=.4mm,segment length=2mm, | |

pre=moveto,pre length=5mm,post length=8mm}] | |

(m) -- (m2); | |

\end{tikzpicture} | |

$$ | |

Currently, {\tt isl} also eliminates equalities involving only parameters | |

in a similar way, provided at least one of the coefficients is equal to one. | |

The application of parameter compression (see below) | |

would obviate the need for removing parametric equalities. | |

\subsubsection{Offline Symmetry Detection}\label{s:offline} | |

Some problems, notably those of \shortciteN{Bygde2010licentiate}, | |

have a collection of constraints, say | |

$b_i(\vec p) + \sp {\vec a} {\vec x} \ge 0$, | |

that only differ in their (parametric) constant terms. | |

These constant terms will be non-negative on different parts | |

of the context and this context may have to be split for each | |

of the constraints. In the worst case, the basic algorithm may | |

have to consider all possible orderings of the constant terms. | |

Instead, {\tt isl} introduces a new parameter, say $u$, and | |

replaces the collection of constraints by the single | |

constraint $u + \sp {\vec a} {\vec x} \ge 0$ along with | |

context constraints $u \le b_i(\vec p)$. | |

Any solution to the new system is also a solution | |

to the original system since | |

$\sp {\vec a} {\vec x} \ge -u \ge -b_i(\vec p)$. | |

Conversely, $m = \min_i b_i(\vec p)$ satisfies the constraints | |

on $u$ and therefore extends a solution to the new system. | |

It can also be plugged into a new solution. | |

See \autoref{s:post} for how this substitution is currently performed | |

in {\tt isl}. | |

The method described in this section can only detect symmetries | |

that are explicitly available in the input. | |

See \autoref{s:online} for the detection | |

and exploitation of symmetries that appear during the course of | |

the dual simplex method. | |

\subsubsection{Parameter Compression}\label{s:compression} | |

It may in some cases be apparent from the equalities in the problem | |

description that there can only be a solution for a sublattice | |

of the parameters. In such cases ``parameter compression'' | |

\shortcite{Meister2004PhD,Meister2008} can be used to replace | |

the parameters by alternative ``dense'' parameters. | |

For example, if there is a constraint $2x = n$, then the system | |

will only have solutions for even values of $n$ and $n$ can be replaced | |

by $2n'$. Similarly, the parameters $n$ and $m$ in a system with | |

the constraint $2n = 3m$ can be replaced by a single parameter $n'$ | |

with $n=3n'$ and $m=2n'$. | |

It is also possible to perform a similar compression on the unknowns, | |

but it would be more complicated as the compression would have to | |

preserve the lexicographical order. Moreover, due to our handling | |

of equalities described above there should be | |

no need for such variable compression. | |

Although parameter compression has been implemented in {\tt isl}, | |

it is currently not yet used during parametric integer programming. | |

\subsection{Postprocessing}\label{s:post} | |

The output of {\tt PipLib} is a quast (quasi-affine selection tree). | |

Each internal node in this tree corresponds to a split of the context | |

based on a parametric constant term in the main tableau with indeterminate | |

sign. Each of these nodes may introduce extra variables in the context | |

corresponding to integer divisions. Each leaf of the tree prescribes | |

the solution in that part of the context that satisfies all the conditions | |

on the path leading to the leaf. | |

Such a quast is a very economical way of representing the solution, but | |

it would not be suitable as the (only) internal representation of | |

sets and relations in {\tt isl}. Instead, {\tt isl} represents | |

the constraints of a set or relation in disjunctive normal form. | |

The result of a parametric integer programming problem is then also | |

converted to this internal representation. Unfortunately, the conversion | |

to disjunctive normal form can lead to an explosion of the size | |

of the representation. | |

In some cases, this overhead would have to be paid anyway in subsequent | |

operations, but in other cases, especially for outside users that just | |

want to solve parametric integer programming problems, we would like | |

to avoid this overhead in future. That is, we are planning on introducing | |

quasts or a related representation as one of several possible internal | |

representations and on allowing the output of {\tt isl\_pip} to optionally | |

be printed as a quast. | |

Currently, {\tt isl} also does not have an internal representation | |

for expressions such as $\min_i b_i(\vec p)$ from the offline | |

symmetry detection of \autoref{s:offline}. | |

Assume that one of these expressions has $n$ bounds $b_i(\vec p)$. | |

If the expression | |

does not appear in the affine expression describing the solution, | |

but only in the constraints, and if moreover, the expression | |

only appears with a positive coefficient, i.e., | |

$\min_i b_i(\vec p) \ge f_j(\vec p)$, then each of these constraints | |

can simply be reduplicated $n$ times, once for each of the bounds. | |

Otherwise, a conversion to disjunctive normal form | |

leads to $n$ cases, each described as $u = b_i(\vec p)$ with constraints | |

$b_i(\vec p) \le b_j(\vec p)$ for $j > i$ | |

and | |

$b_i(\vec p) < b_j(\vec p)$ for $j < i$. | |

Note that even though this conversion leads to a size increase | |

by a factor of $n$, not detecting the symmetry could lead to | |

an increase by a factor of $n!$ if all possible orderings end up being | |

considered. | |

\subsection{Context Tableau}\label{s:context} | |

The main operation that a context tableau needs to provide is a test | |

on the sign of an affine expression over the elements of the context. | |

This sign can be determined by solving two integer linear feasibility | |

problems, one with a constraint added to the context that enforces | |

the expression to be non-negative and one where the expression is | |

negative. As already mentioned by \shortciteN{Feautrier88parametric}, | |

any integer linear feasibility solver could be used, but the {\tt PipLib} | |

implementation uses a recursive call to the dual simplex with Gomory | |

cuts algorithm to determine the feasibility of a context. | |

In {\tt isl}, two ways of handling the context have been implemented, | |

one that performs the recursive call and one, used by default, that | |

uses generalized basis reduction. | |

We start with some optimizations that are shared between the two | |

implementations and then discuss additional details of each of them. | |

\subsubsection{Maintaining Witnesses}\label{s:witness} | |

A common feature of both integer linear feasibility solvers is that | |

they will not only say whether a set is empty or not, but if the set | |

is non-empty, they will also provide a \emph{witness} for this result, | |

i.e., a point that belongs to the set. By maintaining a list of such | |

witnesses, we can avoid many feasibility tests during the determination | |

of the signs of affine expressions. In particular, if the expression | |

evaluates to a positive number on some of these points and to a negative | |

number on some others, then no feasibility test needs to be performed. | |

If all the evaluations are non-negative, we only need to check for the | |

possibility of a negative value and similarly in case of all | |

non-positive evaluations. Finally, in the rare case that all points | |

evaluate to zero or at the start, when no points have been collected yet, | |

one or two feasibility tests need to be performed depending on the result | |

of the first test. | |

When a new constraint is added to the context, the points that | |

violate the constraint are temporarily removed. They are reconsidered | |

when we backtrack over the addition of the constraint, as they will | |

satisfy the negation of the constraint. It is only when we backtrack | |

over the addition of the points that they are finally removed completely. | |

When an extra integer division is added to the context, | |

the new coordinates of the | |

witnesses can easily be computed by evaluating the integer division. | |

The idea of keeping track of witnesses was first used in {\tt barvinok}. | |

\subsubsection{Choice of Constant Term on which to Split} | |

Recall that if there are no rows with a non-positive constant term, | |

but there are rows with an indeterminate sign, then the context | |

needs to be split along the constant term of one of these rows. | |

If there is more than one such row, then we need to choose which row | |

to split on first. {\tt PipLib} uses a heuristic based on the (absolute) | |

sizes of the coefficients. In particular, it takes the largest coefficient | |

of each row and then selects the row where this largest coefficient is smaller | |

than those of the other rows. | |

In {\tt isl}, we take that row for which non-negativity of its constant | |

term implies non-negativity of as many of the constant terms of the other | |

rows as possible. The intuition behind this heuristic is that on the | |

positive side, we will have fewer negative and indeterminate signs, | |

while on the negative side, we need to perform a pivot, which may | |

affect any number of rows meaning that the effect on the signs | |

is difficult to predict. This heuristic is of course much more | |

expensive to evaluate than the heuristic used by {\tt PipLib}. | |

More extensive tests are needed to evaluate whether the heuristic is worthwhile. | |

\subsubsection{Dual Simplex + Gomory Cuts} | |

When a new constraint is added to the context, the first steps | |

of the dual simplex method applied to this new context will be the same | |

or at least very similar to those taken on the original context, i.e., | |

before the constraint was added. In {\tt isl}, we therefore apply | |

the dual simplex method incrementally on the context and backtrack | |

to a previous state when a constraint is removed again. | |

An initial implementation that was never made public would also | |

keep the Gomory cuts, but the current implementation backtracks | |

to before the point where Gomory cuts are added before adding | |

an extra constraint to the context. | |

Keeping the Gomory cuts has the advantage that the sample value | |

is always an integer point and that this point may also satisfy | |

the new constraint. However, due to the technique of maintaining | |

witnesses explained above, | |

we would not perform a feasibility test in such cases and then | |

the previously added cuts may be redundant, possibly resulting | |

in an accumulation of a large number of cuts. | |

If the parameters may be negative, then the same big parameter trick | |

used in the main tableau is applied to the context. This big parameter | |

is of course unrelated to the big parameter from the main tableau. | |

Note that it is not a requirement for this parameter to be ``big'', | |

but it does allow for some code reuse in {\tt isl}. | |

In {\tt PipLib}, the extra parameter is not ``big'', but this may be because | |

the big parameter of the main tableau also appears | |

in the context tableau. | |

Finally, it was reported by \shortciteN{Galea2009personal}, who | |

worked on a parametric integer programming implementation | |

in {\tt PPL} \shortcite{PPL}, | |

that it is beneficial to add cuts for \emph{all} rational coordinates | |

in the context tableau. Based on this report, | |

the initial {\tt isl} implementation was adapted accordingly. | |

\subsubsection{Generalized Basis Reduction}\label{s:GBR} | |

The default algorithm used in {\tt isl} for feasibility checking | |

is generalized basis reduction \shortcite{Cook1991implementation}. | |

This algorithm is also used in the {\tt barvinok} implementation. | |

The algorithm is fairly robust, but it has some overhead. | |

We therefore try to avoid calling the algorithm in easy cases. | |

In particular, we incrementally keep track of points for which | |

the entire unit hypercube positioned at that point lies in the context. | |

This set is described by translates of the constraints of the context | |

and if (rationally) non-empty, any rational point | |

in the set can be rounded up to yield an integer point in the context. | |

A restriction of the algorithm is that it only works on bounded sets. | |

The affine hull of the recession cone therefore needs to be projected | |

out first. As soon as the algorithm is invoked, we then also | |

incrementally keep track of this recession cone. The reduced basis | |

found by one call of the algorithm is also reused as initial basis | |

for the next call. | |

Some problems lead to the | |

introduction of many integer divisions. Within a given context, | |

some of these integer divisions may be equal to each other, even | |

if the expressions are not identical, or they may be equal to some | |

affine combination of other variables. | |

To detect such cases, we compute the affine hull of the context | |

each time a new integer division is added. The algorithm used | |

for computing this affine hull is that of \shortciteN{Karr1976affine}, | |

while the points used in this algorithm are obtained by performing | |

integer feasibility checks on that part of the context outside | |

the current approximation of the affine hull. | |

The list of witnesses is used to construct an initial approximation | |

of the hull, while any extra points found during the construction | |

of the hull is added to this list. | |

Any equality found in this way that expresses an integer division | |

as an \emph{integer} affine combination of other variables is | |

propagated to the main tableau, where it is used to eliminate that | |

integer division. | |

\subsection{Experiments} | |

\autoref{t:comparison} compares the execution times of {\tt isl} | |

(with both types of context tableau) | |

on some more difficult instances to those of other tools, | |

run on an Intel Xeon W3520 @ 2.66GHz. | |

Easier problems such as the | |

test cases distributed with {\tt Pip\-Lib} can be solved so quickly | |

that we would only be measuring overhead such as input/output and conversions | |

and not the running time of the actual algorithm. | |

We compare the following versions: | |

{\tt piplib-1.4.0-5-g0132fd9}, | |

{\tt barvinok-0.32.1-73-gc5d7751}, | |

{\tt isl-0.05.1-82-g3a37260} | |

and {\tt PPL} version 0.11.2. | |

The first test case is the following dependence analysis problem | |

originating from the Phideo project \shortcite{Verhaegh1995PhD} | |

that was communicated to us by Bart Kienhuis: | |

\begin{lstlisting}[flexiblecolumns=true,breaklines=true]{} | |

lexmax { [j1,j2] -> [i1,i2,i3,i4,i5,i6,i7,i8,i9,i10] : 1 <= i1,j1 <= 8 and 1 <= i2,i3,i4,i5,i6,i7,i8,i9,i10 <= 2 and 1 <= j2 <= 128 and i1-1 = j1-1 and i2-1+2*i3-2+4*i4-4+8*i5-8+16*i6-16+32*i7-32+64*i8-64+128*i9-128+256*i10-256=3*j2-3+66 }; | |

\end{lstlisting} | |

This problem was the main inspiration | |

for some of the optimizations in \autoref{s:GBR}. | |

The second group of test cases are projections used during counting. | |

The first nine of these come from \shortciteN{Seghir2006minimizing}. | |

The remaining two come from \shortciteN{Verdoolaege2005experiences} and | |

were used to drive the first, Gomory cuts based, implementation | |

in {\tt isl}. | |

The third and final group of test cases are borrowed from | |

\shortciteN{Bygde2010licentiate} and inspired the offline symmetry detection | |

of \autoref{s:offline}. Without symmetry detection, the running times | |

are 11s and 5.9s. | |

All running times of {\tt barvinok} and {\tt isl} include a conversion | |

to disjunctive normal form. Without this conversion, the final two | |

cases can be solved in 0.07s and 0.21s. | |

The {\tt PipLib} implementation has some fixed limits and will | |

sometimes report the problem to be too complex (TC), while on some other | |

problems it will run out of memory (OOM). | |

The {\tt barvinok} implementation does not support problems | |

with a non-trivial lineality space (line) nor maximization problems (max). | |

The Gomory cuts based {\tt isl} implementation was terminated after 1000 | |

minutes on the first problem. The gbr version introduces some | |

overhead on some of the easier problems, but is overall the clear winner. | |

\begin{table} | |

\begin{center} | |

\begin{tabular}{lrrrrr} | |

& {\tt PipLib} & {\tt barvinok} & {\tt isl} cut & {\tt isl} gbr & {\tt PPL} \\ | |

\hline | |

\hline | |

% bart.pip | |

Phideo & TC & 793m & $>$999m & 2.7s & 372m \\ | |

\hline | |

e1 & 0.33s & 3.5s & 0.08s & 0.11s & 0.18s \\ | |

e3 & 0.14s & 0.13s & 0.10s & 0.10s & 0.17s \\ | |

e4 & 0.24s & 9.1s & 0.09s & 0.11s & 0.70s \\ | |

e5 & 0.12s & 6.0s & 0.06s & 0.14s & 0.17s \\ | |

e6 & 0.10s & 6.8s & 0.17s & 0.08s & 0.21s \\ | |

e7 & 0.03s & 0.27s & 0.04s & 0.04s & 0.03s \\ | |

e8 & 0.03s & 0.18s & 0.03s & 0.04s & 0.01s \\ | |

e9 & OOM & 70m & 2.6s & 0.94s & 22s \\ | |

vd & 0.04s & 0.10s & 0.03s & 0.03s & 0.03s \\ | |

bouleti & 0.25s & line & 0.06s & 0.06s & 0.15s \\ | |

difficult & OOM & 1.3s & 1.7s & 0.33s & 1.4s \\ | |

\hline | |

cnt/sum & TC & max & 2.2s & 2.2s & OOM \\ | |

jcomplex & TC & max & 3.7s & 3.9s & OOM \\ | |

\end{tabular} | |

\caption{Comparison of Execution Times} | |

\label{t:comparison} | |

\end{center} | |

\end{table} | |

\subsection{Online Symmetry Detection}\label{s:online} | |

Manual experiments on small instances of the problems of | |

\shortciteN{Bygde2010licentiate} and an analysis of the results | |

by the approximate MPA method developed by \shortciteN{Bygde2010licentiate} | |

have revealed that these problems contain many more symmetries | |

than can be detected using the offline method of \autoref{s:offline}. | |

In this section, we present an online detection mechanism that has | |

not been implemented yet, but that has shown promising results | |

in manual applications. | |

Let us first consider what happens when we do not perform offline | |

symmetry detection. At some point, one of the | |

$b_i(\vec p) + \sp {\vec a} {\vec x} \ge 0$ constraints, | |

say the $j$th constraint, appears as a column | |

variable, say $c_1$, while the other constraints are represented | |

as rows of the form $b_i(\vec p) - b_j(\vec p) + c$. | |

The context is then split according to the relative order of | |

$b_j(\vec p)$ and one of the remaining $b_i(\vec p)$. | |

The offline method avoids this split by replacing all $b_i(\vec p)$ | |

by a single newly introduced parameter that represents the minimum | |

of these $b_i(\vec p)$. | |

In the online method the split is similarly avoided by the introduction | |

of a new parameter. In particular, a new parameter is introduced | |

that represents | |

$\left| b_j(\vec p) - b_i(\vec p) \right|_+ = | |

\max(b_j(\vec p) - b_i(\vec p), 0)$. | |

In general, let $r = b(\vec p) + \sp {\vec a} {\vec c}$ be a row | |

of the tableau such that the sign of $b(\vec p)$ is indeterminate | |

and such that exactly one of the elements of $\vec a$ is a $1$, | |

while all remaining elements are non-positive. | |

That is, $r = b(\vec p) + c_j - f$ with $f = -\sum_{i\ne j} a_i c_i \ge 0$. | |

We introduce a new parameter $t$ with | |

context constraints $t \ge -b(\vec p)$ and $t \ge 0$ and replace | |

the column variable $c_j$ by $c' + t$. The row $r$ is now equal | |

to $b(\vec p) + t + c' - f$. The constant term of this row is always | |

non-negative because any negative value of $b(\vec p)$ is compensated | |

by $t \ge -b(\vec p)$ while and non-negative value remains non-negative | |

because $t \ge 0$. | |

We need to show that this transformation does not eliminate any valid | |

solutions and that it does not introduce any spurious solutions. | |

Given a valid solution for the original problem, we need to find | |

a non-negative value of $c'$ satisfying the constraints. | |

If $b(\vec p) \ge 0$, we can take $t = 0$ so that | |

$c' = c_j - t = c_j \ge 0$. | |

If $b(\vec p) < 0$, we can take $t = -b(\vec p)$. | |

Since $r = b(\vec p) + c_j - f \ge 0$ and $f \ge 0$, we have | |

$c' = c_j + b(\vec p) \ge 0$. | |

Note that these choices amount to plugging in | |

$t = \left|-b(\vec p)\right|_+ = \max(-b(\vec p), 0)$. | |

Conversely, given a solution to the new problem, we need to find | |

a non-negative value of $c_j$, but this is easy since $c_j = c' + t$ | |

and both of these are non-negative. | |

Plugging in $t = \max(-b(\vec p), 0)$ can be performed as in | |

\autoref{s:post}, but, as in the case of offline symmetry detection, | |

it may be better to provide a direct representation for such | |

expressions in the internal representation of sets and relations | |

or at least in a quast-like output format. | |

\section{Coalescing}\label{s:coalescing} | |

See \shortciteN{Verdoolaege2009isl}, for now. | |

More details will be added later. | |

\section{Transitive Closure} | |

\subsection{Introduction} | |

\begin{definition}[Power of a Relation] | |

Let $R \in \Z^n \to 2^{\Z^{d+d}}$ be a relation and | |

$k \in \Z_{\ge 1}$ | |

a positive number, then power $k$ of relation $R$ is defined as | |

\begin{equation} | |

\label{eq:transitive:power} | |

R^k \coloneqq | |

\begin{cases} | |

R & \text{if $k = 1$} | |

\\ | |

R \circ R^{k-1} & \text{if $k \ge 2$} | |

. | |

\end{cases} | |

\end{equation} | |

\end{definition} | |

\begin{definition}[Transitive Closure of a Relation] | |

Let $R \in \Z^n \to 2^{\Z^{d+d}}$ be a relation, | |

then the transitive closure $R^+$ of $R$ is the union | |

of all positive powers of $R$, | |

$$ | |

R^+ \coloneqq \bigcup_{k \ge 1} R^k | |

. | |

$$ | |

\end{definition} | |

Alternatively, the transitive closure may be defined | |

inductively as | |

\begin{equation} | |

\label{eq:transitive:inductive} | |

R^+ \coloneqq R \cup \left(R \circ R^+\right) | |

. | |

\end{equation} | |

Since the transitive closure of a polyhedral relation | |

may no longer be a polyhedral relation \shortcite{Kelly1996closure}, | |

we can, in the general case, only compute an approximation | |

of the transitive closure. | |

Whereas \shortciteN{Kelly1996closure} compute underapproximations, | |

we, like \shortciteN{Beletska2009}, compute overapproximations. | |

That is, given a relation $R$, we will compute a relation $T$ | |

such that $R^+ \subseteq T$. Of course, we want this approximation | |

to be as close as possible to the actual transitive closure | |

$R^+$ and we want to detect the cases where the approximation is | |

exact, i.e., where $T = R^+$. | |

For computing an approximation of the transitive closure of $R$, | |

we follow the same general strategy as \shortciteN{Beletska2009} | |

and first compute an approximation of $R^k$ for $k \ge 1$ and then project | |

out the parameter $k$ from the resulting relation. | |

\begin{example} | |

As a trivial example, consider the relation | |

$R = \{\, x \to x + 1 \,\}$. The $k$th power of this map | |

for arbitrary $k$ is | |

$$ | |

R^k = k \mapsto \{\, x \to x + k \mid k \ge 1 \,\} | |

. | |

$$ | |

The transitive closure is then | |

$$ | |

\begin{aligned} | |

R^+ & = \{\, x \to y \mid \exists k \in \Z_{\ge 1} : y = x + k \,\} | |

\\ | |

& = \{\, x \to y \mid y \ge x + 1 \,\} | |

. | |

\end{aligned} | |

$$ | |

\end{example} | |

\subsection{Computing an Approximation of $R^k$} | |

\label{s:power} | |

There are some special cases where the computation of $R^k$ is very easy. | |

One such case is that where $R$ does not compose with itself, | |

i.e., $R \circ R = \emptyset$ or $\domain R \cap \range R = \emptyset$. | |

In this case, $R^k$ is only non-empty for $k=1$ where it is equal | |

to $R$ itself. | |

In general, it is impossible to construct a closed form | |

of $R^k$ as a polyhedral relation. | |

We will therefore need to make some approximations. | |

As a first approximations, we will consider each of the basic | |

relations in $R$ as simply adding one or more offsets to a domain element | |

to arrive at an image element and ignore the fact that some of these | |

offsets may only be applied to some of the domain elements. | |

That is, we will only consider the difference set $\Delta\,R$ of the relation. | |

In particular, we will first construct a collection $P$ of paths | |

that move through | |

a total of $k$ offsets and then intersect domain and range of this | |

collection with those of $R$. | |

That is, | |

\begin{equation} | |

\label{eq:transitive:approx} | |

K = P \cap \left(\domain R \to \range R\right) | |

, | |

\end{equation} | |

with | |

\begin{equation} | |

\label{eq:transitive:path} | |

P = \vec s \mapsto \{\, \vec x \to \vec y \mid | |

\exists k_i \in \Z_{\ge 0}, \vec\delta_i \in k_i \, \Delta_i(\vec s) : | |

\vec y = \vec x + \sum_i \vec\delta_i | |

\wedge | |

\sum_i k_i = k > 0 | |

\,\} | |

\end{equation} | |

and with $\Delta_i$ the basic sets that compose | |

the difference set $\Delta\,R$. | |

Note that the number of basic sets $\Delta_i$ need not be | |

the same as the number of basic relations in $R$. | |

Also note that since addition is commutative, it does not | |

matter in which order we add the offsets and so we are allowed | |

to group them as we did in \eqref{eq:transitive:path}. | |

If all the $\Delta_i$s are singleton sets | |

$\Delta_i = \{\, \vec \delta_i \,\}$ with $\vec \delta_i \in \Z^d$, | |

then \eqref{eq:transitive:path} simplifies to | |

\begin{equation} | |

\label{eq:transitive:singleton} | |

P = \{\, \vec x \to \vec y \mid | |

\exists k_i \in \Z_{\ge 0} : | |

\vec y = \vec x + \sum_i k_i \, \vec \delta_i | |

\wedge | |

\sum_i k_i = k > 0 | |

\,\} | |

\end{equation} | |

and then the approximation computed in \eqref{eq:transitive:approx} | |

is essentially the same as that of \shortciteN{Beletska2009}. | |

If some of the $\Delta_i$s are not singleton sets or if | |

some of $\vec \delta_i$s are parametric, then we need | |

to resort to further approximations. | |

To ease both the exposition and the implementation, we will for | |

the remainder of this section work with extended offsets | |

$\Delta_i' = \Delta_i \times \{\, 1 \,\}$. | |

That is, each offset is extended with an extra coordinate that is | |

set equal to one. The paths constructed by summing such extended | |

offsets have the length encoded as the difference of their | |

final coordinates. The path $P'$ can then be decomposed into | |

paths $P_i'$, one for each $\Delta_i$, | |

\begin{equation} | |

\label{eq:transitive:decompose} | |

P' = \left( | |

(P_m' \cup \identity) \circ \cdots \circ | |

(P_2' \cup \identity) \circ | |

(P_1' \cup \identity) | |

\right) \cap | |

\{\, | |

\vec x' \to \vec y' \mid y_{d+1} - x_{d+1} = k > 0 | |

\,\} | |

, | |

\end{equation} | |

with | |

$$ | |

P_i' = \vec s \mapsto \{\, \vec x' \to \vec y' \mid | |

\exists k \in \Z_{\ge 1}, \vec \delta \in k \, \Delta_i'(\vec s) : | |

\vec y' = \vec x' + \vec \delta | |

\,\} | |

. | |

$$ | |

Note that each $P_i'$ contains paths of length at least one. | |

We therefore need to take the union with the identity relation | |

when composing the $P_i'$s to allow for paths that do not contain | |

any offsets from one or more $\Delta_i'$. | |

The path that consists of only identity relations is removed | |

by imposing the constraint $y_{d+1} - x_{d+1} > 0$. | |

Taking the union with the identity relation means that | |

that the relations we compose in \eqref{eq:transitive:decompose} | |

each consist of two basic relations. If there are $m$ | |

disjuncts in the input relation, then a direct application | |

of the composition operation may therefore result in a relation | |

with $2^m$ disjuncts, which is prohibitively expensive. | |

It is therefore crucial to apply coalescing (\autoref{s:coalescing}) | |

after each composition. | |

Let us now consider how to compute an overapproximation of $P_i'$. | |

Those that correspond to singleton $\Delta_i$s are grouped together | |

and handled as in \eqref{eq:transitive:singleton}. | |

Note that this is just an optimization. The procedure described | |

below would produce results that are at least as accurate. | |

For simplicity, we first assume that no constraint in $\Delta_i'$ | |

involves any existentially quantified variables. | |

We will return to existentially quantified variables at the end | |

of this section. | |

Without existentially quantified variables, we can classify | |

the constraints of $\Delta_i'$ as follows | |

\begin{enumerate} | |

\item non-parametric constraints | |

\begin{equation} | |

\label{eq:transitive:non-parametric} | |

A_1 \vec x + \vec c_1 \geq \vec 0 | |

\end{equation} | |

\item purely parametric constraints | |

\begin{equation} | |

\label{eq:transitive:parametric} | |

B_2 \vec s + \vec c_2 \geq \vec 0 | |

\end{equation} | |

\item negative mixed constraints | |

\begin{equation} | |

\label{eq:transitive:mixed} | |

A_3 \vec x + B_3 \vec s + \vec c_3 \geq \vec 0 | |

\end{equation} | |

such that for each row $j$ and for all $\vec s$, | |

$$ | |

\Delta_i'(\vec s) \cap | |

\{\, \vec \delta' \mid B_{3,j} \vec s + c_{3,j} > 0 \,\} | |

= \emptyset | |

$$ | |

\item positive mixed constraints | |

$$ | |

A_4 \vec x + B_4 \vec s + \vec c_4 \geq \vec 0 | |

$$ | |

such that for each row $j$, there is at least one $\vec s$ such that | |

$$ | |

\Delta_i'(\vec s) \cap | |

\{\, \vec \delta' \mid B_{4,j} \vec s + c_{4,j} > 0 \,\} | |

\ne \emptyset | |

$$ | |

\end{enumerate} | |

We will use the following approximation $Q_i$ for $P_i'$: | |

\begin{equation} | |

\label{eq:transitive:Q} | |

\begin{aligned} | |

Q_i = \vec s \mapsto | |

\{\, | |

\vec x' \to \vec y' | |

\mid {} & \exists k \in \Z_{\ge 1}, \vec f \in \Z^d : | |

\vec y' = \vec x' + (\vec f, k) | |

\wedge {} | |

\\ | |

& | |

A_1 \vec f + k \vec c_1 \geq \vec 0 | |

\wedge | |

B_2 \vec s + \vec c_2 \geq \vec 0 | |

\wedge | |

A_3 \vec f + B_3 \vec s + \vec c_3 \geq \vec 0 | |

\,\} | |

. | |

\end{aligned} | |

\end{equation} | |

To prove that $Q_i$ is indeed an overapproximation of $P_i'$, | |

we need to show that for every $\vec s \in \Z^n$, for every | |

$k \in \Z_{\ge 1}$ and for every $\vec f \in k \, \Delta_i(\vec s)$ | |

we have that | |

$(\vec f, k)$ satisfies the constraints in \eqref{eq:transitive:Q}. | |

If $\Delta_i(\vec s)$ is non-empty, then $\vec s$ must satisfy | |

the constraints in \eqref{eq:transitive:parametric}. | |

Each element $(\vec f, k) \in k \, \Delta_i'(\vec s)$ is a sum | |

of $k$ elements $(\vec f_j, 1)$ in $\Delta_i'(\vec s)$. | |

Each of these elements satisfies the constraints in | |

\eqref{eq:transitive:non-parametric}, i.e., | |

$$ | |

\left[ | |

\begin{matrix} | |

A_1 & \vec c_1 | |

\end{matrix} | |

\right] | |

\left[ | |

\begin{matrix} | |

\vec f_j \\ 1 | |

\end{matrix} | |

\right] | |

\ge \vec 0 | |

. | |

$$ | |

The sum of these elements therefore satisfies the same set of inequalities, | |

i.e., $A_1 \vec f + k \vec c_1 \geq \vec 0$. | |

Finally, the constraints in \eqref{eq:transitive:mixed} are such | |

that for any $\vec s$ in the parameter domain of $\Delta$, | |

we have $-\vec r(\vec s) \coloneqq B_3 \vec s + \vec c_3 \le \vec 0$, | |

i.e., $A_3 \vec f_j \ge \vec r(\vec s) \ge \vec 0$ | |

and therefore also $A_3 \vec f \ge \vec r(\vec s)$. | |

Note that if there are no mixed constraints and if the | |

rational relaxation of $\Delta_i(\vec s)$, i.e., | |

$\{\, \vec x \in \Q^d \mid A_1 \vec x + \vec c_1 \ge \vec 0\,\}$, | |

has integer vertices, then the approximation is exact, i.e., | |

$Q_i = P_i'$. In this case, the vertices of $\Delta'_i(\vec s)$ | |

generate the rational cone | |

$\{\, \vec x' \in \Q^{d+1} \mid \left[ | |

\begin{matrix} | |

A_1 & \vec c_1 | |

\end{matrix} | |

\right] \vec x' \,\}$ and therefore $\Delta'_i(\vec s)$ is | |

a Hilbert basis of this cone \shortcite[Theorem~16.4]{Schrijver1986}. | |

Note however that, as pointed out by \shortciteN{DeSmet2010personal}, | |

if there \emph{are} any mixed constraints, then the above procedure may | |

not compute the most accurate affine approximation of | |

$k \, \Delta_i(\vec s)$ with $k \ge 1$. | |

In particular, we only consider the negative mixed constraints that | |

happen to appear in the description of $\Delta_i(\vec s)$, while we | |

should instead consider \emph{all} valid such constraints. | |

It is also sufficient to consider those constraints because any | |

constraint that is valid for $k \, \Delta_i(\vec s)$ is also | |

valid for $1 \, \Delta_i(\vec s) = \Delta_i(\vec s)$. | |

Take therefore any constraint | |

$\spv a x + \spv b s + c \ge 0$ valid for $\Delta_i(\vec s)$. | |

This constraint is also valid for $k \, \Delta_i(\vec s)$ iff | |

$k \, \spv a x + \spv b s + c \ge 0$. | |

If $\spv b s + c$ can attain any positive value, then $\spv a x$ | |

may be negative for some elements of $\Delta_i(\vec s)$. | |

We then have $k \, \spv a x < \spv a x$ for $k > 1$ and so the constraint | |

is not valid for $k \, \Delta_i(\vec s)$. | |

We therefore need to impose $\spv b s + c \le 0$ for all values | |

of $\vec s$ such that $\Delta_i(\vec s)$ is non-empty, i.e., | |

$\vec b$ and $c$ need to be such that $- \spv b s - c \ge 0$ is a valid | |

constraint of $\Delta_i(\vec s)$. That is, $(\vec b, c)$ are the opposites | |

of the coefficients of a valid constraint of $\Delta_i(\vec s)$. | |

The approximation of $k \, \Delta_i(\vec s)$ can therefore be obtained | |

using three applications of Farkas' lemma. The first obtains the coefficients | |

of constraints valid for $\Delta_i(\vec s)$. The second obtains | |

the coefficients of constraints valid for the projection of $\Delta_i(\vec s)$ | |

onto the parameters. The opposite of the second set is then computed | |

and intersected with the first set. The result is the set of coefficients | |

of constraints valid for $k \, \Delta_i(\vec s)$. A final application | |

of Farkas' lemma is needed to obtain the approximation of | |

$k \, \Delta_i(\vec s)$ itself. | |

\begin{example} | |

Consider the relation | |

$$ | |

n \to \{\, (x, y) \to (1 + x, 1 - n + y) \mid n \ge 2 \,\} | |

. | |

$$ | |

The difference set of this relation is | |

$$ | |

\Delta = n \to \{\, (1, 1 - n) \mid n \ge 2 \,\} | |

. | |

$$ | |

Using our approach, we would only consider the mixed constraint | |

$y - 1 + n \ge 0$, leading to the following approximation of the | |

transitive closure: | |

$$ | |

n \to \{\, (x, y) \to (o_0, o_1) \mid n \ge 2 \wedge o_1 \le 1 - n + y \wedge o_0 \ge 1 + x \,\} | |

. | |

$$ | |

If, instead, we apply Farkas's lemma to $\Delta$, i.e., | |

\begin{verbatim} | |

D := [n] -> { [1, 1 - n] : n >= 2 }; | |

CD := coefficients D; | |

CD; | |

\end{verbatim} | |

we obtain | |

\begin{verbatim} | |

{ rat: coefficients[[c_cst, c_n] -> [i2, i3]] : i3 <= c_n and | |

i3 <= c_cst + 2c_n + i2 } | |

\end{verbatim} | |

The pure-parametric constraints valid for $\Delta$, | |

\begin{verbatim} | |

P := { [a,b] -> [] }(D); | |

CP := coefficients P; | |

CP; | |

\end{verbatim} | |

are | |

\begin{verbatim} | |

{ rat: coefficients[[c_cst, c_n] -> []] : c_n >= 0 and 2c_n >= -c_cst } | |

\end{verbatim} | |

Negating these coefficients and intersecting with \verb+CD+, | |

\begin{verbatim} | |

NCP := { rat: coefficients[[a,b] -> []] | |

-> coefficients[[-a,-b] -> []] }(CP); | |

CK := wrap((unwrap CD) * (dom (unwrap NCP))); | |

CK; | |

\end{verbatim} | |

we obtain | |

\begin{verbatim} | |

{ rat: [[c_cst, c_n] -> [i2, i3]] : i3 <= c_n and | |

i3 <= c_cst + 2c_n + i2 and c_n <= 0 and 2c_n <= -c_cst } | |

\end{verbatim} | |

The approximation for $k\,\Delta$, | |

\begin{verbatim} | |

K := solutions CK; | |

K; | |

\end{verbatim} | |

is then | |

\begin{verbatim} | |

[n] -> { rat: [i0, i1] : i1 <= -i0 and i0 >= 1 and i1 <= 2 - n - i0 } | |

\end{verbatim} | |

Finally, the computed approximation for $R^+$, | |

\begin{verbatim} | |

T := unwrap({ [dx,dy] -> [[x,y] -> [x+dx,y+dy]] }(K)); | |

R := [n] -> { [x,y] -> [x+1,y+1-n] : n >= 2 }; | |

T := T * ((dom R) -> (ran R)); | |

T; | |

\end{verbatim} | |

is | |

\begin{verbatim} | |

[n] -> { [x, y] -> [o0, o1] : o1 <= x + y - o0 and | |

o0 >= 1 + x and o1 <= 2 - n + x + y - o0 and n >= 2 } | |

\end{verbatim} | |

\end{example} | |

Existentially quantified variables can be handled by | |

classifying them into variables that are uniquely | |

determined by the parameters, variables that are independent | |

of the parameters and others. The first set can be treated | |

as parameters and the second as variables. Constraints involving | |

the other existentially quantified variables are removed. | |

\begin{example} | |

Consider the relation | |

$$ | |

R = | |

n \to \{\, x \to y \mid \exists \, \alpha_0, \alpha_1: 7\alpha_0 = -2 + n \wedge 5\alpha_1 = -1 - x + y \wedge y \ge 6 + x \,\} | |

. | |

$$ | |

The difference set of this relation is | |

$$ | |

\Delta = \Delta \, R = | |

n \to \{\, x \mid \exists \, \alpha_0, \alpha_1: 7\alpha_0 = -2 + n \wedge 5\alpha_1 = -1 + x \wedge x \ge 6 \,\} | |

. | |

$$ | |

The existentially quantified variables can be defined in terms | |

of the parameters and variables as | |

$$ | |

\alpha_0 = \floor{\frac{-2 + n}7} | |

\qquad | |

\text{and} | |

\qquad | |

\alpha_1 = \floor{\frac{-1 + x}5} | |

. | |

$$ | |

$\alpha_0$ can therefore be treated as a parameter, | |

while $\alpha_1$ can be treated as a variable. | |

This in turn means that $7\alpha_0 = -2 + n$ can be treated as | |

a purely parametric constraint, while the other two constraints are | |

non-parametric. | |

The corresponding $Q$~\eqref{eq:transitive:Q} is therefore | |

$$ | |

\begin{aligned} | |

n \to \{\, (x,z) \to (y,w) \mid | |

\exists\, \alpha_0, \alpha_1, k, f : {} & | |

k \ge 1 \wedge | |

y = x + f \wedge | |

w = z + k \wedge {} \\ | |

& | |

7\alpha_0 = -2 + n \wedge | |

5\alpha_1 = -k + x \wedge | |

x \ge 6 k | |

\,\} | |

. | |

\end{aligned} | |

$$ | |

Projecting out the final coordinates encoding the length of the paths, | |

results in the exact transitive closure | |

$$ | |

R^+ = | |

n \to \{\, x \to y \mid \exists \, \alpha_0, \alpha_1: 7\alpha_1 = -2 + n \wedge 6\alpha_0 \ge -x + y \wedge 5\alpha_0 \le -1 - x + y \,\} | |

. | |

$$ | |

\end{example} | |

The fact that we ignore some impure constraints clearly leads | |

to a loss of accuracy. In some cases, some of this loss can be recovered | |

by not considering the parameters in a special way. | |

That is, instead of considering the set | |

$$ | |

\Delta = \diff R = | |

\vec s \mapsto | |

\{\, \vec \delta \in \Z^{d} \mid \exists \vec x \to \vec y \in R : | |

\vec \delta = \vec y - \vec x | |

\,\} | |

$$ | |

we consider the set | |

$$ | |

\Delta' = \diff R' = | |

\{\, \vec \delta \in \Z^{n+d} \mid \exists | |

(\vec s, \vec x) \to (\vec s, \vec y) \in R' : | |

\vec \delta = (\vec s - \vec s, \vec y - \vec x) | |

\,\} | |

. | |

$$ | |

The first $n$ coordinates of every element in $\Delta'$ are zero. | |

Projecting out these zero coordinates from $\Delta'$ is equivalent | |

to projecting out the parameters in $\Delta$. | |

The result is obviously a superset of $\Delta$, but all its constraints | |

are of type \eqref{eq:transitive:non-parametric} and they can therefore | |

all be used in the construction of $Q_i$. | |

\begin{example} | |

Consider the relation | |

$$ | |

% [n] -> { [x, y] -> [1 + x, 1 - n + y] | n >= 2 } | |

R = n \to \{\, (x, y) \to (1 + x, 1 - n + y) \mid n \ge 2 \,\} | |

. | |

$$ | |

We have | |

$$ | |

\diff R = n \to \{\, (1, 1 - n) \mid n \ge 2 \,\} | |

$$ | |

and so, by treating the parameters in a special way, we obtain | |

the following approximation for $R^+$: | |

$$ | |

n \to \{\, (x, y) \to (x', y') \mid n \ge 2 \wedge y' \le 1 - n + y \wedge x' \ge 1 + x \,\} | |

. | |

$$ | |

If we consider instead | |

$$ | |

R' = \{\, (n, x, y) \to (n, 1 + x, 1 - n + y) \mid n \ge 2 \,\} | |

$$ | |

then | |

$$ | |

\diff R' = \{\, (0, 1, y) \mid y \le -1 \,\} | |

$$ | |

and we obtain the approximation | |

$$ | |

n \to \{\, (x, y) \to (x', y') \mid n \ge 2 \wedge x' \ge 1 + x \wedge y' \le x + y - x' \,\} | |

. | |

$$ | |

If we consider both $\diff R$ and $\diff R'$, then we obtain | |

$$ | |

n \to \{\, (x, y) \to (x', y') \mid n \ge 2 \wedge y' \le 1 - n + y \wedge x' \ge 1 + x \wedge y' \le x + y - x' \,\} | |

. | |

$$ | |

Note, however, that this is not the most accurate affine approximation that | |

can be obtained. That would be | |

$$ | |

n \to \{\, (x, y) \to (x', y') \mid y' \le 2 - n + x + y - x' \wedge n \ge 2 \wedge x' \ge 1 + x \,\} | |

. | |

$$ | |

\end{example} | |

\subsection{Checking Exactness} | |

The approximation $T$ for the transitive closure $R^+$ can be obtained | |

by projecting out the parameter $k$ from the approximation $K$ | |

\eqref{eq:transitive:approx} of the power $R^k$. | |

Since $K$ is an overapproximation of $R^k$, $T$ will also be an | |

overapproximation of $R^+$. | |

To check whether the results are exact, we need to consider two | |

cases depending on whether $R$ is {\em cyclic}, where $R$ is defined | |

to be cyclic if $R^+$ maps any element to itself, i.e., | |

$R^+ \cap \identity \ne \emptyset$. | |

If $R$ is acyclic, then the inductive definition of | |

\eqref{eq:transitive:inductive} is equivalent to its completion, | |

i.e., | |

$$ | |

R^+ = R \cup \left(R \circ R^+\right) | |

$$ | |

is a defining property. | |

Since $T$ is known to be an overapproximation, we only need to check | |

whether | |

$$ | |

T \subseteq R \cup \left(R \circ T\right) | |

. | |

$$ | |

This is essentially Theorem~5 of \shortciteN{Kelly1996closure}. | |

The only difference is that they only consider lexicographically | |

forward relations, a special case of acyclic relations. | |

If, on the other hand, $R$ is cyclic, then we have to resort | |

to checking whether the approximation $K$ of the power is exact. | |

Note that $T$ may be exact even if $K$ is not exact, so the check | |

is sound, but incomplete. | |

To check exactness of the power, we simply need to check | |

\eqref{eq:transitive:power}. Since again $K$ is known | |

to be an overapproximation, we only need to check whether | |

$$ | |

\begin{aligned} | |

K'|_{y_{d+1} - x_{d+1} = 1} & \subseteq R' | |

\\ | |

K'|_{y_{d+1} - x_{d+1} \ge 2} & \subseteq R' \circ K'|_{y_{d+1} - x_{d+1} \ge 1} | |

, | |

\end{aligned} | |

$$ | |

where $R' = \{\, \vec x' \to \vec y' \mid \vec x \to \vec y \in R | |

\wedge y_{d+1} - x_{d+1} = 1\,\}$, i.e., $R$ extended with path | |

lengths equal to 1. | |

All that remains is to explain how to check the cyclicity of $R$. | |

Note that the exactness on the power is always sound, even | |

in the acyclic case, so we only need to be careful that we find | |

all cyclic cases. Now, if $R$ is cyclic, i.e., | |

$R^+ \cap \identity \ne \emptyset$, then, since $T$ is | |

an overapproximation of $R^+$, also | |

$T \cap \identity \ne \emptyset$. This in turn means | |

that $\Delta \, K'$ contains a point whose first $d$ coordinates | |

are zero and whose final coordinate is positive. | |

In the implementation we currently perform this test on $P'$ instead of $K'$. | |

Note that if $R^+$ is acyclic and $T$ is not, then the approximation | |

is clearly not exact and the approximation of the power $K$ | |

will not be exact either. | |

\subsection{Decomposing $R$ into strongly connected components} | |

If the input relation $R$ is a union of several basic relations | |

that can be partially ordered | |

then the accuracy of the approximation may be improved by computing | |

an approximation of each strongly connected components separately. | |

For example, if $R = R_1 \cup R_2$ and $R_1 \circ R_2 = \emptyset$, | |

then we know that any path that passes through $R_2$ cannot later | |

pass through $R_1$, i.e., | |

\begin{equation} | |

\label{eq:transitive:components} | |

R^+ = R_1^+ \cup R_2^+ \cup \left(R_2^+ \circ R_1^+\right) | |

. | |

\end{equation} | |

We can therefore compute (approximations of) transitive closures | |

of $R_1$ and $R_2$ separately. | |

Note, however, that the condition $R_1 \circ R_2 = \emptyset$ | |

is actually too strong. | |

If $R_1 \circ R_2$ is a subset of $R_2 \circ R_1$ | |

then we can reorder the segments | |

in any path that moves through both $R_1$ and $R_2$ to | |

first move through $R_1$ and then through $R_2$. | |

This idea can be generalized to relations that are unions | |

of more than two basic relations by constructing the | |

strongly connected components in the graph with as vertices | |

the basic relations and an edge between two basic relations | |

$R_i$ and $R_j$ if $R_i$ needs to follow $R_j$ in some paths. | |

That is, there is an edge from $R_i$ to $R_j$ iff | |

\begin{equation} | |

\label{eq:transitive:edge} | |

R_i \circ R_j | |

\not\subseteq | |

R_j \circ R_i | |

. | |

\end{equation} | |

The components can be obtained from the graph by applying | |

Tarjan's algorithm \shortcite{Tarjan1972}. | |

In practice, we compute the (extended) powers $K_i'$ of each component | |

separately and then compose them as in \eqref{eq:transitive:decompose}. | |

Note, however, that in this case the order in which we apply them is | |

important and should correspond to a topological ordering of the | |

strongly connected components. Simply applying Tarjan's | |

algorithm will produce topologically sorted strongly connected components. | |

The graph on which Tarjan's algorithm is applied is constructed on-the-fly. | |

That is, whenever the algorithm checks if there is an edge between | |

two vertices, we evaluate \eqref{eq:transitive:edge}. | |

The exactness check is performed on each component separately. | |

If the approximation turns out to be inexact for any of the components, | |

then the entire result is marked inexact and the exactness check | |

is skipped on the components that still need to be handled. | |

It should be noted that \eqref{eq:transitive:components} | |

is only valid for exact transitive closures. | |

If overapproximations are computed in the right hand side, then the result will | |

still be an overapproximation of the left hand side, but this result | |

may not be transitively closed. If we only separate components based | |

on the condition $R_i \circ R_j = \emptyset$, then there is no problem, | |

as this condition will still hold on the computed approximations | |

of the transitive closures. If, however, we have exploited | |

\eqref{eq:transitive:edge} during the decomposition and if the | |

result turns out not to be exact, then we check whether | |

the result is transitively closed. If not, we recompute | |

the transitive closure, skipping the decomposition. | |

Note that testing for transitive closedness on the result may | |

be fairly expensive, so we may want to make this check | |

configurable. | |

\begin{figure} | |

\begin{center} | |

\begin{tikzpicture}[x=0.5cm,y=0.5cm,>=stealth,shorten >=1pt] | |

\foreach \x in {1,...,10}{ | |

\foreach \y in {1,...,10}{ | |

\draw[->] (\x,\y) -- (\x,\y+1); | |

} | |

} | |

\foreach \x in {1,...,20}{ | |

\foreach \y in {5,...,15}{ | |

\draw[->] (\x,\y) -- (\x+1,\y); | |

} | |

} | |

\end{tikzpicture} | |

\end{center} | |

\caption{The relation from \autoref{ex:closure4}} | |

\label{f:closure4} | |

\end{figure} | |

\begin{example} | |

\label{ex:closure4} | |

Consider the relation in example {\tt closure4} that comes with | |

the Omega calculator~\shortcite{Omega_calc}, $R = R_1 \cup R_2$, | |

with | |

$$ | |

\begin{aligned} | |

R_1 & = \{\, (x,y) \to (x,y+1) \mid 1 \le x,y \le 10 \,\} | |

\\ | |

R_2 & = \{\, (x,y) \to (x+1,y) \mid 1 \le x \le 20 \wedge 5 \le y \le 15 \,\} | |

. | |

\end{aligned} | |

$$ | |

This relation is shown graphically in \autoref{f:closure4}. | |

We have | |

$$ | |

\begin{aligned} | |

R_1 \circ R_2 &= | |

\{\, (x,y) \to (x+1,y+1) \mid 1 \le x \le 9 \wedge 5 \le y \le 10 \,\} | |

\\ | |

R_2 \circ R_1 &= | |

\{\, (x,y) \to (x+1,y+1) \mid 1 \le x \le 10 \wedge 4 \le y \le 10 \,\} | |

. | |

\end{aligned} | |

$$ | |

Clearly, $R_1 \circ R_2 \subseteq R_2 \circ R_1$ and so | |

$$ | |

\left( | |

R_1 \cup R_2 | |

\right)^+ | |

= | |

\left(R_2^+ \circ R_1^+\right) | |

\cup R_1^+ | |

\cup R_2^+ | |

. | |

$$ | |

\end{example} | |

\begin{figure} | |

\newcounter{n} | |

\newcounter{t1} | |

\newcounter{t2} | |

\newcounter{t3} | |

\newcounter{t4} | |

\begin{center} | |

\begin{tikzpicture}[>=stealth,shorten >=1pt] | |

\setcounter{n}{7} | |

\foreach \i in {1,...,\value{n}}{ | |

\foreach \j in {1,...,\value{n}}{ | |

\setcounter{t1}{2 * \j - 4 - \i + 1} | |

\setcounter{t2}{\value{n} - 3 - \i + 1} | |

\setcounter{t3}{2 * \i - 1 - \j + 1} | |

\setcounter{t4}{\value{n} - \j + 1} | |

\ifnum\value{t1}>0\ifnum\value{t2}>0 | |

\ifnum\value{t3}>0\ifnum\value{t4}>0 | |

\draw[thick,->] (\i,\j) to[out=20] (\i+3,\j); | |

\fi\fi\fi\fi | |

\setcounter{t1}{2 * \j - 1 - \i + 1} | |

\setcounter{t2}{\value{n} - \i + 1} | |

\setcounter{t3}{2 * \i - 4 - \j + 1} | |

\setcounter{t4}{\value{n} - 3 - \j + 1} | |

\ifnum\value{t1}>0\ifnum\value{t2}>0 | |

\ifnum\value{t3}>0\ifnum\value{t4}>0 | |

\draw[thick,->] (\i,\j) to[in=-20,out=20] (\i,\j+3); | |

\fi\fi\fi\fi | |

\setcounter{t1}{2 * \j - 1 - \i + 1} | |

\setcounter{t2}{\value{n} - 1 - \i + 1} | |

\setcounter{t3}{2 * \i - 1 - \j + 1} | |

\setcounter{t4}{\value{n} - 1 - \j + 1} | |

\ifnum\value{t1}>0\ifnum\value{t2}>0 | |

\ifnum\value{t3}>0\ifnum\value{t4}>0 | |

\draw[thick,->] (\i,\j) to (\i+1,\j+1); | |

\fi\fi\fi\fi | |

} | |

} | |

\end{tikzpicture} | |

\end{center} | |

\caption{The relation from \autoref{ex:decomposition}} | |

\label{f:decomposition} | |

\end{figure} | |

\begin{example} | |

\label{ex:decomposition} | |

Consider the relation on the right of \shortciteN[Figure~2]{Beletska2009}, | |

reproduced in \autoref{f:decomposition}. | |

The relation can be described as $R = R_1 \cup R_2 \cup R_3$, | |

with | |

$$ | |

\begin{aligned} | |

R_1 &= n \mapsto \{\, (i,j) \to (i+3,j) \mid | |

i \le 2 j - 4 \wedge | |

i \le n - 3 \wedge | |

j \le 2 i - 1 \wedge | |

j \le n \,\} | |

\\ | |

R_2 &= n \mapsto \{\, (i,j) \to (i,j+3) \mid | |

i \le 2 j - 1 \wedge | |

i \le n \wedge | |

j \le 2 i - 4 \wedge | |

j \le n - 3 \,\} | |

\\ | |

R_3 &= n \mapsto \{\, (i,j) \to (i+1,j+1) \mid | |

i \le 2 j - 1 \wedge | |

i \le n - 1 \wedge | |

j \le 2 i - 1 \wedge | |

j \le n - 1\,\} | |

. | |

\end{aligned} | |

$$ | |

The figure shows this relation for $n = 7$. | |

Both | |

$R_3 \circ R_1 \subseteq R_1 \circ R_3$ | |

and | |

$R_3 \circ R_2 \subseteq R_2 \circ R_3$, | |

which the reader can verify using the {\tt iscc} calculator: | |

\begin{verbatim} | |

R1 := [n] -> { [i,j] -> [i+3,j] : i <= 2 j - 4 and i <= n - 3 and | |

j <= 2 i - 1 and j <= n }; | |

R2 := [n] -> { [i,j] -> [i,j+3] : i <= 2 j - 1 and i <= n and | |

j <= 2 i - 4 and j <= n - 3 }; | |

R3 := [n] -> { [i,j] -> [i+1,j+1] : i <= 2 j - 1 and i <= n - 1 and | |

j <= 2 i - 1 and j <= n - 1 }; | |

(R1 . R3) - (R3 . R1); | |

(R2 . R3) - (R3 . R2); | |

\end{verbatim} | |

$R_3$ can therefore be moved forward in any path. | |

For the other two basic relations, we have both | |

$R_2 \circ R_1 \not\subseteq R_1 \circ R_2$ | |

and | |

$R_1 \circ R_2 \not\subseteq R_2 \circ R_1$ | |

and so $R_1$ and $R_2$ form a strongly connected component. | |

By computing the power of $R_3$ and $R_1 \cup R_2$ separately | |

and composing the results, the power of $R$ can be computed exactly | |

using \eqref{eq:transitive:singleton}. | |

As explained by \shortciteN{Beletska2009}, applying the same formula | |

to $R$ directly, without a decomposition, would result in | |

an overapproximation of the power. | |

\end{example} | |

\subsection{Partitioning the domains and ranges of $R$} | |

The algorithm of \autoref{s:power} assumes that the input relation $R$ | |

can be treated as a union of translations. | |

This is a reasonable assumption if $R$ maps elements of a given | |

abstract domain to the same domain. | |

However, if $R$ is a union of relations that map between different | |

domains, then this assumption no longer holds. | |

In particular, when an entire dependence graph is encoded | |

in a single relation, as is done by, e.g., | |

\shortciteN[Section~6.1]{Barthou2000MSE}, then it does not make | |

sense to look at differences between iterations of different domains. | |

Now, arguably, a modified Floyd-Warshall algorithm should | |

be applied to the dependence graph, as advocated by | |

\shortciteN{Kelly1996closure}, with the transitive closure operation | |

only being applied to relations from a given domain to itself. | |

However, it is also possible to detect disjoint domains and ranges | |

and to apply Floyd-Warshall internally. | |

\linesnumbered | |

\begin{algorithm} | |

\caption{The modified Floyd-Warshall algorithm of | |

\protect\shortciteN{Kelly1996closure}} | |

\label{a:Floyd} | |

\SetKwInput{Input}{Input} | |

\SetKwInput{Output}{Output} | |

\Input{Relations $R_{pq}$, $0 \le p, q < n$} | |

\Output{Updated relations $R_{pq}$ such that each relation | |

$R_{pq}$ contains all indirect paths from $p$ to $q$ in the input graph} | |

% | |

\BlankLine | |

\SetVline | |

\dontprintsemicolon | |

% | |

\For{$r \in [0, n-1]$}{ | |

$R_{rr} \coloneqq R_{rr}^+$ \nllabel{l:Floyd:closure}\; | |

\For{$p \in [0, n-1]$}{ | |

\For{$q \in [0, n-1]$}{ | |

\If{$p \ne r$ or $q \ne r$}{ | |

$R_{pq} \coloneqq R_{pq} \cup \left(R_{rq} \circ R_{pr}\right) | |

\cup \left(R_{rq} \circ R_{rr} \circ R_{pr}\right)$ | |

\nllabel{l:Floyd:update} | |

} | |

} | |

} | |

} | |

\end{algorithm} | |

Let the input relation $R$ be a union of $m$ basic relations $R_i$. | |

Let $D_{2i}$ be the domains of $R_i$ and $D_{2i+1}$ the ranges of $R_i$. | |

The first step is to group overlapping $D_j$ until a partition is | |

obtained. If the resulting partition consists of a single part, | |

then we continue with the algorithm of \autoref{s:power}. | |

Otherwise, we apply Floyd-Warshall on the graph with as vertices | |

the parts of the partition and as edges the $R_i$ attached to | |

the appropriate pairs of vertices. | |

In particular, let there be $n$ parts $P_k$ in the partition. | |

We construct $n^2$ relations | |

$$ | |

R_{pq} \coloneqq \bigcup_{i \text{ s.t. } \domain R_i \subseteq P_p \wedge | |

\range R_i \subseteq P_q} R_i | |

, | |

$$ | |

apply \autoref{a:Floyd} and return the union of all resulting | |

$R_{pq}$ as the transitive closure of $R$. | |

Each iteration of the $r$-loop in \autoref{a:Floyd} updates | |

all relations $R_{pq}$ to include paths that go from $p$ to $r$, | |

possibly stay there for a while, and then go from $r$ to $q$. | |

Note that paths that ``stay in $r$'' include all paths that | |

pass through earlier vertices since $R_{rr}$ itself has been updated | |

accordingly in previous iterations of the outer loop. | |

In principle, it would be sufficient to use the $R_{pr}$ | |

and $R_{rq}$ computed in the previous iteration of the | |

$r$-loop in Line~\ref{l:Floyd:update}. | |

However, from an implementation perspective, it is easier | |

to allow either or both of these to have been updated | |

in the same iteration of the $r$-loop. | |

This may result in duplicate paths, but these can usually | |

be removed by coalescing (\autoref{s:coalescing}) the result of the union | |

in Line~\ref{l:Floyd:update}, which should be done in any case. | |

The transitive closure in Line~\ref{l:Floyd:closure} | |

is performed using a recursive call. This recursive call | |

includes the partitioning step, but the resulting partition will | |

usually be a singleton. | |

The result of the recursive call will either be exact or an | |

overapproximation. The final result of Floyd-Warshall is therefore | |

also exact or an overapproximation. | |

\begin{figure} | |

\begin{center} | |

\begin{tikzpicture}[x=1cm,y=1cm,>=stealth,shorten >=3pt] | |

\foreach \x/\y in {0/0,1/1,3/2} { | |

\fill (\x,\y) circle (2pt); | |

} | |

\foreach \x/\y in {0/1,2/2,3/3} { | |

\draw (\x,\y) circle (2pt); | |

} | |

\draw[->] (0,0) -- (0,1); | |

\draw[->] (0,1) -- (1,1); | |

\draw[->] (2,2) -- (3,2); | |

\draw[->] (3,2) -- (3,3); | |

\draw[->,dashed] (2,2) -- (3,3); | |

\draw[->,dotted] (0,0) -- (1,1); | |

\end{tikzpicture} | |

\end{center} | |

\caption{The relation (solid arrows) on the right of Figure~1 of | |

\protect\shortciteN{Beletska2009} and its transitive closure} | |

\label{f:COCOA:1} | |

\end{figure} | |

\begin{example} | |

Consider the relation on the right of Figure~1 of | |

\shortciteN{Beletska2009}, | |

reproduced in \autoref{f:COCOA:1}. | |

This relation can be described as | |

$$ | |

\begin{aligned} | |

\{\, (x, y) \to (x_2, y_2) \mid {} & (3y = 2x \wedge x_2 = x \wedge 3y_2 = 3 + 2x \wedge x \ge 0 \wedge x \le 3) \vee {} \\ | |

& (x_2 = 1 + x \wedge y_2 = y \wedge x \ge 0 \wedge 3y \ge 2 + 2x \wedge x \le 2 \wedge 3y \le 3 + 2x) \,\} | |

. | |

\end{aligned} | |

$$ | |

Note that the domain of the upward relation overlaps with the range | |

of the rightward relation and vice versa, but that the domain | |

of neither relation overlaps with its own range or the domain of | |

the other relation. | |

The domains and ranges can therefore be partitioned into two parts, | |

$P_0$ and $P_1$, shown as the white and black dots in \autoref{f:COCOA:1}, | |

respectively. | |

Initially, we have | |

$$ | |

\begin{aligned} | |

R_{00} & = \emptyset | |

\\ | |

R_{01} & = | |

\{\, (x, y) \to (x+1, y) \mid | |

(x \ge 0 \wedge 3y \ge 2 + 2x \wedge x \le 2 \wedge 3y \le 3 + 2x) \,\} | |

\\ | |

R_{10} & = | |

\{\, (x, y) \to (x_2, y_2) \mid (3y = 2x \wedge x_2 = x \wedge 3y_2 = 3 + 2x \wedge x \ge 0 \wedge x \le 3) \,\} | |

\\ | |

R_{11} & = \emptyset | |

. | |

\end{aligned} | |

$$ | |

In the first iteration, $R_{00}$ remains the same ($\emptyset^+ = \emptyset$). | |

$R_{01}$ and $R_{10}$ are therefore also unaffected, but | |

$R_{11}$ is updated to include $R_{01} \circ R_{10}$, i.e., | |

the dashed arrow in the figure. | |

This new $R_{11}$ is obviously transitively closed, so it is not | |

changed in the second iteration and it does not have an effect | |

on $R_{01}$ and $R_{10}$. However, $R_{00}$ is updated to | |

include $R_{10} \circ R_{01}$, i.e., the dotted arrow in the figure. | |

The transitive closure of the original relation is then equal to | |

$R_{00} \cup R_{01} \cup R_{10} \cup R_{11}$. | |

\end{example} | |

\subsection{Incremental Computation} | |

\label{s:incremental} | |

In some cases it is possible and useful to compute the transitive closure | |

of union of basic relations incrementally. In particular, | |

if $R$ is a union of $m$ basic maps, | |

$$ | |

R = \bigcup_j R_j | |

, | |

$$ | |

then we can pick some $R_i$ and compute the transitive closure of $R$ as | |

\begin{equation} | |

\label{eq:transitive:incremental} | |

R^+ = R_i^+ \cup | |

\left( | |

\bigcup_{j \ne i} | |

R_i^* \circ R_j \circ R_i^* | |

\right)^+ | |

. | |

\end{equation} | |

For this approach to be successful, it is crucial that each | |

of the disjuncts in the argument of the second transitive | |

closure in \eqref{eq:transitive:incremental} be representable | |

as a single basic relation, i.e., without a union. | |

If this condition holds, then by using \eqref{eq:transitive:incremental}, | |

the number of disjuncts in the argument of the transitive closure | |

can be reduced by one. | |

Now, $R_i^* = R_i^+ \cup \identity$, but in some cases it is possible | |

to relax the constraints of $R_i^+$ to include part of the identity relation, | |

say on domain $D$. We will use the notation | |

${\cal C}(R_i,D) = R_i^+ \cup \identity_D$ to represent | |

this relaxed version of $R^+$. | |

\shortciteN{Kelly1996closure} use the notation $R_i^?$. | |

${\cal C}(R_i,D)$ can be computed by allowing $k$ to attain | |

the value $0$ in \eqref{eq:transitive:Q} and by using | |

$$ | |

P \cap \left(D \to D\right) | |

$$ | |

instead of \eqref{eq:transitive:approx}. | |

Typically, $D$ will be a strict superset of both $\domain R_i$ | |

and $\range R_i$. We therefore need to check that domain | |

and range of the transitive closure are part of ${\cal C}(R_i,D)$, | |

i.e., the part that results from the paths of positive length ($k \ge 1$), | |

are equal to the domain and range of $R_i$. | |

If not, then the incremental approach cannot be applied for | |

the given choice of $R_i$ and $D$. | |

In order to be able to replace $R^*$ by ${\cal C}(R_i,D)$ | |

in \eqref{eq:transitive:incremental}, $D$ should be chosen | |

to include both $\domain R$ and $\range R$, i.e., such | |

that $\identity_D \circ R_j \circ \identity_D = R_j$ for all $j\ne i$. | |

\shortciteN{Kelly1996closure} say that they use | |

$D = \domain R_i \cup \range R_i$, but presumably they mean that | |

they use $D = \domain R \cup \range R$. | |

Now, this expression of $D$ contains a union, so it not directly usable. | |

\shortciteN{Kelly1996closure} do not explain how they avoid this union. | |

Apparently, in their implementation, | |

they are using the convex hull of $\domain R \cup \range R$ | |

or at least an approximation of this convex hull. | |

We use the simple hull (\autoref{s:simple hull}) of $\domain R \cup \range R$. | |

It is also possible to use a domain $D$ that does {\em not\/} | |

include $\domain R \cup \range R$, but then we have to | |

compose with ${\cal C}(R_i,D)$ more selectively. | |

In particular, if we have | |

\begin{equation} | |

\label{eq:transitive:right} | |

\text{for each $j \ne i$ either } | |

\domain R_j \subseteq D \text{ or } \domain R_j \cap \range R_i = \emptyset | |

\end{equation} | |

and, similarly, | |

\begin{equation} | |

\label{eq:transitive:left} | |

\text{for each $j \ne i$ either } | |

\range R_j \subseteq D \text{ or } \range R_j \cap \domain R_i = \emptyset | |

\end{equation} | |

then we can refine \eqref{eq:transitive:incremental} to | |

$$ | |

R_i^+ \cup | |

\left( | |

\left( | |

\bigcup_{\shortstack{$\scriptstyle\domain R_j \subseteq D $\\ | |

$\scriptstyle\range R_j \subseteq D$}} | |

{\cal C} \circ R_j \circ {\cal C} | |

\right) | |

\cup | |

\left( | |

\bigcup_{\shortstack{$\scriptstyle\domain R_j \cap \range R_i = \emptyset$\\ | |

$\scriptstyle\range R_j \subseteq D$}} | |

\!\!\!\!\! | |

{\cal C} \circ R_j | |

\right) | |

\cup | |

\left( | |

\bigcup_{\shortstack{$\scriptstyle\domain R_j \subseteq D $\\ | |

$\scriptstyle\range R_j \cap \domain R_i = \emptyset$}} | |

\!\!\!\!\! | |

R_j \circ {\cal C} | |

\right) | |

\cup | |

\left( | |

\bigcup_{\shortstack{$\scriptstyle\domain R_j \cap \range R_i = \emptyset$\\ | |

$\scriptstyle\range R_j \cap \domain R_i = \emptyset$}} | |

\!\!\!\!\! | |

R_j | |

\right) | |

\right)^+ | |

. | |

$$ | |

If only property~\eqref{eq:transitive:right} holds, | |

we can use | |

$$ | |

R_i^+ \cup | |

\left( | |

\left( | |

R_i^+ \cup \identity | |

\right) | |

\circ | |

\left( | |

\left( | |

\bigcup_{\shortstack{$\scriptstyle\domain R_j \subseteq D $}} | |

R_j \circ {\cal C} | |

\right) | |

\cup | |

\left( | |

\bigcup_{\shortstack{$\scriptstyle\domain R_j \cap \range R_i = \emptyset$}} | |

\!\!\!\!\! | |

R_j | |

\right) | |

\right)^+ | |

\right) | |

, | |

$$ | |

while if only property~\eqref{eq:transitive:left} holds, | |

we can use | |

$$ | |

R_i^+ \cup | |

\left( | |

\left( | |

\left( | |

\bigcup_{\shortstack{$\scriptstyle\range R_j \subseteq D $}} | |

{\cal C} \circ R_j | |

\right) | |

\cup | |

\left( | |

\bigcup_{\shortstack{$\scriptstyle\range R_j \cap \domain R_i = \emptyset$}} | |

\!\!\!\!\! | |

R_j | |

\right) | |

\right)^+ | |

\circ | |

\left( | |

R_i^+ \cup \identity | |

\right) | |

\right) | |

. | |

$$ | |

It should be noted that if we want the result of the incremental | |

approach to be transitively closed, then we can only apply it | |

if all of the transitive closure operations involved are exact. | |

If, say, the second transitive closure in \eqref{eq:transitive:incremental} | |

contains extra elements, then the result does not necessarily contain | |

the composition of these extra elements with powers of $R_i$. | |

\subsection{An {\tt Omega}-like implementation} | |

While the main algorithm of \shortciteN{Kelly1996closure} is | |

designed to compute and underapproximation of the transitive closure, | |

the authors mention that they could also compute overapproximations. | |

In this section, we describe our implementation of an algorithm | |

that is based on their ideas. | |

Note that the {\tt Omega} library computes underapproximations | |

\shortcite[Section 6.4]{Omega_lib}. | |

The main tool is Equation~(2) of \shortciteN{Kelly1996closure}. | |

The input relation $R$ is first overapproximated by a ``d-form'' relation | |

$$ | |

\{\, \vec i \to \vec j \mid \exists \vec \alpha : | |

\vec L \le \vec j - \vec i \le \vec U | |

\wedge | |

(\forall p : j_p - i_p = M_p \alpha_p) | |

\,\} | |

, | |

$$ | |

where $p$ ranges over the dimensions and $\vec L$, $\vec U$ and | |

$\vec M$ are constant integer vectors. The elements of $\vec U$ | |

may be $\infty$, meaning that there is no upper bound corresponding | |

to that element, and similarly for $\vec L$. | |

Such an overapproximation can be obtained by computing strides, | |

lower and upper bounds on the difference set $\Delta \, R$. | |

The transitive closure of such a ``d-form'' relation is | |

\begin{equation} | |

\label{eq:omega} | |

\{\, \vec i \to \vec j \mid \exists \vec \alpha, k : | |

k \ge 1 \wedge | |

k \, \vec L \le \vec j - \vec i \le k \, \vec U | |

\wedge | |

(\forall p : j_p - i_p = M_p \alpha_p) | |

\,\} | |

. | |

\end{equation} | |

The domain and range of this transitive closure are then | |

intersected with those of the input relation. | |

This is a special case of the algorithm in \autoref{s:power}. | |

In their algorithm for computing lower bounds, the authors | |

use the above algorithm as a substep on the disjuncts in the relation. | |

At the end, they say | |

\begin{quote} | |

If an upper bound is required, it can be calculated in a manner | |

similar to that of a single conjunct [sic] relation. | |

\end{quote} | |

Presumably, the authors mean that a ``d-form'' approximation | |

of the whole input relation should be used. | |

However, the accuracy can be improved by also trying to | |

apply the incremental technique from the same paper, | |

which is explained in more detail in \autoref{s:incremental}. | |

In this case, ${\cal C}(R_i,D)$ can be obtained by | |

allowing the value zero for $k$ in \eqref{eq:omega}, | |

i.e., by computing | |

$$ | |

\{\, \vec i \to \vec j \mid \exists \vec \alpha, k : | |

k \ge 0 \wedge | |

k \, \vec L \le \vec j - \vec i \le k \, \vec U | |

\wedge | |

(\forall p : j_p - i_p = M_p \alpha_p) | |

\,\} | |

. | |

$$ | |

In our implementation we take as $D$ the simple hull | |

(\autoref{s:simple hull}) of $\domain R \cup \range R$. | |

To determine whether it is safe to use ${\cal C}(R_i,D)$, | |

we check the following conditions, as proposed by | |

\shortciteN{Kelly1996closure}: | |

${\cal C}(R_i,D) - R_i^+$ is not a union and for each $j \ne i$ | |

the condition | |

$$ | |

\left({\cal C}(R_i,D) - R_i^+\right) | |

\circ | |

R_j | |

\circ | |

\left({\cal C}(R_i,D) - R_i^+\right) | |

= | |

R_j | |

$$ | |

holds. |