Chebyshev3DCalc.h

Go to the documentation of this file.

// Copyright 2019-2020 CERN and copyright holders of ALICE O2.
// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders.
// All rights not expressly granted are reserved.
//
// This software is distributed under the terms of the GNU General Public
// License v3 (GPL Version 3), copied verbatim in the file "COPYING".
//
// In applying this license CERN does not waive the privileges and immunities
// granted to it by virtue of its status as an Intergovernmental Organization
// or submit itself to any jurisdiction.
 
 
#ifndef ALICEO2_MATHUTILS_CHEBYSHEV3DCALC_H_
#define ALICEO2_MATHUTILS_CHEBYSHEV3DCALC_H_
 
#include <TNamed.h> // for TNamed
#include <cstdio>   // for FILE, stdout
#include <cmath>    // for std::fma
#include <algorithm>
#include <vector>
#include "Rtypes.h" // for Float_t, UShort_t, Int_t, Double_t, etc
 
class TString;
 
// To decrease the compilable code size comment this define. This will exclude the routines
// used for the calculation and saving of the coefficients.
#define _INC_CREATION_Chebyshev3D_
 
// When _BRING_TO_BOUNDARY_ is defined, the point outside of the fitted folume is assumed to be on the surface
// #define _BRING_TO_BOUNDARY_
 
namespace o2
{
namespace math_utils
{
class Chebyshev3DCalc : public TNamed
{
 
 public:
  Chebyshev3DCalc();
 
  Chebyshev3DCalc(const Chebyshev3DCalc& src);
 
  Chebyshev3DCalc(FILE* stream);
 
  ~Chebyshev3DCalc() override
  {
    Clear();
  }
 
  Chebyshev3DCalc& operator=(const Chebyshev3DCalc& rhs);
 
  void Print(const Option_t* opt = "") const override;
 
  void loadData(FILE* stream);
 
  Float_t evaluateDerivative(int dim, const Float_t* par) const;
 
  Float_t evaluateDerivative2(int dim1, int dim2, const Float_t* par) const;
 
#ifdef _INC_CREATION_Chebyshev3D_
 
  void saveData(const char* outfile, Bool_t append = kFALSE) const;
 
  // Writes coefficients data to existing output stream
  // Note: mNumberOfColumns, mNumberOfElementsBound2D and mColumnAtRowBeginning are not stored, will be computed on fly
  // during the loading of this file
  void saveData(FILE* stream = stdout) const;
 
#endif
 
  void initializeRows(int nr);
 
  void initializeColumns(int nc);
 
  Int_t getNumberOfCoefficients() const
  {
    return mNumberOfCoefficients;
  }
 
  Int_t getNumberOfColumns() const
  {
    return (Int_t)mNumberOfColumns;
  }
 
  Int_t getNumberOfRows() const
  {
    return (Int_t)mNumberOfRows;
  }
 
  Int_t getNumberOfElementsBound2D() const
  {
    return (Int_t)mNumberOfElementsBound2D;
  }
 
  Int_t getMaxColumnsAtRow() const;
 
  UShort_t* getNumberOfColumnsAtRow() const
  {
    return mNumberOfColumnsAtRow;
  }
 
  UShort_t* getColAtRowBg() const
  {
    return mColumnAtRowBeginning;
  }
 
  Float_t getPrecision() const
  {
    return mPrecision;
  }
 
  void setPrecision(Float_t prc = 1e-6)
  {
    mPrecision = prc;
  }
 
  void initializeElementBound2D(int ne);
 
  UShort_t* getCoefficientBound2D0() const
  {
    return mCoefficientBound2D0;
  }
 
  UShort_t* getCoefficientBound2D1() const
  {
    return mCoefficientBound2D1;
  }
 
  void Clear(const Option_t* option = "") override;
 
  static Float_t chebyshevEvaluation1D(Float_t x, const Float_t* array, int ncf);
 
  static Float_t chebyshevEvaluation1Derivative(Float_t x, const Float_t* array, int ncf);
 
  static Float_t chebyshevEvaluation1Derivative2(Float_t x, const Float_t* array, int ncf);
 
  void initializeCoefficients(int nc);
 
  Float_t* getCoefficients() const
  {
    return mCoefficients;
  }
 
  static void readLine(TString& str, FILE* stream);
 
  Float_t Eval(const Float_t* par) const;
 
  Double_t Eval(const Double_t* par) const;
 
 private:
  class EvalScratch
  {
   public:
    EvalScratch(int nrows, int ncols)
    {
      if (nrows <= kStackSize && ncols <= kStackSize) {
        mRows = mStackRows;
        mCols = mStackCols;
      } else {
        thread_local std::vector<Float_t> rows, cols;
        rows.resize(std::max<size_t>(nrows, rows.size()));
        cols.resize(std::max<size_t>(ncols, cols.size()));
        mRows = rows.data();
        mCols = cols.data();
      }
    }
    EvalScratch(const EvalScratch&) = delete;
    EvalScratch& operator=(const EvalScratch&) = delete;
    Float_t* rows() { return mRows; }
    Float_t* cols() { return mCols; }
 
   private:
    static constexpr int kStackSize = 64;
    Float_t mStackRows[kStackSize];
    Float_t mStackCols[kStackSize];
    Float_t* mRows;
    Float_t* mCols;
  };
 
  Int_t mNumberOfCoefficients;    
  Int_t mNumberOfRows;            
  Int_t mNumberOfColumns;         
  Int_t mNumberOfElementsBound2D; 
  Float_t mPrecision;             
  UShort_t*
    mNumberOfColumnsAtRow;         //[mNumberOfRows] number of sighificant columns (2nd dim) at each row of 3D coefs matrix
  UShort_t* mColumnAtRowBeginning; //[mNumberOfRows] beginning of significant columns (2nd dim) for row in the 2D
  // boundary matrix
  UShort_t* mCoefficientBound2D0; //[mNumberOfElementsBound2D] 2D matrix defining the boundary of significance for 3D
  // coeffs.matrix
  //(Ncoefs for col/row)
  UShort_t* mCoefficientBound2D1; //[mNumberOfElementsBound2D] 2D matrix defining the start beginning of significant
  // coeffs for col/row
  Float_t* mCoefficients; //[mNumberOfCoefficients] array of Chebyshev coefficients
 
  Float_t* mTemporaryCoefficients2D; //[mNumberOfColumns] no longer used for evaluation
  Float_t* mTemporaryCoefficients1D; //[mNumberOfRows] no longer used for evaluation
 
  ClassDefOverride(o2::math_utils::Chebyshev3DCalc,
                   2) // Class for interpolation of 3D->1 function by Chebyshev parametrization
};
 
inline Float_t Chebyshev3DCalc::chebyshevEvaluation1D(Float_t x, const Float_t* array, int ncf)
{
  if (ncf <= 0) {
    return 0;
  }
 
  Float_t b0, b1, b2, x2 = x + x;
  b0 = array[--ncf];
  b1 = b2 = 0;
 
  for (int i = ncf; i--;) {
    b2 = b1;
    b1 = b0;
    // Clenshaw recurrence, grouped as fma(x2, b1, array[i] - b2). Mathematically
    // identical to `array[i] + x2 * b1 - b2`, but `array[i] - b2` does not depend
    // on the just-updated b1, so the loop-carried chain collapses to a single FMA
    // latency instead of a dependent multiply+add+subtract. This kernel dominates
    // magnetic-field evaluation in (e.g.) muon track extrapolation.
    b0 = std::fma(x2, b1, array[i] - b2);
  }
  return std::fma(-x, b1, b0);
}
 
inline Float_t Chebyshev3DCalc::Eval(const Float_t* par) const
{
  EvalScratch scratch(mNumberOfRows, mNumberOfColumns);
  Float_t* coefs2D = scratch.cols();
  Float_t* coefs1D = scratch.rows();
  for (int id0 = mNumberOfRows; id0--;) {
    int nCLoc = mNumberOfColumnsAtRow[id0]; // number of significant coefs on this row
    int col0 = mColumnAtRowBeginning[id0];  // beginning of local column in the 2D boundary matrix
    for (int id1 = nCLoc; id1--;) {
      int id = id1 + col0;
      coefs2D[id1] = chebyshevEvaluation1D(par[2], mCoefficients + mCoefficientBound2D1[id], mCoefficientBound2D0[id]);
    }
    coefs1D[id0] = chebyshevEvaluation1D(par[1], coefs2D, nCLoc);
  }
  return chebyshevEvaluation1D(par[0], coefs1D, mNumberOfRows);
}
 
inline Double_t Chebyshev3DCalc::Eval(const Double_t* par) const
{
  EvalScratch scratch(mNumberOfRows, mNumberOfColumns);
  Float_t* coefs2D = scratch.cols();
  Float_t* coefs1D = scratch.rows();
  for (int id0 = mNumberOfRows; id0--;) {
    int nCLoc = mNumberOfColumnsAtRow[id0]; // number of significant coefs on this row
    int col0 = mColumnAtRowBeginning[id0];  // beginning of local column in the 2D boundary matrix
    for (int id1 = nCLoc; id1--;) {
      int id = id1 + col0;
      coefs2D[id1] = chebyshevEvaluation1D(par[2], mCoefficients + mCoefficientBound2D1[id], mCoefficientBound2D0[id]);
    }
    coefs1D[id0] = chebyshevEvaluation1D(par[1], coefs2D, nCLoc);
  }
  return chebyshevEvaluation1D(par[0], coefs1D, mNumberOfRows);
}
} // namespace math_utils
} // namespace o2
 
#endif