eic-opticks/sysrap/xfold.py

0001 #!/usr/bin/env python
0002
0003 import numpy as np
0004 import os, re, logging
0005 log = logging.getLogger(__name__)
0006
0007 from opticks.ana.fold import Fold
0008 from opticks.sysrap.stag import stag
0009 from opticks.u4.U4Stack import U4Stack
0010 from opticks.ana.p import *
0011
0012
0013 class XFold(object):
0014     tag = stag()
0015     stack = U4Stack()
0016
0017     @classmethod
0018     def BaseSymbol(cls, xf):
0019         """
0020         :param xf: Fold instance
0021         :return symbol: "A" or "B" : A for Opticks, B for Geant4
0022         """
0023         CX = xf.base.find("CX") > -1
0024         U4 = xf.base.find("U4") > -1
0025         assert CX ^ U4  # exclusive-OR
0026         symbol = "A" if CX else "B"
0027         return symbol
0028
0029     @classmethod
0030     def Ident(cls, x):
0031         bsym = cls.BaseSymbol(x)
0032         ident = cls.tag if bsym == "A" else cls.stack
0033         return ident
0034
0035     def __init__(self, x, symbol=None):
0036         """
0037         :param x: Fold instance, eg either a or b
0038         """
0039         t = stag.Unpack(x.tag) if hasattr(x,"tag") else None
0040         f = getattr(x, "flat", None)
0041         n = stag.NumStarts(t) if not t is None else None
0042         log.info("XFold before stag.StepSplit")
0043         ts,fs = stag.StepSplit(t,x.flat) if not t is None else None
0044         log.info("XFold after stag.StepSplit")
0045         ident = self.Ident(x)
0046
0047         xsymbol = self.BaseSymbol(x)
0048         if xsymbol == "B":
0049             t2 = self.stack.make_stack2tag_mapped(t)
0050             ts2 = self.stack.make_stack2tag_mapped(ts)
0051         elif xsymbol == "A":
0052             t2 = self.stack.make_tag2stack_mapped(t)
0053             ts2 = self.stack.make_tag2stack_mapped(ts)
0054         else:
0055             assert 0
0056         pass
0057         if symbol != xsymbol:
0058             log.error("using unconventional symbol %s xsymbol %s " % (symbol, xsymbol))
0059         pass
0060         #assert symbol == xsymbol
0061
0062         self.x = x       # Fold instance, called x because it is usually "a" or "b"
0063         self.t = t       # (num_photon, SLOTS)  : unpacked consumption tag/stack enumeration integers
0064         self.f = f       # (num_photon, SLOTS)  : uniform rand
0065         self.n = n       # (num_photon,)        : number of steps
0066         self.t2 = t2     # (num_photon, SLOTS)  : native enumeration mapped to the other enumeration
0067
0068         self.ts = ts     # (num_photon, 7, 10)  : 7 is example max_steps which should match A.n.max(), 10 is example max rand per step input
0069         self.fs = fs     # (num_photon, 7, 10)  : flat random values split by steps
0070         self.ts2 = ts2   # (num_photon, 7, 10)  : native enumeration mapped to the other enumeration
0071
0072         self.ident = ident  # A:stag OR B:U4Stack instance for providing names to enumeration codes
0073         self.symbol = symbol
0074         self.xsymbol = xsymbol
0075         self.idx = 0
0076         self.flavor = ""
0077
0078     def __call__(self, idx):
0079         self.idx = idx
0080         self.flavor = "call"
0081         return self
0082
0083     def __getitem__(self, idx):
0084         self.idx = idx
0085         self.flavor = "getitem"
0086         return self
0087
0088     def header(self):
0089         idx = self.idx
0090         seqhis = seqhis_(self.x.seq[idx,0])
0091         return "%s(%d) : %s" % (self.symbol, idx, seqhis)
0092
0093     def rbnd__(self):
0094         return boundary___(self.x.record[self.idx])
0095     def rori__(self):
0096         return orient___(self.x.record[self.idx])
0097     def rpri__(self):
0098         return primIdx___(self.x.record[self.idx])
0099     def rins__(self):
0100         return instanceId__(self.x.record[self.idx])
0101
0102
0103     def rbnd_(self):
0104         bb_ = self.rbnd__()
0105         oo_ = self.rori__()
0106         pp_ = self.rpri__()
0107         ii_ = self.rins__()
0108         assert bb_.shape == oo_.shape
0109         assert bb_.shape == pp_.shape
0110         assert bb_.shape == ii_.shape
0111
0112         wp = np.where( bb_ > 0 )  # mask boundary zero to skip unsets
0113         bb = bb_[wp]
0114         oo = oo_[wp]
0115         pp = pp_[wp]
0116         ii = ii_[wp]
0117
0118         pm = ["+","-"]
0119         lines = [ "%s %-40s %-50s"%(pm[oo[i]],cf.sim.bndnamedict.get(bb[i]), cf.primIdx_meshname_dict.get(pp[i])) for i in range(len(bb))]
0120
0121         lines += [""]
0122         lines += ["- : against the normal (ie inwards from omat to imat)"]
0123         lines += ["+ : with the normal (ie outwards from imat to omat)"]
0124         return lines
0125
0126
0127     def rbnd(self):
0128         return "\n".join(self.rbnd_())
0129
0130     def content(self):
0131         lines = []
0132         members = "t t2 n ts fs ts2".split()
0133         for mem in members:
0134             arr = getattr(self, mem)
0135             name = "%s.%s" % ( self.symbol, mem)
0136             line = "%10s : %s " % (name, str(arr.shape))
0137             lines.append(line)
0138         pass
0139         return "\n".join(lines)
0140
0141     def body(self):
0142         return self.ident.label(self.t[self.idx],self.f[self.idx])
0143
0144     def identification(self):
0145         return "%s : %s" % (self.symbol, self.x.base)
0146
0147     def call_repr(self):
0148         return "\n".join([self.header(), self.content(), self.body()])
0149
0150     def getitem_repr(self):
0151         return "\n".join([self.header(), self.rbnd()])
0152
0153     def __repr__(self):
0154         if self.flavor == "call":
0155             rep = self.call_repr()
0156         elif self.flavor == "getitem":
0157             rep = self.getitem_repr()
0158         else:
0159             rep = self.call_repr()
0160         pass
0161         return rep
0162
0163
0164
0165
0166 if __name__ == '__main__':
0167     logging.basicConfig(level=logging.INFO)
0168
0169     a = Fold.Load("$A_FOLD", symbol="a")
0170     b = Fold.Load("$B_FOLD", symbol="b")
0171     A = XFold(a, symbol="A")
0172     B = XFold(b, symbol="B")
0173
0174