eic-opticks/sysrap/SGenstep.h

0001 #pragma once
0002 /**
0003 SGenstep.h : genstep static utilities
0004 ========================================
0005
0006 TODO: incorporate all/most of the genstep methods from SEvent.{hh,cc}
0007       here to avoid near duplicated impls
0008
0009
0010 Used by:
0011
0012 sysrap/SGenerate.h
0013    GetGencode
0014
0015 sysrap/SFrameGenstep.cc
0016    GenstepID, ConfigureGenstep, MakeArray, GridAxes, GridAxesName
0017
0018 sysrap/SCenterExtentGenstep.cc
0019    [ON WAY OUT] GridAxes, GridAxesName
0020
0021 sysrap/SEvent.cc
0022    ConfigureGenstep, MakeArray used by SEvent::MakeCountGenstep
0023
0024 qudarap/QEvt.cc
0025    Check, Desc, GetGencode used by QEvt::setGenstepUpload
0026
0027
0028 +---------+----------+-------------+
0029 | qty     | NumPy    |  Note       |
0030 +=========+==========+=============+
0031 | q0.i.x  | [:,0,0]  |  gencode    |
0032 +---------+----------+-------------+
0033 | q0.i.w  | [:,0,3]  | num_photon  |
0034 +---------+----------+-------------+
0035
0036
0037 **/
0038
0039 #include <string>
0040 #include <vector>
0041 #include <limits>
0042 #include <csignal>
0043
0044 struct quad6 ;
0045 struct NP ;
0046 struct sslice ;
0047
0048 #include "sxyz.h"
0049 #include "SYSRAP_API_EXPORT.hh"
0050
0051 struct SYSRAP_API SGenstep
0052 {
0053     static constexpr const size_t MAX_SLOT_PER_SLICE = std::numeric_limits<size_t>::max();
0054     static constexpr const char* XYZ_ = "XYZ" ;
0055     static constexpr const char* YZ_  = "YZ" ;
0056     static constexpr const char* XZ_  = "XZ" ;
0057     static constexpr const char* XY_  = "XY" ;
0058
0059     static const char* GridAxesName( int gridaxes );
0060     static int GridAxes(int nx, int ny, int nz);
0061     static unsigned GenstepID( int ix, int iy, int iz, int iw=0 ) ;
0062
0063     static void ConfigureGenstep( quad6& gs,  int gencode, int gridaxes, int gsid, int photons_per_genstep );
0064     static int GetGencode( const quad6& gs ) ;
0065     static size_t GetNumPhoton( const quad6& gs ) ;
0066     static void SetNumPhoton( quad6& gs, size_t num );
0067
0068     static const quad6& GetGenstep(const NP* gs, unsigned gs_idx );
0069     static size_t GetGenstepSlices(std::vector<sslice>& slice, const NP* gs, size_t max_slot );
0070     static void CheckGenstepSlices(const std::vector<sslice>& slice, const NP* gs, size_t max_slot );
0071
0072     static int GetGencode( const quad6* qq, unsigned gs_idx  );
0073     static int GetGencode(    const NP* gs, unsigned gs_idx  );
0074
0075     static size_t GetNumPhoton( const quad6* qq, unsigned gs_idx  );
0076     static size_t GetNumPhoton( const NP* gs, unsigned gs_idx  );
0077
0078     static size_t GetPhotonTotal( const NP* gs );
0079     static size_t GetPhotonTotal( const NP* gs, int gs_start, int gs_stop  );
0080
0081     static void Check(const NP* gs);
0082     static NP* MakeArray(const std::vector<quad6>& gs );
0083
0084     template<typename T>
0085     static NP* MakeTestArray(const std::vector<T>& num_ph);
0086
0087     template<typename T>
0088     static std::string DescNum(const std::vector<T>& num_ph);
0089
0090     static std::string Desc(const NP* gs, int edgeitems);
0091
0092 };
0093
0094
0095 #include <cassert>
0096
0097 #include "scuda.h"
0098 #include "squad.h"
0099 #include "sphoton.h"
0100 #include "sc4u.h"
0101 #include "sslice.h"
0102
0103 #include "OpticksGenstep.h"
0104 #include "NP.hh"
0105
0106
0107
0108 inline const char* SGenstep::GridAxesName( int gridaxes)  // static
0109 {
0110     const char* s = nullptr ;
0111     switch( gridaxes )
0112     {
0113         case XYZ: s = XYZ_ ; break ;
0114         case YZ:  s = YZ_  ; break ;
0115         case XZ:  s = XZ_  ; break ;
0116         case XY:  s = XY_  ; break ;
0117     }
0118     return s ;
0119 }
0120
0121 /**
0122 SGenstep::GridAxes
0123 --------------------------
0124
0125 The nx:ny:nz dimensions of the grid are used to classify it into::
0126
0127     YZ
0128     XZ
0129     XY
0130     XYZ
0131
0132 For a planar grid one of the nx:ny:nz grid dimensions is zero.
0133 XYZ is a catch all for non-planar grids.
0134
0135 **/
0136
0137 inline int SGenstep::GridAxes(int nx, int ny, int nz)  // static
0138 {
0139     int gridaxes = XYZ ;
0140     if( nx == 0 && ny > 0 && nz > 0 )
0141     {
0142         gridaxes = YZ ;
0143     }
0144     else if( nx > 0 && ny == 0 && nz > 0 )
0145     {
0146         gridaxes = XZ ;
0147     }
0148     else if( nx > 0 && ny > 0 && nz == 0 )
0149     {
0150         gridaxes = XY ;
0151     }
0152     return gridaxes ;
0153 }
0154
0155 /**
0156 SGenstep::GenstepID
0157 -------------------
0158
0159 Pack four signed integers (assumed to be in char range -128 to 127)
0160 into a 32 bit unsigtned char using C4U uniform.
0161
0162 **/
0163
0164 inline unsigned SGenstep::GenstepID( int ix, int iy, int iz, int iw )
0165 {
0166     C4U gsid ;   // sc4u.h
0167
0168     gsid.c4.x = ix ;
0169     gsid.c4.y = iy ;
0170     gsid.c4.z = iz ;
0171     gsid.c4.w = iw ;
0172
0173     return gsid.u ;
0174 }
0175
0176
0177
0178
0179 /**
0180 SGenstep::ConfigureGenstep
0181 ---------------------------
0182
0183 TODO: pack enums to make room for a photon_offset
0184
0185 * gsid was MOVED from (1,3) to (0,2) when changing genstep to carry transform
0186
0187 **/
0188
0189 inline void SGenstep::ConfigureGenstep( quad6& gs,  int gencode, int gridaxes, int gsid, int photons_per_genstep )
0190 {
0191     assert( gencode == OpticksGenstep_TORCH || gencode == OpticksGenstep_FRAME );
0192     assert( gridaxes == XYZ ||  gridaxes == YZ || gridaxes == XZ || gridaxes == XY );
0193
0194     gs.q0.i.x = gencode ;
0195     gs.q0.i.y = gridaxes ;
0196     gs.q0.u.z = gsid ;
0197     gs.q0.i.w = photons_per_genstep ;
0198 }
0199
0200 inline int SGenstep::GetGencode( const quad6& gs )
0201 {
0202     return gs.q0.i.x ;
0203 }
0204 inline size_t SGenstep::GetNumPhoton( const quad6& gs )
0205 {
0206     return gs.q0.u.w ;
0207 }
0208
0209
0210
0211
0212
0213 inline void SGenstep::SetNumPhoton( quad6& gs, size_t num )
0214 {
0215     bool num_allowed = num <= MAX_SLOT_PER_SLICE ;
0216     assert( num_allowed );
0217     if(!num_allowed) std::raise(SIGINT);
0218
0219     gs.q0.i.w = num ;
0220 }
0221
0222
0223
0224
0225 inline const quad6& SGenstep::GetGenstep(const NP* gs, unsigned gs_idx )
0226 {
0227     Check(gs);
0228     assert( int(gs_idx) < gs->shape[0] );
0229     quad6* qq = (quad6*)gs->bytes() ;
0230     const quad6& q = qq[gs_idx] ;
0231     return q ;
0232 }
0233
0234
0235 /**
0236 SGenstep::GetGenstepSlices
0237 --------------------------
0238
0239 Populates the genstep slice vector with gs[start:stop] slices
0240 that do not exceed max_photon for the total number of photons.
0241
0242 The third argument *max_slot* was formerly incorrectly named *max_photon*,
0243 which is wrong because its the VRAM constrained *max_slot* which is relevant
0244 to how the work needed to be sliced, not the arbitrary and very large
0245 *max_photon* which is not VRAM constrained.
0246
0247 **/
0248
0249 inline size_t SGenstep::GetGenstepSlices(std::vector<sslice>& slice, const NP* gs, size_t max_slot )
0250 {
0251     Check(gs);
0252     int num_gs = gs ? gs->shape[0] : 0 ;
0253     const quad6* qq = gs ? (quad6*)gs->cvalues<float>() : nullptr ;
0254
0255     sslice sl = {};
0256
0257     sl.gs_start = 0 ;   // gs index starting the slice
0258     sl.gs_stop = 0 ;    // gs index stopping the slice, ie one beyond the last index, python style gs[start:stop]
0259     sl.ph_offset = 0 ;  // total photons before the start of this slice
0260     sl.ph_count = 0 ;   // total photons within this slice
0261
0262
0263     size_t tot_ph = 0 ;
0264
0265     bool extend_slice = false ;
0266
0267     for(int i=0 ; i < num_gs ; i++)
0268     {
0269         bool last_gs = i == num_gs - 1;
0270
0271         const quad6& q = qq[i] ;
0272         size_t num_ph = GetNumPhoton(q);
0273
0274         bool slice_too_big = num_ph > MAX_SLOT_PER_SLICE ;
0275         if(slice_too_big ) std::cerr
0276             << "SGenstep::GetGenstepSlices FATAL "
0277             << " slice_too_big " << ( slice_too_big ? "YES" : "NO " )
0278             << " i " << i
0279             << " num_gs " << num_gs
0280             << " num_ph " << num_ph
0281             << " MAX_SLOT_PER_SLICE " << MAX_SLOT_PER_SLICE
0282             << "\n"
0283             ;
0284         assert( !slice_too_big );
0285
0286
0287         tot_ph += num_ph ;
0288
0289         size_t cand = sl.ph_count + num_ph ;
0290         extend_slice = cand <= max_slot ;
0291
0292         if(0) std::cout
0293             << "SGenstep::GetGenstepSlices"
0294             << " i " << std::setw(4) << i
0295             << " sl.ph_count " << std::setw(7) << sl.ph_count
0296             << " num_ph " << std::setw(7) << num_ph
0297             << " max_slot " << std::setw(7) << max_slot
0298             << " cand  " << std::setw(7) << num_ph
0299             << " extend_slice " << ( extend_slice ? "YES" : "NO " )
0300             << " last_gs " << ( last_gs ? "YES" : "NO " )
0301             << "\n"
0302             ;
0303
0304         if(extend_slice)  // genstep fits within limit
0305         {
0306             sl.gs_stop = i+1 ;
0307             sl.ph_count += num_ph ;
0308         }
0309         else  // genstep does not fit, so close slice and start another
0310         {
0311             sl.gs_stop = i ;
0312             slice.push_back(sl) ;
0313
0314             sl.ph_count = 0 ;          // back to zero
0315             sl.ph_count += num_ph ;
0316             sl.gs_start = i ;
0317             sl.gs_stop = i+1 ;
0318         }
0319         if(last_gs) slice.push_back(sl) ;
0320     }
0321     sslice::SetOffset(slice);
0322
0323     CheckGenstepSlices(slice, gs, max_slot);
0324
0325     return tot_ph ;
0326 }
0327
0328 inline void SGenstep::CheckGenstepSlices(const std::vector<sslice>& slice, const NP* gs, size_t max_slot )
0329 {
0330     size_t gs_tot = GetPhotonTotal(gs);
0331     size_t sl_tot = sslice::TotalPhoton(slice);
0332
0333     bool tot_consistent = gs_tot == sl_tot ;
0334
0335     if(!tot_consistent) std::cerr
0336        << "SGenstep::CheckGenstepSlices"
0337        << " gs_tot " << gs_tot
0338        << " sl_tot " << sl_tot
0339        << " tot_consistent " << ( tot_consistent ? "YES" : "NO " )
0340        << "\n"
0341        ;
0342     assert(tot_consistent);
0343
0344
0345     const quad6* qq = gs ? (quad6*)gs->cvalues<float>() : nullptr ;
0346     int num_sl = slice.size();
0347     for(int i=0 ; i < num_sl ; i++)
0348     {
0349         const sslice& sl = slice[i];
0350
0351         size_t sl_ph_count = 0 ;
0352         for(size_t j=sl.gs_start ; j < sl.gs_stop ; j++ )
0353         {
0354             const quad6& q = qq[j] ;
0355             size_t num_ph = GetNumPhoton(q);
0356             sl_ph_count += num_ph ;
0357         }
0358         bool ph_count_expected = sl_ph_count == sl.ph_count ;
0359
0360         if(!ph_count_expected) std::cerr
0361             << "SGenstep::CheckGenstepSlices"
0362             << " i " << i
0363             << " sl.ph_count " << sl.ph_count
0364             << " sl_ph_count " << sl_ph_count
0365             << "\n"
0366             ;
0367
0368         assert( ph_count_expected );
0369         assert( sl_ph_count <= max_slot ); // total for slice should never exceed max_slot
0370     }
0371 }
0372
0373
0374
0375 inline int SGenstep::GetGencode( const quad6* qq, unsigned gs_idx  ) // static
0376 {
0377     if( qq == nullptr ) return OpticksGenstep_INVALID ;
0378     const quad6& q = qq[gs_idx] ;
0379     return GetGencode(q) ;
0380 }
0381 inline int SGenstep::GetGencode( const NP* gs, unsigned gs_idx  ) // static
0382 {
0383     const quad6& q = GetGenstep(gs, gs_idx);
0384     return GetGencode(q) ;
0385 }
0386
0387
0388 inline size_t SGenstep::GetNumPhoton( const quad6* qq, unsigned gs_idx  ) // static
0389 {
0390     const quad6& q = qq[gs_idx] ;
0391     return GetNumPhoton(q) ;
0392 }
0393
0394 inline size_t SGenstep::GetNumPhoton( const NP* gs, unsigned gs_idx  ) // static
0395 {
0396     if( gs == nullptr) return 0 ;
0397     const quad6& q = GetGenstep(gs, gs_idx);
0398     return GetNumPhoton(q) ;
0399 }
0400
0401
0402 inline size_t SGenstep::GetPhotonTotal( const NP* gs ) // static
0403 {
0404     int num_gs = gs ? gs->shape[0] : 0 ;
0405     return GetPhotonTotal( gs, 0, num_gs );
0406 }
0407
0408 /**
0409 Genstep::GetPhotonTotal
0410 ------------------------
0411
0412 Returns total photon in genstep slice gs[start:stop]
0413
0414 **/
0415
0416 inline size_t SGenstep::GetPhotonTotal( const NP* gs, int gs_start, int gs_stop ) // static
0417 {
0418     int num_gs = gs ? gs->shape[0] : 0 ;
0419     assert( gs_start <= num_gs );
0420     assert( gs_stop  <= num_gs );
0421     size_t tot = 0 ;
0422     for(int i=gs_start ; i < gs_stop ; i++ ) tot += GetNumPhoton(gs, i );
0423     return tot ;
0424 }
0425
0426
0427
0428
0429
0430 inline NP* SGenstep::MakeArray(const std::vector<quad6>& gs ) // static
0431 {
0432     assert( gs.size() > 0);
0433     NP* a = NP::Make<float>( gs.size(), 6, 4 );
0434     a->read2<float>( (float*)gs.data() );
0435     Check(a);
0436     return a ;
0437 }
0438
0439 /**
0440 SGenstep::MakeTestArray
0441 ------------------------
0442
0443 Make test array filled with gensteps containing only num_photon
0444
0445 **/
0446
0447 template<typename T>
0448 inline NP* SGenstep::MakeTestArray(const std::vector<T>& num_ph) // static
0449 {
0450     int num_gs = num_ph.size();
0451     NP* gs = NP::Make<float>(num_gs, 6, 4);
0452     quad6* qq = (quad6*)gs->values<float>() ;
0453     for(int i=0 ; i < num_gs ; i++) SGenstep::SetNumPhoton(qq[i], num_ph[i]);
0454     Check(gs);
0455     return gs ;
0456 }
0457
0458
0459 template<typename T>
0460 inline std::string SGenstep::DescNum(const std::vector<T>& num_ph)
0461 {
0462     T tot = 0 ;
0463     std::stringstream ss ;
0464     ss << "SGenstep::DescNum\n" ;
0465     for(int i=0 ; i < int(num_ph.size()) ; i++ )
0466     {
0467         T num = num_ph[i] ;
0468         ss
0469             << std::setw(3) << i << " : "
0470             << " num : " << std::setw(7) << num
0471             << " tot : " << std::setw(7) << tot
0472             << "\n"
0473             ;
0474
0475         tot += num ;
0476     }
0477
0478     ss << " Grand total : " << tot << "\n\n" ;
0479
0480     std::string str = ss.str() ;
0481     return str ;
0482 }
0483
0484
0485
0486 inline void SGenstep::Check(const NP* gs)  // static
0487 {
0488     if( gs == nullptr ) return ;
0489     assert( gs->uifc == 'f' && gs->ebyte == 4 );
0490     assert( gs->has_shape(-1, 6, 4) );
0491 }
0492
0493 inline std::string SGenstep::Desc(const NP* gs, int edgeitems) // static
0494 {
0495     int num_genstep = gs ? gs->shape[0] : 0 ;
0496
0497     quad6* gs_v = gs ? (quad6*)gs->cvalues<float>() : nullptr ;
0498     std::stringstream ss ;
0499     ss << "SGenstep::DescGensteps num_genstep " << num_genstep << " (" ;
0500
0501     int total = 0 ;
0502     for(int i=0 ; i < num_genstep ; i++)
0503     {
0504         const quad6& _gs = gs_v[i];
0505         unsigned gs_pho = _gs.q0.u.w  ;
0506
0507         if( i < edgeitems || i > num_genstep - edgeitems ) ss << gs_pho << " " ;
0508         else if( i == edgeitems )  ss << "... " ;
0509
0510         total += gs_pho ;
0511     }
0512     ss << ") total " << total  ;
0513     std::string s = ss.str();
0514     return s ;
0515 }
0516
0517