Changeset 0656250 in sasmodels for sasmodels/bumps_model.py


Ignore:
Timestamp:
Mar 24, 2015 4:24:43 AM (9 years ago)
Author:
jhbakker
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:
c210c94
Parents:
33d38d5 (diff), 0a33675 (diff)
Note: this is a merge changeset, the changes displayed below correspond to the merge itself.
Use the (diff) links above to see all the changes relative to each parent.
Message:

Add SESANS example data

File:
1 edited

Legend:

Unmodified
Added
Removed

  • sasmodels/bumps_model.py

    r33d38d5 r0656250  
    11""" 
    2 Sasmodels core. 
     2Wrap sasmodels for direct use by bumps. 
    33""" 
    4 import sys, os 
     4 
    55import datetime 
    66 
    7 from sasmodels import sesans 
     7import numpy as np 
     8 
     9from . import sesans 
    810 
    911# CRUFT python 2.6 
    1012if not hasattr(datetime.timedelta, 'total_seconds'): 
    11     def delay(dt): return dt.days*86400 + dt.seconds + 1e-6*dt.microseconds 
     13    def delay(dt): 
     14        """Return number date-time delta as number seconds""" 
     15        return dt.days * 86400 + dt.seconds + 1e-6 * dt.microseconds 
    1216else: 
    13     def delay(dt): return dt.total_seconds() 
    14  
    15 import numpy as np 
    16  
    17 try: 
    18     from .kernelcl import load_model as _loader 
    19 except RuntimeError,exc: 
    20     import warnings 
    21     warnings.warn(str(exc)) 
    22     warnings.warn("OpenCL not available --- using ctypes instead") 
    23     from .kerneldll import load_model as _loader 
    24  
    25 def load_model(modelname, dtype='single'): 
    26     """ 
    27     Load model by name. 
    28     """ 
    29     sasmodels = __import__('sasmodels.models.'+modelname) 
    30     module = getattr(sasmodels.models, modelname, None) 
    31     model = _loader(module, dtype=dtype) 
    32     return model 
     17    def delay(dt): 
     18        """Return number date-time delta as number seconds""" 
     19        return dt.total_seconds() 
    3320 
    3421 
    … …  
    4128    """ 
    4229    then = datetime.datetime.now() 
    43     return lambda: delay(datetime.datetime.now()-then) 
     30    return lambda: delay(datetime.datetime.now() - then) 
    4431 
    4532 
    … …  
    5239    data = loader.load(filename) 
    5340    if data is None: 
    54         raise IOError("Data %r could not be loaded"%filename) 
     41        raise IOError("Data %r could not be loaded" % filename) 
    5542    return data 
    5643 
    … …  
    6552    from sas.dataloader.data_info import Data1D 
    6653 
    67     Iq = 100*np.ones_like(q) 
     54    Iq = 100 * np.ones_like(q) 
    6855    dIq = np.sqrt(Iq) 
    69     data = Data1D(q, Iq, dx=0.05*q, dy=dIq) 
     56    data = Data1D(q, Iq, dx=0.05 * q, dy=dIq) 
    7057    data.filename = "fake data" 
    7158    data.qmin, data.qmax = q.min(), q.max() 
    … …  
    8572    if qy is None: 
    8673        qy = qx 
    87     Qx,Qy = np.meshgrid(qx,qy) 
    88     Qx,Qy = Qx.flatten(), Qy.flatten() 
    89     Iq = 100*np.ones_like(Qx) 
     74    Qx, Qy = np.meshgrid(qx, qy) 
     75    Qx, Qy = Qx.flatten(), Qy.flatten() 
     76    Iq = 100 * np.ones_like(Qx) 
    9077    dIq = np.sqrt(Iq) 
    9178    mask = np.ones(len(Iq), dtype='bool') 
    … …  
    10087 
    10188    # 5% dQ/Q resolution 
    102     data.dqx_data = 0.05*Qx 
    103     data.dqy_data = 0.05*Qy 
     89    data.dqx_data = 0.05 * Qx 
     90    data.dqy_data = 0.05 * Qy 
    10491 
    10592    detector = Detector() 
    … …  
    114101    data.Q_unit = "1/A" 
    115102    data.I_unit = "1/cm" 
    116     data.q_data = np.sqrt(Qx**2 + Qy**2) 
     103    data.q_data = np.sqrt(Qx ** 2 + Qy ** 2) 
    117104    data.xaxis("Q_x", "A^{-1}") 
    118105    data.yaxis("Q_y", "A^{-1}") 
    … …  
    129116        data.mask = Ringcut(0, radius)(data) 
    130117        if outer is not None: 
    131             data.mask += Ringcut(outer,np.inf)(data) 
     118            data.mask += Ringcut(outer, np.inf)(data) 
    132119    else: 
    133         data.mask = (data.x>=radius) 
     120        data.mask = (data.x >= radius) 
    134121        if outer is not None: 
    135             data.mask &= (data.x<outer) 
     122            data.mask &= (data.x < outer) 
    136123 
    137124 
    … …  
    142129    from sas.dataloader.manipulations import Boxcut 
    143130    if half == 'right': 
    144         data.mask += Boxcut(x_min=-np.inf, x_max=0.0, y_min=-np.inf, y_max=np.inf)(data) 
     131        data.mask += \ 
     132            Boxcut(x_min=-np.inf, x_max=0.0, y_min=-np.inf, y_max=np.inf)(data) 
    145133    if half == 'left': 
    146         data.mask += Boxcut(x_min=0.0, x_max=np.inf, y_min=-np.inf, y_max=np.inf)(data) 
    147  
    148  
    149 def set_top(data, max): 
    150     """ 
    151     Chop the top off the data, above *max*. 
     134        data.mask += \ 
     135            Boxcut(x_min=0.0, x_max=np.inf, y_min=-np.inf, y_max=np.inf)(data) 
     136 
     137 
     138def set_top(data, cutoff): 
     139    """ 
     140    Chop the top off the data, above *cutoff*. 
    152141    """ 
    153142    from sas.dataloader.manipulations import Boxcut 
    154     data.mask += Boxcut(x_min=-np.inf, x_max=np.inf, y_min=-np.inf, y_max=max)(data) 
    155  
    156  
    157 def plot_data(data, iq, vmin=None, vmax=None, scale='log'): 
     143    data.mask += \ 
     144        Boxcut(x_min=-np.inf, x_max=np.inf, y_min=-np.inf, y_max=cutoff)(data) 
     145 
     146 
     147def plot_data(data, Iq, vmin=None, vmax=None, view='log'): 
    158148    """ 
    159149    Plot the target value for the data.  This could be the data itself, 
    … …  
    162152    *scale* can be 'log' for log scale data, or 'linear'. 
    163153    """ 
    164     from numpy.ma import masked_array, masked 
     154    from numpy.ma import masked_array 
    165155    import matplotlib.pyplot as plt 
    166156    if hasattr(data, 'qx_data'): 
    167         iq = iq[:] 
    168         valid = np.isfinite(iq) 
    169         if scale == 'log': 
    170             valid[valid] = (iq[valid] > 0) 
    171             iq[valid] = np.log10(iq[valid]) 
    172         iq[~valid|data.mask] = 0 
    173         #plottable = iq 
    174         plottable = masked_array(iq, ~valid|data.mask) 
     157        Iq = Iq + 0 
     158        valid = np.isfinite(Iq) 
     159        if view == 'log': 
     160            valid[valid] = (Iq[valid] > 0) 
     161            Iq[valid] = np.log10(Iq[valid]) 
     162        elif view == 'q4': 
     163            Iq[valid] = Iq*(data.qx_data[valid]**2+data.qy_data[valid]**2)**2 
     164        Iq[~valid | data.mask] = 0 
     165        #plottable = Iq 
     166        plottable = masked_array(Iq, ~valid | data.mask) 
    175167        xmin, xmax = min(data.qx_data), max(data.qx_data) 
    176168        ymin, ymax = min(data.qy_data), max(data.qy_data) 
    177         if vmin is None: vmin = iq[valid&~data.mask].min() 
    178         if vmax is None: vmax = iq[valid&~data.mask].max() 
    179         plt.imshow(plottable.reshape(128,128), 
     169        try: 
     170            if vmin is None: vmin = Iq[valid & ~data.mask].min() 
     171            if vmax is None: vmax = Iq[valid & ~data.mask].max() 
     172        except: 
     173            vmin, vmax = 0, 1 
     174        plt.imshow(plottable.reshape(128, 128), 
    180175                   interpolation='nearest', aspect=1, origin='upper', 
    181176                   extent=[xmin, xmax, ymin, ymax], vmin=vmin, vmax=vmax) 
    182177    else: # 1D data 
    183         if scale == 'linear': 
    184             idx = np.isfinite(iq) 
    185             plt.plot(data.x[idx], iq[idx]) 
    186         else: 
    187             idx = np.isfinite(iq) 
    188             idx[idx] = (iq[idx]>0) 
    189             plt.loglog(data.x[idx], iq[idx]) 
     178        if view == 'linear' or view == 'q4': 
     179            #idx = np.isfinite(Iq) 
     180            scale = data.x**4 if view == 'q4' else 1.0 
     181            plt.plot(data.x, scale*Iq) #, '.') 
     182        else: 
     183            # Find the values that are finite and positive 
     184            idx = np.isfinite(Iq) 
     185            idx[idx] = (Iq[idx] > 0) 
     186            Iq[~idx] = np.nan 
     187            plt.loglog(data.x, Iq) 
    190188 
    191189 
    … …  
    203201        mdata[mdata <= 0] = masked 
    204202    mtheory = masked_array(theory, mdata.mask) 
    205     mresid = masked_array((theory-data.y)/data.dy, mdata.mask) 
    206  
     203    mresid = masked_array((theory - data.y) / data.dy, mdata.mask) 
     204 
     205    scale = data.x**4 if view == 'q4' else 1.0 
    207206    plt.subplot(121) 
    208     plt.errorbar(data.x, mdata, yerr=data.dy) 
    209     plt.plot(data.x, mtheory, '-', hold=True) 
    210     plt.yscale(view) 
     207    plt.errorbar(data.x, scale*mdata, yerr=data.dy) 
     208    plt.plot(data.x, scale*mtheory, '-', hold=True) 
     209    plt.yscale('linear' if view == 'q4' else view) 
    211210    plt.subplot(122) 
    212211    plt.plot(data.x, mresid, 'x') 
    213     #plt.axhline(1, color='black', ls='--',lw=1, hold=True) 
    214     #plt.axhline(0, color='black', lw=1, hold=True) 
    215     #plt.axhline(-1, color='black', ls='--',lw=1, hold=True) 
    216  
     212 
     213# pylint: disable=unused-argument 
    217214def _plot_sesans(data, theory, view): 
    218215    import matplotlib.pyplot as plt 
    219     resid = (theory - data.y)/data.dy 
     216    resid = (theory - data.y) / data.dy 
    220217    plt.subplot(121) 
    221218    plt.errorbar(data.x, data.y, yerr=data.dy) 
    … …  
    233230    """ 
    234231    import matplotlib.pyplot as plt 
    235     resid = (theory-data.data)/data.err_data 
     232    resid = (theory - data.data) / data.err_data 
    236233    plt.subplot(131) 
    237     plot_data(data, data.data, scale=view) 
     234    plot_data(data, data.data, view=view) 
    238235    plt.colorbar() 
    239236    plt.subplot(132) 
    240     plot_data(data, theory, scale=view) 
     237    plot_data(data, theory, view=view) 
    241238    plt.colorbar() 
    242239    plt.subplot(133) 
    243     plot_data(data, resid, scale='linear') 
     240    plot_data(data, resid, view='linear') 
    244241    plt.colorbar() 
    245242 
    … …  
    250247    *data* is the data to be fitted. 
    251248 
    252     *model* is the SAS model, e.g., from :func:`gen.opencl_model`. 
     249    *model* is the SAS model from :func:`core.load_model`. 
    253250 
    254251    *cutoff* is the integration cutoff, which avoids computing the 
    … …  
    267264        self.model = model 
    268265        self.cutoff = cutoff 
    269 # TODO       if  isinstance(data,SESANSData1D)         
    270266        if hasattr(data, 'lam'): 
    271267            self.data_type = 'sesans' 
    … …  
    280276        if self.data_type == 'sesans': 
    281277            q = sesans.make_q(data.sample.zacceptance, data.Rmax) 
    282             self.index = slice(None,None) 
    283             self.iq = data.y 
    284             self.diq = data.dy 
     278            self.index = slice(None, None) 
     279            self.Iq = data.y 
     280            self.dIq = data.dy 
    285281            self._theory = np.zeros_like(q) 
    286282            q_vectors = [q] 
    287283        elif self.data_type == 'Iqxy': 
    288             self.index = (data.mask==0) & (~np.isnan(data.data)) 
    289             self.iq = data.data[self.index] 
    290             self.diq = data.err_data[self.index] 
     284            self.index = (data.mask == 0) & (~np.isnan(data.data)) 
     285            self.Iq = data.data[self.index] 
     286            self.dIq = data.err_data[self.index] 
    291287            self._theory = np.zeros_like(data.data) 
    292288            if not partype['orientation'] and not partype['magnetic']: 
    293                 q_vectors = [np.sqrt(data.qx_data**2+data.qy_data**2)] 
     289                q_vectors = [np.sqrt(data.qx_data ** 2 + data.qy_data ** 2)] 
    294290            else: 
    295291                q_vectors = [data.qx_data, data.qy_data] 
    296292        elif self.data_type == 'Iq': 
    297             self.index = (data.x>=data.qmin) & (data.x<=data.qmax) & ~np.isnan(data.y) 
    298             self.iq = data.y[self.index] 
    299             self.diq = data.dy[self.index] 
     293            self.index = (data.x >= data.qmin) & (data.x <= data.qmax) & ~np.isnan(data.y) 
     294            self.Iq = data.y[self.index] 
     295            self.dIq = data.dy[self.index] 
    300296            self._theory = np.zeros_like(data.y) 
    301297            q_vectors = [data.x] 
    … …  
    311307        pars = [] 
    312308        for p in model.info['parameters']: 
    313             name, default, limits, ptype = p[0], p[2], p[3], p[4] 
     309            name, default, limits = p[0], p[2], p[3] 
    314310            value = kw.pop(name, default) 
    315311            setattr(self, name, Parameter.default(value, name=name, limits=limits)) 
    316312            pars.append(name) 
    317313        for name in partype['pd-2d']: 
    318             for xpart,xdefault,xlimits in [ 
     314            for xpart, xdefault, xlimits in [ 
    319315                    ('_pd', 0, limits), 
    320                     ('_pd_n', 35, (0,1000)), 
     316                    ('_pd_n', 35, (0, 1000)), 
    321317                    ('_pd_nsigma', 3, (0, 10)), 
    322318                    ('_pd_type', 'gaussian', None), 
    323319                ]: 
    324                 xname = name+xpart 
     320                xname = name + xpart 
    325321                xvalue = kw.pop(xname, xdefault) 
    326322                if xlimits is not None: 
    … …  
    330326        self._parameter_names = pars 
    331327        if kw: 
    332             raise TypeError("unexpected parameters: %s"%(", ".join(sorted(kw.keys())))) 
     328            raise TypeError("unexpected parameters: %s" 
     329                            % (", ".join(sorted(kw.keys())))) 
    333330        self.update() 
    334331 
    … …  
    337334 
    338335    def numpoints(self): 
    339         return len(self.iq) 
     336        """ 
     337            Return the number of points 
     338        """ 
     339        return len(self.Iq) 
    340340 
    341341    def parameters(self): 
    342         return dict((k,getattr(self,k)) for k in self._parameter_names) 
     342        """ 
     343            Return a dictionary of parameters 
     344        """ 
     345        return dict((k, getattr(self, k)) for k in self._parameter_names) 
    343346 
    344347    def theory(self): 
    345348        if 'theory' not in self._cache: 
    346349            if self._fn is None: 
    347                 input = self.model.make_input(self._fn_inputs) 
    348                 self._fn = self.model(input) 
    349  
    350             pars = [getattr(self,p).value for p in self._fn.fixed_pars] 
     350                q_input = self.model.make_input(self._fn_inputs) 
     351                self._fn = self.model(q_input) 
     352 
     353            fixed_pars = [getattr(self, p).value for p in self._fn.fixed_pars] 
    351354            pd_pars = [self._get_weights(p) for p in self._fn.pd_pars] 
    352             #print pars 
    353             self._theory[self.index] = self._fn(pars, pd_pars, self.cutoff) 
     355            #print fixed_pars,pd_pars 
     356            self._theory[self.index] = self._fn(fixed_pars, pd_pars, self.cutoff) 
    354357            #self._theory[:] = self._fn.eval(pars, pd_pars) 
    355358            if self.data_type == 'sesans': 
    356                 P = sesans.hankel(self.data.x, self.data.lam*1e-9, 
    357                                   self.data.sample.thickness/10, self._fn_inputs[0], 
    358                                   self._theory) 
    359                 self._cache['theory'] = P 
     359                result = sesans.hankel(self.data.x, self.data.lam * 1e-9, 
     360                                       self.data.sample.thickness / 10, 
     361                                       self._fn_inputs[0], self._theory) 
     362                self._cache['theory'] = result 
    360363            else: 
    361364                self._cache['theory'] = self._theory 
    … …  
    364367    def residuals(self): 
    365368        #if np.any(self.err ==0): print "zeros in err" 
    366         return (self.theory()[self.index]-self.iq)/self.diq 
     369        return (self.theory()[self.index] - self.Iq) / self.dIq 
    367370 
    368371    def nllf(self): 
    369         R = self.residuals() 
     372        delta = self.residuals() 
    370373        #if np.any(np.isnan(R)): print "NaN in residuals" 
    371         return 0.5*np.sum(R**2) 
    372  
    373     def __call__(self): 
    374         return 2*self.nllf()/self.dof 
     374        return 0.5 * np.sum(delta ** 2) 
     375 
     376    #def __call__(self): 
     377    #    return 2 * self.nllf() / self.dof 
    375378 
    376379    def plot(self, view='log'): 
    … …  
    391394        print "noise", noise 
    392395        if noise is None: 
    393             noise = self.diq[self.index] 
    394         else: 
    395             noise = 0.01*noise 
    396             self.diq[self.index] = noise 
     396            noise = self.dIq[self.index] 
     397        else: 
     398            noise = 0.01 * noise 
     399            self.dIq[self.index] = noise 
    397400        y = self.theory() 
    398         y += y*np.random.randn(*y.shape)*noise 
     401        y += y * np.random.randn(*y.shape) * noise 
    399402        if self.data_type == 'Iq': 
    400403            self.data.y[self.index] = y 
    … …  
    410413 
    411414    def _get_weights(self, par): 
     415        """ 
     416        Get parameter dispersion weights 
     417        """ 
    412418        from . import weights 
    413419 
    414420        relative = self.model.info['partype']['pd-rel'] 
    415421        limits = self.model.info['limits'] 
    416         disperser,value,npts,width,nsigma = [getattr(self, par+ext) 
    417                 for ext in ('_pd_type','','_pd_n','_pd','_pd_nsigma')] 
    418         v,w = weights.get_weights( 
     422        disperser, value, npts, width, nsigma = [ 
     423            getattr(self, par + ext) 
     424            for ext in ('_pd_type', '', '_pd_n', '_pd', '_pd_nsigma')] 
     425        value, weight = weights.get_weights( 
    419426            disperser, int(npts.value), width.value, nsigma.value, 
    420427            value.value, limits[par], par in relative) 
    421         return v,w/w.max() 
     428        return value, weight / np.sum(weight) 
    422429 
    423430    def __getstate__(self): 
    … …  
    428435 
    429436    def __setstate__(self, state): 
     437        # pylint: disable=attribute-defined-outside-init 
    430438        self.__dict__ = state 
    431439 
    … …  
    434442    data = load_data('DEC07086.DAT') 
    435443    set_beam_stop(data, 0.004) 
    436     plot_data(data) 
     444    plot_data(data, data.data) 
    437445    import matplotlib.pyplot as plt; plt.show() 
    438446 
Note: See TracChangeset for help on using the changeset viewer.