g4cx/tests/ab.py

0001 #!/usr/bin/env python
0002 """
0003 ab.py (formerly G4CXSimulateTest_ab.py) : random aligned comparison
0004 ======================================================================
0005
0006 Usage::
0007
0008    gx
0009    ./ab.sh
0010
0011
0012     In [15]: np.all( A.ts == B.ts2[:,:12,:29] )
0013     Out[15]: False
0014
0015     In [17]: A.ts.shape
0016     Out[17]: (10000, 12, 29)
0017
0018     In [18]: B.ts2.shape
0019     Out[18]: (10000, 13, 29)
0020
0021
0022 Warning from comparing different shapes::
0023
0024     In [16]: np.all( A.ts == B.ts2 )
0025     /Users/blyth/miniconda3/bin/ipython:1: DeprecationWarning: elementwise comparison failed; this will raise an error in the future.
0026       #!/Users/blyth/miniconda3/bin/python
0027     Out[16]: False
0028
0029 ::
0030
0031     In [19]: np.where( A.ts != B.ts2[:,:12,:29] )
0032     Out[19]:
0033     (array([9964, 9964, 9964, 9964, 9964, 9964]),
0034      array([ 2,  4,  6,  8, 10, 11]),
0035      array([5, 5, 5, 5, 5, 3]))
0036
0037
0038 Mis-aligned from truncation ?::
0039
0040     In [20]: seqhis_(a.seq[9964,0])
0041     Out[20]: 'TO BT SC BR BR BR BR BR BR BR'
0042
0043     In [21]: seqhis_(b.seq[9964,0])
0044     Out[21]: 'TO BT SC BR BR BR BR BR BR BR'
0045
0046
0047
0048 """
0049 import numpy as np
0050 import logging
0051 log = logging.getLogger(__name__)
0052
0053 from opticks.ana.fold import Fold
0054 from opticks.ana.p import *
0055 from opticks.ana.eprint import eprint, epr, edv
0056
0057 from opticks.sysrap.xfold import XFold
0058 from opticks.sysrap.ABR import ABR
0059 from opticks.sysrap.stag import stag
0060 from opticks.u4.U4Stack import U4Stack
0061
0062 np.set_printoptions(edgeitems=16)
0063
0064 tag = stag()
0065 stack = U4Stack()
0066 SLOTS = stag.SLOTS
0067
0068
0069 if __name__ == '__main__':
0070     logging.basicConfig(level=logging.INFO)
0071     quiet = True
0072     a = Fold.Load("$A_FOLD", symbol="a", quiet=quiet) if "A_FOLD" in os.environ else None
0073     b = Fold.Load("$B_FOLD", symbol="b", quiet=quiet) if "B_FOLD" in os.environ else None
0074
0075     ## SUSPECT photon element (1,3) is starting uninitialized in input photon running
0076     ## forcing the below kludge to be able to compare
0077     np.all( b.record[:,0,1,3] == 0. )
0078     np.all( a.record[:,0,1,3] == 1. )
0079     a.record[:,0,1,3] = 1.   # kludge
0080     b.record[:,0,1,3] = 1.   # kludge
0081
0082     ab = (not a is None) and (not b is None)
0083
0084     print("-------- after Fold.Load" )
0085
0086     A = XFold(a, symbol="A") if not a is None else None
0087
0088     print("-------- after XFold.A" )
0089
0090     B = XFold(b, symbol="B") if not b is None else None
0091
0092     print("-------- after XFold.B" )
0093
0094     AB = ABR(A,B) if ab else None
0095
0096     print("-------- after ABR " )
0097
0098     if ab:
0099         im = epr("im = np.abs(a.inphoton - b.inphoton).max()", globals(), locals() )
0100         pm = epr("pm = np.abs(a.photon - b.photon).max()",     globals(), locals() )
0101         rm = epr("rm = np.abs(a.record - b.record).max()",     globals(), locals() )
0102         sm = epr("sm = np.all( a.seq[:,0] == b.seq[:,0] )",    globals(), locals() )
0103         ww = epr("ww = np.where( a.seq[:,0] != b.seq[:,0] )[0]",    globals(), locals() )
0104
0105         we_ = "we = np.where( A.t.view('|S%(SLOTS)s') != B.t2.view('|S%(SLOTS)s') )[0]" % locals()  # eg stag.h/stag::SLOTS = 64
0106         we = epr(we_,  globals(), locals() )
0107
0108         wm_ = "wm = np.where( A.t.view('|S%(SLOTS)s') == B.t2.view('|S%(SLOTS)s') )[0]" % locals()  # eg stag.h/stag::SLOTS = 64
0109         wm = epr(wm_,  globals(), locals() )
0110
0111
0112         wa = epr("wa = np.unique(np.where( np.abs(a.record - b.record ) > 0.05 )[0])", globals(), locals() )
0113
0114         print("---0---")
0115         eprint("np.all( A.ts == B.ts2 )", globals(), locals() )
0116         print("---1---")
0117         eprint("np.all( A.ts2 == B.ts )", globals(), locals() )
0118         print("---2---")
0119
0120         assert (a.inphoton - b.inphoton).max() < 1e-3
0121         print("---3---")
0122
0123
0124         epr("o = cuss(a.seq[:,0])",  globals(), locals(), rprefix="\n" )
0125         edv("a.record[w0,1,0,2] - b.record[w0,1,0,2] # point 1 z", globals(), locals())
0126         edv("a.record[w3,1,0,2] - b.record[w3,1,0,2] # point 1 z", globals(), locals())
0127     pass
0128
0129     print("./U4RecorderTest_ab.sh ## u4t ")
0130     w = epr("w = np.unique(np.where( np.abs(a.photon - b.photon) > 0.1 )[0])", globals(), locals() )
0131     s = epr("s = a.seq[w,0]", globals(), locals() )
0132     o = epr("o = cuss(s,w)", globals(), locals() , rprefix="\n")
0133
0134     #epr("w1", globals(), locals() )
0135     #abw0 = epr("abw0 = a.photon[w0,:4] - b.photon[w0,:4]", globals(), locals(), rprefix="\n" )
0136
0137     epr("a.base", globals(), locals() )
0138     epr("b.base", globals(), locals() )
0139 pass
0140
0141