eic-opticks/sysrap/SRecord.h

0001 #pragma once
0002 /**
0003 SRecord.h
0004 ==============
0005
0006 Used from SGLFW_Evt.h
0007
0008 **/
0009
0010 #include "NP.hh"
0011
0012
0013 #include "ssys.h"
0014 #include "spath.h"
0015 #include "sphoton.h"
0016 #include "sseq_record.h"
0017 #include "sstr.h"
0018 #include "scuda.h"
0019
0020
0021
0022 struct SRecord
0023 {
0024     static constexpr const char* NAME = "record.npy" ;
0025     static constexpr const char* RPOS_SPEC = "4,GL_FLOAT,GL_FALSE,64,0,false";
0026     static constexpr const char* _level = "SRecord__level" ;
0027     static int level ;
0028
0029     NP* record;
0030     int record_first ;
0031     int record_count ;  // all step points across all photon
0032
0033     float4 mn = {} ;
0034     float4 mx = {} ;
0035     float4 ce = {} ;
0036
0037     static NP*      LoadArray( const char* _fold, const char* _slice );
0038     static SRecord* Load(      const char* _fold, const char* _slice=nullptr );
0039
0040     SRecord(NP* record);
0041     void init() ;
0042
0043     const float* get_mn() const ;
0044     const float* get_mx() const ;
0045     const float* get_ce() const ;
0046
0047     const float get_t0() const ;
0048     const float get_t1() const ;
0049     std::string desc() const ;
0050
0051     void getPhotonAtTime( std::vector<sphoton>* pp, std::vector<quad4>* qq, const char* iprt ) const ;
0052     static bool FindInterpolatedPhotonFromRecordAtTime( sphoton& p, const std::vector<sphoton>& point,  float t);
0053
0054     NP*  getPhotonAtTime(   const char* iprt ) const;
0055     NP*  getSimtraceAtTime( const char* iprt ) const;
0056
0057
0058 };
0059
0060
0061 int SRecord::level = ssys::getenvint(_level,0) ;
0062
0063
0064 /**
0065 SRecord::LoadArray
0066 -------------------
0067
0068 Several forms of slice selection are handled.
0069
0070 The slice argument supports envvar tokens like "$AFOLD_RECORD_SLICE"
0071 that yield a selection string or direct such strings.
0072
0073 0. seqhis history string, eg::
0074
0075     TO BT BT BT BT BR BT BT BT BT BT BT SC BT BT BT BT SD
0076
0077 1. "where" selection, eg pick photon records with y coordinate
0078    of the first step point less than a value::
0079
0080      "[:,0,0,1] < -1.0"
0081
0082 2. "indexSlice" selection, eg every 1000th or the first 10.
0083    This uses partial loading of items from files which enables working
0084    with very large record files, eg 66GB. Example slices::
0085
0086      "[::1000]"
0087      "[:10]"
0088
0089 3. path to index array eg "/tmp/w54.npy" : create that file from ipython with::
0090
0091     In [39]: w54
0092     Out[39]: array([71464, 71485, 71663, 71678, 71690, 71780, 71782, 71798, 71799, 71803, 71810, 71914, 71918, 71919, 71975, 71978, 72019, 72026, 72032, 73149, 73212, 73344, 73372])
0093     In [40]: np.save("/tmp/w54.npy", w54)
0094
0095
0096 **/
0097
0098 inline NP* SRecord::LoadArray(const char* _fold, const char* _slice )
0099 {
0100     const char* path = spath::Resolve(_fold, NAME);
0101     bool looks_unresolved = spath::LooksUnresolved(path, _fold);
0102     if(looks_unresolved)
0103     {
0104         std::cout
0105             << "SRecord::LoadArray"
0106             << " FAILED : DUE TO MISSING ENVVAR\n"
0107             << " _fold [" << ( _fold ? _fold : "-" ) << "]\n"
0108             << " path ["  << (  path ?  path : "-" ) << "]\n"
0109             << " looks_unresolved " << ( looks_unresolved ? "YES" : "NO " )
0110             << "\n"
0111             ;
0112         return nullptr ;
0113     }
0114
0115     NP* a = nullptr ;
0116
0117     if( _slice == nullptr )
0118     {
0119         a = NP::Load(path);
0120         a->set_meta<std::string>("SRecord__LoadArray", "NP::Load");
0121     }
0122     else if(sseq_record::LooksLikeRecordSeqSelection(_slice)) // can be indirect via envvar, once resolved look for starting with "CK/SI/TO"
0123     {
0124         a = sseq_record::LoadRecordSeqSelection(_fold, _slice);
0125         a->set_meta<std::string>("SRecord__LoadArray_METHOD", "sseq_record::LoadRecordSeqSelection");
0126         a->set_meta<std::string>("SRecord__LoadArray_SLICE", _slice );
0127     }
0128     else if(NP::LooksLikeWhereSelection(_slice)) // can be indirect via envvar, once resolved look for "<" or ">"
0129     {
0130         a = NP::LoadThenSlice<float>(path, _slice);
0131         a->set_meta<std::string>("SRecord__LoadArray_METHOD", "NP::LoadThenSlice");
0132         a->set_meta<std::string>("SRecord__LoadArray_SLICE", _slice );
0133     }
0134     else    // eg "[0:100]" OR "/tmp/w54.npy" OR "/tmp/w54.npy[0:10]"
0135     {
0136         a = NP::LoadSlice(path, _slice);
0137         a->set_meta<std::string>("SRecord__LoadArray_METHOD", "NP::LoadSlice");
0138         a->set_meta<std::string>("SRecord__LoadArray_SLICE", _slice );
0139     }
0140     return a ;
0141 }
0142
0143
0144 inline SRecord* SRecord::Load(const char* _fold, const char* _slice )
0145 {
0146     NP* _record = LoadArray(_fold, _slice);
0147     return _record ? new SRecord(_record) : nullptr ;
0148 }
0149
0150
0151
0152
0153 inline SRecord::SRecord(NP* _record)
0154     :
0155     record(_record),
0156     record_first(0),
0157     record_count(0)
0158 {
0159     init();
0160 }
0161
0162 /**
0163 SRecord::init
0164 -------------------
0165
0166 Expected shape of record array like (10000, 10, 4, 4)
0167
0168 **/
0169
0170
0171 inline void SRecord::init()
0172 {
0173     assert(record->shape.size() == 4);
0174     bool is_compressed = record->ebyte == 2 ;
0175     assert( is_compressed == false );
0176
0177     record_first = 0 ;
0178     record_count = record->shape[0]*record->shape[1] ;   // all step points across all photon
0179
0180     bool skip_flagmask_zero = true ;
0181     sphoton::MinMaxPost(&mn.x, &mx.x, record, skip_flagmask_zero );
0182     ce = scuda::center_extent( mn, mx );
0183
0184     record->set_meta<float>("x0", mn.x );
0185     record->set_meta<float>("x1", mx.x );
0186
0187     record->set_meta<float>("y0", mn.y );
0188     record->set_meta<float>("y1", mx.y );
0189
0190     record->set_meta<float>("z0", mn.z );
0191     record->set_meta<float>("z1", mx.z );
0192
0193     record->set_meta<float>("t0", mn.w );
0194     record->set_meta<float>("t1", mx.w );
0195
0196     record->set_meta<float>("cx", ce.x );
0197     record->set_meta<float>("cy", ce.y );
0198     record->set_meta<float>("cz", ce.z );
0199     record->set_meta<float>("ce", ce.w );
0200 }
0201
0202
0203
0204 inline const float* SRecord::get_mn() const
0205 {
0206     return &mn.x ;
0207 }
0208 inline const float* SRecord::get_mx() const
0209 {
0210     return &mx.x ;
0211 }
0212 inline const float* SRecord::get_ce() const
0213 {
0214     return &ce.x ;
0215 }
0216
0217 inline const float SRecord::get_t0() const
0218 {
0219     return mn.w ;
0220 }
0221 inline const float SRecord::get_t1() const
0222 {
0223     return mx.w ;
0224 }
0225
0226
0227 inline std::string SRecord::desc() const
0228 {
0229     const char* lpath = record ? record->lpath.c_str() : nullptr ;
0230
0231     std::stringstream ss ;
0232     ss
0233         << "[SRecord.desc\n"
0234         << " lpath [" << ( lpath ? lpath : "-" ) << "]\n"
0235         << std::setw(20) << " mn " << mn
0236         << std::endl
0237         << std::setw(20) << " mx " << mx
0238         << std::endl
0239         << std::setw(20) << " ce " << ce
0240         << std::endl
0241         << std::setw(20) << " record.sstr " << record->sstr()
0242         << std::endl
0243         << std::setw(20) << " record_first " << record_first
0244         << std::endl
0245         << std::setw(20) << " record_count " << record_count
0246         << std::endl
0247         << ( record ? record->descMeta() : "-" )
0248         << std::endl
0249         << "]SRecord.desc\n"
0250         ;
0251     std::string str = ss.str() ;
0252     return str ;
0253 }
0254
0255 /**
0256 SRecord::getPhotonAtTime
0257 --------------------------
0258
0259 Provide photons corresponding to input simulation time(s)
0260 by interpolation between positions in the record array. Momentum and
0261 polarization are taken from the first of the straddling step points.
0262 When not straddling the photon is skipped, corresponding
0263 to it not being alive at the simulation time provided.
0264
0265 Config example from ~/o/CSGOptiX/cxt_precision.sh with multiple times::
0266
0267      export OPTICKS_INPUT_PHOTON=$AFOLD/record.npy
0268      export OPTICKS_INPUT_PHOTON_RECORD_SLICE="TO BT BT BT SA"
0269      export OPTICKS_INPUT_PHOTON_RECORD_TIME=[10:80:10]  # ns
0270
0271 Python analysis is simplified by arranging RECORD_SLICE
0272 and RECORD_TIME such that all photons result in the same number of
0273 staddles by ensuring that the lowest and highest times are
0274 during the lifetime of the photon. Going further, arranging times
0275 to all be inbetween two step points eg "TO->BT" makes the analysis
0276 of intersects from different distances as simple as possible.
0277
0278
0279 Q : Simtrace input uses q0.f and q1.f for pos and mom, like sphoton.h
0280     Simtrace output uses q2.f and q3.f for origin pos and mom.
0281     Why different layout from input and output, as
0282     demonstrated here and in sevent::add_simtrace ?
0283
0284 A : Perhaps to match sphoton for input. When get chance rationalize to use one layout.
0285     This is not urgent, and will require widespread changes across cpp and py.
0286
0287 **/
0288
0289 inline void SRecord::getPhotonAtTime( std::vector<sphoton>* pp, std::vector<quad4>* qq, const char* _iprt ) const
0290 {
0291     NP* iprt = NP::ARange_FromString<float>(_iprt);
0292     const float* tt = iprt ? iprt->cvalues<float>() : nullptr ;
0293
0294     int ni = record->shape[0] ;
0295     int nj = iprt ? iprt->num_items() : 0 ;
0296     int nl = record->shape[1] ;
0297     int item_bytes = record->item_bytes() ; // encompasses multiple step point sphoton
0298
0299     assert( nj > 0 && tt );
0300     int count0 = -1 ;
0301
0302     for(int i=0 ; i < ni ; i++)
0303     {
0304         // populate step point vector
0305         std::vector<sphoton> point(nl) ;
0306         memcpy( point.data(), record->bytes() + i*item_bytes,  item_bytes );
0307
0308         int count = 0 ;
0309         for(int j=0 ; j < nj ; j++) // time loop
0310         {
0311             float t = tt[j] ;
0312             sphoton p = {} ;
0313             bool found = FindInterpolatedPhotonFromRecordAtTime( p, point, t);
0314             if(!found) continue ;
0315
0316             count += 1 ;   // count times at which an interpolated photon was found
0317
0318             quad4 q = {} ;
0319             p.populate_simtrace_input(q);
0320
0321             if(pp) pp->push_back(p);
0322             if(qq) qq->push_back(q);
0323         }
0324
0325         bool consistent_count = count0 == -1 || count0 == count ;
0326         if( count0 == -1 ) count0 = count ;
0327
0328         if(!consistent_count || level > 0) std::cout
0329             << "SRecord::getPhotonAtTime"
0330             << " i " << std::setw(6) << i
0331             << " count0 " << count0
0332             << " count " << count
0333             << " consistent_count " << ( consistent_count ? "YES" : "NO " )
0334             << ( consistent_count ? "" : "(adjust time range to achieve consistency)" )
0335             << " level " << level
0336             << "\n"
0337             ;
0338     }   // i: photon loop
0339
0340
0341     if(level > 0) std::cout
0342          << "SRecord::getPhotonAtTime"
0343          << " [" << _level << "] " << level
0344          << " _iprt " << ( _iprt ? _iprt : "-" )
0345          << " iprt " << ( iprt ? iprt->sstr() : "-" )
0346          << " ni " << ni
0347          << " nj " << nj
0348          << " nl " << nl
0349          << " tt " << ( tt ? "YES" : "NO " )
0350          << " record " << ( record ? record->sstr() : "-" )
0351          << " count0 " << count0
0352          << " pp.size " << ( pp ? int(pp->size()) : -1 )
0353          << " qq.size " << ( qq ? int(qq->size()) : -1 )
0354          << "\n"
0355          ;
0356
0357 }
0358
0359
0360
0361 /**
0362 SRecord::FindInterpolatedPhotonFromRecordAtTime
0363 ------------------------------------------------
0364
0365 Look for record step points with times that straddle the provided time *t*.
0366 If found interpolate the straddling step points to populate *p*.
0367 When a straddle is found returns true.
0368
0369 **/
0370
0371
0372 inline bool SRecord::FindInterpolatedPhotonFromRecordAtTime( sphoton& p, const std::vector<sphoton>& point,  float t)
0373 {
0374     int nl = point.size();
0375     for(int l=1 ; l < nl ; l++)
0376     {
0377         const sphoton& p0 = point[l-1] ;
0378         const sphoton& p1 = point[l] ;
0379         bool straddle = p0.time <= t  && t < p1.time ;
0380
0381         if(level > 3) std::cout
0382             << "SRecord::FindInterpolatedPhotonFromRecordAtTime"
0383             << " nl " << nl
0384             << " l " << std::setw(3) << l
0385             << " t " << std::fixed << std::setw(7) << std::setprecision(3) << t
0386             << " p0.time " << std::fixed << std::setw(7) << std::setprecision(3) << p0.time
0387             << " p1.time " << std::fixed << std::setw(7) << std::setprecision(3) << p1.time
0388             << " straddle " << ( straddle ? "YES" : "NO " )
0389             << "\n"
0390             ;
0391
0392         if(straddle)
0393         {
0394             p = p0 ;  // (pol,mom) from p0
0395             float frac = (t - p0.time)/(p1.time - p0.time) ;
0396             p.pos = lerp( p0.pos, p1.pos, frac ) ; // interpolated position
0397             p.time = t ;
0398             return true ;
0399         }
0400     }
0401     return false ;
0402 }
0403
0404
0405 inline NP* SRecord::getPhotonAtTime( const char* iprt ) const
0406 {
0407     std::vector<sphoton> vpp ;
0408     getPhotonAtTime(&vpp, nullptr, iprt);
0409     int num_pp = vpp.size();
0410     NP* pp = num_pp > 0 ?  NPX::ArrayFromVec<float,sphoton>(vpp, 4, 4 ) : nullptr ;
0411     return pp ;
0412 }
0413
0414 inline NP* SRecord::getSimtraceAtTime( const char* iprt ) const
0415 {
0416     std::vector<quad4> vqq ;
0417     getPhotonAtTime(nullptr, &vqq, iprt);
0418     int num_qq = vqq.size();
0419     NP* qq = num_qq > 0 ?  NPX::ArrayFromVec<float,quad4>(vqq, 4, 4 ) : nullptr ;
0420     return qq ;
0421 }
0422
0423