From ee838c2ba1631ba1db4ba2dfba525f28d9c09823 Mon Sep 17 00:00:00 2001 From: Dmitri Smirnov Date: Thu, 7 May 2026 18:42:37 -0400 Subject: [PATCH 1/5] chore: remove sysrap/*.py --- sysrap/ABR.py | 109 --- sysrap/OpticksCSG.py | 80 --- sysrap/SABTest.py | 191 ----- sysrap/SCenterExtentGenstep.py | 51 -- sysrap/SProp.py | 100 --- sysrap/__init__.py | 19 - sysrap/dv.py | 106 --- sysrap/sevt.py | 1211 -------------------------------- sysrap/sevt_load.py | 17 - sysrap/sevt_tt.py | 576 --------------- sysrap/sframe.py | 463 ------------ sysrap/sfreq.py | 75 -- sysrap/smath.py | 111 --- sysrap/smonitor.py | 108 --- sysrap/sn_check.py | 121 ---- sysrap/sphotonlite.py | 258 ------- sysrap/stag.py | 275 -------- sysrap/stree.py | 594 ---------------- sysrap/xfold.py | 174 ----- 19 files changed, 4639 deletions(-) delete mode 100644 sysrap/ABR.py delete mode 100644 sysrap/OpticksCSG.py delete mode 100644 sysrap/SABTest.py delete mode 100644 sysrap/SCenterExtentGenstep.py delete mode 100755 sysrap/SProp.py delete mode 100644 sysrap/__init__.py delete mode 100644 sysrap/dv.py delete mode 100644 sysrap/sevt.py delete mode 100644 sysrap/sevt_load.py delete mode 100644 sysrap/sevt_tt.py delete mode 100755 sysrap/sframe.py delete mode 100644 sysrap/sfreq.py delete mode 100644 sysrap/smath.py delete mode 100644 sysrap/smonitor.py delete mode 100644 sysrap/sn_check.py delete mode 100644 sysrap/sphotonlite.py delete mode 100755 sysrap/stag.py delete mode 100644 sysrap/stree.py delete mode 100644 sysrap/xfold.py diff --git a/sysrap/ABR.py b/sysrap/ABR.py deleted file mode 100644 index 5a066da100..0000000000 --- a/sysrap/ABR.py +++ /dev/null @@ -1,109 +0,0 @@ -#!/usr/bin/env python - -class ABR(object): - """ - As this takes a general side by side repr approach it belongs in its own module. - """ - def __init__(self, A, B, symbol="AB"): - self.A = A - self.B = B - self.symbol = symbol - - def _get_idx(self): - return self.A.idx - def _set_idx(self, idx): - self.A.idx = idx - self.B.idx = idx - idx = property(_get_idx, _set_idx) - - def _get_flavor(self): - return self.A.flavor - def _set_flavor(self, flavor): - self.A.flavor = flavor - self.B.flavor = flavor - flavor = property(_get_flavor, _set_flavor) - - - def __call__(self, idx): - self.idx = idx - self.flavor = "call" - return self - def __getitem__(self, idx): - self.idx = idx - self.flavor = "getitem" - return self - - - @classmethod - def MaxLen(cls, lines): - mx = 0 - for i in range(len(lines)): - if len(lines[i]) > mx: mx = len(lines[i]) - pass - return mx - - @classmethod - def LeftRightPaddedSplit(cls, lhs, rhs): - """ - :param lhs: text delimited by newlines - :param rhs: text delimited by newlines - :return l,r: lhs and rhs lines padded to be the same length - """ - l = lhs.split("\n") - r = rhs.split("\n") - nl = len(l) - nr = len(r) - if nl == nr: - pass - elif nl > nr: - r.extend( ["" for _ in range(nl-nr)] ) - nr = len(r) - elif nr > nl: - l.extend( ["" for _ in range(nr-nl)] ) - nl = len(l) - else: - pass - pass - assert nl == nr - return l, r - - @classmethod - def SideBySide(cls, lhs, rhs, margin=10): - l, r = cls.LeftRightPaddedSplit(lhs, rhs) - nl = len(l) - nr = len(r) - assert nl == nr - - lx_ = cls.MaxLen(l) - rx_ = cls.MaxLen(r) - lim = 100 - - lx = min(lim,lx_) - rx = min(lim,rx_) - - pass - fmt_ = "{:%d}{:%d}" - fmt = fmt_ % (lx+margin, rx+margin) - - lines = [] - for i in range(nl): - line = fmt.format(l[i], r[i]) - lines.append(line.rstrip()) - pass - return "\n".join(lines) - - def identification(self): - ia = self.A.identification() if hasattr(self.A, "identification") else "" - ib = self.B.identification() if hasattr(self.B, "identification") else "" - return "\n".join(filter(None,[ia,ib])) - - def __repr__(self): - return "\n".join([self.identification(), self.SideBySide(repr(self.A), repr(self.B))]) - - -if __name__ == '__main__': - pass - - - - diff --git a/sysrap/OpticksCSG.py b/sysrap/OpticksCSG.py deleted file mode 100644 index 66cb1d4dac..0000000000 --- a/sysrap/OpticksCSG.py +++ /dev/null @@ -1,80 +0,0 @@ -# generated Tue Feb 7 17:01:08 2023 -# from /Users/blyth/opticks/sysrap -# base OpticksCSG.h stem OpticksCSG -# with command : ../bin/c_enums_to_python.py OpticksCSG.h -#0 -class CSG_(object): - ZERO = 0 - OFFSET_LIST = 4 - OFFSET_LEAF = 7 - TREE = 1 - UNION = 1 - INTERSECTION = 2 - DIFFERENCE = 3 - NODE = 11 - LIST = 11 - CONTIGUOUS = 11 - DISCONTIGUOUS = 12 - OVERLAP = 13 - LEAF = 101 - SPHERE = 101 - BOX = 102 - ZSPHERE = 103 - TUBS = 104 - CYLINDER = 105 - SLAB = 106 - PLANE = 107 - CONE = 108 - EXBB = 109 - BOX3 = 110 - TRAPEZOID = 111 - CONVEXPOLYHEDRON = 112 - DISC = 113 - SEGMENT = 114 - ELLIPSOID = 115 - TORUS = 116 - HYPERBOLOID = 117 - CUBIC = 118 - INFCYLINDER = 119 - OLDCYLINDER = 120 - PHICUT = 121 - THETACUT = 122 - OLDCONE = 123 - UNDEFINED = 124 - OBSOLETE = 1000 - PARTLIST = 1001 - FLAGPARTLIST = 1002 - FLAGNODETREE = 1003 - FLAGINVISIBLE = 1004 - PMT = 1005 - ZLENS = 1006 - PRISM = 1007 - LAST = 1008 - D2V={'zero': 0, 'intersection': 2, 'union': 1, 'difference': 3, 'contiguous': 11, 'discontiguous': 12, 'overlap': 13, 'sphere': 101, 'box': 102, 'zsphere': 103, 'tubs': 104, 'cylinder': 105, 'slab': 106, 'plane': 107, 'cone': 108, 'oldcone': 123, 'box3': 110, 'trapezoid': 111, 'convexpolyhedron': 112, 'disc': 113, 'segment': 114, 'ellipsoid': 115, 'torus': 116, 'hyperboloid': 117, 'cubic': 118, 'infcylinder': 119, 'oldcylinder': 120, 'phicut': 121, 'thetacut': 122, 'undefined': 124, 'externalbb': 109, 'obsolete': 1000, 'partlist': 1001, 'flagpartlist': 1002, 'flagnodetree': 1003, 'flaginvisible': 1004, 'pmt': 1005, 'zlens': 1006, 'prism': 1007, 'last': 1008} - - - @classmethod - def raw_enum(cls): - return list(filter(lambda kv:type(kv[1]) is int,cls.__dict__.items())) - - @classmethod - def enum(cls): - return cls.D2V.items() if len(cls.D2V) > 0 else cls.raw_enum() - - @classmethod - def desc(cls, typ): - kvs = list(filter(lambda kv:kv[1] == typ, cls.enum())) - return kvs[0][0] if len(kvs) == 1 else "UNKNOWN" - - @classmethod - def descmask(cls, typ): - kvs = list(filter(lambda kv:kv[1] & typ, cls.enum())) - return ",".join(map(lambda kv:kv[0], kvs)) - - @classmethod - def fromdesc(cls, label): - kvs = list(filter(lambda kv:kv[0] == label, cls.enum())) - return kvs[0][1] if len(kvs) == 1 else -1 - - - diff --git a/sysrap/SABTest.py b/sysrap/SABTest.py deleted file mode 100644 index 6588d4b2d5..0000000000 --- a/sysrap/SABTest.py +++ /dev/null @@ -1,191 +0,0 @@ -#!/usr/bin/env python -""" -SABTest.py : analysis script used from SABTest.sh -================================================== - -Started from ~/j/InputPhotonsCheck/InputPhotonsCheck.py - -""" - -import os, logging, numpy as np -log = logging.getLogger(__name__) -os.environ["MODE"] = "3" - -from opticks.sysrap.sevt import SEvt, SAB, SABHit -from opticks.ana.p import cf - -MODE = int(os.environ.get("MODE","3")) -MODE0 = MODE -print("SABTest.py initial MODE:%d " % MODE) - - -GLOBAL = int(os.environ.get("GLOBAL","0")) == 1 -PICK = os.environ.get("PICK","A") -INCPOI = float(os.environ.get("INCPOI","0")) ## increment pyvista presentation point and line size, often need -5 - - -if MODE in [2,3]: - import opticks.ana.pvplt as pvp - print(pvp.__doc__) - print("after import MODE:%d " % MODE) - MODE = MODE0 - # HMM this import overrides MODE, so need to keep defaults the same -pass - - -if __name__ == '__main__': - logging.basicConfig(level=logging.INFO) - - INST_FRAME = int(os.environ.get("INST_FRAME","-1")) - if INST_FRAME == -1: - print("INST_FRAME %d : using default sframe saved with the sevt" % INST_FRAME ) - M2W_OVERRIDE = None - W2M_OVERRIDE = None - else: - print("INST_FRAME %d : W2M_OVERRIDE obtained from cf.inst controlled by envvar INST_FRAME " % INST_FRAME ) - M2W_OVERRIDE = np.eye(4) - M2W_OVERRIDE[:,:3] = cf.inst[INST_FRAME][:,:3] - W2M_OVERRIDE = np.linalg.inv(M2W_OVERRIDE) - pass - - - print(__doc__) - - a = SEvt.Load("$AFOLD", symbol="a", W2M=W2M_OVERRIDE) - b = SEvt.Load("$BFOLD", symbol="b", W2M=W2M_OVERRIDE) - print(repr(a)) - print(repr(b)) - - - if not a is None and b is None: - PICK = "A" - elif not b is None and a is None: - PICK = "B" - else: - pass - pass - - if not a is None and not b is None: - if "SAB" in os.environ: - print("[--- ab = SAB(a,b) ----") - ab = SAB(a,b) - print("]--- ab = SAB(a,b) ----") - - print("[----- repr(ab) ") - print(repr(ab)) - print("]----- repr(ab) ") - else: - print(" NOT doing ab = SAB(a,b) # do that with : export SAB=1" ) - pass - - print("[--- abh = SABHit(a,b) ----") - abh = SABHit(a,b) - print("]--- abh = SABHit(a,b) ----") - print("[----- repr(abh) ") - print(repr(abh)) - print("]----- repr(abh) ") - pass - - assert PICK in ["A","B","AB","BA", "CF"] - if PICK == "A": - ee = [a,] - elif PICK == "B": - ee = [b,] - elif PICK == "AB": - ee = [a,b,] - elif PICK == "BA": - ee = [b,a,] - elif PICK == "CF": - ee = [] - pass - - context = "PICK=%s MODE=%d ~/j/jtds/jtds.sh " % (PICK, MODE ) - print(context) - - pl = None - for e in ee: - if e is None:continue - elabel = "%s : %s " % ( e.symbol.upper(), e.f.base ) - - if hasattr(e.f, 'photon'): - pos = e.f.photon[:,0,:3] - elif hasattr(e.f, 'hit'): - pos = e.f.hit[:,0,:3] - else: - log.info("%s:sevt lacks photon or hit" % e.symbol) - pos = None - pass - - if hasattr(e.f, 'record'): - sel = np.where(e.f.record[:,:,2,3] > 0) # select on wavelength to avoid unfilled zeros - poi = e.f.record[:,:,0,:3][sel] - else: - log.info("%s:sevt lacks record" % e.symbol) - poi = None - pass - - if hasattr(e.f, 'sframe'): - _W2M = e.f.sframe.w2m - log.info("%s:sevt has sframe giving _W2M " % e.symbol) - else: - _W2M = np.eye(4) - log.info("%s:sevt lacks sframe" % e.symbol) - pass - - W2M = _W2M if W2M_OVERRIDE is None else W2M_OVERRIDE - - - print("W2M_OVERRIDE (from manual INST_FRAME, no longer typical)\n",W2M_OVERRIDE) - print("_W2M(from %s.f.sframe.w2m, now standard)\n" % e.symbol,_W2M) - print("W2M\n",W2M) - - if not pos is None: - print("transform gpos (from photon or hit array) to lpos using W2M : ie into target frame : GLOBAL %d " % GLOBAL) - gpos = np.ones( [len(pos), 4 ] ) - gpos[:,:3] = pos - lpos = np.dot( gpos, W2M ) # hmm unfilled global zeros getting transformed somewhere - upos = gpos if GLOBAL else lpos - else: - upos = None - pass - - if not poi is None: - print("transform gpoi (from record array) to lpoi using W2M : ie into target frame : GLOBAL %d " % GLOBAL) - gpoi = np.ones( [len(poi), 4 ] ) - gpoi[:,:3] = poi - lpoi = np.dot( gpoi, W2M ) - upoi = gpoi if GLOBAL else lpoi - else: - upoi = None - pass - - - label = context + " ## " + elabel - - if MODE == 3 and not "SAB" in os.environ: - if pvp.pv is None: - print("MODE is 3 but pvp.pv is None, get into ana environment first eg with hookup_conda_ok") - else: - pl = pvp.pvplt_plotter(label) - pvp.pvplt_viewpoint(pl) # sensitive EYE, LOOK, UP, ZOOM envvars eg EYE=0,-3,0 - - if not upoi is None: - print("add_points green for the upoi : record points ") - pl.add_points( upoi[:,:3], color="green", point_size=3.0 ) - pass - if not upos is None: - print("add_points red for the upos : photon/hit points ") - pl.add_points( upos[:,:3], color="red", point_size=3.0 ) - pass - pl.show_grid() - pl.increment_point_size_and_line_width(INCPOI) - cp = pl.show() - pass - pass - pass - - - -pass - - diff --git a/sysrap/SCenterExtentGenstep.py b/sysrap/SCenterExtentGenstep.py deleted file mode 100644 index cab0f2e5c4..0000000000 --- a/sysrap/SCenterExtentGenstep.py +++ /dev/null @@ -1,51 +0,0 @@ -#!/usr/bin/env python -import numpy as np -from opticks.ana.fold import Fold - -X,Y,Z,W=0,1,2,3 - -class SCenterExtentGenstep(object): - - @classmethod - def Lim(cls, xyz): - lim = np.zeros( (3,2), dtype=np.float32 ) - for A in [X,Y,Z]: lim[A] = [xyz[:,A].min(), xyz[:,A].max() ] - return lim - - - def __init__(self, path): - if path is None: - path = "/tmp/$USER/opticks/SCenterExtentGenstep" - pass - - fold = Fold.Load(path) - - photons = fold.photons - gs = fold.gs - - self.fold = fold - self.peta = fold.peta - self.gs = gs - self.photons = photons - self.plim = self.Lim( photons[:,0] ) - self.glim = self.Lim( gs[:,5] ) - - -if __name__ == '__main__': - - - #path = os.path.expandvars("$CFBASE/CSGIntersectSolidTest/$GEOM") - path = None - - cegs = SCenterExtentGenstep(path) - - peta = cegs.peta - gs = cegs.gs - photons = cegs.photons - - print("cegs.plim\n", cegs.plim) - print("cegs.glim\n", cegs.glim) - - - - diff --git a/sysrap/SProp.py b/sysrap/SProp.py deleted file mode 100755 index dfe018ba0b..0000000000 --- a/sysrap/SProp.py +++ /dev/null @@ -1,100 +0,0 @@ -#!/usr/bin/env python -""" -SProp.py -========== - - -""" -import os, logging, numpy as np -from collections import OrderedDict as odict -log = logging.getLogger(__name__) - -class SProp(object): - BASE = os.environ.get("NP_PROP_BASE", "/tmp") - - @classmethod - def Resolve_(cls, spec): - relp = spec.replace(".","/") - return os.path.join(cls.BASE, relp) - - @classmethod - def Resolve(cls, spec): - return cls.Resolve_(spec) if spec.count(".") > 1 else spec - - @classmethod - def Names(cls, spec): - rawnames = os.listdir(cls.Resolve(spec)) - names = [] - for name in rawnames: - if name.count(".") == 0: - names.append(name) - pass - pass - return names - - @classmethod - def Load(cls, spec): - path = cls.Resolve(spec) - log.debug(" spec %s path %s " % (spec, path)) - try: - a = np.loadtxt(path, usecols=(0,2)) # skip "*eV" column - except (ValueError, IndexError): - a = np.loadtxt(path, dtype=np.object) - pass - if os.path.exists(path+".npy"): - b = np.load(path+".npy") - assert np.all( a == b ) - log.debug(" path %s matches %s " % (path, path+".npy")) - pass - return a - - def __init__(self, specbase="PMTProperty.R12860.", symbol="hama"): - props = odict() - specs = [] - paths = [] - for name in self.Names(specbase): - spec = "%s%s" % ( specbase, name) - path = self.Resolve(spec) - props[name] = self.Load(path) - specs.append(spec) - paths.append(path) - pass - self.props = props - self.specs = specs - self.paths = paths - self.specbase = specbase - self.symbol = symbol - - def __repr__(self): - lines = [] - lines.append("Prop %s %s " % (self.specbase, self.symbol) ) - for k, v in self.props.items(): - lines.append("%35s : %10s : %10s " % ("%s.%s" % (self.symbol, k), str(v.shape), str(v.dtype) )) - pass - return "\n".join(lines) - - def __getattr__(self, spec): - if spec in self.props: - return self.props[spec] - else: - raise AttributeError - pass - -if __name__ == '__main__': - logging.basicConfig(level=logging.INFO) - - hama = SProp("PMTProperty.R12860.", symbol="hama") - nnvt = SProp("PMTProperty.NNVTMCP.", symbol="nnvt") - nnvtq = SProp("PMTProperty.NNVTMCP_HiQE.", symbol="nnvtq") - - print(repr(hama)) - print(repr(nnvt)) - print(repr(nnvtq)) - - for symbol in ["hama", "nnvt", "nnvtq"]: - expr = "np.c_[%(symbol)s.PHC_RINDEX, %(symbol)s.PHC_KINDEX, %(symbol)s.ARC_RINDEX ] " % locals() - print(expr) - print(eval(expr)) - pass - # energy domain consistent within but not between them - diff --git a/sysrap/__init__.py b/sysrap/__init__.py deleted file mode 100644 index a11c1d2a07..0000000000 --- a/sysrap/__init__.py +++ /dev/null @@ -1,19 +0,0 @@ -# -# Copyright (c) 2019 Opticks Team. All Rights Reserved. -# -# This file is part of Opticks -# (see https://bitbucket.org/simoncblyth/opticks). -# -# Licensed under the Apache License, Version 2.0 (the "License"); -# you may not use this file except in compliance with the License. -# You may obtain a copy of the License at -# -# http://www.apache.org/licenses/LICENSE-2.0 -# -# Unless required by applicable law or agreed to in writing, software -# distributed under the License is distributed on an "AS IS" BASIS, -# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. -# See the License for the specific language governing permissions and -# limitations under the License. -# - diff --git a/sysrap/dv.py b/sysrap/dv.py deleted file mode 100644 index 1a74a80e99..0000000000 --- a/sysrap/dv.py +++ /dev/null @@ -1,106 +0,0 @@ -#!/usr/bin/env python - -import os, numpy as np -from opticks.ana.fold import Fold -from opticks.ana.array_repr_mixin import ArrayReprMixin - -class AlignedDV(ArrayReprMixin, object): - """ - Simple deviation comparisons of random aligned arrays - - In [2]: pdv.dv - Out[2]: - array([[ 47, 117, 1732, 4412, 2710, 965, 16, 1, 0, 0], - [ 2746, 5430, 1724, 96, 4, 0, 0, 0, 0, 0], - [ 6404, 2937, 647, 11, 1, 0, 0, 0, 0, 0], - [ 9995, 1, 1, 3, 0, 0, 0, 0, 0, 0], - [10000, 0, 0, 0, 0, 0, 0, 0, 0, 0]], dtype=uint32) - - """ - EDGES = np.array( [0.,1e-6,1e-5,1e-4,1e-3, 1e-2, 1e-1, 1, 10, 100, 1000], dtype=np.float32 ) - COLUMN_LABEL = ["%g" % _ for _ in EDGES[1:]] - - POS = 0 - TIME = 1 - MOM = 2 - POL = 3 - WL = 4 - - ITEMS = (POS, TIME, MOM, POL, WL) - ROW_LABEL = ("pos", "time", "mom", "pol", "wl") - - def __init__(self, a, b, arr="photon", symbol=None): - """ - :param a: opticks.ana.fold Fold instance - :param b: opticks.ana.fold Fold instance - - In [12]: adv[...,0,3].shape - Out[12]: (10000,) - - In [13]: bdv[...,0,3].shape - Out[13]: (10000, 10) - """ - - wseq = np.where( a.seq[:,0] == b.seq[:,0] ) - adv = np.abs( getattr(a,arr)[wseq] - getattr(b,arr)[wseq] ) ## for deviations to be meaningful needs to be same history - - - ## although ... ellipsis indexing can generalize access between photon and record - ## need to use different np.amax axis args to keep the deviation at photon level, not record level - if arr == "photon": - time = adv[:,0,3] - pos = np.amax( adv[:,0,:3], axis=1 ) ## amax of the 3 position deviations, so can operate at photon position level, not x,y,z level - mom = np.amax( adv[:,1,:3], axis=1 ) - pol = np.amax( adv[:,2,:3], axis=1 ) - wl = adv[:,2,3] - elif arr == "record": - time = np.amax( adv[:,:,0,3], axis=1 ) - pos = np.amax( adv[:,:,0,:3], axis=(1,2) ) ## amax of the 3 position deviations, so can operate at photon position level, not x,y,z level - mom = np.amax( adv[:,:,1,:3], axis=(1,2) ) - pol = np.amax( adv[:,:,2,:3], axis=(1,2) ) - wl = np.amax( adv[:,:,2,3], axis=1 ) - else: - assert 0, "expecting photon or record not %s " % arr - pass - - dv = np.zeros( (len(self.ITEMS),len(self.EDGES)-1), dtype=np.uint32 ) - dv[self.POS], bins = np.histogram( pos, bins=self.EDGES ) - dv[self.TIME], bins = np.histogram( time, bins=self.EDGES ) - dv[self.MOM], bins = np.histogram( mom, bins=self.EDGES ) - dv[self.POL], bins = np.histogram( pol, bins=self.EDGES ) - dv[self.WL], bins = np.histogram( wl, bins=self.EDGES ) - - self.symbol = symbol if not symbol is None else "%s.dv" % arr - self.pos = pos - self.time = time - self.mom = mom - self.pol = pol - self.wl = wl - self.dv = dv - - def __repr__(self): - return self.MakeRepr(self.dv, symbol=self.symbol) - - -class HeaderAB(object): - def __init__(self, a, b, cmdline): - self.lines = ["A_FOLD : %s " % a.base, "B_FOLD : %s " % b.base, cmdline ] - def __repr__(self): - return "\n".join(self.lines) - - -if __name__ == '__main__': - a = Fold.Load("$A_FOLD", symbol="a") if "A_FOLD" in os.environ else None - b = Fold.Load("$B_FOLD", symbol="b") if "B_FOLD" in os.environ else None - - hdr = HeaderAB(a,b, "./dv.sh # cd ~/opticks/sysrap") - print(hdr) - - pdv = AlignedDV(a,b, arr="photon", symbol="pdv") - print(pdv) - - rdv = AlignedDV(a,b, arr="record", symbol="rdv") - print(rdv) - - - diff --git a/sysrap/sevt.py b/sysrap/sevt.py deleted file mode 100644 index 3472fd655b..0000000000 --- a/sysrap/sevt.py +++ /dev/null @@ -1,1211 +0,0 @@ -#!/usr/bin/env python -""" -sevt.py (formerly opticks.u4.tests.U4SimulateTest) -==================================================== - - -""" -import os, re, logging, textwrap, numpy as np -log = logging.getLogger(__name__) - -from opticks.ana.fold import Fold - -print("[from opticks.ana.p import * ") -from opticks.ana.p import * -print("]from opticks.ana.p import * ") - -from opticks.ana.eget import eslice_ -from opticks.ana.base import PhotonCodeFlags -from opticks.ana.qcf import QU,QCF,QCFZero -from opticks.ana.npmeta import NPMeta -from opticks.ana.hismask import HisMask - -hm = HisMask() - -pcf = PhotonCodeFlags() -fln = pcf.fln -fla = pcf.fla - - -MODE = int(os.environ.get("MODE","2")) -if MODE > 0: - from opticks.ana.pvplt import * -else: - pass -pass - -QLIM=eslice_("QLIM", "0:10") - - -axes = 0, 2 # X,Z -H,V = axes - - -class SEvt(object): - """ - Higher level wrapper for an Opticks SEvt folder of arrays - - WIP: Concatenation of multiple SEvt, starting from Fold concatenation - """ - @classmethod - def Load(cls, *args, **kwa): - NEVT = int(kwa.get("NEVT",0)) - log.info("SEvt.Load NEVT:%d " % NEVT) - if NEVT > 0: - f = cls.LoadConcat(*args,**kwa) - else: - f = Fold.Load(*args, **kwa ) - pass - return None if f is None else cls(f, **kwa) - - @classmethod - def LoadConcat(cls, *args, **kwa): - """ - HMM maybe should load the separate SEvt and then concat those, - rather than trying to concat at the lower Fold level - """ - NEVT = int(kwa.get("NEVT",0)) - assert NEVT > 0 - f = Fold.LoadConcat(*args, **kwa) - assert hasattr(f, 'ff') - return f - - def __init__(self, f, **kwa): - """ - :param f: Fold instance - :param kwa: override param, only W2M accepted - """ - - self.f = f - self.symbol = f.symbol - self._r = None - self.r = None - self._a = None - self.a = None - self._pid = -1 - self.W2M = kwa.get("W2M", None) - print("sevt.init W2M\n", self.W2M ) - - self.init_run(f) - self.init_meta(f) - self.init_fold_meta_timestamp(f) - self.init_fold_meta(f) - self.init_U4R_names(f) - self.init_photon(f) - self.init_record(f) - self.init_record_sframe(f) - self.init_SEventConfig_meta(f) - self.init_seq(f) - self.init_aux(f) - self.init_ee(f) ## handling of concatenated SEvt is done currently only for ee - self.init_aux_t(f) - self.init_sup(f) - self.init_junoSD_PMT_v2_SProfile(f) - self.init_env(f) # setting pid needs frames so has to be late - - @classmethod - def CommonRunBase(cls, f): - """ - :param f: Fold instance - :return urun_base: str - """ - run_bases = [] - if type(f.base) is str: - run_base = os.path.dirname(f.base) - run_bases.append(run_base) - elif type(f.base) is list: - for base in f.base: - run_base = os.path.dirname(base) - if not run_base in run_bases: - run_bases.append(run_base) - pass - else: - pass - pass - assert len(run_bases) == 1 - urun_base = run_bases[0] - return urun_base - - def init_run(self, f): - """ - HMM the run meta should be same for all of concatenated SEvts - """ - urun_base = self.CommonRunBase(f) - urun_path = os.path.join( urun_base, "run.npy" ) - - # THIS WAS RECURSING AND LOADING THE n001 and p001 folders - #print("[ sevt.init_run Fold.Load urun_path %s loading urun_base %s " % (urun_path,urun_base) ) - #fp = Fold.Load(urun_base) if os.path.exists(urun_path) else None - #print("] sevt.init_run Fold.Load ") - #run_meta = not fp is None and getattr(fp,'run_meta', None) != None - #self.fp = fp - - urun_meta = os.path.join( urun_base, "run_meta.txt" ) - run_meta = NPMeta.Load(urun_meta) if os.path.exists(urun_meta) else None - has_run_meta = not run_meta is None - - prof2stamp_ = lambda prof:prof.split(",")[0] - - rr_ = [prof2stamp_(run_meta.SEvt__BeginOfRun), - prof2stamp_(run_meta.SEvt__EndOfRun )] if has_run_meta else [0,0] - - rr = np.array(rr_, dtype=np.uint64 ) - - self.rr = rr - self.run_meta = run_meta - - qwns = ["FAKES_SKIP",] - for qwn in qwns: - mval = list(map(int,getattr(run_meta, qwn, ["-1"] ) if has_run_meta else ["-2"])) - setattr(self, qwn, mval[0] ) - pass - - - def init_env(self, f): - """ - """ - pid = int(os.environ.get("%sPID" % f.symbol.upper(), "-1")) - opt = os.environ.get("%sOPT" % f.symbol.upper(), "") - off_ = os.environ.get("%sOFF" % f.symbol.upper(), "0,0,0") - off = np.array(list(map(float,off_.split(",")))) - log.info("SEvt.__init__ symbol %s pid %d opt %s off %s " % (f.symbol, pid, opt, str(off)) ) - - self.pid = pid - self.opt = opt - self.off = off - - def init_meta(self, f): - """ - """ - metakey = os.environ.get("METAKEY", "junoSD_PMT_v2_Opticks_meta" ) - meta = getattr(f, metakey, None) - self.metakey = metakey - self.meta = meta - - def init_fold_meta_timestamp(self, f): - """ - """ - if getattr(f, "NPFold_meta", None) == None: return - msrc = f.NPFold_meta - if msrc is None: return - - def mts_(key): - val = msrc.get_value(key) - if val == '': val = "0" - return np.uint64(val) - - bor = mts_("T_BeginOfRun") - boe = mts_("t_BeginOfEvent") - gs0 = mts_("t_setGenstep_0") - gs1 = mts_("t_setGenstep_1") - gs2 = mts_("t_setGenstep_2") - gs3 = mts_("t_setGenstep_3") - gs4 = mts_("t_setGenstep_4") - gs5 = mts_("t_setGenstep_5") - gs6 = mts_("t_setGenstep_6") - gs7 = mts_("t_setGenstep_7") - gs8 = mts_("t_setGenstep_8") - pre = mts_("t_PreLaunch") - pos = mts_("t_PostLaunch") - eoe = mts_("t_EndOfEvent") - - b2e = eoe - boe - ee = np.array([boe, eoe, b2e], dtype=np.uint64 ) - ett = np.array([bor, boe, gs0, gs1, gs2, gs3, gs4, gs5, gs6, gs7, gs8, pre, pos, eoe], dtype=np.uint64 ) - ettl = "bor boe gs0 gs1 gs2 gs3 gs4 gs5 gs6 gs7 gs8 pre pos eoe".split() - ettd = np.zeros( len(ett), dtype=np.uint64 ) - ettd[1:] = np.diff(ett) - - ettv = ett.view("datetime64[us]") - ettc = np.c_[ettl, ett, ettv.astype("|S26"), ettd/1e6, ettd] - - self.ee = ee - self.ett = ett - self.ettd = ettd - self.ettv = ettv - self.ettl = ettl - self.ettc = ettc - - - def init_fold_meta(self, f): - """ - """ - if getattr(f, "NPFold_meta", None) == None: return - msrc = f.NPFold_meta - if msrc is None: return - - GEOM = msrc.get_value("GEOM") - CHECK = msrc.get_value("CHECK") - LAYOUT = msrc.get_value("LAYOUT") - VERSION = msrc.get_value("VERSION") - if LAYOUT.find(" ") > -1: LAYOUT = "" - - SCRIPT = os.environ.get("SCRIPT", "") - - GEOMList = msrc.get_value("${GEOM}_GEOMList", "") - if GEOMList.endswith("GEOMList"): GEOMList = "" - - TITLE = "N=%s %s # %s/%s " % (VERSION, SCRIPT, LAYOUT, CHECK ) - IDE = list(filter(None,[f.symbol.upper(), GEOM, GEOMList, "N%s"%VERSION, LAYOUT, CHECK])) - ID = "_".join(IDE) - - self.CHECK = CHECK - self.LAYOUT = LAYOUT - self.VERSION = VERSION - self.SCRIPT = SCRIPT - self.GEOM = GEOM - self.GEOMList = GEOMList - self.TITLE = TITLE - self.IDE = IDE - self.ID = ID - - def init_U4R_names(self, f): - """ - Now its an array with UNSET - """ - U4R_names = getattr(f, "U4R_names", None) - #SPECS = np.array(U4R_names.lines) if not U4R_names is None else None - SPECS = U4R_names if not U4R_names is None else None - self.SPECS = SPECS - - @classmethod - def Make_flagmask_table(cls, flagmask): - """ - - In [30]: np.c_[ label, fmtab[:,1] ][:20] - Out[30]: - array([['TO|BT|SA', '326185'], - ['EC|TO|BT|SD', '312177'], - ['TO|BT|SR|SA', '56805'], - ['TO|BT|BR|SR|SA', '53334'], - ['TO|BT|BR|SA', '42888'], - ['TO|BT|BR|AB', '37716'], - ['EC|TO|BT|BR|SD', '21709'], - ['TO|BT|BR|SA|SC', '20683'], - ['TO|BT|BR|SC|AB', '15776'], - ['EC|TO|BT|BR|SD|SC', '13679'], - ['TO|BT|AB', '12490'], - ['TO|BT|BR|SR|SA|SC', '8668'], - - """ - _fmtabraw = np.c_[np.unique(flagmask,return_counts=True)] - fmtabraw = _fmtabraw[np.argsort(_fmtabraw[:,1])][::-1] - label = hm.label(fmtabraw[:,0]) - fmtab = np.c_[ label, fmtabraw[:,1] ] - return fmtab - - def init_photon(self, f): - if hasattr(f,'photon') and not f.photon is None: - iix = f.photon[:,1,3].view(np.int32) - flagmask = f.photon[:,3,3].view(np.int32) - fmtab = self.Make_flagmask_table(flagmask) - else: - iix = None - flagmask = None - fmtab = None - pass - self.iix = iix - self.flagmask = flagmask - self.fmtab = fmtab - - - def init_record(self, f): - """ - :param f: fold instance - """ - if hasattr(f,'record') and not f.record is None: - - if f.record.dtype.name == 'float32': - iir = f.record[:,:,1,3].view(np.int32) - idr = f.record[:,:,1,3].view(np.int32) - elif f.record.dtype.name == 'float64': - iir = f.record[:,:,1,3].view(np.int64) - idr = f.record[:,:,1,3].view(np.int64) - else: - assert False - pass - else: - iir = None - idr = None - pass - self.iir = iir - self.idr = idr - - def init_record_sframe(self, f): - - with_sframe = not getattr(f, "sframe", None) is None - with_record = not getattr(f, "record", None) is None - - W2M = self.W2M - - if with_sframe and with_record: - gpos = np.ones( f.record.shape[:-1] ) ## trim last dimension reducing eg (10000,32,4,4) to (10000, 32, 4) - gpos[:,:,:3] = f.record[:,:,0,:3] ## point positions of all photons - w2m = W2M if not W2M is None else f.sframe.w2m - lpos = np.dot( gpos, w2m ) ## transform all points from global to local frame - else: - w2m = W2M if not W2M is None else None - gpos = None - lpos = None - pass - self.w2m = w2m - self.gpos = gpos - self.lpos = lpos - - - - def init_SEventConfig_meta(self, f): - meta = getattr(f, 'SEventConfig_meta', {}) - - ipl = getattr(meta, "InputPhoton", []) - ip = ipl[0] if len(ipl)==1 else None - - ipfl = getattr(meta, "InputPhotonFrame", []) - ipf = ipfl[0] if len(ipfl)==1 else None - - ipcl = getattr(meta, "InputPhotonCheck", []) - ipc = ipc[0] if len(ipcl)==1 else None - - self.ip = ip - self.ipf = ipf - self.ipc = ipc - - - def q_find(self, txt="EC", encoding="utf-8"): - """ - :param txt: - :return a: array of indices of self.q array elements that contain *txt* - - np.char.find returns -1 when the txt is not found in array elements - so adding one makes that zero which is excluded with np.flatnonzero - - Note that because the self.q array elements typically end with blanks - it is not useful to implement q_endswith, instead use this q_find - to find array indices with elements containing the txt. - """ - return np.flatnonzero(1+np.char.find(self.q, txt.encode(encoding) )) - - def q_find_or(self, txt1="EC", txt2="EX", encoding="utf-8"): - return np.flatnonzero(np.logical_or( - 1+np.char.find(self.q,txt1.encode(encoding)), - 1+np.char.find(self.q,txt2.encode(encoding)) - )) - - def q_find_or_(self, txts_="EC,EX", delim=",", encoding="utf-8"): - txts= txts_.split(delim) - aa = [] - for txt in txts: - aa.append( 1+np.char.find(self.q, txt.encode(encoding))) - pass - return np.flatnonzero(np.logical_or.reduce(aa)) - - def q__(self, prefix="TO BT SD", encoding="utf-8"): - return self.q_startswith_or_(prefix) if "," in prefix else self.q_startswith(prefix=prefix) - - def q_startswith(self, prefix="TO BT SD", encoding="utf-8"): - return np.flatnonzero(np.char.startswith(self.q, prefix.encode(encoding) )) - - def q_startswith_or(self, pfx1="TO BT BT BT SA", pfx2="TO BT BT BT SD", encoding="utf-8"): - return np.flatnonzero(np.logical_or( - np.char.startswith(self.q,pfx1.encode(encoding)), - np.char.startswith(self.q,pfx2.encode(encoding)) - )) - - def q_startswith_or_(self, pfxs_="TO BT BT BT SA,TO BT BT BT SD", delim=",", encoding="utf-8"): - """ - :param pfxs_: delimited string list with potentially multiple history strings - :param delim: - - Another way to implement this is with np.any - """ - pfxs = pfxs_.split(delim) - aa = [] - for pfx in pfxs: - aa.append( np.char.startswith(self.q, pfx.encode(encoding))) - pass - return np.flatnonzero(np.logical_or.reduce(aa)) - - def init_seq(self, f): - """ - - +---------------------+----------------------------------------------------------+ - | a.f.seq.shape | (1000000, 2, 2) | - +---------------------+----------------------------------------------------------+ - | a.f.seq[:,0].shape | (1000000, 2) | - +---------------------+----------------------------------------------------------+ - | a.f.seq[0,0] | array([14748335283153718477, 37762157772], dtype=uint64) | - +---------------------+----------------------------------------------------------+ - - """ - symbol = f.symbol - qlim = QLIM - qtab_ = "np.c_[qn,qi,qu][quo][qlim]" - - if hasattr(f,'seq') and not f.seq is None: - q_ = f.seq[:,0] - q = ht.seqhis(q_) # ht from opticks.ana.p - qq = ht.Convert(q_) # (n,32) int8 : for easy access to nibbles - n = np.sum( seqnib_(q_), axis=1 ) # occupied nibbles, ie number of history points - - qu, qi, qn = np.unique(q, return_index=True, return_counts=True) - quo = np.argsort(qn)[::-1] - - qtab = eval(qtab_) - qtab_ = qtab_.replace("q","%s.q" % symbol) - - nosc = np.ones(len(qq), dtype=np.bool_ ) # start all true - nosc[np.where(qq == pcf.SC)[0]] = 0 # knock out photons with scatter in their histories - - nore = np.ones(len(qq), dtype=np.bool_ ) # start all true - nore[np.where(qq == pcf.RE)[0]] = 0 # knock out photons with reemission in their histories - - noscab = np.ones(len(qq), dtype=np.bool_ ) # start all true - noscab[np.where(qq == pcf.SC)[0]] = 0 # knock out photons with scatter in their histories - noscab[np.where(qq == pcf.AB)[0]] = 0 # knock out photons with bulk absorb in their histories - else: - q_ = None - q = None - qq = None - n = None - qu, qi, qn = None, None, None - quo = None - qtab = None - - nosc = None - nore = None - noscab = None - pass - - self.q_ = q_ - self.q = q - self.qq = qq - self.n = n - - self.qu = qu - self.qi = qi - self.qn = qn - - self.quo = quo - self.qtab = qtab - - self.qlim = qlim - self.qtab_ = qtab_ - - self.nosc = nosc # mask of photons without SC in their histories - self.nore = nore # mask of photons without RE in their histories - self.noscab = noscab # mask of photons without SC or AB in their histories - - - - def minimal_qtab(self, sli="[:10]", dump=False): - """ - :return qtab: history table in descending count order - """ - e = self - uq,iq,nq = np.unique(e.q, return_index=True, return_counts=True) - oq = np.argsort(nq)[::-1] - expr ="np.c_[nq,iq,uq][oq]%(sli)s" % locals() - qtab = eval(expr) - if dump: - log.info("minimal_qtab : %s " % expr) - print(qtab) - pass - return qtab - - def init_aux(self, f): - """ - Note aux is currently CPU only - """ - if hasattr(f, 'aux') and not f.aux is None: - fk = f.aux[:,:,2,2].view(np.uint32) ## fakemask : for investigating fakes when FAKES_SKIP is disabled - spec_ = f.aux[:,:,2,3].view(np.int32) ## step spec - max_point = f.aux.shape[1] # instead of hardcoding typical values like 32 or 10, use array shape - uc4 = f.aux[:,:,2,2].copy().view(np.uint8).reshape(-1,max_point,4) ## see sysrap/spho.h c4/C4Pho.h - eph = uc4[:,:,1] # .y ProcessHits enum at step point level - ep = np.max(eph, axis=1 ) # ProcessHits enum at photon level - else: - fk = None - spec_ = None - uc4 = None - eph = None - ep = None - t = None - pass - self.fk = fk - self.spec_ = spec_ - self.uc4 = uc4 - self.eph = eph # ProcessHits EPH enum at step point level - self.ep = ep # ProcessHits EPH enum at photon level - - def init_ee(self, f): - """ - microsecond[us] timestamps (UTC epoch) [BeginOfEvent,EndOfEvent,Begin2End] - - For concatenated SEvt the total of all the EndOfEvent-BeginOfEvent - is placed into ee[-1] - - TIMING METADATA HAS BEEN MOVED FROM THE PHOTON TO THE FOLD - with_ff = not getattr(f, 'ff', None) is None - log.info("init_ee with_ff:%d" % (with_ff)) - if with_ff: - kk = f.ff.keys() - boe = np.uint64(0) - eoe = np.uint64(0) - b2e = np.uint64(0) - for k in kk: - fk = f.ff[k] - k_boe = np.uint64(fk.photon_meta.t_BeginOfEvent[0]) - k_eoe = np.uint64(fk.photon_meta.t_EndOfEvent[0]) - k_b2e = np.uint64(k_eoe-k_boe) - b2e += k_b2e - pass - else: - boe, eoe, b2e = 0,0,0 - pass - self.ee = np.array([boe, eoe, b2e], dtype=np.uint64 ) - """ - - def init_junoSD_PMT_v2_SProfile(self, f): - """ - The timestamps come from sysrap/sstamp.h and are datetime64[us] (UTC) compliant - - pf - uint64_t microsecond timestamps collected by SProfile.h - pfr - uint64_t (last - first) timestamp difference for each SProfile (currently ProcessHits call) - - - More than 10% of time spent in ProcessHits:: - - In [16]: np.sum(a.pfr)/a.ee[-1] - Out[16]: 0.10450881266649384 - - In [17]: np.sum(b.pfr)/b.ee[-1] - Out[17]: 0.11134881257006096 - - """ - if hasattr(f, 'junoSD_PMT_v2_SProfile') and not f.junoSD_PMT_v2_SProfile is None: - pf = f.junoSD_PMT_v2_SProfile - pfmx = np.max(pf[:,1:], axis=1 ) - pfmi = pf[:,1] - pfr = pfmx - pfmi - else: - pf = None - pfmx = None - pfmi = None - pfr = None - pass - self.pf = pf ## CAUTION: multiple ProcessHits calls per photon, so not in photon index order - self.pfmx = pfmx - self.pfmi = pfmi - self.pfr = pfr - - - - def init_aux_t(self, f): - if hasattr(f, 'aux') and not f.aux is None: - t = f.aux[:,:,3,:2].copy().view(np.uint64).squeeze() # step point timestamps - else: - t = None - pass - self.t = t # array of photon step point time stamps (UTC epoch) - - def init_sup(self, f): - """ - sup is CPU only - - photon level beginPhoton endPhoton time stamps from the sup quad4 - """ - if hasattr(f,'sup') and not f.sup is None: - s0 = f.sup[:,0,:2].copy().view(np.uint64).squeeze() # SEvt::beginPhoton (xsup.q0.w.x) - s1 = f.sup[:,0,2:].copy().view(np.uint64).squeeze() # SEvt::finalPhoton (xsup.q0.w.y) - ss = s1 - s0 # endPhoton - beginPhoton - - f0 = f.sup[:,1,:2].copy().view(np.uint64).squeeze() # SEvt::finalPhoton (xsup.q1.w.x) t_PenultimatePoint - f1 = f.sup[:,1,2:].copy().view(np.uint64).squeeze() # SEvt::finalPhoton (xsup.q1.w.y) t_LastPoint - ff = f1 - f0 # LastPoint - PenultimatePoint - - h0 = f.sup[:,2,:2].copy().view(np.uint64).squeeze() # SEvt::addProcessHitsStamp(0) (xsup.q2.w.x) - h1 = f.sup[:,2,2:].copy().view(np.uint64).squeeze() # SEvt::addProcessHitsStamp(0) (xsup.q2.w.y) - hh = h1 - h0 # timestamp range of SEvt::AddProcessHitsStamp(0) calls - hc = f.sup[:,3,0].view(np.int32) - - i0 = f.sup[:,4,:2].copy().view(np.uint64).squeeze() # SEvt::addProcessHitsStamp(1) (xsup.q4.w.x) - i1 = f.sup[:,4,2:].copy().view(np.uint64).squeeze() # SEvt::addProcessHitsStamp(1) (xsup.q4.w.y) - ii = i1 - i0 # timestamp range of SEvt::AddProcessHitsStamp(1) calls - ic = f.sup[:,5,0].view(np.int32) - - hi0 = i0 - h0 - hi1 = i1 - h1 - else: - s0 = None - s1 = None - ss = None - - f0 = None - f1 = None - ff = None - - h0 = None - h1 = None - hh = None - hc = None - - i0 = None - i1 = None - ii = None - ic = None - - hi0 = None - hi1 = None - pass - if not getattr(self, 't', None) is None and not s0 is None: - t0 = s0.min() - tt = self.t.copy() - tt[np.where( tt != 0 )] -= t0 # subtract earliest time "pedestal", but leave the zeros - else: - t0 = None - tt = None - pass - - self.t0 = t0 # scalar : event minimum time stamp (which is minimum s0:beginPhoton timestamp) - self.tt = tt # array of photon step point time stamps relative to t0 - - self.s0 = s0 # array of beginPhoton time stamps (UTC epoch) - self.s1 = s1 # array of endPhoton time stamps (UTC epoch) - self.ss = ss # array of endPhoton-beginPhoton timestamp differences - - self.f0 = f0 # array of PenultimatePoint timestamps (UTC epoch) - self.f1 = f1 # array of LastPoint timestamps (UTC epoch) - self.ff = ff # array of LastPoint-PenultimatePoint timestamp differences - - self.h0 = h0 # array of SEvt::AddProcessHitsStamp(0) range begin (UTC epoch) - self.h1 = h1 # array of SEvt::AddProcessHitsStamp(0) range end (UTC epoch) - self.hh = hh # array of SEvt::AddProcessHitsStamp(0) ranges (microseconds [us]) - self.hc = hc # array of SEvt::AddProcessHitsStamp(0) call counts for each photon - - self.i0 = i0 # array of SEvt::AddProcessHitsStamp(1) range begin (UTC epoch) - self.i1 = i1 # array of SEvt::AddProcessHitsStamp(1) range end (UTC epoch) - self.ii = ii # array of SEvt::AddProcessHitsStamp(1) ranges (microseconds [us]) - self.ic = ic # array of SEvt::AddProcessHitsStamp(1) call counts for each photon - - self.hi0 = hi0 # array of range begin SEvt::AddProcessHitsStamp(1)-SEvt::AddProcessHitsStamp(0) - self.hi1 = hi1 # array of range end SEvt::AddProcessHitsStamp(1)-SEvt::AddProcessHitsStamp(0) - - def __repr__(self): - fmt = "SEvt symbol %s pid %s opt %s off %s %s.f.base %s " - return fmt % ( self.symbol, self.pid, self.opt, str(self.off), self.symbol, self.f.base ) - - def _get_pid(self): - return self._pid - def _set_pid(self, pid): - """ - """ - f = self.f - q = self.q - W2M = self.W2M - symbol = self.f.symbol - if pid > -1 and hasattr(f,'record') and pid < len(f.record): - _r = f.record[pid] - wl = _r[:,2,3] - r = _r[wl > 0] ## select wl>0 to avoid lots of record buffer zeros, example shape (13, 4, 4) - - g = np.ones([len(r),4]) - g[:,:3] = r[:,0,:3] - - w2m = W2M if not W2M is None else f.sframe.w2m - l = np.dot( g, w2m ) - ld = np.diff(l,axis=0) - - pass - _aux = getattr(f, 'aux', None) - _a = None if _aux is None else _aux[pid] - if _a is None: - a = None - ast = "no-ast" - else: - a = _a[wl > 0] - ast = a[:len(a),1,3].copy().view(np.int8)[::4] - pass - qpid = q[pid][0].decode("utf-8").strip() - - npoi = (len(qpid)+1)//3 - qkey = " ".join(map(lambda _:"%0.2d"%_, range(npoi))) - - #label = "%s %sPID:%d\n%s " % (qpid, symbol.upper(), pid, qkey) - label = "%s\n%s" % (qpid, qkey) - else: - _r = None - r = None - _a = None - a = None - ast = None - qpid = None - label = "" - - g = None - l = None - ld = None - pass - self._pid = pid - self._r = _r - self.r = r - self._a = _a - self.a = a - self.ast = ast - self.qpid = qpid - self.label = label - self.g = g - self.l = l - self.ld = ld - - pid = property(_get_pid, _set_pid) - - - #def _get_w2m(self): - # return self._w2m - # - #def _set_w2m(self, w2m=None): - # self._w2m = w2m if not w2m is None else f.sframe.w2m - - - - - def spec_(self, i): - """ - ## need np.abs as detected fakes that are not skipped are negated - In [10]: np.c_[a.spec[5,:a.n[5]],a.SPECS[np.abs(a.spec[5,:a.n[5]])]] - Out[10]: - array([['0', 'UNSET'], - ['1', 'Water/Water:pInnerWater/pLPMT_Hamamatsu_R12860'], - ['2', 'Water/AcrylicMask:pLPMT_Hamamatsu_R12860/HamamatsuR12860pMask'], - ['3', 'AcrylicMask/Water:HamamatsuR12860pMask/pLPMT_Hamamatsu_R12860'], - ['4', 'Water/Pyrex:pLPMT_Hamamatsu_R12860/HamamatsuR12860_PMT_20inch_log_phys'], - ['-5', 'Pyrex/Pyrex:HamamatsuR12860_PMT_20inch_log_phys/HamamatsuR12860_PMT_20inch_body_phys'], - ['6', 'Pyrex/Pyrex:HamamatsuR12860_PMT_20inch_body_phys/HamamatsuR12860_PMT_20inch_body_phys']], dtype=' 0 ) - - for i in range(NEVT): - - afold = AFOLD % i - bfold = BFOLD % i - a = SEvt.Load(afold,symbol="a", quiet=True) - b = SEvt.Load(bfold,symbol="b", quiet=True) - if a is None: - log.fatal("FAILED TO LOAD SEVT A FROM AFOLD: %s " % afold ) - pass - if b is None: - log.fatal("FAILED TO LOAD SEVT B FROM BFOLD: %s " % bfold ) - pass - - ab = SAB(a,b) - mab[i] = ab - - - kk = ab.meta.kk - #print("ab.meta.kk (A,B common keys):%s" % str(kk)) - - skk = ab.meta.skk - tab = ab.meta.tab - - # metadata fields that look like floats, integers, timestamps - ffield = np.unique(np.where(np.char.find(tab,'.') > -1 )[0]) - ifield = np.unique(np.where( np.logical_and( np.char.str_len( tab ) < 15, np.char.find(tab, '.') == -1 ))[0]) - tfield = np.unique(np.where( np.logical_and( np.char.str_len( tab ) > 15, np.char.find(tab, '.') == -1 ))[0]) - - if first: - first = False - skk0 = skk - kk0 = kk - tab0 = tab - ffield0 = ffield - ifield0 = ifield - tfield0 = tfield - else: - assert np.all( skk == skk0 ) - assert kk == kk0 - assert tab.shape == tab0.shape - assert np.all( ffield == ffield0 ) - assert np.all( ifield == ifield0 ) - assert np.all( tfield == tfield0 ) - pass - - itab = tab[ifield] - ftab = tab[ffield] - ttab = tab[tfield] - - itabs.append(itab) - ftabs.append(ftab) - ttabs.append(ttab) - - if not ab.qcf is None: - c2tab[0,i] = ab.qcf.c2sum - c2tab[1,i] = ab.qcf.c2n - c2tab[2,i] = ab.qcf.c2per - pass - - # hmm how to present together with c2sum, c2n, c2per ? - if "c2desc" in os.environ: - fmt = " %s : %%s " % efmt - print( fmt % ( i, ab.qcf.c2desc)) - else: - if "DUMP" in os.environ:print(repr(ab)) - pass - pass - - - self.c2tab = c2tab - self.c2tab_zero = np.all( c2tab == 0. ) - - self.mab = mab - akk = np.array( list(kk0)) - ikk = akk[ifield] - fkk = akk[ffield] - tkk = akk[tfield] - - ik = KeyIdx(ikk, symbol="%s.ik" % self.symbol ) - fk = KeyIdx(fkk, symbol="%s.fk" % self.symbol ) - tk = KeyIdx(tkk, symbol="%s.tk" % self.symbol ) - - c_itab = np.zeros( (itab.shape[0], itab.shape[1], NEVT ), dtype=np.uint64 ) - for i in range(len(itabs)): - c_itab[:,:,i] = itabs[i] - pass - c_ftab = np.zeros( (ftab.shape[0], ftab.shape[1], NEVT ), dtype=np.float64 ) - for i in range(len(ftabs)): - c_ftab[:,:,i] = ftabs[i] - pass - c_ttab = np.zeros( (ttab.shape[0], ttab.shape[1], NEVT ), dtype=np.uint64 ) - for i in range(len(ttabs)): - c_ttab[:,:,i] = ttabs[i] - pass - - self.ikk = ikk - self.fkk = fkk - self.tkk = tkk - - self.ik = ik - self.fk = fk - self.tk = tk - - self.c_itab = c_itab - self.c_ftab = c_ftab - self.c_ttab = c_ttab - - - - def annotab(self, tabsym, keys, nline=3): - """ - :param tabsym: attribute symbol of table - :param keys: array of names for each table item - :param nline: number of lines for each table item - :return str: annoted table string - - Annotate a numpy array repr with keys - TODO: split off into reusable AnnoTab? object - """ - expr = "%%(sym)s.%s" % tabsym - lines = [] - lines.append("\n%s\n" % ( expr % {'sym':self.symbol} )) - rawlines = repr(eval(expr % {'sym':"self" })).split("\n") - for i, line in enumerate(rawlines): - j = i//nline - anno = "%2d:%s" % (j,keys[j]) if i % nline == 0 else "" - lines.append("%s %s " % (line, anno )) - pass - return "\n".join(lines) - - def __repr__(self): - lines = [] - lines.append(self.annotab("c_itab", self.ikk, 3)) - lines.append(self.annotab("c_ftab", self.fkk, 3)) - pass - fmt_ptn = re.compile("FMT:(.*)\\s*") - - for expr in self.exprs: - fmt_match = fmt_ptn.search(expr) - fmt = fmt_match.groups()[0].replace("%%","%") if not fmt_match is None else None - if "c2tab" in expr and self.c2tab_zero: continue - lines.append("\n%s\n" % ( expr % {'sym':self.symbol} )) - value = eval( expr % {'sym':"self"} ) - urep = repr(value) if fmt is None else fmt % value - lines.append(urep) - pass - return "\n".join(lines) - - - -if __name__ == '__main__': - logging.basicConfig(level=logging.INFO) - - FOLD = os.environ.get("FOLD", None) - log.info(" -- SEvt.Load FOLD" ) - a = SEvt.Load(FOLD, symbol="a") # optional photon histories - print(a) - - beg_ = "a.f.record[:,0,0,:3]" - beg = eval(beg_) - - end_ = "a.f.photon[:,0,:3]" - end = eval(end_) - - - #label0, ppos0 = None, None - label0, ppos0 = "b:%s" % beg_ , beg - - #label0, ppos0 = None, None - label1, ppos1 = "r:%s" % end_ , end - - - HEADLINE = "%s %s" % ( a.LAYOUT, a.CHECK ) - label = "\n".join( filter(None, [HEADLINE, label0, label1])) - print(label) - - if MODE == 0: - print("not plotting as MODE 0 in environ") - elif MODE == 2: - fig, ax = mpplt_plotter(label=label) - - ax.set_ylim(-250,250) - ax.set_xlim(-500,500) - - if not ppos0 is None: ax.scatter( ppos0[:,H], ppos0[:,V], s=1, c="b" ) - if not ppos1 is None: ax.scatter( ppos1[:,H], ppos1[:,V], s=1, c="r" ) - - fig.show() - elif MODE == 3: - pl = pvplt_plotter(label) - os.environ["EYE"] = "0,100,165" - os.environ["LOOK"] = "0,0,165" - pvplt_viewpoint(pl) - pl.add_points(ppos0) - pl.show() - pass -pass - diff --git a/sysrap/sevt_load.py b/sysrap/sevt_load.py deleted file mode 100644 index ea717c3811..0000000000 --- a/sysrap/sevt_load.py +++ /dev/null @@ -1,17 +0,0 @@ -#!/usr/bin/env python - -import os, numpy as np -log = logging.getLogger(__name__) -from opticks.sysrap.sevt import SEvt - -NEVT = int(os.environ.get("NEVT", 0)) # when NEVT>0 SEvt.LoadConcat loads and concatenates the SEvt - - -if __name__ == '__main__': - logging.basicConfig(level=logging.INFO) - log.info("FOLD:%s" % os.environ.get("FOLD", "-")) - t = SEvt.Load(symbol="t", NEVT=NEVT) - print(repr(t)) - - - diff --git a/sysrap/sevt_tt.py b/sysrap/sevt_tt.py deleted file mode 100644 index 70322cbdb5..0000000000 --- a/sysrap/sevt_tt.py +++ /dev/null @@ -1,576 +0,0 @@ -#!/usr/bin/env python -""" -sevt_tt.py (formerly u4/tests/U4SimulateTest_tt.py) -====================================================== - -Python script (not module) usage from two bash scripts --------------------------------------------------------- - -:: - - PLOT=STAMP STAMP_ANNO=1 STAMP_LEGEND=1 ~/opticks/u4/tests/tt.sh - - PLOT=STAMP STAMP_TT=90000,5000 STAMP_ANNO=1 ~/j/ntds/ntds.sh tt - - PLOT=PHO_AVG ~/j/ntds/ntds.sh tt - - -PLOT envvar selects what to plot ----------------------------------- - -PLOT=STAMP - time stamp illustration - - STAMP_TT=200000,5000 - STAMP_TT=200k,5k - control time window of STAMP plot in microseconds [us] - - STAMP_ANNO=1 - enable photon index and history annotation - -PLOT=PHO_AVG - average beginPhoton->endPhoton CPU time vs photon point count - -PLOT=PHO_HIS - histogram of beginPhoton->endPhoton CPU times - - PLOT=PHO_HIS ~/j/ntds/ntds.sh tt - -PLOT=PHO_N - A, B photon step point counts (fakes must be skipped in A to get a match) - - PLOT=PHO_N ~/j/ntds/ntds.sh tt - - -Dev Notes ------------- - -HMM: for reusability, its tempting to -just relocate this into sysrap/tests/tt.py -as a runnable python script which is used from -wherever with bash level envvars controlling the events -to be loaded rather than doing python dev to -bring in functionality at python level. - -When is this bash in control approach appropriate -as opposed to moving into sevt.py ? - -The level of generality determines what is appropriate -If the functionality really is general and likely -to be usable from all over the place then it belongs -into sevt.py with python control and flexibility of use. - -However if only likely to use functionality from a few places -pointing at differnt events that doing -at bash level seems appropriate and avoids python -module development. - -:: - - u4t - ./U4SimulateTest.sh tt - - PLOT=STAMP ~/opticks/u4/tests/tt.sh - PLOT=PHO_HIS ~/opticks/u4/tests/tt.sh - PLOT=PHO_AVG ~/opticks/u4/tests/tt.sh - - -""" -import os, numpy as np, textwrap, logging -log = logging.getLogger(__name__) -from opticks.ana.fold import Fold -from opticks.ana.p import * -from opticks.ana.eget import efloatarray_ -from opticks.sysrap.sevt import SEvt - - -SCRIPT = "./U4SimulateTest.sh tt" -os.environ["SCRIPT"] = SCRIPT -LABEL = os.environ.get("LABEL", "U4SimulateTest_tt.py" ) - - -N = int(os.environ.get("N", "-1")) -MODE = int(os.environ.get("MODE", "2")) -assert MODE in [0,2,-2,3,-3] -FLIP = int(os.environ.get("FLIP", "0")) == 1 -TIGHT = int(os.environ.get("TIGHT", "0")) == 1 -STAMP_TT = efloatarray_("STAMP_TT", "162307,3000") # "t0,dt" microseconds 1M=1s - -PLOT = os.environ.get("PLOT", None) -NEVT = int(os.environ.get("NEVT", 0)) # when NEVT>0 SEvt.LoadConcat loads and concatenates the SEvt - -ENVOUT = os.environ.get("ENVOUT", None) -CAP_STEM = PLOT -CAP_BASE = None # set below to a.f.base or b.f.base - -if MODE != 0: - from opticks.ana.pvplt import * - WM = mp.pyplot.get_current_fig_manager() -else: - WM = None -pass - - - -def GetAnnotation(title, EXPRS_, syms): - EXPRS = list(filter(None, textwrap.dedent(EXPRS_).split("\n"))) - - rlines = [] - rlines.append(title) - - for i,sym in enumerate(syms): - for expr_ in EXPRS: - label = expr_[expr_.find("#"):] if expr_.find("#") > -1 else "" - expr_ = expr_.split("#")[0].strip() - expr = expr_ % locals() - print(expr) - try: - val = eval(expr) - except ValueError: - log.fatal("FAILED EVAL : %s " % expr ) - val = "?" - pass - fmt_ = {} - fmt_["str"] = "%30s : %s : %s " - fmt_["int"] = "%30s : %8d : %s " - fmt_["float"] = "%30s : %8.3f : %s " - fmt_[""] = fmt_["float"] - - if type(val) is list: val = " ".join(val) - fmt = fmt_.get(type(val).__name__, fmt_["float"] ) - rlines.append(fmt % ( expr, val, label )) - pass - rlines.append("") - pass - ranno = "\n".join(rlines) - return ranno - - - -if __name__ == '__main__': - logging.basicConfig(level=logging.INFO) - - print("AFOLD:%s" % os.environ.get("AFOLD", "-")) - print("BFOLD:%s" % os.environ.get("BFOLD", "-")) - print("N:%d " % N ) - - if N == -1: - a = SEvt.Load("$AFOLD",symbol="a", NEVT=NEVT) - b = SEvt.Load("$BFOLD",symbol="b", NEVT=NEVT) - syms = ['a','b'] - evts = [a,b] - elif N == 0: - a = SEvt.Load("$AFOLD",symbol="a", NEVT=NEVT) - b = None - syms = ['a'] - evts = [a,] - elif N == 1: - a = None - b = SEvt.Load("$BFOLD",symbol="b", NEVT=NEVT) - syms = ['b'] - evts = [b,] - else: - assert(0) - pass - - ## CAP_BASE is passed via ENVOUT to the invoking bash script - ## to control where the figs folder with screen captures should be - if not a is None: - CAP_BASE = a.f.base - elif not b is None: - CAP_BASE = b.f.base - pass - - # concatenated SEvt have list of bases - if type(CAP_BASE) is list: CAP_BASE = CAP_BASE[0] - - if not a is None:print(repr(a)) - if not b is None:print(repr(b)) - - - if PLOT == "EXPRS": - EXPRS = r""" - %(sym)s - %(sym)s.symbol - w - %(sym)s.q[w] - %(sym)s.n[w] - np.diff(%(sym)s.tt[w]) - np.c_[%(sym)s.t[w].view("datetime64[us]")] - """ - - for sym in syms: - - n = eval("%(sym)s.n" % locals() ) - ww = np.where( n > 24)[0] - - for w in ww: - for e_ in textwrap.dedent(EXPRS).split("\n"): - e = e_ % locals() - print(e) - if len(e) == 0 or e[0] == "#": continue - print(eval(e)) - pass - pass - pass - pass - - - # needs: syms, a, b - - if PLOT == "PHO_HIS": - msg = "histogram of beginPhoton->endPhoton durations" - label = "%s : %s" % (PLOT, msg) - fig, axs = mpplt_plotter(nrows=1, ncols=1, label=label, equal=False) - ax = axs[0] - e_ = "%(sym)s.s1 - %(sym)s.s0 # PHO_HIS " - - #bins = np.logspace( np.log10(0.1),np.log10(500.0), 25 ) - bins = np.linspace( 0, 500, 50 ) - - h = {} - for sym in syms: - e = e_ % locals() - h[sym] = np.histogram( eval(e) , bins=bins ) - ax.plot( h[sym][1][:-1], h[sym][0], label=e ) - pass - ax.legend() - fig.show() - pass - - - if PLOT == "TABLE": - - HEAD = "np.c_[" - TAIL = "]" - FIELDS = list(filter(None,textwrap.dedent(r""" - %(sym)s.s1[w] - %(sym)s.s0[w] # beginPhoton->endPhoton - %(sym)s.t[w,0] - %(sym)s.s0[w] # beginPhoton->firstPoint - %(sym)s.t[w,%(n_1)s] - %(sym)s.t[w,0] # firstPoint->lastPoint - %(sym)s.s1[w] - %(sym)s.t[w,%(n_1)s] # lastPoint->endPhoton - """).split("\n"))) - - LABELS = list(map(lambda _:_[_.find("#")+1:], FIELDS)) - FIELDS = list(map(lambda _:_[:_.find("#")].strip(), FIELDS)) - - print("compare first and last point stamp range with beginPhoton endPhoton range") - - nn = np.arange(2,33) - for n in nn: - n_1 = int(n - 1) - for sym in syms: - w_ = "np.where( %(sym)s.n == %(n)s )[0]" % locals() - w = eval(w_) - expr__ = "".join([HEAD, ",".join(FIELDS), TAIL ]) - expr_ = expr__ % locals() - expr = eval(expr_) - label = " ".join(LABELS) - if len(expr) == 0: continue - - print(expr_) - print(w_) - print(label) - print(expr) - pass - pass - pass - - - print(" average photonBegin->photonEnd for different step counts ") - ssa = np.zeros((33,len(syms),4), dtype=np.float64 ) - for n in range(33): - for i,sym in enumerate(syms): - w_ = "np.where( %(sym)s.n == %(n)s )[0]" % locals() - w = eval(w_) - nw = len(w) - ss_ = "%(sym)s.ss[%(sym)s.n == %(n)s ]" % locals() - ss = eval(ss_) - - ssa[n,i,0] = n - ssa[n,i,1] = nw - - if len(ss) == 0: continue - - ssa[n,i,2] = np.average(ss) - ssa[n,i,3] = np.std(ss, ddof=1)/np.sqrt(np.size(ss)) - pass - pass - - - expr_ = "np.c_[ssa.reshape(-1,8)]" - print(expr_) - print(eval(expr_)) - - if PLOT == "PHO_AVG": - EXPRS_ = r""" - %(sym)s.f.baselabel # - %(sym)s.FAKES_SKIP # from run_meta - np.diff(%(sym)s.rr)[0]/1e6 # Run - %(sym)s.ee[-1]/1e6 # Evt - np.sum(%(sym)s.ss)/1e6 # Pho - np.sum(%(sym)s.ss)/%(sym)s.ee[-1] # Pho/Evt - """ - title = "sevt_tt.PHO_AVG.table NEVT:%(NEVT)d ## " % locals() - ranno = GetAnnotation(title, EXPRS_, syms ) - os.environ["RHSANNO"] = "\n".join(filter(lambda line:line.find("##")==-1, ranno.split("\n"))) - os.environ["RHSANNO_POS"]="0.45,-0.02" - - cut = 10 - labelx = "stat_cut:%(cut)d NEVT:%(NEVT)d " % locals() - label = "PHO_AVG : average beginPhoton->endPhoton duration [us] vs (step point count - 2) # " + labelx - print(label) - - fig, axs = mpplt_plotter(nrows=1, ncols=1, label=label, equal=False) - ax = axs[0] - - #sl = slice(2,16) - sl = slice(None) - style = {'a':"ro-", 'b':"bo-" } - fitstyle = {'a':"r:", 'b':"b:" } - - - mc = np.zeros( (len(syms), 2), dtype=np.float64 ) - - - for i,sym in enumerate(syms): - x_ = ssa[sl,i,0] - 2 # step point count - 2 (so starts from 0) - n_ = ssa[sl,i,1] - y_ = ssa[sl,i,2] - ye_ = ssa[sl,i,3] - - w = np.where( n_ > cut)[0] - - x = x_[w] - y = y_[w] - ye = ye_[w] - - label = "%(sym)s : avg(ss) vs n-2 " % locals() - - #ax.plot( x, y , style[sym], label=label ) - - ax.errorbar( x, y,yerr=ye, fmt=style[sym], label=label ) - - - - wf = np.where( n_ > 10)[0] ## only fit points with some stats - x = x_[wf] - y = y_[wf] - - A = np.vstack([x, np.ones(len(x))]).T - m, c = np.linalg.lstsq(A, y, rcond=None)[0] - mc[i,0] = m - mc[i,1] = c - ax.plot(x, m*x + c, fitstyle[sym], label='Fit: m %10.3f c %10.3f ' % (m,c) ) - pass - ax.legend() - fig.show() - - print(label) - expr = "np.c_[mc]" - print(expr) - print(eval(expr)) - - print("ranno") - print(ranno) - pass - - - if PLOT == "PHO_N": - - EXPRS_ = r""" - %(sym)s.f.baselabel # - %(sym)s.FAKES_SKIP # from run_meta - """ - title = "sevt_tt.PHO_N.table NEVT:%(NEVT)d ## " % locals() - ranno0 = GetAnnotation(title, EXPRS_, syms ) - ranno = "\n".join(filter(lambda line:line.find("##")==-1, ranno0.split("\n"))) - os.environ["RHSANNO"] = ranno - os.environ["RHSANNO_POS"]="0.1,-0.02" - - a_n, a_nc = np.unique(a.n, return_counts=True) - b_n, b_nc = np.unique(b.n, return_counts=True) - - #msg = "consistent when fakes are skipped" - msg = "much more in A when fakes not skipped : NEVT:%(NEVT)d" % locals() - label = "PHO_N : A(N=0),B(N=1) photon step counts : %s " % msg - fig, axs = mpplt_plotter(nrows=1, ncols=1, label=label, equal=False) - ax = axs[0] - ax.set_yscale("log") - #ax.plot( a_n, a_nc, "ro", label="A" ) - #ax.plot( b_n, b_nc, "bo", label="B" ) - - ax.errorbar( a_n, a_nc,yerr=np.sqrt(a_nc.astype(np.float64)),fmt="ro-", label="A" ) - ax.errorbar( b_n, b_nc,yerr=np.sqrt(b_nc.astype(np.float64)),fmt="bo-", label="B" ) - - ax.legend() - fig.show() - pass - - - - - if PLOT == "PHO_SCAT": - msg = "scatter plot of beginPhoton->endPhoton CPU time vs step point count" - label = "%s : %s" % (PLOT, msg) - - fig, axs = mpplt_plotter(nrows=1, ncols=1, label=label, equal=False) - ax = axs[0] - - ax.set_ylim(0,1000) - - for i,sym in enumerate(syms): - x_ = "%(sym)s.n.astype(np.float64)" % locals() - y_ = "%(sym)s.ss" % locals() - x = eval(x_) - y = eval(y_) - if i == 1: x += 0.2 - - ax.scatter( x[x<16], y[x<16], s=2 ) - pass - ax.legend() - fig.show() - pass - - - """ - STAMP NOTES - - HMM: how to annotate the stamp plot with photon indices ? - - For comparibility between A and B it makes more - sense to select based on time stamps rather than photon indices. - - To find time range with big bouncers:: - - a.s0[np.where(a.n > 20)]-a.ee[0] - - """ - if PLOT == "STAMP": - - - t0 = STAMP_TT[0] - t1 = STAMP_TT[0]+STAMP_TT[1] - stt = os.environ.get("STAMP_TT","-") - - label = "A:lhs, B:rhs : STAMP_TT=%s # (t0,t1):(%d,%d) " % (stt,t0,t1) - - #os.environ["SUBTITLE"] = subtitle - #os.environ["THIRDLINE"] = subtitle - - fig, axs = mpplt_plotter(nrows=1, ncols=1, label=label, equal=False) - ax = axs[0] - - #sl = slice(-100,None,1) # last 100 photons, which are actually processed first - #sl = slice(1040,1550,1) - sl = slice(None) - - s0_ = {} - s1_ = {} - tt_ = {} - wt_ = {} - h0_ = {} - h1_ = {} - i0_ = {} - i1_ = {} - - s0 = {} - s1 = {} - tt = {} - wt = {} - h0 = {} - h1 = {} - i0 = {} - i1 = {} - - sz = 0.01 - xx = {'a':[-0.5, -sz], 'b':[sz,0.5] } - zz = {'a':[-0.2, -sz], 'b':[sz,0.2] } - yy = {'a':-0.5, 'b':0.20 } - - qwns = "s0 s1 t h0 h1 i0 i1".split() - q_expr = "%(sym)s.%(q)s[sl] - %(sym)s.ee[0]" - - print(" (t0,t1) (%(t0)d,%(t1)d) " % locals() ) - - for i, sym in enumerate(syms): - print("sym:%s" % sym) - for j in range(len(qwns)): - q = qwns[j] - expr = q_expr % locals() - t = eval(expr) - w = eval("np.where(np.logical_and(t > t0, t < t1 ))[0]") - ws = str(w.shape) - qmn = t.min() - qmx = t.max() - fmt = " %(q)2s %(expr)22s : (%(qmn)8d,%(qmx)8d) : ws:%(ws)s " - print( fmt % locals() ) - pass - pass - - print(" when all the ws are 0: adjust the time range to find some stamps") - - - for i, sym in enumerate(syms): - r0 = eval("%(sym)s.rr[0]" % locals()) - #if i == 1:continue - print( "sym:%s r0:%10.3f " % (sym, r0)) - - ## where selecting on times greater than zero messes up the indices, - ## so instead just rely on exclusion of crazy times for unfilled cases - - cols = "r b g c m c m".split() - assert len(cols) == len(qwns) - - for j in range(len(qwns)): - q = qwns[j] - expr = q_expr % locals() - t = eval(expr) - w = eval("np.where(np.logical_and(t > t0, t < t1 ))[0]") - # photon indices of times within time window - kk = list(range(len(w))) - xmin,xmax = zz[sym] if q == "t" else xx[sym] - - ax.hlines( t[np.logical_and(t > t0, t < t1 )], xmin, xmax, cols[j], label=expr ) - - if "STAMP_ANNO" in os.environ: - for k in kk: - idx = w[k] - anno = "(%s) %s : %d : " % (q, sym.upper(), idx) - if q == "s0": - his = eval("%(sym)s.q[%(idx)s][0].decode(\"utf-8\").strip()" % locals()) - anno += his - elif q == "h0": - hc = eval("%(sym)s.hc[%(idx)s]" % locals()) - anno += " hc:%d " % (hc) - elif q == "i0": - ic = eval("%(sym)s.ic[%(idx)s]" % locals()) - hi0 = eval("%(sym)s.hi0[%(idx)s]" % locals()) - anno += "ic:%d hi0:%d "% (ic, hi0) - else: - anno = None - pass - if not anno is None: - ax.text( yy[sym], t[idx], anno ) - pass - pass - pass - pass - pass - - if "STAMP_LEGEND" in os.environ: - ax.legend() - pass - fig.show() - pass - if not ENVOUT is None and not CAP_STEM is None and not CAP_BASE is None: - envout = "\n".join([ - "export ENVOUT_PATH=%s" % ENVOUT, - "export ENVOUT_CAP_STEM=%s" % CAP_STEM, - "export ENVOUT_CAP_BASE=%s" % CAP_BASE, - "" - ]) - open(ENVOUT, "w").write(envout) - print(envout) - pass - - diff --git a/sysrap/sframe.py b/sysrap/sframe.py deleted file mode 100755 index 8c97d18899..0000000000 --- a/sysrap/sframe.py +++ /dev/null @@ -1,463 +0,0 @@ -#!/usr/bin/env python - -import logging -log = logging.getLogger(__name__) - -from opticks.ana.npmeta import NPMeta -from opticks.ana.axes import * -from opticks.ana.pvplt import mpplt_focus, SIZE - -import os, numpy as np -eary_ = lambda ekey, edef:np.array( list(map(float, os.environ.get(ekey,edef).split(","))) ) - -from opticks.ana.eget import efloat_ - - -class sframe(object): - @classmethod - def Load(cls, fold=None, name="sframe.npy"): - if fold is None: - fold = os.environ.get("FOLD", "") - pass - if fold.endswith(name): - path = fold - else: - path = os.path.join(fold, name) - pass - return cls(path) - - - @classmethod - def DetermineAxes(cls, nx, ny, nz): - """ - :param nx: - :param nx: - :param nx: - - With planar axes the order is arranged to make the longer axis the first horizontal one - followed by the shorter axis as the vertical one. - - +------+ - | | nz -> ( Y, Z ) ny_over_nz > 1 - +------+ - ny - - +------+ - | | ny -> ( Z, Y ) ny_over_nz < 1 - +------+ - nz - - """ - if nx == 0 and ny > 0 and nz > 0: - ny_over_nz = float(ny)/float(nz) - axes = (Y,Z) if ny_over_nz > 1 else (Z,Y) - elif nx > 0 and ny == 0 and nz > 0: - nx_over_nz = float(nx)/float(nz) - axes = (X,Z) if nx_over_nz > 1 else (Z,X) - elif nx > 0 and ny > 0 and nz == 0: - nx_over_ny = float(nx)/float(ny) - axes = (X,Y) if nx_over_ny > 1 else (Y,X) - else: - axes = (X,Y,Z) - pass - return axes - - - @classmethod - def CombineFrame(cls, framelist_): - """ - HMM: probably should do more consistency checks between all the frames - """ - framelist = list(filter(None, framelist_)) - axes = None - for f in framelist: - if axes is None: - axes = f.axes - else: - assert axes == f.axes - pass - pass - return framelist[0] if len(framelist) > 0 else None - - - def __init__(self, path, clear_identity=True ): - """ - Whether clear_identity makes any material difference depends on the identity values. - But it should be done anyhow. For some identity values they will appear as nan in float. - """ - if path is None: - return - pass - - metapath = path.replace(".npy", "_meta.txt") - if os.path.exists(metapath): - meta = NPMeta.Load(metapath) - else: - meta = None - pass - self.meta = meta - - a = np.load(path) - i = a.view(np.int32) - - self.shape = a.shape # for fold to treat like np.array - self.path = path - self.a = a - self.i = i - - ce = a[0,0] - ix0,ix1,iy0,iy1 = i[0,1,:4] # q1.i.xyzw - iz0,iz1,num_photon = i[0,2,:3] # q2.i.xyz - gridscale = a[0,2,3] # q2.f.w - - - midx, mord, gord = i[0,3,:3] # q3.i.xyz - inst = i[0,3,3] # q3.i.w - - - - self.midx = midx - self.mord = mord - self.gord = gord - self.inst = inst - - ## replacing ana/gridspec.py:GridSpec - - - - grid = "".join(["ix0 %(ix0)4d ix1 %(ix1)4d ", - "iy0 %(iy0)4d iy1 %(iy1)4d ", - "iz0 %(iz0)4d iz1 %(iz1)4d ", - "num_photon %(num_photon)4d gridscale %(gridscale)10.4f"]) % locals() - - - - - - - self.grid = grid - self.ce = ce - self.ix0 = ix0 - self.ix1 = ix1 - self.iy0 = iy0 - self.iy1 = iy1 - self.iz0 = iz0 - self.iz1 = iz1 - self.num_photon = num_photon - self.gridscale = gridscale - - - - target = "midx %(midx)6d mord %(mord)6d gord %(gord)6d inst %(inst)7d " % locals() - self.target = target - - - ins_idx_1 = i[1,0,3] - 1 # q0.u.w-1 see sqat4.h qat4::getIdentity qat4::setIdentity - gas_idx_1 = i[1,1,3] - 1 # q1.u.w-1 - ias_idx_1 = i[1,2,3] - 1 # q2.u.w-1 - - ins_idx_2 = i[2,0,3] - 1 # q0.u.w-1 see sqat4.h qat4::getIdentity qat4::setIdentity - gas_idx_2 = i[2,1,3] - 1 # q1.u.w-1 - ias_idx_2 = i[2,2,3] - 1 # q2.u.w-1 - - ins_idx_match = ins_idx_1 == ins_idx_2 - gas_idx_match = gas_idx_1 == gas_idx_2 - ias_idx_match = ias_idx_1 == ias_idx_2 - - if not ins_idx_match: log.warning(f" ins_idx_match {ins_idx_match} ins_idx_1 {ins_idx_1} ins_idx_2 {ins_idx_2} ") - if not gas_idx_match: log.warning(f" gas_idx_match {gas_idx_match} gas_idx_1 {gas_idx_1} gas_idx_2 {gas_idx_2} ") - if not ias_idx_match: log.warning(f" ias_idx_match {ias_idx_match} ias_idx_1 {ias_idx_1} ias_idx_2 {ias_idx_2} ") - - #assert ins_idx_match - #assert gas_idx_match - #assert ias_idx_match - - ins_idx = ins_idx_1 - gas_idx = gas_idx_1 - ias_idx = ias_idx_1 - qat4id = "ins_idx %(ins_idx)6d gas_idx %(gas_idx)4d %(ias_idx)4d " % locals() - - self.ins_idx = ins_idx - self.gas_idx = gas_idx - self.ias_idx = ias_idx - self.qat4id = qat4id - - m2w = a[1].copy() - w2m = a[2].copy() - - if clear_identity: - m2w[0,3] = 0. - m2w[1,3] = 0. - m2w[2,3] = 0. - m2w[3,3] = 1. - - w2m[0,3] = 0. - w2m[1,3] = 0. - w2m[2,3] = 0. - w2m[3,3] = 1. - pass - - self.m2w = m2w - self.w2m = w2m - self.id = np.dot( m2w, w2m ) - - - ins = i[3,0,0] # aux.q0.i.x - gas = i[3,0,1] - ias = i[3,0,2] - - propagate_epsilon = a[3,1,0] # aux.q1.f.x - - - ins_gas_ias = " ins %(ins)6d gas %(gas)4d ias %(ias)4d " % locals() - - self.ins = ins - self.gas = gas - self.ias = ias - self.ins_gas_ias = ins_gas_ias - self.propagate_epsilon = propagate_epsilon - - self.init_grid() - self.init_view() - - - def init_grid(self): - """ - replacing ana/gridspec.py - """ - gord = self.gord - ix0 = self.ix0 - ix1 = self.ix1 - iy0 = self.iy0 - iy1 = self.iy1 - iz0 = self.iz0 - iz1 = self.iz1 - gridscale = self.gridscale - ce = self.ce - extent = ce[3] - s_extent = "e:%7.3f" % extent - - nx = (ix1 - ix0)//2 - ny = (iy1 - iy0)//2 - nz = (iz1 - iz0)//2 - coords = "RTP" if gord == -3 else "XYZ" ## NB RTP IS CORRECT ORDERING radiusUnitVec:thetaUnitVec:phiUnitVec - axes = self.DetermineAxes(nx, ny, nz) - other_axis = Axes.OtherAxis(axes) - planar = len(axes) == 2 - - bbox = np.zeros( (2,3), dtype=np.float32 ) - bbox[0] = list(map(float,[ix0,iy0,iz0])) - bbox[1] = list(map(float,[ix1,iy1,iz1])) - bbox *= gridscale*ce[3] - - - if planar: - H, V = axes - axlabels = coords[H], coords[V] - s_bbox = "bb " - s_bbox += " %7.4f %7.4f " % tuple(bbox[0,axes]) - s_bbox += " %7.4f %7.4f " % tuple(bbox[1,axes]) - else: - H, V, D = axes - axlabels = coords[H], coords[V], coords[D] - s_bbox = "bb %7.4f %7.4f %7.4f %7.4f %7.4f %7.4f " % tuple(bbox.ravel()) - pass - - - self.nx = nx - self.ny = ny - self.nz = nz - self.coords = coords - self.axes = axes - self.other_axis = other_axis - self.axlabels = axlabels - - self.bbox = bbox # without ce_offset - self.s_bbox = s_bbox - - self.extent = extent - self.s_extent = s_extent - - - - def init_view(self): - """ - HMM: input EYE is ignored for planar - """ - ce = self.ce - axes = self.axes - coords = self.coords - planar = len(axes) == 2 - - # below default from envvars are overridden for planar data - eye = eary_("EYE","1.,1.,1.") - up = eary_("UP","0.,0.,1.") - off = eary_("OFF","0.,0.,1.") - look = eary_("LOOK","0.,0.,0.") - #look = ce[:3] - - EYES = efloat_("EYES", "6.") # TODO: why 6 ? how to control FOV to gain more control of this - extent = ce[3] - - if planar: - H, V = axes - up = Axes.Up(H,V) - off = Axes.Off(H,V) - eye = look + extent*off*EYES - else: - H, V, D = axes - axlabels = coords[H], coords[V], coords[D] - - up = XYZ.up - off = XYZ.off - ## hmm in 3D case makes less sense : better to just use the input EYE - - eye = ce[3]*eye*EYES - pass - - self.look = look - self.up = up - self.off = off - self.eye = eye - self.thirdline = self.s_bbox + self.s_extent - - - def pv_compose(self, pl, local=True): - """ - #pl.view_xz() ## TODO: see if view_xz is doing anything when subsequently set_focus/viewup/position - """ - - RESET = "RESET" in os.environ - PARA = "PARA" in os.environ - ZOOM = efloat_("ZOOM", "1.") - PVGRID = "PVGRID" in os.environ - - if PARA: - pl.camera.ParallelProjectionOn() - pass - - look = self.look if local else self.ce[:3] ## HMM: same ? - eye = look + self.off - up = self.up - - ## for reset=True to succeed to auto-set the view, must do this after add_points etc.. - - #eye = look + self.off - - look = self.look - eye = self.eye - up = self.up - - log.info("frame.pv_compose look:%s eye: %s up:%s PARA:%s RESET:%d ZOOM:%s " % (str(look), str(eye), str(up), RESET, PARA, ZOOM )) - - pl.set_focus( look ) - pl.set_viewup( up ) - pl.set_position( eye, reset=RESET ) ## for reset=True to succeed to auto-set the view, must do this after add_points etc.. - pl.camera.Zoom(ZOOM) - - if PVGRID: - pl.show_grid() - pass - - - - def __repr__(self): - - l_ = lambda k,v:"%-12s : %s" % (k, v) - - return "\n".join( - [ - l_("sframe",""), - l_("path",self.path), - l_("meta",repr(self.meta)), - l_("ce", repr(self.ce)), - l_("grid", self.grid), - l_("bbox", repr(self.bbox)), - l_("target", self.target), - l_("qat4id", self.qat4id), - l_("m2w",""), repr(self.m2w), "", - l_("w2m",""), repr(self.w2m), "", - l_("id",""), repr(self.id) , - l_("ins_gas_ias",self.ins_gas_ias) - ]) - - - @classmethod - def FakeXZ(cls, e=10): - path = None - fr = cls(path) - fr.axes = X,Z - fr.axlabels = "X","Z" - fr.bbox = np.array( [[-e,-e,-e],[e,e,e]] ) - return fr - - - def mp_subplots(self, mp): - pass - H,V = self.axes - _H,_V = self.axlabels - - xlim = self.bbox[:,H] - ylim = self.bbox[:,V] - - log.info("mp_subplots xlim %s ylim %s " % (str(xlim), str(ylim))) - #xlim, ylim = mpplt_focus(xlim, ylim) - - topline = os.environ.get("TOPLINE", "sframe.py:mp_subplots:TOPLINE") - botline = os.environ.get("BOTLINE", "sframe.py:mp_subplots:BOTLINE") - thirdline = getattr(self, "thirdline", "thirdline" ) - - title = [topline, botline, thirdline] - - fig, ax = mp.subplots(figsize=SIZE/100.) - fig.suptitle("\n".join(title)) - - ax.set_aspect('equal') - ax.set_xlim(xlim) - ax.set_ylim(ylim) - ax.set_xlabel(_H) - ax.set_ylabel(_V) - - self.fig = fig - self.ax = ax - - return fig, ax - - def mp_scatter(self, pos, **kwa): - assert pos.ndim == 2 and pos.shape[1] == 3 - H,V = self.axes - self.ax.scatter( pos[:,H], pos[:,V], **kwa) - - def mp_arrow(self, pos, mom, **kwa): - assert pos.ndim == 2 and pos.shape[1] == 3 - assert mom.ndim == 2 and mom.shape[1] == 3 - assert len(pos) == len(mom) - H,V = self.axes - - for i in range(len(pos)): - self.ax.arrow( pos[i,H], pos[i,V], mom[i,H], mom[i,V] ) - pass - - - def mp_legend(self): - NOLEGEND = "NOLEGEND" in os.environ - LSIZ = float(os.environ.get("LSIZ",50)) - LOC = os.environ.get("LOC", "upper right") - log.info("mp_legend LSIZ %s LOC %s NOLEGEND %s " % (LSIZ,LOC,NOLEGEND) ) - if not NOLEGEND: - lgnd = self.ax.legend(loc=LOC) - for handle in lgnd.legendHandles: - handle.set_sizes([LSIZ]) - pass - pass - - - - -if __name__ == '__main__': - - fr = sframe.Load() - print(fr) - - - diff --git a/sysrap/sfreq.py b/sysrap/sfreq.py deleted file mode 100644 index 85d156f6da..0000000000 --- a/sysrap/sfreq.py +++ /dev/null @@ -1,75 +0,0 @@ -#!/usr/bin/env python - -import os, numpy as np -from opticks.ana.fold import Fold - -class fold(object): - pass - -class sfreq(object): - @classmethod - def CreateFromArray(cls, a): - usub, nsub = np.unique(a, return_counts=True ) - f = fold() - f.key = usub - f.val = nsub - return cls(f, sort=True) - - @classmethod - def ArrayFromFile(cls, base=None, name="subs.txt"): - if base is None: - base = os.environ.get("FOLD", None) - pass - path = os.path.join(base, name) - a = np.loadtxt(path, dtype="|S32" ) - return a - - @classmethod - def CreateFromFile(cls, base=None, name="subs.txt"): - a = cls.ArrayFromFile(base=base, name=name) - return cls.CreateFromArray(a) - - def __init__(self, f, sort=True): - order = np.argsort(f.val)[::-1] if sort else slice(None) - okey = f.key[order] - oval = f.val[order] - subs = list(map(lambda _:_.decode("utf-8"), okey.view("|S32").ravel() )) - vals = list(map(int, oval)) - - self.f = f # without the sort - self.okey = okey - self.oval = oval - self.subs = subs - self.vals = vals - self.order = order - - - def find_index(self, key): - key = key.encode() if type(key) is str else key - assert type(key) is bytes - ii = np.where( self.okey.view("|S32") == key )[0] ## hmm must use okey to feel the sort - assert len(ii) == 1 - return int(ii[0]) - - def desc_key(self, key): - idx = self.find_index(key) - return self.desc_idx(idx) - - def desc_idx(self, idx): - return "sf %3d : %7d : %s." % (idx, self.vals[idx], self.subs[idx]) - - def __repr__(self): - return "\n".join(self.desc_idx(idx) for idx in range(len(self.subs))) - - def __str__(self): - return str(self.f) - - -if __name__ == '__main__': - sf = sfreq.CreateFromFile() - print(repr(sf)) - - - - - diff --git a/sysrap/smath.py b/sysrap/smath.py deleted file mode 100644 index e66a961e68..0000000000 --- a/sysrap/smath.py +++ /dev/null @@ -1,111 +0,0 @@ -#!/usr/bin/env python - -import numpy as np - -def rotateUz(d, u): - """ - :param d: array of vectors with shape (n,3) - :param u: array of single vector with shape (3,) - :return d1: array of vectors with shape (n,3) - - NumPy translation of smath.h rotateUz:: - - inline SMATH_METHOD void smath::rotateUz(float3& d, const float3& u ) - { - float up = u.x*u.x + u.y*u.y ; - if (up>0.f) - { - up = sqrt(up); - float px = d.x ; - float py = d.y ; - float pz = d.z ; - d.x = (u.x*u.z*px - u.y*py)/up + u.x*pz; - d.y = (u.y*u.z*px + u.x*py)/up + u.y*pz; - d.z = -up*px + u.z*pz; - } - else if (u.z < 0.f ) - { - d.x = -d.x; - d.z = -d.z; - } - } - - """ - assert u.shape == (3,) - assert d.shape == (len(d),3) - d1 = np.zeros( (len(d),3 ), dtype=d.dtype ) - - up = u[0]*u[0] + u[1]*u[1] - if up > 0.: - up = np.sqrt(up) - px = d[:,0] - py = d[:,1] - pz = d[:,2] - d1[:,0] = (u[0]*u[2]*px - u[1]*py)/up + u[0]*pz - d1[:,1] = (u[1]*u[2]*px + u[0]*py)/up + u[1]*pz - d1[:,2] = -up*px + u[2]*pz - elif u[2] < 0.: - d1[:,0] = -d[:,0] - d1[:,1] = d[:,1] # diff as returning not inplace changing - d1[:,2] = -d[:,2] - pass - return d1 - - -def rotateUz_(d, u): - """ - :param d: array of single vector with shape (3,) - :param u: array of vectors with shape (n,3) - :return d1: array of vectors with shape (n,3) - - THIS NEEDS MORE TESTING - """ - assert u.shape == (len(u),3) - assert d.shape == (3,) - d1 = np.zeros( (len(u),3 ), dtype=d.dtype ) - - _up = u[:,0]*u[:,0] + u[:,1]*u[:,1] - wp = np.where( _up > 0. ) - wn = np.where( np.logical_and( _up <= 0., u[:,2] < 0. )) - - up = np.zeros( (len(u),), dtype=d.dtype ) - up[wp] = np.sqrt(_up) - px = d[0] - py = d[1] - pz = d[2] - - d1[wp,0] = (u[wp,0]*u[wp,2]*px - u[wp,1]*py)/up + u[wp,0]*pz - d1[wp,1] = (u[wp,1]*u[wp,2]*px + u[wp,0]*py)/up + u[wp,1]*pz - d1[wp,2] = -up[wp]*px + u[wp,2]*pz - - d1[wn,0] = -d[0] - d1[wn,1] = d[1] - d1[wn,2] = -d[2] - - return d1 - - -def test_rotateUz(): - s = np.sqrt(0.5, dtype=np.float64) - u = np.array([s,0,-s], dtype=s.dtype) - d = np.array([[1,0,0],[s,s,0],[0,1,0],[-s,s,0],[-1,0,0],[-s,-s,0],[0,-1,0],[s,-s,0],[1,0,0]], dtype=s.dtype ) - d1 = rotateUz( d, u ) - - print("u\n",u) - print("d\n",d) - print("d1\n",d1) - - - -if __name__ == '__main__': - - s = np.sqrt(0.5, dtype=np.float64) - u = np.array([[1,0,0],[s,s,0],[0,1,0],[-s,s,0],[-1,0,0],[-s,-s,0],[0,-1,0],[s,-s,0],[1,0,0]], dtype=s.dtype ) - d = np.array([s,0,-s], dtype=s.dtype) - d1 = rotateUz_( d, u ) - - print("u\n",u) - print("d\n",d) - print("d1\n",d1) - - diff --git a/sysrap/smonitor.py b/sysrap/smonitor.py deleted file mode 100644 index 7a2c1db68f..0000000000 --- a/sysrap/smonitor.py +++ /dev/null @@ -1,108 +0,0 @@ -#!/usr/bin/env python -""" -smonitor.py -============ - -~/o/sysrap/smonitor.sh ana - -""" - -import os, numpy as np -import matplotlib.pyplot as mp -SIZE=np.array([1280, 720]) - -if __name__ == '__main__': - mon = np.load("smonitor.npy").astype(np.int64) - - stamp = mon[:,0] - device = mon[:,1] - free = mon[:,2] - total = mon[:,3] - used = mon[:,4] - pid = mon[:,5] - usedGpuMemory = mon[:,6] - proc_count = mon[:,7] - - t = (stamp - stamp[0])/1e6 - usedGpuMemory_GB = usedGpuMemory/1e9 - total_GB = total/1e9 - free_GB = free/1e9 - used_GB = used/1e9 - - u_device = np.unique(device) - u_total = np.unique(total) - u_pid = np.unique(pid) - u_proc_count = np.unique(proc_count) - - assert len(u_device) == 1 - assert len(u_total) == 1 - assert len(u_pid) == 1 - - _device = u_device[0] - _total = u_total[0]/1e9 - _pid = u_pid[0] - _proc_count = u_proc_count[0] - - assert _proc_count == 1 - - delta_usedGpuMemory_GB = np.diff( usedGpuMemory_GB ) - w_delta_usedGpuMemory_GB = np.where( delta_usedGpuMemory_GB > 0.001 )[0] - - auto_start = w_delta_usedGpuMemory_GB[2] if len(w_delta_usedGpuMemory_GB) > 2 else 0 - env_start = int(os.environ["START"]) if "START" in os.environ else 0 - start = env_start if env_start > 0 else auto_start - message = "auto_start:%(auto_start)2d START:%(env_start)2d start:%(start)2d " % locals() - - - idx = np.arange( len(usedGpuMemory_GB) ) - sel = slice(start, None) - - expr = "np.c_[idx[sel], t[sel], usedGpuMemory_GB[sel], usedGpuMemory[sel] ]" - print(expr) - print(eval(expr)) - - _dgb = "(usedGpuMemory_GB[sel][-1]-usedGpuMemory_GB[sel][0])" - dgb = eval(_dgb) - print("dgb %10.3f %s " % ( dgb, _dgb )) - - _db = "(usedGpuMemory[sel][-1]-usedGpuMemory[sel][0])" - db = eval(_db) - print("db %10.3f %s " % ( db, _db )) - - - - _dt = "(t[sel][-1]-t[sel][0])" - dt = eval(_dt) - print("dt %10.3f %s " % ( dt, _dt )) - print("dgb/dt %10.3f " % (dgb/dt)) - print("db/dt %10.3f " % (db/dt)) - - - deg = 1 # linear - pfit = np.polyfit(t[sel], usedGpuMemory_GB[sel], deg) - linefit = np.poly1d(pfit) - linefit_label = "line fit: slope %10.3f [GB/s] intercept %10.3f " % (linefit.coef[0], linefit.coef[1]) - - headline = "smonitor.sh device %(_device)s total_GB %(_total)4.1f pid %(_pid)s " % locals() - publine = "%s : %s " % (message, os.environ.get("PUB", "no-PUB" )) - title = "\n".join([headline, linefit_label, publine]) - print(title) - - fig, ax = mp.subplots(figsize=SIZE/100.) - fig.suptitle(title) - - ax.set_yscale('log') - - ax.plot( t, total_GB , label="total_GB" ) - ax.plot( t, free_GB , label="free_GB" ) - ax.plot( t, used_GB , label="used_GB" ) - - #ax.plot( t, usedGpuMemory_GB , label="proc.usedGpuMemory_GB" ) - ax.scatter( t, usedGpuMemory_GB, label="proc.usedGpuMemory_GB" ) - ax.plot( t[sel], linefit(t[sel]), label=linefit_label ) - - - ax.legend() - fig.show() - - diff --git a/sysrap/sn_check.py b/sysrap/sn_check.py deleted file mode 100644 index ed4342a502..0000000000 --- a/sysrap/sn_check.py +++ /dev/null @@ -1,121 +0,0 @@ -#!/usr/bin/env python -""" -sn_check.py -============ - - -""" - -import numpy as np - -class snd_check(object): - def __init__(self, f, symbol="c"): - lvn = f.soname_names - snd = f.csg.node - snd_fields = snd.shape[1] - assert snd_fields == 17 - self.snd = snd - - -class sn_check(object): - def __init__(self, f, symbol="c"): - lvn = f.soname_names - sn = f._csg.sn - sn_fields = sn.shape[1] - assert sn_fields in [17,19] - with_child = sn_fields == 19 - - s_pa = f._csg.s_pa - s_tv = f._csg.s_tv - s_bb = f._csg.s_bb - - self.sn = sn - self.sn_fields = sn_fields - self.with_child = with_child - - self.s_pa = s_pa - self.s_tv = s_tv - self.s_bb = s_bb - - - lv = sn[:,2] - tv = sn[:,3] - pa = sn[:,4] - bb = sn[:,5] - parent = sn[:,6] - - if with_child: - pass - else: - left = sn[:,7] - right = sn[:,8] - - self.left = left - self.right = right - self.ucheck("left") - self.ucheck("right") - self.check_left() - self.check_right() - pass - - self.lvn = lvn - self.symbol = symbol - - self.lv = lv - self.tv = tv - self.pa = pa - self.bb = bb - self.parent = parent - - for field in "lv tv pa bb parent".split(): - self.ucheck(field) - pass - self.check_ulv() - self.check_utv() - self.check_upa() - self.check_ubb() - self.check_parent() - - def ucheck(self, field="lv"): - uchk = np.unique(getattr(self,field), return_counts=True) - setattr(self, "u%(field)s" % locals(), uchk[0] ) - setattr(self, "n%(field)s" % locals(), uchk[1] ) - - def check_ulv(self): - assert np.all( self.ulv == np.arange(len(self.ulv)) ) - def check_utv(self): - assert self.utv.max() < len(self.s_tv), "integrity of transform refs" - def check_upa(self): - assert self.upa.max() < len(self.s_pa), "integrity of param refs" - def check_ubb(self): - assert self.ubb.max() < len(self.s_bb), "integrity of bbox refs" - def check_parent(self): - assert self.parent.max() < len(self.sn), "integrity of parent node refs" - # assert np.all( self.nparent[1:] == 1 ) - # not fulfilled by the virtual mask polycone - - def check_left(self): - assert self.left.max() < len(self.sn), "integrity of left node refs" - assert np.all( self.nleft[1:] == 1 ) - def check_right(self): - assert self.right.max() < len(self.sn), "integrity of right node refs" - assert np.all( self.nright[1:] == 1 ) - - - def __repr__(self): - sym = self.symbol - #locals()[sym] = self ## curious locals worked with older python, now need globals - globals()[sym] = self - expr = "np.c_[%(sym)s.ulv, %(sym)s.nlv, %(sym)s.lvn[%(sym)s.ulv]][np.where(%(sym)s.nlv>5)]" % locals() - lines = [] - lines.append("opticks.sysrap.sn_check") - lines.append(expr) - lines.append(str(eval(expr))) - lines.append("# NB raw nodes from G4VSolid, not complete binary tree counts") - return "\n".join(lines) - -if __name__ == '__main__': - print("klop") -pass - - diff --git a/sysrap/sphotonlite.py b/sysrap/sphotonlite.py deleted file mode 100644 index 149d49f0f4..0000000000 --- a/sysrap/sphotonlite.py +++ /dev/null @@ -1,258 +0,0 @@ -""" -opticks/sysrap/sphotonlite.py -============================== - -:: - - f = Fold.Load(symbol="f") - from opticks.sysrap.sphotonlite import SPhotonLite - p = SPhotonLite.view(f.photonlite) # photonlite shape (N,4) uint32 - -""" -from __future__ import annotations -import logging -log = logging.getLogger(__name__) - -import numpy as np -from opticks.ana.fold import append_fields - -EFFICIENCY_COLLECT = 0x1 << 13 - - -class SPhotonLite: - """ - Structured view of the 16-byte ``sphotonlite`` record. - - The binary layout (little-endian, 16 B total) is: - - 0-1 : identity (u2) - 2-3 : hitcount (u2) - 4-7 : time (f4) - 8-9 : lposfphi (u2) → normalised to lposfphi_f (f4) - 10-11 : lposcost (u2) → normalised to lposcost_f (f4) - 12-15 : flagmask (u4) - - All public API is provided as **class-methods** – no instance needed. - """ - - # ------------------------------------------------------------------ - # 2. The dtype – defined once, reused everywhere - # ------------------------------------------------------------------ - DTYPE = np.dtype( - [ - ("identity", " np.ndarray: - """0xffff → 1.0, 0 → 0.0 (exact, no rounding error)""" - return u16.astype(np.float32) / 0xFFFF - - # ------------------------------------------------------------------ - # 4. Core conversion: uint32[N,4] → structured array - # ------------------------------------------------------------------ - @classmethod - def view_(cls, buf: np.ndarray) -> np.ndarray: - """ - Convert a ``uint32`` buffer of shape ``(..., 4)`` into a flat - structured array with the extra float fields ``lposfphi_f`` - and ``lposcost_f``. - - Parameters - ---------- - buf : np.ndarray - Must be ``dtype=uint32`` and ``buf.shape[-1] == 4``. - - Returns - ------- - np.ndarray - Flat (1-D) structured array. - """ - if buf.dtype != np.uint32: - raise TypeError(f"buf must be uint32, got {buf.dtype}") - if buf.shape[-1] != 4: - raise ValueError(f"last dimension must be 4, got {buf.shape[-1]}") - - rec = buf.view(cls.DTYPE).reshape(-1) - - - fphi = cls.unpack_uint16_to_float(rec["lposfphi"]) - cost = cls.unpack_uint16_to_float(rec["lposcost"]) - - phi = ( fphi * 2.0 - 1. ) * np.pi # scuda.h:phi_from_fphi - cos = np.arccos(cost) - - rec = append_fields( - rec - , - [ - "lposfphi_f", - "lposcost_f", - "phi", - "cost", - "cos" - ] - , - [ - fphi, - cost, - phi, - cost, - cos - ] - , - [ - np.float32, - np.float32, - np.float32, - np.float32, - np.float32, - ] - , - usemask=False, - ) - return rec - - @classmethod - def view(cls, buf: np.ndarray) -> np.recarray: - rec = cls.view_(buf) - return rec.view(np.recarray) - - - # ------------------------------------------------------------------ - # 6. Pretty-print a few records (useful in notebooks / REPL) - # ------------------------------------------------------------------ - - @classmethod - def pprint(cls, rec: np.ndarray, n: int = 5) -> None: - """Print the first *n* records in a readable table.""" - print(rec[:n]) - - # ------------------------------------------------------------------ - # 7. Optional: create an *empty* array of a given length - # ------------------------------------------------------------------ - @classmethod - def empty(cls, size: int) -> np.ndarray: - """Return a zero-filled structured array of length *size*.""" - return np.zeros(size, dtype=cls.DTYPE) - - - - - - -class SPhotonLite_Merge: - """ - Select hits and merge them using a NumPy implementation - that follows the thrust::reduce_by_key based approach - of sysrap/SPM.cu SPM::merge_partial_select - - For a small number of photons (eg M1) - this precisely duplicates in python the hit selection and merging - done with the CUDA impl - - """ - - def __init__(self, tw=1.0, select_mask = EFFICIENCY_COLLECT ): - self.tw = tw - self.select_mask = select_mask - - def __call__(self, _pl): - """ - :param _pl: raw photonlite input array - :return hlm: raw hitlitemerged array - - Q: If this is applied twice or more, does the hlm array stay the same ? - A: It should do. YES IT DOES STAY THE SAME - """ - log.info("[ photonlite.shape %s " % (repr(_pl.shape),) ) - - tw = self.tw - select_mask = self.select_mask - - # Step 1: Filter hits - flagmask = _pl[:,3] - - select = (flagmask & select_mask) != 0 - if not np.any(select): return np.zeros(0, dtype=_pl.dtype ) - - _hl = _pl[select] - # _hl = _pl[ np.where( _pl[:,3] & (0x1 << 13) )] - - # Step 2: Build 64-bit key: (pmt_id << 48) | bucket - identity = _hl[:,0] & 0xFFFF # lower 16 bits = PMT ID - time = _hl[:,1].view(np.float32) - bucket = np.floor(time / tw).astype(np.uint64) - key = (identity.astype(np.uint64) << 48) | bucket - - # key = ( ( _hl[:,0] & 0xFFFF ).astype(np.uint64) << 48 ) | np.floor( _hl[:,1].view(np.float32)/1. ).astype(np.uint64) - - # Step 3: Sort by key (exactly what Thrust does internally) - - log.info("-argsort[") - sort_idx = np.argsort(key, kind='stable') # stable = same as thrust::stable_sort_by_key - log.info("-argsort]") - key_sorted = key[sort_idx] - _hl_sorted = _hl[sort_idx] - - - # Step 4: Find group boundaries - log.info("-diff[") - key_diff = np.diff(key_sorted, prepend=key_sorted[0]-1) # force first group start - log.info("-diff]") - group_start = np.where(key_diff != 0)[0] - - # Number of output groups - n_groups = len(group_start) - - # Step 5: Reduce each group (vectorized version of sphotonlite_reduce_op) - hlm = np.zeros( (n_groups,4), dtype=_pl.dtype ) - - # Take first hit in group as base (like your CUDA reduce does with 'a') - first_in_group = group_start - - hlm[:] = _hl_sorted[first_in_group] - - # Now reduce the rest of each group - - log.info("-reducing n_groups %d [" % (n_groups,) ) - for i in range(n_groups): - start = group_start[i] - end = group_start[i+1] if i+1 < n_groups else len(key) - - if end - start == 1: - continue - - group_slice = slice(start, end) - - # min time - hlm[i,1] = np.min(_hl_sorted[group_slice,1].view(np.float32)).view(np.uint32) - - # OR all flagmasks - hlm[i,3] |= np.bitwise_or.reduce(_hl_sorted[group_slice,3]) - - all_identity = _hl_sorted[group_slice,0] & 0xffff - - # sum hitcounts - sum_hitcount = np.sum(_hl_sorted[group_slice,0] >> 16, dtype=np.uint32) - hlm[i,0] = sum_hitcount << 16 | all_identity[0] - pass - log.info("-reducing n_groups %d ]" % (n_groups,) ) - log.info("] hlm.shape %s " % ( repr(hlm.shape), )) - return hlm - pass -pass - - - diff --git a/sysrap/stag.py b/sysrap/stag.py deleted file mode 100755 index 185c949261..0000000000 --- a/sysrap/stag.py +++ /dev/null @@ -1,275 +0,0 @@ -#!/usr/bin/env python - -import numpy as np -import os, re, logging -log = logging.getLogger(__name__) -from collections import OrderedDict as odict - - -class stag_item(object): - @classmethod - def Placeholder(cls): - return cls(-1,"placeholder","ERROR" ) - - def __init__(self, code, name, note=""): - self.code = code - self.name = name - self.note = note - - def __str__(self): - return "%2d : %10s : %s " % (self.code, self.name, self.note) - def __repr__(self): - return "%2d : %10s" % (self.code, self.name) - - -class stag(object): - """ - # the below NSEQ, BITS, ... param need to correspond to stag.h static constexpr - """ - enum_ptn = re.compile("^\s*(\w+)\s*=\s*(.*?),*\s*?$") - note_ptn = re.compile("^\s*static constexpr const char\* (\w+)_note = \"(.*)\" ;\s*$") - - PATH = "$OPTICKS_PREFIX/include/sysrap/stag.h" - - NSEQ = 4 ## must match stag.h:NSEQ - BITS = 4 ## must match stag.h:BITS - MASK = ( 0x1 << BITS ) - 1 - SLOTMAX = 64//BITS - SLOTS = SLOTMAX*NSEQ - - - @classmethod - def NumStarts(cls, tg): - ns = np.zeros( (len(tg)), dtype=np.uint8 ) - for i in range(len(tg)): - starts = np.where( tg[i] == tg[0,0] )[0] - ns[i] = len(starts) - pass - return ns - - - @classmethod - def StepSplit(cls, tg, fl=None): - """ - Hmm maybe StepFold is clearer name - - :param tg: unpacked tag array of shape (n, SLOTS) - :param fl: None or flat array of shape (n, SLOTS) - :return tgs OR (tgs,fls): step split arrays of shape (n, max_starts, max_slots) - - In [4]: at[0] - Out[4]: array([ 1, 2, 9, 10, 1, 2, 9, 10, 1, 2, 11, 12, 0, 0, 0, 0], dtype=uint8) - - In [8]: ats[0] - Out[8]: - array([[ 1, 2, 9, 10, 0, 0, 0, 0, 0, 0], - [ 1, 2, 9, 10, 0, 0, 0, 0, 0, 0], - [ 1, 2, 11, 12, 0, 0, 0, 0, 0, 0], - [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]], dtype=uint8) - - """ - if not fl is None: - assert fl.shape == tg.shape - pass - - max_starts = 0 ## corresponds to the maximum number of steps of all photons - max_slots = 0 - for i in range(len(tg)): - starts = np.where( tg[i] == tg[0,0] )[0] # indices of where the very first tag appears in this photons tags - if len(starts) > max_starts: max_starts = len(starts) - ends = np.where( tg[i] == 0 )[0] - end = ends[0] if len(ends) > 0 else len(tg[i]) - mkr = np.concatenate( (starts, np.array((end,))) ) - mkr_slots = np.diff(mkr).max() - if mkr_slots > max_slots: max_slots = mkr_slots - pass - print("max_starts:%s max_slots:%d" % (max_starts, max_slots)) - - tgs = np.zeros((len(tg), max_starts, max_slots), dtype=np.uint8) - fls = np.zeros((len(tg), max_starts, max_slots), dtype=np.float32) if not fl is None else None - - for i in range(len(tg)): - starts = np.where( tg[i] == tg[0,0] )[0] - ends = np.where( tg[i] == 0 )[0] - end = ends[0] if len(ends) > 0 else len(tg[i]) - ## above handles when the tags do not get to zero due to collection truncation - - for j in range(len(starts)): - st = starts[j] - en = starts[j+1] if j+1 < len(starts) else end - tgs[i, j,0:en-st] = tg[i,st:en] - if not fls is None: - fls[i, j,0:en-st] = fl[i,st:en] - pass - pass - pass - return tgs if fls is None else tgs,fls - - - - @classmethod - def Unpack(cls, tag): - """ - :param tag: (n, NSEQ) array of bitpacked tag enumerations - :return tg: (n, SLOTS) array of unpacked tag enumerations - - Usage:: - - # apply stag.Unpack to both as same stag.h bitpacking is used - at = stag.Unpack(a.tag) if hasattr(a,"tag") else None - bt = stag.Unpack(b.tag) if hasattr(b,"tag") else None - - """ - assert tag.shape == (len(tag), cls.NSEQ) - - st = np.zeros( (len(tag), cls.SLOTS), dtype=np.uint8 ) - for i in range(cls.NSEQ): - for j in range(cls.SLOTMAX): - st[:,i*cls.SLOTMAX+j] = (tag[:,i] >> (cls.BITS*j)) & cls.MASK - pass - pass - return st - - def __init__(self, path=PATH): - path = os.path.expandvars(path) - lines = open(path, "r").read().splitlines() - self.path = path - self.lines = lines - self.items = [] - self.parse() - - def find_item(self, name): - for item in self.items: - if item.name == name: return item - pass - return None - - def parse(self): - self.code2item = odict() - self.name2item = odict() - for line in self.lines: - enum_match = self.enum_ptn.match(line) - note_match = self.note_ptn.match(line) - if enum_match: - name, val = enum_match.groups() - pfx = "stag_" - assert name.startswith(pfx) - sname = name[len(pfx):] - code = int(val) - item = stag_item(code, sname, "") - self.items.append(item) - self.code2item[code] = item - self.name2item[sname] = item - log.debug("%40s : name:%20s sname:%10s val:%10s code:%d " % (line,name,sname,val, code) ) - elif note_match: - name, note = note_match.groups() - item = self.find_item(name) - assert not item is None - item.note = note - log.debug(" note %10s : %s " % (name, note)) - pass - pass - pass - - def old_label(self, st): - d = self.d - label_ = lambda _:repr(d.get(_,stag_item.Placeholder())) - ilabel_ = lambda _:"%2d : %s" % ( _, label_(st[_])) - return "\n".join(map(ilabel_, range(len(st)))) - - - def __call__(self, code): - return self.code2item.get(code,stag_item.Placeholder()) - - def label(self, st, fl=None): - """ - :param st: array of stag enumeration codes - :param fl: None or array of flat uniform rands of shape shape as st - """ - if not fl is None: - assert st.shape == fl.shape - pass - - lines = [] - num_zero = 0 - - for i in range(len(st)): - code = st[i] - flat = fl[i] if not fl is None else None - item = self(code) - it = repr(item) - assert item.code == code - if code == st[0] and i > 0: - lines.append("") - pass - label = "%2d : %s " % (i, it) if fl is None else "%2d : %10.4f : %s" % (i, flat, it) - lines.append(label) - if code == 0: - num_zero += 1 - if num_zero == 2: break - pass - pass - return "\n".join(lines) - - def __str__(self): - return "\n".join(self.lines) - - def __repr__(self): - return "\n".join(list(map(repr,self.items))) - - - - -def test_label(): - tag = stag() - #print(tag) - print(repr(tag)) - - st = np.array([[ 1, 2, 9, 10, 1, 2, 11, 12, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], - [ 1, 2, 9, 10, 1, 2, 11, 12, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], - [ 1, 2, 9, 10, 1, 2, 11, 12, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]], dtype=np.uint8) - - print(tag.label(st[0,:10])) - - -def test_StepSplit(): - from numpy import array, uint8 - at = array( - [[ 1, 2, 9, 10, 1, 2, 9, 10, 1, 2, 11, 12, 0, 0, 0, 0], - [ 1, 2, 9, 10, 1, 2, 9, 10, 1, 2, 11, 12, 0, 0, 0, 0], - [ 1, 2, 9, 10, 1, 2, 9, 10, 1, 2, 11, 12, 0, 0, 0, 0]], dtype=uint8) - - x_ats = array( - [[[ 1, 2, 9, 10, 0, 0, 0, 0, 0, 0], - [ 1, 2, 9, 10, 0, 0, 0, 0, 0, 0], - [ 1, 2, 11, 12, 0, 0, 0, 0, 0, 0], - [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]], - - [[ 1, 2, 9, 10, 0, 0, 0, 0, 0, 0], - [ 1, 2, 9, 10, 0, 0, 0, 0, 0, 0], - [ 1, 2, 11, 12, 0, 0, 0, 0, 0, 0], - [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]], - - [[ 1, 2, 9, 10, 0, 0, 0, 0, 0, 0], - [ 1, 2, 9, 10, 0, 0, 0, 0, 0, 0], - [ 1, 2, 11, 12, 0, 0, 0, 0, 0, 0], - [ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]]], dtype=uint8) - - ats = stag.StepSplit(at) - assert np.all( ats == x_ats ) - - - - - - - -if __name__ == '__main__': - logging.basicConfig(level=logging.INFO) - - #test_label() - #test_StepSplit() - - #test_PFold() - - diff --git a/sysrap/stree.py b/sysrap/stree.py deleted file mode 100644 index 45b6440b37..0000000000 --- a/sysrap/stree.py +++ /dev/null @@ -1,594 +0,0 @@ -#!/usr/bin/env python -""" -stree.py -========== - -TODO: improve stree.py repr - -""" - -import os, numpy as np, logging, builtins -log = logging.getLogger(__name__) -from opticks.ana.fold import Fold -from opticks.sysrap.sfreq import sfreq -from opticks.sysrap.OpticksCSG import CSG_ - - -class STR(str): - """STR inherits from str and changes the repr to provide the str : useful for interactive ipython""" - def __repr__(self): - return str(self) - - -class sobject(object): - @classmethod - def Label(cls, spc=5, pfx=10): - prefix = " " * pfx - spacer = " " * spc - return prefix + spacer.join(cls.FIELD) - - @classmethod - def Fields(cls, bi=False): - kls = cls.__name__ - for i, field in enumerate(cls.FIELD): - setattr(cls, field, i) - if bi:setattr(builtins, field, i) - pass - - @classmethod - def Type(cls): - cls.Fields() - kls = cls.__name__ - print("%s.Type()" % kls ) - for i, field in enumerate(cls.FIELD): - name = cls.DTYPE[i][0] - fieldname = "%s.%s" % (kls, field) - print(" %2d : %20s : %s " % (i, fieldname, name)) - pass - print("%s.Label() : " % cls.Label() ) - - @classmethod - def RecordsFromArrays(cls, a): - """ - :param a: ndarray - :return: np.recarray - """ - ra = np.core.records.fromarrays(a.T, dtype=cls.DTYPE ) - return ra - - @classmethod - def Doc(cls): - assert len(cls.FIELD) == len(cls.DTYPE) - lines = [] - for i in range(len(cls.FIELD)): - line = "%2d : %2s : %15s : %s " % (i, cls.FIELD[i], cls.DTYPE[i][0], cls.DTYPE[i][1] ) - lines.append(line) - pass - return STR("\n".join(lines)) - - - - -class sn(sobject): - """ - sn.h CSG constituent node, assuming the default WITH_CHILD - """ - DTYPE = [ - ('typecode', ' 0 - lvid = int(ii[0]) - if len(ii) > 1: - log.warning("multiple lvid %s correspond to arg:%s using first %s " % (str(ii), arg, lvid)) - pass - else: - print("arg %s unhandled %s " % (arg, type(arg))) - assert(0) - pass - return np.where( self.nds.lvid == lvid )[0] - - def find_lvid_node( self, arg, ordinal=0 ): - nn = self.find_lvid_nodes(arg) - return nn[ordinal] - - - def get_children(self, nidx, prepend_arg=False): - """ - This is a direct translation of the C++ stree::get_children - which will inevitably be slow in python. - HMM: how to do this in a more numpy way ? - """ - children = [] - if prepend_arg: - children.append(nidx) - pass - nd = self.nds[nidx] - assert nd.index == nidx - ch = nd.first_child - while ch > -1: - child = self.nds[ch] - assert child.parent == nd.index - children.append(child.index) - ch = child.next_sibling - pass - return children - - def get_progeny_r(self, nidx, progeny): - """ - :param nidx: node index - :param progeny: in/out list of node indices - - Recursively collects node indices of progeny beneath nidx, - not including nidx itself. - """ - children = self.get_children(nidx) - progeny.extend(children) - for nix in children: - self.get_progeny_r(nix, progeny) - pass - - def get_progeny(self, nidx): - progeny = [] - self.get_progeny_r(nidx, progeny) - return np.array(progeny) - - def get_soname(self, nidx): - lvid = self.get_lvid(nidx) - return self.f.soname[lvid] if lvid > -1 and lvid < len(self.f.soname) else None ; - - def get_sub(self, nidx): - return self.f.subs[nidx] if nidx > -1 else None ; - - def get_depth(self, nidx): - return self.nds[nidx].depth if nidx > -1 else -1 ; - - def get_parent(self, nidx): - return self.nds[nidx].parent if nidx > -1 else -1 ; - - def get_lvid(self, nidx): - return self.nds[nidx].lvid if nidx > -1 else -1 ; - - def get_transform(self, nidx): - return self.f.trs[nidx] if nidx > -1 and nidx < len(self.f.trs) else None - - def get_ancestors(self, nidx): - ancestors = [] - parent = self.get_parent(nidx) - while parent > -1: - ancestors.append(parent) - parent = self.get_parent(parent) - pass - ancestors.reverse() - return ancestors - - def __repr__(self): - return repr(self.f) - def __str__(self): - return str(self.f) - - @classmethod - def DepthSpacer(cls, dep, depmax=15): - fmt = "|S%d" % depmax - spc = np.zeros( (depmax,), dtype=np.int8 ) - spc[:] = ord(" ") - spc[dep] = ord("+") - return spc.view(fmt)[0].decode() - - def desc_node(self, nidx, brief=False): - return self.desc_node_(nidx, self.sf, brief=brief) - - def desc_nodes(self, nodes, brief=False): - return "\n".join([self.desc_node(nix, brief=brief) for nix in nodes]) - - def desc_node_(self, nidx, sf, brief=False): - - nd = self.nds[nidx] - dep = self.get_depth(nidx) - spc = self.DepthSpacer(dep) - ndd = snode.Brief(nd) if brief else snode.Desc(nd) - sub = "skip-sub" # self.get_sub(nidx) - sfd = "skip-sfd" # "" if sf is None else sf.desc_key(sub) - son = "skip-son" # self.get_soname(nidx) - return " ".join([spc,ndd,sfd,son]) - - def desc_nodes_(self, nodes, sf): - return "\n".join( [self.desc_node_(nix, sf) for nix in nodes]) - - def make_freq(self, nodes): - assert type(nodes) is np.ndarray - if self.raw_subs is None: - path = os.path.join(self.f.base, "subs.txt") - self.raw_subs = np.loadtxt( path, dtype="|S32") - pass - ssub = self.raw_subs[nodes] - ssf = sfreq.CreateFromArray(ssub) - return ssf - - - def desc_boundary_stats(self): - """ - Counts of volumes with each boundary - - desc_boundary_stats - u_bd, n_bd = np.unique( st.nds.boundary, return_counts=True ) - 0 : 0 : 1 : Galactic///Galactic - 1 : 1 : 3 : Galactic///Rock - 2 : 2 : 1 : Rock///Galactic - 3 : 3 : 1 : Rock//Implicit_RINDEX_NoRINDEX_pDomeAir_pDomeRock/Air - 4 : 4 : 1 : Rock///Rock - 5 : 5 : 1 : Rock//Implicit_RINDEX_NoRINDEX_pExpHall_pExpRockBox/Air - 6 : 6 : 1 : Air/Implicit_RINDEX_NoRINDEX_pExpHall_pPoolCover//Steel - ... - 108 : 108 : 4997 : Pyrex/HamamatsuR12860_PMT_20inch_photocathode_mirror_logsurf/HamamatsuR12860_PMT_20inch_photocathode_mirror_logsurf/Vacuum - 109 : 109 : 4997 : Vacuum/HamamatsuR12860_PMT_20inch_dynode_plate_opsurface//Steel - 110 : 110 : 4997 : Vacuum/HamamatsuR12860_PMT_20inch_outer_edge_opsurface//Steel - 111 : 111 : 4997 : Vacuum/HamamatsuR12860_PMT_20inch_inner_edge_opsurface//Steel - - - """ - lines = [] - lines.append("desc_boundary_stats") - st = self - lines.append("u_bd, n_bd = np.unique( st.nds.boundary, return_counts=True )") - u_bd, n_bd = np.unique( st.nds.boundary, return_counts=True ) - - for i in range(len(u_bd)): - u = u_bd[i] - n = n_bd[i] - line = " %3d : %4d : %6d : %s " % (i, u, n, st.f.standard.bd_names[u] ) - lines.append(line) - pass - return STR("\n".join(lines)) - - - def desc_remainder(self): - lines = [] - lines.append("desc_remainder") - st = self - u_lv, n_lv = np.unique( st.rem.lvid, return_counts=True ) - lines.append("u_lv, n_lv = np.unique( st.rem.lvid, return_counts=True )") - - for i in range(len(u_lv)): - u = u_lv[i] - n = n_lv[i] - line = " %3d : %4d : %6d : %s " % (i, u, n, st.f.soname[u] ) - lines.append(line) - pass - return "\n".join(lines) - - - def desc_csg(self, lvid=105): - """ - In [8]: snd.Label(3,8), n[n[:,snd.lv] == 105] - Out[8]: - (' ix dp sx pt nc fc sx lv tc pm bb xf', - array([[483, 4, 0, 484, 0, -1, -1, 105, 105, 294, 294, -1], - [484, 3, 0, 486, 1, 483, 485, 105, 11, -1, -1, -1], - [485, 3, 1, 486, 0, -1, -1, 105, 103, 295, 295, 183], - [486, 2, 0, 488, 2, 484, 487, 105, 1, -1, -1, -1], - [487, 2, 1, 488, 0, -1, -1, 105, 105, 296, 296, 184], - [488, 1, 0, 495, 2, 486, 494, 105, 1, -1, -1, -1], - [489, 4, 0, 490, 0, -1, -1, 105, 105, 297, 297, -1], - [490, 3, 0, 492, 1, 489, 491, 105, 11, -1, -1, -1], - [491, 3, 1, 492, 0, -1, -1, 105, 103, 298, 298, 186], - [492, 2, 0, 494, 2, 490, 493, 105, 1, -1, -1, -1], - [493, 2, 1, 494, 0, -1, -1, 105, 105, 299, 299, 187], - [494, 1, 1, 495, 2, 492, -1, 105, 1, -1, -1, 188], - [495, 0, -1, -1, 2, 488, -1, 105, 3, -1, -1, -1]], dtype=int32)) - - In [9]: st.f.soname[105] - Out[9]: 'HamamatsuR12860Tail0x61b5500' - """ - csg = self.get_csg(lvid) - lines = [] - lines.append("desc_csg lvid:%d st.f.soname[%d]:%s " % (lvid,lvid,self.get_lv_soname(lvid))) - lines.append(snd.Label(3,8)) - lines.append("%s" % repr(csg)) - return "\n".join(lines) - - def get_csg(self, lvid): - st = self - n = st.f._csg.sn - return n[n[:,sn.lv] == lvid] - - def get_csg_typecode(self, lvid): - n = self.get_csg(lvid) - return n[:,sn.tc] - - - def get_numSolid(self): - return 1+len(self.factor) - - def get_numPrim(self, ridx): - return len(self.rem) if ridx == 0 else self.factor[ridx-1].subtree - - def get_lvid(self, ridx): - """ - :param ridx: - :return lvid: array of meshIdx - - In [3]: ridx = 1 ; st.nds.lvid[st.nds.repeat_index == ridx][:st.factor[ridx-1].subtree] - Out[3]: array([133, 131, 129, 130, 132], dtype=int32) - - In [4]: ridx = 2 ; st.nds.lvid[st.nds.repeat_index == ridx][:st.factor[ridx-1].subtree] - Out[4]: array([128, 118, 119, 127, 126, 120, 125, 121, 122, 123, 124], dtype=int32) - - In [9]: st.nds.lvid[st.nds.repeat_index == 0].shape - Out[9]: (3089,) - """ - st = self - 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] - - def get_lv_soname(self, lv): - return str(self.soname_[lv]) # former .decode("utf-8") - - def descSolid(self, ridx, detail=False): - """ - cf with CSGFoundry.descSolid - """ - numPrim = self.get_numPrim(ridx) - lvid = self.get_lvid(ridx) - u_lvid, n_lvid = np.unique(lvid, return_counts=True ) - n_lvid_one = np.all( n_lvid == 1 ) - p_lvid = lvid if n_lvid_one else u_lvid # present in original order when not "unique" summarizing - pass - lines = [] - lines.append("stree.descSolid ridx %3d numPrim %5d lvid %s n_lvid_one %d" % (ridx, numPrim, str(lvid), n_lvid_one)) - if detail: - lines.append("") - for pass_ in [0,1]: - for i in range(len(p_lvid)): - ulv = p_lvid[i] - nlv = n_lvid[i] - lvn = self.get_lv_soname(ulv) - csg = self.get_csg(ulv) - tc = self.get_csg_typecode(ulv) - assert len(csg) == len(tc) - tcn = " ".join(list(map(lambda _:"%d:%s"%(_,CSG_.desc(_)), tc))) - if pass_ == 1: - lines.append("") - pass - lines.append(" lv:%3d nlv:%2d %50s csg %2d tcn %s " % (ulv, nlv, lvn, len(csg), tcn )) - if pass_ == 1: - lines.append(self.desc_csg(ulv)) - pass - pass - pass - pass - return "\n".join(lines) - - def descSolids(self, detail=False): - lines = [] - numSolid = self.get_numSolid() - lines.append("stree.descSolids numSolid:%d detail:%d " % (numSolid,detail) ) - - q_ridx = int(os.environ.get("RIDX","-1")) - for ridx in range(numSolid): - if q_ridx > -1 and ridx != q_ridx: continue - lines.append(self.descSolid(ridx, detail=detail)) - if detail: - lines.append("") - pass - pass - return "\n".join(lines) - - - - -if __name__ == '__main__': - logging.basicConfig(level=logging.INFO) - f = Fold.Load(symbol="f") - - st = stree(f) - print(repr(st)) - - nidx = os.environ.get("NIDX", 13) - - ancestors = st.get_ancestors(nidx) - print("NIDX %d ancestors: %s " % (nidx, repr(ancestors)) ) - print("st.desc_nodes(ancestors)") - print(st.desc_nodes(ancestors)) - - so = os.environ.get("SO", "sBar") - lvs = st.find_lvid_(so) - print("SO %s lvs %s" % (so, str(lvs))) - - for lv in lvs: - bb = st.find_lvid_nodes(lv) - b = int(bb[0]) - print("lv:%d bb=st.find_lvid_nodes(lv) bb:%s b:%s " % (lv, str(bb),b)) - - anc = st.get_ancestors(b) - print("b:%d anc=st.get_ancestors(b) anc:%s " % (b, str(anc))) - - print("st.desc_nodes(anc, brief=True))") - print(st.desc_nodes(anc, brief=True)) - print("st.desc_nodes([b], brief=True))") - print(st.desc_nodes([b], brief=True)) - pass - - - -if 0: - progeny = st.get_progeny(nidx) - print("NIDX %d progeny: %s " % (nidx, repr(progeny)) ) - #print("st.desc_nodes(progeny)") - #print(st.desc_nodes(progeny)) - - psf = st.make_freq(progeny) - print("st.desc_nodes_(progeny, psf)") - print(st.desc_nodes_(progeny, psf)) - - np.set_printoptions(edgeitems=600) - - print("st.f.trs[progeny].reshape(-1,16)") - print(st.f.trs[progeny].reshape(-1,16)) diff --git a/sysrap/xfold.py b/sysrap/xfold.py deleted file mode 100644 index e987d3a208..0000000000 --- a/sysrap/xfold.py +++ /dev/null @@ -1,174 +0,0 @@ -#!/usr/bin/env python - -import numpy as np -import os, re, logging -log = logging.getLogger(__name__) - -from opticks.ana.fold import Fold -from opticks.sysrap.stag import stag -from opticks.u4.U4Stack import U4Stack -from opticks.ana.p import * - - -class XFold(object): - tag = stag() - stack = U4Stack() - - @classmethod - def BaseSymbol(cls, xf): - """ - :param xf: Fold instance - :return symbol: "A" or "B" : A for Opticks, B for Geant4 - """ - CX = xf.base.find("CX") > -1 - U4 = xf.base.find("U4") > -1 - assert CX ^ U4 # exclusive-OR - symbol = "A" if CX else "B" - return symbol - - @classmethod - def Ident(cls, x): - bsym = cls.BaseSymbol(x) - ident = cls.tag if bsym == "A" else cls.stack - return ident - - def __init__(self, x, symbol=None): - """ - :param x: Fold instance, eg either a or b - """ - t = stag.Unpack(x.tag) if hasattr(x,"tag") else None - f = getattr(x, "flat", None) - n = stag.NumStarts(t) if not t is None else None - log.info("XFold before stag.StepSplit") - ts,fs = stag.StepSplit(t,x.flat) if not t is None else None - log.info("XFold after stag.StepSplit") - ident = self.Ident(x) - - xsymbol = self.BaseSymbol(x) - if xsymbol == "B": - t2 = self.stack.make_stack2tag_mapped(t) - ts2 = self.stack.make_stack2tag_mapped(ts) - elif xsymbol == "A": - t2 = self.stack.make_tag2stack_mapped(t) - ts2 = self.stack.make_tag2stack_mapped(ts) - else: - assert 0 - pass - if symbol != xsymbol: - log.error("using unconventional symbol %s xsymbol %s " % (symbol, xsymbol)) - pass - #assert symbol == xsymbol - - self.x = x # Fold instance, called x because it is usually "a" or "b" - self.t = t # (num_photon, SLOTS) : unpacked consumption tag/stack enumeration integers - self.f = f # (num_photon, SLOTS) : uniform rand - self.n = n # (num_photon,) : number of steps - self.t2 = t2 # (num_photon, SLOTS) : native enumeration mapped to the other enumeration - - self.ts = ts # (num_photon, 7, 10) : 7 is example max_steps which should match A.n.max(), 10 is example max rand per step input - self.fs = fs # (num_photon, 7, 10) : flat random values split by steps - self.ts2 = ts2 # (num_photon, 7, 10) : native enumeration mapped to the other enumeration - - self.ident = ident # A:stag OR B:U4Stack instance for providing names to enumeration codes - self.symbol = symbol - self.xsymbol = xsymbol - self.idx = 0 - self.flavor = "" - - def __call__(self, idx): - self.idx = idx - self.flavor = "call" - return self - - def __getitem__(self, idx): - self.idx = idx - self.flavor = "getitem" - return self - - def header(self): - idx = self.idx - seqhis = seqhis_(self.x.seq[idx,0]) - return "%s(%d) : %s" % (self.symbol, idx, seqhis) - - def rbnd__(self): - return boundary___(self.x.record[self.idx]) - def rori__(self): - return orient___(self.x.record[self.idx]) - def rpri__(self): - return primIdx___(self.x.record[self.idx]) - def rins__(self): - return instanceId__(self.x.record[self.idx]) - - - def rbnd_(self): - bb_ = self.rbnd__() - oo_ = self.rori__() - pp_ = self.rpri__() - ii_ = self.rins__() - assert bb_.shape == oo_.shape - assert bb_.shape == pp_.shape - assert bb_.shape == ii_.shape - - wp = np.where( bb_ > 0 ) # mask boundary zero to skip unsets - bb = bb_[wp] - oo = oo_[wp] - pp = pp_[wp] - ii = ii_[wp] - - pm = ["+","-"] - lines = [ "%s %-40s %-50s"%(pm[oo[i]],cf.sim.bndnamedict.get(bb[i]), cf.primIdx_meshname_dict.get(pp[i])) for i in range(len(bb))] - - lines += [""] - lines += ["- : against the normal (ie inwards from omat to imat)"] - lines += ["+ : with the normal (ie outwards from imat to omat)"] - return lines - - - def rbnd(self): - return "\n".join(self.rbnd_()) - - def content(self): - lines = [] - members = "t t2 n ts fs ts2".split() - for mem in members: - arr = getattr(self, mem) - name = "%s.%s" % ( self.symbol, mem) - line = "%10s : %s " % (name, str(arr.shape)) - lines.append(line) - pass - return "\n".join(lines) - - def body(self): - return self.ident.label(self.t[self.idx],self.f[self.idx]) - - def identification(self): - return "%s : %s" % (self.symbol, self.x.base) - - def call_repr(self): - return "\n".join([self.header(), self.content(), self.body()]) - - def getitem_repr(self): - return "\n".join([self.header(), self.rbnd()]) - - def __repr__(self): - if self.flavor == "call": - rep = self.call_repr() - elif self.flavor == "getitem": - rep = self.getitem_repr() - else: - rep = self.call_repr() - pass - return rep - - - - -if __name__ == '__main__': - logging.basicConfig(level=logging.INFO) - - a = Fold.Load("$A_FOLD", symbol="a") - b = Fold.Load("$B_FOLD", symbol="b") - A = XFold(a, symbol="A") - B = XFold(b, symbol="B") - - From ceb31e3bc3d991add88d44cca1f15047b3897b8d Mon Sep 17 00:00:00 2001 From: Dmitri Smirnov Date: Thu, 7 May 2026 18:42:53 -0400 Subject: [PATCH 2/5] chore: remove sysrap/*.sh --- sysrap/SABTest.sh | 58 ------------------ sysrap/SCenterExtentGenstep.sh | 13 ---- sysrap/SNameTest.sh | 65 -------------------- sysrap/SProp.sh | 5 -- sysrap/dv.sh | 7 --- sysrap/go.sh | 41 ------------- sysrap/npd.sh | 28 --------- sysrap/populate_gl.sh | 8 --- sysrap/sc4u.sh | 14 ----- sysrap/sevt_load.sh | 19 ------ sysrap/sframe.sh | 6 -- sysrap/sfreq.sh | 6 -- sysrap/smonitor.sh | 108 --------------------------------- sysrap/stag.sh | 7 --- sysrap/stree.sh | 10 --- sysrap/xfold.sh | 7 --- 16 files changed, 402 deletions(-) delete mode 100644 sysrap/SABTest.sh delete mode 100755 sysrap/SCenterExtentGenstep.sh delete mode 100644 sysrap/SNameTest.sh delete mode 100755 sysrap/SProp.sh delete mode 100755 sysrap/dv.sh delete mode 100755 sysrap/go.sh delete mode 100755 sysrap/npd.sh delete mode 100755 sysrap/populate_gl.sh delete mode 100755 sysrap/sc4u.sh delete mode 100755 sysrap/sevt_load.sh delete mode 100755 sysrap/sframe.sh delete mode 100755 sysrap/sfreq.sh delete mode 100755 sysrap/smonitor.sh delete mode 100755 sysrap/stag.sh delete mode 100755 sysrap/stree.sh delete mode 100755 sysrap/xfold.sh diff --git a/sysrap/SABTest.sh b/sysrap/SABTest.sh deleted file mode 100644 index 3df3af8973..0000000000 --- a/sysrap/SABTest.sh +++ /dev/null @@ -1,58 +0,0 @@ -#!/bin/bash -usage(){ cat << EOU -SABTest.sh : A-B comparison of persisted Opticks SEvt -======================================================== - -A:B comparison of SEvt and plotting is currently done by very -many Opticks scripts. Aims of SABTest.sh: - -1. replace most of the other A-B comparisons to decrease duplication -2. work from release, ie do not depend on having source tree - -Started from ~/j/InputPhotonsCheck/InputPhotonsCheck.sh - -EOU -} - -vv="BASH_SOURCE" -source $HOME/.opticks/GEOM/GEOM.sh -vv="$vv GEOM" -source $HOME/.opticks/GEOM/EVT.sh -vv="$vv AFOLD BFOLD" - -anascript=SABTest.py - -defarg="info_ls_pdb" -arg=${1:-$defarg} - -if [ "${arg/info}" != "$arg" ]; then - for v in $vv ; do printf "%30s : %s\n" "$v" "${!v}" ; done -fi - -if [ "${arg/ls}" != "$arg" ]; then - ff="AFOLD BFOLD" - for f in $ff ; do echo $f ${!f} && ls -alst ${!f} ; done -fi - -if [ "${arg/pdb}" != "$arg" ]; then - ${IPYTHON:-ipython} --pdb -i $(which $anascript) - [ $? -ne 0 ] && echo $BASH_SOURCE pdb error -fi - -if [ "${arg/ana}" != "$arg" ]; then - ${PYTHON:-python} $(which $anascript) - [ $? -ne 0 ] && echo $BASH_SOURCE ana error -fi - -if [ "${arg/chi2}" != "$arg" ] ; then - sseq_index_test.sh info_run_ana -fi - -if [ "${arg/dhi2}" != "$arg" ] ; then - # development version of "chi2" - unset sseq_index_test__DEBUG - #export sseq_index_test__DEBUG=1 - DEV=1 $HOME/opticks/sysrap/tests/sseq_index_test.sh info_build_run_ana -fi - - diff --git a/sysrap/SCenterExtentGenstep.sh b/sysrap/SCenterExtentGenstep.sh deleted file mode 100755 index 6bf7798b52..0000000000 --- a/sysrap/SCenterExtentGenstep.sh +++ /dev/null @@ -1,13 +0,0 @@ -#!/bin/bash -l - - -export CE=1000,0,0,1 - -SCenterExtentGenstepTest - -${IPYTHON:-ipython} -i SCenterExtentGenstep.py - - - - - diff --git a/sysrap/SNameTest.sh b/sysrap/SNameTest.sh deleted file mode 100644 index a01550729a..0000000000 --- a/sysrap/SNameTest.sh +++ /dev/null @@ -1,65 +0,0 @@ -#!/bin/bash -l -usage(){ cat << EOU - -QTYPE=C SNameTest virtual - - -Default match type is QTYPE=S meaning START:: - - epsilon:sysrap blyth$ SNameTest NNVTMCPPMTsMask_virtual0x - id.desc() - SName::desc numName 141 name[0] sTopRock_domeAir0x578ca70 name[-1] uni10x5832ff0 - findIndex NNVTMCPPMTsMask_virtual0x count 1 idx 110 - findIndices NNVTMCPPMTsMask_virtual0x idxs.size 1 SName::QTypeLabel START - descIndices - 110 : NNVTMCPPMTsMask_virtual0x5f5f900 - - -QTYPE=C uses CONTAIN name matching:: - - epsilon:sysrap blyth$ QTYPE=C SNameTest virtual - id.desc() - SName::desc numName 141 name[0] sTopRock_domeAir0x578ca70 name[-1] uni10x5832ff0 - findIndex virtual count 0 idx -1 - findIndices virtual idxs.size 3 SName::QTypeLabel CONTAIN - descIndices - 110 : NNVTMCPPMTsMask_virtual0x5f5f900 - 117 : HamamatsuR12860sMask_virtual0x5f50d40 - 134 : mask_PMT_20inch_vetosMask_virtual0x5f62e40 - - -QTYPE=E uses EXACT name matching:: - - epsilon:sysrap blyth$ QTYPE=E SNameTest virtual - id.desc() - SName::desc numName 141 name[0] sTopRock_domeAir0x578ca70 name[-1] uni10x5832ff0 - findIndex virtual count 0 idx -1 - findIndices virtual idxs.size 0 SName::QTypeLabel EXACT - descIndices - - epsilon:sysrap blyth$ QTYPE=E SNameTest NNVTMCPPMTsMask_virtual0x5f5f900 - id.desc() - SName::desc numName 141 name[0] sTopRock_domeAir0x578ca70 name[-1] uni10x5832ff0 - findIndex NNVTMCPPMTsMask_virtual0x5f5f900 count 1 idx 110 - findIndices NNVTMCPPMTsMask_virtual0x5f5f900 idxs.size 1 SName::QTypeLabel EXACT - descIndices - 110 : NNVTMCPPMTsMask_virtual0x5f5f900 - - epsilon:sysrap blyth$ QTYPE=E SNameTest NNVTMCPPMTsMask_virtual0x5f5f90 - id.desc() - SName::desc numName 141 name[0] sTopRock_domeAir0x578ca70 name[-1] uni10x5832ff0 - findIndex NNVTMCPPMTsMask_virtual0x5f5f90 count 1 idx 110 - findIndices NNVTMCPPMTsMask_virtual0x5f5f90 idxs.size 0 SName::QTypeLabel EXACT - descIndices - - -EOU -} - - -SNameTest $* - - - - - diff --git a/sysrap/SProp.sh b/sysrap/SProp.sh deleted file mode 100755 index 0549110877..0000000000 --- a/sysrap/SProp.sh +++ /dev/null @@ -1,5 +0,0 @@ -#!/bin/bash -l - -${IPYTHON:-ipython} --pdb -i SProp.py - - diff --git a/sysrap/dv.sh b/sysrap/dv.sh deleted file mode 100755 index fea7fe752e..0000000000 --- a/sysrap/dv.sh +++ /dev/null @@ -1,7 +0,0 @@ -#!/bin/bash -l - -source ../bin/AB_FOLD.sh - -${IPYTHON:-ipython} --pdb -i dv.py $* - - diff --git a/sysrap/go.sh b/sysrap/go.sh deleted file mode 100755 index 43207ad518..0000000000 --- a/sysrap/go.sh +++ /dev/null @@ -1,41 +0,0 @@ -#!/bin/bash -l -## -## Copyright (c) 2019 Opticks Team. All Rights Reserved. -## -## This file is part of Opticks -## (see https://bitbucket.org/simoncblyth/opticks). -## -## Licensed under the Apache License, Version 2.0 (the "License"); -## you may not use this file except in compliance with the License. -## You may obtain a copy of the License at -## -## http://www.apache.org/licenses/LICENSE-2.0 -## -## Unless required by applicable law or agreed to in writing, software -## distributed under the License is distributed on an "AS IS" BASIS, -## WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. -## See the License for the specific language governing permissions and -## limitations under the License. -## - - -opticks- - -sdir=$(pwd) -name=$(basename $sdir) -#bdir=/tmp/$USER/opticks/$name/build -bdir=$(opticks-prefix)/build/$name - -#rm -rf $bdir -mkdir -p $bdir && cd $bdir && pwd - -cmake $sdir \ - -DCMAKE_BUILD_TYPE=Debug \ - -DCMAKE_PREFIX_PATH=$(opticks-prefix)/externals \ - -DCMAKE_INSTALL_PREFIX=$(opticks-prefix) \ - -DCMAKE_MODULE_PATH=$(opticks-home)/cmake/Modules - -make -make install - - diff --git a/sysrap/npd.sh b/sysrap/npd.sh deleted file mode 100755 index 59091884b8..0000000000 --- a/sysrap/npd.sh +++ /dev/null @@ -1,28 +0,0 @@ -#!/bin/bash -l - -msg="=== $BASH_SOURCE :" -arg=${1:-diff} -echo $msg arg [$arg] AVAILABLE CMDS : diff, cp -hhs="NP.hh NPU.hh" - -if [ "$arg" == "diff" ]; then - for hh in $hhs ; do - cmd="diff $HOME/np/$hh $hh" - echo $msg $cmd - eval $cmd - [ $? -ne 0 ] && echo $msg DIFF DETECTED && exit 1 - done -elif [ "$arg" == "cp" ]; then - for hh in $hhs ; do - cmd="cp $HOME/np/$hh $hh" - echo $msg $cmd - eval $cmd - [ $? -ne 0 ] && echo $msg COPY ERROR && exit 2 - done -fi - - -echo $msg suggested commit message : update NP.hh from https://github.com/simoncblyth/np - -exit 0 - diff --git a/sysrap/populate_gl.sh b/sysrap/populate_gl.sh deleted file mode 100755 index 48f2de9a17..0000000000 --- a/sysrap/populate_gl.sh +++ /dev/null @@ -1,8 +0,0 @@ -#!/bin/bash - -cd ~/opticks/sysrap -#cp -r ../examples/UseShaderSGLFW_SScene_encapsulated/gl/* gl/ -#cp -r ../examples/UseGeometryShader/rec_flying_point_persist gl/ - - - diff --git a/sysrap/sc4u.sh b/sysrap/sc4u.sh deleted file mode 100755 index 4303a6fa73..0000000000 --- a/sysrap/sc4u.sh +++ /dev/null @@ -1,14 +0,0 @@ -#!/bin/bash - -name=sc4u -gcc $name.cc -std=c++11 -lstdc++ -I. -I/usr/local/cuda/include -o /tmp/$name -[ $? -ne 0 ] && echo compile error && exit 1 - -/tmp/$name -[ $? -ne 0 ] && echo run error && exit 2 - -ipython -i -c "import numpy as np ; a = np.load('/tmp/p.npy') ; print(a.view(np.int8)) " -[ $? -ne 0 ] && echo ana error && exit 3 - -exit 0 - diff --git a/sysrap/sevt_load.sh b/sysrap/sevt_load.sh deleted file mode 100755 index 8f01a627c9..0000000000 --- a/sysrap/sevt_load.sh +++ /dev/null @@ -1,19 +0,0 @@ -#!/bin/bash -l - - -#nevt=0 # 0 for default single Fold loading -nevt=3 # gt 0 for concat loading -export NEVT=${NEVT:-$nevt} - -evt="000" -[ $NEVT -gt 0 ] && evt="%0.3d" - -export GEOM=V0J008 -export BASE=/tmp/$USER/opticks/GEOM/$GEOM/ntds2 -export FOLD=$BASE/ALL0/${EVT:-$evt} - -echo $BASH_FOLD : FOLD $FOLD NEVT $NEVT - -${IPYTHON:-ipython} --pdb -i sevt_load.py - - diff --git a/sysrap/sframe.sh b/sysrap/sframe.sh deleted file mode 100755 index a3641756d4..0000000000 --- a/sysrap/sframe.sh +++ /dev/null @@ -1,6 +0,0 @@ -#!/bin/bash -l - -export FOLD=/tmp/$USER/opticks/CSG/CSGFoundry_MakeCenterExtentGensteps_Test -${IPYTHON:-ipython} --pdb -i sframe.py - - diff --git a/sysrap/sfreq.sh b/sysrap/sfreq.sh deleted file mode 100755 index 87add41c62..0000000000 --- a/sysrap/sfreq.sh +++ /dev/null @@ -1,6 +0,0 @@ -#!/bin/bash -l - -export FOLD=/tmp/$USER/opticks/U4TreeTest - -${IPYTHON:-ipython} --pdb -i $(dirname $BASH_SOURCE)/sfreq.py - diff --git a/sysrap/smonitor.sh b/sysrap/smonitor.sh deleted file mode 100755 index 07688dcf17..0000000000 --- a/sysrap/smonitor.sh +++ /dev/null @@ -1,108 +0,0 @@ -#!/bin/bash -usage(){ cat << EOU -smonitor.sh : NVML based GPU memory monitor -============================================= - -Usage: - -1. start the monitor process in one tab:: - - ~/o/sysrap/smonitor.sh # builds and starts runloop - -2. start the GPU memory using job in another tab - -3. once the job has completed, ctrl-C the monitor process - which catches the SIGINT signal and saves the memory - profile into smonitor.npy array - -4. make a plot:: - - ~/o/sysrap/smonitor.sh grab - ~/o/sysrap/smonitor.sh ana - - -5. screen capture the plot with annotation - - ~/o/sysrap/smonitor.sh mpcap - PUB=cxs_min_igs_with_rg_dummy ~/o/sysrap/smonitor.sh mppub - - -EOU -} - -SDIR=$(dirname $(realpath $BASH_SOURCE)) - -name=smonitor - -vars="BASH_SOURCE SDIR PWD name TMP FOLD CUDA_PREFIX CUDA_TARGET bin" - -cuda_prefix=/usr/local/cuda-11.7 -CUDA_PREFIX=${CUDA_PREFIX:-$cuda_prefix} -CUDA_TARGET=$CUDA_PREFIX/targets/x86_64-linux - -tmp=/tmp/$USER/opticksGG -TMP=${TMP:-$tmp} -export FOLD=$TMP/$name -mkdir -p $FOLD - -bin=$FOLD/$name -script=$SDIR/$name.py - -cd $FOLD -LOGDIR=$PWD - -#stem=smonitor_okjob -stem=smonitor_cxs_min_igs -pub=no-PUB -export STEM=${STEM:-$stem} -export PUB=${PUB:-$pub} - -export smonitor__SLEEP_US=100000 # 0.1s - - -defarg="info_build_run" -arg=${1:-$defarg} - -if [ "${arg/info}" != "$arg" ]; then - for var in $vars ; do printf "%25s : %s \n" "$var" "${!var}" ; done -fi - -if [ "${arg/build}" != "$arg" ]; then - gcc $SDIR/$name.cc -std=c++11 -lstdc++ \ - -I$SDIR \ - -I$CUDA_TARGET/include \ - -L$CUDA_TARGET/lib \ - -lnvidia-ml \ - -o $bin - [ $? -ne 0 ] && echo $BASH_SOURCE build error && exit 1 -fi - -if [ "${arg/run}" != "$arg" ]; then - $bin - [ $? -ne 0 ] && echo $BASH_SOURCE run error && exit 2 -fi - -if [ "${arg/grab}" != "$arg" ]; then - source $OPTICKS_HOME/bin/rsync.sh $LOGDIR - [ $? -ne 0 ] && echo $BASH_SOURCE grab error && exit 3 -fi - -if [ "${arg/ana}" != "$arg" ]; then - ${IPYTHON:-ipython} --pdb -i $script - [ $? -ne 0 ] && echo $BASH_SOURCE ana error && exit 4 -fi - -if [ "$arg" == "mpcap" -o "$arg" == "mppub" ]; then - export CAP_BASE=$FOLD/figs - export CAP_REL=smonitor - export CAP_STEM=${STEM}_${PUB} - case $arg in - mpcap) source mpcap.sh cap ;; - mppub) source mpcap.sh env ;; - esac - if [ "$arg" == "mppub" ]; then - source epub.sh - fi -fi - -exit 0 diff --git a/sysrap/stag.sh b/sysrap/stag.sh deleted file mode 100755 index e55e79cb3a..0000000000 --- a/sysrap/stag.sh +++ /dev/null @@ -1,7 +0,0 @@ -#!/bin/bash -l - -export A_FOLD=/tmp/$USER/opticks/GeoChain/BoxedSphere/CXRaindropTest -export B_FOLD=/tmp/$USER/opticks/U4RecorderTest -export FOLD=$A_FOLD - -${IPYTHON:-ipython} --pdb -i stag.py diff --git a/sysrap/stree.sh b/sysrap/stree.sh deleted file mode 100755 index ccd36451d9..0000000000 --- a/sysrap/stree.sh +++ /dev/null @@ -1,10 +0,0 @@ -#!/bin/bash -l - - -export FOLD=/tmp/$USER/opticks/U4TreeTest - -${IPYTHON:-ipython} --pdb -i $(dirname $BASH_SOURCE)/stree.py - - - - diff --git a/sysrap/xfold.sh b/sysrap/xfold.sh deleted file mode 100755 index b980d6f265..0000000000 --- a/sysrap/xfold.sh +++ /dev/null @@ -1,7 +0,0 @@ -#!/bin/bash -l - -export A_FOLD=/tmp/$USER/opticks/GeoChain/BoxedSphere/CXRaindropTest -export B_FOLD=/tmp/$USER/opticks/U4RecorderTest -export FOLD=$A_FOLD - -${IPYTHON:-ipython} --pdb -i xfold.py From b82a2c8100bfbc7c274e469fb31a4abe788d16c4 Mon Sep 17 00:00:00 2001 From: Dmitri Smirnov Date: Thu, 7 May 2026 18:43:49 -0400 Subject: [PATCH 3/5] chore: remove qudarap/*.py --- qudarap/__init__.py | 0 1 file changed, 0 insertions(+), 0 deletions(-) delete mode 100644 qudarap/__init__.py diff --git a/qudarap/__init__.py b/qudarap/__init__.py deleted file mode 100644 index e69de29bb2..0000000000 From 11741e5c7b0779c09edd8f4834f9e4d5166ccb26 Mon Sep 17 00:00:00 2001 From: Dmitri Smirnov Date: Thu, 7 May 2026 18:43:58 -0400 Subject: [PATCH 4/5] chore: remove qudarap/*.sh --- qudarap/QBndTest.sh | 3 --- qudarap/QSimTest.sh | 0 qudarap/qgs.sh | 17 ----------------- 3 files changed, 20 deletions(-) delete mode 100755 qudarap/QBndTest.sh delete mode 100644 qudarap/QSimTest.sh delete mode 100755 qudarap/qgs.sh diff --git a/qudarap/QBndTest.sh b/qudarap/QBndTest.sh deleted file mode 100755 index 68fadf8a16..0000000000 --- a/qudarap/QBndTest.sh +++ /dev/null @@ -1,3 +0,0 @@ -#!/bin/bash -l - -${IPYTHON:-ipython} --pdb -i tests/QBndTest.py diff --git a/qudarap/QSimTest.sh b/qudarap/QSimTest.sh deleted file mode 100644 index e69de29bb2..0000000000 diff --git a/qudarap/qgs.sh b/qudarap/qgs.sh deleted file mode 100755 index bc43fb067e..0000000000 --- a/qudarap/qgs.sh +++ /dev/null @@ -1,17 +0,0 @@ -#!/bin/bash -l - -name=qgs -gcc $name.cc \ - -std=c++11 \ - -I/usr/local/cuda/include \ - -I/usr/local/opticks/include/SysRap \ - -lstdc++ -o /tmp/$name - -[ $? -ne 0 ] && echo compile error && exit 1 - - -/tmp/$name -[ $? -ne 0 ] && echo run error && exit 2 - -exit 0 - From 315f14b2b48e88c9173efa1c85f4ab3789afc17c Mon Sep 17 00:00:00 2001 From: Dmitri Smirnov Date: Thu, 4 Jun 2026 14:49:47 -0400 Subject: [PATCH 5/5] chore(tests): remove obsolete analysis and helper scripts - delete legacy Python and shell utilities across CSG, CSGOptiX, g4cx, sysrap, u4, and qudarap - remove launcher wrappers that only referenced deleted scripts - update QSimTest.sh to stop defaulting to fake_propagate.py --- CSG/CSGFoundry.py | 1170 ----------------- CSG/CSGFoundryTest.sh | 18 - CSG/csg.sh | 96 -- CSG/csg_geochain.sh | 368 ------ CSG/foundry.py | 259 ---- CSG/foundry.sh | 5 - CSG/tests/CSGFoundryAB.py | 626 --------- CSG/tests/CSGFoundryLoadTest.py | 184 --- ...SGFoundry_MakeCenterExtentGensteps_Test.py | 56 - CSG/tests/CSGIntersectComparisonTest.py | 45 - CSG/tests/CSGIntersectSolidTest.py | 650 --------- CSG/tests/CSGNodeScanTest.py | 36 - CSG/tests/CSGSimtraceTest.py | 153 --- CSGOptiX/cxs_min.py | 177 --- CSGOptiX/cxs_min.sh | 739 ----------- CSGOptiX/cxs_min_AB.py | 53 - CSGOptiX/cxs_min_igs.sh | 227 ---- CSGOptiX/cxs_min_scan.sh | 58 - CSGOptiX/cxt_min.py | 601 --------- CSGOptiX/cxt_min.sh | 320 ----- CSGOptiX/cxt_precision.py | 172 --- CSGOptiX/cxt_precision.sh | 144 -- CSGOptiX/tests/CSGOptiXSimtraceTest.py | 262 ---- CSGOptiX/tests/CXRaindropTest.py | 80 -- g4cx/tests/G4CXSimtraceMinTest.py | 136 -- g4cx/tests/G4CXSimulateTest.py | 16 - g4cx/tests/G4CXTest.py | 266 ---- g4cx/tests/G4CXTest.sh | 231 ---- g4cx/tests/G4CXTest_GEOM.py | 229 ---- g4cx/tests/G4CXTest_GEOM.sh | 489 ------- g4cx/tests/G4CXTest_raindrop.py | 158 --- g4cx/tests/G4CXTest_raindrop.sh | 274 ---- g4cx/tests/G4CXTest_raindrop_CPU.sh | 16 - g4cx/tests/G4CXTest_raindrop_simtrace.py | 117 -- g4cx/tests/G4CXTest_raindrop_simtrace.sh | 28 - g4cx/tests/ab.py | 141 -- g4cx/tests/cf_G4CXSimtraceTest.py | 131 -- g4cx/tests/gx.py | 73 - g4cx/tests/gx.sh | 32 - g4cx/tests/sev.py | 60 - g4cx/tests/sev.sh | 33 - qudarap/tests/QEvtTest.sh | 88 -- qudarap/tests/QEvtTest_ALL.sh | 24 - qudarap/tests/QSimTest.sh | 4 +- qudarap/tests/fake_propagate.py | 143 -- sysrap/tests/CSG_Test.py | 48 - sysrap/tests/GetXF.py | 13 - sysrap/tests/SEvt_Lifecycle_Test.py | 13 - sysrap/tests/cfcf.py | 95 -- sysrap/tests/cfcf.sh | 19 - sysrap/tests/sframe_test.py | 12 - sysrap/tests/smath_test.py | 33 - sysrap/tests/stree_create_test.py | 22 - sysrap/tests/stree_load_test.py | 228 ---- sysrap/tests/stree_load_test_csg.py | 66 - sysrap/tests/stree_py_test.py | 20 - sysrap/tests/stree_sensor_test.py | 47 - sysrap/tests/stree_test.py | 109 -- u4/tests/PIDX.sh | 11 - u4/tests/U4RecorderTest.py | 207 --- u4/tests/U4RecorderTest_ab.py | 71 - u4/tests/U4SimtraceSimpleTest.py | 475 ------- u4/tests/U4SimtraceSimpleTest.sh | 125 -- u4/tests/U4SimtraceTest.py | 475 ------- u4/tests/U4SimtraceTest.sh | 251 ---- u4/tests/U4SimtraceTest_one_pmt.sh | 12 - u4/tests/U4SimtraceTest_two_pmt.sh | 34 - u4/tests/U4SimtraceTest_two_pmt_cf.sh | 38 - u4/tests/U4SimulateTest.sh | 323 ----- u4/tests/U4SimulateTest_cf.py | 126 -- u4/tests/U4SimulateTest_ck.py | 72 - u4/tests/U4SimulateTest_fk.py | 102 -- u4/tests/U4SimulateTest_mt.py | 4 +- u4/tests/U4SimulateTest_ph.py | 5 +- u4/tests/U4SimulateTest_pr.py | 188 --- u4/tests/U4TreeCreateSSimTest.py | 36 - u4/tests/U4TreeCreateSSimTest.sh | 107 -- u4/tests/U4TreeCreateTest.py | 74 -- u4/tests/U4TreeCreateTest.sh | 104 -- u4/tests/cf.sh | 41 - u4/tests/ck.sh | 3 - u4/tests/fk.sh | 5 - u4/tests/mt.sh | 8 - u4/tests/ph.sh | 11 - u4/tests/pr.sh | 33 - u4/tests/tt.sh | 17 - u4/tests/viz.sh | 83 -- 87 files changed, 5 insertions(+), 12949 deletions(-) delete mode 100644 CSG/CSGFoundry.py delete mode 100755 CSG/CSGFoundryTest.sh delete mode 100755 CSG/csg.sh delete mode 100755 CSG/csg_geochain.sh delete mode 100755 CSG/foundry.py delete mode 100755 CSG/foundry.sh delete mode 100644 CSG/tests/CSGFoundryAB.py delete mode 100644 CSG/tests/CSGFoundryLoadTest.py delete mode 100644 CSG/tests/CSGFoundry_MakeCenterExtentGensteps_Test.py delete mode 100644 CSG/tests/CSGIntersectComparisonTest.py delete mode 100644 CSG/tests/CSGIntersectSolidTest.py delete mode 100644 CSG/tests/CSGNodeScanTest.py delete mode 100644 CSG/tests/CSGSimtraceTest.py delete mode 100644 CSGOptiX/cxs_min.py delete mode 100755 CSGOptiX/cxs_min.sh delete mode 100644 CSGOptiX/cxs_min_AB.py delete mode 100755 CSGOptiX/cxs_min_igs.sh delete mode 100755 CSGOptiX/cxs_min_scan.sh delete mode 100644 CSGOptiX/cxt_min.py delete mode 100755 CSGOptiX/cxt_min.sh delete mode 100644 CSGOptiX/cxt_precision.py delete mode 100755 CSGOptiX/cxt_precision.sh delete mode 100644 CSGOptiX/tests/CSGOptiXSimtraceTest.py delete mode 100644 CSGOptiX/tests/CXRaindropTest.py delete mode 100644 g4cx/tests/G4CXSimtraceMinTest.py delete mode 100644 g4cx/tests/G4CXSimulateTest.py delete mode 100644 g4cx/tests/G4CXTest.py delete mode 100755 g4cx/tests/G4CXTest.sh delete mode 100644 g4cx/tests/G4CXTest_GEOM.py delete mode 100755 g4cx/tests/G4CXTest_GEOM.sh delete mode 100644 g4cx/tests/G4CXTest_raindrop.py delete mode 100755 g4cx/tests/G4CXTest_raindrop.sh delete mode 100755 g4cx/tests/G4CXTest_raindrop_CPU.sh delete mode 100644 g4cx/tests/G4CXTest_raindrop_simtrace.py delete mode 100755 g4cx/tests/G4CXTest_raindrop_simtrace.sh delete mode 100644 g4cx/tests/ab.py delete mode 100644 g4cx/tests/cf_G4CXSimtraceTest.py delete mode 100644 g4cx/tests/gx.py delete mode 100755 g4cx/tests/gx.sh delete mode 100644 g4cx/tests/sev.py delete mode 100755 g4cx/tests/sev.sh delete mode 100755 qudarap/tests/QEvtTest.sh delete mode 100755 qudarap/tests/QEvtTest_ALL.sh delete mode 100755 qudarap/tests/fake_propagate.py delete mode 100644 sysrap/tests/CSG_Test.py delete mode 100644 sysrap/tests/GetXF.py delete mode 100644 sysrap/tests/SEvt_Lifecycle_Test.py delete mode 100644 sysrap/tests/cfcf.py delete mode 100755 sysrap/tests/cfcf.sh delete mode 100644 sysrap/tests/sframe_test.py delete mode 100644 sysrap/tests/smath_test.py delete mode 100644 sysrap/tests/stree_create_test.py delete mode 100644 sysrap/tests/stree_load_test.py delete mode 100644 sysrap/tests/stree_load_test_csg.py delete mode 100644 sysrap/tests/stree_py_test.py delete mode 100644 sysrap/tests/stree_sensor_test.py delete mode 100644 sysrap/tests/stree_test.py delete mode 100755 u4/tests/PIDX.sh delete mode 100644 u4/tests/U4RecorderTest.py delete mode 100644 u4/tests/U4RecorderTest_ab.py delete mode 100644 u4/tests/U4SimtraceSimpleTest.py delete mode 100755 u4/tests/U4SimtraceSimpleTest.sh delete mode 100644 u4/tests/U4SimtraceTest.py delete mode 100755 u4/tests/U4SimtraceTest.sh delete mode 100755 u4/tests/U4SimtraceTest_one_pmt.sh delete mode 100755 u4/tests/U4SimtraceTest_two_pmt.sh delete mode 100755 u4/tests/U4SimtraceTest_two_pmt_cf.sh delete mode 100755 u4/tests/U4SimulateTest.sh delete mode 100644 u4/tests/U4SimulateTest_cf.py delete mode 100644 u4/tests/U4SimulateTest_ck.py delete mode 100644 u4/tests/U4SimulateTest_fk.py delete mode 100644 u4/tests/U4SimulateTest_pr.py delete mode 100644 u4/tests/U4TreeCreateSSimTest.py delete mode 100755 u4/tests/U4TreeCreateSSimTest.sh delete mode 100644 u4/tests/U4TreeCreateTest.py delete mode 100755 u4/tests/U4TreeCreateTest.sh delete mode 100755 u4/tests/cf.sh delete mode 100755 u4/tests/ck.sh delete mode 100755 u4/tests/fk.sh delete mode 100755 u4/tests/mt.sh delete mode 100755 u4/tests/ph.sh delete mode 100755 u4/tests/pr.sh delete mode 100755 u4/tests/tt.sh delete mode 100755 u4/tests/viz.sh diff --git a/CSG/CSGFoundry.py b/CSG/CSGFoundry.py deleted file mode 100644 index e1ea5c94c9..0000000000 --- a/CSG/CSGFoundry.py +++ /dev/null @@ -1,1170 +0,0 @@ -#!/usr/bin/env python -import os, sys, re, numpy as np, logging, datetime -from configparser import ConfigParser -from collections import OrderedDict, Counter -log = logging.getLogger(__name__) - -#from opticks.ana.key import keydir -from opticks.ana.fold import Fold, STR -from opticks.sysrap.OpticksCSG import CSG_ - - - -class KeyNameConfig(object): - """ - Parses ini file of the below form:: - - [key_boundary_regexp] - red = (?PWater/.*/Water$) - blue = Water///Acrylic$ - magenta = Water/Implicit_RINDEX_NoRINDEX_pOuterWaterInCD_pInnerReflector//Tyvek$ - yellow = DeadWater/Implicit_RINDEX_NoRINDEX_pDeadWater_pTyvekFilm//Tyvek$ - pink = Air/CDTyvekSurface//Tyvek$ - cyan = Tyvek//Implicit_RINDEX_NoRINDEX_pOuterWaterInCD_pCentralDetector/Water$ - orange = Water/Steel_surface//Steel$ - grey = - # empty value is special cased to mean all other boundary names - - Defaults path is $HOME/.opticks/GEOM/cxt_min.ini - """ - - - def __init__(self, _path ): - path = os.path.expandvars(_path) - cfg = ConfigParser() - cfg.read(path) - - self._path = _path - self.path = path - self.cfg = cfg - - - def __call__(self, _section): - sect = self.cfg[_section] - bdict = OrderedDict(sect) - - counts = Counter(bdict.values()) - duplicates = [val for val, count in counts.items() if count > 1] - - if duplicates: - log.fatal(f"CSGFoundry.py/KeyNameConfig : Found duplicated values: {duplicates}") - sys.exit(1) - pass - - return bdict - - - - - -class NameTable(object): - def __init__(self, symbol, ix, ixn, d_qix={}, KEY=None): - """ - :param symbol: eg btab or ptab - :param ix: large array of indices (eg boundaries or globalPrimIdx) - :param ixn: smallish array of names - :param d_qix: dict keyed on colors with arrays of ixn indices as values - """ - u, x, c = np.unique(ix, return_index=True, return_counts=True) - - if "KLUDGE" in os.environ: - kludge = u != 0xffff - u = u[kludge] - x = x[kludge] - c = c[kludge] - pass - - # form array with the keys that each boundary index is grouped into - akk = np.repeat("????????????????", len(u)) - for i in range(len(u)): - kk = [] - for k, qix in d_qix.items(): - if len(np.where(np.isin(qix, u[i]))[0]) == 1: - kk.append(k) - pass - pass - akk[i] = ",".join(kk) - pass - n = ixn[u] # IF THIS FAILS FROM 0xffff TRY KLUDGE=1 - o = np.argsort(c)[::-1] - - _tab = "np.c_[akk, u, c, x, n][o]" - tab = eval(_tab) - - self.symbol = symbol - self.akk = akk[o] - self.u = u[o] - self.x = x[o] - self.c = c[o] - self.n = n[o] - - self.KEY = KEY - - self.d_qix = d_qix - self._tab = _tab - self.tab = tab - self.ixn = ixn - - def get_names(self, k): - """ - print("\n".join(cf.primname[cf.cxtp.d_qix['thistle']])) - """ - return self.ixn[self.d_qix[k]] - - def dump_names(self, k): - names = self.get_names(k) - print("\n".join(names)) - - def __repr__(self): - """ - """ - vfmt = " %20s : %4d %8d %8d : %s " - sfmt = " %20s : %4s %8s %8s : %s " - tfmt = " %20s : %4s %8d %8s : %s " - - label = sfmt % ( ".akk", ".u", ".c", ".x", ".n" ) - lines = [] - lines.append("[cf.%s KEY:%s" % (self.symbol, self.KEY) ) - lines.append(label) - - if not self.KEY is None: - if self.KEY.startswith("~"): - KEY = self.KEY[1:] - invert = True - else: - KEY = self.KEY - invert = False - pass - KK = np.array(self.KEY.split(",")) - sel = np.where(np.isin(self.tab[:,0], KK, invert=invert))[0] - else: - sel = slice(None) - pass - - ctot = 0 - for row in self.tab[sel]: - akk = row[0] - u = int(row[1]) - c = int(row[2]) - x = int(row[3]) - n = row[4] - - line = vfmt % ( akk, u, c, x, n ) - ctot += c - lines.append(line) - pass - lines.append(label) - lines.append(tfmt % ("", "",ctot, "",".ctot")) - lines.append("]cf.%s" % self.symbol ) - self.ctot = ctot - return "\n".join(lines) - - - -class GroupedNameTable(object): - def __init__(self, symbol, ix, d_qix, d_anno, cf_ixn, lwid=100): - """ - :param symbol: identifier - :param ix: large array of indices - :param d_qix: dict keyed on colors with indices array values - :param d_anno: dict keyed on colors with label values - :param cf_ixn: small array of all names - :param lwid: width of the label field - """ - self.symbol = symbol - self.ix = ix - self.d_qix = d_qix - self.d_anno = d_anno - self.cf_ixn = cf_ixn - self.lwid = lwid - - wdict = {} - for k,qix in self.d_qix.items(): - _w = np.isin( ix, qix ) # bool array indicating which elem of ix are in the qix array - w = np.where(_w)[0] # indices of ix array with ix names matching the regexp - wdict[k] = w - pass - self.wdict = wdict - - def __repr__(self): - """ - """ - lines = [] - lines.append("[cf.%s___________________________________" % self.symbol ) - fmt = " %%15s %%%(lwid)ds nnn:%%4d %%s " % dict(lwid=self.lwid) - for k,qix in self.d_qix.items(): - _w = np.isin( self.ix, qix ) # bool array indicating which elem of bn are in the qbn array - w = np.where(_w)[0] # indices of bn array with boundaries matching the regexp - label = self.d_anno[k] - nn = self.cf_ixn[qix] - line = fmt % (k, label, len(nn), str(qix[:10])) - lines.append(line) - pass - lines.append("]cf.%s_______ cf.%s.d_qix ____________________________" % ( self.symbol, self.symbol) ) - return "\n".join(lines) - - - - - - -class CSGObject(object): - @classmethod - def Label(cls, spc=5, pfx=10): - prefix = " " * pfx - spacer = " " * spc - return prefix + spacer.join(cls.FIELD) - - @classmethod - def Fields(cls, bi=False): - kls = cls.__name__ - for i, field in enumerate(cls.FIELD): - setattr(cls, field, i) - if bi:setattr(builtins, field, i) - pass - - @classmethod - def Type(cls): - cls.Fields() - kls = cls.__name__ - print("%s.Type()" % kls ) - for i, field in enumerate(cls.FIELD): - name = cls.DTYPE[i][0] - fieldname = "%s.%s" % (kls, field) - print(" %2d : %20s : %s " % (i, fieldname, name)) - pass - print("%s.Label() : " % cls.Label() ) - - @classmethod - def RecordsFromArrays(cls, a): - """ - :param a: ndarray - :return: np.recarray - """ - ra = np.core.records.fromarrays(a.T, dtype=cls.DTYPE ) - return ra - - - -class CSGPrim(CSGObject): - DTYPE = [ - ('numNode', ' $TMP/mm.txt - - """ - PTN = re.compile("\\d+") - def __init__(self, path ): - mm = os.path.expandvars(path) - mm = open(mm, "r").read().splitlines() if os.path.exists(mm) else None - self.mm = mm - if mm is None: - log.fatal("missing %s, which is now a standard part of CSGFoundry " % path ) - sys.exit(1) - pass - - def imm(self, emm): - return list(map(int, self.PTN.findall(emm))) - - def label(self, emm): - imm = self.imm(emm) - labs = [self.mm[i] for i in imm] - lab = " ".join(labs) - - tilde = emm[0] == "t" or emm[0] == "~" - pfx = ( "NOT: " if tilde else " " ) - - if emm == "~0" or emm == "t0": - return "ALL" - elif imm == [1,2,3,4]: - return "ONLY PMT" - elif "," in emm: - return ( "EXCL: " if tilde else "ONLY: " ) + lab - else: - return lab - pass - - def __repr__(self): - return STR("\n".join(self.mm)) - - -class LV(object): - PTN = re.compile("\\d+") - def __init__(self, path): - lv = os.path.expandvars(path) - lv = open(lv, "r").read().splitlines() if os.path.exists(lv) else None - self.lv = lv - if lv is None: - log.fatal("missing %s, which is now a standard part of CSGFoundry " % path ) - sys.exit(1) - pass - - def ilv(self, elv): - return list(map(int, self.PTN.findall(elv))) - - def label(self, elv): - ilv = self.ilv(elv) - mns = [self.lv[i] for i in ilv] - mn = " ".join(mns) - tilde = elv[0] == "t" - lab = "" - if elv == "t": - lab = "ALL" - else: - lab = ( "EXCL: " if tilde else "ONLY: " ) + mn - pass - return lab - - def __str__(self): - return "\n".join(self.lv) - - def __repr__(self): - return "\n".join(["%3d:%s " % (i, self.lv[i]) for i in range(len(self.lv))]) - - - -class Deprecated_NPFold(object): - """ - HMM opticks.ana.fold.Fold looks more developed than this - TODO: eliminate all use of this - """ - INDEX = "NPFold_index.txt" - - @classmethod - def IndexPath(cls, fold): - return os.path.join(fold, cls.INDEX) - - @classmethod - def HasIndex(cls, fold): - return os.path.exists(cls.IndexPath(fold)) - - @classmethod - def Load(cls, *args, **kwa): - return cls(*args, **kwa) - - def __init__(self, *args, **kwa): - fold = os.path.join(*args) - self.fold = fold - self.has_index = self.HasIndex(fold) - - if self.has_index: - self.load_idx(fold) - else: - self.load_fts(fold) - pass - - def load_fts(self, fold): - assert 0 - - def load_idx(self, fold): - keys = open(self.IndexPath(fold),"r").read().splitlines() - aa = [] - kk = [] - ff = [] - subfold = [] - - for k in keys: - path = os.path.join(fold, k) - if k.endswith(".npy"): - assert os.path.exists(path) - log.info("loading path %s k %s " % (path, k)) - a = np.load(path) - aa.append(a) - kk.append(k) - else: - log.info("skip non .npy path %s k %s " % (path, k)) - pass - pass - self.kk = kk - self.aa = aa - self.ff = ff - self.subfold = subfold - self.keys = keys - - - def find(self, k): - return self.kk.index(k) if k in self.kk else -1 - - def has_key(self, k): - return self.find(k) > -1 - - def __repr__(self): - lines = [] - lines.append("NPFold %s " % self.fold ) - for i in range(len(self.kk)): - k = self.kk[i] - a = self.aa[i] - stem, ext = os.path.splitext(os.path.basename(k)) - path = os.path.join(self.fold, k) - lines.append("%10s : %20s : %s " % (stem, str(a.shape), k )) - pass - return "\n".join(lines) - - -class Deprecated_SSim(object): - """ - HMM: this adds little on top of ana.fold.Fold - TODO: get rid of it - """ - BND = "bnd.npy" - - - @classmethod - def Load(cls, simbase): - ssdir = os.path.join(simbase, "SSim") - log.info("SSim.Load simbase %s ssdir %s " % (simbase,ssdir)) - sim = Fold.Load(ssdir) if os.path.isdir(ssdir) else None - return sim - - def __init__(self, fold): - if self.has_key(self.BND): - bnpath = os.path.join(fold, "bnd_names.txt") - assert os.path.exists(bnpath) - bnd_names = open(bnpath,"r").read().splitlines() - self.bnd_names = bnd_names - pass - if hasattr(self, 'bnd_names'): # names list from NP bnd.names metadata - bndnamedict = SSim.namelist_to_namedict(self.bnd_names) - else: - bndnamedict = {} - pass - self.bndnamedict = bndnamedict - pass - - - -class CSGFoundry(object): - FOLD = os.path.expandvars("$TMP/CSG_GGeo/CSGFoundry") - FMT = " %10s : %20s : %s " - - - @classmethod - def namelist_to_namedict(cls, namelist): - nd = {} - if not namelist is None: - nd = dict(zip(range(len(namelist)),list(map(str,namelist)))) - pass - return nd - - - @classmethod - def CFBase(cls): - """ - Precedence order for which envvars are used to derive the CFBase folder is: - - 1. CFBASE_LOCAL - 2. CFBASE_ALT - 3. CFBASE - 4. GEOM - - This allows scripts such as CSGOptiX/cxt_min.sh to set envvars to control - the geometry arrays and meshname.txt etc that is loaded. - For example this could be used such that local analysis on laptop can use - a geometry rsynced from remote. - - """ - if "CFBASE_LOCAL" in os.environ: - cfbase= os.environ["CFBASE_LOCAL"] - note = "via CFBASE_LOCAL" - elif "CFBASE_ALT" in os.environ: - cfbase= os.environ["CFBASE_ALT"] - note = "via CFBASE_ALT" - elif "CFBASE" in os.environ: - cfbase= os.environ["CFBASE"] - note = "via CFBASE" - elif "GEOM" in os.environ: - KEY = os.path.expandvars("${GEOM}_CFBaseFromGEOM") - if KEY in os.environ: - cfbase = os.environ[KEY] - note = "via GEOM,%s" % KEY - else: - cfbase = os.path.expandvars("$HOME/.opticks/GEOM/$GEOM") - note = "via GEOM" - pass - else: - cfbase = None - note = "via NO envvars" - pass - if cfbase is None: - print("CSGFoundry.CFBase returning None, note:%s " % note ) - else: - print("CSGFoundry.CFBase returning [%s], note:[%s] " % (cfbase,note) ) - pass - return cfbase - - @classmethod - def Load(cls, cfbase_=None, symbol=None): - """ - :param cfbase_: typically None, but available when debugging eg comparing two geometries - """ - if cfbase_ is None: - cfbase = cls.CFBase() - else: - cfbase = os.path.expandvars(cfbase_) - pass - log.info("cfbase:%s " % cfbase) - if cfbase is None or not os.path.exists(os.path.join(cfbase, "CSGFoundry")): - print("ERROR CSGFoundry.CFBase returned None OR non-existing CSGFoundry dir so cannot CSGFoundry.Load" ) - return None - pass - assert not cfbase is None - cf = cls(fold=os.path.join(cfbase, "CSGFoundry")) - cf.symbol = symbol - return cf - - @classmethod - def FindDirUpTree(cls, origpath, name="CSGFoundry"): - """ - :param origpath: to traverse upwards looking for sibling dirs with *name* - :param name: sibling directory name to look for - :return full path to *name* directory - """ - elem = origpath.split("/") - found = None - for i in range(len(elem),0,-1): - path = "/".join(elem[:i]) - cand = os.path.join(path, name) - log.debug(cand) - if os.path.isdir(cand): - found = cand - break - pass - pass - return found - - @classmethod - def FindDigest(cls, path): - if path.find("*") > -1: # eg with a glob pattern filename element - base = os.path.dirname(path) - else: - base = path - pass - base = os.path.expandvars(base) - return cls.FindDigest_(base) - - @classmethod - def FindDigest_(cls, path): - """ - :param path: Directory path in which to look for a 32 character digest path element - :return 32 character digest or None: - """ - hexchar = "0123456789abcdef" - digest = None - for elem in path.split("/"): - if len(elem) == 32 and set(elem).issubset(hexchar): - digest = elem - pass - pass - return digest - - def __init__(self, fold): - self.load(fold) - self.meshnamedict = self.namelist_to_namedict(self.meshname) - self.primIdx_meshname_dict = self.make_primIdx_meshname_dict() - - self.mokname = "zero one two three four five six seven eight nine".split() - self.moknamedict = self.namelist_to_namedict(self.mokname) - self.insnamedict = {} - - self.lv = LV(os.path.join(fold, "meshname.txt")) - self.mm = MM(os.path.join(fold, "mmlabel.txt")) - - sim = Fold.Load(fold, "SSim") - - try: - bdn = sim.stree.standard.bnd_names - except AttributeError: - bdn = None - pass - if bdn is None: log.fatal("CSGFoundry fail to access sim.stree.standard.bnd_names : geometry incomplete" ) - if type(bdn) is np.ndarray: sim.bndnamedict = self.namelist_to_namedict(bdn) - pass - - self.kncfg = KeyNameConfig("$HOME/.opticks/GEOM/cxt_min.ini") - - self.bdn_config = self.kncfg("key_boundary_regexp") - self.bdn_config_note = self.kncfg("key_boundary_regexp_note") - - self.prn_config = self.kncfg("key_prim_regexp") - self.prn_config_note = self.kncfg("key_prim_regexp_note") - - self.bdn = bdn - self.sim = sim - - self.npa = self.node.reshape(-1,16)[:,0:6] - self.nbd = self.node.view(np.int32)[:,1,2] - self.nix = self.node.view(np.int32)[:,1,3] - self.nbb = self.node.reshape(-1,16)[:,8:14] - self.ntc = self.node.view(np.int32)[:,3,2] - self.ncm = self.node.view(np.uint32)[:,3,3] >> 31 # node complement - self.ntr = self.node.view(np.uint32)[:,3,3] & 0x7fffffff # node tranIdx+1 - - self.pnn = self.prim.view(np.int32)[:,0,0] # prim num node - self.pno = self.prim.view(np.int32)[:,0,1] # prim node offset - - - - self.pto = self.prim.view(np.int32)[:,0,2] # prim tran offset - self.ppo = self.prim.view(np.int32)[:,0,3] # prim plan offset - - self.crn = self.pno[self.pnn>1] # node indices of compound roots - self.crn_subnum = self.node.view(np.int32)[self.crn,0,0] - self.crn_suboff = self.node.view(np.int32)[self.crn,0,1] - - - self.psb = self.prim.view(np.int32)[:,1,0] # prim sbtIndexOffset - self.plv = self.prim.view(np.int32)[:,1,1] # prim lvid/meshIdx - self.prx = self.prim.view(np.int32)[:,1,2] # prim ridx/repeatIdx - self.pix = self.prim.view(np.int32)[:,1,3] # prim idx - - self.pbb = self.prim.reshape(-1,16)[:,8:14] - - - self.snp = self.solid[:,1,0].view(np.int32) # solid numPrim - self.spo = self.solid[:,1,1].view(np.int32) # solid primOffset - self.sce = self.solid[:,2].view(np.float32) - - - - def simtrace_boundary_analysis(self, bn, KEY=os.environ.get("KEY", None) ): - """ - Group simtrace intersects according to their boundary indices - using groups specified by regexp matching boundary names. - - Typically would have a few hundred boundary names - in cf.bdn but potentially millions of simtrace intersects - in the bn array. - - :param bn: large array of boundary indices - :param KEY: color key OR None - :return cxtb,btab: - """ - d_qbn, d_anno = self.Dict_find_name_indices_re_match(self.bdn, self.bdn_config, self.bdn_config_note ) - btab = NameTable("btab", bn, self.bdn, d_qbn, KEY) - cxtb = GroupedNameTable("cxtb", bn, d_qbn, d_anno, self.bdn, lwid=100) - self.btab = btab - self.cxtb = cxtb - return cxtb, btab - - def simtrace_prim_analysis(self, pr, KEY=os.environ.get("KEY", None) ): - """ - Group simtrace intersects according to their primitive indices - using groups specified by regexp matching solid names of the prim. - - :param pr: large array of primitive indices - :param KEY: - :return cxtp, ptab: GroupedNameTable, NameTable - """ - d_qpr, d_anno = self.Dict_find_name_indices_re_match(self.primname, self.prn_config, self.prn_config_note ) - ptab = NameTable("ptab", pr, self.primname, d_qpr, KEY) - cxtp = GroupedNameTable("cxtp", pr, d_qpr, d_anno, self.primname, lwid=50) - self.ptab = ptab - self.cxtp = cxtp - return cxtp, ptab - - @classmethod - def Dict_find_name_indices_re_match(cls, names, _nameptn_dict, _namenote_dict ): - """ - :param names: array of names - :param _nameptn_dict: label keys, regexp pattern string values - :return d: dict with same label keys and arrays of matching names indices - - Names unmatched by the provided regexp values are - grouped within the key with a blank value if provided. - """ - d = {} - anno = {} - matched = np.array([], dtype=np.int64 ) - unmatched_k = None - for k, v in _nameptn_dict.items(): - if v == '': - unmatched_k = k # usaully grey - continue - pass - qnm, nn, label = cls.Find_name_indices_re_match(names, v) - matched = np.concatenate( (matched, qnm) ) - k_note = _namenote_dict.get(k,"") - lab = "%100s ## %s" % ( label, k_note ) - d[k] = qnm - anno[k] = lab - pass - - c = dict(matched=matched, all=np.arange(len(names))) - c['unmatched'] = np.where(np.isin(c['all'], c['matched'], invert=True ))[0] - if not unmatched_k is None: - d[unmatched_k] = c['unmatched'] - anno[unmatched_k] = "OTHER" - pass - return d, anno - - - @classmethod - def Find_name_indices_re_match(cls, names, _ptn): - """ - :param _ptn: name regexp string - :return qnm,nn,label: - - qnm - array of names array indices that match the pattern - nn - array of names matching the pattern - label - match key when the pattern string includes one eg (?