Add Poisson Bootstrapping Capabilities

AmpTools/GPUManager/GPUManager.cc

@@ -345,8 +345,7 @@ GPUManager::copyDataToGPU( const AmpVecs& a, bool use4Vectors )
   #endif
 
   // copy the weights to the GPU
-  gpuErrChk( cudaMemcpy( m_pfDevWeights, a.m_pdWeights,
-                         m_iGDoubleDataArrSize, cudaMemcpyHostToDevice ) );
+  copyWeightsToGPU( a );
 
   // we only need to copy the four vectors to the GPU if the
   // amplitude evaulation kernels need them -- this is done by
@@ -426,6 +425,12 @@ GPUManager::copyUserVarsToGPU( const AmpVecs& a )
   #endif
 }
 
+void
+GPUManager::copyWeightsToGPU( const AmpVecs& a )
+{
+  gpuErrChk( cudaMemcpy( m_pfDevWeights, a.m_pdWeights,
+                         m_iGDoubleDataArrSize, cudaMemcpyHostToDevice ) );
+}
 
 void
 GPUManager::copyAmpsFromGPU( AmpVecs& a )

AmpTools/GPUManager/GPUManager.h

@@ -80,6 +80,7 @@ class GPUManager
 
   void copyDataToGPU( const AmpVecs& a, bool use4Vectors = true  );
   void copyUserVarsToGPU( const AmpVecs& a );
+  void copyWeightsToGPU( const AmpVecs& a );
 
   void calcAmplitudeAll( const Amplitude* amp, size_t uAmpFactOffset,
                          const vector< vector< int > >* pvPermutations,

AmpTools/IUAmpTools/AmpToolsInterface.cc

@@ -33,7 +33,8 @@
 // held liable for any liability with respect to any claim by the user or
 // any other party arising from use of the program.
 //******************************************************************************
-
+#include <ctime>
+#include <map>
 
 #include "IUAmpTools/AmpToolsInterface.h"
 #include "MinuitInterface/MinuitMinimizationManager.h"
@@ -400,6 +401,41 @@ AmpToolsInterface::randomizeParameter( const string& parName, float min, float m
   invalidateAmps();
 }
 
+void
+AmpToolsInterface::bootstrap( unsigned int seed ){
+
+  if( m_functionality != kFull ) return;
+  if ( seed == 0 ) seed = (unsigned int) time( NULL );
+
+  // give shared acc MC datasets across reactions the same bootstrap seed. Each will
+  // still draw and reassign its own set of weights, but they will be identical. Gen MC 
+  // datasets are never bootstrapped.
+  map< DataReader*, unsigned int > mcSeedForReader;
+  unsigned int nGroups = 0;
+
+  for (unsigned int irct = 0; irct < m_configurationInfo->reactionList().size(); irct++){
+
+    ReactionInfo* reaction = m_configurationInfo->reactionList()[irct];
+    string reactionName( reaction->reactionName() );
+
+    unsigned int mcSeed = 0;
+    if( DataReader* accMC = accMCReader( reactionName ) ){
+      map< DataReader*, unsigned int >::iterator it = mcSeedForReader.find( accMC );
+      if( it == mcSeedForReader.end() ){
+        mcSeed = seed + 1000000u * (++nGroups);
+        mcSeedForReader[accMC] = mcSeed;
+      }
+      else mcSeed = it->second;
+    }
+
+    // data/bkgnd files do not share common seeds across reactions
+    unsigned int dataSeed = seed + irct;
+
+    if (LikelihoodCalculator* likCalc = likelihoodCalculator( reactionName ) )
+      likCalc->applyPoissonWeights( dataSeed, mcSeed ); 
+  }
+}
+
 void
 AmpToolsInterface::finalizeFit( const string& tag ){
 

AmpTools/IUAmpTools/AmpToolsInterface.h

@@ -259,6 +259,15 @@ class AmpToolsInterface{
    */
 
   void randomizeParameter( const string& parName, float min = 0, float max = 1 );
+
+  /** This function will perform Poisson bootstrapping by re-weighting the 
+   * signal, accMC, and background (if available) events, using the provided 
+   * seed. Data and background are always given the same seed per reaction, while
+   * accMC has its own unique, but still traceable, seed. Repeated accMC samples
+   * use the same seed.
+   */
+
+  void bootstrap( unsigned int seed = 0 );
 
   /** Print final fit results to a file.  The tag can be used to
    *  generate a unique name in the case that multiple results are

AmpTools/IUAmpTools/AmpVecs.cc

@@ -59,8 +59,9 @@ AmpVecs::AmpVecs(){
   m_maxFactPerEvent = 0 ;
   m_userVarsPerEvent = 0;
 
-  m_pdData      = 0 ;
-  m_pdWeights   = 0 ;
+  m_pdData            = 0 ;
+  m_pdWeights         = 0 ;
+  m_pdOriginalWeights = 0 ;
 
   m_pdAmps       = 0 ;
   m_pdAmpFactors = 0 ;
@@ -104,6 +105,10 @@ AmpVecs::deallocAmpVecs()
   if(m_pdWeights)
     delete[] m_pdWeights;
   m_pdWeights=0;
+
+  if(m_pdOriginalWeights)
+    delete[] m_pdOriginalWeights;
+  m_pdOriginalWeights=0;
 
   deallocTerms( true );
 
@@ -206,6 +211,7 @@ AmpVecs::loadEvent( const Kinematics* pKinematics, size_t iEvent,
 
     m_pdData = new GDouble[4*m_iNParticles*m_iNEvents];
     m_pdWeights = new GDouble[m_iNEvents];        
+    m_pdOriginalWeights = new GDouble[m_iNEvents];
   }
 
   loadDataArrayElement( pKinematics, iEvent );
@@ -222,7 +228,7 @@ AmpVecs::loadData( DataReader* pDataReader, bool needsUserVarsOnly, size_t chunk
 
   //  Make sure no data is already loaded
 
-  if( m_pdData != NULL || m_pdWeights != NULL ){
+  if( m_pdData != NULL || m_pdWeights != NULL || m_pdOriginalWeights != NULL ){
 
     report( ERROR, kModule ) << "Trying to load data into a non-empty AmpVecs object\n"<<flush;
     assert(false);
@@ -276,6 +282,7 @@ AmpVecs::loadData( DataReader* pDataReader, bool needsUserVarsOnly, size_t chunk
       m_iNParticles = pKinematics->particleList().size();
       m_pdData = new GDouble[4*m_iNParticles*m_iNEvents];
       m_pdWeights = new GDouble[m_iNEvents]; 
+      m_pdOriginalWeights = new GDouble[m_iNEvents];
 
       #ifdef GPU_ACCELERATION
 
@@ -337,6 +344,26 @@ AmpVecs::loadDataArrayElement( const Kinematics* pKinematics, size_t iEvent ){
   }
 
   m_pdWeights[iEvent] = pKinematics->weight();
+  m_pdOriginalWeights[iEvent] = pKinematics->weight();
+}
+
+void
+AmpVecs::applyPoissonWeights( const int* weights ){
+
+  m_dSumWeights = 0;
+
+  for(size_t iEvent = 0; iEvent < m_iNTrueEvents; iEvent++){
+
+    m_pdWeights[iEvent] = m_pdOriginalWeights[iEvent] * weights[iEvent];
+    m_dSumWeights += m_pdWeights[iEvent];
+  }
+
+  m_integralValid = false;
+
+  #ifdef GPU_ACCELERATION
+    m_gpuMan.copyWeightsToGPU( *this );
+  #endif
+
 }
 
 void
@@ -522,6 +549,9 @@ AmpVecs::shareDataWith( AmpVecs* targetAmpVecs, bool needsUserVarsOnly ){
   targetAmpVecs->m_pdWeights = new GDouble[m_iNEvents];
   memcpy( targetAmpVecs->m_pdWeights, m_pdWeights, 
 	  sizeof(GDouble)*m_iNEvents );
+  targetAmpVecs->m_pdOriginalWeights = new GDouble[m_iNEvents];
+  memcpy( targetAmpVecs->m_pdOriginalWeights, m_pdOriginalWeights,
+    sizeof(GDouble)*m_iNEvents );
 
   targetAmpVecs->m_dataLoaded = true;
   targetAmpVecs->m_hasNonUnityWeights = m_hasNonUnityWeights;

AmpTools/IUAmpTools/AmpVecs.h

@@ -81,7 +81,7 @@ struct AmpVecs
   size_t m_iNTrueEvents;
 
   /**
-   * A double that stores the absolute value of the sum of the weights.  (For
+   * A double that stores the absolute value of the sum of the current weights.  (For
    * cases where all weights are unity, this is simply the number of true 
    * events.)
    */
@@ -118,10 +118,16 @@ struct AmpVecs
   GDouble* m_pdData;
 
   /**
-   * An array of length iNEvents that stores the event weights for each event.
+
+  * An array of length iNEvents that stores the current event weights for each event.
    */
   GDouble* m_pdWeights;
 
+  /**
+   * An array of length iNEvents that stores the original event weights for each event.
+   */
+  GDouble* m_pdOriginalWeights;
+
   /**
    * An array of length 2 * iNAmps * iNEvents that stores the real and imaginary
    * parts of the complete decay amplitude (product of factors) for each event.
@@ -296,6 +302,18 @@ struct AmpVecs
    */
   void loadEvent( const Kinematics* pKinematics, size_t iEvent = 0,
                   size_t iNTrueEvents = 1, bool needsUserVarsOnly = false );
+
+  /**
+   * This will overwrite the m_pdWeights array with the original weights multiplied by 
+   * an integer randomly drawn from a Poisson distribution of mean 1. For any given 
+   * weighted event, this will effectively
+   * - (integer == 0) remove the event 
+   * - (integer == 1) keep the event 
+   * - (integer == N) repeat the event N times, for N > 1
+   * 
+   * \param[in] weights an array of integers of length m_iNTrueEvents
+   */
+  void applyPoissonWeights(const int* weights);
 
   /**
    * A helper routine to get an event i from the array of data and weights.

AmpTools/IUAmpTools/AmplitudeManager.cc

@@ -744,14 +744,16 @@ SCOREP_USER_REGION_BEGIN( calcSumLogIntensity, "calcSumLogIntensity", SCOREP_USE
   calcIntensities( a );
 
   for( int iEvent=0; iEvent < a.m_iNTrueEvents; iEvent++ ){
+
+    double weight = a.m_pdWeights[iEvent];
+    if( weight == 0 ) continue; // avoid divide by zero
 
     // here divide out the weight that was put into the intensity calculation
     // and weight the log -- in practice this just contributes an extra constant
     // term in the likelihood equal to sum -w_i * log( w_i ), but the division
     // helps avoid problems with negative weights, which may be used
     // in background subtraction
-    dSumLogI += a.m_pdWeights[iEvent] *
-    G_LOG( a.m_pdIntensity[iEvent] / a.m_pdWeights[iEvent] );
+    dSumLogI += weight * G_LOG( a.m_pdIntensity[iEvent] / weight );
   }
 
 #else

AmpTools/IUAmpTools/LikelihoodCalculator.cc

@@ -38,6 +38,8 @@
 #include <cassert>
 #include <string>
 
+#include "TRandom3.h"
+
 #include "IUAmpTools/LikelihoodCalculator.h"
 #include "IUAmpTools/IntensityManager.h"
 #include "IUAmpTools/DataReader.h"
@@ -305,6 +307,28 @@ SCOREP_USER_REGION_END( dataTerm )
   return sumLnI;
 }
 
+void
+LikelihoodCalculator::applyPoissonWeights( unsigned int dataSeed, unsigned int mcSeed ){
+  if( m_firstDataCalc ) dataTerm();
+
+  TRandom3 rndGen( dataSeed );
+
+  vector<int> sigW( m_ampVecsSignal.m_iNTrueEvents);
+  for( size_t i = 0; i < sigW.size(); ++i) sigW[i] = rndGen.Poisson( 1.0 );
+  m_ampVecsSignal.applyPoissonWeights( &sigW[0] );
+  m_sumDataWeights = m_ampVecsSignal.m_dSumWeights;
+
+  if( m_hasBackground ){
+    vector<int> bkgW( m_ampVecsBkgnd.m_iNTrueEvents);
+    for( size_t i = 0; i < bkgW.size(); ++i) bkgW[i] = rndGen.Poisson( 1.0 );
+    m_ampVecsBkgnd.applyPoissonWeights( &bkgW[0] );
+    m_sumBkgWeights = m_ampVecsBkgnd.m_dSumWeights;
+  }
+
+  m_normInt.applyPoissonWeights( mcSeed );
+
+}
+
 void
 LikelihoodCalculator::invalidateTerms(){
 

AmpTools/IUAmpTools/LikelihoodCalculator.h

@@ -85,6 +85,7 @@ class LikelihoodCalculator : public MIFunctionContribution
   double operator()();
 
   virtual double numSignalEvents();
+  virtual void applyPoissonWeights( unsigned int dataSeed, unsigned int mcSeed );
 
   void invalidateTerms();
 
@@ -108,13 +109,14 @@ class LikelihoodCalculator : public MIFunctionContribution
   void setNumBkgEvents (  double num ) { m_numBkgEvents   = num; }
   void setSumDataWeights(  double sum ) { m_sumDataWeights = sum; }
   void setNumDataEvents(  double num ) { m_numDataEvents  = num; }
+
+  const NormIntInterface& m_normInt;
 
 private:
 
   bool m_hasBackground;
 
   const IntensityManager& m_intenManager;
-  const NormIntInterface& m_normInt;
 
   DataReader* m_dataReaderSignal;
   DataReader* m_dataReaderBkgnd;

AmpTools/IUAmpTools/NormIntInterface.cc

@@ -44,6 +44,8 @@
 #include <cstring>
 #include <iomanip>
 
+#include "TRandom3.h"
+
 #include "IUAmpTools/NormIntInterface.h"
 #include "IUAmpTools/report.h"
 const char* NormIntInterface::kModule = "NormIntInterface";
@@ -603,6 +605,20 @@ NormIntInterface::loadMC() const {
   m_sumGenWeights = m_genMCVecs.m_dSumWeights;
 }
 
+void
+NormIntInterface::applyPoissonWeights( unsigned int seed ) const {
+  if( !m_accMCVecs.m_dataLoaded ) loadMC();
+
+  TRandom3 rndGen(seed);
+
+  vector<int> accW( m_accMCVecs.m_iNTrueEvents );
+  for( size_t i = 0; i < accW.size(); ++i ) accW[i] = rndGen.Poisson(1.0);  
+  m_accMCVecs.applyPoissonWeights( &accW[0] );
+  m_sumAccWeights = m_accMCVecs.m_dSumWeights;
+
+  m_emptyNormIntCache = true;
+}
+
 void
 NormIntInterface::invalidateTerms(){
 

AmpTools/IUAmpTools/NormIntInterface.h

@@ -88,6 +88,7 @@ class NormIntInterface
   // needs to be virtual so parallel implementations can properly
   // override this function
   virtual void forceCacheUpdate( bool normIntOnly = false ) const;
+  virtual void applyPoissonWeights( unsigned int seed ) const;
 
   void invalidateTerms();
 
@@ -103,8 +104,8 @@ class NormIntInterface
   const double* ampIntMatrix() const  { return m_ampIntCache;  }
   const double* normIntMatrix() const { return m_normIntCache; }
 
-  void setGenEvents( double sumWeights ) { m_sumGenWeights = sumWeights; }
-  void setAccEvents( double sumWeights ) { m_sumAccWeights = sumWeights; }
+  void setGenEvents( double sumWeights ) const { m_sumGenWeights = sumWeights; }
+  void setAccEvents( double sumWeights ) const { m_sumAccWeights = sumWeights; }
 
 protected:
 

AmpTools/IUAmpToolsMPI/LikelihoodCalculatorMPI.cc

@@ -276,6 +276,35 @@ LikelihoodCalculatorMPI::computeLikelihood()
   MPI_Send( (double*)&data, 5, MPI_DOUBLE, 0, MPITag::kDoubleSend, MPI_COMM_WORLD );
 }
 
+void
+LikelihoodCalculatorMPI::applyPoissonWeights( unsigned int dataSeed, unsigned int mcSeed ){
+  assert( m_isLeader );
+
+  int cmnd[2] = { m_thisId, LikelihoodManagerMPI::kApplyPoissonWeights };
+  MPI_Bcast( cmnd, 2, MPI_INT, 0, MPI_COMM_WORLD );
+
+  unsigned int seeds[2] = { dataSeed, mcSeed };
+  MPI_Bcast( seeds, 2, MPI_UNSIGNED, 0, MPI_COMM_WORLD );
+
+  // despite the leader holding no sig/bkgd data for reweighting, we need to join 
+  // the normIntInterface MPI handshakes to avoid a desync.
+  m_normInt.applyPoissonWeights( mcSeed );
+
+}
+
+void
+LikelihoodCalculatorMPI::applyBootstrapWeights(){
+  assert( !m_isLeader );
+
+  unsigned int seeds[2];
+  MPI_Bcast( seeds, 2, MPI_UNSIGNED, 0, MPI_COMM_WORLD );
+
+  // give each follower an independent generator
+  unsigned int dataSeed = seeds[0] + 1000003u * (unsigned int) m_rank;
+  unsigned int mcSeed = seeds[1] + 1000003u * (unsigned int) m_rank;
+  LikelihoodCalculator::applyPoissonWeights( dataSeed, mcSeed );
+}
+
 void
 LikelihoodCalculatorMPI::setupMPI()
 {