Changeset fdc538a in sasmodels
 Timestamp:
 Jan 27, 2016 3:47:17 PM (8 years ago)
 Branches:
 master, core_shell_microgels, costrafo411, magnetic_model, release_v0.94, release_v0.95, ticket1257vesicleproduct, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests
 Children:
 09ebe8c
 Parents:
 d15a908
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

sasmodels/resolution.py
r5258859 rfdc538a 135 135 """ 136 136 #print("apply shapes", theory.shape, weight_matrix.shape) 137 Iq = np.dot(theory[None, :], weight_matrix)137 Iq = np.dot(theory[None, :], weight_matrix) 138 138 #print("result shape",Iq.shape) 139 139 return Iq.flatten() … … 153 153 # neither trapezoid nor Simpson's rule improved the accuracy. 154 154 edges = bin_edges(q_calc) 155 edges[edges <0.0] = 0.0 # clip edges below zero156 G = erf( (edges[:,None]  q[None,:]) / (sqrt(2.0)*q_width)[None,:])155 edges[edges < 0.0] = 0.0 # clip edges below zero 156 G = erf((edges[:, None]  q[None, :]) / (sqrt(2.0)*q_width)[None, :]) 157 157 weights = G[1:]  G[:1] 158 weights /= np.sum(weights, axis=0)[None, :]158 weights /= np.sum(weights, axis=0)[None, :] 159 159 return weights 160 160 … … 287 287 # The current algorithm is a midpoint rectangle rule. 288 288 q_edges = bin_edges(q_calc) # Note: requires q > 0 289 q_edges[q_edges <0.0] = 0.0 # clip edges below zero289 q_edges[q_edges < 0.0] = 0.0 # clip edges below zero 290 290 weights = np.zeros((len(q), len(q_calc)), 'd') 291 291 … … 306 306 #print(qi  h, qi + h) 307 307 #print(in_x + abs_x) 308 weights[i, :] = (in_x + abs_x) * np.diff(q_edges) / (2*h)308 weights[i, :] = (in_x + abs_x) * np.diff(q_edges) / (2*h) 309 309 else: 310 310 L = n_height 311 311 for k in range(L, L+1): 312 weights[i, :] += _q_perp_weights(q_edges, qi+k*h/L, w)313 weights[i, :] /= 2*L + 1312 weights[i, :] += _q_perp_weights(q_edges, qi+k*h/L, w) 313 weights[i, :] /= 2*L + 1 314 314 315 315 return weights.T … … 358 358 logscaled data. 359 359 """ 360 if len(x) < 2 or (np.diff(x) <0).any():360 if len(x) < 2 or (np.diff(x) < 0).any(): 361 361 raise ValueError("Expected bins to be an increasing set") 362 362 edges = np.hstack([ … … 373 373 """ 374 374 step = np.diff(q) 375 index = step >max_step375 index = step > max_step 376 376 if np.any(index): 377 377 inserts = [] 378 for q_i, step_i in zip(q[:1][index],step[index]):378 for q_i, step_i in zip(q[:1][index], step[index]): 379 379 n = np.ceil(step_i/max_step) 380 inserts.extend(q_i + np.arange(1, n)*(step_i/n))380 inserts.extend(q_i + np.arange(1, n)*(step_i/n)) 381 381 # Extend a couple of fringes beyond the end of the data 382 inserts.extend(q[1] + np.arange(1, 8)*max_step)383 q_calc = np.sort(np.hstack((q, inserts)))382 inserts.extend(q[1] + np.arange(1, 8)*max_step) 383 q_calc = np.sort(np.hstack((q, inserts))) 384 384 else: 385 385 q_calc = q … … 399 399 if q_min < q[0]: 400 400 if q_min <= 0: q_min = q_min*MIN_Q_SCALE_FOR_NEGATIVE_Q_EXTRAPOLATION 401 n_low = np.ceil((q[0]q_min) / (q[1]q[0])) if q[1] >q[0] else 15401 n_low = np.ceil((q[0]q_min) / (q[1]q[0])) if q[1] > q[0] else 15 402 402 q_low = np.linspace(q_min, q[0], n_low+1)[:1] 403 403 else: 404 404 q_low = [] 405 405 if q_max > q[1]: 406 n_high = np.ceil((q_maxq[1]) / (q[1]q[2])) if q[1] >q[2] else 15406 n_high = np.ceil((q_maxq[1]) / (q[1]q[2])) if q[1] > q[2] else 15 407 407 q_high = np.linspace(q[1], q_max, n_high+1)[1:] 408 408 else: … … 452 452 if q_min < 0: q_min = q[0]*MIN_Q_SCALE_FOR_NEGATIVE_Q_EXTRAPOLATION 453 453 n_low = log_delta_q * (log(q[0])log(q_min)) 454 q_low 454 q_low = np.logspace(log10(q_min), log10(q[0]), np.ceil(n_low)+1)[:1] 455 455 else: 456 456 q_low = [] … … 489 489 that make it slow to evaluate but give it good accuracy. 490 490 """ 491 from scipy.integrate import romberg , dblquad491 from scipy.integrate import romberg 492 492 493 493 if any(k not in form.info['defaults'] for k in pars.keys()): … … 520 520 result[i] = r/(2*h) 521 521 else: 522 w_grid = np.linspace(0, w, 21)[None, :]523 h_grid = np.linspace(h, h, 23)[:, None]522 w_grid = np.linspace(0, w, 21)[None, :] 523 h_grid = np.linspace(h, h, 23)[:, None] 524 524 u = sqrt((qi+h_grid)**2 + w_grid**2) 525 525 Iu = np.interp(u, q_calc, Iq) 526 526 #print(np.trapz(Iu, w_grid, axis=1)) 527 Is = np.trapz(np.trapz(Iu, w_grid, axis=1), h_grid[:, 0])527 Is = np.trapz(np.trapz(Iu, w_grid, axis=1), h_grid[:, 0]) 528 528 result[i] = Is / (2*h*w) 529 """ 530 r, err = dblquad(_int_wh, h, h, lambda h: 0., lambda h: w, 531 args=(qi,)) 532 result[i] = r/(w*2*h) 533 """ 529 # from scipy.integrate import dblquad 530 # r, err = dblquad(_int_wh, h, h, lambda h: 0., lambda h: w, 531 # args=(qi,)) 532 # result[i] = r/(w*2*h) 534 533 535 534 # r should be [float, ...], but it is [array([float]), array([float]),...] … … 553 552 554 553 _fn = lambda q, q0, dq: eval_form(q, form, pars)*gaussian(q, q0, dq) 555 r = [romberg(_fn, max(qinsigma*dqi, 1e10*q[0]), qi+nsigma*dqi, args=(qi, dqi),556 divmax=100, vec_func=True, tol=0, rtol=1e8)557 for qi, dqi in zip(q,q_width)]554 r = [romberg(_fn, max(qinsigma*dqi, 1e10*q[0]), qi+nsigma*dqi, 555 args=(qi, dqi), divmax=100, vec_func=True, tol=0, rtol=1e8) 556 for qi, dqi in zip(q, q_width)] 558 557 return np.asarray(r).flatten() 559 558 … … 595 594 theory = self.Iq(resolution.q_calc) 596 595 output = resolution.apply(theory) 597 answer = [ 9.0618, 8.6402, 8.1187, 7.1392, 6.1528, 598 5.5555, 4.5584, 3.5606, 2.5623, 2.0000 ] 596 answer = [ 597 9.0618, 8.6402, 8.1187, 7.1392, 6.1528, 598 5.5555, 4.5584, 3.5606, 2.5623, 2.0000, 599 ] 599 600 np.testing.assert_allclose(output, answer, atol=1e4) 600 601 … … 608 609 theory = 1000*self.Iq(resolution.q_calc**4) 609 610 output = resolution.apply(theory) 610 answer = [ 8.85785, 8.43012, 7.92687, 6.94566, 6.03660, 611 5.40363, 4.40655, 3.40880, 2.41058, 2.00000 ] 611 answer = [ 612 8.85785, 8.43012, 7.92687, 6.94566, 6.03660, 613 5.40363, 4.40655, 3.40880, 2.41058, 2.00000, 614 ] 612 615 np.testing.assert_allclose(output, answer, atol=1e4) 613 616 … … 620 623 theory = self.Iq(resolution.q_calc) 621 624 output = resolution.apply(theory) 622 answer = [ 11.0, 10.0, 9.0, 8.0, 7.0, 6.0, 5.0, 4.0, 3.0, 2.0 ] 625 answer = [ 626 11.0, 10.0, 9.0, 8.0, 7.0, 6.0, 5.0, 4.0, 3.0, 2.0, 627 ] 623 628 np.testing.assert_allclose(output, answer, atol=1e4) 624 629 … … 632 637 theory = self.Iq(resolution.q_calc) 633 638 output = resolution.apply(theory) 634 answer = [ 11.0, 10.0, 9.0, 8.0, 7.0, 6.0, 5.0, 4.0, 3.0, 2.0 ] 639 answer = [ 640 11.0, 10.0, 9.0, 8.0, 7.0, 6.0, 5.0, 4.0, 3.0, 2.0, 641 ] 635 642 np.testing.assert_allclose(output, answer, atol=1e4) 636 643 … … 652 659 theory = 12.01000.0*resolution.q_calc 653 660 output = resolution.apply(theory) 654 answer = [ 10.44785079, 9.84991299, 8.98101708, 655 7.99906585, 6.99998311, 6.00001689, 656 5.00093415, 4.01898292, 3.15008701, 2.55214921] 661 answer = [ 662 10.44785079, 9.84991299, 8.98101708, 663 7.99906585, 6.99998311, 6.00001689, 664 5.00093415, 4.01898292, 3.15008701, 2.55214921, 665 ] 657 666 np.testing.assert_allclose(output, answer, atol=1e8) 658 667 … … 783 792 } 784 793 form = load_model('ellipsoid', dtype='double') 785 q = np.logspace(log10(4e5), log10(2.5e2), 68)794 q = np.logspace(log10(4e5), log10(2.5e2), 68) 786 795 width, height = 0.117, 0. 787 796 resolution = Slit1D(q, width=width, height=height) … … 1071 1080 #h = 0.0277790 1072 1081 resolution = Slit1D(q, w, h) 1073 _eval_demo_1d(resolution, title="(%g,%g) Slit Resolution"%(w, h))1082 _eval_demo_1d(resolution, title="(%g,%g) Slit Resolution"%(w, h)) 1074 1083 1075 1084 def demo():
Note: See TracChangeset
for help on using the changeset viewer.