eic-opticks/u4/U4Random.cc

0001 #include <iostream>
0002 #include <iomanip>
0003 #include <sstream>
0004 #include <cstdlib>
0005 #include <csignal>
0006
0007 #include "Randomize.hh"
0008 #include "G4Types.hh"
0009 #include "U4Random.hh"
0010 #include "U4Stack.h"
0011 #include "U4StackAuto.h"
0012
0013 #include "NP.hh"
0014 #include "spath.h"
0015 #include "ssys.h"
0016
0017 #include "SEvt.hh"
0018 #include "SBacktrace.h"
0019 #include "SDBG.h"
0020
0021 #include "SLOG.hh"
0022
0023 const plog::Severity U4Random::LEVEL = SLOG::EnvLevel("U4Random", "DEBUG") ;
0024
0025 U4Random* U4Random::INSTANCE = nullptr ;
0026 U4Random* U4Random::Get(){ return INSTANCE ; }
0027
0028 const char* U4Random::NAME = "U4Random" ;
0029
0030
0031 void U4Random::SetSequenceIndex(int index)
0032 {
0033     if(U4Random::Get() == nullptr) return ;
0034
0035     U4Random* rnd = U4Random::Get();
0036     if(rnd->isReady() == false)
0037     {
0038         LOG(LEVEL) << " index " << index << " NOT READY " ;
0039         return ;
0040     }
0041     rnd->setSequenceIndex(index);
0042 }
0043
0044 int U4Random::GetSequenceIndex()
0045 {
0046     if(U4Random::Get() == nullptr) return -1 ;
0047     U4Random* rnd = U4Random::Get();
0048     return rnd->isReady() ? rnd->getSequenceIndex() : -1 ;
0049 }
0050
0051
0052
0053
0054 const char* U4Random::SeqPath(){ return ssys::getenvvar(OPTICKS_RANDOM_SEQPATH, DEFAULT_SEQPATH ) ; } // static
0055
0056
0057 U4Random* U4Random::Create(const char* seq, const char* seqmask)
0058 {
0059     U4Random* rdm = new U4Random(seq, seqmask);
0060     bool ready = rdm->isReady();
0061     LOG_IF(fatal, !ready) << " NOT READY : \n" << rdm->detail() ;
0062     return ready ? rdm : nullptr  ;
0063 }
0064
0065
0066
0067 /**
0068 U4Random::U4Random
0069 -------------------------------
0070
0071 When no seq path argument is provided the envvar OPTICKS_RANDOM_SEQPATH
0072 is consulted to provide the path.
0073
0074 The optional seqmask (a list of size_t or "unsigned long long" indices)
0075 allows working with sub-selections of the full set of streams of randoms.
0076 This allows reproducible running within photon selections
0077 by arranging the same random stream to be consumed in
0078 full-sample and sub-sample running.
0079
0080 Not that *seq* can either be the path to an .npy file
0081 or the path to a directory containing .npy files which
0082 are concatenated using NP::Load/NP::Concatenate.
0083
0084 TODO: when next need to use seqmask, change the interface to enabling it
0085 to be envvar OR SEventConfig based and eliminate the U4Random arguments
0086
0087 **/
0088
0089 U4Random::U4Random(const char* seq, const char* seqmask)
0090     :
0091     m_seqpath_( seq ? seq : SeqPath() ),
0092     m_seqpath(spath::Resolve(m_seqpath_)),
0093     m_seqpath_exists( NP::Exists( m_seqpath ) || NP::ExistsArrayFolder(m_seqpath) ),
0094     m_seq(m_seqpath_exists ? NP::Load(m_seqpath) : nullptr),
0095     m_seq_values(m_seq ? m_seq->cvalues<float>() : nullptr ),
0096     m_seq_ni(m_seq ? m_seq->shape[0] : 0 ),                        // num items
0097     m_seq_nv(m_seq ? m_seq->shape[1]*m_seq->shape[2] : 0 ),        // num values in each item
0098     m_seq_index(-1),
0099
0100     m_cur(m_seq_ni > 0 ? NP::Make<int>(m_seq_ni) : nullptr),
0101     m_cur_values(m_cur ? m_cur->values<int>() : nullptr),
0102     m_recycle(true),
0103     m_default(CLHEP::HepRandom::getTheEngine()),
0104
0105     m_seqmask(seqmask ? NP::Load(seqmask) : nullptr),
0106     m_seqmask_ni( m_seqmask ? m_seqmask->shape[0] : 0 ),
0107     m_seqmask_values(m_seqmask ? m_seqmask->cvalues<size_t>() : nullptr),
0108     //m_flat_debug(ssys::getenvbool("U4Random_flat_debug")),
0109     m_flat_prior(0.),
0110     m_ready(false),
0111     m_select(ssys::getenvintvec("U4Random_select")),
0112     m_select_action(SDBG::Action(ssys::getenvvar("U4Random_select_action", "backtrace")))   // "backtrace" "caller" "interrupt" "summary"
0113 {
0114     init();
0115 }
0116
0117
0118 /**
0119 U4Random::isSelect
0120 ---------------------
0121
0122 Returns true when the (photon index, random index "flat cursor") tuple
0123 provided in the arguments matches one of the pairs provided in the
0124 comma delimited U4Random_select envvar.
0125
0126 An U4Random_select envvar of (-1,-1) is special cased to always match,
0127 and hence may generate very large amounts of output dumping.
0128 A single -1 corresponds to a wildcard:
0129
0130 * (0,-1) will match all cursors for photon idx 0
0131 * (-1,0) will match all cursor 0 for all photon idx
0132
0133 This allows selection of one or more U4Random::flat calls.
0134 For example std::raise(SIGINT) could be called when a random draw of
0135 interest is done in order to examine the call stack for that random
0136 consumption before resuming processing.
0137
0138 **/
0139
0140 bool U4Random::isSelect(int photon_idx, int flat_cursor) const
0141 {
0142    if(m_select == nullptr) return false ;
0143    assert( m_select->size() % 2 == 0 );
0144    for(unsigned p=0 ; p < m_select->size()/2 ; p++)
0145    {
0146        int _pidx   = (*m_select)[p*2+0] ;
0147        int _cursor = (*m_select)[p*2+1] ;
0148        bool match = (( _pidx == photon_idx || _pidx == -1)  && (_cursor == flat_cursor || _cursor == -1)) ;
0149        if(match) return true ;
0150    }
0151    return false ;
0152 }
0153
0154 std::string U4Random::descSelect(int photon_idx, int flat_cursor ) const
0155 {
0156     std::stringstream ss ;
0157     ss << "U4Random_select " << ssys::getenvvar("U4Random_select", "-")
0158        << " m_select->size " << ( m_select ? m_select->size() : 0 )
0159        ;
0160
0161     if(m_select)
0162     {
0163         for(unsigned p=0 ; p < m_select->size()/2 ; p++)
0164         {
0165            int _pidx   = (*m_select)[p*2+0] ;
0166            int _cursor = (*m_select)[p*2+1] ;
0167            bool match = ( ( _pidx == photon_idx || _pidx == -1)  && (_cursor == flat_cursor || _cursor == -1) ) ;
0168            ss << " (" << _pidx << "," << _cursor << ") " << ( match ? "YES" : "NO" )  << " " ;
0169         }
0170     }
0171     std::string s = ss.str();
0172     return s ;
0173 }
0174
0175
0176
0177 void U4Random::init()
0178 {
0179     INSTANCE = this ;
0180     m_ready = m_seq != nullptr ;
0181     LOG_IF(error, m_ready == false)
0182         << desc()
0183         << std::endl
0184         << NOTES
0185         << std::endl
0186         ;
0187
0188     LOG(LEVEL) << desc() ;
0189 }
0190
0191 bool U4Random::isReady() const { return m_ready ; }
0192
0193 std::string U4Random::desc() const
0194 {
0195     std::stringstream ss ;
0196     ss
0197          << "U4Random::desc"
0198          << " U4Random::isReady() " << ( m_ready ? "YES" : "NO" )
0199          << " m_seqpath " << ( m_seqpath ? m_seqpath : "-" )
0200          << " m_seq_ni " << m_seq_ni
0201          << " m_seq_nv " << m_seq_nv
0202          ;
0203     return ss.str();
0204 }
0205
0206 std::string U4Random::detail() const
0207 {
0208     std::stringstream ss ;
0209     ss << "U4Random::detail\n"
0210
0211        << std::setw(20) << "m_seqpath_"        << " : " << ( m_seqpath_ ? m_seqpath_ : "-" ) << "\n"
0212        << std::setw(20) << "m_seqpath"         << " : " << ( m_seqpath  ? m_seqpath : "-" ) << "\n"
0213        << std::setw(20) << "m_seqpath_exists"  << " : " << ( m_seqpath_exists  ? "YES" : "NO " ) << "\n"
0214        << std::setw(20) << "m_seq"             << " : " << ( m_seq ? m_seq->desc() : "-" ) << "\n"
0215        << std::setw(20) << "m_seq_index"       << " : " << m_seq_index  << "\n"
0216        << std::setw(20) << "m_cur"             << " : " << ( m_cur ? m_cur->desc() : "-" ) << "\n"
0217        << std::setw(20) << "m_seqmask "        << " : " << ( m_seqmask ? m_seqmask->desc() : "-" ) << "\n"
0218        << std::setw(20) << "desc \n" << desc() << "\n"
0219        ;
0220     std::string s = ss.str();
0221     return s ;
0222 }
0223
0224
0225 /**
0226 U4Random::getNumIndices
0227 -----------------------------
0228
0229 With seqmask running returned the number of seqmask indices otherwise returns the total number of indices.
0230 This corresponds to the total number of available streams of randoms.
0231
0232 **/
0233
0234 size_t U4Random::getNumIndices() const
0235 {
0236    return m_seq && m_seqmask ? m_seqmask_ni : ( m_seq ? m_seq_ni : 0 ) ;
0237 }
0238
0239 /**
0240 U4Random::SetSeed
0241 -----------------------
0242
0243 static control of the seed, NB calling this while enabled will assert
0244 as there is no role for a seed with pre-cooked randoms
0245
0246 TODO: THIS IS RATHER OUT OF PLACE AS NOT MUCH RELATED TO ALIGNED RUNNING, SO RELOCATE TO SEPARATE STRUCT "U4HepRandomEngine::SetSeed" ?
0247
0248 **/
0249
0250 void U4Random::SetSeed(long seed)  // static
0251 {
0252     CLHEP::HepRandomEngine* engine = CLHEP::HepRandom::getTheEngine();
0253     int dummy = 0 ;
0254     engine->setSeed(seed, dummy);
0255 }
0256
0257 /**
0258 U4Random::getMaskedIndex
0259 ------------------------------
0260
0261 When no seqmask is active this just returns the argument.
0262 When a seqmask selection is active indices from the mask are returned.
0263
0264 Masked running allows to reproduce running on subsets of photons from
0265 a larger sample.
0266
0267 **/
0268
0269 size_t U4Random::getMaskedIndex(int index_)
0270 {
0271     if( m_seqmask == nullptr  ) return index_ ;
0272     assert( index_ < m_seqmask_ni );
0273     size_t idx = m_seqmask_values[index_] ;
0274     return idx ;
0275 }
0276
0277
0278 int U4Random::getSequenceIndex() const
0279 {
0280     return m_seq_index ;
0281 }
0282
0283 /**
0284 U4Random::setSequenceIndex
0285 --------------------------------
0286
0287 Switches random stream when index is not negative.
0288 This is used for example to switch between the separate streams
0289 used for each photon.
0290
0291 A negative index disables the control of the Geant4 random engine.
0292
0293 **/
0294
0295 void U4Random::setSequenceIndex(int index_)
0296 {
0297     LOG(LEVEL) << " index " << index_ ;
0298
0299     if( index_ < 0 )
0300     {
0301 #ifndef PRODUCTION
0302 #ifdef DEBUG_TAG
0303         check_cursor_vs_tagslot() ;
0304 #endif
0305 #endif
0306         m_seq_index = index_ ;
0307         disable() ;
0308     }
0309     else
0310     {
0311         size_t idx = getMaskedIndex(index_);
0312         bool idx_in_range = int(idx) < m_seq_ni ;
0313
0314         if(!idx_in_range)
0315             std::cout
0316                 << "FATAL : OUT OF RANGE : "
0317                 << " m_seq_ni " << m_seq_ni
0318                 << " index_ " << index_
0319                 << " idx " << idx << " (must be < m_seq_ni ) "
0320                 << " desc "  << desc()
0321                 ;
0322         assert( idx_in_range );
0323         m_seq_index = idx ;
0324         enable();
0325     }
0326 }
0327
0328
0329 U4Random::~U4Random()
0330 {
0331 }
0332
0333 /**
0334 U4Random::enable
0335 ----------------------
0336
0337 Invokes CLHEP::HepRandom::setTheEngine to *this* U4Random instance
0338 which means that all subsequent calls to G4UniformRand will provide pre-cooked
0339 randoms from the stream controlled by *U4Random::setSequenceIndex*
0340
0341 **/
0342
0343 void U4Random::enable()
0344 {
0345     CLHEP::HepRandom::setTheEngine(this);
0346 }
0347
0348 /**
0349 U4Random::disable
0350 -----------------------
0351
0352 Returns Geant4 to using to the default engine.
0353
0354 **/
0355
0356 void U4Random::disable()
0357 {
0358     CLHEP::HepRandom::setTheEngine(m_default);
0359 }
0360
0361
0362 /**
0363 U4Random::dump
0364 ----------------------
0365
0366 Invokes G4UniformRand *n* times dumping the values.
0367
0368 **/
0369
0370 void U4Random::dump(unsigned n)
0371 {
0372     for(unsigned i=0 ; i < n ; i++)
0373     {
0374         G4double u = G4UniformRand() ;
0375         std::cout
0376             << " i " << std::setw(5) << i
0377             << " u " << std::fixed << std::setw(10) << std::setprecision(5) << u
0378             << std::endl
0379             ;
0380     }
0381 }
0382
0383
0384 /**
0385 U4Random::flat
0386 --------------------
0387
0388 This is the engine method that gets invoked by G4UniformRand calls
0389 and which returns pre-cooked randoms.
0390 The *m_cur_values* cursor is updated to maintain the place in the sequence.
0391
0392 **/
0393
0394 double U4Random::flat()
0395 {
0396     assert(m_seq_index > -1) ;  // must not call when disabled, use G4UniformRand to use standard engine
0397
0398     int cursor = *(m_cur_values + m_seq_index) ;  // get the cursor value to use for this generation, starting from 0
0399
0400     if( cursor >= m_seq_nv )
0401     {
0402         if(m_recycle == false)
0403         {
0404             std::cout
0405                 << "U4Random::flat"
0406                 << " FATAL : not enough precooked randoms and recycle not enabled "
0407                 << " m_seq_index " << m_seq_index
0408                 << " m_seq_nv " << m_seq_nv
0409                 << " cursor " << cursor
0410                 << std::endl
0411                 ;
0412             assert(0);
0413         }
0414         else
0415         {
0416             /*
0417             std::cout
0418                 << "U4Random::flat"
0419                 << " WARNING : not enough precooked randoms are recycling randoms "
0420                 << " m_seq_index " << m_seq_index
0421                 << " m_seq_nv " << m_seq_nv
0422                 << " cursor " << cursor
0423                 << std::endl
0424                 ;
0425             */
0426             cursor = cursor % m_seq_nv ;
0427         }
0428     }
0429
0430
0431
0432     int idx = m_seq_index*m_seq_nv + cursor ;
0433
0434     float  f = m_seq_values[idx] ;
0435     double d = f ;     // promote random float to double
0436     m_flat_prior = d ;
0437
0438     *(m_cur_values + m_seq_index) += 1 ;          // increment the cursor in the array, for the next generation
0439
0440     if(m_seq_index == SEvt::PIDX)
0441     {
0442         std::cout << "-------U4Random::flat" << std::endl << std::endl ;
0443         LOG(info)
0444             << " SEvt::PIDX " << SEvt::PIDX
0445             << " m_seq_index " << std::setw(4) << m_seq_index
0446             << " m_seq_nv " << std::setw(4) << m_seq_nv
0447             << " cursor " << std::setw(4) << cursor
0448             << " idx " << std::setw(4) << idx
0449             << " d " <<  std::setw(10 ) << std::fixed << std::setprecision(5) << d
0450             ;
0451
0452         char* summary = SBacktrace::Summary();
0453
0454         LOG(info) << std::endl << summary ;
0455
0456     }
0457
0458
0459     bool auto_tag = false ;   // unfortunately stack summaries lack vital lines on Linux
0460
0461     if( auto_tag )
0462     {
0463         char* summary = SBacktrace::Summary();
0464         unsigned stack = U4StackAuto::Classify(summary);
0465         bool is_classified = U4StackAuto::IsClassified(stack) ;
0466
0467         //LOG(info) << " stack " << std::setw(2) << stack << " " << U4Stack::Name(stack) ;
0468
0469         LOG_IF(error, !is_classified) << std::endl << summary ;
0470
0471         bool select = isSelect(m_seq_index, cursor ) || is_classified == false ;
0472         if( select )
0473         {
0474             LOG(info) << descSelect(m_seq_index, cursor) ;
0475             switch(m_select_action)
0476             {
0477                 case SDBG::INTERRUPT: std::raise(SIGINT)        ; break ;
0478                 case SDBG::BACKTRACE: SBacktrace::Dump()        ; break ;
0479                 case SDBG::CALLER:    SBacktrace::DumpCaller()  ; break ;
0480                 case SDBG::SUMMARY:   SBacktrace::DumpSummary() ; break ;
0481             }
0482         }
0483
0484 #ifndef PRODUCTION
0485         SEvt::AddTag(1, stack, f );
0486 #endif
0487     }
0488
0489
0490     return d ;
0491 }
0492
0493
0494 #ifndef PRODUCTION
0495 #ifdef DEBUG_TAG
0496 /**
0497 U4Random::check_cursor_vs_tagslot
0498 ----------------------------------
0499
0500 This is called by setSequenceIndex with index -1 signalling the end
0501 of the index. A comparison between the below counts is made:
0502
0503 * number of randoms provided by U4Random::flat for the last m_seq_index as indicated by the cursor
0504 * random consumption tags added with SEvt::AddTag
0505
0506 **/
0507
0508 void U4Random::check_cursor_vs_tagslot()
0509 {
0510     assert(m_seq_index > -1) ;  // must not call when disabled, use G4UniformRand to use standard engine
0511     int cursor = *(m_cur_values + m_seq_index) ;  // get the cursor value to use for this generation, starting from 0
0512
0513     if(SEvt::Exists(1) == false)
0514     {
0515         LOG(error) << " SEvt::Exists(1) false : cannot make this check " ;
0516         return ;
0517     }
0518
0519     int slot = SEvt::GetTagSlot(1);
0520     bool cursor_slot_match = cursor == slot ;
0521
0522     //LOG(info) << " m_seq_index " << m_seq_index << " cursor " << cursor << " slot " << slot << " cursor_slot_match " << cursor_slot_match ;
0523
0524     if(!cursor_slot_match)
0525     {
0526         m_problem_idx.push_back(m_seq_index);
0527         LOG(error)
0528             << " m_seq_index " << m_seq_index
0529             << " cursor " << cursor
0530             << " slot " << slot
0531             << " cursor_slot_match " << cursor_slot_match
0532             << std::endl
0533             << " PROBABLY SOME RANDOM CONSUMPTION LACKS SEvt::AddTag CALLS "
0534             ;
0535     }
0536     assert( cursor_slot_match );
0537 }
0538 #endif
0539 #endif
0540
0541
0542 void U4Random::saveProblemIdx(const char* fold) const
0543 {
0544     std::cout << "U4Random::saveProblemIdx m_problem_idx.size " << m_problem_idx.size() << " (" ;
0545     for(unsigned i=0 ; i <  m_problem_idx.size() ; i++ ) std::cout << m_problem_idx[i] << " " ;
0546     std::cout << ")" << std::endl ;
0547
0548     NP::Write<int>( fold, "problem_idx.npy", m_problem_idx );
0549 }
0550
0551
0552
0553
0554 /**
0555 U4Random::getFlatTag
0556 ----------------------------
0557
0558 By examination of the SBacktrace and Geant4 process state etc..
0559 need to determine what is the random consumer and assign this consumption
0560 of a random number with a standardized tag enumeration from stag.h
0561 to facilitate random alignment with the GPU simulation.
0562
0563 ::
0564
0565    u4t ; U4Random_select=-1,0 U4Random_select_action=backtrace ./U4RecorderTest.sh run
0566        dump full backtraces for first flat consumption of all photons
0567
0568    u4t ; U4Random_select=-1,0 U4Random_select_action=summary ./U4RecorderTest.sh run
0569        dump summary backtraces for first flat consumption of all photons
0570
0571
0572 unsigned U4Random::getFlatTag()
0573 {
0574     char* summary = SBacktrace::Summary();
0575     unsigned stk = U4Stack::Classify(summary);
0576     LOG(info) << " stk " << std::setw(2) << stk << " U4Stack::Name(stk) " << U4Stack::Name(stk) ;
0577     if(!U4Stack::IsClassified(stk)) LOG(error) << std::endl << summary ;
0578     return stk  ; // stk needs translation into the stag.h enumeration
0579 }
0580 **/
0581
0582
0583
0584 int U4Random::getFlatCursor() const
0585 {
0586     if(m_seq_index < 0) return -1 ;
0587     int cursor = *(m_cur_values + m_seq_index) ;  // get the cursor value to use for this generation, starting from 0
0588     return cursor ;
0589 }
0590
0591 double U4Random::getFlatPrior() const
0592 {
0593     return m_flat_prior ;
0594 }
0595
0596
0597 /**
0598 U4Random::flatArray
0599 --------------------------
0600
0601 This method and several others are required as U4Random ISA CLHEP::HepRandomEngine
0602
0603 **/
0604
0605 void U4Random::flatArray(const int size, double* vect)
0606 {
0607      assert(0);
0608 }
0609 void U4Random::setSeed(long seed, int)
0610 {
0611     assert(0);
0612 }
0613 void U4Random::setSeeds(const long * seeds, int)
0614 {
0615     assert(0);
0616 }
0617 void U4Random::saveStatus( const char filename[]) const
0618 {
0619     assert(0);
0620 }
0621 void U4Random::restoreStatus( const char filename[])
0622 {
0623     assert(0);
0624 }
0625 void U4Random::showStatus() const
0626 {
0627     assert(0);
0628 }
0629 std::string U4Random::name() const
0630 {
0631     return NAME ;
0632 }
0633