eic-opticks/ana/feature.py

0001 #!/usr/bin/env python
0002
0003 import os, logging, numpy as np
0004 log = logging.getLogger(__name__)
0005
0006 PV = not "NOPV" in os.environ
0007 if PV:
0008     try:
0009         from pyvista.plotting.colors import hexcolors
0010     except ImportError:
0011         hexcolors = None
0012     pass
0013 else:
0014     hexcolors = None
0015 pass
0016
0017 def shorten_surfname(name=None, cut=20):
0018     if name is None: name="HamamatsuR12860_PMT_20inch_photocathode_logsurf2"
0019     name = name[:cut//2] + "..." + name[-cut//2:]
0020     return name
0021
0022 def shorten_bndname(name=None, cut=20):
0023     if name is None: name = "Pyrex/HamamatsuR12860_PMT_20inch_photocathode_logsurf2/HamamatsuR12860_PMT_20inch_photocathode_logsurf1/Vacuum"
0024     elem = name.split("/")
0025     if len(elem) == 4:
0026         if len(elem[1]) > cut: elem[1] = shorten_surfname(elem[1], cut)
0027         if len(elem[2]) > cut: elem[2] = shorten_surfname(elem[2], cut)
0028         ret = "/".join(elem)
0029     else:
0030         ret = name
0031     pass
0032     return ret
0033
0034
0035 def make_colors():
0036     """
0037     :return colors: large list of color names with easily recognisable ones first
0038     """
0039     #colors = ["red","green","blue","cyan","magenta","yellow","pink","purple"]
0040     all_colors = list(hexcolors.keys()) if not hexcolors is None else []
0041     easy_colors = "red green blue cyan magenta yellow pink".split()
0042     skip_colors = "bisque beige white aliceblue antiquewhite aqua azure black".split()  # skip colors that look too alike or difficult to see against white or black
0043     skip_colors += "raw_sienna blanchedalmond".split()   # perhaps not available in some matplotlib version or problem from underscore
0044
0045     colors = easy_colors
0046     for c in all_colors:
0047         if c in skip_colors:
0048             continue
0049         if not c in colors:
0050             colors.append(c)
0051         pass
0052     pass
0053     return colors
0054
0055 COLORS = make_colors()
0056
0057
0058 class Feature(object):
0059     """
0060     Trying to generalize feature handling
0061     """
0062     def __init__(self, name, val, vname={}, symbol="ft" ):
0063         """
0064         :param name: string eg "bnd" or "primIdx"
0065         :param val: large array of integer feature values
0066         :param namedict: dict relating feature integers to string names
0067
0068         Canonical usage is from the below SimtraceFeatures
0069
0070
0071         The is an implicit assumption that the number of unique feature values is not enormous,
0072         for example boundary values or prim identity values.
0073
0074         If this comes up with feature names of form "ERR-3105" "ERR-32" etc.. then
0075         probably the CF geometry loaded into python does not match the geometry
0076         used by the C++ run that created the SEvt arrays.
0077
0078         In which case you need to arrange for the CF geometry to be saved alongside
0079         the event arrays and transferred alongside them when grabbed and also
0080         ensure that the CFBASE envvar that controls which geometry is
0081         loaded by CSG/CSGFoundry.py matches the geometry used to create the SEvt arrays.
0082
0083         TODO: automated way to check consistency via metadata matching between geometry and evt
0084         """
0085         uval, ucount = np.unique(val, return_counts=True)
0086
0087         if len(vname) == 0:
0088             nn = ["%s%d" % (name,i) for i in uval]
0089             vname = dict(zip(uval,nn))
0090         pass
0091         pass
0092         idxdesc = np.argsort(ucount)[::-1]
0093         # indices of uval and ucount that reorder those arrays into descending count order
0094
0095         ocount = [ucount[j]       for j in idxdesc]
0096         ouval  = [uval[j]         for j in idxdesc]
0097
0098         # vname needs absolutes to get the names
0099         onames = [vname.get(uval[j], "ERR-%d" % uval[j])  for j in idxdesc]
0100
0101         self.is_bnd = name == "bnd"
0102         self.name = name
0103         self.val = val
0104         self.vname = vname
0105
0106         self.uval = uval
0107         self.unum = len(uval)
0108         self.ucount = ucount
0109         self.idxdesc = idxdesc
0110         self.onames = onames
0111         self.ocount = ocount
0112         self.ouval = ouval
0113
0114         ISEL = os.environ.get("ISEL","")
0115         isel = self.parse_ISEL(ISEL, onames)
0116         sisel = ",".join(map(str, isel))
0117
0118         print( "Feature name %s ISEL: [%s] isel: [%s] sisel [%s] " % (name, ISEL, str(isel), sisel))
0119
0120         self.isel = isel
0121         self.sisel = sisel
0122
0123     @classmethod
0124     def parse_ISEL(cls, ISEL, onames):
0125         """
0126         :param ISEL: comma delimited list of strings or integers
0127         :param onames: names ordered in descending frequency order
0128         :return isels: list of frequency order indices
0129
0130         Integers in the ISEL are interpreted as frequency order indices.
0131
0132         Strings are interpreted as fragments to look for in the ordered names,
0133         (which could be boundary names or prim names for example)
0134         eg use Hama or NNVT to yield the list of frequency order indices
0135         with corresponding names containing those strings.
0136         """
0137         ISELS = list(filter(None,ISEL.split(",")))
0138         isels = []
0139         for i in ISELS:
0140             if i.isnumeric():
0141                 isels.append(int(i))
0142             else:
0143                 for idesc, nam in enumerate(onames):
0144                     if i in nam:
0145                         isels.append(idesc)
0146                     pass
0147                 pass
0148             pass
0149         pass
0150         return isels
0151
0152     def __call__(self, idesc):
0153         """
0154         :param idesc: zero based index less than unum
0155
0156         for frame photons, empty pixels give zero : so not including 0 in ISEL allows to skip
0157         if uval==0 and not 0 in isel: continue
0158
0159         """
0160         assert idesc > -1 and idesc < self.unum
0161         fname = self.onames[idesc]
0162         if self.is_bnd: fname = shorten_bndname(fname)
0163
0164         uval = self.ouval[idesc]
0165         count = self.ocount[idesc]
0166         isel = self.isel
0167
0168         if fname[0] == "_":
0169             fname = fname[1:]
0170         pass
0171         color = COLORS[idesc % len(COLORS)]  # gives the more frequent boundary the easy_color names
0172         msg = " %2d : %5d : %6d : %20s : %80s " % (idesc, uval, count, color, fname  )
0173         selector = self.val == uval
0174
0175         if len(isel) == 0:
0176             skip = False
0177         else:
0178             skip = idesc not in isel
0179         pass
0180         return uval, selector, fname, color, skip, msg
0181
0182     def __str__(self):
0183         lines = []
0184         lines.append(self.desc)
0185         for idesc in range(self.unum):
0186             uval, selector, fname, color, skip, msg = self(idesc)
0187             lines.append(msg)
0188         pass
0189         return "\n".join(lines)
0190
0191     desc = property(lambda self:"ph.%sfeat : %s " % (self.name, str(self.val.shape)))
0192
0193     def __repr__(self):
0194         return "\n".join([
0195             "Feature name %s val %s" % (self.name, str(self.val.shape)),
0196             "uval %s " % str(self.uval),
0197             "ucount %s " % str(self.ucount),
0198             "idxdesc %s " % str(self.idxdesc),
0199             "onames %s " % " ".join(self.onames),
0200             "ocount %s " % str(self.ocount),
0201             "ouval %s " % " ".join(map(str,self.ouval)),
0202             ])
0203
0204
0205 class SimtraceFeatures(object):
0206     """
0207     feat contriols how to select positions, eg  via boundary or identity
0208     allow plotting of subsets with different colors
0209     """
0210     @classmethod
0211     def SubMock(cls, i, num):
0212         p = np.zeros([num, 4, 4], dtype=np.float32)
0213         offset = i*100
0214         for j in range(10):
0215             for k in range(10):
0216                 idx = j*10+k
0217                 if idx < num:
0218                     p[idx,0,0] = float(offset+j*10)
0219                     p[idx,0,1] = 0
0220                     p[idx,0,2] = float(offset+k*10)
0221                     p.view(np.int32)[idx,3,3] = i << 16
0222                 pass
0223             pass
0224         pass
0225         return p
0226
0227     @classmethod
0228     def Mock(cls):
0229         """
0230         Random number of items between 50 and 100 for each of 10 categories
0231         """
0232         aa = []
0233         for i in range(10):
0234             aa.append(cls.SubMock(i, np.random.randint(0,100)))
0235         pass
0236         return np.concatenate(tuple(aa))
0237
0238     def __init__(self, pos, cf=None, featname="pid", do_mok=False, symbol="pf" ):
0239         """
0240         :param pos: SimtracePositions instance
0241         :param cf: CSGFoundry instance
0242
0243         Although at first sight it looks like could use photons array
0244         argument rather than Positions instance, that is not the case when masks are applied.
0245         The pos.p is changed by the mask as is needed such that feature
0246         selectors can be used within the masked arrays.
0247
0248         HMM: identity access is only fully applicable to simtrace, not photons
0249
0250         * ACTUALLY THE SAME INFO IS PRESENT IN PHOTON ARRAYS BUT IN DIFFERENT POSITIONS
0251         * TODO: accomodate the photon layout as well as the simtrace one by using
0252           some common method names with different imps for SimtracePositions and PhotonPositions
0253
0254           * OR: standardize the flag/identity layout between photons and simtrace ?
0255
0256
0257         bnd = p[:,2,3].view(np.int32)
0258         ids = p[:,3,3].view(np.int32)
0259         pid = ids >> 16          # prim_idx
0260         ins = ids & 0xffff       # instance_id
0261
0262              +  +  +  +
0263              +  +  +  +
0264              +  +  + bnd
0265              +  +  + ids
0266
0267         qudarap/qevent.h::
0268
0269             232 QEVENT_METHOD void qevent::add_simtrace( unsigned idx, const quad4& p, const quad2* prd, float tmin )
0270             233 {
0271             234     float t = prd->distance() ;
0272             235     quad4 a ;
0273             236
0274             237     a.q0.f  = prd->q0.f ;
0275             238
0276             239     a.q1.f.x = p.q0.f.x + t*p.q1.f.x ;
0277             240     a.q1.f.y = p.q0.f.y + t*p.q1.f.y ;
0278             241     a.q1.f.z = p.q0.f.z + t*p.q1.f.z ;
0279             242     a.q1.i.w = 0.f ;
0280             243
0281             244     a.q2.f.x = p.q0.f.x ;
0282             245     a.q2.f.y = p.q0.f.y ;
0283             246     a.q2.f.z = p.q0.f.z ;
0284             247     a.q2.u.w = prd->boundary() ; // used by ana/feature.py from CSGOpitXSimtraceTest,py
0285             248
0286             249     a.q3.f.x = p.q1.f.x ;
0287             250     a.q3.f.y = p.q1.f.y ;
0288             251     a.q3.f.z = p.q1.f.z ;
0289             252     a.q3.u.w = prd->identity() ;  // identity from __closesthit__ch (( prim_idx & 0xffff ) << 16 ) | ( instance_id & 0xffff )
0290             253
0291             254     simtrace[idx] = a ;
0292             255 }
0293
0294         TODO: try to pass along the gas_idx in the prd ?
0295         HMM: but maybe that would not distinguish between HighQE and ordinary probably ?
0296         suggests will need to do instance id lookup
0297
0298         """
0299         p = pos.simtrace  ## CAUTION : changed simtrace layout might not be accomodated yet
0300
0301         log.info("[Photons p.ndim %d p.shape %s " % (int(p.ndim), str(p.shape)) )
0302         assert featname in ["pid", "bnd", "ins", "mok"]
0303         if p.ndim == 3:  # eg (121000, 4, 4)
0304             bnd = p[:,2,3].view(np.int32)
0305             ids = p[:,3,3].view(np.int32)   ## see sevent::add_simtrace a.q3.u.w = prd->identity()
0306         elif p.ndim == 4:
0307             bnd = p.view(np.int32)[:,:,2,3]
0308             ids = p.view(np.int32)[:,:,3,3]
0309         else:
0310             log.info("unexpected p.shape %s " % str(p.shape))
0311         pass
0312         pid = ids >> 16       ## HMM: what about MISS ?
0313         ins = ids & 0xffff
0314
0315         log.debug("[ Photons.bndfeat ")
0316         bnd_namedict = {} if cf is None or cf.sim is None else cf.sim.bndnamedict
0317         bnd_namedict[0xffff] = "0xffff"
0318
0319
0320         bndfeat = Feature("bnd", bnd, bnd_namedict)
0321         log.debug("] Photons.bndfeat ")
0322
0323         log.debug("[ Photons.pidfeat ")
0324         pid_namedict = {} if cf is None else cf.primIdx_meshname_dict
0325         pid_namedict[-1] = "-1"
0326
0327
0328         log.info(" pid_namedict: %d  " % len(pid_namedict))
0329         pidfeat = Feature("pid", pid, pid_namedict)
0330         log.debug("] Photons.pidfeat ")
0331
0332         log.debug("[ Photons.insfeat ")
0333         ins_namedict = {} if cf is None else cf.insnamedict
0334         log.info(" ins_namedict: %d  " % len(ins_namedict))
0335         insfeat = Feature("ins", ins, ins_namedict)
0336         log.debug("] Photons.insfeat ")
0337
0338         if do_mok:
0339             log.info("[ Photons.mokfeat ")
0340             mok_namedict = {} if cf is None else cf.moknamedict
0341             mokfeat = Feature("mok", pid, mok_namedict)
0342             log.info("] Photons.mokfeat ")
0343         else:
0344             mokfeat = None
0345         pass
0346
0347         if featname=="pid":
0348             feat = pidfeat
0349         elif featname == "bnd":
0350             feat = bndfeat
0351         elif featname == "ins":
0352             feat = insfeat
0353         elif featname == "mok":
0354             feat = mokfeat
0355         else:
0356             feat = None
0357         pass
0358
0359         self.cf = cf
0360         self.p = p
0361         self.bnd = bnd
0362         self.ids = ids
0363         self.bndfeat = bndfeat
0364         self.pidfeat = pidfeat
0365         self.insfeat = insfeat
0366         self.mokfeat = mokfeat
0367         self.feat = feat
0368         log.info("]Photons")
0369
0370         print(bndfeat)
0371         print(pidfeat)
0372         print(insfeat)
0373
0374
0375     def __repr__(self):
0376         return "\n".join([
0377                "p %s" % str(self.p.shape),
0378                ])
0379
0380
0381 def test_mok(cf):
0382     mock_photons = PhotonFeatures.Mock()
0383     ph = PhotonFeatures(mock_photons, cf, featname="mok", do_mok=True)
0384     print(ph.mokfeat)
0385
0386
0387 if __name__ == '__main__':
0388     from opticks.CSG.CSGFoundry import CSGFoundry
0389     cf = CSGFoundry.Load()
0390     test_mok(cf)
0391