BornAgain ≻ 24 ≻ Documentation ≻ Reference ≻ Fitting ≻ Examples ≻ Custom objective function

Custom objective function

BornAgain fitting scripts define the objective function directly in Python.

In this example we are going to construct a vector of residuals calculated between the experimental and simulated intensity values after applying an additional $sqrt$ function to the amplitudes.

$$ residuals = [r_{0}, r_{1}, … , r_{n-1}], ~~~ r_{i} = \sqrt{e_{i}} - \sqrt{s_{i}} $$

The length of vector n corresponds to the flattened simulation result. If a simulation uses a detector bitmap mask, masked pixels appear as NaN; custom residual functions should ignore or replace them explicitly.

This is done by defining the residual function residuals. It runs the simulation, extracts simulate().intensities(), applies sqrt to both arrays, and returns the flattened difference. The function is then passed directly to lmfit.minimize.

 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

#!/usr/bin/env python3
"""
Using custom objective function to fit GISAS data.

In this example objective function returns vector of residuals computed from
the data and simulation after applying sqrt() to intensity values.
"""
import bornagain as ba
from bornagain import deg, nm
import numpy as np
import lmfit


def get_sample(P):
    """
    Spheres on a hexagonal lattice, parameterized for fitting.
    """
    substrate_mat = ba.RefractiveMaterial(
        "Substrate", (0.28, 0.57, 0.82), 6e-6, 2e-8)
    particle_mat = ba.RefractiveMaterial(
        "Particle", (0.86, 0.24, 0.18), 6e-4, 2e-8)
    particle = ba.Particle(particle_mat, ba.Sphere(P["radius"]))

    layout = ba.Crystal2D(
        particle, ba.HexagonalLattice2D(P["length"], 0))
    layout.setDecayFunction(ba.Profile2DCauchy(100*nm, 100*nm, 0))

    particle_layer = ba.Layer(ba.Vacuum())
    particle_layer.deposit2D(layout)

    sample = ba.Sample()
    sample.addLayer(particle_layer)
    sample.addLayer(ba.Layer(substrate_mat))
    return sample


def get_simulation(P):
    """
    GISAS simulation for the parameterized hexagonal lattice.
    """
    n = 100
    beam = ba.Beam(1e8, 0.1*nm, 0.2*deg)
    detector = ba.SphericalDetector(n, -1*deg, 1*deg, n, 0, 2*deg)
    simulation = ba.ScatteringSimulation(beam, get_sample(P), detector)
    simulation.options().setUseAvgMaterials(False)
    return simulation


def fake_data():
    """
    Noisy synthetic data for a known hexagonal lattice.
    """
    P = {"radius": 6*nm, "length": 12*nm}
    return get_simulation(P).simulate().noisy(0.1, 0.1)


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
        of sqrt-transformed intensities.
        """
        simulation = get_simulation(P.valuesdict())
        flat_sim_values = simulation.simulate().intensities().ravel()
        return np.sqrt(flat_sim_values) - np.sqrt(flat_exp_values)

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

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