IrregularSpline1D.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_IRREGULARSPLINE1D_H
#define ALICEO2_GPUCOMMON_TPCFASTTRANSFORMATION_IRREGULARSPLINE1D_H
 
#include "GPUCommonDef.h"
#include "FlatObject.h"
#include "GPUCommonDef.h"
 
#ifndef __OPENCL__
#include <cstddef>
#include <memory>
#include <cstring>
#endif
 
namespace o2
{
namespace gpu
{
 
class IrregularSpline1D : public FlatObject
{
 public:
  struct Knot {
    float u; 
 
 
    float scale;   
    float scaleL0; 
    float scaleL2; 
    float scaleR2; 
    float scaleR3; 
  };
 
 
  IrregularSpline1D();
 
  IrregularSpline1D(const IrregularSpline1D&) = delete;
 
  IrregularSpline1D& operator=(const IrregularSpline1D&) = delete;
 
  ~IrregularSpline1D() = default;
 
 
 
  using FlatObject::getBufferAlignmentBytes;
  using FlatObject::getClassAlignmentBytes;
 
 
  void cloneFromObject(const IrregularSpline1D& obj, char* newFlatBufferPtr);
  void destroy();
 
 
  using FlatObject::releaseInternalBuffer;
#ifndef GPUCA_GPUCODE
  using FlatObject::moveBufferTo;
#endif
 
 
  using FlatObject::setActualBufferAddress;
  using FlatObject::setFutureBufferAddress;
 
 
  void construct(int32_t numberOfKnots, const float knots[], int32_t numberOfAxisBins);
 
  void constructRegular(int32_t numberOfKnotsU);
 
 
  template <typename T>
  GPUd() void correctEdges(T* data) const;
 
  template <typename T>
  GPUd() static T getSpline(const IrregularSpline1D::Knot& knot1, T f0, T f1, T f2, T f3, float u);
 
  template <typename T>
  GPUd() T getSpline(const T correctedData[], float u) const;
 
  GPUd() int32_t getNumberOfKnots() const { return mNumberOfKnots; }
 
  GPUd() int32_t getKnotIndex(float u) const;
 
  GPUd() const IrregularSpline1D::Knot& getKnot(int32_t i) const { return getKnots()[i]; }
 
  GPUd() const IrregularSpline1D::Knot* getKnots() const { return reinterpret_cast<const IrregularSpline1D::Knot*>(mFlatBufferPtr); }
 
  static constexpr size_t getClassAlignmentBytes() { return 8; }
 
  static constexpr size_t getBufferAlignmentBytes() { return 8; }
 
  static constexpr size_t getDataAlignmentBytes() { return 8; }
 
 
  GPUd() const int32_t* getBin2KnotMap() const { return reinterpret_cast<const int32_t*>(mFlatBufferPtr + mBin2KnotMapOffset); }
 
  int32_t getNumberOfAxisBins() const { return mNumberOfAxisBins; }
 
  GPUd() static void getEdgeCorrectionCoefficients(double u0, double u1, double u2, double u3, double& c0, double& c1, double& c2, double& c3);
 
  void print() const;
 
 private:
  IrregularSpline1D::Knot* getKnotsNonConst() { return reinterpret_cast<IrregularSpline1D::Knot*>(mFlatBufferPtr); }
 
  int32_t* getBin2KnotMapNonConst() { return reinterpret_cast<int32_t*>(mFlatBufferPtr + mBin2KnotMapOffset); }
 
 
  int32_t mNumberOfKnots;      
  int32_t mNumberOfAxisBins;   
  uint32_t mBin2KnotMapOffset; 
 
  ClassDefNV(IrregularSpline1D, 1);
};
 
 
template <typename T>
GPUdi() T IrregularSpline1D::getSpline(const IrregularSpline1D::Knot& knot1, T f0, T f1, T f2, T f3, float u)
{
 
  f0 -= f1;
  f2 -= f1;
  f3 -= f1;
 
  T x = T((u - knot1.u) * knot1.scale);              // scaled u
  T z1 = T(f0 * knot1.scaleL0 + f2 * knot1.scaleL2); // scaled u derivative at the knot 1
  T z2 = T(f2 * knot1.scaleR2 + f3 * knot1.scaleR3); // scaled u derivative at the knot 2
 
  T x2 = x * x;
 
  // f(x) = ax^3 + bx^2 + cx + d
  //
  // f(0) = f1
  // f(1) = f2
  // f'(0) = z1
  // f'(1) = z2
  //
 
  // T d = f1;
  // T c = z1;
  T a = -f2 - f2 + z1 + z2;
  T b = f2 - z1 - a;
  return a * x * x2 + b * x2 + z1 * x + f1;
}
 
template <typename T>
GPUdi() T IrregularSpline1D::getSpline(const T correctedData[], float u) const
{
  int32_t iknot = getKnotIndex(u);
  const IrregularSpline1D::Knot& knot = getKnot(iknot);
  const T* f = correctedData + iknot - 1;
  return getSpline(knot, f[0], f[1], f[2], f[3], u);
}
 
GPUdi() int32_t IrregularSpline1D::getKnotIndex(float u) const
{
  int32_t ibin = (int32_t)(u * mNumberOfAxisBins);
  if (ibin < 0) {
    ibin = 0;
  }
  if (ibin > mNumberOfAxisBins - 1) {
    ibin = mNumberOfAxisBins - 1;
  }
  return getBin2KnotMap()[ibin];
}
 
GPUdi() void IrregularSpline1D::getEdgeCorrectionCoefficients(double u0, double u1, double u2, double u3, double& c0, double& c1, double& c2, double& c3)
{
 
  double du = u2 - u1;
  double x0 = (u0 - u1) / du;
  // double x1 = ( u1 - u1 )/du = 0
  // double x2 = ( u2 - u1 )/du = 1
  double x3 = (u3 - u1) / du;
 
  double cL0 = -1. / (x0 * (x0 - 1.));
  double cL2 = x0 / (x0 - 1.);
  double cR2 = (x3 - 2.) / (x3 - 1.);
  double cR3 = 1. / (x3 * (x3 - 1.));
 
  // Fi = fi - f1
 
  // z1 = F0corr*cL0 + F2*cL2; // scaled u derivative at the knot 1
  // z2 = F2*cR2 + F3*cR3; // scaled u derivative at the knot 2
 
  // a =  -2F2 +   z1 + z2
  // b =   3F2 - 2*z1 - z2 ;
  // c =           z1
  // d = 0
 
  // F(x) = ax^3 + bx^2 + cx + d
  // F(x) = (-2*x^3+3*x^2)*F2 + (x^3-2*x^2+x)*z1 + (x^3-x^2)*z2
  // F(x) = (-2*x^3+3*x^2)*F2 + (x^3-2*x^2+x)*z1 + (x^3-x^2)*(cR2*F2 + cR3*F3)
  // F(x) = (-2*x^3+3*x^2+(x^3-x^2)*cR2)*F2 + (x^3-2*x^2+x)*z1 + (x^3-x^2)*cR3*F3
  // F(x) = x^2(-2*x+3+(x-1)*cR2)*F2 + x*(x-1)^2*z1 + x^2(x-1)*cR3*F3
 
  // z1*x*(x-1)^2 = F(x) - x^2(-2*x+3+(x-1)*cR2)*F2 - x^2(x-1)*cR3*F3
  // z1*x0*(x0-1)^2 = F0 - x0^2(-2*x0+3+(x0-1)*cR2)*F2 - x0^2(x0-1)*cR3*F3
 
  double x01 = x0 - 1.;
 
  // coeff. for z1 at F0, F1, F2, F3:
 
  c0 = 1. / (x0 * x01 * x01);
  c1 = 0;
  c2 = -(-2 * x0 + 3. + x01 * cR2) * x0 / (x01 * x01);
  c3 = -x0 * cR3 / x01;
 
  // F0corr = (z1 - F2*cL2)/cL0;
 
  c0 = c0 / cL0;
  c1 = c1 / cL0;
  c2 = (c2 - cL2) / cL0;
  c3 = c3 / cL0;
 
  // F0corr = c0*F0 + c1*F1 + C2*F2 + c3*F3;
  // f0corr-f1 = c0*(f0-f1) + c1*(f1-f1) + c2*(f2-f1) + c3*(f3-f1);
  // f0corr-f1 = c0*(f0-f1) + c2*(f2-f1) + c3*(f3-f1);
  // f0corr = c0*f0 + c2*f2 + c3*f3 + (1-c0-c2-c3)*f1;
 
  c1 = 1. - c0 - c2 - c3;
}
 
template <typename T>
GPUdi() void IrregularSpline1D::correctEdges(T* data) const
{
  const IrregularSpline1D::Knot* s = getKnots();
  double c0, c1, c2, c3;
  getEdgeCorrectionCoefficients(s[0].u, s[1].u, s[2].u, s[3].u, c0, c1, c2, c3);
  data[0] = c0 * data[0] + c1 * data[1] + c2 * data[2] + c3 * data[3];
  int32_t i = mNumberOfKnots - 1;
  getEdgeCorrectionCoefficients(s[i - 0].u, s[i - 1].u, s[i - 2].u, s[i - 3].u, c0, c1, c2, c3);
  data[i] = c0 * data[i - 0] + c1 * data[i - 1] + c2 * data[i - 2] + c3 * data[i - 3];
}
} // namespace gpu
} // namespace o2
 
#endif