Add Capability to Bootstrap Signal and/or Background Datasets

AmpTools/IUAmpTools/AmpToolsInterface.cc

@@ -170,7 +170,7 @@ AmpToolsInterface::resetConfigurationInfo(ConfigurationInfo* configurationInfo){
       for (unsigned int i = 0; i < m_userDataReaders.size(); i++){
         if (reaction->data().first == m_userDataReaders[i]->name())
           m_dataReaderMap[reactionName]
-	  = m_userDataReaders[i]->newDataReader(reaction->data().second);
+          = m_userDataReaders[i]->newDataReader(reaction->data().second);
         if (reaction->bkgnd().first == m_userDataReaders[i]->name())
           m_bkgndReaderMap[reactionName]
           = m_userDataReaders[i]->newDataReader(reaction->bkgnd().second);
@@ -199,7 +199,7 @@ AmpToolsInterface::resetConfigurationInfo(ConfigurationInfo* configurationInfo){
         report( WARNING, kModule ) << "not creating a DataReader for accepted MC associated with reaction " << reactionName << endl;
 
       if( m_functionality == kFull ){
- 
+
         // ************************
         // create a NormIntInterface
         // ************************
@@ -400,6 +400,34 @@ AmpToolsInterface::randomizeParameter( const string& parName, float min, float m
   invalidateAmps();
 }
 
+void
+AmpToolsInterface::bootstrapSignalData(const string& reactionName, unsigned int seed){
+
+  LikelihoodCalculator* likCalc = likelihoodCalculator(reactionName);
+  if( likCalc == NULL ){
+    report( ERROR, kModule ) << "no LikelihoodCalculator for reaction: " << reactionName << endl;
+    return;
+  }
+  likCalc->bootstrapSignalData( seed );
+
+}
+
+void
+AmpToolsInterface::bootstrapBackgroundData(const string& reactionName, unsigned int seed){
+
+  LikelihoodCalculator* likCalc = likelihoodCalculator(reactionName);
+  DataReader* backgroundReader = bkgndReader(reactionName);
+  if( likCalc == NULL ){
+    report( ERROR, kModule ) << "no LikelihoodCalculator for reaction: " << reactionName << endl;
+    return;
+  }
+  if( backgroundReader == NULL ){
+    report( ERROR, kModule ) << "no background data for reaction: " << reactionName << endl;
+    return;
+  }
+  likCalc->bootstrapBackgroundData( seed );
+}
+
 void
 AmpToolsInterface::finalizeFit( const string& tag ){
 

AmpTools/IUAmpTools/AmpToolsInterface.h

@@ -44,6 +44,7 @@
 #include <map>
 #include <complex>
 #include <fstream>
+#include <ctime>
 
 #include "MinuitInterface/MinuitMinimizationManager.h"
 #include "IUAmpTools/IntensityManager.h"
@@ -259,18 +260,34 @@ class AmpToolsInterface{
    */
 
   void randomizeParameter( const string& parName, float min = 0, float max = 1 );
+
+  /** This function will call the dataReader.resample() method intended to 
+   *  randomly sample, with replacement, the signal events before getEvent()
+   *  is called. This is particularly useful for fits that contain free 
+   *  parameters in the amplitudes themselves, as the bootstrap distributions 
+   *  will provide better uncertainty estimates. If no seed value is given, a
+   *  random one is chosen.
+   */
+
+  void bootstrapSignalData(const string& reactionName, unsigned int seed = static_cast<unsigned int>( time(NULL) ) );
+
+  /** This function will call the dataReader.resample() method intended to 
+   *  randomly sample, with replacement, the background events before getEvent()
+   *  is called. If no seed value is given, a random one is chosen.
+   */
+  void bootstrapBackgroundData(const string& reactionName, unsigned int seed = static_cast<unsigned int>( time(NULL) ) );
 
   /** Print final fit results to a file.  The tag can be used to
    *  generate a unique name in the case that multiple results are
-   *  written for a singele fit job.
+   *  written for a single fit job.
    */
 
   virtual void finalizeFit( const string& tag = "" );
 
 
   /** For manual calculations:  clear all events and calculations.
    *  Call this before loading events from a new reaction or to start
-   *  new calcuations.
+   *  new calculations.
    *
    * \param[in] iDataSet used to index simultaneous manual calculations
    *

AmpTools/IUAmpTools/DataReader.cc

@@ -39,6 +39,7 @@
 #include <vector>
 
 #include "IUAmpTools/DataReader.h"
+const char* DataReader::kModule = "DataReader";
 
 using namespace std;
 

AmpTools/IUAmpTools/DataReader.h

@@ -39,6 +39,9 @@
 
 #include <string>
 #include <vector>
+#include <cassert>
+
+#include "IUAmpTools/report.h"
 
 using namespace std;
 
@@ -105,6 +108,24 @@ class DataReader
    * \see DataReaderMPI
    */
   virtual void resetSource() = 0; 
+
+  /**
+   * Reserved for DataReaders that support redrawing samples in place (e.g. 
+   * bootstrap resampling with replacement) without reconstructing the object.
+   * The default implementation reports that no resampling is supported. 
+   * Overridden by a reader that does.
+   * 
+   * This method should be virtual in the user's class if the DataReaderMPI 
+   * template is to be used.   
+   *    
+   * 
+   * \see DataReaderMPI
+   */
+  virtual void resample( unsigned int seed ){
+    report( ERROR, kModule ) << name() << " does not implement resample(). "
+    << "Cannot bootstrap this data source." << endl;
+    assert( false );
+  }
 
   /**
    * The user should override this function with one that returns the number
@@ -166,7 +187,7 @@ class DataReader
 
   vector<string> m_args;
 
-
+  static const char* kModule;
 };
 
 #endif

AmpTools/IUAmpTools/LikelihoodCalculator.cc

@@ -67,6 +67,7 @@ m_normInt( normInt ),
 m_dataReaderSignal( dataReaderSignal ),
 m_dataReaderBkgnd( dataReaderBkgnd ),
 m_firstDataCalc( true ),
+m_firstBkgndCalc( true ),
 m_firstNormIntCalc( true ),
 m_sumBkgWeights( 0 ),
 m_numBkgEvents( 0 ),
@@ -106,7 +107,7 @@ LikelihoodCalculator::operator()(){
 double
 LikelihoodCalculator::numSignalEvents(){
 
-  if( m_firstDataCalc ) dataTerm();
+  if( m_firstDataCalc || m_firstBkgndCalc ) dataTerm();
   return numDataEvents() - sumBkgWeights();
 }
 
@@ -230,7 +231,7 @@ SCOREP_USER_REGION_BEGIN( dataTerm, "dataTerm", SCOREP_USER_REGION_TYPE_COMMON )
 
     m_numDataEvents = m_ampVecsSignal.m_iNTrueEvents;
     m_sumDataWeights = m_ampVecsSignal.m_dSumWeights;
- 
+
     if( m_ampVecsSignal.m_hasNonUnityWeights && m_hasBackground ){
 
       report( WARNING, kModule ) << "\n"
@@ -244,44 +245,45 @@ SCOREP_USER_REGION_BEGIN( dataTerm, "dataTerm", SCOREP_USER_REGION_TYPE_COMMON )
       << "*   contribution to the likelihood.                            *\n"
       << "****************************************************************\n" << endl;
     }
+  }
 
-    if( m_hasBackground ){
+  if( m_firstBkgndCalc && m_hasBackground ){
 
-      m_ampVecsBkgnd.loadData( m_dataReaderBkgnd, m_intenManager.needsUserVarsOnly() );
-      m_ampVecsBkgnd.allocateTerms( m_intenManager, true );
+    m_ampVecsBkgnd.loadData( m_dataReaderBkgnd, m_intenManager.needsUserVarsOnly() );
+    m_ampVecsBkgnd.allocateTerms( m_intenManager, true );
 
-      if( m_ampVecsBkgnd.m_hasMixedSignWeights ){
-        report( NOTICE, kModule ) << "***************************************************************" << endl;
-        report( NOTICE, kModule ) << "* NOTICE:  Weights with both positive and negative signs were *" << endl;
-        report( NOTICE, kModule ) << "* detected in the background file.  This may be desirable for *" << endl;
-        report( NOTICE, kModule ) << "* some applications.  Older versions of AmpTools (v0.10.x and *" << endl;
-        report( NOTICE, kModule ) << "* prior) will not properly handle this case and will also not *" << endl;
-        report( NOTICE, kModule ) << "* print this notification to the screen.                      *" << endl;
-        report( NOTICE, kModule ) << "***************************************************************" << endl;
-      }
-
-      m_sumBkgWeights = m_ampVecsBkgnd.m_dSumWeights;
-
-      // the extra boolean allows MPI jobs to suppress this check which
-      // may fail on one of the follower nodes if the background sample
-      // is sparse
-      if( m_sumBkgWeights < 0 && !suppressError ){
-        report( ERROR, kModule ) << "****************************************************************" << endl;
-        report( ERROR, kModule ) << "* ERROR: The sum of all background weights is negative.  This  *" << endl;
-        report( ERROR, kModule ) << "*   implies a negative background in the signal region, which  *" << endl;
-        report( ERROR, kModule ) << "*   unphysical.  The weighted sum of the background events     *" << endl;
-        report( ERROR, kModule ) << "*   should represent the background contribution to the signal *" << endl;
-        report( ERROR, kModule ) << "*   region.                                                    *" << endl;
-        report( ERROR, kModule ) << "****************************************************************" << endl;
-
-        assert( false );
-      }
+    if( m_ampVecsBkgnd.m_hasMixedSignWeights ){
+      report( NOTICE, kModule ) << "***************************************************************" << endl;
+      report( NOTICE, kModule ) << "* NOTICE:  Weights with both positive and negative signs were *" << endl;
+      report( NOTICE, kModule ) << "* detected in the background file.  This may be desirable for *" << endl;
+      report( NOTICE, kModule ) << "* some applications.  Older versions of AmpTools (v0.10.x and *" << endl;
+      report( NOTICE, kModule ) << "* prior) will not properly handle this case and will also not *" << endl;
+      report( NOTICE, kModule ) << "* print this notification to the screen.                      *" << endl;
+      report( NOTICE, kModule ) << "***************************************************************" << endl;
+    }
+
+    m_sumBkgWeights = m_ampVecsBkgnd.m_dSumWeights;
+
+    // the extra boolean allows MPI jobs to suppress this check which
+    // may fail on one of the follower nodes if the background sample
+    // is sparse
+    if( m_sumBkgWeights < 0 && !suppressError ){
+      report( ERROR, kModule ) << "****************************************************************" << endl;
+      report( ERROR, kModule ) << "* ERROR: The sum of all background weights is negative.  This  *" << endl;
+      report( ERROR, kModule ) << "*   implies a negative background in the signal region, which  *" << endl;
+      report( ERROR, kModule ) << "*   unphysical.  The weighted sum of the background events     *" << endl;
+      report( ERROR, kModule ) << "*   should represent the background contribution to the signal *" << endl;
+      report( ERROR, kModule ) << "*   region.                                                    *" << endl;
+      report( ERROR, kModule ) << "****************************************************************" << endl;
 
-      m_numBkgEvents = m_ampVecsBkgnd.m_iNTrueEvents;
+      assert( false );
     }
 
-    report( DEBUG, kModule ) << "\tDone." << endl;
+    m_numBkgEvents = m_ampVecsBkgnd.m_iNTrueEvents;
+    m_firstBkgndCalc = false;
   }
+
+  report( DEBUG, kModule ) << "\tDone." << endl;
 
   double sumLnI = m_intenManager.calcSumLogIntensity( m_ampVecsSignal );
 
@@ -305,6 +307,38 @@ SCOREP_USER_REGION_END( dataTerm )
   return sumLnI;
 }
 
+void
+LikelihoodCalculator::bootstrapSignalData( unsigned int seed ){
+#ifdef SCOREP
+SCOREP_USER_REGION_DEFINE( bootstrapSignalData )                                                                                    
+SCOREP_USER_REGION_BEGIN( bootstrapSignalData, "bootstrapSignalData", SCOREP_USER_REGION_TYPE_COMMON )
+#endif
+  m_ampVecsSignal.deallocAmpVecs();
+  m_dataReaderSignal->resample( seed );
+  m_firstDataCalc = true;
+#ifdef SCOREP
+SCOREP_USER_REGION_END( bootstrapSignalData )
+#endif
+}
+
+void
+LikelihoodCalculator::bootstrapBackgroundData( unsigned int seed ){
+#ifdef SCOREP
+SCOREP_USER_REGION_DEFINE( bootstrapBackgroundData )                                                                                    
+SCOREP_USER_REGION_BEGIN( bootstrapBackgroundData, "bootstrapBackgroundData", SCOREP_USER_REGION_TYPE_COMMON )
+#endif
+  if( !m_hasBackground ){
+    report( ERROR, kModule ) << "Requested to bootstrap sample non-existant background dataset" << endl;
+    assert( false );
+  }
+  m_ampVecsBkgnd.deallocAmpVecs();
+  m_dataReaderBkgnd->resample( seed );
+  m_firstBkgndCalc = true;
+#ifdef SCOREP
+SCOREP_USER_REGION_END( bootstrapBackgroundData )
+#endif
+}
+
 void
 LikelihoodCalculator::invalidateTerms(){
 

AmpTools/IUAmpTools/LikelihoodCalculator.h

@@ -87,6 +87,9 @@ class LikelihoodCalculator : public MIFunctionContribution
   virtual double numSignalEvents();
 
   void invalidateTerms();
+
+  virtual void bootstrapSignalData( unsigned int seed );
+  virtual void bootstrapBackgroundData( unsigned int seed);
 
 protected:
 
@@ -121,6 +124,7 @@ class LikelihoodCalculator : public MIFunctionContribution
 
   bool m_firstNormIntCalc;
   bool m_firstDataCalc;
+  bool m_firstBkgndCalc;
 
   double* m_prodFactorArray;
   const double* m_normIntArray;

AmpTools/IUAmpToolsMPI/DataReaderMPI.h

@@ -31,8 +31,8 @@ struct KinStruct {
  * of the instances of this class on the follower nodes will behave as if
  * they have only a subset of the data.  The instance of it on the leader
  * node will behave as if it has all of the data.  Be sure that the getEvent, 
- * resetSource, and numEvents methods in the user-defined class are declared
- * virtual.
+ * resetSource, numEvents, and (if used for bootstrapping) resample methods in
+ * the user-defined class are declared virtual.
  *
  * \ingroup IUAmpToolsMPI
  */
@@ -66,6 +66,7 @@ class DataReaderMPI : public T
   Kinematics* getEvent();
 
   void resetSource();
+  void resample( unsigned int seed );
 
   unsigned int numEvents() const;
 
@@ -229,6 +230,35 @@ void DataReaderMPI<T>::resetSource()
   }
 }
 
+template< class T >
+void DataReaderMPI<T>::resample( unsigned int seed )
+{
+  if( m_isLeader ){
+
+    report( DEBUG, kDRModule ) << "Resampling leader data source " << m_rank << endl;
+
+    // redraw the bootstrap sample using the user's definition in the reader, then
+    // push the selection out to the follower's exactly as done in startup
+    T::resample( seed );
+    distributeData();
+  }
+  else{
+
+    report( DEBUG, kDRModule ) << "Receiving resampled data on process with rank " << m_rank << endl;
+
+    //discard the cached partition
+    for( vector<Kinematics*>::iterator ptrItr = m_ptrCache.begin();
+        ptrItr != m_ptrCache.end();
+        ++ptrItr ){
+
+      delete *ptrItr;
+    }
+    m_ptrCache.clear();
+
+    receiveData();
+  }
+}
+
 template< class T >
 void DataReaderMPI<T>::provideData()
 {