From b909dcb3d44fadcadd245e3b03b65dd822e9edff Mon Sep 17 00:00:00 2001 From: John Halloran Date: Sat, 23 Aug 2025 21:58:48 -0700 Subject: [PATCH 1/5] fix: guard against zero/NaN stretches in apply_interpolation_matrix --- src/diffpy/snmf/snmf_class.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/diffpy/snmf/snmf_class.py b/src/diffpy/snmf/snmf_class.py index a62d021..da1b128 100644 --- a/src/diffpy/snmf/snmf_class.py +++ b/src/diffpy/snmf/snmf_class.py @@ -384,6 +384,9 @@ def apply_interpolation_matrix(self, components=None, weights=None, stretch=None if stretch is None: stretch = self.stretch_ + eps = 1e-8 # guard against divide by zero/NaN stretches + stretch = np.maximum(stretch, eps) + # Compute scaled indices stretch_flat = stretch.reshape(1, self.n_signals * self.n_components) ** -1 stretch_tiled = np.tile(stretch_flat, (self.signal_length, 1)) From a7adeedb513b178db86501c77e0febcbb9043178 Mon Sep 17 00:00:00 2001 From: John Halloran Date: Sat, 23 Aug 2025 22:09:53 -0700 Subject: [PATCH 2/5] refactor: use broadcasting instead of np.tile in apply_interpolation_matrix --- src/diffpy/snmf/snmf_class.py | 17 ++++++----------- 1 file changed, 6 insertions(+), 11 deletions(-) diff --git a/src/diffpy/snmf/snmf_class.py b/src/diffpy/snmf/snmf_class.py index da1b128..c7c0516 100644 --- a/src/diffpy/snmf/snmf_class.py +++ b/src/diffpy/snmf/snmf_class.py @@ -389,17 +389,12 @@ def apply_interpolation_matrix(self, components=None, weights=None, stretch=None # Compute scaled indices stretch_flat = stretch.reshape(1, self.n_signals * self.n_components) ** -1 - stretch_tiled = np.tile(stretch_flat, (self.signal_length, 1)) # Compute `fractional_indices` - fractional_indices = ( - np.tile(np.arange(self.signal_length)[:, None], (1, self.n_signals * self.n_components)) - * stretch_tiled - ) + fractional_indices = np.arange(self.signal_length)[:, None] * stretch_flat # Weighting matrix weights_flat = weights.reshape(1, self.n_signals * self.n_components) - weights_tiled = np.tile(weights_flat, (self.signal_length, 1)) # Bias for indexing into reshaped components # TODO break this up or describe what it does better @@ -439,17 +434,17 @@ def apply_interpolation_matrix(self, components=None, weights=None, stretch=None unweighted_stretched_comps = ( comp_values_1 * (1 - fractional_floor_indices) + comp_values_2 * fractional_floor_indices ) - stretched_components = unweighted_stretched_comps * weights_tiled # Apply weighting + stretched_components = unweighted_stretched_comps * weights_flat # Apply weighting # Compute first derivative - di = -fractional_indices * stretch_tiled + di = -fractional_indices * stretch_flat d_comps_unweighted = comp_values_1 * (-di) + comp_values_2 * di - d_stretched_components = d_comps_unweighted * weights_tiled + d_stretched_components = d_comps_unweighted * weights_flat # Compute second derivative - ddi = -di * stretch_tiled * 2 + ddi = -di * stretch_flat * 2 dd_comps_unweighted = comp_values_1 * (-ddi) + comp_values_2 * ddi - dd_stretched_components = dd_comps_unweighted * weights_tiled + dd_stretched_components = dd_comps_unweighted * weights_flat return stretched_components, d_stretched_components, dd_stretched_components From 6c1f8ec81f7559eaebb3c8e2b33a5569edb0dce4 Mon Sep 17 00:00:00 2001 From: John Halloran Date: Sun, 24 Aug 2025 21:10:46 -0700 Subject: [PATCH 3/5] refactor: flatten from a single buffer in apply_interpolation_matrix() --- src/diffpy/snmf/snmf_class.py | 16 +++++++++++----- 1 file changed, 11 insertions(+), 5 deletions(-) diff --git a/src/diffpy/snmf/snmf_class.py b/src/diffpy/snmf/snmf_class.py index c7c0516..2ef40bf 100644 --- a/src/diffpy/snmf/snmf_class.py +++ b/src/diffpy/snmf/snmf_class.py @@ -421,12 +421,18 @@ def apply_interpolation_matrix(self, components=None, weights=None, stretch=None offset_indices_1 = floor_indices_1 + bias offset_indices_2 = floor_indices_2 + bias - # Extract values - # Note: this "-1" corrects an off-by-one error that may have originated in an earlier line - comp_values_1 = components_bounded.flatten(order="F")[(offset_indices_1 - 1).ravel(order="F")].reshape( + # Flatten components once (Fortran order, column-major) + components_bounded_flat = components_bounded.ravel(order="F") + + # Pre-compute flattened indices + offset_indices_1_flat = (offset_indices_1 - 1).ravel(order="F") + offset_indices_2_flat = (offset_indices_2 - 1).ravel(order="F") + + # Extract values using pre-flattened arrays + comp_values_1 = components_bounded_flat[offset_indices_1_flat].reshape( self.signal_length, self.n_components * self.n_signals, order="F" - ) # order = F uses FORTRAN, column major order - comp_values_2 = components_bounded.flatten(order="F")[(offset_indices_2 - 1).ravel(order="F")].reshape( + ) + comp_values_2 = components_bounded_flat[offset_indices_2_flat].reshape( self.signal_length, self.n_components * self.n_signals, order="F" ) From 52ab77c9326c808cd94a6eba33cea2c7673a7cde Mon Sep 17 00:00:00 2001 From: John Halloran Date: Sun, 24 Aug 2025 22:32:36 -0700 Subject: [PATCH 4/5] refactor: drastically simplify indexing in apply_interpolation_matrix() and remove legacy MATLAB terminology --- src/diffpy/snmf/snmf_class.py | 141 +++++++++++++++++----------------- 1 file changed, 70 insertions(+), 71 deletions(-) diff --git a/src/diffpy/snmf/snmf_class.py b/src/diffpy/snmf/snmf_class.py index 2ef40bf..ee83e06 100644 --- a/src/diffpy/snmf/snmf_class.py +++ b/src/diffpy/snmf/snmf_class.py @@ -372,11 +372,31 @@ def get_objective_function(self, residuals=None, stretch=None): def apply_interpolation_matrix(self, components=None, weights=None, stretch=None): """ - Applies an interpolation-based transformation to the 'components' using `stretch`, - weighted by `weights`. Optionally computes first (`d_stretched_components`) and - second (`dd_stretched_components`) derivatives. + Interpolates each component along its sample axis according to per-(component, signal) + stretch factors, then applies per-(component, signal) weights. Also computes the + first and second derivatives with respect to stretch. Left and right, respectively, + refer to the sample prior to and subsequent to the interpolated sample's position. + + Inputs + ------ + components : array, shape (signal_len, n_components) + Each column is a component with signal_len samples. + weights : array, shape (n_components, n_signals) + Per-(component, signal) weights. + stretch : array, shape (n_components, n_signals) + Per-(component, signal) stretch factors. + + Outputs + ------- + stretched_components : array, shape (signal_len, n_components * n_signals) + Interpolated and weighted components. + d_stretched_components : array, shape (signal_len, n_components * n_signals) + First derivatives with respect to stretch. + dd_stretched_components : array, shape (signal_len, n_components * n_signals) + Second derivatives with respect to stretch. """ + # --- Defaults --- if components is None: components = self.components_ if weights is None: @@ -384,76 +404,55 @@ def apply_interpolation_matrix(self, components=None, weights=None, stretch=None if stretch is None: stretch = self.stretch_ - eps = 1e-8 # guard against divide by zero/NaN stretches - stretch = np.maximum(stretch, eps) - - # Compute scaled indices - stretch_flat = stretch.reshape(1, self.n_signals * self.n_components) ** -1 - - # Compute `fractional_indices` - fractional_indices = np.arange(self.signal_length)[:, None] * stretch_flat - - # Weighting matrix - weights_flat = weights.reshape(1, self.n_signals * self.n_components) - - # Bias for indexing into reshaped components - # TODO break this up or describe what it does better - bias = np.kron( - np.arange(self.n_components) * (self.signal_length + 1), - np.ones((self.signal_length, self.n_signals), dtype=int), - ).reshape(self.signal_length, self.n_components * self.n_signals) - - # Handle boundary conditions for interpolation - components_bounded = np.vstack( - [components, components[-1, :]] - ) # Duplicate last row (like MATLAB, not sure why) - - # Compute floor indices - floor_indices = np.floor(fractional_indices).astype(int) - - floor_indices_1 = np.minimum(floor_indices + 1, self.signal_length) - floor_indices_2 = np.minimum(floor_indices_1 + 1, self.signal_length) - - # Compute fractional part - fractional_floor_indices = fractional_indices - floor_indices - - # Compute offset indices - offset_indices_1 = floor_indices_1 + bias - offset_indices_2 = floor_indices_2 + bias - - # Flatten components once (Fortran order, column-major) - components_bounded_flat = components_bounded.ravel(order="F") - - # Pre-compute flattened indices - offset_indices_1_flat = (offset_indices_1 - 1).ravel(order="F") - offset_indices_2_flat = (offset_indices_2 - 1).ravel(order="F") - - # Extract values using pre-flattened arrays - comp_values_1 = components_bounded_flat[offset_indices_1_flat].reshape( - self.signal_length, self.n_components * self.n_signals, order="F" - ) - comp_values_2 = components_bounded_flat[offset_indices_2_flat].reshape( - self.signal_length, self.n_components * self.n_signals, order="F" + # Dimensions + signal_len = components.shape[0] # number of samples + n_components = components.shape[1] # number of components + n_signals = weights.shape[1] # number of signals + + # Guard stretches + eps = 1e-8 + stretch = np.clip(stretch, eps, None) + stretch_inv = 1.0 / stretch + + # Apply stretching to the original sample indices, represented as a "time-stretch" + t = np.arange(signal_len, dtype=float)[:, None, None] * stretch_inv[None, :, :] + # has shape (signal_len, n_components, n_signals) + + # For each stretched coordinate, find its prior integer (original) index and their difference + i0 = np.floor(t).astype(np.int64) # prior original index + alpha = t - i0.astype(float) # fractional distance between left/right + + # Clip indices to valid range (0, signal_len - 1) to maintain original size + max_idx = signal_len - 1 + i0 = np.clip(i0, 0, max_idx) + i1 = np.clip(i0 + 1, 0, max_idx) + + # Gather sample values + comps_3d = components[:, :, None] # expand components by a dimension for broadcasting across n_signals + c0 = np.take_along_axis(comps_3d, i0, axis=0) # left sample values + c1 = np.take_along_axis(comps_3d, i1, axis=0) # right sample values + + # Linear interpolation to determine stretched sample values + interp = c0 * (1.0 - alpha) + c1 * alpha + interp_weighted = interp * weights[None, :, :] + + # Derivatives + di = -t * stretch_inv[None, :, :] # first-derivative coefficient + ddi = -di * stretch_inv[None, :, :] * 2.0 # second-derivative coefficient + + d_unweighted = c0 * (-di) + c1 * di + dd_unweighted = c0 * (-ddi) + c1 * ddi + + d_weighted = d_unweighted * weights[None, :, :] + dd_weighted = dd_unweighted * weights[None, :, :] + + # Flatten back to expected shape (signal_len, n_components * n_signals) + return ( + interp_weighted.reshape(signal_len, n_components * n_signals), + d_weighted.reshape(signal_len, n_components * n_signals), + dd_weighted.reshape(signal_len, n_components * n_signals), ) - # Interpolation - unweighted_stretched_comps = ( - comp_values_1 * (1 - fractional_floor_indices) + comp_values_2 * fractional_floor_indices - ) - stretched_components = unweighted_stretched_comps * weights_flat # Apply weighting - - # Compute first derivative - di = -fractional_indices * stretch_flat - d_comps_unweighted = comp_values_1 * (-di) + comp_values_2 * di - d_stretched_components = d_comps_unweighted * weights_flat - - # Compute second derivative - ddi = -di * stretch_flat * 2 - dd_comps_unweighted = comp_values_1 * (-ddi) + comp_values_2 * ddi - dd_stretched_components = dd_comps_unweighted * weights_flat - - return stretched_components, d_stretched_components, dd_stretched_components - def apply_transformation_matrix(self, stretch=None, weights=None, residuals=None): """ Computes the transformation matrix `stretch_transformed` for residuals, From cdb3269eec6d524e2c90f85486661c67ad9b8554 Mon Sep 17 00:00:00 2001 From: John Halloran Date: Sun, 24 Aug 2025 22:39:25 -0700 Subject: [PATCH 5/5] style: rename apply_interpolation_matrix() to compute_stretched_components() --- src/diffpy/snmf/snmf_class.py | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/src/diffpy/snmf/snmf_class.py b/src/diffpy/snmf/snmf_class.py index ee83e06..658eeeb 100644 --- a/src/diffpy/snmf/snmf_class.py +++ b/src/diffpy/snmf/snmf_class.py @@ -370,7 +370,7 @@ def get_objective_function(self, residuals=None, stretch=None): function = residual_term + regularization_term + sparsity_term return function - def apply_interpolation_matrix(self, components=None, weights=None, stretch=None): + def compute_stretched_components(self, components=None, weights=None, stretch=None): """ Interpolates each component along its sample axis according to per-(component, signal) stretch factors, then applies per-(component, signal) weights. Also computes the @@ -563,7 +563,7 @@ def update_components(self): Updates `components` using gradient-based optimization with adaptive step size. """ # Compute stretched components using the interpolation function - stretched_components, _, _ = self.apply_interpolation_matrix() # Discard the derivatives + stretched_components, _, _ = self.compute_stretched_components() # Discard the derivatives # Compute reshaped_stretched_components and component_residuals intermediate_reshaped = stretched_components.flatten(order="F").reshape( (self.signal_length * self.n_signals, self.n_components), order="F" @@ -651,7 +651,9 @@ def regularize_function(self, stretch=None): if stretch is None: stretch = self.stretch_ - stretched_components, d_stretch_comps, dd_stretch_comps = self.apply_interpolation_matrix(stretch=stretch) + stretched_components, d_stretch_comps, dd_stretch_comps = self.compute_stretched_components( + stretch=stretch + ) intermediate = stretched_components.flatten(order="F").reshape( (self.signal_length * self.n_signals, self.n_components), order="F" ) @@ -754,8 +756,8 @@ def reconstruct_matrix(components, weights, stretch): """ signal_len = components.shape[0] - n_signals = weights.shape[1] n_components = components.shape[1] + n_signals = weights.shape[1] reconstructed_matrix = np.zeros((signal_len, n_signals)) sample_indices = np.arange(signal_len)