eic-opticks/ana/p.py

0001 #!/usr/bin/env python
0002 """
0003 p.py : HMM perhaps rename to sysrap/sphoton.py beside sphoton.h
0004 ===========================================================================
0005
0006 NB locations and packing used here must match sysrap/sphoton.h::
0007
0008     050 struct sphoton
0009      51 {
0010      52     float3 pos ;
0011      53     float  time ;
0012      54
0013      55     float3 mom ;
0014      56     float  weight ;
0015      57
0016      58     float3 pol ;
0017      59     float  wavelength ;
0018      60
0019      61     unsigned boundary_flag ;  // p.view(np.uint32)[3,0]
0020      62     unsigned identity ;       // p.view(np.uint32)[3,1]
0021      63     unsigned orient_idx ;     // p.view(np.uint32)[3,2]
0022      64     unsigned flagmask ;       // p.view(np.uint32)[3,3]
0023      65
0024      66
0025      67     SPHOTON_METHOD void set_prd( unsigned  boundary, unsigned  identity, float  orient );
0026
0027
0028      72     SPHOTON_METHOD void set_orient(float orient){ orient_idx = ( orient_idx & 0x7fffffffu ) | (( orient < 0.f ? 0x1 : 0x0 ) << 31 ) ; } // clear orient bit and then set it
0029      73     SPHOTON_METHOD void set_idx( unsigned idx ){  orient_idx = ( orient_idx & 0x80000000u ) | ( 0x7fffffffu & idx ) ; }   // retain bit 31 asis
0030      74
0031     105 SPHOTON_METHOD void sphoton::set_prd( unsigned  boundary_, unsigned  identity_, float  orient_ )
0032     106 {
0033     107     set_boundary(boundary_);
0034     108     identity = identity_ ;
0035     109     set_orient( orient_ );
0036     110 }
0037
0038
0039
0040 """
0041 import numpy as np
0042 import hashlib, builtins
0043 from opticks.ana.hismask import HisMask
0044 from opticks.ana.histype import HisType
0045 from opticks.ana.nibble import count_nibbles
0046 from opticks.ana.nbase import cus
0047
0048
0049 hm = HisMask()
0050 ht = HisType()
0051
0052 seqhis_ = lambda s:ht.label(s)
0053 cseqhis_ = lambda l:ht.code(l)
0054 seqnib_ = lambda s:count_nibbles(s)
0055
0056 #seqdesc_ = lambda s:"seqdesc_ %20s " %
0057
0058
0059
0060 class CUSS(np.ndarray): pass
0061
0062 def cuss(s,w=None):
0063     """
0064     CountUniqueSubSelections
0065
0066     :param s: seqhis array with same shape as w that provides the histories of the selection
0067     :param w: None OR array of indices expressing a selection of eg deviant photons created for example using np.unique np.where
0068     :return o: CUSS instance view of cus count-unique-sorted array listing unique seqhis
0069                and their counts in descending frequency order
0070
0071     This aims to make it easier to investigate seqhis sub-selections of selections (eg deviants)
0072     CAUTION: using this auto-populates via builtins the calling scope with np.where arrays
0073     for each of the sub-selections::
0074
0075           w0, w1, w2, ...
0076
0077     Usage example::
0078
0079         w = np.unique(np.where( np.abs(a.photon - b.photon) > 0.1 )[0])  ## unique indices selection
0080         s = a.seq[w,0]     # seqhis histories of the selection
0081         o = cuss(s,w)
0082         print(o)
0083         print(w1)
0084
0085     """
0086     if w is None:
0087         w = np.arange(len(s))
0088     pass
0089     assert w.shape == s.shape
0090     cu = cus(s)   ## count unique sorted
0091     for i in range(len(cu)): setattr( builtins, "w%d"%i, w[s == cu[i,0]] )
0092
0093     o = np.zeros( (cu.shape[0], cu.shape[1]+2), dtype=np.object )
0094     o[:,0] = list(map(lambda _:"w%d" % _ , list(range(len(cu))) ))
0095     o[:,1] = list(map(lambda _:"%30s" % _, seqhis_(cu[:,0]) ))
0096     o[:,2] = list(map(lambda _:"%16s" % _, cu[:,0] ))
0097     o[:,3] = list(map(lambda _:"%16s" % _, cu[:,1] ))
0098
0099     ret = o.view(CUSS)
0100     ret.cu = cu
0101     return ret
0102
0103
0104
0105
0106 digest_ = lambda a:hashlib.md5(a.data.tobytes()).hexdigest()
0107
0108 ## using ellipsis avoids having to duplicate for photons and records
0109 #3,0
0110 boundary___ = lambda p:p.view(np.uint32)[...,3,0] >> 16
0111 boundary__ = lambda p:p.view(np.uint32)[:,3,0] >> 16
0112 boundary_  = lambda p:p.view(np.uint32)[3,0] >> 16
0113
0114 flag___    = lambda p:p.view(np.uint32)[...,3,0] & 0xffff
0115 flag__    = lambda p:p.view(np.uint32)[:,3,0] & 0xffff
0116 flag_     = lambda p:p.view(np.uint32)[3,0] & 0xffff
0117
0118 #3,1
0119 identity___   = lambda p:p.view(np.uint32)[...,3,1]
0120 identity__    = lambda p:p.view(np.uint32)[:,3,1]
0121 identity_     = lambda p:p.view(np.uint32)[3,1]
0122
0123 primIdx___    = lambda p:p.view(np.uint32)[...,3,1] >> 16
0124 primIdx__     = lambda p:p.view(np.uint32)[:,3,1] >> 16
0125
0126
0127 isectPrimIdx_ = lambda isect:isect.view(np.uint32)
0128
0129
0130
0131 instanceId___ = lambda p:p.view(np.uint32)[...,3,1] & 0xffff
0132 instanceId__  = lambda p:p.view(np.uint32)[:,3,1] & 0xffff
0133
0134 primIdx_      = lambda p:identity_(p) >> 16
0135 instanceId_   = lambda p:identity_(p) & 0xffff
0136
0137
0138
0139 ident_ = lambda p:p.view(np.uint32)[...,3,1]
0140 prim_  = lambda p:ident_(p) >> 16
0141 inst_  = lambda p:ident_(p) & 0xffff
0142
0143 iindex_ = lambda p:p.view(np.uint32)[...,1,3]
0144
0145 #3,2
0146 idx_      = lambda p:p.view(np.uint32)[3,2] & 0x7fffffff
0147 orient___  = lambda p:p.view(np.uint32)[...,3,2] >> 31
0148 orient_   = lambda p:p.view(np.uint32)[3,2] >> 31
0149
0150 #3,3
0151 flagmask_ = lambda p:p.view(np.uint32)[3,3]
0152
0153 flagdesc_ = lambda p:" idx(%6d) prd(b%3d p%4d i%5d o%1d ii:%5d) %3s  %15s " % ( idx_(p),  boundary_(p),primIdx_(p),instanceId_(p), orient_(p), iindex_(p), hm.label(flag_(p)),hm.label( flagmask_(p) ))
0154
0155
0156 flagmask__ = lambda p:p.view(np.uint32)[:,3,3]
0157 hit__      = lambda p,msk:p[np.where( ( flagmask__(p) & msk ) == msk)]
0158
0159
0160
0161 ridiff_ = lambda ri:ri[1:,:3] - ri[:-1,:3]
0162
0163 # rdist_(a,i)/rtime_(a,i)/rspeed_(a,i) : distance/time/speed between record point i and i+1
0164 rdist_ = lambda a,i:np.sqrt(np.sum( (a.record[:,i+1,0,:3]-a.record[:,i,0,:3])*(a.record[:,i+1,0,:3]-a.record[:,i,0,:3]) , axis=1 ))
0165 rtime_ = lambda a,i:a.record[:,i+1,0,3] - a.record[:,i,0,3]
0166 rspeed_ = lambda a,i:rdist_(a,i)/rtime_(a,i)
0167
0168
0169 ## HMM r1time_ works when i is an array
0170 # r1dist_ when i is array gives single scalar does not
0171 r1time_ = lambda r,i:r[i+1,0,3] - r[i,0,3]
0172 r1dist_ = lambda r,i:np.sqrt(np.sum( (r[i+1,0,:3]-r[i,0,:3])*(r[i+1,0,:3]-r[i,0,:3]) ))
0173 r1speed_ = lambda r,i:r1dist_(r,i)/r1time_(r,i)
0174
0175
0176 ## rv: do that in a more vectorized way
0177 rvtime_ = lambda r:np.diff(r[:,0,3])
0178 rvstep_ = lambda r:np.diff(r[:,0,:3],axis=0 )
0179 rvdist_ = lambda r:np.sqrt(np.sum(rvstep_(r)*rvstep_(r),axis=1))
0180 rvspeed_ = lambda r:rvdist_(r)/rvtime_(r)
0181
0182 ## alternate avoiding intermediary lambda funcs
0183 #rvdist_ = lambda r:np.sqrt(np.sum(np.diff(r[:,0,:3],axis=0 )*np.diff(r[:,0,:3],axis=0 ),axis=1))
0184 #rvspeed_ = lambda r:np.sqrt(np.sum(np.diff(r[:,0,:3],axis=0 )*np.diff(r[:,0,:3],axis=0 ),axis=1))/np.diff(r[:,0,3])
0185
0186
0187
0188
0189 ## TO PICK THE GEOMETRY APPROPRIATE TO THE RESULT ARRAYS SET CFBASE envvar TO DIRECTORY CONTAINING CSGFoundry dir
0190 print("[ana/p.py:from opticks.CSG.CSGFoundry import CSGFoundry ")
0191 from opticks.CSG.CSGFoundry import CSGFoundry
0192 print("]ana/p.py:from opticks.CSG.CSGFoundry import CSGFoundry ")
0193
0194
0195 print("[ana/p.py:cf = CSGFoundry.Load()")
0196 cf = CSGFoundry.Load()
0197 print("]ana/p.py:cf = CSGFoundry.Load()")
0198
0199 ## HMM: NASTY LOADING cf HERE IN THE MODULE
0200 ## DONE TO AVOID cf ARG FOR THE BELOW : MAYBE CURRYING CAN AVOID ?
0201
0202 if not cf is None:
0203     cf_bnd_  = lambda p:cf.sim.bndnamedict.get(boundary_(p),"cf_bnd_ERR")
0204     cf_prim_ = lambda p:cf.primIdx_meshname_dict.get(prim_(p),"cf_prim_ERR")
0205 else:
0206     cf_bnd_  = lambda p:"no_cf"
0207     cf_prim_  = lambda p:"no_cf"
0208 pass
0209
0210 bflagdesc_ = lambda p:"%s : %60s : %s : %s " % ( flagdesc_(p), cf_prim_(p) , digest_(p[:3])[:8], cf_bnd_(p) )
0211
0212
0213