Project
Loading...
Searching...
No Matches
poissonEM.cxx
Go to the documentation of this file.
1// Copyright 2019-2020 CERN and copyright holders of ALICE O2.
2// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders.
3// All rights not expressly granted are reserved.
4//
5// This software is distributed under the terms of the GNU General Public
6// License v3 (GPL Version 3), copied verbatim in the file "COPYING".
7//
8// In applying this license CERN does not waive the privileges and immunities
9// granted to it by virtue of its status as an Intergovernmental Organization
10// or submit itself to any jurisdiction.
11
12#include <cfloat>
13#include <cmath>
14
15#include <gsl/gsl_blas.h>
16#include <gsl/gsl_vector.h>
17#include <gsl/gsl_matrix.h>
18
19#include "MCHClustering/ClusterConfig.h"
20#include "MCHClustering/PadsPEM.h"
21#include "mathUtil.h"
22#include "mathieson.h"
23
24//
25// TODO : Optimization, generateMixedGaussians2D computed twice.
26//
27
28namespace o2
29{
30namespace mch
31{
32
33extern ClusterConfig clusterConfig;
34
35void iterateEMPoisson(const double* Cij, const double* Ci,
36 const Mask_t* maskCij, const double* qPixels,
37 const double* qPad, double* qPadPrediction, int nPixels,
38 int nPads, double* newQPixels)
39{
40
41 double residu = 0.0;
42
43 // Compute charge prediction on pad j based on pixel charges
44 // qPadPrediction[j] = Sum_i{ Cij[i,j].qPixels[i] }
45 for (int j = 0; j < nPads; j++) {
46 qPadPrediction[j] = 0;
47 for (int i = 0; i < nPixels; i++) {
48 qPadPrediction[j] +=
49 // maskCij[nPads * i + j] * Cij[nPads * i + j] * qPixels[i];
50 Cij[nPads * i + j] * qPixels[i];
51 }
52 // Prevent zero division
53 if (qPadPrediction[j] < 1.0e-6) {
54 if (qPad[j] < 1.0e-6) {
55 qPadPrediction[j] = 1.0;
56 } else {
57 qPadPrediction[j] = 0.1 * qPad[j];
58 }
59 }
60 residu += fabs(qPadPrediction[j] - qPad[j]);
61 }
62
63 // Update the pixel charge qPixels with the
64 // the new predicted charges qPadPrediction
65 // qPixels[i] = qPixels[i] / Ci[i] * Sum_j { Cij[i,j]*qPad[j] /
66 // qPadPrediction[j] }
67 for (int i = 0; i < nPixels; i++) {
68 // Normalization term
69 // r = np.sum( Cij[:,j]*qPad[0:nPads] / qPadPrediction[0:nPads] )
70 if (Ci[i] > 1.0e-10) {
71 double s_i = 0;
72 for (int j = 0; j < nPads; j++) {
73 // s_i += maskCij[nPads * i + j] * Cij[nPads * i + j] * qPad[j] /
74 s_i += Cij[nPads * i + j] * qPad[j] /
75 qPadPrediction[j];
76 }
77 newQPixels[i] = s_i * qPixels[i] / Ci[i];
78 } else {
79 newQPixels[i] = 0;
80 }
81 }
82}
83
84void fastIterateEMPoissonV0(const double* Cij, const double* Ci,
85 const double* qPixels, const double* qPad,
86 double* qPadPrediction, int nPixels, int nPads,
87 double* newQPixels)
88{
89
90 double residu = 0.0;
91 double* qRatio = new double[nPads];
92
93 // Compute charge prediction on pad j based on pixel charges
94 // qPadPrediction[j] = Sum_i{ Cij[i,j].qPixels[i] }
95 for (int j = 0; j < nPads; j++) {
96 qPadPrediction[j] = 0;
97 for (int i = 0; i < nPixels; i++) {
98 qPadPrediction[j] += Cij[nPads * i + j] * qPixels[i];
99 }
100 // Prevent zero division
101
102 if (qPadPrediction[j] < 1.0e-10) {
103 qPadPrediction[j] = qPadPrediction[j] + 1.0e-10;
104 }
105 /*
106 if (qPad[j] < 1.0e-6) {
107 qPadPrediction[j] = 1.0;
108 } else {
109 qPadPrediction[j] = 0.1 * qPad[j];
110 }
111 }
112 */
113 qRatio[j] = qPad[j] / qPadPrediction[j];
114 residu += fabs(qPadPrediction[j] - qPad[j]);
115 }
116
117 // Update the pixel charge qPixels with the
118 // the new predicted charges qPadPrediction
119 // qPixels[i] = qPixels[i] / Ci[i] * Sum_j { Cij[i,j]*qPad[j] /
120 // qPadPrediction[j] }
121 for (int i = 0; i < nPixels; i++) {
122 // Normalization term
123 // r = np.sum( Cij[:,j]*qPad[0:nPads] / qPadPrediction[0:nPads] )
124 if (Ci[i] > 1.0e-10) {
125 double s_i = 0;
126 for (int j = 0; j < nPads; j++) {
127 s_i += Cij[nPads * i + j] * qRatio[j];
128 }
129 newQPixels[i] = s_i * qPixels[i] / Ci[i];
130 } else {
131 newQPixels[i] = 0;
132 }
133 }
134 delete[] qRatio;
135}
136
137void fastIterateEMPoisson(const double* Cij, const double* Ci,
138 const double* qPixels, const double* qPad,
139 double* qPadPrediction, int nPixels, int nPads,
140 double* newQPixels)
141{
142
143 double residu = 0.0;
144 double* qRatio = new double[nPads];
145
146 // Compute charge prediction on pad j based on pixel charges
147 // qPadPrediction[j] = Sum_i{ Cij[i,j].qPixels[i] }
148 // printf("fastIterateEMPoisson Cij=%p, Ci=%p, qPixels=%p, qPad=%p, qPadPrediction=%p, nPixels=%d, nPads=%d, newQPixels=%p\n",
149 // Cij, Ci, qPixels, qPad, qPadPrediction, nPixels, nPads,newQPixels);
150 gsl_matrix_const_view Cij_gsl = gsl_matrix_const_view_array(Cij, nPixels, nPads);
151 gsl_vector_const_view qPixels_gsl = gsl_vector_const_view_array(qPixels, nPixels);
152 gsl_vector_view qPadPrediction_gsl = gsl_vector_view_array(qPadPrediction, nPads);
153
154 gsl_blas_dgemv(CblasTrans, 1.0, &Cij_gsl.matrix, &qPixels_gsl.vector, 0.0, &qPadPrediction_gsl.vector);
155 for (int j = 0; j < nPads; j++) {
156 // Prevent division by zero
157 if (qPadPrediction[j] < 1.0e-10) {
158 qPadPrediction[j] = qPadPrediction[j] + 1.0e-10;
159 }
160 qRatio[j] = qPad[j] / qPadPrediction[j];
161 residu += fabs(qPadPrediction[j] - qPad[j]);
162 }
163
164 // Update the pixel charge qPixels with the
165 // the new predicted charges qPadPrediction
166 // qPixels[i] = qPixels[i] / Ci[i] * Sum_j { Cij[i,j]*qPad[j] /
167 // qPadPrediction[j] }
168 gsl_vector_view qRatio_gsl = gsl_vector_view_array(qRatio, nPads);
169 gsl_vector_view newQPixels_gsl = gsl_vector_view_array(newQPixels, nPixels);
170 if (1) {
171 gsl_blas_dgemv(CblasNoTrans, 1.0, &Cij_gsl.matrix, &qRatio_gsl.vector, 0.0, &newQPixels_gsl.vector);
172
173 for (int i = 0; i < nPixels; i++) {
174 if (Ci[i] > 1.0e-10) {
175 newQPixels[i] = newQPixels[i] * qPixels[i] / Ci[i];
176 } else {
177 newQPixels[i] = 0;
178 }
179 }
180 } else if (0) {
181 for (int i = 0; i < nPixels; i++) {
182 // Normalization term
183 // r = np.sum( Cij[:,j]*qPad[0:nPads] / qPadPrediction[0:nPads] )
184 double s_i = 0;
185 newQPixels[i] = 0.;
186 for (int j = 0; j < nPads; j++) {
187 s_i = s_i + Cij[nPads * i + j] * qRatio[j];
188 }
189 if (Ci[i] > 1.0e-10) {
190 newQPixels[i] = s_i * qPixels[i] / Ci[i];
191 } else {
192 newQPixels[i] = 0;
193 }
194 }
195 } else {
196 for (int i = 0; i < nPixels; i++) {
197 // Normalization term
198 // r = np.sum( Cij[:,j]*qPad[0:nPads] / qPadPrediction[0:nPads] )
199 double s_i = 0;
200 // newQPixels[i] = 0.;
201
202 for (int j = 0; j < nPads; j++) {
203 s_i = s_i + Cij[nPads * i + j] * qRatio[j];
204 }
205 if (Ci[i] > 1.0e-10) {
206 newQPixels[i] = s_i * qPixels[i] / Ci[i];
207 } else {
208 newQPixels[i] = 0;
209 }
210
211 /*
212 double s_i = 0;
213 if (Ci[i] > 1.0e-10) {
214 for (int j = 0; j < nPads; j++) {
215 s_i += Cij[nPads * i + j] * qRatio[j];
216 }
217 newQPixels[i] = s_i * qPixels[i] / Ci[i];
218 } else {
219 newQPixels[i] = 0;
220 }
221 */
222 }
223 }
224
225 delete[] qRatio;
226}
227
228double computeChiSquare(const Pads& pads, const double* qPredictedPads,
229 int iStart, int iEnd)
230{
231 // Compute Chi2 on unsaturated pads
232 double chi2 = 0.0;
233 const double* q = pads.getCharges();
234 const Mask_t* sat = pads.getSaturates();
235 for (int i = iStart; i < iEnd; i++) {
236 double var = (1 - sat[i]) * (q[i] - qPredictedPads[i]);
237 chi2 += var * var;
238 }
239 return chi2;
240}
241
242std::pair<double, double> computeChiSquare(const Pads& pads, const double* qPredictedPads,
243 int N)
244{
245 // Compute Chi2 on unsaturated pads
246 double chi20 = 0.0;
247 double chi21 = 0.0;
248 const double* q = pads.getCharges();
249 const Mask_t* cath = pads.getCathodes();
250 const Mask_t* sat = pads.getSaturates();
251 for (int i = 0; i < N; i++) {
252 double var = (1 - sat[i]) * (q[i] - qPredictedPads[i]);
253 if (cath[i] == 0) {
254 chi20 += var * var;
255 } else {
256 chi21 += var * var;
257 }
258 }
259 return std::make_pair(chi20, chi21);
260}
261
262std::pair<double, double> PoissonEMLoop(const Pads& pads, Pads& pixels,
263 const double* Cij, Mask_t* maskCij,
264 int qCutMode, double minPadError,
265 int nItMax)
266{
267 // The array pixels return the last state
268 //
269 // Init.
270 //
271 int nPads = pads.getNbrOfPads();
272 int nPixels = pixels.getNbrOfPads();
273 //
274 const double* x = pixels.getX();
275 const double* y = pixels.getY();
276 const double* dx = pixels.getDX();
277 const double* dy = pixels.getDY();
278 double* qPixels = new double[pixels.getNbrOfPads()];
279 vectorCopy(pixels.getCharges(), pixels.getNbrOfPads(), qPixels);
280 //
281 const double* qPads = pads.getCharges();
282 //
283 // MaskCij: Used to disable Cij contribution (disable pixels)
284 vectorSetShort(maskCij, 1, nPads * nPixels);
285 // Ci, Cj: sum of Cij by row or by column
286 double Ci[nPixels];
287 double qPadPrediction[nPads];
288 double previousQPixels[nPixels];
289 double qPixCut;
290 // Init convergence criteria
291 bool converge = false;
292 int it = 0;
293 if (clusterConfig.EMLocalMaxLog > clusterConfig.info) {
294 printf("Poisson EM\n");
295 printf(
296 " it. <Pixels_residu> <Pad_residu> max(Pad_residu) "
297 "sum(Pad_residu)/sum(qPad)\n");
298 }
299 double meanPixelsResidu = 0.0;
300 double maxPixelsResidu = 0.0;
301 double maxRelResidu;
302 double pixelVariation;
303 double padRelError;
304 //
305
306 while (!converge) {
307 //
308 // Filter pixels
309 //
310 if (qCutMode == -1) {
311 // Percent of the min charge
312 // qPixCut = 1.02 * vectorMin(qPixels, nPixels);
313 qPixCut = 1.0e-14;
314 } else {
315 // qCutMode = 0
316 // No filtering
317 qPixCut = 0.0;
318 }
319
320 // Disable pixels
321 if (qPixCut > 0.0) {
322 for (int i = 0; i < (nPixels); i++) {
323 if (qPixels[i] < qPixCut) {
324 // old version with mask vectorSetShort(&maskCij[nPads * i], 0, nPads);
325 }
326 }
327 }
328 // Update Ci
329 for (int i = 0; i < nPixels; i++) {
330 Ci[i] = 0;
331 int start = nPads * i;
332 int end = start + nPads;
333 // for (int j = 0; j < nPads; j++) {
334 for (int l = start; l < end; l++) {
335 // Ci[i] += Cij[nPads * i + j] * maskCij[nPads * i + j];
336 // Ci[i] += Cij[nPads * i + j];
337 Ci[i] += Cij[l];
338 }
339 }
340 // Not needed
341 /*
342 double Cj[nPads];
343 for (int j = 0; j < nPads; j++) {
344 Cj[j] = 0;
345 for (int i = 0; i < nPixels; i++) {
346 // Cj[j] += Cij[nPads * i + j] * maskCij[nPads * i + j];
347 Cj[j] += Cij[nPads * i + j];
348 }
349 }
350 */
351 // Store previous qPixels state
352 vectorCopy(qPixels, nPixels, previousQPixels);
353 //
354 //
355 // Poisson EM Iterations
356 //
357
358 // iterateEMPoisson( Cij, Ci, maskCij, qPixels, qPads, qPadPrediction,
359 // nPixels, nPads, qPixels);
360
361 fastIterateEMPoisson(Cij, Ci, previousQPixels, qPads, qPadPrediction, nPixels,
362 nPads, qPixels);
363
364 // Measure of pixel variation to stop
365 // the iteration if required
366 double pixResidu[nPixels];
367 vectorAddVector(previousQPixels, -1.0, qPixels, nPixels, pixResidu);
368 vectorAbs(pixResidu, nPixels, pixResidu);
369 // Pixel variation
370 pixelVariation = vectorSum(pixResidu, nPixels) / vectorSum(qPixels, nPixels);
371 int iMaxResidu = vectorArgMax(pixResidu, nPixels);
372 maxRelResidu = pixResidu[iMaxResidu] / previousQPixels[iMaxResidu];
373
374 // Relative error on pad prediction
375 // Compute a stdError normalized with the pad Esperance
376 double padResidu[nPads];
377 vectorAddVector(qPads, -1.0, qPadPrediction, nPads, padResidu);
378 double var = vectorDotProd(padResidu, padResidu, nPads) / nPads;
379 double E = vectorSum(qPads, nPads) / nPads;
380 padRelError = std::sqrt(var) / E;
381
382 if (clusterConfig.EMLocalMaxLog > clusterConfig.info) {
383 printf(" EM it=%d <pixelResidu>=%10.6f, dQPixel/qPixel=%10.6f, max(dQPix)/qPix=%10.6f, relPadError=%10.6f\n",
384 it, vectorSum(pixResidu, nPixels) / nPixels, pixelVariation, maxRelResidu, padRelError);
385 }
386 // maxPixelVariation = 1 / 20 * minPadError;
387 converge = (pixelVariation < minPadError * 0.03) && (padRelError < minPadError) ||
388 (it > nItMax);
389 it += 1;
390 }
391 if (clusterConfig.EMLocalMaxLog > clusterConfig.info) {
392 printf(" Exit criterom pixelVariation=%d padRelError=%d itend=%d \n", (pixelVariation < minPadError * 0.03), (padRelError < minPadError), (it > nItMax));
393 }
394 // Update pixels charge
395 // Remove small charged pixels (<qPixCut)
396 int oldValueNPads = pixels.getNbrOfPads();
397 pixels.setCharges(qPixels, nPixels);
398 int k = 0;
399 if (qPixCut > 0.0) {
400 k = pixels.removePads(qPixCut);
401 }
402 // Chi2 on cathodes 0/1
403 double chi20, chi21;
404 // Chi2 on cathode1
405 std::pair<double, double> chi = computeChiSquare(pads, qPadPrediction, pads.getNbrOfPads());
406 std::pair<double, double> chiObs = computeChiSquare(pads, qPadPrediction, pads.getNbrOfObsPads());
407 if (clusterConfig.EMLocalMaxLog > clusterConfig.info) {
408 printf(" ??? Chi2 over NbrPads = (%f, %f); Chi2 over NbrObsPads = (%f, %f) \n",
409 chi.first, chi.second, chiObs.first, chiObs.second);
410 }
411 // ??? Must chose a method. A the moment over pads is better
412 if (1) {
413 // Chi2 over Pads
414 chi20 = chi.first;
415 chi21 = chi.second;
416 } else {
417 // Chi2 over ObsPads
418 chi20 = chiObs.first;
419 chi21 = chiObs.second;
420 }
421 // Take care to the leadind dimension is getNbrOfPads()
422 if (clusterConfig.EMLocalMaxLog >= clusterConfig.info) {
423 printf("End poisson EM :\n");
424 printf(" Total Pad Charge = %14.6f\n", vectorSum(qPads, nPads));
425 printf(" Total Predicted Charge = %14.6f\n",
426 vectorSum(qPadPrediction, nPads));
427 printf(" Total Pixel Charge = %14.6f\n", vectorSum(qPixels, nPixels));
428 printf(" EMPoisson qPixCut = %14.6f\n", qPixCut);
429 printf(" Chi20 = %14.6f\n", chi20);
430 printf(" Chi21 = %14.6f\n", chi21);
431 printf(" # of iterations = %d\n", it);
432 printf(" Max Pixel variation = %14.6f\n", maxPixelsResidu);
433 printf(" Nbr of removed pixels = %d/%d\n", oldValueNPads - k,
434 oldValueNPads);
435 }
436 delete[] qPixels;
437 return std::make_pair(chi20, chi21);
438}
439
440} // namespace mch
441} // namespace o2
ClusterConfig.h
Clustering and fifting parameters.
i
int32_t i
Definition GPUCommonAlgorithm.h:443
j
uint32_t j
Definition RawData.h:0
o2::mch::Pads::getCathodes
const Mask_t * getCathodes() const
Definition PadsPEM.h:101
o2::mch::Pads::getNbrOfPads
static int getNbrOfPads(const Pads *pads)
Definition PadsPEM.h:63
o2::mch::Pads::getCharges
const double * getCharges() const
Definition PadsPEM.h:99
o2::mch::Pads::getNbrOfObsPads
int getNbrOfObsPads() const
Definition PadsPEM.h:90
o2::mch::Pads::getSaturates
const Mask_t * getSaturates() const
Definition PadsPEM.h:100
x
GLint GLenum GLint x
Definition glcorearb.h:403
end
GLuint GLuint end
Definition glcorearb.h:469
pixels
GLint GLint GLsizei GLint GLenum GLenum const void * pixels
Definition glcorearb.h:275
start
GLuint start
Definition glcorearb.h:469
o2::mch::fastIterateEMPoissonV0
void fastIterateEMPoissonV0(const double *Cij, const double *Ci, const double *qPixels, const double *qPad, double *qPadPrediction, int nPixels, int nPads, double *newQPixels)
Definition poissonEM.cxx:84
o2::mch::computeChiSquare
double computeChiSquare(const Pads &pads, const double *qPredictedPads, int iStart, int iEnd)
Definition poissonEM.cxx:228
o2::mch::fastIterateEMPoisson
void fastIterateEMPoisson(const double *Cij, const double *Ci, const double *qPixels, const double *qPad, double *qPadPrediction, int nPixels, int nPads, double *newQPixels)
Definition poissonEM.cxx:137
o2::mch::iterateEMPoisson
void iterateEMPoisson(const double *Cij, const double *Ci, const Mask_t *maskCij, const double *qPixels, const double *qPad, double *qPadPrediction, int nPixels, int nPads, double *newQPixels)
Definition poissonEM.cxx:35
o2::mch::clusterConfig
ClusterConfig clusterConfig
Definition ClusterConfig.cxx:34
o2::mch::PoissonEMLoop
std::pair< double, double > PoissonEMLoop(const Pads &pads, Pads &pixels, const double *Cij, Mask_t *maskCij, int qCutMode, double minPadError, int nItMax)
Definition poissonEM.cxx:262
o2
a couple of static helper functions to create timestamp values for CCDB queries or override obsolete ...
Definition BitstreamReader.h:24
o2::mch::ClusterConfig::info
@ info
Describes main steps and high level behaviors.
Definition ClusterConfig.h:62