Merge branch 'main' of github.com:diffpy/diffpy.snmf into create_main

Author: aajayi-21

.codecov.yml

@@ -0,0 +1,34 @@
+# codecov can find this file anywhere in the repo, so we don't need to clutter
+# the root folder.
+#comment: false
+
+codecov:
+  notify:
+    require_ci_to_pass: no
+
+coverage:
+  status:
+    patch:
+      default:
+        target: '70'
+        if_no_uploads: error
+        if_not_found: success
+        if_ci_failed: failure
+    project:
+      default: false
+      library:
+        target: auto
+        if_no_uploads: error
+        if_not_found: success
+        if_ci_failed: failure
+        paths: '!*/tests/.*'
+
+      tests:
+        target: 97.9%
+        paths: '*/tests/.*'
+        if_not_found: success
+
+flags:
+  tests:
+    paths:
+      - tests/

.coveragerc

@@ -0,0 +1,13 @@
+[run]
+source = 
+    diffpy/snmf/
+[report]
+omit =
+    */python?.?/*
+    */site-packages/nose/*
+    # ignore _version.py and versioneer.py
+    *version.*
+    *_version.py
+
+exclude_lines =
+    if __name__ == '__main__':

.github/workflows/main.yml

@@ -0,0 +1,59 @@
+name: CI
+
+on:
+  push:
+    branches:
+      - main
+  pull_request:
+  workflow_dispatch:
+
+jobs:
+  miniconda:
+    name: Miniconda ${{ matrix.os }}
+    runs-on: ${{ matrix.os }}
+    strategy:
+        matrix:
+            os: ["ubuntu-latest"]
+    steps:
+      - name: check out diffpy.snmf
+        uses: actions/checkout@v3
+        with:
+          repository: diffpy/diffpy.snmf
+          # for bookkeeping have diffpy.snmf at the same level as everything else in the
+          # directory tree
+          path: .
+
+      - name: initialize miniconda
+        # this uses a marketplace action that sets up miniconda in a way that makes
+        # it easier to use.  I tried setting it up without this and it was a pain
+        uses: conda-incubator/setup-miniconda@v2
+        with:
+          activate-environment: test
+          # environment.yml file is needed by this action.  Because I don't want
+          # maintain this but rather maintain the requirements files it just has
+          # basic things in it like conda and pip
+          environment-file: ./environment.yml
+          python-version: 3
+          auto-activate-base: false
+
+      - name: install diffpy.snmf requirements
+        shell: bash -l {0}
+        run: |
+          conda config --set always_yes yes --set changeps1 no
+          conda config --add channels conda-forge
+          conda activate test
+          conda install --file requirements/run.txt
+          conda install --file requirements/test.txt
+          pip install .
+
+      - name: Validate diffpy.snmf
+        shell: bash -l {0}
+        run: |
+          conda activate test
+          coverage run run_tests.py
+          coverage report -m
+
+      - name: Upload coverage reports to Codecov
+        uses: codecov/codecov-action@v3
+#        with:
+#          token: ${{ secrets.CODECOV_TOKEN }}

README.md

@@ -1 +1,4 @@
-# diffpy.snmf
+# diffpy.snmf
+
+[![test](https://github.com/diffpy/diffpy.snmf/actions/workflows/main.yml/badge.svg)](https://github.com/diffpy/diffpy.snmf/actions/workflows/main.yml)
+[![codecov](https://codecov.io/gh/diffpy/diffpy.snmf/branch/main/graph/badge.svg)](https://codecov.io/gh/diffpy/diffpy.snmf)

diffpy/snmf/factorizers.py

@@ -0,0 +1,32 @@
+import numpy as np
+import scipy.optimize
+
+
+def lsqnonneg(stretched_component_matrix, target_signal):
+    """Finds the weights of stretched component signals under one-sided constraint.
+
+    Solves ``argmin_x || Ax - b ||_2`` for ``x>=0`` where A is the stretched_component_matrix and b is the target_signal
+    vector. Finds the weights of component signals given undecomposed signal data and stretched components under a
+    one-sided constraint on the weights.
+
+    Parameters
+    ----------
+    stretched_component_matrix: 2d array like
+      The component matrix where each column contains a stretched component signal. Has dimensions R x C where R is
+      the length of the signal and C is the number of components. Does not need to be nonnegative. Corresponds with 'A'
+      from the objective function.
+
+    target_signal: 1d array like
+      The signal that is used as reference against which weight factors will be determined. Any column from the matrix
+      of the entire, unfactorized input data could be used. Has length R. Does not need to be nonnegative. Corresponds
+      with 'b' from the objective function.
+
+    Returns
+    -------
+    1d array like
+      The vector containing component signal weights at a moment. Has length C.
+
+    """
+    stretched_component_matrix = np.asarray(stretched_component_matrix)
+    target_signal = np.asarray(target_signal)
+    return scipy.optimize.nnls(stretched_component_matrix, target_signal)[0]

diffpy/snmf/io.py

@@ -0,0 +1,116 @@
+import numpy as np
+import scipy.sparse
+from pathlib import Path
+from diffpy.utils.parsers.loaddata import loadData
+
+
+def initialize_variables(data_input, component_amount, data_type, sparsity=1, smoothness=1e18):
+    """Determines the variables and initial values used in the SNMF algorithm.
+
+    Parameters
+    ----------
+    data_input: 2d array like
+      The observed or simulated PDF or XRD data provided by the user. Has dimensions R x N where R is the signal length
+      and N is the number of PDF/XRD signals.
+
+    component_amount: int
+      The number of component signals the user would like to decompose 'data_input' into.
+
+    data_type: str
+      The type of data the user has passed into the program. Can assume the value of 'PDF' or 'XRD.'
+
+    sparsity: float, optional
+      The regularization parameter that behaves as the coefficient of a "sparseness" regularization term that enhances
+      the ability to decompose signals in the case of sparse data e.g. X-ray Diffraction data. A non-zero value
+      indicates sparsity in the data; greater magnitudes indicate greater amounts of sparsity.
+
+    smoothness: float, optional
+      The regularization parameter that behaves as the coefficient of a "smoothness" term that ensures that component
+      signal weightings change smoothly with time. Assumes a default value of 1e18.
+
+    Returns
+    -------
+    dictionary
+      The collection of the names and values of the constants used in the algorithm. Contains the number of observed PDF
+      /XRD patterns, the length of each pattern, the type of the data, the number of components the user would like to
+      decompose the data into, an initial guess for the component matrix, and initial guess for the weight factor matrix
+      ,an initial guess for the stretching factor matrix, a parameter controlling smoothness of the solution, a
+      parameter controlling sparseness of the solution, the matrix representing the smoothness term, and a matrix used
+      to construct a hessian matrix.
+
+    """
+    signal_length = data_input.shape[0]
+    moment_amount = data_input.shape[1]
+
+    component_matrix_guess = np.random.rand(signal_length, component_amount)
+    weight_matrix_guess = np.random.rand(component_amount, moment_amount)
+    stretching_matrix_guess = np.ones(component_amount, moment_amount) + np.random.randn(component_amount,
+                                                                                         moment_amount) * 1e-3
+
+    diagonals = [np.ones(moment_amount - 2), -2 * np.ones(moment_amount - 2), np.ones(moment_amount - 2)]
+    smoothness_term = .25 * scipy.sparse.diags(diagonals, [0, 1, 2], shape=(moment_amount - 2, moment_amount))
+
+    hessian_helper_matrix = scipy.sparse.block_diag([smoothness_term.T @ smoothness_term] * component_amount)
+    sequence = np.arange(moment_amount * component_amount).reshape(component_amount, moment_amount).T.flatten()
+    hessian_helper_matrix = hessian_helper_matrix[sequence, :][:, sequence]
+
+    return {
+        "signal_length": signal_length,
+        "moment_amount": moment_amount,
+        "component_matrix_guess": component_matrix_guess,
+        "weight_matrix_guess": weight_matrix_guess,
+        "stretching_matrix_guess": stretching_matrix_guess,
+        "component_amount": component_amount,
+        "data_type": data_type,
+        "smoothness": smoothness,
+        "sparsity": sparsity,
+        "smoothness_term": smoothness_term,
+        "hessian_helper_matrix": hessian_helper_matrix
+    }
+
+
+def load_input_signals(file_path=None):
+    """Processes a directory of a series of PDF/XRD patterns into a usable format.
+
+    Constructs a 2d array out of a directory of PDF/XRD patterns containing each files dependent variable column in a
+    new column. Constructs a 1d array containing the grid values.
+
+    Parameters
+    ----------
+    file_path: str or Path object, optional
+      The path to the directory containing the input XRD/PDF data. If no path is specified, defaults to the current
+      working directory. Accepts a string or a pathlib.Path object. Input data not on the same grid as the first file
+      read will be ignored.
+
+    Returns
+    -------
+    tuple
+      The tuple whose first element is an R x M 2d array made of PDF/XRD patterns as each column; R is the length of the
+      signal and M is the number of patterns. The tuple contains a 1d array containing the values of the grid points as
+      its second element; Has length R.
+
+    """
+
+    if file_path is None:
+        directory_path = Path.cwd()
+    else:
+        directory_path = Path(file_path)
+
+    values_list = []
+    grid_list = []
+    current_grid = []
+    for item in directory_path.iterdir():
+        if item.is_file():
+            data = loadData(item.resolve())
+            if current_grid and current_grid != data[:, 0]:
+                print(f"{item.name} was ignored as it is not on a compatible grid.")
+                continue
+            else:
+                grid_list.append(data[:, 0])
+                current_grid = grid_list[-1]
+                values_list.append(data[:, 1])
+
+    grid_array = np.column_stack(grid_list)
+    grid_vector = np.unique(grid_array, axis=1)
+    values_array = np.column_stack(values_list)
+    return grid_vector, values_array

diffpy/snmf/optimizers.py

@@ -0,0 +1,52 @@
+import numpy as np
+import cvxpy
+
+
+def get_weights(stretched_component_gram_matrix, linear_coefficient, lower_bound, upper_bound):
+    """Finds the weights of stretched component signals under a two-sided constraint
+
+    Solves min J(y) = (linear_coefficient)' * y + (1/2) * y' * (quadratic coefficient) * y where lower_bound <= y <=
+    upper_bound and stretched_component_gram_matrix is symmetric positive definite. Finds the weightings of stretched
+    component signals under a two-sided constraint.
+
+    Parameters
+    ----------
+    stretched_component_gram_matrix: 2d array like
+      The Gram matrix constructed from the stretched component matrix. It is a square positive definite matrix. It has
+      dimensions C x C where C is the number of component signals. Must be symmetric positive definite.
+
+    linear_coefficient: 1d array like
+      The vector containing the product of the stretched component matrix and the transpose of the observed data matrix.
+      Has length C.
+
+    lower_bound: 1d array like
+      The lower bound on the values of the output weights. Has the same dimensions of the function output. Each 
+      element in 'lower_bound' determines the minimum value the corresponding element in the function output may take.
+
+    upper_bound: 1d array like
+      The upper bound on the values of the output weights. Has the same dimensions of the function output. Each element
+      in 'upper_bound' determines the maximum value the corresponding element in the function output may take.
+
+    Returns
+    -------
+    1d array like
+      The vector containing the weightings of the components needed to reconstruct a given input signal from the input
+      set. Has length C
+
+    """
+    stretched_component_gram_matrix = np.asarray(stretched_component_gram_matrix)
+    linear_coefficient = np.asarray(linear_coefficient)
+    upper_bound = np.asarray(upper_bound)
+    lower_bound = np.asarray(lower_bound)
+
+    problem_size = max(linear_coefficient.shape)
+    solution_variable = cvxpy.Variable(problem_size)
+
+    objective = cvxpy.Minimize(
+        linear_coefficient.T @ solution_variable + 0.5 * cvxpy.quad_form(solution_variable,
+                                                                         stretched_component_gram_matrix))
+    constraints = [lower_bound <= solution_variable, solution_variable <= upper_bound]
+
+    cvxpy.Problem(objective, constraints).solve()
+
+    return solution_variable.value

diffpy/snmf/polynomials.py

@@ -0,0 +1,33 @@
+import numpy as np
+
+
+def rooth(linear_coefficient, constant_term):
+    """
+    Returns the largest real root of x^3+(linear_coefficient) * x + constant_term. If there are no real roots return 0.
+
+    Parameters
+    ----------
+    linear_coefficient: nd array like of floats
+        The matrix coefficient of the linear term
+    constant_term: 0d array like, 1d array like of floats or scalar
+        The constant scalar term of the problem
+
+    Returns
+    -------
+    ndarray of floats
+        The largest real root of x^3+(linear_coefficient) * x + constant_term if roots are real, else return 0 array
+
+
+    """
+    linear_coefficient = np.asarray(linear_coefficient)
+    constant_term = np.asarray(constant_term)
+    solution = np.empty_like(linear_coefficient, dtype=np.float64)
+
+    for index, value in np.ndenumerate(linear_coefficient):
+        inputs = [1, 0, value, constant_term]
+        roots = np.roots(inputs)
+        if ((constant_term / 2) ** 2 + (value / 3) ** 3) < 0:  # Discriminant of depressed cubic equation
+            solution[index] = max(np.real(roots))
+        else:
+            solution[index] = 0
+    return solution