Changes in / [fbaef04:ba51e00] in sasmodels


Ignore:
Files:
11 edited

Legend:

Unmodified
Added
Removed

  • doc/guide/scripting.rst

    r4aa5dce rbd7630d  
    1010The key functions are :func:`sasmodels.core.load_model` for loading the 
    1111model definition and compiling the kernel and 
    12 :func:`sasmodels.data.load_data` for calling sasview to load the data. Need 
    13 the data because that defines the resolution function and the q values to 
    14 evaluate. If there is no data, then use :func:`sasmodels.data.empty_data1D` 
    15 or :func:`sasmodels.data.empty_data2D` to create some data with a given $q$. 
    16  
    17 Using sasmodels through bumps 
    18 ============================= 
    19  
    20 With the data and the model, you can wrap it in a *bumps* model with 
     12:func:`sasmodels.data.load_data` for calling sasview to load the data. 
     13 
     14Preparing data 
     15============== 
     16 
     17Usually you will load data via the sasview loader, with the 
     18:func:`sasmodels.data.load_data` function.  For example:: 
     19 
     20    from sasmodels.data import load_data 
     21    data = load_data("sasmodels/example/093191_201.dat") 
     22 
     23You may want to apply a data mask, such a beam stop, and trim high $q$:: 
     24 
     25    from sasmodels.data import set_beam_stop 
     26    set_beam_stop(data, qmin, qmax) 
     27 
     28The :func:`sasmodels.data.set_beam_stop` method simply sets the *mask* 
     29attribute for the data. 
     30 
     31The data defines the resolution function and the q values to evaluate, so 
     32even if you simulating experiments prior to making measurements, you still 
     33need a data object for reference. Use :func:`sasmodels.data.empty_data1D` 
     34or :func:`sasmodels.data.empty_data2D` to create a container with a 
     35given $q$ and $\Delta q/q$.  For example:: 
     36 
     37    import numpy as np 
     38    from sasmodels.data import empty_data1D 
     39 
     40    # 120 points logarithmically spaced from 0.005 to 0.2, with dq/q = 5% 
     41    q = np.logspace(np.log10(5e-3), np.log10(2e-1), 120) 
     42    data = empty_data1D(q, resolution=0.05) 
     43 
     44To use a more realistic model of resolution, or to load data from a file 
     45format not understood by SasView, you can use :class:`sasmodels.data.Data1D` 
     46or :class:`sasmodels.data.Data2D` directly.  The 1D data uses 
     47*x*, *y*, *dx* and *dy* for $x = q$ and $y = I(q)$, and 2D data uses 
     48*x*, *y*, *z*, *dx*, *dy*, *dz* for $x, y = qx, qy$ and $z = I(qx, qy)$. 
     49[Note: internally, the Data2D object uses SasView conventions, 
     50*qx_data*, *qy_data*, *data*, *dqx_data*, *dqy_data*, and *err_data*.] 
     51 
     52For USANS data, use 1D data, but set *dxl* and *dxw* attributes to 
     53indicate slit resolution:: 
     54 
     55    data.dxl = 0.117 
     56 
     57See :func:`sasmodels.resolution.slit_resolution` for details. 
     58 
     59SESANS data is more complicated; if your SESANS format is not supported by 
     60SasView you need to define a number of attributes beyond *x*, *y*.  For 
     61example:: 
     62 
     63    SElength = np.linspace(0, 2400, 61) # [A] 
     64    data = np.ones_like(SElength) 
     65    err_data = np.ones_like(SElength)*0.03 
     66 
     67    class Source: 
     68        wavelength = 6 # [A] 
     69        wavelength_unit = "A" 
     70    class Sample: 
     71        zacceptance = 0.1 # [A^-1] 
     72        thickness = 0.2 # [cm] 
     73 
     74    class SESANSData1D: 
     75        #q_zmax = 0.23 # [A^-1] 
     76        lam = 0.2 # [nm] 
     77        x = SElength 
     78        y = data 
     79        dy = err_data 
     80        sample = Sample() 
     81    data = SESANSData1D() 
     82 
     83    x, y = ... # create or load sesans 
     84    data = smd.Data 
     85 
     86The *data* module defines various data plotters as well. 
     87 
     88Using sasmodels directly 
     89======================== 
     90 
     91Once you have a computational kernel and a data object, you can evaluate 
     92the model for various parameters using 
     93:class:`sasmodels.direct_model.DirectModel`.  The resulting object *f* 
     94will be callable as *f(par=value, ...)*, returning the $I(q)$ for the $q$ 
     95values in the data.  For example:: 
     96 
     97    import numpy as np 
     98    from sasmodels.data import empty_data1D 
     99    from sasmodels.core import load_model 
     100    from sasmodels.direct_model import DirectModel 
     101 
     102    # 120 points logarithmically spaced from 0.005 to 0.2, with dq/q = 5% 
     103    q = np.logspace(np.log10(5e-3), np.log10(2e-1), 120) 
     104    data = empty_data1D(q, resolution=0.05) 
     105    kernel = load_model("ellipsoid) 
     106    f = DirectModel(data, kernel) 
     107    Iq = f(radius_polar=100) 
     108 
     109Polydispersity information is set with special parameter names: 
     110 
     111    * *par_pd* for polydispersity width, $\Delta p/p$, 
     112    * *par_pd_n* for the number of points in the distribution, 
     113    * *par_pd_type* for the distribution type (as a string), and 
     114    * *par_pd_nsigmas* for the limits of the distribution. 
     115 
     116Using sasmodels through the bumps optimizer 
     117=========================================== 
     118 
     119Like DirectModel, you can wrap data and a kernel in a *bumps* model with 
    21120class:`sasmodels.bumps_model.Model` and create an 
    22 class:`sasmodels.bump_model.Experiment` that you can fit with the *bumps* 
     121class:`sasmodels.bumps_model.Experiment` that you can fit with the *bumps* 
    23122interface. Here is an example from the *example* directory such as 
    24123*example/model.py*:: 
    … …  
    75174    SasViewCom bumps.cli example/model.py --preview 
    76175 
    77 Using sasmodels directly 
    78 ======================== 
    79  
    80 Bumps has a notion of parameter boxes in which you can set and retrieve 
    81 values.  Instead of using bumps, you can create a directly callable function 
    82 with :class:`sasmodels.direct_model.DirectModel`.  The resulting object *f* 
    83 will be callable as *f(par=value, ...)*, returning the $I(q)$ for the $q$ 
    84 values in the data.  Polydisperse parameters use the same naming conventions 
    85 as in the bumps model, with e.g., radius_pd being the polydispersity associated 
    86 with radius. 
     176Calling the computation kernel 
     177============================== 
    87178 
    88179Getting a simple function that you can call on a set of q values and return 

  • sasmodels/compare.py

    raf7a97c rbd7630d  
    12881288 
    12891289    if opts['datafile'] is not None: 
    1290         data = load_data(os.path.expanduser(opts['datafile'])) 
     1290        data0 = load_data(os.path.expanduser(opts['datafile'])) 
     1291        data = data0, data0 
    12911292    else: 
    12921293        # Hack around the fact that make_data doesn't take a pair of resolutions 

  • sasmodels/data.py

    rc036ddb rbd7630d  
    183183    *x* is spin echo length and *y* is polarization (P/P0). 
    184184    """ 
     185    isSesans = True 
    185186    def __init__(self, **kw): 
    186187        Data1D.__init__(self, **kw) 
    … …  
    301302        self.wavelength_unit = "A" 
    302303 
     304class Sample(object): 
     305    """ 
     306    Sample attributes. 
     307    """ 
     308    def __init__(self): 
     309        # type: () -> None 
     310        pass 
    303311 
    304312def empty_data1D(q, resolution=0.0, L=0., dL=0.): 

  • sasmodels/direct_model.py

    r1a8c11c r7b9e4dd  
    3131from . import resolution2d 
    3232from .details import make_kernel_args, dispersion_mesh 
     33from .modelinfo import DEFAULT_BACKGROUND 
    3334 
    3435# pylint: disable=unused-import 
    … …  
    349350 
    350351        # Need to pull background out of resolution for multiple scattering 
    351         background = pars.get('background', 0.) 
     352        background = pars.get('background', DEFAULT_BACKGROUND) 
    352353        pars = pars.copy() 
    353354        pars['background'] = 0. 

  • sasmodels/generate.py

    r5399809 r5399809  
    991991    docs = model_info.docs if model_info.docs is not None else "" 
    992992    docs = convert_section_titles_to_boldface(docs) 
    993     pars = make_partable(model_info.parameters.COMMON 
    994                          + model_info.parameters.kernel_parameters) 
     993    if model_info.structure_factor: 
     994        pars = model_info.parameters.kernel_parameters 
     995    else: 
     996        pars = model_info.parameters.COMMON + model_info.parameters.kernel_parameters 
     997    partable = make_partable(pars) 
    995998    subst = dict(id=model_info.id.replace('_', '-'), 
    996999                 name=model_info.name, 
    9971000                 title=model_info.title, 
    998                  parameters=pars, 
     1001                 parameters=partable, 
    9991002                 returns=Sq_units if model_info.structure_factor else Iq_units, 
    10001003                 docs=docs) 

  • sasmodels/kernel_iq.c

    rc57ee9e r70530778  
    192192    QACRotation *rotation, 
    193193    double qx, double qy, 
    194     double *qa_out, double *qc_out) 
     194    double *qab_out, double *qc_out) 
    195195{ 
     196    // Indirect calculation of qab, from qab^2 = |q|^2 - qc^2 
    196197    const double dqc = rotation->R31*qx + rotation->R32*qy; 
    197     // Indirect calculation of qab, from qab^2 = |q|^2 - qc^2 
    198     const double dqa = sqrt(-dqc*dqc + qx*qx + qy*qy); 
    199  
    200     *qa_out = dqa; 
     198    const double dqab_sq = -dqc*dqc + qx*qx + qy*qy; 
     199    //*qab_out = sqrt(fabs(dqab_sq)); 
     200    *qab_out = dqab_sq > 0.0 ? sqrt(dqab_sq) : 0.0; 
    201201    *qc_out = dqc; 
    202202} 

  • sasmodels/modelinfo.py

    r5399809 r5399809  
    4545# Note that scale and background cannot be coordinated parameters whose value 
    4646# depends on the some polydisperse parameter with the current implementation 
     47DEFAULT_BACKGROUND = 1e-3 
    4748COMMON_PARAMETERS = [ 
    4849    ("scale", "", 1, (0.0, np.inf), "", "Source intensity"), 
    49     ("background", "1/cm", 1e-3, (-np.inf, np.inf), "", "Source background"), 
     50    ("background", "1/cm", DEFAULT_BACKGROUND, (-np.inf, np.inf), "", "Source background"), 
    5051] 
    5152assert (len(COMMON_PARAMETERS) == 2 

  • sasmodels/models/ellipsoid.py

    rd277229 rd277229  
    125125import numpy as np 
    126126from numpy import inf, sin, cos, pi 
     127 
     128try: 
     129    from numpy import cbrt 
     130except ImportError: 
     131    def cbrt(x): return x ** (1.0/3.0) 
    127132 
    128133name = "ellipsoid" 

  • sasmodels/models/spinodal.py

    r71b751d r475ff58  
    33---------- 
    44 
    5 This model calculates the SAS signal of a phase separating solution 
    6 under spinodal decomposition. The scattering intensity $I(q)$ is calculated as 
     5This model calculates the SAS signal of a phase separating system  
     6undergoing spinodal decomposition. The scattering intensity $I(q)$ is calculated  
     7as  
    78 
    89.. math:: 
    910    I(q) = I_{max}\frac{(1+\gamma/2)x^2}{\gamma/2+x^{2+\gamma}}+B 
    1011 
    11 where $x=q/q_0$ and $B$ is a flat background. The characteristic structure 
    12 length scales with the correlation peak at $q_0$. The exponent $\gamma$ is 
    13 equal to $d+1$ with d the dimensionality of the off-critical concentration 
    14 mixtures. A transition to $\gamma=2d$ is seen near the percolation threshold 
    15 into the critical concentration regime. 
     12where $x=q/q_0$, $q_0$ is the peak position, $I_{max}$ is the intensity  
     13at $q_0$ (parameterised as the $scale$ parameter), and $B$ is a flat  
     14background. The spinodal wavelength is given by $2\pi/q_0$.  
     15 
     16The exponent $\gamma$ is equal to $d+1$ for off-critical concentration  
     17mixtures (smooth interfaces) and $2d$ for critical concentration mixtures  
     18(entangled interfaces), where $d$ is the dimensionality (ie, 1, 2, 3) of the  
     19system. Thus 2 <= $\gamma$ <= 6. A transition from $\gamma=d+1$ to $\gamma=2d$  
     20is expected near the percolation threshold.  
     21 
     22As this function tends to zero as $q$ tends to zero, in practice it may be  
     23necessary to combine it with another function describing the low-angle  
     24scattering, or to simply omit the low-angle scattering from the fit. 
    1625 
    1726References 
    … …  
    2231Physica A 123,497 (1984). 
    2332 
    24 Authorship and Verification 
    25 ---------------------------- 
     33Revision History 
     34---------------- 
    2635 
    27 * **Author:** Dirk Honecker **Date:** Oct 7, 2016 
     36* **Author:**  Dirk Honecker **Date:** Oct 7, 2016 
     37* **Revised:** Steve King    **Date:** Sep 7, 2018 
    2838""" 
    2939 
    … …  
    3444title = "Spinodal decomposition model" 
    3545description = """\ 
    36       I(q) = scale ((1+gamma/2)x^2)/(gamma/2+x^(2+gamma))+background 
     46      I(q) = Imax ((1+gamma/2)x^2)/(gamma/2+x^(2+gamma)) + background 
    3747 
    3848      List of default parameters: 
    39       scale = scaling 
    40       gamma = exponent 
    41       x = q/q_0 
     49       
     50      Imax = correlation peak intensity at q_0 
     51      background = incoherent background 
     52      gamma = exponent (see model documentation) 
    4253      q_0 = correlation peak position [1/A] 
    43       background = Incoherent background""" 
     54      x = q/q_0""" 
     55       
    4456category = "shape-independent" 
    4557 

  • sasmodels/product.py

    r33d7be3 r33d7be3  
    4242    pars = [] 
    4343    if p_info.have_Fq: 
    44         par = parse_parameter("structure_factor_mode", "", 0, [["P*S","P*(1+beta*(S-1))"]], "", 
    45                         "Structure factor calculation") 
     44        par = parse_parameter( 
     45                "structure_factor_mode", 
     46                "", 
     47                0, 
     48                [["P*S","P*(1+beta*(S-1))"]], 
     49                "", 
     50                "Structure factor calculation") 
    4651        pars.append(par) 
    4752    if p_info.effective_radius_type is not None: 
    48         par = parse_parameter("radius_effective_mode", "", 0, [["unconstrained"] + p_info.effective_radius_type], "", 
    49                         "Effective radius calculation") 
     53        par = parse_parameter( 
     54                "radius_effective_mode", 
     55                "", 
     56                0, 
     57                [["unconstrained"] + p_info.effective_radius_type], 
     58                "", 
     59                "Effective radius calculation") 
    5060        pars.append(par) 
    5161    return pars 
    … …  
    169179    )) 
    170180 
    171 def _intermediates_beta(F1, F2, S, scale, bg, effective_radius): 
    172     # type: (np.ndarray, np.ndarray, np.ndarray, np.ndarray, np.ndarray, float) -> OrderedDict[str, Union[np.ndarray, float]] 
     181def _intermediates_beta(F1,              # type: np.ndarray 
     182                        F2,              # type: np.ndarray 
     183                        S,               # type: np.ndarray 
     184                        scale,           # type: np.ndarray 
     185                        bg,              # type: np.ndarray 
     186                        effective_radius # type: float 
     187                        ): 
     188    # type: (...) -> OrderedDict[str, Union[np.ndarray, float]] 
    173189    """ 
    174190    Returns intermediate results for beta approximation-enabled product. The result may be an array or a float. 

  • sasmodels/sasview_model.py

    rd32de68 r84f2962  
    543543            # pylint: disable=unpacking-non-sequence 
    544544            q, phi = x 
    545             return self.calculate_Iq([q*math.cos(phi)], [q*math.sin(phi)])[0] 
    546         else: 
    547             return self.calculate_Iq([x])[0] 
     545            result, _ = self.calculate_Iq([q*math.cos(phi)], [q*math.sin(phi)]) 
     546            return result[0] 
     547        else: 
     548            result, _ = self.calculate_Iq([x]) 
     549            return result[0] 
    548550 
    549551 
    … …  
    560562        """ 
    561563        if isinstance(x, (list, tuple)): 
    562             return self.calculate_Iq([x[0]], [x[1]])[0] 
    563         else: 
    564             return self.calculate_Iq([x])[0] 
     564            result, _ = self.calculate_Iq([x[0]], [x[1]]) 
     565            return result[0] 
     566        else: 
     567            result, _ = self.calculate_Iq([x]) 
     568            return result[0] 
    565569 
    566570    def evalDistribution(self, qdist): 
    … …  
    596600            # Check whether we have a list of ndarrays [qx,qy] 
    597601            qx, qy = qdist 
    598             return self.calculate_Iq(qx, qy) 
     602            result, _ = self.calculate_Iq(qx, qy) 
     603            return result 
    599604 
    600605        elif isinstance(qdist, np.ndarray): 
    601606            # We have a simple 1D distribution of q-values 
    602             return self.calculate_Iq(qdist) 
     607            result, _ = self.calculate_Iq(qdist) 
     608            return result 
    603609 
    604610        else: 
    … …  
    640646        if composition and composition[0] == 'product': # only P*S for now 
    641647            with calculation_lock: 
    642                 self._calculate_Iq(qx) 
    643                 # for compatibility with sasview 4.3 
    644                 results = self._intermediate_results() 
    645                 return results["P(Q)"], results["S(Q)"] 
     648                _, lazy_results = self._calculate_Iq(qx) 
     649                # for compatibility with sasview 4.x 
     650                results = lazy_results() 
     651                pq_data = results.get("P(Q)") 
     652                sq_data = results.get("S(Q)") 
     653                return pq_data, sq_data 
    646654        else: 
    647655            return None 
    648656 
    649     def calculate_Iq(self, qx, qy=None): 
    650         # type: (Sequence[float], Optional[Sequence[float]]) -> np.ndarray 
     657    def calculate_Iq(self, 
     658                     qx,     # type: Sequence[float] 
     659                     qy=None # type: Optional[Sequence[float]] 
     660                     ): 
     661        # type: (...) -> Tuple[np.ndarray, Callable[[], collections.OrderedDict[str, np.ndarray]]] 
    651662        """ 
    652663        Calculate Iq for one set of q with the current parameters. 
    … …  
    656667        This should NOT be used for fitting since it copies the *q* vectors 
    657668        to the card for each evaluation. 
     669 
     670        The returned tuple contains the scattering intensity followed by a 
     671        callable which returns a dictionary of intermediate data from 
     672        ProductKernel. 
    658673        """ 
    659674        ## uncomment the following when trying to debug the uncoordinated calls 
    … …  
    688703        result = calculator(call_details, values, cutoff=self.cutoff, 
    689704                            magnetic=is_magnetic) 
     705        lazy_results = getattr(calculator, 'results', 
     706                               lambda: collections.OrderedDict()) 
    690707        #print("result", result) 
    691         self._intermediate_results = getattr(calculator, 'results', None) 
     708 
    692709        calculator.release() 
    693710        #self._model.release() 
    694         return result 
     711 
     712        return result, lazy_results 
    695713 
    696714    def calculate_ER(self): 
    … …  
    881899    q = np.linspace(-0.35, 0.35, 500) 
    882900    qx, qy = np.meshgrid(q, q) 
    883     result = model.calculate_Iq(qx.flatten(), qy.flatten()) 
     901    result, _ = model.calculate_Iq(qx.flatten(), qy.flatten()) 
    884902    result = result.reshape(qx.shape) 
    885903 
Note: See TracChangeset for help on using the changeset viewer.