File indexing completed on 2026-04-09 07:48:56
0001
0002 import os, sys, re, numpy as np, logging, datetime
0003 from configparser import ConfigParser
0004 from collections import OrderedDict, Counter
0005 log = logging.getLogger(__name__)
0006
0007
0008 from opticks.ana.fold import Fold, STR
0009 from opticks.sysrap.OpticksCSG import CSG_
0010
0011
0012
0013 class KeyNameConfig(object):
0014 """
0015 Parses ini file of the below form::
0016
0017 [key_boundary_regexp]
0018 red = (?P<Water_Virtuals>Water/.*/Water$)
0019 blue = Water///Acrylic$
0020 magenta = Water/Implicit_RINDEX_NoRINDEX_pOuterWaterInCD_pInnerReflector//Tyvek$
0021 yellow = DeadWater/Implicit_RINDEX_NoRINDEX_pDeadWater_pTyvekFilm//Tyvek$
0022 pink = Air/CDTyvekSurface//Tyvek$
0023 cyan = Tyvek//Implicit_RINDEX_NoRINDEX_pOuterWaterInCD_pCentralDetector/Water$
0024 orange = Water/Steel_surface//Steel$
0025 grey =
0026 # empty value is special cased to mean all other boundary names
0027
0028 Defaults path is $HOME/.opticks/GEOM/cxt_min.ini
0029 """
0030
0031
0032 def __init__(self, _path ):
0033 path = os.path.expandvars(_path)
0034 cfg = ConfigParser()
0035 cfg.read(path)
0036
0037 self._path = _path
0038 self.path = path
0039 self.cfg = cfg
0040
0041
0042 def __call__(self, _section):
0043 sect = self.cfg[_section]
0044 bdict = OrderedDict(sect)
0045
0046 counts = Counter(bdict.values())
0047 duplicates = [val for val, count in counts.items() if count > 1]
0048
0049 if duplicates:
0050 log.fatal(f"CSGFoundry.py/KeyNameConfig : Found duplicated values: {duplicates}")
0051 sys.exit(1)
0052 pass
0053
0054 return bdict
0055
0056
0057
0058
0059
0060 class NameTable(object):
0061 def __init__(self, symbol, ix, ixn, d_qix={}, KEY=None):
0062 """
0063 :param symbol: eg btab or ptab
0064 :param ix: large array of indices (eg boundaries or globalPrimIdx)
0065 :param ixn: smallish array of names
0066 :param d_qix: dict keyed on colors with arrays of ixn indices as values
0067 """
0068 u, x, c = np.unique(ix, return_index=True, return_counts=True)
0069
0070 if "KLUDGE" in os.environ:
0071 kludge = u != 0xffff
0072 u = u[kludge]
0073 x = x[kludge]
0074 c = c[kludge]
0075 pass
0076
0077
0078 akk = np.repeat("????????????????", len(u))
0079 for i in range(len(u)):
0080 kk = []
0081 for k, qix in d_qix.items():
0082 if len(np.where(np.isin(qix, u[i]))[0]) == 1:
0083 kk.append(k)
0084 pass
0085 pass
0086 akk[i] = ",".join(kk)
0087 pass
0088 n = ixn[u]
0089 o = np.argsort(c)[::-1]
0090
0091 _tab = "np.c_[akk, u, c, x, n][o]"
0092 tab = eval(_tab)
0093
0094 self.symbol = symbol
0095 self.akk = akk[o]
0096 self.u = u[o]
0097 self.x = x[o]
0098 self.c = c[o]
0099 self.n = n[o]
0100
0101 self.KEY = KEY
0102
0103 self.d_qix = d_qix
0104 self._tab = _tab
0105 self.tab = tab
0106 self.ixn = ixn
0107
0108 def get_names(self, k):
0109 """
0110 print("\n".join(cf.primname[cf.cxtp.d_qix['thistle']]))
0111 """
0112 return self.ixn[self.d_qix[k]]
0113
0114 def dump_names(self, k):
0115 names = self.get_names(k)
0116 print("\n".join(names))
0117
0118 def __repr__(self):
0119 """
0120 """
0121 vfmt = " %20s : %4d %8d %8d : %s "
0122 sfmt = " %20s : %4s %8s %8s : %s "
0123 tfmt = " %20s : %4s %8d %8s : %s "
0124
0125 label = sfmt % ( ".akk", ".u", ".c", ".x", ".n" )
0126 lines = []
0127 lines.append("[cf.%s KEY:%s" % (self.symbol, self.KEY) )
0128 lines.append(label)
0129
0130 if not self.KEY is None:
0131 if self.KEY.startswith("~"):
0132 KEY = self.KEY[1:]
0133 invert = True
0134 else:
0135 KEY = self.KEY
0136 invert = False
0137 pass
0138 KK = np.array(self.KEY.split(","))
0139 sel = np.where(np.isin(self.tab[:,0], KK, invert=invert))[0]
0140 else:
0141 sel = slice(None)
0142 pass
0143
0144 ctot = 0
0145 for row in self.tab[sel]:
0146 akk = row[0]
0147 u = int(row[1])
0148 c = int(row[2])
0149 x = int(row[3])
0150 n = row[4]
0151
0152 line = vfmt % ( akk, u, c, x, n )
0153 ctot += c
0154 lines.append(line)
0155 pass
0156 lines.append(label)
0157 lines.append(tfmt % ("", "",ctot, "",".ctot"))
0158 lines.append("]cf.%s" % self.symbol )
0159 self.ctot = ctot
0160 return "\n".join(lines)
0161
0162
0163
0164 class GroupedNameTable(object):
0165 def __init__(self, symbol, ix, d_qix, d_anno, cf_ixn, lwid=100):
0166 """
0167 :param symbol: identifier
0168 :param ix: large array of indices
0169 :param d_qix: dict keyed on colors with indices array values
0170 :param d_anno: dict keyed on colors with label values
0171 :param cf_ixn: small array of all names
0172 :param lwid: width of the label field
0173 """
0174 self.symbol = symbol
0175 self.ix = ix
0176 self.d_qix = d_qix
0177 self.d_anno = d_anno
0178 self.cf_ixn = cf_ixn
0179 self.lwid = lwid
0180
0181 wdict = {}
0182 for k,qix in self.d_qix.items():
0183 _w = np.isin( ix, qix )
0184 w = np.where(_w)[0]
0185 wdict[k] = w
0186 pass
0187 self.wdict = wdict
0188
0189 def __repr__(self):
0190 """
0191 """
0192 lines = []
0193 lines.append("[cf.%s___________________________________" % self.symbol )
0194 fmt = " %%15s %%%(lwid)ds nnn:%%4d %%s " % dict(lwid=self.lwid)
0195 for k,qix in self.d_qix.items():
0196 _w = np.isin( self.ix, qix )
0197 w = np.where(_w)[0]
0198 label = self.d_anno[k]
0199 nn = self.cf_ixn[qix]
0200 line = fmt % (k, label, len(nn), str(qix[:10]))
0201 lines.append(line)
0202 pass
0203 lines.append("]cf.%s_______ cf.%s.d_qix ____________________________" % ( self.symbol, self.symbol) )
0204 return "\n".join(lines)
0205
0206
0207
0208
0209
0210
0211 class CSGObject(object):
0212 @classmethod
0213 def Label(cls, spc=5, pfx=10):
0214 prefix = " " * pfx
0215 spacer = " " * spc
0216 return prefix + spacer.join(cls.FIELD)
0217
0218 @classmethod
0219 def Fields(cls, bi=False):
0220 kls = cls.__name__
0221 for i, field in enumerate(cls.FIELD):
0222 setattr(cls, field, i)
0223 if bi:setattr(builtins, field, i)
0224 pass
0225
0226 @classmethod
0227 def Type(cls):
0228 cls.Fields()
0229 kls = cls.__name__
0230 print("%s.Type()" % kls )
0231 for i, field in enumerate(cls.FIELD):
0232 name = cls.DTYPE[i][0]
0233 fieldname = "%s.%s" % (kls, field)
0234 print(" %2d : %20s : %s " % (i, fieldname, name))
0235 pass
0236 print("%s.Label() : " % cls.Label() )
0237
0238 @classmethod
0239 def RecordsFromArrays(cls, a):
0240 """
0241 :param a: ndarray
0242 :return: np.recarray
0243 """
0244 ra = np.core.records.fromarrays(a.T, dtype=cls.DTYPE )
0245 return ra
0246
0247
0248
0249 class CSGPrim(CSGObject):
0250 DTYPE = [
0251 ('numNode', '<i4'),
0252 ('nodeOffset', '<i4'),
0253 ('tranOffset', '<i4'),
0254 ('planOffset', '<i4'),
0255 ('sbtIndexOffset', '<i4'),
0256 ('meshIdx', '<i4'),
0257 ('repeatIdx', '<i4'),
0258 ('primIdx', '<i4'),
0259 ]
0260
0261 EXTRA = [
0262 ('BBMin_x', '<f4'),
0263 ('BBMin_y', '<f4'),
0264 ('BBMin_z', '<f4'),
0265 ('BBMax_x', '<f4'),
0266 ('BBMax_y', '<f4'),
0267 ('BBMax_z', '<f4'),
0268 ('spare32', '<f4'),
0269 ('spare33', '<f4'),
0270 ]
0271
0272 FIELD = "nn no to po sb lv ri pi".split()
0273 XFIELD = "nx ny nz mx my mz s2 s3".split()
0274
0275
0276
0277 class CSGNode(CSGObject):
0278 DTYPE = [
0279 ('p0', '<f4'),
0280 ('p1', '<f4'),
0281 ('p2', '<f4'),
0282 ('p3', '<f4'),
0283 ('p4', '<f4'),
0284 ('p5', '<f4'),
0285 ('boundary', '<i4'),
0286 ('index', '<i4'),
0287 ('BBMin_x', '<f4'),
0288 ('BBMin_y', '<f4'),
0289 ('BBMin_z', '<f4'),
0290 ('BBMax_x', '<f4'),
0291 ('BBMax_y', '<f4'),
0292 ('BBMax_z', '<f4'),
0293 ('typecode','<i4'),
0294 ('comptran','<i4'),
0295 ]
0296 FIELD = "p0 p1 p2 p3 p4 p5 bn ix nx ny nz mx my mz tc ct".split()
0297
0298
0299
0300
0301 class MM(object):
0302 """
0303 The input file is now a standard resident of the CSGFoundry directory
0304 Note that the input file was formerly created using::
0305
0306 ggeo.py --mmtrim > $TMP/mm.txt
0307
0308 """
0309 PTN = re.compile("\\d+")
0310 def __init__(self, path ):
0311 mm = os.path.expandvars(path)
0312 mm = open(mm, "r").read().splitlines() if os.path.exists(mm) else None
0313 self.mm = mm
0314 if mm is None:
0315 log.fatal("missing %s, which is now a standard part of CSGFoundry " % path )
0316 sys.exit(1)
0317 pass
0318
0319 def imm(self, emm):
0320 return list(map(int, self.PTN.findall(emm)))
0321
0322 def label(self, emm):
0323 imm = self.imm(emm)
0324 labs = [self.mm[i] for i in imm]
0325 lab = " ".join(labs)
0326
0327 tilde = emm[0] == "t" or emm[0] == "~"
0328 pfx = ( "NOT: " if tilde else " " )
0329
0330 if emm == "~0" or emm == "t0":
0331 return "ALL"
0332 elif imm == [1,2,3,4]:
0333 return "ONLY PMT"
0334 elif "," in emm:
0335 return ( "EXCL: " if tilde else "ONLY: " ) + lab
0336 else:
0337 return lab
0338 pass
0339
0340 def __repr__(self):
0341 return STR("\n".join(self.mm))
0342
0343
0344 class LV(object):
0345 PTN = re.compile("\\d+")
0346 def __init__(self, path):
0347 lv = os.path.expandvars(path)
0348 lv = open(lv, "r").read().splitlines() if os.path.exists(lv) else None
0349 self.lv = lv
0350 if lv is None:
0351 log.fatal("missing %s, which is now a standard part of CSGFoundry " % path )
0352 sys.exit(1)
0353 pass
0354
0355 def ilv(self, elv):
0356 return list(map(int, self.PTN.findall(elv)))
0357
0358 def label(self, elv):
0359 ilv = self.ilv(elv)
0360 mns = [self.lv[i] for i in ilv]
0361 mn = " ".join(mns)
0362 tilde = elv[0] == "t"
0363 lab = ""
0364 if elv == "t":
0365 lab = "ALL"
0366 else:
0367 lab = ( "EXCL: " if tilde else "ONLY: " ) + mn
0368 pass
0369 return lab
0370
0371 def __str__(self):
0372 return "\n".join(self.lv)
0373
0374 def __repr__(self):
0375 return "\n".join(["%3d:%s " % (i, self.lv[i]) for i in range(len(self.lv))])
0376
0377
0378
0379 class Deprecated_NPFold(object):
0380 """
0381 HMM opticks.ana.fold.Fold looks more developed than this
0382 TODO: eliminate all use of this
0383 """
0384 INDEX = "NPFold_index.txt"
0385
0386 @classmethod
0387 def IndexPath(cls, fold):
0388 return os.path.join(fold, cls.INDEX)
0389
0390 @classmethod
0391 def HasIndex(cls, fold):
0392 return os.path.exists(cls.IndexPath(fold))
0393
0394 @classmethod
0395 def Load(cls, *args, **kwa):
0396 return cls(*args, **kwa)
0397
0398 def __init__(self, *args, **kwa):
0399 fold = os.path.join(*args)
0400 self.fold = fold
0401 self.has_index = self.HasIndex(fold)
0402
0403 if self.has_index:
0404 self.load_idx(fold)
0405 else:
0406 self.load_fts(fold)
0407 pass
0408
0409 def load_fts(self, fold):
0410 assert 0
0411
0412 def load_idx(self, fold):
0413 keys = open(self.IndexPath(fold),"r").read().splitlines()
0414 aa = []
0415 kk = []
0416 ff = []
0417 subfold = []
0418
0419 for k in keys:
0420 path = os.path.join(fold, k)
0421 if k.endswith(".npy"):
0422 assert os.path.exists(path)
0423 log.info("loading path %s k %s " % (path, k))
0424 a = np.load(path)
0425 aa.append(a)
0426 kk.append(k)
0427 else:
0428 log.info("skip non .npy path %s k %s " % (path, k))
0429 pass
0430 pass
0431 self.kk = kk
0432 self.aa = aa
0433 self.ff = ff
0434 self.subfold = subfold
0435 self.keys = keys
0436
0437
0438 def find(self, k):
0439 return self.kk.index(k) if k in self.kk else -1
0440
0441 def has_key(self, k):
0442 return self.find(k) > -1
0443
0444 def __repr__(self):
0445 lines = []
0446 lines.append("NPFold %s " % self.fold )
0447 for i in range(len(self.kk)):
0448 k = self.kk[i]
0449 a = self.aa[i]
0450 stem, ext = os.path.splitext(os.path.basename(k))
0451 path = os.path.join(self.fold, k)
0452 lines.append("%10s : %20s : %s " % (stem, str(a.shape), k ))
0453 pass
0454 return "\n".join(lines)
0455
0456
0457 class Deprecated_SSim(object):
0458 """
0459 HMM: this adds little on top of ana.fold.Fold
0460 TODO: get rid of it
0461 """
0462 BND = "bnd.npy"
0463
0464
0465 @classmethod
0466 def Load(cls, simbase):
0467 ssdir = os.path.join(simbase, "SSim")
0468 log.info("SSim.Load simbase %s ssdir %s " % (simbase,ssdir))
0469 sim = Fold.Load(ssdir) if os.path.isdir(ssdir) else None
0470 return sim
0471
0472 def __init__(self, fold):
0473 if self.has_key(self.BND):
0474 bnpath = os.path.join(fold, "bnd_names.txt")
0475 assert os.path.exists(bnpath)
0476 bnd_names = open(bnpath,"r").read().splitlines()
0477 self.bnd_names = bnd_names
0478 pass
0479 if hasattr(self, 'bnd_names'):
0480 bndnamedict = SSim.namelist_to_namedict(self.bnd_names)
0481 else:
0482 bndnamedict = {}
0483 pass
0484 self.bndnamedict = bndnamedict
0485 pass
0486
0487
0488
0489 class CSGFoundry(object):
0490 FOLD = os.path.expandvars("$TMP/CSG_GGeo/CSGFoundry")
0491 FMT = " %10s : %20s : %s "
0492
0493
0494 @classmethod
0495 def namelist_to_namedict(cls, namelist):
0496 nd = {}
0497 if not namelist is None:
0498 nd = dict(zip(range(len(namelist)),list(map(str,namelist))))
0499 pass
0500 return nd
0501
0502
0503 @classmethod
0504 def CFBase(cls):
0505 """
0506 Precedence order for which envvars are used to derive the CFBase folder is:
0507
0508 1. CFBASE_LOCAL
0509 2. CFBASE_ALT
0510 3. CFBASE
0511 4. GEOM
0512
0513 This allows scripts such as CSGOptiX/cxt_min.sh to set envvars to control
0514 the geometry arrays and meshname.txt etc that is loaded.
0515 For example this could be used such that local analysis on laptop can use
0516 a geometry rsynced from remote.
0517
0518 """
0519 if "CFBASE_LOCAL" in os.environ:
0520 cfbase= os.environ["CFBASE_LOCAL"]
0521 note = "via CFBASE_LOCAL"
0522 elif "CFBASE_ALT" in os.environ:
0523 cfbase= os.environ["CFBASE_ALT"]
0524 note = "via CFBASE_ALT"
0525 elif "CFBASE" in os.environ:
0526 cfbase= os.environ["CFBASE"]
0527 note = "via CFBASE"
0528 elif "GEOM" in os.environ:
0529 KEY = os.path.expandvars("${GEOM}_CFBaseFromGEOM")
0530 if KEY in os.environ:
0531 cfbase = os.environ[KEY]
0532 note = "via GEOM,%s" % KEY
0533 else:
0534 cfbase = os.path.expandvars("$HOME/.opticks/GEOM/$GEOM")
0535 note = "via GEOM"
0536 pass
0537 else:
0538 cfbase = None
0539 note = "via NO envvars"
0540 pass
0541 if cfbase is None:
0542 print("CSGFoundry.CFBase returning None, note:%s " % note )
0543 else:
0544 print("CSGFoundry.CFBase returning [%s], note:[%s] " % (cfbase,note) )
0545 pass
0546 return cfbase
0547
0548 @classmethod
0549 def Load(cls, cfbase_=None, symbol=None):
0550 """
0551 :param cfbase_: typically None, but available when debugging eg comparing two geometries
0552 """
0553 if cfbase_ is None:
0554 cfbase = cls.CFBase()
0555 else:
0556 cfbase = os.path.expandvars(cfbase_)
0557 pass
0558 log.info("cfbase:%s " % cfbase)
0559 if cfbase is None or not os.path.exists(os.path.join(cfbase, "CSGFoundry")):
0560 print("ERROR CSGFoundry.CFBase returned None OR non-existing CSGFoundry dir so cannot CSGFoundry.Load" )
0561 return None
0562 pass
0563 assert not cfbase is None
0564 cf = cls(fold=os.path.join(cfbase, "CSGFoundry"))
0565 cf.symbol = symbol
0566 return cf
0567
0568 @classmethod
0569 def FindDirUpTree(cls, origpath, name="CSGFoundry"):
0570 """
0571 :param origpath: to traverse upwards looking for sibling dirs with *name*
0572 :param name: sibling directory name to look for
0573 :return full path to *name* directory
0574 """
0575 elem = origpath.split("/")
0576 found = None
0577 for i in range(len(elem),0,-1):
0578 path = "/".join(elem[:i])
0579 cand = os.path.join(path, name)
0580 log.debug(cand)
0581 if os.path.isdir(cand):
0582 found = cand
0583 break
0584 pass
0585 pass
0586 return found
0587
0588 @classmethod
0589 def FindDigest(cls, path):
0590 if path.find("*") > -1:
0591 base = os.path.dirname(path)
0592 else:
0593 base = path
0594 pass
0595 base = os.path.expandvars(base)
0596 return cls.FindDigest_(base)
0597
0598 @classmethod
0599 def FindDigest_(cls, path):
0600 """
0601 :param path: Directory path in which to look for a 32 character digest path element
0602 :return 32 character digest or None:
0603 """
0604 hexchar = "0123456789abcdef"
0605 digest = None
0606 for elem in path.split("/"):
0607 if len(elem) == 32 and set(elem).issubset(hexchar):
0608 digest = elem
0609 pass
0610 pass
0611 return digest
0612
0613 def __init__(self, fold):
0614 self.load(fold)
0615 self.meshnamedict = self.namelist_to_namedict(self.meshname)
0616 self.primIdx_meshname_dict = self.make_primIdx_meshname_dict()
0617
0618 self.mokname = "zero one two three four five six seven eight nine".split()
0619 self.moknamedict = self.namelist_to_namedict(self.mokname)
0620 self.insnamedict = {}
0621
0622 self.lv = LV(os.path.join(fold, "meshname.txt"))
0623 self.mm = MM(os.path.join(fold, "mmlabel.txt"))
0624
0625 sim = Fold.Load(fold, "SSim")
0626
0627 try:
0628 bdn = sim.stree.standard.bnd_names
0629 except AttributeError:
0630 bdn = None
0631 pass
0632 if bdn is None: log.fatal("CSGFoundry fail to access sim.stree.standard.bnd_names : geometry incomplete" )
0633 if type(bdn) is np.ndarray: sim.bndnamedict = self.namelist_to_namedict(bdn)
0634 pass
0635
0636 self.kncfg = KeyNameConfig("$HOME/.opticks/GEOM/cxt_min.ini")
0637
0638 self.bdn_config = self.kncfg("key_boundary_regexp")
0639 self.bdn_config_note = self.kncfg("key_boundary_regexp_note")
0640
0641 self.prn_config = self.kncfg("key_prim_regexp")
0642 self.prn_config_note = self.kncfg("key_prim_regexp_note")
0643
0644 self.bdn = bdn
0645 self.sim = sim
0646
0647 self.npa = self.node.reshape(-1,16)[:,0:6]
0648 self.nbd = self.node.view(np.int32)[:,1,2]
0649 self.nix = self.node.view(np.int32)[:,1,3]
0650 self.nbb = self.node.reshape(-1,16)[:,8:14]
0651 self.ntc = self.node.view(np.int32)[:,3,2]
0652 self.ncm = self.node.view(np.uint32)[:,3,3] >> 31
0653 self.ntr = self.node.view(np.uint32)[:,3,3] & 0x7fffffff
0654
0655 self.pnn = self.prim.view(np.int32)[:,0,0]
0656 self.pno = self.prim.view(np.int32)[:,0,1]
0657
0658
0659
0660 self.pto = self.prim.view(np.int32)[:,0,2]
0661 self.ppo = self.prim.view(np.int32)[:,0,3]
0662
0663 self.crn = self.pno[self.pnn>1]
0664 self.crn_subnum = self.node.view(np.int32)[self.crn,0,0]
0665 self.crn_suboff = self.node.view(np.int32)[self.crn,0,1]
0666
0667
0668 self.psb = self.prim.view(np.int32)[:,1,0]
0669 self.plv = self.prim.view(np.int32)[:,1,1]
0670 self.prx = self.prim.view(np.int32)[:,1,2]
0671 self.pix = self.prim.view(np.int32)[:,1,3]
0672
0673 self.pbb = self.prim.reshape(-1,16)[:,8:14]
0674
0675
0676 self.snp = self.solid[:,1,0].view(np.int32)
0677 self.spo = self.solid[:,1,1].view(np.int32)
0678 self.sce = self.solid[:,2].view(np.float32)
0679
0680
0681
0682 def simtrace_boundary_analysis(self, bn, KEY=os.environ.get("KEY", None) ):
0683 """
0684 Group simtrace intersects according to their boundary indices
0685 using groups specified by regexp matching boundary names.
0686
0687 Typically would have a few hundred boundary names
0688 in cf.bdn but potentially millions of simtrace intersects
0689 in the bn array.
0690
0691 :param bn: large array of boundary indices
0692 :param KEY: color key OR None
0693 :return cxtb,btab:
0694 """
0695 d_qbn, d_anno = self.Dict_find_name_indices_re_match(self.bdn, self.bdn_config, self.bdn_config_note )
0696 btab = NameTable("btab", bn, self.bdn, d_qbn, KEY)
0697 cxtb = GroupedNameTable("cxtb", bn, d_qbn, d_anno, self.bdn, lwid=100)
0698 self.btab = btab
0699 self.cxtb = cxtb
0700 return cxtb, btab
0701
0702 def simtrace_prim_analysis(self, pr, KEY=os.environ.get("KEY", None) ):
0703 """
0704 Group simtrace intersects according to their primitive indices
0705 using groups specified by regexp matching solid names of the prim.
0706
0707 :param pr: large array of primitive indices
0708 :param KEY:
0709 :return cxtp, ptab: GroupedNameTable, NameTable
0710 """
0711 d_qpr, d_anno = self.Dict_find_name_indices_re_match(self.primname, self.prn_config, self.prn_config_note )
0712 ptab = NameTable("ptab", pr, self.primname, d_qpr, KEY)
0713 cxtp = GroupedNameTable("cxtp", pr, d_qpr, d_anno, self.primname, lwid=50)
0714 self.ptab = ptab
0715 self.cxtp = cxtp
0716 return cxtp, ptab
0717
0718 @classmethod
0719 def Dict_find_name_indices_re_match(cls, names, _nameptn_dict, _namenote_dict ):
0720 """
0721 :param names: array of names
0722 :param _nameptn_dict: label keys, regexp pattern string values
0723 :return d: dict with same label keys and arrays of matching names indices
0724
0725 Names unmatched by the provided regexp values are
0726 grouped within the key with a blank value if provided.
0727 """
0728 d = {}
0729 anno = {}
0730 matched = np.array([], dtype=np.int64 )
0731 unmatched_k = None
0732 for k, v in _nameptn_dict.items():
0733 if v == '':
0734 unmatched_k = k
0735 continue
0736 pass
0737 qnm, nn, label = cls.Find_name_indices_re_match(names, v)
0738 matched = np.concatenate( (matched, qnm) )
0739 k_note = _namenote_dict.get(k,"")
0740 lab = "%100s ## %s" % ( label, k_note )
0741 d[k] = qnm
0742 anno[k] = lab
0743 pass
0744
0745 c = dict(matched=matched, all=np.arange(len(names)))
0746 c['unmatched'] = np.where(np.isin(c['all'], c['matched'], invert=True ))[0]
0747 if not unmatched_k is None:
0748 d[unmatched_k] = c['unmatched']
0749 anno[unmatched_k] = "OTHER"
0750 pass
0751 return d, anno
0752
0753
0754 @classmethod
0755 def Find_name_indices_re_match(cls, names, _ptn):
0756 """
0757 :param _ptn: name regexp string
0758 :return qnm,nn,label:
0759
0760 qnm
0761 array of names array indices that match the pattern
0762 nn
0763 array of names matching the pattern
0764 label
0765 match key when the pattern string includes one eg (?<label>.*$)
0766 otherwise the pattern string itself
0767
0768 Typically would have a few hundred to a few thousand names for a geometry
0769 so the slowness of np.vectorize is not a problem.
0770
0771 Example _ptn::
0772
0773 .*
0774 .*/Vacuum$
0775 .*/LatticedShellSteel$
0776 .*/Steel$
0777 Water/.*/Water$
0778 (?P<Water_Virtuals>:Water/.*/Water$)
0779 .*/Tyvek$
0780
0781 """
0782 ptn = re.compile(_ptn)
0783 _re_match = lambda s:bool(ptn.match(s))
0784 re_match = np.vectorize(_re_match)
0785
0786 def _re_key(s):
0787 keys = list(ptn.match(s).groupdict().keys())
0788 return keys[0] if len(keys) == 1 else _ptn
0789 pass
0790 re_key = np.vectorize(_re_key)
0791
0792 _qnm = re_match(names)
0793 qnm = np.where(_qnm)[0]
0794 nn = names[qnm]
0795
0796 if len(qnm) == 0:
0797 log.info("_ptn:[%s] did not match any names" % _ptn )
0798 label = "NO_MATCH"
0799 else:
0800 _kk = re_key(nn)
0801 kk = np.unique(_kk)
0802 assert len(kk) == 1
0803 label = kk[0]
0804 pass
0805 return qnm, nn, label
0806
0807
0808 def find_primIdx_from_nodeIdx(self, nodeIdx_):
0809 """
0810 If nodeIdx is valid there should always be exactly
0811 one prim in which it appears.
0812
0813 For example use to lookup primname that contain a selection
0814 of nodeIdx::
0815
0816 In [4]: np.c_[b.primname[np.unique(a.find_primIdx_from_nodeIdx(w))]]
0817 Out[4]:
0818 array([['NNVTMCPPMTsMask_virtual'],
0819 ['HamamatsuR12860sMask_virtual'],
0820 ['HamamatsuR12860_PMT_20inch_pmt_solid_1_4'],
0821 ['HamamatsuR12860_PMT_20inch_inner_solid_1_4'],
0822 ['base_steel'],
0823 ['uni_acrylic1']], dtype=object)
0824
0825
0826 """
0827 a = self
0828
0829 numNode = a.prim[:,0,0].view(np.int32)
0830 nodeOffset = a.prim[:,0,1].view(np.int32)
0831
0832 primIdx = np.zeros(len(nodeIdx_), dtype=np.int32)
0833 for i, nodeIdx in enumerate(nodeIdx_):
0834 primIdx[i] = np.where(np.logical_and( nodeIdx >= nodeOffset, nodeIdx < nodeOffset+numNode ))[0]
0835 pass
0836 return primIdx
0837
0838
0839 def find_lvid_from_nodeIdx(self, nodeIdx_):
0840 """
0841 """
0842 primIdx = self.find_primIdx_from_nodeIdx(nodeIdx_)
0843 return self.prim[primIdx].view(np.int32)[:,1,1]
0844
0845 def find_lvname_from_nodeIdx(self, nodeIdx_):
0846 lvid = self.find_lvid_from_nodeIdx(nodeIdx_)
0847 return self.meshname[lvid]
0848
0849
0850
0851 def meshIdx(self, primIdx):
0852 """
0853 Lookup the midx of primIdx prim
0854
0855 :param primIdx:
0856 :return midx:
0857 """
0858 assert primIdx < len(self.prim)
0859 midx = self.prim[primIdx].view(np.uint32)[1,1]
0860 return midx
0861
0862 def make_primIdx_meshname_dict(self):
0863 """
0864 See notes/issues/cxs_2d_plotting_labels_suggest_meshname_order_inconsistency.rst
0865 this method resolved an early naming bug
0866
0867 CSG/CSGPrim.h::
0868
0869 95 PRIM_METHOD unsigned meshIdx() const { return q1.u.y ; } // aka lvIdx
0870 96 PRIM_METHOD void setMeshIdx(unsigned midx){ q1.u.y = midx ; }
0871
0872 """
0873 d = {}
0874 for primIdx in range(len(self.prim)):
0875 midx = self.meshIdx (primIdx)
0876 assert midx < len(self.meshname)
0877 mnam = self.meshname[midx]
0878 d[primIdx] = mnam
0879
0880 pass
0881 return d
0882
0883 def loadtxt(self, path):
0884 """
0885 When the path contains only a single line using dtype np.object or |S100
0886 both result in an object with no length::
0887
0888 array('solidXJfixture', dtype=object)
0889
0890 """
0891
0892
0893
0894
0895
0896 txt = open(path).read().splitlines()
0897 a_txt = np.array(txt, dtype=np.str_ )
0898 return a_txt
0899
0900
0901 def brief(self):
0902 symbol = getattr(self, "symbol", "no-symbol")
0903 return "CSGFoundry %s cfbase %s " % (symbol, self.cfbase)
0904
0905 def load(self, fold):
0906 cfbase = os.path.dirname(fold)
0907 self.cfbase = cfbase
0908 self.base = cfbase
0909 log.info("load : fold %s cfbase: %s " % (fold, cfbase) )
0910
0911 if not os.path.isdir(fold):
0912 log.fatal("CSGFoundry folder %s does not exist " % fold)
0913 log.fatal("create foundry folder from OPTICKS_KEY geocache with CSG_GGeo/run.sh " )
0914 assert 0
0915 pass
0916
0917 names = os.listdir(fold)
0918 stems = []
0919 stamps = []
0920 for name in filter(lambda name:name.endswith(".npy") or name.endswith(".txt"), names):
0921 path = os.path.join(fold, name)
0922 st = os.stat(path)
0923 stamp = datetime.datetime.fromtimestamp(st.st_ctime)
0924 stamps.append(stamp)
0925 stem = name[:-4]
0926 a = np.load(path) if name.endswith(".npy") else self.loadtxt(path)
0927 setattr(self, stem, a)
0928 stems.append(stem)
0929 pass
0930
0931 min_stamp = min(stamps)
0932 max_stamp = max(stamps)
0933 now_stamp = datetime.datetime.now()
0934 dif_stamp = max_stamp - min_stamp
0935 age_stamp = now_stamp - max_stamp
0936
0937
0938
0939
0940
0941
0942 assert dif_stamp.microseconds < 1e6, "stamp divergence detected microseconds %d : so are seeing mixed up results from multiple runs " % dif_stamp.microseconds
0943
0944 self.stems = stems
0945 self.min_stamp = min_stamp
0946 self.max_stamp = max_stamp
0947 self.age_stamp = age_stamp
0948 self.stamps = stamps
0949 self.fold = fold
0950 self.pr = CSGPrim.RecordsFromArrays(self.prim[:,:2].reshape(-1,8).view(np.int32))
0951 self.nf = CSGNode.RecordsFromArrays(self.node.reshape(-1,16).view(np.float32))
0952 self.ni = CSGNode.RecordsFromArrays(self.node.reshape(-1,16).view(np.int32))
0953
0954
0955 def desc(self, stem):
0956 a = getattr(self, stem)
0957 ext = ".txt" if a.dtype == 'O' else ".npy"
0958 pstem = "bnd" if stem == "bndname" else stem
0959 path = os.path.join(self.fold, "%s%s" % (pstem, ext))
0960 return self.FMT % (stem, str(a.shape), path)
0961
0962 def head(self):
0963 return "\n".join(map(str, [self.fold, "min_stamp:%s" % self.min_stamp, "max_stamp:%s" % self.max_stamp, "age_stamp:%s" % self.age_stamp]))
0964
0965 def body(self):
0966 return "\n".join(map(lambda stem:self.desc(stem),self.stems))
0967
0968 def __repr__(self):
0969 return "\n".join([self.head(),self.body()])
0970
0971 def dump_node_boundary(self):
0972 logging.info("dump_node_boundary")
0973 node = self.node
0974
0975 node_boundary = node.view(np.uint32)[:,1,2]
0976 ubs, ubs_count = np.unique(node_boundary, return_counts=True)
0977
0978 for ub, ub_count in zip(ubs, ubs_count):
0979 bn = self.bdn[ub] if not self.bdn is None else "-"
0980 print(" %4d : %6d : %s " % (ub, ub_count, bn))
0981 pass
0982
0983 def descSolid(self, ridx, detail=True):
0984 """
0985 After CSGFoundry::dumpSolid
0986 """
0987 label = self.solid[ridx,0,:4].copy().view("|S16")[0].decode("utf8")
0988 numPrim = self.solid[ridx,1,0]
0989 primOffset = self.solid[ridx,1,1]
0990 prs = self.prim[primOffset:primOffset+numPrim]
0991 iprs = prs.view(np.int32)
0992
0993 lvid = iprs[:,1,1]
0994 u_lvid, x_lvid, n_lvid = np.unique( lvid, return_index=True, return_counts = True )
0995 lv_one = np.all( n_lvid == 1 )
0996
0997 lines = []
0998 lines.append("CSGFoundry.descSolid ridx %2d label %16s numPrim %6d primOffset %6d lv_one %d " % (ridx,label, numPrim, primOffset, lv_one))
0999
1000 if detail:
1001 if lv_one:
1002 for pidx in range(primOffset, primOffset+numPrim):
1003 lines.append(self.descPrim(pidx))
1004 pass
1005 else:
1006 x_order = np.argsort(x_lvid)
1007 for i in x_order:
1008 ulv = u_lvid[i]
1009 nlv = n_lvid[i]
1010 xlv = x_lvid[i]
1011 pidx = primOffset + xlv
1012 lines.append(" i %3d ulv %3d xlv %4d nlv %3d : %s " % (i, ulv, xlv, nlv, self.descPrim(pidx) ))
1013 pass
1014 pass
1015 pass
1016 return "\n".join(lines)
1017
1018 def descLV(self, lvid, detail=False):
1019 lv = self.prim.view(np.int32)[:,1,1]
1020 pidxs = np.where(lv==lvid)[0]
1021
1022 lines = []
1023 lines.append("descLV lvid:%d meshname:%s pidxs:%s" % (lvid, self.meshname[lvid],str(pidxs)))
1024 for pidx in pidxs:
1025 lines.append(self.descPrim(pidx, detail=detail))
1026 pass
1027 return STR("\n".join(lines))
1028
1029 def descLVDetail(self, lvid):
1030 return self.descLV(lvid, detail=True)
1031
1032 def descPrim(self, pidx, detail=False):
1033
1034 pr = self.prim[pidx]
1035 ipr = pr.view(np.int32)
1036
1037 numNode = ipr[0,0]
1038 nodeOffset = ipr[0,1]
1039 repeatIdx = ipr[1,2]
1040 primIdxLocal = ipr[1,3]
1041
1042 lvid = ipr[1,1]
1043 lvn = self.meshname[lvid]
1044 tcn = self.descNodeTC(nodeOffset, numNode)
1045
1046 bnd_ = self.node[nodeOffset:nodeOffset+numNode,1,2].view(np.int32)
1047 ubnd_ = np.unique(bnd_)
1048 assert len(ubnd_) == 1, "all nodes of prim should have same boundary "
1049 ubnd = ubnd_[0]
1050 line = "pidx %4d lv %3d pxl %4d : %50s : %s : bnd %s : %s " % (pidx, lvid, primIdxLocal, lvn, tcn, ubnd, self.bdn[ubnd] )
1051
1052 lines = []
1053 lines.append(line)
1054 if detail:
1055 lines.append(self.descNodeParam(nodeOffset,numNode))
1056 lines.append(self.descNodeBB(nodeOffset,numNode))
1057 lines.append(self.descNodeBoundaryIndex(nodeOffset,numNode))
1058 lines.append(self.descNodeTCTran(nodeOffset,numNode, mask_signbit=False))
1059 lines.append(self.descNodeTCTran(nodeOffset,numNode, mask_signbit=True))
1060 pass
1061 return STR("\n".join(lines))
1062
1063 def descPrimDetail(self, pidx):
1064 return self.descPrim(pidx,detail=True)
1065
1066 def descNodeFloat(self, nodeOffset, numNode, sli="[:,:6]", label=""):
1067 symbol = self.symbol
1068 locals()[symbol] = self
1069 expr = "%(symbol)s.node[%(nodeOffset)d:%(nodeOffset)d+%(numNode)d].reshape(-1,16)"
1070 expr += sli
1071 expr = expr % locals()
1072 lines = []
1073 lines.append("%s # %s " % (expr,label))
1074 lines.append(str(eval(expr)))
1075 return STR("\n".join(lines))
1076
1077 def descNodeParam(self, nodeOffset, numNode):
1078 return self.descNodeFloat(nodeOffset,numNode,"[:,:6]", "descNodeParam" )
1079 def descNodeBB(self, nodeOffset, numNode):
1080 return self.descNodeFloat(nodeOffset,numNode,"[:,8:14]", "descNodeBB" )
1081
1082 def descNodeInt(self, nodeOffset, numNode, sli="[:,6:8]", label="", mask_signbit=False):
1083 symbol = self.symbol
1084 locals()[symbol] = self
1085 expr = "%(symbol)s.node[%(nodeOffset)d:%(nodeOffset)d+%(numNode)d].reshape(-1,16).view(np.int32)"
1086 expr += sli
1087 if mask_signbit:
1088 expr += " & 0x7ffffff "
1089 pass
1090 expr = expr % locals()
1091 lines = []
1092 lines.append("%s # %s " % (expr,label))
1093 lines.append(str(eval(expr)))
1094 return STR("\n".join(lines))
1095
1096 def descNodeBoundaryIndex(self, nodeOffset, numNode):
1097 return self.descNodeInt(nodeOffset,numNode,"[:,6:8]", "descNodeBoundaryIndex" )
1098 def descNodeTCTran(self, nodeOffset, numNode, mask_signbit=False):
1099 return self.descNodeInt(nodeOffset,numNode,"[:,14:16]", "descNodeTCTran", mask_signbit=mask_signbit )
1100
1101
1102
1103
1104 def descNodeTC(self, nodeOffset, numNode, sumcut=7):
1105 """
1106 :param nodeOffset: absolute index within self.node array
1107 :param numNode: number of contiguous nodes
1108 :param sumcut: number of nodes above which a summary output is returned
1109 :return str: representing CSG typecodes
1110
1111 (nodeOffset,numNode) will normally correspond to the node range of a single prim pidx
1112
1113 Output examples::
1114
1115 tcn 1:union 1:union 108:cone 105:cylinder 105:cylinder 0:zero 0:zero
1116 tcn 2:intersection 1:union 2:intersection 103:zsphere 105:cylinder 103:!zsphere 105:!cylinder
1117 tcn 3(2:intersection) 2(1:union) 4(105:cylinder) 2(103:zsphere) 4(0:zero)
1118
1119 """
1120 nds = self.node[nodeOffset:nodeOffset+numNode].view(np.uint32)
1121 neg = ( nds[:,3,3] & 0x80000000 ) >> 31
1122 tcs = nds[:,3,2]
1123
1124 if len(tcs) > sumcut:
1125 u_tc, x_tc, n_tc = np.unique(tcs, return_index=True, return_counts=True)
1126 x_order = np.argsort(x_tc)
1127 tce = []
1128 for i in x_order:
1129 utc = u_tc[i]
1130 ntc = n_tc[i]
1131 xtc = x_tc[i]
1132 xng = neg[xtc]
1133 tce.append( "%d(%d:%s%s)" % (ntc,utc,"!" if xng else "", CSG_.desc(utc) ))
1134 pass
1135 tcn = " ".join(tce)
1136 else:
1137 tcn = " ".join(list(map(lambda _:"%d:%s%s" % (tcs[_],"!" if neg[_] else "",CSG_.desc(tcs[_])),range(len(tcs)))))
1138 pass
1139
1140 return "no %5d nn %4d tcn %s tcs %s" % (nodeOffset, numNode, tcn, str(tcs) )
1141
1142
1143
1144 def descSolids(self, detail=True):
1145 num_solid = len(self.solid)
1146 lines = []
1147 q_ridx = int(os.environ.get("RIDX", "-1"))
1148 for ridx in range(num_solid):
1149 if q_ridx > -1 and ridx != q_ridx: continue
1150 if detail:
1151 lines.append("")
1152 pass
1153 lines.append(self.descSolid(ridx, detail=detail))
1154 pass
1155 return STR("\n".join(lines))
1156
1157
1158 if __name__ == '__main__':
1159 logging.basicConfig(level=logging.INFO)
1160
1161 cf = CSGFoundry.Load()
1162 print(cf)
1163
1164
1165
1166
1167
1168
1169
1170
