redo Hermite interpolant computation to hopefully be more efficient

(in practice it doesn't seem to make much difference, though :( :(),
and to be a bit easier to read/understand/modify the code


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

1 files changed, 105 insertions, 93 deletions

diff --git a/src/GeneralizedPolynomial-Uniform/interpolate.maple b/src/GeneralizedPolynomial-Uniform/interpolate.maple
index 3f627db..58a2c97 100644
--- a/src/GeneralizedPolynomial-Uniform/interpolate.maple
+++ b/src/GeneralizedPolynomial-Uniform/interpolate.maple
@@ -86,7 +86,7 @@ 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
+# * 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
@@ -105,50 +105,45 @@ end proc;
 #		+ 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_list = A set of the interpolation coefficients (coefficients in
+# 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_coord_list = A list of lists of coordinates specifying which
-#		     derivatives are computed by the deriv_proc_list[]
-#		     procedures, for example
-#			[[x], [y], [x,y]]
-# deriv_proc_list = A list of procedures for computing finite-difference
-#		    derivatives.  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
-#		    example,
-#			[
-#			  proc(i::integer, j::integer)
-#			  - 1/2*DATA(i-1,j) + 1/2*DATA(i+1,j)
-#			  end proc
+# 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
 #			,
-#			  proc(i::integer, j::integer)
-#			  - 1/2*DATA(i,j-1) + 1/2*DATA(i,j+1)
-#			  end proc
+#			  {y}   = proc(i::integer, j::integer)
+#				  - 1/2*DATA(i,j-1) + 1/2*DATA(i,j+1)
+#				  end proc
 #			,
-#			  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_list = A list 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_list = A list of positions (each a list of numeric values)
-#		    where the interpolant is to match each possible product
-#		    of 1st derivatives"1st"
-#		    derivatives (actually all derivatives the given data DATA(),
-#		    for example
-#			[[0,0], [0,1], [1,0], [1,1]]
+#			  {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
@@ -156,71 +151,88 @@ end proc;
 #
 Hermite_polynomial_interpolant :=
 proc(
-      fn::procedure, coeff_list::list(name), coord_list::list(name),
-      deriv_coord_list::list(list(name)), fd_deriv_proc_list::list(procedure),
-      fn_posn_list::list(list(numeric)), deriv_posn_list::list(list(numeric))
+      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_eset,
-      fn_expr, deriv_eset;
-
-# set of equations {fn(posn) = DATA(posn)}
-fn_eset := map(
-		  # return equation that fn(this point) = DATA(this point)
-		  proc(posn::list(integer))
-		  fn(op(posn)) = 'DATA'(op(posn));
-		  end proc
-		,
-		  {op(fn_posn_list)}
-	      );
-
-# set of sets of equations
-#	{
-#	  { deriv1(posn1) = fd_deriv1(posn1),
-#	    deriv2(posn1) = fd_deriv2(posn1),
-#	    ... },
-#	  { deriv1(posn2) = fd_deriv1(posn2),
-#	    deriv2(posn2) = fd_deriv2(posn2),
-#	    ... },
-#	  ...
-#	}
-fn_expr := fn(op(coord_list));
+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
-       #	{deriv(posn) = fd_deriv(posn)}
-       # for this point
-       proc(posn::list(integer))
-       {
-	op(
-	 zip(
-		# return equation that
-		#	deriv(posn) = fd_deriv(posn)
-		# for this deriv and this point
-		proc(deriv_coords::list(name), fd_deriv_proc::procedure)
-		local fn_deriv_proc;
-		fn_deriv_proc := unapply( diff(fn_expr,deriv_coords),
-					  op(coord_list) );
-		fn_deriv_proc(op(posn)) = fd_deriv_proc(op(posn));
-		end proc
-	      ,
-		[op(deriv_coord_list)]
-	      ,
-		[op(fd_deriv_proc_list)]
-	    )
-	  )
-       }
+       # 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
      ,
-       {op(deriv_posn_list)}
+       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);
+
+
+#
+# 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
+		  );
 
-# set of equations {deriv-i(posn-j) = fd_deriv-i(posn-j)}
-deriv_eset := `union`(op(%));
 
-# solve equations for coeffs
-solve(fn_eset union deriv_eset, {op(coeff_list)});
+#
+# compute the coefficients in terms of the DATA() values
+#
+subs(subs_eqnset, coeff_eqns);
+eval(%);
 
-# interpolant as a polynomial in the coordinates
-return subs(%, eval(fn))(op(coord_list));
+#
+# compute the interpolant as a polynomial in the coordinates
+#
+subs(%, fn(op(coord_list)));
 end proc;
 
 ################################################################################