Changeset bc873053 in sasview


Ignore:
Timestamp:
Nov 10, 2015 5:39:44 PM (9 years ago)
Author:
Paul Kienzle <pkienzle@…>
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, release-4.1.1, release-4.1.2, release-4.2.2, release_4.0.1, ticket-1009, ticket-1094-headless, ticket-1242-2d-resolution, ticket-1243, ticket-1249, ticket885, unittest-saveload
Children:
80ba1a2
Parents:
40c69b3
Message:

Fix slit height resolution problem. Still issues with slit width. Refs #472

Files:
4 edited

Legend:

Unmodified
Added
Removed
  • .gitignore

    r98f3234 rbc873053  
    2323**/test/logs 
    2424*.pyc 
     25*.so 
    2526default_categories.json 
    2627 
  • src/sas/models/resolution.py

    rbe0c318 rbc873053  
    1111MINIMUM_RESOLUTION = 1e-8 
    1212 
    13 class Resolution1D(object): 
     13class Resolution(object): 
    1414    """ 
    1515    Abstract base class defining a 1D resolution function. 
     
    3232 
    3333 
    34 class Perfect1D(Resolution1D): 
     34class Perfect1D(Resolution): 
    3535    """ 
    3636    Resolution function to use when there is no actual resolution smearing 
     
    4545 
    4646 
    47 class Pinhole1D(Resolution1D): 
     47class Pinhole1D(Resolution): 
    4848    r""" 
    4949    Pinhole aperture with q-dependent gaussian resolution. 
     
    7777 
    7878 
    79 class Slit1D(Resolution1D): 
     79class Slit1D(Resolution): 
    8080    """ 
    8181    Slit aperture with a complicated resolution function. 
     
    9292    The *weight_matrix* is computed by :func:`slit1d_resolution` 
    9393    """ 
    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): 
    10595        # Remember what width/height was used even though we won't need them 
    10696        # after the weight matrix is constructed 
    10797        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) 
    108109 
    109110        self.q = q.flatten() 
     
    169170            \int_{-\Delta q_v}^{\Delta q_v} I(u) du 
    170171    """ 
    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. 
    175175    q_edges = bin_edges(q_calc) # Note: requires q > 0 
    176176    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 >= q-width) * (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 >= qi-h) * (q_calc <= qi+h) 
     199            weights[i,:] = in_x * np.diff(q_edges) / (2*h) 
     200 
     201    return weights.T 
    196202 
    197203 
     
    212218    function. 
    213219    """ 
    214     height # keep lint happy 
    215220    q_min, q_max = np.min(q), np.max(np.sqrt(q**2 + width**2)) 
    216221    return geometric_extrapolation(q, q_min, q_max) 
     
    233238        ]) 
    234239    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 zero 
    242     qpsq = q**2 - q0**2 
    243     qpsq[qpsq<0] = 0 
    244     return sqrt(qpsq) 
    245240 
    246241 
     
    297292    *points_per_decade* sets the ratio between consecutive steps such 
    298293    that there will be $n$ points used for every factor of 10 increase 
    299     in $q$. 
     294    in *q*. 
    300295 
    301296    If *points_per_decade* is not given, it will be estimated as follows. 
     
    316311    Substituting: 
    317312 
    318     .. math:: 
    319  
    320313        n_\text{extend} = (n-1) (\log q_\text{max} - \log q_n) 
    321             / (\log q_n - \log q_1) 
     314            / (\log q_n - log q_1) 
    322315    """ 
    323316    q = np.sort(q) 
     
    359352 
    360353 
    361 def romberg_slit_1d(q, delta_qv, form, pars): 
     354def romberg_slit_1d(q, width, height, form, pars): 
    362355    """ 
    363356    Romberg integration for slit resolution. 
     
    366359    that make it slow to evaluate but give it good accuracy. 
    367360    """ 
    368     from scipy.integrate import romberg 
     361    from scipy.integrate import romberg, dblquad 
    369362 
    370363    if any(k not in form.info['defaults'] for k in pars.keys()): 
     
    374367                         (", ".join(sorted(extra)), ", ".join(sorted(keys)))) 
    375368 
    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=1e-8) 
    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=1e-8) 
     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=1e-8) 
     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 
    380391    # r should be [float, ...], but it is [array([float]), array([float]),...] 
    381     return np.asarray(r).flatten()/delta_qv[0] 
     392    return result 
    382393 
    383394 
     
    520531 
    521532    def compare(self, q, output, answer, tolerance): 
    522         err = (output - answer)/answer 
    523         idx = abs(err) >= tolerance 
    524         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) 
    526537        np.testing.assert_allclose(output, answer, rtol=tolerance) 
    527538 
     
    609620        data = np.loadtxt(data_string.split('\n')).T 
    610621        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) 
    612623        q_calc = slit_extend_q(interpolate(q, 2*np.pi/radius/20), 
    613                                delta_qv[0], delta_qv[0]) 
     624                               delta_qv[0], 0.) 
    614625        resolution = Slit1D(q, width=delta_qv, height=0, q_calc=q_calc) 
    615626        output = self.Iq_sphere(pars, resolution) 
     
    629640        form = load_model('ellipsoid', dtype='double') 
    630641        q = np.logspace(log10(4e-5),log10(2.5e-2), 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) 
    634645        output = resolution.apply(eval_form(resolution.q_calc, form, pars)) 
    635646        # TODO: 10% is too much error; use better algorithm 
     647        #print np.max(abs(answer-output)/answer) 
    636648        self.compare(q, output, answer, 0.1) 
    637649 
     
    860872 
    861873def _eval_demo_1d(resolution, title): 
     874    import sys 
    862875    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) 
    867894 
    868895    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) 
    870897    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) 
    871905 
    872906    import matplotlib.pyplot as plt 
    873907    plt.loglog(resolution.q_calc, theory, label='unsmeared') 
    874908    plt.loglog(resolution.q, Iq, label='smeared', hold=True) 
     909    plt.loglog(resolution.q, Iq_romb, label='romberg smeared', hold=True) 
    875910    plt.legend() 
    876911    plt.title(title) 
     
    879914 
    880915def demo_pinhole_1d(): 
    881     q = np.logspace(-3, -1, 400) 
     916    q = np.logspace(-4, np.log10(0.2), 400) 
    882917    q_width = 0.1*q 
    883918    resolution = Pinhole1D(q, q_width) 
     
    885920 
    886921def 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)) 
    892929 
    893930def demo(): 
     
    895932    plt.subplot(121) 
    896933    demo_pinhole_1d() 
     934    #plt.yscale('linear') 
    897935    plt.subplot(122) 
    898936    demo_slit_1d() 
     937    #plt.yscale('linear') 
    899938    plt.show() 
    900939 
  • test/sasfit/test/utest_fit_line.py

    racf8e4a5 rbc873053  
    1111from sas.models.LineModel import LineModel 
    1212from sas.models.Constant import Constant 
     13 
     14from bumps import fitters 
     15try: 
     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=[]) 
     20except: 
     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 
    1328 
    1429class testFitModule(unittest.TestCase): 
     
    6782 
    6883    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) 
    7685        self.fit_single(isdream=(alg=='dream')) 
    7786 
  • test/sasfit/test/utest_fit_smeared.py

    racf8e4a5 rbc873053  
    1313from sas.models.CylinderModel import CylinderModel 
    1414from sas.models.SphereModel import SphereModel 
     15 
     16 
     17from bumps import fitters 
     18try: 
     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=[]) 
     23except: 
     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 
    1531 
    1632class testFitModule(unittest.TestCase): 
     
    5773            Cylinder fit with dispersion 
    5874        """ 
    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', {}) 
    6476        self._dispersion(fitter = Fit()) 
    6577 
Note: See TracChangeset for help on using the changeset viewer.