KnotArray.h

// -*- C++ -*-
//
// This file is part of LHAPDF
// Copyright (C) 2012-2024 The LHAPDF collaboration (see AUTHORS for details)
//
#pragma once
#ifndef LHAPDF_KnotArray_H
#define LHAPDF_KnotArray_H
 
#include "LHAPDF/Exceptions.h"
#include "LHAPDF/Utils.h"
#include <algorithm>  // for std::upper_bound, std::find
 
namespace {
 
 
  // Hide some internal functions from outside API view
 
  // General function to find the knot below a given value
  // PERFORMANCE-CRITICAL: This function is very hot. Keep optimized for speed. When upgrading to C++20, add [[unlikely]] to the edge case branch.
  inline size_t indexbelow(double value, const std::vector<double>& knots) {
    const size_t n = knots.size();
    const double* b = knots.data();
    const double* e = b + n;
 
    size_t i = std::upper_bound(b, e, value) - b;
    if (i >= n) i = n - 1;
    return i - 1;
  }
 
  int findPidInPids(int pid, const std::vector<int>& pids) {
    std::vector<int>::const_iterator it = std::find(pids.begin(), pids.end(), pid);
    if (it == pids.end())
      return -1;
    else
      // save to cast, given small number of pid's
      return static_cast<int>(std::distance(pids.begin(), it));
  }
 
 
}
 
namespace LHAPDF {
 
 
  /// @brief Internal storage class for PDF data point grids
  ///
  /// We use "array" to refer to the "raw" knot grid, while "grid" means a grid-based PDF.
  class KnotArray{
  public:
 
    /// How many flavours are stored in the grid
    size_t size() const { return _shape.back(); }
    /// How many flavours are stored in the grid
    size_t numflavs() const { return _shape.back(); }
 
    /// How many x knots are there
    size_t xsize() const { return _shape[0]; }
 
    /// How many q2 knots are there
    size_t q2size() const { return _shape[1]; }
 
    /// Is this container empty?
    bool empty() const { return _grid.empty(); }
 
    /// Find the largest grid index below given x, such that xknots[index] < x
    inline size_t ixbelow(double x) const { return indexbelow(x, _xs); }
 
    /// Find the largest grid index below given q2, such that q2knots[index] < q2
    inline size_t iq2below(double q2) const { return indexbelow(q2, _q2s); }
 
    /// Convenient accessor to the grid values
    double xf(int ix, int iq2, int ipid) const {
      return _grid[ix*_shape[2]*_shape[1] + iq2*_shape[2] + ipid];
    }
 
    /// @brief Convenient accessor to the polynomial coefficients
    ///
    /// @note Returns reference rather than value, to be able to read multiple adjacent coeffs at once.
    const double& coeff(int ix, int iq2, int pid, int in) const {
      return _coeffs[ix*(_shape[1])*_shape[2]*4 + iq2*_shape[2]*4 + pid*4 + in];
    }
 
    /// Accessor to the internal 'lookup table' for the pid's
    int lookUpPid(size_t id) const { return _lookup[id]; }
 
    /// Set of x knots
    double xs(size_t id) const { return _xs[id]; }
 
    /// Set of log10(x) knots
    double logxs(size_t id) const { return _logxs[id]; }
 
    /// Set of Q2 knots
    double q2s(size_t id) const { return _q2s[id]; }
 
    /// Set of log10(Q2) knots
    double logq2s(size_t id) const { return _logq2s[id]; }
 
    /// @brief Shape of the interpolation grid
    ///
    /// @todo Document/explain better!
    size_t shape(size_t id) const { return _shape[id]; }
 
    /// Check if value within the boundaries of xknots
    bool inRangeX(double x) const {
      if (x < _xs.front()) return false;
      if (x > _xs.back())  return false;
      return true;
    }
 
    /// Check if value within the boundaries of q2knots
    bool inRangeQ2(double q2) const {
      if (q2 < _q2s.front()) return false;
      if (q2 > _q2s.back())  return false;
      return true;
    }
 
    /// @todo Document and normalise naming style
    inline int get_pid(int id) const {
      // hardcoded lookup table for particle ids
      // -6,...,-1,21/0,1,...,6,22
      // if id outside of this range, search in list of ids
      if (-6 <= id && id <= 6) return _lookup[id + 6];
      else if (id == 21) return _lookup[0 + 6];
      else if (id == 22) return _lookup[13];
      else return findPidInPids(id, _pids);
    }
 
    /// @todo Document and normalise naming style
    bool has_pid(int id) const {
      return get_pid(id) != -1;
    }
 
    /// @todo Document and normalise naming style
    void initPidLookup();
 
    /// @todo Document
    void fillLogKnots();
 
 
    /// @name Const accessors to the internal data container
    /// @{
 
    const std::vector<double>& xs() const { return _xs; }
 
    const std::vector<double>& logxs() const { return _logxs; }
 
    const std::vector<double>& q2s() const { return _q2s; }
 
    const std::vector<double>& logq2s() const { return _logq2s; }
 
    /// @}
 
 
    /// @name Non const accessors for programmatic filling
    /// @{
 
    std::vector<double>& setCoeffs() { return _coeffs; }
 
    std::vector<double>& setGrid() { return _grid; }
 
    std::vector<double>& setxknots() { return _xs; }
 
    std::vector<double>& setq2knots() { return _q2s; }
 
    std::vector<size_t>& setShape() { return _shape; }
 
    std::vector<int>& setPids() { return _pids; }
 
    /// @}
 
 
  private:
 
    /// @brief Shape of the interpolation grid
    ///
    /// @todo Explain/document better!
    std::vector<size_t> _shape;
 
    /// Grid values
    std::vector<double> _grid;
 
    /// Storage for the precomputed polynomial coefficients
    std::vector<double> _coeffs;
 
    /// Order the PIDs are filled in
    std::vector<int> _pids;
    std::vector<int> _lookup;
 
    /// @name Knots
    /// @{
    std::vector<double> _xs;
    std::vector<double> _q2s;
    std::vector<double> _logxs;
    std::vector<double> _logq2s;
    /// @}
 
  };
 
 
  /// Internal storage class for alpha_s interpolation grids
  class AlphaSArray {
  public:
 
    /// @name Construction etc.
    ///@{
 
    /// Default constructor just for std::map insertability
    AlphaSArray() {}
 
    /// Constructor from Q2 knot values and alpha_s values
    AlphaSArray(const std::vector<double>& q2knots, const std::vector<double>& as)
      : _q2s(q2knots), _as(as)
    {
      _syncq2s();
    }
 
    ///@}
 
 
    /// @name Q2 stuff
    ///@{
 
    /// Q2 knot vector accessor
    const std::vector<double>& q2s() const { return _q2s; }
 
    /// log(Q2) knot vector accessor
    const std::vector<double>& logq2s() const { return _logq2s; }
 
    /// Get the index of the closest Q2 knot row <= q2
    ///
    /// If the value is >= q2_max, return i_max-1 (for polynomial spine construction)
    size_t iq2below(double q2) const {
      // Test that Q2 is in the grid range
      if (q2 < q2s().front()) throw AlphaSError("Q2 value " + to_str(q2) + " is lower than lowest-Q2 grid point at " + to_str(q2s().front()));
      if (q2 > q2s().back()) throw AlphaSError("Q2 value " + to_str(q2) + " is higher than highest-Q2 grid point at " + to_str(q2s().back()));
      /// Find the closest knot below the requested value
      size_t i = upper_bound(q2s().begin(), q2s().end(), q2) - q2s().begin();
      if (i == q2s().size()) i -= 1; // can't return the last knot index
      i -= 1; // have to step back to get the knot <= q2 behaviour
      return i;
    }
 
    /// Get the index of the closest logQ2 knot row <= logq2
    ///
    /// If the value is >= q2_max, return i_max-1 (for polynomial spine construction)
    size_t ilogq2below(double logq2) const {
      // Test that log(Q2) is in the grid range
      if (logq2 < logq2s().front()) throw GridError("logQ2 value " + to_str(logq2) + " is lower than lowest-logQ2 grid point at " + to_str(logq2s().front()));
      if (logq2 > logq2s().back()) throw GridError("logQ2 value " + to_str(logq2) + " is higher than highest-logQ2 grid point at " + to_str(logq2s().back()));
      /// Find the closest knot below the requested value
      size_t i = upper_bound(logq2s().begin(), logq2s().end(), logq2) - logq2s().begin();
      if (i == logq2s().size()) i -= 1; // can't return the last knot index
      i -= 1; // have to step back to get the knot <= q2 behaviour
      return i;
    }
 
    ///@}
 
 
    /// @name alpha_s values at Q2 points
    ///@{
 
    /// alpha_s value accessor (const)
    const std::vector<double>& alphas() const { return _as; }
    // /// alpha_s value accessor (non-const)
    // std::vector<double>& alphas() { return _as; }
    // /// alpha_s value setter
    // void setalphas(const valarray& xfs) { _as = as; }
 
    ///@}
 
 
    /// @name alpha_s derivatives vs (log)Q2, useful for interpolation
    ///@{
 
    /// Forward derivative w.r.t. logQ2
    double ddlogq_forward(size_t i) const {
      return (alphas()[i+1] - alphas()[i]) / (logq2s()[i+1] - logq2s()[i]);
    }
 
    /// Backward derivative w.r.t. logQ2
    double ddlogq_backward(size_t i) const {
      return (alphas()[i] - alphas()[i-1]) / (logq2s()[i] - logq2s()[i-1]);
    }
 
    /// Central (avg of forward and backward) derivative w.r.t. logQ2
    double ddlogq_central(size_t i) const {
      return 0.5 * (ddlogq_forward(i) + ddlogq_backward(i));
    }
 
    ///@}
 
 
  private:
 
    /// Synchronise the log(Q2) array from the Q2 one
    void _syncq2s() {
      _logq2s.resize(_q2s.size());
      for (size_t i = 0; i < _q2s.size(); ++i) _logq2s[i] = log(_q2s[i]);
    }
 
    /// List of Q2 knots
    std::vector<double> _q2s;
    /// List of log(Q2) knots
    std::vector<double> _logq2s;
    /// List of alpha_s values across the knot array
    std::vector<double> _as;
 
  };
}
#endif