eic-opticks/sysrap/dv.py

0001 #!/usr/bin/env python
0002
0003 import os, numpy as np
0004 from opticks.ana.fold import Fold
0005 from opticks.ana.array_repr_mixin import ArrayReprMixin
0006
0007 class AlignedDV(ArrayReprMixin, object):
0008     """
0009     Simple deviation comparisons of random aligned arrays
0010
0011     In [2]: pdv.dv
0012     Out[2]:
0013     array([[   47,   117,  1732,  4412,  2710,   965,    16,     1,     0,     0],
0014            [ 2746,  5430,  1724,    96,     4,     0,     0,     0,     0,     0],
0015            [ 6404,  2937,   647,    11,     1,     0,     0,     0,     0,     0],
0016            [ 9995,     1,     1,     3,     0,     0,     0,     0,     0,     0],
0017            [10000,     0,     0,     0,     0,     0,     0,     0,     0,     0]], dtype=uint32)
0018
0019     """
0020     EDGES = np.array( [0.,1e-6,1e-5,1e-4,1e-3, 1e-2, 1e-1, 1, 10, 100, 1000], dtype=np.float32 )
0021     COLUMN_LABEL = ["%g" % _ for _ in EDGES[1:]]
0022
0023     POS = 0
0024     TIME = 1
0025     MOM = 2
0026     POL = 3
0027     WL = 4
0028
0029     ITEMS = (POS, TIME, MOM, POL, WL)
0030     ROW_LABEL = ("pos", "time", "mom", "pol", "wl")
0031
0032     def __init__(self, a, b, arr="photon", symbol=None):
0033         """
0034         :param a: opticks.ana.fold Fold instance
0035         :param b: opticks.ana.fold Fold instance
0036
0037         In [12]: adv[...,0,3].shape
0038         Out[12]: (10000,)
0039
0040         In [13]: bdv[...,0,3].shape
0041         Out[13]: (10000, 10)
0042         """
0043
0044         wseq = np.where( a.seq[:,0] == b.seq[:,0] )
0045         adv = np.abs( getattr(a,arr)[wseq] - getattr(b,arr)[wseq] )  ## for deviations to be meaningful needs to be same history
0046
0047
0048         ## although ... ellipsis indexing can generalize access between photon and record
0049         ## need to use different np.amax axis args to keep the deviation at photon level, not record level
0050         if arr == "photon":
0051             time = adv[:,0,3]
0052             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
0053             mom  = np.amax( adv[:,1,:3], axis=1 )
0054             pol  = np.amax( adv[:,2,:3], axis=1 )
0055             wl   = adv[:,2,3]
0056         elif arr == "record":
0057             time = np.amax( adv[:,:,0,3], axis=1 )
0058             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
0059             mom  = np.amax( adv[:,:,1,:3], axis=(1,2) )
0060             pol  = np.amax( adv[:,:,2,:3], axis=(1,2) )
0061             wl   = np.amax( adv[:,:,2,3], axis=1 )
0062         else:
0063             assert 0, "expecting photon or record not %s " % arr
0064         pass
0065
0066         dv = np.zeros( (len(self.ITEMS),len(self.EDGES)-1), dtype=np.uint32 )
0067         dv[self.POS], bins = np.histogram( pos, bins=self.EDGES )
0068         dv[self.TIME], bins = np.histogram( time, bins=self.EDGES )
0069         dv[self.MOM], bins = np.histogram( mom, bins=self.EDGES )
0070         dv[self.POL], bins = np.histogram( pol, bins=self.EDGES )
0071         dv[self.WL], bins = np.histogram( wl, bins=self.EDGES )
0072
0073         self.symbol = symbol if not symbol is None else "%s.dv" % arr
0074         self.pos = pos
0075         self.time = time
0076         self.mom = mom
0077         self.pol = pol
0078         self.wl = wl
0079         self.dv = dv
0080
0081     def __repr__(self):
0082         return self.MakeRepr(self.dv, symbol=self.symbol)
0083
0084
0085 class HeaderAB(object):
0086     def __init__(self, a, b, cmdline):
0087         self.lines = ["A_FOLD : %s " % a.base, "B_FOLD : %s " % b.base, cmdline ]
0088     def __repr__(self):
0089         return "\n".join(self.lines)
0090
0091
0092 if __name__ == '__main__':
0093     a = Fold.Load("$A_FOLD", symbol="a") if "A_FOLD" in os.environ else None
0094     b = Fold.Load("$B_FOLD", symbol="b") if "B_FOLD" in os.environ else None
0095
0096     hdr = HeaderAB(a,b, "./dv.sh   # cd ~/opticks/sysrap")
0097     print(hdr)
0098
0099     pdv = AlignedDV(a,b, arr="photon", symbol="pdv")
0100     print(pdv)
0101
0102     rdv = AlignedDV(a,b, arr="record", symbol="rdv")
0103     print(rdv)
0104
0105
0106