Changeset ff7119b in sasmodels


Ignore:
Timestamp:
Aug 26, 2014 10:27:06 PM (10 years ago)
Author:
Paul Kienzle <pkienzle@…>
Branches:
master, core_shell_microgels, costrafo411, magnetic_model, release_v0.94, release_v0.95, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests
Children:
5d4777d
Parents:
a7684e5
Message:

docu update

Files:
1 deleted
11 edited
1 moved

Legend:

Unmodified
Added
Removed

  • compare-new.py

    r13d86bc rff7119b  
    66import numpy as np 
    77 
    8 from sasmodels.core import BumpsModel, plot_data, tic, opencl_model, dll_model 
     8from sasmodels.bumps_model import BumpsModel, plot_data, tic 
     9from sasmodels.gen import opencl_model, dll_model 
    910 
    1011def sasview_model(modelname, **pars): 

  • sasmodels/alignment.py

    r14de349 rff7119b  
    33 
    44Some web sites say that maximizing performance for OpenCL code requires 
    5 aligning data on certain memory boundaries. 
     5aligning data on certain memory boundaries.  The following functions 
     6provide this service: 
    67 
    7 :func:`data_alignment` queries all devices in the OpenCL context, returning 
    8 the most restrictive alignment. 
     8:func:`align_data` aligns an existing array, returning a new array of the 
     9correct alignment. 
    910 
    10 :func:`align_data` aligns an existing array. 
     11:func:`align_empty` to create an empty array of the correct alignment. 
    1112 
    12 :func:`align_empty` to create a new array of the correct alignment. 
     13Set alignment to :func:`gpu.environment()` attribute *boundary*. 
    1314 
    1415Note:  This code is unused. So far, tests have not demonstrated any 
    … …  
    1718to decide whether it is really required. 
    1819""" 
    19  
    2020import numpy as np 
    21 import pyopencl as cl 
    22  
    23 def data_alignment(context): 
    24     """ 
    25     Return the desired byte alignment across all devices. 
    26     """ 
    27     # Note: rely on the fact that alignment will be 2^k 
    28     return max(d.min_data_type_align_size for d in context.devices) 
    2921 
    3022def align_empty(shape, dtype, alignment=128): 
     23    """ 
     24    Return an empty array aligned on the alignment boundary. 
     25    """ 
    3126    size = np.prod(shape) 
    3227    dtype = np.dtype(dtype) 
    … …  
    4035 
    4136def align_data(v, dtype, alignment=128): 
     37    """ 
     38    Return a copy of an array on the alignment boundary. 
     39    """ 
    4240    # if v is contiguous, aligned, and of the correct type then just return v 
    4341    view = align_empty(v.shape, dtype, alignment=alignment) 
    … …  
    4543    return view 
    4644 
    47 def work_group_size(context, kernel): 
    48     """ 
    49     Return the desired work group size for the context. 
    50     """ 
    51     max(kernel.get_work_group_info(cl.kernel_work_group_info.PREFERRED_WORK_GROUP_SIZE_MULTIPLE, d) 
    52         for d in context.devices) 
    53  
    54  

  • sasmodels/bumps_model.py

    r13d86bc rff7119b  
    1 #!/usr/bin/env python 
    2 # -*- coding: utf-8 -*- 
    3  
     1""" 
     2Sasmodels core. 
     3""" 
    44import sys, os 
    55import datetime 
    6 import warnings 
    76 
    87import numpy as np 
    98 
    109 
    11 def opencl_model(kernel_module, dtype="single"): 
    12     from sasmodels import gen, gpu 
    13  
    14     source, info = gen.make(kernel_module) 
    15     ## for debugging, save source to a .cl file, edit it, and reload as model 
    16     #open(modelname+'.cl','w').write(source) 
    17     #source = open(modelname+'.cl','r').read() 
    18     return gpu.GpuModel(source, info, dtype) 
    19  
    20  
    21 if sys.platform == 'darwin': 
    22     COMPILE = "gcc-mp-4.7 -shared -fPIC -std=c99 -fopenmp -O2 -Wall %s -o %s -lm -lgomp" 
    23 elif os.name == 'nt': 
    24     COMPILE = "gcc -shared -fPIC -std=c99 -fopenmp -O2 -Wall %s -o %s -lm" 
    25 else: 
    26     COMPILE = "cc -shared -fPIC -std=c99 -fopenmp -O2 -Wall %s -o %s -lm" 
    27 DLL_PATH = "/tmp" 
    28  
    29  
    30 def dll_path(info): 
    31     from os.path import join as joinpath, split as splitpath, splitext 
    32     basename = splitext(splitpath(info['filename'])[1])[0] 
    33     return joinpath(DLL_PATH, basename+'.so') 
    34  
    35  
    36 def dll_model(kernel_module): 
    37     import os 
    38     import tempfile 
    39     from sasmodels import gen, dll 
    40  
    41     source, info = gen.make(kernel_module) 
    42     dllpath = dll_path(info) 
    43     if not os.path.exists(dllpath) \ 
    44             or (os.path.getmtime(dllpath) < os.path.getmtime(info['filename'])): 
    45         # Replace with a proper temp file 
    46         srcfile = tempfile.mkstemp(suffix=".c",prefix="sas_"+info['name']) 
    47         open(srcfile, 'w').write(source) 
    48         os.system(COMPILE%(srcfile, dllpath)) 
    49         ## comment the following to keep the generated c file 
    50         #os.unlink(srcfile) 
    51     return dll.DllModel(dllpath, info) 
    52  
    53  
    54 TIC = None 
    5510def tic(): 
    56     global TIC 
     11    """ 
     12    Timer function. 
     13 
     14    Use "toc=tic()" to start the clock and "toc()" to measure 
     15    a time interval. 
     16    """ 
    5717    then = datetime.datetime.now() 
    58     TIC = lambda: (datetime.datetime.now()-then).total_seconds() 
    59     return TIC 
    60  
    61  
    62 def toc(): 
    63     return TIC() 
     18    return lambda: (datetime.datetime.now()-then).total_seconds() 
    6419 
    6520 
    6621def load_data(filename): 
     22    """ 
     23    Load data using a sasview loader. 
     24    """ 
    6725    from sans.dataloader.loader import Loader 
    6826    loader = Loader() 
    … …  
    7331 
    7432 
     33def empty_data1D(q): 
     34    """ 
     35    Create empty 1D data using the given *q* as the x value. 
     36 
     37    Resolutions dq/q is 5%. 
     38    """ 
     39 
     40    from sans.dataloader.data_info import Data1D 
     41 
     42    Iq = 100*np.ones_like(q) 
     43    dIq = np.sqrt(Iq) 
     44    data = Data1D(q, Iq, dx=0.05*q, dy=dIq) 
     45    data.filename = "fake data" 
     46    data.qmin, data.qmax = q.min(), q.max() 
     47    return data 
     48 
     49 
    7550def empty_data2D(qx, qy=None): 
     51    """ 
     52    Create empty 2D data using the given mesh. 
     53 
     54    If *qy* is missing, create a square mesh with *qy=qx*. 
     55 
     56    Resolution dq/q is 5%. 
     57    """ 
    7658    from sans.dataloader.data_info import Data2D, Detector 
    7759 
    … …  
    11496 
    11597 
    116 def empty_data1D(q): 
    117     from sans.dataloader.data_info import Data1D 
    118  
    119     Iq = 100*np.ones_like(q) 
    120     dIq = np.sqrt(Iq) 
    121     data = Data1D(q, Iq, dx=0.05*q, dy=dIq) 
    122     data.filename = "fake data" 
    123     data.qmin, data.qmax = q.min(), q.max() 
    124     return data 
    125  
    126  
    12798def set_beam_stop(data, radius, outer=None): 
     99    """ 
     100    Add a beam stop of the given *radius*.  If *outer*, make an annulus. 
     101    """ 
    128102    from sans.dataloader.manipulations import Ringcut 
    129103    if hasattr(data, 'qx_data'): 
    … …  
    138112 
    139113def set_half(data, half): 
     114    """ 
     115    Select half of the data, either "right" or "left". 
     116    """ 
    140117    from sans.dataloader.manipulations import Boxcut 
    141118    if half == 'right': 
    … …  
    146123 
    147124def set_top(data, max): 
     125    """ 
     126    Chop the top off the data, above *max*. 
     127    """ 
    148128    from sans.dataloader.manipulations import Boxcut 
    149129    data.mask += Boxcut(x_min=-np.inf, x_max=np.inf, y_min=-np.inf, y_max=max)(data) 
    … …  
    151131 
    152132def plot_data(data, iq, vmin=None, vmax=None, scale='log'): 
     133    """ 
     134    Plot the target value for the data.  This could be the data itself, 
     135    the theory calculation, or the residuals. 
     136 
     137    *scale* can be 'log' for log scale data, or 'linear'. 
     138    """ 
    153139    from numpy.ma import masked_array, masked 
    154140    import matplotlib.pyplot as plt 
    … …  
    172158 
    173159 
    174 def plot_result2D(data, theory, view='log'): 
     160def _plot_result1D(data, theory, view): 
     161    """ 
     162    Plot the data and residuals for 1D data. 
     163    """ 
     164    import matplotlib.pyplot as plt 
     165    from numpy.ma import masked_array, masked 
     166    #print "not a number",sum(np.isnan(data.y)) 
     167    #data.y[data.y<0.05] = 0.5 
     168    mdata = masked_array(data.y, data.mask) 
     169    mdata[np.isnan(mdata)] = masked 
     170    if view is 'log': 
     171        mdata[mdata <= 0] = masked 
     172    mtheory = masked_array(theory, mdata.mask) 
     173    mresid = masked_array((theory-data.y)/data.dy, mdata.mask) 
     174 
     175    plt.subplot(121) 
     176    plt.errorbar(data.x, mdata, yerr=data.dy) 
     177    plt.plot(data.x, mtheory, '-', hold=True) 
     178    plt.yscale(view) 
     179    plt.subplot(122) 
     180    plt.plot(data.x, mresid, 'x') 
     181    #plt.axhline(1, color='black', ls='--',lw=1, hold=True) 
     182    #plt.axhline(0, color='black', lw=1, hold=True) 
     183    #plt.axhline(-1, color='black', ls='--',lw=1, hold=True) 
     184 
     185 
     186def _plot_result2D(data, theory, view): 
     187    """ 
     188    Plot the data and residuals for 2D data. 
     189    """ 
    175190    import matplotlib.pyplot as plt 
    176191    resid = (theory-data.data)/data.err_data 
    … …  
    185200    plt.colorbar() 
    186201 
    187  
    188 def plot_result1D(data, theory, view='log'): 
    189     import matplotlib.pyplot as plt 
    190     from numpy.ma import masked_array, masked 
    191     #print "not a number",sum(np.isnan(data.y)) 
    192     #data.y[data.y<0.05] = 0.5 
    193     mdata = masked_array(data.y, data.mask) 
    194     mdata[np.isnan(mdata)] = masked 
    195     if view is 'log': 
    196         mdata[mdata <= 0] = masked 
    197     mtheory = masked_array(theory, mdata.mask) 
    198     mresid = masked_array((theory-data.y)/data.dy, mdata.mask) 
    199  
    200     plt.subplot(121) 
    201     plt.errorbar(data.x, mdata, yerr=data.dy) 
    202     plt.plot(data.x, mtheory, '-', hold=True) 
    203     plt.yscale(view) 
    204     plt.subplot(122) 
    205     plt.plot(data.x, mresid, 'x') 
    206     #plt.axhline(1, color='black', ls='--',lw=1, hold=True) 
    207     #plt.axhline(0, color='black', lw=1, hold=True) 
    208     #plt.axhline(-1, color='black', ls='--',lw=1, hold=True) 
     202def plot_result(data, theory, view='log'): 
     203    """ 
     204    Plot the data and residuals. 
     205    """ 
     206    if hasattr(data, 'qx_data'): 
     207        _plot_result2D(data, theory, view) 
     208    else: 
     209        _plot_result1D(data, theory, view) 
    209210 
    210211 
    211212class BumpsModel(object): 
     213    """ 
     214    Return a bumps wrapper for a SAS model. 
     215 
     216    *data* is the data to be fitted. 
     217 
     218    *model* is the SAS model, e.g., from :func:`gen.opencl_model`. 
     219 
     220    *cutoff* is the integration cutoff, which avoids computing the 
     221    the SAS model where the polydispersity weight is low. 
     222 
     223    Model parameters can be initialized with additional keyword 
     224    arguments, or by assigning to model.parameter_name.value. 
     225 
     226    The resulting bumps model can be used directly in a FitProblem call. 
     227    """ 
    212228    def __init__(self, data, model, cutoff=1e-5, **kw): 
    213229        from bumps.names import Parameter 
    214         from . import gpu 
    215230 
    216231        # interpret data 
    217         self.is2D = hasattr(data,'qx_data') 
    218232        self.data = data 
    219         if self.is2D: 
     233        if hasattr(data, 'qx_data'): 
    220234            self.index = (data.mask==0) & (~np.isnan(data.data)) 
    221235            self.iq = data.data[self.index] 
    … …  
    294308 
    295309    def plot(self, view='log'): 
    296         if self.is2D: 
    297             plot_result2D(self.data, self.theory(), view=view) 
    298         else: 
    299             plot_result1D(self.data, self.theory(), view=view) 
     310        plot_result(self.data, self.theory(), view=view) 
    300311 
    301312    def save(self, basename): 

  • sasmodels/dll.py

    rce27e21 rff7119b  
    2626    for single and 'd', 'float64' or 'double' for double.  Double precision 
    2727    is an optional extension which may not be available on all devices. 
     28 
     29    Call :meth:`release` when done with the kernel. 
    2830    """ 
    2931    def __init__(self, dllpath, info): 
    … …  
    6971        return DllInput(q_vectors) 
    7072 
     73    def release(self): 
     74        pass # TODO: should release the dll 
     75 
    7176 
    7277class DllInput(object): 
    … …  
    100105 
    101106class DllKernel(object): 
     107    """ 
     108    Callable SAS kernel. 
     109 
     110    *kernel* is the DllKernel object to call. 
     111 
     112    *info* is the module information 
     113 
     114    *input* is the DllInput q vectors at which the kernel should be 
     115    evaluated. 
     116 
     117    The resulting call method takes the *pars*, a list of values for 
     118    the fixed parameters to the kernel, and *pd_pars*, a list of (value,weight) 
     119    vectors for the polydisperse parameters.  *cutoff* determines the 
     120    integration limits: any points with combined weight less than *cutoff* 
     121    will not be calculated. 
     122 
     123    Call :meth:`release` when done with the kernel instance. 
     124    """ 
    102125    def __init__(self, kernel, info, input): 
    103126        self.input = input 
    104127        self.kernel = kernel 
    105         self.info = info 
    106128        self.res = np.empty(input.nq, input.dtype) 
    107129        dim = '2d' if input.is_2D else '1d' 

  • sasmodels/gen.py

    ra7684e5 rff7119b  
    5555them with the desired polydispersity. 
    5656 
     57The available kernel parameters are defined as a list, with each parameter 
     58defined as a sublist with the following elements: 
     59 
     60    *name* is the name that will be used in the call to the kernel 
     61    function and the name that will be displayed to the user.  Names 
     62    should be lower case, with words separated by underscore.  If 
     63    acronyms are used, the whole acronym should be upper case. 
     64 
     65    *units* should be one of *degrees* for angles, *Ang* for lengths, 
     66    *1e-6/Ang^2* for SLDs. 
     67 
     68    *default value* will be the initial value for  the model when it 
     69    is selected, or when an initial value is not otherwise specified. 
     70 
     71    [*lb*, *ub*] are the hard limits on the parameter value, used to limit 
     72    the polydispersity density function.  In the fit, the parameter limits 
     73    given to the fit are the limits  on the central value of the parameter. 
     74    If there is polydispersity, it will evaluate parameter values outside 
     75    the fit limits, but not outside the hard limits specified in the model. 
     76    If there are no limits, use +/-inf imported from numpy. 
     77 
     78    *type* indicates how the parameter will be used.  "volume" parameters 
     79    will be used in all functions.  "orientation" parameters will be used 
     80    in *Iqxy* and *Imagnetic*.  "magnetic* parameters will be used in 
     81    *Imagnetic* only.  If *type* is the empty string, the parameter will 
     82    be used in all of *Iq*, *Iqxy* and *Imagnetic*. 
     83 
     84    *description* is a short description of the parameter.  This will 
     85    be displayed in the parameter table and used as a tool tip for the 
     86    parameter value in the user interface. 
     87 
    5788The kernel module must set variables defining the kernel meta data: 
    5889 
    … …  
    6596    while the model parameters are being edited. 
    6697 
    67     *parameters* is a list of parameters.  Each parameter is itself a 
    68     list containing *name*, *units*, *default value*, 
    69     [*lower bound*, *upper bound*], *type* and *description*. 
     98    *parameters* is the list of parameters.  Parameters in the kernel 
     99    functions must appear in the same order as they appear in the 
     100    parameters list.  Two additional parameters, *scale* and *background* 
     101    are added to the beginning of the parameter list.  They will show up 
     102    in the documentation as model parameters, but they are never sent to 
     103    the kernel functions. 
    70104 
    71105    *source* is the list of C-99 source files that must be joined to 
    … …  
    78112    *VR* is a python function defining the volume ratio.  If it is not 
    79113    present, the volume ratio is 1. 
     114 
     115An *info* dictionary is constructed from the kernel meta data and 
     116returned to the caller.  It includes the additional fields 
     117 
     118 
     119The model evaluator, function call sequence consists of q inputs and the return vector, 
     120followed by the loop value/weight vector, followed by the values for 
     121the non-polydisperse parameters, followed by the lengths of the 
     122polydispersity loops.  To construct the call for 1D models, the 
     123categories *fixed-1d* and *pd-1d* list the names of the parameters 
     124of the non-polydisperse and the polydisperse parameters respectively. 
     125Similarly, *fixed-2d* and *pd-2d* provide parameter names for 2D models. 
     126The *pd-rel* category is a set of those parameters which give 
     127polydispersitiy as a portion of the value (so a 10% length dispersity 
     128would use a polydispersity value of 0.1) rather than absolute 
     129dispersity such as an angle plus or minus 15 degrees. 
     130 
     131The *volume* category lists the volume parameters in order for calls 
     132to volume within the kernel (used for volume normalization) and for 
     133calls to ER and VR for effective radius and volume ratio respectively. 
     134 
     135The *orientation* and *magnetic* categories list the orientation and 
     136magnetic parameters.  These are used by the sasview interface.  The 
     137blank category is for parameters such as scale which don't have any 
     138other marking. 
    80139 
    81140The doc string at the start of the kernel module will be used to 
    … …  
    86145file names and extensions. 
    87146 
    88 Parameters are defined as follows: 
    89  
    90     *name* is the name that will be used in the call to the kernel 
    91     function and the name that will be displayed to the user.  Names 
    92     should be lower case, with words separated by underscore.  If 
    93     acronyms are used, the whole acronym should be upper case. 
    94  
    95     *units* should be one of *degrees* for angles, *Ang* for lengths, 
    96     *1e-6/Ang^2* for SLDs. 
    97  
    98     *default value* will be the initial value for  the model when it 
    99     is selected, or when an initial value is not otherwise specified. 
    100  
    101     *limits* are the valid bounds of the parameter, used to limit the 
    102     polydispersity density function.   In the fit, the parameter limits 
    103     given to the fit are the limits  on the central value of the parameter. 
    104     If there is polydispersity, it will evaluate parameter values outside 
    105     the fit limits, but not outside the hard limits specified in the model. 
    106     If there are no limits, use +/-inf imported from numpy. 
    107  
    108     *type* indicates how the parameter will be used.  "volume" parameters 
    109     will be used in all functions.  "orientation" parameters will be used 
    110     in *Iqxy* and *Imagnetic*.  "magnetic* parameters will be used in 
    111     *Imagnetic* only.  If *type* is none, the parameter will be used in 
    112     all of *Iq*, *Iqxy* and *Imagnetic*.  This will probably be a 
    113     is the empty string if the parameter is used in all model calculations, 
    114     it is "volu 
    115  
    116     *description* is a short description of the parameter.  This will 
    117     be displayed in the parameter table and used as a tool tip for the 
    118     parameter value in the user interface. 
    119147 
    120148The function :func:`make` loads the metadata from the module and returns 
    121 the kernel source.  The function :func:`doc` extracts 
     149the kernel source.  The function :func:`doc` extracts the doc string 
     150and adds the parameter table to the top.  The function :func:`sources` 
     151returns a list of files required by the model. 
    122152""" 
    123153 
    124154# TODO: identify model files which have changed since loading and reload them. 
    125155 
    126 __all__ = ["make, doc"] 
     156__all__ = ["make, doc", "sources"] 
    127157 
    128158import os.path 
    … …  
    158188    ] 
    159189 
     190# Minimum width for a default value (this is shorter than the column header 
     191# width, so will be ignored). 
    160192PARTABLE_VALUE_WIDTH = 10 
    161193 
    162194# Header included before every kernel. 
    163 # This makes sure that the appropriate math constants are defined, and 
     195# This makes sure that the appropriate math constants are defined, and does 
     196# whatever is required to make the kernel compile as pure C rather than 
     197# as an OpenCL kernel. 
    164198KERNEL_HEADER = """\ 
    165199// GENERATED CODE --- DO NOT EDIT --- 
    … …  
    231265    } 
    232266 
    233 # Generic kernel template for opencl/openmp. 
     267# Generic kernel template for the polydispersity loop. 
    234268# This defines the opencl kernel that is available to the host.  The same 
    235269# structure is used for Iq and Iqxy kernels, so extra flexibility is needed 
    … …  
    333367 
    334368def kernel_name(info, is_2D): 
     369    """ 
     370    Name of the exported kernel symbol. 
     371    """ 
    335372    return info['name'] + "_" + ("Iqxy" if is_2D else "Iq") 
    336373 
    … …  
    431468 
    432469def make_partable(info): 
     470    """ 
     471    Generate the parameter table to include in the sphinx documentation. 
     472    """ 
    433473    pars = info['parameters'] 
    434474    column_widths = [ 
    … …  
    467507    raise ValueError("%r not found in %s"%(filename, search_path)) 
    468508 
    469 def make_model(search_path, info): 
     509def sources(info): 
     510    """ 
     511    Return a list of the sources file paths for the module. 
     512    """ 
     513    from os.path import abspath, dirname, join as joinpath 
     514    search_path = [ dirname(info['filename']), 
     515                    abspath(joinpath(dirname(__file__),'models')) ] 
     516    return [_search(search_path) for f in info['source']] 
     517 
     518def make_model(info): 
     519    """ 
     520    Generate the code for the kernel defined by info, using source files 
     521    found in the given search path. 
     522    """ 
    470523    kernel_Iq = make_kernel(info, is_2D=False) 
    471524    kernel_Iqxy = make_kernel(info, is_2D=True) 
    472     source = [open(_search(search_path, f)).read() for f in info['source']] 
     525    source = [open(f).read() for f in sources(info)] 
    473526    kernel = "\n\n".join([KERNEL_HEADER]+source+[kernel_Iq, kernel_Iqxy]) 
    474527    return kernel 
    … …  
    479532 
    480533    Returns a dictionary of categories. 
    481  
    482     The function call sequence consists of q inputs and the return vector, 
    483     followed by the loop value/weight vector, followed by the values for 
    484     the non-polydisperse parameters, followed by the lengths of the 
    485     polydispersity loops.  To construct the call for 1D models, the 
    486     categories *fixed-1d* and *pd-1d* list the names of the parameters 
    487     of the non-polydisperse and the polydisperse parameters respectively. 
    488     Similarly, *fixed-2d* and *pd-2d* provide parameter names for 2D models. 
    489     The *pd-rel* category is a set of those parameters which give 
    490     polydispersitiy as a portion of the value (so a 10% length dispersity 
    491     would use a polydispersity value of 0.1) rather than absolute 
    492     dispersity such as an angle plus or minus 15 degrees. 
    493  
    494     The *volume* category lists the volume parameters in order for calls 
    495     to volume within the kernel (used for volume normalization) and for 
    496     calls to ER and VR for effective radius and volume ratio respectively. 
    497  
    498     The *orientation* and *magnetic* categories list the orientation and 
    499     magnetic parameters.  These are used by the sasview interface.  The 
    500     blank category is for parameters such as scale which don't have any 
    501     other marking. 
    502534    """ 
    503535    partype = { 
    … …  
    535567    """ 
    536568    # TODO: allow Iq and Iqxy to be defined in python 
    537     from os.path import abspath, dirname, join as joinpath 
     569    from os.path import abspath 
    538570    #print kernelfile 
    539571    info = dict( 
    … …  
    550582    info['partype'] = categorize_parameters(info['parameters']) 
    551583 
    552     search_path = [ dirname(info['filename']), 
    553                     abspath(joinpath(dirname(__file__),'models')) ] 
    554     source = make_model(search_path, info) 
     584    source = make_model(info) 
    555585 
    556586    return source, info 
    … …  
    565595                 doc=kernel_module.__doc__) 
    566596    return DOC_HEADER%subst 
     597 
     598# Compiler platform details 
     599if sys.platform == 'darwin': 
     600    COMPILE = "gcc-mp-4.7 -shared -fPIC -std=c99 -fopenmp -O2 -Wall %s -o %s -lm -lgomp" 
     601elif os.name == 'nt': 
     602    COMPILE = "gcc -shared -fPIC -std=c99 -fopenmp -O2 -Wall %s -o %s -lm" 
     603else: 
     604    COMPILE = "cc -shared -fPIC -std=c99 -fopenmp -O2 -Wall %s -o %s -lm" 
     605DLL_PATH = "/tmp" 
     606 
     607 
     608def dll_path(info): 
     609    """ 
     610    Path to the compiled model defined by *info*. 
     611    """ 
     612    from os.path import join as joinpath, split as splitpath, splitext 
     613    basename = splitext(splitpath(info['filename'])[1])[0] 
     614    return joinpath(DLL_PATH, basename+'.so') 
     615 
     616 
     617def dll_model(kernel_module, dtype=None): 
     618    """ 
     619    Load the compiled model defined by *kernel_module*. 
     620 
     621    Recompile if any files are newer than the model file. 
     622 
     623    *dtype* is ignored.  Compiled files are always double. 
     624 
     625    The DLL is not loaded until the kernel is called so models an 
     626    be defined without using too many resources. 
     627    """ 
     628    import tempfile 
     629    from sasmodels import dll 
     630 
     631    source, info = make(kernel_module) 
     632    source_files = sources(info) + [info['filename']] 
     633    newest = max(os.path.getmtime(f) for f in source_files) 
     634    dllpath = dll_path(info) 
     635    if not os.path.exists(dllpath) or os.path.getmtime(dllpath)<newest: 
     636        # Replace with a proper temp file 
     637        srcfile = tempfile.mkstemp(suffix=".c",prefix="sas_"+info['name']) 
     638        open(srcfile, 'w').write(source) 
     639        os.system(COMPILE%(srcfile, dllpath)) 
     640        ## comment the following to keep the generated c file 
     641        os.unlink(srcfile) 
     642    return dll.DllModel(dllpath, info) 
     643 
     644 
     645def opencl_model(kernel_module, dtype="single"): 
     646    """ 
     647    Load the OpenCL model defined by *kernel_module*. 
     648 
     649    Access to the OpenCL device is delayed until the kernel is called 
     650    so models can be defined without using too many resources. 
     651    """ 
     652    from sasmodels import gpu 
     653 
     654    source, info = make(kernel_module) 
     655    ## for debugging, save source to a .cl file, edit it, and reload as model 
     656    #open(modelname+'.cl','w').write(source) 
     657    #source = open(modelname+'.cl','r').read() 
     658    return gpu.GpuModel(source, info, dtype) 
    567659 
    568660 

  • sasmodels/gpu.py

    r1780d59 rff7119b  
    242242 
    243243class GpuKernel(object): 
     244    """ 
     245    Callable SAS kernel. 
     246 
     247    *kernel* is the GpuKernel object to call. 
     248 
     249    *info* is the module information 
     250 
     251    *input* is the DllInput q vectors at which the kernel should be 
     252    evaluated. 
     253 
     254    The resulting call method takes the *pars*, a list of values for 
     255    the fixed parameters to the kernel, and *pd_pars*, a list of (value,weight) 
     256    vectors for the polydisperse parameters.  *cutoff* determines the 
     257    integration limits: any points with combined weight less than *cutoff* 
     258    will not be calculated. 
     259 
     260    Call :meth:`release` when done with the kernel instance. 
     261    """ 
    244262    def __init__(self, kernel, info, input): 
    245263        self.input = input 

  • sasmodels/models/README.rst

    ra7684e5 rff7119b  
    3333will mean extending the data handler to support multiple cross sections 
    3434in the same data set. 
     35 
     36Need to write code to generate the polydispersity loops in python for 
     37kernels that are only implemented in python. 

  • sasmodels/models/cylinder.c

    ra7684e5 rff7119b  
    1414    real length) 
    1515{ 
    16     const real h = REAL(0.5)*length; 
     16    const real halflength = REAL(0.5)*length; 
    1717    real summ = REAL(0.0); 
     18    // real lower=0, upper=M_PI_2; 
    1819    for (int i=0; i<76 ;i++) { 
    19         //const real zi = ( Gauss76Z[i]*(uplim-lolim) + uplim + lolim )/2.0; 
     20        // translate a point in [-1,1] to a point in [lower,upper] 
     21        //const real zi = ( Gauss76Z[i]*(upper-lower) + upper + lower )/2.0; 
    2022        const real zi = REAL(0.5)*(Gauss76Z[i]*M_PI_2 + M_PI_2); 
    21         summ += Gauss76Wt[i] * CylKernel(q, radius, h, zi); 
     23        summ += Gauss76Wt[i] * CylKernel(q, radius, halflength, zi); 
    2224    } 
    23     //const real form = (uplim-lolim)/2.0*summ; 
     25    // translate dx in [-1,1] to dx in [lower,upper] 
     26    //const real form = (upper-lower)/2.0*summ; 
    2427    const real form = summ * M_PI_4; 
    2528 

  • sasmodels/models/cylinder_clone.c

    r32c160a rff7119b  
    4545    // # The following correction factor exists in sasview, but it can't be 
    4646    // # right, so we are leaving it out for now. 
    47     // const real correction = fabs(cn)*M_PI_2; 
     47    const real spherical_integration = fabs(cn)*M_PI_2; 
    4848    const real q = sqrt(qx*qx+qy*qy); 
    4949    const real cos_val = cn*cos(cyl_phi*M_PI_180)*(qx/q) + sn*(qy/q); 
    … …  
    6565    // The additional volume factor is for polydisperse volume normalization. 
    6666    const real s = (sldCyl - sldSolv) * form_volume(radius, length); 
    67     return REAL(1.0e8) * form * s * s; // * correction; 
     67    return REAL(1.0e8) * form * s * s * spherical_integration; 
    6868} 

  • sasmodels/models/lib/cylkernel.c

    r14de349 rff7119b  
    11// Compute the form factor for a cylinder 
    2 // qq is the q-value for the calculation (1/A) 
     2//     F^2(q) sin a 
     3// given 
     4//     F(q) = sin(Q L/2 cos a)/(Q L/2 cos a) * 2 J1(Q R sin a) / (Q R sin a) 
     5// q is the q-value for the calculation (1/A) 
    36// radius is the radius of the cylinder (A) 
    4 // h is the HALF-LENGTH of the cylinder = L/2 (A) 
    5 real CylKernel(real q, real radius, real h, real theta); 
    6 real CylKernel(real q, real radius, real h, real theta) 
     7// halflength is the HALF-LENGTH of the cylinder = L/2 (A) 
     8real CylKernel(real q, real radius, real halflength, real theta); 
     9real CylKernel(real q, real radius, real halflength, real theta) 
    710{ 
    811    real sn, cn; 
    912    SINCOS(theta, sn, cn); 
    1013    const real besarg = q*radius*sn; 
    11     const real siarg = q*h*cn; 
     14    const real siarg = q*halflength*cn; 
    1215    // lim_{x->0} J1(x)/x = 1/2,  lim_{x->0} sin(x)/x = 1 
    1316    const real bj = (besarg == REAL(0.0) ? REAL(0.5) : J1(besarg)/besarg); 
    1417    const real si = (siarg == REAL(0.0) ? REAL(1.0) : sin(siarg)/siarg); 
    15     return REAL(4.0)*sin(theta)*bj*bj*si*si; 
     18    return REAL(4.0)*sn*bj*bj*si*si; 
    1619} 

  • sasmodels/sasview_model.py

    ra7684e5 rff7119b  
    11import math 
    22from copy import deepcopy 
     3import warnings 
    34 
    45import numpy as np 
    56 
     7try: 
     8    import pyopencl 
     9    from .gen import opencl_model as load_model 
     10except ImportError: 
     11    warnings.warn("OpenCL not available --- using ctypes instead") 
     12    from .gen import dll_model as load_model 
     13 
     14 
    615def make_class(kernel_module, dtype='single'): 
    7     from .core import opencl_model 
    8     model =  opencl_model(kernel_module, dtype=dtype) 
     16    """ 
     17    Load the sasview model defined in *kernel_module*. 
     18    :param kernel_module: 
     19    :param dtype: 
     20    :return: 
     21    """ 
     22    model =  load_model(kernel_module, dtype=dtype) 
    923    def __init__(self, multfactor=1): 
    1024        SasviewModel.__init__(self, model) 
    … …  
    242256 
    243257    def calculate_Iq(self, *args): 
     258        """ 
     259        Calculate Iq for one set of q with the current parameters. 
     260 
     261        If the model is 1D, use *q*.  If 2D, use *qx*, *qy*. 
     262 
     263        This should NOT be used for fitting since it copies the *q* vectors 
     264        to the card for each evaluation. 
     265        """ 
    244266        q_vectors = [np.asarray(q) for q in args] 
    245267        fn = self._model(self._model.make_input(q_vectors)) 

  • sasmodels/weights.py

    r1780d59 rff7119b  
     1""" 
     2SAS distributions for polydispersity. 
     3""" 
    14# TODO: include dispersion docs with the disperser models 
    25from math import sqrt 
    … …  
    116119 
    117120 
     121# dispersion name -> disperser lookup table. 
    118122models = dict((d.type,d) for d in ( 
    119123    GaussianDispersion, RectangleDispersion, 
    … …  
    123127 
    124128def get_weights(disperser, n, width, nsigmas, value, limits, relative): 
     129    """ 
     130    Return the set of values and weights for a polydisperse parameter. 
     131 
     132    *disperser* is the name of the disperser. 
     133 
     134    *n* is the number of points in the weight vector. 
     135 
     136    *width* is the width of the disperser distribution. 
     137 
     138    *nsigmas* is the number of sigmas to span for the dispersion convolution. 
     139 
     140    *value* is the value of the parameter in the model. 
     141 
     142    *limits* is [lb, ub], the lower and upper bound of the weight value. 
     143 
     144    *relative* is true if *width* is defined in proportion to the value 
     145    of the parameter, and false if it is an absolute width. 
     146 
     147    Returns *(value,weight)*, where *value* and *weight* are vectors. 
     148    """ 
    125149    cls = models[disperser] 
    126150    obj = cls(n, width, nsigmas) 
Note: See TracChangeset for help on using the changeset viewer.