Changeset ea75043 in sasmodels


Ignore:
Timestamp:
Mar 29, 2016 6:22:13 PM (9 years ago)
Author:
Paul Kienzle <pkienzle@…>
Branches:
master, core_shell_microgels, costrafo411, magnetic_model, release_v0.94, release_v0.95, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests
Children:
d6f5da6
Parents:
1d61d07
Message:

add support for oriented usans using direct 2d resolution integral

Files:
1 added
6 edited

Legend:

Unmodified
Added
Removed

  • sasmodels/bumps_model.py

    r9217ef8 rea75043  
    223223        """ 
    224224        data, theory, resid = self._data, self.theory(), self.residuals() 
    225         plot_theory(data, theory, resid, view) 
     225        plot_theory(data, theory, resid, view, Iq_calc = self.Iq_calc) 
    226226 
    227227    def simulate_data(self, noise=None): 

  • sasmodels/compare.py

    rfa1582e rea75043  
    473473        data = empty_data2D(np.linspace(-qmax, qmax, nq), resolution=res) 
    474474        data.accuracy = opts['accuracy'] 
    475         set_beam_stop(data, 0.004) 
     475        set_beam_stop(data, 0.0004) 
    476476        index = ~data.mask 
    477477    else: 

  • sasmodels/data.py

    r1d61d07 rea75043  
    322322 
    323323def plot_theory(data, theory, resid=None, view='log', 
    324                 use_data=True, limits=None): 
     324                use_data=True, limits=None, Iq_calc=None): 
    325325    """ 
    326326    Plot theory calculation. 
    … …  
    343343        _plot_result2D(data, theory, resid, view, use_data, limits=limits) 
    344344    else: 
    345         _plot_result1D(data, theory, resid, view, use_data, limits=limits) 
     345        _plot_result1D(data, theory, resid, view, use_data, 
     346                       limits=limits, Iq_calc=Iq_calc) 
    346347 
    347348 
    … …  
    366367 
    367368@protect 
    368 def _plot_result1D(data, theory, resid, view, use_data, limits=None): 
     369def _plot_result1D(data, theory, resid, view, use_data, 
     370                   limits=None, Iq_calc=None): 
    369371    """ 
    370372    Plot the data and residuals for 1D data. 
    … …  
    376378    use_theory = theory is not None 
    377379    use_resid = resid is not None 
    378     num_plots = (use_data or use_theory) + use_resid 
     380    use_calc = use_theory and Iq_calc is not None 
     381    num_plots = (use_data or use_theory) + use_calc + use_resid 
     382 
    379383 
    380384    scale = data.x**4 if view == 'q4' else 1.0 
    … …  
    416420        plt.ylabel('$I(q)$') 
    417421 
     422    if use_calc: 
     423        # Only have use_calc if have use_theory 
     424        plt.subplot(1, num_plots, 2) 
     425        qx, qy, Iqxy = Iq_calc 
     426        plt.pcolormesh(qx, qy, np.log10(Iqxy)) 
     427        plt.xlabel("$q_x$/A$^{-1}$") 
     428        plt.xlabel("$q_y$/A$^{-1}$") 
     429        #plt.axis('equal') 
     430 
    418431    if use_resid: 
    419432        mresid = masked_array(resid, data.mask.copy()) 
    … …  
    422435 
    423436        if num_plots > 1: 
    424             plt.subplot(1, num_plots, (use_data or use_theory) + 1) 
     437            plt.subplot(1, num_plots, use_calc + 2) 
    425438        plt.plot(data.x, mresid, '-') 
    426439        plt.xlabel("$q$/A$^{-1}$") 
    … …  
    562575    plottable = masked_array(image, ~valid | data.mask) 
    563576    # Divide range by 10 to convert from angstroms to nanometers 
    564     xmin, xmax = min(data.qx_data)/10, max(data.qx_data)/10 
    565     ymin, ymax = min(data.qy_data)/10, max(data.qy_data)/10 
     577    xmin, xmax = min(data.qx_data), max(data.qx_data) 
     578    ymin, ymax = min(data.qy_data), max(data.qy_data) 
    566579    if view == 'log': 
    567580        vmin, vmax = np.log10(vmin), np.log10(vmax) 
    568581    plt.imshow(plottable.reshape(len(data.x_bins), len(data.y_bins)), 
    569                interpolation='nearest', aspect=1, origin='upper', 
     582               interpolation='nearest', aspect=1, origin='lower', 
    570583               extent=[xmin, xmax, ymin, ymax], vmin=vmin, vmax=vmax) 
    571     plt.xlabel("$q_x$/nm$^{-1}$") 
    572     plt.ylabel("$q_y$/nm$^{-1}$") 
     584    plt.xlabel("$q_x$/A$^{-1}$") 
     585    plt.ylabel("$q_y$/A$^{-1}$") 
    573586    return vmin, vmax 
    574587 

  • sasmodels/direct_model.py

    r02e70ff rea75043  
    6363        elif hasattr(data, 'qx_data'): 
    6464            self.data_type = 'Iqxy' 
     65        elif getattr(data, 'oriented', False): 
     66            self.data_type = 'Iq-oriented' 
    6567        else: 
    6668            self.data_type = 'Iq' 
    … …  
    116118                  and getattr(data, 'dxw', None) is not None): 
    117119                res = resolution.Slit1D(data.x[index], 
    118                                         width=data.dxh[index], 
    119                                         height=data.dxw[index]) 
     120                                        qx_width=data.dxw[index], 
     121                                        qy_width=data.dxl[index]) 
    120122            else: 
    121123                res = resolution.Perfect1D(data.x[index]) 
    … …  
    123125            #self._theory = np.zeros_like(self.Iq) 
    124126            q_vectors = [res.q_calc] 
     127            q_mono = [] 
     128        elif self.data_type == 'Iq-oriented': 
     129            index = (data.x >= data.qmin) & (data.x <= data.qmax) 
     130            if data.y is not None: 
     131                index &= ~np.isnan(data.y) 
     132                Iq = data.y[index] 
     133                dIq = data.dy[index] 
     134            else: 
     135                Iq, dIq = None, None 
     136            if (getattr(data, 'dxl', None) is None 
     137                or getattr(data, 'dxw', None) is None): 
     138                raise ValueError("oriented sample with 1D data needs slit resolution") 
     139 
     140            res = resolution2d.Slit2D(data.x[index], 
     141                                      qx_width=data.dxw[index], 
     142                                      qy_width=data.dxl[index]) 
     143            q_vectors = res.q_calc 
    125144            q_mono = [] 
    126145        else: 
    … …  
    142161        y = Iq + np.random.randn(*dy.shape) * dy 
    143162        self.Iq = y 
    144         if self.data_type == 'Iq': 
     163        if self.data_type in ('Iq', 'Iq-oriented'): 
    145164            self._data.dy[self.index] = dy 
    146165            self._data.y[self.index] = y 
    … …  
    155174        if self._kernel is None: 
    156175            self._kernel = make_kernel(self._model, self._kernel_inputs)  # pylint: disable=attribute-dedata_type 
    157             self._kernel_mono = make_kernel(self._model, self._kernel_mono_inputs) if self._kernel_mono_inputs else None 
     176            self._kernel_mono = (make_kernel(self._model, self._kernel_mono_inputs) 
     177                                 if self._kernel_mono_inputs else None) 
    158178 
    159179        Iq_calc = call_kernel(self._kernel, pars, cutoff=cutoff) 
    160         Iq_mono = call_kernel(self._kernel_mono, pars, mono=True) if self._kernel_mono_inputs else None 
     180        # TODO: may want to plot the raw Iq for other than oriented usans 
     181        self.Iq_calc = None 
    161182        if self.data_type == 'sesans': 
     183            Iq_mono = (call_kernel(self._kernel_mono, pars, mono=True) 
     184                       if self._kernel_mono_inputs else None) 
    162185            result = sesans.transform(self._data, 
    163186                                   self._kernel_inputs[0], Iq_calc,  
    … …  
    165188        else: 
    166189            result = self.resolution.apply(Iq_calc) 
     190            if hasattr(self.resolution, 'nx'): 
     191                self.Iq_calc = ( 
     192                    self.resolution.qx_calc, self.resolution.qy_calc, 
     193                    np.reshape(Iq_calc, (self.resolution.ny, self.resolution.nx)) 
     194                ) 
    167195        return result         
    168196 

  • sasmodels/resolution.py

    r8b935d1 rea75043  
    8282        self.q_calc = (pinhole_extend_q(q, q_width, nsigma=nsigma) 
    8383                       if q_calc is None else np.sort(q_calc)) 
    84         self.weight_matrix = pinhole_resolution( 
    85             self.q_calc, self.q, np.maximum(q_width, MINIMUM_RESOLUTION)) 
     84        self.weight_matrix = pinhole_resolution(self.q_calc, self.q, 
     85                                np.maximum(q_width, MINIMUM_RESOLUTION)) 
    8686 
    8787    def apply(self, theory): 
    … …  
    9191class Slit1D(Resolution): 
    9292    """ 
    93     Slit aperture with a complicated resolution function. 
     93    Slit aperture with resolution function. 
    9494 
    9595    *q* points at which the data is measured. 
    9696 
    97     *qx_width* slit width 
    98  
    99     *qy_height* slit height 
     97    *dqx* slit width in qx 
     98 
     99    *dqy* slit height in qy 
    100100 
    101101    *q_calc* is the list of points to calculate, or None if this should 
    … …  
    104104    The *weight_matrix* is computed by :func:`slit1d_resolution` 
    105105    """ 
    106     def __init__(self, q, width, height=0., q_calc=None): 
    107         # Remember what width/height was used even though we won't need them 
     106    def __init__(self, q, qx_width, qy_width=0., q_calc=None): 
     107        # Remember what width/dqy was used even though we won't need them 
    108108        # after the weight matrix is constructed 
    109         self.width, self.height = width, height 
     109        self.qx_width, self.qy_width = qx_width, qy_width 
    110110 
    111111        # Allow independent resolution on each point even though it is not 
    112112        # needed in practice. 
    113         if np.isscalar(width): 
    114             width = np.ones(len(q))*width 
     113        if np.isscalar(qx_width): 
     114            qx_width = np.ones(len(q))*qx_width 
    115115        else: 
    116             width = np.asarray(width) 
    117         if np.isscalar(height): 
    118             height = np.ones(len(q))*height 
     116            qx_width = np.asarray(qx_width) 
     117        if np.isscalar(qy_width): 
     118            qy_width = np.ones(len(q))*qy_width 
    119119        else: 
    120             height = np.asarray(height) 
     120            qy_width = np.asarray(qy_width) 
    121121 
    122122        self.q = q.flatten() 
    123         self.q_calc = slit_extend_q(q, width, height) \ 
     123        self.q_calc = slit_extend_q(q, qx_width, qy_width) \ 
    124124            if q_calc is None else np.sort(q_calc) 
    125125        self.weight_matrix = \ 
    126             slit_resolution(self.q_calc, self.q, width, height) 
     126            slit_resolution(self.q_calc, self.q, qx_width, qy_width) 
    127127 
    128128    def apply(self, theory): 
    … …  
    396396    """ 
    397397    q = np.sort(q) 
    398     if q_min < q[0]: 
     398    if q_min + 2*MINIMUM_RESOLUTION < q[0]: 
    399399        if q_min <= 0: q_min = q_min*MIN_Q_SCALE_FOR_NEGATIVE_Q_EXTRAPOLATION 
    400400        n_low = np.ceil((q[0]-q_min) / (q[1]-q[0])) if q[1] > q[0] else 15 
    … …  
    402402    else: 
    403403        q_low = [] 
    404     if q_max > q[-1]: 
     404    if q_max - 2*MINIMUM_RESOLUTION > q[-1]: 
    405405        n_high = np.ceil((q_max-q[-1]) / (q[-1]-q[-2])) if q[-1] > q[-2] else 15 
    406406        q_high = np.linspace(q[-1], q_max, n_high+1)[1:] 
    … …  
    595595        Slit smearing with perfect resolution. 
    596596        """ 
    597         resolution = Slit1D(self.x, width=0, height=0, q_calc=self.x) 
     597        resolution = Slit1D(self.x, qx_width=0, qy_width=0, q_calc=self.x) 
    598598        theory = self.Iq(resolution.q_calc) 
    599599        output = resolution.apply(theory) 
    … …  
    605605        Slit smearing with height 0.005 
    606606        """ 
    607         resolution = Slit1D(self.x, width=0, height=0.005, q_calc=self.x) 
     607        resolution = Slit1D(self.x, qx_width=0, qy_width=0.005, q_calc=self.x) 
    608608        theory = self.Iq(resolution.q_calc) 
    609609        output = resolution.apply(theory) 
    … …  
    620620        """ 
    621621        q = np.logspace(-4, -1, 10) 
    622         resolution = Slit1D(q, width=0.2, height=np.inf) 
     622        resolution = Slit1D(q, qx_width=0.2, qy_width=np.inf) 
    623623        theory = 1000*self.Iq(resolution.q_calc**4) 
    624624        output = resolution.apply(theory) 
    … …  
    634634        Slit smearing with width 0.0002 
    635635        """ 
    636         resolution = Slit1D(self.x, width=0.0002, height=0, q_calc=self.x) 
     636        resolution = Slit1D(self.x, qx_width=0.0002, qy_width=0, q_calc=self.x) 
    637637        theory = self.Iq(resolution.q_calc) 
    638638        output = resolution.apply(theory) 
    … …  
    647647        Slit smearing with width > 100*height. 
    648648        """ 
    649         resolution = Slit1D(self.x, width=0.0002, height=0.000001, 
     649        resolution = Slit1D(self.x, qx_width=0.0002, qy_width=0.000001, 
    650650                            q_calc=self.x) 
    651651        theory = self.Iq(resolution.q_calc) 
    … …  
    772772        data = np.loadtxt(data_string.split('\n')).T 
    773773        q, delta_qv, _, answer = data 
    774         resolution = Slit1D(q, width=delta_qv, height=0) 
     774        resolution = Slit1D(q, qx_width=delta_qv, qy_width=0) 
    775775        output = self._eval_sphere(pars, resolution) 
    776776        # TODO: eliminate Igor test since it is too inaccurate to be useful. 
    … …  
    792792        q_calc = slit_extend_q(interpolate(q, 2*np.pi/radius/20), 
    793793                               delta_qv[0], 0.) 
    794         resolution = Slit1D(q, width=delta_qv, height=0, q_calc=q_calc) 
     794        resolution = Slit1D(q, qx_width=delta_qv, qy_width=0, q_calc=q_calc) 
    795795        output = self._eval_sphere(pars, resolution) 
    796796        # TODO: relative error should be lower 
    … …  
    810810        q = np.logspace(log10(4e-5), log10(2.5e-2), 68) 
    811811        width, height = 0.117, 0. 
    812         resolution = Slit1D(q, width=width, height=height) 
     812        resolution = Slit1D(q, qx_width=width, qy_width=height) 
    813813        answer = romberg_slit_1d(q, width, height, form, pars) 
    814814        output = resolution.apply(eval_form(resolution.q_calc, form, pars)) 
    … …  
    10671067 
    10681068    if isinstance(resolution, Slit1D): 
    1069         width, height = resolution.width, resolution.height 
     1069        width, height = resolution.dqx, resolution.dqy 
    10701070        Iq_romb = romberg_slit_1d(resolution.q, width, height, model, pars) 
    10711071    else: 

  • sasmodels/resolution2d.py

    r823e620 rea75043  
    1010from numpy import pi, cos, sin, sqrt 
    1111 
     12from . import resolution 
    1213from .resolution import Resolution 
    1314 
    … …  
    2021NR = {'xhigh':10, 'high':5, 'med':5, 'low':3} 
    2122NPHI = {'xhigh':20, 'high':12, 'med':6, 'low':4} 
     23 
     24## Defaults 
     25N_SLIT_PERP = {'xhigh':2000, 'high':1000, 'med':500, 'low':250} 
     26N_SLIT_PERP_DOC = ", ".join("%s=%d"%(name,value) for value,name in 
     27                            sorted((v,k) for k,v in N_SLIT_PERP.items())) 
    2228 
    2329class Pinhole2D(Resolution): 
    … …  
    167173        else: 
    168174            return theory 
     175 
     176 
     177class Slit2D(Resolution): 
     178    """ 
     179    Slit aperture with resolution function on an oriented sample. 
     180 
     181    *q* points at which the data is measured. 
     182 
     183    *qx_width* slit width in qx 
     184 
     185    *qy_width* slit height in qy; current implementation requires a fixed 
     186    qy_width for all q points. 
     187 
     188    *q_calc* is the list of q points to calculate, or None if this 
     189    should be estimated from the *q* and *qx_width*. 
     190 
     191    *accuracy* determines the number of *qy* points to compute for each *q*. 
     192    The values are stored in sasmodels.resolution2d.N_SLIT_PERP.  The default 
     193    values are: %s 
     194    """ 
     195    __doc__ = __doc__%N_SLIT_PERP_DOC 
     196    def __init__(self, q, qx_width, qy_width=0., q_calc=None, accuracy='low'): 
     197        # Remember what q and width was used even though we won't need them 
     198        # after the weight matrix is constructed 
     199        self.q, self.qx_width, self.qy_width = q, qx_width, qy_width 
     200 
     201        # Allow independent resolution on each qx point even though it is not 
     202        # needed in practice.  Set qy_width to the maximum qy width. 
     203        if np.isscalar(qx_width): 
     204            qx_width = np.ones(len(q))*qx_width 
     205        else: 
     206            qx_width = np.asarray(qx_width) 
     207        if not np.isscalar(qy_width): 
     208            qy_width = np.max(qy_width) 
     209 
     210        # Build grid of qx, qy points 
     211        if q_calc is not None: 
     212            qx_calc = np.sort(q_calc) 
     213        else: 
     214            qx_calc = resolution.pinhole_extend_q(q, qx_width, nsigma=3) 
     215        qy_calc = np.linspace(-qy_width, qy_width, N_SLIT_PERP[accuracy]*10) 
     216        self.q_calc = [v.flatten() for v in np.meshgrid(qx_calc, qy_calc)] 
     217        self.qx_calc, self.qy_calc = qx_calc, qy_calc 
     218        self.nx, self.ny = len(qx_calc), len(qy_calc) 
     219        self.dy = 2*qy_width/self.ny 
     220 
     221        # Build weight matrix for resolution integration 
     222        if np.any(qx_width > 0): 
     223            self.weights = resolution.pinhole_resolution(qx_calc, q, 
     224                    np.maximum(qx_width, resolution.MINIMUM_RESOLUTION)) 
     225        elif len(qx_calc)==len(q) and np.all(qx_calc == q): 
     226            self.weights = None 
     227        else: 
     228            raise ValueError("Slit2D fails with q_calc != q") 
     229 
     230    def apply(self, theory): 
     231        Iq = np.sum(theory.reshape(self.ny, self.nx), axis=0)*self.dy 
     232        if self.weights is not None: 
     233            Iq = resolution.apply_resolution_matrix(self.weights, Iq) 
     234        return Iq 
     235 
Note: See TracChangeset for help on using the changeset viewer.