Honeycomb_fit.py 14.5 KB
Newer Older
1
#!/usr/bin/env python3
2
3
4
5
"""
This example demonstrates how to fit a complex experimental setup using BornAgain.
It is based on real data published in  https://doi.org/10.1002/advs.201700856
by A. Glavic et al.
Wuttke, Joachim's avatar
Wuttke, Joachim committed
6
In this example we utilize the scalar reflectometry  engine to fit polarized
7
8
9
data without spin-flip for performance reasons.
"""

10
11
import os, sys
import numpy
12
13
import matplotlib.pyplot as plt
from scipy.optimize import differential_evolution
14
import bornagain as ba
15
from bornagain import angstrom, sample_tools as st
16
17

# number of points on which the computed result is plotted
Beerwerth, Randolf's avatar
Beerwerth, Randolf committed
18
scan_size = 1500
19
20
21
22
23

# restrict the Q-range of the data used for fitting
qmin = 0.08
qmax = 1.4

Wuttke, Joachim's avatar
Wuttke, Joachim committed
24
datadir = os.getenv('BA_EXAMPLE_DATA_DIR', '')
25

26
27
28
####################################################################
#                  Create Sample and Simulation                    #
####################################################################
Beerwerth, Randolf's avatar
Beerwerth, Randolf committed
29
30


31
def get_sample(parameters, sign, ms150=1):
Beerwerth, Randolf's avatar
Beerwerth, Randolf committed
32
33
34

    m_Air = ba.MaterialBySLD("Air", 0, 0)
    m_PyOx  = ba.MaterialBySLD("PyOx",
35
                               (parameters["sld_PyOx_real"] + \
Beerwerth, Randolf's avatar
Beerwerth, Randolf committed
36
                                 sign * ms150 * parameters["msld_PyOx"] )* 1e-6,
37
                               parameters["sld_PyOx_imag"] * 1e-6)
Beerwerth, Randolf's avatar
Beerwerth, Randolf committed
38
    m_Py2   = ba.MaterialBySLD("Py2",
39
40
41
                               ( parameters["sld_Py2_real"] + \
                                 sign * ms150 * parameters["msld_Py2"] ) * 1e-6,
                               parameters["sld_Py2_imag"] * 1e-6)
Beerwerth, Randolf's avatar
Beerwerth, Randolf committed
42
    m_Py1   = ba.MaterialBySLD("Py1",
43
44
45
                               ( parameters["sld_Py1_real"] + \
                                 sign * ms150 * parameters["msld_Py1"] ) * 1e-6,
                               parameters["sld_Py1_imag"] * 1e-6)
Beerwerth, Randolf's avatar
Beerwerth, Randolf committed
46
47
48
49
50
51
52
53
54
55
56
57
    m_SiO2 = ba.MaterialBySLD("SiO2", parameters["sld_SiO2_real"]*1e-6,
                              parameters["sld_SiO2_imag"]*1e-6)
    m_Si = ba.MaterialBySLD("Substrate", parameters["sld_Si_real"]*1e-6,
                            parameters["sld_Si_imag"]*1e-6)

    l_Air = ba.Layer(m_Air)
    l_PyOx = ba.Layer(m_PyOx, parameters["t_PyOx"]*angstrom)
    l_Py2 = ba.Layer(m_Py2, parameters["t_Py2"]*angstrom)
    l_Py1 = ba.Layer(m_Py1, parameters["t_Py1"]*angstrom)
    l_SiO2 = ba.Layer(m_SiO2, parameters["t_SiO2"]*angstrom)
    l_Si = ba.Layer(m_Si)

58
59
60
61
62
    r_PyOx = ba.LayerRoughness(parameters["r_PyOx"]*angstrom)
    r_Py2 = ba.LayerRoughness(parameters["r_Py2"]*angstrom)
    r_Py1 = ba.LayerRoughness(parameters["r_Py1"]*angstrom)
    r_SiO2 = ba.LayerRoughness(parameters["r_SiO2"]*angstrom)
    r_Si = ba.LayerRoughness(parameters["r_Si"]*angstrom)
Beerwerth, Randolf's avatar
Beerwerth, Randolf committed
63

64
    sample = ba.MultiLayer()
Beerwerth, Randolf's avatar
Beerwerth, Randolf committed
65

66
67
68
69
70
71
    sample.addLayer(l_Air)
    sample.addLayerWithTopRoughness(l_PyOx, r_PyOx)
    sample.addLayerWithTopRoughness(l_Py2, r_Py2)
    sample.addLayerWithTopRoughness(l_Py1, r_Py1)
    sample.addLayerWithTopRoughness(l_SiO2, r_SiO2)
    sample.addLayerWithTopRoughness(l_Si, r_Si)
Beerwerth, Randolf's avatar
Beerwerth, Randolf committed
72
73
74

    sample.setRoughnessModel(ba.RoughnessModel.NEVOT_CROCE)

75
    return sample
Beerwerth, Randolf's avatar
Beerwerth, Randolf committed
76
77


78
def get_simulation(q_axis, parameters, sign, ms150=False):
Beerwerth, Randolf's avatar
Beerwerth, Randolf committed
79

80
    q_distr = ba.DistributionGaussian(0., 1., 25, 3.)
Beerwerth, Randolf's avatar
Beerwerth, Randolf committed
81
82

    dq = parameters["dq"]*q_axis
83
    scan = ba.QzScan(q_axis)
Wuttke, Joachim's avatar
Wuttke, Joachim committed
84
    scan.setVectorResolution(q_distr, dq)
Beerwerth, Randolf's avatar
Beerwerth, Randolf committed
85
86


87
    if ms150:
Beerwerth, Randolf's avatar
Beerwerth, Randolf committed
88
89
90
        sample = get_sample(parameters=parameters,
                            sign=sign,
                            ms150=parameters["ms150"])
91
    else:
Beerwerth, Randolf's avatar
Beerwerth, Randolf committed
92
93
        sample = get_sample(parameters=parameters, sign=sign, ms150=1)

94
    simulation = ba.SpecularSimulation(scan, sample)
95
    simulation.setBackground(ba.ConstantBackground(5e-7))
Beerwerth, Randolf's avatar
Beerwerth, Randolf committed
96

97
    return simulation
Beerwerth, Randolf's avatar
Beerwerth, Randolf committed
98
99


100
def run_simulation(q_axis, fitParams, *, sign, ms150=False):
101
102
103
104
105
106
    parameters = dict(fitParams, **fixedParams)

    simulation = get_simulation(q_axis, parameters, sign, ms150)
    result = simulation.simulate()
    result.data_field().scale(parameters["intensity"])
    return result
Beerwerth, Randolf's avatar
Beerwerth, Randolf committed
107
108


109
110
def qr(result):
    """
111
    Return q and reflectivity arrays from simulation result.
112
    """
113
114
    q = numpy.array(result.convertedBinCenters(ba.Coords_QSPACE))
    r = numpy.array(result.array(ba.Coords_QSPACE))
Beerwerth, Randolf's avatar
Beerwerth, Randolf committed
115

116
    return q, r
Beerwerth, Randolf's avatar
Beerwerth, Randolf committed
117
118


119
120
121
####################################################################
#                         Plot Handling                            #
####################################################################
Beerwerth, Randolf's avatar
Beerwerth, Randolf committed
122
123


124
125
def plot(qs, rs, exps, shifts, labels, filename):
    """
126
    Plot the simulated result together with the experimental data.
127
128
129
    """
    fig = plt.figure()
    ax = fig.add_subplot(111)
Beerwerth, Randolf's avatar
Beerwerth, Randolf committed
130

131
    for q, r, exp, shift, l in zip(qs, rs, exps, shifts, labels):
Beerwerth, Randolf's avatar
Beerwerth, Randolf committed
132
133
134
135
136
137
138
139
140
141

        ax.errorbar(exp[0],
                    exp[1]/shift,
                    yerr=exp[2]/shift,
                    fmt='.',
                    markersize=0.75,
                    linewidth=0.5)

        ax.plot(q, r/shift, label=l)

142
143
    ax.set_yscale('log')
    plt.legend()
Beerwerth, Randolf's avatar
Beerwerth, Randolf committed
144

145
146
    plt.xlabel("Q [nm${}^{-1}$]")
    plt.ylabel("R")
Beerwerth, Randolf's avatar
Beerwerth, Randolf committed
147

148
149
    plt.tight_layout()
    plt.savefig(filename)
Beerwerth, Randolf's avatar
Beerwerth, Randolf committed
150
151


152
def plot_sld_profile(fitParams, filename):
Beerwerth, Randolf's avatar
Beerwerth, Randolf committed
153

154
155
    plt.figure()
    parameters = dict(fitParams, **fixedParams)
Beerwerth, Randolf's avatar
Beerwerth, Randolf committed
156

157
158
    z_300_p, sld_300_p = st.materialProfile(get_sample(parameters, 1))
    z_300_m, sld_300_m = st.materialProfile(get_sample(parameters, -1))
Beerwerth, Randolf's avatar
Beerwerth, Randolf committed
159

160
    z_150_p, sld_150_p = st.materialProfile(
Beerwerth, Randolf's avatar
Beerwerth, Randolf committed
161
        get_sample(parameters, 1, ms150=parameters["ms150"]))
162
    z_150_m, sld_150_m = st.materialProfile(
Beerwerth, Randolf's avatar
Beerwerth, Randolf committed
163
164
        get_sample(parameters, -1, ms150=parameters["ms150"]))

165
    plt.figure()
166
167
    plt.plot(z_300_p, numpy.real(sld_300_p)*1e6, label=r"300K $+$")
    plt.plot(z_300_m, numpy.real(sld_300_m)*1e6, label=r"300K $-$")
Beerwerth, Randolf's avatar
Beerwerth, Randolf committed
168

169
170
    plt.plot(z_150_p, numpy.real(sld_150_p)*1e6, label=r"150K $+$")
    plt.plot(z_150_m, numpy.real(sld_150_m)*1e6, label=r"150K $-$")
Beerwerth, Randolf's avatar
Beerwerth, Randolf committed
171
172
173
174

    plt.xlabel(r"$z$ [A]")
    plt.ylabel(r"$\delta(z) \cdot 10^6$")

175
176
177
    plt.legend()
    plt.tight_layout()
    plt.savefig(filename)
178
#    plt.close()
179
180
181
182
183
184


####################################################################
#                          Data Handling                           #
####################################################################

Beerwerth, Randolf's avatar
Beerwerth, Randolf committed
185

186
187
188
189
190
191
192
def normalizeData(data):
    """
    Removes duplicate q values from the input data,
    normalizes it such that the maximum of the reflectivity is
    unity and rescales the q-axis to inverse nm
    """
    # delete repeated data
Beerwerth, Randolf's avatar
Beerwerth, Randolf committed
193
    r0 = numpy.where(data[0] - numpy.roll(data[0], 1) == 0)
194
    data = numpy.delete(data, r0, 1)
Beerwerth, Randolf's avatar
Beerwerth, Randolf committed
195
196
197
198
199
200

    data[0] = data[0]/angstrom
    norm = numpy.max(data[1])
    data[1] = data[1]/norm
    data[2] = data[2]/norm

201
    # sort by q axis
Beerwerth, Randolf's avatar
Beerwerth, Randolf committed
202
203
204
    so = numpy.argsort(data[0])
    data = data[:, so]

205
206
    return data

Beerwerth, Randolf's avatar
Beerwerth, Randolf committed
207

208
def get_Experimental_data(filename, qmin, qmax):
Beerwerth, Randolf's avatar
Beerwerth, Randolf committed
209

210
211
212
    filepath = os.path.join(datadir, filename)
    with open(filepath, 'r') as f:
        input_Data = numpy.genfromtxt(f, unpack=True, usecols=(0, 2, 3))
213
    data = normalizeData(input_Data)
Beerwerth, Randolf's avatar
Beerwerth, Randolf committed
214
215
216
217
218

    minIndex = numpy.argmin(numpy.abs(data[0] - qmin))
    maxIndex = numpy.argmin(numpy.abs(data[0] - qmax))

    return data[:, minIndex:maxIndex + 1]
219
220
221
222
223
224


####################################################################
#                          Fit Function                            #
####################################################################

Beerwerth, Randolf's avatar
Beerwerth, Randolf committed
225

226
def relative_difference(sim, exp):
Beerwerth, Randolf's avatar
Beerwerth, Randolf committed
227
228
229
    result = (exp - sim)/(exp + sim)
    return numpy.sum(result*result)/len(sim)

230
231

def create_Parameter_dictionary(parameterNames, *args):
Beerwerth, Randolf's avatar
Beerwerth, Randolf committed
232
233
234
    return {name: value for name, value in zip(parameterNames, *args)}


235
class FitObjective:
236

237
    def __init__(self, q_axis, rdata, simulationFactory, parameterNames):
Beerwerth, Randolf's avatar
Beerwerth, Randolf committed
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
        if isinstance(q_axis, list) and isinstance(rdata, list) and \
                                    isinstance(simulationFactory, list):
            self._q = q_axis
            self._r = rdata
            self._simulationFactory = simulationFactory

        elif not isinstance(q_axis, list) and not isinstance(rdata, list) \
                              and not isinstance(simulationFactory, list):
            self._q = [q_axis]
            self._r = [rdata]
            self._simulationFactory = [simulationFactory]

        else:
            raise Exception("Inconsistent parameters")

        self._parameterNames = parameterNames

255
    def __call__(self, *args):
Wuttke, Joachim's avatar
Wuttke, Joachim committed
256
257
        fitParameters = create_Parameter_dictionary(self._parameterNames,
                                                    *args)
258
        print(f"FitParamters = {fitParameters}")
Beerwerth, Randolf's avatar
Beerwerth, Randolf committed
259

260
        result_metric = 0
Beerwerth, Randolf's avatar
Beerwerth, Randolf committed
261

262
        for q, r, sim in zip(self._q, self._r, self._simulationFactory):
263
            sim_result = sim(q, fitParameters).array()
Beerwerth, Randolf's avatar
Beerwerth, Randolf committed
264
265
            result_metric += relative_difference(sim_result, r)

266
        return result_metric
Beerwerth, Randolf's avatar
Beerwerth, Randolf committed
267
268


Wuttke, Joachim's avatar
Wuttke, Joachim committed
269
270
def run_fit_differential_evolution(q_axis, rdata, simulationFactory,
                                   startParams):
Beerwerth, Randolf's avatar
Beerwerth, Randolf committed
271
272
273
274

    bounds = [(par[1], par[2]) for n, par in startParams.items()]
    parameters = [par[0] for n, par in startParams.items()]
    parameterNames = [n for n, par in startParams.items()]
275
276
    print(f"Bounds = {bounds}")

Wuttke, Joachim's avatar
Wuttke, Joachim committed
277
278
    objective = FitObjective(q_axis, rdata, simulationFactory,
                             parameterNames)
Beerwerth, Randolf's avatar
Beerwerth, Randolf committed
279

280
    chi2_initial = objective(parameters)
Beerwerth, Randolf's avatar
Beerwerth, Randolf committed
281
282
283
284
285
286
287
288
289

    result = differential_evolution(objective,
                                    bounds,
                                    maxiter=200,
                                    popsize=len(bounds)*10,
                                    mutation=(0.5, 1.5),
                                    disp=True,
                                    tol=1e-2)

Wuttke, Joachim's avatar
Wuttke, Joachim committed
290
291
    resultParameters = create_Parameter_dictionary(parameterNames,
                                                   result.x)
292
293
294
295
296
297
298
299
300
301
302
303
    chi2_final = objective(resultParameters.values())

    print(f"Initial chi2: {chi2_initial}")
    print(f"Final chi2: {chi2_final}")
    return resultParameters


####################################################################
#                          Main Function                           #
####################################################################

if __name__ == '__main__':
Beerwerth, Randolf's avatar
Beerwerth, Randolf committed
304

305
    fixedParams = {
Beerwerth, Randolf's avatar
Beerwerth, Randolf committed
306
307
308
309
310
311
312
        "sld_PyOx_imag": (0, 0, 0),
        "sld_Py2_imag": (0, 0, 0),
        "sld_Py1_imag": (0, 0, 0),
        "sld_SiO2_imag": (0, 0, 0),
        "sld_Si_imag": (0, 0, 0),
        "sld_SiO2_real": (3.47, 3, 4),
        "sld_Si_real": (2.0704, 2, 3),
Wuttke, Joachim's avatar
Wuttke, Joachim committed
313
        "dq": (0.018, 0, 0.1),
Beerwerth, Randolf's avatar
Beerwerth, Randolf committed
314
315
    }

316
    if len(sys.argv) > 1 and sys.argv[1] == "fit":
Beerwerth, Randolf's avatar
Beerwerth, Randolf committed
317

318
319
        # some sensible start parameters for fitting
        startParams = {
Beerwerth, Randolf's avatar
Beerwerth, Randolf committed
320
321
322
323
324
325
            "intensity": (1.04, 0, 3),
            "t_PyOx": (77, 60, 100),
            "t_Py2": (56, 40, 70),
            "t_Py1": (64, 50, 80),
            "t_SiO2": (16, 10, 30),
            "sld_PyOx_real": (1.915, 1.6, 2.2),
Wuttke, Joachim's avatar
Wuttke, Joachim committed
326
            "sld_Py2_real": (5, 3, 6),
Beerwerth, Randolf's avatar
Beerwerth, Randolf committed
327
328
329
330
331
332
333
334
335
336
337
338
            "sld_Py1_real": (4.62, 3, 6),
            "r_PyOx": (27, 5, 35),
            "r_Py2": (12, 5, 20),
            "r_Py1": (12, 5, 20),
            "r_SiO2": (17, 2, 25),
            "r_Si": (18, 2, 25),
            "msld_PyOx": (0.25, 0, 1),
            "msld_Py2": (0.63, 0, 1),
            "msld_Py1": (0.64, 0, 1),
            "ms150": (1, 0.9, 1.1),
        }

339
        fit = True
Beerwerth, Randolf's avatar
Beerwerth, Randolf committed
340

341
342
    else:
        # result from our own fitting
Beerwerth, Randolf's avatar
Beerwerth, Randolf committed
343
344
345
346
347
348
349
350
351
352
353
        startParams = {
            'intensity': 0.9482344993285265,
            't_PyOx': 74.97056190221168,
            't_Py2': 61.75823766477464,
            't_Py1': 54.058310970786316,
            't_SiO2': 23.127048588278402,
            'sld_PyOx_real': 2.199791584033569,
            'sld_Py2_real': 4.980316982224387,
            'sld_Py1_real': 4.612135848532186,
            'r_PyOx': 31.323366207013787,
            'r_Py2': 9.083768897940645,
Wuttke, Joachim's avatar
Wuttke, Joachim committed
354
            'r_Py1': 5,
Beerwerth, Randolf's avatar
Beerwerth, Randolf committed
355
356
357
358
359
360
361
362
363
            'r_SiO2': 14.43455709065263,
            'r_Si': 14.948233893986075,
            'msld_PyOx': 0.292684104601585,
            'msld_Py2': 0.5979447434271339,
            'msld_Py1': 0.56376339230351,
            'ms150': 1.083311702077648
        }

        startParams = {d: (v, ) for d, v in startParams.items()}
364
365
        fit = False

Beerwerth, Randolf's avatar
Beerwerth, Randolf committed
366
367
368
369
    fixedParams = {d: v[0] for d, v in fixedParams.items()}
    paramsInitial = {d: v[0] for d, v in startParams.items()}

    def run_Simulation_300_p(qzs, params):
370
        return run_simulation(qzs, params, sign=1)
Beerwerth, Randolf's avatar
Beerwerth, Randolf committed
371
372

    def run_Simulation_300_m(qzs, params):
373
374
        return run_simulation(qzs, params, sign=-1)

Beerwerth, Randolf's avatar
Beerwerth, Randolf committed
375
    def run_Simulation_150_p(qzs, params):
376
        return run_simulation(qzs, params, sign=1, ms150=True)
Beerwerth, Randolf's avatar
Beerwerth, Randolf committed
377
378

    def run_Simulation_150_m(qzs, params):
379
        return run_simulation(qzs, params, sign=-1, ms150=True)
Beerwerth, Randolf's avatar
Beerwerth, Randolf committed
380

381
    qzs = numpy.linspace(qmin, qmax, scan_size)
Beerwerth, Randolf's avatar
Beerwerth, Randolf committed
382
383
    q_300_p, r_300_p = qr(run_Simulation_300_p(qzs, paramsInitial))
    q_300_m, r_300_m = qr(run_Simulation_300_m(qzs, paramsInitial))
384

Beerwerth, Randolf's avatar
Beerwerth, Randolf committed
385
386
    q_150_p, r_150_p = qr(run_Simulation_150_p(qzs, paramsInitial))
    q_150_m, r_150_m = qr(run_Simulation_150_m(qzs, paramsInitial))
387

388
389
    data_300_p = get_Experimental_data("honeycomb/300_p.dat", qmin, qmax)
    data_300_m = get_Experimental_data("honeycomb/300_m.dat", qmin, qmax)
390

391
392
    data_150_p = get_Experimental_data("honeycomb/150_p.dat", qmin, qmax)
    data_150_m = get_Experimental_data("honeycomb/150_m.dat", qmin, qmax)
393

Wuttke, Joachim's avatar
Wuttke, Joachim committed
394
    plot_sld_profile(paramsInitial,
395
                     "Honeycomb_Fit_sld_profile_initial.pdf")
Wuttke, Joachim's avatar
Wuttke, Joachim committed
396
397
    plot([q_300_p, q_300_m, q_150_p, q_150_m],
         [r_300_p, r_300_m, r_150_p, r_150_m],
Beerwerth, Randolf's avatar
Beerwerth, Randolf committed
398
         [data_300_p, data_300_m, data_150_p, data_150_m], [1, 1, 10, 10],
399
         ["300K $+$", "300K $-$", "150K $+$", "150K $-$"],
400
         "Honeycomb_Fit_reflectivity_initial.pdf")
401
402
403

    # fit and plot fit
    if fit:
Wuttke, Joachim's avatar
Wuttke, Joachim committed
404
405
406
        dataSimTuple = [[
            data_300_p[0], data_300_m[0], data_150_p[0], data_150_m[0]
        ], [data_300_p[1], data_300_m[1], data_150_p[1], data_150_m[1]],
Beerwerth, Randolf's avatar
Beerwerth, Randolf committed
407
408
409
410
411
                        [
                            run_Simulation_300_p, run_Simulation_300_m,
                            run_Simulation_150_p, run_Simulation_150_m
                        ]]

Wuttke, Joachim's avatar
Wuttke, Joachim committed
412
413
        fitResult = run_fit_differential_evolution(*dataSimTuple,
                                                   startParams)
Beerwerth, Randolf's avatar
Beerwerth, Randolf committed
414

415
416
        print("Fit Result:")
        print(fitResult)
Beerwerth, Randolf's avatar
Beerwerth, Randolf committed
417
418
419
420
421
422
423
424
425

        q_300_p, r_300_p = qr(run_Simulation_300_p(qzs, fitResult))
        q_300_m, r_300_m = qr(run_Simulation_300_m(qzs, fitResult))

        q_150_p, r_150_p = qr(run_Simulation_150_p(qzs, fitResult))
        q_150_m, r_150_m = qr(run_Simulation_150_m(qzs, fitResult))

        plot([q_300_p, q_300_m, q_150_p, q_150_m],
             [r_300_p, r_300_m, r_150_p, r_150_m],
Wuttke, Joachim's avatar
Wuttke, Joachim committed
426
427
             [data_300_p, data_300_m, data_150_p, data_150_m],
             [1, 1, 10, 10],
428
             ["300K $+$", "300K $-$", "150K $+$", "150K $-$"],
429
             "Honeycomb_Fit_reflectivity_fit.pdf")
Beerwerth, Randolf's avatar
Beerwerth, Randolf committed
430

431
        plot_sld_profile(fitResult, "Honeycomb_Fit_sld_profile_fit.pdf")
432
433

    plt.show()