SplineSpec.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_GPUCOMMON_TPCFASTTRANSFORMATION_SPLINESPEC_H
#define ALICEO2_GPUCOMMON_TPCFASTTRANSFORMATION_SPLINESPEC_H
 
#include "Spline1D.h"
#include "FlatObject.h"
#include "GPUCommonDef.h"
#include "SplineUtil.h"
 
#if !defined(__CLING__) && !defined(G__ROOT) && !defined(GPUCA_GPUCODE) && !defined(GPUCA_NO_VC)
#include <Vc/Vc>
#include <Vc/SimdArray>
#endif
 
class TFile;
 
namespace o2
{
namespace gpu
{
 
template <typename DataT>
class SplineContainer : public FlatObject
{
 public:
  typedef typename Spline1D<DataT>::SafetyLevel SafetyLevel;
  typedef typename Spline1D<DataT>::Knot Knot;
 
 
  GPUd() static constexpr int32_t getVersion() { return (1 << 16) + Spline1D<DataT>::getVersion(); }
 
 
  SplineContainer() = default;
 
  SplineContainer(const SplineContainer&) = delete;
 
  ~SplineContainer() = default;
 
 
#if !defined(GPUCA_GPUCODE) && !defined(GPUCA_STANDALONE)
  void approximateFunction(const double xMin[/* mXdim */], const double xMax[/* mXdim */],
                           std::function<void(const double x[/* mXdim */], double f[/*mYdim*/])> F,
                           const int32_t nAuxiliaryDataPoints[/* mXdim */] = nullptr);
#endif
 
 
#if !defined(GPUCA_GPUCODE) && !defined(GPUCA_STANDALONE)
  int32_t writeToFile(TFile& outf, const char* name);
 
  static SplineContainer* readFromFile(TFile& inpf, const char* name);
#endif
 
 
  GPUd() int32_t getXdimensions() const { return mXdim; }
 
  GPUd() int32_t getYdimensions() const { return mYdim; }
 
  GPUd() static constexpr size_t getParameterAlignmentBytes() { return 16; }
 
  GPUd() int32_t getNumberOfParameters() const { return this->calcNumberOfParameters(mYdim); }
 
  GPUd() size_t getSizeOfParameters() const { return sizeof(DataT) * this->getNumberOfParameters(); }
 
  GPUd() int32_t getNumberOfKnots() const { return mNknots; }
 
  GPUd() int32_t getNumberOfParametersPerKnot() const { return calcNumberOfParametersPerKnot(mYdim); }
 
  GPUd() const Spline1D<DataT>& getGrid(int32_t dimX) const { return mGrid[dimX]; }
 
  GPUd() void getKnotU(int32_t iKnot, int32_t u[/* mXdim */]) const;
 
  GPUd() int32_t getKnotIndex(const int32_t iKnot[/* mXdim */]) const;
 
  GPUd() DataT* getParameters() { return mParameters; }
 
  GPUd() const DataT* getParameters() const { return mParameters; }
 
 
  GPUd() size_t getGridOffset(int32_t dimX) const { return mGrid[dimX].getFlatBufferPtr() - mFlatBufferPtr; }
 
  GPUd() void setXrange(const DataT xMin[/* mXdim */], const DataT xMax[/* mXdim */]);
 
  void print() const;
 
 
  GPUd() int32_t calcNumberOfParameters(int32_t nYdim) const
  {
    return calcNumberOfParametersPerKnot(nYdim) * getNumberOfKnots();
  }
 
  GPUd() int32_t calcNumberOfParametersPerKnot(int32_t nYdim) const
  {
    return (1 << mXdim) * nYdim; // 2^mXdim parameters per Y dimension
  }
 
 
#if !defined(GPUCA_GPUCODE) && !defined(GPUCA_STANDALONE) // code invisible on GPU and in the standalone compilation
  static int32_t test(const bool draw = 0, const bool drawDataPoints = 1);
#endif
 
 
  using FlatObject::getBufferAlignmentBytes;
  using FlatObject::getClassAlignmentBytes;
 
#if !defined(GPUCA_GPUCODE)
  void cloneFromObject(const SplineContainer& obj, char* newFlatBufferPtr);
  void moveBufferTo(char* newBufferPtr);
#endif
 
  using FlatObject::releaseInternalBuffer;
 
  void destroy();
  void setActualBufferAddress(char* actualFlatBufferPtr);
  void setFutureBufferAddress(char* futureFlatBufferPtr);
 
 protected:
#if !defined(GPUCA_GPUCODE)
  void recreate(int32_t nXdim, int32_t nYdim, const int32_t nKnots[/* nXdim */]);
 
  void recreate(int32_t nXdim, int32_t nYdim, const int32_t nKnots[/* nXdim */], const int32_t* const knotU[/* nXdim */]);
#endif
 
 
  int32_t mXdim = 0;   
  int32_t mYdim = 0;   
  int32_t mNknots = 0; 
 
  Spline1D<DataT>* mGrid; 
  DataT* mParameters;     
 
  ClassDefNV(SplineContainer, 1);
};
 
template <typename DataT>
GPUdi() void SplineContainer<DataT>::getKnotU(int32_t iKnot, int32_t u[/* mXdim */]) const
{
  for (int32_t dim = 0; dim < mXdim; dim++) {
    int32_t n = mGrid[dim].getNumberOfKnots();
    u[dim] = mGrid[dim].getKnot(iKnot % n).getU();
    iKnot /= n;
  }
}
 
template <typename DataT>
GPUdi() int32_t SplineContainer<DataT>::getKnotIndex(const int32_t iKnot[/* mXdim */]) const
{
  int32_t ind = iKnot[0];
  int32_t n = 1;
  for (int32_t dim = 1; dim < mXdim; dim++) {
    n *= mGrid[dim - 1].getNumberOfKnots();
    ind += n * iKnot[dim];
  }
  return ind;
}
 
template <typename DataT>
GPUdi() void SplineContainer<DataT>::
  setXrange(const DataT xMin[/* mXdim */], const DataT xMax[/* mXdim */])
{
  for (int32_t i = 0; i < mXdim; i++) {
    mGrid[i].setXrange(xMin[i], xMax[i]);
  }
}
 
template <typename DataT, int32_t XdimT, int32_t YdimT, int32_t SpecT>
class SplineSpec;
 
template <typename DataT, int32_t XdimT, int32_t YdimT>
class SplineSpec<DataT, XdimT, YdimT, 0> : public SplineContainer<DataT>
{
  typedef SplineContainer<DataT> TBase;
 
 public:
  typedef typename TBase::SafetyLevel SafetyLevel;
  typedef typename TBase::Knot Knot;
 
 
  GPUd() void interpolate(const DataT x[/*mXdim*/], GPUgeneric() DataT S[/*mYdim*/]) const
  {
    const auto nXdimTmp = SplineUtil::getNdim<XdimT>(mXdim);
    const auto nXdim = nXdimTmp.get();
    const auto maxXdimTmp = SplineUtil::getMaxNdim<XdimT>(mXdim);
    DataT u[maxXdimTmp.get()];
    for (int32_t i = 0; i < nXdim; i++) {
      u[i] = mGrid[i].convXtoU(x[i]);
    }
    interpolateU<SafetyLevel::kSafe>(mXdim, mYdim, mParameters, u, S);
  }
 
  template <SafetyLevel SafeT = SafetyLevel::kSafe>
  GPUd() void interpolateU(int32_t inpXdim, int32_t inpYdim, GPUgeneric() const DataT Parameters[],
                           const DataT u[/*inpXdim*/], GPUgeneric() DataT S[/*inpYdim*/]) const
  {
    const auto nXdimTmp = SplineUtil::getNdim<XdimT>(mXdim);
    const auto nXdim = nXdimTmp.get();
    const auto maxXdimTmp = SplineUtil::getMaxNdim<XdimT>(mXdim);
    const auto maxXdim = maxXdimTmp.get();
    const auto nYdimTmp = SplineUtil::getNdim<YdimT>(mYdim);
    const auto nYdim = nYdimTmp.get();
    const auto maxYdimTmp = SplineUtil::getMaxNdim<XdimT>(mYdim);
    const auto maxYdim = maxYdimTmp.get();
 
    // const auto nParameters = 1 << (2 * nXdim);         //total Nr of Parameters necessary for one interpolation
    const auto nKnotParametersPerY = 1 << nXdim;       // Nr of Parameters per Knot per Y dimension
    const auto nKnotParameters = (1 << nXdim) * nYdim; // Nr of Parameters per Knot
 
    DataT iParameters[(1 << (2 * maxXdim)) * maxYdim]; // Array for all parameters
 
    // get the indices of the "most left" Knot:
 
    int32_t indices[maxXdim]; // indices of the 'most left' knot
    for (int32_t i = 0; i < nXdim; i++) {
      indices[i] = mGrid[i].getLeftKnotIndexForU(u[i]);
    }
    // get all the needed parameters into one array iParameters[nParameters]:
    int32_t indicestmp[maxXdim];
    for (int32_t i = 0; i < nKnotParametersPerY; i++) { // for every necessary Knot
      for (int32_t k = 0; k < nXdim; k++) {
        indicestmp[k] = indices[k] + (i / (1 << k)) % 2; // get the knot-indices in every dimension (mirrored order binary counting)
      }
      int32_t index = TBase::getKnotIndex(indicestmp); // get index of the current Knot
 
      for (int32_t j = 0; j < nKnotParameters; j++) { // and fill the iparameter array with according parameters
        iParameters[i * nKnotParameters + j] = Parameters[index * nKnotParameters + j];
      }
    }
    // now start with the interpolation loop:
 
    constexpr auto maxInterpolations = (1 << (2 * maxXdim - 2)) * maxYdim;
 
    DataT S0[maxInterpolations];
    DataT D0[maxInterpolations];
    DataT S1[maxInterpolations];
    DataT D1[maxInterpolations];
 
    int32_t nInterpolations = (1 << (2 * nXdim - 2)) * nYdim;
    int32_t nKnots = 1 << (nXdim);
 
    for (int32_t d = 0; d < nXdim; d++) {            // for every dimension
      DataT* pointer[4] = {S0, D0, S1, D1};          // pointers for interpolation arrays S0, D0, S1, D1 point to Arraystart
      for (int32_t i = 0; i < nKnots; i++) {         // for every knot
        for (int32_t j = 0; j < nKnots; j++) {       // for every parametertype
          int32_t pointernr = 2 * (i % 2) + (j % 2); // to which array should it be delivered
          for (int32_t k = 0; k < nYdim; k++) {
            pointer[pointernr][0] = iParameters[(i * nKnots + j) * nYdim + k];
            pointer[pointernr]++;
          }
        } // end for j (every parametertype)
      } // end for i (every knot)
 
      const typename Spline1D<DataT>::Knot& knotL = mGrid[d].getKnot(indices[d]);
      DataT coordinate = u[d];
      typedef Spline1DSpec<DataT, 0, 0> TGridX;
      const TGridX& gridX = *((const TGridX*)&(mGrid[d]));
      gridX.interpolateU(nInterpolations, knotL, S0, D0, S1, D1, coordinate, iParameters);
      nInterpolations /= 4;
      nKnots /= 2;
    } // end d (every dimension)
 
    for (int32_t i = 0; i < nYdim; i++) {
      S[i] = iParameters[i]; // write into result-array
      // LOG(info)<<iParameters[i] <<", ";
    }
  } // end interpolateU
 
 protected:
  using TBase::mGrid;
  using TBase::mParameters;
  using TBase::mXdim;
  using TBase::mYdim;
  using TBase::TBase; // inherit constructors and hide them
};
 
template <typename DataT, int32_t XdimT, int32_t YdimT>
class SplineSpec<DataT, XdimT, YdimT, 1>
  : public SplineSpec<DataT, XdimT, YdimT, 0>
{
  typedef SplineContainer<DataT> TVeryBase;
  typedef SplineSpec<DataT, XdimT, YdimT, 0> TBase;
 
 public:
  typedef typename TVeryBase::SafetyLevel SafetyLevel;
 
#if !defined(GPUCA_GPUCODE)
  SplineSpec() : SplineSpec(nullptr) {}
 
  SplineSpec(const int32_t nKnots[/*XdimT*/]) : TBase()
  {
    recreate(nKnots);
  }
  SplineSpec(const int32_t nKnots[/*XdimT*/], const int32_t* const knotU[/*XdimT*/])
    : TBase()
  {
    recreate(nKnots, knotU);
  }
  SplineSpec(const SplineSpec& v) : TBase()
  {
    TBase::cloneFromObject(v, nullptr);
  }
  void recreate(const int32_t nKnots[/*XdimT*/])
  {
    TBase::recreate(XdimT, YdimT, nKnots);
  }
 
  void recreate(const int32_t nKnots[/*XdimT*/], const int32_t* const knotU[/*XdimT*/])
  {
    TBase::recreate(XdimT, YdimT, nKnots, knotU);
  }
#endif
 
  GPUd() constexpr int32_t getXdimensions() const { return XdimT; }
 
  GPUd() constexpr int32_t getYdimensions() const { return YdimT; }
 
 
  template <SafetyLevel SafeT = SafetyLevel::kSafe>
  GPUd() void interpolateU(GPUgeneric() const DataT Parameters[],
                           const DataT u[/*XdimT*/], GPUgeneric() DataT S[/*YdimT*/]) const
  {
    TBase::template interpolateU<SafeT>(XdimT, YdimT, Parameters, u, S);
  }
 
 private:
#if !defined(GPUCA_GPUCODE)
  using TBase::recreate;
#endif
  using TBase::interpolateU;
};
 
template <typename DataT, int32_t XdimT, int32_t YdimT>
class SplineSpec<DataT, XdimT, YdimT, 2>
  : public SplineSpec<DataT, XdimT, YdimT, 0>
{
  typedef SplineContainer<DataT> TVeryBase;
  typedef SplineSpec<DataT, XdimT, YdimT, 0> TBase;
 
 public:
  typedef typename TVeryBase::SafetyLevel SafetyLevel;
 
#if !defined(GPUCA_GPUCODE)
  SplineSpec() : SplineSpec((XdimT > 0 ? XdimT : 0), (YdimT > 0 ? YdimT : 0), nullptr) {}
 
  SplineSpec(int32_t nXdim, int32_t nYdim, const int32_t nKnots[/* nXdim */]) : TBase()
  {
    this->recreate(nXdim, nYdim, nKnots);
  }
 
  SplineSpec(int32_t nXdim, int32_t nYdim, const int32_t nKnots[/* nXdim */], const int32_t* const knotU[/* nXdim */])
    : TBase()
  {
    this->recreate(nXdim, nYdim, nKnots, knotU);
  }
 
  SplineSpec(const SplineSpec& v) : TBase()
  {
    cloneFromObject(v, nullptr);
  }
 
  void recreate(int32_t nXdim, int32_t nYdim, const int32_t nKnots[/* nXdim */])
  {
    checkDimensions(nXdim, nYdim);
    TBase::recreate(nXdim, nYdim, nKnots);
  }
 
  void recreate(int32_t nXdim, int32_t nYdim, const int32_t nKnots[/* nXdim */], const int32_t* const knotU[/* nXdim */])
  {
    checkDimensions(nXdim, nYdim);
    TBase::recreate(nXdim, nYdim, nKnots, knotU);
  }
 
#endif
 
 
  using TBase::interpolateU;
 
  void checkDimensions(int32_t& nXdim, int32_t& nYdim)
  {
    if (XdimT > 0 && nXdim != XdimT) {
      assert(0);
      nXdim = XdimT;
    }
    if (XdimT < 0 && nXdim > abs(XdimT)) {
      assert(0);
      nXdim = abs(XdimT);
    }
    if (nXdim < 0) {
      assert(0);
      nXdim = 0;
    }
    if (YdimT > 0 && nYdim != YdimT) {
      assert(0);
      nYdim = YdimT;
    }
    if (YdimT < 0 && nYdim > abs(YdimT)) {
      assert(0);
      nYdim = abs(YdimT);
    }
    if (nYdim < 0) {
      assert(0);
      nYdim = 0;
    }
  }
};
 
template <typename DataT, int32_t XdimT>
class SplineSpec<DataT, XdimT, 1, 3>
  : public SplineSpec<DataT, XdimT, 1, SplineUtil::getSpec(XdimT, 999)>
{
  typedef SplineSpec<DataT, XdimT, 1, SplineUtil::getSpec(XdimT, 999)> TBase;
 
 public:
  using TBase::TBase; // inherit constructors
 
  GPUd() DataT interpolate(const DataT x[]) const
  {
    DataT S = 0.;
    TBase::interpolate(x, &S);
    return S;
  }
 
  // this parent method should be public anyhow,
  // but the compiler gets confused w/o this extra declaration
  using TBase::interpolate;
};
 
} // namespace gpu
} // namespace o2
 
#endif