BornAgain ≻ 24 ≻ Documentation ≻ Python scripting ≻ Fitting

Overview

In fitting, we estimate the optimum parameters in a numerical model, specifically a scattering simulation, by minimizing the difference between simulated and reference data.

BornAgain provides the simulation result. Experimental data import, parameter management, and the objective function are ordinary Python code. This keeps fitting scripts independent of any BornAgain-specific data format and lets users combine BornAgain with standard packages such as numpy, fabio, lmfit, and scipy.optimize.

In the following we will show how to fit using the BornAgain Python API.

Introductory example

In the following, a very simple example shows how to fit a parametric model to given data.

The model is a specular reflectometry scan, with a sample consisting of 20 alternating Ti and Ni layers on a Si substrate. Using this model, synthetic data have been generated with GenX. These data are part of the BornAgain sources, testdata/specular/genx_alternating_layers.dat.gz. To make them findable by the script, the environment variable BA_DATA_DIR must point to the testdata/ directory in your local BornAgain repository.

From the build directory, run as:

BA_DATA_DIR=../testdata python3 ../auto/Examples/fit/...

The fit model is identical to the model used for generating the data. There is just one fit parameter, namely the thickness of the Ti layers. The resulting fit is indistinguishable from the data:

Script

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96

#!/usr/bin/env python3
"""
Basic example how to fit specular data.
The sample consists of twenty alternating Ti and Ni layers.
Reference data was generated with GenX.
We fit just one parameter, the thickness of the Ti layers,
which has an original value of 30 angstroms.
"""

import bornagain as ba, numpy as np, os
from bornagain import angstrom, ba_fitmonitor, ba_io, nm
import lmfit


def load_data():
    datadir = os.getenv('BA_DATA_DIR', '')
    if not datadir:
        raise Exception("Environment variable BA_DATA_DIR not set")
    fname = os.path.join(datadir, "specular/genx_alternating_layers.dat.gz")

    two_alpha, intensity = ba_io.read_columns(fname, usecols=(0, 1))
    alpha = 0.5*two_alpha*ba.deg
    return alpha, intensity


def datafield_from_arrays(x, y):
    frame = ba.Frame(ba.ListScan("alpha_i (rad)", x))
    return ba.Datafield(frame, y.tolist())


def get_sample(P):
    # Materials
    vacuum = ba.Vacuum()
    ti_color = (0.05, 0.62, 0.55)
    ti_mat = ba.SLDMaterial("Ti", ti_color, -1.9493e-06, 0)
    ni_color = (0.93, 0.48, 0.14)
    ni_mat = ba.SLDMaterial("Ni", ni_color, 9.4245e-06, 0)
    si_color = (0.30, 0.62, 0.86)
    si_mat = ba.SLDMaterial("Si", si_color, 2.0704e-06, 0)

    # Layers
    layer_Ti = ba.Layer(ti_mat, P["thickness_Ti"])
    layer_Ni = ba.Layer(ni_mat, 70*angstrom)

    # Periodic stack
    n_repetitions = 10
    stack = ba.LayerStack(n_repetitions)
    stack.addLayer(layer_Ti)
    stack.addLayer(layer_Ni)

    # Sample
    sample = ba.Sample()
    sample.addLayer(ba.Layer(vacuum))
    sample.addStack(stack)
    sample.addLayer(ba.Layer(si_mat))

    return sample


def get_simulation(alpha_axis, P):
    scan = ba.AlphaScan(alpha_axis)
    scan.setWavelength(1.54*angstrom)
    sample = get_sample(P)

    return ba.SpecularSimulation(scan, sample)


if __name__ == '__main__':
    alpha, exp_values = load_data()
    exp_data = datafield_from_arrays(alpha, exp_values)
    sim_result = None  # latest simulation, shared with the plot callback

    def residuals(P):
        """
        Runs a simulation for given parameters P, and returns residuals
        vector.
        """
        global sim_result
        sim_result = get_simulation(alpha, P.valuesdict()).simulate()
        return exp_values - sim_result.intensities()

    plot_observer = ba_fitmonitor.PlotterSpecular(pause=0.5)
    def plot_iteration(P, iteration, resid):
        if iteration % 10 == 0 and sim_result is not None:
            plot_observer.plot(exp_data, sim_result, P,
                               float(np.sum(resid*resid)))

    P = lmfit.Parameters()
    P.add("thickness_Ti", 50*angstrom, min=10*angstrom, max=60*angstrom)

    result = lmfit.minimize(residuals, P, iter_cb=plot_iteration)
    print(lmfit.fit_report(result))
    finalP = result.params.valuesdict()
    ba.showSample3D(get_sample(finalP), sample_size=120*nm, seed=0)

    plot_observer.show()
auto/Examples/fit/specular/Specular1Par.py

Explanations

The arrays alpha and exp_values contain the data to be fitted. They are read with ba_io.read_columns and converted explicitly to the units required by the simulation; see Data import.

The lmfit.Parameters object P contains the fit parameters. The function residuals(P) is the objective passed to the minimizer: it converts the lmfit parameters to a plain dictionary, runs the BornAgain simulation, extracts simulated intensities as a NumPy array, and returns the residual vector.

The actual fit is controlled by the external minimizer:

result = lmfit.minimize(residuals, P)

For masked detector pixels, BornAgain simulations return NaN. Residual functions that use detector masks should replace those entries explicitly, for example with np.where(np.isfinite(r), r, 0.0). The helper ba_fit.valid_pixel_residual(exp, sim) does this, and additionally excludes non-finite or negative experimental values, the dead-pixel markers of real detectors.