Changeset bc873053 in sasview
 Timestamp:
 Nov 10, 2015 5:39:44 PM (8 years ago)
 Branches:
 master, ESS_GUI, ESS_GUI_Docs, ESS_GUI_batch_fitting, ESS_GUI_bumps_abstraction, ESS_GUI_iss1116, ESS_GUI_iss879, ESS_GUI_iss959, ESS_GUI_opencl, ESS_GUI_ordering, ESS_GUI_sync_sascalc, costrafo411, magnetic_scatt, release4.1.1, release4.1.2, release4.2.2, release_4.0.1, ticket1009, ticket1094headless, ticket12422dresolution, ticket1243, ticket1249, ticket885, unittestsaveload
 Children:
 80ba1a2
 Parents:
 40c69b3
 Files:

 4 edited
Legend:
 Unmodified
 Added
 Removed

.gitignore
r98f3234 rbc873053 23 23 **/test/logs 24 24 *.pyc 25 *.so 25 26 default_categories.json 26 27 
src/sas/models/resolution.py
rbe0c318 rbc873053 11 11 MINIMUM_RESOLUTION = 1e8 12 12 13 class Resolution 1D(object):13 class Resolution(object): 14 14 """ 15 15 Abstract base class defining a 1D resolution function. … … 32 32 33 33 34 class Perfect1D(Resolution 1D):34 class Perfect1D(Resolution): 35 35 """ 36 36 Resolution function to use when there is no actual resolution smearing … … 45 45 46 46 47 class Pinhole1D(Resolution 1D):47 class Pinhole1D(Resolution): 48 48 r""" 49 49 Pinhole aperture with qdependent gaussian resolution. … … 77 77 78 78 79 class Slit1D(Resolution 1D):79 class Slit1D(Resolution): 80 80 """ 81 81 Slit aperture with a complicated resolution function. … … 92 92 The *weight_matrix* is computed by :func:`slit1d_resolution` 93 93 """ 94 def __init__(self, q, width, height, q_calc=None): 95 # TODO: maybe issue warnings rather than raising errors 96 if not np.isscalar(width): 97 if np.any(np.diff(width) > 0.0): 98 raise ValueError("Slit resolution requires fixed width slits") 99 width = width[0] 100 if not np.isscalar(height): 101 if np.any(np.diff(height) > 0.0): 102 raise ValueError("Slit resolution requires fixed height slits") 103 height = height[0] 104 94 def __init__(self, q, width, height=0., q_calc=None): 105 95 # Remember what width/height was used even though we won't need them 106 96 # after the weight matrix is constructed 107 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) 108 109 109 110 self.q = q.flatten() … … 169 170 \int_{\Delta q_v}^{\Delta q_v} I(u) du 170 171 """ 171 if width == 0.0 and height == 0.0: 172 #print "condition zero" 173 return 1 174 172 #np.set_printoptions(precision=6, linewidth=10000) 173 174 # The current algorithm is a midpoint rectangle rule. 175 175 q_edges = bin_edges(q_calc) # Note: requires q > 0 176 176 q_edges[q_edges<0.0] = 0.0 # clip edges below zero 177 178 #np.set_printoptions(linewidth=10000) 179 if width <= 100.0 * height or height == 0: 180 # The current algorithm is a midpoint rectangle rule. In the test case, 181 # neither trapezoid nor Simpson's rule improved the accuracy. 182 #print "condition h", q_edges.shape, q.shape, q_calc.shape 183 weights = np.zeros((len(q), len(q_calc)), 'd') 184 for i, qi in enumerate(q): 185 weights[i, :] = np.diff(q_to_u(q_edges, qi)) 186 weights /= width 187 weights = weights.T 188 else: 189 #print "condition w" 190 # Make q_calc into a row vector, and q into a column vector 191 q, q_calc = q[None,:], q_calc[:,None] 192 in_x = (q_calc >= qwidth) * (q_calc <= q+width) 193 weights = np.diff(q_edges)[:,None] * in_x 194 195 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 196 202 197 203 … … 212 218 function. 213 219 """ 214 height # keep lint happy215 220 q_min, q_max = np.min(q), np.max(np.sqrt(q**2 + width**2)) 216 221 return geometric_extrapolation(q, q_min, q_max) … … 233 238 ]) 234 239 return edges 235 236 237 def q_to_u(q, q0):238 """239 Convert *q* values to *u* values for the integral computed at *q0*.240 """241 # array([value])**2  value**2 is not always zero242 qpsq = q**2  q0**2243 qpsq[qpsq<0] = 0244 return sqrt(qpsq)245 240 246 241 … … 297 292 *points_per_decade* sets the ratio between consecutive steps such 298 293 that there will be $n$ points used for every factor of 10 increase 299 in $q$.294 in *q*. 300 295 301 296 If *points_per_decade* is not given, it will be estimated as follows. … … 316 311 Substituting: 317 312 318 .. math::319 320 313 n_\text{extend} = (n1) (\log q_\text{max}  \log q_n) 321 / (\log q_n  \log q_1)314 / (\log q_n  log q_1) 322 315 """ 323 316 q = np.sort(q) … … 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 
test/sasfit/test/utest_fit_line.py
racf8e4a5 rbc873053 11 11 from sas.models.LineModel import LineModel 12 12 from sas.models.Constant import Constant 13 14 from bumps import fitters 15 try: 16 from bumps.options import FIT_CONFIG 17 def set_fitter(alg, opts): 18 FIT_CONFIG.selected_id = alg 19 FIT_CONFIG.values[alg].update(opts, monitors=[]) 20 except: 21 # CRUFT: Bumps changed its handling of fit options around 0.7.5.6 22 def set_fitter(alg, opts): 23 #print "fitting",alg,opts 24 #print "options",fitters.FIT_OPTIONS[alg].__dict__ 25 fitters.FIT_DEFAULT = alg 26 fitters.FIT_OPTIONS[alg].options.update(opts, monitors=[]) 27 13 28 14 29 class testFitModule(unittest.TestCase): … … 67 82 68 83 def fit_bumps(self, alg, **opts): 69 #Importing the Fit module 70 from bumps import fitters 71 fitters.FIT_DEFAULT = alg 72 fitters.FIT_OPTIONS[alg].options.update(opts) 73 fitters.FIT_OPTIONS[alg].options.update(monitors=[]) 74 #print "fitting",alg,opts 75 #kprint "options",fitters.FIT_OPTIONS[alg].__dict__ 84 set_fitter(alg, opts) 76 85 self.fit_single(isdream=(alg=='dream')) 77 86 
test/sasfit/test/utest_fit_smeared.py
racf8e4a5 rbc873053 13 13 from sas.models.CylinderModel import CylinderModel 14 14 from sas.models.SphereModel import SphereModel 15 16 17 from bumps import fitters 18 try: 19 from bumps.options import FIT_CONFIG 20 def set_fitter(alg, opts): 21 FIT_CONFIG.selected_id= alg 22 FIT_CONFIG.values[alg].update(opts, monitors=[]) 23 except: 24 # CRUFT: Bumps changed its handling of fit options around 0.7.5.6 25 def set_fitter(alg, opts): 26 #print "fitting",alg,opts 27 #print "options",fitters.FIT_OPTIONS[alg].__dict__ 28 fitters.FIT_DEFAULT = alg 29 fitters.FIT_OPTIONS[alg].options.update(opts, monitors=[]) 30 15 31 16 32 class testFitModule(unittest.TestCase): … … 57 73 Cylinder fit with dispersion 58 74 """ 59 alg = 'lm' 60 from bumps import fitters 61 fitters.FIT_DEFAULT = alg 62 #fitters.FIT_OPTIONS[alg].options.update(opts) 63 fitters.FIT_OPTIONS[alg].options.update(monitors=[]) 75 set_fitter('lm', {}) 64 76 self._dispersion(fitter = Fit()) 65 77
Note: See TracChangeset
for help on using the changeset viewer.