%\documentclass[12pt]{article}
%\setlength{\oddsidemargin}{0.25in} 
%\setlength{\evensidemargin}{0.25in}
%\setlength{\topmargin}{-0.5in} 
%\setlength{\headheight}{0.25in}
%\setlength{\headsep}{0.25in}
%\setlength{\textheight}{8.75in}
%\setlength{\textwidth}{6in} 
%\setlength{\footskip}{0.5in}
%%\renewcommand{\baselinestretch}{1.5}
%\setlength{\parindent}{0dd}
%%\pagestyle{myheadings}
%\markright{\today}
%\setlength{\parskip}{1.5ex plus 0.5ex minus 0.5ex}
%\usepackage{amsfonts}
%\usepackage{epsfig}
%\newcommand{\field}[1]{\mathbb{#1}}

%\begin{document}

\chapter{Regions and Lattices}
\section{Creating {\tt Region} Objects}

The {\tt Region} class in Opt provides a convenient way of describing
regions in space.
Regions are built from a selection of primitive shapes and using the
various set operations.

A new {\tt Region} object is created by calling
\begin{verbatim}
reg = Region(type, parameters)
\end{verbatim}
{\tt type} is one of the primitive shape names, which at this point
include {\tt sphere}, {\tt halfspace} and {\tt convpoly}. {\tt
parameters} is a structure whose fields are the various parameters
that describe the shape. Each shape requires the field {\tt dim} which
indicates the dimension of the region. 

In addition, {\tt sphere} requires the fields {\tt center} and {\tt
radius}. {\tt center} is a vector of real numbers giving the coordinates
of the center of the sphere object. The radius of the sphere is a real
number given in {\tt radius}.

A halfspace is the set described by
\[ \{ \mathbf{x} \in \field{R}^n \: | \: \mathbf{a'x} \geq b \} \]
and thus requires the fields {\tt a} and {\tt b}.

{\tt convpoly} is short for convex polytope. A {\tt convpoly} object
is created as the convex hull of the list of points given in {\tt
points}. {\tt points} is a matrix with {\tt dim} columns, with each row
giving the coordinates for a separate point. The {\tt convpoly} is more
easily represented internally as an intersection of halfspaces, or the
set
\[ \{ \mathbf{x} \in \field{R}^n \: | \: \mathbf{Ax} \geq \mathbf{b} \} \]
which is determined as follows.

First the convex hull is  determined using {\sc Matlab}'s {\tt convhulln}
routine, which returns a list of the facets which make up the surface of
the polytope. A facet in $N$ dimensions is an ($N$-1)-dimensional
simplex, each of which lies on a hyperplane. A hyperplane is the set
\[ \{ \mathbf{x} \in \field{R}^n \: | \: \mathbf{a'x} = b \} \]
or 
\[ \{ \mathbf{x} \in \field{R}^n \: | \: \sum_{i=1}^N a_i x_i = b \} \]
Let $\mathbf{x}^1 \ldots \mathbf{x}^N$ be the points which determine a
facet, and thus its hyperplane. $x_k^n$ is the $k$th coordinate of the
$n$th point. Then the coefficients of the hyperplane
are given by the following determinants:
%
\[ a_1 = \left|
\begin{array}{llll}
1 & x^1_2 & \cdots & x^1_N \\
\vdots & \vdots & \ddots & \vdots \\
1 & x^N_2 & \cdots & x^N_N 
\end{array}
\right|
\qquad 
a_2 = \left|
\begin{array}{lllll}
x^1_1 & 1 & x^1_3 & \cdots & x^1_N \\
\vdots & \vdots & \vdots & \ddots & \vdots \\
x^N_1 & 1 & x^N_3 & \cdots & x^N_N 
\end{array}
\right|
\qquad \cdots \]
%
\[
a_N = \left|
\begin{array}{llll}
 x^1_1 & \cdots & x^1_{N-1} & 1 \\
\vdots & \ddots & \vdots & \vdots \\
 x^N_1 & \cdots & x^N_{N-1} & 1
\end{array}
\right| \qquad
b = \left| 
\begin{array}{lll}
x^1_1 & \cdots & x^1_N \\
\vdots & \ddots & \vdots \\
x^N_1 & \cdots & x^N_N 
\end{array}
\right|
\]
This hyperplane is the boundary of one of the halfspaces whose
intersection describes the polytope. The direction of the inequality of
the corresponding halfspace is determined by checking points in other
facets of the polytope. It is possible that there are multiple coplanar
facets of the polytope. Points of a coplanar facet should give an
inconclusive result, but because of numerical error can give an
incorrect answer. Because of this, all of the points on the surface of
the convex hull are checked, and direction is chosen as that indicated
by the majority of the points.

\subsection{Forming Composite Regions}
Composite {\tt Region} objects may be formed by using the familiar set
operations such as {\tt union}, {\tt intersect} and {\tt not}. For ease
of expressing unions and intersections of a series of objects, the
operators {\tt +} and {\tt *} may be substituted for {\tt union} and
{\tt intersect}, respectively. {\tt setdiff(a,b)} returns the region
which is the set difference $a / b$, or $a \cap \neg b$. 

A composite region object is always formed as a result of a binary
operation, and therefore contains the data indicating the type of
operation and the two operands. Thus, a complicated composite {\tt
Region} can be thought of as a binary operation tree whose leaves are
primitive shapes.

The operation {\tt not(a)} or $\sim${\tt a} where {\tt a} is a primitive
{\tt Region} returns a primitive region object identical to {\tt a} but
with the complement flag set. If, however, {\tt a} is a {\tt composite}
region, the complement operation is applied recursively down through the
operation tree.

\subsection{Simple Region Transformations}
Two functions are available to perform simple transformations on a {\tt
Region} object. 

\subsubsection{Offset}
A region offset is done using the {\tt or} function. This is more
conveniently called using the equivalent operator {\tt |}. The syntax,
for an $n$-dimensional region, is
\begin{verbatim}
oreg = reg | offset
\end{verbatim}
where {\tt offset} is a vector of length $n$, giving the offset amount
in each coordinate direction. 

Shifting a {\tt sphere} is straightforward; its center is simply
adjusted by the offset amount.

A {\tt halfspace}, as mentioned earlier, is described by
\[ \{ \mathbf{x} \in \field{R}^n \: | \: \mathbf{a'x} \geq b \} \]
and has the hyperplane $\mathbf{a}'\mathbf{x} = b$ as its boundary.
$\mathbf{a}$ can be thought of as the hyperplanes's normal vector, and
$b$ as the distance of the hyperplane from the origin along the normal
vector $\mathbf{a}$. Shifting a hyperplane has the effect of sliding it
along its normal vector. The distance slid is the length of the
projection of the offset vector along the normal direction, or
$\mathbf{a} \cdot \mathbf{offset}$. Therefore the $b$ parameter must be
adjusted as follows:
\[ b_{\mbox{off}} = b + \mathbf{a} \cdot \mathbf{offset} \]

Since a {\tt convpoly} is simply an intersection of several halfspaces,
the procedure for shifting it is the same; each element of the {\tt
convpoly}'s $\mathbf{b}$ vector is adjusted as follows:
\[ b_{k, \mbox{off}} = b_k + \mathbf{a_k} \cdot \mathbf{offset} \]
where each vector $\mathbf{a_k}$ is a row of the $\mathbf{A}$ matrix. In
addition, the extreme points of the {\tt convpoly} are shifted by the
same vector, like the center point of a {\tt sphere}.

\subsubsection{Rotation}
A region can be rotated using the {\tt rot} function. The syntax is
\begin{verbatim}
rreg = rot(reg, ang, dims)
\end{verbatim}
This function rotates {\tt reg} by {\tt ang} radians in the
2-dimensional subspace indicated by the 2-vector {\tt dims}. For
example, {\tt rot(reg, pi/4, [1 3])} will rotate {\tt reg} 45 degrees in
the $xz$-subspace; that is, about the $y$-axis.

For a {\tt sphere}, the center point is rotated by an angle $\theta$ by
premultiplying the vector of coordinates in the specified 2-dimensional
subspace by the rotation matrix
\[ \mathbf{Q}_{\theta} = \left[
\begin{array}{rr}
\cos\theta & -\sin\theta \\
\sin\theta & \cos\theta
\end{array}
\right] . \]

For a {\tt halfspace} with boundary hyperplane $\mathbf{a}'\mathbf{x} =
b$, the rotation is accomplished similarly by premultiplying the components
$\mathbf{a}$ vector corresponding to the specified 2-dimensional
subspace by $\mathbf{Q}_{\theta}$. Since $\mathbf{a}$ is the normal
vector of the hyperplane, it can be rotated like any other vector, while
$b$, the distance of the hyperplane from the origin along the normal
vector, remains constant.

For a {\tt convpoly}, which is the intersection of halfspaces, each
component halfspace is treated the same as above. In addition, the
extreme points of the {\tt convpoly} are rotated as well, like the
center point of the {\tt sphere}.


If either of these functions is called on a {\tt composite} region, the
operation is applied recursively down through the operation tree so that
each primitive object is offset or rotated by the same amount.

\section{Creating {\tt Lattice} Objects}
The {\tt Lattice} class in Opt is used to specify a lattice, or pattern
of points. A lattice in $n$ dimensions is described by its basis matrix
$\mathbf{M} \in \field{R}^{n \times m}$. $\mathbf{M}$ is post-multiplied
by an integer vector $\mathbf{\xi} \in \field{Z}^n$ to generate the
coordinates of the lattice points in cartesian space. While $m$ can be
any integer greater than $n$, it is most convenient to require that
$m=n$. A lattice can be scaled and offset as well. The full expression
for a point on a lattice is then
\[ \mathbf{x} = \alpha\mathbf{M\xi} + \mathbf{c} \]
where $\alpha \in \field{R}$ and $\mathbf{c} \in \field{R}^n$.

A {\tt Lattice} object is created by calling
\begin{verbatim}
lat = Lattice(M)
\end{verbatim}
where {\tt M} is the $n \times n$ basis matrix. The {\tt Lattice} object
will also contain a scale factor, preset to 1, and an offset, preset to
the zero-vector, both of which can be adjusted as follows.

A {\tt Lattice} can be scaled by using the {\tt *} operator. The syntax
\begin{verbatim}
slat = a * lat}
\end{verbatim}
scales {\tt lat} by a scalar {\tt a} and returns the {\tt Lattice}
object {\tt slat}. The object {\tt slat} will contain the same basis
matrix as well as the scale factor.

A {\tt Lattice} can be offset by using the {\tt +} operator. The syntax
\begin{verbatim}
olat = lat + off
\end{verbatim}
shifts {\tt lat} by the amounts in the vector {\tt off} and returns the
result in {\tt olat}. As with the {\tt *} operator, {\tt olat} will have
the same basis matrix as {\tt lat}, with only the offset data adjusted.

\section{{\tt Region} and {\tt Lattice} Operations}
\subsection{{\tt isin} Function}
The function {\tt isin} takes a list of points and returns those located
inside a {\tt Region}. The syntax of the function is
\begin{verbatim}
result = isin(reg, points) 
\end{verbatim}
Both {\tt result} and {\tt points} are matrices in which each row gives
the coordinates of a point.

If {\tt reg} is a region primitive, the result is determined by checking
each input point against the rules for that type of shape.

 A point
$\mathbf{x} = (x_1, \ldots, x_N)$ is inside a {\tt sphere} of radius
$\rho$ and center $\mathbf{c}$ if
\[ \sum_{i=1}^N (x_i - c_i)^2 \leq \rho^2 \]
or 
\[ || \mathbf{x} - \mathbf{c} || \leq \rho \]

A point $\mathbf{x}$ is in a {\tt halfspace} if 
\[ \mathbf{a'x} \geq b \]

A point $\mathbf{x}$ is in a {\tt convpoly} if
\[ \mathbf{Ax} \geq \mathbf{b} \]
where $\mathbf{A}$ and $\mathbf{b}$ are determined as described above.

If {\tt reg} is a composite region, then it was formed as a result of
set operations on primitive regions. Each of the two child regions is
queried, and the result for a given point is computed. For example, if
the two child regions were combined by an {\tt intersect} operation, the
a point is in the region only if the point is both of the child regions.
If a child region is itself a composite, the process continues until a
primitive region is queried. The results for the component primitive
regions are then combined using set operations to produce the result for
the composite region.

\subsection{{\tt msop} Function}
\label{msop}
The function {\tt msop} takes a composite {\tt Region} as input and
determines a simplified sum-of-products (or union-of-intersections)
equivalent. This is useful in determining whether a region is bounded,
and if so, finding a bounding parallelpiped guaranteed to contain the
entire region.

Any function on primitive regions $R_1, \ldots, R_m$ using the operations
$\cup$, $\cap$ and $\neg$ can be reduced to a sum-of-products form, for
example
\[ f(R_1, \ldots, R_m) = \bigcup_{i=1}^p \bigcap_{j \in [1,m]}
	\delta_{i,j} R_m \]

Each composite {\tt Region} object stores the identity of each unique
primitive object that it is composed of. The function {\tt probett} is
then used to assemble the entire truth table of the Boolean function
that the composite {\tt Region} represents. The function {\tt qm} is
then used to find a reduced sum-of-products form for the truth table
using the Quine-McCluskey method of finding prime implicants.

The reduced function is returned as a cell array, of which each element
represents a product term. Each product term is represented as a vector
of numbers, each of which is the index number of a primitive region in
the composite. An example output is
\begin{verbatim}
fcn =
    
    [2 3 5] [1 2]
\end{verbatim}
which indicates that the simplified representation of the composite
region is
\[ (R_2 \cap R_3 \cap R_5) \cup (R_1 \cap R_2) \]

{\tt msop} is used by the function {\tt bb}.

\subsection{{\tt bb} Function}

The function {\tt bb} determines whether a {\tt Region} object is
bounded, and if so, returns its bounding box. It is necessary for
determing the lattice points which lie in a given region. Therefore, the
function takes two arguments: a {\tt Region} object and a {\tt Lattice}
object. The bounding box is returned not as Cartesian coordinates, but
in terms of intervals along the lattice basis vectors. Thus the bounding
box is in fact a parallelepiped which is guaranteed to circumscribe the
{\tt Region} in question.

For a {\tt convpoly} region of dimension $n$, the bounding
parallelepiped is determined by testing each of the region's extreme
points.

A bounding interval must be found for each lattice basis vector. For
each basis vector, any of the extreme points could contribute to the
bounding interval, and as such each extreme point must be tested. 

Once the bounding interval is found for a given basis vector, two faces
of the bounding parallelepiped can be constructed. Each face of the
bounding parallelepiped lies in the hyperplane spanned by all of the
other basis vectors, and is itself a parallelepiped of dimension $n-1$.
Since we require that the lattice basis matrix be full-rank, any $n-1$
basis vectors will span an $(n-1)$-dimensional subspace, and will have a
1-dimensional nullspace. This nullspace can be seen as the normal vector
of the $(n-1)$-dimensional subspace. Together, the nullspace and the
current basis vector lie in a 2-dimensional subspace (although they may
be collinear, in which case they would not span the plane).

Each extreme point $\mathbf{p}$ is projected onto the orthonormal basis
to the nullspace, $\mathbf{u}$, and the minimum and maximum projections
determine the bounding interval. The projection for a given extreme
point and lattice basis vector is
\[ \mbox{proj}_{\mathbf{u}} \mathbf{p} = \frac{ \mathbf{p} \cdot
\mathbf{u} }{|\mathbf{u}|^2} \mathbf{u} = (\mathbf{p} \cdot \mathbf{u})
\mathbf{u} \]
The angle $\theta$ between this projection and the current basis vector
$\mathbf{v}$ is given by
\[ \cos \theta = \frac{  \mbox{proj}_{\mathbf{u}} \mathbf{p}  \cdot \mathbf{v} }{ |\mbox{proj}_{\mathbf{u}} \mathbf{p}|
|\mathbf{v}|} \]
 
Together, $\mbox{proj}_{\mathbf{u}} \mathbf{p}$ and $\theta$ determine
the distance $r'$ along the direction of the basis vector
$\mathbf{v}$. By the definition of cosine,
\[ \cos \theta = \frac{|\mbox{proj}_{\mathbf{u}} \mathbf{p}|}{r'} \]
and rearranging gives
\begin{eqnarray*}
 r' 
&=&
 \frac{ |\mbox{proj}_{\mathbf{u}} \mathbf{p}|^2 |\mathbf{v}| }
      { (\mbox{proj}_{\mathbf{u}} \mathbf{p} \cdot \mathbf{v}) } \\ 
&=& 
\frac{ (\mathbf{p} \cdot \mathbf{u})^2 |\mathbf{v}|}
     { [(\mathbf{p} \cdot \mathbf{u}) \mathbf{u}] \cdot \mathbf{v} }
\end{eqnarray*}

$r'$ must be scaled by the norm of the basis vector to determine its
distance along that vector, giving
\[ r = \frac{r'}{|\mathbf{v}|} = \frac{  (\mathbf{p} \cdot \mathbf{u})^2}{ [(\mathbf{p} \cdot \mathbf{u}) \mathbf{u}] \cdot \mathbf{v}
} \]

If it turns out that the projection onto the nullspace is the
zero-vector, the distance along the lattice basis vector zero, and the
computation is not done, or a division-by-zero error will result.

If the basis for the nullspace and the lattice basis vector are
collinear, the angle between them is 0, making 
\[ r = \frac{ | \mbox{proj}_{\mathbf{u}} \mathbf{p} | }{ |\mathbf{v}| } \]


%\begin{figure}[h]
%\begin{center}
%\epsfxsize=4in\epsfbox{art/bb.eps}
%\end{center}
%\caption{Determining the bounding box for a polyhedron.}
%\label{bb}
%\end{figure} 

Because lattice points are determined by multiplying an integer vector
by the basis matrix, the interval is rounded away from zero.

A {\tt sphere} region's bounding box is most easily determined by
inscribing it inside a hypercube, which has $2^n$ extreme points. The
hypercube's bounding parallelepiped is then determined the same way as
the {\tt convpoly}'s.

A {\tt halfspace}, as well as a complemented {\tt sphere} and {\tt
convpoly} are always unbounded. However, the intersection of at least
$n+1$ halfspaces is possible bounded. This possibility is covered by the
treatment of {\tt composite} regions.

The bounding parallelepiped of a {\tt composite} region is found by
properly combining the bounding parallelepipeds of the constituent
regions. This is accomplished by first obtaining a sum-of-products (or
union-of-intersections) representation of the {\tt composite} region,
using the {\tt msop} function described in section \ref{msop}.

Each product term is the intersection of  primitive regions. If
any of the product terms is unbounded, the entire region, which is the
union of all of the product terms, will be unbounded. Otherwise, the
bounding intervals are found by taking the minimum of the axis minima of
each product term and the maximum of the axis maxima of each product
term.

For an individual product term, the bounding intervals are found by
taking the maximum of the axis minima of each primitive region, and the
minimum of the axis maxima of each primitive region. If at least $n+1$
{\tt halfspaces} are present in a product term, it is possible that
their intersection is bounded.

A halfspace, and intersections of halfspaces, are convex polytopes. To
determine whether a convex polytope is bounded, it is necessary to solve
two linear programs, as follows.

For a convex polytope 
\[ P = \{ \mathbf{x} \in \field{R}^n \: | \:
\mathbf{Ax} \geq \mathbf{b} \} \] 
the {\em recession cone}  is the set of all directions $\mathbf{d}$ along
which we can move indefinitely away from a point inside $P$. The
recession cone is the set
\[ \{ \mathbf{d} \in \field{R}^n \: | \: \mathbf{Ad} \geq \mathbf{0} \} \]

If $P$ is bounded, the recession cone will consist of only the
zero-vector, which is always a solution to $\mathbf{Ad} \geq
\mathbf{0}$. Therefore, if either of the two linear programs
\begin{eqnarray*}
\min & & \pm \mathbf{c}'\mathbf{x} \\
\mbox{subject to} & & \mathbf{Ax} \geq \mathbf{0}
\end{eqnarray*}
has a solution other than the zero-vector, a non-trivial recession cone
exists, and the set $P$ is unbounded. Conversely, if the only solution
to both linear programs is the zero-vector, no non-trivial recession
cone exists, and the set $P$ is bounded. 

It is possible for one choice of the cost vector $\mathbf{c}$ for a
false indication of boundedness to result. If the unbounded direction is
along the hyperplane $\mathbf{c}'\mathbf{y} = \mathbf{0}$ it is possible
for a linear program solver to return the zero-vector as an optimal
solution for both linear programs, as any solution along the hyperplane
has minimum cost. However, if $\mathbf{c}$ is chosen randomly, the
probability of this occurring is extremely small. If not for the finite
precision of numbers on a computer, the probability would in fact be zero.

The linear program solvers recommended for use with {\tt Opt} require
that the linear program be given in {\em standard form}; that is,
\[ \begin{array}{rcl}
\min & & 
\begin{array}{rcl} \mathbf{c}'\mathbf{x} & & \end{array}
	\\
\mbox{subject to} & & 
	\begin{array}[t]{rcl} 
		\mathbf{Ax} &=& \mathbf{b} \\
		\mathbf{x} &\geq& \mathbf{0} 
	\end{array}
\end{array} \]

The duals of our linear programs are
\[ \begin{array}{rcl}
\min & & 
\begin{array}{rcl} \mathbf{0} \cdot \mathbf{y} & & \end{array}
	\\
\mbox{subject to} & & 
	\begin{array}[t]{rcl} 
		\mathbf{A'}\mathbf{y} &=& \pm \mathbf{c} \\
		\mathbf{y} &\geq& \mathbf{0} 
	\end{array}
\end{array} \]
which are in standard form.
By duality theory, if the dual of a linear program is infeasible, the
primal is either infeasible or unbounded. As we know that the primal in
our case is feasible, it must be unbounded. If the dual is feasible and
has a finite optimum, then the primal must also be feasible and have a
finite optimum. 

Therefore, to test the boundedness of the intersection of halfspaces, we
concatenate their $\mathbf{a}$ vectors to form the matrix $\mathbf{A}$,
and we concatenate their $b$'s to form the vector $\mathbf{b}$. We then
run both of the dual linear programs. If either is infeasible, the
polytope $P$ formed by the intersection of halfspaces is unbounded. If
both are feasible, and have finite optima, $P$ is bounded. If it is
determined that $P$ is bounded, the extreme points of $P$ are found, and
we follow the procedure detailed above for finding the bounding
parallelepiped of a convex polytope. 

A vector $\mathbf{p}$ is an extreme point of $P$ if
\begin{enumerate}
\item All equality constraints are active.
\item At least $n$  linearly independent constraints are active at
$\mathbf{p}$.
\item All constraints are satisfied.
\end{enumerate}
Our representation of the convex polytope contains no equality
constraints, so item 1 is automatically satisfied. We test all possible
combinations of $n$ linearly independent constraints, and solve the
linear system
\[ \mathbf{a_i}'\mathbf{x} = b_i, \qquad i \in I_k \]
where the $I_k$ is an index set over the $k$th combination of $n$
linearly independent constraints. If the solution $\mathbf{x}$
satistfies all of the other constraints as well, it is an extreme point
of the polytope.

One final point: If the {\tt Lattice} given is offset, the offset will
be applied to the extreme points before the bounding intervals are
found. Therefore the bounding intervals which are returned by the
function describe a parallelepiped which circumscribes a {\tt Region}
offset by the same amount. 

\subsection{{\tt inbox} function}
The function {\tt inbox} returns a list of points on a given lattice
that lie within a bounding box (parallelepiped). The syntax is
\begin{verbatim}
locs = inbox(lat, box)
\end{verbatim}
where {\tt lat} is a {\tt Lattice} object and {\tt box} is a vector of
the form
\begin{verbatim}
box = [1min 1max 2min 2max ... ]
\end{verbatim}
with entries which are the minimum and maximum extent of the
parallelepiped along each of the lattice basis vectors. The points in
{\tt locs} are returned as a matrix in which each row lists the
Cartesian coordinates of a single point. 

To generate the list of points, we produce an integer grid within the
bounding intervals. Each point on this integer grid is used as
$\mathbf{\xi}$ in the expression
\[ \mathbf{x} = \alpha\mathbf{M\xi} + \mathbf{c} \]
If the {\tt Lattice} is offset, the offset $\mathbf{c}$ will be applied
after the points are calculated, and so the bounding intervals must be
given in terms of an offset region in space. This, as described above,
is automatically taken into account by the {\tt bb} function.

\subsection{{\tt points} function}
The {\tt points} function return a list of the points on a given lattice
that lie within a given region. The syntax is
\begin{verbatim}
locs = points(reg, lat)
\end{verbatim}
The function works by calling three previously-described functions as
follows:
\begin{enumerate}
\item {\tt bb} is called to determine the bounding parallelepiped for the
region;
\item {\tt inbox} is called to find the lattice points within the
bounding parallelepiped;
\item {\tt isin} is called to determine which of the points returned by
{\tt inbox} are in fact within the region
\end{enumerate}

\section{Usage of {\tt Region}s and {\tt Lattice}s in {\tt Opt}}
Two useful applications of {\tt Region} and {\tt Lattice} to {\tt Opt}
are
\begin{itemize}
\item Easy formation of optimizable sequences in a variety of shapes and
patterns
\item Easy specification of frequency responses
\end{itemize}

An $n$-dimensional optimizable sequence can be created using the command
\begin{verbatim}
seq = optArray(lattice, region, pool)
\end{verbatim}
where {\tt pool} is a handle to a pool of optimizables. This command
creates a sequence with support on lattice points within the given
region. If the region and lattice given are not of the same
dimensionality, or if the region given is unbounded, an error is
returned. Otherwise, we compute the bounding parallelepiped, and then
determine which points within that parallelepiped are inside {\tt
region}.

In addition, if {\tt seq} has already been specified as an \optarr\:,
the command
\begin{verbatim}
sseq = seq(reg)
\end{verbatim}
creates a subsequence {\tt sseq} which contains only those points
located within {\tt reg}.

For frequency gridding, a convenient procedure is to define a region of
the desired shape and a lattice of sufficient tightness, and use the
{\tt points} function to output a grid of point locations.

\subsection{Examples}
\label{sect.circimpulse}

As an example, let us create a 2-D filter with its impulse response on a
circular region of support. Such an impulse response pattern is useful
when a circularly-symmetric frequency response is desired. In addition,
the impulse response is to be real and symmetric about the origin so
that the frequency response will be symmetric as well.

To produce an impulse response with the desired symmetry properties,
we must take care that the coefficient at each point $(n_1, n_2)$ depend on
the same optimization variables as its reflection at $(-n_1,
-n_2)$. First we define a square {\tt Lattice} of points from which the
impulse response locations will be chosen.
\begin{verbatim}
sqlat = Lattice(eye(2)); 
\end{verbatim}
Then, we define a Region enclosing a single quadrant
\begin{verbatim}
p1.dim = 2;
p1.points = [0 0; 0 N; N 0; N N];
quadrantBoundary = Region('convpoly', p1);
\end{verbatim}
where integer parameter {\tt N} can be increased if better performance
is desired. The following code defines a circular region centered at the
origin with radius {\tt N}.
\begin{verbatim}
p2.dim = 2;
p2.center = [0 0];
p2.radius = N;
circleReg = Region('sphere', p2);
\end{verbatim}
Next, we create the quarter-circle region by intersecting region {\tt
circleReg} with {\tt quadrantBoundary}.
\begin{verbatim}
qCircle = quadrantBoundary * circleReg;
\end{verbatim}
The following creates two \optarr{}s with support on the points on the
lattice {\tt sqlat} within the region {\tt qCircle}, with {\tt X} a
space of optimization variables:
\begin{verbatim}
aa1 = optArray(sqlat, qCircle, X);
aa2 = optArray(sqlat, qCircle, X);
\end{verbatim}
Finally, we combine various reflections of the two impulse responses to
create the whole circular array:
\begin{verbatim}
aa = aa1 + aa1' + flip(aa2, 1) + flip(aa2, 2);
\end{verbatim}
The command {\tt flip(a, dim)} flips an \optarr\ in dimension {\tt dim};
that is, all point locations are changed from $x= \{\ldots, x_{\mbox{dim}},
\ldots\}$ to $\hat{x} = \{ \ldots, -x_{\mbox{dim}}, \ldots \}$. So {\tt
aa1} is in the first quadrant, {\tt aa1'} in the third, {\tt flip(aa2,
1)} in the second and {\tt flip(aa2, 2)} in the fourth. Although points
along the axes will depend on two optimization variables each, 

The array is shown in \rfig{art.circArray}.

\begin{figure}[h]
%\centerline{\epsfxsize=4in\epsfbox{art/circArray.eps}}
\begin{center}
\includegraphics[width=4.0in]{art/circArray.eps}
\end{center}
\caption{Circular Array.}
\label{art.circArray}
\end{figure}


%\end{document}
