Spline1DSpec.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_SPLINE1DSPEC_H
#define ALICEO2_GPUCOMMON_TPCFASTTRANSFORMATION_SPLINE1DSPEC_H
 
#include "GPUCommonDef.h"
#include "FlatObject.h"
#include "SplineUtil.h"
 
#if !defined(GPUCA_GPUCODE)
#include <functional>
#endif
 
class TFile;
 
namespace o2
{
namespace gpu
{
 
template <typename DataT>
class Spline1DContainer : public FlatObject
{
 public:
  enum SafetyLevel { kNotSafe,
                     kSafe };
 
  struct Knot {
    DataT u;  
    DataT Li; 
    GPUd() int32_t getU() const { return (int32_t)(u + 0.1f); }
  };
 
 
  GPUd() static constexpr int32_t getVersion() { return 1; }
 
 
  Spline1DContainer() = default;
 
  Spline1DContainer(const Spline1DContainer&) = delete;
 
  ~Spline1DContainer() = default;
 
 
#if !defined(GPUCA_GPUCODE) && !defined(GPUCA_STANDALONE)
  void approximateFunction(double xMin, double xMax,
                           std::function<void(double x, double f[/*mYdim*/])> F,
                           int32_t nAuxiliaryDataPoints = 4);
#endif
 
 
#if !defined(GPUCA_GPUCODE) && !defined(GPUCA_STANDALONE)
  int32_t writeToFile(TFile& outf, const char* name);
 
  static Spline1DContainer* readFromFile(TFile& inpf, const char* name);
#endif
 
 
  GPUd() int32_t getUmax() const { return mUmax; }
 
  GPUd() int32_t getYdimensions() const { return mYdim; }
 
  GPUd() size_t getParameterAlignmentBytes() const
  {
    size_t s = 2 * sizeof(DataT) * mYdim;
    return (s < 16) ? s : 16;
  }
 
  GPUd() int32_t getNumberOfParameters() const { return calcNumberOfParameters(mYdim); }
 
  GPUd() size_t getSizeOfParameters() const { return sizeof(DataT) * getNumberOfParameters(); }
 
  GPUd() int32_t getNumberOfKnots() const { return mNumberOfKnots; }
 
  GPUd() const Knot* getKnots() const { return reinterpret_cast<const Knot*>(mFlatBufferPtr); }
 
  template <SafetyLevel SafeT = SafetyLevel::kSafe>
  GPUd() const Knot& getKnot(int32_t i) const
  {
    if (SafeT == SafetyLevel::kSafe) {
      i = (i < 0) ? 0 : (i >= mNumberOfKnots ? mNumberOfKnots - 1 : i);
    }
    return getKnots()[i];
  }
 
  template <SafetyLevel SafeT = SafetyLevel::kSafe>
  GPUd() int32_t getLeftKnotIndexForU(DataT u) const;
 
  GPUd() DataT* getParameters() { return mParameters; }
 
  GPUd() const DataT* getParameters() const { return mParameters; }
 
 
  GPUd() const int32_t* getUtoKnotMap() const { return mUtoKnotMap; }
 
  GPUd() DataT convXtoU(DataT x) const { return (x - mXmin) * mXtoUscale; }
 
  GPUd() DataT convUtoX(DataT u) const { return mXmin + u / mXtoUscale; }
 
  GPUd() DataT getXmin() const { return mXmin; }
 
  GPUd() DataT getXmax() const { return mXmin + mUmax / mXtoUscale; }
 
  GPUd() DataT getXtoUscale() const { return mXtoUscale; }
 
  GPUd() void setXrange(DataT xMin, DataT xMax);
 
  void print() const;
 
 
  GPUd() int32_t calcNumberOfParameters(int32_t nYdim) const { return (2 * nYdim) * getNumberOfKnots(); }
 
 
#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 Spline1DContainer& obj, char* newFlatBufferPtr);
  void moveBufferTo(char* newBufferPtr);
#endif
 
  using FlatObject::releaseInternalBuffer;
 
  void destroy();
  void setActualBufferAddress(char* actualFlatBufferPtr);
  void setFutureBufferAddress(char* futureFlatBufferPtr);
 
 protected:
  Knot* getKnots() { return reinterpret_cast<Knot*>(mFlatBufferPtr); }
 
  int32_t* getUtoKnotMap() { return mUtoKnotMap; }
 
#if !defined(GPUCA_GPUCODE)
  void recreate(int32_t nYdim, int32_t numberOfKnots);
 
  void recreate(int32_t nYdim, int32_t numberOfKnots, const int32_t knotU[]);
#endif
 
 
  int32_t mYdim = 0;            
  int32_t mNumberOfKnots = 0;   
  int32_t mUmax = 0;            
  DataT mXmin = 0;              
  DataT mXtoUscale = 0;         
  int32_t* mUtoKnotMap = nullptr; 
  DataT* mParameters = nullptr; 
 
  ClassDefNV(Spline1DContainer, 1);
};
 
template <typename DataT>
template <typename Spline1DContainer<DataT>::SafetyLevel SafeT>
GPUdi() int32_t Spline1DContainer<DataT>::getLeftKnotIndexForU(DataT u) const
{
  int32_t iu = u < 0 ? 0 : (u > (float)mUmax ? mUmax : (int32_t)u);
  if (SafeT == SafetyLevel::kSafe) {
    iu = (iu < 0) ? 0 : (iu > mUmax ? mUmax : iu);
  }
  return getUtoKnotMap()[iu];
}
 
template <typename DataT>
GPUdi() void Spline1DContainer<DataT>::setXrange(DataT xMin, DataT xMax)
{
  mXmin = xMin;
  double l = ((double)xMax) - xMin;
  if (l < 1.e-8) {
    l = 1.e-8;
  }
  mXtoUscale = mUmax / l;
}
 
template <typename DataT, int32_t YdimT, int32_t SpecT>
class Spline1DSpec;
 
template <typename DataT, int32_t YdimT>
class Spline1DSpec<DataT, YdimT, 0> : public Spline1DContainer<DataT>
{
  typedef Spline1DContainer<DataT> TBase;
 
 public:
  typedef typename TBase::SafetyLevel SafetyLevel;
  typedef typename TBase::Knot Knot;
 
 
  GPUd() void interpolate(DataT x, GPUgeneric() DataT S[/*mYdim*/]) const
  {
    interpolateU<SafetyLevel::kSafe>(mYdim, mParameters, convXtoU(x), S);
  }
 
  template <SafetyLevel SafeT = SafetyLevel::kSafe>
  GPUd() void interpolateU(int32_t inpYdim, GPUgeneric() const DataT Parameters[],
                           DataT u, GPUgeneric() DataT S[/*nYdim*/]) const
  {
    const auto nYdimTmp = SplineUtil::getNdim<YdimT>(inpYdim);
    const auto nYdim = nYdimTmp.get();
    int32_t iknot = TBase::template getLeftKnotIndexForU<SafeT>(u);
    const DataT* d = Parameters + (2 * nYdim) * iknot;
    interpolateU(nYdim, getKnots()[iknot], &(d[0]), &(d[nYdim]), &(d[2 * nYdim]), &(d[3 * nYdim]), u, S);
  }
 
  template <typename T>
  GPUd() void interpolateU(int32_t inpYdim, const Knot& knotL,
                           GPUgeneric() const T Sl[/*mYdim*/], GPUgeneric() const T Dl[/*mYdim*/],
                           GPUgeneric() const T Sr[/*mYdim*/], GPUgeneric() const T Dr[/*mYdim*/],
                           DataT u, GPUgeneric() T S[/*mYdim*/]) const
  {
    const auto nYdimTmp = SplineUtil::getNdim<YdimT>(inpYdim);
    const auto nYdim = nYdimTmp.get();
    T uu = T(u - knotL.u);
    T li = T(knotL.Li);
    T v = uu * li; // scaled u
    for (int32_t dim = 0; dim < nYdim; ++dim) {
      T df = (Sr[dim] - Sl[dim]) * li;
      T a = Dl[dim] + Dr[dim] - df - df;
      T b = df - Dl[dim] - a;
      S[dim] = ((a * v + b) * v + Dl[dim]) * uu + Sl[dim];
    }
    /*
     another way to calculate f(u):
     T uu = T(u - knotL.u);
     T v = uu * T(knotL.Li); // scaled u
     T vm1 = v-1;
     T v2 = v * v;
     float cSr = v2*(3-2*v);
     float cSl = 1-cSr;
     float cDl = v*vm1*vm1*knotL.L;
     float cDr = v2*vm1*knotL.L;
     return cSl*Sl + cSr*Sr + cDl*Dl + cDr*Dr;
    */
  }
 
  template <typename T>
  GPUd() static void getUderivatives(const Knot& knotL, DataT u,
                                     T& dSl, T& dDl, T& dSr, T& dDr)
  {
    u = u - knotL.u;
    T v = u * T(knotL.Li); // scaled u
    T vm1 = v - 1.;
    T a = u * vm1;
    T v2 = v * v;
    dSr = v2 * (3. - 2 * v);
    dSl = 1. - dSr;
    dDl = vm1 * a;
    dDr = v * a;
    // F(u) = dSl * Sl + dSr * Sr + dDl * Dl + dDr * Dr;
  }
 
  using TBase::convXtoU;
  using TBase::getKnot;
  using TBase::getKnots;
  using TBase::getNumberOfKnots;
 
 protected:
  using TBase::mParameters;
  using TBase::mYdim;
  using TBase::TBase; // inherit constructors and hide them
  ClassDefNV(Spline1DSpec, 0);
};
 
template <typename DataT, int32_t YdimT>
class Spline1DSpec<DataT, YdimT, 1>
  : public Spline1DSpec<DataT, YdimT, 0>
{
  typedef Spline1DContainer<DataT> TVeryBase;
  typedef Spline1DSpec<DataT, YdimT, 0> TBase;
 
 public:
  typedef typename TVeryBase::SafetyLevel SafetyLevel;
 
#if !defined(GPUCA_GPUCODE)
  Spline1DSpec() : Spline1DSpec(2) {}
 
  Spline1DSpec(int32_t numberOfKnots) : TBase()
  {
    recreate(numberOfKnots);
  }
  Spline1DSpec(int32_t numberOfKnots, const int32_t knotU[])
    : TBase()
  {
    recreate(numberOfKnots, knotU);
  }
  Spline1DSpec(const Spline1DSpec& v) : TBase()
  {
    TBase::cloneFromObject(v, nullptr);
  }
  void recreate(int32_t numberOfKnots) { TBase::recreate(YdimT, numberOfKnots); }
 
  void recreate(int32_t numberOfKnots, const int32_t knotU[])
  {
    TBase::recreate(YdimT, numberOfKnots, knotU);
  }
#endif
 
  GPUd() constexpr int32_t getYdimensions() const { return YdimT; }
 
  GPUd() constexpr size_t getParameterAlignmentBytes() const
  {
    size_t s = 2 * sizeof(DataT) * YdimT;
    return (s < 16) ? s : 16;
  }
 
  GPUd() int32_t getNumberOfParameters() const { return (2 * YdimT) * getNumberOfKnots(); }
 
  GPUd() size_t getSizeOfParameters() const { return (sizeof(DataT) * 2 * YdimT) * getNumberOfKnots(); }
 
 
  template <SafetyLevel SafeT = SafetyLevel::kSafe>
  GPUd() void interpolateU(GPUgeneric() const DataT Parameters[],
                           DataT u, GPUgeneric() DataT S[/*nYdim*/]) const
  {
    TBase::template interpolateU<SafeT>(YdimT, Parameters, u, S);
  }
 
  template <typename T>
  GPUd() void interpolateU(const typename TBase::Knot& knotL,
                           GPUgeneric() const T Sl[/*mYdim*/], GPUgeneric() const T Dl[/*mYdim*/],
                           GPUgeneric() const T Sr[/*mYdim*/], GPUgeneric() const T Dr[/*mYdim*/],
                           DataT u, GPUgeneric() T S[/*mYdim*/]) const
  {
    TBase::interpolateU(YdimT, knotL, Sl, Dl, Sr, Dr, u, S);
  }
 
  using TBase::getNumberOfKnots;
 
 private:
#if !defined(GPUCA_GPUCODE)
  using TBase::recreate;
#endif
  using TBase::interpolateU;
};
 
template <typename DataT, int32_t YdimT>
class Spline1DSpec<DataT, YdimT, 2>
  : public Spline1DSpec<DataT, YdimT, 0>
{
  typedef Spline1DContainer<DataT> TVeryBase;
  typedef Spline1DSpec<DataT, YdimT, 0> TBase;
 
 public:
  typedef typename TVeryBase::SafetyLevel SafetyLevel;
 
#if !defined(GPUCA_GPUCODE)
  Spline1DSpec() : Spline1DSpec(0, 2) {}
 
  Spline1DSpec(int32_t nYdim, int32_t numberOfKnots) : TBase()
  {
    TBase::recreate(nYdim, numberOfKnots);
  }
  Spline1DSpec(int32_t nYdim, int32_t numberOfKnots, const int32_t knotU[]) : TBase()
  {
    TBase::recreate(nYdim, numberOfKnots, knotU);
  }
  Spline1DSpec(const Spline1DSpec& v) : TBase()
  {
    TVeryBase::cloneFromObject(v, nullptr);
  }
  void recreate(int32_t nYdim, int32_t numberOfKnots) { TBase::recreate(nYdim, numberOfKnots); }
 
  void recreate(int32_t nYdim, int32_t numberOfKnots, const int32_t knotU[])
  {
    TBase::recreate(nYdim, numberOfKnots, knotU);
  }
#endif
 
 
  using TBase::interpolateU;
  ClassDefNV(Spline1DSpec, 0);
};
 
template <typename DataT>
class Spline1DSpec<DataT, 1, 3>
  : public Spline1DSpec<DataT, 1, SplineUtil::getSpec(999)>
{
  typedef Spline1DSpec<DataT, 1, SplineUtil::getSpec(999)> TBase;
 
 public:
  using TBase::TBase; // inherit constructors
 
  GPUd() DataT interpolate(DataT x) const
  {
    DataT S = 0.;
    TBase::interpolate(x, &S);
    return S;
  }
};
 
} // namespace gpu
} // namespace o2
 
#endif