Changeset 03cac08 in sasmodels


Ignore:
Timestamp:
Mar 20, 2016 9:44:11 PM (9 years ago)
Author:
Paul Kienzle <pkienzle@…>
Branches:
master, core_shell_microgels, costrafo411, magnetic_model, release_v0.94, release_v0.95, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests
Children:
303d8d6
Parents:
d5ac45f
Message:

new generator produces code that compiles

Files:
7 edited

Legend:

Unmodified
Added
Removed

  • doc/developer/calculator.rst

    rd5ac45f r03cac08  
    77- KERNEL declares a function to be available externally 
    88- KERNEL_NAME is the name of the function being declared 
     9- MAX_PD is the maximum depth of the polydispersity loop 
    910- NPARS is the number of parameters in the kernel 
    1011- PARAMETER_TABLE is the declaration of the parameters to the kernel:: 
    … …  
    2930        double sld_solvent 
    3031 
    31 - CALL_IQ(q, nq, i, pars) is the declaration of a call to the kernel:: 
     32- CALL_IQ(q, i, pars) is the declaration of a call to the kernel:: 
    3233 
    3334    Cylinder: 
    3435 
    35         #define CALL_IQ(q, nq, i, var) \ 
    36         Iq(q[i], \ 
     36        #define CALL_IQ(q, i, var) Iq(q[i], \ 
    3737        var.length, \ 
    3838        var.radius, \ 
    … …  
    4242    Multi-shell cylinder: 
    4343 
    44         #define CALL_IQ(q, nq, i, var) \ 
    45         Iq(q[i], \ 
     44        #define CALL_IQ(q, i, var) Iq(q[i], \ 
    4645        var.num_shells, \ 
    4746        var.length, \ 
    … …  
    5049        var.sld_solvent) 
    5150 
     51    Cylinder2D: 
     52 
     53        #define CALL_IQ(q, i, var) Iqxy(q[2*i], q[2*i+1], \ 
     54        var.length, \ 
     55        var.radius, \ 
     56        var.sld, \ 
     57        var.sld_solvent, \ 
     58        var.theta, \ 
     59        var.phi) 
     60 
    5261- CALL_VOLUME(var) is similar, but for calling the form volume:: 
    5362 
    … …  
    6978        inline bool constrained(p1, p2, p3) { return expression; } 
    7079        #define INVALID(var) constrained(var.p1, var.p2, var.p3) 
    71  
    72 - IQ_FUNC could be Iq or Iqxy 
    73 - IQ_PARS could be q[i] or q[2*i],q[2*i+1] 
    7480 
    7581Our design supports a limited number of polydispersity loops, wherein 
    … …  
    200206 
    201207TODO: cutoff 
     208 
     209For accuracy we may want to introduce Kahan summation into the integration:: 
     210 
     211 
     212    double accumulated_error = 0.0; 
     213    ... 
     214    #if USE_KAHAN_SUMMATION 
     215        const double y = next - accumulated_error; 
     216        const double t = ret + y; 
     217        accumulated_error = (t - ret) - y; 
     218        ret = t; 
     219    #else 
     220        ret += next; 
     221    #endif 

  • sasmodels/generate.py

    rd5ac45f r03cac08  
    236236 
    237237TEMPLATE_ROOT = dirname(__file__) 
     238 
     239MAX_PD = 4 
    238240 
    239241F16 = np.dtype('float16') 
    … …  
    420422    return _template_cache[filename][1] 
    421423 
     424_FN_TEMPLATE = """\ 
     425double %(name)s(%(pars)s); 
     426double %(name)s(%(pars)s) { 
     427    %(body)s 
     428} 
     429 
     430 
     431""" 
     432 
    422433def _gen_fn(name, pars, body): 
    423434    """ 
    … …  
    431442         } 
    432443    """ 
    433     template = """\ 
    434 double %(name)s(%(pars)s); 
    435 double %(name)s(%(pars)s) { 
    436     %(body)s 
    437 } 
    438  
    439  
    440 """ 
    441444    par_decl = ', '.join('double ' + p for p in pars) if pars else 'void' 
    442     return template % {'name': name, 'body': body, 'pars': par_decl} 
    443  
    444 def _gen_call_pars(name, pars): 
    445     name += "." 
    446     return ",".join(name+p for p in pars) 
     445    return _FN_TEMPLATE % {'name': name, 'body': body, 'pars': par_decl} 
     446 
     447def _call_pars(prefix, pars): 
     448    """ 
     449    Return a list of *prefix.parameter* from parameter items. 
     450    """ 
     451    prefix += "." 
     452    return [prefix+p.name for p in pars] 
     453 
     454_IQXY_PATTERN = re.compile("^((inline|static) )? *(double )? *Iqxy *([(]|$)", 
     455                           flags=re.MULTILINE) 
     456def _have_Iqxy(sources): 
     457    """ 
     458    Return true if any file defines Iqxy. 
     459 
     460    Note this is not a C parser, and so can be easily confused by 
     461    non-standard syntax.  Also, it will incorrectly identify the following 
     462    as having Iqxy:: 
     463 
     464        /* 
     465        double Iqxy(qx, qy, ...) { ... fill this in later ... } 
     466        */ 
     467 
     468    If you want to comment out an Iqxy function, use // on the front of the 
     469    line instead. 
     470    """ 
     471    for code in sources: 
     472        if _IQXY_PATTERN.search(code): 
     473            return True 
     474    else: 
     475        return False 
    447476 
    448477def make_source(model_info): 
    … …  
    464493    # for computing volume even if we allow non-disperse volume parameters. 
    465494 
    466     # Load template 
    467     source = [load_template('kernel_header.c')] 
    468  
    469     # Load additional sources 
    470     source += [open(f).read() for f in model_sources(model_info)] 
    471  
    472     # Prepare defines 
    473     defines = [] 
    474  
    475     iq_parameters = [p.name 
    476                      for p in model_info['parameters'][2:]  # skip scale, background 
    477                      if p.name in model_info['par_set']['1d']] 
    478     iqxy_parameters = [p.name 
    479                        for p in model_info['parameters'][2:]  # skip scale, background 
    480                        if p.name in model_info['par_set']['2d']] 
    481     volume_parameters = model_info['par_type']['volume'] 
     495    # kernel_iq assumes scale and background are the first parameters; 
     496    # they should be first for 1d and 2d parameter lists as well. 
     497    assert model_info['parameters'][0].name == 'scale' 
     498    assert model_info['parameters'][1].name == 'background' 
     499 
     500    # Identify parameter types 
     501    iq_parameters = model_info['par_type']['1d'][2:] 
     502    iqxy_parameters = model_info['par_type']['2d'][2:] 
     503    vol_parameters = model_info['par_type']['volume'] 
     504 
     505    # Load templates and user code 
     506    kernel_header = load_template('kernel_header.c') 
     507    kernel_code = load_template('kernel_iq.c') 
     508    user_code = [open(f).read() for f in model_sources(model_info)] 
     509 
     510    # Build initial sources 
     511    source = [kernel_header] + user_code 
    482512 
    483513    # Generate form_volume function, etc. from body only 
    484514    if model_info['form_volume'] is not None: 
    485         pnames = [p.name for p in volume_parameters] 
     515        pnames = [p.name for p in vol_parameters] 
    486516        source.append(_gen_fn('form_volume', pnames, model_info['form_volume'])) 
    487517    if model_info['Iq'] is not None: 
    … …  
    492522        source.append(_gen_fn('Iqxy', pnames, model_info['Iqxy'])) 
    493523 
    494     # Fill in definitions for volume parameters 
    495     if volume_parameters: 
    496         deref_vol = ",".join("v."+p.name for p in volume_parameters) 
    497         defines.append(('CALL_VOLUME(v)', 'form_volume(%s)\n'%deref_vol)) 
     524    # Define the parameter table 
     525    source.append("#define PARAMETER_TABLE \\") 
     526    source.append("\\\n    ".join("double %s;"%p.name 
     527                                   for p in model_info['parameters'][2:])) 
     528 
     529    # Define the function calls 
     530    if vol_parameters: 
     531        refs = ",".join(_call_pars("v", vol_parameters)) 
     532        call_volume = "#define CALL_VOLUME(v) form_volume(%s)"%refs 
    498533    else: 
    499534        # Model doesn't have volume.  We could make the kernel run a little 
    500535        # faster by not using/transferring the volume normalizations, but 
    501536        # the ifdef's reduce readability more than is worthwhile. 
    502         defines.append(('CALL_VOLUME(v)', '0.0')) 
    503  
    504     # Fill in definitions for Iq parameters 
    505     defines.append(('KERNEL_NAME', model_info['name'])) 
    506     defines.append(('IQ_PARAMETERS', ', '.join(iq_parameters))) 
    507     if fixed_1d: 
    508         defines.append(('IQ_FIXED_PARAMETER_DECLARATIONS', 
    509                         ', \\\n    '.join('const double %s' % p for p in fixed_1d))) 
    510     # Fill in definitions for Iqxy parameters 
    511     defines.append(('IQXY_KERNEL_NAME', model_info['name'] + '_Iqxy')) 
    512     defines.append(('IQXY_PARAMETERS', ', '.join(iqxy_parameters))) 
    513     if fixed_2d: 
    514         defines.append(('IQXY_FIXED_PARAMETER_DECLARATIONS', 
    515                         ', \\\n    '.join('const double %s' % p for p in fixed_2d))) 
    516     if pd_2d: 
    517         defines.append(('IQXY_WEIGHT_PRODUCT', 
    518                         '*'.join(p + '_w' for p in pd_2d))) 
    519         defines.append(('IQXY_DISPERSION_LENGTH_DECLARATIONS', 
    520                         ', \\\n    '.join('const int N%s' % p for p in pd_2d))) 
    521         defines.append(('IQXY_DISPERSION_LENGTH_SUM', 
    522                         '+'.join('N' + p for p in pd_2d))) 
    523         open_loops, close_loops = build_polydispersity_loops(pd_2d) 
    524         defines.append(('IQXY_OPEN_LOOPS', 
    525                         open_loops.replace('\n', ' \\\n'))) 
    526         defines.append(('IQXY_CLOSE_LOOPS', 
    527                         close_loops.replace('\n', ' \\\n'))) 
    528     # Need to know if we have a theta parameter for Iqxy; it is not there 
    529     # for the magnetic sphere model, for example, which has a magnetic 
    530     # orientation but no shape orientation. 
    531     if 'theta' in pd_2d: 
    532         defines.append(('IQXY_HAS_THETA', '1')) 
    533  
    534     #for d in defines: print(d) 
    535     defines = '\n'.join('#define %s %s' % (k, v) for k, v in defines) 
    536     sources = '\n\n'.join(source) 
    537     return C_KERNEL_TEMPLATE % { 
    538         'DEFINES': defines, 
    539         'SOURCES': sources, 
    540         } 
     537        call_volume = "#define CALL_VOLUME(v) 0.0" 
     538    source.append(call_volume) 
     539 
     540    refs = ["q[i]"] + _call_pars("v", iq_parameters) 
     541    call_iq = "#define CALL_IQ(q,i,v) Iq(%s)" % (",".join(refs)) 
     542    if _have_Iqxy(user_code): 
     543        # Call 2D model 
     544        refs = ["q[2*i]", "q[2*i+1]"] + _call_pars("v", iqxy_parameters) 
     545        call_iqxy = "#define CALL_IQ(q,i,v) Iqxy(%s)" % (",".join(refs)) 
     546    else: 
     547        # Call 1D model with sqrt(qx^2 + qy^2) 
     548        warnings.warn("Creating Iqxy = Iq(sqrt(qx^2 + qy^2))") 
     549        # still defined:: refs = ["q[i]"] + _call_pars("v", iq_parameters) 
     550        pars_sqrt = ["sqrt(q[2*i]*q[2*i]+q[2*i+1]*q[2*i+1])"] + refs[1:] 
     551        call_iqxy = "#define CALL_IQ(q,i,v) Iq(%s)" % (",".join(pars_sqrt)) 
     552 
     553    # Fill in definitions for numbers of parameters 
     554    source.append("#define MAX_PD %s"%MAX_PD) 
     555    source.append("#define NPARS %d"%(len(model_info['parameters'])-2)) 
     556 
     557    # TODO: allow mixed python/opencl kernels? 
     558 
     559    # define the Iq kernel 
     560    source.append("#define KERNEL_NAME %s_Iq"%model_info['name']) 
     561    source.append(call_iq) 
     562    source.append(kernel_code) 
     563    source.append("#undef CALL_IQ") 
     564    source.append("#undef KERNEL_NAME") 
     565 
     566    # define the Iqxy kernel from the same source with different #defines 
     567    source.append("#define KERNEL_NAME %s_Iqxy"%model_info['name']) 
     568    source.append(call_iqxy) 
     569    source.append(kernel_code) 
     570    source.append("#undef CALL_IQ") 
     571    source.append("#undef KERNEL_NAME") 
     572 
     573    return '\n'.join(source) 
    541574 
    542575def categorize_parameters(pars): 
    … …  
    588621    } 
    589622    for p in pars: 
    590         par_type[p.type if p.type else 'other'].append(p.name) 
     623        par_type[p.type if p.type else 'other'].append(p) 
    591624    return  par_type 
    592625 
    … …  
    605638    pars = [Parameter(*p) for p in model_info['parameters']] 
    606639    # Fill in the derived attributes 
     640    par_type = collect_types(pars) 
     641    par_type.update(categorize_parameters(pars)) 
    607642    model_info['parameters'] = pars 
    608     partype = categorize_parameters(pars) 
    609643    model_info['limits'] = dict((p.name, p.limits) for p in pars) 
    610     model_info['par_type'] = collect_types(pars) 
    611     model_info['par_set'] = categorize_parameters(pars) 
     644    model_info['par_type'] = par_type 
    612645    model_info['defaults'] = dict((p.name, p.default) for p in pars) 
    613646    if model_info.get('demo', None) is None: 
    614647        model_info['demo'] = model_info['defaults'] 
    615     model_info['has_2d'] = partype['orientation'] or partype['magnetic'] 
     648    model_info['has_2d'] = par_type['orientation'] or par_type['magnetic'] 
     649 
     650def create_default_functions(model_info): 
     651    """ 
     652    Autogenerate missing functions, such as Iqxy from Iq. 
     653 
     654    This only works for Iqxy when Iq is written in python. :func:`make_source` 
     655    performs a similar role for Iq written in C. 
     656    """ 
     657    if model_info['Iq'] is not None and model_info['Iqxy'] is None: 
     658        if model_info['par_type']['1d'] != model_info['par_type']['2d']: 
     659            raise ValueError("Iqxy model is missing") 
     660        Iq = model_info['Iq'] 
     661        def Iqxy(qx, qy, **kw): 
     662            return Iq(np.sqrt(qx**2 + qy**2), **kw) 
     663        model_info['Iqxy'] = Iqxy 
    616664 
    617665def make_model_info(kernel_module): 
    … …  
    693741    functions = "ER VR form_volume Iq Iqxy shape sesans".split() 
    694742    model_info.update((k, getattr(kernel_module, k, None)) for k in functions) 
     743    create_default_functions(model_info) 
    695744    return model_info 
    696745 

  • sasmodels/kernel_header.c

    r2e44ac7 r03cac08  
    11#ifdef __OPENCL_VERSION__ 
    22# define USE_OPENCL 
     3#elif defined(_OPENMP) 
     4# define USE_OPENMP 
    35#endif 
    4  
    5 #define USE_KAHAN_SUMMATION 0 
    66 
    77// If opencl is not available, then we are compiling a C function 
    … …  
    1616         #include <float.h> 
    1717         #define kernel extern "C" __declspec( dllexport ) 
    18          inline double trunc(double x) { return x>=0?floor(x):-floor(-x); } 
    19          inline double fmin(double x, double y) { return x>y ? y : x; } 
    20          inline double fmax(double x, double y) { return x<y ? y : x; } 
    21          inline double isnan(double x) { return _isnan(x); } 
     18         static double trunc(double x) { return x>=0?floor(x):-floor(-x); } 
     19         static double fmin(double x, double y) { return x>y ? y : x; } 
     20         static double fmax(double x, double y) { return x<y ? y : x; } 
     21         static double isnan(double x) { return _isnan(x); } 
    2222         #define NAN (std::numeric_limits<double>::quiet_NaN()) // non-signalling NaN 
    2323         static double cephes_expm1(double x) { 
    … …  
    4848         #define kernel extern "C" 
    4949     #endif 
    50      inline void SINCOS(double angle, double &svar, double &cvar) { svar=sin(angle); cvar=cos(angle); } 
     50     static void SINCOS(double angle, double &svar, double &cvar) { svar=sin(angle); cvar=cos(angle); } 
    5151#  else 
    5252     #include <stdio.h> 
    … …  
    9999#  define M_4PI_3 4.18879020478639 
    100100#endif 
    101 //inline double square(double x) { return pow(x,2.0); } 
    102 //inline double square(double x) { return pown(x,2); } 
    103 inline double square(double x) { return x*x; } 
    104 inline double cube(double x) { return x*x*x; } 
    105 inline double sinc(double x) { return x==0 ? 1.0 : sin(x)/x; } 
     101static double square(double x) { return x*x; } 
     102static double cube(double x) { return x*x*x; } 
     103static double sinc(double x) { return x==0 ? 1.0 : sin(x)/x; } 
    106104 

  • sasmodels/kernel_iq.c

    rd5ac45f r03cac08  
    1212*/ 
    1313 
     14#ifndef _PAR_BLOCK_ // protected block so we can include this code twice. 
     15#define _PAR_BLOCK_ 
    1416 
    1517#define MAX_PD 4  // MAX_PD is the max number of polydisperse parameters 
    16 #define PD_2N 16  // PD_2N is the size of the coordination step table 
    1718 
    1819typedef struct { 
    … …  
    3132 
    3233typedef struct { 
    33     PARAMETER_DECL; 
     34    PARAMETER_TABLE; 
    3435} ParameterBlock; 
    35  
    36 #define FULL_KERNEL_NAME KERNEL_NAME ## _ ## IQ_FUNC 
    37 KERNEL 
    38 void FULL_KERNEL_NAME( 
     36#endif 
     37 
     38 
     39kernel 
     40void KERNEL_NAME( 
    3941    int nq,                 // number of q values 
    4042    global const ProblemDetails *problem, 
    … …  
    4244    global const double *pars, 
    4345    global const double *q, // nq q values, with padding to boundary 
    44     global double *result,  // nq return values, again with padding 
     46    global double *result,  // nq+3 return values, again with padding 
    4547    const double cutoff,    // cutoff in the polydispersity weight product 
    4648    const int pd_start,     // where we are in the polydispersity loop 
    47     const int pd_stop,      // where we are stopping in the polydispersity loop 
     49    const int pd_stop       // where we are stopping in the polydispersity loop 
    4850    ) 
    4951{ 
    … …  
    5254  // walk the polydispersity cube. 
    5355  local ParameterBlock local_pars;  // current parameter values 
    54   const double *parvec = &local_pars;  // Alias named parameters with a vector 
     56  double *pvec = (double *)(&local_pars);  // Alias named parameters with a vector 
    5557 
    5658  local int offset[NPARS-2]; 
    … …  
    6062    // Shouldn't need to copy!! 
    6163    for (int k=0; k < NPARS; k++) { 
    62       parvec[k] = pars[k+2];  // skip scale and background 
     64      pvec[k] = pars[k+2];  // skip scale and background 
    6365    } 
    6466 
    … …  
    6870    for (int i=0; i < nq; i++) { 
    6971    { 
    70       const double scattering = IQ_FUNC(IQ_PARS, IQ_PARAMETERS); 
     72      const double scattering = CALL_IQ(q, i, local_pars); 
    7173      result[i] += pars[0]*scattering + pars[1]; 
    7274    } 
    … …  
    99101    norm = result[nq]; 
    100102    vol = result[nq+1]; 
    101     norm_vol = results[nq+2]; 
     103    norm_vol = result[nq+2]; 
    102104  } 
    103105 
    104106  // Location in the polydispersity hypercube, one index per dimension. 
    105   local int pd_index[PD_MAX]; 
     107  local int pd_index[MAX_PD]; 
    106108 
    107109  // Trigger the reset behaviour that happens at the end the fast loop 
    108110  // by setting the initial index >= weight vector length. 
    109   pd_index[0] = pd_length[0]; 
     111  pd_index[0] = problem->pd_length[0]; 
     112 
    110113 
    111114  // need product of weights at every Iq calc, so keep product of 
    … …  
    120123  for (int loop_index=pd_start; loop_index < pd_stop; loop_index++) { 
    121124    // check if indices need to be updated 
    122     if (pd_index[0] >= pd_length[0]) { 
     125    if (pd_index[0] >= problem->pd_length[0]) { 
    123126 
    124127      // RESET INDICES 
    125       pd_index[0] = loop_index%pd_length[0]; 
     128      pd_index[0] = loop_index%problem->pd_length[0]; 
    126129      partial_weight = 1.0; 
    127130      partial_volweight = 1.0; 
    128131      for (int k=1; k < MAX_PD; k++) { 
    129         pd_index[k] = (loop_index%pd_length[k])/pd_stride[k]; 
    130         const double wi = weights[pd_offset[0]+pd_index[0]]; 
     132        pd_index[k] = (loop_index%problem->pd_length[k])/problem->pd_stride[k]; 
     133        const double wi = weights[problem->pd_offset[0]+pd_index[0]]; 
    131134        partial_weight *= wi; 
    132         if (pd_isvol[k]) partial_volweight *= wi; 
     135        if (problem->pd_isvol[k]) partial_volweight *= wi; 
    133136      } 
    134137      for (int k=0; k < NPARS; k++) { 
    135         int coord = par_coord[k]; 
    136         int this_offset = par_offset[k]; 
     138        int coord = problem->par_coord[k]; 
     139        int this_offset = problem->par_offset[k]; 
    137140        int block_size = 1; 
    138141        for (int bit=0; bit < MAX_PD && coord != 0; bit++) { 
    139142          if (coord&1) { 
    140143              this_offset += block_size * pd_index[bit]; 
    141               block_size *= pd_length[bit]; 
     144              block_size *= problem->pd_length[bit]; 
    142145          } 
    143146          coord /= 2; 
    144147        } 
    145148        offset[k] = this_offset; 
    146         parvec[k] = pars[this_offset]; 
    147       } 
    148       weight = partial_weight * weights[pd_offset[0]+pd_index[0]] 
    149       if (theta_var >= 0) { 
    150         spherical_correction = fabs(cos(M_PI_180*parvec[theta_var])); 
    151       } 
    152       if (!fast_theta) weight *= spherical_correction; 
     149        pvec[k] = pars[this_offset]; 
     150      } 
     151      weight = partial_weight * weights[problem->pd_offset[0]+pd_index[0]]; 
     152      if (problem->theta_var >= 0) { 
     153        spherical_correction = fabs(cos(M_PI_180*pvec[problem->theta_var])); 
     154      } 
     155      if (!problem->fast_theta) { 
     156        weight *= spherical_correction; 
     157      } 
    153158 
    154159    } else { 
    … …  
    156161      // INCREMENT INDICES 
    157162      pd_index[0] += 1; 
    158       const double wi = weights[pd_offset[0]+pd_index[0]]; 
     163      const double wi = weights[problem->pd_offset[0]+pd_index[0]]; 
    159164      weight = partial_weight*wi; 
    160       if (pd_isvol[0]) vol_weight *= wi; 
     165      if (problem->pd_isvol[0]) vol_weight *= wi; 
    161166      for (int k=0; k < problem->fast_coord_count; k++) { 
    162         parvec[ fast_coord_index[k]] 
    163             = pars[offset[fast_coord_index[k]] + pd_index[0]]; 
    164       } 
    165       if (fast_theta) weight *= fabs(cos(M_PI_180*parvec[theta_var])); 
    166  
     167        pvec[problem->fast_coord_index[k]] 
     168            = pars[offset[problem->fast_coord_index[k]]++]; 
     169      } 
     170      if (problem->fast_theta) { 
     171        weight *= fabs(cos(M_PI_180*pvec[problem->theta_var])); 
     172      } 
    167173    } 
    168174 
    … …  
    180186      #endif 
    181187      for (int i=0; i < nq; i++) { 
    182       { 
    183         const double scattering = CALL_IQ(q, nq, i, local_pars); 
     188        const double scattering = CALL_IQ(q, i, local_pars); 
    184189        result[i] += weight*scattering; 
    185190      } 
     191    } 
    186192  } 
    187193  //Makes a normalization avialable for the next round 
    … …  
    191197 
    192198  //End of the PD loop we can normalize 
    193   if (pd_stop == pd_stride[MAX_PD-1]) {} 
     199  if (pd_stop >= problem->pd_stride[MAX_PD-1]) { 
    194200    #ifdef USE_OPENMP 
    195201    #pragma omp parallel for 

  • sasmodels/kerneldll.py

    r17bbadd r03cac08  
    6161if sys.platform == 'darwin': 
    6262    #COMPILE = "gcc-mp-4.7 -shared -fPIC -std=c99 -fopenmp -O2 -Wall %s -o %s -lm -lgomp" 
    63     COMPILE = "gcc -shared -fPIC -std=c99 -O2 -Wall %(source)s -o %(output)s -lm" 
     63    COMPILE = "gcc -shared -fPIC -std=c99 -O2 -Wall %(source)s -o %(output)s -lm -fno-unknown-pragmas" 
    6464elif os.name == 'nt': 
    6565    # call vcvarsall.bat before compiling to set path, headers, libs, etc. 
    … …  
    8080        COMPILE = "gcc -shared -fPIC -std=c99 -O2 -Wall %(source)s -o %(output)s -lm" 
    8181        if "SAS_OPENMP" in os.environ: 
    82             COMPILE = COMPILE + " -fopenmp" 
     82            COMPILE += " -fopenmp" 
     83        else: 
     84            COMPILE += " -fWno-unknown-pragmas" 
    8385else: 
    8486    COMPILE = "cc -shared -fPIC -fopenmp -std=c99 -O2 -Wall %(source)s -o %(output)s -lm" 

  • sasmodels/models/cylinder.c

    r50e1e40 r03cac08  
    33double Iqxy(double qx, double qy, double sld, double solvent_sld, 
    44    double radius, double length, double theta, double phi); 
     5 
     6#define INVALID(v) (v.radius<0 || v.length<0) 
    57 
    68double form_volume(double radius, double length) 
    … …  
    1517    double length) 
    1618{ 
    17     // TODO: return NaN if radius<0 or length<0? 
    1819    // precompute qr and qh to save time in the loop 
    1920    const double qr = q*radius; 
    … …  
    4748    double phi) 
    4849{ 
    49     // TODO: return NaN if radius<0 or length<0? 
    5050    double sn, cn; // slots to hold sincos function output 
    5151 

  • sasmodels/models/gel_fit.c

    r30b4ddf r03cac08  
    1 double form_volume(void); 
    2  
    3 double Iq(double q, 
    4           double guinier_scale, 
    5           double lorentzian_scale, 
    6           double gyration_radius, 
    7           double fractal_exp, 
    8           double cor_length); 
    9  
    10 double Iqxy(double qx, double qy, 
    11           double guinier_scale, 
    12           double lorentzian_scale, 
    13           double gyration_radius, 
    14           double fractal_exp, 
    15           double cor_length); 
    16  
    17 static double _gel_fit_kernel(double q, 
     1static double Iq(double q, 
    182          double guinier_scale, 
    193          double lorentzian_scale, 
    … …  
    248    // Lorentzian Term 
    259    ////////////////////////double a(x[i]*x[i]*zeta*zeta); 
    26     double lorentzian_term = q*q*cor_length*cor_length; 
     10    double lorentzian_term = square(q*cor_length); 
    2711    lorentzian_term = 1.0 + ((fractal_exp + 1.0)/3.0)*lorentzian_term; 
    2812    lorentzian_term = pow(lorentzian_term, (fractal_exp/2.0) ); 
    … …  
    3014    // Exponential Term 
    3115    ////////////////////////double d(x[i]*x[i]*rg*rg); 
    32     double exp_term = q*q*gyration_radius*gyration_radius; 
     16    double exp_term = square(q*gyration_radius); 
    3317    exp_term = exp(-1.0*(exp_term/3.0)); 
    3418 
    … …  
    3721    return result; 
    3822} 
    39 double form_volume(void){ 
    40     // Unused, so free to return garbage. 
    41     return NAN; 
    42 } 
    43  
    44 double Iq(double q, 
    45           double guinier_scale, 
    46           double lorentzian_scale, 
    47           double gyration_radius, 
    48           double fractal_exp, 
    49           double cor_length) 
    50 { 
    51     return _gel_fit_kernel(q, 
    52                           guinier_scale, 
    53                           lorentzian_scale, 
    54                           gyration_radius, 
    55                           fractal_exp, 
    56                           cor_length); 
    57 } 
    58  
    59 // Iqxy is never called since no orientation or magnetic parameters. 
    60 double Iqxy(double qx, double qy, 
    61           double guinier_scale, 
    62           double lorentzian_scale, 
    63           double gyration_radius, 
    64           double fractal_exp, 
    65           double cor_length) 
    66 { 
    67     double q = sqrt(qx*qx + qy*qy); 
    68     return _gel_fit_kernel(q, 
    69                           guinier_scale, 
    70                           lorentzian_scale, 
    71                           gyration_radius, 
    72                           fractal_exp, 
    73                           cor_length); 
    74 } 
    75  
Note: See TracChangeset for help on using the changeset viewer.