SpECTRE Documentation Coverage Report
Current view: top level - NumericalAlgorithms/Interpolation - MultiLinearSpanInterpolation.hpp Hit Total Coverage
Commit: 697db44127ce57424c079f629b3283267ec6bd78 Lines: 15 33 45.5 %
Date: 2024-06-21 21:13:41
Legend: Lines: hit not hit

          Line data    Source code
       1           0 : // Distributed under the MIT License.
       2             : // See LICENSE.txt for details.
       3             : 
       4             : #pragma once
       5             : 
       6             : #include <algorithm>
       7             : #include <array>
       8             : #include <cmath>
       9             : #include <memory>
      10             : 
      11             : #include "DataStructures/Index.hpp"
      12             : #include "Utilities/ConstantExpressions.hpp"
      13             : #include "Utilities/ErrorHandling/Assert.hpp"
      14             : #include "Utilities/GenerateInstantiations.hpp"
      15             : #include "Utilities/Gsl.hpp"
      16             : #include "Utilities/Requires.hpp"
      17             : #include "Utilities/TMPL.hpp"
      18             : 
      19             : namespace intrp {
      20             : 
      21             : /*!
      22             :  * \brief Performs linear interpolation in arbitrary dimensions.
      23             :  *  The class is non-owning and expects a C-ordered array, (n, x, y, z).
      24             :  *  The variable index, n, varies fastest in memory.
      25             :  *  Note that this class is intentionally non-pupable.
      26             :  *
      27             :  *  \tparam Dimension dimensionality of the table
      28             :  *  \tparam NumberOfVariables number of variables stored in the table
      29             :  *  \tparam UniformSpacing indicated whether the table has uniform _spacing or
      30             :  * not. This is useful for performance reasons in finding the appropriate table
      31             :  * index.
      32             :  */
      33             : 
      34             : template <size_t Dimension, size_t NumberOfVariables, bool UniformSpacing>
      35           1 : class MultiLinearSpanInterpolation {
      36             : 
      37             :  public:
      38             :   template <size_t ThisDimension>
      39           0 :   struct Weight {
      40           0 :     std::array<double, two_to_the(ThisDimension)> weights;
      41           0 :     Index<two_to_the(ThisDimension)> index;
      42             :   };
      43             : 
      44           0 :   size_t extents(const size_t which_dimension) const {
      45             :     return number_of_points_[which_dimension];
      46             :   }
      47             : 
      48           0 :   void extrapolate_above_data(const size_t which_dimension, const bool value) {
      49             :     allow_extrapolation_abov_data_[which_dimension] = value;
      50             :   }
      51             : 
      52           0 :   void extrapolate_below_data(const size_t which_dimension, const bool value) {
      53             :     allow_extrapolation_below_data_[which_dimension] = value;
      54             :   }
      55             : 
      56             :   /// Compute interpolation weights for 1D tables
      57           1 :   Weight<1> get_weights(const double x1) const;
      58             : 
      59             :   /// Compute interpolation weights for 2D tables
      60           1 :   Weight<2> get_weights(const double x1, const double x2) const;
      61             : 
      62             :   /// Compute interpolation weights for 3D tables
      63           1 :   Weight<3> get_weights(const double x1, const double x2,
      64             :                         const double x3) const;
      65             : 
      66           0 :   double interpolate(const Weight<Dimension>& weights,
      67             :                      const size_t which_variable = 0) const {
      68             :     double result = 0;
      69             : 
      70             :     for (size_t nn = 0; nn < weights.weights.size(); ++nn) {
      71             :       result += weights.weights[nn] *
      72             :                 y_[which_variable + NumberOfVariables * weights.index[nn]];
      73             :     }
      74             : 
      75             :     return result;
      76             :   }
      77             : 
      78             :   template <size_t... variables_to_interpolate>
      79           0 :   std::array<double, sizeof...(variables_to_interpolate)> interpolate(
      80             :       const Weight<Dimension>& weights) const {
      81             :     static_assert(sizeof...(variables_to_interpolate) <= NumberOfVariables,
      82             :                   "You are trying to interpolate more variables than this "
      83             :                   "container holds.");
      84             : 
      85             :     return std::array<double, sizeof...(variables_to_interpolate)>{
      86             :         interpolate(weights, variables_to_interpolate)...};
      87             :   }
      88             : 
      89             :   template <size_t NumberOfVariablesToInterpolate, typename... T,
      90             :             Requires<(std::is_floating_point_v<typename std::remove_cv_t<T>> and
      91             :                       ...)> = nullptr>
      92           0 :   std::array<double, NumberOfVariablesToInterpolate> interpolate(
      93             :       std::array<size_t, NumberOfVariablesToInterpolate>&
      94             :           variables_to_interpolate,
      95             :       const T&... target_points) const {
      96             :     static_assert(
      97             :         sizeof...(T) == Dimension,
      98             :         "You need to provide the correct number of interpolation points");
      99             :     static_assert(NumberOfVariablesToInterpolate <= NumberOfVariables,
     100             :                   "You are trying to interpolate more variables than this "
     101             :                   "container holds.");
     102             : 
     103             :     auto weights = get_weights(target_points...);
     104             : 
     105             :     std::array<double, NumberOfVariablesToInterpolate> interpolatedValues;
     106             : 
     107             :     for (size_t nn = 0; nn < NumberOfVariablesToInterpolate; ++nn) {
     108             :       interpolatedValues[nn] =
     109             :           interpolate(weights, variables_to_interpolate[nn]);
     110             :     }
     111             : 
     112             :     return interpolatedValues;
     113             :   }
     114             : 
     115             :   template <size_t NumberOfVariablesToInterpolate, size_t... I>
     116           0 :   std::array<double, NumberOfVariablesToInterpolate> interpolate(
     117             :       std::array<size_t, NumberOfVariablesToInterpolate>&
     118             :           variables_to_interpolate,
     119             :       std::array<double, Dimension>& target_points,
     120             :       std::index_sequence<I...>) const {
     121             :     return interpolate(variables_to_interpolate, target_points[I]...);
     122             :   }
     123             : 
     124             :   template <size_t NumberOfVariablesToInterpolate>
     125           0 :   std::array<double, NumberOfVariablesToInterpolate> interpolate(
     126             :       std::array<size_t, NumberOfVariablesToInterpolate>&
     127             :           variables_to_interpolate,
     128             :       std::array<double, Dimension>& target_points) const {
     129             :     return interpolate(variables_to_interpolate, target_points,
     130             :                        std::make_index_sequence<Dimension>{});
     131             :   }
     132             : 
     133           0 :   MultiLinearSpanInterpolation() = default;
     134             : 
     135           0 :   MultiLinearSpanInterpolation(
     136             :       std::array<gsl::span<const double>, Dimension> x_,
     137             :       gsl::span<const double> y_, Index<Dimension> number_of_points__);
     138             : 
     139           0 :   double lower_bound(const size_t which_dimension) const {
     140             :     return x_[which_dimension][0];
     141             :   }
     142             : 
     143           0 :   double upper_bound(const size_t which_dimension) const {
     144             :     return x_[which_dimension][number_of_points_[which_dimension] - 1];
     145             :   }
     146             : 
     147             :  private:
     148             :   /// Allow extrapolation above table bounds.
     149           1 :   std::array<bool, Dimension> allow_extrapolation_below_data_;
     150             :   /// Allow extrapolation below table bounds.
     151           1 :   std::array<bool, Dimension> allow_extrapolation_abov_data_;
     152             :   /// Inverse pacing of the table. Only used for uniform grids
     153           1 :   std::array<double, Dimension> inverse_spacing_;
     154             :   /// Spacing of the table. Only used for uniform grids
     155           1 :   std::array<double, Dimension> spacing_;
     156             :   /// Number of points per dimension
     157           1 :   Index<Dimension> number_of_points_;
     158             : 
     159           0 :   using DataPointer = gsl::span<const double>;
     160             :   /// X values of the table. Only used if allocated
     161           1 :   std::array<DataPointer, Dimension> x_;
     162             :   /// Y values of the table. Only used if allocated
     163           1 :   DataPointer y_;
     164             : 
     165             :   /// Determine relative index bracketing for non-uniform _spacing
     166           1 :   size_t find_index_general(const size_t which_dimension,
     167             :                             const double& target_points) const;
     168             : 
     169             :   /// Determine relative index bracketing for non-uniform _spacing
     170           1 :   size_t find_index_uniform(const size_t which_dimension,
     171             :                             const double& target_points) const;
     172             : 
     173           0 :   size_t find_index(const size_t which_dimension,
     174             :                     const double& target_points) const {
     175             :     if constexpr (UniformSpacing) {
     176             :       return find_index_uniform(which_dimension, target_points);
     177             :     } else {
     178             :       return find_index_general(which_dimension, target_points);
     179             :     }
     180             :   }
     181             : };
     182             : 
     183             : template <size_t Dimension, size_t NumberOfVariables, bool UniformSpacing>
     184             : size_t MultiLinearSpanInterpolation<
     185             :     Dimension, NumberOfVariables,
     186             :     UniformSpacing>::find_index_general(const size_t which_dimension,
     187             :                                         const double& target_points) const {
     188             :   size_t lower_index_bound = 0;
     189             :   size_t upper_index_bound = number_of_points_[which_dimension] - 1;
     190             : 
     191             :   ASSERT((target_points > x_[which_dimension][lower_index_bound]) or
     192             :              allow_extrapolation_below_data_[which_dimension],
     193             :          "Interpolation exceeds lower table bounds");
     194             :   ASSERT((target_points < x_[which_dimension][upper_index_bound]) or
     195             :              allow_extrapolation_abov_data_[which_dimension],
     196             :          "Interpolation exceeds upper table bounds");
     197             : 
     198             :   while (upper_index_bound > 1 + lower_index_bound) {
     199             :     size_t current_index =
     200             :         lower_index_bound + (upper_index_bound - lower_index_bound) / 2;
     201             :     if (target_points < x_[which_dimension][current_index]) {
     202             :       upper_index_bound = current_index;
     203             :     } else {
     204             :       lower_index_bound = current_index;
     205             :     }
     206             :   }
     207             :   return lower_index_bound;
     208             : }
     209             : 
     210             : /// Determine relative index bracketing for non-uniform _spacing
     211             : template <size_t Dimension, size_t NumberOfVariables, bool UniformSpacing>
     212             : size_t MultiLinearSpanInterpolation<
     213             :     Dimension, NumberOfVariables,
     214             :     UniformSpacing>::find_index_uniform(const size_t which_dimension,
     215             :                                         const double& target_points) const {
     216             :   // Compute coordinate relative to lowest table bound
     217             :   const auto relative_coordinate = (target_points - x_[which_dimension][0]);
     218             : 
     219             :   auto current_index = static_cast<size_t>(relative_coordinate *
     220             :                                            inverse_spacing_[which_dimension]);
     221             : 
     222             :   // We are exceeding the table bounds:
     223             :   // Use linear extrapolation based of the lowest
     224             :   // two points in the table
     225             :   ASSERT(allow_extrapolation_below_data_[which_dimension] or
     226             :              UNLIKELY(relative_coordinate >= 0.),
     227             :          "Interpolation exceeds lower table bounds.");
     228             : 
     229             :   // We are exceeding the table bounds:
     230             :   // Use linear extrapolation based of the highest
     231             :   // two points in the table
     232             : 
     233             :   ASSERT(allow_extrapolation_abov_data_[which_dimension] or
     234             :              UNLIKELY(current_index + 1 < number_of_points_[which_dimension]),
     235             :          "Interpolation exceeds upper table bounds.");
     236             : 
     237             :   // Enforce index ranges
     238             :   current_index = std::min(number_of_points_[which_dimension] - 2,
     239             :                            std::max(0ul, current_index));
     240             : 
     241             :   return current_index;
     242             : }
     243             : 
     244             : /// Compute interpolation weights for 1D tables
     245             : template <size_t Dimension, size_t NumberOfVariables, bool UniformSpacing>
     246             : auto MultiLinearSpanInterpolation<Dimension, NumberOfVariables,
     247             :                                   UniformSpacing>::get_weights(const double x1)
     248             :     const -> Weight<1> {
     249             :   // Relative normalized coordinate in x-direction
     250             :   double xx;
     251             : 
     252             :   Weight<1> weights;
     253             : 
     254             :   auto index = Index<1>(find_index(0, x1));
     255             : 
     256             :   if constexpr (UniformSpacing) {
     257             :     xx = (x1 - x_[0][index[0]]) * inverse_spacing_[0];
     258             :   } else {
     259             :     xx = (x1 - x_[0][index[0]]) / (x_[0][index[0] + 1] - x_[0][index[0]]);
     260             :   }
     261             : 
     262             :   // Note: first index varies fastest
     263             :   weights.weights = std::array<double, 2>({
     264             :       (1. - xx),  // 0
     265             :       xx,         // 1
     266             :   });
     267             : 
     268             :   // Compute indices
     269             : 
     270             :   weights.index[0] = index[0];
     271             :   weights.index[1] = index[0] + 1;
     272             : 
     273             :   return weights;
     274             : }
     275             : 
     276             : /// Compute interpolation weights for 2D tables
     277             : 
     278             : template <size_t Dimension, size_t NumberOfVariables, bool UniformSpacing>
     279             : auto MultiLinearSpanInterpolation<Dimension, NumberOfVariables,
     280             :                                   UniformSpacing>::get_weights(const double x1,
     281             :                                                                const double x2)
     282             :     const -> Weight<2> {
     283             :   Weight<2> weights;
     284             : 
     285             :   auto index = Index<2>(find_index(0, x1), find_index(1, x2));
     286             : 
     287             :   // Relative normalized coordinate in x,y-direction
     288             :   double xx, yy;
     289             :   if constexpr (UniformSpacing) {
     290             :     xx = (x1 - x_[0][index[0]]) * inverse_spacing_[0];
     291             :     yy = (x2 - x_[1][index[1]]) * inverse_spacing_[1];
     292             :   } else {
     293             :     xx = (x1 - x_[0][index[0]]) / (x_[0][index[0] + 1] - x_[0][index[0]]);
     294             :     yy = (x2 - x_[1][index[1]]) / (x_[1][index[1] + 1] - x_[1][index[1]]);
     295             :   }
     296             : 
     297             :   // Note: first index varies fastest
     298             :   weights.weights = std::array<double, 4>({
     299             :       (1. - xx) * (1. - yy),  // 00
     300             :       xx * (1. - yy),         // 10
     301             :       (1. - xx) * yy,         // 01
     302             :       xx * yy,                // 11
     303             :   });
     304             : 
     305             :   // Compute indices
     306             :   //
     307             :   for (size_t j = 0; j < 2; ++j) {
     308             :     for (size_t i = 0; i < 2; ++i) {
     309             :       auto tmp_index = Index<2>(index[0] + i, index[1] + j);
     310             :       weights.index[i + 2 * j] = collapsed_index(tmp_index, number_of_points_);
     311             :     }
     312             :   }
     313             : 
     314             :   return weights;
     315             : }
     316             : 
     317             : /// Compute interpolation weights for 3D tables
     318             : template <size_t Dimension, size_t NumberOfVariables, bool UniformSpacing>
     319             : auto MultiLinearSpanInterpolation<Dimension, NumberOfVariables,
     320             :                                   UniformSpacing>::get_weights(const double x1,
     321             :                                                                const double x2,
     322             :                                                                const double x3)
     323             :     const -> Weight<3> {
     324             :   // Relative normalized coordinate in x,y,z-direction
     325             :   double xx, yy, zz;
     326             : 
     327             :   Weight<3> weights;
     328             : 
     329             :   auto index =
     330             :       Index<3>(find_index(0, x1), find_index(1, x2), find_index(2, x3));
     331             : 
     332             :   if constexpr (UniformSpacing) {
     333             :     xx = (x1 - x_[0][index[0]]) * inverse_spacing_[0];
     334             :     yy = (x2 - x_[1][index[1]]) * inverse_spacing_[1];
     335             :     zz = (x3 - x_[2][index[2]]) * inverse_spacing_[2];
     336             :   } else {
     337             :     xx = (x1 - x_[0][index[0]]) / (x_[0][index[0] + 1] - x_[0][index[0]]);
     338             :     yy = (x2 - x_[1][index[1]]) / (x_[1][index[1] + 1] - x_[1][index[1]]);
     339             :     zz = (x3 - x_[2][index[2]]) / (x_[2][index[2] + 1] - x_[2][index[2]]);
     340             :   }
     341             : 
     342             :   // Note: first index varies fastest
     343             :   weights.weights = std::array<double, 8>({
     344             :       (1. - xx) * (1. - yy) * (1. - zz),  // 000
     345             :       xx * (1. - yy) * (1. - zz),         // 100
     346             :       (1. - xx) * yy * (1. - zz),         // 010
     347             :       xx * yy * (1. - zz),                // 110
     348             :       (1. - xx) * (1. - yy) * zz,         // 001
     349             :       xx * (1. - yy) * zz,                // 101
     350             :       (1. - xx) * yy * zz,                // 011
     351             :       xx * yy * zz,                       // 111
     352             :   });
     353             : 
     354             :   // Compute indices
     355             :   //
     356             :   for (size_t k = 0; k < 2; ++k) {
     357             :     for (size_t j = 0; j < 2; ++j) {
     358             :       for (size_t i = 0; i < 2; ++i) {
     359             :         auto tmp_index = Index<3>(index[0] + i, index[1] + j, index[2] + k);
     360             :         weights.index[i + 2 * (j + 2 * k)] =
     361             :             collapsed_index(tmp_index, number_of_points_);
     362             :       }
     363             :     }
     364             :   }
     365             : 
     366             :   return weights;
     367             : }
     368             : 
     369             : template <size_t Dimension, size_t NumberOfVariables, bool UniformSpacing>
     370             : MultiLinearSpanInterpolation<Dimension, NumberOfVariables, UniformSpacing>::
     371             :     MultiLinearSpanInterpolation(
     372             :         std::array<gsl::span<double const>, Dimension> x,
     373             :         gsl::span<double const> y, Index<Dimension> number_of_points)
     374             :     : number_of_points_(number_of_points), x_(x), y_(y) {
     375             :   for (size_t i = 0; i < Dimension; ++i) {
     376             :     spacing_[i] = x_[i][1] - x_[i][0];
     377             :     inverse_spacing_[i] = 1. / spacing_[i];
     378             :     allow_extrapolation_below_data_[i] = false;
     379             :     allow_extrapolation_abov_data_[i] = false;
     380             :   }
     381             : }
     382             : 
     383             : /// Multilinear span interpolation with uniform grid spacing
     384             : template <size_t Dimension, size_t NumberOfVariables>
     385           1 : using UniformMultiLinearSpanInterpolation =
     386             :     MultiLinearSpanInterpolation<Dimension, NumberOfVariables, true>;
     387             : 
     388             : /// Multilinear span interpolation with non-uniform grid spacing
     389             : template <size_t Dimension, size_t NumberOfVariables>
     390           1 : using GeneralMultiLinearSpanInterpolation =
     391             :     MultiLinearSpanInterpolation<Dimension, NumberOfVariables, false>;
     392             : 
     393             : }  // namespace intrp

Generated by: LCOV version 1.14