Skip to content

Fixes for regularization due to scipy depreciations #60

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
florisvb merged 2 commits into from
Feb 13, 2025
Merged

Fixes for regularization due to scipy depreciations #60

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
6a36b4d
updated regularization for scipy >=1.12
asteppke Feb 13, 2025
c2cbdd0
fixed cumtrapz depreciation
asteppke Feb 13, 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
2 changes: 1 addition & 1 deletion pynumdiff/linear_model/_linear_model.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -414,7 +414,7 @@ def __integrate_dxdt_hat_matrix__(dxdt_hat, dt):
#assert isinstance(dxdt_hat, np.matrix)
if len(dxdt_hat.shape) == 1:
dxdt_hat = np.reshape(dxdt_hat, [1, len(dxdt_hat)])
x = np.array(scipy.integrate.cumtrapz(dxdt_hat, axis=1))
x = np.array(scipy.integrate.cumulative_trapezoid(dxdt_hat, axis=1))
first_value = x[:, 0:1] - np.mean(dxdt_hat[:, 0:1], axis=1).reshape(dxdt_hat.shape[0], 1)
x = np.hstack((first_value, x))*dt
return x
Expand Down
20 changes: 10 additions & 10 deletions pynumdiff/total_variation_regularization/__chartrand_tvregdiff__.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -228,7 +228,7 @@ def AT(w): return (sum(w) * np.ones(n + 1) -
g = AT(A(u)) + ATb + alph * L * u

# Prepare to solve linear equation.
tol = 1e-4
rtol = 1e-4
# Simple preconditioner.
P = alph * sparse.spdiags(L.diagonal() + 1, 0, n + 1, n + 1)

Expand All @@ -237,8 +237,8 @@ def linop(v): return (alph * L * v + AT(A(v)))

if diagflag:
[s, info_i] = sparse.linalg.cg(
linop, g, x0=None, tol=tol, maxiter=maxit, callback=None,
M=P, atol='legacy')
linop, g, x0=None, rtol=rtol, maxiter=maxit, callback=None,
M=P, atol=0)
#print('iteration {0:4d}: relative change = {1:.3e}, '
# 'gradient norm = {2:.3e}\n'.format(ii,
# np.linalg.norm(
Expand All @@ -251,8 +251,8 @@ def linop(v): return (alph * L * v + AT(A(v)))
# print("WARNING - illegal input or breakdown")
else:
[s, info_i] = sparse.linalg.cg(
linop, g, x0=None, tol=tol, maxiter=maxit, callback=None,
M=P, atol='legacy')
linop, g, x0=None, rtol=rtol, maxiter=maxit, callback=None,
M=P, atol=0)
# Update solution.
u = u - s
# Display plot.
Expand Down Expand Up @@ -300,16 +300,16 @@ def AT(w): return (sum(w) * np.ones(len(w)) -
# droptol = 1.0e-2
R = sparse.dia_matrix(np.linalg.cholesky(B.todense()))
# Prepare to solve linear equation.
tol = 1.0e-4
rtol = 1.0e-4

def linop(v): return (alph * L * v + AT(A(v)))
linop = splin.LinearOperator((n, n), linop)

print(maxit)
if diagflag:
[s, info_i] = sparse.linalg.cg(
linop, -g, x0=None, tol=tol, maxiter=maxit, callback=None,
M=np.dot(R.transpose(), R), atol='legacy')
linop, -g, x0=None, rtol=rtol, maxiter=maxit, callback=None,
M=np.dot(R.transpose(), R), atol=0)
print('iteration {0:4d}: relative change = {1:.3e}, '
'gradient norm = {2:.3e}\n'.format(ii,
np.linalg.norm(s[0]) /
Expand All @@ -322,8 +322,8 @@ def linop(v): return (alph * L * v + AT(A(v)))

else:
[s, info_i] = sparse.linalg.cg(
linop, -g, x0=None, tol=tol, maxiter=maxit, callback=None,
M=np.dot(R.transpose(), R), atol='legacy')
linop, -g, x0=None, rtol=rtol, maxiter=maxit, callback=None,
M=np.dot(R.transpose(), R), atol=0)
# Update current solution
u = u + s
# Display plot.
Expand Down
4 changes: 2 additions & 2 deletions pynumdiff/utils/utility.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -213,7 +213,7 @@ def finite_difference(x, dt):
# Trapazoidal integration, with interpolated final point so that the lengths match.
def integrate_dxdt_hat(dxdt_hat, dt):
"""
Wrapper for scipy.integrate.cumtrapz to integrate dxdt_hat that ensures the integral has the same length
Wrapper for scipy.integrate.cumulative_trapezoid to integrate dxdt_hat that ensures the integral has the same length

:param dxdt_hat: estimate derivative of timeseries
:type dxdt_hat: np.array
Expand All @@ -224,7 +224,7 @@ def integrate_dxdt_hat(dxdt_hat, dt):
:return: integral of dxdt_hat
:rtype: np.array
"""
x = scipy.integrate.cumtrapz(dxdt_hat)
x = scipy.integrate.cumulative_trapezoid(dxdt_hat)
first_value = x[0] - np.mean(dxdt_hat[0:1])
x = np.hstack((first_value, x))*dt
return x
Expand Down
Loading