doc/implementation.tex

   1 \section{Implementation details}
   2
   3 \subsection{An interior point of a polyhedron}
   4 \label{s:interior}
   5
   6 We often need a point that lies in the interior of a polyhedron.
   7 The function \ai[\tt]{inner\_point} implements the following algorithm.
   8 Each polyhedron $P$ can be written as the sum of a polytope $P'$ and a cone $C$
   9 (the \ai{recession cone} or \ai{characteristic cone} of $P$).
  10 Adding a positive multiple of the sum of the extremal rays of $C$ to
  11 the \ai{barycenter}
  12 $$
  13 \frac 1 N \sum_i \vec v_i(\vec p)
  14 $$
  15 of $P'$, where $N$ is the number of vertices, results in a point
  16 in the interior of $P$.
  17
  18 \subsection{The integer points in the fundamental parallelepiped of a simple cone}
  19
  20 \label{s:fundamental}
  21
  22 This section is based on \shortciteN[Lemma 5.1]{Barvinok1992volume} and
  23 \shortciteN{Koeppe2006experiments}.
  24
  25 \sindex{simple}{cone}
  26 \sindex{open}{facet}
  27 \sindex{open}{ray}
  28 \sindex{explicit}{representation}
  29 In this section we will deal exclusively with \ai{simple cone}s,
  30 i.e. $d$-dimensional cones with $d$ extremal rays and $d$ facets.
  31 \index{open facet}%
  32 Some of the facets of these cones may be open.
  33 Since we will mostly be dealing with cones in their
  34 \ai{explicit representation}, we will have occasion to speak of
  35 ``\ai{open ray}s'', by which we will mean that the facet not
  36 containing the ray is open.  (There is only one such facet because the cone
  37 is simple.)
  38
  39 \sindex{fundamental}{parallelepiped}
  40 \begin{definition}[Fundamental parallelepiped]
  41 Let $K = \vec v + \poshull \lb\, \vec u_i \,\rb$ be
  42 a closed (shifted) cone, then the \defindex{fundamental parallelepiped} $\Pi$
  43 of $K$ is
  44 $$
  45 \Pi = \vec v +
  46 \lb\, \sum_i \alpha_i \vec u_i \mid 0 \leq \alpha_i < 1 \,\rb
  47 .
  48 $$
  49 If some of the rays $\vec u_i$ of $K$ are open, then the constraints on
  50 the corresponding coefficient $\alpha_i$ are such that $0 < \alpha_i \le 1$.
  51 \end{definition}
  52
  53 \begin{lemma}[Integer points in the fundamental parallelepiped of a simple cone]
  54 Let $K = \vec v + \poshull \lb\, \vec u_i \,\rb$ be a closed simple cone
  55 and let $A$ be the matrix with the generators $\vec u_i$ of $K$
  56 as rows.
  57 Furthermore let $V A W^{-1} = S = \diag \vec s$ be the \indac{SNF} of $A$.
  58 Then the integer points in the fundamental parallelepiped of $K$ are given
  59 by
  60 \begin{eqnarray}
  61 \label{eq:parallelepiped}
  62 \vec w^\T & = & \vec v^\T + \fractional{(\vec k^\T W - \vec v^\T) A^{-1}} A
  63 \\
  64 \nonumber
  65 & = &
  66 \vec v^T +
  67 \sum_{i=1}^d
  68     \fractional{\sps{\sum_{j=1}^d k_j \vec w^T_j - \vec v^\T}{\vec u^*_i}} \vec u_i,
  69 \end{eqnarray}
  70 where $\vec u^*_i$ are the columns of $A^{-1}$ and $k_j \in \ZZ$ ranges
  71 over $0 \le k_j < s_j$.
  72 \end{lemma}
  73
  74 \begin{proof}
  75 Since $0 \le \fractional{x} < 1$, it is clear that each such $\vec w$
  76 lies inside the fundamental parallelepiped.
  77 Furthermore,
  78 \begin{eqnarray*}
  79 \vec w^\T & = & \vec v^\T + \fractional{(\vec k^\T W - \vec v^\T) A^{-1}} A
  80 \\
  81 & = &
  82 \vec v^T +
  83 \left(
  84 (\vec k^\T W - \vec v^\T) A^{-1} - \floor{(\vec k^\T W - \vec v^\T) A^{-1}}
  85 \right) A
  86 \\
  87 & = &
  88 \underbrace{\vec k^\T W\mathstrut}_{\in \ZZ^{1\times d}}
  89 +
  90 \underbrace{\floor{(\vec k^\T W - \vec v^\T) A^{-1}}}_{\in \ZZ^{1\times d}}
  91 \underbrace{A\mathstrut}_{\in \ZZ^{d\times d}} \in \ZZ^{1\times d}.
  92 \end{eqnarray*}
  93 Finally, if two such $\vec w$ are equal, i.e., $\vec w_1 = \vec w_2$,
  94 then
  95 \begin{eqnarray*}
  96 \vec 0^\T = \vec w_1^\T - \vec w_2^\T
  97 & = & \vec k_1^\T W - \vec k_2^\T W + \vec p^\T A
  98 \\
  99 & = & \left(\vec k_1^\T - \vec k_2^\T \right) W + \vec p^\T V^{-1} S W,
 100 \end{eqnarray*}
 101 with $\vec p \in \ZZ^d$,
 102 or $\vec k_1 \equiv \vec k_2 \mod \vec s$, i.e., $\vec k_1 = \vec k_2$.
 103 Since $\det S = \det A$, we obtain all points in the fundamental parallelepiped
 104 by taking all $\vec k \in \ZZ^d$ satisfying $0 \le k_j < s_j$.
 105 \end{proof}
 106
 107 If the cone $K$ is not closed then the coefficients of the open rays
 108 should be in $(0,1]$ rather than in $[0,1)$.
 109 In (\ref{eq:parallelepiped}),
 110 we therefore need to replace the fractional part $\fractional{x} = x - \floor{x}$
 111 by $\cractional{x} = x - \ceil{x-1}$ for the open rays.
 112
 113 \begin{figure}
 114 \intercol=1.2cm
 115 \begin{xy}
 116 <\intercol,0pt>:<0pt,\intercol>::
 117 \POS@i@={(0,-3),(0,0),(4,2),(4,-3)},{0*[grey]\xypolyline{*}}
 118 \POS@i@={(0,-3),(0,0),(4,2)},{0*[|(2)]\xypolyline{}}
 119 \POS(-1,0)\ar(4.5,0)
 120 \POS(0,-3)\ar(0,3)
 121 \POS(0,0)\ar@[|(3)](0,-1)
 122 \POS(0,0)\ar@[|(3)](2,1)
 123 \POS(0,-1)\ar@{--}@[|(2)](2,0)
 124 \POS(2,1)\ar@{--}@[|(2)](2,0)
 125 \POS(0,0)*{\bullet}
 126 \POS(1,0)*{\bullet}
 127 \end{xy}
 128 \caption{The integer points in the fundamental parallelepiped of $K$}
 129 \label{f:parallelepiped}
 130 \end{figure}
 131
 132 \begin{example}
 133 Let $K$ be the cone
 134 $$
 135 K = \sm{0 \\ 0} + \poshull \lb\, \sm{2 \\ 1}, \sm{0 \\ -1} \,\rb
 136 ,
 137 $$
 138 shown in Figure~\ref{f:parallelepiped}.
 139 Then
 140 $$
 141 A = \sm{2 & 1\\0 & -1} \qquad A^{-1} = \sm{1/2 & 1/2 \\ 0 & -1 }
 142 $$
 143 and
 144 $$
 145 \sm{1 & 0 \\ 1 & 1 } \sm{2 & 1\\0 & -1} = \sm{1 & 0 \\ 0 & 2} \sm{2 & 1 \\ 1 & 0}.
 146 $$
 147 We have $\det A = \det S = 2$ and
 148 $\vec k_1^\T = \sm{0 & 0}$ and $\vec k_2^\T = \sm{0 & 1}$.
 149 Therefore,
 150 $$
 151 \vec w_1^\T = \fractional{\sm{0 & 0} \sm{2 & 1 \\ 1 & 0} \sm{1/2 & 1/2 \\ 0 & -1 }}
 152 \sm{2 & 1\\0 & -1} = \sm{0 & 0}
 153 $$
 154 and
 155 \begin{eqnarray*}
 156 \vec w_2^\T & = &
 157 \fractional{\sm{0 & 1} \sm{2 & 1 \\ 1 & 0} \sm{1/2 & 1/2 \\ 0 & -1 }}
 158 \sm{2 & 1\\0 & -1}
 159 \\
 160 & = &
 161 \sm{1/2 & 1/2} \sm{2 & 1\\0 & -1} = \sm{1 & 0}.
 162 \end{eqnarray*}
 163 \end{example}
 164
 165
 166
 167
 168 \subsection{Barvinok's decomposition of simple cones in primal space}
 169 \label{s:decomposition}
 170
 171 As described by \shortciteN{DeLoera2003effective}, the first
 172 implementation of Barvinok's counting algorithm applied
 173 \ai{Barvinok's decomposition} \shortcite{Barvinok1994} in the \ai{dual space}.
 174 \ai{Brion's polarization trick} \shortcite{Brion88} then ensures that you
 175 do not need to worry about lower-dimensional faces in the decomposition.
 176 Another way of avoiding the lower-dimensional faces, in the \ai{primal space},
 177 is to perturb the vertex of the cone such that none of the lower-dimensional
 178 face encountered contain any integer points \shortcite{Koeppe2006primal}.
 179 In this section, we describe another technique that is based on allowing
 180 some of the facets of the cone to be open.
 181
 182 The basic step in Barvinok's decomposition is to replace a
 183 $d$-dimensional simple cone
 184 $K = \poshull \lb\, \vec u_i \,\rb_{i=1}^d \subset \QQ^d$
 185 by a signed sum of (at most) $d$ cones $K_j$
 186 with a smaller determinant (in absolute value).
 187 The cones are obtained by successively replacing each generator
 188 of $K$ by an appropriately chosen
 189 $\vec w = \sum_{i=1}^d \alpha_i \vec u_i$, i.e.,
 190 \begin{equation}
 191 \label{eq:K_j}
 192 K_j =
 193 \poshull \left(\lb\, \vec u_i \,\rb_{i=1}^d
 194 \setminus \{\, \vec u_j \,\} \cup \{\, \vec w \,\}\right)
 195 .
 196 \end{equation}
 197 To see that we can use these $K_j$ to perform a decomposition,
 198 rearrange the $\vec u_i$ such that for all $1 \le i \le k$ we have
 199 $\alpha_i < 0$ and for all $k+1 \le i \le d'$ we have $\alpha_i > 0$,
 200 with $d - d'$ the number of zero $\alpha_i$.
 201 We may assume $k < d'$; otherwise replace $\vec w \in B$ by
 202 $-\vec w \in B$.  We have
 203 $$
 204 \vec w + \sum_{i=1}^k (-\alpha_i) \vec u_i =
 205 \sum_{i=k+1}^{d'} \alpha_i \vec u_i
 206 $$
 207 or
 208 \begin{equation}
 209 \label{eq:sub}
 210 \sum_{i=0}^k \beta_i \vec u_i =
 211 \sum_{i=k+1}^{d'} \alpha_i \vec u_i
 212 ,
 213 \end{equation}
 214 with $\vec u_0 = \vec w$, $\beta_0 = 1$ and $\beta_i = -\alpha_i > 0$
 215 for $1 \le i \le k$.  Any two $\vec u_j$ and $\vec u_l$ on the same side
 216 of the equality are on opposite sides of the linear hull $H$ of
 217 the other $\vec u_i$s since there exists a convex combination
 218 of $\vec u_j$ and $\vec u_l$ on this hyperplane.
 219 In particular, since $\alpha_j$ and $\alpha_l$ have the same sign,
 220 we have
 221 \begin{equation}
 222 \label{eq:opposite}
 223 \frac {\alpha_j}{\alpha_j+\alpha_l} \vec u_j
 224 +
 225 \frac {\alpha_l}{\alpha_j+\alpha_l} \vec u_l
 226 \in H
 227 \qquad\text{for $\alpha_i \alpha_l > 0$}
 228 .
 229 \end{equation}
 230 The corresponding cones $K_j$ and $K_l$ (with $K_0 = K$)
 231 therefore intersect in a common face $F \subset H$.
 232 Let
 233 $$
 234 K' :=
 235 \poshull \left(\lb\, \vec u_i \,\rb_{i=1}^d \cup \{\, \vec w \,\}\right)
 236 ,
 237 $$
 238 then any $\vec x \in K'$ lies both in some cone $K_i$ with
 239 $0 \le i \le k$ and in some cone $K_i$ with $k+1 \le i \le d'$.
 240 (Just subtract an appropriate multiple of Equation~(\ref{eq:sub}).)
 241 The cones
 242 $\{\, K_i \,\}_{i=0}^k$
 243 and
 244 $\{\, K_i \,\}_{i=k+1}^{d'}$
 245 therefore both form a triangulation of $K'$ and hence
 246 \begin{equation}
 247 \label{eq:triangulations}
 248 \indf{K'}
 249 =
 250 \indf{K} + \sum_{i=1}^k \indf{K_i} - \sum_{j\in J_1} \indf{F_j}
 251 =
 252 \sum_{i=k+1}^{d'} \indf{K_i} - \sum_{j\in J_2} \indf{F_j}
 253 \end{equation}
 254 or
 255 \begin{equation}
 256 \label{eq:decomposition}
 257 \indf{K} = \sum_{i=1}^{d'} \varepsilon_i \indf{K_i} + \sum_j \delta_j \indf{F_j}
 258 ,
 259 \end{equation}
 260 with $\varepsilon_i = -1$ for $1 \le i \le k$,
 261 $\varepsilon_i = 1$ for $k+1 \le i \le d'$,
 262 $\delta_j \in \{ -1, 1 \}$ and $F_j$ some lower-dimensional faces.
 263 Figure~\ref{fig:w} shows the possible configurations
 264 in the case of a $3$-dimensional cone.
 265
 266 \begin{figure}
 267 \intercol=0.48cm
 268 \begin{center}
 269 \begin{minipage}{0cm}
 270 \begin{xy}
 271 <\intercol,0pt>:<0pt,\intercol>::
 272 *
 273 \xybox{
 274 \POS(-2,-1)="a"*+!U{+}
 275 \POS(2,0)="b"*+!L{+}
 276 \POS(0,2)="c"*+!D{+}
 277 \POS(0,0)="w"*+!DR{\vec w}
 278 \POS"a"\ar@{-}"b"
 279 \POS"b"\ar@{-}"c"
 280 \POS"c"\ar@{-}"a"
 281 \POS"a"\ar@{--}"w"
 282 \POS"b"\ar@{--}"w"
 283 \POS"c"\ar@{--}"w"
 284 }="a"
 285 +R+(2,0)*!L
 286 \xybox{
 287 \POS(-2,-1)="a"*+!U{+}
 288 \POS(2,0)="b"*+!L{-}
 289 \POS(0,2)="c"*+!D{+}
 290 \POS(-3,1)="w"*+!DR{\vec w}
 291 \POS"a"\ar@{-}"b"
 292 \POS"b"\ar@{-}"c"
 293 \POS"c"\ar@{-}"a"
 294 \POS"a"\ar@{--}"w"
 295 \POS"b"\ar@{--}"w"
 296 \POS"c"\ar@{--}"w"
 297 }="b"
 298 +R+(2,0)*!L
 299 \xybox{
 300 \POS(-2,-1)="a"*+!U{-}
 301 \POS(2,0)="b"*+!U{+}
 302 \POS(0,2)="c"*+!D{-}
 303 \POS(5,-1)="w"*+!L{\vec w}
 304 \POS"a"\ar@{-}"b"
 305 \POS"b"\ar@{-}"c"
 306 \POS"c"\ar@{-}"a"
 307 \POS"a"\ar@{--}"w"
 308 \POS"b"\ar@{--}"w"
 309 \POS"c"\ar@{--}"w"
 310 }
 311 \POS"a"
 312 +D-(0,1)*!U
 313 \xybox{
 314 \POS(-2,-1)="a"*+!U{0}
 315 \POS(2,0)="b"*+!L{+}
 316 \POS(0,2)="c"*+!D{+}
 317 \POS(1,1)="w"*+!DL{\vec w}
 318 \POS"a"\ar@{-}"b"
 319 \POS"b"\ar@{-}"c"
 320 \POS"c"\ar@{-}"a"
 321 \POS"a"\ar@{--}"w"
 322 }
 323 \POS"b"
 324 +DL-(0,1)*!UL
 325 \xybox{
 326 \POS(-2,-1)="a"*+!U{0}
 327 \POS(2,0)="b"*+!U{+}
 328 \POS(0,2)="c"*+!D{-}
 329 \POS(4,-2)="w"*+!L{\vec w}
 330 \POS"a"\ar@{-}"b"
 331 \POS"b"\ar@{-}"c"
 332 \POS"c"\ar@{-}"a"
 333 \POS"a"\ar@{--}"w"
 334 \POS"b"\ar@{--}"w"
 335 }
 336 \end{xy}
 337 \end{minipage}
 338 \end{center}
 339 \caption[Possible locations of the vector $\vec w$ with respect to the rays
 340 of a $3$-dimensional cone.]
 341 {Possible locations of $\vec w$ with respect to the rays
 342 of a $3$-dimensional cone.  The figure shows a section of the cones.}
 343 \label{fig:w}
 344 \end{figure}
 345
 346 As explained above there are several ways of avoiding the lower-dimensional
 347 faces in (\ref{eq:decomposition}).  Here we will apply the following proposition.
 348 \begin{proposition}[\shortciteN{Koeppe2007parametric}]
 349 \label{p:inclusion-exclusion}
 350   Let
 351   \begin{equation}
 352     \label{eq:full-source-identity}
 353     \sum_{i\in {I_1}} \epsilon_i [P_i] + \sum_{i\in {I_2}} \delta_k [P_i] = 0
 354   \end{equation}
 355   be a (finite) linear identity of indicator functions of closed
 356   polyhedra~$P_i\subseteq\QQ^d$, where the
 357   polyhedra~$P_i$ with $i \in I_1$ are full-dimensional and those with $i \in I_2$
 358   lower-dimensional.  Let each closed polyhedron be given as
 359 $$
 360     P_i = \left\{\, \vec x \mid \sp{b^*_{i,j}}{x} \ge \beta_{i,j} \text{
 361       for $j\in J_i$}\,\right\}
 362   .
 363 $$
 364   Let $\vec y\in\QQ^d$ be a vector such that $\langle \vec b^*_{i,j}, \vec
 365   y\rangle \neq 0$ for all $i\in I_1\cup I_2$, $j\in J_i$.
 366   For each $i\in I_1$, we define the half-open polyhedron
 367   \begin{equation}
 368     \label{eq:half-open-by-y}
 369     \begin{aligned}
 370       \tilde P_i = \Bigl\{\, \vec x\in\QQ^d \mid {}&
 371             \sp{b^*_{i,j}}{x} \ge \beta_{i,j}
 372       \text{ for $j\in J_i$ with $\sp{b^*_{i,j}}{y} > 0$,} \\
 373       & \sp{b^*_{i,j}}{x} > \beta_{i,j}
 374       \text{ for $j\in J_i$ with $\sp{b^*_{i,j}}{y} < 0$} \,\Bigr\}.
 375     \end{aligned}
 376   \end{equation}
 377   Then
 378   \begin{equation}
 379     \label{eq:target-identity}
 380     \sum_{i\in I_1} \epsilon_i [\tilde P_i] = 0.
 381   \end{equation}
 382 \end{proposition}
 383 When applying this proposition to (\ref{eq:decomposition}), we obtain
 384 \begin{equation}
 385 \label{eq:decomposition:2}
 386 \indf{\tilde K} = \sum_{i=1}^{d'} \varepsilon_i \indf{\tilde K_i}
 387 ,
 388 \end{equation}
 389 where we start out
 390 from a given $\tilde K$, which may be $K$ itself, i.e., a fully closed cone,
 391 or the result of a previous application of the proposition, either through
 392 a triangulation (Section~\ref{s:triangulation}) or a previous decomposition.
 393 In either case, a suitable $\vec y$ is available, either as an interior
 394 point of the cone or as the vector used in the previous application
 395 (which may require a slight perturbation if it happens to lie on one of
 396 the new facets of the cones $K_i$).
 397 We are, however, free to construct a new $\vec y$ on each application
 398 of the proposition.
 399 In fact, we will not even construct such a vector explicitly, but
 400 rather apply a set of rules that is equivalent to a valid choice of $\vec y$.
 401 Below, we will present an ``intuitive'' motivation for these rules.
 402 For a more algebraic, shorter, and arguably simpler motivation we
 403 refer to \shortciteN{Koeppe2007parametric}.
 404
 405 The vector $\vec y$ has to satisfy $\sp{b^*_j}y > 0$ for normals $\vec b^*_j$
 406 of closed facets and $\sp{b^*_j}y < 0$ for normals $\vec b^*_j$ of open facets of
 407 $\tilde K$.
 408 These constraints delineate a non-empty open cone $R$ from which
 409 $\vec y$ should be selected.  For some of the new facets of the cones
 410 $\tilde K_j$, the cone $R$ will not be cut by the affine hull of the facet.
 411 The closedness of these facets is therefore predetermined by $\tilde K$.
 412 For the other facets, a choice will have to be made.
 413 To be able to make the choice based on local information and without
 414 computing an explicit vector $\vec y$, we use the following convention.
 415 We first assign an arbitrary total order to the rays.
 416 If (the affine hull of) a facet separates the two rays not on the facet $\vec u_i$
 417 and $\vec u_j$, i.e., $\alpha_i \alpha_j > 0$ (\ref{eq:opposite}), then
 418 we choose $\vec y$ to lie on the side of the smallest ray, according
 419 to the chosen order.
 420 That is, $\sp{{\tilde n}_{ij}}y > 0$, for
 421 $\vec {\tilde n}_{ij}$ the normal of the facet pointing towards this smallest ray.
 422 Otherwise, i.e., if $\alpha_i \alpha_j < 0$,
 423 the interior of $K$ will lie on one side
 424 of the facet and then we choose $\vec y$ to lie on the other side.
 425 That is, $\sp{{\tilde n}_{ij}}y > 0$, for
 426 $\vec {\tilde n}_{ij}$ the normal of the facet pointing away from the cone $K$.
 427 Figure~\ref{fig:primal:examples} shows some example decompositions with
 428 an explicitly marked $\vec y$.
 429
 430 \begin{figure}
 431 \begin{align*}
 432 \intercol=0.48cm
 433 \begin{xy}
 434 <\intercol,0pt>:<0pt,\intercol>::
 435 \POS(-2,-1)="a"*+!U{+}
 436 \POS(2,0)="b"*+!L{+}
 437 \POS(0,2)="c"*+!D{+}
 438 \POS"a"\ar@{-}@[|(3)]"b"
 439 \POS"b"\ar@{-}@[|(3)]"c"
 440 \POS"c"\ar@{-}@[|(3)]"a"
 441 \POS(0.3,0.6)*{\bullet},*+!L{\vec y}
 442 \end{xy}
 443 & =
 444 \intercol=0.48cm
 445 \begin{xy}
 446 <\intercol,0pt>:<0pt,\intercol>::
 447 \POS(2,0)="b"*+!L{+}
 448 \POS(0,2)="c"*+!D{+}
 449 \POS(0,0)="w"*+!DR{\vec w}
 450 \POS"b"\ar@{-}@[|(3)]"c"
 451 \POS"b"\ar@{-}@[|(3)]"w"
 452 \POS"c"\ar@{-}@[|(3)]"w"
 453 \POS(0.3,0.6)*{\bullet},*+!L{\vec y}
 454 \end{xy}
 455 +
 456 \begin{xy}
 457 <\intercol,0pt>:<0pt,\intercol>::
 458 \POS(-2,-1)="a"*+!U{+}
 459 \POS(0,2)="c"*+!D{+}
 460 \POS(0,0)="w"*+!DR{\vec w}
 461 \POS"c"\ar@{-}@[|(3)]"a"
 462 \POS"a"\ar@{-}@[|(3)]"w"
 463 \POS"c"\ar@{--}@[|(3)]"w"
 464 \POS(0.3,0.6)*{\bullet},*+!L{\vec y}
 465 \end{xy}
 466 +
 467 \begin{xy}
 468 <\intercol,0pt>:<0pt,\intercol>::
 469 \POS(-2,-1)="a"*+!U{+}
 470 \POS(2,0)="b"*+!L{+}
 471 \POS(0,0)="w"*+!DR{\vec w}
 472 \POS"a"\ar@{-}@[|(3)]"b"
 473 \POS"a"\ar@{--}@[|(3)]"w"
 474 \POS"b"\ar@{--}@[|(3)]"w"
 475 \POS(0.3,0.6)*{\bullet},*+!L{\vec y}
 476 \end{xy}
 477 \\
 478 \intercol=0.48cm
 479 \begin{xy}
 480 <\intercol,0pt>:<0pt,\intercol>::
 481 \POS(-2,-1)="a"*+!U{+}
 482 \POS(2,0)="b"*+!L{+}
 483 \POS(0,2)="c"*+!D{+}
 484 \POS"a"\ar@{--}@[|(3)]"b"
 485 \POS"b"\ar@{-}@[|(3)]"c"
 486 \POS"c"\ar@{--}@[|(3)]"a"
 487 \POS(-2.5,-1.5)*{\bullet},*+!U{\vec y}
 488 \end{xy}
 489 & =
 490 \intercol=0.48cm
 491 \begin{xy}
 492 <\intercol,0pt>:<0pt,\intercol>::
 493 \POS(2,0)="b"*+!L{+}
 494 \POS(0,2)="c"*+!D{+}
 495 \POS(0,0)="w"*+!DR{\vec w}
 496 \POS"b"\ar@{-}@[|(3)]"c"
 497 \POS"b"\ar@{--}@[|(3)]"w"
 498 \POS"c"\ar@{--}@[|(3)]"w"
 499 \POS(-2.5,-1.5)*{\bullet},*+!U{\vec y}
 500 \end{xy}
 501 +
 502 \begin{xy}
 503 <\intercol,0pt>:<0pt,\intercol>::
 504 \POS(-2,-1)="a"*+!U{+}
 505 \POS(0,2)="c"*+!D{+}
 506 \POS(0,0)="w"*+!DR{\vec w}
 507 \POS"c"\ar@{--}@[|(3)]"a"
 508 \POS"a"\ar@{--}@[|(3)]"w"
 509 \POS"c"\ar@{-}@[|(3)]"w"
 510 \POS(-2.5,-1.5)*{\bullet},*+!U{\vec y}
 511 \end{xy}
 512 +
 513 \begin{xy}
 514 <\intercol,0pt>:<0pt,\intercol>::
 515 \POS(-2,-1)="a"*+!U{+}
 516 \POS(2,0)="b"*+!L{+}
 517 \POS(0,0)="w"*+!DR{\vec w}
 518 \POS"a"\ar@{--}@[|(3)]"b"
 519 \POS"a"\ar@{-}@[|(3)]"w"
 520 \POS"b"\ar@{-}@[|(3)]"w"
 521 \POS(-2.5,-1.5)*{\bullet},*+!U{\vec y}
 522 \end{xy}
 523 \\
 524 \intercol=0.48cm
 525 \begin{xy}
 526 <\intercol,0pt>:<0pt,\intercol>::
 527 \POS(-2,-1)="a"*+!U{+}
 528 \POS(2,0)="b"*+!L{+}
 529 \POS(0,2)="c"*+!D{+}
 530 \POS"a"\ar@{--}@[|(3)]"b"
 531 \POS"b"\ar@{-}@[|(3)]"c"
 532 \POS"c"\ar@{-}@[|(3)]"a"
 533 \POS(1,-1.5)*{\bullet},*+!L{\vec y}
 534 \end{xy}
 535 & =
 536 \intercol=0.48cm
 537 \begin{xy}
 538 <\intercol,0pt>:<0pt,\intercol>::
 539 \POS(2,0)="b"*+!L{+}
 540 \POS(0,2)="c"*+!D{+}
 541 \POS(0,0)="w"*+!DR{\vec w}
 542 \POS"b"\ar@{-}@[|(3)]"c"
 543 \POS"b"\ar@{--}@[|(3)]"w"
 544 \POS"c"\ar@{-}@[|(3)]"w"
 545 \POS(1,-1.5)*{\bullet},*+!L{\vec y}
 546 \end{xy}
 547 +
 548 \begin{xy}
 549 <\intercol,0pt>:<0pt,\intercol>::
 550 \POS(-2,-1)="a"*+!U{+}
 551 \POS(0,2)="c"*+!D{+}
 552 \POS(0,0)="w"*+!DR{\vec w}
 553 \POS"c"\ar@{-}@[|(3)]"a"
 554 \POS"a"\ar@{--}@[|(3)]"w"
 555 \POS"c"\ar@{--}@[|(3)]"w"
 556 \POS(1,-1.5)*{\bullet},*+!L{\vec y}
 557 \end{xy}
 558 +
 559 \begin{xy}
 560 <\intercol,0pt>:<0pt,\intercol>::
 561 \POS(-2,-1)="a"*+!U{+}
 562 \POS(2,0)="b"*+!L{+}
 563 \POS(0,0)="w"*+!DR{\vec w}
 564 \POS"a"\ar@{--}@[|(3)]"b"
 565 \POS"a"\ar@{-}@[|(3)]"w"
 566 \POS"b"\ar@{-}@[|(3)]"w"
 567 \POS(1,-1.5)*{\bullet},*+!L{\vec y}
 568 \end{xy}
 569 \\
 570 \intercol=0.48cm
 571 \begin{xy}
 572 <\intercol,0pt>:<0pt,\intercol>::
 573 \POS(-2,-1)="a"*+!U{0}
 574 \POS(2,0)="b"*+!L{+}
 575 \POS(0,2)="c"*+!D{+}
 576 \POS"a"\ar@{-}@[|(3)]"b"
 577 \POS"b"\ar@{-}@[|(3)]"c"
 578 \POS"c"\ar@{-}@[|(3)]"a"
 579 \POS(1,0.2)*{\bullet},*+!R{\vec y}
 580 \end{xy}
 581 & =
 582 \intercol=0.48cm
 583 \begin{xy}
 584 <\intercol,0pt>:<0pt,\intercol>::
 585 \POS(-2,-1)="a"*+!U{0}
 586 \POS(2,0)="b"*+!L{+}
 587 \POS(1,1)="w"*+!DL{\vec w}
 588 \POS"a"\ar@{-}@[|(3)]"b"
 589 \POS"a"\ar@{-}@[|(3)]"w"
 590 \POS"b"\ar@{-}@[|(3)]"w"
 591 \POS(1,0.2)*{\bullet},*+!R{\vec y}
 592 \end{xy}
 593 +
 594 \begin{xy}
 595 <\intercol,0pt>:<0pt,\intercol>::
 596 \POS(-2,-1)="a"*+!U{0}
 597 \POS(0,2)="c"*+!D{+}
 598 \POS(1,1)="w"*+!DL{\vec w}
 599 \POS"c"\ar@{-}@[|(3)]"a"
 600 \POS"a"\ar@{--}@[|(3)]"w"
 601 \POS"c"\ar@{-}@[|(3)]"w"
 602 \POS(1,0.2)*{\bullet},*+!R{\vec y}
 603 \end{xy}
 604 \\
 605 \intercol=0.48cm
 606 \begin{xy}
 607 <\intercol,0pt>:<0pt,\intercol>::
 608 \POS(-2,-1)="a"*+!U{0}
 609 \POS(2,0)="b"*+!U{+}
 610 \POS(0,2)="c"*+!D{-}
 611 \POS"a"\ar@{-}@[|(3)]"b"
 612 \POS"b"\ar@{--}@[|(3)]"c"
 613 \POS"c"\ar@{-}@[|(3)]"a"
 614 \POS(1.5,1.5)*{\bullet},*+!D{\vec y}
 615 \end{xy}
 616 & =
 617 \intercol=0.48cm
 618 \begin{xy}
 619 <\intercol,0pt>:<0pt,\intercol>::
 620 \POS(-2,-1)="a"*+!U{0}
 621 \POS(0,2)="c"*+!D{-}
 622 \POS(4,-2)="w"*+!L{\vec w}
 623 \POS"c"\ar@{-}@[|(3)]"a"
 624 \POS"a"\ar@{-}@[|(3)]"w"
 625 \POS"c"\ar@{--}@[|(3)]"w"
 626 \POS(1.5,1.5)*{\bullet},*+!D{\vec y}
 627 \end{xy}
 628 -
 629 \begin{xy}
 630 <\intercol,0pt>:<0pt,\intercol>::
 631 \POS(-2,-1)="a"*+!U{0}
 632 \POS(2,0)="b"*+!U{+}
 633 \POS(4,-2)="w"*+!L{\vec w}
 634 \POS"a"\ar@{--}@[|(3)]"b"
 635 \POS"a"\ar@{-}@[|(3)]"w"
 636 \POS"b"\ar@{--}@[|(3)]"w"
 637 \POS(1.5,1.5)*{\bullet},*+!D{\vec y}
 638 \end{xy}
 639 \\
 640 \intercol=0.48cm
 641 \begin{xy}
 642 <\intercol,0pt>:<0pt,\intercol>::
 643 \POS(-2,-1)="a"*+!U{0}
 644 \POS(2,0)="b"*+!U{+}
 645 \POS(0,2)="c"*+!D{-}
 646 \POS"a"\ar@{--}@[|(3)]"b"
 647 \POS"b"\ar@{--}@[|(3)]"c"
 648 \POS"c"\ar@{-}@[|(3)]"a"
 649 \POS(4.7,-2.5)*{\bullet},*+!R{\vec y}
 650 \end{xy}
 651 & =
 652 \intercol=0.48cm
 653 \begin{xy}
 654 <\intercol,0pt>:<0pt,\intercol>::
 655 \POS(-2,-1)="a"*+!U{0}
 656 \POS(0,2)="c"*+!D{-}
 657 \POS(4,-2)="w"*+!L{\vec w}
 658 \POS"c"\ar@{-}@[|(3)]"a"
 659 \POS"a"\ar@{--}@[|(3)]"w"
 660 \POS"c"\ar@{--}@[|(3)]"w"
 661 \POS(4.7,-2.5)*{\bullet},*+!R{\vec y}
 662 \end{xy}
 663 -
 664 \begin{xy}
 665 <\intercol,0pt>:<0pt,\intercol>::
 666 \POS(-2,-1)="a"*+!U{0}
 667 \POS(2,0)="b"*+!U{+}
 668 \POS(4,-2)="w"*+!L{\vec w}
 669 \POS"a"\ar@{-}@[|(3)]"b"
 670 \POS"a"\ar@{--}@[|(3)]"w"
 671 \POS"b"\ar@{--}@[|(3)]"w"
 672 \POS(4.7,-2.5)*{\bullet},*+!R{\vec y}
 673 \end{xy}
 674 \end{align*}
 675 \caption{Examples of decompositions in primal space.}
 676 \label{fig:primal:examples}
 677 \end{figure}
 678
 679 To see that there is a $\vec y$ satisfying the above constraints,
 680 we need to show that $R \cap S$ is non-empty, with
 681 $S = \{ \vec y \mid \sp{{\tilde n}_{i_kj_k}}y > 0 \text{ for all $k$}\}$.
 682 It will be easier to show this set is non-empty when the $\vec u_i$ form
 683 an orthogonal basis.  Applying a non-singular linear transformation $T$
 684 does not change the decomposition of $\vec w$ in terms of the $\vec u_i$ (i.e., the
 685 $\alpha_i$ remain unchanged), nor does this change
 686 any of the scalar products in the constraints that define $R \cap S$
 687 (the normals are transformed by $\left(T^{-1}\right)^\T$).
 688 Finding a vector $\vec y \in T(R \cap S)$ ensures that
 689 $T^{-1}(\vec y) \in R \cap S$.
 690 Without loss of generality, we can therefore assume for the purpose of
 691 showing that $R \cap S$ is non-empty that
 692 the $\vec u_i$ indeed form an orthogonal basis.
 693
 694 In the orthogonal basis, we have $\vec b_i^* = \vec u_i$
 695 and the corresponding inward normal $\vec N_i$ is either
 696 $\vec u_i$ or $-\vec u_i$.
 697 Furthermore, each normal of a facet of $S$ of the first type is of the
 698 form $\vec {\tilde n}_{i_kj_k} = a_k \vec u_{i_k} - b_k \vec u_{j_k}$, with
 699 $a_k, b_k > 0$ and ${i_k} < {j_k}$,
 700 while for the second type each normal is of the form
 701 $\vec {\tilde n}_{i_kj_k} = -a_k \vec u_{i_k} - b_k \vec u_{j_k}$, with
 702 $a_k, b_k > 0$.
 703 If $\vec {\tilde n}_{i_kj_k} = a_k \vec u_{i_k} - b_k \vec u_{j_k}$
 704 is the normal of a facet of $S$
 705 then either
 706 $(\vec N_{i_k}, \vec N_{j_k}) = (\vec u_{i_k}, \vec u_{j_k})$
 707 or
 708 $(\vec N_{i_k}, \vec N_{j_k}) = (-\vec u_{i_k}, -\vec u_{j_k})$.
 709 Otherwise, the facet would not cut $R$.
 710 Similarly,
 711 if $\vec {\tilde n}_{i_kj_k} = -a_k \vec u_{i_k} - b_k \vec u_{j_k}$
 712 is the normal of a facet of $S$
 713 then either
 714 $(\vec N_{i_k}, \vec N_{j_k}) = (\vec u_{i_k}, -\vec u_{j_k})$
 715 or
 716 $(\vec N_{i_k}, \vec N_{j_k}) = (-\vec u_{i_k}, \vec u_{j_k})$.
 717 Assume now that $R \cap S$ is empty, then there exist
 718 $\lambda_k, \mu_i \ge 0$ not all zero
 719 such that
 720 $\sum_k \lambda_k \vec {\tilde n}_{i_kj_k} + \sum_l \mu_i \vec N_i = \vec 0$.
 721 Assume $\lambda_k > 0$ for some facet of the first type.
 722 If $\vec N_{j_k} = -\vec u_{j_k}$, then $-b_k$ can only be canceled
 723 by another facet $k'$ of the first type with $j_k = i_{k'}$, but then
 724 also $\vec N_{j_{k'}} = -\vec u_{j_{k'}}$.  Since the $j_k$ are strictly
 725 increasing, this sequence has to stop with a strictly positive coefficient
 726 for the largest $\vec u_{j_k}$ in this sequence.
 727 If, on the other hand, $\vec N_{i_k} = \vec u_{i_k}$, then $a_k$ can only
 728 be canceled by the normal of a facet $k'$ of the second kind
 729 with $i_k = j_{k'}$, but then
 730 $\vec N_{i_{k'}} = -\vec u_{i_{k'}}$ and we return to the first case.
 731 Finally, if $\lambda_k > 0$ only for normals of facets of the second type,
 732 then either $\vec N_{i_k} = -\vec u_{i_k}$ or $\vec N_{j_k} = -\vec u_{j_k}$
 733 and so the coefficient of one of these basis vectors will be strictly
 734 negative.
 735 That is, the sum of the normals will never be zero and
 736 the set $R \cap S$ is non-empty.
 737
 738 For each ray $\vec u_j$ of cone $K_i$, i.e., the cone with $\vec u_i$ replaced
 739 by $\vec w$, we now need to determine whether the facet not containing this
 740 ray is closed or not.  We denote the (inward) normal of this cone by
 741 $\vec n_{ij}$.  Note that cone $K_j$ (if it appears in (\ref{eq:triangulations}),
 742 i.e., $\alpha_j \ne 0$) has the same facet opposite $\vec u_i$
 743 and its normal $\vec n_{ji}$ will be equal to either $\vec n_{ij}$ or
 744 $-\vec n_{ij}$, depending on whether we are dealing with an ``external'' facet,
 745 i.e., a facet of $K'$, or an ``internal'' facet.
 746 If, on the other hand, $\alpha_j = 0$, then $\vec n_{ij} = \vec n_{0j}$.
 747 If $\sp{n_{ij}}y > 0$, then the facet is closed.
 748 Otherwise it is open.
 749 It follows that the two (or more) occurrences of external facets are either all open
 750 or all closed, while for internal facets, exactly one is closed.
 751
 752 First consider the facet not containing $\vec u_0 = \vec w$.
 753 If $\alpha_i > 0$, then $\vec u_i$ and $\vec w$ are on the same side of the facet
 754 and so $\vec n_{i0} = \vec n_{0i}$.  Otherwise, $\vec n_{i0} = -\vec n_{i0}$.
 755 Second, if $\alpha_j = 0$, then replacing $\vec u_i$ by $\vec w$ does not
 756 change the affine hull of the facet and so $\vec n_{ij} = \vec n_{0j}$.
 757 Now consider the case that $\alpha_i \alpha_j < 0$, i.e., $\vec u_i$
 758 and $\vec u_j$ are on the same side of the hyperplane through the other rays.
 759 If we project $\vec u_i$, $\vec u_j$ and $\vec w$ onto a plane orthogonal
 760 to the ridge through the other rays, then the possible locations of $\vec w$
 761 with respect to $\vec u_i$ and $\vec u_j$ are shown in Figure~\ref{fig:w:same}.
 762 If both $\vec n_{0i}$ and $\vec n_{0j}$ are closed then $\vec y$ lies in region~1
 763 and therefore $\vec n_{ij}$ (as well as $\vec n_{ji}$) is closed too.
 764 Similarly, if both $\vec n_{0i}$ and $\vec n_{0j}$ are open then so is
 765 $\vec n_{ij}$.  If only one of the facets is closed, then, as explained above,
 766 we choose $\vec n_{ij}$ to be open, i.e., we take $\vec y$ to lie in region~3
 767 or~5.
 768 Figure~\ref{fig:w:opposite} shows the possible configurations
 769 for the case that $\alpha_i \alpha_j > 0$.
 770 If exactly one of $\vec n_{0i}$ and $\vec n_{0j}$ is closed, then
 771 $\vec y$ lies in region~3 or region~5 and therefore $\vec n_{ij}$ is closed iff
 772 $\vec n_{0j}$ is closed.
 773 Otherwise, as explained above, we choose $\vec n_{ij}$ to be closed if $i < j$.
 774
 775 \begin{figure}
 776 \intercol=0.7cm
 777 \begin{minipage}{0cm}
 778 \begin{xy}
 779 <\intercol,0pt>:<0pt,\intercol>::
 780 \POS(-4,0)="a"\ar@[|(2)]@{-}(4,0)
 781 \POS?(0)/\intercol/="b"\POS(2,0)*\xybox{"b"-"a":(0,0)\ar^{\vec n_{0i}}(0,0.75)}
 782 \POS(2,0)*{\bullet},*+!U{\vec u_j}
 783 \POS(-2,-3)="a"\ar@[|(2)]@{-}(2,3)
 784 \POS?(0)/\intercol/="b"\POS(1,1.5)*\xybox{"b"-"a":(0,0)\ar^{\vec n_{0j}}(0,-0.75)}
 785 \POS(1,1.5)*{\bullet},*+!R{\vec u_i}
 786 \POS(-2,3)="a"\ar@{-}(2,-3)
 787 \POS?(0)/\intercol/="b"\POS(1.5,-2.25)
 788 *\xybox{"b"-"a":(0,0)\ar_{\vec n_{ji}}^{\vec n_{ij}}(0,+0.75)}
 789 \POS(1.5,-2.25)*{\bullet},*+!R{\vec u_0 = \vec w}
 790 \POS(0,0)*{\bullet}
 791 \POS(3,1.5)*+[o][F]{\scriptstyle 1}
 792 \POS(0,2.5)*+[o][F]{\scriptstyle 2}
 793 \POS(-3,1.5)*+[o][F]{\scriptstyle 3}
 794 \POS(-3,-1.5)*+[o][F]{\scriptstyle 4}
 795 \POS(0,-3)*+[o][F]{\scriptstyle 5}
 796 \POS(3,-1.5)*+[o][F]{\scriptstyle 6}
 797 \end{xy}
 798 \end{minipage}
 799 \begin{minipage}{0cm}
 800 \begin{xy}
 801 <\intercol,0pt>:<0pt,\intercol>::
 802 \POS(-4,0)="a"\ar@[|(2)]@{-}(4,0)
 803 \POS?(0)/\intercol/="b"\POS(2,0)*\xybox{"b"-"a":(0,0)\ar^{\vec n_{0i}}(0,0.75)}
 804 \POS(2,0)*{\bullet},*+!U{\vec u_j}
 805 \POS(-2,-3)="a"\ar@[|(2)]@{-}(2,3)
 806 \POS?(0)/\intercol/="b"\POS(1,1.5)*\xybox{"b"-"a":(0,0)\ar^{\vec n_{0j}}(0,-0.75)}
 807 \POS(1,1.5)*{\bullet},*+!R{\vec u_i}
 808 \POS(-2,3)="a"\ar@{-}(2,-3)
 809 \POS?(0)/\intercol/="b"\POS(-1.5,2.25)
 810 *\xybox{"b"-"a":(0,0)\ar_{\vec n_{ji}}^{\vec n_{ij}}(0,+0.75)}
 811 \POS(-1.5,2.25)*{\bullet},*+!R{\vec u_0 = \vec w}
 812 \POS(0,0)*{\bullet}
 813 \POS(3,1.5)*+[o][F]{\scriptstyle 1}
 814 \POS(0,2.5)*+[o][F]{\scriptstyle 2}
 815 \POS(-3,1.5)*+[o][F]{\scriptstyle 3}
 816 \POS(-3,-1.5)*+[o][F]{\scriptstyle 4}
 817 \POS(0,-3)*+[o][F]{\scriptstyle 5}
 818 \POS(3,-1.5)*+[o][F]{\scriptstyle 6}
 819 \end{xy}
 820 \end{minipage}
 821 \caption{Possible locations of $\vec w$ with respect to $\vec u_i$ and
 822 $\vec u_j$, projected onto a plane orthogonal to the other rays, when
 823 $\alpha_i \alpha_j < 0$.}
 824 \label{fig:w:same}
 825 \end{figure}
 826
 827 \begin{figure}
 828 \intercol=0.7cm
 829 \begin{minipage}{0cm}
 830 \begin{xy}
 831 <\intercol,0pt>:<0pt,\intercol>::
 832 \POS(-4,0)="a"\ar@[|(2)]@{-}(4,0)
 833 \POS?(0)/\intercol/="b"\POS(2,0)*\xybox{"b"-"a":(0,0)\ar^{\vec n_{0i}}(0,0.75)}
 834 \POS(2,0)*{\bullet},*+!U{\vec u_j}
 835 \POS(-2,3)="a"\ar@[|(2)]@{-}(2,-3)
 836 \POS?(0)/\intercol/="b"\POS(-1,1.5)*\xybox{"b"-"a":(0,0)\ar^{\vec n_{0j}}(0,+0.75)}
 837 \POS(-1,1.5)*{\bullet},*+!R{\vec u_i}
 838 \POS(-2,-3)="a"\ar@{-}(2,3)
 839 \POS?(0)/\intercol/="b"\POS(1.5,2.25)
 840 *\xybox{"b"-"a":(0,0)\ar^{\vec n_{ji}}(0,+0.75)
 841 \POS(0,0)\ar_{\vec n_{ij}}(0,-0.75)}
 842 \POS(1.5,2.25)*{\bullet},*+!L{\vec u_0 = \vec w}
 843 \POS(0,0)*{\bullet}
 844 \POS(3,1.5)*+[o][F]{\scriptstyle 1}
 845 \POS(0,2.5)*+[o][F]{\scriptstyle 2}
 846 \POS(-3,1.5)*+[o][F]{\scriptstyle 3}
 847 \POS(-3,-1.5)*+[o][F]{\scriptstyle 4}
 848 \POS(0,-3)*+[o][F]{\scriptstyle 5}
 849 \POS(3,-1.5)*+[o][F]{\scriptstyle 6}
 850 \end{xy}
 851 \end{minipage}
 852 \begin{minipage}{0cm}
 853 \begin{xy}
 854 <\intercol,0pt>:<0pt,\intercol>::
 855 \POS(-4,0)="a"\ar@[|(2)]@{-}(4,0)
 856 \POS?(0)/\intercol/="b"\POS(2,0)*\xybox{"b"-"a":(0,0)\ar^{\vec n_{0i}}(0,0.75)}
 857 \POS(2,0)*{\bullet},*+!U{\vec u_j}
 858 \POS(-2,3)="a"\ar@[|(2)]@{-}(2,-3)
 859 \POS?(0)/\intercol/="b"\POS(-1,1.5)*\xybox{"b"-"a":(0,0)\ar^{\vec n_{0j}}(0,+0.75)}
 860 \POS(-1,1.5)*{\bullet},*+!R{\vec u_i}
 861 \POS(-2,-3)="a"\ar@{-}(2,3)
 862 \POS?(0)/\intercol/="b"\POS(-1.5,-2.25)
 863 *\xybox{"b"-"a":(0,0)\ar^{\vec n_{ji}}(0,+0.75)
 864 \POS(0,0)\ar_{\vec n_{ij}}(0,-0.75)}
 865 \POS(-1.5,-2.25)*{\bullet},*+!L{\vec u_0 = \vec w}
 866 \POS(0,0)*{\bullet}
 867 \POS(3,1.5)*+[o][F]{\scriptstyle 1}
 868 \POS(0,2.5)*+[o][F]{\scriptstyle 2}
 869 \POS(-3,1.5)*+[o][F]{\scriptstyle 3}
 870 \POS(-3,-1.5)*+[o][F]{\scriptstyle 4}
 871 \POS(0,-3)*+[o][F]{\scriptstyle 5}
 872 \POS(3,-1.5)*+[o][F]{\scriptstyle 6}
 873 \end{xy}
 874 \end{minipage}
 875 \caption{Possible locations of $\vec w$ with respect to $\vec u_i$ and
 876 $\vec u_j$, projected onto a plane orthogonal to the other rays, when
 877 $\alpha_i \alpha_j > 0$.}
 878 \label{fig:w:opposite}
 879 \end{figure}
 880
 881 The algorithm is summarized in Algorithm~\ref{alg:closed}, where
 882 we use the convention that in cone $K_i$, $\vec u_i$ refers to
 883 $\vec u_0 = \vec w$.
 884 Note that we do not need any of the rays or normals in this code.
 885 The only information we need is the closedness of the facets in the
 886 original cone and the signs of the $\alpha_i$.
 887
 888 \begin{algorithm}
 889 \begin{tabbing}
 890 next \= next \= next \= \kill
 891 if $\alpha_j = 0$ \\
 892 \> closed[$K_i$][$\vec u_j$] := closed[$\tilde K$][$\vec u_j$] \\
 893 else if $i = j$ \\
 894 \> if $\alpha_j > 0$ \\
 895 \> \> closed[$K_i$][$\vec u_j$] := closed[$\tilde K$][$\vec u_j$] \\
 896 \> else \\
 897 \> \> closed[$K_i$][$\vec u_j$] := $\lnot$closed[$\tilde K$][$\vec u_j$] \\
 898 else if $\alpha_i \alpha_j > 0$ \\
 899 \> if closed[$\tilde K$][$\vec u_i$] = closed[$\tilde K$][$\vec u_j$] \\
 900 \> \> closed[$K_i$][$\vec u_j$] := $i < j$ \\
 901 \> else \\
 902 \> \> closed[$K_i$][$\vec u_j$] := closed[$\tilde K$][$\vec u_j$] \\
 903 else \\
 904 \> closed[$K_i$][$\vec u_j$] := closed[$\tilde K$][$\vec u_i$] and
 905 closed[$\tilde K$][$\vec u_j$]
 906 \end{tabbing}
 907 \caption{Determine whether the facet opposite $\vec u_j$ is closed in $K_i$.}
 908 \label{alg:closed}
 909 \end{algorithm}
 910
 911 \subsection{Triangulation in primal space}
 912 \label{s:triangulation}
 913
 914 As in the case for Barvinok's decomposition (Section~\ref{s:decomposition}),
 915 we can transform a triangulation of a (closed) cone into closed simple cones
 916 into a triangulation of half-open simple cones that fully partitions the
 917 original cone, i.e., such that the half-open simple cones do not intersect at their
 918 facets.
 919 Again, we apply Proposition~\ref{p:inclusion-exclusion} with $\vec y$
 920 an interior point of the cone (Section~\ref{s:interior}).