Merge 9f67e54 into 4158f08

skxray/core/correlation.py

@@ -51,6 +51,7 @@
 import numpy as np
 
 import skxray.core.utils as core
+from lmfit import minimize, Parameters
 
 logger = logging.getLogger(__name__)
 
@@ -84,7 +85,7 @@ def multi_tau_auto_corr(num_levels, num_bufs, labels, images):
     Returns
     -------
     g2 : array
-        matrix of one-time correlation
+        matrix of normalized intensity-intensity autocorrelation
         shape (num_levels, number of labels(ROI))
 
     lag_steps : array
@@ -94,6 +95,18 @@ def multi_tau_auto_corr(num_levels, num_bufs, labels, images):
     Note
     ----
 
+    The normalized intensity-intensity time-autocorrelation function
+    is defined as
+
+    :math ::
+    g2(Q, t) = \\frac{<I(Q, t)I(Q, t + delay)> }{<I(Q, t)>^2}
+
+    ; delay > 0
+
+    Here, I(Q, t) refers to the scattering strength at the momentum
+    transfer vector Q in reciprocal space at time t, and the brackets
+    <...> refer to averages over time t.
+
     This implementation is based on code in the language Yorick
     by Mark Sutton, based on published work. [1]_
 
@@ -367,3 +380,130 @@ def extract_label_indices(labels):
     label_mask = labels[labels > 0]
 
     return label_mask, pixel_list
+
+
+def auto_corr_scat_factor(lags, beta, relaxation_rate, baseline=1):
+    """
+    This model will provide normalized intensity-intensity time
+    correlation data to be minimized.
+
+    Parameters
+    ----------
+    lags : array
+        delay time
+
+    beta : float
+        optical contrast (speckle contrast), a sample-independent
+        beamline parameter
+
+    relaxation_rate : float
+        relaxation time associated with the samples dynamics.
+
+    baseline : float, optional
+        baseline of one time correlation
+        equal to one for ergodic samples
+
+    Returns
+    -------
+    g2 : array
+        normalized intensity-intensity time autocorreltion
+
+    Note :
+    The intensity-intensity autocorrelation g2 is connected to the intermediate
+    scattering factor(ISF) g1
+
+    :math ::
+        g2(q, t) = \\beta[g1(q, t)]^{2} + g_\\infty
+
+    For a system undergoing  diffusive dynamics,
+
+    :math ::
+        g1(q, t) = e^{-\\Gamma t}
+
+    :math ::
+        g2(q, t) = \\beta e^{-2\\Gamma t} + g_\\infty
+
+    These implementation are based on published work. [1]_
+
+     References
+    ----------
+
+    .. [1] L. Li, P. Kwasniewski, D. Orsi, L. Wiegart, L. Cristofolini, C. Caronna
+       and A. Fluerasu, " Photon statistics and speckle visibility spectroscopy with
+       partially coherent X-rays," J. Synchrotron Rad. vol 21, p 1288-1295, 2014
+
+    """
+    return beta*np.exp(-2*relaxation_rate*lags) + baseline
+
+
+def _residual_auto_corr(params, lags, g2_data, *eps_data):
+    """
+    This will provide difference between experiment data and the model
+    Parameters
+    ----------
+    params : dict or Parameters
+        beta - float, optical contrast (speckle contrast),
+        relaxation_rate - float, relaxation time associated with the
+        samples dynamics,
+        baseline -float, optional baseline of one time
+        correlation equal to one for ergodic samples
+        have to give as either dictionary or Parameters
+        eg: One of the following
+        #create a dictionary
+        {'beta': 0.2, 'relaxation_rate':10, 'baseline':1}
+
+        # create a set of Parameters
+        from lmfit import Parameters
+        params = Parameters()
+        params.add('beta',  value=0.2,  min=0, max=0.3)
+        params.add('relaxation_rate', value=10, min=0, max=11)
+        params.add('baseline', value=1)
+
+    lags : array
+        delay time
+
+    g2_data : array
+        normalized intensity-intensity time autocorreltion
+
+    eps_data : array, optional
+        standard error of the normalized intensity-intensity
+        time autocorreltion
+
+    Returns
+    -------
+    residual : array
+        difference between experimental result and the model
+    """
+    # create set of parameters
+    beta = params['beta'].value
+    relaxation_rate = params['relaxation_rate'].value
+    baseline = params['baseline'].value
+
+    return (g2_data - auto_corr_scat_factor(lags, beta, relaxation_rate,
+                                            baseline=1))/eps_data
+
+
+def fit_auto_corr(params, x, data, eps_data):
+    """
+    Minimize the function and calculate final result
+
+    Parameters
+    ----------
+    params: dict
+        parameters dictionary
+
+    x : array
+       x
+
+    data : array
+
+    eps_data : array
+
+    Returns
+    -------
+    final_result : array
+        minimized fitting to results using least square model
+
+    """
+    result = minimize(_residual_auto_corr, params, args=(x, data, eps_data))
+    return data + result.residual

skxray/core/tests/test_correlation.py

@@ -38,11 +38,18 @@
 import numpy as np
 from numpy.testing import (assert_array_almost_equal,
                            assert_almost_equal)
-from skimage import data
+
+from nose.tools import assert_equal, assert_true, assert_raises
 
 import skxray.core.correlation as corr
 import skxray.core.roi as roi
-from skxray.testing.decorators import skip_if
+import skxray.core.utils as core
+
+from skxray.testing.decorators import known_fail_if, skip_if
+import numpy.testing as npt
+
+from skimage import data
+from lmfit import Parameters
 
 logger = logging.getLogger(__name__)
 
@@ -90,3 +97,46 @@ def test_image_stack_correlation():
 
     assert_almost_equal(True, np.all(g2[:, 0], axis=0))
     assert_almost_equal(True, np.all(g2[:, 1], axis=0))
+
+    num_buf = 5
+
+    # check the number of buffers are even
+    assert_raises(ValueError,
+                  lambda: corr.multi_tau_auto_corr(num_levels, num_buf,
+                                                   coins_mesh, coins_stack))
+    # check image shape and labels shape are equal
+    #assert_raises(ValueError,
+    #            lambda : corr.multi_tau_auto_corr(num_levels, num_bufs,
+    #                                                indices, coins_stack))
+
+    # check the number of pixels is zero
+    mesh = np.zeros_like(coins)
+    assert_raises(ValueError,
+                  lambda: corr.multi_tau_auto_corr(num_levels, num_bufs,
+                                                   mesh, coins_stack))
+
+
+def test_auto_corr_scat_factor():
+    num_levels, num_bufs = 3, 4
+    tot_channels, lags = core.multi_tau_lags(num_levels, num_bufs)
+    beta = 0.5
+    relaxation_rate = 10.0
+    baseline = 1.0
+
+    g2 = corr.auto_corr_scat_factor(lags, beta, relaxation_rate, baseline)
+
+    assert_array_almost_equal(g2, np.array([1.5, 1.0, 1.0, 1.0, 1.0,
+                                            1.0, 1.0, 1.0]), decimal=8)
+
+
+def fit_auto_corr():
+    params = Parameters()
+    params.add('beta',  value=0.1699, min=0.167, max=0.21556)
+    params.add('relaxation_rate', value=6.159, min=6.158, max=6.197)
+    params.add('baseline', value=1, min=0.8, max=1.0)
+
+    num_levels, num_bufs = 2, 4
+    tot_channels, lags = core.multi_tau_lags(num_levels, num_bufs)
+    data = np.array([1.216, 1.212, 1.208, 1.204, 1.196])
+
+    fit_result = corr.fit_auto_corr(params, lags[1:], data, eps_data=1)