Changeset 380e8c9 in sasmodels


Ignore:
Timestamp:
Mar 24, 2016 10:56:44 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:
151f3bc
Parents:
60eab2a
Message:

progress on new polydispersity loop, but still broken

Location:
sasmodels
Files:
6 edited

Legend:

Unmodified
Added
Removed

  • sasmodels/compare.py

    re9b1663d r380e8c9  
    306306    Format the parameter list for printing. 
    307307    """ 
    308     if is2d: 
    309         exclude = lambda n: False 
    310     else: 
    311         par1d = model_info['par_type']['1d'] 
    312         exclude = lambda n: n not in par1d 
    313308    lines = [] 
    314     for p in model_info['parameters']: 
    315         if exclude(p.name): continue 
     309    for p in model_info['parameters'].type['2d' if is2d else '1d']: 
    316310        fields = dict( 
    317311            value=pars.get(p.name, p.default), 

  • sasmodels/core.py

    r60eab2a r380e8c9  
    177177    value = values.get(parameter.name, parameter.default) 
    178178    if parameter.type not in ('volume', 'orientation'): 
    179         return np.array([value]), np.array([1.0]) 
     179        return [value], [1.0] 
    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])) 
     228    vw_pairs = [(get_weights(p, pars) if active(p.name) else ([p.default], [1.0])) 
    229229                for p in kernel.info['parameters']] 
    230230    values, weights = zip(*vw_pairs) 

  • sasmodels/generate.py

    r60eab2a r380e8c9  
    614614    # Note: the reversing view, x[::-1], does not require a copy 
    615615    pd_length = np.array([len(w) for w in weights]) 
    616     print (pd_length) 
    617     print (weights) 
    618616    pd_offset = np.cumsum(np.hstack((0, pd_length))) 
    619617    pd_isvol = np.array([p.type=='volume' for p in pars]) 
    620618    idx = np.argsort(pd_length)[::-1][:max_pd] 
    621     print (idx) 
    622619    pd_stride = np.cumprod(np.hstack((1, pd_length[idx][:-1]))) 
    623620    par_offsets = np.cumsum(np.hstack((2, pd_length)))[:-1] 
     621    coord_offset = par_offset+npars 
     622    fast_coord_offset = par_offset+2*npars 
    624623 
    625624    theta_par = -1 
    … …  
    637636    details[par_offset+0*npars:par_offset+1*npars] = par_offsets 
    638637    details[par_offset+1*npars:par_offset+2*npars] = 0  # no coordination for most 
    639     details[par_offset+2*npars:par_offset+3*npars] = 0  # no fast coord with 0 
    640     coord_offset = par_offset+1*npars 
    641638    for k,parameter_num in enumerate(idx): 
    642639        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 
    643642    details[constants_offset] = 1   # fast_coord_count: one fast index 
    644643    details[constants_offset+1] = theta_par 
    645     print ("details",details) 
     644    print("polydispersity details") 
     645    print_details(model_info, details) 
    646646    return details 
     647 
     648def 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]) 
    647665 
    648666def constrained_poly_details(model_info, weights, constraints): 

  • sasmodels/kernel_header.c

    r5cf3c33 r380e8c9  
    103103#  define M_4PI_3 4.18879020478639 
    104104#endif 
    105 static double square(double x) { return x*x; } 
    106 static double cube(double x) { return x*x*x; } 
    107 static double sinc(double x) { return x==0 ? 1.0 : sin(x)/x; } 
     105static inline double square(double x) { return x*x; } 
     106static inline double cube(double x) { return x*x*x; } 
     107static inline double sinc(double x) { return x==0 ? 1.0 : sin(x)/x; } 
    108108 

  • sasmodels/kernel_iq.c

    r60eab2a r380e8c9  
    5454  double *pvec = (double *)(&local_values);  // Alias named parameters with a vector 
    5555 
    56 #if MAX_PD > 0 
    57   if (problem->pd_length[0] == 1) { 
    58 #endif // MAX_PD > 0 
     56  // Monodisperse computation 
     57  if (pd_stop == 1) { 
    5958    // Shouldn't need to copy!! 
    6059    for (int k=0; k < NPARS; k++) { 
    … …  
    7271    } 
    7372    return; 
     73  } 
     74 
    7475#if MAX_PD > 0 
    75   } 
    76  
    77   // polydispersity loop index positions 
    78   local int offset[NPARS];  // NPARS excludes scale/background 
    79  
    80   printf("Entering polydispersity\n"); 
     76 
     77printf("Entering polydispersity\n"); 
    8178  // Since we are no longer looping over the entire polydispersity hypercube 
    8279  // for each q, we need to track the normalization values for each q in a 
    … …  
    107104  // Location in the polydispersity hypercube, one index per dimension. 
    108105  local int pd_index[MAX_PD]; 
     106 
     107  // polydispersity loop index positions 
     108  local int offset[NPARS];  // NPARS excludes scale/background 
    109109 
    110110  // Trigger the reset behaviour that happens at the end the fast loop 
    … …  
    132132      for (int k=1; k < MAX_PD; k++) { 
    133133        pd_index[k] = (loop_index%problem->pd_length[k])/problem->pd_stride[k]; 
    134         const double wi = weights[problem->pd_offset[0]+pd_index[0]]; 
     134        const double wi = weights[problem->pd_offset[k]+pd_index[k]]; 
    135135        partial_weight *= wi; 
    136136        if (problem->pd_isvol[k]) partial_volweight *= wi; 
    137137      } 
     138      printf("slow %d: ", loop_index); 
    138139      for (int k=0; k < NPARS; k++) { 
    139140        int coord = problem->par_coord[k]; 
    … …  
    149150        offset[k] = this_offset; 
    150151        pvec[k] = values[this_offset]; 
    151       } 
     152        printf("p[%d]=v[%d]=%g ", k, offset[k], pvec[k]); 
     153      } 
     154      printf("\n"); 
    152155      weight = partial_weight * weights[problem->pd_offset[0]+pd_index[0]]; 
    153156      if (problem->theta_par >= 0) { 
    … …  
    157160        weight *= spherical_correction; 
    158161      } 
     162      pd_index[0] += 1; 
    159163 
    160164    } else { 
    161165 
    162166      // INCREMENT INDICES 
    163       pd_index[0] += 1; 
    164167      const double wi = weights[problem->pd_offset[0]+pd_index[0]]; 
    165168      weight = partial_weight*wi; 
    166169      if (problem->pd_isvol[0]) vol_weight *= wi; 
     170      printf("fast %d: ", loop_index); 
    167171      for (int k=0; k < problem->fast_coord_count; k++) { 
    168          printf("fast loop %d coord %d idx %d ?%d offset %d %g\n", 
    169          k,problem->fast_coord_pars[k],pd_index[0], 
    170         problem->fast_coord_pars[k], 
    171             offset[problem->fast_coord_pars[k]], 
    172             values[offset[problem->fast_coord_pars[k]]] 
    173           ); 
    174         pvec[problem->fast_coord_pars[k]] 
    175             = values[offset[problem->fast_coord_pars[k]]++]; 
    176  
    177       } 
    178       if (problem->theta_par ==problem->pd_par[0]) { 
     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]) { 
    179178        weight *= fabs(cos(M_PI_180*pvec[problem->theta_par])); 
    180179      } 
     180      pd_index[0] += 1; 
    181181    } 
    182182    #ifdef INVALID 
    … …  
    201201  } 
    202202 
    203   //Makes a normalization available for the next round 
     203  // Make normalization available for the next round 
    204204  result[nq] = norm; 
    205205  result[nq+1] = vol; 
    206206  result[nq+2] = norm_vol; 
    207207 
    208   //End of the PD loop we can normalize 
     208  // End of the PD loop we can normalize 
    209209  if (pd_stop >= problem->pd_stride[MAX_PD-1]) { 
    210210    #ifdef USE_OPENMP 

  • sasmodels/kerneldll.py

    re69b36f r380e8c9  
    270270        max_pd = self.info['max_pd'] 
    271271        start, stop = 0, details[4*max_pd-1] 
    272         print(details) 
     272        print("in kerneldll") 
     273        print("details", details) 
     274        print("weights", weights) 
     275        print("values", values) 
    273276        args = [ 
    274277            self.q_input.nq, # nq 
Note: See TracChangeset for help on using the changeset viewer.