Co-refinement of multiple contrast datasets

A demonstration of how to do co-refinement of several datasets with refnx.

[1]:
%matplotlib inline
import os.path

from refnx.dataset import ReflectDataset
from refnx.analysis import Transform, CurveFitter, Objective, GlobalObjective, Parameter
from refnx.reflect import SLD, ReflectModel
import refnx

import scipy
import numpy as np
import matplotlib.pyplot as plt
[2]:
print('refnx: %s\nscipy: %s\nnumpy: %s' % (refnx.version.version, scipy.version.version, np.version.version))
refnx: 0.1.12.dev0+497f93a
scipy: 1.5.0.dev0+f614064
numpy: 1.18.2

These are datasets used in refnx testing, distributed with every refnx install.

[3]:
pth = os.path.dirname(refnx.__file__)

e361 = ReflectDataset(os.path.join(pth, 'analysis', 'test', 'e361r.txt'))
e365 = ReflectDataset(os.path.join(pth, 'analysis', 'test', 'e365r.txt'))
e366 = ReflectDataset(os.path.join(pth, 'analysis', 'test', 'e366r.txt'))

Make some SLD objects to represent all the materials.

[4]:
si = SLD(2.07, 'Si')
sio2 = SLD(3.47, 'SiO2')
polymer = SLD(2.0, 'polymer')
d2o = SLD(6.36, 'D2O')
h2o = SLD(-0.56, 'H2O')
cm3 = SLD(3.5, 'cm3.5')

The SLDs are used to create Slabs.

[5]:
sio2_l = sio2(30, 3)
polymer_l = polymer(250, 3)

# we're going to share the water/polymer roughness across all 3 datasets
water_poly_rough = Parameter(3, 'water_poly_rough')
d2o_l = d2o(0, water_poly_rough)
h2o_l = h2o(0, water_poly_rough)
cm3_l = cm3(0, water_poly_rough)

Set the limits for the parameters we wish to vary. Each contrast uses the same polymer SLD. We account for contrast change using the volume fraction of solvent.

[6]:
sio2_l.thick.setp(vary=True, bounds=(1, 50))

polymer_l.thick.setp(vary=True, bounds=(200, 300))
polymer_l.sld.real.setp(vary=True, bounds=(0.1, 2))
polymer_l.vfsolv.setp(vary=True, bounds=(0, 1))

We create a different Structure for each contrast of interest. It’s important to note here that the Structures all share the same si, sio2_l, polymer_l objects. This means the Structures all share the same parameters. The only thing that’s different is the solvent contrast. By default the Structure object solvates with the SLD of the last slab. This behaviour can be modified by changing the Structure.solvent attribute.

[7]:
structure361 = si | sio2_l | polymer_l | d2o_l
structure365 = si | sio2_l | polymer_l | cm3_l
structure366 = si | sio2_l | polymer_l | h2o_l

Create a ReflectModel from the Structure. These are responsible for calculating the generative model, doing resolution smearing, applying a scale factor, and adding a Q-independent constant background.

[8]:
model361 = ReflectModel(structure361)
model365 = ReflectModel(structure365)
model366 = ReflectModel(structure366)

model361.scale.setp(vary=True, bounds=(0.9, 1.1))
model361.bkg.setp(vary=True, bounds=(0.9e-8, 3e-5))
model365.scale.setp(vary=True, bounds=(0.9, 1.1))
model365.bkg.setp(vary=True, bounds=(0.9e-8, 3e-5))
model366.bkg.setp(vary=True, bounds=(0.9e-8, 3e-5))

Objectives are created from the datasets and the model. Here we also add a Transform to fit as logR vs Q.

[9]:
objective361 = Objective(model361, e361, transform=Transform('logY'))
objective365 = Objective(model365, e365, transform=Transform('logY'))
objective366 = Objective(model366, e366, transform=Transform('logY'))

A GlobalObjective is formed from the individual Objectives. This means that they’re all analysed together.

[10]:
global_objective = GlobalObjective([objective361, objective365, objective366])

Create the CurveFitter and fit with differential evolution.

[11]:
# create the fit instance
fitter = CurveFitter(global_objective)
fitter.fit('differential_evolution');
[12]:
global_objective.plot()
plt.legend()
plt.xlabel('Q')
plt.ylabel('logR');
_images/reflectometry_global_22_0.png
[13]:
print(global_objective)
_______________________________________________________________________________

--Global Objective--
________________________________________________________________________________
Objective - 4667369904
Dataset = e361r
datapoints = 99
chi2 = 563.0497329906885
Weighted = True
Transform = Transform('logY')
________________________________________________________________________________
Parameters:       ''
________________________________________________________________________________
Parameters: 'instrument parameters'
<Parameter:    'scale'    , value=1.01633 +/- 0.00239, bounds=[0.9, 1.1]>
<Parameter:     'bkg'     , value=1.51658e-05 +/- 3.06e-07, bounds=[9e-09, 3e-05]>
<Parameter:'dq - resolution', value=5 (fixed)  , bounds=[-inf, inf]>
________________________________________________________________________________
Parameters: 'Structure - '
________________________________________________________________________________
Parameters:      'Si'
<Parameter: 'Si - thick'  , value=0 (fixed)  , bounds=[-inf, inf]>
<Parameter:  'Si - sld'   , value=2.07 (fixed)  , bounds=[-inf, inf]>
<Parameter:  'Si - isld'  , value=0 (fixed)  , bounds=[-inf, inf]>
<Parameter: 'Si - rough'  , value=0 (fixed)  , bounds=[-inf, inf]>
<Parameter:'Si - volfrac solvent', value=0 (fixed)  , bounds=[0.0, 1.0]>
________________________________________________________________________________
Parameters:     'SiO2'
<Parameter:'SiO2 - thick' , value=7.90349 +/- 0.284, bounds=[1, 50]>
<Parameter: 'SiO2 - sld'  , value=3.47 (fixed)  , bounds=[-inf, inf]>
<Parameter: 'SiO2 - isld' , value=0 (fixed)  , bounds=[-inf, inf]>
<Parameter:'SiO2 - rough' , value=3 (fixed)  , bounds=[-inf, inf]>
<Parameter:'SiO2 - volfrac solvent', value=0 (fixed)  , bounds=[0.0, 1.0]>
________________________________________________________________________________
Parameters:    'polymer'
<Parameter:'polymer - thick', value=210.987 +/- 0.218, bounds=[200, 300]>
<Parameter:'polymer - sld', value=0.248507 +/- 0.00631, bounds=[0.1, 2]>
<Parameter:'polymer - isld', value=0 (fixed)  , bounds=[-inf, inf]>
<Parameter:'polymer - rough', value=3 (fixed)  , bounds=[-inf, inf]>
<Parameter:'polymer - volfrac solvent', value=0.0372458 +/- 0.00216, bounds=[0, 1]>
________________________________________________________________________________
Parameters:      'D2O'
<Parameter: 'D2O - thick' , value=0 (fixed)  , bounds=[-inf, inf]>
<Parameter:  'D2O - sld'  , value=6.36 (fixed)  , bounds=[-inf, inf]>
<Parameter: 'D2O - isld'  , value=0 (fixed)  , bounds=[-inf, inf]>
<Parameter:'water_poly_rough', value=3 (fixed)  , bounds=[-inf, inf]>
<Parameter:'D2O - volfrac solvent', value=0 (fixed)  , bounds=[0.0, 1.0]>


________________________________________________________________________________
Objective - 4667369960
Dataset = e365r
datapoints = 99
chi2 = 406.84293694513104
Weighted = True
Transform = Transform('logY')
________________________________________________________________________________
Parameters:       ''
________________________________________________________________________________
Parameters: 'instrument parameters'
<Parameter:    'scale'    , value=0.955107 +/- 0.00384, bounds=[0.9, 1.1]>
<Parameter:     'bkg'     , value=2.13006e-05 +/- 2.47e-07, bounds=[9e-09, 3e-05]>
<Parameter:'dq - resolution', value=5 (fixed)  , bounds=[-inf, inf]>
________________________________________________________________________________
Parameters: 'Structure - '
________________________________________________________________________________
Parameters:      'Si'
<Parameter: 'Si - thick'  , value=0 (fixed)  , bounds=[-inf, inf]>
<Parameter:  'Si - sld'   , value=2.07 (fixed)  , bounds=[-inf, inf]>
<Parameter:  'Si - isld'  , value=0 (fixed)  , bounds=[-inf, inf]>
<Parameter: 'Si - rough'  , value=0 (fixed)  , bounds=[-inf, inf]>
<Parameter:'Si - volfrac solvent', value=0 (fixed)  , bounds=[0.0, 1.0]>
________________________________________________________________________________
Parameters:     'SiO2'
<Parameter:'SiO2 - thick' , value=7.90349 +/- 0.284, bounds=[1, 50]>
<Parameter: 'SiO2 - sld'  , value=3.47 (fixed)  , bounds=[-inf, inf]>
<Parameter: 'SiO2 - isld' , value=0 (fixed)  , bounds=[-inf, inf]>
<Parameter:'SiO2 - rough' , value=3 (fixed)  , bounds=[-inf, inf]>
<Parameter:'SiO2 - volfrac solvent', value=0 (fixed)  , bounds=[0.0, 1.0]>
________________________________________________________________________________
Parameters:    'polymer'
<Parameter:'polymer - thick', value=210.987 +/- 0.218, bounds=[200, 300]>
<Parameter:'polymer - sld', value=0.248507 +/- 0.00631, bounds=[0.1, 2]>
<Parameter:'polymer - isld', value=0 (fixed)  , bounds=[-inf, inf]>
<Parameter:'polymer - rough', value=3 (fixed)  , bounds=[-inf, inf]>
<Parameter:'polymer - volfrac solvent', value=0.0372458 +/- 0.00216, bounds=[0, 1]>
________________________________________________________________________________
Parameters:     'cm3.5'
<Parameter:'cm3.5 - thick', value=0 (fixed)  , bounds=[-inf, inf]>
<Parameter: 'cm3.5 - sld' , value=3.5 (fixed)  , bounds=[-inf, inf]>
<Parameter:'cm3.5 - isld' , value=0 (fixed)  , bounds=[-inf, inf]>
<Parameter:'water_poly_rough', value=3 (fixed)  , bounds=[-inf, inf]>
<Parameter:'cm3.5 - volfrac solvent', value=0 (fixed)  , bounds=[0.0, 1.0]>


________________________________________________________________________________
Objective - 4667369624
Dataset = e366r
datapoints = 99
chi2 = 361.5822051208103
Weighted = True
Transform = Transform('logY')
________________________________________________________________________________
Parameters:       ''
________________________________________________________________________________
Parameters: 'instrument parameters'
<Parameter:    'scale'    , value=1 (fixed)  , bounds=[-inf, inf]>
<Parameter:     'bkg'     , value=2.17328e-05 +/- 2.52e-07, bounds=[9e-09, 3e-05]>
<Parameter:'dq - resolution', value=5 (fixed)  , bounds=[-inf, inf]>
________________________________________________________________________________
Parameters: 'Structure - '
________________________________________________________________________________
Parameters:      'Si'
<Parameter: 'Si - thick'  , value=0 (fixed)  , bounds=[-inf, inf]>
<Parameter:  'Si - sld'   , value=2.07 (fixed)  , bounds=[-inf, inf]>
<Parameter:  'Si - isld'  , value=0 (fixed)  , bounds=[-inf, inf]>
<Parameter: 'Si - rough'  , value=0 (fixed)  , bounds=[-inf, inf]>
<Parameter:'Si - volfrac solvent', value=0 (fixed)  , bounds=[0.0, 1.0]>
________________________________________________________________________________
Parameters:     'SiO2'
<Parameter:'SiO2 - thick' , value=7.90349 +/- 0.284, bounds=[1, 50]>
<Parameter: 'SiO2 - sld'  , value=3.47 (fixed)  , bounds=[-inf, inf]>
<Parameter: 'SiO2 - isld' , value=0 (fixed)  , bounds=[-inf, inf]>
<Parameter:'SiO2 - rough' , value=3 (fixed)  , bounds=[-inf, inf]>
<Parameter:'SiO2 - volfrac solvent', value=0 (fixed)  , bounds=[0.0, 1.0]>
________________________________________________________________________________
Parameters:    'polymer'
<Parameter:'polymer - thick', value=210.987 +/- 0.218, bounds=[200, 300]>
<Parameter:'polymer - sld', value=0.248507 +/- 0.00631, bounds=[0.1, 2]>
<Parameter:'polymer - isld', value=0 (fixed)  , bounds=[-inf, inf]>
<Parameter:'polymer - rough', value=3 (fixed)  , bounds=[-inf, inf]>
<Parameter:'polymer - volfrac solvent', value=0.0372458 +/- 0.00216, bounds=[0, 1]>
________________________________________________________________________________
Parameters:      'H2O'
<Parameter: 'H2O - thick' , value=0 (fixed)  , bounds=[-inf, inf]>
<Parameter:  'H2O - sld'  , value=-0.56 (fixed)  , bounds=[-inf, inf]>
<Parameter: 'H2O - isld'  , value=0 (fixed)  , bounds=[-inf, inf]>
<Parameter:'water_poly_rough', value=3 (fixed)  , bounds=[-inf, inf]>
<Parameter:'H2O - volfrac solvent', value=0 (fixed)  , bounds=[0.0, 1.0]>


Now we’re going to do some MCMC sampling. We discard the first 400 steps, then save 1 in every 100 steps, for a total of 30 saved steps.

[14]:
fitter.sample(400, random_state=1)
fitter.sampler.reset()
fitter.sample(30, nthin=100, random_state=1);
100%|██████████| 400/400 [02:05<00:00,  3.31it/s]
100%|██████████| 3000/3000 [15:17<00:00,  3.34it/s]
[15]:
global_objective.corner();
_images/reflectometry_global_26_0.png