diff --git a/.github/workflows/tests-on-pr.yml b/.github/workflows/tests-on-pr.yml index fabc1066..3c428bee 100644 --- a/.github/workflows/tests-on-pr.yml +++ b/.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 diff --git a/news/declass-obj.rst b/news/declass-obj.rst new file mode 100644 index 00000000..5573f45c --- /dev/null +++ b/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:** + +* + +**Removed:** + +* + +**Fixed:** + +* + +**Security:** + +* diff --git a/requirements/conda.txt b/requirements/conda.txt index 7f992d10..60eeedae 100644 --- a/requirements/conda.txt +++ b/requirements/conda.txt @@ -1,4 +1,6 @@ numpy +matplotlib scipy +cvxpy diffpy.utils numdifftools diff --git a/requirements/test.txt b/requirements/tests.txt similarity index 100% rename from requirements/test.txt rename to requirements/tests.txt diff --git a/src/diffpy/snmf/snmf_class.py b/src/diffpy/snmf/snmf_class.py index 658eeebf..4999168f 100644 --- a/src/diffpy/snmf/snmf_class.py +++ b/src/diffpy/snmf/snmf_class.py @@ -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): """ diff --git a/tests/test_snmf_optimizer.py b/tests/test_snmf_optimizer.py index 42f6ae91..48c0a27f 100644 --- a/tests/test_snmf_optimizer.py +++ b/tests/test_snmf_optimizer.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)