Changeset 839283a in sasmodels

Timestamp:

Mar 21, 2016 10:35:21 AM (8 years ago)

Author:

wojciech

Branches:

master, core_shell_microgels, costrafo411, magnetic_model, release_v0.94, release_v0.95, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests

Children:

1edf610

Parents:

0880966 (diff), c072f83 (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 ​https://github.com/SasView/sasmodels into polydisp

Location:

sasmodels

Files:

• generate.py (modified) (2 diffs)
• kernel_iq.c (modified) (9 diffs)
• kernelcl.py (modified) (6 diffs)
• core.py (modified) (3 diffs)
• kerneldll.py (modified) (1 diff)

sasmodels/generate.py

@@ -669 +669 @@
     c = p + 3*npars
 
-    details = np.zeros(c + 3, 'int32')
+    details = np.zeros(c + 2, 'int32')
     details[0*max_pd:1*max_pd] = range(max_pd)       # pd_par: arbitrary order; use first
     details[1*max_pd:2*max_pd] = [1]*max_pd          # pd_length: only one element
@@ -707 +707 @@
             theta_par = -1
 
-    details = np.empty(constants_offset + 3, 'int32')
+    details = np.empty(constants_offset + 2, 'int32')
     details[0*max_pd:1*max_pd] = idx             # pd_par
     details[1*max_pd:2*max_pd] = pd_length[idx]

sasmodels/kernel_iq.c

@@ -20 +20 @@
     int32_t pd_par[MAX_PD];     // id of the nth polydispersity variable
     int32_t pd_length[MAX_PD];  // length of the nth polydispersity weight vector
-    int32_t pd_offset[MAX_PD];  // offset of pd weights in the par & weight vector
+    int32_t pd_offset[MAX_PD];  // offset of pd weights in the value & weight vector
     int32_t pd_stride[MAX_PD];  // stride to move to the next index at this level
     int32_t pd_isvol[MAX_PD];   // True if parameter is a volume weighting parameter
-    int32_t par_offset[NPARS];  // offset of par values in the par & weight vector
+    int32_t par_offset[NPARS];  // offset of par values in the value & weight vector
     int32_t par_coord[NPARS];   // polydispersity coordination bitvector
     int32_t fast_coord_pars[NPARS]; // ids of the fast coordination parameters
     int32_t fast_coord_count;   // number of parameters coordinated with pd 1
-    int32_t theta_var;          // id of spherical correction variable
+    int32_t theta_par;          // id of spherical correction variable
 } ProblemDetails;
 
@@ -43 +43 @@
     global const ProblemDetails *problem,
     global const double *weights,
-    global const double *pars,
+    global const double *values,
     global const double *q, // nq q values, with padding to boundary
     global double *result,  // nq+3 return values, again with padding
@@ -51 +51 @@
   // Storage for the current parameter values.  These will be updated as we
   // walk the polydispersity cube.
-  local ParameterBlock local_pars;  // current parameter values
-  double *pvec = (double *)(&local_pars);  // Alias named parameters with a vector
+  local ParameterBlock local_values;  // current parameter values
+  double *pvec = (double *)(&local_values);  // Alias named parameters with a vector
 
   local int offset[NPARS];  // NPARS excludes scale/background
@@ -61 +61 @@
 
     for (int k=0; k < NPARS; k++) {
-      pvec[k] = pars[k+2];  // skip scale and background
-    }
-
-    const double volume = CALL_VOLUME(local_pars);
+      pvec[k] = values[k+2];  // skip scale and background
+    }
+
+    const double volume = CALL_VOLUME(local_values);
     #ifdef USE_OPENMP
     #pragma omp parallel for
     #endif
     for (int i=0; i < nq; i++) {
-      const double scattering = CALL_IQ(q, i, local_pars);
-      result[i] = pars[0]*scattering/volume + pars[1];
+      const double scattering = CALL_IQ(q, i, local_values);
+      result[i] = values[0]*scattering/volume + values[1];
     }
     return;
@@ -147 +147 @@
         }
         offset[k] = this_offset;
-        pvec[k] = pars[this_offset];
+        pvec[k] = values[this_offset];
       }
       weight = partial_weight * weights[problem->pd_offset[0]+pd_index[0]];
-      if (problem->theta_var >= 0) {
-        spherical_correction = fabs(cos(M_PI_180*pvec[problem->theta_var]));
-      }
-      if (problem->theta_var == problem->pd_par[0]) {
+      if (problem->theta_par >= 0) {
+        spherical_correction = fabs(cos(M_PI_180*pvec[problem->theta_par]));
+      }
+      if (problem->theta_par == problem->pd_par[0]) {
         weight *= spherical_correction;
       }
@@ -166 +166 @@
       for (int k=0; k < problem->fast_coord_count; k++) {
         pvec[problem->fast_coord_pars[k]]
-            = pars[offset[problem->fast_coord_pars[k]]++];
-      }
-      if (problem->theta_var ==problem->pd_par[0]) {
-        weight *= fabs(cos(M_PI_180*pvec[problem->theta_var]));
+            = values[offset[problem->fast_coord_pars[k]]++];
+      }
+      if (problem->theta_par ==problem->pd_par[0]) {
+        weight *= fabs(cos(M_PI_180*pvec[problem->theta_par]));
       }
     }
 
     #ifdef INVALID
-    if (INVALID(local_pars)) continue;
+    if (INVALID(local_values)) continue;
     #endif
 
@@ -181 +181 @@
     if (weight > cutoff) {
       norm += weight;
-      vol += vol_weight * CALL_VOLUME(local_pars);
+      vol += vol_weight * CALL_VOLUME(local_values);
       norm_vol += vol_weight;
 
@@ -188 +188 @@
       #endif
       for (int i=0; i < nq; i++) {
-        const double scattering = CALL_IQ(q, i, local_pars);
+        const double scattering = CALL_IQ(q, i, local_values);
         result[i] += weight*scattering;
       }
@@ -207 +207 @@
         result[i] *= norm_vol/vol;
       }
-      result[i] = pars[0]*result[i]/norm + pars[1];
+      result[i] = values[0]*result[i]/norm + values[1];
     }
   }

sasmodels/kernelcl.py

@@ -336 +336 @@
         self.program = None
 
-    def __call__(self, q_vectors):
+    def make_calculator(self, q_vectors, details):
         if self.program is None:
             compiler = environment().compile_program
-            self.program = compiler(self.info['name'], self.source, self.dtype,
-                                    self.fast)
+            self.program = compiler(self.info['name'], self.source,
+                                    self.dtype, self.fast)
         is_2d = len(q_vectors) == 2
         kernel_name = generate.kernel_name(self.info, is_2d)
         kernel = getattr(self.program, kernel_name)
-        return GpuKernel(kernel, self.info, q_vectors, self.dtype)
+        return GpuKernel(kernel, self.info, q_vectors, details, self.dtype)
 
     def release(self):
@@ -387 +387 @@
         # at this point, so instead using 32, which is good on the set of
         # architectures tested so far.
-        self.q_vectors = [_stretch_input(q, self.dtype, 32) for q in q_vectors]
+        if self.is_2d:
+            # Note: 17 rather than 15 because results is 2 elements
+            # longer than input.
+            width = ((self.nq+17)//16)*16
+            self.q = np.empty((width, 2), dtype=dtype)
+            self.q[:self.nq, 0] = q_vectors[0]
+            self.q[:self.nq, 1] = q_vectors[1]
+        else:
+            # Note: 33 rather than 31 because results is 2 elements
+            # longer than input.
+            width = ((self.nq+33)//32)*32
+            self.q = np.empty(width, dtype=dtype)
+            self.q[:self.nq] = q_vectors[0]
+        self.global_size = [self.q.shape[0]]
         context = env.get_context(self.dtype)
-        self.global_size = [self.q_vectors[0].size]
         #print("creating inputs of size", self.global_size)
-        self.q_buffers = [
-            cl.Buffer(context, mf.READ_ONLY | mf.COPY_HOST_PTR, hostbuf=q)
-            for q in self.q_vectors
-        ]
+        # COPY_HOST_PTR initiates transfer as necessary?
+        self.q_b = cl.Buffer(context, mf.READ_ONLY | mf.COPY_HOST_PTR,
+                             hostbuf=self.q)
 
     def release(self):
@@ -400 +411 @@
         Free the memory.
         """
-        for b in self.q_buffers:
-            b.release()
-        self.q_buffers = []
+        if self.q is not None:
+            self.q.release()
+            self.q = None
 
     def __del__(self):
@@ -427 +438 @@
     Call :meth:`release` when done with the kernel instance.
     """
-    def __init__(self, kernel, model_info, q_vectors, dtype):
+    def __init__(self, kernel, model_info, q_vectors, details, dtype):
+        if details.dtype != np.int32:
+            raise TypeError("numeric type does not match the kernel type")
+
+        max_pd = self.info['max_pd']
+        npars = len(model_info['parameters'])-2
         q_input = GpuInput(q_vectors, dtype)
+        self.dtype = dtype
         self.kernel = kernel
         self.info = model_info
-        self.res = np.empty(q_input.nq, q_input.dtype)
-        dim = '2d' if q_input.is_2d else '1d'
-        self.fixed_pars = model_info['partype']['fixed-' + dim]
-        self.pd_pars = model_info['partype']['pd-' + dim]
+        self.details = details
+        self.pd_stop_index = 4*max_pd-1
+        # plus three for the normalization values
+        self.result = np.empty(q_input.nq+3, q_input.dtype)
+        #self.dim = '2d' if q_input.is_2d else '1d'
 
         # Inputs and outputs for each kernel call
@@ -440 +458 @@
         env = environment()
         self.queue = env.get_queue(dtype)
-        self.loops_b = cl.Buffer(self.queue.context, mf.READ_WRITE,
-                                 2 * MAX_LOOPS * q_input.dtype.itemsize)
-        self.res_b = cl.Buffer(self.queue.context, mf.READ_WRITE,
+
+        # details is int32 data, padded to a 32 integer boundary
+        size = 4*((self.info['mono'].size+7)//8)*8 # padded to 32 byte boundary
+        self.details_b = cl.Buffer(self.queue.context,
+                                   mf.READ_ONLY | mf.COPY_HOST_PTR,
+                                   hostbuf=details)
+        size = np.sum(details[max_pd:2*max_pd])
+        self.weights_b = cl.Buffer(self.queue.context, mf.READ_ONLY, size)
+        size = np.sum(details[max_pd:2*max_pd])+npars
+        self.values_b = cl.Buffer(self.queue.context, mf.READ_ONLY, size)
+        self.result_b = cl.Buffer(self.queue.context, mf.READ_WRITE,
                                q_input.global_size[0] * q_input.dtype.itemsize)
-        self.q_input = q_input
-
-        self._need_release = [self.loops_b, self.res_b, self.q_input]
-
-    def __call__(self, details, weights, values, cutoff):
+        self.q_input = q_input # allocated by GpuInput above
+
+        self._need_release = [
+            self.details_b, self.weights_b, self.values_b, self.result_b,
+            self.q_input,
+        ]
+
+    def __call__(self, weights, values, cutoff):
         real = (np.float32 if self.q_input.dtype == generate.F32
                 else np.float64 if self.q_input.dtype == generate.F64
@@ -454 +483 @@
                 else np.float32)  # will never get here, so use np.float32
 
-        #print "pars", fixed_pars, pd_pars
-        res_bi = self.res_b
-        nq = np.uint32(self.q_input.nq)
-        if pd_pars:
-            cutoff = real(cutoff)
-            loops_N = [np.uint32(len(p[0])) for p in pd_pars]
-            loops = np.hstack(pd_pars) \
-                if pd_pars else np.empty(0, dtype=self.q_input.dtype)
-            loops = np.ascontiguousarray(loops.T, self.q_input.dtype).flatten()
-            #print("loops",Nloops, loops)
-
-            #import sys; print("opencl eval",pars)
-            #print("opencl eval",pars)
-            if len(loops) > 2 * MAX_LOOPS:
-                raise ValueError("too many polydispersity points")
-
-            loops_bi = self.loops_b
-            cl.enqueue_copy(self.queue, loops_bi, loops)
-            loops_l = cl.LocalMemory(len(loops.data))
-            #ctx = environment().context
-            #loops_bi = cl.Buffer(ctx, mf.READ_ONLY|mf.COPY_HOST_PTR, hostbuf=loops)
-            dispersed = [loops_bi, loops_l, cutoff] + loops_N
-        else:
-            dispersed = []
-        fixed = [real(p) for p in fixed_pars]
-        args = self.q_input.q_buffers + [res_bi, nq] + dispersed + fixed
+        if weights.dtype != real or values.dtype != real:
+            raise TypeError("numeric type does not match the kernel type")
+
+        cl.enqueue_copy(self.queue, self.weights_b, weights)
+        cl.enqueue_copy(self.queue, self.values_b, values)
+
+        args = [
+            np.uint32(self.q_input.nq),
+            np.uint32(0),
+            np.uint32(self.details[self.pd_stop_index]),
+            self.details_b,
+            self.weights_b,
+            self.values_b,
+            self.q_input.q_b,
+            self.result_b,
+            real(cutoff),
+        ]
         self.kernel(self.queue, self.q_input.global_size, None, *args)
-        cl.enqueue_copy(self.queue, self.res, res_bi)
-
-        return self.res
+        cl.enqueue_copy(self.queue, self.result, self.result_b)
+
+        return self.result[:self.nq]
 
     def release(self):

sasmodels/core.py

@@ -183 +183 @@
     value = values.get(parameter.name, parameter.default)
     if parameter.type not in ('volume', 'orientation'):
-        return [value], []
+        return np.array([value]), np.array([1.0])
     relative = parameter.type == 'volume'
     limits = parameter.limits
@@ -208 +208 @@
     return value, weight
 
-def call_kernel(kernel, values, cutoff=0, mono=False):
+def call_kernel(kernel, pars, cutoff=0, mono=False):
     """
     Call *kernel* returned from :func:`make_kernel` with parameters *pars*.
@@ -221 +221 @@
     *mono* is True if polydispersity should be set to none on all parameters.
     """
-    if mono or True:
-        pars = np.array([values.get(p.name, p.default)
-                         for p in kernel.info['parameters']],
-                        dtype=kernel.dtype)
-        weights = np.array([1.0], dtype=kernel.dtype)
+    if mono:
+        values = [pars.get(p.name, p.default) for p in kernel.info['parameters']]
+        weights = [1.0]*len(values)
+    else:
+        wv_pairs = [get_weights(p, pars) for p in kernel.info['parameters']]
+        weights, values = [v for v in zip(*wv_pairs)]
+
+    #TODO: This is what we thought to do if max([len(w) for w in weights]) > 1:
+    if max([w for w in weights]) > 1:
+        details = generate.poly_details(kernel.info, weights)
+    else:
         details = kernel.info['mono_details']
-        return kernel(details,  weights, pars, cutoff)
-    else:
-        pairs = [get_weights(p, values) for p in kernel.info['parameters']]
-        weights, pars = [v for v in zip(*pairs)]
-        details = generate.poly_details(kernel.info, weights, pars)
-        weights, pars = [np.hstack(v) for v in (weights, pars)]
-        return kernel(details, weights, pars, cutoff)
+
+    weights, values = [np.hstack(v) for v in (weights, values)]
+
+    weights = weights.astype(dtype=kernel.dtype)
+    values = values.astype(dtype=kernel.dtype)
+
+    return kernel(details, weights, values, cutoff)
 
 def call_ER_VR(model_info, vol_pars):

sasmodels/kerneldll.py

@@ -271 +271 @@
         #weights = np.asarray(weights, dtype=real)
         #values = np.asarray(values, dtype=real)
+        #TODO: How can I access max_pd and is this the way to do it?
+        #max_pd = model_info['max_pd']
+        max_pd = 1
         args = [
             self.q_input.nq, # nq
+            #TODO: pd_start will need to be changed
             0, # pd_start
-            1, # pd_stop
+            details[3*max_pd:4*max_pd], # pd_stop pd_stride[MAX_PD]
             details.ctypes.data, # problem
             weights.ctypes.data,  # weights