reflectometry.reduction.rebin

1 """ 2 1-D and 2-D rebinning code. 3 """ 4 import numpy 5 6 from reflectometry.reduction import _reduction 7

8 -def rebin(x,I,xo,Io=None,dtype=numpy.float64):

9 """ 10 Rebin a vector. 11 12 x are the existing bin edges 13 xo are the new bin edges 14 I are the existing counts (one fewer than edges) 15 16 Io will be used if present, but be sure that it is a contiguous 17 array of the correct shape and size. 18 19 dtype is the type to use for the intensity vectors. This can be 20 integer (uint8, uint16, uint32) or real (float32 or f, float64 or d). 21 The edge vectors are all coerced to doubles. 22 23 Note that total intensity is not preserved for integer rebinning. 24 The algorithm uses truncation so total intensity will be down on 25 average by half the total number of bins. 26 """ 27 # Coerce axes to float arrays 28 x,xo = _input(x,dtype='d'),_input(xo,dtype='d') 29 shape_in = numpy.array([x.shape[0]-1]) 30 shape_out = numpy.array([xo.shape[0]-1]) 31 32 # Coerce counts to correct type and check shape 33 if dtype is None: 34 try: 35 dtype = I.dtype 36 except AttributeError: 37 dtype = numpy.float64 38 I = _input(I, dtype=dtype) 39 if shape_in != I.shape: 40 raise TypeError("input array incorrect shape %s"%I.shape) 41 42 # Create output vector 43 Io = _output(Io,shape_out,dtype=dtype) 44 45 # Call rebin on type if it is available 46 try: 47 rebincore = getattr(_reduction,'rebin_'+I.dtype.name) 48 except AttributeError: 49 raise TypeError("rebin supports uint8 uint16 uint32 float32 float64, not " 50 + I.dtype.name) 51 rebincore(x,I,xo,Io) 52 return Io

53

54 -def rebin2d(x,y,I,xo,yo,Io=None,dtype=None):

55 """ 56 Rebin a matrix. 57 58 x,y are the existing bin edges 59 xo,yo are the new bin edges 60 I is the existing counts (one fewer than edges in each direction) 61 62 For example, with x representing the column edges in each row and 63 y representing the row edges in each column, the following 64 represents a uniform field:: 65 66 >>> from reflectometry.reduction.rebin import rebin2d 67 >>> x,y = [0,2,4,5], [0,1,3] 68 >>> z = [[2,2,1],[4,4,2]] 69 70 We can check this by rebinning with uniform size bins:: 71 72 >>> xo,yo = range(6), range(4) 73 >>> rebin2d(y,x,z,yo,xo) 74 array([[ 1., 1., 1., 1., 1.], 75 [ 1., 1., 1., 1., 1.], 76 [ 1., 1., 1., 1., 1.]]) 77 78 dtype is the type to use for the intensity vectors. This can be 79 integer (uint8, uint16, uint32) or real (float32 or f, float64 or d). 80 The edge vectors are all coerced to doubles. 81 82 Note that total intensity is not preserved for integer rebinning. 83 The algorithm uses truncation so total intensity will be down on 84 average by half the total number of bins. 85 86 Io will be used if present, if it is contiguous and if it has the 87 correct shape and type for the input. Otherwise it will raise a 88 TypeError. This will allow you to rebin the slices of an appropriately 89 ordered matrix without making copies. 90 """ 91 # Coerce axes to float arrays 92 x,y,xo,yo = [_input(v, dtype='d') for v in (x,y,xo,yo)] 93 shape_in = numpy.array([x.shape[0]-1,y.shape[0]-1]) 94 shape_out = numpy.array([xo.shape[0]-1,yo.shape[0]-1]) 95 96 # Coerce counts to correct type and check shape 97 if dtype is None: 98 try: dtype = I.dtype 99 except AttributeError: dtype = numpy.float64 100 I = _input(I, dtype=dtype) 101 if (shape_in != I.shape).any(): 102 raise TypeError("input array incorrect shape %s"%str(I.shape)) 103 104 # Create output vector 105 Io = _output(Io,shape_out,dtype=dtype) 106 107 # Call rebin on type if it is available 108 try: 109 rebincore = getattr(_reduction,'rebin2d_'+I.dtype.name) 110 except AttributeError: 111 raise TypeError("rebin2d supports uint8 uint16 uint32 float32 float64, not " 112 + I.dtype.name) 113 rebincore(x,y,I,xo,yo,Io) 114 return Io

115

116 -def _input(v, dtype='d'):

117 """ 118 Force v to be a contiguous array of the correct type, avoiding copies 119 if possible. 120 """ 121 return numpy.ascontiguousarray(v,dtype=dtype)

122

123 -def _output(v, shape, dtype=numpy.float64):

124 """ 125 Create a contiguous array of the correct shape and type to hold a 126 returned array, reusing an existing array if possible. 127 """ 128 if v is None: 129 return numpy.empty(shape,dtype=dtype) 130 if not (isinstance(v,numpy.ndarray) 131 and v.dtype == numpy.dtype(dtype) 132 and (v.shape == shape).all() 133 and v.flags.contiguous): 134 raise TypeError("output vector must be contiguous %s of size %s" 135 %(dtype,shape)) 136 return v

137 138 139 # ================ Test code ================== 140 # TODO: move test code to its own file

141 -def _check1d(from_bins,val,to_bins,target):

142 target = _input(target) 143 for (f,F) in [(from_bins,val), (from_bins[::-1],val[::-1])]: 144 for (t,T) in [(to_bins,target), (to_bins[::-1],target[::-1])]: 145 result = rebin(f,F,t) 146 assert numpy.linalg.norm(T-result) < 1e-14, \ 147 "rebin failed for %s->%s %s"%(f,t,result)

148

149 -def _test1d():

150 # Split a value 151 _check1d([1,2,3,4],[10,20,30],[1,2.5,4],[20,40]) 152 153 # bin is a superset of rebin 154 _check1d([0,1,2,3,4],[5,10,20,30],[1,2.5,3],[20,10]); 155 156 # bin is a subset of rebin 157 _check1d([ 1, 2, 3, 4, 5, 6], 158 [ 10, 20, 30, 40, 50], 159 [ 2.5, 3.5], [25]) 160 161 # one bin to many 162 _check1d([1, 2, 3, 4, 5, 6], 163 [10, 20, 30, 40, 50], 164 [2.1, 2.2, 2.3, 2.4 ], 165 [2, 2, 2]); 166 167 # many bins to one 168 _check1d([1, 2, 3, 4, 5, 6], 169 [10, 20, 30, 40, 50], 170 [ 2.5, 4.5 ], 171 [60])

172

173 -def _check2d(x,y,z,xo,yo,zo):

174 result = rebin2d(x,y,z,xo,yo) 175 target = numpy.array(zo,dtype=result.dtype) 176 assert numpy.linalg.norm(target-result) < 1e-14, \ 177 "rebin2d failed for %s,%s->%s,%s\n%s\n%s\n%s"%(x,y,z,xo,yo,zo)

178

179 -def _uniform_test(x,y):

180 z = numpy.array([y],'d') * numpy.array([x],'d').T 181 xedges = numpy.concatenate([(0,),numpy.cumsum(x)]) 182 yedges = numpy.concatenate([(0,),numpy.cumsum(y)]) 183 nx = numpy.round(xedges[-1]) 184 ny = numpy.round(yedges[-1]) 185 ox = numpy.arange(nx+1) 186 oy = numpy.arange(ny+1) 187 target = numpy.ones([nx,ny],'d') 188 _check2d(xedges,yedges,z,ox,oy,target)

189

190 -def _test2d():

191 x,y,I = [0,3,5,7], [0,1,3], [[3,6],[2,4],[2,4]] 192 xo,yo,Io = range(8), range(4), [[1]*3]*7 193 x,y,I,xo,yo,Io = [numpy.array(A,'d') for A in [x,y,I,xo,yo,Io]] 194 195 # Try various types and orders on a non-square matrix 196 _check2d(x,y,I,xo,yo,Io) 197 _check2d(x[::-1],y,I[::-1,:],xo,yo,Io) 198 _check2d(x,y[::-1],I[:,::-1],xo,yo,Io) 199 _check2d(x,y,I,[7,3,0],yo,[[4]*3,[3]*3]) 200 _check2d(x,y,I,xo,[3,2,0],[[1,2]]*7) 201 _check2d(y,x,I.T,yo,xo,Io.T) # C vs. Fortran ordering 202 203 # Test smallest possible result 204 _check2d([-1,2,4], [0,1,3], [[3,6],[2,4]], 205 [1,2], [1,2], [1]) 206 # subset/superset 207 _check2d([0,1,2,3], [0,1,2,3], [[1]*3]*3, 208 [0.5,1.5,2.5], [0.5,1.5,2.5], [[1]*2]*2) 209 for dtype in ['uint8','uint16','uint32','float32','float64']: 210 _check2d([0,1,2,3,4], [0,1,2,3,4], 211 numpy.array([[1]*4]*4, dtype=dtype), 212 [-2,-1,2,5,6], [-2,-1,2,5,6], 213 numpy.array([[0,0,0,0],[0,4,4,0],[0,4,4,0],[0,0,0,0]], 214 dtype=dtype) 215 ) 216 # non-square test 217 _uniform_test([1,2.5,4,0.5],[3,1,2.5,1,3.5]) 218 _uniform_test([3,2],[1,2])

219

220 -def test():

221 _test1d() 222 _test2d()

223 224 if __name__ == "__main__": test() 225

Source Code for Module reflectometry.reduction.rebin