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:
  GPUd() static constexpr uint32_t getNParametersAllTerms(const uint32_t degree, const uint32_t dim) { return (degree == 0) ? binomialCoeff(dim - 1, 0) : binomialCoeff(dim - 1 + degree, degree) + getNParametersAllTerms(degree - 1, dim); }
 
  GPUd() static constexpr uint32_t getNParametersInteractionOnly(const uint32_t degree, const uint32_t dim) { return (degree == 0) ? binomialCoeff(dim - 1, 0) : binomialCoeff(dim, degree) + getNParametersInteractionOnly(degree - 1, dim); }
 
  GPUd() static constexpr uint32_t getNParameters(const uint32_t degree, const uint32_t dim, const bool interactionOnly)
  {
    if (interactionOnly) {
      return getNParametersInteractionOnly(degree, dim);
    } else {
      return getNParametersAllTerms(degree, dim);
    }
  }
 
 private:
  GPUd() static constexpr uint32_t factorial(const uint32_t n) { return (n == 0) || (n == 1) ? 1 : n * factorial(n - 1); }
 
  GPUd() static constexpr uint32_t binomialCoeff(const uint32_t n, const 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); }
 
  GPUd() static constexpr uint32_t mod10(const uint32_t a, const uint32_t b) { 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);
 
  GPUd() static constexpr uint32_t getNewPos(const uint32_t degreePol, const uint32_t pos, const uint32_t digitPos);
 
  template <uint32_t DegreePol>
  GPUd() static constexpr float prodTerm(const float x[], const uint32_t pos);
 
  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>
GPUd() constexpr uint32_t MultivariatePolynomialHelper<Dim, Degree, InteractionOnly>::getNewPos(const uint32_t degreePol, const uint32_t pos, const uint32_t digitPos)
{
  if (degreePol > digitPos) {
    // check if digit of current position is at is max position
    if (mod10(pos, pow10(digitPos)) == Dim) {
      // increase digit of left position
      const uint32_t leftDigit = digitPos + 1;
      const uint32_t posTmp = pos + pow10(leftDigit);
      const uint32_t refDigit = mod10(posTmp, pow10(digitPos + 1));
 
      // resetting digits to the right if digit exceeds number of dimensions
      const uint32_t posReset = resetIndices(degreePol, posTmp, leftDigit - 1, degreePol - digitPos, refDigit);
 
      // check next digit
      return getNewPos(degreePol, posReset, digitPos + 1);
    }
    return getNewPos(degreePol, pos, digitPos + 1);
  }
  return pos;
}
 
template <uint32_t Dim, uint32_t Degree, bool InteractionOnly>
template <uint32_t DegreePol>
GPUd() constexpr float MultivariatePolynomialHelper<Dim, Degree, InteractionOnly>::prodTerm(const float x[], const uint32_t pos)
{
  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>(x, pos);
  }
  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);
    }
 
    // sum up the term for corrent term and set posotion for next combination
    return par[Index] * prodTerm<DegreePol>(x, posNew) + 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