Changes in / [d60b433:c9e31e2] in sasmodels


Ignore:
Files:
23 added
20 deleted
43 edited

Legend:

Unmodified
Added
Removed

  • .gitignore

    r946cdc8e rd6adfbe  
    99/doc/api/ 
    1010/doc/model/ 
     11/doc/ref/models 
    1112.mplconfig 
     13/pylint_violations.txt 
     14/coverage.xml 
     15/.coverage 

  • compare.py

    r29f5536 rb89f519  
    8181        return np.random.rand() 
    8282    else: 
    83         # length, scale, background in [0,200] 
    84         return 200*np.random.rand() 
     83        # values from 0 to 2*x for all other parameters 
     84        return 2*np.random.rand()*(v if v != 0 else 1) 
    8585 
    8686def randomize_model(name, pars, seed=None): 
    … …  
    145145    return value, average_time 
    146146 
    147 def make_data(qmax, is2D, Nq=128): 
     147def make_data(qmax, is2D, Nq=128, view='log'): 
    148148    if is2D: 
    149149        from sasmodels.bumps_model import empty_data2D, set_beam_stop 
    … …  
    153153    else: 
    154154        from sasmodels.bumps_model import empty_data1D 
    155         qmax = math.log10(qmax) 
    156         data = empty_data1D(np.logspace(qmax-3, qmax, Nq)) 
     155        if view == 'log': 
     156            qmax = math.log10(qmax) 
     157            q = np.logspace(qmax-3, qmax, Nq) 
     158        else: 
     159            q = np.linspace(0.001*qmax, qmax, Nq) 
     160        data = empty_data1D(q) 
    157161        index = slice(None, None) 
    158162    return data, index 
    159163 
    160164def compare(name, pars, Ncpu, Nocl, opts, set_pars): 
     165    view = 'linear' if '-linear' in opts else 'log' if '-log' in opts else 'q4' if '-q4' in opts else 'log' 
     166 
    161167    opt_values = dict(split 
    162168                      for s in opts for split in ((s.split('='),)) 
    … …  
    166172    Nq = int(opt_values.get('-Nq', '128')) 
    167173    is2D = not "-1d" in opts 
    168     data, index = make_data(qmax, is2D, Nq) 
     174    data, index = make_data(qmax, is2D, Nq, view=view) 
    169175 
    170176 
    … …  
    174180 
    175181    # randomize parameters 
     182    pars.update(set_pars) 
    176183    if '-random' in opts or '-random' in opt_values: 
    177184        seed = int(opt_values['-random']) if '-random' in opt_values else None 
    178185        pars, seed = randomize_model(name, pars, seed=seed) 
    179186        print "Randomize using -random=%i"%seed 
    180     pars.update(set_pars) 
    181187 
    182188    # parameter selection 
    … …  
    223229    if Ncpu > 0: 
    224230        if Nocl > 0: plt.subplot(131) 
    225         plot_data(data, cpu, scale='log') 
     231        plot_data(data, cpu, view=view) 
    226232        plt.title("%s t=%.1f ms"%(comp,cpu_time)) 
    227233        cbar_title = "log I" 
    228234    if Nocl > 0: 
    229235        if Ncpu > 0: plt.subplot(132) 
    230         plot_data(data, ocl, scale='log') 
     236        plot_data(data, ocl, view=view) 
    231237        plt.title("opencl t=%.1f ms"%ocl_time) 
    232238        cbar_title = "log I" 
    … …  
    234240        plt.subplot(133) 
    235241        if '-abs' in opts: 
    236             err,errstr = resid, "abs err" 
    237         else: 
    238             err,errstr = relerr, "rel err" 
     242            err,errstr,errview = resid, "abs err", "linear" 
     243        else: 
     244            err,errstr,errview = abs(relerr), "rel err", "log" 
    239245        #err,errstr = ocl/cpu,"ratio" 
    240         plot_data(data, err, scale='log') #'linear') 
     246        plot_data(data, err, view=errview) 
    241247        plt.title("max %s = %.3g"%(errstr, max(abs(err[index])))) 
    242         cbar_title = "log "+errstr 
     248        cbar_title = errstr if errview=="linear" else "log "+errstr 
    243249    if is2D: 
    244250        h = plt.colorbar() 
    … …  
    281287    -pars/-nopars* prints the parameter set or not 
    282288    -abs/-rel* plot relative or absolute error 
     289    -linear/-log/-q4 intensity scaling 
    283290    -hist/-nohist* plot histogram of relative error 
    284291 
    … …  
    301308    'nopars','pars', 
    302309    'rel','abs', 
     310    'linear', 'log', 'q4', 
    303311    'hist','nohist', 
    304312    ]) 

  • compare.sh

    • Property mode changed from 100644 to 100755
    r8a3e0af r5134b2c  
    55#PYOPENCL_CTX=${CTX:-1} 
    66PYTHONPATH=../bumps:../periodictable:$SASVIEW 
    7 export PYOPENCL_CTX PYTHONPATH 
     7export PYTHONPATH 
    88 
    9 echo PYTHONPATH=$PYTHONPATH 
    109set -x 
    1110 

  • compare_all.sh

    r34d49af r5134b2c  
    22 
    33SASVIEW=$(ls -d ../sasview/build/lib.*) 
    4 PYOPENCL_CTX=1 
    54PYTHONPATH=../bumps:../periodictable:$SASVIEW 
    6 export PYOPENCL_CTX PYTHONPATH 
     5export PYTHONPATH 
    76 
    87echo PYTHONPATH=$PYTHONPATH 

  • doc/Makefile

    r19dcb933 r61ba623  
    3232MODELS_RST = $(patsubst ../sasmodels/models/%.py,model/%.rst,$(MODELS_PY)) 
    3333 
    34 model/img/%: ../sasmodels/models/img/% 
    35         $(COPY) $< $@ 
    36  
    37 model/%.rst: ../sasmodels/models/%.py 
    38         $(PYTHON) genmodel.py $< $@ 
    39  
    40 .PHONY: help clean html dirhtml pickle json htmlhelp qthelp latex changes linkcheck doctest build 
    41  
    4234help: 
    4335        @echo "Please use \`make <target>' where <target> is one of" 
    … …  
    5345        @echo "  doctest   to run all doctests embedded in the documentation (if enabled)" 
    5446 
     47model/img/%: ../sasmodels/models/img/% 
     48        $(COPY) $< $@ 
     49 
     50model/%.rst: ../sasmodels/models/%.py 
     51        $(PYTHON) genmodel.py $< $@ 
     52 
     53ref/models/index.rst: gentoc.py $(MODELS_PY) 
     54        $(PYTHON) gentoc.py $(MODELS_PY) 
     55 
     56.PHONY: help clean html dirhtml pickle json htmlhelp qthelp latex changes linkcheck doctest build 
     57 
    5558api: genapi.py 
    5659        $(RMDIR) api 
    … …  
    6164        $(MKDIR) model/img 
    6265 
    63 build: model $(MODELS_RST) $(IMG_DEST) api 
     66build: model $(MODELS_RST) $(IMG_DEST) api ref/models/index.rst 
    6467        #cd ../.. && python setup.py build 
    6568 

  • doc/genmodel.py

    rcb6ecf4 r61ba623  
    33sys.path.insert(0,'..') 
    44 
    5 # Convert ../sasmodels/models.name.py to sasmodels.models.name 
     5# Convert ../sasmodels/models/name.py to sasmodels.models.name 
    66module_name = sys.argv[1][3:-3].replace('/','.').replace('\\','.') 
    77print module_name 

  • doc/index.rst

    r19dcb933 r61ba623  
    1515 
    1616   guide/index.rst 
    17    models/index.rst 
     17   ref/index.rst 
    1818   api/index.rst 
    1919 

  • sasmodels/alignment.py

    rff7119b r3c56da87  
    3434    return view 
    3535 
    36 def align_data(v, dtype, alignment=128): 
     36def align_data(x, dtype, alignment=128): 
    3737    """ 
    3838    Return a copy of an array on the alignment boundary. 
    3939    """ 
    40     # if v is contiguous, aligned, and of the correct type then just return v 
    41     view = align_empty(v.shape, dtype, alignment=alignment) 
    42     view[:] = v 
     40    # if x is contiguous, aligned, and of the correct type then just return x 
     41    view = align_empty(x.shape, dtype, alignment=alignment) 
     42    view[:] = x 
    4343    return view 
    4444 

  • sasmodels/bumps_model.py

    r29f5536 r3c56da87  
    22Sasmodels core. 
    33""" 
    4 import sys, os 
    54import datetime 
    65 
    … …  
    98# CRUFT python 2.6 
    109if not hasattr(datetime.timedelta, 'total_seconds'): 
    11     def delay(dt): return dt.days*86400 + dt.seconds + 1e-6*dt.microseconds 
     10    def delay(dt): 
     11        """Return number date-time delta as number seconds""" 
     12        return dt.days * 86400 + dt.seconds + 1e-6 * dt.microseconds 
    1213else: 
    13     def delay(dt): return dt.total_seconds() 
     14    def delay(dt): 
     15        """Return number date-time delta as number seconds""" 
     16        return dt.total_seconds() 
    1417 
    1518import numpy as np 
    … …  
    1720try: 
    1821    from .kernelcl import load_model as _loader 
    19 except RuntimeError,exc: 
     22except RuntimeError, exc: 
    2023    import warnings 
    2124    warnings.warn(str(exc)) 
    … …  
    2730    Load model by name. 
    2831    """ 
    29     sasmodels = __import__('sasmodels.models.'+modelname) 
     32    sasmodels = __import__('sasmodels.models.' + modelname) 
    3033    module = getattr(sasmodels.models, modelname, None) 
    3134    model = _loader(module, dtype=dtype) 
    … …  
    4144    """ 
    4245    then = datetime.datetime.now() 
    43     return lambda: delay(datetime.datetime.now()-then) 
     46    return lambda: delay(datetime.datetime.now() - then) 
    4447 
    4548 
    … …  
    5255    data = loader.load(filename) 
    5356    if data is None: 
    54         raise IOError("Data %r could not be loaded"%filename) 
     57        raise IOError("Data %r could not be loaded" % filename) 
    5558    return data 
    5659 
    … …  
    6568    from sas.dataloader.data_info import Data1D 
    6669 
    67     Iq = 100*np.ones_like(q) 
     70    Iq = 100 * np.ones_like(q) 
    6871    dIq = np.sqrt(Iq) 
    69     data = Data1D(q, Iq, dx=0.05*q, dy=dIq) 
     72    data = Data1D(q, Iq, dx=0.05 * q, dy=dIq) 
    7073    data.filename = "fake data" 
    7174    data.qmin, data.qmax = q.min(), q.max() 
    … …  
    8588    if qy is None: 
    8689        qy = qx 
    87     Qx,Qy = np.meshgrid(qx,qy) 
    88     Qx,Qy = Qx.flatten(), Qy.flatten() 
    89     Iq = 100*np.ones_like(Qx) 
     90    Qx, Qy = np.meshgrid(qx, qy) 
     91    Qx, Qy = Qx.flatten(), Qy.flatten() 
     92    Iq = 100 * np.ones_like(Qx) 
    9093    dIq = np.sqrt(Iq) 
    9194    mask = np.ones(len(Iq), dtype='bool') 
    … …  
    100103 
    101104    # 5% dQ/Q resolution 
    102     data.dqx_data = 0.05*Qx 
    103     data.dqy_data = 0.05*Qy 
     105    data.dqx_data = 0.05 * Qx 
     106    data.dqy_data = 0.05 * Qy 
    104107 
    105108    detector = Detector() 
    … …  
    114117    data.Q_unit = "1/A" 
    115118    data.I_unit = "1/cm" 
    116     data.q_data = np.sqrt(Qx**2 + Qy**2) 
     119    data.q_data = np.sqrt(Qx ** 2 + Qy ** 2) 
    117120    data.xaxis("Q_x", "A^{-1}") 
    118121    data.yaxis("Q_y", "A^{-1}") 
    … …  
    129132        data.mask = Ringcut(0, radius)(data) 
    130133        if outer is not None: 
    131             data.mask += Ringcut(outer,np.inf)(data) 
     134            data.mask += Ringcut(outer, np.inf)(data) 
    132135    else: 
    133         data.mask = (data.x>=radius) 
     136        data.mask = (data.x >= radius) 
    134137        if outer is not None: 
    135             data.mask &= (data.x<outer) 
     138            data.mask &= (data.x < outer) 
    136139 
    137140 
    … …  
    142145    from sas.dataloader.manipulations import Boxcut 
    143146    if half == 'right': 
    144         data.mask += Boxcut(x_min=-np.inf, x_max=0.0, y_min=-np.inf, y_max=np.inf)(data) 
     147        data.mask += \ 
     148            Boxcut(x_min=-np.inf, x_max=0.0, y_min=-np.inf, y_max=np.inf)(data) 
    145149    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*. 
     150        data.mask += \ 
     151            Boxcut(x_min=0.0, x_max=np.inf, y_min=-np.inf, y_max=np.inf)(data) 
     152 
     153 
     154def set_top(data, cutoff): 
     155    """ 
     156    Chop the top off the data, above *cutoff*. 
    152157    """ 
    153158    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'): 
     159    data.mask += \ 
     160        Boxcut(x_min=-np.inf, x_max=np.inf, y_min=-np.inf, y_max=cutoff)(data) 
     161 
     162 
     163def plot_data(data, Iq, vmin=None, vmax=None, view='log'): 
    158164    """ 
    159165    Plot the target value for the data.  This could be the data itself, 
    … …  
    162168    *scale* can be 'log' for log scale data, or 'linear'. 
    163169    """ 
    164     from numpy.ma import masked_array, masked 
     170    from numpy.ma import masked_array 
    165171    import matplotlib.pyplot as plt 
    166172    if hasattr(data, 'qx_data'): 
    167         iq = iq+0 
    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) 
     173        Iq = Iq + 0 
     174        valid = np.isfinite(Iq) 
     175        if view == 'log': 
     176            valid[valid] = (Iq[valid] > 0) 
     177            Iq[valid] = np.log10(Iq[valid]) 
     178        elif view == 'q4': 
     179            Iq[valid] = Iq*(data.qx_data[valid]**2+data.qy_data[valid]**2)**2 
     180        Iq[~valid | data.mask] = 0 
     181        #plottable = Iq 
     182        plottable = masked_array(Iq, ~valid | data.mask) 
    175183        xmin, xmax = min(data.qx_data), max(data.qx_data) 
    176184        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), 
     185        try: 
     186            if vmin is None: vmin = Iq[valid & ~data.mask].min() 
     187            if vmax is None: vmax = Iq[valid & ~data.mask].max() 
     188        except: 
     189            vmin, vmax = 0, 1 
     190        plt.imshow(plottable.reshape(128, 128), 
    180191                   interpolation='nearest', aspect=1, origin='upper', 
    181192                   extent=[xmin, xmax, ymin, ymax], vmin=vmin, vmax=vmax) 
    182193    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]) 
     194        if view == 'linear' or view == 'q4': 
     195            #idx = np.isfinite(Iq) 
     196            scale = data.x**4 if view == 'q4' else 1.0 
     197            plt.plot(data.x, scale*Iq) #, '.') 
     198        else: 
     199            # Find the values that are finite and positive 
     200            idx = np.isfinite(Iq) 
     201            idx[idx] = (Iq[idx] > 0) 
     202            Iq[~idx] = np.nan 
     203            plt.loglog(data.x, Iq) 
    190204 
    191205 
    … …  
    203217        mdata[mdata <= 0] = masked 
    204218    mtheory = masked_array(theory, mdata.mask) 
    205     mresid = masked_array((theory-data.y)/data.dy, mdata.mask) 
    206  
     219    mresid = masked_array((theory - data.y) / data.dy, mdata.mask) 
     220 
     221    scale = data.x**4 if view == 'q4' else 1.0 
    207222    plt.subplot(121) 
    208     plt.errorbar(data.x, mdata, yerr=data.dy) 
    209     plt.plot(data.x, mtheory, '-', hold=True) 
    210     plt.yscale(view) 
     223    plt.errorbar(data.x, scale*mdata, yerr=data.dy) 
     224    plt.plot(data.x, scale*mtheory, '-', hold=True) 
     225    plt.yscale('linear' if view == 'q4' else view) 
    211226    plt.subplot(122) 
    212227    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  
     228 
     229# pylint: disable=unused-argument 
    217230def _plot_sesans(data, theory, view): 
    218231    import matplotlib.pyplot as plt 
    219     resid = (theory - data.y)/data.dy 
     232    resid = (theory - data.y) / data.dy 
    220233    plt.subplot(121) 
    221234    plt.errorbar(data.x, data.y, yerr=data.dy) 
    … …  
    233246    """ 
    234247    import matplotlib.pyplot as plt 
    235     resid = (theory-data.data)/data.err_data 
     248    resid = (theory - data.data) / data.err_data 
    236249    plt.subplot(131) 
    237     plot_data(data, data.data, scale=view) 
     250    plot_data(data, data.data, view=view) 
    238251    plt.colorbar() 
    239252    plt.subplot(132) 
    240     plot_data(data, theory, scale=view) 
     253    plot_data(data, theory, view=view) 
    241254    plt.colorbar() 
    242255    plt.subplot(133) 
    243     plot_data(data, resid, scale='linear') 
     256    plot_data(data, resid, view='linear') 
    244257    plt.colorbar() 
    245258 
    … …  
    267280        self.model = model 
    268281        self.cutoff = cutoff 
    269 # TODO       if  isinstance(data,SESANSData1D)         
     282# TODO       if  isinstance(data,SESANSData1D) 
    270283        if hasattr(data, 'lam'): 
    271284            self.data_type = 'sesans' 
    … …  
    280293        if self.data_type == 'sesans': 
    281294            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 
     295            self.index = slice(None, None) 
     296            self.Iq = data.y 
     297            self.dIq = data.dy 
    285298            self._theory = np.zeros_like(q) 
    286299            q_vectors = [q] 
    287300        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] 
     301            self.index = (data.mask == 0) & (~np.isnan(data.data)) 
     302            self.Iq = data.data[self.index] 
     303            self.dIq = data.err_data[self.index] 
    291304            self._theory = np.zeros_like(data.data) 
    292305            if not partype['orientation'] and not partype['magnetic']: 
    293                 q_vectors = [np.sqrt(data.qx_data**2+data.qy_data**2)] 
     306                q_vectors = [np.sqrt(data.qx_data ** 2 + data.qy_data ** 2)] 
    294307            else: 
    295308                q_vectors = [data.qx_data, data.qy_data] 
    296309        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] 
     310            self.index = (data.x >= data.qmin) & (data.x <= data.qmax) & ~np.isnan(data.y) 
     311            self.Iq = data.y[self.index] 
     312            self.dIq = data.dy[self.index] 
    300313            self._theory = np.zeros_like(data.y) 
    301314            q_vectors = [data.x] 
    … …  
    311324        pars = [] 
    312325        for p in model.info['parameters']: 
    313             name, default, limits, ptype = p[0], p[2], p[3], p[4] 
     326            name, default, limits = p[0], p[2], p[3] 
    314327            value = kw.pop(name, default) 
    315328            setattr(self, name, Parameter.default(value, name=name, limits=limits)) 
    316329            pars.append(name) 
    317330        for name in partype['pd-2d']: 
    318             for xpart,xdefault,xlimits in [ 
     331            for xpart, xdefault, xlimits in [ 
    319332                    ('_pd', 0, limits), 
    320                     ('_pd_n', 35, (0,1000)), 
     333                    ('_pd_n', 35, (0, 1000)), 
    321334                    ('_pd_nsigma', 3, (0, 10)), 
    322335                    ('_pd_type', 'gaussian', None), 
    323336                ]: 
    324                 xname = name+xpart 
     337                xname = name + xpart 
    325338                xvalue = kw.pop(xname, xdefault) 
    326339                if xlimits is not None: 
    … …  
    330343        self._parameter_names = pars 
    331344        if kw: 
    332             raise TypeError("unexpected parameters: %s"%(", ".join(sorted(kw.keys())))) 
     345            raise TypeError("unexpected parameters: %s" 
     346                            % (", ".join(sorted(kw.keys())))) 
    333347        self.update() 
    334348 
    … …  
    337351 
    338352    def numpoints(self): 
    339         return len(self.iq) 
     353        """ 
     354            Return the number of points 
     355        """ 
     356        return len(self.Iq) 
    340357 
    341358    def parameters(self): 
    342         return dict((k,getattr(self,k)) for k in self._parameter_names) 
     359        """ 
     360            Return a dictionary of parameters 
     361        """ 
     362        return dict((k, getattr(self, k)) for k in self._parameter_names) 
    343363 
    344364    def theory(self): 
    345365        if 'theory' not in self._cache: 
    346366            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] 
     367                input_value = self.model.make_input(self._fn_inputs) 
     368                self._fn = self.model(input_value) 
     369 
     370            fixed_pars = [getattr(self, p).value for p in self._fn.fixed_pars] 
    351371            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) 
     372            #print fixed_pars,pd_pars 
     373            self._theory[self.index] = self._fn(fixed_pars, pd_pars, self.cutoff) 
    354374            #self._theory[:] = self._fn.eval(pars, pd_pars) 
    355375            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 
     376                result = sesans.hankel(self.data.x, self.data.lam * 1e-9, 
     377                                       self.data.sample.thickness / 10, 
     378                                       self._fn_inputs[0], self._theory) 
     379                self._cache['theory'] = result 
    360380            else: 
    361381                self._cache['theory'] = self._theory 
    … …  
    364384    def residuals(self): 
    365385        #if np.any(self.err ==0): print "zeros in err" 
    366         return (self.theory()[self.index]-self.iq)/self.diq 
     386        return (self.theory()[self.index] - self.Iq) / self.dIq 
    367387 
    368388    def nllf(self): 
    369         R = self.residuals() 
     389        delta = self.residuals() 
    370390        #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 
     391        return 0.5 * np.sum(delta ** 2) 
     392 
     393    #def __call__(self): 
     394    #    return 2 * self.nllf() / self.dof 
    375395 
    376396    def plot(self, view='log'): 
    … …  
    391411        print "noise", noise 
    392412        if noise is None: 
    393             noise = self.diq[self.index] 
    394         else: 
    395             noise = 0.01*noise 
    396             self.diq[self.index] = noise 
     413            noise = self.dIq[self.index] 
     414        else: 
     415            noise = 0.01 * noise 
     416            self.dIq[self.index] = noise 
    397417        y = self.theory() 
    398         y += y*np.random.randn(*y.shape)*noise 
     418        y += y * np.random.randn(*y.shape) * noise 
    399419        if self.data_type == 'Iq': 
    400420            self.data.y[self.index] = y 
    … …  
    410430 
    411431    def _get_weights(self, par): 
     432        """ 
     433            Get parameter dispersion weights 
     434        """ 
    412435        from . import weights 
    413436 
    414437        relative = self.model.info['partype']['pd-rel'] 
    415438        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( 
     439        disperser, value, npts, width, nsigma = [ 
     440            getattr(self, par + ext) 
     441            for ext in ('_pd_type', '', '_pd_n', '_pd', '_pd_nsigma')] 
     442        value, weight = weights.get_weights( 
    419443            disperser, int(npts.value), width.value, nsigma.value, 
    420444            value.value, limits[par], par in relative) 
    421         return v,w/w.max() 
     445        return value, weight / np.sum(weight) 
    422446 
    423447    def __getstate__(self): 
    … …  
    428452 
    429453    def __setstate__(self, state): 
     454        # pylint: disable=attribute-defined-outside-init 
    430455        self.__dict__ = state 
    431456 
    … …  
    434459    data = load_data('DEC07086.DAT') 
    435460    set_beam_stop(data, 0.004) 
    436     plot_data(data) 
     461    plot_data(data, data.data) 
    437462    import matplotlib.pyplot as plt; plt.show() 
    438463 

  • sasmodels/convert.py

    rfd1c792 r3c56da87  
    3737    Convert model from old style parameter names to new style. 
    3838    """ 
     39    _,_ = name,pars # lint 
    3940    raise NotImplementedError 
    4041    # need to load all new models in order to determine old=>new 
    … …  
    5455    Convert model from new style parameter names to old style. 
    5556    """ 
    56     sasmodels = __import__('sasmodels.models.'+name) 
    57     newmodel = getattr(sasmodels.models, name, None) 
    58     mapping = newmodel.oldpars 
    59     oldname = newmodel.oldname 
     57    __import__('sasmodels.models.'+name) 
     58    from . import models 
     59    model = getattr(models, name, None) 
     60    mapping = model.oldpars 
     61    oldname = model.oldname 
    6062    oldpars = _rename_pars(_unscale_sld(pars), mapping) 
    6163    return oldname, oldpars 

  • sasmodels/direct_model.py

    r16bc3fc rf734e7d  
    33import numpy as np 
    44 
    5 from . import models 
    6 from . import weights 
    7  
    8 try: 
    9     from .kernelcl import load_model 
    10 except ImportError,exc: 
    11     warnings.warn(str(exc)) 
    12     warnings.warn("using ctypes instead") 
    13     from .kerneldll import load_model 
    14  
    15 def load_model_definition(model_name): 
    16     __import__('sasmodels.models.'+model_name) 
    17     model_definition = getattr(models, model_name, None) 
    18     return model_definition 
    19  
    20 # load_model is imported above.  It looks like the following 
    21 #def load_model(model_definition, dtype='single): 
    22 #    if kerneldll: 
    23 #        if source is newer than compiled: compile 
    24 #        load dll 
    25 #        return kernel 
    26 #    elif kernelcl: 
    27 #        compile source on context 
    28 #        return kernel 
    29  
    30  
    31 def make_kernel(model, q_vectors): 
    32     """ 
    33     Return a computation kernel from the model definition and the q input. 
    34     """ 
    35     input = model.make_input(q_vectors) 
    36     return model(input) 
    37  
    38 def get_weights(kernel, pars, name): 
    39     """ 
    40     Generate the distribution for parameter *name* given the parameter values 
    41     in *pars*. 
    42  
    43     Searches for "name", "name_pd", "name_pd_type", "name_pd_n", "name_pd_sigma" 
    44     """ 
    45     relative = name in kernel.info['partype']['pd-rel'] 
    46     limits = kernel.info['limits'] 
    47     disperser = pars.get(name+'_pd_type', 'gaussian') 
    48     value = pars.get(name, kernel.info['defaults'][name]) 
    49     npts = pars.get(name+'_pd_n', 0) 
    50     width = pars.get(name+'_pd', 0.0) 
    51     nsigma = pars.get(name+'_pd_nsigma', 3.0) 
    52     v,w = weights.get_weights( 
    53         disperser, npts, width, nsigma, 
    54         value, limits[name], relative) 
    55     return v,w/np.sum(w) 
    56  
    57 def call_kernel(kernel, pars): 
    58     fixed_pars = [pars.get(name, kernel.info['defaults'][name]) 
    59                   for name in kernel.fixed_pars] 
    60     pd_pars = [get_weights(kernel, pars, name) for name in kernel.pd_pars] 
    61     return kernel(fixed_pars, pd_pars) 
     5from .core import load_model_definition, make_kernel, call_kernel 
     6from .core import load_model_cl as load_model 
     7if load_model is None: 
     8    warnings.warn("unable to load opencl; using ctypes instead") 
     9    from .core import load_model_dll as load_model 
    6210 
    6311class DirectModel: 

  • sasmodels/generate.py

    rae7b97b r3c56da87  
    187187# TODO: identify model files which have changed since loading and reload them. 
    188188 
    189 __all__ = ["make, doc", "sources", "use_single"] 
     189__all__ = ["make", "doc", "sources", "use_single"] 
    190190 
    191191import sys 
    192 import os 
    193 import os.path 
     192from os.path import abspath, dirname, join as joinpath, exists 
    194193import re 
    195194 
    196195import numpy as np 
     196C_KERNEL_TEMPLATE_PATH = joinpath(dirname(__file__), 'kernel_template.c') 
    197197 
    198198F64 = np.dtype('float64') 
    … …  
    201201# Scale and background, which are parameters common to every form factor 
    202202COMMON_PARAMETERS = [ 
    203     [ "scale", "", 1, [0, np.inf], "", "Source intensity" ], 
    204     [ "background", "1/cm", 0, [0, np.inf], "", "Source background" ], 
     203    ["scale", "", 1, [0, np.inf], "", "Source intensity"], 
     204    ["background", "1/cm", 0, [0, np.inf], "", "Source background"], 
    205205    ] 
    206206 
    … …  
    210210RST_UNITS = { 
    211211    "Ang": "|Ang|", 
     212    "1/Ang": "|Ang^-1|", 
    212213    "1/Ang^2": "|Ang^-2|", 
    213214    "1e-6/Ang^2": "|1e-6Ang^-2|", 
    … …  
    229230PARTABLE_VALUE_WIDTH = 10 
    230231 
    231 # Header included before every kernel. 
    232 # This makes sure that the appropriate math constants are defined, and does 
    233 # whatever is required to make the kernel compile as pure C rather than 
    234 # as an OpenCL kernel. 
    235 KERNEL_HEADER = """\ 
    236 // GENERATED CODE --- DO NOT EDIT --- 
    237 // Code is produced by sasmodels.gen from sasmodels/models/MODEL.c 
    238  
    239 #ifdef __OPENCL_VERSION__ 
    240 # define USE_OPENCL 
    241 #endif 
    242  
    243 // If opencl is not available, then we are compiling a C function 
    244 // Note: if using a C++ compiler, then define kernel as extern "C" 
    245 #ifndef USE_OPENCL 
    246 #  ifdef __cplusplus 
    247      #include <cmath> 
    248      #if defined(_MSC_VER) 
    249      #define kernel extern "C" __declspec( dllexport ) 
    250      #else 
    251      #define kernel extern "C" 
    252      #endif 
    253      using namespace std; 
    254      inline void SINCOS(double angle, double &svar, double &cvar) 
    255        { svar=sin(angle); cvar=cos(angle); } 
    256 #  else 
    257      #include <math.h> 
    258      #if defined(_MSC_VER) 
    259      #define kernel __declspec( dllexport ) 
    260      #else 
    261      #define kernel 
    262      #endif 
    263      #define SINCOS(angle,svar,cvar) do {svar=sin(angle);cvar=cos(angle);} while (0) 
    264 #  endif 
    265 #  define global 
    266 #  define local 
    267 #  define constant const 
    268 #  define powr(a,b) pow(a,b) 
    269 #else 
    270 #  ifdef USE_SINCOS 
    271 #    define SINCOS(angle,svar,cvar) svar=sincos(angle,&cvar) 
    272 #  else 
    273 #    define SINCOS(angle,svar,cvar) do {svar=sin(angle);cvar=cos(angle);} while (0) 
    274 #  endif 
    275 #endif 
    276  
    277 // Standard mathematical constants: 
    278 //   M_E, M_LOG2E, M_LOG10E, M_LN2, M_LN10, M_PI, M_PI_2=pi/2, M_PI_4=pi/4, 
    279 //   M_1_PI=1/pi, M_2_PI=2/pi, M_2_SQRTPI=2/sqrt(pi), SQRT2, SQRT1_2=sqrt(1/2) 
    280 // OpenCL defines M_constant_F for float constants, and nothing if double 
    281 // is not enabled on the card, which is why these constants may be missing 
    282 #ifndef M_PI 
    283 #  define M_PI 3.141592653589793 
    284 #endif 
    285 #ifndef M_PI_2 
    286 #  define M_PI_2 1.570796326794897 
    287 #endif 
    288 #ifndef M_PI_4 
    289 #  define M_PI_4 0.7853981633974483 
    290 #endif 
    291  
    292 // Non-standard pi/180, used for converting between degrees and radians 
    293 #ifndef M_PI_180 
    294 #  define M_PI_180 0.017453292519943295 
    295 #endif 
    296 """ 
    297  
    298  
    299 # The I(q) kernel and the I(qx, qy) kernel have one and two q parameters 
    300 # respectively, so the template builder will need to do extra work to 
    301 # declare, initialize and pass the q parameters. 
    302 KERNEL_1D = { 
    303     'fn': "Iq", 
    304     'q_par_decl': "global const double *q,", 
    305     'qinit': "const double qi = q[i];", 
    306     'qcall': "qi", 
    307     'qwork': ["q"], 
    308     } 
    309  
    310 KERNEL_2D = { 
    311     'fn': "Iqxy", 
    312     'q_par_decl': "global const double *qx,\n    global const double *qy,", 
    313     'qinit': "const double qxi = qx[i];\n    const double qyi = qy[i];", 
    314     'qcall': "qxi, qyi", 
    315     'qwork': ["qx", "qy"], 
    316     } 
    317  
    318 # Generic kernel template for the polydispersity loop. 
    319 # This defines the opencl kernel that is available to the host.  The same 
    320 # structure is used for Iq and Iqxy kernels, so extra flexibility is needed 
    321 # for q parameters.  The polydispersity loop is built elsewhere and 
    322 # substituted into this template. 
    323 KERNEL_TEMPLATE = """\ 
    324 kernel void %(name)s( 
    325     %(q_par_decl)s 
    326     global double *result, 
    327 #ifdef USE_OPENCL 
    328     global double *loops_g, 
    329 #else 
    330     const int Nq, 
    331 #endif 
    332     local double *loops, 
    333     const double cutoff, 
    334     %(par_decl)s 
    335     ) 
    336 { 
    337 #ifdef USE_OPENCL 
    338   // copy loops info to local memory 
    339   event_t e = async_work_group_copy(loops, loops_g, (%(pd_length)s)*2, 0); 
    340   wait_group_events(1, &e); 
    341  
    342   int i = get_global_id(0); 
    343   int Nq = get_global_size(0); 
    344 #endif 
    345  
    346 #ifdef USE_OPENCL 
    347   if (i < Nq) 
    348 #else 
    349   #pragma omp parallel for 
    350   for (int i=0; i < Nq; i++) 
    351 #endif 
    352   { 
    353     %(qinit)s 
    354     double ret=0.0, norm=0.0; 
    355     double vol=0.0, norm_vol=0.0; 
    356 %(loops)s 
    357     if (vol*norm_vol != 0.0) { 
    358       ret *= norm_vol/vol; 
    359     } 
    360     result[i] = scale*ret/norm+background; 
    361   } 
    362 } 
    363 """ 
    364  
    365 # Polydispersity loop level. 
    366 # This pulls the parameter value and weight from the looping vector in order 
    367 # in preperation for a nested loop. 
    368 LOOP_OPEN="""\ 
    369 for (int %(name)s_i=0; %(name)s_i < N%(name)s; %(name)s_i++) { 
    370   const double %(name)s = loops[2*(%(name)s_i%(offset)s)]; 
    371   const double %(name)s_w = loops[2*(%(name)s_i%(offset)s)+1];\ 
    372 """ 
    373  
    374  
    375  
    376 ########################################################## 
    377 #                                                        # 
    378 #   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!   # 
    379 #   !!                                              !!   # 
    380 #   !!  KEEP THIS CODE CONSISTENT WITH PYKERNEL.PY  !!   # 
    381 #   !!                                              !!   # 
    382 #   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!   # 
    383 #                                                        # 
    384 ########################################################## 
    385  
    386  
    387 # Polydispersity loop body. 
    388 # This computes the weight, and if it is sufficient, calls the scattering 
    389 # function and adds it to the total.  If there is a volume normalization, 
    390 # it will also be added here. 
    391 LOOP_BODY="""\ 
    392 const double weight = %(weight_product)s; 
    393 if (weight > cutoff) { 
    394   const double I = %(fn)s(%(qcall)s, %(pcall)s); 
    395   if (I>=0.0) { // scattering cannot be negative 
    396     ret += weight*I%(sasview_spherical)s; 
    397     norm += weight; 
    398     %(volume_norm)s 
    399   } 
    400   //else { printf("exclude qx,qy,I:%%g,%%g,%%g\\n",%(qcall)s,I); } 
    401 } 
    402 //else { printf("exclude weight:%%g\\n",weight); }\ 
    403 """ 
    404  
    405 # Use this when integrating over orientation 
    406 SPHERICAL_CORRECTION="""\ 
    407 // Correction factor for spherical integration p(theta) I(q) sin(theta) dtheta 
    408 double spherical_correction = (Ntheta>1 ? fabs(sin(M_PI_180*theta)) : 1.0);\ 
    409 """ 
    410 # Use this to reproduce sasview behaviour 
    411 SASVIEW_SPHERICAL_CORRECTION="""\ 
    412 // Correction factor for spherical integration p(theta) I(q) sin(theta) dtheta 
    413 double spherical_correction = (Ntheta>1 ? fabs(cos(M_PI_180*theta))*M_PI_2 : 1.0);\ 
    414 """ 
    415  
    416 # Volume normalization. 
    417 # If there are "volume" polydispersity parameters, then these will be used 
    418 # to call the form_volume function from the user supplied kernel, and accumulate 
    419 # a normalized weight. 
    420 VOLUME_NORM="""const double vol_weight = %(vol_weight)s; 
    421     vol += vol_weight*form_volume(%(vol_pars)s); 
    422     norm_vol += vol_weight;\ 
    423 """ 
    424  
    425 # functions defined as strings in the .py module 
    426 WORK_FUNCTION="""\ 
    427 double %(name)s(%(pars)s); 
    428 double %(name)s(%(pars)s) 
    429 { 
    430 %(body)s 
    431 }\ 
    432 """ 
    433  
    434232# Documentation header for the module, giving the model name, its short 
    435233# description and its parameter table.  The remainder of the doc comes 
    436234# from the module docstring. 
    437 DOC_HEADER=""".. _%(name)s: 
     235DOC_HEADER = """.. _%(name)s: 
    438236 
    439237%(label)s 
    … …  
    449247""" 
    450248 
    451 def indent(s, depth): 
    452     """ 
    453     Indent a string of text with *depth* additional spaces on each line. 
    454     """ 
    455     spaces = " "*depth 
    456     sep = "\n"+spaces 
    457     return spaces + sep.join(s.split("\n")) 
    458  
    459  
    460 def kernel_name(info, is_2D): 
    461     """ 
    462     Name of the exported kernel symbol. 
    463     """ 
    464     return info['name'] + "_" + ("Iqxy" if is_2D else "Iq") 
    465  
     249def make_partable(pars): 
     250    """ 
     251    Generate the parameter table to include in the sphinx documentation. 
     252    """ 
     253    pars = COMMON_PARAMETERS + pars 
     254    column_widths = [ 
     255        max(len(p[0]) for p in pars), 
     256        max(len(p[-1]) for p in pars), 
     257        max(len(RST_UNITS[p[1]]) for p in pars), 
     258        PARTABLE_VALUE_WIDTH, 
     259        ] 
     260    column_widths = [max(w, len(h)) 
     261                     for w, h in zip(column_widths, PARTABLE_HEADERS)] 
     262 
     263    sep = " ".join("="*w for w in column_widths) 
     264    lines = [ 
     265        sep, 
     266        " ".join("%-*s" % (w, h) for w, h in zip(column_widths, PARTABLE_HEADERS)), 
     267        sep, 
     268        ] 
     269    for p in pars: 
     270        lines.append(" ".join([ 
     271            "%-*s" % (column_widths[0], p[0]), 
     272            "%-*s" % (column_widths[1], p[-1]), 
     273            "%-*s" % (column_widths[2], RST_UNITS[p[1]]), 
     274            "%*g" % (column_widths[3], p[2]), 
     275            ])) 
     276    lines.append(sep) 
     277    return "\n".join(lines) 
     278 
     279def _search(search_path, filename): 
     280    """ 
     281    Find *filename* in *search_path*. 
     282 
     283    Raises ValueError if file does not exist. 
     284    """ 
     285    for path in search_path: 
     286        target = joinpath(path, filename) 
     287        if exists(target): 
     288            return target 
     289    raise ValueError("%r not found in %s" % (filename, search_path)) 
     290 
     291def sources(info): 
     292    """ 
     293    Return a list of the sources file paths for the module. 
     294    """ 
     295    search_path = [dirname(info['filename']), 
     296                   abspath(joinpath(dirname(__file__), 'models'))] 
     297    return [_search(search_path, f) for f in info['source']] 
    466298 
    467299def use_single(source): 
    … …  
    481313 
    482314 
    483 def make_kernel(info, is_2D): 
    484     """ 
    485     Build a kernel call from metadata supplied by the user. 
    486  
    487     *info* is the json object defined in the kernel file. 
    488  
    489     *form* is either "Iq" or "Iqxy". 
    490  
    491     This does not create a complete OpenCL kernel source, only the top 
    492     level kernel call with polydispersity and a call to the appropriate 
    493     Iq or Iqxy function. 
    494     """ 
    495  
    496     # If we are building the Iqxy kernel, we need to propagate qx,qy 
    497     # parameters, otherwise we can 
    498     dim = "2d" if is_2D else "1d" 
    499     fixed_pars = info['partype']['fixed-'+dim] 
    500     pd_pars = info['partype']['pd-'+dim] 
    501     vol_pars = info['partype']['volume'] 
    502     q_pars = KERNEL_2D if is_2D else KERNEL_1D 
    503     fn = q_pars['fn'] 
    504  
    505     # Build polydispersity loops 
    506     depth = 4 
    507     offset = "" 
    508     loop_head = [] 
    509     loop_end = [] 
    510     for name in pd_pars: 
    511         subst = { 'name': name, 'offset': offset } 
    512         loop_head.append(indent(LOOP_OPEN%subst, depth)) 
    513         loop_end.insert(0, (" "*depth) + "}") 
    514         offset += '+N'+name 
    515         depth += 2 
    516  
    517     # The volume parameters in the inner loop are used to call the volume() 
    518     # function in the kernel, with the parameters defined in vol_pars and the 
    519     # weight product defined in weight.  If there are no volume parameters, 
    520     # then there will be no volume normalization. 
    521     if vol_pars: 
    522         subst = { 
    523             'vol_weight': "*".join(p+"_w" for p in vol_pars), 
    524             'vol_pars': ", ".join(vol_pars), 
    525             } 
    526         volume_norm = VOLUME_NORM%subst 
    527     else: 
    528         volume_norm = "" 
    529  
    530     # Define the inner loop function call 
    531     # The parameters to the f(q,p1,p2...) call should occur in the same 
    532     # order as given in the parameter info structure.  This may be different 
    533     # from the parameter order in the call to the kernel since the kernel 
    534     # call places all fixed parameters before all polydisperse parameters. 
    535     fq_pars = [p[0] for p in info['parameters'][len(COMMON_PARAMETERS):] 
    536                if p[0] in set(fixed_pars+pd_pars)] 
    537     if False and "theta" in pd_pars: 
    538         spherical_correction = [indent(SPHERICAL_CORRECTION, depth)] 
    539         weights = [p+"_w" for p in pd_pars]+['spherical_correction'] 
    540         sasview_spherical = "" 
    541     elif True and "theta" in pd_pars: 
    542         spherical_correction = [indent(SASVIEW_SPHERICAL_CORRECTION,depth)] 
    543         weights = [p+"_w" for p in pd_pars] 
    544         sasview_spherical = "*spherical_correction" 
    545     else: 
    546         spherical_correction = [] 
    547         weights = [p+"_w" for p in pd_pars] 
    548         sasview_spherical = "" 
    549     subst = { 
    550         'weight_product': "*".join(weights), 
    551         'volume_norm': volume_norm, 
    552         'fn': fn, 
    553         'qcall': q_pars['qcall'], 
    554         'pcall': ", ".join(fq_pars), # skip scale and background 
    555         'sasview_spherical': sasview_spherical, 
    556         } 
    557     loop_body = [indent(LOOP_BODY%subst, depth)] 
    558     loops = "\n".join(loop_head+spherical_correction+loop_body+loop_end) 
    559  
    560     # declarations for non-pd followed by pd pars 
    561     # e.g., 
    562     #     const double sld, 
    563     #     const int Nradius 
    564     fixed_par_decl = ",\n    ".join("const double %s"%p for p in fixed_pars) 
    565     pd_par_decl = ",\n    ".join("const int N%s"%p for p in pd_pars) 
    566     if fixed_par_decl and pd_par_decl: 
    567         par_decl = ",\n    ".join((fixed_par_decl, pd_par_decl)) 
    568     elif fixed_par_decl: 
    569         par_decl = fixed_par_decl 
    570     else: 
    571         par_decl = pd_par_decl 
    572  
    573     # Finally, put the pieces together in the kernel. 
    574     subst = { 
    575         # kernel name is, e.g., cylinder_Iq 
    576         'name': kernel_name(info, is_2D), 
    577         # to declare, e.g., global double q[], 
    578         'q_par_decl': q_pars['q_par_decl'], 
    579         # to declare, e.g., double sld, int Nradius, int Nlength 
    580         'par_decl': par_decl, 
    581         # to copy global to local pd pars we need, e.g., Nradius+Nlength 
    582         'pd_length': "+".join('N'+p for p in pd_pars), 
    583         # the q initializers, e.g., double qi = q[i]; 
    584         'qinit': q_pars['qinit'], 
    585         # the actual polydispersity loop 
    586         'loops': loops, 
    587         } 
    588     kernel = KERNEL_TEMPLATE%subst 
    589  
    590     # If the working function is defined in the kernel metadata as a 
    591     # string, translate the string to an actual function definition 
    592     # and put it before the kernel. 
    593     if info[fn]: 
    594         subst = { 
    595             'name': fn, 
    596             'pars': ", ".join("double "+p for p in q_pars['qwork']+fq_pars), 
    597             'body': info[fn], 
    598             } 
    599         kernel = "\n".join((WORK_FUNCTION%subst, kernel)) 
    600     return kernel 
    601  
    602 def make_partable(pars): 
    603     """ 
    604     Generate the parameter table to include in the sphinx documentation. 
    605     """ 
    606     pars = COMMON_PARAMETERS + pars 
    607     column_widths = [ 
    608         max(len(p[0]) for p in pars), 
    609         max(len(p[-1]) for p in pars), 
    610         max(len(RST_UNITS[p[1]]) for p in pars), 
    611         PARTABLE_VALUE_WIDTH, 
    612         ] 
    613     column_widths = [max(w, len(h)) 
    614                      for w,h in zip(column_widths, PARTABLE_HEADERS)] 
    615  
    616     sep = " ".join("="*w for w in column_widths) 
    617     lines = [ 
    618         sep, 
    619         " ".join("%-*s"%(w,h) for w,h in zip(column_widths, PARTABLE_HEADERS)), 
    620         sep, 
    621         ] 
    622     for p in pars: 
    623         lines.append(" ".join([ 
    624             "%-*s"%(column_widths[0],p[0]), 
    625             "%-*s"%(column_widths[1],p[-1]), 
    626             "%-*s"%(column_widths[2],RST_UNITS[p[1]]), 
    627             "%*g"%(column_widths[3],p[2]), 
    628             ])) 
    629     lines.append(sep) 
    630     return "\n".join(lines) 
    631  
    632 def _search(search_path, filename): 
    633     """ 
    634     Find *filename* in *search_path*. 
    635  
    636     Raises ValueError if file does not exist. 
    637     """ 
    638     for path in search_path: 
    639         target = os.path.join(path, filename) 
    640         if os.path.exists(target): 
    641             return target 
    642     raise ValueError("%r not found in %s"%(filename, search_path)) 
    643  
    644 def sources(info): 
    645     """ 
    646     Return a list of the sources file paths for the module. 
    647     """ 
    648     from os.path import abspath, dirname, join as joinpath 
    649     search_path = [ dirname(info['filename']), 
    650                     abspath(joinpath(dirname(__file__),'models')) ] 
    651     return [_search(search_path, f) for f in info['source']] 
    652  
    653 def make_model(info): 
    654     """ 
    655     Generate the code for the kernel defined by info, using source files 
    656     found in the given search path. 
    657     """ 
    658     source = [open(f).read() for f in sources(info)] 
    659     # If the form volume is defined as a string, then wrap it in a 
    660     # function definition and place it after the external sources but 
    661     # before the kernel functions.  If the kernel functions are strings, 
    662     # they will be translated in the make_kernel call. 
    663     if info['form_volume']: 
    664         subst = { 
    665             'name': "form_volume", 
    666             'pars': ", ".join("double "+p for p in info['partype']['volume']), 
    667             'body': info['form_volume'], 
    668             } 
    669         source.append(WORK_FUNCTION%subst) 
    670     kernel_Iq = make_kernel(info, is_2D=False) if not callable(info['Iq']) else "" 
    671     kernel_Iqxy = make_kernel(info, is_2D=True) if not callable(info['Iqxy']) else "" 
    672     kernel = "\n\n".join([KERNEL_HEADER]+source+[kernel_Iq, kernel_Iqxy]) 
    673     return kernel 
     315def kernel_name(info, is_2D): 
     316    """ 
     317    Name of the exported kernel symbol. 
     318    """ 
     319    return info['name'] + "_" + ("Iqxy" if is_2D else "Iq") 
     320 
    674321 
    675322def categorize_parameters(pars): 
    … …  
    686333 
    687334    for p in pars: 
    688         name,ptype = p[0],p[4] 
     335        name, ptype = p[0], p[4] 
    689336        if ptype == 'volume': 
    690337            partype['pd-1d'].append(name) 
    … …  
    699346            partype['fixed-2d'].append(name) 
    700347        else: 
    701             raise ValueError("unknown parameter type %r"%ptype) 
     348            raise ValueError("unknown parameter type %r" % ptype) 
    702349        partype[ptype].append(name) 
    703350 
    704351    return partype 
     352 
     353def indent(s, depth): 
     354    """ 
     355    Indent a string of text with *depth* additional spaces on each line. 
     356    """ 
     357    spaces = " "*depth 
     358    sep = "\n" + spaces 
     359    return spaces + sep.join(s.split("\n")) 
     360 
     361 
     362def build_polydispersity_loops(pd_pars): 
     363    """ 
     364    Build polydispersity loops 
     365 
     366    Returns loop opening and loop closing 
     367    """ 
     368    LOOP_OPEN = """\ 
     369for (int %(name)s_i=0; %(name)s_i < N%(name)s; %(name)s_i++) { 
     370  const double %(name)s = loops[2*(%(name)s_i%(offset)s)]; 
     371  const double %(name)s_w = loops[2*(%(name)s_i%(offset)s)+1];\ 
     372""" 
     373    depth = 4 
     374    offset = "" 
     375    loop_head = [] 
     376    loop_end = [] 
     377    for name in pd_pars: 
     378        subst = {'name': name, 'offset': offset} 
     379        loop_head.append(indent(LOOP_OPEN % subst, depth)) 
     380        loop_end.insert(0, (" "*depth) + "}") 
     381        offset += '+N' + name 
     382        depth += 2 
     383    return "\n".join(loop_head), "\n".join(loop_end) 
     384 
     385C_KERNEL_TEMPLATE = None 
     386def make_model(info): 
     387    """ 
     388    Generate the code for the kernel defined by info, using source files 
     389    found in the given search path. 
     390    """ 
     391    # TODO: need something other than volume to indicate dispersion parameters 
     392    # No volume normalization despite having a volume parameter. 
     393    # Thickness is labelled a volume in order to trigger polydispersity. 
     394    # May want a separate dispersion flag, or perhaps a separate category for 
     395    # disperse, but not volume.  Volume parameters also use relative values 
     396    # for the distribution rather than the absolute values used by angular 
     397    # dispersion.  Need to be careful that necessary parameters are available 
     398    # for computing volume even if we allow non-disperse volume parameters. 
     399 
     400    # Load template 
     401    global C_KERNEL_TEMPLATE 
     402    if C_KERNEL_TEMPLATE is None: 
     403        with open(C_KERNEL_TEMPLATE_PATH) as fid: 
     404            C_KERNEL_TEMPLATE = fid.read() 
     405 
     406    # Load additional sources 
     407    source = [open(f).read() for f in sources(info)] 
     408 
     409    # Prepare defines 
     410    defines = [] 
     411    partype = info['partype'] 
     412    pd_1d = partype['pd-1d'] 
     413    pd_2d = partype['pd-2d'] 
     414    fixed_1d = partype['fixed-1d'] 
     415    fixed_2d = partype['fixed-1d'] 
     416 
     417    iq_parameters = [p[0] 
     418                     for p in info['parameters'][2:] # skip scale, background 
     419                     if p[0] in set(fixed_1d + pd_1d)] 
     420    iqxy_parameters = [p[0] 
     421                       for p in info['parameters'][2:] # skip scale, background 
     422                       if p[0] in set(fixed_2d + pd_2d)] 
     423    volume_parameters = [p[0] 
     424                         for p in info['parameters'] 
     425                         if p[4] == 'volume'] 
     426 
     427    # Fill in defintions for volume parameters 
     428    if volume_parameters: 
     429        defines.append(('VOLUME_PARAMETERS', 
     430                        ','.join(volume_parameters))) 
     431        defines.append(('VOLUME_WEIGHT_PRODUCT', 
     432                        '*'.join(p + '_w' for p in volume_parameters))) 
     433 
     434    # Generate form_volume function from body only 
     435    if info['form_volume'] is not None: 
     436        if volume_parameters: 
     437            vol_par_decl = ', '.join('double ' + p for p in volume_parameters) 
     438        else: 
     439            vol_par_decl = 'void' 
     440        defines.append(('VOLUME_PARAMETER_DECLARATIONS', 
     441                        vol_par_decl)) 
     442        fn = """\ 
     443double form_volume(VOLUME_PARAMETER_DECLARATIONS); 
     444double form_volume(VOLUME_PARAMETER_DECLARATIONS) { 
     445    %(body)s 
     446} 
     447""" % {'body':info['form_volume']} 
     448        source.append(fn) 
     449 
     450    # Fill in definitions for Iq parameters 
     451    defines.append(('IQ_KERNEL_NAME', info['name'] + '_Iq')) 
     452    defines.append(('IQ_PARAMETERS', ', '.join(iq_parameters))) 
     453    if fixed_1d: 
     454        defines.append(('IQ_FIXED_PARAMETER_DECLARATIONS', 
     455                        ', \\\n    '.join('const double %s' % p for p in fixed_1d))) 
     456    if pd_1d: 
     457        defines.append(('IQ_WEIGHT_PRODUCT', 
     458                        '*'.join(p + '_w' for p in pd_1d))) 
     459        defines.append(('IQ_DISPERSION_LENGTH_DECLARATIONS', 
     460                        ', \\\n    '.join('const int N%s' % p for p in pd_1d))) 
     461        defines.append(('IQ_DISPERSION_LENGTH_SUM', 
     462                        '+'.join('N' + p for p in pd_1d))) 
     463        open_loops, close_loops = build_polydispersity_loops(pd_1d) 
     464        defines.append(('IQ_OPEN_LOOPS', 
     465                        open_loops.replace('\n', ' \\\n'))) 
     466        defines.append(('IQ_CLOSE_LOOPS', 
     467                        close_loops.replace('\n', ' \\\n'))) 
     468    if info['Iq'] is not None: 
     469        defines.append(('IQ_PARAMETER_DECLARATIONS', 
     470                        ', '.join('double ' + p for p in iq_parameters))) 
     471        fn = """\ 
     472double Iq(double q, IQ_PARAMETER_DECLARATIONS); 
     473double Iq(double q, IQ_PARAMETER_DECLARATIONS) { 
     474    %(body)s 
     475} 
     476""" % {'body':info['Iq']} 
     477        source.append(fn) 
     478 
     479    # Fill in definitions for Iqxy parameters 
     480    defines.append(('IQXY_KERNEL_NAME', info['name'] + '_Iqxy')) 
     481    defines.append(('IQXY_PARAMETERS', ', '.join(iqxy_parameters))) 
     482    if fixed_2d: 
     483        defines.append(('IQXY_FIXED_PARAMETER_DECLARATIONS', 
     484                        ', \\\n    '.join('const double %s' % p for p in fixed_2d))) 
     485    if pd_2d: 
     486        defines.append(('IQXY_WEIGHT_PRODUCT', 
     487                        '*'.join(p + '_w' for p in pd_2d))) 
     488        defines.append(('IQXY_DISPERSION_LENGTH_DECLARATIONS', 
     489                        ', \\\n    '.join('const int N%s' % p for p in pd_2d))) 
     490        defines.append(('IQXY_DISPERSION_LENGTH_SUM', 
     491                        '+'.join('N' + p for p in pd_2d))) 
     492        open_loops, close_loops = build_polydispersity_loops(pd_2d) 
     493        defines.append(('IQXY_OPEN_LOOPS', 
     494                        open_loops.replace('\n', ' \\\n'))) 
     495        defines.append(('IQXY_CLOSE_LOOPS', 
     496                        close_loops.replace('\n', ' \\\n'))) 
     497    if info['Iqxy'] is not None: 
     498        defines.append(('IQXY_PARAMETER_DECLARATIONS', 
     499                        ', '.join('double ' + p for p in iqxy_parameters))) 
     500        fn = """\ 
     501double Iqxy(double qx, double qy, IQXY_PARAMETER_DECLARATIONS); 
     502double Iqxy(double qx, double qy, IQXY_PARAMETER_DECLARATIONS) { 
     503    %(body)s 
     504} 
     505""" % {'body':info['Iqxy']} 
     506        source.append(fn) 
     507 
     508    # Need to know if we have a theta parameter for Iqxy; it is not there 
     509    # for the magnetic sphere model, for example, which has a magnetic 
     510    # orientation but no shape orientation. 
     511    if 'theta' in pd_2d: 
     512        defines.append(('IQXY_HAS_THETA', '1')) 
     513 
     514    #for d in defines: print d 
     515    DEFINES = '\n'.join('#define %s %s' % (k, v) for k, v in defines) 
     516    SOURCES = '\n\n'.join(source) 
     517    return C_KERNEL_TEMPLATE % { 
     518        'DEFINES':DEFINES, 
     519        'SOURCES':SOURCES, 
     520        } 
    705521 
    706522def make(kernel_module): 
    … …  
    713529    """ 
    714530    # TODO: allow Iq and Iqxy to be defined in python 
    715     from os.path import abspath 
    716531    #print kernelfile 
    717532    info = dict( 
    718         filename = abspath(kernel_module.__file__), 
    719         name = kernel_module.name, 
    720         title = kernel_module.title, 
    721         description = kernel_module.description, 
    722         parameters = COMMON_PARAMETERS + kernel_module.parameters, 
    723         source = getattr(kernel_module, 'source', []), 
    724         oldname = kernel_module.oldname, 
    725         oldpars = kernel_module.oldpars, 
     533        filename=abspath(kernel_module.__file__), 
     534        name=kernel_module.name, 
     535        title=kernel_module.title, 
     536        description=kernel_module.description, 
     537        parameters=COMMON_PARAMETERS + kernel_module.parameters, 
     538        source=getattr(kernel_module, 'source', []), 
     539        oldname=kernel_module.oldname, 
     540        oldpars=kernel_module.oldpars, 
    726541        ) 
    727542    # Fill in attributes which default to None 
    728     info.update((k,getattr(kernel_module, k, None)) 
     543    info.update((k, getattr(kernel_module, k, None)) 
    729544                for k in ('ER', 'VR', 'form_volume', 'Iq', 'Iqxy')) 
    730545    # Fill in the derived attributes 
    731     info['limits'] = dict((p[0],p[3]) for p in info['parameters']) 
     546    info['limits'] = dict((p[0], p[3]) for p in info['parameters']) 
    732547    info['partype'] = categorize_parameters(info['parameters']) 
    733     info['defaults'] = dict((p[0],p[2]) for p in info['parameters']) 
    734  
    735     source = make_model(info) 
    736  
     548    info['defaults'] = dict((p[0], p[2]) for p in info['parameters']) 
     549 
     550    # Assume if one part of the kernel is python then all parts are. 
     551    source = make_model(info) if not callable(info['Iq']) else None 
    737552    return source, info 
    738553 
    … …  
    741556    Return the documentation for the model. 
    742557    """ 
    743     subst = dict(name=kernel_module.name.replace('_','-'), 
     558    subst = dict(name=kernel_module.name.replace('_', '-'), 
    744559                 label=" ".join(kernel_module.name.split('_')).capitalize(), 
    745560                 title=kernel_module.title, 
    746561                 parameters=make_partable(kernel_module.parameters), 
    747562                 docs=kernel_module.__doc__) 
    748     return DOC_HEADER%subst 
     563    return DOC_HEADER % subst 
    749564 
    750565 
    751566 
    752567def demo_time(): 
     568    from .models import cylinder 
    753569    import datetime 
    754     from .models import cylinder 
    755     toc = lambda: (datetime.datetime.now()-tic).total_seconds() 
    756570    tic = datetime.datetime.now() 
    757     source, info = make(cylinder) 
    758     print "time:",toc() 
     571    make(cylinder) 
     572    toc = (datetime.datetime.now() - tic).total_seconds() 
     573    print "time:", toc 
    759574 
    760575def main(): 
    … …  
    764579        name = sys.argv[1] 
    765580        import sasmodels.models 
    766         __import__('sasmodels.models.'+name) 
     581        __import__('sasmodels.models.' + name) 
    767582        model = getattr(sasmodels.models, name) 
    768         source, info = make(model); print source 
    769         #print doc(model) 
     583        source, _ = make(model) 
     584        print source 
    770585 
    771586if __name__ == "__main__": 

  • sasmodels/kernelcl.py

    rda32ec3 r3c56da87  
    3030try: 
    3131    import pyopencl as cl 
    32 except ImportError,exc: 
     32    # Ask OpenCL for the default context so that we know that one exists 
     33    cl.create_some_context(interactive=False) 
     34except Exception, exc: 
    3335    warnings.warn(str(exc)) 
    3436    raise RuntimeError("OpenCL not available") 
    3537 
    36 try: 
    37     context = cl.create_some_context(interactive=False) 
    38     del context 
    39 except cl.RuntimeError, exc: 
    40     warnings.warn(str(exc)) 
    41     raise RuntimeError("OpenCl not available") 
    42  
    4338from pyopencl import mem_flags as mf 
    4439 
    4540from . import generate 
    46 from .kernelpy import PyInput, PyKernel 
     41from .kernelpy import PyModel 
    4742 
    4843F64_DEFS = """\ 
    … …  
    6863    """ 
    6964    source, info = generate.make(kernel_module) 
     65    if callable(info.get('Iq', None)): 
     66        return PyModel(info) 
    7067    ## for debugging, save source to a .cl file, edit it, and reload as model 
    7168    #open(info['name']+'.cl','w').write(source) 
    … …  
    113110    device.min_data_type_align_size//4. 
    114111    """ 
    115     remainder = vector.size%boundary 
     112    remainder = vector.size % boundary 
    116113    if remainder != 0: 
    117114        size = vector.size + (boundary - remainder) 
    118         vector = np.hstack((vector, [extra]*(size-vector.size))) 
     115        vector = np.hstack((vector, [extra] * (size - vector.size))) 
    119116    return np.ascontiguousarray(vector, dtype=dtype) 
    120117 
    … …  
    129126    """ 
    130127    dtype = np.dtype(dtype) 
    131     if dtype==generate.F64 and not all(has_double(d) for d in context.devices): 
     128    if dtype == generate.F64 and not all(has_double(d) for d in context.devices): 
    132129        raise RuntimeError("Double precision not supported for devices") 
    133130 
    … …  
    138135    if context.devices[0].type == cl.device_type.GPU: 
    139136        header += "#define USE_SINCOS\n" 
    140     program  = cl.Program(context, header+source).build() 
     137    program = cl.Program(context, header + source).build() 
    141138    return program 
    142139 
    … …  
    163160 
    164161        if not self.context: 
    165             self.context = self._find_context() 
     162            self.context = _get_default_context() 
    166163 
    167164        # Byte boundary for data alignment 
    … …  
    176173        try: 
    177174            self.context = cl.create_some_context(interactive=False) 
    178         except Exception,exc: 
     175        except Exception, exc: 
    179176            warnings.warn(str(exc)) 
    180177            warnings.warn("pyopencl.create_some_context() failed") 
    181178            warnings.warn("the environment variable 'PYOPENCL_CTX' might not be set correctly") 
    182  
    183     def _find_context(self): 
    184         default = None 
    185         for platform in cl.get_platforms(): 
    186             for device in platform.get_devices(): 
    187                 if device.type == cl.device_type.GPU: 
    188                     return cl.Context([device]) 
    189                 if default is None: 
    190                     default = device 
    191  
    192         if not default: 
    193             raise RuntimeError("OpenCL device not found") 
    194  
    195         return cl.Context([default]) 
    196179 
    197180    def compile_program(self, name, source, dtype): 
    … …  
    206189            del self.compiled[name] 
    207190 
     191def _get_default_context(): 
     192    default = None 
     193    for platform in cl.get_platforms(): 
     194        for device in platform.get_devices(): 
     195            if device.type == cl.device_type.GPU: 
     196                return cl.Context([device]) 
     197            if default is None: 
     198                default = device 
     199 
     200    if not default: 
     201        raise RuntimeError("OpenCL device not found") 
     202 
     203    return cl.Context([default]) 
     204 
    208205 
    209206class GpuModel(object): 
    … …  
    233230        self.__dict__ = state.copy() 
    234231 
    235     def __call__(self, input): 
    236         # Support pure python kernel call 
    237         if input.is_2D and callable(self.info['Iqxy']): 
    238             return PyKernel(self.info['Iqxy'], self.info, input) 
    239         elif not input.is_2D and callable(self.info['Iq']): 
    240             return PyKernel(self.info['Iq'], self.info, input) 
    241  
    242         if self.dtype != input.dtype: 
     232    def __call__(self, input_value): 
     233        if self.dtype != input_value.dtype: 
    243234            raise TypeError("data and kernel have different types") 
    244235        if self.program is None: 
    245             self.program = environment().compile_program(self.info['name'],self.source, self.dtype) 
    246         kernel_name = generate.kernel_name(self.info, input.is_2D) 
     236            compiler = environment().compile_program 
     237            self.program = compiler(self.info['name'], self.source, self.dtype) 
     238        kernel_name = generate.kernel_name(self.info, input_value.is_2D) 
    247239        kernel = getattr(self.program, kernel_name) 
    248         return GpuKernel(kernel, self.info, input) 
     240        return GpuKernel(kernel, self.info, input_value) 
    249241 
    250242    def release(self): 
    … …  
    261253        ctypes and some may be pure python. 
    262254        """ 
    263         # Support pure python kernel call 
    264         if len(q_vectors) == 1 and callable(self.info['Iq']): 
    265             return PyInput(q_vectors, dtype=self.dtype) 
    266         elif callable(self.info['Iqxy']): 
    267             return PyInput(q_vectors, dtype=self.dtype) 
    268         else: 
    269             return GpuInput(q_vectors, dtype=self.dtype) 
     255        return GpuInput(q_vectors, dtype=self.dtype) 
    270256 
    271257# TODO: check that we don't need a destructor for buffers which go out of scope 
    … …  
    300286        self.q_vectors = [_stretch_input(q, self.dtype, 32) for q in q_vectors] 
    301287        self.q_buffers = [ 
    302             cl.Buffer(env.context,  mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=q) 
     288            cl.Buffer(env.context, mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=q) 
    303289            for q in self.q_vectors 
    304290        ] 
    … …  
    318304    *info* is the module information 
    319305 
    320     *input* is the DllInput q vectors at which the kernel should be 
     306    *q_input* is the DllInput q vectors at which the kernel should be 
    321307    evaluated. 
    322308 
    … …  
    329315    Call :meth:`release` when done with the kernel instance. 
    330316    """ 
    331     def __init__(self, kernel, info, input): 
    332         self.input = input 
     317    def __init__(self, kernel, info, q_input): 
     318        self.q_input = q_input 
    333319        self.kernel = kernel 
    334320        self.info = info 
    335         self.res = np.empty(input.nq, input.dtype) 
    336         dim = '2d' if input.is_2D else '1d' 
    337         self.fixed_pars = info['partype']['fixed-'+dim] 
    338         self.pd_pars = info['partype']['pd-'+dim] 
     321        self.res = np.empty(q_input.nq, q_input.dtype) 
     322        dim = '2d' if q_input.is_2D else '1d' 
     323        self.fixed_pars = info['partype']['fixed-' + dim] 
     324        self.pd_pars = info['partype']['pd-' + dim] 
    339325 
    340326        # Inputs and outputs for each kernel call 
    … …  
    342328        env = environment() 
    343329        self.loops_b = [cl.Buffer(env.context, mf.READ_WRITE, 
    344                                   2*MAX_LOOPS*input.dtype.itemsize) 
     330                                  2 * MAX_LOOPS * q_input.dtype.itemsize) 
    345331                        for _ in env.queues] 
    346332        self.res_b = [cl.Buffer(env.context, mf.READ_WRITE, 
    347                                 input.global_size[0]*input.dtype.itemsize) 
     333                                q_input.global_size[0] * q_input.dtype.itemsize) 
    348334                      for _ in env.queues] 
    349335 
    350336 
    351     def __call__(self, pars, pd_pars, cutoff=1e-5): 
    352         real = np.float32 if self.input.dtype == generate.F32 else np.float64 
    353         fixed = [real(p) for p in pars] 
    354         cutoff = real(cutoff) 
    355         loops = np.hstack(pd_pars) 
    356         loops = np.ascontiguousarray(loops.T, self.input.dtype).flatten() 
    357         Nloops = [np.uint32(len(p[0])) for p in pd_pars] 
    358         #print "loops",Nloops, loops 
    359  
    360         #import sys; print >>sys.stderr,"opencl eval",pars 
    361         #print "opencl eval",pars 
    362         if len(loops) > 2*MAX_LOOPS: 
    363             raise ValueError("too many polydispersity points") 
     337    def __call__(self, fixed_pars, pd_pars, cutoff=1e-5): 
     338        real = np.float32 if self.q_input.dtype == generate.F32 else np.float64 
     339 
    364340        device_num = 0 
     341        queuei = environment().queues[device_num] 
    365342        res_bi = self.res_b[device_num] 
    366         queuei = environment().queues[device_num] 
    367         loops_bi = self.loops_b[device_num] 
    368         loops_l = cl.LocalMemory(len(loops.data)) 
    369         cl.enqueue_copy(queuei, loops_bi, loops) 
    370         #ctx = environment().context 
    371         #loops_bi = cl.Buffer(ctx, mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=loops) 
    372         args = self.input.q_buffers + [res_bi,loops_bi,loops_l,cutoff] + fixed + Nloops 
    373         self.kernel(queuei, self.input.global_size, None, *args) 
     343        nq = np.uint32(self.q_input.nq) 
     344        if pd_pars: 
     345            cutoff = real(cutoff) 
     346            loops_N = [np.uint32(len(p[0])) for p in pd_pars] 
     347            loops = np.hstack(pd_pars) \ 
     348                if pd_pars else np.empty(0, dtype=self.q_input.dtype) 
     349            loops = np.ascontiguousarray(loops.T, self.q_input.dtype).flatten() 
     350            #print "loops",Nloops, loops 
     351 
     352            #import sys; print >>sys.stderr,"opencl eval",pars 
     353            #print "opencl eval",pars 
     354            if len(loops) > 2 * MAX_LOOPS: 
     355                raise ValueError("too many polydispersity points") 
     356 
     357            loops_bi = self.loops_b[device_num] 
     358            cl.enqueue_copy(queuei, loops_bi, loops) 
     359            loops_l = cl.LocalMemory(len(loops.data)) 
     360            #ctx = environment().context 
     361            #loops_bi = cl.Buffer(ctx, mf.READ_ONLY|mf.COPY_HOST_PTR, hostbuf=loops) 
     362            dispersed = [loops_bi, loops_l, cutoff] + loops_N 
     363        else: 
     364            dispersed = [] 
     365        fixed = [real(p) for p in fixed_pars] 
     366        args = self.q_input.q_buffers + [res_bi, nq] + dispersed + fixed 
     367        self.kernel(queuei, self.q_input.global_size, None, *args) 
    374368        cl.enqueue_copy(queuei, self.res, res_bi) 
    375369 

  • sasmodels/kerneldll.py

    r68d3c1b r3c56da87  
    1111 
    1212from . import generate 
    13 from .kernelpy import PyInput, PyKernel 
     13from .kernelpy import PyInput, PyModel 
    1414 
    1515from .generate import F32, F64 
    … …  
    1919    COMPILE = "gcc -shared -fPIC -std=c99 -O2 -Wall %(source)s -o %(output)s -lm" 
    2020elif os.name == 'nt': 
    21     # make sure vcvarsall.bat is called first in order to set compiler, headers, lib paths, etc. 
     21    # call vcvarsall.bat before compiling to set path, headers, libs, etc. 
    2222    if "VCINSTALLDIR" in os.environ: 
    2323        # MSVC compiler is available, so use it. 
     24        # TODO: remove intermediate OBJ file created in the directory 
     25        # TODO: maybe don't use randomized name for the c file 
    2426        COMPILE = "cl /nologo /Ox /MD /W3 /GS- /DNDEBUG /Tp%(source)s /openmp /link /DLL /INCREMENTAL:NO /MANIFEST /OUT:%(output)s" 
    25         # Can't find VCOMP90.DLL (don't know why), so remove openmp support from windows compiler build 
     27        # Can't find VCOMP90.DLL (don't know why), so remove openmp support 
     28        # from windows compiler build 
    2629        #COMPILE = "cl /nologo /Ox /MD /W3 /GS- /DNDEBUG /Tp%(source)s /link /DLL /INCREMENTAL:NO /MANIFEST /OUT:%(output)s" 
    2730    else: 
    … …  
    5457    be defined without using too many resources. 
    5558    """ 
    56     import tempfile 
    57  
    5859    source, info = generate.make(kernel_module) 
     60    if callable(info.get('Iq',None)): 
     61        return PyModel(info) 
    5962    source_files = generate.sources(info) + [info['filename']] 
    6063    newest = max(os.path.getmtime(f) for f in source_files) 
    … …  
    6871        status = os.system(command) 
    6972        if status != 0: 
    70             print "compile failed.  File is in %r"%filename 
     73            raise RuntimeError("compile failed.  File is in %r"%filename) 
    7174        else: 
    7275            ## uncomment the following to keep the generated c file 
    … …  
    7679 
    7780 
    78 IQ_ARGS = [c_void_p, c_void_p, c_int, c_void_p, c_double] 
    79 IQXY_ARGS = [c_void_p, c_void_p, c_void_p, c_int, c_void_p, c_double] 
     81IQ_ARGS = [c_void_p, c_void_p, c_int] 
     82IQXY_ARGS = [c_void_p, c_void_p, c_void_p, c_int] 
    8083 
    8184class DllModel(object): 
    … …  
    107110        self.dll = ct.CDLL(self.dllpath) 
    108111 
     112        pd_args_1d = [c_void_p, c_double] + [c_int]*Npd1d if Npd1d else [] 
     113        pd_args_2d= [c_void_p, c_double] + [c_int]*Npd2d if Npd2d else [] 
    109114        self.Iq = self.dll[generate.kernel_name(self.info, False)] 
    110         self.Iq.argtypes = IQ_ARGS + [c_double]*Nfixed1d + [c_int]*Npd1d 
     115        self.Iq.argtypes = IQ_ARGS + pd_args_1d + [c_double]*Nfixed1d 
    111116 
    112117        self.Iqxy = self.dll[generate.kernel_name(self.info, True)] 
    113         self.Iqxy.argtypes = IQXY_ARGS + [c_double]*Nfixed2d + [c_int]*Npd2d 
     118        self.Iqxy.argtypes = IQXY_ARGS + pd_args_2d + [c_double]*Nfixed2d 
    114119 
    115120    def __getstate__(self): 
    … …  
    119124        self.__dict__ = state 
    120125 
    121     def __call__(self, input): 
    122         # Support pure python kernel call 
    123         if input.is_2D and callable(self.info['Iqxy']): 
    124             return PyKernel(self.info['Iqxy'], self.info, input) 
    125         elif not input.is_2D and callable(self.info['Iq']): 
    126             return PyKernel(self.info['Iq'], self.info, input) 
     126    def __call__(self, q_input): 
     127        if self.dll is None: self._load_dll() 
     128        kernel = self.Iqxy if q_input.is_2D else self.Iq 
     129        return DllKernel(kernel, self.info, q_input) 
    127130 
    128         if self.dll is None: self._load_dll() 
    129         kernel = self.Iqxy if input.is_2D else self.Iq 
    130         return DllKernel(kernel, self.info, input) 
    131  
     131    # pylint: disable=no-self-use 
    132132    def make_input(self, q_vectors): 
    133133        """ 
    … …  
    152152    *info* is the module information 
    153153 
    154     *input* is the DllInput q vectors at which the kernel should be 
     154    *q_input* is the DllInput q vectors at which the kernel should be 
    155155    evaluated. 
    156156 
    … …  
    163163    Call :meth:`release` when done with the kernel instance. 
    164164    """ 
    165     def __init__(self, kernel, info, input): 
     165    def __init__(self, kernel, info, q_input): 
    166166        self.info = info 
    167         self.input = input 
     167        self.q_input = q_input 
    168168        self.kernel = kernel 
    169         self.res = np.empty(input.nq, input.dtype) 
    170         dim = '2d' if input.is_2D else '1d' 
     169        self.res = np.empty(q_input.nq, q_input.dtype) 
     170        dim = '2d' if q_input.is_2D else '1d' 
    171171        self.fixed_pars = info['partype']['fixed-'+dim] 
    172172        self.pd_pars = info['partype']['pd-'+dim] 
    … …  
    175175        self.p_res = self.res.ctypes.data 
    176176 
    177     def __call__(self, pars, pd_pars, cutoff): 
    178         real = np.float32 if self.input.dtype == F32 else np.float64 
    179         fixed = [real(p) for p in pars] 
    180         cutoff = real(cutoff) 
    181         loops = np.hstack(pd_pars) 
    182         loops = np.ascontiguousarray(loops.T, self.input.dtype).flatten() 
    183         loops_N = [np.uint32(len(p[0])) for p in pd_pars] 
     177    def __call__(self, fixed_pars, pd_pars, cutoff): 
     178        real = np.float32 if self.q_input.dtype == F32 else np.float64 
    184179 
    185         nq = c_int(self.input.nq) 
    186         p_loops = loops.ctypes.data 
    187         args = self.input.q_pointers + [self.p_res, nq, p_loops, cutoff] + fixed + loops_N 
     180        nq = c_int(self.q_input.nq) 
     181        if pd_pars: 
     182            cutoff = real(cutoff) 
     183            loops_N = [np.uint32(len(p[0])) for p in pd_pars] 
     184            loops = np.hstack(pd_pars) 
     185            loops = np.ascontiguousarray(loops.T, self.q_input.dtype).flatten() 
     186            p_loops = loops.ctypes.data 
     187            dispersed = [p_loops, cutoff] + loops_N 
     188        else: 
     189            dispersed = [] 
     190        fixed = [real(p) for p in fixed_pars] 
     191        args = self.q_input.q_pointers + [self.p_res, nq] + dispersed + fixed 
    188192        #print pars 
    189193        self.kernel(*args) 

  • sasmodels/kernelpy.py

    r6edb74a r3c56da87  
    11import numpy as np 
    2 from numpy import pi, sin, cos, sqrt 
    3  
    4 from .generate import F32, F64 
     2from numpy import pi, cos 
     3 
     4from .generate import F64 
     5 
     6class PyModel(object): 
     7    def __init__(self, info): 
     8        self.info = info 
     9 
     10    def __call__(self, q_input): 
     11        kernel = self.info['Iqxy'] if q_input.is_2D else self.info['Iq'] 
     12        return PyKernel(kernel, self.info, q_input) 
     13 
     14    # pylint: disable=no-self-use 
     15    def make_input(self, q_vectors): 
     16        return PyInput(q_vectors, dtype=F64) 
     17 
     18    def release(self): 
     19        pass 
    520 
    621class PyInput(object): 
    … …  
    2742        self.dtype = dtype 
    2843        self.is_2D = (len(q_vectors) == 2) 
    29         self.q_vectors = [np.ascontiguousarray(q,self.dtype) for q in q_vectors] 
     44        self.q_vectors = [np.ascontiguousarray(q, self.dtype) for q in q_vectors] 
    3045        self.q_pointers = [q.ctypes.data for q in q_vectors] 
    3146 
    … …  
    4156    *info* is the module information 
    4257 
    43     *input* is the DllInput q vectors at which the kernel should be 
     58    *q_input* is the DllInput q vectors at which the kernel should be 
    4459    evaluated. 
    4560 
    … …  
    5267    Call :meth:`release` when done with the kernel instance. 
    5368    """ 
    54     def __init__(self, kernel, info, input): 
     69    def __init__(self, kernel, info, q_input): 
    5570        self.info = info 
    56         self.input = input 
    57         self.res = np.empty(input.nq, input.dtype) 
    58         dim = '2d' if input.is_2D else '1d' 
     71        self.q_input = q_input 
     72        self.res = np.empty(q_input.nq, q_input.dtype) 
     73        dim = '2d' if q_input.is_2D else '1d' 
    5974        # Loop over q unless user promises that the kernel is vectorized by 
    6075        # taggining it with vectorized=True 
    … …  
    6277            if dim == '2d': 
    6378                def vector_kernel(qx, qy, *args): 
    64                     return np.array([kernel(qxi,qyi,*args) for qxi,qyi in zip(qx,qy)]) 
     79                    return np.array([kernel(qxi, qyi, *args) 
     80                                     for qxi, qyi in zip(qx, qy)]) 
    6581            else: 
    6682                def vector_kernel(q, *args): 
    67                     return np.array([kernel(qi,*args) for qi in q]) 
     83                    return np.array([kernel(qi, *args) 
     84                                     for qi in q]) 
    6885            self.kernel = vector_kernel 
    6986        else: 
    7087            self.kernel = kernel 
    71         fixed_pars = info['partype']['fixed-'+dim] 
    72         pd_pars = info['partype']['pd-'+dim] 
     88        fixed_pars = info['partype']['fixed-' + dim] 
     89        pd_pars = info['partype']['pd-' + dim] 
    7390        vol_pars = info['partype']['volume'] 
    7491 
    7592        # First two fixed pars are scale and background 
    7693        pars = [p[0] for p in info['parameters'][2:]] 
    77         offset = len(self.input.q_vectors) 
    78         self.args = self.input.q_vectors + [None]*len(pars) 
    79         self.fixed_index = np.array([pars.index(p)+offset for p in fixed_pars[2:] ]) 
    80         self.pd_index = np.array([pars.index(p)+offset for p in pd_pars]) 
    81         self.vol_index = np.array([pars.index(p)+offset for p in vol_pars]) 
    82         try: self.theta_index = pars.index('theta')+offset 
     94        offset = len(self.q_input.q_vectors) 
     95        self.args = self.q_input.q_vectors + [None] * len(pars) 
     96        self.fixed_index = np.array([pars.index(p) + offset 
     97                                     for p in fixed_pars[2:]]) 
     98        self.pd_index = np.array([pars.index(p) + offset 
     99                                  for p in pd_pars]) 
     100        self.vol_index = np.array([pars.index(p) + offset 
     101                                   for p in vol_pars]) 
     102        try: self.theta_index = pars.index('theta') + offset 
    83103        except ValueError: self.theta_index = -1 
    84104 
    … …  
    94114        # First two fixed 
    95115        scale, background = fixed[:2] 
    96         for index,value in zip(self.fixed_index, fixed[2:]): 
     116        for index, value in zip(self.fixed_index, fixed[2:]): 
    97117            args[index] = float(value) 
    98         res = _loops(form, form_volume, cutoff, scale, background,  args, 
     118        res = _loops(form, form_volume, cutoff, scale, background, args, 
    99119                     pd, self.pd_index, self.vol_index, self.theta_index) 
    100120 
    … …  
    102122 
    103123    def release(self): 
    104         self.input = None 
     124        self.q_input = None 
    105125 
    106126def _loops(form, form_volume, cutoff, scale, background, 
    … …  
    134154    """ 
    135155 
    136     ########################################################## 
    137     #                                                        # 
    138     #   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!   # 
    139     #   !!                                              !!   # 
    140     #   !!  KEEP THIS CODE CONSISTENT WITH GENERATE.PY  !!   # 
    141     #   !!                                              !!   # 
    142     #   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!   # 
    143     #                                                        # 
    144     ########################################################## 
     156    ################################################################ 
     157    #                                                              # 
     158    #   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!   # 
     159    #   !!                                                    !!   # 
     160    #   !!  KEEP THIS CODE CONSISTENT WITH KERNEL_TEMPLATE.C  !!   # 
     161    #   !!                                                    !!   # 
     162    #   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!   # 
     163    #                                                              # 
     164    ################################################################ 
    145165 
    146166    weight = np.empty(len(pd), 'd') 
    … …  
    164184        stride = np.array([1]) 
    165185        vol_weight_index = slice(None, None) 
     186        # keep lint happy 
     187        fast_value = [None] 
     188        fast_weight = [None] 
    166189 
    167190    ret = np.zeros_like(args[0]) 
    … …  
    171194    for k in range(stride[-1]): 
    172195        # update polydispersity parameter values 
    173         fast_index = k%stride[0] 
     196        fast_index = k % stride[0] 
    174197        if fast_index == 0:  # bottom loop complete ... check all other loops 
    175198            if weight.size > 0: 
    176                 for i,index, in enumerate(k%stride): 
     199                for i, index, in enumerate(k % stride): 
    177200                    args[pd_index[i]] = pd[i][0][index] 
    178201                    weight[i] = pd[i][1][index] 
    … …  
    180203            args[pd_index[0]] = fast_value[fast_index] 
    181204            weight[0] = fast_weight[fast_index] 
    182         # This computes the weight, and if it is sufficient, calls the scattering 
    183         # function and adds it to the total.  If there is a volume normalization, 
    184         # it will also be added here. 
     205        # This computes the weight, and if it is sufficient, calls the 
     206        # scattering function and adds it to the total.  If there is a volume 
     207        # normalization, it will also be added here. 
    185208        # Note: make sure this is consistent with the code in PY_LOOP_BODY!! 
    186209        # Note: can precompute w1*w2*...*wn 
    … …  
    188211        if w > cutoff: 
    189212            I = form(*args) 
    190             positive = (I>=0.0) 
    191  
    192             # Note: can precompute spherical correction if theta_index is not the fast index 
    193             # Correction factor for spherical integration p(theta) I(q) sin(theta) dtheta 
    194             #spherical_correction = abs(sin(pi*args[theta_index])) if theta_index>=0 else 1.0 
    195             #spherical_correction = abs(cos(pi*args[theta_index]))*pi/2 if theta_index>=0 else 1.0 
    196             spherical_correction = 1.0 
    197             ret += w*I*spherical_correction*positive 
    198             norm += w*positive 
     213            positive = (I >= 0.0) 
     214 
     215            # Note: can precompute spherical correction if theta_index is not 
     216            # the fast index. Correction factor for spherical integration 
     217            #spherical_correction = abs(cos(pi*args[phi_index])) if phi_index>=0 else 1.0 
     218            spherical_correction = (abs(cos(pi * args[theta_index])) * pi / 2 
     219                                    if theta_index >= 0 else 1.0) 
     220            #spherical_correction = 1.0 
     221            ret += w * I * spherical_correction * positive 
     222            norm += w * positive 
    199223 
    200224            # Volume normalization. 
    201             # If there are "volume" polydispersity parameters, then these will be used 
    202             # to call the form_volume function from the user supplied kernel, and accumulate 
    203             # a normalized weight. 
    204             # Note: can precompute volume norm if the fast index is not a volume index 
     225            # If there are "volume" polydispersity parameters, then these 
     226            # will be used to call the form_volume function from the user 
     227            # supplied kernel, and accumulate a normalized weight. 
     228            # Note: can precompute volume norm if fast index is not a volume 
    205229            if form_volume: 
    206230                vol_args = [args[index] for index in vol_index] 
    207231                vol_weight = np.prod(weight[vol_weight_index]) 
    208                 vol += vol_weight*form_volume(*vol_args)*positive 
    209                 vol_norm += vol_weight*positive 
    210  
    211     positive = (vol*vol_norm != 0.0) 
     232                vol += vol_weight * form_volume(*vol_args) * positive 
     233                vol_norm += vol_weight * positive 
     234 
     235    positive = (vol * vol_norm != 0.0) 
    212236    ret[positive] *= vol_norm[positive] / vol[positive] 
    213     result = scale*ret/norm+background 
     237    result = scale * ret / norm + background 
    214238    return result 

  • sasmodels/model_test.py

    r5428233 r3c56da87  
    11# -*- coding: utf-8 -*- 
    22""" 
    3 Created on Tue Feb 17 11:43:56 2015 
    4  
    5 @author: David 
     3Run model unit tests. 
     4 
     5Usage:: 
     6 
     7    python -m sasmodels.model_test [opencl|dll|opencl_and_dll] model1 model2 ... 
     8 
     9    if model1 is 'all', then all except the remaining models will be tested 
     10 
     11Each model is tested using the default parameters at q=0.1, (qx,qy)=(0.1,0.1), 
     12and the ER and VR are computed.  The return values at these points are not 
     13considered.  The test is only to verify that the models run to completion, 
     14and do not produce inf or NaN. 
     15 
     16Tests are defined with the *tests* attribute in the model.py file.  *tests* 
     17is a list of individual tests to run, where each test consists of the 
     18parameter values for the test, the q-values and the expected results.  For 
     19the effective radius test, the q-value should be 'ER'.  For the VR test, 
     20the q-value should be 'VR'.  For 1-D tests, either specify the q value or 
     21a list of q-values, and the corresponding I(q) value, or list of I(q) values. 
     22 
     23That is:: 
     24 
     25    tests = [ 
     26        [ {parameters}, q, I(q)], 
     27        [ {parameters}, [q], [I(q)] ], 
     28        [ {parameters}, [q1, q2, ...], [I(q1), I(q2), ...]], 
     29 
     30        [ {parameters}, (qx, qy), I(qx, Iqy)], 
     31        [ {parameters}, [(qx1, qy1), (qx2, qy2), ...], 
     32                        [I(qx1,qy1), I(qx2,qy2), ...]], 
     33 
     34        [ {parameters}, 'ER', ER(pars) ], 
     35        [ {parameters}, 'VR', VR(pars) ], 
     36        ... 
     37    ] 
     38 
     39Parameters are *key:value* pairs, where key is one of the parameters of the 
     40model and value is the value to use for the test.  Any parameters not given 
     41in the parameter list will take on the default parameter value. 
     42 
     43Precision defaults to 5 digits (relative). 
    644""" 
    745 
     46import sys 
    847import unittest 
    9 import warnings 
     48 
    1049import numpy as np 
    1150 
    12 from os.path import basename, dirname, join as joinpath 
    13 from glob import glob 
    14  
    15 try: 
    16     from .kernelcl import load_model 
    17 except ImportError,exc: 
    18     warnings.warn(str(exc)) 
    19     warnings.warn("using ctypes instead") 
    20     from .kerneldll import load_model 
    21  
    22 def load_kernel(model, dtype='single'):    
    23     kernel = load_model(model, dtype=dtype) 
    24     kernel.info['defaults'] = dict((p[0],p[2]) for p in kernel.info['parameters']) 
    25     return kernel 
    26  
    27 def get_weights(model, pars, name): 
    28     from . import weights 
    29      
    30     relative = name in model.info['partype']['pd-rel'] 
    31     disperser = pars.get(name+"_pd_type", "gaussian") 
    32     value = pars.get(name, model.info['defaults'][name]) 
    33     width = pars.get(name+"_pd", 0.0) 
    34     npts = pars.get(name+"_pd_n", 30) 
    35     nsigma = pars.get(name+"_pd_nsigma", 3.0) 
    36     v,w = weights.get_weights( 
    37             disperser, npts, width, nsigma, 
    38             value, model.info['limits'][name], relative) 
    39     return v,w/np.sum(w) 
    40  
    41 def eval_kernel(kernel, q, pars, cutoff=1e-5): 
    42     input = kernel.make_input(q) 
    43     finput = kernel(input) 
    44  
    45     fixed_pars = [pars.get(name, finput.info['defaults'][name]) 
    46                   for name in finput.fixed_pars] 
    47     pd_pars = [get_weights(finput, pars, p) for p in finput.pd_pars] 
    48     return finput(fixed_pars, pd_pars, cutoff) 
     51from .core import list_models, load_model_definition 
     52from .core import load_model_cl, load_model_dll 
     53from .core import make_kernel, call_kernel, call_ER, call_VR 
    4954 
    5055def annotate_exception(exc, msg): 
    … …  
    5560    Example:: 
    5661        >>> D = {} 
    57         >>> try:  
     62        >>> try: 
    5863        ...    print D['hello'] 
    59         ... except Exception,exc:  
     64        ... except Exception,exc: 
    6065        ...    annotate_exception(exc, "while accessing 'D'") 
    6166        ...    raise 
    … …  
    6671    args = exc.args 
    6772    if not args: 
    68         arg0 = msg 
     73        exc.args = (msg,) 
    6974    else: 
    70         arg0 = " ".join((args[0],msg)) 
    71     exc.args = tuple([arg0] + list(args[1:])) 
    72      
    73 def suite(): 
    74     root = dirname(__file__) 
    75     files = sorted(glob(joinpath(root, 'models', "[a-zA-Z]*.py"))) 
    76     models_names = [basename(f)[:-3] for f in files] 
    77      
     75        try: 
     76            arg0 = " ".join((args[0],msg)) 
     77            exc.args = tuple([arg0] + list(args[1:])) 
     78        except: 
     79            exc.args = (" ".join((str(exc),msg)),) 
     80 
     81def make_suite(loaders, models): 
     82 
     83    ModelTestCase = _hide_model_case_from_nosetests() 
    7884    suite = unittest.TestSuite() 
    79      
    80     for model_name in models_names: 
    81         module = __import__('sasmodels.models.' + model_name) 
    82         module = getattr(module, 'models', None) 
    83  
    84         model = getattr(module, model_name, None) 
    85         tests = getattr(model, 'tests', []) 
    86          
    87         if tests: 
     85 
     86    if models[0] == 'all': 
     87        skip = models[1:] 
     88        models = list_models() 
     89    else: 
     90        skip = [] 
     91    for model_name in models: 
     92        if model_name in skip: continue 
     93        model_definition = load_model_definition(model_name) 
     94 
     95        smoke_tests = [ 
     96            [{},0.1,None], 
     97            [{},(0.1,0.1),None], 
     98            [{},'ER',None], 
     99            [{},'VR',None], 
     100            ] 
     101        tests = smoke_tests + getattr(model_definition, 'tests', []) 
     102 
     103        if tests: # in case there are no smoke tests... 
    88104            #print '------' 
    89105            #print 'found tests in', model_name 
    90106            #print '------' 
    91      
    92             kernel = load_kernel(model) 
    93             suite.addTest(ModelTestCase(model_name, kernel, tests)) 
     107 
     108            # if ispy then use the dll loader to call pykernel 
     109            # don't try to call cl kernel since it will not be 
     110            # available in some environmentes. 
     111            ispy = callable(getattr(model_definition,'Iq', None)) 
     112 
     113            # test using opencl if desired 
     114            if not ispy and ('opencl' in loaders and load_model_cl): 
     115                test_name = "Model: %s, Kernel: OpenCL"%model_name 
     116                test = ModelTestCase(test_name, model_definition, 
     117                                     load_model_cl, tests) 
     118                #print "defining", test_name 
     119                suite.addTest(test) 
     120 
     121            # test using dll if desired 
     122            if ispy or ('dll' in loaders and load_model_dll): 
     123                test_name = "Model: %s, Kernel: dll"%model_name 
     124                test = ModelTestCase(test_name, model_definition, 
     125                                     load_model_dll, tests) 
     126                #print "defining", test_name 
     127                suite.addTest(test) 
    94128 
    95129    return suite 
    96130 
    97 class ModelTestCase(unittest.TestCase): 
    98      
    99     def __init__(self, model_name, kernel, tests): 
    100         unittest.TestCase.__init__(self) 
    101          
    102         self.model_name = model_name 
    103         self.kernel = kernel 
    104         self.tests = tests 
    105  
    106     def runTest(self): 
    107         #print '------' 
    108         #print self.model_name 
    109         #print '------' 
    110         try: 
    111             for test in self.tests: 
    112                 params = test[0] 
    113                 Q = test[1] 
    114                 I = test[2] 
    115                        
    116                 if not isinstance(Q, list): 
    117                     Q = [Q] 
    118                 if not isinstance(I, list): 
    119                     I = [I] 
    120                      
    121                 if isinstance(Q[0], tuple): 
    122                     npQ = [np.array([Qi[d] for Qi in Q]) for d in xrange(len(Q[0]))] 
     131def _hide_model_case_from_nosetests(): 
     132    class ModelTestCase(unittest.TestCase): 
     133        def __init__(self, test_name, definition, loader, tests): 
     134            unittest.TestCase.__init__(self) 
     135 
     136            self.test_name = test_name 
     137            self.definition = definition 
     138            self.loader = loader 
     139            self.tests = tests 
     140 
     141        def runTest(self): 
     142            try: 
     143                model = self.loader(self.definition) 
     144                for test in self.tests: 
     145                    self._run_one_test(model, test) 
     146 
     147            except Exception,exc: 
     148                annotate_exception(exc, self.test_name) 
     149                raise 
     150 
     151        def _run_one_test(self, model, test): 
     152            pars, x, y = test 
     153 
     154            if not isinstance(y, list): 
     155                y = [y] 
     156            if not isinstance(x, list): 
     157                x = [x] 
     158 
     159            self.assertEqual(len(y), len(x)) 
     160 
     161            if x[0] == 'ER': 
     162                actual = [call_ER(model.info, pars)] 
     163            elif x[0] == 'VR': 
     164                actual = [call_VR(model.info, pars)] 
     165            elif isinstance(x[0], tuple): 
     166                Qx,Qy = zip(*x) 
     167                q_vectors = [np.array(Qx), np.array(Qy)] 
     168                kernel = make_kernel(model, q_vectors) 
     169                actual = call_kernel(kernel, pars) 
     170            else: 
     171                q_vectors = [np.array(x)] 
     172                kernel = make_kernel(model, q_vectors) 
     173                actual = call_kernel(kernel, pars) 
     174 
     175            self.assertGreater(len(actual), 0) 
     176            self.assertEqual(len(y), len(actual)) 
     177 
     178            for xi, yi, actual_yi in zip(x, y, actual): 
     179                if yi is None: 
     180                    # smoke test --- make sure it runs and produces a value 
     181                    self.assertTrue(np.isfinite(actual_yi), 
     182                        'invalid f(%s): %s' % (xi, actual_yi)) 
    123183                else: 
    124                     npQ = [np.array(Q)] 
    125  
    126                 self.assertTrue(Q) 
    127                 self.assertEqual(len(I), len(Q))     
    128              
    129                 Iq = eval_kernel(self.kernel, npQ, params) 
    130              
    131                 self.assertGreater(len(Iq), 0)     
    132                 self.assertEqual(len(I), len(Iq))               
    133                  
    134                 for q, i, iq in zip(Q, I, Iq): 
    135                     err = np.abs(i - iq) 
    136                     nrm = np.abs(i) 
    137              
    138                     self.assertLess(err * 10**5, nrm, 'q:%s; expected:%s; actual:%s' % (q, i, iq)) 
    139                      
    140         except Exception,exc:  
    141             annotate_exception(exc, '\r\nModel: %s' % self.model_name) 
    142             raise 
     184                    err = abs(yi - actual_yi) 
     185                    nrm = abs(yi) 
     186                    self.assertLess(err * 10**5, nrm, 
     187                        'f(%s); expected:%s; actual:%s' % (xi, yi, actual_yi)) 
     188 
     189    return ModelTestCase 
     190 
     191 
     192# let nosetests sniff out the tests 
     193def model_tests(): 
     194    tests = make_suite(['opencl','dll'],['all']) 
     195    for test_i in tests: 
     196        yield test_i.runTest 
    143197 
    144198def main(): 
    145     #unittest.main() 
     199    models = sys.argv[1:] 
     200    if models and models[0] == 'opencl': 
     201        if load_model_cl is None: 
     202            print >>sys.stderr, "opencl is not available" 
     203            return 1 
     204        loaders = ['opencl'] 
     205        models = models[1:] 
     206    elif models and models[0] == 'dll': 
     207        # TODO: test if compiler is available? 
     208        loaders = ['dll'] 
     209        models = models[1:] 
     210    elif models and models[0] == 'opencl_and_dll': 
     211        if load_model_cl is None: 
     212            print >>sys.stderr, "opencl is not available" 
     213            return 1 
     214        loaders = ['opencl', 'dll'] 
     215        models = models[1:] 
     216    else: 
     217        loaders = ['opencl', 'dll'] 
     218    if not models: 
     219        print >>sys.stderr, "usage: python -m sasmodels.model_test [opencl|dll|opencl_and_dll] model1 model2 ..." 
     220        print >>sys.stderr, "if model1 is 'all', then all except the remaining models will be tested" 
     221        return 1 
     222 
     223    #run_tests(loaders, models) 
    146224    runner = unittest.TextTestRunner() 
    147     runner.run(suite()) 
     225    result = runner.run(make_suite(loaders, models)) 
     226    return 1 if result.failures or result.errors else 0 
    148227 
    149228if __name__ == "__main__": 
    150     main() 
     229    sys.exit(main()) 

  • sasmodels/models/barbell.py

    r0706431 r3c56da87  
    44r""" 
    55 
    6 Calculates the scattering from a barbell-shaped cylinder (This model simply becomes the DumBellModel when the length of 
    7 the cylinder, *L*, is set to zero). That is, a sphereocylinder with spherical end caps that have a radius larger than 
    8 that of the cylinder and the center of the end cap radius lies outside of the cylinder. All dimensions of the BarBell 
    9 are considered to be monodisperse. See the diagram for the details of the geometry and restrictions on parameter values. 
     6Calculates the scattering from a barbell-shaped cylinder (This model simply 
     7becomes the DumBellModel when the length of the cylinder, *L*, is set to zero). 
     8That is, a sphereocylinder with spherical end caps that have a radius larger 
     9than that of the cylinder and the center of the end cap radius lies outside 
     10of the cylinder. All dimensions of the BarBell are considered to be 
     11monodisperse. See the diagram for the details of the geometry and restrictions 
     12on parameter values. 
    1013 
    1114Definition 
    … …  
    1619The barbell geometry is defined as 
    1720 
    18 .. image:: img/image105.jpg 
     21.. image:: img/barbell_geometry.jpg 
    1922 
    20 where *r* is the radius of the cylinder. All other parameters are as defined in the diagram. 
     23where *r* is the radius of the cylinder. All other parameters are as defined 
     24in the diagram. 
    2125 
    2226Since the end cap radius 
    23 *R* >= *r* and by definition for this geometry *h* < 0, *h* is then defined by *r* and *R* as 
     27*R* >= *r* and by definition for this geometry *h* < 0, *h* is then 
     28defined by *r* and *R* as 
    2429 
    2530*h* = -1 \* sqrt(*R*\ :sup:`2` - *r*\ :sup:`2`) 
    … …  
    4651             \over QR\sin\theta \left(1-t^2\right)^{1/2}} 
    4752 
    48 The < > brackets denote an average of the structure over all orientations. <*A* :sup:`2`\ *(q)*> is then the form 
    49 factor, *P(q)*. The scale factor is equivalent to the volume fraction of cylinders, each of volume, *V*. Contrast is 
    50 the difference of scattering length densities of the cylinder and the surrounding solvent. 
     53The < > brackets denote an average of the structure over all orientations. 
     54<*A* :sup:`2`\ *(q)*> is then the form factor, *P(q)*. The scale factor is 
     55equivalent to the volume fraction of cylinders, each of volume, *V*. Contrast 
     56is the difference of scattering length densities of the cylinder and the 
     57surrounding solvent. 
    5158 
    5259The volume of the barbell is 
    … …  
    6976        \left( 4R^3 6R^2h - 2h^3 + 3r^2L \right)^{-1} 
    7077 
    71 **The requirement that** *R* >= *r* **is not enforced in the model!** It is up to you to restrict this during analysis. 
     78**The requirement that** *R* >= *r* **is not enforced in the model!** It is 
     79up to you to restrict this during analysis. 
    7280 
    73 This example dataset is produced by running the Macro PlotBarbell(), using 200 data points, *qmin* = 0.001 |Ang^-1|, 
    74 *qmax* = 0.7 |Ang^-1| and the following default values 
     81This example dataset is produced by running the Macro PlotBarbell(), 
     82using 200 data points, *qmin* = 0.001 |Ang^-1|, *qmax* = 0.7 |Ang^-1|, 
     83*sld* = 4e-6 |Ang^-2| and the default model values. 
    7584 
    76 ==============  ========  ============= 
    77 Parameter name  Units     Default value 
    78 ==============  ========  ============= 
    79 scale           None      1.0 
    80 len_bar         |Ang|     400.0 
    81 rad_bar         |Ang|     20.0 
    82 rad_bell        |Ang|     40.0 
    83 sld_barbell     |Ang^-2|  1.0e-006 
    84 sld_solv        |Ang^-2|  6.3e-006 
    85 background      |cm^-1|   0 
    86 ==============  ========  ============= 
    87  
    88 .. image:: img/image110.jpg 
     85.. image:: img/barbell_1d.jpg 
    8986 
    9087*Figure. 1D plot using the default values (w/256 data point).* 
    9188 
    92 For 2D data: The 2D scattering intensity is calculated similar to the 2D cylinder model. For example, for 
    93 |theta| = 45 deg and |phi| = 0 deg with default values for other parameters 
     89For 2D data: The 2D scattering intensity is calculated similar to the 2D 
     90cylinder model. For example, for |theta| = 45 deg and |phi| = 0 deg with 
     91default values for other parameters 
    9492 
    95 .. image:: img/image111.jpg 
     93.. image:: img/barbell_2d.jpg 
    9694 
    9795*Figure. 2D plot (w/(256X265) data points).* 
    … …  
    106104 
    107105REFERENCE 
     106--------- 
    108107 
    109108H Kaya, *J. Appl. Cryst.*, 37 (2004) 37 223-230 
    … …  
    112111 
    113112""" 
    114 from numpy import pi, inf 
     113from numpy import inf 
    115114 
    116115name = "barbell" 
    … …  
    122121                It must be that rad_bar <(=) rad_bell. 
    123122""" 
     123category = "shape:cylinder" 
    124124 
    125125parameters = [ 
    … …  
    135135 
    136136source = [ "lib/J1.c", "lib/gauss76.c", "barbell.c" ] 
    137  
    138 def ER(radius, length): 
    139     return 1.0 
    140137 
    141138# parameters for demo 

  • sasmodels/models/bcc.c

    rc95dc908 r6137124  
    112112 
    113113  const double Da = d_factor*dnn; 
    114   const double qDa_2 = q*q*Da*Da; 
    115114  const double s1 = dnn/sqrt(0.75); 
    116115 

  • sasmodels/models/bcc.py

    re166cb9 r3c56da87  
    22#note model title and parameter table are automatically inserted 
    33#note - calculation requires double precision 
    4 """ 
    5 Calculates the scattering from a **body-centered cubic lattice** with paracrystalline distortion. Thermal vibrations 
    6 are considered to be negligible, and the size of the paracrystal is infinitely large. Paracrystalline distortion is 
    7 assumed to be isotropic and characterized by a Gaussian distribution. 
     4r""" 
     5Calculates the scattering from a **body-centered cubic lattice** with 
     6paracrystalline distortion. Thermal vibrations are considered to be negligible, 
     7and the size of the paracrystal is infinitely large. Paracrystalline distortion 
     8is assumed to be isotropic and characterized by a Gaussian distribution. 
    89 
    910The returned value is scaled to units of |cm^-1|\ |sr^-1|, absolute scale. 
    … …  
    1213---------- 
    1314 
    14 The scattering intensity *I(q)* is calculated as 
     15The scattering intensity $I(q)$ is calculated as 
    1516 
    16 .. image:: img/image167.jpg 
     17.. math: 
    1718 
    18 where *scale* is the volume fraction of spheres, *Vp* is the volume of the primary particle, *V(lattice)* is a volume 
    19 correction for the crystal structure, *P(q)* is the form factor of the sphere (normalized), and *Z(q)* is the 
    20 paracrystalline structure factor for a body-centered cubic structure. 
     19    I(q) = \frac{\text{scale}}{V_P} V_\text{lattice} P(q) Z(q) 
    2120 
    22 Equation (1) of the 1990 reference is used to calculate *Z(q)*, using equations (29)-(31) from the 1987 paper for 
    23 *Z1*\ , *Z2*\ , and *Z3*\ . 
    2421 
    25 The lattice correction (the occupied volume of the lattice) for a body-centered cubic structure of particles of radius 
    26 *R* and nearest neighbor separation *D* is 
     22where *scale* is the volume fraction of spheres, *Vp* is the volume of the 
     23primary particle, *V(lattice)* is a volume correction for the crystal 
     24structure, $P(q)$ is the form factor of the sphere (normalized), and $Z(q)$ 
     25is the paracrystalline structure factor for a body-centered cubic structure. 
    2726 
    28 .. image:: img/image159.jpg 
     27Equation (1) of the 1990 reference is used to calculate $Z(q)$, using 
     28equations (29)-(31) from the 1987 paper for *Z1*\ , *Z2*\ , and *Z3*\ . 
    2929 
    30 The distortion factor (one standard deviation) of the paracrystal is included in the calculation of *Z(q)* 
     30The lattice correction (the occupied volume of the lattice) for a 
     31body-centered cubic structure of particles of radius $R$ and nearest neighbor 
     32separation $D$ is 
    3133 
    32 .. image:: img/image160.jpg 
     34.. math: 
    3335 
    34 where *g* is a fractional distortion based on the nearest neighbor distance. 
     36    V_\text{lattice} = \frac{16\pi}{3} \frac{R^3}{\left(D\sqrt{2}\right)^3} 
     37 
     38 
     39The distortion factor (one standard deviation) of the paracrystal is included 
     40in the calculation of $Z(q)$ 
     41 
     42.. math: 
     43 
     44    \Delta a = g D 
     45 
     46where $g$ is a fractional distortion based on the nearest neighbor distance. 
    3547 
    3648The body-centered cubic lattice is 
    3749 
    38 .. image:: img/image168.jpg 
     50.. image:: img/bcc_lattice.jpg 
    3951 
    4052For a crystal, diffraction peaks appear at reduced q-values given by 
    4153 
    42 .. image:: img/image162.jpg 
     54.. math: 
    4355 
    44 where for a body-centered cubic lattice, only reflections where (\ *h* + *k* + *l*\ ) = even are allowed and 
    45 reflections where (\ *h* + *k* + *l*\ ) = odd are forbidden. Thus the peak positions correspond to (just the first 5) 
     56    \frac{qD}{2\pi} = \sqrt{h^2 + k^2 + l^2} 
    4657 
    47 .. image:: img/image169.jpg 
     58where for a body-centered cubic lattice, only reflections where 
     59$(h + k + l) = \text{even}$ are allowed and reflections where 
     60$(h + k + l) = \text{odd}$ are forbidden. Thus the peak positions 
     61correspond to (just the first 5) 
    4862 
    49 **NB: The calculation of** *Z(q)* **is a double numerical integral that must be carried out with a high density of** 
    50 **points to properly capture the sharp peaks of the paracrystalline scattering.** So be warned that the calculation is 
    51 SLOW. Go get some coffee. Fitting of any experimental data must be resolution smeared for any meaningful fit. This 
    52 makes a triple integral. Very, very slow. Go get lunch! 
     63.. math: 
    5364 
    54 This example dataset is produced using 200 data points, *qmin* = 0.001 |Ang^-1|, *qmax* = 0.1 |Ang^-1| and the above 
    55 default values. 
     65    \begin{eqnarray} 
     66    &q/q_o&&\quad 1&& \ \sqrt{2} && \ \sqrt{3} && \ \sqrt{4} && \ \sqrt{5} \\ 
     67    &\text{Indices}&& (110) && (200) && (211) && (220) && (310) 
     68    \end{eqnarray} 
    5669 
    57 .. image:: img/image170.jpg 
     70**NB: The calculation of $Z(q)$ is a double numerical integral that must 
     71be carried out with a high density of points to properly capture the sharp 
     72peaks of the paracrystalline scattering.** So be warned that the calculation 
     73is SLOW. Go get some coffee. Fitting of any experimental data must be 
     74resolution smeared for any meaningful fit. This makes a triple integral. 
     75Very, very slow. Go get lunch! 
    5876 
    59 *Figure. 1D plot in the linear scale using the default values (w/200 data point).* 
     77This example dataset is produced using 200 data points, 
     78*qmin* = 0.001 |Ang^-1|, *qmax* = 0.1 |Ang^-1| and the above default values. 
    6079 
    61 The 2D (Anisotropic model) is based on the reference below where *I(q)* is approximated for 1d scattering. Thus the 
    62 scattering pattern for 2D may not be accurate. Note that we are not responsible for any incorrectness of the 2D model 
    63 computation. 
     80.. image:: img/bcc_1d.jpg 
    6481 
    65 .. image:: img/image165.gif 
     82*Figure. 1D plot in the linear scale using the default values 
     83(w/200 data point).* 
    6684 
    67 .. image:: img/image171.jpg 
     85The 2D (Anisotropic model) is based on the reference below where $I(q)$ is 
     86approximated for 1d scattering. Thus the scattering pattern for 2D may not 
     87be accurate. Note that we are not responsible for any incorrectness of the 2D 
     88model computation. 
     89 
     90.. image:: img/bcc_orientation.gif 
     91 
     92.. image:: img/bcc_2d.jpg 
    6893 
    6994*Figure. 2D plot using the default values (w/200X200 pixels).* 
    7095 
    7196REFERENCE 
     97--------- 
    7298 
    7399Hideki Matsuoka et. al. *Physical Review B*, 36 (1987) 1754-1765 
    … …  
    78104""" 
    79105 
    80 from numpy import pi, inf 
     106from numpy import inf 
    81107 
    82108name = "bcc_paracrystal" 
    … …  
    87113    assumed to be isotropic and characterized by a Gaussian distribution. 
    88114    """ 
     115category="shape:paracrystal" 
    89116 
    90117parameters = [ 
    … …  
    94121    [ "radius", "Ang",  40, [0, inf], "volume","Particle radius" ], 
    95122    [ "sld", "1e-6/Ang^2", 4, [-inf,inf], "", "Particle scattering length density" ], 
    96     [ "solvent_sld", "1e-6/Ang^2", 1, [-inf,inf], "","Solvent scattering length density" ], 
     123    [ "solvent_sld", "1e-6/Ang^2", 1, [-inf,inf], "", "Solvent scattering length density" ], 
    97124    [ "theta", "degrees", 60, [-inf, inf], "orientation","In plane angle" ], 
    98125    [ "phi", "degrees", 60, [-inf, inf], "orientation","Out of plane angle" ], 
    … …  
    101128 
    102129source = [ "lib/J1.c", "lib/gauss150.c", "bcc.c" ] 
    103  
    104 def ER(radius, length): 
    105     return 0 
    106130 
    107131# parameters for demo 

  • sasmodels/models/broad_peak.py

    rf57d123 r3c56da87  
    66layered structures, etc. 
    77 
    8 The d-spacing corresponding to the broad peak is a characteristic distance  
    9 between the scattering inhomogeneities (such as in lamellar, cylindrical, or  
     8The d-spacing corresponding to the broad peak is a characteristic distance 
     9between the scattering inhomogeneities (such as in lamellar, cylindrical, or 
    1010spherical morphologies, or for bicontinuous structures). 
    1111 
    … …  
    3030    q = \sqrt{q_x^2 + q_y^2} 
    3131 
    32 =====================  =========  ============= 
    33 Parameter name          Units     Default value 
    34 =====================  =========  ============= 
    35 lorentz_scale(=C)       None      10 
    36 porod_scale  (=A)       None      1e-05 
    37 lorentz_length (= |xi| )  |Ang|     50 
    38 peak_pos  (=Q0)         |Ang^-1|  0.1 
    39 porod_exp  (=n)         None      2 
    40 lorentz_exp (=m)        None      3 
    41 Background (=B)        |cm^-1|   0.1 
    42 ==================  ========  ============= 
    4332 
    4433.. image:: img/image175.jpg 
    … …  
    4736 
    4837REFERENCE 
     38--------- 
    4939 
    5040None. 
    … …  
    5444""" 
    5545 
    56 import numpy as np 
    57 from numpy import pi, inf, sin, cos, sqrt, exp, log, fabs 
     46from numpy import inf, sqrt 
    5847 
    5948name = "broad_peak" 
    … …  
    7160      lorentz_exp = Lorentzian exponent 
    7261      background = Incoherent background""" 
     62category="shape-independent" 
    7363 
    7464parameters = [ 
    … …  
    9181 
    9282 
    93 def form_volume(): 
    94     return 1 
     83def Iq(q, porod_scale, porod_exp, lorentz_scale, lorentz_length, peak_pos, lorentz_exp): 
     84    inten = (porod_scale/q**porod_exp + lorentz_scale 
     85        / (1.0 + (abs(q-peak_pos)*lorentz_length)**lorentz_exp)) 
     86    return inten 
     87Iq.vectorized = True  # Iq accepts an array of Q values 
    9588 
    96 def Iq(q, porod_scale, porod_exp, lorentz_scale, lorentz_length, peak_pos, lorentz_exp): 
    97     inten = porod_scale/pow(q,porod_exp) + lorentz_scale/(1.0 \ 
    98         + pow((fabs(q-peak_pos)*lorentz_length),lorentz_exp)) 
    99     return inten   
    100  
    101 # FOR VECTORIZED VERSION, UNCOMMENT THE NEXT LINE 
    102 Iq.vectorized = True 
    103  
    104 def Iqxy(qx, qy, porod_scale, porod_exp, lorentz_scale, lorentz_length, peak_pos, lorentz_exp): 
    105     return Iq(sqrt(qx**2 + qy**2), porod_scale, porod_exp, lorentz_scale, lorentz_length, peak_pos, lorentz_exp) 
    106  
    107 # FOR VECTORIZED VERSION, UNCOMMENT THE NEXT LINE 
    108 Iqxy.vectorized = True 
     89def Iqxy(qx, qy, *args): 
     90    return Iq(sqrt(qx**2 + qy**2), *args) 
     91Iqxy.vectorized = True # Iqxy accepts an array of Qx, Qy values 
    10992 
    11093 
    … …  
    11497    lorentz_scale=10,lorentz_length=50, peak_pos=0.1, lorentz_exp=2, 
    11598    ) 
     99 
    116100oldname = "BroadPeakModel" 
    117 oldpars = dict(porod_scale='scale_p', porod_exp='exponent_p',  
    118         lorentz_scale='scale_l', lorentz_length='length_l', peak_pos='q_peak',  
     101oldpars = dict(porod_scale='scale_p', porod_exp='exponent_p', 
     102        lorentz_scale='scale_l', lorentz_length='length_l', peak_pos='q_peak', 
    119103        lorentz_exp='exponent_l', scale=None) 

  • sasmodels/models/capped_cylinder.py

    ra503bfd ra5d0d00  
    130130        solvent_sld: SLD of the solvent. 
    131131""" 
     132category = "shape:cylinder" 
    132133 
    133134parameters = [ 

  • sasmodels/models/core_shell_cylinder.py

    ra503bfd ra5d0d00  
    133133        phi: the axis_phi of the cylinder 
    134134""" 
     135category = "shape:cylinder" 
    135136 
    136137parameters = [ 

  • sasmodels/models/cylinder.py

    r143e2f7 r9890053  
    128128            f(q,alpha)^(2)*sin(alpha)*dalpha + background 
    129129""" 
     130category = "shape:cylinder" 
    130131 
    131132parameters = [ 
    … …  
    157158    sld=6, solvent_sld=1, 
    158159    #radius=5, length=20, 
    159     radius=260, length=290, 
    160     theta=30, phi=0, 
     160    radius=20, length=300, 
     161    theta=60, phi=60, 
    161162    radius_pd=.2, radius_pd_n=9, 
    162163    length_pd=.2,length_pd_n=10, 
    163     theta_pd=15, theta_pd_n=45, 
    164     phi_pd=15, phi_pd_n=1, 
     164    theta_pd=10, theta_pd_n=5, 
     165    phi_pd=10, phi_pd_n=5, 
    165166    ) 
    166167 
    … …  
    171172 
    172173 
    173 # test values: 
    174 # [ 
    175 #   [ {parameters}, q, I(q)], 
    176 #   [ {parameters}, [q], [I(q)] ], 
    177 #   [ {parameters}, [q1, q2, ...], [I(q1), I(q2), ...]], 
    178  
    179 #   [ {parameters}, (qx, qy), I(qx, Iqy)], 
    180 #   [ {parameters}, [(qx1, qy1), (qx2, qy2), ...], [I(qx1,qy1), I(qx2,qy2), ...], 
    181 #   ... 
    182 # ] 
    183 # Precision defaults to 7 digits (relative) for single, 14 for double 
    184 # May want a limited precision version that checks for 8-n or 15-n digits respectively 
    185174qx,qy = 0.2*np.cos(2.5), 0.2*np.sin(2.5) 
    186175tests = [ 

  • sasmodels/models/ellipsoid.py

    ra503bfd r3c56da87  
    116116""" 
    117117 
    118 from numpy import pi, inf 
     118from numpy import inf 
    119119 
    120120name = "ellipsoid" 
    … …  
    136136                Re: equatorial radius of the ellipsoid 
    137137""" 
     138category = "shape:ellipsoid" 
    138139 
    139140parameters = [ 

  • sasmodels/models/fcc.c

    re7b3d7b r6137124  
    112112 
    113113  const double Da = d_factor*dnn; 
    114   const double qDa_2 = q*q*Da*Da; 
    115114  const double s1 = dnn/sqrt(0.75); 
    116115 

  • sasmodels/models/fcc.py

    re7b3d7b r3c56da87  
    33#note - calculation requires double precision 
    44r""" 
    5 Calculates the scattering from a **face-centered cubic lattice** with paracrystalline distortion. Thermal vibrations 
    6 are considered to be negligible, and the size of the paracrystal is infinitely large. Paracrystalline distortion is 
    7 assumed to be isotropic and characterized by a Gaussian distribution. 
     5Calculates the scattering from a **face-centered cubic lattice** with 
     6paracrystalline distortion. Thermal vibrations are considered to be 
     7negligible, and the size of the paracrystal is infinitely large. 
     8Paracrystalline distortion is assumed to be isotropic and characterized by 
     9a Gaussian distribution. 
    810 
    911The returned value is scaled to units of |cm^-1|\ |sr^-1|, absolute scale. 
    … …  
    1618.. image:: img/image158.jpg 
    1719 
    18 where *scale* is the volume fraction of spheres, *Vp* is the volume of the primary particle, *V(lattice)* is a volume 
    19 correction for the crystal structure, *P(q)* is the form factor of the sphere (normalized), and *Z(q)* is the 
    20 paracrystalline structure factor for a face-centered cubic structure. 
     20where *scale* is the volume fraction of spheres, *Vp* is the volume of 
     21the primary particle, *V(lattice)* is a volume correction for the crystal 
     22structure, *P(q)* is the form factor of the sphere (normalized), and *Z(q)* 
     23is the paracrystalline structure factor for a face-centered cubic structure. 
    2124 
    22 Equation (1) of the 1990 reference is used to calculate *Z(q)*, using equations (23)-(25) from the 1987 paper for 
    23 *Z1*\ , *Z2*\ , and *Z3*\ . 
     25Equation (1) of the 1990 reference is used to calculate *Z(q)*, using 
     26equations (23)-(25) from the 1987 paper for *Z1*\ , *Z2*\ , and *Z3*\ . 
    2427 
    25 The lattice correction (the occupied volume of the lattice) for a face-centered cubic structure of particles of radius 
    26 *R* and nearest neighbor separation *D* is 
     28The lattice correction (the occupied volume of the lattice) for a 
     29face-centered cubic structure of particles of radius *R* and nearest 
     30neighbor separation *D* is 
    2731 
    2832.. image:: img/image159.jpg 
    2933 
    30 The distortion factor (one standard deviation) of the paracrystal is included in the calculation of *Z(q)* 
     34The distortion factor (one standard deviation) of the paracrystal is 
     35included in the calculation of *Z(q)* 
    3136 
    3237.. image:: img/image160.jpg 
    … …  
    4247.. image:: img/image162.jpg 
    4348 
    44 where for a face-centered cubic lattice *h*\ , *k*\ , *l* all odd or all even are allowed and reflections where 
    45 *h*\ , *k*\ , *l* are mixed odd/even are forbidden. Thus the peak positions correspond to (just the first 5) 
     49where for a face-centered cubic lattice *h*\ , *k*\ , *l* all odd or all 
     50even are allowed and reflections where *h*\ , *k*\ , *l* are mixed odd/even 
     51are forbidden. Thus the peak positions correspond to (just the first 5) 
    4652 
    4753.. image:: img/image163.jpg 
    4854 
    49 **NB: The calculation of** *Z(q)* **is a double numerical integral that must be carried out with a high density of** 
    50 **points to properly capture the sharp peaks of the paracrystalline scattering.** So be warned that the calculation is 
    51 SLOW. Go get some coffee. Fitting of any experimental data must be resolution smeared for any meaningful fit. This 
    52 makes a triple integral. Very, very slow. Go get lunch! 
     55**NB: The calculation of** *Z(q)* **is a double numerical integral that 
     56must be carried out with a high density of** **points to properly capture 
     57the sharp peaks of the paracrystalline scattering.** So be warned that the 
     58calculation is SLOW. Go get some coffee. Fitting of any experimental data 
     59must be resolution smeared for any meaningful fit. This makes a triple 
     60integral. Very, very slow. Go get lunch! 
    5361 
    54 This example dataset is produced using 200 data points, *qmin* = 0.01 |Ang^-1|, *qmax* = 0.1 |Ang^-1| and the above 
    55 default values. 
     62This example dataset is produced using 200 data points, *qmin* = 0.01 |Ang^-1|, 
     63*qmax* = 0.1 |Ang^-1| and the above default values. 
    5664 
    5765.. image:: img/image164.jpg 
    … …  
    5967*Figure. 1D plot in the linear scale using the default values (w/200 data point).* 
    6068 
    61 The 2D (Anisotropic model) is based on the reference below where *I(q)* is approximated for 1d scattering. Thus the 
    62 scattering pattern for 2D may not be accurate. Note that we are not responsible for any incorrectness of the 2D model 
    63 computation. 
     69The 2D (Anisotropic model) is based on the reference below where *I(q)* is 
     70approximated for 1d scattering. Thus the scattering pattern for 2D may not 
     71be accurate. Note that we are not responsible for any incorrectness of the 
     722D model computation. 
    6473 
    6574.. image:: img/image165.gif 
    … …  
    7887""" 
    7988 
    80 from numpy import pi, inf 
     89from numpy import inf 
    8190 
    8291name = "fcc_paracrystal" 
    … …  
    8796    assumed to be isotropic and characterized by a Gaussian distribution. 
    8897    """ 
     98category = "shape:paracrystal" 
    8999 
    90100parameters = [ 
    … …  
    101111 
    102112source = [ "lib/J1.c", "lib/gauss150.c", "fcc.c" ] 
    103  
    104 def ER(radius, length): 
    105     return 0 
    106113 
    107114# parameters for demo 
    … …  
    120127# names and the target sasview model name. 
    121128oldname='FCCrystalModel' 
    122 oldpars=dict(sld='sldSph', 
    123              solvent_sld='sldSolv') 
     129oldpars=dict(sld='sldSph', solvent_sld='sldSolv') 

  • sasmodels/models/hardsphere.py

    r301e096 r3c56da87  
    11# Note: model title and parameter table are inserted automatically 
    2 r"""This calculates the interparticle structure factor for monodisperse spherical particles interacting through hard 
    3 sphere (excluded volume) interactions. 
     2r"""Calculate the interparticle structure factor for monodisperse 
     3spherical particles interacting through hard sphere (excluded volume) 
     4interactions. 
    45 
    5 The calculation uses the Percus-Yevick closure where the interparticle potential is 
     6The calculation uses the Percus-Yevick closure where the interparticle 
     7potential is 
    68 
    7 .. image:: img/HardSphere_223.PNG 
     9.. math: 
     10 
     11    U(r) = \begin{cases} 
     12    \infty & r < 2R \\ 
     13    0 & r \geq 2R 
     14    \end{cases} 
    815 
    916where *r* is the distance from the center of the sphere of a radius *R*. 
    … …  
    1320.. math:: 
    1421 
    15     Q = \sqrt{Q_x^2 + Q_y^2} 
     22    q = \sqrt{q_x^2 + q_y^2} 
    1623 
    1724 
    18 ==============  ========  ============= 
    19 Parameter name  Units     Default value 
    20 ==============  ========  ============= 
    21 effect_radius   |Ang|     50.0 
    22 volfraction     None      0.2 
    23 ==============  ========  ============= 
    24  
    25 .. image:: img/HardSphere_224.jpg 
     25.. image:: img/HardSphere_1d.jpg 
    2626 
    2727*Figure. 1D plot using the default values (in linear scale).* 
    … …  
    3232""" 
    3333 
    34 from numpy import pi, inf 
     34from numpy import inf 
    3535 
    3636name = "hardsphere" 
    … …  
    3939        [Hard sphere structure factor, with Percus-Yevick closure] 
    4040        Interparticle S(Q) for random, non-interacting spheres. 
    41                 May be a reasonable approximation for other shapes of 
    42                 particles that freely rotate, and for moderately polydisperse 
     41        May be a reasonable approximation for other shapes of 
     42        particles that freely rotate, and for moderately polydisperse 
    4343        systems. Though strictly the maths needs to be modified -  
    44                 which sasview does not do yet. 
    45                 effect_radius is the hard sphere radius 
    46                 volfraction is the volume fraction occupied by the spheres. 
     44        which sasview does not do yet. 
     45        effect_radius is the hard sphere radius 
     46        volfraction is the volume fraction occupied by the spheres. 
    4747""" 
     48category = "structure-factor" 
    4849 
    4950parameters = [ 
    … …  
    9697Iqxy = """ 
    9798    // never called since no orientation or magnetic parameters. 
    98     return -1.0; 
     99    return Iq(sqrt(qx*qx+qy*qy), IQ_PARAMETERS); 
    99100    """ 
    100101 

  • sasmodels/models/lamellar.py

    ra503bfd r3c56da87  
    4747""" 
    4848 
    49 from numpy import pi, inf 
     49from numpy import inf 
    5050 
    5151name = "lamellar" 
    … …  
    6161                scale = scale factor 
    6262""" 
     63category = "shape:lamellae" 
    6364 
    6465parameters = [ 
    … …  
    8889 
    8990Iqxy = """ 
    90     // never called since no orientation or magnetic parameters. 
    91     return -1.0; 
     91    return Iq(sqrt(qx*qx+qy*qy), IQ_PARAMETERS); 
    9292    """ 
    9393 

  • sasmodels/models/lamellarCaille.py

    rdc02af0 r3c56da87  
    11# Note: model title and parameter table are inserted automatically 
    22r""" 
    3 This model provides the scattering intensity, *I(q)* = *P(q)* \* *S(q)*, for a lamellar phase where a random 
    4 distribution in solution are assumed. Here a Caille S(Q) is used for the lamellar stacks. 
     3This model provides the scattering intensity, $I(q) = P(q) S(q)$, for a 
     4lamellar phase where a random distribution in solution are assumed. 
     5Here a Caille $S(Q)$ is used for the lamellar stacks. 
    56 
    6 The scattering intensity *I(q)* is 
     7The scattering intensity $I(q)$ is 
    78 
    8 .. image:: img/lamellarCaille_139.PNG 
     9.. math: 
     10 
     11    I(q) = 2\pi \frac{P(q)S(q)}{\delta q^2} 
    912 
    1013The form factor is 
    1114 
    12 .. image:: img/lamellarCaille_134.PNG 
     15.. math: 
     16 
     17    P(q) = \frac{2\Delta\rho^2}{q^2}\left(1-\cos q\delta \right) 
    1318 
    1419and the structure factor is 
    1520 
    16 .. image:: img/lamellarCaille_.PNG 
     21.. math: 
     22 
     23    S(q) = 1 + 2 \sum_1^{N-1}\left(1-\frac{n}{N}\right) 
     24           \cos(qdn)\exp\left(-\frac{2q^2d^2\alpha(n)}{2}\right) 
    1725 
    1826where 
    1927 
    20 .. image:: img/lamellarCaille_.PNG 
     28.. math: 
    2129 
    22 Here *d* = (repeat) spacing, |delta| = bilayer thickness, the contrast |drho| = SLD(headgroup) - SLD(solvent), 
    23 K = smectic bending elasticity, B = compression modulus, and N = number of lamellar plates (*n_plates*). 
     30    \begin{eqnarray} 
     31    \alpha(n) &=& \frac{\eta_{cp}}{4\pi^2} \left(\ln(\pi n)+\gamma_E\right)  \\ 
     32    \gamma_E &=& 0.5772156649 && \text{Euler's constant} \\ 
     33    \eta_{cp} &=& \frac{q_o^2k_B T}{8\pi\sqrt{K\overline{B}}} && \text{Caille constant} 
     34    \end{eqnarray} 
    2435 
    25 NB: **When the Caille parameter is greater than approximately 0.8 to 1.0, the assumptions of the model are incorrect.** 
    26 And due to a complication of the model function, users are responsible for making sure that all the assumptions are 
    27 handled accurately (see the original reference below for more details). 
     36Here $d$ = (repeat) spacing, $\delta$ = bilayer thickness, 
     37the contrast $\Delta\rho$ = SLD(headgroup) - SLD(solvent), 
     38$K$ = smectic bending elasticity, $B$ = compression modulus, and 
     39$N$ = number of lamellar plates (*n_plates*). 
    2840 
    29 Non-integer numbers of stacks are calculated as a linear combination of results for the next lower and higher values. 
     41NB: **When the Caille parameter is greater than approximately 0.8 to 1.0, the 
     42assumptions of the model are incorrect.** And due to a complication of the 
     43model function, users are responsible for making sure that all the assumptions 
     44are handled accurately (see the original reference below for more details). 
    3045 
    31 The 2D scattering intensity is calculated in the same way as 1D, where the *q* vector is defined as 
     46Non-integer numbers of stacks are calculated as a linear combination of 
     47results for the next lower and higher values. 
     48 
     49The 2D scattering intensity is calculated in the same way as 1D, where the 
     50$q$ vector is defined as 
    3251 
    3352.. math:: 
    3453 
    35     Q = \sqrt{Q_x^2 + Q_y^2} 
     54    q = \sqrt{q_x^2 + q_y^2} 
    3655 
    3756The returned value is in units of |cm^-1|, on absolute scale. 
    3857 
    39 ==============  ========  ============= 
    40 Parameter name  Units     Default value 
    41 ==============  ========  ============= 
    42 background      |cm^-1|   0.0 
    43 contrast        |Ang^-2|  5e-06 
    44 scale           None      1 
    45 delta           |Ang|     30 
    46 n_plates        None      20 
    47 spacing         |Ang|     400 
    48 caille          |Ang^-2|  0.1 
    49 ==============  ========  ============= 
    50  
    51 .. image:: img/lamellarPS_142.jpg 
     58.. image:: img/lamellarCaille_1d.jpg 
    5259 
    5360*Figure. 1D plot using the default values (w/6000 data point).* 
    5461 
    55 Our model uses the form factor calculations implemented in a c-library provided by the NIST Center for Neutron Research 
    56 (Kline, 2006). 
     62Our model uses the form factor calculations as implemented in a c library 
     63provided by the NIST Center for Neutron Research (Kline, 2006). 
    5764 
    5865REFERENCE 
     66--------- 
    5967 
    6068F Nallet, R Laversanne, and D Roux, J. Phys. II France, 3, (1993) 487-502 
    … …  
    6270also in J. Phys. Chem. B, 105, (2001) 11081-11088 
    6371""" 
    64 from numpy import pi, inf 
     72from numpy import inf 
    6573 
    6674name = "lamellarPS" 
    … …  
    7583                scale = scale factor 
    7684""" 
     85category = "shape:lamellae" 
    7786 
    7887parameters = [ 
    … …  
    8594    [ "spacing", "Ang", 400., [0.0,inf], "volume", 
    8695      "d-spacing of Caille S(Q)" ], 
    87     [ "Caille_parameter", "Ang^-2", 0.1, [0.0,0.8], "", 
     96    [ "Caille_parameter", "1/Ang^2", 0.1, [0.0,0.8], "", 
    8897      "Caille parameter" ], 
    8998    [ "sld", "1e-6/Ang^2", 6.3, [-inf,inf], "", 
    … …  
    102111 
    103112Iqxy = """ 
    104     // never called since no orientation or magnetic parameters. 
    105     return -1.0; 
     113    return Iq(sqrt(qx*qx+qy*qy), IQ_PARAMETERS); 
    106114    """ 
    107115 

  • sasmodels/models/lamellarCailleHG.py

    rcd55ac3 r3c56da87  
    11# Note: model title and parameter table are inserted automatically 
    22r""" 
    3 This model provides the scattering intensity, *I(q)* = *P(q)* \* *S(q)*, for a lamellar phase where a random 
    4 distribution in solution are assumed. Here a Caille S(Q) is used for the lamellar stacks. 
     3This model provides the scattering intensity, $I(q) = P(q)S(q)$, for a lamellar 
     4phase where a random distribution in solution are assumed. Here a Caille $S(Q)$ 
     5is used for the lamellar stacks. 
    56 
    6 The scattering intensity *I(q)* is 
     7The scattering intensity $I(q)$ is 
    78 
    8 .. image:: img/lamellarCailleHG_139.PNG 
     9.. math:: 
    910 
    10 The form factor is 
     11    I(q) = 2 \pi \frac{P(q)S(q)}{\delta q^2} 
    1112 
    12 .. image:: img/lamellarCailleHG_143.PNG 
    1313 
    14 and the structure factor is 
     14The form factor $P(q)$ is 
    1515 
    16 .. image:: img/lamellarCailleHG_140.PNG 
     16.. math:: 
     17 
     18        P(q) = \frac{4}{q^2}\big\{ 
     19        \Delta\rho_H \left[\sin[q(\delta_H + \delta_T)] - \sin(q\delta_T)\right] 
     20            + \Delta\rho_T\sin(q\delta_T)\big\}^2 
     21 
     22and the structure factor $S(q)$ is 
     23 
     24.. math:: 
     25 
     26    S(q) = 1 + 2 \sum_1^{N-1}\left(1-\frac{n}{N}\right) 
     27        \cos(qdn)\exp\left(-\frac{2q^2d^2\alpha(n)}{2}\right) 
    1728 
    1829where 
    1930 
    20 .. image:: img/lamellarCailleHG_141.PNG 
     31.. math:: 
    2132 
    22 where |delta|\ T = tail length (or *tail_length*), |delta|\ H = head thickness (or *h_thickness*), 
    23 |drho|\ H = SLD(headgroup) - SLD(solvent), and |drho|\ T = SLD(tail) - SLD(headgroup). 
    24 Here *d* = (repeat) spacing, *K* = smectic bending elasticity, *B* = compression modulus, and N = number of lamellar 
    25 plates (*n_plates*). 
     33    \begin{eqnarray} 
     34    \alpha(n) &=& \frac{\eta_{cp}}{4\pi^2} \left(\ln(\pi n)+\gamma_E\right)  \\ 
     35    \gamma_E &=& 0.5772156649&&\text{Euler's constant} \\ 
     36    \eta_{cp} &=& \frac{q_o^2k_B T}{8\pi\sqrt{K\overline{B}}} && \text{Caille constant} 
     37    \end{eqnarray} 
    2638 
    27 NB: **When the Caille parameter is greater than approximately 0.8 to 1.0, the assumptions of the model are incorrect.** 
    28 And due to a complication of the model function, users are responsible for making sure that all the assumptions are 
    29 handled accurately (see the original reference below for more details). 
    3039 
    31 Non-integer numbers of stacks are calculated as a linear combination of results for the next lower and higher values. 
     40$\delta_T$ is the tail length (or *tail_length*), $\delta_H$ is the head 
     41thickness (or *head_length*), $\Delta\rho_H$ is SLD(headgroup) - SLD(solvent), 
     42and $\Delta\rho_T$ is SLD(tail) - SLD(headgroup). Here $d$ is (repeat) spacing, 
     43$K$ is smectic bending elasticity, $B$ is compression modulus, and $N$ is the 
     44number of lamellar plates (*Nlayers*). 
    3245 
    33 The 2D scattering intensity is calculated in the same way as 1D, where the *q* vector is defined as 
     46NB: **When the Caille parameter is greater than approximately 0.8 to 1.0, the 
     47assumptions of the model are incorrect.**  And due to a complication of the 
     48model function, users are responsible for making sure that all the assumptions 
     49are handled accurately (see the original reference below for more details). 
     50 
     51Non-integer numbers of stacks are calculated as a linear combination of 
     52results for the next lower and higher values. 
     53 
     54The 2D scattering intensity is calculated in the same way as 1D, where 
     55the $q$ vector is defined as 
    3456 
    3557.. math:: 
    3658 
    37     Q = \sqrt{Q_x^2 + Q_y^2} 
     59    q = \sqrt{q_x^2 + q_y^2} 
    3860 
    3961The returned value is in units of |cm^-1|, on absolute scale. 
    4062 
    41 ==============  ========  ============= 
    42 Parameter name  Units     Default value 
    43 ==============  ========  ============= 
    44 background      |cm^-1|   0.001 
    45 sld_head        |Ang^-2|  2e-06 
    46 scale           None      1 
    47 sld_solvent     |Ang^-2|  6e-06 
    48 deltaH          |Ang|     2 
    49 deltaT          |Ang|     10 
    50 sld_tail        |Ang^-2|  0 
    51 n_plates        None      30 
    52 spacing         |Ang|     40 
    53 caille          |Ang^-2|  0.001 
    54 ==============  ========  ============= 
    55  
    56 .. image:: img/lamellarCailleHG_142.jpg 
     63.. image:: img/lamellarCailleHG_1d.jpg 
    5764 
    5865*Figure. 1D plot using the default values (w/6000 data point).* 
    5966 
    60 Our model uses the form factor calculations implemented in a c-library provided by the NIST Center for Neutron Research 
    61 (Kline, 2006). 
     67Our model uses the form factor calculations implemented in a C library provided 
     68by the NIST Center for Neutron Research (Kline, 2006). 
    6269 
    6370REFERENCE 
    … …  
    6774also in J. Phys. Chem. B, 105, (2001) 11081-11088 
    6875""" 
    69 from numpy import pi, inf 
     76from numpy import inf 
    7077 
    7178name = "lamellarCailleHG" 
    … …  
    8289                scale = scale factor 
    8390""" 
     91category = "shape:lamellae" 
    8492 
    8593parameters = [ 
    … …  
    113121 
    114122Iqxy = """ 
    115     // never called since no orientation or magnetic parameters. 
    116     return -1.0; 
     123    return Iq(sqrt(qx*qx+qy*qy), IQ_PARAMETERS); 
    117124    """ 
    118125 
    … …  
    122129demo = dict( 
    123130        scale=1, background=0, 
    124                 Nlayers=20,spacing=200., 
    125         Caille_parameter=0.05, 
    126                 tail_length=15,head_length=10, 
    127         sld=-1, head_sld=4.0, solvent_sld=6.0, 
    128                 tail_length_pd= 0.1, tail_length_pd_n=20, 
    129                 head_length_pd= 0.05, head_length_pd_n=30, 
    130                 spacing_pd= 0.2, spacing_pd_n=40 
    131          ) 
     131        Nlayers=20, 
     132        spacing=200., Caille_parameter=0.05, 
     133        tail_length=15,head_length=10, 
     134        #sld=-1, head_sld=4.0, solvent_sld=6.0, 
     135        sld=-1, head_sld=4.1, solvent_sld=6.0, 
     136        tail_length_pd= 0.1, tail_length_pd_n=20, 
     137        head_length_pd= 0.05, head_length_pd_n=30, 
     138        spacing_pd= 0.2, spacing_pd_n=40 
     139   ) 
    132140 
    133141oldname = 'LamellarPSHGModel' 

  • sasmodels/models/lamellarCailleHG_kernel.c

    rcd55ac3 r12c810f  
    88      double Nlayers,  
    99      double dd, 
    10           double Cp,  
     10      double Cp, 
    1111      double tail_sld, 
    1212      double head_sld, 
    … …  
    1818      double Nlayers,  
    1919      double dd, 
    20           double Cp,  
     20      double Cp, 
    2121      double tail_sld, 
    2222      double head_sld, 
    … …  
    3434  NN = trunc(Nlayers);    //be sure that NN is an integer 
    3535   
    36   Pq = (head_sld-solvent_sld)*(sin(qval*(head_length+tail_length))-sin(qval*tail_length)) +  
     36  Pq = (head_sld-solvent_sld)*(sin(qval*(head_length+tail_length))-sin(qval*tail_length)) + 
    3737              (tail_sld-solvent_sld)*sin(qval*tail_length); 
    3838  Pq *= Pq; 
    … …  
    4242  ii=0; 
    4343  Sq = 0.0; 
    44   for(ii=1;ii<(NNint-1);ii+=1) { 
     44  for(ii=1;ii<=(NNint-1);ii+=1) { 
    4545 
    4646    //fii = (double)ii;   //do I really need to do this? - unused variable, removed 18Feb2015 
    … …  
    5555    //temp *= cos(dd*qval*ii/(1.0+t1)); 
    5656    temp *= cos(dd*qval*ii); 
     57    //if (temp < 0) printf("q=%g: ii=%d, cos(dd*q*ii)=cos(%g) < 0\n",qval,ii,dd*qval*ii); 
    5758    //temp *= exp(-1.0*(t2 + t3)/(2.0*(1.0+t1)) ); 
    5859    temp *= exp(-t2/2.0 ); 
    … …  
    6566  Sq += 1.0; 
    6667 
     68  //if (Sq < 0) printf("q=%g: S(q) =%g\n", qval, Sq); 
     69 
    6770  inten = 2.0*M_PI*Pq*Sq/(dd*qval*qval); 
    6871 
    6972  inten *= 1.0e-04;   // 1/A to 1/cm 
    70  
    7173  return(inten); 
    7274} 

  • sasmodels/models/lamellarFFHG.py

    rdc02af0 r3c56da87  
    11# Note: model title and parameter table are inserted automatically 
    22r""" 
    3 This model provides the scattering intensity, *I(q)*, for a lyotropic lamellar phase where a random distribution in 
    4 solution are assumed. The SLD of the head region is taken to be different from the SLD of the tail region. 
     3This model provides the scattering intensity, *I(q)*, for a lyotropic lamellar 
     4phase where a random distribution in solution are assumed. The SLD of the head 
     5region is taken to be different from the SLD of the tail region. 
    56 
    67*2.1.31.1. Definition* 
    … …  
    1617.. image:: img/lamellarFFHG_.jpg 
    1718 
    18 where |delta|\ T = tail length (or *tail_length*), |delta|\ H = head thickness (or *h_thickness*), 
    19 |drho|\ H = SLD(headgroup) - SLD(solvent), and |drho|\ T = SLD(tail) - SLD(solvent). 
     19where |delta|\ T = tail length (or *tail_length*), |delta|\ H = head thickness 
     20(or *h_thickness*), |drho|\ H = SLD(headgroup) - SLD(solvent), 
     21and |drho|\ T = SLD(tail) - SLD(solvent). 
    2022 
    21 The 2D scattering intensity is calculated in the same way as 1D, where the *q* vector is defined as 
     23The 2D scattering intensity is calculated in the same way as 1D, where 
     24the *q* vector is defined as 
    2225 
    2326.. math:: 
    … …  
    2528    Q = \sqrt{Q_x^2 + Q_y^2} 
    2629 
    27 The returned value is in units of |cm^-1|, on absolute scale. In the parameters, *sld_tail* = SLD of the tail group, 
    28 and *sld_head* = SLD of the head group. 
    29  
    30 ==============  ========  ============= 
    31 Parameter name  Units     Default value 
    32 ==============  ========  ============= 
    33 background      |cm^-1|   0.0 
    34 head_sld        |Ang^-2|  3e-06 
    35 scale            None      1 
    36 solvent_sld     |Ang^-2|  6e-06 
    37 head_length     |Ang|     10 
    38 tail_length     |Ang|     15 
    39 sld (tail)      |Ang^-2|  0 
    40 ==============  ========  ============= 
     30The returned value is in units of |cm^-1|, on absolute scale. In the 
     31parameters, *sld_tail* = SLD of the tail group, and *sld_head* = SLD 
     32of the head group. 
    4133 
    4234.. image:: img/lamellarFFHG_138.jpg 
    … …  
    4436*Figure. 1D plot using the default values (w/1000 data point).* 
    4537 
    46 Our model uses the form factor calculations implemented in a c-library provided by the NIST Center for Neutron Research 
    47 (Kline, 2006). 
     38Our model uses the form factor calculations implemented in a C library 
     39provided by the NIST Center for Neutron Research (Kline, 2006). 
    4840 
    4941REFERENCE 
    … …  
    5648""" 
    5749 
    58 from numpy import pi, inf 
     50from numpy import inf 
    5951 
    6052name = "lamellar_FFHG" 
    … …  
    7163                scale = scale factor 
    7264""" 
     65category = "shape:lamellae" 
    7366 
    7467parameters = [ 
    … …  
    9588Iq = """ 
    9689    const double qsq = q*q; 
    97         const double drh = head_sld - solvent_sld; 
    98         const double drt = sld  - solvent_sld;          //correction 13FEB06 by L.Porcar 
    99         const double qT = q*tail_length; 
    100         double Pq, inten; 
    101         Pq = drh*(sin(q*(head_length+tail_length))-sin(qT)) + drt*sin(qT); 
    102         Pq *= Pq; 
    103         Pq *= 4.0/(qsq); 
    104          
    105         inten = 2.0e-4*M_PI*Pq/qsq;      
    106          
    107         return inten /= 2.0*(head_length+tail_length);  //normalize by the bilayer thickness     
    108                  
     90    const double drh = head_sld - solvent_sld; 
     91    const double drt = sld - solvent_sld;    //correction 13FEB06 by L.Porcar 
     92    const double qT = q*tail_length; 
     93    double Pq, inten; 
     94    Pq = drh*(sin(q*(head_length+tail_length))-sin(qT)) + drt*sin(qT); 
     95    Pq *= Pq; 
     96    Pq *= 4.0/(qsq); 
     97 
     98    inten = 2.0e-4*M_PI*Pq/qsq; 
     99 
     100    // normalize by the bilayer thickness 
     101    inten /= 2.0*(head_length+tail_length); 
     102 
     103    return inten; 
    109104    """ 
    110          
     105 
    111106Iqxy = """ 
    112     // never called since no orientation or magnetic parameters. 
    113     return -1.0; 
     107    return Iq(sqrt(qx*qx+qy*qy), IQ_PARAMETERS); 
    114108    """ 
    115109 
    … …  
    118112 
    119113demo = dict( 
    120         scale=1, background=0, 
    121                 tail_length=15,head_length=10, 
    122         sld=0.4, head_sld=3.0, solvent_sld=6.0, 
    123                 tail_length_pd= 0.2, tail_length_pd_n=40, 
    124                 head_length_pd= 0.01, head_length_pd_n=40, 
    125           ) 
     114    scale=1, background=0, 
     115    tail_length=15,head_length=10, 
     116    sld=0.4, head_sld=3.0, solvent_sld=6.0, 
     117    tail_length_pd= 0.2, tail_length_pd_n=40, 
     118    head_length_pd= 0.01, head_length_pd_n=40, 
     119    ) 
     120 
    126121oldname = 'LamellarFFHGModel' 
    127 oldpars = dict(head_length='h_thickness', sld='sld_tail', head_sld='sld_head', solvent_sld='sld_solvent') 
     122oldpars = dict(head_length='h_thickness', sld='sld_tail', 
     123               head_sld='sld_head', solvent_sld='sld_solvent') 
    128124 

  • sasmodels/models/lamellarPC.py

    rdc02af0 r3c56da87  
    11# Note: model title and parameter table are inserted automatically 
    22r""" 
    3 This model calculates the scattering from a stack of repeating lamellar structures. The stacks of lamellae (infinite 
    4 in lateral dimension) are treated as a paracrystal to account for the repeating spacing. The repeat distance is further 
    5 characterized by a Gaussian polydispersity. **This model can be used for large multilamellar vesicles.** 
     3This model calculates the scattering from a stack of repeating lamellar 
     4structures. The stacks of lamellae (infinite in lateral dimension) are 
     5treated as a paracrystal to account for the repeating spacing. The repeat 
     6distance is further characterized by a Gaussian polydispersity. **This model 
     7can be used for large multilamellar vesicles.** 
    68 
    79*2.1.33.1. Definition* 
    … …  
    1113.. image:: img/image145.jpg 
    1214 
    13 The form factor of the bilayer is approximated as the cross section of an infinite, planar bilayer of thickness *t* 
     15The form factor of the bilayer is approximated as the cross section of an 
     16infinite, planar bilayer of thickness *t* 
    1417 
    1518.. image:: img/image146.jpg 
    1619 
    17 Here, the scale factor is used instead of the mass per area of the bilayer (*G*). The scale factor is the volume 
    18 fraction of the material in the bilayer, *not* the total excluded volume of the paracrystal. *Z*\ :sub:`N`\ *(q)* 
    19 describes the interference effects for aggregates consisting of more than one bilayer. The equations used are (3-5) 
     20Here, the scale factor is used instead of the mass per area of the 
     21bilayer (*G*). The scale factor is the volume fraction of the material in 
     22the bilayer, *not* the total excluded volume of the paracrystal. 
     23*Z*\ :sub:`N`\ *(q)* describes the interference effects for aggregates 
     24consisting of more than one bilayer. The equations used are (3-5) 
    2025from the Bergstrom reference below. 
    2126 
    22 Non-integer numbers of stacks are calculated as a linear combination of the lower and higher values 
     27Non-integer numbers of stacks are calculated as a linear combination of 
     28the lower and higher values 
    2329 
    2430.. image:: img/image147.jpg 
    2531 
    26 The 2D scattering intensity is the same as 1D, regardless of the orientation of the *q* vector which is defined as 
     32The 2D scattering intensity is the same as 1D, regardless of the orientation 
     33of the *q* vector which is defined as 
    2734 
    2835.. math:: 
    … …  
    3037    Q = \sqrt{Q_x^2 + Q_y^2} 
    3138 
    32 The parameters of the model are *Nlayers* = no. of layers, and *pd_spacing* = polydispersity of spacing. 
     39The parameters of the model are *Nlayers* = no. of layers, and 
     40*pd_spacing* = polydispersity of spacing. 
    3341 
    3442==============  ========  ============= 
    … …  
    4957*Figure. 1D plot using the default values above (w/20000 data point).* 
    5058 
    51 Our model uses the form factor calculations implemented in a c-library provided by the NIST Center for Neutron Research 
    52 (Kline, 2006). 
     59Our model uses the form factor calculations implemented in a C library 
     60provided by the NIST Center for Neutron Research (Kline, 2006). 
    5361 
    5462REFERENCE 
    5563 
    56 M Bergstrom, J S Pedersen, P Schurtenberger, S U Egelhaaf, *J. Phys. Chem. B*, 103 (1999) 9888-9897 
     64M Bergstrom, J S Pedersen, P Schurtenberger, S U Egelhaaf, 
     65*J. Phys. Chem. B*, 103 (1999) 9888-9897 
    5766 
    5867""" 
    5968 
    60 from numpy import pi, inf 
     69from numpy import inf 
    6170 
    6271name = "lamellarPC" 
    … …  
    7180                scale = scale factor 
    7281""" 
     82category = "shape:lamellae" 
    7383 
    7484parameters = [ 
    … …  
    8999    ] 
    90100 
    91          
     101 
    92102source = [ "lamellarPC_kernel.c"] 
    93103 
    94 # No volume normalization despite having a volume parameter 
    95 # This should perhaps be volume normalized? 
    96104form_volume = """ 
    97105    return 1.0; 
    … …  
    99107 
    100108Iqxy = """ 
    101     // never called since no orientation or magnetic parameters. 
    102     return -1.0; 
     109    return Iq(sqrt(qx*qx+qy*qy), IQ_PARAMETERS); 
    103110    """ 
    104111 
    … …  
    107114 
    108115demo = dict( 
    109         scale=1, background=0, 
    110                 thickness=33,Nlayers=20,spacing=250,spacing_polydisp=0.2, 
    111         sld=1.0, solvent_sld=6.34, 
    112                 thickness_pd= 0.2, thickness_pd_n=40 
    113          ) 
     116    scale=1, background=0, 
     117    thickness=33, Nlayers=20, spacing=250, spacing_polydisp=0.2, 
     118    sld=1.0, solvent_sld=6.34, 
     119    thickness_pd= 0.2, thickness_pd_n=40 
     120    ) 
     121 
    114122oldname = 'LamellarPCrystalModel' 
    115 oldpars = dict(spacing_polydisp='pd_spacing', sld='sld_layer',solvent_sld='sld_solvent') 
     123oldpars = dict( 
     124    spacing_polydisp='pd_spacing', sld='sld_layer', 
     125    solvent_sld='sld_solvent' 
     126    ) 
    116127 
    117128 

  • sasmodels/models/lamellarPC_kernel.c

    rdc02af0 rf734e7d  
    2828        //get the fractional part of Nlayers, to determine the "mixing" of N's 
    2929         
    30         n1 = trunc(Nlayers);            //rounds towards zero 
     30        n1 = (long)trunc(Nlayers);              //rounds towards zero 
    3131        n2 = n1 + 1; 
    3232        xn = (double)n2 - Nlayers;                      //fractional contribution of n1 

  • sasmodels/models/parallelepiped.py

    rc5b7d07 r9890053  
    99---------- 
    1010 
    11 This model provides the form factor, *P(q)*, for a rectangular parallelepiped  
    12 (below) where the form factor is normalized by the volume of the  
    13 parallelepiped. If you need to apply polydispersity, see also the  
     11This model provides the form factor, *P(q)*, for a rectangular parallelepiped 
     12(below) where the form factor is normalized by the volume of the 
     13parallelepiped. If you need to apply polydispersity, see also the 
    1414RectangularPrismModel_. 
    1515 
    … …  
    2020    P(Q) = {\text{scale} \over V} F^2(Q) + \text{background} 
    2121 
    22 where the volume *V* = *A B C* and the averaging < > is applied over all  
     22where the volume *V* = *A B C* and the averaging < > is applied over all 
    2323orientations for 1D. 
    2424 
    … …  
    2727*Figure. Parallelepiped with the corresponding Definition of sides. 
    2828 
    29 The edge of the solid must satisfy the condition that** *A* < *B* < *C*.  
    30 Then, assuming *a* = *A* / *B* < 1, *b* = *B* / *B* = 1, and  
     29The edge of the solid must satisfy the condition that** *A* < *B* < *C*. 
     30Then, assuming *a* = *A* / *B* < 1, *b* = *B* / *B* = 1, and 
    3131*c* = *C* / *B* > 1, the form factor is 
    3232 
    3333.. math:: 
    3434 
    35     P(q) = \frac{\textstyle{scale}}{V}\int_0^1 \phi(\mu \sqrt{1-\sigma^2},a)  
     35    P(q) = \frac{\textstyle{scale}}{V}\int_0^1 \phi(\mu \sqrt{1-\sigma^2},a) 
    3636    [S(\mu c \sigma/2)]^2 d\sigma 
    3737 
    … …  
    4040.. math:: 
    4141 
    42     \phi(\mu,a) = \int_0^1 \{S[\frac{\mu}{2}\cos(\frac{\pi}{2}u)]  
     42    \phi(\mu,a) = \int_0^1 \{S[\frac{\mu}{2}\cos(\frac{\pi}{2}u)] 
    4343    S[\frac{\mu a}{2}\sin(\frac{\pi}{2}u)]\}^2 du 
    4444 
    4545    S(x) = \frac{\sin x}{x} 
    46      
     46 
    4747    \mu = qB 
    4848 
    … …  
    5353    \Delta\rho = \rho_{\textstyle p} - \rho_{\textstyle solvent} 
    5454 
    55 The scattering intensity per unit volume is returned in units of |cm^-1|;  
     55The scattering intensity per unit volume is returned in units of |cm^-1|; 
    5656ie, *I(q)* = |phi| *P(q)*\ . 
    5757 
    58 NB: The 2nd virial coefficient of the parallelpiped is calculated based on  
    59 the averaged effective radius (= sqrt(*short_a* \* *short_b* / |pi|)) and  
     58NB: The 2nd virial coefficient of the parallelpiped is calculated based on 
     59the averaged effective radius (= sqrt(*short_a* \* *short_b* / |pi|)) and 
    6060length(= *long_c*) values, and used as the effective radius for 
    6161*S(Q)* when *P(Q)* \* *S(Q)* is applied. 
    6262 
    63 To provide easy access to the orientation of the parallelepiped, we define  
     63To provide easy access to the orientation of the parallelepiped, we define 
    6464three angles |theta|, |phi| and |bigpsi|. The definition of |theta| and |phi| 
    65 is the same as for the cylinder model (see also figures below).  
    66 The angle |bigpsi| is the rotational angle around the *long_c* axis against  
    67 the *q* plane. For example, |bigpsi| = 0 when the *short_b* axis is parallel  
     65is the same as for the cylinder model (see also figures below). 
     66The angle |bigpsi| is the rotational angle around the *long_c* axis against 
     67the *q* plane. For example, |bigpsi| = 0 when the *short_b* axis is parallel 
    6868to the *x*-axis of the detector. 
    6969 
    … …  
    8383---------- 
    8484 
    85 Validation of the code was done by comparing the output of the 1D calculation  
    86 to the angular average of the output of a 2D calculation over all possible  
    87 angles. The Figure below shows the comparison where the solid dot refers to  
    88 averaged 2D while the line represents the result of the 1D calculation (for  
    89 the averaging, 76, 180, 76 points are taken for the angles of |theta|, |phi|,  
     85Validation of the code was done by comparing the output of the 1D calculation 
     86to the angular average of the output of a 2D calculation over all possible 
     87angles. The Figure below shows the comparison where the solid dot refers to 
     88averaged 2D while the line represents the result of the 1D calculation (for 
     89the averaging, 76, 180, 76 points are taken for the angles of |theta|, |phi|, 
    9090and |psi| respectively). 
    9191 
    … …  
    9696*Figure. Comparison between 1D and averaged 2D.* 
    9797 
    98 This model reimplements the form factor calculations implemented in a c-library  
     98This model reimplements the form factor calculations implemented in a c-library 
    9999provided by the NIST Center for Neutron Research (Kline, 2006). 
    100100 
    101101""" 
    102102 
    103 from numpy import pi, inf 
     103from numpy import pi, inf, sqrt 
    104104 
    105105name = "parallelepiped" 
    … …  
    108108     P(q)= scale/V*integral from 0 to 1 of ... 
    109109           phi(mu*sqrt(1-sigma^2),a) * S(mu*c*sigma/2)^2 * dsigma 
    110       
     110 
    111111            phi(mu,a) = integral from 0 to 1 of .. 
    112             (S((mu/2)*cos(pi*u/2))*S((mu*a/2)*sin(pi*u/2)))^2 * du 
     112        (S((mu/2)*cos(pi*u/2))*S((mu*a/2)*sin(pi*u/2)))^2 * du 
    113113            S(x) = sin(x)/x 
    114             mu = q*B 
     114        mu = q*B 
    115115""" 
     116category = "shape:parallelpiped" 
    116117 
    117118parameters = [ 
    118 #   [ "name", "units", default, [lower, upper], "type", 
    119 #     "description" ], 
    120     [ "sld", "6e-6/Ang^2", 4, [-inf,inf], "", 
    121       "Parallelepiped scattering length density" ], 
    122     [ "solvent_sld", "1e-6/Ang^2", 1, [-inf,inf], "", 
    123       "Solvent scattering length density" ], 
    124     [ "a_side", "Ang", 35, [0, inf], "volume", 
    125       "Shorter side of the parallelepiped" ], 
    126     [ "b_side", "Ang", 75, [0, inf], "volume", 
    127       "Second side of the parallelepiped" ], 
    128     [ "c_side", "Ang", 400, [0, inf], "volume", 
    129       "Larger side of the parallelepiped" ], 
    130     [ "theta", "degrees", 60, [-inf, inf], "orientation", 
    131       "In plane angle" ], 
    132     [ "phi", "degrees", 60, [-inf, inf], "orientation", 
    133       "Out of plane angle" ], 
    134     [ "psi", "degrees", 60, [-inf, inf], "orientation", 
    135       "Rotation angle around its own c axis against q plane" ], 
     119    #   [ "name", "units", default, [lower, upper], "type", 
     120    #     "description" ], 
     121    ["sld", "1e-6/Ang^2", 4, [-inf, inf], "", 
     122     "Parallelepiped scattering length density"], 
     123    ["solvent_sld", "1e-6/Ang^2", 1, [-inf, inf], "", 
     124     "Solvent scattering length density"], 
     125    ["a_side", "Ang", 35, [0, inf], "volume", 
     126     "Shorter side of the parallelepiped"], 
     127    ["b_side", "Ang", 75, [0, inf], "volume", 
     128     "Second side of the parallelepiped"], 
     129    ["c_side", "Ang", 400, [0, inf], "volume", 
     130     "Larger side of the parallelepiped"], 
     131    ["theta", "degrees", 60, [-inf, inf], "orientation", 
     132     "In plane angle"], 
     133    ["phi", "degrees", 60, [-inf, inf], "orientation", 
     134     "Out of plane angle"], 
     135    ["psi", "degrees", 60, [-inf, inf], "orientation", 
     136     "Rotation angle around its own c axis against q plane"], 
    136137    ] 
    137138 
    138 source = [ "lib/J1.c", "lib/gauss76.c", "parallelepiped.c" ] 
     139source = ["lib/J1.c", "lib/gauss76.c", "parallelepiped.c"] 
    139140 
    140141def ER(a_side, b_side, c_side): 
    … …  
    147148    b = 0.5 * c_side 
    148149    t1 = a * a * b 
    149     t2 = 1.0 + (b/a)*(1.0+a/b/2.0)*(1.0+pi*a/b/2.0) 
     150    t2 = 1.0 + (b / a) * (1.0 + a / b / 2.0) * (1.0 + pi * a / b / 2.0) 
    150151    ddd = 3.0 * t1 * t2 
    151152 
    152     return 0.5 * (ddd)**(1./3.) 
     153    return 0.5 * (ddd) ** (1. / 3.) 
    153154 
    154155# parameters for demo 
    … …  
    168169# For testing against the old sasview models, include the converted parameter 
    169170# names and the target sasview model name. 
    170 oldname='ParallelepipedModel' 
    171 oldpars=dict(theta='parallel_theta', phi='parallel_phi', psi='parallel_psi', 
    172              a_side='short_a', b_side='short_b', c_side='long_c', 
    173              sld='sldPipe', solvent_sld='sldSolv') 
     171oldname = 'ParallelepipedModel' 
     172oldpars = dict(theta='parallel_theta', phi='parallel_phi', psi='parallel_psi', 
     173               a_side='short_a', b_side='short_b', c_side='long_c', 
     174               sld='sldPipe', solvent_sld='sldSolv') 
    174175 

  • sasmodels/models/sphere.py

    ra503bfd r3c56da87  
    5858""" 
    5959 
    60 from numpy import pi, inf 
     60from numpy import inf 
    6161 
    6262name = "sphere" 
    … …  
    7070    solvent_sld: the SLD of the solvent 
    7171""" 
     72category = "shape:sphere" 
    7273 
    7374parameters = [ 

  • sasmodels/models/spherepy.py

    ra503bfd r3c56da87  
    5959 
    6060import numpy as np 
    61 from numpy import pi, inf, sin, cos, sqrt, exp, log 
     61from numpy import pi, inf, sin, cos, sqrt, log 
    6262 
    6363name = "sphere" 
    … …  
    7171    solvent_sld: the SLD of the solvent 
    7272""" 
     73category = "shape:sphere" 
    7374 
    7475parameters = [ 

  • sasmodels/models/stickyhardsphere.py

    r9cb1415 r3c56da87  
    11# Note: model title and parameter table are inserted automatically 
    2 r"""This calculates the interparticle structure factor for a hard sphere fluid with a narrow attractive well. A perturbative 
    3 solution of the Percus-Yevick closure is used. The strength of the attractive well is described in terms of "stickiness" 
    4 as defined below. The returned value is a dimensionless structure factor, *S(q)*. 
     2r""" 
     3This calculates the interparticle structure factor for a hard sphere fluid 
     4with a narrow attractive well. A perturbative solution of the Percus-Yevick 
     5closure is used. The strength of the attractive well is described in terms 
     6of "stickiness" as defined below. The returned value is a dimensionless 
     7structure factor, *S(q)*. 
    58 
    6 The perturb (perturbation parameter), |epsilon|, should be held between 0.01 and 0.1. It is best to hold the 
    7 perturbation parameter fixed and let the "stickiness" vary to adjust the interaction strength. The stickiness, |tau|, 
    8 is defined in the equation below and is a function of both the perturbation parameter and the interaction strength. 
    9 |tau| and |epsilon| are defined in terms of the hard sphere diameter (|sigma| = 2\*\ *R*\ ), the width of the square 
    10 well, |bigdelta| (same units as *R*), and the depth of the well, *Uo*, in units of kT. From the definition, it is clear 
    11 that smaller |tau| means stronger attraction. 
     9The perturb (perturbation parameter), |epsilon|, should be held between 0.01 
     10and 0.1. It is best to hold the perturbation parameter fixed and let 
     11the "stickiness" vary to adjust the interaction strength. The stickiness, 
     12|tau|, is defined in the equation below and is a function of both the 
     13perturbation parameter and the interaction strength. |tau| and |epsilon| 
     14are defined in terms of the hard sphere diameter (|sigma| = 2\*\ *R*\ ), the 
     15width of the square well, |bigdelta| (same units as *R*), and the depth of 
     16the well, *Uo*, in units of kT. From the definition, it is clear that 
     17smaller |tau| means stronger attraction. 
    1218 
    1319.. image:: img/stickyhardsphere_228.PNG 
    … …  
    1723.. image:: img/stickyhardsphere_229.PNG 
    1824 
    19 The Percus-Yevick (PY) closure was used for this calculation, and is an adequate closure for an attractive interparticle 
    20 potential. This solution has been compared to Monte Carlo simulations for a square well fluid, with good agreement. 
     25The Percus-Yevick (PY) closure was used for this calculation, and is an 
     26adequate closure for an attractive interparticle potential. This solution 
     27has been compared to Monte Carlo simulations for a square well fluid, with 
     28good agreement. 
    2129 
    22 The true particle volume fraction, |phi|, is not equal to *h*, which appears in most of the reference. The two are 
    23 related in equation (24) of the reference. The reference also describes the relationship between this perturbation 
    24 solution and the original sticky hard sphere (or adhesive sphere) model by Baxter. 
     30The true particle volume fraction, |phi|, is not equal to *h*, which appears 
     31in most of the reference. The two are related in equation (24) of the 
     32reference. The reference also describes the relationship between this 
     33perturbation solution and the original sticky hard sphere (or adhesive 
     34sphere) model by Baxter. 
    2535 
    26 NB: The calculation can go haywire for certain combinations of the input parameters, producing unphysical solutions - in 
    27 this case errors are reported to the command window and the *S(q)* is set to -1 (so it will disappear on a log-log 
    28 plot). Use tight bounds to keep the parameters to values that you know are physical (test them) and keep nudging them 
    29 until the optimization does not hit the constraints. 
     36NB: The calculation can go haywire for certain combinations of the input 
     37parameters, producing unphysical solutions - in this case errors are 
     38reported to the command window and the *S(q)* is set to -1 (so it will 
     39disappear on a log-log plot). Use tight bounds to keep the parameters to 
     40values that you know are physical (test them) and keep nudging them until 
     41the optimization does not hit the constraints. 
    3042 
    31 In sasview the effective radius will be calculated from the parameters used in the form factor P(Q) that this  
    32 S(Q) is combined with. 
     43In sasview the effective radius will be calculated from the parameters 
     44used in the form factor P(Q) that this S(Q) is combined with. 
    3345 
    34 For 2D data: The 2D scattering intensity is calculated in the same way as 1D, where the *q* vector is defined as 
     46For 2D data: The 2D scattering intensity is calculated in the same way 
     47as 1D, where the *q* vector is defined as 
    3548 
    3649.. math:: 
    … …  
    5770 
    5871# TODO: refactor so that we pull in the old sansmodels.c_extensions 
    59   
    60 from numpy import pi, inf 
     72 
     73from numpy import inf 
    6174 
    6275name = "stickyhardsphere" 
    … …  
    6477description = """\ 
    6578        [Sticky hard sphere structure factor, with Percus-Yevick closure] 
    66         Interparticle structure factor S(Q)for a hard sphere fluid with  
     79        Interparticle structure factor S(Q)for a hard sphere fluid with 
    6780                a narrow attractive well. Fits are prone to deliver non-physical 
    68                 parameters, use with care and read the references in the full manual.  
    69                 In sasview the effective radius will be calculated from the  
     81                parameters, use with care and read the references in the full manual. 
     82                In sasview the effective radius will be calculated from the 
    7083                parameters used in P(Q). 
    7184""" 
     85category = "structure-factor" 
    7286 
    7387parameters = [ 
    74 #   [ "name", "units", default, [lower, upper], "type", 
    75 #     "description" ], 
    76     [ "effect_radius", "Ang", 50.0, [0, inf], "volume", 
    77       "effective radius of hard sphere" ], 
    78     [ "volfraction", "", 0.2, [0, 0.74], "", 
    79       "volume fraction of hard spheres" ], 
    80     [ "perturb", "", 0.05, [0.01, 0.1], "", 
    81       "perturbation parameter, epsilon" ], 
    82     [ "stickiness", "",  0.20, [-inf,inf], "", 
    83       "stickiness, tau" ], 
     88    #   [ "name", "units", default, [lower, upper], "type", 
     89    #     "description" ], 
     90    ["effect_radius", "Ang", 50.0, [0, inf], "volume", 
     91     "effective radius of hard sphere"], 
     92    ["volfraction", "", 0.2, [0, 0.74], "", 
     93     "volume fraction of hard spheres"], 
     94    ["perturb", "", 0.05, [0.01, 0.1], "", 
     95     "perturbation parameter, epsilon"], 
     96    ["stickiness", "", 0.20, [-inf, inf], "", 
     97     "stickiness, tau"], 
    8498    ] 
    85          
     99 
    86100# No volume normalization despite having a volume parameter 
    87101# This should perhaps be volume normalized? 
    … …  
    91105 
    92106Iq = """ 
    93         double onemineps,eta; 
    94         double sig,aa,etam1,etam1sq,qa,qb,qc,radic; 
    95         double lam,lam2,test,mu,alpha,beta; 
    96         double kk,k2,k3,ds,dc,aq1,aq2,aq3,aq,bq1,bq2,bq3,bq,sq; 
    97          
    98         onemineps = 1.0-perturb; 
    99         eta = volfraction/onemineps/onemineps/onemineps; 
    100          
    101         sig = 2.0 * effect_radius; 
    102         aa = sig/onemineps; 
    103         etam1 = 1.0 - eta; 
    104         etam1sq=etam1*etam1; 
    105         //C 
    106         //C  SOLVE QUADRATIC FOR LAMBDA 
    107         //C 
    108         qa = eta/12.0; 
    109         qb = -1.0*(stickiness + eta/etam1); 
    110         qc = (1.0 + eta/2.0)/etam1sq; 
    111         radic = qb*qb - 4.0*qa*qc; 
    112         if(radic<0) { 
    113                 //if(x>0.01 && x<0.015) 
    114                 //      Print "Lambda unphysical - both roots imaginary" 
    115                 //endif 
    116                 return(-1.0); 
    117         } 
    118         //C   KEEP THE SMALLER ROOT, THE LARGER ONE IS UNPHYSICAL 
    119         lam = (-1.0*qb-sqrt(radic))/(2.0*qa); 
    120         lam2 = (-1.0*qb+sqrt(radic))/(2.0*qa); 
    121         if(lam2<lam) { 
    122                 lam = lam2; 
    123         } 
    124         test = 1.0 + 2.0*eta; 
    125         mu = lam*eta*etam1; 
    126         if(mu>test) { 
    127                 //if(x>0.01 && x<0.015) 
    128                 // Print "Lambda unphysical mu>test" 
    129                 //endif 
    130                 return(-1.0); 
    131         } 
    132         alpha = (1.0 + 2.0*eta - mu)/etam1sq; 
    133         beta = (mu - 3.0*eta)/(2.0*etam1sq); 
    134         //C 
    135         //C   CALCULATE THE STRUCTURE FACTOR 
    136         //C 
    137         kk = q*aa; 
    138         k2 = kk*kk; 
    139         k3 = kk*k2; 
    140         SINCOS(kk,ds,dc); 
    141         //ds = sin(kk); 
    142         //dc = cos(kk); 
    143         aq1 = ((ds - kk*dc)*alpha)/k3; 
    144         aq2 = (beta*(1.0-dc))/k2; 
    145         aq3 = (lam*ds)/(12.0*kk); 
    146         aq = 1.0 + 12.0*eta*(aq1+aq2-aq3); 
    147         // 
    148         bq1 = alpha*(0.5/kk - ds/k2 + (1.0 - dc)/k3); 
    149         bq2 = beta*(1.0/kk - ds/k2); 
    150         bq3 = (lam/12.0)*((1.0 - dc)/kk); 
    151         bq = 12.0*eta*(bq1+bq2-bq3); 
    152         // 
    153         sq = 1.0/(aq*aq +bq*bq); 
    154          
    155         return(sq); 
     107    double onemineps,eta; 
     108    double sig,aa,etam1,etam1sq,qa,qb,qc,radic; 
     109    double lam,lam2,test,mu,alpha,beta; 
     110    double kk,k2,k3,ds,dc,aq1,aq2,aq3,aq,bq1,bq2,bq3,bq,sq; 
     111 
     112    onemineps = 1.0-perturb; 
     113    eta = volfraction/onemineps/onemineps/onemineps; 
     114 
     115    sig = 2.0 * effect_radius; 
     116    aa = sig/onemineps; 
     117    etam1 = 1.0 - eta; 
     118    etam1sq=etam1*etam1; 
     119    //C 
     120    //C  SOLVE QUADRATIC FOR LAMBDA 
     121    //C 
     122    qa = eta/12.0; 
     123    qb = -1.0*(stickiness + eta/etam1); 
     124    qc = (1.0 + eta/2.0)/etam1sq; 
     125    radic = qb*qb - 4.0*qa*qc; 
     126    if(radic<0) { 
     127        //if(x>0.01 && x<0.015) 
     128        //      Print "Lambda unphysical - both roots imaginary" 
     129        //endif 
     130        return(-1.0); 
     131    } 
     132    //C   KEEP THE SMALLER ROOT, THE LARGER ONE IS UNPHYSICAL 
     133    lam = (-1.0*qb-sqrt(radic))/(2.0*qa); 
     134    lam2 = (-1.0*qb+sqrt(radic))/(2.0*qa); 
     135    if(lam2<lam) { 
     136        lam = lam2; 
     137    } 
     138    test = 1.0 + 2.0*eta; 
     139    mu = lam*eta*etam1; 
     140    if(mu>test) { 
     141        //if(x>0.01 && x<0.015) 
     142        // Print "Lambda unphysical mu>test" 
     143        //endif 
     144        return(-1.0); 
     145    } 
     146    alpha = (1.0 + 2.0*eta - mu)/etam1sq; 
     147    beta = (mu - 3.0*eta)/(2.0*etam1sq); 
     148    //C 
     149    //C   CALCULATE THE STRUCTURE FACTOR 
     150    //C 
     151    kk = q*aa; 
     152    k2 = kk*kk; 
     153    k3 = kk*k2; 
     154    SINCOS(kk,ds,dc); 
     155    //ds = sin(kk); 
     156    //dc = cos(kk); 
     157    aq1 = ((ds - kk*dc)*alpha)/k3; 
     158    aq2 = (beta*(1.0-dc))/k2; 
     159    aq3 = (lam*ds)/(12.0*kk); 
     160    aq = 1.0 + 12.0*eta*(aq1+aq2-aq3); 
     161    // 
     162    bq1 = alpha*(0.5/kk - ds/k2 + (1.0 - dc)/k3); 
     163    bq2 = beta*(1.0/kk - ds/k2); 
     164    bq3 = (lam/12.0)*((1.0 - dc)/kk); 
     165    bq = 12.0*eta*(bq1+bq2-bq3); 
     166    // 
     167    sq = 1.0/(aq*aq +bq*bq); 
     168 
     169    return(sq); 
    156170""" 
     171 
    157172Iqxy = """ 
    158     // never called since no orientation or magnetic parameters. 
    159     return -1.0; 
     173    return Iq(sqrt(qx*qx+qy*qy), IQ_PARAMETERS); 
    160174    """ 
    161175 
    … …  
    165179oldname = 'StickyHSStructure' 
    166180oldpars = dict() 
    167 demo = dict(effect_radius = 200,volfraction = 0.2,perturb=0.05,stickiness=0.2,effect_radius_pd = 0.1,effect_radius_pd_n = 40) 
     181demo = dict(effect_radius=200, volfraction=0.2, perturb=0.05, 
     182            stickiness=0.2, effect_radius_pd=0.1, effect_radius_pd_n=40) 
    168183 
    169  
    170  

  • sasmodels/models/triaxial_ellipsoid.py

    ra503bfd r3c56da87  
    9090""" 
    9191 
    92 from numpy import pi, inf 
     92from numpy import inf 
    9393 
    9494name = "triaxial_ellipsoid" 
    … …  
    100100        not be correct. 
    101101""" 
     102category = "shape:ellipsoid" 
    102103 
    103104parameters = [ 

  • sasmodels/sasview_model.py

    r0a82216 r3c56da87  
    2626try: 
    2727    from .kernelcl import load_model 
    28 except ImportError,exc: 
     28except ImportError, exc: 
    2929    warnings.warn(str(exc)) 
    3030    warnings.warn("using ctypes instead") 
    … …  
    4141    will produce a class with a name compatible with SasView 
    4242    """ 
    43     model =  load_model(kernel_module, dtype=dtype) 
     43    model = load_model(kernel_module, dtype=dtype) 
    4444    def __init__(self, multfactor=1): 
    4545        SasviewModel.__init__(self, model) 
    4646    attrs = dict(__init__=__init__) 
    47     ConstructedModel = type(model.info[namestyle],  (SasviewModel,), attrs) 
     47    ConstructedModel = type(model.info[namestyle], (SasviewModel,), attrs) 
    4848    return ConstructedModel 
    4949 
    … …  
    6969        self.dispersion = dict() 
    7070        partype = model.info['partype'] 
    71         for name,units,default,limits,ptype,description in model.info['parameters']: 
     71        for name, units, default, limits, _, _ in model.info['parameters']: 
    7272            self.params[name] = default 
    73             self.details[name] = [units]+limits 
     73            self.details[name] = [units] + limits 
    7474 
    7575        for name in partype['pd-2d']: 
    … …  
    8383        self.orientation_params = ( 
    8484            partype['orientation'] 
    85             + [n+'.width' for n in partype['orientation']] 
     85            + [n + '.width' for n in partype['orientation']] 
    8686            + partype['magnetic']) 
    8787        self.magnetic_params = partype['magnetic'] 
    88         self.fixed = [n+'.width' for n in partype['pd-2d']] 
     88        self.fixed = [n + '.width' for n in partype['pd-2d']] 
    8989        self.non_fittable = [] 
    9090 
    9191        ## independent parameter name and unit [string] 
    92         self.input_name = model.info.get("input_name","Q") 
    93         self.input_unit = model.info.get("input_unit","A^{-1}") 
    94         self.output_name = model.info.get("output_name","Intensity") 
    95         self.output_unit = model.info.get("output_unit","cm^{-1}") 
     92        self.input_name = model.info.get("input_name", "Q") 
     93        self.input_unit = model.info.get("input_unit", "A^{-1}") 
     94        self.output_name = model.info.get("output_name", "Intensity") 
     95        self.output_unit = model.info.get("output_unit", "cm^{-1}") 
    9696 
    9797        ## _persistency_dict is used by sas.perspectives.fitting.basepage 
    … …  
    120120 
    121121 
     122    # pylint: disable=no-self-use 
    122123    def getProfile(self): 
    123124        """ 
    … …  
    139140        # Look for dispersion parameters 
    140141        toks = name.split('.') 
    141         if len(toks)==2: 
     142        if len(toks) == 2: 
    142143            for item in self.dispersion.keys(): 
    143                 if item.lower()==toks[0].lower(): 
     144                if item.lower() == toks[0].lower(): 
    144145                    for par in self.dispersion[item]: 
    145146                        if par.lower() == toks[1].lower(): 
    … …  
    149150            # Look for standard parameter 
    150151            for item in self.params.keys(): 
    151                 if item.lower()==name.lower(): 
     152                if item.lower() == name.lower(): 
    152153                    self.params[item] = value 
    153154                    return 
    … …  
    164165        # Look for dispersion parameters 
    165166        toks = name.split('.') 
    166         if len(toks)==2: 
     167        if len(toks) == 2: 
    167168            for item in self.dispersion.keys(): 
    168                 if item.lower()==toks[0].lower(): 
     169                if item.lower() == toks[0].lower(): 
    169170                    for par in self.dispersion[item]: 
    170171                        if par.lower() == toks[1].lower(): 
    … …  
    173174            # Look for standard parameter 
    174175            for item in self.params.keys(): 
    175                 if item.lower()==name.lower(): 
     176                if item.lower() == name.lower(): 
    176177                    return self.params[item] 
    177178 
    … …  
    182183        Return a list of all available parameters for the model 
    183184        """ 
    184         list = self.params.keys() 
     185        param_list = self.params.keys() 
    185186        # WARNING: Extending the list with the dispersion parameters 
    186         list.extend(self.getDispParamList()) 
    187         return list 
     187        param_list.extend(self.getDispParamList()) 
     188        return param_list 
    188189 
    189190    def getDispParamList(self): 
    … …  
    192193        """ 
    193194        # TODO: fix test so that parameter order doesn't matter 
    194         ret = ['%s.%s'%(d.lower(), p) 
     195        ret = ['%s.%s' % (d.lower(), p) 
    195196               for d in self._model.info['partype']['pd-2d'] 
    196197               for p in ('npts', 'nsigmas', 'width')] 
    … …  
    212213        **DEPRECATED**: use calculate_Iq instead 
    213214        """ 
    214         if isinstance(x, (list,tuple)): 
     215        if isinstance(x, (list, tuple)): 
     216            # pylint: disable=unpacking-non-sequence 
    215217            q, phi = x 
    216218            return self.calculate_Iq([q * math.cos(phi)], 
    … …  
    230232        **DEPRECATED**: use calculate_Iq instead 
    231233        """ 
    232         if isinstance(x, (list,tuple)): 
    233             return self.calculate_Iq([float(x[0])],[float(x[1])])[0] 
     234        if isinstance(x, (list, tuple)): 
     235            return self.calculate_Iq([float(x[0])], [float(x[1])])[0] 
    234236        else: 
    235237            return self.calculate_Iq([float(x)])[0] 
    … …  
    263265 
    264266 
    265         :param qdist: ndarray of scalar q-values or list [qx,qy] where qx,qy are 1D ndarrays 
    266         """ 
    267         if isinstance(qdist, (list,tuple)): 
     267        :param qdist: ndarray of scalar q-values or list [qx,qy] 
     268        where qx,qy are 1D ndarrays 
     269        """ 
     270        if isinstance(qdist, (list, tuple)): 
    268271            # Check whether we have a list of ndarrays [qx,qy] 
    269272            qx, qy = qdist 
    270273            partype = self._model.info['partype'] 
    271274            if not partype['orientation'] and not partype['magnetic']: 
    272                 return self.calculate_Iq(np.sqrt(qx**2+qy**2)) 
     275                return self.calculate_Iq(np.sqrt(qx ** 2 + qy ** 2)) 
    273276            else: 
    274277                return self.calculate_Iq(qx, qy) 
    … …  
    279282 
    280283        else: 
    281             raise TypeError("evalDistribution expects q or [qx, qy], not %r"%type(qdist)) 
     284            raise TypeError("evalDistribution expects q or [qx, qy], not %r" 
     285                            % type(qdist)) 
    282286 
    283287    def calculate_Iq(self, *args): 
    … …  
    313317            fv = ER(*values) 
    314318            #print values[0].shape, weights.shape, fv.shape 
    315             return np.sum(weights*fv) / np.sum(weights) 
     319            return np.sum(weights * fv) / np.sum(weights) 
    316320 
    317321    def calculate_VR(self): 
    … …  
    327331            vol_pars = self._model.info['partype']['volume'] 
    328332            values, weights = self._dispersion_mesh(vol_pars) 
    329             whole,part = VR(*values) 
    330             return np.sum(weights*part)/np.sum(weights*whole) 
     333            whole, part = VR(*values) 
     334            return np.sum(weights * part) / np.sum(weights * whole) 
    331335 
    332336    def set_dispersion(self, parameter, dispersion): 
    … …  
    373377 
    374378    def _get_weights(self, par): 
     379        """ 
     380            Return dispersion weights 
     381            :param par parameter name 
     382        """ 
    375383        from . import weights 
    376384 
    … …  
    378386        limits = self._model.info['limits'] 
    379387        dis = self.dispersion[par] 
    380         v,w = weights.get_weights( 
     388        value, weight = weights.get_weights( 
    381389            dis['type'], dis['npts'], dis['width'], dis['nsigmas'], 
    382390            self.params[par], limits[par], par in relative) 
    383         return v,w/w.max() 
    384  
     391        return value, weight / np.sum(weight) 
     392 

  • sasmodels/sesans.py

    rc97724e r3c56da87  
    2222# TODO: dead code; for now the call to the hankel transform happens in BumpsModel 
    2323class SesansCalculator: 
    24     def __init__(self, sans_kernel, q_zmax, Rmax, SElength, wavelength, thickness): 
    25         self._set_kernel(sans_kernel, q_zmax, Rmax) 
     24    def __init__(self, kernel, q_zmax, Rmax, SElength, wavelength, thickness): 
     25        self._set_kernel(kernel, q_zmax, Rmax) 
    2626        self.SElength = SElength 
    2727        self.wavelength = wavelength 
    2828        self.thickness = thickness 
    2929 
    30     def _set_kernel(self, sans_kernel, q_zmax, Rmax): 
    31         input = sans_kernel.make_input([make_q(q_zmax, Rmax)]) 
    32         self.sans_calculator = sans_kernel(input) 
     30    def _set_kernel(self, kernel, q_zmax, Rmax): 
     31        kernel_input = kernel.make_input([make_q(q_zmax, Rmax)]) 
     32        self.sans_calculator = kernel(kernel_input) 
    3333 
    3434    def __call__(self, pars, pd_pars, cutoff=1e-5): 

  • sasmodels/weights.py

    r5d4777d r3c56da87  
    2626        return pars 
    2727 
     28    # pylint: disable=no-self-use 
    2829    def set_weights(self, values, weights): 
    2930        raise RuntimeError("set_weights is only available for ArrayDispersion") 

  • setup.py

    rda32ec3 r040575f  
     1from setuptools import setup,find_packages 
     2 
     3packages = find_packages(exclude=['contrib', 'docs', 'tests*']) 
     4package_data = { 
     5    'sasmodels.models': ['*.c'], 
     6} 
     7required = [] 
     8 
     9setup( 
     10    name="sasmodels", 
     11    version = "1.0.0a", 
     12    description = "sasmodels package", 
     13    long_description=open('README.rst').read(), 
     14    author = "SasView Collaboration", 
     15    author_email = "management@sasview.org", 
     16    url = "http://www.sasview.org", 
     17    keywords = "small-angle x-ray and neutron scattering analysis", 
     18    download_url = "https://github.com/SasView/sasmodels", 
     19    classifiers=[ 
     20        'Development Status :: 4 - Beta', 
     21        'Environment :: Console', 
     22        'Intended Audience :: Science/Research', 
     23        'License :: Public Domain', 
     24        'Operating System :: OS Independent', 
     25        'Programming Language :: Python', 
     26        'Topic :: Scientific/Engineering', 
     27        'Topic :: Scientific/Engineering :: Chemistry', 
     28        'Topic :: Scientific/Engineering :: Physics', 
     29    ], 
     30    packages=packages, 
     31    package_data=package_data, 
     32    install_requires = required, 
     33    extras_require = { 
     34        'OpenCL': ["pyopencl"], 
     35        'Bumps': ["bumps"], 
     36        } 
     37 
     38    ) 
Note: See TracChangeset for help on using the changeset viewer.