MultivariatePolynomialHelper.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_TPC_MULTIVARIATEPOLYNOMIALHELPER
#define ALICEO2_TPC_MULTIVARIATEPOLYNOMIALHELPER
 
#include "GPUCommonDef.h"
 
#if !defined(GPUCA_GPUCODE)
#include <vector>
#include <string>
#include <cassert>
#if !defined(GPUCA_NO_FMT)
#include <fmt/format.h>
#endif
#endif
 
class TLinearFitter;
 
namespace o2::gpu
{
 
#if !defined(GPUCA_GPUCODE)
 
struct MultivariatePolynomialContainer {
 
  MultivariatePolynomialContainer(const uint32_t dim, const uint32_t degree, const uint32_t nParameters, const float params[/* nParameters*/], const bool interactionOnly) : mDim{dim}, mDegree{degree}, mParams{params, params + nParameters}, mInteractionOnly{interactionOnly} {};
 
  MultivariatePolynomialContainer() = default;
 
  const uint32_t mDim{};              
  const uint32_t mDegree{};           
  const std::vector<float> mParams{}; 
  const bool mInteractionOnly{};      
};
#endif
 
class MultivariatePolynomialParametersHelper
{
 public:
  template <uint32_t Degree, uint32_t Dim>
  GPUd() static constexpr uint32_t getNParametersAllTerms()
  {
    if constexpr (Degree == 0) {
      return binomialCoeff<Dim - 1, 0>();
    } else {
      return binomialCoeff<Dim - 1 + Degree, Degree>() + getNParametersAllTerms<Degree - 1, Dim>();
    }
  }
 
  GPUd() static constexpr uint32_t getNParametersAllTerms(uint32_t degree, uint32_t dim)
  {
    if (degree == 0) {
      return binomialCoeff(dim - 1, 0);
    } else {
      return binomialCoeff(dim - 1 + degree, degree) + getNParametersAllTerms(degree - 1, dim);
    }
  }
 
  template <uint32_t Degree, uint32_t Dim>
  GPUd() static constexpr uint32_t getNParametersInteractionOnly()
  {
    if constexpr (Degree == 0) {
      return binomialCoeff<Dim - 1, 0>();
    } else {
      return binomialCoeff<Dim, Degree>() + getNParametersInteractionOnly<Degree - 1, Dim>();
    }
  }
 
  GPUd() static constexpr uint32_t getNParametersInteractionOnly(uint32_t degree, uint32_t dim)
  {
    if (degree == 0) {
      return binomialCoeff(dim - 1, 0);
    } else {
      return binomialCoeff(dim, degree) + getNParametersInteractionOnly(degree - 1, dim);
    }
  }
 
  template <uint32_t Degree, uint32_t Dim, bool InteractionOnly>
  GPUd() static constexpr uint32_t getNParameters()
  {
    if constexpr (InteractionOnly) {
      return getNParametersInteractionOnly<Degree, Dim>();
    } else {
      return getNParametersAllTerms<Degree, Dim>();
    }
  }
 
  GPUd() static constexpr uint32_t getNParameters(uint32_t degree, uint32_t dim, bool interactionOnly)
  {
    if (interactionOnly) {
      return getNParametersInteractionOnly(degree, dim);
    } else {
      return getNParametersAllTerms(degree, dim);
    }
  }
 
 private:
  template <uint32_t N>
  GPUd() static constexpr uint32_t factorial()
  {
    if constexpr (N == 0 || N == 1) {
      return 1;
    } else {
      return N * factorial<N - 1>();
    }
  }
 
  GPUd() static constexpr uint32_t factorial(uint32_t n) { return n == 0 || n == 1 ? 1 : n * factorial(n - 1); }
 
  template <uint32_t N, uint32_t K>
  GPUd() static constexpr uint32_t binomialCoeff()
  {
    return factorial<N>() / (factorial<K>() * factorial<N - K>());
  }
 
  GPUd() static constexpr uint32_t binomialCoeff(uint32_t n, uint32_t k)
  {
    return factorial(n) / (factorial(k) * factorial(n - k));
  }
};
 
template <uint32_t Dim, uint32_t Degree, bool InteractionOnly>
class MultivariatePolynomialHelper : public MultivariatePolynomialParametersHelper
{
  static constexpr uint16_t FMaxdim = 10;   
  static constexpr uint16_t FMaxdegree = 9; 
 
#if !defined(GPUCA_GPUCODE)
  static_assert(Dim <= MultivariatePolynomialHelper<Dim, Degree, InteractionOnly>::FMaxdim && Degree <= MultivariatePolynomialHelper<Dim, Degree, InteractionOnly>::FMaxdegree, "Max. number of dimensions or degrees exceeded!");
#endif
 
 public:
  GPUd() static constexpr float evalPol(GPUgeneric() const float par[/*number of parameters*/], const float x[/*number of dimensions*/])
  {
    return par[0] + loopDegrees<1>(par, x);
  }
 
  GPUd() static constexpr uint32_t getDim() { return Dim; }
 
  GPUd() static constexpr uint32_t getDegree() { return Degree; }
 
  GPUd() static constexpr bool isInteractionOnly() { return InteractionOnly; }
 
 private:
  GPUd() static constexpr uint32_t pow10(const uint32_t n) { return n == 0 ? 1 : 10 * pow10(n - 1); }
 
  template <uint32_t N>
  GPUd() static constexpr uint32_t pow10()
  {
    if constexpr (N == 0) {
      return 1;
    } else {
      return 10 * pow10<N - 1>();
    }
  }
 
  GPUd() static constexpr uint32_t mod10(const uint32_t a, const uint32_t b) { return (a / b) % 10; }
 
  template <uint32_t A, uint32_t B>
  GPUd() static constexpr uint32_t mod10()
  {
    return (A / B) % 10;
  }
 
  GPUd() static constexpr uint32_t resetIndices(const uint32_t degreePol, const uint32_t pos, const uint32_t leftDigit, const uint32_t iter, const uint32_t refDigit);
 
  template <uint32_t DegreePol, uint32_t Pos, uint32_t DigitPos>
  GPUd() static constexpr uint32_t getNewPos();
 
  template <uint32_t DegreePol, uint32_t Pos>
  GPUd() static constexpr float prodTerm(const float x[]);
 
  template <uint32_t DegreePol, uint32_t posNew>
  static constexpr bool checkInteraction();
 
  template <uint32_t DegreePol, uint32_t Pos, uint32_t Index>
  GPUd() static constexpr float sumTerms(GPUgeneric() const float par[], const float x[]);
 
  template <uint32_t DegreePol>
  GPUd() static constexpr float loopDegrees(GPUgeneric() const float par[], const float x[]);
};
 
#if !defined(GPUCA_GPUCODE) // disable specialized class for GPU, as it should be only used for the fitting!
template <>
class MultivariatePolynomialHelper<0, 0, false> : public MultivariatePolynomialParametersHelper
{
 public:
  MultivariatePolynomialHelper(const uint32_t nDim, const uint32_t degree, const bool interactionOnly) : mDim{nDim}, mDegree{degree}, mInteractionOnly{interactionOnly} { assert(mDegree <= FMaxdegree); };
 
  MultivariatePolynomialHelper() = default;
 
  ~MultivariatePolynomialHelper() = default;
 
  void print() const;
 
  std::string getFormula() const;
 
  std::string getTLinearFitterFormula() const;
 
  std::vector<std::string> getTerms() const;
 
  TLinearFitter getTLinearFitter() const;
 
  static std::vector<float> fit(TLinearFitter& fitter, std::vector<double>& x, std::vector<double>& y, std::vector<double>& error, const bool clearPoints);
 
  std::vector<float> fit(std::vector<double>& x, std::vector<double>& y, std::vector<double>& error, const bool clearPoints) const;
 
  float evalPol(const float par[/*number of parameters*/], const float x[/*number of dimensions*/]) const
  {
    return evalPol(par, x, mDegree, mDim, mInteractionOnly);
  }
 
  float evalPol(const float par[], const float x[], const uint32_t degree, const uint32_t dim, const bool interactionOnly) const;
 
  uint32_t getDim() const { return mDim; }
 
  uint32_t getDegree() const { return mDegree; }
 
  bool isInteractionOnly() const { return mInteractionOnly; }
 
 protected:
  uint32_t mDim{};         
  uint32_t mDegree{};      
  bool mInteractionOnly{}; 
 
 private:
  static constexpr uint16_t FMaxdegree = 9; 
 
  template <class Type>
  Type combination_with_repetiton(const uint32_t degree, const uint32_t dim, const float par[], int32_t& indexPar, const float x[], const bool interactionOnly) const;
};
#endif
 
//=================================================================================
//============================ inline implementations =============================
//=================================================================================
 
template <uint32_t Dim, uint32_t Degree, bool InteractionOnly>
GPUd() constexpr uint32_t MultivariatePolynomialHelper<Dim, Degree, InteractionOnly>::resetIndices(const uint32_t degreePol, const uint32_t pos, const uint32_t leftDigit, const uint32_t iter, const uint32_t refDigit)
{
  if (iter <= degreePol) {
    const int32_t powTmp = pow10(leftDigit);
    const int32_t rightDigit = mod10(pos, powTmp);
    const int32_t posTmp = pos - (rightDigit - refDigit) * powTmp;
    return resetIndices(degreePol, posTmp, leftDigit - 1, iter + 1, refDigit);
  }
  return pos;
}
 
template <uint32_t Dim, uint32_t Degree, bool InteractionOnly>
template <uint32_t DegreePol, uint32_t Pos, uint32_t DigitPos>
GPUd() constexpr uint32_t MultivariatePolynomialHelper<Dim, Degree, InteractionOnly>::getNewPos()
{
  if constexpr (DegreePol > DigitPos) {
    // check if digit of current position is at is max position
    if constexpr (mod10<Pos, pow10<DigitPos>()>() == Dim) {
      // increase digit of left position
      constexpr uint32_t LeftDigit = DigitPos + 1;
      constexpr uint32_t PowLeftDigit = pow10<LeftDigit>();
      constexpr uint32_t PosTmp = Pos + PowLeftDigit;
      constexpr uint32_t RefDigit = mod10<PosTmp, PowLeftDigit>();
 
      // resetting digits to the right if digit exceeds number of dimensions
      constexpr uint32_t PosReset = resetIndices(DegreePol, PosTmp, LeftDigit - 1, DegreePol - DigitPos, RefDigit);
 
      // check next digit
      return getNewPos<DegreePol, PosReset, DigitPos + 1>();
    } else {
      return getNewPos<DegreePol, Pos, DigitPos + 1>();
    }
  } else {
    return Pos;
  }
}
 
template <uint32_t Dim, uint32_t Degree, bool InteractionOnly>
template <uint32_t DegreePol, uint32_t Pos>
GPUd() constexpr float MultivariatePolynomialHelper<Dim, Degree, InteractionOnly>::prodTerm(const float x[])
{
  if constexpr (DegreePol > 0) {
    // extract index of the dimension which is decoded in the digit
    const uint32_t index = mod10<Pos, pow10<DegreePol - 1>()>();
    return x[index] * prodTerm<DegreePol - 1, Pos>(x);
  }
  return 1;
}
 
template <uint32_t Dim, uint32_t Degree, bool InteractionOnly>
template <uint32_t DegreePol, uint32_t posNew>
constexpr bool MultivariatePolynomialHelper<Dim, Degree, InteractionOnly>::checkInteraction()
{
  if constexpr (DegreePol > 1) {
    constexpr bool isInteraction = mod10<posNew, pow10<DegreePol - 1>()>() == mod10<posNew, pow10<DegreePol - 2>()>();
    if constexpr (isInteraction) {
      return true;
    }
    return checkInteraction<DegreePol - 1, posNew>();
  }
  return false;
}
 
template <uint32_t Dim, uint32_t Degree, bool InteractionOnly>
template <uint32_t DegreePol, uint32_t Pos, uint32_t Index>
GPUd() constexpr float MultivariatePolynomialHelper<Dim, Degree, InteractionOnly>::sumTerms(GPUgeneric() const float par[], const float x[])
{
  // checking if the current position is reasonable e.g. if the max dimension is x[4]: for Pos=15 -> x[1]*x[5] the position is set to 22 -> x[2]*x[2]
  constexpr uint32_t PosNew = getNewPos<DegreePol, Pos, 0>();
  if constexpr (mod10<PosNew, pow10<DegreePol>()>() != 1) {
 
    // check if all digits in posNew are unequal: For interaction_only terms with x[Dim]*x[Dim]... etc. can be skipped
    if constexpr (InteractionOnly && checkInteraction<DegreePol, PosNew>()) {
      return sumTerms<DegreePol, PosNew + 1, Index>(par, x);
    } else {
      // sum up the term for corrent term and set posotion for next combination
      return par[Index] * prodTerm<DegreePol, PosNew>(x) + sumTerms<DegreePol, PosNew + 1, Index + 1>(par, x);
    }
  }
  return 0;
}
 
template <uint32_t Dim, uint32_t Degree, bool InteractionOnly>
template <uint32_t DegreePol>
GPUd() constexpr float MultivariatePolynomialHelper<Dim, Degree, InteractionOnly>::loopDegrees(GPUgeneric() const float par[], const float x[])
{
  if constexpr (DegreePol <= Degree) {
    constexpr uint32_t index{getNParameters<DegreePol - 1, Dim, InteractionOnly>()}; // offset of the index for accessing the parameters
    return sumTerms<DegreePol, 0, index>(par, x) + loopDegrees<DegreePol + 1>(par, x);
  }
  return 0;
}
 
} // namespace o2::gpu
 
#endif