Changes in / [168052c:ed82794] in sasmodels


Ignore:
Files:
14 edited

Legend:

Unmodified
Added
Removed

  • extra/pylint.rc

    r37a7252 rd4666ca  
    77# pygtk.require(). 
    88#init-hook='import os, sys; sys.path.extend([os.getcwd(),os.path.join(os.getcwd(),"extra")])' 
     9init-hook='import sys; sys.path.append('extra') 
    910 
    1011# Profiled execution. 
    … …  
    6465    star-args, 
    6566    unbalanced-tuple-unpacking, 
     67    locally-enabled, 
    6668    locally-disabled, 
    6769 
    … …  
    107109good-names=_, 
    108110    input, 
    109     i,j,k,n,x,y,z, 
     111    h,i,j,k,n,w,x,y,z, 
    110112    q,qx,qy,qz, 
    111113    dt,dx,dy,dz,id, 
    112     Iq,dIq,Qx,Qy,Qz, 
     114    Iq,dIq,Qx,Qy,Qz,Iqxy, 
    113115    p, 
    114116    ER, call_ER, VR, call_VR, 
     117    Rmax, SElength, 
    115118 
    116119# Bad variable names which should always be refused, separated by a comma 
    … …  
    335338 
    336339# Maximum number of arguments for function / method 
    337 max-args=5 
     340max-args=15 
    338341 
    339342# Argument names that match this expression will be ignored. Default to name 
    … …  
    342345 
    343346# Maximum number of locals for function / method body 
    344 max-locals=15 
     347max-locals=25 
    345348 
    346349# Maximum number of return / yield for function / method body 

  • sasmodels/bumps_model.py

    r37a7252 r190fc2b  
    1212""" 
    1313 
     14import warnings 
     15 
     16import numpy as np 
     17 
     18from .data import plot_theory 
     19from .direct_model import DataMixin 
     20 
    1421__all__ = [ 
    1522    "Model", "Experiment", 
    1623    ] 
    17  
    18 import warnings 
    19  
    20 import numpy as np 
    21  
    22 from .data import plot_theory 
    23 from .direct_model import DataMixin 
    2424 
    2525# CRUFT: old style bumps wrapper which doesn't separate data and model 

  • sasmodels/compare.py

    rd15a908 r190fc2b  
    2828 
    2929from __future__ import print_function 
     30 
     31import sys 
     32import math 
     33from os.path import basename, dirname, join as joinpath 
     34import glob 
     35import datetime 
     36import traceback 
     37 
     38import numpy as np 
     39 
     40from . import core 
     41from . import kerneldll 
     42from . import generate 
     43from .data import plot_theory, empty_data1D, empty_data2D 
     44from .direct_model import DirectModel 
     45from .convert import revert_model, constrain_new_to_old 
    3046 
    3147USAGE = """ 
    … …  
    8096           + USAGE) 
    8197 
    82  
    83  
    84 import sys 
    85 import math 
    86 from os.path import basename, dirname, join as joinpath 
    87 import glob 
    88 import datetime 
    89 import traceback 
    90  
    91 import numpy as np 
    92  
     98kerneldll.ALLOW_SINGLE_PRECISION_DLLS = True 
     99 
     100# List of available models 
    93101ROOT = dirname(__file__) 
    94 sys.path.insert(0, ROOT)  # Make sure sasmodels is first on the path 
    95  
    96  
    97 from . import core 
    98 from . import kerneldll 
    99 from . import generate 
    100 from .data import plot_theory, empty_data1D, empty_data2D 
    101 from .direct_model import DirectModel 
    102 from .convert import revert_model, constrain_new_to_old 
    103 kerneldll.ALLOW_SINGLE_PRECISION_DLLS = True 
    104  
    105 # List of available models 
    106102MODELS = [basename(f)[:-3] 
    107103          for f in sorted(glob.glob(joinpath(ROOT, "models", "[a-zA-Z]*.py")))] 
    … …  
    117113        return dt.total_seconds() 
    118114 
     115 
     116class push_seed(object): 
     117    """ 
     118    Set the seed value for the random number generator. 
     119 
     120    When used in a with statement, the random number generator state is 
     121    restored after the with statement is complete. 
     122 
     123    :Parameters: 
     124 
     125    *seed* : int or array_like, optional 
     126        Seed for RandomState 
     127 
     128    :Example: 
     129 
     130    Seed can be used directly to set the seed:: 
     131 
     132        >>> from numpy.random import randint 
     133        >>> push_seed(24) 
     134        <...push_seed object at...> 
     135        >>> print(randint(0,1000000,3)) 
     136        [242082    899 211136] 
     137 
     138    Seed can also be used in a with statement, which sets the random 
     139    number generator state for the enclosed computations and restores 
     140    it to the previous state on completion:: 
     141 
     142        >>> with push_seed(24): 
     143        ...    print(randint(0,1000000,3)) 
     144        [242082    899 211136] 
     145 
     146    Using nested contexts, we can demonstrate that state is indeed 
     147    restored after the block completes:: 
     148 
     149        >>> with push_seed(24): 
     150        ...    print(randint(0,1000000)) 
     151        ...    with push_seed(24): 
     152        ...        print(randint(0,1000000,3)) 
     153        ...    print(randint(0,1000000)) 
     154        242082 
     155        [242082    899 211136] 
     156        899 
     157 
     158    The restore step is protected against exceptions in the block:: 
     159 
     160        >>> with push_seed(24): 
     161        ...    print(randint(0,1000000)) 
     162        ...    try: 
     163        ...        with push_seed(24): 
     164        ...            print(randint(0,1000000,3)) 
     165        ...            raise Exception() 
     166        ...    except: 
     167        ...        print("Exception raised") 
     168        ...    print(randint(0,1000000)) 
     169        242082 
     170        [242082    899 211136] 
     171        Exception raised 
     172        899 
     173    """ 
     174    def __init__(self, seed=None): 
     175        self._state = np.random.get_state() 
     176        np.random.seed(seed) 
     177 
     178    def __enter__(self): 
     179        return None 
     180 
     181    def __exit__(self, *args): 
     182        np.random.set_state(self._state) 
    119183 
    120184def tic(): 
    … …  
    179243        return [0, (2*v if v > 0 else 1)] 
    180244 
     245 
    181246def _randomize_one(p, v): 
    182247    """ 
    … …  
    188253        return np.random.uniform(*parameter_range(p, v)) 
    189254 
     255 
    190256def randomize_pars(pars, seed=None): 
    191257    """ 
    … …  
    197263    :func:`constrain_pars` needs to be called afterward.. 
    198264    """ 
    199     np.random.seed(seed) 
    200     # Note: the sort guarantees order `of calls to random number generator 
    201     pars = dict((p, _randomize_one(p, v)) 
    202                 for p, v in sorted(pars.items())) 
     265    with push_seed(seed): 
     266        # Note: the sort guarantees order `of calls to random number generator 
     267        pars = dict((p, _randomize_one(p, v)) 
     268                    for p, v in sorted(pars.items())) 
    203269    return pars 
    204270 
    … …  
    463529    if Nbase > 0: 
    464530        if Ncomp > 0: plt.subplot(131) 
    465         plot_theory(data, base_value, view=view, plot_data=False, limits=limits) 
     531        plot_theory(data, base_value, view=view, use_data=False, limits=limits) 
    466532        plt.title("%s t=%.1f ms"%(base.engine, base_time)) 
    467533        #cbar_title = "log I" 
    468534    if Ncomp > 0: 
    469535        if Nbase > 0: plt.subplot(132) 
    470         plot_theory(data, comp_value, view=view, plot_data=False, limits=limits) 
     536        plot_theory(data, comp_value, view=view, use_data=False, limits=limits) 
    471537        plt.title("%s t=%.1f ms"%(comp.engine, comp_time)) 
    472538        #cbar_title = "log I" 
    … …  
    478544            err, errstr, errview = abs(relerr), "rel err", "log" 
    479545        #err,errstr = base/comp,"ratio" 
    480         plot_theory(data, None, resid=err, view=errview, plot_data=False) 
     546        plot_theory(data, None, resid=err, view=errview, use_data=False) 
    481547        plt.title("max %s = %.3g"%(errstr, max(abs(err)))) 
    482548        #cbar_title = errstr if errview=="linear" else "log "+errstr 

  • sasmodels/compare_many.py

    rd15a908 r4f2478e  
    261261 
    262262if __name__ == "__main__": 
     263    #from .compare import push_seed 
     264    #with push_seed(1): main() 
    263265    main() 

  • sasmodels/core.py

    rd15a908 r190fc2b  
    22Core model handling routines. 
    33""" 
    4 __all__ = [ 
    5     "list_models", "load_model_definition", "precompile_dll", 
    6     "load_model", "make_kernel", "call_kernel", "call_ER", "call_VR", 
    7     ] 
    84 
    95from os.path import basename, dirname, join as joinpath 
    … …  
    2420    HAVE_OPENCL = False 
    2521 
     22__all__ = [ 
     23    "list_models", "load_model_definition", "precompile_dll", 
     24    "load_model", "make_kernel", "call_kernel", "call_ER", "call_VR", 
     25] 
    2626 
    2727def list_models(): 

  • sasmodels/data.py

    rd15a908 r69ec80f  
    9090        self.x, self.y, self.dx, self.dy = x, y, dx, dy 
    9191        self.dxl = None 
     92        self.filename = None 
     93        self.qmin = x.min() if x is not None else np.NaN 
     94        self.qmax = x.max() if x is not None else np.NaN 
     95        self.mask = np.isnan(y) if y is not None else None 
     96        self._xaxis, self._xunit = "x", "" 
     97        self._yaxis, self._yunit = "y", "" 
    9298 
    9399    def xaxis(self, label, unit): 
    … …  
    108114 
    109115class Data2D(object): 
    110     def __init__(self): 
     116    def __init__(self, x=None, y=None, z=None, dx=None, dy=None, dz=None): 
     117        self.qx_data, self.dqx_data = x, dx 
     118        self.qy_data, self.dqy_data = y, dy 
     119        self.data, self.err_data = z, dz 
     120        self.mask = ~np.isnan(z) if z is not None else None 
     121        self.q_data = np.sqrt(x**2 + y**2) 
     122        self.qmin = 1e-16 
     123        self.qmax = np.inf 
    111124        self.detector = [] 
    112125        self.source = Source() 
     126        self.Q_unit = "1/A" 
     127        self.I_unit = "1/cm" 
     128        self.xaxis("Q_x", "A^{-1}") 
     129        self.yaxis("Q_y", "A^{-1}") 
     130        self.zaxis("Intensity", r"\text{cm}^{-1}") 
     131        self._xaxis, self._xunit = "x", "" 
     132        self._yaxis, self._yunit = "y", "" 
     133        self._zaxis, self._zunit = "z", "" 
     134        self.x_bins, self.y_bins = None, None 
    113135 
    114136    def xaxis(self, label, unit): 
    … …  
    139161 
    140162class Detector(object): 
     163    """ 
     164    Detector attributes. 
     165    """ 
     166    def __init__(self, pixel_size=(None, None), distance=None): 
     167        self.pixel_size = Vector(*pixel_size) 
     168        self.distance = distance 
     169 
     170class Source(object): 
     171    """ 
     172    Beam attributes. 
     173    """ 
    141174    def __init__(self): 
    142         self.pixel_size = Vector() 
    143  
    144 class Source(object): 
    145     pass 
     175        self.wavelength = np.NaN 
     176        self.wavelength_unit = "A" 
    146177 
    147178 
    … …  
    158189    data = Data1D(q, Iq, dx=resolution * q, dy=dIq) 
    159190    data.filename = "fake data" 
    160     data.qmin, data.qmax = q.min(), q.max() 
    161     data.mask = np.zeros(len(q), dtype='bool') 
    162191    return data 
    163192 
    … …  
    173202    if qy is None: 
    174203        qy = qx 
     204    # 5% dQ/Q resolution 
    175205    Qx, Qy = np.meshgrid(qx, qy) 
    176206    Qx, Qy = Qx.flatten(), Qy.flatten() 
    177207    Iq = 100 * np.ones_like(Qx) 
    178208    dIq = np.sqrt(Iq) 
    179     mask = np.ones(len(Iq), dtype='bool') 
    180  
    181     data = Data2D() 
    182     data.filename = "fake data" 
    183     data.qx_data = Qx 
    184     data.qy_data = Qy 
    185     data.data = Iq 
    186     data.err_data = dIq 
    187     data.mask = mask 
    188     data.qmin = 1e-16 
    189     data.qmax = np.inf 
    190  
    191     # 5% dQ/Q resolution 
    192209    if resolution != 0: 
    193210        # https://www.ncnr.nist.gov/staff/hammouda/distance_learning/chapter_15.pdf 
    … …  
    197214        # radial (which instead it should be inverse). 
    198215        Q = np.sqrt(Qx**2 + Qy**2) 
    199         data.dqx_data = resolution * Q 
    200         data.dqy_data = resolution * Q 
     216        dqx = resolution * Q 
     217        dqy = resolution * Q 
    201218    else: 
    202         data.dqx_data = data.dqy_data = None 
    203  
    204     detector = Detector() 
    205     detector.pixel_size.x = 5 # mm 
    206     detector.pixel_size.y = 5 # mm 
    207     detector.distance = 4 # m 
    208     data.detector.append(detector) 
     219        dqx = dqy = None 
     220 
     221    data = Data2D(x=Qx, y=Qy, z=Iq, dx=dqx, dy=dqy, dz=dIq) 
    209222    data.x_bins = qx 
    210223    data.y_bins = qy 
     224    data.filename = "fake data" 
     225 
     226    # pixel_size in mm, distance in m 
     227    detector = Detector(pixel_size=(5, 5), distance=4) 
     228    data.detector.append(detector) 
    211229    data.source.wavelength = 5 # angstroms 
    212230    data.source.wavelength_unit = "A" 
    213     data.Q_unit = "1/A" 
    214     data.I_unit = "1/cm" 
    215     data.q_data = np.sqrt(Qx ** 2 + Qy ** 2) 
    216     data.xaxis("Q_x", "A^{-1}") 
    217     data.yaxis("Q_y", "A^{-1}") 
    218     data.zaxis("Intensity", r"\text{cm}^{-1}") 
    219231    return data 
    220232 
    … …  
    228240    # do not repeat. 
    229241    if hasattr(data, 'lam'): 
    230         _plot_result_sesans(data, None, None, plot_data=True, limits=limits) 
     242        _plot_result_sesans(data, None, None, use_data=True, limits=limits) 
    231243    elif hasattr(data, 'qx_data'): 
    232         _plot_result2D(data, None, None, view, plot_data=True, limits=limits) 
     244        _plot_result2D(data, None, None, view, use_data=True, limits=limits) 
    233245    else: 
    234         _plot_result1D(data, None, None, view, plot_data=True, limits=limits) 
     246        _plot_result1D(data, None, None, view, use_data=True, limits=limits) 
    235247 
    236248 
    237249def plot_theory(data, theory, resid=None, view='log', 
    238                 plot_data=True, limits=None): 
     250                use_data=True, limits=None): 
    239251    if hasattr(data, 'lam'): 
    240         _plot_result_sesans(data, theory, resid, plot_data=True, limits=limits) 
     252        _plot_result_sesans(data, theory, resid, use_data=True, limits=limits) 
    241253    elif hasattr(data, 'qx_data'): 
    242         _plot_result2D(data, theory, resid, view, plot_data, limits=limits) 
     254        _plot_result2D(data, theory, resid, view, use_data, limits=limits) 
    243255    else: 
    244         _plot_result1D(data, theory, resid, view, plot_data, limits=limits) 
     256        _plot_result1D(data, theory, resid, view, use_data, limits=limits) 
    245257 
    246258 
    … …  
    251263        except: 
    252264            traceback.print_exc() 
    253             pass 
    254265 
    255266    return wrapper 
    … …  
    257268 
    258269@protect 
    259 def _plot_result1D(data, theory, resid, view, plot_data, limits=None): 
     270def _plot_result1D(data, theory, resid, view, use_data, limits=None): 
    260271    """ 
    261272    Plot the data and residuals for 1D data. 
    … …  
    264275    from numpy.ma import masked_array, masked 
    265276 
    266     plot_theory = theory is not None 
    267     plot_resid = resid is not None 
    268  
    269     if data.y is None: 
    270         plot_data = False 
     277    use_data = use_data and data.y is not None 
     278    use_theory = theory is not None 
     279    use_resid = resid is not None 
     280    num_plots = (use_data or use_theory) + use_resid 
     281 
    271282    scale = data.x**4 if view == 'q4' else 1.0 
    272283 
    273     if plot_data or plot_theory: 
    274         if plot_resid: 
    275             plt.subplot(121) 
    276  
     284    if use_data or use_theory: 
    277285        #print(vmin, vmax) 
    278286        all_positive = True 
    279287        some_present = False 
    280         if plot_data: 
     288        if use_data: 
    281289            mdata = masked_array(data.y, data.mask.copy()) 
    282290            mdata[~np.isfinite(mdata)] = masked 
    … …  
    288296 
    289297 
    290         if plot_theory: 
     298        if use_theory: 
    291299            mtheory = masked_array(theory, data.mask.copy()) 
    292300            mtheory[~np.isfinite(mtheory)] = masked 
    … …  
    299307        if limits is not None: 
    300308            plt.ylim(*limits) 
     309 
     310        if num_plots > 1: 
     311            plt.subplot(1, num_plots, 1) 
    301312        plt.xscale('linear' if not some_present else view) 
    302313        plt.yscale('linear' 
    … …  
    306317        plt.ylabel('$I(q)$') 
    307318 
    308     if plot_resid: 
    309         if plot_data or plot_theory: 
    310             plt.subplot(122) 
    311  
     319    if use_resid: 
    312320        mresid = masked_array(resid, data.mask.copy()) 
    313321        mresid[~np.isfinite(mresid)] = masked 
    314322        some_present = (mresid.count() > 0) 
     323 
     324        if num_plots > 1: 
     325            plt.subplot(1, num_plots, (use_data or use_theory) + 1) 
    315326        plt.plot(data.x/10, mresid, '-') 
    316327        plt.xlabel("$q$/nm$^{-1}$") 
    … …  
    320331 
    321332@protect 
    322 def _plot_result_sesans(data, theory, resid, plot_data, limits=None): 
     333def _plot_result_sesans(data, theory, resid, use_data, limits=None): 
    323334    import matplotlib.pyplot as plt 
    324     if data.y is None: 
    325         plot_data = False 
    326     plot_theory = theory is not None 
    327     plot_resid = resid is not None 
    328  
    329     if plot_data or plot_theory: 
    330         if plot_resid: 
    331             plt.subplot(121) 
    332         if plot_data: 
     335    use_data = use_data and data.y is not None 
     336    use_theory = theory is not None 
     337    use_resid = resid is not None 
     338    num_plots = (use_data or use_theory) + use_resid 
     339 
     340    if use_data or use_theory: 
     341        if num_plots > 1: 
     342            plt.subplot(1, num_plots, 1) 
     343        if use_data: 
    333344            plt.errorbar(data.x, data.y, yerr=data.dy) 
    334345        if theory is not None: 
    … …  
    340351 
    341352    if resid is not None: 
    342         if plot_data or plot_theory: 
    343             plt.subplot(122) 
    344  
     353        if num_plots > 1: 
     354            plt.subplot(1, num_plots, (use_data or use_theory) + 1) 
    345355        plt.plot(data.x, resid, 'x') 
    346356        plt.xlabel('spin echo length (nm)') 
    … …  
    349359 
    350360@protect 
    351 def _plot_result2D(data, theory, resid, view, plot_data, limits=None): 
     361def _plot_result2D(data, theory, resid, view, use_data, limits=None): 
    352362    """ 
    353363    Plot the data and residuals for 2D data. 
    354364    """ 
    355365    import matplotlib.pyplot as plt 
    356     if data.data is None: 
    357         plot_data = False 
    358     plot_theory = theory is not None 
    359     plot_resid = resid is not None 
     366    use_data = use_data and data.data is not None 
     367    use_theory = theory is not None 
     368    use_resid = resid is not None 
     369    num_plots = use_data + use_theory + use_resid 
    360370 
    361371    # Put theory and data on a common colormap scale 
    362     if limits is None: 
    363         vmin, vmax = np.inf, -np.inf 
    364         if plot_data: 
    365             target = data.data[~data.mask] 
    366             datamin = target[target>0].min() if view == 'log' else target.min() 
    367             datamax = target.max() 
    368             vmin = min(vmin, datamin) 
    369             vmax = max(vmax, datamax) 
    370         if plot_theory: 
    371             theorymin = theory[theory>0].min() if view=='log' else theory.min() 
    372             theorymax = theory.max() 
    373             vmin = min(vmin, theorymin) 
    374             vmax = max(vmax, theorymax) 
    375     else: 
     372    vmin, vmax = np.inf, -np.inf 
     373    if use_data: 
     374        target = data.data[~data.mask] 
     375        datamin = target[target > 0].min() if view == 'log' else target.min() 
     376        datamax = target.max() 
     377        vmin = min(vmin, datamin) 
     378        vmax = max(vmax, datamax) 
     379    if use_theory: 
     380        theorymin = theory[theory > 0].min() if view == 'log' else theory.min() 
     381        theorymax = theory.max() 
     382        vmin = min(vmin, theorymin) 
     383        vmax = max(vmax, theorymax) 
     384 
     385    # Override data limits from the caller 
     386    if limits is not None: 
    376387        vmin, vmax = limits 
    377388 
    378     if plot_data: 
    379         if plot_theory and plot_resid: 
    380             plt.subplot(131) 
    381         elif plot_theory or plot_resid: 
    382             plt.subplot(121) 
     389    # Plot data 
     390    if use_data: 
     391        if num_plots > 1: 
     392            plt.subplot(1, num_plots, 1) 
    383393        _plot_2d_signal(data, target, view=view, vmin=vmin, vmax=vmax) 
    384394        plt.title('data') 
    … …  
    386396        h.set_label('$I(q)$') 
    387397 
    388     if plot_theory: 
    389         if plot_data and plot_resid: 
    390             plt.subplot(132) 
    391         elif plot_data: 
    392             plt.subplot(122) 
    393         elif plot_resid: 
    394             plt.subplot(121) 
     398    # plot theory 
     399    if use_theory: 
     400        if num_plots > 1: 
     401            plt.subplot(1, num_plots, use_data+1) 
    395402        _plot_2d_signal(data, theory, view=view, vmin=vmin, vmax=vmax) 
    396403        plt.title('theory') 
    … …  
    400407                    else '$I(q)$') 
    401408 
    402     #if plot_data or plot_theory: 
    403     #    plt.colorbar() 
    404  
    405     if plot_resid: 
    406         if plot_data and plot_theory: 
    407             plt.subplot(133) 
    408         elif plot_data or plot_theory: 
    409             plt.subplot(122) 
     409    # plot resid 
     410    if use_resid: 
     411        if num_plots > 1: 
     412            plt.subplot(1, num_plots, use_data+use_theory+1) 
    410413        _plot_2d_signal(data, resid, view='linear') 
    411414        plt.title('residuals') 

  • sasmodels/generate.py

    re66c9f9 r190fc2b  
    197197# TODO: identify model files which have changed since loading and reload them. 
    198198 
    199 __all__ = ["make", "doc", "sources", "convert_type"] 
    200  
    201199import sys 
    202200from os.path import abspath, dirname, join as joinpath, exists, basename, \ 
    … …  
    206204 
    207205import numpy as np 
     206 
     207__all__ = ["make", "doc", "sources", "convert_type"] 
     208 
    208209C_KERNEL_TEMPLATE_PATH = joinpath(dirname(__file__), 'kernel_template.c') 
    209210 
    … …  
    216217    F128 = None 
    217218 
    218  
    219219# Scale and background, which are parameters common to every form factor 
    220220COMMON_PARAMETERS = [ 
    … …  
    222222    ["background", "1/cm", 0, [0, np.inf], "", "Source background"], 
    223223    ] 
    224  
    225224 
    226225# Conversion from units defined in the parameter table for each model 
    … …  
    266265 
    267266def format_units(units): 
     267    """ 
     268    Convert units into ReStructured Text format. 
     269    """ 
    268270    return "string" if isinstance(units, list) else RST_UNITS.get(units, units) 
    269271 
    … …  
    335337    """ 
    336338    Convert code from double precision to the desired type. 
     339 
     340    Floating point constants are tagged with 'f' for single precision or 'L' 
     341    for long double precision. 
    337342    """ 
    338343    if dtype == F16: 
    … …  
    350355 
    351356def _convert_type(source, type_name, constant_flag): 
     357    """ 
     358    Replace 'double' with *type_name* in *source*, tagging floating point 
     359    constants with *constant_flag*. 
     360    """ 
    352361    # Convert double keyword to float/long double/half. 
    353362    # Accept an 'n' # parameter for vector # values, where n is 2, 4, 8 or 16. 
    … …  
    625634                            %re.escape(string.punctuation)) 
    626635def _convert_section_titles_to_boldface(lines): 
     636    """ 
     637    Do the actual work of identifying and converting section headings. 
     638    """ 
    627639    prior = None 
    628640    for line in lines: 
    … …  
    642654        yield prior 
    643655 
    644 def convert_section_titles_to_boldface(string): 
    645     return "\n".join(_convert_section_titles_to_boldface(string.split('\n'))) 
     656def convert_section_titles_to_boldface(s): 
     657    """ 
     658    Use explicit bold-face rather than section headings so that the table of 
     659    contents is not polluted with section names from the model documentation. 
     660 
     661    Sections are identified as the title line followed by a line of punctuation 
     662    at least as long as the title line. 
     663    """ 
     664    return "\n".join(_convert_section_titles_to_boldface(s.split('\n'))) 
    646665 
    647666def doc(kernel_module): 
    … …  
    666685 
    667686def demo_time(): 
     687    """ 
     688    Show how long it takes to process a model. 
     689    """ 
    668690    from .models import cylinder 
    669691    import datetime 
    … …  
    674696 
    675697def main(): 
     698    """ 
     699    Program which prints the source produced by the model. 
     700    """ 
    676701    if len(sys.argv) <= 1: 
    677702        print("usage: python -m sasmodels.generate modelname") 

  • sasmodels/models/lamellarCailleHG.py

    rd18f8a8 rd4666ca  
    9292 
    9393parameters = [ 
    94               #   [ "name", "units", default, [lower, upper], "type", 
    95               #     "description" ], 
    96               [ "tail_length", "Ang", 10, [0, inf], "volume", 
    97                 "Tail thickness" ], 
    98               [ "head_length", "Ang", 2, [0, inf], "volume", 
    99                 "head thickness" ], 
    100               [ "Nlayers", "", 30, [0, inf], "", 
    101                 "Number of layers" ], 
    102               [ "spacing", "Ang", 40., [0.0,inf], "volume", 
    103                 "d-spacing of Caille S(Q)" ], 
    104               [ "Caille_parameter", "", 0.001, [0.0,0.8], "", 
    105                 "Caille parameter" ], 
    106               [ "sld", "1e-6/Ang^2", 0.4, [-inf,inf], "", 
    107                 "Tail scattering length density" ], 
    108               [ "head_sld", "1e-6/Ang^2", 2.0, [-inf,inf], "", 
    109                 "Head scattering length density" ], 
    110               [ "solvent_sld", "1e-6/Ang^2", 6, [-inf,inf], "", 
    111                 "Solvent scattering length density" ], 
     94    #   [ "name", "units", default, [lower, upper], "type", 
     95    #     "description" ], 
     96    ["tail_length", "Ang", 10, [0, inf], "volume", 
     97     "Tail thickness"], 
     98    ["head_length", "Ang", 2, [0, inf], "volume", 
     99     "head thickness"], 
     100    ["Nlayers", "", 30, [0, inf], "", 
     101     "Number of layers"], 
     102    ["spacing", "Ang", 40., [0.0, inf], "volume", 
     103     "d-spacing of Caille S(Q)"], 
     104    ["Caille_parameter", "", 0.001, [0.0, 0.8], "", 
     105     "Caille parameter"], 
     106    ["sld", "1e-6/Ang^2", 0.4, [-inf, inf], "", 
     107     "Tail scattering length density"], 
     108    ["head_sld", "1e-6/Ang^2", 2.0, [-inf, inf], "", 
     109     "Head scattering length density"], 
     110    ["solvent_sld", "1e-6/Ang^2", 6, [-inf, inf], "", 
     111     "Solvent scattering length density"], 
    112112    ] 
    113113 
    114 source = [ "lamellarCailleHG_kernel.c"] 
     114source = ["lamellarCailleHG_kernel.c"] 
    115115 
    116116# No volume normalization despite having a volume parameter 
    … …  
    128128 
    129129demo = dict( 
    130             scale=1, background=0, 
    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            ) 
     130    scale=1, background=0, 
     131    Nlayers=20, spacing=200., Caille_parameter=0.05, 
     132    tail_length=15, head_length=10, 
     133    #sld=-1, head_sld=4.0, solvent_sld=6.0, 
     134    sld=-1, head_sld=4.1, solvent_sld=6.0, 
     135    tail_length_pd=0.1, tail_length_pd_n=20, 
     136    head_length_pd=0.05, head_length_pd_n=30, 
     137    spacing_pd=0.2, spacing_pd_n=40, 
     138    ) 
    140139 
    141140oldname = 'LamellarPSHGModel' 
    142 oldpars = dict(tail_length='deltaT',head_length='deltaH',Nlayers='n_plates',Caille_parameter='caille', sld='sld_tail', head_sld='sld_head',solvent_sld='sld_solvent') 
     141oldpars = dict( 
     142    tail_length='deltaT', head_length='deltaH', Nlayers='n_plates', 
     143    Caille_parameter='caille', sld='sld_tail', head_sld='sld_head', 
     144    solvent_sld='sld_solvent') 

  • sasmodels/models/lib/sph_j1c.c

    r9c461c7 ra3e78c3  
    11/** 
    2 * Spherical Bessel function j1(x)/x 
     2* Spherical Bessel function 3*j1(x)/x 
    33* 
    44* Used for low q to avoid cancellation error. 
    … …  
    88* requires the first 3 terms.  Double precision requires the 4th term. 
    99* The fifth term is not needed, and is commented out. 
     10* Taylor expansion: 
     11*      1.0 + q2*(-3./30. + q2*(3./840.))+ q2*(-3./45360. + q2*(3./3991680.)))) 
     12* Expression returned from Herbie (herbie.uwpise.org/demo): 
     13*      const double t = ((1. + 3.*q2*q2/5600.) - q2/20.); 
     14*      return t*t; 
    1015*/ 
    1116 
    … …  
    1823    SINCOS(q, sin_q, cos_q); 
    1924 
    20     const double bessel = (q < 1.e-1) ? 
    21         1.0 + q2*(-3./30. + q2*(3./840. + q2*(-3./45360.))) // + q2*(3./3991680.))) 
     25    const double bessel = (q < 0.384038453352533) 
     26        ? (1.0 + q2*(-3./30. + q2*(3./840.))) 
    2227        : 3.0*(sin_q/q - cos_q)/q2; 
    2328 
    2429    return bessel; 
     30 
     31 /* 
     32    // Code to test various expressions 
     33    if (sizeof(q2) > 4) { 
     34        return 3.0*(sin_q/q - cos_q)/q2; 
     35    } else if (q < 0.384038453352533) { 
     36        //const double t = ((1. + 3.*q2*q2/5600.) - q2/20.); return t*t; 
     37        return 1.0 + q2*q2*(3./840.) - q2*(3./30.); 
     38        //return 1.0 + q2*(-3./30. + q2*(3./840.)); 
     39        //return 1.0 + q2*(-3./30. + q2*(3./840. + q2*(-3./45360.))); 
     40        //return 1.0 + q2*(-3./30. + q2*(3./840. + q2*(-3./45360. + q2*(3./3991680.)))); 
     41    } else { 
     42        return 3.0*(sin_q/q - cos_q)/q2; 
     43    } 
     44*/ 
    2545} 

  • sasmodels/models/pearl_necklace.py

    rcf404cb r841753c  
    11r""" 
    2 This model provides the form factor for a pearl necklace composed of two  
    3 elements: *N* pearls (homogeneous spheres of radius *R*) freely jointed by *M*  
     2This model provides the form factor for a pearl necklace composed of two 
     3elements: *N* pearls (homogeneous spheres of radius *R*) freely jointed by *M* 
    44rods (like strings - with a total mass *Mw* = *M* \* *m*\ :sub:`r` + *N* \* *m*\ 
    5 :sub:`s`, and the string segment length (or edge separation) *l*  
     5:sub:`s`, and the string segment length (or edge separation) *l* 
    66(= *A* - 2\ *R*)). *A* is the center-to-center pearl separation distance. 
    77 
    … …  
    1313---------- 
    1414 
    15 The output of the scattering intensity function for the PearlNecklaceModel is  
     15The output of the scattering intensity function for the PearlNecklaceModel is 
    1616given by (Schweins, 2004) 
    1717 
    … …  
    3737    \beta(q) &= \frac{\int_{qR}^{q(A-R)}\frac{sin(t)}{t}dt}{ql} 
    3838 
    39 where the mass *m*\ :sub:`i` is (SLD\ :sub:`i` - SLD\ :sub:`solvent`) \*  
     39where the mass *m*\ :sub:`i` is (SLD\ :sub:`i` - SLD\ :sub:`solvent`) \* 
    4040(volume of the *N* pearls/rods). *V* is the total volume of the necklace. 
    4141 
    42 The 2D scattering intensity is the same as $P(q)$ above, regardless of the  
     42The 2D scattering intensity is the same as $P(q)$ above, regardless of the 
    4343orientation of the *q* vector. 
    4444 
    … …  
    5454REFERENCE 
    5555 
    56 R Schweins and K Huber, *Particle Scattering Factor of Pearl Necklace Chains*,  
     56R Schweins and K Huber, *Particle Scattering Factor of Pearl Necklace Chains*, 
    5757*Macromol. Symp.* 211 (2004) 25-42 2004 
    5858""" 
    … …  
    7979 
    8080#             ["name", "units", default, [lower, upper], "type","description"], 
    81 parameters = [["radius", "Angstrom", 80.0, [0, inf], "volume",  
     81parameters = [["radius", "Angstrom", 80.0, [0, inf], "volume", 
    8282               "Mean radius of the chained spheres"], 
    83               ["edge_separation", "Angstrom", 350.0, [0, inf], "volume",  
     83              ["edge_separation", "Angstrom", 350.0, [0, inf], "volume", 
    8484               "Mean separation of chained particles"], 
    85               ["string_thickness", "Angstrom", 2.5, [0, inf], "volume",  
     85              ["string_thickness", "Angstrom", 2.5, [0, inf], "volume", 
    8686               "Thickness of the chain linkage"], 
    87               ["number_of_pearls", "none", 3, [0, inf], "volume",  
     87              ["number_of_pearls", "none", 3, [0, inf], "volume", 
    8888               "Mean number of pearls in each necklace"], 
    89               ["sld", "Angstrom^2", 1.0, [-inf, inf], "",  
     89              ["sld", "Angstrom^2", 1.0, [-inf, inf], "", 
    9090               "Scattering length density of the chained spheres"], 
    91               ["string_sld", "Angstrom^2", 1.0, [-inf, inf], "",  
     91              ["string_sld", "Angstrom^2", 1.0, [-inf, inf], "", 
    9292               "Scattering length density of the chain linkage"], 
    93               ["solvent_sld", "Angstrom^2", 6.3, [-inf, inf], "",  
     93              ["solvent_sld", "Angstrom^2", 6.3, [-inf, inf], "", 
    9494               "Scattering length density of the solvent"], 
    95               ] 
     95             ] 
    9696 
    9797source = ["lib/Si.c", "pearl_necklace.c"] 
    … …  
    110110# names and the target sasview model name. 
    111111oldname = 'PearlNecklaceModel' 
    112 oldpars = dict(scale='scale',background='background',radius='radius', 
     112oldpars = dict(scale='scale', background='background', radius='radius', 
    113113               number_of_pearls='num_pearls', solvent_sld='sld_solv', 
    114114               string_thickness='thick_string', sld='sld_pearl', 

  • sasmodels/models/power_law.py

    reb69cce r841753c  
    1212    I(q) = \text{scale} \cdot q^{-\text{power}} + \text{background} 
    1313 
    14 Note the minus sign in front of the exponent. The exponent *power*  
     14Note the minus sign in front of the exponent. The exponent *power* 
    1515should therefore be entered as a **positive** number for fitting. 
    1616 
    17 Also note that unlike many other models, *scale* in this model  
    18 is NOT explicitly related to a volume fraction. Be careful if  
     17Also note that unlike many other models, *scale* in this model 
     18is NOT explicitly related to a volume fraction. Be careful if 
    1919combining this model with other models. 
    2020 
    … …  
    3131from numpy import inf, sqrt 
    3232 
    33 name =  "power_law" 
     33name = "power_law" 
    3434title = "Simple power law with a flat background" 
    3535 
    36 description = """\ 
    37         Evaluates the function 
    38         I(q) = scale * q^(-power) + background 
    39         NB: enter power as a positive number! 
    40         """ 
     36description = """ 
     37    Evaluates the function 
     38    I(q) = scale * q^(-power) + background 
     39    NB: enter power as a positive number! 
     40    """ 
    4141category = "shape-independent" 
    4242 
    … …  
    4545 
    4646# NB: Scale and Background are implicit parameters on every model 
    47 def Iq(q,power): 
     47def Iq(q, power): 
     48    # pylint: disable=missing-docstring 
    4849    inten = (q**-power) 
    4950    return inten 
    … …  
    5152 
    5253def Iqxy(qx, qy, *args): 
     54    # pylint: disable=missing-docstring 
    5355    return Iq(sqrt(qx ** 2 + qy ** 2), *args) 
    5456Iqxy.vectorized = True # Iqxy accepts an array of qx, qy values 
    … …  
    6466 
    6567tests = [ 
    66         [ {'scale': 1.0, 'power': 4.0, 'background' : 0.0}, [0.0106939, 0.469418], [7.64644e+07, 20.5949]] 
    67         ] 
     68    [{'scale': 1.0, 'power': 4.0, 'background' : 0.0}, 
     69     [0.0106939, 0.469418], [7.64644e+07, 20.5949]], 
     70    ] 

  • sasmodels/resolution.py

    rfdc538a r190fc2b  
    55""" 
    66from __future__ import division 
     7 
     8from scipy.special import erf 
     9from numpy import sqrt, log, log10 
     10import numpy as np 
    711 
    812__all__ = ["Resolution", "Perfect1D", "Pinhole1D", "Slit1D", 
    … …  
    1115           "interpolate", "linear_extrapolation", "geometric_extrapolation", 
    1216           ] 
    13  
    14 from scipy.special import erf 
    15 from numpy import sqrt, log, log10 
    16 import numpy as np 
    1717 
    1818MINIMUM_RESOLUTION = 1e-8 

  • sasmodels/resolution2d.py

    r9404dd3 r841753c  
    22#This software was developed by the University of Tennessee as part of the 
    33#Distributed Data Analysis of Neutron Scattering Experiments (DANSE) 
    4 #project funded by the US National Science Foundation.  
     4#project funded by the US National Science Foundation. 
    55#See the license text in license.txt 
    66""" 
    … …  
    1919## Defaults 
    2020NR = {'xhigh':10, 'high':5, 'med':5, 'low':3} 
    21 NPHI ={'xhigh':20, 'high':12, 'med':6, 'low':4} 
     21NPHI = {'xhigh':20, 'high':12, 'med':6, 'low':4} 
    2222 
    2323class Pinhole2D(Resolution): 
    … …  
    2525    Gaussian Q smearing class for SAS 2d data 
    2626    """ 
    27       
     27 
    2828    def __init__(self, data=None, index=None, 
    2929                 nsigma=NSIGMA, accuracy='Low', coords='polar'): 
    3030        """ 
    3131        Assumption: equally spaced bins in dq_r, dq_phi space. 
    32          
     32 
    3333        :param data: 2d data used to set the smearing parameters 
    3434        :param index: 1d array with len(data) to define the range 
    … …  
    4242        ## number of bins in r axis for over-sampling 
    4343        self.nr = NR[accuracy.lower()] 
    44         ## number of bins in phi axis for over-sampling  
     44        ## number of bins in phi axis for over-sampling 
    4545        self.nphi = NPHI[accuracy.lower()] 
    4646        ## maximum nsigmas 
    47         self.nsigma= nsigma 
     47        self.nsigma = nsigma 
    4848        self.coords = coords 
    4949        self._init_data(data, index) 
    … …  
    8585    def _calc_res(self): 
    8686        """ 
    87         Over sampling of r_nbins times phi_nbins, calculate Gaussian weights,  
     87        Over sampling of r_nbins times phi_nbins, calculate Gaussian weights, 
    8888        then find smeared intensity 
    89         """     
     89        """ 
    9090        nr, nphi = self.nr, self.nphi 
    9191        # Total number of bins = # of bins 
    … …  
    143143        if self.coords == 'polar': 
    144144            q_r = sqrt(qx**2 + qy**2) 
    145             qx_res = ( (dqx*cos(dphi) + q_r) * cos(-q_phi) + 
     145            qx_res = ((dqx*cos(dphi) + q_r) * cos(-q_phi) + 
    146146                           dqy*sin(dphi) * sin(-q_phi)) 
    147147            qy_res = (-(dqx*cos(dphi) + q_r) * sin(-q_phi) + 
    … …  
    167167        else: 
    168168            return theory 
    169  
    170 """ 
    171 if __name__ == '__main__': 
    172     ## Test w/ 2D linear function 
    173     x = 0.001*np.arange(1, 11) 
    174     dx = np.ones(len(x))*0.0003 
    175     y = 0.001*np.arange(1, 11) 
    176     dy = np.ones(len(x))*0.001 
    177     z = np.ones(10) 
    178     dz = sqrt(z) 
    179  
    180     from sas.dataloader import Data2D 
    181     #for i in range(10): print(i, 0.001 + i*0.008/9.0) 
    182     #for i in range(100): print(i, int(math.floor( (i/ (100/9.0)) ))) 
    183     out = Data2D() 
    184     out.data = z 
    185     out.qx_data = x 
    186     out.qy_data = y 
    187     out.dqx_data = dx 
    188     out.dqy_data = dy 
    189     out.q_data = sqrt(dx * dx + dy * dy) 
    190     index = np.ones(len(x), dtype = bool) 
    191     out.mask = index 
    192     from sas.models.LineModel import LineModel 
    193     model = LineModel() 
    194     model.setParam("A", 0) 
    195  
    196     smear = Smearer2D(out, model, index) 
    197     #smear.set_accuracy('Xhigh') 
    198     value = smear.get_value() 
    199     ## All data are ones, so the smeared should also be ones. 
    200     print("Data length =", len(value)) 
    201     print(" 2D linear function, I = 0 + 1*qy") 
    202     text = " Gaussian weighted averaging on a 2D linear function will " 
    203     text += "provides the results same as without the averaging." 
    204     print(text) 
    205     print("qx_data", "qy_data", "I_nonsmear", "I_smeared") 
    206     for ind in range(len(value)): 
    207         print(x[ind], y[ind], model.evalDistribution([x, y])[ind], value[ind]) 
    208  
    209  
    210 if __name__ == '__main__': 
    211     ## Another Test w/ constant function 
    212     x = 0.001*np.arange(1,11) 
    213     dx = np.ones(len(x))*0.001 
    214     y = 0.001*np.arange(1,11) 
    215     dy = np.ones(len(x))*0.001 
    216     z = np.ones(10) 
    217     dz = sqrt(z) 
    218  
    219     from DataLoader import Data2D 
    220     #for i in range(10): print(i, 0.001 + i*0.008/9.0) 
    221     #for i in range(100): print(i, int(math.floor( (i/ (100/9.0)) ))) 
    222     out = Data2D() 
    223     out.data = z 
    224     out.qx_data = x 
    225     out.qy_data = y 
    226     out.dqx_data = dx 
    227     out.dqy_data = dy 
    228     index = np.ones(len(x), dtype = bool) 
    229     out.mask = index 
    230     from sas.models.Constant import Constant 
    231     model = Constant() 
    232  
    233     value = Smearer2D(out,model,index).get_value() 
    234     ## All data are ones, so the smeared values should also be ones. 
    235     print("Data length =",len(value), ", Data=",value) 
    236 """ 

  • sasmodels/sesans.py

    r384d114 r190fc2b  
    4343    *I* [cm$^{-1}$] is the value of the SANS model at *q* 
    4444    """ 
    45     G = np.zeros(len(SElength), 'd') 
    46     for i in range(len(SElength)): 
    47         integr = besselj(0, q*SElength[i])*Iq*q 
    48         G[i] = np.sum(integr) 
     45    G = np.zeros_like(SElength, 'd') 
     46    for i, SElength_i in enumerate(SElength): 
     47        integral = besselj(0, q*SElength_i)*Iq*q 
     48        G[i] = np.sum(integral) 
    4949 
    5050    # [m^-1] step size in q, needed for integration 
    51     dq=(q[1]-q[0])*1e10 
     51    dq = (q[1]-q[0])*1e10 
    5252 
    5353    # integration step, convert q into [m**-1] and 2 pi circle integration 
Note: See TracChangeset for help on using the changeset viewer.