# PyLops and OccamyPy together

Authors
Affiliations
Francesco Picetti
Image and Sound Processing Lab, Politecnico di Milano
Ettore Biondi
SeismoLab, California Institute of Technology

# Using PyLops and OccamyPy together¶

@Author: Francesco Picetti - picettifrancesco@gmail.com

In this notebook, we will combine PyLops and OccamyPy.

First, let me express my gratitude to Matteo Ravasi; PyLops has strongly inspired us to publish OccamyPy, as our contribution to the scientific community.

In the following, we leverage PyLops rich library of operators to be used with OccamyPy, which we found to be agile for distributed computation and nonlinear problems.

import numpy as np
import occamypy as o
import matplotlib.pyplot as plt
import pylops

## Mechanism of the interface¶

Let's say we have a VectorNumpy

x = o.VectorNumpy((3,5)).set(1.)

and we want to use a PyLops operator like

Dp = pylops.Diagonal(np.arange(x.size))

we can cast such operator to OccamyPy with FromPylops

D = o.pylops_interface.FromPylops(domain=x, range=x, op=Dp)

Now we have a OccamyPy operator that can be used in our inversion

D.dotTest(True)

At the same time, we can cast a Numpy-based OccamyPy operator to PyLops with ToPylops

S = o.Scaling(x, 2.)
Sp = o.pylops_interface.ToPylops(S)

Such operator can be used inside PyLops

pylops.utils.dottest(Sp, x.size, x.size, tol=1e-6, complexflag=0, verb=True)

print("\t", np.isclose((Sp * x.arr.ravel()).reshape(x.shape), 2. * x.arr).all())

## Example 1: CT Scan Imaging¶

Let's replicate this example. The modeling equation is $\mathbf{d} = \mathbf{R} \, \mathbf{m}$

def radoncurve(x, r, theta):
curve = (r - ny // 2) / (np.sin(np.deg2rad(theta)) + 1e-15) \
+ np.tan(np.deg2rad(90 - theta)) * x + ny // 2
return curve

Load the image and create the 2D Radon transform operator as done in PyLops

x = np.load('data/shepp_logan_phantom.npy').T
x = x / x.max()
nx, ny = x.shape

ntheta = 150
theta = np.linspace(0., 180., ntheta, endpoint=False)

centeredh=True, interp=False,
engine='numpy', dtype='float64')

Let's see the operator at work

y = RLop.H * x.ravel()
y = y.reshape(ntheta, ny)

xrec = RLop * y.ravel()
xrec = xrec.reshape(nx, ny)

fig, axs = plt.subplots(1, 3, figsize=(10, 4))
axs[0].imshow(x.T, vmin=0, vmax=1, cmap='gray')
axs[0].set_title('Model')
axs[0].axis('tight')
axs[1].imshow(y.T, cmap='gray')
axs[1].set_title('Data')
axs[1].axis('tight')
axs[2].imshow(xrec.T, cmap='gray')
axs[2].axis('tight')
fig.tight_layout()
plt.show()

Now let's build the OccamyPy equivalent

x_ = o.VectorNumpy(x)
y_ = o.VectorNumpy((ntheta, ny))
R_ = o.pylops_interface.FromPylops(x_, y_, RLop.H)
y_ = R_ * x_
xrec_ = R_.H * y_
fig, axs = plt.subplots(1, 3, figsize=(10, 4))
axs[0].imshow(x_.plot().T, vmin=0, vmax=1, cmap='gray')
axs[0].set_title('Model')
axs[0].axis('tight')
axs[1].imshow(y_.plot().T, cmap='gray')
axs[1].set_title('Data')
axs[1].axis('tight')
axs[2].imshow(xrec_.plot().T, cmap='gray')
axs[2].axis('tight')
fig.tight_layout()
plt.show()

Set up the TV-regularized inversion as done in PyLops

mu = 1.5
lamda = 1.
niter = 3
niterinner = 4

problemTV = o.GeneralizedLasso(x_.clone().zero(), y_, R_, eps=lamda / mu,

SB = o.SplitBregman(o.BasicStopper(niter=niter),
niter_inner=niterinner,
niter_solver=20,
linear_solver="LSQR",
warm_start=True)

SB.setDefaults(save_obj=True)
SB.run(problemTV, verbose=True, inner_verbose=False)
fig, axs = plt.subplots(1, 2, figsize=(7, 4))
axs[0].imshow(x_.plot().T, vmin=0, vmax=1, cmap='gray')
axs[0].set_title('Model')
axs[0].axis('tight')
axs[1].imshow(problemTV.model.plot().T, vmin=0, vmax=1, cmap='gray')
axs[1].set_title('Inverted domain')
axs[1].axis('tight')
fig.tight_layout()
plt.show()