Changeset 5ff1b03 in sasmodels


Ignore:
Timestamp:
Mar 25, 2016 11:44:37 AM (9 years ago)
Author:
Paul Kienzle <pkienzle@…>
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:
c499331
Parents:
ba32cdd
Message:

working kerneldll

Location:
sasmodels
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • sasmodels/core.py

    r380e8c9 r5ff1b03  
    177177    value = values.get(parameter.name, parameter.default) 
    178178    if parameter.type not in ('volume', 'orientation'): 
    179         return [value], [1.0] 
     179        return [value], [] 
    180180    relative = parameter.type == 'volume' 
    181181    limits = parameter.limits 
    … …  
    226226        active = lambda name: True 
    227227 
    228     vw_pairs = [(get_weights(p, pars) if active(p.name) else ([p.default], [1.0])) 
     228    vw_pairs = [(get_weights(p, pars) if active(p.name) else ([p.default], [])) 
    229229                for p in kernel.info['parameters']] 
    230230    values, weights = zip(*vw_pairs) 

  • sasmodels/generate.py

    r380e8c9 r5ff1b03  
    508508        # faster by not using/transferring the volume normalizations, but 
    509509        # the ifdef's reduce readability more than is worthwhile. 
    510         call_volume = "#define CALL_VOLUME(v) 0.0" 
     510        call_volume = "#define CALL_VOLUME(v) 1.0" 
    511511    source.append(call_volume) 
    512512 
    … …  
    579579    model_info['max_pd'] = min(partable.num_pd, MAX_PD) 
    580580 
     581class CoordinationDetails(object): 
     582    def __init__(self, model_info): 
     583        max_pd = model_info['max_pd'] 
     584        npars = len(model_info['parameters'].kernel_pars()) 
     585        par_offset = 4*max_pd 
     586        self.details = np.zeros(par_offset + 3*npars + 4, 'i4') 
     587 
     588        # generate views on different parts of the array 
     589        self._pd_par     = self.details[0*max_pd:1*max_pd] 
     590        self._pd_length  = self.details[1*max_pd:2*max_pd] 
     591        self._pd_offset  = self.details[2*max_pd:3*max_pd] 
     592        self._pd_stride  = self.details[3*max_pd:4*max_pd] 
     593        self._par_offset = self.details[par_offset+0*npars:par_offset+1*npars] 
     594        self._par_coord  = self.details[par_offset+1*npars:par_offset+2*npars] 
     595        self._pd_coord   = self.details[par_offset+2*npars:par_offset+3*npars] 
     596 
     597        # theta_par is fixed 
     598        self.details[-1] = model_info['parameters'].theta_par 
     599 
     600    @property 
     601    def ctypes(self): return self.details.ctypes 
     602    @property 
     603    def pd_par(self): return self._pd_par 
     604    @property 
     605    def pd_length(self): return self._pd_length 
     606    @property 
     607    def pd_offset(self): return self._pd_offset 
     608    @property 
     609    def pd_stride(self): return self._pd_stride 
     610    @property 
     611    def pd_coord(self): return self._pd_coord 
     612    @property 
     613    def par_coord(self): return self._par_coord 
     614    @property 
     615    def par_offset(self): return self._par_offset 
     616    @property 
     617    def num_coord(self): return self.details[-2] 
     618    @num_coord.setter 
     619    def num_coord(self, v): self.details[-2] = v 
     620    @property 
     621    def total_pd(self): return self.details[-3] 
     622    @total_pd.setter 
     623    def total_pd(self, v): self.details[-3] = v 
     624    @property 
     625    def num_active(self): return self.details[-4] 
     626    @num_active.setter 
     627    def num_active(self, v): self.details[-4] = v 
     628 
     629    def show(self): 
     630        print("total_pd", self.total_pd) 
     631        print("num_active", self.num_active) 
     632        print("pd_par", self.pd_par) 
     633        print("pd_length", self.pd_length) 
     634        print("pd_offset", self.pd_offset) 
     635        print("pd_stride", self.pd_stride) 
     636        print("par_offsets", self.par_offset) 
     637        print("num_coord", self.num_coord) 
     638        print("par_coord", self.par_coord) 
     639        print("pd_coord", self.pd_coord) 
     640        print("theta par", self.details[-1]) 
     641 
    581642def mono_details(model_info): 
    582     # TODO: move max_pd into ParameterTable? 
    583     max_pd = model_info['max_pd'] 
    584     pars = model_info['parameters'].kernel_pars() 
    585     npars = len(pars) 
    586     par_offset = 5*max_pd 
    587     constants_offset = par_offset + 3*npars 
    588  
    589     details = np.zeros(constants_offset + 2, 'int32') 
    590     details[0*max_pd:1*max_pd] = range(max_pd)       # pd_par: arbitrary order; use first 
    591     details[1*max_pd:2*max_pd] = [1]*max_pd          # pd_length: only one element 
    592     details[2*max_pd:3*max_pd] = range(max_pd)       # pd_offset: consecutive 1.0 weights 
    593     details[3*max_pd:4*max_pd] = [1]*max_pd          # pd_stride: vectors of length 1 
    594     details[4*max_pd:5*max_pd] = [0]*max_pd          # pd_isvol: doens't matter if no norm 
    595     details[par_offset+0*npars:par_offset+1*npars] = range(2, npars+2) # par_offset: skip scale and background 
    596     details[par_offset+1*npars:par_offset+2*npars] = [0]*npars         # no coordination 
    597     #details[p+npars] = 1 # par_coord[0] is coordinated with the first par? 
    598     details[par_offset+2*npars:par_offset+3*npars] = 0 # fast coord with 0 
    599     details[constants_offset]   = 1     # fast_coord_count: one fast index 
    600     details[constants_offset+1] = -1    # theta_par: None 
     643    details = CoordinationDetails(model_info) 
     644    # The zero defaults for monodisperse systems are mostly fine 
     645    details.par_offset[:] = np.arange(2, len(details.par_offset)+2) 
    601646    return details 
    602647 
    603648def poly_details(model_info, weights): 
    604649    weights = weights[2:] 
    605  
    606     # TODO: move max_pd into ParameterTable? 
    607650    max_pd = model_info['max_pd'] 
    608     pars = model_info['parameters'].kernel_pars() 
    609     npars = len(pars) 
    610     par_offset = 5*max_pd 
    611     constants_offset = par_offset + 3*npars 
    612651 
    613652    # Decreasing list of polydispersity lengths 
    614653    # Note: the reversing view, x[::-1], does not require a copy 
    615654    pd_length = np.array([len(w) for w in weights]) 
     655    num_active = np.sum(pd_length>1) 
     656    if num_active > max_pd: 
     657        raise ValueError("Too many polydisperse parameters") 
     658 
    616659    pd_offset = np.cumsum(np.hstack((0, pd_length))) 
    617     pd_isvol = np.array([p.type=='volume' for p in pars]) 
    618     idx = np.argsort(pd_length)[::-1][:max_pd] 
    619     pd_stride = np.cumprod(np.hstack((1, pd_length[idx][:-1]))) 
    620     par_offsets = np.cumsum(np.hstack((2, pd_length)))[:-1] 
    621     coord_offset = par_offset+npars 
    622     fast_coord_offset = par_offset+2*npars 
    623  
    624     theta_par = -1 
    625     if 'theta_par' in model_info: 
    626         theta_par = model_info['theta_par'] 
    627         if theta_par >= 0 and pd_length[theta_par] <= 1: 
    628             theta_par = -1 
    629  
    630     details = np.empty(constants_offset + 2, 'int32') 
    631     details[0*max_pd:1*max_pd] = idx             # pd_par 
    632     details[1*max_pd:2*max_pd] = pd_length[idx] 
    633     details[2*max_pd:3*max_pd] = pd_offset[idx] 
    634     details[3*max_pd:4*max_pd] = pd_stride 
    635     details[4*max_pd:5*max_pd] = pd_isvol[idx] 
    636     details[par_offset+0*npars:par_offset+1*npars] = par_offsets 
    637     details[par_offset+1*npars:par_offset+2*npars] = 0  # no coordination for most 
    638     for k,parameter_num in enumerate(idx): 
    639         details[coord_offset+parameter_num] = 2**k 
    640     details[fast_coord_offset] = idx[0] 
    641     details[fast_coord_offset+1:fast_coord_offset+npars] = 0  # no fast coord with 0 
    642     details[constants_offset] = 1   # fast_coord_count: one fast index 
    643     details[constants_offset+1] = theta_par 
    644     print("polydispersity details") 
    645     print_details(model_info, details) 
     660    idx = np.argsort(pd_length)[::-1][:num_active] 
     661    par_length = np.array([max(len(w),1) for w in weights]) 
     662    pd_stride = np.cumprod(np.hstack((1, par_length[idx]))) 
     663    par_offsets = np.cumsum(np.hstack((2, par_length))) 
     664 
     665    details = CoordinationDetails(model_info) 
     666    details.pd_par[:num_active] = idx 
     667    details.pd_length[:num_active] = pd_length[idx] 
     668    details.pd_offset[:num_active] = pd_offset[idx] 
     669    details.pd_stride[:num_active] = pd_stride[:-1] 
     670    details.par_offset[:] = par_offsets[:-1] 
     671    details.total_pd = pd_stride[-1] 
     672    details.num_active = num_active 
     673    # Without constraints coordinated parameters are just the pd parameters 
     674    details.par_coord[:num_active] = idx 
     675    details.pd_coord[:num_active] = 2**np.arange(num_active) 
     676    details.num_coord = num_active 
     677    #details.show() 
    646678    return details 
    647  
    648 def print_details(model_info, details): 
    649     max_pd = model_info['max_pd'] 
    650     pars = model_info['parameters'].kernel_pars() 
    651     npars = len(pars) 
    652     par_offset = 5*max_pd 
    653     constants_offset = par_offset + 3*npars 
    654  
    655     print("pd_par", details[0*max_pd:1*max_pd]) 
    656     print("pd_length", details[1*max_pd:2*max_pd]) 
    657     print("pd_offset", details[2*max_pd:3*max_pd]) 
    658     print("pd_stride", details[3*max_pd:4*max_pd]) 
    659     print("pd_isvol", details[4*max_pd:5*max_pd]) 
    660     print("par_offsets", details[par_offset+0*npars:par_offset+1*npars]) 
    661     print("par_coord", details[par_offset+1*npars:par_offset+2*npars]) 
    662     print("fast_coord_pars", details[par_offset+2*npars:par_offset+3*npars]) 
    663     print("fast_coord_count", details[constants_offset]) 
    664     print("theta par", details[constants_offset+1]) 
    665679 
    666680def constrained_poly_details(model_info, weights, constraints): 

  • sasmodels/kernel_iq.c

    rba32cdd r5ff1b03  
    2121    int32_t pd_offset[MAX_PD];  // offset of pd weights in the value & weight vector 
    2222    int32_t pd_stride[MAX_PD];  // stride to move to the next index at this level 
    23     int32_t pd_isvol[MAX_PD];   // True if parameter is a volume weighting parameter 
    2423#endif // MAX_PD > 0 
    25     int32_t par_offset[NPARS];  // offset of par values in the value & weight vector 
    26     int32_t par_coord[NPARS];   // polydispersity coordination bitvector 
    27     int32_t fast_coord_pars[NPARS]; // ids of the fast coordination parameters 
    28     int32_t fast_coord_count;   // number of parameters coordinated with pd 1 
     24    int32_t par_offset[NPARS];  // offset of par value blocks in the value & weight vector 
     25    int32_t par_coord[NPARS];   // ids of the coordination parameters 
     26    int32_t pd_coord[NPARS];    // polydispersity coordination bitvector 
     27    int32_t num_active;         // number of non-trivial pd loops 
     28    int32_t total_pd;           // total number of voxels in hypercube 
     29    int32_t num_coord;          // number of coordinated parameters 
    2930    int32_t theta_par;          // id of spherical correction variable 
    3031} ProblemDetails; 
    … …  
    5455  double *pvec = (double *)(&local_values);  // Alias named parameters with a vector 
    5556 
    56   // Monodisperse computation 
    57   if (pd_stop == 1) { 
    58     // Shouldn't need to copy!! 
    59     for (int k=0; k < NPARS; k++) { 
    60       pvec[k] = values[k+2];  // skip scale and background 
    61     } 
    62  
    63     const double volume = CALL_VOLUME(local_values); 
     57  // Fill in the initial variables 
     58  #ifdef USE_OPENMP 
     59  #pragma omp parallel for 
     60  #endif 
     61  for (int k=0; k < NPARS; k++) { 
     62    pvec[k] = values[problem->par_offset[k]]; 
     63  } 
     64 
     65  // If it is the first round initialize the result to zero, otherwise 
     66  // assume that the previous result has been passed back. 
     67  // Note: doing this even in the monodisperse case in order to handle the 
     68  // rare case where the model parameters are invalid and zero is returned. 
     69  // So slightly increased cost for slightly smaller code size. 
     70  if (pd_start == 0) { 
    6471    #ifdef USE_OPENMP 
    6572    #pragma omp parallel for 
    6673    #endif 
     74    for (int i=0; i < nq+1; i++) { 
     75      result[i] = 0.0; 
     76    } 
     77  } 
     78 
     79  // Monodisperse computation 
     80  if (problem->num_active == 0) { 
     81    #ifdef INVALID 
     82    if (INVALID(local_values)) { return; } 
     83    #endif 
     84 
     85    const double norm = CALL_VOLUME(local_values); 
     86    #ifdef USE_OPENMP 
     87    #pragma omp parallel for 
     88    #endif 
     89    result[nq] = norm; // Total volume normalization 
    6790    for (int i=0; i < nq; i++) { 
    6891      double scattering = CALL_IQ(q, i, local_values); 
    69       if (volume != 0.0) scattering /= volume; 
    70       result[i] = values[0]*scattering + values[1]; 
     92      result[i] = values[0]*scattering/norm + values[1]; 
    7193    } 
    7294    return; 
    … …  
    7496 
    7597#if MAX_PD > 0 
    76   //printf("Entering polydispersity\n"); 
    77  
     98  //printf("Entering polydispersity from %d to %d\n", pd_start, pd_stop); 
    7899  // Since we are no longer looping over the entire polydispersity hypercube 
    79   // for each q, we need to track the normalization values for each q in a 
    80   // separate work vector. 
    81   double norm;   // contains sum over weights 
    82   double vol; // contains sum over volume 
    83   double norm_vol; // contains weights over volume 
    84  
    85   // Initialize the results to zero 
    86   if (pd_start == 0) { 
    87     norm_vol = 0.0; 
    88     norm = 0.0; 
    89     vol = 0.0; 
    90  
    91     #ifdef USE_OPENMP 
    92     #pragma omp parallel for 
    93     #endif 
    94     for (int i=0; i < nq; i++) { 
    95       result[i] = 0.0; 
    96     } 
    97   } else { 
    98     //Pulling values from previous segment 
    99     norm = result[nq]; 
    100     vol = result[nq+1]; 
    101     norm_vol = result[nq+2]; 
    102   } 
    103  
    104   // Location in the polydispersity hypercube, one index per dimension. 
    105   local int pd_index[MAX_PD]; 
    106  
    107   // polydispersity loop index positions 
    108   local int offset[NPARS];  // NPARS excludes scale/background 
    109  
    110   // Trigger the reset behaviour that happens at the end the fast loop 
    111   // by setting the initial index >= weight vector length. 
    112   pd_index[0] = problem->pd_length[0]; 
    113  
     100  // for each q, we need to track the normalization values between calls. 
     101  double norm = 0.0; 
    114102 
    115103  // need product of weights at every Iq calc, so keep product of 
    116104  // weights from the outer loops so that weight = partial_weight * fast_weight 
    117105  double partial_weight = NAN; // product of weight w4*w3*w2 but not w1 
    118   double partial_volweight = NAN; 
    119   double weight = 1.0;        // set to 1 in case there are no weights 
    120   double vol_weight = 1.0;    // set to 1 in case there are no vol weights 
    121   double spherical_correction = 1.0;  // correction for latitude variation 
     106  double spherical_correction = 1.0;  // cosine correction for latitude variation 
     107 
     108  // Location in the polydispersity hypercube, one index per dimension. 
     109  local int pd_index[MAX_PD]; 
     110 
     111  // Location of the coordinated parameters in their own sub-cubes. 
     112  local int offset[NPARS]; 
     113 
     114  // Trigger the reset behaviour that happens at the end the fast loop 
     115  // by setting the initial index >= weight vector length. 
     116  const int fast_length = problem->pd_length[0]; 
     117  pd_index[0] = fast_length; 
    122118 
    123119  // Loop over the weights then loop over q, accumulating values 
    124120  for (int loop_index=pd_start; loop_index < pd_stop; loop_index++) { 
    125121    // check if indices need to be updated 
    126     if (pd_index[0] >= problem->pd_length[0]) { 
    127  
    128       // RESET INDICES 
    129       pd_index[0] = loop_index%problem->pd_length[0]; 
     122    if (pd_index[0] == fast_length) { 
     123      //printf("should be here with %d active\n", problem->num_active); 
     124 
     125      // Compute position in polydispersity hypercube 
     126      for (int k=0; k < problem->num_active; k++) { 
     127        pd_index[k] = (loop_index/problem->pd_stride[k])%problem->pd_length[k]; 
     128        //printf("pd_index[%d] = %d\n",k,pd_index[k]); 
     129      } 
     130 
     131      // Compute partial weights 
    130132      partial_weight = 1.0; 
    131       partial_volweight = 1.0; 
    132       for (int k=1; k < MAX_PD; k++) { 
    133         pd_index[k] = (loop_index%problem->pd_length[k])/problem->pd_stride[k]; 
    134         const double wi = weights[problem->pd_offset[k]+pd_index[k]]; 
     133      //printf("partial weight %d: ", loop_index); 
     134      for (int k=1; k < problem->num_active; k++) { 
     135        double wi = weights[problem->pd_offset[k] + pd_index[k]]; 
     136        //printf("pd[%d]=par[%d]=%g ", k, problem->pd_par[k], wi); 
    135137        partial_weight *= wi; 
    136         if (problem->pd_isvol[k]) partial_volweight *= wi; 
    137       } 
     138      } 
     139      //printf("\n"); 
     140 
     141      // Update parameter offsets in weight vector 
    138142      //printf("slow %d: ", loop_index); 
    139       for (int k=0; k < NPARS; k++) { 
    140         int coord = problem->par_coord[k]; 
    141         int this_offset = problem->par_offset[k]; 
     143      for (int k=0; k < problem->num_coord; k++) { 
     144        int par = problem->par_coord[k]; 
     145        int coord = problem->pd_coord[k]; 
     146        int this_offset = problem->par_offset[par]; 
    142147        int block_size = 1; 
    143         for (int bit=0; bit < MAX_PD && coord != 0; bit++) { 
     148        for (int bit=0; coord != 0; bit++) { 
    144149          if (coord&1) { 
    145150              this_offset += block_size * pd_index[bit]; 
    146151              block_size *= problem->pd_length[bit]; 
    147152          } 
    148           coord /= 2; 
     153          coord >>= 1; 
    149154        } 
    150         offset[k] = this_offset; 
    151         pvec[k] = values[this_offset]; 
    152         //printf("p[%d]=v[%d]=%g ", k, offset[k], pvec[k]); 
     155        offset[par] = this_offset; 
     156        pvec[par] = values[this_offset]; 
     157        //printf("par[%d]=v[%d]=%g \n", k, offset[k], pvec[k]); 
     158        // if theta is not coordinated with fast index, precompute spherical correction 
     159        if (par == problem->theta_par && !(problem->par_coord[k]&1)) { 
     160          spherical_correction = fmax(fabs(cos(M_PI_180*pvec[problem->theta_par])), 1e-6); 
     161        } 
    153162      } 
    154163      //printf("\n"); 
    155       weight = partial_weight * weights[problem->pd_offset[0]+pd_index[0]]; 
    156       if (problem->theta_par >= 0) { 
    157         spherical_correction = fabs(cos(M_PI_180*pvec[problem->theta_par])); 
    158       } 
    159       if (problem->theta_par == problem->pd_par[0]) { 
    160         weight *= spherical_correction; 
    161       } 
    162       pd_index[0] += 1; 
    163  
    164     } else { 
    165  
    166       // INCREMENT INDICES 
    167       const double wi = weights[problem->pd_offset[0]+pd_index[0]]; 
    168       weight = partial_weight*wi; 
    169       if (problem->pd_isvol[0]) vol_weight *= wi; 
    170       //printf("fast %d: ", loop_index); 
    171       for (int k=0; k < problem->fast_coord_count; k++) { 
    172         const int pindex = problem->fast_coord_pars[k]; 
    173         pvec[pindex] = values[++offset[pindex]]; 
    174         //printf("p[%d]=v[%d]=%g ", pindex, offset[pindex], pvec[pindex]); 
    175       } 
    176       //printf("\n"); 
    177       if (problem->theta_par == problem->pd_par[0]) { 
    178         weight *= fabs(cos(M_PI_180*pvec[problem->theta_par])); 
    179       } 
    180       pd_index[0] += 1; 
    181     } 
     164    } 
     165 
     166    // Increment fast index 
     167    const double wi = weights[problem->pd_offset[0] + pd_index[0]++]; 
     168    double weight = partial_weight*wi; 
     169    //printf("fast %d: ", loop_index); 
     170    for (int k=0; k < problem->num_coord; k++) { 
     171      if (problem->pd_coord[k]&1) { 
     172        const int par = problem->par_coord[k]; 
     173        pvec[par] = values[offset[par]++]; 
     174        //printf("p[%d]=v[%d]=%g ", par, offset[par]-1, pvec[par]); 
     175        // if theta is coordinated with fast index, compute spherical correction each time 
     176        if (par == problem->theta_par) { 
     177          spherical_correction = fmax(fabs(cos(M_PI_180*pvec[problem->theta_par])), 1e-6); 
     178        } 
     179      } 
     180    } 
     181    //printf("\n"); 
     182 
    182183    #ifdef INVALID 
    183184    if (INVALID(local_values)) continue; 
    … …  
    187188    // Note: weight==0 must always be excluded 
    188189    if (weight > cutoff) { 
    189       norm += weight; 
    190       vol += vol_weight * CALL_VOLUME(local_values); 
    191       norm_vol += vol_weight; 
     190      // spherical correction has some nasty effects when theta is +90 or -90 
     191      // where it becomes zero.  If the entirety of the correction 
     192      weight *= spherical_correction; 
     193      norm += weight * CALL_VOLUME(local_values); 
    192194 
    193195      #ifdef USE_OPENMP 
    … …  
    202204 
    203205  // Make normalization available for the next round 
    204   result[nq] = norm; 
    205   result[nq+1] = vol; 
    206   result[nq+2] = norm_vol; 
     206  result[nq] += norm; 
    207207 
    208208  // End of the PD loop we can normalize 
    209   if (pd_stop >= problem->pd_stride[MAX_PD-1]) { 
     209  if (pd_stop >= problem->total_pd) { 
    210210    #ifdef USE_OPENMP 
    211211    #pragma omp parallel for 
    212212    #endif 
    213213    for (int i=0; i < nq; i++) { 
    214       if (vol*norm_vol != 0.0) { 
    215         result[i] *= norm_vol/vol; 
    216       } 
    217214      result[i] = values[0]*result[i]/norm + values[1]; 
    218215    } 

  • sasmodels/kerneldll.py

    r151f3bc r5ff1b03  
    263263                else np.float64 if self.q_input.dtype == generate.F64 
    264264                else np.float128) 
    265         assert details.dtype == np.int32 
     265        assert isinstance(details, generate.CoordinationDetails) 
    266266        assert weights.dtype == real and values.dtype == real 
    267267 
    268268        max_pd = self.info['max_pd'] 
    269         start, stop = 0, details[4*max_pd-1] 
    270         print("in kerneldll") 
    271         print("details", details) 
    272         print("weights", weights) 
    273         print("values", values) 
     269        start, stop = 0, details.total_pd 
     270        #print("in kerneldll") 
     271        #print("weights", weights) 
     272        #print("values", values) 
    274273        args = [ 
    275274            self.q_input.nq, # nq 
Note: See TracChangeset for help on using the changeset viewer.