qudarap/tests/QSim_MockTest.cc

0001 /**
0002 QSim_MockTest.cc : CPU tests of QSim.hh/qsim.h CUDA code using MOCK_CURAND mocking
0003 =======================================================================================
0004
0005 Testing GPU code on CPU requires mocking of CUDA API including:
0006
0007 1. tex2D lookups
0008 2. curand random generation
0009 3. erfcinvf : inverse complementary error function
0010
0011 There are now lots of examples of curand mocking,
0012 search for MOCK_CURAND, MOCK_CUDA. See::
0013
0014     sysrap/srngcpu.h
0015     sysrap/scurand.h
0016
0017 Mocking tex2D lookups is not so common. See::
0018
0019     sysrap/s_mock_texture.h
0020     sysrap/stexture.h
0021
0022 and search for MOCK_TEXTURE, MOCK_CUDA.
0023
0024
0025 HMM: QSim.hh not very amenable to standalone use because the boatload of
0026 headers that come with it : so operating at lower level.
0027
0028 Standalone compile and run with::
0029
0030    ~/o/qudarap/tests/QSim_MockTest.sh
0031
0032 **/
0033
0034 #include "NPFold.h"
0035
0036 #include "ssys.h"
0037 #include "spath.h"
0038 #include "scuda.h"
0039 #include "smath.h"    // includes s_mock_erfinvf.h when MOCK_CUDA is defined
0040 #include "squad.h"
0041 #include "srec.h"
0042 #include "stag.h"
0043 #include "sflow.h"
0044 #include "sphoton.h"
0045 #include "sstate.h"
0046 #include "scerenkov.h"
0047
0048 #include "srngcpu.h"
0049 using RNG = srngcpu ;
0050
0051
0052 #include "stexture.h"   // includes s_mock_texture.h when MOCK_TEXTURE OR MOCK_CUDA defined
0053
0054 #include "SPMT.h"
0055 #include "SBnd.h"
0056 #include "OpticksPhoton.hh"
0057
0058 #include "QBase.hh"
0059 #include "QPMT.hh"
0060 #include "QBnd.hh"
0061 #include "QOptical.hh"
0062
0063 #include "qpmt.h"
0064 #include "qbnd.h"
0065 #include "qsim.h"
0066
0067 struct QSim_MockTest
0068 {
0069     static const char* TEST ;
0070
0071     static constexpr const char* BASE = "$HOME/.opticks/GEOM/$GEOM/CSGFoundry/SSim/stree/standard" ;
0072     static constexpr const char* BND =
0073     "Pyrex/HamamatsuR12860_PMT_20inch_photocathode_mirror_logsurf/HamamatsuR12860_PMT_20inch_photocathode_mirror_logsurf/Vacuum" ;
0074
0075     const char* FOLD ;
0076     const char* CHECK ;
0077     RNG         rng ;
0078
0079     const NP* optical ;
0080     const NP* bnd ;
0081
0082     const QBase*    q_base ;
0083     const QOptical* q_optical ;
0084     const QBnd*     q_bnd ;
0085
0086     const SBnd*     s_bnd ;
0087     int   boundary ;
0088
0089     const NPFold* jpmt ;
0090     const QPMT<float>* q_pmt ;
0091     qsim* sim ;
0092     const int num ;
0093     const int PIDX ;
0094     NP*  a ;
0095
0096     QSim_MockTest();
0097
0098     void init();
0099     void init_bnd();
0100     std::string desc() const;
0101
0102     int generate_photon_dummy();
0103     int uniform_sphere();
0104
0105     int propagate_at_boundary_manual();
0106
0107     void setup_prd(    quad2& prd );
0108     void setup_photon( sphoton& p );
0109
0110 #ifdef WITH_CUSTOM4
0111     int propagate_at_surface_CustomART_manual();
0112 #endif
0113
0114     int fill_state();
0115     int propagate_at_boundary();
0116     int propagate();
0117     void SmearNormal(int chk, float value);
0118     int SmearNormal_SigmaAlpha_one();
0119     int run();
0120
0121     int main();
0122 };
0123
0124 const char* QSim_MockTest::TEST = ssys::getenvvar("TEST", "ALL") ;
0125
0126 inline QSim_MockTest::QSim_MockTest()
0127     :
0128     FOLD(ssys::getenvvar("FOLD")),
0129     CHECK(spath::Basename(FOLD)),
0130     rng(),
0131     optical(NP::Load(BASE, "optical.npy")),
0132     bnd(    NP::Load(BASE, "bnd.npy")),
0133     q_base( new QBase ),
0134     q_optical(optical ? new QOptical(optical) : nullptr),
0135     q_bnd(    bnd     ? new QBnd(bnd)         : nullptr),
0136     s_bnd(    bnd     ? new SBnd(bnd)         : nullptr),
0137     boundary(s_bnd ? s_bnd->getBoundaryIndex(BND) : -1),
0138     jpmt(SPMT::Serialize()),
0139     q_pmt( jpmt ? new QPMT<float>( jpmt ) : nullptr),
0140     sim(new qsim),
0141     num(ssys::getenvint("NUM",1000)),
0142     PIDX(ssys::getenvint("PIDX",-100)),
0143     a(nullptr)
0144 {
0145     init();
0146 }
0147
0148
0149 inline void QSim_MockTest::init()
0150 {
0151     rng.seed = 1 ;
0152     init_bnd();
0153
0154     assert( q_pmt );
0155     //rng.set_fake(0.); // 0/1:forces transmit/reflect
0156
0157     sim->base = q_base ? q_base->d_base : nullptr ;
0158     sim->pmt = q_pmt->d_pmt ;
0159 }
0160
0161
0162 inline void QSim_MockTest::init_bnd()
0163 {
0164     assert( bnd );
0165     bool have_boundary = boundary > -1 ;
0166     if(!have_boundary) std::cerr
0167         << "QSim_MockTest::init_bnd"
0168         << std::endl
0169         << " FATAL FAILED TO LOOKUP boundary " << boundary
0170         << std::endl
0171         << " BND " << BND
0172         << std::endl
0173         << " NO BOUNDARY WITH THAT SPEC PRESENT WITHIN LOADED GEOMETRY "
0174         << std::endl
0175         << " BASE " << BASE
0176         << std::endl
0177         << " PROBABLY THE GEOM ENVVAR IS MIS-CONFIGURED : CHANGE IT USING GEOM BASH FUNCTION "
0178         << std::endl
0179         ;
0180     assert( have_boundary );
0181     sim->bnd = q_bnd->d_qb ;
0182 }
0183
0184
0185 inline std::string QSim_MockTest::desc() const
0186 {
0187     std::stringstream ss ;
0188     ss << "QSim_MockTest::desc" << std::endl
0189        << " bnd " << ( bnd ? bnd->sstr() : "-" )
0190        << " optical " << ( optical ? optical->sstr() : "-" )
0191        << " boundary " << boundary << std::endl
0192        << " s_bnd.getBoundarySpec " << ( s_bnd ? s_bnd->getBoundarySpec(boundary) : "-" )
0193        << std::endl
0194        ;
0195     std::string str = ss.str();
0196     return str ;
0197 }
0198
0199 inline int QSim_MockTest::generate_photon_dummy()
0200 {
0201     sphoton p ;
0202     p.zero();
0203     quad6 gs ;
0204     gs.zero();
0205     unsigned photon_id = 0 ;
0206     unsigned genstep_id = 0 ;
0207     sim->generate_photon_dummy(p, rng, gs, photon_id, genstep_id );
0208     return 0 ;
0209 }
0210
0211 inline int QSim_MockTest::uniform_sphere()
0212 {
0213     for(int i=0 ; i < 10 ; i++)
0214     {
0215         float3 dir = sim->uniform_sphere(rng);
0216         printf("//test_uniform_sphere dir (%10.4f %10.4f %10.4f) \n", dir.x, dir.y, dir.z );
0217     }
0218     return 0 ;
0219 }
0220
0221
0222 /**
0223 QSim_MockTest::propagate_at_boundary_manual
0224 ---------------------------------------------
0225
0226                n
0227            i   :   r
0228             \  :  /
0229              \ : /
0230               \:/
0231           -----+-----------
0232                 \
0233                  \
0234                   \
0235                    t
0236
0237 Consider mom is some direction, say +Z::
0238
0239    (0, 0, 1)
0240
0241 There is a circle of vectors that are perpendicular
0242 to that mom, all in the XY plane::
0243
0244    ( cos(phi), sin(phi), 0 )    phi 0->2pi
0245
0246 Clearly the dot product if that and +Z is zero.
0247
0248 **/
0249 inline int QSim_MockTest::propagate_at_boundary_manual()
0250 {
0251     float3 mom = normalize(make_float3(1.f, 0.f, -1.f)) ;
0252
0253     std::cout
0254         << " QSim_MockTest::propagate_at_boundary "
0255         << " mom " << mom
0256         << std::endl
0257         ;
0258
0259     quad2 prd ;
0260     setup_prd(prd) ;
0261
0262     sctx ctx ;
0263     ctx.prd = &prd ;
0264
0265     sstate& s = ctx.s ;
0266     s.material1.x = 1.0f ;
0267     s.material2.x = 1.5f ;
0268     // The two RINDEX are the only state needed for propagate_at_boundary
0269     // TODO: get these in mocked up way from the bnd texture using the wavelength
0270
0271
0272     sphoton& p = ctx.p ;
0273     p.zero();
0274     p.wavelength = 440.f ;
0275
0276     unsigned flag = 0 ;
0277     const int N = 16 ;
0278
0279     std::vector<sphoton> pp(N*2) ;
0280
0281     for(int i=0 ; i < N ; i++)
0282     {
0283         float frac_twopi = float(i)/float(N)  ;
0284
0285         p.mom = mom ;
0286         p.set_polarization(frac_twopi) ;
0287
0288         sphoton p0(p) ;
0289         int ctrl = sim->propagate_at_boundary(flag, rng, ctx) ;
0290
0291         pp[i*2+0 ] = p0 ;
0292         pp[i*2+1 ] = p ;
0293
0294         std::cout
0295             << " flag " << OpticksPhoton::Flag(flag)
0296             << " ctrl " <<  sflow::desc(ctrl)
0297             << std::endl
0298             << " p0 " << p0.descDir()
0299             << std::endl
0300             << " p  " << p.descDir()
0301             << std::endl
0302             ;
0303     }
0304
0305     NP* a = NP::Make<float>(N,2,4,4) ;
0306     a->read2<float>( (float*)pp.data() );
0307     a->save("$FOLD/pp.npy");
0308     return 0 ;
0309 }
0310
0311
0312 inline int QSim_MockTest::propagate_at_boundary()
0313 {
0314     std::cout
0315         << "QSim_MockTest::propagate_at_boundary"
0316         << std::endl
0317         ;
0318
0319     quad2 prd ;
0320     setup_prd(prd) ;
0321
0322     sctx ctx ;
0323     ctx.idx = 0 ;
0324     ctx.prd = &prd ;
0325
0326     sphoton& p = ctx.p ;
0327     setup_photon(p);
0328
0329     std::cout << "p0 " << p << std::endl ;
0330
0331     float cosTheta = -dot( p.mom, *prd.normal() );
0332
0333     sim->bnd->fill_state(ctx.s, boundary, ctx.p.wavelength, cosTheta, ctx.idx );
0334
0335
0336     unsigned flag = 0 ;
0337     int ctrl = sim->propagate_at_boundary(flag, rng, ctx) ;
0338
0339     std::cout
0340         << " flag " << OpticksPhoton::Flag(flag)
0341         << " ctrl " <<  sflow::desc(ctrl)
0342         << std::endl
0343         ;
0344
0345     std::cout << "p1 " << p << std::endl ;
0346
0347     return 0 ;
0348 }
0349
0350
0351 inline void QSim_MockTest::setup_prd( quad2& prd )
0352 {
0353     float3 nrm = make_float3(0.f, 0.f, 1.f ); // surface normal in +z direction
0354
0355     prd.q0.f.x = nrm.x ;
0356     prd.q0.f.y = nrm.y ;
0357     prd.q0.f.z = nrm.z ;
0358     prd.q0.f.w = 0.1f ;    // distance
0359
0360     prd.q1.f.x = 0.5f ;  // lposcost : local position cosTheta of intersect
0361     prd.q1.u.y = 0u ;
0362     prd.q1.u.z = 1001 ;   // lpmtid : sensor_identifier (< 17612 )
0363     prd.q1.u.w = boundary ;
0364 }
0365 inline void QSim_MockTest::setup_photon( sphoton& p)
0366 {
0367     p.zero();
0368     p.pos = make_float3( 0.f, 0.f, 0.f );
0369     p.mom = normalize(make_float3(1.f, 0.f, -1.f));
0370     p.pol = normalize(make_float3(0.f, 1.f,  0.f));
0371     p.wavelength = 440.f ;
0372 }
0373
0374
0375 #ifdef WITH_CUSTOM4
0376
0377 /**
0378 QSim_MockTest::propagate_at_surface_CustomART_manual
0379 -----------------------------------------------------
0380
0381 Lower level version with manual sstate filling
0382
0383 **/
0384
0385 inline int QSim_MockTest::propagate_at_surface_CustomART_manual()
0386 {
0387     quad2 prd ;
0388     setup_prd(prd) ;
0389
0390     sctx ctx ;
0391     ctx.prd = &prd ;
0392
0393     sphoton& p = ctx.p ;
0394     setup_photon(p);
0395
0396     // lpmtcat doesnt matter as Pyrex and Vacuum are same for all of them
0397     float n0 = sim->pmt->get_lpmtcat_rindex_wl( 0, 0, 0, p.wavelength );
0398     float n3 = sim->pmt->get_lpmtcat_rindex_wl( 0, 3, 0, p.wavelength );
0399
0400     std::cout
0401         << "QSim_MockTest::propagate_at_surface_CustomART_manual "
0402         << " boundary " << boundary
0403         << " prd.q1.u.w " << prd.q1.u.w
0404         << " prd.q0.f:nrm " << prd.q0.f
0405         << " p.mom " << p.mom
0406         << " n0 " << std::fixed << std::setw(10) << std::setprecision(4) << n0
0407         << " n3 " << std::fixed << std::setw(10) << std::setprecision(4) << n3
0408         << std::endl
0409         ;
0410
0411     ctx.s.material1.x = n0 ;  // Pyrex RINDEX
0412     ctx.s.material2.x = n3 ;  // Vacuum RINDEX
0413
0414     unsigned flag = 0 ;
0415     int ctrl = sim->propagate_at_surface_CustomART(flag, rng, ctx) ;
0416
0417     std::cout
0418         << "QSim_MockTest::propagate_at_surface_CustomART_manual "
0419         << " flag " << flag << " : " << OpticksPhoton::Flag(flag)
0420         << " ctrl " << ctrl << " : " << sflow::desc(ctrl)
0421         << std::endl
0422         ;
0423     return 0 ;
0424 }
0425 #endif
0426
0427
0428
0429 inline int QSim_MockTest::fill_state()
0430 {
0431     std::cout
0432         << "QSim_MockTest::fill_state"
0433         << " boundary " << boundary
0434         << std::endl
0435         ;
0436
0437     sctx ctx ;
0438     ctx.p.wavelength = 440.f  ; // nm
0439     ctx.idx = 0 ;
0440
0441     for(float cosTheta=-1.f ; cosTheta <= 1.f ; cosTheta+= 2.f )
0442     {
0443         sim->bnd->fill_state(ctx.s, boundary, ctx.p.wavelength, cosTheta, ctx.idx );
0444         std::cout
0445             << " cosTheta " << cosTheta
0446             << std::endl
0447             << " ctx.s "
0448             << std::endl
0449             << ctx.s
0450             << std::endl
0451             ;
0452     }
0453     return 0 ;
0454 }
0455
0456
0457 /**
0458 QSim_MockTest::propagate
0459 --------------------------
0460
0461 Does qsim::propagate_to_boundary and branches to the
0462 appropriate qsim::propagate_at_boundary method.
0463
0464 **/
0465
0466 inline int QSim_MockTest::propagate()
0467 {
0468     quad2 prd ;
0469     setup_prd(prd) ;
0470
0471     sctx ctx ;
0472     ctx.idx = 0 ;
0473     ctx.prd = &prd ;
0474
0475     sphoton& p = ctx.p ;
0476     setup_photon(p);
0477
0478     int bounce = 0 ;
0479
0480     sphoton p0(p) ;
0481     sim->propagate(bounce, rng, ctx );
0482     sphoton p1(p) ;
0483
0484     std::cout << "p0 " << p0 << std::endl ;
0485     std::cout << "p1 " << p1 << std::endl ;
0486
0487     return 0 ;
0488 }
0489
0490 /**
0491 QSim_MockTest::SmearNormal
0492 ---------------------------
0493
0494 This MOCK_CUDA test of qsim::SmearNormal_SigmaAlpha qsim::SmearNormal_Polish
0495 is similar to sysrap/tests/S4OpBoundaryProcessTest.sh
0496
0497 +------+-------------+
0498 | chk  |  value      |
0499 +======+=============+
0500 |  0   | sigma_alpha |
0501 +------+-------------+
0502 |  1   | polish      |
0503 +------+-------------+
0504 **/
0505
0506 inline void QSim_MockTest::SmearNormal(int chk, float value)
0507 {
0508     float3 direct = make_float3(0.f, 0.f, -1.f );
0509     float3 normal = make_float3(0.f, 0.f,  1.f );
0510
0511     int ni = num ;
0512     int nj = 4 ;
0513
0514     a = NP::Make<float>( ni, nj );
0515
0516     a->set_meta<std::string>("source", "QSim_MockTest.sh" );
0517     a->set_meta<std::string>("normal", scuda::serialize(normal) );
0518     a->set_meta<std::string>("direct", scuda::serialize(direct) );
0519     a->set_meta<float>("value", value );
0520     a->set_meta<std::string>("valuename", chk == 0 ? "sigma_alpha" : "polish"  );
0521     a->names.push_back( chk == 0 ? "SmearNormal_SigmaAlpha" : "SmearNormal_Polish" );
0522
0523
0524     sctx ctx ;
0525
0526
0527     float* aa = a->values<float>();
0528     for(int i=0 ; i < ni ; i++)
0529     {
0530         rng.setSequenceIndex(i);
0531         ctx.idx = i ;
0532         float3* smeared_normal = (float3*)(aa + i*nj + 0) ;
0533         switch(chk)
0534         {
0535             case 0: sim->SmearNormal_SigmaAlpha(rng, smeared_normal, &direct, &normal, value, ctx ); break ;
0536             case 1: sim->SmearNormal_Polish(    rng, smeared_normal, &direct, &normal, value, ctx ); break ;
0537         }
0538     }
0539
0540     a->save("$FOLD/q.npy");
0541 }
0542
0543
0544
0545 inline int QSim_MockTest::SmearNormal_SigmaAlpha_one()
0546 {
0547     rng.setSequenceIndex(0);
0548
0549     float sigma_alpha = 0.1f ;
0550     float3 direct = make_float3(0.f, 0.f, -1.f );
0551     float3 normal = make_float3(0.f, 0.f,  1.f );
0552
0553     float3 smeared_normal = make_float3( 0.f, 0.f, 0.f ) ;
0554     sctx ctx ;
0555     ctx.idx = 0 ;
0556
0557     sim->SmearNormal_SigmaAlpha(rng, &smeared_normal, &direct, &normal, sigma_alpha, ctx );
0558     return 0 ;
0559 }
0560
0561 inline int QSim_MockTest::run()
0562 {
0563     if(     strcmp(CHECK,"smear_normal_sigma_alpha")==0) SmearNormal(0, 0.1f) ;
0564     else if(strcmp(CHECK,"smear_normal_polish")==0)     SmearNormal(1, 0.8f) ;
0565     else
0566     {
0567         std::cerr
0568             << "QSim_MockTest::run"
0569             << " CHECK " << ( CHECK ? CHECK : "-" )
0570             << " UNHANDLED "
0571             << std::endl
0572             ;
0573     }
0574     return 0 ;
0575 }
0576
0577
0578 inline int QSim_MockTest::main()
0579 {
0580     bool ALL = strcmp(TEST, "ALL") == 0 ;
0581     int rc = 0 ;
0582
0583     if(ALL||0==strcmp(TEST, "propagate_at_surface_CustomART_manual"))
0584     {
0585 #ifdef WITH_CUSTOM4
0586        rc += propagate_at_surface_CustomART_manual() ;
0587 #else
0588        std::cout << "NOT-WITH_CUSTOM4 " << std::endl ;
0589        rc += 1 ;
0590 #endif
0591     }
0592
0593     if(ALL||0==strcmp(TEST,"propagate_at_boundary_manual")) rc += propagate_at_boundary_manual() ;
0594     if(ALL||0==strcmp(TEST,"fill_state"))                   rc += fill_state() ;
0595     if(ALL||0==strcmp(TEST,"propagate_at_boundary"))        rc += propagate_at_boundary() ;
0596     if(ALL||0==strcmp(TEST,"propagate"))                    rc += propagate() ;
0597     if(ALL||0==strcmp(TEST,"SmearNormal_SigmaAlpha_one"))   rc += SmearNormal_SigmaAlpha_one() ;
0598     if(ALL||0==strcmp(TEST,"run"))                          rc += run();
0599
0600     return rc ;
0601 }
0602
0603
0604
0605 int main(int argc, char** argv)
0606 {
0607     QSim_MockTest t ;
0608     std::cout << t.desc() ;
0609     return t.main() ;
0610 }
0611