Co-refinement of multiple contrast datasets

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

[1]:
%matplotlib inline
from __future__ import print_function, division

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.0.16.dev0+e7fc2a0
scipy: 1.1.0
numpy: 1.14.5

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 - 112202351224
Dataset = e361r
datapoints = 99
chi2 = 544.5574581564464
Weighted = True
Transform = <refnx.analysis.objective.Transform object at 0x1a1fc80160>
________________________________________________________________________________
Parameters:       ''
________________________________________________________________________________
Parameters: 'instrument parameters'
<Parameter:    'scale'    value=    1.01649     +/- 0.00239, bounds=[0.9, 1.1]>
<Parameter:     'bkg'     value=  1.51192e-05   +/- 3.05e-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=[-inf, inf]>
________________________________________________________________________________
Parameters:     'SiO2'
<Parameter:'SiO2 - thick' value=    8.22988     +/- 0.277, 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=[-inf, inf]>
________________________________________________________________________________
Parameters:    'polymer'
<Parameter:'polymer - thick'value=    211.078     +/- 0.211, bounds=[200, 300]>
<Parameter:'polymer - sld'value=   0.250823     +/- 0.00632, 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.0381575    +/- 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=[-inf, inf]>


________________________________________________________________________________
Objective - 112202350888
Dataset = e365r
datapoints = 99
chi2 = 408.3378781046983
Weighted = True
Transform = <refnx.analysis.objective.Transform object at 0x1a1fc80198>
________________________________________________________________________________
Parameters:       ''
________________________________________________________________________________
Parameters: 'instrument parameters'
<Parameter:    'scale'    value=   0.955047     +/- 0.00384, bounds=[0.9, 1.1]>
<Parameter:     'bkg'     value=  2.12472e-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=[-inf, inf]>
________________________________________________________________________________
Parameters:     'SiO2'
<Parameter:'SiO2 - thick' value=    8.22988     +/- 0.277, 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=[-inf, inf]>
________________________________________________________________________________
Parameters:    'polymer'
<Parameter:'polymer - thick'value=    211.078     +/- 0.211, bounds=[200, 300]>
<Parameter:'polymer - sld'value=   0.250823     +/- 0.00632, 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.0381575    +/- 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=[-inf, inf]>


________________________________________________________________________________
Objective - 112202351280
Dataset = e366r
datapoints = 99
chi2 = 378.02526887284046
Weighted = True
Transform = <refnx.analysis.objective.Transform object at 0x1a1fc802e8>
________________________________________________________________________________
Parameters:       ''
________________________________________________________________________________
Parameters: 'instrument parameters'
<Parameter:    'scale'    value=       1        (fixed)  , bounds=[-inf, inf]>
<Parameter:     'bkg'     value=  2.16416e-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=[-inf, inf]>
________________________________________________________________________________
Parameters:     'SiO2'
<Parameter:'SiO2 - thick' value=    8.22988     +/- 0.277, 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=[-inf, inf]>
________________________________________________________________________________
Parameters:    'polymer'
<Parameter:'polymer - thick'value=    211.078     +/- 0.211, bounds=[200, 300]>
<Parameter:'polymer - sld'value=   0.250823     +/- 0.00632, 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.0381575    +/- 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=[-inf, inf]>


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 [01:34<00:00,  4.40it/s]
100%|██████████| 3000/3000 [11:47<00:00,  4.43it/s]
[15]:
global_objective.corner();
_images/reflectometry_global_26_0.png