File indexing completed on 2026-04-09 07:49:55
0001
0002 """
0003 stree.py
0004 ==========
0005
0006 TODO: improve stree.py repr
0007
0008 """
0009
0010 import os, numpy as np, logging, builtins
0011 log = logging.getLogger(__name__)
0012 from opticks.ana.fold import Fold
0013 from opticks.sysrap.sfreq import sfreq
0014 from opticks.sysrap.OpticksCSG import CSG_
0015
0016
0017 class STR(str):
0018 """STR inherits from str and changes the repr to provide the str : useful for interactive ipython"""
0019 def __repr__(self):
0020 return str(self)
0021
0022
0023 class sobject(object):
0024 @classmethod
0025 def Label(cls, spc=5, pfx=10):
0026 prefix = " " * pfx
0027 spacer = " " * spc
0028 return prefix + spacer.join(cls.FIELD)
0029
0030 @classmethod
0031 def Fields(cls, bi=False):
0032 kls = cls.__name__
0033 for i, field in enumerate(cls.FIELD):
0034 setattr(cls, field, i)
0035 if bi:setattr(builtins, field, i)
0036 pass
0037
0038 @classmethod
0039 def Type(cls):
0040 cls.Fields()
0041 kls = cls.__name__
0042 print("%s.Type()" % kls )
0043 for i, field in enumerate(cls.FIELD):
0044 name = cls.DTYPE[i][0]
0045 fieldname = "%s.%s" % (kls, field)
0046 print(" %2d : %20s : %s " % (i, fieldname, name))
0047 pass
0048 print("%s.Label() : " % cls.Label() )
0049
0050 @classmethod
0051 def RecordsFromArrays(cls, a):
0052 """
0053 :param a: ndarray
0054 :return: np.recarray
0055 """
0056 ra = np.core.records.fromarrays(a.T, dtype=cls.DTYPE )
0057 return ra
0058
0059 @classmethod
0060 def Doc(cls):
0061 assert len(cls.FIELD) == len(cls.DTYPE)
0062 lines = []
0063 for i in range(len(cls.FIELD)):
0064 line = "%2d : %2s : %15s : %s " % (i, cls.FIELD[i], cls.DTYPE[i][0], cls.DTYPE[i][1] )
0065 lines.append(line)
0066 pass
0067 return STR("\n".join(lines))
0068
0069
0070
0071
0072 class sn(sobject):
0073 """
0074 sn.h CSG constituent node, assuming the default WITH_CHILD
0075 """
0076 DTYPE = [
0077 ('typecode', '<i4'),
0078 ('complement', '<i4'),
0079 ('lvid', '<i4'),
0080 ('xform', '<i4'),
0081 ('param', '<i4'),
0082 ('aabb', '<i4'),
0083 ('parent', '<i4'),
0084 ('sibdex', '<i4'),
0085 ('num_child', '<i4'),
0086 ('first_child', '<i4'),
0087 ('next_sibling', '<i4'),
0088 ('index', '<i4'),
0089 ('depth', '<i4'),
0090 ('note', '<i4'),
0091 ('coincide', '<i4'),
0092 ('label0', '<i4'),
0093 ('label1', '<i4'),
0094 ('label2', '<i4'),
0095 ('label3', '<i4'),
0096 ]
0097
0098 FIELD = "tc cm lv xf pa bb pr sx nc fc ns ix dp nt l0 l1 l2 l3".split()
0099
0100
0101
0102
0103
0104
0105
0106 class snd(sobject):
0107 """
0108 snd.hh CSG constituent node
0109 """
0110 DTYPE = [
0111 ('index', '<i4'),
0112 ('depth', '<i4'),
0113 ('sibdex', '<i4'),
0114 ('parent', '<i4'),
0115 ('num_child', '<i4'),
0116 ('first_child', '<i4'),
0117 ('next_sibling', '<i4'),
0118 ('lvid', '<i4'),
0119 ('typecode', '<i4'),
0120 ('param', '<i4'),
0121 ('aabb', '<i4'),
0122 ('xform', '<i4'),
0123 ]
0124
0125 FIELD = "ix dp sx pt nc fc sx lv tc pm bb xf".split()
0126
0127
0128 class snode(sobject):
0129 """
0130 Volume structure node
0131 """
0132 DTYPE = [
0133 ('index', '<i4'),
0134 ('depth', '<i4'),
0135 ('sibdex', '<i4'),
0136 ('parent', '<i4'),
0137 ('num_child', '<i4'),
0138 ('first_child', '<i4'),
0139 ('next_sibling', '<i4'),
0140 ('lvid', '<i4'),
0141 ('copyno', '<i4'),
0142 ('sensor_id', '<i4'),
0143 ('sensor_index', '<i4'),
0144 ('repeat_index', '<i4'),
0145 ('repeat_ordinal', '<i4'),
0146 ('boundary', '<i4'),
0147 ('sensor_name', '<i4')
0148 ]
0149
0150 FIELD = "ix dp sx pt nc fc sx lv cp se sx ri ro bd sn".split()
0151
0152 @classmethod
0153 def Desc(cls, rec):
0154 """
0155 :param rec: single rec eg obtained by nds[0]
0156 :return str: formatted description
0157 """
0158 return "snode ix:%7d dh:%2d sx:%5d pt:%7d nc:%5d fc:%7d ns:%7d lv:%3d cp:%7d se:%7d." % tuple(rec)
0159
0160 @classmethod
0161 def Brief(cls, rec):
0162 return "snode ix:%7d dh:%2s nc:%5d lv:%3d se:%7d." % (rec.index, rec.depth, rec.num_child, rec.lvid, rec.sensor_id)
0163
0164
0165 class sfactor(sobject):
0166 """
0167 Just handles the first four fields, not the |S32 digest
0168 """
0169 DTYPE = [
0170 ( 'index', '<i4'),
0171 ( 'freq', '<i4'),
0172 ( 'sensors', '<i4'),
0173 ( 'subtree', '<i4')
0174 ]
0175 FIELD = "ix fr se su".split()
0176
0177
0178 class stree(object):
0179 @classmethod
0180 def MakeTxtArray(cls, lines):
0181 maxlen = max(list(map(len, lines)))
0182 ta = np.array( lines, dtype="|S%d" % maxlen )
0183 return ta
0184
0185 def __init__(self, f, symbol="st"):
0186
0187
0188 sff = f.subs_freq if getattr( f, "subs_freq", None ) != None else None
0189
0190 sf = None if sff is None else sfreq(sff)
0191 nds = None if f.nds is None else snode.RecordsFromArrays(f.nds)
0192 rem = None if f.rem is None else snode.RecordsFromArrays(f.rem)
0193 csg = None if f._csg is None else sn.RecordsFromArrays(f._csg.sn)
0194
0195
0196 factor = None if f.factor is None else sfactor.RecordsFromArrays(f.factor[:,:4])
0197
0198 soname_ = f.soname_names
0199
0200 self.sf = sf
0201 self.f = f
0202 self.symbol = symbol
0203 self.nds = nds
0204 self.rem = rem
0205 self.csg = csg
0206 self.factor = factor
0207 self.raw_subs = None
0208 self.soname_ = soname_
0209
0210 def find_lvid(self, q_soname, starting=True):
0211 """
0212 Pedestrian way to find string in a list
0213 """
0214 lines = self.f.soname.lines
0215 lvid = -1
0216 for i,soname in enumerate(lines):
0217 match = soname.startswith(q_soname) if starting else soname == q_soname
0218 if match:
0219 lvid = i
0220 break
0221 pass
0222 pass
0223 return lvid
0224
0225 def find_lvid_(self, q_soname_, starting=True ):
0226 """
0227 find array indices starting or exactly matching q_soname
0228 unlike the above this returns an array with all matching indices
0229 """
0230 q_soname = q_soname_.encode() if type(q_soname_) is str else q_soname_
0231 assert type(q_soname) is bytes
0232 if starting:
0233 ii = np.where( np.core.defchararray.find(self.soname_, q_soname) == 0 )[0]
0234 else:
0235 ii = np.where( self.soname_ == q_soname )[0]
0236 pass
0237 return ii
0238
0239 def find_lvid_nodes(self, arg):
0240 """
0241 :param arg: integer lvid or string soname used to obtain lvid
0242 :return: array of indices of all nodes with the specified lvid index or solid name
0243 """
0244 if type(arg) in [int, np.int32, np.int64, np.uint32, np.uint64]:
0245 lvid = int(arg)
0246 elif type(arg) is str or type(arg) is bytes:
0247 ii = self.find_lvid_(arg)
0248 assert len(ii) > 0
0249 lvid = int(ii[0])
0250 if len(ii) > 1:
0251 log.warning("multiple lvid %s correspond to arg:%s using first %s " % (str(ii), arg, lvid))
0252 pass
0253 else:
0254 print("arg %s unhandled %s " % (arg, type(arg)))
0255 assert(0)
0256 pass
0257 return np.where( self.nds.lvid == lvid )[0]
0258
0259 def find_lvid_node( self, arg, ordinal=0 ):
0260 nn = self.find_lvid_nodes(arg)
0261 return nn[ordinal]
0262
0263
0264 def get_children(self, nidx, prepend_arg=False):
0265 """
0266 This is a direct translation of the C++ stree::get_children
0267 which will inevitably be slow in python.
0268 HMM: how to do this in a more numpy way ?
0269 """
0270 children = []
0271 if prepend_arg:
0272 children.append(nidx)
0273 pass
0274 nd = self.nds[nidx]
0275 assert nd.index == nidx
0276 ch = nd.first_child
0277 while ch > -1:
0278 child = self.nds[ch]
0279 assert child.parent == nd.index
0280 children.append(child.index)
0281 ch = child.next_sibling
0282 pass
0283 return children
0284
0285 def get_progeny_r(self, nidx, progeny):
0286 """
0287 :param nidx: node index
0288 :param progeny: in/out list of node indices
0289
0290 Recursively collects node indices of progeny beneath nidx,
0291 not including nidx itself.
0292 """
0293 children = self.get_children(nidx)
0294 progeny.extend(children)
0295 for nix in children:
0296 self.get_progeny_r(nix, progeny)
0297 pass
0298
0299 def get_progeny(self, nidx):
0300 progeny = []
0301 self.get_progeny_r(nidx, progeny)
0302 return np.array(progeny)
0303
0304 def get_soname(self, nidx):
0305 lvid = self.get_lvid(nidx)
0306 return self.f.soname[lvid] if lvid > -1 and lvid < len(self.f.soname) else None ;
0307
0308 def get_sub(self, nidx):
0309 return self.f.subs[nidx] if nidx > -1 else None ;
0310
0311 def get_depth(self, nidx):
0312 return self.nds[nidx].depth if nidx > -1 else -1 ;
0313
0314 def get_parent(self, nidx):
0315 return self.nds[nidx].parent if nidx > -1 else -1 ;
0316
0317 def get_lvid(self, nidx):
0318 return self.nds[nidx].lvid if nidx > -1 else -1 ;
0319
0320 def get_transform(self, nidx):
0321 return self.f.trs[nidx] if nidx > -1 and nidx < len(self.f.trs) else None
0322
0323 def get_ancestors(self, nidx):
0324 ancestors = []
0325 parent = self.get_parent(nidx)
0326 while parent > -1:
0327 ancestors.append(parent)
0328 parent = self.get_parent(parent)
0329 pass
0330 ancestors.reverse()
0331 return ancestors
0332
0333 def __repr__(self):
0334 return repr(self.f)
0335 def __str__(self):
0336 return str(self.f)
0337
0338 @classmethod
0339 def DepthSpacer(cls, dep, depmax=15):
0340 fmt = "|S%d" % depmax
0341 spc = np.zeros( (depmax,), dtype=np.int8 )
0342 spc[:] = ord(" ")
0343 spc[dep] = ord("+")
0344 return spc.view(fmt)[0].decode()
0345
0346 def desc_node(self, nidx, brief=False):
0347 return self.desc_node_(nidx, self.sf, brief=brief)
0348
0349 def desc_nodes(self, nodes, brief=False):
0350 return "\n".join([self.desc_node(nix, brief=brief) for nix in nodes])
0351
0352 def desc_node_(self, nidx, sf, brief=False):
0353
0354 nd = self.nds[nidx]
0355 dep = self.get_depth(nidx)
0356 spc = self.DepthSpacer(dep)
0357 ndd = snode.Brief(nd) if brief else snode.Desc(nd)
0358 sub = "skip-sub"
0359 sfd = "skip-sfd"
0360 son = "skip-son"
0361 return " ".join([spc,ndd,sfd,son])
0362
0363 def desc_nodes_(self, nodes, sf):
0364 return "\n".join( [self.desc_node_(nix, sf) for nix in nodes])
0365
0366 def make_freq(self, nodes):
0367 assert type(nodes) is np.ndarray
0368 if self.raw_subs is None:
0369 path = os.path.join(self.f.base, "subs.txt")
0370 self.raw_subs = np.loadtxt( path, dtype="|S32")
0371 pass
0372 ssub = self.raw_subs[nodes]
0373 ssf = sfreq.CreateFromArray(ssub)
0374 return ssf
0375
0376
0377 def desc_boundary_stats(self):
0378 """
0379 Counts of volumes with each boundary
0380
0381 desc_boundary_stats
0382 u_bd, n_bd = np.unique( st.nds.boundary, return_counts=True )
0383 0 : 0 : 1 : Galactic///Galactic
0384 1 : 1 : 3 : Galactic///Rock
0385 2 : 2 : 1 : Rock///Galactic
0386 3 : 3 : 1 : Rock//Implicit_RINDEX_NoRINDEX_pDomeAir_pDomeRock/Air
0387 4 : 4 : 1 : Rock///Rock
0388 5 : 5 : 1 : Rock//Implicit_RINDEX_NoRINDEX_pExpHall_pExpRockBox/Air
0389 6 : 6 : 1 : Air/Implicit_RINDEX_NoRINDEX_pExpHall_pPoolCover//Steel
0390 ...
0391 108 : 108 : 4997 : Pyrex/HamamatsuR12860_PMT_20inch_photocathode_mirror_logsurf/HamamatsuR12860_PMT_20inch_photocathode_mirror_logsurf/Vacuum
0392 109 : 109 : 4997 : Vacuum/HamamatsuR12860_PMT_20inch_dynode_plate_opsurface//Steel
0393 110 : 110 : 4997 : Vacuum/HamamatsuR12860_PMT_20inch_outer_edge_opsurface//Steel
0394 111 : 111 : 4997 : Vacuum/HamamatsuR12860_PMT_20inch_inner_edge_opsurface//Steel
0395
0396
0397 """
0398 lines = []
0399 lines.append("desc_boundary_stats")
0400 st = self
0401 lines.append("u_bd, n_bd = np.unique( st.nds.boundary, return_counts=True )")
0402 u_bd, n_bd = np.unique( st.nds.boundary, return_counts=True )
0403
0404 for i in range(len(u_bd)):
0405 u = u_bd[i]
0406 n = n_bd[i]
0407 line = " %3d : %4d : %6d : %s " % (i, u, n, st.f.standard.bd_names[u] )
0408 lines.append(line)
0409 pass
0410 return STR("\n".join(lines))
0411
0412
0413 def desc_remainder(self):
0414 lines = []
0415 lines.append("desc_remainder")
0416 st = self
0417 u_lv, n_lv = np.unique( st.rem.lvid, return_counts=True )
0418 lines.append("u_lv, n_lv = np.unique( st.rem.lvid, return_counts=True )")
0419
0420 for i in range(len(u_lv)):
0421 u = u_lv[i]
0422 n = n_lv[i]
0423 line = " %3d : %4d : %6d : %s " % (i, u, n, st.f.soname[u] )
0424 lines.append(line)
0425 pass
0426 return "\n".join(lines)
0427
0428
0429 def desc_csg(self, lvid=105):
0430 """
0431 In [8]: snd.Label(3,8), n[n[:,snd.lv] == 105]
0432 Out[8]:
0433 (' ix dp sx pt nc fc sx lv tc pm bb xf',
0434 array([[483, 4, 0, 484, 0, -1, -1, 105, 105, 294, 294, -1],
0435 [484, 3, 0, 486, 1, 483, 485, 105, 11, -1, -1, -1],
0436 [485, 3, 1, 486, 0, -1, -1, 105, 103, 295, 295, 183],
0437 [486, 2, 0, 488, 2, 484, 487, 105, 1, -1, -1, -1],
0438 [487, 2, 1, 488, 0, -1, -1, 105, 105, 296, 296, 184],
0439 [488, 1, 0, 495, 2, 486, 494, 105, 1, -1, -1, -1],
0440 [489, 4, 0, 490, 0, -1, -1, 105, 105, 297, 297, -1],
0441 [490, 3, 0, 492, 1, 489, 491, 105, 11, -1, -1, -1],
0442 [491, 3, 1, 492, 0, -1, -1, 105, 103, 298, 298, 186],
0443 [492, 2, 0, 494, 2, 490, 493, 105, 1, -1, -1, -1],
0444 [493, 2, 1, 494, 0, -1, -1, 105, 105, 299, 299, 187],
0445 [494, 1, 1, 495, 2, 492, -1, 105, 1, -1, -1, 188],
0446 [495, 0, -1, -1, 2, 488, -1, 105, 3, -1, -1, -1]], dtype=int32))
0447
0448 In [9]: st.f.soname[105]
0449 Out[9]: 'HamamatsuR12860Tail0x61b5500'
0450 """
0451 csg = self.get_csg(lvid)
0452 lines = []
0453 lines.append("desc_csg lvid:%d st.f.soname[%d]:%s " % (lvid,lvid,self.get_lv_soname(lvid)))
0454 lines.append(snd.Label(3,8))
0455 lines.append("%s" % repr(csg))
0456 return "\n".join(lines)
0457
0458 def get_csg(self, lvid):
0459 st = self
0460 n = st.f._csg.sn
0461 return n[n[:,sn.lv] == lvid]
0462
0463 def get_csg_typecode(self, lvid):
0464 n = self.get_csg(lvid)
0465 return n[:,sn.tc]
0466
0467
0468 def get_numSolid(self):
0469 return 1+len(self.factor)
0470
0471 def get_numPrim(self, ridx):
0472 return len(self.rem) if ridx == 0 else self.factor[ridx-1].subtree
0473
0474 def get_lvid(self, ridx):
0475 """
0476 :param ridx:
0477 :return lvid: array of meshIdx
0478
0479 In [3]: ridx = 1 ; st.nds.lvid[st.nds.repeat_index == ridx][:st.factor[ridx-1].subtree]
0480 Out[3]: array([133, 131, 129, 130, 132], dtype=int32)
0481
0482 In [4]: ridx = 2 ; st.nds.lvid[st.nds.repeat_index == ridx][:st.factor[ridx-1].subtree]
0483 Out[4]: array([128, 118, 119, 127, 126, 120, 125, 121, 122, 123, 124], dtype=int32)
0484
0485 In [9]: st.nds.lvid[st.nds.repeat_index == 0].shape
0486 Out[9]: (3089,)
0487 """
0488 st = self
0489 return st.nds.lvid[st.nds.repeat_index == ridx] if ridx == 0 else st.nds.lvid[st.nds.repeat_index == ridx][:st.factor[ridx-1].subtree]
0490
0491 def get_lv_soname(self, lv):
0492 return str(self.soname_[lv])
0493
0494 def descSolid(self, ridx, detail=False):
0495 """
0496 cf with CSGFoundry.descSolid
0497 """
0498 numPrim = self.get_numPrim(ridx)
0499 lvid = self.get_lvid(ridx)
0500 u_lvid, n_lvid = np.unique(lvid, return_counts=True )
0501 n_lvid_one = np.all( n_lvid == 1 )
0502 p_lvid = lvid if n_lvid_one else u_lvid
0503 pass
0504 lines = []
0505 lines.append("stree.descSolid ridx %3d numPrim %5d lvid %s n_lvid_one %d" % (ridx, numPrim, str(lvid), n_lvid_one))
0506 if detail:
0507 lines.append("")
0508 for pass_ in [0,1]:
0509 for i in range(len(p_lvid)):
0510 ulv = p_lvid[i]
0511 nlv = n_lvid[i]
0512 lvn = self.get_lv_soname(ulv)
0513 csg = self.get_csg(ulv)
0514 tc = self.get_csg_typecode(ulv)
0515 assert len(csg) == len(tc)
0516 tcn = " ".join(list(map(lambda _:"%d:%s"%(_,CSG_.desc(_)), tc)))
0517 if pass_ == 1:
0518 lines.append("")
0519 pass
0520 lines.append(" lv:%3d nlv:%2d %50s csg %2d tcn %s " % (ulv, nlv, lvn, len(csg), tcn ))
0521 if pass_ == 1:
0522 lines.append(self.desc_csg(ulv))
0523 pass
0524 pass
0525 pass
0526 pass
0527 return "\n".join(lines)
0528
0529 def descSolids(self, detail=False):
0530 lines = []
0531 numSolid = self.get_numSolid()
0532 lines.append("stree.descSolids numSolid:%d detail:%d " % (numSolid,detail) )
0533
0534 q_ridx = int(os.environ.get("RIDX","-1"))
0535 for ridx in range(numSolid):
0536 if q_ridx > -1 and ridx != q_ridx: continue
0537 lines.append(self.descSolid(ridx, detail=detail))
0538 if detail:
0539 lines.append("")
0540 pass
0541 pass
0542 return "\n".join(lines)
0543
0544
0545
0546
0547 if __name__ == '__main__':
0548 logging.basicConfig(level=logging.INFO)
0549 f = Fold.Load(symbol="f")
0550
0551 st = stree(f)
0552 print(repr(st))
0553
0554 nidx = os.environ.get("NIDX", 13)
0555
0556 ancestors = st.get_ancestors(nidx)
0557 print("NIDX %d ancestors: %s " % (nidx, repr(ancestors)) )
0558 print("st.desc_nodes(ancestors)")
0559 print(st.desc_nodes(ancestors))
0560
0561 so = os.environ.get("SO", "sBar")
0562 lvs = st.find_lvid_(so)
0563 print("SO %s lvs %s" % (so, str(lvs)))
0564
0565 for lv in lvs:
0566 bb = st.find_lvid_nodes(lv)
0567 b = int(bb[0])
0568 print("lv:%d bb=st.find_lvid_nodes(lv) bb:%s b:%s " % (lv, str(bb),b))
0569
0570 anc = st.get_ancestors(b)
0571 print("b:%d anc=st.get_ancestors(b) anc:%s " % (b, str(anc)))
0572
0573 print("st.desc_nodes(anc, brief=True))")
0574 print(st.desc_nodes(anc, brief=True))
0575 print("st.desc_nodes([b], brief=True))")
0576 print(st.desc_nodes([b], brief=True))
0577 pass
0578
0579
0580
0581 if 0:
0582 progeny = st.get_progeny(nidx)
0583 print("NIDX %d progeny: %s " % (nidx, repr(progeny)) )
0584
0585
0586
0587 psf = st.make_freq(progeny)
0588 print("st.desc_nodes_(progeny, psf)")
0589 print(st.desc_nodes_(progeny, psf))
0590
0591 np.set_printoptions(edgeitems=600)
0592
0593 print("st.f.trs[progeny].reshape(-1,16)")
0594 print(st.f.trs[progeny].reshape(-1,16))
