qudarap/tests/QCKTest.cc

0001 /**
0002 QCKTest.cc
0003 =============
0004
0005 Loads ICDF created and persisted by QCerenkovIntegralTest
0006 and uses it to create Cerenkov energy samples for a variety of BetaInverse
0007 by both lookup sampling and for comparison by the traditional rejection sampling approach.
0008
0009 **/
0010
0011 #include <iostream>
0012 #include <iomanip>
0013 #include <random>
0014
0015 #include "SStr.hh"
0016 #include "SRng.hh"
0017 #include "SPath.hh"
0018 #include "QCK.hh"
0019 #include "NP.hh"
0020 #include "OPTICKS_LOG.hh"
0021
0022 void test_energy_lookup_one( const QCK<double>* qck, double BetaInverse, double u  )
0023 {
0024     double dt ;
0025     bool icdf_ = true ;
0026     double en = qck->energy_lookup_( BetaInverse, u, dt, icdf_ );
0027     int p = 7 ;
0028     int w = p + 10 ;
0029     LOG(info)
0030         << " BetaInverse " << std::fixed << std::setw(w) << std::setprecision(p) << BetaInverse
0031         << " u " << std::fixed << std::setw(w) << std::setprecision(p) << u
0032         << " en " << std::fixed << std::setw(w) << std::setprecision(p) << en
0033         << " dt " << std::fixed << std::setw(w) << std::setprecision(p) << dt
0034         ;
0035 }
0036
0037 void test_energy_lookup_many( const QCK<double>* qck, double BetaInverse, unsigned ni, const char* base, bool icdf_ )
0038 {
0039     bool biok = qck->is_permissable(BetaInverse);
0040     if(biok == false)
0041     {
0042         LOG(fatal)
0043             << " BetaInverse not permitted as no photons " << std::fixed << std::setw(10) << std::setprecision(4) << BetaInverse
0044             ;
0045         return ;
0046     }
0047
0048     NP* tt = NP::Make<double>( ni );
0049     NP* uu = NP::MakeUniform<double>( ni ) ;
0050
0051     NP* en = qck->energy_lookup(BetaInverse, uu, tt, icdf_ ) ;
0052     const char* path = SPath::MakePath<double>(base, nullptr, BetaInverse, icdf_ ? "test_energy_lookup_many_icdf.npy" :  "test_energy_lookup_many_s2cn.npy" );
0053     en->save(path);
0054
0055     const char* tt_path = SPath::MakePath<double>(base, nullptr, BetaInverse, "test_energy_lookup_many_tt.npy" );
0056     tt->save(tt_path);
0057 }
0058
0059
0060 void test_energy_sample_one( const QCK<double>* qck, double BetaInverse )
0061 {
0062     double dt ;
0063     unsigned seed = 0u ;
0064     SRng<double> rng(seed) ;
0065
0066     double emin = 1.55 ;
0067     double emax = 15.5 ;
0068     unsigned count ;
0069     double en = qck->energy_sample_( emin, emax, BetaInverse, rng, dt, count );
0070
0071     int p = 7 ;
0072     int w = p + 10 ;
0073
0074     LOG(info)
0075         << " BetaInverse " << std::fixed << std::setw(w) << std::setprecision(p) << BetaInverse
0076         << " en " << std::fixed << std::setw(w) << std::setprecision(p) << en
0077         << " dt " << std::fixed << std::setw(w) << std::setprecision(p) << dt
0078         << " count " << std::setw(7) << count
0079         ;
0080 }
0081
0082 void test_energy_sample_many( const QCK<double>* qck, double BetaInverse, unsigned ni, const char* base )
0083 {
0084     bool biok = qck->is_permissable(BetaInverse);
0085     if(biok == false)
0086     {
0087         LOG(fatal)
0088             << " BetaInverse not permitted as no photons " << std::fixed << std::setw(10) << std::setprecision(4) << BetaInverse
0089             ;
0090         return ;
0091     }
0092
0093     unsigned seed = 0u ;
0094     SRng<double> rng(seed) ;
0095
0096     NP* tt = NP::Make<double>( ni );
0097     NP* en = qck->energy_sample( BetaInverse, rng, ni, tt ) ;
0098     const char* path = SPath::MakePath<double>( base, nullptr, BetaInverse, "test_energy_sample_many.npy" );
0099     en->save(path) ;
0100
0101     const char* tt_path = SPath::MakePath<double>( base, nullptr, BetaInverse, "test_energy_sample_many_tt.npy" );
0102     tt->save(tt_path) ;
0103 }
0104
0105
0106 int main(int argc, char** argv)
0107 {
0108     OPTICKS_LOG(argc, argv);
0109
0110     char d = 'S' ;
0111     char t = argc > 1 ? argv[1][0] : d ;
0112
0113     const char* reldir = nullptr ;
0114     switch(t)
0115     {
0116         case 'S': reldir = "test_makeICDF_SplitBin" ; break ;
0117         case 'U': reldir = "test_makeICDF_UpperCut" ; break ;
0118     }
0119     assert( reldir && "require S or U argument to pick ICDF reldir" );
0120     const char* icdf_base = SPath::Resolve("$TMP/QCerenkovIntegralTest", false);
0121     LOG(info) << " t " << t << " icdf_base " << icdf_base << " reldir " << reldir ;
0122
0123     // hmm: save samples into ICDF subdir to explicitly record the vital connection
0124
0125     QCK<double>* qck = QCK<double>::Load(icdf_base, reldir);
0126     if( qck == nullptr ) return 0 ;
0127
0128     qck->init();
0129     const char* lfold = qck->lfold.c_str() ;
0130
0131     LOG(info) << " qck.lfold " << lfold ;
0132     LOG(info) << " qck.bis  " << qck->bis->desc();
0133     LOG(info) << " qck.bis->meta  " << qck->bis->meta;
0134
0135     LOG(info) << " qck.s2c  " << qck->s2c->desc();
0136     LOG(info) << " qck.s2cn " << qck->s2cn->desc();
0137
0138     int create_dirs = 2 ;  //2:dirpath
0139     const char* sample_base = SPath::Resolve(lfold , "QCKTest", create_dirs ) ;
0140     LOG(info) << " sample_base " << sample_base ;
0141
0142     std::vector<double> vbis ;
0143     for( double bi=1.0 ; bi < qck->rmx ; bi+=0.05 ) vbis.push_back(bi);
0144     vbis.push_back(1.792);  // extreme peak : some tiny fraction of a photon
0145     //bis.push_back(1.45);   // pdomain assert, from going slightly non-monotonic for _UpperCut not _SplitBin
0146
0147     NP* bis = NP::Make<double>(vbis);
0148
0149     std::stringstream ss ;
0150     ss << "qck.lfold:" << lfold
0151        << std::endl
0152        << "qck.bis.meta:" << qck->bis->meta
0153        << std::endl
0154        ;
0155     bis->meta = ss.str();   // metadata about this test, used for annotations from python
0156
0157
0158     int nbis =  bis->shape[0] ;
0159     unsigned num_gen = 1000000 ;
0160
0161     LOG(info)
0162         << " num_gen " << num_gen
0163         << " nbis " << nbis
0164         ;
0165
0166
0167     double emin, emax ;
0168     for(int i=0 ; i < nbis ; i++)
0169     {
0170         double BetaInverse = bis->get<double>(i) ;
0171         //qck->energy_range_s2cn( emin, emax, BetaInverse, true );
0172         qck->energy_range_avph( emin, emax, BetaInverse, true );
0173     }
0174
0175
0176     for(int i=0 ; i < nbis ; i++)
0177     {
0178         double BetaInverse = bis->get<double>(i) ;
0179         LOG(info) << " lookup BetaInverse icdf:1 " << std::fixed << std::setw(10) << std::setprecision(4) << BetaInverse ;
0180         test_energy_lookup_many( qck, BetaInverse, num_gen, sample_base, true );
0181     }
0182
0183     for(int i=0 ; i < nbis ; i++)
0184     {
0185         double BetaInverse = bis->get<double>(i) ;
0186         LOG(info) << " lookup BetaInverse icdf:0 " << std::fixed << std::setw(10) << std::setprecision(4) << BetaInverse ;
0187         test_energy_lookup_many( qck, BetaInverse, num_gen, sample_base, false );
0188     }
0189
0190     for(int i=0 ; i < nbis ; i++)
0191     {
0192         double BetaInverse = bis->get<double>(i) ;
0193         LOG(info) << " sampling BetaInverse " << std::fixed << std::setw(10) << std::setprecision(4) << BetaInverse ;
0194         test_energy_sample_many( qck, BetaInverse, num_gen, sample_base );
0195     }
0196
0197     bis->save(sample_base, "bis.npy" );
0198
0199     return 0 ;
0200 }
0201