BornAgain ≻ 1.19 ≻ Documentation ≻ Python scripting ≻ Fitting ≻ GISAS 2D fits

Fitting GISAS in 2D

This example demonstrates fitting of 2D GISAS data.

Faked experimenta data

The sample and instrument model is basically the same as in our basic GISAS simulation example, except that it depends on four external parameters: beam intensity, detector background, radius and height of the cylindrical nanoparticles.

Faked experimental data have been generated by the script devtools/fakeData/fake-gisas1.py.

The intensity for each detector pixel has been drawn at random from a Poisson distribution, with rate basically given by the fit model except for some imperfections that have been introduced on purpose to make the subsequent fitting more difficult and more realistic:

  • Refractive indices changed, and absorption strongly enhanced;
  • Cylinders replaced by a mixture of cones and segmented spheroids;
  • Beam wavelength and alpha_i slightly changed;
  • Detector x-axis skewed.

The faked data can be found at Examples/fit/data/faked-gisas1.txt.gz.

Fit script

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

Fit window

Parametric sample and simulation model

 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

"""
Parametric model of a GISAS simulation.
The idealized sample model consists of dilute cylinders on a substrate.
"""

import bornagain as ba
from bornagain import deg, nm

mat_vacuum = ba.HomogeneousMaterial("Vacuum", 0, 0)
mat_substrate = ba.HomogeneousMaterial("Substrate", 6e-6, 2e-8)
mat_particle = ba.HomogeneousMaterial("Particle", 6e-4, 2e-8)


def get_sample(params):
    cylinder_height = params["cylinder_height"]
    cylinder_radius = params["cylinder_radius"]

    ff = ba.FormFactorCylinder(cylinder_radius, cylinder_height)
    cylinder = ba.Particle(mat_particle, ff)
    layout = ba.ParticleLayout()
    layout.addParticle(cylinder)

    layer_1 = ba.Layer(mat_vacuum)
    layer_1.addLayout(layout)
    layer_2 = ba.Layer(mat_substrate)

    sample = ba.MultiLayer()
    sample.addLayer(layer_1)
    sample.addLayer(layer_2)
    return sample


def get_simulation(params):
    beam = ba.Beam(10**params['lg(intensity)'], 0.1*nm,
                   ba.Direction(0.2*deg, 0))
    det = ba.SphericalDetector(100, -1.5*deg, 1.5*deg, 100, 0, 3*deg)
    sample = get_sample(params)

    simulation = ba.GISASSimulation(beam, sample, det)
    simulation.setBackground(
        ba.ConstantBackground(10**params['lg(background)']))

    return simulation

def start_parameters_1():
    params = ba.Parameters()
    params.add("lg(intensity)", 5)
    params.add("lg(background)", 1)
    params.add("cylinder_height", 6.*nm, min=0.01)
    params.add("cylinder_radius", 6.*nm, min=0.01)
    return params
gisas_model1.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 the function argument params which can be either a Python dict, or an instance of the BornAgain class Parameters .

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 function start_parameters_1 is supplied along with the model in order to facilitate the benchmarking of different fith methods.

Main fit 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

#!/usr/bin/env python3
"""
Basic GISAS 2D fitting example.
The model is in gisas_model1.
Fake experimental data are generated by gisas_fake1.
"""

import gisas_model1 as model
import bornagain as ba
import numpy as np
from matplotlib import pyplot as plt


def run_fitting():
    real_data = np.loadtxt("../data/faked-gisas1.txt.gz", dtype=float)

    fit_objective = ba.FitObjective()
    fit_objective.addSimulationAndData(model.get_simulation, real_data)

    fit_objective.initPrint(10)  # Print on every 10th iteration.
    fit_objective.initPlot(10)  # Plot on every 10th iteration. Slow!

    minimizer = ba.Minimizer()
    params = model.start_parameters_1()
    result = minimizer.minimize(fit_objective.evaluate, params)

    fit_objective.finalize(result)

    print("Fit completed.")
    print("chi2:", result.minValue())
    for fitPar in result.parameters():
        print(fitPar.name(), fitPar.value, fitPar.error)

    # Save simulation image corresponding to the best fit parameters
    np.savetxt("fit.txt", fit_objective.simulationResult().array())


if __name__ == '__main__':
    run_fitting()
    plt.show()
fit_gisas.py

Explanation

This script loads the (faked) experimental data, and performs a fit.

The FitObjective created here is used to put into correspondence the real data, represented by a numpy array, and the simulation, represented by the get_simulation callable. On every fit iteration FitObjective

  • will generate a new simulation object for a given set of fit parameter values using the get_simulation callable,
  • run it to obtain simulated detector intensities,
  • calculate $\chi^2$ between simulated and real data.

The method fit_objective.evaluate provided by the FitObjective class interface is used here as an objective function. The method accepts fit parameters and returns the $\chi^2$ value, calculated between experimental and simulated images for the given values of the fit parameters.

The method is passed to the minimizer together with the initial fit parameter values. The minimizer.minimize starts a fit that will continue further without user intervention until the minimum is found or the minimizer failed to converge. The rest of the code demonstrates how to access the fit results.