1D Post-stack inversion with proximal solvers¶
In [1]:
%load_ext autoreload
%autoreload 2
%matplotlib inline
import os
import warnings
warnings.filterwarnings('ignore')
import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import filtfilt
from scipy.linalg import solve
from scipy.sparse.linalg import lsqr
from matplotlib.colors import LinearSegmentedColormap
from pylops.basicoperators import *
from pylops.signalprocessing import *
from pylops.utils.wavelets import *
from pylops.avo.poststack import *
from pylops.optimization.sparsity import *
from pyproximal.proximal import *
from pyproximal.proximal import VStack as VStackprox
from pyproximal import ProxOperator
from pyproximal.optimization.primal import *
from pyproximal.optimization.primaldual import *
from pyproximal.optimization.bregman import *
from pyproximal.optimization.segmentation import *
In [2]:
os.environ['NUMBA_NUM_THREADS']
Out[2]:
'4'
In [3]:
def RRE(x, xinv):
return np.linalg.norm(x-xinv) / np.linalg.norm(x)
In [4]:
def PSNR(x, xinv):
return 10 * np.log10(len(xinv) * np.max(xinv) / np.linalg.norm(x-xinv))
We consider the following problem
$$ \mathbf{x} = arg min_\mathbf{x} \frac{\sigma}{2} ||\mathbf{W}\mathbf{D}\mathbf{x} - \mathbf{y} ||_2^2 + ||\mathbf{D} \mathbf{x}||_1 $$True proximal of convolution operator¶
In [5]:
nt = 101
dt = 0.004
t = np.arange(nt)*dt
wav, th, wavc = ricker(t[:11], f0=20)
x = 1800 * np.ones(nt)
x[nt//4:nt//2] = 2000
x[nt//2:2*nt//3] = 1700
x = np.log(x)
# Add gradient at the end
#istart = nt-35
#iend = nt-20
#xgrad = np.arange(iend - istart) * (x[0] - x[istart-1]) / (iend - istart) + x[istart-1]
#x[istart:iend] = xgrad
#x[iend:] = x[0]
In [6]:
# Background model
x0 = filtfilt(np.ones(30)/30, 1, x)
# Data
Lop = PoststackLinearModelling(wav, nt)
y = Lop * x
y = y + np.random.normal(0, 0e-3, x.shape)
In [7]:
sigma=0.01
l1 = VStackprox([L1(sigma=sigma), L2(sigma=sigma)], [60, nt-60])#L1(sigma=sigma)
l2 = L2(Op=Lop, b=y, niter=10)
Dop = FirstDerivative(nt, edge=True, kind='forward')
# Split-Bregman
mu = 1.3
lamda = 0.15
niter = 100
niterinner = 3
xsb = splitbregman(Lop, y, [Dop], x0=x0, niter_outer=niter, niter_inner=niterinner, mu=mu, epsRL1s=[lamda],
tol=1e-4, tau=1, **dict(iter_lim=30, damp=1e-3))[0]
# L-ADMM
#tau = 1. / np.real((Lop.H*Lop).eigs(neigs=1, which='LM')[0]) # optimal tau=1/maxeig(Lop^H Lop)
#mu = 0.99 * tau / np.real((Dop.H*Dop).eigs(neigs=1, which='LM')[0]) # optimal mu=tau/maxeig(Dop^H Dop)
L = 8 #np.real((Dop.H*Dop).eigs(neigs=1, which='LM')[0])
tau = 1.
mu = 0.99 * tau / L # optimal mu<=tau/maxeig(Dop^H Dop)
xladmm, zladmm = LinearizedADMM(l2, l1, Dop, tau=tau, mu=mu, x0=x0, niter=500)
fig, axs = plt.subplots(2, 1, figsize=(10, 7))
fig.suptitle('post-stack inversion', fontsize=14, fontweight='bold')
axs[0].plot(t, y, 'k', lw=2, label=r'$y$')
axs[1].plot(t, x, 'k', lw=2, label=r'$x$')
axs[1].plot(t, x0, 'k', lw=1, label=r'$x0$')
axs[1].plot(t, xladmm, 'b', lw=2, label=r'$x_{L-ADMM}$')
axs[1].plot(t, xsb, 'r', lw=2, label=r'$x_{Split-Bregman}$')
axs[1].legend();
In [8]:
sigma = 0.05
l1 = VStackprox([L1(), L2()], [60, nt-60]) #L1()
l2 = L2(Op=Lop, b=y, niter=20, warm=True)
Dop = FirstDerivative(nt, edge=True, kind='forward')
# approach with norms of 2 operators (CHECK!)
#tau = 1. / np.real((Lop.H*Lop).eigs(neigs=1, which='LM')[0]) # optimal tau=1/maxeig(Lop^H Lop)
#mu = 0.99 * tau / np.real((Dop.H*Dop).eigs(neigs=1, which='LM')[0]) # optimal mu=tau/maxeig(Dop^H Dop)
# most common approach
L = np.real((Dop.H*Dop).eigs(neigs=1, which='LM')[0])
tau = 1.
mu = 0.99 * tau / L # optimal mu<=tau/maxeig(Dop^H Dop)
xb = Bregman(l2, l1, x0, LinearizedADMM, A=Dop, alpha=sigma,
niterouter=10, show=True, **dict(tau=tau, mu=mu, niter=400))
fig, ax = plt.subplots(1, 1, figsize=(10, 3))
fig.suptitle('post-stack inversion', fontsize=14, fontweight='bold')
ax.plot(t, x, 'k', lw=2, label=r'$x$')
ax.plot(t, x0, 'k', lw=1, label=r'$x0$')
ax.plot(t, xladmm, 'b', lw=2, label=r'$x_{L-ADMM}$')
ax.plot(t, xb, 'r', lw=2, label=r'$x_{Bregman/L-ADMM}$')
ax.legend();
Bregman
---------------------------------------------------------
Proximal operator (f): <class 'pyproximal.proximal.L2.L2'>
Proximal operator (g): <class 'pyproximal.proximal.VStack.VStack'>
Linear operator (A): <class 'pylops.basicoperators.firstderivative.FirstDerivative'>
Inner Solver: <function LinearizedADMM at 0x7f784addac10>
alpha = 5.000000e-02 tolf = 1.000000e-10 tolx = 1.000000e-10
niter = 10
Itn x[0] f g J = f + g
1 7.50200e+00 7.573e-04 1.202e-02 1.277e-02
2 7.49609e+00 8.639e-07 1.350e-02 1.350e-02
3 7.49589e+00 4.446e-07 1.350e-02 1.350e-02
4 7.49583e+00 3.161e-07 1.350e-02 1.350e-02
5 7.49579e+00 2.548e-07 1.350e-02 1.350e-02
6 7.49576e+00 2.189e-07 1.350e-02 1.350e-02
7 7.49573e+00 1.950e-07 1.350e-02 1.350e-02
8 7.49570e+00 1.776e-07 1.350e-02 1.350e-02
9 7.49568e+00 1.642e-07 1.350e-02 1.350e-02
10 7.49565e+00 1.535e-07 1.350e-02 1.350e-02
Total time (s) = 6.71
---------------------------------------------------------
And PrimalDual solver
In [9]:
sigma=0.01
l1 = VStackprox([L1(sigma=sigma), L2(sigma=sigma)], [60, nt-60]) #L1(sigma=sigma)
l2 = L2(Op=Lop, b=y, niter=10)
Dop = FirstDerivative(nt, edge=True, kind='forward')
# most common approach
L = 8. # np.real((Dop.H*Dop).eigs(neigs=1, which='LM')[0])
#tau = 1.
#mu = 0.99 / (tau * L)
tau = 0.95 / np.sqrt(L)
mu = 0.95 / (tau * L)
xpd = PrimalDual(l2, l1, Dop, x0, tau=tau, mu=mu, niter=200, show=True)
xpd_ada, steps = AdaptivePrimalDual(l2, l1, Dop, x0, tau=tau, mu=mu, s=0.1,
niter=80, show=True)
fig, ax = plt.subplots(1, 1, figsize=(10, 3))
fig.suptitle('post-stack inversion', fontsize=14, fontweight='bold')
ax.plot(t, x, 'k', lw=2, label=r'$x$')
ax.plot(t, x0, 'k', lw=1, label=r'$x0$')
ax.plot(t, xladmm, 'b', lw=2, label=r'$x_{L-ADMM}$')
ax.plot(t, xpd, 'r', lw=2, label=r'$x_{Primal-Dual}$')
ax.plot(t, xpd_ada, 'g', lw=2, label=r'$x_{Ada-Primal-Dual}$')
ax.legend();
Primal-dual: min_x f(Ax) + x^T z + g(x)
---------------------------------------------------------
Proximal operator (f): <class 'pyproximal.proximal.L2.L2'>
Proximal operator (g): <class 'pyproximal.proximal.VStack.VStack'>
Linear operator (A): <class 'pylops.basicoperators.firstderivative.FirstDerivative'>
Additional vector (z): None
tau = 0.33587572106361 mu = 0.3535533905932738
theta = 1.00 niter = 200
Itn x[0] f g z^x J = f + g + z^x
1 7.49553e+00 1.370e-02 2.890e-03 0.000e+00 1.659e-02
2 7.49518e+00 5.442e-03 3.668e-03 0.000e+00 9.109e-03
3 7.49471e+00 3.333e-03 3.860e-03 0.000e+00 7.193e-03
4 7.49426e+00 2.504e-03 3.837e-03 0.000e+00 6.342e-03
5 7.49398e+00 2.078e-03 3.926e-03 0.000e+00 6.004e-03
6 7.49393e+00 1.816e-03 4.082e-03 0.000e+00 5.898e-03
7 7.49415e+00 1.636e-03 4.155e-03 0.000e+00 5.791e-03
8 7.49463e+00 1.502e-03 4.128e-03 0.000e+00 5.630e-03
9 7.49533e+00 1.397e-03 4.044e-03 0.000e+00 5.441e-03
10 7.49619e+00 1.312e-03 3.918e-03 0.000e+00 5.231e-03
21 7.50734e+00 7.678e-04 3.477e-03 0.000e+00 4.245e-03
41 7.52025e+00 3.957e-04 3.194e-03 0.000e+00 3.590e-03
61 7.52078e+00 2.487e-04 3.124e-03 0.000e+00 3.373e-03
81 7.51280e+00 1.883e-04 2.986e-03 0.000e+00 3.175e-03
101 7.50558e+00 1.608e-04 2.892e-03 0.000e+00 3.053e-03
121 7.50107e+00 1.406e-04 2.826e-03 0.000e+00 2.967e-03
141 7.49821e+00 1.237e-04 2.772e-03 0.000e+00 2.896e-03
161 7.49668e+00 1.078e-04 2.742e-03 0.000e+00 2.850e-03
181 7.49586e+00 8.930e-05 2.716e-03 0.000e+00 2.805e-03
192 7.49558e+00 7.918e-05 2.699e-03 0.000e+00 2.778e-03
193 7.49555e+00 7.831e-05 2.697e-03 0.000e+00 2.776e-03
194 7.49553e+00 7.745e-05 2.696e-03 0.000e+00 2.774e-03
195 7.49551e+00 7.659e-05 2.696e-03 0.000e+00 2.772e-03
196 7.49548e+00 7.574e-05 2.695e-03 0.000e+00 2.771e-03
197 7.49546e+00 7.490e-05 2.694e-03 0.000e+00 2.769e-03
198 7.49544e+00 7.407e-05 2.692e-03 0.000e+00 2.767e-03
199 7.49542e+00 7.325e-05 2.691e-03 0.000e+00 2.764e-03
200 7.49539e+00 7.244e-05 2.689e-03 0.000e+00 2.762e-03
Total time (s) = 0.37
---------------------------------------------------------
Adaptive Primal-dual: min_x f(Ax) + x^T z + g(x)
---------------------------------------------------------
Proximal operator (f): <class 'pyproximal.proximal.L2.L2'>
Proximal operator (g): <class 'pyproximal.proximal.VStack.VStack'>
Linear operator (A): <class 'pylops.basicoperators.firstderivative.FirstDerivative'>
Additional vector (z): None
tau0 = 3.358757e-01 mu0 = 3.535534e-01
alpha0 = 5.000000e-01 eta = 9.500000e-01
s = 1.000000e-01 delta = 1.500000e+00
niter = 80 tol = 1.000000e-10
Itn x[0] f g z^x J = f + g + z^x
1 7.49532e+00 1.370e-02 2.895e-03 0.000e+00 1.660e-02
2 7.49449e+00 3.815e-03 3.887e-03 0.000e+00 7.702e-03
3 7.49311e+00 1.940e-03 4.146e-03 0.000e+00 6.086e-03
4 7.49251e+00 1.351e-03 4.287e-03 0.000e+00 5.639e-03
5 7.49241e+00 9.318e-04 4.754e-03 0.000e+00 5.686e-03
6 7.49266e+00 7.886e-04 4.973e-03 0.000e+00 5.762e-03
7 7.49412e+00 7.429e-04 5.030e-03 0.000e+00 5.773e-03
8 7.49583e+00 6.301e-04 5.286e-03 0.000e+00 5.916e-03
9 7.49773e+00 5.929e-04 5.309e-03 0.000e+00 5.902e-03
16 7.51903e+00 2.539e-04 4.526e-03 0.000e+00 4.780e-03
24 7.52814e+00 1.279e-04 3.770e-03 0.000e+00 3.898e-03
32 7.53140e+00 1.110e-04 3.520e-03 0.000e+00 3.631e-03
40 7.53010e+00 9.756e-05 3.644e-03 0.000e+00 3.741e-03
48 7.52560e+00 8.656e-05 3.285e-03 0.000e+00 3.372e-03
56 7.51918e+00 8.027e-05 3.181e-03 0.000e+00 3.261e-03
64 7.51355e+00 7.259e-05 3.134e-03 0.000e+00 3.207e-03
71 7.50884e+00 6.571e-05 2.968e-03 0.000e+00 3.034e-03
72 7.50752e+00 6.477e-05 2.935e-03 0.000e+00 3.000e-03
73 7.50730e+00 6.442e-05 2.925e-03 0.000e+00 2.990e-03
74 7.50617e+00 6.421e-05 2.903e-03 0.000e+00 2.967e-03
75 7.50595e+00 6.326e-05 2.897e-03 0.000e+00 2.960e-03
76 7.50502e+00 6.319e-05 2.879e-03 0.000e+00 2.943e-03
77 7.50484e+00 6.236e-05 2.875e-03 0.000e+00 2.937e-03
78 7.50370e+00 6.219e-05 2.855e-03 0.000e+00 2.917e-03
79 7.50345e+00 6.152e-05 2.850e-03 0.000e+00 2.911e-03
80 7.50263e+00 6.080e-05 2.838e-03 0.000e+00 2.899e-03
Total time (s) = 0.15
In [10]:
plt.figure()
plt.plot(steps[0], 'k')
plt.plot(steps[1], 'r');
We can also choose the niter in L2 to be adaptive, start with fewer and then increase
In [11]:
def niter(count, niter0, k, nitermax):
niter = min(int(niter0 + k * count), nitermax)
return niter
plt.figure()
plt.plot(np.array([niter(i, 5, 0.1, 10) for i in np.arange(100)]), 'k');
In [12]:
sigma=0.01
l1 = VStackprox([L1(sigma=sigma), L2(sigma=sigma)], [60, nt-60]) #L1(sigma=sigma)
l2 = L2(Op=Lop, b=y, niter=lambda x: niter(x, 5, 0.1, 10))
Dop = FirstDerivative(nt, edge=True, kind='forward')
# most common approach
L = 8. # np.real((Dop.H*Dop).eigs(neigs=1, which='LM')[0])
tau = 1.
mu = 0.99 / (tau * L)
xpd_ada = PrimalDual(l2, l1, Dop, x0, tau=tau, mu=mu, niter=200, show=True)
fig, ax = plt.subplots(1, 1, figsize=(10, 3))
fig.suptitle('post-stack inversion', fontsize=14, fontweight='bold')
ax.plot(t, x, 'k', lw=2, label=r'$x$')
ax.plot(t, x0, 'k', lw=1, label=r'$x0$')
ax.plot(t, xpd, 'b', lw=2, label=r'$x_{Primal-Dual}$')
ax.plot(t, xpd_ada, 'r', lw=2, label=r'$x_{Primal-Dual}$ Ada')
ax.legend();
Primal-dual: min_x f(Ax) + x^T z + g(x)
---------------------------------------------------------
Proximal operator (f): <class 'pyproximal.proximal.L2.L2'>
Proximal operator (g): <class 'pyproximal.proximal.VStack.VStack'>
Linear operator (A): <class 'pylops.basicoperators.firstderivative.FirstDerivative'>
Additional vector (z): None
tau = 1.0 mu = 0.12375
theta = 1.00 niter = 200
Itn x[0] f g z^x J = f + g + z^x
1 7.48278e+00 2.534e-02 8.576e-03 0.000e+00 3.391e-02
2 7.48455e+00 3.822e-03 5.657e-03 0.000e+00 9.480e-03
3 7.48640e+00 1.859e-03 4.846e-03 0.000e+00 6.705e-03
4 7.48751e+00 1.576e-03 4.484e-03 0.000e+00 6.060e-03
5 7.48810e+00 1.470e-03 4.197e-03 0.000e+00 5.667e-03
6 7.48871e+00 1.306e-03 4.144e-03 0.000e+00 5.449e-03
7 7.49000e+00 1.133e-03 4.335e-03 0.000e+00 5.468e-03
8 7.49212e+00 1.018e-03 4.411e-03 0.000e+00 5.429e-03
9 7.49455e+00 9.055e-04 4.353e-03 0.000e+00 5.258e-03
10 7.49702e+00 8.201e-04 4.175e-03 0.000e+00 4.995e-03
21 7.51420e+00 3.169e-04 3.339e-03 0.000e+00 3.656e-03
41 7.52226e+00 1.428e-04 3.241e-03 0.000e+00 3.383e-03
61 7.50948e+00 8.959e-05 2.928e-03 0.000e+00 3.018e-03
81 7.49261e+00 7.542e-05 2.756e-03 0.000e+00 2.832e-03
101 7.48143e+00 7.064e-05 2.952e-03 0.000e+00 3.022e-03
121 7.47836e+00 6.040e-05 3.035e-03 0.000e+00 3.096e-03
141 7.48180e+00 4.583e-05 3.028e-03 0.000e+00 3.074e-03
161 7.48810e+00 3.359e-05 2.933e-03 0.000e+00 2.967e-03
181 7.49388e+00 2.572e-05 2.859e-03 0.000e+00 2.885e-03
192 7.49619e+00 2.389e-05 2.842e-03 0.000e+00 2.865e-03
193 7.49636e+00 2.382e-05 2.840e-03 0.000e+00 2.864e-03
194 7.49653e+00 2.377e-05 2.839e-03 0.000e+00 2.862e-03
195 7.49670e+00 2.374e-05 2.837e-03 0.000e+00 2.861e-03
196 7.49686e+00 2.373e-05 2.835e-03 0.000e+00 2.859e-03
197 7.49701e+00 2.373e-05 2.833e-03 0.000e+00 2.857e-03
198 7.49716e+00 2.376e-05 2.831e-03 0.000e+00 2.855e-03
199 7.49731e+00 2.379e-05 2.828e-03 0.000e+00 2.852e-03
200 7.49744e+00 2.385e-05 2.826e-03 0.000e+00 2.849e-03
Total time (s) = 0.36
---------------------------------------------------------
Approximated proximal of convolution operator¶
Now we use the fact that the modelling operator is a convolution to speed up the solution of the proximal operator.
First we ensure that we can rewrite our modelling operator as pure convolution and use that filter as input to our fast L2Convolve
In [13]:
derivative = np.array([0.5, 0, -0.5])
h = np.convolve(derivative, wav, mode='same')
plt.figure(figsize=(10, 3))
plt.plot(h, 'k')
# data
nfft = nt
y1 = np.fft.ifft(np.fft.fft(x, nfft) * np.fft.fft(h, nfft))
plt.figure(figsize=(10, 3))
plt.plot(y, 'k')
plt.plot(y1, '--r')
# LinearizedADMM inversion
sigma=0.01
l1 = VStackprox([L1(sigma=sigma), L2(sigma=sigma)], [50, nt-50]) #L1(sigma=sigma)
l2 = L2Convolve(h=h, b=y, nfft=nt)
Dop = FirstDerivative(nt, edge=True, kind='forward')
L = np.real((Dop.H*Dop).eigs(neigs=1, which='LM')[0])
tau = 1.
mu = 0.99 * tau / L # optimal mu<=tau/maxeig(Dop^H Dop)
xladmm, zladmm = LinearizedADMM(l2, l1, Dop, tau=tau, mu=mu, x0=x0, niter=1000, show=True)
# Bregman inversion
sigma=0.1
l1 = VStackprox([L1(), L2()], [50, nt-50]) #L1()
xb = Bregman(l2, l1, x0, LinearizedADMM, A=Dop, niterouter=5, alpha=sigma,
show=True, **dict(tau=tau, mu=mu, niter=400))
# recenter signal
xladmm = np.pad(xladmm[:-wavc], (wavc, 0), constant_values=x0[0])
xb = np.pad(xb[:-wavc], (wavc, 0), constant_values=x0[0])
fig, ax = plt.subplots(1, 1, figsize=(10, 3))
fig.suptitle('post-stack inversion', fontsize=14, fontweight='bold')
ax.plot(t, x, 'k', lw=2, label=r'$x$')
ax.plot(t, x0, 'k', lw=1, label=r'$x0$')
ax.plot(t, xladmm, 'b', lw=2, label=r'$x_{L-ADMM}$ Fastconv')
ax.plot(t, xb, 'r', lw=2, label=r'$x_{Bregman}$ Fastconv')
ax.legend();
Linearized-ADMM
---------------------------------------------------------
Proximal operator (f): <class 'pyproximal.proximal.L2.L2Convolve'>
Proximal operator (g): <class 'pyproximal.proximal.VStack.VStack'>
Linear operator (A): <class 'pylops.basicoperators.firstderivative.FirstDerivative'>
tau = 1.000000e+00 mu = 2.475599e-01 niter = 1000
Itn x[0] f g J = f + g
1 7.49707e+00 1.883e-02 2.467e-03 2.130e-02
2 7.49919e+00 8.626e-03 3.115e-03 1.174e-02
3 7.50121e+00 5.155e-03 3.226e-03 8.381e-03
4 7.50285e+00 3.685e-03 3.569e-03 7.254e-03
5 7.50404e+00 2.952e-03 3.738e-03 6.689e-03
6 7.50481e+00 2.527e-03 3.784e-03 6.311e-03
7 7.50522e+00 2.259e-03 3.774e-03 6.033e-03
8 7.50537e+00 2.077e-03 3.706e-03 5.784e-03
9 7.50536e+00 1.943e-03 3.645e-03 5.589e-03
10 7.50531e+00 1.836e-03 3.642e-03 5.477e-03
101 7.49348e+00 2.868e-04 3.038e-03 3.325e-03
201 7.48695e+00 1.295e-04 2.814e-03 2.944e-03
301 7.48549e+00 7.956e-05 2.747e-03 2.826e-03
401 7.48932e+00 6.839e-05 2.709e-03 2.777e-03
501 7.49089e+00 6.954e-05 2.679e-03 2.749e-03
601 7.49261e+00 7.034e-05 2.662e-03 2.732e-03
701 7.49384e+00 6.761e-05 2.651e-03 2.718e-03
801 7.49488e+00 6.629e-05 2.642e-03 2.709e-03
901 7.49583e+00 6.633e-05 2.635e-03 2.702e-03
992 7.49663e+00 6.682e-05 2.630e-03 2.697e-03
993 7.49664e+00 6.682e-05 2.630e-03 2.697e-03
994 7.49665e+00 6.683e-05 2.630e-03 2.697e-03
995 7.49666e+00 6.683e-05 2.630e-03 2.697e-03
996 7.49666e+00 6.684e-05 2.630e-03 2.697e-03
997 7.49667e+00 6.685e-05 2.630e-03 2.697e-03
998 7.49668e+00 6.685e-05 2.630e-03 2.697e-03
999 7.49669e+00 6.686e-05 2.630e-03 2.697e-03
1000 7.49670e+00 6.686e-05 2.630e-03 2.697e-03
Total time (s) = 0.09
---------------------------------------------------------
Bregman
---------------------------------------------------------
Proximal operator (f): <class 'pyproximal.proximal.L2.L2Convolve'>
Proximal operator (g): <class 'pyproximal.proximal.VStack.VStack'>
Linear operator (A): <class 'pylops.basicoperators.firstderivative.FirstDerivative'>
Inner Solver: <function LinearizedADMM at 0x7f784addac10>
alpha = 1.000000e-01 tolf = 1.000000e-10 tolx = 1.000000e-10
niter = 5
Itn x[0] f g J = f + g
1 7.51210e+00 2.982e-03 2.094e-02 2.392e-02
2 7.49974e+00 3.810e-05 2.680e-02 2.684e-02
3 7.49933e+00 3.687e-05 2.683e-02 2.686e-02
4 7.49903e+00 3.640e-05 2.685e-02 2.689e-02
5 7.49882e+00 3.618e-05 2.687e-02 2.691e-02
Total time (s) = 0.19
---------------------------------------------------------
This seems to work under the following conditions:
- we need to assume that our filter has a certain phase (not centered to zero)... and then recenter things ourselves
- the model at the top and base should have same AI, otherwise the modelling artificially creates a discontinuity and a new event in the data... if not the case, we need pad the base with an artificial gradient
- we need to separate the regularization between the model part and the gradient