NDPiecewisePolynomials.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_NDPIECEWISEPOLYNOMIALS
#define ALICEO2_TPC_NDPIECEWISEPOLYNOMIALS
 
#include "GPUCommonLogger.h"
#include "FlatObject.h"
#include "MultivariatePolynomialHelper.h"
#include "GPUCommonMath.h"
 
#if !defined(GPUCA_GPUCODE) && !defined(GPUCA_STANDALONE)
#include <vector>
#include <functional>
#endif
 
class TFile;
 
namespace o2::gpu
{
 
#if !defined(GPUCA_GPUCODE)
struct NDPiecewisePolynomialContainer {
 
  NDPiecewisePolynomialContainer(const uint32_t dim, const uint32_t degree, const uint32_t nParameters, const float params[/* nParameters*/], const bool interactionOnly, const float xMin[/* Dim*/], const float xMax[/* Dim*/], const uint32_t nVertices[/* Dim*/]) : mDim{dim}, mDegree{degree}, mParams{params, params + nParameters}, mInteractionOnly{interactionOnly}, mMin{xMin, xMin + dim}, mMax{xMax, xMax + dim}, mN{nVertices, nVertices + dim} {};
 
  NDPiecewisePolynomialContainer() = default;
 
  const uint32_t mDim{};              
  const uint32_t mDegree{};           
  const std::vector<float> mParams{}; 
  const bool mInteractionOnly{};      
  const std::vector<float> mMin{};    
  const std::vector<float> mMax{};    
  const std::vector<uint32_t> mN{};   
};
#endif
 
template <uint32_t Dim, uint32_t Degree, bool InteractionOnly>
class NDPiecewisePolynomials : public FlatObject
{
 public:
#if !defined(GPUCA_GPUCODE) && !defined(GPUCA_STANDALONE)
  NDPiecewisePolynomials(const float min[/* Dim */], const float max[/* Dim */], const uint32_t n[/* Dim */]) { init(min, max, n); }
  NDPiecewisePolynomials(const char* fileName, const char* name)
  {
    loadFromFile(fileName, name);
  };
#endif // !defined(GPUCA_GPUCODE) && !defined(GPUCA_STANDALONE)
  NDPiecewisePolynomials() = default;
 
  ~NDPiecewisePolynomials() = default;
 
  NDPiecewisePolynomials(const NDPiecewisePolynomials& obj) { cloneFromObject(obj, nullptr); }
 
#if !defined(GPUCA_GPUCODE)
  void cloneFromObject(const NDPiecewisePolynomials& obj, char* newFlatBufferPtr);
 
  void moveBufferTo(char* newBufferPtr);
#endif // !defined(GPUCA_GPUCODE)
 
  void destroy();
 
  void setActualBufferAddress(char* actualFlatBufferPtr);
 
  void setFutureBufferAddress(char* futureFlatBufferPtr);
 
  GPUd() float eval(float x[/* Dim */]) const
  {
    int32_t index[Dim];
    setIndex<Dim - 1>(x, index);
    clamp<Dim - 1>(x, index);
    return MultivariatePolynomialHelper<Dim, Degree, InteractionOnly>::evalPol(getParameters(index), x);
  }
 
  GPUd() float evalUnsafe(const float x[/* Dim */]) const
  {
    int32_t index[Dim];
    setIndex<Dim - 1>(x, index);
    return MultivariatePolynomialHelper<Dim, Degree, InteractionOnly>::evalPol(getParameters(index), x);
  }
 
  GPUd() float evalPol(const float x[/* Dim */], const int32_t index[/* Dim */]) const
  {
    return MultivariatePolynomialHelper<Dim, Degree, InteractionOnly>::evalPol(getParameters(index), x);
  }
 
  GPUd() float getXMin(const uint32_t dim) const { return mMin[dim]; }
 
  GPUd() float getXMax(const uint32_t dim) const { return mMax[dim]; }
 
  GPUd() float getInvSpacing(const uint32_t dim) const { return mInvSpacing[dim]; }
 
  GPUd() uint32_t getNVertices(const uint32_t dim) const { return mN[dim]; }
 
  GPUd() uint32_t getNPolynomials(const uint32_t dim) const { return mN[dim] - 1; }
 
  GPUd() const float* getParams() const { return mParams; }
 
  void init(const float min[/* Dim */], const float max[/* Dim */], const uint32_t n[/* Dim */]);
 
  void setDefault();
 
#ifndef GPUCA_GPUCODE
  void setParams(const float params[/* getNParameters() */]) { std::copy(params, params + getNParameters(), mParams); }
#endif
 
#if !defined(GPUCA_GPUCODE) && !defined(GPUCA_STANDALONE)
  void performFits(const std::function<double(const double x[/* Dim */])>& func, const uint32_t nAuxiliaryPoints[/* Dim */]);
 
  void performFits(const std::vector<float>& x, const std::vector<float>& y);
 
  void loadFromFile(TFile& inpf, const char* name);
 
  void loadFromFile(const char* fileName, const char* name);
 
  void writeToFile(TFile& outf, const char* name) const;
 
  void dumpToTree(const uint32_t nSamplingPoints[/* Dim */], const char* outName = "debug.root", const char* treeName = "tree", const bool recreateFile = true) const;
 
#endif // !defined(GPUCA_GPUCODE) && !defined(GPUCA_STANDALONE)
#ifndef GPUCA_GPUCODE
 
  void setFromContainer(const NDPiecewisePolynomialContainer& container);
 
  uint32_t getNPolynomials() const;
 
  NDPiecewisePolynomialContainer getContainer() const { return NDPiecewisePolynomialContainer{Dim, Degree, getNParameters(), mParams, InteractionOnly, mMin, mMax, mN}; }
#endif
 
  uint32_t getNParameters() const { return getNPolynomials() * MultivariatePolynomialParametersHelper::getNParameters<Degree, Dim, InteractionOnly>(); }
 
  GPUd() static constexpr uint32_t getDim() { return Dim; }
 
  GPUd() static constexpr uint32_t getDegree() { return Degree; }
 
  GPUd() static constexpr bool isInteractionOnly() { return InteractionOnly; }
 
 private:
  using DataTParams = float;     
  float mMin[Dim]{};             
  float mMax[Dim]{};             
  float mInvSpacing[Dim]{};      
  uint32_t mN[Dim]{};            
  DataTParams* mParams{nullptr}; 
 
  GPUd() int32_t getVertex(const float x, const uint32_t dim) const { return ((x - mMin[dim]) * mInvSpacing[dim]); }
 
  template <uint32_t TermDim>
  GPUd() uint32_t getTerms() const
  {
    if constexpr (TermDim == 0) {
      return 1;
    } else {
      return (mN[TermDim - 1] - 1) * getTerms<TermDim - 1>();
    }
  }
 
  GPUd() uint32_t getDataIndex(const int32_t ix[/* Dim */]) const { return getDataIndex<Dim - 1>(ix) * MultivariatePolynomialParametersHelper::getNParameters<Degree, Dim, InteractionOnly>(); }
 
  template <uint32_t DimTmp>
  GPUd() uint32_t getDataIndex(const int32_t ix[/* Dim */]) const;
 
  GPUd() const float* getParameters(const int32_t index[/* Dim */]) const { return &mParams[getDataIndex(index)]; }
 
  template <uint32_t DimTmp>
  GPUd() void setIndex(const float x[/* Dim */], int32_t index[/* Dim */]) const;
 
  template <uint32_t DimTmp>
  GPUd() void clamp(float x[/* Dim */], int32_t index[/* Dim */]) const;
 
#if !defined(GPUCA_GPUCODE) && !defined(GPUCA_STANDALONE)
  void fitInnerGrid(const std::function<double(const double x[/* Dim */])>& func, const uint32_t nAuxiliaryPoints[/* Dim */], const int32_t currentIndex[/* Dim */], TLinearFitter& fitter, std::vector<double>& xCords, std::vector<double>& response);
#endif
 
  void checkPos(const uint32_t iMax[/* Dim */], int32_t pos[/* Dim */]) const;
 
  double getStepWidth(const uint32_t dim, const int32_t nAuxiliaryPoints) const { return 1 / (static_cast<double>(mInvSpacing[dim]) * (nAuxiliaryPoints - 1)); }
 
  double getVertexPosition(const uint32_t ix, const int32_t dim) const { return ix / static_cast<double>(mInvSpacing[dim]) + mMin[dim]; }
 
#if !defined(GPUCA_GPUCODE)
  std::size_t sizeOfParameters() const { return getNParameters() * sizeof(DataTParams); }
#endif // #if !defined(GPUCA_GPUCODE)
 
  // construct the object (flatbuffer)
  void construct();
 
  ClassDefNV(NDPiecewisePolynomials, 1);
};
 
//=================================================================================
//============================ inline implementations =============================
//=================================================================================
 
#if !defined(GPUCA_GPUCODE)
template <uint32_t Dim, uint32_t Degree, bool InteractionOnly>
void NDPiecewisePolynomials<Dim, Degree, InteractionOnly>::setFromContainer(const NDPiecewisePolynomialContainer& container)
{
  if (Dim != container.mDim) {
    LOGP(info, "wrong number of dimensions! this {} container {}", Dim, container.mDim);
    return;
  }
  if (Degree != container.mDegree) {
    LOGP(info, "wrong number of degrees! this {} container {}", Degree, container.mDegree);
    return;
  }
  if (InteractionOnly != container.mInteractionOnly) {
    LOGP(info, "InteractionOnly is set for this object to {}, but stored as {} in the container", InteractionOnly, container.mInteractionOnly);
    return;
  }
  init(container.mMin.data(), container.mMax.data(), container.mN.data());
  setParams(container.mParams.data());
}
 
template <uint32_t Dim, uint32_t Degree, bool InteractionOnly>
void NDPiecewisePolynomials<Dim, Degree, InteractionOnly>::setDefault()
{
  const auto nParamsPerPol = MultivariatePolynomialParametersHelper::getNParameters<Degree, Dim, InteractionOnly>();
  const auto nPols = getNPolynomials();
  std::vector<float> params(nParamsPerPol);
  params.front() = 1;
  for (unsigned int i = 0; i < nPols; ++i) {
    std::copy(params.begin(), params.end(), &mParams[i * nParamsPerPol]);
  }
}
 
template <uint32_t Dim, uint32_t Degree, bool InteractionOnly>
uint32_t NDPiecewisePolynomials<Dim, Degree, InteractionOnly>::getNPolynomials() const
{
  uint32_t nP = getNPolynomials(0);
  for (uint32_t i = 1; i < Dim; ++i) {
    nP *= getNPolynomials(i);
  }
  return nP;
}
 
template <uint32_t Dim, uint32_t Degree, bool InteractionOnly>
void NDPiecewisePolynomials<Dim, Degree, InteractionOnly>::checkPos(const uint32_t iMax[/* Dim */], int32_t pos[/* Dim */]) const
{
  for (uint32_t i = 0; i < Dim; ++i) {
    if (pos[i] == int32_t(iMax[i])) {
      ++pos[i + 1];
      std::fill_n(pos, i + 1, 0);
    }
  }
}
 
#endif // !defined(GPUCA_GPUCODE) && !defined(GPUCA_STANDALONE)
 
#ifndef GPUCA_GPUCODE
template <uint32_t Dim, uint32_t Degree, bool InteractionOnly>
void NDPiecewisePolynomials<Dim, Degree, InteractionOnly>::cloneFromObject(const NDPiecewisePolynomials<Dim, Degree, InteractionOnly>& obj, char* newFlatBufferPtr)
{
  const char* oldFlatBufferPtr = obj.mFlatBufferPtr;
  FlatObject::cloneFromObject(obj, newFlatBufferPtr);
  for (uint32_t i = 0; i < Dim; ++i) {
    mMin[i] = obj.mMin[i];
    mMax[i] = obj.mMax[i];
    mInvSpacing[i] = obj.mInvSpacing[i];
    mN[i] = obj.mN[i];
  }
  if (obj.mParams) {
    mParams = FlatObject::relocatePointer(oldFlatBufferPtr, mFlatBufferPtr, obj.mParams);
  }
}
 
template <uint32_t Dim, uint32_t Degree, bool InteractionOnly>
void NDPiecewisePolynomials<Dim, Degree, InteractionOnly>::moveBufferTo(char* newFlatBufferPtr)
{
  char* oldFlatBufferPtr = mFlatBufferPtr;
  FlatObject::moveBufferTo(newFlatBufferPtr);
  char* currFlatBufferPtr = mFlatBufferPtr;
  mFlatBufferPtr = oldFlatBufferPtr;
  setActualBufferAddress(currFlatBufferPtr);
}
 
template <uint32_t Dim, uint32_t Degree, bool InteractionOnly>
void NDPiecewisePolynomials<Dim, Degree, InteractionOnly>::construct()
{
  FlatObject::startConstruction();
  const std::size_t flatbufferSize = sizeOfParameters();
  FlatObject::finishConstruction(flatbufferSize);
  mParams = reinterpret_cast<DataTParams*>(mFlatBufferPtr);
}
 
template <uint32_t Dim, uint32_t Degree, bool InteractionOnly>
void NDPiecewisePolynomials<Dim, Degree, InteractionOnly>::init(const float min[], const float max[], const uint32_t n[])
{
  for (uint32_t i = 0; i < Dim; ++i) {
    mMin[i] = min[i];
    mMax[i] = max[i];
    mN[i] = n[i];
    mInvSpacing[i] = (mN[i] - 1) / (mMax[i] - mMin[i]);
  }
  construct();
}
#endif // !GPUCA_GPUCODE
 
template <uint32_t Dim, uint32_t Degree, bool InteractionOnly>
void NDPiecewisePolynomials<Dim, Degree, InteractionOnly>::destroy()
{
  mParams = nullptr;
  FlatObject::destroy();
}
 
template <uint32_t Dim, uint32_t Degree, bool InteractionOnly>
void NDPiecewisePolynomials<Dim, Degree, InteractionOnly>::setActualBufferAddress(char* actualFlatBufferPtr)
{
  FlatObject::setActualBufferAddress(actualFlatBufferPtr);
  mParams = reinterpret_cast<DataTParams*>(mFlatBufferPtr);
}
 
template <uint32_t Dim, uint32_t Degree, bool InteractionOnly>
void NDPiecewisePolynomials<Dim, Degree, InteractionOnly>::setFutureBufferAddress(char* futureFlatBufferPtr)
{
  mParams = FlatObject::relocatePointer(mFlatBufferPtr, futureFlatBufferPtr, mParams);
  FlatObject::setFutureBufferAddress(futureFlatBufferPtr);
}
 
template <uint32_t Dim, uint32_t Degree, bool InteractionOnly>
template <uint32_t DimTmp>
GPUd() uint32_t NDPiecewisePolynomials<Dim, Degree, InteractionOnly>::getDataIndex(const int32_t ix[/* Dim */]) const
{
  if constexpr (DimTmp > 0) {
    return ix[DimTmp] * getTerms<DimTmp>() + getDataIndex<DimTmp - 1>(ix);
  }
  return ix[DimTmp];
}
 
template <uint32_t Dim, uint32_t Degree, bool InteractionOnly>
template <uint32_t DimTmp>
GPUdi() void NDPiecewisePolynomials<Dim, Degree, InteractionOnly>::setIndex(const float x[/* Dim */], int32_t index[/* Dim */]) const
{
  index[DimTmp] = getVertex(x[DimTmp], DimTmp);
  if constexpr (DimTmp > 0) {
    return setIndex<DimTmp - 1>(x, index);
  }
  return;
}
 
template <uint32_t Dim, uint32_t Degree, bool InteractionOnly>
template <uint32_t DimTmp>
GPUdi() void NDPiecewisePolynomials<Dim, Degree, InteractionOnly>::clamp(float x[/* Dim */], int32_t index[/* Dim */]) const
{
  if (index[DimTmp] <= 0) {
    index[DimTmp] = 0;
 
    if (x[DimTmp] < mMin[DimTmp]) {
      x[DimTmp] = mMin[DimTmp];
    }
 
  } else {
    if (index[DimTmp] >= int32_t(mN[DimTmp] - 1)) {
      index[DimTmp] = mN[DimTmp] - 2;
      x[DimTmp] = mMax[DimTmp];
    }
  }
 
  if constexpr (DimTmp > 0) {
    return clamp<DimTmp - 1>(x, index);
  }
}
 
} // namespace o2::gpu
 
#endif