Scripting Interface

Need some basic details here of how to load models and data via script, evaluate them at given parameter values and run bumps fits.

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

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

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.

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 an 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)

Assume that bumps has been installed and the bumps command is available. Maybe need to set the path to sasmodels/sasview using PYTHONPATH=path/to/sasmodels:path/to/sasview/src. To run the model use the bumps program:

$ bumps example/model.py --preview

Note that bumps and sasmodels are included as part of the SasView distribution. On windows, bumps can be called from the cmd prompt as follows:

SasViewCom bumps.cli example/model.py --preview

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()
../../../../_images/direct_call.png

Fig. 138 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:

SasViewCom example/cylinder_eval.py