From c59cc7048d9c87ab77d7bf0f2796175e4144d08d Mon Sep 17 00:00:00 2001 From: John Halloran Date: Tue, 19 Aug 2025 22:06:07 -0700 Subject: [PATCH 1/3] refactor: reconstruct separate from residuals --- src/diffpy/snmf/snmf_class.py | 65 ++++++++++++++++++++++++----------- 1 file changed, 44 insertions(+), 21 deletions(-) diff --git a/src/diffpy/snmf/snmf_class.py b/src/diffpy/snmf/snmf_class.py index 2983e57..68ec23a 100644 --- a/src/diffpy/snmf/snmf_class.py +++ b/src/diffpy/snmf/snmf_class.py @@ -304,12 +304,7 @@ def outer_loop(self): def get_residual_matrix(self, components=None, weights=None, stretch=None): """ - Return the residuals (difference) between the source matrix and its reconstruction - from the given components, weights, and stretch factors. - - Each component profile is stretched, interpolated to fractional positions, - weighted per signal, and summed to form the reconstruction. The residuals - are the source matrix minus this reconstruction. + Return the residuals (difference) between the source matrix and its reconstruction. Parameters ---------- @@ -329,21 +324,8 @@ def get_residual_matrix(self, components=None, weights=None, stretch=None): if stretch is None: stretch = self.stretch - residuals = -self.source_matrix.copy() - sample_indices = np.arange(components.shape[0]) # (signal_len,) - - for comp in range(components.shape[1]): # loop over components - residuals += ( - np.interp( - sample_indices[:, None] - / stretch[comp][None, :], # fractional positions (signal_len, n_signals) - sample_indices, # (signal_len,) - components[:, comp], # component profile (signal_len,) - left=components[0, comp], - right=components[-1, comp], - ) - * weights[comp][None, :] # broadcast (n_signals,) over rows - ) + reconstructed_matrix = reconstruct_matrix(components, weights, stretch) + residuals = reconstructed_matrix - self.source_matrix return residuals @@ -718,3 +700,44 @@ def cubic_largest_real_root(p, q): y = np.max(real_roots, axis=0) * (delta < 0) # Keep only real roots when delta < 0 return y + + +def reconstruct_matrix(components=None, weights=None, stretch=None): + """ + Construct the approximation of the source matrix corresponding to the + given components, weights, and stretch factors. + + Each component profile is stretched, interpolated to fractional positions, + weighted per signal, and summed to form the reconstruction. + + Parameters + ---------- + components : (signal_len, n_components) array + weights : (n_components, n_signals) array + stretch : (n_components, n_signals) array + + Returns + ------- + reconstructed_matrix : (signal_len, n_signals) array + """ + + signal_len = components.shape[0] + n_signals = weights.shape[1] + n_components = components.shape[1] + + reconstructed_matrix = np.zeros((signal_len, n_signals)) + sample_indices = np.arange(signal_len) + + for comp in range(n_components): # loop over components + reconstructed_matrix += ( + np.interp( + sample_indices[:, None] / stretch[comp][None, :], # fractional positions (signal_len, n_signals) + sample_indices, # (signal_len,) + components[:, comp], # component profile (signal_len,) + left=components[0, comp], + right=components[-1, comp], + ) + * weights[comp][None, :] # broadcast (n_signals,) over rows + ) + + return reconstructed_matrix From 98db94c74ee1e60b98f00fe577e1bad843b41fca Mon Sep 17 00:00:00 2001 From: John Halloran Date: Tue, 19 Aug 2025 22:52:46 -0700 Subject: [PATCH 2/3] fix: inputs to reconstruct_matrix should not be optional --- src/diffpy/snmf/snmf_class.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/diffpy/snmf/snmf_class.py b/src/diffpy/snmf/snmf_class.py index 68ec23a..3934267 100644 --- a/src/diffpy/snmf/snmf_class.py +++ b/src/diffpy/snmf/snmf_class.py @@ -702,7 +702,7 @@ def cubic_largest_real_root(p, q): return y -def reconstruct_matrix(components=None, weights=None, stretch=None): +def reconstruct_matrix(components, weights, stretch): """ Construct the approximation of the source matrix corresponding to the given components, weights, and stretch factors. From 4c058886cf9ca0b813900549c5b6751d1fbc2224 Mon Sep 17 00:00:00 2001 From: John Halloran Date: Tue, 19 Aug 2025 23:02:53 -0700 Subject: [PATCH 3/3] docs: updated default documentation for init_weights --- src/diffpy/snmf/snmf_class.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/diffpy/snmf/snmf_class.py b/src/diffpy/snmf/snmf_class.py index 3934267..e0e1e8c 100644 --- a/src/diffpy/snmf/snmf_class.py +++ b/src/diffpy/snmf/snmf_class.py @@ -83,7 +83,7 @@ def __init__( ---------- source_matrix : ndarray The data to be decomposed. Shape is (length_of_signal, number_of_conditions). - init_weights : ndarray Optional Default = rng.beta(a=2.5, b=1.5, size=(n_components, n_signals)) + init_weights : ndarray Optional Default = rng.beta(a=2.0, b=2.0, size=(n_components, n_signals)) The initial guesses for the component weights at each stretching condition. Shape is (number_of_components, number_of_signals) Must provide exactly one of this or n_components.