1 // Author: ruben.shahoyan@cern.ch 09/09/2006
3 ////////////////////////////////////////////////////////////////////////////////
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. //
10 // The user-callable methods are: //
11 // To create the interpolation use: //
12 // AliCheb3D(const char* funName, // name of the file with user function //
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 //
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. //
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. //
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); //
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. //
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. //
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 //
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! //
58 // For the properties of Chebyshev parameterization see: //
59 // H.Wind, CERN EP Internal Report, 81-12/Rev. //
61 ////////////////////////////////////////////////////////////////////////////////
68 #include <TMethodCall.h>
71 #include <TObjArray.h>
72 #include "AliCheb3DCalc.h"
79 class AliCheb3D: public TNamed
83 AliCheb3D(const AliCheb3D& src);
84 AliCheb3D(const char* inpFile);
85 AliCheb3D(FILE* stream);
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);
94 ~AliCheb3D() {Clear();}
96 AliCheb3D& operator=(const AliCheb3D& rhs);
97 void Eval(Float_t *par,Float_t *res);
98 Float_t Eval(Float_t *par,int idim);
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);
117 void LoadData(const char* inpFile);
118 void LoadData(FILE* stream);
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;
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);
134 void Clear(Option_t* option = "");
135 void SetDimOut(int d);
136 void PrepareBoundaries(Float_t *bmin,Float_t *bmax);
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);
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
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
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
167 ClassDef(AliCheb3D,1) // Chebyshev parametrization for 3D->N function
170 //__________________________________________________________________________________________
171 inline Bool_t AliCheb3D::IsInside(Float_t *par) const
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;
180 //__________________________________________________________________________________________
181 inline Bool_t AliCheb3D::IsInside(Double_t *par) const
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;
190 //__________________________________________________________________________________________
191 inline void AliCheb3D::Eval(Float_t *par, Float_t *res)
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);
199 //__________________________________________________________________________________________
200 inline Float_t AliCheb3D::Eval(Float_t *par, int idim)
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);
208 //__________________________________________________________________________________________
209 inline void AliCheb3D::EvalDeriv3D(Float_t *par, Float_t dbdr[3][3])
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];
216 //__________________________________________________________________________________________
217 inline void AliCheb3D::EvalDeriv3D2(Float_t *par, Float_t dbdrdr[3][3][3])
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];
225 //__________________________________________________________________________________________
226 inline void AliCheb3D::EvalDeriv(int dimd,Float_t *par, Float_t *res)
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];;
234 //__________________________________________________________________________________________
235 inline void AliCheb3D::EvalDeriv2(int dimd1,int dimd2,Float_t *par, Float_t *res)
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];
243 //__________________________________________________________________________________________
244 inline Float_t AliCheb3D::EvalDeriv(int dimd,Float_t *par, int idim)
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];
252 //__________________________________________________________________________________________
253 inline Float_t AliCheb3D::EvalDeriv2(int dimd1,int dimd2,Float_t *par, int idim)
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];
261 //__________________________________________________________________________________________
262 inline Float_t AliCheb3D::MapToInternal(Float_t x,Int_t d) const
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;
271 return (x-fBOffset[d])*fBScale[d];