STEER/AliCheb3D.h

   1 // Author: ruben.shahoyan@cern.ch   09/09/2006
   2
   3 ////////////////////////////////////////////////////////////////////////////////
   4 //                                                                            //
   5 // AliCheb3D produces the interpolation of the user 3D->NDimOut arbitrary     //
   6 // function supplied in "void (*fcn)(float* inp,float* out)" format           //
   7 // either in a separate macro file or as a function pointer.                  //
   8 // Only coefficients needed to guarantee the requested precision are kept.    //
   9 //                                                                            //
  10 // The user-callable methods are:                                             //
  11 // To create the interpolation use:                                           //
  12 // AliCheb3D(const char* funName,  // name of the file with user function     //
  13 //          or                                                                //
  14 // AliCheb3D(void (*ptr)(float*,float*),// pointer on the  user function      //
  15 //        Int_t     DimOut,     // dimensionality of the function's output    //
  16 //        Float_t  *bmin,       // lower 3D bounds of interpolation domain    //
  17 //        Float_t  *bmax,       // upper 3D bounds of interpolation domain    //
  18 //        Int_t    *npoints,    // number of points in each of 3 input        //
  19 //                              // dimension, defining the interpolation grid //
  20 //        Float_t   prec=1E-6); // requested max.absolute difference between  //
  21 //                              // the interpolation and any point on grid    //
  22 //                                                                            //
  23 // To test obtained parameterization use the method                           //
  24 // TH1* TestRMS(int idim,int npoints = 1000,TH1* histo=0);                    //
  25 // it will compare the user output of the user function and interpolation     //
  26 // for idim-th output dimension and fill the difference in the supplied       //
  27 // histogram. If no histogram is supplied, it will be created.                //
  28 //                                                                            //
  29 // To save the interpolation data:                                            //
  30 // SaveData(const char* filename, Bool_t append )                             //
  31 // write text file with data. If append is kTRUE and the output file already  //
  32 // exists, data will be added in the end of the file.                         //
  33 // Alternatively, SaveData(FILE* stream) will write the data to               //
  34 // already existing stream.                                                   //
  35 //                                                                            //
  36 // To read back already stored interpolation use either the constructor       //
  37 // AliCheb3D(const char* inpFile);                                            //
  38 // or the default constructor AliCheb3D() followed by                         //
  39 // AliCheb3D::LoadData(const char* inpFile);                                  //
  40 //                                                                            //
  41 // To compute the interpolation use Eval(float* par,float *res) method, with  //
  42 // par being 3D vector of arguments (inside the validity region) and res is   //
  43 // the array of DimOut elements for the output.                               //
  44 //                                                                            //
  45 // If only one component (say, idim-th) of the output is needed, use faster   //
  46 // Float_t Eval(Float_t *par,int idim) method.                                //
  47 //                                                                            //
  48 // void Print(option="") will print the name, the ranges of validity and      //
  49 // the absolute precision of the parameterization. Option "l" will also print //
  50 // the information about the number of coefficients for each output           //
  51 // dimension.                                                                 //
  52 //                                                                            //
  53 // NOTE: during the evaluation no check is done for parameter vector being    //
  54 // outside the interpolation region. If there is such a risk, use             //
  55 // Bool_t IsInside(float *par) method. Chebyshev parameterization is not      //
  56 // good for extrapolation!                                                    //
  57 //                                                                            //
  58 // For the properties of Chebyshev parameterization see:                      //
  59 // H.Wind, CERN EP Internal Report, 81-12/Rev.                                //
  60 //                                                                            //
  61 ////////////////////////////////////////////////////////////////////////////////
  62
  63
  64 #ifndef _ALICHEB3D_
  65 #define _ALICHEB3D_
  66 #include <stdio.h>
  67 #include <TNamed.h>
  68 #include <TMethodCall.h>
  69 #include <TMath.h>
  70 #include <TH1.h>
  71 #include <TObjArray.h>
  72 #include "AliCheb3DCalc.h"
  73
  74 class TString;
  75 class TSystem;
  76 class TRandom;
  77
  78
  79 class AliCheb3D: public TNamed
  80 {
  81  public:
  82   AliCheb3D();
  83   AliCheb3D(const AliCheb3D& src);
  84   AliCheb3D(const char* inpFile);
  85   AliCheb3D(FILE* stream);
  86   //
  87 #ifdef _INC_CREATION_ALICHEB3D_
  88   AliCheb3D(const char* funName, Int_t DimOut, Float_t  *bmin,Float_t  *bmax, Int_t *npoints, Float_t  prec=1E-6);
  89   AliCheb3D(void (*ptr)(float*,float*), Int_t DimOut, Float_t  *bmin,Float_t  *bmax, Int_t *npoints, Float_t  prec=1E-6);
  90   AliCheb3D(void (*ptr)(float*,float*), int DimOut, Float_t  *bmin,Float_t  *bmax, Int_t *npX,Int_t *npY,Int_t *npZ, Float_t prec=1E-6);
  91   AliCheb3D(void (*ptr)(float*,float*), int DimOut, Float_t  *bmin,Float_t  *bmax, Float_t prec=1E-6);
  92 #endif
  93   //
  94   ~AliCheb3D()                                                                 {Clear();}
  95   //
  96   AliCheb3D&   operator=(const AliCheb3D& rhs);
  97   void         Eval(Float_t  *par,Float_t  *res);
  98   Float_t      Eval(Float_t  *par,int idim);
  99   //
 100   void         EvalDeriv(int dimd, Float_t  *par,Float_t  *res);
 101   void         EvalDeriv2(int dimd1, int dimd2,Float_t  *par,Float_t  *res);
 102   Float_t      EvalDeriv(int dimd,Float_t  *par, int idim);
 103   Float_t      EvalDeriv2(int dimd1,int dimd2, Float_t  *par, int idim);
 104   void         EvalDeriv3D(Float_t *par, Float_t dbdr[3][3]);
 105   void         EvalDeriv3D2(Float_t *par, Float_t dbdrdr[3][3][3]);
 106   void         Print(Option_t* opt="")                                   const;
 107   Bool_t       IsInside(Float_t  *par)                                   const;
 108   Bool_t       IsInside(Double_t  *par)                                  const;
 109   AliCheb3DCalc*  GetChebCalc(int i)                                     const {return (AliCheb3DCalc*)fChebCalc.UncheckedAt(i);}
 110   Float_t      GetBoundMin(int i)                                        const {return fBMin[i];}
 111   Float_t      GetBoundMax(int i)                                        const {return fBMax[i];}
 112   Float_t*     GetBoundMin()                                             const {return (float*)fBMin;}
 113   Float_t*     GetBoundMax()                                             const {return (float*)fBMax;}
 114   Float_t      GetPrecision()                                            const {return fPrec;}
 115   void         ShiftBound(int id,float dif);
 116   //
 117   void         LoadData(const char* inpFile);
 118   void         LoadData(FILE* stream);
 119   //
 120 #ifdef _INC_CREATION_ALICHEB3D_
 121   int*         GetNCNeeded(float xyz[3],int DimVar, float mn,float mx, float prec);
 122   void         EstimateNPoints(float Prec, int gridBC[3][3]);
 123   void         SaveData(const char* outfile,Bool_t append=kFALSE)        const;
 124   void         SaveData(FILE* stream=stdout)                             const;
 125   //
 126   void         SetUsrFunction(const char* name);
 127   void         SetUsrFunction(void (*ptr)(float*,float*));
 128   void         EvalUsrFunction(Float_t  *x, Float_t  *res);
 129   TH1*         TestRMS(int idim,int npoints = 1000,TH1* histo=0);
 130   static Int_t CalcChebCoefs(Float_t  *funval,int np, Float_t  *outCoefs, Float_t  prec=-1);
 131 #endif
 132   //
 133  protected:
 134   void         Clear(Option_t* option = "");
 135   void         SetDimOut(int d);
 136   void         PrepareBoundaries(Float_t  *bmin,Float_t  *bmax);
 137   //
 138 #ifdef _INC_CREATION_ALICHEB3D_
 139   void         EvalUsrFunction();
 140   void         DefineGrid(Int_t* npoints);
 141   Int_t        ChebFit();                                                                 // fit all output dimensions
 142   Int_t        ChebFit(int dmOut);
 143 #endif
 144   //
 145   Float_t      MapToInternal(Float_t  x,Int_t d) const; // map x to [-1:1]
 146   Float_t      MapToExternal(Float_t  x,Int_t d) const {return x/fBScale[d]+fBOffset[d];}   // map from [-1:1] to x
 147   //
 148  protected:
 149   Int_t        fDimOut;            // dimension of the ouput array
 150   Float_t      fPrec;              // requested precision
 151   Float_t      fBMin[3];           // min boundaries in each dimension
 152   Float_t      fBMax[3];           // max boundaries in each dimension
 153   Float_t      fBScale[3];         // scale for boundary mapping to [-1:1] interval
 154   Float_t      fBOffset[3];        // offset for boundary mapping to [-1:1] interval
 155   TObjArray    fChebCalc;          // Chebyshev parameterization for each output dimension
 156   //
 157   Int_t        fMaxCoefs;          //! max possible number of coefs per parameterization
 158   Int_t        fNPoints[3];        //! number of used points in each dimension
 159   Float_t      fArgsTmp[3];        //! temporary vector for coefs caluclation
 160   Float_t      fBuff[6];           //! buffer for coordinate transformations
 161   Float_t *    fResTmp;            //! temporary vector for results of user function caluclation
 162   Float_t *    fGrid;              //! temporary buffer for Chebyshef roots grid
 163   Int_t        fGridOffs[3];       //! start of grid for each dimension
 164   TString      fUsrFunName;        //! name of user macro containing the function of  "void (*fcn)(float*,float*)" format
 165   TMethodCall* fUsrMacro;          //! Pointer to MethodCall for function from user macro
 166   //
 167   ClassDef(AliCheb3D,1)  // Chebyshev parametrization for 3D->N function
 168 };
 169
 170 //__________________________________________________________________________________________
 171 inline Bool_t  AliCheb3D::IsInside(Float_t  *par) const
 172 {
 173   // check if the point is inside of the fitted box
 174   const float kTol = 1.e-4;
 175   for (int i=3;i--;) if (fBMin[i]-par[i]>kTol || par[i]-fBMax[i]>kTol) return kFALSE;
 176   //if(par[i]<fBMin[i]||par[i]>fBMax[i]) return kFALSE;
 177   return kTRUE;
 178 }
 179
 180 //__________________________________________________________________________________________
 181 inline Bool_t  AliCheb3D::IsInside(Double_t  *par) const
 182 {
 183   // check if the point is inside of the fitted box
 184   const float kTol = 1.e-4;
 185   for (int i=3;i--;) if (fBMin[i]-par[i]>kTol || par[i]-fBMax[i]>kTol) return kFALSE;
 186   //if(par[i]<fBMin[i]||par[i]>fBMax[i]) return kFALSE;
 187   return kTRUE;
 188 }
 189
 190 //__________________________________________________________________________________________
 191 inline void AliCheb3D::Eval(Float_t  *par, Float_t  *res)
 192 {
 193   // evaluate Chebyshev parameterization for 3d->DimOut function
 194   for (int i=3;i--;) fArgsTmp[i] = MapToInternal(par[i],i);
 195   for (int i=fDimOut;i--;) res[i] = GetChebCalc(i)->Eval(fArgsTmp);
 196   //
 197 }
 198
 199 //__________________________________________________________________________________________
 200 inline Float_t AliCheb3D::Eval(Float_t  *par, int idim)
 201 {
 202   // evaluate Chebyshev parameterization for idim-th output dimension of 3d->DimOut function
 203   for (int i=3;i--;) fArgsTmp[i] = MapToInternal(par[i],i);
 204   return GetChebCalc(idim)->Eval(fArgsTmp);
 205   //
 206 }
 207
 208 //__________________________________________________________________________________________
 209 inline void AliCheb3D::EvalDeriv3D(Float_t *par, Float_t dbdr[3][3])
 210 {
 211   // return gradient matrix
 212   for (int i=3;i--;) fArgsTmp[i] = MapToInternal(par[i],i);
 213   for (int ib=3;ib--;) for (int id=3;id--;) dbdr[ib][id] = GetChebCalc(ib)->EvalDeriv(id,fArgsTmp)*fBScale[id];
 214 }
 215
 216 //__________________________________________________________________________________________
 217 inline void AliCheb3D::EvalDeriv3D2(Float_t *par, Float_t dbdrdr[3][3][3])
 218 {
 219   // return gradient matrix
 220   for (int i=3;i--;) fArgsTmp[i] = MapToInternal(par[i],i);
 221   for (int ib=3;ib--;) for (int id=3;id--;)for (int id1=3;id1--;)
 222     dbdrdr[ib][id][id1] = GetChebCalc(ib)->EvalDeriv2(id,id1,fArgsTmp)*fBScale[id]*fBScale[id1];
 223 }
 224
 225 //__________________________________________________________________________________________
 226 inline void AliCheb3D::EvalDeriv(int dimd,Float_t  *par, Float_t  *res)
 227 {
 228   // evaluate Chebyshev parameterization derivative for 3d->DimOut function
 229   for (int i=3;i--;) fArgsTmp[i] = MapToInternal(par[i],i);
 230   for (int i=fDimOut;i--;) res[i] = GetChebCalc(i)->EvalDeriv(dimd,fArgsTmp)*fBScale[dimd];;
 231   //
 232 }
 233
 234 //__________________________________________________________________________________________
 235 inline void AliCheb3D::EvalDeriv2(int dimd1,int dimd2,Float_t  *par, Float_t  *res)
 236 {
 237   // evaluate Chebyshev parameterization 2nd derivative over dimd1 and dimd2 dimensions for 3d->DimOut function
 238   for (int i=3;i--;) fArgsTmp[i] = MapToInternal(par[i],i);
 239   for (int i=fDimOut;i--;) res[i] = GetChebCalc(i)->EvalDeriv2(dimd1,dimd2,fArgsTmp)*fBScale[dimd1]*fBScale[dimd2];
 240   //
 241 }
 242
 243 //__________________________________________________________________________________________
 244 inline Float_t AliCheb3D::EvalDeriv(int dimd,Float_t  *par, int idim)
 245 {
 246   // evaluate Chebyshev parameterization derivative over dimd dimention for idim-th output dimension of 3d->DimOut function
 247   for (int i=3;i--;) fArgsTmp[i] = MapToInternal(par[i],i);
 248   return GetChebCalc(idim)->EvalDeriv(dimd,fArgsTmp)*fBScale[dimd];
 249   //
 250 }
 251
 252 //__________________________________________________________________________________________
 253 inline Float_t AliCheb3D::EvalDeriv2(int dimd1,int dimd2,Float_t  *par, int idim)
 254 {
 255   // evaluate Chebyshev parameterization 2ns derivative over dimd1 and dimd2 dimensions for idim-th output dimension of 3d->DimOut function
 256   for (int i=3;i--;) fArgsTmp[i] = MapToInternal(par[i],i);
 257   return GetChebCalc(idim)->EvalDeriv2(dimd1,dimd2,fArgsTmp)*fBScale[dimd1]*fBScale[dimd2];
 258   //
 259 }
 260
 261 //__________________________________________________________________________________________
 262 inline Float_t AliCheb3D::MapToInternal(Float_t  x,Int_t d) const
 263 {
 264   // map x to [-1:1]
 265 #ifdef _BRING_TO_BOUNDARY_
 266   float res = (x-fBOffset[d])*fBScale[d];
 267   if (res<-1) return -1;
 268   if (res> 1) return 1;
 269   return res;
 270 #else
 271   return (x-fBOffset[d])*fBScale[d];
 272 #endif
 273 }
 274
 275 #endif