Changes in / [d5ac45f:31d22de] in sasmodels


Ignore:
Location:
sasmodels
Files:
2 added
1 edited

Legend:

Unmodified
Added
Removed
  • sasmodels/kernelpy.py

    rd5ac45f rd5ac45f  
    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.