<reponame>wStockhausen/admb<gh_stars>10-100
% $Id$
%
% Author: <NAME>
% Copyright (c) 2008 Regents of the University of California
%

\section{The necessity for adjoint code}

When you write code for variable objects in \ADM, all the
derivatives are calculated for you in a transparent manner. 
To accomplish this, \ADM\ must save certain information for
later use. We shall refer to this as derivative information.
Each arithmetic operation generates about 32~bytes of
derivative information.  If you have some simple function that
has 20 arithmetic operations, it will therefore
generate 640 bytes of derivative information every time
it is called. The purpose of writing adjoint code is to
reduce the amount of derivative information that must
be calculated. For a function that is called many times,
this can greatly reduce the amount of derivative information that
must be stored.

In this chapter, we investigate how to write and debug
adjoint code. To begin, we investigate how to write adjoint code for a
simple function that takes 1--4 independent variables and returns 1~dependent variable. 
The adjoint code for such functions is
simpler to write than that for a general function---such as the
singular value decomposition of a matrix, which we will consider later.


\section{Writing adjoint code: a simple case}

Consider a simple function $f$ which takes 1~independent variable $x$
and returns a dependent variable $y$, i.e.,
$$y=f(x)$$
where $f(x)=\exp(-x^2/2)$. The code for this example can be written like
\begin{lstlisting}
  dvariable errf(const prevariable& x)
  {
    return exp(-0.5*square(x));
  }
\end{lstlisting}
There are three arithmetic operations here: square, multiplication, and exponentiation,
so 96 bytes of derivative information will be generated. (Actually, the \texttt{return} operation
also generates 32 bytes of derivative information, but we will ignore that for now.)
A less efficient way to write the code (but more useful for 
showing adjoint code for this simple example) would be:
\begin{lstlisting}
  dvariable errf(const prevariable& x)
  {
    dvariable y;
    y=exp(-0.5*square(x));
    return y;
  }
\end{lstlisting}
Here is the same code with the derivative calculated by the one line of adjoint
code \texttt{double dfx=-value(x)*value(y);}:
\begin{lstlisting}
  dvariable errf(const prevariable& x)
  {
    dvariable y;
    value(y)=exp(-0.5*square(value(x)));
    double dfx=-value(x)*value(y);
    AD_SET_DERIVATIVES1(y,x,dfx);  // 1 dependent variable 
    return y;
  }
\end{lstlisting}
So what is going on here? Consider the line
\begin{lstlisting}
  value(y)=exp(-0.5*square(value(x))); 
\end{lstlisting}
The \texttt{value} function returns a constant type, that is, a double, 
that has the same value as the corresponding \texttt{dvariable} or prevariable.
In fact, it is the {\it same} object.  That is, it shares the same address,
but the type has been changed to double. So, the above line of
code assigns the value \texttt{exp(-0.5*square(value(x)))} to \texttt{y},
but without generating any derivative code. Similarly, since
the calculations are made on \texttt{value(x)}, these calculations
will not generate any derivative code. So, it is the responsibility of the
programmer to calculate the derivative code and store it where it
can be used later. The line
\begin{lstlisting}
  double dfx=-value(x)*value(y);
\end{lstlisting}
calculates the derivative $f^\prime(x)$
of $y$ with respect to $x$ and stores it with the line of code
\begin{lstlisting}
  AD_SET_DERIVATIVES1(y,x,dfx);
\end{lstlisting}
This code will only generate 32 bytes of derivative information.


\section{Debugging adjoint code: a simple case}

The simplest way to debug the adjoint code is to put your new
function into an \ADM\ template file and use the \texttt{-dd 1}
command line option to call the derivative checker.
\XX{\fontindexentry{tt}{-dd}}{command line argument}
\XX{command line arguments}{\fontindexentry{tt}{-dd}}
\begin{lstlisting}
DATA_SECTION
PARAMETER_SECTION
  init_number x
 !! x=2;
  objective_function_value f;
PROCEDURE_SECTION
  f=square(errf(x));
  
GLOBALS_SECTION
  #include <admodel.h>
  dvariable errf(const prevariable& x)
  {
    dvariable y;
    value(y)=exp(-0.5*square(value(x)));
    double dfx=-value(x)*value(y);
    AD_SET_DERIVATIVES1(y,x,dfx);  // 1 dependent variable 
    return y;
  }
\end{lstlisting}


\section{Adjoint code for more than\br one independent variable}

The following code shows how to write the adjoint code for a function with two
independent variables:
\begin{lstlisting}
DATA_SECTION
  vector lengths(1,10)
  vector ages(1,10)
  !! lengths.fill_seqadd(1,1);
  !! ages.fill_seqadd(1,1);
  !! lengths=sqrt(lengths);

PARAMETER_SECTION
  init_bounded_number linf(0,10) 
  init_bounded_number rho(0,1) 
  objective_function_value f;
PROCEDURE_SECTION
  for (int i=1;i<=10;i++)
    f+=square(lengths(i)-vb_growth(linf,rho,ages(i)));
  
GLOBALS_SECTION
  #include <admodel.h>

  dvariable vb_growth(const prevariable& linf, const prevariable& rho, 
    double t)
  {
    double clinf=value(linf);
    double crho=value(rho);
    dvariable len;
    value(len)=clinf*(1-pow(crho,t));
    double dflinf=1-pow(crho,t);
    double dfrho=-clinf*t*pow(crho,t-1);
    AD_SET_DERIVATIVES2(len,rho,dfrho,linf,dflinf);  // 2 dependent variable 
    return len;
  }
\end{lstlisting}
This approach to writing adjoint code has been implemented for functions of
up to four independent variables.
\begin{lstlisting}
DATA_SECTION
  vector lengths(1,10)
  vector ages(1,10)
  !! lengths.fill_seqadd(1,1);
  !! ages.fill_seqadd(1,1);
  !! lengths=sqrt(lengths);

PARAMETER_SECTION
  init_bounded_number linf(0,10) 
  init_bounded_number rho(0,1) 
  init_number t0 
  init_bounded_number gamma(.1,1.9) 
  objective_function_value f;
PROCEDURE_SECTION
  for (int i=1;i<=10;i++)
    f+=square(lengths(i)-vb_growth(linf,rho,t0,gamma,ages(i)));
  
GLOBALS_SECTION
  #include <admodel.h>

  dvariable vb_growth(const prevariable& linf, const prevariable& rho, 
    const prevariable& t0, const prevariable gamma,double t)
  {
    double clinf=value(linf);
    double ct0=value(t0);
    double crho=value(rho);
    double cgamma=value(gamma);
    dvariable len;
    value(len)=pow(clinf*(1-pow(crho,t-ct0)),cgamma);
    double tmp=cgamma*pow(clinf*(1-pow(crho,t-ct0)),cgamma-1);
    double dflinf=tmp*(1-pow(crho,t-ct0));
    double dft0=tmp*(clinf*log(crho)*pow(crho,t-ct0));
    double dfrho=-tmp*clinf*(t-ct0)*pow(crho,t-ct0-1);
    double dfgamma=value(len)*log(clinf*(1-pow(crho,t-ct0)));
    AD_SET_DERIVATIVES4(len,t0,dft0,rho,dfrho,linf,dflinf,gamma,dfgamma); 
       // 4 dependent variable 
    return len;
  }
\end{lstlisting}


\section{Structured calculation of derivatives in adjoint code}

Until now, we have deliberately calculated the derivatives with respect to the independent
variables in an {\it ad-hoc} fashion. While this approach works for simple functions,
it rapidly becomes untenable when the function is more complicated. In the following example,
we have calculated the derivatives in a more structured fashion. Notice that to calculate the
derivatives, every line of code in the function is repeated in the opposite
order (commented out, of course) and the corresponding derivatives are calculated.
\begin{lstlisting}
DATA_SECTION
  vector lengths(1,10)
  vector ages(1,10)
  !! lengths.fill_seqadd(1,1);
  !! ages.fill_seqadd(1,1);
  !! lengths=sqrt(lengths);

PARAMETER_SECTION
  init_bounded_number linf(0,10) 
  init_bounded_number rho(0,1) 
  init_number t0 
  init_bounded_number gamma(.1,1.9) 
  objective_function_value f;
PROCEDURE_SECTION
  for (int i=1;i<=10;i++)
    f+=square(lengths(i)-vb_growth(linf,rho,t0,gamma,ages(i)));
  
GLOBALS_SECTION
  #include <admodel.h>

  dvariable vb_growth(const prevariable& linf, const prevariable& rho, 
    const prevariable& t0, const prevariable gamma,double t)
  {
    double clinf=value(linf);
    double ct0=value(t0);
    double crho=value(rho);
    double cgamma=value(gamma);
    dvariable len;
    double u1=pow(crho,t-ct0);
    double u2=clinf*(1-u1); 
    value(len)=pow(u2,cgamma);
    double dflen=1.0;
    //value(len)=pow(u2,cgamma);
    double dfu2=dflen*cgamma*pow(u2,cgamma-1.0);
    double dfgamma=dflen*value(len)*log(u2);
    //double u2=clinf*(1-u1); 
    double dflinf=dfu2*(1-u1);
    double dfu1=-dfu2*clinf;
    //double u1=pow(crho,t-ct0);
    double dfrho=dfu1*(t-ct0)*pow(crho,t-ct0-1.0);
    double dft0=-dfu1*u1*log(crho);

    AD_SET_DERIVATIVES4(len,t0,dft0,rho,dfrho,linf,dflinf,gamma,dfgamma);  
      // 4 dependent variable 
    return len;
  }
\end{lstlisting}


\section{General adjoint code}

So far, the adjoint code has been for a simple function which has from~1 to~4
independent variables and returns 1~dependent variable. 
Now we consider the general case where the function can take
any number of dependent variables and return any number of 
dependent variables, and these variables can be in the form of numbers, 
vectors, or matrices.
\begin{lstlisting}
 dvar_vector  operator * (_CONST dvar_matrix& m,_CONST dvar_vector& x )
 {
   if (x.indexmin() != m.colmin() || x.indexmax() != m.colmax())
   {
     cerr << " Incompatible array bounds in dvar_vector  operator * "
            << "(_CONST dvar_matrix& m,_CONST dvar_vector& x)\n";
     ad_exit(21);
   }

   dvar_vector tmp(m.rowmin(),m.rowmax());
   double sum;

   for (int i=m.rowmin(); i<=m.rowmax(); i++)
   {
     sum=0.0;
     for (int j=x.indexmin(); j<=x.indexmax(); j++)
     {
       sum+=(m.elem(i)).elem_value(j)*x.elem_value(j);
     }
     tmp.elem_value(i)=sum;
   }
  save_identifier_string("PLACE4");
  x.save_dvar_vector_value();
  x.save_dvar_vector_position();
  save_identifier_string("PLACE3");
  m.save_dvar_matrix_value();
  m.save_dvar_matrix_position();
  save_identifier_string("PLACE2");
  tmp.save_dvar_vector_position();
  save_identifier_string("PLACE1");
  
  ADJOINT_CODE(dmdv_prod);
   return(tmp);
 }
\end{lstlisting}
To calculate the adjoint code, it will be necessary to have the
values of the matrix \texttt{m} and the vector \texttt{x}.
This is accomplished with the instructions
\begin{lstlisting}
  x.save_dvar_vector_value();
  m.save_dvar_matrix_value();
\end{lstlisting}
Also, to calculate the derivatives, it will be necessary to know where the
derivatives with respect to the independent and dependent variables are 
located. This information is saved with the instructions
\begin{lstlisting}
  x.save_dvar_vector_position();
  m.save_dvar_matrix_position();
  tmp.save_dvar_vector_position();
\end{lstlisting}
Finally, we need to save the name of the routine that calculates
the adjoint code, so that it can be called at the appropriate time. To write the
code for the adjoint calculations, keep in mind that everything must be recovered
from the stack in the reverse order from which it was put on the stack.
This process can be a bit confusing and if you don't do it properly, the 
stack will become corrupted and nothing will work.  To help diagnose
problems, function \texttt{save\_identifier\_string} can be
used to put a string on the stack. This string value can be checked in the
adjoint code with the \texttt{verify\_identifier\_string} function. At least two of
these functions should be left in any adjoint code, so that stack
integrity can be monitored if problems show up later. For optimized code,
they are not used, and so contribute almost nothing to the overhead.

The adjoint code begins by reading the information that was saved on the stack.
An object of type \texttt{dvar\_vector\_position} contains both the size and
address information associated with a \texttt{dvar\_vector}---which are needed to
recover or store derivative values, or to build a \texttt{dvector} with the
same shape as the {\t dvar\_vector}.
The function \texttt{restore\_dvar\_vector\_der\-ivatives} gets the values of the
derivatives with respect to the dependent variables, so they can be used in the
adjoint code. The functions
\begin{lstlisting}
   dfx.save_dvector_derivatives(x_pos);
   dfm.save_dmatrix_derivatives(m_pos);
\end{lstlisting}
use the position information to save the derivatives with respect to the
independent variables in the appropriate places.
\begin{lstlisting}
void dmdv_prod(void)
{
  verify_identifier_string("PLACE1");
  dvar_vector_position tmp_pos=restore_dvar_vector_position();
  verify_identifier_string("PLACE2");
  dvar_matrix_position m_pos=restore_dvar_matrix_position();
  dmatrix m=restore_dvar_matrix_value(m_pos);
  verify_identifier_string("PLACE3");
  dvar_vector_position x_pos=restore_dvar_vector_position();
  dvector x=restore_dvar_vector_value(x_pos);
  verify_identifier_string("PLACE4");
  dvector dftmp=restore_dvar_vector_derivatives(tmp_pos);
  
  dmatrix dfm(m_pos);
  dvector dfx(x_pos.indexmin(),x_pos.indexmax());
  dfm.initialize();
  dfx.initialize();
  
   double dfsum;
   for (int i=m.rowmax(); i>=m.rowmin(); i--)
   {
     // tmp.elem_value(i)=sum;
     dfsum=dftmp.elem(i);
     for (int j=x.indexmax(); j>=x.indexmin(); j--)
     {
       //sum+=(m.elem(i)).elem_value(j)*x.elem_value(j);
       dfm.elem(i,j)+=dfsum*x.elem(j);
       dfx.elem(j)+=dfsum*m.elem(i,j);
     }
     //sum=0.0;
     dfsum=0.0;
   }
   dfx.save_dvector_derivatives(x_pos);
   dfm.save_dmatrix_derivatives(m_pos);
}
\end{lstlisting}
