diff --git a/AmpTools/IUAmpTools/AmpToolsInterface.cc b/AmpTools/IUAmpTools/AmpToolsInterface.cc index d4cb3c4..5a3a21b 100644 --- a/AmpTools/IUAmpTools/AmpToolsInterface.cc +++ b/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 ){ diff --git a/AmpTools/IUAmpTools/AmpToolsInterface.h b/AmpTools/IUAmpTools/AmpToolsInterface.h index 24c7b8a..914131a 100644 --- a/AmpTools/IUAmpTools/AmpToolsInterface.h +++ b/AmpTools/IUAmpTools/AmpToolsInterface.h @@ -44,6 +44,7 @@ #include #include #include +#include #include "MinuitInterface/MinuitMinimizationManager.h" #include "IUAmpTools/IntensityManager.h" @@ -259,10 +260,26 @@ 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( 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( 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 = "" ); @@ -270,7 +287,7 @@ class AmpToolsInterface{ /** 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 * diff --git a/AmpTools/IUAmpTools/DataReader.cc b/AmpTools/IUAmpTools/DataReader.cc index 0e0a240..8b3a44a 100644 --- a/AmpTools/IUAmpTools/DataReader.cc +++ b/AmpTools/IUAmpTools/DataReader.cc @@ -39,6 +39,7 @@ #include #include "IUAmpTools/DataReader.h" +const char* DataReader::kModule = "DataReader"; using namespace std; diff --git a/AmpTools/IUAmpTools/DataReader.h b/AmpTools/IUAmpTools/DataReader.h index 19000fc..72be429 100644 --- a/AmpTools/IUAmpTools/DataReader.h +++ b/AmpTools/IUAmpTools/DataReader.h @@ -39,6 +39,9 @@ #include #include +#include + +#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 m_args; - + static const char* kModule; }; #endif diff --git a/AmpTools/IUAmpTools/LikelihoodCalculator.cc b/AmpTools/IUAmpTools/LikelihoodCalculator.cc index af65939..af6f328 100644 --- a/AmpTools/IUAmpTools/LikelihoodCalculator.cc +++ b/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(){ diff --git a/AmpTools/IUAmpTools/LikelihoodCalculator.h b/AmpTools/IUAmpTools/LikelihoodCalculator.h index 62df8c0..bc9cd0f 100644 --- a/AmpTools/IUAmpTools/LikelihoodCalculator.h +++ b/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; diff --git a/AmpTools/IUAmpToolsMPI/DataReaderMPI.h b/AmpTools/IUAmpToolsMPI/DataReaderMPI.h index a98b88e..23e01d6 100644 --- a/AmpTools/IUAmpToolsMPI/DataReaderMPI.h +++ b/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::resetSource() } } +template< class T > +void DataReaderMPI::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::iterator ptrItr = m_ptrCache.begin(); + ptrItr != m_ptrCache.end(); + ++ptrItr ){ + + delete *ptrItr; + } + m_ptrCache.clear(); + + receiveData(); + } +} + template< class T > void DataReaderMPI::provideData() { diff --git a/AmpTools/IUAmpToolsMPI/LikelihoodCalculatorMPI.cc b/AmpTools/IUAmpToolsMPI/LikelihoodCalculatorMPI.cc index cfcacf8..f71ad54 100644 --- a/AmpTools/IUAmpToolsMPI/LikelihoodCalculatorMPI.cc +++ b/AmpTools/IUAmpToolsMPI/LikelihoodCalculatorMPI.cc @@ -215,6 +215,46 @@ LikelihoodCalculatorMPI::operator()() return -2 * lnL; } +void +LikelihoodCalculatorMPI::bootstrapSignalData( unsigned int seed ){ + + if( m_isLeader ){ + + int cmnd[2]; + cmnd[0] = m_thisId; + cmnd[1] = LikelihoodManagerMPI::kBootstrapSignalData; + MPI_Bcast( cmnd, 2, MPI_INT, 0, MPI_COMM_WORLD); + + // second, command-specific broadcast carrying the randomized seed for bootstrapping + unsigned int seedBuf = seed; + MPI_Bcast( &seedBuf, 1, MPI_UNSIGNED, 0, MPI_COMM_WORLD ); + } + + // calls resample() on DataReaderMPI for the leader, and redistributes to the + // followers, discarding the old cached partition + LikelihoodCalculator::bootstrapSignalData( seed ); +} + +void +LikelihoodCalculatorMPI::bootstrapBackgroundData( unsigned int seed ){ + + if( m_isLeader ){ + + int cmnd[2]; + cmnd[0] = m_thisId; + cmnd[1] = LikelihoodManagerMPI::kBootstrapBackgroundData; + MPI_Bcast( cmnd, 2, MPI_INT, 0, MPI_COMM_WORLD); + + // second, command-specific broadcast carrying the randomized seed for bootstrapping + unsigned int seedBuf = seed; + MPI_Bcast( &seedBuf, 1, MPI_UNSIGNED, 0, MPI_COMM_WORLD ); + } + + // calls resample() on DataReaderMPI for the leader, and redistributes to the + // followers, discarding the old cached partition + LikelihoodCalculator::bootstrapBackgroundData( seed ); +} + double LikelihoodCalculatorMPI::numSignalEvents(){ diff --git a/AmpTools/IUAmpToolsMPI/LikelihoodCalculatorMPI.h b/AmpTools/IUAmpToolsMPI/LikelihoodCalculatorMPI.h index eb88cd6..3a4dab5 100644 --- a/AmpTools/IUAmpToolsMPI/LikelihoodCalculatorMPI.h +++ b/AmpTools/IUAmpToolsMPI/LikelihoodCalculatorMPI.h @@ -123,6 +123,20 @@ class LikelihoodCalculatorMPI : public LikelihoodCalculator */ double operator()(); + /** + * Triggers a bootstrap resampling of the signal data. On a leader it + * broadcasts the command to all followers and resamples the leader's own + * DataReaderMPI, pushing newly-drawn events to the followers. + */ + void bootstrapSignalData( unsigned int seed ); + + /** + * Triggers a bootstrap resampling of the background data. On a leader it + * broadcasts the command to all followers and resamples the leader's own + * DataReaderMPI, pushing newly-drawn events to the followers. + */ + void bootstrapBackgroundData( unsigned int seed ); + /** * This sends the finalize fit flag to all follower jobs which breaks them * out of the deliverLikelihood method of the LikelihoodManagerMPI. diff --git a/AmpTools/IUAmpToolsMPI/LikelihoodManagerMPI.cc b/AmpTools/IUAmpToolsMPI/LikelihoodManagerMPI.cc index 3582a22..dd3aab2 100644 --- a/AmpTools/IUAmpToolsMPI/LikelihoodManagerMPI.cc +++ b/AmpTools/IUAmpToolsMPI/LikelihoodManagerMPI.cc @@ -72,6 +72,7 @@ LikelihoodManagerMPI::deliverLikelihood() int* fitFlag = &(cmnd[1]); LikelihoodCalculatorMPI* likCalc; + unsigned int seed; map< int, LikelihoodCalculatorMPI* >::iterator mapItr; MPI_Bcast( cmnd, 2, MPI_INT, 0, MPI_COMM_WORLD ); @@ -106,6 +107,20 @@ LikelihoodManagerMPI::deliverLikelihood() likCalc->computeLikelihood(); break; + + case kBootstrapSignalData: + + MPI_Bcast( &seed, 1, MPI_UNSIGNED, 0, MPI_COMM_WORLD ); + + likCalc->bootstrapSignalData( seed ); + break; + + case kBootstrapBackgroundData: + + MPI_Bcast( &seed, 1, MPI_UNSIGNED, 0, MPI_COMM_WORLD ); + + likCalc->bootstrapBackgroundData( seed ); + break; default: diff --git a/AmpTools/IUAmpToolsMPI/LikelihoodManagerMPI.h b/AmpTools/IUAmpToolsMPI/LikelihoodManagerMPI.h index 5d1c79e..045f9df 100644 --- a/AmpTools/IUAmpToolsMPI/LikelihoodManagerMPI.h +++ b/AmpTools/IUAmpToolsMPI/LikelihoodManagerMPI.h @@ -60,6 +60,8 @@ class LikelihoodManagerMPI kComputeIntegrals, kUpdateParameters, kUpdateAmpParameter, + kBootstrapSignalData, + kBootstrapBackgroundData, kFinalizeFit, kExit };