text:
\documentclass{article} % For LaTeX2e
\usepackage{multicol}
\usepackage{nips13submit_e,times}
\usepackage{hyperref}
\usepackage{url}
\usepackage{amsmath}
\usepackage{amsfonts}
\usepackage{breakurl}
\usepackage{graphicx}
\usepackage{caption}
\usepackage{subcaption}
\usepackage{float}
\usepackage{fancyvrb}
\usepackage{booktabs}
\usepackage[font=footnotesize,labelfont=bf]{caption}
\usepackage[margin=1in]{geometry}
% \documentstyle[nips13submit_09,times,art10]{article} % For LaTeX 2.09
\bibliographystyle{ieeetr}

\title{Predicting water pump failure in Tanzania and optimising maintenance routes}

\author{
Sam Kim \\
Harvard College\\
Cambridge, MA 02138 \\
\texttt{samuelkim@college.harvard.edu} \\
\And
Gareth Haslam \\
Harvard Extension School\\
Cambridge, MA 02138 \\
\texttt{haslam.gareth@gmail.com} \\
}

% The \author macro works with any number of authors. There are two commands
% used to separate the names and addresses of multiple authors: \And and \AND.
%
% Using \And between authors leaves it to \LaTeX{} to determine where to break
% the lines. Using \AND forces a linebreak at that point. So, if \LaTeX{}
% puts 3 of 4 authors names on the first line, and the last on the second
% line, try using \AND instead of \And before the third author name.

\newcommand{\fix}{\marginpar{FIX}}
\newcommand{\new}{\marginpar{NEW}}

\nipsfinalcopy % Uncomment for camera-ready version

\begin{document}

\maketitle

\begin{abstract}
We describe MCMC techniques to model the functionality of water pumps in Tanzania and suggest an optimised route for their maintenance. Using field data describing various attributes of each pump, such as construction year, installer, and location, we describe a method to predict the current status of the pump from a set of three possibilities (functioning, functioning needs repair, and not functioning). We also optimize the route that a maintenance crew could take to repair damaged pumps in a version of the well--known, NP--hard, traveling salesman problem. We find acceptable solutions using stochastic metaheuristics. Code and more information are available at \url{http://ghaslam.github.io/AM207/}.
\end{abstract}

\begin{multicols}{2}

\section{Introduction}

Access to clean water is a fundamental need for people all over the world. Around 11\% of the world's population currently lack such access due to either lack of natural resources or lack of infrastructure \cite{UN2013}. In Tanzania, despite the availability of water, only 33.5\% have access to a piped source \cite{Morisset2012}. Many others depend on pumps for their supply. When these pumps fail, that can create serious problems for the people who depend on them. Previous work on predicting the failure of mechanical devices such as pumps has benefitted from the use of both statistical analysis as well as data obtained from direct observation of the machines, for example diagnostic vibration sensors \cite{Nakamura2007}. As the popularity of machine learning and 'big data' continues to grow, researchers are also studying how failure can be predicted based on training data from large historical datasets, for example for rod pumps in the oil industry \cite{Liu2013}. Liu's research used subject matter experts to first identify whether a pump was operating normally or about to fail and then trained their model to recognise the features in the data such as daily usage, that predicted the label.

We aim to use Bayesian Methods to assign one of three categorical labels indicating the functional status each pump (functioning, functioning needs repair, and not functioning). The Bayesion models will be built on training data, then used to predict the functional status for wells with unknown status. Then, we treat the repair of the damaged pumps as a constrained optimisation problem over the space of all possible routes that the maintenance crew could take to reach the damaged wells. We ignore the existing road infrastructure, terrain, height etc. and assume that all locations are equally accessible.

\section{The Data}

The data for this project was provided as part of a data science challenge run by DrivenData.org \cite{DrivenData2015}. DrivenData aims to encourage data scientists from around the world to take part in online competitions to develop predictive statistical models that can address important problems in the fields of health, international development, and the environment. The data for the Tanzanian water pump competition was aggregated from the Tanzanian Ministry of Water by Taarifa, an open source platform for tracking infrastructure related issues \cite{Taarifa2015}.  The dataset contains records of 59,400 water pumps with information on a range of 39 parameters, such as pump construction year, installer, type of pump, water quality, population around the well, cost of water, quantity and others. These include a mix of categorical and numerical data, and there are also many missing values. Each pump is assigned a unique ID as well as precise coordinates for its location. A subset of the parameters are shown in Table 1. The final column in Table 1 indicates the current status of the pump which the model will attempt to predict. The optimisation algorithm will then attempt to find the best route that reaches the damaged pumps in the second part.

\begin{table}[H]
\caption{Sample of the DrivenData Tanzania water wells dataset}
\centering
\resizebox{\columnwidth}{!}{%
\begin{tabular}{lrlrlrrrr}
\toprule
{} &  Installer & GPS Height &  Longitude &                     Latitude &  Construction\_Year &  Source\_Type &  Quantity\_Type &  Status\_Group \\
\midrule
69572 &   Roman         & 1390             &       34.93 &  -9.86 &   1983 &                         Spring &          Enough  &           Functioning \\
8776 &   Grumeti          & 1399             &      34.69 &  -2.15 &   2002 &                          Rainwater &           Insufficient &          Functioning Needs Repair  \\
34310 &   World Vision & 686               &    37.46 &  -3.82 &   1967 &                           Dam &           Enough &           Functioning Needs Repair \\
67743 &   UNICEF       & 263                &   38.49 &  -11.15  &   1997 &                         Machine &           Dry &         Non-functioning \\
19728 &   Artisan          &  0                  &   31.13 &     -1.83 &   2007 &                         Rainwater &           Seasonal &           Functioning  \\
\bottomrule
\end{tabular}
}
\end{table}

In addition to the parameters given in the original dataset, we calculate some additional parameters such as age (based on the construction year), and relative remoteness (based on the average distance of the 5 nearest neighbours). Figure \ref{fig:map-phys} shows a map of the physical geography of Tanzania and another showing the locations of the pumps. The high density of pumps appears as a solid blue colour but actually represents many individual points.

\begin{figure*}
  \centering
  \begin{subfigure}[b]{0.5\textwidth}
    \centering
    \includegraphics[width=\textwidth]{figures/Tanzania}
    \caption{Physical geography of Tanzania.}
    \label{fig:map-phys}
  \end{subfigure}~\begin{subfigure}[b]{0.5\textwidth}
    \centering
    \includegraphics[width=\textwidth]{figures/PumpMap}
    \caption{Distribution of water pumps.}
    \label{fig:map-points}
  \end{subfigure}
  \caption{Geography of Tanzania and distribution of water pumps.}
  \label{fig:map}
\end{figure*}

\section{Modeling water well functionality}

\subsection{Model Classification}

Bayesian methods are only compatible with numerical data and do not work directly with classification problems, so we need to convert the categorical data and classes into numerical data and labels.

For the classification, we can take two approaches:
\begin{enumerate}
\item Assume that the 3 labels, "non functional", "functional needs repair", and "functional" lie on a spectrum where the distance between "non functional" and "functional" is greater than the difference between "non functional" and "functional needs repair" or "functional needs repair" and "functional. We can then have the Bayesian model predict a number $y$, and then convert this to a label using a range, where the exact limits can be adjusted:
\begin{align*}
\text{functional: } & 0.67\leq y_i \\
\text{functional needs repair: } & 0.33<y_i<0.67 \\
\text{non functional: } &  y_i \leq 0.33
\end{align*}
\item The assumption that the 3 labels lie on a linear spectrum is not necessarily a safe assumption, and the limits are set rather arbitrarily. A more natural and common method in machine learning is to build a classifier that decides between 2 classes. In this case, we can either build 3 one-versus-one classifiers or 3 one-versus-rest classifiers and use majority vote (or probabilities in the case of a tie) to make the final classification.
\end{enumerate}

We use the first approach for its simplicity, although further work would investigate the performance of the second approach, which is much more versatile at the cost of computational complexity.

During the training stage of our model, we also need to convert the given labels into numbers. Assuming that we use the first approach above for prediction, we use a similar scheme:
\begin{equation*}
y_i =
\begin{cases}
0 & \text{ non functional}\\
0.5 & \text{ functional needs repair}\\
1 & \text{ functional}
\end{cases}
\end{equation*}

\subsection{Categorical data}

Many of the data features are categorical data that do not have any numerical interpretation, so we cannot convert them the same way that we converted our label. For example, the "installer" feature has labels such as "UNICEF," "Roman," and "Artisan." To deal with these, we use one-of-k representation in which we add a new feature column for every unique value. For example, we would add the columns "installer=UNICEF," "installer=Roman," and "installer=Artisan." The "installer=UNICEF" would be 1 if that data's "installer" was "UNICEF," and 0 otherwise. So if there are $k$ unique values for a particular categorical feature, then for each data point, we add on a vector of length $k$ which has a single 1 and $k-1$ 0s.

\subsection{The model}

We model $y_i$ using a logistic function that constrains the range to $0 < y_i < 1$, and the parameter for the logistic function is controlled by the features and weights, which are also parameters that we need to calculate. Our model comes from the assumption that there are certain factors impacting decay rate, and so the functionality is dependent on this decay rate.
We then select features from the data which we have a prior belief to have a strong influence on functionality, choose priors on their distributions and construct a Metropolis-Hastings Sampler to select the optimum value for each parameter. New predictions of functionality can then be made from this model.

We also add on a term for noise, $\epsilon_i$, which we assume is Gaussian noise controlled by the standard deviation $\sigma$.

\begin{gather*}
y_i = \frac{1}{1+e^{\alpha_i}}+\sigma\epsilon_i \\
\alpha_i = \beta_0 + \beta_1 x_{i,1} + \beta_2 x_{i,2} + ... + \beta_n x_{i,n}
\end{gather*}

There are $n$ features and $n+2$ parameters ($n+1$ $\beta$s and $\sigma$). The parameters are found by sampling from the posterior:

\begin{equation*}
p(Y,\Theta)=p(Y|\Theta)p(\Theta)
\end{equation*}

\section{Prediction Results}

\subsection{Results}

\begin{figure*}
	\begin{subfigure}{0.5\textwidth}
		\includegraphics[width=\textwidth]{figures/trace_block}
		\caption{Blockwise updating}
		\label{fig:trace_block}
	\end{subfigure}
	\begin{subfigure}{0.5\textwidth}
		\includegraphics[width=\textwidth]{figures/trace_component}
		\caption{Componentwise updating}
		\label{fig:trace_comp}
	\end{subfigure}
	\caption{Traceplot of model parameters, using MH sampler.}
	\label{fig:trace}
\end{figure*}

Figure~\ref{fig:trace} shows the trace plot of the model parameters with blockwise updating and componentwise updating. These are after the burn-in period. Blockwise updating has an acceptance rate of nearly 0 because we are updating 8 parameters at a time. Componentwise updating has a much higher acceptance rate, although at a huge time cost.

The parameters chosen are "longitude," "latitude," "age" (2015-"construction\_year"), "gps\_height," "quantity=dry," and "population."

We see that $\beta_5$ has significant predictive power since this is the farthest from 0. This corresponds to the "quantity=dry" feature. $\beta_0$ is the constant offset, which we would expect to be non-zero.

By fine-tuning the step parameters and carefully choosing features, we reach a cross-validation prediction accuracy of 60.6\%.

\subsection{Comparison to machine learning approaches}

Here we compare our methods with other common machine learning methods.

The k-Nearest Neighbors algorithm (also known as the k-NN algorithm) is a machine learning algorithm that can be applied to classification problems. There is no training aspect to the algorithm. Predictions are made by taking the $k$ nearest data points to the desired data point in feature space, and taking the majority of their labels. The labels can also be weighted by $1/d$, where $d$ is the distance from the desired data point to the labeled data point. We apply the k-NN algorithm using only the longitude and latitude data, so we are picking the geographically nearest neighbors. Accuracy in cross-validation is shown in the below figure as a function of $k$. We see that it reaches 69\% prediction accuracy.

\includegraphics[width=8cm]{figures/kNN_accuracy_vs_K.png}

Decision trees are a machine learning method in which a tree of rules is built. One traverses down the tree, and each node represents a rule based on the data's features. The leaves are the final decision. The nodes are built based on which features have the most predictive power for classification. Random forests are an extension of decision trees in which many trees are trained, and stochasticity is introduced to make the trees different, which decreases bias. We show in randomforest\_explore.ipynb that random forests on a small subset of features can reach an accuracy of 79\%. Further experimentation and more sophisticated feature extraction can reach prediction accuracies of 83\% on this data.

Bayesian modeling is probably not the best method to use for this type of problem. It assumes that the classification can be described by some sort of model, and we would have to guess this model ahead of time. Decision trees and random forests do not assume any sort of relation between the variable features and the labels, and merely seek out the variable features with the highest predictive power. Second, there are a lot of parameters in the MH sampler that need to be tuned, including initial guesses for the parameters and step sizes. Additionally, Bayesian modeling takes extremely long to train, on the order of hours on a laptop, assuming that the initial guesses for the parameters and step sizes are set correctly. Training k-NN and random forests take on the order of minutes.

\section{Optimizing maintenance crew route}

In the second part of this project, we aim to suggest an optimised route by which a maintenance crew could visit the damaged pumps. The choice of locations can either be taken from the originally labeled pumps or based on the outputs of our model's labels.

\subsection{The traveling salesman problem}

Having hopefully identified a method for accurately predicting which pumps will fail, we can imagine that the Tanzanian Ministry of Water will want to task a maintenance crew with visiting each of these locations. We may also assume that the maintenance crew are based in the commercial capital of Dar Es Salaam. The crew wishes to visit all $N$ locations in as efficient manner as possible. This problem is equivalent to the traveling salesman problem (TSP), first described in 1932 by Karl Menger \cite{Menger}. A compact mathematical formulation of the problem is given by Miller et. al. \cite{Miller1960} as follows:

\begin{align*}
\textrm{Find variables }& X_{ij}  \textrm{ and } U_i i,j = 1, 2, \cdots , n \textrm{ that minimize}:\\
Z &= \min \sum_{i=0}^n \sum_{j\ne i,j=0}^n C_{ij} X_{ij}\\
\end{align*}

\begin{align*}
\textrm{subject to} & \\
	& \sum_{i=1}^n X_{ij} = 1 && j = 1, 2 , \cdots, n \\
	& \sum_{j=1}^n X_{ij} = 1 && i = 1, 2, \cdots, n \\
	&U_i - U_j + nX_{ij} \le n-1 && i, j = 2, \cdots, n i\ne j \\
	& x_{ij} \in \{0, 1\} && i,j=0, \cdots, n \\
\end{align*}

The optimal solution is subject to the following constraints:

\begin{itemize}
  \item $X_{ij}$ is a binary decision variable indicating whether to accept a move from location $i$ to location $j$.
  \item $C_{ij}$ is the distance between location $i$ and location $j$.
  \item The objective function is the sum of the distances for routes that we decide to take.
  \item Each location is visited once and only once.
\end{itemize}

The difficulty with solving the TSP is that for large numbers of locations, it becomes impossible to test all possible routes as the problem has $\mathcal{O}$$(n!)$ complexity and is classified as NP-hard. Therefore, it can be sufficient to find an answer that is 'good enough' even if it is not globally optimal. To find  this 'good enough' solution we use Simulated Annealing (SA).

Figure \ref{fig:TSP} shows the effect of simulated annealing on reducing the total distance travelled by the maintenance crew and suggest a near-optimal route between the pumps in need of repair. Due to the large number of locations, only a subset are shown in order to make viewing easier. The depot in Dar Es Salaam is highlighted in red. The trace plot shown in \ref{fig:Plain-TSP-Trace} shows that the SA algorithm has approached a stable route distance which is significantly shorter than the original random route.

\begin{figure*}
	\begin{subfigure}{0.5\textwidth}
		\includegraphics[width=\textwidth]{figures/Plain-TSP-Trace}
		\caption{Decrease in distance due to annealing}
		\label{fig:Plain-TSP-Trace}
	\end{subfigure}
	\begin{subfigure}{0.5\textwidth}
		\includegraphics[width=\textwidth]{figures/Optimised-Route}
		\caption{Optimised maintenance route with Dar Es Saleem in red.}
		\label{fig:trace_comp}
	\end{subfigure}
	\caption{Simulated annealing for travelling salesman.}
	\label{fig:TSP}
\end{figure*}

\subsection{Modifications}

With one maintenance crew, the total distance travelled is 41,000 km, which is obviously too much for one crew to travel. It is also important to note that we use the Haversine formula for calculating the distance between locations as the straight-line, Euclidean distance is imprecise for large distances on a sphere \cite{Sinnot1983}. An obvious extension of the problem is to imagine that instead of making a single visiting every pump, the crew make multiple trips, returning to the depot in Dar Es Salaam each time (or equivalently, that multiple crews leave and return to the depot each day). This leads to a decrease in total distance travelled from 41,000 km to 34,500 km. This problem is commonly referred to as the vehicle routing problem (VRP). Gavish et. al. describe a range of other interesting formulations of the TSP which would be relevant to our situation \cite{Gavish1978}. Finally, we looked at prioritising the damaged pumps that were identified as being 'remote', with the assumption that people depending on these pumps would have fewer alternative sources nearby that they could visit if their own pump failed. We first identified the average distance from each pump to its 5 nearest neighbours using the k-dimensional tree method. We then added this weight to the cost function used by the SA algorithm. As expected, this led to a significant reorganisation of the route and an increase in the total distance travelled but with the benefit of reaching those in need more quickly. The effect of these modifications are show in \ref{fig:Modifications}.

\begin{figure*}
	\begin{subfigure}{0.5\textwidth}
		\includegraphics[width=\textwidth]{figures/Multiple-salesmen}
		\caption{Effect of adding multiple maintenance crews.}
		\label{fig:Multiple-crews}
	\end{subfigure}
	\begin{subfigure}{0.5\textwidth}
		\includegraphics[width=\textwidth]{figures/Prioritising-remote-pumps}
		\caption{Optimised maintenance route for prioritising more remotes pumps.}
		\label{fig:Priority-remote}
	\end{subfigure}
	\caption{Adding multiple maintenance crews and prioritising remote pumps.}
	\label{fig:Modifications}
\end{figure*}

It would be interesting in future work on this data to look at solving the Knapsack Problem to simulate the maintenance crew having a constrained amount of materials to allocate to each pump repair and also the selection of the optimal position of the depot. Another interesting exercise would be to look at how alternative solutions for the TSP compare in terms of both processing time and distance -- for example Genetic Algorithms have been used successfully to solve the TSP.

\section*{Conclusions}

We have shown that Bayesian modelling can produce predictions for the operating status of water pumps that are significantly better than the naive baseline prediction. The approach assumes an a priori model of how the pump functions based on parameters that we believe to be significant and then updates the relative strength of each of these parameters probabilistically to converge towards to optimal values. However, we have also shown that other Machine Learning Methods can produce superior predictions based on simpler assumptions. Having developed methods for identifying faulty pumps we have demonstrated the use of simulated annealing to optimise the route that a maintenance crew can take and added constraints such as a fixed depot, and multiple crews. Interesting avenues for future work have also been suggested.

\bibliography{Write-up_Ref}

\end{multicols}
\end{document}

ext: tex


text:
% Conclusions

\chapter{Conclusions}

Check Final Year Project Guide for the content of Conclusions chapter.

\clearpage %force the next chapter to start on a new page. Keep that as the last line of your chapter!

max_stars_repo_name: visuve/metropolia-thesis-latex

max_issues_repo_name: visuve/metropolia-thesis-latex

max_forks_repo_name: visuve/metropolia-thesis-latex


text:
%!TEX root = report.tex

\section{Methods}\label{sec:methods}

\subsection{Environment}\label{sec:environment}
\begin{figure}[h]
    \centering
    \includegraphics[width=\linewidth]{img/Simulation.png}
    \caption{A visual representation of the environment. The agent is shown in white, leaving behind a trail of inflammable dug cells (shown in brown). The trees (green) can ignite to become burning cells (red), which heat up the neighbouring cells (yellow). When a burning cell runs out of fuel it dies (black)}
    \label{fig:simulation}
\end{figure}
The forest fire is simulated using a grid of cells, each of which has a number of attributes. The main ones are heat potential, temperature, ignition threshold and fuel. Heat potential is the amount of heat, once ignited, the cell radiates to its neighbours to increase their temperature. As soon as the temperature of a cell reaches the ignition threshold, the cell ignites and burns. A burning cell consumes one fuel per iteration, and stops burning and radiating heat when the fuel is empty, becoming a dead cell. The heat from a burning cell reaches the cells directly north, south, east and west of it. If that cell is flammable its temperature is increased by the heat potential of the burning cell, otherwise nothing happens.

For a visual representation please see Figure \ref{fig:simulation}. The shape of the grid is always square and has a size of either 10-by-10 or 14-by-14 cells. The green cells represent trees that can be ignited. The agent (or bulldozer) is represented by the white tile. Wherever it moves it destroys the trees and an empty, inflammable (brown) cell is formed. A line of these dug cells forms a fire line over which the fire cannot spread. Dead cells are represented in black. The agent has to move each time step, it is not allowed to idle on the same cell, and it can only die if it moves into an actively burning (red) cell. The only possible actions for the agent to take are the 4 valid movements (north, south, east, west), and the agent is always digging as it moves. The environment reaches a terminal state and restarts if the agent dies, or if there are no more burning cells.

\subsubsection{The Reward Function}\label{sec:reward_function}
Any reinforcement learning algorithm will rely on the two things that constitute the input to the agent. The quality of the state representation, or in other words how much of and how well the agent can perceive the environment, and the quality of the reward function, which defines how well the agent's notion of success corresponds with the desired behaviour. Because of this, the performance of the agent may depend heavily on the design of the reward function.

The reward function also determines the speed at which the agent will be able to learn the optimal policy. To take the gradient descent analogy of a problem landscape, if the reward function produces a smooth gradient to the optimal solution, the agent will be able to find a path to that solution more easily than if the reward resembles a flat landscape with sparse spikes in which the value jumps from almost always 0 or negative to a positive reward \citep{sutton_barto_2018}. In other words, the agent should be provided gradual feedback instead of sparse and delayed rewards in order to facilitate fast and efficient learning.

We defined the reward function as
\begin{equation}\label{eq:reward_function}
    R_t = \left\{\begin{array}{lr}
        1000, & \text{ Fire is contained }\\
        1000 * (p), & \text{ Fire burns out }\\
        -1000, & \text{ Agent dies }\\
        -1, & \text{ Otherwise }

        \end{array}\right\},
\end{equation}
where $p$ is the percent of the map undamaged by either fire or digging. Note that this reward function does not produce a smooth gradient, and the performance of the agent will likely suffer as a result.

Containment of the fire is defined as
\begin{equation}\label{eq:containment}
    \sum_{f \in \mathcal{F}} \sum_{b \in \mathcal{B}} astar(f, b) = 0,
\end{equation}
where $f$ is a burning cell from the set of currently burning cells $\mathcal{F}$ and $b$ is a cell on the border of the map from the fixed set of border cells $\mathcal{B}$. The function $astar$ is defined as
\begin{equation}\label{eq:astar}
    astar(f, b) = \left\{
        \begin{array}{lr}
            1, & \text{ if A* path exists }\\
            0, & \text{ Otherwise }
        \end{array}\right\},
\end{equation}
where a path is a sequence of directly connecting cells starting with cell $f$ and ending with cell $b$, determined using the A* path-finding algorithm and not allowing diagonal steps. The intuition is that if there exists a path between any burning cell and any cell on the border of the map, then there exists a way for the fire to spread beyond control and so it is not contained.

\subsubsection{The State Representation}\label{sec:state_rep}
The state of the environment, as it is visible to the agent, consists of 3 matrices of size $N^2$ with a boolean domain resulting in, after flattening, an array of $3 * N^2$ boolean inputs where $N$ is the size of the square map. The environment is thus represented three times: One layer contains only the agent position, translating to a matrix of zeros except for a single one representing the position of the agent. The second layer consists of the positions of the fire. Cells that are on fire are represented by a 1, the rest are set to 0. The third layer represents the fire lines cut by the agent in a similar boolean fashion, resulting in a total of 300 inputs to the agent when $N=10$. This is shown graphically in Figure \ref{fig:visiongrid}.

\begin{figure}[h]
    \centering
    \includegraphics[width=1\linewidth]{img/State_Representation.png}
    \caption{An example of the state representation. Each layer shows an important aspect of the world: The location of the fire, inflammable cells and the agent itself.}
    \label{fig:visiongrid}
\end{figure}

This state representation can speed up the learning process as well as increase performance \citep{knegt2018opponent}. Indeed it had a noticeable effect on the performance and learning speed of our implementation compared to a single matrix representing the gray-scaled map as input, likely because the agent can more easily differentiate between important attributes due to the boolean domain, and because only the relevant information is presented. Furthermore, the agent can now see whether the cell it is occupying is already dug or not.

\subsection{RL Algorithm}\label{sec:agent}

\subsubsection{Known Problems with Connectionist RL}\label{sec:problems}
The combination of function approximation (the neural network), bootstrapping (TD methods that update $Q$-values using estimated return values) and off-policy training (the $Q$-Learning algorithm) make up the deadly triad \citep{sutton_barto_2018}, which is known to cause instabilities and even divergence in the learning process. This is at least partly due to correlations in the sequence of consecutive observations, the fact that even small updates to $Q$ can change the policy and therefore the data distribution, and the correlations between the $Q$-values $Q(s_t, a_t; \theta_t)$ and the target values $Y_t = r_{t+1} + \gamma \max_a Q(s_{t+1}, a; \theta_t)$ \citep{mnih2015human}.

As will be clear by the end of this section, all three of these elements of the deadly triad are present in some of the algorithms. There are, however, two well known strategies to reduce this effect and stabilise learning: Using experience replay and a target network.

\subsubsection{Network Architecture}\label{sec:architecture}
First we define the neural network architecture\footnote{The source code is available at https://github.com/dashdeckers/Wildfire-Control-Python and can be accessed and used for non-commercial uses only.} that will approximate $Q$, or the action-value function. We use only a shallow network with one hidden layer of 50 units using a sigmoidal activation function, and a simple output layer with a linear activation. The network takes as an input a state $s$, consisting of an array of size $3*N^2$, and outputs an array of size 4 corresponding to the $Q$-values of taking any of the 4 possible actions $a$ in that state. We used the Adam optimizer to train this network.

\subsubsection{Experience Replay}\label{sec:exp_replay}
\begin{figure}[h]
    \centering
    \includegraphics[width=1\linewidth]{img/Experience_Replay.png}
    \caption{Schematic structure of the experience replay process}
    \label{fig:expreplay}
\end{figure}
Instead of training the network on incoming experiences $e_t=(s_t,a_t,r_t,s_{t+1})$ at each time step $t$ directly, the experiences can be stored in a memory buffer $B_t=\{e_1,...e_t\}$ and the network can be trained on a mini-batch of memories randomly sampled from this buffer $(s,a,r,s') \sim U(B)$ (see Figure \ref{fig:expreplay}). The buffer $B$ allows for $M$ experiences to be stored until it is full, at which point the oldest memories will be discarded to make room for new ones.

This helps because neural networks assume that each training sample is identically and independently distributed from the population, which is in this case the set of $Q$-values that the network is set up to approximate, and training directly on incoming experiences violates this assumption in two regards.

For one, consecutive incoming experiences are obviously highly correlated because the environment does not radically change after any action (unless it is a terminal state). This means that any experience differs from the previous experience only to the extent to which the environment can change in one update step from a single action from the agent.

Secondly, the distribution of incoming experiences is also dependent on the agent's current policy. If the current policy determines that the agent should head east, then the next experiences recorded by the agent will involve the agent headed east. Apart from the obvious correlation, this can also lead to feedback loops and the network parameters could get stuck in local optima or even diverge \citep{tsitsiklis1997analysis,mnih2015human}.

In addition to stabilising the learning process, this approach also increases sample efficiency by allowing memories to be re-used multiple times until they are discarded when they reach the end of the queue. We can also pre-initialize this buffer with memories before learning to provide the agent with demonstration data, as will be explained in Section \ref{sec:demo_data}.

\subsubsection{Target Network}\label{sec:target_network}
Another source of instability arises when we use the predictions of the network to generate the target values which directly update the weights of that same network, as in the standard $Q$-learning update
\begin{equation}\label{eq:update_rule}
    \theta_{t+1} = \theta_t + \alpha(Y_t - Q(s_t, a_t; \theta_t))\nabla_{\theta_t}Q(s_t, a_t; \theta_t),
\end{equation}
where $\alpha$ is the learning rate and the target $Y_t$ is defined as
\begin{equation}\label{eq:target}
    Y_t \equiv r_t + \gamma \max_a Q(s_{t+1}, a; \theta_t).
\end{equation}

The update rule in Equation \eqref{eq:update_rule} resembles stochastic gradient descent, updating the value $Q(s_t, a_t; \theta_t)$ towards the target value $Y_t$. Using the same network can lead to unwanted feedback loops. A delay to the loop can be added to reduce these effects by using a periodically updated, frozen copy of the network with weights $\theta^-$ to generate the target values instead (as shown schematically in Figure \ref{fig:targetnet})
\begin{equation}\label{eq:target_target}
    Y^{Target}_t \equiv r_{t+1} + \gamma \max_a Q(s_{t+1}, a; \theta_t^-),
\end{equation}
where every $C$ time steps, the weights of the network are copied into the target network $\theta_t^- = \theta_t$. This modification was also used by \cite{mnih2015human} alongside experience replay to make the algorithm more stable and dramatically increase its performance.

\begin{figure}[h]
    \centering
    \includegraphics[width=1\linewidth]{img/Target_Network.png}
    \caption{Using a target network to add a delay to the feedback loop. The main network affects (copies itself into) the target network only every $C$ iterations.}
    \label{fig:targetnet}
\end{figure}

\subsubsection{Q-Learning versus SARSA}\label{sec:ql_sarsa}
We implemented both $Q$-Learning \citep{watkins1989learning} and SARSA \citep{rummery1994line} to investigate the difference in performance when using an off-policy versus on-policy algorithm.

There are two commonly used policies, the greedy and the $\epsilon$-greedy policy, and both depend on $Q$. The greedy policy simply always chooses the best action to take in the current state at any given time step: $\argmax_a Q(s_t, a)$.

The $\epsilon$-greedy policy is similar, but it regulates the exploration/exploitation trade-off with the parameter $\epsilon \in [0,1]$. With a probability of $\epsilon$, this policy chooses a random action and otherwise chooses the best action. By annealing this probability towards 0 over time, we can ensure a high exploration rate in the beginning and then slowly decrease it to switch focus to the exploitation stage.

In $Q$-learning the agent follows the $\epsilon$-greedy policy while optimizing the greedy policy, making it an off-policy algorithm. Putting together everything so far, we end up with the following loss function for the base $Q$-Network algorithm:
\begin{equation}
    L_i(\theta_i)= \sum_{(s,a,r,s') \sim U(B)}[(Y^{Target}-Q(s,a,\theta_i))^2].
\end{equation}
Where $Y^{Target}$ is the target value defined in Equation \eqref{eq:target_target}.

In SARSA however, as the acronym suggests, we also save the next action to the memory buffer with the tuple $(s,a,r,s',a')$. This way the agent can both follow and optimize for the same $\epsilon$-greedy policy, making this an on-policy algorithm with the loss function:
\begin{equation}
    L_i(\theta_i)= \sum_{(s,a,r,s',a') \sim U(B)}[(Y^{SARSA}-Q(s,a,\theta_i))^2].
\end{equation}
with
\begin{equation}\label{eq:target_sarsa}
    Y^{SARSA}_t \equiv r_{t+1} + \gamma Q(s_{t+1}, a_{t+1}; \theta_t^-).
\end{equation}
This should lead to an increase in learning stability due to one less element of the deadly triad present in the system. Because SARSA updates $Q$ in a way that takes into account the random actions the agent sometimes takes it is expected to find a safer, but less optimal policy \citep{sutton_barto_2018} and report a higher average performance during learning.

\subsubsection{Dueling Networks}\label{sec:dueling}
One limitation of the $Q$-Network approach is that it is not able to estimate the value of a state and the action separately \citep{wang2015dueling}. This ability can be very useful, since often in states an action has no relevant consequence. The dueling network architecture achieves this by having two streams each predicting different things, instead of one stream predicting both. It is implemented by configuring two hidden layers in parallel, replacing the single hidden layer that is standard in the Q-Network. One of these two layers will estimate the state value $V$, while the other layer will estimate the action advantages $A$ for all possible actions. The two are merged together according to Equation \eqref{eq:dueling}. This equation differs from the original paper, since here we do not use any convolutional layers.

\begin{equation} \label{eq:dueling}
  \begin{array}{l}
    Q(s,a; \alpha, \beta) = V(s; \beta) + \\
    \Big(A(s, a; \alpha) - \frac{1}{c} \sum\limits_{a'} A(s,a'; \alpha)\Big)
  \end{array}
\end{equation}

Here, $\alpha$ and $\beta$ denote the weights of the two fully connected layers and $c$ the number of actions. This equation automatically prevents the layer for the state value from estimating anything related to the action advantages, since the sum of the advantages is kept at zero. The equation results in the $Q$-values which can be used in the same way as the single stream Q-values. All learning techniques and code can be recycled, and the target network and experience replay is also used. Given the results from the original paper, this modification is expected to boost the maximum performance and especially the speed of learning.

\subsubsection{Dueling SARSA}\label{sec:duelingsarsa}
As well as combining dueling networks with $Q$-Learning, we also combine on-policy SARSA with the dueling network architecture to introduce a new RL algorithm, Dueling SARSA. We expect this algorithm to benefit from both modifications for a twice improved performance and both increased learning speed and learning stability. Just as the other three algorithms, Dueling SARSA also uses the target network and experience replay.

\subsubsection{Demonstration Data}\label{sec:demo_data}
Because the reward function defined in Section \ref{sec:reward_function} provides only sparse and delayed rewards, the agent might require additional guidance to learn to contain the fire in a reasonable training time. To point it in the right direction, we can fill the memory buffer with demonstration data before learning as mentioned in Section \ref{sec:exp_replay}. The agent can then use these memories during the learning process to its advantage.

The demonstration data only needs to show the agent how it might be able to collect the containment reward, to this effect we developed a simple algorithm that moves the agent around the fire in a clockwise direction. It does this by choosing randomly from one of two possible actions, unless one of the actions would lead to the death of the agent in which case it chooses the safe action. The two possible actions depend on the position of the agent relative to the fire (This is shown schematically in Figure \ref{fig:demodata}).

The environment is reset upon containment as defined in Equation \eqref{eq:containment}. Only memories (transitions) leading to successful containments are stored. This results in an average of 35 memories per episode, or containments of the fire, on the 10x10 map, and around 48 memories per episode on the 14x14 map. The number of episodes of demonstration data required can be specified.

\begin{figure}[h]
    \centering
    \includegraphics[width=1\linewidth]{img/Demo-data_Baseline.png}
    \caption{The possible actions for the agent to take, based on the position (quadrant) of the agent relative to the origin of the fire}
    \label{fig:demodata}
\end{figure}

\subsection{Experimental Setup}\label{sec:experiment}

\subsubsection{Baseline Algorithm}\label{sec:baseline}
To be able to reliably compare the performance of the different algorithms, we need a stable baseline. We build on the ideas discussed in Section \ref{sec:demo_data} to define the algorithm shown in Algorithm \ref{alg:baseline}. This algorithm is identical to the demonstration data generation except that it continues until the fire has burnt out, so it does not stop as soon as the fire is contained.

\begin{algorithm}
  \caption{Baseline algorithm to contain the fire}
  \label{alg:baseline}
  \begin{algorithmic}[1]
    \Procedure{RunBaseline}{}
    \State {totalreward = 0}
    \While {\textbf{not} done}
    \State {action = random(possible actions)}
    \If {action is dangerous}
    \State {action = other possible action}
    \EndIf
    \State {reward, done = execute(action)}
    \State {totalreward = totalreward + reward}
    \EndWhile
    \State \Return {totalreward}
    \EndProcedure
  \end{algorithmic}
\end{algorithm}

\subsubsection{Collection of Results}\label{sec:datacollection}
We investigate the performance of both $Q$-Learning and SARSA, both with and without the dueling network architecture modification, for a total of 4 different algorithms. The hyper-parameters we used throughout all experiments are shown in Table \ref{tab:hyperparameters}.

Each algorithm and parameter combination was run 10 times for 10,000 episodes per run. We compared the baseline performance to 2 map sizes ($N=10$ and $N=14$), 3 different amounts of demonstration data (0, 100 and 1000 episodes) and 4 algorithms for a total of 240 runs (260 including the baselines). At the start of each run, the environment is initialized with trees and a single cell at the center of the map is ignited. The agent starts at a random location on a circle centered around the fire with a radius of either 1, 2 or 3 (also randomly selected). Since each run took approximately 4.5 hours on average, they were calculated on the Peregrine computing cluster provided by the University of Groningen.

\begin{table}[H]
\centering
\begin{tabular}{|l|l|}
\hline
Memory size & 20000 \\ \hline
Batch size & 32 \\ \hline
Target update ($C$) & 20 \\ \hline
Gamma ($\gamma$) & 0.999 \\ \hline
Alpha ($\alpha$) & 0.005 \\ \hline
Epsilon decay ($\epsilon$) & 0.01 \\ \hline
Epsilon maximum & 1 \\ \hline
Epsilon minimum & 0.005 \\ \hline
\end{tabular}
\caption{All relevant hyper-parameters used in the training process. These values were selected by performing an informal search using the Q-Learning algorithm without dueling networks. The target is updated every $C$ episodes. The epsilon value is decayed after every episode.}
\label{tab:hyperparameters}
\end{table}

max_stars_repo_name: dashdeckers/Wildfire-Control-Python

max_issues_repo_name: dashdeckers/Wildfire-Control-Python

max_forks_repo_name: dashdeckers/Wildfire-Control-Python


text:
% $Id$
%
% Author: David Fournier
% Copyright (c) 2008 Regents of the University of California
%

This chapter attempts to list and document all the functions
available in \ADM. It will always be incomplete, since functions are
continually being added. If you are aware of a function that is
not documented please, contact me at \href{mailto:users@admb-project.org}{ADMB Users Group} and
let me know.

\section{Naming conventions for documenting functions}

Wherever applicable, the \texttt{name} function has been supplied for
constant and variable objects (such as \texttt{double} and \texttt{dvariable}).
Instead of repeating the description for both kinds of objects,
the convention of referring to both types as ``number,'' ``vector,''
``matrix,'' etc., will be observed.

\section{Mathematical functions}

 The following functions have been included in \scAD, by overloading the
\cplus\ library functions or adding additional functions where necessary:
\begin{lstlisting}
  acos atan  cos cosh cube exp (mfexp) fabs gammln (sfabs) log log_comb
  log10 log_density_poisson pow square sqrt sin sinh tan tanh
\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}}
\X{\fontindexentry{tt}{gammln}}\X{\fontindexentry{tt}{log\_comb}}
\XX{functions}{\fontindexentry{tt}{sin}} \XX{functions}{\fontindexentry{tt}{cos}} \XX{functions}{\fontindexentry{tt}{tan}}
\XX{functions}{\fontindexentry{tt}{asin}} \XX{functions}{\fontindexentry{tt}{atan}} \XX{functions}{\fontindexentry{tt}{acos}}
\XX{functions}{\fontindexentry{tt}{sinh}} \XX{functions}{\fontindexentry{tt}{cosh}} \XX{functions}{\fontindexentry{tt}{tanh}}
\XX{functions}{\fontindexentry{tt}{fabs}} \XX{functions}{\fontindexentry{tt}{sfabs}} \XX{functions}{\fontindexentry{tt}{exp}}
\XX{functions}{\fontindexentry{tt}{log}} \XX{functions}{\fontindexentry{tt}{log10}} \XX{functions}{\fontindexentry{tt}{sqrt}}
\XX{functions}{\fontindexentry{tt}{pow}}
\XX{functions}{\fontindexentry{tt}{gammln}}\XX{functions}{\fontindexentry{tt}{log\_comb}}
These functions can be used on \texttt{number}s or \texttt{vector\_object}s
in the form
\begin{lstlisting}
number = function(number);
vector_object = function(vector_object);
\end{lstlisting}
When operating on \texttt{vector\_object}s, the functions operate
elemen-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}\XX{\fontindexentry{tt}{max}}{operation on a vector}
\XX{vector}{minimum element}\XX{\fontindexentry{tt}{min}}{operation on a vector}
The functions \texttt{min} and \texttt{max}, when applied to a \texttt{vector\_object},
return a \texttt{number} that is equal to the minimum or maximum element of the
\texttt{vector\_object}.

The function \texttt{gammln} is the logarithm of the gamma function.
\XX{gamma function}{logarithm}

The function \texttt{log\_comb(n,k)} is the logarithm of the  function,
the combination of \texttt{n} things taken \texttt{k} at a time.
It is defined via the logarithm of the gamma function
for non-integer values,
and is differentiable.

\X{\fontindexentry{tt}{log\_comb} function}

\section{Operations on arrays}

\subsection{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. Spreadsheet users might find it much more
natural to define the product of matrices as an element-wise operation,
such as
$$z_{ij}= x_{ij}*y_{ij}$$
The ``classical'' mathematical definition for the
matrix product has been assigned to the overloaded 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{functions}{element-wise product of vectors}
\XX{vector}{element-wise product}\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}
\XX{functions}{element-wise division of vectors}
\begin{lstlisting}
vector_object = elem_div(vector_object,vector_object) // element-wise divide
\end{lstlisting}
\afterlistingthing{$z_i= x_i/y_i $}

\XX{functions}{element-wise product of matrices}
\XX{matrix}{element-wise product}
\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}
\XX{functions}{element-wise division of matrices}
\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{functions}{identity matrix function}
\XX{matrix}{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{Probability densities and related functions:\BR \texttt{poisson negative binomial cauchy}}
\XX{functions}{Cauchy density}
\XX{probability densities}{Cauchy}

\begin{lstlisting}
number log_density_cauchy(number x);
\end{lstlisting}
\afterlistingthing{returns the logarithm of the Cauchy density function
at \texttt{x}.}

\XX{functions}{Poisson density}
\XX{probability densities}{Poisson}
\begin{lstlisting}
number log_density_poisson(number x,number mu);
\end{lstlisting}
\afterlistingthing{returns the logarithm of the Poisson density function
at \texttt{x}
with mean \texttt{mu}.}

\XX{functions}{negative binomial density}
\XX{probability densities}{negative binomial density}
\begin{lstlisting}
number log_negbinomial_density(number x,number mu,number tau);
\end{lstlisting}
\afterlistingthing{returns the logarithm of the negative  binomial density function
with mean \texttt{mu} and over-dispersion (variance/mean) \texttt{tau}.
\texttt{tau} must be greater than~1.}

\section{The operations \texttt{det inv norm norm2 min max sum}}

\subsubsection{The determinant of a matrix object}
\XX{matrix}{determinant}\X{determinant}
\XX{functions}{determinant of a matrix}
\X{\fontindexentry{tt}{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}
\XX{matrix}{inverse}\X{inverse} \X{\fontindexentry{tt}{inv}}
\XX{functions}{inverse of a matrix}

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

\subsubsection{The norm of a vector object}
\XX{functions}{norm of a vector}
\XX{vector}{\fontindexentry{tt}{norm}}

\begin{lstlisting}
number = norm(vector_object)
\end{lstlisting}
\afterlistingthing{$z=\sqrt{\sum_{i} x_{i}^2} $}
\bigskip %zzz

\subsubsection{The norm squared of a vector object}
\XX{vector}{norm squared}\X{\fontindexentry{tt}{norm2}}
\XX{matrix}{norm squared}\X{\fontindexentry{tt}{norm2}}
\XX{functions}{norm squared of a vector}

\begin{lstlisting}
number = norm2(vector_object)
\end{lstlisting}
\afterlistingthing{$z=\sum_{i} x_{i}^2 $}
\bigskip

\subsubsection{The norm of a matrix object}
\XX{matrix}{\fontindexentry{tt}{norm}}\X{\fontindexentry{tt}{norm}}
\XX{functions}{norm of a matrix}

\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}
\XX{functions}{norm squared of a matrix}
\XX{matrix}{norm squared}\X{\fontindexentry{tt}{norm2}}

\begin{lstlisting}
number = norm2(matrix_object)
\end{lstlisting}
\afterlistingthing{$z_{ij}=x_{ji} $}
\bigskip

\subsubsection{The transpose of a matrix object}
\XX{functions}{transpose of a matrix}
\XX{matrix}{transpose}

\begin{lstlisting}
matrix_object = trans(matrix_object)
\end{lstlisting}
\afterlistingthing{$z=\sum_{ij} x_{ij}^2 $}
\bigskip

\subsubsection{The sum over the elements of a vector object}
\XX{functions}{sum over the elements of a vector}
\XX{vector}{sum over the elements}\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}}\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}}\XX{\fontindexentry{tt}{colsum}}{operation on a matrix}
\XX{functions}{sum over the columns of 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}\XX{\fontindexentry{tt}{min}}{operation on a vector}
\XX{functions}{minimum element of a vector}

\begin{lstlisting}
number = min(vector_object)
\end{lstlisting}
\bigskip

\subsubsection{The maximum element of a vector object}
\XX{vector}{maximum element}\XX{\fontindexentry{tt}{max}}{operation on a vector}
\XX{functions}{maximum element of a vector}

\begin{lstlisting}
number = max(vector_object)
\end{lstlisting}
\bigskip

\section{Eigenvalues and eigenvectors of a symmetric matrix}

\XX{eigenvalues}{not differentiable}
\XX{eigenvectors}{not differentiable}
\XX{functions}{eigenvalues of symmetric matrix}
\XX{functions}{eigenvectors of symmetric matrix}

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.

\XX{eigenvectors} {of a symmetric matrix}
\begin{lstlisting}
matrix_object = eigenvectors(matrix_object)
\end{lstlisting}
\afterlistingthing{returns in a matrix the eigenvectors of a symmetric matrix.}

 It is the user's responsibility to
ensure that the matrix is actually symmetric. The routine symmetrizes
the matrix, so 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{eigenvectors}.

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

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{solve@\texttt{solve} function}
\XX{functions}{solving linear system of equations}
\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} while solving the system of linear
equations.  Since this is useful in multivariate
analysis, we have also included a function that returns the
determinant when 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}

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 either 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{functions}{filling arrays and matrices}
\X{filling arrays and matrices}
\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, 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, base+offset, 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, base+offset, 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, base+offset, 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 $j^\textrm{th}$ column of a matrix with a vector.}

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

\section{Methods for filling arrays and matrices\br with random numbers}

This method of filling containers with random numbers is becoming
obsolete. The preferred method is to use the
\texttt{random\_number\_generator} class.
See Section~\ref{rng} for instructions on using this class.
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 {\it 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.
{\it This function is now obsolete. You should use the}
\texttt{random\_number\_generator} \textit{class to generate 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 of 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.
{\it This function is now obsolete. You should use the}
\texttt{random\_number\_generator} \textit{class to generate 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~$p$ is a \texttt{dvector} such that
\texttt{p[i]} is the probability of returning~$i$. The elements of \texttt{p}
must sum to~$1$.}

\section{Methods for obtaining shape information\br from containers}

When this code was first written, the maximum dimension of arrays
was about four. At this level, it perhaps make sense to think of
a 1-dimensional array as a vector, a 2-dimensional array as
a matrix, etc. For a matrix, one thinks in terms of rows and columns.
However, with the adoption of ragged container objects up to
eight dimensions (at present), a more generic method of obtaining
shape information of these objects was called for.

If \texttt{v} is a vector object, then
\begin{lstlisting}
int v.indexmin()
int v.indexmax()
\end{lstlisting}
return the minimum and maximum valid indices for \texttt{v}.
If \texttt{m} is a matrix object, then
\begin{lstlisting}
int v.rowmin()
int v.rowmax()
int v.colmin()
int v.colmax()
\end{lstlisting}
return the minimum and maximum valid row and column indices
for \texttt{m}. These functions make sense for a matrix where every
row is a vector with the same minimum and maximum valid
indices.  For a ragged matrix, this is no longer the case, so
the \texttt{rowmin()} and \texttt{rowmax()} functions don't make sense.
To deal with a ragged matrix, one may need to calculate
the minimum and maximum valid indices for each row of the ragged matrix.
To facilitate this approach, the functions \texttt{indexmin} and
\texttt{indexmax} have been defined for all container classes.  So,
for example, if \texttt{w} is a 6-dimensional array, then
\begin{lstlisting}
int w.indexmin()
int w.indexmax()
\end{lstlisting}
return the minimum and maximum valid indices for the first
index of \texttt{w}. For a matrix object \texttt{m},
\texttt{m.indexmin()} and \texttt{m.colmin()} are the same
and as long as \texttt{m} is not ragged,
\begin{lstlisting}
  m(m.indexmin()).indexmin()
\end{lstlisting}
is the same as
\begin{lstlisting}
  m.colmin()
\end{lstlisting}
and
\begin{lstlisting}
  m(m.indexmin()).indexmax()
\end{lstlisting}
is the same as
\begin{lstlisting}
  m.colmax()
\end{lstlisting}

\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}
\XX{matrix}{\fontindexentry{tt}{column}}
\begin{lstlisting}
vector column(matrix& M,int& j)
\end{lstlisting}
\afterlistingthing{extracts the $j^\textrm{th}$ 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.
Note that the operation \texttt{M(i)} has the same effect.}

\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 \fontindexentry{tt}{()}}
The function call operator \texttt{()} has been overloaded in two ways to
provide for the extraction of a subvector.
\XX{\fontindexentry{tt}{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 \fontindexentry{tt}{()} 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},
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 \fontindexentry{tt}{++}}
\X{operator\fontindexentry{tt}{-{}-}}
\XX{operator \fontindexentry{tt}{++}}{use with subvectors}
\XX{operator \fontindexentry{tt}{-{}-}}{use with subvectors}
\XX{operator \fontindexentry{tt}{++}}{for \fontindexentry{tt}{dvectors}}
\XX{operator \fontindexentry{tt}{-{}-}}{for \fontindexentry{tt}{dvectors}}
\XX{functions}{\fontindexentry{tt}{++} use with vectors}
\XX{functions}{\fontindexentry{tt}{-{}-} use with vectors}
\XX{functions}{\fontindexentry{tt}{shift} use with vectors}
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{-{}-}
increment and decrement the index ranges by~1.
\begin{lstlisting}
dvector  u(1,20);
dvector  w(1,19);
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
               // probably not what you want
w=v;   // error different index ranges
\end{lstlisting}
It is important to realize that from the point of view of the vector
\texttt{u}, both of the above assignments have the same effect. It will have
elements~2 through~20 set equal to~2. The difference is in the side
effects on the vector \texttt{v}. The operation \texttt{++v}
will increase the minimum and maximum valid index range of the vector
\texttt{v} by~1. This increase is permanent. On the other hand,
the operation \texttt{-{}-u(2,20)} decrements the valid index bounds of
the {\it subvector} \texttt{u(2,20)}. This is a distinct object from
the vector \texttt{u}, although both objects share a common area for their
components. Thus, the valid index bounds of \texttt{u} are not
affected by this process.
The use of subvectors, along with 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}

\XX{functions}{\fontindexentry{tt}{shift} use with vectors}
The \texttt{shift} function can be used to set the minimum (and maximum)
valid index for a vector.
\begin{lstlisting}
dvector u(10,100);  // minimum valid index is 10
                    // maximum valid index is 100
u.shift(25);        // minimum valid index is 25
                    // maximum valid index is 115
\end{lstlisting}
In particular, the operators \texttt{-{}-} and \texttt{++}
are just convenient shorthand for using the \texttt{shift}
function to change the minimum valid index by~1.
\begin{lstlisting}
dvector u(10,100);     // minimum valid index is 10
                       // maximum valid index is 100
u.shift(u.indexmin()-1); // minimum valid index is 9
--u;                   // same result as u.shift(u.indexmin()-1)
u.shift(u.indexmin()+1); // minimum valid index is 11
++u;                   // same result as u.shift(u.indexmin()+1)
\end{lstlisting}

\XX{subobjects}{accessing with\fontindexentry{tt}{sub}}
\XX{\fontindexentry{tt}{sub}}{accessing subobjects}

\section{Accessing subobjects of higher-dimensional arrays}

The  \texttt{()} operator cannot be used to access subobjects of
arrays of dimension~2 or greater, because this operator has
already been defined to do something else. For example,
for a \texttt{dmatrix M}, \texttt{M(1,2)} is an {\it element}
of \texttt{M}.  To access subobjects of higher-dimensional arrays,
use the \texttt{sub} member function. If \texttt{M} is a matrix object,
then \texttt{M.sub(2,6)} is a matrix object with minimum valid
index~2 and maximum valid index~6 (provided, of course, that the
minimum valid index for \texttt{M} is less than or equal to~2 and
the maximum valid index is greater than or equal to~6).
If \texttt{T} is a 3-dimensional object, then
\texttt{T.sub(2,5)} is a 3-dimensional object, provided that
the index bounds are legal.

\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 a part of methods for calculating the
derivatives of differentiable functions (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
by using the \texttt{()} operator.
\X{operator \fontindexentry{tt}{()}}
\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 order.

\section{Statistical functions}

\begin{lstlisting}
cumd_norm
inv_cumd_norm
cumd_cauchy
inv_cumd_cauchy
\end{lstlisting}
These are the cumulative distribution function and the inverse cumulative distribution
function for the normal and Cauchy distributions.

\section{The random number generator class}\label{rng}
\XX{cumulative distribution function}{normal}
\XX{cumulative distribution function}{Cauchy}
\XX{\fontindexentry{tt}{random\_number\_generator}}{class}
\XX{function}{\fontindexentry{tt}{random\_number\_generator}}

The random number generator class is used to manage the generation of
random numbers. A random number generator object is created with the
declaration
\begin{lstlisting}
  random_number_generator r(n);
\end{lstlisting}
where \texttt{n} is the seed that initializes the random number generator.
Any number of random number generators may be declared.
This class can be used to manage random number generation with the following functions:
\begin{lstlisting}
randpoisson(lambda,r);  // generate a Poisson with mean lambda
randnegbinomiual(mu,tau,r);  // generate a negative binomial with mean mu
                            // and over-dispersion tau (tau>1)
randn(r);     // generate a normally distributed random number
randu(r);    // generate a uniformly distributed random number
v.fill_randu(r)  // fill a vector v
v.fill_randn(r)  // fill a vector v
v.fill_randpoisson(mu,r)  // fill a vector v with Poisson distributed
                          // random variables with mean mu
v.fill_rand(mu,tau,r)  // fill a vector v with negative binomial distributed
            // random variables with mean mu and over-dispersion var/mu = tau
v.fill_multinomial(r,p)  // fill a vector v
                      // p is a vector of probabilities
m.fill_randu(r)  // fill a matrix m
m.fill_randn(r)  // fill matrix m
m.fill_randpoisson(lambda,r)  // fill a matrix m
\end{lstlisting}
The incomplete beta function $I_x(a,b)$ is defined by
\begin{equation}
  I_x(a,b)=\frac{1}{B(a,b)}\int_0^x t^{a-1}(1-t)^{b-1}\,\textrm{d}t\qquad (a,b,>0)
\end{equation}
This is also the cumulative distribution function for the beta family
of probability distributions.
The function is named \texttt{betai} and is invoked by
\XX{cumulative distribution function}{\fontindexentry{tt}{beta}}
\XX{functions}{\fontindexentry{tt}{beta} function, incomplete}
\begin{lstlisting}
\\ .....
dvariable p=betai(a,b,x);
\end{lstlisting}

\section{The \texttt{adstring} class operations}
\XX{\fontindexentry{tt}{adstring} class}{operations on}
\X{operations on strings}
\X{string operations}

The \texttt{adstring} class was defined before there was a standardized
\cplus\ string class. It does not contain all the features
that a full string class should have. It is, however, easier to use
in many cases than the standard C string operations.
\begin{lstlisting}
adstring s;
adstring t;
s="first_part";
t="second_part";
adstring u = s + " ___ " + t;
cout << u << endl;
\end{lstlisting}
should print out
\begin{lstlisting}
  first_part ___ second_part
\end{lstlisting}
The operation \texttt{+} concatenates two \texttt{adstring} objects.
It can be used to concatenate C-style strings by first
turning them into \texttt{adstring} objects, as in
\begin{lstlisting}
  adstring u = adstring("xxx") + adstring("yyy");
\end{lstlisting}
One can also append to a string with the \texttt{+=}
operator, as in
\begin{lstlisting}
adstring u = "abc";
u += v;
adstring w = "abc";
w += 'f';
\end{lstlisting}
which adds the \texttt{adstring} object \texttt{v} to \texttt{u} and the
character `f' to \texttt{w}.
It is also possible to cast an \texttt{adstring} object to
a C-like \texttt{char *} string, as in
\begin{lstlisting}
adstring u = "abc"
char * c = (char*)(u);
\end{lstlisting}
Then it may be used as you would use a C-like string.

\section{Miscellaneous functions}

\begin{lstlisting}
posfun(x,eps,pen)
\end{lstlisting}
\afterlistingthing{The \texttt{posfun} function constrains the argument \texttt{x} to be positive.
For $\texttt{x} > \texttt{eps}$, it is the identity function.}

The current source code for the \texttt{posfun} function appears below:
\XX{functions}{\fontindexentry{tt}{posfun}}
\X{\fontindexentry{tt}{posfun} function}
\begin{lstlisting}
dvariable posfun(const dvariable&x,const double eps,dvariable& pen)
{
  if (x>=eps) {
    return x;
  } else {
    pen+=.01*square(x-eps);
    return eps/(2-x/eps);
  }
}
\end{lstlisting}

\begin{lstlisting}
mfexp(_CONST prevariable& x)
\end{lstlisting}
\afterlistingthing{The \texttt{mfexp} function is the exponential function
that is modified for large values of its argument, to prevent
floating point overflows.}
The current source code for the \texttt{mfexp} function appears below:
\XX{functions}{\fontindexentry{tt}{mfexp}}
\X{\fontindexentry{tt}{mfexp} function}
\begin{lstlisting}
  dvariable mfexp(_CONST prevariable& x)
  {
    double b=60;
    if (x<b)
    {
      return exp(x);
    }
    else
    {
      return exp(b)*(1.+2.*(x-b))/(1.+x-b);
    }
  }
\end{lstlisting}

max_stars_repo_name: wStockhausen/admb

max_issues_repo_name: wStockhausen/admb

max_forks_repo_name: wStockhausen/admb