Remove old code for piece-wise linear periodic interpolation

The new approach is better.

src/util/iceModelVec2T.cc

@@ -747,8 +747,7 @@ void IceModelVec2T::init_interpolation(const std::vector<double> &ts) {
                                          &m_data->time[m_data->first],
                                          m_data->n_records,
                                          times_requested.data(),
-                                         times_requested.size(),
-                                         m_data->period));
+                                         times_requested.size()));
 }
 
 /**

src/util/interpolation.cc

@@ -28,17 +28,15 @@ namespace pism {
 
 Interpolation::Interpolation(InterpolationType type,
                              const std::vector<double> &input_x,
-                             const std::vector<double> &output_x,
-                             double period)
+                             const std::vector<double> &output_x)
   : Interpolation(type, input_x.data(), input_x.size(),
-                  output_x.data(), output_x.size(), period) {
+                  output_x.data(), output_x.size()) {
   // empty
 }
 
 Interpolation::Interpolation(InterpolationType type,
                              const double *input_x, unsigned int input_x_size,
-                             const double *output_x, unsigned int output_x_size,
-                             double period) {
+                             const double *output_x, unsigned int output_x_size) {
 
   // the trivial case (the code below requires input_x_size >= 2)
   if (input_x_size < 2) {
@@ -74,9 +72,6 @@ Interpolation::Interpolation(InterpolationType type,
   case PIECEWISE_CONSTANT:
     init_piecewise_constant(input_x, input_x_size, output_x, output_x_size);
     break;
-  case LINEAR_PERIODIC:
-    init_linear_periodic(input_x, input_x_size, output_x, output_x_size, period);
-    break;
     // LCOV_EXCL_START
   default:
     throw RuntimeError(PISM_ERROR_LOCATION, "invalid interpolation type");
@@ -210,63 +205,6 @@ void Interpolation::init_piecewise_constant(const double *input_x, unsigned int
   }
 }
 
-void Interpolation::init_linear_periodic(const double *input_x, unsigned int input_x_size,
-                                         const double *output_x, unsigned int output_x_size,
-                                         double period) {
-
-  assert(period > 0);
-  assert(input_x_size >= 2);
-
-  m_left.resize(output_x_size);
-  m_right.resize(output_x_size);
-  m_alpha.resize(output_x_size);
-
-  // compute indexes and weights
-  for (unsigned int i = 0; i < output_x_size; ++i) {
-    double x = output_x[i];
-
-    unsigned int L = 0, R = 0;
-    if (x < input_x[0]) {
-      L = input_x_size - 1;
-      R = 0.0;
-    } else {
-      // note: use "input_x_size" instead of "input_x_size - 1" to support extrapolation on
-      // the right
-      L = gsl_interp_bsearch(input_x, x, 0, input_x_size);
-      R = L + 1 < input_x_size ? L + 1 : 0;
-    }
-
-    double
-      x_l = input_x[L],
-      x_r = input_x[R],
-      alpha = 0.0;
-    if (L < R) {
-      // regular case
-      alpha = (x - x_l) / (x_r - x_l);
-    } else {
-      double
-        x0 = input_x[0],
-        dx = (period - x_l) + x0;
-      assert(dx > 0);
-      if (x > x0) {
-        // interval from the last point of the input grid to the period
-        alpha = (x - x_l) / dx;
-      } else {
-        // interval from 0 to the first point of the input grid
-        alpha = 1.0 - (x_r - x) / dx;
-      }
-    }
-
-    assert(L < input_x_size);
-    assert(R < input_x_size);
-    assert(alpha >= 0.0 and alpha <= 1.0);
-
-    m_left[i]  = L;
-    m_right[i] = R;
-    m_alpha[i] = alpha;
-  }
-}
-
 static std::map<size_t, double> weights_piecewise_constant(const double *x,
                                                            size_t x_size,
                                                            double a,

src/util/interpolation.hh

@@ -25,7 +25,7 @@
 
 namespace pism {
 
-enum InterpolationType {LINEAR, NEAREST, PIECEWISE_CONSTANT, LINEAR_PERIODIC};
+enum InterpolationType {LINEAR, NEAREST, PIECEWISE_CONSTANT};
 
 /**
  * Class encapsulating linear and piece-wise constant interpolation indexes and weights.
@@ -76,18 +76,13 @@ enum InterpolationType {LINEAR, NEAREST, PIECEWISE_CONSTANT, LINEAR_PERIODIC};
  * intervals: [0, 1) and [1, x_e). Here the value x_e is irrelevant because we use
  * constant extrapolation for points outside the interval covered by input data.
  *
- *
- * Linear interpolation for periodic data (annual temperature cycles, etc).
- *
- * Input data are assumed to be periodic on the interval from 0 to `period`.
- *
  */
 class Interpolation {
 public:
   Interpolation(InterpolationType type, const std::vector<double> &input_x,
-                const std::vector<double> &output_x, double period = 0.0);
+                const std::vector<double> &output_x);
   Interpolation(InterpolationType type, const double *input_x, unsigned int input_x_size,
-                const double *output_x, unsigned int output_x_size, double period = 0.0);
+                const double *output_x, unsigned int output_x_size);
 
   const std::vector<int>& left() const;
   const std::vector<int>& right() const;
@@ -119,9 +114,6 @@ private:
                     const double *output_x, unsigned int output_x_size);
   void init_piecewise_constant(const double *input_x, unsigned int input_x_size,
                                const double *output_x, unsigned int output_x_size);
-  void init_linear_periodic(const double *input_x, unsigned int input_x_size,
-                            const double *output_x, unsigned int output_x_size,
-                            double period);
 };
 
 std::map<size_t, double> integration_weights(const double *x,

test/miscellaneous.py

@@ -1077,7 +1077,7 @@ def test_interpolation_other():
 
     try:
         PISM.Interpolation(PISM.PIECEWISE_CONSTANT, [2, 1], [2, 1])
-        assert False, "failed to detect non-increasing data"
+        assert False, "failed to detect non-increasing times"
     except RuntimeError as e:
         print(e)
         pass
@@ -1101,27 +1101,6 @@ def test_interpolation_other():
     np.testing.assert_almost_equal(I.right(), [1, 2])
     np.testing.assert_almost_equal(I.alpha(), [0.5, 0.5])
 
-
-def test_linear_periodic():
-    "Linear (periodic) interpolation"
-
-    period = 1.0
-    x_p = np.linspace(0, 0.9, 10) + 0.05
-    y_p = np.sin(2 * np.pi * x_p)
-
-    # grid with end points added (this makes periodic interpolation unnecessary)
-    x = np.r_[0, x_p, 1]
-    y = np.sin(2 * np.pi * x)
-
-    # target grid
-    xx = np.linspace(0, 1, 21)
-
-    yy_p = PISM.Interpolation(PISM.LINEAR_PERIODIC, x_p, xx, period).interpolate(y_p)
-
-    yy = PISM.Interpolation(PISM.LINEAR, x, xx).interpolate(y)
-
-    np.testing.assert_almost_equal(yy, yy_p)
-
 def test_nearest_neighbor():
     "Nearest neighbor interpolation"
     x = [-1, 1]