Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 14:14:24

0001 #ifndef MinL3AlgoUnivErr_H
0002 #define MinL3AlgoUnivErr_H
0003 
0004 //=============================================================================
0005 
0006 /** class MinL3AlgoUnivErr
0007 
0008  * \author R.Ofierzynski, CERN, 2007/08/23
0009  *                              under class name MinL3AlgoUniv
0010  *  Modified by A.Fedotov :
0011  *           24.07.09: a calculation of statistical errors implemented
0012  *                     on top of revision 1.2 2007/08/23 12:38:02 ;
0013  *                     the code remains backward compatible
0014  *           20.10.09: class name changed to MinL3AlgoUnivErr in order to
0015  *                     exclude any effect on older applications that make use
0016  *                     of the MinL3AlgoUniv class
0017  *
0018  *=============================================================================
0019  *
0020  *  General purpose
0021  *  ~~~~~~~~~~~~~~~
0022  *  Implementation of the L3 Collaboration algorithm to solve a system
0023  *                        Ax = B
0024  *  by minimization of |Ax-B| using an iterative linear approach
0025  *
0026  *  This class should be universal, i.e. working with DetIds or whatever else 
0027  *  will be invented to identify Subdetector parts
0028  *
0029  *  The bookkeeping of the cluster size and its elements
0030  *  has to be done by the user.
0031  *
0032  *  Calculation of statistical errors
0033  *  ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
0034  *  The solution whose errors have to be obtained, is found by a call to
0035  *  the `iterate' function. The errors are also found within that procedure
0036  *  in a general phenomenological approach consisting in
0037  *    - splitting the full sample to a certain number n of equal subsamples
0038  *                           (an optional argument nSubsamples of `iterate'
0039  *                            determines n; n = 10  by default)    
0040  *    - solving the problem for every part separately
0041  *    - then the spread of partial solutions is a measure of the error:
0042  *                    error = rms / sqrt (n - 1).                          (1)
0043  *  The relative precision of such estimate is believed to be 
0044  *          1 / sqrt[2(n - 1]            which yields 24% for n = 10.
0045  *  The user can fetch the errors by calling a `getError' function, and
0046  *  the average partial solution -- via `getMeanPartialSolution'.
0047  *
0048  *  Known PROBLEMS:
0049  *     1. If the event statistics for a particular cell is low enough, then
0050  *  the number of subsamples where the cell is present, n_cell, can be less
0051  *  than n (e.g, it is always so if a cell is active in 5 events while we
0052  *  split the full sample into n=10 parts). Then the error of this cell
0053  *  becomes wrong because a part of the full statistics gets lost for the cell
0054  *  (b.t.w., n_cell is actually used in eq.(1) ).
0055  *  The user can check the presence of such cells via function
0056  *  `numberOfWrongErrors' and check which cells have wrong errors with the
0057  *  functions `getErrorQuality' which gives the ratio n_cell/n -- the fraction
0058  *  of the full statistics used for the error estimation.
0059  *    2. Cases have been observed where the total solution converged nicely 
0060  *  with the increasing no. of iterations, while some of partial solutions 
0061  *  did not. Apparently, the errors can explode in such cases and do not
0062  *  reflect the real stat. errors of the total solution. This seems to be
0063  *  a problem of instabilities of the L3 method itself.
0064  *
0065  */
0066 
0067 //=============================================================================
0068 
0069 #include <vector>
0070 #include <iostream>
0071 #include <map>
0072 #include <cmath>
0073 
0074 //=============================================================================
0075 template <class IDdet>
0076 class MinL3AlgoUnivErr {
0077 public:
0078   typedef std::map<IDdet, float> IDmapF;
0079   typedef typename IDmapF::value_type IDmapFvalue;
0080   typedef typename IDmapF::iterator iter_IDmapF;
0081   typedef typename IDmapF::const_iterator citer_IDmapF;
0082   typedef std::map<IDdet, int> IDmapI;
0083   typedef typename IDmapI::value_type IDmapIvalue;
0084   typedef typename IDmapI::iterator iter_IDmapI;
0085   typedef typename IDmapI::const_iterator citer_IDmapI;
0086 
0087   //----------------------------------------------
0088   /// Default constructor
0089   /// kweight_ = event weight
0090 
0091   MinL3AlgoUnivErr(float kweight_ = 0.);
0092 
0093   //----------------------------------------------
0094   /// Destructor
0095 
0096   ~MinL3AlgoUnivErr();
0097 
0098   //----------------------------------------------
0099   ///  method doing the full calibration running nIter number of times,
0100   ///  recalibrating the event matrix after each iteration with the
0101   ///  new solution
0102   ///  returns the vector of calibration coefficients built
0103   ///  from all iteration solutions
0104   ///
0105   ///  The calibration is also done for nSubsamples sub-samples in order
0106   ///  to be able to estimate statistical errors of the main solution.
0107   ///
0108   ///     >> also to be used also as recipe on how to use the calibration
0109   ///     >> methods one-by-one with a re-selection of the events in between
0110   ///  the iterations<<
0111 
0112   IDmapF iterate(const std::vector<std::vector<float> >& eventMatrix,
0113                  const std::vector<std::vector<IDdet> >& idMatrix,
0114                  const std::vector<float>& energyVector,
0115                  const int& nIter,
0116                  const bool& normalizeFlag = false,
0117                  const int& nSubsamples = 10);
0118   //----------------------------------------------
0119   /// method to get the stat. error on the correction factor for cell id
0120   /// (to be called after the `iterate')
0121   ///
0122   ///    special values: getError = -2. : no info for the cell
0123   ///                             = -1. : the cell was met in one partial
0124   ///                                     solution only => the error equals to
0125   ///                                     INFINITY
0126 
0127   float getError(IDdet id) const;
0128 
0129   //----------------------------------------------
0130   /// method to get the stat. errors on the correction factors
0131   /// for all cells together (to be called after the `iterate').
0132   /// A map (id,error) is returned containing all the cells
0133   /// for which the information is available
0134   ///
0135   ///    special value: error = -1. : the cell was met in one partial
0136   ///                                 solution only => the error equals to
0137   ///                                 INFINITY
0138 
0139   IDmapF getError() const;
0140 
0141   //----------------------------------------------
0142   /// method to get the number of cells where the errors are incorrect
0143   /// due to nSamples_cell < nSubsamples
0144   /// (to be called after the `iterate') .
0145   ///
0146   /// this can happen for a cell if the number of events nev_cell
0147   /// where the cell is active (i.e. energy > 0),  is small enough,
0148   /// e.g. it is always so if nev_cell < nSubsamples.
0149 
0150   int numberOfWrongErrors() const;
0151 
0152   //----------------------------------------------
0153   /// two methods to get the fraction of full statistics used in calculation
0154   /// of stat. errors (to be called after the `iterate').
0155   ///
0156   /// a fraction < 1 indicates that the error is incorrect.
0157   ///
0158   ///     version with one argument: the fraction for a particular cells
0159   ///                                (special returned value: = -1. if no
0160   ///                                info for the cell is available =>
0161   ///                                a wrong id has been given)
0162   ///             w/o arguments    : the fractions for all cells together
0163 
0164   float getErrorQuality(IDdet id) const;
0165 
0166   IDmapF getErrorQuality() const;
0167 
0168   //----------------------------------------------
0169   /// method to get
0170   /// the mean partial solution for the correction factor for cell id
0171   /// (to be called after the `iterate')
0172   ///
0173   ///    special return value: 999. : no info for the cell
0174 
0175   float getMeanPartialSolution(IDdet id) const;
0176 
0177   //----------------------------------------------
0178   /// method to get the mean partial solution for all cells together
0179   /// (to be called after the `iterate')
0180   /// A map (id,mean) is returned containing all the cells
0181   /// for which the information is available
0182 
0183   IDmapF getMeanPartialSolution() const;
0184 
0185   //----------------------------------------------
0186   /// add event to the calculation of the calibration vector
0187 
0188   void addEvent(const std::vector<float>& myCluster, const std::vector<IDdet>& idCluster, const float& energy);
0189 
0190   //----------------------------------------------
0191   /// recalibrate before next iteration:
0192   /// give previous solution vector as argument
0193 
0194   std::vector<float> recalibrateEvent(const std::vector<float>& myCluster,
0195                                       const std::vector<IDdet>& idCluster,
0196                                       const IDmapF& newCalibration,
0197                                       const int& isol = -1,
0198                                       const int& iter = -1);
0199 
0200   //----------------------------------------------
0201   /// get the solution at the end of the calibration as a map between
0202   /// DetIds and calibration constant
0203 
0204   IDmapF getSolution(const bool resetsolution = true);
0205 
0206   //----------------------------------------------
0207   /// reset for new iteration
0208 
0209   void resetSolution();
0210 
0211   //----------------------------------------------
0212 
0213 private:
0214   float kweight;
0215   int countEvents;
0216   IDmapF wsum;
0217   IDmapF Ewsum;
0218   IDmapI idToIndex;                 // map: cell id -> index of cell info
0219                                     // in sumPartSolu... vectors
0220   std::vector<int> sumPartSolu0;    // number of partial solutions
0221   std::vector<float> sumPartSolu1;  // sum of partial solutions
0222   std::vector<float> sumPartSolu2;  // sum of squared partial solutions
0223 
0224   int nOfSubsamples;  //  a copy of an input argument of `iterate' function
0225 
0226   //----------------------------------------------
0227   // register a partial solution in data members
0228   void registerPartialSolution(const IDmapF& partialSolution);
0229 
0230   //----------------------------------------------
0231 
0232 };  //end of class MinL3AlgoUnivErr prototype
0233 
0234 //=============================================================================
0235 
0236 template <class IDdet>
0237 MinL3AlgoUnivErr<IDdet>::MinL3AlgoUnivErr(float kweight_) : kweight(kweight_), countEvents(0) {
0238   std::cout << "MinL3AlgoUnivErr : L3 algo with a calculation"
0239             << " of stat.errors working..." << std::endl;
0240   resetSolution();
0241 }
0242 
0243 //=============================================================================
0244 
0245 template <class IDdet>
0246 MinL3AlgoUnivErr<IDdet>::~MinL3AlgoUnivErr() {}
0247 
0248 //=============================================================================
0249 
0250 template <class IDdet>
0251 typename MinL3AlgoUnivErr<IDdet>::IDmapF MinL3AlgoUnivErr<IDdet>::iterate(
0252     const std::vector<std::vector<float> >& eventMatrix,
0253     const std::vector<std::vector<IDdet> >& idMatrix,
0254     const std::vector<float>& energyVector,
0255     const int& nIter,
0256     const bool& normalizeFlag,
0257     const int& nSubsamples) {
0258   // clear the data members which are filled inside the function
0259   //                               (in registerPartialSolution called from here)
0260 
0261   nOfSubsamples = nSubsamples;  // keep the input argument for use in other
0262                                 // functions
0263 
0264   idToIndex.clear();
0265   sumPartSolu0.clear();
0266   sumPartSolu1.clear();
0267   sumPartSolu2.clear();
0268 
0269   IDmapF totalSolution;
0270 
0271   // Loop over samples/solutions:
0272   //    isol = 0                 : all events with the solution stored in
0273   //                               totalSolution
0274   //    isol = 1,...,nSubsamples : partial solutions are found for sub-samples
0275   //                               with the info on the solutions stored in the
0276   //                               data members
0277   //                                    idToIndex, sumPartSolu...
0278   //                               in order to be able to estimate the stat.
0279   //                               errors later
0280 
0281   for (int isol = 0; isol <= nSubsamples; isol++) {
0282     IDmapF sampleSolution;  // solution for the sample
0283     IDmapF iterSolution;    // intermediate solution after an iteration
0284     std::vector<std::vector<float> > myEventMatrix;
0285     std::vector<std::vector<IDdet> > myIdMatrix;
0286     std::vector<float> myEnergyVector;
0287 
0288     // Select the sample.
0289     // Fill myEventMatrix, myIdMatrix and myEnergyVector
0290     // either with all evs or with independent event subsamples
0291 
0292     if (isol == 0)  // total solution
0293     {
0294       myEventMatrix = eventMatrix;
0295       myIdMatrix = idMatrix;
0296       myEnergyVector = energyVector;
0297     } else  // partial solution # isol
0298     {
0299       // clear containers filled for the previous sample
0300       sampleSolution.clear();
0301       myEventMatrix.clear();
0302       myIdMatrix.clear();
0303       myEnergyVector.clear();
0304 
0305       for (int i = 0; i < static_cast<int>(eventMatrix.size()); i++) {
0306         // select every nSubsamples'th event to the subsample
0307         if (i % nSubsamples + 1 == isol) {
0308           myEventMatrix.push_back(eventMatrix[i]);
0309           myIdMatrix.push_back(idMatrix[i]);
0310           myEnergyVector.push_back(energyVector[i]);
0311         }
0312       }
0313     }
0314 
0315     int Nevents = myEventMatrix.size();  // Number of events to calibrate with
0316     countEvents = 0;
0317 
0318     int i;
0319 
0320     // Iterate the correction
0321     for (int iter = 1; iter <= nIter; iter++) {
0322       // if normalization flag is set, normalize energies
0323       float sumOverEnergy;
0324       if (normalizeFlag) {
0325         float scale = 0.;
0326 
0327         for (i = 0; i < Nevents; i++) {
0328           sumOverEnergy = 0.;
0329           for (unsigned j = 0; j < myEventMatrix[i].size(); j++) {
0330             sumOverEnergy += myEventMatrix[i][j];
0331           }
0332           sumOverEnergy /= myEnergyVector[i];
0333           scale += sumOverEnergy;
0334         }
0335         scale /= Nevents;
0336 
0337         for (i = 0; i < Nevents; i++) {
0338           myEnergyVector[i] *= scale;
0339         }
0340       }  // end normalize energies
0341 
0342       // now the real work starts:
0343       for (int iEvt = 0; iEvt < Nevents; iEvt++) {
0344         addEvent(myEventMatrix[iEvt], myIdMatrix[iEvt], myEnergyVector[iEvt]);
0345       }
0346       iterSolution = getSolution();
0347       if (iterSolution.empty()) {
0348         sampleSolution.clear();
0349         break;  // exit the iteration loop leaving sampleSolution empty
0350       }
0351 
0352       // re-calibrate eventMatrix with solution
0353       for (int ievent = 0; ievent < Nevents; ievent++) {
0354         myEventMatrix[ievent] = recalibrateEvent(myEventMatrix[ievent], myIdMatrix[ievent], iterSolution, isol, iter);
0355       }
0356 
0357       // save solution into the sampleSolution map
0358       for (iter_IDmapF i = iterSolution.begin(); i != iterSolution.end(); i++) {
0359         iter_IDmapF itotal = sampleSolution.find(i->first);
0360         if (itotal == sampleSolution.end()) {
0361           sampleSolution.insert(IDmapFvalue(i->first, i->second));
0362         } else {
0363           itotal->second *= i->second;
0364         }
0365       }
0366 
0367       //      resetSolution(); // reset for new iteration,
0368       //               now: getSolution does it automatically if not vetoed
0369     }  // end iterate correction
0370 
0371     if (isol == 0)  // total solution
0372     {
0373       totalSolution = sampleSolution;
0374     } else  // partial solution => register it in sumPartSolu...
0375     {
0376       registerPartialSolution(sampleSolution);
0377     }
0378 
0379   }  // end of the loop over solutions/samples
0380 
0381   return totalSolution;
0382 }
0383 
0384 //=============================================================================
0385 
0386 template <class IDdet>
0387 void MinL3AlgoUnivErr<IDdet>::addEvent(const std::vector<float>& myCluster,
0388                                        const std::vector<IDdet>& idCluster,
0389                                        const float& energy) {
0390   countEvents++;
0391 
0392   float w, invsumXmatrix;
0393   float eventw;
0394 
0395   // Loop over the crystal matrix to find the sum
0396   float sumXmatrix = 0.;
0397   for (unsigned i = 0; i < myCluster.size(); i++) {
0398     sumXmatrix += myCluster[i];
0399   }
0400 
0401   // event weighting
0402   eventw = 1 - fabs(1 - sumXmatrix / energy);
0403   eventw = pow(eventw, kweight);
0404 
0405   if (sumXmatrix != 0.) {
0406     invsumXmatrix = 1 / sumXmatrix;
0407     // Loop over the crystal matrix (3x3,5x5,7x7) again
0408     // and calculate the weights for each xtal
0409     for (unsigned i = 0; i < myCluster.size(); i++) {
0410       w = myCluster[i] * invsumXmatrix;
0411 
0412       // include the weights into wsum, Ewsum
0413       iter_IDmapF iwsum = wsum.find(idCluster[i]);
0414       if (iwsum == wsum.end())
0415         wsum.insert(IDmapFvalue(idCluster[i], w * eventw));
0416       else
0417         iwsum->second += w * eventw;
0418 
0419       iter_IDmapF iEwsum = Ewsum.find(idCluster[i]);
0420       if (iEwsum == Ewsum.end())
0421         Ewsum.insert(IDmapFvalue(idCluster[i], (w * eventw * energy * invsumXmatrix)));
0422       else
0423         iEwsum->second += (w * eventw * energy * invsumXmatrix);
0424     }
0425   }
0426   //  else {std::cout << " Debug: dropping null event: " << countEvents << std::endl;}
0427 }
0428 
0429 //=============================================================================
0430 
0431 template <class IDdet>
0432 typename MinL3AlgoUnivErr<IDdet>::IDmapF MinL3AlgoUnivErr<IDdet>::getSolution(const bool resetsolution) {
0433   IDmapF solution;
0434 
0435   for (iter_IDmapF i = wsum.begin(); i != wsum.end(); i++) {
0436     iter_IDmapF iEwsum = Ewsum.find(i->first);
0437     float myValue = 1;
0438     if (i->second != 0)
0439       myValue = iEwsum->second / i->second;
0440 
0441     solution.insert(IDmapFvalue(i->first, myValue));
0442   }
0443 
0444   if (resetsolution)
0445     resetSolution();
0446 
0447   return solution;
0448 }
0449 
0450 //=============================================================================
0451 
0452 template <class IDdet>
0453 void MinL3AlgoUnivErr<IDdet>::resetSolution() {
0454   wsum.clear();
0455   Ewsum.clear();
0456 }
0457 
0458 //=============================================================================
0459 
0460 template <class IDdet>
0461 std::vector<float> MinL3AlgoUnivErr<IDdet>::recalibrateEvent(const std::vector<float>& myCluster,
0462                                                              const std::vector<IDdet>& idCluster,
0463                                                              const IDmapF& newCalibration,
0464                                                              const int& isol,  // for a printout only
0465                                                              const int& iter   // for a printout only
0466 ) {
0467   std::vector<float> newCluster(myCluster);
0468 
0469   for (unsigned i = 0; i < myCluster.size(); i++) {
0470     citer_IDmapF icalib = newCalibration.find(idCluster[i]);
0471     if (icalib != newCalibration.end()) {
0472       newCluster[i] *= icalib->second;
0473     } else {
0474       std::cout << "No calibration available for this element." << std::endl;
0475       std::cout << "   isol = " << isol << "   iter = " << iter << "   idCluster[i] = " << idCluster[i] << "\n";
0476     }
0477   }
0478 
0479   return newCluster;
0480 }
0481 
0482 //=============================================================================
0483 
0484 template <class IDdet>
0485 void MinL3AlgoUnivErr<IDdet>::registerPartialSolution(const IDmapF& partialSolution)
0486 
0487 {
0488   int lastIndex = sumPartSolu0.size() - 1;  // index of the last element
0489                                             // of the parallel vectors
0490 
0491   for (citer_IDmapF cell = partialSolution.begin(); cell != partialSolution.end(); ++cell) {
0492     IDdet id = cell->first;
0493     float corr = cell->second;
0494 
0495     // where is the cell in another map?
0496     iter_IDmapI cell2 = idToIndex.find(id);
0497 
0498     if (cell2 == idToIndex.end()) {
0499       // the cell is met for the first time in patial solutions
0500       // => insert the info to the end of the vectors
0501 
0502       sumPartSolu0.push_back(1);
0503       sumPartSolu1.push_back(corr);
0504       sumPartSolu2.push_back(corr * corr);
0505       idToIndex.insert(IDmapIvalue(id, ++lastIndex));
0506     } else {
0507       // add the info to the already registered cell
0508 
0509       int index = cell2->second;
0510       sumPartSolu0[index] += 1;
0511       sumPartSolu1[index] += corr;
0512       sumPartSolu2[index] += corr * corr;
0513     }
0514   }
0515 }
0516 
0517 //=============================================================================
0518 
0519 template <class IDdet>
0520 float MinL3AlgoUnivErr<IDdet>::getError(IDdet id) const
0521 
0522 {
0523   float error;
0524   citer_IDmapI cell = idToIndex.find(id);
0525   if (cell == idToIndex.end())
0526     error = -2.;  // no info for the cell
0527   else {
0528     int i = cell->second;
0529     int n = sumPartSolu0[i];
0530     if (n <= 1)
0531       error = -1.;  // 1 entry => error estimate impossible
0532     else {
0533       float meanX = sumPartSolu1[i] / n;
0534       float meanX2 = sumPartSolu2[i] / n;
0535 
0536       error = sqrt(fabs(meanX2 - meanX * meanX) / (n - 1.));
0537     }
0538   }
0539   return error;
0540 }
0541 
0542 //=============================================================================
0543 
0544 template <class IDdet>
0545 typename MinL3AlgoUnivErr<IDdet>::IDmapF MinL3AlgoUnivErr<IDdet>::getError() const
0546 
0547 {
0548   IDmapF errors;
0549   float error;
0550 
0551   for (citer_IDmapI cell = idToIndex.begin(); cell != idToIndex.end(); ++cell) {
0552     int i = cell->second;
0553     int n = sumPartSolu0[i];
0554     if (n <= 1)
0555       error = -1.;  // 1 entry => error estimate impossible
0556     else {
0557       float meanX = sumPartSolu1[i] / n;
0558       float meanX2 = sumPartSolu2[i] / n;
0559 
0560       error = sqrt(fabs(meanX2 - meanX * meanX) / (n - 1));
0561     }
0562 
0563     errors.insert(IDmapFvalue(cell->first, error));
0564   }
0565   return errors;
0566 }
0567 //=============================================================================
0568 
0569 template <class IDdet>
0570 float MinL3AlgoUnivErr<IDdet>::getErrorQuality(IDdet id) const
0571 
0572 {
0573   float fraction;
0574   citer_IDmapI cell = idToIndex.find(id);
0575   if (cell == idToIndex.end())
0576     fraction = -1.;  // no info for the cell
0577                      // => return a special value
0578   else {
0579     int i = cell->second;
0580     int n = sumPartSolu0[i];
0581     if (n < nOfSubsamples)
0582       fraction = float(n) / nOfSubsamples;
0583     else
0584       fraction = 1.;
0585   }
0586   return fraction;
0587 }
0588 
0589 //=============================================================================
0590 
0591 template <class IDdet>
0592 typename MinL3AlgoUnivErr<IDdet>::IDmapF MinL3AlgoUnivErr<IDdet>::getErrorQuality() const
0593 
0594 {
0595   IDmapF fractions;
0596   float fraction;
0597 
0598   for (citer_IDmapI cell = idToIndex.begin(); cell != idToIndex.end(); ++cell) {
0599     int i = cell->second;
0600     int n = sumPartSolu0[i];
0601 
0602     if (n < nOfSubsamples)
0603       fraction = float(n) / nOfSubsamples;
0604     else
0605       fraction = 1.;
0606 
0607     fractions.insert(IDmapFvalue(cell->first, fraction));
0608   }
0609 
0610   return fractions;
0611 }
0612 
0613 //=============================================================================
0614 
0615 template <class IDdet>
0616 int MinL3AlgoUnivErr<IDdet>::numberOfWrongErrors() const
0617 
0618 {
0619   int nWrong = 0;
0620 
0621   for (citer_IDmapI cell = idToIndex.begin(); cell != idToIndex.end(); ++cell) {
0622     int i = cell->second;
0623     int n = sumPartSolu0[i];
0624 
0625     if (n < nOfSubsamples)
0626       nWrong++;
0627   }
0628 
0629   return nWrong;
0630 }
0631 
0632 //=============================================================================
0633 
0634 template <class IDdet>
0635 float MinL3AlgoUnivErr<IDdet>::getMeanPartialSolution(IDdet id) const
0636 
0637 {
0638   float meanX;
0639   citer_IDmapI cell = idToIndex.find(id);
0640   if (cell == idToIndex.end())
0641     meanX = 999.;  // no info for the cell
0642   else {
0643     int i = cell->second;
0644     int n = sumPartSolu0[i];
0645     meanX = sumPartSolu1[i] / n;
0646   }
0647   return meanX;
0648 }
0649 
0650 //=============================================================================
0651 
0652 template <class IDdet>
0653 typename MinL3AlgoUnivErr<IDdet>::IDmapF MinL3AlgoUnivErr<IDdet>::getMeanPartialSolution() const
0654 
0655 {
0656   IDmapF solution;
0657 
0658   for (citer_IDmapI cell = idToIndex.begin(); cell != idToIndex.end(); ++cell) {
0659     int i = cell->second;
0660     int n = sumPartSolu0[i];
0661     float meanX = sumPartSolu1[i] / n;
0662     solution.insert(IDmapFvalue(cell->first, meanX));
0663   }
0664   return solution;
0665 }
0666 
0667 //=============================================================================
0668 
0669 #endif  // MinL3AlgoUnivErr_H