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 {


  class KnotArray{
  public:

    size_t size() const { return _shape.back(); }
    size_t numflavs() const { return _shape.back(); }

    size_t xsize() const { return _shape[0]; }

    size_t q2size() const { return _shape[1]; }

    bool empty() const { return _grid.empty(); }

    inline size_t ixbelow(double x) const { return indexbelow(x, _xs); }

    inline size_t iq2below(double q2) const { return indexbelow(q2, _q2s); }

    double xf(int ix, int iq2, int ipid) const {
      return _grid[ix*_shape[2]*_shape[1] + iq2*_shape[2] + ipid];
    }

    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];
    }

    int lookUpPid(size_t id) const { return _lookup[id]; }

    double xs(size_t id) const { return _xs[id]; }

    double logxs(size_t id) const { return _logxs[id]; }

    double q2s(size_t id) const { return _q2s[id]; }

    double logq2s(size_t id) const { return _logq2s[id]; }

    size_t shape(size_t id) const { return _shape[id]; }

    bool inRangeX(double x) const {
      if (x < _xs.front()) return false;
      if (x > _xs.back())  return false;
      return true;
    }

    bool inRangeQ2(double q2) const {
      if (q2 < _q2s.front()) return false;
      if (q2 > _q2s.back())  return false;
      return true;
    }

    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);
    }

    bool has_pid(int id) const {
      return get_pid(id) != -1;
    }

    void initPidLookup();

    void fillLogKnots();



    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; }




    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:

    std::vector<size_t> _shape;

    std::vector<double> _grid;

    std::vector<double> _coeffs;

    std::vector<int> _pids;
    std::vector<int> _lookup;

    std::vector<double> _xs;
    std::vector<double> _q2s;
    std::vector<double> _logxs;
    std::vector<double> _logq2s;

  };


  class AlphaSArray {
  public:


    AlphaSArray() {}

    AlphaSArray(const std::vector<double>& q2knots, const std::vector<double>& as)
      : _q2s(q2knots), _as(as)
    {
      _syncq2s();
    }




    const std::vector<double>& q2s() const { return _q2s; }

    const std::vector<double>& logq2s() const { return _logq2s; }

    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()));
      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;
    }

    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()));
      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;
    }




    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; }




    double ddlogq_forward(size_t i) const {
      return (alphas()[i+1] - alphas()[i]) / (logq2s()[i+1] - logq2s()[i]);
    }

    double ddlogq_backward(size_t i) const {
      return (alphas()[i] - alphas()[i-1]) / (logq2s()[i] - logq2s()[i-1]);
    }

    double ddlogq_central(size_t i) const {
      return 0.5 * (ddlogq_forward(i) + ddlogq_backward(i));
    }



  private:

    void _syncq2s() {
      _logq2s.resize(_q2s.size());
      for (size_t i = 0; i < _q2s.size(); ++i) _logq2s[i] = log(_q2s[i]);
    }

    std::vector<double> _q2s;
    std::vector<double> _logq2s;
    std::vector<double> _as;

  };
}
#endif