Changeset 7e923c2 in sasmodels for 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

Location:
sasmodels
Files:
1 added
18 edited

Legend:

Unmodified
Added
Removed

  • 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.