Changeset 303d8d6 in sasmodels


Ignore:
Timestamp:
Mar 21, 2016 2:49:21 AM (8 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:
4a72d1a, 3a45c2c
Parents:
03cac08
Message:

new calculator says hello before crashing

Location:
sasmodels
Files:
9 edited

Legend:

Unmodified
Added
Removed

  • sasmodels/bumps_model.py

    ra84a0ca r303d8d6  
    8181    from bumps.names import Parameter 
    8282 
    83     pars = {} 
     83    pars = {} # => floating point parameters 
    8484    for p in model_info['parameters']: 
    8585        value = kwargs.pop(p.name, p.default) 
    8686        pars[p.name] = Parameter.default(value, name=p.name, limits=p.limits) 
    87     for name in model_info['partype']['pd-2d']: 
     87    for name in model_info['par_type']['pd']: 
    8888        for xpart, xdefault, xlimits in [ 
    8989                ('_pd', 0., pars[name].limits), 
    … …  
    9595            pars[xname] = Parameter.default(xvalue, name=xname, limits=xlimits) 
    9696 
    97     pd_types = {} 
    98     for name in model_info['partype']['pd-2d']: 
     97    pd_types = {}  # => distribution names 
     98    for name in model_info['par_type']['pd']: 
    9999        xname = name + '_pd_type' 
    100100        xvalue = kwargs.pop(xname, 'gaussian') 

  • sasmodels/compare.py

    r72a081d r303d8d6  
    308308        exclude = lambda n: False 
    309309    else: 
    310         partype = model_info['partype'] 
    311         par1d = set(partype['fixed-1d']+partype['pd-1d']) 
     310        par1d = model_info['par_type']['1d'] 
    312311        exclude = lambda n: n not in par1d 
    313312    lines = [] 
    … …  
    870869        pars, pd_types = bumps_model.create_parameters(model_info, **opts['pars']) 
    871870        if not opts['is2d']: 
    872             active = [base + ext 
    873                       for base in model_info['partype']['pd-1d'] 
    874                       for ext in ['', '_pd', '_pd_n', '_pd_nsigma']] 
    875             active.extend(model_info['partype']['fixed-1d']) 
    876             for k in active: 
    877                 v = pars[k] 
    878                 v.range(*parameter_range(k, v.value)) 
     871            for name in model_info['par_type']['1d']: 
     872                for ext in ['', '_pd', '_pd_n', '_pd_nsigma']: 
     873                    k = name+ext 
     874                    v = pars.get(k, None) 
     875                    if v is not None: 
     876                        v.range(*parameter_range(k, v.value)) 
    879877        else: 
    880878            for k, v in pars.items(): 

  • sasmodels/core.py

    rcf52f9c r303d8d6  
    173173    return model(q_vectors) 
    174174 
    175 def get_weights(model_info, pars, name): 
     175def get_weights(parameter, values): 
    176176    """ 
    177177    Generate the distribution for parameter *name* given the parameter values 
    … …  
    181181    from the *pars* dictionary for parameter value and parameter dispersion. 
    182182    """ 
    183     relative = name in model_info['partype']['pd-rel'] 
    184     limits = model_info['limits'][name] 
    185     disperser = pars.get(name+'_pd_type', 'gaussian') 
    186     value = pars.get(name, model_info['defaults'][name]) 
    187     npts = pars.get(name+'_pd_n', 0) 
    188     width = pars.get(name+'_pd', 0.0) 
    189     nsigma = pars.get(name+'_pd_nsigma', 3.0) 
     183    value = values.get(parameter.name, parameter.default) 
     184    if parameter.type not in ('volume', 'orientation'): 
     185        return [value], [] 
     186    relative = parameter.type == 'volume' 
     187    limits = parameter.limits 
     188    disperser = values.get(parameter.name+'_pd_type', 'gaussian') 
     189    npts = values.get(parameter.name+'_pd_n', 0) 
     190    width = values.get(parameter.name+'_pd', 0.0) 
     191    nsigma = values.get(parameter.name+'_pd_nsigma', 3.0) 
    190192    value, weight = weights.get_weights( 
    191193        disperser, npts, width, nsigma, value, limits, relative) 
    … …  
    206208    return value, weight 
    207209 
    208 def call_kernel(kernel, pars, cutoff=0, mono=False): 
     210def call_kernel(kernel, values, cutoff=0, mono=False): 
    209211    """ 
    210212    Call *kernel* returned from :func:`make_kernel` with parameters *pars*. 
    … …  
    219221    *mono* is True if polydispersity should be set to none on all parameters. 
    220222    """ 
    221     fixed_pars = [pars.get(name, kernel.info['defaults'][name]) 
    222                   for name in kernel.fixed_pars] 
    223     if mono: 
    224         pd_pars = [( np.array([pars[name]]), np.array([1.0]) ) 
    225                    for name in kernel.pd_pars] 
    226     else: 
    227         pd_pars = [get_weights(kernel.info, pars, name) for name in kernel.pd_pars] 
    228     return kernel(fixed_pars, pd_pars, cutoff=cutoff) 
     223    if mono or True: 
     224        pars = np.array([values.get(p.name, p.default) 
     225                         for p in kernel.info['parameters']]) 
     226        weights = np.array([1.0]) 
     227        details = kernel.info['mono_details'] 
     228        return kernel(pars, weights, details, cutoff) 
     229    else: 
     230        pairs = [get_weights(p, values) for p in kernel.info['parameters']] 
     231        weights, pars = [v for v in zip(*pairs)] 
     232        details = generate.poly_details(kernel.info, weights, pars) 
     233        weights, pars = [np.hstack(v) for v in (weights, pars)] 
     234        return kernel(pars, weights, details, cutoff) 
    229235 
    230236def call_ER_VR(model_info, vol_pars): 
    … …  
    249255 
    250256 
    251 def call_ER(info, pars): 
    252     """ 
    253     Call the model ER function using *pars*. 
    254     *info* is either *model.info* if you have a loaded model, or *kernel.info* 
    255     if you have a model kernel prepared for evaluation. 
    256     """ 
    257     ER = info.get('ER', None) 
     257def call_ER(model_info, values): 
     258    """ 
     259    Call the model ER function using *values*. *model_info* is either 
     260    *model.info* if you have a loaded model, or *kernel.info* if you 
     261    have a model kernel prepared for evaluation. 
     262    """ 
     263    ER = model_info.get('ER', None) 
    258264    if ER is None: 
    259265        return 1.0 
    260266    else: 
    261         vol_pars = [get_weights(info, pars, name) 
    262                     for name in info['partype']['volume']] 
     267        vol_pars = [get_weights(parameter, values) 
     268                    for parameter in model_info['parameters'] 
     269                    if parameter.type == 'volume'] 
    263270        value, weight = dispersion_mesh(vol_pars) 
    264271        individual_radii = ER(*value) 
    … …  
    266273        return np.sum(weight*individual_radii) / np.sum(weight) 
    267274 
    268 def call_VR(info, pars): 
     275def call_VR(model_info, values): 
    269276    """ 
    270277    Call the model VR function using *pars*. 
    … …  
    272279    if you have a model kernel prepared for evaluation. 
    273280    """ 
    274     VR = info.get('VR', None) 
     281    VR = model_info.get('VR', None) 
    275282    if VR is None: 
    276283        return 1.0 
    277284    else: 
    278         vol_pars = [get_weights(info, pars, name) 
    279                     for name in info['partype']['volume']] 
     285        vol_pars = [get_weights(parameter, values) 
     286                    for parameter in model_info['parameters'] 
     287                    if parameter.type == 'volume'] 
    280288        value, weight = dispersion_mesh(vol_pars) 
    281289        whole, part = VR(*value) 

  • sasmodels/direct_model.py

    r02e70ff r303d8d6  
    6666            self.data_type = 'Iq' 
    6767 
    68         partype = model.info['partype'] 
    69  
    7068        if self.data_type == 'sesans': 
    7169             
    … …  
    8179            q_mono = sesans.make_all_q(data) 
    8280        elif self.data_type == 'Iqxy': 
     81            partype = model.info['par_type'] 
    8382            if not partype['orientation'] and not partype['magnetic']: 
    8483                raise ValueError("not 2D without orientation or magnetic parameters") 

  • sasmodels/generate.py

    r03cac08 r303d8d6  
    354354    return [_search(search_path, f) for f in model_info['source']] 
    355355 
    356  
    357356def convert_type(source, dtype): 
    358357    """ 
    … …  
    419418    if filename not in _template_cache or mtime > _template_cache[filename][0]: 
    420419        with open(path) as fid: 
    421             _template_cache[filename] = (mtime, fid.read()) 
     420            _template_cache[filename] = (mtime, fid.read(), path) 
    422421    return _template_cache[filename][1] 
     422 
     423def model_templates(): 
     424    # TODO: fails DRY; templates are listed in two places. 
     425    # should instead have model_info contain a list of paths 
     426    return [joinpath(TEMPLATE_ROOT, filename) 
     427            for filename in ('kernel_header.c', 'kernel_iq.c')] 
     428 
    423429 
    424430_FN_TEMPLATE = """\ 
    … …  
    450456    """ 
    451457    prefix += "." 
    452     return [prefix+p.name for p in pars] 
     458    return [prefix+p for p in pars] 
    453459 
    454460_IQXY_PATTERN = re.compile("^((inline|static) )? *(double )? *Iqxy *([(]|$)", 
    … …  
    513519    # Generate form_volume function, etc. from body only 
    514520    if model_info['form_volume'] is not None: 
    515         pnames = [p.name for p in vol_parameters] 
     521        pnames = [p for p in vol_parameters] 
    516522        source.append(_gen_fn('form_volume', pnames, model_info['form_volume'])) 
    517523    if model_info['Iq'] is not None: 
    518         pnames = ['q'] + [p.name for p in iq_parameters] 
     524        pnames = ['q'] + [p for p in iq_parameters] 
    519525        source.append(_gen_fn('Iq', pnames, model_info['Iq'])) 
    520526    if model_info['Iqxy'] is not None: 
    521         pnames = ['qx', 'qy'] + [p.name for p in iqxy_parameters] 
     527        pnames = ['qx', 'qy'] + [p for p in iqxy_parameters] 
    522528        source.append(_gen_fn('Iqxy', pnames, model_info['Iqxy'])) 
    523529 
    … …  
    552558 
    553559    # Fill in definitions for numbers of parameters 
    554     source.append("#define MAX_PD %s"%MAX_PD) 
     560    source.append("#define MAX_PD %s"%model_info['max_pd']) 
    555561    source.append("#define NPARS %d"%(len(model_info['parameters'])-2)) 
    556562 
    … …  
    584590    * *magnetic* set of parameters that are used to compute magnetic 
    585591      patterns (which includes all 1D and 2D parameters) 
    586     * *sesans* set of parameters that are used to compute sesans patterns 
    587      (which is just 1D without background) 
    588     * *pd-relative* is the set of parameters with relative distribution 
     592    * *pd_relative* is the set of parameters with relative distribution 
    589593      width (e.g., radius +/- 10%) rather than absolute distribution 
    590594      width (e.g., theta +/- 6 degrees). 
    591595    """ 
    592596    par_set = {} 
    593     par_set['1d'] = [p for p in pars if p.type not in ('orientation', 'magnetic')] 
    594     par_set['2d'] = [p for p in pars if p.type != 'magnetic'] 
    595     par_set['magnetic'] = [p for p in pars] 
    596     par_set['pd'] = [p for p in pars if p.type in ('volume', 'orientation')] 
    597     par_set['pd_relative'] = [p for p in pars if p.type == 'volume'] 
     597    par_set['1d'] = [p.name for p in pars if p.type not in ('orientation', 'magnetic')] 
     598    par_set['2d'] = [p.name for p in pars if p.type != 'magnetic'] 
     599    par_set['magnetic'] = [p.name for p in pars] 
     600    par_set['pd'] = [p.name for p in pars if p.type in ('volume', 'orientation')] 
     601    par_set['pd_relative'] = [p.name for p in pars if p.type == 'volume'] 
    598602    return par_set 
    599603 
    … …  
    621625    } 
    622626    for p in pars: 
    623         par_type[p.type if p.type else 'other'].append(p) 
     627        par_type[p.type if p.type else 'other'].append(p.name) 
    624628    return  par_type 
    625629 
    … …  
    641645    par_type.update(categorize_parameters(pars)) 
    642646    model_info['parameters'] = pars 
    643     model_info['limits'] = dict((p.name, p.limits) for p in pars) 
    644647    model_info['par_type'] = par_type 
    645     model_info['defaults'] = dict((p.name, p.default) for p in pars) 
    646648    if model_info.get('demo', None) is None: 
    647         model_info['demo'] = model_info['defaults'] 
     649        model_info['demo'] = dict((p.name, p.default) for p in pars) 
    648650    model_info['has_2d'] = par_type['orientation'] or par_type['magnetic'] 
     651 
     652def mono_details(max_pd, npars): 
     653    par_offset = 5*max_pd 
     654    const_offset = par_offset + 3*npars 
     655 
     656    mono = np.zeros(const_offset + 3, 'i') 
     657    mono[0] = 0                # pd_par: arbitrary order; use first 
     658    mono[1*max_pd] = 1         # pd_length: only one element 
     659    mono[2*max_pd] = 2         # pd_offset: skip scale and background 
     660    mono[3*max_pd] = 1         # pd_stride: vectors of length 1 
     661    mono[4*max_pd-1] = 1       # pd_stride[-1]: only one element in set 
     662    mono[4*max_pd] = 0         # pd_isvol: doens't matter if no norm 
     663    mono[par_offset:par_offset+npars] = np.arange(2, npars+2, dtype='i') 
     664    # par_offset: copied in order 
     665    mono[par_offset+npars:par_offset+2*npars] = 0 
     666    # par_coord: no coordinated parameters 
     667    mono[par_offset+npars] = 1 # par_coord[0]: except par 0 
     668    mono[par_offset+2*npars:par_offset+3*npars] = 0 
     669    # fast coord with 0 
     670    mono[const_offset] = 1     # fast_coord_count: one fast index 
     671    mono[const_offset+1] = -1  # theta_var: None 
     672    mono[const_offset+2] = 0   # fast_theta: False 
     673    return mono 
     674 
     675def poly_details(model_info, weights, pars, constraints=None): 
     676    if constraints is not None: 
     677        # Need to find the independently varying pars and sort them 
     678        # Need to build a coordination list for the dependent variables 
     679        # Need to generate a constraints function which takes values 
     680        # and weights, returning par blocks 
     681        raise ValueError("Can't handle constraints yet") 
     682 
     683    raise ValueError("polydispersity not supported") 
     684 
    649685 
    650686def create_default_functions(model_info): 
    … …  
    688724      model variants (e.g., the list of cases) or is None if there are no 
    689725      model variants 
    690     * *defaults* is the *{parameter: value}* table built from the parameter 
    691       description table. 
    692     * *limits* is the *{parameter: [min, max]}* table built from the 
    693       parameter description table. 
    694     * *partypes* categorizes the model parameters. See 
     726    * *par_type* categorizes the model parameters. See 
    695727      :func:`categorize_parameters` for details. 
    696728    * *demo* contains the *{parameter: value}* map used in compare (and maybe 
    … …  
    709741      *model_info* blocks for the composition objects.  This allows us to 
    710742      build complete product and mixture models from just the info. 
     743    * *max_pd* is the max polydispersity dimension.  This is constant and 
     744      should not be reset.  You may be able to change it when the program 
     745      starts by setting *sasmodels.generate.MAX_PD*. 
    711746 
    712747    """ 
    … …  
    719754        name = " ".join(w.capitalize() for w in kernel_id.split('_')) 
    720755    model_info = dict( 
     756        max_pd=MAX_PD, 
    721757        id=kernel_id,  # string used to load the kernel 
    722758        filename=abspath(kernel_module.__file__), 
    … …  
    736772        oldpars=getattr(kernel_module, 'oldpars', {}), 
    737773        tests=getattr(kernel_module, 'tests', []), 
     774        mono_details = mono_details(MAX_PD, len(kernel_module.parameters)) 
    738775        ) 
    739776    process_parameters(model_info) 

  • sasmodels/kernel_iq.c

    r03cac08 r303d8d6  
    4040void KERNEL_NAME( 
    4141    int nq,                 // number of q values 
     42    const int pd_start,     // where we are in the polydispersity loop 
     43    const int pd_stop,      // where we are stopping in the polydispersity loop 
    4244    global const ProblemDetails *problem, 
    4345    global const double *weights, 
    … …  
    4547    global const double *q, // nq q values, with padding to boundary 
    4648    global double *result,  // nq+3 return values, again with padding 
    47     const double cutoff,    // cutoff in the polydispersity weight product 
    48     const int pd_start,     // where we are in the polydispersity loop 
    49     const int pd_stop       // where we are stopping in the polydispersity loop 
     49    const double cutoff     // cutoff in the polydispersity weight product 
    5050    ) 
    5151{ 
    52  
     52printf("hello\n"); 
    5353  // Storage for the current parameter values.  These will be updated as we 
    5454  // walk the polydispersity cube. 
    … …  
    177177    #endif 
    178178 
     179    // Accumulate I(q) 
     180    // Note: weight==0 must always be excluded 
    179181    if (weight > cutoff) { 
    180182      norm += weight; 

  • sasmodels/kerneldll.py

    r03cac08 r303d8d6  
    142142 
    143143    source = generate.convert_type(source, dtype) 
    144     source_files = generate.model_sources(model_info) + [model_info['filename']] 
     144    source_files = (generate.model_sources(model_info) 
     145                    + [model_info['filename']] 
     146                    + generate.model_templates()) 
    145147    dll = dll_path(model_info, dtype) 
    146148    newest = max(os.path.getmtime(f) for f in source_files) 
    … …  
    172174    return DllModel(filename, model_info, dtype=dtype) 
    173175 
    174  
    175 IQ_ARGS = [c_void_p, c_void_p, c_int] 
    176 IQXY_ARGS = [c_void_p, c_void_p, c_void_p, c_int] 
    177  
    178176class DllModel(object): 
    179177    """ 
    … …  
    197195 
    198196    def _load_dll(self): 
    199         Nfixed1d = len(self.info['partype']['fixed-1d']) 
    200         Nfixed2d = len(self.info['partype']['fixed-2d']) 
    201         Npd1d = len(self.info['partype']['pd-1d']) 
    202         Npd2d = len(self.info['partype']['pd-2d']) 
    203  
    204197        #print("dll", self.dllpath) 
    205198        try: 
    … …  
    212205              else c_double if self.dtype == generate.F64 
    213206              else c_longdouble) 
    214         pd_args_1d = [c_void_p, fp] + [c_int]*Npd1d if Npd1d else [] 
    215         pd_args_2d = [c_void_p, fp] + [c_int]*Npd2d if Npd2d else [] 
     207 
     208        # int, int, int, int*, double*, double*, double*, double*, double*, double 
     209        argtypes = [c_int]*3 + [c_void_p]*5 + [fp] 
    216210        self.Iq = self.dll[generate.kernel_name(self.info, False)] 
    217         self.Iq.argtypes = IQ_ARGS + pd_args_1d + [fp]*Nfixed1d 
    218  
    219211        self.Iqxy = self.dll[generate.kernel_name(self.info, True)] 
    220         self.Iqxy.argtypes = IQXY_ARGS + pd_args_2d + [fp]*Nfixed2d 
     212        self.Iq.argtypes = argtypes 
     213        self.Iqxy.argtypes = argtypes 
    221214 
    222215    def __getstate__(self): 
    … …  
    263256        self.q_input = q_input 
    264257        self.kernel = kernel 
    265         self.res = np.empty(q_input.nq, q_input.dtype) 
     258        self.res = np.empty(q_input.nq+3, q_input.dtype) 
    266259        dim = '2d' if q_input.is_2d else '1d' 
    267         self.fixed_pars = model_info['partype']['fixed-' + dim] 
    268         self.pd_pars = model_info['partype']['pd-' + dim] 
     260        self.parameters = model_info['par_type'][dim] 
    269261 
    270262        # In dll kernel, but not in opencl kernel 
    271263        self.p_res = self.res.ctypes.data 
    272264 
    273     def __call__(self, fixed_pars, pd_pars, cutoff): 
     265    def __call__(self, details, values, weights, cutoff): 
    274266        real = (np.float32 if self.q_input.dtype == generate.F32 
    275267                else np.float64 if self.q_input.dtype == generate.F64 
    276268                else np.float128) 
    277  
    278         nq = c_int(self.q_input.nq) 
    279         if pd_pars: 
    280             cutoff = real(cutoff) 
    281             loops_N = [np.uint32(len(p[0])) for p in pd_pars] 
    282             loops = np.hstack(pd_pars) 
    283             loops = np.ascontiguousarray(loops.T, self.q_input.dtype).flatten() 
    284             p_loops = loops.ctypes.data 
    285             dispersed = [p_loops, cutoff] + loops_N 
    286         else: 
    287             dispersed = [] 
    288         fixed = [real(p) for p in fixed_pars] 
    289         args = self.q_input.q_pointers + [self.p_res, nq] + dispersed + fixed 
    290         #print(pars) 
     269        args = [ 
     270            self.q_input.nq, # nq 
     271            0, # pd_start 
     272            1, # pd_stop 
     273            details.ctypes.data, # problem 
     274            weights.ctypes.data,  # weights 
     275            values.ctypes.data,  #pars 
     276            self.q_input.q_pointers[0], #q 
     277            self.p_res,   # results 
     278            real(cutoff), # cutoff 
     279            ] 
    291280        self.kernel(*args) 
    292281 
    293         return self.res 
     282        return self.res[:-3] 
    294283 
    295284    def release(self): 

  • sasmodels/product.py

    r72a081d r303d8d6  
    9999        # a parameter map. 
    100100        par_map = {} 
    101         p_info = p_kernel.info['partype'] 
    102         s_info = s_kernel.info['partype'] 
     101        p_info = p_kernel.info['par_type'] 
     102        s_info = s_kernel.info['par_type'] 
    103103        vol_pars = set(p_info['volume']) 
    104104        if dim == '2d': 

  • sasmodels/resolution.py

    ra146eaa r303d8d6  
    504504    from scipy.integrate import romberg 
    505505 
    506     if any(k not in form.info['defaults'] for k in pars.keys()): 
    507         keys = set(form.info['defaults'].keys()) 
    508         extra = set(pars.keys()) - keys 
     506    par_set = set([p.name for p in form.info['parameters']]) 
     507    if any(k not in par_set for k in pars.keys()): 
     508        extra = set(pars.keys()) - par_set 
    509509        raise ValueError("bad parameters: [%s] not in [%s]"% 
    510510                         (", ".join(sorted(extra)), ", ".join(sorted(keys)))) 
    … …  
    558558    from scipy.integrate import romberg 
    559559 
    560     if any(k not in form.info['defaults'] for k in pars.keys()): 
    561         keys = set(form.info['defaults'].keys()) 
    562         extra = set(pars.keys()) - keys 
     560    par_set = set([p.name for p in form.info['parameters']]) 
     561    if any(k not in par_set for k in pars.keys()): 
     562        extra = set(pars.keys()) - par_set 
    563563        raise ValueError("bad parameters: [%s] not in [%s]"% 
    564564                         (", ".join(sorted(extra)), ", ".join(sorted(keys)))) 
Note: See TracChangeset for help on using the changeset viewer.