# interpolate.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_coeff_var_store - print C code to store coeff vars "somewhere"
# 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);

fclose(file_name);

NULL;
end proc;

################################################################################

#
# This function prints a sequence of C expression to assign the data-value
# variables, eg
#	data_m1_p1 = DATA(-1,1);
#
# 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
   );
fclose(file_name);

NULL;
end proc;

################################################################################

#
# This function prints a sequence of C expression to assign the interpolation
# coefficients to  COEFF(...)  expressions, eg
#	COEFF(1,-1) = coeff_dx_p1_m1
#
# Arguments:
# posn_list = The same list of positions as was used to compute the
#	      interpolating polynomial.
# coeff_name_prefix = A prefix string for the coefficient names.
# file_name = The file name to write the coefficients to.  This is
#	      truncated before writing.
#
print_interp_coeff_var_store :=
proc(
      posn_list::list(list(numeric)),
      coeff_name_prefix::string,
      file_name::string
    )

ftruncate(file_name);
map(
       proc(posn::list(numeric))
       fprintf(file_name,
	       "%a = %s;\n",
	       'COEFF'(op(posn)),
	       coeff_name(posn,coeff_name_prefix));
       end proc
     ,
       posn_list
   );
fclose(file_name);

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", %);

fclose(file_name);

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;