13416593565

Committed 19 Feb 2025 03:44PM UTC coverage: 94.408%. Remained the same

Build # 13416593565

Build Type

push

github

Committed by

web-flow

Commit Message

gh-478: glass.core.array -> glass.arraytools (#524)

Run Details

167 of 171 branches covered (97.66%)

Branch coverage included in aggregate %.

13 of 13 new or added lines in 4 files covered. (100.0%)

25 existing lines in 3 files now uncovered.

1133 of 1206 relevant lines covered (93.95%)

7.51 hits per line

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

96.84

/glass/algorithm.py

"""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 [Lawson95]_ as described by
    [Bro97]_.

    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.

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

glass-dev / glass / 13416593565

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