Fitting GISAS in 2D

This example demonstrates fitting of 2D GISAS data.

Fit script

The fit script shown and discussed below displays its progress visually, and terminates with the following result:

Fit window

 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
#!/usr/bin/env python3
"""
Minimal working fit examples: finds radius of sphere in Born approximation.
"""
import bornagain as ba
from bornagain import deg, nm, nm2
import lmfit


def get_sample(P):
    """
    Returns sample for given set of parameters.
    """
    radius = P["radius"]

    particle_color = (0.86, 0.24, 0.18)
    particle_mat = ba.RefractiveMaterial("Particle", particle_color, 6e-4, 2e-8)
    particle = ba.Particle(particle_mat, ba.Sphere(radius))

    layer_top = ba.Layer(ba.Vacuum())
    layer_bottom = ba.Layer(ba.Vacuum())
    layer_top.deposit2D(ba.Dilute2D(0.01/nm2, particle))

    sample = ba.Sample()
    sample.addLayer(layer_top)
    sample.addLayer(layer_bottom)

    return sample


def get_simulation(P):
    """
    Returns GISAS simulation for given set of parameters.
    """
    sample = get_sample(P)

    n = 100
    beam = ba.Beam(1, 0.1*nm, 0.2*deg)
    detector = ba.SphericalDetector(n, -1*deg, 1*deg, n, 0., 2*deg)
    simulation = ba.ScatteringSimulation(beam, sample, detector)

    return simulation


def fake_data():
    """
    Generating "experimental" data by running simulation with default parameters.
    """
    simulation = get_simulation({'radius': 5 * nm})
    result = simulation.simulate()
    return result


if __name__ == '__main__':
    data = fake_data()
    flat_exp_values = data.intensities().ravel()

    def residuals(P):
        """
        Runs a simulation for given parameters P, and returns residuals
        vector.
        """
        sim_values = get_simulation(P.valuesdict()).simulate().intensities()
        flat_sim_values = sim_values.ravel()
        return flat_exp_values - flat_sim_values

    P = lmfit.Parameters()
    P.add("radius", 4. * nm, min=0.01)

    result = lmfit.minimize(residuals, P)
    print(lmfit.fit_report(result))
    finalP = result.params.valuesdict()
    ba.showSample3D(get_sample(finalP), sample_size=120*nm, seed=0)
auto/Examples/fit/scatter2d/fit2d.py

Explanation

The function get_sample is basically the same as in our basic GISAS simulation example, except that radius and height of the cylindrical disks are now supplied as external parameters. These parameters are passed through a Python dictionary.

The function get_simulation has the same function argument params, from which it takes beam intensity and background level. These two parameters are on a logarithmic scale.

The residual function is defined directly in the script. It converts lmfit.Parameters to a dictionary, runs the simulation, extracts simulated intensities as a NumPy array, and returns the flattened difference to the experimental data. This residual function is passed to lmfit.minimize.