Changeset dbde9f8 in sasmodels
 Timestamp:
 Nov 10, 2015 5:01:46 PM (9 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:
 3a45c2c, eb588ef
 Parents:
 2781b2e
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

sasmodels/resolution.py
r346bc88 rdbde9f8 74 74 75 75 def apply(self, theory): 76 print "q calc", self.q_calc77 print "weights", self.weight_matrix.shape78 76 return apply_resolution_matrix(self.weight_matrix, theory) 79 77 … … 94 92 The *weight_matrix* is computed by :func:`slit1d_resolution` 95 93 """ 96 def __init__(self, q, width, height, q_calc=None): 97 # TODO: maybe issue warnings rather than raising errors 98 if not np.isscalar(width): 99 if np.any(np.diff(width) > 0.0): 100 raise ValueError("Slit resolution requires fixed width slits") 101 width = width[0] 102 if not np.isscalar(height): 103 if np.any(np.diff(height) > 0.0): 104 raise ValueError("Slit resolution requires fixed height slits") 105 height = height[0] 106 94 def __init__(self, q, width, height=0., q_calc=None): 107 95 # Remember what width/height was used even though we won't need them 108 96 # after the weight matrix is constructed 109 97 self.width, self.height = width, height 98 99 # Allow independent resolution on each point even though it is not 100 # needed in practice. 101 if np.isscalar(width): 102 width = np.ones(len(q))*width 103 else: 104 width = np.asarray(width) 105 if np.isscalar(height): 106 height = np.ones(len(q))*height 107 else: 108 height = np.asarray(height) 110 109 111 110 self.q = q.flatten() … … 171 170 \int_{\Delta q_v}^{\Delta q_v} I(u) du 172 171 """ 173 if width == 0.0 and height == 0.0: 174 #print "condition zero" 175 return 1 176 172 #np.set_printoptions(precision=6, linewidth=10000) 173 174 # The current algorithm is a midpoint rectangle rule. 177 175 q_edges = bin_edges(q_calc) # Note: requires q > 0 178 176 q_edges[q_edges<0.0] = 0.0 # clip edges below zero 179 180 #np.set_printoptions(linewidth=10000) 181 if width <= 100.0 * height or height == 0: 182 # The current algorithm is a midpoint rectangle rule. In the test case, 183 # neither trapezoid nor Simpson's rule improved the accuracy. 184 #print "condition h", q_edges.shape, q.shape, q_calc.shape 185 weights = np.zeros((len(q), len(q_calc)), 'd') 186 for i, qi in enumerate(q): 187 weights[i, :] = np.diff(q_to_u(q_edges, qi)) 188 weights /= width 189 weights = weights.T 190 else: 191 #print "condition w" 192 # Make q_calc into a row vector, and q into a column vector 193 q, q_calc = q[None,:], q_calc[:,None] 194 in_x = (q_calc >= qwidth) * (q_calc <= q+width) 195 weights = np.diff(q_edges)[:,None] * in_x 196 197 return weights 177 weights = np.zeros((len(q), len(q_calc)), 'd') 178 179 for i, (qi, w, h) in enumerate(zip(q, width, height)): 180 if w == 0. and h == 0.: 181 # Perfect resolution, so return the theory value directly. 182 # Note: assumes that q is a subset of q_calc. If qi need not be 183 # in q_calc, then we can do a weighted interpolation by looking 184 # up qi in q_calc, then weighting the result by the relative 185 # distance to the neighbouring points. 186 weights[i, :] = (q_calc == qi) 187 elif h == 0 or w > 100.0 * h: 188 # Convert bin edges from q to u 189 u_limit = np.sqrt(qi**2 + w**2) 190 u_edges = q_edges**2  qi**2 191 u_edges[q_edges < qi] = 0. 192 u_edges[q_edges > u_limit] = u_limit**2  qi**2 193 weights[i, :] = np.diff(np.sqrt(u_edges))/w 194 #print "i, qi",i,qi,qi+width 195 #print q_calc 196 #print weights[i,:] 197 elif w == 0: 198 in_x = (q_calc >= qih) * (q_calc <= qi+h) 199 weights[i,:] = in_x * np.diff(q_edges) / (2*h) 200 201 return weights.T 198 202 199 203 … … 214 218 function. 215 219 """ 216 height # keep lint happy217 220 q_min, q_max = np.min(q), np.max(np.sqrt(q**2 + width**2)) 218 221 return geometric_extrapolation(q, q_min, q_max) … … 235 238 ]) 236 239 return edges 237 238 239 def q_to_u(q, q0):240 """241 Convert *q* values to *u* values for the integral computed at *q0*.242 """243 # array([value])**2  value**2 is not always zero244 qpsq = q**2  q0**2245 qpsq[qpsq<0] = 0246 return sqrt(qpsq)247 240 248 241 … … 359 352 360 353 361 def romberg_slit_1d(q, delta_qv, form, pars):354 def romberg_slit_1d(q, width, height, form, pars): 362 355 """ 363 356 Romberg integration for slit resolution. … … 366 359 that make it slow to evaluate but give it good accuracy. 367 360 """ 368 from scipy.integrate import romberg 361 from scipy.integrate import romberg, dblquad 369 362 370 363 if any(k not in form.info['defaults'] for k in pars.keys()): … … 374 367 (", ".join(sorted(extra)), ", ".join(sorted(keys)))) 375 368 376 _fn = lambda u, q0: eval_form(sqrt(q0**2 + u**2), form, pars) 377 r = [romberg(_fn, 0, delta_qv[0], args=(qi,), 378 divmax=100, vec_func=True, tol=0, rtol=1e8) 379 for qi in q] 369 if np.isscalar(width): 370 width = [width]*len(q) 371 if np.isscalar(height): 372 height = [height]*len(q) 373 _int_w = lambda w, qi: eval_form(sqrt(qi**2 + w**2), form, pars) 374 _int_h = lambda h, qi: eval_form(abs(qi+h), form, pars) 375 _int_wh = lambda w, h, qi: eval_form(sqrt((qi+h)**2 + w**2), form, pars) 376 result = np.empty(len(q)) 377 for i, (qi, w, h) in enumerate(zip(q, width, height)): 378 if h == 0.: 379 r = romberg(_int_w, 0, w, args=(qi,), 380 divmax=100, vec_func=True, tol=0, rtol=1e8) 381 result[i] = r/w 382 elif w == 0.: 383 r = romberg(_int_h, h, h, args=(qi,), 384 divmax=100, vec_func=True, tol=0, rtol=1e8) 385 result[i] = r/(2*h) 386 else: 387 r, err = dblquad(_int_wh, h, h, lambda h: 0., lambda h: width, 388 args=(qi,)) 389 result[i] = r/(w*2*h) 390 380 391 # r should be [float, ...], but it is [array([float]), array([float]),...] 381 return np.asarray(r).flatten()/delta_qv[0]392 return result 382 393 383 394 … … 520 531 521 532 def compare(self, q, output, answer, tolerance): 522 err = (output  answer)/answer523 idx = abs(err) >= tolerance524 problem = zip(q[idx], output[idx], answer[idx], err[idx])525 print "\n".join(str(v) for v in problem)533 #err = (output  answer)/answer 534 #idx = abs(err) >= tolerance 535 #problem = zip(q[idx], output[idx], answer[idx], err[idx]) 536 #print "\n".join(str(v) for v in problem) 526 537 np.testing.assert_allclose(output, answer, rtol=tolerance) 527 538 … … 609 620 data = np.loadtxt(data_string.split('\n')).T 610 621 q, delta_qv, _, answer = data 611 answer = romberg_slit_1d(q, delta_qv, self.model, pars)622 answer = romberg_slit_1d(q, delta_qv, 0., self.model, pars) 612 623 q_calc = slit_extend_q(interpolate(q, 2*np.pi/radius/20), 613 delta_qv[0], delta_qv[0])624 delta_qv[0], 0.) 614 625 resolution = Slit1D(q, width=delta_qv, height=0, q_calc=q_calc) 615 626 output = self.Iq_sphere(pars, resolution) … … 629 640 form = load_model('ellipsoid', dtype='double') 630 641 q = np.logspace(log10(4e5),log10(2.5e2), 68) 631 delta_qv = [0.117]632 resolution = Slit1D(q, width= delta_qv, height=0)633 answer = romberg_slit_1d(q, delta_qv, form, pars)642 width, height = 0.117, 0. 643 resolution = Slit1D(q, width=width, height=height) 644 answer = romberg_slit_1d(q, width, height, form, pars) 634 645 output = resolution.apply(eval_form(resolution.q_calc, form, pars)) 635 646 # TODO: 10% is too much error; use better algorithm 647 #print np.max(abs(answeroutput)/answer) 636 648 self.compare(q, output, answer, 0.1) 637 649 … … 860 872 861 873 def _eval_demo_1d(resolution, title): 874 import sys 862 875 from sasmodels import core 863 from sasmodels.models import cylinder 864 ## or alternatively: 865 # cylinder = core.load_model_definition('cylinder') 866 model = core.load_model(cylinder) 876 name = sys.argv[1] if len(sys.argv) > 1 else 'cylinder' 877 878 if name == 'cylinder': 879 pars = {'length':210, 'radius':500} 880 elif name == 'teubner_strey': 881 pars = {'a2':0.003, 'c1':1e4, 'c2':1e10, 'background':0.312643} 882 elif name == 'sphere' or name == 'spherepy': 883 pars = TEST_PARS_SLIT_SPHERE 884 elif name == 'ellipsoid': 885 pars = { 886 'scale':0.05, 887 'rpolar':500, 'requatorial':15000, 888 'sld':6, 'solvent_sld': 1, 889 } 890 else: 891 pars = {} 892 defn = core.load_model_definition(name) 893 model = core.load_model(defn) 867 894 868 895 kernel = core.make_kernel(model, [resolution.q_calc]) 869 theory = core.call_kernel(kernel, {'length':210, 'radius':500})896 theory = core.call_kernel(kernel, pars) 870 897 Iq = resolution.apply(theory) 898 899 if isinstance(resolution, Slit1D): 900 width, height = resolution.width, resolution.height 901 Iq_romb = romberg_slit_1d(resolution.q, width, height, model, pars) 902 else: 903 dq = resolution.q_width 904 Iq_romb = romberg_pinhole_1d(resolution.q, dq, model, pars) 871 905 872 906 import matplotlib.pyplot as plt 873 907 plt.loglog(resolution.q_calc, theory, label='unsmeared') 874 908 plt.loglog(resolution.q, Iq, label='smeared', hold=True) 909 plt.loglog(resolution.q, Iq_romb, label='romberg smeared', hold=True) 875 910 plt.legend() 876 911 plt.title(title) … … 879 914 880 915 def demo_pinhole_1d(): 881 q = np.logspace( 3, 1, 400)916 q = np.logspace(4, np.log10(0.2), 400) 882 917 q_width = 0.1*q 883 918 resolution = Pinhole1D(q, q_width) … … 885 920 886 921 def demo_slit_1d(): 887 q = np.logspace(3, 1, 400) 888 qx_width = 0.005 889 qy_width = 0.0 890 resolution = Slit1D(q, qx_width, qy_width) 891 _eval_demo_1d(resolution, title="0.005 Qx Slit Resolution") 922 q = np.logspace(4, np.log10(0.2), 400) 923 w = h = 0. 924 #w = 0.005 925 w = 0.0277790 926 #h = 0.00277790 927 resolution = Slit1D(q, w, h) 928 _eval_demo_1d(resolution, title="(%g,%g) Slit Resolution"%(w,h)) 892 929 893 930 def demo(): … … 895 932 plt.subplot(121) 896 933 demo_pinhole_1d() 934 #plt.yscale('linear') 897 935 plt.subplot(122) 898 936 demo_slit_1d() 937 #plt.yscale('linear') 899 938 plt.show() 900 939
Note: See TracChangeset
for help on using the changeset viewer.