eic-opticks/ana/dbgseed.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 This script is invoked from oks-dbgseed, which
0024 first creates the seeds with commands like::
0025
0026    GGeoViewTest --dbgseed --trivial --cerenkov --compute
0027    GGeoViewTest --dbgseed --trivial --cerenkov
0028
0029 Looks like on Linux the interop photon buffer is being seeded but is
0030 then getting overwritten ??
0031
0032 ::
0033
0034     In [22]: cc[:10]
0035     Out[22]:
0036     array([[  0,  80],
0037            [  1, 108],
0038            [  2,  77],
0039            [  3,  30],
0040            [  4,  99],
0041            [  5, 105],
0042            [  6, 106],
0043            [  7,  85],
0044            [  8,  94],
0045            [  9,  29]])
0046
0047     In [23]: ii[42:42+10]
0048     Out[23]:
0049     array([[  0,  80],
0050            [  1, 108],
0051            [  2,  77],
0052            [  3,  30],
0053            [  4,  99],
0054            [  5, 105],
0055            [  6, 106],
0056            [  7,  85],
0057            [  8,  94],
0058            [  9,  29]])
0059
0060
0061     n [11]: ii[:43]
0062     Out[11]:
0063     array([[-965867012,          5],
0064            [-965852569,     192504],
0065            [-961979442,          4],
0066            [-961978135,         54],
0067            [-961921942,         49],
0068            [-961761495,        286],
0069            ...
0070            [-918184732,         18],
0071            [-918183179,         21],
0072            [-918181729,        108],
0073            [-918180539,       7062],
0074            [-918178532,         12],
0075            [-918178359,     368294],
0076            [         0,         80]])
0077
0078
0079 Things go wrong from item 11883::
0080
0081     In [41]: i[11880:12000,0,0].view(np.int32)
0082     Out[41]:
0083     array([       160,        160,        160, -965867012, -965867012,
0084            -965867012, -965867012, -965867012, -965852569, -965852569,
0085            -965852569, -965852569, -965852569, -965852569, -965852569,
0086            -965852569, -965852569, -965852569, -965852569, -965852569,
0087            -965852569, -965852569, -965852569, -965852569, -965852569,
0088            -965852569, -965852569, -965852569, -965852569, -965852569,
0089            -965852569, -965852569, -965852569, -965852569, -965852569,
0090            -965852569, -965852569, -965852569, -965852569, -965852569,
0091
0092
0093     In [42]: i.shape
0094     Out[42]: (612841, 4, 4)
0095
0096 ::
0097
0098     In [49]: i[11882]
0099     Out[49]:
0100     array([[      0.   ,       0.   ,       0.   ,       0.   ],
0101            [-806767.188,   -8289.5  ,       0.   ,       0.   ],
0102            [ -15244.033, -806768.125,   -8315.324,  -15232.615],
0103            [-806767.188,   -8289.5  ,       0.   ,       0.   ]], dtype=float32)
0104
0105     In [50]: i[11883]
0106     Out[50]:
0107     array([[ -15234.496, -806782.125,   -8417.453,  -15215.982],
0108            [-806768.062,   -8395.   ,         nan,       0.   ],
0109            [ -15234.496, -806782.125,   -8290.   ,  -15215.982],
0110            [-806768.062,   -8267.547,         nan,       0.   ]], dtype=float32)
0111
0112
0113     In [47]: i[11882].view(np.int32)
0114     Out[47]:
0115     array([[       160,          0,          0,          0],
0116            [-918227213, -972978688,     227185,     227189],
0117            [-965857246, -918227198, -972952244, -965868938],
0118            [-918227213, -972978688,     227189,     227191]], dtype=int32)
0119
0120     In [48]: i[11883].view(np.int32)
0121     Out[48]:
0122     array([[-965867012, -918226974, -972847664, -965885970],
0123            [-918227199, -972870656,         -1,     317936],
0124            [-965867012, -918226974, -972978176, -965885970],
0125            [-918227199, -973001168,         -1,     317944]], dtype=int32)
0126
0127
0128
0129
0130 """
0131
0132 import os, sys, datetime, logging, numpy as np
0133 from opticks.ana.base import opticks_main
0134 from opticks.ana.nbase import count_unique
0135
0136 log = logging.getLogger(__name__)
0137
0138
0139 def stmp_(st, fmt="%Y%m%d-%H%M"):
0140     return datetime.datetime.fromtimestamp(st.st_ctime).strftime(fmt)
0141
0142 def stamp_(path, fmt="%Y%m%d-%H%M"):
0143     try:
0144         st = os.stat(path)
0145     except OSError:
0146         return "FILE-DOES-NOT-EXIST"
0147     return stmp_(st, fmt=fmt)
0148
0149 def x_(_):
0150     p = os.path.expandvars(_)
0151     st = stamp_(p)
0152     log.info( " %s -> %s (%s) " % (_, p, st))
0153     return p
0154
0155 def check_dbgseed(a,g):
0156     """
0157     The seeds should be genstep_id from 0:num_genstep-1
0158     """
0159     aa = count_unique(a[:,0,0].view(np.int32))
0160     assert np.all(aa[:,0] == np.arange(0,len(aa)))
0161
0162     xx = g[:,0,3].view(np.int32)   ## photons per genstep
0163
0164     assert len(aa) == len(xx)
0165     assert np.all(aa[:,1] == xx)
0166
0167
0168 if __name__ == '__main__':
0169     args = opticks_main(src="torch", tag="1", det="dayabay")
0170
0171     np.set_printoptions(suppress=True, precision=3)
0172
0173     cpath = x_("$TMP/dbgseed_compute.npy")
0174     ipath = x_("$TMP/dbgseed_interop.npy")
0175
0176     log.info("cpath : %s " % cpath)
0177     log.info("ipath : %s " % ipath)
0178
0179     if not(os.path.exists(cpath) and os.path.exists(ipath)):
0180         log.warning("SKIP due to missing path")
0181         sys.exit(0)  ## very particular test, not a standardized test yet so dont treat as a fail
0182
0183     c = np.load(cpath)
0184     i = np.load(ipath)
0185
0186     log.info(" c : %s " % repr(c.shape) )
0187     log.info(" i : %s " % repr(i.shape) )
0188
0189     cj = c[:,0,0].view(np.int32)
0190     ij = i[:,0,0].view(np.int32)
0191
0192
0193     if args.src in ("cerenkov", "scintillation"):
0194         g = np.load(x_("$OPTICKS_DATA_DIR/gensteps/%s/%s/%s.npy" % (args.det,args.src,args.tag) ))
0195     elif args.src == "torch":
0196         g = np.load(x_("$TMP/torchdbg.npy"))
0197     else:
0198         assert 0, args.src
0199
0200     if args.src == "torch":
0201         print "cj", cj[99500:]
0202         print "ij", ij[99500:]
0203     pass
0204
0205     check_dbgseed(c,g)
0206     check_dbgseed(i,g)
0207
0208     ii = count_unique(i[:,0,0].view(np.int32))
0209     cc = count_unique(c[:,0,0].view(np.int32))
0210
0211
0212
0213
0214
0215
0216