Skip to content

Commit 2e9fdf6

Browse files
mrava87mrava87
committed
Added Convolve1D, FirstDerivative, SecondDerivative, Laplacian, and Smoothing1D
1 parent 7ef1be1 commit 2e9fdf6

File tree

17 files changed

+1112
-46
lines changed

17 files changed

+1112
-46
lines changed

‎docs/source/api/index.rst‎

Lines changed: 12 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -27,6 +27,18 @@ Basic operators
2727
BlockDiag
2828

2929

30+
Smoothing and derivatives
31+
~~~~~~~~~~~~~~~~~~~~~~~~~
32+
33+
.. autosummary::
34+
:toctree: generated/
35+
36+
Smoothing1D
37+
FirstDerivative
38+
SecondDerivative
39+
Laplacian
40+
41+
3042
Signal processing
3143
~~~~~~~~~~~~~~~~~
3244

‎examples/plot_convolve.py‎

Lines changed: 79 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,79 @@
1+
"""
2+
Convolution
3+
===========
4+
This example shows how to use the
5+
:py:class:`pylops_distributed.signalprocessing.Convolve1D` operator to perform
6+
convolution between two signals.
7+
8+
Note that when the input dataset is distributed across multiple nodes,
9+
additional care should be taken when applying convolution.
10+
In PyLops-distributed, we leverage the :func:`dask.array.map_block` and
11+
:func:`dask.array.map_overlap` functionalities of dask when chunking is
12+
performed along the dimension of application of the convolution and along any
13+
other dimension, respectively. This is however handled by the inner working of
14+
:py:class:`pylops_distributed.signalprocessing.Convolve1D`.
15+
16+
"""
17+
import numpy as np
18+
import dask.array as da
19+
import matplotlib.pyplot as plt
20+
from scipy.sparse.linalg import lsqr
21+
22+
import pylops_distributed
23+
from pylops.utils.wavelets import ricker
24+
25+
plt.close('all')
26+
27+
###############################################################################
28+
# We will start by creating a zero signal of lenght :math:`nt` and we will
29+
# place a unitary spike at its center. We also create our filter to be
30+
# applied by means of
31+
# :py:class:`pylops-distributed.signalprocessing.Convolve1D` operator.
32+
# Following the seismic example mentioned above, the filter is a
33+
# `Ricker wavelet <http://subsurfwiki.org/wiki/Ricker_wavelet>`_
34+
# with dominant frequency :math:`f_0 = 30 Hz`.
35+
nt = 1001
36+
dt = 0.004
37+
t = np.arange(nt)*dt
38+
x = np.zeros(nt)
39+
x[int(nt/2)] = 1
40+
x = da.from_array(x, chunks=nt // 2 + 1)
41+
h, th, hcenter = ricker(t[:101], f0=30)
42+
43+
Cop = pylops_distributed.signalprocessing.Convolve1D(nt, h=h, offset=hcenter,
44+
chunks=(nt // 2 + 1,
45+
nt // 2 + 1),
46+
dtype='float32')
47+
y = Cop * x
48+
xinv = Cop / y
49+
50+
fig, ax = plt.subplots(1, 1, figsize=(10, 3))
51+
ax.plot(t, x, 'k', lw=2, label=r'$x$')
52+
ax.plot(t, y, 'r', lw=2, label=r'$y=Ax$')
53+
ax.plot(t, xinv, '--g', lw=2, label=r'$x_{ext}$')
54+
ax.set_title('Convolve 1d data', fontsize=14, fontweight='bold')
55+
ax.legend()
56+
ax.set_xlim(1.9, 2.1)
57+
58+
###############################################################################
59+
# We show now that also a filter with mixed phase (i.e., not centered
60+
# around zero) can be applied and inverted for using the
61+
# :py:class:`pylops.signalprocessing.Convolve1D`
62+
# operator.
63+
Cop = pylops_distributed.signalprocessing.Convolve1D(nt, h=h, offset=hcenter - 3,
64+
chunks=(nt // 2 + 1,
65+
nt // 2 + 1),
66+
dtype='float32')
67+
y = Cop * x
68+
y1 = Cop.H * x
69+
xinv = Cop / y
70+
71+
fig, ax = plt.subplots(1, 1, figsize=(10, 3))
72+
ax.plot(t, x, 'k', lw=2, label=r'$x$')
73+
ax.plot(t, y, 'r', lw=2, label=r'$y=Ax$')
74+
ax.plot(t, y1, 'b', lw=2, label=r'$y=A^Hx$')
75+
ax.plot(t, xinv, '--g', lw=2, label=r'$x_{ext}$')
76+
ax.set_title('Convolve 1d data with non-zero phase filter', fontsize=14,
77+
fontweight='bold')
78+
ax.set_xlim(1.9, 2.1)
79+
ax.legend()

‎examples/plot_derivative.py‎

Lines changed: 184 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,184 @@
1+
"""
2+
Derivatives
3+
===========
4+
This example shows how to use the suite of derivative operators, namely
5+
:py:class:`pylops_distributed.FirstDerivative`,
6+
:py:class:`pylops_distributed.SecondDerivative`, and
7+
:py:class:`pylops_distributed.Laplacian`.
8+
9+
The derivative operators can be applied along any dimension of a N-dimensional
10+
input. When the input is chuncked along any other direction(s) than the one
11+
the derivative is applied, the derivative is efficiently performed without
12+
neither communication between workers nor duplication of part of the input
13+
array. On the other hand when the input is chuncked along the direction where
14+
the derivative is applied, the chunks are partially overlapping such that no
15+
communication is required between the workers when applying the derivative.
16+
17+
In some applications, the user cannot avoid this second scenario to happen
18+
(e.g, when the derivative should be applied to all the dimensions of the
19+
dataset). Nevertheless our implementation makes this possible in a fully
20+
transparent way and with very little additional overhead.
21+
22+
"""
23+
import numpy as np
24+
import dask.array as da
25+
import matplotlib.pyplot as plt
26+
27+
import pylops
28+
import pylops_distributed
29+
30+
plt.close('all')
31+
np.random.seed(0)
32+
33+
###############################################################################
34+
# Let's start by looking at a simple first-order centered derivative. We
35+
# chunck the vector in 3 chunks.
36+
nx = 100
37+
nchunks = 3
38+
39+
x = np.zeros(nx)
40+
x[int(nx/2)] = 1
41+
xd = da.from_array(x, chunks=nx // nchunks + 1)
42+
print('x:', xd)
43+
44+
Dop = pylops.FirstDerivative(nx, dtype='float32')
45+
dDop = pylops_distributed.FirstDerivative(nx, compute=(True, True),
46+
dtype='float32')
47+
48+
y = Dop * x
49+
xadj = Dop.H * y
50+
51+
yd = Dop * xd
52+
xadjd = Dop.H * yd
53+
54+
55+
fig, axs = plt.subplots(3, 1, figsize=(13, 8))
56+
axs[0].stem(np.arange(nx), x, linefmt='k', markerfmt='ko',
57+
use_line_collection=True)
58+
axs[0].set_title('Input', size=20, fontweight='bold')
59+
axs[1].stem(np.arange(nx), y, linefmt='--r', markerfmt='ro',
60+
label='Numpy', use_line_collection=True)
61+
axs[1].stem(np.arange(nx), yd, linefmt='--r', markerfmt='ro',
62+
label='Dask', use_line_collection=True)
63+
axs[1].set_title('Forward', size=20, fontweight='bold')
64+
axs[1].legend()
65+
axs[2].stem(np.arange(nx), xadj, linefmt='k', markerfmt='ko',
66+
label='Numpy', use_line_collection=True)
67+
axs[2].stem(np.arange(nx), xadjd, linefmt='--r', markerfmt='ro',
68+
label='Dask', use_line_collection=True)
69+
axs[2].set_title('Adjoint', size=20, fontweight='bold')
70+
axs[2].legend()
71+
plt.tight_layout()
72+
73+
###############################################################################
74+
# As expected we obtain the same result, with the only difference that
75+
# in the second case we did not need to explicitly create a matrix,
76+
# saving memory and computational time.
77+
#
78+
# Let's move onto applying the same first derivative to a 2d array in
79+
# the first direction. Now we consider two cases, first when the data is
80+
# chunked along the first direction and second when its chunked along the
81+
# second direction.
82+
nx, ny = 11, 21
83+
nchunks = 2
84+
85+
A = np.zeros((nx, ny))
86+
A[nx//2, ny//2] = 1.
87+
A1d = da.from_array(A, chunks= (nx // nchunks + 1, ny))
88+
A2d = da.from_array(A, chunks= (nx , ny // nchunks + 1))
89+
print('A1d:', A1d)
90+
print('A2d:', A2d)
91+
92+
Dop = pylops_distributed.FirstDerivative(nx * ny, dims=(nx, ny),
93+
compute=(True, True),
94+
dir=0, dtype='float64')
95+
96+
B1d = np.reshape(Dop * A1d.flatten(), (nx, ny))
97+
B2d = np.reshape(Dop * A2d.flatten(), (nx, ny))
98+
99+
fig, axs = plt.subplots(1, 3, figsize=(12, 3))
100+
fig.suptitle('First Derivative in 1st direction', fontsize=12,
101+
fontweight='bold', y=0.95)
102+
im = axs[0].imshow(A, interpolation='nearest', cmap='rainbow')
103+
axs[0].axis('tight')
104+
axs[0].set_title('x')
105+
plt.colorbar(im, ax=axs[0])
106+
im = axs[1].imshow(B1d, interpolation='nearest', cmap='rainbow')
107+
axs[1].axis('tight')
108+
axs[1].set_title('y (1st dir chunks)')
109+
plt.colorbar(im, ax=axs[1])
110+
im = axs[2].imshow(B2d, interpolation='nearest', cmap='rainbow')
111+
axs[2].axis('tight')
112+
axs[2].set_title('y (2nd dir chunks)')
113+
plt.colorbar(im, ax=axs[2])
114+
plt.tight_layout()
115+
plt.subplots_adjust(top=0.8)
116+
117+
###############################################################################
118+
# We can now do the same for the second derivative
119+
A = np.zeros((nx, ny))
120+
A[nx//2, ny//2] = 1.
121+
A1d = da.from_array(A, chunks= (nx // nchunks + 1, ny))
122+
A2d = da.from_array(A, chunks= (nx , ny // nchunks + 1))
123+
print('A1d:', A1d)
124+
print('A2d:', A2d)
125+
126+
Dop = pylops_distributed.SecondDerivative(nx * ny, dims=(nx, ny),
127+
compute=(True, True),
128+
dir=0, dtype='float64')
129+
130+
B1d = np.reshape(Dop * A1d.flatten(), (nx, ny))
131+
B2d = np.reshape(Dop * A2d.flatten(), (nx, ny))
132+
133+
fig, axs = plt.subplots(1, 3, figsize=(12, 3))
134+
fig.suptitle('Second Derivative in 1st direction', fontsize=12,
135+
fontweight='bold', y=0.95)
136+
im = axs[0].imshow(A, interpolation='nearest', cmap='rainbow')
137+
axs[0].axis('tight')
138+
axs[0].set_title('x')
139+
plt.colorbar(im, ax=axs[0])
140+
im = axs[1].imshow(B1d, interpolation='nearest', cmap='rainbow')
141+
axs[1].axis('tight')
142+
axs[1].set_title('y (1st dir chunks)')
143+
plt.colorbar(im, ax=axs[1])
144+
im = axs[2].imshow(B2d, interpolation='nearest', cmap='rainbow')
145+
axs[2].axis('tight')
146+
axs[2].set_title('y (2nd dir chunks)')
147+
plt.colorbar(im, ax=axs[2])
148+
plt.tight_layout()
149+
plt.subplots_adjust(top=0.8)
150+
151+
152+
###############################################################################
153+
# We use the symmetrical Laplacian operator as well
154+
# as a asymmetrical version of it (by adding more weight to the
155+
# derivative along one direction)
156+
157+
# symmetrical
158+
L2symop = pylops_distributed.Laplacian(dims=(nx, ny), weights=(1, 1),
159+
compute=(True, True), dtype='float64')
160+
161+
# asymmetrical
162+
L2asymop = pylops_distributed.Laplacian(dims=(nx, ny), weights=(3, 1),
163+
compute=(True, True), dtype='float64')
164+
165+
Bsym = np.reshape(L2symop * A1d.flatten(), (nx, ny))
166+
Basym = np.reshape(L2asymop * A2d.flatten(), (nx, ny))
167+
168+
fig, axs = plt.subplots(1, 3, figsize=(10, 3))
169+
fig.suptitle('Laplacian', fontsize=12,
170+
fontweight='bold', y=0.95)
171+
im = axs[0].imshow(A, interpolation='nearest', cmap='rainbow')
172+
axs[0].axis('tight')
173+
axs[0].set_title('x')
174+
plt.colorbar(im, ax=axs[0])
175+
im = axs[1].imshow(Bsym, interpolation='nearest', cmap='rainbow')
176+
axs[1].axis('tight')
177+
axs[1].set_title('y sym')
178+
plt.colorbar(im, ax=axs[1])
179+
im = axs[2].imshow(Basym, interpolation='nearest', cmap='rainbow')
180+
axs[2].axis('tight')
181+
axs[2].set_title('y asym')
182+
plt.colorbar(im, ax=axs[2])
183+
plt.tight_layout()
184+
plt.subplots_adjust(top=0.8)

‎examples/plot_smoothing1d.py‎

Lines changed: 92 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,92 @@
1+
r"""
2+
1D Smoothing
3+
============
4+
5+
This example shows how to use the :py:class:`pylops_distributed.Smoothing1D`
6+
operator to smooth an input signal along a given axis.
7+
8+
A smoothing operator is a simple compact filter on lenght :math:`n_{smooth}`
9+
and each elements is equal to :math:`1/n_{smooth}`.
10+
"""
11+
12+
import numpy as np
13+
import dask.array as da
14+
import matplotlib.pyplot as plt
15+
16+
import pylops_distributed
17+
18+
plt.close('all')
19+
20+
###############################################################################
21+
# Define the input parameters: number of samples of input signal (``N``) and
22+
# lenght of the smoothing filter regression coefficients (:math:`n_{smooth}`).
23+
# In this first case the input signal is one at the center and zero elsewhere.
24+
N = 31
25+
nchunks = 3
26+
nsmooth = 7
27+
x = np.zeros(N)
28+
x[int(N/2)] = 1
29+
x = da.from_array(x, chunks=N // nchunks + 1)
30+
print('x:', x)
31+
32+
Sop = pylops_distributed.Smoothing1D(nsmooth=nsmooth, dims=[N],
33+
dtype='float32')
34+
35+
y = Sop * x
36+
xadj = Sop.H * y
37+
38+
fig, ax = plt.subplots(1, 1, figsize=(10, 3))
39+
ax.plot(x, 'k', lw=2, label=r'$x$')
40+
ax.plot(y, 'r', lw=2, label=r'$y=Ax$')
41+
ax.set_title('Smoothing in 1st direction', fontsize=14, fontweight='bold')
42+
ax.legend()
43+
44+
###############################################################################
45+
# Let's repeat the same exercise with a random signal as input. After applying
46+
# smoothing, we will also try to invert it.
47+
N = 120
48+
nchunks = 3
49+
nsmooth = 13
50+
x = np.random.normal(0, 1, N)
51+
x = da.from_array(x, chunks=N // nchunks + 1)
52+
print('x:', x)
53+
Sop = pylops_distributed.Smoothing1D(nsmooth=13, dims=(N), dtype='float32')
54+
55+
y = Sop * x
56+
xest = Sop / y
57+
58+
fig, ax = plt.subplots(1, 1, figsize=(10, 3))
59+
ax.plot(x, 'k', lw=2, label=r'$x$')
60+
ax.plot(y, 'r', lw=2, label=r'$y=Ax$')
61+
ax.plot(xest, '--g', lw=2, label=r'$x_{ext}$')
62+
ax.set_title('Smoothing in 1st direction',
63+
fontsize=14, fontweight='bold')
64+
ax.legend()
65+
66+
###############################################################################
67+
# Finally we show that the same operator can be applied to multi-dimensional
68+
# data along a chosen axis.
69+
nx, ny = 21, 31
70+
nchunks = 2, 3
71+
A = np.zeros((nx, ny))
72+
A[5, 10] = 1
73+
A = da.from_array(A, chunks= (nx // nchunks[0] + 1, ny // nchunks[1] + 1))
74+
print('A:', A)
75+
76+
Sop = pylops_distributed.Smoothing1D(nsmooth=5, dims=(nx, ny),
77+
dir=0, dtype='float64')
78+
B = da.reshape(Sop * A.ravel(), (nx, ny))
79+
80+
fig, axs = plt.subplots(1, 2, figsize=(10, 3))
81+
fig.suptitle('Smoothing in 1st direction for 2d data', fontsize=14,
82+
fontweight='bold', y=0.95)
83+
im = axs[0].imshow(A, interpolation='nearest', vmin=0, vmax=1)
84+
axs[0].axis('tight')
85+
axs[0].set_title('Model')
86+
plt.colorbar(im, ax=axs[0])
87+
im = axs[1].imshow(B, interpolation='nearest', vmin=0, vmax=1)
88+
axs[1].axis('tight')
89+
axs[1].set_title('Data')
90+
plt.colorbar(im, ax=axs[1])
91+
plt.tight_layout()
92+
plt.subplots_adjust(top=0.8)

‎pylops_distributed/__init__.py‎

Lines changed: 4 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -10,6 +10,10 @@
1010
from .basicoperators import HStack
1111
from .basicoperators import Block
1212
from .basicoperators import BlockDiag
13+
from .basicoperators import FirstDerivative
14+
from .basicoperators import SecondDerivative
15+
from .basicoperators import Laplacian
16+
from .basicoperators import Smoothing1D
1317

1418
from .optimization.cg import cg, cgls
1519

0 commit comments

Comments
 (0)