Line data Source code
1 : /**************************************************************************
2 : * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 : * *
4 : * Author: The ALICE Off-line Project. *
5 : * Contributors are mentioned in the code where appropriate. *
6 : * *
7 : * Permission to use, copy, modify and distribute this software and its *
8 : * documentation strictly for non-commercial purposes is hereby granted *
9 : * without fee, provided that the above copyright notice appears in all *
10 : * copies and that both the copyright notice and this permission notice *
11 : * appear in the supporting documentation. The authors make no claims *
12 : * about the suitability of this software for any purpose. It is *
13 : * provided "as is" without express or implied warranty. *
14 : **************************************************************************/
15 :
16 :
17 : #include "AliCheb2DStackS.h"
18 : #include "AliLog.h"
19 : #include <TMath.h>
20 : #include <limits.h>
21 :
22 176 : ClassImp(AliCheb2DStackS)
23 :
24 : //____________________________________________________________________
25 : AliCheb2DStackS::AliCheb2DStackS()
26 0 : : AliCheb2DStack()
27 0 : ,fParScale(0)
28 0 : ,fParHVar(0)
29 0 : ,fCoeffs(0)
30 0 : {
31 : // Default constructor
32 0 : }
33 :
34 : //____________________________________________________________________
35 : AliCheb2DStackS::AliCheb2DStackS(stFun_t fun, int nSlices, int dimOut,
36 : const float bmin[2],const float bmax[2],
37 : const int np[2], const float* dead, const float *rowXI,
38 : const float* precD)
39 0 : : AliCheb2DStack(nSlices,dimOut,bmin,bmax,dead,rowXI)
40 0 : ,fParScale(new Float_t[nSlices*dimOut])
41 0 : ,fParHVar(new Float_t[nSlices*dimOut])
42 0 : ,fCoeffs(0)
43 0 : {
44 : // create stack of 2D->dimOut Chebyshev parameterizations debined in 2 dimensions between bmin and bmax,
45 : // and trained with function fun on 2D grid on np points.
46 : // Truncate each precision of each output dimension parameterization i to precD[i] if precD!=0, or prec
47 : //
48 0 : int *npd = new int[2*dimOut];
49 0 : int maxcoefs=np[0]*np[1]*fNSlices, maxrows=np[0]*fNSlices;
50 0 : for (int i=0;i<dimOut;i++) {
51 0 : npd[2*i] = np[0];
52 0 : npd[2*i+1] = np[1];
53 : }
54 0 : CheckDimensions(npd); // basic check
55 0 : fNCols = new UChar_t[maxrows];
56 0 : fCoeffs = new Short_t[maxcoefs];
57 0 : CreateParams(fun, npd, precD);
58 0 : delete[] npd;
59 : //
60 0 : Print();
61 : //
62 0 : }
63 :
64 : //____________________________________________________________________
65 : AliCheb2DStackS::AliCheb2DStackS(stFun_t fun, int nSlices, int dimOut,
66 : const float bmin[2],const float bmax[2],
67 : const int np[][2], const float* dead, const float *rowXI,
68 : const float* precD)
69 0 : :AliCheb2DStack(nSlices,dimOut,bmin,bmax,dead,rowXI)
70 0 : ,fParScale(new Float_t[nSlices*dimOut])
71 0 : ,fParHVar(new Float_t[nSlices*dimOut])
72 0 : ,fCoeffs(0)
73 0 : {
74 : // create stack of 2D->dimOut Chebyshev parameterizations debined in 2 dimensions between bmin and bmax,
75 : // and trained with function fun on 2D grid on np[i][2] points for i-th output dimension.
76 : // Truncate each precision of each output dimension parameterization i to precD[i] if precD!=0, or prec
77 : //
78 : int maxcoefs=0, maxrows=0;
79 0 : for (int id=fDimOut;id--;) {
80 0 : maxcoefs += np[id][0]*np[id][1];
81 0 : maxrows += np[id][0];
82 : }
83 0 : fNCols = new UChar_t[maxrows*fNSlices];
84 0 : fCoeffs = new Short_t[maxcoefs*fNSlices];
85 0 : CreateParams(fun, (const int*)np, precD);
86 0 : Print();
87 : //
88 0 : }
89 :
90 : //____________________________________________________________________
91 : AliCheb2DStackS::~AliCheb2DStackS()
92 0 : {
93 : // D-tor
94 0 : delete[] fCoeffs;
95 0 : delete[] fParScale;
96 0 : delete[] fParHVar;
97 0 : }
98 :
99 : //____________________________________________________________________
100 : void AliCheb2DStackS::Eval(int sliceID, const float *par, float *res) const
101 : {
102 : // evaluate Chebyshev parameterization for 2d->DimOut function at sliceID
103 0 : float p0,p1;
104 0 : MapToInternal(sliceID,par,p0,p1);
105 0 : int pid = sliceID*fDimOut;
106 0 : const UChar_t *rows = &fNRows[pid]; // array of fDimOut rows for current slice
107 0 : const UChar_t *cols = &fNCols[fColEntry[sliceID]]; // array of columns per row for current slice
108 0 : const Short_t *cfs = &fCoeffs[fCoeffsEntry[sliceID]]; // array of coefficients for current slice
109 0 : const Float_t *scl = &fParScale[pid];
110 0 : const Float_t *hvr = &fParHVar[pid];
111 0 : for (int id=0;id<fDimOut;id++) {
112 0 : int nr = *rows++; // N rows in the matrix of coeffs for given dimension
113 0 : for (int ir=0;ir<nr;ir++) {
114 0 : int nc = *cols++; // N of significant colums at this row
115 0 : fWSpace[ir] = ChebEval1D(p1,cfs,nc); // interpolation of Cheb. coefs along row
116 0 : cfs += nc; // prepare coefs for the next row
117 : }
118 0 : res[id] = ChebEval1D(p0,fWSpace,nr) * (*scl++) + (*hvr++);
119 : }
120 : //
121 0 : }
122 :
123 : //____________________________________________________________________
124 : Float_t AliCheb2DStackS::Eval(int sliceID, int dimOut, const float *par) const
125 : {
126 : // evaluate Chebyshev parameterization for requested output dimension only at requested sliceID
127 0 : float p0,p1;
128 0 : MapToInternal(sliceID,par,p0,p1);
129 0 : int pid = sliceID*fDimOut;
130 0 : const UChar_t *rows = &fNRows[pid]; // array of fDimOut rows for current slice
131 0 : const UChar_t *cols = &fNCols[fColEntry[sliceID]]; // array of columns per row for current slice
132 0 : const Short_t *cfs = &fCoeffs[fCoeffsEntry[sliceID]]; // array of coefficients for current slice
133 0 : const Float_t *scl = &fParScale[pid];
134 0 : const Float_t *hvr = &fParHVar[pid];
135 0 : while (dimOut) {
136 0 : for (int ir=*rows++;ir--;) cfs += *cols++; // go to the matrix of needed row
137 0 : dimOut--;
138 0 : scl++;
139 0 : hvr++;
140 : }
141 0 : int nr = *rows++; // N rows in the matrix of coeffs for given dimension
142 0 : for (int ir=0;ir<nr;ir++) {
143 0 : int nc = *cols++; // N of significant colums at this row
144 0 : fWSpace[ir] = ChebEval1D(p1,cfs,nc); // interpolation of Cheb. coefs along row
145 0 : cfs += nc; // prepare coefs for the next row
146 : }
147 0 : return ChebEval1D(p0,fWSpace,nr) * (*scl) + (*hvr);
148 : //
149 0 : }
150 :
151 : //____________________________________________________________________
152 : void AliCheb2DStackS::CreateParams(stFun_t fun, const int *np, const float* prc)
153 : {
154 : // create parameterizations
155 : //
156 : // temporary space for max possible coeffs, rows etc
157 0 : float **grids = new float*[fDimOut]; // Chebyshev grids for each output dimension
158 : int maxSpace = 1, totSpace = 1;
159 : Bool_t sameGrid = kTRUE;
160 0 : int ref0=np[0],ref1=np[1];
161 0 : for (int id=fDimOut;id--;) {
162 0 : int nsp = np[2*id]*np[2*id+1];
163 0 : if (maxSpace<nsp) maxSpace = nsp;
164 0 : totSpace += nsp;
165 0 : if (ref0!=np[2*id] || ref1!=np[2*id+1]) sameGrid = kFALSE;
166 : }
167 : //
168 : // save pointers to recover the beggining of arrays later
169 0 : UChar_t* nRows0 = fNRows;
170 0 : UChar_t* nCols0 = fNCols;
171 0 : Short_t* coeffs0 = fCoeffs;
172 : //
173 0 : float *tmpCoef2D = new float[maxSpace]; // temporary workspace for coeffs
174 0 : float *tmpVals = new float[totSpace]; // temporary workspace for function values
175 : //
176 0 : for (int isl=0;isl<fNSlices;isl++) {
177 0 : for (int id=fDimOut;id--;) grids[id] = DefineGrid(isl, id, &np[2*id]);
178 0 : fCoeffsEntry[isl] = fCoeffs - coeffs0; // offset of the given slice coeffs
179 0 : fColEntry[isl] = fNCols - nCols0; // offset of the given slice columns dimensions
180 : //
181 0 : if (sameGrid) FillFunValues(fun, isl, grids[0], &np[0], tmpVals);
182 : else {
183 : int slot = 0;
184 0 : for (int id=0;id<fDimOut;id++) {
185 0 : FillFunValues(fun, isl, id, grids[id], &np[2*id], tmpVals+slot); // get values for single dimensions
186 0 : slot += np[2*id]*np[2*id+1];
187 : }
188 : }
189 : int slot = 0;
190 0 : for (int id=0;id<fDimOut;id++) {
191 0 : float prcs = prc ? prc[id] : fgkDefPrec;
192 0 : fParScale[isl*fDimOut+id] = ChebFit(&np[2*id], tmpVals+slot, tmpCoef2D, prcs);
193 0 : slot += np[2*id]*np[2*id+1];
194 : }
195 0 : for (int id=fDimOut;id--;) delete[] grids[id];
196 : }
197 0 : delete[] grids;
198 0 : delete[] tmpCoef2D;
199 0 : delete[] tmpVals;
200 : //
201 0 : fNCoefsTot = fCoeffs-coeffs0; // size of final coeffs array
202 : //
203 0 : fNRows = nRows0;
204 : // redefine arrays in compressed form, clean temp. stuff
205 0 : fNCols = new UChar_t[fNRowsTot];
206 0 : memcpy(fNCols, nCols0, fNRowsTot*sizeof(UChar_t));
207 0 : delete[] nCols0;
208 0 : fCoeffs = new Short_t[fNCoefsTot];
209 0 : memcpy(fCoeffs, coeffs0, fNCoefsTot*sizeof(Short_t));
210 0 : delete[] coeffs0;
211 : //
212 0 : }
213 :
214 : //____________________________________________________________________
215 : Float_t AliCheb2DStackS::ChebFit(const int np[2], const float* tmpVals, float* tmpCoef2D, float prec)
216 : {
217 : // prepare Cheb.fit for 2D -> single output dimension
218 : //
219 0 : int maxDim = TMath::Max(np[0],np[1]);
220 0 : memset(tmpCoef2D,0,np[0]*np[1]*sizeof(float));
221 : //
222 0 : for (int id1=np[1];id1--;) { // create Cheb.param for each node of 2nd input dimension
223 0 : int nc = CalcChebCoefs(&tmpVals[id1*np[0]], np[0], fWSpace, -1);
224 0 : for (int id0=nc;id0--;) tmpCoef2D[id1 + id0*np[1]] = fWSpace[id0];
225 : }
226 : //
227 : // once each 1d slice of given 2d slice is parametrized, parametrize the Cheb.coeffs
228 : float mxAbs = -1; // we need largest coeff for rescaling to +-MaxShort range
229 0 : for (int id0=np[0];id0--;) {
230 0 : CalcChebCoefs(&tmpCoef2D[id0*np[1]], np[1], fWSpace, -1);
231 0 : for (int id1=np[1];id1--;) {
232 0 : tmpCoef2D[id1+id0*np[1]] = fWSpace[id1];
233 0 : if (TMath::Abs(fWSpace[id1])>mxAbs) mxAbs = TMath::Abs(fWSpace[id1]);
234 : }
235 : }
236 : //
237 0 : float scl = mxAbs>0 ? SHRT_MAX/mxAbs : 1;
238 : // rescale/truncate to +-MaxShort range
239 0 : for (int id0=np[0];id0--;) {
240 0 : for (int id1=np[1];id1--;) {
241 0 : int id = id1 + id0*np[1];
242 0 : Short_t cfs = Short_t(tmpCoef2D[id1+id0*np[1]]*scl);
243 0 : tmpCoef2D[id] = float(cfs);
244 : }
245 : }
246 : //
247 0 : float rTiny = 0.5*prec*scl/maxDim; // neglect coefficient below this threshold
248 : //
249 : // now find 1D curve which separates significant coefficients of 2D matrix
250 : // from nonsignificant ones (up to prec)
251 : // double resid = 0;
252 0 : int ncNZero=0, nRows = np[0]; // find max significant row
253 0 : for (int id0=np[0];id0--;) {
254 0 : fNCols[id0]=0;
255 : double resid = 0;
256 0 : for (int id1=np[1];id1--;) {
257 0 : int id = id1 + id0*np[1];
258 0 : float cfa = TMath::Abs(tmpCoef2D[id]);
259 0 : if (cfa < rTiny) {tmpCoef2D[id] = 0; continue;} // neglect coefs below the threshold
260 0 : resid += cfa;
261 : // if (resid<prec) continue; // this coeff is negligible
262 0 : if (resid<rTiny) continue; // this coeff is negligible
263 : // resid -= cfa; // otherwise go back 1 step
264 0 : fNCols[id0] = id1+1; // how many coefs to keep
265 0 : break;
266 : }
267 0 : if (fNCols[id0]) ncNZero++;
268 0 : else if (!ncNZero) nRows--; // decrease N significant rows untile 1st non-0 column met
269 : }
270 : //
271 : // find max significant column and fill the coeffs in the storage
272 0 : for (int id0=0;id0<nRows;id0++) {
273 0 : int nc = *fNCols++;
274 0 : for (int id1=0;id1<nc;id1++) {
275 0 : *fCoeffs++ = Short_t(tmpCoef2D[id1+id0*np[1]]);
276 : }
277 : }
278 0 : *fNRows++ = nRows;
279 0 : fNRowsTot += nRows;
280 : //
281 0 : return 1./scl;
282 : }
283 :
284 : //____________________________________________________________________
285 : void AliCheb2DStackS::FillFunValues(stFun_t fun, int slice, int dim, const float *grid,
286 : const int np[2], float* vals)
287 : {
288 : // fill function values on the grid
289 0 : float args[2];
290 : float minv=1e15,maxv=-1e15;
291 0 : for (int id1=np[1];id1--;) {
292 0 : args[1] = grid[id1];
293 0 : for (int id0=np[0];id0--;) { //
294 0 : args[0] = grid[np[1]+id0];
295 0 : fun(slice, args, fWSpace);
296 0 : vals[id1*np[0] + id0] = fWSpace[dim];
297 0 : if (minv>fWSpace[dim]) minv = fWSpace[dim];
298 0 : if (maxv<fWSpace[dim]) maxv = fWSpace[dim];
299 : }
300 : }
301 0 : fParHVar[slice*fDimOut+dim] = (maxv+minv)/2.;
302 : //
303 : // map values to +- fWSpace
304 0 : float offs = fParHVar[slice*fDimOut+dim];
305 0 : for (int id1=np[1];id1--;) for (int id0=np[0];id0--;) vals[id1*np[0] + id0] -= offs;
306 : //
307 0 : }
308 :
309 : //____________________________________________________________________
310 : void AliCheb2DStackS::FillFunValues(stFun_t fun, int slice, const float *grid,
311 : const int np[2], float* vals)
312 : {
313 : // fill function values on the grid for all output dimension at same time (same grid)
314 0 : float args[2];
315 0 : float *minv = new float[fDimOut], *maxv = new float[fDimOut];
316 0 : for (int dim=fDimOut;dim--;) {
317 0 : minv[dim] = 1e15;
318 0 : maxv[dim] =-1e15;
319 : }
320 0 : int stepDim = np[0]*np[1];
321 0 : for (int id1=np[1];id1--;) {
322 0 : args[1] = grid[id1];
323 0 : for (int id0=np[0];id0--;) { //
324 0 : args[0] = grid[np[1]+id0];
325 0 : fun(slice, args, fWSpace);
326 0 : for (int dim=fDimOut;dim--;) {
327 0 : vals[stepDim*dim + id1*np[0] + id0] = fWSpace[dim];
328 0 : if (minv[dim]>fWSpace[dim]) minv[dim] = fWSpace[dim];
329 0 : if (maxv[dim]<fWSpace[dim]) maxv[dim] = fWSpace[dim];
330 : }
331 : }
332 : }
333 0 : for (int dim=fDimOut;dim--;) {
334 0 : float offs = fParHVar[slice*fDimOut+dim] = (maxv[dim]+minv[dim])/2.;
335 : // map values to +- fWSpace
336 0 : float* valsdim = vals+stepDim*dim;
337 0 : for (int id1=np[1];id1--;) for (int id0=np[0];id0--;) valsdim[id1*np[0] + id0] -= offs;
338 : }
339 : //
340 0 : delete[] minv;
341 0 : delete[] maxv;
342 : //
343 0 : }
344 :
345 : //__________________________________________________________________________________________
346 : void AliCheb2DStackS::Print(const Option_t* opt) const
347 : {
348 0 : printf("S*");
349 0 : AliCheb2DStack::Print(opt);
350 0 : }
351 :
352 : //__________________________________________________________________________________________
353 : void AliCheb2DStackS::PrintSlice(int isl, const Option_t* opt) const
354 : {
355 : // print slice info
356 : //
357 0 : TString opts = opt; opts.ToLower();
358 0 : Bool_t showcf = opts.Contains("c");
359 0 : int par0 = isl*fDimOut;
360 0 : const UChar_t* rows = &fNRows[par0];
361 0 : const UChar_t* cols = &fNCols[fColEntry[isl]];
362 0 : const Short_t* cfs = &fCoeffs[fCoeffsEntry[isl]];
363 0 : printf("#%3d ",isl);
364 0 : if (showcf) printf("\n");
365 0 : for (int id=0;id<fDimOut;id++) {
366 0 : int nr = *rows++;
367 : int ncmax=0,ncf=0;
368 0 : for (int ir=0;ir<nr;ir++) {
369 0 : ncf += cols[ir];
370 0 : if (cols[ir]>ncmax) ncmax = cols[ir];
371 : }
372 0 : printf("D%d: %4d coefs in %3dx%3d (Scl:%.2e Hv:%.2e)| ",
373 0 : id,ncf,nr,ncmax,fParScale[isl*fDimOut+id],fParHVar[isl*fDimOut+id]);
374 0 : if (showcf) {
375 0 : printf("\n");
376 0 : for (int ir=0;ir<nr;ir++) {
377 0 : for (int ic=0;ic<cols[ir];ic++) printf("%+6d ",*cfs++); printf("\n");
378 : }
379 0 : }
380 0 : cols += nr; // cols entry for next dimension
381 : //
382 : }
383 0 : if (!showcf) printf("\n");
384 : //
385 0 : }
|
