eic-opticks/ana/nibble.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
0024     In [14]: np.where( seqhis & 0xf << (4*5) != 0 )
0025     Out[14]: (array([5, 6, 7, 8, 9]),)
0026
0027 """
0028 from collections import OrderedDict as odict
0029 from functools import reduce
0030 import numpy as np
0031
0032
0033 msk_ = lambda n:(1 << 4*n) - 1  # msk_(0)=0x0 msk_(1)=0xf msk_(2)=0xff msk_(3)=0xfff
0034
0035 desc_ = lambda _:"\n".join( "%0.16x" % m for m in _ )
0036
0037
0038 def make_msk(n=16):
0039     """
0040     :param n: normally 16
0041     :return msk: array of shape (n,) with nibble masks
0042
0043     In [11]: print( "\n".join( "%0.16x" % m for m in msk ) )
0044     0000000000000000
0045     000000000000000f
0046     00000000000000ff
0047     0000000000000fff
0048     000000000000ffff
0049     00000000000fffff
0050     0000000000ffffff
0051     000000000fffffff
0052     00000000ffffffff
0053     0000000fffffffff
0054     000000ffffffffff
0055     00000fffffffffff
0056     0000ffffffffffff
0057     000fffffffffffff
0058     00ffffffffffffff
0059     0fffffffffffffff
0060     """
0061     msk = np.zeros( n, dtype=np.uint64  )
0062     for i in range(n):
0063         msk[i] = msk_(i)
0064     pass
0065     return msk
0066
0067
0068 nib_ = lambda n:( 0xf << 4*(n-1) ) if n > 0 else 0    # nib_(0)=0xf nib_(1)=0xf0 nib_(2)=0xf00 nib_(3)=0xf000
0069
0070 def make_nib(n):
0071     """
0072     :return nib: array of n uint64 with nibble masks of 4 bits each for slots from 0 to n-1
0073
0074     HMM: because of the zero its missing top nibble, should be using n=17 ?
0075
0076     In [13]: print("\n".join( "%0.16x" % m for m in nib ))
0077     0000000000000000
0078     000000000000000f
0079     00000000000000f0
0080     0000000000000f00
0081     000000000000f000
0082     00000000000f0000
0083     0000000000f00000
0084     000000000f000000
0085     00000000f0000000
0086     0000000f00000000
0087     000000f000000000
0088     00000f0000000000
0089     0000f00000000000
0090     000f000000000000
0091     00f0000000000000
0092     0f00000000000000
0093     """
0094     nib = np.zeros( n, dtype=np.uint64  )
0095     for i in range(n):
0096         nib[i] = nib_(i)
0097     pass
0098     return nib
0099
0100
0101 def count_nibbles(x_):
0102     """
0103     :param x: array of uint64, eq seqhis
0104     :return count: number of filled nibbles
0105
0106     HMM works for arrays but not single values::
0107
0108          ufunc 'right_shift' not supported for the input types,
0109          and the inputs could not be safely coerced to any supported types according to the casting rule ''safe''
0110
0111
0112     https://stackoverflow.com/questions/38225571/count-number-of-zero-nibbles-in-an-unsigned-64-bit-integer
0113
0114     In [14]: count_nibbles(msk)
0115     Out[14]: array([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15], dtype=uint64)
0116
0117     In [15]: count_nibbles(nib)
0118     Out[15]: array([0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1], dtype=uint64)
0119     """
0120     x = x_.copy()
0121
0122     ## gather the zeroness (the OR of all 4 bits)
0123     x |= x >> 1               # or-with-1bit-right-shift-self is or-of-each-consequtive-2-bits
0124     x |= x >> 2               # or-with-2bit-right-shift-self is or-of-each-consequtive-4-bits in the lowest slot
0125     x &= 0x1111111111111111   # pluck the zeroth bit of each of the 16 nibbles
0126
0127     x = (x + (x >> 4)) & 0xF0F0F0F0F0F0F0F    # sum occupied counts of consequtive nibbles, and pluck them
0128     count = (x * 0x101010101010101) >> 56     #  add up byte totals into top byte,  and shift that down to pole 64-8 = 56
0129
0130     return count
0131
0132
0133
0134 tst_ = lambda i:reduce(lambda a,b:a|b, map(lambda n:0x1 << (4*n), range(i) ),0)  # marching 1s to mimic seqhis nibbles
0135 tst0_ = lambda n:(1 << 4*n)   # leaves unoccupied nibbles to the right, so not valid standin for seqhis nibbles
0136
0137 def test_tst(n):
0138     """
0139      0 : 0000000000000000
0140      1 : 0000000000000001
0141      2 : 0000000000000011
0142      3 : 0000000000000111
0143      4 : 0000000000001111
0144      5 : 0000000000011111
0145      6 : 0000000000111111
0146      7 : 0000000001111111
0147      8 : 0000000011111111
0148      9 : 0000000111111111
0149     10 : 0000001111111111
0150     11 : 0000011111111111
0151     12 : 0000111111111111
0152     13 : 0001111111111111
0153     14 : 0011111111111111
0154     15 : 0111111111111111
0155     """
0156     print("\n".join(map(lambda i:"%2d : %0.16x" % (i, tst_(i) ), range(n))))
0157
0158
0159 def make_tst(n):
0160    """
0161    Create array of 64bit uints as standin for seqhis
0162    photon histories with 4 bits per recorded step point
0163    """
0164    tst = np.zeros( n, dtype=np.uint64 )
0165    for i in range(n):
0166        tst[i] = tst_(i)
0167    pass
0168    xpo = np.arange( len(tst), dtype=np.int32 )
0169    return tst, xpo
0170
0171
0172 def make_tst2(n):
0173     """
0174     Create array of 64bit uints as standin for seqhis
0175     photon histories with 4 bits per recorded step point
0176
0177     In [3]: print(desc_(tst))
0178     0000000000000000
0179     0000000000000001
0180     0000000000000011
0181     0000000000000111
0182     0000000000001111
0183     0000000000011111
0184     0000000000111111
0185     0000000001111111
0186     0000000011111111
0187     0000000111111111
0188     0000001111111111
0189     0000011111111111
0190     0000111111111111
0191     0001111111111111
0192     0011111111111111
0193     0111111111111111
0194     0111111111111111
0195     0011111111111111
0196     0001111111111111
0197     0000111111111111
0198     0000011111111111
0199     0000001111111111
0200     0000000111111111
0201     0000000011111111
0202     0000000001111111
0203     0000000000111111
0204     0000000000011111
0205     0000000000001111
0206     0000000000000111
0207     0000000000000011
0208     0000000000000001
0209     0000000000000000
0210     """
0211     tst = np.zeros( 2*n, dtype=np.uint64 )
0212     xpo = np.zeros( 2*n, dtype=np.uint64 )
0213     for i in range(n):
0214         tst[i] = tst_(i)
0215         xpo[i] = i
0216     for i in range(n):
0217         tst[n+i] = tst_(n-1-i)
0218         xpo[n+i] = n-1-i
0219     pass
0220     return tst, xpo
0221
0222
0223 def dump_msk_nib(msk, nib):
0224     print(" %16s : %16s : %16s " % ( "nib", "msk", "~msk" ))
0225     for i in range(len(msk)):
0226         print(" %16x : %16x : %16x " % ( nib[i], msk[i], ~msk[i]  ))
0227     pass
0228
0229
0230
0231
0232 if __name__ == '__main__':
0233     n = 16
0234     msk = make_msk(n)
0235     nib = make_nib(n)
0236     dump_msk_nib(msk, nib)
0237
0238     #test_tst(n)
0239     tst,xpo = make_tst2(n)
0240
0241     ## number of occupied nibbles in the seqhis
0242     npo = np.zeros(len(tst), dtype=np.int32)
0243     qpo = odict()
0244     wpo = odict()
0245
0246     for i in range(16):
0247         wpo[i] = np.where( tst >> (4*i) != 0 )
0248         npo[wpo[i]] = i+1
0249         """
0250         wpo[i]
0251            right shift by i nibbles and see if any bits remain
0252
0253         Notice the overwriting into npo, this works because
0254         are using ascending order
0255
0256         """
0257         #qpo[i] = np.where(np.logical_and(np.logical_and( tst & ~msk[i] == 0, tst & msk[i] == tst ), tst & nib[i] != 0  ))[0]
0258         #qpo[i] = np.where(reduce(np.logical_and, (tst & ~msk[i] == 0, tst & msk[i] == tst, tst & nib[i] != 0  )))
0259         #qpo[i] = np.where(reduce(np.logical_and, (tst & ~msk[i] == 0, tst & nib[i] != 0  )))
0260         #qpo[i] = np.where(np.logical_and(tst & ~msk[i] == 0, tst & nib[i] != 0  ))
0261         qpo[i] = np.where( (tst & ~msk[i] == 0) & (tst & nib[i] != 0) )
0262     pass
0263     """
0264     qpo[i] (i in 1..16)  gives tst indices having each of the possible number of occupied nibbles
0265
0266     tst & ~msk[i] == 0
0267         and-ing with anti-mask[i] removes high side nibbles,
0268         but leaves all those below
0269
0270     tst & nib[i] != 0
0271         requiring something in nibble i cuts out all those below
0272
0273     All the indices selected have the same number of occupied nibbles
0274
0275     This may be specialized to the seqhis case where never get
0276     zero nibbles on the low side.
0277     """
0278
0279     print(" %2s : %16s : %2s : %2s " % ( "i", "tst", "xp", "np" ))
0280     for i in range(len(tst)):
0281         print(" %2d : %16x : %2d : %2d " % ( i, tst[i],  xpo[i], npo[i]  ))
0282     pass
0283
0284     assert np.all( npo == xpo )
0285
0286
0287     print(" %2s : %16s : %s " % ( "i", "tst", "qpo" ))
0288     for i in range(16):
0289         print(" %2d : %16x : %s " % (  i, tst[i],  repr(qpo[i]) ))
0290     pass
0291
0292