Simple walkthrough#

The heart of the DDM code (found in the ddm_calc.py file) is the computation of the image structure function. This is found by taking the average of the Fourier transforms of all image differences. By “image differences,” I mean the result of subtracting two images separated by a given lag time \(\Delta t\).

To describe the process mathematically, we find the difference between images separated by some lag time \(\Delta t\):

\[\Delta I = I(x,y;t) - I(x,y;t + \Delta t)\]

For a given \(\Delta t\) all such image differences are calculated. We then Fourier transform each \(\Delta I\) and average all of the same \(\Delta t\).

This results in the image structure function \(D(q_x,q_y,\Delta t)\).

Importing the necessary modules#

Modules for plotting#

We use matplotlib for creating figures and plots. Note that we use:

%matplotlib inline

This sets the backend of matplotlib to inline which means the plots are included in the notebook. If you want the plots to also be interactive (e.g., having the ability to zoom, scroll, and save) then use:

%matplotlib notebook

[1]:

%matplotlib inline
import matplotlib
import matplotlib.pyplot as plt

Modules for numerical work#

Here, we import numpy and xarray.

Note that xarray might not have come with your Anaconda Python distribution (or whichever other distribution you installed). If that is the case, you’ll need to install this package.

[2]:

import numpy as np #numerical python used for working with arrays, mathematical operations
import xarray as xr #package for labeling and adding metadata to multi-dimensional arrays

Import the PyDDM package#

Make sure you append to sys.path the directory containing the PyDDM code.

[3]:

import sys
sys.path.append("../../PyDDM") #must point to the PyDDM folder
import ddm_analysis_and_fitting as ddm

Initiazing DDM_Fit class and fitting our data to a model#

Next, we initialize the DDM_Fit class which will help us fit the calculated DDM matrix to a model of our choosing.

When the DDM_Fit class is initialized, it will display a table of the parameters, with the initial guesses and bounds, that were in the parameters of the yaml file. It will also indicate if, for the chosen model, initial guesses for parameters are missing.

Note: The intial guesses and bounds for the parameters are given in the yaml file. But you may realize that you need to adjust the bounds. You could change the yaml file and reload it. But you can also change the initial guesses and bounds for the parameters as shown in this code snippet:

ddm_fit.set_parameter_initial_guess('Tau', 0.5)
ddm_fit.set_parameter_bounds('Tau', [0.001, 1000])

Also, note that to initialize the DDM_Fit class, we need to pass it the YAML file which contains the metadata and parameters. Below, we use the fact that the path to the YAML file is stored in the variable ddm_calc.data_yaml. But we could also use the path to the YAML file just like we did when initializing the DDM_Analysis class. That is, this would work as well:

ddm_fit = ddm.DDM_Fit("../../Examples/example_data_silica_beads.yml")

[6]:

ddm_fit = ddm.DDM_Fit(ddm_calc.data_yaml)

	Initial guess	Minimum	Maximum
Amplitude	100.0	1.000	1000000.0
Tau	1.0	0.001	10.0
Background	25000.0	0.000	10000000.0
StretchingExp	1.0	0.500	1.1

Loading file C:/Users/rmcgorty/Documents/GitHub/PyDDM/Examples/images_nobin_40x_128x128_8bit_ddmmatrix.nc ...

Below, the fit is performed and saved into an xarray Dataset. We also have a table of the parameters for each q. Since this table may be long (especially if you have a large image size), you might want to not display it. If that’s the case, then run this command with the display_table parmeter set to False, i.e.:

fit01 = ddm_fit.fit(name_fit = 'fit01', display_table=False)

[7]:

#Note that the method `fit` has many optional parameters, but here we
# allow those optional parameter to take their default values

fit01 = ddm_fit.fit(name_fit = 'fit01')

In function 'get_tau_vs_q_fit', using new tau...
Fit is saved in fittings dictionary with key 'fit01'.

	q	Amplitude	Tau	Background	StretchingExp
0	0.000000	100.000000	1.000000	25000.000000	1.000000
1	0.202840	2221.053486	10.000000	105.857478	0.696463
2	0.405681	3006.600450	10.000000	8.836142	0.572435
3	0.608521	3718.728783	4.503455	0.000000	0.687558
4	0.811362	5092.262467	3.367577	0.000000	0.758814
5	1.014202	5298.076962	1.495877	0.000000	0.875560
6	1.217043	5017.826391	0.923533	0.000000	0.919762
7	1.419883	5018.286724	0.707998	0.000000	0.933032
8	1.622723	4804.352014	0.486938	0.000000	0.973702
9	1.825564	5181.793843	0.446446	0.000000	0.965262
10	2.028404	5611.243988	0.371388	0.000000	0.979897
11	2.231245	6403.155522	0.308168	0.000000	0.996226
12	2.434085	7812.097726	0.273778	0.000000	1.024677
13	2.636926	8842.276777	0.231335	0.000000	1.049213
14	2.839766	9033.580795	0.205506	0.000000	1.034606
15	3.042607	9228.078528	0.190861	0.000000	1.036096
16	3.245447	8215.724811	0.166399	0.000000	1.037051
17	3.448287	7064.034828	0.153635	0.000000	1.040212
18	3.651128	5159.959987	0.139709	0.000000	1.032205
19	3.853968	3564.060652	0.127278	0.000000	1.037030
20	4.056809	2320.133466	0.118597	0.000000	1.026376
21	4.259649	1422.696297	0.111623	0.000000	0.998286
22	4.462490	866.841353	0.102422	0.000000	0.964719
23	4.665330	514.012003	0.098157	6.102412	0.915522
24	4.868170	301.203055	0.094404	18.732899	0.890644
25	5.071011	194.349793	0.101695	31.757336	0.885361
26	5.273851	152.163959	0.094730	28.431701	0.801239
27	5.476692	118.809414	0.095692	32.667078	0.789463
28	5.679532	96.618516	0.098254	34.735010	0.754149
29	5.882373	77.994077	0.109526	38.194460	0.746305
30	6.085213	73.573291	0.109010	37.473100	0.710879
31	6.288053	66.603954	0.101589	36.948208	0.686781

Trying a different mathematical model#

Above, we fit the DDM matrix to a model where the intermediate scattering function (ISF or \(f(q,\Delta t)\)) is an exponential. That is, we fit our data to:

\[D(q,\Delta t) = A(q) \left[ 1 - \exp \left(-\Delta t / \tau(q) \right) ^{s(q)} \right] + B(q)\]

where \(A\) is the amplitude, \(B\) is the background, \(\tau\) is the decay time, and \(s\) is the stretching exponent.

An alternative is to use other information to determine \(A(q)\) and \(B\) and then fit only the ISF. We will try this below.

[9]:

ddm_fit.reload_fit_model_by_name("ISF - Single Exponential")

	Initial guess	Minimum	Maximum
Tau	1.0	0.001	10.0
StretchingExp	1.0	0.500	1.1

[10]:

fit02 = ddm_fit.fit(name_fit = 'fit02', display_table=False)

In function 'get_tau_vs_q_fit', using new tau...
Fit is saved in fittings dictionary with key 'fit02'.

[11]:

ddm.fit_report(fit02, q_indices=[3,6,9,22], forced_qs=[4,16], use_new_tau=True, show=True)

In function 'get_tau_vs_q_fit', using new tau...
In hf.plot_one_tau_vs_q function, using new tau...

Switching method for calculating background#

[12]:

ddm_fit.ddm_dataset = ddm.recalculate_ISF_with_new_background(ddm_fit.ddm_dataset, background_method=5)

[13]:

ddm_fit.reload_fit_model_by_name("ISF - Single Exponential")

	Initial guess	Minimum	Maximum
Tau	1.0	0.001	10.0
StretchingExp	1.0	0.500	1.1

[14]:

fit03 = ddm_fit.fit(name_fit = 'fit03', display_table=False)

In function 'get_tau_vs_q_fit', using new tau...
Fit is saved in fittings dictionary with key 'fit03'.

[15]:

ddm.fit_report(fit03, q_indices=[3,6,9,22], forced_qs=[4,16], use_new_tau=True, show=True)

In function 'get_tau_vs_q_fit', using new tau...
In hf.plot_one_tau_vs_q function, using new tau...

Comparing fit methods#

[16]:

fig, ax = plt.subplots(figsize=(8,8/1.618))

q_lo, q_hi = 1,28
ax.loglog(fit01.q[q_lo:q_hi], fit01.parameters.loc['Tau'][q_lo:q_hi], 'o', c='tab:red', label="DDM Fit")
ax.loglog(fit02.q[q_lo:q_hi], fit02.parameters.loc['Tau'][q_lo:q_hi], 'o', c='tab:blue', label="ISF Fit")
ax.loglog(fit03.q[q_lo:q_hi], fit03.parameters.loc['Tau'][q_lo:q_hi], 'o', c='tab:green', label="ISF Fit (different B est)")

plt.legend(loc=0)

ax.set_xlabel("$q$ ($\mu$m$^{-1}$)")
ax.set_ylabel("Decay time, $\\tau$ (s)")

[16]:

Text(0, 0.5, 'Decay time, $\\tau$ (s)')

Interactive with matplotlib#

Using the Browse_DDM_Fits class as shown below, one can interactively inspect the fits (to either the DDM matrix or the ISF) by selecting the appropriate point on the \(\tau\) vs \(q\) plot or by pressing ‘N’ or ‘P’ for next or previous.

[17]:

%matplotlib notebook
fig, (ax, ax2) = plt.subplots(2, 1, figsize=(10,10/1.618))
browser = ddm.Browse_DDM_Fits(fig, ax, ax2, fit02)

fig.canvas.mpl_connect('pick_event', browser.on_pick)
fig.canvas.mpl_connect('key_press_event', browser.on_press)

ax.set_title('Decay time vs wavevector')
ax.set_xlabel("q")
ax.set_ylabel("tau (s)")

Click on a point in the tau vs q plot to see a fit.
Or press 'N' or 'P' to display next or previous fit.

[17]:

Text(0, 0.5, 'tau (s)')

[ ]:

Saving the results#

We have a few options for saving the results of the fit. We can save the information contained in the xarray Dataset fit01 as an Excel file. Or we can save it as a netCDF file using the xarray function `to_netcdf <https://xarray.pydata.org/en/stable/generated/xarray.Dataset.to_netcdf.html>`__. If we have the fit results as a netCDF file, then we can reload it easily with teh xarray function `open_dataset <https://xarray.pydata.org/en/stable/generated/xarray.open_dataset.html>`__.

[13]:

ddm.save_fit_results_to_excel(fit01)

[14]:

fit01.to_netcdf("fit_of_ddmmatrix.nc")

Simple walkthrough

Contents