Changeset 3707eee in sasmodels


Ignore:
Timestamp:
Apr 7, 2016 12:26:04 PM (9 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:
6d6508e
Parents:
a8a7f08 (diff), cb3a824 (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:

Merge branch 'polydisp' of github.com:sasview/sasmodels into polydisp

Location:
sasmodels
Files:
11 edited

Legend:

Unmodified
Added
Removed

  • sasmodels/compare.py

    r21b116f r8bd7b77  
    209209    Choose a parameter range based on parameter name and initial value. 
    210210    """ 
     211    # process the polydispersity options 
    211212    if p.endswith('_pd_n'): 
    212213        return [0, 100] 
    … …  
    221222        else: 
    222223            return [-180, 180] 
     224    elif p.endswith('_pd'): 
     225        return [0, 1] 
    223226    elif 'sld' in p: 
    224227        return [-0.5, 10] 
    225     elif p.endswith('_pd'): 
    226         return [0, 1] 
    227228    elif p == 'background': 
    228229        return [0, 10] 
    229230    elif p == 'scale': 
    230231        return [0, 1e3] 
    231     elif p == 'case_num': 
    232         # RPA hack 
    233         return [0, 10] 
    234232    elif v < 0: 
    235         # Kxy parameters in rpa model can be negative 
    236233        return [2*v, -2*v] 
    237234    else: 
    … …  
    239236 
    240237 
    241 def _randomize_one(p, v): 
     238def _randomize_one(model_info, p, v): 
    242239    """ 
    243240    Randomize a single parameter. 
    … …  
    245242    if any(p.endswith(s) for s in ('_pd_n', '_pd_nsigma', '_pd_type')): 
    246243        return v 
     244 
     245    # Find the parameter definition 
     246    for par in model_info['parameters'].call_parameters: 
     247        if par.name == p: 
     248            break 
    247249    else: 
    248         return np.random.uniform(*parameter_range(p, v)) 
    249  
    250  
    251 def randomize_pars(pars, seed=None): 
     250        raise ValueError("unknown parameter %r"%p) 
     251 
     252    if len(par.limits) > 2:  # choice list 
     253        return np.random.randint(len(par.limits)) 
     254 
     255    limits = par.limits 
     256    if not np.isfinite(limits).all(): 
     257        low, high = parameter_range(p, v) 
     258        limits = (max(limits[0], low), min(limits[1], high)) 
     259 
     260    return np.random.uniform(*limits) 
     261 
     262 
     263def randomize_pars(model_info, pars, seed=None): 
    252264    """ 
    253265    Generate random values for all of the parameters. 
    … …  
    260272    with push_seed(seed): 
    261273        # Note: the sort guarantees order `of calls to random number generator 
    262         pars = dict((p, _randomize_one(p, v)) 
     274        pars = dict((p, _randomize_one(model_info, p, v)) 
    263275                    for p, v in sorted(pars.items())) 
    264276    return pars 
    … …  
    293305        # Make sure phi sums to 1.0 
    294306        if pars['case_num'] < 2: 
    295             pars['Phia'] = 0. 
    296             pars['Phib'] = 0. 
     307            pars['Phi1'] = 0. 
     308            pars['Phi2'] = 0. 
    297309        elif pars['case_num'] < 5: 
    298             pars['Phia'] = 0. 
    299         total = sum(pars['Phi'+c] for c in 'abcd') 
    300         for c in 'abcd': 
     310            pars['Phi1'] = 0. 
     311        total = sum(pars['Phi'+c] for c in '1234') 
     312        for c in '1234': 
    301313            pars['Phi'+c] /= total 
    302314 
    … …  
    805817    #pars.update(set_pars)  # set value before random to control range 
    806818    if opts['seed'] > -1: 
    807         pars = randomize_pars(pars, seed=opts['seed']) 
     819        pars = randomize_pars(model_info, pars, seed=opts['seed']) 
    808820        print("Randomize using -random=%i"%opts['seed']) 
    809821    if opts['mono']: 

  • sasmodels/compare_many.py

    ree8f734 r8bd7b77  
    108108 
    109109    if is_2d: 
    110         if not model_info['has_2d']: 
     110        if not model_info['parameters'].has_2d: 
    111111            print(',"1-D only"') 
    112112            return 

  • sasmodels/convert.json

    rd2bb604 r6e7ff6d  
    602602    "RPAModel",  
    603603    { 
     604      "N1": "Na", "N2": "Nb", "N3": "Nc", "N4": "Nd", 
     605      "L1": "La", "L2": "Lb", "L3": "Lc", "L4": "Ld", 
     606      "v1": "va", "v2": "vb", "v3": "vc", "v4": "vd", 
     607      "b1": "ba", "b2": "bb", "b3": "bc", "b4": "bd", 
     608      "Phi1": "Phia", "Phi2": "Phib", "Phi3": "Phic", "Phi4": "Phid", 
    604609      "case_num": "lcase_n" 
    605610    } 
    … …  
    620625    } 
    621626  ],  
     627  "_spherepy": [ 
     628    "SphereModel", 
     629    { 
     630      "sld": "sldSph", 
     631      "radius": "radius", 
     632      "sld_solvent": "sldSolv" 
     633    } 
     634  ], 
    622635  "spherical_sld": [ 
    623636    "SphereSLDModel",  

  • sasmodels/convert.py

    rf247314 r8bd7b77  
    206206        elif name == 'rpa': 
    207207            # convert scattering lengths from femtometers to centimeters 
    208             for p in "La", "Lb", "Lc", "Ld": 
     208            for p in "L1", "L2", "L3", "L4": 
    209209                if p in oldpars: oldpars[p] *= 1e-13 
    210210        elif name == 'core_shell_parallelepiped': 

  • sasmodels/generate.py

    r4937980 r6e7ff6d  
    125125An *model_info* dictionary is constructed from the kernel meta data and 
    126126returned to the caller. 
    127  
    128 The model evaluator, function call sequence consists of q inputs and the return vector, 
    129 followed by the loop value/weight vector, followed by the values for 
    130 the non-polydisperse parameters, followed by the lengths of the 
    131 polydispersity loops.  To construct the call for 1D models, the 
    132 categories *fixed-1d* and *pd-1d* list the names of the parameters 
    133 of the non-polydisperse and the polydisperse parameters respectively. 
    134 Similarly, *fixed-2d* and *pd-2d* provide parameter names for 2D models. 
    135 The *pd-rel* category is a set of those parameters which give 
    136 polydispersitiy as a portion of the value (so a 10% length dispersity 
    137 would use a polydispersity value of 0.1) rather than absolute 
    138 dispersity such as an angle plus or minus 15 degrees. 
    139  
    140 The *volume* category lists the volume parameters in order for calls 
    141 to volume within the kernel (used for volume normalization) and for 
    142 calls to ER and VR for effective radius and volume ratio respectively. 
    143  
    144 The *orientation* and *magnetic* categories list the orientation and 
    145 magnetic parameters.  These are used by the sasview interface.  The 
    146 blank category is for parameters such as scale which don't have any 
    147 other marking. 
    148127 
    149128The doc string at the start of the kernel module will be used to 
    … …  
    552531    @property 
    553532    def ctypes(self): return self.details.ctypes 
     533 
    554534    @property 
    555535    def pd_par(self): return self._pd_par 
     536 
    556537    @property 
    557538    def pd_length(self): return self._pd_length 
     539 
    558540    @property 
    559541    def pd_offset(self): return self._pd_offset 
     542 
    560543    @property 
    561544    def pd_stride(self): return self._pd_stride 
     545 
    562546    @property 
    563547    def pd_coord(self): return self._pd_coord 
     548 
    564549    @property 
    565550    def par_coord(self): return self._par_coord 
     551 
    566552    @property 
    567553    def par_offset(self): return self._par_offset 
     554 
     555    @property 
     556    def num_active(self): return self.details[-4] 
     557    @num_active.setter 
     558    def num_active(self, v): self.details[-4] = v 
     559 
     560    @property 
     561    def total_pd(self): return self.details[-3] 
     562    @total_pd.setter 
     563    def total_pd(self, v): self.details[-3] = v 
     564 
    568565    @property 
    569566    def num_coord(self): return self.details[-2] 
    570567    @num_coord.setter 
    571568    def num_coord(self, v): self.details[-2] = v 
    572     @property 
    573     def total_pd(self): return self.details[-3] 
    574     @total_pd.setter 
    575     def total_pd(self, v): self.details[-3] = v 
    576     @property 
    577     def num_active(self): return self.details[-4] 
    578     @num_active.setter 
    579     def num_active(self, v): self.details[-4] = v 
     569 
     570    @property 
     571    def theta_par(self): return self.details[-1] 
    580572 
    581573    def show(self): 
    … …  
    638630 
    639631 
     632def create_vector_Iq(model_info): 
     633    """ 
     634    Define Iq as a vector function if it exists. 
     635    """ 
     636    Iq = model_info['Iq'] 
     637    if callable(Iq) and not getattr(Iq, 'vectorized', False): 
     638        def vector_Iq(q, *args): 
     639            """ 
     640            Vectorized 1D kernel. 
     641            """ 
     642            return np.array([Iq(qi, *args) for qi in q]) 
     643        model_info['Iq'] = vector_Iq 
     644 
     645def create_vector_Iqxy(model_info): 
     646    """ 
     647    Define Iqxy as a vector function if it exists, or default it from Iq(). 
     648    """ 
     649    Iq, Iqxy = model_info['Iq'], model_info['Iqxy'] 
     650    if callable(Iqxy) and not getattr(Iqxy, 'vectorized', False): 
     651        def vector_Iqxy(qx, qy, *args): 
     652            """ 
     653            Vectorized 2D kernel. 
     654            """ 
     655            return np.array([Iqxy(qxi, qyi, *args) for qxi, qyi in zip(qx, qy)]) 
     656        model_info['Iqxy'] = vector_Iqxy 
     657    elif callable(Iq): 
     658        # Iq is vectorized because create_vector_Iq was already called. 
     659        def default_Iqxy(qx, qy, *args): 
     660            """ 
     661            Default 2D kernel. 
     662            """ 
     663            return Iq(np.sqrt(qx**2 + qy**2), *args) 
     664        model_info['Iqxy'] = default_Iqxy 
     665 
    640666def create_default_functions(model_info): 
    641667    """ 
    … …  
    643669 
    644670    This only works for Iqxy when Iq is written in python. :func:`make_source` 
    645     performs a similar role for Iq written in C. 
    646     """ 
    647     if callable(model_info['Iq']) and model_info['Iqxy'] is None: 
    648         partable = model_info['parameters'] 
    649         if partable.has_2d: 
    650             raise ValueError("Iqxy model is missing") 
    651         Iq = model_info['Iq'] 
    652         def Iqxy(qx, qy, **kw): 
    653             return Iq(np.sqrt(qx**2 + qy**2), **kw) 
    654         model_info['Iqxy'] = Iqxy 
    655  
     671    performs a similar role for Iq written in C.  This also vectorizes 
     672    any functions that are not already marked as vectorized. 
     673    """ 
     674    create_vector_Iq(model_info) 
     675    create_vector_Iqxy(model_info)  # call create_vector_Iq() first 
    656676 
    657677def load_kernel_module(model_name): 

  • sasmodels/kernel_iq.c

    rcb3a824 r6e7ff6d  
    4242    const int32_t pd_start,     // where we are in the polydispersity loop 
    4343    const int32_t pd_stop,      // where we are stopping in the polydispersity loop 
    44     global const ProblemDetails *problem, 
     44    global const ProblemDetails *details, 
    4545    global const double *weights, 
    4646    global const double *values, 
    … …  
    6060  #endif 
    6161  for (int k=0; k < NPARS; k++) { 
    62     pvec[k] = values[problem->par_offset[k]]; 
     62    pvec[k] = values[details->par_offset[k]]; 
    6363  } 
    6464 
    … …  
    7878 
    7979  // Monodisperse computation 
    80   if (problem->num_active == 0) { 
     80  if (details->num_active == 0) { 
    8181    #ifdef INVALID 
    8282    if (INVALID(local_values)) { return; } 
    … …  
    116116  // Trigger the reset behaviour that happens at the end the fast loop 
    117117  // by setting the initial index >= weight vector length. 
    118   const int fast_length = problem->pd_length[0]; 
     118  const int fast_length = details->pd_length[0]; 
    119119  pd_index[0] = fast_length; 
    120120 
    … …  
    123123    // check if indices need to be updated 
    124124    if (pd_index[0] == fast_length) { 
    125       //printf("should be here with %d active\n", problem->num_active); 
     125      //printf("should be here with %d active\n", details->num_active); 
    126126 
    127127      // Compute position in polydispersity hypercube 
    128       for (int k=0; k < problem->num_active; k++) { 
    129         pd_index[k] = (loop_index/problem->pd_stride[k])%problem->pd_length[k]; 
     128      for (int k=0; k < details->num_active; k++) { 
     129        pd_index[k] = (loop_index/details->pd_stride[k])%details->pd_length[k]; 
    130130        //printf("pd_index[%d] = %d\n",k,pd_index[k]); 
    131131      } 
    … …  
    134134      partial_weight = 1.0; 
    135135      //printf("partial weight %d: ", loop_index); 
    136       for (int k=1; k < problem->num_active; k++) { 
    137         double wi = weights[problem->pd_offset[k] + pd_index[k]]; 
    138         //printf("pd[%d]=par[%d]=%g ", k, problem->pd_par[k], wi); 
     136      for (int k=1; k < details->num_active; k++) { 
     137        double wi = weights[details->pd_offset[k] + pd_index[k]]; 
     138        //printf("pd[%d]=par[%d]=%g ", k, details->pd_par[k], wi); 
    139139        partial_weight *= wi; 
    140140      } 
    … …  
    143143      // Update parameter offsets in weight vector 
    144144      //printf("slow %d: ", loop_index); 
    145       for (int k=0; k < problem->num_coord; k++) { 
    146         int par = problem->par_coord[k]; 
    147         int coord = problem->pd_coord[k]; 
    148         int this_offset = problem->par_offset[par]; 
     145      for (int k=0; k < details->num_coord; k++) { 
     146        int par = details->par_coord[k]; 
     147        int coord = details->pd_coord[k]; 
     148        int this_offset = details->par_offset[par]; 
    149149        int block_size = 1; 
    150150        for (int bit=0; coord != 0; bit++) { 
    151151          if (coord&1) { 
    152152              this_offset += block_size * pd_index[bit]; 
    153               block_size *= problem->pd_length[bit]; 
     153              block_size *= details->pd_length[bit]; 
    154154          } 
    155155          coord >>= 1; 
    … …  
    159159        //printf("par[%d]=v[%d]=%g \n", k, offset[k], pvec[k]); 
    160160        // if theta is not coordinated with fast index, precompute spherical correction 
    161         if (par == problem->theta_par && !(problem->par_coord[k]&1)) { 
    162           spherical_correction = fmax(fabs(cos(M_PI_180*pvec[problem->theta_par])), 1e-6); 
     161        if (par == details->theta_par && !(details->par_coord[k]&1)) { 
     162          spherical_correction = fmax(fabs(cos(M_PI_180*pvec[details->theta_par])), 1e-6); 
    163163        } 
    164164      } 
    … …  
    167167 
    168168    // Increment fast index 
    169     const double wi = weights[problem->pd_offset[0] + pd_index[0]++]; 
     169    const double wi = weights[details->pd_offset[0] + pd_index[0]++]; 
    170170    double weight = partial_weight*wi; 
    171171    //printf("fast %d: ", loop_index); 
    172     for (int k=0; k < problem->num_coord; k++) { 
    173       if (problem->pd_coord[k]&1) { 
    174         const int par = problem->par_coord[k]; 
     172    for (int k=0; k < details->num_coord; k++) { 
     173      if (details->pd_coord[k]&1) { 
     174        const int par = details->par_coord[k]; 
    175175        pvec[par] = values[offset[par]++]; 
    176176        //printf("p[%d]=v[%d]=%g ", par, offset[par]-1, pvec[par]); 
    177177        // if theta is coordinated with fast index, compute spherical correction each time 
    178         if (par == problem->theta_par) { 
    179           spherical_correction = fmax(fabs(cos(M_PI_180*pvec[problem->theta_par])), 1e-6); 
     178        if (par == details->theta_par) { 
     179          spherical_correction = fmax(fabs(cos(M_PI_180*pvec[details->theta_par])), 1e-6); 
    180180        } 
    181181      } 
    … …  
    209209 
    210210  // End of the PD loop we can normalize 
    211   if (pd_stop >= problem->total_pd) { 
     211  if (pd_stop >= details->total_pd) { 
    212212    const double scale = values[0]; 
    213213    const double background = values[1]; 

  • sasmodels/kernelpy.py

    r1e2a1ba ra8a7f08  
    8787    """ 
    8888    def __init__(self, kernel, model_info, q_input): 
     89        self.dtype = np.dtype('d') 
    8990        self.info = model_info 
    9091        self.q_input = q_input 
    9192        self.res = np.empty(q_input.nq, q_input.dtype) 
    92         dim = '2d' if q_input.is_2d else '1d' 
    93         # Loop over q unless user promises that the kernel is vectorized by 
    94         # taggining it with vectorized=True 
    95         if not getattr(kernel, 'vectorized', False): 
    96             if dim == '2d': 
    97                 def vector_kernel(qx, qy, *args): 
    98                     """ 
    99                     Vectorized 2D kernel. 
    100                     """ 
    101                     return np.array([kernel(qxi, qyi, *args) 
    102                                      for qxi, qyi in zip(qx, qy)]) 
     93        self.kernel = kernel 
     94        self.dim = '2d' if q_input.is_2d else '1d' 
     95 
     96        partable = model_info['parameters'] 
     97        kernel_parameters = (partable.iqxy_parameters if q_input.is_2d 
     98                             else partable.iq_parameters) 
     99        volume_parameters = partable.form_volume_parameters 
     100 
     101        # Create an array to hold the parameter values.  There will be a 
     102        # single array whose values are updated as the calculator goes 
     103        # through the loop.  Arguments to the kernel and volume functions 
     104        # will use views into this vector, relying on the fact that a 
     105        # an array of no dimensions acts like a scalar. 
     106        parameter_vector = np.empty(len(partable.call_parameters)-2, 'd') 
     107 
     108        # Create views into the array to hold the arguments 
     109        offset = 0 
     110        kernel_args, volume_args = [], [] 
     111        for p in partable.kernel_parameters: 
     112            if p.length == 1: 
     113                # Scalar values are length 1 vectors with no dimensions. 
     114                v = parameter_vector[offset:offset+1].reshape(()) 
    103115            else: 
    104                 def vector_kernel(q, *args): 
    105                     """ 
    106                     Vectorized 1D kernel. 
    107                     """ 
    108                     return np.array([kernel(qi, *args) 
    109                                      for qi in q]) 
    110             self.kernel = vector_kernel 
     116                # Vector values are simple views. 
     117                v = parameter_vector[offset:offset+p.length] 
     118            offset += p.length 
     119            if p in kernel_parameters: 
     120                kernel_args.append(v) 
     121            if p in volume_parameters: 
     122                volume_args.append(v) 
     123 
     124        # Hold on to the parameter vector so we can use it to call kernel later. 
     125        # This may also be required to preserve the views into the vector. 
     126        self._parameter_vector = parameter_vector 
     127 
     128        # Generate a closure which calls the kernel with the views into the 
     129        # parameter array. 
     130        if q_input.is_2d: 
     131            form = model_info['Iqxy'] 
     132            qx, qy = q_input.q[:,0], q_input.q[:,1] 
     133            self._form = lambda: form(qx, qy, *kernel_args) 
    111134        else: 
    112             self.kernel = kernel 
    113         fixed_pars = model_info['partype']['fixed-' + dim] 
    114         pd_pars = model_info['partype']['pd-' + dim] 
    115         vol_pars = model_info['partype']['volume'] 
    116  
    117         # First two fixed pars are scale and background 
    118         pars = [p.name for p in model_info['parameters'][2:]] 
    119         offset = len(self.q_input.q_vectors) 
    120         self.args = self.q_input.q_vectors + [None] * len(pars) 
    121         self.fixed_index = np.array([pars.index(p) + offset 
    122                                      for p in fixed_pars[2:]]) 
    123         self.pd_index = np.array([pars.index(p) + offset 
    124                                   for p in pd_pars]) 
    125         self.vol_index = np.array([pars.index(p) + offset 
    126                                    for p in vol_pars]) 
    127         try: self.theta_index = pars.index('theta') + offset 
    128         except ValueError: self.theta_index = -1 
    129  
    130         # Caller needs fixed_pars and pd_pars 
    131         self.fixed_pars = fixed_pars 
    132         self.pd_pars = pd_pars 
    133  
    134     def __call__(self, fixed, pd, cutoff=1e-5): 
    135         #print("fixed",fixed) 
    136         #print("pd", pd) 
    137         args = self.args[:]  # grab a copy of the args 
    138         form, form_volume = self.kernel, self.info['form_volume'] 
    139         # First two fixed 
    140         scale, background = fixed[:2] 
    141         for index, value in zip(self.fixed_index, fixed[2:]): 
    142             args[index] = float(value) 
    143         res = _loops(form, form_volume, cutoff, scale, background, args, 
    144                      pd, self.pd_index, self.vol_index, self.theta_index) 
    145  
     135            form = model_info['Iq'] 
     136            q = q_input.q 
     137            self._form = lambda: form(q, *kernel_args) 
     138 
     139        # Generate a closure which calls the form_volume if it exists. 
     140        form_volume = model_info['form_volume'] 
     141        self._volume = ((lambda: form_volume(*volume_args)) if form_volume 
     142                        else (lambda: 1.0)) 
     143 
     144    def __call__(self, details, weights, values, cutoff): 
     145        # type: (.generate.CoordinationDetails, np.ndarray, np.ndarray, float) -> np.ndarray 
     146        res = _loops(self._parameter_vector, self._form, self._volume, 
     147                     self.q_input.nq, details, weights, values, cutoff) 
    146148        return res 
    147149 
    … …  
    152154        self.q_input = None 
    153155 
    154 def _loops(form, form_volume, cutoff, scale, background, 
    155            args, pd, pd_index, vol_index, theta_index): 
    156     """ 
    157     *form* is the name of the form function, which should be vectorized over 
    158     q, but otherwise have an interface like the opencl kernels, with the 
    159     q parameters first followed by the individual parameters in the order 
    160     given in model.parameters (see :mod:`sasmodels.generate`). 
    161  
    162     *form_volume* calculates the volume of the shape.  *vol_index* gives 
    163     the list of volume parameters 
    164  
    165     *cutoff* ignores the corners of the dispersion hypercube 
    166  
    167     *scale*, *background* multiplies the resulting form and adds an offset 
    168  
    169     *args* is the prepopulated set of arguments to the form function, starting 
    170     with the q vectors, and including slots for all the parameters.  The 
    171     values for the parameters will be substituted with values from the 
    172     dispersion functions. 
    173  
    174     *pd* is the list of dispersion parameters 
    175  
    176     *pd_index* are the indices of the dispersion parameters in *args* 
    177  
    178     *vol_index* are the indices of the volume parameters in *args* 
    179  
    180     *theta_index* is the index of the theta parameter for the sasview 
    181     spherical correction, or -1 if there is no angular dispersion 
    182     """ 
    183  
     156def _loops(parameters,  # type: np.ndarray 
     157           form,        # type: Callable[[], np.ndarray] 
     158           form_volume, # type: Callable[[], float] 
     159           nq,          # type: int 
     160           details,     # type: .generate.CoordinationDetails 
     161           weights,     # type: np.ndarray 
     162           values,      # type: np.ndarray 
     163           cutoff,      # type: float 
     164           ):           # type: (...) -> None 
    184165    ################################################################ 
    185166    #                                                              # 
    … …  
    191172    #                                                              # 
    192173    ################################################################ 
    193  
    194     weight = np.empty(len(pd), 'd') 
    195     if weight.size > 0: 
    196         # weight vector, to be populated by polydispersity loops 
    197  
    198         # identify which pd parameters are volume parameters 
    199         vol_weight_index = np.array([(index in vol_index) for index in pd_index]) 
    200  
    201         # Sort parameters in decreasing order of pd length 
    202         Npd = np.array([len(pdi[0]) for pdi in pd], 'i') 
    203         order = np.argsort(Npd)[::-1] 
    204         stride = np.cumprod(Npd[order]) 
    205         pd = [pd[index] for index in order] 
    206         pd_index = pd_index[order] 
    207         vol_weight_index = vol_weight_index[order] 
    208  
    209         fast_value = pd[0][0] 
    210         fast_weight = pd[0][1] 
     174    parameters[:] = values[details.par_offset] 
     175    scale, background = values[0], values[1] 
     176    if details.num_active == 0: 
     177        norm = float(form_volume()) 
     178        if norm > 0.0: 
     179            return (scale/norm)*form() + background 
     180        else: 
     181            return np.ones(nq, 'd')*background 
     182 
     183    partial_weight = np.NaN 
     184    spherical_correction = 1.0 
     185    pd_stride = details.pd_stride[:details.num_active] 
     186    pd_length = details.pd_length[:details.num_active] 
     187    pd_offset = details.pd_offset[:details.num_active] 
     188    pd_index = np.empty_like(pd_offset) 
     189    offset = np.empty_like(details.par_offset) 
     190    theta = details.theta_par 
     191    fast_length = pd_length[0] 
     192    pd_index[0] = fast_length 
     193    total = np.zeros(nq, 'd') 
     194    norm = 0.0 
     195    for loop_index in range(details.total_pd): 
     196        # update polydispersity parameter values 
     197        if pd_index[0] == fast_length: 
     198            pd_index[:] = (loop_index/pd_stride)%pd_length 
     199            partial_weight = np.prod(weights[pd_offset+pd_index][1:]) 
     200            for k in range(details.num_coord): 
     201                par = details.par_coord[k] 
     202                coord = details.pd_coord[k] 
     203                this_offset = details.par_offset[par] 
     204                block_size = 1 
     205                for bit in xrange(32): 
     206                    if coord&1: 
     207                        this_offset += block_size * pd_index[bit] 
     208                        block_size *= pd_length[bit] 
     209                    coord >>= 1 
     210                    if coord == 0: break 
     211                offset[par] = this_offset 
     212                parameters[par] = values[this_offset] 
     213                if par == theta and not (details.par_coord[k]&1): 
     214                    spherical_correction = max(abs(cos(pi/180 * parameters[theta])), 1e-6) 
     215        for k in range(details.num_coord): 
     216            if details.pd_coord[k]&1: 
     217                #par = details.par_coord[k] 
     218                parameters[par] = values[offset[par]] 
     219                #print "par",par,offset[par],parameters[par+2] 
     220                offset[par] += 1 
     221                if par == theta: 
     222                    spherical_correction = max(abs(cos(pi/180 * parameters[theta])), 1e-6) 
     223 
     224        weight = partial_weight * weights[pd_offset[0] + pd_index[0]] 
     225        pd_index[0] += 1 
     226        if weight > cutoff: 
     227            # Call the scattering function 
     228            # Assume that NaNs are only generated if the parameters are bad; 
     229            # exclude all q for that NaN.  Even better would be to have an 
     230            # INVALID expression like the C models, but that is too expensive. 
     231            I = form() 
     232            if np.isnan(I).any(): continue 
     233 
     234            # update value and norm 
     235            weight *= spherical_correction 
     236            total += weight * I 
     237            norm += weight * form_volume() 
     238 
     239    if norm > 0.0: 
     240        return (scale/norm)*total + background 
    211241    else: 
    212         stride = np.array([1]) 
    213         vol_weight_index = slice(None, None) 
    214         # keep lint happy 
    215         fast_value = [None] 
    216         fast_weight = [None] 
    217  
    218     ret = np.zeros_like(args[0]) 
    219     norm = np.zeros_like(ret) 
    220     vol = np.zeros_like(ret) 
    221     vol_norm = np.zeros_like(ret) 
    222     for k in range(stride[-1]): 
    223         # update polydispersity parameter values 
    224         fast_index = k % stride[0] 
    225         if fast_index == 0:  # bottom loop complete ... check all other loops 
    226             if weight.size > 0: 
    227                 for i, index, in enumerate(k % stride): 
    228                     args[pd_index[i]] = pd[i][0][index] 
    229                     weight[i] = pd[i][1][index] 
    230         else: 
    231             args[pd_index[0]] = fast_value[fast_index] 
    232             weight[0] = fast_weight[fast_index] 
    233  
    234         # Computes the weight, and if it is not sufficient then ignore this 
    235         # parameter set. 
    236         # Note: could precompute w1*...*wn so we only need to multiply by w0 
    237         w = np.prod(weight) 
    238         if w > cutoff: 
    239             # Note: can precompute spherical correction if theta_index is not 
    240             # the fast index. Correction factor for spherical integration 
    241             #spherical_correction = abs(cos(pi*args[phi_index])) if phi_index>=0 else 1.0 
    242             spherical_correction = (abs(cos(pi * args[theta_index])) * pi / 2 
    243                                     if theta_index >= 0 else 1.0) 
    244             #spherical_correction = 1.0 
    245  
    246             # Call the scattering function and adds it to the total. 
    247             I = form(*args) 
    248             if np.isnan(I).any(): continue 
    249             ret += w * I * spherical_correction 
    250             norm += w 
    251  
    252             # Volume normalization. 
    253             # If there are "volume" polydispersity parameters, then these 
    254             # will be used to call the form_volume function from the user 
    255             # supplied kernel, and accumulate a normalized weight. 
    256             # Note: can precompute volume norm if fast index is not a volume 
    257             if form_volume: 
    258                 vol_args = [args[index] for index in vol_index] 
    259                 vol_weight = np.prod(weight[vol_weight_index]) 
    260                 vol += vol_weight * form_volume(*vol_args) 
    261                 vol_norm += vol_weight 
    262  
    263     positive = (vol * vol_norm != 0.0) 
    264     ret[positive] *= vol_norm[positive] / vol[positive] 
    265     result = scale * ret / norm + background 
    266     return result 
     242        return np.ones(nq, 'd')*background 

  • sasmodels/model_test.py

    rd2bb604 r6e7ff6d  
    8989 
    9090        if is_py:  # kernel implemented in python 
    91             continue # TODO: re-enable python tests 
    9291            test_name = "Model: %s, Kernel: python"%model_name 
    9392            test_method_name = "test_%s_python" % model_name 

  • sasmodels/modelinfo.py

    ree8f734 r6e7ff6d  
    189189    can automatically be promoted to magnetic parameters, each of which 
    190190    will have a magnitude and a direction, which may be different from 
    191     other sld parameters. 
     191    other sld parameters. The volume parameters are used for calls 
     192    to form_volume within the kernel (required for volume normalization) 
     193    and for calls to ER and VR for effective radius and volume ratio 
     194    respectively. 
    192195 
    193196    *description* is a short description of the parameter.  This will 
    … …  
    197200    Additional values can be set after the parameter is created: 
    198201 
    199     *length* is the length of the field if it is a vector field 
    200     *length_control* is the parameter which sets the vector length 
    201     *is_control* is True if the parameter is a control parameter for a vector 
    202     *polydisperse* is true if the parameter accepts a polydispersity 
    203     *relative_pd* is true if that polydispersity is relative 
     202    * *length* is the length of the field if it is a vector field 
     203 
     204    * *length_control* is the parameter which sets the vector length 
     205 
     206    * *is_control* is True if the parameter is a control parameter for a vector 
     207 
     208    * *polydisperse* is true if the parameter accepts a polydispersity 
     209 
     210    * *relative_pd* is true if that polydispersity is a portion of the 
     211    value (so a 10% length dipsersity would use a polydispersity value of 0.1) 
     212    rather than absolute dispersisity (such as an angle plus or minus 
     213    15 degrees). 
    204214 
    205215    In the usual process these values are set by :func:`make_parameter_table` 

  • sasmodels/models/broad_peak.py

    rec45c4f r65279d8  
    101101    """ 
    102102    return Iq(sqrt(qx ** 2 + qy ** 2), *args) 
    103  
    104103Iqxy.vectorized = True # Iqxy accepts an array of qx, qy values 
    105104 

  • sasmodels/models/poly_gauss_coil.py

    rec45c4f r65279d8  
    5050""" 
    5151 
    52 from numpy import inf, sqrt, exp, power 
     52import numpy as np 
     53from numpy import inf, exp, power, sqrt 
    5354 
    5455name =  "poly_gauss_coil" 
    … …  
    7071    # pylint: disable = missing-docstring 
    7172    u = polydispersity - 1.0 
    72     z = ((q * radius_gyration) * (q * radius_gyration)) / (1.0 + 2.0 * u) 
    73     if (q == 0).any(): 
    74         inten = i_zero 
     73    z = (q*radius_gyration)**2 / (1.0 + 2.0*u) 
     74    # need to trap the case of the polydispersity being 1 (ie, monodispersity!) 
     75    if polydispersity == 1.0: 
     76        inten = i_zero * 2.0 * (exp(-z) + z - 1.0) 
    7577    else: 
    76     # need to trap the case of the polydispersity being 1 (ie, monodispersity!) 
    77         if polydispersity == 1: 
    78             inten = i_zero * 2.0 * (exp(-z) + z - 1.0 ) / (z * z) 
    79         else: 
    80             minusoneonu = -1.0 / u 
    81             inten = i_zero * 2.0 * (power((1.0 + u * z),minusoneonu) + z - 1.0 ) / ((1.0 + u) * (z * z)) 
     78        inten = i_zero * 2.0 * (power(1.0 + u*z, -1.0/u) + z - 1.0) / (1.0 + u) 
     79    index = q != 0. 
     80    inten[~index] = i_zero 
     81    inten[index] /= z[index]**2 
    8282    return inten 
    83 #Iq.vectorized =  True  # Iq accepts an array of q values 
     83Iq.vectorized =  True  # Iq accepts an array of q values 
    8484 
    8585def Iqxy(qx, qy, *args): 
    8686    # pylint: disable = missing-docstring 
    8787    return Iq(sqrt(qx ** 2 + qy ** 2), *args) 
    88 #Iqxy.vectorized = True # Iqxy accepts an array of qx, qy values 
     88Iqxy.vectorized = True # Iqxy accepts an array of qx, qy values 
    8989 
    9090demo =  dict(scale = 1.0, 
Note: See TracChangeset for help on using the changeset viewer.