eic-opticks/sysrap/s_seq.h

0001 #pragma once
0002 /**
0003 s_seq.h : Loading and serving up precooked randoms, usaully GPU generated by curand
0004 =====================================================================================
0005
0006 +-----------------------+------------------------------------------------+
0007 |  Two defaults         |   Notes                                        |
0008 +=======================+================================================+
0009 |   DEFAULT_SEQPATH_H1  |  1 file with (100k,16,16) precooked randoms    |
0010 +-----------------------+------------------------------------------------+
0011 |   DEFAULT_SEQPATH_M1  | 10 files with (100k,16,16) precooked randoms   |
0012 +-----------------------+------------------------------------------------+
0013
0014 As loading the (M1,16,16) randoms from 10 files and concatenating them
0015 takes a few too many seconds the defualt is the first 100k option.
0016
0017 Switch to the larger 1M option with::
0018
0019    export s_seq__SeqPath_DEFAULT_LARGE=1
0020
0021 Or alternatively define a path to your own precooked random files::
0022
0023    export OPTICKS_RANDOM_SEQPATH=...
0024
0025 **/
0026
0027 #include "NP.hh"
0028 #include "ssys.h"
0029
0030 struct s_seq
0031 {
0032     static constexpr const char* OPTICKS_RANDOM_SEQPATH = "OPTICKS_RANDOM_SEQPATH" ;
0033     static constexpr const char* EKEY = "s_seq__SeqPath_DEFAULT_LARGE" ;
0034     static constexpr const char* DEFAULT_SEQPATH_H1 =
0035     "$HOME/.opticks/precooked/QSimTest/rng_sequence/rng_sequence_f_ni1000000_nj16_nk16_tranche100000/rng_sequence_f_ni100000_nj16_nk16_ioffset000000.npy" ;
0036     static constexpr const char* DEFAULT_SEQPATH_M1 =
0037     "$HOME/.opticks/precooked/QSimTest/rng_sequence/rng_sequence_f_ni1000000_nj16_nk16_tranche100000" ;
0038
0039     static const char* SeqPath() ;
0040
0041     s_seq();
0042
0043     std::string desc() const  ;
0044     double flat();
0045     std::string demo(int n) ;
0046     int  getSequenceIndex() const ;
0047     void setSequenceIndex(int index_);   // -ve to disable, must be less than ni
0048     bool is_enabled() const ;
0049
0050 private:
0051     const char*   m_seqpath ;
0052     const NP*     m_seq ;
0053     const float*  m_seq_values ;
0054     int           m_seq_ni ;
0055     int           m_seq_nv ;
0056     int           m_seq_index ;
0057     int           m_pidx ;
0058     NP*           m_cur ;
0059     int*          m_cur_values ;
0060     bool          m_recycle ;
0061     double        m_flat_prior ;
0062 };
0063
0064 inline const char* s_seq::SeqPath() // static
0065 {
0066     bool DEFAULT_LARGE = ssys::getenvbool(EKEY);
0067     const char* default_seqpath = DEFAULT_LARGE ? DEFAULT_SEQPATH_M1 : DEFAULT_SEQPATH_H1 ;
0068     return ssys::getenvvar(OPTICKS_RANDOM_SEQPATH, default_seqpath ) ;
0069 }
0070
0071 inline s_seq::s_seq()
0072     :
0073     m_seqpath(U::Resolve(SeqPath())),
0074     m_seq(m_seqpath ? NP::LoadIfExists(m_seqpath) : nullptr),
0075     m_seq_values(m_seq ? m_seq->cvalues<float>() : nullptr ),
0076     m_seq_ni(m_seq ? m_seq->shape[0] : 0 ),                        // num items
0077     m_seq_nv(m_seq ? m_seq->shape[1]*m_seq->shape[2] : 0 ),        // num values in each item
0078     m_seq_index(-1),
0079     m_pidx(ssys::getenvint("PIDX",-100)),
0080     m_cur(NP::Make<int>(m_seq_ni)),
0081     m_cur_values(m_cur->values<int>()),
0082     m_recycle(true)
0083 {
0084 }
0085
0086 inline std::string s_seq::desc() const
0087 {
0088     std::stringstream ss ;
0089     ss << "s_seq::desc"
0090        << std::endl
0091        << " m_seqpath " << ( m_seqpath ? m_seqpath : "-" )
0092        << std::endl
0093        << " m_seq " << ( m_seq ? m_seq->sstr() : "-" )
0094        << std::endl
0095        ;
0096     std::string str = ss.str();
0097     return str ;
0098 }
0099
0100 inline double s_seq::flat()
0101 {
0102     assert(m_seq_index > -1) ;  // must not call when disabled, use G4UniformRand to use standard engine
0103     int cursor = *(m_cur_values + m_seq_index) ;  // get the cursor value to use for this generation, starting from 0
0104
0105     if( cursor >= m_seq_nv )
0106     {
0107         if(m_recycle == false)
0108         {
0109             std::cout
0110                 << "s_seq::flat"
0111                 << " FATAL : not enough precooked randoms and recycle not enabled "
0112                 << " m_seq_index " << m_seq_index
0113                 << " m_seq_nv " << m_seq_nv
0114                 << " cursor " << cursor
0115                 << std::endl
0116                 ;
0117             assert(0);
0118         }
0119         else
0120         {
0121             std::cout
0122                 << "s_seq::flat"
0123                 << " WARNING : not enough precooked randoms are recycling randoms "
0124                 << " m_seq_index " << m_seq_index
0125                 << " m_seq_nv " << m_seq_nv
0126                 << " cursor " << cursor
0127                 << std::endl
0128                 ;
0129             cursor = cursor % m_seq_nv ;
0130         }
0131     }
0132
0133
0134     int idx = m_seq_index*m_seq_nv + cursor ;
0135
0136     float  f = m_seq_values[idx] ;
0137     double d = f ;               // promote random float to double
0138     m_flat_prior = d ;
0139
0140 #ifdef MOCK_CUDA_DEBUG
0141     if(m_seq_index == m_pidx) printf("//s_seq::flat.MOCK_CUDA_DEBUG m_seq_index %5d cursor %4d d %10.5f \n", m_seq_index, cursor, d ) ;
0142 #endif
0143
0144     *(m_cur_values + m_seq_index) += 1 ;   // increment the cursor in the array, for the next generation
0145
0146     return d ;
0147 }
0148
0149 inline std::string s_seq::demo(int n)
0150 {
0151     std::stringstream ss ;
0152     ss << "s_seq::demo m_seq_index " << m_seq_index << std::endl ;
0153     for(int i=0 ; i < n ; i++) ss
0154          << std::setw(4) << i
0155          << " : "
0156          << std::fixed << std::setw(10) << std::setprecision(5) << flat()
0157          << std::endl
0158          ;
0159     std::string str = ss.str();
0160     return str ;
0161 }
0162
0163 inline int  s_seq::getSequenceIndex() const
0164 {
0165     return m_seq_index ;
0166 }
0167 inline void s_seq::setSequenceIndex(int index_)
0168 {
0169     if( index_ < 0 )
0170     {
0171         m_seq_index = index_ ;
0172     }
0173     else
0174     {
0175         int idx = index_ ; // ASSUME NO MASKS
0176         bool idx_in_range = int(idx) < m_seq_ni ;
0177         if(!idx_in_range)
0178             std::cout
0179                 << "s_seq::setSequenceIndex"
0180                 << "FATAL : OUT OF RANGE : "
0181                 << " m_seq_ni " << m_seq_ni
0182                 << " index_ " << index_
0183                 << " idx " << idx << " (must be < m_seq_ni ) "
0184                 << " desc "  << desc()
0185                 ;
0186         assert( idx_in_range );
0187         m_seq_index = idx ;
0188     }
0189 }
0190 inline bool s_seq::is_enabled() const
0191 {
0192     return m_seq_index > -1 ;
0193 }
0194
0195