eic-opticks/qudarap/QCK.cc

0001
0002 #include <iostream>
0003 #include <iomanip>
0004 #include <chrono>
0005 #include <csignal>
0006
0007
0008 #include "NP.hh"
0009 #include "QCK.hh"
0010 #include "SLOG.hh"
0011
0012 template<typename T>
0013 void QCK<T>::save(const char* base, const char* reldir) const
0014 {
0015     rindex->save(base, reldir, "rindex.npy");
0016     bis->save(base, reldir, "bis.npy");
0017     avph->save(base, reldir, "avph.npy");
0018     s2c->save(base, reldir, "s2c.npy");
0019     s2cn->save(base, reldir, "s2cn.npy");
0020     icdf->save(base, reldir, "icdf.npy");
0021     icdf_prop->save(base, reldir, "icdf_prop.npy");
0022 }
0023
0024 template<typename T>
0025 QCK<T>* QCK<T>::Load(const char* base, const char* reldir)   // static
0026 {
0027     NP* rindex = NP::Load(base, reldir, "rindex.npy");
0028     NP* bis = NP::Load(base, reldir, "bis.npy");
0029     NP* avph = NP::Load(base, reldir, "avph.npy");
0030     NP* s2c = NP::Load(base, reldir, "s2c.npy");
0031     NP* s2cn = NP::Load(base, reldir, "s2cn.npy");
0032     NP* icdf = NP::Load(base, reldir, "icdf.npy");
0033     NP* icdf_prop = NP::Load(base, reldir, "icdf_prop.npy");
0034
0035     QCK* qck = nullptr ;
0036     if(rindex == nullptr || bis == nullptr || s2c == nullptr || s2cn == nullptr || avph == nullptr || icdf == nullptr || icdf_prop == nullptr)
0037     {
0038         LOG(fatal)
0039             << " QCK::Load FAILURE "
0040             << " base " << base
0041             << " reldir " << reldir
0042             << std::endl
0043             << " rindex " << rindex
0044             << " bis " << bis
0045             << " s2c " << s2c
0046             << " s2cn " << s2cn
0047             << " icdf " << icdf
0048             << " icdf_prop " << icdf_prop
0049             ;
0050         LOG(error)
0051             << " use QCerenkovIntegralTest to create the ICDF arrays "
0052             ;
0053     }
0054     else
0055     {
0056         qck = new QCK<T> ;
0057         qck->lfold = rindex->lfold ;
0058         qck->rindex = rindex ;
0059         qck->bis = bis ;
0060         qck->avph = avph ;
0061         qck->s2c = s2c ;
0062         qck->s2cn = s2cn ;
0063         qck->icdf = icdf ;
0064         qck->icdf_prop = icdf_prop ;
0065     }
0066     return qck ;
0067 }
0068
0069 template<typename T>
0070 void QCK<T>::init()
0071 {
0072
0073     assert( rindex ) ;
0074     assert( rindex->shape.size() == 2 );
0075     rindex->minmax<T>(emn, emx, 0u );
0076     rindex->minmax<T>(rmn, rmx, 1u );
0077
0078     int w = 10 ;
0079     int p = 5 ;
0080
0081     LOG(info)
0082         << " emn " << std::setw(w) << std::fixed << std::setprecision(p) << emn
0083         << " emx " << std::setw(w) << std::fixed << std::setprecision(p) << emx
0084         << " rmn " << std::setw(w) << std::fixed << std::setprecision(p) << rmn
0085         << " rmx " << std::setw(w) << std::fixed << std::setprecision(p) << rmx
0086         ;
0087
0088
0089 }
0090
0091
0092 template<typename T> int QCK<T>::find_bisIdx( const T BetaInverse ) const
0093 {
0094     bool in_range ;
0095     unsigned column = 0u ; // placeholder, as bis is 1d
0096     int ibin = bis->pfindbin<T>(BetaInverse, column, in_range );
0097     assert( in_range == true );
0098     int item = ibin - 1 ;  // ibin is 1-based for in_range bins, so item is 0-based
0099     assert( item > -1 );
0100
0101 /*
0102     int w = 10 ;
0103     int p = 6 ;
0104     std::cout
0105         << "QCK::find_bisIdx"
0106         << " BetaInverse "      << std::fixed << std::setw(w) << std::setprecision(p) << BetaInverse
0107         << " in_range "         << std::setw(5) << in_range
0108         << " ibin "             << std::setw(5) << ibin
0109         << " item "             << std::setw(5) << item
0110         << std::endl
0111         ;
0112 */
0113
0114     return item ;
0115 }
0116
0117
0118
0119 /*
0120 QCK::energy_lookup
0121 ----------------------
0122
0123 For BetaInverse such as 1.0 which is less than RINDEX_min
0124 Cerenkov can happen anywhere across the full energy range.
0125
0126 For BetaInverse approaching the limits of where Cerenkov
0127 can happen, like when there is only a small RINDEX island
0128 peaking above the BetaInverse sea the ICDF will be shaped
0129 more and more like a step function, so no matter the
0130 input random number will get close to the same energy.
0131
0132 For BetaInverse that is greater than RINDEX_max the "island" has
0133 been submerged and Cerenkov does not happen, **so this should not be called**
0134 as it will trip an assert.
0135
0136 Avoid tripping the assert by only calling QCK::energy_looking with BetaInverse
0137 that returns true from QCK::is_permissable.
0138
0139 see ~/np/tests/NPget_edgesTest.cc
0140
0141 NB this is currently using NP::pdomain binary bin search and lookups on the CDF rather
0142 than NP::interp/NP::interpHD direct readoffs from the ICDF.
0143 However for usage on GPU the ICDF is typically required for the GPU texture array,
0144 so forming the ICDF is normally necessary.
0145
0146 Actually a GPU equivalent of NP::pdomain (qprop?) could also be used, but
0147 the performance of that is expected to be dramatically less than texture lookups,
0148 although this needs to be measured.
0149
0150 */
0151
0152 template<typename T> T QCK<T>::energy_lookup_( const T BetaInverse, const T u, double& dt, bool icdf_ ) const
0153 {
0154     typedef std::chrono::high_resolution_clock::time_point TP ;
0155     TP t0 = std::chrono::high_resolution_clock::now() ;
0156
0157     int ibis = find_bisIdx(BetaInverse);
0158     T en ;
0159
0160     if(icdf_)
0161     {
0162         // *NP::interpHD* asserts with GPU ready icdf of shape (1000, 1000, 4)
0163         // need to use icdf_prop of shape eg (1000, 1000, 4, 2 ) with the interleaved domain
0164         // for CPU testing
0165
0166         unsigned hd_factor = icdf_prop->get_meta<unsigned>("hd_factor", 0u);
0167         en = icdf_prop->interpHD<T>(u, hd_factor, ibis );
0168     }
0169     else
0170     {
0171         bool dump = false ;
0172         en = s2cn->pdomain<T>(u, ibis, dump );
0173     }
0174
0175
0176     TP t1 = std::chrono::high_resolution_clock::now() ;
0177     std::chrono::duration<double> t10 = t1 - t0;
0178     dt = t10.count()*DT_SCALE ;
0179
0180     return en  ;
0181 }
0182
0183
0184
0185 /**
0186 QCK::energy_range_s2cn
0187 ------------------------
0188
0189 For _SplitBin the energy range obtained from s2cn is the full energy range.
0190
0191 **/
0192
0193 template<typename T> void QCK<T>::energy_range_s2cn( T& emin, T& emax, const T BetaInverse, bool dump ) const
0194 {
0195     int ibis = find_bisIdx(BetaInverse);
0196     unsigned k = 0 ; // 1st payload slot is domain, the energy
0197     emin = s2cn->get<T>(ibis,  0, k) ;   // first energy
0198     emax = s2cn->get<T>(ibis, -1, k) ;   // last energy
0199
0200     T edom = emax - emin ;
0201     if(dump) std::cout
0202         << "QCK::energy_range_s2cn"
0203         << " BetaInverse " << std::setw(10) << std::fixed << std::setprecision(4) << BetaInverse
0204         << " ibis " << std::setw(10) << ibis
0205         << " bis[ibis] " <<  std::setw(10) << std::fixed << std::setprecision(4) << bis->get<T>(ibis)
0206         << " emin " << std::setw(10) << std::fixed << std::setprecision(4) << emin
0207         << " emax " << std::setw(10) << std::fixed << std::setprecision(4) << emax
0208         << " edom " << std::setw(10) << std::fixed << std::setprecision(4) << edom
0209         << std::endl
0210         ;
0211 }
0212
0213 template<typename T> void QCK<T>::energy_range_avph( T& emin, T& emax, const T BetaInverse, bool dump ) const
0214 {
0215     int ibis = find_bisIdx(BetaInverse);
0216     T BetaInverse2 = avph->get<T>(ibis,0) ;
0217     emin = avph->get<T>(ibis,1) ;
0218     emax = avph->get<T>(ibis,2) ;
0219     T averageNumPhotons = avph->get<T>(ibis,3) ;
0220
0221     T edom = emax - emin ;
0222     if(dump) std::cout
0223         << "QCK::energy_range_avph"
0224         << " BetaInverse " << std::setw(10) << std::fixed << std::setprecision(4) << BetaInverse
0225         << " BetaInverse2 " << std::setw(10) << std::fixed << std::setprecision(4) << BetaInverse2
0226         << " ibis " << std::setw(10) << ibis
0227         << " bis[ibis] " <<  std::setw(10) << std::fixed << std::setprecision(4) << bis->get<T>(ibis)
0228         << " emin " << std::setw(10) << std::fixed << std::setprecision(4) << emin
0229         << " emax " << std::setw(10) << std::fixed << std::setprecision(4) << emax
0230         << " edom " << std::setw(10) << std::fixed << std::setprecision(4) << edom
0231         << " averageNumPhotons " << std::setw(10) << std::fixed << std::setprecision(4) << averageNumPhotons
0232         << std::endl
0233         ;
0234 }
0235
0236
0237
0238 template<typename T> const double QCK<T>::DT_SCALE = 1e6 ;
0239
0240 /**
0241 QCK::energy_sample_
0242 ----------------------
0243
0244 *u1_s2max* needs to be less than s2 for the energy to be accepted
0245
0246 Using full emn to emx range is problematic for BetaInverse close to the
0247 RINDEX peak as most of the randoms will be wasted because the vast majority
0248 of the range is not Cerenkov permitted.
0249
0250
0251 **/
0252
0253 template<typename T> T QCK<T>::energy_sample_( const T emin, const T emax, const T BetaInverse, const std::function<T()>& rng, double& dt, unsigned& count ) const
0254 {
0255     typedef std::chrono::high_resolution_clock::time_point TP ;
0256     TP t0 = std::chrono::high_resolution_clock::now() ;
0257
0258     const T one(1.) ;
0259     T u0, u1, en, ri, ct, s2, u1_s2max ;
0260
0261     const T ctmax = BetaInverse/rmx ;
0262     const T s2max = ( one - ctmax )*( one + ctmax ) ;
0263     count = 0 ;
0264
0265     do {
0266         u0 = rng() ;
0267         u1 = rng() ;
0268
0269         en = emin + u0*(emax-emin) ;
0270         ri = rindex->interp<T>(en) ;
0271         ct = BetaInverse/ri ;
0272         s2 = ( one - ct )*( one + ct ) ;
0273         u1_s2max = u1*s2max ;
0274
0275         count += 1 ;
0276
0277     } while ( u1_s2max > s2 ) ;
0278
0279
0280     TP t1 = std::chrono::high_resolution_clock::now() ;
0281     std::chrono::duration<double> t10 = t1 - t0;
0282     dt = t10.count()*DT_SCALE ;
0283
0284     return en ;
0285 }
0286
0287
0288 template<typename T> NP* QCK<T>::energy_lookup( const T BetaInverse, const NP* uu, NP* tt, bool icdf_ ) const
0289 {
0290     unsigned ndim = uu->shape.size() ;
0291     bool ndim_expect = ndim == 1  ;
0292     assert( ndim_expect );
0293     if(!ndim_expect) std::raise(SIGINT);
0294
0295     unsigned ni = uu->shape[0] ;
0296     const T* uu_v = uu->cvalues<T>();
0297
0298
0299     if(tt) assert( tt->shape.size() == 1u && unsigned(tt->shape[0]) == ni );
0300     T* tt_v = tt ? tt->values<T>() : nullptr ;
0301
0302
0303     NP* en = NP::MakeLike(uu);
0304     T* en_v = en->values<T>();
0305     double dt ;
0306
0307     for(unsigned i=0 ; i < ni ; i++)
0308     {
0309         en_v[i] = energy_lookup_( BetaInverse, uu_v[i],  dt, icdf_ ) ;
0310         if(tt_v) tt_v[i] = dt ;
0311     }
0312
0313     en->set_meta<std::string>("qck_lfold", lfold );
0314
0315     return en ;
0316 }
0317
0318 /**
0319 QCK::energy_sample
0320 -------------------
0321
0322 Returns array of energy samples of shape (ni,) and if argument tt is non-null and of shape (ni,) sets cputimes
0323
0324 **/
0325
0326 template<typename T> NP* QCK<T>::energy_sample( const T BetaInverse, const std::function<T()>& rng, unsigned ni, NP* tt ) const
0327 {
0328     NP* en = NP::Make<T>(ni);
0329     T* en_v = en->values<T>();
0330
0331     if(tt) assert( tt->shape.size() == 1u && unsigned(tt->shape[0]) == ni );
0332     T* tt_v = tt ? tt->values<T>() : nullptr ;
0333
0334     double dt ;
0335     unsigned count_max(0);
0336     unsigned count(0) ;
0337
0338     T emin ;
0339     T emax ;
0340     bool dump = false ;
0341     energy_range_avph( emin, emax, BetaInverse, dump);
0342
0343
0344     for(unsigned i=0 ; i < ni ; i++)
0345     {
0346         en_v[i] = energy_sample_( emin, emax, BetaInverse, rng, dt, count  ) ;
0347         if(tt_v) tt_v[i] = dt ;
0348
0349         if( count > count_max ) count_max = count ;
0350
0351         if(count > 1000)
0352             std::cout
0353                 << "QCK::energy_sample "
0354                 << " i " << std::setw(7) << i
0355                 << " count " << std::setw(7) << count
0356                 << std::endl
0357                 ;
0358
0359     }
0360
0361     LOG(info) << " count_max " << count_max ;
0362     en->set_meta<std::string>("qck_lfold", lfold );
0363
0364     return en ;
0365 }
0366
0367
0368
0369 /**
0370 QCK::is_permissable
0371 ----------------------
0372
0373 Look for the BetaInverse value within the bis array, returning true when in range.
0374
0375 **/
0376
0377 template<typename T> bool QCK<T>::is_permissable( const T BetaInverse) const
0378 {
0379     assert( bis->shape.size() == 1 );
0380
0381     bool in_range ;
0382     unsigned column = 0u ; // placeholder, as bis is 1d
0383     int ibin = bis->pfindbin<T>(BetaInverse, column, in_range );
0384
0385     T lo, hi ;
0386     bis->get_edges<T>(lo,hi,column,ibin);    // when not in_range the  hi and lo are edge values
0387
0388     int w = 10 ;
0389     int p = 6 ;
0390
0391     std::cout
0392         << "QCK::is_permissable "
0393         << " BetaInverse "      << std::fixed << std::setw(w) << std::setprecision(p) << BetaInverse
0394         << " bis.ibin "         << std::setw(5) << ibin
0395         << " bis.in_range "     << std::setw(5) << in_range
0396         << " bis.get_edges ["
0397         << " " << std::fixed << std::setw(w) << std::setprecision(p) << lo
0398         << " " << std::fixed << std::setw(w) << std::setprecision(p) << hi
0399         << "]"
0400         << std::endl
0401         ;
0402
0403     assert( BetaInverse >= lo && BetaInverse <= hi );
0404     return in_range ;
0405 }
0406
0407
0408 template struct QCK<float> ;
0409 template struct QCK<double> ;
0410