eic-opticks/ana/dv.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 dv.py
0024 ======
0025
0026
0027 Non-aligned deviation checking
0028 ---------------------------------
0029
0030 ::
0031
0032     tp() { tboolean-;PROXYLV=18 tboolean-proxy-ip $* ; }
0033
0034     ab.rpost_dv maxdvmax: 0.02283 maxdv: 0.02283        0  0.02283  skip:SC AB RE
0035       idx        msg :                            sel :    lcu1     lcu2  :       nitem     nelem/    ndisc: fdisc  mx/mn/av        mx/       mn/      avg  eps:eps
0036      0000            :                    TO BT BT SA :    8794     8794  :        7710    123360/       17: 0.000  mx/mn/av   0.02283/        0/1.705e-06  eps:0.0002
0037      0001            :                       TO BR SA :     580      617  :          33       396/        0: 0.000  mx/mn/av         0/        0/        0  eps:0.0002
0038      0002            :                 TO BT BR BT SA :     561      527  :          27       540/        1: 0.002  mx/mn/av   0.02283/        0/4.227e-05  eps:0.0002
0039     ab.rpol_dv maxdvmax:       0 maxdv:       0        0        0  skip:SC AB RE
0040       idx        msg :                            sel :    lcu1     lcu2  :       nitem     nelem/    ndisc: fdisc  mx/mn/av        mx/       mn/      avg  eps:eps
0041      0000            :                    TO BT BT SA :    8794     8794  :        7710     92520/        0: 0.000  mx/mn/av         0/        0/        0  eps:0.0002
0042      0001            :                       TO BR SA :     580      617  :          33       297/        0: 0.000  mx/mn/av         0/        0/        0  eps:0.0002
0043      0002            :                 TO BT BR BT SA :     561      527  :          27       405/        0: 0.000  mx/mn/av         0/        0/        0  eps:0.0002
0044     ab.ox_dv maxdvmax: 0.00238 maxdv:0.0001221        0  0.00238  skip:SC AB RE
0045       idx        msg :                            sel :    lcu1     lcu2  :       nitem     nelem/    ndisc: fdisc  mx/mn/av        mx/       mn/      avg  eps:eps
0046      0000            :                    TO BT BT SA :    8794     8794  :        7710     92520/        0: 0.000  mx/mn/av 0.0001221/        0/3.652e-06  eps:0.0002
0047      0001            :                       TO BR SA :     580      617  :          33       396/        0: 0.000  mx/mn/av         0/        0/        0  eps:0.0002
0048      0002            :                 TO BT BR BT SA :     561      527  :          27       324/       14: 0.043  mx/mn/av   0.00238/        0/3.823e-05  eps:0.0002
0049
0050
0051
0052 Observations:
0053
0054 * nitem much lower than photon counts when you have a BR in history
0055 * "--reflectcheat" is not enabled by default : enabling it can increase the stats of the comparison
0056
0057
0058 Non-random-aligned deviation checking relies on
0059
0060 1. same input photons for both simulations A and B
0061 2. "accidental" history alignment between A and B : helped by BR reflectcheat
0062
0063 Need to find some corroboration:
0064
0065     For histories like "TO BT BT SA" what happens is purely geometry, so
0066     for those photons in A and B that follow this history can directly compare.
0067     But some fraction of BT in simulation A will be BR in the simulation B as
0068     the random numbers are not aligned.  To reduce this the reflectcheat technique
0069     is used to decide which fraction of sample to reflect and which to transmit
0070     based on the ratio of record_id to total
0071
0072
0073 Examining the deviation::
0074
0075     In [15]: ab.rpost_dv
0076     Out[15]:
0077     rpost_dv
0078      0000            :                          TO SA :   55321    55303  :     55249/      0: 0.000  mx/mn/av 0.0000/0.0000/     0
0079      0001            :                    TO BT BT SA :   39222    39231  :     34492/      8: 0.000  mx/mn/av 0.0138/0.0000/3.192e-06
0080      0002            :                       TO BR SA :    2768     2814  :       188/      0: 0.000  mx/mn/av 0.0000/0.0000/     0
0081      0003            :                 TO BT BR BT SA :    2425     2369  :       125/      0: 0.000  mx/mn/av 0.0000/0.0000/     0
0082      0004            :              TO BT BR BR BT SA :     151      142  :         1/      0: 0.000  mx/mn/av 0.0000/0.0000/     0
0083
0084     In [12]: dv = ab.rpost_dv.dvs[1].dv
0085     In [16]: av = ab.rpost_dv.dvs[1].av
0086     In [17]: bv = ab.rpost_dv.dvs[1].bv
0087
0088     In [14]: np.where( dv > 0 )
0089     Out[14]:
0090     (
0091         A([ 8019,  8019,  8019,  8019, 13879, 13879, 13879, 13879]),
0092         A([0, 1, 2, 3, 0, 1, 2, 3]),
0093         A([1, 1, 1, 1, 0, 0, 0, 0])
0094     )
0095
0096     In [18]: wdv = np.where( dv > 0 )
0097
0098     In [19]: av[wdv]
0099     Out[19]:
0100     A([  30.4181,   30.4181,   30.4181,   30.4181,  116.2219,  116.2219,  116.2219,  116.2219])
0101
0102     In [20]: bv[wdv]
0103     Out[20]:
0104     A([  30.4319,   30.4319,   30.4319,   30.4319,  116.2357,  116.2357,  116.2357,  116.2357])
0105
0106     In [21]: av[wdv] - bv[wdv]
0107     Out[21]:
0108     A([-0.0138, -0.0138, -0.0138, -0.0138, -0.0138, -0.0138, -0.0138, -0.0138])
0109
0110
0111 """
0112 import os, sys, logging, numpy as np
0113 from opticks.ana.log import fatal_, error_, warning_, info_, debug_
0114 from opticks.ana.log import underline_, blink_
0115 log = logging.getLogger(__name__)
0116
0117 assert 0, "dont use this : use qdv.py : its much faster "
0118
0119 class Level(object):
0120     FATAL = 20
0121     ERROR = 10
0122     WARNING = 0
0123     INFO = -10
0124     DEBUG = -20
0125
0126     level2name = { FATAL:"FATAL", ERROR:"ERROR", WARNING:"WARNING", INFO:"INFO", DEBUG:"DEBUG" }
0127     name2level = { "FATAL":FATAL, "ERROR":ERROR, "WARNING":WARNING, "INFO":INFO, "DEBUG":DEBUG  }
0128     level2func = { FATAL:fatal_, ERROR:error_, WARNING:warning_, INFO:info_, DEBUG:debug_ }
0129
0130
0131     @classmethod
0132     def FromName(cls, name):
0133         level = cls.name2level[name]
0134         return cls(name, level)
0135     @classmethod
0136     def FromLevel(cls, level):
0137         name = cls.level2name[level]
0138         return cls(name, level)
0139
0140     def __init__(self, name, level):
0141         self.name = name
0142         self.level = level
0143         self.fn_ = self.level2func[level]
0144
0145
0146 class Dv(object):
0147
0148    FMT  =       "  %9d %9d : %5d %5d %5d : %6.4f %6.4f %6.4f : %9.4f %9.4f %9.4f  "
0149    CFMT =       "  %9s %9s : %5s %5s %5s : %6s %6s %6s : %9s %9s %9s    "
0150    CFMT_CUTS =  "  %9s %9s : %5s %5s %5s : %6.4f %6.4f %6.4f : %9s %9s %9s    "
0151    CFMT_COLS =  "nitem nelem nwar nerr nfat fwar ferr ffat mx mn avg".split()
0152
0153    LMT  = " %0.4d %10s : %30s : %7d  %7d "   # labels seqhis line
0154    CLMT = " %4s %10s : %30s : %7s  %7s "
0155
0156    clabel = CLMT % ( "idx", "msg", "sel", "lcu1", "lcu2" )
0157    cblank = CLMT % ( "", "", "", "", "" )
0158
0159
0160    def __init__(self, tab, idx, sel, av, bv, lcu, dvmax, msg=""):
0161        """
0162        :param tab: DvTab instance
0163        :param idx: unskipped orignal seqhis line index
0164        :param sel: single line selection eg 'TO BT BT SA'
0165        :param av: evt a values array within selection
0166        :param bv: evt b values array within selection
0167        :param lcu: list of length 3 with (seqhis-bigint, a-count, b-count)
0168        :param dvmax: triplet of floats for warn/error/fatal deviation levels
0169
0170        Access an Dv instance in ipython::
0171
0172             In [12]: ab.ox_dv.dvs[2]
0173             Out[12]:  0002            :                 TO BT BR BT SA :     561      527  :          27       324/       14: 0.043  mx/mn/av   0.00238/        0/3.823e-05  eps:0.0002
0174
0175        Get at the values::
0176
0177            In [16]: av,bv = ab.ox_dv.dvs[2].av, ab.ox_dv.dvs[2].bv
0178
0179        """
0180        label = self.LMT % ( idx, msg, sel, lcu[1], lcu[2] )
0181        assert len(dvmax) == 3
0182
0183        dv = np.abs( av - bv )
0184        nitem = len(dv)
0185        nelem = dv.size
0186
0187        if nelem>0:
0188
0189            mx = dv.max()
0190            mn = dv.min()
0191            avg = dv.sum()/float(nelem)
0192
0193            disc=[ dv[dv>dvmax[0]], dv[dv>dvmax[1]], dv[dv>dvmax[2]] ]
0194            ndisc = map(len, disc)     # elements, not items
0195            fdisc =  map(lambda _:float(_)/float(nelem), ndisc )
0196        else:
0197            mx = None
0198            mn = None
0199            avg = None
0200            ndisc = None
0201            fdisc = None
0202        pass
0203
0204        self.tab = tab
0205        self.label = label
0206        self.nitem = nitem
0207        self.nelem = nelem
0208        self.ndisc = ndisc
0209        self.fdisc = fdisc
0210        self.mx = mx
0211        self.mn = mn
0212        self.avg = avg
0213        self.ismax = False # set from DvTab
0214
0215        self.av = av
0216        self.bv = bv
0217        self.dv = dv
0218        self.lcu  = lcu
0219        self.dvmax = dvmax
0220        self.msg = msg
0221
0222        if self.mx > self.dvmax[2]:
0223            lev = Level.FromName("FATAL")
0224            lmsg = "  > dvmax[2] %.4f " % self.dvmax[2]
0225        elif self.mx > self.dvmax[1]:
0226            lev = Level.FromName("ERROR")
0227            lmsg = "  > dvmax[1] %.4f " % self.dvmax[1]
0228        elif self.mx > self.dvmax[0]:
0229            lev = Level.FromName("WARNING")
0230            lmsg = "  > dvmax[0] %.4f " % self.dvmax[0]
0231        else:
0232            lev = Level.FromName("INFO")
0233            lmsg = ""
0234        pass
0235        self.fn_ = lev.fn_
0236        self.lev = lev
0237        self.lmsg = lmsg
0238
0239
0240    @classmethod
0241    def columns(cls):
0242        cdesc = cls.CFMT % tuple(cls.CFMT_COLS)
0243        clabel = cls.clabel
0244        return "%s : %s  " % (clabel, cdesc )
0245
0246    def columns2(self):
0247        cdesc2 = self.CFMT_CUTS % ("","","","","", self.dvmax[0], self.dvmax[1], self.dvmax[2], "", "", "")
0248        return "%s : %s  " % (self.cblank, cdesc2 )
0249
0250
0251    def __repr__(self):
0252        if self.nelem>0:
0253            desc =  self.FMT % ( self.nitem, self.nelem, self.ndisc[0], self.ndisc[1], self.ndisc[2], self.fdisc[0], self.fdisc[1], self.fdisc[2], self.mx, self.mn, self.avg )
0254        else:
0255            desc = ""
0256        pass
0257
0258        if self.ismax:
0259            #pdesc = self.fn_(desc)
0260            pdesc = underline_(self.fn_(desc))
0261        else:
0262            pdesc = self.fn_(desc)
0263        pass
0264        return "%s : %s : %s : %s " % (self.label, pdesc, self.fn_("%20s"%self.lev.name), self.fn_(self.lmsg)  )
0265
0266
0267 class DvTab(object):
0268     def is_skip(self, sel):
0269         """
0270         SC AB RE
0271           see notes/issues/sc_ab_re_alignment.rst
0272
0273         When comparing non-aligned events it is necessary to
0274         avoid accidental history alignment causing deviation fails
0275         by skipping any selection that includes SC AB or RE, assuming
0276         reflect cheat is used.
0277         """
0278         sqs = sel.split()
0279         skip = False
0280         for skp in self.skips:
0281              with_sk = skp in sqs
0282              if with_sk:
0283                  skip = True
0284              pass
0285         pass
0286         return skip
0287
0288     def __init__(self, name, seqtab, ab, skips="SC AB RE", selbase=None ):
0289         """
0290         :param name: ox_dv, rpost_dv, rpol_dv
0291         :param seqtab: ab.ahis SeqTable
0292         :param ab:
0293         :param skips:
0294         :param selbase: either None or "ALIGN"
0295
0296         """
0297         log.info("[ sel %s selbase %s  " % (name, selbase) )
0298         self.name = name
0299         self.seqtab = seqtab
0300         self.ab = ab
0301         self.dirty = False
0302         self.skips = skips.split()
0303         #self.sli = slice(None)
0304         self.sli = slice(0,10)
0305         self.selbase = selbase
0306
0307         labels = self.seqtab.labels       # eg list of length 17 : ['TO BT BT SA', 'TO BR SA', ... ]
0308
0309         cu = self.seqtab.cu               # eg with shape (17,3)  the 3 columns being (seqhis, a-count, b-count )
0310         assert len(labels) == len(cu)
0311         nsel = len(labels)
0312
0313         self.labels = labels
0314         self.cu = cu
0315         self.nsel = nsel
0316
0317         ab.aselhis = selbase
0318
0319         self.dvs = self.make_dvs_slowly()
0320
0321         self.findmax()
0322         log.info("] %s " % name )
0323
0324
0325     def make_dvs_slowly(self):
0326         """
0327         Changing selection for every line like this
0328         is too slow for more than 100k photons
0329         """
0330         cu = self.cu
0331         ab = self.ab
0332         dvs = []
0333         for i in range(self.nsel)[self.sli]:
0334             sel = self.labels[i]
0335
0336             if self.is_skip(sel):
0337                 continue
0338             pass
0339
0340             lcu = cu[i]
0341             assert len(lcu) == 3
0342             _, na, nb = lcu
0343
0344             ab.aselhis = sel              # set selection to just this history sequence
0345             assert True == ab.checkrec()  # checking that both evt a and b give same number of record points
0346
0347             dv = self.dv_(i, sel, lcu)    # Dv instance comparing values within single line selection eg all 'TO BT BT SA'
0348
0349             if dv is None:
0350                 log.debug("dv None for i:%d sel:%s lcu:%r " % ( i, sel, lcu))
0351             else:
0352                 dvs.append(dv)
0353             pass
0354
0355             ab.aselhis = self.selbase
0356         pass
0357         return dvs
0358
0359
0360
0361     def _get_dvmax(self):
0362         """
0363         :return dvmax: list of three values corresponding to warn/error/fatal deviation cuts
0364
0365         rpost
0366             extent*2/65535 is the size of shortnorm compression bins
0367             so its the closest can expect values to get.
0368
0369             TODO: currently the time cuts are not correct, a for "t"
0370             should be using different domain
0371
0372         rpol
0373              is extremely compressed, so bins are big : 1*2/255
0374         ox
0375              not compressed stored as floats : other things will limit deviations
0376
0377         """
0378         ab = self.ab
0379         if self.name == "rpost_dv":
0380             eps = ab.fdom[0,0,3]*2.0/((0x1 << 16) - 1)
0381             dvmax = [eps, 1.5*eps, 2.0*eps]
0382         elif self.name == "rpol_dv":
0383             eps = 1.0*2.0/((0x1 << 8) - 1)
0384             dvmax = [eps, 1.5*eps, 2.0*eps]
0385         elif self.name == "ox_dv":
0386             dvmax = ab.ok.pdvmax    ## adhoc input guess at appropriate cuts
0387         else:
0388             assert self.name
0389         pass
0390         return dvmax
0391     dvmax = property(_get_dvmax)
0392
0393
0394     def _get_dvmaxt(self):
0395         """
0396         ta 19::
0397
0398             In [1]: ab.rpost_dv.dvmax
0399             Out[1]: [0.022827496757457846, 0.03424124513618677, 0.04565499351491569]
0400
0401             In [2]: ab.rpost_dv.dvmaxt
0402             Out[2]: [0.0003043666242088853, 0.00045654993631332797, 0.0006087332484177706]
0403
0404         Depends on the rule-of-thumb Opticks::setupTimeDomains but time cuts should be much
0405         more stringent than position ones::
0406
0407             In [3]: ab.fdom[0,0,3]
0408             Out[3]: 748.0
0409
0410             In [4]: ab.fdom[1,0,1]
0411             Out[4]: 9.973333
0412
0413             In [5]: ab.fdom[0,0,3]/ab.fdom[1,0,1]
0414             Out[5]: 75.0
0415
0416         """
0417         ab = self.ab
0418         if self.name == "rpost_dv":
0419             eps_t = ab.fdom[1,0,1]*2.0/((0x1 << 16) - 1)
0420             dvmaxt = [eps_t, 1.5*eps_t, 2.0*eps_t]
0421         else:
0422             assert self.name
0423         pass
0424         return dvmaxt
0425     dvmaxt = property(_get_dvmaxt)
0426
0427
0428     def dv_(self, i, sel, lcu):
0429         """
0430         :param i: unskipped orignal seqhis line index
0431         :param sel: seqhis label eg 'TO BT BT SA'
0432         :param lcu: count unique for single line eg [36045,  8794, 8794]
0433
0434         Accessing av and bv values within the active single line selection ab.aselhis,
0435         and creating a Dv instance from them.
0436         """
0437         log.info("[  %s " % sel )
0438
0439         ab = self.ab
0440         if self.name == "rpost_dv":
0441             av = ab.a.rpost()
0442             bv = ab.b.rpost()
0443         elif self.name == "rpol_dv":
0444             av = ab.a.rpol()
0445             bv = ab.b.rpol()
0446         elif self.name == "ox_dv":
0447             av = ab.a.ox[:,:3,:]
0448             bv = ab.b.ox[:,:3,:]
0449         else:
0450             assert self.name
0451         pass
0452         dvmax = self.dvmax
0453
0454         assert ab.a.sel == ab.b.sel
0455         sel = ab.a.sel
0456         dv = Dv(self, i, sel, av, bv, lcu, dvmax )
0457         ret = dv if len(dv.dv) > 0 else None
0458         log.info("]")
0459         return ret
0460
0461     def _get_float(self, att):
0462         return map(lambda dv:float(getattr(dv, att)), filter(None,self.dvs))
0463
0464     maxdv = property(lambda self:self._get_float("mx"))
0465     def _get_maxdvmax(self):
0466         maxdv_ = self.maxdv
0467         return max(maxdv_) if len(maxdv_) > 0 else -1
0468     maxdvmax   = property(_get_maxdvmax)
0469
0470     def findmax(self):
0471         maxdv = map(lambda _:float(_.mx), self.dvs)
0472         mmaxdv = max(maxdv) if len(maxdv) > 0 else -1
0473         for dv in self.dvs:
0474             if dv.mx == mmaxdv and dv.lev.level > Level.INFO:
0475                 dv.ismax = True
0476             pass
0477         pass
0478
0479     def _get_maxlevel(self):
0480         """
0481         Overall level of the table : INFO, WARNING, ERROR or FATAL
0482         based on the maximum level of the lines
0483         """
0484         levs = map(lambda dv:dv.lev.level, self.dvs)
0485         mxl = max(levs) if len(levs) > 0 else None
0486         return Level.FromLevel(mxl) if mxl is not None else None
0487     maxlevel = property(_get_maxlevel)
0488
0489     def _get_RC(self):
0490         maxlevel = self.maxlevel
0491         return 1 if maxlevel is not None and  maxlevel.level > Level.WARNING else 0
0492     RC = property(_get_RC)
0493
0494     fdiscreps = property(lambda self:self._get_float("fdiscrep"))
0495     def _get_fdiscmax(self):
0496         fdiscreps_ = self.fdiscreps
0497         return max(fdiscreps_) if len(fdiscreps_) > 0 else -1
0498     fdiscmax   = property(_get_fdiscmax)
0499
0500
0501     def _get_smry(self):
0502         return "%s fdiscmax:%s fdiscreps:%r maxdvmax:%s " % ( self.name, self.fdiscmax, self.fdiscreps, self.maxdvmax  )
0503     smry = property(_get_smry)
0504
0505     def _get_brief(self):
0506         skips = " ".join(self.skips)
0507         gfmt_ = lambda _:"%.4f" % float(_)
0508
0509         if self.maxlevel is None:
0510             return "maxlevel None"
0511         else:
0512             return "maxdvmax:%s  level:%s  RC:%d       skip:%s" % ( gfmt_(self.maxdvmax), self.maxlevel.fn_(self.maxlevel.name), self.RC,  skips )
0513         pass
0514
0515
0516     brief = property(_get_brief)
0517
0518     def __repr__(self):
0519         if len(self.dvs) == 0:
0520             return "\n".join(["ab.%s" % self.name, "no dvs" ])
0521         else:
0522             return "\n".join( ["ab.%s" % self.name, self.brief, self.dvs[0].columns2(), Dv.columns()] + map(repr, filter(None,self.dvs[self.sli]) ) + ["."] )
0523         pass
0524
0525     def __getitem__(self, sli):
0526          self.sli = sli
0527          return self
0528
0529
0530
0531 if __name__ == '__main__':
0532     pass
0533 pass
0534
0535
0536
0537
0538