interpolation.hpp

Go to the documentation of this file.

/*
 * interpolation.hpp
 *
 * Groups functions to do linear interpolation
 *
 * Copyright (c) 2019, The Chancellor, Masters and Scholars of the University
 * of Oxford, VITO nv, and the 'Slide' Developers.
 * See the licence file LICENCE.txt for more information.
 */
 
#pragma once
 
#include "../settings/settings.hpp"
 
#include <vector>
#include <array>
#include <iostream>
#include <cassert>
#include <algorithm>
#include <utility>
#include <cmath>
#include <cstdlib>
 
namespace slide {
template <typename Tx, typename Ty>
auto linInt_noexcept(bool bound, Tx &xdat, Ty &ydat, int nin, double x, bool is_fixed = false)
{
  /*
   * function for linear interpolation with the data points provided as two arrays
   *
   * IN
   * bound  boolean deciding what to do if the value of x is out of range of xdat
   *            if true, an error is thrown
   *            if false, the value closest to x is returned
   * xdat   x data points in strictly increasing order
   * ydat   y data points
   * nin        number of data points
   * x      x point at which value is needed
   * is_fixed If the difference between values are fixed.
   *
   * OUT
   * y      y value corresponding to x
   * status     0 if successful, 1 if x>first and -1 if x < last
   */
 
  double yy{ 0.0 }; 
  int status = 0;   
 
  if (bound) {
    if (x < xdat[0]) {
      status = -1;
      return std::make_pair(yy, status);
    } else if (x > xdat[nin - 1]) {
      status = 1;
      return std::make_pair(yy, status);
    }
  }
 
  if (x <= xdat[0])
    yy = ydat[0]; 
  else if (x >= xdat[nin - 1])
    yy = ydat[nin - 1]; 
  else {
    int i_low{ 0 };
 
    if (is_fixed) {
      double dt = xdat[1] - xdat[0];
      i_low = static_cast<int>((x - xdat[0]) / dt) + 1;
    } else {
      const auto it = std::lower_bound(xdat.begin(), xdat.begin() + nin, x);
      i_low = static_cast<int>(it - xdat.begin());
    }
 
    const double xr = xdat[i_low]; 
    const double yr = ydat[i_low];
    const double xl = xdat[i_low - 1]; 
    const double yl = ydat[i_low - 1];
    yy = yl + (yr - yl) * (x - xl) / (xr - xl);
  }
 
  return std::make_pair(yy, status);
}
 
template <typename Tx, typename Ty>
double linInt(bool verbose, bool bound, Tx &xdat, Ty &ydat, int nin, double x, bool is_fixed = false)
{
  /*
   * function for linear interpolation with the data points provided as two arrays
   *
   * IN
   * verbose    if false, no error message is written (but the error is still thrown)
   *            if true, an error message is written
   * bound  boolean deciding what to do if the value of x is out of range of xdat
   *            if true, an error is thrown
   *            if false, the value closest to x is returned
   * xdat   x data points in strictly increasing order
   * ydat   y data points
   * nin        number of data points
   * x      x point at which value is needed
   * is_fixed If the difference between values are fixed.
   *
   * OUT
   * y      y value corresponding to x
   *
   * THROWS
   * 1      if bound = true && if x is out of bounds, i.e. smaller than the smallest value of xdat or larger than the largest value of xdat
   *                i.e. bound AND (x < xdat [0] OR x > xdat[end])
   */
 
  auto [yy, status] = linInt_noexcept(bound, xdat, ydat, nin, x, is_fixed);
 
  if (status) {
    if (verbose)
      std::cerr << "ERROR in Interpolation::linInt: x is out of bounds. x = " << x << " while xmin = " << xdat[0] << " and xmax is " << xdat[nin - 1] << ".\n";
    throw 1;
  }
 
  return yy;
}
 
} // namespace slide