design changes to support Cactus-style optimizations,

start implementation of interpolate()


git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinAnalysis/AHFinderDirect/trunk@117 f88db872-0e4f-0410-b76b-b9085cfa78c5

2 files changed, 136 insertions, 38 deletions

diff --git a/src/jtutil/interpolate.cc b/src/jtutil/interpolate.cc
index b1aef2f..661f81d 100644
--- a/src/jtutil/interpolate.cc
+++ b/src/jtutil/interpolate.cc
@@ -4,6 +4,7 @@
 // interpolation::Lagrange_uniform_1d::Lagrange_uniform
 // interpolation::Lagrange_uniform_1d::setup_data_x
 // interpolation::Lagrange_uniform_1d::setup_data_y
+// interpolation::Lagrange_uniform_1d::setup_interp_x
 // interpolation::Lagrange_uniform_1d::bdod
 // interpolation::Lagrange_uniform_1d::interpolate
 // interpolation::Lagrange_uniform_1d::Jacobian
@@ -30,8 +31,26 @@ namespace interpolation
 template <typename fpt>
   Lagrange_uniform_1d<fpt>::Lagrange_uniform_1d(int order_in)
 	: order(order_in), N_interp(order_in+1),
+	  x_interp_valid(false), x_interp(0.0), end_action('?'),
 	  x_map_ptr(NULL), y_array(NULL)
 {
+switch	(order)
+	{
+case 1:
+case 2:
+case 3:
+case 4:
+case 5:
+case 6:
+	// yes, we support this ==> no-op here
+	break;
+default:
+	error_exit(ERROR_EXIT,
+"***** Lagrange_uniform_1d<fpt>::Lagrange_uniform_1d():\n"
+"        order=%d isn't supported!\n",
+,
+		   order);					/*NOTREACHED*/
+	}	
 }
 	  }
 
@@ -47,6 +66,7 @@ template <typename fpt>
   Lagrange_uniform_1d<fpt>::setup_data_x(const linear_map<fpt>& x_map_in)
 {
 x_map_ptr = &x_map_in;
+
 if (x_map_ptr->N_points() < N_interp)
    then error_exit(ERROR_EXIT,
 "***** Lagrange_uniform_1d<fpt>::setup_data_x():\n"
@@ -83,18 +103,58 @@ y_array = y_array_in;
 //*****************************************************************************
 
 //
-// This function does the actual interpolation.
+// This function is used by client code to set up
+// * the x coordinate at which interpolation is to be done
+// * what sort of off-the-end-of-the-data handling is desired
 //
 namespace interpolation
 	  {
 template <typename fpt>
-  fpt Lagrange_uniform_1d<fpt>::interpolate(fpt x, char end_action = 'o')
-	const
+  void Lagrange_uniform_1d<fpt>::setup_interp_x(fpt x_interp_in,
+						char end_action_in = 'o')
+{
+x_interp = x_interp_in,
+end_action = end_action_in;
+x_interp_valid = true;
+}
+	  }
+
+//*****************************************************************************
+
+//
+// This function is the user-level driver to do the actual interpolation.
+//
+namespace interpolation
+	  {
+template <typename fpt>
+  fpt Lagrange_uniform_1d<fpt>::interpolate() const
 {
 assert(x_map_ptr != NULL);
 assert(y_array != NULL);
+assert(x_interp_valid);
 
-// FIXME FIXME FIXME
+int i_lo, i_hi;
+const int i_center = bdod(i_lo, i_hi);
+
+const fpt x_relative = x - x_map_ptr->fp_of_int(i_center);
+const fpt x_scaled = x_map_ptr->inverse_delta_fp() * x_relative;
+
+switch	(order)
+	{
+case 1:	return interpolate_local_order1(x_scaled, & y_array[i_center]);
+case 2:	return interpolate_local_order2(x_scaled, & y_array[i_center]);
+case 3:	return interpolate_local_order3(x_scaled, & y_array[i_center]);
+case 4:	return interpolate_local_order4(x_scaled, & y_array[i_center]);
+case 5:	return interpolate_local_order5(x_scaled, & y_array[i_center]);
+case 6:	return interpolate_local_order6(x_scaled, & y_array[i_center]);
+default:
+	error_exit(PANIC_EXIT,
+"***** Lagrange_uniform_1d<fpt>::intgerpolate():\n"
+"        order=%d isn't supported!\n",
+"        (this should never happen!)\n"
+,
+		   order);					/*NOTREACHED*/
+	}
 }
 	  }
 
@@ -110,9 +170,7 @@ assert(y_array != NULL);
 namespace interpolation
 	  {
 template <typename fpt>
-  void Lagrange_uniform_1d<fpt>::bdod(fpt x, char end_action = 'o',
-				      int& i_lo, int& i_hi)
-	const
+  int Lagrange_uniform_1d<fpt>::bdod(int& i_lo, int& i_hi) const
 {
 //
 // Suppose we have a data array indexed by [i], from which we wish to
@@ -140,6 +198,7 @@ template <typename fpt>
 //
 
 assert(x_map_ptr != NULL);
+assert(x_valid_flag);
 
 const enum linear_map<fpt>::noninteger_action floor_or_round
 	= jtutil::is_even(N_interp) ? linear_map<fpt>::floor
@@ -172,6 +231,7 @@ case 'e':
 			   double(x_map_ptr->max_fp()),
 			   i_center, i_lo, i_hi);		/*NOTREACHED*/
 case 'o':
+case '*':
 	if (i_lo < x_map_ptr->min_int())
 	   then {
 		const int shift = x_map_ptr->min_int() - i_lo;
@@ -189,11 +249,11 @@ case 'o':
 		break;
 		}
 	break;
-	
-case '*':
+
 default:
 	error_exit(ERROR_EXIT,
 "***** Lagrange_uniform<fpt>::bdod(): invalid end_action='%c'='0x%02x'!\n",
+"                                     (this should never happen!)\n"
 		   int(end_action), unsigned int(end_action));	/*NOTREACHED*/
 	}
 
diff --git a/src/jtutil/interpolate.hh b/src/jtutil/interpolate.hh
index c636fa5..715ac65 100644
--- a/src/jtutil/interpolate.hh
+++ b/src/jtutil/interpolate.hh
@@ -68,7 +68,7 @@
 // interpolation position lies (fuzzily) outside the data range.
 // '*' ==> Like 'o', but also allow the interpolation position to lie
 //	   (fuzzily) outside the data range.  In other words, this case
-//	   case accepts *any* interpolation position.
+//	   accepts *any* interpolation position.
 //
 // For Lagrange interpolation all three possibilities are valid;
 // for spline interpolation the "off-centering as needed" is intrinsic
@@ -89,15 +89,11 @@
 // 4. Set up the x coordinates where interpolated data is desired.
 // 5. Do the actual interpolation.
 //
-// The routines here all allow 2-5, 3-5, and/or 4-5 to be repeated as needed.
 // The idea is that earlier steps can preallocate storage arrays and/or
 // precompute coefficients to be used by later stages, so repeating the
 // later stages may have lower overheads than redoing the entire process.
-//
-// One could also imagine doing these steps in the order 12435, i.e.
-// interpolating many different sets of date to the same position(s).
-// The Cactus interpolator interface supports this nicely, but alas
-// the interface here doesn't support this at all. :( :(
+// The interface here requires that 1 be done first, then allows 2-4
+// to be done in any order before 5.
 //
 
 //
@@ -109,15 +105,17 @@
 
 //******************************************************************************
 
+namespace interpolation
+	  {
+
+//**************************************
+
 //
 // This class provides 1D Lagrange global polynomial interpolation
 // for uniformly spaced data.
 //
 // ...  const  qualifiers refer to results of  interpolate()
 //
-
-namespace interpolation
-	  {
 template <typename fpt>
 class	Lagrange_uniform_1d
 	{
@@ -127,12 +125,21 @@ public:
 	void setup_data_x(const linear_map<fpt>& x_map_in);
 
 	// this function sets up the y data values
+	// n.b. indexing off this pointer is [i], i.e. it need not be 0-origin
 	void setup_data_y(const fpt y_array_in[]);
 
+	// this function sets up
+	// * the x coordinate where interpolation is to be done
+	// * what sort of off-the-end-of-the-data handling is desired
+	void setup_interp_x(fpt x_interp_in, char end_action_in = 'o');
+
 	// the actual interpolation
-	// ... needs  setup_data_x()  and  setup_data_y()
-	//     to have already been called
-	fpt interpolate(fpt x, char end_action = 'o') const;
+	// ... requires that
+	//	setup_data_x()
+	//	setup_data_y()
+	//	setup_interp_x()
+	//     have already been called
+	fpt interpolate();
 
 	// this function computes the backwards domain of dependence
 	// of the interpolation, i.e. the interval of [i] such that
@@ -141,20 +148,24 @@ public:
 	// for  interpolate()
 	//
 	// ... returns domain-of-dependence interval in [i_lo,i_hi]
-	// ... returns nominal center point (i_center in source code)
-	//     as function value
-	// ... needs  setup_data_x()  to have already been called,
-	//     but does *not* need  setup_data_y()  to have already been called
-	int bdod(fpt x, char end_action = 'o', int& i_lo, int& i_hi) const;
-
-	// Jacobian  d interpolate(x,end_action) / d y[i]
-	// ... needs setup_data_x() to have already been called,
-	//     but does *not* need  setup_data_y()  to have already been called
-	fpt Jacobian(fpt x, char end_action = 'o', int i) const;
+	// ... returns nominal center point of underlying local interpolant 
+	//     (i_center in source code) as function value
+	// ... requires that
+	//	setup_data_x()
+	//	setup_interp_x()
+	//     have already been called
+	int bdod(int& i_lo, int& i_hi) const;
+
+	// Jacobian  d (value returned by interpolate()) / d y[i]
+	// ... requires that
+	//	setup_data_x()
+	//	setup_interp_x()
+	//     have already been called
+	fpt Jacobian(int i) const;
 
 	// constructor, destructor
+	// ... compiler-generated no-op destructor is ok
 	Lagrange_uniform_1d(int order_in);
-	~Lagrange_uniform_1d() { }
 
 private:
 	// we forbid copying and passing by value
@@ -166,13 +177,40 @@ private:
 
 private:
 	const int order;
-	const int N_interp;			// number of points in local
-						// interpolator, i.e.
-						// linear = 2 points,
-						// quadratic = 3 points, etc
+	const int N_interp;		// number of points in local
+					// interpolator, i.e.
+					// linear = 2 points,
+					// quadratic = 3 points, etc
+
+	bool x_interp_valid;		// valid flag for x_interp, end_action
+	fpt x_interp;
+	char end_action;
 
 	// n.b. non-const pointers to const data!
 	const linear_map<fpt> *x_map_ptr;	// --> client-owned data
 	const fpt *y_array;			// --> client-owned data
 	};
-	  };
+
+//**************************************
+
+//
+// ***** underlying local interpolants of varying orders *****
+//
+
+// ... x here is in units of the grid spacing,
+//     relative to the center point
+// ... these make no checks for data being in range
+template <typename fpt>
+  fpt interpolate_local_order1(const fpt x, const fpt *y_center);
+template <typename fpt>
+  fpt interpolate_local_order2(const fpt x, const fpt *y_center);
+template <typename fpt>
+  fpt interpolate_local_order3(const fpt x, const fpt *y_center);
+template <typename fpt>
+  fpt interpolate_local_order4(const fpt x, const fpt *y_center);
+fpt interpolate_local_order5(const fpt x, const fpt *y_center);
+fpt interpolate_local_order6(const fpt x, const fpt *y_center);
+
+//**************************************
+
+	  };	// close namespace interpolation::