Changeset 69ec80f in sasmodels
- Timestamp:
- Jan 28, 2016 5:42:09 PM (9 years ago)
- Branches:
- master, core_shell_microgels, costrafo411, magnetic_model, release_v0.94, release_v0.95, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests
- Children:
- 841753c
- Parents:
- a3e78c3
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
sasmodels/data.py
rd15a908 r69ec80f 90 90 self.x, self.y, self.dx, self.dy = x, y, dx, dy 91 91 self.dxl = None 92 self.filename = None 93 self.qmin = x.min() if x is not None else np.NaN 94 self.qmax = x.max() if x is not None else np.NaN 95 self.mask = np.isnan(y) if y is not None else None 96 self._xaxis, self._xunit = "x", "" 97 self._yaxis, self._yunit = "y", "" 92 98 93 99 def xaxis(self, label, unit): … … 108 114 109 115 class Data2D(object): 110 def __init__(self): 116 def __init__(self, x=None, y=None, z=None, dx=None, dy=None, dz=None): 117 self.qx_data, self.dqx_data = x, dx 118 self.qy_data, self.dqy_data = y, dy 119 self.data, self.err_data = z, dz 120 self.mask = ~np.isnan(z) if z is not None else None 121 self.q_data = np.sqrt(x**2 + y**2) 122 self.qmin = 1e-16 123 self.qmax = np.inf 111 124 self.detector = [] 112 125 self.source = Source() 126 self.Q_unit = "1/A" 127 self.I_unit = "1/cm" 128 self.xaxis("Q_x", "A^{-1}") 129 self.yaxis("Q_y", "A^{-1}") 130 self.zaxis("Intensity", r"\text{cm}^{-1}") 131 self._xaxis, self._xunit = "x", "" 132 self._yaxis, self._yunit = "y", "" 133 self._zaxis, self._zunit = "z", "" 134 self.x_bins, self.y_bins = None, None 113 135 114 136 def xaxis(self, label, unit): … … 139 161 140 162 class Detector(object): 163 """ 164 Detector attributes. 165 """ 166 def __init__(self, pixel_size=(None, None), distance=None): 167 self.pixel_size = Vector(*pixel_size) 168 self.distance = distance 169 170 class Source(object): 171 """ 172 Beam attributes. 173 """ 141 174 def __init__(self): 142 self.pixel_size = Vector() 143 144 class Source(object): 145 pass 175 self.wavelength = np.NaN 176 self.wavelength_unit = "A" 146 177 147 178 … … 158 189 data = Data1D(q, Iq, dx=resolution * q, dy=dIq) 159 190 data.filename = "fake data" 160 data.qmin, data.qmax = q.min(), q.max()161 data.mask = np.zeros(len(q), dtype='bool')162 191 return data 163 192 … … 173 202 if qy is None: 174 203 qy = qx 204 # 5% dQ/Q resolution 175 205 Qx, Qy = np.meshgrid(qx, qy) 176 206 Qx, Qy = Qx.flatten(), Qy.flatten() 177 207 Iq = 100 * np.ones_like(Qx) 178 208 dIq = np.sqrt(Iq) 179 mask = np.ones(len(Iq), dtype='bool')180 181 data = Data2D()182 data.filename = "fake data"183 data.qx_data = Qx184 data.qy_data = Qy185 data.data = Iq186 data.err_data = dIq187 data.mask = mask188 data.qmin = 1e-16189 data.qmax = np.inf190 191 # 5% dQ/Q resolution192 209 if resolution != 0: 193 210 # https://www.ncnr.nist.gov/staff/hammouda/distance_learning/chapter_15.pdf … … 197 214 # radial (which instead it should be inverse). 198 215 Q = np.sqrt(Qx**2 + Qy**2) 199 d ata.dqx_data= resolution * Q200 d ata.dqy_data= resolution * Q216 dqx = resolution * Q 217 dqy = resolution * Q 201 218 else: 202 data.dqx_data = data.dqy_data = None 203 204 detector = Detector() 205 detector.pixel_size.x = 5 # mm 206 detector.pixel_size.y = 5 # mm 207 detector.distance = 4 # m 208 data.detector.append(detector) 219 dqx = dqy = None 220 221 data = Data2D(x=Qx, y=Qy, z=Iq, dx=dqx, dy=dqy, dz=dIq) 209 222 data.x_bins = qx 210 223 data.y_bins = qy 224 data.filename = "fake data" 225 226 # pixel_size in mm, distance in m 227 detector = Detector(pixel_size=(5, 5), distance=4) 228 data.detector.append(detector) 211 229 data.source.wavelength = 5 # angstroms 212 230 data.source.wavelength_unit = "A" 213 data.Q_unit = "1/A"214 data.I_unit = "1/cm"215 data.q_data = np.sqrt(Qx ** 2 + Qy ** 2)216 data.xaxis("Q_x", "A^{-1}")217 data.yaxis("Q_y", "A^{-1}")218 data.zaxis("Intensity", r"\text{cm}^{-1}")219 231 return data 220 232 … … 228 240 # do not repeat. 229 241 if hasattr(data, 'lam'): 230 _plot_result_sesans(data, None, None, plot_data=True, limits=limits)242 _plot_result_sesans(data, None, None, use_data=True, limits=limits) 231 243 elif hasattr(data, 'qx_data'): 232 _plot_result2D(data, None, None, view, plot_data=True, limits=limits)244 _plot_result2D(data, None, None, view, use_data=True, limits=limits) 233 245 else: 234 _plot_result1D(data, None, None, view, plot_data=True, limits=limits)246 _plot_result1D(data, None, None, view, use_data=True, limits=limits) 235 247 236 248 237 249 def plot_theory(data, theory, resid=None, view='log', 238 plot_data=True, limits=None):250 use_data=True, limits=None): 239 251 if hasattr(data, 'lam'): 240 _plot_result_sesans(data, theory, resid, plot_data=True, limits=limits)252 _plot_result_sesans(data, theory, resid, use_data=True, limits=limits) 241 253 elif hasattr(data, 'qx_data'): 242 _plot_result2D(data, theory, resid, view, plot_data, limits=limits)254 _plot_result2D(data, theory, resid, view, use_data, limits=limits) 243 255 else: 244 _plot_result1D(data, theory, resid, view, plot_data, limits=limits)256 _plot_result1D(data, theory, resid, view, use_data, limits=limits) 245 257 246 258 … … 251 263 except: 252 264 traceback.print_exc() 253 pass254 265 255 266 return wrapper … … 257 268 258 269 @protect 259 def _plot_result1D(data, theory, resid, view, plot_data, limits=None):270 def _plot_result1D(data, theory, resid, view, use_data, limits=None): 260 271 """ 261 272 Plot the data and residuals for 1D data. … … 264 275 from numpy.ma import masked_array, masked 265 276 266 plot_theory = theory is not None267 plot_resid = residis not None268 269 if data.y is None:270 plot_data = False 277 use_data = use_data and data.y is not None 278 use_theory = theory is not None 279 use_resid = resid is not None 280 num_plots = (use_data or use_theory) + use_resid 281 271 282 scale = data.x**4 if view == 'q4' else 1.0 272 283 273 if plot_data or plot_theory: 274 if plot_resid: 275 plt.subplot(121) 276 284 if use_data or use_theory: 277 285 #print(vmin, vmax) 278 286 all_positive = True 279 287 some_present = False 280 if plot_data:288 if use_data: 281 289 mdata = masked_array(data.y, data.mask.copy()) 282 290 mdata[~np.isfinite(mdata)] = masked … … 288 296 289 297 290 if plot_theory:298 if use_theory: 291 299 mtheory = masked_array(theory, data.mask.copy()) 292 300 mtheory[~np.isfinite(mtheory)] = masked … … 299 307 if limits is not None: 300 308 plt.ylim(*limits) 309 310 if num_plots > 1: 311 plt.subplot(1, num_plots, 1) 301 312 plt.xscale('linear' if not some_present else view) 302 313 plt.yscale('linear' … … 306 317 plt.ylabel('$I(q)$') 307 318 308 if plot_resid: 309 if plot_data or plot_theory: 310 plt.subplot(122) 311 319 if use_resid: 312 320 mresid = masked_array(resid, data.mask.copy()) 313 321 mresid[~np.isfinite(mresid)] = masked 314 322 some_present = (mresid.count() > 0) 323 324 if num_plots > 1: 325 plt.subplot(1, num_plots, (use_data or use_theory) + 1) 315 326 plt.plot(data.x/10, mresid, '-') 316 327 plt.xlabel("$q$/nm$^{-1}$") … … 320 331 321 332 @protect 322 def _plot_result_sesans(data, theory, resid, plot_data, limits=None):333 def _plot_result_sesans(data, theory, resid, use_data, limits=None): 323 334 import matplotlib.pyplot as plt 324 if data.y is None:325 plot_data = False326 plot_theory = theoryis not None327 plot_resid = resid is not None328 329 if plot_data or plot_theory:330 if plot_resid:331 plt.subplot(1 21)332 if plot_data:335 use_data = use_data and data.y is not None 336 use_theory = theory is not None 337 use_resid = resid is not None 338 num_plots = (use_data or use_theory) + use_resid 339 340 if use_data or use_theory: 341 if num_plots > 1: 342 plt.subplot(1, num_plots, 1) 343 if use_data: 333 344 plt.errorbar(data.x, data.y, yerr=data.dy) 334 345 if theory is not None: … … 340 351 341 352 if resid is not None: 342 if plot_data or plot_theory: 343 plt.subplot(122) 344 353 if num_plots > 1: 354 plt.subplot(1, num_plots, (use_data or use_theory) + 1) 345 355 plt.plot(data.x, resid, 'x') 346 356 plt.xlabel('spin echo length (nm)') … … 349 359 350 360 @protect 351 def _plot_result2D(data, theory, resid, view, plot_data, limits=None):361 def _plot_result2D(data, theory, resid, view, use_data, limits=None): 352 362 """ 353 363 Plot the data and residuals for 2D data. 354 364 """ 355 365 import matplotlib.pyplot as plt 356 if data.data is None:357 plot_data = False358 plot_theory = theoryis not None359 plot_resid = resid is not None366 use_data = use_data and data.data is not None 367 use_theory = theory is not None 368 use_resid = resid is not None 369 num_plots = use_data + use_theory + use_resid 360 370 361 371 # Put theory and data on a common colormap scale 362 if limits is None: 363 vmin, vmax = np.inf, -np.inf 364 if plot_data: 365 target = data.data[~data.mask] 366 datamin = target[target>0].min() if view == 'log' else target.min() 367 datamax = target.max() 368 vmin = min(vmin, datamin) 369 vmax = max(vmax, datamax) 370 if plot_theory: 371 theorymin = theory[theory>0].min() if view=='log' else theory.min() 372 theorymax = theory.max() 373 vmin = min(vmin, theorymin) 374 vmax = max(vmax, theorymax) 375 else: 372 vmin, vmax = np.inf, -np.inf 373 if use_data: 374 target = data.data[~data.mask] 375 datamin = target[target > 0].min() if view == 'log' else target.min() 376 datamax = target.max() 377 vmin = min(vmin, datamin) 378 vmax = max(vmax, datamax) 379 if use_theory: 380 theorymin = theory[theory > 0].min() if view == 'log' else theory.min() 381 theorymax = theory.max() 382 vmin = min(vmin, theorymin) 383 vmax = max(vmax, theorymax) 384 385 # Override data limits from the caller 386 if limits is not None: 376 387 vmin, vmax = limits 377 388 378 if plot_data: 379 if plot_theory and plot_resid: 380 plt.subplot(131) 381 elif plot_theory or plot_resid: 382 plt.subplot(121) 389 # Plot data 390 if use_data: 391 if num_plots > 1: 392 plt.subplot(1, num_plots, 1) 383 393 _plot_2d_signal(data, target, view=view, vmin=vmin, vmax=vmax) 384 394 plt.title('data') … … 386 396 h.set_label('$I(q)$') 387 397 388 if plot_theory: 389 if plot_data and plot_resid: 390 plt.subplot(132) 391 elif plot_data: 392 plt.subplot(122) 393 elif plot_resid: 394 plt.subplot(121) 398 # plot theory 399 if use_theory: 400 if num_plots > 1: 401 plt.subplot(1, num_plots, use_data+1) 395 402 _plot_2d_signal(data, theory, view=view, vmin=vmin, vmax=vmax) 396 403 plt.title('theory') … … 400 407 else '$I(q)$') 401 408 402 #if plot_data or plot_theory: 403 # plt.colorbar() 404 405 if plot_resid: 406 if plot_data and plot_theory: 407 plt.subplot(133) 408 elif plot_data or plot_theory: 409 plt.subplot(122) 409 # plot resid 410 if use_resid: 411 if num_plots > 1: 412 plt.subplot(1, num_plots, use_data+use_theory+1) 410 413 _plot_2d_signal(data, resid, view='linear') 411 414 plt.title('residuals')
Note: See TracChangeset
for help on using the changeset viewer.