Changeset a3f125f0 in sasview for src


Ignore:
Timestamp:
Mar 30, 2015 10:51:41 AM (10 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:
edfc8ac
Parents:
d5419f7f
Message:

refactor resolution calculation to enable sasmodels resolution calcuator for pinhole smearing, but don't enable it yet

Location:
src/sas
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • src/sas/models/qsmearing.py

    r9f7fbd9 ra3f125f0  
    5959    if _found_resolution == True: 
    6060        return QSmearer(data1D, model) 
     61        #return pinhole_smear(data1D, model) 
    6162 
    6263    # Look for slit smearing data 
    … …  
    8182    # If we found slit smearing data, return a slit smearer 
    8283    if _found_slit == True: 
    83         return PySlitSmearer(data1D, model) 
     84        #return SlitSmearer(data1D, model) 
     85        return slit_smear(data1D, model) 
    8486    return None 
    8587             
    … …  
    198200                iq_in[len(iq_in) - 1] = iq_in_high[0] 
    199201            # Append the extrapolated points to the data points 
    200             if self.nbins_low > 0:                              
     202            if self.nbins_low > 0: 
    201203                iq_in_temp = numpy.append(iq_in_low, iq_in) 
    202204            if self.nbins_high > 0: 
    … …  
    586588 
    587589 
    588 from .resolution import Slit1D 
    589 class PySlitSmearer(object): 
    590     def __init__(self, data1D, model = None): 
     590from .resolution import Slit1D, Pinhole1D 
     591class PySmear(object): 
     592    """ 
     593    Wrapper for pure python sasmodels resolution functions. 
     594    """ 
     595    def __init__(self, resolution, model): 
    591596        self.model = model 
    592  
    593         q = data1D.x 
    594         width = data1D.dxw if data1D.dxw is not None else 0 
    595         height = data1D.dxl if data1D.dxl is not None else 0 
    596         # TODO: width and height seem to be reversed 
    597         self.resolution = Slit1D(q, height, width) 
    598  
    599     def __call__(self, iq_in, first_bin=0, last_bin=None): 
    600         if last_bin is None or last_bin >= len(iq_in): 
    601             last_bin = len(iq_in) - 1 
     597        self.resolution = resolution 
     598        self.offset = numpy.searchsorted(self.resolution.q_calc, self.resolution.q[0]) 
     599 
     600    def apply(self, iq_in, first_bin=0, last_bin=None): 
     601        """ 
     602        Apply the resolution function to the data. 
     603 
     604        Note that this is called with iq_in matching data.x, but with 
     605        iq_in[first_bin:last_bin] set to theory values for these bins, 
     606        and the remainder left undefined.  The first_bin, last_bin values 
     607        should be those returned from get_bin_range. 
     608 
     609        The returned value is of the same length as iq_in, with the range 
     610        first_bin:last_bin set to the resolution smeared values. 
     611        """ 
     612        if last_bin is None: last_bin = len(iq_in) 
     613        start, end = first_bin + self.offset, last_bin + self.offset 
    602614        q_calc = self.resolution.q_calc 
    603615        iq_calc = numpy.empty_like(q_calc) 
    604         iq_calc[:first_bin] = 0 
    605         iq_calc[first_bin:last_bin+1] = iq_in 
    606         if last_bin < len(q_calc)-1: 
    607             iq_calc[last_bin:] = self.model.evalDistribution(q_calc[last_bin:]) 
    608         iq_out = self.resolution.apply(iq_calc) 
    609         return iq_out[first_bin:last_bin+1] 
     616        if start > 0: 
     617            iq_calc[:start] = self.model.evalDistribution(q_calc[:start]) 
     618        if end+1 < len(q_calc): 
     619            iq_calc[end+1:] = self.model.evalDistribution(q_calc[end+1:]) 
     620        iq_calc[start:end+1] = iq_in[first_bin:last_bin+1] 
     621        smeared = self.resolution.apply(iq_calc) 
     622        return smeared 
     623    __call__ = apply 
    610624 
    611625    def get_bin_range(self, q_min=None, q_max=None): 
    612626        """ 
    613  
    614         :param q_min: minimum q-value to smear 
    615         :param q_max: maximum q-value to smear 
    616  
    617         """ 
    618         # assume the data values are sorted 
    619         first = numpy.searchsorted(self.resolution.q, q_min) 
    620         # assume that the resolution is much larger than the q range 
    621         last = len(self.resolution.q)-1 
    622         return first, last 
     627        For a given q_min, q_max, find the corresponding indices in the data. 
     628 
     629        Returns first, last. 
     630 
     631        Note that these are indexes into q from the data, not the q_calc 
     632        needed by the resolution function.  Note also that these are the 
     633        indices, not the range limits.  That is, the complete range will be 
     634        q[first:last+1]. 
     635        """ 
     636        q = self.resolution.q 
     637        first = numpy.searchsorted(q, q_min) 
     638        last = numpy.searchsorted(q, q_max) 
     639        return first, min(last,len(q)-1) 
     640 
     641def slit_smear(data, model=None): 
     642    q = data.x 
     643    width = data.dxw if data.dxw is not None else 0 
     644    height = data.dxl if data.dxl is not None else 0 
     645    # TODO: width and height seem to be reversed 
     646    return PySmear(Slit1D(q, height, width), model) 
     647 
     648def pinhole_smear(data, model=None): 
     649    q = data.x 
     650    width = data.dx if data.dx is not None else 0 
     651    return PySmear(Pinhole1D(q, width), model) 

  • src/sas/models/resolution.py

    r9f7fbd9 ra3f125f0  
    5656    be estimated from the *q* and *q_width*. 
    5757    """ 
    58     def __init__(self, q, q_width, q_calc=None): 
     58    def __init__(self, q, q_width, q_calc=None, nsigma=3): 
    5959        #*min_step* is the minimum point spacing to use when computing the 
    6060        #underlying model.  It should be on the order of 
    … …  
    6868        # default to Perfect1D if the pinhole geometry is not defined. 
    6969        self.q, self.q_width = q, q_width 
    70         self.q_calc = pinhole_extend_q(q, q_width) \ 
     70        self.q_calc = pinhole_extend_q(q, q_width, nsigma=nsigma) \ 
    7171            if q_calc is None else np.sort(q_calc) 
    7272        self.weight_matrix = pinhole_resolution(self.q_calc, 
    … …  
    203203    """ 
    204204    q_min, q_max = np.min(q - nsigma*q_width), np.max(q + nsigma*q_width) 
    205     return geometric_extrapolation(q, q_min, q_max) 
     205    return linear_extrapolation(q, q_min, q_max) 
    206206 
    207207 
    … …  
    267267    """ 
    268268    Extrapolate *q* out to [*q_min*, *q_max*] using the step size in *q* as 
    269     a guide.  Extrapolation below uses steps the same size as the first 
    270     interval.  Extrapolation above uses steps the same size as the final 
     269    a guide.  Extrapolation below uses about the same size as the first 
     270    interval.  Extrapolation above uses about the same size as the final 
    271271    interval. 
    272272 
    … …  
    276276    if q_min < q[0]: 
    277277        if q_min <= 0: q_min = q[0]/10 
    278         delta = q[1] - q[0] 
    279         q_low = np.arange(q_min, q[0], delta) 
     278        n_low = np.ceil((q[0]-q_min) / (q[1]-q[0])) if q[1]>q[0] else 15 
     279        q_low = np.linspace(q_min, q[0], n_low+1)[:-1] 
    280280    else: 
    281281        q_low = [] 
    282282    if q_max > q[-1]: 
    283         delta = q[-1] - q[-2] 
    284         q_high = np.arange(q[-1]+delta, q_max, delta) 
     283        n_high = np.ceil((q_max-q[-1]) / (q[-1]-q[-2])) if q[-1]>q[-2] else 15 
     284        q_high = np.linspace(q[-1], q_max, n_high+1)[1:] 
    285285    else: 
    286286        q_high = [] 
    … …  
    288288 
    289289 
    290 def geometric_extrapolation(q, q_min, q_max): 
     290def geometric_extrapolation(q, q_min, q_max, points_per_decade=None): 
    291291    r""" 
    292292    Extrapolate *q* to [*q_min*, *q_max*] using geometric steps, with the 
    … …  
    295295    if *q_min* is zero or less then *q[0]/10* is used instead. 
    296296 
    297     Starting at $q_1$ and stepping geometrically by $\Delta q$ 
    298     to $q_n$ in $n$ points gives a geometric average of: 
     297    *points_per_decade* sets the ratio between consecutive steps such 
     298    that there will be $n$ points used for every factor of 10 increase 
     299    in *q*. 
     300 
     301    If *points_per_decade* is not given, it will be estimated as follows. 
     302    Starting at $q_1$ and stepping geometrically by $\Delta q$ to $q_n$ 
     303    in $n$ points gives a geometric average of: 
    299304 
    300305    .. math:: 
    … …  
    315320    """ 
    316321    q = np.sort(q) 
    317     delta_q = (len(q)-1)/(log(q[-1]) - log(q[0])) 
     322    if points_per_decade is None: 
     323        log_delta_q = (len(q) - 1) / (log(q[-1]) - log(q[0])) 
     324    else: 
     325        log_delta_q = log(10.) / points_per_decade 
    318326    if q_min < q[0]: 
    319327        if q_min < 0: q_min = q[0]/10 
    320         n_low = delta_q * (log(q[0])-log(q_min)) 
     328        n_low = log_delta_q * (log(q[0])-log(q_min)) 
    321329        q_low  = np.logspace(log10(q_min), log10(q[0]), np.ceil(n_low)+1)[:-1] 
    322330    else: 
    323331        q_low = [] 
    324332    if q_max > q[-1]: 
    325         n_high = delta_q * (log(q_max)-log(q[-1])) 
     333        n_high = log_delta_q * (log(q_max)-log(q[-1])) 
    326334        q_high = np.logspace(log10(q[-1]), log10(q_max), np.ceil(n_high)+1)[1:] 
    327335    else: 
    … …  
    329337    return np.concatenate([q_low, q, q_high]) 
    330338 
     339 
     340############################################################################ 
     341# unit tests 
     342############################################################################ 
     343import unittest 
     344 
     345 
     346def eval_form(q, form, pars): 
     347    from sasmodels import core 
     348    kernel = core.make_kernel(form, [q]) 
     349    theory = core.call_kernel(kernel, pars) 
     350    kernel.release() 
     351    return theory 
     352 
     353 
     354def gaussian(q, q0, dq): 
     355    from numpy import exp, pi 
     356    return exp(-0.5*((q-q0)/dq)**2)/(sqrt(2*pi)*dq) 
     357 
     358 
     359def romberg_slit_1d(q, delta_qv, form, pars): 
     360    """ 
     361    Romberg integration for slit resolution. 
     362 
     363    This is an adaptive integration technique.  It is called with settings 
     364    that make it slow to evaluate but give it good accuracy. 
     365    """ 
     366    from scipy.integrate import romberg 
     367 
     368    if any(k not in form.info['defaults'] for k in pars.keys()): 
     369        keys = set(form.info['defaults'].keys()) 
     370        extra = set(pars.keys()) - keys 
     371        raise ValueError("bad parameters: [%s] not in [%s]"% 
     372                         (", ".join(sorted(extra)), ", ".join(sorted(keys)))) 
     373 
     374    _fn = lambda u, q0: eval_form(sqrt(q0**2 + u**2), form, pars) 
     375    r = [romberg(_fn, 0, delta_qv[0], args=(qi,), 
     376                 divmax=100, vec_func=True, tol=0, rtol=1e-8) 
     377         for qi in q] 
     378    # r should be [float, ...], but it is [array([float]), array([float]),...] 
     379    return np.asarray(r).flatten()/delta_qv[0] 
     380 
     381 
     382def romberg_pinhole_1d(q, q_width, form, pars, nsigma=5): 
     383    """ 
     384    Romberg integration for pinhole resolution. 
     385 
     386    This is an adaptive integration technique.  It is called with settings 
     387    that make it slow to evaluate but give it good accuracy. 
     388    """ 
     389    from scipy.integrate import romberg 
     390 
     391    if any(k not in form.info['defaults'] for k in pars.keys()): 
     392        keys = set(form.info['defaults'].keys()) 
     393        extra = set(pars.keys()) - keys 
     394        raise ValueError("bad parameters: [%s] not in [%s]"% 
     395                         (", ".join(sorted(extra)), ", ".join(sorted(keys)))) 
     396 
     397    _fn = lambda q, q0, dq: eval_form(q, form, pars)*gaussian(q, q0, dq) 
     398    r = [romberg(_fn, max(qi-nsigma*dqi,1e-10*q[0]), qi+nsigma*dqi, args=(qi, dqi), 
     399                 divmax=100, vec_func=True, tol=0, rtol=1e-8) 
     400         for qi,dqi in zip(q,q_width)] 
     401    return np.asarray(r).flatten() 
     402 
     403 
     404class ResolutionTest(unittest.TestCase): 
     405 
     406    def setUp(self): 
     407        self.x = 0.001*np.arange(1, 11) 
     408        self.y = self.Iq(self.x) 
     409 
     410    def Iq(self, q): 
     411        "Linear function for resolution unit test" 
     412        return 12.0 - 1000.0*q 
     413 
     414    def test_perfect(self): 
     415        """ 
     416        Perfect resolution and no smearing. 
     417        """ 
     418        resolution = Perfect1D(self.x) 
     419        theory = self.Iq(resolution.q_calc) 
     420        output = resolution.apply(theory) 
     421        np.testing.assert_equal(output, self.y) 
     422 
     423    def test_slit_zero(self): 
     424        """ 
     425        Slit smearing with perfect resolution. 
     426        """ 
     427        resolution = Slit1D(self.x, width=0, height=0, q_calc=self.x) 
     428        theory = self.Iq(resolution.q_calc) 
     429        output = resolution.apply(theory) 
     430        np.testing.assert_equal(output, self.y) 
     431 
     432    @unittest.skip("not yet supported") 
     433    def test_slit_high(self): 
     434        """ 
     435        Slit smearing with height 0.005 
     436        """ 
     437        resolution = Slit1D(self.x, width=0, height=0.005, q_calc=self.x) 
     438        theory = self.Iq(resolution.q_calc) 
     439        output = resolution.apply(theory) 
     440        answer = [ 9.0618, 8.6402, 8.1187, 7.1392, 6.1528, 
     441                   5.5555, 4.5584, 3.5606, 2.5623, 2.0000 ] 
     442        np.testing.assert_allclose(output, answer, atol=1e-4) 
     443 
     444    @unittest.skip("not yet supported") 
     445    def test_slit_both_high(self): 
     446        """ 
     447        Slit smearing with width < 100*height. 
     448        """ 
     449        q = np.logspace(-4, -1, 10) 
     450        resolution = Slit1D(q, width=0.2, height=np.inf) 
     451        theory = 1000*self.Iq(resolution.q_calc**4) 
     452        output = resolution.apply(theory) 
     453        answer = [ 8.85785, 8.43012, 7.92687, 6.94566, 6.03660, 
     454                   5.40363, 4.40655, 3.40880, 2.41058, 2.00000 ] 
     455        np.testing.assert_allclose(output, answer, atol=1e-4) 
     456 
     457    @unittest.skip("not yet supported") 
     458    def test_slit_wide(self): 
     459        """ 
     460        Slit smearing with width 0.0002 
     461        """ 
     462        resolution = Slit1D(self.x, width=0.0002, height=0, q_calc=self.x) 
     463        theory = self.Iq(resolution.q_calc) 
     464        output = resolution.apply(theory) 
     465        answer = [ 11.0, 10.0, 9.0, 8.0, 7.0, 6.0, 5.0, 4.0, 3.0, 2.0 ] 
     466        np.testing.assert_allclose(output, answer, atol=1e-4) 
     467 
     468    @unittest.skip("not yet supported") 
     469    def test_slit_both_wide(self): 
     470        """ 
     471        Slit smearing with width > 100*height. 
     472        """ 
     473        resolution = Slit1D(self.x, width=0.0002, height=0.000001, 
     474                            q_calc=self.x) 
     475        theory = self.Iq(resolution.q_calc) 
     476        output = resolution.apply(theory) 
     477        answer = [ 11.0, 10.0, 9.0, 8.0, 7.0, 6.0, 5.0, 4.0, 3.0, 2.0 ] 
     478        np.testing.assert_allclose(output, answer, atol=1e-4) 
     479 
     480    def test_pinhole_zero(self): 
     481        """ 
     482        Pinhole smearing with perfect resolution 
     483        """ 
     484        resolution = Pinhole1D(self.x, 0.0*self.x) 
     485        theory = self.Iq(resolution.q_calc) 
     486        output = resolution.apply(theory) 
     487        np.testing.assert_equal(output, self.y) 
     488 
     489    def test_pinhole(self): 
     490        """ 
     491        Pinhole smearing with dQ = 0.001 [Note: not dQ/Q = 0.001] 
     492        """ 
     493        resolution = Pinhole1D(self.x, 0.001*np.ones_like(self.x), 
     494                               q_calc=self.x) 
     495        theory = 12.0-1000.0*resolution.q_calc 
     496        output = resolution.apply(theory) 
     497        answer = [ 10.44785079, 9.84991299, 8.98101708, 
     498                  7.99906585, 6.99998311, 6.00001689, 
     499                  5.00093415, 4.01898292, 3.15008701, 2.55214921] 
     500        np.testing.assert_allclose(output, answer, atol=1e-8) 
     501 
     502 
     503class IgorComparisonTest(unittest.TestCase): 
     504 
     505    def setUp(self): 
     506        self.pars = TEST_PARS_PINHOLE_SPHERE 
     507        from sasmodels import core 
     508        from sasmodels.models import sphere 
     509        self.model = core.load_model(sphere, dtype='double') 
     510 
     511    def Iq_sphere(self, pars, resolution): 
     512        from sasmodels import core 
     513        kernel = core.make_kernel(self.model, [resolution.q_calc]) 
     514        theory = core.call_kernel(kernel, pars) 
     515        result = resolution.apply(theory) 
     516        kernel.release() 
     517        return result 
     518 
     519    def compare(self, q, output, answer, tolerance): 
     520        err = (output - answer)/answer 
     521        idx = abs(err) >= tolerance 
     522        problem = zip(q[idx], output[idx], answer[idx], err[idx]) 
     523        print "\n".join(str(v) for v in problem) 
     524        np.testing.assert_allclose(output, answer, rtol=tolerance) 
     525 
     526    def test_perfect(self): 
     527        """ 
     528        Compare sphere model with NIST Igor SANS, no resolution smearing. 
     529        """ 
     530        pars = TEST_PARS_SLIT_SPHERE 
     531        data_string = TEST_DATA_SLIT_SPHERE 
     532 
     533        data = np.loadtxt(data_string.split('\n')).T 
     534        q, width, answer, _ = data 
     535        resolution = Perfect1D(q) 
     536        output = self.Iq_sphere(pars, resolution) 
     537        self.compare(q, output, answer, 1e-6) 
     538 
     539    def test_pinhole(self): 
     540        """ 
     541        Compare pinhole resolution smearing with NIST Igor SANS 
     542        """ 
     543        pars = TEST_PARS_PINHOLE_SPHERE 
     544        data_string = TEST_DATA_PINHOLE_SPHERE 
     545 
     546        data = np.loadtxt(data_string.split('\n')).T 
     547        q, q_width, answer = data 
     548        resolution = Pinhole1D(q, q_width) 
     549        output = self.Iq_sphere(pars, resolution) 
     550        # TODO: relative error should be lower 
     551        self.compare(q, output, answer, 3e-4) 
     552 
     553    def test_pinhole_romberg(self): 
     554        """ 
     555        Compare pinhole resolution smearing with romberg integration result. 
     556        """ 
     557        pars = TEST_PARS_PINHOLE_SPHERE 
     558        data_string = TEST_DATA_PINHOLE_SPHERE 
     559        pars['radius'] *= 5 
     560        radius = pars['radius'] 
     561 
     562        data = np.loadtxt(data_string.split('\n')).T 
     563        q, q_width, answer = data 
     564        answer = romberg_pinhole_1d(q, q_width, self.model, pars) 
     565        ## Getting 0.1% requires 5 sigma and 200 points per fringe 
     566        #q_calc = interpolate(pinhole_extend_q(q, q_width, nsigma=5), 
     567        #                     2*np.pi/radius/200) 
     568        #tol = 0.001 
     569        ## The default 3 sigma and no extra points gets 1% 
     570        q_calc, tol = None, 0.01 
     571        resolution = Pinhole1D(q, q_width, q_calc=q_calc) 
     572        output = self.Iq_sphere(pars, resolution) 
     573        if 0: # debug plot 
     574            import matplotlib.pyplot as plt 
     575            resolution = Perfect1D(q) 
     576            source = self.Iq_sphere(pars, resolution) 
     577            plt.loglog(q, source, '.') 
     578            plt.loglog(q, answer, '-', hold=True) 
     579            plt.loglog(q, output, '-', hold=True) 
     580            plt.show() 
     581        self.compare(q, output, answer, tol) 
     582 
     583    def test_slit(self): 
     584        """ 
     585        Compare slit resolution smearing with NIST Igor SANS 
     586        """ 
     587        pars = TEST_PARS_SLIT_SPHERE 
     588        data_string = TEST_DATA_SLIT_SPHERE 
     589 
     590        data = np.loadtxt(data_string.split('\n')).T 
     591        q, delta_qv, _, answer = data 
     592        resolution = Slit1D(q, width=delta_qv, height=0) 
     593        output = self.Iq_sphere(pars, resolution) 
     594        # TODO: eliminate Igor test since it is too inaccurate to be useful. 
     595        # This means we can eliminate the test data as well, and instead 
     596        # use a generated q vector. 
     597        self.compare(q, output, answer, 0.5) 
     598 
     599    def test_slit_romberg(self): 
     600        """ 
     601        Compare slit resolution smearing with romberg integration result. 
     602        """ 
     603        pars = TEST_PARS_SLIT_SPHERE 
     604        data_string = TEST_DATA_SLIT_SPHERE 
     605        radius = pars['radius'] 
     606 
     607        data = np.loadtxt(data_string.split('\n')).T 
     608        q, delta_qv, _, answer = data 
     609        answer = romberg_slit_1d(q, delta_qv, self.model, pars) 
     610        q_calc = slit_extend_q(interpolate(q, 2*np.pi/radius/20), 
     611                               delta_qv[0], delta_qv[0]) 
     612        resolution = Slit1D(q, width=delta_qv, height=0, q_calc=q_calc) 
     613        output = self.Iq_sphere(pars, resolution) 
     614        # TODO: relative error should be lower 
     615        self.compare(q, output, answer, 0.025) 
     616 
     617    def test_ellipsoid(self): 
     618        """ 
     619        Compare romberg integration for ellipsoid model. 
     620        """ 
     621        from .core import load_model 
     622        pars = { 
     623            'scale':0.05, 
     624            'rpolar':500, 'requatorial':15000, 
     625            'sld':6, 'solvent_sld': 1, 
     626            } 
     627        form = load_model('ellipsoid', dtype='double') 
     628        q = np.logspace(log10(4e-5),log10(2.5e-2), 68) 
     629        delta_qv = [0.117] 
     630        resolution = Slit1D(q, width=delta_qv, height=0) 
     631        answer = romberg_slit_1d(q, delta_qv, form, pars) 
     632        output = resolution.apply(eval_form(resolution.q_calc, form, pars)) 
     633        # TODO: 10% is too much error; use better algorithm 
     634        self.compare(q, output, answer, 0.1) 
     635 
     636    #TODO: can sas q spacing be too sparse for the resolution calculation? 
     637    @unittest.skip("suppress sparse data test; not supported by current code") 
     638    def test_pinhole_sparse(self): 
     639        """ 
     640        Compare pinhole resolution smearing with NIST Igor SANS on sparse data 
     641        """ 
     642        pars = TEST_PARS_PINHOLE_SPHERE 
     643        data_string = TEST_DATA_PINHOLE_SPHERE 
     644 
     645        data = np.loadtxt(data_string.split('\n')).T 
     646        q, q_width, answer = data[:, ::20] # Take every nth point 
     647        resolution = Pinhole1D(q, q_width) 
     648        output = self.Iq_sphere(pars, resolution) 
     649        self.compare(q, output, answer, 1e-6) 
     650 
     651 
     652# pinhole sphere parameters 
     653TEST_PARS_PINHOLE_SPHERE = { 
     654    'scale': 1.0, 'background': 0.01, 
     655    'radius': 60.0, 'sld': 1, 'solvent_sld': 6.3, 
     656    } 
     657# Q, dQ, I(Q) calculated by NIST Igor SANS package 
     658TEST_DATA_PINHOLE_SPHERE = """\ 
     6590.001278 0.0002847 2538.41176383 
     6600.001562 0.0002905 2536.91820405 
     6610.001846 0.0002956 2535.13182479 
     6620.002130 0.0003017 2533.06217813 
     6630.002414 0.0003087 2530.70378586 
     6640.002698 0.0003165 2528.05024192 
     6650.002982 0.0003249 2525.10408349 
     6660.003266 0.0003340 2521.86667499 
     6670.003550 0.0003437 2518.33907750 
     6680.003834 0.0003539 2514.52246995 
     6690.004118 0.0003646 2510.41798319 
     6700.004402 0.0003757 2506.02690988 
     6710.004686 0.0003872 2501.35067884 
     6720.004970 0.0003990 2496.38678318 
     6730.005253 0.0004112 2491.16237596 
     6740.005537 0.0004237 2485.63911673 
     6750.005821 0.0004365 2479.83657083 
     6760.006105 0.0004495 2473.75676948 
     6770.006389 0.0004628 2467.40145990 
     6780.006673 0.0004762 2460.77293372 
     6790.006957 0.0004899 2453.86724627 
     6800.007241 0.0005037 2446.69623838 
     6810.007525 0.0005177 2439.25775219 
     6820.007809 0.0005318 2431.55421398 
     6830.008093 0.0005461 2423.58785521 
     6840.008377 0.0005605 2415.36158137 
     6850.008661 0.0005750 2406.87009473 
     6860.008945 0.0005896 2398.12841186 
     6870.009229 0.0006044 2389.13360806 
     6880.009513 0.0006192 2379.88958042 
     6890.009797 0.0006341 2370.39776774 
     6900.010080 0.0006491 2360.69528793 
     6910.010360 0.0006641 2350.85169027 
     6920.010650 0.0006793 2340.42023633 
     6930.010930 0.0006945 2330.11206013 
     6940.011220 0.0007097 2319.20109972 
     6950.011500 0.0007251 2308.43503981 
     6960.011780 0.0007404 2297.44820179 
     6970.012070 0.0007558 2285.83853677 
     6980.012350 0.0007713 2274.41290746 
     6990.012640 0.0007868 2262.36219581 
     7000.012920 0.0008024 2250.51169731 
     7010.013200 0.0008180 2238.45596231 
     7020.013490 0.0008336 2225.76495666 
     7030.013770 0.0008493 2213.29618391 
     7040.014060 0.0008650 2200.19110751 
     7050.014340 0.0008807 2187.34050325 
     7060.014620 0.0008965 2174.30529864 
     7070.014910 0.0009123 2160.61632548 
     7080.015190 0.0009281 2147.21038112 
     7090.015470 0.0009440 2133.62023580 
     7100.015760 0.0009598 2119.37907426 
     7110.016040 0.0009757 2105.45234903 
     7120.016330 0.0009916 2090.86319102 
     7130.016610 0.0010080 2076.60576032 
     7140.016890 0.0010240 2062.19214565 
     7150.017180 0.0010390 2047.10550219 
     7160.017460 0.0010550 2032.38715621 
     7170.017740 0.0010710 2017.52560123 
     7180.018030 0.0010880 2001.99124318 
     7190.018310 0.0011040 1986.84662060 
     7200.018600 0.0011200 1971.03389745 
     7210.018880 0.0011360 1955.61395119 
     7220.019160 0.0011520 1940.08291563 
     7230.019450 0.0011680 1923.87672225 
     7240.019730 0.0011840 1908.10656374 
     7250.020020 0.0012000 1891.66297192 
     7260.020300 0.0012160 1875.66789021 
     7270.020580 0.0012320 1859.56357196 
     7280.020870 0.0012490 1842.79468290 
     7290.021150 0.0012650 1826.50064489 
     7300.021430 0.0012810 1810.11533702 
     7310.021720 0.0012970 1793.06840882 
     7320.022000 0.0013130 1776.51153580 
     7330.022280 0.0013290 1759.87201249 
     7340.022570 0.0013460 1742.57354412 
     7350.022850 0.0013620 1725.79397319 
     7360.023140 0.0013780 1708.35831550 
     7370.023420 0.0013940 1691.45256069 
     7380.023700 0.0014110 1674.48561783 
     7390.023990 0.0014270 1656.86525366 
     7400.024270 0.0014430 1639.79847285 
     7410.024550 0.0014590 1622.68887088 
     7420.024840 0.0014760 1604.96421100 
     7430.025120 0.0014920 1587.85768129 
     7440.025410 0.0015080 1569.99297335 
     7450.025690 0.0015240 1552.84580279 
     7460.025970 0.0015410 1535.54074115 
     7470.026260 0.0015570 1517.75249337 
     7480.026540 0.0015730 1500.40115023 
     7490.026820 0.0015900 1483.03632237 
     7500.027110 0.0016060 1465.05942429 
     7510.027390 0.0016220 1447.67682181 
     7520.027670 0.0016390 1430.46495191 
     7530.027960 0.0016550 1412.49232282 
     7540.028240 0.0016710 1395.13182318 
     7550.028520 0.0016880 1377.93439837 
     7560.028810 0.0017040 1359.99528971 
     7570.029090 0.0017200 1342.67274512 
     7580.029370 0.0017370 1325.55375609 
     759""" 
     760 
     761# Slit sphere parameters 
     762TEST_PARS_SLIT_SPHERE = { 
     763    'scale': 0.01, 'background': 0.01, 
     764    'radius': 60000, 'sld': 1, 'solvent_sld': 4, 
     765    } 
     766# Q dQ I(Q) I_smeared(Q) 
     767TEST_DATA_SLIT_SPHERE = """\ 
     7682.26097e-05 0.117 5.5781372896e+09 1.4626077708e+06 
     7692.53847e-05 0.117 5.0363141458e+09 1.3117318023e+06 
     7702.81597e-05 0.117 4.4850108103e+09 1.1594863713e+06 
     7713.09347e-05 0.117 3.9364658459e+09 1.0094881630e+06 
     7723.37097e-05 0.117 3.4019975074e+09 8.6518941303e+05 
     7733.92597e-05 0.117 2.4139519814e+09 6.0232158311e+05 
     7744.48097e-05 0.117 1.5816877820e+09 3.8739994090e+05 
     7755.03597e-05 0.117 9.3715407224e+08 2.2745304775e+05 
     7765.59097e-05 0.117 4.8387917428e+08 1.2101295768e+05 
     7776.14597e-05 0.117 2.0193586928e+08 6.0055107771e+04 
     7786.70097e-05 0.117 5.5886110911e+07 3.2749521065e+04 
     7797.25597e-05 0.117 3.7782348010e+06 2.6350963616e+04 
     7807.81097e-05 0.117 5.3407817904e+06 2.9624963314e+04 
     7818.36597e-05 0.117 2.7975485523e+07 3.4403962254e+04 
     7828.92097e-05 0.117 4.9845448282e+07 3.6130017903e+04 
     7839.47597e-05 0.117 6.0092588905e+07 3.3495107849e+04 
     7841.00310e-04 0.117 5.6823430831e+07 2.7475726279e+04 
     7851.05860e-04 0.117 4.3857024036e+07 2.0144282226e+04 
     7861.11410e-04 0.117 2.7277144760e+07 1.3647403260e+04 
     7871.22510e-04 0.117 3.3119334113e+06 6.6519711526e+03 
     7881.33610e-04 0.117 1.4412859402e+06 6.9726212813e+03 
     7891.44710e-04 0.117 8.5620162463e+06 8.1441335775e+03 
     7901.55810e-04 0.117 9.6957429033e+06 6.4559996521e+03 
     7911.66910e-04 0.117 4.3818341914e+06 3.6252154156e+03 
     7921.78010e-04 0.117 2.7448997387e+05 2.4006505342e+03 
     7931.89110e-04 0.117 8.0472009936e+05 2.8187789089e+03 
     7942.00210e-04 0.117 2.8149552834e+06 3.0915662855e+03 
     7952.11310e-04 0.117 2.7510907861e+06 2.3722530293e+03 
     7962.22410e-04 0.117 1.0053133293e+06 1.4473468311e+03 
     7972.33510e-04 0.117 5.8428305052e+03 1.2048540556e+03 
     7982.44610e-04 0.117 5.1699305004e+05 1.4625670042e+03 
     7992.55710e-04 0.117 1.2120227268e+06 1.5010705973e+03 
     8002.66810e-04 0.117 9.7896842846e+05 1.1336343426e+03 
     8012.77910e-04 0.117 2.5507264791e+05 8.1848818080e+02 
     8023.05660e-04 0.117 5.2403101181e+05 7.4913374357e+02 
     8033.33410e-04 0.117 5.8699343809e+04 4.4669964560e+02 
     8043.61160e-04 0.117 3.0844327150e+05 4.6774007542e+02 
     8053.88910e-04 0.117 8.3360142970e+03 2.7169550220e+02 
     8064.16660e-04 0.117 1.8630080583e+05 3.0710983679e+02 
     8074.44410e-04 0.117 3.1616804732e-01 1.7959006831e+02 
     8084.72160e-04 0.117 1.1299016314e+05 2.0763952339e+02 
     8094.99910e-04 0.117 2.9952522747e+03 1.2536542765e+02 
     8105.27660e-04 0.117 6.7625695649e+04 1.4013969777e+02 
     8115.55410e-04 0.117 7.6927460089e+03 8.2145593180e+01 
     8126.10910e-04 0.117 1.1229057779e+04 8.4519745643e+01 
     8136.66410e-04 0.117 1.3035567943e+04 8.1554625609e+01 
     8147.21910e-04 0.117 1.3309931343e+04 7.4437319172e+01 
     8157.77410e-04 0.117 1.2462626212e+04 6.4697088261e+01 
     8168.32910e-04 0.117 1.0912927143e+04 5.3773301044e+01 
     8178.88410e-04 0.117 9.0172597469e+03 4.2843375753e+01 
     8189.43910e-04 0.117 7.0496495917e+03 3.2771032724e+01 
     8199.99410e-04 0.117 5.2030483682e+03 2.4113557144e+01 
     8201.05491e-03 0.117 3.5988976711e+03 1.7160773658e+01 
     8211.11041e-03 0.117 2.2996060652e+03 1.2016626459e+01 
     8221.22141e-03 0.117 6.4766590598e+02 6.0373017740e+00 
     8231.33241e-03 0.117 4.1963483264e+01 4.5215452974e+00 
     8241.44341e-03 0.117 6.3370708246e+01 5.1054681903e+00 
     8251.55441e-03 0.117 3.0736750577e+02 5.9176165298e+00 
     8261.66541e-03 0.117 5.0327682399e+02 5.9815000189e+00 
     8271.77641e-03 0.117 5.4084331454e+02 5.1634639625e+00 
     8281.88741e-03 0.117 4.3488671756e+02 3.8535158148e+00 
     8291.99841e-03 0.117 2.6322287860e+02 2.5824997753e+00 
     8302.10941e-03 0.117 1.0793633150e+02 1.7315517194e+00 
     8312.22041e-03 0.117 1.8474448850e+01 1.4077213604e+00 
     8322.33141e-03 0.117 1.5864062279e+00 1.4771560682e+00 
     8332.44241e-03 0.117 3.2267213848e+01 1.6916253448e+00 
     8342.55341e-03 0.117 7.4289116207e+01 1.8274751193e+00 
     8352.66441e-03 0.117 9.9000521929e+01 1.7706812289e+00 
     836""" 
     837 
     838def main(): 
     839    """ 
     840    Run tests given is sys.argv. 
     841 
     842    Returns 0 if success or 1 if any tests fail. 
     843    """ 
     844    import sys 
     845    import xmlrunner 
     846 
     847    suite = unittest.TestSuite() 
     848    suite.addTest(unittest.defaultTestLoader.loadTestsFromModule(sys.modules[__name__])) 
     849 
     850    runner = xmlrunner.XMLTestRunner(output='logs') 
     851    result = runner.run(suite) 
     852    return 1 if result.failures or result.errors else 0 
     853 
     854 
     855############################################################################ 
     856# usage demo 
     857############################################################################ 
     858 
     859def _eval_demo_1d(resolution, title): 
     860    from sasmodels import core 
     861    from sasmodels.models import cylinder 
     862    ## or alternatively: 
     863    # cylinder = core.load_model_definition('cylinder') 
     864    model = core.load_model(cylinder) 
     865 
     866    kernel = core.make_kernel(model, [resolution.q_calc]) 
     867    theory = core.call_kernel(kernel, {'length':210, 'radius':500}) 
     868    Iq = resolution.apply(theory) 
     869 
     870    import matplotlib.pyplot as plt 
     871    plt.loglog(resolution.q_calc, theory, label='unsmeared') 
     872    plt.loglog(resolution.q, Iq, label='smeared', hold=True) 
     873    plt.legend() 
     874    plt.title(title) 
     875    plt.xlabel("Q (1/Ang)") 
     876    plt.ylabel("I(Q) (1/cm)") 
     877 
     878def demo_pinhole_1d(): 
     879    q = np.logspace(-3, -1, 400) 
     880    q_width = 0.1*q 
     881    resolution = Pinhole1D(q, q_width) 
     882    _eval_demo_1d(resolution, title="10% dQ/Q Pinhole Resolution") 
     883 
     884def demo_slit_1d(): 
     885    q = np.logspace(-3, -1, 400) 
     886    qx_width = 0.005 
     887    qy_width = 0.0 
     888    resolution = Slit1D(q, qx_width, qy_width) 
     889    _eval_demo_1d(resolution, title="0.005 Qx Slit Resolution") 
     890 
     891def demo(): 
     892    import matplotlib.pyplot as plt 
     893    plt.subplot(121) 
     894    demo_pinhole_1d() 
     895    plt.subplot(122) 
     896    demo_slit_1d() 
     897    plt.show() 
     898 
     899 
     900if __name__ == "__main__": 
     901    #demo() 
     902    main() 

  • src/sas/perspectives/fitting/fitting.py

    rf76bf17 ra3f125f0  
    954954                return True 
    955955            except: 
     956                raise 
    956957                msg = "Fitting error: %s" % str(sys.exc_value) 
    957958                wx.PostEvent(self.parent, StatusEvent(status=msg, info="error", 

  • src/sas/perspectives/fitting/model_thread.py

    r2f4b430 ra3f125f0  
    169169            first_bin, last_bin = self.smearer.get_bin_range(self.qmin, 
    170170                                                             self.qmax) 
    171             mask = self.data.x[first_bin:last_bin] 
    172             output[first_bin:last_bin] = self.model.evalDistribution(mask) 
     171            mask = self.data.x[first_bin:last_bin+1] 
     172            output[first_bin:last_bin+1] = self.model.evalDistribution(mask) 
    173173            output = self.smearer(output, first_bin, last_bin) 
    174174        else: 
Note: See TracChangeset for help on using the changeset viewer.