interp1D.hpp

Go to the documentation of this file.

// Copyright (c) 2024 Andrea Enrico Franzoni (andreaenrico.franzoni@gmail.com)
//
// This file is part of PPCKO
//
// Permission is hereby granted, free of charge, to any person obtaining a copy
// of PPCKO and associated documentation files (the PPCKO software), to deal
// PPCKO without restriction, including without limitation the rights
// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
// copies of PPCKO, and to permit persons to whom PPCKO is
// furnished to do so, subject to the following conditions:
//
// PPCKO IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
// OUT OF OR IN CONNECTION WITH PPCKO OR THE USE OR OTHER DEALINGS IN
// PPCKO.
 
#ifndef EXAMPLES_SRC_INTERP1D_INTERP1D_HPP_
#define EXAMPLES_SRC_INTERP1D_INTERP1D_HPP_
#include <algorithm>
#include <array>
#include <exception>
#include <functional>
#include <iterator>
#include <limits>
#include <type_traits>

 
 
namespace apsc
{
template <typename RAIterator, typename Key, typename ExtractKey,
          typename ExtractValue, typename CompareKey = std::less<Key>>
auto
interp1D(RAIterator const &begin, RAIterator const &end, Key const &keyVal,
         ExtractKey const &extractKey, ExtractValue const &extractValue,
         CompareKey const &comp = std::less<Key>())
{
  // I avoid users using wrong iterators
  // A nice use of iterator_traits and iterator_tags
  using category = typename std::iterator_traits<RAIterator>::iterator_category;
  static_assert(std::is_same_v<category, std::bidirectional_iterator_tag> ||
                  std::is_same_v<category, std::random_access_iterator_tag>,
                "Iterators must be (at least) bidirectional");
  // I need at least two point for interpolating anything. This checks also that
  // end is after begin!
  if(std::distance(begin, end) < 1)
    throw std::runtime_error(
      "Interp1D: I need at least 2 points to interpolate!");
 
  RAIterator a{begin};
  RAIterator b{end};
  // bisection
  for(auto dis = std::distance(a, b); dis > 1;)
    {
      RAIterator c = std::next(a, dis / 2); // midpoint
      if(comp(keyVal, extractKey(*c)))      // keyVal on the left of c
        b = c;
      else // keyVal on the right of c or on c
        a = c;
      dis = std::distance(a, b); // recompute distance
    }
  b = std::next(a, 1); // get other end of interval
  if(b == end)         // complex situation I need to go back of 1
    {
      b = a;
      std::advance(a, -1); // here I need bi-directionality!
    }
  const auto valueLeft = extractValue(*a);
  const Key  keyLeft = extractKey(*a);
  const auto valueRight = extractValue(*b);
  const Key  keyRight = extractKey(*b);
  const auto len = keyRight - keyLeft;
  // I assume no nodes are repeated
  const auto coeffRight = (keyVal - keyLeft) / len;
  const auto coeffLeft = 1.0 - coeffRight;
  return valueLeft * coeffLeft + valueRight * coeffRight;
}
 
} // namespace apsc
 
#endif /* EXAMPLES_SRC_INTERP1D_INTERP1D_HPP_ */