Changes in / [6b7e2f7:524e5c4] in sasmodels


Ignore:
Location:
sasmodels
Files:
11 edited

Legend:

Unmodified
Added
Removed

  • sasmodels/core.py

    r725ee36 r3d9001f  
    3333except ImportError: 
    3434    pass 
     35 
     36try: 
     37    np.meshgrid([]) 
     38    meshgrid = np.meshgrid 
     39except Exception: 
     40    # CRUFT: np.meshgrid requires multiple vectors 
     41    def meshgrid(*args): 
     42        """Allow meshgrid with a single argument""" 
     43        if len(args) > 1: 
     44            return np.meshgrid(*args) 
     45        else: 
     46            return [np.asarray(v) for v in args] 
    3547 
    3648# TODO: refactor composite model support 

  • sasmodels/details.py

    r725ee36 r40a87fa  
    2020    np.meshgrid([]) 
    2121    meshgrid = np.meshgrid 
    22 except Exception: 
     22except ValueError: 
    2323    # CRUFT: np.meshgrid requires multiple vectors 
    2424    def meshgrid(*args): 
    … …  
    7171        #   pd_offset[max_pd]  offset of pd values in parameter array 
    7272        #   pd_stride[max_pd]  index of pd value in loop = n//stride[k] 
    73         #   num_eval           total length of pd loop 
    74         #   num_weights        total length of the weight vector 
     73        #   pd_prod            total length of pd loop 
     74        #   pd_sum             total length of the weight vector 
    7575        #   num_active         number of pd params 
    7676        #   theta_par          parameter number for theta parameter 
    77         self.buffer = np.empty(4*max_pd + 4, 'i4') 
     77        self.buffer = np.zeros(4*max_pd + 4, 'i4') 
    7878 
    7979        # generate views on different parts of the array 
    … …  
    8686        self.theta_par = parameters.theta_offset 
    8787 
    88         # offset and length are for all parameters, not just pd parameters 
    89         # They are not sent to the kernel function, though they could be. 
    90         # They are used by the composite models (sum and product) to 
    91         # figure out offsets into the combined value list. 
    92         self.offset = None  # type: np.ndarray 
    93         self.length = None  # type: np.ndarray 
    94  
    95         # keep hold of ifno show() so we can break a values vector 
    96         # into the individual components 
    97         self.info = model_info 
    98  
    9988    @property 
    10089    def pd_par(self): 
    … …  
    118107 
    119108    @property 
    120     def num_eval(self): 
     109    def pd_prod(self): 
    121110        """Total size of the pd mesh""" 
    122111        return self.buffer[-4] 
    123112 
    124     @num_eval.setter 
    125     def num_eval(self, v): 
     113    @pd_prod.setter 
     114    def pd_prod(self, v): 
    126115        """Total size of the pd mesh""" 
    127116        self.buffer[-4] = v 
    128117 
    129118    @property 
    130     def num_weights(self): 
     119    def pd_sum(self): 
    131120        """Total length of all the weight vectors""" 
    132121        return self.buffer[-3] 
    133122 
    134     @num_weights.setter 
    135     def num_weights(self, v): 
     123    @pd_sum.setter 
     124    def pd_sum(self, v): 
    136125        """Total length of all the weight vectors""" 
    137126        self.buffer[-3] = v 
    … …  
    157146        self.buffer[-1] = v 
    158147 
    159     def show(self, values=None): 
     148    def show(self): 
    160149        """Print the polydispersity call details to the console""" 
    161         print("===== %s details ===="%self.info.name) 
    162         print("num_active:%d  num_eval:%d  num_weights:%d  theta=%d" 
    163               % (self.num_active, self.num_eval, self.num_weights, self.theta_par)) 
    164         if self.pd_par.size: 
    165             print("pd_par", self.pd_par) 
    166             print("pd_length", self.pd_length) 
    167             print("pd_offset", self.pd_offset) 
    168             print("pd_stride", self.pd_stride) 
    169         if values is not None: 
    170             nvalues = self.info.parameters.nvalues 
    171             print("scale, background", values[:2]) 
    172             print("val", values[2:nvalues]) 
    173             print("pd", values[nvalues:nvalues+self.num_weights]) 
    174             print("wt", values[nvalues+self.num_weights:nvalues+2*self.num_weights]) 
    175             print("offsets", self.offset) 
    176  
    177  
    178 def make_details(model_info, length, offset, num_weights): 
     150        print("num_active", self.num_active) 
     151        print("pd_prod", self.pd_prod) 
     152        print("pd_sum", self.pd_sum) 
     153        print("theta par", self.theta_par) 
     154        print("pd_par", self.pd_par) 
     155        print("pd_length", self.pd_length) 
     156        print("pd_offset", self.pd_offset) 
     157        print("pd_stride", self.pd_stride) 
     158 
     159 
     160def mono_details(model_info): 
     161    # type: (ModelInfo) -> CallDetails 
     162    """ 
     163    Return a :class:`CallDetails` object for a monodisperse calculation 
     164    of the model defined by *model_info*. 
     165    """ 
     166    call_details = CallDetails(model_info) 
     167    call_details.pd_prod = 1 
     168    call_details.pd_sum = model_info.parameters.nvalues 
     169    call_details.pd_par[:] = np.arange(0, model_info.parameters.max_pd) 
     170    call_details.pd_length[:] = 1 
     171    call_details.pd_offset[:] = np.arange(0, model_info.parameters.max_pd) 
     172    call_details.pd_stride[:] = 1 
     173    return call_details 
     174 
     175 
     176def poly_details(model_info, weights): 
    179177    """ 
    180178    Return a :class:`CallDetails` object for a polydisperse calculation 
    181     of the model defined by *model_info*.  Polydispersity is defined by 
    182     the *length* of the polydispersity distribution for each parameter 
    183     and the *offset* of the distribution in the polydispersity array. 
    184     Monodisperse parameters should use a polydispersity length of one 
    185     with weight 1.0. *num_weights* is the total length of the polydispersity 
    186     array. 
    187     """ 
    188     # type: (ModelInfo, np.ndarray, np.ndarray, int) -> CallDetails 
    189     #pars = model_info.parameters.call_parameters[2:model_info.parameters.npars+2] 
    190     #print(", ".join(str(i)+"-"+p.id for i,p in enumerate(pars))) 
    191     #print("len:",length) 
    192     #print("off:",offset) 
    193  
    194     # Check that we arn't using too many polydispersity loops 
    195     num_active = np.sum(length > 1) 
     179    of the model defined by *model_info* for the given set of *weights*. 
     180    *weights* is a list of vectors, one for each parameter.  Monodisperse 
     181    parameters should provide a weight vector containing [1.0]. 
     182    """ 
     183    # type: (ModelInfo) -> CallDetails 
     184    #print("weights",weights) 
     185    #weights = weights[2:] # Skip scale and background 
     186 
     187    # Decreasing list of polydispersity lengths 
     188    #print([p.id for p in model_info.parameters.call_parameters]) 
     189    pd_length = np.array([len(w) 
     190                          for w in weights[2:2+model_info.parameters.npars]]) 
     191    num_active = np.sum(pd_length > 1) 
    196192    max_pd = model_info.parameters.max_pd 
    197193    if num_active > max_pd: 
    198194        raise ValueError("Too many polydisperse parameters") 
    199195 
    200     # Decreasing list of polydpersity lengths 
     196    pd_offset = np.cumsum(np.hstack((0, pd_length))) 
     197    #print(", ".join(str(i)+"-"+p.id for i,p in enumerate(model_info.parameters.call_parameters))) 
     198    #print("len:",pd_length) 
     199    #print("off:",pd_offset) 
    201200    # Note: the reversing view, x[::-1], does not require a copy 
    202     idx = np.argsort(length)[::-1][:max_pd] 
    203     pd_stride = np.cumprod(np.hstack((1, length[idx]))) 
     201    idx = np.argsort(pd_length)[::-1][:max_pd] 
     202    pd_stride = np.cumprod(np.hstack((1, pd_length[idx]))) 
    204203 
    205204    call_details = CallDetails(model_info) 
    206205    call_details.pd_par[:max_pd] = idx 
    207     call_details.pd_length[:max_pd] = length[idx] 
    208     call_details.pd_offset[:max_pd] = offset[idx] 
     206    call_details.pd_length[:max_pd] = pd_length[idx] 
     207    call_details.pd_offset[:max_pd] = pd_offset[idx] 
    209208    call_details.pd_stride[:max_pd] = pd_stride[:-1] 
    210     call_details.num_eval = pd_stride[-1] 
    211     call_details.num_weights = num_weights 
     209    call_details.pd_prod = pd_stride[-1] 
     210    call_details.pd_sum = sum(len(w) for w in weights) 
    212211    call_details.num_active = num_active 
    213     call_details.length = length 
    214     call_details.offset = offset 
    215212    #call_details.show() 
    216213    return call_details 
    217  
    218  
    219 ZEROS = tuple([0.]*31) 
    220 def make_kernel_args(kernel, pairs): 
    221     # type: (Kernel, Tuple[List[np.ndarray], List[np.ndarray]]) -> Tuple[CallDetails, np.ndarray, bool] 
    222     """ 
    223     Converts (value, weight) pairs into parameters for the kernel call. 
    224  
    225     Returns a CallDetails object indicating the polydispersity, a data object 
    226     containing the different values, and the magnetic flag indicating whether 
    227     any magnetic magnitudes are non-zero. Magnetic vectors (M0, phi, theta) are 
    228     converted to rectangular coordinates (mx, my, mz). 
    229     """ 
    230     npars = kernel.info.parameters.npars 
    231     nvalues = kernel.info.parameters.nvalues 
    232     scalars = [v[0][0] for v in pairs] 
    233     values, weights = zip(*pairs[2:npars+2]) if npars else ((),()) 
    234     length = np.array([len(w) for w in weights]) 
    235     offset = np.cumsum(np.hstack((0, length))) 
    236     call_details = make_details(kernel.info, length, offset[:-1], offset[-1]) 
    237     # Pad value array to a 32 value boundaryd 
    238     data_len = nvalues + 2*sum(len(v) for v in values) 
    239     extra = (32 - data_len%32)%32 
    240     data = np.hstack((scalars,) + values + weights + ZEROS[:extra]) 
    241     data = data.astype(kernel.dtype) 
    242     is_magnetic = convert_magnetism(kernel.info.parameters, data) 
    243     #call_details.show() 
    244     return call_details, data, is_magnetic 
    245  
    246  
    247 def convert_magnetism(parameters, values): 
    248     """ 
    249     Convert magnetism values from polar to rectangular coordinates. 
    250  
    251     Returns True if any magnetism is present. 
    252     """ 
    253     mag = values[parameters.nvalues-3*parameters.nmagnetic:parameters.nvalues] 
    254     mag = mag.reshape(-1, 3) 
    255     scale = mag[:,0] 
    256     if np.any(scale): 
    257         theta, phi = mag[:, 1]*pi/180., mag[:, 2]*pi/180. 
    258         cos_theta = cos(theta) 
    259         mag[:, 0] = scale*cos_theta*cos(phi)  # mx 
    260         mag[:, 1] = scale*sin(theta)  # my 
    261         mag[:, 2] = -scale*cos_theta*sin(phi)  # mz 
    262         return True 
    263     else: 
    264         return False 
    265214 
    266215 
    … …  
    290239        value = pars 
    291240    return value, weight 
     241 
     242 
     243 
     244def build_details(kernel, pairs): 
     245    # type: (Kernel, Tuple[List[np.ndarray], List[np.ndarray]]) -> Tuple[CallDetails, np.ndarray, bool] 
     246    """ 
     247    Converts (value, weight) pairs into parameters for the kernel call. 
     248 
     249    Returns a CallDetails object indicating the polydispersity, a data object 
     250    containing the different values, and the magnetic flag indicating whether 
     251    any magnetic magnitudes are non-zero. Magnetic vectors (M0, phi, theta) are 
     252    converted to rectangular coordinates (mx, my, mz). 
     253    """ 
     254    values, weights = zip(*pairs) 
     255    scalars = [v[0] for v in values] 
     256    if all(len(w) == 1 for w in weights): 
     257        call_details = mono_details(kernel.info) 
     258        # Pad value array to a 32 value boundary 
     259        data_len = 3*len(scalars) 
     260        extra = ((data_len+31)//32)*32 - data_len 
     261        data = np.array(scalars+scalars+[1.]*len(scalars)+[0.]*extra, 
     262                        dtype=kernel.dtype) 
     263    else: 
     264        call_details = poly_details(kernel.info, weights) 
     265        # Pad value array to a 32 value boundary 
     266        data_len = len(scalars) + 2*sum(len(v) for v in values) 
     267        extra = ((data_len+31)//32)*32 - data_len 
     268        data = np.hstack(scalars+list(values)+list(weights)+[0.]*extra) 
     269        data = data.astype(kernel.dtype) 
     270    is_magnetic = convert_magnetism(kernel.info.parameters, data) 
     271    #call_details.show() 
     272    return call_details, data, is_magnetic 
     273 
     274def convert_magnetism(parameters, values): 
     275    """ 
     276    Convert magnetism values from polar to rectangular coordinates. 
     277 
     278    Returns True if any magnetism is present. 
     279    """ 
     280    mag = values[parameters.nvalues-3*parameters.nmagnetic:parameters.nvalues] 
     281    mag = mag.reshape(-1, 3) 
     282    scale = mag[:,0] 
     283    if np.any(scale): 
     284        theta, phi = mag[:, 1]*pi/180., mag[:, 2]*pi/180. 
     285        cos_theta = cos(theta) 
     286        mag[:, 0] = scale*cos_theta*cos(phi)  # mx 
     287        mag[:, 1] = scale*sin(theta)  # my 
     288        mag[:, 2] = -scale*cos_theta*sin(phi)  # mz 
     289        return True 
     290    else: 
     291        return False 

  • sasmodels/direct_model.py

    rbde38b5 r40a87fa  
    3030from . import resolution 
    3131from . import resolution2d 
    32 from .details import make_kernel_args, dispersion_mesh 
     32from .details import build_details, dispersion_mesh 
    3333 
    3434try: 
    … …  
    7070                for p in parameters.call_parameters] 
    7171 
    72     call_details, values, is_magnetic = make_kernel_args(calculator, vw_pairs) 
     72    call_details, values, is_magnetic = build_details(calculator, vw_pairs) 
    7373    #print("values:", values) 
    7474    return calculator(call_details, values, cutoff, is_magnetic) 

  • sasmodels/kernel.py

    rbde38b5 r9eb3632  
    3939    info = None  # type: ModelInfo 
    4040    results = None # type: List[np.ndarray] 
    41     dtype = None  # type: np.dtype 
    4241 
    4342    def __call__(self, call_details, values, cutoff, magnetic): 

  • sasmodels/kernel_iq.c

    rbde38b5 r0f00d95  
    2222    int32_t pd_stride[MAX_PD];  // stride to move to the next index at this level 
    2323#endif // MAX_PD > 0 
    24     int32_t num_eval;           // total number of voxels in hypercube 
    25     int32_t num_weights;        // total length of the weights vector 
     24    int32_t pd_prod;            // total number of voxels in hypercube 
     25    int32_t pd_sum;             // total length of the weights vector 
    2626    int32_t num_active;         // number of non-trivial pd loops 
    2727    int32_t theta_par;          // id of spherical correction variable 
    2828} ProblemDetails; 
    2929 
    30 // Intel HD 4000 needs private arrays to be a multiple of 4 long 
    3130typedef struct { 
    3231    PARAMETER_TABLE 
    33 } ParameterTable; 
    34 typedef union { 
    35     ParameterTable table; 
    36     double vector[4*((NUM_PARS+3)/4)]; 
    3732} ParameterBlock; 
    3833#endif // _PAR_BLOCK_ 
    … …  
    8479    ) 
    8580{ 
     81 
    8682  // Storage for the current parameter values.  These will be updated as we 
    87   // walk the polydispersity cube. 
     83  // walk the polydispersity cube.  local_values will be aliased to pvec. 
    8884  ParameterBlock local_values; 
     85  double *pvec = (double *)&local_values; 
    8986 
    9087#if defined(MAGNETIC) && NUM_MAGNETIC>0 
    91   // Location of the sld parameters in the parameter vector. 
     88  // Location of the sld parameters in the parameter pvec. 
    9289  // These parameters are updated with the effective sld due to magnetism. 
    9390  #if NUM_MAGNETIC > 3 
    … …  
    113110  #endif 
    114111  for (int i=0; i < NUM_PARS; i++) { 
    115     local_values.vector[i] = values[2+i]; 
    116 //printf("p%d = %g\n",i, local_values.vector[i]); 
    117   } 
    118 //printf("NUM_VALUES:%d  NUM_PARS:%d  MAX_PD:%d\n", NUM_VALUES, NUM_PARS, MAX_PD); 
    119 //printf("start:%d  stop:%d\n", pd_start, pd_stop); 
    120  
    121   double pd_norm = (pd_start == 0 ? 0.0 : result[nq]); 
     112    pvec[i] = values[2+i]; 
     113//printf("p%d = %g\n",i, pvec[i]); 
     114  } 
     115 
     116  double pd_norm; 
     117//printf("start: %d %d\n",pd_start, pd_stop); 
    122118  if (pd_start == 0) { 
     119    pd_norm = 0.0; 
    123120    #ifdef USE_OPENMP 
    124121    #pragma omp parallel for 
    125122    #endif 
    126123    for (int q_index=0; q_index < nq; q_index++) result[q_index] = 0.0; 
     124//printf("initializing %d\n", nq); 
     125  } else { 
     126    pd_norm = result[nq]; 
    127127  } 
    128128//printf("start %d %g %g\n", pd_start, pd_norm, result[0]); 
    129129 
    130130#if MAX_PD>0 
    131   global const double *pd_value = values + NUM_VALUES; 
    132   global const double *pd_weight = pd_value + details->num_weights; 
     131  global const double *pd_value = values + NUM_VALUES + 2; 
     132  global const double *pd_weight = pd_value + details->pd_sum; 
    133133#endif 
    134134 
    … …  
    169169  global const double *v0 = pd_value + details->pd_offset[0]; 
    170170  global const double *w0 = pd_weight + details->pd_offset[0]; 
    171 //printf("w0:%p, values:%p, diff:%ld, %d\n",w0,values,(w0-values), NUM_VALUES); 
     171//printf("w0:%p, values:%p, diff:%d, %d\n",w0,values,(w0-values),NUM_VALUES); 
    172172#endif 
    173173 
    … …  
    186186  int step = pd_start; 
    187187 
     188 
    188189#if MAX_PD>4 
    189190  const double weight5 = 1.0; 
    190191  while (i4 < n4) { 
    191     local_values.vector[p4] = v4[i4]; 
     192    pvec[p4] = v4[i4]; 
    192193    double weight4 = w4[i4] * weight5; 
    193 //printf("step:%d level %d: p:%d i:%d n:%d value:%g weight:%g\n", step, 4, p4, i4, n4, local_values.vector[p4], weight4); 
     194//printf("step:%d level %d: p:%d i:%d n:%d value:%g weight:%g\n", step, 4, p4, i4, n4, pvec[p4], weight4); 
    194195#elif MAX_PD>3 
    195196    const double weight4 = 1.0; 
    … …  
    197198#if MAX_PD>3 
    198199  while (i3 < n3) { 
    199     local_values.vector[p3] = v3[i3]; 
     200    pvec[p3] = v3[i3]; 
    200201    double weight3 = w3[i3] * weight4; 
    201 //printf("step:%d level %d: p:%d i:%d n:%d value:%g weight:%g\n", step, 3, p3, i3, n3, local_values.vector[p3], weight3); 
     202//printf("step:%d level %d: p:%d i:%d n:%d value:%g weight:%g\n", step, 3, p3, i3, n3, pvec[p3], weight3); 
    202203#elif MAX_PD>2 
    203204    const double weight3 = 1.0; 
    … …  
    205206#if MAX_PD>2 
    206207  while (i2 < n2) { 
    207     local_values.vector[p2] = v2[i2]; 
     208    pvec[p2] = v2[i2]; 
    208209    double weight2 = w2[i2] * weight3; 
    209 //printf("step:%d level %d: p:%d i:%d n:%d value:%g weight:%g\n", step, 2, p2, i2, n2, local_values.vector[p2], weight2); 
     210//printf("step:%d level %d: p:%d i:%d n:%d value:%g weight:%g\n", step, 2, p2, i2, n2, pvec[p2], weight2); 
    210211#elif MAX_PD>1 
    211212    const double weight2 = 1.0; 
    … …  
    213214#if MAX_PD>1 
    214215  while (i1 < n1) { 
    215     local_values.vector[p1] = v1[i1]; 
     216    pvec[p1] = v1[i1]; 
    216217    double weight1 = w1[i1] * weight2; 
    217 //printf("step:%d level %d: p:%d i:%d n:%d value:%g weight:%g\n", step, 1, p1, i1, n1, local_values.vector[p1], weight1); 
     218//printf("step:%d level %d: p:%d i:%d n:%d value:%g weight:%g\n", step, 1, p1, i1, n1, pvec[p1], weight1); 
    218219#elif MAX_PD>0 
    219220    const double weight1 = 1.0; 
    … …  
    221222#if MAX_PD>0 
    222223  if (slow_theta) { // Theta is not in inner loop 
    223     spherical_correction = fmax(fabs(cos(M_PI_180*local_values.vector[theta_par])), 1.e-6); 
     224    spherical_correction = fmax(fabs(cos(M_PI_180*pvec[theta_par])), 1.e-6); 
    224225  } 
    225226  while(i0 < n0) { 
    226     local_values.vector[p0] = v0[i0]; 
     227    pvec[p0] = v0[i0]; 
    227228    double weight0 = w0[i0] * weight1; 
    228 //printf("step:%d level %d: p:%d i:%d n:%d value:%g weight:%g\n", step, 0, p0, i0, n0, local_values.vector[p0], weight0); 
     229//printf("step:%d level %d: p:%d i:%d n:%d value:%g weight:%g\n", step, 0, p0, i0, n0, pvec[p0], weight0); 
    229230    if (fast_theta) { // Theta is in inner loop 
    230       spherical_correction = fmax(fabs(cos(M_PI_180*local_values.vector[p0])), 1.e-6); 
     231      spherical_correction = fmax(fabs(cos(M_PI_180*pvec[p0])), 1.e-6); 
    231232    } 
    232233#else 
    … …  
    234235#endif 
    235236 
    236 //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"); 
     237//printf("step:%d of %d, pars:",step,pd_stop); for (int i=0; i < NUM_PARS; i++) printf("p%d=%g ",i, pvec[i]); printf("\n"); 
    237238//printf("sphcor: %g\n", spherical_correction); 
    238239 
    239240    #ifdef INVALID 
    240     if (!INVALID(local_values.table)) 
     241    if (!INVALID(local_values)) 
    241242    #endif 
    242243    { 
    … …  
    247248        // would be problems looking at models with theta=90. 
    248249        const double weight = weight0 * spherical_correction; 
    249         pd_norm += weight * CALL_VOLUME(local_values.table); 
     250        pd_norm += weight * CALL_VOLUME(local_values); 
    250251 
    251252        #ifdef USE_OPENMP 
    … …  
    262263          // TODO: what is the magnetic scattering at q=0 
    263264          if (qsq > 1.e-16) { 
    264             double p[4];  // dd, du, ud, uu 
     265            double p[4]; 
    265266            p[0] = (qy*cos_mspin + qx*sin_mspin)/qsq; 
    266267            p[3] = -p[0]; 
    … …  
    277278                  #define M3 NUM_PARS+13 
    278279                  #define SLD(_M_offset, _sld_offset) \ 
    279                       local_values.vector[_sld_offset] = xs * (axis \ 
     280                      pvec[_sld_offset] = xs * (axis \ 
    280281                      ? (index==1 ? -values[_M_offset+2] : values[_M_offset+2]) \ 
    281282                      : mag_sld(qx, qy, pk, values[_M_offset], values[_M_offset+1], \ 
    … …  
    295296                  } 
    296297                  #endif 
    297                   scattering += CALL_IQ(q, q_index, local_values.table); 
     298                  scattering += CALL_IQ(q, q_index, local_values); 
    298299                } 
    299300              } 
    … …  
    301302          } 
    302303#else  // !MAGNETIC 
    303           const double scattering = CALL_IQ(q, q_index, local_values.table); 
     304          const double scattering = CALL_IQ(q, q_index, local_values); 
    304305#endif // !MAGNETIC 
    305306//printf("q_index:%d %g %g %g %g\n",q_index, scattering, weight, spherical_correction, weight0); 

  • sasmodels/kernel_iq.cl

    rbde38b5 r4f1f876  
    2222    int32_t pd_stride[MAX_PD];  // stride to move to the next index at this level 
    2323#endif // MAX_PD > 0 
    24     int32_t num_eval;           // total number of voxels in hypercube 
    25     int32_t num_weights;        // total length of the weights vector 
     24    int32_t pd_prod;            // total number of voxels in hypercube 
     25    int32_t pd_sum;             // total length of the weights vector 
    2626    int32_t num_active;         // number of non-trivial pd loops 
    2727    int32_t theta_par;          // id of spherical correction variable 
    … …  
    9090 
    9191  // Storage for the current parameter values.  These will be updated as we 
    92   // walk the polydispersity cube. 
     92  // walk the polydispersity cube.  local_values will be aliased to pvec. 
    9393  ParameterBlock local_values; 
     94 
     95  // Fill in the initial variables 
     96  for (int i=0; i < NUM_PARS; i++) { 
     97    local_values.vector[i] = values[2+i]; 
     98//if (q_index==0) printf("p%d = %g\n",i, local_values.vector[i]); 
     99  } 
    94100 
    95101#if defined(MAGNETIC) && NUM_MAGNETIC>0 
    … …  
    111117#endif // MAGNETIC 
    112118 
    113   // Fill in the initial variables 
    114   //   values[0] is scale 
    115   //   values[1] is background 
    116   for (int i=0; i < NUM_PARS; i++) { 
    117     local_values.vector[i] = values[2+i]; 
    118 //if (q_index==0) printf("p%d = %g\n",i, local_values.vector[i]); 
    119   } 
    120 //if (q_index==0) printf("NUM_VALUES:%d  NUM_PARS:%d  MAX_PD:%d\n", NUM_VALUES, NUM_PARS, MAX_PD); 
    121 //if (q_index==0) printf("start:%d stop:%d\n", pd_start, pd_stop); 
    122  
    123   double pd_norm = (pd_start == 0 ? 0.0 : result[nq]); 
    124   double this_result = (pd_start == 0 ? 0.0 : result[q_index]); 
     119  double pd_norm, this_result; 
     120  if (pd_start == 0) { 
     121    pd_norm = this_result = 0.0; 
     122  } else { 
     123    pd_norm = result[nq]; 
     124    this_result = result[q_index]; 
     125  } 
    125126//if (q_index==0) printf("start %d %g %g\n", pd_start, pd_norm, this_result); 
    126127 
    127128#if MAX_PD>0 
    128   global const double *pd_value = values + NUM_VALUES; 
    129   global const double *pd_weight = pd_value + details->num_weights; 
     129  global const double *pd_value = values + NUM_VALUES + 2; 
     130  global const double *pd_weight = pd_value + details->pd_sum; 
    130131#endif 
    131132 
    … …  
    255256        // TODO: what is the magnetic scattering at q=0 
    256257        if (qsq > 1.e-16) { 
    257           double p[4];  // dd, du, ud, uu 
     258          double p[4];  // spin_i, spin_f 
    258259          p[0] = (qy*cos_mspin + qx*sin_mspin)/qsq; 
    259260          p[3] = -p[0]; 

  • sasmodels/kernelcl.py

    rbde38b5 r40a87fa  
    552552        ] 
    553553        #print("Calling OpenCL") 
    554         #call_details.show(values) 
     554        #print("values",values) 
    555555        # Call kernel and retrieve results 
    556556        last_call = None 
    557557        step = 100 
    558         for start in range(0, call_details.num_eval, step): 
    559             stop = min(start + step, call_details.num_eval) 
     558        for start in range(0, call_details.pd_prod, step): 
     559            stop = min(start+step, call_details.pd_prod) 
    560560            #print("queuing",start,stop) 
    561561            args[1:3] = [np.int32(start), np.int32(stop)] 
    … …  
    563563                                None, *args, wait_for=last_call)] 
    564564        cl.enqueue_copy(self.queue, self.result, self.result_b) 
    565         #print("result", self.result) 
    566565 
    567566        # Free buffers 
    … …  
    571570        scale = values[0]/self.result[self.q_input.nq] 
    572571        background = values[1] 
    573         #print("scale",scale,values[0],self.result[self.q_input.nq]) 
     572        #print("scale",scale,background) 
    574573        return scale*self.result[:self.q_input.nq] + background 
    575574        # return self.result[:self.q_input.nq] 

  • sasmodels/kerneldll.py

    rbde38b5 r40a87fa  
    380380            self.real(cutoff), # cutoff 
    381381        ] 
    382         #print("Calling DLL") 
    383         #call_details.show(values) 
     382        #print("kerneldll values", values) 
    384383        step = 100 
    385         for start in range(0, call_details.num_eval, step): 
    386             stop = min(start + step, call_details.num_eval) 
     384        for start in range(0, call_details.pd_prod, step): 
     385            stop = min(start+step, call_details.pd_prod) 
    387386            args[1:3] = [start, stop] 
     387            #print("calling DLL") 
    388388            kernel(*args) # type: ignore 
    389389 

  • sasmodels/kernelpy.py

    rbde38b5 r40a87fa  
    100100    """ 
    101101    def __init__(self, kernel, model_info, q_input): 
    102         # type: (callable, ModelInfo, List[np.ndarray]) -> None 
    103102        self.dtype = np.dtype('d') 
    104103        self.info = model_info 
    … …  
    157156 
    158157    def __call__(self, call_details, values, cutoff, magnetic): 
    159         # type: (CallDetails, np.ndarray, np.ndarray, float, bool) -> np.ndarray 
    160         if magnetic: 
    161             raise NotImplementedError("Magnetism not implemented for pure python models") 
    162         #print("Calling python kernel") 
    163         #call_details.show(values) 
     158        assert isinstance(call_details, details.CallDetails) 
    164159        res = _loops(self._parameter_vector, self._form, self._volume, 
    165160                     self.q_input.nq, call_details, values, cutoff) 
    … …  
    167162 
    168163    def release(self): 
    169         # type: () -> None 
    170164        """ 
    171165        Free resources associated with the kernel. 
    … …  
    195189            return np.ones(nq, 'd')*background 
    196190 
    197     pd_value = values[2+n_pars:2+n_pars + call_details.num_weights] 
    198     pd_weight = values[2+n_pars + call_details.num_weights:] 
     191    pd_value = values[2+n_pars:2+n_pars + call_details.pd_sum] 
     192    pd_weight = values[2+n_pars + call_details.pd_sum:] 
    199193 
    200194    pd_norm = 0.0 
    … …  
    215209 
    216210    total = np.zeros(nq, 'd') 
    217     for loop_index in range(call_details.num_eval): 
     211    for loop_index in range(call_details.pd_prod): 
    218212        # update polydispersity parameter values 
    219213        if p0_index == p0_length: 

  • sasmodels/mixture.py

    rfe496dd r7ae2b7f  
    1111*ProductModel(P, S)*. 
    1212""" 
    13 from __future__ import print_function 
    14  
    1513from copy import copy 
    1614import numpy as np  # type: ignore 
    … …  
    1816from .modelinfo import Parameter, ParameterTable, ModelInfo 
    1917from .kernel import KernelModel, Kernel 
    20 from .details import make_details 
    2118 
    2219try: 
    2320    from typing import List 
     21    from .details import CallDetails 
    2422except ImportError: 
    2523    pass 
    … …  
    2826    # type: (List[ModelInfo]) -> ModelInfo 
    2927    """ 
    30     Create info block for mixture model. 
     28    Create info block for product model. 
    3129    """ 
    3230    flatten = [] 
    … …  
    4038    # Build new parameter list 
    4139    combined_pars = [] 
    42     demo = {} 
    4340    for k, part in enumerate(parts): 
    4441        # Parameter prefix per model, A_, B_, ... 
    … …  
    4643        # to support vector parameters 
    4744        prefix = chr(ord('A')+k) + '_' 
    48         scale =  Parameter(prefix+'scale', default=1.0, 
    49                            description="model intensity for " + part.name) 
    50         combined_pars.append(scale) 
     45        combined_pars.append(Parameter(prefix+'scale')) 
    5146        for p in part.parameters.kernel_parameters: 
    5247            p = copy(p) 
    … …  
    5651                p.length_control = prefix + p.length_control 
    5752            combined_pars.append(p) 
    58         demo.update((prefix+k, v) for k, v in part.demo.items() 
    59                     if k != "background") 
    60     #print("pars",combined_pars) 
    6153    parameters = ParameterTable(combined_pars) 
    62     parameters.max_pd = sum(part.parameters.max_pd for part in parts) 
    6354 
    6455    model_info = ModelInfo() 
    … …  
    7970    # Remember the component info blocks so we can build the model 
    8071    model_info.composition = ('mixture', parts) 
    81     model_info.demo = demo 
    82     return model_info 
    8372 
    8473 
    … …  
    8978        self.parts = parts 
    9079 
    91     def make_kernel(self, q_vectors): 
     80    def __call__(self, q_vectors): 
    9281        # type: (List[np.ndarray]) -> MixtureKernel 
    9382        # Note: may be sending the q_vectors to the n times even though they 
    … …  
    115104        self.info =  model_info 
    116105        self.kernels = kernels 
    117         self.dtype = self.kernels[0].dtype 
    118106 
    119     def __call__(self, call_details, values, cutoff, magnetic): 
    120         # type: (CallDetails, np.ndarray, np.ndarry, float, bool) -> np.ndarray 
    121         scale, background = values[0:2] 
     107    def __call__(self, call_details, value, weight, cutoff): 
     108        # type: (CallDetails, np.ndarray, np.ndarry, float) -> np.ndarray 
     109        scale, background = value[0:2] 
    122110        total = 0.0 
    123111        # remember the parts for plotting later 
    124112        self.results = [] 
    125         offset = 2 # skip scale & background 
    126         parts = MixtureParts(self.info, self.kernels, call_details, values) 
    127         for kernel, kernel_details, kernel_values in parts: 
    128             #print("calling kernel", kernel.info.name) 
    129             result = kernel(kernel_details, kernel_values, cutoff, magnetic) 
    130             #print(kernel.info.name, result) 
    131             total += result 
    132             self.results.append(result) 
     113        for kernel, kernel_details in zip(self.kernels, call_details.parts): 
     114            part_result = kernel(kernel_details, value, weight, cutoff) 
     115            total += part_result 
     116            self.results.append(part_result) 
    133117 
    134118        return scale*total + background 
    … …  
    139123            k.release() 
    140124 
    141  
    142 class MixtureParts(object): 
    143     def __init__(self, model_info, kernels, call_details, values): 
    144         # type: (ModelInfo, List[Kernel], CallDetails, np.ndarray) -> None 
    145         self.model_info = model_info 
    146         self.parts = model_info.composition[1] 
    147         self.kernels = kernels 
    148         self.call_details = call_details 
    149         self.values = values 
    150         self.spin_index = model_info.parameters.npars + 2 
    151         #call_details.show(values) 
    152  
    153     def __iter__(self): 
    154         # type: () -> PartIterable 
    155         self.part_num = 0 
    156         self.par_index = 2 
    157         self.mag_index = self.spin_index + 3 
    158         return self 
    159  
    160     def next(self): 
    161         # type: () -> Tuple[List[Callable], CallDetails, np.ndarray] 
    162         if self.part_num >= len(self.parts): 
    163             raise StopIteration() 
    164         info = self.parts[self.part_num] 
    165         kernel = self.kernels[self.part_num] 
    166         call_details = self._part_details(info, self.par_index) 
    167         values = self._part_values(info, self.par_index, self.mag_index) 
    168         values = values.astype(kernel.dtype) 
    169         #call_details.show(values) 
    170  
    171         self.part_num += 1 
    172         self.par_index += info.parameters.npars + 1 
    173         self.mag_index += 3 * len(info.parameters.magnetism_index) 
    174  
    175         return kernel, call_details, values 
    176  
    177     def _part_details(self, info, par_index): 
    178         # type: (ModelInfo, int) -> CallDetails 
    179         full = self.call_details 
    180         # par_index is index into values array of the current parameter, 
    181         # which includes the initial scale and background parameters. 
    182         # We want the index into the weight length/offset for each parameter. 
    183         # Exclude the initial scale and background, so subtract two, but each 
    184         # component has its own scale factor which we need to skip when 
    185         # constructing the details for the kernel, so add one, giving a 
    186         # net subtract one. 
    187         index = slice(par_index - 1, par_index - 1 + info.parameters.npars) 
    188         length = full.length[index] 
    189         offset = full.offset[index] 
    190         # The complete weight vector is being sent to each part so that 
    191         # offsets don't need to be adjusted. 
    192         part = make_details(info, length, offset, full.num_weights) 
    193         return part 
    194  
    195     def _part_values(self, info, par_index, mag_index): 
    196         # type: (ModelInfo, int, int) -> np.ndarray 
    197         #print(info.name, par_index, self.values[par_index:par_index + info.parameters.npars + 1]) 
    198         scale = self.values[par_index] 
    199         pars = self.values[par_index + 1:par_index + info.parameters.npars + 1] 
    200         nmagnetic = len(info.parameters.magnetism_index) 
    201         if nmagnetic: 
    202             spin_state = self.values[self.spin_index:self.spin_index + 3] 
    203             mag_index = self.values[mag_index:mag_index + 3 * nmagnetic] 
    204         else: 
    205             spin_state = [] 
    206             mag_index = [] 
    207         nvalues = self.model_info.parameters.nvalues 
    208         nweights = self.call_details.num_weights 
    209         weights = self.values[nvalues:nvalues+2*nweights] 
    210         zero = self.values.dtype.type(0.) 
    211         values = [[scale, zero], pars, spin_state, mag_index, weights] 
    212         # Pad value array to a 32 value boundary 
    213         spacer = (32 - sum(len(v) for v in values)%32)%32 
    214         values.append([zero]*spacer) 
    215         values = np.hstack(values).astype(self.kernels[0].dtype) 
    216         return values 

  • sasmodels/sasview_model.py

    rbde38b5 rfd19811  
    2424from . import weights 
    2525from . import modelinfo 
    26 from .details import make_kernel_args, dispersion_mesh 
     26from .details import build_details, dispersion_mesh 
    2727 
    2828try: 
    … …  
    565565        parameters = self._model_info.parameters 
    566566        pairs = [self._get_weights(p) for p in parameters.call_parameters] 
    567         call_details, values, is_magnetic = make_kernel_args(calculator, pairs) 
     567        call_details, values, is_magnetic = build_details(calculator, pairs) 
    568568        #call_details.show() 
    569569        #print("pairs", pairs) 
Note: See TracChangeset for help on using the changeset viewer.