Analysis¶

It is time to analyse some data. We here show how to set up an Analysis object and use it to first fit an artificial vanadium measurement. Next, we use the fitted resolution to fit an artificial measurement of a model with diffusion and some elastic scattering.

We extract and plot the relevant parameters. Finally, we show how to fit directly to the diffusion model.

In the near future, it will be possible to fit the width and area of the Lorentzian to the diffusion model as well.

In [1]:

# Imports
import pooch

from easydynamics.analysis.analysis import Analysis
from easydynamics.experiment import Experiment
from easydynamics.sample_model import BrownianTranslationalDiffusion
from easydynamics.sample_model import ComponentCollection
from easydynamics.sample_model import DeltaFunction
from easydynamics.sample_model import Gaussian
from easydynamics.sample_model import Lorentzian
from easydynamics.sample_model import Polynomial
from easydynamics.sample_model.background_model import BackgroundModel
from easydynamics.sample_model.instrument_model import InstrumentModel
from easydynamics.sample_model.resolution_model import ResolutionModel
from easydynamics.sample_model.sample_model import SampleModel

%matplotlib widget

# Imports import pooch from easydynamics.analysis.analysis import Analysis from easydynamics.experiment import Experiment from easydynamics.sample_model import BrownianTranslationalDiffusion from easydynamics.sample_model import ComponentCollection from easydynamics.sample_model import DeltaFunction from easydynamics.sample_model import Gaussian from easydynamics.sample_model import Lorentzian from easydynamics.sample_model import Polynomial from easydynamics.sample_model.background_model import BackgroundModel from easydynamics.sample_model.instrument_model import InstrumentModel from easydynamics.sample_model.resolution_model import ResolutionModel from easydynamics.sample_model.sample_model import SampleModel %matplotlib widget

In [2]:

# Load the vanadium data
vanadium_experiment = Experiment('Vanadium')
file_path = pooch.retrieve(
    url='https://github.com/easyscience/dynamics-lib/raw/refs/heads/docs/docs/docs/tutorials/data/vanadium_data_example.h5',
    known_hash='16cc1b327c303feeb88fb9dda5390dc4880b62396b1793f98c6fef0b27c7b873',
)

vanadium_experiment.load_hdf5(filename=file_path)

# Load the vanadium data vanadium_experiment = Experiment('Vanadium') file_path = pooch.retrieve( url='https://github.com/easyscience/dynamics-lib/raw/refs/heads/docs/docs/docs/tutorials/data/vanadium_data_example.h5', known_hash='16cc1b327c303feeb88fb9dda5390dc4880b62396b1793f98c6fef0b27c7b873', ) vanadium_experiment.load_hdf5(filename=file_path)

Downloading data from 'https://github.com/easyscience/dynamics-lib/raw/refs/heads/docs/docs/docs/tutorials/data/vanadium_data_example.h5' to file '/home/runner/.cache/pooch/67f654c1ef73b5bbe61d7cbd3ffd0657-vanadium_data_example.h5'.

/home/runner/work/dynamics-lib/dynamics-lib/.pixi/envs/default/lib/python3.12/site-packages/requests/__init__.py:113: RequestsDependencyWarning: urllib3 (2.6.3) or chardet (7.0.1)/charset_normalizer (3.4.5) doesn't match a supported version!
  warnings.warn(

In [3]:

# Example of Analysis with a simple sample model and instrument model
# The scattering from vanadium is purely elastic, so we model it with a
# delta function
delta_function = DeltaFunction(display_name='DeltaFunction', area=1)
sample_model = SampleModel(
    components=delta_function,
)

# The resolution is in this case modeled as a Gaussian. However, we can
# add as many components as we like to the resolution model
res_gauss = Gaussian(width=0.1)
res_gauss.area.fixed = True
resolution_components = ComponentCollection()
resolution_components.append_component(res_gauss)
resolution_model = ResolutionModel(components=resolution_components)

# The background model is created in the same way. In this case, we use
# a flat background
background_model = BackgroundModel(components=Polynomial(coefficients=[0.001]))

# We combine the resolution abd background model into an instrument
# model. This model also contains a small energy offset to account for
# instrument misalignment.

instrument_model = InstrumentModel(
    resolution_model=resolution_model,
    background_model=background_model,
)

# Collect everything into an analysis object.
vanadium_analysis = Analysis(
    display_name='Vanadium Full Analysis',
    experiment=vanadium_experiment,
    sample_model=sample_model,
    instrument_model=instrument_model,
)

# Let us first fit a single Q index and plot the data and model to see
# how it looks
fit_result_independent_single_Q = vanadium_analysis.fit(fit_method='independent', Q_index=5)
vanadium_analysis.plot_data_and_model(Q_index=5)

# Example of Analysis with a simple sample model and instrument model # The scattering from vanadium is purely elastic, so we model it with a # delta function delta_function = DeltaFunction(display_name='DeltaFunction', area=1) sample_model = SampleModel( components=delta_function, ) # The resolution is in this case modeled as a Gaussian. However, we can # add as many components as we like to the resolution model res_gauss = Gaussian(width=0.1) res_gauss.area.fixed = True resolution_components = ComponentCollection() resolution_components.append_component(res_gauss) resolution_model = ResolutionModel(components=resolution_components) # The background model is created in the same way. In this case, we use # a flat background background_model = BackgroundModel(components=Polynomial(coefficients=[0.001])) # We combine the resolution abd background model into an instrument # model. This model also contains a small energy offset to account for # instrument misalignment. instrument_model = InstrumentModel( resolution_model=resolution_model, background_model=background_model, ) # Collect everything into an analysis object. vanadium_analysis = Analysis( display_name='Vanadium Full Analysis', experiment=vanadium_experiment, sample_model=sample_model, instrument_model=instrument_model, ) # Let us first fit a single Q index and plot the data and model to see # how it looks fit_result_independent_single_Q = vanadium_analysis.fit(fit_method='independent', Q_index=5) vanadium_analysis.plot_data_and_model(Q_index=5)

Out[3]:

In [4]:

# It looks good, so let us fit all Q indices independently and plot the
# results
fit_result_independent_all_Q = vanadium_analysis.fit(fit_method='independent')
vanadium_analysis.plot_data_and_model()

# It looks good, so let us fit all Q indices independently and plot the # results fit_result_independent_all_Q = vanadium_analysis.fit(fit_method='independent') vanadium_analysis.plot_data_and_model()

Out[4]:

In [5]:

# Inspect the Parameters as a scipp Dataset
vanadium_analysis.parameters_to_dataset()

# Inspect the Parameters as a scipp Dataset vanadium_analysis.parameters_to_dataset()

Out[5]:

scipp.Dataset (6.66 KB)

• Dimensions:
  • Q: 16
• Coordinates: (1)
  • Q
    (Q)
    float64
    1/Å
    0.1, 0.227, ..., 1.873, 2.0

    Values:
    array([0.1       , 0.22666667, 0.35333333, 0.48      , 0.60666667,
           0.73333333, 0.86      , 0.98666667, 1.11333333, 1.24      ,
           1.36666667, 1.49333333, 1.62      , 1.74666667, 1.87333333,
           2.        ])

• Data: (7)
  • DeltaFunction area
    (Q)
    float64
    meV
    0.522, 0.529, ..., 0.538, 0.533

    σ = 0.018, 0.018, ..., 0.017, 0.017

    Values:
    array([0.52209408, 0.5285229 , 0.51741732, 0.53239306, 0.52764215,
           0.5299243 , 0.52765768, 0.51477418, 0.52913757, 0.53258438,
           0.5324246 , 0.53238945, 0.52271039, 0.52427006, 0.53758902,
           0.53293296])
    
    Variances (σ²):
    array([0.00030823, 0.00033506, 0.00025861, 0.00028047, 0.00027407,
           0.0002244 , 0.00044408, 0.00024497, 0.00029737, 0.00034951,
           0.00030663, 0.00047018, 0.00037023, 0.00033697, 0.00028746,
           0.00028007])

  • DeltaFunction center
    (Q)
    float64
    meV
    0.0, 0.0, ..., 0.0, 0.0

    σ = 0.0, 0.0, ..., 0.0, 0.0

    Values:
    array([0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.])
    
    Variances (σ²):
    array([0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.])

  • Gaussian area
    (Q)
    float64
    meV
    1.0, 1.0, ..., 1.0, 1.0

    σ = 0.0, 0.0, ..., 0.0, 0.0

    Values:
    array([1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.])
    
    Variances (σ²):
    array([0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.])

  • Gaussian center
    (Q)
    float64
    meV
    0.0, 0.0, ..., 0.0, 0.0

    σ = 0.0, 0.0, ..., 0.0, 0.0

    Values:
    array([0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.])
    
    Variances (σ²):
    array([0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.])

  • Gaussian width
    (Q)
    float64
    meV
    0.102, 0.100, ..., 0.103, 0.102

    σ = 0.003, 0.003, ..., 0.003, 0.003

    Values:
    array([0.10222849, 0.09980219, 0.10094395, 0.10308798, 0.10108169,
           0.10126272, 0.10063122, 0.09902007, 0.10372129, 0.10254205,
           0.10142844, 0.10315841, 0.10076283, 0.10044881, 0.10346177,
           0.10187953])
    
    Variances (σ²):
    array([9.25121501e-06, 8.81462902e-06, 7.61444384e-06, 8.02505303e-06,
           7.63407023e-06, 6.32632931e-06, 1.23202070e-05, 6.87361362e-06,
           8.78013502e-06, 9.96635847e-06, 8.46568221e-06, 1.33366848e-05,
           1.06131452e-05, 9.28744245e-06, 8.02343194e-06, 7.89725747e-06])

  • Polynomial_c0
    (Q)
    float64
    𝟙
    0.099, 0.095, ..., 0.098, 0.101

    σ = 0.002, 0.002, ..., 0.002, 0.002

    Values:
    array([0.09942061, 0.09482336, 0.10250515, 0.09830358, 0.10212681,
           0.1023704 , 0.09957118, 0.1020655 , 0.09643967, 0.10180256,
           0.10191825, 0.09198844, 0.09427405, 0.09806923, 0.09764543,
           0.10051938])
    
    Variances (σ²):
    array([5.53666458e-06, 5.69780006e-06, 4.81032728e-06, 4.89897712e-06,
           4.99355959e-06, 4.08027374e-06, 7.91272766e-06, 4.56234563e-06,
           5.13312615e-06, 6.29261176e-06, 5.53177991e-06, 7.73519893e-06,
           6.33375780e-06, 5.95523123e-06, 4.94577573e-06, 4.98191241e-06])

  • energy_offset
    (Q)
    float64
    meV
    0.001, -0.001, ..., -0.001, 0.001

    σ = 0.004, 0.004, ..., 0.003, 0.003

    Values:
    array([ 0.00135435, -0.00052522,  0.00093028, -0.00192069, -0.00366331,
            0.00128853,  0.000512  ,  0.00090986, -0.00100558,  0.00036124,
            0.00054189, -0.00216345, -0.00137654,  0.00129853, -0.001399  ,
            0.00055349])
    
    Variances (σ²):
    array([1.33891037e-05, 1.37163920e-05, 1.12418484e-05, 1.18955355e-05,
           1.14728274e-05, 9.35171956e-06, 1.83656779e-05, 1.03814854e-05,
           1.29102630e-05, 1.47729559e-05, 1.26723273e-05, 1.99321937e-05,
           1.56095496e-05, 1.41010062e-05, 1.20573093e-05, 1.16915411e-05])

In [6]:

# Plot some of fitted parameters as a function of Q
vanadium_analysis.plot_parameters(names=['DeltaFunction area'])

# Plot some of fitted parameters as a function of Q vanadium_analysis.plot_parameters(names=['DeltaFunction area'])

Out[6]:

In [7]:

vanadium_analysis.plot_parameters(names=['Gaussian width'])

vanadium_analysis.plot_parameters(names=['Gaussian width'])

Out[7]:

In [8]:

vanadium_analysis.plot_parameters(names=['energy_offset'])

vanadium_analysis.plot_parameters(names=['energy_offset'])

Out[8]:

In [9]:

# Now it's time to look at the data we want to fit. We first load the
# data
diffusion_experiment = Experiment('Diffusion')

file_path = pooch.retrieve(
    url='https://github.com/easyscience/dynamics-lib/raw/refs/heads/docs/docs/docs/tutorials/data/diffusion_data_example.h5',
)

diffusion_experiment.load_hdf5(filename=file_path)

# Now it's time to look at the data we want to fit. We first load the # data diffusion_experiment = Experiment('Diffusion') file_path = pooch.retrieve( url='https://github.com/easyscience/dynamics-lib/raw/refs/heads/docs/docs/docs/tutorials/data/diffusion_data_example.h5', ) diffusion_experiment.load_hdf5(filename=file_path)

Downloading data from 'https://github.com/easyscience/dynamics-lib/raw/refs/heads/docs/docs/docs/tutorials/data/diffusion_data_example.h5' to file '/home/runner/.cache/pooch/23ec0ed9b23916e9513b147db8a588d9-diffusion_data_example.h5'.

SHA256 hash of downloaded file: 5fe846b19aacbda8b8b936eb2e5310d025dc56c25b0b353521e7d6b921f229ab
Use this value as the 'known_hash' argument of 'pooch.retrieve' to ensure that the file hasn't changed if it is downloaded again in the future.

In [10]:

# Now we set up the model, similarly to how we set up the model for the
# vanadium data.

delta_function = DeltaFunction(display_name='DeltaFunction', area=0.2)
lorentzian = Lorentzian(display_name='Lorentzian', area=0.5, width=0.3)
component_collection = ComponentCollection(
    components=[delta_function, lorentzian],
)

sample_model = SampleModel(
    components=component_collection,
)

background_model = BackgroundModel(components=Polynomial(coefficients=[0.001]))

instrument_model = InstrumentModel(
    background_model=background_model,
)

diffusion_analysis = Analysis(
    display_name='Diffusion Full Analysis',
    experiment=diffusion_experiment,
    sample_model=sample_model,
    instrument_model=instrument_model,
)

# We need to hack in the resolution model from the vanadium analysis,
# since the setters and getters overwrite the model. This will be fixed
# asap.
diffusion_analysis.instrument_model._resolution_model = (
    vanadium_analysis.instrument_model.resolution_model
)

# We fix all parameters of the resolution model.
diffusion_analysis.instrument_model.resolution_model.fix_all_parameters()

# Now we set up the model, similarly to how we set up the model for the # vanadium data. delta_function = DeltaFunction(display_name='DeltaFunction', area=0.2) lorentzian = Lorentzian(display_name='Lorentzian', area=0.5, width=0.3) component_collection = ComponentCollection( components=[delta_function, lorentzian], ) sample_model = SampleModel( components=component_collection, ) background_model = BackgroundModel(components=Polynomial(coefficients=[0.001])) instrument_model = InstrumentModel( background_model=background_model, ) diffusion_analysis = Analysis( display_name='Diffusion Full Analysis', experiment=diffusion_experiment, sample_model=sample_model, instrument_model=instrument_model, ) # We need to hack in the resolution model from the vanadium analysis, # since the setters and getters overwrite the model. This will be fixed # asap. diffusion_analysis.instrument_model._resolution_model = ( vanadium_analysis.instrument_model.resolution_model ) # We fix all parameters of the resolution model. diffusion_analysis.instrument_model.resolution_model.fix_all_parameters()

In [11]:

# Let us see how good the starting parameters are
diffusion_analysis.plot_data_and_model()

# Let us see how good the starting parameters are diffusion_analysis.plot_data_and_model()

Out[11]:

In [12]:

# Now we fit the data and plot the result. Looks good!
diffusion_analysis.fit(fit_method='independent')
diffusion_analysis.plot_data_and_model()

# Now we fit the data and plot the result. Looks good! diffusion_analysis.fit(fit_method='independent') diffusion_analysis.plot_data_and_model()

Out[12]:

In [13]:

# Let us look at the most interesting fit parameters
diffusion_analysis.plot_parameters(names=['Lorentzian width', 'Lorentzian area'])

# Let us look at the most interesting fit parameters diffusion_analysis.plot_parameters(names=['Lorentzian width', 'Lorentzian area'])

Out[13]:

In [14]:

# It will be possible to fit this to a DiffusionModel, but that will
# come later.

# It will be possible to fit this to a DiffusionModel, but that will # come later.

In [15]:

# Let us now fit directly to a diffusion model. We replace the
# Lorentzian with a Brownian translational diffusion model and keep the
# other parameters the same.
delta_function = DeltaFunction(display_name='DeltaFunction', area=0.2)
component_collection = ComponentCollection(
    components=[delta_function],
)
diffusion_model = BrownianTranslationalDiffusion(
    display_name='Brownian Translational Diffusion', diffusion_coefficient=2.4e-9, scale=0.5
)

sample_model = SampleModel(
    components=component_collection,
    diffusion_models=diffusion_model,
)

background_model = BackgroundModel(components=Polynomial(coefficients=[0.001]))

instrument_model = InstrumentModel(
    background_model=background_model,
)

diffusion_model_analysis = Analysis(
    display_name='Diffusion Full Analysis',
    experiment=diffusion_experiment,
    sample_model=sample_model,
    instrument_model=instrument_model,
)

# We again need to hack in the resolution model from the vanadium
# analysis, since the setters and getters overwrite the model. This will
# be fixed asap.
diffusion_model_analysis.instrument_model._resolution_model = (
    vanadium_analysis.instrument_model.resolution_model
)
diffusion_model_analysis.instrument_model.resolution_model.fix_all_parameters()

# Let us see how good the starting parameters are
diffusion_model_analysis.plot_data_and_model()

# Let us now fit directly to a diffusion model. We replace the # Lorentzian with a Brownian translational diffusion model and keep the # other parameters the same. delta_function = DeltaFunction(display_name='DeltaFunction', area=0.2) component_collection = ComponentCollection( components=[delta_function], ) diffusion_model = BrownianTranslationalDiffusion( display_name='Brownian Translational Diffusion', diffusion_coefficient=2.4e-9, scale=0.5 ) sample_model = SampleModel( components=component_collection, diffusion_models=diffusion_model, ) background_model = BackgroundModel(components=Polynomial(coefficients=[0.001])) instrument_model = InstrumentModel( background_model=background_model, ) diffusion_model_analysis = Analysis( display_name='Diffusion Full Analysis', experiment=diffusion_experiment, sample_model=sample_model, instrument_model=instrument_model, ) # We again need to hack in the resolution model from the vanadium # analysis, since the setters and getters overwrite the model. This will # be fixed asap. diffusion_model_analysis.instrument_model._resolution_model = ( vanadium_analysis.instrument_model.resolution_model ) diffusion_model_analysis.instrument_model.resolution_model.fix_all_parameters() # Let us see how good the starting parameters are diffusion_model_analysis.plot_data_and_model()

Out[15]:

In [16]:

# We now fit all the data simultaneously to the diffusion model, then
# plot the result. Looks good.
diffusion_model_analysis.fit(fit_method='simultaneous')
diffusion_model_analysis.plot_data_and_model()

# We now fit all the data simultaneously to the diffusion model, then # plot the result. Looks good. diffusion_model_analysis.fit(fit_method='simultaneous') diffusion_model_analysis.plot_data_and_model()

Out[16]:

In [17]:

# Let us look at the fitted diffusion coefficient
diffusion_model.get_all_parameters()

# Let us look at the fitted diffusion coefficient diffusion_model.get_all_parameters()

Out[17]:

[<Parameter 'diffusion_coefficient': 1.126e-08 ± 9.799e-11 m^2/s, bounds=[0.0:inf]>,
 <Parameter 'scale': 0.6938 ± 0.0043, bounds=[0.0:inf]>]