Merge branch 'MultiComponentStateModeCalculator_MT' into 'master'

ATLASRECTS-4612 MultiComponentStateModeCalculator : Remove mutable state, counters atomic, pass the array of vectors of mixture as arguement

See merge request atlas/athena!14567

Tracking/TrkFitter/TrkGaussianSumFilter/TrkGaussianSumFilter/MultiComponentStateModeCalculator.h

 /*
-  Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+  Copyright (C) 2002-2018 CERN for the benefit of the ATLAS collaboration
 */
 
 /***********************************************************************************
@@ -8,7 +8,7 @@
 begin                : Thursday 6th July 2006
 author               : atkinson
 email                : Tom.Atkinson@cern.ch
-description          : Class to calculate the mode (q/p) of a gaussian mixture
+description          : Class to calculate the mode (q/p) of a gaussian mixtureArray
 ***********************************************************************************/
 
 #ifndef Trk_MultiComponentStateModeCalculator_H
@@ -16,6 +16,9 @@ description          : Class to calculate the mode (q/p) of a gaussian mixture
 
 #include "AthenaBaseComps/AthAlgTool.h"
 #include "TrkGaussianSumFilter/IMultiComponentStateModeCalculator.h"
+#include <atomic>
+#include <array>
+#include <vector>
 
 namespace Trk{
 
@@ -38,72 +41,75 @@ class MultiComponentStateModeCalculator : public AthAlgTool, virtual public IMul
   StatusCode finalize();
 
   //!< IMultiComponentStateModeCalculator interface method to calculate mode
-  Amg::VectorX calculateMode( const MultiComponentState& ) const;
-
- private:
-  
-  //!< Private method to extract the weight, mean and sigma values from the multi-component state
-  void fillMixture( const MultiComponentState& ) const;
-
-  //!< Private method to find the mode using the Newton-Raphson method based on a starting guess
-  double findMode( double, int ) const;
-
-  //!< Private method to determine the pdf of the cashed mixture at a given value
-  double pdf( double, int ) const;
-
-  //!< Private method to determine the first order derivative of the pdf at a given value
-  double d1pdf( double, int ) const;
-
-  //!< Private method to determine the second order derivative of the pdf at a given value
-  double d2pdf( double, int) const;
-
-  //!< Private method to determine the value of the a gaussian distribution at a given value
-  double gaus( double x, double mean, double sigma ) const;
-
-  double findModeGlobal(double, int) const;
-  
-  double width( int i)  const;
-  
-  double findRoot(double &result, double xlo, double xhi, double value, double i) const;
-
+  virtual Amg::VectorX calculateMode( const MultiComponentState& ) const override;
 
  private:
 
-  //!< Private mixture structure
+  //!< Private Mixture structure
   struct Mixture{
 
     // Default constructor
-    Mixture()
-      :
+    Mixture():
       weight(0.),
       mean(0.),
       sigma(0.)
     {}
 
     // Constructor with arguments
-    Mixture( double aWeight, double aMean, double aSigma )
-      :
+    Mixture( double aWeight, double aMean, double aSigma ):
       weight( aWeight ),
       mean( aMean ),
       sigma( aSigma )
     {}
-
     double weight;
     double mean;
     double sigma;
   };
 
+ 
+  //!< Private method to extract the weight, mean and sigma values from the multi-component state
+  void fillMixture( std::array<std::vector< Mixture >,5>& mixtureArray, 
+                    const MultiComponentState& ) const;
+
+  //!< Private method to find the mode using the Newton-Raphson method based on a starting guess
+  double findMode( const std::array<std::vector< Mixture >,5>& mixtureArray,
+                   double, int ) const;
+
+  //!< Private method to determine the pdf of the cashed mixtureArray at a given value
+  double pdf( const std::array<std::vector< Mixture >,5>& mixtureArray,
+              double, int ) const;
+
+  //!< Private method to determine the first order derivative of the pdf at a given value
+  double d1pdf( const std::array<std::vector< Mixture >,5>& mixtureArray,
+                double, int ) const;
+
+  //!< Private method to determine the second order derivative of the pdf at a given value
+  double d2pdf( const std::array<std::vector< Mixture >,5>& mixtureArray,
+                double, int) const;
+
+  //!< Private method to determine the value of the a gaussian distribution at a given value
+  double gaus( double x, double mean, double sigma ) const;
+
+  double findModeGlobal(const std::array<std::vector< Mixture >,5>& mixtureArray,
+                        double, int) const;
+  
+  double width( const std::array<std::vector< Mixture >,5>& mixtureArray,
+                int i)  const;
+  
+  double findRoot(const std::array<std::vector< Mixture >,5>& mixtureArray,
+                  double &result, double xlo, double xhi, double value, double i) const;
+
+
  private:
 
   int                              m_outputlevel;                      //!< to cache current output level
   
-  mutable std::vector< Mixture >   m_mixture[5];    //!< Mixture
   
   //ModeCalualtor Stats
-  mutable int                      m_NumberOfCalls;
-  mutable int                      m_ConverganceFilures;
-  mutable int                      m_NoErrorMatrix;
-  mutable int                      m_MixtureSizeZero;
+  mutable std::atomic<int>                      m_NumberOfCalls;
+  mutable std::atomic<int>                      m_ConverganceFilures;
+  mutable std::atomic<int>                      m_NoErrorMatrix;
+  mutable std::atomic<int>                      m_MixtureSizeZero;
   
   
 };

Tracking/TrkFitter/TrkGaussianSumFilter/src/MultiComponentStateModeCalculator.cxx

 /*
-  Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+  Copyright (C) 2002-2018 CERN for the benefit of the ATLAS collaboration
 */
 
 /****************************************************************************************
@@ -17,8 +17,6 @@ description          : Implementation code for MultiComponentStateModeCalculator
 //#include "EventPrimitives/AmgMatrixPlugin.h"
 
 #include "TrkParameters/TrackParameters.h"
-
-
 #include <map>
 
 Trk::MultiComponentStateModeCalculator::MultiComponentStateModeCalculator( const std::string& type, const std::string& name, const IInterface* parent )
@@ -70,7 +68,7 @@ StatusCode Trk::MultiComponentStateModeCalculator::finalize()
 
 }
 
-Amg::VectorX Trk::MultiComponentStateModeCalculator::calculateMode( const Trk::MultiComponentState& multiComponentState ) const
+Amg::VectorX Trk::MultiComponentStateModeCalculator::calculateMode(const Trk::MultiComponentState& multiComponentState ) const
 {
   ++m_NumberOfCalls;
   
@@ -83,63 +81,65 @@ Amg::VectorX Trk::MultiComponentStateModeCalculator::calculateMode( const Trk::M
     ++m_NoErrorMatrix;
     return modes;
   }
-
-  fillMixture( multiComponentState );
-  
+  std::array<std::vector< Mixture >,5> mixtureArray{};
+  fillMixture( mixtureArray,multiComponentState );
+ 
+  int tmp_MixtureSizeZero{0}; 
   for (int i=0 ; i<5 ; i++){
     std::multimap< double, double, std::greater<double> > para_startMap;
-
-    // Loop over the mixture and find the starting point
-    std::vector<Mixture>::const_iterator component = m_mixture[i].begin();
+    // Loop over the mixtureArray and find the starting point
+    std::vector<Mixture>::const_iterator component = mixtureArray[i].begin();
  
-    if (m_mixture[i].size() == 0) ++m_MixtureSizeZero;
+    if (mixtureArray[i].size() == 0) {
+      ++tmp_MixtureSizeZero;
+    }
  
-    for( ; component != m_mixture[i].end() ; ++component )
-      para_startMap.insert( std::pair<double, double>( pdf( component->mean, i), component->mean ) );
-    
+    for( ; component != mixtureArray[i].end() ; ++component ){
+      para_startMap.insert( std::pair<double, double>( pdf( mixtureArray,component->mean, i), component->mean ) );
+    }
+  
     double para_start = para_startMap.begin()->second;
-
-    modes[i] =  findMode( para_start, i ) ;
+    modes[i] =  findMode( mixtureArray,para_start, i ) ;
 
     // Calculate the FWHM and return this back so that it can be used to correct the covariance matrix 
     if( para_start != modes[i] ){
       // mode calculation was successful now calulate FWHM
-      double currentWidth =  width(i);
+      double currentWidth =  width(mixtureArray,i);
       modes[i+5] = -1; // Failure is flagged with a value less than 0;
       
-      double pdfVal = pdf( modes[i], i );      
+      double pdfVal = pdf( mixtureArray,modes[i], i );      
       double highX(0); 
       double lowX(0);      
 
       double upperbound =modes[i] + 1.5 * currentWidth;
       while(true){
-        if( pdf( upperbound, i ) > pdfVal*0.5 ){
+        if( pdf( mixtureArray,upperbound, i ) > pdfVal*0.5 ){
           upperbound += currentWidth;
         } else {
           break;
         }
       }
 
-      ATH_MSG_VERBOSE ( "HighX  PDFval, high val, low val [ " << pdfVal << ", " << pdf( modes[i], i ) << ", " << pdf( upperbound, i ) << "]") ;    
+      ATH_MSG_VERBOSE ( "HighX  PDFval, high val, low val [ " << pdfVal << ", " << pdf( mixtureArray,modes[i], i ) << ", " << pdf( mixtureArray,upperbound, i ) << "]") ;    
 
-      bool highXFound = findRoot( highX, modes[i], upperbound, pdfVal*0.5, i);
+      bool highXFound = findRoot( mixtureArray,highX, modes[i], upperbound, pdfVal*0.5, i);
       
       double lowerbound =modes[i] - 1.5 * currentWidth;
       while(true){
-        if( pdf( lowerbound, i ) > pdfVal*0.5 ){
+        if( pdf( mixtureArray,lowerbound, i ) > pdfVal*0.5 ){
           lowerbound -= currentWidth;
         } else {
           break;
         }
       }
 
-      ATH_MSG_VERBOSE ( "LowX   PDFval, high val, low val [ " << pdfVal << ", " << pdf( lowerbound, i ) << ", " << pdf( modes[i], i ) << "]") ;    
+      ATH_MSG_VERBOSE ( "LowX   PDFval, high val, low val [ " << pdfVal << ", " << pdf( mixtureArray,lowerbound, i ) << ", " << pdf( mixtureArray,modes[i], i ) << "]") ;    
 
-      bool lowXFound  = findRoot( lowX , lowerbound, modes[i], pdfVal*0.5, i);
+      bool lowXFound  = findRoot( mixtureArray,lowX , lowerbound, modes[i], pdfVal*0.5, i);
       
       if (highXFound && lowXFound ){
         double FWHM = highX - lowX;
-        ATH_MSG_DEBUG ( "PDFval, high val, low val [ " << pdfVal << ", " << pdf( highX, i ) << ", " << pdf( lowX, i ) << "]") ;    
+        ATH_MSG_DEBUG ( "PDFval, high val, low val [ " << pdfVal << ", " << pdf( mixtureArray,highX, i ) << ", " << pdf( mixtureArray,lowX, i ) << "]") ;    
         
         if( FWHM <= 0 ) {
           ATH_MSG_DEBUG(i << " Width is neagtive? " << highX << " " << lowX << " " <<  modes[i] );
@@ -161,19 +161,19 @@ Amg::VectorX Trk::MultiComponentStateModeCalculator::calculateMode( const Trk::M
         }
       }
     }
-    
-
     para_startMap.clear();
   }
+  m_MixtureSizeZero+=tmp_MixtureSizeZero;
   return modes;
 
 }
 
-void Trk::MultiComponentStateModeCalculator::fillMixture( const Trk::MultiComponentState& multiComponentState ) const
+void Trk::MultiComponentStateModeCalculator::fillMixture( std::array<std::vector< Mixture >,5>& mixtureArray,
+                                                          const Trk::MultiComponentState& multiComponentState ) const
 {
   
   for (int i=0; i<5; i++){
-    m_mixture[i].clear();
+    mixtureArray[i].clear();
   }
   
   // Loop over all the components in the multi-component state
@@ -205,11 +205,8 @@ void Trk::MultiComponentStateModeCalculator::fillMixture( const Trk::MultiCompon
       }
     }
 
-
-
-    Mixture mixture(weight, mean, sigma );
-
-    m_mixture[i].push_back( mixture );
+    Mixture newmixture(weight, mean, sigma );
+    mixtureArray[i].push_back( newmixture);
    }
   }
 
@@ -217,7 +214,8 @@ void Trk::MultiComponentStateModeCalculator::fillMixture( const Trk::MultiCompon
 
 }
 
-double Trk::MultiComponentStateModeCalculator::findMode( double xStart, int i ) const
+double Trk::MultiComponentStateModeCalculator::findMode( const std::array<std::vector< Mixture >,5>& mixtureArray,
+                                                         double xStart, int i ) const
 {
 
   int iteration(0);
@@ -230,16 +228,17 @@ double Trk::MultiComponentStateModeCalculator::findMode( double xStart, int i )
 
       double previousMode(mode);
       
-      if (d2pdf( previousMode, i) != 0.0) {
-	      mode = previousMode - d1pdf( previousMode, i ) / d2pdf( previousMode, i );
+      if (d2pdf( mixtureArray,previousMode, i) != 0.0) {
+	      mode = previousMode - d1pdf(mixtureArray,previousMode, i ) / d2pdf( mixtureArray,previousMode, i );
       }
       else {
         ATH_MSG_DEBUG( "Second derivative is zero ... Returning the original value" );
         return xStart;
       }
 
-      if ( ( pdf(mode, i) + pdf(previousMode, i)) != 0.0 ) {
-        tolerance = fabs( pdf(mode, i) - pdf(previousMode, i) ) / ( pdf(mode, i) + pdf(previousMode, i) );
+      if ( ( pdf(mixtureArray,mode, i) + pdf(mixtureArray,previousMode, i)) != 0.0 ) {
+        tolerance = fabs( pdf(mixtureArray,mode, i) - pdf(mixtureArray,previousMode, i) ) 
+          / ( pdf(mixtureArray,mode, i) + pdf(mixtureArray,previousMode, i) );
       }
       else {
         ATH_MSG_DEBUG( "Dividing by zero ... Returning the original value" );
@@ -260,7 +259,8 @@ double Trk::MultiComponentStateModeCalculator::findMode( double xStart, int i )
 
 }
 
-double Trk::MultiComponentStateModeCalculator::findModeGlobal(double mean,int i) const
+double Trk::MultiComponentStateModeCalculator::findModeGlobal( const std::array<std::vector< Mixture >,5>& mixtureArray,
+                                                               double mean,int i) const
 {
   
   double start(-1);
@@ -279,7 +279,7 @@ double Trk::MultiComponentStateModeCalculator::findModeGlobal(double mean,int i)
   double iterate(fabs(mean/1000));
  
   for (double counter(start); counter < end; counter+=iterate) {
-    double value(pdf(counter,i));
+    double value(pdf(mixtureArray,counter,i));
     if (value > maximum) {
       maximum = value;
       mode = counter;
@@ -288,28 +288,30 @@ double Trk::MultiComponentStateModeCalculator::findModeGlobal(double mean,int i)
   return mode;
 }
 
-double Trk::MultiComponentStateModeCalculator::pdf( double x, int i ) const
+double Trk::MultiComponentStateModeCalculator::pdf( const std::array<std::vector< Mixture >,5>& mixtureArray, 
+                                                    double x, int i ) const
 {
 
   double pdf(0.);
 
-  std::vector< Mixture >::const_iterator component = m_mixture[i].begin();
+  std::vector< Mixture >::const_iterator component = mixtureArray[i].begin();
 
-  for ( ; component != m_mixture[i].end() ; ++component )
+  for ( ; component != mixtureArray[i].end() ; ++component )
     pdf += component->weight * gaus( x, component->mean, component->sigma );
 
   return pdf;
 
 }
 
-double Trk::MultiComponentStateModeCalculator::d1pdf( double x, int i  ) const
+double Trk::MultiComponentStateModeCalculator::d1pdf( const std::array<std::vector< Mixture >,5>& mixtureArray,
+                                                      double x, int i  ) const
 {
 
   double result(0.);
 
-  std::vector< Mixture >::const_iterator component = m_mixture[i].begin();
+  std::vector< Mixture >::const_iterator component = mixtureArray[i].begin();
 
-  for( ; component != m_mixture[i].end() ; ++component ){
+  for( ; component != mixtureArray[i].end() ; ++component ){
 
     double z = ( x - component->mean ) / component->sigma;
 
@@ -321,15 +323,16 @@ double Trk::MultiComponentStateModeCalculator::d1pdf( double x, int i  ) const
 
 }
 
-double Trk::MultiComponentStateModeCalculator::d2pdf( double x, int i ) const
+double Trk::MultiComponentStateModeCalculator::d2pdf( const std::array<std::vector< Mixture >,5>& mixtureArray,
+                                                      double x, int i ) const
 {
 
 
   double result(0.);
 
-  std::vector< Mixture >::const_iterator component = m_mixture[i].begin();
+  std::vector< Mixture >::const_iterator component = mixtureArray[i].begin();
 
-  for( ; component != m_mixture[i].end() ; ++component ){
+  for( ; component != mixtureArray[i].end() ; ++component ){
 
     double z = ( x - component->mean ) / component->sigma;
 
@@ -356,14 +359,15 @@ double Trk::MultiComponentStateModeCalculator::gaus( double x, double mean, doub
 }
 
 
-double Trk::MultiComponentStateModeCalculator::width( int i ) const
+double Trk::MultiComponentStateModeCalculator::width( const std::array<std::vector< Mixture >,5>& mixtureArray ,
+                                                      int i ) const
 {
 
   double pdf(0.);
 
-  std::vector< Mixture >::const_iterator component = m_mixture[i].begin();
+  std::vector< Mixture >::const_iterator component = mixtureArray[i].begin();
 
-  for ( ; component != m_mixture[i].end() ; ++component )
+  for ( ; component != mixtureArray[i].end() ; ++component )
     pdf += component->weight * component->sigma;
 
   return pdf;
@@ -372,7 +376,10 @@ double Trk::MultiComponentStateModeCalculator::width( int i ) const
 
 
 
-double Trk::MultiComponentStateModeCalculator::findRoot(double &result, double xlo, double xhi, double value, double i) const
+double Trk::MultiComponentStateModeCalculator::findRoot(const std::array<std::vector< Mixture >,5>& mixtureArray,
+                                                        double &result, double xlo, 
+                                                        double xhi, double value, 
+                                                        double i) const
 {
  // Do the root finding using the Brent-Decker method. Returns a boolean status and
  // loads 'result' with our best guess at the root if true.
@@ -381,8 +388,8 @@ double Trk::MultiComponentStateModeCalculator::findRoot(double &result, double x
 
  
   double a(xlo),b(xhi);
-  double fa= pdf(a, i) - value;
-  double fb= pdf(b, i) - value;
+  double fa= pdf(mixtureArray,a,i) - value;
+  double fb= pdf(mixtureArray,b,i) - value;
   
   if(fb*fa > 0) {
     ATH_MSG_DEBUG( "BrentRootFinder::findRoot: initial interval does not bracket a root: ("
@@ -479,7 +486,7 @@ double Trk::MultiComponentStateModeCalculator::findRoot(double &result, double x
     else {
       b += (m > 0 ? +tol : -tol);
     }
-    fb= pdf(b, i) - value;
+    fb= pdf(mixtureArray,b, i) - value;
 
   }
   // Return our best guess if we run out of iterations