Line data    Source code

       1           0 : // Distributed under the MIT License.
       2             : // See LICENSE.txt for details.
       3             : 
       4             : #pragma once
       5             : 
       6             : #include <array>
       7             : #include <cstddef>
       8             : #include <gsl/gsl_multifit.h>
       9             : #include <vector>
      10             : 
      11             : /// \cond
      12             : namespace PUP {
      13             : class er;
      14             : }  // namespace PUP
      15             : /// \endcond
      16             : 
      17             : namespace intrp {
      18             : /*!
      19             :  * \ingroup NumericalAlgorithmsGroup
      20             :  * \brief A linear least squares solver class
      21             :  *
      22             :  * A wrapper class for the gsl linear least squares solver which determines
      23             :  * the coefficients of best fit for a polynomial of order `Order` given a set
      24             :  * of data points `x_values` and `y_values` representing some function y(x).
      25             :  * Note that the `interpolate` function requires a set of fit coefficients,
      26             :  * which can be obtained by first calling the `fit_coefficients` function.
      27             :  * The parameter `num_observations` refers to the number of entries in
      28             :  * `x_values`.
      29             :  *
      30             :  * The details of the linear least squares solver can be seen here:
      31             :  * [GSL documentation](https://www.gnu.org/software/gsl/doc/html/lls.html#).
      32             :  *
      33             :  * \warning This class is not thread-safe, the reason for this is because
      34             :  * making the buffers used in the class thread-local would be a large use
      35             :  * of memory. This class should be treated as a function that internally
      36             :  * stores workspace variables to reduce memory allocations.
      37             :  */
      38             : 
      39             : template <size_t Order>
      40           1 : class LinearLeastSquares {
      41             :  public:
      42           0 :   LinearLeastSquares(const size_t num_observations);
      43             : 
      44           0 :   LinearLeastSquares() = default;
      45           0 :   LinearLeastSquares(const LinearLeastSquares& /*rhs*/) = delete;
      46           0 :   LinearLeastSquares& operator=(const LinearLeastSquares& /*rhs*/) = delete;
      47           0 :   LinearLeastSquares(LinearLeastSquares&& /*rhs*/) = default;
      48           0 :   LinearLeastSquares& operator=(LinearLeastSquares&& rhs) = default;
      49           0 :   ~LinearLeastSquares();
      50             : 
      51           0 :   double interpolate(const std::array<double, Order + 1> coefficients,
      52             :                      const double x_to_interp_to);
      53             :   template <typename T>
      54           0 :   std::array<double, Order + 1> fit_coefficients(const T& x_values,
      55             :                                                  const T& y_values);
      56             : 
      57             :   /*!
      58             :    * The `x_values` are a sequence of times, positions, or other abscissa.
      59             :    * The `y_values` are a std::vector of sequences of ordinates corresponding
      60             :    * to these abscissa, with each element of the
      61             :    * std::vector containing one sequence. Therefore, one set of fit
      62             :    * coefficients is found for each entry in the std::vector `y_values`.
      63             :    * In other words, the data in `y_values` is stored as:
      64             :    * `y_values[fit_index][x_value_index]`, where `fit_index` runs over the
      65             :    * number of different sequences being fit to, and `x_value_index` runs over
      66             :    * the different entries corrsponding to those in `x_values`.
      67             :    */
      68             :   template <typename T>
      69           1 :   std::vector<std::array<double, Order + 1>> fit_coefficients(
      70             :       const T& x_values, const std::vector<T>& y_values);
      71             : 
      72             :   // NOLINTNEXTLINE(google-runtime-references)
      73           0 :   void pup(PUP::er& p);
      74             : 
      75             :  private:
      76           0 :   size_t num_observations_;
      77           0 :   gsl_matrix* X = nullptr;
      78           0 :   gsl_matrix* covariance_matrix_ = nullptr;
      79           0 :   gsl_vector* y = nullptr;
      80           0 :   gsl_vector* c = nullptr;
      81             : };
      82             : }  // namespace intrp