eic-opticks/sysrap/SSim.cc

0001 #include <map>
0002 #include <csignal>
0003
0004 #include "NPFold.h"
0005 #include "scuda.h"
0006 #include "squad.h"
0007 #include "sdigest.h"
0008 #include "stree.h"
0009
0010 #include "SLOG.hh"
0011 #include "SStr.hh"
0012 #include "ssys.h"
0013 #include "spath.h"
0014 #include "sstamp.h"
0015 #include "SScene.h"
0016 #include "SBitSet.h"
0017
0018 #include "SEvt.hh"
0019 #include "SSim.hh"
0020 #include "SBnd.h"
0021 #include "SPrd.h"
0022 #include "SPMT.h"
0023
0024 const plog::Severity SSim::LEVEL = SLOG::EnvLevel("SSim", "DEBUG");
0025 const int SSim::stree_level = ssys::getenvint("SSim__stree_level", 0) ;
0026
0027 SSim* SSim::INSTANCE = nullptr ;
0028
0029 SSim* SSim::Get(){ return INSTANCE ; }
0030
0031
0032 SSim* SSim::CreateOrReuse()
0033 {
0034     return INSTANCE ? INSTANCE : Create() ;
0035 }
0036
0037 void SSim::AddExtraSubfold(const char* k, const char* _dir) // static
0038 {
0039     const char* dir = spath::Resolve(_dir);
0040     LOG(LEVEL)
0041         << " k " << k
0042         << " _dir " << _dir
0043         << " dir " << dir
0044          ;
0045     if(NPFold::Exists(dir))
0046     {
0047         NPFold* fold = NPFold::Load(dir) ;
0048         LOG(LEVEL) << " fold " << ( fold ? "YES" : "NO " ) ;
0049         AddExtraSubfold(k, fold );
0050     }
0051     else
0052     {
0053         LOG(LEVEL) << " DOESNT EXIST : SKIP " ;
0054     }
0055 }
0056 void SSim::AddExtraSubfold(const char* k, NPFold* f) // static
0057 {
0058     SSim* ss = CreateOrReuse();
0059     LOG_IF(error, ss == nullptr ) << " SSim::INSTANCE not instanciated yet " ;
0060     if(ss == nullptr) return ;
0061     ss->add_extra_subfold(k, f);
0062 }
0063
0064
0065 int SSim::Compare( const SSim* a , const SSim* b )
0066 {
0067     return ( a && b ) ? NPFold::Compare(a->top, b->top) : -1 ;
0068 }
0069
0070 std::string SSim::DescCompare( const SSim* a , const SSim* b )
0071 {
0072     std::stringstream ss ;
0073     ss << "SSim::DescCompare"
0074        << " a " << ( a ? "Y" : "N" )
0075        << " b " << ( b ? "Y" : "N" )
0076        << std::endl
0077        << ( ( a && b ) ? NPFold::DescCompare( a->top, b->top ) : "-" )
0078        << std::endl
0079        ;
0080     std::string s = ss.str();
0081     return s ;
0082 }
0083
0084
0085 SSim* SSim::Create()
0086 {
0087     LOG_IF(fatal, INSTANCE) << "replacing SSim::INSTANCE" ;
0088     new SSim ;
0089     return INSTANCE ;
0090 }
0091
0092
0093 /**
0094 SSim::Load from persisted geometry  : used for testing
0095 -------------------------------------------------------
0096
0097 **/
0098
0099
0100 SSim* SSim::Load(){ return Load_(nullptr) ; }
0101
0102 SSim* SSim::Load_(const char* base_)
0103 {
0104     LOG(LEVEL) << "[" ;
0105     const char* ssb = spath::Resolve("$CFBaseFromGEOM/CSGFoundry") ;
0106     const char* base = spath::Resolve(base_ ? base_ : ssb );
0107     LOG(LEVEL)
0108        << " ssb [" << ( ssb ? ssb : "-" ) << "]"
0109        << " base_ [" << ( base_ ? base_ : "-" ) << "]"
0110        << " base [" << ( base ? base : "-" ) << "]"
0111        ;
0112
0113     SSim* sim = new SSim ;
0114     sim->load(base);    // reldir defaults to "SSim"
0115
0116     LOG(LEVEL) << "]" ;
0117     return sim ;
0118 }
0119
0120
0121 /**
0122 SSim::Load
0123 ------------
0124
0125 Default reldir is "SSim" so the base directory is for example::
0126
0127     $HOME/.opticks/GEOM/$GEOM/CSGFoundry
0128     /tmp/GEOM/$GEOM/CSGFoundry
0129
0130
0131 **/
0132
0133 SSim* SSim::Load(const char* base, const char* reldir)
0134 {
0135     LOG_IF(fatal, base==nullptr)
0136         << " base NULL "
0137         << " base " << ( base ? base : "-" )
0138         << " reldir " << ( reldir ? reldir : "-" )
0139         ;
0140
0141     SSim* sim = new SSim ;
0142     sim->load(base, reldir);
0143     return sim ;
0144 }
0145
0146
0147
0148
0149
0150 SSim::SSim()
0151     :
0152     relp(ssys::getenvvar("SSim__RELP", RELP_DEFAULT )),
0153     top(nullptr),
0154     extra(nullptr),
0155     tree(new stree),
0156     scene(new SScene),
0157     toploadtime(0)
0158 {
0159     init(); // just sets tree level
0160 }
0161
0162 /**
0163 SSim::init
0164 ------------
0165
0166 **/
0167
0168 void SSim::init()
0169 {
0170     INSTANCE = this ;
0171     tree->set_level(stree_level);
0172 }
0173
0174 /**
0175 SSim::AnnotateFrame
0176 ----------------------
0177
0178
0179 **/
0180
0181
0182 void SSim::AnnotateFrame( sframe& fr, const SBitSet* elv, const char* caller ) // static
0183 {
0184     const stree* tree = INSTANCE ? INSTANCE->tree : nullptr ;
0185     const char* tree_digest = tree->get_tree_digest() ;
0186
0187     std::vector<unsigned char> extra ;
0188     std::string _dynamic0 = tree->make_tree_digest( extra );
0189     assert( strcmp( _dynamic0.c_str(), tree_digest ) == 0 ); // with empty extra should get same as standard tree digest
0190
0191     if(elv) elv->serialize(extra);
0192     std::string _dynamic1 = tree->make_tree_digest( extra );
0193     const char* dynamic = _dynamic1.c_str();
0194
0195     fr.set_tree( tree_digest );
0196     fr.set_dynamic( dynamic );
0197
0198     LOG(info)
0199        << " caller " << ( caller ? caller : "-" )
0200        << " tree " << ( tree ? "YES" : "NO " )
0201        << " elv " << ( elv ? "YES" : "NO " )
0202        << " extra.size " << extra.size()
0203        << " tree_digest " << ( tree_digest ? tree_digest : "-" )
0204        << " dynamic " << ( dynamic ? dynamic : "-" )
0205        ;
0206 }
0207
0208
0209
0210 stree* SSim::get_tree() const { return tree ; }
0211 SScene* SSim::get_scene() const { return scene ; }
0212
0213 // TRANSITIONAL WHILST HAVING DIFFICULTIES WITH FULL GEOM CONVERSION
0214 void SSim::set_override_scene(SScene* _scene){ scene = _scene ; }
0215
0216
0217
0218 /**
0219 SSim::initSceneFromTree
0220 ------------------------
0221
0222 This needs to be invoked after the tree has been populated by U4Tree
0223
0224 **/
0225
0226 void SSim::initSceneFromTree()
0227 {
0228     scene->initFromTree(tree);
0229 }
0230
0231
0232
0233
0234
0235 /**
0236 SSim::lookup_mtline
0237 ---------------------
0238
0239 Lookup matline for bnd texture or array access
0240 from an original Geant4 material creation index
0241 as obtained by G4Material::GetIndex
0242
0243 NB this original mtindex is NOT GENERALLY THE SAME
0244 as the Opticks material index.
0245
0246 **/
0247
0248 int SSim::lookup_mtline( int mtindex ) const
0249 {
0250     return tree->lookup_mtline(mtindex);
0251 }
0252
0253 std::string SSim::desc_mt() const
0254 {
0255     return tree->desc_mt() ;
0256 }
0257
0258
0259
0260
0261
0262
0263 /**
0264 SSim::desc
0265 ------------
0266
0267 This and the following accessors require serialization to populate top
0268
0269 **/
0270
0271 std::string SSim::desc() const { return top ? top->desc() : "-" ; }
0272 std::string SSim::brief() const { return top ? top->brief() : "-" ; }
0273
0274 const NP* SSim::get(const char* k) const
0275 {
0276     assert( top );
0277     const NPFold* f = top->find_subfold( relp );
0278     if( f == nullptr ) std::cerr
0279         << "SSim::get"
0280         << " relp[" << ( relp ? relp : "-" ) << "]"
0281         << " k[" << ( k ? k : "-" ) << "]"
0282         << " f null "
0283         << std::endl
0284         ;
0285
0286     return f ? f->get(k) : nullptr ;
0287 }
0288 void SSim::set(const char* k, const NP* a)
0289 {
0290     assert( top );
0291     NPFold* f = top->find_subfold_( relp );
0292     f->set( k, a );
0293 }
0294
0295
0296 const NP* SSim::get_bnd() const {  return get(snam::BND);  }
0297
0298 /**
0299 SSim::getBndName
0300 -------------------
0301
0302 Return the bnd name for boundary index bidx using the
0303 metadata names list associated with the bnd.npy array.
0304
0305 **/
0306
0307 const char* SSim::getBndName(unsigned bidx) const
0308 {
0309     const NP* bnd = get_bnd();
0310     bool valid = bnd && bidx < bnd->names.size() ;
0311     if(!valid) return nullptr ;
0312     const std::string& name = bnd->names[bidx] ;
0313     return name.c_str()  ;  // no need for strdup as it lives in NP vector
0314 }
0315 int SSim::getBndIndex(const char* bname) const
0316 {
0317     unsigned count = 0 ;
0318     const NP* bnd = get_bnd();
0319     int bidx = bnd->get_name_index(bname, count );
0320     bool bname_found = count == 1 && bidx > -1  ;
0321
0322     LOG_IF(fatal, !bname_found)
0323         << " bname " << bname
0324         << " bidx " << bidx
0325         << " count " << count
0326         << " bname_found " << bname_found
0327         ;
0328
0329     assert( bname_found );
0330     return bidx ;
0331 }
0332
0333 const SBnd* SSim::get_sbnd() const
0334 {
0335     const NP* bnd = get_bnd();
0336     return bnd ? new SBnd(bnd) : nullptr  ;
0337 }
0338 const SPrd* SSim::get_sprd() const
0339 {
0340     const NP* bnd = get_bnd();
0341     return bnd ? new SPrd(bnd) : nullptr  ;
0342 }
0343
0344 /**
0345 SSim::get_jpmt
0346 ---------------
0347
0348 Note that if the top fold does not have the JPMT_RELP "extra/jpmt"  subfold
0349 then this returns nullptr.
0350
0351 So that means must first call SSim::AddExtraSubfold
0352
0353 Returning a deepcopy avoids subsequent parent changing assert
0354 from NPFold::add_subfold. But causes other problems.
0355
0356 **/
0357
0358 const NPFold* SSim::get_jpmt_nocopy() const
0359 {
0360     NPFold* f = top ? top->find_subfold_(JPMT_RELP) : nullptr ;
0361     //if(f) f->set_verbose_r();
0362     LOG(LEVEL)
0363         << " (NPFold)top " << ( top ? "YES" : "NO " )
0364         << " (NPFold)f " << ( f ? "YES" : "NO " )
0365         ;
0366     return f ;
0367 }
0368
0369 const NPFold* SSim::get_jpmt() const
0370 {
0371     const NPFold* f = get_jpmt_nocopy();
0372     const NPFold* fc  = f ? f->deepcopy() : nullptr ;
0373
0374     LOG(LEVEL)
0375         << " (NPFold)f " << ( f ? "YES" : "NO " )
0376         << " (NPFold)fc " << ( fc ? "YES" : "NO " )
0377         ;
0378
0379     return fc ;
0380 }
0381 const SPMT* SSim::get_spmt() const
0382 {
0383     const NPFold* jpmt = get_jpmt();
0384     const SPMT* spmt = jpmt ? new SPMT(jpmt) : nullptr ;
0385
0386     LOG(LEVEL)
0387         << " (NPFold)jpmt " << ( jpmt ? "YES" : "NO " )
0388         << " (SPMT)spmt " << ( spmt ? "YES" : "NO " )
0389         ;
0390
0391     return spmt ;
0392 }
0393
0394 /**
0395 SSim::get_spmt_f
0396 ------------------
0397
0398 This is invoked from QSim::UploadComponents and
0399 used to instanciate QPMT for GPU uploading.
0400
0401 TODO: tidy up this chain of four methods by
0402 moving the details into statics within SPMT.h
0403
0404 **/
0405
0406 const NPFold* SSim::get_spmt_f() const
0407 {
0408     const SPMT* spmt = get_spmt() ;
0409     const NPFold* spmt_f = spmt ? spmt->serialize() : nullptr ;
0410
0411     LOG(LEVEL)
0412         << " (SPMT)spmt " << ( spmt ? "YES" : "NO " )
0413         << " (NPFold)spmt_f " << ( spmt_f ? "YES" : "NO " )
0414         ;
0415
0416     return spmt_f ;
0417 }
0418
0419
0420
0421 void SSim::add_extra_subfold(const char* k, NPFold* f )
0422 {
0423     assert(k);
0424     LOG_IF(LEVEL, f == nullptr) << "k:" << k  << " f null " ;
0425     if(f == nullptr) return ;
0426
0427     if( extra == nullptr ) extra = new NPFold ;
0428     extra->add_subfold(k,f);
0429 }
0430
0431 /**
0432 SSim::AddMultiFilm
0433
0434 */
0435
0436 void SSim::AddMultiFilm(const char*k, const NP* f){
0437
0438
0439     SSim* ss = CreateOrReuse();
0440     LOG_IF(error, ss == nullptr ) << " SSim::INSTANCE not instanciated yet " ;
0441     if(ss == nullptr) return ;
0442     ss->set_extra(k, f);
0443 }
0444 void SSim::set_extra(const char* k,  const NP* f){
0445     assert(f);
0446     if( extra == nullptr ) extra = new NPFold ;
0447     extra->add(k, f);
0448
0449 }
0450
0451 const NP* SSim::get_extra(const char* k) const{
0452
0453     return extra ? extra->get(k) : nullptr ;
0454
0455 }
0456 /**
0457 SSim::save
0458 ------------
0459
0460 Canonical usage from CSGFoundry::save_ with::
0461
0462     sim->save(dir, SSim::RELDIR)
0463
0464 **/
0465
0466 void SSim::save(const char* base, const char* reldir)
0467 {
0468     if(top == nullptr) serialize() ;
0469     LOG_IF(fatal, top == nullptr) << " top null : MUST serialize before save, serialize failed ? " ;
0470     assert( top != nullptr ) ;
0471
0472     const char* dir = spath::Resolve(base, reldir) ;   // default reldir "SSim"
0473     top->save(dir);
0474
0475     tree->save_desc(dir, stree::RELDIR );  // implicit "desc" last element
0476 }
0477
0478
0479
0480 /**
0481 SSim::load
0482 ------------
0483
0484 **/
0485
0486 void SSim::load(const char* base, const char* reldir)
0487 {
0488     const char* dir = spath::Resolve(base, reldir) ;
0489     load_(dir);
0490 }
0491
0492 void SSim::load_(const char* dir)
0493 {
0494     LOG(LEVEL) << "[" ;
0495     LOG_IF(fatal, top != nullptr)  << " top is NOT nullptr : cannot SSim::load into pre-serialized instance " ;
0496     top = new NPFold ;
0497
0498     LOG(LEVEL) << "[ top.load [" << dir << "]" ;
0499
0500     int64_t t0 = sstamp::Now();
0501     top->load(dir) ;
0502     toploadtime = sstamp::Now() - t0 ;
0503
0504     LOG(LEVEL) << "] top.load [" << dir << "] toploadtime/1e6 " << std::fixed << std::setw(9) << std::setprecision(6) << toploadtime/1e6 ;
0505
0506     NPFold* f_tree = top->get_subfold( stree::RELDIR ) ;
0507     tree->import_( f_tree );
0508
0509
0510
0511     NPFold* f_scene = top->get_subfold( SScene::RELDIR ) ;
0512     scene->import_( f_scene );
0513
0514
0515     std::stringstream ss ;
0516     ss << dir << "/" << stree::RELDIR ;
0517     std::string treedir = ss.str();
0518     tree->loaddir = strdup(treedir.c_str());
0519
0520
0521     afterLoadOrCreate();
0522
0523     LOG(LEVEL) << "]" ;
0524 }
0525
0526
0527 /**
0528 SSim::afterLoadOrCreate
0529 ------------------------
0530
0531 Trying to progress with FRAME_TRANSITION by replacing the old
0532 CSGFoundry::AfterLoadOrCreate with SSim::afterLoadOrCreate
0533
0534 **/
0535
0536 void SSim::afterLoadOrCreate()
0537 {
0538     SEvt::CreateOrReuse() ;   // creates 1/2 SEvt depending on OPTICKS_INTEGRATION_MODE
0539     assert(tree);
0540
0541
0542     sfr fr = tree->get_frame_moi() ;
0543     SEvt::SetFr(fr);
0544 }
0545
0546
0547 /**
0548 SSim::serialize
0549 -----------------
0550
0551 NPFold layout::
0552
0553    SSim/extra
0554    SSim/stree
0555
0556
0557 This is invoked by::
0558
0559     CSGOptiX::InitSim
0560     SSim::save
0561
0562 Q: Is this needed other than for saving ? How so ?
0563
0564
0565 **/
0566
0567
0568 void SSim::serialize()
0569 {
0570     bool has_top = hasTop();
0571
0572     LOG(LEVEL) << "[" ;
0573     LOG_IF(fatal, has_top )  << " has_top : cannot serialize twice : DONT SERIALIZE AFTER LOADING SSim " ;
0574     assert( !has_top );
0575
0576     top = new NPFold ;
0577     NPFold* f_tree = tree->serialize() ;
0578     top->add_subfold( stree::RELDIR, f_tree );
0579
0580     NPFold* f_scene = scene->serialize() ;
0581     top->add_subfold( SScene::RELDIR, f_scene );
0582
0583     if( extra )
0584     {
0585         top->add_subfold( EXTRA, extra );
0586     }
0587
0588     LOG(LEVEL) << "]" ;
0589 }
0590
0591 bool SSim::hasTop() const
0592 {
0593    return top != nullptr ;
0594 }
0595
0596
0597
0598
0599
0600 template<typename ... Args>
0601 void SSim::addFake( Args ... args )
0602 {
0603     std::vector<std::string> specs = {args...};
0604     LOG(LEVEL) << "specs.size " << specs.size()  ;
0605     addFake_(specs);
0606
0607 }
0608 template void SSim::addFake( const char* );
0609 template void SSim::addFake( const char*, const char* );
0610 template void SSim::addFake( const char*, const char*, const char* );
0611
0612 /**
0613 SSim::addFake_
0614 ----------------
0615
0616 Fabricates boundaries and appends them to the bnd and optical arrays
0617
0618 **/
0619
0620 void SSim::addFake_( const std::vector<std::string>& specs )
0621 {
0622     bool has_optical = hasOptical();
0623     LOG_IF(fatal, !has_optical) << " optical+bnd are required " ;
0624     assert(has_optical);
0625
0626     const NP* optical = get(snam::OPTICAL);
0627     const NP* bnd     = get(snam::BND);
0628
0629     NP* opticalplus = nullptr ;
0630     NP* bndplus = nullptr ;
0631
0632     Add( &opticalplus, &bndplus, optical, bnd, specs );
0633
0634     //NOTE: are leaking the old ones
0635     set(snam::OPTICAL, opticalplus);
0636     set(snam::BND,     bndplus);
0637 }
0638
0639
0640
0641 /**
0642 SSim::Add
0643 -----------
0644
0645 Coordinates addition of boundaries to the optical and bnd buffers using the boundary string
0646 specification.
0647
0648 **/
0649
0650 void SSim::Add(
0651     NP** opticalplus,
0652     NP** bndplus,
0653     const NP* optical,
0654     const NP* bnd,
0655     const std::vector<std::string>& specs ) // static
0656 {
0657     *opticalplus = AddOptical(optical, bnd->names, specs );
0658     *bndplus = AddBoundary( bnd, specs );
0659 }
0660
0661 /**
0662 SSim::AddOptical
0663 ------------------
0664
0665 Used from SSim::Add in coordination with SSim::AddBoundary.
0666 Using this alone would break optical:bnd consistency.
0667
0668 optical buffer has 4 uint for each species and 4 species for each boundary
0669
0670 Water/Steel_surface/Steel_surface/Steel
0671   19    0    0    0
0672   21    0    3   20
0673   21    0    3   20
0674    4    0    0    0
0675
0676 The .x is the 1-based material or surface index with 0 signifying none
0677 which shoild only ever happen for surfaces.
0678
0679 **/
0680
0681 NP* SSim::AddOptical(
0682     const NP* optical,
0683     const std::vector<std::string>& bnames,
0684     const std::vector<std::string>& specs )
0685 {
0686     unsigned ndim = optical->shape.size() ;
0687     bool ndim_expect = ndim == 2 ;
0688     assert( ndim_expect );
0689     if(!ndim_expect) std::raise(SIGINT);
0690
0691     unsigned num_bnd = bnames.size() ;
0692     unsigned num_add = specs.size()  ;
0693     unsigned ni = optical->shape[0] ;
0694     unsigned nj = optical->shape[1] ;
0695
0696     bool num_bnd_expect = 4*num_bnd == ni ;
0697     bool nj_expect = nj == 4  ;
0698     bool optical_expect = optical->ebyte == 4 && optical->uifc == 'u' ;
0699
0700     assert( num_bnd_expect );
0701     assert( nj_expect );
0702     assert( optical_expect ) ;
0703
0704     if(!num_bnd_expect) std::raise(SIGINT);
0705     if(!nj_expect) std::raise(SIGINT);
0706     if(!optical_expect) std::raise(SIGINT);
0707
0708     unsigned item_bytes = optical->ebyte*optical->itemsize_(0);
0709     bool item_bytes_expect = item_bytes == 16u ;
0710     assert( item_bytes_expect );
0711     if(!item_bytes_expect ) std::raise(SIGINT);
0712
0713     NP* opticalplus = new NP(optical->dtype);
0714     std::vector<NP::INT> opticalplus_shape(optical->shape);
0715     opticalplus_shape[0] += 4*num_add ;
0716     opticalplus->set_shape(opticalplus_shape) ;
0717
0718     unsigned offset = 0 ;
0719     unsigned optical_arr_bytes = optical->arr_bytes() ;
0720     memcpy( opticalplus->bytes() + offset, optical->bytes(), optical_arr_bytes );
0721     offset += optical_arr_bytes ;
0722
0723     uint4 item = make_uint4( 0u, 0u, 0u, 0u );
0724
0725     for(unsigned b=0 ; b < num_add ; b++)
0726     {
0727         const char* spec = SStr::Trim(specs[b].c_str());
0728         std::vector<std::string> elem ;
0729         SStr::Split(spec, '/', elem );
0730
0731         bool four_elem = elem.size() == 4 ;
0732         LOG_IF(fatal, four_elem == false) << " expecting four elem spec [" << spec << "] elem.size " << elem.size() ;
0733         assert(four_elem);
0734
0735         for(unsigned s=0 ; s < 4 ; s++ )
0736         {
0737             const char* qname = elem[s].c_str();
0738             int i, j ;
0739             bool found = SBnd::FindName(i, j, qname, bnames );
0740
0741             NP::INT idx = i*4 + j ;
0742
0743             const char* ibytes = nullptr ;
0744             NP::INT num_bytes = 0 ;
0745
0746             if(found)
0747             {
0748                 optical->itembytes_( &ibytes, num_bytes, idx );
0749             }
0750             else if(strstr(qname, "perfect"))
0751             {
0752                 assert( s == 1 || s == 2 );  // only expecting "perfect" surfaces not materials
0753                  //
0754                 // NB: when found==false (i,j) will be stale or undefined SO DO NOT USE THEM HERE
0755                 //
0756                 // NB: the only way this item is currently used is checking of item.x (aka s.optical.x)  > 0
0757                 //     to indicate that propagate_at_surface should be used and not propagate_at_boundary
0758                 //     while the value for .x has traditionally been a 1-based surface index
0759                 //     that index is at qsim.h level just metadata : it is never used to lookup anything
0760                 //
0761                 // TODO: mint a new index to use for added surfaces, rather than here just using 99u
0762                 //
0763                 item.x = 99u ;
0764                 item.y = 99u ;
0765                 item.z = 99u ;
0766                 item.w = 99u ;
0767
0768                 ibytes = (const char*)&item;
0769                 num_bytes = sizeof(uint4);
0770             }
0771             else
0772             {
0773                 LOG(error) << "SBin::FindName failed to find qname [" << qname << "] from within the bnames.size " << bnames.size() ;
0774                 for(unsigned z=0 ; z < bnames.size() ; z++) LOG(error) << " z " << z << " bnames[z] " << bnames[z] ;
0775                 assert( 0 );
0776             }
0777             assert( ibytes != nullptr );
0778             assert( num_bytes == item_bytes );
0779             memcpy( opticalplus->bytes() + offset,  ibytes, item_bytes );
0780             offset += item_bytes ;
0781         }
0782     }
0783     return opticalplus ;
0784 }
0785
0786
0787
0788
0789 /**
0790 SSim::AddBoundary
0791 ------------------------
0792
0793 Canonically invoked from SSim::Add in coordination with SSim::AddOptical to maintain consistency.
0794
0795 Creates new array containing the src array with extra boundaries constructed
0796 from materials and surfaces already present in the src array as configured by the
0797 specs argument.
0798
0799 **/
0800
0801 NP* SSim::AddBoundary( const NP* dsrc, const std::vector<std::string>& specs ) // static
0802 {
0803     const NP* src = NP::MakeNarrowIfWide(dsrc) ;
0804
0805     unsigned ndim = src->shape.size() ;
0806     bool ndim_expect = ndim == 5  ;
0807     assert( ndim_expect );
0808     if(!ndim_expect) std::raise(SIGINT);
0809
0810     unsigned ni = src->shape[0] ;
0811     unsigned nj = src->shape[1] ;
0812     unsigned nk = src->shape[2] ;
0813     unsigned nl = src->shape[3] ;
0814     unsigned nm = src->shape[4] ;
0815
0816     LOG(LEVEL)
0817         << " src.ebyte " << src->ebyte
0818         << " src.desc " << src->desc()
0819         ;
0820
0821     bool src_expect = src->ebyte == 4 && ni > 0  && nj == 4 ;
0822     // expecting 2nd dimension to be 4: omat/osur/isur/imat
0823     assert( src_expect ) ;
0824     if(!src_expect) std::raise(SIGINT);
0825
0826
0827     unsigned src_bytes = src->arr_bytes() ;
0828     unsigned bnd_bytes = src->ebyte*src->itemsize_(0) ;
0829     unsigned sub_bytes = src->ebyte*src->itemsize_(0,0) ;
0830
0831     bool bnd_bytes_expect =  bnd_bytes == 4*sub_bytes ;
0832     assert( bnd_bytes_expect );
0833     if(!bnd_bytes_expect) std::raise(SIGINT);
0834
0835     NP* dst = new NP(src->dtype);
0836
0837     std::vector<NP::INT> dst_shape(src->shape);
0838     dst_shape[0] += specs.size() ;
0839     dst->set_shape(dst_shape) ;
0840
0841     std::vector<std::string> names ;
0842     src->get_names(names);
0843
0844     std::vector<std::string> dst_names(names);
0845     LOG(LEVEL)
0846         << " dst_names.size before " << dst_names.size()
0847         << " specs.size " << specs.size()
0848         ;
0849
0850     unsigned offset = 0 ;
0851     memcpy( dst->bytes() + offset, src->bytes(), src_bytes );
0852     offset += src_bytes ;
0853
0854     for(unsigned b=0 ; b < specs.size() ; b++)
0855     {
0856         const char* spec = SStr::Trim(specs[b].c_str());
0857         dst_names.push_back(spec);
0858
0859         std::vector<std::string> elem ;
0860         SStr::Split(spec, '/', elem );
0861
0862         bool four_elem = elem.size() == 4 ;
0863         LOG_IF(fatal, four_elem == false) << " expecting four elem spec [" << spec << "] elem.size " << elem.size() ;
0864         assert(four_elem);
0865
0866         for(unsigned s=0 ; s < 4 ; s++ )
0867         {
0868             const char* qname = elem[s].c_str();
0869             int i, j ;
0870             bool found = SBnd::FindName(i, j, qname, names );
0871
0872             const char* ibytes = nullptr ;
0873             NP::INT num_bytes = 0 ;
0874
0875             if(found)
0876             {
0877                 src->itembytes_( &ibytes, num_bytes, i, j );
0878             }
0879             else if(strstr(qname, "perfect"))
0880             {
0881                 std::vector<float> values ;
0882                 GetPerfectValues( values, nk, nl, nm, qname );
0883                 ibytes = (const char*)values.data();
0884                 num_bytes = sizeof(float)*values.size();
0885             }
0886             else
0887             {
0888                 LOG(fatal) << " FAILED to find qname " << qname ;
0889                 assert( 0 );
0890             }
0891
0892             assert( ibytes != nullptr );
0893             assert( num_bytes == sub_bytes );
0894
0895             memcpy( dst->bytes() + offset,  ibytes, num_bytes );
0896             offset += sub_bytes ;
0897         }
0898     }
0899
0900     LOG(LEVEL) << " dst_names.size after " << dst_names.size() ;
0901
0902     dst->set_names( dst_names );
0903     dst->meta = src->meta ;    // need to pass along the domain metadata
0904
0905     std::vector<std::string> dst_names_check ;
0906     dst->get_names(dst_names_check);
0907
0908     LOG(LEVEL) << " dst_names_check.size after " << dst_names_check.size() ;
0909
0910     return dst ;
0911 }
0912
0913
0914
0915
0916 /**
0917 SSim::GetPerfectValues
0918 -------------------------
0919
0920 bnd with shape (44, 4, 2, 761, 4, )::
0921
0922    ni : boundaries
0923    nj : 0:omat/1:osur/2:isur/3:imat
0924    nk : 0 or 1 property group
0925    nl : wavelengths
0926    nm : payload
0927
0928 **/
0929
0930 void SSim::GetPerfectValues(
0931     std::vector<float>& values,
0932     unsigned nk, unsigned nl, unsigned nm, const char* name ) // static
0933 {
0934     LOG(LEVEL) << name << " nk " << nk << " nl " << nl << " nm " << nm ;
0935
0936     assert( nk == 2 );
0937     assert( nl > 0 );
0938     assert( nm == 4 );
0939
0940     float4 payload[2] ;
0941     if(     strstr(name, "perfectDetectSurface"))   payload[0] = make_float4(  1.f, 0.f, 0.f, 0.f );
0942     else if(strstr(name, "perfectAbsorbSurface"))   payload[0] = make_float4(  0.f, 1.f, 0.f, 0.f );
0943     else if(strstr(name, "perfectSpecularSurface")) payload[0] = make_float4(  0.f, 0.f, 1.f, 0.f );
0944     else if(strstr(name, "perfectDiffuseSurface"))  payload[0] = make_float4(  0.f, 0.f, 0.f, 1.f );
0945     else                                            payload[0] = make_float4( -1.f, -1.f, -1.f, -1.f );
0946
0947     payload[1] = make_float4( -1.f, -1.f, -1.f, -1.f );
0948
0949     values.resize( nk*nl*nm );
0950     unsigned idx = 0 ;
0951     unsigned count = 0 ;
0952     for(unsigned k=0 ; k < nk ; k++)          // over payload groups
0953     {
0954         const float4& pay = payload[k] ;
0955         for(unsigned l=0 ; l < nl ; l++)      // payload repeated over wavelength samples
0956         {
0957             for(unsigned m=0 ; m < nm ; m++)  // over payload values
0958             {
0959                  idx = k*nl*nm + l*nm + m ;
0960                  assert( idx == count );
0961                  count += 1 ;
0962                  switch(m)
0963                  {
0964                      case 0: values[idx] = pay.x ; break ;
0965                      case 1: values[idx] = pay.y ; break ;
0966                      case 2: values[idx] = pay.z ; break ;
0967                      case 3: values[idx] = pay.w ; break ;
0968                  }
0969             }
0970         }
0971     }
0972 }
0973
0974 bool SSim::hasOptical() const
0975 {
0976     const NP* optical = get(snam::OPTICAL);
0977     const NP* bnd = get(snam::BND);
0978     bool has_optical = optical != nullptr && bnd != nullptr ;
0979     return has_optical ;
0980 }
0981
0982
0983 std::string SSim::descOptical() const
0984 {
0985     const NP* optical = get(snam::OPTICAL);
0986     const NP* bnd = get(snam::BND);
0987
0988     if(optical == nullptr && bnd == nullptr) return "SSim::descOptical null" ;
0989     return DescOptical(optical, bnd);
0990 }
0991
0992 std::string SSim::DescOptical(const NP* optical, const NP* bnd )
0993 {
0994     int num_bnd = bnd->shape[0] ;
0995     int num_bnd_names = bnd->names.size() ;
0996     assert( num_bnd == num_bnd_names );
0997
0998     int num_optical = optical->shape[0] ;
0999     bool consistent = num_optical == num_bnd*4  ;
1000
1001     typedef std::map<unsigned, std::string> MUS ;
1002     MUS surf ;
1003
1004     std::stringstream ss ;
1005     ss << "SSim::DescOptical"
1006        << " optical " << optical->sstr()
1007        << " bnd " << bnd->sstr()
1008        << " num_bnd_names " << num_bnd_names
1009        << " consistent " << ( consistent ? "YES" : "NO:ERROR" )
1010        << std::endl
1011        ;
1012
1013     assert(consistent);
1014     assert( optical->shape.size() == 2 );
1015
1016     unsigned ni = optical->shape[0] ;
1017     unsigned nj = optical->shape[1] ;
1018     assert( nj == 4 );
1019
1020     const unsigned* oo = optical->cvalues<unsigned>();
1021
1022     std::vector<std::string> elem ;
1023
1024     for(unsigned i=0 ; i < ni ; i++)
1025     {
1026         unsigned b = i/4 ;
1027         unsigned ii = i % 4 ;
1028         if( ii == 0 )
1029         {
1030             elem.clear();
1031             const std::string& spec = bnd->names[b] ;
1032             SStr::Split( spec.c_str(), '/', elem );
1033             ss << std::setw(4) << b << " " << spec<< std::endl ;
1034         }
1035
1036         const std::string& name = elem[ii] ;
1037
1038         ss << std::setw(4) << i << " : " << std::setw(4) << ii << " : " ;
1039
1040         for(unsigned j=0 ; j < nj ; j++)
1041         {
1042             unsigned idx = i*nj + j ;
1043             unsigned val = oo[idx] ;
1044             ss << std::setw(4) << val << " " ;
1045
1046             if( j == 0 && val > 0 && ( ii == 1 || ii == 2)  )
1047             {
1048                 surf[val] = name ;
1049             }
1050         }
1051         ss << " " << name << std::endl ;
1052     }
1053
1054     ss << " surfaces ....... " << std::endl ;
1055     for(MUS::const_iterator it=surf.begin() ; it != surf.end() ; it++)
1056     {
1057         ss << std::setw(5) << it->first << " : " << it->second << std::endl ;
1058     }
1059
1060     std::string s = ss.str() ;
1061     return s ;
1062 }
1063
1064 std::string SSim::GetItemDigest( const NP* bnd, int i, int j, int w )
1065 {
1066     //std::string dig = SDigestNP::Item(bnd, i, j ) ;
1067     std::string dig = sdigest::Item(bnd, i, j ) ;
1068     std::string sdig = dig.substr(0, w);
1069     return sdig ;
1070 }
1071
1072 bool SSim::findName( int& i, int& j, const char* qname ) const
1073 {
1074     const NP* bnd = get_bnd();
1075     return bnd ? SBnd::FindName(i, j, qname, bnd->names) : false ;
1076 }
1077
1078
1079
1080
1081
1082