initial checkin of new LocalInterp thorn

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

1 files changed, 313 insertions, 0 deletions

diff --git a/src/GeneralizedPolynomial-Uniform/interpolate.maple b/src/GeneralizedPolynomial-Uniform/interpolate.maple
new file mode 100644
index 0000000..4938f86
--- /dev/null
+++ b/src/GeneralizedPolynomial-Uniform/interpolate.maple
@@ -0,0 +1,313 @@
+# GInterpolate.maple -- compute generalized interpolation formulas/coefficients
+# $Id$
+
+#
+# <<<representation of numbers, data values, etc>>>
+# polynomial_interpolant - compute polynomial interpolant
+# coeff_as_lc_of_data - coefficients of ... (linear combination of data)
+#
+# print_coeff__lc_of_data - print C code to compute coefficients
+# print_data_var_assign - print C code to assign data-value variables
+# print_interp_cmpt__lc_of_data - print C code for computation of interpolant
+#
+# coeff_name - name of coefficient of data at a given [m] coordinate
+# data_var_name - name of variable storing data value at a given [m] coordinate
+#
+
+################################################################################
+
+#
+# ***** representation of numbers, data values, etc *****
+#
+# We use RATIONAL(p.0,q.0) to denote the rational number p/q.
+#
+# We use DATA(...) to represent the data values being interpolated at a
+# specified [m] coordinate, where the arguments are the [m] coordinates.
+#
+# We use COEFF(...) to represent the molecule coefficient at a specified
+# [m] coordinate, where the arguments are the [m] coordinates.
+#
+# For example, the usual 1-D centered 2nd order 1st derivative molecule
+# would be written
+#	RATIONAL(-1.0,2.0)*DATA(-1) + RATIONA(1.0,2.0)*DATA(1)
+# and its coefficients as
+#	COEFF(-1) = RATIONAL(-1.0,2.0)
+#	COEFF(1) = RATIONAL(1.0,2.0)
+#
+
+################################################################################
+################################################################################
+################################################################################
+
+#
+# This function computes a polynomial interpolant in any number of dimensions.
+#
+# Arguments:
+# fn = The interpolation function.  This should be a procedure in the
+#      coordinates, having the coefficients as global variables.  For
+#      example,
+#	  proc(x,y) c00 + c10*x + c01*y end proc
+# coeff_list = A set of the interpolation coefficients (coefficients in
+#	       the interpolation function), for example [c00, c10, c01].
+# coord_list = A list of the coordinates (independent variables in the
+#	       interpolation function), for example [x,y].
+# posn_list = A list of positions (each a list of numeric values) where the
+#	      interpolant is to use data, for example  hypercube([0,0], [1,1]).
+#	      Any positions may be used; if they're redundant (as in the
+#	      example) the least-squares interpolant is computed.
+#
+# Results:
+# This function returns the interpolating polynomial, in the form of
+# an algebraic expression in the coordinates and the data values.
+#
+polynomial_interpolant :=
+proc(
+      fn::procedure, coeff_list::list(name),
+      coord_list::list(name), posn_list::list(list(numeric))
+    )
+local posn, data_eqns, coeff_eqns;
+
+# coefficients of interpolating polynomial
+data_eqns := {  seq( fn(op(posn))='DATA'(op(posn)) , posn=posn_list )  };
+coeff_eqns := linalg[leastsqrs](data_eqns, {op(coeff_list)});
+if (has(coeff_eqns, '_t'))
+   then error "interpolation coefficients aren't uniquely determined!";
+end if;
+
+# interpolant as a polynomial in the coordinates
+return subs(coeff_eqns, eval(fn))(op(coord_list));
+end proc;
+
+################################################################################
+
+#
+# This function takes as input an interpolating polynomial, expresses
+# it as a linear combination of the data values, and returns the coefficeints
+# of that form.
+# 
+# Arguments:
+# interpolant = The interpolating polynomial (an algebraic expression
+#		in the coordinates and the data values).
+# posn_list = The same list of positions as was used to compute the
+#	      interpolating polynomial.
+#
+# Results:
+# This function returns the coefficients, as a list of equations of the
+# form   COEFF(...) = value , where each  value  is a polynomial in the
+# coordinates.  The order of the list matches that of  posn_list.
+#
+coeff_as_lc_of_data :=
+proc(
+      interpolant::algebraic,
+      posn_list::list(list(numeric))
+    )
+local data_list, interpolant_as_lc_of_data;
+
+# interpolant as a linear combination of the data values
+data_list := [ seq( 'DATA'(op(posn)) , posn=posn_list ) ];
+interpolant_as_lc_of_data := collect(interpolant, data_list);
+
+# coefficients of the data values in the linear combination
+return map(
+	      proc(posn::list(numeric))
+	      coeff(interpolant_as_lc_of_data, DATA(op(posn)));
+	      'COEFF'(op(posn)) = %;
+	      end proc
+	    ,
+	      posn_list
+	  );
+end proc;
+
+################################################################################
+################################################################################
+################################################################################
+
+#
+# This function prints C expressions for the coefficients of an
+# interpolating polynomial.  (The polynomial is expressed as linear
+# combinations of the data values with coefficients which are
+# RATIONAL(p,q) calls.)
+#
+# Arguments:
+# coeff_list = A list of the coefficients, as returned from
+#	       coeff_as_lc_of_data() .
+# coeff_name_prefix = A prefix string for the coefficient names.
+# temp_name_type = The C type to be used for Maple-introduced temporary
+#		   names, eg. "double".
+# file_name = The file name to write the coefficients to.  This is
+#	      truncated before writing.
+#
+print_coeff__lc_of_data :=
+proc( coeff_list::list(specfunc(numeric,COEFF) = algebraic),
+      coeff_name_prefix::string,
+      temp_name_type::string,
+      file_name::string )
+global `codegen/C/function/informed`;
+local coeff_list2, cmpt_list, temp_name_list;
+
+# convert LHS of each equation from a COEFF() call (eg COEFF(-1,+1))
+# to a Maple/C variable name (eg coeff_I_m1_p1)
+coeff_list2 := map(
+		      proc(coeff_eqn::specfunc(numeric,COEFF) = algebraic)
+		      local posn;
+		      posn := [op(lhs(coeff_eqn))];
+		      coeff_name(posn,coeff_name_prefix);
+		      convert(%, name);	# codegen[C] wants LHS
+					# to be an actual Maple *name*
+		      % = fix_rationals(rhs(coeff_eqn));
+		      end proc
+		    ,
+		      coeff_list
+		  );
+
+#
+# generate the C code
+#
+
+# tell codegen[C] not to warn about unknown RATIONAL() and DATA() "fn calls"
+# via undocumented :( global table
+`codegen/C/function/informed`['RATIONAL'] := true;
+`codegen/C/function/informed`['DATA'] := true;
+
+ftruncate(file_name);
+
+# optimized computation sequence for all the coefficients
+# (may use local variables t0,t1,t2,...)
+cmpt_list := [codegen[optimize](coeff_list2, tryhard)];
+
+# list of the t0,t1,t2,... local variables
+temp_name_list := nonmatching_names(map(lhs,cmpt_list), coeff_name_prefix);
+
+# declare the t0,t1,t2,... local variables (if there are any)
+if (nops(temp_name_list) > 0)
+   then print_name_list_dcl(%, temp_name_type, file_name);
+fi;
+
+# now print the optimized computation sequence
+codegen[C](cmpt_list, filename=file_name);
+
+NULL;
+end proc;
+
+################################################################################
+
+#
+# This function prints a sequence of C expression to assign the data-value
+# variables.
+#
+# Arguments:
+# posn_list = The same list of positions as was used to compute the
+#	      interpolating polynomial.
+# data_var_name_prefix = A prefix string for the data variable names.
+# file_name = The file name to write the coefficients to.  This is
+#	      truncated before writing.
+#
+print_data_var_assign :=
+proc(
+      posn_list::list(list(numeric)),
+      data_var_name_prefix::string,
+      file_name::string
+    )
+
+ftruncate(file_name);
+
+map(
+       proc(posn::list(numeric))
+       fprintf(file_name,
+	       "%s = %a;\n",
+	       data_var_name(posn,data_var_name_prefix),
+	       DATA(op(posn)));
+       end proc
+     ,
+       posn_list
+   );
+
+NULL;
+end proc;
+
+################################################################################
+
+#
+# This function prints a C expression to compute the interpolant,
+# using the coefficients computed by  print_coeff__lc_of_data()
+# (i.e. expressing the interpolant as a linear combination of the
+# data values).
+#
+# Arguments:
+# posn_list = The same list of positions as was used to compute the
+#	      interpolating polynomial.
+# result_var_name = The (string) name of the variable to which the
+#		    result is to be assigned.
+# coeff_name_prefix = A prefix string for the coefficient names.
+# data_var_name_prefix = A prefix string for the data variable names.
+# file_name = The file name to write the coefficients to.  This is
+#	      truncated before writing.
+#
+print_interp_cmpt__lc_of_data :=
+proc(
+      posn_list::list(list(numeric)),
+      result_var_name::string,
+      coeff_name_prefix::string,
+      data_var_name_prefix::string,
+      file_name::string
+    )
+
+ftruncate(file_name);
+fprintf(file_name, "%s =\n", result_var_name);
+
+# list of "coeff*data_var" terms
+map(
+       proc(posn::list(numeric))
+       sprintf("%s*%s",
+	       coeff_name(posn,coeff_name_prefix),
+	       data_var_name(posn,data_var_name_prefix));
+       end proc
+     ,
+       posn_list
+   );
+
+ListTools[Join](%, "\n\t+ ");
+cat(op(%));
+
+fprintf(file_name, "\t%s;\n", %);
+NULL;
+end proc;
+
+################################################################################
+################################################################################
+################################################################################
+
+#
+# This function computes the name of the coefficient of the data at a
+# given [m] position, i.e. it encapsulates our naming convention for this.
+#
+# Arguments:
+# posn = (in) The [m] coordinates.
+# name_prefix = A prefix string for the coefficient name.
+#
+# Results:
+# The function returns the coefficient, as a Maple string.
+#
+coeff_name :=
+proc(posn::list(numeric), name_prefix::string)
+cat(name_prefix, sprint_numeric_list(posn));
+end proc;
+
+################################################################################
+
+#
+# This function computes the name of the variable in which the C code
+# will store the input data at a given [m] position, i.e. it encapsulates
+# our naming convention for this.
+#
+# Arguments:
+# posn = (in) The [m] coordinates.
+# name_prefix = A prefix string for the variable name.
+#
+# Results:
+# The function returns the variable name, as a Maple string.
+#
+data_var_name :=
+proc(posn::list(numeric), name_prefix::string)
+cat(name_prefix, sprint_numeric_list(posn));
+end proc;