aboutsummaryrefslogtreecommitdiff
path: root/src/GeneralizedPolynomial-Uniform/interpolate.maple
diff options
context:
space:
mode:
Diffstat (limited to 'src/GeneralizedPolynomial-Uniform/interpolate.maple')
-rw-r--r--src/GeneralizedPolynomial-Uniform/interpolate.maple509
1 files changed, 0 insertions, 509 deletions
diff --git a/src/GeneralizedPolynomial-Uniform/interpolate.maple b/src/GeneralizedPolynomial-Uniform/interpolate.maple
deleted file mode 100644
index 7337648..0000000
--- a/src/GeneralizedPolynomial-Uniform/interpolate.maple
+++ /dev/null
@@ -1,509 +0,0 @@
-# interpolate.maple -- compute interpolation formulas/coefficients
-# $Header$
-
-#
-# <<<representation of numbers, data values, etc>>>
-# Lagrange_polynomial_interpolant - compute Lagrange polynomial interpolant
-# Hermite_polynomial_interpolant - compute Hermite polynomial interpolant
-# coeffs_as_lc_of_data - coefficients of ... (linear combination of data)
-#
-# print_coeffs__lc_of_data - print C code to compute coefficients
-# print_fetch_data - print C code to fetch input array chunk into struct data
-# print_store_coeffs - print C code to store struct coeffs "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 Lagrange 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.
-#
-Lagrange_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 computes a Hermite polynomial interpolant in any
-# number of dimensions. This is a polynomial which
-# * has values which match the given data DATA() at a specified set of
-# points, and
-# * has derivatives which match the specified finite-difference derivatives
-# of the given data DATA() at a specified set of points
-#
-# For the derivative matching, we actually match all possible products
-# of 1st derivatives, i.e. in 2-D we match dx, dy, and dxy, in 3-D we
-# match dx, dy, dz, dxy, dxz, dyz, and dxyz, etc etc.
-#
-# Arguments:
-# fn = The interpolation function. This should be a procedure in the
-# coordinates, having the coefficients as global variables. For
-# example,
-# proc(x,y)
-# + c03*y^3 + c13*x*y^3 + c23*x^2*y^3 + c33*x^3*y^3
-# + c02*y^2 + c12*x*y^2 + c22*x^2*y^2 + c32*x^3*y^2
-# + c01*y + c11*x*y + c21*x^2*y + c31*x^3*y
-# + c00 + c10*x + c20*x^2 + c30*x^3
-# end proc;
-# coeff_set = A set of the interpolation coefficients (coefficients in
-# the interpolation function), for example
-# {
-# c03, c13, c23, c33,
-# c02, c12, c22, c32,
-# c01, c11, c21, c31,
-# c00, c10, c20, c30
-# }
-# coord_list = A list of the coordinates (independent variables in the
-# interpolation function), for example [x,y].
-# deriv_set = A set of equations of the form
-# {coords} = proc
-# giving the derivatives which are to be matched, and the
-# procedures to compute their finite-difference approximations.
-# Each procedure should take N_dims integer arguments specifying
-# an evaluation point, and return a suitable linear combination
-# of the DATA() for the derivative at that point. For example
-# {
-# {x} = proc(i::integer, j::integer)
-# - 1/2*DATA(i-1,j) + 1/2*DATA(i+1,j)
-# end proc
-# ,
-# {y} = proc(i::integer, j::integer)
-# - 1/2*DATA(i,j-1) + 1/2*DATA(i,j+1)
-# end proc
-# ,
-# {x,y} = proc(i::integer, j::integer)
-# - 1/4*DATA(i-1,j+1) + 1/4*DATA(i+1,j+1)
-# + 1/4*DATA(i-1,j-1) - 1/4*DATA(i+1,j-1)
-# end proc
-# }
-# fn_posn_set = A set of positions (each a list of numeric values)
-# where the interpolant is to match the given data DATA(),
-# for example
-# {[0,0], [0,1], [1,0], [1,1]}
-# deriv_posn_set = A list of positions (each a list of numeric values)
-# where the interpolant is to match the derivatives
-# specified by deriv_set , for example
-# {[0,0], [0,1], [1,0], [1,1]}
-#
-# Results:
-# This function returns the interpolating polynomial, in the form of
-# an algebraic expression in the coordinates and the data values.
-#
-Hermite_polynomial_interpolant :=
-proc(
- fn::procedure,
- coeff_set::set(name),
- coord_list::list(name),
- deriv_set::set(set(name) = procedure),
- fn_posn_set::set(list(numeric)),
- deriv_posn_set::set(list(numeric))
- )
-local fn_eqnset, deriv_eqnset, coeff_eqns, subs_eqnset;
-
-
-#
-# compute a set of equations
-# {fn(posn) = DATA(posn)}
-# giving the function values to be matched
-#
-fn_eqnset := map(
- # return equation that fn(posn) = DATA(posn)
- proc(posn::list(integer))
- fn(op(posn)) = 'DATA'(op(posn));
- end proc
- ,
- fn_posn_set
- );
-
-
-#
-# compute a set of equations
-# { diff(fn,coords)(posn) = DERIV(coords)(posn) }
-# giving the derivative values to be matched, where DERIV(coords)
-# is a placeholder for the appropriate derivative
-#
-map(
- # return set of equations for this particular derivative
- proc(deriv_coords::set(name))
- local deriv_fn;
- fn(op(coord_list));
- diff(%, op(deriv_coords));
- deriv_fn := unapply(%, op(coord_list));
- map(
- proc(posn::list(integer))
- deriv_fn(op(posn)) = 'DERIV'(op(deriv_coords))(op(posn));
- end proc
- ,
- deriv_posn_set
- );
- end proc
- ,
- map(lhs, deriv_set)
- );
-deriv_eqnset := `union`(op(%));
-
-
-#
-# solve overall set of equations for coefficients
-# in terms of DATA() and DERIV() values
-#
-coeff_eqns := solve[linear](fn_eqnset union deriv_eqnset, coeff_set);
-if (indets(map(rhs,%)) <> {})
- then error "no unique solution for coefficients -- %1 eqns for %2 coeffs",
- nops(fn_eqnset union deriv_eqnset),
- nops(coeff_set);
-fi;
-
-
-#
-# compute a set of substitution equations
-# {'DERIV'(coords) = procedure}
-#
-subs_eqnset := map(
- proc(eqn::set(name) = procedure)
- 'DERIV'(op(lhs(eqn))) = rhs(eqn);
- end proc
- ,
- deriv_set
- );
-
-
-#
-# compute the coefficients in terms of the DATA() values
-#
-subs(subs_eqnset, coeff_eqns);
-eval(%);
-
-#
-# compute the interpolant as a polynomial in the coordinates
-#
-subs(%, 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 data positions used in the interpolant.
-#
-# 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.
-#
-coeffs_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
-# coeffs_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_coeffs__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->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_fetch_data :=
-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 store the interpolation
-# coefficients in COEFF(...) expressions, eg
-# COEFF(1,-1) = factor * coeffs->coeff_p1_m1;
-#
-# Arguments:
-# posn_list = The list of positions in the molecule.
-# coeff_name_prefix = A prefix string for the coefficient names,
-# eg "factor * coeffs->coeff_"
-# file_name = The file name to write the coefficients to. This is
-# truncated before writing.
-#
-print_store_coeffs :=
-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 evaluate a molecule, i.e.
-# to compute the molecule as a linear combination of the data values.
-#
-# Arguments:
-# posn_list = The list of positions in the molecule.
-# 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_evaluate_molecule :=
-proc(
- posn_list::list(list(numeric)),
- coeff_name_prefix::string,
- data_var_name_prefix::string,
- file_name::string
- )
-
-ftruncate(file_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 + ");
-cat(op(%));
-fprintf(file_name, " %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;