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>
#include <type_traits>
#endif
 
class TFile;
 
namespace o2::gpu
{
 
template <typename DataT>
struct Knot {
  DataT u;  
  DataT Li; 
  GPUd() int32_t getU() const { return (int32_t)(u + 0.1f); }
};
 
enum SafetyLevel { kNotSafe,
                   kSafe };
 
template <typename DataT, class FlatBase = FlatObject>
class Spline1DContainerBase : public FlatBase
{
 public:
 
  GPUd() static constexpr int32_t getVersion() { return 1; }
 
 
  Spline1DContainerBase() = default;
 
  Spline1DContainerBase(const Spline1DContainerBase&) = delete;
 
  ~Spline1DContainerBase() = default;
 
 
  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() 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);
 
 
  GPUd() int32_t calcNumberOfParameters(int32_t nYdim) const { return (2 * nYdim) * getNumberOfKnots(); }
 
  int32_t mYdim = 0;          
  int32_t mNumberOfKnots = 0; 
  int32_t mUmax = 0;          
  DataT mXmin = 0;            
  DataT mXtoUscale = 0;       
};
 
template <typename DataT, typename FlatBase = FlatObject>
class Spline1DContainer; // forward declaration
 
template <typename DataT>
class Spline1DContainer<DataT, FlatObject> : public Spline1DContainerBase<DataT, FlatObject>
{
 public:
  GPUd() const int32_t* getUtoKnotMap() const { return mUtoKnotMap; }
 
  GPUd() const Knot<DataT>* getKnots() const { return reinterpret_cast<const Knot<DataT>*>(this->mFlatBufferPtr); }
 
  template <SafetyLevel SafeT = SafetyLevel::kSafe>
  GPUd() const Knot<DataT>& getKnot(int32_t i) const
  {
    if (SafeT == SafetyLevel::kSafe) {
      i = (i < 0) ? 0 : (i >= this->getNumberOfKnots() ? this->getNumberOfKnots() - 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; }
 
#if !defined(GPUCA_GPUCODE) && !defined(GPUCA_STANDALONE)
  void approximateFunction(double xMin, double xMax, std::function<void(double x, double f[])> F, int32_t nAuxiliaryDataPoints = 4);
 
  int32_t writeToFile(TFile& outf, const char* name);
 
  static Spline1DContainer* readFromFile(TFile& inpf, const char* name);
 
  static int32_t test(const bool draw = 0, const bool drawDataPoints = 1);
#endif
 
  void print() const;
 
#if !defined(GPUCA_GPUCODE)
  void cloneFromObject(const Spline1DContainer& obj, char* newFlatBufferPtr);
  void moveBufferTo(char* newBufferPtr);
 
  template <class OtherFlatBase>
  void importFrom(const Spline1DContainerBase<DataT, OtherFlatBase>& src);
#endif
 
  void destroy();
  void setActualBufferAddress(char* actualFlatBufferPtr);
  void setFutureBufferAddress(char* futureFlatBufferPtr);
 
 protected:
#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* getUtoKnotMap() { return mUtoKnotMap; }
 
  Knot<DataT>* getKnots() { return reinterpret_cast<Knot<DataT>*>(this->mFlatBufferPtr); }
 
  int32_t* mUtoKnotMap = nullptr; 
  DataT* mParameters = nullptr;   
 
  ClassDefNV(Spline1DContainer, 1);
};
 
template <typename DataT>
class Spline1DContainer<DataT, NoFlatObject> : public Spline1DContainerBase<DataT, NoFlatObject>
{
 public:
  GPUd() const int32_t* getUtoKnotMapFromBuffer(const char* flatBuf) const
  {
    return reinterpret_cast<const int32_t*>(flatBuf + this->getNumberOfKnots() * sizeof(Knot<DataT>));
  }
 
  template <SafetyLevel SafeT = SafetyLevel::kSafe>
  GPUd() int32_t getLeftKnotIndexForUFromBuffer(const char* flatBuf, DataT u) const
  {
    int32_t iu = u < 0 ? 0 : (u > (float)this->mUmax ? this->mUmax : (int32_t)u);
    if (SafeT == SafetyLevel::kSafe) {
      iu = (iu < 0) ? 0 : (iu > this->mUmax ? this->mUmax : iu);
    }
    return getUtoKnotMapFromBuffer(flatBuf)[iu];
  }
 
  GPUd() const Knot<DataT>* getKnotsFromBuffer(const char* flatBuf) const { return reinterpret_cast<const Knot<DataT>*>(flatBuf); }
 
  template <SafetyLevel SafeT = SafetyLevel::kSafe>
  GPUd() const Knot<DataT>& getKnotFromBuffer(const char* flatBuf, int32_t i) const
  {
    if (SafeT == SafetyLevel::kSafe) {
      i = (i < 0) ? 0 : (i >= this->getNumberOfKnots() ? this->getNumberOfKnots() - 1 : i);
    }
    return getKnotsFromBuffer(flatBuf)[i];
  }
 
  // Lifecycle no-ops: NoFlatObject splines have no owned buffer.
  void destroy()
  {
    this->mNumberOfKnots = 0;
    this->mUmax = 0;
    this->mYdim = 0;
    this->mXmin = 0.;
    this->mXtoUscale = 1.;
    this->mFlatBufferSize = 0;
  }
  GPUdi() void setActualBufferAddress(char*) {}
  GPUdi() void setFutureBufferAddress(char*) {}
 
#if !defined(GPUCA_GPUCODE)
  template <class OtherFlatBase>
  void importFrom(const Spline1DContainerBase<DataT, OtherFlatBase>& src)
  {
    this->mYdim = src.getYdimensions();
    this->mNumberOfKnots = src.getNumberOfKnots();
    this->mUmax = src.getUmax();
    this->mXmin = src.getXmin();
    this->mXtoUscale = src.getXtoUscale();
    this->mFlatBufferSize = src.getFlatBufferSize();
  }
#endif
};
 
template <typename DataT>
template <SafetyLevel SafeT>
GPUdi() int32_t Spline1DContainer<DataT, FlatObject>::getLeftKnotIndexForU(DataT u) const
{
  int32_t iu = u < 0 ? 0 : (u > (float)this->mUmax ? this->mUmax : (int32_t)u);
  if (SafeT == SafetyLevel::kSafe) {
    iu = (iu < 0) ? 0 : (iu > this->mUmax ? this->mUmax : iu);
  }
  return getUtoKnotMap()[iu];
}
 
template <typename DataT, class FlatBase>
GPUdi() void Spline1DContainerBase<DataT, FlatBase>::setXrange(DataT xMin, DataT xMax)
{
  mXmin = xMin;
  double l = ((double)xMax) - xMin;
  if (l < 1.e-8) {
    l = 1.e-8;
  }
  mXtoUscale = this->mUmax / l;
}
 
template <typename DataT, int32_t YdimT, int32_t SpecT, class FlatBase = FlatObject>
class Spline1DSpec;
 
template <typename DataT, int32_t YdimT, class FlatBase>
class Spline1DSpec<DataT, YdimT, 0, FlatBase> : public Spline1DContainer<DataT, FlatBase>
{
  using Container = Spline1DContainer<DataT, FlatBase>;
 
 public:
  using KnotType = Knot<DataT>;
 
 
  GPUd() void interpolate(DataT x, GPUgeneric() DataT S[/*mYdim*/]) const
  {
    if constexpr (std::is_same_v<FlatBase, FlatObject>) {
      interpolateAtU<SafetyLevel::kSafe>(this->mYdim, this->mParameters, this->convXtoU(x), S);
    }
  }
 
  template <SafetyLevel SafeT = SafetyLevel::kSafe>
  GPUd() void interpolateAtU(int32_t inpYdim, GPUgeneric() const DataT Parameters[], DataT u, GPUgeneric() DataT S[/*nYdim*/]) const
  {
    if constexpr (std::is_same_v<FlatBase, FlatObject>) {
      const auto nYdimTmp = SplineUtil::getNdim<YdimT>(inpYdim);
      const auto nYdim = nYdimTmp.get();
      int32_t iknot = this->template getLeftKnotIndexForU<SafeT>(u);
      const DataT* d = Parameters + (2 * nYdim) * iknot;
      interpolateAtU(nYdim, this->getKnots()[iknot], &(d[0]), &(d[nYdim]), &(d[2 * nYdim]), &(d[3 * nYdim]), u, S);
    }
  }
 
  template <typename T>
  GPUd() void interpolateAtU(int32_t inpYdim, const KnotType& 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 dSdSl, dSdDl, dSdSr, dSdDr;
    getSderivativesOverParsAtU<T>(knotL, u, dSdSl, dSdDl, dSdSr, dSdDr);
 
    for (int32_t dim = 0; dim < nYdim; ++dim) {
      S[dim] = dSdSr * Sr[dim] + dSdSl * Sl[dim] + dSdDl * Dl[dim] + dSdDr * Dr[dim];
    }
 
    /*
    another way to calculate f(u):
 
    if (u < (DataT)0) {
      u = (DataT)0;
    }
    if (u > (DataT)ParentSpec::getUmax()) {
      u = (DataT)ParentSpec::getUmax();
    }
 
    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];
    }
    */
  }
 
  template <typename T>
  GPUd() void getSderivativesOverParsAtU(const KnotType& knotL, DataT u, T& dSdSl, T& dSdDl, T& dSdSr, T& dSdDr) const
  {
 
    if (u < (DataT)0) {
      u = (DataT)0;
    }
    if (u > (DataT)Container::getUmax()) {
      u = (DataT)Container::getUmax();
    }
 
    u = u - knotL.u;
    T v = u * T(knotL.Li); // scaled u
    T vm1 = v - T(1.);
    T a = u * vm1;
    T v2 = v * v;
    dSdSr = v2 * (T(3.) - v - v);
    dSdSl = T(1.) - dSdSr;
    dSdDl = vm1 * a;
    dSdDr = v * a;
    // S(u) = dSdSl * Sl + dSdSr * Sr + dSdDl * Dl + dSdDr * Dr;
  }
 
  template <typename T>
  GPUd() void getSDderivativesOverParsAtU(const KnotType& knotL, DataT u, T& dSdSl, T& dSdDl, T& dSdSr, T& dSdDr, T& dDdSl, T& dDdDl, T& dDdSr, T& dDdDr) const
  {
 
    if (u < (DataT)0) {
      u = (DataT)0;
    }
    if (u > (DataT)Container::getUmax()) {
      u = (DataT)Container::getUmax();
    }
 
    u = u - knotL.u;
    T v = u * T(knotL.Li); // scaled u
    T vm1 = v - T(1.);
    T a = u * vm1;
    T v2 = v * v;
    dSdSr = v2 * (T(3.) - v - v);
    dSdSl = T(1.) - dSdSr;
    dSdDl = vm1 * a;
    dSdDr = v * a;
 
    T dv = T(knotL.Li);
    dDdSr = T(6.) * v * (T(1.) - v) * dv;
    dDdSl = -dDdSr;
    dDdDl = vm1 * (v + v + vm1);
    dDdDr = v * (v + vm1 + vm1);
    // S(u) = dSdSl * Sl + dSdSr * Sr + dSdDl * Dl + dSdDr * Dr;
    // D(u) = dS(u)/du = dDdSl * Sl + dDdSr * Sr + dDdDl * Dl + dDdDr * Dr;
  }
};
 
template <typename DataT, int32_t YdimT, class FlatBase>
class Spline1DSpec<DataT, YdimT, 1, FlatBase> : public Spline1DSpec<DataT, YdimT, 0, FlatBase>
{
  using ParentSpec = Spline1DSpec<DataT, YdimT, 0, FlatBase>;
 
 public:
#if !defined(GPUCA_GPUCODE)
  Spline1DSpec() : ParentSpec()
  {
    if constexpr (!std::is_same_v<FlatBase, NoFlatObject>) {
      recreate(2);
    }
  }
 
  Spline1DSpec(int32_t numberOfKnots) : ParentSpec() { recreate(numberOfKnots); }
 
  Spline1DSpec(int32_t numberOfKnots, const int32_t knotU[]) : ParentSpec() { recreate(numberOfKnots, knotU); }
 
  Spline1DSpec(const Spline1DSpec& v) : ParentSpec() { ParentSpec::cloneFromObject(v, nullptr); }
 
  void recreate(int32_t numberOfKnots) { ParentSpec::recreate(YdimT, numberOfKnots); }
 
  void recreate(int32_t numberOfKnots, const int32_t knotU[]) { ParentSpec::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) * this->getNumberOfKnots(); }
 
  GPUd() size_t getSizeOfParameters() const { return (sizeof(DataT) * 2 * YdimT) * this->getNumberOfKnots(); }
 
 
  template <SafetyLevel SafeT = SafetyLevel::kSafe>
  GPUd() void interpolateAtU(GPUgeneric() const DataT Parameters[], DataT u, GPUgeneric() DataT S[/*nYdim*/]) const
  {
    ParentSpec::template interpolateAtU<SafeT>(YdimT, Parameters, u, S);
  }
 
  template <typename T>
  GPUd() void interpolateAtU(const typename ParentSpec::KnotType& 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
  {
    ParentSpec::interpolateAtU(YdimT, knotL, Sl, Dl, Sr, Dr, u, S);
  }
};
 
template <typename DataT, int32_t YdimT, class FlatBase>
class Spline1DSpec<DataT, YdimT, 2, FlatBase> : public Spline1DSpec<DataT, YdimT, 0, FlatBase>
{
  using ParentSpec = Spline1DSpec<DataT, YdimT, 0, FlatBase>;
  using Container = Spline1DContainer<DataT, FlatBase>;
 
 public:
#if !defined(GPUCA_GPUCODE)
  Spline1DSpec() : ParentSpec()
  {
    if constexpr (!std::is_same_v<FlatBase, NoFlatObject>) {
      ParentSpec::recreate(0, 2);
    }
  }
 
  Spline1DSpec(int32_t nYdim, int32_t numberOfKnots) : ParentSpec()
  {
    ParentSpec::recreate(nYdim, numberOfKnots);
  }
 
  Spline1DSpec(int32_t nYdim, int32_t numberOfKnots, const int32_t knotU[]) : ParentSpec()
  {
    ParentSpec::recreate(nYdim, numberOfKnots, knotU);
  }
 
  Spline1DSpec(const Spline1DSpec& v) : ParentSpec()
  {
    Container::cloneFromObject(v, nullptr);
  }
 
  void recreate(int32_t nYdim, int32_t numberOfKnots) { ParentSpec::recreate(nYdim, numberOfKnots); }
 
  void recreate(int32_t nYdim, int32_t numberOfKnots, const int32_t knotU[]) { ParentSpec::recreate(nYdim, numberOfKnots, knotU); }
#endif
};
 
template <typename DataT, class FlatBase>
class Spline1DSpec<DataT, 1, 3, FlatBase> : public Spline1DSpec<DataT, 1, SplineUtil::getSpec(999), FlatBase>
{
  using ParentSpec = Spline1DSpec<DataT, 1, SplineUtil::getSpec(999), FlatBase>;
 
 public:
  GPUd() DataT interpolate(DataT x) const
  {
    DataT S = 0;
    ParentSpec::interpolate(x, &S);
    return S;
  }
};
 
} // namespace o2::gpu
 
#endif