Merge pull request #46 from aajayi-21/update_weights2

function update_weights

Author: sbillinge

diffpy/snmf/subroutines.py

@@ -134,11 +134,51 @@ def construct_weight_matrix(components):
         raise ValueError(f"Number of components = {number_of_components}. Number of components must be >= 1")
     if number_of_signals == 0:
         raise ValueError(f"Number of signals = {number_of_signals}. Number_of_signals must be >= 1.")
-    weights_matrix = np.zeros((number_of_components,number_of_signals))
+    weights_matrix = np.zeros((number_of_components, number_of_signals))
     for i, component in enumerate(components):
         weights_matrix[i] = component.weights
     return weights_matrix
 
+
+def update_weights(components, data_input, method=None):
+    """Updates the weights matrix.
+
+    Updates the weights matrix and the weights vector for each ComponentSignal object.
+
+    Parameters
+    ----------
+    components: tuple of ComponentSignal objects
+      The tuple containing the component signals.
+    method: str
+      The string specifying which method should be used to find a new weight matrix: non-negative least squares or a
+      quadratic program.
+    data_input: 2d array
+      The 2d array containing the user-provided signals.
+
+    Returns
+    -------
+    2d array
+      The 2d array containing the weight factors for each component for each signal from `data_input`. Has dimensions
+      K x M where K is the number of components and M is the number of signals in `data_input.`
+    """
+    data_input = np.asarray(data_input)
+    weight_matrix = construct_weight_matrix(components)
+    number_of_signals = len(components[0].weights)
+    number_of_components = len(components)
+    signal_length = len(components[0].grid)
+    for signal in range(number_of_signals):
+        stretched_components = np.zeros((signal_length, number_of_components))
+        for i, component in enumerate(components):
+            stretched_components[:, i] = component.apply_stretch(signal)[0]
+        if method == 'align':
+            weights = lsqnonneg(stretched_components, data_input[:,signal])
+        else:
+            weights = get_weights(stretched_components.T @ stretched_components,
+                                  -stretched_components.T @ data_input[:, signal], 0, 1)
+            weight_matrix[:, signal] = weights
+    return weight_matrix
+
+
 def initialize_arrays(number_of_components, number_of_moments, signal_length):
     """Generates the initial guesses for the weight, stretching, and component matrices
 

diffpy/snmf/tests/test_subroutines.py

@@ -2,7 +2,8 @@
 import numpy as np
 from diffpy.snmf.containers import ComponentSignal
 from diffpy.snmf.subroutines import objective_function, get_stretched_component, reconstruct_data, get_residual_matrix, \
-    update_weights_matrix, initialize_arrays, lift_data, initialize_components, construct_stretching_matrix, construct_component_matrix, construct_weight_matrix
+    update_weights_matrix, initialize_arrays, lift_data, initialize_components, construct_stretching_matrix, \
+    construct_component_matrix, construct_weight_matrix, update_weights
 
 to = [
     ([[[1, 2], [3, 4]], [[5, 6], [7, 8]], 1e11, [[1, 2], [3, 4]], [[1, 2], [3, 4]], 1], 2.574e14),
@@ -207,6 +208,7 @@ def test_construct_component_matrix(tccm):
     for component in tccm:
         np.testing.assert_allclose(actual[component.id], component.iq)
 
+
 tcwm = [
         ([ComponentSignal([0,.25,.5,.75,1],20,0)]),
         # ([ComponentSignal([0,.25,.5,.75,1],0,0)]), # 0 signal length. Failure expected
@@ -225,3 +227,28 @@ def test_construct_weight_matrix(tcwm):
     actual = construct_weight_matrix(tcwm)
     for component in tcwm:
         np.testing.assert_allclose(actual[component.id], component.weights)
+
+
+tuw = [([ComponentSignal([0, .25, .5, .75, 1], 2, 0), ComponentSignal([0, .25, .5, .75, 1], 2, 1),
+         ComponentSignal([0, .25, .5, .75, 1], 2, 2)], [[1, 1], [1.2, 1.3], [1.3, 1.4], [1.4, 1.5], [2, 2.1]], None),
+       ([ComponentSignal([0, .25, .5, .75, 1], 2, 0), ComponentSignal([0, .25, .5, .75, 1], 2, 1),
+         ComponentSignal([0, .25, .5, .75, 1], 2, 2)], [[1, 1], [1.2, 1.3], [1.3, 1.4], [1.4, 1.5], [2, 2.1]], "align"),
+       ([ComponentSignal([0, .25, .5, .75, 1], 2, 0), ComponentSignal([0, .25, .5, .75, 1], 2, 1),
+         ComponentSignal([0, .25, .5, .75, 1], 2, 2)], [[0, 0], [0, 0], [0, 0], [0, 0], [0, 0]], None),
+       ([ComponentSignal([0, .25, .5, .75, 1], 2, 0), ComponentSignal([0, .25, .5, .75, 1], 2, 1),
+         ComponentSignal([0, .25, .5, .75, 1], 2, 2)], [[0, 0], [0, 0], [0, 0], [0, 0], [0, 0]], "align"),
+       ([ComponentSignal([0, .25, .5, .75, 1], 2, 0), ComponentSignal([0, .25, .5, .75, 1], 2, 1),
+         ComponentSignal([0, .25, .5, .75, 1], 2, 2)], [[-.5, 1], [1.2, -1.3], [1.1, -1], [0, -1.5], [0, .1]], None),
+       ([ComponentSignal([0, .25, .5, .75, 1], 2, 0), ComponentSignal([0, .25, .5, .75, 1], 2, 1),
+         ComponentSignal([0, .25, .5, .75, 1], 2, 2)], [[-.5, 1], [1.2, -1.3], [1.1, -1], [0, -1.5], [0, .1]], "align"),
+       # ([ComponentSignal([0, .25, .5, .75, 1], 0, 0), ComponentSignal([0, .25, .5, .75, 1], 0, 1),
+         # ComponentSignal([0, .25, .5, .75, 1], 0, 2)], [[1, 1], [1.2, 1.3], [1.3, 1.4], [1.4, 1.5], [2, 2.1]], None),
+       # ([ComponentSignal([0, .25, .5, .75, 1], 0, 0), ComponentSignal([0, .25, .5, .75, 1], 0, 1),
+         # ComponentSignal([0, .25, .5, .75, 1], 0, 2)], [], None),
+       # ([ComponentSignal([0, .25, .5, .75, 1], 2, 0), ComponentSignal([0, .25, .5, .75, 1], 2, 1),
+         # ComponentSignal([0, .25, .5, .75, 1], 2, 2)], [], 170),
+       ]
+@pytest.mark.parametrize('tuw', tuw)
+def test_update_weights(tuw):
+    actual = update_weights(tuw[0], tuw[1], tuw[2])
+    assert np.shape(actual) == (len(tuw[0]), len(tuw[0][0].weights))