X-Git-Url: http://review.tizen.org/git/?a=blobdiff_plain;f=doc%2Fimplementation.tex;h=d5ece80e4c8ba7876ed2769f38434b4a72848dfe;hb=63fb8a7f484648c3caa25351c8c94ac2395ec563;hp=cec087439ebb41d18c9f8a82384a2cf0c9e7fa31;hpb=94b16637d9adbc3bdb6466b675a0d1b328354899;p=platform%2Fupstream%2Fisl.git diff --git a/doc/implementation.tex b/doc/implementation.tex index cec0874..d5ece80 100644 --- a/doc/implementation.tex +++ b/doc/implementation.tex @@ -94,7 +94,7 @@ The difference set ($\Delta \, R$) of $R$ is the set of differences between image elements and the corresponding domain elements, $$ -\Delta \, R \coloneqq +\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 @@ -151,6 +151,669 @@ 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. @@ -288,7 +951,7 @@ P = \{\, \vec x \to \vec y \mid \end{equation} and then the approximation computed in \eqref{eq:transitive:approx} is essentially the same as that of \shortciteN{Beletska2009}. -If some of $\Delta_i$s are not singleton sets or if +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. @@ -445,6 +1108,113 @@ A_1 & \vec c_1 \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 @@ -505,6 +1275,76 @@ n \to \{\, x \to y \mid \exists \, \alpha_0, \alpha_1: 7\alpha_1 = -2 + n \wedge $$ \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 @@ -988,7 +1828,7 @@ $$ 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 part of ${\cal C}(R_i,D)$, +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 @@ -1017,7 +1857,7 @@ In particular, if we have \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, either +and, similarly, \begin{equation} \label{eq:transitive:left} \text{for each $j \ne i$ either }