Changeset ec87470 in sasmodels


Ignore:
Timestamp:
Aug 3, 2016 10:26:20 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:
3a45c2c, 3d9001f
Parents:
edf06e1 (diff), d119f34 (diff)
Note: this is a merge changeset, the changes displayed below correspond to the merge itself.
Use the (diff) links above to see all the changes relative to each parent.
Message:

Merge remote-tracking branch 'origin/master' into polydisp

Location:
sasmodels
Files:
30 edited

Legend:

Unmodified
Added
Removed

  • sasmodels/compare.py

    r51ec7e8 rf67f26c  
    451451        # paying for parameter conversion each time to keep life simple, if not fast 
    452452        oldpars = revert_pars(model_info, pars) 
     453        #print("sasview pars",oldpars) 
    453454        for k, v in oldpars.items(): 
    454455            name_attr = k.split('.')  # polydispersity components 
    … …  
    667668    #    h = plt.colorbar() 
    668669    #    h.ax.set_title(cbar_title) 
     670    fig = plt.gcf() 
     671    fig.suptitle(opts['name']) 
    669672 
    670673    if Ncomp > 0 and Nbase > 0 and '-hist' in opts: 

  • sasmodels/conversion_table.py

    rf1765a2 rf67f26c  
    509509        "MultiShellModel", 
    510510        { 
     511            "radius": "core_radius", 
     512            "sld_solvent": "core_sld", 
     513            "n_pairs": "n_pairs", 
     514            "thick_shell": "s_thickness", 
    511515            "sld": "shell_sld", 
    512516            "thick_solvent": "w_thickness", 
    513             "thick_shell": "s_thickness", 
    514             "radius": "core_radius", 
    515             "sld_solvent": "core_sld" 
    516517        } 
    517518    ], 

  • sasmodels/convert.py

    r2c74c11 rd119f34  
    316316        elif name == 'core_multi_shell': 
    317317            pars['n'] = min(math.ceil(pars['n']), 4) 
     318        elif name == 'onion': 
     319            pars['n_shells'] = math.ceil(pars['n_shells']) 
    318320        elif name == 'spherical_sld': 
    319321            for k in range(1, 11): 

  • sasmodels/core.py

    ra4280bd r2547694  
    66__all__ = [ 
    77    "list_models", "load_model", "load_model_info", 
    8     "build_model", "precompile_dll", 
     8    "build_model", "precompile_dlls", 
    99    ] 
    1010 
    1111import os 
    12 from os.path import basename, dirname, join as joinpath, splitext 
     12from os.path import basename, dirname, join as joinpath 
    1313from glob import glob 
    1414 
    … …  
    5757#    build_model 
    5858 
    59 def list_models(): 
     59KINDS = ("all", "py", "c", "double", "single", "1d", "2d", 
     60         "nonmagnetic", "magnetic") 
     61def list_models(kind=None): 
    6062    # type: () -> List[str] 
    6163    """ 
    6264    Return the list of available models on the model path. 
    6365    """ 
     66    if kind and kind not in KINDS: 
     67        raise ValueError("kind not in " + ", ".join(KINDS)) 
    6468    root = dirname(__file__) 
    6569    files = sorted(glob(joinpath(root, 'models', "[a-zA-Z]*.py"))) 
    6670    available_models = [basename(f)[:-3] for f in files] 
    67     return available_models 
     71    selected = [name for name in available_models if _matches(name, kind)] 
     72 
     73    return selected 
     74 
     75def _matches(name, kind): 
     76    if kind is None or kind == "all": 
     77        return True 
     78    info = load_model_info(name) 
     79    pars = info.parameters.kernel_parameters 
     80    if kind == "py" and callable(info.Iq): 
     81        return True 
     82    elif kind == "c" and not callable(info.Iq): 
     83        return True 
     84    elif kind == "double" and not info.single: 
     85        return True 
     86    elif kind == "single" and info.single: 
     87        return True 
     88    elif kind == "2d" and any(p.type == 'orientation' for p in pars): 
     89        return True 
     90    elif kind == "1d" and any(p.type != 'orientation' for p in pars): 
     91        return True 
     92    elif kind == "magnetic" and any(p.type == 'sld' for p in pars): 
     93        return True 
     94    elif kind == "nonmagnetic" and any(p.type != 'sld' for p in pars): 
     95        return True 
     96    return False 
    6897 
    6998def load_model(model_name, dtype=None, platform='ocl'): 
    … …  
    129158            return mixture.MixtureModel(model_info, models) 
    130159        elif composition_type == 'product': 
    131             from . import product 
    132160            P, S = models 
    133161            return product.ProductModel(model_info, P, S) 
    … …  
    189217    if platform is None: 
    190218        platform = "ocl" 
    191     if platform=="ocl" and not HAVE_OPENCL: 
     219    if platform == "ocl" and not HAVE_OPENCL: 
    192220        platform = "dll" 
    193221 
    … …  
    198226 
    199227    # Convert special type names "half", "fast", and "quad" 
    200     fast = (dtype=="fast") 
     228    fast = (dtype == "fast") 
    201229    if fast: 
    202230        dtype = "single" 
    203     elif dtype=="quad": 
     231    elif dtype == "quad": 
    204232        dtype = "longdouble" 
    205     elif dtype=="half": 
     233    elif dtype == "half": 
    206234        dtype = "f16" 
    207235 
    208236    # Convert dtype string to numpy dtype. 
    209237    if dtype is None: 
    210         numpy_dtype = generate.F32 if platform=="ocl" and model_info.single else generate.F64 
     238        numpy_dtype = (generate.F32 if platform == "ocl" and model_info.single 
     239                       else generate.F64) 
    211240    else: 
    212241        numpy_dtype = np.dtype(dtype) 
    213242 
    214243    # Make sure that the type is supported by opencl, otherwise use dll 
    215     if platform=="ocl": 
     244    if platform == "ocl": 
    216245        env = kernelcl.environment() 
    217246        if not env.has_type(numpy_dtype): 
    … …  
    221250 
    222251    return numpy_dtype, fast, platform 
     252 
     253def list_models_main(): 
     254    import sys 
     255    kind = sys.argv[1] if len(sys.argv) > 1 else "all" 
     256    print("\n".join(list_models(kind))) 
     257 
     258if __name__ == "__main__": 
     259    list_models_main() 

  • sasmodels/details.py

    r9eb3632 r4f1f876  
    176176    if all(len(w)==1 for w in weights): 
    177177        call_details = mono_details(kernel.info) 
    178         data = np.array(scalars+scalars+[1]*len(scalars), dtype=kernel.dtype) 
     178        # Pad value array to a 32 value boundary 
     179        data_len = 3*len(scalars) 
     180        extra = ((data_len+31)//32)*32 - data_len 
     181        data = np.array(scalars+scalars+[1.]*len(scalars)+[0.]*extra, dtype=kernel.dtype) 
    179182    else: 
    180183        call_details = poly_details(kernel.info, weights) 
    181         data = np.hstack(scalars+list(values)+list(weights)).astype(kernel.dtype) 
     184        # Pad value array to a 32 value boundary 
     185        data_len = len(scalars) + 2*sum(len(v) for v in values) 
     186        extra = ((data_len+31)//32)*32 - data_len 
     187        data = np.hstack(scalars+list(values)+list(weights)+[0.]*extra).astype(kernel.dtype) 
    182188    is_magnetic = convert_magnetism(kernel.info.parameters, data) 
    183189    #call_details.show() 

  • sasmodels/generate.py

    ra4280bd r0f00d95  
    166166import re 
    167167import string 
    168 import warnings 
    169168 
    170169import numpy as np  # type: ignore 
    … …  
    497496    if callable(model_info.Iq): 
    498497        raise ValueError("can't compile python model") 
     498        #return None 
    499499 
    500500    # TODO: need something other than volume to indicate dispersion parameters 

  • sasmodels/kernel_iq.c

    r56547a8 r0f00d95  
    173173 
    174174 
    175   double spherical_correction=1.0; 
     175#if MAX_PD>0 
    176176  const int theta_par = details->theta_par; 
    177 #if MAX_PD>0 
    178177  const int fast_theta = (theta_par == p0); 
    179178  const int slow_theta = (theta_par >= 0 && !fast_theta); 
     179  double spherical_correction = 1.0; 
    180180#else 
    181   const int slow_theta = (theta_par >= 0); 
     181  // Note: if not polydisperse the weights cancel and we don't need the 
     182  // spherical correction. 
     183  const double spherical_correction = 1.0; 
    182184#endif 
    183185 
    … …  
    218220    const double weight1 = 1.0; 
    219221#endif 
    220     if (slow_theta) { // Theta is not in inner loop 
    221       spherical_correction = fmax(fabs(cos(M_PI_180*pvec[theta_par])), 1.e-6); 
    222     } 
    223 #if MAX_PD>0 
     222#if MAX_PD>0 
     223  if (slow_theta) { // Theta is not in inner loop 
     224    spherical_correction = fmax(fabs(cos(M_PI_180*pvec[theta_par])), 1.e-6); 
     225  } 
    224226  while(i0 < n0) { 
    225227    pvec[p0] = v0[i0]; 

  • sasmodels/kernel_iq.cl

    r56547a8 r4f1f876  
    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_ 
    … …  
    8792  // walk the polydispersity cube.  local_values will be aliased to pvec. 
    8893  ParameterBlock local_values; 
    89   double *pvec = (double *)&local_values; 
    9094 
    9195  // Fill in the initial variables 
    9296  for (int i=0; i < NUM_PARS; i++) { 
    93     pvec[i] = values[2+i]; 
    94 //if (q_index==0) printf("p%d = %g\n",i, pvec[i]); 
     97    local_values.vector[i] = values[2+i]; 
     98//if (q_index==0) printf("p%d = %g\n",i, local_values.vector[i]); 
    9599  } 
    96100 
    97101#if defined(MAGNETIC) && NUM_MAGNETIC>0 
    98   // Location of the sld parameters in the parameter pvec. 
     102  // Location of the sld parameters in the parameter vector. 
    99103  // These parameters are updated with the effective sld due to magnetism. 
    100104  #if NUM_MAGNETIC > 3 
    … …  
    166170 
    167171 
    168   double spherical_correction=1.0; 
     172#if MAX_PD>0 
    169173  const int theta_par = details->theta_par; 
    170 #if MAX_PD>0 
    171   const int fast_theta = (theta_par == p0); 
    172   const int slow_theta = (theta_par >= 0 && !fast_theta); 
     174  const bool fast_theta = (theta_par == p0); 
     175  const bool slow_theta = (theta_par >= 0 && !fast_theta); 
     176  double spherical_correction = 1.0; 
    173177#else 
    174   const int slow_theta = (theta_par >= 0); 
     178  // Note: if not polydisperse the weights cancel and we don't need the 
     179  // spherical correction. 
     180  const double spherical_correction = 1.0; 
    175181#endif 
    176182 
    … …  
    181187  const double weight5 = 1.0; 
    182188  while (i4 < n4) { 
    183     pvec[p4] = v4[i4]; 
     189    local_values.vector[p4] = v4[i4]; 
    184190    double weight4 = w4[i4] * weight5; 
    185 //if (q_index == 0) printf("step:%d level %d: p:%d i:%d n:%d value:%g weight:%g\n", step, 4, p4, i4, n4, pvec[p4], weight4); 
     191//if (q_index == 0) 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); 
    186192#elif MAX_PD>3 
    187193    const double weight4 = 1.0; 
    … …  
    189195#if MAX_PD>3 
    190196  while (i3 < n3) { 
    191     pvec[p3] = v3[i3]; 
     197    local_values.vector[p3] = v3[i3]; 
    192198    double weight3 = w3[i3] * weight4; 
    193 //if (q_index == 0) printf("step:%d level %d: p:%d i:%d n:%d value:%g weight:%g\n", step, 3, p3, i3, n3, pvec[p3], weight3); 
     199//if (q_index == 0) 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); 
    194200#elif MAX_PD>2 
    195201    const double weight3 = 1.0; 
    … …  
    197203#if MAX_PD>2 
    198204  while (i2 < n2) { 
    199     pvec[p2] = v2[i2]; 
     205    local_values.vector[p2] = v2[i2]; 
    200206    double weight2 = w2[i2] * weight3; 
    201 //if (q_index == 0) printf("step:%d level %d: p:%d i:%d n:%d value:%g weight:%g\n", step, 2, p2, i2, n2, pvec[p2], weight2); 
     207//if (q_index == 0) 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); 
    202208#elif MAX_PD>1 
    203209    const double weight2 = 1.0; 
    … …  
    205211#if MAX_PD>1 
    206212  while (i1 < n1) { 
    207     pvec[p1] = v1[i1]; 
     213    local_values.vector[p1] = v1[i1]; 
    208214    double weight1 = w1[i1] * weight2; 
    209 //if (q_index == 0) printf("step:%d level %d: p:%d i:%d n:%d value:%g weight:%g\n", step, 1, p1, i1, n1, pvec[p1], weight1); 
     215//if (q_index == 0) 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); 
    210216#elif MAX_PD>0 
    211217    const double weight1 = 1.0; 
    212218#endif 
    213     if (slow_theta) { // Theta is not in inner loop 
    214       spherical_correction = fmax(fabs(cos(M_PI_180*pvec[theta_par])), 1.e-6); 
    215     } 
    216 #if MAX_PD>0 
     219#if MAX_PD>0 
     220  if (slow_theta) { // Theta is not in inner loop 
     221    spherical_correction = fmax(fabs(cos(M_PI_180*local_values.vector[theta_par])), 1.e-6); 
     222  } 
    217223  while(i0 < n0) { 
    218     pvec[p0] = v0[i0]; 
     224    local_values.vector[p0] = v0[i0]; 
    219225    double weight0 = w0[i0] * weight1; 
    220 //if (q_index == 0) printf("step:%d level %d: p:%d i:%d n:%d value:%g weight:%g\n", step, 0, p0, i0, n0, pvec[p0], weight0); 
     226//if (q_index == 0) 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); 
    221227    if (fast_theta) { // Theta is in inner loop 
    222       spherical_correction = fmax(fabs(cos(M_PI_180*pvec[p0])), 1.e-6); 
     228      spherical_correction = fmax(fabs(cos(M_PI_180*local_values.vector[p0])), 1.e-6); 
    223229    } 
    224230#else 
    … …  
    226232#endif 
    227233 
    228 //if (q_index == 0) {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"); } 
     234//if (q_index == 0) {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"); } 
    229235//if (q_index == 0) printf("sphcor: %g\n", spherical_correction); 
    230236 
    231237    #ifdef INVALID 
    232     if (!INVALID(local_values)) 
     238    if (!INVALID(local_values.table)) 
    233239    #endif 
    234240    { 
    … …  
    239245        // would be problems looking at models with theta=90. 
    240246        const double weight = weight0 * spherical_correction; 
    241         pd_norm += weight * CALL_VOLUME(local_values); 
     247        pd_norm += weight * CALL_VOLUME(local_values.table); 
    242248 
    243249#if defined(MAGNETIC) && NUM_MAGNETIC > 0 
    … …  
    265271                #define M3 NUM_PARS+13 
    266272                #define SLD(_M_offset, _sld_offset) \ 
    267                     pvec[_sld_offset] = xs * (axis \ 
     273                    local_values.vector[_sld_offset] = xs * (axis \ 
    268274                    ? (index==1 ? -values[_M_offset+2] : values[_M_offset+2]) \ 
    269275                    : mag_sld(qx, qy, pk, values[_M_offset], values[_M_offset+1], \ 
    … …  
    283289                } 
    284290                #endif 
    285                 scattering += CALL_IQ(q, q_index, local_values); 
     291                scattering += CALL_IQ(q, q_index, local_values.table); 
    286292              } 
    287293            } 
    … …  
    289295        } 
    290296#else  // !MAGNETIC 
    291         const double scattering = CALL_IQ(q, q_index, local_values); 
     297        const double scattering = CALL_IQ(q, q_index, local_values.table); 
    292298#endif // !MAGNETIC 
    293299        this_result += weight * scattering; 

  • sasmodels/kernelcl.py

    ra4280bd r20317b3  
    6565    raise RuntimeError("OpenCL not available") 
    6666 
    67 from pyopencl import mem_flags as mf  # type: ignore 
    68 from pyopencl.characterize import get_fast_inaccurate_build_options  # type: ignore 
    69 # CRUFT: pyopencl < 2017.1  (as of June 2016 needs quotes around include path) 
    70 def _quote_path(v): 
    71     """ 
    72     Quote the path if it is not already quoted. 
    73  
    74     If v starts with '-', then assume that it is a -I option or similar 
    75     and do not quote it.  This is fragile:  -Ipath with space needs to 
    76     be quoted. 
    77     """ 
    78     return '"'+v+'"' if v and ' ' in v and not v[0] in "\"'-" else v 
    79  
    80 if hasattr(cl, '_DEFAULT_INCLUDE_OPTIONS'): 
    81     cl._DEFAULT_INCLUDE_OPTIONS = [_quote_path(v) for v in cl._DEFAULT_INCLUDE_OPTIONS] 
    82  
    8367from pyopencl import mem_flags as mf 
    8468from pyopencl.characterize import get_fast_inaccurate_build_options 
    … …  
    9377except ImportError: 
    9478    pass 
     79 
     80# CRUFT: pyopencl < 2017.1  (as of June 2016 needs quotes around include path) 
     81def quote_path(v): 
     82    """ 
     83    Quote the path if it is not already quoted. 
     84 
     85    If v starts with '-', then assume that it is a -I option or similar 
     86    and do not quote it.  This is fragile:  -Ipath with space needs to 
     87    be quoted. 
     88    """ 
     89    return '"'+v+'"' if v and ' ' in v and not v[0] in "\"'-" else v 
     90 
     91def fix_pyopencl_include(): 
     92    import pyopencl as cl 
     93    if hasattr(cl, '_DEFAULT_INCLUDE_OPTIONS'): 
     94        cl._DEFAULT_INCLUDE_OPTIONS = [quote_path(v) for v in cl._DEFAULT_INCLUDE_OPTIONS] 
     95 
     96fix_pyopencl_include() 
     97 
    9598 
    9699# The max loops number is limited by the amount of local memory available 
    … …  
    256259        Return a OpenCL context for the kernels of type dtype. 
    257260        """ 
    258         for context, queue in zip(self.context, self.queues): 
     261        for context in self.context: 
    259262            if all(has_type(d, dtype) for d in context.devices): 
    260263                return context 
    … …  
    282285        if key not in self.compiled: 
    283286            context = self.get_context(dtype) 
    284             logging.info("building %s for OpenCL %s" 
    285                          % (key, context.devices[0].name.strip())) 
     287            logging.info("building %s for OpenCL %s", key, 
     288                         context.devices[0].name.strip()) 
    286289            program = compile_model(self.get_context(dtype), 
    287290                                    str(source), dtype, fast) 
    … …  
    299302 
    300303def _get_default_context(): 
    301     # type: () -> cl.Context 
     304    # type: () -> List[cl.Context] 
    302305    """ 
    303306    Get an OpenCL context, preferring GPU over CPU, and preferring Intel 
    … …  
    318321    for platform in cl.get_platforms(): 
    319322        # AMD provides a much weaker CPU driver than Intel/Apple, so avoid it. 
    320         # If someone has bothered to install the AMD/NVIDIA drivers, prefer them over the integrated 
    321         # graphics driver that may have been supplied with the CPU chipset. 
    322         preferred_cpu = platform.vendor.startswith('Intel') or platform.vendor.startswith('Apple') 
    323         preferred_gpu = platform.vendor.startswith('Advanced') or platform.vendor.startswith('NVIDIA') 
     323        # If someone has bothered to install the AMD/NVIDIA drivers, prefer 
     324        # them over the integrated graphics driver that may have been supplied 
     325        # with the CPU chipset. 
     326        preferred_cpu = (platform.vendor.startswith('Intel') 
     327                         or platform.vendor.startswith('Apple')) 
     328        preferred_gpu = (platform.vendor.startswith('Advanced') 
     329                         or platform.vendor.startswith('NVIDIA')) 
    324330        for device in platform.get_devices(): 
    325331            if device.type == cl.device_type.GPU: 
    326                 # If the existing type is not GPU then it will be CUSTOM or ACCELERATOR, 
    327                 # so don't override it. 
     332                # If the existing type is not GPU then it will be CUSTOM 
     333                # or ACCELERATOR so don't override it. 
    328334                if gpu is None or (preferred_gpu and gpu.type == cl.device_type.GPU): 
    329335                    gpu = device 
    … …  
    334340                # System has cl.device_type.ACCELERATOR or cl.device_type.CUSTOM 
    335341                # Intel Phi for example registers as an accelerator 
    336                 # Since the user installed a custom device on their system and went through the 
    337                 # pain of sorting out OpenCL drivers for it, lets assume they really do want to 
    338                 # use it as their primary compute device. 
     342                # Since the user installed a custom device on their system 
     343                # and went through the pain of sorting out OpenCL drivers for 
     344                # it, lets assume they really do want to use it as their 
     345                # primary compute device. 
    339346                gpu = device 
    340347 
    341     # order the devices by gpu then by cpu; when searching for an available device by data type they 
    342     # will be checked in this order, which means that if the gpu supports double then the cpu will never 
    343     # be used (though we may make it possible to explicitly request the cpu at some point). 
     348    # order the devices by gpu then by cpu; when searching for an available 
     349    # device by data type they will be checked in this order, which means 
     350    # that if the gpu supports double then the cpu will never be used (though 
     351    # we may make it possible to explicitly request the cpu at some point). 
    344352    devices = [] 
    345353    if gpu is not None: 
    … …  
    372380        self.fast = fast 
    373381        self.program = None # delay program creation 
     382        self._kernels = None 
    374383 
    375384    def __getstate__(self): 
    … …  
    394403            names = [generate.kernel_name(self.info, k) for k in variants] 
    395404            kernels = [getattr(self.program, k) for k in names] 
    396             self._kernels = dict((k,v) for k,v in zip(variants, kernels)) 
     405            self._kernels = dict((k, v) for k, v in zip(variants, kernels)) 
    397406        is_2d = len(q_vectors) == 2 
    398407        if is_2d: 
    … …  
    500509    def __init__(self, kernel, dtype, model_info, q_vectors): 
    501510        # type: (cl.Kernel, np.dtype, ModelInfo, List[np.ndarray]) -> None 
    502         max_pd = model_info.parameters.max_pd 
    503         npars = len(model_info.parameters.kernel_parameters)-2 
    504511        q_input = GpuInput(q_vectors, dtype) 
    505512        self.kernel = kernel 
    … …  
    516523 
    517524        self.result_b = cl.Buffer(self.queue.context, mf.READ_WRITE, 
    518                                q_input.global_size[0] * dtype.itemsize) 
     525                                  q_input.global_size[0] * dtype.itemsize) 
    519526        self.q_input = q_input # allocated by GpuInput above 
    520527 
    521         self._need_release = [ self.result_b, self.q_input ] 
     528        self._need_release = [self.result_b, self.q_input] 
    522529        self.real = (np.float32 if dtype == generate.F32 
    523530                     else np.float64 if dtype == generate.F64 

  • sasmodels/kerneldll.py

    r9eb3632 r739aad4  
    5757import numpy as np  # type: ignore 
    5858 
     59try: 
     60    import tinycc 
     61except ImportError: 
     62    tinycc = None 
     63 
    5964from . import generate 
    6065from .kernel import KernelModel, Kernel 
    … …  
    7075    pass 
    7176 
    72 if os.name == 'nt': 
    73     ARCH = "" if sys.maxint > 2**32 else "x86"  # maxint=2**31-1 on 32 bit 
    74     # Windows compiler; check if TinyCC is available 
    75     try: 
    76         import tinycc 
    77     except ImportError: 
    78         tinycc = None 
    79     # call vcvarsall.bat before compiling to set path, headers, libs, etc. 
     77if "SAS_COMPILER" in os.environ: 
     78    compiler = os.environ["SAS_COMPILER"] 
     79elif os.name == 'nt': 
     80    # If vcvarsall.bat has been called, then VCINSTALLDIR is in the environment 
     81    # and we can use the MSVC compiler.  Otherwise, if tinycc is available 
     82    # the use it.  Otherwise, hope that mingw is available. 
    8083    if "VCINSTALLDIR" in os.environ: 
    81         # MSVC compiler is available, so use it.  OpenMP requires a copy of 
    82         # vcomp90.dll on the path.  One may be found here: 
    83         #       C:/Windows/winsxs/x86_microsoft.vc90.openmp*/vcomp90.dll 
    84         # Copy this to the python directory and uncomment the OpenMP COMPILE 
    85         # TODO: remove intermediate OBJ file created in the directory 
    86         # TODO: maybe don't use randomized name for the c file 
    87         # TODO: maybe ask distutils to find MSVC 
    88         CC = "cl /nologo /Ox /MD /W3 /GS- /DNDEBUG".split() 
    89         if "SAS_OPENMP" in os.environ: 
    90             CC.append("/openmp") 
    91         LN = "/link /DLL /INCREMENTAL:NO /MANIFEST".split() 
    92         def compile_command(source, output): 
    93             return CC + ["/Tp%s"%source] + LN + ["/OUT:%s"%output] 
     84        compiler = "msvc" 
    9485    elif tinycc: 
    95         # TinyCC compiler. 
    96         CC = [tinycc.TCC] + "-shared -rdynamic -Wall".split() 
    97         def compile_command(source, output): 
    98             return CC + [source, "-o", output] 
     86        compiler = "tinycc" 
    9987    else: 
    100         # MinGW compiler. 
    101         CC = "gcc -shared -std=c99 -O2 -Wall".split() 
    102         if "SAS_OPENMP" in os.environ: 
    103             CC.append("-fopenmp") 
    104         def compile_command(source, output): 
    105             return CC + [source, "-o", output, "-lm"] 
     88        compiler = "mingw" 
    10689else: 
    107     ARCH = "" 
     90    compiler = "unix" 
     91 
     92ARCH = "" if sys.maxint > 2**32 else "x86"  # maxint=2**31-1 on 32 bit 
     93if compiler == "unix": 
    10894    # Generic unix compile 
    10995    # On mac users will need the X code command line tools installed 
    … …  
    11298    # add openmp support if not running on a mac 
    11399    if sys.platform != "darwin": 
     100        CC.append("-fopenmp") 
     101    def compile_command(source, output): 
     102        return CC + [source, "-o", output, "-lm"] 
     103elif compiler == "msvc": 
     104    # Call vcvarsall.bat before compiling to set path, headers, libs, etc. 
     105    # MSVC compiler is available, so use it.  OpenMP requires a copy of 
     106    # vcomp90.dll on the path.  One may be found here: 
     107    #       C:/Windows/winsxs/x86_microsoft.vc90.openmp*/vcomp90.dll 
     108    # Copy this to the python directory and uncomment the OpenMP COMPILE 
     109    # TODO: remove intermediate OBJ file created in the directory 
     110    # TODO: maybe don't use randomized name for the c file 
     111    # TODO: maybe ask distutils to find MSVC 
     112    CC = "cl /nologo /Ox /MD /W3 /GS- /DNDEBUG".split() 
     113    if "SAS_OPENMP" in os.environ: 
     114        CC.append("/openmp") 
     115    LN = "/link /DLL /INCREMENTAL:NO /MANIFEST".split() 
     116    def compile_command(source, output): 
     117        return CC + ["/Tp%s"%source] + LN + ["/OUT:%s"%output] 
     118elif compiler == "tinycc": 
     119    # TinyCC compiler. 
     120    CC = [tinycc.TCC] + "-shared -rdynamic -Wall".split() 
     121    def compile_command(source, output): 
     122        return CC + [source, "-o", output] 
     123elif compiler == "mingw": 
     124    # MinGW compiler. 
     125    CC = "gcc -shared -std=c99 -O2 -Wall".split() 
     126    if "SAS_OPENMP" in os.environ: 
    114127        CC.append("-fopenmp") 
    115128    def compile_command(source, output): 

  • sasmodels/models/core_shell_ellipsoid.py

    rec45c4f r7f1ee79  
    9999category = "shape:ellipsoid" 
    100100 
    101 single = False  # TODO: maybe using sph_j1c inside gfn would help? 
    102101# pylint: disable=bad-whitespace, line-too-long 
    103102#             ["name", "units", default, [lower, upper], "type", "description"], 

  • sasmodels/models/gaussian_peak.py

    r2c74c11 r7f1ee79  
    3838category = "shape-independent" 
    3939 
    40 single = False 
    4140#             ["name", "units", default, [lower, upper], "type","description"], 
    4241parameters = [["q0", "1/Ang", 0.05, [-inf, inf], "", "Peak position"], 

  • sasmodels/models/hardsphere.py

    r9eb3632 r7f1ee79  
    5858category = "structure-factor" 
    5959structure_factor = True 
    60 single = False 
     60single = False # TODO: check 
    6161 
    6262#             ["name", "units", default, [lower, upper], "type","description"], 

  • sasmodels/models/hollow_cylinder.c

    r2f5c6d4 rf95556f  
    66        double solvent_sld, double theta, double phi); 
    77 
    8 #define INVALID(v) (v.core_radius >= v.radius || v.radius >= v.length) 
     8#define INVALID(v) (v.core_radius >= v.radius) 
    99 
    1010// From Igor library 
    11 static double hollow_cylinder_scaling(double integrand, double delrho, double volume) 
     11static double hollow_cylinder_scaling( 
     12    double integrand, double delrho, double volume) 
    1213{ 
    13         double answer; 
    14         // Multiply by contrast^2 
    15         answer = integrand*delrho*delrho; 
     14    double answer; 
     15    // Multiply by contrast^2 
     16    answer = integrand*delrho*delrho; 
    1617 
    17         //normalize by cylinder volume 
    18         answer *= volume*volume; 
     18    //normalize by cylinder volume 
     19    answer *= volume*volume; 
    1920 
    20         //convert to [cm-1] 
    21         answer *= 1.0e-4; 
     21    //convert to [cm-1] 
     22    answer *= 1.0e-4; 
    2223 
    23         return answer; 
     24    return answer; 
    2425} 
    2526 
    26 static double _hollow_cylinder_kernel(double q, double core_radius, double radius, 
    27         double length, double dum) 
     27 
     28static double _hollow_cylinder_kernel( 
     29    double q, double core_radius, double radius, double length, double dum) 
    2830{ 
    29     double gamma,arg1,arg2,lam1,lam2,psi,sinarg,t2,retval;              //local variables 
    30      
    31     gamma = core_radius/radius; 
    32     arg1 = q*radius*sqrt(1.0-dum*dum);          //1=shell (outer radius) 
    33     arg2 = q*core_radius*sqrt(1.0-dum*dum);                     //2=core (inner radius) 
    34     if (arg1 == 0.0){ 
    35         lam1 = 1.0; 
    36     }else{ 
    37         lam1 = sas_J1c(arg1); 
    38     } 
    39     if (arg2 == 0.0){ 
    40         lam2 = 1.0; 
    41     }else{ 
    42         lam2 = sas_J1c(arg2); 
    43     } 
    44     //Todo: Need to check psi behavior as gamma goes to 1. 
    45     psi = (lam1 -  gamma*gamma*lam2)/(1.0-gamma*gamma);         //SRK 10/19/00 
    46     sinarg = q*length*dum/2.0; 
    47     if (sinarg == 0.0){ 
    48         t2 = 1.0; 
    49     }else{ 
    50         t2 = sin(sinarg)/sinarg; 
    51     } 
     31    const double qs = q*sqrt(1.0-dum*dum); 
     32    const double lam1 = sas_J1c(radius*qs); 
     33    const double lam2 = sas_J1c(core_radius*qs); 
     34    const double gamma_sq = square(core_radius/radius); 
     35    //Note: lim_{r -> r_c} psi = J0(core_radius*qs) 
     36    const double psi = (lam1 - gamma_sq*lam2)/(1.0 - gamma_sq); //SRK 10/19/00 
     37    const double t2 = sinc(q*length*dum/2.0); 
     38    return square(psi*t2); 
     39} 
    5240 
    53     retval = psi*psi*t2*t2; 
    54      
    55     return(retval); 
    56 } 
    57 static double hollow_cylinder_analytical_2D_scaled(double q, double q_x, double q_y, double radius, double core_radius, double length, double sld, 
    58         double solvent_sld, double theta, double phi) { 
    59         double cyl_x, cyl_y; //, cyl_z 
    60         //double q_z; 
    61         double vol, cos_val, delrho; 
    62         double answer; 
    63         //convert angle degree to radian 
    64         double pi = 4.0*atan(1.0); 
    65         theta = theta * pi/180.0; 
    66         phi = phi * pi/180.0; 
    67         delrho = solvent_sld - sld; 
    6841 
    69         // Cylinder orientation 
    70         cyl_x = cos(theta) * cos(phi); 
    71         cyl_y = sin(theta); 
    72         //cyl_z = -cos(theta) * sin(phi); 
     42static double hollow_cylinder_analytical_2D_scaled( 
     43    double q, double q_x, double q_y, double radius, double core_radius, 
     44    double length, double sld, double solvent_sld, double theta, double phi) 
     45{ 
     46    double cyl_x, cyl_y; //, cyl_z 
     47    //double q_z; 
     48    double vol, cos_val, delrho; 
     49    double answer; 
     50    //convert angle degree to radian 
     51    theta = theta * M_PI_180; 
     52    phi = phi * M_PI_180; 
     53    delrho = solvent_sld - sld; 
    7354 
    74         // q vector 
    75         //q_z = 0; 
     55    // Cylinder orientation 
     56    cyl_x = cos(theta) * cos(phi); 
     57    cyl_y = sin(theta); 
     58    //cyl_z = -cos(theta) * sin(phi); 
    7659 
    77         // Compute the angle btw vector q and the 
    78         // axis of the cylinder 
    79         cos_val = cyl_x*q_x + cyl_y*q_y;// + cyl_z*q_z; 
     60    // q vector 
     61    //q_z = 0; 
    8062 
    81         answer = _hollow_cylinder_kernel(q, core_radius, radius, length, cos_val); 
     63    // Compute the angle btw vector q and the 
     64    // axis of the cylinder 
     65    cos_val = cyl_x*q_x + cyl_y*q_y;// + cyl_z*q_z; 
    8266 
    83         vol = form_volume(radius, core_radius, length); 
    84         answer = hollow_cylinder_scaling(answer, delrho, vol); 
     67    answer = _hollow_cylinder_kernel(q, core_radius, radius, length, cos_val); 
    8568 
    86         return answer; 
     69    vol = form_volume(radius, core_radius, length); 
     70    answer = hollow_cylinder_scaling(answer, delrho, vol); 
     71 
     72    return answer; 
    8773} 
    8874 
    … …  
    9076double form_volume(double radius, double core_radius, double length) 
    9177{ 
    92         double pi = 4.0*atan(1.0); 
    93         double v_shell = pi*length*(radius*radius-core_radius*core_radius); 
    94         return(v_shell); 
     78    double v_shell = M_PI*length*(radius*radius-core_radius*core_radius); 
     79    return(v_shell); 
    9580} 
    9681 
    9782 
    98 double Iq(double q, double radius, double core_radius, double length, double sld, 
    99         double solvent_sld) 
     83double Iq(double q, double radius, double core_radius, double length, 
     84    double sld, double solvent_sld) 
    10085{ 
    10186    int i; 
    102         int nord=76;                    //order of integration 
    103         double lower,upper,zi, inter;           //upper and lower integration limits 
    104         double summ,answer,delrho;                      //running tally of integration 
    105         double norm,volume;     //final calculation variables 
    106          
    107         delrho = solvent_sld - sld; 
    108         lower = 0.0; 
    109         upper = 1.0;            //limits of numerical integral 
     87    double lower,upper,zi, inter;               //upper and lower integration limits 
     88    double summ,answer,delrho;                  //running tally of integration 
     89    double norm,volume; //final calculation variables 
    11090 
    111         summ = 0.0;                     //initialize intergral 
    112         for(i=0;i<nord;i++) { 
    113                 zi = ( Gauss76Z[i] * (upper-lower) + lower + upper )/2.0; 
    114                 inter = Gauss76Wt[i] * _hollow_cylinder_kernel(q, core_radius, radius, length, zi); 
    115                 summ += inter; 
    116         } 
    117          
    118         norm = summ*(upper-lower)/2.0; 
    119         volume = form_volume(radius, core_radius, length); 
    120         answer = hollow_cylinder_scaling(norm, delrho, volume); 
    121          
    122         return(answer); 
     91    lower = 0.0; 
     92    upper = 1.0;                //limits of numerical integral 
     93 
     94    summ = 0.0;                 //initialize intergral 
     95    for (i=0;i<76;i++) { 
     96        zi = ( Gauss76Z[i] * (upper-lower) + lower + upper )/2.0; 
     97        inter = Gauss76Wt[i] * _hollow_cylinder_kernel(q, core_radius, radius, length, zi); 
     98        summ += inter; 
     99    } 
     100 
     101    norm = summ*(upper-lower)/2.0; 
     102    volume = form_volume(radius, core_radius, length); 
     103    delrho = solvent_sld - sld; 
     104    answer = hollow_cylinder_scaling(norm, delrho, volume); 
     105 
     106    return(answer); 
    123107} 
    124108 
    125 double Iqxy(double qx, double qy, double radius, double core_radius, double length, double sld, 
    126         double solvent_sld, double theta, double phi) 
     109 
     110double Iqxy(double qx, double qy, double radius, double core_radius, 
     111    double length, double sld, double solvent_sld, double theta, double phi) 
    127112{ 
    128         double q; 
    129         q = sqrt(qx*qx+qy*qy); 
    130         return hollow_cylinder_analytical_2D_scaled(q, qx/q, qy/q, radius, core_radius, length, sld, solvent_sld, theta, phi); 
     113    const double q = sqrt(qx*qx+qy*qy); 
     114    return hollow_cylinder_analytical_2D_scaled(q, qx/q, qy/q, radius, core_radius, length, sld, solvent_sld, theta, phi); 
    131115} 

  • sasmodels/models/hollow_cylinder.py

    r42356c8 r58210db  
    7979# pylint: enable=bad-whitespace, line-too-long 
    8080 
    81 source = ["lib/polevl.c","lib/sas_J1.c", "lib/gauss76.c", "hollow_cylinder.c"] 
     81source = ["lib/polevl.c", "lib/sas_J1.c", "lib/gauss76.c", "hollow_cylinder.c"] 
    8282 
    8383# pylint: disable=W0613 

  • sasmodels/models/lamellar_hg_stack_caille.py

    r42356c8 r7f1ee79  
    9090category = "shape:lamellae" 
    9191 
    92 single = False 
     92single = False  # TODO: check 
    9393parameters = [ 
    9494    #   [ "name", "units", default, [lower, upper], "type", 

  • sasmodels/models/lamellar_stack_caille.py

    r42356c8 r7f1ee79  
    8383category = "shape:lamellae" 
    8484 
    85 single = False 
     85single = False  # TODO: check 
    8686# pylint: disable=bad-whitespace, line-too-long 
    8787#             ["name", "units", default, [lower, upper], "type","description"], 

  • sasmodels/models/linear_pearls.py

    r42356c8 rd2d6100  
    6161    ] 
    6262# pylint: enable=bad-whitespace, line-too-long 
     63single=False 
    6364 
    6465source = ["linear_pearls.c"] 

  • sasmodels/models/multilayer_vesicle.py

    r42356c8 rf67f26c  
    7474    ["volfraction", "",  0.05, [0.0, 1],  "", "volume fraction of vesicles"], 
    7575    ["radius", "Ang", 60.0, [0.0, inf],  "", "Core radius of the multishell"], 
    76     ["thick_shell", "Ang",        10.0, [0.0, inf],  "sld", "Shell thickness"], 
     76    ["thick_shell", "Ang",        10.0, [0.0, inf],  "", "Shell thickness"], 
    7777    ["thick_solvent", "Ang",        10.0, [0.0, inf],  "", "Water thickness"], 
    7878    ["sld_solvent",    "1e-6/Ang^2",  6.4, [-inf, inf], "sld", "Core scattering length density"], 

  • sasmodels/models/onion.c

    r639c4e3 rd119f34  
    22static double 
    33f_exp(double q, double r, double sld_in, double sld_out, 
    4     double thickness, double A) 
    5 { 
    6   const double vol = M_4PI_3 * cube(r); 
    7   const double qr = q * r; 
    8   const double alpha = A * r/thickness; 
    9   const double bes = sph_j1c(qr); 
    10   const double B = (sld_out - sld_in)/expm1(A); 
    11   const double C = sld_in - B; 
    12   double fun; 
    13   if (qr == 0.0) { 
    14     fun = 1.0; 
    15   } else if (fabs(A) > 0.0) { 
    16     const double qrsq = qr*qr; 
    17     const double alphasq = alpha*alpha; 
    18     const double sumsq = alphasq + qrsq; 
    19     double sinqr, cosqr; 
    20     SINCOS(qr, sinqr, cosqr); 
    21     fun = -3.0*( 
    22             ((alphasq - qrsq)*sinqr/qr - 2.0*alpha*cosqr) / sumsq 
    23                 - (alpha*sinqr/qr - cosqr) 
    24         ) / sumsq; 
    25   } else { 
    26     fun = bes; 
    27   } 
    28   return vol * (B*fun + C*bes); 
    29 } 
    30  
    31 static double 
    32 f_linear(double q, double r, double sld, double slope) 
     4    double thickness, double A, double side) 
    335{ 
    346  const double vol = M_4PI_3 * cube(r); 
    357  const double qr = q * r; 
    368  const double bes = sph_j1c(qr); 
    37   double fun = 0.0; 
    38   if (qr > 0.0) { 
    39     const double qrsq = qr*qr; 
     9  const double alpha = A * r/thickness; 
     10  double result; 
     11  if (qr == 0.0) { 
     12    result = 1.0; 
     13  } else if (fabs(A) > 0.0) { 
     14    const double qrsq = qr * qr; 
     15    const double alphasq = alpha * alpha; 
     16    const double sumsq = alphasq + qrsq; 
    4017    double sinqr, cosqr; 
    4118    SINCOS(qr, sinqr, cosqr); 
    42     // Jae-He's code seems to simplify to this 
    43     //     fun = 3.0 * slope * r * (2.0*qr*sinqr - (qrsq-2.0)*cosqr)/(qrsq*qrsq); 
    44     // Rederiving the math, we get the following instead: 
    45     fun = 3.0 * slope * r * (2.0*cosqr + qr*sinqr)/(qrsq*qrsq); 
     19    const double t1 = (alphasq - qrsq)*sinqr/qr - 2.0*alpha*cosqr; 
     20    const double t2 = alpha*sinqr/qr - cosqr; 
     21    const double fun = -3.0*(t1/sumsq - t2)/sumsq; 
     22    const double slope = (sld_out - sld_in)/expm1(A); 
     23    const double contrast = slope*exp(A*side); 
     24    result = contrast*fun + (sld_in-slope)*bes; 
     25  } else { 
     26    result = sld_in*bes; 
    4627  } 
    47   return vol * (sld*bes + fun); 
    48 } 
    49  
    50 static double 
    51 f_constant(double q, double r, double sld) 
    52 { 
    53   const double bes = sph_j1c(q * r); 
    54   const double vol = M_4PI_3 * cube(r); 
    55   return sld * vol * bes; 
     28  return vol * result; 
    5629} 
    5730 
    … …  
    6942static double 
    7043Iq(double q, double sld_core, double core_radius, double sld_solvent, 
    71     double n, double sld_in[], double sld_out[], double thickness[], 
     44    double n_shells, double sld_in[], double sld_out[], double thickness[], 
    7245    double A[]) 
    7346{ 
    74   int i; 
    75   double r = core_radius; 
    76   double f = f_constant(q, r, sld_core); 
    77   for (i=0; i<n; i++){ 
    78     const double r0 = r; 
    79     r += thickness[i]; 
    80     if (r == r0) { 
    81       // no thickness, so nothing to add 
    82     } else if (fabs(A[i]) < 1.0e-16 || sld_out[i] == sld_in[i]) { 
    83       f -= f_constant(q, r0, sld_in[i]); 
    84       f += f_constant(q, r, sld_in[i]); 
    85     } else if (fabs(A[i]) < 1.0e-4) { 
    86       const double slope = (sld_out[i] - sld_in[i])/thickness[i]; 
    87       f -= f_linear(q, r0, sld_in[i], slope); 
    88       f += f_linear(q, r, sld_out[i], slope); 
    89     } else { 
    90       f -= f_exp(q, r0, sld_in[i], sld_out[i], thickness[i], A[i]); 
    91       f += f_exp(q, r, sld_in[i], sld_out[i], thickness[i], A[i]); 
    92     } 
     47  int n = (int)(n_shells+0.5); 
     48  double r_out = core_radius; 
     49  double f = f_exp(q, r_out, sld_core, 0.0, 0.0, 0.0, 0.0); 
     50  for (int i=0; i < n; i++){ 
     51    const double r_in = r_out; 
     52    r_out += thickness[i]; 
     53    f -= f_exp(q, r_in, sld_in[i], sld_out[i], thickness[i], A[i], 0.0); 
     54    f += f_exp(q, r_out, sld_in[i], sld_out[i], thickness[i], A[i], 1.0); 
    9355  } 
    94   f -= f_constant(q, r, sld_solvent); 
     56  f -= f_exp(q, r_out, sld_solvent, 0.0, 0.0, 0.0, 0.0); 
    9557  const double f2 = f * f * 1.0e-4; 
    9658 

  • sasmodels/models/onion.py

    r63c6a08 rd119f34  
    394394    # Could also specify them individually as 
    395395    # "A1": 0, "A2": -1, "A3": 1e-4, "A4": 1, 
     396    #"core_radius_pd_n": 10, 
     397    #"core_radius_pd": 0.4, 
     398    #"thickness4_pd_n": 10, 
     399    #"thickness4_pd": 0.4, 
    396400    } 
    397401 

  • sasmodels/models/polymer_micelle.py

    r42356c8 rd2d6100  
    5151# pylint: enable=bad-whitespace, line-too-long 
    5252 
     53single = False 
     54 
    5355source = ["lib/sph_j1c.c", "polymer_micelle_kernel.c"] 
    5456 

  • sasmodels/models/polymer_micelle_kernel.c

    r2c74c11 rc915373  
    3636    // Self-correlation term of the chains 
    3737    const double qrg2 = q*radius_gyr*q*radius_gyr; 
    38     const double debye_chain = (qrg2 == 0.0) ? 1.0 : 2.0*(exp(-qrg2)-1+qrg2)/(qrg2*qrg2); 
     38    const double debye_chain = (qrg2 == 0.0) ? 1.0 : 2.0*(expm1(-qrg2)+qrg2)/(qrg2*qrg2); 
    3939    const double term2 = n_aggreg * beta_corona * beta_corona * debye_chain; 
    4040 
    4141    // Interference cross-term between core and chains 
    42     const double chain_ampl = (qrg2 == 0.0) ? 1.0 : (1-exp(-qrg2))/qrg2; 
     42    const double chain_ampl = (qrg2 == 0.0) ? 1.0 : -expm1(-qrg2)/qrg2; 
    4343    const double bes_corona = sinc(q*(radius_core + d_penetration * radius_gyr)); 
    4444    const double term3 = 2 * n_aggreg * n_aggreg * beta_core * beta_corona * 

  • sasmodels/models/rpa.c

    rd680a2b r6351bfa  
    2525  double S0ba,Pbb,S0bb,Pbc,S0bc,Pbd,S0bd; 
    2626  double S0ca,S0cb,Pcc,S0cc,Pcd,S0cd; 
    27   double S0da,S0db,S0dc,Pdd,S0dd; 
    28   double Kaa,Kbb,Kcc,Kdd; 
    29   double Kba,Kca,Kda,Kcb,Kdb,Kdc; 
     27  double S0da,S0db,S0dc; 
     28  double Pdd,S0dd; 
     29  double Kaa,Kbb,Kcc; 
     30  double Kba,Kca,Kcb; 
     31  double Kda,Kdb,Kdc,Kdd; 
    3032  double Zaa,Zab,Zac,Zba,Zbb,Zbc,Zca,Zcb,Zcc; 
    3133  double DenT,T11,T12,T13,T21,T22,T23,T31,T32,T33; 

  • sasmodels/models/surface_fractal.py

    rec45c4f r33875e3  
    8080parameters = [["radius",        "Ang", 10.0, [0, inf],   "", 
    8181               "Particle radius"], 
    82               ["surface_dim",   "",    2.0,  [0, inf],   "", 
     82              ["surface_dim",   "",    2.0,  [1, 3],   "", 
    8383               "Surface fractal dimension"], 
    8484              ["cutoff_length", "Ang", 500., [0.0, inf], "", 

  • sasmodels/models/unified_power_Rg.py

    r2c74c11 rec77322  
    8080 
    8181    with errstate(divide='ignore', invalid='ignore'): 
    82         result = np.empty(q.shape, 'd') 
     82        result = np.zeros(q.shape, 'd') 
    8383        for i in range(ilevel): 
    8484            exp_now = exp(-(q*rg[i])**2/3.) 

  • sasmodels/sasview_model.py

    r9eb3632 r4edec6f  
    514514        pairs = [self._get_weights(p) for p in parameters.call_parameters] 
    515515        call_details, values, is_magnetic = build_details(calculator, pairs) 
     516        #call_details.show() 
     517        #print("pairs", pairs) 
     518        #print("params", self.params) 
     519        #print("values", values) 
     520        #print("is_mag", is_magnetic) 
    516521        result = calculator(call_details, values, cutoff=self.cutoff, 
    517522                            magnetic=is_magnetic) 
    … …  
    598603        if par.name not in self.params: 
    599604            if par.name == self.multiplicity_info.control: 
    600                 return [self.multiplicity], [] 
     605                return [self.multiplicity], [1.0] 
    601606            else: 
    602                 return [np.NaN], [] 
     607                return [np.NaN], [1.0] 
    603608        elif par.polydisperse: 
    604609            dis = self.dispersion[par.name] 
    … …  
    608613            return value, weight / np.sum(weight) 
    609614        else: 
    610             return [self.params[par.name]], [] 
     615            return [self.params[par.name]], [1.0] 
    611616 
    612617def test_model(): 

  • sasmodels/kernel_header.c

    r1557a1e redf06e1  
    6161         #define NEED_EXPM1 
    6262         #define NEED_TGAMMA 
     63         // expf missing from windows? 
     64         #define expf exp 
    6365     #else 
    6466         #include <tgmath.h> // C99 type-generic math, so sin(float) => sinf 

  • sasmodels/models/lib/sas_erf.c

    r1557a1e redf06e1  
    280280    else 
    281281        x = a;*/ 
    282  
    283     x = fabsf(a); 
     282    // TODO: tinycc does not support fabsf 
     283    x = fabs(a); 
    284284 
    285285 
    … …  
    302302            return (0.0); 
    303303    } 
    304  
    305304    z = expf(z); 
    306305 
    … …  
    332331    float y, z; 
    333332 
    334     if (fabsf(x) > 1.0) 
     333    // TODO: tinycc does not support fabsf 
     334    if (fabs(x) > 1.0) 
    335335        return (1.0 - erfcf(x)); 
    336336 

  • sasmodels/models/spherical_sld.py

    r63c6a08 redf06e1  
    217217              ] 
    218218# pylint: enable=bad-whitespace, line-too-long 
    219 source = ["lib/sas_erf.c", "lib/librefl.c",  "lib/sph_j1c.c", "spherical_sld.c"] 
     219source = ["lib/polevl.c", "lib/sas_erf.c", "lib/librefl.c",  "lib/sph_j1c.c", "spherical_sld.c"] 
    220220single = False  # TODO: fix low q behaviour 
    221221 
Note: See TracChangeset for help on using the changeset viewer.