Line data Source code
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_H
65 : #define ALICHEB3D_H
66 :
67 : #include <TNamed.h>
68 : #include <TObjArray.h>
69 : #include "AliCheb3DCalc.h"
70 :
71 : class TString;
72 : class TSystem;
73 : class TRandom;
74 : class TH1;
75 : class TMethodCall;
76 : class TRandom;
77 : class TROOT;
78 : class stdio;
79 :
80 :
81 :
82 : class AliCheb3D: public TNamed
83 : {
84 : public:
85 : AliCheb3D();
86 : AliCheb3D(const AliCheb3D& src);
87 : AliCheb3D(const char* inpFile);
88 : AliCheb3D(FILE* stream);
89 : //
90 : #ifdef _INC_CREATION_ALICHEB3D_
91 : AliCheb3D(const char* funName, Int_t DimOut, const Float_t *bmin, const Float_t *bmax, Int_t *npoints, Float_t prec=1E-6, const Float_t* precD=0);
92 : AliCheb3D(void (*ptr)(float*,float*), Int_t DimOut, Float_t *bmin,Float_t *bmax, Int_t *npoints, Float_t prec=1E-6, const Float_t* precD=0);
93 : 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, const Float_t* precD=0);
94 : AliCheb3D(void (*ptr)(float*,float*), int DimOut, Float_t *bmin,Float_t *bmax, Float_t prec=1E-6, Bool_t run=kTRUE, const Float_t* precD=0);
95 : #endif
96 : //
97 160440 : ~AliCheb3D() {Clear();}
98 : //
99 : AliCheb3D& operator=(const AliCheb3D& rhs);
100 : void Eval(const Float_t *par, Float_t *res);
101 : Float_t Eval(const Float_t *par,int idim);
102 : void Eval(const Double_t *par, Double_t *res);
103 : Double_t Eval(const Double_t *par,int idim);
104 : //
105 : void EvalDeriv(int dimd, const Float_t *par, Float_t *res);
106 : void EvalDeriv2(int dimd1, int dimd2, const Float_t *par,Float_t *res);
107 : Float_t EvalDeriv(int dimd, const Float_t *par, int idim);
108 : Float_t EvalDeriv2(int dimd1,int dimd2, const Float_t *par, int idim);
109 : void EvalDeriv3D(const Float_t *par, Float_t dbdr[3][3]);
110 : void EvalDeriv3D2(const Float_t *par, Float_t dbdrdr[3][3][3]);
111 : void Print(const Option_t* opt="") const;
112 : Bool_t IsInside(const Float_t *par) const;
113 : Bool_t IsInside(const Double_t *par) const;
114 : //
115 17637716 : AliCheb3DCalc* GetChebCalc(int i) const {return (AliCheb3DCalc*)fChebCalc.UncheckedAt(i);}
116 0 : Float_t GetBoundMin(int i) const {return fBMin[i];}
117 0 : Float_t GetBoundMax(int i) const {return fBMax[i];}
118 0 : Float_t* GetBoundMin() const {return (float*)fBMin;}
119 0 : Float_t* GetBoundMax() const {return (float*)fBMax;}
120 0 : Float_t GetPrecision() const {return fPrec;}
121 : void ShiftBound(int id,float dif);
122 : //
123 : void LoadData(const char* inpFile);
124 : void LoadData(FILE* stream);
125 : //
126 : #ifdef _INC_CREATION_ALICHEB3D_
127 : void InvertSign();
128 : int* GetNCNeeded(float xyz[3],int DimVar, float mn,float mx, float prec, Int_t npCheck=30);
129 : void EstimateNPoints(float prec, int gridBC[3][3],Int_t npd1=30,Int_t npd2=30,Int_t npd3=30);
130 : void SaveData(const char* outfile,Bool_t append=kFALSE) const;
131 : void SaveData(FILE* stream=stdout) const;
132 : //
133 : void SetUsrFunction(const char* name);
134 : void SetUsrFunction(void (*ptr)(float*,float*));
135 : void EvalUsrFunction(const Float_t *x, Float_t *res);
136 : TH1* TestRMS(int idim,int npoints = 1000,TH1* histo=0);
137 : static Int_t CalcChebCoefs(const Float_t *funval,int np, Float_t *outCoefs, Float_t prec=-1);
138 : #endif
139 : //
140 : protected:
141 : void Clear(const Option_t* option = "");
142 : void SetDimOut(const int d, const float* prec=0);
143 : void PrepareBoundaries(const Float_t *bmin,const Float_t *bmax);
144 : //
145 : #ifdef _INC_CREATION_ALICHEB3D_
146 : void EvalUsrFunction();
147 : void DefineGrid(Int_t* npoints);
148 : Int_t ChebFit(); // fit all output dimensions
149 : Int_t ChebFit(int dmOut);
150 : void SetPrecision(float prec) {fPrec = prec;}
151 : #endif
152 : //
153 : Float_t MapToInternal(Float_t x,Int_t d) const; // map x to [-1:1]
154 0 : Float_t MapToExternal(Float_t x,Int_t d) const {return x/fBScale[d]+fBOffset[d];} // map from [-1:1] to x
155 : Double_t MapToInternal(Double_t x,Int_t d) const; // map x to [-1:1]
156 : Double_t MapToExternal(Double_t x,Int_t d) const {return x/fBScale[d]+fBOffset[d];} // map from [-1:1] to x
157 : //
158 : protected:
159 : Int_t fDimOut; // dimension of the ouput array
160 : Float_t fPrec; // requested precision
161 : Float_t fBMin[3]; // min boundaries in each dimension
162 : Float_t fBMax[3]; // max boundaries in each dimension
163 : Float_t fBScale[3]; // scale for boundary mapping to [-1:1] interval
164 : Float_t fBOffset[3]; // offset for boundary mapping to [-1:1] interval
165 : TObjArray fChebCalc; // Chebyshev parameterization for each output dimension
166 : //
167 : Int_t fMaxCoefs; //! max possible number of coefs per parameterization
168 : Int_t fNPoints[3]; //! number of used points in each dimension
169 : Float_t fArgsTmp[3]; //! temporary vector for coefs caluclation
170 : Float_t * fResTmp; //! temporary vector for results of user function caluclation
171 : Float_t * fGrid; //! temporary buffer for Chebyshef roots grid
172 : Int_t fGridOffs[3]; //! start of grid for each dimension
173 : TString fUsrFunName; //! name of user macro containing the function of "void (*fcn)(float*,float*)" format
174 : TMethodCall* fUsrMacro; //! Pointer to MethodCall for function from user macro
175 : //
176 : static const Float_t fgkMinPrec; // smallest precision
177 : //
178 46026 : ClassDef(AliCheb3D,2) // Chebyshev parametrization for 3D->N function
179 : };
180 :
181 : //__________________________________________________________________________________________
182 : inline Bool_t AliCheb3D::IsInside(const Float_t *par) const
183 : {
184 : // check if the point is inside of the fitted box
185 0 : for (int i=3;i--;) if (fBMin[i]>par[i] || par[i]>fBMax[i]) return kFALSE;
186 0 : return kTRUE;
187 0 : }
188 :
189 : //__________________________________________________________________________________________
190 : inline Bool_t AliCheb3D::IsInside(const Double_t *par) const
191 : {
192 : // check if the point is inside of the fitted box
193 42506027 : for (int i=3;i--;) if (fBMin[i]>par[i] || par[i]>fBMax[i]) return kFALSE;
194 3146283 : return kTRUE;
195 3331565 : }
196 :
197 : //__________________________________________________________________________________________
198 : inline void AliCheb3D::Eval(const Float_t *par, Float_t *res)
199 : {
200 : // evaluate Chebyshev parameterization for 3d->DimOut function
201 0 : for (int i=3;i--;) fArgsTmp[i] = MapToInternal(par[i],i);
202 0 : for (int i=fDimOut;i--;) res[i] = GetChebCalc(i)->Eval(fArgsTmp);
203 : //
204 0 : }
205 : //__________________________________________________________________________________________
206 : inline void AliCheb3D::Eval(const Double_t *par, Double_t *res)
207 : {
208 : // evaluate Chebyshev parameterization for 3d->DimOut function
209 25564302 : for (int i=3;i--;) fArgsTmp[i] = MapToInternal(par[i],i);
210 22723824 : for (int i=fDimOut;i--;) res[i] = GetChebCalc(i)->Eval(fArgsTmp);
211 : //
212 2840478 : }
213 :
214 : //__________________________________________________________________________________________
215 : inline Double_t AliCheb3D::Eval(const Double_t *par, int idim)
216 : {
217 : // evaluate Chebyshev parameterization for idim-th output dimension of 3d->DimOut function
218 2676816 : for (int i=3;i--;) fArgsTmp[i] = MapToInternal(par[i],i);
219 297424 : return GetChebCalc(idim)->Eval(fArgsTmp);
220 : //
221 : }
222 :
223 : //__________________________________________________________________________________________
224 : inline Float_t AliCheb3D::Eval(const Float_t *par, int idim)
225 : {
226 : // evaluate Chebyshev parameterization for idim-th output dimension of 3d->DimOut function
227 0 : for (int i=3;i--;) fArgsTmp[i] = MapToInternal(par[i],i);
228 0 : return GetChebCalc(idim)->Eval(fArgsTmp);
229 : //
230 : }
231 :
232 : //__________________________________________________________________________________________
233 : inline void AliCheb3D::EvalDeriv3D(const Float_t *par, Float_t dbdr[3][3])
234 : {
235 : // return gradient matrix
236 0 : for (int i=3;i--;) fArgsTmp[i] = MapToInternal(par[i],i);
237 0 : for (int ib=3;ib--;) for (int id=3;id--;) dbdr[ib][id] = GetChebCalc(ib)->EvalDeriv(id,fArgsTmp)*fBScale[id];
238 0 : }
239 :
240 : //__________________________________________________________________________________________
241 : inline void AliCheb3D::EvalDeriv3D2(const Float_t *par, Float_t dbdrdr[3][3][3])
242 : {
243 : // return gradient matrix
244 0 : for (int i=3;i--;) fArgsTmp[i] = MapToInternal(par[i],i);
245 0 : for (int ib=3;ib--;) for (int id=3;id--;)for (int id1=3;id1--;)
246 0 : dbdrdr[ib][id][id1] = GetChebCalc(ib)->EvalDeriv2(id,id1,fArgsTmp)*fBScale[id]*fBScale[id1];
247 0 : }
248 :
249 : //__________________________________________________________________________________________
250 : inline void AliCheb3D::EvalDeriv(int dimd, const Float_t *par, Float_t *res)
251 : {
252 : // evaluate Chebyshev parameterization derivative for 3d->DimOut function
253 0 : for (int i=3;i--;) fArgsTmp[i] = MapToInternal(par[i],i);
254 0 : for (int i=fDimOut;i--;) res[i] = GetChebCalc(i)->EvalDeriv(dimd,fArgsTmp)*fBScale[dimd];;
255 : //
256 0 : }
257 :
258 : //__________________________________________________________________________________________
259 : inline void AliCheb3D::EvalDeriv2(int dimd1,int dimd2, const Float_t *par, Float_t *res)
260 : {
261 : // evaluate Chebyshev parameterization 2nd derivative over dimd1 and dimd2 dimensions for 3d->DimOut function
262 0 : for (int i=3;i--;) fArgsTmp[i] = MapToInternal(par[i],i);
263 0 : for (int i=fDimOut;i--;) res[i] = GetChebCalc(i)->EvalDeriv2(dimd1,dimd2,fArgsTmp)*fBScale[dimd1]*fBScale[dimd2];
264 : //
265 0 : }
266 :
267 : //__________________________________________________________________________________________
268 : inline Float_t AliCheb3D::EvalDeriv(int dimd, const Float_t *par, int idim)
269 : {
270 : // evaluate Chebyshev parameterization derivative over dimd dimention for idim-th output dimension of 3d->DimOut function
271 0 : for (int i=3;i--;) fArgsTmp[i] = MapToInternal(par[i],i);
272 0 : return GetChebCalc(idim)->EvalDeriv(dimd,fArgsTmp)*fBScale[dimd];
273 : //
274 : }
275 :
276 : //__________________________________________________________________________________________
277 : inline Float_t AliCheb3D::EvalDeriv2(int dimd1,int dimd2, const Float_t *par, int idim)
278 : {
279 : // evaluate Chebyshev parameterization 2ns derivative over dimd1 and dimd2 dimensions for idim-th output dimension of 3d->DimOut function
280 0 : for (int i=3;i--;) fArgsTmp[i] = MapToInternal(par[i],i);
281 0 : return GetChebCalc(idim)->EvalDeriv2(dimd1,dimd2,fArgsTmp)*fBScale[dimd1]*fBScale[dimd2];
282 : //
283 : }
284 :
285 : //__________________________________________________________________________________________
286 : inline Float_t AliCheb3D::MapToInternal(Float_t x,Int_t d) const
287 : {
288 : // map x to [-1:1]
289 : #ifdef _BRING_TO_BOUNDARY_
290 : float res = (x-fBOffset[d])*fBScale[d];
291 : if (res<-1) return -1;
292 : if (res> 1) return 1;
293 : return res;
294 : #else
295 0 : return (x-fBOffset[d])*fBScale[d];
296 : #endif
297 : }
298 :
299 : //__________________________________________________________________________________________
300 : inline Double_t AliCheb3D::MapToInternal(Double_t x,Int_t d) const
301 : {
302 : // map x to [-1:1]
303 : #ifdef _BRING_TO_BOUNDARY_
304 : double res = (x-fBOffset[d])*fBScale[d];
305 : if (res<-1) return -1;
306 : if (res> 1) return 1;
307 : return res;
308 : #else
309 18827412 : return (x-fBOffset[d])*fBScale[d];
310 : #endif
311 : }
312 :
313 : #endif
|
