* now documents the Hermite interpolator

* add table describing derivative codes
* clarify wording for Jacobian queries
* clarify wording for spurious compiler warnings
* other minur wording tweaks


git-svn-id: http://svn.cactuscode.org/arrangements/CactusBase/LocalInterp/trunk@130 df1f8a13-aa1d-4dd4-9681-27ded5b42416

1 files changed, 130 insertions, 77 deletions

diff --git a/doc/documentation.tex b/doc/documentation.tex
index fa59fec..6845291 100644
--- a/doc/documentation.tex
+++ b/doc/documentation.tex
@@ -99,6 +99,7 @@
 \def\thorn#1{{\bf #1}}
 
 \def\Cplusplus{\hbox{C\hspace{-.05em}\raisebox{.4ex}{\tiny\bf ++}}}
+\def\cf{cf.\hbox{}}
 \def\ie{i.e.\hbox{}}
 \def\eg{e.g.\hbox{}}
 
@@ -209,9 +210,11 @@ guarantees a smooth ($C^1$) interpolant and (for smooth input data)
 interpolation error.  Unfortunately, this comes at the cost of a
 larger molecule size than the same-order Lagrange interplator, and
 a much larger interpolation error if the interpolation molecules are
-off-centered.  See
-section~\ref{localinterp/sect-Hermite}
-for details.
+off-centered.
+
+For most purposes Lagrange polynomial interpolation is better;
+only use Hermite polynomial interpolation if you need a smooth
+interpolant.
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
@@ -372,13 +375,13 @@ order~2, cubic is order~3, etc.
 
 If no grid boundaries or excised points are nearby, the interpolator
 centers the molecules around the interpolation point as much as possible.
-Figure~\ref{localinterp/fig-molecule-size+centering} gives the molecule size
-and details of the centering policy for each interpolation operator
+Table~\ref{localinterp/tab-molecule-size+centering} gives the molecule
+size and details of the centering policy for each interpolation operator
 and order.
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
-\begin{figure}
+\begin{table}
 
 \def\Linterval{\lower0.60ex\hbox{\hskip-0.20em$\big[$}}
 \def\Rinterval{\lower0.60ex\hbox{\hskip-0.25em$\big)$}}
@@ -500,8 +503,8 @@ Interpolation Operator				& Order	& Size
 	use a molecule containing the grid points $\{ -2, -1, 0, 1, 2 \}$
 	for all interpolation points in the range $[-0.5, 0.5)$.
 	}
-\label{localinterp/fig-molecule-size+centering}
-\end{figure}
+\label{localinterp/tab-molecule-size+centering}
+\end{table}
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
@@ -517,7 +520,7 @@ off-center the interpolation molecules, or refuse to interpolate
 	 {\tt CCTK\_INTERP\_ERROR\_POINT\_OUTSIDE}.
 	 }%%%
 {} or \verb|CCTK_INTERP_ERROR_POINT_EXCISED| error code).  This behavior
-is controlled by the parameter-table entries
+is controlled by the (optional) parameter-table entries
 \begin{verbatim}
 const CCTK_INT N_boundary_points_to_omit[2*N_dims]
 const CCTK_REAL boundary_off_centering_tolerance[2*N_dims]
@@ -560,7 +563,7 @@ zones.
 \subsubsection{Off-Centering Molecules near Grid Boundaries
 	       and/or Excision Regions}
 
-The parameter-table entries
+The (optional) parameter-table entries
 \begin{verbatim}
 const CCTK_REAL boundary_off_centering_tolerance[2*N_dims]
 const CCTK_REAL boundary_extrapolation_tolerance[2*N_dims]
@@ -581,7 +584,7 @@ regions in the data coordinate space:
 \item[The ``default-centering region'']
 	is that region in interpolation-point space where the
 	interpolator can use the default molecule centering
-	(described in figure~\ref{localinterp/fig-molecule-size+centering}),
+	(described in table~\ref{localinterp/tab-molecule-size+centering}),
 	\ie{} where the default-centering molecules require data
 	only from the data-valid region.
 \end{description}
@@ -828,19 +831,21 @@ molecules.  Hypercube-shaped molecules are the simplest and most common
 case, but one could also imagine (say) octagon-shaped molecules in 2-D,
 or some generalization of this in higher numbers of dimensions.
 
-The parameter
+The (optional) parameter
 \begin{verbatim}
 /* set with Util_TableSetString() */
 const char molecule_family[];
 \end{verbatim}
 may be used to set or query the molecule family.
 
-If this key is present in the parameter table, this sets the
-molecule family/shape based on that string.
-
-If this key isn't present in the parameter table, then the
+If this key is present in the parameter table, the interpolator sets
+the molecule family/shape based on the value specified.
+If this key {\em isn't\/} present in the parameter table, then the
 interpolator sets it to the molecule family being used.
 
+At present only hypercubed-shaped molecules are implemented; these
+are specified by \verb|molecule_family| set to \verb|"cube"|.
+
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
 \subsection{Non-Contiguous Input Arrays}
@@ -898,8 +903,7 @@ then instead of estimating values of that function at the interpolation
 points, the interpolator can instead or additionally estimate values
 of various derivatives of that function at the interpolation points.
 (We don't currently implement it, but one could also imagine
-interpolating more general molecule operations such as Laplacians,
-running integrals, etc.)
+interpolating more general molecule operations such as Laplacians.)
 
 The following (optional) parameter-table entries are used to specify
 this:
@@ -928,21 +932,43 @@ use.  An \verb|operation_codes[out]| value which is $\ge 0$ is taken
 as a decimal integer encoding a coordinate partial derivative: each
 decimal digit means to take the coordinate partial derivative along
 that (1-origin) axis; the order of the digits in a number is ignored.
-For example:
-\begin{flushleft}
+Table~\ref{localinterp/tab-derivative-codes} summarizes these resulting
+derivative codes.
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+\begin{table}
+\begin{center}
 \begin{tabular}{cl}
 \verb|operation_codes[out]|
-	& What it means								\\
-\hline %-----------------------------------------------------------------------
-0	& no derivative, ``just'' interpolate the input array			\\
-1	& interpolate $\partial \big/ \partial x^1$ of the input array		\\
-2	& interpolate $\partial \big/ \partial x^2$ of the input array		\\
-12 or 21
-	& interpolate $\partial^2 \big/ \partial x^1 \partial x^2$
-	  of the input array							\\
-33	& interpolate $\partial^2 \big/ \partial (x^3)^2$ of the input array	%%%\\
+	& What it means							\\
+\hline %---------------------------------------------------------------
+0	& interpolate the input array itself (no derivative)		\\[1ex]
+1	& interpolate $\partial \big/ \partial x^1$ of the input array	\\
+2	& interpolate $\partial \big/ \partial x^2$ of the input array	\\
+3	& interpolate $\partial \big/ \partial x^3$ of the input array	\\[1ex]
+12 or 21& interpolate $\partial^2 \big/ \partial x^1 \partial x^2$
+	  of the input array						\\
+13 or 31& interpolate $\partial^2 \big/ \partial x^1 \partial x^3$
+	  of the input array						\\
+23 or 32& interpolate $\partial^2 \big/ \partial x^2 \partial x^3$
+	  of the input array						\\
+11	& interpolate $\partial^2 \big/ \partial (x^1)^2$ of the input array \\
+22	& interpolate $\partial^2 \big/ \partial (x^2)^2$ of the input array \\
+33	& interpolate $\partial^2 \big/ \partial (x^3)^2$ of the input array \\
+\hline %---------------------------------------------------------------
 \end{tabular}
-\end{flushleft}
+\end{center}
+\caption[Derivative Codes]
+	{
+	This table gives the codes in {\tt operation\_codes[out]}
+	for each possible 1st or 2nd~derivative in 3-D; for 1-D or
+	2-D codes referring to nonexistent coordinates are invalid.
+	}
+\label{localinterp/tab-derivative-codes}
+\end{table}
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
 At present we do {\em not\/} have any \verb|#define|s for the
 operation codes in the Cactus header files.  However, you can avoid
@@ -960,16 +986,17 @@ section~\ref{localinterp/sect-example-derivatives}.
 \subsection{Jacobian and Domain of Dependence}
 \label{localinterp/sect-Jacobian}
 
-Sometimes we want to know the Jacobian of the interpolator,
+The Jacobian of the interpolator is defined as
 \begin{equation}
 \frac{\partial \hbox{\tt output\_array[out][pt]}}
      {\partial \hbox{\tt input\_array[in][(i,j,k)]}}
 						\label{localinterp/eqn-Jacobian}
 \end{equation}
-and/or ``just'' the set of \verb|(i,j,k)| for which this is nonzero
-(\ie{} the sparsity structure of the Jacobian, equivalently the domain
-of dependence of the output result, or equivalently the interpolation
-molecule size and shape) for a given \verb|out|, \verb|in|, and \verb|pt|.%%%
+We may want to know the Jacobian itself, and/or ``just'' the set of
+\verb|(i,j,k)| for which this is nonzero (\ie{} the Jacobian's sparsity
+structure, or equivalently the domain of dependence of the output result,
+or equivalently the interpolation molecule size and shape) for a given
+\verb|out|, \verb|in|, and \verb|pt|.%%%
 \footnote{%%%
 	 For something like a spline interpolator the Jacobian
 	 is generally nonzero for all {\tt (i,j,k)}, but
@@ -1006,9 +1033,14 @@ questions:
 	on the actual floating-point values being interpolated?
 	Equivalently, is the interpolation nonlinear?
 \end{itemize}
-Because the different cases differ so much in their complexity,
-we define several distinct APIs for obtaining the interpolator's
-Jacobian and/or domain of dependence.
+
+Because the different cases differ so much in the amount of information
+required to describe the Jacobian, it's hard to define a single API
+which covers all cases without burdening the simpler cases with
+excessive complexity.  Instead, the interpolator defines a
+Jacobian-structure--query API to determine which case applies,
+together with (conceptually) several different APIs for the
+different cases.  (At the moment only a single case is implemented.)
 
 %%%%%%%%%%%%%%%%%%%%
 
@@ -1031,19 +1063,18 @@ to 1.  Otherwise (\ie{} if the interpolation molecule size and shape
 are independent of the interpolation coordinates) it should set this
 parameter to 0.
 
-If the interpolator supports computing derivatives as described
-in section~\ref{localinterp/sect-derivatives}, {\em and\/}
+If the interpolator supports computing derivatives as described in
+section~\ref{localinterp/sect-derivatives}, {\em and\/}
 if the interpolation molecule's size and/or shape varies with
-\verb|operation_codes[]|, the interpolator should sets the
-parameter
+\verb|operation_codes[]|, the interpolator sets the parameter
 \begin{verbatim}
 CCTK_INT MSS_is_fn_of_which_operation;
 \end{verbatim}
 to 1.  Otherwise (\ie{} if the interpolator doesn't support computing
 derivatives, or if the interpolator does support computing derivatives
 but the interpolation molecule size and shape are independent of the
-\verb|operation_code[]| values), it should set this parameter to 0.
-Note that this query tests whether the molecule size and/or shape
+\verb|operation_code[]| values), the interpolator sets this parameter
+to 0.  Note that this query tests whether the molecule size and/or shape
 depend on \verb|operation_codes[]| in general, independent of whether
 there are in fact any distinct values (or even any values at all) passed
 in \verb|operation_codes[]| in this particular interpolator call.  In
@@ -1052,24 +1083,33 @@ not about this particular call.
 
 If the interpolation molecule's size and/or shape varies with the
 actual floating-point values of the input arrays, the interpolator
-should set the parameter
+sets the parameter
 \begin{verbatim}
 CCTK_INT MSS_is_fn_of_input_array_values;
 \end{verbatim}
 to 1.  Otherwise (\ie{} if the interpolation molecule size and shape
-are independent of the input array values; this is a necessary, but not
-sufficient, condition for the interpolation to be linear), it should
-set this parameter to 0.
+are independent of the input array values; this is a necessary, but
+not sufficient, condition for the interpolation to be linear), the
+interpolator sets this parameter to 0.
 
-If the actual floating-point values of the Jacobian~$(\ref{localinterp/eqn-Jacobian})$
+If the actual floating-point values of the
+Jacobian~$(\ref{localinterp/eqn-Jacobian})$
 (for a given \verb|out|, \verb|in|, and \verb|pt|) depend on the actual
 floating-point values of the input arrays (\ie{} if the interpolation
-is nonlinear), the interpolator should set the parameter
+is nonlinear), the interpolator sets the parameter
 \begin{verbatim}
 CCTK_INT Jacobian_is_fn_of_input_array_values;
 \end{verbatim}
-to 1.  Otherwise (\ie{} if the interpolation is linear) it should set
-this parameter to 0.
+to 1.  Otherwise (\ie{} if the interpolation is linear) the interpolator
+sets this parameter to 0.
+
+The current implementation always sets
+\begin{verbatim}
+MSS_is_fn_of_interp_coords = 0
+MSS_is_fn_of_which_operation = 0
+MSS_is_fn_of_input_array_values = 0
+Jacobian_is_fn_of_input_array_values = 0
+\end{verbatim}
 
 %%%%%%%%%%%%%%%%%%%%
 
@@ -1091,6 +1131,8 @@ MSS_is_fn_of_which_operation = 0
 MSS_is_fn_of_input_array_values = 0
 Jacobian_is_fn_of_input_array_values = 0
 \end{verbatim}
+(These are precisely the values set by the current implementation.)
+In the rest of this section we describe the query API for this case.
 
 The following parameters may be used to query the molecule size:
 \begin{verbatim}
@@ -1157,11 +1199,10 @@ An example may help to clarify this:  Suppose we have a 1-D grid
 with 11~grid points, with integer subscripts~0 through~10 inclusive,
 and interpolation coordinates given by \verb|coord_origin = 0.0|
 and \verb|coord_delta = 0.1|.
-Suppose further that we're doing quadratic interpolation, using
-3-point (vertex-centered) molecules, and that the interpolator
-centers the interpolation molecules as close to the interpolation
-point as possible subject to the constraint that the molecules
-never require data from outside the grid.
+Suppose further that we're doing Lagrange polynomial interpolation,
+with \verb|order = 2| and hence
+(by table~\ref{localinterp/tab-molecule-size+centering})
+using 3-point molecules.
 Finally, suppose that we query the Jacobian molecule positions for
 the \verb|N_interp_points=14| interpolation points 0.0, 0.04, 0.06,
 0.10, 0.14, 0.16, 0.20, 0.80, 0.84, 0.86, 0.90, 0.94, 0.96, 1.00.
@@ -1194,14 +1235,16 @@ interp_molecule_positions[0][13] =  9    /*   1.00     [0.8, 0.9, 1.0] */
 The way the generalized polynomial interpolator is implemented it's
 easy to also do Savitzky-Golay smoothing.%%%
 \footnote{%%%
-	 See section 14.8 of Numerical Recipes
-	 2nd edition for a general discussion of
+	 See section~14.8 of Numerical Recipes
+	 (2nd edition) for a general discussion of
 	 Savitzky-Golay smoothing.
 	 }%%%
 {}  This is best described by way of an example:  Suppose we're doing
-1-D cubic interpolation, \ie{} at each interpolation point we're using
-a cubic polynomial fitted to 4~surrounding data points.  For
-Savitzky-Golay smoothing, we would instead {\em least-squares fit\/}
+1-D cubic interpolation, using
+(by table~\ref{localinterp/tab-molecule-size+centering})
+4-point molecules.  In other words, at each interpolation point we
+use a cubic interpolation polynomial fitted to 4~surrounding data points.
+For Savitzky-Golay smoothing, we would instead {\em least-squares fit\/}
 a cubic polynomial to some {\em larger\/} number of surrounding data
 points.  This combines interpolation with smoothing, so there's less
 amplification of noise in the input data in the interpolation outputs.
@@ -1212,10 +1255,10 @@ const CCTK_INT smoothing;
 \end{verbatim}
 specifies how much (how many points) to enlarge the interpolation
 molecule for this.  The default is 0 (no smoothing).  1 would mean to
-enlarge the molecule by 1 point, \eg{}  to use a 5-point molecule
-instead of the usual 4-point one for cubic interpolation.  2 would mean
-to enlarge by 2 points, \eg{}  to use a 6-point molecule for cubic
-interpolation.  Etc etc.
+enlarge the molecule by 1 point (\eg{} to use a 5-point molecule instead
+of the usual 4-point one for cubic interpolation).  2 would mean to
+enlarge by 2 points, (\eg{} to use a 6-point molecule for cubic
+interpolation).  Etc etc.
 
 Note that in $>1$~dimension, the default hypercube-shaped molecules
 already use more data points than the number of free parameters in
@@ -1229,8 +1272,7 @@ Savitzky-Golay smoothing is basically free apart from the increase in
 the molecule size, e.g. a \verb|smoothing|=2 cubic interpolation has
 exactly the same cost as any other 6-point--molecule interpolation.
 
-Alas, at the moment only the trivial case \verb|smoothing|=0 is
-implemented, but the framework is all there for more general cases.
+The current implementation but only supports \verb|smoothing=0|.
 
 %%%%%%%%%%%%%%%%%%%%
 
@@ -1257,15 +1299,26 @@ though, \eg{} the current coefficients for 3-D for orders~1-4
 take about 8~cpu minutes to generate using Maple~7 on a 1.7~GHz~P4.
 
 Note that when compiling the code in the directory
-\verb|LocalInterp/src/UniformCartesian/|, you may get compiler
-warnings about casts discarding \verb|const| qualifiers from pointers
-in (the \verb|#include|-ed file) \verb|template.c|.  Don't worry --
-the code is actually ok, the problem is that C's \verb|const|-qualifier
-rules are overly strict in some cases.  See question 11.10 in the
-(excellent!) \verb|comp.lang.c| online FAQ list
-(\verb|http://www.eskimo.com/~scs/C-faq/top.html|) for details.
-(I think \Cplusplus{} has at least partially fixed the
-\verb|const|-qualifier rules.)
+\verb|src/GeneralizedPolynomial-Uniform| of this thorn, you may get
+compiler warnings about casts discarding \verb|const| qualifiers from
+pointers in (the \verb|#include|-ed file) \verb|template.c|.  Don't
+worry -- the code is actually ok.%%%
+\footnote{%%%
+	 The warnings are really a design bug in C's
+	 const-qualifier rules when pointers point to
+	 other pointers.  See question~11.10 in the
+	 (excellent!) {\tt comp.lang.c} online FAQ
+	 ({\tt http://www.eskimo.com/~scs/C-faq/top.html})
+	 for further discussion.  gcc~2.95.3 gives the
+	 spurious warnings, but gcc~3.2.1 doesn't, nor
+	 does the Intel C compiler in any version I've
+	 tried.  Also, \Cplusplus{} has fixed the problems
+	 in C's const-qualifier rules, so the same code
+	 compiled as \Cplusplus{} shouldn't give any
+	 warnings (and doesn't on the compilers I've tried).
+	 }%%%
+{}  You may also get compiler warnings about unused variables in
+this same file; again don't worry, the code is ok.
 
 See the \verb|README| file in the source code directory
 \verb|LocalInterp/src/UniformCartesian/| for further details on the