Skip to content

gh-599: swap all xp vars to xq (and similar) #600

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 2 commits into from
Apr 11, 2025
Merged

gh-599: swap all xp vars to xq (and similar) #600

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
4027bfd
gh-599: swap all `xp` vars to `xq` (and similar)
paddyroddy Apr 8, 2025
7ebe5df
Merge branch 'main' into paddy/issue-599
paddyroddy Apr 10, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
26 changes: 13 additions & 13 deletions glass/algorithm.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -67,28 +67,28 @@ def nnls(
maxiter = 3 * n

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


Expand Down
18 changes: 9 additions & 9 deletions glass/arraytools.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -91,8 +91,8 @@ def broadcast_leading_axes(

def ndinterp( # noqa: PLR0913
x: float | NDArray[np.float64],
xp: Sequence[float] | NDArray[np.float64],
fp: Sequence[float] | NDArray[np.float64],
xq: Sequence[float] | NDArray[np.float64],
fq: Sequence[float] | NDArray[np.float64],
axis: int = -1,
left: float | None = None,
right: float | None = None,
Expand All @@ -105,16 +105,16 @@ def ndinterp( # noqa: PLR0913
----------
x
The x-coordinates.
xp
xq
The x-coordinates of the data points.
fp
The function values corresponding to the x-coordinates in *xp*.
fq
The function values corresponding to the x-coordinates in *xq*.
axis
The axis to interpolate over.
left
The value to return for x < xp[0].
The value to return for x < xq[0].
right
The value to return for x > xp[-1].
The value to return for x > xq[-1].
period
The period of the function, used for interpolating periodic data.

Expand All @@ -124,9 +124,9 @@ def ndinterp( # noqa: PLR0913

"""
return np.apply_along_axis(
partial(np.interp, x, xp),
partial(np.interp, x, xq),
axis,
fp,
fq,
left=left,
right=right,
period=period,
Expand Down
26 changes: 13 additions & 13 deletions tests/test_arraytools.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -57,40 +57,40 @@ def test_broadcast_leading_axes() -> None:
def test_ndinterp() -> None:
# test 1d interpolation

xp = np.array([0, 1, 2, 3, 4])
yp = np.array([1.1, 1.2, 1.3, 1.4, 1.5])
xq = np.array([0, 1, 2, 3, 4])
yq = np.array([1.1, 1.2, 1.3, 1.4, 1.5])

x: float | NDArray[np.float64] = 0.5
y = glass.arraytools.ndinterp(x, xp, yp)
y = glass.arraytools.ndinterp(x, xq, yq)
assert np.shape(y) == ()
np.testing.assert_allclose(y, 1.15, atol=1e-15)

x = np.array([0.5, 1.5, 2.5])
y = glass.arraytools.ndinterp(x, xp, yp)
y = glass.arraytools.ndinterp(x, xq, yq)
assert np.shape(y) == (3,)
np.testing.assert_allclose(y, [1.15, 1.25, 1.35], atol=1e-15)

x = np.array([[0.5, 1.5], [2.5, 3.5]])
y = glass.arraytools.ndinterp(x, xp, yp)
y = glass.arraytools.ndinterp(x, xq, yq)
assert np.shape(y) == (2, 2)
np.testing.assert_allclose(y, [[1.15, 1.25], [1.35, 1.45]], atol=1e-15)

# test nd interpolation in final axis

yp = np.array([[1.1, 1.2, 1.3, 1.4, 1.5], [2.1, 2.2, 2.3, 2.4, 2.5]])
yq = np.array([[1.1, 1.2, 1.3, 1.4, 1.5], [2.1, 2.2, 2.3, 2.4, 2.5]])

x = 0.5
y = glass.arraytools.ndinterp(x, xp, yp)
y = glass.arraytools.ndinterp(x, xq, yq)
assert np.shape(y) == (2,)
np.testing.assert_allclose(y, [1.15, 2.15], atol=1e-15)

x = np.array([0.5, 1.5, 2.5])
y = glass.arraytools.ndinterp(x, xp, yp)
y = glass.arraytools.ndinterp(x, xq, yq)
assert np.shape(y) == (2, 3)
np.testing.assert_allclose(y, [[1.15, 1.25, 1.35], [2.15, 2.25, 2.35]], atol=1e-15)

x = np.array([[0.5, 1.5], [2.5, 3.5]])
y = glass.arraytools.ndinterp(x, xp, yp)
y = glass.arraytools.ndinterp(x, xq, yq)
assert np.shape(y) == (2, 2, 2)
np.testing.assert_allclose(
y,
Expand All @@ -100,17 +100,17 @@ def test_ndinterp() -> None:

# test nd interpolation in middle axis

yp = np.array(
yq = np.array(
[[[1.1], [1.2], [1.3], [1.4], [1.5]], [[2.1], [2.2], [2.3], [2.4], [2.5]]],
)

x = 0.5
y = glass.arraytools.ndinterp(x, xp, yp, axis=1)
y = glass.arraytools.ndinterp(x, xq, yq, axis=1)
assert np.shape(y) == (2, 1)
np.testing.assert_allclose(y, [[1.15], [2.15]], atol=1e-15)

x = np.array([0.5, 1.5, 2.5])
y = glass.arraytools.ndinterp(x, xp, yp, axis=1)
y = glass.arraytools.ndinterp(x, xq, yq, axis=1)
assert np.shape(y) == (2, 3, 1)
np.testing.assert_allclose(
y,
Expand All @@ -119,7 +119,7 @@ def test_ndinterp() -> None:
)

x = np.array([[0.5, 1.5, 2.5, 3.5], [3.5, 2.5, 1.5, 0.5], [0.5, 3.5, 1.5, 2.5]])
y = glass.arraytools.ndinterp(x, xp, yp, axis=1)
y = glass.arraytools.ndinterp(x, xq, yq, axis=1)
assert np.shape(y) == (2, 3, 4, 1)
np.testing.assert_allclose(
y,
Expand Down