# Ticket #312: slit_test.py

File slit_test.py, 18.7 KB (added by pkienzle, 9 years ago) |
---|

Line | |
---|---|

1 | r""" |

2 | Explore different integration schemes for slit smearing. |

3 | |

4 | Introduction |

5 | ============ |

6 | |

7 | Slit smearing evaluates the following integral: |

8 | |

9 | .. math:: |

10 | |

11 | I_s(q_o) = \int_0^{\Delta q_v} I(\sqrt{q_o^2 + u^2}) du |

12 | |

13 | where $\Delta q_v$ is the vertical slit contribution. See |

14 | `http://www.ncnr.nist.gov/staff/hammouda/distance_learning/chapter_61.pdf`_ |

15 | |

16 | Procedure |

17 | ========= |

18 | |

19 | This program is not intended as an end user application, but rather as a |

20 | platform for exploring implementations of the resolution integral. |

21 | Different experiments can be performed by modifying the :func:`demo` function. |

22 | The integration functions have the name smear_METHOD, for a variety of |

23 | different methods. Each one takes a form factor with parameters, a set of |

24 | q values, a resolution, and maybe a set of q_calc values at which to evaluate |

25 | the theory function. The *ms* and *mt* variables in :func:`demo`, determine |

26 | which methods are being compared. :func:`interpolate` allows you to fill |

27 | in the integral between q points since the demo function is rapidly |

28 | oscillating. *TEST_DATA_SLIT_SMEAR* contains the values calculated in |

29 | Igor for a sphere model with a radius of 6 microns and a contrast of 3e-6 |

30 | inverse Angstroms. |

31 | |

32 | :func:`smear_romberg` is a direct implementation of the integral using |

33 | adaptive romberg integration from scipy. This method would add considerable |

34 | complexity to sasmodels since we could no longer rely on evaluating the |

35 | models at a fixed set of q points for all parameter sets in a fit. A |

36 | "good enough" method with fixed q, although less accurate and requiring |

37 | more q points on average to evaluate, will be a win in complexity and perhaps |

38 | performance. The romberg integration parameters are set to produce a |

39 | pretty good value for the integral. |

40 | |

41 | The remaining methods used a fixed set of *q* points, *q_calc*, to evaluate the |

42 | integral. They are :func:`smear_simpsons`, :func:`smear_trapezoid`, |

43 | and :func:`smear_rectangular`. Except for :func:`smear_simpsons_pow`, |

44 | *q_calc* is extrapolated to the full width of the resolution function. |

45 | The extrapolation uses exponential spacing between the points, much like *q* |

46 | spacing within the measured data range, though :func:`smear_simpsons_lin` |

47 | uses linear spacing. |

48 | |

49 | :func:`smear_rectangular` is implemented as a non-square weight matrix which |

50 | can be multiplied by a set of computed values *I(q_calc)* to produce |

51 | *I_smeared(q)*. The same style could be used for trapezoid and Simpson's |

52 | rule with the appropriate bookkeeping. This was not done because at least |

53 | for the demo problem, rectangular performed better than the other |

54 | non-adaptive methods, and so wins on both simplicity and accuracy. |

55 | |

56 | :func:`smear_simpsons_pow` uses simpson's rule to estimate the integral within |

57 | the range of the data, and a power law approximation for the tail above the |

58 | maximum q value. This is akin to the Igor implementation of slit smearing, |

59 | though it uses the theory function itself to estimate the power law rather |

60 | than data. |

61 | |

62 | The power law approximation has two problems: |

63 | |

64 | 1. the power is estimated from the data, which means that this will only |

65 | work if the data has a flat region at the tail, and |

66 | 2. any features such as lamellar peaks that occur beyond the end of the |

67 | data will not be correctly modeled. |

68 | |

69 | The demo shows that this method performs poorly for high q when oscillations |

70 | in the theory have not completely decayed by the end of the q range, which |

71 | is not really a surprise. It is not a problem in practice for usans data on |

72 | igor since polydispersity will likely dampen any ringing that a monodisperse |

73 | model might exhibit. The demo does not simulate a model containing a peak at |

74 | q beyond the end of the data. |

75 | |

76 | |

77 | Conclusions |

78 | =========== |

79 | |

80 | At least for this theory function, the simple rectangular integral is as |

81 | good as the more complicated trapezoid and Simpson's rule integration. |

82 | Due to the oscillations in the theory function, the relative error is |

83 | huge (50%) at high q. Oversampling by a factor of 10 (requiring 1020 |

84 | evaluations rather than 122), reduces the error to 2.5%. Another factor |

85 | of 10 (9130 evaluations) reduces error to 0.15%. Getting error down to |

86 | 5e-5 requires about 90,000 evaluations. |

87 | |

88 | In comparison, romberg integration requires only 1360 evaluations. |

89 | |

90 | """ |

91 | import numpy as np |

92 | from numpy import sqrt, cos, sin, pi, log10, log, exp |

93 | |

94 | def smear_romberg(q, form, pars, delta_qv, q_calc=None): |

95 | """ |

96 | Romberg integration. |

97 | |

98 | This is an adaptive integration technique. It is called with settings |

99 | that make it slow to evaluate but give it good accuracy. |

100 | """ |

101 | from scipy.integrate import romberg |

102 | evaluations = [0] |

103 | def _fn(u, q0): |

104 | evaluations[0] += 1 if np.isscalar(q0) else len(q0) |

105 | return form(sqrt(q0**2 + u**2), *pars) |

106 | r = [romberg(_fn, 0, delta_qv, args=(qi,), |

107 | divmax=100, vec_func=True, tol=0, rtol=1e-14) |

108 | for qi in q] |

109 | print "romberg evaluations:", evaluations[0] |

110 | return r |

111 | |

112 | |

113 | def smear_simpsons_pow(q, form, pars, delta_qv, q_calc=None): |

114 | """ |

115 | Use Simpson's rule with power low approximation to compute the integral. |

116 | |

117 | The *q* range between the end of the data and *delta_qv* is not sampled. |

118 | Instead it is approximated by a power law computed from the (I(delta_qv)) and |

119 | the final 25\% of the theory function. Igor fits the data to the a power |

120 | law and assumes that it will be an adequate estimate for the fitted |

121 | theory function. |

122 | """ |

123 | from scipy.integrate import simps |

124 | if q_calc is None: q_calc = q |

125 | x = np.hstack((q_calc, q_calc[-1]+delta_qv)) |

126 | print "#q:",q.size, "#q_calc:",x.size |

127 | y = form(x, *pars) |

128 | A, m = fit_power_law(x, y, portion=0.25) |

129 | return [simps(y[i:-1], q_to_u(x[i:-1], q[i])) |

130 | + power_law_residual(delta_qv, q[i], q[-1], A, m) |

131 | for i in range(len(q))] |

132 | |

133 | |

134 | def smear_simpsons_lin(q, form, pars, delta_qv, q_calc=None): |

135 | """ |

136 | Use Simpson's rule with linear extrapolation to compute the integral. |

137 | |

138 | The evaluation points *q* are extrapolated to *delta_qv+q[-1]* using |

139 | the final *q* step as the step size. |

140 | |

141 | For a typical NCNR USANS dataset this requires over 10000 extra evaluations |

142 | of the function. |

143 | """ |

144 | from scipy.integrate import simps |

145 | if q_calc is None: q_calc = q |

146 | x = linear_extrapolation(q_calc, q_calc[-1]+delta_qv) |

147 | print "#q:",q.size, "#q_calc:",x.size |

148 | y = form(x, *pars) |

149 | return [simps(y[j:], q_to_u(x[j:], qi)) |

150 | for qi in q for j in [np.searchsorted(x,qi)]] |

151 | |

152 | |

153 | def smear_simpsons(q, form, pars, delta_qv, q_calc=None): |

154 | """ |

155 | Use Simpson's rule with geometric extrapolation to compute the integral. |

156 | |

157 | The evaluation points *q* are extrapolated to *delta_qv+q[-1]* using |

158 | the geometric average of *q* as the step size. |

159 | |

160 | For a typical NCNR USANS dataset the function is evaluated at twice as |

161 | many points as there are points in q. |

162 | """ |

163 | from scipy.integrate import simps |

164 | if q_calc is None: q_calc = q |

165 | x = geometric_extrapolation(q_calc, q_calc[-1]+delta_qv) |

166 | print "#q:",q.size, "#q_calc:",x.size |

167 | y = form(x, *pars) |

168 | return [simps(y[j:], q_to_u(x[j:], qi)) |

169 | for qi in q for j in [np.searchsorted(x,qi)]] |

170 | |

171 | |

172 | def smear_trapezoid(q, form, pars, delta_qv, q_calc=None): |

173 | """ |

174 | Use the trapezoid rule with geometric extrapolation to compute the integral. |

175 | """ |

176 | from scipy.integrate import trapz |

177 | |

178 | if q_calc is None: q_calc = q |

179 | x = geometric_extrapolation(q_calc, q_calc[-1]+delta_qv) |

180 | print "#q:",q.size, "#q_calc:",x.size |

181 | y = form(x, *pars) |

182 | return [trapz(y[j:], q_to_u(x[j:], qi)) |

183 | for qi in q for j in [np.searchsorted(x,qi)]] |

184 | |

185 | |

186 | def smear_rectangular(q, form, pars, delta_qv, q_calc=None): |

187 | """ |

188 | Use the trapezoid rule with geometric extrapolation to compute the integral. |

189 | """ |

190 | from scipy.integrate import trapz |

191 | if q_calc is None: q_calc = q |

192 | x = geometric_extrapolation(q_calc, q_calc[-1]+delta_qv) |

193 | print "#q:",q.size, "#q_calc:",x.size |

194 | edges = bin_edges(x) |

195 | W = np.zeros((len(q), len(x)), 'd') |

196 | for i, qi in enumerate(q): |

197 | W[i,:] = np.diff(q_to_u(edges, qi)) |

198 | |

199 | y = form(x, *pars) |

200 | return np.dot(W,y) |

201 | |

202 | |

203 | ########################################################################### |

204 | # Helper functions |

205 | |

206 | def sphere(q, radius, contrast): |

207 | """ |

208 | Sphere form factor with no polydispersity |

209 | """ |

210 | qr = q * radius |

211 | sn, cn = sin(qr), cos(qr) |

212 | bes = 3 * (sn - qr * cn) / qr ** 3 # may be 0/0 but we fix that next line |

213 | #bes[qr == 0] = 1 |

214 | fq = bes * contrast |

215 | return 4*pi*radius**3/3 * 1.0e-4 * fq ** 2 |

216 | |

217 | |

218 | |

219 | def q_to_u(q, q0): |

220 | """ |

221 | Convert *q* values to *u* values for the integral computed at *q0*. |

222 | """ |

223 | # array([value])**2 - value**2 is not always zero |

224 | qpsq = q**2 - q0**2 |

225 | qpsq[qpsq<0] = 0 |

226 | return sqrt(qpsq) |

227 | |

228 | |

229 | def interpolate(q, max_step): |

230 | """ |

231 | Returns *q_calc* with points spaced at most max_step apart. |

232 | """ |

233 | step = np.diff(q) |

234 | index = step>max_step |

235 | if np.any(index): |

236 | inserts = [] |

237 | for q_i,step_i in zip(q[:-1][index],step[index]): |

238 | n = np.ceil(step_i/max_step) |

239 | inserts.extend(q_i + np.arange(1,n)*(step_i/n)) |

240 | # Extend a couple of fringes beyond the end of the data |

241 | inserts.extend(q[-1] + np.arange(1,16)*max_step) |

242 | q_calc = np.sort(np.hstack((q,inserts))) |

243 | else: |

244 | q_calc = q |

245 | return q_calc |

246 | |

247 | def linear_extrapolation(q, qmax): |

248 | """ |

249 | Extrapolate *q* out to *qmax* using the final step in *q* as the step |

250 | size. |

251 | """ |

252 | extrapolate = np.arange(q[-1], qmax, q[-1]-q[-2]) |

253 | x = np.concatenate((q, extrapolate[1:])) |

254 | |

255 | |

256 | def geometric_extrapolation(q, qmax): |

257 | r""" |

258 | Extrapolate *q* out to *qmax* using geometric steps, with the average |

259 | geometric step size in *q* as the step size. |

260 | |

261 | Starting at $q_1$ and stepping geometrically by $\Delta q$ |

262 | to $q_n$ in $n$ points gives a geometric average of: |

263 | |

264 | .. math:: |

265 | |

266 | \log \Delta q = (\log q_n - log q_1) / (n - 1) |

267 | |

268 | From this we can compute the number of steps required to extend $q$ |

269 | from $q_n$ to $q_\text{max}$ by $\Delta q$ as: |

270 | |

271 | .. math:: |

272 | |

273 | n_\text{extend} = (\log q_\text{max} - \log q_n) / \log \Delta q |

274 | |

275 | Substituting: |

276 | |

277 | n_\text{extend} = (n-1) (\log q_\text{max} - \log q_n) |

278 | / (\log q_n - log q_1) |

279 | """ |

280 | n_extend = (len(q)-1)*(log(qmax) - log(q[-1]))/(log(q[-1]) - log(q[0])) |

281 | extrapolate = np.logspace(log10(q[-1]), log10(qmax), np.ceil(n_extend)+1) |

282 | x = np.concatenate((q, extrapolate[1:])) |

283 | return x |

284 | |

285 | |

286 | def fit_power_law(x, y, portion=1.0): |

287 | """ |

288 | Fit a power law $y=Ax^{-m}$ against *x*, *y*. If *portion* is given, use |

289 | only the given portion from the tail of the data in the fit. |

290 | """ |

291 | fitting_points = max(2, np.ceil(len(x)*portion)) |

292 | p = np.polyfit(log(x[-fitting_points:]), log(y[-fitting_points:]), 1) |

293 | A,m = exp(p[1]), -p[0] |

294 | #print "A,m",A,m |

295 | |

296 | |

297 | def power_law_residual(delta_qv, q0, qmax, A, m): |

298 | r""" |

299 | Residual term for the integration, assuming I(q) follows a power law for |

300 | high q. |

301 | |

302 | This integrates $I_p(\sqrt(q_o^2 + u^2})$ by $u$ for $u$ corresponding to |

303 | $q_\text{max}$ to $u=\Delta q_v$, where $I_p$ is: |

304 | |

305 | .. math:: |

306 | |

307 | I_p(q) = A q^{-m} |

308 | """ |

309 | from scipy.integrate import romberg |

310 | a, b = sqrt(qmax**2 - q0**2), delta_qv |

311 | _residual_fn = lambda u, q0, A, m: A*(q0**2 + u**2)**(-m/2) |

312 | r = romberg(_residual_fn, a, b, args=(q0, A, m), |

313 | divmax=100, vec_func=True, tol=0, rtol=1e-8) |

314 | return r |

315 | |

316 | def bin_edges(x): |

317 | """ |

318 | Determine bin edges from bin centers, assuming that edges are centered |

319 | between the bins. |

320 | |

321 | Note: this uses the arithmetic mean, which may not be appropriate for |

322 | log-scaled data. |

323 | """ |

324 | if len(x) < 2 or (np.diff(x)<0).any(): |

325 | raise ValueError("Expected bins to be an increasing set") |

326 | edges = np.hstack([ |

327 | x[0] - 0.5*(x[1] - x[0]), # first point minus half first interval |

328 | 0.5*(x[1:] + x[:-1]), # mid points of all central intervals |

329 | x[-1] + 0.5*(x[-1] - x[-2]), # last point plus half last interval |

330 | ]) |

331 | return edges |

332 | |

333 | def Iq_smeared(q, form, pars, scale=1, background=0., delta_qv=0.117, |

334 | method='trapezoid', max_step=None): |

335 | """ |

336 | Compute slit-smeared version of *form(q;pars)* |

337 | """ |

338 | q_calc = interpolate(q, max_step) if max_step is not None else q |

339 | integrator = globals()['smear_'+method] |

340 | res = integrator(q, form, pars, delta_qv, q_calc) |

341 | return scale/delta_qv*np.asarray(res) + background |

342 | |

343 | |

344 | def Iq(q, form, pars, scale=1, background=0): |

345 | """ |

346 | Compute unsmeared version of *form(q;pars)* |

347 | """ |

348 | return scale*form(q, *pars) + background |

349 | |

350 | |

351 | TEST_DATA_SLIT_SPHERE = """\ |

352 | 2.26097e-05 0.117 5.5781372896e+09 1.4626077708e+06 |

353 | 2.53847e-05 0.117 5.0363141458e+09 1.3117318023e+06 |

354 | 2.81597e-05 0.117 4.4850108103e+09 1.1594863713e+06 |

355 | 3.09347e-05 0.117 3.9364658459e+09 1.0094881630e+06 |

356 | 3.37097e-05 0.117 3.4019975074e+09 8.6518941303e+05 |

357 | 3.92597e-05 0.117 2.4139519814e+09 6.0232158311e+05 |

358 | 4.48097e-05 0.117 1.5816877820e+09 3.8739994090e+05 |

359 | 5.03597e-05 0.117 9.3715407224e+08 2.2745304775e+05 |

360 | 5.59097e-05 0.117 4.8387917428e+08 1.2101295768e+05 |

361 | 6.14597e-05 0.117 2.0193586928e+08 6.0055107771e+04 |

362 | 6.70097e-05 0.117 5.5886110911e+07 3.2749521065e+04 |

363 | 7.25597e-05 0.117 3.7782348010e+06 2.6350963616e+04 |

364 | 7.81097e-05 0.117 5.3407817904e+06 2.9624963314e+04 |

365 | 8.36597e-05 0.117 2.7975485523e+07 3.4403962254e+04 |

366 | 8.92097e-05 0.117 4.9845448282e+07 3.6130017903e+04 |

367 | 9.47597e-05 0.117 6.0092588905e+07 3.3495107849e+04 |

368 | 1.00310e-04 0.117 5.6823430831e+07 2.7475726279e+04 |

369 | 1.05860e-04 0.117 4.3857024036e+07 2.0144282226e+04 |

370 | 1.11410e-04 0.117 2.7277144760e+07 1.3647403260e+04 |

371 | 1.22510e-04 0.117 3.3119334113e+06 6.6519711526e+03 |

372 | 1.33610e-04 0.117 1.4412859402e+06 6.9726212813e+03 |

373 | 1.44710e-04 0.117 8.5620162463e+06 8.1441335775e+03 |

374 | 1.55810e-04 0.117 9.6957429033e+06 6.4559996521e+03 |

375 | 1.66910e-04 0.117 4.3818341914e+06 3.6252154156e+03 |

376 | 1.78010e-04 0.117 2.7448997387e+05 2.4006505342e+03 |

377 | 1.89110e-04 0.117 8.0472009936e+05 2.8187789089e+03 |

378 | 2.00210e-04 0.117 2.8149552834e+06 3.0915662855e+03 |

379 | 2.11310e-04 0.117 2.7510907861e+06 2.3722530293e+03 |

380 | 2.22410e-04 0.117 1.0053133293e+06 1.4473468311e+03 |

381 | 2.33510e-04 0.117 5.8428305052e+03 1.2048540556e+03 |

382 | 2.44610e-04 0.117 5.1699305004e+05 1.4625670042e+03 |

383 | 2.55710e-04 0.117 1.2120227268e+06 1.5010705973e+03 |

384 | 2.66810e-04 0.117 9.7896842846e+05 1.1336343426e+03 |

385 | 2.77910e-04 0.117 2.5507264791e+05 8.1848818080e+02 |

386 | 3.05660e-04 0.117 5.2403101181e+05 7.4913374357e+02 |

387 | 3.33410e-04 0.117 5.8699343809e+04 4.4669964560e+02 |

388 | 3.61160e-04 0.117 3.0844327150e+05 4.6774007542e+02 |

389 | 3.88910e-04 0.117 8.3360142970e+03 2.7169550220e+02 |

390 | 4.16660e-04 0.117 1.8630080583e+05 3.0710983679e+02 |

391 | 4.44410e-04 0.117 3.1616804732e-01 1.7959006831e+02 |

392 | 4.72160e-04 0.117 1.1299016314e+05 2.0763952339e+02 |

393 | 4.99910e-04 0.117 2.9952522747e+03 1.2536542765e+02 |

394 | 5.27660e-04 0.117 6.7625695649e+04 1.4013969777e+02 |

395 | 5.55410e-04 0.117 7.6927460089e+03 8.2145593180e+01 |

396 | 6.10910e-04 0.117 1.1229057779e+04 8.4519745643e+01 |

397 | 6.66410e-04 0.117 1.3035567943e+04 8.1554625609e+01 |

398 | 7.21910e-04 0.117 1.3309931343e+04 7.4437319172e+01 |

399 | 7.77410e-04 0.117 1.2462626212e+04 6.4697088261e+01 |

400 | 8.32910e-04 0.117 1.0912927143e+04 5.3773301044e+01 |

401 | 8.88410e-04 0.117 9.0172597469e+03 4.2843375753e+01 |

402 | 9.43910e-04 0.117 7.0496495917e+03 3.2771032724e+01 |

403 | 9.99410e-04 0.117 5.2030483682e+03 2.4113557144e+01 |

404 | 1.05491e-03 0.117 3.5988976711e+03 1.7160773658e+01 |

405 | 1.11041e-03 0.117 2.2996060652e+03 1.2016626459e+01 |

406 | 1.22141e-03 0.117 6.4766590598e+02 6.0373017740e+00 |

407 | 1.33241e-03 0.117 4.1963483264e+01 4.5215452974e+00 |

408 | 1.44341e-03 0.117 6.3370708246e+01 5.1054681903e+00 |

409 | 1.55441e-03 0.117 3.0736750577e+02 5.9176165298e+00 |

410 | 1.66541e-03 0.117 5.0327682399e+02 5.9815000189e+00 |

411 | 1.77641e-03 0.117 5.4084331454e+02 5.1634639625e+00 |

412 | 1.88741e-03 0.117 4.3488671756e+02 3.8535158148e+00 |

413 | 1.99841e-03 0.117 2.6322287860e+02 2.5824997753e+00 |

414 | 2.10941e-03 0.117 1.0793633150e+02 1.7315517194e+00 |

415 | 2.22041e-03 0.117 1.8474448850e+01 1.4077213604e+00 |

416 | 2.33141e-03 0.117 1.5864062279e+00 1.4771560682e+00 |

417 | 2.44241e-03 0.117 3.2267213848e+01 1.6916253448e+00 |

418 | 2.55341e-03 0.117 7.4289116207e+01 1.8274751193e+00 |

419 | 2.66441e-03 0.117 9.9000521929e+01 1.7706812289e+00 |

420 | """ |

421 | |

422 | def demo(): |

423 | from matplotlib import pyplot as plt |

424 | |

425 | ## ==== Set data and model parameters |

426 | delta_qv = 0.117 |

427 | q = np.logspace(-5,-3,400) |

428 | #q = np.logspace(log10(2.3e-5),log10(2.7e-3),68*10) |

429 | #q = np.logspace(log10(2.3e-5),log10(2.7e-3),68) |

430 | if 1: # use values from Igor |

431 | data = np.loadtxt(TEST_DATA_SLIT_SPHERE.split('\n')).T |

432 | q, dq, _, answer = data |

433 | delta_qv = dq[0] |

434 | else: |

435 | answer = None |

436 | scale, background = 0.01, 0. #0.01 |

437 | pars = (60000, 3) # radius, contrast |

438 | |

439 | ## Set the maximum step size between q values for numerical integration |

440 | ## This can be based on the biggest dimension in the model, the visible |

441 | ## oscillation in the data, or None if the data is assumed to be j |

442 | d_max = pars[0] |

443 | points_per_fringe = 2000 |

444 | max_step = 2*pi/(d_max*points_per_fringe) |

445 | #max_step = None |

446 | #print "max_step",max_step |

447 | |

448 | ## ==== Select which models are being compared |

449 | ## Models are smear_MODEL functions above. |

450 | ## romberg is the slow accurate model computed with adaptive integration |

451 | ## rectangular is the simplest model to implement, and hopefully good |

452 | ## enough. It is the only one that uses max_step at the time of this |

453 | ## writing, and is what we will use in sasmodels since it seems to be |

454 | ## good enough. |

455 | #ms,mt = "simpsons","rectangular" |

456 | ms,mt = "romberg","rectangular" |

457 | #ms,mt = "romberg","trapezoid" |

458 | #ms,mt = "romberg","simpsons" |

459 | |

460 | ## ==== Evaluate the model and its comparison |

461 | q_calc = interpolate(q, max_step) if max_step is not None else q |

462 | y = Iq(q_calc, sphere, pars, scale=scale, background=background) |

463 | ys = Iq_smeared(q, sphere, pars, delta_qv=delta_qv, max_step=max_step, |

464 | scale=scale, background=background, method=ms) |

465 | yt = Iq_smeared(q, sphere, pars, delta_qv=delta_qv, max_step=max_step, |

466 | scale=scale, background=background, method=mt) |

467 | |

468 | ## ==== plots |

469 | plt.subplot(211) |

470 | plt.loglog(q_calc, y, '-', label="unsmeared") |

471 | plt.loglog(q, ys, '.', label=ms, hold=True) |

472 | plt.loglog(q, yt, '.', label=mt, hold=True) |

473 | if answer is not None: |

474 | plt.loglog(q, answer, '-', label='igor', hold=True) |

475 | plt.ylabel("y") |

476 | ## Uncomment to see linear fringing in the model; this works particularly |

477 | ## well when q_calc is set to extrapolated and interpolated points. |

478 | #plt.xscale('linear') |

479 | plt.grid(True) |

480 | plt.legend() |

481 | plt.subplot(212) |

482 | err = (ys-yt)/ys |

483 | print mt, "err max:", max(abs(err)), " median:", np.median(abs(err)) |

484 | plt.semilogx(q, err, '-', label=mt) |

485 | ## Uncomment to compare test model to igor |

486 | #yt=answer # use igor output for comparison |

487 | #plt.semilogx(q, (ys-yt)/ys, '-', label='Igor', hold=True) |

488 | plt.xlabel("q") |

489 | plt.ylabel("(%s - %s)/%s"%(ms,mt,ms)) |

490 | plt.legend() |

491 | |

492 | |

493 | if __name__ == "__main__": |

494 | demo() |

495 | import pylab; pylab.show() |