Line data    Source code

       1             : /*
       2             :  * interpolation.hpp
       3             :  *
       4             :  * Groups functions to do linear interpolation
       5             :  *
       6             :  * Copyright (c) 2019, The Chancellor, Masters and Scholars of the University
       7             :  * of Oxford, VITO nv, and the 'Slide' Developers.
       8             :  * See the licence file LICENCE.txt for more information.
       9             :  */
      10             : 
      11             : #pragma once
      12             : 
      13             : #include "../settings/settings.hpp"
      14             : 
      15             : #include <vector>
      16             : #include <array>
      17             : #include <iostream>
      18             : #include <cassert>
      19             : #include <algorithm>
      20             : #include <utility>
      21             : #include <cmath>
      22             : #include <cstdlib>
      23             : 
      24             : namespace slide {
      25             : template <typename Tx, typename Ty>
      26     2888602 : auto linInt_noexcept(bool bound, Tx &xdat, Ty &ydat, int nin, double x, bool is_fixed = false)
      27             : {
      28             :   /*
      29             :    * function for linear interpolation with the data points provided as two arrays
      30             :    *
      31             :    * IN
      32             :    * bound      boolean deciding what to do if the value of x is out of range of xdat
      33             :    *                    if true, an error is thrown
      34             :    *                    if false, the value closest to x is returned
      35             :    * xdat       x data points in strictly increasing order
      36             :    * ydat       y data points
      37             :    * nin                number of data points
      38             :    * x          x point at which value is needed
      39             :    * is_fixed If the difference between values are fixed.
      40             :    *
      41             :    * OUT
      42             :    * y          y value corresponding to x
      43             :    * status     0 if successful, 1 if x>first and -1 if x < last
      44             :    */
      45             : 
      46     2888602 :   double yy{ 0.0 }; //!< Some programs depend on 0.0 initial condition when status != 0, do not change.
      47     2888602 :   int status = 0;   //!< Set the status as inverse of bound. So that first two branches of if are invalid if bound is true.
      48             :   //!< check that x is within the limits of the data points
      49             : 
      50     2888602 :   if (bound) {
      51     2888602 :     if (x < xdat[0]) {
      52           0 :       status = -1;
      53           0 :       return std::make_pair(yy, status);
      54     2888602 :     } else if (x > xdat[nin - 1]) {
      55           0 :       status = 1;
      56           0 :       return std::make_pair(yy, status);
      57             :     }
      58             :   }
      59             : 
      60     2888602 :   if (x <= xdat[0])
      61           0 :     yy = ydat[0]; //!< if x is below the minimum value, return the y value of the first data point
      62     2888602 :   else if (x >= xdat[nin - 1])
      63           1 :     yy = ydat[nin - 1]; //!< if x is above the maximum value, return the y value of the last data point
      64             :   else {
      65             :     //!< scan the data points
      66     2888601 :     int i_low{ 0 };
      67             :     //!< For fixed step no need to iterate:
      68             : 
      69     2888601 :     if (is_fixed) {
      70     2888601 :       double dt = xdat[1] - xdat[0];
      71     2888601 :       i_low = static_cast<int>((x - xdat[0]) / dt) + 1;
      72             :     } else {
      73             :       //!< binary search algorithm:
      74             :       //!< i_low will be the first index which compares greater than x;
      75             :       //!< const auto it = std::find_if(std::begin(xdat), std::begin(xdat) + nin, [x](double element) { return (x < element); }); -> Linear search if needed.
      76           0 :       const auto it = std::lower_bound(xdat.begin(), xdat.begin() + nin, x);
      77           0 :       i_low = static_cast<int>(it - xdat.begin());
      78             :     }
      79             : 
      80     2888601 :     const double xr = xdat[i_low]; //!< then that point is the point 'to the right' of x
      81     2888601 :     const double yr = ydat[i_low];
      82     2888601 :     const double xl = xdat[i_low - 1]; //!< while the previous point is the point 'to the left' of x
      83     2888601 :     const double yl = ydat[i_low - 1];
      84     2888601 :     yy = yl + (yr - yl) * (x - xl) / (xr - xl);
      85             :   }
      86             : 
      87     2888602 :   return std::make_pair(yy, status);
      88             : }
      89             : 
      90             : template <typename Tx, typename Ty>
      91     2888602 : double linInt(bool verbose, bool bound, Tx &xdat, Ty &ydat, int nin, double x, bool is_fixed = false)
      92             : {
      93             :   /*
      94             :    * function for linear interpolation with the data points provided as two arrays
      95             :    *
      96             :    * IN
      97             :    * verbose    if false, no error message is written (but the error is still thrown)
      98             :    *                    if true, an error message is written
      99             :    * bound      boolean deciding what to do if the value of x is out of range of xdat
     100             :    *                    if true, an error is thrown
     101             :    *                    if false, the value closest to x is returned
     102             :    * xdat       x data points in strictly increasing order
     103             :    * ydat       y data points
     104             :    * nin                number of data points
     105             :    * x          x point at which value is needed
     106             :    * is_fixed If the difference between values are fixed.
     107             :    *
     108             :    * OUT
     109             :    * y          y value corresponding to x
     110             :    *
     111             :    * THROWS
     112             :    * 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
     113             :    *                            i.e. bound AND (x < xdat [0] OR x > xdat[end])
     114             :    */
     115             : 
     116     2888602 :   auto [yy, status] = linInt_noexcept(bound, xdat, ydat, nin, x, is_fixed);
     117             : 
     118     2888602 :   if (status) {
     119           0 :     if (verbose)
     120           0 :       std::cerr << "ERROR in Interpolation::linInt: x is out of bounds. x = " << x << " while xmin = " << xdat[0] << " and xmax is " << xdat[nin - 1] << ".\n";
     121           0 :     throw 1;
     122             :   }
     123             : 
     124     2888602 :   return yy;
     125             : }
     126             : 
     127             : } // namespace slide