eic-opticks/ana/edges.py

0001 #!/usr/bin/env python
0002
0003 import logging
0004 log = logging.getLogger(__name__)
0005 import numpy as np
0006 import matplotlib.pyplot as plt
0007
0008 ri = np.array([
0009        [ 1.55 ,  1.478],
0010        [ 1.795,  1.48 ],
0011        [ 2.105,  1.484],
0012        [ 2.271,  1.486],
0013        [ 2.551,  1.492],
0014        [ 2.845,  1.496],
0015        [ 3.064,  1.499],
0016        [ 4.133,  1.526],
0017        [ 6.2  ,  1.619],
0018        [ 6.526,  1.618],
0019        [ 6.889,  1.527],
0020        [ 7.294,  1.554],
0021        [ 7.75 ,  1.793],
0022        [ 8.267,  1.783],
0023        [ 8.857,  1.664],
0024        [ 9.538,  1.554],
0025        [10.33 ,  1.454],
0026        [15.5  ,  1.454]
0027       ])
0028
0029
0030 def divide_bins( e, mul ):
0031     """
0032     :param e: 1d array of monotonic edges
0033     :param mul: integer multiplicity with which to divide edges, eg 2 or 3 splits each edge into 2 or 3 etc..
0034     :return ee: array with extra edges obtained by splitting the input edges
0035
0036
0037          +--------+--------+     3 values, 2 bins
0038
0039          +----+---+---+----+     5 values, 4 bins    (mul 2)
0040
0041          +--+-+-+-+-+-+--+-+     9 values, 8 bins
0042
0043
0044     """
0045     nb = len(e)-1           # number of bins is one less than number of values
0046     nbb = nb*mul            # bins multiply
0047     ee = np.zeros(nbb+1)    # add one to to give number of values
0048
0049     print(" divide_bins mul %d len(e) %d len(ee) %d " % ( mul, len(e), len(ee) ))
0050
0051     for i in range(len(e)-1):
0052         a = np.linspace( e[i], e[i+1], 1+mul )
0053         i0 = i*mul
0054         i1 = (i+1)*mul
0055         print( " %2d %7.4f %7.4f i0 %2d i1 %2d     %s " % (i, e[i], e[i+1], i0, i1, a ) )
0056         ee[i0:i1+1] = a
0057     pass
0058     if mul == 1: assert np.all( ee == e )
0059
0060     non_monotonic = np.any(ee[:-1] > ee[1:])
0061     if non_monotonic:
0062         log.error("non_monotonic")
0063         print(ee)
0064     pass
0065     return ee
0066
0067
0068
0069 def edges_plot(e):
0070     fig, ax = plt.subplots(figsize=[12.8, 7.2])
0071     for m in range(1,10):
0072         ee = divide_bins(e, mul=m)
0073         yy = np.repeat( m, len(ee) )
0074         ax.scatter( ee, yy, label="mul %d " % m )
0075     pass
0076     ax.legend()
0077     fig.show()
0078
0079
0080
0081 if __name__ == '__main__':
0082     logging.basicConfig(level=logging.INFO)
0083     e = ri[:,0]
0084     #ee = divide_bins(e, mul=2 )
0085     edges_plot(e)
0086
0087
0088