refactor: change get_objective_function into a static method and getter (#174)

* refactor: compute objective function in a static method and retrieve via getter

* refactor: change get_objective_function into a static method and getter

* fix: update location for _tests-on-pr.yml

* fix: use expected name for requirements/tests.txt

* fix: add cvxpy to requirements/conda.txt

* fix: add matplotlib to requirements/conda.txt

---------

Co-authored-by: John Halloran <jhalloran@oxy.edu>

Author: john-halloran

.github/workflows/tests-on-pr.yml

@@ -1,15 +1,12 @@
 name: Tests on PR
 
 on:
-  push:
-    branches:
-      - main
   pull_request:
   workflow_dispatch:
 
 jobs:
-  validate:
-    uses: Billingegroup/release-scripts/.github/workflows/_tests-on-pr.yml@v0
+  tests-on-pr:
+    uses: scikit-package/release-scripts/.github/workflows/_tests-on-pr.yml@v0
     with:
       project: diffpy.snmf
       c_extension: false

news/declass-obj.rst

@@ -0,0 +1,23 @@
+**Added:**
+
+* Implement tests for ``compute_objective_function()``
+
+**Changed:**
+
+* Refactor ``get_objective_function()`` into a static method and getter
+
+**Deprecated:**
+
+* <news item>
+
+**Removed:**
+
+* <news item>
+
+**Fixed:**
+
+* <news item>
+
+**Security:**
+
+* <news item>

requirements/conda.txt

@@ -1,4 +1,6 @@
 numpy
+matplotlib
 scipy
+cvxpy
 diffpy.utils
 numdifftools

requirements/test.txt renamed to requirements/tests.txt

@@ -359,16 +359,30 @@ def get_residual_matrix(self, components=None, weights=None, stretch=None):
         return residuals
 
     def get_objective_function(self, residuals=None, stretch=None):
-        if residuals is None:
-            residuals = self.residuals
-        if stretch is None:
-            stretch = self.stretch_
-        residual_term = 0.5 * np.linalg.norm(residuals, "fro") ** 2
-        regularization_term = 0.5 * self.rho * np.linalg.norm(self._spline_smooth_operator @ stretch.T, "fro") ** 2
-        sparsity_term = self.eta * np.sum(np.sqrt(self.components_))  # Square root penalty
-        # Final objective function value
-        function = residual_term + regularization_term + sparsity_term
-        return function
+        """
+        Return the objective value, passing stored attributes or overrides
+        to _compute_objective_function().
+
+        Parameters
+        ----------
+        residuals : ndarray, optional
+            Residual matrix to use instead of self.residuals.
+        stretch : ndarray, optional
+            Stretch matrix to use instead of self.stretch_.
+
+        Returns
+        -------
+        float
+            Current objective function value.
+        """
+        return SNMFOptimizer._compute_objective_function(
+            components=self.components_,
+            residuals=self.residuals if residuals is None else residuals,
+            stretch=self.stretch_ if stretch is None else stretch,
+            rho=self.rho,
+            eta=self.eta,
+            spline_smooth_operator=self._spline_smooth_operator,
+        )
 
     def compute_stretched_components(self, components=None, weights=None, stretch=None):
         """
@@ -702,6 +716,59 @@ def objective(stretch_vec):
         # Update stretch with the optimized values
         self.stretch_ = result.x.reshape(self.stretch_.shape)
 
+    @staticmethod
+    def _compute_objective_function(components, residuals, stretch, rho, eta, spline_smooth_operator):
+        r"""
+        Computes the objective function used in stretched non-negative matrix factorization.
+
+        Parameters
+        ----------
+        components : ndarray
+            Non-negative matrix of component signals :math:`X`.
+        residuals : ndarray
+            Difference between reconstructed and observed data.
+        stretch : ndarray
+            Stretching factors :math:`A` applied to each component across samples.
+        rho : float
+            Regularization parameter enforcing smooth variation in :math:`A`.
+        eta : float
+            Sparsity-promoting regularization parameter applied to :math:`X`.
+        spline_smooth_operator : ndarray
+            Linear operator :math:`L` penalizing non-smooth changes in :math:`A`.
+
+        Returns
+        -------
+        float
+            Value of the stretched-NMF objective function.
+
+        Notes
+        -----
+        The stretched-NMF objective function :math:`J` is
+
+        .. math::
+
+           J(X, Y, A) =
+              \tfrac{1}{2} \lVert Z - Y\,S(A)X \rVert_F^2
+              + \tfrac{\rho}{2} \lVert L A \rVert_F^2
+              + \eta \sum_{i,j} \sqrt{X_{ij}} \,,
+
+        where :math:`Z` is the data matrix, :math:`Y` contains the non-negative
+        weights, :math:`S(A)` denotes the spline-interpolated stretching operator,
+        and :math:`\lVert \cdot \rVert_F` is the Frobenius norm.
+
+        Special cases
+        -------------
+        - :math:`\rho = 0` — no smoothness regularization on stretching factors.
+        - :math:`\eta = 0` — no sparsity promotion on components.
+        - :math:`\rho = \eta = 0` — reduces to the classical NMF least-squares
+          objective :math:`\tfrac{1}{2} \lVert Z - YX \rVert_F^2`.
+
+        """
+        residual_term = 0.5 * np.linalg.norm(residuals, "fro") ** 2
+        regularization_term = 0.5 * rho * np.linalg.norm(spline_smooth_operator @ stretch.T, "fro") ** 2
+        sparsity_term = eta * np.sum(np.sqrt(components))
+        return residual_term + regularization_term + sparsity_term
+
 
 def cubic_largest_real_root(p, q):
     """

src/diffpy/snmf/snmf_class.py

@@ -45,3 +45,75 @@ def test_final_objective_below_threshold(inputs):
     # Basic sanity check and the actual assertion
     assert np.isfinite(model.objective_function)
     assert model.objective_function < 5e6
+
+
+@pytest.mark.parametrize(
+    "inputs, expected",
+    # inputs tuple: (components, residuals, stretch, rho, eta, spline smoothness operator)
+    [
+        # Case 0: No smoothness or sparsity penalty, reduces to standard NMF objective
+        # residual Frobenius norm^2 = 3^2 + 4^2 = 25 -> 0.5 * 25 = 12.5
+        (
+            (
+                np.array([[0.0, 0.0], [3.0, 4.0]]),
+                np.array([[0.0, 0.0], [3.0, 4.0]]),
+                np.ones((2, 2)),
+                0.0,
+                0.0,
+                np.zeros((2, 2)),
+            ),
+            12.5,
+        ),
+        # Case 1: rho = 0, sparsity penalty only
+        # sqrt components sum = 1 + 2 + 3 + 4 = 10 -> eta * 10 = 5
+        # residual term remains 12.5 -> total = 17.5
+        (
+            (
+                np.array([[1.0, 4.0], [9.0, 16.0]]),
+                np.array([[3.0, 4.0], [0.0, 0.0]]),
+                np.ones((2, 2)),
+                0.0,
+                0.5,
+                np.zeros((2, 2)),
+            ),
+            17.5,
+        ),
+        # Case 2: eta = 0, smoothness penalty only
+        # residual = 12.5, smoothing = 0.5 * 1 * 1 = 0.5 -> total = 13.0
+        (
+            (
+                np.array([[1.0, 2.0], [3.0, 4.0]]),
+                np.array([[3.0, 4.0], [0.0, 0.0]]),
+                np.array([[1.0, 2.0]]),
+                1.0,
+                0.0,
+                np.array([[1.0, -1.0]]),
+            ),
+            13.0,
+        ),
+        # Case 3: penalty for smoothness and sparsity
+        # residual = 2.5, sparsity = 1.5, smoothing = 9 -> total = 13.0
+        (
+            (
+                np.array([[1.0, 4.0]]),
+                np.array([[1.0, 2.0]]),
+                np.array([[1.0, 4.0]]),
+                2.0,
+                0.5,
+                np.array([[3.0, 0.0]]),
+            ),
+            13.0,
+        ),
+    ],
+)
+def test_compute_objective_function(inputs, expected):
+    components, residuals, stretch, rho, eta, operator = inputs
+    result = SNMFOptimizer._compute_objective_function(
+        components=components,
+        residuals=residuals,
+        stretch=stretch,
+        rho=rho,
+        eta=eta,
+        spline_smooth_operator=operator,
+    )
+    assert np.isclose(result, expected)