• Home
  • Features
  • Pricing
  • Docs
  • Announcements
  • Sign In

glass-dev / glass / 13416593565

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

push

github

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

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
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✔
STATUS · Troubleshooting · Open an Issue · Sales · Support · CAREERS · ENTERPRISE · START FREE · SCHEDULE DEMO
ANNOUNCEMENTS · TWITTER · TOS & SLA · Supported CI Services · What's a CI service? · Automated Testing

© 2026 Coveralls, Inc