-                      rba32cdd
+                      r5ff1b03
     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
 #endif // MAX_PD > 0
+    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 par_offset[NPARS];  // offset of par value blocks in the value & weight vector
+    int32_t par_coord[NPARS];   // ids of the coordination parameters
+    int32_t pd_coord[NPARS];    // polydispersity coordination bitvector
+    int32_t num_active;         // number of non-trivial pd loops
+    int32_t total_pd;           // total number of voxels in hypercube
+    int32_t num_coord;          // number of coordinated parameters
     int32_t theta_par;          // id of spherical correction variable
 } ProblemDetails;
 …
   double *pvec = (double *)(&local_values);  // Alias named parameters with a vector
+  // Monodisperse computation
+  if (pd_stop == 1) {
+    // Shouldn't need to copy!!
+    for (int k=0; k < NPARS; k++) {
+      pvec[k] = values[k+2];  // skip scale and background
+    }
+    const double volume = CALL_VOLUME(local_values);
+  // Fill in the initial variables
+  #ifdef USE_OPENMP
+  #pragma omp parallel for
+  #endif
+  for (int k=0; k < NPARS; k++) {
+    pvec[k] = values[problem->par_offset[k]];
+  }
+  // If it is the first round initialize the result to zero, otherwise
+  // assume that the previous result has been passed back.
+  // Note: doing this even in the monodisperse case in order to handle the
+  // rare case where the model parameters are invalid and zero is returned.
+  // So slightly increased cost for slightly smaller code size.
+  if (pd_start == 0) {
     #ifdef USE_OPENMP
     #pragma omp parallel for
     #endif
+    for (int i=0; i < nq+1; i++) {
+      result[i] = 0.0;
+    }
+  }
+  // Monodisperse computation
+  if (problem->num_active == 0) {
+    #ifdef INVALID
+    if (INVALID(local_values)) { return; }
+    #endif
+    const double norm = CALL_VOLUME(local_values);
+    #ifdef USE_OPENMP
+    #pragma omp parallel for
+    #endif
+    result[nq] = norm; // Total volume normalization
     for (int i=0; i < nq; i++) {
       double scattering = CALL_IQ(q, i, local_values);
+      if (volume != 0.0) scattering /= volume;
+      result[i] = values[0]*scattering + values[1];
+      result[i] = values[0]*scattering/norm + values[1];
+    }
     return;
 …
 #if MAX_PD > 0
+  //printf("Entering polydispersity\n");
+  //printf("Entering polydispersity from %d to %d\n", pd_start, pd_stop);
   // Since we are no longer looping over the entire polydispersity hypercube
+  // for each q, we need to track the normalization values for each q in a
+  // separate work vector.
+  double norm;   // contains sum over weights
+  double vol; // contains sum over volume
+  double norm_vol; // contains weights over volume
+  // Initialize the results to zero
+  if (pd_start == 0) {
+    norm_vol = 0.0;
+    norm = 0.0;
+    vol = 0.0;
+    #ifdef USE_OPENMP
+    #pragma omp parallel for
+    #endif
+    for (int i=0; i < nq; i++) {
+      result[i] = 0.0;
+    }
+  } else {
+    //Pulling values from previous segment
+    norm = result[nq];
+    vol = result[nq+1];
+    norm_vol = result[nq+2];
+  }
+  // Location in the polydispersity hypercube, one index per dimension.
+  local int pd_index[MAX_PD];
+  // polydispersity loop index positions
+  local int offset[NPARS];  // NPARS excludes scale/background
+  // Trigger the reset behaviour that happens at the end the fast loop
+  // by setting the initial index >= weight vector length.
+  pd_index[0] = problem->pd_length[0];
+  // for each q, we need to track the normalization values between calls.
+  double norm = 0.0;
   // need product of weights at every Iq calc, so keep product of
   // weights from the outer loops so that weight = partial_weight * fast_weight
   double partial_weight = NAN; // product of weight w4*w3*w2 but not w1
+  double partial_volweight = NAN;
+  double weight = 1.0;        // set to 1 in case there are no weights
+  double vol_weight = 1.0;    // set to 1 in case there are no vol weights
+  double spherical_correction = 1.0;  // correction for latitude variation
+  double spherical_correction = 1.0;  // cosine correction for latitude variation
+  // Location in the polydispersity hypercube, one index per dimension.
+  local int pd_index[MAX_PD];
+  // Location of the coordinated parameters in their own sub-cubes.
+  local int offset[NPARS];
+  // Trigger the reset behaviour that happens at the end the fast loop
+  // by setting the initial index >= weight vector length.
+  const int fast_length = problem->pd_length[0];
+  pd_index[0] = fast_length;
   // Loop over the weights then loop over q, accumulating values
   for (int loop_index=pd_start; loop_index < pd_stop; loop_index++) {
     // check if indices need to be updated
+    if (pd_index[0] >= problem->pd_length[0]) {
+      // RESET INDICES
+      pd_index[0] = loop_index%problem->pd_length[0];
+    if (pd_index[0] == fast_length) {
+      //printf("should be here with %d active\n", problem->num_active);
+      // Compute position in polydispersity hypercube
+      for (int k=0; k < problem->num_active; k++) {
+        pd_index[k] = (loop_index/problem->pd_stride[k])%problem->pd_length[k];
+        //printf("pd_index[%d] = %d\n",k,pd_index[k]);
+      }
+      // Compute partial weights
       partial_weight = 1.0;
       partial_volweight = 1.0;
       for (int k=1; k < MAX_PD; k++) {
         pd_index[k] = (loop_index%problem->pd_length[k])/problem->pd_stride[k];
         const double wi = weights[problem->pd_offset[k]+pd_index[k]];
+      //printf("partial weight %d: ", loop_index);
+      for (int k=1; k < problem->num_active; k++) {
+        double wi = weights[problem->pd_offset[k] + pd_index[k]];
+        //printf("pd[%d]=par[%d]=%g ", k, problem->pd_par[k], wi);
         partial_weight *= wi;
+        if (problem->pd_isvol[k]) partial_volweight *= wi;
+      }
+      }
+      //printf("\n");
+      // Update parameter offsets in weight vector
       //printf("slow %d: ", loop_index);
+      for (int k=0; k < NPARS; k++) {
+        int coord = problem->par_coord[k];
+        int this_offset = problem->par_offset[k];
+      for (int k=0; k < problem->num_coord; k++) {
+        int par = problem->par_coord[k];
+        int coord = problem->pd_coord[k];
+        int this_offset = problem->par_offset[par];
         int block_size = 1;
         for (int bit=0; bit < MAX_PD && coord != 0; bit++) {
+        for (int bit=0; coord != 0; bit++) {
           if (coord&1) {
               this_offset += block_size * pd_index[bit];
               block_size *= problem->pd_length[bit];
+          }
           coord /= 2;
+          coord >>= 1;
+        }
+        offset[k] = this_offset;
+        pvec[k] = values[this_offset];
+        //printf("p[%d]=v[%d]=%g ", k, offset[k], pvec[k]);
+        offset[par] = this_offset;
+        pvec[par] = values[this_offset];
+        //printf("par[%d]=v[%d]=%g \n", k, offset[k], pvec[k]);
+        // if theta is not coordinated with fast index, precompute spherical correction
+        if (par == problem->theta_par && !(problem->par_coord[k]&1)) {
+          spherical_correction = fmax(fabs(cos(M_PI_180*pvec[problem->theta_par])), 1e-6);
+        }
+      }
       //printf("\n");
+      weight = partial_weight * weights[problem->pd_offset[0]+pd_index[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;
+      }
+      pd_index[0] += 1;
+    } else {
+      // INCREMENT INDICES
+      const double wi = weights[problem->pd_offset[0]+pd_index[0]];
+      weight = partial_weight*wi;
+      if (problem->pd_isvol[0]) vol_weight *= wi;
+      //printf("fast %d: ", loop_index);
+      for (int k=0; k < problem->fast_coord_count; k++) {
+        const int pindex = problem->fast_coord_pars[k];
+        pvec[pindex] = values[++offset[pindex]];
+        //printf("p[%d]=v[%d]=%g ", pindex, offset[pindex], pvec[pindex]);
+      }
+      //printf("\n");
+      if (problem->theta_par == problem->pd_par[0]) {
+        weight *= fabs(cos(M_PI_180*pvec[problem->theta_par]));
+      }
+      pd_index[0] += 1;
+    }
+    }
+    // Increment fast index
+    const double wi = weights[problem->pd_offset[0] + pd_index[0]++];
+    double weight = partial_weight*wi;
+    //printf("fast %d: ", loop_index);
+    for (int k=0; k < problem->num_coord; k++) {
+      if (problem->pd_coord[k]&1) {
+        const int par = problem->par_coord[k];
+        pvec[par] = values[offset[par]++];
+        //printf("p[%d]=v[%d]=%g ", par, offset[par]-1, pvec[par]);
+        // if theta is coordinated with fast index, compute spherical correction each time
+        if (par == problem->theta_par) {
+          spherical_correction = fmax(fabs(cos(M_PI_180*pvec[problem->theta_par])), 1e-6);
+        }
+      }
+    }
+    //printf("\n");
     #ifdef INVALID
     if (INVALID(local_values)) continue;
 …
     // Note: weight==0 must always be excluded
     if (weight > cutoff) {
+      norm += weight;
+      vol += vol_weight * CALL_VOLUME(local_values);
+      norm_vol += vol_weight;
+      // spherical correction has some nasty effects when theta is +90 or -90
+      // where it becomes zero.  If the entirety of the correction
+      weight *= spherical_correction;
+      norm += weight * CALL_VOLUME(local_values);
       #ifdef USE_OPENMP
 …
   // Make normalization available for the next round
+  result[nq] = norm;
+  result[nq+1] = vol;
+  result[nq+2] = norm_vol;
+  result[nq] += norm;
   // End of the PD loop we can normalize
   if (pd_stop >= problem->pd_stride[MAX_PD-1]) {
+  if (pd_stop >= problem->total_pd) {
     #ifdef USE_OPENMP
     #pragma omp parallel for
     #endif
     for (int i=0; i < nq; i++) {
-      if (vol*norm_vol != 0.0) {
-        result[i] *= norm_vol/vol;
+      }
       result[i] = values[0]*result[i]/norm + values[1];
+    }

Note: See TracChangeset for help on using the changeset viewer.

SasView

Changeset 5ff1b03 in sasmodels for sasmodels/kernel_iq.c

Legend:

sasmodels/kernel_iq.c

Download in other formats: