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
|
#!/usr/bin/env python3
"""
External minimize: using lmfit minimizers for BornAgain fits.
"""
import bornagain as ba
from bornagain import ba_fitmonitor, deg, nm
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
vector.
"""
simulation = get_simulation(P.valuesdict())
flat_sim_values = simulation.simulate().intensities().ravel()
return flat_exp_values - flat_sim_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,
iter_cb=ba_fitmonitor.Printer(10))
print(lmfit.fit_report(result))
finalP = result.params.valuesdict()
ba.showSample3D(get_sample(finalP), sample_size=300*nm, seed=0)
|