Changeset 32c160a in sasmodels


Ignore:
Timestamp:
Aug 25, 2014 12:55:08 AM (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:
13d86bc
Parents:
1f21edf
Message:

support ER/VR python kernels; move metadata to python

Files:
3 added
5 edited

Legend:

Unmodified
Added
Removed

  • CylinderModel.py

    rce27e21 r32c160a  
    55""" 
    66from sasmodels.sasview_model import make_class 
    7 CylinderModel = make_class('cylinder_clone','CylinderModel') 
     7from sasmodels.models import cylinder_clone as cylinder 
     8CylinderModel = make_class(cylinder, dtype='single') 

  • sasmodels/core.py

    rce27e21 r32c160a  
    88import numpy as np 
    99 
    10 def load_model(modelname): 
    11     from os.path import abspath, join as joinpath, dirname 
    12     from sasmodels import gen 
    13     modelpath = abspath(joinpath(dirname(gen.__file__), 'models', 
    14                                  modelname+'.c')) 
    15     return gen.make(modelpath) 
    16  
    17  
    18  
    19 def opencl_model(modelname, dtype="single"): 
    20     from sasmodels import gpu 
    21  
    22     source, info, _ = load_model(modelname) 
    23     # for debugging, save source to a .cl file, edit it, and reload as model 
     10 
     11def 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 
    2416    #open(modelname+'.cl','w').write(source) 
    2517    #source = open(modelname+'.cl','r').read() 
    … …  
    4234 
    4335 
    44 def dll_model(modelname): 
     36def dll_model(kernel_module): 
    4537    import os 
    46     from sasmodels import dll 
    47  
    48     source, info, _ = load_model(modelname) 
     38    import tempfile 
     39    from sasmodels import gen, dll 
     40 
     41    source, info = gen.make(kernel_module) 
    4942    dllpath = dll_path(info) 
    5043    if not os.path.exists(dllpath) \ 
    5144            or (os.path.getmtime(dllpath) < os.path.getmtime(info['filename'])): 
    5245        # Replace with a proper temp file 
    53         srcfile = '/tmp/%s.c'%modelname 
     46        srcfile = tempfile.mkstemp(suffix=".c",prefix="sas_"+info['name']) 
    5447        open(srcfile, 'w').write(source) 
    5548        os.system(COMPILE%(srcfile, dllpath)) 
     49        ## comment the following to keep the generated c file 
     50        #os.unlink(srcfile) 
    5651    return dll.DllModel(dllpath, info) 
    5752 

  • sasmodels/gen.py

    rce27e21 r32c160a  
    316316    return "\n".join(lines) 
    317317 
    318 def make_doc(kernelfile, info, doc): 
    319     doc = doc%{'parameters': make_partable(info)} 
    320     return doc 
    321  
    322 def make_model(kernelfile, info, source): 
     318def _search(search_path, filename): 
     319    """ 
     320    Find *filename* in *search_path*. 
     321 
     322    Raises ValueError if file does not exist. 
     323    """ 
     324    for path in search_path: 
     325        target = os.path.join(path, filename) 
     326        if os.path.exists(target): 
     327            return target 
     328    raise ValueError("%r not found in %s"%(filename, search_path)) 
     329 
     330def make_model(search_path, info): 
    323331    kernel_Iq = make_kernel(info, is_2D=False) 
    324332    kernel_Iqxy = make_kernel(info, is_2D=True) 
    325     path = os.path.dirname(kernelfile) 
    326     extra = [open("%s/%s"%(path,f)).read() for f in info['include']] 
    327     kernel = "\n\n".join([KERNEL_HEADER]+extra+[source, kernel_Iq, kernel_Iqxy]) 
     333    source = [open(_search(search_path, f)).read() for f in info['source']] 
     334    kernel = "\n\n".join([KERNEL_HEADER]+source+[kernel_Iq, kernel_Iqxy]) 
    328335    return kernel 
    329  
    330 def parse_file(kernelfile): 
    331     source = open(kernelfile).read() 
    332  
    333     # select parameters out of the source file 
    334     parts = source.split("PARAMETERS") 
    335     if len(parts) != 3: 
    336         raise ValueError("PARAMETERS block missing from %r"%kernelfile) 
    337     info_source = parts[1].strip() 
    338     try: 
    339         info = relaxed_loads(info_source) 
    340     except: 
    341         print "in json text:" 
    342         print "\n".join("%2d: %s"%(i+1,s) 
    343                         for i,s in enumerate(info_source.split('\n'))) 
    344         raise 
    345         #raise ValueError("PARAMETERS block could not be parsed from %r"%kernelfile) 
    346  
    347     # select documentation out of the source file 
    348     parts = source.split("DOCUMENTATION") 
    349     if len(parts) == 3: 
    350         doc = make_doc(kernelfile, info, parts[1].strip()) 
    351     elif len(parts) == 1: 
    352         raise ValueError("DOCUMENTATION block is missing from %r"%kernelfile) 
    353     else: 
    354         raise ValueError("DOCUMENTATION block incorrect from %r"%kernelfile) 
    355  
    356     return source, info, doc 
    357336 
    358337def categorize_parameters(pars): 
    … …  
    408387    return partype 
    409388 
    410 def make(kernelfile): 
     389def make(kernel_module): 
    411390    """ 
    412391    Build an OpenCL function from the source in *kernelfile*. 
    … …  
    415394    will be a JSON definition containing 
    416395    """ 
     396    # TODO: allow Iq and Iqxy to be defined in python 
     397    from os.path import abspath, dirname, join as joinpath 
    417398    #print kernelfile 
    418     source, info, doc = parse_file(kernelfile) 
    419     info['filename'] = kernelfile 
    420     info['parameters'] = COMMON_PARAMETERS + info['parameters'] 
     399    info = dict( 
     400        filename = abspath(kernel_module.__file__), 
     401        name = kernel_module.name, 
     402        title = kernel_module.title, 
     403        source = kernel_module.source, 
     404        description = kernel_module.description, 
     405        parameters = COMMON_PARAMETERS + kernel_module.parameters, 
     406        ER = getattr(kernel_module, 'ER', None), 
     407        VR = getattr(kernel_module, 'VR', None), 
     408        ) 
     409    info['limits'] = dict((p[0],p[3]) for p in info['parameters']) 
    421410    info['partype'] = categorize_parameters(info['parameters']) 
    422     info['limits'] = dict((p[0],p[3]) for p in info['parameters']) 
    423     doc = make_doc(kernelfile, info, doc) 
    424     model = make_model(kernelfile, info, source) 
    425     return model, info, doc 
     411 
     412    search_path = [ dirname(info['filename']), 
     413                    abspath(joinpath(dirname(__file__),'models')) ] 
     414    source = make_model(search_path, info) 
     415 
     416    return source, info 
    426417 
    427418 

  • sasmodels/models/cylinder_clone.c

    r1f21edf r32c160a  
    1 /* PARAMETERS 
    2 { 
    3 name: "CylinderModel", 
    4 title: "Cylinder with uniform scattering length density", 
    5 include: [ "lib/J1.c", "lib/gauss76.c", "lib/cylkernel.c"], 
    6 parameters: [ 
    7    // [ "name", "units", default, [lower, upper], "type", "description" ], 
    8    [ "sldCyl", "1e-6/Ang^2", 4e-6, [-Infinity,Infinity], "", 
    9      "Cylinder scattering length density" ], 
    10    [ "sldSolv", "1e-6/Ang^2", 1e-6, [-Infinity,Infinity], "", 
    11      "Solvent scattering length density" ], 
    12    [ "radius", "Ang",  20, [0, Infinity], "volume", 
    13      "Cylinder radius" ], 
    14    [ "length", "Ang",  400, [0, Infinity], "volume", 
    15      "Cylinder length" ], 
    16    [ "cyl_theta", "degrees", 60, [-Infinity, Infinity], "orientation", 
    17      "In plane angle" ], 
    18    [ "cyl_phi", "degrees", 60, [-Infinity, Infinity], "orientation", 
    19      "Out of plane angle" ], 
    20 ], 
    21 description: "f(q)= 2*(sldCyl - sldSolv)*V*sin(qLcos(alpha/2))/[qLcos(alpha/2)]*J1(qRsin(alpha/2))/[qRsin(alpha)]", 
    22 } 
    23 PARAMETERS END 
    24  
    25 DOCUMENTATION 
    26 .. _CylinderModel: 
    27  
    28 CylinderModel 
    29 ============= 
    30  
    31 This model provides the form factor for a right circular cylinder with uniform 
    32 scattering length density. The form factor is normalized by the particle volume. 
    33  
    34 For information about polarised and magnetic scattering, click here_. 
    35  
    36 Definition 
    37 ---------- 
    38  
    39 The output of the 2D scattering intensity function for oriented cylinders is 
    40 given by (Guinier, 1955) 
    41  
    42 .. math:: 
    43  
    44     P(q,\alpha) = \frac{\text{scale}}{V}f^2(q) + \text{bkg} 
    45  
    46 where 
    47  
    48 .. math:: 
    49  
    50     f(q) = 2 (\Delta \rho) V 
    51            \frac{\sin (q L/2 \cos \alpha)}{q L/2 \cos \alpha} 
    52            \frac{J_1 (q r \sin \alpha)}{q r \sin \alpha} 
    53  
    54 and $\alpha$ is the angle between the axis of the cylinder and $\vec q$, $V$ 
    55 is the volume of the cylinder, $L$ is the length of the cylinder, $r$ is the 
    56 radius of the cylinder, and $d\rho$ (contrast) is the scattering length density 
    57 difference between the scatterer and the solvent. $J_1$ is the first order 
    58 Bessel function. 
    59  
    60 To provide easy access to the orientation of the cylinder, we define the 
    61 axis of the cylinder using two angles $\theta$ and $\phi$. Those angles 
    62 are defined in Figure :num:`figure #CylinderModel-orientation`. 
    63  
    64 .. _CylinderModel-orientation: 
    65  
    66 .. figure:: img/image061.JPG 
    67  
    68     Definition of the angles for oriented cylinders. 
    69  
    70 .. figure:: img/image062.JPG 
    71  
    72     Examples of the angles for oriented pp against the detector plane. 
    73  
    74 NB: The 2nd virial coefficient of the cylinder is calculated based on the 
    75 radius and length values, and used as the effective radius for $S(Q)$ 
    76 when $P(Q) \cdot S(Q)$ is applied. 
    77  
    78 The returned value is scaled to units of |cm^-1| and the parameters of 
    79 the CylinderModel are the following: 
    80  
    81 %(parameters)s 
    82  
    83 The output of the 1D scattering intensity function for randomly oriented 
    84 cylinders is then given by 
    85  
    86 .. math:: 
    87  
    88     P(q) = \frac{\text{scale}}{V} 
    89         \int_0^{\pi/2} f^2(q,\alpha) \sin \alpha d\alpha + \text{background} 
    90  
    91 The *theta* and *phi* parameters are not used for the 1D output. Our 
    92 implementation of the scattering kernel and the 1D scattering intensity 
    93 use the c-library from NIST. 
    94  
    95 Validation of the CylinderModel 
    96 ------------------------------- 
    97  
    98 Validation of our code was done by comparing the output of the 1D model 
    99 to the output of the software provided by the NIST (Kline, 2006). 
    100 Figure :num:`figure #CylinderModel-compare` shows a comparison of 
    101 the 1D output of our model and the output of the NIST software. 
    102  
    103 .. _CylinderModel-compare: 
    104  
    105 .. figure:: img/image065.JPG 
    106  
    107     Comparison of the SasView scattering intensity for a cylinder with the 
    108     output of the NIST SANS analysis software. 
    109     The parameters were set to: *Scale* = 1.0, *Radius* = 20 |Ang|, 
    110     *Length* = 400 |Ang|, *Contrast* = 3e-6 |Ang^-2|, and 
    111     *Background* = 0.01 |cm^-1|. 
    112  
    113 In general, averaging over a distribution of orientations is done by 
    114 evaluating the following 
    115  
    116 .. math:: 
    117  
    118     P(q) = \int_0^{\pi/2} d\phi 
    119         \int_0^\pi p(\theta, \phi) P_0(q,\alpha) \sin \theta d\theta 
    120  
    121  
    122 where $p(\theta,\phi)$ is the probability distribution for the orientation 
    123 and $P_0(q,\alpha)$ is the scattering intensity for the fully oriented 
    124 system. Since we have no other software to compare the implementation of 
    125 the intensity for fully oriented cylinders, we can compare the result of 
    126 averaging our 2D output using a uniform distribution $p(\theta, \phi) = 1.0$. 
    127 Figure :num:`figure #CylinderModel-crosscheck` shows the result of 
    128 such a cross-check. 
    129  
    130 .. _CylinderModel-crosscheck: 
    131  
    132 .. figure:: img/image066.JPG 
    133  
    134     Comparison of the intensity for uniformly distributed cylinders 
    135     calculated from our 2D model and the intensity from the NIST SANS 
    136     analysis software. 
    137     The parameters used were: *Scale* = 1.0, *Radius* = 20 |Ang|, 
    138     *Length* = 400 |Ang|, *Contrast* = 3e-6 |Ang^-2|, and 
    139     *Background* = 0.0 |cm^-1|. 
    140  
    141 DOCUMENTATION END 
    142 */ 
    1431real form_volume(real radius, real length); 
    1442real Iq(real q, real sld, real solvent_sld, real radius, real length); 
    1453real Iqxy(real qx, real qy, real sld, real solvent_sld, real radius, real length, real theta, real phi); 
    146  
    1474 
    1485real form_volume(real radius, real length) 
    … …  
    17330    return REAL(1.0e8) * form * s * s; 
    17431} 
    175  
    17632 
    17733real Iqxy(real qx, real qy, 

  • sasmodels/sasview_model.py

    r1780d59 r32c160a  
    44import numpy as np 
    55 
    6 def make_class(name, class_name=None, dtype='single'): 
     6def make_class(kernel_module, dtype='single'): 
    77    from .core import opencl_model 
    8     if class_name is None: 
    9         class_name = "".join(part.capitalize() for part in name.split('_')+['model']) 
    10     model =  opencl_model(name, dtype=dtype) 
    11     class ConstructedModel(SasviewModel): 
    12         kernel = model 
    13         def __init__(self, multfactor=1): 
    14             SasviewModel.__init__(self, self.kernel) 
    15     ConstructedModel.__name__ = class_name 
     8    model =  opencl_model(kernel_module, dtype=dtype) 
     9    def __init__(self, multfactor=1): 
     10        SasviewModel.__init__(self, self.kernel) 
     11    attrs = dict(__init__=__init__, kernel=model) 
     12    ConstructedModel = type(model.info['name'],  (SasviewModel,), attrs) 
     13    #class ConstructedModel(SasviewModel): 
     14    #    kernel = model 
     15    #    def __init__(self, multfactor=1): 
     16    #        SasviewModel.__init__(self, self.kernel) 
     17    #ConstructedModel.__name__ = model.info['name'] 
    1618    return ConstructedModel 
    1719 
    … …  
    267269            values, weights = self._dispersion_mesh(vol_pars) 
    268270            fv = ER(*values) 
     271            #print values[0].shape, weights.shape, fv.shape 
    269272            return np.sum(weights*fv) / np.sum(weights) 
    270273 
    … …  
    275278        :return: the value of the volf ratio 
    276279        """ 
    277         VR = self._model.info.get('ER', None) 
     280        VR = self._model.info.get('VR', None) 
    278281        if VR is None: 
    279282            return 1.0 
    … …  
    323326        values = [v.flatten() for v in np.meshgrid(*values)] 
    324327        weights = np.vstack([v.flatten() for v in np.meshgrid(*weights)]) 
    325         weights = np.prod(weights, axis=1) 
     328        weights = np.prod(weights, axis=0) 
    326329        return values, weights 
    327330 
Note: See TracChangeset for help on using the changeset viewer.