% $Id: autodif.tex 2762 2014-12-05 21:45:49Z johnoel $
\documentclass{admbmanual}
\makeatletter\@twosidefalse\makeatother
\hypersetup{urlcolor=black}

% Some local definitions.
%
% Do not use a blank line between the question and the answer.
\newcommand{\question}[1]{\bigskip\noindent{\bf #1?}\par\medskip\noindent}

\newcommand\nr{{Newton-Raphson}}

% Some legacy shortcuts.
\newcommand\pd[2]{\frac{\partial #1}{\partial #2}}
\newcommand\bs{$\backslash$}
\newcommand\adtypes{\texttt{dvariable}, \texttt{dvar\_array},
  \texttt{dvar\_matrix}, \texttt{dvector}, and \texttt{dmatrix}}
\newcommand\admbversion{11.2}

% Set the plot symbol size for PicTeX pictures.
\renewcommand\setplotsymbol{\fiverm .}

% Make a table with commands followed by their descriptions.
\newenvironment{commandtable}
  {
    \par
    \begin{tabular}%
      {@{\vrule height 12pt depth7pt width0pt}@{\tt-}l @{\rm\quad} l}
  }
  % Use a & to separate command from description.  End line with a \\
  {
    \end{tabular}
    \medskip
   }

% Make a table with days of the week followed by up to 10 columns of numbers
% vertically aligned on their decimal points.  Wanted it in place, so don't use
% the table environment.
\newenvironment{daysdecimaltable}
  {
  %\begin{table}[!h]
    \par
    \medskip
    \begin{tabular}{l @{\qquad} *{10}{r@{}l} }
  }
  % Usage:  (have to put decimal in, because some numbers don't have them.
  % Monday  &  25&.2 &  -41 & & 433 & .12 &  10 & .11 & 23& \\
  % Tuesday & 34&.21 &  356 & &  23 &     &  23 & .1  &  5& \\
  {
    \end{tabular}
    \bigskip
    \par
 % \end{table}
  }

% Title definition.
\title{%
  \largetitlepart{AUTODIF}
  \smalltitlepart{%
    A \cplus\ Array Language Extension\\
    with Automatic Differentiation\\
    for Use in Nonlinear Modeling and Statistics}
  \vspace{4ex}\textsf{\textit{Version \admbversion~~(2013-05-01)\\[3pt]
      Revised manual~~(2013-06-24)}}\vspace{3ex}
}
 % Author definition.
\author{\textsf{\textit{David Fournier}}}
% Should be U/lc.
\manualname{Autodif}

% We want an index.
\makeindex

% Start the thing.
\begin{document}

% We want a title page.
\maketitle

~\vfill
\noindent ADMB Foundation, Honolulu.\\\\
\noindent This is the manual for AUTODIF version \admbversion.\\\\
\noindent Copyright \copyright\ 1991, 1992, 1993, 1994, 1996, 2000, 2008, 2011,
2013 David Fournier\\\\
\noindent The latest edition of the manual is available at:\\
\url{http://admb-project.org/documentation/manuals}

% We want a table of contents.
\tableofcontents

\chapter{Introduction}
\label{ch:introduction}

\X{dvariable@\texttt{dvariable}}
\X{dvar\_vector@\texttt{dvar\_vector}}
\question{What is \scAD}
\scAD\ is an array language extension to \cplus\ which allows the user to
automatically calculate the partial derivatives of a function of many
independent variables by simply redeclaring the usual C floating point types
doubles and arrays of doubles to be the corresponding \scAD\ type:
\texttt{dvariable}, and \texttt{dvar\_vector}. The user can write the code for
the function using any valid \cplus\ syntax. \scAD\ does not impose any
restrictions on allowable code constructions. Using \scAD, the programmer can
manipulate vectors and matrices as single objects in complicated mathematical
calculations, and at the same time, obtain derivatives of the resulting function
with respect to the independent variables.

When calculating derivatives, there is no restriction on calling user-supplied
functions, which may themselves call other functions---so long as the types have
been redefined to be the corresponding \scAD\ types. Since \cplus\ is upwardly
compatible with \textsc{ansi C}, existing numerical routines written in C can be
easily modified to use \scAD.

\question{Who should use \scAD}
Anyone who would like to have the derivatives of a complicated function computed
automatically without any additional programming effort will find \scAD\ useful.
It is particularly useful for optimization problems involving a differentiable
function of many independent variables. Such problems occur, for example, in
nonlinear statistical modeling, i.e., in nonlinear parameter estimation, in
sensitivity analysis, and in finding roots of systems of equations.

\question{Why is automatic differentiation better than estimating the
  derivatives by finite differences}
\X{automatic differentiation}%
There are two reasons why the \scAD\ method is better than finite differences
for calculating partial derivatives.
\XX{accuracy}{of derivative calculations}
\begin{enumerate}

  \item With \scAD, the derivatives are calculated to the limit of machine
  accuracy, that is, to the same degree of accuracy as the function itself. The
  derivatives are obtained as accurately as they would be if an analytical
  expression for the derivatives were evaluated on the computer. With finite
  differences, the accuracy of the derivative calculations seldom approaches
  that of the limit of the machine accuracy for calculating the derivatives.
  Consequently, optimization schemes based on automatic differentiation tend to
  be more stable and perform better than those that use finite difference
  approximations.

  \item To obtain estimates of the partial derivatives by finite difference
  approximations, the amount of calculation necessary is proportional to the
  number of independent variables times the amount of time needed to calculate
  the function itself. With \scAD, the amount of calculation is less than 5
  times the amount of calculation needed to calculate the function itself
  \cite{griewankcorliss1991}. In particular, the computing time does not depend
  on the number of independent variables. This makes \scAD\ particularly
  attractive for calculating derivatives for functions with many independent
  variables.
\end{enumerate}
\XX{speed}{of calculations}

\question{Is it necessary to learn \cplus\ before one can use \scAD}
Any experienced C programmer can learn to use \scAD\ in a few hours. While it is
necessary to have a \cplus\ compiler to use \scAD, only a small number of extra
concepts besides those used in the C language need to be learned. This is due to
the fact that \cplus\ is essentially a superset of C, so valid
\textsc{ansi}-compliant C code needs only a very small amount of modification to
be converted to valid \cplus\ code.

We feel that the best way to get a feel for how \scAD\ is used is to study a
series of examples. These examples have been chosen to illustrate the use of
\scAD\ on problems in different fields from engineering to statistical analysis.

\question{How does automatic differentiation work}
Every computer program for calculating a function value, no matter how
complicated, can be broken down into a series of binary and unary mathematical
operations. Automatic differentiation is simply the chain rule from elementary
calculus applied to this series. The following simple example illustrates the
chain rule:
\begin{align}\label{eq:introduction-01}
   \nonumber u&=x^2+y^2\\
             v&=x+3y\\
   \nonumber f&=\sin(u)+\exp(v)
\end{align}

We want to calculate the derivatives of the function $f$ with respect to the
variables~$x$ and~$y$. By the chain rule,
\begin{equation*}
  \pd fx=\pd fu\pd ux + \pd fv\pd vx\qquad\qquad
          \pd fy=\pd fu\pd uy + \pd fv\pd vy
\end{equation*}
Now
 \begin{gather}\label{eq:introduction-02}
  \pd fu=\cos(u) \qquad\qquad\qquad\qquad  \pd fv=\exp(v)\\[6pt]
           \pd ux=2x \qquad \pd vx=1 \qquad\qquad    \pd uy=2y\qquad\pd vy=3
 \end{gather}
\noindent such that
\begin{equation*}
  \pd fx=\cos(u)*2x+\exp(v)*1\qquad\qquad\pd fy=\cos(u)*2y+\exp(v)*3
\end{equation*}
That is really all there is to it! Your program is broken down into a series of
simple operations like these so that the derivatives can be calculated. The
important point is that no matter how complicated the mathematical expressions
are, they can be thought of as consisting of a long series of binary and unary
operations.

\question{What is adjoint (precompiled) derivative code and why is it so
  important} While the reverse mode of automatic differentiation is extremely
efficient in terms of the number of calculations required to compute
derivatives, the method requires a large amount of temporary storage. Each
arithmetic operation used to calculate the original function generates 28--30
bytes of temporary storage. Consider the fact that inverting a 100-by-100 matrix
requires $10^6$ operations. This would generate 28--30 megabytes of temporary
storage requirements. Clearly, it would be impractical to carry out large
calculations under these conditions. \X{precompiled derivative code}

The use of precompiled derivative code combines the best of two worlds, so to
speak. The small number of arithmetic operations required by the reverse mode
for calculating derivatives is retained while, at the same time, the amount of
temporary storage required is greatly reduced. The amount of storage required
for inverting a 100-by-100 matrix inverse is reduced to about 500K. Even better,
the amount of temporary storage requirement generated by adding two matrices is
about 32 bytes, no matter how large the matrices are.

Precompiled derivative code means that the actual derivative calculations are
written in \cplus\ functions. Following object oriented design considerations,
the precompiled code is completely ``encapsulated'' so that it is invisible to
the user. This means that it makes no difference to the user whether or not a
particular function has a precompiled derivative code component written for it.
At the user level, everything looks the same---except that performance is
greatly enhanced. The use of precompiled derivative code produces the ultimate
in reuseable code. For example, writing good code for the derivative of the
inverse of a matrix is not a trivial procedure. However, once it has been
written, any user of \scAD\ can take advantage of it simply by taking the
inverse of a matrix object with code like
\begin{lstlisting}
  bbM1=inv(M);  // M1 and M are dvar_matrix objects
\end{lstlisting}

\scAD\ combines array and matrix classes with the use of precompiled derivative
code for common array and matrix calculations and assignments to produce a
superior environment for efficient derivative calculations.

\section{A simple example of a derivative calculation}

While \scAD\ has been supplied with all the computational power described above,
it has been designed so that, at the user level, it appears like enhanced C
code, or perhaps as enhanced C code with an additional array language. Before
illustrating some of \scAD's array language facilities, we will begin with a
simple C~code example.

The following C code for the function \texttt{fcomp}
calculates the sum of squares
\begin{equation*}
  \sum_{i=0}^{n-1} (x_i-1)^2
\end{equation*}
This code has been written in C without using any \cplus\ or \scAD\ extensions,
to show how C code can be modified to use \scAD\ for calculating derivatives.

\begin{lstlisting}
#include <math.h>

double fcomp(int n,double * x)
{
  double z;
  double tmp;
  int i;
  z=0;
  for (i=0;i<n;i++)
  {
    tmp=x[i]-1;
    z=z+tmp*tmp;
  }
  return(z);
}
\end{lstlisting}
\texttt{fcomp} computes the sum of squares of $n$ double precision floating
point numbers, which are contained in an array of length \texttt{n} pointed to
by \texttt{x}. The expression \texttt{x[i]} refers to the
$\texttt{i}^\textrm{th}$ element of the array. (Note that the default
construction for C begins indexing an array with zero.)

\section{Modifying \texttt{fcomp} to use \scAD}

\X{dvariable@\texttt{dvariable}}
The arithmetic in the function \texttt{fcomp} is carried out using objects of
type \texttt{double}, that is, eight byte floating point numbers. In order to
calculate derivatives at the same time as the function is being calculated, the
class \texttt{double} is replaced with a new \scAD\ class, the
\texttt{dvariable}. An object of type \texttt{dvariable} ``appears'' to the user
to be an object of type \texttt{double} in that it can be used in arithmetic and
mathematical expressions, and it has a numerical value that can be displayed.
The difference between an object of type \texttt{double} and an object of type
\texttt{dvariable} is that when arithmetic is done using a \texttt{dvariable}
object, information necessary for the calculation of derivatives is
automatically generated as well.

\X{dvar\_vector@\texttt{dvar\_vector}}
In addition to objects of type \texttt{dvariable}, we shall need a class of
``vector'' object, corresponding to \texttt{x}, that generates derivative
information. In \scAD, this class is the \texttt{dvar\_vector} class.

The following listing is the \texttt{fcomp} function modified to use automatic
differentiation.
\input fcomp_s.cpp
\noindent There are five changes.
\begin{enumerate}
  \item The line \texttt{\#include <fvar.hpp>} has been included. The file
  \texttt{fvar.hpp} contains information about \scAD\ for the \cplus\ compiler.

  \item The array of variables \texttt{x} has been redefined to be of type
  \texttt{dvar\_vector}. This \scAD\ type is like an array of doubles for which
  derivative information is collected. Besides providing for the calculation of
  derivatives, the use of the \texttt{dvar\_vector} class for \texttt{x}
  provides an additional improvement. With the original C code, there was no way
  of determining in the subroutine \texttt{fcomp} whether or not enough memory
  had actually been allocated to the array~\texttt{x}. \scAD\ provides optional
  index bounds checking for all container classes, such as
  \texttt{dvar\_arrays}. This greatly speeds up program development.

  \item The variables \texttt{z} and \texttt{tmp} have been declared to be of
  type \texttt{dvariable}. The \texttt{dvariable} type is the fundamental type
  for automatic differentiation. It corresponds to the C floating point type
  \texttt{double}.

  \item The function \texttt{value} has been included in the line
  \texttt{return(value(z));} The function \texttt{value()} returns the numerical
  part of an \texttt{dvariable}. This is necessary to convert the
  \texttt{dvariable} to a \texttt{double} so that the value can be returned and
  used in subsequent calculations (such as in a function minimizer) without
  generating any more derivative calculations.

  \item The range of the index of the \texttt{for} loop has been changed to
  begin at 1 and include \texttt{n}. While this change is not strictly necessary
  for this example, \scAD\ allows the user to declare arrays that have arbitrary
  valid index ranges. This facility often contributes to the production of
  simpler looking code.
\end{enumerate}

\section{Interfacing with the rest of the program}

The form of the interface of the function \texttt{fcomp} with the rest of the
program will depend on the user's application. This section illustrates two
applications. The first example shows how to simply calculate the derivatives so
they can be used for any other purpose whatsoever. The second example shows how
to invoke a function minimizer, supplied with \scAD, to find the minimum of the
function \texttt{fcomp}. In any event, it is necessary to identify the
independent variables for \scAD. (These are the variables for which the partial
derivatives are calculated.) The independent variables are identified by
declaring them to be of type \texttt{independent\_variables}.
\X{independent variables}
The declaration \texttt{independent\_variables x(1,nvar)} creates an array of 20
independent variables and identifies the independent variables for \scAD. An
object of type \texttt{gradient\_structure} must be declared and remain in scope
until all the derivative calculations have been finished.
\XX{minimization}{of a function}
\X{gradient\_structure@\texttt{gradient\_structure}}
\begin{lstlisting}
#include <math.h>
#include <fvar.hpp>

double fcomp(dvar_vector); // Function prototype declaration
void main()
{
  double f;
  int i;
  int nvar=20;
  independent_variables x(1,nvar);  // Identify the independent variables
  dvector g(1,nvar);  // Holds the vector of partial derivatives (gradient)
  gradient_structure gs;    // Must declare this structure to manage
                            // derivative calculations
  f=fcomp(nvar,x);
  gradcalc(nvar,g);        // The derivatives are calculated
  cout <<" The gradient vector is\n"<<g<<"\n"; // Print out the derivatives
}                                                   // on the screen
\end{lstlisting}

\XX{constructor}{from \texttt{independent\_variables} to \texttt{dvar\_vector}}

Notice that although \texttt{x} has been declared to be an object of type
\texttt{independent\_variables} in the main routine, it is declared to be an
object of type \texttt{dvar\_vector} in \texttt{fcomp}. While this might appear
to be an error, it is actually an important part of the interface. It is
necessary to identify both what the independent variables are and when one
should begin to collect derivative information for calculations involving them.
When the vector \texttt{x} is passed to the function \texttt{fcomp}, the
compiler will invoke a constructor to convert the type
\texttt{independent\_variables} to the type \texttt{dvar\_vector}. At this time,
the steps necessary for identifying the independent variables and initializing
the collection of information needed for calculating the derivatives are carried
out. With the Borland \cplus\ compiler, a warning message
\XX{warning messages}%
{temporary used for type \texttt{xxx} in function \texttt{yyy}}
\begin{lstlisting}
  Temporary used for type xxx in function yyy
\end{lstlisting}
\noindent is printed out, notifying the user that a constructor has been invoked
to perform the necessary conversion. Such a message is usually benign and we
tend to ignore them, since we have tried to design \scAD\ so that unwanted type
conversions do not occur. However, it may indicate that some unintended
conversion has taken place. This is one of the subtler pitfalls of \cplus\
programming. If you want to avoid such warning messages, it is possible to
explicitly ``cast'' the object \texttt{x} to the desired type as
\XX{cast}{from \texttt{independent\_variables} to \texttt{dvar\_vector}}
\begin{lstlisting}
  f=fcomp(nvar,dvar_vector(x));
\end{lstlisting}
\X{BEGIN\_MINIMIZATION@\texttt{BEGIN\_MINIMIZATION}}
\X{END\_MINIMIZATION@\texttt{END\_MINIMIZATION}}
\XX{macros}{BEGIN\_MINIMIZATION@\texttt{BEGIN\_MINIMIZATION}}
\XX{macros}{END\_MINIMIZATION@\texttt{END\_MINIMIZATION}}

The next example shows how the above code is modified to invoke the function
minimizer routine. The function to be minimized and the function minimizer are
embedded in a loop so that it is not necessary for the user's function to call
the function minimizer or for the function minimizer to call the user's
function. This eliminates the need to communicate with either routine by passing
global variables. Since the code that sets up the loop is independent of the
user's function to be minimized, we have defined two macros,
\texttt{BEGIN\_MINIMIZATION} and \texttt{END\_MINIMIZATION}, which will expand
to produce the loop structure. The first macro, \texttt{BEGIN\_MINIMIZATION}
takes five arguments: the number of variables in the function to be minimized,
the name of the vector of independent variables, the name of the dependent
variable, the name of the vector of partial derivatives, and the name of the
\texttt{fmm} control structure. It is the user's responsibility to declare these
five objects. The second macro, \texttt{END\_MINIMIZATION} takes two arguments:
the number of variables and the name of the vector of partial derivatives.
\input simple.cpp

It is not necessary to use the macros. They have simply been provided for
convenience. If you prefer, their expansion can be written out, as in the
following example:
\begin{lstlisting}
#include <math.h>
#include <fvar.hpp>

double fcomp(dvar_vector); // Function prototype declaration
void main()
{
  int nvar=20; // This is the number of independent variables
  independent_variables x(1,nvar); // these are the independent variables
  double f;  // This is the dependent variable
  dvector g(1,nvar);  // Holds the vector of partial derivatives (gradient)
  fmm fmc(nvar);      // Create structure to manage minimization
  // Expansion of BEGIN_MINIMIZATION(nvar,x,f,g,fmc) macro
  gradient_structure gs;
  while (fmc.ireturn >= 0)   // Begin loop for minimization
  {
    fmc.fmin(f,x,g);         // Calls the minimization routine
    if (fmc.ireturn > 0)     // Loop for evaluating the function and
    {                        // derivatives
    // End of the Expansion of BEGIN_MINIMIZATION macro
      f=fcomp(nvar,x);
    // Expansion of END_MINIMIZATION(nvar,g) macro
      gradcalc(nvar,g);
    }
  }
  // End of the Expansion of END_MINIMIZATION macro
  cout << " The minimizing values are\n" << x << "\n"; // Print out answer
}
\end{lstlisting}

The minimization routine is controlled by the object \texttt{fmc} of type
\texttt{fmm}. Different aspects of the minimization---such as the maximum number
of function evaluations before stopping, convergence criteria, and printing
options---can be adjusted by the user through this object. See
Chapter~\ref{ch:function-minimizing}, on the function minimizers, for more
details. The flow through the loops is controlled by the
\texttt{integer fmc.ireturn}.

If you wish to use the conjugate gradient function minimizer instead of the
quasi-Newton function minimizer, you should declare the object \texttt{fmc} to
be of type \texttt{fmmc} rather than type \texttt{fmm}. Nothing else needs to be
changed.
\XX{function minimizer}{conjugate gradient}
\X{conjugate gradient function minimizer}

\section{Using \scAD's built in\br vector and matrix calculations}

\XX{overloading}{of functions in \cplus}
An important feature of \cplus\ is the ability to overload of arithmetic
operations, i.e., extending their definition to user-defined types. For example,
the plus operation (\texttt{+}) can be extended to user-defined string objects,
so if \texttt{string1} and \texttt{string2} are two string objects, then
\begin{lstlisting}
  string1 + string2
\end{lstlisting}
produces a string object that contains the concatenation of \texttt{string1} and
\texttt{string2}. \scAD\ uses this feature to extend the operations of addition,
subtraction, multiplication, and division to various combinations of matrices,
vectors, and numbers. As a result, the matrices and vectors can be employed as
objects in complex mathematical calculations while obtaining derivatives. This
greatly simplifies and speeds up code development. Addition and multiplication
of vectors and matrices appear in many engineering and statistical applications.

\XX{bounds checking}{for container objects}
\X{safe arrays}

This example introduces two new types: the \texttt{dvector} and the
\texttt{dmatrix}. These are the constant analogues of the \texttt{dvar\_vector}
and \texttt{dvar\_matrix} types. They are ``safe array'' objects that contain
their own size information, and for which optional array bounds checking is
provided. The memory for all these objects is managed by \scAD\ in a fashion
that is transparent to the user. \scAD\ can combine these constant and variable
objects automatically so, for example, if \texttt{A} is a \texttt{dmatrix} and
\texttt{M} is a \texttt{dvar\_matrix}, the expression \texttt{A*M} in the user
code will produce an \texttt{dvar\_matrix} that is the product of \texttt{A} and
\texttt{M}. At the same time, the array bounds on \texttt{A} and \texttt{M} will
be checked to ensure that the operation is legal for these two objects.
\X{dvector@\texttt{dvector}}
\X{dvar\_vector@\texttt{dvar\_vector}}
\X{dmatrix@\texttt{dmatrix}}
\X{dvar\_matrix@\texttt{dvar\_matrix}}

To illustrate these vector and matrix operations, we have taken the following
example from~\cite{griewank1988}. The use of a cubic equation of state yields
the Hemholtz energy of a mixed fluid of unit volume at the absolute temperature
$T$ as
\begin{equation*}
  f(x)=RT\sum_{i=1}^n x_i \log\left(\frac{x_i}{1-b^Tx}\right)
-\frac{x^TAx}{\sqrt{8}b^Tx}
   \,\log\left(%
       \frac{1+\big(1+\sqrt{2}\big)b^Tx}%
               {1+\big(1-\sqrt{2}\big)b^Tx}\right)
\end{equation*}
where $R$ is the universal gas constant and
\begin{equation*}
  0\le x,b \in {\mbox{\bf R}}^n\qquad
                      A=A^T\in {\mbox{\bf R}}^{n\times n}
\end{equation*}
During the simulation of an oil reservoir, this function and its derivatives
must be calculated at thousands of points in space and time.

In this example, there are two new applications of the overloaded multiplication
operator (\texttt{*}). $b^Tx$ denotes the dot product of the $n$-dimensional
vector~$b$ with the $n$-dimensional vector~$x$, $Ax$~denotes the multiplication
of the $n$-dimensional vector~$x$ by the $n\times n$~matrix~$A$, and $x^T Ax$ is
the result of the dot product of the vector~$x$ with the vector~$Ax$. Here is
the code for the function that calculates the Hemholtz energy using the
overloaded operators.
\input hem_func.cpp

\X{overloaded operators}
All the vector and matrix multiplications are carried out using the overloaded
operators \texttt{*}, \texttt{+}, and \texttt{-}, which enable the code to be
developed extremely quickly. Notice that multiplication~(using \texttt{*}) can
be carried out for any combination of \texttt{dvector}, \texttt{dmatrix},
\texttt{dvar\_vector}, and \texttt{dvar\_matrix}---so long as the operation
makes sense mathematically. \scAD\ can combine the constant types
\texttt{double}, \texttt{dmatrix}, and \texttt{dvector} with the variable types
\texttt{dvariable}, \texttt{dvar\_vector}, and \texttt{dvar\_matrix} in a manner
transparent to the user.

Here is the code that calls the function \texttt{hemholtz\_energy}.
\input hemholtz.cpp

\section{The log-likelihood function for a\br multivariate normal distribution}

As a further example of how to use \scAD's vector and matrix calculations,
consider the problem of estimating the parameters for a multivariate normal
model when the means and covariance matrix depend on a vector of common
parameters $\Theta$. Let $Y$ be a vector of observations $(y_1, y_2,\ldots,y_n)$
that satisfies the relationship
\begin{equation*}
  Y = A(\Theta)X + \epsilon
\end{equation*}
where $X$ is a vector $(x_1, x_2,\ldots,x_m)$ of controls, $A(\Theta)$ is a
$n\times m$ matrix depending on a parameter vector
$\Theta=(\theta_1,\ldots,\theta_p)$, and $\epsilon=(\epsilon_1,
\epsilon_2,\ldots,\epsilon_n)$ is a random vector with a multivariate normal
distribution, with mean vector $0$ and covariance matrix $\Sigma(\Theta)$.

The log-likelihood function for the parameters $\Theta$ is given by
\begin{equation*}
  -0.5 \ln\Big(\mbox{\rm det}\big(\Sigma(\Theta)\big)\Big)-
   0.5\big(Y-A(\Theta)X\big)^T
   \,\Sigma(\Theta)^{-1}
   \big(Y-A(\Theta)X\big)
\end{equation*}
The maximum likelihood estimates for the parameters are found by maximizing the
log-likelihood function with respect to $\Theta$. This is a nonlinear
optimization problem. The best methods for solving such a problem require the
derivatives of the log-likelihood function with respect to the parameter vector
$\Theta$. Using \scAD's vector and matrix classes, the code for calculating the
log-likelihood function and obtaining the derivatives can be written as follows.
(Actually, there are even more efficient methods of coding these calculations
with \scAD.)
\XX{multivariate normal distribution}{log-likelihood function}

\begin{lstlisting}
dvariable log_likelihood(dvector x,dvar_vector Y,dvar_matrix Sigma,
  dvar_matrix A)
 // It is assumed that Sigma and A have already been calculated by
 // other routines in terms of the parameters Theta
{
  dvariable f;
  dvar_vector diff=Y-A*x;
  f=-0.5*log(det(Sigma))-.5*diff*inv(Sigma)*diff;
  return f;
}
\end{lstlisting}

When this code is executed, the information necessary for calculating the
derivatives is automatically generated without further intervention by the user.
In addition, optional array bounds checking is provided to ensure that all
vector and matrix objects have the right dimensions for the operations
performed.

\question{What are ragged arrays} Traditional programming languages often used
for numerically intensive applications---such as \textsc{fortran}---only support
``flat'' data structures such as matrices. The natural data structures for many
applications simply do not come in this simple form. \scAD\ has extended these
simple forms to ragged matrices and ragged 3-dimensional arrays. A ragged matrix
is a 2-dimensional array where the rows of the array can contain a different
number of elements and have different valid index ranges. A ragged 3-dimensional
array is made up of matrices of different sizes, and these matrices can
themselves be ragged matrices.

Such ragged objects can be used to create many useful data structures. For
example, the ``weights'' in a feed-forward neural network can be best described
by a ragged 3-dimensional array. This can be used to provide an extremely
compact description of the neural network. See Chapter~\ref{ch:neural-networks},
on neural networks, for more details.

Since \cplus\ and \scAD\ are extensible, there is in principle no limit to the
complexity of the data structures that can be developed. This means that you
will never hit a ``dead end'' when developing an application.

%xx \htmlnewfile

\chapter{Getting Started}

\section{System requirements: PC implementations}
\X{system requirements}
\X{hardware}

\scAD\ is designed to operate on 80286, 80386, and 80486-based microcomputers
with a full 640 Kb main memory running under the \textsc{ms-dos} operating
system version 3.2 or higher. A hard disk is not required. However, \scAD\ may
generate several Mb of gradient information to be stored on an external device.
A hard disk and a \textsc{ram} disk are highly recommended for greatest speed.
The \scAD\ libraries have been compiled with either in-line code for a numeric
coprocessor or with code to emulate a numeric coprocessor. A numeric coprocessor
is also highly recommended for greatest speed. The optimum configuration would
include a numeric coprocessor, plus a hard disk and \textsc{ram} disk with at
least 2 Mb free memory, for storing the temporary files created by \scAD. (Large
applications may require even more disk memory.) \scAD\ appears to be compatible
with most popular extended memory managers and multitasking systems.

\X{libraries}
The \scAD\ system consists of a header file \texttt{fvar.hpp} and two or more
libraries. There is a safe library, which provides bounds checking and
initialization of arrays. This library should be used for code development. In
addition, there is an optimized library for faster execution of production
programs. Check the \texttt{read.me} file on your \scAD\ distribution disk for
the names and number of libraries supplied with your version of \scAD.

\section{System requirements: other implementations}

\scAD\ has been ported to a number of other platforms. The necessary information
for each specific platform should be contained in a \texttt{read.me} file.

\section{Installation}

The simplest way to install the \scAD\ system is to place the header file in the
directory where your compiler normally searches for header files, and the
libraries in the directory where your linker normally searches for libraries. If
you choose this option, there is no need to specify additional paths for the
compiler and linker to search for header files and libraries. Alternatively, you
may wish to install the \scAD\ system in its own directory. In this case, you
will need to tell the compiler and linker where to look for the header file and
libraries.

An installation program is included on the distribution disk. This program is
invoked by typing \texttt{a:install} at the \textsc{dos} prompt (assuming the
distribution disk is in drive \texttt{a:}). You will then be asked for the drive
of the source disk (in this case, \texttt{a:}) and the drive where you want to
install \scAD. Next, you will be asked for the directories in which to install
the \scAD\ header file, libraries, and examples. If your \cplus\ compiler is
located in \texttt{c:\bs cc}, you might want to specify \texttt{c:} as the
destination drive and \texttt{\bs cc\bs include} and \texttt{\bs cc\bs lib} as
the directories into which to install the include file and library files,
respectively. If you press \texttt{Enter} for any of these directories, the
corresponding files will be installed in a directory called \texttt{AUTODIF} on
the destination drive.

\textit{Be sure to read the} \texttt{read.me} \textit{file for last minute
  changes and compiler-specific information.}

\section{Setting up the Borland compilers for \scAD}
\X{Borland}
\XX{IDE@\textsc{ide}}{Borland}
\X{memory model}

To compile \scAD\ applications with the integrated development environment
(\textsc{ide}), the following options must be selected within \textsc{ide}:
\begin{itemize}
  \item Select {\bf O}ptions {\bf C}ompiler {\bf C}ode Generation. Choose {\bf
    L}arge memory model. (Optional: select Pre-c{\bf o}mpiled headers if you
  wish to speed up compilation.) Select Mo{\bf r}e for Advanced Code generation.
  Select 8028{\bf 7} to generate direct 80287 inline code or {\bf E}mulation to
  use numeric coprocessor emulation. (\textit{Optional:} Select {\bf D}ebug info
  in \textsc{obj}s if you intend to use the debugger.)

  \item (\textit{Optional:} If you intend to generate an overlay executable,
  select {\bf O}ptions {\bf C}ompiler {\bf E}ntry/Exit Code \textsc{dos}
  {\bf o}verlay.)

  \item Select {\bf O}ptions {\bf D}irectories. If you have not stored your
  \scAD\ header file and library files in the default directories, enter under
  {\bf I}nclude Directories the path where you have stored the \scAD\ header
  file (\texttt{fvar.hpp}), and under {\bf L}ibrary Directories the path where
  you have stored the \scAD\ libraries.

  \item Select {\bf P}roject {\bf O}pen to create a project file. Include your
  source files and the \scAD\ library that you wish to use.
\end{itemize}

\X{make}
To compile \scAD\ applications using a makefile, the \texttt{bcc} or
\texttt{bccx} command-line compilers and \texttt{tlink} or \texttt{tlinkx}
linkers must be used.

\XX{compiler options}{Borland}
\X{memory model}
You should invoke the following compiler and linker options (see the file
\texttt{makefile} supplied with the examples):
\begin{commandtable}
     I  &  followed by the path of your \scAD\ include file.\\
     ml  &  to specify the large memory model. \\
     L  &  followed by the path of your \scAD\ libraries.\\
    c  &  to suppress linking after compilation.\\
    v  &  if you want to use the Borland debugger.\\
    f287  &   to use inline 80287 instructions.\\
\end{commandtable}

\XX{linker options}{Borland}
You should invoke the \texttt{/c} and \texttt{/v} linker options if you are
using debugging.

\XX{linking}{Borland linkers}
The following files must be linked in the following order (see the file
\texttt{makefile} supplied with the examples):
\begin{enumerate}
  \item \texttt{c0l.obj} startup module for the large memory model.

  \item Your source code \texttt{.obj} files.

  \item \texttt{graphics.lib}, if you have a graphic application.

  \item \texttt{ado7.lib} or \texttt{ade.lib}.

  \item \texttt{fp87.lib} for the inline floating-point library or
  \texttt{emu.lib} for the emulation library.

  \item \texttt{mathl.lib}, the math library for the large memory model.

  \item \texttt{cl.lib}, the run-time library for the large memory model.
\end{enumerate}

\section{Setting up the Zortech compilers for \scAD}

\X{Zortech}
\XX{IDE@\textsc{ide}}{Zortech}
To compile \scAD\ applications with the integrated development environment
\textsc{ide}), the following options must be selected within \textsc{ide}:
\begin{itemize}
  \item Select {\bf C}ompile Com{\bf p}ile Options. Select {\bf M}emory Models
  L{\bf a}rge to select the large memory model. Select Object {\bf C}ode Inline
  808{\bf 7} to generate direct inline 8087 code.

  \item Append the directories where you have stored the \scAD\ header file and
  the \scAD\ libraries to the environment strings \texttt{INCLUDE} and
  \texttt{LIB}, as discussed in the compiler documentation.
\end{itemize}

\XX{compiler options}{Zortech}
\X{memory model}
You should invoke the following compiler and linker options (see the file
\texttt{makefile} supplied with the examples):

\begin{commandtable}
 I  &  followed by the path of your \scAD\ include file.\\
 ml  &  to specify the large memory model.\\
 c  &  to suppress linking after compilation.\\
 g  &  if you want to use the Zortech debugger.\\
 f  &  if you want to generate in-line 8087 code.\\
\end{commandtable}

\section{Examples}

\X{examples}
Many of the examples in this manual are included on the distribution
%disketted fixed
as a ``self-expanding'' archive program \texttt{examples.exe}. When this program
is executed, it expands into several files. The installation program will
install the examples in the directory you specify. A \texttt{makefile} is
included to simplify compilation of the examples. Read the file
\texttt{examples.txt} for a more complete description of the examples.

%xx \htmlnewfile

\chapter{The \scAD\ Classes}

\X{\fontindexentry{tt}{double}}
\X{\fontindexentry{tt}{float}}
\X{\fontindexentry{tt}{dvariable}}
\X{\fontindexentry{tt}{dvector}}
\X{\fontindexentry{tt}{dvar\_vector}}
\X{constant objects}
\X{container classes}
\X{variable objects}
\X{\fontindexentry{tt}{dmatrix}}
\X{\fontindexentry{tt}{dvar\_matrix}}
\X{\fontindexentry{tt}{d3array}}
\X{\fontindexentry{tt}{dvar3\_array}}
\X{\fontindexentry{tt}{long int}}
\X{\fontindexentry{tt}{ivector}}
\X{\fontindexentry{tt}{lvector}}
\X{\fontindexentry{tt}{imatrix}}
\X{\fontindexentry{tt}{lmatrix}}
\X{\fontindexentry{tt}{int}}

The \scAD\ floating point container classes have been designed to facilitate the
calculation of derivatives. As in differential calculus, they reflect the main
dichotomy between constant objects for which no derivative calculations are
required, and variable objects for which derivative information is generated.
Both constant and variable objects can be further classified either as single
objects (a number), or 1-dimensional arrays, 2-dimensional arrays, and so on.
The constant and variable versions of each array type (container class) have
been designed to appear as similar as possible, so that identical code can be
written for calculations involving either constant or variable objects. See
Table~\ref{tab:types}.
\begin{table}[!h]
\begin{center}
\begin{tabular}{@{\vrule height 12pt depth7pt width0pt} l  |  c  |  c }
 \hline
 & \bf Constant Type & \bf Variable Type \\
 \hline
Single element&\texttt{double}, \texttt{float}&\texttt{dvariable} \\
1-dimensional array&\texttt{dvector}&\texttt{dvar\_vector} \\
2-dimensional array&\texttt{dmatrix}&\texttt{dvar\_matrix} \\
3-dimensional array&\texttt{d3array}&\texttt{dvar3\_array} \\[3pt]
\hline
\end{tabular}
\end{center}
\emptycaption
\label{tab:types}
\end{table}

\scAD\ also has 1- and 2-dimensional container classes for integers. See
Table~\ref{tab:integer-container-classes}.
\begin{table}[!h]
\begin{center}
\begin{tabular}{@{\vrule height 12pt depth7pt width0pt}  l  |  l }
\hline
Single element&\texttt{int}, \texttt{long int} \\
1-dimensional array&\texttt{ivector}, \texttt{lvector} \\
2-dimensional array&\texttt{imatrix}, \texttt{lmatrix} \\
\hline
\end{tabular}
\end{center}
\emptycaption
\label{tab:integer-container-classes}
\end{table}

The \texttt{ivector} and \texttt{lvector} classes are arrays of integers of type
\texttt{int} and \texttt{long int}. The \texttt{imatrix} and \texttt{lmatrix}
classes are 2-dimensional arrays of integers of type \texttt{int} and
\texttt{long~int}.

\section{The \texttt{dvariable} class}

Objects of the \texttt{dvariable} class are like numbers (variables). When they
are used in calculations, derivative information is generated. The value of a
\texttt{dvariable} is an object of type \texttt{double}. All the standard C
arithmetic and logical functions have been defined for type \texttt{dvariable}.
A \texttt{dvariable~z} is created by this declaration:
\begin{lstlisting}
  dvariable z;
\end{lstlisting}

Since the \texttt{dvariable} class is designed to appear to the user like a
\texttt{double}, it is not necessary to discuss many of the operations that can
be performed. They appear to be simply the same as those performed on an object
of type \texttt{double}. There is, however, an important difference between the
behavior of a \texttt{dvariable} and that of a \texttt{double} under a copy
operation. A copy occurs, for example, when an object is passed as an argument
to a function (so long as the function has not been declared to take a
reference). Suppose that the function \texttt{func} has been declared
\X{copy operation}
\begin{lstlisting}
  void func(dvariable)
\end{lstlisting}
so that \texttt{func} is a function that takes one argument of type
\texttt{dvariable}, and a portion of the user's code looks like
\begin{lstlisting}
dvariable v;
// ...
func(v);
\end{lstlisting}
\XX{constructor}{copy initializer}
\X{\fontindexentry{tt}{dvariable(dvariable\&)}}

When the function \texttt{func(v)} is called, the copy initializer constructor
will be called first, in order to construct a copy of \texttt{v} and place that
on the stack before passing control to \texttt{func}. For previous versions of
\scAD\ (before Version 1.05), the copy initializer performed a shallow
copy---that is it passed a pointer to the component of the \texttt{dvariable}
containing the \texttt{double}. As a result, the \texttt{dvariable} in the
function \texttt{func} and the \texttt{dvariable v} in the calling routine both
point to the same \texttt{double}, so changes to the value of the
\texttt{dvariable} in \texttt{func} will affect \texttt{v} even though
\texttt{v} was passed by value. To get better compatibility with the behavior of
C code (and the behavior of a \texttt{double}) for Version 1.05 or later, the
\texttt{dvariable} copy initializer constructor has been modified so that it
produces a distinct copy. This could cause a problem for old \scAD\ code. If a
\texttt{dvariable} has been passed by value to a function and its value modified
in the body of the function, it used to be the case that the value of the
\texttt{dvariable} in the calling function would be changed as well. This is no
longer the case. Care must be taken to modify old code where necessary.
\XX{copy operation}{shallow}

%
% \XX{distinct copy}{using assignment operator = to create}
% \XX{distinct copy}{of a dvariable}
% This section applies only to versions of \scAD\ with
% version number less than 1.05.
% For older versions of \scAD\ (version number less than 1.05)
% when it is necessary to create a distinct copy of a \texttt{dvariable}
% the assignment operator~``\texttt{=}'' should be used.
% To create a distinct copy of \texttt{v} whose value is not affected by
% the function \texttt{func} the preceding code could be modified to
% \begin{lstlisting}
% dvariable v;
% dvariable vcopy;
% vcopy=v; // make a distinct copy of v to pass to the function func
% // ...
% func(vcopy);
% \end{lstlisting}
% \noindent
% {\bf It is important to realize that the following simpler code will not
% accomplish the objective of creating a distinct copy of \texttt{v}.}
% \begin{lstlisting}
% dvariable v;        // Not the right way to make a distinct copy
% dvariable vcopy=v;  // This will NOT make a distinct copy
% // ...
% func(vcopy);
% \end{lstlisting}
% \X{initialization}\XX{copy initializer constructor}{for a dvariable}
% \noindent The problem is that while the code
% \texttt{dvariable vcopy=v;} appears to be an assignment, it is in fact
% interpreted by the compiler as an initialization
% and will cause the copy initializer constructor to be invoked.
% Thus a distinct copy of v is not created and the problem persists.
%

Questions about distinct copies and the form of the copy initializer also arise
with container classes, such as the \scAD\ vector and matrix classes. This is
discussed below for the \texttt{dvar\_vector} class.

There is one other operation on \texttt{dvariable}s that has no counterpart for
\texttt{double}s: the
\X{\fontindexentry{tt}{value} function}
\texttt{value} function:
\begin{lstlisting}
  double value(dvariable&);
\end{lstlisting}
\noindent This function returns a \texttt{double}, which is the value of the
\texttt{dvariable}. The main use of this function is to pass to some other part
of the user's program the value of a some computation for which derivative
information has been calculated. That way, the value can be used without
creating a derivative calculation associated either with the passing of the
value or its subsequent use in other parts of the program. It is not necessary
to use \texttt{value}s for stream output of \texttt{dvariable}s, because the
stream operators \texttt{<{}<} and \texttt{>{}>} have been overloaded to accept
these objects. The value function should be used sparingly. In particular, it
must not be used in the middle of calculations for which the derivatives are
desired. The following example will produce an error:
\XX{stream operators}{\fontindexentry{tt}{<{}<}, \texttt{>{}>}}
\XX{operators}{\fontindexentry{tt}{<{}<}, \texttt{>{}>}}
\XX{stream input and output}{of \texttt{dvariables}}
\begin{lstlisting}
dvariable x,y,z;
double a;
// ...             This code is incorrect
y=x*x;
a=value(y);
cout << " y = " << a << "\n";
z=sin(a); // This is where an error is introduced into the derivatives
cout << " z = " << z << "\n"; This will print out the value of z
\end{lstlisting}

The following correct version of the code will produce exactly the same output
and the derivative calculations will be correct:
\begin{lstlisting}
dvariable x,y,z;
double a;
// ...             This code is correct
y=x*x;
cout << " y = " << y << "\n";
z=sin(y); // Use y instead of a at this point
cout << " z = " << z << "\n"; This will print out the value of z
\end{lstlisting}

It is possible to create an array \texttt{u} of 10 \texttt{dvariables} by using
the declaration \XX{array}{of \texttt{dvariables}}
\begin{lstlisting}
  dvariable u[10]; // creates an array of 10 dvariables indexed from 0..9
\end{lstlisting}

However, this is not the recommended construction. Instead, the
\texttt{dvar\_vector} class should be used as in the declaration
\begin{lstlisting}
  dvar\_vector u(0,9); // creates a dvar_vector of size 10 indexed from 0..9
\end{lstlisting}
This construction generates more efficient code and produces a ``safe'' object
with optional bounds checking. The \texttt{dvariable u[10]} construction may be
useful for converting pre-existing C code where the declaration \texttt{double}
or \texttt{float} is modified to \texttt{dvariable}.

\XX{\fontindexentry{tt}{float}}{no corresponding \scAD\ type}
There is at present no \scAD\ type that corresponds to the C~type
\texttt{float}. All \scAD\ variable types contain a \texttt{double}.

\section{The \texttt{dvector} and \texttt{dvar\_vector} classes}

Since the \texttt{dvector} and \texttt{dvar\_vector} classes are intended to
appear identical to the user (with the exception that derivative information is
generated for the \texttt{dvar\_vector} class), we shall discuss the use of
these classes simultaneously.

\X{safe arrays}
\XX{maximum size}{of a \texttt{dvar\_vector}}
The \texttt{dvar\_vector} class appears to the user as a vector (array) of
objects of type \texttt{dvariable}. Bounds checking for array and matrix objects
can be implemented---or not---by linking the user's program with the appropriate
\scAD\ library. The index for a \texttt{dvector} or a \texttt{dvar\_vector} is
of type \texttt{int} (a signed integer), so permissible array bounds are
$-32,768$ to $32,767$ for $16$-bit \textsc{dos} versions, where an object of
type \texttt{int} occupies two bytes. The maximum number of elements that can be
contained in a \texttt{dvar\_vector} is $8191$ in $16$-bit \textsc{dos} versions
of \scAD.

\section{Creating \texttt{dvector}s and \texttt{dvar\_vector}s}

Several declarations can be used to create a \texttt{dvector} or
\texttt{dvar\_vector} object:
\begin{lstlisting}
dvector v(1,n); // Creates a dvector with array bounds 1..n
dvector v=w;    // Creates a dvector and initializes it with w
dvector w("{1.2,4,-5.23e-1}");  // Creates a dvector with array bounds 1..3
                     // and initializes it with 1.2, 4, -5.23e-1.
dvector w("mystuff.dat"); //Creates a dvector whose contents are
                     // read in from the file mystuff.dat. The data in
                     // mystuff.dat must be invalid numeric format
\end{lstlisting}

The code
\begin{lstlisting}
  dvector v=w;
\end{lstlisting}
creates a \texttt{dvector v} by using the ``copy initializer''
\texttt{dvector(dvector\&)} constructor. This constructor performs a shallow
copy, so after this line of code is executed, \texttt{v} and \texttt{w} will
both point to the same area of memory where the array is contained. The same
type of constructor is invoked when a \texttt{dvector} or \texttt{dvar\_vector}
object is passed by value to a function. For example, suppose that a function
\texttt{users\_function} has been defined where the function prototype is
\begin{lstlisting}
  void users_function(dvar_vector v);
\end{lstlisting}
and that the following statement appears somewhere in the user's code:
\begin{lstlisting}
dvar_vector w(1,n);
. . .
users_function(w);
. . .
\end{lstlisting}

To pass the \texttt{dvar\_vector w} to \texttt{users\_function}, a copy of
\texttt{w} is made and put on the stack. The copy initializer constructor is
called by the compiler to do this. The copy initializer uses a ``shallow'' copy
to create the copy of \texttt{w}, which is passed to the function.

A \texttt{dvector} contains a pointer that points to the area of memory where
the elements of the \texttt{dvector} are contained. With a shallow copy, the
value of this pointer is simply passed to the new copy of the \texttt{dvector
  w}. As a result, both \texttt{w} in the calling routine and the copy of
\texttt{w} in the function \texttt{users\_function} point to the same area of
memory. So, changes made to this \texttt{dvector} in \texttt{users\_function}
will change the values of the corresponding entries of~\texttt{w} in the calling
routine, even though the \texttt{dvector} was passed by value. This can create a
problem if the user wants to operate on the \texttt{dvector} in
\texttt{users\_function} without changing the value of \texttt{w} in the calling
routine.

For ordinary C language objects, changes within a function to an argument passed
by value to the function do not affect the object passed. This is not the case
for the \scAD\ container classes \texttt{dvariable}, \texttt{dvector},
\texttt{dvar\_vector}, \texttt{dmatrix}, or \texttt{dvar\_matrix}, or
higher-dimensional arrays. If you wish to pass a ``copy'' of an object to a
function such that the original object is not affected by changes to the passed
object, you must explicitly make a distinct copy of the object by using the
assignment operator (\texttt{=}), as in this example.
\XX{AVOID}{using a shallow copy when a distinct object is desired}
\begin{lstlisting}
dvar_vector v(1,n);  // Want to pass a copy of v to the function func
dvar_vector w(1,n);  // This will be the distinct copy of v
w=v;                 // This creates the distinct copy
users_function(w);   // Can now call the function
\end{lstlisting}

The user must avoid the simpler \cplus\ construction of declaring an initialized
object. Such a declaration will cause the copy initializer to be invoked, so a
shallow copy will be performed.
\begin{lstlisting}
dvar_vector v(1,n);
dvar_vector w=v;     // This uses the copy initializer so that w is
                     // a shallow copy of v
users_function(w);   // Changes to w in the function will cause changes
                     // in v
\end{lstlisting}

It is possible to create a \texttt{dvector} by explicitly naming the data it
will contain: \XX{constructor}{\fontindexentry{tt}{dvector}}
\begin{lstlisting}
  dvector w("{12,24,56 78 91.4 23.455e+2}");
\end{lstlisting}

The data must be enclosed in double quotes and begin and end with braces
(\texttt{\{{}\}}). The numbers must be in valid numeric format and be separated
by blanks or commas. It is also possible to create a \texttt{dvector} using the
contents of a file. The declaration
\begin{lstlisting}
  dvector w("mystuff.dat");
\end{lstlisting}
will create a \texttt{dvector} \texttt{w} that contains the contents of the file
\texttt{mystuff.dat}. It is not necessary to know how many numbers are contained
in the file to use this declaration. An error message will be generated if the
file does not exist, the file is empty, or there is non-numeric data in it. The
data in the file must be separated by spaces. This constructor will create a
\texttt{dvector} whose minimum valid index is~1.

It is also possible to create a \texttt{dvector} or a \texttt{dvar\_vector} by
reading in a column of data from a file. Suppose the file \texttt{``work.prn''}
contains the data
\begin{daysdecimaltable}
  Monday     & 25&.12 & -41&       & 433&.12 & 10&.11 & 23& \\
  Tuesday    & 34&.21 & 356&       &  23&    & 23&.1  &  5& \\
  Wednesday  & 10&    &   3&.13e-4 &  43&.11 &  3&.23 &  4& \\
  Thursday   & 12&.1  &  53&.13    & 453&.1  & -5&.13 &  4&.1 \\
\end{daysdecimaltable}
The declaration
\begin{lstlisting}
  dvar_vector u("work.prn",3);
\end{lstlisting}
will create a \texttt{dvar\_vector} whose elements are the third column of the
file \texttt{work.prn}. That~is,
\begin{lstlisting}
  u = (\ts-41,\ts356,\ts3.13e-4,\ts53.13\ts)
\end{lstlisting}
The entries of the file must be separated by blanks and all the entries in the
column being read must be in valid numeric format. This constructor will create
a \texttt{dvector} whose minimum valid index is~1.

% \section{Avoiding arithmetic statements in the initialization of
% \texttt{dvariable} objects}
% {\bf This section only applies to version of \scAD\ with
% version number less than 1.05}
%
% Including statements like
% \begin{lstlisting}
% dvariable x=a*b;
% \end{lstlisting}
% \noindent in your code can introduce errors.
% This is due to the use of shallow copies in \scAD\ combined with the
% the method used to allocate memory for objects the results containing
% %use of a circular queue of temporary objects  for containing
% intermediate calculations. If \texttt{a} and \texttt{b} are
% \texttt{dvariables} then the operation \texttt{a*b} is stored in a temporary
% variable. The shallow copy which is then used to initialize
% the variable \texttt{x} will leave \texttt{x} pointing to the area
% of memory allocated to the temporary variable. Later on this
% memory location may be overwritten leading to incorrect
% results in the program. To avoid this type of error the proper
% construction is to split up this statement into two statements.
% \begin{lstlisting}
% dvariable x;
% x=a*b;
% \end{lstlisting}
% \XX{AVOID} {Initializing a \texttt{dvariable} with an arithmetic expression}
% \noindent The statement \texttt{x=a*b} performs a real copy so that
% \texttt{x} does not share the same memory location as the temporary
% variable.
% This restriction applies only to objects of type \texttt{dvariable}.
% It is perfectly all right to use arithmetic statements in
% the initialization of variable array objects. For example
% \begin{lstlisting}
% dvar_matrix M1=M*trans(M);
% \end{lstlisting}
% \noindent where \texttt{M} is a \texttt{dvar\_matrix} or
% \begin{lstlisting}
% dvar_vector v= 10*u;
% \end{lstlisting}
% \noindent where \texttt{u} is a \texttt{dvar\_vector}.
% resizable

Since this version of \scAD\ does not feature resizable arrays, the only way
that you can place the product of two \texttt{dvar\_matrices}~\texttt{M1}
and~\texttt{M2} in a \texttt{dvar\_matrix}~\texttt{M3} without knowing or
calculating the size of the product is to use the construction
\begin{lstlisting}
  dvar_matrix M3=M1*M2;
\end{lstlisting}
so the copy initializer constructor will create an object with the correct size.
Of course, it is possible to calculate the size of the product and explicitly
declare the \texttt{dvar\_vector~M3} by using the \texttt{dvar\_matrix} access
functions to obtain the row and column dimensions of the matrices being
multiplied, as in the code
\begin{lstlisting}
dvar\_matrix M3(M1.rowmin(),M1.rowmax(),M2.colmin(),M2.colmax());
M3=M1*M2;
\end{lstlisting}
but this solution lacks elegance.

\section{The \texttt{dvector} and \texttt{dvar\_vector} access functions}

Most of the components of the \scAD\ classes are private, so they cannot be
accessed directly by the user. To allow the user to access the relevant elements
of the classes, public class member access functions are provided:
\begin{lstlisting}
int indexmin();     //Returns the minimum allowable index  for an array
int indexmax();     //Returns the maximum allowable index  for an array
int size();         //Returns the size of an array
int shift(int min); //Changes the allowable indices to min and min+size-1
\end{lstlisting}

These functions can be used to determine the size and valid index ranges of a
vector object that has been passed to a user's function. Thus, it is not
necessary to explicitly pass such information when designing a function that
uses a vector object. The information is already passed with the object itself.
If \texttt{v} is a \texttt{dvector}, the function \texttt{v.indexmin()} will
return the minimum valid index for \texttt{v}. The \texttt{shift} function
allows the valid index range for a \texttt{dvector} or \texttt{dvar\_vector} to
be changed. As an example of a function that acts on a vector objects, here is a
listing for a function that overloads the multiplication operator (\texttt{*});
\texttt{v * w} will be the dot product (also sometimes called the ``scalar
product'') of the \texttt{dvectors}~\texttt{v} and~\texttt{w}:
\begin{lstlisting}
double operator * (dvector& t1, dvector& t2)
{
  // Check to see if vectors have the same valid index bounds
  if (t1.indexmin() != t2.indexmin() ||  t1.indexmax() != t2.indexmax())
  {
    cerr << "Index bounds do not match in double operator * "
            "(dvector&,dvector&)\n";
    exit(1);  // Stop if vectors are not compatible
  }
  double tmp;
  tmp=0;
  for (int i=t1.indexmin(); i<=t1.indexmax(); i++)
  {
    tmp+=t1[i]*t2[i];
  }
  return(tmp);
}
\end{lstlisting}
\XX{AVOID}{only creating a version of a function that takes a variable argument}

If you write your own functions, it is a good idea to provide both a constant
and variable version at the time of writing. This is especially true if you
write a variable version of the routine. Suppose that the user has written a
function \texttt{void userfun(dvar\_matrix)} that takes a \texttt{dvar\_matrix}
argument, but has neglected to make a version \texttt{void userfun(dmatrix)}
that takes a \texttt{dmatrix} argument. At some later time, the user wants to
use the function \texttt{userfun} with an argument of type \texttt{dmatrix},
such as in the following code segment.
\begin{lstlisting}
{
  dmatrix M(1,5,1,5);
  userfun(M); // M will be converted to a dvar_matrix
  // ...
}
\end{lstlisting}
Since there is no version of \texttt{userfun} that accepts an argument of type
\texttt{dmatrix}, the \cplus\ compiler will convert \texttt{M} to an object of
type \texttt{dvar\_matrix}, by invoking the constructor for a
\texttt{dvar\_matrix(dmatrix\&)}. This is probably not what the user intends and
will be disastrous if an object of type \texttt{gradient\_structure} is not in
scope at the time.

\XX{dot product}{of vectors}
The corresponding variable version of the dot product for \texttt{dvar\_vectors}
would look like
\begin{lstlisting}
dvariable operator * (dvar_vector& t1, dvar_vector& t2)
{
  if (t1.indexmin() != t2.indexmin() ||  t1.indexmax() != t2.indexmax())
  {
    cerr << "Index bounds do not match in dvariable operator * "
            "(dvar_vector&,dvar_vector&)\n";
    exit(1);
  }
  RETURN_ARRAYS_INCREMENT();
  dvariable tmp;
  tmp=0;

  for (int i=t1.indexmin(); i<=t1.indexmax(); i++)
  {
    tmp+=t1[i]*t2[i];
  }
  RETURN_ARRAYS_DECREMENT();
  return(tmp);
}
\end{lstlisting}
\X{\fontindexentry{tt}{RETURN\_ARRAYS\_INCREMENT}}
\X{\fontindexentry{tt}{RETURN\_ARRAYS\_DECREMENT}}

The appearance of the statements
\begin{lstlisting}
  RETURN_ARRAYS_INCREMENT();
\end{lstlisting}
and
\begin{lstlisting}
  RETURN_ARRAYS_DECREMENT();
\end{lstlisting}
requires some explanation. \scAD\ employs a stack of queues to hold intermediate
arithmetic results for variable objects. To ensure that these queues are not
overflowed, you should always increment the queue structure when entering a
function that returns a variable object, and decrement the queue structure when
leaving the function. If you fail to increment the queue structure and a queue
overflow occurs, the program will continue to operate, but the arithmetic will
be corrupted. If you fail to decrement the queue structure when leaving the
function, you may eventually increment past the end of the stack of queues,
producing an error message
\begin{lstlisting}
Overflow in RETURN_ARRAYS stack---Increase NUM_RETURN_ARRAYS
There may be a RETURN_ARRAYS_INCREMENT()
which is not matched by a RETURN_ARRAYS_DECREMENT()
\end{lstlisting}

Simply include the increment and decrement statements in every function you
create that returns a variable object and everything will be fine. It is not
necessary to include the increment and decrement statements in a function that
does not return a variable type, because such a function cannot be used in an
arithmetic statement. \textit{You should never include the}
\texttt{RETURN\_ARRAYS\_INCREMENT()} \textit{and}
\texttt{RETURN\_ARRAYS\_DECREMENT()}
\textit{functions in a function that returns a void or constant type, such as}
\texttt{double}, \texttt{dvector}, \textit{or} \texttt{dmatrix}. The reason is
that such a function might be called when there is no object from the class
\texttt{gradient\_structure} in scope. Calling
\texttt{RETURN\_ARRAYS\_INCREMENT()} or \texttt{RETURN\_ARRAYS\_DE\-CRE\-MENT()}
when there is no object of type \texttt{gradient\_structure} in scope could
cause a serious program error.
\XX{AVOID}{calling \texttt{RETURN\_ARRAYS\_INCREMENT()}}
\XX{AVOID}{calling \texttt{RETURN\_ARRAYS\_DECREMENT()}}

\section{Accessing array elements}
\X{operators \texttt{()} and \texttt{[]}}
\XX{access functions}{array elements}

Both the \texttt{()} and the \texttt{[]} operators can be used to access
elements of arrays. There is absolutely no difference between them. If
\texttt{v} is a \texttt{dvector}, the $\texttt{i}^\textrm{th}$ element of
\texttt{v} can be accessed by writing either \texttt{v[i]} or \texttt{v(i)}.
Optional array bounds checking is provided for both operators. Array bounds
checking is implemented or turned off by linking with the appropriate library.
No changes in the user's code are required. Linking with the safe libraries
enables array bounds checking, while linking with the optimized libraries
disables array bounds checking.
\XX{libraries}{safe}
\XX{libraries}{optimized}

\section{Creating column vectors and row vectors}

It is important to recognize that there is a difference in \scAD\ between a
vector object and a matrix object that has either 1 row or 1 column. This means
that a vector object is neither a row vector nor a column vector, it is simply a
vector object. We have adopted this treatment of vector objects so that it is
possible to multiply vector objects such as \texttt{u} and \texttt{v} simply by
writing \texttt{u*v}. If vector objects were the same as matrices with 1~column,
it would be necessary to write \texttt{trans(u)*v} to multiply two vectors. Here
\texttt{trans} denotes the transpose of a \texttt{matrix} object.

\X{\fontindexentry{tt}{column\_vector function}}
\X{\fontindexentry{tt}{row\_vector} function}
\XX{constructor}{\fontindexentry{tt}{dmatrix}}
Sometimes it is useful to be able to treat vector objects as though they were
matrices with one row or column. The functions \texttt{column\_vector} and
\texttt{row\_vector} enable this to be done. The function
\begin{lstlisting}
  dmatrix column_vector(dvector& v);
\end{lstlisting}
will convert \texttt{v} into a matrix with 1 column. The valid column index is
from~1 to~1. The valid row index bounds are the same as the valid index bounds
for \texttt{v}. The function
\begin{lstlisting}
  dmatrix row_vector(dvector& v);
\end{lstlisting}
will convert \texttt{v} into a matrix with 1 row. The valid row index is from~1
to~1. The valid column index bounds are the same as the valid index bounds for
\texttt{v}.

As an example of how to use these functions, consider the problem of computing
the outer product of two vectors. The outer product of two vectors $(u_i)$ and
$(v_j)$ is the matrix $w_{ij}$, where $w_{ij}=u_iv_j$. The statement
\begin{lstlisting}
  dmatrix w=column_matrix(v)*row_matrix(u);
\end{lstlisting}
will compute the outer product of the \texttt{dvectors} \texttt{u} and
\texttt{v}. Of course, you can also use the function \texttt{outer\_prod},
supplied with the \scAD\ libraries, to accomplish the same thing:
\begin{lstlisting}
  dmatrix w=outer_prod(u,v);
\end{lstlisting}

\XX{operators}{\fontindexentry{tt}{<{}<}}
Another use for the \texttt{column\_vector} function is to print out a vector
one element to a line. The overloaded operator \texttt{<{}<} will write out a
vector object all on one line. This is useful if, for example, you want to put a
vector object into a row in a spreadsheet, but it can be a problem for some
editors that cannot read such long lines in a file. The code
\begin{lstlisting}
  cout << column_vector(v);
\end{lstlisting}
will print out the vector \texttt{v} with one entry per line.

\section{The \texttt{dmatrix} and \texttt{dvar\_matrix} classes}

The \texttt{dmatrix} and \texttt{dvar\_matrix} classes are implemented,
respectively, as arrays of \texttt{dvector} and \texttt{dvar\_vectors}. If
\texttt{M} is a \texttt{dmatrix}, then \texttt{M[j]} or \texttt{M(j)} is a
\texttt{dvector} whose elements are the $\texttt{j}^\textrm{th}$ row of
\texttt{M}. Since each row is a \texttt{dvector} that contains its own size
(shape) information, it is a simple matter to construct \texttt{dmatrix}
objects---so-called ``ragged'' matrices---that have a different numbers of rows
in each column. Of course, these are not matrices in the usual sense in which
the word is used, but they can be useful container objects for some
applications. Therefore, we have extended the \texttt{dmatrix} and
\texttt{dvar\_matrix} classes to include them. As with the \texttt{dvector and
  dvar\_vector} classes, \texttt{dmatrix and dvar\_matrix} classes appear
identical to the user, with the exception that derivative information is
collected for the \texttt{dvar\_matrix class}. An object of type
\texttt{dvar\_matrix} is created by using the declarations
\begin{lstlisting}
dmatrix M(lbr,lur,lbc,luc); // Creates a (lur-lbr+1)$\times$(luc-lbc+1) matrix
dmatrix M1=M;               // M has already been defined. M1 is initialized
                            // by M.
dmatrix M(lbr,lur,ivec_lb,ivec_ub); // Creates a ``ragged'' dmatrix with a
                            // variable number of rows in each column
dmatrix M("filename");      // Creates a dmatrix whose elements are the
                            // contents of the file ``filename''
dmatrix M( "{ 1.5,-2.0,3.0,3.5 }"
         "{ 2.0,14.0,5.1,2.2 }" );// creates and initializes a dmatrix
\end{lstlisting}

The constructor
\begin{lstlisting}
  dmatrix(lbr,lur,ivec_lb,ivec_ub)
\end{lstlisting}
creates a \texttt{dmatrix} with \texttt{lur-lbr+1} rows. The
$\texttt{i}^\textrm{th}$ row is a \texttt{dvector} of size $\texttt{ivec\_ub[i]}
- \texttt{ivec\_lb[i]+1}$ with a range of valid indices running from
\texttt{ivec\_lb[i]} to \texttt{ivec\_ub[i]}. \texttt{ivec\_lb} and
\texttt{ivec\_ub} are instances of the class \texttt{ivector}, that is, vectors
of integers. They must have valid indices ranging from \texttt{lbr} to
\texttt{lur}. For example, suppose we wish to make a ragged
\texttt{dvar\_matrix} with 4 rows for which the valid row index bounds should be
from~2 to~5. The declaration should be
\begin{lstlisting}
  dvar_matrix M(2,5,ivec_lb,ivec_ub);
\end{lstlisting}

Suppose the desired valid row indices are 1 to 6 for row 1, 0 to 3 for row 2, 1
to 2 for row~4, and~3 to~9 for row~5. The legal bounds on the \texttt{ivector}s
\texttt{ivec\_lb} and \texttt{ivec\_ub} must be equal to the legal column bounds
on \texttt{M}. We can create such \texttt{ivector} objects with the declaration
\begin{lstlisting}
ivector ivec_lb(2,5);
ivector ivec_ub(2,5);
\end{lstlisting}

The correct values for \texttt{ivec\_lb} and \texttt{ivec\_ub} can be inserted
for example by using the ``fill'' function
\begin{lstlisting}
ivec_lb.fill("{1,0,1,3}");
ivec_ub.fill("{6,3,2,9}");
\end{lstlisting}

Alternatively, the desired \texttt{ivector} objects can be created with the
declarations
\begin{lstlisting}
ivector ivec_lb("{1,0,1,3}");
ivector ivec_ub("{6,3,2,9}");
\end{lstlisting}

These declarations will create \texttt{ivector} objects whose legal bounds go
from~1 to~4. To change the valid index bounds to the desired values from~2 to~5,
the class member function \texttt{shift()} can be used:
\begin{lstlisting}
ivec_lb.shift(2); // legal bounds will be from 2 to 5
ivec_ub.shift(2); // legal bounds will be from 2 to 5
\end{lstlisting}

As is the case for vector objects, the copy initializer for matrix objects
performs a light copy. To create a distinct new matrix object, you must use the
assignment operator (\texttt{=}).

The declaration
\begin{lstlisting}
  dmatrix M("filename");
\end{lstlisting}
will create a \texttt{dmatrix} \texttt{M} whose elements are the contents of the
file \texttt{filename}. All entries in the file must be in valid numeric format.
All entries must be separated by blanks or commas. Each line in the file becomes
a line of the matrix. Ragged matrices are supported, so the files may have
different number of entries in each line. For the \textsc{dos} version of \scAD,
no line may contain more than 6550 elements.

To create a \texttt{dmatrix} and initialize it with data, a declaration of the
form
\begin{lstlisting}
  dmatrix W( "{ 1.5,-2.0,3.0,3.5 }"
           "{ 2.0,14.0,5.1,2.2 }" );
\end{lstlisting}
can be used.

This declaration will create a dmatrix with~2~rows and~4~columns. The minimum
valid index for the rows and columns is set equal to~1.

\section{The \texttt{dmatrix} access functions}

\XX{\fontindexentry{tt}{dmatrix}}{access functions}
\XX{access functions}{for \texttt{dmatrix} and \texttt{dvar\_matrix} class}
The components of all the \scAD\ classes are private, which means that they
cannot be accessed directly by the user. Public class member access functions
are provided to allow the user to access the relevant elements of the classes:
\begin{lstlisting}
int rowmin();   //Returns the minimum allowable row index  for a matrix
int rowmax();   //Returns the maximum allowable row index  for a matrix
int colmin();   //Returns the minimum allowable column index  for a matrix
int colmax();   //Returns the maximum allowable column index  for a matrix
int rowsize();  //Returns the number of rows in an array
int colsize();  //Returns the number of columns in an array
void colshift(int&);  //Changes the range of valid indices for the columns
void rowshift(int&);  //Changes the range of valid indices for the rows
\end{lstlisting}
\X{\fontindexentry{tt}{rowmax()}}
\X{\fontindexentry{tt}{colmin()}}
\X{\fontindexentry{tt}{colmax()}}
\X{\fontindexentry{tt}{rowsize()}}
\X{\fontindexentry{tt}{colsize()}}
\X{\fontindexentry{tt}{colshift(int\&)}}
\X{\fontindexentry{tt}{rowshift(int\&)}}

\section{The 3-dimensional arrays\br \texttt{d3array} and \texttt{dvar3\_array}}
\X{\fontindexentry{tt}{d3array} and \texttt{dvar3\_array} classes}

Three-dimensional arrays have been implemented in this version of \scAD\ mainly
as containers for collections of 2-dimensional arrays. If \texttt{u} is a
\texttt{d3array} then \texttt{u[1]} or \texttt{u(1)} is a dmatrix (provided
that~1 is a valid index for \texttt{u}). We shall refer to the object obtained
by fixing the first index of a 3-dimensional array as a ``slice.'' An element of
a 3-dimensional array is determined by picking a slice, a row, and a column.
\XX{constructor}{\fontindexentry{tt}{d3array}} \XX{slice}{of a \texttt{d3array}}

An object of type \texttt{dvar3\_array} is created by using the declarations
\begin{lstlisting}
dvar_3darray M(ls,us,lr,ur,lc,uc); // Creates a (us-ls+1)*(ur-lr+1)*(uc-lc+1)
                                   // dvar_3darray
dvar_3darray M1=M;                 // M has already been defined. M1 is
                                   // initialized by M.
d3array(ls,us,ivec_lr,ivec_ur,ivec_lc,ivec_uc); // Creates a ragged d3array
                            // with a variable size dmatrices in each slice.
\end{lstlisting}

\section{The \texttt{d3array} access functions}
\XX{access functions}{for a \texttt{d3array}}
\XX{\fontindexentry{tt}{d3array}}{access functions}

\begin{lstlisting}
int slicemin();   //Returns the minimum allowable slice index  for a d3array
int slicemax();   //Returns the maximum allowable slice index  for a d3array
int rowmin();   //Returns the minimum allowable row index  for a d3array
int rowmax();   //Returns the maximum allowable row index  for a d3array
int colmin();   //Returns the minimum allowable column index  for a d3array
int colmax();   //Returns the maximum allowable column index  for a d3array
int slicesize();  //Returns the number of slices in a d3array
int rowsize();  //Returns the number of rows in an d3array
int colsize();  //Returns the number of columns in an array
void sliceshift(int&);  //Changes the range of valid indices for the slices
void colshift(int&);  //Changes the range of valid indices for the columns
void rowshift(int&);  //Changes the range of valid indices for the rows
\end{lstlisting}
\X{\fontindexentry{tt}{slicemin()}}
\X{\fontindexentry{tt}{slicemax()}}
\X{\fontindexentry{tt}{slicesize()}}
\X{\fontindexentry{tt}{rowmin()}}
\X{\fontindexentry{tt}{rowmax()}}
\X{\fontindexentry{tt}{colmin()}}
\X{\fontindexentry{tt}{colmax()}}
\X{\fontindexentry{tt}{rowsize()}}
\X{\fontindexentry{tt}{colsize()}}
\X{\fontindexentry{tt}{colshift(int\&)}}
\X{\fontindexentry{tt}{rowshift(int\&)}}
\X{\fontindexentry{tt}{slicehift(int\&)}}
\XX{conversion}{between between integer and floating point objects}

\section{Converting between the integer\br and floating point container objects}

It is often useful to be able to convert a floating point object into an integer
object, or an integer object into a floating point object. (See, for instance,
the example on creating a bootstrap sample, where a \texttt{dvector} object of
random numbers is converted into an \texttt{ivector} object of integer indices.)
It is possible to convert between vector objects as follows:
\begin{lstlisting}
ivector <--> dvector
lvector <--> dvector
lvector <--> ivector
\end{lstlisting}

These conversions may be made explicit, as in the following code fragment:
\begin{lstlisting}
dvector u(1,10);
ivector iv(1,10);
// ...
iv = ivector(u); // iv will be filled with the integer parts of the
                 // components of u
\end{lstlisting}
or the conversion may be ``implicit,'' that is, supplied by the compiler.
\X{implicit conversion}
\XX{conversion}{implicit}
\begin{lstlisting}
dvector u(1,10);
// ...
ivector iv = u; // iv will be filled with the integer parts of the
                // components of u
\end{lstlisting}

In the second example, the \texttt{ivector} \texttt{iv} will have minimum and
maximum valid indices determined by the minimum and maximum valid indices of the
\texttt{dvector}~\texttt{u}.

%xx \htmlnewfile

\chapter{Operations and Functions}

\section{The operators  $+$ and  $-$}
\X{operators $+$ and $-$}
\XX{vector}{component-wise sum of}
\X{addition of vectors}
\XX{optimizing performance}{using the best operators for a calculation}

Since \scAD\ supplies a large number of operators and functions, there is often
more than one way to carry out a given calculation. While all equivalent methods
will produce the same answer, one method may produce more efficient code. For
example, the two lines of~code
\begin{lstlisting}
A=A+B;
//
A+=B;
\end{lstlisting}
will both add the object \texttt{A} to the object \texttt{B} and store the
result in \texttt{A}. (This will work whether \texttt{A} and \texttt{B} are
numbers, vectors, or matrices.) However, the code produced by \texttt{A+=B} is
more efficient in all cases. It is both faster and generates less temporary data
to be stored.

While the type of object produced by an arithmetic operation depends on whether
the arguments are constants (\texttt{double}, \texttt{dvector} or
\texttt{dmatrix}) or variables (\texttt{dvariable}, \texttt{dvar\_vector}, or
\texttt{dvar\_matrix}), the sort of operation defined by the operation does not.
The symbol~\texttt{+} between two vector objects will denote the component-wise
sum of these vectors and produce a vector object corresponding to their sum.
This is true regardless of whether one or both objects are constant or variable
objects. We denote this by writing
\begin{lstlisting}
  vector_object = vector_object + vector_object // vector sum
\end{lstlisting}

We shall refer to the ``output'' object on the left side of the equals sign
($=$) by $z$ if it is a number, $z_i$ if it is a vector and $z_{ij}$ if it is a
matrix---in order to express the components of the operation. We shall refer to
the first input object (the one on the right side of the equals sign) by $x$,
$x_i$, and $x_{ij}$, and the second input object---if there is one---by $y$,
$y_i$, and $y_{ij}$. For any operation, the property of being a variable is
dominant in that if either of the input objects is a variable, then the output
object is a variable.
\begin{lstlisting}
vector_object = vector_object + vector_object // vector sum
\end{lstlisting}
\afterlistingthing{$z_i=x_i+y_i$}
\begin{lstlisting}
vector_object = number + vector_object // add number to vector
\end{lstlisting}
\afterlistingthing{$z_i=x+y_i$}
\begin{lstlisting}
vector_object = vector_object + number // add number to vector
\end{lstlisting}
\afterlistingthing{$z_i=x_i+y$}
\begin{lstlisting}
matrix_object = matrix_object + matrix_object // matrix sum
\end{lstlisting}
\afterlistingthing{$z_{ij}=x_{ij}+y_{ij}$}
\begin{lstlisting}
matrix_object = matrix_object + number // matrix sum
\end{lstlisting}
\afterlistingthing{$z_{ij}=x_{ij}+y$}
\XX{vector}{component-wise difference of}\X{subtraction of vectors}
\begin{lstlisting}
vector_object = vector_object - vector_object // vector difference
\end{lstlisting}
\afterlistingthing{$z_i=z_i-y_i$}
\XX{vector}{subtracting a number from a}
\begin{lstlisting}
vector_object = vector_object - number // subtract number from a vector
\end{lstlisting}
\afterlistingthing{$z_i=x_i-y$}
\begin{lstlisting}
matrix_object = matrix_object - number // matrix difference
\end{lstlisting}
\afterlistingthing{$z_{ij}=x_{ij}-y$}
\XX{matrix}{subtracting a number from a}
\begin{lstlisting}
matrix_object = matrix_object - matrix_object // matrix difference
\end{lstlisting}
\afterlistingthing{$z_{ij}=x_{ij}-y_{ij}$}

\section{The operators \ttpluseq,  \ttminuseq, \ttmultiplyeq, and \ttdivideeq}
\X{operators \ttpluseq, \ttminuseq, \ttmultiplyeq, and \ttdivideeq}

These operators generate faster code than the corresponding operators defined
above, since their use avoids both the necessity for creating a temporary object
to hold the result, and the extra overhead associated with assigning the values
contained in this temporary object to the left-hand side of the expression.
\begin{lstlisting}
vector_object += vector_object // vector sum
\end{lstlisting}
\afterlistingthing{$x_i \pluseq y_i$}
\begin{lstlisting}
vector_object += number // add number to vector
\end{lstlisting}
\afterlistingthing{$x_i \pluseq d$}
\begin{lstlisting}
vector_object -= vector_object // vector difference
\end{lstlisting}
\afterlistingthing{$x_i \minuseq y_i$}
\begin{lstlisting}
vector_object -= number // subtract number from vector
\end{lstlisting}
\afterlistingthing{$x_i \minuseq d$}
\begin{lstlisting}
vector_object /= number // divide vector by number
\end{lstlisting}
\afterlistingthing{$x_i \divideeq d$}
\begin{lstlisting}
vector_object *= number // multiply vector by number
\end{lstlisting}
\afterlistingthing{$x_i \multiplyeq d$}
\begin{lstlisting}
matrix_object += matrix_object // matrix sum
\end{lstlisting}
\afterlistingthing{$x_{ij} \pluseq y_{ij}$}
\begin{lstlisting}
matrix_object += number // add number to matrix
\end{lstlisting}
\afterlistingthing{$x_{ij} \pluseq d$}
\begin{lstlisting}
matrix_object -= matrix_object // matrix subtraction
\end{lstlisting}
\afterlistingthing{$x_{ij} \minuseq y_{ij}$}
\begin{lstlisting}
matrix_object -= number // subtract number from matrix
\end{lstlisting}
\afterlistingthing{$x_{ij} \minuseq d$}

\section{The operators $*$ and  $/$}
\X{$/$}
\X{operators $*$ and $/$}

\XX{vector}{dot product of two}
\X{dot product}
\begin{lstlisting}
number = vector_object * vector_object // vector dot product
\end{lstlisting}
\afterlistingthing{$z=\sum_i x_iy_i$}
\XX{vector}{multiplying a number by a}\X{scalar product}
\begin{lstlisting}
vector_object = number * vector_object // scalar product
\end{lstlisting}
\afterlistingthing{$z_i=x*y_i$}
\begin{lstlisting}
vector_object = vector_object * number // scalar product
\end{lstlisting}
\afterlistingthing{$z_i=x_i*y$}
\XX{vector}{multiplying a matrix by a}
\XX{matrix}{multiplying a vector times a}
\begin{lstlisting}
vector_object = vector_object * matrix_object // vector times matrix
\end{lstlisting}
\afterlistingthing{$z_j=\sum_i x_i*y_{ij}$}
\XX{vector}{multiplying a matrix times a}
\XX{matrix}{multiplying a vector by a}
\begin{lstlisting}
vector_object = matrix_object * vector_object // matrix times vector
\end{lstlisting}
\afterlistingthing{$z_i=\sum_j x_{ij}*y_{j}$}
\XX{matrix}{multiplication of two}
\X{matrix multiplication}
\begin{lstlisting}
matrix_object = matrix_object * matrix_object // matrix times matrix
\end{lstlisting}
\afterlistingthing{$z_{ij}=\sum_k x_{ik}*y_{kj}$}
\XX{matrix}{multiplication by a number}
\begin{lstlisting}
matrix_object = number * matrix_object // scalar times matrix
\end{lstlisting}
\afterlistingthing{$z_{ij}= x*y_{ij}$}
An additional vector-matrix operation is the outer product of two vectors.
\XX{vector}{outer product of two}
\XX{outer product}{of two vectors}
\begin{lstlisting}
matrix_object = outer_prod(vector_object,vector_object)
  // outer product of two vectors
\end{lstlisting}
\afterlistingthing{$z_{ij}= x_i*y_j$}

\XX{vector}{dividing by a number}
\begin{lstlisting}
vector_object = vector_object / number // divide a vector by a number
\end{lstlisting}
\afterlistingthing{$z_i=x_i/y$}
\XX{vector}{dividing  a number by}
\begin{lstlisting}
vector_object = number / vector_object // divide a number by a vector elements
\end{lstlisting}
\afterlistingthing{$z_i=x/y_i$}
\XX{matrix}{dividing by a number}
\begin{lstlisting}
matrix_object = matrix_object / number // divide a matrix by a number
\end{lstlisting}
\afterlistingthing{$z_{ij}=x_{ij}/y$}

\section{The concatenation of vector objects}

\XX{concatenation}{of vector objects} \XX{vector objects}{concatenation}
The \texttt{\&} operator has been overloaded to produce the concatenation of two
vector objects:
\begin{lstlisting}
  vector_object=vector_object & vector_object
\end{lstlisting}
If \texttt{x}=$(x_r,\ldots,x_n)$ and \texttt{y}=$(y_s,\ldots,y_m)$, then
\texttt{x \& y}= $(x_r,\ldots,x_n,y_s,\ldots,y_m)$ and the minimum valid index
for \texttt{x \& y} is equal to the minimum valid index for~\texttt{x}.

Since in the present version of \scAD\ arrays are not resizable, it is not
possible to write something like \X{resizable arrays}
\begin{lstlisting}
  x = x & y;
\end{lstlisting}
where \texttt{x} and \texttt{y} are \texttt{dvectors}. Instead, you must do
something like
\begin{lstlisting}
  dvector z = x & y;
\end{lstlisting}

Future versions of \scAD\ will incorporate resizable arrays.

\section{Element-wise operations}

There are several operations familiar to users of spreadsheets that do not
appear as often in classical mathematical calculations. For example, spreadsheet
users often wish to multiply one column in a spreadsheet by the corresponding
elements of another column. Spread sheet users might find it much more natural
to define the product of matrices as an element-wise operation such as
\begin{equation*}
  z_{ij}= x_{ij}*y_{ij}
\end{equation*}
The ``classical'' mathematical definition for the matrix product has been
assigned to the overloaded multiplication operator (\texttt{*}), so large
mathematical formulas involving vector and matrix operations can be written in a
concise notation. Typically, spreadsheet-type calculations are not so
complicated and do not suffer so much from being forced to adopt a
``function-style'' of notation.

Since addition and subtraction are already defined in an element-wise manner, it
is only necessary to define element-wise operations for multiplication and
division. We have named these functions \texttt{elem\_prod} and
\texttt{elem\_div}.

\XX{vector}{element-wise product of}
\XX{\fontindexentry{tt}{elem\_prod}}{element-wise product}
\begin{lstlisting}
vector_object = elem_prod(vector_object,vector_object) // element-wise multiply
\end{lstlisting}
\afterlistingthing{$z_i= x_i*y_i$}

\XX{vector}{element-wise division}
\XX{\fontindexentry{tt}{elem\_div}}{element-wise division}
\begin{lstlisting}
vector_object = elem_div(vector_object,vector_object) // element-wise divide
\end{lstlisting}
\afterlistingthing{$z_i= x_i/y_i$}

\XX{matrix}{element-wise product of}
\begin{lstlisting}
matrix_object = elem_prod(matrix_object,matrix_object) // element-wise multiply
\end{lstlisting}
\afterlistingthing{$z_{ij}= x_{ij}*y_{ij}$}

\XX{matrix}{element-wise division}
\begin{lstlisting}
matrix_object = elem_div(matrix_object,matrix_object) // element-wise divide
\end{lstlisting}
\afterlistingthing{$z_{ij}= x_{ij}/y_{ij}$}

\section{The identity matrix function \texttt{identity\_matrix}}
\XX{matrix}{the identity matrix function}

\begin{lstlisting}
matrix_object = identity_matrix(int min,int max)
\end{lstlisting}
\afterlistingthing{creates a square identity matrix with minimum valid index
  \texttt{min} and maximum valid index~\texttt{max}.}

\section{The operations \texttt{det}, \texttt{inv}, \texttt{norm},
  \texttt{norm2}, \br \texttt{min}, \texttt{max}, \texttt{sum}, \texttt{rowsum},
  and \texttt{colsum}}

\subsubsection{The determinant of a matrix object}
\XX{matrix}{the determinant of}
\X{determinant}
\X{det}

(The matrix must be square: that is, the number of rows must equal the number of
columns.)
\begin{lstlisting}
  matrix_object = det(matrix_object)
\end{lstlisting}

\subsubsection{The inverse of a matrix object}
(The matrix must be square: that is, the number of rows must equal the number of
columns.)
\XX{matrix}{the inverse of}
\X{inverse}
\X{inv@\texttt{inv}}
\begin{lstlisting}
  matrix_object = inv(matrix_object)
\end{lstlisting}

\subsubsection{The norm of a vector object}
\XX{vector}{the norm of}
\XX{vector}{the norm squared of}
\XX{matrix}{the norm of}
\X{\fontindexentry{tt}{norm}}
\begin{lstlisting}
number = norm(vector_object)
\end{lstlisting}
\afterlistingthing{$z=\sqrt{\sum_{i} x_{i}^2} $}
\bigskip

\subsubsection{The norm squared of a vector\_object}
\XX{matrix}{the norm squared of}\X{\fontindexentry{tt}{norm2}}
\begin{lstlisting}
number = norm2(vector_object)
\end{lstlisting}
\afterlistingthing{$z=\sum_{i} x_{i}^2 $}
\bigskip

\subsubsection{The norm of a matrix\_object}
\begin{lstlisting}
number = norm(matrix_object)
\end{lstlisting}
\afterlistingthing{$z=\sqrt{\sum_{ij} x_{ij}^2} $}
\bigskip

\subsubsection{The norm squared of a matrix\_object}
\begin{lstlisting}
number = norm2(matrix_object)
\end{lstlisting}
\afterlistingthing{$z=\sum_{ij} x_{ij}^2 $}
\bigskip

\subsubsection{The sum over the elements of a vector object}
\XX{vector}{sum over the elements of}
\XX{\fontindexentry{tt}{sum}}{operation on a vector}
\begin{lstlisting}
number = sum(vector_object)
\end{lstlisting}
\afterlistingthing{$z=\sum_i x_i $}
\bigskip

\subsubsection{The row sums of a matrix object}
\XX{matrix}{\fontindexentry{tt}{rowsum of}}
\XX{\fontindexentry{tt}{rowsum}}{operation on a matrix}
\begin{lstlisting}
vector = rowsum(matrix_object)
\end{lstlisting}
\afterlistingthing{$z_i=\sum_j x_{ij}$}
\bigskip

\subsubsection{The column sums of a matrix object}
\XX{matrix}{\fontindexentry{tt}{colsum} of}
\XX{\fontindexentry{tt}{colsum}}{operation on a matrix}
\begin{lstlisting}
vector = colsum(matrix_object)
\end{lstlisting}
\afterlistingthing{$z_j=\sum_i x_{ij}$}
\bigskip

\subsubsection{The minimum element of a vector object}
\XX{vector}{minimum element of}
\XX{\fontindexentry{tt}{min}}{operation on a vector}
\begin{lstlisting}
number = min(vector_object)
\end{lstlisting}
\XX{vector}{maximum element of}
\XX{\fontindexentry{tt}{max}}{operation on a vector}
\bigskip

\subsubsection{The maximum element of a vector object}
\begin{lstlisting}
number = max(vector_object)
\end{lstlisting}
\bigskip

\section{Eigenvalues and eigenvectors\br of a symmetric matrix}
\XX{eigenvalues}{not differentiable}
\XX{eigenvectors}{not differentiable}

While we have included eigenvalue and eigenvector routines for both constant and
variable matrix objects, you should be aware that, in general, the eigenvectors
and eigenvalues are not differentiable functions of the variables determining
the matrix.

The eigenvalues of a symmetric matrix
\XX{eigenvalues}{of a symmetric matrix}
\begin{lstlisting}
  vector_object = eigenvalues(matrix_object)
\end{lstlisting}
are returned in a vector. It is the user's responsibility to ensure that the
matrix is actually symmetric. The routine symmetrizes the matrix so that the
eigenvalues returned are actually those for the symmetrized matrix.

The eigenvectors of a symmetric matrix:
\XX{eigenvectors}{of a symmetric matrix}
\begin{lstlisting}
  matrix_object = eigenvectors(matrix_object)
\end{lstlisting}
are returned in a matrix. It is the user's responsibility to ensure that the
matrix is actually symmetric. The routine symmetrizes the matrix so that the
eigenvectors returned are actually those for the symmetrized matrix. The
eigenvectors are located in the columns of the matrix. The $i^\textrm{th}$
eigenvalue returned by the function \texttt{eigenvalues} corresponds to the
$i^\textrm{th}$ eigenvector returned by the function \texttt{eigenvector}.

\X{Choleski decomposition of a symmetric matrix}
\XX{symmetric matrix}{Choleski decomposition}

\section{The Choleski decomposition of a\br positive definite symmetric matrix}

For a positive definite symmetric matrix \texttt{S}, the Choleski decomposition
of \texttt{S} is a lower triangular matrix \texttt{T} satisfying the
relationship \texttt{S=T*trans(T)}. If \texttt{S} is a (positive definite
symmetric) matrix object and \texttt{T} is a matrix object, the line of code
\begin{lstlisting}
  T=choleski_decomp(S);
\end{lstlisting}
will calculate the Choleski decomposition of \texttt{S} and put it into
\texttt{T}.
\X{solving a system of linear equations}
\X{\fontindexentry{tt}{solve} function}
\XX{optimizing performance}{using the best operators for a calculation}

\section{Solving a system of linear equations}

If \texttt{y} is a vector and \texttt{M} is an invertible matrix, then finding a
vector \texttt{x} such that
\begin{lstlisting}
  x=inv(M)*y
\end{lstlisting}
will be referred to as solving the system of linear equations determined by
\texttt{y} and \texttt{M}. Of course, it is possible to use the \texttt{inv}
function to accomplish this task, but it is much more efficient to use the
\texttt{solve} function.
\begin{lstlisting}
  vector x=solve(M,y); // x will satisfy x=inv(M)*y;
\end{lstlisting}

It turns out that it is a simple matter to calculate the determinant of the
matrix \texttt{M} at the same time as the system of linear equations is solved.
Since this is useful in multivariate analysis, we have also included a function
that returns the determinant at the same time as the system of equations is
solved. To avoid floating point overflow or underflow when working with large
matrices, the logarithm of the absolute value of the determinant, together with
the sign of the determinant, are returned. The constant form of the
\texttt{solve} function is
\begin{lstlisting}
double ln_det;
double sign;
dvector x=solve(M,y,ln_det,sign);
\end{lstlisting}
while the variable form is
\begin{lstlisting}
dvariable ln_det;
dvariable sign;
dvar_vector x=solve(M,y,ln_det,sign);
\end{lstlisting}

The \texttt{solve} function is useful for calculating the log-likelihood
function for a multivariate normal distribution. Such a log-likelihood function
involves a calculation similar to
\begin{lstlisting}
  l = -.5*log(det(S)) -.5*y*inv(S)*y
\end{lstlisting}
where \texttt{S} is a matrix object and \texttt{y} is a vector object. It is
much more efficient to carry out this calculation using the \texttt{solve}
function. The following code illustrates the calculations for variable objects:
\XX{multivariate normal distribution}%
{calculation of the log-likelihood function for}
\begin{lstlisting}
dvariable ln_det;
dvariable sign;
dvariable l;
dvar_vector tmp=solve(M,y,ln_det,sign);
l=-.5*ln_det-y*tmp;
\end{lstlisting}

\section{Methods for filling arrays and matrices}
\X{filling arrays and matrices}

While it is always possible to fill vectors and matrices by using loops and
filling them element by element, this is tedious and prone to error. To simplify
this task, a selection of methods for filling vectors and matrices with random
numbers, or a specified sequence of numbers, is available. There are also
methods for filling row and columns of matrices with vectors. In this section
the symbol \texttt{vector} can refer to either a \texttt{dvector} or a
\texttt{dvar\_vector}, while the symbol \texttt{matrix} can refer to either a
\texttt{dmatrix} or a \texttt{dvar\_matrix}.

\XX{vector}{\fontindexentry{tt}{fill}}
\XX{\fontindexentry{tt}{fill}}{filling a vector}
\begin{lstlisting}
void vector::fill("{m,n,...,}")
\end{lstlisting}
\afterlistingthing{fills a vector with a sequence of the form \texttt{n},
  \texttt{m},$\ldots$ The number of elements in the string must match the size
  of the vector.}

\XX{vector}{\fontindexentry{tt}{fill\_seqadd}}
\XX{\fontindexentry{tt}{fill\_seqadd}}{filling a vector}
\begin{lstlisting}
void vector::fill_seqadd(double& base, double& offset)
\end{lstlisting}
\afterlistingthing{fills a vector with a sequence of the form
\texttt{base}, \texttt{base+offset}, \texttt{base+2*offset},$\ldots$}

For example, if \texttt{v} is a \texttt{dvector} created by the statement
\begin{lstlisting}
  dvector v(0,4);
\end{lstlisting}
then the statement
\begin{lstlisting}
  v.fill_seqadd(-1,.5);
\end{lstlisting}
will fill \texttt{v} with the numbers $(-1.0,-0.5,0.0,0.5,1.0)$.

\XX{matrix}{\fontindexentry{tt}{rowfill\_seqadd}}
\XX{\fontindexentry{tt}{rowfill\_seqadd}}{filling a matrix}
\begin{lstlisting}
void matrix::rowfill_seqadd(int& i,double& base, double& offset)
\end{lstlisting}
\afterlistingthing{fills row \texttt{i} of a matrix with a sequence of the form
  \texttt{base}, \texttt{base+offset}, \texttt{base+2*offset},$\ldots$}

\XX{matrix}{\fontindexentry{tt}{colfill\_seqadd}}
\XX{\fontindexentry{tt}{colfill\_seqadd}}{filling a matrix}
\begin{lstlisting}
void matrix::colfill_seqadd(int& j,double& base, double& offset)
\end{lstlisting}
\afterlistingthing{fills column \texttt{j} of a matrix with a sequence of the
  form \texttt{base}, \texttt{base+offset}, \texttt{base+2*offset},$\ldots$}

\XX{matrix}{\fontindexentry{tt}{colfill}}
\XX{\fontindexentry{tt}{colfill}}{filling a matrix column with a vector}
\XX{matrix}{\fontindexentry{tt}{rowfill}}
\XX{\fontindexentry{tt}{rowfill}}{filling a matrix row with a vector}
\begin{lstlisting}
void matrix::colfill(int& j,vector&)
\end{lstlisting}
\afterlistingthing{fills the $\texttt{j}^\textrm{th}$ column of a matrix with a
  vector.}

\begin{lstlisting}
void matrix::rowfill(int& i,vector&)
\end{lstlisting}
\afterlistingthing{fills the $\texttt{i}^\textrm{th}$ row of a matrix with a
  vector.}

In this section, a uniformly distributed random number is assumed to have a
uniform distribution on $[0,1]$. A normally distributed random number is assumed
to have mean $0$ and variance~$1$. A binomially distributed random number is
assumed to have a parameter~$p$, where $1$ is returned with probability~$p$,
and~$0$ is returned with probability $1-p$. A multinomially distributed random
variable is assumed to have a vector of parameters $P$, where $i$ is returned
with probability~$p_i$. If the components of $P$ do not sum to $1$, the vector
will be normalized so that the components \textit{do} sum to~$1$.

\XX{vector}{\fontindexentry{tt}{fill\_randu}}
\XX{\fontindexentry{tt}{fill\_randu}}{filling a vector with random numbers}
\begin{lstlisting}
void vector::fill_randu(long int& n)
\end{lstlisting}
\afterlistingthing{fills a vector with a sequence of uniformly distributed
  random numbers. The \texttt{long int n} is a seed for the random number
  generator. Changing \texttt{n} will produce a different sequence of random
  numbers.}

\XX{matrix}{\fontindexentry{tt}{colfill\_randu}}
\XX{\fontindexentry{tt}{colfill\_randu}}{filling a matrix with random numbers}
\begin{lstlisting}
void matrix::colfill_randu(int& j,long int& n)
\end{lstlisting}
\afterlistingthing{fills column \texttt{j} of a matrix with a sequence of
  uniformly distributed random numbers. The \texttt{long int n} is a seed for
  the random number generator. Changing \texttt{n} will produce a different
  sequence of random numbers.}

\XX{matrix}{\fontindexentry{tt}{rowfill\_randu}}
\XX{\fontindexentry{tt}{rowfill\_randu}}{filling a matrix with random numbers}
\begin{lstlisting}
void matrix::rowfill_randu(int& i,long int& n)
\end{lstlisting}
\afterlistingthing{fills row \texttt{i} of a matrix with a sequence of uniformly
  distributed random numbers.}

\XX{vector}{\fontindexentry{tt}{fill\_randbi}}
\XX{\fontindexentry{tt}{fill\_randbi}}{filling a vector with random numbers}
\begin{lstlisting}
void vector::fill_randbi(long int& n, double& p)
\end{lstlisting}
\afterlistingthing{fills a vector with a sequence random numbers from a binomial
  distribution.}

\XX{vector}{\fontindexentry{tt}{fill\_randn}}
\XX{\fontindexentry{tt}{fill\_randn}}{filling a vector with random numbers}
\begin{lstlisting}
void vector::fill_randn(long int& n)
\end{lstlisting}
\afterlistingthing{fills a vector with a sequence of normally distributed random
  numbers.}

\XX{matrix}{\fontindexentry{tt}{rowfill\_randn}}
\XX{\fontindexentry{tt}{rowfill\_randn}}{filling a matrix with random numbers}
\begin{lstlisting}
void matrix::colfill_randn(int& j,long int& n)
\end{lstlisting}
\afterlistingthing{fills column \texttt{j} of a matrix with a sequence of
  normally distributed random numbers.}

\begin{lstlisting}
void matrix::rowfill_randn(int& i,long int& n)
\end{lstlisting}
\afterlistingthing{fills row \texttt{i} of a matrix with a sequence of normally
  distributed random numbers.}

\XX{vector}{\fontindexentry{tt}{fill\_multinomial}}
\XX{\fontindexentry{tt}{fill\_multinomial}}%
{filling a vector with random numbers}
\begin{lstlisting}
void vector::fill_multinomial(long int& n, dvector& p)
\end{lstlisting}
\afterlistingthing{fills a vector with a sequence random numbers from a
  multinomial distribution. The parameter~\texttt{p} is a \texttt{dvector} such
  that \texttt{p[}$i$\texttt{]} is the probability of returning $i$. The
  elements of \texttt{p} must sum to $1$.}

\section{Methods for extracting from arrays and matrices}
\X{extracting data from arrays and matrices}
\XX{matrix}{\fontindexentry{tt}{extract\_column}}
\XX{\fontindexentry{tt}{extract\_column}}{from a matrix}

\begin{lstlisting}
vector extract_column(matrix& M,int& j)
\end{lstlisting}
\afterlistingthing{extracts a column from a matrix and puts it into a vector.}

\begin{lstlisting}
vector extract_row(matrix& M,int& i)
\end{lstlisting}
\afterlistingthing{extracts a row from a matrix and puts it into a vector.}

\XX{\fontindexentry{tt}{extract\_diagonal}}{from a matrix}
\begin{lstlisting}
vector extract_diagonal(matrix& M)
\end{lstlisting}
\afterlistingthing{extracts the diagonal elements from a matrix and puts them
  into a vector.}

\X{operator \texttt{()}}
The function call operator (\texttt{()}) has been overloaded in two ways to
provide for the extraction of a subvector. X{vector}{extracting a subvector}
\X{extracting a subvector}
\XX{matrix}{\fontindexentry{tt}{extract\_row}}
\XX{\fontindexentry{tt}{extract\_row}}{from a matrix}
\XX{vector}{function call \texttt{()} to extract subvector}
\begin{lstlisting}
  vector(ivector&)
\end{lstlisting}
An \texttt{ivector} object is used to specify the elements of the vector to be
chosen. If \texttt{u} and \texttt{v} are \texttt{dvectors} and \texttt{i} is an
\texttt{ivector}, then the construction
\begin{lstlisting}
  dvector u = v(i);
\end{lstlisting}
will extract the members of \texttt{v} indexed by \texttt{i} and put them in the
\texttt{dvector u}. The size of~\texttt{u} is equal to the size of~\texttt{i}.
The \texttt{dvector u} will have minimum valid index and maximum valid index
equal to the minimum valid index and maximum valid index of \texttt{i}. The size
of~\texttt{i} can be larger than the size of~\texttt{v}, in which case some
elements of \texttt{v} must be repeated. The elements of the \texttt{ivector i}
must lie in the valid index range for \texttt{v}.

If \texttt{v} is a \texttt{dvector} and \texttt{i1} and \texttt{i2} are two
integers,
\begin{lstlisting}
  u(i1,i2)
\end{lstlisting}
is a \texttt{dvector}, which is a subvector of \texttt{v} (provided, of course,
that \texttt{i1} and \texttt{i2} are valid indices for \texttt{v}). Subvectors
can appear on both the left and right-hand side of an assignment:
\begin{lstlisting}
dvector  u(1,20);
dvector v(1,19);
v = 2.0;   // assigns the value 2 to all elements of v
u(1,19) = v; // assigns the value 2 to elements 1 through 19 of u
\end{lstlisting}
\X{operator $++$}
\X{operator $-{}-$}
\XX{operator $++$}{use with subvectors}
\XX{operator $-{}-$}{use with subvectors}
\XX{operator $++$}{for \texttt{dvectors}}
\XX{operator $-{}-$}{for \texttt{dvectors}}

In the above example, suppose that we wanted to assign the vector \texttt{v} to
elements 2 through 20 of the vector \texttt{u}. To do this, we must first ensure
that they have the same valid index ranges. The operators \texttt{++} and
\texttt{-{}-}, when applied to vector objects, increment and decrement the index
ranges by~1. See the code fragment below:
\begin{lstlisting}
dvector u(1,20);
dvector v(1,19);
v = 2.0;   // assigns the value 2 to all elements of v
--u(2,20) = v; // assigns the value 2 to elements 2 through 20 of u
u(2,20) = ++v; // assigns the value 2 to elements 2 through 20 of u
\end{lstlisting}
It is important to realize that from the point of view of the vector \texttt{v},
both of the above assignments have the same effect. It will have elements 2
through 20 set equal to~2. The difference is the side effects on the vector
\texttt{v}. The use of subvectors, and increment and decrement operations, can
be used to remove loops from the code. Note that
\XX{subvectors}{using to remove loops from code}
\begin{lstlisting}
dvector x(1,n)
dvector y(1,n)
dvector z(1,n)
for (int i=2;i<=n;i++)
{
  x(i)=y(i-1)*z(i-1);
}
\end{lstlisting}
can be written as
\begin{lstlisting}
dvector x(1,n)
dvector y(1,n)
dvector z(1,n)
x(2,n)=++elem_prod(y(1,n-1),z(1,n-1));  // elem_prod is element-wise
                                       // multiplication of vectors
\end{lstlisting}

\section{Making a bootstrap sample}
X{bootstrap sample}{code for}

As an example of how one can use some of these operations, consider the problem
of creating bootstrap samples from a set of data. Suppose we wish to create a
random sample (with replacement) of the original data set. Assume that the
original data is contained in a \texttt{dvector} object named \texttt{data}. It
is not necessary that the new samples be the same size as the original sample,
so we shall allow the size of the bootstrap sample to be a variable. The code
for the function \texttt{bootstrap} follows:
\begin{lstlisting}
dvector bootstrap(dvector& data, int& sample_size, long int seed)
{
  // data contains the data which is used to make the bootstrap sample
  // sample_size is the size of the desired bootstrap sample
  // seed is a seed for the random number generator
  dvector tmp(1,sample_size);
  tmp.fill_randu(seed);  // 0<tmp(j)<1             for 1<j<sample_size
  tmp=tmp*data.size()+1; // 1<tmp(j)<data.size()+1 for 1<j<sample_size
  ivector iselect(tmp);  // Take the integer parts of tmp and put them
                         // into iselect
                         // 1<=iselect(j)<=data.size() for 1<=j<=sample_size
  tmp=data(iselect);     // Select the bootstrap sample and
                         // store it in tmp
  return(tmp);
}
\end{lstlisting}

The above code could be employed in the user's code to make bootstrap samples by
using something like the following code fragment:
\begin{lstlisting}
// ...
dvector data(1,nobs)
// ... somehow the data is put into data
// ...
dvector boot_sample=bootstrap(data,200,1231); // Make a sample of size 200
                          // and create the dvector boot_sample to hold it
// ...
\end{lstlisting}

\section{Sorting vectors and matrices}
\X{sorting}
\XX{\fontindexentry{tt}{dvector}}{sorting a}
\XX{\fontindexentry{tt}{dmatrix}}{sorting a}
\XX{sorting}{\fontindexentry{tt}{dvector}}
\XX{sorting}{\fontindexentry{tt}{dmatrix}}

While sorting is not strictly part of the methods for calculating the
derivatives of differentiable functions (as it is a highly non-differentiable
operation), it is so useful for pre- and post-processing data that we have
included some functions for sorting \texttt{dvector} and \texttt{dmatrix}
objects. If \texttt{v} is a \texttt{dvector} the statement
\begin{lstlisting}
  dvector w=sort(v);
\end{lstlisting}
will sort the elements of \texttt{v} in ascending order and put them in the
\texttt{dvector} object \texttt{w}. The minimum and maximum valid indices of
\texttt{w} will be the same as those of \texttt{v}. If desired, an index table
for the sort can be constructed by passing an \texttt{ivector} along with the
\texttt{dvector}. This index table can be used to sort other vectors in the same
order as the original vector using the \texttt{()} operator.
\X{operator \texttt{()}}
\begin{lstlisting}
dvector u={4,2,1};
dvector v={1,6,5}
ivector ind(1,3);
dvector w=sort(u,ind); // ind will contain an index table for the sort
// Now w=(1,2,4) and ind=(3,2,1)
dvector ww=v(ind);     // This is the use of the ( ) operator for subset
                       // selection.
// Now ww=(5,6,1)
\end{lstlisting}

The \texttt{sort} function for a \texttt{dmatrix} object sorts the columns of
the \texttt{dmatrix} into ascending order, using the column specified to do the
sorting. For example,
\begin{lstlisting}
  dmatrix MM = sort(M,3);
\end{lstlisting}
will put the sorted matrix into \texttt{MM}, and the third column of \texttt{MM}
will be sorted in ascending~\mbox{order}.

\section{Mathematical functions}

The following functions have been included in \scAD\ by overloading the \cplus\
library functions:
\begin{lstlisting}
  sin cos tan asin atan acos sinh cosh tanh fabs (sfabs) exp log log10 sqrt pow
\end{lstlisting}
\X{\fontindexentry{tt}{sin}}
\X{\fontindexentry{tt}{cos}}
\X{\fontindexentry{tt}{tan}}
\X{\fontindexentry{tt}{asin}}
\X{\fontindexentry{tt}{atan}}
\X{\fontindexentry{tt}{acos}}
\X{\fontindexentry{tt}{sinh}}
\X{\fontindexentry{tt}{cosh}}
\X{\fontindexentry{tt}{tanh}}
\X{\fontindexentry{tt}{fabs}}
\X{\fontindexentry{tt}{sfabs}}
\X{\fontindexentry{tt}{exp}}
\X{\fontindexentry{tt}{log}}
\X{\fontindexentry{tt}{log10}}
\X{\fontindexentry{tt}{sqrt}}
\X{\fontindexentry{tt}{pow}}
These functions can be used on numbers or \texttt{vector\_objects} in the form
\begin{lstlisting}
number = function(number);
vector_object = function(vector_object);
\end{lstlisting}
When operating on \texttt{vector\_objects}, the functions operate element by
element. So, if \texttt{y} is a \texttt{dvector} whose elements are
$(y_1,\ldots,y_n)$, then \texttt{exp(y)} is a \texttt{dvector} whose elements
are $(\exp(y_1),\ldots,\exp(y_n))$.

\XX{vector}{maximum element of}
\XX{\fontindexentry{tt}{max}}{operation on a vector}
\XX{vector}{minimum element of}
\XX{\fontindexentry{tt}{min}}{operation on a vector}
The functions \texttt{min} and \texttt{max}, when applied to a
\texttt{vector\_object}, return a number that is equal to the minimum or maximum
element of the \texttt{vector\_object}.

\section{Using the smoothed absolute value function \texttt{sfabs}}

The function \texttt{sfabs} is a ``smoothed'' absolute value function that
agrees with the standard absolution function for all values of its argument with
absolute value greater than 0.001. For values of its argument less than 0.001 in
absolute value, the function uses cubic splines to produce a differentiable
approximation to the absolute value function. This function is intended to be
used when a differentiable function that approximates the absolute value
function is needed.

While the function \texttt{sfabs} has been set to agree with the absolute value
function \texttt{fabs} for all values greater than 0.001, it is a simple matter
to rescale the function to agree with the absolute value function over any
desired range, by multiplying and dividing it by the desired scale factor:
\begin{lstlisting}
double scale
double x
// ...
double v=scale*sfabs(x/scale); // This function will agree with fabs
                            // for all components of w with absolute value
                            // greater than scale*.001
\end{lstlisting}
\begin{figure}[h]
  \centering\hskip1pt
 % \vbox{\hsize=.95\textwidth
  \hskip1pt\input robgrph3.tex%}
  \caption{Plot of the smoothed absolute value function
    \texttt{y=1000*sfabs(x/1000)} and the normal  absolute value function
    \texttt{y=fabs(x)}.}
  \label{fig:operations-functions-01}
\end{figure}
In Figure~\ref{fig:operations-functions-01}, the function \texttt{sfabs} has
been ``magnified'' by a scale factor of 1000, to show how it differs from the
ordinary absolute value function \texttt{fabs}.

%xx \htmlnewfile

\chapter{Advanced Concepts}
\label{ch:advanced-concepts}

This chapter is devoted to various matters that enable the \scAD\ user to create
more advanced applications, and to obtain faster performance from the \scAD\
code created. Beginning users can skip or briefly skim the topics contained
here.
\XX{optimization considerations}{temporary files}

\section{Reducing the size of the temporary files}

A disadvantage of the reverse mode of automatic differentiation is the large
amount of temporary storage required for keeping the intermediate data for the
derivative calculations. Every arithmetic operation on a \texttt{dvariable}
object generates 28--30 bytes of temporary storage. To see what this means,
consider that to invert an $n\times n$ matrix requires about $n^3$ such
operations. So, the temporary storage required for inverting a $100\times100$
matrix would be about 28--30~megabytes.

To reduce the requirements for temporary storage, \scAD\ employs precompiled
derivative code for common operations on variable container classes. Depending
on the operation, this can greatly reduce the temporary storage used. For the
matrix inverse, the storage for a $100\times 100$ matrix is about 500K, a
reduction of about~60 to~1. For a vector dot product, the reduction in temporary
storage requirements is about~3.5 to~1. While this is much less than the matrix
inverse, it can still be a considerable saving in large problems.

To achieve the best performance, you should use \scAD's vector and array
operations on variable objects whenever possible, rather than writing your own
code to manipulate elements of these container classes. Consider the following
code fragment which ``normalizes'' a \texttt{dvar\_vector} so that its
components sum to \texttt{1}:
\begin{lstlisting}
dvar_vector v(1,n);
// get some data into v somehow
// ...
// now normalize v so that its components sum to 1
dvariable sum=0.;
for (int i=1;i<=n;i++)
{
  sum+=v(i);
}
for (i=1;i<=n;i++)
{
  v(i)=v(i)/sum;
}
\end{lstlisting}
\X{sum@\texttt{sum}}
\XX{optimizing performance}{using the best operators for a calculation}

The same normalization can be achieved by the following code. It uses the \scAD\
function \texttt{sum}, which sums the components of a vector object. It also
uses the divide operator, for dividing a vector by a number. It will produce
less temporary derivative information.
\begin{lstlisting}
dvar_vector v(1,n);
// get some data into v somehow
// ...
// now normalize v so that its components sum to 1
v=v/sum(v);
\end{lstlisting}

Even more efficient code is generated by employing the \texttt{/=} operator,
which divides a vector by a number:
\begin{lstlisting}
  v/=sum(v);
\end{lstlisting}

\XX{optimization considerations}{accessing container class elements}

\section{Creating efficient code when accessing arrays}

Consider the following code fragment for adding two matrix objects:
\begin{lstlisting}
dvar_matrix a(1,n,1,n);
dvar_matrix b(1,n,1,n);
dvar_matrix c(1,n,1,n);
// ...
for (int i=1;i<=n;i++)
{
  for (int j=1;j<=n;j++)
  {
    a(i,j)=b(i,j)+c(i,j);
  }
}
\end{lstlisting}

A \textsc{fortran} compiler that only has to deal with a small number of array
types can be expected to produce well-optimized code from the above source code.
\cplus\ compilers are less likely to be able to exploit that fact that the first
index \texttt{i} is constant in the inner loop. It is to be expected that, in
the future, the ability of \cplus\ compilers to optimize code will improve. In
the meantime, code segments such as the above can be optimized by the user
through the use of ``references.'' If \texttt{a} is a \texttt{dvar\_matrix},
then \texttt{a(i)} (or \texttt{a[i]}) is a \texttt{dvar\_vector}. This fact can
be used to reduce the redundant addressing operations on the first index in the
above example, by defining reference objects for these \texttt{dvar\_vectors}.
If you are not familiar with the concept of references, which are extremely
useful in \cplus, consult your \cplus\ reference manual.
\begin{lstlisting}
dvar_matrix a(1,n,1,n);
dvar_matrix b(1,n,1,n);
dvar_matrix c(1,n,1,n);
// ...
for (int i=1;i<=n;i++)
{
  dvar_vector& ai=a(i);
  dvar_vector& bi=b(i);
  dvar_vector& ci=c(i);
  for (int j=1;j<=n;j++)
  {
    ai(j)=bi(j)+ci(j);
  }
}
\end{lstlisting}
Note, however, that the most efficient way to add two matrices is to use the
code
\begin{lstlisting}
  a=b+c;
\end{lstlisting}

Of course, there is some overhead involved in establishing the references, so
this would not be worthwhile for very small values of \texttt{n}.

\section{Ragged matrices and 3-dimensional arrays}

\XX{data structures}{ragged matrices}
\X{ragged matrices}
\XX{data structures}{ragged 3-dimensional arrays}
\X{ragged 3-dimensional arrays}
A big advantage of \cplus\ over traditional scientific programming languages,
such as \textsc{fortran}, is the increased facility for defining and
manipulating complex data structures. Many real-world situations cannot be
conveniently described in terms of the simple ``flat'' data structures
represented by vectors, matrices, and regular 3-dimensional arrays.

As an aid to creating more complex data structures, \scAD's 2- and 3-dimensional
container classes have been extended to what we call ``ragged arrays.'' In
short, a ragged matrix is a matrix whose rows consist of vectors of different
lengths and varying valid index ranges. A ragged 3-dimensional array can be
viewed as an array of matrices of different sizes, but such matrices may or may
not be ragged themselves.

\section{Specifying ragged matrices}

To specify a ragged matrix requires two integers for the minimum and maximum
valid row indices, a vector of integers for the minimum valid index of the
vector forming each row of the matrix, and a vector of integers for the maximum
valid index of the vector forming each row of the matrix.
\begin{lstlisting}
int min=0;
int max=4;
ivector minind(0,4);
ivector maxind(0,4);
// get the desired values into the ivectors minind and maxind
// ...
dmatrix M(min,max,minind,maxind);
\end{lstlisting}

Notice that the minimum and maximum valid indices for \texttt{minind} and
\texttt{maxind} must be the same as the minimum and maximum valid row indices
for the matrix \texttt{M}.

Suppose that $\texttt{minind(3)}=-1$ and $\texttt{maxind(3)}=5$. Then the fourth
row of \texttt{M} will be a vector with seven components whose minimum valid
index is $-1$ and whose maximum valid index is~$5$.

Often one desires to make a ragged matrix where the minimum valid index for all
the rows is the same. This can be accomplished by filling the \texttt{ivector}
\texttt{minind} with the same number in each component. For convenience, a
matrix constructor that takes an \texttt{int} for its third argument has been
supplied.
\begin{lstlisting}
int min=0;
int max=4;
int minind=1;
ivector maxind(0,4);
// get the desired values into the ivector maxind
dmatrix M(min,max,minind,maxind);
\end{lstlisting}
This will make a ragged matrix \texttt{M}. The minimum valid indices for vectors
making up the rows of \texttt{M} will all be equal to~1.

\section{Specifying ragged 3-dimensional arrays}

The range of constructors available for ragged 3-dimensional arrays is
considerably more complicated than for ragged matrices. The most general
constructor for a ragged 3-dimensional array has the following prototype:
\begin{lstlisting}
  d3_array(int sl,int su,ivector& rl,ivector& ru,imatrix& cl,imatrix& cu);
\end{lstlisting}

The $i^\textrm{th}$ slice of the \texttt{d3\_array} is a ragged \texttt{dmatrix}
that is constructed using the data \texttt{rl($i$)}, \texttt{ru($i$)},
\texttt{cl($i$)}, and \texttt{cu($i$)}, as described above. Notice that
\texttt{rl($i$)} and \texttt{ru($i$)} are of type \texttt{int}, while
\texttt{cl($i$)} and \texttt{cu($i$)} are objects of type \texttt{ivector}. So,
this all makes sense. This is the most ``ragged'' of 3-dimensional objects.
There is also a constructor with the following prototype:
\begin{lstlisting}
  d3_array(int sl,int su,ivector& rl,ivector& ru,ivector& cl,ivector& cu);
\end{lstlisting}
The $i^\textrm{th}$ slice of the \texttt{d3\_array} will use the data
%\texttt{rl(i)\ ru(i)\ cl(i)\ cu(i)}
\texttt{rl($i$)}, \texttt{ru($i$)}, \texttt{cl($i$)}, and \texttt{cu($i$)} to
make a \texttt{dmatrix}. Since these arguments are all of type \texttt{int},
each slice of the ragged 3-dimensional array consists of a ``non-ragged'' or
regular \texttt{matrix} object. An example of the use of a ragged 3-dimensional
array for a data structure is given in Chapter~\ref{ch:neural-networks}. The
complete list of constructors provided for specifying 3-dimensional arrays is:
\begin{lstlisting}
d3_array( int& sl,  int& sh,  int& nrl,
       int& nrh, int& ncl,  int& nch);

d3_array(int& sl,int& sh,ivector& nrl,ivector& nrh,
      imatrix& ncl,imatrix& nch);

d3_array(int& sl,int& sh,ivector& nrl,ivector& nrh,
      int ncl,imatrix& nch);

d3_array(int& sl,int& sh,ivector& nrl,ivector& nrh,
      ivector& ncl,ivector& nch);

d3_array(int& sl,int& sh,int& nrl,ivector& nrh,
      int& ncl,ivector& nch);

d3_array(int& sl,int& sh,int& nrl,ivector& nrh,
      int& ncl,int& nch);

d3_array(int sl,int sh,int nrl,ivector& nrh,
      int ncl,imatrix& nch);
\end{lstlisting}

\section{Complex data structures:\br making order out of chaos}

As your model becomes more complex, the number of \scAD\ objects required will
increase until the code again starts to become unmanageable. One solution is to
group these objects themselves into structures or classes, so they can be
manipulated as single objects.

Suppose that your model has five vector objects whose minimum valid indices
are~1 and whose maximum valid indices are all equal to an integer $n$, whose
value is known at run time. Denote these vector objects by \texttt{a1},
\texttt{a2},$\ldots$, \texttt{a5}. Let the model have five more vector objects
whose minimum valid indices are~1 and whose maximum valid indices are~$2n$.
Denote these vector objects by \texttt{b1}, \texttt{b2},$\ldots$, \texttt{b5}.
Finally, suppose the model has a (regular) matrix object whose valid row indices
go from~1 to~$n$ and whose valid column indices go from~0 to~$2n$. Denote this
matrix by \texttt{M}. We can encapsulate all these objects in a class as
follows.
\begin{lstlisting}
// class definition for constant_model_parameters
class constant_model_parameters
{
public:
  int size;
  dvector a1;
  dvector a2;
  dvector a3;
  dvector a4;
  dvector a5;
  dvector b1;
  dvector b2;
  dvector b3;
  dvector b4;
  dvector b5;
  dmatrix M;
  constant_model_parameters(int n); // prototype for constructor
};
\end{lstlisting}

At the point in your program where you want to ``create'' this object, you will
need to invoke the constructor described below. The invocation is simply a
declaration of the object with its name, such as:
\begin{lstlisting}
int n=10;
// ...
constant_model_parameters cmp(n); //creates the object with name "cmp"
\end{lstlisting}

The main part of the implementation of the class is in the definition of the
constructor.
\begin{lstlisting}
constant_model_parameters::constant_model_parameters(int n) :
   a1(1,n),a2(1,n),a3(1,n),a4(1,n),a5(1,n),
   b1(1,2*n),b2(1,2*n),b3(1,2*n),b4(1,2*n),b5(1,2*n),
   M(1,n,0,2*n) {size=n;}
\end{lstlisting}
If you are not familiar with the form of the constructor above, you should
consult a \cplus\ reference book. The idea is that the constructors for all the
classes that are members of the class \texttt{constant\_model\_parameters} are
invoked after the colon (\texttt{:}) and before the body of the constructor.

Note that while for the ``safe'' library all container classes are initialized
to zero, for the optimized library, they are not. If you want to initialize some
of these objects, such as \texttt{a1} and \texttt{M}, this should be done in the
body of the constructor.
\begin{lstlisting}
constant_model_parameters::constant_model_parameters(int n) :
   a1(1,n),a2(1,n),a3(1,n),a4(1,n),a5(1,n),
   b1(1,2*n),b2(1,2*n),b3(1,2*n),b4(1,2*n),b5(1,2*n),
   M(1,n,0,2*n)
{
   size=n;
   a1.initialize();  // set components of a1 and M to 0
   M.initialize();
}
\end{lstlisting}

All the objects in the class \texttt{constant\_model\_parameters} are as the
name implies: constants. This class could be useful for reading in the
parameters from a file, for example, and otherwise manipulating them. At some
point, presumably the user wants the parameters transformed into variable
objects so that derivative information can be generated. This can be done by
defining a corresponding variable class:

\begin{lstlisting}
// class definition for variable_model_parameters
class variable_model_parameters
{
public:
  dvar_vector a1;
  dvar_vector a2;
  dvar_vector a3;
  dvar_vector a4;
  dvar_vector a5;
  dvar_vector b1;
  dvar_vector b2;
  dvar_vector b3;
  dvar_vector b4;
  dvar_vector b5;
  dvar_matrix M;
  // This constructor transforms the constant parameters into
  // variable parameters
  variable_model_parameters(constant_model_parameters&); // prototype
                                                     // for constructor
};
\end{lstlisting}

The constructor for the class \texttt{variable\_model\_parameters} takes the
constant objects in the class \texttt{constant\_model\_parameters} and converts
them to variable objects.
\begin{lstlisting}
variable_model_parameters::variable_model_parameters
   (constant_model_parameters& cms) :
   a1(cms.a1),a2(cms.a2),a3(cms.a3),a4(cms.a4),a5(cms.a5),
   b1(cms.b1),b2(cms.b2),b3(cms.b3),b4(cms.b4),b5(cms.b5),
   M(cms.M) {size=cms.size;}
\end{lstlisting}

After all this has been done, the constant and variable classes could be used,
as~in:
\begin{lstlisting}
void fcomp(variable_model_parameters&); // function prototype for fcomp
int n=10;
constant_model_parameters cms(n);
// get some data into cms
//
// call a function which takes an object of type variable_model_parameters
fcomp(cms); // cms will be converted into an object of type
            // variable_model_parameters
\end{lstlisting}
Since in this example, the only constructor that has been defined for an object
of type \texttt{variable\_model\_parameters} takes an argument of type
\texttt{constant\_model\_parameters}, the only way to create the variable object
will be to first create the constant object. This may or may not be what is
required. Of course, a constructor can be defined to create a variable object
from the \texttt{int n}.

\section{Default constructors for complex data structures}

You should be aware that the \cplus\ compiler attempts to create default
constructors for complex classes if they are not provided by the user. To do
this, it employs the constructors provided for the members of the complex class.
For example, if you employ the lines of code
\begin{lstlisting}
constant_model_parameters cmp(10);
constant_model_parameters cmp1=cmp;
\end{lstlisting}
the compiler use a default copy constructor for the class
\texttt{constant\_model\_parameters} (assuming that you, the user, have not
already defined one). In order to create the default copy constructor, the
compiler will invoke the copy constructors provided by \scAD\ for
\texttt{dvector} and \texttt{dmatrix} classes. Since these constructors perform
shallow copies, the \texttt{dvector} and \texttt{dmatrix} objects in
\texttt{cmp1} will not be distinct from the corresponding objects \texttt{cmp}.
If you want to define a copy constructor for the class
\texttt{constant\_model\_parameters} that will create a distinct object, it
could be done as follows:
\begin{lstlisting}
constant_model_parameters::
  constant_model_parameters(constant_model_parameters cms) :
  a1(1,cms.size),a2(1,cms.size),a3(1,cms.size),a4(1,cms.size),
  a5(1,cms.size),b1(1,2*cms.n),b2(1,2*cms.n),b3(1,2*cms.n),
  b4(1,2*cms.n),b5(1,2*cms.n),
  M(1,cms.size,0,2*cms.size)
{
  size=cms.size;
  a1=cms.a1;
  a2=cms.a2;
  a3=cms.a3;
  a4=cms.a4;
  a5=cms.a5;
  b1=cms.b1;
  b2=cms.b2;
  b3=cms.b3;
  b4=cms.b4;
  b5=cms.b5;
  M=cms.M;
}
\end{lstlisting}

%xx \htmlnewfile

\chapter{The \scAD\ Libraries}

\XX{\scAD}{libraries}
\XX{libraries}{safe}
\XX{libraries}{optimized}
\scAD\ is supplied with safe and optimized versions of its object module
libraries. The safe version is intended to be used for code development. It
provides for bounds checking of arrays and initialization of all declared
objects. This initialization sets the value of all array members to zero and
initializes the derivative structure associated with each variable type object.
The initialization of the derivative structures can involve a fair amount of
overhead, so after a project has been debugged, the user may wish to use the
optimized version of the library. The optimized version does not provide array
bounds checking or initialization of objects. It is recommended that during the
development stage of the program, the user initialize all objects that require
initialization---even though this is not necessary when using the safe library.
The project code can be easily linked with the optimized library later if
desired. If an application runs successfully using the safe library, but does
not run properly using the optimized library, there is probably an uninitialized
object being used somewhere in the code.

The names and number of libraries supplied with \scAD\ vary from compiler to
compiler, so check the \texttt{read.me} file on the distribution disk for the
names of the libraries supplied with your version of \scAD.

\X{initializing \scAD\ classes}
The class member function \texttt{initialize} is provided for initializing all
appropriate \scAD\ classes. For example, to declare and initialize a
\texttt{dvar\_matrix M}, you could use the following code fragment:
\begin{lstlisting}
dvar_matrix M(1,20,1,30) // Creates a 20 by 30 matrix
M.initialize();          // Sets all the elements of M to 0 and initializes
                         // the derivative structures for M
\end{lstlisting}

An example of code where an object requires initialization is:
\begin{lstlisting}
{
  // . . .
  dvariable u; //u is declared but not initialized
  for (int i=1;i<=n;i++)
  {
    u+=x[i];  // x is a vector object which was previously declared
  }            // This is an error for the optimized library
 }
\end{lstlisting}

In this example, the \texttt{dvariable u} is declared, but not initialized. The
use of the uninitialized \texttt{u} in the expression \texttt{u+=x[i]} is an
error if the code is linked with the optimized library. It will perform fine if
the code is linked with the safe library. To use the optimized library, this
code could be rewritten as:
\begin{lstlisting}
{
  // . . .
  dvariable u=0; //u is declared and initialized
  for (int i=1;i<=n;i++)
  {
    u+=x[i];  // x is a vector object which was previously declared
  }
}
\end{lstlisting}
The code could also be written as:
\begin{lstlisting}
{
  // . . .
  dvariable u; //u is declared but not initialized
  u=0;         //u is initialized
  for (int i=1;i<=n;i++)
  {
    u+=x[i];  // x is a vector object which was previously declared
  }
}
\end{lstlisting}
It is slightly less efficient to initialize \texttt{x} in this way.

For the \textsc{borland} \cplus\ compiler, both libraries have been compiled to
enable overlays.

\XX{overlays}{Borland \cplus\ compiler and}
\X{Borland}

%xx \htmlnewfile

\chapter{Input and Output}
\X{input and output}
\X{iostream}
\X{operators \texttt{<{}<} and \texttt{>{}>}}
\XX{vector}{I/O operations}

Input and output methods for the \scAD\ class objects have been implemented by
overloading of the \cplus\ stream I/O operators \texttt{<{}<} and \texttt{>{}>}
for all \scAD\ classes. The so-called \textit{iostream} package, introduced in
AT\&T~Release 2.0 of \cplus, is fully supported. For a detailed discussion of
\cplus\ stream I/O, you should consult your \cplus\ manual (and the appropriate
header files). In this section, we discuss extension of stream I/O to \scAD\
class objects. Formatted input and output uses the standard \cplus\ stream
classes \texttt{istream}, \texttt{ostream}, \texttt{ifstream}, and
\texttt{ofstream}. In addition, unformatted (``binary'') file I/O is implemented
through two additional stream classes \texttt{uostream} and \texttt{uistream}.

\section{Formatted stream I/O}
\XX{input and output}{formatted}
\X{formatting}
\XX{manipulators}{\fontindexentry{tt}{setw}}
\XX{manipulators}{\fontindexentry{tt}{setfixex}}
\XX{manipulators}{\fontindexentry{tt}{setscientific}}
\XX{manipulators}{\fontindexentry{tt}{setprecision}}

Variables of type \adtypes\ may be input or output using the standard stream
operators, and formatted using the standard manipulators. A two additional
manipulators are provided to simplify selection of either fixed point or
``scientific'' (e-format) notation. The use of the manipulators \texttt{setw()}
(set width), \texttt{setprecision()}, \texttt{setfixed}, \texttt{setscientific},
and \texttt{endl} is illustrated below.
\begin{lstlisting}
#include <iomanip.h>
#include <fvar.hpp>

main()
{
  int n = 30;
  int m = 5;
  gradient_structure gs;
  dvar_matrix x(1, n, 1, m);
  .
  .   // code to calculate x
  .
  // e format with default width and precision
  cout << setscientific << x << endl;

  // fixed point format
  cout << setw(13) << setprecision(3) << setfixed << x << endl;
  .
  .
  .
}
\end{lstlisting}
The first output statement will produce 30 lines with five numbers per line in
``scientific'' (e-format) separated by one or more blanks. The second output
statement will produce a table with 30 rows and 5 columns, with each member of
\texttt{x} printed with a fixed decimal place, with 3 significant figures in a
field 13 bytes wide. If the manipulators are omitted, the output will appear
with the default width, precision, and format specified by the compiler, and
with each number separated by one space.

The following standard manipulators are also supported:
\XX{manipulators}{\fontindexentry{tt}{setf}}
\begin{lstlisting}
setf(f, ios::floatfield);  // to select fixed or scientific
setf(b, ios::basefield);   // to base of number
setf(a, ios::adjustfield); // to select justification of fields
// where
//   f  may be ios::fixed or ios::scientific
//   b  may be ios::dec, ios::oct or ios::hex
//   a  may be ios::left, ios::right or ios::internal
\end{lstlisting}

%Some compilers to not yet support \texttt{iostream} package. The \scAD\ system
%includes the overloaded function \texttt{fform( )} which mimics the
%older style \cplus\ formatting function \texttt{form(const char* ... )}.
%The function prototypes for \texttt{fform( )} are:
%\begin{lstlisting}
%char* fform(const char*, dvariable&);
%char* fform(const char*, dvar_vector&);
%char* fform(const char*, dvar_matrix&);
%char* fform(const char*, dvector&);
%char* fform(const char*, dmatrix&);
%\end{lstlisting}
%\noindent The first argument is any valid C format string as described
%in the documentation for \texttt{printf}. The second
%argument is the object to be printed.
%The function \texttt{fform( )} is noticeably slower than the \texttt{iostream}
%package and the size of the object it can print is limited by
%by the size of an internal buffer. It offers no known advantages
%over the \texttt{iostream} package and is only provided to enable some
%control over the format of the output.
%{\bf There are incompatibilities between the new iostream package and the
%older stream package. Users are cautioned against trying to use both.}

\section{Error checking}
\X{I/O error checking}
\X{operator \texttt{!}}
\XX{operator \texttt{!}}{I/O error checking and}

Although the stream operators \texttt{>{}>} and \texttt{<{}<} make it very easy
to perform input or output on the \scAD\ container classes, it is important to
remember that they do not perform any error checking. While a complete
discussion of error states for the stream classes is beyond the scope of this
manual, a brief discussion of the use of the operator~\texttt{!} for I/O error
checking is included. If \texttt{infile} is a stream object, then after any I/O
operation involving \texttt{infile}, \texttt{!infile} will be true if the
operation has failed and false if the operation has succeeded. For example:
\begin{lstlisting}
dvector x(1,10);
ifstream infile("data"); // associate the file data with the ifstream
                         // object infile
if (!infile) // Check if the file has been successfully opened
{
  cerr << "Error trying to open file data\n";
  exit(1); // User has decided to stop on this error condition
}

infile >> x ;            // read data into the dvector x;
if (!infile) // Check if the input operation was successful
{
  cerr << "Error trying to read dvector x from file data\n";
  exit(1); // User has decided to stop on this error condition
}
\end{lstlisting}

If you neglect to put in any error checking for I/O operations, your program may
very well carry on blissfully, even though the I/O operation has failed. The
symptoms of this failure can be quite confusing and waste a lot of valuable
program development time.

\section{Unformatted stream I/O}
\XX{input and output}{unformatted}
\XX{input and output}{binary}

It is often useful to store intermediate calculations in an external file.
Unformatted (``binary'') file I/O is the most rapid and space-efficient method
for such operations. It is supported through two additional stream classes
\texttt{uostream} and \texttt{uistream}. These two classes are derived classes
from the \texttt{ofstream} and \texttt{ifstream} classes in the
\textit{iostream} package.

\section{An example of input and output\br for \scAD\ classes}
\XX{input and output}{example}

The following listing illustrates the main features of the \scAD\ system I/O
support. It also illustrates the fill operations and \cplus\ scoping rules to
automatically destroy some variables.
\input io_examp.cpp

%xx \htmlnewfile

\chapter{Temporary Files}
\X{\fontindexentry{tt}{gradfil1.tmp}}
\X{\fontindexentry{tt}{gradfil2.tmp}}
\X{\fontindexentry{tt}{cmpdiff.tmp}}
\XX{temporary files}{\fontindexentry{tt}{gradfil1.tmp}}
\XX{temporary files}{\fontindexentry{tt}{gradfil2.tmp}}
\XX{temporary files}{\fontindexentry{tt}{cmpdiff.tmp}}
\XX{temporary files}{for saving derivative information}

The \scAD\ system saves on a stack information required to calculate the
derivatives during the evaluation of the function. The amount of information
depends on the number of variables for which derivative information is required,
and on the specific nature of the computations. It is difficult to predict
generally how much information will be generated. If the length of this stack
exceeds a preset amount, it is stored on the disk in one or two unformatted
temporary gradient files, \texttt{gradfil1.tmp} and \texttt{gradfil2.tmp}. There
is also derivative information from ``precompiled'' derivative calculations,
which is stored in the file \texttt{cmpdiff.tmp}.

\XX{\fontindexentry{sc}{ram} disk}{and temporary files}
\XX{temporary files}{use of the \textsc{dos} environment string \textsc{tmp}}
\XX{temporary files}{use of the \textsc{dos} environment string \textsc{tmp}1}
\XX{temporary files}{default directory for}
\XX{environment string}{specifying directory for temporary files}
The location of these files is controlled by the \textsc{dos} environment
strings, \texttt{TMP} and \texttt{TMP1}. The file \texttt{gradfil1.tmp} will be
created in the directory indicated by \texttt{TMP}, while the file
\texttt{cmpdiff.tmp} will be created in the directory indicated by
\texttt{TMP1}. If either or both of these strings is undefined, the current
directory will be used to create the files. Greatest speed will be achieved if
the temporary files are located on a \textsc{ram} disk. Assuming that your
\textsc{ram} disk is drive \texttt{E:}, type \texttt{set tmp=e:} at the
\textsc{dos} prompt. If \texttt{gradfil1.tmp} becomes too large for your
\textsc{ram} disk (or whatever device is specified by the \texttt{TMP}
environment string), the rest of the file will be stored in
\texttt{gradfil2.tmp} on the device referred to by the \textsc{dos} environment
string \texttt{TMP1}. For example, if you have typed \texttt{set tmp=e:} and
\texttt{set tmp1=c:\bs temp}, the first part of the gradient information will be
stored on \texttt{e:} (presumably a \textsc{ram} disk), while the second part
will be stored in the directory \texttt{c:\bs temp}. If the directory does not
exist, \scAD\ will display an error message and retire. If either the
\texttt{TMP} or \texttt{TMP1} environment string does not exist, \scAD\ will
create \texttt{gradfil1.tmp} or \texttt{gradfil2.tmp} in the current directory.

\XX{temporary files}{effect of abnormal program termination on}
The temporary gradient files are normally deleted when objects of class
\texttt{gradient\_structure} go out of scope. Abnormal termination of the
program (such as by pressing \texttt{Ctrl-Break} or use of the derivative
checker) will cause some temporary gradient files to remain in your file system.
These files may become quite large---several~Mb or more, so it is a good idea to
remove them.

%xx \htmlnewfile

\chapter{Global Variables}

\section{Adjusting the \scAD\ system global variables}

\XX{global variables}{access functions for}
\X{\fontindexentry{tt}{gradient\_structure}}
\X{controlling memory allocation}
The \scAD\ system declares a number of global variables that control memory
allocation for the \scAD\ structures. These global variables can be changed by
the use of ``access functions'' prior to declaring an object of type
\texttt{gradient\_structure}. For most applications, the default values of these
global variables should suffice. You will not need to worry about them until you
start building larger applications.

\bigskip
\XX{\fontindexentry{tt}{NUM\_RETURN\_ARRAYS}}{default value}
\begin{lstlisting}
int NUM_RETURN_ARRAYS = 10;
\end{lstlisting}
This global variable determines the maximum allowable depth of nesting of
functions that return \scAD\ variable types (see the discussion about the
\texttt{RETURN\_ARRAYS\_INCREMENT()} and \texttt{RETURN\_ARRAYS\_DECREMENT()}
instructions). The default value of 10 will suffice for almost all applications
and most users will not need to concern themselves with this variable.

To increase the number to, say, 15, you should put the instruction
\X{\fontindexentry{tt}{RETURN\_ARRAYS\_INCREMENT}}
\X{\fontindexentry{tt}{RETURN\_ARRAYS\_DECREMENT}}
\XX{error messages}{\fontindexentry{tt}%
{Overflow in RETURN\_ARRAYS stack. Increase NUM\_RETURN\_ARRAYS}}
\begin{lstlisting}
  gradient_structure::set_NUM_RETURN_ARRAYS(15);
\end{lstlisting}
into your program before an object of type \texttt{gradient\_structure} is
declared to manage the derivative calculations. If you exceed the maximum
allowable depth of nesting, the message
\begin{lstlisting}
  Overflow in RETURN_ARRAYS stack---Increase NUM_RETURN_ARRAYS
\end{lstlisting}
will appear. This message can also occur if you have put a
\texttt{RETURN\_ARRAYS\_INCREMENT} that is not matched by a
\texttt{RETURN\_ARRAYS\_DECREMENT} into one of your functions.

\bigskip
\X{\fontindexentry{tt}{GRADSTACK\_BUFFER\_SIZE}}
\begin{lstlisting}
long int GRADSTACK_BUFFER_SIZE = 2200;
\end{lstlisting}
This global variable determines the number of entries contained in the buffer
that, in turn, contains the information necessary for calculating derivatives.
For historical reasons, the actual amount of memory (in bytes) reserved for the
buffer is equal to the value of \texttt{GRADSTACK\_BUFFER\_SIZE} multiplied by
the size (in bytes) of an \scAD\ structure, \texttt{grad\_\-stack\_\-entry}. For
16-bit \textsc{dos} compilers, the size of \texttt{grad\_stack\_entry} is
28~bytes. The default value of \texttt{GRADSTACK\_BUFFER\_SIZE} is 2200 bytes,
which reserves a buffer of~$28*2200=61600$~bytes. Since for 16-bit \textsc{dos},
versions of \scAD\ buffers are limited in size to a maximum of 64K, 2200 is
about the largest size that can be used.

If you make the value smaller (say, 500), it will free up some memory (about
50K), but the program will execute more slowly, because it must store data on
the hard disk more often. So, this is a desperation move in a situation where
you must find some more memory for your program. To decrease the number to 500,
you should put the instruction
\begin{lstlisting}
  gradient_structure::set_GRAD_STACK_BUFFER_SIZE(500);
\end{lstlisting}
into your program before declaring a \texttt{gradient\_structure} object.

\bigskip
\X{\fontindexentry{tt}{CMPDIF\_BUFFER\_SIZE}}
\X{\fontindexentry{tt}{set\_CMPDIF\_BUFFER\_SIZE}}
\begin{lstlisting}
long int CMPDIF_BUFFER_SIZE = 32000L;
\end{lstlisting}
This global variable determines the size in bytes of the buffer used to contain
the information generated by the ``precompiled'' derivative code. To change this
variable, the instruction
\begin{lstlisting}
  static void gradient_structure::set_CMPDIF_BUFFER_SIZE(long int i);
\end{lstlisting}
is used.

\bigskip
\X{\fontindexentry{tt}{RETURN\_ARRAYS\_SIZE}}
\begin{lstlisting}
int RETURN_ARRAYS_SIZE = 50;
\end{lstlisting}
This global variable controls the amount of complexity that one line of
arithmetic can have. The present default value should be large enough for any
conceivable purpose. It is recommended that the user leave it alone.

\bigskip
\X{\fontindexentry{tt}{MAX\_NVAR\_OFFSET}}
\XX{independent variables}{maximum number of}
\begin{lstlisting}
int MAX_NVAR_OFFSET = 200;
\end{lstlisting}
This global variable determines the maximum number of independent variables that
can be used. In can be increased, for example, to 500, by including the
instruction
\begin{lstlisting}
  gradient_structure::set_MAX_NVAR_OFFSET(500);
\end{lstlisting}
into your program before declaring a \texttt{gradient\_structure} object. If you
have more independent variables than the value of \texttt{MAX\_NVAR\_OFFSET},
the error message
\begin{lstlisting}
  You need to increase the global variable MAX_NVAR_OFFSET to xxx
\end{lstlisting}
will appear and execution will be halted.
\XX{error messages} {\fontindexentry{tt}%
{You need to increase the global variable  MAX\_NVAR\_OFFSET to xxx}}

\bigskip
\X{\fontindexentry{tt}{ARRAY\_MEMBLOCK\_SIZE}}
\begin{lstlisting}
unsigned long ARRAY_MEMBLOCK_SIZE=100000L;
\end{lstlisting}
This global variable determines the maximum amount of memory (in bytes)
available for the \scAD\ variable type container class objects
(\texttt{dvar\_vector}, \texttt{dvar\_matrix}, etc.). It should be set to the
largest single amount of memory required when an object of type
\texttt{gradient\_structure} is declared. An \scAD\ variable array object
requires about 8 bytes for each element of the array. A \texttt{dvar\_vector}
object of size 500 will require slightly more than 4000 bytes of memory. A
\texttt{dvar\_matrix} object with 20~rows and 30~columns would require $ 20
\times 30 \times 8 = 4800$ bytes of memory. To provide the maximum amount of
memory for other purposes, you can reduce the size of
\texttt{ARRAY\_MEMBLOCK\_SIZE} until an error message telling you to increase
\texttt{ARRAY\_MEMBLOCK\_SIZE} appears.
\XX{error messages}{\fontindexentry{tt}%
{Need to increase ARRAY\_MEMBLOCK\_SIZE parameter}}

\X{\fontindexentry{tt}{MAX\_DLINKS}}
\begin{lstlisting}
unsigned MAX_DLINKS  = 1000;
\end{lstlisting}
The access function
\XX{\fontindexentry{tt}{dvariable}}%
{controlling the maximum number of with \texttt{set\_MAX\_DLINKS}}
\XX{\fontindexentry{tt}%
{set\_MAX\_DLINKS}}{setting the maximum number of dvariable objects}
\begin{lstlisting}
  gradient_structure::set_MAX_DLINKS(int i);
\end{lstlisting}
determines the maximum number of \texttt{dvariable} objects that can be in scope
at one time. The default value is 1000, which should be enough for most
purposes. If the value is exceeded, the message
\XX{error messages}{\fontindexentry{tt}%
{Need to increase the maximum number of dlinks}}
\begin{lstlisting}
  Need to increase the maximum number of dlinks
\end{lstlisting}
will appear and program execution will be terminated.

\subsubsection{Other access functions}

The complete list of access functions for the \scAD\ system global variables is:
\XX{access functions}{for global objects}
\begin{lstlisting}
static void set_NUM_RETURN_ARRAYS(int i);
static void set_ARRAY_MEMBLOCK_SIZE(unsigned long i);
static void set_GRADSTACK_BUFFER_SIZE(long int i);
static void set_MAX_NVAR_OFFSET(unsigned int i);
static void set_MAX_DLINKS(int i);
static void set_CMPDIF_BUFFER_SIZE(long int i);
\end{lstlisting}
These functions may not be used while an object of type
\texttt{gradient\_structure} is in scope. Doing so will produce an error message
and program execution will be halted.

\section{Stack size}
\XX{stack}{allocating sufficient space for the}

The size of the stack is usually controlled by a compiler-specific global
variable. The default stack size is generally too small for \scAD. It is
recommended that the stack size be increased to 15000 or 20000 bytes. The \scAD\
library routines have been compiled without the option of checking for stack
overflow, so if you do not reserve enough area for the stack, your program will
probably crash the system without any warning, or simply behave in some strange
manner.
\XX{stack}{checking for overflow of the}

Stack size is set by inserting an appropriate line in the code prior to the
\texttt{main()} block. The following examples set the stack length to 15000
bytes. Refer to your compiler documentation for more details.

\X{\fontindexentry{tt}{\_stklen}}
\X{Borland}
\XX{stack}{setting the stack size for Borland \string\cplus}
For the Borland compiler, the appropriate syntax is:
\begin{lstlisting}
extern unsigned _stklen = 15000;
main()
{
  // ... your code
}
\end{lstlisting}
\X{\fontindexentry{tt}{\_stack}}
\X{Zortech}
\XX{stack}{setting the stack size for Zortech \string\cplus}
\noindent For the Zortech compiler, the appropriate syntax is:
\begin{lstlisting}
unsigned int _stack = 15000;
main()
{
  // ... your code
}
\end{lstlisting}

%xx \htmlnewfile

\chapter{The \scAD\ Function\br Minimizing Routines}
\label{ch:function-minimizing}
\XX{minimization of functions}{quasi-Newton method}
\XX{minimization of functions}{conjugate gradient method}
\X{conjugate gradient method}
\X{quasi-Newton method}

\scAD\ is supplied with both a quasi-Newton function minimizer and a conjugate
gradient function minimizer. A quasi-Newton function minimizer is an effective
method for minimizing smooth functions of up to one or two hundred parameters.
At some point, the calculations involved in updating the inverse Hessian
approximation (on the order of $n^3$) and the memory requirements for the
inverse Hessian approximation (with $n(n+1)/2$ floating point elements), where
$n$ is the number of parameters, become prohibitive. The conjugate gradient
method, supplied with \scAD, should then be used.

\XX{minimization of functions}{user interface}
We have provided the user with two types of interfaces to the function minimizer
routines. The first approach provides the user with the maximum amount of
encapsulation of the minimization routine. It effectively allows the user to
say, ``Here is my function. Minimize it.'' The disadvantage of this approach is
some loss of flexibility in managing memory and controlling the minimization
process. To provide maximum flexibility, we also provide a lower-level interface
that allows the user to exercise complete control over the minimization process.

\X{minimize function}
The fully encapsulated minimization procedure is invoked by calling the function
\texttt{minimize}. The following code invokes the \texttt{minimize} routine to
find the minimum for the user's function \texttt{userfun}.
\input minimize.cpp

\XX{\fontindexentry{tt}{independent\_variables}}{declaration of}
\XX{\fontindexentry{tt}{independent\_variables}}{initialization of}
In this example, it is the user's responsibility to determine the number of
independent variables \texttt{nvar} and to declare the array of type
\texttt{independent\_variables} that holds the independent variables. The
elements of the \texttt{independent\_variables} array \texttt{x} are initialized
to zero. If other initial values are desired, the user must initialize the array
as well. An example of a typical simple function that could be minimized by this
code is the code for the function
\begin{equation*}
  \big(x[1]-1\big)^2 + \sum_{i=1}^{i<n} \Big(x[i+1]-x[i]\Big)^2
\end{equation*}
whose minimum is at the point $(1,1,\ldots,1)$.
\input userfun.cpp

\X{\fontindexentry{tt}{dvar\_vector}}
\X{\fontindexentry{tt}{dvariable}}
While in the main routine, the vector of independent variables is declared to be
an object of type \texttt{independent\_variables}, in \texttt{userfun}, it must
be declared to be an object of type \texttt{dvar\_vector}. This is the \scAD\
``array'' type for which derivative information is calculated. It corresponds to
an array or vector of doubles. Also, the quantities \texttt{z} and \texttt{tmp},
which correspond to objects of type \texttt{double} for which derivative
information is required, must be declared to be of type \texttt{dvariable}. The
function \texttt{value} takes a \texttt{dvariable} as its argument and returns a
\texttt{double}, the value of the \texttt{dvariable}.

This simple implementation of the function minimizer is useful for the beginning
user who doesn't wish to get too involved in the intricacies of \scAD\ or who
wants a quick solution to a minimization problem. What are its disadvantages?
The main disadvantage is that the type of arguments taken by the function
\texttt{userfun} are limited to a predefined set. In the present case, the
function \texttt{userfun} is declared to be a function that takes one argument
of type \texttt{dvar\_vector\&}. The line of code
\begin{lstlisting}
  double minimum_value= minimize(x,userfun);
\end{lstlisting}
causes a pointer to the function \texttt{userfun} to be passed to the function
\texttt{minimize}. (When a function is named without the accompanying
parentheses, a pointer to the function is generated.) Thus, the second argument
of \texttt{minimize} is a pointer to a function that takes one argument of type
\texttt{dvar\_vector\&} and returns a~\texttt{double}.

To enable the user to retain maximum control over memory management and the
minimization process, the function minimizer and the function being minimized
should be put on the same level in the program. This is accomplished by putting
them into a loop. The general form of this control structure is
\begin{lstlisting}
double f;
independent_variables x(1,nvar);
dvector g(1,nvar);
fmm fmc(nvar);  // Use this declaration for the quasi-Newton function minimizer
// fmmc fmc(nvar);  // Use this declaration for the conjugate-gradient function
                    //minimizer
gradient_structure gs;
while (fmc.ireturn >=0)
{
  fmc.fmin(f,x,g);
  if (fmc.ireturn > 0)
  {
    f=userfun(x,...) // The users function
    gradcalc(nvar,g);
  }
}
\end{lstlisting}

The \texttt{dvector} objects \texttt{g} and \texttt{x} must have minimum valid
indices equal to~1 and maximum valid indices equal to the number of independent
variables in the minimization.

\X{\fontindexentry{tt}{BEGIN\_MINIMIZATION}}
\X{\fontindexentry{tt}{END\_MINIMIZATION}}
Since many parts of the code sequence are independent of the particular
application, the code segments for calling the function minimizer and evaluating
the derivatives has been made into a macro. Using the macros,
\texttt{BEGIN\_MINIMIZATION(nvar,f,x,g,fmc)} and
\texttt{END\_MINIMIZATION(nvar,g)}, the above code can be written as
\begin{lstlisting}
double f;
independent_variables x(1,nvar);
dvector g(1,nvar);
fmm fmc(nvar); // Use this declaration for the quasi-Newton function minimizer
// fmmc fmc(nvar);  // Use this declaration for the conjugate-gradient function
                    //minimizer
BEGIN_MINIMIZATION(nvar,f,x,g,fmc)
  f=userfun(x,...) // The users function
END_MINIMIZATION(nvar,g)
\end{lstlisting}

The macro \texttt{BEGIN\_MINIMIZATION} takes five arguments: the number of
independent variables, the name of the object of type \texttt{double} that will
contain the value of the function being minimized, the independent variables,
the gradient vector, and the name of the function-minimizing control structure.
The macro \texttt{END\_MINIMIZATION} takes two arguments: the number of
independent variables and the gradient vector.
\XX{\fontindexentry{tt}{BEGIN\_MINIMIZATION}}{number of arguments}
\XX{\fontindexentry{tt}{END\_MINIMIZATION}}{number of arguments}

\XX{\fontindexentry{tt}{fmm} class}{controlling minimization with}
\XX{\fontindexentry{tt}{fmmc} class}{controlling minimization with}
Notice that even when using the macros, it is the responsibility of the user to
declare the name of the object of type \texttt{double}, \texttt{f}, which will
contain the value of the function being minimized, the
\texttt{independent\_variables x} and the \texttt{dvector g}, as well as the
\texttt{fmm} object \texttt{fmc}. The quasi-Newton function minimizer
\texttt{fmin} is controlled by, and is a member function of, the class
\texttt{fmm}. The constructor for the class \texttt{fmm} takes \texttt{nvar},
the number of variables, as an argument and creates the structure \texttt{fmc},
which controls the minimization. The user can control the operation of the
function minimizer by changing certain values of the \texttt{fmm} structure
\texttt{fmc}. The conjugate gradient function minimizer \texttt{fmin} is a
member function of the class \texttt{fmmc}. The only difference in calling the
quasi-Newton or conjugate-gradient function minimizers is in the declaration of
an object of type \texttt{fmm} for the quasi-Newton method, and of type
\texttt{fmmc} for the conjugate-gradient method. The public members of this
structure are given below.
\XX{\fontindexentry{tt}{fmm} class}{\fontindexentry{tt}{maxfn}}
\XX{\fontindexentry{tt}{fmm} class}{\fontindexentry{tt}{min\_improve}}
\XX{\fontindexentry{tt}{fmm} class}{\fontindexentry{tt}{crit}}
\XX{\fontindexentry{tt}{fmm} class}{\fontindexentry{tt}{imax}}
\XX{\fontindexentry{tt}{fmm} class}{\fontindexentry{tt}{iprint}}
\XX{\fontindexentry{tt}{fmm} class}{\fontindexentry{tt}{scroll\_flag}}
\XX{\fontindexentry{tt}{fmm} class}{\fontindexentry{tt}{ireturn}}

\begin{lstlisting}
class fmm
{
public:
                        // control variables:

  long maxfn;        // maximum number of function evaluations
                     // default value: maxfn = 500
  double min_improve;  // if the function doesn't decrease by at least this
                       // amount over 10 iterations then quit.
                       // default value: min_improve=1.e-6
  double crit;         // gradient convergence criterion; fmin terminates if
                       //   maximum absolute value of the
                       //   gradient component is less than crit.
                       //   default value: crit = 1e-4
  long imax;         // maximum number of function evaluations in linear
                     // search without improving the current estimate
                     // before fmin terminates  default value : imax = 30
  long iprint;       // number of iterations between screen displays of
                     // current parameter estimates default value: iprint = 1
  long scroll_flag;  // set to 1 to scroll display on screen;
                       //        0 to overwrite screen
                       // default value: scroll_flag = 1
  int ireturn;      // used to control the loop containing fmin and the
                     // user's function
}
\end{lstlisting}

There are five criteria by which the function minimizer may decide that
convergence has been achieved.
\X{convergence criteria for function minimizer}
\XX{function minimizer}{convergence criteria for}
\begin{enumerate}
  \item If the magnitude of all the gradient components is less than
  \texttt{crit} (default value 1.0\e{-4}).

  \item If the maximum number of function evaluations, \texttt{maxfn}, is
  exceeded (default value 500).

  \item If the function value does not decrease by more than
  \texttt{min\_improve} over 10 iterations (default value 1.0\e{-6}).

  \item If the function minimizer tries to evaluate the function more than
  \texttt{imax} times within one linear search (default value 30).

  \item If the function minimizer detects that the user has pressed the `Q' or
  `q' key.
\end{enumerate}

Here is an example of the use of these switches:
\begin{lstlisting}
dvector g(1,nvar);
dvector x(1,nvar);
double f;
fmm fmc(nvar);
fmc.maxfn=1000;       // Maximum number of function evaluations is 1000
fmc.iprint=10;        // Current best estimate will be printed out every
                      // ten iterations
fmc.crit=.01;         // Convergence if derivatives are all smaller in
                      // magnitude than .01
fmc.min_improve=0.0;  // This convergence criterion won't be used
fmc.imax=0;           // This convergence criterion won't be used
gradient_structure gs;
while (fmc.ireturn >=0)
{
  fmc.fmin(f,x,g);
  if (fmc.ireturn > 0)
  {
    f=userfun(x,...) // The users function
    gradcalc(nvar,g);
  }
}
\end{lstlisting}

The function prototype for \texttt{fmin} is:
\begin{lstlisting}[escapechar=\%]
void fmm::fmin(double& f, independent_variables & x, dvector & g);
/*--------------------------------------------------------------------------*
 * Purpose: quasi-Newton numerical minimization                             *
 *  f       : current value of function to be minimized                     *
 *  x       : current values of the variables over which the minimization   *
 *            occurs; instance of class independent variables with elements *
 *            from 1 to n                                                   *
 *  g       : current value of the gradient as returned by the function     *
 *            gradcalc( ... ); instance of class dvector with elements      *
 *            from 1 to n                                                   *
 *--------------------------------------------------------------------------*/
\end{lstlisting}

The execution of \texttt{fmin} may be terminated by pressing the `Q' key on the
keyboard. There may be a slight delay between the keypress and the termination
of \texttt{fmin}.

\section{Putting bounds on parameters}

It is often necessary to put bounds on the values that a parameter can take. For
example, if the expression $1/y$ is to be evaluated in a program, an arithmetic
error will occur if $y=0$. Similarly, $e^x$ will probably produce an error if
$x=100000$.

\X{\fontindexentry{tt}{boundp}}
\XX{bounding functions}{\fontindexentry{tt}{boundp}}
The function \texttt{boundp} is supplied with \scAD\ to simplify the placing of
bounds on parameters.

\begin{lstlisting}
dvariable boundp( dvariable x, double fmin, double fmax, dvariable& fpen)
{
  dvariable t,y;
  t=fmin + (fmax-fmin)*(sin(x*1.570795)+1)/2;
  if (x < 0.00001)
  {
    fpen+=(x-0.00001)*(x-0.00001);
  }
  if (x > 0.9999)
  {
    fpen+=(x-0.9999)*(x-0.9999);
  }
  if (x < -1)
  {
    fpen+=1000*(x+1)*(x+1);
  }
  if (x > 1)
  {
    fpen+=1000*(x-1)*(x-1);
  }
  return(t);
}
\end{lstlisting}

The function \texttt{boundp} constrains the variable to lie between the values
\texttt{fmin} and \texttt{fmax}. Both the values \texttt{fmin} and \texttt{fmax}
can actually be attained, so if you want a variable to be greater than zero, the
code
\begin{lstlisting}
  y=boundp(x,0.0,1.0,fpen);
\end{lstlisting}
will not do the job. Instead, you can accomplish this in one of two ways. You
can use a quadratic transformation if you really don't care how large \texttt{x}
can become.
\begin{lstlisting}
  y=x*x;
\end{lstlisting}
You can make a decision about what is meant by ``close to zero'' for your
application---in this case, 0.00001:
\begin{lstlisting}
  y=boundp(x,0.00001,1.0,fpen);
\end{lstlisting}

\XX{penalty function}{use wish bounding functions}
The penalty function \texttt{fpen} is included in the \texttt{boundp} routine,
because otherwise, the function would have a zero derivative at its lower and
upper bounds. This causes the derivative with respect to the bounded parameter
to be equal to zero. The result is that if the minimization routine is started
with a parameter value already at the upper or lower bound, then this parameter
will never be changed. The penalty shifts the position of the critical point
slightly within the bounds. Of course, if the initial parameter estimate happens
to coincide exactly with this critical point, the derivative will be zero, but
this is very unlikely.

If you neglect to add the penalty function to the function being minimized, the
minimization routine will still work. The only problem is that you run the risk
of one of the bounded parameters getting ``stuck'' at its upper or lower bound.

\XX{penalty function}{getting the sign correct}
The penalty should be added to the function that you wish to minimize. If you
wish to maximize the function (in which case you are going to change the sign of
the function before passing it to the function minimizer), then you should
subtract the penalty. \textit{If you get the sign of the penalty wrong, the
  routines will not work properly.}

\X{\fontindexentry{tt}{boundpin}}
\XX{bounding functions}{\fontindexentry{tt}{boundpin}}
\XX{bounding functions}{inverse function for}
To begin the minimization, you generally know the initial value of the
constrained parameter~\texttt{y}, and you need to find the value of~\texttt{x}
that will produce this value, so you can initialize the proper component of the
vector of \texttt{independent\_variables} with it. The function
\texttt{boundpin}, which is the inverse function of \texttt{boundp}, will yield
the desired value for~\texttt{x}.
\begin{lstlisting}
// double boundpin(double y, double fmin, double fmax)
x=boundpin(y,fmin,fmax);
\end{lstlisting}

\section{The derivative checker}
\XX{derivative checker}{checking correctness of derivatives with}

It may happen that the derivatives being passed to \texttt{fmin} become
incorrect for some values of the parameters. Two possible reasons for this are
1) an error in the internal derivative calculations of \scAD\ itself or 2) the
use of a user's function that is not differentiable at the point in question.
The former cause becomes increasingly less likely as experience with \scAD\
continues and any remaining bugs are removed. The latter can occur because when
constructing a large complicated function, it is always possible (and not
uncommon, in our experience) to add some non-differentiable components as an
oversight.

If the function becomes nondifferentiable at some point or the derivative
calculations are incorrect, the function minimizer will begin to experience
difficulty in making progress. To help identify derivative errors or cases of
non-differentiability when they occur, the \scAD\ system includes a interactive
derivative checker, which compares the value of the derivatives calculated by
\texttt{gradcalc()} with finite difference approximations. The derivative
checker can be invoked while the function minimizer is in progress by pressing
the `C' key on the keyboard.
\XX{derivative checker}{invoking from the function minimizer}

The derivative checker operates in two modes. In the default mode, the
derivative is checked for one independent variable at a time. In this mode, you
will be asked to supply the index of the variable for which you want to check
the derivative. Alternatively, if you want to check all derivatives, enter $0$
when prompted for the index of the variable to check. Next, you will be asked to
enter the step size to use for the finite difference approximation. Enter a
number which will alter the value of the of the variables to be checked in the
fifth or sixth significant digit. The derivative checker will then tabulate the
values of the function, the variable being checked, the analytical derivative
(from \texttt{gradcalc()}), the finite difference approximation, and the index
of the variable being checked. The derivative checker may be interrupted by
pressing the `X' key on the keyboard. On completion, the derivative checker
exits to the operating system.

For many models, the finite approximations to the derivatives may change
radically as the step size used in the derivative approximation is changed. How
then can one be sure that the analytical derivatives are correct or incorrect?
The following ``rule of thumb'' seems to work fairly well. First, identify a
gradient component that seems to be incorrect. Pick this component as the
independent variable whose derivative is to be checked. Evaluate the finite
difference approximation for different step sizes, with the step sizes differing
by a factor of about 3. Typical step sizes might be \texttt{1.e-4},
\texttt{3.e-5}, \texttt{1.e-5}, \texttt{3.e-6},$\ldots$ Now, if the finite
difference approximations stay almost the same for two consecutive step sizes,
but are very different from the analytical approximations, this is a good
indication that the analytical derivatives are incorrect.

Frequent use of the derivative checker is advisable during development of a new
application, particularly after changes that modify the objective function.

\XX{derivative checker}{invoking explicitly}
The derivative checker may also be invoked explicitly. Its function prototype is
\begin{lstlisting}[escapechar=\%]
void derch(double & f, independent_variables & x, dvector & g,
           const int n, int & ireturn);
/*--------------------------------------------------------------------------*
 * Purpose: compare analytical derivative calculation with finite           *
 *          difference approximation                                        *
 * Input:                                                                   *
 *  f       : current value of function                                     *
 *  x       : current values of the parameters of the function; instance of *
 *            class independent_variables with elements ranging from 1 to n *
 *  g       : current value of the gradient as returned by the function     *
 *            gradcalc( ... ); instance of class dvector with elements from *
 *            1 to n                                                        *
 *  n       : number of variables                                           *
 *  ireturn : return status; set ireturn = 3 prior to first call to derch() *
 * Output:                                                                  *
 *  ireturn : return status; ireturn = 4 or 5 on return from derch()        *
 *--------------------------------------------------------------------------*/
\end{lstlisting}

The following code fragment will display the values of the function, the
variable being checked, the analytical derivative (from \texttt{gradcalc()}),
the finite difference approximation, and the index of the variable being
checked.
\begin{lstlisting}
#include <math.h>
#include <fvar.hpp>

// prototype for user supplied code to compute f as a function of x
void f_comp( double & f, independent_variables & x, ... );
void main()
{
  int n = ...;        // number of variables in function
  independent_variables x(n);
  set_x( ... );      //  user supplied code to set the values of x
  dvector g(1,n);
  gradient_structure gs;
  double f = 0;
  int derch_return = 3;

  f_comp(f, x, ... );   // first call to get value of function for current x
  while (derch_return >=3)
  {
    derch(f, x, g, n, derch_return);
    {
      f_comp(f, x, ... );    // subsequent calls for with values of x
                             // varied by the step size within derch()
      gradcalc(n, g);
    }
  }
}
\end{lstlisting}

%xx \htmlnewfile

\chapter{Statistical Functions}

The principal purpose of \scAD\ is for constructing computer programs that use
mathematical derivatives. One such application is in nonlinear parameter
estimation, that is, the construction of nonlinear statistical models. It is
convenient to use simulations that produce data with known statistical
properties, in order to analyze the behavior of such models. To facilitate this
process, \scAD\ is furnished with a small number of statistical functions that
can be applied to both constant and variable vector objects.
\X{statistical functions}
\X{simulation models}
\X{nonlinear-parameter estimation}

\section{Filling vectors with random numbers}
\XX{random numbers}{filling vectors with}
\XX{filling vectors} {with random numbers}

In this section, a uniformly distributed random number is assumed to have a
uniform distribution on~$[0,1]$. A normally distributed random number is assumed
to have mean~$0$ and variance~$1$. A binomially distributed random number is
assumed to have a parameter~$p$, where~1 is returned with probability $p$, and
$0$ is returned with probability $1-p$. The \texttt{long~int\&~n} is the random
number seed specified by the user.

\X{\fontindexentry{tt}{fill\_randu}}
\X{\fontindexentry{tt}{fill\_randn}}
\X{\fontindexentry{tt}{fill\_randbi}}
\X{\fontindexentry{tt}{fill\_multinomial}}
\begin{lstlisting}
dvector::fill_randu(long int& n) // Fill a dvector from a uniform
                                 // distribution on 0,1
dvector::fill_randn(long int& n) // Fill a dvector from a standard normal
                                 // distribution
dvector::fill_randbi(int& n,double p) // Fill a dvector from a binomial
                        //distribution A 1 is returned with probability p
                        // A 0 is returned with probability 1-p
dvector::fill_multinomial(int& n,dvector p)// Fill a dvector from a
             // multinomial  distribution The distribution is determined by
            // the dvector p. An i is returned with probability p(i) if the
            // components of p sum to 1. Otherwise p is normalized so that
            // its components sum to 1.
ivector::fill_multinomial(int& n,dvector p)// Fill a dvector from a
             // multinomial distribution The distribution is determined by
            // the dvector p. An i is returned with probability p(i) if the
            // components of p sum to 1. Otherwise p is normalized so that
            // its components sum to 1.
\end{lstlisting}

The \texttt{fill\_multinomial} routine expects to receive a \texttt{dvector} of
parameters, that is, $P=(p_1,\ldots,p_n)$. The routine fills the elements of a
\texttt{dvector} with integers whose values range from $1$ to~$n$, where $i$
occurs with probability~$p_i$. If the components of $P$ do not sum to $1$, they
will be normalized so that the components do sum to~$1$.

As an example of how to use these functions, consider the following code
segment, which generates random numbers from a normal distribution contaminated
with a small proportion of large errors. This kind of distribution is useful for
studying ``robust'' parameter estimation procedures. Robust estimation
procedures are intended to be insensitive to such deviations from the main
distribution.

\XX{random number generator}{example of the use of}
\X{robust estimation}
Assume that the main distribution has mean~20 and standard deviation~2, while
the contaminating distribution has mean~20 and standard deviation~6. We want the
probability that the random numbers belong to the main distribution to be 0.90
and the probability that they belong to the contaminating distribution to
be~0.10. We shall generate a random vector of 100 such random numbers.
\begin{lstlisting}
// ...
dvector rn(1,100);
{
  dvector w(1,100);
  rn.fill_randn(212); // fills rn with numbers drawn from
                          //  a normal distribution
  w.fill_randbi(1521,0.90); // fills w with numbers drawn from
                            // a binomial distribution
  for (int i=1; i<=100;i++)
  {
    if (w(i)==1)  // This condition will occur with probability 0.90
    {
      rn(i)=2*rn(i)+20;  // standard deviation 2 and mean 20
    }
    else// This condition will occur with probability 0.10
    {
      rn(i)=6*rn(i)+20;  // standard deviation 6 and mean 20
    }
  }
} // Now the dvector w goes out of scope and the memory it used
  // is released...
\end{lstlisting}

\section{Generating a sample from a mixture\br of $n^2$ normal distributions}
\X{mixture of normal distributions}
\XX{multinomial fill}{of an ivector object}

This example extends the previous example to a mixture of more than two
distributions. Two interesting aspects of this example are the use of the
multinomial fill of an \texttt{ivector} object and the use of the resulting
\texttt{ivector} object in a general selection of the form
\texttt{dvector(ivector)}.
\X{operator \texttt{()}}
\XX{vector}{extracting a subvector}
\X{extracting a subvector}
\XX{vector}{function call \texttt{()} to extract subvector}

\begin{lstlisting}
int nsample;  // nsample is the sample size to be generated
int n;        // n is the number of groups in the mixture
dvector p(1,ngroups); // p contains the mixture proportions;
dvector mean(1,ngroups); // mean contains the mean for each group
dvector sd(1,ngroups); // sd contains the standard deviations for each group
dvector sample(1,nsample); // sample will contain the simulated data
ivector choose(1,nsample); // this will determine which group in the mixture
                           // each observation came from
//...
// somehow get the mixture proportions into p, the means into mean,
// and the standard deviations into sd
// ...
choose.fill_multinomial(1011,p);

sample.fill_randn(211); // fill the sample with standard normal deviates
sample=elem_prod(sd(choose),sample)+mean(choose); //This
            // creates the mixture of normal distributions
\end{lstlisting}

Notice that the $i^\textrm{th}$ element of \texttt{choose} will be equal to $j$
with probability $p_j$, so each element in the sample will have a probability
$p_j$ of having \texttt{mean($j$)} as its mean and \texttt{sd($j$)} as its
standard deviation. This means that each element of \texttt{sample} will have a
probability $p_j$ of belonging to the $j^\textrm{th}$ group of normal
distributions. It follows that the elements of \texttt{sample} have the desired
mixture distribution.

\section{Functions that provide summary statistics}
\XX{\fontindexentry{tt}{mean}}{mean of a vector object}
\XX{\fontindexentry{tt}{std\_dev}}{standard deviation a vector object}

The functions
\begin{lstlisting}
number mean(vector_object);
number std_dev(vector_object);
\end{lstlisting}
\noindent return the mean and standard deviation of a \texttt{dvector} or
\texttt{dvar\_vector}.

%xx \htmlnewfile

\chapter{Robust Nonlinear Regression}
\label{ch:robust-nonlinear-regression}

Many authors have indicated that the usual method of least-squares estimation of
parameters is too sensitive to the existence of a small number of large errors,
so-called outliers, in the data being analyzed. Although an extensive theory of
robust estimators that are not so sensitive to such outliers has developed, many
people still use standard statistical packages based on least-squares estimation
for statistical problems such as multiple regression. With the \scAD\ system, it
is a simple matter to use superior robust-regression estimation procedures on
both linear and nonlinear statistical models. We should emphasize that we are
not experts in the field of robust regression and we make no claims of
optimality or near optimality for these methods. The methods in this chapter
seem to work well and they have been designed so that they can be easily
extended to large nonlinear regression problems with thousands of observations,
and at least hundreds (and perhaps thousands) of parameters. (In 16-bit
\textsc{dos} versions of \scAD, the 64K barrier may limit the size of the
problem that can be considered.)
\XX{least-squares}{parameter estimation}
\XX{robust estimation}{in nonlinear parameter estimation}
\XX{parameter estimation}{least-squares}
\XX{parameter estimation}{robust-nonlinear}

Our purpose in this chapter is to present an estimation procedure that performs
almost as well as the usual least-squares estimator when the model errors are
normally distributed, and that performs much better than least squares when the
model errors contain some large outliers. With such an estimator, you don't have
to worry about whether the use of robust techniques is justified. If your model
errors are exactly normally distributed, you have lost very little, and if they
are not normally distributed, you may have gained a lot.

The reader who is not interested in the derivation of the formulas in this
section can skip to
Section~\ref{sec:using-routines}, %xx \number\mychapno.2 check this!!
where the methods for using the routines are discussed. The necessary functions
are supplied with the \scAD\ libraries.

\section{Statistical theory for robust regression}

Assume that we have $n$ data points, $Y_i$, $i=1,\ldots,n$. The $Y_i$ are
assumed to have been produced by the process
\begin{equation}
  \label{eq:robust-nonlinear-regression-00}% \eqno(\number\mychapno.0)
  Y_i=f({\mbox{\boldmath $x$}}_i;\Theta)+\epsilon_i
\end{equation}
where for each $i$, ${\mbox{\boldmath $x$}}_i$ is a known vector, $\Theta$ is an
unknown vector of parameters, and the~$\epsilon_i$ are random variables. We want
to get estimates for the parameters~$\Theta$. Some common special cases are:
\begin{description}

\item[Multiple regression:]
for each $i$,
${\mbox{\boldmath $x$}}_i =(x_{i1},\ldots,x_{ir})$ is an $r$-vector and
$\Theta =(\theta_{0},\theta_{1},\ldots,\theta_{r})$ is an $(r+1)$-vector. The
functional relationship is
\begin{equation*}
  f({\mbox{\boldmath $x$}}_i;\Theta)
  =\theta_0+x_{i1}\theta_1+x_{i2}\theta_2+\cdots+x_{ir}\theta_r
\end{equation*}

\item[Polynomial regression:]
for each $i$, ${\mbox{\boldmath $x$}}_i =(x_i)$  is a $1$-vector and
$\Theta =(\theta_{0},\theta_{1},\ldots,\theta_{r})$ is an $(r+1)$-vector. The
functional relationship is
\begin{equation*}
  f({\mbox{\boldmath $x$}}_i;\Theta)
  =\theta_0+x_i\theta_1+x_i^2\theta_2+\cdots+x_i^r\theta_r
\end{equation*}

\item[A nonlinear model:]

For a simple nonlinear model, we shall consider the case where
${\mbox{\boldmath $x$}}_i$ is a 1-dimensional vector and $\Theta$ is a
1-dimensional vector whose single component is denoted by $b$. The functional
relationship is
\begin{equation}
  \label{eq:robust-nonlinear-regression-01} % \eqno(\number\mychapno.1)
  f({\mbox{\boldmath $x$}}_i;b)=bx_i+b^2
\end{equation}
\end{description}

The form of the estimation procedure that should be used to estimate the
parameter vector~$\Theta$ depends on the assumptions that are made about the
nature of the random variables~$\epsilon_i$. Assume that the~$\epsilon_i$ are
independent identically distributed random variables with a probability density
function~$\phi(u;\Lambda)$, where~$\Lambda$ is a vector of parameters whose
components may be known or unknown.

The joint probability density function for the $(Y_1,\ldots,Y_n)$ is given by
\begin{equation}
  \label{eq:robust-nonlinear-regression-02} %\eqno(\number\mychapno.2)
  \prod_{i=1}^n \phi\big(Y_i-f({\mbox{\boldmath $x$}}_i;\Theta);\Lambda\big)
\end{equation}
Taking the logarithm of
equation~(\ref{eq:robust-nonlinear-regression-02}) %(\number\mychapno.2)
gives the expression
\begin{equation}
  \label{eq:robust-nonlinear-regression-03} %\eqno(\number\mychapno.3)
  \sum_{i=1}^n \log
  \Big(\phi\big(Y_i-f({\mbox{\boldmath $x$}}_i;\Theta);\Lambda\big)\,\Big)
\end{equation}
Considered as a function of the parameters $\Theta$ and $\Lambda$,
expression~(\ref{eq:robust-nonlinear-regression-03}) % $(\number\mychapno.3)$
is known as the log-likelihood function. The maximum likelihood estimates for
the parameters $\Theta$ and $\Lambda$ are the values of parameters that maximize
the likelihood or its equivalent, the log-likelihood function.

If we assume that the $\epsilon_i$ are normally distributed with mean~$0$ and
standard deviation~$\sigma$, then
\begin{equation*}
  \phi(u;\Lambda)= \frac{1}{\sqrt{2\pi}\sigma}
    \,\exp\left(-\frac{u^2}{2\sigma^2}\right)
\end{equation*}
and the log-likelihood function becomes
\begin{equation}
  \label{eq:robust-nonlinear-regression-04} % \eqno(\number\mychapno.4)
  -n\log(\sigma)-\frac{1}{2\sigma^2}
  \sum_{i=1}^n\big(Y_i-f({\mbox{\boldmath $x$}}_i;\Theta)\big)^2
\end{equation}

In this case, $\Lambda=(\sigma)$ is a 1-dimensional vector.

The maximum likelihood estimates $\widehat\Theta$ for the parameters $\Theta$
that maximize
equation~(\ref{eq:robust-nonlinear-regression-04}) %(\number\mychapno.4)
are the same as the those values of the parameters that minimize
\begin{equation}
  \label{eq:robust-nonlinear-regression-05} % \eqno(\number\mychapno.5)
  \sum_{i=1}^n\big(Y_i-f({\mbox{\boldmath $x$}}_i;\Theta)\big)^2
\end{equation}
Since these parameter estimates minimize the squared difference between the
predicted and observed values, they are known as the least-squares estimates.
For normally distributed model errors, the maximum likelihood estimates and the
least-squares estimates coincide. Since the
expression~(\ref{eq:robust-nonlinear-regression-05}) % (\number\mychapno.5)
does not depend on~$\sigma$, the maximum-likelihood estimate for~$\Theta$ does
not depend on the estimate for~$\sigma$. This is a special property shared by
the normal distribution and some other distributions. It is does not hold in
general. The maximum likelihood estimate~$\hat\sigma$ for the standard
deviation~$\sigma$ is given by
\begin{equation}
  \label{eq:robust-nonlinear-regression-06} % \eqno(\number\mychapno.6)
  \hat\sigma=\sqrt{\frac{1}{n}
    \sum_{i=1}^n\big(Y_i-f({\mbox{\boldmath $x$}}_i;\Theta)\big)^2}
\end{equation}
\XX{maximum likelihood estimation}{for robust regression}

The reason that the normal least-squares estimation is so sensitive to outliers
is that the normal distribution has a small ``tail,'' so the probability of
occurrence of an event drops off very quickly as one moves away from the mean a
distance of a few standard deviations. One approach to develop more robust
estimation methods is to use distributions with fatter tails. A distribution
with a fatter tail than the normal distribution's is the double exponential
distribution whose probability density function (assuming the mean is~0) is
given by:
\X{two-sided exponential distribution}
\begin{equation}
  \label{eq:robust-nonlinear-regression-07} % \eqno(\number\mychapno.7)
  \frac{1}{2\sigma}
    \exp\left|-\frac{u}{\sigma}\right|
\end{equation}
The corresponding log-likelihood function is
\begin{equation}
  \label{eq:robust-nonlinear-regression-08} % \eqno(\number\mychapno.8)
  -n\log(\sigma)-\frac{1}{\sigma}
  \sum_{i=1}^n\big|Y_i-f({\mbox{\boldmath $x$}}_i;\Theta)\big|
\end{equation}
The values of the parameters $\Theta$ that maximize
equation~(\ref{eq:robust-nonlinear-regression-08}) % (\number\mychapno.8)
are the same as the values that minimize the sum
\begin{equation}
  \label{eq:robust-nonlinear-regression-09} % \eqno(\number\mychapno.9)
 \sum_{i=1}^n\big|Y_i-f({\mbox{\boldmath $x$}}_i;\Theta)\big|
\end{equation}
Since equation~(\ref{eq:robust-nonlinear-regression-09}) % (\number\mychapno.9)
is the total absolute deviation between the predicted and observed values, these
estimates are known as the minimum total absolute deviation estimates for the
parameters.

While the minimum total absolute deviation estimates are fairly insensitive to a
small number of large model errors, they do not perform very well when the model
errors are normally distributed. To develop a maximum likelihood estimation
procedure that is less sensitive to outliers, but that also performs well when
the model errors are normally distributed, we have employed a distribution that
is a mixture of a normal distribution and another distribution with a fatter
tail. The probability density functions of the fat-tailed distribution is given
by:
\begin{equation}
  \label{eq:robust-nonlinear-regression-10} % \eqno(\number\mychapno.10)
  \frac{\sqrt{2}}{\pi\sigma e}\left\{1+\frac{x^4}{{(\sigma e)}^4}\right\}^{-1}
\end{equation}
The extra parameter $e$ has been added to the log-likelihood
function~(\ref{eq:robust-nonlinear-regression-10}) % $(\number\mychapno.10)
to adjust the ``spread'' of the fat-tailed distribution. The value of $e$ has
been set equal to~3. A mixture of two random variables with mixing proportion
$p$ is formed by choosing one of the random variables with probability~$p$ and
the other random variable with probability~$1-p$. If the probability density
functions of the mixture components are $\psi(u)$ and $\chi(u)$, then the
probability density function of the mixture is $p\psi(u)+(1-p)\chi(u)$. We shall
assume that the model errors $\epsilon_i$ are a mixture of the
distribution~(\ref{eq:robust-nonlinear-regression-10}) % (\number\mychapno.10)
with probability~0.05 and a normal distribution with probability~0.95. The
quantity $0.05$ can be interpreted as the proportion of observations
contaminated by large errors.

The corresponding log-likelihood function for the observations is
\begin{align}
  \label{eq:robust-nonlinear-regression-11} % \eqno(\number\mychapno.11)
  \nonumber -n\ln(\sigma)+\sum_{i=1}^n \log n
  & \Bigg[
  (1-p)\exp\bigg\{-\frac{\big(Y_i-f({\mbox{\boldmath $x$}}_i;\Theta)\big)^2}
  {2\sigma^2}\bigg\}\\
  & \quad + \frac{2p}{\sqrt{\pi} e}
  \bigg\{1+\frac{\big(Y_i-f({\mbox{\boldmath $x$}}_i;\Theta)\big)^4}
  {(e\sigma)^4}\bigg\}^{-1}
  \Bigg]
\end{align}

The parameter estimates $\widetilde\Theta$ and $\tilde\sigma$ for $\Theta$ and
$\sigma$, obtained by maximizing
equation~(\ref{eq:robust-nonlinear-regression-05}), % $(\number\mychapno.5)$
will be referred to as the ``robust-mixture estimator'' for the parameters. For
the robust-mixture estimator, the estimate for $\Theta$ is not independent of
the estimate for $\sigma$, so to estimate the parameters, the entire
expression~(\ref{eq:robust-nonlinear-regression-11}) % (\number\mychapno.11)
must be maximized for $\Theta$ and $\sigma$ simultaneously.

While it is true that the robust parameter estimates are the values of the
parameters that maximize equation~(\ref{eq:robust-nonlinear-regression-11}), it
is possible to reformulate the problem slightly so it is posed in terms of
parameters more easily estimated. Such reparameterizations are often useful in
nonlinear optimization. Consider the log-likelihood
function~(\ref{eq:robust-nonlinear-regression-11}). % (\number\mychapno.11)
If~$p$ is set equal to zero, the log-likelihood reduces to the normal case,
expression~(\ref{eq:robust-nonlinear-regression-04}). So, the value of the
parameter~$\sigma$ that would maximize
equation~(\ref{eq:robust-nonlinear-regression-11}) when~$p=0$ is~$\hat\sigma$,
given by
equation~(\ref{eq:robust-nonlinear-regression-06}). % (\number\mychapno.6)

Define a new parameter $\alpha$ by $\alpha=\sigma/\hat\sigma$ such that
$\sigma=\alpha\hat\sigma$. Replace the parameter~$\sigma$ in
equation~(\ref{eq:robust-nonlinear-regression-11}) % (\number\mychapno.11)
by the parameter~$\alpha\hat\sigma$. The log-likelihood function now takes the
form
\begin{align}
  \label{eq:robust-nonlinear-regression-12} % \eqno(\number\mychapno.12)
  \nonumber
  n\ln(\alpha\hat\sigma)-\sum_{i=1}^n \log
  & \Bigg[
  (1-p)\exp\bigg\{-\frac{\big(Y_i-f({\mbox{\boldmath $x$}}_i;\Theta)\big)^2}
  {2(\alpha\hat\sigma)^2}\bigg\}\\
  &\quad + \frac{2p}{\sqrt{\pi}e}
  \bigg\{1+\frac{\big(Y_i-f({\mbox{\boldmath $x$}}_i;\Theta)\big)^4}
  {(e\alpha\hat\sigma)^4}\bigg\}^{-1}
  \Bigg]
\end{align}
The parameter~$\alpha$ can be interpreted as an estimate of how much of the
residuals in the model fit are contributed by the large outliers generated by
the contaminating distribution. If~$\alpha$ is close to~$1$, most of the
contribution to the residuals comes from the ``ordinary'' errors associated with
the normal distribution. If~$\alpha$ is close to~$0$, most of the contribution
to the residuals comes from the large outliers associated with the fat-tailed
distribution. We expect that the maximum likelihood estimate~$\tilde \alpha$
for~$\alpha$ will lie between~$0$ and~$1$. The parameter~$\alpha$ can also be
interpreted as a ``cutoff'' switch. As~$\alpha$ is made smaller, the influence
of the larger residuals on the parameter estimates will be reduced. This is much
like looking at the data and removing the data points that fit the model badly.
\X{outliers}{in robust regression}

The code for the log-likelihood
function~(\ref{eq:robust-nonlinear-regression-12}) % \number\mychapno.12
can be written down very concisely using \scAD's vector operations. Suppose that
the observations $Y_i$ are contained in a \texttt{dvector} called \texttt{obs},
and the predicted values $f({\mbox{\boldmath $x$}}_i;\Theta)$ are contained in a
\texttt{dvar\_vector} called \texttt{pred}. The vector of residuals $Y_i
-f({\mbox{\boldmath $x$}}_i;\Theta)$ and the vector of squared residuals can be
calculated in two statements:
\begin{lstlisting}
dvar_vector diff = obs-pred;     // These are the residuals
dvar_vector diff2 = pow(diff,2); // These are the squared residuals
v_hat = mean(diff2);
\end{lstlisting}
The estimate of the variance is given by the mean of the vector of squared
residuals. If $\texttt{b}=2p/\sqrt{\pi}e$ and $\texttt{a2}=2\alpha$, the code
for equation~(\ref{eq:robust-nonlinear-regression-12}) % \number\mychapno.12
can be written as
\begin{lstlisting}
  dvariable log_likelihood = 0.5*diff.size()*log(a2*v_hat)
    -sum(log((1.-pcon)*exp(-diff2/(2*a2*v_hat))
    + b/(1.+pow(diff2/(a2*v_hat),2))));
\end{lstlisting}
\XX{vector operations}{example of}
\XX{vector operations}{to reduce the amount of temporary storage required}

While the combination of scalars and vectors in statements like \texttt{1.+v}
can seem a bit strange at first, one quickly gets used ot it. The main thing to
consider is that some combinations of operations are not associative, so it may
be necessary to include parentheses to ensure that the operations are carried
out in the correct order. As is often the case, using vector operations instead
of writing the code in terms of the components of the vector objects will
greatly reduce the amount of temporary storage required.

As the examples below illustrate, the parameter $\alpha$ is often not
well-determined by the data. Rather than estimating this parameter directly, we
feel that the best way to do robust regression is probably to use a fixed value
for $\alpha$. We recommend the value~$0.7$ for fitting the model, although we
have no strong reason for preferring this particular value.

\section{Using the \scAD\ robust\br nonlinear regression routines}
\label{sec:using-routines}

The \scAD\ system permits the user to obtain these robust regression estimates
in a simple fashion. Suppose that the observations $Y_i$ are contained in a
\texttt{dvector} called \texttt{OBS\_Y} and that the values predicted by the
model for these observations have been computed and are contained in a
\texttt{dvar\_array} called \texttt{PRED\_Y}. To compute the value of the
log-likelihood
function~(\ref{eq:robust-nonlinear-regression-05}) % (\number\mychapno.5)
with a fixed value of $\alpha$, the routine \texttt{robust\_regression} is used.
The function prototype for this routine is:
\begin{lstlisting}
  dvariable robust_regression(dvector& OBS_Y,dvar_vector& PRED_Y,double& alpha);
\end{lstlisting}
The function \texttt{robust\_regression} returns minus the value of the
log-likelihood function. (The sign has been changed for the function minimizer.)
\begin{figure}
  \centering\hskip1pt
   \input graph.tex
   \emptycaption{}
   \label{fig:robust-nonlinear-regression-01} %  \number\mychapno.1
\end{figure}

To estimate the parameter $\alpha$ as well, the parameter must be passed as a
\texttt{dvariable}. The function prototype is
\begin{lstlisting}[escapechar=\%]
  dvariable %
  %robust_regression(dvector& OBS_Y,dvar_vector& PRED_Y,dvariable& alpha);
\end{lstlisting}

\XX{bounding functions}{\fontindexentry{tt}{boundp}}
\XX{bounding functions}{\fontindexentry{tt}{boundpin}}
For the routine \texttt{robust\_regression} to work properly when the parameter
$\alpha$ is being estimated, it is necessary to put bounds on the values that
$\alpha$ can assume. This is done by using the \texttt{boundp} function supplied
with \scAD. The parameter $\alpha$ must be restricted such that $0<a<=1$.
(Actually, the upper bound is not necessary, but it will make the routine run a
bit more efficiently.) The code
\begin{lstlisting}
  a=boundp(x(i),.001,1.0,zpen);
\end{lstlisting}
restricts the value of $\alpha$ to lie within $0.001$ and $1.0$. To find the
initial value that will give a desired initial value for~$\alpha$, the inverse
function \texttt{boundpin} can be used. To obtain the initial value that gives a
value of $0.7$ for~$\alpha$, and to put the initial value in the $i$th component
of the vector of independent variables \texttt{x}, the code
\begin{lstlisting}
  x(i)=boundpin(0.7,001,1.00);
\end{lstlisting}
would be employed.

\section{A simple example of robust regression}

\XX{robust regression}{a simple example} We shall fit a simple 1-parameter model
using the robust regression routine \texttt{robust\_regres\-sion}, which
estimates the value of the parameter $\alpha$. The form of the model is
\begin{equation*}
  Y_i=bZ_i+b^2
\end{equation*}
where the true value of~$b$ is~5. The data for this example were produced by
simulating 18 real data points with normally distributed errors with mean $0$
and standard deviation 40, then adding two outliers that were 3.1 standard
deviations and 6 standard deviations removed from the true value. Of course,
such large outliers could be identified by eye in this simple example. For large
nonlinear models, however, there is often no simple graphical representation of
the fit that can enable such a simple inspection of the residuals. A small
circle ($\circ$) was used to mark the outliers. The results of the fitting
procedure are shown in
Figure~\ref{fig:robust-nonlinear-regression-01} % \number\mychapno.1.
The least-squares fit with the outliers removed (estimated $b=4.45$) is the
closest to the true relationship. It represents the best estimates that can be
obtained if the outliers are identified with certainty. The robust fit including
the outliers (estimated $b=4.06$) did considerably better than did the
least-squares fit including the outliers (estimated~$b=3.51$).

For the robust fit, the estimate of $\hat\sigma$ is $70.1$, while the estimate
of $\tilde\sigma$ is $55.6$. So, the estimate of $\alpha$ is $0.79$. Since the
real value of~$\sigma$ was $40.0$, the robust regression routine has
underestimated the contribution to the residuals made by the outliers.

The \scAD\ code for this model is given below:
\input robust.cpp

\input fcomp_r.cpp

\input  newreg2.cpp

\section{Setting the cutoff point in the\br robust regression routines}
\label{sec:setting-cutoff}

It is important to be able to specify the ratio $\tilde\sigma/\hat\sigma$,
because this parameter is often not well determined by the data. In the previous
example, the robust estimation procedure has somewhat underestimated the
contribution to the overall variance of the residuals made by the outliers. Now
consider the analysis for the same data as above---with one change. The value of
the second outlier has been set equal to $-75.0$, rather than $-125.0$. Although
this represents an improvement in the quality of the data set, the parameter
estimates (estimated $b=3.80$) for the robust-mixture estimation procedure when
the parameter $\alpha$ is estimated are actually worse than before. The reason
for this is that the estimation of the ratio $\alpha=\tilde\sigma/\hat\sigma$
has increased to $0.94$, which means that the robust estimation procedure no
longer recognizes that there are outliers in the data. The robust-mixture
estimates to the data for which the ratio $\tilde\sigma/\hat\sigma$ has been
fixed at $0.7$ (estimated $b=4.20$) is included. This is considerably better
than the least-squares fit (estimated $b=3.70$).

\begin{figure}
  \centering\hskip1pt
   \input graph1.tex
   \emptycaption{}
   \label{fig:robust-nonlinear-regression-02}  %\number\mychapno.2
\end{figure}

We suspect that a little ``robustness'' is much better than none at all, and
that overall best performance of the model may be obtained by simply fixing the
value for the ratio~$\tilde\sigma/\hat\sigma$ at some nominal value, such
as~$0.7$.

\section{A Monte Carlo evaluation\br of the robust regression routine}
\label{sec:monte-carlo}

While the above example is indicative of the superiority of the robust-mixture
estimator, it only indicates how the estimators perform in one special case.
From a statistical point of view, one is interested in the average performance
of the estimators for many different realizations of the data. This section
describes the results of such an investigation. The data were fit by the
least-squares estimator, the minimum total absolute deviation estimator, and the
robust-mixture estimator (with fixed $a=0.7$).

For each simulation, a set of 30 pairs $(Y_i,z_i)$ were generated, where
$z_i=i$, \mbox{for $i=1,\ldots,30$}, and
\begin{equation*}
  Y_i=bz_i+b^2+\epsilon_i
\end{equation*}
The parameter $b$ was set equal to~5. For the case of normally distributed model
errors, the~$\epsilon_i$ are normally distributed random variables with standard
deviation~40. For model errors exhibiting moderate deviations from normality,
the $\epsilon_i$ were a mixture of normally distributed random variables with
standard deviation 40, with probability 0.90, and normally distributed random
variables with standard deviation 120, with probability~0.10. For model errors
exhibiting large deviations from normality, $\epsilon_i$ were a mixture of
normally distributed random variables with standard deviation 40, with
probability 0.90, and normally distributed random variables with standard
deviation 200, with probability~0.10.

Five hundred replications were made for each of the error distributions. For all
the simulations, the value of ratio $\alpha$ was kept fixed at $0.7$ for the
robust-mixture estimator.

\begin{wrapfigure}{r}{.45\textwidth}%
  \centering
  % {\includegraphics[height=2in]{foo.eps}}%
  \vspace*{-.75\baselineskip}%  Had to do this to fit on page.
  {\input robgrph2.tex \par}%
  \caption{Plot of Robust-mixture estimates absolute error vs. Least-squares
    estimates absolute error.}%
  \label{fig:monte-carlo-03}
\end{wrapfigure}

Now it is known that for normally distributed errors, the usual least-squares
estimates are the minimum variance unbiased estimates of the model parameters.
It follows that if the robust estimates are unbiased, they must have greater
variance than the least-squares estimates. This greater variance is the penalty
we must pay for using robust estimates in the case where they are not necessary.
For the problem studied here, however, the standard deviation of the robust
estimates was only about 8\% greater than the least-squares estimates for the
case of normally distributed errors, while the robust-mixture estimator had a
standard deviation about 20\% smaller than that of the least-squares estimates
for moderate deviations from normality, and about 35\% smaller for large
deviations from normality.

\begin{table}[!h]
  % Recursive: three sets of three columns.
  \begin{tabular}{ l   *{3}{ c *{3}{r@{.}l } }     }
    \hline\\[-11pt]
    \hline\\[-9.0pt]
              && \multicolumn{6}{c}{Normal errors}
              && \multicolumn{6}{c}{Moderate deviations}
              && \multicolumn{6}{c}{Large deviations}  %&\multicolumn{1}{c}{}
    \\[3pt]
    \cline{3-8}\cline{10-15}\cline{17-22} \\[-6pt]
    Mean
    &&\multicolumn{2}{l}{LS} & \multicolumn{2}{l}{RM} & \multicolumn{2}{l}{MD}
    &&\multicolumn{2}{l}{LS} & \multicolumn{2}{l}{RM} & \multicolumn{2}{l}{MD}
    &&\multicolumn{2}{l}{LS} & \multicolumn{2}{l}{RM} & \multicolumn{2}{l}{MD}\\
    Mean
    & \multicolumn{1}{c}{\qquad} & 5&01 & 5&01 & 5&00
    & \multicolumn{1}{c}{\qquad} & 4&99 & 5&01 & 5&01
    & \multicolumn{1}{c}{\qquad} & 4&98 & 4&99 & 5&01\\
    Std Error
    && 0&012 & 0&013 & 0&015 && 0&016 & 0&014 & 0&016 && 0&023 & 0&015 & 0&016\\
    Std Dev
    && 0&26  & 0&28  & 0&33  && 0&36  & 0&30  & 0&36  &&  0&51 & 0&33  & 0&36\\
    \multicolumn{3}{l}{\% improvement}
    & \multicolumn{1}{c}{} & -0&077 & -0&27
    &&\multicolumn{2}{c}{} & 20&0   &  0&00
    &&\multicolumn{2}{c}{} & 35&3   &  29&4\\[3pt]
    \hline
  \end{tabular}
  \caption{Comparison of the least squares (LS),  minimum total absolute
    deviation (MD), and robust-mixture (RM) estimates for $b$ in the model
    \mbox{$Y=bz+b^2$.} The true value was $b=5.0$}
  \label{tab:robust-nonlinear-regression-01}
\end{table}

All three estimators appear to be unbiased. The standard deviation of the
least-squares estimates almost doubles as we go from the normal errors to errors
with large deviations from normality. The robust-mixture estimator exhibits much
greater stability. The minimum total absolute deviation estimator performs
almost as well as the robust-mixture estimator when there are large outliers
present, but does not perform very well for normally distributed errors or for
medium-sized outliers present.

Another way to compare the performance of the robust-mixture estimator over the
least-squares estimator is illustrated by Figure~\ref{fig:monte-carlo-03}.
% {fig:robust-nonlinear-regression-02}. % \number\mychapno.2.
The absolute value of the difference between the robust-mixture estimator
for $b$ and the actual value of $b$ is plotted against the
absolute value of the difference between the least-squares estimator for
$b$ and the actual value of $b$. For points that lie below the
diagonal line, the robust-mixture estimator has smaller absolute error
 than the least-squares estimator.
The plot is for the large-deviations set of simulations.

\section{Investigating extremely bad data\br with the robust-mixture
estimator}
\label{sec:investigating-bad-data}

While the errors in the data in the previous examples exhibited moderate and
large deviations from normality, the data were not really ``bad.'' A much more
difficult problem for parameter estimation is created if all the large errors
are concentrated on one side of the true line. To compare the estimation methods
for this kind of ``bad'' data, we generated a data set where the majority (30
data points) of the data are distributed around the line $y=x$, with the
outliers (20 data points) distributed around the point (6,2). The data were fit
with the 2-parameter linear model
\begin{equation*}
  Y=a+bx
\end{equation*}
where the true values were $a=0.0$ and $b=1.0$.

Various fits to the data are shown in
Figure~\ref{fig:robust-nonlinear-regression-03}. %\number\mychapno.3.
When the value of $\alpha$ is set equal to~0.1, the estimation procedure
completely ignores the 20 outliers. We anticipate that varying the value of
$\alpha$, noting whether the residuals suddenly shift, could be a useful way to
detect these types of outliers in large nonlinear models. Of course, in a large
model, it is not always possible to view the results in this simple way. One way
to examine how the vector of residuals is changing is to calculate the Euclidean
distance between the vectors. If the residuals are changing slowly, the distance
between the residuals should be small. A sudden shift in the residuals should
manifest itself in a large distance between the vectors.
\begin{figure}
  \centering\hskip1pt
   \input grph3.tex
   \emptycaption{}
   \label{fig:robust-nonlinear-regression-03}
\end{figure}

Table~\ref{tab:robust-nonlinear-regression-02} % \number\mychapno.2
shows the distance between the vectors of residuals for different values of
$\alpha$. It is clear that there are two groups, one group for values
$0.9>\alpha>0.4$ and another group for $0.3>\alpha>0.1$. Constructing such a
table could be a possible way to search for ``interesting'' structure in the
outliers in a large nonlinear model.

\begin{table}[!h]
  \begin{tabular}{c @{\qquad}*{9}{r}}
    \hline\\[-11pt]
    \hline\\[-9.0pt]
    $\alpha$&0.9&0.8&0.7&0.6&0.5&0.4&0.3&0.2&0.1 \\
    \cline{2-10}\\[-9.0pt]
    0.9&0.0&\\%&0.3&1.0&1.9&2.7&3.4&11.6&14.0&15.9\cr
    0.8&0.3&0.0&\\%0.6&1.6&2.4&3.1&11.9&14.3&16.2\cr
    0.7&1.0&0.6&0.0&\\%0.9&1.7&2.4&12.5&14.9&16.9\cr
    0.6&1.9&1.6&0.9&0.0&\\%0.7&1.5&13.4&15.8&17.7\cr
    0.5&2.7&2.4&1.7&0.7&0.0&\\%0.7&14.1&16.5&18.4\cr
    0.4&3.4&3.1&2.4&1.5&0.7&0.0&\\%14.7&17.1&19.0\cr
    0.3&11.6&11.9&12.5&13.4&14.1&14.7&0.0&\\%2.4&4.3\cr
    0.2&14.0&14.3&14.9&15.8&16.5&17.1&2.4&0.0&\\%1.9\cr
    0.1&15.9&16.2&16.9&17.7&18.4&19.0&4.3&1.9&0.0\\
    %0.9&0.000&0.311&1.007&1.940&2.733&3.478&11.614&14.041&15.985\cr
    %0.&0.311&0.000&0.696&1.629&2.423&3.171&11.905&14.332&16.277\cr
    %0.&1.007&0.696&0.000&0.934&1.729&2.487&12.557&14.984&16.927\cr
    %0.&1.940&1.629&0.934&0.000&0.799&1.574&13.424&15.848&17.790\cr
    %0.&2.733&2.423&1.729&0.799&0.000&0.790&14.132&16.553&18.492\cr
    %0.&3.478&3.171&2.487&1.574&0.790&0.000&14.723&17.136&19.070\cr
    %0.&11.614&11.905&12.557&13.424&14.132&14.723&0.000&2.427&4.371\cr
    %0.&14.041&14.332&14.984&15.848&16.553&17.136&2.427&0.000&1.944\cr
    %0.&15.985&16.277&16.927&17.790&18.492&19.070&4.371&1.944&0.000\cr
    \hline
  \end{tabular}
  \caption{Euclidean distance between the vectors of residuals for the
    robust-mixture estimator for different values of $\alpha$.}
  \label{tab:robust-nonlinear-regression-02}
\end{table}

%xx \htmlnewfile

\chapter{Problems with More than\br One Dependent Variable}
\label{ch:problems-more-dependent-variable}
\X{several dependent variables}
\XX{several dependent variables}{use of the operator \texttt{<{}<}}

\section{Using more than one dependent variable}

'Till now, all the problems considered have had only one dependent variable.
Many problems of interest, such as the problem of solving a system of nonlinear
equations, require the computation of derivatives with respect to more than one
dependent variable.

The operator \texttt{<{}<} has been overloaded to implement more than one
dependent variable. The user simply writes something like
\begin{lstlisting}
  u << -exp(x(1))+y; // u is a dependent variable
\end{lstlisting}
where \texttt{u} is a \texttt{dvariable} and will become a dependent variable
for which derivatives can be calculated. The operator \texttt{<{}<} acts like
the equals sign (\texttt{=}), with the additional property that \texttt{u} is
made a dependent variable. Of course, it is possible to use \texttt{u} in later
calculations, such~as
\begin{lstlisting}
u << -exp(x(1))+y; // u is a dependent variable
// ...
v << u*u; // v is a dependent variable
\end{lstlisting}

The process of declaring dependent variables in this way is global---that is, it
does not have to be done in one function. The dependent variables are ordered in
the order in which they are encountered when the code is executed.

The default maximum number of dependent variables (which is at least 10) depends
on the implementation of \scAD\ you are using. If desired, it can be increased
by using the static member function
\XX{several dependent variables}{default number of}
\XX{dependent variables}{\fontindexentry{tt}%
{set\_NUM\_DEPENDENT\_VARIABLES(int i);}}
\X{\fontindexentry{tt}{set\_NUM\_DEPENDENT\_VARIABLES(int i);}}
\begin{lstlisting}
  gradient_structure::set_NUM_DEPENDENT_VARIABLES(int i);
\end{lstlisting}

For example, to set the maximum number of dependent variables to~25, you would
put the lines
\begin{lstlisting}
gradient_structure::set_NUM_DEPENDENT_VARIABLES(25);
// ...  must occur before the gradient_structure declaration
gradient_structure gs;
\end{lstlisting}
into your code. The derivatives are calculated by calling the function
\begin{lstlisting}
  jacobcalc(int\& nvar,dmatrix\& jac)
\end{lstlisting}

The derivatives are put into the \texttt{dmatrix} called \texttt{jac}. The
minimum valid row index of \texttt{jac} must be~1, while the maximum valid row
index must be greater than or equal to the number of dependent variables. Each
row of \texttt{jac} must have a minimum valid index of~1 and a maximum valid
index greater than or equal to \texttt{nvar}, where \texttt{nvar} is the number
of independent variables.

\X{\fontindexentry{tt}{gradient\_structure}}
\X{\fontindexentry{tt}{jacobcalc}}
\XX{finding roots of equations}{Newton-Raphson method}
\XX{Newton-Raphson method}{roots of equations}
\XX{Newton-Raphson method}{example of}

\section {Finding roots of systems of\br equations with Newton-Raphson}
\label{sec:finding-roots}

Suppose that you wish to solve the nonlinear system of equations
\begin{align}\label{eq:problems-more-dependent-variable-01}
  \nonumber %\eqno(\number\mychapno.1)
f_1(x_1,\ldots,x_n)&=0\\ \nonumber
f_2(x_1,\ldots,x_n)&=0\\
f_2(x_1,\ldots,x_n)&=0\\  \nonumber
& \,\,\,\vdots\\ \nonumber
f_n(x_1,\ldots,x_n)&=0\\ \nonumber
\end{align}
In general, this is a very difficult problem. However, under suitable regularity
conditions, if you can find an initial estimate for $(x_1,\ldots,x_n)$ that is
close enough, then the iterative Newton-Raphson technique will efficiently
converge to the solution. Let $X$ denote the vector $(x_1,\ldots,x_n)$. Let
\begin{equation*}
  J(X)=\frac{\partial f_i(X)}{\partial x_j}
\end{equation*}
be the Jacobian matrix. Then to first order,
\begin{equation}
  \label{eq:problems-more-dependent-variable-02} % \eqno(\number\mychapno.2)
  f(X+\delta)=f(X)+J(X)\delta,
\end{equation}
where $f(X)$ denotes the vector $f_1(X),\ldots,f_n(X)$ and $J(X)\delta$ denotes
the multiplication of the vector $\delta$ by the matrix $J(X)$. Equating
(\ref{eq:problems-more-dependent-variable-02}) %(\number\mychapno.2)
to zero and solving for $\delta$ yields
\begin{equation}
  \label{eq:problems-more-dependent-variable-03} %\eqno(\number\mychapno.3)
  \delta=-J(X)^{-1}f(X)
\end{equation}
Now we proceed iteratively. If $X^{(r)}$ is the current estimate for $X$, let
\begin{equation*}
  X^{(r+1)}=X^{(r)} +\delta
\end{equation*}
be the next estimate where $\delta$ is defined
by~(\ref{eq:problems-more-dependent-variable-03}). % (\number\mychapno.3).

\section{Implementation of the Newton-Raphson technique}

This example shows the \nr\ technique for two equations.
\input newt.cpp

We have chosen a rather ``degenerate'' example where the Jacobian matrix is
almost singular at the solution. The two functions are
\begin{align*}
  f_1(x_1,x_2)&=100x_1^2+x_2^2-1 \\[6pt]
  f_2(x_1,x_2)&=100x_1^2+x_2^2-1+.0001-.0001x_1^4+.01x_2^2
\end{align*}

The calculation of the Jacobian and the Newton-Raphson ``update'' step can be
accomplished with two lines of code:
\begin{lstlisting}
jacobcalc(nvar,jacobian)
x-=inv(jacobian)*f
\end{lstlisting}
The iteration continues until the norm of the \texttt{dvector} called \texttt{f}
is less than 1.0\e{-12}.

The only variable types occur in the user's function \texttt{function}. All the
arithmetic for the Newton-Raphson update is done with vectors and matrices that
are of constant type. The \texttt{value} function is used to turn the
\texttt{dvar\_vector} called \texttt{f} in \texttt{function} into a
\texttt{dvector} object, so it can be returned to the calling routine.

%xx \htmlnewfile

\chapter{A Complete Nonlinear Parameter Estimation Program}
\label{nonlinear-estimation}

A major difference between linear and nonlinear statistical modeling is the
inherent instability of nonlinear models. It seems to be an empirical fact that
for any reasonably complicated model, if you simply write correct naive computer
code for estimating the parameters, the program will not work. By ``work,'' we
mean you need to get the information you want out of the data you are trying to
analyze. In general, there are two reasons for this behavior:
\begin{enumerate}
  \item The function being calculated is numerically unstable, so a zero divide
  or a floating point overflow will occur. For example, calculate~$1/x$ when~$x$
  has the value~$0$, or calculate~$\exp(x)$ when~$x$ has the value~$90000$.

  \item The ``correct'' values of the parameter estimates occur at a local
  minimum of the objective function, not at a global minimum.
\end{enumerate}

The problem we consider in this chapter exhibits both kinds of pathology
mentioned above. We will describe the techniques necessary for producing a
stable, controllable nonlinear parameter estimation routine for solving this
problem. The code for this example is included on the \scAD\ distribution disk.

\section{Finite mixture problems}
\label{sec:finite-mixture}

Consider a set of~$m$~random variables~$X_j$, indexed by~$j$, with probability
density functions~$F_j(x)$. A new random variable~$Y$ can be created from
the~$X_j$ as follows. A realization of~$Y$ is produced by first picking one of
the~$X_j$ at random, with probability~$p_j$, and then taking a realization of
the resulting~$X_j$. The random variable~$Y$ is a finite mixture of the~$X_j$,
with mixing proportions $p_j$. The probability density function for~$Y_j$ is
given by~$\sum_{j=1}^m p_jF_j(x)$. The particular finite mixture that we shall
address is a mixture of normal distributions. Effective techniques for
estimating parameters for this problem were first considered by
Hasselblad~\cite{hasselblad1966}.

Assume that $X_j$ is normally distributed with mean~$\mu_j$ and standard
deviation~$\sigma_j$. The probability density function for~$X_j$ is given by
\begin{equation*}
  \frac{1}{\sqrt{2\pi}\sigma_j}\,
  \exp\left(\, \frac{-(x-\mu_j)^2}{2\sigma_j^2}\right)
\end{equation*}
while the probability density function for the mixture~$Y$ is given by:
\begin{equation*}
  \sum_{j=1}^m \frac{p_j}{\sqrt{2\pi}\sigma_j}
                        \,\exp\left(\frac{-(x-\mu_j)^2}{2\sigma_j^2}\right)
\end{equation*}
Now suppose that we are given a vector $(y_1,\ldots,y_n)$ of $n$ realizations of
the random variable~$Y$. We wish to use the information in the~$y_i$ (and any
other auxiliary information that may be at our disposal) to estimate the
parameters~$p_j$,~$\mu_j$, and~$\sigma_j$.

In general, there are many kinds of auxiliary information that may be present in
such mixture problems. For example, we may or may not know the number of groups
present in the mixture. We may know some or all of the means,~$\mu_j$, or we may
have some knowledge of the approximate value of some or all of the~$\mu_j$. We
may know the value of some or all of the~$\sigma_j$, or we may know that all
the~$\sigma_j$ have the same value. A general routine for estimating the
parameters in mixture problems should be able, as much as possible, to
incorporate auxiliary information about the parameters into the estimation
process.

The maximum likelihood estimates for the parameters are found by finding the
correct local maximum of the log-likelihood function
\begin{equation}\label{eq:nonlinear-estimation-01} % \eqno(\number\mychapno.1
  \max_{p,\mu,\sigma}
   \Bigg\{
       \sum_{i=1}^n \,\log
         \Bigg[\,
             \sum_{j=1}^m \frac{p_j}{\sqrt{2\pi}\sigma_j}
               \,\exp
                   \left(
                       \frac{-(y_i-\mu_j)^2}{2\sigma_j^2}
                   \right)\
          \Bigg]
    \Bigg\}
\end{equation}
subject to any constraints on the parameters that may be appropriate for the
problem.

For mixtures of normal distributions, the correct solution never occurs at the
global maximum of the log-likelihood function when there is more than one group
in the mixture. This is due to the fact that it is always possible to make the
log-likelihood function as large as desired, by setting the mean of one of the
groups equal to the value of an observation and letting the corresponding
standard deviation parameter tend to zero. The value of the probability density
function for this group will then tend to infinity at the point in question. It
follows that to find the correct solution, the values of the standard deviations
should be bounded away from zero.

The following code calculates the likelihood function corresponding to
equation~(\ref{eq:nonlinear-estimation-01}). Notice that none of the constraints
that we intend to place on permissible parameter values appears in the code. The
constraint conditions will be implemented elsewhere. Avoiding them at this stage
produces a clean piece of code that can easily be seen to correspond to
equation~(\ref{eq:nonlinear-estimation-01}). The constant terms $1/\sqrt{2\pi}$
have been ignored.

\input like.cpp
The term \texttt{1.e-20} is put in to improve the numerical stability of the
program. The routine \texttt{log\_likelihood} is intended to be used with a
function-minimizing routine to estimate the parameters. In the course of such a
process, the function minimizer invariably will pick some parameters that fit
certain observations extremely badly. An extremely bad fit is characterized by
the value of \texttt{diff} being large for all groups. As a result, due to
floating point underflow,
$\exp(-\texttt{diff}*\texttt{diff}/(2.*\texttt{v(j)}))$ can be equal to zero for
all groups, so \texttt{sum} has the value zero. Taking the logarithm of~$0$
would cause a floating point exception or error. Adding the small value
\texttt{1.e-20} allows the program to carry on gracefully.

It should be noted that Hasselblad grouped the observed data into size classes
and analyzed the grouped data. This can easily be done---and should be done---if
the number of observations is large, because the analysis of the grouped data
can be much faster.

\section{Putting bounds on parameter values}
\X{\fontindexentry{tt}{boundp}}
\XX{bounding functions}{\fontindexentry{tt}{boundp}}
\X{penalty function}

The function minimizer routine deals with a vector \texttt{x} of numbers. The
likelihood calculation routine \texttt{log\_likelihood} deals with vectors of
proportions, means, and standard deviations. We need a routine to ``put'' the
$x$s into the model.
\input reset.cpp

\X{\fontindexentry{tt}{boundp}}
The function \texttt{reset} not only puts the $x$s into the parameters of the
function \texttt{like\_fun}, it also determines what parameters are active
(being estimated at the time), by examining the \texttt{ivector control}, and
puts constraints on the values of the parameters through the function
\texttt{boundp}. The function
\begin{lstlisting}
  dvariable boundp(dvariable x,double fmin,double fmax,dvariable& fpen)
\end{lstlisting}
\X{penalty function}
bounds the variable \texttt{x} to the closed interval
$[$\texttt{fmin, fmax}$]$.
The term \texttt{fpen} is a penalty term that modifies the behavior of the
function. It must be added to the function being minimized. \textit{Note that
  this means that if you want to maximize a function, you must change the sign
  of the function before adding} \texttt{fpen} \textit{to it!} The function
\texttt{reset} supplies complete control over what parameters are active,
through the vector \texttt{control}. For example, if \texttt{control(1)} is
equal to~1 while \texttt{control(2)} and \texttt{control(3)} are equal to~0,
then only the proportions at age will be estimated, while the means and standard
deviations will be kept fixed. The \texttt{dvector}s called \texttt{mumin} and
\texttt{mumax} determine lower and upper bounds on the means, while the
\texttt{dvector}s \texttt{sdmin} and \texttt{sdmax} determine lower and upper
bounds on the standard deviations.

Suppose that you want the standard deviations to be the same for all the groups.
This can be accomplished easily by modifying the last part of \texttt{reset} to
\begin{lstlisting}
if (control(3)>0)  // control(3) determines whether the sd's are active
{
  for (int j=1;j<=ngroups;j++)
  {
    sd(j)=boundp(x(ii),sdmin(j),sdmax(j),zpen); //Note that ii is not
                                                // incremented
  }
  ii++;   // Now increment ii
}
\end{lstlisting}

Another modification would be to enable individual mean or standard deviations
to be fixed or active.

\section{Getting the initial \texttt{x} vector}
\X{\fontindexentry{tt}{boundpin}}
\XX{bounding functions}{\fontindexentry{tt}{boundpin}}

It should be clear that the vector \texttt{x} is a transitory object associated
with the minimization routine. Without the function \texttt{reset}, it is not
even clear to what model parameter a particular component of \texttt{x}
corresponds or what value of the parameter it produces. For this reason, it does
not makes sense to save the values of \texttt{x} during or after the analysis.
It is the model parameters \texttt{p}, \texttt{mu}, and \texttt{sd} that should
be saved. However, this raises a question. Given the values of the model
parameters, the values of the control vector, and the values of the bounds on
the parameters, what should the initial values of the \texttt{x} vector be? This
problem is solved by the functions \texttt{boundpin} and \texttt{xinit}. \input
xinit.cpp

The control structure of \texttt{xinit} exactly parallels that of
\texttt{reset}, so no matter what switch combinations or bounds are used, the
correct values will be put into the components of~\texttt{x}. The function
\begin{lstlisting}
  double boundpin(double x,double fmin,double fmax)
\end{lstlisting}
is the inverse of the function \texttt{boundp}. Notice that the control vector,
\texttt{control}, determines how many active parameters there are.

The job of the function \texttt{nvarcal} is to use this information to calculate
the number of active parameters.
\input nvarcal.cpp

The three functions \texttt{reset}, \texttt{xinit}, and \texttt{nvarcal} all
share a common control structure and make it possible to coordinate the various
user options in the estimation routine.

Here is the code for the main routine in the mixture analysis program. The data
are assumed to be in a file named \texttt{mixture.dat}. The initial values for
the parameters are assumed to be in a file named \texttt{mixture.par}. \input
mixture.cpp

The job of the function \texttt{fcomp} is to put the \texttt{x} values into the
model by calling the function \texttt{reset}, and then to evaluate the
log-likelihood by calling the function \texttt{log\_likelihood}.
\input fcomp_m.cpp

Notice that while the parameterization of the mixture proportions,~\texttt{p},
has restricted the proportions to lie between~$0$ and~$1$, the proportions have
not been restricted, so they sum to~1. The function \texttt{normalize\_p}
restricts the~\texttt{p} so that they sum to~$1$.
\input normaliz.cpp
\X{penalty function}
The penalty $1000.*\log(\texttt{psum})*log(\texttt{psum})$ has been added to the
function being minimized, to remove the degeneracy in the parameterization of
the~\texttt{p}. Without this penalty, the minimizing value would be independent
of the value of \texttt{psum}. Including the penalty ensures that the minimum
value will occur for the value $\texttt{psum}=1$.

\section{Saving the parameter estimates}

\scAD\ uses the same memory locations for the derivative calculations as are
used for holding the variable objects. As a consequence, after the derivatives
have been calculated, the current value of the parameter estimates has been
lost. In order to save the results after the minimization has been carried out,
it is necessary to calculate the parameters one more time and print them into a
file. The simplest way to do this is to use the code for \texttt{fcomp} that
already exists, and to put a switch into \texttt{fcomp} so that saving the
parameters is enabled. The parameters are written into the file
\texttt{mixture.par} by the routine \texttt{save\_pars}.
\input savepar.cpp

%xx \htmlnewfile

\chapter{A Neural Network}\label{ch:neural-networks}
\X{neural net}
\XX{neural net}{feed-forward}

This example illustrates the use of \scAD\thinspace's 2- and 3-dimensional
ragged arrays.

\section{Description of the neural network}
\label{sec:description-nn}

An $n$-layer feed-forward neural network can be viewed as a process that, given
a set of $m_1$ inputs $(x_{11},\ldots,x_{1m_1})$ at layer $1$, produces a set of
$m_n$ outputs $(x_{n1},\ldots,x_{nm_n})$ at layer $n$. Level~$1$ is referred to
as the ``input layer.'' Level~$n$ is referred to as the ``output layer.''
Level~$2$ through Level~$n-1$ are referred to as ``hidden layers.'' The
quantities $x_{ij}$ are referred to as ``nodes.''
\XX{neural net}{input layer}
\XX{neural net}{output layer}
\XX{neural net}{nodes}

In a feed-forward neural network, the values of the nodes in each layer are
determined by the values in the previous layer via the relationship
\begin{equation}\label{eq:neural-network-01} %  \eqno(\number\mychapno.1)
  x_{i+1,j}=\Phi\left(\sum_{k=1}^{m_{i-1}}w_{ijk}x_{ik}+b_{i+1,j}\right)
\end{equation}
The function $\Phi$ is referred to as a ``squashing function.'' To the best of
our knowledge, the particular form of $\Phi$ does not seem to be especially
important, so long as it is a continuously differentiable monotone increasing
function that maps $(-\infty,\infty)$ into some desired bounded set $(a,b)$. The
bounded set is often of the form $(-0.5,0.5)$ or~$(0,1)$. The $w_{ijk}$ are
referred to as ``weights,'' while the $b_{ij}$ are called ``bias terms.''
\XX{neural net}{squashing function}
\XX{neural net}{bias terms}

The shape or topology of the neural net is completely determined by the number
of levels, $n$, and the number of nodes in each level, $m_i$. Assume that the
$m_i$ are contained in an \texttt{ivector} object \texttt{num\_nodes}.

\section{Implementation of the neural network}

\X{2-dimensional ragged array}
The $x_{ij}$ can be viewed as a 2-dimensional ragged array with $n$~rows and
$m_i$~elements in the $i^\textrm{th}$~row. Such an array can be created by the
declaration
\X{ragged matrix}
\begin{lstlisting}
  dmatrix x(1,num_levels,1,num_nodes); // these are the nodes
\end{lstlisting}
The bias terms $b_{ij}$ can also be viewed as a ragged matrix \texttt{B} with
$n-1$ rows. The declaration for~\texttt{B} requires a slight modification of the
\texttt{ivector} called \texttt{num\_nodes} to reflect the fact that \texttt{B}
does not need to contain row~1. The legal row index values for \texttt{B} should
be from~2 to~$n$.
\begin{lstlisting}
ivector iv(2,num_levels);
for (int i=2;i<=num_levels;i++)
{
  iv(i)=num_nodes(i); // Put the desired elements of num_nodes into iv
}
// Now declare B
dmatrix B(2,num_levels,1,iv);
\end{lstlisting}

For each $i$, $1\le i<n$, the weights, $w_{ijk}$, can be viewed as an $m_{i+1}
\times m_i$~matrix. Thus, the $w_{ijk}$ can be viewed as a ragged 3-dimensional
array \texttt{W}. To create \texttt{W}, we need some \texttt{ivector} objects to
contain the variable number of rows and columns.
\X{3-dimensional ragged array}
\begin{lstlisting}
ivector nrows(1,num_levels-1);
ivector ncols(1,num_levels-1);
for (i=1;i<num_levels;i++)
{
  nrows(i)=num_nodes(i+1);
  ncols(i)=num_nodes(i);
}
// Now declare W
d3_array W(1,num_levels-1,1,nrows,1,ncols);
\end{lstlisting}

If the squashing function $\Phi$ is extended to vector objects, by acting
elementwise on the elements of the vector, the code segment to carry out the
calculations in equation~(\ref{eq:neural-network-01}) %(\number\mychapno.1)
can be written as
\begin{lstlisting}
for (int i=1;i<num_levels; i++)
{
  x[i+1]=Phi(W[i]*x[i]+B[i+1]);
}
\end{lstlisting}

\section{ Initializing the vector of active parameters\br and resetting
active parameter values\br with the \texttt{set\_value} functions}

To set up the initial vector of parameter estimates needed by the function
minimizer, it is necessary to copy all the values of the active parameters. To
copy values into or from 2- and 3-dimensional objects, one must use a series of
nested loops. This process is tedious and prone to error. A collection of
functions has been included in \scAD\ to carry out these sort of tasks. The
function
\begin{lstlisting}
  set_value_inv(A,y,ii); // this is the ``inverse'' of the set_value function
\end{lstlisting}
will copy the elements of \texttt{A} into the \texttt{dvector} called
\texttt{y}, starting at the offset specified by the \texttt{int ii}. \texttt{ii}
will be incremented so that it points at the correct offset in \texttt{y} to add
the values of other objects. \texttt{A} may be a \texttt{dvector}, a
\texttt{dmatrix}, or a \texttt{d3\_array}. For a \texttt{dvector} called
\texttt{v}, the statements
\begin{lstlisting}
int ii=6;
set_value_inv(v,y,ii);
\end{lstlisting}
will have the same effect as
\begin{lstlisting}
int ii=6;
for (int i=v.indexmin();i<=v.indexmax();i++)
{
  y(ii++)=v(i);
}
\end{lstlisting}
\X{\fontindexentry{tt}{set\_value\_inv}}
\X{\fontindexentry{tt}{set\_value}}
\X{\fontindexentry{tt}{size\_count}}
\X{counting the number of elements in container classes}

The function
\begin{lstlisting}
  set_value(A,y,ii);
\end{lstlisting}
takes the values from the \texttt{dvar\_vector} \texttt{y}, beginning at the
offset \texttt{ii}, and puts them into the object \texttt{A}, which may be a
\texttt{dvar\_vector}, a \texttt{dvar\_dmatrix}, or a \texttt{dvar3\_array}. Of
course, it is the user's responsibility to declare the vector object~\texttt{y}
to have the appropriate size. To count up the number of elements in the
container objects, the function \texttt{size\_count()} has been defined. The
functions \texttt{set\_value()}, \texttt{set\_value\_inv()}, and
\texttt{size\_count()} have been used in the neural net code to simplify the
construction of the following functions:
\begin{itemize}
  \item \texttt{nvar\_calc}, which calculates the number of active parameters to
  be estimated,

  \item \texttt{yinit}, which puts the intial value of all the active parameters
  into the vector \texttt{y} for the function minimizer, and

  \item \texttt{reset}, which calculates the current estimate of the active
  parameters from the vector~\texttt{y}.
\end{itemize}
The code for these functions follows.
\begin{lstlisting}
// Count up the number of active parameters
int nvar_calc(dmatrix& B, d3_array& W)
{
  int nvar=0;
  nvar+=size_count(B);
  nvar+=size_count(W);
  return nvar;
}

// Put the parameter values from W and B into the vector y
void yinit(d3_array& W,dmatrix& B,independent_variables& y,int& num_levels,
  ivector& num_nodes)
{
  int ii=1;
  set_value_inv(W,y,ii);
  set_value_inv(B,y,ii);
}

// Put the values from the dvar_vector y into W and B
void reset(dvar3_array W,dvar_matrix& B,dvar_vector& y,int& num_levels,
  ivector& num_nodes)
{
  int ii=1;
  set_value(W,y,ii);
  set_value(B,y,ii);
}
\end{lstlisting}

\section{Training the network}
\XX{neural net}{training set} \XX{neural net}{weights}
\XX{neural net}{bias terms}

A training set for a neural net is a collection of input layer vectors together
with a corresponding set of output layer vectors. For each input layer vector,
the corresponding output layer vector can be thought of as the desired correct
output response of the network for the given input stimulus. The weights and
bias terms of the network are adjusted so that the difference between desired
response and the actual response is minimized. We shall assume that the training
set is contained in a file called \texttt{learning.smp}. The training set will
be read in every time the response of the neural net is calculated. This
provides for very large training sets that could not be kept in memory.
\begin{lstlisting}
ifstream infile("learning.smp");  // Open the file with the training data
infile >> num_learn;  // read in the number of examples to be used
// Loop through the training data
for (int k=1;k<= num_learn; k++)
{
  infile >> x[1];
  if (!infile)
  {
    cerr << "Error reading x[1] from file learning.smp\n";
    exit(1);
  }
  infile >> correct_response;
  if (!infile)
  {
    cerr << "Error reading correct_response from file learning.smp\n";
    exit(1);
  }
  for (int i=1;i<num_levels; i++)
  {
    x[i+1]=Phi(W[i]*x[i]+B[i+1]);
  }
  dvar_vector vdiff=correct_response-x[num_levels];
  z+=vdiff*vdiff; // z sums up the squared difference between the
                  // desired response and the actual response

  // Add the penalty terms
  dvariable v;
  v=norm(W);
  z+=.01*v*v;
  v=norm(B);
  z+=.01*v*v;
}
\end{lstlisting}

The parameters for \texttt{W} and \texttt{B} are chosen such that the total
squared error is minimized. A small quadratic penalty function of the model
parameters is added to keep the values of the parameters close to zero. This is
much like the use of ridge regression \cite{hoerlkennard1970} in linear models.
We don't know whether or not the use of such penalty functions really improves
the performance of the neural net, but it seemed to produce good results for the
problem considered here.
\X{ridge regression}

\section{The squashing function}

\XX{squashing function}{use of inverse tangent for}
For the squashing function, $\Phi$, we have employed the inverse tangent
function~\texttt{atan}.
\begin{lstlisting}
dvar_vector Phi(dvar_vector& v)
{
  // This squashing function maps into the interval (-.5,.5)
  RETURN_ARRAYS_INCREMENT;  // Need this because the function returns a
                            // variable type
  int min=v.minindex();
  int max=v.maxindex();
  dvar_vector tmp(min,max)
  for (int i=min;i<=max;i++)
  {
    tmp(i)=0.31831*atan(v(i)); //0.31831 is 1/PI
    // tmp(i)=0.31831*atan(v(i))+0.5; // This maps into (0,1)
  }
  RETURN_ARRAYS_DECREMENT;
  return(tmp);
}
\end{lstlisting}

\section{Specifying the structure of the neural net\br
  and getting initial parameter estimates}

The structure of the neural net is determined by specifying the number of layers
and the number of nodes in each layer. This information is put into the file
\texttt{struct.nrl}.

To begin training the neural net, the \texttt{d3\_array} called \texttt{W} and
the \texttt{dmatrix} called \texttt{B} were filled with normally distributed
random numbers, with mean 0 and standard deviation~1.0. The value of~1.0 for the
standard deviation was picked by trial and error. If the standard deviation is
too small, the minimization routine doesn't seem to be able to find a good
minimum and the procedure ``stalls'' at a bad local minimum.

The program \texttt{WINIT.EXE} reads the file \texttt{struct.nrl} and writes the
initial parameter values into the file \texttt{WEIGHTS.PAR}.

\section{A three spiral problem}

Touretsky and Pomerlau \cite{touretzkypomerlau1989} discussed the problem of
training a neural network to recognize points that lie on or near a double
spiral. We thought it would be interesting to extend this problem to points that
lie on or near a triple spiral.

\begin{wrapfigure}{r}{.33\textwidth}%
      \centering
        %{\includegraphics[height=2in]{foo.eps}}%
        {\input grph5.tex \par}%
      \emptycaption{}%
      \label{fig:description-nn-01}
\end{wrapfigure}

The training set for the neural net is shown in
Figure~\ref{fig:description-nn-01}. % \number\mychapno.1.
The net is being trained to determine whether a point in a disk lies on or near
one of three spirals. The input level to the net consists of two nodes: the
Cartesian coordinates of the point. The output level of the net has three nodes.
The node corresponding to the spiral on which a particular point is located is
intended to be either turned on (have value $+.5$), or turned off (have value
$-.5$). Actually, any value for a node that is less than zero is considered
``off'' while any value which is greater than $0$ is considered ``on.'' Due to
the nature of the inverse tangent squashing function, the maximum value produced
by the neural net is about $0.42$, while the minimum value is about $-0.42$. If
the neural net has $0$, $2$, or $3$ output nodes turned on ($>0$) for a point,
then that point is considered unclassified.

A typical input consists of the pair of coordinates of a point on the spiral
$(0.0618,0.1901)$ together with the triple $(0.5,-0.5,-0.5)$, which is the
desired output for this point, since it lies on the first spiral arm.

We trained neural nets of increasing complexity until a design was obtained that
had enough parameters to classify all the points in the learning sets correctly.
The neural net we chose has 6 layers with 2 input nodes, 8 nodes in each of the
4 hidden layers, and 3 output nodes. This neural net has 267 parameters to be
estimated.

\X{\fontindexentry{tt}{GRADFIL1.TMP}}
Training the neural net required about 20 hours on a 20 MHz 386 computer. The
file \texttt{GRADFIL1.TMP}, which stores the temporary information needed to
calculate the gradient, was about 3 MB in size. It was stored on a \textsc{ram}
disk. At any time, it is possible to stop the program by typing `q.' The current
parameter estimates are printed into the file \texttt{EST.PAR}. To restart the
learning process, you must first copy the file \texttt{EST.PAR} into the file
\texttt{WEIGHTS.PAR}. Then it is possible to carry out training sessions during
periods when the computer is not required.
\begin{figure}[h]
  \centering\hskip1pt
   \input grph6.tex
   \emptycaption{}
   \label{fig:description-nn-02} % \number\mychapno.2
\end{figure}
Figure~\ref{fig:description-nn-02} % \number\mychapno.2
illustrates how the neural net has classified the points in the disk. The points
in \texttt{A}, \texttt{B}, and \texttt{C} have been classified as being
associated with the corresponding spiral. The small points in \texttt{D} are
unclassified. The large points in \texttt{D} are the training set. The
unclassified points appear to mark fairly regular boundaries between the
spirals. Our solution does not seem to exhibit the ``bumps and gaps'' reported
by Touretsky and Pomerlau for the double spiral solution. Perhaps this is due to
the superiority of our minimization procedure over the method they employed. It
may also be due to the use of more parameters, combined with a small quadratic
penalty function on the parameters to introduce stability into the model. In any
event, we seem to have obtained a nice solution to the problem.

\bibliographystyle{plain}
\bibliography{autodif}

\printindex  % Make the index here.

\end{document}
