# Changeset 6b7e2f7 in sasmodels

Ignore:
Timestamp:
Aug 17, 2016 3:59:36 PM (5 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:
6dc78e4
Parents:
fe496dd (diff), 524e5c4 (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 'master' into polydisp

Files:
19 edited
2 moved

Unmodified
Removed
• ## doc/ref/magnetism/magnetism.rst

 rb1c40bee For magnetic scattering, only the magnetization component $\mathbf{M_\perp}$ perpendicular to the scattering vector $\mathbf{Q}$ contributes to the the magnetic perpendicular to the scattering vector $\mathbf{Q}$ contributes to the magnetic scattering length. .. figure:: img/mag_vector.bmp mag_img/mag_vector.bmp The magnetic scattering length density is then .. figure:: img/M_angles_pic.bmp mag_img/M_angles_pic.bmp If the angles of the $Q$ vector and the spin-axis $x'$ to the $x$ - axis are $\phi$ and $\theta_{up}$, respectively, then, depending on the spin state of the neutrons, the scattering length densities, including the nuclear scattering length density ($\beta_N$) are length density ($\beta{_N}$) are .. math:: ===========   ================================================================ M0_sld        = $D_M M_0$ Up_theta      = $\theta_up$ Up_theta      = $\theta_{up}$ M_theta       = $\theta_M$ M_phi         = $\phi_M$
• ## sasmodels/models/core_multi_shell.py

 rb1c40bee For information about polarised and magnetic scattering, see the :doc:magnetic help <../sasgui/perspectives/fitting/mag_help> documentation. the :ref:magnetism documentation. Our model uses the form factor calculations implemented in a c-library provided
• ## sasmodels/models/core_shell_sphere.py

 rb1c40bee For information about polarised and magnetic scattering, see the :doc:magnetic help <../sasgui/perspectives/fitting/mag_help> documentation. the :ref:magnetism documentation. Definition
• ## sasmodels/models/cylinder.py

 rb1c40bee The form factor is normalized by the particle volume V = \piR^2L. For information about polarised and magnetic scattering, see the :doc:magnetic help <../sasgui/perspectives/fitting/mag_help> documentation. the :ref:magnetism documentation. Definition
• ## sasmodels/models/fuzzy_sphere.py

 rb1c40bee r""" For information about polarised and magnetic scattering, see the :doc:magnetic help <../sasgui/perspectives/fitting/mag_help> documentation. the :ref:magnetism documentation. Definition
• ## sasmodels/models/multilayer_vesicle.py

 rb1c40bee For information about polarised and magnetic scattering, see the :doc:magnetic help <../sasgui/perspectives/fitting/mag_help> documentation. the :ref:magnetism documentation. This code is based on the form factor calculations implemented in the NIST
• ## sasmodels/models/parallelepiped.py

 rb1c40bee The form factor is normalized by the particle volume. For information about polarised and magnetic scattering, see the :doc:magnetic help <../sasgui/perspectives/fitting/mag_help> documentation. the :ref:magnetism documentation. Definition
• ## sasmodels/models/sphere.py

 rb1c40bee r""" For information about polarised and magnetic scattering, see the :doc:magnetic help <../sasgui/perspectives/fitting/mag_help> documentation. the :ref:magnetism documentation. documentation.
• ## sasmodels/core.py

 r672978c except ImportError: pass try: np.meshgrid([]) meshgrid = np.meshgrid except Exception: # CRUFT: np.meshgrid requires multiple vectors def meshgrid(*args): """Allow meshgrid with a single argument""" if len(args) > 1: return np.meshgrid(*args) else: return [np.asarray(v) for v in args] # TODO: refactor composite model support
• ## sasmodels/details.py

 r40a87fa np.meshgrid([]) meshgrid = np.meshgrid except ValueError: except Exception: # CRUFT: np.meshgrid requires multiple vectors def meshgrid(*args): #   pd_offset[max_pd]  offset of pd values in parameter array #   pd_stride[max_pd]  index of pd value in loop = n//stride[k] #   pd_prod            total length of pd loop #   pd_sum             total length of the weight vector #   num_eval           total length of pd loop #   num_weights        total length of the weight vector #   num_active         number of pd params #   theta_par          parameter number for theta parameter self.buffer = np.zeros(4*max_pd + 4, 'i4') self.buffer = np.empty(4*max_pd + 4, 'i4') # generate views on different parts of the array self.theta_par = parameters.theta_offset # offset and length are for all parameters, not just pd parameters # They are not sent to the kernel function, though they could be. # They are used by the composite models (sum and product) to # figure out offsets into the combined value list. self.offset = None  # type: np.ndarray self.length = None  # type: np.ndarray # keep hold of ifno show() so we can break a values vector # into the individual components self.info = model_info @property def pd_par(self): @property def pd_prod(self): def num_eval(self): """Total size of the pd mesh""" return self.buffer[-4] @pd_prod.setter def pd_prod(self, v): @num_eval.setter def num_eval(self, v): """Total size of the pd mesh""" self.buffer[-4] = v @property def pd_sum(self): def num_weights(self): """Total length of all the weight vectors""" return self.buffer[-3] @pd_sum.setter def pd_sum(self, v): @num_weights.setter def num_weights(self, v): """Total length of all the weight vectors""" self.buffer[-3] = v self.buffer[-1] = v def show(self): def show(self, values=None): """Print the polydispersity call details to the console""" print("num_active", self.num_active) print("pd_prod", self.pd_prod) print("pd_sum", self.pd_sum) print("theta par", self.theta_par) print("pd_par", self.pd_par) print("pd_length", self.pd_length) print("pd_offset", self.pd_offset) print("pd_stride", self.pd_stride) def mono_details(model_info): # type: (ModelInfo) -> CallDetails """ Return a :class:CallDetails object for a monodisperse calculation of the model defined by *model_info*. """ call_details = CallDetails(model_info) call_details.pd_prod = 1 call_details.pd_sum = model_info.parameters.nvalues call_details.pd_par[:] = np.arange(0, model_info.parameters.max_pd) call_details.pd_length[:] = 1 call_details.pd_offset[:] = np.arange(0, model_info.parameters.max_pd) call_details.pd_stride[:] = 1 return call_details def poly_details(model_info, weights): print("===== %s details ===="%self.info.name) print("num_active:%d  num_eval:%d  num_weights:%d  theta=%d" % (self.num_active, self.num_eval, self.num_weights, self.theta_par)) if self.pd_par.size: print("pd_par", self.pd_par) print("pd_length", self.pd_length) print("pd_offset", self.pd_offset) print("pd_stride", self.pd_stride) if values is not None: nvalues = self.info.parameters.nvalues print("scale, background", values[:2]) print("val", values[2:nvalues]) print("pd", values[nvalues:nvalues+self.num_weights]) print("wt", values[nvalues+self.num_weights:nvalues+2*self.num_weights]) print("offsets", self.offset) def make_details(model_info, length, offset, num_weights): """ Return a :class:CallDetails object for a polydisperse calculation of the model defined by *model_info* for the given set of *weights*. *weights* is a list of vectors, one for each parameter.  Monodisperse parameters should provide a weight vector containing [1.0]. """ # type: (ModelInfo) -> CallDetails #print("weights",weights) #weights = weights[2:] # Skip scale and background # Decreasing list of polydispersity lengths #print([p.id for p in model_info.parameters.call_parameters]) pd_length = np.array([len(w) for w in weights[2:2+model_info.parameters.npars]]) num_active = np.sum(pd_length > 1) of the model defined by *model_info*.  Polydispersity is defined by the *length* of the polydispersity distribution for each parameter and the *offset* of the distribution in the polydispersity array. Monodisperse parameters should use a polydispersity length of one with weight 1.0. *num_weights* is the total length of the polydispersity array. """ # type: (ModelInfo, np.ndarray, np.ndarray, int) -> CallDetails #pars = model_info.parameters.call_parameters[2:model_info.parameters.npars+2] #print(", ".join(str(i)+"-"+p.id for i,p in enumerate(pars))) #print("len:",length) #print("off:",offset) # Check that we arn't using too many polydispersity loops num_active = np.sum(length > 1) max_pd = model_info.parameters.max_pd if num_active > max_pd: raise ValueError("Too many polydisperse parameters") pd_offset = np.cumsum(np.hstack((0, pd_length))) #print(", ".join(str(i)+"-"+p.id for i,p in enumerate(model_info.parameters.call_parameters))) #print("len:",pd_length) #print("off:",pd_offset) # Decreasing list of polydpersity lengths # Note: the reversing view, x[::-1], does not require a copy idx = np.argsort(pd_length)[::-1][:max_pd] pd_stride = np.cumprod(np.hstack((1, pd_length[idx]))) idx = np.argsort(length)[::-1][:max_pd] pd_stride = np.cumprod(np.hstack((1, length[idx]))) call_details = CallDetails(model_info) call_details.pd_par[:max_pd] = idx call_details.pd_length[:max_pd] = pd_length[idx] call_details.pd_offset[:max_pd] = pd_offset[idx] call_details.pd_length[:max_pd] = length[idx] call_details.pd_offset[:max_pd] = offset[idx] call_details.pd_stride[:max_pd] = pd_stride[:-1] call_details.pd_prod = pd_stride[-1] call_details.pd_sum = sum(len(w) for w in weights) call_details.num_eval = pd_stride[-1] call_details.num_weights = num_weights call_details.num_active = num_active call_details.length = length call_details.offset = offset #call_details.show() return call_details ZEROS = tuple([0.]*31) def make_kernel_args(kernel, pairs): # type: (Kernel, Tuple[List[np.ndarray], List[np.ndarray]]) -> Tuple[CallDetails, np.ndarray, bool] """ Converts (value, weight) pairs into parameters for the kernel call. Returns a CallDetails object indicating the polydispersity, a data object containing the different values, and the magnetic flag indicating whether any magnetic magnitudes are non-zero. Magnetic vectors (M0, phi, theta) are converted to rectangular coordinates (mx, my, mz). """ npars = kernel.info.parameters.npars nvalues = kernel.info.parameters.nvalues scalars = [v[0][0] for v in pairs] values, weights = zip(*pairs[2:npars+2]) if npars else ((),()) length = np.array([len(w) for w in weights]) offset = np.cumsum(np.hstack((0, length))) call_details = make_details(kernel.info, length, offset[:-1], offset[-1]) # Pad value array to a 32 value boundaryd data_len = nvalues + 2*sum(len(v) for v in values) extra = (32 - data_len%32)%32 data = np.hstack((scalars,) + values + weights + ZEROS[:extra]) data = data.astype(kernel.dtype) is_magnetic = convert_magnetism(kernel.info.parameters, data) #call_details.show() return call_details, data, is_magnetic def convert_magnetism(parameters, values): """ Convert magnetism values from polar to rectangular coordinates. Returns True if any magnetism is present. """ mag = values[parameters.nvalues-3*parameters.nmagnetic:parameters.nvalues] mag = mag.reshape(-1, 3) scale = mag[:,0] if np.any(scale): theta, phi = mag[:, 1]*pi/180., mag[:, 2]*pi/180. cos_theta = cos(theta) mag[:, 0] = scale*cos_theta*cos(phi)  # mx mag[:, 1] = scale*sin(theta)  # my mag[:, 2] = -scale*cos_theta*sin(phi)  # mz return True else: return False value = pars return value, weight def build_details(kernel, pairs): # type: (Kernel, Tuple[List[np.ndarray], List[np.ndarray]]) -> Tuple[CallDetails, np.ndarray, bool] """ Converts (value, weight) pairs into parameters for the kernel call. Returns a CallDetails object indicating the polydispersity, a data object containing the different values, and the magnetic flag indicating whether any magnetic magnitudes are non-zero. Magnetic vectors (M0, phi, theta) are converted to rectangular coordinates (mx, my, mz). """ values, weights = zip(*pairs) scalars = [v[0] for v in values] if all(len(w) == 1 for w in weights): call_details = mono_details(kernel.info) # Pad value array to a 32 value boundary data_len = 3*len(scalars) extra = ((data_len+31)//32)*32 - data_len data = np.array(scalars+scalars+[1.]*len(scalars)+[0.]*extra, dtype=kernel.dtype) else: call_details = poly_details(kernel.info, weights) # Pad value array to a 32 value boundary data_len = len(scalars) + 2*sum(len(v) for v in values) extra = ((data_len+31)//32)*32 - data_len data = np.hstack(scalars+list(values)+list(weights)+[0.]*extra) data = data.astype(kernel.dtype) is_magnetic = convert_magnetism(kernel.info.parameters, data) #call_details.show() return call_details, data, is_magnetic def convert_magnetism(parameters, values): """ Convert magnetism values from polar to rectangular coordinates. Returns True if any magnetism is present. """ mag = values[parameters.nvalues-3*parameters.nmagnetic:parameters.nvalues] mag = mag.reshape(-1, 3) scale = mag[:,0] if np.any(scale): theta, phi = mag[:, 1]*pi/180., mag[:, 2]*pi/180. cos_theta = cos(theta) mag[:, 0] = scale*cos_theta*cos(phi)  # mx mag[:, 1] = scale*sin(theta)  # my mag[:, 2] = -scale*cos_theta*sin(phi)  # mz return True else: return False
• ## sasmodels/direct_model.py

 r40a87fa from . import resolution from . import resolution2d from .details import build_details, dispersion_mesh from .details import make_kernel_args, dispersion_mesh try: for p in parameters.call_parameters] call_details, values, is_magnetic = build_details(calculator, vw_pairs) call_details, values, is_magnetic = make_kernel_args(calculator, vw_pairs) #print("values:", values) return calculator(call_details, values, cutoff, is_magnetic)
• ## sasmodels/kernel.py

 r9eb3632 info = None  # type: ModelInfo results = None # type: List[np.ndarray] dtype = None  # type: np.dtype def __call__(self, call_details, values, cutoff, magnetic):
• ## sasmodels/kernel_iq.c

 r0f00d95 int32_t pd_stride[MAX_PD];  // stride to move to the next index at this level #endif // MAX_PD > 0 int32_t pd_prod;            // total number of voxels in hypercube int32_t pd_sum;             // total length of the weights vector int32_t num_eval;           // total number of voxels in hypercube int32_t num_weights;        // total length of the weights vector int32_t num_active;         // number of non-trivial pd loops int32_t theta_par;          // id of spherical correction variable } ProblemDetails; // Intel HD 4000 needs private arrays to be a multiple of 4 long typedef struct { PARAMETER_TABLE } ParameterTable; typedef union { ParameterTable table; double vector[4*((NUM_PARS+3)/4)]; } ParameterBlock; #endif // _PAR_BLOCK_ ) { // Storage for the current parameter values.  These will be updated as we // walk the polydispersity cube.  local_values will be aliased to pvec. // walk the polydispersity cube. ParameterBlock local_values; double *pvec = (double *)&local_values; #if defined(MAGNETIC) && NUM_MAGNETIC>0 // Location of the sld parameters in the parameter pvec. // Location of the sld parameters in the parameter vector. // These parameters are updated with the effective sld due to magnetism. #if NUM_MAGNETIC > 3 #endif for (int i=0; i < NUM_PARS; i++) { pvec[i] = values[2+i]; //printf("p%d = %g\n",i, pvec[i]); } double pd_norm; //printf("start: %d %d\n",pd_start, pd_stop); local_values.vector[i] = values[2+i]; //printf("p%d = %g\n",i, local_values.vector[i]); } //printf("NUM_VALUES:%d  NUM_PARS:%d  MAX_PD:%d\n", NUM_VALUES, NUM_PARS, MAX_PD); //printf("start:%d  stop:%d\n", pd_start, pd_stop); double pd_norm = (pd_start == 0 ? 0.0 : result[nq]); if (pd_start == 0) { pd_norm = 0.0; #ifdef USE_OPENMP #pragma omp parallel for #endif for (int q_index=0; q_index < nq; q_index++) result[q_index] = 0.0; //printf("initializing %d\n", nq); } else { pd_norm = result[nq]; } //printf("start %d %g %g\n", pd_start, pd_norm, result[0]); #if MAX_PD>0 global const double *pd_value = values + NUM_VALUES + 2; global const double *pd_weight = pd_value + details->pd_sum; global const double *pd_value = values + NUM_VALUES; global const double *pd_weight = pd_value + details->num_weights; #endif global const double *v0 = pd_value + details->pd_offset[0]; global const double *w0 = pd_weight + details->pd_offset[0]; //printf("w0:%p, values:%p, diff:%d, %d\n",w0,values,(w0-values),NUM_VALUES); //printf("w0:%p, values:%p, diff:%ld, %d\n",w0,values,(w0-values), NUM_VALUES); #endif int step = pd_start; #if MAX_PD>4 const double weight5 = 1.0; while (i4 < n4) { pvec[p4] = v4[i4]; local_values.vector[p4] = v4[i4]; double weight4 = w4[i4] * weight5; //printf("step:%d level %d: p:%d i:%d n:%d value:%g weight:%g\n", step, 4, p4, i4, n4, pvec[p4], weight4); //printf("step:%d level %d: p:%d i:%d n:%d value:%g weight:%g\n", step, 4, p4, i4, n4, local_values.vector[p4], weight4); #elif MAX_PD>3 const double weight4 = 1.0; #if MAX_PD>3 while (i3 < n3) { pvec[p3] = v3[i3]; local_values.vector[p3] = v3[i3]; double weight3 = w3[i3] * weight4; //printf("step:%d level %d: p:%d i:%d n:%d value:%g weight:%g\n", step, 3, p3, i3, n3, pvec[p3], weight3); //printf("step:%d level %d: p:%d i:%d n:%d value:%g weight:%g\n", step, 3, p3, i3, n3, local_values.vector[p3], weight3); #elif MAX_PD>2 const double weight3 = 1.0; #if MAX_PD>2 while (i2 < n2) { pvec[p2] = v2[i2]; local_values.vector[p2] = v2[i2]; double weight2 = w2[i2] * weight3; //printf("step:%d level %d: p:%d i:%d n:%d value:%g weight:%g\n", step, 2, p2, i2, n2, pvec[p2], weight2); //printf("step:%d level %d: p:%d i:%d n:%d value:%g weight:%g\n", step, 2, p2, i2, n2, local_values.vector[p2], weight2); #elif MAX_PD>1 const double weight2 = 1.0; #if MAX_PD>1 while (i1 < n1) { pvec[p1] = v1[i1]; local_values.vector[p1] = v1[i1]; double weight1 = w1[i1] * weight2; //printf("step:%d level %d: p:%d i:%d n:%d value:%g weight:%g\n", step, 1, p1, i1, n1, pvec[p1], weight1); //printf("step:%d level %d: p:%d i:%d n:%d value:%g weight:%g\n", step, 1, p1, i1, n1, local_values.vector[p1], weight1); #elif MAX_PD>0 const double weight1 = 1.0; #if MAX_PD>0 if (slow_theta) { // Theta is not in inner loop spherical_correction = fmax(fabs(cos(M_PI_180*pvec[theta_par])), 1.e-6); spherical_correction = fmax(fabs(cos(M_PI_180*local_values.vector[theta_par])), 1.e-6); } while(i0 < n0) { pvec[p0] = v0[i0]; local_values.vector[p0] = v0[i0]; double weight0 = w0[i0] * weight1; //printf("step:%d level %d: p:%d i:%d n:%d value:%g weight:%g\n", step, 0, p0, i0, n0, pvec[p0], weight0); //printf("step:%d level %d: p:%d i:%d n:%d value:%g weight:%g\n", step, 0, p0, i0, n0, local_values.vector[p0], weight0); if (fast_theta) { // Theta is in inner loop spherical_correction = fmax(fabs(cos(M_PI_180*pvec[p0])), 1.e-6); spherical_correction = fmax(fabs(cos(M_PI_180*local_values.vector[p0])), 1.e-6); } #else #endif //printf("step:%d of %d, pars:",step,pd_stop); for (int i=0; i < NUM_PARS; i++) printf("p%d=%g ",i, pvec[i]); printf("\n"); //printf("step:%d of %d, pars:",step,pd_stop); for (int i=0; i < NUM_PARS; i++) printf("p%d=%g ",i, local_values.vector[i]); printf("\n"); //printf("sphcor: %g\n", spherical_correction); #ifdef INVALID if (!INVALID(local_values)) if (!INVALID(local_values.table)) #endif { // would be problems looking at models with theta=90. const double weight = weight0 * spherical_correction; pd_norm += weight * CALL_VOLUME(local_values); pd_norm += weight * CALL_VOLUME(local_values.table); #ifdef USE_OPENMP // TODO: what is the magnetic scattering at q=0 if (qsq > 1.e-16) { double p[4]; double p[4];  // dd, du, ud, uu p[0] = (qy*cos_mspin + qx*sin_mspin)/qsq; p[3] = -p[0]; #define M3 NUM_PARS+13 #define SLD(_M_offset, _sld_offset) \ pvec[_sld_offset] = xs * (axis \ local_values.vector[_sld_offset] = xs * (axis \ ? (index==1 ? -values[_M_offset+2] : values[_M_offset+2]) \ : mag_sld(qx, qy, pk, values[_M_offset], values[_M_offset+1], \ } #endif scattering += CALL_IQ(q, q_index, local_values); scattering += CALL_IQ(q, q_index, local_values.table); } } } #else  // !MAGNETIC const double scattering = CALL_IQ(q, q_index, local_values); const double scattering = CALL_IQ(q, q_index, local_values.table); #endif // !MAGNETIC //printf("q_index:%d %g %g %g %g\n",q_index, scattering, weight, spherical_correction, weight0);
• ## sasmodels/kernel_iq.cl

 r4f1f876 int32_t pd_stride[MAX_PD];  // stride to move to the next index at this level #endif // MAX_PD > 0 int32_t pd_prod;            // total number of voxels in hypercube int32_t pd_sum;             // total length of the weights vector int32_t num_eval;           // total number of voxels in hypercube int32_t num_weights;        // total length of the weights vector int32_t num_active;         // number of non-trivial pd loops int32_t theta_par;          // id of spherical correction variable // Storage for the current parameter values.  These will be updated as we // walk the polydispersity cube.  local_values will be aliased to pvec. // walk the polydispersity cube. ParameterBlock local_values; // Fill in the initial variables for (int i=0; i < NUM_PARS; i++) { local_values.vector[i] = values[2+i]; //if (q_index==0) printf("p%d = %g\n",i, local_values.vector[i]); } #if defined(MAGNETIC) && NUM_MAGNETIC>0 #endif // MAGNETIC double pd_norm, this_result; if (pd_start == 0) { pd_norm = this_result = 0.0; } else { pd_norm = result[nq]; this_result = result[q_index]; } // Fill in the initial variables //   values[0] is scale //   values[1] is background for (int i=0; i < NUM_PARS; i++) { local_values.vector[i] = values[2+i]; //if (q_index==0) printf("p%d = %g\n",i, local_values.vector[i]); } //if (q_index==0) printf("NUM_VALUES:%d  NUM_PARS:%d  MAX_PD:%d\n", NUM_VALUES, NUM_PARS, MAX_PD); //if (q_index==0) printf("start:%d stop:%d\n", pd_start, pd_stop); double pd_norm = (pd_start == 0 ? 0.0 : result[nq]); double this_result = (pd_start == 0 ? 0.0 : result[q_index]); //if (q_index==0) printf("start %d %g %g\n", pd_start, pd_norm, this_result); #if MAX_PD>0 global const double *pd_value = values + NUM_VALUES + 2; global const double *pd_weight = pd_value + details->pd_sum; global const double *pd_value = values + NUM_VALUES; global const double *pd_weight = pd_value + details->num_weights; #endif // TODO: what is the magnetic scattering at q=0 if (qsq > 1.e-16) { double p[4];  // spin_i, spin_f double p[4];  // dd, du, ud, uu p[0] = (qy*cos_mspin + qx*sin_mspin)/qsq; p[3] = -p[0];
• ## sasmodels/kernelcl.py

 r40a87fa ] #print("Calling OpenCL") #print("values",values) #call_details.show(values) # Call kernel and retrieve results last_call = None step = 100 for start in range(0, call_details.pd_prod, step): stop = min(start+step, call_details.pd_prod) for start in range(0, call_details.num_eval, step): stop = min(start + step, call_details.num_eval) #print("queuing",start,stop) args[1:3] = [np.int32(start), np.int32(stop)] None, *args, wait_for=last_call)] cl.enqueue_copy(self.queue, self.result, self.result_b) #print("result", self.result) # Free buffers scale = values[0]/self.result[self.q_input.nq] background = values[1] #print("scale",scale,background) #print("scale",scale,values[0],self.result[self.q_input.nq]) return scale*self.result[:self.q_input.nq] + background # return self.result[:self.q_input.nq]
• ## sasmodels/kerneldll.py

 r40a87fa self.real(cutoff), # cutoff ] #print("kerneldll values", values) #print("Calling DLL") #call_details.show(values) step = 100 for start in range(0, call_details.pd_prod, step): stop = min(start+step, call_details.pd_prod) for start in range(0, call_details.num_eval, step): stop = min(start + step, call_details.num_eval) args[1:3] = [start, stop] #print("calling DLL") kernel(*args) # type: ignore
• ## sasmodels/kernelpy.py

 r40a87fa """ def __init__(self, kernel, model_info, q_input): # type: (callable, ModelInfo, List[np.ndarray]) -> None self.dtype = np.dtype('d') self.info = model_info def __call__(self, call_details, values, cutoff, magnetic): assert isinstance(call_details, details.CallDetails) # type: (CallDetails, np.ndarray, np.ndarray, float, bool) -> np.ndarray if magnetic: raise NotImplementedError("Magnetism not implemented for pure python models") #print("Calling python kernel") #call_details.show(values) res = _loops(self._parameter_vector, self._form, self._volume, self.q_input.nq, call_details, values, cutoff) def release(self): # type: () -> None """ Free resources associated with the kernel. return np.ones(nq, 'd')*background pd_value = values[2+n_pars:2+n_pars + call_details.pd_sum] pd_weight = values[2+n_pars + call_details.pd_sum:] pd_value = values[2+n_pars:2+n_pars + call_details.num_weights] pd_weight = values[2+n_pars + call_details.num_weights:] pd_norm = 0.0 total = np.zeros(nq, 'd') for loop_index in range(call_details.pd_prod): for loop_index in range(call_details.num_eval): # update polydispersity parameter values if p0_index == p0_length:
• ## sasmodels/mixture.py

 r7ae2b7f *ProductModel(P, S)*. """ from __future__ import print_function from copy import copy import numpy as np  # type: ignore from .modelinfo import Parameter, ParameterTable, ModelInfo from .kernel import KernelModel, Kernel from .details import make_details try: from typing import List from .details import CallDetails except ImportError: pass # type: (List[ModelInfo]) -> ModelInfo """ Create info block for product model. Create info block for mixture model. """ flatten = [] # Build new parameter list combined_pars = [] demo = {} for k, part in enumerate(parts): # Parameter prefix per model, A_, B_, ... # to support vector parameters prefix = chr(ord('A')+k) + '_' combined_pars.append(Parameter(prefix+'scale')) scale =  Parameter(prefix+'scale', default=1.0, description="model intensity for " + part.name) combined_pars.append(scale) for p in part.parameters.kernel_parameters: p = copy(p) p.length_control = prefix + p.length_control combined_pars.append(p) demo.update((prefix+k, v) for k, v in part.demo.items() if k != "background") #print("pars",combined_pars) parameters = ParameterTable(combined_pars) parameters.max_pd = sum(part.parameters.max_pd for part in parts) model_info = ModelInfo() # Remember the component info blocks so we can build the model model_info.composition = ('mixture', parts) model_info.demo = demo return model_info self.parts = parts def __call__(self, q_vectors): def make_kernel(self, q_vectors): # type: (List[np.ndarray]) -> MixtureKernel # Note: may be sending the q_vectors to the n times even though they self.info =  model_info self.kernels = kernels def __call__(self, call_details, value, weight, cutoff): # type: (CallDetails, np.ndarray, np.ndarry, float) -> np.ndarray scale, background = value[0:2] self.dtype = self.kernels[0].dtype def __call__(self, call_details, values, cutoff, magnetic): # type: (CallDetails, np.ndarray, np.ndarry, float, bool) -> np.ndarray scale, background = values[0:2] total = 0.0 # remember the parts for plotting later self.results = [] for kernel, kernel_details in zip(self.kernels, call_details.parts): part_result = kernel(kernel_details, value, weight, cutoff) total += part_result self.results.append(part_result) offset = 2 # skip scale & background parts = MixtureParts(self.info, self.kernels, call_details, values) for kernel, kernel_details, kernel_values in parts: #print("calling kernel", kernel.info.name) result = kernel(kernel_details, kernel_values, cutoff, magnetic) #print(kernel.info.name, result) total += result self.results.append(result) return scale*total + background k.release() class MixtureParts(object): def __init__(self, model_info, kernels, call_details, values): # type: (ModelInfo, List[Kernel], CallDetails, np.ndarray) -> None self.model_info = model_info self.parts = model_info.composition[1] self.kernels = kernels self.call_details = call_details self.values = values self.spin_index = model_info.parameters.npars + 2 #call_details.show(values) def __iter__(self): # type: () -> PartIterable self.part_num = 0 self.par_index = 2 self.mag_index = self.spin_index + 3 return self def next(self): # type: () -> Tuple[List[Callable], CallDetails, np.ndarray] if self.part_num >= len(self.parts): raise StopIteration() info = self.parts[self.part_num] kernel = self.kernels[self.part_num] call_details = self._part_details(info, self.par_index) values = self._part_values(info, self.par_index, self.mag_index) values = values.astype(kernel.dtype) #call_details.show(values) self.part_num += 1 self.par_index += info.parameters.npars + 1 self.mag_index += 3 * len(info.parameters.magnetism_index) return kernel, call_details, values def _part_details(self, info, par_index): # type: (ModelInfo, int) -> CallDetails full = self.call_details # par_index is index into values array of the current parameter, # which includes the initial scale and background parameters. # We want the index into the weight length/offset for each parameter. # Exclude the initial scale and background, so subtract two, but each # component has its own scale factor which we need to skip when # constructing the details for the kernel, so add one, giving a # net subtract one. index = slice(par_index - 1, par_index - 1 + info.parameters.npars) length = full.length[index] offset = full.offset[index] # The complete weight vector is being sent to each part so that # offsets don't need to be adjusted. part = make_details(info, length, offset, full.num_weights) return part def _part_values(self, info, par_index, mag_index): # type: (ModelInfo, int, int) -> np.ndarray #print(info.name, par_index, self.values[par_index:par_index + info.parameters.npars + 1]) scale = self.values[par_index] pars = self.values[par_index + 1:par_index + info.parameters.npars + 1] nmagnetic = len(info.parameters.magnetism_index) if nmagnetic: spin_state = self.values[self.spin_index:self.spin_index + 3] mag_index = self.values[mag_index:mag_index + 3 * nmagnetic] else: spin_state = [] mag_index = [] nvalues = self.model_info.parameters.nvalues nweights = self.call_details.num_weights weights = self.values[nvalues:nvalues+2*nweights] zero = self.values.dtype.type(0.) values = [[scale, zero], pars, spin_state, mag_index, weights] # Pad value array to a 32 value boundary spacer = (32 - sum(len(v) for v in values)%32)%32 values.append([zero]*spacer) values = np.hstack(values).astype(self.kernels[0].dtype) return values
• ## sasmodels/sasview_model.py

 rfd19811 from . import weights from . import modelinfo from .details import build_details, dispersion_mesh from .details import make_kernel_args, dispersion_mesh try: parameters = self._model_info.parameters pairs = [self._get_weights(p) for p in parameters.call_parameters] call_details, values, is_magnetic = build_details(calculator, pairs) call_details, values, is_magnetic = make_kernel_args(calculator, pairs) #call_details.show() #print("pairs", pairs)
Note: See TracChangeset for help on using the changeset viewer.