Usage and Examples¶
Some examples are also available in Jupyter Notebook format in the Examples/ directory of the repository.
Basic Usage¶
from pymcr.mcr import McrAR
mcrar = McrAR()
# MCR assumes a system of the form: D = CS^T
#
# Data that you will provide (hyperspectral context):
# D [n_pixels, n_frequencies] # Hyperspectral image unraveled in space (2D)
#
# initial_spectra [n_components, n_frequencies] ## S^T in the literature
# OR
# initial_conc [n_pixels, n_components] ## C in the literature
# If you have an initial estimate of the spectra
mcrar.fit(D, ST=initial_spectra)
# Otherwise, if you have an initial estimate of the concentrations
mcrar.fit(D, C=initial_conc)
Default: Ordinary Least-Squares (OLS) with Non-Negativity Constraints¶
Using OLS and post-iterative non-negativity constraints is a classic way of performing MCR-ALS. OLS algorithms are often faster in practice, but NNLS method may converge MCR faster.
from pymcr.mcr import McrAR
mcrar = McrAR()
# Equivalent to
mcrar = McrAR(c_regr='OLS', st_regr='OLS')
# Equivalent to
from pymcr.mcr import McrAR
from pymcr.regressors import OLS
from pymcr.constraints import ConstraintNonneg
mcrar = McrAR(c_regr=OLS(), st_regr=OLS(), c_constraints=[ConstraintNonneg()],
st_constraints=[ConstraintNonneg()])
# Then use the fit method with initial guesses of ST or C.
Non-Negative Least-Squares (NNLS), No Constraints¶
Using OLS and post-iterative non-negativity constraints is a classic way of performing MCR-ALS. OLS algorithms are often faster in practice, but NNLS method may converge MCR faster.
from pymcr.mcr import McrAR
from pymcr.regressors import NNLS
mcrar = McrAR(c_regr=NNLS(), st_regr=NNLS(), c_constraints=[], st_constraints=[])
# Then use the fit method with initial guesses of ST or C.
Mix NNLS and OLS Regressors, Non-Negativity and Sum-to-1 (i.e., Norm) Constraints¶
In this example, NNLS is used to regress the signature/spectrum and OLS is used for abundance/concentration. Additionally, the ConstraintNorm imposes that after each iteration the summation of all concentrations for each data point sums-to-one. This is common and physical for many use-cases such as spectroscopy (as long as all analytes have a signature [i.e., no 0-intensity signatures])
from pymcr.mcr import McrAR
from pymcr.regressors import NNLS, OLS
from pymcr.constraints import ConstraintNonneg, ConstraintNorm
# Note constraint order matters
mcrar = McrAR(max_iter=100, st_regr='NNLS', c_regr='OLS',
c_constraints=[ConstraintNonneg(), ConstraintNorm()])
# Equivalent to
# Note constraint order matters
mcrar = McrAR(max_iter=100, st_regr=NNLS(), c_regr=OLS(),
c_constraints=[ConstraintNonneg(), ConstraintNorm()])
# Then use the fit method with initial guesses of ST or C.
Ridge-Regression from Scikit-Learn¶
In this example, NNLS is used to regress the signatures/spectra but a ridge regression regressor is imported from scikit-learn
from pymcr.mcr import McrAR
from pymcr.regressors import NNLS
from pymcr.constraints import ConstraintNonneg, ConstraintNorm
from sklearn.linear_model.ridge import Ridge
# Note that an instance of Ridge can be instantiated (with hyperparameters)
# within the instantiation of McrAR.
mcrar = McrAR(max_iter=100, tol_increase=2, tol_err_change=1e-10,
c_regr=Ridge(alpha=1), st_regr=NNLS(),
c_constraints=[ConstraintNonneg(), ConstraintNorm()])
# Then use the fit method with initial guesses of ST or C.
Single-Gaussian Distribution Spectral Fitting with LMFIT¶
In this example, the LMFIT Library is used to perform Gaussian distribution fitting on each spectral component.
A longer example of this can be found in the Examples/NIST_JRes_Paper Jupyter Notebook.
import numpy as np
from lmfit.models import GaussianModel
from pymcr.mcr import McrAR
from pymcr.constraints import Constraint, ConstraintNorm
class ConstraintSingleGauss(Constraint):
"""
Perform a nonlinear least-squares fitting to enforce a Gaussian.
Parameters
----------
copy : bool
Make copy of input data, A; otherwise, overwrite (if mutable)
axis : int
Axis to perform fitting over
"""
def __init__(self, copy=False, axis=-1):
""" A must be non-negative"""
self.copy = copy
self.axis = axis
def transform(self, A):
""" Fit """
n_components = list(A.shape)
x = np.arange(n_components[self.axis])
n_components.pop(self.axis)
assert len(n_components)==1, "Input must be 2D"
n_components = n_components[0]
A_fit = 0*A
for num in range(n_components):
if (self.axis == -1) | (self.axis == 1):
y = A[num, :]
else:
y = A[:, num]
mod = GaussianModel()
pars = mod.guess(y, x=x)
out = mod.fit(y, pars, x=x)
if (self.axis == -1) | (self.axis == 1):
A_fit[num,:] = 1*out.best_fit
else:
A_fit[:, num] = 1*out.best_fit
if self.copy:
return A_fit
else:
A *= 0
A += A_fit
return A
mcrar = McrAR(max_iter=1000, st_regr='NNLS', c_regr='NNLS',
c_constraints=[ConstraintNorm()],
st_constraints=[ConstraintSingleGauss()])
Semi-Learned Fitting¶
In semi-learned fitting, one has a priori knowledge about some or all of the spectra and/or concentration components. To use this information, the known entities are put into the initial guesses and the “fix” parameter is used to tell the fit method the indices. Effectively the fit method will not alter these values during iterations.
# Example 1: 2 known concentration maps
# D: Mixture Spectra
# C_guess: Concentration guesses, but the 0 and 1 indices are known ahead of time
mcrar.fit(D, C=C_guess, c_fix=[0,1])
# Example 1: 2 known concentration maps, 1 known spectrum
# D: Mixture Spectra
# C_guess: Concentration guesses, but the 0 and 1 indices are known ahead of time
# St_guess: Spectra guesses, but the 0-index is known ahead of time
# c_first: In first iteration, solve for C first (fixed St)?
mcrar.fit(D_known, C=C_guess, ST=St_guess, c_fix=[0,1], st_fix=[0], c_first=True)