Changeset b26d4c8 in sasmodels


Ignore:
Timestamp:
Mar 20, 2016 4:46:08 PM (8 years ago)
Author:
wojciech
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:
025c82d
Parents:
119bd3d (diff), cf52f9c (diff)
Note: this is a merge changeset, the changes displayed below correspond to the merge itself.
Use the (diff) links above to see all the changes relative to each parent.
Message:

Merged with remote

Files:
2 added
4 edited

Legend:

Unmodified
Added
Removed

  • doc/developer/index.rst

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

  • sasmodels/core.py

    r35647ab rcf52f9c  
    216216    with an error of about 1%, which is usually less than the measurement 
    217217    uncertainty. 
    218     """ 
    219     print pars 
     218 
     219    *mono* is True if polydispersity should be set to none on all parameters. 
     220    """ 
    220221    fixed_pars = [pars.get(name, kernel.info['defaults'][name]) 
    221222                  for name in kernel.fixed_pars] 

  • sasmodels/kernel_iq.c

    r119bd3d r4b2972c  
    1212*/ 
    1313 
    14 /* 
    15 The environment needs to provide the following #defines: 
    16  
    17    USE_OPENCL is defined if running in opencl 
    18    KERNEL declares a function to be available externally 
    19    KERNEL_NAME is the name of the function being declared 
    20    NPARS is the number of parameters in the kernel 
    21    PARAMETER_DECL is the declaration of the parameters to the kernel. 
    22  
    23        Cylinder: 
    24  
    25            #define PARAMETER_DECL \ 
    26            double length; \ 
    27            double radius; \ 
    28            double sld; \ 
    29            double sld_solvent 
    30  
    31        Note: scale and background are not included 
    32  
    33        Multi-shell cylinder (10 shell max): 
    34  
    35            #define PARAMETER_DECL \ 
    36            double num_shells; \ 
    37            double length; \ 
    38            double radius[10]; \ 
    39            double sld[10]; \ 
    40            double sld_solvent 
    41  
    42    CALL_IQ(q, nq, i, pars) is the declaration of a call to the kernel. 
    43  
    44        Cylinder: 
    45  
    46            #define CALL_IQ(q, nq, i, var) \ 
    47            Iq(q[i], \ 
    48            var.length, \ 
    49            var.radius, \ 
    50            var.sld, \ 
    51            var.sld_solvent) 
    52  
    53        Multi-shell cylinder: 
    54            #define CALL_IQ(q, nq, i, var) \ 
    55            Iq(q[i], \ 
    56            var.num_shells, \ 
    57            var.length, \ 
    58            var.radius, \ 
    59            var.sld, \ 
    60            var.sld_solvent) 
    61  
    62    CALL_VOLUME(var) is similar, but for calling the form volume. 
    63  
    64    INVALID is a test for model parameters in the correct range 
    65  
    66        Cylinder: 
    67  
    68            #define INVALID(var) 0 
    69  
    70        BarBell: 
    71  
    72            #define INVALID(var) (var.bell_radius > var.radius) 
    73  
    74        Model with complicated constraints: 
    75  
    76            inline bool constrained(p1, p2, p3) { return expression; } 
    77            #define INVALID(var) constrained(var.p1, var.p2, var.p3) 
    78  
    79    IQ_FUNC could be Iq or Iqxy 
    80    IQ_PARS could be q[i] or q[2*i],q[2*i+1] 
    81  
    82 Our design supports a limited number of polydispersity loops, wherein 
    83 we need to cycle through the values of the polydispersity, calculate 
    84 the I(q, p) for each combination of parameters, and perform a normalized 
    85 weighted sum across all the weights.  Parameters may be passed to the 
    86 underlying calculation engine as scalars or vectors, but the polydispersity 
    87 calculator treats the parameter set as one long vector. 
    88  
    89 Let's assume we have 6 parameters in the model, with two polydisperse:: 
    90  
    91     0: scale        {scl = constant} 
    92     1: background   {bkg = constant} 
    93     5: length       {l = vector of 30pts} 
    94     4: radius       {r = vector of 10pts} 
    95     3: sld          {s = constant/(radius**2*length)} 
    96     2: sld_solvent  {s2 = constant} 
    97  
    98 This generates the following call to the kernel (where x stands for an 
    99 arbitrary value that is not used by the kernel evaluator): 
    100  
    101     NPARS = 4  // scale and background are in all models 
    102     problem { 
    103         pd_par = {5, 4, x, x}         // parameters *radius* and *length* vary 
    104         pd_length = {30, 10, 0, 0}    // *length* has more, so it is first 
    105         pd_offset = {10, 0, x, x}     // *length* starts at index 10 in weights 
    106         pd_stride = {1, 30, 300, 300} // cumulative product of pd length 
    107         pd_isvol = {1, 1, x, x}       // true if weight is a volume weight 
    108         par_offset = {2, 3, 303, 313}  // parameter offsets 
    109         par_coord = {0, 3, 2, 1} // bitmap of parameter dependencies 
    110         fast_coord_index = {5, 3, x, x} 
    111         fast_coord_count = 2  // two parameters vary with *length* distribution 
    112         theta_var = -1   // no spherical correction 
    113         fast_theta = 0   // spherical correction angle is not pd 1 
    114     } 
    115  
    116     weight = { l0, .., l29, r0, .., r9} 
    117     pars = { scl, bkg, l0, ..., l29, r0, r1, ..., r9, 
    118              s[l0,r0], ... s[l0,r9], s[l1,r0], ... s[l29,r9] , s2} 
    119  
    120     nq = 130 
    121     q = { q0, q1, ..., q130, x, x }  # pad to 8 element boundary 
    122     result = {r1, ..., r130, norm, vol, vol_norm, x, x, x, x, x, x, x} 
    123  
    124  
    125 The polydisperse parameters are stored in as an array of parameter 
    126 indices, one for each polydisperse parameter, stored in pd_par[n]. 
    127 Non-polydisperse parameters do not appear in this array. Each polydisperse 
    128 parameter has a weight vector whose length is stored in pd_length[n], 
    129 The weights are stored in a contiguous vector of weights for all 
    130 parameters, with the starting position for the each parameter stored 
    131 in pd_offset[n].  The values corresponding to the weights are stored 
    132 together in a separate weights[] vector, with offset stored in 
    133 par_offset[pd_par[n]]. Polydisperse parameters should be stored in 
    134 decreasing order of length for highest efficiency. 
    135  
    136 We limit the number of polydisperse dimensions to MAX_PD (currently 4). 
    137 This cuts the size of the structure in half compared to allowing a 
    138 separate polydispersity for each parameter.  This will help a little 
    139 bit for models with large numbers of parameters, such as the onion model. 
    140  
    141 Parameters may be coordinated.  That is, we may have the value of one 
    142 parameter depend on a set of other parameters, some of which may be 
    143 polydisperse.  For example, if sld is inversely proportional to the 
    144 volume of a cylinder, and the length and radius are independently 
    145 polydisperse, then for each combination of length and radius we need a 
    146 separate value for the sld.  The caller must provide a coordination table 
    147 for each parameter containing the value for each parameter given the 
    148 value of the polydisperse parameters v1, v2, etc.  The tables for each 
    149 parameter are arranged contiguously in a vector, with offset[k] giving the 
    150 starting location of parameter k in the vector.  Each parameter defines 
    151 coord[k] as a bit mask indicating which polydispersity parameters the 
    152 parameter depends upon. Usually this is zero, indicating that the parameter 
    153 is independent, but for the cylinder example given, the bits for the 
    154 radius and length polydispersity parameters would both be set, the result 
    155 being a (#radius x #length) table, or maybe a (#length x #radius) table 
    156 if length comes first in the polydispersity table. 
    157  
    158 NB: If we can guarantee that a compiler and OpenCL driver are available, 
    159 we could instead create the coordination function on the fly for each 
    160 parameter, saving memory and transfer time, but requiring a C compiler 
    161 as part of the environment. 
    162  
    163 In ordering the polydisperse parameters by decreasing length we can 
    164 iterate over the longest dispersion weight vector first.  All parameters 
    165 coordinated with this weight vector (the 'fast' parameters), can be 
    166 updated with a simple increment to the next position in the parameter 
    167 value table.  The indices of these parameters is stored in fast_coord_index[], 
    168 with fast_coord_count being the number of fast parameters.  A total 
    169 of NPARS slots is allocated to allow for the case that all parameters 
    170 are coordinated with the fast index, though this will likely be mostly 
    171 empty.  When the fast increment count reaches the end of the weight 
    172 vector, then the index of the second polydisperse parameter must be 
    173 incremented, and all of its coordinated parameters updated.  Because this 
    174 operation is not in the inner loop, a slower algorithm can be used. 
    175  
    176 If there is no polydispersity we pretend that it is polydisperisty with one 
    177 parameter, pd_start=0 and pd_stop=1.  We may or may not short circuit the 
    178 calculation in this case, depending on how much time it saves. 
    179  
    180 The problem details structure can be allocated and sent in as an integer 
    181 array using the read-only flag.  This allows us to copy it once per fit 
    182 along with the weights vector, since features such as the number of 
    183 polydisperity elements per pd parameter or the coordinated won't change 
    184 between function evaluations.  A new parameter vector is sent for 
    185 each I(q) evaluation. 
    186  
    187 To protect against expensive evaluations taking all the GPU resource 
    188 on large fits, the entire polydispersity will not be computed at once. 
    189 Instead, a start and stop location will be sent, indicating where in the 
    190 polydispersity loop the calculation should start and where it should 
    191 stop.  We can do this for arbitrary start/stop points since we have 
    192 unwound the nested loop.  Instead, we use the same technique as array 
    193 index translation, using div and mod to figure out the i,j,k,... 
    194 indices in the virtual nested loop. 
    195  
    196 The results array will be initialized to zero for polydispersity loop 
    197 entry zero, and preserved between calls to [start, stop] so that the 
    198 results accumulate by the time the loop has completed.  Background and 
    199 scale will be applied when the loop reaches the end.  This does require 
    200 that the results array be allocated read-write, which is less efficient 
    201 for the GPU, but it makes the calling sequence much more manageable. 
    202  
    203 Scale and background cannot be coordinated with other polydisperse parameters 
    204  
    205 Oriented objects in 2-D need a spherical correction on the angular variation 
    206 in order to preserve the 'surface area' of the weight distribution. 
    207  
    208 Cutoff specifies the integration area by discarding regions in polydisperisty 
    209 hypercubue that has no parameters defined 
    210 */ 
    21114 
    21215#define MAX_PD 4  // MAX_PD is the max number of polydisperse parameters 
    … …  
    23134} ParameterBlock; 
    23235 
    233 #define KERNEL_NAME test_Iq 
    234 #define FULL_KERNEL_NAME test_Iq 
    235 #define IQ_FUNC Iq 
    236  
    237  
     36#define FULL_KERNEL_NAME KERNEL_NAME ## _ ## IQ_FUNC 
     37KERNEL 
    23838void FULL_KERNEL_NAME( 
    23939    int nq,                 // number of q values 
    … …  
    404204  } 
    405205} 
    406 } 
    407 } 

  • sasmodels/kernelpy.py

    r352964b r119bd3d  
    128128 
    129129    def __call__(self, fixed, pd, cutoff=1e-5): 
    130         print("fixed",fixed) 
    131         print("pd", pd) 
     130        #print("fixed",fixed) 
     131        #print("pd", pd) 
    132132        args = self.args[:]  # grab a copy of the args 
    133133        form, form_volume = self.kernel, self.info['form_volume'] 
    … …  
    187187    ################################################################ 
    188188 
    189     #TODO: Wojtek's notes 
    190     #TODO: The goal is to restructure polydispersity loop 
    191     #so it allows fitting arbitrary polydispersity parameters 
    192     #and it accepts cases like coupled parameters 
    193189    weight = np.empty(len(pd), 'd') 
    194190    if weight.size > 0: 
    … …  
    264260    result = scale * ret / norm + background 
    265261    return result 
    266 ======= 
    267 """ 
    268 Python driver for python kernels 
    269  
    270 Calls the kernel with a vector of $q$ values for a single parameter set. 
    271 Polydispersity is supported by looping over different parameter sets and 
    272 summing the results.  The interface to :class:`PyModel` matches those for 
    273 :class:`kernelcl.GpuModel` and :class:`kerneldll.DllModel`. 
    274 """ 
    275 import numpy as np 
    276 from numpy import pi, cos 
    277  
    278 from .generate import F64 
    279  
    280 class PyModel(object): 
    281     """ 
    282     Wrapper for pure python models. 
    283     """ 
    284     def __init__(self, model_info): 
    285         self.info = model_info 
    286  
    287     def __call__(self, q_vectors): 
    288         q_input = PyInput(q_vectors, dtype=F64) 
    289         kernel = self.info['Iqxy'] if q_input.is_2d else self.info['Iq'] 
    290         return PyKernel(kernel, self.info, q_input) 
    291  
    292     def release(self): 
    293         """ 
    294         Free resources associated with the model. 
    295         """ 
    296         pass 
    297  
    298 class PyInput(object): 
    299     """ 
    300     Make q data available to the gpu. 
    301  
    302     *q_vectors* is a list of q vectors, which will be *[q]* for 1-D data, 
    303     and *[qx, qy]* for 2-D data.  Internally, the vectors will be reallocated 
    304     to get the best performance on OpenCL, which may involve shifting and 
    305     stretching the array to better match the memory architecture.  Additional 
    306     points will be evaluated with *q=1e-3*. 
    307  
    308     *dtype* is the data type for the q vectors. The data type should be 
    309     set to match that of the kernel, which is an attribute of 
    310     :class:`GpuProgram`.  Note that not all kernels support double 
    311     precision, so even if the program was created for double precision, 
    312     the *GpuProgram.dtype* may be single precision. 
    313  
    314     Call :meth:`release` when complete.  Even if not called directly, the 
    315     buffer will be released when the data object is freed. 
    316     """ 
    317     def __init__(self, q_vectors, dtype): 
    318         self.nq = q_vectors[0].size 
    319         self.dtype = dtype 
    320         self.is_2d = (len(q_vectors) == 2) 
    321         self.q_vectors = [np.ascontiguousarray(q, self.dtype) for q in q_vectors] 
    322         self.q_pointers = [q.ctypes.data for q in self.q_vectors] 
    323  
    324     def release(self): 
    325         """ 
    326         Free resources associated with the model inputs. 
    327         """ 
    328         self.q_vectors = [] 
    329  
    330 class PyKernel(object): 
    331     """ 
    332     Callable SAS kernel. 
    333  
    334     *kernel* is the DllKernel object to call. 
    335  
    336     *model_info* is the module information 
    337  
    338     *q_input* is the DllInput q vectors at which the kernel should be 
    339     evaluated. 
    340  
    341     The resulting call method takes the *pars*, a list of values for 
    342     the fixed parameters to the kernel, and *pd_pars*, a list of (value,weight) 
    343     vectors for the polydisperse parameters.  *cutoff* determines the 
    344     integration limits: any points with combined weight less than *cutoff* 
    345     will not be calculated. 
    346  
    347     Call :meth:`release` when done with the kernel instance. 
    348     """ 
    349     def __init__(self, kernel, model_info, q_input): 
    350         self.info = model_info 
    351         self.q_input = q_input 
    352         self.res = np.empty(q_input.nq, q_input.dtype) 
    353         dim = '2d' if q_input.is_2d else '1d' 
    354         # Loop over q unless user promises that the kernel is vectorized by 
    355         # taggining it with vectorized=True 
    356         if not getattr(kernel, 'vectorized', False): 
    357             if dim == '2d': 
    358                 def vector_kernel(qx, qy, *args): 
    359                     """ 
    360                     Vectorized 2D kernel. 
    361                     """ 
    362                     return np.array([kernel(qxi, qyi, *args) 
    363                                      for qxi, qyi in zip(qx, qy)]) 
    364             else: 
    365                 def vector_kernel(q, *args): 
    366                     """ 
    367                     Vectorized 1D kernel. 
    368                     """ 
    369                     return np.array([kernel(qi, *args) 
    370                                      for qi in q]) 
    371             self.kernel = vector_kernel 
    372         else: 
    373             self.kernel = kernel 
    374         fixed_pars = model_info['partype']['fixed-' + dim] 
    375         pd_pars = model_info['partype']['pd-' + dim] 
    376         vol_pars = model_info['partype']['volume'] 
    377  
    378         # First two fixed pars are scale and background 
    379         pars = [p.name for p in model_info['parameters'][2:]] 
    380         offset = len(self.q_input.q_vectors) 
    381         self.args = self.q_input.q_vectors + [None] * len(pars) 
    382         self.fixed_index = np.array([pars.index(p) + offset 
    383                                      for p in fixed_pars[2:]]) 
    384         self.pd_index = np.array([pars.index(p) + offset 
    385                                   for p in pd_pars]) 
    386         self.vol_index = np.array([pars.index(p) + offset 
    387                                    for p in vol_pars]) 
    388         try: self.theta_index = pars.index('theta') + offset 
    389         except ValueError: self.theta_index = -1 
    390  
    391         # Caller needs fixed_pars and pd_pars 
    392         self.fixed_pars = fixed_pars 
    393         self.pd_pars = pd_pars 
    394  
    395     def __call__(self, fixed, pd, cutoff=1e-5): 
    396         #print("fixed",fixed) 
    397         #print("pd", pd) 
    398         args = self.args[:]  # grab a copy of the args 
    399         form, form_volume = self.kernel, self.info['form_volume'] 
    400         # First two fixed 
    401         scale, background = fixed[:2] 
    402         for index, value in zip(self.fixed_index, fixed[2:]): 
    403             args[index] = float(value) 
    404         res = _loops(form, form_volume, cutoff, scale, background, args, 
    405                      pd, self.pd_index, self.vol_index, self.theta_index) 
    406  
    407         return res 
    408  
    409     def release(self): 
    410         """ 
    411         Free resources associated with the kernel. 
    412         """ 
    413         self.q_input = None 
    414  
    415 def _loops(form, form_volume, cutoff, scale, background, 
    416            args, pd, pd_index, vol_index, theta_index): 
    417     """ 
    418     *form* is the name of the form function, which should be vectorized over 
    419     q, but otherwise have an interface like the opencl kernels, with the 
    420     q parameters first followed by the individual parameters in the order 
    421     given in model.parameters (see :mod:`sasmodels.generate`). 
    422  
    423     *form_volume* calculates the volume of the shape.  *vol_index* gives 
    424     the list of volume parameters 
    425  
    426     *cutoff* ignores the corners of the dispersion hypercube 
    427  
    428     *scale*, *background* multiplies the resulting form and adds an offset 
    429  
    430     *args* is the prepopulated set of arguments to the form function, starting 
    431     with the q vectors, and including slots for all the parameters.  The 
    432     values for the parameters will be substituted with values from the 
    433     dispersion functions. 
    434  
    435     *pd* is the list of dispersion parameters 
    436  
    437     *pd_index* are the indices of the dispersion parameters in *args* 
    438  
    439     *vol_index* are the indices of the volume parameters in *args* 
    440  
    441     *theta_index* is the index of the theta parameter for the sasview 
    442     spherical correction, or -1 if there is no angular dispersion 
    443     """ 
    444  
    445     ################################################################ 
    446     #                                                              # 
    447     #   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!   # 
    448     #   !!                                                    !!   # 
    449     #   !!  KEEP THIS CODE CONSISTENT WITH KERNEL_TEMPLATE.C  !!   # 
    450     #   !!                                                    !!   # 
    451     #   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!   # 
    452     #                                                              # 
    453     ################################################################ 
    454  
    455     weight = np.empty(len(pd), 'd') 
    456     if weight.size > 0: 
    457         # weight vector, to be populated by polydispersity loops 
    458  
    459         # identify which pd parameters are volume parameters 
    460         vol_weight_index = np.array([(index in vol_index) for index in pd_index]) 
    461  
    462         # Sort parameters in decreasing order of pd length 
    463         Npd = np.array([len(pdi[0]) for pdi in pd], 'i') 
    464         order = np.argsort(Npd)[::-1] 
    465         stride = np.cumprod(Npd[order]) 
    466         pd = [pd[index] for index in order] 
    467         pd_index = pd_index[order] 
    468         vol_weight_index = vol_weight_index[order] 
    469  
    470         fast_value = pd[0][0] 
    471         fast_weight = pd[0][1] 
    472     else: 
    473         stride = np.array([1]) 
    474         vol_weight_index = slice(None, None) 
    475         # keep lint happy 
    476         fast_value = [None] 
    477         fast_weight = [None] 
    478  
    479     ret = np.zeros_like(args[0]) 
    480     norm = np.zeros_like(ret) 
    481     vol = np.zeros_like(ret) 
    482     vol_norm = np.zeros_like(ret) 
    483     for k in range(stride[-1]): 
    484         # update polydispersity parameter values 
    485         fast_index = k % stride[0] 
    486         if fast_index == 0:  # bottom loop complete ... check all other loops 
    487             if weight.size > 0: 
    488                 for i, index, in enumerate(k % stride): 
    489                     args[pd_index[i]] = pd[i][0][index] 
    490                     weight[i] = pd[i][1][index] 
    491         else: 
    492             args[pd_index[0]] = fast_value[fast_index] 
    493             weight[0] = fast_weight[fast_index] 
    494  
    495         # Computes the weight, and if it is not sufficient then ignore this 
    496         # parameter set. 
    497         # Note: could precompute w1*...*wn so we only need to multiply by w0 
    498         w = np.prod(weight) 
    499         if w > cutoff: 
    500             # Note: can precompute spherical correction if theta_index is not 
    501             # the fast index. Correction factor for spherical integration 
    502             #spherical_correction = abs(cos(pi*args[phi_index])) if phi_index>=0 else 1.0 
    503             spherical_correction = (abs(cos(pi * args[theta_index])) * pi / 2 
    504                                     if theta_index >= 0 else 1.0) 
    505             #spherical_correction = 1.0 
    506  
    507             # Call the scattering function and adds it to the total. 
    508             I = form(*args) 
    509             if np.isnan(I).any(): continue 
    510             ret += w * I * spherical_correction 
    511             norm += w 
    512  
    513             # Volume normalization. 
    514             # If there are "volume" polydispersity parameters, then these 
    515             # will be used to call the form_volume function from the user 
    516             # supplied kernel, and accumulate a normalized weight. 
    517             # Note: can precompute volume norm if fast index is not a volume 
    518             if form_volume: 
    519                 vol_args = [args[index] for index in vol_index] 
    520                 vol_weight = np.prod(weight[vol_weight_index]) 
    521                 vol += vol_weight * form_volume(*vol_args) 
    522                 vol_norm += vol_weight 
    523  
    524     positive = (vol * vol_norm != 0.0) 
    525     ret[positive] *= vol_norm[positive] / vol[positive] 
    526     result = scale * ret / norm + background 
    527     return result 
Note: See TracChangeset for help on using the changeset viewer.