function update_weights_matrix, get_residual_matrix, reconstruct_data

diffpy/snmf/subroutines.py

@@ -1,4 +1,6 @@
 import numpy as np
+from diffpy.snmf.optimizers import get_weights
+from diffpy.snmf.factorizers import lsqnonneg
 
 
 def objective_function(residual_matrix, stretching_factor_matrix, smoothness, smoothness_term, component_matrix,
@@ -75,3 +77,167 @@ def get_stretched_component(stretching_factor, component, signal_length):
     component = np.asarray(component)
     normalized_grid = np.arange(signal_length)
     return np.interp(normalized_grid / stretching_factor, normalized_grid, component, left=0, right=0)
+
+
+def update_weights_matrix(component_amount, signal_length, stretching_factor_matrix, component_matrix, data_input,
+                          moment_amount, weights_matrix, method):
+    """Update the weight factors matrix.
+
+    Parameters
+    ----------
+    component_amount: int
+      The number of component signals the user would like to determine from the experimental data.
+
+    signal_length: int
+      The length of the experimental signal patterns
+
+    stretching_factor_matrix: 2d array like
+      The matrx containing the stretching factors of the calculated component signals. Has dimensions K x M where K is
+      the number of component signals and M is the number of XRD/PDF patterns.
+
+    component_matrix: 2d array like
+      The matrix containing the unstretched calculated component signals. Has dimensions N x K where N is the length of
+      the signals and K is the number of component signals.
+
+    data_input: 2d array like
+      The experimental series of PDF/XRD patterns. Has dimensions N x M where N is the length of the PDF/XRD signals and
+      M is the number of PDF/XRD patterns.
+
+    moment_amount: int
+      The number of PDF/XRD patterns from the experimental data.
+
+    weights_matrix: 2d array like
+      The matrix containing the weights of the stretched component signals. Has dimensions K x M where K is the number
+      of component signals and M is the number of XRD/PDF patterns.
+
+    method: str
+      The string specifying the method for obtaining individual weights.
+
+    Returns
+    -------
+    2d array like
+      The matrix containing the new weight factors of the stretched component signals.
+
+    """
+    stretching_factor_matrix = np.asarray(stretching_factor_matrix)
+    component_matrix = np.asarray(component_matrix)
+    data_input = np.asarray(data_input)
+    weights_matrix = np.asarray(weights_matrix)
+    weight = np.zeros(component_amount)
+    for i in range(moment_amount):
+        stretched_components = np.zeros((signal_length, component_amount))
+        for n in range(component_amount):
+            stretched_components[:, n] = get_stretched_component(stretching_factor_matrix[n, i], component_matrix[:, n],
+                                                                 signal_length)
+        if method == 'align':
+            weight = lsqnonneg(stretched_components[0:signal_length, :], data_input[0:signal_length, i])
+        else:
+            weight = get_weights(
+                stretched_components[0:signal_length, :].T @ stretched_components[0:signal_length, :],
+                -1 * stretched_components[0:signal_length, :].T @ data_input[0:signal_length, i],
+                0, 1)
+        weights_matrix[:, i] = weight
+    return weights_matrix
+
+
+def get_residual_matrix(component_matrix, weights_matrix, stretching_matrix, data_input, moment_amount,
+                        component_amount, signal_length):
+    """Obtains the residual matrix between the experimental data and calculated data
+
+    Calculates the difference between the experimental data and the reconstructed experimental data created from the
+    calculated components, weights, and stretching factors.
+
+    Parameters
+    ----------
+    component_matrix: 2d array like
+      The matrix containing the calculated component signals. Has dimensions N x K where N is the length of the signal
+      and K is the number of calculated component signals.
+
+    weights_matrix: 2d array like
+      The matrix containing the calculated weights of the stretched component signals. Has dimensions K x M where K is
+      the number of components and M is the number of moments or experimental PDF/XRD patterns.
+
+    stretching_matrix: 2d array like
+      The matrix containing the calculated stretching factors of the calculated component signals. Has dimensions K x M
+      where K is the number of components and M is the number of moments or experimental PDF/XRD patterns.
+
+    data_input: 2d array like
+      The matrix containing the experimental PDF/XRD data. Has dimensions N x M where N is the length of the signals and
+      M is the number of signal patterns.
+
+    moment_amount: int
+      The number of PDF/XRD patterns.
+
+    component_amount: int
+      The number of component signals that user would like to experimental data.
+
+    signal_length: int
+      The length of the signals.
+
+
+    Returns
+    -------
+    2d array like
+
+    """
+    component_matrix = np.asarray(component_matrix)
+    weights_matrix = np.asarray(weights_matrix)
+    stretching_matrix = np.asarray(stretching_matrix)
+    data_input = np.asarray(data_input)
+    residual_matrx = -1 * data_input
+    for m in range(moment_amount):
+        residual = residual_matrx[:, m]
+        for k in range(component_amount):
+            residual = residual + weights_matrix[k, m] * get_stretched_component(stretching_matrix[k, m],
+                                                                                 component_matrix[:, k], signal_length)
+        residual_matrx[:, m] = residual
+    return residual_matrx
+
+
+def reconstruct_data(stretching_factor_matrix, component_matrix, weight_matrix, component_amount,
+                     moment_amount, signal_length):
+    """Reconstructs the experimental data from the component signals, stretching factors, and weights.
+
+    Parameters
+    ----------
+    stretching_factor_matrix: 2d array like
+      The matrix containing the stretching factors of the component signals.
+
+    component_matrix: 2d array like
+      The matrix containing the unstretched component signals
+
+    weight_matrix: 2d array like
+
+    component_amount: int
+
+    moment_amount: int
+
+    signal_length: int
+
+    Returns
+    -------
+    2d array like
+
+    """
+    stretching_factor_matrix = np.asarray(stretching_factor_matrix)
+    component_matrix = np.asarray(component_matrix)
+    weight_matrix = np.asarray(weight_matrix)
+    stretched_component_series = []
+    for moment in range(moment_amount):
+        for component in range(component_amount):
+            stretched_component = get_stretched_component(stretching_factor_matrix[component, moment],
+                                                          component_matrix[:, component], signal_length)
+            stretched_component_series.append(stretched_component)
+            print(stretched_component)
+    stretched_component_series = np.column_stack(stretched_component_series)
+    print(stretched_component_series)
+
+    reconstructed_data = []
+    moment = 0
+    for s_component in range(0, moment_amount * component_amount, component_amount):
+        block = stretched_component_series[:, s_component:s_component + component_amount]
+        for component in range(component_amount):
+            block[:, component] = block[:, component] * weight_matrix[component, moment]
+        reconstructed_data.append(block)
+        moment += 1
+    return np.column_stack(reconstructed_data)

diffpy/snmf/tests/test_subroutines.py

@@ -1,6 +1,7 @@
 import pytest
 import numpy as np
-from diffpy.snmf.subroutines import objective_function, get_stretched_component
+from diffpy.snmf.subroutines import objective_function, get_stretched_component, reconstruct_data, get_residual_matrix, \
+    update_weights_matrix
 
 to = [
     ([[[1, 2], [3, 4]], [[5, 6], [7, 8]], 1e11, [[1, 2], [3, 4]], [[1, 2], [3, 4]], 1], 2.574e14),
@@ -39,3 +40,54 @@ def test_get_stretched_component(tgso):
     actual = get_stretched_component(tgso[0][0], tgso[0][1], tgso[0][2])
     expected = tgso[1]
     np.testing.assert_allclose(actual, expected, rtol=1e-03)
+
+
+tuwm = [([2, 2, [[.5, .6], [.7, .8]], [[1, 2], [4, 8]], [[1.6, 2.8], [5, 8.8]], 2, [[.78, .12], [.5, .5]], None],
+         [[0, 1], [1, 1]])
+
+        ]
+
+
+@pytest.mark.parametrize('tuwm', tuwm)
+def test_update_weights_matrix(tuwm):
+    actual = update_weights_matrix(tuwm[0][0], tuwm[0][1], tuwm[0][2], tuwm[0][3], tuwm[0][4], tuwm[0][5], tuwm[0][6],
+                                   tuwm[0][7])
+    expected = tuwm[1]
+    np.testing.assert_allclose(actual, expected)
+
+
+tgrm = [
+    ([[[1, 2], [3, 4]], [[.25], [.75]], [[.9], [.7]], [[11, 22], [33, 44]], 1, 2, 2], [[-9.25, -22], [-30.6190, -44]]),
+    ([[[1, 2], [3, 4]], [[1], [1]], [[1], [1]], [[11, 22], [33, 44]], 1, 2, 2], [[-8, -22], [-26, -44]])
+    # positive whole numbers
+    # negative whole numbers
+    # single component
+    # positive float
+    # negative floats
+
+]
+
+
+@pytest.mark.parametrize('tgrm', tgrm)
+def test_get_residual_matrix(tgrm):
+    actual = get_residual_matrix(tgrm[0][0], tgrm[0][1], tgrm[0][2], tgrm[0][3], tgrm[0][4], tgrm[0][5], tgrm[0][6])
+    expected = tgrm[1]
+    np.testing.assert_allclose(actual, expected)
+
+
+trd = [
+    ([np.array([[.5], [.5]]), np.array([[1, 2], [3, 4]]), np.array([[.25], [.75]]), 2, 1, 2],
+     np.array([[.25, 1.5], [0, 0]])),
+    ([[[.9], [.7]], [[1, 2], [3, 4], [5, 6], [7, 8], [9, 10]], [[.25], [.75]], 2, 1, 5],
+     [[.25, 1.5], [.8056, 3.6429], [1.3611, 5.7857], [1.9167, 5.3571], [1.25, 0]]),
+    ([[[.5, .7], [.5, .8]], [[1.1, 2.2], [3.3, 4.4]], [[.25, .5], [.75, .5]], 2, 2, 2],
+     [[.275, .55, 1.65, 1.1], [0, .9429, 0, 1.65]])
+
+]
+
+
+@pytest.mark.parametrize('trd', trd)
+def test_reconstruct_data(trd):
+    actual = reconstruct_data(trd[0][0], trd[0][1], trd[0][2], trd[0][3], trd[0][4], trd[0][5])
+    expected = trd[1]
+    np.testing.assert_allclose(actual, expected)