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
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
|
#!/usr/bin/env python3
"""
An example of using the Bayesian sampling library emcee with BornAgain.
Author: Andrew McCluskey (andrew.mccluskey@ess.eu)
"""
import bornagain as ba, corner, emcee, numpy as np, os, scipy
import matplotlib.pyplot as plt
from bornagain import nm
from bornagain.numpyutil import Arrayf64Converter as dac
np.random.seed(1)
datadir = os.getenv('BA_DATA_DIR', '')
if not datadir:
raise Exception("Environment variable BA_DATA_DIR not set")
def get_sample(ni_thickness, ti_thickness):
# pure real scattering-length densities (in angstrom^-2)
si_sld_real = 2.0704e-06 # Si (substrate)
ni_sld_real = 9.4245e-06 # Ni
ti_sld_real = -1.9493e-06 # Ti
# materials
vacuum = ba.MaterialBySLD()
material_ni = ba.MaterialBySLD("Ni", ni_sld_real, 0)
material_ti = ba.MaterialBySLD("Ti", ti_sld_real, 0)
material_substrate = ba.MaterialBySLD("SiSubstrate", si_sld_real, 0)
# layers
vacuum_layer = ba.Layer(vacuum)
ni_layer = ba.Layer(material_ni, ni_thickness)
ti_layer = ba.Layer(material_ti, ti_thickness)
substrate_layer = ba.Layer(material_substrate)
# periodic stack
n_repetitions = 10
stack = ba.LayerStack(n_repetitions)
stack.addLayer(ti_layer)
stack.addLayer(ni_layer)
# sample
sample = ba.Sample()
sample.addLayer(vacuum_layer)
sample.addStack(stack)
sample.addLayer(substrate_layer)
return sample
def get_simulation(sample, points):
scan = ba.AlphaScan(ba.ListScan("alpha_i (rad)", points))
scan.setWavelength(0.154 * nm);
return ba.SpecularSimulation(scan, sample)
def run_simulation(points, ni_thickness, ti_thickness):
sample = get_sample(ni_thickness, ti_thickness)
simulation = get_simulation(sample, points)
result = simulation.simulate()
return dac.npArray(result.dataArray())
if __name__ == '__main__':
filepath = os.path.join(datadir, "specular/genx_alternating_layers.dat.gz")
flags = ba.ImportSettings1D("2alpha (deg)", "#", "", 1, 2)
data = ba.readData1D(filepath, ba.csv1D, flags)
q = dac.npArray(data.xCenters())
y = dac.asNpArray(data.dataArray())
dy = y * 0.1 # arbitrary uncertainties
def log_likelihood(P):
"""
Calculates the log-likelihood for the normal uncertainties
:tuple sim_var: the variable parameters
:array x: the abscissa data (q-values)
:array y: the ordinate data (R-values)
:array yerr: the ordinate uncertainty (dR-values)
:return: log-likelihood
"""
y_sim = run_simulation(q, *P)
sigma2 = dy**2 + y_sim**2
return -0.5*np.sum((y - y_sim)**2/sigma2 + np.log(sigma2))
# Use differential evolution to find the maximum likelihood estimate
objective_function = lambda *args: -log_likelihood(*args) # just change sign
initial = np.array([9.0, 1.0]) + 0.1*np.random.randn(2)
solution = scipy.optimize.differential_evolution(
objective_function,
bounds=((5.0, 9.0), (1.0, 10.0)))
ni_thickness_ml, ti_thickness_ml = solution.x
print('MLE Ni Thickness', ni_thickness_ml, 'nm')
print('MLE Ti Thickness', ti_thickness_ml, 'nm')
# Perform the likelihood sampling
nwalkers = 32
ndim = 2
pos = solution.x + 1e-4*np.random.randn(nwalkers, ndim)
sampler = emcee.EnsembleSampler(nwalkers,
ndim,
log_likelihood)
sampler.run_mcmc(pos, <%= test_mode ? 1 : 1000 %>, progress=True)
# Plot and show corner plot of samples
flat_samples = sampler.get_chain(flat=True)
corner.corner(flat_samples,
labels=['Ni-thickness/nm', 'Ti-thickness/nm'])
<%- if test_mode -%>
<%- elsif figure_mode -%>
<%- else -%>
plt.show()
<%- end -%>
# Plot and show MLE and data of reflectivity
plt.errorbar(q, y, dy, marker='.', ls='')
plt.plot(
q,
run_simulation(q, *flat_samples.mean(axis=0)),
'-')
plt.xlabel('$\\alpha$/rad')
plt.ylabel('$R$')
plt.yscale('log')
<%- if test_mode -%>
<%- elsif figure_mode -%>
plt.close()
<%- else -%>
plt.show()
<%- end -%>
|
