From 38bbb00d37b260d0d259da14ff984d54e3fd49e6 Mon Sep 17 00:00:00 2001 From: John Halloran Date: Tue, 12 Aug 2025 14:36:22 -0700 Subject: [PATCH 1/5] refactor: get residual matrix without a helper --- src/diffpy/snmf/snmf_class.py | 72 +++++++++++++++++++++++++++++++++-- 1 file changed, 68 insertions(+), 4 deletions(-) diff --git a/src/diffpy/snmf/snmf_class.py b/src/diffpy/snmf/snmf_class.py index 489aabd..17cabea 100644 --- a/src/diffpy/snmf/snmf_class.py +++ b/src/diffpy/snmf/snmf_class.py @@ -302,6 +302,51 @@ 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. + + Parameters + ---------- + components : (signal_len, n_components) array, optional + weights : (n_components, n_signals) array, optional + stretch : (n_components, n_signals) array, optional + + Returns + ------- + residuals : (signal_len, n_signals) array + """ + + if components is None: + components = self.components + if weights is None: + weights = self.weights + 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 + ) + + return residuals + + def old_get_residual_matrix(self, components=None, weights=None, stretch=None): # Initialize residual matrix as negative of source_matrix if components is None: components = self.components @@ -310,10 +355,29 @@ def get_residual_matrix(self, components=None, weights=None, stretch=None): if stretch is None: stretch = self.stretch residuals = -self.source_matrix.copy() - # Compute transformed components for all (k, m) pairs - for k in range(weights.shape[0]): # K - stretched_components, _, _ = apply_interpolation(stretch[k, :], components[:, k]) # Only use Ax - residuals += weights[k, :] * stretched_components # Element-wise scaling and sum + + # Discrete sample positions along the component axis + sample_indices = np.arange(components.shape[0]) # (N,) + + for comp in range(components.shape[1]): # loop over components + component_profile = components[:, comp] # (N,) + stretch_factors = stretch[comp, :] # (M,) + + # Compute scaled/fractional positions along component_profile + fractional_positions = sample_indices[:, None] / stretch_factors[None, :] + + # Interpolate component_profile at fractional positions, clamp to ends + interpolated_component = np.interp( + fractional_positions, + sample_indices, + component_profile, + left=component_profile[0], + right=component_profile[-1], + ) + + # Accumulate weighted contribution into residuals + residuals += interpolated_component * weights[comp, None, :] # (M,) broadcast + return residuals def get_objective_function(self, residuals=None, stretch=None): From da3a93894412b4b00530bcfb5b0cde9da7f64452 Mon Sep 17 00:00:00 2001 From: John Halloran Date: Tue, 12 Aug 2025 22:41:03 -0700 Subject: [PATCH 2/5] perf: remove unused derivatives from apply_interpolation --- src/diffpy/snmf/snmf_class.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/src/diffpy/snmf/snmf_class.py b/src/diffpy/snmf/snmf_class.py index 17cabea..0b79556 100644 --- a/src/diffpy/snmf/snmf_class.py +++ b/src/diffpy/snmf/snmf_class.py @@ -653,7 +653,7 @@ def update_weights(self): # Populate t using apply_interpolation for k in range(self.n_components): - t[:, k] = apply_interpolation(self.stretch[k, m], self.components[:, k])[0].squeeze() + t[:, k] = apply_interpolation(self.stretch[k, m], self.components[:, k]).squeeze() # Solve quadratic problem for y y = self.solve_quadratic_program(t=t, m=m) @@ -780,6 +780,7 @@ def apply_interpolation(a, x): intr_x_tail = np.full((x_len - len(idx_int), interpolated_x.shape[1]), interpolated_x[-1, :]) interpolated_x = np.vstack([interpolated_x, intr_x_tail]) + """ # Compute first derivative (d_intr_x) di = -idx_frac / a d_intr_x = x[idx_int] * (-di) + x[np.minimum(idx_int + 1, x_len - 1)] * di @@ -789,5 +790,6 @@ def apply_interpolation(a, x): ddi = -di / a + idx_frac * a**-2 dd_intr_x = x[idx_int] * (-ddi) + x[np.minimum(idx_int + 1, x_len - 1)] * ddi dd_intr_x = np.vstack([dd_intr_x, np.zeros((x_len - len(idx_int), dd_intr_x.shape[1]))]) + """ - return interpolated_x, d_intr_x, dd_intr_x + return interpolated_x # , d_intr_x, dd_intr_x From c18d8682621e7011a78f6278ea1154f215d0271c Mon Sep 17 00:00:00 2001 From: John Halloran Date: Tue, 12 Aug 2025 22:45:46 -0700 Subject: [PATCH 3/5] chore: remove old residual matrix and reference to derivatives --- src/diffpy/snmf/snmf_class.py | 49 +---------------------------------- 1 file changed, 1 insertion(+), 48 deletions(-) diff --git a/src/diffpy/snmf/snmf_class.py b/src/diffpy/snmf/snmf_class.py index 0b79556..424585e 100644 --- a/src/diffpy/snmf/snmf_class.py +++ b/src/diffpy/snmf/snmf_class.py @@ -346,40 +346,6 @@ def get_residual_matrix(self, components=None, weights=None, stretch=None): return residuals - def old_get_residual_matrix(self, components=None, weights=None, stretch=None): - # Initialize residual matrix as negative of source_matrix - if components is None: - components = self.components - if weights is None: - weights = self.weights - if stretch is None: - stretch = self.stretch - residuals = -self.source_matrix.copy() - - # Discrete sample positions along the component axis - sample_indices = np.arange(components.shape[0]) # (N,) - - for comp in range(components.shape[1]): # loop over components - component_profile = components[:, comp] # (N,) - stretch_factors = stretch[comp, :] # (M,) - - # Compute scaled/fractional positions along component_profile - fractional_positions = sample_indices[:, None] / stretch_factors[None, :] - - # Interpolate component_profile at fractional positions, clamp to ends - interpolated_component = np.interp( - fractional_positions, - sample_indices, - component_profile, - left=component_profile[0], - right=component_profile[-1], - ) - - # Accumulate weighted contribution into residuals - residuals += interpolated_component * weights[comp, None, :] # (M,) broadcast - - return residuals - def get_objective_function(self, residuals=None, stretch=None): if residuals is None: residuals = self.residuals @@ -751,7 +717,6 @@ def cubic_largest_real_root(p, q): def apply_interpolation(a, x): """ Applies an interpolation-based transformation to `x` based on scaling `a`. - Also computes first (`d_intr_x`) and second (`dd_intr_x`) derivatives. """ x_len = len(x) @@ -780,16 +745,4 @@ def apply_interpolation(a, x): intr_x_tail = np.full((x_len - len(idx_int), interpolated_x.shape[1]), interpolated_x[-1, :]) interpolated_x = np.vstack([interpolated_x, intr_x_tail]) - """ - # Compute first derivative (d_intr_x) - di = -idx_frac / a - d_intr_x = x[idx_int] * (-di) + x[np.minimum(idx_int + 1, x_len - 1)] * di - d_intr_x = np.vstack([d_intr_x, np.zeros((x_len - len(idx_int), d_intr_x.shape[1]))]) - - # Compute second derivative (dd_intr_x) - ddi = -di / a + idx_frac * a**-2 - dd_intr_x = x[idx_int] * (-ddi) + x[np.minimum(idx_int + 1, x_len - 1)] * ddi - dd_intr_x = np.vstack([dd_intr_x, np.zeros((x_len - len(idx_int), dd_intr_x.shape[1]))]) - """ - - return interpolated_x # , d_intr_x, dd_intr_x + return interpolated_x From 640dacfce6d09079b4a35dcdc51b052e6b9d5fb0 Mon Sep 17 00:00:00 2001 From: John Halloran Date: Tue, 12 Aug 2025 23:54:58 -0700 Subject: [PATCH 4/5] refactor: replace remaining apply_interpolation with np.interp --- src/diffpy/snmf/snmf_class.py | 69 ++++++++++------------------------- 1 file changed, 20 insertions(+), 49 deletions(-) diff --git a/src/diffpy/snmf/snmf_class.py b/src/diffpy/snmf/snmf_class.py index 424585e..3d408ec 100644 --- a/src/diffpy/snmf/snmf_class.py +++ b/src/diffpy/snmf/snmf_class.py @@ -608,24 +608,29 @@ def update_components(self): def update_weights(self): """ - Updates weights using matrix operations, solving a quadratic program to do so. + Updates weights by building the stretched component matrix `stretched_comps` with np.interp + and solving a quadratic program for each signal. """ - signal_length = self.signal_length - n_signals = self.n_signals - - for m in range(n_signals): - t = np.zeros((signal_length, self.n_components)) - - # Populate t using apply_interpolation - for k in range(self.n_components): - t[:, k] = apply_interpolation(self.stretch[k, m], self.components[:, k]).squeeze() - - # Solve quadratic problem for y - y = self.solve_quadratic_program(t=t, m=m) + sample_indices = np.arange(self.signal_length) + for signal in range(self.n_signals): + # Stretch factors for this signal across components: + this_stretch = self.stretch[:, signal] + # Build stretched_comps[:, k] by interpolating component at frac. pos. index / this_stretch[comp] + stretched_comps = np.empty((self.signal_length, self.n_components), dtype=self.components.dtype) + for comp in range(self.n_components): + pos = sample_indices / this_stretch[comp] + stretched_comps[:, comp] = np.interp( + pos, + sample_indices, + self.components[:, comp], + left=self.components[0, comp], + right=self.components[-1, comp], + ) - # Update Y - self.weights[:, m] = y + # Solve quadratic problem for a given signal and update its weight + new_weight = self.solve_quadratic_program(t=stretched_comps, m=signal) + self.weights[:, signal] = new_weight def regularize_function(self, stretch=None): if stretch is None: @@ -712,37 +717,3 @@ 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 apply_interpolation(a, x): - """ - Applies an interpolation-based transformation to `x` based on scaling `a`. - """ - x_len = len(x) - - # Ensure `a` is an array and reshape for broadcasting - a = np.atleast_1d(np.asarray(a)) # Ensures a is at least 1D - - # Compute fractional indices, broadcasting over `a` - fractional_indices = np.arange(x_len)[:, None] / a # Shape (N, M) - - integer_indices = np.floor(fractional_indices).astype(int) # Integer part (still (N, M)) - valid_mask = integer_indices < (x_len - 1) # Ensure indices are within bounds - - # Apply valid_mask to keep correct indices - idx_int = np.where(valid_mask, integer_indices, x_len - 2) # Prevent out-of-bounds indexing (previously "I") - idx_frac = np.where(valid_mask, fractional_indices, integer_indices) # Keep aligned (previously "i") - - # Ensure x is a 1D array - x = np.asarray(x).ravel() - - # Compute interpolated_x (linear interpolation) - interpolated_x = x[idx_int] * (1 - idx_frac + idx_int) + x[np.minimum(idx_int + 1, x_len - 1)] * ( - idx_frac - idx_int - ) - - # Fill the tail with the last valid value - intr_x_tail = np.full((x_len - len(idx_int), interpolated_x.shape[1]), interpolated_x[-1, :]) - interpolated_x = np.vstack([interpolated_x, intr_x_tail]) - - return interpolated_x From f81ac60f78d962e8a2a857988524dd80204bc0d1 Mon Sep 17 00:00:00 2001 From: John Halloran Date: Wed, 13 Aug 2025 00:15:05 -0700 Subject: [PATCH 5/5] style: remove references to old variable names --- src/diffpy/snmf/snmf_class.py | 22 +++++++++++----------- 1 file changed, 11 insertions(+), 11 deletions(-) diff --git a/src/diffpy/snmf/snmf_class.py b/src/diffpy/snmf/snmf_class.py index 3d408ec..06385ef 100644 --- a/src/diffpy/snmf/snmf_class.py +++ b/src/diffpy/snmf/snmf_class.py @@ -636,19 +636,19 @@ def regularize_function(self, stretch=None): if stretch is None: stretch = self.stretch - K = self.n_components - M = self.n_signals - N = self.signal_length - stretched_components, d_stretch_comps, dd_stretch_comps = self.apply_interpolation_matrix(stretch=stretch) - intermediate = stretched_components.flatten(order="F").reshape((N * M, K), order="F") - residuals = intermediate.sum(axis=1).reshape((N, M), order="F") - self.source_matrix + intermediate = stretched_components.flatten(order="F").reshape( + (self.signal_length * self.n_signals, self.n_components), order="F" + ) + residuals = ( + intermediate.sum(axis=1).reshape((self.signal_length, self.n_signals), order="F") - self.source_matrix + ) fun = self.get_objective_function(residuals, stretch) - tiled_res = np.tile(residuals, (1, K)) + tiled_res = np.tile(residuals, (1, self.n_components)) grad_flat = np.sum(d_stretch_comps * tiled_res, axis=0) - gra = grad_flat.reshape((M, K), order="F").T + gra = grad_flat.reshape((self.n_signals, self.n_components), order="F").T gra += self.rho * stretch @ (self._spline_smooth_operator.T @ self._spline_smooth_operator) # Hessian would go here @@ -657,10 +657,10 @@ def regularize_function(self, stretch=None): def update_stretch(self): """ - Updates matrix A using constrained optimization (equivalent to fmincon in MATLAB). + Updates stretching matrix using constrained optimization (equivalent to fmincon in MATLAB). """ - # Flatten A for compatibility with the optimizer (since SciPy expects 1D input) + # Flatten stretch for compatibility with the optimizer (since SciPy expects 1D input) stretch_flat_initial = self.stretch.flatten() # Define the optimization function @@ -682,7 +682,7 @@ def objective(stretch_vec): bounds=bounds, ) - # Update A with the optimized values + # Update stretch with the optimized values self.stretch = result.x.reshape(self.stretch.shape)