diff --git a/src/util/iceModelVec2T.cc b/src/util/iceModelVec2T.cc index 06f53abc7d..766a2720ae 100644 --- a/src/util/iceModelVec2T.cc +++ b/src/util/iceModelVec2T.cc @@ -747,8 +747,7 @@ void IceModelVec2T::init_interpolation(const std::vector &ts) { &m_data->time[m_data->first], m_data->n_records, times_requested.data(), - times_requested.size(), - m_data->period)); + times_requested.size())); } /** diff --git a/src/util/interpolation.cc b/src/util/interpolation.cc index ddc6cdff25..855c7ad4a5 100644 --- a/src/util/interpolation.cc +++ b/src/util/interpolation.cc @@ -28,17 +28,15 @@ namespace pism { Interpolation::Interpolation(InterpolationType type, const std::vector &input_x, - const std::vector &output_x, - double period) + const std::vector &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 weights_piecewise_constant(const double *x, size_t x_size, double a, diff --git a/src/util/interpolation.hh b/src/util/interpolation.hh index 950cdf0432..6335c784f8 100644 --- a/src/util/interpolation.hh +++ b/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 &input_x, - const std::vector &output_x, double period = 0.0); + const std::vector &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& left() const; const std::vector& 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 integration_weights(const double *x, diff --git a/test/miscellaneous.py b/test/miscellaneous.py index f10e247821..bc67398be9 100644 --- a/test/miscellaneous.py +++ b/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]