]>
Commit | Line | Data |
---|---|---|
0eea9d4d | 1 | // Author: ruben.shahoyan@cern.ch 09/09/2006 |
40389866 | 2 | |
0eea9d4d | 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 | #include <TString.h> | |
64 | #include <TSystem.h> | |
0eea9d4d | 65 | #include <TROOT.h> |
40389866 | 66 | #include <TRandom.h> |
0eea9d4d | 67 | #include "AliCheb3D.h" |
68 | ||
69 | ||
70 | ||
0eea9d4d | 71 | ClassImp(AliCheb3D) |
72 | ||
40389866 | 73 | //__________________________________________________________________________________________ |
74 | AliCheb3D::AliCheb3D() : | |
75 | fDimOut(0), | |
76 | fPrec(0), | |
77 | fChebCalc(1), | |
78 | fMaxCoefs(0), | |
79 | fResTmp(0), | |
80 | fGrid(0), | |
81 | fUsrFunName(""), | |
82 | fUsrMacro(0) | |
0eea9d4d | 83 | { |
40389866 | 84 | for (int i=3;i--;) fBMin[i] = fBMax[i] = fBScale[i] = fBOffset[i] = 0; |
0eea9d4d | 85 | } |
86 | ||
40389866 | 87 | //__________________________________________________________________________________________ |
88 | AliCheb3D::AliCheb3D(const AliCheb3D& src) : | |
89 | TNamed(src), | |
90 | fDimOut(src.fDimOut), | |
91 | fPrec(src.fPrec), | |
92 | fChebCalc(1), | |
93 | fMaxCoefs(src.fMaxCoefs), | |
94 | fResTmp(0), | |
95 | fGrid(0), | |
96 | fUsrFunName(src.fUsrFunName), | |
97 | fUsrMacro(0) | |
0eea9d4d | 98 | { |
40389866 | 99 | // read coefs from text file |
100 | for (int i=3;i--;) { | |
101 | fBMin[i] = src.fBMin[i]; | |
102 | fBMax[i] = src.fBMax[i]; | |
103 | fBScale[i] = src.fBScale[i]; | |
104 | fBOffset[i] = src.fBOffset[i]; | |
105 | fNPoints[i] = src.fNPoints[i]; | |
106 | } | |
107 | for (int i=0;i<fDimOut;i++) { | |
108 | AliCheb3DCalc* cbc = src.GetChebCalc(i); | |
109 | if (cbc) fChebCalc.AddAtAndExpand(new AliCheb3DCalc(*cbc),i); | |
110 | } | |
0eea9d4d | 111 | } |
112 | ||
40389866 | 113 | //__________________________________________________________________________________________ |
114 | AliCheb3D::AliCheb3D(const char* inpFile) : | |
115 | fDimOut(0), | |
116 | fPrec(0), | |
117 | fChebCalc(1), | |
118 | fMaxCoefs(0), | |
119 | fResTmp(0), | |
120 | fGrid(0), | |
121 | fUsrFunName(""), | |
122 | fUsrMacro(0) | |
0bc7b414 | 123 | { |
40389866 | 124 | // read coefs from text file |
125 | for (int i=3;i--;) fBMin[i] = fBMax[i] = fBScale[i] = fBOffset[i] = 0; | |
126 | LoadData(inpFile); | |
c437b1a5 | 127 | } |
128 | ||
40389866 | 129 | //__________________________________________________________________________________________ |
130 | AliCheb3D::AliCheb3D(FILE* stream) : | |
131 | fDimOut(0), | |
132 | fPrec(0), | |
133 | fChebCalc(1), | |
134 | fMaxCoefs(0), | |
135 | fResTmp(0), | |
136 | fGrid(0), | |
137 | fUsrFunName(""), | |
138 | fUsrMacro(0) | |
c437b1a5 | 139 | { |
40389866 | 140 | // read coefs from stream |
141 | for (int i=3;i--;) fBMin[i] = fBMax[i] = fBScale[i] = fBOffset[i] = 0; | |
142 | LoadData(stream); | |
0bc7b414 | 143 | } |
0eea9d4d | 144 | |
145 | //__________________________________________________________________________________________ | |
146 | #ifdef _INC_CREATION_ALICHEB3D_ | |
40389866 | 147 | AliCheb3D::AliCheb3D(const char* funName, int DimOut, Float_t *bmin,Float_t *bmax, Int_t *npoints, Float_t prec) : |
148 | TNamed(funName,funName), | |
149 | fDimOut(0), | |
150 | fPrec(TMath::Max(1.E-12f,prec)), | |
151 | fChebCalc(1), | |
152 | fMaxCoefs(0), | |
153 | fResTmp(0), | |
154 | fGrid(0), | |
155 | fUsrFunName("") , | |
156 | fUsrMacro(0) | |
0eea9d4d | 157 | { |
158 | // Construct the parameterization for the function | |
159 | // funName : name of the file containing the function: void funName(Float_t * inp,Float_t * out) | |
160 | // DimOut : dimension of the vector computed by the user function | |
161 | // bmin : array of 3 elements with the lower boundaries of the region where the function is defined | |
162 | // bmax : array of 3 elements with the upper boundaries of the region where the function is defined | |
163 | // npoints : array of 3 elements with the number of points to compute in each of 3 dimension | |
164 | // prec : max allowed absolute difference between the user function and computed parameterization on the requested grid | |
165 | // | |
0eea9d4d | 166 | if (DimOut<1) {Error("AliCheb3D","Requested output dimension is %d\nStop\n",fDimOut); exit(1);} |
167 | SetDimOut(DimOut); | |
168 | PrepareBoundaries(bmin,bmax); | |
169 | DefineGrid(npoints); | |
170 | SetUsrFunction(funName); | |
171 | ChebFit(); | |
172 | // | |
173 | } | |
174 | #endif | |
175 | ||
176 | //__________________________________________________________________________________________ | |
177 | #ifdef _INC_CREATION_ALICHEB3D_ | |
40389866 | 178 | AliCheb3D::AliCheb3D(void (*ptr)(float*,float*), int DimOut, Float_t *bmin,Float_t *bmax, Int_t *npoints, Float_t prec) : |
179 | fDimOut(0), | |
180 | fPrec(TMath::Max(1.E-12f,prec)), | |
181 | fChebCalc(1), | |
182 | fMaxCoefs(0), | |
183 | fResTmp(0), | |
184 | fGrid(0), | |
185 | fUsrFunName(""), | |
186 | fUsrMacro(0) | |
0eea9d4d | 187 | { |
188 | // Construct the parameterization for the function | |
189 | // ptr : pointer on the function: void fun(Float_t * inp,Float_t * out) | |
190 | // DimOut : dimension of the vector computed by the user function | |
191 | // bmin : array of 3 elements with the lower boundaries of the region where the function is defined | |
192 | // bmax : array of 3 elements with the upper boundaries of the region where the function is defined | |
193 | // npoints : array of 3 elements with the number of points to compute in each of 3 dimension | |
194 | // prec : max allowed absolute difference between the user function and computed parameterization on the requested grid | |
195 | // | |
0eea9d4d | 196 | if (DimOut<1) {Error("AliCheb3D","Requested output dimension is %d\nStop\n",fDimOut); exit(1);} |
197 | SetDimOut(DimOut); | |
198 | PrepareBoundaries(bmin,bmax); | |
199 | DefineGrid(npoints); | |
200 | SetUsrFunction(ptr); | |
201 | ChebFit(); | |
202 | // | |
203 | } | |
204 | #endif | |
205 | ||
40389866 | 206 | //__________________________________________________________________________________________ |
207 | #ifdef _INC_CREATION_ALICHEB3D_ | |
208 | AliCheb3D::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) : | |
209 | fDimOut(0), | |
210 | fPrec(TMath::Max(1.E-12f,prec)), | |
211 | fChebCalc(1), | |
212 | fMaxCoefs(0), | |
213 | fResTmp(0), | |
214 | fGrid(0), | |
215 | fUsrFunName(""), | |
216 | fUsrMacro(0) | |
217 | { | |
218 | // Construct very economic parameterization for the function | |
219 | // ptr : pointer on the function: void fun(Float_t * inp,Float_t * out) | |
220 | // DimOut : dimension of the vector computed by the user function | |
221 | // bmin : array of 3 elements with the lower boundaries of the region where the function is defined | |
222 | // bmax : array of 3 elements with the upper boundaries of the region where the function is defined | |
223 | // npX : array of 3 elements with the number of points to compute in each dimension for 1st component | |
224 | // npY : array of 3 elements with the number of points to compute in each dimension for 2nd component | |
225 | // npZ : array of 3 elements with the number of points to compute in each dimension for 3d component | |
226 | // prec : max allowed absolute difference between the user function and computed parameterization on the requested grid | |
227 | // | |
228 | if (DimOut<1) {Error("AliCheb3D","Requested output dimension is %d\nStop\n",fDimOut); exit(1);} | |
229 | SetDimOut(DimOut); | |
230 | PrepareBoundaries(bmin,bmax); | |
231 | SetUsrFunction(ptr); | |
232 | // | |
233 | DefineGrid(npX); | |
234 | ChebFit(0); | |
235 | DefineGrid(npY); | |
236 | ChebFit(1); | |
237 | DefineGrid(npZ); | |
238 | ChebFit(2); | |
239 | // | |
240 | } | |
241 | #endif | |
242 | ||
243 | ||
244 | //__________________________________________________________________________________________ | |
245 | #ifdef _INC_CREATION_ALICHEB3D_ | |
246 | AliCheb3D::AliCheb3D(void (*ptr)(float*,float*), int DimOut, Float_t *bmin,Float_t *bmax, Float_t prec) : | |
247 | fDimOut(0), | |
248 | fPrec(TMath::Max(1.E-12f,prec)), | |
249 | fChebCalc(1), | |
250 | fMaxCoefs(0), | |
251 | fResTmp(0), | |
252 | fGrid(0), | |
253 | fUsrFunName(""), | |
254 | fUsrMacro(0) | |
255 | { | |
256 | // Construct very economic parameterization for the function with automatic calculation of the root's grid | |
257 | // ptr : pointer on the function: void fun(Float_t * inp,Float_t * out) | |
258 | // DimOut : dimension of the vector computed by the user function | |
259 | // bmin : array of 3 elements with the lower boundaries of the region where the function is defined | |
260 | // bmax : array of 3 elements with the upper boundaries of the region where the function is defined | |
261 | // prec : max allowed absolute difference between the user function and computed parameterization on the requested grid | |
262 | // | |
263 | if (DimOut!=3) {Error("AliCheb3D","This constructor works only for 3D fits, %dD fit was requested\n",fDimOut); exit(1);} | |
264 | SetDimOut(DimOut); | |
265 | PrepareBoundaries(bmin,bmax); | |
266 | SetUsrFunction(ptr); | |
267 | // | |
268 | int gridNC[3][3]; | |
269 | EstimateNPoints(prec,gridNC); | |
270 | DefineGrid(gridNC[0]); | |
271 | ChebFit(0); | |
272 | DefineGrid(gridNC[1]); | |
273 | ChebFit(1); | |
274 | DefineGrid(gridNC[2]); | |
275 | ChebFit(2); | |
276 | // | |
277 | } | |
278 | #endif | |
279 | ||
280 | ||
281 | //__________________________________________________________________________________________ | |
282 | AliCheb3D& AliCheb3D::operator=(const AliCheb3D& rhs) | |
283 | { | |
284 | if (this != &rhs) { | |
285 | Clear(); | |
286 | fDimOut = rhs.fDimOut; | |
287 | fPrec = rhs.fPrec; | |
288 | fMaxCoefs = rhs.fMaxCoefs; | |
289 | fUsrFunName = rhs.fUsrFunName; | |
290 | fUsrMacro = 0; | |
291 | for (int i=3;i--;) { | |
292 | fBMin[i] = rhs.fBMin[i]; | |
293 | fBMax[i] = rhs.fBMax[i]; | |
294 | fBScale[i] = rhs.fBScale[i]; | |
295 | fBOffset[i] = rhs.fBOffset[i]; | |
296 | fNPoints[i] = rhs.fNPoints[i]; | |
297 | } | |
298 | for (int i=0;i<fDimOut;i++) { | |
299 | AliCheb3DCalc* cbc = rhs.GetChebCalc(i); | |
300 | if (cbc) fChebCalc.AddAtAndExpand(new AliCheb3DCalc(*cbc),i); | |
301 | } | |
302 | } | |
303 | return *this; | |
304 | // | |
305 | } | |
0bc7b414 | 306 | |
0eea9d4d | 307 | //__________________________________________________________________________________________ |
308 | void AliCheb3D::Clear(Option_t*) | |
309 | { | |
310 | if (fResTmp) { delete[] fResTmp; fResTmp = 0; } | |
311 | if (fGrid) { delete[] fGrid; fGrid = 0; } | |
312 | if (fUsrMacro) { delete fUsrMacro; fUsrMacro = 0;} | |
313 | fChebCalc.Delete(); | |
314 | // | |
315 | } | |
316 | ||
317 | //__________________________________________________________________________________________ | |
318 | void AliCheb3D::Print(Option_t* opt) const | |
319 | { | |
320 | printf("%s: Chebyshev parameterization for 3D->%dD function. Precision: %e\n",GetName(),fDimOut,fPrec); | |
321 | printf("Region of validity: [%+.5e:%+.5e] [%+.5e:%+.5e] [%+.5e:%+.5e]\n",fBMin[0],fBMax[0],fBMin[1],fBMax[1],fBMin[2],fBMax[2]); | |
322 | TString opts = opt; opts.ToLower(); | |
323 | if (opts.Contains("l")) for (int i=0;i<fDimOut;i++) {printf("Output dimension %d:\n",i+1); GetChebCalc(i)->Print();} | |
324 | // | |
325 | } | |
326 | ||
0eea9d4d | 327 | //__________________________________________________________________________________________ |
328 | void AliCheb3D::PrepareBoundaries(Float_t *bmin,Float_t *bmax) | |
329 | { | |
330 | // Set and check boundaries defined by user, prepare coefficients for their conversion to [-1:1] interval | |
331 | // | |
332 | for (int i=3;i--;) { | |
333 | fBMin[i] = bmin[i]; | |
334 | fBMax[i] = bmax[i]; | |
335 | fBScale[i] = bmax[i]-bmin[i]; | |
336 | if (fBScale[i]<=0) { | |
337 | Error("PrepareBoundaries","Boundaries for %d-th dimension are not increasing: %+.4e %+.4e\nStop\n",i,fBMin[i],fBMax[i]); | |
338 | exit(1); | |
339 | } | |
340 | fBOffset[i] = bmin[i] + fBScale[i]/2.0; | |
341 | fBScale[i] = 2./fBScale[i]; | |
342 | } | |
343 | // | |
344 | } | |
345 | ||
40389866 | 346 | |
0eea9d4d | 347 | //__________________________________________________________________________________________ |
348 | #ifdef _INC_CREATION_ALICHEB3D_ | |
40389866 | 349 | |
350 | // Pointer on user function (faster altrnative to TMethodCall) | |
351 | void (*gUsrFunAliCheb3D) (float* ,float* ); | |
352 | ||
353 | void AliCheb3D::EvalUsrFunction() | |
354 | { | |
355 | // call user supplied function | |
356 | if (gUsrFunAliCheb3D) gUsrFunAliCheb3D(fArgsTmp,fResTmp); | |
357 | else fUsrMacro->Execute(); | |
358 | } | |
359 | ||
0eea9d4d | 360 | void AliCheb3D::SetUsrFunction(const char* name) |
361 | { | |
362 | // load user macro with function definition and compile it | |
363 | gUsrFunAliCheb3D = 0; | |
364 | fUsrFunName = name; | |
365 | gSystem->ExpandPathName(fUsrFunName); | |
366 | if (fUsrMacro) delete fUsrMacro; | |
367 | TString tmpst = fUsrFunName; | |
368 | tmpst += "+"; // prepare filename to compile | |
369 | if (gROOT->LoadMacro(tmpst.Data())) {Error("SetUsrFunction","Failed to load user function from %s\nStop\n",name); exit(1);} | |
370 | fUsrMacro = new TMethodCall(); | |
371 | tmpst = tmpst.Data() + tmpst.Last('/')+1; //Strip away any path preceding the macro file name | |
372 | int dot = tmpst.Last('.'); | |
373 | if (dot>0) tmpst.Resize(dot); | |
374 | fUsrMacro->InitWithPrototype(tmpst.Data(),"Float_t *,Float_t *"); | |
375 | long args[2]; | |
376 | args[0] = (long)fArgsTmp; | |
377 | args[1] = (long)fResTmp; | |
378 | fUsrMacro->SetParamPtrs(args); | |
379 | // | |
380 | } | |
381 | #endif | |
382 | ||
383 | //__________________________________________________________________________________________ | |
384 | #ifdef _INC_CREATION_ALICHEB3D_ | |
385 | void AliCheb3D::SetUsrFunction(void (*ptr)(float*,float*)) | |
386 | { | |
387 | if (fUsrMacro) delete fUsrMacro; | |
388 | fUsrMacro = 0; | |
389 | fUsrFunName = ""; | |
390 | gUsrFunAliCheb3D = ptr; | |
391 | } | |
392 | #endif | |
393 | ||
394 | //__________________________________________________________________________________________ | |
395 | #ifdef _INC_CREATION_ALICHEB3D_ | |
396 | void AliCheb3D::EvalUsrFunction(Float_t *x, Float_t *res) { | |
397 | for (int i=3;i--;) fArgsTmp[i] = x[i]; | |
398 | if (gUsrFunAliCheb3D) gUsrFunAliCheb3D(fArgsTmp,fResTmp); | |
399 | else fUsrMacro->Execute(); | |
400 | for (int i=fDimOut;i--;) res[i] = fResTmp[i]; | |
401 | } | |
402 | #endif | |
403 | ||
404 | //__________________________________________________________________________________________ | |
405 | #ifdef _INC_CREATION_ALICHEB3D_ | |
406 | Int_t AliCheb3D::CalcChebCoefs(Float_t *funval,int np, Float_t *outCoefs, Float_t prec) | |
407 | { | |
408 | // Calculate Chebyshev coeffs using precomputed function values at np roots. | |
409 | // If prec>0, estimate the highest coeff number providing the needed precision | |
410 | // | |
411 | double sm; // do summations in double to minimize the roundoff error | |
412 | for (int ic=0;ic<np;ic++) { // compute coeffs | |
413 | sm = 0; | |
414 | for (int ir=0;ir<np;ir++) { | |
415 | float rt = TMath::Cos( ic*(ir+0.5)*TMath::Pi()/np); | |
416 | sm += funval[ir]*rt; | |
417 | } | |
418 | outCoefs[ic] = Float_t( sm * ((ic==0) ? 1./np : 2./np) ); | |
419 | } | |
420 | // | |
421 | if (prec<=0) return np; | |
422 | // | |
423 | sm = 0; | |
424 | int cfMax = 0; | |
425 | for (cfMax=np;cfMax--;) { | |
426 | sm += TMath::Abs(outCoefs[cfMax]); | |
427 | if (sm>=prec) break; | |
428 | } | |
429 | if (++cfMax==0) cfMax=1; | |
430 | return cfMax; | |
431 | // | |
432 | } | |
433 | #endif | |
434 | ||
435 | //__________________________________________________________________________________________ | |
436 | #ifdef _INC_CREATION_ALICHEB3D_ | |
437 | void AliCheb3D::DefineGrid(Int_t* npoints) | |
438 | { | |
439 | // prepare the grid of Chebyshev roots in each dimension | |
440 | const int kMinPoints = 1; | |
441 | int ntot = 0; | |
442 | fMaxCoefs = 1; | |
443 | for (int id=3;id--;) { | |
444 | fNPoints[id] = npoints[id]; | |
445 | if (fNPoints[id]<kMinPoints) { | |
40389866 | 446 | Error("DefineGrid","at %d-th dimension %d point is requested, at least %d is needed\nStop\n",id,fNPoints[id],kMinPoints); |
0eea9d4d | 447 | exit(1); |
448 | } | |
449 | ntot += fNPoints[id]; | |
450 | fMaxCoefs *= fNPoints[id]; | |
451 | } | |
40389866 | 452 | printf("Computing Chebyshev nodes on [%2d/%2d/%2d] grid\n",npoints[0],npoints[1],npoints[2]); |
453 | if (fGrid) delete[] fGrid; | |
0eea9d4d | 454 | fGrid = new Float_t [ntot]; |
455 | // | |
456 | int curp = 0; | |
457 | for (int id=3;id--;) { | |
458 | int np = fNPoints[id]; | |
459 | fGridOffs[id] = curp; | |
460 | for (int ip=0;ip<np;ip++) { | |
461 | Float_t x = TMath::Cos( TMath::Pi()*(ip+0.5)/np ); | |
462 | fGrid[curp++] = MapToExternal(x,id); | |
463 | } | |
464 | } | |
465 | // | |
466 | } | |
467 | #endif | |
468 | ||
469 | //__________________________________________________________________________________________ | |
470 | #ifdef _INC_CREATION_ALICHEB3D_ | |
471 | Int_t AliCheb3D::ChebFit() | |
472 | { | |
473 | // prepare parameterization for all output dimensions | |
474 | int ir=0; | |
475 | for (int i=fDimOut;i--;) ir+=ChebFit(i); | |
476 | return ir; | |
477 | } | |
478 | #endif | |
479 | ||
480 | //__________________________________________________________________________________________ | |
481 | #ifdef _INC_CREATION_ALICHEB3D_ | |
482 | Int_t AliCheb3D::ChebFit(int dmOut) | |
483 | { | |
484 | // prepare paramaterization of 3D function for dmOut-th dimension | |
485 | int maxDim = 0; | |
486 | for (int i=0;i<3;i++) if (maxDim<fNPoints[i]) maxDim = fNPoints[i]; | |
487 | Float_t *fvals = new Float_t [ fNPoints[0] ]; | |
488 | Float_t *tmpCoef3D = new Float_t [ fNPoints[0]*fNPoints[1]*fNPoints[2] ]; | |
489 | Float_t *tmpCoef2D = new Float_t [ fNPoints[0]*fNPoints[1] ]; | |
490 | Float_t *tmpCoef1D = new Float_t [ maxDim ]; | |
491 | // | |
492 | Float_t RTiny = fPrec/Float_t(maxDim); // neglect coefficient below this threshold | |
493 | // | |
494 | // 1D Cheb.fit for 0-th dimension at current steps of remaining dimensions | |
495 | int ncmax = 0; | |
496 | // | |
40389866 | 497 | printf("Dim%d : 00.00%% Done",dmOut);fflush(stdout); |
0eea9d4d | 498 | AliCheb3DCalc* cheb = GetChebCalc(dmOut); |
499 | // | |
40389866 | 500 | float ncals2count = fNPoints[2]*fNPoints[1]*fNPoints[0]; |
501 | float ncals = 0; | |
502 | float frac = 0; | |
503 | float fracStep = 0.001; | |
504 | // | |
0eea9d4d | 505 | for (int id2=fNPoints[2];id2--;) { |
506 | fArgsTmp[2] = fGrid[ fGridOffs[2]+id2 ]; | |
507 | // | |
508 | for (int id1=fNPoints[1];id1--;) { | |
509 | fArgsTmp[1] = fGrid[ fGridOffs[1]+id1 ]; | |
510 | // | |
511 | for (int id0=fNPoints[0];id0--;) { | |
512 | fArgsTmp[0] = fGrid[ fGridOffs[0]+id0 ]; | |
513 | EvalUsrFunction(); // compute function values at Chebyshev roots of 0-th dimension | |
514 | fvals[id0] = fResTmp[dmOut]; | |
40389866 | 515 | float fr = (++ncals)/ncals2count; |
516 | if (fr-frac>=fracStep) { | |
517 | frac = fr; | |
518 | printf("\b\b\b\b\b\b\b\b\b\b\b"); | |
519 | printf("%05.2f%% Done",fr*100); | |
520 | fflush(stdout); | |
521 | } | |
522 | // | |
0eea9d4d | 523 | } |
524 | int nc = CalcChebCoefs(fvals,fNPoints[0], tmpCoef1D, fPrec); | |
525 | for (int id0=fNPoints[0];id0--;) tmpCoef2D[id1 + id0*fNPoints[1]] = tmpCoef1D[id0]; | |
526 | if (ncmax<nc) ncmax = nc; // max coefs to be kept in dim0 to guarantee needed precision | |
527 | } | |
528 | // | |
529 | // once each 1d slice of given 2d slice is parametrized, parametrize the Cheb.coeffs | |
530 | for (int id0=fNPoints[0];id0--;) { | |
531 | CalcChebCoefs( tmpCoef2D+id0*fNPoints[1], fNPoints[1], tmpCoef1D, -1); | |
532 | for (int id1=fNPoints[1];id1--;) tmpCoef3D[id2 + fNPoints[2]*(id1+id0*fNPoints[1])] = tmpCoef1D[id1]; | |
533 | } | |
534 | } | |
535 | // | |
536 | // now fit the last dimensions Cheb.coefs | |
537 | for (int id0=fNPoints[0];id0--;) { | |
538 | for (int id1=fNPoints[1];id1--;) { | |
539 | CalcChebCoefs( tmpCoef3D+ fNPoints[2]*(id1+id0*fNPoints[1]), fNPoints[2], tmpCoef1D, -1); | |
540 | for (int id2=fNPoints[2];id2--;) tmpCoef3D[id2+ fNPoints[2]*(id1+id0*fNPoints[1])] = tmpCoef1D[id2]; // store on place | |
541 | } | |
542 | } | |
543 | // | |
544 | // now find 2D surface which separates significant coefficients of 3D matrix from nonsignificant ones (up to fPrec) | |
545 | int *tmpCoefSurf = new Int_t[ fNPoints[0]*fNPoints[1] ]; | |
546 | for (int id0=fNPoints[0];id0--;) for (int id1=fNPoints[1];id1--;) tmpCoefSurf[id1+id0*fNPoints[1]]=0; | |
547 | Double_t resid = 0; | |
548 | for (int id0=fNPoints[0];id0--;) { | |
549 | for (int id1=fNPoints[1];id1--;) { | |
550 | for (int id2=fNPoints[2];id2--;) { | |
551 | int id = id2 + fNPoints[2]*(id1+id0*fNPoints[1]); | |
552 | Float_t cfa = TMath::Abs(tmpCoef3D[id]); | |
40389866 | 553 | if (cfa < RTiny) {tmpCoef3D[id] = 0; continue;} // neglect coefs below the threshold |
0eea9d4d | 554 | resid += cfa; |
555 | if (resid<fPrec) continue; // this coeff is negligible | |
556 | // otherwise go back 1 step | |
557 | resid -= cfa; | |
558 | tmpCoefSurf[id1+id0*fNPoints[1]] = id2+1; // how many coefs to keep | |
559 | break; | |
560 | } | |
561 | } | |
562 | } | |
563 | /* | |
564 | printf("\n\nCoeffs\n"); | |
565 | int cnt = 0; | |
566 | for (int id0=0;id0<fNPoints[0];id0++) { | |
567 | for (int id1=0;id1<fNPoints[1];id1++) { | |
568 | for (int id2=0;id2<fNPoints[2];id2++) { | |
569 | printf("%2d%2d%2d %+.4e |",id0,id1,id2,tmpCoef3D[cnt++]); | |
570 | } | |
571 | printf("\n"); | |
572 | } | |
573 | printf("\n"); | |
574 | } | |
575 | */ | |
576 | // see if there are rows to reject, find max.significant column at each row | |
577 | int NRows = fNPoints[0]; | |
578 | int *tmpCols = new int[NRows]; | |
579 | for (int id0=fNPoints[0];id0--;) { | |
580 | int id1 = fNPoints[1]; | |
581 | while (id1>0 && tmpCoefSurf[(id1-1)+id0*fNPoints[1]]==0) id1--; | |
582 | tmpCols[id0] = id1; | |
583 | } | |
584 | // find max significant row | |
585 | for (int id0=NRows;id0--;) {if (tmpCols[id0]>0) break; NRows--;} | |
586 | // find max significant column and fill the permanent storage for the max sigificant column of each row | |
587 | cheb->InitRows(NRows); // create needed arrays; | |
588 | int *NColsAtRow = cheb->GetNColsAtRow(); | |
589 | int *ColAtRowBg = cheb->GetColAtRowBg(); | |
590 | int NCols = 0; | |
591 | int NElemBound2D = 0; | |
592 | for (int id0=0;id0<NRows;id0++) { | |
593 | NColsAtRow[id0] = tmpCols[id0]; // number of columns to store for this row | |
594 | ColAtRowBg[id0] = NElemBound2D; // begining of this row in 2D boundary surface | |
595 | NElemBound2D += tmpCols[id0]; | |
596 | if (NCols<NColsAtRow[id0]) NCols = NColsAtRow[id0]; | |
597 | } | |
598 | cheb->InitCols(NCols); | |
599 | delete[] tmpCols; | |
600 | // | |
601 | // create the 2D matrix defining the boundary of significance for 3D coeffs.matrix | |
602 | // and count the number of siginifacnt coefficients | |
603 | // | |
604 | cheb->InitElemBound2D(NElemBound2D); | |
605 | int *CoefBound2D0 = cheb->GetCoefBound2D0(); | |
606 | int *CoefBound2D1 = cheb->GetCoefBound2D1(); | |
607 | fMaxCoefs = 0; // redefine number of coeffs | |
608 | for (int id0=0;id0<NRows;id0++) { | |
609 | int nCLoc = NColsAtRow[id0]; | |
40389866 | 610 | int Col0 = ColAtRowBg[id0]; |
0eea9d4d | 611 | for (int id1=0;id1<nCLoc;id1++) { |
40389866 | 612 | CoefBound2D0[Col0 + id1] = tmpCoefSurf[id1+id0*fNPoints[1]]; // number of coefs to store for 3-d dimension |
613 | CoefBound2D1[Col0 + id1] = fMaxCoefs; | |
614 | fMaxCoefs += CoefBound2D0[Col0 + id1]; | |
0eea9d4d | 615 | } |
616 | } | |
617 | // | |
618 | // create final compressed 3D matrix for significant coeffs | |
619 | cheb->InitCoefs(fMaxCoefs); | |
620 | Float_t *Coefs = cheb->GetCoefs(); | |
621 | int count = 0; | |
622 | for (int id0=0;id0<NRows;id0++) { | |
623 | int ncLoc = NColsAtRow[id0]; | |
40389866 | 624 | int Col0 = ColAtRowBg[id0]; |
0eea9d4d | 625 | for (int id1=0;id1<ncLoc;id1++) { |
40389866 | 626 | int ncf2 = CoefBound2D0[Col0 + id1]; |
0eea9d4d | 627 | for (int id2=0;id2<ncf2;id2++) { |
628 | Coefs[count++] = tmpCoef3D[id2 + fNPoints[2]*(id1+id0*fNPoints[1])]; | |
629 | } | |
630 | } | |
631 | } | |
632 | /* | |
633 | printf("\n\nNewSurf\n"); | |
634 | for (int id0=0;id0<fNPoints[0];id0++) { | |
635 | for (int id1=0;id1<fNPoints[1];id1++) { | |
636 | printf("(%2d %2d) %2d |",id0,id1,tmpCoefSurf[id1+id0*fNPoints[1]]); | |
637 | } | |
638 | printf("\n"); | |
639 | } | |
640 | */ | |
641 | // | |
642 | delete[] tmpCoefSurf; | |
643 | delete[] tmpCoef1D; | |
644 | delete[] tmpCoef2D; | |
645 | delete[] tmpCoef3D; | |
646 | delete[] fvals; | |
647 | // | |
40389866 | 648 | printf("\b\b\b\b\b\b\b\b\b\b\b\b"); |
649 | printf("100.00%% Done\n"); | |
0eea9d4d | 650 | return 1; |
651 | } | |
652 | #endif | |
653 | ||
654 | //_______________________________________________ | |
655 | #ifdef _INC_CREATION_ALICHEB3D_ | |
656 | void AliCheb3D::SaveData(const char* outfile,Bool_t append) const | |
657 | { | |
658 | // writes coefficients data to output text file, optionallt appending on the end of existing file | |
659 | TString strf = outfile; | |
660 | gSystem->ExpandPathName(strf); | |
661 | FILE* stream = fopen(strf,append ? "a":"w"); | |
662 | SaveData(stream); | |
663 | fclose(stream); | |
664 | // | |
665 | } | |
666 | #endif | |
667 | ||
668 | //_______________________________________________ | |
669 | #ifdef _INC_CREATION_ALICHEB3D_ | |
670 | void AliCheb3D::SaveData(FILE* stream) const | |
671 | { | |
672 | // writes coefficients data to existing output stream | |
673 | // | |
674 | fprintf(stream,"\n# These are automatically generated data for the Chebyshev interpolation of 3D->%dD function\n",fDimOut); | |
675 | fprintf(stream,"#\nSTART %s\n",GetName()); | |
676 | fprintf(stream,"# Dimensionality of the output\n%d\n",fDimOut); | |
677 | fprintf(stream,"# Interpolation abs. precision\n%+.8e\n",fPrec); | |
678 | // | |
679 | fprintf(stream,"# Lower boundaries of interpolation region\n"); | |
680 | for (int i=0;i<3;i++) fprintf(stream,"%+.8e\n",fBMin[i]); | |
681 | fprintf(stream,"# Upper boundaries of interpolation region\n"); | |
682 | for (int i=0;i<3;i++) fprintf(stream,"%+.8e\n",fBMax[i]); | |
40389866 | 683 | fprintf(stream,"# Parameterization for each output dimension follows:\n"); |
0eea9d4d | 684 | // |
685 | for (int i=0;i<fDimOut;i++) GetChebCalc(i)->SaveData(stream); | |
686 | fprintf(stream,"#\nEND %s\n#\n",GetName()); | |
687 | // | |
688 | } | |
689 | #endif | |
690 | ||
691 | //_______________________________________________ | |
692 | void AliCheb3D::LoadData(const char* inpFile) | |
693 | { | |
694 | TString strf = inpFile; | |
695 | gSystem->ExpandPathName(strf); | |
696 | FILE* stream = fopen(strf.Data(),"r"); | |
697 | LoadData(stream); | |
698 | fclose(stream); | |
699 | // | |
700 | } | |
701 | ||
702 | //_______________________________________________ | |
703 | void AliCheb3D::LoadData(FILE* stream) | |
704 | { | |
705 | if (!stream) {Error("LoadData","No stream provided.\nStop"); exit(1);} | |
706 | TString buffs; | |
707 | Clear(); | |
708 | AliCheb3DCalc::ReadLine(buffs,stream); | |
709 | if (!buffs.BeginsWith("START")) {Error("LoadData","Expected: \"START <fit_name>\", found \"%s\"\nStop\n",buffs.Data());exit(1);} | |
710 | SetName(buffs.Data()+buffs.First(' ')+1); | |
711 | // | |
712 | AliCheb3DCalc::ReadLine(buffs,stream); // N output dimensions | |
713 | fDimOut = buffs.Atoi(); | |
714 | if (fDimOut<1) {Error("LoadData","Expected: '<number_of_output_dimensions>', found \"%s\"\nStop\n",buffs.Data());exit(1);} | |
715 | // | |
716 | SetDimOut(fDimOut); | |
717 | // | |
718 | AliCheb3DCalc::ReadLine(buffs,stream); // Interpolation abs. precision | |
719 | fPrec = buffs.Atof(); | |
720 | if (fPrec<=0) {Error("LoadData","Expected: '<abs.precision>', found \"%s\"\nStop\n",buffs.Data());exit(1);} | |
721 | // | |
722 | for (int i=0;i<3;i++) { // Lower boundaries of interpolation region | |
723 | AliCheb3DCalc::ReadLine(buffs,stream); | |
724 | fBMin[i] = buffs.Atof(); | |
725 | } | |
726 | for (int i=0;i<3;i++) { // Upper boundaries of interpolation region | |
727 | AliCheb3DCalc::ReadLine(buffs,stream); | |
728 | fBMax[i] = buffs.Atof(); | |
729 | } | |
730 | PrepareBoundaries(fBMin,fBMax); | |
731 | // | |
732 | // data for each output dimension | |
733 | for (int i=0;i<fDimOut;i++) GetChebCalc(i)->LoadData(stream); | |
734 | // | |
735 | // check end_of_data record | |
736 | AliCheb3DCalc::ReadLine(buffs,stream); | |
737 | if (!buffs.BeginsWith("END") || !buffs.Contains(GetName())) { | |
738 | Error("LoadData","Expected \"END %s\", found \"%s\".\nStop\n",GetName(),buffs.Data()); | |
739 | exit(1); | |
740 | } | |
741 | // | |
742 | } | |
743 | ||
744 | //_______________________________________________ | |
745 | void AliCheb3D::SetDimOut(int d) | |
746 | { | |
747 | fDimOut = d; | |
748 | if (fResTmp) delete fResTmp; | |
40389866 | 749 | fResTmp = new Float_t[fDimOut]; |
0eea9d4d | 750 | fChebCalc.Delete(); |
751 | for (int i=0;i<d;i++) fChebCalc.AddAtAndExpand(new AliCheb3DCalc(),i); | |
752 | } | |
753 | ||
754 | //_______________________________________________ | |
755 | void AliCheb3D::ShiftBound(int id,float dif) | |
756 | { | |
757 | if (id<0||id>2) {printf("Maximum 3 dimensions are supported\n"); return;} | |
758 | fBMin[id] += dif; | |
759 | fBMax[id] += dif; | |
760 | fBOffset[id] += dif; | |
761 | } | |
762 | ||
763 | //_______________________________________________ | |
764 | #ifdef _INC_CREATION_ALICHEB3D_ | |
765 | TH1* AliCheb3D::TestRMS(int idim,int npoints,TH1* histo) | |
766 | { | |
767 | // fills the difference between the original function and parameterization (for idim-th component of the output) | |
768 | // to supplied histogram. Calculations are done in npoints random points. | |
769 | // If the hostgram was not supplied, it will be created. It is up to the user to delete it! | |
770 | if (!fUsrMacro) { | |
771 | printf("No user function is set\n"); | |
772 | return 0; | |
773 | } | |
774 | if (!histo) histo = new TH1D(GetName(),"Control: Function - Parametrization",100,-2*fPrec,2*fPrec); | |
775 | for (int ip=npoints;ip--;) { | |
776 | gRandom->RndmArray(3,(Float_t *)fArgsTmp); | |
777 | for (int i=3;i--;) fArgsTmp[i] = fBMin[i] + fArgsTmp[i]*(fBMax[i]-fBMin[i]); | |
778 | EvalUsrFunction(); | |
779 | Float_t valFun = fResTmp[idim]; | |
780 | Eval(fArgsTmp,fResTmp); | |
781 | Float_t valPar = fResTmp[idim]; | |
782 | histo->Fill(valFun - valPar); | |
783 | } | |
784 | return histo; | |
785 | // | |
786 | } | |
787 | #endif | |
40389866 | 788 | |
789 | //_______________________________________________ | |
790 | #ifdef _INC_CREATION_ALICHEB3D_ | |
791 | void AliCheb3D::EstimateNPoints(float Prec, int gridBC[3][3]) | |
792 | { | |
793 | const float sclA[9] = {0.1, 0.5, 0.9, 0.1, 0.5, 0.9, 0.1, 0.5, 0.9} ; | |
794 | const float sclB[9] = {0.1, 0.1, 0.1, 0.5, 0.5, 0.5, 0.9, 0.9, 0.9} ; | |
795 | const float sclDim[2] = {0.01,0.99}; | |
796 | const int compDim[3][2] = { {1,2}, {2,0}, {0,1} }; | |
797 | static float xyz[3]; | |
798 | // | |
799 | for (int i=3;i--;)for (int j=3;j--;) gridBC[i][j] = -1; | |
800 | // | |
801 | for (int idim=0;idim<3;idim++) { | |
802 | float dimMN = fBMin[idim] + sclDim[0]*(fBMax[idim]-fBMin[idim]); | |
803 | float dimMX = fBMin[idim] + sclDim[1]*(fBMax[idim]-fBMin[idim]); | |
804 | // | |
805 | for (int it=0;it<9;it++) { // test in 9 points | |
806 | int id1 = compDim[idim][0]; // 1st fixed dim | |
807 | int id2 = compDim[idim][1]; // 2nd fixed dim | |
808 | xyz[ id1 ] = fBMin[id1] + sclA[it]*( fBMax[id1]-fBMin[id1] ); | |
809 | xyz[ id2 ] = fBMin[id2] + sclB[it]*( fBMax[id2]-fBMin[id2] ); | |
810 | // | |
811 | int* npt = GetNCNeeded(xyz,idim, dimMN,dimMX, Prec); // npoints for Bx,By,Bz | |
812 | for (int ib=0;ib<3;ib++) if (npt[ib]>gridBC[ib][idim]) gridBC[ib][idim] = npt[ib]+2; | |
813 | // | |
814 | } | |
815 | } | |
816 | } | |
817 | ||
818 | ||
819 | int* AliCheb3D::GetNCNeeded(float xyz[3],int DimVar, float mn,float mx, float prec) | |
820 | { | |
821 | // estimate needed number of chebyshev coefs for given function desctiption in DimVar dimension | |
822 | // The values for two other dimensions must be set beforehand | |
823 | // | |
824 | static int curNC[3]; | |
825 | static int retNC[3]; | |
826 | const int kMaxPoint = 400; | |
827 | float* gridVal = new float[3*kMaxPoint]; | |
828 | float* coefs = new float[3*kMaxPoint]; | |
829 | // | |
830 | float scale = mx-mn; | |
831 | float offs = mn + scale/2.0; | |
832 | scale = 2./scale; | |
833 | // | |
834 | int curNP; | |
835 | int maxNC=-1; | |
836 | int maxNCPrev=-1; | |
837 | for (int i=0;i<3;i++) retNC[i] = -1; | |
838 | for (int i=0;i<3;i++) fArgsTmp[i] = xyz[i]; | |
839 | // | |
840 | for (curNP=5; curNP<kMaxPoint; curNP+=5) { | |
841 | maxNCPrev = maxNC; | |
842 | // | |
843 | for (int i=0;i<curNP;i++) { // get function values on Cheb. nodes | |
844 | float x = TMath::Cos( TMath::Pi()*(i+0.5)/curNP ); | |
845 | fArgsTmp[DimVar] = x/scale+offs; // map to requested interval | |
846 | EvalUsrFunction(); | |
847 | for (int ib=3;ib--;) gridVal[ib*kMaxPoint + i] = fResTmp[ib]; | |
848 | } | |
849 | // | |
850 | for (int ib=0;ib<3;ib++) { | |
851 | curNC[ib] = AliCheb3D::CalcChebCoefs(&gridVal[ib*kMaxPoint], curNP, &coefs[ib*kMaxPoint],prec); | |
852 | if (maxNC < curNC[ib]) maxNC = curNC[ib]; | |
853 | if (retNC[ib] < curNC[ib]) retNC[ib] = curNC[ib]; | |
854 | } | |
855 | if ( (curNP-maxNC)>3 && (maxNC-maxNCPrev)<1 ) break; | |
856 | maxNCPrev = maxNC; | |
857 | // | |
858 | } | |
859 | delete[] gridVal; | |
860 | delete[] coefs; | |
861 | return retNC; | |
862 | // | |
863 | } | |
864 | #endif |