Changeset 7e923c2 in sasmodels


Ignore:
Timestamp:
Aug 7, 2018 10:45:52 PM (6 years ago)
Author:
Paul Kienzle <pkienzle@…>
Branches:
master, core_shell_microgels, magnetic_model, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests
Children:
1a7ddc9
Parents:
c036ddb (diff), e9b17b18 (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 branch 'master' into beta_approx

Files:
46 added
20 edited

Legend:

Unmodified
Added
Removed

  • doc/guide/gpu_setup.rst

    r59485a4 r63602b1  
    139139the compiler. 
    140140 
    141 On Windows, set *SASCOMPILER=tinycc* for the tinycc compiler, 
    142 *SASCOMPILER=msvc* for the Microsoft Visual C compiler, 
    143 or *SASCOMPILER=mingw* for the MinGW compiler. If TinyCC is available 
     141On Windows, set *SAS_COMPILER=tinycc* for the tinycc compiler, 
     142*SAS_COMPILER=msvc* for the Microsoft Visual C compiler, 
     143or *SAS_COMPILER=mingw* for the MinGW compiler. If TinyCC is available 
    144144on the python path (it is provided with SasView), that will be the 
    145145default. If you want one of the other compilers, be sure to have it 

  • example/model_ellipsoid_hayter_msa.py

    r8a5f021 r9a99993  
    1616 
    1717# DEFINE THE MODEL 
    18 kernel = load_model('ellipsoid*hayter_msa') 
     18kernel = load_model('ellipsoid@hayter_msa') 
    1919 
    2020pars = dict(scale=6.4, background=0.06, sld=0.33, sld_solvent=2.15, radius_polar=14.0, 

  • sasmodels/compare.py

    r1e7b202a raf7a97c  
    112112    -edit starts the parameter explorer 
    113113    -help/-html shows the model docs instead of running the model 
     114 
     115    === environment variables === 
     116    -DSAS_MODELPATH=path sets directory containing custom models 
     117    -DSAS_OPENCL=vendor:device|none sets the target OpenCL device 
     118    -DXDG_CACHE_HOME=~/.cache sets the pyopencl cache root (linux only) 
     119    -DSAS_COMPILER=tinycc|msvc|mingw|unix sets the DLL compiler 
     120    -DSAS_OPENMP=1 turns on OpenMP for the DLLs 
     121    -DSAS_DLL_PATH=path sets the path to the compiled modules 
    114122 
    115123The interpretation of quad precision depends on architecture, and may 
    … …  
    11131121 
    11141122    invalid = [o[1:] for o in flags 
    1115                if o[1:] not in NAME_OPTIONS 
    1116                and not any(o.startswith('-%s='%t) for t in VALUE_OPTIONS)] 
     1123               if not (o[1:] in NAME_OPTIONS 
     1124                       or any(o.startswith('-%s='%t) for t in VALUE_OPTIONS) 
     1125                       or o.startswith('-D'))] 
    11171126    if invalid: 
    11181127        print("Invalid options: %s"%(", ".join(invalid))) 
    … …  
    12151224        elif arg == '-html':    opts['html'] = True 
    12161225        elif arg == '-help':    opts['html'] = True 
     1226        elif arg.startswith('-D'): 
     1227            var, val = arg[2:].split('=') 
     1228            os.environ[var] = val 
    12171229    # pylint: enable=bad-whitespace,C0321 
    12181230 

  • sasmodels/core.py

    rc036ddb r7e923c2  
    3737if CUSTOM_MODEL_PATH == "": 
    3838    CUSTOM_MODEL_PATH = joinpath(os.path.expanduser("~"), ".sasmodels", "custom_models") 
    39     if not os.path.isdir(CUSTOM_MODEL_PATH): 
    40         os.makedirs(CUSTOM_MODEL_PATH) 
     39    #if not os.path.isdir(CUSTOM_MODEL_PATH): 
     40    #    os.makedirs(CUSTOM_MODEL_PATH) 
    4141 
    4242# TODO: refactor composite model support 
    … …  
    232232        if not callable(model_info.Iq): 
    233233            source = generate.make_source(model_info)['dll'] 
    234             old_path = kerneldll.DLL_PATH 
     234            old_path = kerneldll.SAS_DLL_PATH 
    235235            try: 
    236                 kerneldll.DLL_PATH = path 
     236                kerneldll.SAS_DLL_PATH = path 
    237237                dll = kerneldll.make_dll(source, model_info, dtype=numpy_dtype) 
    238238            finally: 
    239                 kerneldll.DLL_PATH = old_path 
     239                kerneldll.SAS_DLL_PATH = old_path 
    240240            compiled_dlls.append(dll) 
    241241    return compiled_dlls 

  • sasmodels/data.py

    rc036ddb r7e923c2  
    502502            # Note: masks merge, so any masked theory points will stay masked, 
    503503            # and the data mask will be added to it. 
    504             mtheory = masked_array(theory, data.mask.copy()) 
     504            #mtheory = masked_array(theory, data.mask.copy()) 
     505            theory_x = data.x[data.mask == 0] 
     506            mtheory = masked_array(theory) 
    505507            mtheory[~np.isfinite(mtheory)] = masked 
    506508            if view is 'log': 
    507509                mtheory[mtheory <= 0] = masked 
    508             plt.plot(data.x, scale*mtheory, '-') 
     510            plt.plot(theory_x, scale*mtheory, '-') 
    509511            all_positive = all_positive and (mtheory > 0).all() 
    510512            some_present = some_present or (mtheory.count() > 0) 
    … …  
    542544 
    543545    if use_resid: 
    544         mresid = masked_array(resid, data.mask.copy()) 
     546        theory_x = data.x[data.mask == 0] 
     547        mresid = masked_array(resid) 
    545548        mresid[~np.isfinite(mresid)] = masked 
    546549        some_present = (mresid.count() > 0) 
    … …  
    548551        if num_plots > 1: 
    549552            plt.subplot(1, num_plots, use_calc + 2) 
    550         plt.plot(data.x, mresid, '.') 
     553        plt.plot(theory_x, mresid, '.') 
    551554        plt.xlabel("$q$/A$^{-1}$") 
    552555        plt.ylabel('residuals') 

  • sasmodels/direct_model.py

    rd18d6dd r1a8c11c  
    250250            qmax = getattr(data, 'qmax', np.inf) 
    251251            accuracy = getattr(data, 'accuracy', 'Low') 
    252             index = ~data.mask & (q >= qmin) & (q <= qmax) 
     252            index = (data.mask == 0) & (q >= qmin) & (q <= qmax) 
    253253            if data.data is not None: 
    254254                index &= ~np.isnan(data.data) 
    … …  
    263263        elif self.data_type == 'Iq': 
    264264            index = (data.x >= data.qmin) & (data.x <= data.qmax) 
     265            mask = getattr(data, 'mask', None) 
     266            if mask is not None: 
     267                index &= (mask == 0) 
    265268            if data.y is not None: 
    266269                index &= ~np.isnan(data.y) 

  • sasmodels/kerneldll.py

    rc036ddb r7e923c2  
    9999    pass 
    100100# pylint: enable=unused-import 
     101 
     102if "SAS_DLL_PATH" in os.environ: 
     103    SAS_DLL_PATH = os.environ["SAS_DLL_PATH"] 
     104else: 
     105    # Assume the default location of module DLLs is in .sasmodels/compiled_models. 
     106    SAS_DLL_PATH = os.path.join(os.path.expanduser("~"), ".sasmodels", "compiled_models") 
    101107 
    102108if "SAS_COMPILER" in os.environ: 
    … …  
    161167        return CC + [source, "-o", output, "-lm"] 
    162168 
    163 # Assume the default location of module DLLs is in .sasmodels/compiled_models. 
    164 DLL_PATH = os.path.join(os.path.expanduser("~"), ".sasmodels", "compiled_models") 
    165  
    166169ALLOW_SINGLE_PRECISION_DLLS = True 
    167170 
    … …  
    200203        return path 
    201204 
    202     return joinpath(DLL_PATH, basename) 
     205    return joinpath(SAS_DLL_PATH, basename) 
    203206 
    204207 
    … …  
    209212    exist yet if it hasn't been compiled. 
    210213    """ 
    211     return os.path.join(DLL_PATH, dll_name(model_info, dtype)) 
     214    return os.path.join(SAS_DLL_PATH, dll_name(model_info, dtype)) 
    212215 
    213216 
    … …  
    228231    models are not allowed as DLLs. 
    229232 
    230     Set *sasmodels.kerneldll.DLL_PATH* to the compiled dll output path. 
     233    Set *sasmodels.kerneldll.SAS_DLL_PATH* to the compiled dll output path. 
     234    Alternatively, set the environment variable *SAS_DLL_PATH*. 
    231235    The default is in ~/.sasmodels/compiled_models. 
    232236    """ 
    … …  
    247251    if need_recompile: 
    248252        # Make sure the DLL path exists 
    249         if not os.path.exists(DLL_PATH): 
    250             os.makedirs(DLL_PATH) 
     253        if not os.path.exists(SAS_DLL_PATH): 
     254            os.makedirs(SAS_DLL_PATH) 
    251255        basename = splitext(os.path.basename(dll))[0] + "_" 
    252256        system_fd, filename = tempfile.mkstemp(suffix=".c", prefix=basename) 

  • sasmodels/model_test.py

    • Property mode changed from 100644 to 100755
    r3221de0 r012cd34  
    376376        stream.writeln(traceback.format_exc()) 
    377377        return 
    378     # Run the test suite 
    379     suite.run(result) 
    380  
    381     # Print the failures and errors 
    382     for _, tb in result.errors: 
    383         stream.writeln(tb) 
    384     for _, tb in result.failures: 
    385         stream.writeln(tb) 
    386378 
    387379    # Warn if there are no user defined tests. 
    … …  
    393385    # iterator since we don't have direct access to the list of tests in the 
    394386    # test suite. 
     387    # In Qt5 suite.run() will clear all tests in the suite after running 
     388    # with no way of retaining them for the test below, so let's check 
     389    # for user tests before running the suite. 
    395390    for test in suite: 
    396391        if not test.info.tests: 
    … …  
    399394    else: 
    400395        stream.writeln("Note: no test suite created --- this should never happen") 
     396 
     397    # Run the test suite 
     398    suite.run(result) 
     399 
     400    # Print the failures and errors 
     401    for _, tb in result.errors: 
     402        stream.writeln(tb) 
     403    for _, tb in result.failures: 
     404        stream.writeln(tb) 
    401405 
    402406    output = stream.getvalue() 

  • sasmodels/modelinfo.py

    rc036ddb r7e923c2  
    585585                Parameter('up:frac_f', '', 0., [0., 1.], 
    586586                          'magnetic', 'fraction of spin up final'), 
    587                 Parameter('up:angle', 'degress', 0., [0., 360.], 
     587                Parameter('up:angle', 'degrees', 0., [0., 360.], 
    588588                          'magnetic', 'spin up angle'), 
    589589            ]) 

  • sasmodels/models/guinier.py

    r2d81cfe rc9fc873  
    77.. math:: 
    88 
    9     I(q) = \text{scale} \cdot \exp{\left[ \frac{-Q^2R_g^2}{3} \right]} 
     9    I(q) = \text{scale} \cdot \exp{\left[ \frac{-Q^2 R_g^2 }{3} \right]} 
    1010            + \text{background} 
    1111 
    … …  
    1919 
    2020.. math:: q = \sqrt{q_x^2 + q_y^2} 
     21 
     22In scattering, the radius of gyration $R_g$ quantifies the objects's 
     23distribution of SLD (not mass density, as in mechanics) from the objects's 
     24SLD centre of mass. It is defined by 
     25 
     26.. math:: R_g^2 = \frac{\sum_i\rho_i\left(r_i-r_0\right)^2}{\sum_i\rho_i} 
     27 
     28where $r_0$ denotes the object's SLD centre of mass and $\rho_i$ is the SLD at 
     29a point $i$. 
     30 
     31Notice that $R_g^2$ may be negative (since SLD can be negative), which happens 
     32when a form factor $P(Q)$ is increasing with $Q$ rather than decreasing. This 
     33can occur for core/shell particles, hollow particles, or for composite 
     34particles with domains of different SLDs in a solvent with an SLD close to the 
     35average match point. (Alternatively, this might be regarded as there being an 
     36internal inter-domain "structure factor" within a single particle which gives 
     37rise to a peak in the scattering). 
     38 
     39To specify a negative value of $R_g^2$ in SasView, simply give $R_g$ a negative 
     40value ($R_g^2$ will be evaluated as $R_g |R_g|$). Note that the physical radius  
     41of gyration, of the exterior of the particle, will still be large and positive.  
     42It is only the apparent size from the small $Q$ data that will give a small or  
     43negative value of $R_g^2$. 
    2144 
    2245References 
    … …  
    4265 
    4366#             ["name", "units", default, [lower, upper], "type","description"], 
    44 parameters = [["rg", "Ang", 60.0, [0, inf], "", "Radius of Gyration"]] 
     67parameters = [["rg", "Ang", 60.0, [-inf, inf], "", "Radius of Gyration"]] 
    4568 
    4669Iq = """ 
    47     double exponent = rg*rg*q*q/3.0; 
     70    double exponent = fabs(rg)*rg*q*q/3.0; 
    4871    double value = exp(-exponent); 
    4972    return value; 
    … …  
    6689 
    6790# parameters for demo 
    68 demo = dict(scale=1.0, rg=60.0) 
     91demo = dict(scale=1.0,  background=0.001, rg=60.0 ) 
    6992 
    7093# parameters for unit tests 

  • sasmodels/resolution.py

    r0b9c6df r9e7837a  
    2020MINIMUM_RESOLUTION = 1e-8 
    2121MINIMUM_ABSOLUTE_Q = 0.02  # relative to the minimum q in the data 
    22 PINHOLE_N_SIGMA = 2.5 # From: Barker & Pedersen 1995 JAC 
     22# According to (Barker & Pedersen 1995 JAC), 2.5 sigma is a good limit. 
     23# According to simulations with github.com:scattering/sansresolution.git 
     24# it is better to use asymmetric bounds (2.5, 3.0) 
     25PINHOLE_N_SIGMA = (2.5, 3.0) 
    2326 
    2427class Resolution(object): 
    … …  
    9093        # from the geometry, they may appear since we are using a truncated 
    9194        # gaussian to represent resolution rather than a skew distribution. 
    92         cutoff = MINIMUM_ABSOLUTE_Q*np.min(self.q) 
    93         self.q_calc = self.q_calc[self.q_calc >= cutoff] 
     95        #cutoff = MINIMUM_ABSOLUTE_Q*np.min(self.q) 
     96        #self.q_calc = self.q_calc[self.q_calc >= cutoff] 
    9497 
    9598        # Build weight matrix from calculated q values 
    … …  
    188191    cdf = erf((edges[:, None] - q[None, :]) / (sqrt(2.0)*q_width)[None, :]) 
    189192    weights = cdf[1:] - cdf[:-1] 
    190     # Limit q range to +/- 2.5 sigma 
    191     qhigh = q + nsigma*q_width 
    192     #qlow = q - nsigma*q_width  # linear limits 
    193     qlow = q*q/qhigh  # log limits 
     193    # Limit q range to (-2.5,+3) sigma 
     194    try: 
     195        nsigma_low, nsigma_high = nsigma 
     196    except TypeError: 
     197        nsigma_low = nsigma_high = nsigma 
     198    qhigh = q + nsigma_high*q_width 
     199    qlow = q - nsigma_low*q_width  # linear limits 
     200    ##qlow = q*q/qhigh  # log limits 
    194201    weights[q_calc[:, None] < qlow[None, :]] = 0. 
    195202    weights[q_calc[:, None] > qhigh[None, :]] = 0. 
    … …  
    365372 
    366373 
    367 def pinhole_extend_q(q, q_width, nsigma=3): 
     374def pinhole_extend_q(q, q_width, nsigma=PINHOLE_N_SIGMA): 
    368375    """ 
    369376    Given *q* and *q_width*, find a set of sampling points *q_calc* so 
    … …  
    371378    function. 
    372379    """ 
    373     q_min, q_max = np.min(q - nsigma*q_width), np.max(q + nsigma*q_width) 
     380    try: 
     381        nsigma_low, nsigma_high = nsigma 
     382    except TypeError: 
     383        nsigma_low = nsigma_high = nsigma 
     384    q_min, q_max = np.min(q - nsigma_low*q_width), np.max(q + nsigma_high*q_width) 
    374385    return linear_extrapolation(q, q_min, q_max) 
    375386 

  • sasmodels/generate.py

    rd86f0fc rc036ddb  
    671671 
    672672 
    673 # type in IQXY pattern could be single, float, double, long double, ... 
    674 _IQXY_PATTERN = re.compile(r"(^|\s)double\s+I(?P<mode>q(ab?c|xy))\s*[(]", 
     673_IQXY_PATTERN = re.compile(r"(^|\s)double\s+I(?P<mode>q(ac|abc|xy))\s*[(]", 
    675674                           flags=re.MULTILINE) 
    676675def find_xy_mode(source): 
    … …  
    701700    return 'qa' 
    702701 
     702# Note: not presently used.  Need to know whether Fq is available before 
     703# trying to compile the source.  Leave the code here in case we decide that 
     704# define have_Fq for each form factor is too tedious and error prone. 
     705_FQ_PATTERN = re.compile(r"(^|\s)void\s+Fq[(]", flags=re.MULTILINE) 
     706def contains_Fq(source): 
     707    # type: (List[str]) -> bool 
     708    """ 
     709    Return True if C source defines "void Fq(". 
     710    """ 
     711    for code in source: 
     712        m = _FQ_PATTERN.search(code) 
     713        if m is not None: 
     714            return True 
     715    return False 
    703716 
    704717def _add_source(source, code, path, lineno=1): 
    … …  
    730743    # dispersion.  Need to be careful that necessary parameters are available 
    731744    # for computing volume even if we allow non-disperse volume parameters. 
    732  
    733745    partable = model_info.parameters 
    734746 
    … …  
    743755    for path, code in user_code: 
    744756        _add_source(source, code, path) 
    745  
    746757    if model_info.c_code: 
    747758        _add_source(source, model_info.c_code, model_info.filename, 
    … …  
    751762    q, qx, qy, qab, qa, qb, qc \ 
    752763        = [Parameter(name=v) for v in 'q qx qy qab qa qb qc'.split()] 
     764 
    753765    # Generate form_volume function, etc. from body only 
    754766    if isinstance(model_info.form_volume, str): 
    … …  
    789801    source.append("\\\n".join(p.as_definition() 
    790802                              for p in partable.kernel_parameters)) 
    791  
    792803    # Define the function calls 
    793804    if partable.form_volume_parameters: 
    … …  
    800811        call_volume = "#define CALL_VOLUME(v) 1.0" 
    801812    source.append(call_volume) 
    802  
    803813    model_refs = _call_pars("_v.", partable.iq_parameters) 
    804     pars = ",".join(["_q"] + model_refs) 
    805     call_iq = "#define CALL_IQ(_q, _v) Iq(%s)" % pars 
     814 
     815    if model_info.have_Fq: 
     816        pars = ",".join(["_q", "&_F1", "&_F2",] + model_refs) 
     817        call_iq = "#define CALL_FQ(_q, _F1, _F2, _v) Fq(%s)" % pars 
     818        clear_iq = "#undef CALL_FQ" 
     819    else: 
     820        pars = ",".join(["_q"] + model_refs) 
     821        call_iq = "#define CALL_IQ(_q, _v) Iq(%s)" % pars 
     822        clear_iq = "#undef CALL_IQ" 
    806823    if xy_mode == 'qabc': 
    807824        pars = ",".join(["_qa", "_qb", "_qc"] + model_refs) 
    … …  
    812829        call_iqxy = "#define CALL_IQ_AC(_qa,_qc,_v) Iqac(%s)" % pars 
    813830        clear_iqxy = "#undef CALL_IQ_AC" 
    814     elif xy_mode == 'qa': 
     831    elif xy_mode == 'qa' and not model_info.have_Fq: 
    815832        pars = ",".join(["_qa"] + model_refs) 
    816833        call_iqxy = "#define CALL_IQ_A(_qa,_v) Iq(%s)" % pars 
    817834        clear_iqxy = "#undef CALL_IQ_A" 
     835    elif xy_mode == 'qa' and model_info.have_Fq: 
     836        pars = ",".join(["_qa", "&_F1", "&_F2",] + model_refs) 
     837        # Note: uses rare C construction (expr1, expr2) which computes 
     838        # expr1 then expr2 and evaluates to expr2.  This allows us to 
     839        # leave it looking like a function even though it is returning 
     840        # its values by reference. 
     841        call_iqxy = "#define CALL_FQ_A(_qa,_F1,_F2,_v) (Fq(%s),_F2)" % pars 
     842        clear_iqxy = "#undef CALL_FQ_A" 
    818843    elif xy_mode == 'qxy': 
    819844        orientation_refs = _call_pars("_v.", partable.orientation_parameters) 
    … …  
    831856    magpars = [k-2 for k, p in enumerate(partable.call_parameters) 
    832857               if p.type == 'sld'] 
    833  
    834858    # Fill in definitions for numbers of parameters 
    835859    source.append("#define MAX_PD %s"%partable.max_pd) 
    … …  
    839863    source.append("#define MAGNETIC_PARS %s"%",".join(str(k) for k in magpars)) 
    840864    source.append("#define PROJECTION %d"%PROJECTION) 
    841  
    842865    # TODO: allow mixed python/opencl kernels? 
    843  
    844     ocl = _kernels(kernel_code, call_iq, call_iqxy, clear_iqxy, model_info.name) 
    845     dll = _kernels(kernel_code, call_iq, call_iqxy, clear_iqxy, model_info.name) 
     866    ocl = _kernels(kernel_code, call_iq, clear_iq, call_iqxy, clear_iqxy, model_info.name) 
     867    dll = _kernels(kernel_code, call_iq, clear_iq, call_iqxy, clear_iqxy, model_info.name) 
     868 
    846869    result = { 
    847870        'dll': '\n'.join(source+dll[0]+dll[1]+dll[2]), 
    848871        'opencl': '\n'.join(source+ocl[0]+ocl[1]+ocl[2]), 
    849872    } 
    850  
    851873    return result 
    852874 
    853875 
    854 def _kernels(kernel, call_iq, call_iqxy, clear_iqxy, name): 
     876def _kernels(kernel, call_iq, clear_iq, call_iqxy, clear_iqxy, name): 
    855877    # type: ([str,str], str, str, str) -> List[str] 
    856878    code = kernel[0] 
    … …  
    862884        '#line 1 "%s Iq"' % path, 
    863885        code, 
    864         "#undef CALL_IQ", 
     886        clear_iq, 
    865887        "#undef KERNEL_NAME", 
    866888        ] 

  • sasmodels/kernel_iq.c

    r7c35fda rc036ddb  
    2929//  CALL_IQ(q, table) : call the Iq function for 1D calcs. 
    3030//  CALL_IQ_A(q, table) : call the Iq function with |q| for 2D data. 
     31//  CALL_FQ(q, F1, F2, table) : call the Fq function for 1D calcs. 
     32//  CALL_FQ_A(q, F1, F2, table) : call the Iq function with |q| for 2D data. 
    3133//  CALL_IQ_AC(qa, qc, table) : call the Iqxy function for symmetric shapes 
    3234//  CALL_IQ_ABC(qa, qc, table) : call the Iqxy function for asymmetric shapes 
    … …  
    3638//  PROJECTION : equirectangular=1, sinusoidal=2 
    3739//      see explore/jitter.py for definitions. 
     40 
    3841 
    3942#ifndef _PAR_BLOCK_ // protected block so we can include this code twice. 
    … …  
    8487  out_spin = clip(out_spin, 0.0, 1.0); 
    8588  // Previous version of this function took the square root of the weights, 
    86   // under the assumption that  
     89  // under the assumption that 
    8790  // 
    8891  //     w*I(q, rho1, rho2, ...) = I(q, sqrt(w)*rho1, sqrt(w)*rho2, ...) 
    … …  
    270273#endif // _QABC_SECTION 
    271274 
    272  
    273275// ==================== KERNEL CODE ======================== 
    274  
     276#define COMPUTE_F1_F2 defined(CALL_FQ) 
    275277kernel 
    276278void KERNEL_NAME( 
    … …  
    341343  // version must loop over all q. 
    342344  #ifdef USE_OPENCL 
    343     double pd_norm = (pd_start == 0 ? 0.0 : result[nq]); 
    344     double this_result = (pd_start == 0 ? 0.0 : result[q_index]); 
     345    #if COMPUTE_F1_F2 
     346      double pd_norm = (pd_start == 0 ? 0.0 : result[nq]); 
     347      double weight_norm = (pd_start == 0 ? 0.0 : result[2*nq+1]); 
     348      double this_F2 = (pd_start == 0 ? 0.0 : result[q_index]); 
     349      double this_F1 = (pd_start == 0 ? 0.0 : result[q_index+nq+1]); 
     350    #else 
     351      double pd_norm = (pd_start == 0 ? 0.0 : result[nq]); 
     352      double this_result = (pd_start == 0 ? 0.0 : result[q_index]); 
     353    #endif 
    345354  #else // !USE_OPENCL 
    346     double pd_norm = (pd_start == 0 ? 0.0 : result[nq]); 
     355    #if COMPUTE_F1_F2 
     356      double pd_norm = (pd_start == 0 ? 0.0 : result[nq]); 
     357      double weight_norm = (pd_start == 0 ? 0.0 : result[2*nq+1]); 
     358    #else 
     359      double pd_norm = (pd_start == 0 ? 0.0 : result[nq]); 
     360    #endif 
    347361    if (pd_start == 0) { 
    348362      #ifdef USE_OPENMP 
    349363      #pragma omp parallel for 
    350364      #endif 
    351       for (int q_index=0; q_index < nq; q_index++) result[q_index] = 0.0; 
     365      #if COMPUTE_F1_F2 
     366          for (int q_index=0; q_index < 2*nq+2; q_index++) result[q_index] = 0.0; 
     367      #else 
     368          for (int q_index=0; q_index < nq; q_index++) result[q_index] = 0.0; 
     369      #endif 
    352370    } 
    353371    //if (q_index==0) printf("start %d %g %g\n", pd_start, pd_norm, result[0]); 
    … …  
    442460// inner loop and defines the macros that use them. 
    443461 
    444 #if defined(CALL_IQ) 
    445   // unoriented 1D 
     462 
     463#if defined(CALL_FQ) // COMPUTE_F1_F2 is true 
     464  // unoriented 1D returning <F> and <F^2> 
     465  double qk; 
     466  double F1, F2; 
     467  #define FETCH_Q() do { qk = q[q_index]; } while (0) 
     468  #define BUILD_ROTATION() do {} while(0) 
     469  #define APPLY_ROTATION() do {} while(0) 
     470  #define CALL_KERNEL() CALL_FQ(qk,F1,F2,local_values.table) 
     471 
     472#elif defined(CALL_FQ_A) 
     473  // unoriented 2D return <F> and <F^2> 
     474  double qx, qy; 
     475  double F1, F2; 
     476  #define FETCH_Q() do { qx = q[2*q_index]; qy = q[2*q_index+1]; } while (0) 
     477  #define BUILD_ROTATION() do {} while(0) 
     478  #define APPLY_ROTATION() do {} while(0) 
     479  #define CALL_KERNEL() CALL_FQ_A(sqrt(qx*qx+qy*qy),F1,F2,local_values.table) 
     480 
     481#elif defined(CALL_IQ) 
     482  // unoriented 1D return <F^2> 
    446483  double qk; 
    447484  #define FETCH_Q() do { qk = q[q_index]; } while (0) 
    448485  #define BUILD_ROTATION() do {} while(0) 
    449486  #define APPLY_ROTATION() do {} while(0) 
    450   #define CALL_KERNEL() CALL_IQ(qk, local_values.table) 
     487  #define CALL_KERNEL() CALL_IQ(qk,local_values.table) 
    451488 
    452489#elif defined(CALL_IQ_A) 
    … …  
    482519  #define APPLY_ROTATION() qabc_apply(&rotation, qx, qy, &qa, &qb, &qc) 
    483520  #define CALL_KERNEL() CALL_IQ_ABC(qa, qb, qc, local_values.table) 
     521 
    484522#elif defined(CALL_IQ_XY) 
    485523  // direct call to qx,qy calculator 
    … …  
    495533// logic in the IQ_AC and IQ_ABC branches.  This will become more important 
    496534// if we implement more projections, or more complicated projections. 
    497 #if defined(CALL_IQ) || defined(CALL_IQ_A)  // no orientation 
     535#if defined(CALL_IQ) || defined(CALL_IQ_A) || defined(CALL_FQ) || defined(CALL_FQ_A) 
     536  // no orientation 
    498537  #define APPLY_PROJECTION() const double weight=weight0 
    499538#elif defined(CALL_IQ_XY) // pass orientation to the model 
    … …  
    646685    if (weight > cutoff) { 
    647686      pd_norm += weight * CALL_VOLUME(local_values.table); 
     687      #if COMPUTE_F1_F2 
     688      weight_norm += weight; 
     689      #endif 
    648690      BUILD_ROTATION(); 
    649691 
    … …  
    693735          } 
    694736        #else  // !MAGNETIC 
    695           const double scattering = CALL_KERNEL(); 
     737          #if COMPUTE_F1_F2 
     738            CALL_KERNEL(); // sets F1 and F2 by reference 
     739          #else 
     740            const double scattering = CALL_KERNEL(); 
     741          #endif 
    696742        #endif // !MAGNETIC 
    697743//printf("q_index:%d %g %g %g %g\n", q_index, scattering, weight0); 
    698744 
    699745        #ifdef USE_OPENCL 
    700           this_result += weight * scattering; 
     746          #if COMPUTE_F1_F2 
     747            this_F1 += weight * F1; 
     748            this_F2 += weight * F2; 
     749          #else 
     750            this_result += weight * scattering; 
     751          #endif 
    701752        #else // !USE_OPENCL 
    702           result[q_index] += weight * scattering; 
     753          #if COMPUTE_F1_F2 
     754            result[q_index] += weight * F2; 
     755            result[q_index+nq+1] += weight * F1; 
     756          #else 
     757            result[q_index] += weight * scattering; 
     758          #endif 
    703759        #endif // !USE_OPENCL 
    704760      } 
    705761    } 
    706762  } 
    707  
    708763// close nested loops 
    709764++step; 
    … …  
    726781// Remember the current result and the updated norm. 
    727782#ifdef USE_OPENCL 
    728   result[q_index] = this_result; 
    729   if (q_index == 0) result[nq] = pd_norm; 
     783  #if COMPUTE_F1_F2 
     784    result[q_index] = this_F2; 
     785    result[q_index+nq+1] = this_F1; 
     786    if (q_index == 0) { 
     787      result[nq] = pd_norm; 
     788      result[2*nq+1] = weight_norm; 
     789    } 
     790  #else 
     791    result[q_index] = this_result; 
     792    if (q_index == 0) result[nq] = pd_norm; 
     793  #endif 
     794 
    730795//if (q_index == 0) printf("res: %g/%g\n", result[0], pd_norm); 
    731796#else // !USE_OPENCL 
    732   result[nq] = pd_norm; 
     797  #if COMPUTE_F1_F2 
     798    result[nq] = pd_norm; 
     799    result[2*nq+1] = weight_norm; 
     800  #else 
     801    result[nq] = pd_norm; 
     802  #endif 
    733803//printf("res: %g/%g\n", result[0], pd_norm); 
    734804#endif // !USE_OPENCL 
    735805 
    736806// ** clear the macros in preparation for the next kernel ** 
     807#undef COMPUTE_F1_F2 
    737808#undef PD_INIT 
    738809#undef PD_OPEN 

  • sasmodels/kernelcl.py

    rd86f0fc rc036ddb  
    538538        self.dtype = dtype 
    539539        self.dim = '2d' if q_input.is_2d else '1d' 
    540         # plus three for the normalization values 
    541         self.result = np.empty(q_input.nq+1, dtype) 
     540        # leave room for f1/f2 results in case we need to compute beta for 1d models 
     541        num_returns = 1 if self.dim == '2d' else 2  # 
     542        # plus 1 for the normalization value 
     543        self.result = np.empty((q_input.nq+1)*num_returns, dtype) 
    542544 
    543545        # Inputs and outputs for each kernel call 
    … …  
    547549 
    548550        self.result_b = cl.Buffer(self.queue.context, mf.READ_WRITE, 
    549                                   q_input.global_size[0] * dtype.itemsize) 
     551                                  q_input.global_size[0] * num_returns * dtype.itemsize) 
    550552        self.q_input = q_input # allocated by GpuInput above 
    551553 
    … …  
    556558                     else np.float32)  # will never get here, so use np.float32 
    557559 
    558     def __call__(self, call_details, values, cutoff, magnetic): 
     560    def Iq(self, call_details, values, cutoff, magnetic): 
     561        # type: (CallDetails, np.ndarray, np.ndarray, float, bool) -> np.ndarray 
     562        self._call_kernel(call_details, values, cutoff, magnetic) 
     563        #print("returned",self.q_input.q, self.result) 
     564        pd_norm = self.result[self.q_input.nq] 
     565        scale = values[0]/(pd_norm if pd_norm != 0.0 else 1.0) 
     566        background = values[1] 
     567        #print("scale",scale,background) 
     568        return scale*self.result[:self.q_input.nq] + background 
     569    __call__ = Iq 
     570 
     571    def beta(self, call_details, values, cutoff, magnetic): 
     572        # type: (CallDetails, np.ndarray, np.ndarray, float, bool) -> np.ndarray 
     573        if self.dim == '2d': 
     574            raise NotImplementedError("beta not yet supported for 2D") 
     575        self._call_kernel(call_details, values, cutoff, magnetic) 
     576        w_norm = self.result[2*self.q_input.nq + 1] 
     577        pd_norm = self.result[self.q_input.nq] 
     578        if w_norm == 0.: 
     579            w_norm = 1. 
     580        F2 = self.result[:self.q_input.nq]/w_norm 
     581        F1 = self.result[self.q_input.nq+1:2*self.q_input.nq+1]/w_norm 
     582        volume_avg = pd_norm/w_norm 
     583        return F1, F2, volume_avg 
     584 
     585    def _call_kernel(self, call_details, values, cutoff, magnetic): 
    559586        # type: (CallDetails, np.ndarray, np.ndarray, float, bool) -> np.ndarray 
    560587        context = self.queue.context 
    … …  
    573600        #print("Calling OpenCL") 
    574601        #call_details.show(values) 
    575         # Call kernel and retrieve results 
     602        #Call kernel and retrieve results 
    576603        wait_for = None 
    577604        last_nap = time.clock() 
    … …  
    598625                v.release() 
    599626 
    600         pd_norm = self.result[self.q_input.nq] 
    601         scale = values[0]/(pd_norm if pd_norm != 0.0 else 1.0) 
    602         background = values[1] 
    603         #print("scale",scale,values[0],self.result[self.q_input.nq],background) 
    604         return scale*self.result[:self.q_input.nq] + background 
    605         # return self.result[:self.q_input.nq] 
    606  
    607627    def release(self): 
    608628        # type: () -> None 

  • sasmodels/models/core_shell_sphere.py

    rdc76240 rc036ddb  
    5858title = "Form factor for a monodisperse spherical particle with particle with a core-shell structure." 
    5959description = """ 
    60     F^2(q) = 3/V_s [V_c (sld_core-sld_shell) (sin(q*radius)-q*radius*cos(q*radius))/(q*radius)^3 
    61                    + V_s (sld_shell-sld_solvent) (sin(q*r_s)-q*r_s*cos(q*r_s))/(q*r_s)^3] 
     60    F(q) = [V_c (sld_core-sld_shell) 3 (sin(q*radius)-q*radius*cos(q*radius))/(q*radius)^3 
     61            + V_s (sld_shell-sld_solvent) 3 (sin(q*r_s)-q*r_s*cos(q*r_s))/(q*r_s)^3] 
    6262 
    6363            V_s: Volume of the sphere shell 

  • sasmodels/models/ellipsoid.c

    r108e70e rc036ddb  
    55} 
    66 
     7/* Fq overrides Iq 
    78static  double 
    89Iq(double q, 
    … …  
    3738    return 1.0e-4 * s * s * form; 
    3839} 
     40*/ 
     41 
     42static void 
     43Fq(double q, 
     44    double *F1, 
     45    double *F2, 
     46    double sld, 
     47    double sld_solvent, 
     48    double radius_polar, 
     49    double radius_equatorial) 
     50{ 
     51    // Using ratio v = Rp/Re, we can implement the form given in Guinier (1955) 
     52    //     i(h) = int_0^pi/2 Phi^2(h a sqrt(cos^2 + v^2 sin^2) cos dT 
     53    //          = int_0^pi/2 Phi^2(h a sqrt((1-sin^2) + v^2 sin^2) cos dT 
     54    //          = int_0^pi/2 Phi^2(h a sqrt(1 + sin^2(v^2-1)) cos dT 
     55    // u-substitution of 
     56    //     u = sin, du = cos dT 
     57    //     i(h) = int_0^1 Phi^2(h a sqrt(1 + u^2(v^2-1)) du 
     58    const double v_square_minus_one = square(radius_polar/radius_equatorial) - 1.0; 
     59    // translate a point in [-1,1] to a point in [0, 1] 
     60    // const double u = GAUSS_Z[i]*(upper-lower)/2 + (upper+lower)/2; 
     61    const double zm = 0.5; 
     62    const double zb = 0.5; 
     63    double total_F2 = 0.0; 
     64    double total_F1 = 0.0; 
     65    for (int i=0;i<GAUSS_N;i++) { 
     66        const double u = GAUSS_Z[i]*zm + zb; 
     67        const double r = radius_equatorial*sqrt(1.0 + u*u*v_square_minus_one); 
     68        const double f = sas_3j1x_x(q*r); 
     69        total_F2 += GAUSS_W[i] * f * f; 
     70        total_F1 += GAUSS_W[i] * f; 
     71    } 
     72    // translate dx in [-1,1] to dx in [lower,upper] 
     73    const double form_squared_avg = total_F2*zm; 
     74    const double form_avg = total_F1*zm; 
     75    const double s = (sld - sld_solvent) * form_volume(radius_polar, radius_equatorial); 
     76    *F1 = 1e-2 * s * form_avg; 
     77    *F2 = 1e-4 * s * s * form_squared_avg; 
     78} 
    3979 
    4080static double 

  • sasmodels/models/ellipsoid.py

    r2d81cfe rc036ddb  
    161161             ] 
    162162 
     163 
    163164source = ["lib/sas_3j1x_x.c", "lib/gauss76.c", "ellipsoid.c"] 
     165 
     166have_Fq = True 
    164167 
    165168def ER(radius_polar, radius_equatorial): 

  • sasmodels/models/lib/sphere_form.c

    r925ad6e r01c8d9e  
    1313    return 1.0e-4*square(contrast * fq); 
    1414} 
     15 

  • sasmodels/models/sphere.py

    ref07e95 rc036ddb  
    6969source = ["lib/sas_3j1x_x.c", "lib/sphere_form.c"] 
    7070 
    71 # No volume normalization despite having a volume parameter 
    72 # This should perhaps be volume normalized? 
    73 form_volume = """ 
     71c_code = """ 
     72static double form_volume(double radius) 
     73{ 
    7474    return sphere_volume(radius); 
    75     """ 
     75} 
    7676 
    77 Iq = """ 
    78     return sphere_form(q, radius, sld, sld_solvent); 
    79     """ 
     77static void Fq(double q, double *F1,double *F2, double sld, double solvent_sld, double radius) 
     78{ 
     79    const double fq = sas_3j1x_x(q*radius); 
     80    const double contrast = (sld - solvent_sld); 
     81    const double form = 1e-2 * contrast * sphere_volume(radius) * fq; 
     82    *F1 = form; 
     83    *F2 = form*form; 
     84} 
     85""" 
     86 
     87# TODO: figure this out by inspection 
     88have_Fq = True 
    8089 
    8190def ER(radius): 

  • sasmodels/product.py

    r2d81cfe rc036ddb  
    1616import numpy as np  # type: ignore 
    1717 
    18 from .modelinfo import ParameterTable, ModelInfo 
     18from .modelinfo import ParameterTable, ModelInfo, Parameter 
    1919from .kernel import KernelModel, Kernel 
    2020from .details import make_details, dispersion_mesh 
    … …  
    3737ER_ID = "radius_effective" 
    3838VF_ID = "volfraction" 
     39BETA_DEFINITION = ("beta_mode", "", 0, [["P*S"],["P*(1+beta*(S-1))"]], "", 
     40                   "Structure factor dispersion calculation mode") 
    3941 
    4042# TODO: core_shell_sphere model has suppressed the volume ratio calculation 
    … …  
    7476    translate_name = dict((old.id, new.id) for old, new 
    7577                          in zip(s_pars.kernel_parameters[1:], s_list)) 
    76     combined_pars = p_pars.kernel_parameters + s_list 
     78    beta = [Parameter(*BETA_DEFINITION)] if p_info.have_Fq else [] 
     79    combined_pars = p_pars.kernel_parameters + s_list + beta 
    7780    parameters = ParameterTable(combined_pars) 
    7881    parameters.max_pd = p_pars.max_pd + s_pars.max_pd 
    … …  
    151154        #: Structure factor modelling interaction between particles. 
    152155        self.S = S 
     156 
    153157        #: Model precision. This is not really relevant, since it is the 
    154158        #: individual P and S models that control the effective dtype, 
    … …  
    168172        # in opencl; or both in opencl, but one in single precision and the 
    169173        # other in double precision). 
     174 
    170175        p_kernel = self.P.make_kernel(q_vectors) 
    171176        s_kernel = self.S.make_kernel(q_vectors) 
    … …  
    211216        p_offset = call_details.offset[:p_npars] 
    212217        p_details = make_details(p_info, p_length, p_offset, nweights) 
     218 
    213219        # Set p scale to the volume fraction in s, which is the first of the 
    214220        # 'S' parameters in the parameter list, or 2+np in 0-origin. 
    … …  
    254260        s_values = np.hstack(s_values).astype(self.s_kernel.dtype) 
    255261 
     262        # beta mode is the first parameter after the structure factor pars 
     263        beta_index = 2+p_npars+s_npars 
     264        beta_mode = values[beta_index] 
     265 
    256266        # Call the kernels 
    257         p_result = self.p_kernel(p_details, p_values, cutoff, magnetic) 
    258         s_result = self.s_kernel(s_details, s_values, cutoff, False) 
    259  
    260         #print("p_npars",p_npars,s_npars,p_er,s_vr,values[2+p_npars+1:2+p_npars+s_npars]) 
     267        s_result = self.s_kernel.Iq(s_details, s_values, cutoff, False) 
     268        scale, background = values[0], values[1] 
     269        if beta_mode: 
     270            F1, F2, volume_avg = self.p_kernel.beta(p_details, p_values, cutoff, magnetic) 
     271            combined_scale = scale*volfrac/volume_avg 
     272            # Define lazy results based on intermediate values. 
     273            # The return value for the calculation should be an ordered 
     274            # dictionary containing any result the user might want to see 
     275            # at the end of the calculation, including scalars, strings, 
     276            # and plottable data.  Don't want to build this structure during 
     277            # fits, only when displaying the final result (or a one-off 
     278            # computation which asks for it). 
     279            # Do not use the current hack of storing the intermediate values 
     280            # in self.results since that leads to awkward threading issues. 
     281            # Instead return the function along with the bundle of inputs. 
     282            # P and Q may themselves have intermediate results they want to 
     283            # include, such as A and B if P = A + B.  Might use this mechanism 
     284            # to return the computed effective radius as well. 
     285            #def lazy_results(Q, S, F1, F2, scale): 
     286            #    """ 
     287            #    beta = F1**2 / F2  # what about divide by zero errors? 
     288            #    return { 
     289            #        'P' : Data1D(Q, scale*F2), 
     290            #        'beta': Data1D(Q, beta), 
     291            #        'S' : Data1D(Q, S), 
     292            #        'Seff': Data1D(Q, 1 + beta*(S-1)), 
     293            #        'I' : Data1D(Q, scale*(F2 + (F1**2)*(S-1)) + background), 
     294            #    } 
     295            #lazy_pars = s_result, F1, F2, combined_scale 
     296            self.results = [F2, s_result] 
     297            final_result = combined_scale*(F2 + (F1**2)*(s_result - 1)) + background 
     298        else: 
     299            p_result = self.p_kernel.Iq(p_details, p_values, cutoff, magnetic) 
     300            # remember the parts for plotting later 
     301            self.results = [p_result, s_result] 
     302            final_result = scale*(p_result*s_result) + background 
     303 
    261304        #call_details.show(values) 
    262305        #print("values", values) 
    … …  
    265308        #s_details.show(s_values) 
    266309        #print("=>", s_result) 
    267  
    268         # remember the parts for plotting later 
    269         self.results = [p_result, s_result] 
    270  
    271310        #import pylab as plt 
    272311        #plt.subplot(211); plt.loglog(self.p_kernel.q_input.q, p_result, '-') 
    273312        #plt.subplot(212); plt.loglog(self.s_kernel.q_input.q, s_result, '-') 
    274313        #plt.figure() 
    275  
    276         return values[0]*(p_result*s_result) + values[1] 
     314        return final_result 
    277315 
    278316    def release(self): 
Note: See TracChangeset for help on using the changeset viewer.