Index: /trunk/doc/pslib/psLibADD.tex
===================================================================
--- /trunk/doc/pslib/psLibADD.tex	(revision 306)
+++ /trunk/doc/pslib/psLibADD.tex	(revision 307)
@@ -35,7 +35,187 @@
 \section{PanSTARRS Library PSLib}
 
-\subsection{Data Manipulation}
-
-\subsubsection{Matrix Operations}
+\subsection{Math Utilities}
+
+\subsection{Sorting}
+
+A variety of sorting algorithms exist, with a wide range in speed for
+different set sizes.  Three of the best sorting algorithms are
+``heapsort'', ``quicksort'', and ``mergesort'' (see, e.g.,
+\citealt{knuth}, \citealt{sedgewick}, \citealt{press}).  These three
+sorting algorithms all run in a time of $O(n \log n)$ in the best
+case, but have different worse-case run times.  Implementations of all
+three sorting algorithms are described in references above and in
+various places online (e.g.,
+http://www.iti.fh-flensburg.de/lang/algorithmen/sortieren/algoen.htm).
+The linux \code{qsort} function orignially mades use of the quicksort
+algorithm, but now uses \tbd{what algorithm?}.  The GSL function
+\code{gsl_heapsort} uses the heapsort algorithm.  Both of these
+implemenations require a pointer to the comparison function, which may
+result in a slower code.  
+
+For PSLib, the sorting functions shall be implemented via
+\code{qsort}.  The function \code{psSort} shall return the sorted
+result of the input array without over-writing the input array.  The
+function \code{psSortIndex} shall return an integer index to the
+sequence of the input array without overwriting the input array.
+Given the following line of code, \code{out = psSortIndex (NULL,
+&in);}, the elements of the array \code{out} are in the sequence
+\code{in.arr[out->arr[0]]} to \code{in.arr[out->arr[in.n - 1]]}. 
+
+\subsection{Smoothing: Boxcar and Gaussian}
+\label{smooth}
+
+Smoothing may occasionally be perfomed on data.  We present the
+algorithms for two typical versions: boxcar and Gaussian smoothing.
+In both smoothing techniques, given a series of data values $f_i$, the
+smoothed values $g_i$ are determined by calculating a linear
+combination based on the input data point and its nearest $2N$
+neighbors in the form:
+
+\[ g_i = \sum_{n=-N}^N c_n f_i \]
+%
+where the values of $c_n$ determine the filter type.  For boxcar
+smoothing, the values $c_n$ are constant, and must be equal to $1/(2N
++ 1)$ to maintain the zeroth moment of the data (care must be taken at
+the ends of the data range to reduce the value of $c_n$ as fewer input
+data points may be used).  For Gaussian smoothing, the crucial
+parameter is $\sigma$, the standard deviation.  The
+value of $N$ should be chosen to be large enough, $N = 5\sigma$, and
+the values of $c_n$ are then just the Gaussian curve:
+
+\[ c_n = \frac{e^{\frac{-n^2}{2\sigma^2}}}{\sqrt{2\pi\sigma^2}} \]
+
+\subsection{Statistics}
+
+The general statistics function \code{psStats} performs a variety of
+statistical calculations on a collection of floating point numbers.
+These statistics include both sample values, in which the formal
+statistic is calculated for the complete sample, and robust values, in
+which the statistic is estimated on the basis of an assumption of an
+underlying population.  While more reliable in general, the robust
+estimators may be somewhat slower than the sample statisic, so their
+use may vary depending on the context.  In addition, certain
+sigma-clipped values are defined as an intermediate choice between the
+sample and robust estimators.  
+
+\subsubsection{Sample Statistics}
+
+We define the following statistical terms, assuming there is a set of
+data elements $x_i$.
+
+\paragraph{Mean}
+
+The mean is defined as: \[ \bar{x} = \frac{1}{N} \sum_{i = 1}^{N} x_i \]
+
+\paragraph{Median}
+
+The median is defined as the value for which 50\% of the data values
+are larger and 50\% are smaller.  For a sample, the values are sorted
+and the median is given by the value of the middle element, if the
+number of values is odd, and by the average of the two middle values
+if the number of values is even.  This median should be avoided for
+samples which are large (e.g., $N > 10000$ elements) as the basic
+robust median is quicker and more accurate.  
+
+\paragraph{Upper and Lower Quartiles}
+
+The upper and lower quartiles ($U_{\frac{1}{4}}$ and
+$L_{\frac{1}{4}}$) are first and last quarter equivalents to the
+median.  The upper quartile is the value for which 25\% of the data
+are larger and 75\% are smaller.  The lower quartile is the value for
+which 75\% of the data are larger and 25\% are smaller.  For a sample,
+the values are sorted and the quartiles are given by the values of
+elements $N/4$ and $3N/4$.  For the purpose of a sample upper and
+lower quartile, it is sufficient to provide the closest integer entry
+to these values.  The sample quartiles should be avoided for samples
+which are large (e.g., $N > 10000$ elements) as the robust quartiles
+are quicker and more accurate.
+
+\paragraph{Standard Deviation}
+
+The standard deviation of the sample is given by:
+
+\[ \sigma = \sqrt{\sum_{i = 1}^N \frac{(x_i - \bar{x})^2}{N - 1}} \]
+
+To minimize the numerical rounding error, this should be calculated
+numerically as:
+
+\[ \sigma = \sqrt{\frac{1}{N - 1} [ \sum_{i = 1}^{N} (x_i - \bar{x})^2 - \frac{1}{N} (\sum_{i = 1}^{N} (x_i - \bar{x}))^2 ]} \]
+
+\subsubsection{Clipped Statistics}
+
+The clipped statistics are used to determine the mean and $\sigma$ of
+a distribution in the presence of outliers. The clipped statistics are
+quicker than the complete robust statistical estimators and more
+reliable than the sample statistics.  The clipped statistics algorithm
+produces the clipped mean and clipped standard deviation.  The
+algorithm has 2 parameters: $N$, the number of iterations and $k$, the
+multiplying factor of $\sigma$ used to exclude outliers.  Typical
+values for both $N$ and $k$ are 3.  The algorithm is as follows:
+
+\begin{enumerate}
+\item Compute the sample median. 
+\item Compute the sample standard deviation.
+\item Use the sample median as the first estimator of the mean, $\bar{x}$.
+\item Use the sample standard deviation as the first estimator of the
+  true standard deviation, $\sigma$.
+\item Repeat the following N times:
+  \begin{enumerate}
+  \item Exclude all values $x_{i}$ for while $|x_{i} - \bar{x}| > k \sigma$.
+  \item Compute the mean and standard deviation of the sub-sample.
+  \item Use the new mean for $\bar{x}$.
+  \item Use the new standard deviations for $\sigma$.
+\end{enumerate}
+\item The last calulcated value of $\bar{x}$ is the clipped mean.
+\item The last calulcated value of $\sigma$ is the clipped standard deviation.
+\end{enumerate}
+
+\subsubsection{Robust Statistics}
+
+The robust version of the statistics provide estimators of basic
+statistical concepts which are reliable even for data samples with
+significant contamination.  A typical case is the situation in which
+the data of interest represent a primary population of interest with a
+single-valued mean and standard deviation and a secondary population
+of data with a substantially different distribution.  For example, an
+image of an uncrowded night-time field may consist of a sparse
+collection of stars and an overall background level.  The majority of
+pixels have the background value, with some variance due to noise
+sources, but many are significantly higher (contributed by stars) or
+significantly lower (dead pixels).  If we want to measure the mean of
+the background, a robust mean is necessary as the outliers will
+strongly bias the sample statistics.
+
+The robust statistics are calculated by constructing a histogram of
+the values and performing the measurements on the histogram.  The
+choice of a bin size requires some care.  If the data are integer
+valued, the natural bin size is an integer.  Otherwise, the bin should
+be a fraction of an estimate of the standard deviation.  Use the
+sample upper and lower quartiles to determine an estimate of the
+standard deviation: $\sigma_e = (U_{\frac{1}{4}} - L_{\frac{1}{4}}) /
+1.34$.  The bin size shall be set at $\sigma_e / 10$.  The remaining
+steps of the algorithm are as follows:
+
+\begin{itemize}
+\item Construct the histogram with the specified bin size.
+\item Smooth the histogram by a Gaussian with $\sigma_s = \sigma_e /
+  4$.  
+\item Find the bin with the peak value in the range $L_{\frac{1}{4}}$
+  to $U_{\frac{1}{4}}$; this is the robust mode $\mbox{mode}_r$.  
+\item Determine $dL = (U_{\frac{1}{4}} - L_{\frac{1}{4}}) / 8$.
+\item Fit a Gaussian to the bins in the range $\mbox{mode}_r - dL$ to
+  $\mbox{mode}_r + dL$.
+\item The resulting fit parameters are the robust mean,
+$\mbox{mean}_r$ and the robust standard deviation, $\sigma_r$.
+\end{itemize}
+
+To determine the robust median, construct the cumulative histogram
+from the histogram above. Select the bin which contains the 50th
+percentile value and its two neighbors.  Fit a quadratic to these
+three points.  The robust median value is the coordinate of the
+quadratic which returns the 50\% value.
+
+
+\subsection{Matrix Operations}
 
 In this section, we define the linear algebra operations performed on
@@ -56,5 +236,5 @@
 \code{gsl_linalg_LU_decomp}.
 
-\paragraph{LU Decomposition}
+\subsubsection{LU Decomposition}
 \label{LUdecomp}
 
@@ -91,5 +271,5 @@
 \[ \alpha_{ij} = \frac{1}{\beta_{jj}} \left( a_{ij} - \sum_{k=1}^{j-1} \alpha_{ik}\beta_{kj} \right) \]
 
-\paragraph{Calculate a matrix determinant}
+\subsubsection{Calculate a matrix determinant}
 
 The determinant $D$ of a matrix $a_{ij}$ is calculated from the
@@ -105,5 +285,5 @@
 shall be used.
 
-\paragraph{Solving a Linear Equation}
+\subsubsection{Solving a Linear Equation}
 
 The LU decomposition of a matrix may be used to solve the
@@ -121,5 +301,5 @@
 \[ x_i = \frac{1}{\beta_{ii}}\left(y_i - \sum_{j=i+1}^N \beta_{ij} y_ij\right) \]
 
-\paragraph{Invert a matrix}
+\subsubsection{Invert a matrix}
 
 Inversion of a matrix using the LU decomposition is performed by
@@ -130,5 +310,5 @@
 operation shall be implemented using the GSL function \code{gsl_linalg_LU_invert}.
 
-\paragraph{Perform matrix addition, subtraction and multiplication}
+\subsubsection{Perform matrix addition, subtraction and multiplication}
 
 Matrix binary arithmetic operations differ from image binary
@@ -148,5 +328,5 @@
 operators: $+, -, \times, /$
 
-\paragraph{Transpose a matrix}
+\subsubsection{Transpose a matrix}
 
 The transpose of a matrix is simply the reorganization of the matrix
@@ -159,6 +339,8 @@
 where $M_{ji}$ is the matrix to be transposed.
 
-\paragraph{Convert a matrix to a vector}
+\subsubsection{Convert a matrix to a vector}
 \TBD{write me} 
+
+\subsection{Fitting}
 
 \subsubsection{General Polynomial Fitting}
@@ -255,189 +437,9 @@
 \end{center}
 
-\subsection{Sorting}
-
-A variety of sorting algorithms exist, with a wide range in speed for
-different set sizes.  Three of the best sorting algorithms are
-``heapsort'', ``quicksort'', and ``mergesort'' (see, e.g.,
-\citealt{knuth}, \citealt{sedgewick}, \citealt{press}).  These three
-sorting algorithms all run in a time of $O(n \log n)$ in the best
-case, but have different worse-case run times.  Implementations of all
-three sorting algorithms are described in references above and in
-various places online (e.g.,
-http://www.iti.fh-flensburg.de/lang/algorithmen/sortieren/algoen.htm).
-The linux \code{qsort} function orignially mades use of the quicksort
-algorithm, but now uses \tbd{what algorithm?}.  The GSL function
-\code{gsl_heapsort} uses the heapsort algorithm.  Both of these
-implemenations require a pointer to the comparison function, which may
-result in a slower code.  
-
-For PSLib, the sorting functions shall be implemented via
-\code{qsort}.  The function \code{psSort} shall return the sorted
-result of the input array without over-writing the input array.  The
-function \code{psSortIndex} shall return an integer index to the
-sequence of the input array without overwriting the input array.
-Given the following line of code, \code{out = psSortIndex (NULL,
-&in);}, the elements of the array \code{out} are in the sequence
-\code{in.arr[out->arr[0]]} to \code{in.arr[out->arr[in.n - 1]]}. 
-
-\subsection{Smoothing: Boxcar and Gaussian}
-\label{smooth}
-
-Smoothing may occasionally be perfomed on data.  We present the
-algorithms for two typical versions: boxcar and Gaussian smoothing.
-In both smoothing techniques, given a series of data values $f_i$, the
-smoothed values $g_i$ are determined by calculating a linear
-combination based on the input data point and its nearest $2N$
-neighbors in the form:
-
-\[ g_i = \sum_{n=-N}^N c_n f_i \]
-%
-where the values of $c_n$ determine the filter type.  For boxcar
-smoothing, the values $c_n$ are constant, and must be equal to $1/(2N
-+ 1)$ to maintain the zeroth moment of the data (care must be taken at
-the ends of the data range to reduce the value of $c_n$ as fewer input
-data points may be used).  For Gaussian smoothing, the crucial
-parameter is $\sigma$, the standard deviation.  The
-value of $N$ should be chosen to be large enough, $N = 5\sigma$, and
-the values of $c_n$ are then just the Gaussian curve:
-
-\[ c_n = \frac{e^{\frac{-n^2}{2\sigma^2}}}{\sqrt{2\pi\sigma^2}} \]
-
-\subsection{Statistics}
-
-The general statistics function \code{psStats} performs a variety of
-statistical calculations on a collection of floating point numbers.
-These statistics include both sample values, in which the formal
-statistic is calculated for the complete sample, and robust values, in
-which the statistic is estimated on the basis of an assumption of an
-underlying population.  While more reliable in general, the robust
-estimators may be somewhat slower than the sample statisic, so their
-use may vary depending on the context.  In addition, certain
-sigma-clipped values are defined as an intermediate choice between the
-sample and robust estimators.  
-
-\subsubsection{Sample Statistics}
-
-We define the following statistical terms, assuming there is a set of
-data elements $x_i$.
-
-\paragraph{Mean}
-
-The mean is defined as: \[ \bar{x} = \frac{1}{N} \sum_{i = 1}^{N} x_i \]
-
-\paragraph{Median}
-
-The median is defined as the value for which 50\% of the data values
-are larger and 50\% are smaller.  For a sample, the values are sorted
-and the median is given by the value of the middle element, if the
-number of values is odd, and by the average of the two middle values
-if the number of values is even.  This median should be avoided for
-samples which are large (e.g., $N > 10000$ elements) as the basic
-robust median is quicker and more accurate.  
-
-\paragraph{Upper and Lower Quartiles}
-
-The upper and lower quartiles ($U_{\frac{1}{4}}$ and
-$L_{\frac{1}{4}}$) are first and last quarter equivalents to the
-median.  The upper quartile is the value for which 25\% of the data
-are larger and 75\% are smaller.  The lower quartile is the value for
-which 75\% of the data are larger and 25\% are smaller.  For a sample,
-the values are sorted and the quartiles are given by the values of
-elements $N/4$ and $3N/4$.  For the purpose of a sample upper and
-lower quartile, it is sufficient to provide the closest integer entry
-to these values.  The sample quartiles should be avoided for samples
-which are large (e.g., $N > 10000$ elements) as the robust quartiles
-are quicker and more accurate.
-
-\paragraph{Standard Deviation}
-
-The standard deviation of the sample is given by:
-
-\[ \sigma = \sqrt{\sum_{i = 1}^N \frac{(x_i - \bar{x})^2}{N - 1}} \]
-
-To minimize the numerical rounding error, this should be calculated
-numerically as:
-
-\[ \sigma = \sqrt{\frac{1}{N - 1} [ \sum_{i = 1}^{N} (x_i - \bar{x})^2 - \frac{1}{N} (\sum_{i = 1}^{N} (x_i - \bar{x}))^2 ]} \]
-
-\subsubsection{Clipped Statistics}
-
-The clipped statistics are used to determine the mean and $\sigma$ of
-a distribution in the presence of outliers. The clipped statistics are
-quicker than the complete robust statistical estimators and more
-reliable than the sample statistics.  The clipped statistics algorithm
-produces the clipped mean and clipped standard deviation.  The
-algorithm has 2 parameters: $N$, the number of iterations and $k$, the
-multiplying factor of $\sigma$ used to exclude outliers.  Typical
-values for both $N$ and $k$ are 3.  The algorithm is as follows:
-
-\begin{enumerate}
-\item Compute the sample median. 
-\item Compute the sample standard deviation.
-\item Use the sample median as the first estimator of the mean, $\bar{x}$.
-\item Use the sample standard deviation as the first estimator of the
-  true standard deviation, $\sigma$.
-\item Repeat the following N times:
-  \begin{enumerate}
-  \item Exclude all values $x_{i}$ for while $|x_{i} - \bar{x}| > k \sigma$.
-  \item Compute the mean and standard deviation of the sub-sample.
-  \item Use the new mean for $\bar{x}$.
-  \item Use the new standard deviations for $\sigma$.
-\end{enumerate}
-\item The last calulcated value of $\bar{x}$ is the clipped mean.
-\item The last calulcated value of $\sigma$ is the clipped standard deviation.
-\end{enumerate}
-
-\subsubsection{Robust Statistics}
-
-The robust version of the statistics provide estimators of basic
-statistical concepts which are reliable even for data samples with
-significant contamination.  A typical case is the situation in which
-the data of interest represent a primary population of interest with a
-single-valued mean and standard deviation and a secondary population
-of data with a substantially different distribution.  For example, an
-image of an uncrowded night-time field may consist of a sparse
-collection of stars and an overall background level.  The majority of
-pixels have the background value, with some variance due to noise
-sources, but many are significantly higher (contributed by stars) or
-significantly lower (dead pixels).  If we want to measure the mean of
-the background, a robust mean is necessary as the outliers will
-strongly bias the sample statistics.
-
-The robust statistics are calculated by constructing a histogram of
-the values and performing the measurements on the histogram.  The
-choice of a bin size requires some care.  If the data are integer
-valued, the natural bin size is an integer.  Otherwise, the bin should
-be a fraction of an estimate of the standard deviation.  Use the
-sample upper and lower quartiles to determine an estimate of the
-standard deviation: $\sigma_e = (U_{\frac{1}{4}} - L_{\frac{1}{4}}) /
-1.34$.  The bin size shall be set at $\sigma_e / 10$.  The remaining
-steps of the algorithm are as follows:
-
-\begin{itemize}
-\item Construct the histogram with the specified bin size.
-\item Smooth the histogram by a Gaussian with $\sigma_s = \sigma_e /
-  4$.  
-\item Find the bin with the peak value in the range $L_{\frac{1}{4}}$
-  to $U_{\frac{1}{4}}$; this is the robust mode $\mbox{mode}_r$.  
-\item Determine $dL = (U_{\frac{1}{4}} - L_{\frac{1}{4}}) / 8$.
-\item Fit a Gaussian to the bins in the range $\mbox{mode}_r - dL$ to
-  $\mbox{mode}_r + dL$.
-\item The resulting fit parameters are the robust mean,
-$\mbox{mean}_r$ and the robust standard deviation, $\sigma_r$.
-\end{itemize}
-
-To determine the robust median, construct the cumulative histogram
-from the histogram above. Select the bin which contains the 50th
-percentile value and its two neighbors.  Fit a quadratic to these
-three points.  The robust median value is the coordinate of the
-quadratic which returns the 50\% value.
-
-
-\section{Polynomials}
+\subsection{Polynomials}
 
 We will employ Chebyshev polynomials (NR \S 5.8) to approximate functions:
 \begin{equation}
-f(x) = \Sum_{i=0}^{n} c_i T_i(x)
+f(x) = \sum_{i=0}^{n} c_i T_i(x)
 \end{equation}
 These have some desirable features:
@@ -450,15 +452,11 @@
 
 The first few Chebyshev polynomials are:
-\begin{equation}
-T_0(x) = 1
-
-T_1(x) = x
-
-T_2(x) = 2x^2 - 1
-
-T_3(x) = 4x^3 - 3x
-
-T_4(x) = 8x^4 - 8x^2 + 1
-\end{equation}
+\begin{eqnarray}
+T_0(x) & = & 1 \\
+T_1(x) & = & x \\
+T_2(x) & = & 2x^2 - 1 \\
+T_3(x) & = & 4x^3 - 3x \\
+T_4(x) & = & 8x^4 - 8x^2 + 1 \\
+\end{eqnarray}
 Chebyshev polynomials follow the recurrence relation:
 \begin{equation}
@@ -468,17 +466,53 @@
 Practically, Chebyshev polynomials should be evaluated using Clenshaw's recurrence
 formula (NR \S 5.5):
-\begin{equation}
-d_j = 2xd_{j+1} - d_{j+2} + c_j
-
-f(x) = x*d_1 - d_2 + 1/2 c_0
-\end{equation}
+\begin{eqnarray}
+d_j  & = & 2xd_{j+1} - d_{j+2} + c_j \\
+f(x) & = & x*d_1 - d_2 + 1/2 c_0 \\
+\end{eqnarray}
 
 It shall be the responsibility of the user to convert the domain into the range
 $-1 < x < 1$.
 
-\subsection{Multi-dimensional polynomials}
+\subsubsection{Multi-dimensional polynomials}
 
 Multi-dimensional polynomials shall be composed of multiplications of 1D polynomials.
 
 
+\subsection{Astronomy Utilities}
+
+\section{Modules}
+
+\subsection{Image Processing Modules}
+\subsubsection{debias}
+\subsubsection{mask}
+\subsubsection{trim}
+\subsubsection{flatten}
+\subsubsection{sky/fringe subtract}
+\subsubsection{warp}
+\subsubsection{stack}
+\subsubsection{difference}
+\subsubsection{kernel convolution}
+\subsubsection{special stack}
+
+\subsection{Object Detection Modules}
+\subsubsection{find peaks}
+\subsubsection{background }
+\subsubsection{aperture photometry}
+\subsubsection{get shape}
+
+\subsection{Miscellaneous Modules}
+
+\section{Analysis Stages}
+\subsection{Phase 1}
+\subsection{Phase 2}
+\subsection{Phase 3}
+\subsection{Phase 4}
+\subsection{Cal 1}
+\subsection{Cal 2}
+\subsection{Cal 3}
+\subsection{Astrom Ref}
+\subsection{Photom Ref}
+
+\section{Architectual Components}
+
 \end{document}