// -*- C++ -*- // // This file is part of LHAPDF // Copyright (C) 2012-2021 The LHAPDF collaboration (see AUTHORS for details) // #include "LHAPDF/LogBicubicInterpolator.h" #include namespace LHAPDF { namespace { // Unnamed namespace struct shared_data{ // Pre-calculate parameters double logx, logq2, dlogx_1, dlogq_0, dlogq_1, dlogq_2, tlogq; double tlogx; // bools to find out if at grid edges bool q2_lower, q2_upper; }; shared_data fill(const KnotArray& grid, double x, double q2, size_t ix, size_t iq2){ shared_data shared; shared.logx = std::log(x); shared.logq2 = std::log(q2); // only required for fallback mode // Pre-calculate parameters shared.dlogx_1 = grid.logxs(ix+1) - grid.logxs(ix); shared.tlogx = (shared.logx - grid.logxs(ix)) / shared.dlogx_1; shared.dlogq_0 = 1./(grid.logq2s(iq2) - grid.logq2s(iq2-1)); // only ever need the 1/values of the differences shared.dlogq_1 = grid.logq2s(iq2+1) - grid.logq2s(iq2); shared.dlogq_2 = 1./(grid.logq2s(iq2+2) - grid.logq2s(iq2+1)); shared.tlogq = (shared.logq2 - grid.logq2s(iq2)) / shared.dlogq_1; shared.q2_lower = ( (iq2 == 0) || (grid.q2s(iq2) == grid.q2s(iq2-1))); shared.q2_upper = ( (iq2+1 == grid.q2size() - 1 ) || (grid.q2s(iq2+1) == grid.q2s(iq2+2)) ); return shared; } /// One-dimensional linear interpolation for y(x) /// @todo Expose for re-use inline double _interpolateLinear(double x, double xl, double xh, double yl, double yh) { assert(x >= xl); assert(xh >= x); return yl + (x - xl) / (xh - xl) * (yh - yl); } /// One-dimensional cubic interpolation /// /// @arg t is the fractional distance of the evaluation x into the dx /// interval. @arg vl and @arg vh are the function values at the low and /// high edges of the interval. @arg vl and @arg vh are linearly /// extrapolated value changes from the product of dx and the discrete low- /// and high-edge derivative estimates. /// /// @note See Numerical Recipes 3.6: http://www.it.uom.gr/teaching/linearalgebra/NumericalRecipiesInC/c3-6.pdf /// /// @todo Expose for re-use inline double _interpolateCubic(double t, double vl, double vdl, double vh, double vdh) { // Pre-calculate powers of t const double t2 = t*t; const double t3 = t*t2; // Calculate polynomial (grouped by input param rather than powers of t for efficiency) const double p0 = (2*t3 - 3*t2 + 1)*vl; const double m0 = (t3 - 2*t2 + t)*vdl; const double p1 = (-2*t3 + 3*t2)*vh; const double m1 = (t3 - t2)*vdh; return p0 + m0 + p1 + m1; } /// Cubic interpolation using a passed array of coefficients /// /// @todo Expose for re-use inline double _interpolateCubic(double t, const double* coeffs) { const double x = t; const double x2 = x*x; const double x3 = x2*x; return coeffs[0]*x3 + coeffs[1]*x2 + coeffs[2]*x + coeffs[3]; } double _interpolate(const KnotArray& grid, size_t ix, size_t iq2, int id, shared_data& _share) { double vl = _interpolateCubic(_share.tlogx, &grid.coeff(ix,iq2,id,0)); double vh = _interpolateCubic(_share.tlogx, &grid.coeff(ix,iq2+1,id,0)); // Derivatives in Q2 /// @todo changes derivative approximation to take knot distance into account double vdl, vdh; if (_share.q2_lower) { // Forward difference for lower q vdl = (vh - vl); // Central difference for higher q double vhh = _interpolateCubic(_share.tlogx, &grid.coeff(ix,iq2+2,id,0)); vdh = (vdl + (vhh - vh) * _share.dlogq_1 * _share.dlogq_2) * 0.5; } else if (_share.q2_upper) { // Backward difference for higher q vdh = (vh - vl); // Central difference for lower q double vll = _interpolateCubic(_share.tlogx, &grid.coeff(ix,iq2-1,id,0)); vdl = (vdh + (vl - vll) * _share.dlogq_1 * _share.dlogq_0) * 0.5; } else { // Central difference for both q double vll = _interpolateCubic(_share.tlogx, &grid.coeff(ix,iq2-1,id,0)); // replace with better derivative estimate? vdl = ( (vh - vl) + (vl - vll)*_share.dlogq_1 * _share.dlogq_0 ) * 0.5; double vhh = _interpolateCubic(_share.tlogx, &grid.coeff(ix,iq2+2,id,0)); vdh = ( (vh - vl) + (vhh - vh)*_share.dlogq_1 * _share.dlogq_2 ) * 0.5; } return _interpolateCubic(_share.tlogq, vl, vdl, vh, vdh); } double _interpolateFallback(const KnotArray& grid, size_t ix, size_t iq2, int id, shared_data& _share) { // First interpolate in x const double logx0 = grid.logxs(ix); const double logx1 = grid.logxs(ix+1); const double f_ql = _interpolateLinear(_share.logx, logx0, logx1, grid.xf(ix, iq2, id), grid.xf(ix+1, iq2, id)); const double f_qh = _interpolateLinear(_share.logx, logx0, logx1, grid.xf(ix, iq2+1, id), grid.xf(ix+1, iq2+1, id)); // Then interpolate in Q2, using the x-ipol results as anchor points return _interpolateLinear(_share.logq2, grid.logq2s(iq2), grid.logq2s(iq2+1), f_ql, f_qh); } void _checkGridSize(const KnotArray& grid, const size_t ix, const size_t iq2) { // Raise an error if there are too few knots even for a linear fall-back const size_t nxknots = grid.xsize(); const size_t nq2knots = grid.q2size(); /// @todo MK: do you really need different number of knots in the directions? /// Probably should be <2 for both methods, and fall back to linear in both cases. if (nxknots < 4) throw GridError("PDF subgrids are required to have at least 4 x-knots for use with LogBicubicInterpolator"); if (nq2knots < 2) throw GridError("PDF subgrids are required to have at least 2 Q-knots for use with LogBicubicInterpolator"); // Check x and q index ranges -- we always need i and i+1 indices to be valid const size_t ixmax = nxknots - 1; const size_t iq2max = nq2knots - 1; if (ix+1 > ixmax) // also true if ix is off the end throw GridError("Attempting to access an x-knot index past the end of the array, in linear fallback mode"); if (iq2+1 > iq2max) // also true if iq2 is off the end throw GridError("Attempting to access an Q-knot index past the end of the array, in linear fallback mode"); } } double LogBicubicInterpolator::_interpolateXQ2(const KnotArray& grid, double x, size_t ix, double q2, size_t iq2, int id) const { _checkGridSize(grid, ix, iq2); shared_data shared = fill(grid, x, q2, ix, iq2); // if it is an upper and a lower edge, there can only be two notes // use lineare fallback in that case, but have default be the cubic case if (!shared.q2_lower || !shared.q2_upper) return _interpolate(grid, ix, iq2, id, shared); // Fallback mode return _interpolateFallback(grid, ix, iq2, id, shared); } /// Interpolate a whole set of standard parton IDs at once, minimising recomputation void LogBicubicInterpolator::_interpolateXQ2(const KnotArray& grid, double x, size_t ix, double q2, size_t iq2, std::vector& ret) const { _checkGridSize(grid, ix, iq2); shared_data shared = fill(grid, x, q2, ix, iq2); if (!shared.q2_lower || !shared.q2_upper) { for (int pid = -6; pid <= 6; ++pid) { int id = grid.lookUpPid(pid + 6); ret[pid + 6] = (id == -1) ? 0 : _interpolate(grid, ix, iq2, id, shared); } } else { for (int pid = -6; pid <= 6; ++pid) { const int id = grid.lookUpPid(pid + 6); ret[pid + 6] = (id == -1) ? 0 : _interpolateFallback(grid, ix, iq2, id, shared); } } } }