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 "AliCheb2DStackF.h"
18 : #include "AliLog.h"
19 : #include <TMath.h>
20 :
21 176 : ClassImp(AliCheb2DStackF)
22 :
23 : //____________________________________________________________________
24 : AliCheb2DStackF::AliCheb2DStackF()
25 0 : :AliCheb2DStack()
26 0 : ,fCoeffs(0)
27 0 : {
28 : // Default constructor
29 0 : }
30 :
31 : //____________________________________________________________________
32 : AliCheb2DStackF::AliCheb2DStackF(stFun_t fun, int nSlices, int dimOut, const float bmin[2],const float bmax[2],
33 : const int np[2], const float* dead, const float *rowXI, const float* precD)
34 0 : :AliCheb2DStack(nSlices,dimOut,bmin,bmax,dead,rowXI)
35 0 : ,fCoeffs(0)
36 0 : {
37 : // create stack of 2D->dimOut Chebyshev parameterizations debined in 2 dimensions between bmin and bmax,
38 : // and trained with function fun on 2D grid on np points.
39 : // Truncate each precision of each output dimension parameterization i to precD[i] if precD!=0, or prec
40 : //
41 0 : int *npd = new int[2*dimOut];
42 0 : int maxcoefs=np[0]*np[1]*fNSlices, maxrows=np[0]*fNSlices;
43 0 : for (int i=0;i<dimOut;i++) {
44 0 : npd[2*i] = np[0];
45 0 : npd[2*i+1] = np[1];
46 : }
47 0 : CheckDimensions(npd); // basic check
48 0 : fNCols = new UChar_t[maxrows];
49 0 : fCoeffs = new Float_t[maxcoefs];
50 0 : CreateParams(fun, npd, precD);
51 0 : delete[] npd;
52 : //
53 0 : Print();
54 : //
55 0 : }
56 :
57 : //____________________________________________________________________
58 : AliCheb2DStackF::AliCheb2DStackF(stFun_t fun, int nSlices, int dimOut,
59 : const float bmin[2],const float bmax[2],
60 : const int np[][2], const float* dead, const float *rowXI, const float* precD)
61 0 : :AliCheb2DStack(nSlices,dimOut,bmin,bmax,dead,rowXI)
62 0 : ,fCoeffs(0)
63 0 : {
64 : // create stack of 2D->dimOut Chebyshev parameterizations debined in 2 dimensions between bmin and bmax,
65 : // and trained with function fun on 2D grid on np[i][2] points for i-th output dimension.
66 : // Truncate each precision of each output dimension parameterization i to precD[i] if precD!=0, or prec
67 : //
68 : int maxcoefs=0, maxrows=0;
69 0 : for (int id=fDimOut;id--;) {
70 0 : maxcoefs += np[id][0]*np[id][1];
71 0 : maxrows += np[id][0];
72 : }
73 0 : fNCols = new UChar_t[maxrows*fNSlices];
74 0 : fCoeffs = new Float_t[maxcoefs*fNSlices];
75 0 : CreateParams(fun, (const int*)np, precD);
76 : //
77 0 : Print();
78 0 : }
79 :
80 : //____________________________________________________________________
81 : AliCheb2DStackF::~AliCheb2DStackF()
82 0 : {
83 : // D-tor
84 0 : delete[] fCoeffs;
85 0 : }
86 :
87 : //____________________________________________________________________
88 : void AliCheb2DStackF::Eval(int sliceID, const float *par, float *res) const
89 : {
90 : // evaluate Chebyshev parameterization for 2d->DimOut function at sliceID
91 0 : float p0,p1;
92 0 : MapToInternal(sliceID, par,p0,p1);
93 0 : const UChar_t *rows = &fNRows[sliceID*fDimOut]; // array of fDimOut rows for current slice
94 0 : const UChar_t *cols = &fNCols[fColEntry[sliceID]]; // array of columns per row for current slice
95 0 : const Float_t *cfs = &fCoeffs[fCoeffsEntry[sliceID]]; // array of coefficients for current slice
96 0 : for (int id=0;id<fDimOut;id++) {
97 0 : int nr = *rows++; // N rows in the matrix of coeffs for given dimension
98 0 : for (int ir=0;ir<nr;ir++) {
99 0 : int nc = *cols++; // N of significant colums at this row
100 0 : fWSpace[ir] = ChebEval1D(p1,cfs,nc); // interpolation of Cheb. coefs along row
101 0 : cfs += nc; // prepare coefs for the next row
102 : }
103 0 : res[id] = ChebEval1D(p0,fWSpace,nr);
104 : }
105 : //
106 0 : }
107 :
108 : //____________________________________________________________________
109 : Float_t AliCheb2DStackF::Eval(int sliceID, int dimOut, const float *par) const
110 : {
111 : // evaluate Chebyshev parameterization for requested output dimension only at requested sliceID
112 0 : float p0,p1;
113 0 : MapToInternal(sliceID,par,p0,p1);
114 0 : int pid = sliceID*fDimOut;
115 0 : const UChar_t *rows = &fNRows[pid]; // array of fDimOut rows for current slice
116 0 : const UChar_t *cols = &fNCols[fColEntry[sliceID]]; // array of columns per row for current slice
117 0 : const Float_t *cfs = &fCoeffs[fCoeffsEntry[sliceID]]; // array of coefficients for current slice
118 0 : while (dimOut) {
119 0 : for (int ir=*rows++;ir--;) cfs += *cols++; // go to the matrix of needed row
120 0 : dimOut--;
121 : }
122 0 : int nr = *rows++; // N rows in the matrix of coeffs for given dimension
123 0 : for (int ir=0;ir<nr;ir++) {
124 0 : int nc = *cols++; // N of significant colums at this row
125 0 : fWSpace[ir] = ChebEval1D(p1,cfs,nc); // interpolation of Cheb. coefs along row
126 0 : cfs += nc; // prepare coefs for the next row
127 : }
128 0 : return ChebEval1D(p0,fWSpace,nr);
129 : //
130 0 : }
131 :
132 : //____________________________________________________________________
133 : void AliCheb2DStackF::CreateParams(stFun_t fun, const int *np, const float* prc)
134 : {
135 : // create parameterizations
136 : //
137 : // temporary space for max possible coeffs, rows etc
138 0 : float **grids = new float*[fDimOut]; // Chebyshev grids for each output dimension
139 : int maxSpace = 1, totSpace = 1;
140 : Bool_t sameGrid = kTRUE;
141 0 : int ref0=np[0],ref1=np[1];
142 0 : for (int id=fDimOut;id--;) {
143 0 : int nsp = np[2*id]*np[2*id+1];
144 0 : if (maxSpace<nsp) maxSpace = nsp;
145 0 : totSpace += nsp;
146 0 : if (ref0!=np[2*id] || ref1!=np[2*id+1]) sameGrid = kFALSE;
147 : }
148 : //
149 : // save pointers to recover the beggining of arrays later
150 0 : UChar_t* nRows0 = fNRows;
151 0 : UChar_t* nCols0 = fNCols;
152 0 : float* coeffs0 = fCoeffs;
153 : //
154 0 : float *tmpCoef2D = new float[maxSpace]; // temporary workspace
155 0 : float *tmpVals = new float[totSpace]; // temporary workspace for function values
156 : //
157 0 : for (int isl=0;isl<fNSlices;isl++) {
158 0 : for (int id=fDimOut;id--;) grids[id] = DefineGrid(isl, id, &np[2*id]);
159 0 : fCoeffsEntry[isl] = fCoeffs - coeffs0; // offset of the given slice coeffs
160 0 : fColEntry[isl] = fNCols - nCols0; // offset of the given slice columns dimensions
161 : //
162 0 : if (sameGrid) FillFunValues(fun, isl, grids[0], &np[0], tmpVals);
163 : else {
164 : int slot = 0;
165 0 : for (int id=0;id<fDimOut;id++) {
166 0 : FillFunValues(fun, isl, id, grids[id], &np[2*id], tmpVals+slot); // get values for single dimensions
167 0 : slot += np[2*id]*np[2*id+1];
168 : }
169 : }
170 : //
171 : int slot = 0;
172 0 : for (int id=0;id<fDimOut;id++) {
173 0 : ChebFit(&np[2*id], tmpVals+slot, tmpCoef2D, prc ? prc[id] : fgkDefPrec);
174 0 : slot += np[2*id]*np[2*id+1];
175 : }
176 0 : for (int id=fDimOut;id--;) delete[] grids[id];
177 : }
178 0 : delete[] grids;
179 0 : delete[] tmpCoef2D;
180 : //
181 0 : fNCoefsTot = fCoeffs-coeffs0; // size of final coeffs array
182 : //
183 0 : fNRows = nRows0;
184 : // redefine arrays in compressed form, clean temp. stuff
185 0 : fNCols = new UChar_t[fNRowsTot];
186 0 : memcpy(fNCols, nCols0, fNRowsTot*sizeof(UChar_t));
187 0 : delete[] nCols0;
188 0 : fCoeffs = new Float_t[fNCoefsTot];
189 0 : memcpy(fCoeffs, coeffs0, fNCoefsTot*sizeof(Float_t));
190 0 : delete[] coeffs0;
191 : //
192 0 : }
193 :
194 : //____________________________________________________________________
195 : void AliCheb2DStackF::ChebFit(const int np[2], const float* tmpVals, float* tmpCoef2D, float prec)
196 : {
197 : // prepare Cheb.fit for 2D -> single output dimension
198 : //
199 0 : int ncmax=0, maxDim = TMath::Max(np[0],np[1]);
200 0 : memset(tmpCoef2D,0,np[0]*np[1]*sizeof(float));
201 : //
202 0 : float rTiny = 0.5*prec/maxDim; // neglect coefficient below this threshold
203 : //
204 0 : for (int id1=np[1];id1--;) { // create Cheb.param for each node of 2nd input dimension
205 0 : int nc = CalcChebCoefs(&tmpVals[id1*np[0]], np[0], fWSpace, -1);
206 0 : for (int id0=nc;id0--;) tmpCoef2D[id1 + id0*np[1]] = fWSpace[id0];
207 0 : if (ncmax<nc) ncmax = nc; // max coefs to be kept in dim0 to guarantee needed precision
208 : }
209 : //
210 : // once each 1d slice of given 2d slice is parametrized, parametrize the Cheb.coeffs
211 0 : for (int id0=np[0];id0--;) {
212 0 : CalcChebCoefs(&tmpCoef2D[id0*np[1]], np[1], fWSpace, -1);
213 0 : for (int id1=np[1];id1--;) tmpCoef2D[id1+id0*np[1]] = fWSpace[id1];
214 : }
215 : //
216 : // now find 1D curve which separates significant coefficients of 2D matrix from nonsignificant ones (up to prec)
217 : // double resid = 0;
218 0 : int ncNZero=0, nRows = np[0]; // find max significant row
219 0 : for (int id0=np[0];id0--;) {
220 0 : fNCols[id0]=0;
221 : double resid = 0;
222 0 : for (int id1=np[1];id1--;) {
223 0 : int id = id1 + id0*np[1];
224 0 : float cfa = TMath::Abs(tmpCoef2D[id]);
225 0 : if (cfa < rTiny) {tmpCoef2D[id] = 0; continue;} // neglect coefs below the threshold
226 0 : resid += cfa;
227 : // if (resid<prec) continue; // this coeff is negligible
228 0 : if (resid<rTiny) continue; // this coeff is negligible
229 : // resid -= cfa; // otherwise go back 1 step
230 0 : fNCols[id0] = id1+1; // how many coefs to keep
231 0 : break;
232 : }
233 0 : if (fNCols[id0]) ncNZero++;
234 0 : else if (!ncNZero) nRows--; // decrease N significant rows untile 1st non-0 column met
235 : }
236 : //
237 : // find max significant column and fill the storage
238 : // for the max sigificant column of each row
239 0 : for (int id0=0;id0<nRows;id0++) {
240 0 : int nc = *fNCols++;
241 0 : for (int id1=0;id1<nc;id1++) *fCoeffs++ = tmpCoef2D[id1+id0*np[1]];
242 : }
243 0 : *fNRows++ = nRows;
244 0 : fNRowsTot += nRows;
245 : //
246 0 : }
247 :
248 : //____________________________________________________________________
249 : void AliCheb2DStackF::FillFunValues(stFun_t fun, int slice, int dim, const float *grid,
250 : const int np[2], float* vals)
251 : {
252 : // fill function values on the grid
253 0 : float args[2];
254 0 : for (int id1=np[1];id1--;) {
255 0 : args[1] = grid[id1];
256 0 : for (int id0=np[0];id0--;) { //
257 0 : args[0] = grid[np[1]+id0];
258 0 : fun(slice, args, fWSpace);
259 0 : vals[id1*np[0] + id0] = fWSpace[dim];
260 : }
261 : }
262 0 : }
263 :
264 : //____________________________________________________________________
265 : void AliCheb2DStackF::FillFunValues(stFun_t fun, int slice, const float *grid, const int np[2], float* vals)
266 : {
267 : // fill function values on the grid for all dimensions at once (if grids are the same)
268 0 : float args[2];
269 0 : int slotStep = np[0]*np[1];
270 0 : for (int id1=np[1];id1--;) {
271 0 : args[1] = grid[id1];
272 0 : for (int id0=np[0];id0--;) { //
273 0 : args[0] = grid[np[1]+id0];
274 0 : fun(slice, args, fWSpace);
275 : int slotDim = 0;
276 0 : for (int dim=0;dim<fDimOut;dim++) {
277 0 : vals[slotDim + id1*np[0] + id0] = fWSpace[dim];
278 0 : slotDim += slotStep;
279 : }
280 : }
281 : }
282 0 : }
283 :
284 : //__________________________________________________________________________________________
285 : void AliCheb2DStackF::Print(const Option_t* opt) const
286 : {
287 0 : printf("F*");
288 0 : AliCheb2DStack::Print(opt);
289 0 : }
290 :
291 : //__________________________________________________________________________________________
292 : void AliCheb2DStackF::PrintSlice(int isl, const Option_t* opt) const
293 : {
294 : // print slice info
295 : //
296 0 : TString opts = opt; opts.ToLower();
297 0 : Bool_t showcf = opts.Contains("c");
298 0 : int par0 = isl*fDimOut;
299 0 : const UChar_t* rows = &fNRows[par0];
300 0 : const UChar_t* cols = &fNCols[fColEntry[isl]];
301 0 : const Float_t* cfs = &fCoeffs[fCoeffsEntry[isl]];
302 0 : printf("#%3d ",isl);
303 0 : if (showcf) printf("\n");
304 0 : for (int id=0;id<fDimOut;id++) {
305 0 : int nr = *rows++;
306 : int ncmax=0,ncf=0;
307 0 : for (int ir=0;ir<nr;ir++) {
308 0 : ncf += cols[ir];
309 0 : if (cols[ir]>ncmax) ncmax = cols[ir];
310 : }
311 0 : printf("D%d: %4d coefs in %3dx%3d| ",id,ncf,nr,ncmax);
312 0 : if (showcf) {
313 0 : printf("\n");
314 0 : for (int ir=0;ir<nr;ir++) {
315 0 : for (int ic=0;ic<cols[ir];ic++) printf("%+.2e ",*cfs++); printf("\n");
316 : }
317 0 : }
318 0 : cols += nr; // cols entry for next dimension
319 : //
320 : }
321 0 : if (!showcf) printf("\n");
322 : //
323 0 : }
|
