sysrap/tests/SPMT_test.py

0001 #!/usr/bin/env python
0002 """
0003 SPMT_test.py
0004 ==============
0005
0006 """
0007 import os, textwrap, numpy as np
0008 from opticks.ana.fold import Fold
0009 import matplotlib.pyplot as plt
0010 SIZE = np.array([1280,720])
0011 np.set_printoptions(edgeitems=16)
0012
0013
0014 def old_check_test(s):
0015     """
0016     """
0017     qi = s.test.get_pmtid_qe
0018     qc = s.test.get_pmtcat_qe
0019     ct = s.test.get_pmtcat
0020
0021     fig, ax = plt.subplots(1, figsize=[12.8, 7.2] )
0022
0023     for i in range(3):
0024         ax.plot( qc[i,:,0], qc[i,:,1], label="qc[%d]"%i )
0025     pass
0026     ax.legend()
0027
0028     #for pmtid in range(25):
0029     #    ax.plot( qi[pmtid,:,0], qi[pmtid,:,1], label=pmtid )
0030     #pass
0031     fig.show()
0032 pass
0033
0034
0035
0036
0037 '-', '--', '-.', ':', ''
0038
0039 plotmap = {
0040    'As':["r","-"],
0041    'Ap':["r","--"],
0042    'Aa':["r","-."],
0043    'A_':["r",":"],
0044    'Rs':["g","-"],
0045    'Rp':["g","--"],
0046    'Ra':["g","-."],
0047    'R_':["g",":"],
0048    'Ts':["b","-"],
0049    'Tp':["b","--"],
0050    'Ta':["b","-."],
0051    'T_':["b",":"],
0052    }
0053
0054 color_     = lambda _:plotmap.get(_,["k","."])[0]
0055 linestyle_ = lambda _:plotmap.get(_,["k","."])[1]
0056
0057
0058 class ART(object):
0059
0060     @classmethod
0061     def check_nan(cls,f):
0062         print("\ncheck_nan f.base %s " % f.base)
0063         qwns = "args spec extra ARTE stack ll comp art nstack nll ncomp nart".split()
0064         for qwn in qwns:
0065             q = getattr(f, qwn, None)
0066             if q is None: continue
0067             expr = "np.where( np.isnan(f.%s.ravel()) )[0]" % qwn
0068             nan = eval(expr)
0069             print(" %-50s : %s " % (expr,eval(expr)))
0070         pass
0071
0072     @classmethod
0073     def check_domain(cls, f):
0074         lpmtid_domain = f.lpmtid_domain
0075         lpmtcat_domain = f.lpmtcat_domain
0076         assert len(lpmtid_domain) == len(lpmtcat_domain)
0077
0078         mct_domain = f.mct_domain
0079         st_domain = f.st_domain
0080         assert len(mct_domain) == len(st_domain)
0081
0082
0083
0084     def __init__(self, f):
0085         self.f = f
0086
0087         self.check_nan(f)
0088         self.check_domain(f)
0089
0090         args = f.args.squeeze()
0091         spec = f.spec.squeeze()
0092         extra = f.extra.squeeze()
0093         ARTE = f.ARTE.squeeze()
0094
0095         stack = f.stack.squeeze()
0096         ll = f.ll.squeeze()
0097         comp = f.comp.squeeze()
0098         art = f.art.squeeze()
0099
0100         num_lpmtid = len(f.lpmtid_domain)
0101         num_mct = len(f.mct_domain)
0102         art_shape = (num_lpmtid, num_mct, 4, 4) if num_mct > 1 else (num_lpmtid, 4, 4)
0103         assert art.shape == art_shape, (art_shape, art.shape)
0104
0105         PMTIDX = np.fromstring(os.environ.get("PMTIDX","0"),dtype=np.int64, sep=",")
0106         lpmtid = f.lpmtid_domain[PMTIDX]
0107         lpmtcat = f.lpmtcat_domain[PMTIDX]
0108
0109         OPT = "A_,R_,T_,As,Rs,Ts,Ap,Rp,Tp,Aa,Ra,Ta"
0110         opt = os.environ.get("OPT", OPT)
0111
0112
0113         expr = "np.c_[PMTIDX,lpmtid,lpmtcat].T"
0114         etab = eval(expr)
0115         mtab = "PMTIDX %s lpmtid %s lpmtcat %s " % ( str(PMTIDX), str(lpmtid), str(lpmtcat) )
0116         title = "%s : OPT %s \n%s" % (f.base, opt, etab )
0117
0118         print(title)
0119
0120
0121         fig, ax = plt.subplots(1, figsize=SIZE/100.)
0122         fig.suptitle(title)
0123
0124         for i, pmtidx in enumerate(PMTIDX):
0125             As   = art[pmtidx,:,0,0]
0126             Ap   = art[pmtidx,:,0,1]
0127             Aa   = art[pmtidx,:,0,2]
0128             A_   = art[pmtidx,:,0,3]
0129
0130             Rs   = art[pmtidx,:,1,0]
0131             Rp   = art[pmtidx,:,1,1]
0132             Ra   = art[pmtidx,:,1,2]
0133             R_   = art[pmtidx,:,1,3]
0134
0135             Ts   = art[pmtidx,:,2,0]
0136             Tp   = art[pmtidx,:,2,1]
0137             Ta   = art[pmtidx,:,2,2]
0138             T_   = art[pmtidx,:,2,3]
0139
0140             SF     = art[pmtidx,:,3,0]
0141             wl     = art[pmtidx,:,3,1]
0142             ARTa   = art[pmtidx,:,3,2]
0143             mct    = art[pmtidx,:,3,3]
0144
0145             if i == 0:
0146                 label_ = lambda _:_
0147             else:
0148                 label_ = lambda _:None
0149             pass
0150
0151             if "As" in opt:ax.plot(  mct, As, label=label_("As"), color=color_("As"), linestyle=linestyle_("As") )
0152             if "Ap" in opt:ax.plot(  mct, Ap, label=label_("Ap"), color=color_("Ap"), linestyle=linestyle_("Ap"))
0153             if "Aa" in opt:ax.plot(  mct, Aa, label=label_("Aa"), color=color_("Aa"), linestyle=linestyle_("Aa"))
0154             if "A_" in opt:ax.plot(  mct, A_, label=label_("A_"), color=color_("A_"), linestyle=linestyle_("A_"))
0155
0156             if "Rs" in opt:ax.plot(  mct, Rs, label=label_("Rs"), color=color_("Rs"), linestyle=linestyle_("Rs"))
0157             if "Rp" in opt:ax.plot(  mct, Rp, label=label_("Rp"), color=color_("Rp"), linestyle=linestyle_("Rp"))
0158             if "Ra" in opt:ax.plot(  mct, Ra, label=label_("Ra"), color=color_("Ra"), linestyle=linestyle_("Ra"))
0159             if "R_" in opt:ax.plot(  mct, R_, label=label_("R_"), color=color_("R_"), linestyle=linestyle_("R_"))
0160
0161             if "Ts" in opt:ax.plot(  mct, Ts, label=label_("Ts"), color=color_("Ts"), linestyle=linestyle_("Ts"))
0162             if "Tp" in opt:ax.plot(  mct, Tp, label=label_("Tp"), color=color_("Tp"), linestyle=linestyle_("Tp"))
0163             if "Ta" in opt:ax.plot(  mct, Ta, label=label_("Ta"), color=color_("Ta"), linestyle=linestyle_("Ta"))
0164             if "T_" in opt:ax.plot(  mct, T_, label=label_("T_"), color=color_("T_"), linestyle=linestyle_("T_"))
0165
0166             if "SF" in opt:ax.plot(  mct, SF, label=label_("SF"), color=color_("SF"), linestyle=linestyle_("SF"))
0167             if "wl" in opt:ax.plot(  mct, wl, label=label_("wl"), color=color_("wl"), linestyle=linestyle_("wl") )
0168             if "ARTa" in opt:ax.plot(  mct, ARTa, label=label_("ARTa"), color=color_("ARTa"), linestyle=linestyle_("ARTa") )
0169             if "mct" in opt:ax.plot(  mct, mct, label=label_("mct"), color=color_("mct"), linestyle=linestyle_("mct") )
0170         pass
0171         ax.legend()
0172         fig.show()
0173
0174
0175
0176 class TESTFOLD(object):
0177     """
0178     np.c_[np.unique(t.get_lpmtcat_from_lpmtidx, return_counts=True)]
0179     np.c_[np.unique(t.get_lpmtcat_from_lpmtidx, return_counts=True)][:,1].sum()
0180
0181     t.get_qescale_from_lpmtidx.shape
0182     t.get_qescale_from_lpmtidx.min()
0183     t.get_qescale_from_lpmtidx.max()
0184
0185     t.get_qescale_from_lpmtidx[:17612].min()         ## CD_LPMT
0186     t.get_qescale_from_lpmtidx[:17612].max()
0187
0188     t.get_qescale_from_lpmtidx[17612:17612+2400].min()       ## WP_PMT
0189     t.get_qescale_from_lpmtidx[17612:17612+2400].max()
0190
0191     t.get_qescale_from_lpmtidx[17612+2400:17612+2400+324].min()   ## WP_ATM_LPMT
0192     t.get_qescale_from_lpmtidx[17612+2400:17612+2400+324].max()
0193
0194     t.get_qescale_from_lpmtidx[17612+2400+324:17612+2400+324+5].min()  ## WP_WAL_PMT
0195     t.get_qescale_from_lpmtidx[17612+2400+324:17612+2400+324+5].max()
0196
0197
0198     np.all( t.get_lcqs_from_lpmtidx[:,0] == t.get_lpmtcat_from_lpmtidx )
0199     np.all( t.get_lcqs_from_lpmtidx[:,1].view(np.float32) == t.get_qescale_from_lpmtidx )
0200
0201
0202     # S_PMT
0203
0204     np.all( s.s_qescale == t.get_s_qescale_from_spmtid )  # consistency check between SPMT.h input and output
0205
0206     t.get_s_qeshape.shape
0207     t.get_s_qeshape
0208
0209     """
0210
0211     @classmethod
0212     def EXPR(cls):
0213         return list(map(str.strip,textwrap.dedent(cls.__doc__).split("\n")))
0214
0215     def __init__(self, t, s ):
0216         self.t = t
0217         self.s = s
0218
0219     def __repr__(self):
0220         lines = []
0221
0222         t = self.t
0223         s = self.s
0224
0225         for expr in self.EXPR():
0226             lines.append(expr)
0227             if expr == "" or expr[0] == "#": continue
0228             lines.append(repr(eval(expr)))
0229         pass
0230         return "\n".join(lines)
0231
0232
0233 class QESHAPEPLOT(object):
0234     def __init__(self, t, s, title="QESHAPEPLOT" ):
0235         self.t = t
0236         self.s = s
0237
0238         qesh = t.get_s_qeshape[0]    # interpolated values
0239         s_qesh = s.s_qeshape[0,:-1]  # prop values, excluding special last value
0240
0241         fig, ax = plt.subplots(1, figsize=SIZE/100.)
0242         fig.suptitle(title)
0243
0244         ax.plot( qesh[:,0], qesh[:,1], label="qesh" )
0245
0246         sel = np.logical_and( s_qesh[:,0] >= qesh[:,0].min(), s_qesh[:,0] < qesh[:,0].max() )
0247         # select the prop values within the interpolated range
0248         ax.scatter( s_qesh[sel,0], s_qesh[sel,1], label="s_qesh" )
0249
0250         ax.legend()
0251         fig.show()
0252
0253
0254 if __name__ == '__main__':
0255
0256     TEST = os.environ.get("TEST","testfold")
0257
0258     s = Fold.Load("$FOLD/spmt", symbol="s")
0259     print(repr(s))
0260
0261     if TEST == "c4scan":
0262         f = Fold.Load("$FOLD/testfold/c4scan", symbol="f")
0263         print(repr(f))
0264         a = ART(f)
0265     elif TEST == "testfold":
0266
0267         t = Fold.Load("$FOLD/testfold", symbol="t")
0268         print(repr(t))
0269
0270         tf = TESTFOLD(t, s)
0271         print(repr(tf))
0272
0273
0274         QESHAPEPLOT(t, s)
0275
0276
0277     elif TEST == "spmt":
0278         pass
0279     pass
0280
0281
0282