Changeset 6aee3ab in sasmodels


Ignore:
Timestamp:
Nov 2, 2017 3:24:41 PM (6 years ago)
Author:
Paul Kienzle <pkienzle@…>
Branches:
master, core_shell_microgels, magnetic_model, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests
Children:
ce99754
Parents:
9d9fcbd
Message:

make sure that jitter angle defaults to zero if no polydispersity

Location:
sasmodels
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • sasmodels/details.py

    r767dca8 r6aee3ab  
    193193    #print("off:",offset) 
    194194 
    195     # Check that we arn't using too many polydispersity loops 
     195    # Check that we aren't using too many polydispersity loops 
    196196    num_active = np.sum(length > 1) 
    197197    max_pd = model_info.parameters.max_pd 
  • sasmodels/kernel_iq.c

    rff10479 r6aee3ab  
    1313// NOTE: the following macros are defined in generate.py: 
    1414// 
    15 //  MAX_PD : the maximum number of dispersity loops allowed 
     15//  MAX_PD : the maximum number of dispersity loops allowed for this model, 
     16//      which will be at most modelinfo.MAX_PD. 
    1617//  NUM_PARS : the number of parameters in the parameter table 
    1718//  NUM_VALUES : the number of values to skip at the start of the 
     
    282283#endif 
    283284 
    284   // Storage for the current parameter values.  These will be updated as we 
    285   // walk the dispersity mesh. 
     285  // ** Fill in the local values table ** 
     286  // Storage for the current parameter values. 
     287  // These will be updated as we walk the dispersity mesh. 
    286288  ParameterBlock local_values; 
    287  
     289  //   values[0] is scale 
     290  //   values[1] is background 
     291  #ifdef USE_OPENMP 
     292  #pragma omp parallel for 
     293  #endif 
     294  for (int i=0; i < NUM_PARS; i++) { 
     295    local_values.vector[i] = values[2+i]; 
     296    //if (q_index==0) printf("p%d = %g\n",i, local_values.vector[i]); 
     297  } 
     298  //if (q_index==0) printf("NUM_VALUES:%d  NUM_PARS:%d  MAX_PD:%d\n", NUM_VALUES, NUM_PARS, MAX_PD); 
     299  //if (q_index==0) printf("start:%d  stop:%d\n", pd_start, pd_stop); 
     300 
     301  // ** Precompute magnatism values ** 
    288302#if defined(MAGNETIC) && NUM_MAGNETIC>0 
    289303  // Location of the sld parameters in the parameter vector. 
     
    302316#endif // MAGNETIC 
    303317 
    304   // Fill in the initial variables 
    305   //   values[0] is scale 
    306   //   values[1] is background 
    307   #ifdef USE_OPENMP 
    308   #pragma omp parallel for 
    309   #endif 
    310   for (int i=0; i < NUM_PARS; i++) { 
    311     local_values.vector[i] = values[2+i]; 
    312 //if (q_index==0) printf("p%d = %g\n",i, local_values.vector[i]); 
    313   } 
    314 //if (q_index==0) printf("NUM_VALUES:%d  NUM_PARS:%d  MAX_PD:%d\n", NUM_VALUES, NUM_PARS, MAX_PD); 
    315 //if (q_index==0) printf("start:%d  stop:%d\n", pd_start, pd_stop); 
    316  
     318  // ** Fill in the initial results ** 
    317319  // If pd_start is zero that means that we are starting a new calculation, 
    318320  // and must initialize the result to zero.  Otherwise, we are restarting 
     
    336338      for (int q_index=0; q_index < nq; q_index++) result[q_index] = 0.0; 
    337339    } 
    338 //if (q_index==0) printf("start %d %g %g\n", pd_start, pd_norm, result[0]); 
     340    //if (q_index==0) printf("start %d %g %g\n", pd_start, pd_norm, result[0]); 
    339341#endif // !USE_OPENCL 
    340342 
     
    390392      PD_OPEN(0,1) 
    391393 
    392       ... main loop body ... 
     394      // ... main loop body ... 
     395      APPLY_PROJECTION    // convert jitter values to spherical coords 
     396      BUILD_ROTATION      // construct the rotation matrix qxy => qabc 
     397      for each q 
     398          FETCH_Q         // set qx,qy from the q input vector 
     399          APPLY_ROTATION  // convert qx,qy to qa,qb,qc 
     400          CALL_KERNEL     // scattering = Iqxy(qa, qb, qc, p1, p2, ...) 
    393401 
    394402      ++step;  // increment counter representing position in dispersity mesh 
     
    412420*/ 
    413421 
    414 // Define looping variables 
    415 #define PD_INIT(_LOOP) \ 
    416   const int n##_LOOP = details->pd_length[_LOOP]; \ 
    417   const int p##_LOOP = details->pd_par[_LOOP]; \ 
    418   global const double *v##_LOOP = pd_value + details->pd_offset[_LOOP]; \ 
    419   global const double *w##_LOOP = pd_weight + details->pd_offset[_LOOP]; \ 
    420   int i##_LOOP = (pd_start/details->pd_stride[_LOOP])%n##_LOOP; 
    421  
    422 // Jump into the middle of the dispersity loop 
    423 #define PD_OPEN(_LOOP,_OUTER) \ 
    424   while (i##_LOOP < n##_LOOP) { \ 
    425     local_values.vector[p##_LOOP] = v##_LOOP[i##_LOOP]; \ 
    426     const double weight##_LOOP = w##_LOOP[i##_LOOP] * weight##_OUTER; 
    427  
    428 // create the variable "weight#=1.0" where # is the outermost level+1 (=MAX_PD). 
    429 #define _PD_OUTERMOST_WEIGHT(_n) const double weight##_n = 1.0; 
    430 #define PD_OUTERMOST_WEIGHT(_n) _PD_OUTERMOST_WEIGHT(_n) 
    431  
    432 // Close out the loop 
    433 #define PD_CLOSE(_LOOP) \ 
    434     if (step >= pd_stop) break; \ 
    435     ++i##_LOOP; \ 
    436   } \ 
    437   i##_LOOP = 0; 
    438  
    439 // The main loop body is essentially: 
    440 // 
    441 //    BUILD_ROTATION      // construct the rotation matrix qxy => qabc 
    442 //    for each q 
    443 //        FETCH_Q         // set qx,qy from the q input vector 
    444 //        APPLY_ROTATION  // convert qx,qy to qa,qb,qc 
    445 //        CALL_KERNEL     // scattering = Iqxy(qa, qb, qc, p1, p2, ...) 
    446 // 
     422 
     423// ** prepare inner loops ** 
     424 
    447425// Depending on the shape type (radial, axial, triaxial), the variables 
    448 // and calling parameters will be slightly different.  These macros 
    449 // capture the differences in one spot so the rest of the code is easier 
    450 // to read. 
     426// and calling parameters in the loop body will be slightly different. 
     427// Macros capture the differences in one spot so the rest of the code 
     428// is easier to read. The code below both declares variables for the 
     429// inner loop and defines the macros that use them. 
    451430 
    452431#if defined(CALL_IQ) 
     
    486465  // psi and dpsi are only for IQ_ABC, so they are processed here. 
    487466  const double psi = values[details->theta_par+4]; 
     467  local_values.table.psi = 0.; 
    488468  #define BUILD_ROTATION() qabc_rotation(&rotation, theta, phi, psi, dtheta, dphi, local_values.table.psi) 
    489469  #define APPLY_ROTATION() qabc_apply(rotation, qx, qy, &qa, &qb, &qc) 
     
    501481  const double theta = values[details->theta_par+2]; 
    502482  const double phi = values[details->theta_par+3]; 
     483  // Make sure jitter angle defaults to zero if there is no jitter distribution 
     484  local_values.table.theta = 0.; 
     485  local_values.table.phi = 0.; 
    503486  // The "jitter" angles (dtheta, dphi, dpsi) are stored with the 
    504487  // dispersity values and copied to the local parameter table as 
     
    523506#endif // !spherosymmetric projection 
    524507 
     508// ** define looping macros ** 
     509 
     510// Define looping variables 
     511#define PD_INIT(_LOOP) \ 
     512  const int n##_LOOP = details->pd_length[_LOOP]; \ 
     513  const int p##_LOOP = details->pd_par[_LOOP]; \ 
     514  global const double *v##_LOOP = pd_value + details->pd_offset[_LOOP]; \ 
     515  global const double *w##_LOOP = pd_weight + details->pd_offset[_LOOP]; \ 
     516  int i##_LOOP = (pd_start/details->pd_stride[_LOOP])%n##_LOOP; 
     517 
     518// Jump into the middle of the dispersity loop 
     519#define PD_OPEN(_LOOP,_OUTER) \ 
     520  while (i##_LOOP < n##_LOOP) { \ 
     521    local_values.vector[p##_LOOP] = v##_LOOP[i##_LOOP]; \ 
     522    const double weight##_LOOP = w##_LOOP[i##_LOOP] * weight##_OUTER; 
     523 
     524// create the variable "weight#=1.0" where # is the outermost level+1 (=MAX_PD). 
     525#define _PD_OUTERMOST_WEIGHT(_n) const double weight##_n = 1.0; 
     526#define PD_OUTERMOST_WEIGHT(_n) _PD_OUTERMOST_WEIGHT(_n) 
     527 
     528// Close out the loop 
     529#define PD_CLOSE(_LOOP) \ 
     530    if (step >= pd_stop) break; \ 
     531    ++i##_LOOP; \ 
     532  } \ 
     533  i##_LOOP = 0; 
    525534 
    526535// ====== construct the loops ======= 
     
    537546int step = pd_start; 
    538547 
     548// *** define loops for each of 0, 1, 2, ..., modelinfo.MAX_PD-1 *** 
     549 
    539550// define looping variables 
    540551#if MAX_PD>4 
     
    571582  PD_OPEN(0,1) 
    572583#endif 
    573  
    574584 
    575585//if (q_index==0) {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");} 
     
    672682#endif // !USE_OPENCL 
    673683 
    674 // clear the macros in preparation for the next kernel 
     684// ** clear the macros in preparation for the next kernel ** 
    675685#undef PD_INIT 
    676686#undef PD_OPEN 
  • sasmodels/modelinfo.py

    re3571cb r6aee3ab  
    3030    TestCondition = Tuple[ParameterSetUser, TestInput, TestValue] 
    3131 
    32 MAX_PD = 4 #: Maximum number of simultaneously polydisperse parameters 
     32# If MAX_PD changes, need to change the loop macros in kernel_iq.c 
     33MAX_PD = 5 #: Maximum number of simultaneously polydisperse parameters 
    3334 
    3435# assumptions about common parameters exist throughout the code, such as: 
Note: See TracChangeset for help on using the changeset viewer.