Line data    Source code

       1             : //* This file is part of the MOOSE framework
       2             : //* https://mooseframework.inl.gov
       3             : //*
       4             : //* All rights reserved, see COPYRIGHT for full restrictions
       5             : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
       6             : //*
       7             : //* Licensed under LGPL 2.1, please see LICENSE for details
       8             : //* https://www.gnu.org/licenses/lgpl-2.1.html
       9             : 
      10             : #include "MultiDimensionalInterpolation.h"
      11             : 
      12             : #include <cassert>
      13             : #include <fstream>
      14             : #include <stdexcept>
      15             : 
      16             : template <typename T>
      17           4 : MultiDimensionalInterpolationTempl<T>::MultiDimensionalInterpolationTempl(
      18             :     const std::vector<std::vector<Real>> & base_points, const MultiIndex<Real> & data)
      19           4 :   : _base_points(std::vector<std::vector<Real>>()), _data(MultiIndex<Real>({}))
      20             : {
      21           4 :   setData(base_points, data);
      22           4 : }
      23             : 
      24             : template <typename T>
      25           1 : MultiDimensionalInterpolationTempl<T>::MultiDimensionalInterpolationTempl()
      26           1 :   : _base_points(std::vector<std::vector<Real>>()), _data(MultiIndex<Real>({}))
      27             : {
      28           1 : }
      29             : 
      30             : template <typename T>
      31             : void
      32           6 : MultiDimensionalInterpolationTempl<T>::setData(const std::vector<std::vector<Real>> & base_points,
      33             :                                                const MultiIndex<Real> & data)
      34             : {
      35           6 :   _degenerate_interpolation = false;
      36           6 :   _setup_complete = false;
      37             : 
      38             :   // stash away original dimension of data because we might modify it
      39           6 :   _original_dim = data.dim();
      40           6 :   _degenerate_index.resize(_original_dim, false);
      41             : 
      42             :   // need to check this here to make sure the next loop is safe
      43           6 :   if (data.dim() != base_points.size())
      44             :   {
      45           0 :     std::ostringstream oss;
      46           0 :     oss << "Leading dimension of base_points " << base_points.size() << " and data dimensionality "
      47           0 :         << data.dim() << " are inconsistent.";
      48           0 :     throw std::domain_error(oss.str());
      49           0 :   }
      50             : 
      51             :   // first a corner case is addressed; if the size of the MultiIndex in one
      52             :   // dimension is 1, then finding the interpolation "stencil" becomes much
      53             :   // harder. It is easier to slice the data array and keep track of which values
      54             :   // to discard in the interpolation routines
      55           6 :   MultiIndex<Real>::size_type slice_dimensions;
      56           6 :   MultiIndex<Real>::size_type slice_indices;
      57           6 :   MultiIndex<Real>::size_type size = data.size();
      58           6 :   std::vector<std::vector<Real>> modified_base_points;
      59          28 :   for (unsigned int j = 0; j < data.dim(); ++j)
      60          22 :     if (size[j] == 1)
      61             :     {
      62           7 :       slice_dimensions.push_back(j);
      63           7 :       slice_indices.push_back(0);
      64           7 :       _degenerate_index[j] = true;
      65             :     }
      66             :     else
      67          15 :       modified_base_points.push_back(base_points[j]);
      68             : 
      69             :   // if all dimensions have size 1, then there is only one value to return
      70             :   // this is a weird corner case that we should support to make life easier
      71             :   // for whoever uses this object
      72             :   // In this case, we can just use the passed base_points and data objects
      73             :   // and short-cut the computation in the interpolation function
      74           6 :   if (modified_base_points.size() == 0)
      75             :   {
      76           1 :     _degenerate_interpolation = true;
      77           1 :     _base_points = base_points;
      78           1 :     _data = data;
      79             :   }
      80             :   else
      81             :   {
      82             :     // set the data
      83           5 :     _base_points = modified_base_points;
      84           5 :     _data = data.slice(slice_dimensions, slice_indices);
      85             :   }
      86           6 :   errorCheck();
      87           6 : }
      88             : 
      89             : template <typename T>
      90             : void
      91           6 : MultiDimensionalInterpolationTempl<T>::errorCheck()
      92             : {
      93           6 :   if (_base_points.size() != _data.dim())
      94           0 :     throw std::domain_error("Size of the leading dimension of base points must be equal to "
      95             :                             "dimensionality of data array");
      96             : 
      97           6 :   MultiIndex<Real>::size_type shape = _data.size();
      98          26 :   for (unsigned int j = 0; j < _data.dim(); ++j)
      99             :   {
     100          20 :     if (shape[j] != _base_points[j].size())
     101             :     {
     102           0 :       std::ostringstream oss;
     103           0 :       oss << "Dimension " << j << " has " << _base_points[j].size()
     104           0 :           << " base points but dimension of data array is " << shape[j];
     105           0 :       throw std::domain_error(oss.str());
     106           0 :     }
     107             : 
     108          65 :     for (unsigned int i = 1; i < _base_points[j].size(); ++i)
     109             :     {
     110          45 :       if (_base_points[j][i - 1] >= _base_points[j][i])
     111             :       {
     112           0 :         std::ostringstream oss;
     113           0 :         oss << "Base point values for dimension " << j << " are not strictly increasing: bp["
     114           0 :             << i - 1 << "]: " << _base_points[j][i - 1] << " bc[" << i
     115           0 :             << "]: " << _base_points[j][i];
     116           0 :         throw std::domain_error(oss.str());
     117           0 :       }
     118             :     }
     119             :   }
     120             : 
     121             :   // now the setup is complete
     122           6 :   _setup_complete = true;
     123           6 : }
     124             : 
     125             : template <typename T>
     126             : void
     127        2667 : MultiDimensionalInterpolationTempl<T>::linearSearch(std::vector<T> & values,
     128             :                                                     MultiIndex<Real>::size_type & indices) const
     129             : {
     130             :   assert(values.size() == _data.dim());
     131             : 
     132        2667 :   indices.resize(_data.dim());
     133       10668 :   for (unsigned int d = 0; d < _data.dim(); ++d)
     134        8001 :     indices[d] = linearSearchHelper(values[d], _base_points[d]);
     135        2667 : }
     136             : 
     137             : template <typename T>
     138             : unsigned int
     139        8001 : MultiDimensionalInterpolationTempl<T>::linearSearchHelper(T & x,
     140             :                                                           const std::vector<Real> & vector) const
     141             : {
     142             : // check order of vector only in debug
     143             : #ifndef NDEBUG
     144             :   for (unsigned int i = 1; i < vector.size(); ++i)
     145             :   {
     146             :     if (vector[i - 1] >= vector[i])
     147             :     {
     148             :       std::ostringstream oss;
     149             :       oss << "In linearSearchHelper vector is not strictly increasing: vector[" << i - 1
     150             :           << "]: " << vector[i - 1] << " vector[" << i << "]: " << vector[i];
     151             :       throw std::domain_error(oss.str());
     152             :     }
     153             :   }
     154             : #endif
     155             : 
     156             :   // smaller/larger cases for tabulation
     157        8001 :   if (x < vector[0])
     158             :   {
     159           2 :     x = vector[0];
     160           2 :     return 0;
     161             :   }
     162        7999 :   else if (x >= vector.back())
     163             :   {
     164             :     // this requires explanation because it looks like a bug (I swear it isn't, it's a feature)
     165             :     // if x >= largest base point value and the index i = vector.size() - 1 is returned, then during
     166             :     // the interpolation, i + 1 will be accessed, causing a segementation fault. Logic in the
     167             :     // interpolation is also bad because it may slow it down. By returning i = vector.size() - 2,
     168             :     // the last two values in the base point array are used, but x is reset to be equal to the last
     169             :     // value. This will give the right interpolation behavior.
     170         610 :     x = vector.back();
     171         610 :     return vector.size() - 2;
     172             :   }
     173             : 
     174       15499 :   for (unsigned int j = 0; j < vector.size() - 2; ++j)
     175       12230 :     if (x >= vector[j] && x < vector[j + 1])
     176        4120 :       return j;
     177        3269 :   return vector.size() - 2;
     178             : }
     179             : 
     180             : template <typename T>
     181             : T
     182        2664 : MultiDimensionalInterpolationTempl<T>::multiLinearInterpolation(const std::vector<T> & x) const
     183             : {
     184             :   // ensure that object has been set up properly
     185        2664 :   if (!_setup_complete)
     186           0 :     throw std::domain_error(
     187             :         "Object has been constructed properly. This happens when the empty constructor is called "
     188             :         "and setData is not called before interpolation.");
     189             : 
     190             :   // assert does not seem to be enough to check data consistency here
     191             :   // because this function is exposed to user and will be frequently used
     192        2664 :   if (x.size() != _original_dim)
     193             :   {
     194           0 :     std::ostringstream oss;
     195           0 :     oss << "In sample the parameter x has size " << x.size() << " but data has dimension "
     196           0 :         << _data.dim();
     197           0 :     throw std::domain_error(oss.str());
     198           0 :   }
     199             : 
     200             :   // in this case there is only one entry that resides at flat index 0
     201        2664 :   if (_degenerate_interpolation)
     202           1 :     return _data[0];
     203             : 
     204             :   // make a copy of x because linearSearch needs to be able to modify it
     205             :   // also adjust for degenerate indices that were removed
     206        2663 :   std::vector<T> y(_data.dim());
     207        2663 :   unsigned int l = 0;
     208       10654 :   for (unsigned int j = 0; j < _original_dim; ++j)
     209        7991 :     if (!_degenerate_index[j])
     210             :     {
     211        7989 :       y[l] = x[j];
     212        7989 :       ++l;
     213             :     }
     214             : 
     215             :   // obtain the indices using linearSearch & get flat index in _data array
     216        2663 :   MultiIndex<Real>::size_type indices;
     217        2663 :   linearSearch(y, indices);
     218             : 
     219             :   // pre-compute the volume of the hypercube and weights used for the interpolation
     220        2663 :   Real volume = 1;
     221        2663 :   std::vector<std::vector<T>> weights(_data.dim());
     222       10652 :   for (unsigned int j = 0; j < _data.dim(); ++j)
     223             :   {
     224        7989 :     volume *= _base_points[j][indices[j] + 1] - _base_points[j][indices[j]];
     225             : 
     226             :     // now compute weight for each dimension, note that
     227             :     // weight[j][0] = bp[j][high] - y[j]
     228             :     // weight[j][0] = y[j] - bp[j][low]
     229             :     // because in the interpolation the value on the "left" is weighted with the
     230             :     // distance of y to the "right"
     231        7989 :     weights[j].resize(2);
     232        7989 :     weights[j][0] = _base_points[j][indices[j] + 1] - y[j];
     233        7989 :     weights[j][1] = y[j] - _base_points[j][indices[j]];
     234             :   }
     235             : 
     236             :   // get flat index and stride
     237        2663 :   unsigned int flat_index = _data.flatIndex(indices);
     238        2663 :   MultiIndex<Real>::size_type stride = _data.stride();
     239             : 
     240             :   // this is the interpolation routine, evaluation can be expensive if dim is large
     241             : 
     242             :   // first we retrieve the data around the point y
     243             :   // we use a 2^dim MultiIndex hypercube for generating the permulations
     244        2663 :   MultiIndex<Real>::size_type shape(_data.dim(), 2);
     245        2663 :   MultiIndex<Real> hypercube = MultiIndex<Real>(shape);
     246        2663 :   T interpolated_value = 0;
     247       26630 :   for (const auto & p : hypercube)
     248             :   {
     249             :     // this loop computes the flat index for the point in the
     250             :     // hypercube and accumulates the total weight for this point
     251             :     // the selection of the dimensional weight is based on the swapping
     252             :     // done in the pre-computing loop before
     253       21304 :     unsigned int index = flat_index;
     254       21304 :     T w = 1;
     255       85216 :     for (unsigned int j = 0; j < p.first.size(); ++j)
     256             :     {
     257       63912 :       index += stride[j] * p.first[j];
     258       63912 :       w *= weights[j][p.first[j]];
     259             :     }
     260             : 
     261             :     // add the contribution of this base point to the interpolated value
     262       21304 :     interpolated_value += w * _data[index];
     263             :   }
     264             : 
     265        2663 :   return interpolated_value / volume;
     266        2663 : }
     267             : 
     268             : template class MultiDimensionalInterpolationTempl<Real>;
     269             : template class MultiDimensionalInterpolationTempl<ADReal>;