External minimizer

In this example we are demonstrating how to run a typical fitting task in BornAgain using a third party minimizer.

The BornAgain fit parameters and minimizer interface were developed with the idea to simplify the switch between our own minimization engines and other, possibly more advanced minimization libraries. Particularly, we have been inspired by the lmfit Python package.

This makes the switch between the BornAgain and lmfit minimizers very easy.

Using the BornAgain default minimizer

import bornagain as ba

params = ba.Parameters()
params.add('radius', value=7*nm, min=5*nm, max=8*nm)
params.add('length', value=10*nm, min=8*nm, max=14*nm)

result = ba.Minimizer().minimize(fit_objective.evaluate_residuals, params)
fit_objective.finalize(result)

Using the lmfit minimizer

import lmfit

params = lmfit.Parameters()
params.add('radius', value=7*nm, min=5*nm, max=8*nm)
params.add('length', value=10*nm, min=8*nm, max=14*nm)

result = lmfit.minimize(fit_objective.evaluate_residuals, params)
fit_objective.finalize(result)

print(result.params.pretty_print())

The complete script for the lmfit based fitting is shown below.

 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
97
98
99
#!/usr/bin/env python3
"""
External minimize: using lmfit minimizers for BornAgain fits.
"""
import numpy as np
from matplotlib import pyplot as plt
import bornagain as ba
from bornagain import deg, angstrom, nm
import lmfit


def get_sample(params):
    """
    Returns a sample with cylinders and pyramids on a substrate,
    forming a hexagonal lattice.
    """
    radius = params['radius']
    lattice_length = params['length']

    m_vacuum = ba.HomogeneousMaterial("Vacuum", 0, 0)
    m_substrate = ba.HomogeneousMaterial("Substrate", 6e-6, 2e-8)
    m_particle = ba.HomogeneousMaterial("Particle", 6e-4, 2e-8)

    sphere_ff = ba.FormFactorFullSphere(radius)
    sphere = ba.Particle(m_particle, sphere_ff)
    particle_layout = ba.ParticleLayout()
    particle_layout.addParticle(sphere)

    interference = ba.InterferenceFunction2DLattice(
        ba.HexagonalLattice2D(lattice_length))
    pdf = ba.FTDecayFunction2DCauchy(10*nm, 10*nm, 0)
    interference.setDecayFunction(pdf)

    particle_layout.setInterferenceFunction(interference)

    vacuum_layer = ba.Layer(m_vacuum)
    vacuum_layer.addLayout(particle_layout)
    substrate_layer = ba.Layer(m_substrate, 0)
    multi_layer = ba.MultiLayer()
    multi_layer.addLayer(vacuum_layer)
    multi_layer.addLayer(substrate_layer)
    return multi_layer


def get_simulation(params):
    """
    Create and return GISAXS simulation with beam and detector defined
    """
    simulation = ba.GISASSimulation()
    simulation.setDetectorParameters(100, -1*deg, 1*deg, 100, 0, 2*deg)
    simulation.setBeamParameters(1*angstrom, 0.2*deg, 0)
    simulation.beam().setIntensity(1e+08)
    simulation.setSample(get_sample(params))
    return simulation


def create_real_data():
    """
    Generating "real" data by adding noise to the simulated data.
    """
    params = {'radius': 6*nm, 'length': 12*nm}
    simulation = get_simulation(params)
    simulation.runSimulation()

    # retrieving simulated data in the form of numpy array
    real_data = simulation.result().array()

    # spoiling simulated data with noise to produce "real" data
    np.random.seed(0)
    noise_factor = 0.1
    noisy = np.random.normal(real_data, noise_factor*np.sqrt(real_data))
    noisy[noisy < 0.1] = 0.1
    return noisy


def run_fitting():
    """
    main function to run fitting
    """
    real_data = create_real_data()

    fit_objective = ba.FitObjective()
    fit_objective.addSimulationAndData(get_simulation, real_data, 1)
    fit_objective.initPrint(10)

    params = lmfit.Parameters()
    params.add('radius', value=7*nm, min=5*nm, max=8*nm)
    params.add('length', value=10*nm, min=8*nm, max=14*nm)

    result = lmfit.minimize(fit_objective.evaluate_residuals, params)
    fit_objective.finalize(result)

    print(result.params.pretty_print())
    print(lmfit.fit_report(result))


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