Chebyshev3D.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_MATHUTILS_CHEBYSHEV3D_H_
#define ALICEO2_MATHUTILS_CHEBYSHEV3D_H_
 
#include <TNamed.h>                    // for TNamed
#include <TObjArray.h>                 // for TObjArray
#include <cstdio>                      // for FILE, stdout
#include "MathUtils/Chebyshev3DCalc.h" // for Chebyshev3DCalc, etc
#include "Rtypes.h"                    // for Float_t, Int_t, Double_t, Bool_t, etc
#include "TString.h"                   // for TString
 
class TH1;         // lines 15-15
class TMethodCall; // lines 16-16
 
namespace o2
{
namespace math_utils
{
 
class Chebyshev3D : public TNamed
{
 public:
  Chebyshev3D();
 
  Chebyshev3D(const Chebyshev3D& src);
 
  Chebyshev3D(const char* inpFile);
 
  Chebyshev3D(FILE* stream);
 
#ifdef _INC_CREATION_Chebyshev3D_
  Chebyshev3D(const char* funName, Int_t dimOut, const Float_t* bmin, const Float_t* bmax, const Int_t* npoints,
              Float_t prec = 1E-6, const Float_t* precD = nullptr);
  Chebyshev3D(void (*ptr)(float*, float*), Int_t dimOut, const Float_t* bmin, const Float_t* bmax, const Int_t* npoints,
              Float_t prec = 1E-6, const Float_t* precD = nullptr);
  Chebyshev3D(void (*ptr)(float*, float*), int dimOut, const Float_t* bmin, const Float_t* bmax,
              const Int_t* npX, const Int_t* npY, const Int_t* npZ,
              Float_t prec = 1E-6, const Float_t* precD = nullptr);
  Chebyshev3D(void (*ptr)(float*, float*), int DimOut, const Float_t* bmin, const Float_t* bmax, Float_t prec = 1E-6,
              Bool_t run = kTRUE, const Float_t* precD = nullptr);
#endif
 
  ~Chebyshev3D() override
  {
    Clear();
  }
 
  Chebyshev3D& operator=(const Chebyshev3D& rhs);
 
  void Eval(const Float_t* par, Float_t* res);
 
  Float_t Eval(const Float_t* par, int idim);
 
  void Eval(const Double_t* par, Double_t* res);
 
  Double_t Eval(const Double_t* par, int idim);
 
  void evaluateDerivative(int dimd, const Float_t* par, Float_t* res);
 
  void evaluateDerivative2(int dimd1, int dimd2, const Float_t* par, Float_t* res);
 
  Float_t evaluateDerivative(int dimd, const Float_t* par, int idim);
 
  Float_t evaluateDerivative2(int dimd1, int dimd2, const Float_t* par, int idim);
 
  void evaluateDerivative3D(const Float_t* par, Float_t dbdr[3][3]);
 
  void evaluateDerivative3D2(const Float_t* par, Float_t dbdrdr[3][3][3]);
 
  void Print(const Option_t* opt = "") const override;
 
  Bool_t isInside(const Float_t* par) const;
 
  Bool_t isInside(const Double_t* par) const;
 
  Chebyshev3DCalc* getChebyshevCalc(int i) const
  {
    return (Chebyshev3DCalc*)mChebyshevParameter.UncheckedAt(i);
  }
 
  Float_t getBoundMin(int i) const
  {
    return mMinBoundaries[i];
  }
 
  Float_t getBoundMax(int i) const
  {
    return mMaxBoundaries[i];
  }
 
  Float_t* getBoundMin() const
  {
    return (float*)mMinBoundaries;
  }
 
  Float_t* getBoundMax() const
  {
    return (float*)mMaxBoundaries;
  }
 
  Float_t* getBoundaryMappingScale() const
  {
    return (float*)mBoundaryMappingScale;
  }
 
  Float_t* getBoundaryMappingOffset() const
  {
    return (float*)mBoundaryMappingOffset;
  }
 
  Float_t getPrecision() const
  {
    return mPrecision;
  }
 
  void shiftBound(int id, float dif);
 
  void loadData(const char* inpFile);
 
  void loadData(FILE* stream);
 
#ifdef _INC_CREATION_Chebyshev3D_
  void invertSign();
  int* getNcNeeded(float xyz[3], int dimVar, float mn, float mx, float prec, Int_t npCheck = 30);
  void estimateNumberOfPoints(float prec, int gridBC[3][3], Int_t npd1 = 30, Int_t npd2 = 30, Int_t npd3 = 30);
  void saveData(const char* outfile, Bool_t append = kFALSE) const;
  void saveData(FILE* stream = stdout) const;
 
  void setuserFunction(const char* name);
  void setuserFunction(void (*ptr)(float*, float*));
  void evaluateUserFunction(const Float_t* x, Float_t* res);
  TH1* TestRMS(int idim, int npoints = 1000, TH1* histo = nullptr);
  static Int_t calculateChebyshevCoefficients(const Float_t* funval, int np, Float_t* outCoefs, Float_t prec = -1);
#endif
 
 protected:
  void Clear(const Option_t* option = "") override;
 
  void setDimOut(const int d, const float* prec = nullptr);
 
  void prepareBoundaries(const Float_t* bmin, const Float_t* bmax);
 
#ifdef _INC_CREATION_Chebyshev3D_
  void evaluateUserFunction();
  void defineGrid(const Int_t* npoints);
  Int_t chebyshevFit(); // fit all output dimensions
  Int_t chebyshevFit(int dmOut);
  void setPrecision(float prec)
  {
    mPrecision = prec;
  }
#endif
 
  Float_t mapToInternal(Float_t x, Int_t d) const; // map x to [-1:1]
  Float_t mapToExternal(Float_t x, Int_t d) const
  {
    return x / mBoundaryMappingScale[d] + mBoundaryMappingOffset[d];
  }                                                  // map from [-1:1] to x
  Double_t mapToInternal(Double_t x, Int_t d) const; // map x to [-1:1]
  Double_t mapToExternal(Double_t x, Int_t d) const
  {
    return x / mBoundaryMappingScale[d] + mBoundaryMappingOffset[d];
  } // map from [-1:1] to x
 
 private:
  Int_t mOutputArrayDimension;       
  Float_t mPrecision;                
  Float_t mMinBoundaries[3];         
  Float_t mMaxBoundaries[3];         
  Float_t mBoundaryMappingScale[3];  
  Float_t mBoundaryMappingOffset[3]; 
  TObjArray mChebyshevParameter;     
 
  Int_t mMaxCoefficients;               
  Int_t mNumberOfPoints[3];             
  Float_t mTemporaryCoefficient[3];     
  Float_t* mTemporaryUserResults;       
  Float_t* mTemporaryChebyshevGrid;     
  Int_t mTemporaryChebyshevGridOffs[3]; 
  TString mUserFunctionName;            
  TMethodCall* mUserMacro;              
 
  static const Float_t sMinimumPrecision; 
 
  ClassDefOverride(o2::math_utils::Chebyshev3D,
                   2) // Chebyshev parametrization for 3D->N function
};
 
inline Bool_t Chebyshev3D::isInside(const Float_t* par) const
{
  for (int i = 3; i--;) {
    if (mMinBoundaries[i] > par[i] || par[i] > mMaxBoundaries[i]) {
      return kFALSE;
    }
  }
  return kTRUE;
}
 
inline Bool_t Chebyshev3D::isInside(const Double_t* par) const
{
  for (int i = 3; i--;) {
    if (mMinBoundaries[i] > par[i] || par[i] > mMaxBoundaries[i]) {
      return kFALSE;
    }
  }
  return kTRUE;
}
 
inline void Chebyshev3D::Eval(const Float_t* par, Float_t* res)
{
  Float_t coef[3]; // mapped arguments, local so that evaluation is reentrant
  for (int i = 3; i--;) {
    coef[i] = mapToInternal(par[i], i);
  }
  for (int i = mOutputArrayDimension; i--;) {
    res[i] = getChebyshevCalc(i)->Eval(coef);
  }
}
 
inline void Chebyshev3D::Eval(const Double_t* par, Double_t* res)
{
  Float_t coef[3]; // mapped arguments, local so that evaluation is reentrant
  for (int i = 3; i--;) {
    coef[i] = mapToInternal(par[i], i);
  }
  for (int i = mOutputArrayDimension; i--;) {
    res[i] = getChebyshevCalc(i)->Eval(coef);
  }
}
 
inline Double_t Chebyshev3D::Eval(const Double_t* par, int idim)
{
  Float_t coef[3]; // mapped arguments, local so that evaluation is reentrant
  for (int i = 3; i--;) {
    coef[i] = mapToInternal(par[i], i);
  }
  return getChebyshevCalc(idim)->Eval(coef);
}
 
inline Float_t Chebyshev3D::Eval(const Float_t* par, int idim)
{
  Float_t coef[3]; // mapped arguments, local so that evaluation is reentrant
  for (int i = 3; i--;) {
    coef[i] = mapToInternal(par[i], i);
  }
  return getChebyshevCalc(idim)->Eval(coef);
}
 
inline void Chebyshev3D::evaluateDerivative3D(const Float_t* par, Float_t dbdr[3][3])
{
  Float_t coef[3]; // mapped arguments, local so that evaluation is reentrant
  for (int i = 3; i--;) {
    coef[i] = mapToInternal(par[i], i);
  }
  for (int ib = 3; ib--;) {
    for (int id = 3; id--;) {
      dbdr[ib][id] = getChebyshevCalc(ib)->evaluateDerivative(id, coef) * mBoundaryMappingScale[id];
    }
  }
}
 
inline void Chebyshev3D::evaluateDerivative3D2(const Float_t* par, Float_t dbdrdr[3][3][3])
{
  Float_t coef[3]; // mapped arguments, local so that evaluation is reentrant
  for (int i = 3; i--;) {
    coef[i] = mapToInternal(par[i], i);
  }
  for (int ib = 3; ib--;) {
    for (int id = 3; id--;) {
      for (int id1 = 3; id1--;) {
        dbdrdr[ib][id][id1] = getChebyshevCalc(ib)->evaluateDerivative2(id, id1, coef) *
                              mBoundaryMappingScale[id] * mBoundaryMappingScale[id1];
      }
    }
  }
}
 
// Evaluates Chebyshev parameterization derivative for 3d->DimOut function
inline void Chebyshev3D::evaluateDerivative(int dimd, const Float_t* par, Float_t* res)
{
  Float_t coef[3]; // mapped arguments, local so that evaluation is reentrant
  for (int i = 3; i--;) {
    coef[i] = mapToInternal(par[i], i);
  }
  for (int i = mOutputArrayDimension; i--;) {
    res[i] = getChebyshevCalc(i)->evaluateDerivative(dimd, coef) * mBoundaryMappingScale[dimd];
  };
}
 
// Evaluates Chebyshev parameterization 2nd derivative over dimd1 and dimd2 dimensions for 3d->DimOut function
inline void Chebyshev3D::evaluateDerivative2(int dimd1, int dimd2, const Float_t* par, Float_t* res)
{
  Float_t coef[3]; // mapped arguments, local so that evaluation is reentrant
  for (int i = 3; i--;) {
    coef[i] = mapToInternal(par[i], i);
  }
  for (int i = mOutputArrayDimension; i--;) {
    res[i] = getChebyshevCalc(i)->evaluateDerivative2(dimd1, dimd2, coef) *
             mBoundaryMappingScale[dimd1] * mBoundaryMappingScale[dimd2];
  }
}
 
inline Float_t Chebyshev3D::evaluateDerivative(int dimd, const Float_t* par, int idim)
{
  Float_t coef[3]; // mapped arguments, local so that evaluation is reentrant
  for (int i = 3; i--;) {
    coef[i] = mapToInternal(par[i], i);
  }
  return getChebyshevCalc(idim)->evaluateDerivative(dimd, coef) * mBoundaryMappingScale[dimd];
}
 
inline Float_t Chebyshev3D::evaluateDerivative2(int dimd1, int dimd2, const Float_t* par, int idim)
{
  Float_t coef[3]; // mapped arguments, local so that evaluation is reentrant
  for (int i = 3; i--;) {
    coef[i] = mapToInternal(par[i], i);
  }
  return getChebyshevCalc(idim)->evaluateDerivative2(dimd1, dimd2, coef) *
         mBoundaryMappingScale[dimd1] * mBoundaryMappingScale[dimd2];
}
 
inline Float_t Chebyshev3D::mapToInternal(Float_t x, Int_t d) const
{
#ifdef _BRING_TO_BOUNDARY_
  T res = (x - mBoundaryMappingOffset[d]) * mBoundaryMappingScale[d];
  if (res < -1) {
    return -1;
  }
  if (res > 1) {
    return 1;
  }
  return res;
#else
  return (x - mBoundaryMappingOffset[d]) * mBoundaryMappingScale[d];
#endif
}
 
inline Double_t Chebyshev3D::mapToInternal(Double_t x, Int_t d) const
{
#ifdef _BRING_TO_BOUNDARY_
  T res = (x - mBoundaryMappingOffset[d]) * mBoundaryMappingScale[d];
  if (res < -1) {
    return -1;
  }
  if (res > 1) {
    return 1;
  }
  return res;
#else
  return (x - mBoundaryMappingOffset[d]) * mBoundaryMappingScale[d];
#endif
}
} // namespace math_utils
} // namespace o2
 
#endif