eic-opticks/ana/abstat.py

0001 #!/usr/bin/env python
0002 #
0003 # Copyright (c) 2019 Opticks Team. All Rights Reserved.
0004 #
0005 # This file is part of Opticks
0006 # (see https://bitbucket.org/simoncblyth/opticks).
0007 #
0008 # Licensed under the Apache License, Version 2.0 (the "License");
0009 # you may not use this file except in compliance with the License.
0010 # You may obtain a copy of the License at
0011 #
0012 #   http://www.apache.org/licenses/LICENSE-2.0
0013 #
0014 # Unless required by applicable law or agreed to in writing, software
0015 # distributed under the License is distributed on an "AS IS" BASIS,
0016 # WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
0017 # See the License for the specific language governing permissions and
0018 # limitations under the License.
0019 #
0020
0021 """
0022
0023 ABStat slicing::
0024
0025     In [22]: st[10:20]
0026     Out[22]:
0027     ABStat 10 slice(10, 20, None) na,nb,qctx,reclab,X,Y,Z,T,A,B,C,R
0028     ===== ===== =============================== ====================== ===== ===== ===== ===== ===== ===== ===== =====
0029     na    nb    qctx                            reclab                 X     Y     Z     T     A     B     C     R
0030     ===== ===== =============================== ====================== ===== ===== ===== ===== ===== ===== ===== =====
0031     45490 45054 TO_SC_BT_BT_BT_BT_SA/2/XYZTABCR TO SC [BT] BT BT BT SA  1.12  1.29  1.12  0.86  0.70  1.03  1.00  0.89
0032     45490 45054 TO_SC_BT_BT_BT_BT_SA/3/XYZTABCR TO SC BT [BT] BT BT SA  1.03  1.39  1.17  0.87  0.71  0.98  1.01  0.76
0033     45490 45054 TO_SC_BT_BT_BT_BT_SA/4/XYZTABCR TO SC BT BT [BT] BT SA  1.04  0.93  0.85  1.63  0.70  0.99  1.01  1.11
0034     45490 45054 TO_SC_BT_BT_BT_BT_SA/5/XYZTABCR TO SC BT BT BT [BT] SA  0.98  0.91  0.82  1.62  1.02  1.01  0.99  1.21
0035     45490 45054 TO_SC_BT_BT_BT_BT_SA/6/XYZTABCR TO SC BT BT BT BT [SA]  0.80  0.99  1.01  4.78  1.02  1.01  0.99  1.18
0036     28955 28649 TO_BT_BT_BT_BT_AB/0/XYZTABCR    [TO] BT BT BT BT AB     1.63  1.63  1.63  1.63  1.63  1.63  1.63  1.63
0037     28955 28649 TO_BT_BT_BT_BT_AB/1/XYZTABCR    TO [BT] BT BT BT AB     1.63  1.63  1.63  1.63  1.63  1.63  1.63  1.63
0038     28955 28649 TO_BT_BT_BT_BT_AB/2/XYZTABCR    TO BT [BT] BT BT AB     1.63  1.63  1.63  1.63  1.63  1.63  1.63  1.63
0039     28955 28649 TO_BT_BT_BT_BT_AB/3/XYZTABCR    TO BT BT [BT] BT AB     1.63  1.63  1.63  1.63  1.63  1.63  1.63  1.63
0040     28955 28649 TO_BT_BT_BT_BT_AB/4/XYZTABCR    TO BT BT BT [BT] AB     1.63  1.63  1.63  1.63  1.63  1.63  1.63  1.63
0041     ===== ===== =============================== ====================== ===== ===== ===== ===== ===== ===== ===== =====
0042
0043     ABStat 10 slice(10, 20, None) na,nb,qctx,reclab,X,Y,Z,T,A,B,C,R
0044
0045
0046 Lookup a line via reclab::
0047
0048     In [10]: st[st.st.reclab == "TO BT BT BT BT DR BT BT BT BT BT BT BT BT [SA]"]
0049     Out[10]:
0050     ABStat 1 iv,is,na,nb,qctx,reclab,X,Y,Z,T,A,B,C,R
0051     == == ==== ==== ======================================================= ============================================== ===== ===== ===== ===== ===== ===== ===== =====
0052     iv is na   nb   qctx                                                    reclab                                         X     Y     Z     T     A     B     C     R
0053     == == ==== ==== ======================================================= ============================================== ===== ===== ===== ===== ===== ===== ===== =====
0054     78 11 5339 5367 TO_BT_BT_BT_BT_DR_BT_BT_BT_BT_BT_BT_BT_BT_SA/e/XYZTABCR TO BT BT BT BT DR BT BT BT BT BT BT BT BT [SA]  1.05  1.24  1.00 79.14  0.96  1.28  1.14  1.00
0055     == == ==== ==== ======================================================= ============================================== ===== ===== ===== ===== ===== ===== ===== =====
0056
0057
0058 """
0059 import os, logging, numpy as np
0060 import numpy.lib.recfunctions as rf
0061
0062 from opticks.ana.make_rst_table import recarray_as_rst
0063
0064 log = logging.getLogger(__name__)
0065
0066
0067 class ABStat(object):
0068     """
0069     Hmm stats should probably have a standard path within the event tree
0070     """
0071     STATPATH = "$TMP/stat.npy"
0072     SUPTITLE = "  %(det)s/%(src)s/%(tag)s  %(iv)s/%(is)s %(na)d/%(nb)d   %(reclab)-50s  XYZT: %(X)4.2f %(Y)4.2f %(Z)4.2f %(T)4.2f ABCW: %(A)4.2f %(B)4.2f %(C)4.2f %(W)4.2f  seqc2 %(seqc2)4.2f dstc2 %(distc2)4.2f "
0073     SKIP = "qctx".split()
0074
0075     @classmethod
0076     def path_(cls):
0077         return os.path.expandvars(cls.STATPATH)
0078
0079     def __init__(self, ok, st):
0080         self.ok = ok
0081         self.st = st
0082         self.ar = self.ary
0083         self.sli = slice(0,None,1)
0084
0085     def save(self):
0086         np.save(self.path_(),self.st)
0087
0088     @classmethod
0089     def load(cls, ok):
0090         ra = np.load(cls.path_()).view(np.recarray)
0091         return cls(ok, ra)
0092
0093     def __repr__(self):
0094         return "\n".join([self.brief,recarray_as_rst(self.st[self.sli], skip=self.SKIP),self.brief])
0095
0096     def _get_brief(self):
0097         return "ABStat %s %s " % (len(self.st[self.sli]), ",".join(self.names) )
0098     brief = property(_get_brief)
0099
0100     def _get_names(self):
0101         names = map( lambda k:None if k in self.SKIP else k, self.st.dtype.names )
0102         return filter(None, names)
0103     names = property(_get_names)
0104
0105     def __getitem__(self, sli):
0106         if type(sli) is int:
0107             sli = slice(sli,sli+1)
0108         pass
0109         self.sli = sli
0110         return self
0111
0112     def _get_ary(self, qwns=None):
0113         if qwns is None:
0114             qwns = list(self.ok.qwns)
0115         pass
0116         return np.vstack([self.st[q] for q in qwns]).T
0117     ary = property(_get_ary)
0118
0119     def _get_suptitle(self):
0120         """
0121         After making a single line selection this provides a plot title::
0122
0123             In [11]: st[26].suptitle
0124             Out[11]: ' 26/5 20238/20140   TO [RE] BT BT BT BT SA                              XYZT: 0.85 0.00 0.00 1.31 ABCW: 1.12 1.37 1.10 0.78  seqc2 0.24 distc2 1.10 '
0125
0126             In [16]: st[st.st.reclab=="TO RE BT [BT] BT BT DR BT BT BT BT SA"]
0127             Out[16]:
0128             ABStat 1 iv,is,na,nb,reclab,X,Y,Z,T,A,B,C,W,seqc2,distc2
0129             === == === === ===================================== ===== ===== ===== ===== ===== ===== ===== ===== ===== ======
0130             iv  is na  nb  reclab                                X     Y     Z     T     A     B     C     W     seqc2 distc2
0131             === == === === ===================================== ===== ===== ===== ===== ===== ===== ===== ===== ===== ======
0132             868 97 125 126 TO RE BT [BT] BT BT DR BT BT BT BT SA  0.00  0.00  0.00  1.84  0.00  0.00  0.00  1.24  0.00  1.08
0133             === == === === ===================================== ===== ===== ===== ===== ===== ===== ===== ===== ===== ======
0134
0135             ABStat 1 iv,is,na,nb,reclab,X,Y,Z,T,A,B,C,W,seqc2,distc2
0136
0137             In [17]: st.suptitle
0138             Out[17]: ' 868/97 125/126   TO RE BT [BT] BT BT DR BT BT BT BT SA               XYZT: 0.00 0.00 0.00 1.84 ABCW: 0.00 0.00 0.00 1.24  seqc2 0.00 distc2 1.08 '
0139
0140         """
0141         nli = len(self.st[self.sli])
0142         if nli == 1:
0143             st_tup = self.st[self.sli][0]
0144
0145             d = {}
0146             d.update(self.ok.ctx)
0147             d.update(dict(zip(self.st.dtype.names,st_tup)))
0148             return self.SUPTITLE % d
0149         pass
0150         log.warning("st.suptitle requires single line slice")
0151         return None
0152     suptitle = property(_get_suptitle)
0153
0154
0155     def chi2sel(self, chi2cut=None,  statcut=None, style=None):
0156         """
0157         :return indices: propagation record point indices with chi2 sum greater than cut
0158         """
0159         if chi2cut is None:
0160             chi2cut = self.ok.chi2selcut
0161         if statcut is None:
0162             statcut = self.ok.statcut
0163         if style is None:
0164             style = "distc2stat"
0165         pass
0166
0167         if style == "distc2":
0168             s = np.where( self.st.distc2 > chi2cut)
0169         elif style == "distc2stat":
0170             s = np.where( np.logical_and(self.st.distc2 > chi2cut, self.st.na > statcut))
0171         elif style == "seqc2":
0172             s = np.where( self.st.seqc2 > chi2cut)
0173         elif style == "qwnsum":
0174             s = np.where( np.sum(self.ar, axis=1) > chi2cut )
0175         else:
0176             assert 0, style
0177         pass
0178
0179         log.info("style %s chi2cut %s statcut %s nsel %d " % (style, chi2cut, statcut, len(s[0])))
0180
0181         return s
0182
0183     def qctxsel(self, chi2cut=None, statcut=None, style="qwnsum"):
0184         """
0185         :return qctx list: were chi2 sum exceeds the cut
0186         """
0187         return self.st.qctx[self.chi2sel(chi2cut=chi2cut, statcut=statcut, style=style)]
0188
0189     def reclabsel(self, chi2cut=None, statcut=None, style="distc2stat"):
0190         """
0191         :return reclab list: were chi2 sum exceeds the cut
0192         """
0193         return self.st.reclab[self.chi2sel(chi2cut=chi2cut, statcut=statcut, style=style)]
0194
0195
0196     @classmethod
0197     def dump(cls, st=None):
0198         if st is None:
0199             st = cls.load()
0200         pass
0201         print(st)
0202
0203
0204 if __name__ == '__main__':
0205     from opticks.ana.main import opticks_main
0206     ok = opticks_main()
0207     st = ABStat.load(ok)
0208
0209     #print(st)
0210     print(st[st.chi2sel()])
0211
0212     ar = st.ar
0213
0214