From d8d4e1122c0b9ceaf44e79401f38b66350775b8a Mon Sep 17 00:00:00 2001 From: John Halloran Date: Fri, 6 Jun 2025 12:45:48 -0700 Subject: [PATCH 01/11] feat: Add random state feature. --- src/diffpy/snmf/snmf_class.py | 25 ++++++++++++++++++------- 1 file changed, 18 insertions(+), 7 deletions(-) diff --git a/src/diffpy/snmf/snmf_class.py b/src/diffpy/snmf/snmf_class.py index 9d58e5f2..150f6541 100644 --- a/src/diffpy/snmf/snmf_class.py +++ b/src/diffpy/snmf/snmf_class.py @@ -4,8 +4,20 @@ class SNMFOptimizer: - def __init__(self, MM, Y0=None, X0=None, A=None, rho=1e12, eta=610, max_iter=500, tol=5e-7, components=None): - print("Initializing SNMF Optimizer") + def __init__( + self, + MM, + Y0=None, + X0=None, + A=None, + rho=1e12, + eta=610, + max_iter=500, + tol=5e-7, + components=None, + random_state=None, + ): + self.MM = MM self.X0 = X0 self.Y0 = Y0 @@ -15,23 +27,22 @@ def __init__(self, MM, Y0=None, X0=None, A=None, rho=1e12, eta=610, max_iter=500 # Capture matrix dimensions self.N, self.M = MM.shape self.num_updates = 0 + self.rng = np.random.default_rng(random_state) if Y0 is None: if components is None: raise ValueError("Must provide either Y0 or a number of components.") else: self.K = components - self.Y0 = np.random.beta(a=2.5, b=1.5, size=(self.K, self.M)) # This is untested + self.Y0 = self.rng.beta(a=2.5, b=1.5, size=(self.K, self.M)) else: self.K = Y0.shape[0] - # Initialize A, X0 if not provided if self.A is None: - self.A = np.ones((self.K, self.M)) + np.random.randn(self.K, self.M) * 1e-3 # Small perturbation + self.A = np.ones((self.K, self.M)) + self.rng.normal(0, 1e-3, size=(self.K, self.M)) if self.X0 is None: - self.X0 = np.random.rand(self.N, self.K) # Ensures values in [0,1] + self.X0 = self.rng.random((self.N, self.K)) - # Initialize solution matrices to be iterated on self.X = np.maximum(0, self.X0) self.Y = np.maximum(0, self.Y0) From ae45726407c90ff5f7c0121bf9959b96f0ea0f6b Mon Sep 17 00:00:00 2001 From: John Halloran Date: Sat, 7 Jun 2025 23:21:52 -0700 Subject: [PATCH 02/11] Add class docstring --- src/diffpy/snmf/snmf_class.py | 58 +++++++++++++++++++++++++++++++++++ 1 file changed, 58 insertions(+) diff --git a/src/diffpy/snmf/snmf_class.py b/src/diffpy/snmf/snmf_class.py index 150f6541..a40c2512 100644 --- a/src/diffpy/snmf/snmf_class.py +++ b/src/diffpy/snmf/snmf_class.py @@ -17,6 +17,64 @@ def __init__( components=None, random_state=None, ): + """Run sNMF based on an ndarray, parameters, and either a number + of components or a set of initial guess matrices. + + Currently instantiating the SNMFOptimizer class runs all the analysis + immediately. The results can then be accessed as instance attributes + of the class (X, Y, and A). Eventually, this will be changed such + that __init__ only prepares for the optimization, which will can then + be done using fit_transform. + + Parameters + ---------- + MM: ndarray + A numpy array containing the data to be decomposed. Rows correspond + to different samples/angles, while columns correspond to different + conditions with different stretching. Currently, there is no option + to treat the first column (commonly containing 2theta angles, sample + index, etc) differently, so if present it must be stripped in advance. + Y0: ndarray + A numpy array containing initial guesses for the component weights + at each stretching condition, with number of rows equal to the assumed + number of components and number of columns equal to the number of + conditions (same number of columns as MM). Must be provided if components + is not provided. Will override components if both are provided. + X0: ndarray + A numpy array containing initial guesses for the intensities of each + component per row/sample/angle. Has rows equal to the rows of MM and + columns equal to n_components or the number of rows of Y0. + A: ndarray + A numpy array containing initial guesses for the stretching factor for + each component, at each condition. Has number of rows equal to n_components + or the number of rows of Y0, and columns equal to the number of conditions + (columns of MM). + rho: float + A stretching factor that influences the decomposition. Zero corresponds to + no stretching present. Relatively insensitive and typically adjusted in + powers of 10. + eta: float + A sparsity factor than influences the decomposition. Should be set to zero + for non sparse data such as PDF. Can be used to improve results for sparse + data such as XRD, but due to instability, should be used only after first + selecting the best value for rho. + max_iter: int + The maximum number of times to update each of A, X, and Y before stopping + the optimization. + tol: float + The minimum fractional improvement in the objective function to allow + without terminating the optimization. Note that a minimum of 20 updates + are run before this parameter is checked. + components: int + The number of components to attempt to extract from MM. Note that this will + be overridden by Y0 if that is provided, but must be provided if no Y0 is + provided. + random_state: int + Used to set a reproducible seed for the initial matrices used in the + optimization. Due to the non-convex nature of the problem, results may vary + even with the same initial guesses, so this does not make the program + deterministic. + """ self.MM = MM self.X0 = X0 From 3d7c8b6a7a9273b71aa49a2db82f5fb0a4feaeb1 Mon Sep 17 00:00:00 2001 From: John Halloran Date: Sun, 8 Jun 2025 01:13:02 -0700 Subject: [PATCH 03/11] components->n_components --- src/diffpy/snmf/main.py | 21 +-------------------- src/diffpy/snmf/snmf_class.py | 15 ++++++++------- 2 files changed, 9 insertions(+), 27 deletions(-) diff --git a/src/diffpy/snmf/main.py b/src/diffpy/snmf/main.py index 7afacf6b..a512f67e 100644 --- a/src/diffpy/snmf/main.py +++ b/src/diffpy/snmf/main.py @@ -7,27 +7,8 @@ Y0 = np.loadtxt("input/W0.txt", dtype=float) N, M = MM.shape -# Convert to DataFrames for display -# df_X = pd.DataFrame(X, columns=[f"Comp_{i+1}" for i in range(X.shape[1])]) -# df_Y = pd.DataFrame(Y, columns=[f"Sample_{i+1}" for i in range(Y.shape[1])]) -# df_MM = pd.DataFrame(MM, columns=[f"Sample_{i+1}" for i in range(MM.shape[1])]) -# df_Y0 = pd.DataFrame(Y0, columns=[f"Sample_{i+1}" for i in range(Y0.shape[1])]) - -# Print the matrices -""" -print("Feature Matrix (X):\n", df_X, "\n") -print("Coefficient Matrix (Y):\n", df_Y, "\n") -print("Data Matrix (MM):\n", df_MM, "\n") -print("Initial Guess (Y0):\n", df_Y0, "\n") -""" - - -my_model = snmf_class.SNMFOptimizer(MM=MM, Y0=Y0, X0=X0, A=A0, components=2) +my_model = snmf_class.SNMFOptimizer(MM=MM, Y0=Y0, X0=X0, A=A0, n_components=2) print("Done") -# print(f"My final guess for X: {my_model.X}") -# print(f"My final guess for Y: {my_model.Y}") -# print(f"Compare to true X: {X_norm}") -# print(f"Compare to true Y: {Y_norm}") np.savetxt("my_norm_X.txt", my_model.X, fmt="%.6g", delimiter=" ") np.savetxt("my_norm_Y.txt", my_model.Y, fmt="%.6g", delimiter=" ") np.savetxt("my_norm_A.txt", my_model.A, fmt="%.6g", delimiter=" ") diff --git a/src/diffpy/snmf/snmf_class.py b/src/diffpy/snmf/snmf_class.py index a40c2512..da0fcc2d 100644 --- a/src/diffpy/snmf/snmf_class.py +++ b/src/diffpy/snmf/snmf_class.py @@ -14,7 +14,7 @@ def __init__( eta=610, max_iter=500, tol=5e-7, - components=None, + n_components=None, random_state=None, ): """Run sNMF based on an ndarray, parameters, and either a number @@ -38,8 +38,9 @@ def __init__( A numpy array containing initial guesses for the component weights at each stretching condition, with number of rows equal to the assumed number of components and number of columns equal to the number of - conditions (same number of columns as MM). Must be provided if components - is not provided. Will override components if both are provided. + conditions (same number of columns as MM). Must be provided if + n_components is not provided. Will override n_components if both are + provided. X0: ndarray A numpy array containing initial guesses for the intensities of each component per row/sample/angle. Has rows equal to the rows of MM and @@ -65,7 +66,7 @@ def __init__( The minimum fractional improvement in the objective function to allow without terminating the optimization. Note that a minimum of 20 updates are run before this parameter is checked. - components: int + n_components: int The number of components to attempt to extract from MM. Note that this will be overridden by Y0 if that is provided, but must be provided if no Y0 is provided. @@ -88,10 +89,10 @@ def __init__( self.rng = np.random.default_rng(random_state) if Y0 is None: - if components is None: - raise ValueError("Must provide either Y0 or a number of components.") + if n_components is None: + raise ValueError("Must provide either Y0 or n_components.") else: - self.K = components + self.K = n_components self.Y0 = self.rng.beta(a=2.5, b=1.5, size=(self.K, self.M)) else: self.K = Y0.shape[0] From a0483b4b855ebdca2765c2da493edd346101280c Mon Sep 17 00:00:00 2001 From: John Halloran Date: Mon, 9 Jun 2025 15:56:44 -0700 Subject: [PATCH 04/11] Updated docstring --- src/diffpy/snmf/snmf_class.py | 70 +++++++++++++++++------------------ 1 file changed, 33 insertions(+), 37 deletions(-) diff --git a/src/diffpy/snmf/snmf_class.py b/src/diffpy/snmf/snmf_class.py index da0fcc2d..6ab1fb40 100644 --- a/src/diffpy/snmf/snmf_class.py +++ b/src/diffpy/snmf/snmf_class.py @@ -4,6 +4,18 @@ class SNMFOptimizer: + """A self-contained implementation of the stretched NMF algorithm (sNMF), + including sparse stretched NMF. + + Instantiating the SNMFOptimizer class runs all the analysis immediately. + The results matrices can then be accessed as instance attributes + of the class (X, Y, and A). + + For more information on sNMF, please reference: + Gu, R., Rakita, Y., Lan, L. et al. Stretched non-negative matrix factorization. + npj Comput Mater 10, 193 (2024). https://doi.org/10.1038/s41524-024-01377-5 + """ + def __init__( self, MM, @@ -17,48 +29,33 @@ def __init__( n_components=None, random_state=None, ): - """Run sNMF based on an ndarray, parameters, and either a number - of components or a set of initial guess matrices. - - Currently instantiating the SNMFOptimizer class runs all the analysis - immediately. The results can then be accessed as instance attributes - of the class (X, Y, and A). Eventually, this will be changed such - that __init__ only prepares for the optimization, which will can then - be done using fit_transform. + """Initialize an instance of SNMF and run the optimization Parameters ---------- MM: ndarray - A numpy array containing the data to be decomposed. Rows correspond - to different samples/angles, while columns correspond to different - conditions with different stretching. Currently, there is no option - to treat the first column (commonly containing 2theta angles, sample - index, etc) differently, so if present it must be stripped in advance. + The array containing the data to be decomposed. Shape is (length_of_signal, + number_of_conditions). Y0: ndarray - A numpy array containing initial guesses for the component weights - at each stretching condition, with number of rows equal to the assumed - number of components and number of columns equal to the number of - conditions (same number of columns as MM). Must be provided if - n_components is not provided. Will override n_components if both are - provided. + The array containing initial guesses for the component weights + at each stretching condition. Shape is (number of components, number of + conditions) Must be provided if n_components is not provided. Will override + n_components if both are provided. X0: ndarray - A numpy array containing initial guesses for the intensities of each - component per row/sample/angle. Has rows equal to the rows of MM and - columns equal to n_components or the number of rows of Y0. + The array containing initial guesses for the intensities of each component per + row/sample/angle. Shape is (length_of_signal, number_of_components). A: ndarray - A numpy array containing initial guesses for the stretching factor for - each component, at each condition. Has number of rows equal to n_components - or the number of rows of Y0, and columns equal to the number of conditions - (columns of MM). + The array containing initial guesses for the stretching factor for each component, + at each condition. Shape is (number_of_components, number_of_conditions). rho: float - A stretching factor that influences the decomposition. Zero corresponds to - no stretching present. Relatively insensitive and typically adjusted in - powers of 10. + The float which sets a stretching factor that influences the decomposition. + Zero corresponds to no stretching present. Relatively insensitive and typically + adjusted in powers of 10. eta: float - A sparsity factor than influences the decomposition. Should be set to zero - for non sparse data such as PDF. Can be used to improve results for sparse - data such as XRD, but due to instability, should be used only after first - selecting the best value for rho. + The integer which sets a sparsity factor than influences the decomposition. + Should be set to zero for non sparse data such as PDF. Can be used to improve + results for sparse data such as XRD, but due to instability, should be used + only after first selecting the best value for rho. max_iter: int The maximum number of times to update each of A, X, and Y before stopping the optimization. @@ -71,10 +68,9 @@ def __init__( be overridden by Y0 if that is provided, but must be provided if no Y0 is provided. random_state: int - Used to set a reproducible seed for the initial matrices used in the - optimization. Due to the non-convex nature of the problem, results may vary - even with the same initial guesses, so this does not make the program - deterministic. + The integer which acts as a reproducible seed for the initial matrices used in + the optimization. Due to the non-convex nature of the problem, results may vary + even with the same initial guesses, so this does not make the program deterministic. """ self.MM = MM From d39cbe0e9c89ed017b2098aa2e445227fcc16d5f Mon Sep 17 00:00:00 2001 From: John Halloran Date: Thu, 12 Jun 2025 17:12:43 -0700 Subject: [PATCH 05/11] Shorten and reformat docstring --- src/diffpy/snmf/snmf_class.py | 62 ++++++++++++++++------------------- 1 file changed, 28 insertions(+), 34 deletions(-) diff --git a/src/diffpy/snmf/snmf_class.py b/src/diffpy/snmf/snmf_class.py index 6ab1fb40..3650abc0 100644 --- a/src/diffpy/snmf/snmf_class.py +++ b/src/diffpy/snmf/snmf_class.py @@ -4,8 +4,7 @@ class SNMFOptimizer: - """A self-contained implementation of the stretched NMF algorithm (sNMF), - including sparse stretched NMF. + """A implementation of stretched NMF (sNMF), including sparse stretched NMF. Instantiating the SNMFOptimizer class runs all the analysis immediately. The results matrices can then be accessed as instance attributes @@ -33,44 +32,39 @@ def __init__( Parameters ---------- - MM: ndarray - The array containing the data to be decomposed. Shape is (length_of_signal, - number_of_conditions). - Y0: ndarray - The array containing initial guesses for the component weights - at each stretching condition. Shape is (number of components, number of - conditions) Must be provided if n_components is not provided. Will override - n_components if both are provided. - X0: ndarray - The array containing initial guesses for the intensities of each component per + MM : ndarray + The data to be decomposed. Shape is (length_of_signal, number_of_conditions). + Y0 : ndarray + The initial guesses for the component weights at each stretching condition. + Shape is (number of components, number ofconditions) Must be provided if + n_components is not provided. Will override n_components if both are provided. + X0 : ndarray + The initial guesses for the intensities of each component per row/sample/angle. Shape is (length_of_signal, number_of_components). - A: ndarray - The array containing initial guesses for the stretching factor for each component, - at each condition. Shape is (number_of_components, number_of_conditions). - rho: float - The float which sets a stretching factor that influences the decomposition. - Zero corresponds to no stretching present. Relatively insensitive and typically - adjusted in powers of 10. - eta: float - The integer which sets a sparsity factor than influences the decomposition. - Should be set to zero for non sparse data such as PDF. Can be used to improve - results for sparse data such as XRD, but due to instability, should be used - only after first selecting the best value for rho. - max_iter: int + A : ndarray + The initial guesses for the stretching factor for each component, at each + condition. Shape is (number_of_components, number_of_conditions). + rho : float + The stretching factor that influences the decomposition. Zero corresponds to no + stretching present. Relatively insensitive and typically adjusted in powers of 10. + eta : float + The sparsity factor than influences the decomposition. Should be set to zero for + non-sparse data such as PDF. Can be used to improve results for sparse data such + as XRD, but due to instability, should be used only after first selecting the + best value for rho. + max_iter : int The maximum number of times to update each of A, X, and Y before stopping the optimization. - tol: float - The minimum fractional improvement in the objective function to allow - without terminating the optimization. Note that a minimum of 20 updates - are run before this parameter is checked. - n_components: int + tol : float + The convergence threshold. This is the minimum fractional improvement in the + objective function to allow without terminating the optimization. Note that + a minimum of 20 updates are run before this parameter is checked. + n_components : int The number of components to attempt to extract from MM. Note that this will be overridden by Y0 if that is provided, but must be provided if no Y0 is provided. - random_state: int - The integer which acts as a reproducible seed for the initial matrices used in - the optimization. Due to the non-convex nature of the problem, results may vary - even with the same initial guesses, so this does not make the program deterministic. + random_state : int + The seed for the initial matrices used in the optimization. """ self.MM = MM From c783a029aafbfb50cb203479612aa600953b6956 Mon Sep 17 00:00:00 2001 From: John Halloran Date: Fri, 13 Jun 2025 14:36:19 -0700 Subject: [PATCH 06/11] docstring typo --- src/diffpy/snmf/snmf_class.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/diffpy/snmf/snmf_class.py b/src/diffpy/snmf/snmf_class.py index 3650abc0..784bc015 100644 --- a/src/diffpy/snmf/snmf_class.py +++ b/src/diffpy/snmf/snmf_class.py @@ -48,10 +48,10 @@ def __init__( The stretching factor that influences the decomposition. Zero corresponds to no stretching present. Relatively insensitive and typically adjusted in powers of 10. eta : float - The sparsity factor than influences the decomposition. Should be set to zero for + The sparsity factor that influences the decomposition. Should be set to zero for non-sparse data such as PDF. Can be used to improve results for sparse data such as XRD, but due to instability, should be used only after first selecting the - best value for rho. + best value for rho. Suggested adjustment is by powers of 2. max_iter : int The maximum number of times to update each of A, X, and Y before stopping the optimization. From 7f8e33d5f05a57c6682421424cb4b34c816f194a Mon Sep 17 00:00:00 2001 From: John Halloran Date: Sun, 15 Jun 2025 12:47:30 -0700 Subject: [PATCH 07/11] Flag self.rng as private --- src/diffpy/snmf/snmf_class.py | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/src/diffpy/snmf/snmf_class.py b/src/diffpy/snmf/snmf_class.py index 784bc015..14812552 100644 --- a/src/diffpy/snmf/snmf_class.py +++ b/src/diffpy/snmf/snmf_class.py @@ -60,11 +60,12 @@ def __init__( objective function to allow without terminating the optimization. Note that a minimum of 20 updates are run before this parameter is checked. n_components : int - The number of components to attempt to extract from MM. Note that this will + The number of components to extract from MM. Note that this will be overridden by Y0 if that is provided, but must be provided if no Y0 is provided. random_state : int - The seed for the initial matrices used in the optimization. + The seed for the initial guesses at the matrices (A, X, and Y) created by + the decomposition. """ self.MM = MM @@ -76,21 +77,21 @@ def __init__( # Capture matrix dimensions self.N, self.M = MM.shape self.num_updates = 0 - self.rng = np.random.default_rng(random_state) + self._rng = np.random.default_rng(random_state) if Y0 is None: if n_components is None: raise ValueError("Must provide either Y0 or n_components.") else: self.K = n_components - self.Y0 = self.rng.beta(a=2.5, b=1.5, size=(self.K, self.M)) + self.Y0 = self._rng.beta(a=2.5, b=1.5, size=(self.K, self.M)) else: self.K = Y0.shape[0] if self.A is None: - self.A = np.ones((self.K, self.M)) + self.rng.normal(0, 1e-3, size=(self.K, self.M)) + self.A = np.ones((self.K, self.M)) + self._rng.normal(0, 1e-3, size=(self.K, self.M)) if self.X0 is None: - self.X0 = self.rng.random((self.N, self.K)) + self.X0 = self._rng.random((self.N, self.K)) self.X = np.maximum(0, self.X0) self.Y = np.maximum(0, self.Y0) From 6ba837f0479ff0865a3a375466951d0b2dfb35fb Mon Sep 17 00:00:00 2001 From: John Halloran Date: Sun, 15 Jun 2025 13:26:02 -0700 Subject: [PATCH 08/11] Make logic for n_components and Y0 more rigid --- src/diffpy/snmf/main.py | 2 +- src/diffpy/snmf/snmf_class.py | 45 ++++++++++++++++++++--------------- 2 files changed, 27 insertions(+), 20 deletions(-) diff --git a/src/diffpy/snmf/main.py b/src/diffpy/snmf/main.py index a512f67e..9d9721a5 100644 --- a/src/diffpy/snmf/main.py +++ b/src/diffpy/snmf/main.py @@ -7,7 +7,7 @@ Y0 = np.loadtxt("input/W0.txt", dtype=float) N, M = MM.shape -my_model = snmf_class.SNMFOptimizer(MM=MM, Y0=Y0, X0=X0, A=A0, n_components=2) +my_model = snmf_class.SNMFOptimizer(MM=MM, Y0=Y0, X0=X0, A0=A0) print("Done") np.savetxt("my_norm_X.txt", my_model.X, fmt="%.6g", delimiter=" ") np.savetxt("my_norm_Y.txt", my_model.Y, fmt="%.6g", delimiter=" ") diff --git a/src/diffpy/snmf/snmf_class.py b/src/diffpy/snmf/snmf_class.py index 14812552..94291d22 100644 --- a/src/diffpy/snmf/snmf_class.py +++ b/src/diffpy/snmf/snmf_class.py @@ -20,7 +20,7 @@ def __init__( MM, Y0=None, X0=None, - A=None, + A0=None, rho=1e12, eta=610, max_iter=500, @@ -36,12 +36,12 @@ def __init__( The data to be decomposed. Shape is (length_of_signal, number_of_conditions). Y0 : ndarray The initial guesses for the component weights at each stretching condition. - Shape is (number of components, number ofconditions) Must be provided if - n_components is not provided. Will override n_components if both are provided. + Shape is (number_of_components, number_of_conditions) Must provide exactly one + of this or n_components. X0 : ndarray The initial guesses for the intensities of each component per row/sample/angle. Shape is (length_of_signal, number_of_components). - A : ndarray + A0 : ndarray The initial guesses for the stretching factor for each component, at each condition. Shape is (number_of_components, number_of_conditions). rho : float @@ -60,18 +60,14 @@ def __init__( objective function to allow without terminating the optimization. Note that a minimum of 20 updates are run before this parameter is checked. n_components : int - The number of components to extract from MM. Note that this will - be overridden by Y0 if that is provided, but must be provided if no Y0 is - provided. + The number of components to extract from MM. Must be provided when and only when + Y0 is not provided. random_state : int The seed for the initial guesses at the matrices (A, X, and Y) created by the decomposition. """ self.MM = MM - self.X0 = X0 - self.Y0 = Y0 - self.A = A self.rho = rho self.eta = eta # Capture matrix dimensions @@ -79,22 +75,33 @@ def __init__( self.num_updates = 0 self._rng = np.random.default_rng(random_state) + # Enforce exclusive specification of n_components or Y0 + if (n_components is None) == (Y0 is not None): + raise ValueError("Must provide exactly one of Y0 or n_components, but not both.") + + # Initialize Y0 and determine number of components if Y0 is None: - if n_components is None: - raise ValueError("Must provide either Y0 or n_components.") - else: - self.K = n_components - self.Y0 = self._rng.beta(a=2.5, b=1.5, size=(self.K, self.M)) + self.K = n_components + self.Y = self._rng.beta(a=2.5, b=1.5, size=(self.K, self.M)) else: self.K = Y0.shape[0] + self.Y = Y0 + # Initialize A if not provided if self.A is None: self.A = np.ones((self.K, self.M)) + self._rng.normal(0, 1e-3, size=(self.K, self.M)) - if self.X0 is None: - self.X0 = self._rng.random((self.N, self.K)) + else: + self.A = A0 + + # Initialize X0 if not provided + if self.X is None: + self.X = self._rng.random((self.N, self.K)) + else: + self.X = X0 - self.X = np.maximum(0, self.X0) - self.Y = np.maximum(0, self.Y0) + # Enforce non-negativity + self.X = np.maximum(0, self.X) + self.Y = np.maximum(0, self.Y) # Second-order spline: Tridiagonal (-2 on diagonal, 1 on sub/superdiagonals) self.P = 0.25 * diags([1, -2, 1], offsets=[0, 1, 2], shape=(self.M - 2, self.M)) From f45b4bb3fe57f0735a7192d21e15b8bb3f81f460 Mon Sep 17 00:00:00 2001 From: John Halloran Date: Sun, 15 Jun 2025 14:04:14 -0700 Subject: [PATCH 09/11] added class attributes to docstring --- src/diffpy/snmf/snmf_class.py | 126 +++++++++++++++++++++++----------- 1 file changed, 87 insertions(+), 39 deletions(-) diff --git a/src/diffpy/snmf/snmf_class.py b/src/diffpy/snmf/snmf_class.py index 94291d22..950da7cd 100644 --- a/src/diffpy/snmf/snmf_class.py +++ b/src/diffpy/snmf/snmf_class.py @@ -13,6 +13,52 @@ class SNMFOptimizer: For more information on sNMF, please reference: Gu, R., Rakita, Y., Lan, L. et al. Stretched non-negative matrix factorization. npj Comput Mater 10, 193 (2024). https://doi.org/10.1038/s41524-024-01377-5 + + Attributes + ---------- + MM : ndarray + The original, unmodified data to be decomposed and later, compared against. + Shape is (length_of_signal, number_of_conditions). + Y : ndarray + The best guess (or while running, the current guess) for the stretching + factor matrix. + X : ndarray + The best guess (or while running, the current guess) for the matrix of + component intensities. + A : ndarray + The best guess (or while running, the current guess) for the matrix of + component weights. + rho : float + The stretching factor that influences the decomposition. Zero corresponds to no + stretching present. Relatively insensitive and typically adjusted in powers of 10. + eta : float + The sparsity factor that influences the decomposition. Should be set to zero for + non-sparse data such as PDF. Can be used to improve results for sparse data such + as XRD, but due to instability, should be used only after first selecting the + best value for rho. Suggested adjustment is by powers of 2. + max_iter : int + The maximum number of times to update each of A, X, and Y before stopping + the optimization. + tol : float + The convergence threshold. This is the minimum fractional improvement in the + objective function to allow without terminating the optimization. Note that + a minimum of 20 updates are run before this parameter is checked. + n_components : int + The number of components to extract from MM. Must be provided when and only when + Y0 is not provided. + random_state : int + The seed for the initial guesses at the matrices (A, X, and Y) created by + the decomposition. + num_updates : int + The total number of times that any of (A, X, and Y) have had their values changed. + If not terminated by other means, this value is used to stop when reaching max_iter. + objective_function: float + The value corresponding to the minimization of the difference between the MM and the + products of A, X, and Y. For full details see the sNMF paper. Smaller corresponds to + better agreement and is desirable. + objective_difference : float + The change in the objective function value since the last update. A negative value + means that the result improved. """ def __init__( @@ -28,7 +74,7 @@ def __init__( n_components=None, random_state=None, ): - """Initialize an instance of SNMF and run the optimization + """Initialize an instance of SNMF and run the optimization. Parameters ---------- @@ -71,7 +117,7 @@ def __init__( self.rho = rho self.eta = eta # Capture matrix dimensions - self.N, self.M = MM.shape + self._N, self._M = MM.shape self.num_updates = 0 self._rng = np.random.default_rng(random_state) @@ -81,21 +127,21 @@ def __init__( # Initialize Y0 and determine number of components if Y0 is None: - self.K = n_components - self.Y = self._rng.beta(a=2.5, b=1.5, size=(self.K, self.M)) + self._K = n_components + self.Y = self._rng.beta(a=2.5, b=1.5, size=(self._K, self._M)) else: - self.K = Y0.shape[0] + self._K = Y0.shape[0] self.Y = Y0 # Initialize A if not provided if self.A is None: - self.A = np.ones((self.K, self.M)) + self._rng.normal(0, 1e-3, size=(self.K, self.M)) + self.A = np.ones((self._K, self._M)) + self._rng.normal(0, 1e-3, size=(self._K, self._M)) else: self.A = A0 # Initialize X0 if not provided if self.X is None: - self.X = self._rng.random((self.N, self.K)) + self.X = self._rng.random((self._N, self._K)) else: self.X = X0 @@ -104,19 +150,19 @@ def __init__( self.Y = np.maximum(0, self.Y) # Second-order spline: Tridiagonal (-2 on diagonal, 1 on sub/superdiagonals) - self.P = 0.25 * diags([1, -2, 1], offsets=[0, 1, 2], shape=(self.M - 2, self.M)) + self.P = 0.25 * diags([1, -2, 1], offsets=[0, 1, 2], shape=(self._M - 2, self._M)) self.PP = self.P.T @ self.P # Set up residual matrix, objective function, and history self.R = self.get_residual_matrix() self.objective_function = self.get_objective_function() self.objective_difference = None - self.objective_history = [self.objective_function] + self._objective_history = [self.objective_function] # Set up tracking variables for updateX() - self.preX = None - self.GraX = np.zeros_like(self.X) # Gradient of X (zeros for now) - self.preGraX = np.zeros_like(self.X) # Previous gradient of X (zeros for now) + self._preX = None + self._GraX = np.zeros_like(self.X) # Gradient of X (zeros for now) + self._preGraX = np.zeros_like(self.X) # Previous gradient of X (zeros for now) regularization_term = 0.5 * rho * np.linalg.norm(self.P @ self.A.T, "fro") ** 2 sparsity_term = eta * np.sum(np.sqrt(self.X)) # Square root penalty @@ -151,20 +197,20 @@ def __init__( # loop to normalize X # effectively just re-running class with non-normalized X, normalized Y/A as inputs, then only update X # reset difference trackers and initialize - self.preX = None - self.GraX = np.zeros_like(self.X) # Gradient of X (zeros for now) - self.preGraX = np.zeros_like(self.X) # Previous gradient of X (zeros for now) + self._preX = None + self._GraX = np.zeros_like(self.X) # Gradient of X (zeros for now) + self._preGraX = np.zeros_like(self.X) # Previous gradient of X (zeros for now) self.R = self.get_residual_matrix() self.objective_function = self.get_objective_function() self.objective_difference = None - self.objective_history = [self.objective_function] + self._objective_history = [self.objective_function] for norm_iter in range(100): self.updateX() self.R = self.get_residual_matrix() self.objective_function = self.get_objective_function() print(f"Objective function after normX: {self.objective_function:.5e}") - self.objective_history.append(self.objective_function) - self.objective_difference = self.objective_history[-2] - self.objective_history[-1] + self._objective_history.append(self.objective_function) + self.objective_difference = self._objective_history[-2] - self._objective_history[-1] if self.objective_difference < self.objective_function * tol and norm_iter >= 20: break # end of normalization (and program) @@ -172,15 +218,15 @@ def __init__( print("Finished optimization.") def optimize_loop(self): - self.preGraX = self.GraX.copy() + self._preGraX = self._GraX.copy() self.updateX() self.num_updates += 1 self.R = self.get_residual_matrix() self.objective_function = self.get_objective_function() print(f"Objective function after updateX: {self.objective_function:.5e}") - self.objective_history.append(self.objective_function) + self._objective_history.append(self.objective_function) if self.objective_difference is None: - self.objective_difference = self.objective_history[-1] - self.objective_function + self.objective_difference = self._objective_history[-1] - self.objective_function # Now we update Y self.updateY2() @@ -188,7 +234,7 @@ def optimize_loop(self): self.R = self.get_residual_matrix() self.objective_function = self.get_objective_function() print(f"Objective function after updateY2: {self.objective_function:.5e}") - self.objective_history.append(self.objective_function) + self._objective_history.append(self.objective_function) self.updateA2() @@ -196,8 +242,8 @@ def optimize_loop(self): self.R = self.get_residual_matrix() self.objective_function = self.get_objective_function() print(f"Objective function after updateA2: {self.objective_function:.5e}") - self.objective_history.append(self.objective_function) - self.objective_difference = self.objective_history[-2] - self.objective_history[-1] + self._objective_history.append(self.objective_function) + self.objective_difference = self._objective_history[-2] - self._objective_history[-1] def apply_interpolation(self, a, x, return_derivatives=False): """ @@ -469,36 +515,38 @@ def updateX(self): # Compute `AX` using the interpolation function AX, _, _ = self.apply_interpolation_matrix() # Skip the other two outputs # Compute RA and RR - intermediate_RA = AX.flatten(order="F").reshape((self.N * self.M, self.K), order="F") - RA = intermediate_RA.sum(axis=1).reshape((self.N, self.M), order="F") + intermediate_RA = AX.flatten(order="F").reshape((self._N * self._M, self._K), order="F") + RA = intermediate_RA.sum(axis=1).reshape((self._N, self._M), order="F") RR = RA - self.MM # Compute gradient `GraX` - self.GraX = self.apply_transformation_matrix(R=RR).toarray() # toarray equivalent of full, make non-sparse + self._GraX = self.apply_transformation_matrix( + R=RR + ).toarray() # toarray equivalent of full, make non-sparse # Compute initial step size `L0` L0 = np.linalg.eigvalsh(self.Y.T @ self.Y).max() * np.max([self.A.max(), 1 / self.A.min()]) # Compute adaptive step size `L` - if self.preX is None: + if self._preX is None: L = L0 else: - num = np.sum((self.GraX - self.preGraX) * (self.X - self.preX)) # Element-wise multiplication - denom = np.linalg.norm(self.X - self.preX, "fro") ** 2 # Frobenius norm squared + num = np.sum((self._GraX - self._preGraX) * (self.X - self._preX)) # Element-wise multiplication + denom = np.linalg.norm(self.X - self._preX, "fro") ** 2 # Frobenius norm squared L = num / denom if denom > 0 else L0 if L <= 0: L = L0 # Store our old X before updating because it is used in step selection - self.preX = self.X.copy() + self._preX = self.X.copy() while True: # iterate updating X - x_step = self.preX - self.GraX / L + x_step = self._preX - self._GraX / L # Solve x^3 + p*x + q = 0 for the largest real root self.X = np.square(cubic_largest_real_root(-x_step, self.eta / (2 * L))) # Mask values that should be set to zero mask = self.X**2 * L / 2 - L * self.X * x_step + self.eta * np.sqrt(self.X) < 0 self.X = mask * self.X - objective_improvement = self.objective_history[-1] - self.get_objective_function( + objective_improvement = self._objective_history[-1] - self.get_objective_function( R=self.get_residual_matrix() ) @@ -515,9 +563,9 @@ def updateY2(self): Updates Y using matrix operations, solving a quadratic program via `solve_mkr_box`. """ - K = self.K - N = self.N - M = self.M + K = self._K + N = self._N + M = self._M for m in range(M): T = np.zeros((N, K)) # Initialize T as an (N, K) zero matrix @@ -544,9 +592,9 @@ def regularize_function(self, A=None): if A is None: A = self.A - K = self.K - M = self.M - N = self.N + K = self._K + M = self._M + N = self._N # Compute interpolated matrices AX, TX, HX = self.apply_interpolation_matrix(A=A, return_derivatives=True) From bdfab77e003ec8fdc7cf6b038947aa86987f1901 Mon Sep 17 00:00:00 2001 From: Simon Billinge Date: Sun, 15 Jun 2025 17:09:45 -0400 Subject: [PATCH 10/11] fix: cleaner import of SNMFOptimizer --- src/diffpy/snmf/main.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/diffpy/snmf/main.py b/src/diffpy/snmf/main.py index 9d9721a5..ac6cfd46 100644 --- a/src/diffpy/snmf/main.py +++ b/src/diffpy/snmf/main.py @@ -1,5 +1,5 @@ import numpy as np -import snmf_class +from snmf_class import SNMFOptimizer X0 = np.loadtxt("input/X0.txt", dtype=float) MM = np.loadtxt("input/MM.txt", dtype=float) From 4e13df11020b65fb75d8f71ad9a99890fb4f9247 Mon Sep 17 00:00:00 2001 From: Simon Billinge Date: Sun, 15 Jun 2025 17:10:56 -0400 Subject: [PATCH 11/11] fix: correct class instantiation after change in import --- src/diffpy/snmf/main.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/diffpy/snmf/main.py b/src/diffpy/snmf/main.py index ac6cfd46..e38b5ab9 100644 --- a/src/diffpy/snmf/main.py +++ b/src/diffpy/snmf/main.py @@ -7,7 +7,7 @@ Y0 = np.loadtxt("input/W0.txt", dtype=float) N, M = MM.shape -my_model = snmf_class.SNMFOptimizer(MM=MM, Y0=Y0, X0=X0, A0=A0) +my_model = SNMFOptimizer(MM=MM, Y0=Y0, X0=X0, A0=A0) print("Done") np.savetxt("my_norm_X.txt", my_model.X, fmt="%.6g", delimiter=" ") np.savetxt("my_norm_Y.txt", my_model.Y, fmt="%.6g", delimiter=" ")