Changes in / [0ce5710:62cf915] in sasmodels


Ignore:
Files:
7 deleted
49 edited

Legend:

Unmodified
Added
Removed

  • doc/developer/index.rst

    rb85be2d r56fc97a  
    33*************************** 
    44 
    5 .. toctree:: 
    6    :numbered: 4 
    7    :maxdepth: 4 
    85 
    9    calculator.rst 

  • doc/genapi.py

    ra5b8477 r3d5c6f8  
    5959    #('alignment', 'GPU data alignment [unused]'), 
    6060    ('bumps_model', 'Bumps interface'), 
    61     ('compare', 'Compare models on different compute engines'), 
    6261    ('convert', 'Sasview to sasmodel converter'), 
    6362    ('core', 'Model access'), 
    64     ('data', 'Data layout and plotting routines'), 
    65     ('details', 'Parameter packing for kernel calls'), 
    6663    ('direct_model', 'Simple interface'), 
    6764    ('exception', 'Annotate exceptions'), 
    68     #('frozendict', 'Freeze a dictionary to make it immutable'), 
    6965    ('generate', 'Model parser'), 
    70     ('kernel', 'Evaluator type definitions'), 
    7166    ('kernelcl', 'OpenCL model evaluator'), 
    7267    ('kerneldll', 'Ctypes model evaluator'), 
    7368    ('kernelpy', 'Python model evaluator'), 
    74     ('list_pars', 'Identify all parameters in all models'), 
    75     ('mixture', 'Mixture model evaluator'), 
    7669    ('model_test', 'Unit test support'), 
    77     ('modelinfo', 'Parameter and model definitions'), 
    78     ('product', 'Product model evaluator'), 
    7970    ('resolution', '1-D resolution functions'), 
    8071    ('resolution2d', '2-D resolution functions'), 
    8172    ('sasview_model', 'Sasview interface'), 
    82     ('sesans', 'SESANS calculation routines'), 
     73    ('sesans', 'SESANS model evaluator'), 
     74    ('weights', 'Distribution functions'), 
    8375    #('transition', 'Model stepper for automatic model selection'), 
    84     ('weights', 'Distribution functions'), 
    8576] 
    8677package='sasmodels' 

  • doc/genmodel.py

    ra5b8477 r89f4163  
    1 from __future__ import print_function 
    2  
    31import sys, os, math, re 
    42import numpy as np 
    53import matplotlib.pyplot as plt 
     4import pylab 
    65sys.path.insert(0, os.path.abspath('..')) 
    76from sasmodels import generate, core 
    … …  
    98from sasmodels.data import empty_data1D, empty_data2D 
    109 
    11 try: 
    12     from typing import Dict, Any 
    13 except ImportError: 
    14     pass 
    15 else: 
    16     from matplotlib.axes import Axes 
    17     from sasmodels.kernel import KernelModel 
    18     from sasmodels.modelinfo import ModelInfo 
     10 
     11# Convert ../sasmodels/models/name.py to name 
     12model_name = os.path.basename(sys.argv[1])[:-3] 
     13model_info = core.load_model_info(model_name) 
     14model = core.build_model(model_info) 
     15 
     16# Load the doc string from the module definition file and store it in rst 
     17docstr = generate.make_doc(model_info) 
     18 
     19 
     20# Calculate 1D curve for default parameters 
     21pars = dict((p.name, p.default) for p in model_info['parameters']) 
     22 
     23# Plotting ranges and options 
     24opts = { 
     25    'xscale'    : 'log', 
     26    'yscale'    : 'log' if not model_info['structure_factor'] else 'linear', 
     27    'zscale'    : 'log' if not model_info['structure_factor'] else 'linear', 
     28    'q_min'     : 0.001, 
     29    'q_max'     : 1.0, 
     30    'nq'        : 1000, 
     31    'nq2d'      : 1000, 
     32    'vmin'      : 1e-3,  # floor for the 2D data results 
     33    'qx_max'    : 0.5, 
     34    #'colormap'  : 'gist_ncar', 
     35    'colormap'  : 'nipy_spectral', 
     36    #'colormap'  : 'jet', 
     37} 
     38 
    1939 
    2040def plot_1d(model, opts, ax): 
    21     # type: (KernelModel, Dict[str, Any], Axes) -> None 
    22     """ 
    23     Create a 1-D image. 
    24     """ 
    2541    q_min, q_max, nq = opts['q_min'], opts['q_max'], opts['nq'] 
    2642    q_min = math.log10(q_min) 
    … …  
    3147    Iq1D = calculator() 
    3248 
    33     ax.plot(q, Iq1D, color='blue', lw=2, label=model.info.name) 
     49    ax.plot(q, Iq1D, color='blue', lw=2, label=model_info['name']) 
    3450    ax.set_xlabel(r'$Q \/(\AA^{-1})$') 
    3551    ax.set_ylabel(r'$I(Q) \/(\mathrm{cm}^{-1})$') 
    … …  
    3955 
    4056def plot_2d(model, opts, ax): 
    41     # type: (KernelModel, Dict[str, Any], Axes) -> None 
    42     """ 
    43     Create a 2-D image. 
    44     """ 
    4557    qx_max, nq2d = opts['qx_max'], opts['nq2d'] 
    46     q = np.linspace(-qx_max, qx_max, nq2d) # type: np.ndarray 
     58    q = np.linspace(-qx_max, qx_max, nq2d) 
    4759    data2d = empty_data2D(q, resolution=0.0) 
    4860    calculator = DirectModel(data2d, model) 
    … …  
    5264        Iq2D = np.log(np.clip(Iq2D, opts['vmin'], np.inf)) 
    5365    ax.imshow(Iq2D, interpolation='nearest', aspect=1, origin='lower', 
    54               extent=[-qx_max, qx_max, -qx_max, qx_max], cmap=opts['colormap']) 
     66        extent=[-qx_max, qx_max, -qx_max, qx_max], cmap=opts['colormap']) 
    5567    ax.set_xlabel(r'$Q_x \/(\AA^{-1})$') 
    5668    ax.set_ylabel(r'$Q_y \/(\AA^{-1})$') 
    5769 
    58 def figfile(model_info): 
    59     # type: (ModelInfo) -> str 
    60     return model_info.id + '_autogenfig.png' 
     70# Generate image  
     71fig_height = 3.0 # in 
     72fig_left = 0.6 # in 
     73fig_right = 0.5 # in 
     74fig_top = 0.6*0.25 # in 
     75fig_bottom = 0.6*0.75 
     76if model_info['has_2d']: 
     77    plot_height = fig_height - (fig_top+fig_bottom) 
     78    plot_width = plot_height 
     79    fig_width = 2*(plot_width + fig_left + fig_right) 
     80    aspect = (fig_width, fig_height) 
     81    ratio = aspect[0]/aspect[1] 
     82    ax_left = fig_left/fig_width 
     83    ax_bottom = fig_bottom/fig_height 
     84    ax_height = plot_height/fig_height 
     85    ax_width = ax_height/ratio # square axes 
     86    fig = plt.figure(figsize=aspect) 
     87    ax2d = fig.add_axes([0.5+ax_left, ax_bottom, ax_width, ax_height]) 
     88    plot_2d(model, opts, ax2d) 
     89    ax1d = fig.add_axes([ax_left, ax_bottom, ax_width, ax_height]) 
     90    plot_1d(model, opts, ax1d) 
     91    #ax.set_aspect('square') 
     92else: 
     93    plot_height = fig_height - (fig_top+fig_bottom) 
     94    plot_width = (1+np.sqrt(5))/2*fig_height 
     95    fig_width = plot_width + fig_left + fig_right 
     96    ax_left = fig_left/fig_width 
     97    ax_bottom = fig_bottom/fig_height 
     98    ax_width = plot_width/fig_width 
     99    ax_height = plot_height/fig_height 
     100    aspect = (fig_width, fig_height) 
     101    fig = plt.figure(figsize=aspect) 
     102    ax1d = fig.add_axes([ax_left, ax_bottom, ax_width, ax_height]) 
     103    plot_1d(model, opts, ax1d) 
    61104 
    62 def make_figure(model_info, opts): 
    63     # type: (ModelInfo, Dict[str, Any]) -> None 
    64     """ 
    65     Generate the figure file to include in the docs. 
    66     """ 
    67     model = core.build_model(model_info) 
     105# Save image in model/img 
     106figname = model_name + '_autogenfig.png' 
     107filename = os.path.join('model', 'img', figname) 
     108plt.savefig(filename, bbox_inches='tight') 
     109#print "figure saved in",filename 
    68110 
    69     fig_height = 3.0 # in 
    70     fig_left = 0.6 # in 
    71     fig_right = 0.5 # in 
    72     fig_top = 0.6*0.25 # in 
    73     fig_bottom = 0.6*0.75 
    74     if model_info.parameters.has_2d: 
    75         plot_height = fig_height - (fig_top+fig_bottom) 
    76         plot_width = plot_height 
    77         fig_width = 2*(plot_width + fig_left + fig_right) 
    78         aspect = (fig_width, fig_height) 
    79         ratio = aspect[0]/aspect[1] 
    80         ax_left = fig_left/fig_width 
    81         ax_bottom = fig_bottom/fig_height 
    82         ax_height = plot_height/fig_height 
    83         ax_width = ax_height/ratio # square axes 
    84         fig = plt.figure(figsize=aspect) 
    85         ax2d = fig.add_axes([0.5+ax_left, ax_bottom, ax_width, ax_height]) 
    86         plot_2d(model, opts, ax2d) 
    87         ax1d = fig.add_axes([ax_left, ax_bottom, ax_width, ax_height]) 
    88         plot_1d(model, opts, ax1d) 
    89         #ax.set_aspect('square') 
    90     else: 
    91         plot_height = fig_height - (fig_top+fig_bottom) 
    92         plot_width = (1+np.sqrt(5))/2*fig_height 
    93         fig_width = plot_width + fig_left + fig_right 
    94         ax_left = fig_left/fig_width 
    95         ax_bottom = fig_bottom/fig_height 
    96         ax_width = plot_width/fig_width 
    97         ax_height = plot_height/fig_height 
    98         aspect = (fig_width, fig_height) 
    99         fig = plt.figure(figsize=aspect) 
    100         ax1d = fig.add_axes([ax_left, ax_bottom, ax_width, ax_height]) 
    101         plot_1d(model, opts, ax1d) 
     111# Auto caption for figure 
     112captionstr = '\n' 
     113captionstr += '.. figure:: img/' + model_info['id'] + '_autogenfig.png\n' 
     114captionstr += '\n' 
     115if model_info['has_2d']: 
     116    captionstr += '    1D and 2D plots corresponding to the default parameters of the model.\n' 
     117else: 
     118    captionstr += '    1D plot corresponding to the default parameters of the model.\n' 
     119captionstr += '\n' 
    102120 
    103     # Save image in model/img 
    104     path = os.path.join('model', 'img', figfile(model_info)) 
    105     plt.savefig(path, bbox_inches='tight') 
    106     #print("figure saved in",path) 
     121# Add figure reference and caption to documentation (at end, before References) 
     122pattern = '\*\*REFERENCE' 
     123m = re.search(pattern, docstr.upper()) 
    107124 
    108 def gen_docs(model_info): 
    109     # type: (ModelInfo) -> None 
    110     """ 
    111     Generate the doc string with the figure inserted before the references. 
    112     """ 
     125if m: 
     126    docstr1 = docstr[:m.start()] 
     127    docstr2 = docstr[m.start():] 
     128    docstr = docstr1 + captionstr + docstr2 
     129else: 
     130    print '------------------------------------------------------------------' 
     131    print 'References NOT FOUND for model: ', model_info['id'] 
     132    print '------------------------------------------------------------------' 
     133    docstr = docstr + captionstr 
    113134 
    114     # Load the doc string from the module definition file and store it in rst 
    115     docstr = generate.make_doc(model_info) 
    116  
    117     # Auto caption for figure 
    118     captionstr = '\n' 
    119     captionstr += '.. figure:: img/' + figfile(model_info) + '\n' 
    120     captionstr += '\n' 
    121     if model_info.parameters.has_2d: 
    122         captionstr += '    1D and 2D plots corresponding to the default parameters of the model.\n' 
    123     else: 
    124         captionstr += '    1D plot corresponding to the default parameters of the model.\n' 
    125     captionstr += '\n' 
    126  
    127     # Add figure reference and caption to documentation (at end, before References) 
    128     pattern = '\*\*REFERENCE' 
    129     match = re.search(pattern, docstr.upper()) 
    130  
    131     if match: 
    132         docstr1 = docstr[:match.start()] 
    133         docstr2 = docstr[match.start():] 
    134         docstr = docstr1 + captionstr + docstr2 
    135     else: 
    136         print('------------------------------------------------------------------') 
    137         print('References NOT FOUND for model: ', model_info.id) 
    138         print('------------------------------------------------------------------') 
    139         docstr += captionstr 
    140  
    141     open(sys.argv[2],'w').write(docstr) 
    142  
    143 def process_model(path): 
    144     # type: (str) -> None 
    145     """ 
    146     Generate doc file and image file for the given model definition file. 
    147     """ 
    148  
    149     # Load the model file 
    150     model_name = os.path.basename(path)[:-3] 
    151     model_info = core.load_model_info(model_name) 
    152  
    153     # Plotting ranges and options 
    154     opts = { 
    155         'xscale'    : 'log', 
    156         'yscale'    : 'log' if not model_info.structure_factor else 'linear', 
    157         'zscale'    : 'log' if not model_info.structure_factor else 'linear', 
    158         'q_min'     : 0.001, 
    159         'q_max'     : 1.0, 
    160         'nq'        : 1000, 
    161         'nq2d'      : 1000, 
    162         'vmin'      : 1e-3,  # floor for the 2D data results 
    163         'qx_max'    : 0.5, 
    164         #'colormap'  : 'gist_ncar', 
    165         'colormap'  : 'nipy_spectral', 
    166         #'colormap'  : 'jet', 
    167     } 
    168  
    169     # Generate the RST file and the figure.  Order doesn't matter. 
    170     gen_docs(model_info) 
    171     make_figure(model_info, opts) 
    172  
    173 if __name__ == "__main__": 
    174     process_model(sys.argv[1]) 
     135open(sys.argv[2],'w').write(docstr) 

  • doc/gentoc.py

    ra5b8477 r5041682  
    99from sasmodels.core import load_model_info 
    1010 
    11 try: 
    12     from typing import Optional, BinaryIO, List, Dict 
    13 except ImportError: 
    14     pass 
    15 else: 
    16     from sasmodels.modelinfo import ModelInfo 
    1711 
    1812TEMPLATE="""\ 
    … …  
    3327 
    3428def _make_category(category_name, label, title, parent=None): 
    35     # type: (str, str, str, Optional[BinaryIO]) -> BinaryIO 
    3629    file = open(joinpath(MODEL_TOC_PATH, category_name+".rst"), "w") 
    3730    file.write(TEMPLATE%{'label':label, 'title':title, 'bar':'*'*len(title)}) 
    … …  
    4134 
    4235def _add_subcategory(category_name, parent): 
    43     # type: (str, BinaryIO) -> None 
    4436    parent.write("    %s.rst\n"%category_name) 
    4537 
    4638def _add_model(file, model_name): 
    47     # type: (IO[str], str) -> None 
    4839    file.write("    ../../model/%s.rst\n"%model_name) 
    4940 
    5041def _maybe_make_category(category, models, cat_files, model_toc): 
    51     # type: (str, List[str], Dict[str, BinaryIO], BinaryIO) -> None 
    5242    if category not in cat_files: 
    5343        print("Unexpected category %s containing"%category, models, file=sys.stderr) 
    … …  
    5646 
    5747def generate_toc(model_files): 
    58     # type: (List[str]) -> None 
    5948    if not model_files: 
    6049        print("gentoc needs a list of model files", file=sys.stderr) 
    6150 
    6251    # find all categories 
    63     category = {} # type: Dict[str, List[str]] 
     52    category = {} 
    6453    for item in model_files: 
    6554        # assume model is in sasmodels/models/name.py, and ignore the full path 
    … …  
    6756        if model_name.startswith('_'): continue 
    6857        model_info = load_model_info(model_name) 
    69         if model_info.category is None: 
     58        if model_info['category'] is None: 
    7059            print("Missing category for", item, file=sys.stderr) 
    7160        else: 
    72             category.setdefault(model_info.category,[]).append(model_name) 
     61            category.setdefault(model_info['category'],[]).append(model_name) 
    7362 
    7463    # Check category names 

  • sasmodels/alignment.py

    r7ae2b7f r3c56da87  
    1818to decide whether it is really required. 
    1919""" 
    20 import numpy as np  # type: ignore 
     20import numpy as np 
    2121 
    2222def align_empty(shape, dtype, alignment=128): 

  • sasmodels/bumps_model.py

    r04dc697 rea75043  
    1111 
    1212""" 
    13 from __future__ import print_function 
    14  
    15 __all__ = [ "Model", "Experiment" ] 
    16  
    17 import numpy as np  # type: ignore 
     13 
     14import warnings 
     15 
     16import numpy as np 
    1817 
    1918from .data import plot_theory 
    2019from .direct_model import DataMixin 
    2120 
    22 try: 
    23     from typing import Dict, Union, Tuple, Any 
    24     from .data import Data1D, Data2D 
    25     from .kernel import KernelModel 
    26     from .modelinfo import ModelInfo 
    27     Data = Union[Data1D, Data2D] 
    28 except ImportError: 
    29     pass 
    30  
    31 try: 
    32     # Optional import. This allows the doc builder and nosetests to run even 
    33     # when bumps is not on the path. 
    34     from bumps.names import Parameter # type: ignore 
    35 except ImportError: 
    36     pass 
     21__all__ = [ 
     22    "Model", "Experiment", 
     23    ] 
     24 
     25# CRUFT: old style bumps wrapper which doesn't separate data and model 
     26# pylint: disable=invalid-name 
     27def BumpsModel(data, model, cutoff=1e-5, **kw): 
     28    r""" 
     29    Bind a model to data, along with a polydispersity cutoff. 
     30 
     31    *data* is a :class:`data.Data1D`, :class:`data.Data2D` or 
     32    :class:`data.Sesans` object.  Use :func:`data.empty_data1D` or 
     33    :func:`data.empty_data2D` to define $q, \Delta q$ calculation 
     34    points for displaying the SANS curve when there is no measured data. 
     35 
     36    *model* is a runnable module as returned from :func:`core.load_model`. 
     37 
     38    *cutoff* is the polydispersity weight cutoff. 
     39 
     40    Any additional *key=value* pairs are model dependent parameters. 
     41 
     42    Returns an :class:`Experiment` object. 
     43 
     44    Note that the usual Bumps semantics is not fully supported, since 
     45    assigning *M.name = parameter* on the returned experiment object 
     46    does not set that parameter in the model.  Range setting will still 
     47    work as expected though. 
     48 
     49    .. deprecated:: 0.1 
     50        Use :class:`Experiment` instead. 
     51    """ 
     52    warnings.warn("Use of BumpsModel is deprecated.  Use bumps_model.Experiment instead.") 
     53 
     54    # Create the model and experiment 
     55    model = Model(model, **kw) 
     56    experiment = Experiment(data=data, model=model, cutoff=cutoff) 
     57 
     58    # Copy the model parameters up to the experiment object. 
     59    for k, v in model.parameters().items(): 
     60        setattr(experiment, k, v) 
     61    return experiment 
    3762 
    3863 
    3964def create_parameters(model_info, **kwargs): 
    40     # type: (ModelInfo, **Union[float, str, Parameter]) -> Tuple[Dict[str, Parameter], Dict[str, str]] 
    4165    """ 
    4266    Generate Bumps parameters from the model info. 
    … …  
    4771    Any additional *key=value* pairs are initial values for the parameters 
    4872    to the models.  Uninitialized parameters will use the model default 
    49     value.  The value can be a float, a bumps parameter, or in the case 
    50     of the distribution type parameter, a string. 
     73    value. 
    5174 
    5275    Returns a dictionary of *{name: Parameter}* containing the bumps 
    … …  
    5477    *{name: str}* containing the polydispersity distribution types. 
    5578    """ 
    56     pars = {}     # type: Dict[str, Parameter] 
    57     pd_types = {} # type: Dict[str, str] 
    58     for p in model_info.parameters.call_parameters: 
     79    # lazy import; this allows the doc builder and nosetests to run even 
     80    # when bumps is not on the path. 
     81    from bumps.names import Parameter 
     82 
     83    pars = {} 
     84    for p in model_info['parameters']: 
    5985        value = kwargs.pop(p.name, p.default) 
    6086        pars[p.name] = Parameter.default(value, name=p.name, limits=p.limits) 
    61         if p.polydisperse: 
    62             for part, default, limits in [ 
    63                     ('_pd', 0., pars[p.name].limits), 
    64                     ('_pd_n', 35., (0, 1000)), 
    65                     ('_pd_nsigma', 3., (0, 10)), 
    66                 ]: 
    67                 name = p.name + part 
    68                 value = kwargs.pop(name, default) 
    69                 pars[name] = Parameter.default(value, name=name, limits=limits) 
    70             name = p.name + '_pd_type' 
    71             pd_types[name] = str(kwargs.pop(name, 'gaussian')) 
     87    for name in model_info['partype']['pd-2d']: 
     88        for xpart, xdefault, xlimits in [ 
     89                ('_pd', 0., pars[name].limits), 
     90                ('_pd_n', 35., (0, 1000)), 
     91                ('_pd_nsigma', 3., (0, 10)), 
     92            ]: 
     93            xname = name + xpart 
     94            xvalue = kwargs.pop(xname, xdefault) 
     95            pars[xname] = Parameter.default(xvalue, name=xname, limits=xlimits) 
     96 
     97    pd_types = {} 
     98    for name in model_info['partype']['pd-2d']: 
     99        xname = name + '_pd_type' 
     100        xvalue = kwargs.pop(xname, 'gaussian') 
     101        pd_types[xname] = xvalue 
    72102 
    73103    if kwargs:  # args not corresponding to parameters 
    … …  
    88118    """ 
    89119    def __init__(self, model, **kwargs): 
    90         # type: (KernelModel, **Dict[str, Union[float, Parameter]]) -> None 
    91         self.sasmodel = model 
     120        self._sasmodel = model 
    92121        pars, pd_types = create_parameters(model.info, **kwargs) 
    93122        for k, v in pars.items(): 
    … …  
    99128 
    100129    def parameters(self): 
    101         # type: () -> Dict[str, Parameter] 
    102130        """ 
    103131        Return a dictionary of parameters objects for the parameters, 
    … …  
    107135 
    108136    def state(self): 
    109         # type: () -> Dict[str, Union[Parameter, str]] 
    110137        """ 
    111138        Return a dictionary of current values for all the parameters, 
    … …  
    132159    The resulting model can be used directly in a Bumps FitProblem call. 
    133160    """ 
    134     _cache = None # type: Dict[str, np.ndarray] 
    135161    def __init__(self, data, model, cutoff=1e-5): 
    136         # type: (Data, Model, float) -> None 
     162 
    137163        # remember inputs so we can inspect from outside 
    138164        self.model = model 
    139165        self.cutoff = cutoff 
    140         self._interpret_data(data, model.sasmodel) 
    141         self._cache = {} 
     166        self._interpret_data(data, model._sasmodel) 
     167        self.update() 
    142168 
    143169    def update(self): 
    144         # type: () -> None 
    145170        """ 
    146171        Call when model parameters have changed and theory needs to be 
    147172        recalculated. 
    148173        """ 
    149         self._cache.clear() 
     174        self._cache = {} 
    150175 
    151176    def numpoints(self): 
    152         # type: () -> float 
    153177        """ 
    154178        Return the number of data points 
    … …  
    157181 
    158182    def parameters(self): 
    159         # type: () -> Dict[str, Parameter] 
    160183        """ 
    161184        Return a dictionary of parameters 
    … …  
    164187 
    165188    def theory(self): 
    166         # type: () -> np.ndarray 
    167189        """ 
    168190        Return the theory corresponding to the model parameters. 
    … …  
    177199 
    178200    def residuals(self): 
    179         # type: () -> np.ndarray 
    180201        """ 
    181202        Return theory minus data normalized by uncertainty. 
    … …  
    185206 
    186207    def nllf(self): 
    187         # type: () -> float 
    188208        """ 
    189209        Return the negative log likelihood of seeing data given the model 
    … …  
    193213        delta = self.residuals() 
    194214        #if np.any(np.isnan(R)): print("NaN in residuals") 
    195         return 0.5 * np.sum(delta**2) 
     215        return 0.5 * np.sum(delta ** 2) 
    196216 
    197217    #def __call__(self): 
    … …  
    199219 
    200220    def plot(self, view='log'): 
    201         # type: (str) -> None 
    202221        """ 
    203222        Plot the data and residuals. 
    … …  
    207226 
    208227    def simulate_data(self, noise=None): 
    209         # type: (float) -> None 
    210228        """ 
    211229        Generate simulated data. 
    … …  
    215233 
    216234    def save(self, basename): 
    217         # type: (str) -> None 
    218235        """ 
    219236        Save the model parameters and data into a file. 
    … …  
    226243 
    227244    def __getstate__(self): 
    228         # type: () -> Dict[str, Any] 
    229245        # Can't pickle gpu functions, so instead make them lazy 
    230246        state = self.__dict__.copy() 
    … …  
    233249 
    234250    def __setstate__(self, state): 
    235         # type: (Dict[str, Any]) -> None 
    236251        # pylint: disable=attribute-defined-outside-init 
    237252        self.__dict__ = state 
    238  

  • sasmodels/compare.py

    r7ae2b7f rf247314  
    3434import traceback 
    3535 
    36 import numpy as np  # type: ignore 
     36import numpy as np 
    3737 
    3838from . import core 
    3939from . import kerneldll 
     40from . import product 
    4041from .data import plot_theory, empty_data1D, empty_data2D 
    4142from .direct_model import DirectModel 
    … …  
    209210    Choose a parameter range based on parameter name and initial value. 
    210211    """ 
    211     # process the polydispersity options 
    212212    if p.endswith('_pd_n'): 
    213213        return [0, 100] 
    … …  
    222222        else: 
    223223            return [-180, 180] 
     224    elif 'sld' in p: 
     225        return [-0.5, 10] 
    224226    elif p.endswith('_pd'): 
    225227        return [0, 1] 
    226     elif 'sld' in p: 
    227         return [-0.5, 10] 
    228228    elif p == 'background': 
    229229        return [0, 10] 
    230230    elif p == 'scale': 
    231231        return [0, 1e3] 
     232    elif p == 'case_num': 
     233        # RPA hack 
     234        return [0, 10] 
    232235    elif v < 0: 
     236        # Kxy parameters in rpa model can be negative 
    233237        return [2*v, -2*v] 
    234238    else: 
    … …  
    236240 
    237241 
    238 def _randomize_one(model_info, p, v): 
     242def _randomize_one(p, v): 
    239243    """ 
    240244    Randomize a single parameter. 
    … …  
    242246    if any(p.endswith(s) for s in ('_pd_n', '_pd_nsigma', '_pd_type')): 
    243247        return v 
    244  
    245     # Find the parameter definition 
    246     for par in model_info.parameters.call_parameters: 
    247         if par.name == p: 
    248             break 
    249     else: 
    250         raise ValueError("unknown parameter %r"%p) 
    251  
    252     if len(par.limits) > 2:  # choice list 
    253         return np.random.randint(len(par.limits)) 
    254  
    255     limits = par.limits 
    256     if not np.isfinite(limits).all(): 
    257         low, high = parameter_range(p, v) 
    258         limits = (max(limits[0], low), min(limits[1], high)) 
    259  
    260     return np.random.uniform(*limits) 
    261  
    262  
    263 def randomize_pars(model_info, pars, seed=None): 
     248    else: 
     249        return np.random.uniform(*parameter_range(p, v)) 
     250 
     251 
     252def randomize_pars(pars, seed=None): 
    264253    """ 
    265254    Generate random values for all of the parameters. 
    … …  
    272261    with push_seed(seed): 
    273262        # Note: the sort guarantees order `of calls to random number generator 
    274         pars = dict((p, _randomize_one(model_info, p, v)) 
     263        pars = dict((p, _randomize_one(p, v)) 
    275264                    for p, v in sorted(pars.items())) 
    276265    return pars 
    … …  
    284273    cylinder radius in this case). 
    285274    """ 
    286     name = model_info.id 
     275    name = model_info['id'] 
    287276    # if it is a product model, then just look at the form factor since 
    288277    # none of the structure factors need any constraints. 
    … …  
    305294        # Make sure phi sums to 1.0 
    306295        if pars['case_num'] < 2: 
    307             pars['Phi1'] = 0. 
    308             pars['Phi2'] = 0. 
     296            pars['Phia'] = 0. 
     297            pars['Phib'] = 0. 
    309298        elif pars['case_num'] < 5: 
    310             pars['Phi1'] = 0. 
    311         total = sum(pars['Phi'+c] for c in '1234') 
    312         for c in '1234': 
     299            pars['Phia'] = 0. 
     300        total = sum(pars['Phi'+c] for c in 'abcd') 
     301        for c in 'abcd': 
    313302            pars['Phi'+c] /= total 
    314303 
    … …  
    317306    Format the parameter list for printing. 
    318307    """ 
     308    if is2d: 
     309        exclude = lambda n: False 
     310    else: 
     311        partype = model_info['partype'] 
     312        par1d = set(partype['fixed-1d']+partype['pd-1d']) 
     313        exclude = lambda n: n not in par1d 
    319314    lines = [] 
    320     parameters = model_info.parameters 
    321     for p in parameters.user_parameters(pars, is2d): 
     315    for p in model_info['parameters']: 
     316        if exclude(p.name): continue 
    322317        fields = dict( 
    323             value=pars.get(p.id, p.default), 
    324             pd=pars.get(p.id+"_pd", 0.), 
    325             n=int(pars.get(p.id+"_pd_n", 0)), 
    326             nsigma=pars.get(p.id+"_pd_nsgima", 3.), 
    327             type=pars.get(p.id+"_pd_type", 'gaussian')) 
     318            value=pars.get(p.name, p.default), 
     319            pd=pars.get(p.name+"_pd", 0.), 
     320            n=int(pars.get(p.name+"_pd_n", 0)), 
     321            nsigma=pars.get(p.name+"_pd_nsgima", 3.), 
     322            type=pars.get(p.name+"_pd_type", 'gaussian')) 
    328323        lines.append(_format_par(p.name, **fields)) 
    329324    return "\n".join(lines) 
    … …  
    368363 
    369364    # grab the sasview model, or create it if it is a product model 
    370     if model_info.composition: 
    371         composition_type, parts = model_info.composition 
     365    if model_info['composition']: 
     366        composition_type, parts = model_info['composition'] 
    372367        if composition_type == 'product': 
    373368            from sas.models.MultiplicationModel import MultiplicationModel 
    … …  
    459454    """ 
    460455    # initialize the code so time is more accurate 
    461     if Nevals > 1: 
    462         value = calculator(**suppress_pd(pars)) 
     456    value = calculator(**suppress_pd(pars)) 
    463457    toc = tic() 
    464458    for _ in range(max(Nevals, 1)):  # make sure there is at least one eval 
    … …  
    667661    """ 
    668662    # Get the default values for the parameters 
    669     pars = {} 
    670     for p in model_info.parameters.call_parameters: 
    671         parts = [('', p.default)] 
    672         if p.polydisperse: 
    673             parts.append(('_pd', 0.0)) 
    674             parts.append(('_pd_n', 0)) 
    675             parts.append(('_pd_nsigma', 3.0)) 
    676             parts.append(('_pd_type', "gaussian")) 
    677         for ext, val in parts: 
    678             if p.length > 1: 
    679                 dict(("%s%d%s"%(p.id,k,ext), val) for k in range(p.length)) 
    680             else: 
    681                 pars[p.id+ext] = val 
     663    pars = dict((p.name, p.default) for p in model_info['parameters']) 
     664 
     665    # Fill in default values for the polydispersity parameters 
     666    for p in model_info['parameters']: 
     667        if p.type in ('volume', 'orientation'): 
     668            pars[p.name+'_pd'] = 0.0 
     669            pars[p.name+'_pd_n'] = 0 
     670            pars[p.name+'_pd_nsigma'] = 3.0 
     671            pars[p.name+'_pd_type'] = "gaussian" 
    682672 
    683673    # Plug in values given in demo 
    684674    if use_demo: 
    685         pars.update(model_info.demo) 
     675        pars.update(model_info['demo']) 
    686676    return pars 
    687677 
    … …  
    710700    try: 
    711701        model_info = core.load_model_info(name) 
    712     except ImportError as exc: 
     702    except ImportError, exc: 
    713703        print(str(exc)) 
    714704        print("Could not find model; use one of:\n    " + models) 
    … …  
    784774 
    785775    if len(engines) == 0: 
    786         engines.extend(['single', 'double']) 
     776        engines.extend(['single', 'sasview']) 
    787777    elif len(engines) == 1: 
    788         if engines[0][0] != 'double': 
    789             engines.append('double') 
     778        if engines[0][0] != 'sasview': 
     779            engines.append('sasview') 
    790780        else: 
    791781            engines.append('single') 
    … …  
    817807    #pars.update(set_pars)  # set value before random to control range 
    818808    if opts['seed'] > -1: 
    819         pars = randomize_pars(model_info, pars, seed=opts['seed']) 
     809        pars = randomize_pars(pars, seed=opts['seed']) 
    820810        print("Randomize using -random=%i"%opts['seed']) 
    821811    if opts['mono']: 
    … …  
    860850    Explore the model using the Bumps GUI. 
    861851    """ 
    862     import wx  # type: ignore 
    863     from bumps.names import FitProblem  # type: ignore 
    864     from bumps.gui.app_frame import AppFrame  # type: ignore 
     852    import wx 
     853    from bumps.names import FitProblem 
     854    from bumps.gui.app_frame import AppFrame 
    865855 
    866856    problem = FitProblem(Explore(opts)) 
    … …  
    883873    """ 
    884874    def __init__(self, opts): 
    885         from bumps.cli import config_matplotlib  # type: ignore 
     875        from bumps.cli import config_matplotlib 
    886876        from . import bumps_model 
    887877        config_matplotlib() 
    … …  
    889879        model_info = opts['def'] 
    890880        pars, pd_types = bumps_model.create_parameters(model_info, **opts['pars']) 
    891         # Initialize parameter ranges, fixing the 2D parameters for 1D data. 
    892881        if not opts['is2d']: 
    893             for p in model_info.parameters.user_parameters(is2d=False): 
    894                 for ext in ['', '_pd', '_pd_n', '_pd_nsigma']: 
    895                     k = p.name+ext 
    896                     v = pars.get(k, None) 
    897                     if v is not None: 
    898                         v.range(*parameter_range(k, v.value)) 
     882            active = [base + ext 
     883                      for base in model_info['partype']['pd-1d'] 
     884                      for ext in ['', '_pd', '_pd_n', '_pd_nsigma']] 
     885            active.extend(model_info['partype']['fixed-1d']) 
     886            for k in active: 
     887                v = pars[k] 
     888                v.range(*parameter_range(k, v.value)) 
    899889        else: 
    900890            for k, v in pars.items(): 

  • sasmodels/compare_many.py

    r7ae2b7f rce346b6  
    1717import traceback 
    1818 
    19 import numpy as np  # type: ignore 
     19import numpy as np 
    2020 
    2121from . import core 
     22from . import generate 
    2223from .compare import (MODELS, randomize_pars, suppress_pd, make_data, 
    2324                      make_engine, get_pars, columnize, 
    … …  
    108109 
    109110    if is_2d: 
    110         if not model_info['parameters'].has_2d: 
     111        if not model_info['has_2d']: 
    111112            print(',"1-D only"') 
    112113            return 
    … …  
    124125        try: 
    125126            result = fn(**pars) 
    126         except Exception: 
     127        except KeyboardInterrupt: 
     128            raise 
     129        except: 
    127130            traceback.print_exc() 
    128131            print("when comparing %s for %d"%(name, seed)) 
    … …  
    243246        base = sys.argv[5] 
    244247        comp = sys.argv[6] if len(sys.argv) > 6 else "sasview" 
    245     except Exception: 
     248    except: 
    246249        traceback.print_exc() 
    247250        print_usage() 

  • sasmodels/convert.json

    r6e7ff6d rec45c4f  
    602602    "RPAModel",  
    603603    { 
    604       "N1": "Na", "N2": "Nb", "N3": "Nc", "N4": "Nd", 
    605       "L1": "La", "L2": "Lb", "L3": "Lc", "L4": "Ld", 
    606       "v1": "va", "v2": "vb", "v3": "vc", "v4": "vd", 
    607       "b1": "ba", "b2": "bb", "b3": "bc", "b4": "bd", 
    608       "Phi1": "Phia", "Phi2": "Phib", "Phi3": "Phic", "Phi4": "Phid", 
    609604      "case_num": "lcase_n" 
    610605    } 
    … …  
    625620    } 
    626621  ],  
    627   "_spherepy": [ 
    628     "SphereModel", 
    629     { 
    630       "sld": "sldSph", 
    631       "radius": "radius", 
    632       "sld_solvent": "sldSolv" 
    633     } 
    634   ], 
    635622  "spherical_sld": [ 
    636623    "SphereSLDModel",  
    637624    { 
    638       "n": "n_shells", 
    639625      "radius_core": "rad_core",  
    640626      "sld_solvent": "sld_solv" 

  • sasmodels/convert.py

    r7ae2b7f rf247314  
    6262    # but only if we haven't already byteified it 
    6363    if isinstance(data, dict) and not ignore_dicts: 
    64         return dict((_byteify(key, ignore_dicts=True), 
    65                      _byteify(value, ignore_dicts=True)) 
    66                     for key, value in data.items()) 
     64        return { 
     65            _byteify(key, ignore_dicts=True): _byteify(value, ignore_dicts=True) 
     66            for key, value in data.iteritems() 
     67        } 
    6768    # if it's anything else, return it in its original form 
    6869    return data 
    … …  
    152153def revert_name(model_info): 
    153154    _read_conversion_table() 
    154     oldname, oldpars = CONVERSION_TABLE.get(model_info.id, [None, {}]) 
     155    oldname, oldpars = CONVERSION_TABLE.get(model_info['id'], [None, {}]) 
    155156    return oldname 
    156157 
    157158def _get_old_pars(model_info): 
    158159    _read_conversion_table() 
    159     oldname, oldpars = CONVERSION_TABLE.get(model_info.id, [None, {}]) 
     160    oldname, oldpars = CONVERSION_TABLE.get(model_info['id'], [None, {}]) 
    160161    return oldpars 
    161162 
    … …  
    164165    Convert model from new style parameter names to old style. 
    165166    """ 
    166     if model_info.composition is not None: 
    167         composition_type, parts = model_info.composition 
     167    if model_info['composition'] is not None: 
     168        composition_type, parts = model_info['composition'] 
    168169        if composition_type == 'product': 
    169170            P, S = parts 
    … …  
    179180 
    180181    # Note: update compare.constrain_pars to match 
    181     name = model_info.id 
    182     if name in MODELS_WITHOUT_SCALE or model_info.structure_factor: 
     182    name = model_info['id'] 
     183    if name in MODELS_WITHOUT_SCALE or model_info['structure_factor']: 
    183184        if oldpars.pop('scale', 1.0) != 1.0: 
    184185            warnings.warn("parameter scale not used in sasview %s"%name) 
    185     if name in MODELS_WITHOUT_BACKGROUND or model_info.structure_factor: 
     186    if name in MODELS_WITHOUT_BACKGROUND or model_info['structure_factor']: 
    186187        if oldpars.pop('background', 0.0) != 0.0: 
    187188            warnings.warn("parameter background not used in sasview %s"%name) 
    … …  
    205206        elif name == 'rpa': 
    206207            # convert scattering lengths from femtometers to centimeters 
    207             for p in "L1", "L2", "L3", "L4": 
     208            for p in "La", "Lb", "Lc", "Ld": 
    208209                if p in oldpars: oldpars[p] *= 1e-13 
    209210        elif name == 'core_shell_parallelepiped': 
    … …  
    224225    Restrict parameter values to those that will match sasview. 
    225226    """ 
    226     name = model_info.id 
     227    name = model_info['id'] 
    227228    # Note: update convert.revert_model to match 
    228     if name in MODELS_WITHOUT_SCALE or model_info.structure_factor: 
     229    if name in MODELS_WITHOUT_SCALE or model_info['structure_factor']: 
    229230        pars['scale'] = 1 
    230     if name in MODELS_WITHOUT_BACKGROUND or model_info.structure_factor: 
     231    if name in MODELS_WITHOUT_BACKGROUND or model_info['structure_factor']: 
    231232        pars['background'] = 0 
    232233    # sasview multiplies background by structure factor 

  • sasmodels/core.py

    rfa5fd8d r4d76711  
    22Core model handling routines. 
    33""" 
    4 from __future__ import print_function 
    5  
    64__all__ = [ 
    7     "list_models", "load_model", "load_model_info", 
    8     "build_model", "precompile_dll", 
     5    "list_models", "load_model_info", "precompile_dll", 
     6    "build_model", "make_kernel", "call_kernel", "call_ER_VR", 
    97    ] 
    108 
    11 from os.path import basename, dirname, join as joinpath 
     9from os.path import basename, dirname, join as joinpath, splitext 
    1210from glob import glob 
    1311 
    14 import numpy as np # type: ignore 
    15  
     12import numpy as np 
     13 
     14from . import models 
     15from . import weights 
    1616from . import generate 
    17 from . import modelinfo 
    18 from . import product 
     17# TODO: remove circular references between product and core 
     18# product uses call_ER/call_VR, core uses make_product_info/ProductModel 
     19#from . import product 
    1920from . import mixture 
    2021from . import kernelpy 
    … …  
    2324    from . import kernelcl 
    2425    HAVE_OPENCL = True 
    25 except Exception: 
     26except: 
    2627    HAVE_OPENCL = False 
    2728 
    2829try: 
    29     from typing import List, Union, Optional, Any 
    30     DType = Union[None, str, np.dtype] 
    31     from .kernel import KernelModel 
    32 except ImportError: 
    33     pass 
    34  
     30    np.meshgrid([]) 
     31    meshgrid = np.meshgrid 
     32except ValueError: 
     33    # CRUFT: np.meshgrid requires multiple vectors 
     34    def meshgrid(*args): 
     35        if len(args) > 1: 
     36            return np.meshgrid(*args) 
     37        else: 
     38            return [np.asarray(v) for v in args] 
    3539 
    3640# TODO: refactor composite model support 
    … …  
    4751 
    4852def list_models(): 
    49     # type: () -> List[str] 
    5053    """ 
    5154    Return the list of available models on the model path. 
    … …  
    5760 
    5861def isstr(s): 
    59     # type: (Any) -> bool 
    6062    """ 
    6163    Return True if *s* is a string-like object. 
    6264    """ 
    6365    try: s + '' 
    64     except Exception: return False 
     66    except: return False 
    6567    return True 
    6668 
    67 def load_model(model_name, dtype=None, platform='ocl'): 
    68     # type: (str, DType, str) -> KernelModel 
     69def load_model(model_name, **kw): 
    6970    """ 
    7071    Load model info and build model. 
    71  
    72     *model_name* is the name of the model as used by :func:`load_model_info`. 
    73     Additional keyword arguments are passed directly to :func:`build_model`. 
    74     """ 
    75     return build_model(load_model_info(model_name), 
    76                        dtype=dtype, platform=platform) 
     72    """ 
     73    return build_model(load_model_info(model_name), **kw) 
    7774 
    7875 
    7976def load_model_info(model_name): 
    80     # type: (str) -> modelinfo.ModelInfo 
    8177    """ 
    8278    Load a model definition given the model name. 
    … …  
    9288    parts = model_name.split('*') 
    9389    if len(parts) > 1: 
     90        from . import product 
     91        # Note: currently have circular reference 
    9492        if len(parts) > 2: 
    9593            raise ValueError("use P*S to apply structure factor S to model P") 
    … …  
    9896 
    9997    kernel_module = generate.load_kernel_module(model_name) 
    100     return modelinfo.make_model_info(kernel_module) 
     98    return generate.make_model_info(kernel_module) 
    10199 
    102100 
    103101def build_model(model_info, dtype=None, platform="ocl"): 
    104     # type: (modelinfo.ModelInfo, np.dtype, str) -> KernelModel 
    105102    """ 
    106103    Prepare the model for the default execution platform. 
    … …  
    121118    otherwise it uses the default "ocl". 
    122119    """ 
    123     composition = model_info.composition 
     120    composition = model_info.get('composition', None) 
    124121    if composition is not None: 
    125122        composition_type, parts = composition 
    … …  
    134131            raise ValueError('unknown mixture type %s'%composition_type) 
    135132 
    136     # If it is a python model, return it immediately 
    137     if callable(model_info.Iq): 
    138         return kernelpy.PyModel(model_info) 
    139  
    140133    ## for debugging: 
    141134    ##  1. uncomment open().write so that the source will be saved next time 
    … …  
    144137    ##  4. rerun "python -m sasmodels.direct_model $MODELNAME" 
    145138    ##  5. uncomment open().read() so that source will be regenerated from model 
    146     # open(model_info.name+'.c','w').write(source) 
    147     # source = open(model_info.name+'.cl','r').read() 
     139    # open(model_info['name']+'.c','w').write(source) 
     140    # source = open(model_info['name']+'.cl','r').read() 
    148141    source = generate.make_source(model_info) 
    149142    if dtype is None: 
    150         dtype = generate.F32 if model_info.single else generate.F64 
     143        dtype = 'single' if model_info['single'] else 'double' 
     144    if callable(model_info.get('Iq', None)): 
     145        return kernelpy.PyModel(model_info) 
    151146    if (platform == "dll" 
    152147            or not HAVE_OPENCL 
    … …  
    157152 
    158153def precompile_dll(model_name, dtype="double"): 
    159     # type: (str, DType) -> Optional[str] 
    160154    """ 
    161155    Precompile the dll for a model. 
    … …  
    174168    source = generate.make_source(model_info) 
    175169    return kerneldll.make_dll(source, model_info, dtype=dtype) if source else None 
     170 
     171 
     172def get_weights(model_info, pars, name): 
     173    """ 
     174    Generate the distribution for parameter *name* given the parameter values 
     175    in *pars*. 
     176 
     177    Uses "name", "name_pd", "name_pd_type", "name_pd_n", "name_pd_sigma" 
     178    from the *pars* dictionary for parameter value and parameter dispersion. 
     179    """ 
     180    relative = name in model_info['partype']['pd-rel'] 
     181    limits = model_info['limits'][name] 
     182    disperser = pars.get(name+'_pd_type', 'gaussian') 
     183    value = pars.get(name, model_info['defaults'][name]) 
     184    npts = pars.get(name+'_pd_n', 0) 
     185    width = pars.get(name+'_pd', 0.0) 
     186    nsigma = pars.get(name+'_pd_nsigma', 3.0) 
     187    value, weight = weights.get_weights( 
     188        disperser, npts, width, nsigma, value, limits, relative) 
     189    return value, weight / np.sum(weight) 
     190 
     191def dispersion_mesh(pars): 
     192    """ 
     193    Create a mesh grid of dispersion parameters and weights. 
     194 
     195    Returns [p1,p2,...],w where pj is a vector of values for parameter j 
     196    and w is a vector containing the products for weights for each 
     197    parameter set in the vector. 
     198    """ 
     199    value, weight = zip(*pars) 
     200    value = [v.flatten() for v in meshgrid(*value)] 
     201    weight = np.vstack([v.flatten() for v in meshgrid(*weight)]) 
     202    weight = np.prod(weight, axis=0) 
     203    return value, weight 
     204 
     205def call_kernel(kernel, pars, cutoff=0, mono=False): 
     206    """ 
     207    Call *kernel* returned from *model.make_kernel* with parameters *pars*. 
     208 
     209    *cutoff* is the limiting value for the product of dispersion weights used 
     210    to perform the multidimensional dispersion calculation more quickly at a 
     211    slight cost to accuracy. The default value of *cutoff=0* integrates over 
     212    the entire dispersion cube.  Using *cutoff=1e-5* can be 50% faster, but 
     213    with an error of about 1%, which is usually less than the measurement 
     214    uncertainty. 
     215 
     216    *mono* is True if polydispersity should be set to none on all parameters. 
     217    """ 
     218    fixed_pars = [pars.get(name, kernel.info['defaults'][name]) 
     219                  for name in kernel.fixed_pars] 
     220    if mono: 
     221        pd_pars = [( np.array([pars[name]]), np.array([1.0]) ) 
     222                   for name in kernel.pd_pars] 
     223    else: 
     224        pd_pars = [get_weights(kernel.info, pars, name) for name in kernel.pd_pars] 
     225    return kernel(fixed_pars, pd_pars, cutoff=cutoff) 
     226 
     227def call_ER_VR(model_info, vol_pars): 
     228    """ 
     229    Return effect radius and volume ratio for the model. 
     230 
     231    *info* is either *kernel.info* for *kernel=make_kernel(model,q)* 
     232    or *model.info*. 
     233 
     234    *pars* are the parameters as expected by :func:`call_kernel`. 
     235    """ 
     236    ER = model_info.get('ER', None) 
     237    VR = model_info.get('VR', None) 
     238    value, weight = dispersion_mesh(vol_pars) 
     239 
     240    individual_radii = ER(*value) if ER else 1.0 
     241    whole, part = VR(*value) if VR else (1.0, 1.0) 
     242 
     243    effect_radius = np.sum(weight*individual_radii) / np.sum(weight) 
     244    volume_ratio = np.sum(weight*part)/np.sum(weight*whole) 
     245    return effect_radius, volume_ratio 
     246 
     247 
     248def call_ER(model_info, values): 
     249    """ 
     250    Call the model ER function using *values*. *model_info* is either 
     251    *model.info* if you have a loaded model, or *kernel.info* if you 
     252    have a model kernel prepared for evaluation. 
     253    """ 
     254    ER = model_info.get('ER', None) 
     255    if ER is None: 
     256        return 1.0 
     257    else: 
     258        vol_pars = [get_weights(model_info, values, name) 
     259                    for name in model_info['partype']['volume']] 
     260        value, weight = dispersion_mesh(vol_pars) 
     261        individual_radii = ER(*value) 
     262        #print(values[0].shape, weights.shape, fv.shape) 
     263        return np.sum(weight*individual_radii) / np.sum(weight) 
     264 
     265def call_VR(model_info, values): 
     266    """ 
     267    Call the model VR function using *pars*. 
     268    *info* is either *model.info* if you have a loaded model, or *kernel.info* 
     269    if you have a model kernel prepared for evaluation. 
     270    """ 
     271    VR = model_info.get('VR', None) 
     272    if VR is None: 
     273        return 1.0 
     274    else: 
     275        vol_pars = [get_weights(model_info, values, name) 
     276                    for name in model_info['partype']['volume']] 
     277        value, weight = dispersion_mesh(vol_pars) 
     278        whole, part = VR(*value) 
     279        return np.sum(weight*part)/np.sum(weight*whole) 
     280 
     281# TODO: remove call_ER, call_VR 
     282 

  • sasmodels/custom/__init__.py

    r7ae2b7f rcd0a808  
    1414try: 
    1515    # Python 3.5 and up 
    16     from importlib.util import spec_from_file_location, module_from_spec  # type: ignore 
     16    from importlib.util import spec_from_file_location, module_from_spec 
    1717    def load_module_from_path(fullname, path): 
    1818        spec = spec_from_file_location(fullname, path) 

  • sasmodels/data.py

    ra5b8477 rd6f5da6  
    3535import traceback 
    3636 
    37 import numpy as np  # type: ignore 
    38  
    39 try: 
    40     from typing import Union, Dict, List, Optional 
    41 except ImportError: 
    42     pass 
    43 else: 
    44     Data = Union["Data1D", "Data2D", "SesansData"] 
     37import numpy as np 
    4538 
    4639def load_data(filename): 
    47     # type: (str) -> Data 
    4840    """ 
    4941    Load data using a sasview loader. 
    5042    """ 
    51     from sas.sascalc.dataloader.loader import Loader  # type: ignore 
     43    from sas.sascalc.dataloader.loader import Loader 
    5244    loader = Loader() 
    5345    data = loader.load(filename) 
    … …  
    5850 
    5951def set_beam_stop(data, radius, outer=None): 
    60     # type: (Data, float, Optional[float]) -> None 
    6152    """ 
    6253    Add a beam stop of the given *radius*.  If *outer*, make an annulus. 
    6354    """ 
    64     from sas.dataloader.manipulations import Ringcut  # type: ignore 
     55    from sas.dataloader.manipulations import Ringcut 
    6556    if hasattr(data, 'qx_data'): 
    6657        data.mask = Ringcut(0, radius)(data) 
    … …  
    7465 
    7566def set_half(data, half): 
    76     # type: (Data, str) -> None 
    7767    """ 
    7868    Select half of the data, either "right" or "left". 
    7969    """ 
    80     from sas.dataloader.manipulations import Boxcut  # type: ignore 
     70    from sas.dataloader.manipulations import Boxcut 
    8171    if half == 'right': 
    8272        data.mask += \ 
    … …  
    8878 
    8979def set_top(data, cutoff): 
    90     # type: (Data, float) -> None 
    9180    """ 
    9281    Chop the top off the data, above *cutoff*. 
    9382    """ 
    94     from sas.dataloader.manipulations import Boxcut  # type: ignore 
     83    from sas.dataloader.manipulations import Boxcut 
    9584    data.mask += \ 
    9685        Boxcut(x_min=-np.inf, x_max=np.inf, y_min=-np.inf, y_max=cutoff)(data) 
    … …  
    125114    """ 
    126115    def __init__(self, x=None, y=None, dx=None, dy=None): 
    127         # type: (Optional[np.ndarray], Optional[np.ndarray], Optional[np.ndarray], Optional[np.ndarray]) -> None 
    128116        self.x, self.y, self.dx, self.dy = x, y, dx, dy 
    129117        self.dxl = None 
    … …  
    139127 
    140128    def xaxis(self, label, unit): 
    141         # type: (str, str) -> None 
    142129        """ 
    143130        set the x axis label and unit 
    … …  
    147134 
    148135    def yaxis(self, label, unit): 
    149         # type: (str, str) -> None 
    150136        """ 
    151137        set the y axis label and unit 
    … …  
    154140        self._yunit = unit 
    155141 
    156 class SesansData(Data1D): 
    157     def __init__(self, **kw): 
    158         Data1D.__init__(self, **kw) 
    159         self.lam = None # type: Optional[np.ndarray] 
     142 
    160143 
    161144class Data2D(object): 
    … …  
    192175    """ 
    193176    def __init__(self, x=None, y=None, z=None, dx=None, dy=None, dz=None): 
    194         # type: (Optional[np.ndarray], Optional[np.ndarray], Optional[np.ndarray], Optional[np.ndarray], Optional[np.ndarray], Optional[np.ndarray]) -> None 
    195177        self.qx_data, self.dqx_data = x, dx 
    196178        self.qy_data, self.dqy_data = y, dy 
    … …  
    215197 
    216198    def xaxis(self, label, unit): 
    217         # type: (str, str) -> None 
    218199        """ 
    219200        set the x axis label and unit 
    … …  
    223204 
    224205    def yaxis(self, label, unit): 
    225         # type: (str, str) -> None 
    226206        """ 
    227207        set the y axis label and unit 
    … …  
    231211 
    232212    def zaxis(self, label, unit): 
    233         # type: (str, str) -> None 
    234213        """ 
    235214        set the y axis label and unit 
    … …  
    244223    """ 
    245224    def __init__(self, x=None, y=None, z=None): 
    246         # type: (float, float, Optional[float]) -> None 
    247225        self.x, self.y, self.z = x, y, z 
    248226 
    … …  
    252230    """ 
    253231    def __init__(self, pixel_size=(None, None), distance=None): 
    254         # type: (Tuple[float, float], float) -> None 
    255232        self.pixel_size = Vector(*pixel_size) 
    256233        self.distance = distance 
    … …  
    261238    """ 
    262239    def __init__(self): 
    263         # type: () -> None 
    264240        self.wavelength = np.NaN 
    265241        self.wavelength_unit = "A" 
    … …  
    267243 
    268244def empty_data1D(q, resolution=0.0): 
    269     # type: (np.ndarray, float) -> Data1D 
    270245    """ 
    271246    Create empty 1D data using the given *q* as the x value. 
    … …  
    284259 
    285260def empty_data2D(qx, qy=None, resolution=0.0): 
    286     # type: (np.ndarray, Optional[np.ndarray], float) -> Data2D 
    287261    """ 
    288262    Create empty 2D data using the given mesh. 
    … …  
    298272    Qx, Qy = np.meshgrid(qx, qy) 
    299273    Qx, Qy = Qx.flatten(), Qy.flatten() 
    300     Iq = 100 * np.ones_like(Qx)  # type: np.ndarray 
     274    Iq = 100 * np.ones_like(Qx) 
    301275    dIq = np.sqrt(Iq) 
    302276    if resolution != 0: 
    … …  
    326300 
    327301def plot_data(data, view='log', limits=None): 
    328     # type: (Data, str, Optional[Tuple[float, float]]) -> None 
    329302    """ 
    330303    Plot data loaded by the sasview loader. 
    … …  
    350323def plot_theory(data, theory, resid=None, view='log', 
    351324                use_data=True, limits=None, Iq_calc=None): 
    352     # type: (Data, Optional[np.ndarray], Optional[np.ndarray], str, bool, Optional[Tuple[float,float]], Optional[np.ndarray]) -> None 
    353325    """ 
    354326    Plot theory calculation. 
    … …  
    365337    *limits* sets the intensity limits on the plot; if None then the limits 
    366338    are inferred from the data. 
    367  
    368     *Iq_calc* is the raw theory values without resolution smearing 
    369339    """ 
    370340    if hasattr(data, 'lam'): 
    … …  
    378348 
    379349def protect(fn): 
    380     # type: (Callable) -> Callable 
    381350    """ 
    382351    Decorator to wrap calls in an exception trapper which prints the 
    … …  
    389358        try: 
    390359            return fn(*args, **kw) 
    391         except Exception: 
     360        except KeyboardInterrupt: 
     361            raise 
     362        except: 
    392363            traceback.print_exc() 
    393364 
    … …  
    398369def _plot_result1D(data, theory, resid, view, use_data, 
    399370                   limits=None, Iq_calc=None): 
    400     # type: (Data1D, Optional[np.ndarray], Optional[np.ndarray], str, bool, Optional[Tuple[float, float]], Optional[np.ndarray]) -> None 
    401371    """ 
    402372    Plot the data and residuals for 1D data. 
    403373    """ 
    404     import matplotlib.pyplot as plt  # type: ignore 
    405     from numpy.ma import masked_array, masked  # type: ignore 
     374    import matplotlib.pyplot as plt 
     375    from numpy.ma import masked_array, masked 
    406376 
    407377    use_data = use_data and data.y is not None 
    … …  
    476446@protect 
    477447def _plot_result_sesans(data, theory, resid, use_data, limits=None): 
    478     # type: (SesansData, Optional[np.ndarray], Optional[np.ndarray], bool, Optional[Tuple[float, float]]) -> None 
    479448    """ 
    480449    Plot SESANS results. 
    481450    """ 
    482     import matplotlib.pyplot as plt  # type: ignore 
     451    import matplotlib.pyplot as plt 
    483452    use_data = use_data and data.y is not None 
    484453    use_theory = theory is not None 
    … …  
    487456 
    488457    if use_data or use_theory: 
    489         is_tof = (data.lam != data.lam[0]).any() 
     458        is_tof = np.any(data.lam!=data.lam[0]) 
    490459        if num_plots > 1: 
    491460            plt.subplot(1, num_plots, 1) 
    492461        if use_data: 
    493462            if is_tof: 
    494                 plt.errorbar(data.x, np.log(data.y)/(data.lam*data.lam), 
    495                              yerr=data.dy/data.y/(data.lam*data.lam)) 
     463                plt.errorbar(data.x, np.log(data.y)/(data.lam*data.lam), yerr=data.dy/data.y/(data.lam*data.lam)) 
    496464            else: 
    497465                plt.errorbar(data.x, data.y, yerr=data.dy) 
    … …  
    521489@protect 
    522490def _plot_result2D(data, theory, resid, view, use_data, limits=None): 
    523     # type: (Data2D, Optional[np.ndarray], Optional[np.ndarray], str, bool, Optional[Tuple[float,float]]) -> None 
    524491    """ 
    525492    Plot the data and residuals for 2D data. 
    526493    """ 
    527     import matplotlib.pyplot as plt  # type: ignore 
     494    import matplotlib.pyplot as plt 
    528495    use_data = use_data and data.data is not None 
    529496    use_theory = theory is not None 
    … …  
    533500    # Put theory and data on a common colormap scale 
    534501    vmin, vmax = np.inf, -np.inf 
    535     target = None # type: Optional[np.ndarray] 
    536502    if use_data: 
    537503        target = data.data[~data.mask] 
    … …  
    582548@protect 
    583549def _plot_2d_signal(data, signal, vmin=None, vmax=None, view='log'): 
    584     # type: (Data2D, np.ndarray, Optional[float], Optional[float], str) -> Tuple[float, float] 
    585550    """ 
    586551    Plot the target value for the data.  This could be the data itself, 
    … …  
    589554    *scale* can be 'log' for log scale data, or 'linear'. 
    590555    """ 
    591     import matplotlib.pyplot as plt  # type: ignore 
    592     from numpy.ma import masked_array  # type: ignore 
     556    import matplotlib.pyplot as plt 
     557    from numpy.ma import masked_array 
    593558 
    594559    image = np.zeros_like(data.qx_data) 
    … …  
    624589 
    625590def demo(): 
    626     # type: () -> None 
    627591    """ 
    628592    Load and plot a SAS dataset. 
    … …  
    631595    set_beam_stop(data, 0.004) 
    632596    plot_data(data) 
    633     import matplotlib.pyplot as plt  # type: ignore 
    634     plt.show() 
     597    import matplotlib.pyplot as plt; plt.show() 
    635598 
    636599 

  • sasmodels/direct_model.py

    ra5b8477 r4d76711  
    2323from __future__ import print_function 
    2424 
    25 import numpy as np  # type: ignore 
    26  
    27 # TODO: fix sesans module 
    28 from . import sesans  # type: ignore 
    29 from . import weights 
     25import numpy as np 
     26 
     27from .core import call_kernel, call_ER_VR 
     28from . import sesans 
    3029from . import resolution 
    3130from . import resolution2d 
    32 from . import details 
    33  
    34 try: 
    35     from typing import Optional, Dict, Tuple 
    36 except ImportError: 
    37     pass 
    38 else: 
    39     from .data import Data 
    40     from .kernel import Kernel, KernelModel 
    41     from .modelinfo import Parameter, ParameterSet 
    42  
    43 def call_kernel(kernel, pars, cutoff=0., mono=False): 
    44     # type: (Kernel, ParameterSet, float, bool) -> np.ndarray 
    45     """ 
    46     Call *kernel* returned from *model.make_kernel* with parameters *pars*. 
    47  
    48     *cutoff* is the limiting value for the product of dispersion weights used 
    49     to perform the multidimensional dispersion calculation more quickly at a 
    50     slight cost to accuracy. The default value of *cutoff=0* integrates over 
    51     the entire dispersion cube.  Using *cutoff=1e-5* can be 50% faster, but 
    52     with an error of about 1%, which is usually less than the measurement 
    53     uncertainty. 
    54  
    55     *mono* is True if polydispersity should be set to none on all parameters. 
    56     """ 
    57     parameters = kernel.info.parameters 
    58     if mono: 
    59         active = lambda name: False 
    60     elif kernel.dim == '1d': 
    61         active = lambda name: name in parameters.pd_1d 
    62     elif kernel.dim == '2d': 
    63         active = lambda name: name in parameters.pd_2d 
    64     else: 
    65         active = lambda name: True 
    66  
    67     vw_pairs = [(get_weights(p, pars) if active(p.name) 
    68                  else ([pars.get(p.name, p.default)], [])) 
    69                 for p in parameters.call_parameters] 
    70  
    71     call_details, weights, values = details.build_details(kernel, vw_pairs) 
    72     return kernel(call_details, weights, values, cutoff) 
    73  
    74 def get_weights(parameter, values): 
    75     # type: (Parameter, Dict[str, float]) -> Tuple[np.ndarray, np.ndarray] 
    76     """ 
    77     Generate the distribution for parameter *name* given the parameter values 
    78     in *pars*. 
    79  
    80     Uses "name", "name_pd", "name_pd_type", "name_pd_n", "name_pd_sigma" 
    81     from the *pars* dictionary for parameter value and parameter dispersion. 
    82     """ 
    83     value = float(values.get(parameter.name, parameter.default)) 
    84     relative = parameter.relative_pd 
    85     limits = parameter.limits 
    86     disperser = values.get(parameter.name+'_pd_type', 'gaussian') 
    87     npts = values.get(parameter.name+'_pd_n', 0) 
    88     width = values.get(parameter.name+'_pd', 0.0) 
    89     nsigma = values.get(parameter.name+'_pd_nsigma', 3.0) 
    90     if npts == 0 or width == 0: 
    91         return [value], [] 
    92     value, weight = weights.get_weights( 
    93         disperser, npts, width, nsigma, value, limits, relative) 
    94     return value, weight / np.sum(weight) 
    9531 
    9632class DataMixin(object): 
    … …  
    11652    """ 
    11753    def _interpret_data(self, data, model): 
    118         # type: (Data, KernelModel) -> None 
    11954        # pylint: disable=attribute-defined-outside-init 
    12055 
    … …  
    13267            self.data_type = 'Iq' 
    13368 
     69        partype = model.info['partype'] 
     70 
    13471        if self.data_type == 'sesans': 
    13572             
    … …  
    14582            q_mono = sesans.make_all_q(data) 
    14683        elif self.data_type == 'Iqxy': 
    147             #if not model.info.parameters.has_2d: 
    148             #    raise ValueError("not 2D without orientation or magnetic parameters") 
     84            if not partype['orientation'] and not partype['magnetic']: 
     85                raise ValueError("not 2D without orientation or magnetic parameters") 
    14986            q = np.sqrt(data.qx_data**2 + data.qy_data**2) 
    15087            qmin = getattr(data, 'qmin', 1e-16) 
    … …  
    217154 
    218155    def _set_data(self, Iq, noise=None): 
    219         # type: (np.ndarray, Optional[float]) -> None 
    220156        # pylint: disable=attribute-defined-outside-init 
    221157        if noise is not None: 
    … …  
    235171 
    236172    def _calc_theory(self, pars, cutoff=0.0): 
    237         # type: (ParameterSet, float) -> np.ndarray 
    238173        if self._kernel is None: 
    239             self._kernel = self._model.make_kernel(self._kernel_inputs) 
     174            self._kernel = self._model.make_kernel(self._kernel_inputs)  # pylint: disable=attribute-dedata_type 
    240175            self._kernel_mono = (self._model.make_kernel(self._kernel_mono_inputs) 
    241176                                 if self._kernel_mono_inputs else None) 
    … …  
    271206    """ 
    272207    def __init__(self, data, model, cutoff=1e-5): 
    273         # type: (Data, KernelModel, float) -> None 
    274208        self.model = model 
    275209        self.cutoff = cutoff 
    … …  
    278212 
    279213    def __call__(self, **pars): 
    280         # type: (**float) -> np.ndarray 
    281214        return self._calc_theory(pars, cutoff=self.cutoff) 
    282215 
     216    def ER_VR(self, **pars): 
     217        """ 
     218        Compute the equivalent radius and volume ratio for the model. 
     219        """ 
     220        return call_ER_VR(self.model.info, pars) 
     221 
    283222    def simulate_data(self, noise=None, **pars): 
    284         # type: (Optional[float], **float) -> None 
    285223        """ 
    286224        Generate simulated data for the model. 
    … …  
    290228 
    291229def main(): 
    292     # type: () -> None 
    293230    """ 
    294231    Program to evaluate a particular model at a set of q values. 
    … …  
    306243        try: 
    307244            values = [float(v) for v in call.split(',')] 
    308         except Exception: 
     245        except: 
    309246            values = [] 
    310247        if len(values) == 1: 

  • sasmodels/generate.py

    ra5b8477 r81ec7c8  
    2121 
    2222    *VR(p1, p2, ...)* returns the volume ratio for core-shell style forms. 
    23  
    24     #define INVALID(v) (expr)  returns False if v.parameter is invalid 
    25     for some parameter or other (e.g., v.bell_radius < v.radius).  If 
    26     necessary, the expression can call a function. 
    2723 
    2824These functions are defined in a kernel module .py script and an associated 
    … …  
    6965The constructor code will generate the necessary vectors for computing 
    7066them with the desired polydispersity. 
     67 
     68The available kernel parameters are defined as a list, with each parameter 
     69defined as a sublist with the following elements: 
     70 
     71    *name* is the name that will be used in the call to the kernel 
     72    function and the name that will be displayed to the user.  Names 
     73    should be lower case, with words separated by underscore.  If 
     74    acronyms are used, the whole acronym should be upper case. 
     75 
     76    *units* should be one of *degrees* for angles, *Ang* for lengths, 
     77    *1e-6/Ang^2* for SLDs. 
     78 
     79    *default value* will be the initial value for  the model when it 
     80    is selected, or when an initial value is not otherwise specified. 
     81 
     82    *limits = [lb, ub]* are the hard limits on the parameter value, used to 
     83    limit the polydispersity density function.  In the fit, the parameter limits 
     84    given to the fit are the limits  on the central value of the parameter. 
     85    If there is polydispersity, it will evaluate parameter values outside 
     86    the fit limits, but not outside the hard limits specified in the model. 
     87    If there are no limits, use +/-inf imported from numpy. 
     88 
     89    *type* indicates how the parameter will be used.  "volume" parameters 
     90    will be used in all functions.  "orientation" parameters will be used 
     91    in *Iqxy* and *Imagnetic*.  "magnetic* parameters will be used in 
     92    *Imagnetic* only.  If *type* is the empty string, the parameter will 
     93    be used in all of *Iq*, *Iqxy* and *Imagnetic*.  "sld" parameters 
     94    can automatically be promoted to magnetic parameters, each of which 
     95    will have a magnitude and a direction, which may be different from 
     96    other sld parameters. 
     97 
     98    *description* is a short description of the parameter.  This will 
     99    be displayed in the parameter table and used as a tool tip for the 
     100    parameter value in the user interface. 
     101 
    71102The kernel module must set variables defining the kernel meta data: 
    72103 
    … …  
    123154    polydispersity values for tests. 
    124155 
    125 A :class:`modelinfo.ModelInfo` structure is constructed from the kernel meta 
    126 data and returned to the caller. 
     156An *model_info* dictionary is constructed from the kernel meta data and 
     157returned to the caller. 
     158 
     159The model evaluator, function call sequence consists of q inputs and the return vector, 
     160followed by the loop value/weight vector, followed by the values for 
     161the non-polydisperse parameters, followed by the lengths of the 
     162polydispersity loops.  To construct the call for 1D models, the 
     163categories *fixed-1d* and *pd-1d* list the names of the parameters 
     164of the non-polydisperse and the polydisperse parameters respectively. 
     165Similarly, *fixed-2d* and *pd-2d* provide parameter names for 2D models. 
     166The *pd-rel* category is a set of those parameters which give 
     167polydispersitiy as a portion of the value (so a 10% length dispersity 
     168would use a polydispersity value of 0.1) rather than absolute 
     169dispersity such as an angle plus or minus 15 degrees. 
     170 
     171The *volume* category lists the volume parameters in order for calls 
     172to volume within the kernel (used for volume normalization) and for 
     173calls to ER and VR for effective radius and volume ratio respectively. 
     174 
     175The *orientation* and *magnetic* categories list the orientation and 
     176magnetic parameters.  These are used by the sasview interface.  The 
     177blank category is for parameters such as scale which don't have any 
     178other marking. 
    127179 
    128180The doc string at the start of the kernel module will be used to 
    … …  
    132184file systems are case-sensitive, so only use lower case characters for 
    133185file names and extensions. 
     186 
     187 
     188The function :func:`make` loads the metadata from the module and returns 
     189the kernel source.  The function :func:`make_doc` extracts the doc string 
     190and adds the parameter table to the top.  The function :func:`model_sources` 
     191returns a list of files required by the model. 
    134192 
    135193Code follows the C99 standard with the following extensions and conditions:: 
    … …  
    143201    all double declarations may be converted to half, float, or long double 
    144202    FLOAT_SIZE is the number of bytes in the converted variables 
    145  
    146 :func:`load_kernel_module` loads the model definition file and 
    147 :modelinfo:`make_model_info` parses it. :func:`make_source` 
    148 converts C-based model definitions to C source code, including the 
    149 polydispersity integral.  :func:`model_sources` returns the list of 
    150 source files the model depends on, and :func:`timestamp` returns 
    151 the latest time stamp amongst the source files (so you can check if 
    152 the model needs to be rebuilt). 
    153  
    154 The function :func:`make_doc` extracts the doc string and adds the 
    155 parameter table to the top.  *make_figure* in *sasmodels/doc/genmodel* 
    156 creates the default figure for the model.  [These two sets of code 
    157 should mignrate into docs.py so docs can be updated in one place]. 
    158203""" 
    159204from __future__ import print_function 
    160205 
     206#TODO: identify model files which have changed since loading and reload them. 
    161207#TODO: determine which functions are useful outside of generate 
    162208#__all__ = ["model_info", "make_doc", "make_source", "convert_type"] 
    163209 
    164 from os.path import abspath, dirname, join as joinpath, exists, getmtime 
     210import sys 
     211from os.path import abspath, dirname, join as joinpath, exists, basename, \ 
     212    splitext 
    165213import re 
    166214import string 
    167215import warnings 
    168  
    169 import numpy as np  # type: ignore 
    170  
    171 from .modelinfo import Parameter 
     216from collections import namedtuple 
     217 
     218import numpy as np 
     219 
    172220from .custom import load_custom_kernel_module 
    173221 
    174 try: 
    175     from typing import Tuple, Sequence, Iterator 
    176     from .modelinfo import ModelInfo 
    177 except ImportError: 
    178     pass 
    179  
    180 TEMPLATE_ROOT = dirname(__file__) 
     222PARAMETER_FIELDS = ['name', 'units', 'default', 'limits', 'type', 'description'] 
     223Parameter = namedtuple('Parameter', PARAMETER_FIELDS) 
     224 
     225C_KERNEL_TEMPLATE_PATH = joinpath(dirname(__file__), 'kernel_template.c') 
    181226 
    182227F16 = np.dtype('float16') 
    … …  
    187232except TypeError: 
    188233    F128 = None 
     234 
     235# Scale and background, which are parameters common to every form factor 
     236COMMON_PARAMETERS = [ 
     237    ["scale", "", 1, [0, np.inf], "", "Source intensity"], 
     238    ["background", "1/cm", 1e-3, [0, np.inf], "", "Source background"], 
     239    ] 
    189240 
    190241# Conversion from units defined in the parameter table for each model 
    … …  
    233284 
    234285def format_units(units): 
    235     # type: (str) -> str 
    236286    """ 
    237287    Convert units into ReStructured Text format. 
    … …  
    240290 
    241291def make_partable(pars): 
    242     # type: (List[Parameter]) -> str 
    243292    """ 
    244293    Generate the parameter table to include in the sphinx documentation. 
    … …  
    271320 
    272321def _search(search_path, filename): 
    273     # type: (List[str], str) -> str 
    274322    """ 
    275323    Find *filename* in *search_path*. 
    … …  
    283331    raise ValueError("%r not found in %s" % (filename, search_path)) 
    284332 
    285  
    286333def model_sources(model_info): 
    287     # type: (ModelInfo) -> List[str] 
    288334    """ 
    289335    Return a list of the sources file paths for the module. 
    290336    """ 
    291     search_path = [dirname(model_info.filename), 
     337    search_path = [dirname(model_info['filename']), 
    292338                   abspath(joinpath(dirname(__file__), 'models'))] 
    293     return [_search(search_path, f) for f in model_info.source] 
    294  
    295 def timestamp(model_info): 
    296     # type: (ModelInfo) -> int 
    297     """ 
    298     Return a timestamp for the model corresponding to the most recently 
    299     changed file or dependency. 
    300     """ 
    301     source_files = (model_sources(model_info) 
    302                     + model_templates() 
    303                     + [model_info.filename]) 
    304     newest = max(getmtime(f) for f in source_files) 
    305     return newest 
     339    return [_search(search_path, f) for f in model_info['source']] 
     340 
     341# Pragmas for enable OpenCL features.  Be sure to protect them so that they 
     342# still compile even if OpenCL is not present. 
     343_F16_PRAGMA = """\ 
     344#if defined(__OPENCL_VERSION__) // && !defined(cl_khr_fp16) 
     345#  pragma OPENCL EXTENSION cl_khr_fp16: enable 
     346#endif 
     347""" 
     348 
     349_F64_PRAGMA = """\ 
     350#if defined(__OPENCL_VERSION__) // && !defined(cl_khr_fp64) 
     351#  pragma OPENCL EXTENSION cl_khr_fp64: enable 
     352#endif 
     353""" 
    306354 
    307355def convert_type(source, dtype): 
    308     # type: (str, np.dtype) -> str 
    309356    """ 
    310357    Convert code from double precision to the desired type. 
    … …  
    315362    if dtype == F16: 
    316363        fbytes = 2 
    317         source = _convert_type(source, "half", "f") 
     364        source = _F16_PRAGMA + _convert_type(source, "half", "f") 
    318365    elif dtype == F32: 
    319366        fbytes = 4 
    … …  
    321368    elif dtype == F64: 
    322369        fbytes = 8 
    323         # no need to convert if it is already double 
     370        source = _F64_PRAGMA + source  # Source is already double 
    324371    elif dtype == F128: 
    325372        fbytes = 16 
    … …  
    331378 
    332379def _convert_type(source, type_name, constant_flag): 
    333     # type: (str, str, str) -> str 
    334380    """ 
    335381    Replace 'double' with *type_name* in *source*, tagging floating point 
    … …  
    350396 
    351397def kernel_name(model_info, is_2d): 
    352     # type: (ModelInfo, bool) -> str 
    353398    """ 
    354399    Name of the exported kernel symbol. 
    355400    """ 
    356     return model_info.name + "_" + ("Iqxy" if is_2d else "Iq") 
     401    return model_info['name'] + "_" + ("Iqxy" if is_2d else "Iq") 
    357402 
    358403 
    359404def indent(s, depth): 
    360     # type: (str, int) -> str 
    361405    """ 
    362406    Indent a string of text with *depth* additional spaces on each line. 
    … …  
    367411 
    368412 
    369 _template_cache = {}  # type: Dict[str, Tuple[int, str, str]] 
    370 def load_template(filename): 
    371     # type: (str) -> str 
    372     path = joinpath(TEMPLATE_ROOT, filename) 
    373     mtime = getmtime(path) 
    374     if filename not in _template_cache or mtime > _template_cache[filename][0]: 
    375         with open(path) as fid: 
    376             _template_cache[filename] = (mtime, fid.read(), path) 
    377     return _template_cache[filename][1] 
    378  
    379 def model_templates(): 
    380     # type: () -> List[str] 
    381     # TODO: fails DRY; templates are listed in two places. 
    382     # should instead have model_info contain a list of paths 
    383     return [joinpath(TEMPLATE_ROOT, filename) 
    384             for filename in ('kernel_header.c', 'kernel_iq.c')] 
    385  
    386  
    387 _FN_TEMPLATE = """\ 
    388 double %(name)s(%(pars)s); 
    389 double %(name)s(%(pars)s) { 
    390     %(body)s 
    391 } 
    392  
    393  
     413LOOP_OPEN = """\ 
     414for (int %(name)s_i=0; %(name)s_i < N%(name)s; %(name)s_i++) { 
     415  const double %(name)s = loops[2*(%(name)s_i%(offset)s)]; 
     416  const double %(name)s_w = loops[2*(%(name)s_i%(offset)s)+1];\ 
    394417""" 
    395  
    396 def _gen_fn(name, pars, body): 
    397     # type: (str, List[Parameter], str) -> str 
    398     """ 
    399     Generate a function given pars and body. 
    400  
    401     Returns the following string:: 
    402  
    403          double fn(double a, double b, ...); 
    404          double fn(double a, double b, ...) { 
    405              .... 
    406          } 
    407     """ 
    408     par_decl = ', '.join(p.as_function_argument() for p in pars) if pars else 'void' 
    409     return _FN_TEMPLATE % {'name': name, 'body': body, 'pars': par_decl} 
    410  
    411 def _call_pars(prefix, pars): 
    412     # type: (str, List[Parameter]) -> List[str] 
    413     """ 
    414     Return a list of *prefix.parameter* from parameter items. 
    415     """ 
    416     return [p.as_call_reference(prefix) for p in pars] 
    417  
    418 _IQXY_PATTERN = re.compile("^((inline|static) )? *(double )? *Iqxy *([(]|$)", 
    419                            flags=re.MULTILINE) 
    420 def _have_Iqxy(sources): 
    421     # type: (List[str]) -> bool 
    422     """ 
    423     Return true if any file defines Iqxy. 
    424  
    425     Note this is not a C parser, and so can be easily confused by 
    426     non-standard syntax.  Also, it will incorrectly identify the following 
    427     as having Iqxy:: 
    428  
    429         /* 
    430         double Iqxy(qx, qy, ...) { ... fill this in later ... } 
    431         */ 
    432  
    433     If you want to comment out an Iqxy function, use // on the front of the 
    434     line instead. 
    435     """ 
    436     for code in sources: 
    437         if _IQXY_PATTERN.search(code): 
    438             return True 
    439     else: 
    440         return False 
    441  
     418def build_polydispersity_loops(pd_pars): 
     419    """ 
     420    Build polydispersity loops 
     421 
     422    Returns loop opening and loop closing 
     423    """ 
     424    depth = 4 
     425    offset = "" 
     426    loop_head = [] 
     427    loop_end = [] 
     428    for name in pd_pars: 
     429        subst = {'name': name, 'offset': offset} 
     430        loop_head.append(indent(LOOP_OPEN % subst, depth)) 
     431        loop_end.insert(0, (" "*depth) + "}") 
     432        offset += '+N' + name 
     433        depth += 2 
     434    return "\n".join(loop_head), "\n".join(loop_end) 
     435 
     436C_KERNEL_TEMPLATE = None 
    442437def make_source(model_info): 
    443     # type: (ModelInfo) -> str 
    444438    """ 
    445439    Generate the OpenCL/ctypes kernel from the module info. 
    446440 
    447     Uses source files found in the given search path.  Returns None if this 
    448     is a pure python model, with no C source components. 
    449     """ 
    450     if callable(model_info.Iq): 
    451         raise ValueError("can't compile python model") 
     441    Uses source files found in the given search path. 
     442    """ 
     443    if callable(model_info['Iq']): 
     444        return None 
    452445 
    453446    # TODO: need something other than volume to indicate dispersion parameters 
    … …  
    460453    # for computing volume even if we allow non-disperse volume parameters. 
    461454 
    462     partable = model_info.parameters 
    463  
    464     # Identify parameters for Iq, Iqxy, Iq_magnetic and form_volume. 
    465     # Note that scale and volume are not possible types. 
    466  
    467     # Load templates and user code 
    468     kernel_header = load_template('kernel_header.c') 
    469     kernel_code = load_template('kernel_iq.c') 
    470     user_code = [open(f).read() for f in model_sources(model_info)] 
    471  
    472     # Build initial sources 
    473     source = [kernel_header] + user_code 
    474  
    475     # Make parameters for q, qx, qy so that we can use them in declarations 
    476     q, qx, qy = [Parameter(name=v) for v in ('q', 'qx', 'qy')] 
    477     # Generate form_volume function, etc. from body only 
    478     if isinstance(model_info.form_volume, str): 
    479         pars = partable.form_volume_parameters 
    480         source.append(_gen_fn('form_volume', pars, model_info.form_volume)) 
    481     if isinstance(model_info.Iq, str): 
    482         pars = [q] + partable.iq_parameters 
    483         source.append(_gen_fn('Iq', pars, model_info.Iq)) 
    484     if isinstance(model_info.Iqxy, str): 
    485         pars = [qx, qy] + partable.iqxy_parameters 
    486         source.append(_gen_fn('Iqxy', pars, model_info.Iqxy)) 
    487  
    488     # Define the parameter table 
    489     source.append("#define PARAMETER_TABLE \\") 
    490     source.append("\\\n".join(p.as_definition() 
    491                               for p in partable.kernel_parameters)) 
    492  
    493     # Define the function calls 
    494     if partable.form_volume_parameters: 
    495         refs = _call_pars("_v.", partable.form_volume_parameters) 
    496         call_volume = "#define CALL_VOLUME(_v) form_volume(%s)" % (",".join(refs)) 
    497     else: 
    498         # Model doesn't have volume.  We could make the kernel run a little 
    499         # faster by not using/transferring the volume normalizations, but 
    500         # the ifdef's reduce readability more than is worthwhile. 
    501         call_volume = "#define CALL_VOLUME(v) 1.0" 
    502     source.append(call_volume) 
    503  
    504     refs = ["_q[_i]"] + _call_pars("_v.", partable.iq_parameters) 
    505     call_iq = "#define CALL_IQ(_q,_i,_v) Iq(%s)" % (",".join(refs)) 
    506     if _have_Iqxy(user_code): 
    507         # Call 2D model 
    508         refs = ["q[2*i]", "q[2*i+1]"] + _call_pars("_v.", partable.iqxy_parameters) 
    509         call_iqxy = "#define CALL_IQ(_q,_i,_v) Iqxy(%s)" % (",".join(refs)) 
    510     else: 
    511         # Call 1D model with sqrt(qx^2 + qy^2) 
    512         warnings.warn("Creating Iqxy = Iq(sqrt(qx^2 + qy^2))") 
    513         # still defined:: refs = ["q[i]"] + _call_pars("v", iq_parameters) 
    514         pars_sqrt = ["sqrt(_q[2*_i]*_q[2*_i]+_q[2*_i+1]*_q[2*_i+1])"] + refs[1:] 
    515         call_iqxy = "#define CALL_IQ(_q,_i,_v) Iq(%s)" % (",".join(pars_sqrt)) 
    516  
    517     # Fill in definitions for numbers of parameters 
    518     source.append("#define MAX_PD %s"%partable.max_pd) 
    519     source.append("#define NPARS %d"%partable.npars) 
    520  
    521     # TODO: allow mixed python/opencl kernels? 
    522  
    523     # define the Iq kernel 
    524     source.append("#define KERNEL_NAME %s_Iq"%model_info.name) 
    525     source.append(call_iq) 
    526     source.append(kernel_code) 
    527     source.append("#undef CALL_IQ") 
    528     source.append("#undef KERNEL_NAME") 
    529  
    530     # define the Iqxy kernel from the same source with different #defines 
    531     source.append("#define KERNEL_NAME %s_Iqxy"%model_info.name) 
    532     source.append(call_iqxy) 
    533     source.append(kernel_code) 
    534     source.append("#undef CALL_IQ") 
    535     source.append("#undef KERNEL_NAME") 
    536  
    537     return '\n'.join(source) 
     455    # Load template 
     456    global C_KERNEL_TEMPLATE 
     457    if C_KERNEL_TEMPLATE is None: 
     458        with open(C_KERNEL_TEMPLATE_PATH) as fid: 
     459            C_KERNEL_TEMPLATE = fid.read() 
     460 
     461    # Load additional sources 
     462    source = [open(f).read() for f in model_sources(model_info)] 
     463 
     464    # Prepare defines 
     465    defines = [] 
     466    partype = model_info['partype'] 
     467    pd_1d = partype['pd-1d'] 
     468    pd_2d = partype['pd-2d'] 
     469    fixed_1d = partype['fixed-1d'] 
     470    fixed_2d = partype['fixed-1d'] 
     471 
     472    iq_parameters = [p.name 
     473                     for p in model_info['parameters'][2:]  # skip scale, background 
     474                     if p.name in set(fixed_1d + pd_1d)] 
     475    iqxy_parameters = [p.name 
     476                       for p in model_info['parameters'][2:]  # skip scale, background 
     477                       if p.name in set(fixed_2d + pd_2d)] 
     478    volume_parameters = [p.name 
     479                         for p in model_info['parameters'] 
     480                         if p.type == 'volume'] 
     481 
     482    # Fill in defintions for volume parameters 
     483    if volume_parameters: 
     484        defines.append(('VOLUME_PARAMETERS', 
     485                        ','.join(volume_parameters))) 
     486        defines.append(('VOLUME_WEIGHT_PRODUCT', 
     487                        '*'.join(p + '_w' for p in volume_parameters))) 
     488 
     489    # Generate form_volume function from body only 
     490    if model_info['form_volume'] is not None: 
     491        if volume_parameters: 
     492            vol_par_decl = ', '.join('double ' + p for p in volume_parameters) 
     493        else: 
     494            vol_par_decl = 'void' 
     495        defines.append(('VOLUME_PARAMETER_DECLARATIONS', 
     496                        vol_par_decl)) 
     497        fn = """\ 
     498double form_volume(VOLUME_PARAMETER_DECLARATIONS); 
     499double form_volume(VOLUME_PARAMETER_DECLARATIONS) { 
     500    %(body)s 
     501} 
     502""" % {'body':model_info['form_volume']} 
     503        source.append(fn) 
     504 
     505    # Fill in definitions for Iq parameters 
     506    defines.append(('IQ_KERNEL_NAME', model_info['name'] + '_Iq')) 
     507    defines.append(('IQ_PARAMETERS', ', '.join(iq_parameters))) 
     508    if fixed_1d: 
     509        defines.append(('IQ_FIXED_PARAMETER_DECLARATIONS', 
     510                        ', \\\n    '.join('const double %s' % p for p in fixed_1d))) 
     511    if pd_1d: 
     512        defines.append(('IQ_WEIGHT_PRODUCT', 
     513                        '*'.join(p + '_w' for p in pd_1d))) 
     514        defines.append(('IQ_DISPERSION_LENGTH_DECLARATIONS', 
     515                        ', \\\n    '.join('const int N%s' % p for p in pd_1d))) 
     516        defines.append(('IQ_DISPERSION_LENGTH_SUM', 
     517                        '+'.join('N' + p for p in pd_1d))) 
     518        open_loops, close_loops = build_polydispersity_loops(pd_1d) 
     519        defines.append(('IQ_OPEN_LOOPS', 
     520                        open_loops.replace('\n', ' \\\n'))) 
     521        defines.append(('IQ_CLOSE_LOOPS', 
     522                        close_loops.replace('\n', ' \\\n'))) 
     523    if model_info['Iq'] is not None: 
     524        defines.append(('IQ_PARAMETER_DECLARATIONS', 
     525                        ', '.join('double ' + p for p in iq_parameters))) 
     526        fn = """\ 
     527double Iq(double q, IQ_PARAMETER_DECLARATIONS); 
     528double Iq(double q, IQ_PARAMETER_DECLARATIONS) { 
     529    %(body)s 
     530} 
     531""" % {'body':model_info['Iq']} 
     532        source.append(fn) 
     533 
     534    # Fill in definitions for Iqxy parameters 
     535    defines.append(('IQXY_KERNEL_NAME', model_info['name'] + '_Iqxy')) 
     536    defines.append(('IQXY_PARAMETERS', ', '.join(iqxy_parameters))) 
     537    if fixed_2d: 
     538        defines.append(('IQXY_FIXED_PARAMETER_DECLARATIONS', 
     539                        ', \\\n    '.join('const double %s' % p for p in fixed_2d))) 
     540    if pd_2d: 
     541        defines.append(('IQXY_WEIGHT_PRODUCT', 
     542                        '*'.join(p + '_w' for p in pd_2d))) 
     543        defines.append(('IQXY_DISPERSION_LENGTH_DECLARATIONS', 
     544                        ', \\\n    '.join('const int N%s' % p for p in pd_2d))) 
     545        defines.append(('IQXY_DISPERSION_LENGTH_SUM', 
     546                        '+'.join('N' + p for p in pd_2d))) 
     547        open_loops, close_loops = build_polydispersity_loops(pd_2d) 
     548        defines.append(('IQXY_OPEN_LOOPS', 
     549                        open_loops.replace('\n', ' \\\n'))) 
     550        defines.append(('IQXY_CLOSE_LOOPS', 
     551                        close_loops.replace('\n', ' \\\n'))) 
     552    if model_info['Iqxy'] is not None: 
     553        defines.append(('IQXY_PARAMETER_DECLARATIONS', 
     554                        ', '.join('double ' + p for p in iqxy_parameters))) 
     555        fn = """\ 
     556double Iqxy(double qx, double qy, IQXY_PARAMETER_DECLARATIONS); 
     557double Iqxy(double qx, double qy, IQXY_PARAMETER_DECLARATIONS) { 
     558    %(body)s 
     559} 
     560""" % {'body':model_info['Iqxy']} 
     561        source.append(fn) 
     562 
     563    # Need to know if we have a theta parameter for Iqxy; it is not there 
     564    # for the magnetic sphere model, for example, which has a magnetic 
     565    # orientation but no shape orientation. 
     566    if 'theta' in pd_2d: 
     567        defines.append(('IQXY_HAS_THETA', '1')) 
     568 
     569    #for d in defines: print(d) 
     570    defines = '\n'.join('#define %s %s' % (k, v) for k, v in defines) 
     571    sources = '\n\n'.join(source) 
     572    return C_KERNEL_TEMPLATE % { 
     573        'DEFINES': defines, 
     574        'SOURCES': sources, 
     575        } 
     576 
     577def categorize_parameters(pars): 
     578    """ 
     579    Build parameter categories out of the the parameter definitions. 
     580 
     581    Returns a dictionary of categories. 
     582 
     583    Note: these categories are subject to change, depending on the needs of 
     584    the UI and the needs of the kernel calling function. 
     585 
     586    The categories are as follows: 
     587 
     588    * *volume* list of volume parameter names 
     589    * *orientation* list of orientation parameters 
     590    * *magnetic* list of magnetic parameters 
     591    * *<empty string>* list of parameters that have no type info 
     592 
     593    Each parameter is in one and only one category. 
     594 
     595    The following derived categories are created: 
     596 
     597    * *fixed-1d* list of non-polydisperse parameters for 1D models 
     598    * *pd-1d* list of polydisperse parameters for 1D models 
     599    * *fixed-2d* list of non-polydisperse parameters for 2D models 
     600    * *pd-d2* list of polydisperse parameters for 2D models 
     601    """ 
     602    partype = { 
     603        'volume': [], 'orientation': [], 'magnetic': [], 'sld': [], '': [], 
     604        'fixed-1d': [], 'fixed-2d': [], 'pd-1d': [], 'pd-2d': [], 
     605        'pd-rel': set(), 
     606    } 
     607 
     608    for p in pars: 
     609        if p.type == 'volume': 
     610            partype['pd-1d'].append(p.name) 
     611            partype['pd-2d'].append(p.name) 
     612            partype['pd-rel'].add(p.name) 
     613        elif p.type == 'magnetic': 
     614            partype['fixed-2d'].append(p.name) 
     615        elif p.type == 'orientation': 
     616            partype['pd-2d'].append(p.name) 
     617        elif p.type in ('', 'sld'): 
     618            partype['fixed-1d'].append(p.name) 
     619            partype['fixed-2d'].append(p.name) 
     620        else: 
     621            raise ValueError("unknown parameter type %r" % p.type) 
     622        partype[p.type].append(p.name) 
     623 
     624    return partype 
     625 
     626def process_parameters(model_info): 
     627    """ 
     628    Process parameter block, precalculating parameter details. 
     629    """ 
     630    # convert parameters into named tuples 
     631    for p in model_info['parameters']: 
     632        if p[4] == '' and (p[0].startswith('sld') or p[0].endswith('sld')): 
     633            p[4] = 'sld' 
     634            # TODO: make sure all models explicitly label their sld parameters 
     635            #raise ValueError("%s.%s needs to be explicitly set to type 'sld'" %(model_info['id'], p[0])) 
     636 
     637    pars = [Parameter(*p) for p in model_info['parameters']] 
     638    # Fill in the derived attributes 
     639    model_info['parameters'] = pars 
     640    partype = categorize_parameters(pars) 
     641    model_info['limits'] = dict((p.name, p.limits) for p in pars) 
     642    model_info['partype'] = partype 
     643    model_info['defaults'] = dict((p.name, p.default) for p in pars) 
     644    if model_info.get('demo', None) is None: 
     645        model_info['demo'] = model_info['defaults'] 
     646    model_info['has_2d'] = partype['orientation'] or partype['magnetic'] 
     647 
    538648 
    539649def load_kernel_module(model_name): 
    540     # type: (str) -> module 
    541650    if model_name.endswith('.py'): 
    542651        kernel_module = load_custom_kernel_module(model_name) 
    … …  
    548657 
    549658 
     659def make_model_info(kernel_module): 
     660    """ 
     661    Interpret the model definition file, categorizing the parameters. 
     662 
     663    The module can be loaded with a normal python import statement if you 
     664    know which module you need, or with __import__('sasmodels.model.'+name) 
     665    if the name is in a string. 
     666 
     667    The *model_info* structure contains the following fields: 
     668 
     669    * *id* is the id of the kernel 
     670    * *name* is the display name of the kernel 
     671    * *filename* is the full path to the module defining the file (if any) 
     672    * *title* is a short description of the kernel 
     673    * *description* is a long description of the kernel (this doesn't seem 
     674      very useful since the Help button on the model page brings you directly 
     675      to the documentation page) 
     676    * *docs* is the docstring from the module.  Use :func:`make_doc` to 
     677    * *category* specifies the model location in the docs 
     678    * *parameters* is the model parameter table 
     679    * *single* is True if the model allows single precision 
     680    * *structure_factor* is True if the model is useable in a product 
     681    * *defaults* is the *{parameter: value}* table built from the parameter 
     682      description table. 
     683    * *limits* is the *{parameter: [min, max]}* table built from the 
     684      parameter description table. 
     685    * *partypes* categorizes the model parameters. See 
     686      :func:`categorize_parameters` for details. 
     687    * *demo* contains the *{parameter: value}* map used in compare (and maybe 
     688      for the demo plot, if plots aren't set up to use the default values). 
     689      If *demo* is not given in the file, then the default values will be used. 
     690    * *tests* is a set of tests that must pass 
     691    * *source* is the list of library files to include in the C model build 
     692    * *Iq*, *Iqxy*, *form_volume*, *ER*, *VR* and *sesans* are python functions 
     693      implementing the kernel for the module, or None if they are not 
     694      defined in python 
     695    * *composition* is None if the model is independent, otherwise it is a 
     696      tuple with composition type ('product' or 'mixture') and a list of 
     697      *model_info* blocks for the composition objects.  This allows us to 
     698      build complete product and mixture models from just the info. 
     699 
     700    """ 
     701    parameters = COMMON_PARAMETERS + kernel_module.parameters 
     702    filename = abspath(kernel_module.__file__) 
     703    kernel_id = splitext(basename(filename))[0] 
     704    name = getattr(kernel_module, 'name', None) 
     705    if name is None: 
     706        name = " ".join(w.capitalize() for w in kernel_id.split('_')) 
     707    model_info = dict( 
     708        id=kernel_id,  # string used to load the kernel 
     709        filename=abspath(kernel_module.__file__), 
     710        name=name, 
     711        title=kernel_module.title, 
     712        description=kernel_module.description, 
     713        parameters=parameters, 
     714        composition=None, 
     715        docs=kernel_module.__doc__, 
     716        category=getattr(kernel_module, 'category', None), 
     717        single=getattr(kernel_module, 'single', True), 
     718        structure_factor=getattr(kernel_module, 'structure_factor', False), 
     719        control=getattr(kernel_module, 'control', None), 
     720        demo=getattr(kernel_module, 'demo', None), 
     721        source=getattr(kernel_module, 'source', []), 
     722        tests=getattr(kernel_module, 'tests', []), 
     723        ) 
     724    process_parameters(model_info) 
     725    # Check for optional functions 
     726    functions = "ER VR form_volume Iq Iqxy shape sesans".split() 
     727    model_info.update((k, getattr(kernel_module, k, None)) for k in functions) 
     728    return model_info 
    550729 
    551730section_marker = re.compile(r'\A(?P<first>[%s])(?P=first)*\Z' 
    552731                            %re.escape(string.punctuation)) 
    553732def _convert_section_titles_to_boldface(lines): 
    554     # type: (Sequence[str]) -> Iterator[str] 
    555733    """ 
    556734    Do the actual work of identifying and converting section headings. 
    … …  
    574752 
    575753def convert_section_titles_to_boldface(s): 
    576     # type: (str) -> str 
    577754    """ 
    578755    Use explicit bold-face rather than section headings so that the table of 
    … …  
    585762 
    586763def make_doc(model_info): 
    587     # type: (ModelInfo) -> str 
    588764    """ 
    589765    Return the documentation for the model. 
    … …  
    591767    Iq_units = "The returned value is scaled to units of |cm^-1| |sr^-1|, absolute scale." 
    592768    Sq_units = "The returned value is a dimensionless structure factor, $S(q)$." 
    593     docs = convert_section_titles_to_boldface(model_info.docs) 
    594     pars = make_partable(model_info.parameters.COMMON 
    595                          + model_info.parameters.kernel_parameters) 
    596     subst = dict(id=model_info.id.replace('_', '-'), 
    597                  name=model_info.name, 
    598                  title=model_info.title, 
    599                  parameters=pars, 
    600                  returns=Sq_units if model_info.structure_factor else Iq_units, 
     769    docs = convert_section_titles_to_boldface(model_info['docs']) 
     770    subst = dict(id=model_info['id'].replace('_', '-'), 
     771                 name=model_info['name'], 
     772                 title=model_info['title'], 
     773                 parameters=make_partable(model_info['parameters']), 
     774                 returns=Sq_units if model_info['structure_factor'] else Iq_units, 
    601775                 docs=docs) 
    602776    return DOC_HEADER % subst 
    … …  
    604778 
    605779def demo_time(): 
    606     # type: () -> None 
    607780    """ 
    608781    Show how long it takes to process a model. 
    609782    """ 
     783    from .models import cylinder 
    610784    import datetime 
    611     from .modelinfo import make_model_info 
    612     from .models import cylinder 
    613  
    614785    tic = datetime.datetime.now() 
    615786    make_source(make_model_info(cylinder)) 
    … …  
    618789 
    619790def main(): 
    620     # type: () -> None 
    621791    """ 
    622792    Program which prints the source produced by the model. 
    623793    """ 
    624     import sys 
    625     from .modelinfo import make_model_info 
    626  
    627794    if len(sys.argv) <= 1: 
    628795        print("usage: python -m sasmodels.generate modelname") 

  • sasmodels/kernelcl.py

    ra5b8477 r4d76711  
    4848harmless, albeit annoying. 
    4949""" 
    50 from __future__ import print_function 
    51  
    5250import os 
    5351import warnings 
    5452 
    55 import numpy as np  # type: ignore 
     53import numpy as np 
    5654 
    5755try: 
    58     raise NotImplementedError("OpenCL not yet implemented for new kernel template") 
    59     import pyopencl as cl  # type: ignore 
     56    import pyopencl as cl 
    6057    # Ask OpenCL for the default context so that we know that one exists 
    6158    cl.create_some_context(interactive=False) 
    … …  
    6461    raise RuntimeError("OpenCL not available") 
    6562 
    66 from pyopencl import mem_flags as mf  # type: ignore 
    67 from pyopencl.characterize import get_fast_inaccurate_build_options  # type: ignore 
     63from pyopencl import mem_flags as mf 
     64from pyopencl.characterize import get_fast_inaccurate_build_options 
    6865 
    6966from . import generate 
    70 from .kernel import KernelModel, Kernel 
    71  
    72 try: 
    73     from typing import Tuple, Callable, Any 
    74     from .modelinfo import ModelInfo 
    75     from .details import CallDetails 
    76 except ImportError: 
    77     pass 
    7867 
    7968# The max loops number is limited by the amount of local memory available 
    … …  
    8473# of polydisperse parameters. 
    8574MAX_LOOPS = 2048 
    86  
    87  
    88 # Pragmas for enable OpenCL features.  Be sure to protect them so that they 
    89 # still compile even if OpenCL is not present. 
    90 _F16_PRAGMA = """\ 
    91 #if defined(__OPENCL_VERSION__) // && !defined(cl_khr_fp16) 
    92 #  pragma OPENCL EXTENSION cl_khr_fp16: enable 
    93 #endif 
    94 """ 
    95  
    96 _F64_PRAGMA = """\ 
    97 #if defined(__OPENCL_VERSION__) // && !defined(cl_khr_fp64) 
    98 #  pragma OPENCL EXTENSION cl_khr_fp64: enable 
    99 #endif 
    100 """ 
    10175 
    10276 
    … …  
    168142        raise RuntimeError("%s not supported for devices"%dtype) 
    169143 
    170     source_list = [generate.convert_type(source, dtype)] 
    171  
    172     if dtype == generate.F16: 
    173         source_list.insert(0, _F16_PRAGMA) 
    174     elif dtype == generate.F64: 
    175         source_list.insert(0, _F64_PRAGMA) 
    176  
     144    source = generate.convert_type(source, dtype) 
    177145    # Note: USE_SINCOS makes the intel cpu slower under opencl 
    178146    if context.devices[0].type == cl.device_type.GPU: 
    179         source_list.insert(0, "#define USE_SINCOS\n") 
     147        source = "#define USE_SINCOS\n" + source 
    180148    options = (get_fast_inaccurate_build_options(context.devices[0]) 
    181149               if fast else []) 
    182     source = "\n".join(source_list) 
    183150    program = cl.Program(context, source).build(options=options) 
    184151    return program 
    … …  
    253220        key = "%s-%s-%s"%(name, dtype, fast) 
    254221        if key not in self.compiled: 
    255             print("compiling",name) 
     222            #print("compiling",name) 
    256223            dtype = np.dtype(dtype) 
    257             program = compile_model(self.get_context(dtype), 
    258                                     str(source), dtype, fast) 
     224            program = compile_model(self.get_context(dtype), source, dtype, fast) 
    259225            self.compiled[key] = program 
    260226        return self.compiled[key] 
    … …  
    319285 
    320286 
    321 class GpuModel(KernelModel): 
     287class GpuModel(object): 
    322288    """ 
    323289    GPU wrapper for a single model. 
    … …  
    351317        if self.program is None: 
    352318            compiler = environment().compile_program 
    353             self.program = compiler(self.info.name, self.source, 
    354                                     self.dtype, self.fast) 
     319            self.program = compiler(self.info['name'], self.source, self.dtype, 
     320                                    self.fast) 
    355321        is_2d = len(q_vectors) == 2 
    356322        kernel_name = generate.kernel_name(self.info, is_2d) 
    … …  
    363329        """ 
    364330        if self.program is not None: 
    365             environment().release_program(self.info.name) 
     331            environment().release_program(self.info['name']) 
    366332            self.program = None 
    367333 
    … …  
    399365        # at this point, so instead using 32, which is good on the set of 
    400366        # architectures tested so far. 
    401         if self.is_2d: 
    402             # Note: 17 rather than 15 because results is 2 elements 
    403             # longer than input. 
    404             width = ((self.nq+17)//16)*16 
    405             self.q = np.empty((width, 2), dtype=dtype) 
    406             self.q[:self.nq, 0] = q_vectors[0] 
    407             self.q[:self.nq, 1] = q_vectors[1] 
    408         else: 
    409             # Note: 33 rather than 31 because results is 2 elements 
    410             # longer than input. 
    411             width = ((self.nq+33)//32)*32 
    412             self.q = np.empty(width, dtype=dtype) 
    413             self.q[:self.nq] = q_vectors[0] 
    414         self.global_size = [self.q.shape[0]] 
     367        self.q_vectors = [_stretch_input(q, self.dtype, 32) for q in q_vectors] 
    415368        context = env.get_context(self.dtype) 
     369        self.global_size = [self.q_vectors[0].size] 
    416370        #print("creating inputs of size", self.global_size) 
    417         self.q_b = cl.Buffer(context, mf.READ_ONLY | mf.COPY_HOST_PTR, 
    418                              hostbuf=self.q) 
     371        self.q_buffers = [ 
     372            cl.Buffer(context, mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=q) 
     373            for q in self.q_vectors 
     374        ] 
    419375 
    420376    def release(self): 
    … …  
    422378        Free the memory. 
    423379        """ 
    424         if self.q is not None: 
    425             self.q.release() 
    426             self.q = None 
     380        for b in self.q_buffers: 
     381            b.release() 
     382        self.q_buffers = [] 
    427383 
    428384    def __del__(self): 
    429385        self.release() 
    430386 
    431 class GpuKernel(Kernel): 
     387class GpuKernel(object): 
    432388    """ 
    433389    Callable SAS kernel. 
    … …  
    449405    Call :meth:`release` when done with the kernel instance. 
    450406    """ 
    451     def __init__(self, kernel, model_info, q_vectors): 
    452         # type: (KernelModel, ModelInfo, List[np.ndarray]) -> None 
    453         max_pd = model_info.parameters.max_pd 
    454         npars = len(model_info.parameters.kernel_parameters)-2 
    455         q_input = GpuInput(q_vectors, kernel.dtype) 
     407    def __init__(self, kernel, model_info, q_vectors, dtype): 
     408        q_input = GpuInput(q_vectors, dtype) 
    456409        self.kernel = kernel 
    457410        self.info = model_info 
    458         self.dtype = kernel.dtype 
    459         self.dim = '2d' if q_input.is_2d else '1d' 
    460         self.pd_stop_index = 4*max_pd-1 
    461         # plus three for the normalization values 
    462         self.result = np.empty(q_input.nq+3, q_input.dtype) 
     411        self.res = np.empty(q_input.nq, q_input.dtype) 
     412        dim = '2d' if q_input.is_2d else '1d' 
     413        self.fixed_pars = model_info['partype']['fixed-' + dim] 
     414        self.pd_pars = model_info['partype']['pd-' + dim] 
    463415 
    464416        # Inputs and outputs for each kernel call 
    465417        # Note: res may be shorter than res_b if global_size != nq 
    466418        env = environment() 
    467         self.queue = env.get_queue(kernel.dtype) 
    468  
    469         # details is int32 data, padded to an 8 integer boundary 
    470         size = ((max_pd*5 + npars*3 + 2 + 7)//8)*8 
    471         self.result_b = cl.Buffer(self.queue.context, mf.READ_WRITE, 
    472                                q_input.global_size[0] * kernel.dtype.itemsize) 
    473         self.q_input = q_input # allocated by GpuInput above 
    474  
    475         self._need_release = [ self.result_b, self.q_input ] 
    476  
    477     def __call__(self, call_details, weights, values, cutoff): 
    478         # type: (CallDetails, np.ndarray, np.ndarray, float) -> np.ndarray 
     419        self.queue = env.get_queue(dtype) 
     420        self.loops_b = cl.Buffer(self.queue.context, mf.READ_WRITE, 
     421                                 2 * MAX_LOOPS * q_input.dtype.itemsize) 
     422        self.res_b = cl.Buffer(self.queue.context, mf.READ_WRITE, 
     423                               q_input.global_size[0] * q_input.dtype.itemsize) 
     424        self.q_input = q_input 
     425 
     426        self._need_release = [self.loops_b, self.res_b, self.q_input] 
     427 
     428    def __call__(self, fixed_pars, pd_pars, cutoff): 
    479429        real = (np.float32 if self.q_input.dtype == generate.F32 
    480430                else np.float64 if self.q_input.dtype == generate.F64 
    481431                else np.float16 if self.q_input.dtype == generate.F16 
    482432                else np.float32)  # will never get here, so use np.float32 
    483         assert call_details.dtype == np.int32 
    484         assert weights.dtype == real and values.dtype == real 
    485  
    486         context = self.queue.context 
    487         details_b = cl.Buffer(context, mf.READ_ONLY | mf.COPY_HOST_PTR, 
    488                               hostbuf=call_details) 
    489         weights_b = cl.Buffer(context, mf.READ_ONLY | mf.COPY_HOST_PTR, 
    490                               hostbuf=weights) 
    491         values_b = cl.Buffer(context, mf.READ_ONLY | mf.COPY_HOST_PTR, 
    492                              hostbuf=values) 
    493  
    494         start, stop = 0, self.details[self.pd_stop_index] 
    495         args = [ 
    496             np.uint32(self.q_input.nq), np.uint32(start), np.uint32(stop), 
    497             self.details_b, self.weights_b, self.values_b, 
    498             self.q_input.q_b, self.result_b, real(cutoff), 
    499         ] 
     433 
     434        #print "pars", fixed_pars, pd_pars 
     435        res_bi = self.res_b 
     436        nq = np.uint32(self.q_input.nq) 
     437        if pd_pars: 
     438            cutoff = real(cutoff) 
     439            loops_N = [np.uint32(len(p[0])) for p in pd_pars] 
     440            loops = np.hstack(pd_pars) \ 
     441                if pd_pars else np.empty(0, dtype=self.q_input.dtype) 
     442            loops = np.ascontiguousarray(loops.T, self.q_input.dtype).flatten() 
     443            #print("loops",Nloops, loops) 
     444 
     445            #import sys; print("opencl eval",pars) 
     446            #print("opencl eval",pars) 
     447            if len(loops) > 2 * MAX_LOOPS: 
     448                raise ValueError("too many polydispersity points") 
     449 
     450            loops_bi = self.loops_b 
     451            cl.enqueue_copy(self.queue, loops_bi, loops) 
     452            loops_l = cl.LocalMemory(len(loops.data)) 
     453            #ctx = environment().context 
     454            #loops_bi = cl.Buffer(ctx, mf.READ_ONLY|mf.COPY_HOST_PTR, hostbuf=loops) 
     455            dispersed = [loops_bi, loops_l, cutoff] + loops_N 
     456        else: 
     457            dispersed = [] 
     458        fixed = [real(p) for p in fixed_pars] 
     459        args = self.q_input.q_buffers + [res_bi, nq] + dispersed + fixed 
    500460        self.kernel(self.queue, self.q_input.global_size, None, *args) 
    501         cl.enqueue_copy(self.queue, self.result, self.result_b) 
    502         [v.release() for v in (details_b, weights_b, values_b)] 
    503  
    504         return self.result[:self.nq] 
     461        cl.enqueue_copy(self.queue, self.res, res_bi) 
     462 
     463        return self.res 
    505464 
    506465    def release(self): 

  • sasmodels/kerneldll.py

    ra5b8477 r4d76711  
    4949import os 
    5050import tempfile 
    51 import ctypes as ct  # type: ignore 
    52 from ctypes import c_void_p, c_int32, c_longdouble, c_double, c_float  # type: ignore 
    53  
    54 import numpy as np  # type: ignore 
     51import ctypes as ct 
     52from ctypes import c_void_p, c_int, c_longdouble, c_double, c_float 
     53import _ctypes 
     54 
     55import numpy as np 
    5556 
    5657from . import generate 
    57 from .kernel import KernelModel, Kernel 
    58 from .kernelpy import PyInput 
     58from .kernelpy import PyInput, PyModel 
    5959from .exception import annotate_exception 
    60 from .generate import F16, F32, F64 
    61  
    62 try: 
    63     from typing import Tuple, Callable, Any 
    64     from .modelinfo import ModelInfo 
    65     from .details import CallDetails 
    66 except ImportError: 
    67     pass 
    6860 
    6961# Compiler platform details 
    … …  
    8981        COMPILE = "gcc -shared -fPIC -std=c99 -O2 -Wall %(source)s -o %(output)s -lm" 
    9082        if "SAS_OPENMP" in os.environ: 
    91             COMPILE += " -fopenmp" 
     83            COMPILE = COMPILE + " -fopenmp" 
    9284else: 
    9385    COMPILE = "cc -shared -fPIC -fopenmp -std=c99 -O2 -Wall %(source)s -o %(output)s -lm" 
    … …  
    9890 
    9991 
    100 def dll_name(model_info, dtype): 
    101     # type: (ModelInfo, np.dtype) ->  str 
    102     """ 
    103     Name of the dll containing the model.  This is the base file name without 
    104     any path or extension, with a form such as 'sas_sphere32'. 
    105     """ 
    106     bits = 8*dtype.itemsize 
    107     return "sas_%s%d"%(model_info.id, bits) 
    108  
    109  
    110 def dll_path(model_info, dtype): 
    111     # type: (ModelInfo, np.dtype) -> str 
    112     """ 
    113     Complete path to the dll for the model.  Note that the dll may not 
    114     exist yet if it hasn't been compiled. 
    115     """ 
    116     return os.path.join(DLL_PATH, dll_name(model_info, dtype)+".so") 
    117  
    118  
    119 def make_dll(source, model_info, dtype=F64): 
    120     # type: (str, ModelInfo, np.dtype) -> str 
    121     """ 
    122     Returns the path to the compiled model defined by *kernel_module*. 
    123  
    124     If the model has not been compiled, or if the source file(s) are newer 
    125     than the dll, then *make_dll* will compile the model before returning. 
    126     This routine does not load the resulting dll. 
     92def dll_path(model_info, dtype="double"): 
     93    """ 
     94    Path to the compiled model defined by *model_info*. 
     95    """ 
     96    from os.path import join as joinpath, split as splitpath, splitext 
     97    basename = splitext(splitpath(model_info['filename'])[1])[0] 
     98    if np.dtype(dtype) == generate.F32: 
     99        basename += "32" 
     100    elif np.dtype(dtype) == generate.F64: 
     101        basename += "64" 
     102    else: 
     103        basename += "128" 
     104    return joinpath(DLL_PATH, basename+'.so') 
     105 
     106 
     107def make_dll(source, model_info, dtype="double"): 
     108    """ 
     109    Load the compiled model defined by *kernel_module*. 
     110 
     111    Recompile if any files are newer than the model file. 
    127112 
    128113    *dtype* is a numpy floating point precision specifier indicating whether 
    129     the model should be single, double or long double precision.  The default 
    130     is double precision, *np.dtype('d')*. 
    131  
    132     Set *sasmodels.ALLOW_SINGLE_PRECISION_DLLS* to False if single precision 
    133     models are not allowed as DLLs. 
     114    the model should be single or double precision.  The default is double 
     115    precision. 
     116 
     117    The DLL is not loaded until the kernel is called so models can 
     118    be defined without using too many resources. 
    134119 
    135120    Set *sasmodels.kerneldll.DLL_PATH* to the compiled dll output path. 
    136121    The default is the system temporary directory. 
    137     """ 
    138     if dtype == F16: 
     122 
     123    Set *sasmodels.ALLOW_SINGLE_PRECISION_DLLS* to True if single precision 
     124    models are allowed as DLLs. 
     125    """ 
     126    if callable(model_info.get('Iq', None)): 
     127        return PyModel(model_info) 
     128     
     129    dtype = np.dtype(dtype) 
     130    if dtype == generate.F16: 
    139131        raise ValueError("16 bit floats not supported") 
    140     if dtype == F32 and not ALLOW_SINGLE_PRECISION_DLLS: 
    141         dtype = F64  # Force 64-bit dll 
    142     # Note: dtype may be F128 for long double precision 
    143  
    144     newest = generate.timestamp(model_info) 
     132    if dtype == generate.F32 and not ALLOW_SINGLE_PRECISION_DLLS: 
     133        dtype = generate.F64  # Force 64-bit dll 
     134 
     135    if dtype == generate.F32: # 32-bit dll 
     136        tempfile_prefix = 'sas_' + model_info['name'] + '32_' 
     137    elif dtype == generate.F64: 
     138        tempfile_prefix = 'sas_' + model_info['name'] + '64_' 
     139    else: 
     140        tempfile_prefix = 'sas_' + model_info['name'] + '128_' 
     141  
     142    source = generate.convert_type(source, dtype) 
     143    source_files = generate.model_sources(model_info) + [model_info['filename']] 
    145144    dll = dll_path(model_info, dtype) 
     145    newest = max(os.path.getmtime(f) for f in source_files) 
    146146    if not os.path.exists(dll) or os.path.getmtime(dll) < newest: 
    147         basename = dll_name(model_info, dtype) + "_" 
    148         fid, filename = tempfile.mkstemp(suffix=".c", prefix=basename) 
    149         source = generate.convert_type(source, dtype) 
     147        # Replace with a proper temp file 
     148        fid, filename = tempfile.mkstemp(suffix=".c", prefix=tempfile_prefix) 
    150149        os.fdopen(fid, "w").write(source) 
    151150        command = COMPILE%{"source":filename, "output":dll} 
    … …  
    161160 
    162161 
    163 def load_dll(source, model_info, dtype=F64): 
    164     # type: (str, ModelInfo, np.dtype) -> "DllModel" 
     162def load_dll(source, model_info, dtype="double"): 
    165163    """ 
    166164    Create and load a dll corresponding to the source, info pair returned 
    … …  
    174172 
    175173 
    176 class DllModel(KernelModel): 
     174IQ_ARGS = [c_void_p, c_void_p, c_int] 
     175IQXY_ARGS = [c_void_p, c_void_p, c_void_p, c_int] 
     176 
     177class DllModel(object): 
    177178    """ 
    178179    ctypes wrapper for a single model. 
    … …  
    190191     
    191192    def __init__(self, dllpath, model_info, dtype=generate.F32): 
    192         # type: (str, ModelInfo, np.dtype) -> None 
    193193        self.info = model_info 
    194194        self.dllpath = dllpath 
    195         self._dll = None  # type: ct.CDLL 
     195        self.dll = None 
    196196        self.dtype = np.dtype(dtype) 
    197197 
    198198    def _load_dll(self): 
    199         # type: () -> None 
     199        Nfixed1d = len(self.info['partype']['fixed-1d']) 
     200        Nfixed2d = len(self.info['partype']['fixed-2d']) 
     201        Npd1d = len(self.info['partype']['pd-1d']) 
     202        Npd2d = len(self.info['partype']['pd-2d']) 
     203 
    200204        #print("dll", self.dllpath) 
    201205        try: 
    202             self._dll = ct.CDLL(self.dllpath) 
     206            self.dll = ct.CDLL(self.dllpath) 
    203207        except: 
    204208            annotate_exception("while loading "+self.dllpath) 
    … …  
    208212              else c_double if self.dtype == generate.F64 
    209213              else c_longdouble) 
    210  
    211         # int, int, int, int*, double*, double*, double*, double*, double*, double 
    212         argtypes = [c_int32]*3 + [c_void_p]*5 + [fp] 
    213         self._Iq = self._dll[generate.kernel_name(self.info, is_2d=False)] 
    214         self._Iqxy = self._dll[generate.kernel_name(self.info, is_2d=True)] 
    215         self._Iq.argtypes = argtypes 
    216         self._Iqxy.argtypes = argtypes 
     214        pd_args_1d = [c_void_p, fp] + [c_int]*Npd1d if Npd1d else [] 
     215        pd_args_2d = [c_void_p, fp] + [c_int]*Npd2d if Npd2d else [] 
     216        self.Iq = self.dll[generate.kernel_name(self.info, False)] 
     217        self.Iq.argtypes = IQ_ARGS + pd_args_1d + [fp]*Nfixed1d 
     218 
     219        self.Iqxy = self.dll[generate.kernel_name(self.info, True)] 
     220        self.Iqxy.argtypes = IQXY_ARGS + pd_args_2d + [fp]*Nfixed2d 
     221         
     222        self.release() 
    217223 
    218224    def __getstate__(self): 
    219         # type: () -> Tuple[ModelInfo, str] 
    220225        return self.info, self.dllpath 
    221226 
    222227    def __setstate__(self, state): 
    223         # type: (Tuple[ModelInfo, str]) -> None 
    224228        self.info, self.dllpath = state 
    225         self._dll = None 
     229        self.dll = None 
    226230 
    227231    def make_kernel(self, q_vectors): 
    228         # type: (List[np.ndarray]) -> DllKernel 
    229232        q_input = PyInput(q_vectors, self.dtype) 
    230         # Note: pickle not supported for DllKernel 
    231         if self._dll is None: 
    232             self._load_dll() 
    233         kernel = self._Iqxy if q_input.is_2d else self._Iq 
     233        if self.dll is None: self._load_dll() 
     234        kernel = self.Iqxy if q_input.is_2d else self.Iq 
    234235        return DllKernel(kernel, self.info, q_input) 
    235236 
    236237    def release(self): 
    237         # type: () -> None 
    238238        """ 
    239239        Release any resources associated with the model. 
    … …  
    244244            libHandle = dll._handle 
    245245            #libHandle = ct.c_void_p(dll._handle) 
    246             del dll, self._dll 
    247             self._dll = None 
     246            del dll, self.dll 
     247            self.dll = None 
    248248            #_ctypes.FreeLibrary(libHandle) 
    249249            ct.windll.kernel32.FreeLibrary(libHandle) 
    … …  
    252252 
    253253 
    254 class DllKernel(Kernel): 
     254class DllKernel(object): 
    255255    """ 
    256256    Callable SAS kernel. 
    … …  
    272272    """ 
    273273    def __init__(self, kernel, model_info, q_input): 
    274         # type: (Callable[[], np.ndarray], ModelInfo, PyInput) -> None 
    275         self.kernel = kernel 
    276274        self.info = model_info 
    277275        self.q_input = q_input 
    278         self.dtype = q_input.dtype 
    279         self.dim = '2d' if q_input.is_2d else '1d' 
    280         self.result = np.empty(q_input.nq+1, q_input.dtype) 
    281         self.real = (np.float32 if self.q_input.dtype == generate.F32 
    282                      else np.float64 if self.q_input.dtype == generate.F64 
    283                      else np.float128) 
    284  
    285     def __call__(self, call_details, weights, values, cutoff): 
    286         # type: (CallDetails, np.ndarray, np.ndarray, float) -> np.ndarray 
    287  
    288         #print("in kerneldll") 
    289         #print("weights", weights) 
    290         #print("values", values) 
    291         start, stop = 0, call_details.total_pd 
    292         args = [ 
    293             self.q_input.nq, # nq 
    294             start, # pd_start 
    295             stop, # pd_stop pd_stride[MAX_PD] 
    296             call_details.ctypes.data, # problem 
    297             weights.ctypes.data,  # weights 
    298             values.ctypes.data,  #pars 
    299             self.q_input.q.ctypes.data, #q 
    300             self.result.ctypes.data,   # results 
    301             self.real(cutoff), # cutoff 
    302             ] 
    303         self.kernel(*args) # type: ignore 
    304         return self.result[:-1] 
     276        self.kernel = kernel 
     277        self.res = np.empty(q_input.nq, q_input.dtype) 
     278        dim = '2d' if q_input.is_2d else '1d' 
     279        self.fixed_pars = model_info['partype']['fixed-' + dim] 
     280        self.pd_pars = model_info['partype']['pd-' + dim] 
     281 
     282        # In dll kernel, but not in opencl kernel 
     283        self.p_res = self.res.ctypes.data 
     284 
     285    def __call__(self, fixed_pars, pd_pars, cutoff): 
     286        real = (np.float32 if self.q_input.dtype == generate.F32 
     287                else np.float64 if self.q_input.dtype == generate.F64 
     288                else np.float128) 
     289 
     290        nq = c_int(self.q_input.nq) 
     291        if pd_pars: 
     292            cutoff = real(cutoff) 
     293            loops_N = [np.uint32(len(p[0])) for p in pd_pars] 
     294            loops = np.hstack(pd_pars) 
     295            loops = np.ascontiguousarray(loops.T, self.q_input.dtype).flatten() 
     296            p_loops = loops.ctypes.data 
     297            dispersed = [p_loops, cutoff] + loops_N 
     298        else: 
     299            dispersed = [] 
     300        fixed = [real(p) for p in fixed_pars] 
     301        args = self.q_input.q_pointers + [self.p_res, nq] + dispersed + fixed 
     302        #print(pars) 
     303        self.kernel(*args) 
     304 
     305        return self.res 
    305306 
    306307    def release(self): 
    307         # type: () -> None 
    308308        """ 
    309309        Release any resources associated with the kernel. 
    310310        """ 
    311         self.q_input.release() 
     311        pass 

  • sasmodels/kernelpy.py

    r7ae2b7f r4d76711  
    77:class:`kernelcl.GpuModel` and :class:`kerneldll.DllModel`. 
    88""" 
    9 import numpy as np  # type: ignore 
    10 from numpy import pi, cos  #type: ignore 
    11  
    12 from . import details 
     9import numpy as np 
     10from numpy import pi, cos 
     11 
    1312from .generate import F64 
    14 from .kernel import KernelModel, Kernel 
    15  
    16 try: 
    17     from typing import Union, Callable 
    18 except: 
    19     pass 
    20 else: 
    21     DType = Union[None, str, np.dtype] 
    22  
    23 class PyModel(KernelModel): 
     13 
     14class PyModel(object): 
    2415    """ 
    2516    Wrapper for pure python models. 
    2617    """ 
    2718    def __init__(self, model_info): 
    28         # Make sure Iq and Iqxy are available and vectorized 
    29         _create_default_functions(model_info) 
    3019        self.info = model_info 
    3120 
    3221    def make_kernel(self, q_vectors): 
    3322        q_input = PyInput(q_vectors, dtype=F64) 
    34         kernel = self.info.Iqxy if q_input.is_2d else self.info.Iq 
     23        kernel = self.info['Iqxy'] if q_input.is_2d else self.info['Iq'] 
    3524        return PyKernel(kernel, self.info, q_input) 
    3625 
    … …  
    6453        self.dtype = dtype 
    6554        self.is_2d = (len(q_vectors) == 2) 
    66         if self.is_2d: 
    67             self.q = np.empty((self.nq, 2), dtype=dtype) 
    68             self.q[:, 0] = q_vectors[0] 
    69             self.q[:, 1] = q_vectors[1] 
    70         else: 
    71             self.q = np.empty(self.nq, dtype=dtype) 
    72             self.q[:self.nq] = q_vectors[0] 
     55        self.q_vectors = [np.ascontiguousarray(q, self.dtype) for q in q_vectors] 
     56        self.q_pointers = [q.ctypes.data for q in self.q_vectors] 
    7357 
    7458    def release(self): 
    … …  
    7660        Free resources associated with the model inputs. 
    7761        """ 
    78         self.q = None 
    79  
    80 class PyKernel(Kernel): 
     62        self.q_vectors = [] 
     63 
     64class PyKernel(object): 
    8165    """ 
    8266    Callable SAS kernel. 
    … …  
    9882    """ 
    9983    def __init__(self, kernel, model_info, q_input): 
    100         self.dtype = np.dtype('d') 
    10184        self.info = model_info 
    10285        self.q_input = q_input 
    10386        self.res = np.empty(q_input.nq, q_input.dtype) 
    104         self.kernel = kernel 
    105         self.dim = '2d' if q_input.is_2d else '1d' 
    106  
    107         partable = model_info.parameters 
    108         kernel_parameters = (partable.iqxy_parameters if q_input.is_2d 
    109                              else partable.iq_parameters) 
    110         volume_parameters = partable.form_volume_parameters 
    111  
    112         # Create an array to hold the parameter values.  There will be a 
    113         # single array whose values are updated as the calculator goes 
    114         # through the loop.  Arguments to the kernel and volume functions 
    115         # will use views into this vector, relying on the fact that a 
    116         # an array of no dimensions acts like a scalar. 
    117         parameter_vector = np.empty(len(partable.call_parameters)-2, 'd') 
    118  
    119         # Create views into the array to hold the arguments 
    120         offset = 0 
    121         kernel_args, volume_args = [], [] 
    122         for p in partable.kernel_parameters: 
    123             if p.length == 1: 
    124                 # Scalar values are length 1 vectors with no dimensions. 
    125                 v = parameter_vector[offset:offset+1].reshape(()) 
     87        dim = '2d' if q_input.is_2d else '1d' 
     88        # Loop over q unless user promises that the kernel is vectorized by 
     89        # taggining it with vectorized=True 
     90        if not getattr(kernel, 'vectorized', False): 
     91            if dim == '2d': 
     92                def vector_kernel(qx, qy, *args): 
     93                    """ 
     94                    Vectorized 2D kernel. 
     95                    """ 
     96                    return np.array([kernel(qxi, qyi, *args) 
     97                                     for qxi, qyi in zip(qx, qy)]) 
    12698            else: 
    127                 # Vector values are simple views. 
    128                 v = parameter_vector[offset:offset+p.length] 
    129             offset += p.length 
    130             if p in kernel_parameters: 
    131                 kernel_args.append(v) 
    132             if p in volume_parameters: 
    133                 volume_args.append(v) 
    134  
    135         # Hold on to the parameter vector so we can use it to call kernel later. 
    136         # This may also be required to preserve the views into the vector. 
    137         self._parameter_vector = parameter_vector 
    138  
    139         # Generate a closure which calls the kernel with the views into the 
    140         # parameter array. 
    141         if q_input.is_2d: 
    142             form = model_info.Iqxy 
    143             qx, qy = q_input.q[:,0], q_input.q[:,1] 
    144             self._form = lambda: form(qx, qy, *kernel_args) 
     99                def vector_kernel(q, *args): 
     100                    """ 
     101                    Vectorized 1D kernel. 
     102                    """ 
     103                    return np.array([kernel(qi, *args) 
     104                                     for qi in q]) 
     105            self.kernel = vector_kernel 
    145106        else: 
    146             form = model_info.Iq 
    147             q = q_input.q 
    148             self._form = lambda: form(q, *kernel_args) 
    149  
    150         # Generate a closure which calls the form_volume if it exists. 
    151         form_volume = model_info.form_volume 
    152         self._volume = ((lambda: form_volume(*volume_args)) if form_volume 
    153                         else (lambda: 1.0)) 
    154  
    155     def __call__(self, call_details, weights, values, cutoff): 
    156         assert isinstance(call_details, details.CallDetails) 
    157         res = _loops(self._parameter_vector, self._form, self._volume, 
    158                      self.q_input.nq, call_details, weights, values, cutoff) 
     107            self.kernel = kernel 
     108        fixed_pars = model_info['partype']['fixed-' + dim] 
     109        pd_pars = model_info['partype']['pd-' + dim] 
     110        vol_pars = model_info['partype']['volume'] 
     111 
     112        # First two fixed pars are scale and background 
     113        pars = [p.name for p in model_info['parameters'][2:]] 
     114        offset = len(self.q_input.q_vectors) 
     115        self.args = self.q_input.q_vectors + [None] * len(pars) 
     116        self.fixed_index = np.array([pars.index(p) + offset 
     117                                     for p in fixed_pars[2:]]) 
     118        self.pd_index = np.array([pars.index(p) + offset 
     119                                  for p in pd_pars]) 
     120        self.vol_index = np.array([pars.index(p) + offset 
     121                                   for p in vol_pars]) 
     122        try: self.theta_index = pars.index('theta') + offset 
     123        except ValueError: self.theta_index = -1 
     124 
     125        # Caller needs fixed_pars and pd_pars 
     126        self.fixed_pars = fixed_pars 
     127        self.pd_pars = pd_pars 
     128 
     129    def __call__(self, fixed, pd, cutoff=1e-5): 
     130        #print("fixed",fixed) 
     131        #print("pd", pd) 
     132        args = self.args[:]  # grab a copy of the args 
     133        form, form_volume = self.kernel, self.info['form_volume'] 
     134        # First two fixed 
     135        scale, background = fixed[:2] 
     136        for index, value in zip(self.fixed_index, fixed[2:]): 
     137            args[index] = float(value) 
     138        res = _loops(form, form_volume, cutoff, scale, background, args, 
     139                     pd, self.pd_index, self.vol_index, self.theta_index) 
     140 
    159141        return res 
    160142 
    … …  
    163145        Free resources associated with the kernel. 
    164146        """ 
    165         self.q_input.release() 
    166147        self.q_input = None 
    167148 
    168 def _loops(parameters, form, form_volume, nq, call_details, 
    169            weights, values, cutoff): 
    170     # type: (np.ndarray, Callable[[], np.ndarray], Callable[[], float], int, details.CallDetails, np.ndarray, np.ndarray, float) -> None 
     149def _loops(form, form_volume, cutoff, scale, background, 
     150           args, pd, pd_index, vol_index, theta_index): 
     151    """ 
     152    *form* is the name of the form function, which should be vectorized over 
     153    q, but otherwise have an interface like the opencl kernels, with the 
     154    q parameters first followed by the individual parameters in the order 
     155    given in model.parameters (see :mod:`sasmodels.generate`). 
     156 
     157    *form_volume* calculates the volume of the shape.  *vol_index* gives 
     158    the list of volume parameters 
     159 
     160    *cutoff* ignores the corners of the dispersion hypercube 
     161 
     162    *scale*, *background* multiplies the resulting form and adds an offset 
     163 
     164    *args* is the prepopulated set of arguments to the form function, starting 
     165    with the q vectors, and including slots for all the parameters.  The 
     166    values for the parameters will be substituted with values from the 
     167    dispersion functions. 
     168 
     169    *pd* is the list of dispersion parameters 
     170 
     171    *pd_index* are the indices of the dispersion parameters in *args* 
     172 
     173    *vol_index* are the indices of the volume parameters in *args* 
     174 
     175    *theta_index* is the index of the theta parameter for the sasview 
     176    spherical correction, or -1 if there is no angular dispersion 
     177    """ 
     178 
    171179    ################################################################ 
    172180    #                                                              # 
    … …  
    178186    #                                                              # 
    179187    ################################################################ 
    180     parameters[:] = values[call_details.par_offset] 
    181     scale, background = values[0], values[1] 
    182     if call_details.num_active == 0: 
    183         norm = float(form_volume()) 
    184         if norm > 0.0: 
    185             return (scale/norm)*form() + background 
     188 
     189    weight = np.empty(len(pd), 'd') 
     190    if weight.size > 0: 
     191        # weight vector, to be populated by polydispersity loops 
     192 
     193        # identify which pd parameters are volume parameters 
     194        vol_weight_index = np.array([(index in vol_index) for index in pd_index]) 
     195 
     196        # Sort parameters in decreasing order of pd length 
     197        Npd = np.array([len(pdi[0]) for pdi in pd], 'i') 
     198        order = np.argsort(Npd)[::-1] 
     199        stride = np.cumprod(Npd[order]) 
     200        pd = [pd[index] for index in order] 
     201        pd_index = pd_index[order] 
     202        vol_weight_index = vol_weight_index[order] 
     203 
     204        fast_value = pd[0][0] 
     205        fast_weight = pd[0][1] 
     206    else: 
     207        stride = np.array([1]) 
     208        vol_weight_index = slice(None, None) 
     209        # keep lint happy 
     210        fast_value = [None] 
     211        fast_weight = [None] 
     212 
     213    ret = np.zeros_like(args[0]) 
     214    norm = np.zeros_like(ret) 
     215    vol = np.zeros_like(ret) 
     216    vol_norm = np.zeros_like(ret) 
     217    for k in range(stride[-1]): 
     218        # update polydispersity parameter values 
     219        fast_index = k % stride[0] 
     220        if fast_index == 0:  # bottom loop complete ... check all other loops 
     221            if weight.size > 0: 
     222                for i, index, in enumerate(k % stride): 
     223                    args[pd_index[i]] = pd[i][0][index] 
     224                    weight[i] = pd[i][1][index] 
    186225        else: 
    187             return np.ones(nq, 'd')*background 
    188  
    189     partial_weight = np.NaN 
    190     spherical_correction = 1.0 
    191     pd_stride = call_details.pd_stride[:call_details.num_active] 
    192     pd_length = call_details.pd_length[:call_details.num_active] 
    193     pd_offset = call_details.pd_offset[:call_details.num_active] 
    194     pd_index = np.empty_like(pd_offset) 
    195     offset = np.empty_like(call_details.par_offset) 
    196     theta = call_details.theta_par 
    197     fast_length = pd_length[0] 
    198     pd_index[0] = fast_length 
    199     total = np.zeros(nq, 'd') 
    200     norm = 0.0 
    201     for loop_index in range(call_details.total_pd): 
    202         # update polydispersity parameter values 
    203         if pd_index[0] == fast_length: 
    204             pd_index[:] = (loop_index/pd_stride)%pd_length 
    205             partial_weight = np.prod(weights[pd_offset+pd_index][1:]) 
    206             for k in range(call_details.num_coord): 
    207                 par = call_details.par_coord[k] 
    208                 coord = call_details.pd_coord[k] 
    209                 this_offset = call_details.par_offset[par] 
    210                 block_size = 1 
    211                 for bit in range(len(pd_offset)): 
    212                     if coord&1: 
    213                         this_offset += block_size * pd_index[bit] 
    214                         block_size *= pd_length[bit] 
    215                     coord >>= 1 
    216                     if coord == 0: break 
    217                 offset[par] = this_offset 
    218                 parameters[par] = values[this_offset] 
    219                 if par == theta and not (call_details.par_coord[k]&1): 
    220                     spherical_correction = max(abs(cos(pi/180 * parameters[theta])), 1e-6) 
    221         for k in range(call_details.num_coord): 
    222             if call_details.pd_coord[k]&1: 
    223                 #par = call_details.par_coord[k] 
    224                 parameters[par] = values[offset[par]] 
    225                 #print "par",par,offset[par],parameters[par+2] 
    226                 offset[par] += 1 
    227                 if par == theta: 
    228                     spherical_correction = max(abs(cos(pi/180 * parameters[theta])), 1e-6) 
    229  
    230         weight = partial_weight * weights[pd_offset[0] + pd_index[0]] 
    231         pd_index[0] += 1 
    232         if weight > cutoff: 
    233             # Call the scattering function 
    234             # Assume that NaNs are only generated if the parameters are bad; 
    235             # exclude all q for that NaN.  Even better would be to have an 
    236             # INVALID expression like the C models, but that is too expensive. 
    237             I = form() 
     226            args[pd_index[0]] = fast_value[fast_index] 
     227            weight[0] = fast_weight[fast_index] 
     228 
     229        # Computes the weight, and if it is not sufficient then ignore this 
     230        # parameter set. 
     231        # Note: could precompute w1*...*wn so we only need to multiply by w0 
     232        w = np.prod(weight) 
     233        if w > cutoff: 
     234            # Note: can precompute spherical correction if theta_index is not 
     235            # the fast index. Correction factor for spherical integration 
     236            #spherical_correction = abs(cos(pi*args[phi_index])) if phi_index>=0 else 1.0 
     237            spherical_correction = (abs(cos(pi * args[theta_index])) * pi / 2 
     238                                    if theta_index >= 0 else 1.0) 
     239            #spherical_correction = 1.0 
     240 
     241            # Call the scattering function and adds it to the total. 
     242            I = form(*args) 
    238243            if np.isnan(I).any(): continue 
    239  
    240             # update value and norm 
    241             weight *= spherical_correction 
    242             total += weight * I 
    243             norm += weight * form_volume() 
    244  
    245     if norm > 0.0: 
    246         return (scale/norm)*total + background 
    247     else: 
    248         return np.ones(nq, 'd')*background 
    249  
    250  
    251 def _create_default_functions(model_info): 
    252     """ 
    253     Autogenerate missing functions, such as Iqxy from Iq. 
    254  
    255     This only works for Iqxy when Iq is written in python. :func:`make_source` 
    256     performs a similar role for Iq written in C.  This also vectorizes 
    257     any functions that are not already marked as vectorized. 
    258     """ 
    259     _create_vector_Iq(model_info) 
    260     _create_vector_Iqxy(model_info)  # call create_vector_Iq() first 
    261  
    262  
    263 def _create_vector_Iq(model_info): 
    264     """ 
    265     Define Iq as a vector function if it exists. 
    266     """ 
    267     Iq = model_info.Iq 
    268     if callable(Iq) and not getattr(Iq, 'vectorized', False): 
    269         #print("vectorizing Iq") 
    270         def vector_Iq(q, *args): 
    271             """ 
    272             Vectorized 1D kernel. 
    273             """ 
    274             return np.array([Iq(qi, *args) for qi in q]) 
    275         vector_Iq.vectorized = True 
    276         model_info.Iq = vector_Iq 
    277  
    278 def _create_vector_Iqxy(model_info): 
    279     """ 
    280     Define Iqxy as a vector function if it exists, or default it from Iq(). 
    281     """ 
    282     Iq, Iqxy = model_info.Iq, model_info.Iqxy 
    283     if callable(Iqxy): 
    284         if not getattr(Iqxy, 'vectorized', False): 
    285             #print("vectorizing Iqxy") 
    286             def vector_Iqxy(qx, qy, *args): 
    287                 """ 
    288                 Vectorized 2D kernel. 
    289                 """ 
    290                 return np.array([Iqxy(qxi, qyi, *args) for qxi, qyi in zip(qx, qy)]) 
    291             vector_Iqxy.vectorized = True 
    292             model_info.Iqxy = vector_Iqxy 
    293     elif callable(Iq): 
    294         #print("defaulting Iqxy") 
    295         # Iq is vectorized because create_vector_Iq was already called. 
    296         def default_Iqxy(qx, qy, *args): 
    297             """ 
    298             Default 2D kernel. 
    299             """ 
    300             return Iq(np.sqrt(qx**2 + qy**2), *args) 
    301         default_Iqxy.vectorized = True 
    302         model_info.Iqxy = default_Iqxy 
    303  
     244            ret += w * I * spherical_correction 
     245            norm += w 
     246 
     247            # Volume normalization. 
     248            # If there are "volume" polydispersity parameters, then these 
     249            # will be used to call the form_volume function from the user 
     250            # supplied kernel, and accumulate a normalized weight. 
     251            # Note: can precompute volume norm if fast index is not a volume 
     252            if form_volume: 
     253                vol_args = [args[index] for index in vol_index] 
     254                vol_weight = np.prod(weight[vol_weight_index]) 
     255                vol += vol_weight * form_volume(*vol_args) 
     256                vol_norm += vol_weight 
     257 
     258    positive = (vol * vol_norm != 0.0) 
     259    ret[positive] *= vol_norm[positive] / vol[positive] 
     260    result = scale * ret / norm + background 
     261    return result 

  • sasmodels/list_pars.py

    r6d6508e ra84a0ca  
    2525    for name in sorted(MODELS): 
    2626        model_info = load_model_info(name) 
    27         for p in model_info.parameters.kernel_parameters: 
     27        for p in model_info['parameters']: 
    2828            partable.setdefault(p.name, []) 
    2929            partable[p.name].append(name) 

  • sasmodels/mixture.py

    r7ae2b7f r72a081d  
    1212""" 
    1313from copy import copy 
    14 import numpy as np  # type: ignore 
     14import numpy as np 
    1515 
    16 from .modelinfo import Parameter, ParameterTable, ModelInfo 
    17 from .kernel import KernelModel, Kernel 
     16from .generate import process_parameters, COMMON_PARAMETERS, Parameter 
    1817 
    19 try: 
    20     from typing import List 
    21     from .details import CallDetails 
    22 except ImportError: 
    23     pass 
     18SCALE=0 
     19BACKGROUND=1 
     20EFFECT_RADIUS=2 
     21VOLFRACTION=3 
    2422 
    2523def make_mixture_info(parts): 
    26     # type: (List[ModelInfo]) -> ModelInfo 
    2724    """ 
    2825    Create info block for product model. 
    … …  
    3027    flatten = [] 
    3128    for part in parts: 
    32         if part.composition and part.composition[0] == 'mixture': 
    33             flatten.extend(part.composition[1]) 
     29        if part['composition'] and part['composition'][0] == 'mixture': 
     30            flatten.extend(part['compostion'][1]) 
    3431        else: 
    3532            flatten.append(part) 
    … …  
    3734 
    3835    # Build new parameter list 
    39     combined_pars = [] 
     36    pars = COMMON_PARAMETERS + [] 
    4037    for k, part in enumerate(parts): 
    4138        # Parameter prefix per model, A_, B_, ... 
    42         # Note that prefix must also be applied to id and length_control 
    43         # to support vector parameters 
    4439        prefix = chr(ord('A')+k) + '_' 
    45         combined_pars.append(Parameter(prefix+'scale')) 
    46         for p in part.parameters.kernel_parameters: 
    47             p = copy(p) 
    48             p.name = prefix + p.name 
    49             p.id = prefix + p.id 
    50             if p.length_control is not None: 
    51                 p.length_control = prefix + p.length_control 
    52             combined_pars.append(p) 
    53     parameters = ParameterTable(combined_pars) 
     40        for p in part['parameters']: 
     41            # No background on the individual mixture elements 
     42            if p.name == 'background': 
     43                continue 
     44            # TODO: promote Parameter to a full class 
     45            # this code knows too much about the implementation! 
     46            p = list(p) 
     47            if p[0] == 'scale':  # allow model subtraction 
     48                p[3] = [-np.inf, np.inf]  # limits 
     49            p[0] = prefix+p[0]   # relabel parameters with A_, ... 
     50            pars.append(p) 
    5451 
    55     model_info = ModelInfo() 
    56     model_info.id = '+'.join(part.id for part in parts) 
    57     model_info.name = ' + '.join(part.name for part in parts) 
    58     model_info.filename = None 
    59     model_info.title = 'Mixture model with ' + model_info.name 
    60     model_info.description = model_info.title 
    61     model_info.docs = model_info.title 
    62     model_info.category = "custom" 
    63     model_info.parameters = parameters 
    64     #model_info.single = any(part['single'] for part in parts) 
    65     model_info.structure_factor = False 
    66     model_info.variant_info = None 
    67     #model_info.tests = [] 
    68     #model_info.source = [] 
     52    model_info = {} 
     53    model_info['id'] = '+'.join(part['id']) 
     54    model_info['name'] = ' + '.join(part['name']) 
     55    model_info['filename'] = None 
     56    model_info['title'] = 'Mixture model with ' + model_info['name'] 
     57    model_info['description'] = model_info['title'] 
     58    model_info['docs'] = model_info['title'] 
     59    model_info['category'] = "custom" 
     60    model_info['parameters'] = pars 
     61    #model_info['single'] = any(part['single'] for part in parts) 
     62    model_info['structure_factor'] = False 
     63    model_info['variant_info'] = None 
     64    #model_info['tests'] = [] 
     65    #model_info['source'] = [] 
    6966    # Iq, Iqxy, form_volume, ER, VR and sesans 
    7067    # Remember the component info blocks so we can build the model 
    71     model_info.composition = ('mixture', parts) 
     68    model_info['composition'] = ('mixture', parts) 
     69    process_parameters(model_info) 
     70    return model_info 
    7271 
    7372 
    74 class MixtureModel(KernelModel): 
     73class MixtureModel(object): 
    7574    def __init__(self, model_info, parts): 
    76         # type: (ModelInfo, List[KernelModel]) -> None 
    7775        self.info = model_info 
    7876        self.parts = parts 
    7977 
    8078    def __call__(self, q_vectors): 
    81         # type: (List[np.ndarray]) -> MixtureKernel 
    8279        # Note: may be sending the q_vectors to the n times even though they 
    8380        # are only needed once.  It would mess up modularity quite a bit to 
    … …  
    8683        # in opencl; or both in opencl, but one in single precision and the 
    8784        # other in double precision). 
    88         kernels = [part.make_kernel(q_vectors) for part in self.parts] 
     85        kernels = [part(q_vectors) for part in self.parts] 
    8986        return MixtureKernel(self.info, kernels) 
    9087 
    9188    def release(self): 
    92         # type: () -> None 
    9389        """ 
    9490        Free resources associated with the model. 
    … …  
    9894 
    9995 
    100 class MixtureKernel(Kernel): 
     96class MixtureKernel(object): 
    10197    def __init__(self, model_info, kernels): 
    102         # type: (ModelInfo, List[Kernel]) -> None 
    103         self.dim = kernels[0].dim 
    104         self.info =  model_info 
     98        dim = '2d' if kernels[0].q_input.is_2d else '1d' 
     99 
     100        # fixed offsets starts at 2 for scale and background 
     101        fixed_pars, pd_pars = [], [] 
     102        offsets = [[2, 0]] 
     103        #vol_index = [] 
     104        def accumulate(fixed, pd, volume): 
     105            # subtract 1 from fixed since we are removing background 
     106            fixed_offset, pd_offset = offsets[-1] 
     107            #vol_index.extend(k+pd_offset for k,v in pd if v in volume) 
     108            offsets.append([fixed_offset + len(fixed) - 1, pd_offset + len(pd)]) 
     109            pd_pars.append(pd) 
     110        if dim == '2d': 
     111            for p in kernels: 
     112                partype = p.info['partype'] 
     113                accumulate(partype['fixed-2d'], partype['pd-2d'], partype['volume']) 
     114        else: 
     115            for p in kernels: 
     116                partype = p.info['partype'] 
     117                accumulate(partype['fixed-1d'], partype['pd-1d'], partype['volume']) 
     118 
     119        #self.vol_index = vol_index 
     120        self.offsets = offsets 
     121        self.fixed_pars = fixed_pars 
     122        self.pd_pars = pd_pars 
     123        self.info = model_info 
    105124        self.kernels = kernels 
     125        self.results = None 
    106126 
    107     def __call__(self, call_details, value, weight, cutoff): 
    108         # type: (CallDetails, np.ndarray, np.ndarry, float) -> np.ndarray 
    109         scale, background = value[0:2] 
     127    def __call__(self, fixed_pars, pd_pars, cutoff=1e-5): 
     128        scale, background = fixed_pars[0:2] 
    110129        total = 0.0 
    111         # remember the parts for plotting later 
    112         self.results = [] 
    113         for kernel, kernel_details in zip(self.kernels, call_details.parts): 
    114             part_result = kernel(kernel_details, value, weight, cutoff) 
     130        self.results = []  # remember the parts for plotting later 
     131        for k in range(len(self.offsets)-1): 
     132            start_fixed, start_pd = self.offsets[k] 
     133            end_fixed, end_pd = self.offsets[k+1] 
     134            part_fixed = [fixed_pars[start_fixed], 0.0] + fixed_pars[start_fixed+1:end_fixed] 
     135            part_pd = [pd_pars[start_pd], 0.0] + pd_pars[start_pd+1:end_pd] 
     136            part_result = self.kernels[k](part_fixed, part_pd) 
    115137            total += part_result 
    116             self.results.append(part_result) 
     138            self.results.append(scale*sum+background) 
    117139 
    118140        return scale*total + background 
    119141 
    120142    def release(self): 
    121         # type: () -> None 
    122         for k in self.kernels: 
    123             k.release() 
     143        self.p_kernel.release() 
     144        self.q_kernel.release() 
    124145 

  • sasmodels/model_test.py

    r7ae2b7f r4d76711  
    4343Precision defaults to 5 digits (relative). 
    4444""" 
    45 #TODO: rename to tests so that tab completion works better for models directory 
    46  
    4745from __future__ import print_function 
    4846 
    … …  
    5048import unittest 
    5149 
    52 import numpy as np  # type: ignore 
     50import numpy as np 
    5351 
    5452from .core import list_models, load_model_info, build_model, HAVE_OPENCL 
    55 from .details import dispersion_mesh 
    56 from .direct_model import call_kernel, get_weights 
     53from .core import call_kernel, call_ER, call_VR 
    5754from .exception import annotate_exception 
    58 from .modelinfo import expand_pars 
    59  
    60 try: 
    61     from typing import List, Iterator, Callable 
    62 except ImportError: 
    63     pass 
    64 else: 
    65     from .modelinfo import ParameterTable, ParameterSet, TestCondition, ModelInfo 
    66     from .kernelpy import PyModel, PyInput, PyKernel, DType 
    67  
    68 def call_ER(model_info, pars): 
    69     # type: (ModelInfo, ParameterSet) -> float 
    70     """ 
    71     Call the model ER function using *values*. 
    72  
    73     *model_info* is either *model.info* if you have a loaded model, 
    74     or *kernel.info* if you have a model kernel prepared for evaluation. 
    75     """ 
    76     if model_info.ER is None: 
    77         return 1.0 
    78     else: 
    79         value, weight = _vol_pars(model_info, pars) 
    80         individual_radii = model_info.ER(*value) 
    81         return np.sum(weight*individual_radii) / np.sum(weight) 
    82  
    83 def call_VR(model_info, pars): 
    84     # type: (ModelInfo, ParameterSet) -> float 
    85     """ 
    86     Call the model VR function using *pars*. 
    87  
    88     *model_info* is either *model.info* if you have a loaded model, 
    89     or *kernel.info* if you have a model kernel prepared for evaluation. 
    90     """ 
    91     if model_info.VR is None: 
    92         return 1.0 
    93     else: 
    94         value, weight = _vol_pars(model_info, pars) 
    95         whole, part = model_info.VR(*value) 
    96         return np.sum(weight*part)/np.sum(weight*whole) 
    97  
    98 def _vol_pars(model_info, pars): 
    99     # type: (ModelInfo, ParameterSet) -> Tuple[np.ndarray, np.ndarray] 
    100     vol_pars = [get_weights(p, pars) 
    101                 for p in model_info.parameters.call_parameters 
    102                 if p.type == 'volume'] 
    103     value, weight = dispersion_mesh(model_info, vol_pars) 
    104     return value, weight 
    105  
     55 
     56#TODO: rename to tests so that tab completion works better for models directory 
    10657 
    10758def make_suite(loaders, models): 
    108     # type: (List[str], List[str]) -> unittest.TestSuite 
    10959    """ 
    11060    Construct the pyunit test suite. 
    … …  
    11666    *models* is the list of models to test, or *["all"]* to test all models. 
    11767    """ 
    118     ModelTestCase = _hide_model_case_from_nose() 
     68 
     69    ModelTestCase = _hide_model_case_from_nosetests() 
    11970    suite = unittest.TestSuite() 
    12071 
    … …  
    13586        # don't try to call cl kernel since it will not be 
    13687        # available in some environmentes. 
    137         is_py = callable(model_info.Iq) 
     88        is_py = callable(model_info['Iq']) 
    13889 
    13990        if is_py:  # kernel implemented in python 
    … …  
    173124 
    174125 
    175 def _hide_model_case_from_nose(): 
    176     # type: () -> type 
     126def _hide_model_case_from_nosetests(): 
    177127    class ModelTestCase(unittest.TestCase): 
    178128        """ 
    … …  
    185135        def __init__(self, test_name, model_info, test_method_name, 
    186136                     platform, dtype): 
    187             # type: (str, ModelInfo, str, str, DType) -> None 
    188137            self.test_name = test_name 
    189138            self.info = model_info 
    … …  
    191140            self.dtype = dtype 
    192141 
    193             setattr(self, test_method_name, self.run_all) 
     142            setattr(self, test_method_name, self._runTest) 
    194143            unittest.TestCase.__init__(self, test_method_name) 
    195144 
    196         def run_all(self): 
    197             # type: () -> None 
     145        def _runTest(self): 
    198146            smoke_tests = [ 
    199                 ({}, 0.1, None), 
    200                 ({}, (0.1, 0.1), None), 
    201                 ({}, 'ER', None), 
    202                 ({}, 'VR', None), 
     147                [{}, 0.1, None], 
     148                [{}, (0.1, 0.1), None], 
     149                [{}, 'ER', None], 
     150                [{}, 'VR', None], 
    203151                ] 
    204152 
    205             tests = self.info.tests 
     153            tests = self.info['tests'] 
    206154            try: 
    207155                model = build_model(self.info, dtype=self.dtype, 
    208156                                    platform=self.platform) 
    209157                for test in smoke_tests + tests: 
    210                     self.run_one(model, test) 
     158                    self._run_one_test(model, test) 
    211159 
    212160                if not tests and self.platform == "dll": 
    … …  
    222170                raise 
    223171 
    224         def run_one(self, model, test): 
    225             # type: (PyModel, TestCondition) -> None 
    226             user_pars, x, y = test 
    227             pars = expand_pars(self.info.parameters, user_pars) 
     172        def _run_one_test(self, model, test): 
     173            pars, x, y = test 
    228174 
    229175            if not isinstance(y, list): 
    … …  
    239185                actual = [call_VR(model.info, pars)] 
    240186            elif isinstance(x[0], tuple): 
    241                 qx, qy = zip(*x) 
    242                 q_vectors = [np.array(qx), np.array(qy)] 
     187                Qx, Qy = zip(*x) 
     188                q_vectors = [np.array(Qx), np.array(Qy)] 
    243189                kernel = model.make_kernel(q_vectors) 
    244190                actual = call_kernel(kernel, pars) 
    … …  
    248194                actual = call_kernel(kernel, pars) 
    249195 
    250             self.assertTrue(len(actual) > 0) 
     196            self.assertGreater(len(actual), 0) 
    251197            self.assertEqual(len(y), len(actual)) 
    252198 
    … …  
    264210 
    265211def is_near(target, actual, digits=5): 
    266     # type: (float, float, int) -> bool 
    267212    """ 
    268213    Returns true if *actual* is within *digits* significant digits of *target*. 
    … …  
    273218 
    274219def main(): 
    275     # type: () -> int 
    276220    """ 
    277221    Run tests given is sys.argv. 
    … …  
    307251  python -m sasmodels.model_test [-v] [opencl|dll] model1 model2 ... 
    308252 
    309 If -v is included on the command line, then use verbose output. 
     253If -v is included on the command line, then use verboe output. 
    310254 
    311255If neither opencl nor dll is specified, then models will be tested with 
    312 both OpenCL and dll; the compute target is ignored for pure python models. 
     256both opencl and dll; the compute target is ignored for pure python models. 
    313257 
    314258If model1 is 'all', then all except the remaining models will be tested. 
    … …  
    325269 
    326270def model_tests(): 
    327     # type: () -> Iterator[Callable[[], None]] 
    328271    """ 
    329272    Test runner visible to nosetests. 
    … …  
    333276    tests = make_suite(['opencl', 'dll'], ['all']) 
    334277    for test_i in tests: 
    335         yield test_i.run_all 
     278        yield test_i._runTest 
    336279 
    337280 

  • sasmodels/models/cylinder.c

    r03cac08 r26141cb  
    33double Iqxy(double qx, double qy, double sld, double solvent_sld, 
    44    double radius, double length, double theta, double phi); 
    5  
    6 #define INVALID(v) (v.radius<0 || v.length<0) 
    75 
    86double form_volume(double radius, double length) 
    … …  
    1715    double length) 
    1816{ 
     17    // TODO: return NaN if radius<0 or length<0? 
    1918    // precompute qr and qh to save time in the loop 
    2019    const double qr = q*radius; 
    … …  
    4847    double phi) 
    4948{ 
     49    // TODO: return NaN if radius<0 or length<0? 
    5050    double sn, cn; // slots to hold sincos function output 
    5151 

  • sasmodels/models/cylinder.py

    r7ae2b7f rf247314  
    8282""" 
    8383 
    84 import numpy as np  # type: ignore 
    85 from numpy import pi, inf  # type: ignore 
     84import numpy as np 
     85from numpy import pi, inf 
    8686 
    8787name = "cylinder" 

  • sasmodels/models/flexible_cylinder.c

    re6408d0 re6408d0  
    1 static double 
    2 form_volume(length, kuhn_length, radius) 
     1double form_volume(double length, double kuhn_length, double radius); 
     2double Iq(double q, double length, double kuhn_length, double radius, 
     3          double sld, double solvent_sld); 
     4double Iqxy(double qx, double qy, double length, double kuhn_length, 
     5            double radius, double sld, double solvent_sld); 
     6double flexible_cylinder_kernel(double q, double length, double kuhn_length, 
     7                                double radius, double sld, double solvent_sld); 
     8 
     9 
     10double form_volume(double length, double kuhn_length, double radius) 
    311{ 
    412    return 1.0; 
    513} 
    614 
    7 static double 
    8 Iq(double q, 
    9    double length, 
    10    double kuhn_length, 
    11    double radius, 
    12    double sld, 
    13    double solvent_sld) 
     15double flexible_cylinder_kernel(double q, 
     16          double length, 
     17          double kuhn_length, 
     18          double radius, 
     19          double sld, 
     20          double solvent_sld) 
    1421{ 
    15     const double contrast = sld - solvent_sld; 
    16     const double cross_section = sas_J1c(q*radius); 
    17     const double volume = M_PI*radius*radius*length; 
    18     const double flex = Sk_WR(q, length, kuhn_length); 
    19     return 1.0e-4 * volume * square(contrast*cross_section) * flex; 
     22 
     23    const double cont = sld-solvent_sld; 
     24    const double qr = q*radius; 
     25    //const double crossSect = (2.0*J1(qr)/qr)*(2.0*J1(qr)/qr); 
     26    const double crossSect = sas_J1c(qr); 
     27    double flex = Sk_WR(q,length,kuhn_length); 
     28    flex *= crossSect*crossSect; 
     29    flex *= M_PI*radius*radius*length; 
     30    flex *= cont*cont; 
     31    flex *= 1.0e-4; 
     32    return flex; 
    2033} 
     34 
     35double Iq(double q, 
     36          double length, 
     37          double kuhn_length, 
     38          double radius, 
     39          double sld, 
     40          double solvent_sld) 
     41{ 
     42 
     43    double result = flexible_cylinder_kernel(q, length, kuhn_length, radius, sld, solvent_sld); 
     44    return result; 
     45} 
     46 
     47double Iqxy(double qx, double qy, 
     48            double length, 
     49            double kuhn_length, 
     50            double radius, 
     51            double sld, 
     52            double solvent_sld) 
     53{ 
     54    double q; 
     55    q = sqrt(qx*qx+qy*qy); 
     56    double result = flexible_cylinder_kernel(q, length, kuhn_length, radius, sld, solvent_sld); 
     57 
     58    return result; 
     59} 

  • sasmodels/models/gel_fit.c

    r03cac08 r30b4ddf  
    1 static double Iq(double q, 
     1double form_volume(void); 
     2 
     3double Iq(double q, 
     4          double guinier_scale, 
     5          double lorentzian_scale, 
     6          double gyration_radius, 
     7          double fractal_exp, 
     8          double cor_length); 
     9 
     10double Iqxy(double qx, double qy, 
     11          double guinier_scale, 
     12          double lorentzian_scale, 
     13          double gyration_radius, 
     14          double fractal_exp, 
     15          double cor_length); 
     16 
     17static double _gel_fit_kernel(double q, 
    218          double guinier_scale, 
    319          double lorentzian_scale, 
    … …  
    824    // Lorentzian Term 
    925    ////////////////////////double a(x[i]*x[i]*zeta*zeta); 
    10     double lorentzian_term = square(q*cor_length); 
     26    double lorentzian_term = q*q*cor_length*cor_length; 
    1127    lorentzian_term = 1.0 + ((fractal_exp + 1.0)/3.0)*lorentzian_term; 
    1228    lorentzian_term = pow(lorentzian_term, (fractal_exp/2.0) ); 
    … …  
    1430    // Exponential Term 
    1531    ////////////////////////double d(x[i]*x[i]*rg*rg); 
    16     double exp_term = square(q*gyration_radius); 
     32    double exp_term = q*q*gyration_radius*gyration_radius; 
    1733    exp_term = exp(-1.0*(exp_term/3.0)); 
    1834 
    … …  
    2137    return result; 
    2238} 
     39double form_volume(void){ 
     40    // Unused, so free to return garbage. 
     41    return NAN; 
     42} 
     43 
     44double Iq(double q, 
     45          double guinier_scale, 
     46          double lorentzian_scale, 
     47          double gyration_radius, 
     48          double fractal_exp, 
     49          double cor_length) 
     50{ 
     51    return _gel_fit_kernel(q, 
     52                          guinier_scale, 
     53                          lorentzian_scale, 
     54                          gyration_radius, 
     55                          fractal_exp, 
     56                          cor_length); 
     57} 
     58 
     59// Iqxy is never called since no orientation or magnetic parameters. 
     60double Iqxy(double qx, double qy, 
     61          double guinier_scale, 
     62          double lorentzian_scale, 
     63          double gyration_radius, 
     64          double fractal_exp, 
     65          double cor_length) 
     66{ 
     67    double q = sqrt(qx*qx + qy*qy); 
     68    return _gel_fit_kernel(q, 
     69                          guinier_scale, 
     70                          lorentzian_scale, 
     71                          gyration_radius, 
     72                          fractal_exp, 
     73                          cor_length); 
     74} 
     75 

  • sasmodels/models/hardsphere.py

    rd2bb604 rec45c4f  
    149149   """ 
    150150 
     151Iqxy = """ 
     152    // never called since no orientation or magnetic parameters. 
     153    return Iq(sqrt(qx*qx+qy*qy), IQ_PARAMETERS); 
     154    """ 
     155 
    151156# ER defaults to 0.0 
    152157# VR defaults to 1.0 

  • sasmodels/models/hayter_msa.py

    rd2bb604 rec45c4f  
    8787    return 1.0; 
    8888    """ 
     89Iqxy = """ 
     90    // never called since no orientation or magnetic parameters. 
     91    return Iq(sqrt(qx*qx+qy*qy), IQ_PARAMETERS); 
     92    """ 
    8993# ER defaults to 0.0 
    9094# VR defaults to 1.0 

  • sasmodels/models/lamellar.py

    rd2bb604 rec45c4f  
    8282    """ 
    8383 
     84Iqxy = """ 
     85    return Iq(sqrt(qx*qx+qy*qy), IQ_PARAMETERS); 
     86    """ 
     87 
    8488# ER defaults to 0.0 
    8589# VR defaults to 1.0 

  • sasmodels/models/lamellar_hg.py

    rd2bb604 rec45c4f  
    101101    """ 
    102102 
     103Iqxy = """ 
     104    return Iq(sqrt(qx*qx+qy*qy), IQ_PARAMETERS); 
     105    """ 
     106 
    103107# ER defaults to 0.0 
    104108# VR defaults to 1.0 

  • sasmodels/models/lamellar_hg_stack_caille.py

    rd2bb604 rec45c4f  
    120120    """ 
    121121 
     122Iqxy = """ 
     123    return Iq(sqrt(qx*qx+qy*qy), IQ_PARAMETERS); 
     124    """ 
     125 
    122126# ER defaults to 0.0 
    123127# VR defaults to 1.0 

  • sasmodels/models/lamellar_stack_caille.py

    rd2bb604 rec45c4f  
    104104    """ 
    105105 
     106Iqxy = """ 
     107    return Iq(sqrt(qx*qx+qy*qy), IQ_PARAMETERS); 
     108    """ 
     109 
    106110# ER defaults to 0.0 
    107111# VR defaults to 1.0 

  • sasmodels/models/lamellar_stack_paracrystal.py

    rd2bb604 rec45c4f  
    132132    """ 
    133133 
     134Iqxy = """ 
     135    return Iq(sqrt(qx*qx+qy*qy), IQ_PARAMETERS); 
     136    """ 
     137 
    134138# ER defaults to 0.0 
    135139# VR defaults to 1.0 

  • sasmodels/models/lib/sas_JN.c

    re6408d0 re6408d0  
    4848*/ 
    4949 
    50 double sas_JN( int n, double x ); 
    51  
    52 double sas_JN( int n, double x ) { 
     50static double 
     51sas_JN( int n, double x ) { 
    5352 
    5453    const double MACHEP = 1.11022302462515654042E-16; 

  • sasmodels/models/lib/sph_j1c.c

    rba32cdd re6f1410  
    77* using double precision that are the source. 
    88*/ 
    9 double sph_j1c(double q); 
    109 
    1110// The choice of the number of terms in the series and the cutoff value for 
    … …  
    4443#endif 
    4544 
     45double sph_j1c(double q); 
    4646double sph_j1c(double q) 
    4747{ 

  • sasmodels/models/lib/sphere_form.c

    rba32cdd rad90df9  
    1 double sphere_volume(double radius); 
    2 double sphere_form(double q, double radius, double sld, double solvent_sld); 
    3  
    4 double sphere_volume(double radius) 
     1inline double 
     2sphere_volume(double radius) 
    53{ 
    64    return M_4PI_3*cube(radius); 
    75} 
    86 
    9 double sphere_form(double q, double radius, double sld, double solvent_sld) 
     7inline double 
     8sphere_form(double q, double radius, double sld, double solvent_sld) 
    109{ 
    1110    const double fq = sphere_volume(radius) * sph_j1c(q*radius); 

  • sasmodels/models/lib/wrc_cyl.c

    rba32cdd re7678b2  
    22    Functions for WRC implementation of flexible cylinders 
    33*/ 
    4 double Sk_WR(double q, double L, double b); 
    5  
    6  
    74static double 
    85AlphaSquare(double x) 
    … …  
    366363} 
    367364 
     365double Sk_WR(double q, double L, double b); 
    368366double Sk_WR(double q, double L, double b) 
    369367{ 

  • sasmodels/models/onion.c

    rce896fd rfdb1487  
    44    double thickness, double A) 
    55{ 
    6   const double vol = M_4PI_3 * cube(r); 
     6  const double vol = 4.0/3.0 * M_PI * r * r * r; 
    77  const double qr = q * r; 
    88  const double alpha = A * r/thickness; 
    … …  
    1919    double sinqr, cosqr; 
    2020    SINCOS(qr, sinqr, cosqr); 
    21     fun = -3.0*( 
     21    fun = -3.0( 
    2222            ((alphasq - qrsq)*sinqr/qr - 2.0*alpha*cosqr) / sumsq 
    2323                - (alpha*sinqr/qr - cosqr) 
    … …  
    3232f_linear(double q, double r, double sld, double slope) 
    3333{ 
    34   const double vol = M_4PI_3 * cube(r); 
     34  const double vol = 4.0/3.0 * M_PI * r * r * r; 
    3535  const double qr = q * r; 
    3636  const double bes = sph_j1c(qr); 
    … …  
    5252{ 
    5353  const double bes = sph_j1c(q * r); 
    54   const double vol = M_4PI_3 * cube(r); 
     54  const double vol = 4.0/3.0 * M_PI * r * r * r; 
    5555  return sld * vol * bes; 
    5656} 
    … …  
    6464    r += thickness[i]; 
    6565  } 
    66   return M_4PI_3*cube(r); 
     66  return 4.0/3.0 * M_PI * r * r * r; 
    6767} 
    6868 
    6969static double 
    70 Iq(double q, double sld_core, double core_radius, double sld_solvent, 
    71     double n, double sld_in[], double sld_out[], double thickness[], 
     70Iq(double q, double core_sld, double core_radius, double solvent_sld, 
     71    double n, double in_sld[], double out_sld[], double thickness[], 
    7272    double A[]) 
    7373{ 
    7474  int i; 
    75   double r = core_radius; 
    76   double f = f_constant(q, r, sld_core); 
     75  r = core_radius; 
     76  f = f_constant(q, r, core_sld); 
    7777  for (i=0; i<n; i++){ 
    7878    const double r0 = r; 
    … …  
    9292    } 
    9393  } 
    94   f -= f_constant(q, r, sld_solvent); 
    95   const double f2 = f * f * 1.0e-4; 
     94  f -= f_constant(q, r, solvent_sld); 
     95  f2 = f * f * 1.0e-4; 
    9696 
    9797  return f2; 

  • sasmodels/models/onion.py

    rfa5fd8d r416609b  
    293293 
    294294#             ["name", "units", default, [lower, upper], "type","description"], 
    295 parameters = [["sld_core", "1e-6/Ang^2", 1.0, [-inf, inf], "", 
     295parameters = [["core_sld", "1e-6/Ang^2", 1.0, [-inf, inf], "", 
    296296               "Core scattering length density"], 
    297297              ["core_radius", "Ang", 200., [0, inf], "volume", 
    298298               "Radius of the core"], 
    299               ["sld_solvent", "1e-6/Ang^2", 6.4, [-inf, inf], "", 
     299              ["solvent_sld", "1e-6/Ang^2", 6.4, [-inf, inf], "", 
    300300               "Solvent scattering length density"], 
    301               ["n_shells", "", 1, [0, 10], "volume", 
     301              ["n", "", 1, [0, 10], "volume", 
    302302               "number of shells"], 
    303               ["sld_in[n_shells]", "1e-6/Ang^2", 1.7, [-inf, inf], "", 
     303              ["in_sld[n]", "1e-6/Ang^2", 1.7, [-inf, inf], "", 
    304304               "scattering length density at the inner radius of shell k"], 
    305               ["sld_out[n_shells]", "1e-6/Ang^2", 2.0, [-inf, inf], "", 
     305              ["out_sld[n]", "1e-6/Ang^2", 2.0, [-inf, inf], "", 
    306306               "scattering length density at the outer radius of shell k"], 
    307               ["thickness[n_shells]", "Ang", 40., [0, inf], "volume", 
     307              ["thickness[n]", "Ang", 40., [0, inf], "volume", 
    308308               "Thickness of shell k"], 
    309               ["A[n_shells]", "", 1.0, [-inf, inf], "", 
     309              ["A[n]", "", 1.0, [-inf, inf], "", 
    310310               "Decay rate of shell k"], 
    311311              ] 
    312312 
    313 source = ["lib/sph_j1c.c", "onion.c"] 
    314  
    315 #def Iq(q, *args, **kw): 
    316 #    return q 
    317  
    318 profile_axes = ['Radius (A)', 'SLD (1e-6/A^2)'] 
    319 def profile(core_sld, core_radius, solvent_sld, n_shells, 
    320             in_sld, out_sld, thickness, A): 
     313#source = ["lib/sph_j1c.c", "onion.c"] 
     314 
     315def Iq(q, *args, **kw): 
     316    return q 
     317 
     318def Iqxy(qx, *args, **kw): 
     319    return qx 
     320 
     321 
     322def shape(core_sld, core_radius, solvent_sld, n, in_sld, out_sld, thickness, A): 
    321323    """ 
    322     Returns shape profile with x=radius, y=SLD. 
     324    SLD profile 
    323325    """ 
    324326 
    325     total_radius = 1.25*(sum(thickness[:n_shells]) + core_radius + 1) 
     327    total_radius = 1.25*(sum(thickness[:n]) + core_radius + 1) 
    326328    dr = total_radius/400  # 400 points for a smooth plot 
    327329 
    … …  
    336338 
    337339    # add in the shells 
    338     for k in range(n_shells): 
     340    for k in range(n): 
    339341        # Left side of each shells 
    340342        r0 = r[-1] 
    … …  
    363365    beta.append(solvent_sld) 
    364366 
    365     return np.asarray(r), np.asarray(beta)*1e-6 
     367    return np.asarray(r), np.asarray(beta) 
    366368 
    367369def ER(core_radius, n, thickness): 
    … …  
    372374 
    373375demo = { 
    374     "sld_solvent": 2.2, 
    375     "sld_core": 1.0, 
     376    "solvent_sld": 2.2, 
     377    "core_sld": 1.0, 
    376378    "core_radius": 100, 
    377379    "n": 4, 
    378     "sld_in": [0.5, 1.5, 0.9, 2.0], 
    379     "sld_out": [nan, 0.9, 1.2, 1.6], 
     380    "in_sld": [0.5, 1.5, 0.9, 2.0], 
     381    "out_sld": [nan, 0.9, 1.2, 1.6], 
    380382    "thickness": [50, 75, 150, 75], 
    381383    "A": [0, -1, 1e-4, 1], 

  • sasmodels/models/rpa.c

    rd2bb604 r13ed84c  
    11double Iq(double q, double case_num, 
    2     double N[], double Phi[], double v[], double L[], double b[], 
     2    double Na, double Phia, double va, double a_sld, double ba, 
     3    double Nb, double Phib, double vb, double b_sld, double bb, 
     4    double Nc, double Phic, double vc, double c_sld, double bc, 
     5    double Nd, double Phid, double vd, double d_sld, double bd, 
    36    double Kab, double Kac, double Kad, 
    47    double Kbc, double Kbd, double Kcd 
    58    ); 
    69 
     10double Iqxy(double qx, double qy, double case_num, 
     11    double Na, double Phia, double va, double a_sld, double ba, 
     12    double Nb, double Phib, double vb, double b_sld, double bb, 
     13    double Nc, double Phic, double vc, double c_sld, double bc, 
     14    double Nd, double Phid, double vd, double d_sld, double bd, 
     15    double Kab, double Kac, double Kad, 
     16    double Kbc, double Kbd, double Kcd 
     17    ); 
     18 
     19double form_volume(void); 
     20 
     21double form_volume(void) 
     22{ 
     23    return 1.0; 
     24} 
     25 
    726double Iq(double q, double case_num, 
    8     double N[], double Phi[], double v[], double L[], double b[], 
     27    double Na, double Phia, double va, double La, double ba, 
     28    double Nb, double Phib, double vb, double Lb, double bb, 
     29    double Nc, double Phic, double vc, double Lc, double bc, 
     30    double Nd, double Phid, double vd, double Ld, double bd, 
    931    double Kab, double Kac, double Kad, 
    1032    double Kbc, double Kbd, double Kcd 
    … …  
    1436#if 0  // Sasview defaults 
    1537  if (icase <= 1) { 
    16     N[0]=N[1]=1000.0; 
    17     Phi[0]=Phi[1]=0.0000001; 
     38    Na=Nb=1000.0; 
     39    Phia=Phib=0.0000001; 
    1840    Kab=Kac=Kad=Kbc=Kbd=-0.0004; 
    19     L[0]=L[1]=1e-12; 
    20     v[0]=v[1]=100.0; 
    21     b[0]=b[1]=5.0; 
     41    La=Lb=1e-12; 
     42    va=vb=100.0; 
     43    ba=bb=5.0; 
    2244  } else if (icase <= 4) { 
    23     Phi[0]=0.0000001; 
     45    Phia=0.0000001; 
    2446    Kab=Kac=Kad=-0.0004; 
    25     L[0]=1e-12; 
    26     v[0]=100.0; 
    27     b[0]=5.0; 
     47    La=1e-12; 
     48    va=100.0; 
     49    ba=5.0; 
    2850  } 
    2951#else 
    3052  if (icase <= 1) { 
    31     N[0]=N[1]=0.0; 
    32     Phi[0]=Phi[1]=0.0; 
     53    Na=Nb=0.0; 
     54    Phia=Phib=0.0; 
    3355    Kab=Kac=Kad=Kbc=Kbd=0.0; 
    34     L[0]=L[1]=L[3]; 
    35     v[0]=v[1]=v[3]; 
    36     b[0]=b[1]=0.0; 
     56    La=Lb=Ld; 
     57    va=vb=vd; 
     58    ba=bb=0.0; 
    3759  } else if (icase <= 4) { 
    38     N[0] = 0.0; 
    39     Phi[0]=0.0; 
     60    Na = 0.0; 
     61    Phia=0.0; 
    4062    Kab=Kac=Kad=0.0; 
    41     L[0]=L[3]; 
    42     v[0]=v[3]; 
    43     b[0]=0.0; 
     63    La=Ld; 
     64    va=vd; 
     65    ba=0.0; 
    4466  } 
    4567#endif 
    4668 
    47   const double Xa = q*q*b[0]*b[0]*N[0]/6.0; 
    48   const double Xb = q*q*b[1]*b[1]*N[1]/6.0; 
    49   const double Xc = q*q*b[2]*b[2]*N[2]/6.0; 
    50   const double Xd = q*q*b[3]*b[3]*N[3]/6.0; 
     69  const double Xa = q*q*ba*ba*Na/6.0; 
     70  const double Xb = q*q*bb*bb*Nb/6.0; 
     71  const double Xc = q*q*bc*bc*Nc/6.0; 
     72  const double Xd = q*q*bd*bd*Nd/6.0; 
    5173 
    5274  // limit as Xa goes to 0 is 1 
    … …  
    7698#if 0 
    7799  const double S0aa = icase<5 
    78                       ? 1.0 : N[0]*Phi[0]*v[0]*Paa; 
     100                      ? 1.0 : Na*Phia*va*Paa; 
    79101  const double S0bb = icase<2 
    80                       ? 1.0 : N[1]*Phi[1]*v[1]*Pbb; 
    81   const double S0cc = N[2]*Phi[2]*v[2]*Pcc; 
    82   const double S0dd = N[3]*Phi[3]*v[3]*Pdd; 
     102                      ? 1.0 : Nb*Phib*vb*Pbb; 
     103  const double S0cc = Nc*Phic*vc*Pcc; 
     104  const double S0dd = Nd*Phid*vd*Pdd; 
    83105  const double S0ab = icase<8 
    84                       ? 0.0 : sqrt(N[0]*v[0]*Phi[0]*N[1]*v[1]*Phi[1])*Pa*Pb; 
     106                      ? 0.0 : sqrt(Na*va*Phia*Nb*vb*Phib)*Pa*Pb; 
    85107  const double S0ac = icase<9 
    86                       ? 0.0 : sqrt(N[0]*v[0]*Phi[0]*N[2]*v[2]*Phi[2])*Pa*Pc*exp(-Xb); 
     108                      ? 0.0 : sqrt(Na*va*Phia*Nc*vc*Phic)*Pa*Pc*exp(-Xb); 
    87109  const double S0ad = icase<9 
    88                       ? 0.0 : sqrt(N[0]*v[0]*Phi[0]*N[3]*v[3]*Phi[3])*Pa*Pd*exp(-Xb-Xc); 
     110                      ? 0.0 : sqrt(Na*va*Phia*Nd*vd*Phid)*Pa*Pd*exp(-Xb-Xc); 
    89111  const double S0bc = (icase!=4 && icase!=7 && icase!= 9) 
    90                       ? 0.0 : sqrt(N[1]*v[1]*Phi[1]*N[2]*v[2]*Phi[2])*Pb*Pc; 
     112                      ? 0.0 : sqrt(Nb*vb*Phib*Nc*vc*Phic)*Pb*Pc; 
    91113  const double S0bd = (icase!=4 && icase!=7 && icase!= 9) 
    92                       ? 0.0 : sqrt(N[1]*v[1]*Phi[1]*N[3]*v[3]*Phi[3])*Pb*Pd*exp(-Xc); 
     114                      ? 0.0 : sqrt(Nb*vb*Phib*Nd*vd*Phid)*Pb*Pd*exp(-Xc); 
    93115  const double S0cd = (icase==0 || icase==2 || icase==5) 
    94                       ? 0.0 : sqrt(N[2]*v[2]*Phi[2]*N[3]*v[3]*Phi[3])*Pc*Pd; 
     116                      ? 0.0 : sqrt(Nc*vc*Phic*Nd*vd*Phid)*Pc*Pd; 
    95117#else  // sasview equivalent 
    96 //printf("Xc=%g, S0cc=%g*%g*%g*%g\n",Xc,N[2],Phi[2],v[2],Pcc); 
    97   double S0aa = N[0]*Phi[0]*v[0]*Paa; 
    98   double S0bb = N[1]*Phi[1]*v[1]*Pbb; 
    99   double S0cc = N[2]*Phi[2]*v[2]*Pcc; 
    100   double S0dd = N[3]*Phi[3]*v[3]*Pdd; 
    101   double S0ab = sqrt(N[0]*v[0]*Phi[0]*N[1]*v[1]*Phi[1])*Pa*Pb; 
    102   double S0ac = sqrt(N[0]*v[0]*Phi[0]*N[2]*v[2]*Phi[2])*Pa*Pc*exp(-Xb); 
    103   double S0ad = sqrt(N[0]*v[0]*Phi[0]*N[3]*v[3]*Phi[3])*Pa*Pd*exp(-Xb-Xc); 
    104   double S0bc = sqrt(N[1]*v[1]*Phi[1]*N[2]*v[2]*Phi[2])*Pb*Pc; 
    105   double S0bd = sqrt(N[1]*v[1]*Phi[1]*N[3]*v[3]*Phi[3])*Pb*Pd*exp(-Xc); 
    106   double S0cd = sqrt(N[2]*v[2]*Phi[2]*N[3]*v[3]*Phi[3])*Pc*Pd; 
     118//printf("Xc=%g, S0cc=%g*%g*%g*%g\n",Xc,Nc,Phic,vc,Pcc); 
     119  double S0aa = Na*Phia*va*Paa; 
     120  double S0bb = Nb*Phib*vb*Pbb; 
     121  double S0cc = Nc*Phic*vc*Pcc; 
     122  double S0dd = Nd*Phid*vd*Pdd; 
     123  double S0ab = sqrt(Na*va*Phia*Nb*vb*Phib)*Pa*Pb; 
     124  double S0ac = sqrt(Na*va*Phia*Nc*vc*Phic)*Pa*Pc*exp(-Xb); 
     125  double S0ad = sqrt(Na*va*Phia*Nd*vd*Phid)*Pa*Pd*exp(-Xb-Xc); 
     126  double S0bc = sqrt(Nb*vb*Phib*Nc*vc*Phic)*Pb*Pc; 
     127  double S0bd = sqrt(Nb*vb*Phib*Nd*vd*Phid)*Pb*Pd*exp(-Xc); 
     128  double S0cd = sqrt(Nc*vc*Phic*Nd*vd*Phid)*Pc*Pd; 
    107129switch(icase){ 
    108130  case 0: 
    … …  
    289311  // Note: 1e-13 to convert from fm to cm for scattering length 
    290312  const double sqrt_Nav=sqrt(6.022045e+23) * 1.0e-13; 
    291   const double Lad = icase<5 ? 0.0 : (L[0]/v[0] - L[3]/v[3])*sqrt_Nav; 
    292   const double Lbd = icase<2 ? 0.0 : (L[1]/v[1] - L[3]/v[3])*sqrt_Nav; 
    293   const double Lcd = (L[2]/v[2] - L[3]/v[3])*sqrt_Nav; 
     313  const double Lad = icase<5 ? 0.0 : (La/va - Ld/vd)*sqrt_Nav; 
     314  const double Lbd = icase<2 ? 0.0 : (Lb/vb - Ld/vd)*sqrt_Nav; 
     315  const double Lcd = (Lc/vc - Ld/vd)*sqrt_Nav; 
    294316 
    295317  const double result=Lad*Lad*S11 + Lbd*Lbd*S22 + Lcd*Lcd*S33 
    … …  
    299321 
    300322} 
     323 
     324double Iqxy(double qx, double qy, 
     325    double case_num, 
     326    double Na, double Phia, double va, double a_sld, double ba, 
     327    double Nb, double Phib, double vb, double b_sld, double bb, 
     328    double Nc, double Phic, double vc, double c_sld, double bc, 
     329    double Nd, double Phid, double vd, double d_sld, double bd, 
     330    double Kab, double Kac, double Kad, 
     331    double Kbc, double Kbd, double Kcd 
     332    ) 
     333{ 
     334    double q = sqrt(qx*qx + qy*qy); 
     335    return Iq(q, 
     336        case_num, 
     337        Na, Phia, va, a_sld, ba, 
     338        Nb, Phib, vb, b_sld, bb, 
     339        Nc, Phic, vc, c_sld, bc, 
     340        Nd, Phid, vd, d_sld, bd, 
     341        Kab, Kac, Kad, 
     342        Kbc, Kbd, Kcd); 
     343} 

  • sasmodels/models/rpa.py

    ra5b8477 rec45c4f  
    8686#   ["name", "units", default, [lower, upper], "type","description"], 
    8787parameters = [ 
    88     ["case_num", "", 1, [CASES], "", "Component organization"], 
     88    ["case_num", CASES, 0, [0, 10], "", "Component organization"], 
    8989 
    90     ["N[4]", "", 1000.0, [1, inf], "", "Degree of polymerization"], 
    91     ["Phi[4]", "", 0.25, [0, 1], "", "volume fraction"], 
    92     ["v[4]", "mL/mol", 100.0, [0, inf], "", "specific volume"], 
    93     ["L[4]", "fm", 10.0, [-inf, inf], "", "scattering length"], 
    94     ["b[4]", "Ang", 5.0, [0, inf], "", "segment length"], 
     90    ["Na", "", 1000.0, [1, inf], "", "Degree of polymerization"], 
     91    ["Phia", "", 0.25, [0, 1], "", "volume fraction"], 
     92    ["va", "mL/mol", 100.0, [0, inf], "", "specific volume"], 
     93    ["La", "fm", 10.0, [-inf, inf], "", "scattering length"], 
     94    ["ba", "Ang", 5.0, [0, inf], "", "segment length"], 
     95 
     96    ["Nb", "", 1000.0, [1, inf], "", "Degree of polymerization"], 
     97    ["Phib", "", 0.25, [0, 1], "", "volume fraction"], 
     98    ["vb", "mL/mol", 100.0, [0, inf], "", "specific volume"], 
     99    ["Lb", "fm", 10.0, [-inf, inf], "", "scattering length"], 
     100    ["bb", "Ang", 5.0, [0, inf], "", "segment length"], 
     101 
     102    ["Nc", "", 1000.0, [1, inf], "", "Degree of polymerization"], 
     103    ["Phic", "", 0.25, [0, 1], "", "volume fraction"], 
     104    ["vc", "mL/mol", 100.0, [0, inf], "", "specific volume"], 
     105    ["Lc", "fm", 10.0, [-inf, inf], "", "scattering length"], 
     106    ["bc", "Ang", 5.0, [0, inf], "", "segment length"], 
     107 
     108    ["Nd", "", 1000.0, [1, inf], "", "Degree of polymerization"], 
     109    ["Phid", "", 0.25, [0, 1], "", "volume fraction"], 
     110    ["vd", "mL/mol", 100.0, [0, inf], "", "specific volume"], 
     111    ["Ld", "fm", 10.0, [-inf, inf], "", "scattering length"], 
     112    ["bd", "Ang", 5.0, [0, inf], "", "segment length"], 
    95113 
    96114    ["Kab", "", -0.0004, [-inf, inf], "", "Interaction parameter"], 

  • sasmodels/models/spherical_sld.py

    rd2bb604 rec45c4f  
    170170# pylint: disable=bad-whitespace, line-too-long 
    171171#            ["name", "units", default, [lower, upper], "type", "description"], 
    172 parameters = [["n",                "",               1,      [0, 9],         "", "number of shells"], 
     172parameters = [["n_shells",         "",               1,      [0, 9],         "", "number of shells"], 
    173173              ["radius_core",      "Ang",            50.0,   [0, inf],       "", "intern layer thickness"], 
    174174              ["sld_core",         "1e-6/Ang^2",     2.07,   [-inf, inf],    "", "sld function flat"], 
    … …  
    192192 
    193193demo = dict( 
    194     n=4, 
     194    n_shells=4, 
    195195    scale=1.0, 
    196196    solvent_sld=1.0, 

  • sasmodels/models/squarewell.py

    rd2bb604 rec45c4f  
    123123""" 
    124124 
     125Iqxy = """ 
     126    return Iq(sqrt(qx*qx+qy*qy), IQ_PARAMETERS); 
     127    """ 
     128 
    125129# ER defaults to 0.0 
    126130# VR defaults to 1.0 

  • sasmodels/models/stickyhardsphere.py

    rd2bb604 rec45c4f  
    171171""" 
    172172 
     173Iqxy = """ 
     174    return Iq(sqrt(qx*qx+qy*qy), IQ_PARAMETERS); 
     175    """ 
     176 
    173177# ER defaults to 0.0 
    174178# VR defaults to 1.0 

  • sasmodels/product.py

    r7ae2b7f rf247314  
    1111*ProductModel(P, S)*. 
    1212""" 
    13 import numpy as np  # type: ignore 
     13import numpy as np 
    1414 
    15 from .details import dispersion_mesh 
    16 from .modelinfo import suffix_parameter, ParameterTable, ModelInfo 
    17 from .kernel import KernelModel, Kernel 
     15from .core import call_ER_VR 
     16from .generate import process_parameters 
    1817 
    19 try: 
    20     from typing import Tuple 
    21     from .modelinfo import ParameterSet 
    22     from .details import CallDetails 
    23 except ImportError: 
    24     pass 
    25  
    26 # TODO: make estimates available to constraints 
    27 #ESTIMATED_PARAMETERS = [ 
    28 #    ["est_effect_radius", "A", 0.0, [0, np.inf], "", "Estimated effective radius"], 
    29 #    ["est_volume_ratio", "", 1.0, [0, np.inf], "", "Estimated volume ratio"], 
    30 #] 
     18SCALE=0 
     19BACKGROUND=1 
     20RADIUS_EFFECTIVE=2 
     21VOLFRACTION=3 
    3122 
    3223# TODO: core_shell_sphere model has suppressed the volume ratio calculation 
    3324# revert it after making VR and ER available at run time as constraints. 
    3425def make_product_info(p_info, s_info): 
    35     # type: (ModelInfo, ModelInfo) -> ModelInfo 
    3626    """ 
    3727    Create info block for product model. 
    3828    """ 
    39     p_id, p_name, p_pars = p_info.id, p_info.name, p_info.parameters 
    40     s_id, s_name, s_pars = s_info.id, s_info.name, s_info.parameters 
    41     p_set = set(p.id for p in p_pars.call_parameters) 
    42     s_set = set(p.id for p in s_pars.call_parameters) 
    43  
    44     if p_set & s_set: 
    45         # there is some overlap between the parameter names; tag the 
    46         # overlapping S parameters with name_S 
    47         s_list = [(suffix_parameter(par, "_S") if par.id in p_set else par) 
    48                   for par in s_pars.kernel_parameters] 
    49         combined_pars = p_pars.kernel_parameters + s_list 
    50     else: 
    51         combined_pars = p_pars.kernel_parameters + s_pars.kernel_parameters 
    52     parameters = ParameterTable(combined_pars) 
    53  
    54     model_info = ModelInfo() 
    55     model_info.id = '*'.join((p_id, s_id)) 
    56     model_info.name = ' X '.join((p_name, s_name)) 
    57     model_info.filename = None 
    58     model_info.title = 'Product of %s and %s'%(p_name, s_name) 
    59     model_info.description = model_info.title 
    60     model_info.docs = model_info.title 
    61     model_info.category = "custom" 
    62     model_info.parameters = parameters 
    63     #model_info.single = p_info.single and s_info.single 
    64     model_info.structure_factor = False 
    65     model_info.variant_info = None 
    66     #model_info.tests = [] 
    67     #model_info.source = [] 
     29    p_id, p_name, p_pars = p_info['id'], p_info['name'], p_info['parameters'] 
     30    s_id, s_name, s_pars = s_info['id'], s_info['name'], s_info['parameters'] 
     31    # We require models to start with scale and background 
     32    assert s_pars[SCALE].name == 'scale' 
     33    assert s_pars[BACKGROUND].name == 'background' 
     34    # We require structure factors to start with effect radius and volfraction 
     35    assert s_pars[RADIUS_EFFECTIVE].name == 'radius_effective' 
     36    assert s_pars[VOLFRACTION].name == 'volfraction' 
     37    # Combine the parameter sets.  We are skipping the first three 
     38    # parameters of S since scale, background are defined in P and 
     39    # effect_radius is set from P.ER(). 
     40    pars = p_pars + s_pars[3:] 
     41    # check for duplicates; can't use assertion since they may never be checked 
     42    if len(set(p.name for p in pars)) != len(pars): 
     43        raise ValueError("Duplicate parameters in %s and %s"%(p_id)) 
     44    model_info = {} 
     45    model_info['id'] = '*'.join((p_id, s_id)) 
     46    model_info['name'] = ' X '.join((p_name, s_name)) 
     47    model_info['filename'] = None 
     48    model_info['title'] = 'Product of %s and structure factor %s'%(p_name, s_name) 
     49    model_info['description'] = model_info['title'] 
     50    model_info['docs'] = model_info['title'] 
     51    model_info['category'] = "custom" 
     52    model_info['parameters'] = pars 
     53    #model_info['single'] = p_info['single'] and s_info['single'] 
     54    model_info['structure_factor'] = False 
     55    model_info['variant_info'] = None 
     56    #model_info['tests'] = [] 
     57    #model_info['source'] = [] 
    6858    # Iq, Iqxy, form_volume, ER, VR and sesans 
    69     model_info.composition = ('product', [p_info, s_info]) 
     59    model_info['composition'] = ('product', [p_info, s_info]) 
     60    process_parameters(model_info) 
    7061    return model_info 
    7162 
    72 class ProductModel(KernelModel): 
     63class ProductModel(object): 
    7364    def __init__(self, model_info, P, S): 
    74         # type: (ModelInfo, KernelModel, KernelModel) -> None 
    7565        self.info = model_info 
    7666        self.P = P 
    … …  
    7868 
    7969    def __call__(self, q_vectors): 
    80         # type: (List[np.ndarray]) -> Kernel 
    8170        # Note: may be sending the q_vectors to the GPU twice even though they 
    8271        # are only needed once.  It would mess up modularity quite a bit to 
    … …  
    8574        # in opencl; or both in opencl, but one in single precision and the 
    8675        # other in double precision). 
    87         p_kernel = self.P.make_kernel(q_vectors) 
    88         s_kernel = self.S.make_kernel(q_vectors) 
     76        p_kernel = self.P(q_vectors) 
     77        s_kernel = self.S(q_vectors) 
    8978        return ProductKernel(self.info, p_kernel, s_kernel) 
    9079 
    9180    def release(self): 
    92         # type: (None) -> None 
    9381        """ 
    9482        Free resources associated with the model. 
    … …  
    9886 
    9987 
    100 class ProductKernel(Kernel): 
     88class ProductKernel(object): 
    10189    def __init__(self, model_info, p_kernel, s_kernel): 
    102         # type: (ModelInfo, Kernel, Kernel) -> None 
     90        dim = '2d' if p_kernel.q_input.is_2d else '1d' 
     91 
     92        # Need to know if we want 2D and magnetic parameters when constructing 
     93        # a parameter map. 
     94        par_map = {} 
     95        p_info = p_kernel.info['partype'] 
     96        s_info = s_kernel.info['partype'] 
     97        vol_pars = set(p_info['volume']) 
     98        if dim == '2d': 
     99            num_p_fixed = len(p_info['fixed-2d']) 
     100            num_p_pd = len(p_info['pd-2d']) 
     101            num_s_fixed = len(s_info['fixed-2d']) 
     102            num_s_pd = len(s_info['pd-2d']) - 1 # exclude effect_radius 
     103            # volume parameters are amongst the pd pars for P, not S 
     104            vol_par_idx = [k for k,v in enumerate(p_info['pd-2d']) 
     105                           if v in vol_pars] 
     106        else: 
     107            num_p_fixed = len(p_info['fixed-1d']) 
     108            num_p_pd = len(p_info['pd-1d']) 
     109            num_s_fixed = len(s_info['fixed-1d']) 
     110            num_s_pd = len(s_info['pd-1d']) - 1  # exclude effect_radius 
     111            # volume parameters are amongst the pd pars for P, not S 
     112            vol_par_idx = [k for k,v in enumerate(p_info['pd-1d']) 
     113                           if v in vol_pars] 
     114 
     115        start = 0 
     116        par_map['p_fixed'] = np.arange(start, start+num_p_fixed) 
     117        # User doesn't set scale, background or effect_radius for S in P*S, 
     118        # so borrow values from end of p_fixed.  This makes volfraction the 
     119        # first S parameter. 
     120        start += num_p_fixed 
     121        par_map['s_fixed'] = np.hstack(([start,start], 
     122                                        np.arange(start, start+num_s_fixed-2))) 
     123        par_map['volfraction'] = num_p_fixed 
     124        start += num_s_fixed-2 
     125        # vol pars offset from the start of pd pars 
     126        par_map['vol_pars'] = [start+k for k in vol_par_idx] 
     127        par_map['p_pd'] = np.arange(start, start+num_p_pd) 
     128        start += num_p_pd-1 
     129        par_map['s_pd'] = np.hstack((start, 
     130                                     np.arange(start, start+num_s_pd-1))) 
     131 
     132        self.fixed_pars = model_info['partype']['fixed-' + dim] 
     133        self.pd_pars = model_info['partype']['pd-' + dim] 
    103134        self.info = model_info 
    104135        self.p_kernel = p_kernel 
    105136        self.s_kernel = s_kernel 
     137        self.par_map = par_map 
    106138 
    107     def __call__(self, details, weights, values, cutoff): 
    108         # type: (CallDetails, np.ndarray, np.ndarray, float) -> np.ndarray 
     139    def __call__(self, fixed_pars, pd_pars, cutoff=1e-5): 
     140        pars = fixed_pars + pd_pars 
     141        scale = pars[SCALE] 
     142        background = pars[BACKGROUND] 
     143        s_volfraction = pars[self.par_map['volfraction']] 
     144        p_fixed = [pars[k] for k in self.par_map['p_fixed']] 
     145        s_fixed = [pars[k] for k in self.par_map['s_fixed']] 
     146        p_pd = [pars[k] for k in self.par_map['p_pd']] 
     147        s_pd = [pars[k] for k in self.par_map['s_pd']] 
     148        vol_pars = [pars[k] for k in self.par_map['vol_pars']] 
     149 
    109150        effect_radius, vol_ratio = call_ER_VR(self.p_kernel.info, vol_pars) 
    110151 
    111         p_details, s_details = details.parts 
    112         p_res = self.p_kernel(p_details, weights, values, cutoff) 
    113         s_res = self.s_kernel(s_details, weights, values, cutoff) 
     152        p_fixed[SCALE] = s_volfraction 
     153        p_fixed[BACKGROUND] = 0.0 
     154        s_fixed[SCALE] = scale 
     155        s_fixed[BACKGROUND] = 0.0 
     156        s_fixed[2] = s_volfraction/vol_ratio 
     157        s_pd[0] = [effect_radius], [1.0] 
     158 
     159        p_res = self.p_kernel(p_fixed, p_pd, cutoff) 
     160        s_res = self.s_kernel(s_fixed, s_pd, cutoff) 
    114161        #print s_fixed, s_pd, p_fixed, p_pd 
    115162 
    116         return values[0]*(p_res*s_res) + values[1] 
     163        return p_res*s_res + background 
    117164 
    118165    def release(self): 
    119         # type: () -> None 
    120166        self.p_kernel.release() 
    121         self.s_kernel.release() 
     167        self.q_kernel.release() 
    122168 
    123 def call_ER_VR(model_info, pars): 
    124     """ 
    125     Return effect radius and volume ratio for the model. 
    126     """ 
    127     if model_info.ER is None and model_info.VR is None: 
    128         return 1.0, 1.0 
    129  
    130     value, weight = _vol_pars(model_info, pars) 
    131  
    132     if model_info.ER is not None: 
    133         individual_radii = model_info.ER(*value) 
    134         effect_radius = np.sum(weight*individual_radii) / np.sum(weight) 
    135     else: 
    136         effect_radius = 1.0 
    137  
    138     if model_info.VR is not None: 
    139         whole, part = model_info.VR(*value) 
    140         volume_ratio = np.sum(weight*part)/np.sum(weight*whole) 
    141     else: 
    142         volume_ratio = 1.0 
    143  
    144     return effect_radius, volume_ratio 
    145  
    146 def _vol_pars(model_info, pars): 
    147     # type: (ModelInfo, ParameterSet) -> Tuple[np.ndarray, np.ndarray] 
    148     vol_pars = [get_weights(p, pars) 
    149                 for p in model_info.parameters.call_parameters 
    150                 if p.type == 'volume'] 
    151     value, weight = dispersion_mesh(model_info, vol_pars) 
    152     return value, weight 
    153  

  • sasmodels/resolution.py

    r7ae2b7f r4d76711  
    66from __future__ import division 
    77 
    8 from scipy.special import erf  # type: ignore 
    9 from numpy import sqrt, log, log10, exp, pi  # type: ignore 
    10 import numpy as np  # type: ignore 
     8from scipy.special import erf 
     9from numpy import sqrt, log, log10 
     10import numpy as np 
    1111 
    1212__all__ = ["Resolution", "Perfect1D", "Pinhole1D", "Slit1D", 
    … …  
    3535    smeared theory I(q). 
    3636    """ 
    37     q = None  # type: np.ndarray 
    38     q_calc = None  # type: np.ndarray 
     37    q = None 
     38    q_calc = None 
    3939    def apply(self, theory): 
    4040        """ 
    … …  
    476476    *pars* are the parameter values to use when evaluating. 
    477477    """ 
    478     from sasmodels import direct_model 
     478    from sasmodels import core 
    479479    kernel = form.make_kernel([q]) 
    480     theory = direct_model.call_kernel(kernel, pars) 
     480    theory = core.call_kernel(kernel, pars) 
    481481    kernel.release() 
    482482    return theory 
    … …  
    489489    *q0* is the center, *dq* is the width and *q* are the points to evaluate. 
    490490    """ 
     491    from numpy import exp, pi 
    491492    return exp(-0.5*((q-q0)/dq)**2)/(sqrt(2*pi)*dq) 
    492493 
    … …  
    499500    that make it slow to evaluate but give it good accuracy. 
    500501    """ 
    501     from scipy.integrate import romberg  # type: ignore 
    502  
    503     par_set = set([p.name for p in form.info.parameters.call_parameters]) 
    504     if any(k not in par_set for k in pars.keys()): 
    505         extra = set(pars.keys()) - par_set 
    506         raise ValueError("bad parameters: [%s] not in [%s]" 
    507                          % (", ".join(sorted(extra)), 
    508                             ", ".join(sorted(pars.keys())))) 
     502    from scipy.integrate import romberg 
     503 
     504    if any(k not in form.info['defaults'] for k in pars.keys()): 
     505        keys = set(form.info['defaults'].keys()) 
     506        extra = set(pars.keys()) - keys 
     507        raise ValueError("bad parameters: [%s] not in [%s]"% 
     508                         (", ".join(sorted(extra)), ", ".join(sorted(keys)))) 
    509509 
    510510    if np.isscalar(width): 
    … …  
    554554    that make it slow to evaluate but give it good accuracy. 
    555555    """ 
    556     from scipy.integrate import romberg  # type: ignore 
    557  
    558     par_set = set([p.name for p in form.info.parameters.call_parameters]) 
    559     if any(k not in par_set for k in pars.keys()): 
    560         extra = set(pars.keys()) - par_set 
    561         raise ValueError("bad parameters: [%s] not in [%s]" 
    562                          % (", ".join(sorted(extra)), 
    563                             ", ".join(sorted(pars.keys())))) 
     556    from scipy.integrate import romberg 
     557 
     558    if any(k not in form.info['defaults'] for k in pars.keys()): 
     559        keys = set(form.info['defaults'].keys()) 
     560        extra = set(pars.keys()) - keys 
     561        raise ValueError("bad parameters: [%s] not in [%s]"% 
     562                         (", ".join(sorted(extra)), ", ".join(sorted(keys)))) 
    564563 
    565564    _fn = lambda q, q0, dq: eval_form(q, form, pars)*gaussian(q, q0, dq) 
    … …  
    693692 
    694693    def _eval_sphere(self, pars, resolution): 
    695         from sasmodels import direct_model 
     694        from sasmodels import core 
    696695        kernel = self.model.make_kernel([resolution.q_calc]) 
    697         theory = direct_model.call_kernel(kernel, pars) 
     696        theory = core.call_kernel(kernel, pars) 
    698697        result = resolution.apply(theory) 
    699698        kernel.release() 
    … …  
    751750        #tol = 0.001 
    752751        ## The default 3 sigma and no extra points gets 1% 
    753         q_calc = None  # type: np.ndarray 
    754         tol = 0.01 
     752        q_calc, tol = None, 0.01 
    755753        resolution = Pinhole1D(q, q_width, q_calc=q_calc) 
    756754        output = self._eval_sphere(pars, resolution) 
    757755        if 0: # debug plot 
    758             import matplotlib.pyplot as plt  # type: ignore 
     756            import matplotlib.pyplot as plt 
    759757            resolution = Perfect1D(q) 
    760758            source = self._eval_sphere(pars, resolution) 
    … …  
    10281026    """ 
    10291027    import sys 
    1030     import xmlrunner  # type: ignore 
     1028    import xmlrunner 
    10311029 
    10321030    suite = unittest.TestSuite() 
    … …  
    10451043    import sys 
    10461044    from sasmodels import core 
    1047     from sasmodels import direct_model 
    10481045    name = sys.argv[1] if len(sys.argv) > 1 else 'cylinder' 
    10491046 
    … …  
    10661063 
    10671064    kernel = model.make_kernel([resolution.q_calc]) 
    1068     theory = direct_model.call_kernel(kernel, pars) 
     1065    theory = core.call_kernel(kernel, pars) 
    10691066    Iq = resolution.apply(theory) 
    10701067 
    … …  
    10761073        Iq_romb = romberg_pinhole_1d(resolution.q, dq, model, pars) 
    10771074 
    1078     import matplotlib.pyplot as plt  # type: ignore 
     1075    import matplotlib.pyplot as plt 
    10791076    plt.loglog(resolution.q_calc, theory, label='unsmeared') 
    10801077    plt.loglog(resolution.q, Iq, label='smeared', hold=True) 
    … …  
    11111108    Run the resolution demos. 
    11121109    """ 
    1113     import matplotlib.pyplot as plt  # type: ignore 
     1110    import matplotlib.pyplot as plt 
    11141111    plt.subplot(121) 
    11151112    demo_pinhole_1d() 

  • sasmodels/resolution2d.py

    r7ae2b7f rd6f5da6  
    77from __future__ import division 
    88 
    9 import numpy as np  # type: ignore 
    10 from numpy import pi, cos, sin, sqrt  # type: ignore 
     9import numpy as np 
     10from numpy import pi, cos, sin, sqrt 
    1111 
    1212from . import resolution 

  • sasmodels/sasview_model.py

    r60f03de r81ec7c8  
    2121import logging 
    2222 
    23 import numpy as np  # type: ignore 
     23import numpy as np 
    2424 
    2525from . import core 
    2626from . import custom 
    2727from . import generate 
    28 from . import weights 
    29 from . import details 
    30 from . import modelinfo 
    3128 
    3229try: 
    33     from typing import Dict, Mapping, Any, Sequence, Tuple, NamedTuple, List, Optional, Union, Callable 
    34     from .modelinfo import ModelInfo, Parameter 
     30    from typing import Dict, Mapping, Any, Sequence, Tuple, NamedTuple, List, Optional 
    3531    from .kernel import KernelModel 
    3632    MultiplicityInfoType = NamedTuple( 
    … …  
    3834        [("number", int), ("control", str), ("choices", List[str]), 
    3935         ("x_axis_label", str)]) 
    40     SasviewModelType = Callable[[int], "SasviewModel"] 
    4136except ImportError: 
    4237    pass 
    4338 
    4439# TODO: separate x_axis_label from multiplicity info 
    45 # The profile x-axis label belongs with the profile generating function 
     40# The x-axis label belongs with the profile generating function 
    4641MultiplicityInfo = collections.namedtuple( 
    4742    'MultiplicityInfo', 
    … …  
    4944) 
    5045 
    51 # TODO: figure out how to say that the return type is a subclass 
    5246def load_standard_models(): 
    53     # type: () -> List[SasviewModelType] 
    5447    """ 
    5548    Load and return the list of predefined models. 
    … …  
    6255        try: 
    6356            models.append(_make_standard_model(name)) 
    64         except Exception: 
     57        except: 
    6558            logging.error(traceback.format_exc()) 
    6659    return models 
    … …  
    6861 
    6962def load_custom_model(path): 
    70     # type: (str) -> SasviewModelType 
    7163    """ 
    7264    Load a custom model given the model path. 
    7365    """ 
    7466    kernel_module = custom.load_custom_kernel_module(path) 
    75     model_info = modelinfo.make_model_info(kernel_module) 
     67    model_info = generate.make_model_info(kernel_module) 
    7668    return _make_model_from_info(model_info) 
    7769 
    7870 
    7971def _make_standard_model(name): 
    80     # type: (str) -> SasviewModelType 
    8172    """ 
    8273    Load the sasview model defined by *name*. 
    … …  
    8778    """ 
    8879    kernel_module = generate.load_kernel_module(name) 
    89     model_info = modelinfo.make_model_info(kernel_module) 
     80    model_info = generate.make_model_info(kernel_module) 
    9081    return _make_model_from_info(model_info) 
    9182 
    9283 
    9384def _make_model_from_info(model_info): 
    94     # type: (ModelInfo) -> SasviewModelType 
    9585    """ 
    9686    Convert *model_info* into a SasView model wrapper. 
    9787    """ 
    98     def __init__(self, multiplicity=None): 
     88    model_info['variant_info'] = None  # temporary hack for older sasview 
     89    def __init__(self, multiplicity=1): 
    9990        SasviewModel.__init__(self, multiplicity=multiplicity) 
    10091    attrs = _generate_model_attributes(model_info) 
    10192    attrs['__init__'] = __init__ 
    102     ConstructedModel = type(model_info.name, (SasviewModel,), attrs) # type: SasviewModelType 
     93    ConstructedModel = type(model_info['id'], (SasviewModel,), attrs) 
    10394    return ConstructedModel 
    10495 
    … …  
    113104    All the attributes should be immutable to avoid accidents. 
    114105    """ 
    115  
    116     # TODO: allow model to override axis labels input/output name/unit 
    117  
    118     # Process multiplicity 
    119     non_fittable = []  # type: List[str] 
    120     xlabel = model_info.profile_axes[0] if model_info.profile is not None else "" 
    121     variants = MultiplicityInfo(0, "", [], xlabel) 
    122     for p in model_info.parameters.kernel_parameters: 
    123         if p.name == model_info.control: 
    124             non_fittable.append(p.name) 
    125             variants = MultiplicityInfo( 
    126                 len(p.choices), p.name, p.choices, xlabel 
    127             ) 
    128             break 
    129         elif p.is_control: 
    130             non_fittable.append(p.name) 
    131             variants = MultiplicityInfo( 
    132                 int(p.limits[1]), p.name, p.choices, xlabel 
    133             ) 
    134             break 
    135  
    136     # Organize parameter sets 
    137     orientation_params = [] 
    138     magnetic_params = [] 
    139     fixed = [] 
    140     for p in model_info.parameters.user_parameters(): 
    141         if p.type == 'orientation': 
    142             orientation_params.append(p.name) 
    143             orientation_params.append(p.name+".width") 
    144             fixed.append(p.name+".width") 
    145         if p.type == 'magnetic': 
    146             orientation_params.append(p.name) 
    147             magnetic_params.append(p.name) 
    148             fixed.append(p.name+".width") 
    149  
    150     # Build class dictionary 
    151106    attrs = {}  # type: Dict[str, Any] 
    152107    attrs['_model_info'] = model_info 
    153     attrs['name'] = model_info.name 
    154     attrs['id'] = model_info.id 
    155     attrs['description'] = model_info.description 
    156     attrs['category'] = model_info.category 
    157     attrs['is_structure_factor'] = model_info.structure_factor 
    158     attrs['is_form_factor'] = model_info.ER is not None 
     108    attrs['name'] = model_info['name'] 
     109    attrs['id'] = model_info['id'] 
     110    attrs['description'] = model_info['description'] 
     111    attrs['category'] = model_info['category'] 
     112 
     113    # TODO: allow model to override axis labels input/output name/unit 
     114 
     115    #self.is_multifunc = False 
     116    non_fittable = []  # type: List[str] 
     117    variants = MultiplicityInfo(0, "", [], "") 
     118    attrs['is_structure_factor'] = model_info['structure_factor'] 
     119    attrs['is_form_factor'] = model_info['ER'] is not None 
    159120    attrs['is_multiplicity_model'] = variants[0] > 1 
    160121    attrs['multiplicity_info'] = variants 
     122 
     123    partype = model_info['partype'] 
     124    orientation_params = ( 
     125            partype['orientation'] 
     126            + [n + '.width' for n in partype['orientation']] 
     127            + partype['magnetic']) 
     128    magnetic_params = partype['magnetic'] 
     129    fixed = [n + '.width' for n in partype['pd-2d']] 
     130 
    161131    attrs['orientation_params'] = tuple(orientation_params) 
    162132    attrs['magnetic_params'] = tuple(magnetic_params) 
    163133    attrs['fixed'] = tuple(fixed) 
     134 
    164135    attrs['non_fittable'] = tuple(non_fittable) 
    165136 
    … …  
    221192    dispersion = None  # type: Dict[str, Any] 
    222193    #: units and limits for each parameter 
    223     details = None     # type: Dict[str, Sequence[Any]] 
    224     #                  # actual type is Dict[str, List[str, float, float]] 
    225     #: multiplicity value, or None if no multiplicity on the model 
     194    details = None     # type: Mapping[str, Tuple(str, float, float)] 
     195    #: multiplicity used, or None if no multiplicity controls 
    226196    multiplicity = None     # type: Optional[int] 
    227     #: memory for polydispersity array if using ArrayDispersion (used by sasview). 
    228     _persistency_dict = None # type: Dict[str, Tuple[np.ndarray, np.ndarray]] 
    229  
    230     def __init__(self, multiplicity=None): 
    231         # type: (Optional[int]) -> None 
    232  
    233         # TODO: _persistency_dict to persistency_dict throughout sasview 
    234         # TODO: refactor multiplicity to encompass variants 
    235         # TODO: dispersion should be a class 
    236         # TODO: refactor multiplicity info 
    237         # TODO: separate profile view from multiplicity 
    238         # The button label, x and y axis labels and scale need to be under 
    239         # the control of the model, not the fit page.  Maximum flexibility, 
    240         # the fit page would supply the canvas and the profile could plot 
    241         # how it wants, but this assumes matplotlib.  Next level is that 
    242         # we provide some sort of data description including title, labels 
    243         # and lines to plot. 
    244  
    245         # Get the list of hidden parameters given the mulitplicity 
    246         # Don't include multiplicity in the list of parameters 
     197 
     198    def __init__(self, multiplicity): 
     199        # type: () -> None 
     200        print("initializing", self.name) 
     201        #raise Exception("first initialization") 
     202        self._model = None 
     203 
     204        ## _persistency_dict is used by sas.perspectives.fitting.basepage 
     205        ## to store dispersity reference. 
     206        self._persistency_dict = {} 
     207 
    247208        self.multiplicity = multiplicity 
    248         if multiplicity is not None: 
    249             hidden = self._model_info.get_hidden_parameters(multiplicity) 
    250             hidden |= set([self.multiplicity_info.control]) 
    251         else: 
    252             hidden = set() 
    253  
    254         self._persistency_dict = {} 
     209 
    255210        self.params = collections.OrderedDict() 
    256211        self.dispersion = {} 
    257212        self.details = {} 
    258         for p in self._model_info.parameters.user_parameters(): 
    259             if p.name in hidden: 
    260                 continue 
     213 
     214        for p in self._model_info['parameters']: 
    261215            self.params[p.name] = p.default 
    262             self.details[p.id] = [p.units, p.limits[0], p.limits[1]] 
    263             if p.polydisperse: 
    264                 self.details[p.id+".width"] = [ 
    265                     "", 0.0, 1.0 if p.relative_pd else np.inf 
    266                 ] 
    267                 self.dispersion[p.name] = { 
    268                     'width': 0, 
    269                     'npts': 35, 
    270                     'nsigmas': 3, 
    271                     'type': 'gaussian', 
    272                 } 
     216            self.details[p.name] = [p.units] + p.limits 
     217 
     218        for name in self._model_info['partype']['pd-2d']: 
     219            self.dispersion[name] = { 
     220                'width': 0, 
     221                'npts': 35, 
     222                'nsigmas': 3, 
     223                'type': 'gaussian', 
     224            } 
    273225 
    274226    def __get_state__(self): 
    275         # type: () -> Dict[str, Any] 
    276227        state = self.__dict__.copy() 
    277228        state.pop('_model') 
    … …  
    282233 
    283234    def __set_state__(self, state): 
    284         # type: (Dict[str, Any]) -> None 
    285235        self.__dict__ = state 
    286236        self._model = None 
    287237 
    288238    def __str__(self): 
    289         # type: () -> str 
    290239        """ 
    291240        :return: string representation 
    … …  
    294243 
    295244    def is_fittable(self, par_name): 
    296         # type: (str) -> bool 
    297245        """ 
    298246        Check if a given parameter is fittable or not 
    … …  
    305253 
    306254 
     255    # pylint: disable=no-self-use 
    307256    def getProfile(self): 
    308         # type: () -> (np.ndarray, np.ndarray) 
    309257        """ 
    310258        Get SLD profile 
    … …  
    313261                beta is a list of the corresponding SLD values 
    314262        """ 
    315         args = [] # type: List[Union[float, np.ndarray]] 
    316         for p in self._model_info.parameters.kernel_parameters: 
    317             if p.id == self.multiplicity_info.control: 
    318                 args.append(float(self.multiplicity)) 
    319             elif p.length == 1: 
    320                 args.append(self.params.get(p.id, np.NaN)) 
    321             else: 
    322                 args.append([self.params.get(p.id+str(k), np.NaN) 
    323                              for k in range(1,p.length+1)]) 
    324         return self._model_info.profile(*args) 
     263        return None, None 
    325264 
    326265    def setParam(self, name, value): 
    327         # type: (str, float) -> None 
    328266        """ 
    329267        Set the value of a model parameter 
    … …  
    352290 
    353291    def getParam(self, name): 
    354         # type: (str) -> float 
    355292        """ 
    356293        Set the value of a model parameter 
    … …  
    376313 
    377314    def getParamList(self): 
    378         # type: () -> Sequence[str] 
    379315        """ 
    380316        Return a list of all available parameters for the model 
    381317        """ 
    382         param_list = list(self.params.keys()) 
     318        param_list = self.params.keys() 
    383319        # WARNING: Extending the list with the dispersion parameters 
    384320        param_list.extend(self.getDispParamList()) 
    … …  
    386322 
    387323    def getDispParamList(self): 
    388         # type: () -> Sequence[str] 
    389324        """ 
    390325        Return a list of polydispersity parameters for the model 
    391326        """ 
    392327        # TODO: fix test so that parameter order doesn't matter 
    393         ret = ['%s.%s' % (p.name.lower(), ext) 
    394                for p in self._model_info.parameters.user_parameters() 
    395                for ext in ('npts', 'nsigmas', 'width') 
    396                if p.polydisperse] 
     328        ret = ['%s.%s' % (d.lower(), p) 
     329               for d in self._model_info['partype']['pd-2d'] 
     330               for p in ('npts', 'nsigmas', 'width')] 
    397331        #print(ret) 
    398332        return ret 
    399333 
    400334    def clone(self): 
    401         # type: () -> "SasviewModel" 
    402335        """ Return a identical copy of self """ 
    403336        return deepcopy(self) 
    404337 
    405338    def run(self, x=0.0): 
    406         # type: (Union[float, (float, float), List[float]]) -> float 
    407339        """ 
    408340        Evaluate the model 
    … …  
    417349            # pylint: disable=unpacking-non-sequence 
    418350            q, phi = x 
    419             return self.calculate_Iq([q*math.cos(phi)], [q*math.sin(phi)])[0] 
    420         else: 
    421             return self.calculate_Iq([x])[0] 
     351            return self.calculate_Iq([q * math.cos(phi)], 
     352                                     [q * math.sin(phi)])[0] 
     353        else: 
     354            return self.calculate_Iq([float(x)])[0] 
    422355 
    423356 
    424357    def runXY(self, x=0.0): 
    425         # type: (Union[float, (float, float), List[float]]) -> float 
    426358        """ 
    427359        Evaluate the model in cartesian coordinates 
    … …  
    434366        """ 
    435367        if isinstance(x, (list, tuple)): 
    436             return self.calculate_Iq([x[0]], [x[1]])[0] 
    437         else: 
    438             return self.calculate_Iq([x])[0] 
     368            return self.calculate_Iq([float(x[0])], [float(x[1])])[0] 
     369        else: 
     370            return self.calculate_Iq([float(x)])[0] 
    439371 
    440372    def evalDistribution(self, qdist): 
    441         # type: (Union[np.ndarray, Tuple[np.ndarray, np.ndarray], List[np.ndarray]]) -> np.ndarray 
    442373        r""" 
    443374        Evaluate a distribution of q-values. 
    … …  
    470401            # Check whether we have a list of ndarrays [qx,qy] 
    471402            qx, qy = qdist 
    472             if not self._model_info.parameters.has_2d: 
     403            partype = self._model_info['partype'] 
     404            if not partype['orientation'] and not partype['magnetic']: 
    473405                return self.calculate_Iq(np.sqrt(qx ** 2 + qy ** 2)) 
    474406            else: 
    … …  
    483415                            % type(qdist)) 
    484416 
    485     def calculate_Iq(self, qx, qy=None): 
    486         # type: (Sequence[float], Optional[Sequence[float]]) -> np.ndarray 
     417    def calculate_Iq(self, *args): 
    487418        """ 
    488419        Calculate Iq for one set of q with the current parameters. 
    … …  
    495426        if self._model is None: 
    496427            self._model = core.build_model(self._model_info) 
    497         if qy is not None: 
    498             q_vectors = [np.asarray(qx), np.asarray(qy)] 
    499         else: 
    500             q_vectors = [np.asarray(qx)] 
    501         kernel = self._model.make_kernel(q_vectors) 
    502         pairs = [self._get_weights(p) 
    503                  for p in self._model_info.parameters.call_parameters] 
    504         call_details, weight, value = details.build_details(kernel, pairs) 
    505         result = kernel(call_details, weight, value, cutoff=self.cutoff) 
    506         kernel.release() 
     428        q_vectors = [np.asarray(q) for q in args] 
     429        fn = self._model.make_kernel(q_vectors) 
     430        pars = [self.params[v] for v in fn.fixed_pars] 
     431        pd_pars = [self._get_weights(p) for p in fn.pd_pars] 
     432        result = fn(pars, pd_pars, self.cutoff) 
     433        fn.q_input.release() 
     434        fn.release() 
    507435        return result 
    508436 
    509437    def calculate_ER(self): 
    510         # type: () -> float 
    511438        """ 
    512439        Calculate the effective radius for P(q)*S(q) 
    … …  
    514441        :return: the value of the effective radius 
    515442        """ 
    516         if self._model_info.ER is None: 
     443        ER = self._model_info.get('ER', None) 
     444        if ER is None: 
    517445            return 1.0 
    518446        else: 
    519             value, weight = self._dispersion_mesh() 
    520             fv = self._model_info.ER(*value) 
     447            values, weights = self._dispersion_mesh() 
     448            fv = ER(*values) 
    521449            #print(values[0].shape, weights.shape, fv.shape) 
    522             return np.sum(weight * fv) / np.sum(weight) 
     450            return np.sum(weights * fv) / np.sum(weights) 
    523451 
    524452    def calculate_VR(self): 
    525         # type: () -> float 
    526453        """ 
    527454        Calculate the volf ratio for P(q)*S(q) 
    … …  
    529456        :return: the value of the volf ratio 
    530457        """ 
    531         if self._model_info.VR is None: 
     458        VR = self._model_info.get('VR', None) 
     459        if VR is None: 
    532460            return 1.0 
    533461        else: 
    534             value, weight = self._dispersion_mesh() 
    535             whole, part = self._model_info.VR(*value) 
    536             return np.sum(weight * part) / np.sum(weight * whole) 
     462            values, weights = self._dispersion_mesh() 
     463            whole, part = VR(*values) 
     464            return np.sum(weights * part) / np.sum(weights * whole) 
    537465 
    538466    def set_dispersion(self, parameter, dispersion): 
    539         # type: (str, weights.Dispersion) -> Dict[str, Any] 
    540467        """ 
    541468        Set the dispersion object for a model parameter 
    … …  
    560487            from . import weights 
    561488            disperser = weights.dispersers[dispersion.__class__.__name__] 
    562             dispersion = weights.MODELS[disperser]() 
     489            dispersion = weights.models[disperser]() 
    563490            self.dispersion[parameter] = dispersion.get_pars() 
    564491        else: 
    … …  
    566493 
    567494    def _dispersion_mesh(self): 
    568         # type: () -> List[Tuple[np.ndarray, np.ndarray]] 
    569495        """ 
    570496        Create a mesh grid of dispersion parameters and weights. 
    … …  
    574500        parameter set in the vector. 
    575501        """ 
    576         pars = [self._get_weights(p) 
    577                 for p in self._model_info.parameters.call_parameters 
    578                 if p.type == 'volume'] 
    579         return details.dispersion_mesh(self._model_info, pars) 
     502        pars = self._model_info['partype']['volume'] 
     503        return core.dispersion_mesh([self._get_weights(p) for p in pars]) 
    580504 
    581505    def _get_weights(self, par): 
    582         # type: (Parameter) -> Tuple[np.ndarray, np.ndarray] 
    583506        """ 
    584507        Return dispersion weights for parameter 
    585508        """ 
    586         if par.name not in self.params: 
    587             if par.name == self.multiplicity_info.control: 
    588                 return [self.multiplicity], [] 
    589             else: 
    590                 return [np.NaN], [] 
    591         elif par.polydisperse: 
    592             dis = self.dispersion[par.name] 
    593             value, weight = weights.get_weights( 
    594                 dis['type'], dis['npts'], dis['width'], dis['nsigmas'], 
    595                 self.params[par.name], par.limits, par.relative_pd) 
    596             return value, weight / np.sum(weight) 
    597         else: 
    598             return [self.params[par.name]], [] 
     509        from . import weights 
     510        relative = self._model_info['partype']['pd-rel'] 
     511        limits = self._model_info['limits'] 
     512        dis = self.dispersion[par] 
     513        value, weight = weights.get_weights( 
     514            dis['type'], dis['npts'], dis['width'], dis['nsigmas'], 
     515            self.params[par], limits[par], par in relative) 
     516        return value, weight / np.sum(weight) 
     517 
    599518 
    600519def test_model(): 
    601     # type: () -> float 
    602520    """ 
    603521    Test that a sasview model (cylinder) can be run. 
    … …  
    607525    return cylinder.evalDistribution([0.1,0.1]) 
    608526 
    609 def test_rpa(): 
    610     # type: () -> float 
    611     """ 
    612     Test that a sasview model (cylinder) can be run. 
    613     """ 
    614     RPA = _make_standard_model('rpa') 
    615     rpa = RPA(3) 
    616     return rpa.evalDistribution([0.1,0.1]) 
    617  
    618527 
    619528def test_model_list(): 
    620     # type: () -> None 
    621529    """ 
    622530    Make sure that all models build as sasview models. 

  • sasmodels/sesans.py

    r7ae2b7f r2684d45  
    1212from __future__ import division 
    1313 
    14 import numpy as np  # type: ignore 
    15 from numpy import pi, exp  # type: ignore 
     14import numpy as np 
     15from numpy import pi, exp 
    1616from scipy.special import jv as besselj 
    1717#import direct_model.DataMixin as model 

  • sasmodels/weights.py

    ra936688 ra936688  
    33""" 
    44# TODO: include dispersion docs with the disperser models 
    5 from math import sqrt  # type: ignore 
    6 import numpy as np  # type: ignore 
    7 from scipy.special import gammaln  # type: ignore 
     5from math import sqrt 
     6import numpy as np 
     7from scipy.special import gammaln 
    88 
    99class Dispersion(object): 
Note: See TracChangeset for help on using the changeset viewer.