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


Ignore:
Location:
sasmodels
Files:
11 edited

Legend:

Unmodified
Added
Removed

  • sasmodels/core.py

    r3d9001f r725ee36  
    3333except ImportError: 
    3434    pass 
    35  
    36 try: 
    37     np.meshgrid([]) 
    38     meshgrid = np.meshgrid 
    39 except 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] 
    4735 
    4836# TODO: refactor composite model support 

  • sasmodels/details.py

    r40a87fa r725ee36  
    2020    np.meshgrid([]) 
    2121    meshgrid = np.meshgrid 
    22 except ValueError: 
     22except Exception: 
    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         #   pd_prod            total length of pd loop 
    74         #   pd_sum             total length of the weight vector 
     73        #   num_eval           total length of pd loop 
     74        #   num_weights        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.zeros(4*max_pd + 4, 'i4') 
     77        self.buffer = np.empty(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 
    8899    @property 
    89100    def pd_par(self): 
    … …  
    107118 
    108119    @property 
    109     def pd_prod(self): 
     120    def num_eval(self): 
    110121        """Total size of the pd mesh""" 
    111122        return self.buffer[-4] 
    112123 
    113     @pd_prod.setter 
    114     def pd_prod(self, v): 
     124    @num_eval.setter 
     125    def num_eval(self, v): 
    115126        """Total size of the pd mesh""" 
    116127        self.buffer[-4] = v 
    117128 
    118129    @property 
    119     def pd_sum(self): 
     130    def num_weights(self): 
    120131        """Total length of all the weight vectors""" 
    121132        return self.buffer[-3] 
    122133 
    123     @pd_sum.setter 
    124     def pd_sum(self, v): 
     134    @num_weights.setter 
     135    def num_weights(self, v): 
    125136        """Total length of all the weight vectors""" 
    126137        self.buffer[-3] = v 
    … …  
    146157        self.buffer[-1] = v 
    147158 
    148     def show(self): 
     159    def show(self, values=None): 
    149160        """Print the polydispersity call details to the console""" 
    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  
    160 def 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  
    176 def poly_details(model_info, weights): 
     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 
     178def make_details(model_info, length, offset, num_weights): 
    177179    """ 
    178180    Return a :class:`CallDetails` object for a polydisperse calculation 
    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) 
     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) 
    192196    max_pd = model_info.parameters.max_pd 
    193197    if num_active > max_pd: 
    194198        raise ValueError("Too many polydisperse parameters") 
    195199 
    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) 
     200    # Decreasing list of polydpersity lengths 
    200201    # Note: the reversing view, x[::-1], does not require a copy 
    201     idx = np.argsort(pd_length)[::-1][:max_pd] 
    202     pd_stride = np.cumprod(np.hstack((1, pd_length[idx]))) 
     202    idx = np.argsort(length)[::-1][:max_pd] 
     203    pd_stride = np.cumprod(np.hstack((1, length[idx]))) 
    203204 
    204205    call_details = CallDetails(model_info) 
    205206    call_details.pd_par[:max_pd] = idx 
    206     call_details.pd_length[:max_pd] = pd_length[idx] 
    207     call_details.pd_offset[:max_pd] = pd_offset[idx] 
     207    call_details.pd_length[:max_pd] = length[idx] 
     208    call_details.pd_offset[:max_pd] = offset[idx] 
    208209    call_details.pd_stride[:max_pd] = pd_stride[:-1] 
    209     call_details.pd_prod = pd_stride[-1] 
    210     call_details.pd_sum = sum(len(w) for w in weights) 
     210    call_details.num_eval = pd_stride[-1] 
     211    call_details.num_weights = num_weights 
    211212    call_details.num_active = num_active 
     213    call_details.length = length 
     214    call_details.offset = offset 
    212215    #call_details.show() 
    213216    return call_details 
     217 
     218 
     219ZEROS = tuple([0.]*31) 
     220def 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 
     247def 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 
    214265 
    215266 
    … …  
    239290        value = pars 
    240291    return value, weight 
    241  
    242  
    243  
    244 def 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  
    274 def 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

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

  • sasmodels/kernel.py

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

  • sasmodels/kernel_iq.c

    r0f00d95 rbde38b5  
    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 pd_prod;            // total number of voxels in hypercube 
    25     int32_t pd_sum;             // total length of the weights vector 
     24    int32_t num_eval;           // total number of voxels in hypercube 
     25    int32_t num_weights;        // 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 
    3031typedef struct { 
    3132    PARAMETER_TABLE 
     33} ParameterTable; 
     34typedef union { 
     35    ParameterTable table; 
     36    double vector[4*((NUM_PARS+3)/4)]; 
    3237} ParameterBlock; 
    3338#endif // _PAR_BLOCK_ 
    … …  
    7984    ) 
    8085{ 
    81  
    8286  // Storage for the current parameter values.  These will be updated as we 
    83   // walk the polydispersity cube.  local_values will be aliased to pvec. 
     87  // walk the polydispersity cube. 
    8488  ParameterBlock local_values; 
    85   double *pvec = (double *)&local_values; 
    8689 
    8790#if defined(MAGNETIC) && NUM_MAGNETIC>0 
    88   // Location of the sld parameters in the parameter pvec. 
     91  // Location of the sld parameters in the parameter vector. 
    8992  // These parameters are updated with the effective sld due to magnetism. 
    9093  #if NUM_MAGNETIC > 3 
    … …  
    110113  #endif 
    111114  for (int i=0; i < NUM_PARS; i++) { 
    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); 
     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]); 
    118122  if (pd_start == 0) { 
    119     pd_norm = 0.0; 
    120123    #ifdef USE_OPENMP 
    121124    #pragma omp parallel for 
    122125    #endif 
    123126    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 + 2; 
    132   global const double *pd_weight = pd_value + details->pd_sum; 
     131  global const double *pd_value = values + NUM_VALUES; 
     132  global const double *pd_weight = pd_value + details->num_weights; 
    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:%d, %d\n",w0,values,(w0-values),NUM_VALUES); 
     171//printf("w0:%p, values:%p, diff:%ld, %d\n",w0,values,(w0-values), NUM_VALUES); 
    172172#endif 
    173173 
    … …  
    186186  int step = pd_start; 
    187187 
    188  
    189188#if MAX_PD>4 
    190189  const double weight5 = 1.0; 
    191190  while (i4 < n4) { 
    192     pvec[p4] = v4[i4]; 
     191    local_values.vector[p4] = v4[i4]; 
    193192    double weight4 = w4[i4] * weight5; 
    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); 
     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); 
    195194#elif MAX_PD>3 
    196195    const double weight4 = 1.0; 
    … …  
    198197#if MAX_PD>3 
    199198  while (i3 < n3) { 
    200     pvec[p3] = v3[i3]; 
     199    local_values.vector[p3] = v3[i3]; 
    201200    double weight3 = w3[i3] * weight4; 
    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); 
     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); 
    203202#elif MAX_PD>2 
    204203    const double weight3 = 1.0; 
    … …  
    206205#if MAX_PD>2 
    207206  while (i2 < n2) { 
    208     pvec[p2] = v2[i2]; 
     207    local_values.vector[p2] = v2[i2]; 
    209208    double weight2 = w2[i2] * weight3; 
    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); 
     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); 
    211210#elif MAX_PD>1 
    212211    const double weight2 = 1.0; 
    … …  
    214213#if MAX_PD>1 
    215214  while (i1 < n1) { 
    216     pvec[p1] = v1[i1]; 
     215    local_values.vector[p1] = v1[i1]; 
    217216    double weight1 = w1[i1] * weight2; 
    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); 
     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); 
    219218#elif MAX_PD>0 
    220219    const double weight1 = 1.0; 
    … …  
    222221#if MAX_PD>0 
    223222  if (slow_theta) { // Theta is not in inner loop 
    224     spherical_correction = fmax(fabs(cos(M_PI_180*pvec[theta_par])), 1.e-6); 
     223    spherical_correction = fmax(fabs(cos(M_PI_180*local_values.vector[theta_par])), 1.e-6); 
    225224  } 
    226225  while(i0 < n0) { 
    227     pvec[p0] = v0[i0]; 
     226    local_values.vector[p0] = v0[i0]; 
    228227    double weight0 = w0[i0] * weight1; 
    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); 
     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); 
    230229    if (fast_theta) { // Theta is in inner loop 
    231       spherical_correction = fmax(fabs(cos(M_PI_180*pvec[p0])), 1.e-6); 
     230      spherical_correction = fmax(fabs(cos(M_PI_180*local_values.vector[p0])), 1.e-6); 
    232231    } 
    233232#else 
    … …  
    235234#endif 
    236235 
    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"); 
     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"); 
    238237//printf("sphcor: %g\n", spherical_correction); 
    239238 
    240239    #ifdef INVALID 
    241     if (!INVALID(local_values)) 
     240    if (!INVALID(local_values.table)) 
    242241    #endif 
    243242    { 
    … …  
    248247        // would be problems looking at models with theta=90. 
    249248        const double weight = weight0 * spherical_correction; 
    250         pd_norm += weight * CALL_VOLUME(local_values); 
     249        pd_norm += weight * CALL_VOLUME(local_values.table); 
    251250 
    252251        #ifdef USE_OPENMP 
    … …  
    263262          // TODO: what is the magnetic scattering at q=0 
    264263          if (qsq > 1.e-16) { 
    265             double p[4]; 
     264            double p[4];  // dd, du, ud, uu 
    266265            p[0] = (qy*cos_mspin + qx*sin_mspin)/qsq; 
    267266            p[3] = -p[0]; 
    … …  
    278277                  #define M3 NUM_PARS+13 
    279278                  #define SLD(_M_offset, _sld_offset) \ 
    280                       pvec[_sld_offset] = xs * (axis \ 
     279                      local_values.vector[_sld_offset] = xs * (axis \ 
    281280                      ? (index==1 ? -values[_M_offset+2] : values[_M_offset+2]) \ 
    282281                      : mag_sld(qx, qy, pk, values[_M_offset], values[_M_offset+1], \ 
    … …  
    296295                  } 
    297296                  #endif 
    298                   scattering += CALL_IQ(q, q_index, local_values); 
     297                  scattering += CALL_IQ(q, q_index, local_values.table); 
    299298                } 
    300299              } 
    … …  
    302301          } 
    303302#else  // !MAGNETIC 
    304           const double scattering = CALL_IQ(q, q_index, local_values); 
     303          const double scattering = CALL_IQ(q, q_index, local_values.table); 
    305304#endif // !MAGNETIC 
    306305//printf("q_index:%d %g %g %g %g\n",q_index, scattering, weight, spherical_correction, weight0); 

  • sasmodels/kernel_iq.cl

    r4f1f876 rbde38b5  
    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 pd_prod;            // total number of voxels in hypercube 
    25     int32_t pd_sum;             // total length of the weights vector 
     24    int32_t num_eval;           // total number of voxels in hypercube 
     25    int32_t num_weights;        // 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.  local_values will be aliased to pvec. 
     92  // walk the polydispersity cube. 
    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   } 
    10094 
    10195#if defined(MAGNETIC) && NUM_MAGNETIC>0 
    … …  
    117111#endif // MAGNETIC 
    118112 
    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   } 
     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]); 
    126125//if (q_index==0) printf("start %d %g %g\n", pd_start, pd_norm, this_result); 
    127126 
    128127#if MAX_PD>0 
    129   global const double *pd_value = values + NUM_VALUES + 2; 
    130   global const double *pd_weight = pd_value + details->pd_sum; 
     128  global const double *pd_value = values + NUM_VALUES; 
     129  global const double *pd_weight = pd_value + details->num_weights; 
    131130#endif 
    132131 
    … …  
    256255        // TODO: what is the magnetic scattering at q=0 
    257256        if (qsq > 1.e-16) { 
    258           double p[4];  // spin_i, spin_f 
     257          double p[4];  // dd, du, ud, uu 
    259258          p[0] = (qy*cos_mspin + qx*sin_mspin)/qsq; 
    260259          p[3] = -p[0]; 

  • sasmodels/kernelcl.py

    r40a87fa rbde38b5  
    552552        ] 
    553553        #print("Calling OpenCL") 
    554         #print("values",values) 
     554        #call_details.show(values) 
    555555        # Call kernel and retrieve results 
    556556        last_call = None 
    557557        step = 100 
    558         for start in range(0, call_details.pd_prod, step): 
    559             stop = min(start+step, call_details.pd_prod) 
     558        for start in range(0, call_details.num_eval, step): 
     559            stop = min(start + step, call_details.num_eval) 
    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) 
    565566 
    566567        # Free buffers 
    … …  
    570571        scale = values[0]/self.result[self.q_input.nq] 
    571572        background = values[1] 
    572         #print("scale",scale,background) 
     573        #print("scale",scale,values[0],self.result[self.q_input.nq]) 
    573574        return scale*self.result[:self.q_input.nq] + background 
    574575        # return self.result[:self.q_input.nq] 

  • sasmodels/kerneldll.py

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

  • sasmodels/kernelpy.py

    r40a87fa rbde38b5  
    100100    """ 
    101101    def __init__(self, kernel, model_info, q_input): 
     102        # type: (callable, ModelInfo, List[np.ndarray]) -> None 
    102103        self.dtype = np.dtype('d') 
    103104        self.info = model_info 
    … …  
    156157 
    157158    def __call__(self, call_details, values, cutoff, magnetic): 
    158         assert isinstance(call_details, details.CallDetails) 
     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) 
    159164        res = _loops(self._parameter_vector, self._form, self._volume, 
    160165                     self.q_input.nq, call_details, values, cutoff) 
    … …  
    162167 
    163168    def release(self): 
     169        # type: () -> None 
    164170        """ 
    165171        Free resources associated with the kernel. 
    … …  
    189195            return np.ones(nq, 'd')*background 
    190196 
    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:] 
     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:] 
    193199 
    194200    pd_norm = 0.0 
    … …  
    209215 
    210216    total = np.zeros(nq, 'd') 
    211     for loop_index in range(call_details.pd_prod): 
     217    for loop_index in range(call_details.num_eval): 
    212218        # update polydispersity parameter values 
    213219        if p0_index == p0_length: 

  • sasmodels/mixture.py

    r7ae2b7f rfe496dd  
    1111*ProductModel(P, S)*. 
    1212""" 
     13from __future__ import print_function 
     14 
    1315from copy import copy 
    1416import numpy as np  # type: ignore 
    … …  
    1618from .modelinfo import Parameter, ParameterTable, ModelInfo 
    1719from .kernel import KernelModel, Kernel 
     20from .details import make_details 
    1821 
    1922try: 
    2023    from typing import List 
    21     from .details import CallDetails 
    2224except ImportError: 
    2325    pass 
    … …  
    2628    # type: (List[ModelInfo]) -> ModelInfo 
    2729    """ 
    28     Create info block for product model. 
     30    Create info block for mixture model. 
    2931    """ 
    3032    flatten = [] 
    … …  
    3840    # Build new parameter list 
    3941    combined_pars = [] 
     42    demo = {} 
    4043    for k, part in enumerate(parts): 
    4144        # Parameter prefix per model, A_, B_, ... 
    … …  
    4346        # to support vector parameters 
    4447        prefix = chr(ord('A')+k) + '_' 
    45         combined_pars.append(Parameter(prefix+'scale')) 
     48        scale =  Parameter(prefix+'scale', default=1.0, 
     49                           description="model intensity for " + part.name) 
     50        combined_pars.append(scale) 
    4651        for p in part.parameters.kernel_parameters: 
    4752            p = copy(p) 
    … …  
    5156                p.length_control = prefix + p.length_control 
    5257            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) 
    5361    parameters = ParameterTable(combined_pars) 
     62    parameters.max_pd = sum(part.parameters.max_pd for part in parts) 
    5463 
    5564    model_info = ModelInfo() 
    … …  
    7079    # Remember the component info blocks so we can build the model 
    7180    model_info.composition = ('mixture', parts) 
     81    model_info.demo = demo 
     82    return model_info 
    7283 
    7384 
    … …  
    7889        self.parts = parts 
    7990 
    80     def __call__(self, q_vectors): 
     91    def make_kernel(self, q_vectors): 
    8192        # type: (List[np.ndarray]) -> MixtureKernel 
    8293        # Note: may be sending the q_vectors to the n times even though they 
    … …  
    104115        self.info =  model_info 
    105116        self.kernels = kernels 
    106  
    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] 
     117        self.dtype = self.kernels[0].dtype 
     118 
     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] 
    110122        total = 0.0 
    111123        # remember the parts for plotting later 
    112124        self.results = [] 
    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) 
     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) 
    117133 
    118134        return scale*total + background 
    … …  
    123139            k.release() 
    124140 
     141 
     142class 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

    rfd19811 rbde38b5  
    2424from . import weights 
    2525from . import modelinfo 
    26 from .details import build_details, dispersion_mesh 
     26from .details import make_kernel_args, 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 = build_details(calculator, pairs) 
     567        call_details, values, is_magnetic = make_kernel_args(calculator, pairs) 
    568568        #call_details.show() 
    569569        #print("pairs", pairs) 
Note: See TracChangeset for help on using the changeset viewer.