# Home

## 3. Scripting Interface

## 3.1. Preparing your environment

You can use python scripts to load and plot your data, create SAS models and fit parameters. You can save a script to a file such as example/model.py and run it later. However, this requires a scripting environment with the correct packages installed.

You can either use the SasView application itself (versions after 5.0.5), as both bumps and sasmodels are included as part of the distribution, so for example on Windows:

    > sasview model.py

(Note that it may be necessary to first add the folder path to sasmodels/sasview to your Path environment variable for this to work; set PATH=%PATH%;C:\your\path\here\ . The folder path can be found in the Help > About box if you are running the SasView GUI.)

or create a Python environment with pip:

    > pip install sasmodels sasdata matplotlib bumps periodictable
    > python model.py

(You can also create a Python environment using conda, see: https://github.com/SasView/sasview/wiki/DevNotes_CondaDevEnviroment)

The pip command also works within a Jupyter notebook

    %pip install sasmodels sasdata matplotlib bumps periodictable

On a Mac the command for invoking SasView will be something like:

    /Applications/Sasview5.dmg/Contents/MacOS/sasview

depending on where it is actually installed. This can either be used directly or can be symlinked into your path, for example:

    mkdir ~/bin
    ln -s /path/to/Applications/SasView5.dmg/Contents/MacOS/sasview ~/bin

## 3.2. Preparing your data

The key functions are then `core.load_model()` for loading the model definition and compiling the kernel and `data.load_data()` for calling sasview to load the data.

Usually you will load data via the sasview loader, with the `data.load_data()` function. For example:

    from sasmodels.data import load_data
    data = load_data("sasmodels/example/093191_201.dat")

You may want to apply a data mask, such a beam stop, and trim high \(q\):

    from sasmodels.data import set_beam_stop
    set_beam_stop(data, qmin, qmax)

The `data.set_beam_stop()` method simply sets the mask attribute for the data.

The data defines the resolution function and the q values to evaluate, so even if you simulating experiments prior to making measurements, you still need a data object for reference. Use `data.empty_data1D()` or `data.empty_data2D()` to create a container with a given \(q\) and \(\Delta q/q\). For example:

    import numpy as np
    from sasmodels.data import empty_data1D

    # 120 points logarithmically spaced from 0.005 to 0.2, with dq/q = 5%
    q = np.logspace(np.log10(5e-3), np.log10(2e-1), 120)
    data = empty_data1D(q, resolution=0.05)

To use a more realistic model of resolution, or to load data from a file format not understood by SasView, you can use `data.Data1D` or `data.Data2D` directly. The 1D data uses x, y, dx and dy for \(x = q\) and \(y = I(q)\), and 2D data uses x, y, z, dx, dy, dz for \(x, y = qx, qy\) and \(z = I(qx, qy)\). [Note: internally, the Data2D object uses SasView conventions, qx_data, qy_data, data, dqx_data, dqy_data, and err_data.]

For USANS data, use 1D data, but set dxl and dxw attributes to indicate slit resolution:

    data.dxl = 0.117

See `resolution.slit_resolution()` for details.

SESANS data is more complicated; if your SESANS format is not supported by SasView you need to define a number of attributes beyond x, y. For example:

    SElength = np.linspace(0, 2400, 61) # [A]
    data = np.ones_like(SElength)
    err_data = np.ones_like(SElength)*0.03

    class Source:
        wavelength = 6 # [A]
        wavelength_unit = "A"
    class Sample:
        zacceptance = 0.1 # [A^-1]
        thickness = 0.2 # [cm]

    class SESANSData1D:
        #q_zmax = 0.23 # [A^-1]
        lam = 0.2 # [nm]
        x = SElength
        y = data
        dy = err_data
        sample = Sample()
    data = SESANSData1D()

    x, y = ... # create or load sesans
    data = smd.Data

The data module defines various data plotters as well.

## 3.3. Using sasmodels directly

Once you have a computational kernel and a data object, you can evaluate the model for various parameters using `direct_model.DirectModel`. The resulting object f will be callable as f(par=value, …), returning the \(I(q)\) for the \(q\) values in the data. For example:

    import numpy as np
    from sasmodels.data import empty_data1D
    from sasmodels.core import load_model
    from sasmodels.direct_model import DirectModel

    # 120 points logarithmically spaced from 0.005 to 0.2, with dq/q = 5%
    q = np.logspace(np.log10(5e-3), np.log10(2e-1), 120)
    data = empty_data1D(q, resolution=0.05)
    kernel = load_model("ellipsoid)
    f = DirectModel(data, kernel)
    Iq = f(radius_polar=100)

Polydispersity information is set with special parameter names:

> - par_pd for polydispersity width, \(\Delta p/p\),

> - par_pd_n for the number of points in the distribution,

> - par_pd_type for the distribution type (as a string), and

> - par_pd_nsigmas for the limits of the distribution.

## 3.4. Using sasmodels through the bumps optimizer

Like DirectModel, you can wrap data and a kernel in a bumps model with `bumps_model.Model` and create a `bumps_model.Experiment` that you can fit with the bumps interface. Here is an example from the example directory such as example/model.py:

    import sys
    from bumps.names import *
    from sasmodels.core import load_model
    from sasmodels.bumps_model import Model, Experiment
    from sasmodels.data import load_data, set_beam_stop, set_top

    """ IMPORT THE DATA USED """
    radial_data = load_data('DEC07267.DAT')
    set_beam_stop(radial_data, 0.00669, outer=0.025)
    set_top(radial_data, -.0185)

    kernel = load_model("ellipsoid")

    model = Model(kernel,
        scale=0.08,
        radius_polar=15, radius_equatorial=800,
        sld=.291, sld_solvent=7.105,
        background=0,
        theta=90, phi=0,
        theta_pd=15, theta_pd_n=40, theta_pd_nsigma=3,
        radius_polar_pd=0.222296, radius_polar_pd_n=1, radius_polar_pd_nsigma=0,
        radius_equatorial_pd=.000128, radius_equatorial_pd_n=1, radius_equatorial_pd_nsigma=0,
        phi_pd=0, phi_pd_n=20, phi_pd_nsigma=3,
        )

    # SET THE FITTING PARAMETERS
    model.radius_polar.range(15, 1000)
    model.radius_equatorial.range(15, 1000)
    model.theta_pd.range(0, 360)
    model.background.range(0,1000)
    model.scale.range(0, 10)

    #cutoff = 0     # no cutoff on polydisperisity loops
    #cutoff = 1e-5  # default cutoff
    cutoff = 1e-3  # low precision cutoff
    M = Experiment(data=radial_data, model=model, cutoff=cutoff)
    problem = FitProblem(M)

To run the model from your python environment use the installed bumps program:

    > bumps example/model.py --preview

Alternatively, on Windows, bumps can be called from the cmd prompt as follows:

    > sasview -m bumps.cli example/model.py --preview

## 3.5. Calling the computation kernel

Getting a simple function that you can call on a set of q values and return a result is not so simple. Since the time critical use case (fitting) involves calling the function over and over with identical \(q\) values, we chose to optimize the call by only transfering the \(q\) values to the GPU once at the start of the fit. We do this by creating a `kernel.Kernel` object from the `kernel.KernelModel` returned from `core.load_model()` using the `kernel.KernelModel.make_kernel()` method. What it actual does depends on whether it is running as a DLL, as OpenCL or as a pure python kernel. Once the kernel is in hand, we can then marshal a set of parameters into a `details.CallDetails` object and ship it to the kernel using the `direct_model.call_kernel()` function. To accesses the underlying \(<F(q)>\) and \(<F^2(q)>\), use `direct_model.call_Fq()` instead.

The following example should help, example/cylinder_eval.py:

    from numpy import logspace, sqrt
    from matplotlib import pyplot as plt
    from sasmodels.core import load_model
    from sasmodels.direct_model import call_kernel, call_Fq

    model = load_model('cylinder')
    q = logspace(-3, -1, 200)
    kernel = model.make_kernel([q])
    pars = {'radius': 200, 'radius_pd': 0.1, 'scale': 2}
    Iq = call_kernel(kernel, pars)
    F, Fsq, Reff, V, Vratio = call_Fq(kernel, pars)

    plt.loglog(q, Iq, label='2 I(q)')
    plt.loglog(q, F**2/V, label='<F(q)>^2/V')
    plt.loglog(q, Fsq/V, label='<F^2(q)>/V')
    plt.xlabel('q (1/A)')
    plt.ylabel('I(q) (1/cm)')
    plt.title('Cylinder with radius 200.')
    plt.legend()
    plt.show()

Fig. 3.1 Figure 1: Comparison between \(I(q)\), \(<F(q)>\) and \(<F^2(q)>\) for cylinder model.

This compares \(I(q)\) with \(<F(q)>\) and \(<F^2(q)>\) for a cylinder with radius=200 +/- 20 and scale=2. Note that call_Fq does not include scale and background, nor does it normalize by the average volume. The definition of \(F = \rho V \hat F\) scaled by the contrast and volume, compared to the canonical cylinder \(\hat F\), with \(\hat F(0) = 1\). Integrating over polydispersity and orientation, the returned values are \(\sum_{r,w\in N(r_o, r_o/10)} \sum_\theta w F(q,r_o,\theta)\sin\theta\) and \(\sum_{r,w\in N(r_o, r_o/10)} \sum_\theta w F^2(q,r_o,\theta)\sin\theta\).

On Windows, this example can be called from the cmd prompt using sasview as as the python interpreter:

    > sasview example/cylinder_eval.py

## 3.6. Using sasmodels and bumps from a Jupyter notebook

You can also use sasmodels/bumps to fit experimental data from a Jupyter notebook by constructing and computing the model in an analogous manner to that shown above. For an example notebook see:

https://github.com/SasView/documents/blob/master/Notebooks/sasmodels_fitting.ipynb