reflectometry.reduction.data

1 # This program is public domain 2 """ 3 Data was the initial foray into a data representation for reflectometry, 4 and in particular for polarization correction, which was the first task. 5 6 It's main attributes is that it defines standard names for the data 7 axes x,y,z and values v, as well as a holder for labels and units. 8 9 Many of the ideas of data are preserved and live on in refldata. At 10 some point refldata may be made to inherit from data, but for data 11 will act as a platform for developing plotting ideas, and is not part 12 of the main reflectometry code. 13 """ 14 import numpy 15

16 -def edges_from_centers(x):

17 """ 18 Given a series of bin centers, compute the bin edges corresponding 19 to those centers assuming the edeges are mid-way between the centers. 20 Note: assumes the bins are the same width. If this is not the case, 21 there is insufficient information to uniquely reconstruct the edges 22 given the centers. 23 """ 24 z = numpy.zeros(len(x)+1) 25 z[1:-1] = centers_from_edges(x) 26 z[0] = x[0] - (z[1]-x[0]) 27 z[-1] = x[-1] + (x[-1]-z[-2]) 28 return z

29

30 -def centers_from_edges(x):

31 """ 32 Given a series of bin edges, return the bin centers. 33 """ 34 return (x[:-1]+x[1:])/2

35

36 -def dims(a):

37 """ 38 Printable form for matrix dimensions. 39 """ 40 return "x".join([str(v) for v in a.shape])

41

42 -def dict2text(d,indent=0):

43 """ 44 Display the contents of a dictionary with some formatting. 45 """ 46 keys = d.keys() 47 keys.sort() 48 lines = [] 49 for key in keys: 50 val = d[key] 51 if isinstance(val,dict): 52 text = "\n"+dict2text(val,indent+2) 53 elif isinstance(val,numpy.ndarray): 54 text = "array size " + dims(val) 55 else: 56 text = str(val) 57 lines.append(" "*indent + key + ": " + text) 58 return "\n".join(lines)

59 60

61 -class Data(object):

62 """ 63 Data object. 64 65 For one-dimensional data use x,v. 66 For two and three dimensional data use x,y,z,v. 67 68 Multi-dimensional data can be regular or irregular. If it is irregular, 69 the size of x,y,z must match the size of z. If it is regular, 70 the size of each vector x,y,z must match the corresponding dimension of v. 71 72 Note: for reflectivity the natural dimensions are x='Qz',y='Qx',z='Qy' 73 """ 74 _x,_xedges,dx = None,None,None 75 _y,_yedges,dy = None,None,None 76 _z,_zedges,dz = None,None,None 77 v,variance = None,None 78 xlabel,xunits='x',None 79 ylabel,yunits='y',None 80 zlabel,zunits='z',None 81 vlabel,vunits='v',None 82 83 # Store variance but allow 1-sigma uncertainty interface

84 - def _getdv(self): return numpy.sqrt(self.variance)

85 - def _setdv(self,dv): self.variance = dv**2

86 dv = property(_getdv,_setdv,'1-sigma uncertainty') 87 88 # Store centers and retrieve edges

89 - def _getx(self): return self._x

90 - def _setx(self,t):

91 self._x = t 92 self._xedges = edges_from_centers(t)

93 x = property(_getx,_setx,'x centers')

94 - def _gety(self): return self._y

95 - def _sety(self,t):

96 self._y = t 97 self._yedges = edges_from_centers(t)

98 y = property(_gety,_sety,'y centers')

99 - def _getz(self): return self._x

100 - def _setz(self,t):

101 self._z = t 102 self._zedges = edges_from_centers(t)

103 z = property(_getz,_setz,'z centers') 104 105 106 # Store edges and retrieve centers

107 - def _getxedges(self): return self._xedges

108 - def _setxedges(self,t):

109 self._x = centers_from_edges(t) 110 self._xedges = t

111 xedges = property(_getxedges,_setxedges,'x edges')

112 - def _getyedges(self): return self._yedges

113 - def _setyedges(self,t):

114 self._y = centers_from_edges(t) 115 self._yedges = t

116 yedges = property(_getyedges,_setyedges,'y edges')

117 - def _getzedges(self): return self._zedges

118 - def _setzedges(self,t):

119 self._z = centers_from_edges(t) 120 self._zedges = t

121 zedges = property(_getzedges,_setzedges,'z edges') 122 123 # Set attributes on initialization

124 - def __init__(self,**kw):

125 self.messages = [] 126 for k,v in kw.iteritems(): setattr(self,k,v)

127

128 - def summary(self):

129 """Return a text description of the contents of data""" 130 return dict2text(self.prop.__dict__)

131

132 - def log(self,msg):

133 """Record corrections that have been applied to the data""" 134 self.messages.append(msg)

135

136 - def __str__(self):

137 return "Data(%s)"%(dims(self.v))

138

139 -class PolarizedData(object):

140 """ 141 Polarized data object. This holds all four cross sections of a 142 polarized measurement, and provides operations for aligning and 143 transforming the cross sections. 144 """ 145 146 xlabel,xunits='Qz','inv A' 147 ylabel,yunits='Qx','inv A' 148 zlabel,zunits='Qy','inv A' 149 vlabel,vunits='Reflectivity',None 150

151 - def __init__(self, **kw):

152 self.pp,self.pm,self.mp,self.mm = Data(),Data(),Data(),Data() 153 self.set(**kw) 154 self.messages = []

155

156 - def set(self,**kw):

157 """Set a number of attributes simultaneously""" 158 for k,v in kw.iteritems(): 159 if hasattr(self,k): setattr(self,k,v)

160

161 - def __str__(self):

162 return "\n".join(["++"+str(self.pp),"--"+str(self.mm), 163 "+-"+str(self.pm),"-+"+str(self.mp)])

164

165 - def ispolarized(self):

166 """Indicates that the data is polarized data.""" 167 return True

168

169 - def isaligned(self):

170 """Test if all four cross sections have the same Q values.""" 171 # TODO: perform the check (monoref only) 172 return True

173

174 - def log(self,msg):

175 """Record corrections that have been applied to the data""" 176 self.messages.append(msg)

177

178 - def spin_asymmetry(self):

179 """ 180 Return the spin asymmetry for the pp and mm crosssections. 181 This is a Data object with one data line. 182 183 TODO: support multidimensional data 184 """ 185 if self.isaligned(): 186 pp, Vpp = self.pp.v, self.pp.variance 187 mm, Vmm = self.mm.v, self.mm.variance 188 else: 189 # TODO: implement matching and interpolation 190 x = match_ordinal(self.pp,self.mm) 191 pp, Vpp = interp(x,self.pp) 192 mm, Vmm = interp(x,self.mm) 193 v = (pp-mm)/(pp+mm) 194 V = v**2 * ( (1/(pp-mm) - 1/(pp+mm))**2 * Vpp 195 + (1/(pp-mm) + 1/(pp+mm))**2 * Vmm ) 196 result = Data(x,v,variance=V, 197 xlabel=self.xlabel, xunits=self.xunits, 198 vlabel="Spin asymmetry", vunits=None) 199 return result

200

201 -def refl(n=100,noise=0.02):

202 """ 203 Example 1D data: fresnel reflectivity for neutrons off silicon 204 n = number of data points 205 noise = relative error in simulated data points 206 """ 207 Q = numpy.linspace(0.,0.5,n) 208 dQ = Q*0.001 209 f = numpy.sqrt(Q**2 - 16*numpy.pi*2.07e-6 + 0j) 210 R = numpy.abs( (Q-f)/(Q+f) )**2 211 dR = noise*R 212 R = R + dR*numpy.random.randn(n) 213 data = Data(x=Q,dx=dQ,v=R,dv=dR,xlabel='Q',xunits='inv A',vlabel='R') 214 return data

215

216 -def peaks(n=40,noise=0.02):

217 """Example 2D data: peaks function""" 218 nr,nc = n,n+5 219 x = numpy.linspace(-3,3,nr) 220 y = numpy.linspace(-3,3,nc) 221 [X,Y] = numpy.meshgrid(x,y) 222 v = 3*(1-X)**2*numpy.exp(-X**2 - (Y+1)**2) \ 223 - 10*(X/5 - X**3 - Y**5)*numpy.exp(-X**2-Y**2) \ 224 - 1/3*numpy.exp(-(X+1)**2 - Y**2) 225 v = v + noise*numpy.random.randn(nc,nr) 226 data = Data(x=x,y=y,v=v) 227 return data

228

229 -def noise2d(n=40,noise=0.02,bkg=1e-10):

230 """Example 2D data: positive noise""" 231 nr,nc = n,n+1 232 x = N.linspace(0,1,nr) 233 y = N.linspace(0,1,nc) 234 v = N.abs(noise*N.random.randn(nc,nr))+bkg 235 data = Data(x=x,y=y,v=v) 236 return data

237

238 -def noisepol2d(n=40):

239 """Example 2D polarized data: positive noise""" 240 pp = noise2d(n,noise=0.05) 241 pm = noise2d(n,noise=0.3); pm.v *= 0.2; 242 mp = noise2d(n,noise=0.3); mp.v *= 0.2; 243 mm = noise2d(n,noise=0.05) 244 data = PolarizedData() 245 data.pp,data.mm = pp,mm 246 data.pm,data.mp = pm,mp 247 return data

248

249 -def peakspol(n=40):

250 """Example 2D polarized data: peaks function""" 251 pp = peaks(n,noise=0.05) 252 pm = peaks(n,noise=0.3); pm.v *= 0.4; 253 mp = peaks(n,noise=0.3); mp.v *= 0.4; 254 mm = peaks(n,noise=0.05) 255 data = PolarizedData() 256 data.pp,data.mm = pp,mm 257 data.pm,data.mp = pm,mp 258 return data

259

260 -def demo():

261 """Smoke test demo""" 262 d = noisepol2d(n=3) 263 print "polarized 2D noise:",d 264 print "++",d.pp,d.pp.v 265 print "--",d.mm,d.mm.v 266 print "+-",d.pm,d.pm.v 267 print "-+",d.mp,d.mp.v 268 print "x,y",d.mp.x,d.mp.y 269 print "edges x,y",d.mp.xedges,d.mp.yedges

270 271 if __name__ == "__main__": 272 #demo() 273 print refl() 274

Source Code for Module reflectometry.reduction.data