13335643083

Committed 14 Feb 2025 07:00PM UTC coverage: 93.028% (-1.3%) from 94.286%

Build # 13335643083

Build Type

Pull #497

github

Committed by

web-flow

Commit Message

Merge ff2a741bb into 4497645a8

Pull Request Pull Request #497: gh-496: functions for legacy mode

Run Details

165 of 171 branches covered (96.49%)

Branch coverage included in aggregate %.

57 of 77 new or added lines in 2 files covered. (74.03%)

5 existing lines in 1 file now uncovered.

1116 of 1206 relevant lines covered (92.54%)

7.4 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

92.63

/glass/core/algorithm.py

"""Core module for algorithms."""

from __future__ import annotations

from typing import TYPE_CHECKING

import numpy as np

if TYPE_CHECKING:
    from numpy.typing import NDArray


def nnls(
    a: NDArray[np.float64],
    b: NDArray[np.float64],
    *,
    tol: float = 0.0,
    maxiter: int | None = None,
) -> NDArray[np.float64]:
    """
    Compute a non-negative least squares solution.

    Implementation of the algorithm due to [1] as described in [2].

    Parameters
    ----------
    a
        The matrix.
    b
        The vector.
    tol
        The tolerance for convergence.
    maxiter
        The maximum number of iterations.

    Returns
    -------
        The non-negative least squares solution.

    Raises
    ------
    ValueError
        If ``a`` is not a matrix.
    ValueError
        If ``b`` is not a vector.
    ValueError
        If the shapes of ``a`` and ``b`` do not match.

    References
    ----------
    * [1] Lawson, C. L. and Hanson, R. J. (1995), Solving Least Squares
          Problems. doi: 10.1137/1.9781611971217
    * [2] Bro, R. and De Jong, S. (1997), A fast
          non-negativity-constrained least squares algorithm. J.
          Chemometrics, 11, 393-401.

    """
    a = np.asanyarray(a)
    b = np.asanyarray(b)

    if a.ndim != 2:
        msg = "input `a` is not a matrix"
        raise ValueError(msg)
    if b.ndim != 1:
        msg = "input `b` is not a vector"
        raise ValueError(msg)
    if a.shape[0] != b.shape[0]:
        msg = "the shapes of `a` and `b` do not match"
        raise ValueError(msg)

    _, n = a.shape

    if maxiter is None:
        maxiter = 3 * n

    index = np.arange(n)
    p = np.full(n, fill_value=False)
    x = np.zeros(n)
    for _ in range(maxiter):
        if np.all(p):
            break
        w = np.dot(b - a @ x, a)
        m = index[~p][np.argmax(w[~p])]
        if w[m] <= tol:
            break
        p[m] = True
        while True:
            ap = a[:, p]
            xp = x[p]
            sp = np.linalg.solve(ap.T @ ap, b @ ap)
            t = sp <= 0
            if not np.any(t):
                break
            alpha = -np.min(xp[t] / (xp[t] - sp[t]))
            x[p] += alpha * (sp - xp)
            p[x <= 0] = False
        x[p] = sp
        x[~p] = 0
    return x


def cov_clip(
    cov: NDArray[np.float64],
    rtol: float | None = None,
) -> NDArray[np.float64]:
    """
    Covariance matrix from clipping non-positive eigenvalues.

    The relative tolerance *rtol* is defined as for
    :func:`~array_api.linalg.matrix_rank`.

    Parameter
    ---------
    cov
        A symmetric matrix (or a stack of matrices).
    rtol
        An optional relative tolerance for eigenvalues to be considered
        positive.

    Returns
    -------
        Covariance matrix with negative eigenvalues clipped.

    """
    xp = cov.__array_namespace__()

    # Hermitian eigendecomposition
    w, v = xp.linalg.eigh(cov)

    # get tolerance if not given
    if rtol is None:
        rtol = max(v.shape[-2], v.shape[-1]) * xp.finfo(w.dtype).eps

    # clip negative diagonal values
    w = xp.clip(w, rtol * xp.max(w, axis=-1, keepdims=True), None)

    # put matrix back together
    # enforce symmetry
    v = xp.sqrt(w[..., None, :]) * v
    cov_clipped: NDArray[np.float64] = xp.matmul(v, xp.matrix_transpose(v))
    return cov_clipped


def nearcorr(
    a: NDArray[np.float64],
    *,
    tol: float | None = None,
    niter: int = 100,
) -> NDArray[np.float64]:
    """
    Compute the nearest correlation matrix.

    Returns the nearest correlation matrix using the alternating
    projections algorithm of [Higham02]_.

    Parameters
    ----------
    a
        Square matrix (or a stack of square matrices).
    tol
        Tolerance for convergence. Default is 16 times machine epsilon.
    niter
        Maximum number of iterations.

    Returns
    -------
        Nearest correlation matrix.

    """
    xp = a.__array_namespace__()

    # shorthand for Frobenius norm
    frob = xp.linalg.matrix_norm

    # get size of the covariance matrix and flatten leading dimensions
    *dim, m, n = a.shape
    if m != n:
        msg = "non-square matrix"
        raise ValueError(msg)

    # default tolerance
    if tol is None:
        tol = n * xp.finfo(a.dtype).eps

    # current result, flatten leading dimensions
    y = a.reshape(-1, n, n)

    # initial correction is zero
    ds = xp.zeros_like(a)

    # store identity matrix
    diag = xp.eye(n)

    # find the nearest correlation matrix
    for _ in range(niter):
        # apply Dykstra's correction to current result
        r = y - ds

        # project onto positive semi-definite matrices
        x = cov_clip(r)

        # compute Dykstra's correction
        ds = x - r

        # project onto matrices with unit diagonal
        y = (1 - diag) * x + diag

        # check for convergence
        if xp.all(frob(y - x) <= tol * frob(y)):
            break

    # return result in original shape
    return y.reshape(*dim, n, n)


def cov_nearest(
    cov: NDArray[np.float64],
    tol: float | None = None,
    niter: int = 100,
) -> NDArray[np.float64]:
    """
    Covariance matrix from nearest correlation matrix.

    Divides *cov* along rows and columns by the square root of the
    diagonal, then computes the nearest valid correlation matrix using
    :func:`nearcorr`, before scaling rows and columns back.  The
    diagonal of the input is hence unchanged.

    Parameters
    ----------
    cov
        A square matrix (or a stack of matrices).
    tol
        Tolerance for convergence, see :func:`nearcorr`.
    niter
        Maximum number of iterations.

    Returns
    -------
        Covariance matrix from nearest correlation matrix.

    """
    xp = cov.__array_namespace__()

    # get the diagonal
    diag = xp.linalg.diagonal(cov)

    # cannot fix negative diagonal
    if xp.any(diag < 0):
        msg = "negative values on the diagonal"
        raise ValueError(msg)

    # store the normalisation of the matrix
    norm: NDArray[np.float64] = xp.sqrt(diag)
    norm = norm[..., None, :] * norm[..., :, None]

    # find nearest correlation matrix
    corr = cov / xp.where(norm > 0, norm, 1.0)
    return nearcorr(corr, niter=niter, tol=tol) * norm

1	"""Core module for algorithms."""
2
3	from __future__ import annotations	8✔
4
5	from typing import TYPE_CHECKING	8✔
6
7	import numpy as np	8✔
8
9	if TYPE_CHECKING:
10	from numpy.typing import NDArray
11
12
13	def nnls(	8✔
14	a: NDArray[np.float64],
15	b: NDArray[np.float64],
16	*,
17	tol: float = 0.0,
18	maxiter: int \| None = None,
19	) -> NDArray[np.float64]:
20	"""
21	Compute a non-negative least squares solution.
22
23	Implementation of the algorithm due to [1] as described in [2].
24
25	Parameters
26	----------
27	a
28	The matrix.
29	b
30	The vector.
31	tol
32	The tolerance for convergence.
33	maxiter
34	The maximum number of iterations.
35
36	Returns
37	-------
38	The non-negative least squares solution.
39
40	Raises
41	------
42	ValueError
43	If ``a`` is not a matrix.
44	ValueError
45	If ``b`` is not a vector.
46	ValueError
47	If the shapes of ``a`` and ``b`` do not match.
48
49	References
50	----------
51	* [1] Lawson, C. L. and Hanson, R. J. (1995), Solving Least Squares
52	Problems. doi: 10.1137/1.9781611971217
53	* [2] Bro, R. and De Jong, S. (1997), A fast
54	non-negativity-constrained least squares algorithm. J.
55	Chemometrics, 11, 393-401.
56
57	"""
58	a = np.asanyarray(a)	8✔
59	b = np.asanyarray(b)	8✔
60
61	if a.ndim != 2:	8✔
62	msg = "input `a` is not a matrix"	8✔
63	raise ValueError(msg)	8✔
64	if b.ndim != 1:	8✔
65	msg = "input `b` is not a vector"	8✔
66	raise ValueError(msg)	8✔
67	if a.shape[0] != b.shape[0]:	8✔
68	msg = "the shapes of `a` and `b` do not match"	8✔
69	raise ValueError(msg)	8✔
70
71	_, n = a.shape	8✔
72
73	if maxiter is None:	8✔
74	maxiter = 3 * n	8✔
75
76	index = np.arange(n)	8✔
77	p = np.full(n, fill_value=False)	8✔
78	x = np.zeros(n)	8✔
79	for _ in range(maxiter):	8✔
80	if np.all(p):	8✔
81	break	8✔
82	w = np.dot(b - a @ x, a)	8✔
83	m = index[~p][np.argmax(w[~p])]	8✔
84	if w[m] <= tol:	8✔
85	break	8✔
86	p[m] = True	8✔
87	while True:	8✔
88	ap = a[:, p]	8✔
89	xp = x[p]	8✔
90	sp = np.linalg.solve(ap.T @ ap, b @ ap)	8✔
91	t = sp <= 0	8✔
92	if not np.any(t):	8✔
93	break	8✔
94	alpha = -np.min(xp[t] / (xp[t] - sp[t]))	×
95	x[p] += alpha * (sp - xp)	×
96	p[x <= 0] = False	×
97	x[p] = sp	8✔
98	x[~p] = 0	8✔
99	return x	8✔
100
101
102	def cov_clip(	8✔
103	cov: NDArray[np.float64],
104	rtol: float \| None = None,
105	) -> NDArray[np.float64]:
106	"""
107	Covariance matrix from clipping non-positive eigenvalues.
108
109	The relative tolerance rtol is defined as for
110	:func:`~array_api.linalg.matrix_rank`.
111
112	Parameter
113	---------
114	cov
115	A symmetric matrix (or a stack of matrices).
116	rtol
117	An optional relative tolerance for eigenvalues to be considered
118	positive.
119
120	Returns
121	-------
122	Covariance matrix with negative eigenvalues clipped.
123
124	"""
125	xp = cov.__array_namespace__()	8✔
126
127	# Hermitian eigendecomposition
128	w, v = xp.linalg.eigh(cov)	8✔
129
130	# get tolerance if not given
131	if rtol is None:	8✔
132	rtol = max(v.shape[-2], v.shape[-1]) * xp.finfo(w.dtype).eps	8✔
133
134	# clip negative diagonal values
135	w = xp.clip(w, rtol * xp.max(w, axis=-1, keepdims=True), None)	8✔
136
137	# put matrix back together
138	# enforce symmetry
139	v = xp.sqrt(w[..., None, :]) * v	8✔
140	cov_clipped: NDArray[np.float64] = xp.matmul(v, xp.matrix_transpose(v))	8✔
141	return cov_clipped	8✔
142
143
144	def nearcorr(	8✔
145	a: NDArray[np.float64],
146	*,
147	tol: float \| None = None,
148	niter: int = 100,
149	) -> NDArray[np.float64]:
150	"""
151	Compute the nearest correlation matrix.
152
153	Returns the nearest correlation matrix using the alternating
154	projections algorithm of [Higham02]_.
155
156	Parameters
157	----------
158	a
159	Square matrix (or a stack of square matrices).
160	tol
161	Tolerance for convergence. Default is 16 times machine epsilon.
162	niter
163	Maximum number of iterations.
164
165	Returns
166	-------
167	Nearest correlation matrix.
168
169	"""
170	xp = a.__array_namespace__()	8✔
171
172	# shorthand for Frobenius norm
173	frob = xp.linalg.matrix_norm	8✔
174
175	# get size of the covariance matrix and flatten leading dimensions
176	*dim, m, n = a.shape	8✔
177	if m != n:	8✔
NEW 178	msg = "non-square matrix"	×
NEW 179	raise ValueError(msg)	×
180
181	# default tolerance
182	if tol is None:	8✔
183	tol = n * xp.finfo(a.dtype).eps	8✔
184
185	# current result, flatten leading dimensions
186	y = a.reshape(-1, n, n)	8✔
187
188	# initial correction is zero
189	ds = xp.zeros_like(a)	8✔
190
191	# store identity matrix
192	diag = xp.eye(n)	8✔
193
194	# find the nearest correlation matrix
195	for _ in range(niter):	8✔
196	# apply Dykstra's correction to current result
197	r = y - ds	8✔
198
199	# project onto positive semi-definite matrices
200	x = cov_clip(r)	8✔
201
202	# compute Dykstra's correction
203	ds = x - r	8✔
204
205	# project onto matrices with unit diagonal
206	y = (1 - diag) * x + diag	8✔
207
208	# check for convergence
209	if xp.all(frob(y - x) <= tol * frob(y)):	8✔
210	break	8✔
211
212	# return result in original shape
213	return y.reshape(*dim, n, n)	8✔
214
215
216	def cov_nearest(	8✔
217	cov: NDArray[np.float64],
218	tol: float \| None = None,
219	niter: int = 100,
220	) -> NDArray[np.float64]:
221	"""
222	Covariance matrix from nearest correlation matrix.
223
224	Divides cov along rows and columns by the square root of the
225	diagonal, then computes the nearest valid correlation matrix using
226	:func:`nearcorr`, before scaling rows and columns back. The
227	diagonal of the input is hence unchanged.
228
229	Parameters
230	----------
231	cov
232	A square matrix (or a stack of matrices).
233	tol
234	Tolerance for convergence, see :func:`nearcorr`.
235	niter
236	Maximum number of iterations.
237
238	Returns
239	-------
240	Covariance matrix from nearest correlation matrix.
241
242	"""
243	xp = cov.__array_namespace__()	8✔
244
245	# get the diagonal
246	diag = xp.linalg.diagonal(cov)	8✔
247
248	# cannot fix negative diagonal
249	if xp.any(diag < 0):	8✔
NEW 250	msg = "negative values on the diagonal"	×
NEW 251	raise ValueError(msg)	×
252
253	# store the normalisation of the matrix
254	norm: NDArray[np.float64] = xp.sqrt(diag)	8✔
255	norm = norm[..., None, :] * norm[..., :, None]	8✔
256
257	# find nearest correlation matrix
258	corr = cov / xp.where(norm > 0, norm, 1.0)	8✔
259	return nearcorr(corr, niter=niter, tol=tol) * norm	8✔

glass-dev / glass / 13335643083

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous