Changeset 6e7ff6d in sasmodels
- Timestamp:
- Apr 7, 2016 2:01:54 PM (9 years ago)
- 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:
- 3543141
- Parents:
- 8bd7b77
- Location:
- sasmodels
- Files:
-
- 6 edited
Legend:
- Unmodified
- Added
- Removed
-
sasmodels/convert.json
r8bd7b77 r6e7ff6d 625 625 } 626 626 ], 627 "_spherepy": [ 628 "SphereModel", 629 { 630 "sld": "sldSph", 631 "radius": "radius", 632 "sld_solvent": "sldSolv" 633 } 634 ], 627 635 "spherical_sld": [ 628 636 "SphereSLDModel", -
sasmodels/generate.py
r4937980 r6e7ff6d 125 125 An *model_info* dictionary is constructed from the kernel meta data and 126 126 returned 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 for130 the non-polydisperse parameters, followed by the lengths of the131 polydispersity loops. To construct the call for 1D models, the132 categories *fixed-1d* and *pd-1d* list the names of the parameters133 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 give136 polydispersitiy as a portion of the value (so a 10% length dispersity137 would use a polydispersity value of 0.1) rather than absolute138 dispersity such as an angle plus or minus 15 degrees.139 140 The *volume* category lists the volume parameters in order for calls141 to volume within the kernel (used for volume normalization) and for142 calls to ER and VR for effective radius and volume ratio respectively.143 144 The *orientation* and *magnetic* categories list the orientation and145 magnetic parameters. These are used by the sasview interface. The146 blank category is for parameters such as scale which don't have any147 other marking.148 127 149 128 The doc string at the start of the kernel module will be used to … … 552 531 @property 553 532 def ctypes(self): return self.details.ctypes 533 554 534 @property 555 535 def pd_par(self): return self._pd_par 536 556 537 @property 557 538 def pd_length(self): return self._pd_length 539 558 540 @property 559 541 def pd_offset(self): return self._pd_offset 542 560 543 @property 561 544 def pd_stride(self): return self._pd_stride 545 562 546 @property 563 547 def pd_coord(self): return self._pd_coord 548 564 549 @property 565 550 def par_coord(self): return self._par_coord 551 566 552 @property 567 553 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 568 565 @property 569 566 def num_coord(self): return self.details[-2] 570 567 @num_coord.setter 571 568 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] 580 572 581 573 def show(self): … … 638 630 639 631 632 def 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 645 def 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 640 666 def create_default_functions(model_info): 641 667 """ … … 643 669 644 670 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 656 676 657 677 def load_kernel_module(model_name): -
sasmodels/kernel_iq.c
r398aa94 r6e7ff6d 42 42 const int32_t pd_start, // where we are in the polydispersity loop 43 43 const int32_t pd_stop, // where we are stopping in the polydispersity loop 44 global const ProblemDetails * problem,44 global const ProblemDetails *details, 45 45 global const double *weights, 46 46 global const double *values, … … 60 60 #endif 61 61 for (int k=0; k < NPARS; k++) { 62 pvec[k] = values[ problem->par_offset[k]];62 pvec[k] = values[details->par_offset[k]]; 63 63 } 64 64 … … 78 78 79 79 // Monodisperse computation 80 if ( problem->num_active == 0) {80 if (details->num_active == 0) { 81 81 #ifdef INVALID 82 82 if (INVALID(local_values)) { return; } … … 116 116 // Trigger the reset behaviour that happens at the end the fast loop 117 117 // 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]; 119 119 pd_index[0] = fast_length; 120 120 … … 123 123 // check if indices need to be updated 124 124 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); 126 126 127 127 // 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]; 130 130 //printf("pd_index[%d] = %d\n",k,pd_index[k]); 131 131 } … … 134 134 partial_weight = 1.0; 135 135 //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); 139 139 partial_weight *= wi; 140 140 } … … 143 143 // Update parameter offsets in weight vector 144 144 //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]; 149 149 int block_size = 1; 150 150 for (int bit=0; coord != 0; bit++) { 151 151 if (coord&1) { 152 152 this_offset += block_size * pd_index[bit]; 153 block_size *= problem->pd_length[bit];153 block_size *= details->pd_length[bit]; 154 154 } 155 155 coord >>= 1; … … 159 159 //printf("par[%d]=v[%d]=%g \n", k, offset[k], pvec[k]); 160 160 // 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); 163 163 } 164 164 } … … 167 167 168 168 // 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]++]; 170 170 double weight = partial_weight*wi; 171 171 //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]; 175 175 pvec[par] = values[offset[par]++]; 176 176 //printf("p[%d]=v[%d]=%g ", par, offset[par]-1, pvec[par]); 177 177 // 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); 180 180 } 181 181 } … … 209 209 210 210 // End of the PD loop we can normalize 211 if (pd_stop >= problem->total_pd) {211 if (pd_stop >= details->total_pd) { 212 212 const double scale = values[0]; 213 213 const double background = values[1]; … … 216 216 #endif 217 217 for (int i=0; i < nq; i++) { 218 result[i] = (norm>0. ? scale* scattering/norm + background : background);218 result[i] = (norm>0. ? scale*result[i]/norm + background : background); 219 219 } 220 220 } -
sasmodels/kernelpy.py
r1e2a1ba r6e7ff6d 87 87 """ 88 88 def __init__(self, kernel, model_info, q_input): 89 self.dtype = np.dtype('d') 89 90 self.info = model_info 90 91 self.q_input = q_input 91 92 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), 'd') 107 108 # Create views into the array to hold the arguments 109 offset = 2 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(()) 103 115 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(p) 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) 111 134 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) 146 148 return res 147 149 … … 152 154 self.q_input = None 153 155 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 156 def _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 184 165 ################################################################ 185 166 # # … … 191 172 # # 192 173 ################################################################ 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[2:] = 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 offset[par] += 1 220 if par == theta: 221 spherical_correction = max(abs(cos(pi/180 * parameters[theta])), 1e-6) 222 223 weight = partial_weight * weights[pd_offset[0] + pd_index[0]] 224 pd_index[0] += 1 225 if weight > cutoff: 226 # Call the scattering function 227 # Assume that NaNs are only generated if the parameters are bad; 228 # exclude all q for that NaN. Even better would be to have an 229 # INVALID expression like the C models, but that is too expensive. 230 I = form() 231 if np.isnan(I).any(): continue 232 233 # update value and norm 234 weight *= spherical_correction 235 total += weight * I 236 norm += weight * form_volume() 237 238 if norm > 0.0: 239 return (scale/norm)*total + background 211 240 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 241 return np.ones(nq, 'd')*background -
sasmodels/model_test.py
rd2bb604 r6e7ff6d 89 89 90 90 if is_py: # kernel implemented in python 91 continue # TODO: re-enable python tests92 91 test_name = "Model: %s, Kernel: python"%model_name 93 92 test_method_name = "test_%s_python" % model_name -
sasmodels/modelinfo.py
ree8f734 r6e7ff6d 189 189 can automatically be promoted to magnetic parameters, each of which 190 190 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. 192 195 193 196 *description* is a short description of the parameter. This will … … 197 200 Additional values can be set after the parameter is created: 198 201 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). 204 214 205 215 In the usual process these values are set by :func:`make_parameter_table`
Note: See TracChangeset
for help on using the changeset viewer.