The bot created this issue to inform you that pyup.io has been set up on this repo.
Once you have closed it, the bot will open pull requests for updates as soon as they are available.

Take advantage of broadcasting ability of numpy.matmul as a way to speed up integral computation (i.e., inner for loop) in MDC operator.

This requires a slight change in the code as frequency need to move to first dimension. It will make fft and ifft less optimal but as the for loop is bottleneck this should give an overall speed up to the operator.

The "Improve this page" links on read the docs are dead.

Example:

Leads to:

So it links to
https://github.com/pylops/edit/master/docs/installation.rst

instead to
https://github.com/Statoil/pylops/edit/master/docs/source/index.rst
or simply
https://github.com/Statoil/pylops/blob/master/docs/source/index.rst

I tried to find the source of the problem in the conf.py and the _templates/breadcrumbs.html files, but I couldn't locate the problem.

Check first scipy version when oaconvolve was introduced and force it as minimum requirement for pylops

Create NMO operator performing NMO modelling like in http://sepwww.stanford.edu/sep/prof/bei/vela/paper_html/node6.html

Change

Move from function based to class based solvers. This would require creating a base class with the following:

• __init__
• step: run a single step
• callback: run after any step (a user can just overwrite it - currently is passed to the solver as function)
• print: print update
• solve: run entire solver

Motivation

• This allows running multiple inversion with different solver parametrization (niter, tol, etc) without having to rerun the preparation part, which could be time consuming and independent on the solver parameters.
• This allows users to run both the solver as is or run it in a more fine-grained fashion (using step).

Countermotivation

• This is a breaking change, is it worth?
• Scipy has function based solvers, if we move away from them (for also cg, cgls, lsqr) we will need to call them differently inside our apps depending on whether we use numpy or cupy arrays - right now we just call a different solver but with a very similar function signature.

Hi, thank you guys for providing these good algorithms for us to do research.

I am a student majoring in geophysics, I am using your Marchenko redatuming python code. After I changed the virtual source locations, the Green's function the code produced is incorrect.

I think the parameter trav has changed with the virtual source location, but the function MarchenkoWM.apply_onepoint may need to fix. Could you please check and fix it?

Thank you for your help. I am looking forward to your reply.
Stay safe!
Mengli

Is there a way to convert a LinearOperator to an explicit sparse matrix representation.

So a sparse equivalent of the todense() method that returns a scipy.sparse.csr_matrix instead of a numpy array, for example?

If not, do you have a suggestion for a workaround? One possibility is to apply the linear operator to each of the (dense) unit vectors and record the indices and entries of the non zeros of the results. Is there a more elegant solution?

Appreciate the help!

matplotlib is used only within utils.wavelets. Consider removing those plots and free the user version from matplotlib dependency

checkimag seems to be the same as checkreal, but I suspect that the np.reals should be changed to np.imags:

Add code to handle edges where stencil goes out of available signal support and additional input parameter edges=True/False to choose whether to use reduce order derivatives at edges or just ignore (as per today)

Code is actually present and commented out in SecondDerivative at the moment. It is required to:

• validate that both options pass dot test
• validate that the operator can be inverted for both options

Modify all return statement patterns:

if: 
    return ... 
else: 
   return ...

to always return once (and the same number of elements - this change can only be implemented in v2.0.0)

Not exactly a bug, but it should be avoided np.flatten() unless needed (both codebase and documentation examples).

By using Gradient Operator for Total Variation denosing in my project I got some artifacts in the solution. If I define my own sparse version it is fine.
You are using central differences in Gradient/FirstDerivative. With wikipedia 'oscillating functions can yield zero derivative', a working Gradient with backward differences and my desired solution of piecewise constant TV/BV functions, I think of central differences as my problem.

Is there a way of setting the difference scheme? Or would you implement a parameter for this in pylops future versions?

The example on CT scan imaging generates an exception because the pylops.signalprocessing.Radon2D does not appear to accept for the kind argument a function object.

Using version 1.8.0 of PyLops under Python 3.7.6.

---------------------------------------------------------------------------
NotImplementedError                       Traceback (most recent call last)
<ipython-input-506-a87c57c49272> in <module>
     17                                     theta, kind=radoncurve,
     18                                     centeredh=True, interp=False,
---> 19                                     engine='numba', dtype='float64')
     20 
     21 y = RLop.H * x.T.ravel()

/usr/local/lib/python3.7/site-packages/pylops/signalprocessing/Radon2D.py in Radon2D(taxis, haxis, pxaxis, kind, centeredh, interp, onthefly, engine, dtype)
    179         f = _hyperbolic if engine == 'numpy' else _hyperbolic_numba
    180     else:
--> 181         raise NotImplementedError('kind must be linear, '
    182                                   'parabolic, or hyperbolic...')
    183     # make axes unitless

NotImplementedError: kind must be linear, parabolic, or hyperbolic...

Package versions on my system (OS X Catalina)

Package                            Version
---------------------------------- --------------------
absl-py                            0.7.0
alabaster                          0.7.11
algopy                             0.5.7
altgraph                           0.16.1
amqp                               2.5.0
appdirs                            1.4.3
appnope                            0.1.0
astor                              0.7.1
astroid                            1.6.5
Babel                              2.6.0
backports-abc                      0.5
backports.functools-lru-cache      1.5
backports.shutil-get-terminal-size 1.0.0
backports.weakref                  1.0.post1
beautifulsoup4                     4.6.1
billiard                           3.6.0.0
bleach                             1.5.0
bs4                                0.0.1
bspline                            0.1.1
bullseye                           0.1.1
Cartopy                            0.17.0
celery                             4.3.0
certifi                            2018.4.16
chardet                            3.0.4
click                              6.7
cloudpickle                        0.5.3
CommonMark                         0.7.5
configparser                       3.5.0
corner                             2.0.1
CVXcanon                           0.1.1
cvxopt                             1.2.0
cvxpy                              1.0.6
cycler                             0.10.0
cypy                               0.2.0
Cython                             0.29.5
dask                               0.18.2
decon                              0.1.7
decorator                          4.3.0
dicom                              0.9.9.post1
dill                               0.2.8.2
Django                             1.11.20
docutils                           0.14
ecos                               2.0.5
emcee                              2.2.1
entrypoints                        0.2.3
enum34                             1.1.6
et-xmlfile                         1.0.1
fastcache                          1.0.2
Flask                              1.1.1
Flask-Login                        0.4.1
Flask-Redis                        0.3.0
Flask-SSE                          0.2.1
Flask-WTF                          0.14.2
funcsigs                           1.0.2
functools32                        3.2.3.post2
future                             0.16.0
futures                            3.2.0
gast                               0.2.2
geos                               0.2.1
gevent                             1.3.5
google-pasta                       0.1.7
greenlet                           0.4.14
grpcio                             1.18.0
gunicorn                           19.9.0
h5py                               2.8.0
html5lib                           0.9999999
idna                               2.6
imagesize                          1.0.0
ipaddress                          1.0.22
ipykernel                          4.9.0
ipython                            5.8.0
ipython-genutils                   0.2.0
ipywidgets                         7.4.0
ir-fits                            0.6.4
isort                              4.3.4
itsdangerous                       0.24
jdcal                              1.4
jedi                               0.12.1
Jinja2                             2.10.3
joblib                             0.12.4
jsonschema                         2.6.0
jupyter                            1.0.0
jupyter-client                     5.2.3
jupyter-console                    5.2.0
jupyter-core                       4.4.0
Keras-Applications                 1.0.8
Keras-Preprocessing                1.0.6
keyring                            13.2.1
kiwisolver                         1.0.1
kombu                              4.6.0
latest                             0.4.1
lazy-object-proxy                  1.3.1
lxml                               4.2.3
macholib                           1.10
Mandelbrot                         0.0.0
Markdown                           3.0.1
MarkupSafe                         1.0
MASS                               0.4.12
MASS-Report-Reader                 0.4.14
matplotlib                         2.2.4
mccabe                             0.6.1
mistune                            0.8.3
mock                               2.0.0
modulegraph                        0.17
mpld3                              0.3
mpmath                             1.0.0
msgpack-python                     0.5.6
multiprocess                       0.70.6.1
nbconvert                          5.3.1
nbformat                           4.4.0
nbsphinx                           0.3.4
networkx                           2.1
notebook                           5.6.0
numdifftools                       0.9.20
numpy                              1.16.5
numpydoc                           0.8.0
olefile                            0.45.1
openpyxl                           2.5.4
opt-einsum                         2.3.2
osqp                               0.4.0
packaging                          17.1
pandas                             0.24.2
pandocfilters                      1.4.2
parso                              0.3.1
pathlib                            1.0.1
pathlib2                           2.3.2
patsy                              0.5.0
pbr                                5.1.1
pep8                               1.7.1
pexpect                            4.6.0
pickleshare                        0.7.4
Pillow                             5.2.0
pip                                19.3.1
plotMatrix2pdf                     0.3.2
progressbar2                       3.38.0
prometheus-client                  0.3.1
prompt-toolkit                     1.0.15
protobuf                           3.6.1
psutil                             5.4.6
ptyprocess                         0.6.0
pur                                5.2.2
py2app                             0.15
PyAstronomy                        0.13.0
pycodestyle                        2.4.0
pydicom                            1.1.0
pyflakes                           2.0.0
pygam                              0.5.5
Pygments                           2.2.0
pylint                             1.9.3
pymc                               2.3.7
pymc3                              3.6
pyparsing                          2.2.0
pyshp                              2.0.1
pysym                              0.2.3
python-dateutil                    2.7.3
python-utils                       2.3.0
pytz                               2018.5
PyWavelets                         0.5.2
PyYAML                             3.13
pyzmq                              18.0.1
QtAwesome                          0.4.4
qtconsole                          4.4.1
QtPy                               1.4.2
recommonmark                       0.4.0
redis                              3.2.1
regex                              2019.4.14
relaxography                       0.4.7
requests                           2.19.1
rope                               0.10.7
scandir                            1.9.0
scikit-image                       0.14.0
scikit-learn                       0.20.2
scipy                              1.2.2
scs                                2.0.2
seaborn                            0.9.0
Send2Trash                         1.5.0
setuptools                         41.4.0
Shapely                            1.6.4.post2
simplegeneric                      0.8.1
simpy                              3.0.11
singledispatch                     3.4.0.3
six                                1.12.0
snowballstemmer                    1.2.1
sparse                             0.3.1
Sphinx                             1.7.6
sphinxcontrib-websupport           1.1.0
spyder-kernels                     0.2.6
statsmodels                        0.9.0
subprocess32                       3.5.2
sympy                              1.2
tb-nightly                         1.15.0a20190911
tensorboard                        1.13.1
tensorflow                         1.13.1
tensorflow-estimator               1.13.0
tensorflow-estimator-2.0-preview   1.14.0.dev2019091300
tensorflow-tensorboard             1.5.1
termcolor                          1.1.0
terminado                          0.8.1
testpath                           0.3.1
Theano                             1.0.4
toolz                              0.9.0
tornado                            5.1
tqdm                               4.24.0
traitlets                          4.3.2
tslib                              1.6
typing                             3.6.4
untangle                           1.1.1
urllib3                            1.23
version                            0.1.1
vine                               1.3.0
wcwidth                            0.1.7
webencodings                       0.5.1
Werkzeug                           0.16.0
wheel                              0.33.6
widgetsnbextension                 3.4.0
wrapt                              1.11.2
WTForms                            2.2.1
xlrd                               1.2.0
xlutils                            2.0.0
xlwt                               1.3.0
xmltodict                          0.11.0
zerorpc                            0.6.1

Some operators are only R-linear, such as Re, Im, complex conjugation, and RFFT. Is there an interest in including these in PyLops? If so, I think the dot product test would need to change to include a special case for these operators, where the inner product is performed with complex values converted to R^2: <a+ib, c+id> = ac + bd.

I can create a pull request with this if desired.

For speeding up routines involving fft pyFFTW (https://github.com/pyFFTW/pyFFTW) may be a good choice.

It comes with an extra dependency and will leave the user with the manual installation of FFTW as described in pyFFTW manual, but it has wrappers to numpy/scipy fits, so we could fallback to bumpy if FFTW is not available.

After release 1.12.1 importing LinearOperator using: from pylops import LinearOperator throws the following error in Google Colab

ImportError                               Traceback (most recent call last)
<ipython-input-4-70b12ce7c165> in <module>()
      1 import numpy as np
----> 2 from pylops import LinearOperator
      3 
      4 
      5 class UndersampleObesrver(LinearOperator):

4 frames
/usr/local/lib/python3.6/dist-packages/pylops/__init__.py in <module>()
----> 1 from .LinearOperator import LinearOperator
      2 from .basicoperators import Regression
      3 from .basicoperators import LinearRegression
      4 from .basicoperators import MatrixMult
      5 from .basicoperators import Identity

/usr/local/lib/python3.6/dist-packages/pylops/LinearOperator.py in <module>()
     15 from scipy.sparse import csr_matrix
     16 
---> 17 from pylops.utils.backend import get_array_module, get_module, get_sparse_eye
     18 from pylops.optimization.solver import cgls
     19 

/usr/local/lib/python3.6/dist-packages/pylops/utils/__init__.py in <module>()
      1 from .utils import Report
----> 2 from .dottest import dottest
      3 from .deps import *
      4 from .backend import *

/usr/local/lib/python3.6/dist-packages/pylops/utils/dottest.py in <module>()
      1 import numpy as np
----> 2 from pylops.utils.backend import get_module
      3 
      4 
      5 def dottest(Op, nr, nc, tol=1e-6, complexflag=0, raiseerror=True, verb=False,

/usr/local/lib/python3.6/dist-packages/pylops/utils/backend.py in <module>()
      8     import cupy as cp
      9     import cupyx
---> 10     from cupyx.scipy.linalg import block_diag as cp_block_diag
     11     from cupyx.scipy.linalg import toeplitz as cp_toeplitz
     12     from cupyx.scipy.sparse import csc_matrix as cp_csc_matrix

ImportError: cannot import name 'block_diag'

Seems like some dependency error with cupy. I tested locally installed version via pip and it worked. Downgrading to version 1.11.1 fixes the error

Add the ability to combine linear operators with a numpy 2D arrays (or other matrices like the ones from scipy.sparse) using VStack, HStack, Block, etc. like this:

import pylops
import numpy as np

# 2x2 identity operator
A = pylops.Identity(2)

# 2x2 identity matrix
B = np.eye(2)

C = pylops.VStack([A, B])
print(C @ np.array([1, 2]))

Currently this results in
AttributeError: 'numpy.ndarray' object has no attribute 'matvec'
while a numpy 2D array does support the @ operator (as substitute for .matvec()). The same holds for other types of matrices. The implementation of VStack, HStack, Block, etc. should use the @ operator if .matvec is not implemented.

Table creation requires looping over px and py which results in slow initialisation of operator Spread3D. Moreover making tables for 3D case is very memory demanding. Two improvements are required:

• creation of tables via numba
• possibility to use numba with onthefly=True

Create operator that performs acoustic wavefield separation

The implementation of the Kronecker product in pylops.Kronecker is faulty. As a demonstration, I compute the Kronecker product of 2 simple matrices, using numpy.kron, and compare it to the same product with pylops.Kronecker:

import pylops
import numpy

A = numpy.array([[1.0, 0.0],
                 [1.0, 1.0]])

B = numpy.array([[1.0, 0.0],
                 [0.0, 1.0]])

print(numpy.kron(A,B))
print(pylops.Kronecker(pylops.MatrixMult(A),
                       pylops.MatrixMult(B)) @ numpy.eye(4))

The results is are different:

[[1. 0. 0. 0.]
 [0. 1. 0. 0.]
 [1. 0. 1. 0.]
 [0. 1. 0. 1.]]
[[1. 0. 0. 0.]
 [0. 0. 1. 0.]
 [1. 1. 0. 0.]
 [0. 0. 1. 1.]]

To get a correct implementation, the matrix y should be transposed before it is raveled (both in _matvec and _rmatvec), and the role of Op1 and Op2 should be switched in _matvec only. The resulting correct implementation is

def _matvec(self, x):
    x = x.reshape(self.Op1.shape[1], self.Op2.shape[1])
    y = self.Op1.matmat(x).T
    y = self.Op2.matmat(y).T.ravel()
    return y
def _rmatvec(self, x):
    x = x.reshape(self.Op1.shape[0], self.Op2.shape[0])
    y = self.Op1H.matmat(x).T
    y = self.Op2H.matmat(y).T.ravel()
    return y

Indeed, the implementation of _rmatvec should be identical to that of _matvec except for using the adjoint operators, Because taking the adjoint does not switch the order of the Kronecker product: (AxB)^H = A^H x B^H

Create Kroneker operator performing Kroneker product (https://en.wikipedia.org/wiki/Kronecker_product) of a list of linear operators.

Motivation

This Issue is intended to collect pointers to various transforms that we would like to include in PyLops (as well as discuss the level of maturity of the Python libraries providing such transforms).

Contourlets

• https://github.com/mazayux/pycontourlet

and create tutorial based on https://ieeexplore.ieee.org/stamp/stamp.jsp?arnumber=9833541

DCT: done 🚀

- https://docs.scipy.org/doc/scipy/reference/generated/scipy.fftpack.dct.html

and 2D: https://stackoverflow.com/questions/7110899/how-do-i-apply-a-dct-to-an-image-in-python

See this for an example: https://inst.eecs.berkeley.edu/~ee123/sp16/Sections/JPEG_DCT_Demo.html

DTCW

Make operator based on:

• https://dtcwt.readthedocs.io/en/0.12.0/

Shearlet

Make operator for Shearlet transform using either of those libraries:

• https://github.com/stefanloock/pyshearlab and http://na.math.uni-goettingen.de/pyshearlab/

• https://github.com/grlee77/PyShearlets

and create tutorial based on http://citeseerx.ist.psu.edu/viewdoc/download;jsessionid=250093F1970AD5FCFCF9921F97DA9BE1?doi=10.1.1.572.1557&rep=rep1&type=pdf

When creating a block operator, often some coefficients are zero. Right now, you need to use a Zero operator like:

pylops.Block([[A, A               ],
              [A, pylops.Zero(n,m)]])

while in for example scipy.sparse.bmat you can just write

bmat([[A, A   ],
      [A, None]])

which is much more expressive, and the correct dimensions of the zero block are inferred from the other blocks.

Hi, I'm kind of struggling to get to grips with the operator logic. I have the following code:

A = np.array([[1,0],[0,1]])
I1 = MatrixMult(A)
I2 = MatrixMult(A)
I3 = MatrixMult(A)
BOp = BlockDiag(op_list)

According to my understanding I'd assume any matrix that has 3 rows will have the identity matrix applied to every element of it. So if every element is another matrx such as follows:

X = np.array([[[2,3],[3,2]],[ [4,5],[5,5]],[[6,6],[0,1]]])

I'd expect [[2,3],[3,2]] * [[1,0],[0,1]] to be placed into the first row, first column of my result matrix however a dimension error is thrown. Where is the flaw in my logic?

Thanks

So far the convention of prepending _ to highlight protected members is not used. While no operator requires the user to change anything after initialization and its use is limited to * (or matvec) and .H (or rmatvec), it may be better to be consistent with common practice.

Also, in some cases some variables are checked at initialization, it may more natural to convert those into @property. This is not just a pure stylistic change but may make better code (although will raise the bar for new developers?)

Identity is so far implemented with a reference y=x.

This may not always be the right choice, allow both copy and reference (inplace) with an optional flag.

Convolve1D throws an error with scipy<1.1.0 because it uses scipy.signal.fftconvolve with axes parameters which was added in 1.1.0.

Try using scipy.signal.fftconvolve with axes instead of scipy.signal.convolve and scipy.signal.correlate with np.apply_along_axis and compare speed.

As there is no fftcorrelate the adjoint can be done by flipping h.

Just a personal preference. I like when the link to the main page is straight next to the title. Often I end up on the GitHub-site but I actually want the docs, so I appreciate this shortcut.

With the example of empymod:

Which, in this case, would be https://pylops.readthedocs.io.

I suggest to change the module name from lops to pylops, to be the same as the package name to install.

I installed it with pip install pylops, fired-up an IPython session, and run import pylops just to get an error that pylops is apparently not installed. So you have to go to the documentation and examples to find out how to call it. Ah, it's name is lops, so import lops.

I might not be the only one facing that issue. Most modules follow the convention to name the module the same as the package you install (e.g. numpy, scipy, numba, matplotlib, basically most I can think of, although there are exceptions).

I thought I open an issue to discuss it and change it now (v1.0.1) or never, as the package is still pretty new and not too much code is written yet with import lops. This will probably change pretty soon and changing the module name would not be a good idea anymore down the road.

So far a single damping factor (epsR and epsI) can be provided to pylops.avo.prestack.PrestackInversion. However as it is known that some parameters are more constrained than others in any AVO linearisation it is useful to provide different damping factors for the different parameters to invert for.

Does "PoststackLinearModelling" operators, be able to handle "3D wavelet variation" to cope up with change in amplitude spectrum. If yes, can you please add a Tutorial.

Or is it in the list to be added in future ?

Implement a variety of sparse solvers:

• ISTA
• FISTA
• OMP
• SPGL1 (see this discussion drrelyea/spgl1#14)

Hey,

I've created a 6x9 Fredholm1 operator from the following array

G = np.array([ [[2,1,0],[0,5,1]],[[2,3,1],[0,2,1]],[[3,1,2],[2,3,4]] ])
#G.shape => 3x2x3
Gop = Fredholm1(G)
#Gop.shape => 6x9

What should the shape of X be if I want to perform G * X ?

Sys - win10
pip - 19.1.1
python 3.7.3

Installed pylops using pip and when I opened the mdd.py file installed under \lib\site-package\pylops I noticed that the definition of MDD is different to the one in this repo. The example "08. Multi-Dimensional Deconvolution" seems to only work with the version that you get from pip and if you try to use the version of the function on this repo it complains about dimension error. The significant difference is on line 251. The repo version seems to transform the array before calling the MDC whereas the pip version just calls the MDC.

Create Wavelet operator that wraps the wavelet transform from https://pywavelets.readthedocs.io/en/latest/.

As for pyFFTW this will add a dependency but give a new useful functionality. PyWavelets seems to be very mature and stable though, so that should not be a problem.

After reading the paper quoted, I did not find any explanation in regards with F^(H)GF (does the paper itself explain it or one of referenced papers?). I'm assuming the idea is to create an object that holds a copy of one of the data sets and whenever you call MDC.* you perform convolution between the copy held by the object and the argument to the function. (I might be misunderstanding?). Could you provide a quick alternative in pylops to the following matlab code?

A(:,:,1) = [1,2,1;2,3,0];
A(:,:,2) = [3,2,2;1,2,3];

B(:,:,1) = [1,1;1,1];
B(:,:,2) = [1,1;1,1];



convn(A,B)

When trying to write the alternative I wasn't certain what values nt and nv should be. Also, I know the result will be another operator so is there a way to actually display the result of the convolution?

• Add a boolean member called symmetric that defines is the operator is symmetric or not. This will be always False if the operator is not square, while it can be fixed for a square operator or depend to its input parameters like in the case of Block and BlockDiagonal operators
• Add also a 'validation' method to LinearOperator called isymmetric which validates the following equality <x, Ay>==<y, Ax>
• Check symmetric when solving an inverse problem in one of the classes that invoke script solvers and use it to choose the most appropriate solver (e.g., cg for symmetric operators)

The linear operator part of pylops could be interesting and appeal to a wider community. The geophysical part, however, is more for a specific crowd of geophysicist.

I suggest to remove the geophysical parts and create with it a new package. This might attract more interest in both, pylops from the wider scientific community which would otherwise be put-off by the geophysics-focus, and the new package from the geophysics community which is not necessarily interested in low-level stuff which is present in pylops.

Please mention and add a link to the corresponding Arxiv preprint https://arxiv.org/abs/1907.12349 into the readme of this repo.

Using the SciPy lstsq method, I can compute a solution to the problem Ax = b:

from scipy.linalg import lstsq
import numpy as np

"""Solving the system Ax = b:"""

A = np.array([[2.0, 0.0], [-1.0, 1.0], [0.0, 2.0]])
b = np.array([1.0, 0.0, -1.0])

x = lstsq(a=A, b=b)

"""Gives the correct solution, x = (1/3, -1.3)"""
print(x)

However, I have been unable to reproduce this calculation in PyLops. All of the tutorials I've found deal with scenarions where the matrix A is 1D, or it has the same dimensions as the data. Copious trial and error has so far filed to get me out of a dimension mismatch. Here is the sort of solution I'm trying:

import pylops

"""Now want to repeat this calculation using PyLops"""
LRop = pylops.Regression(A, order=1, dtype='float64')
print(LRop.shape)
solver = pylops.NormalEquationsInversion(Op=LRop, Regs=None, data=b)
print(solver)

I have no doubt I'm being very stupid, but if you could give me an example to get my simple problem up and running (or explains what I'm conceptually not grasping). I reckon I'll be able to make progress with PyLops from here.

Hi,

I have a LinearOperator object that is complex and a data vector that is also complex. I want to feed these into the NormalEquationsInversion but I want my output solution to be real.

I was wondering if there is a way to split the "Op" LinearOperator to one that has the real and one on the imag components so that when I want to perform an Inversion I can get the desired outcome by means that I will explain with an example. Assume I have the following

Aop = pylops.MatrixMult(
    sparse.bsr_matrix(matrix), dtype="complex128"
)

FFTop = pylops.signalprocessing.FFT2D(dims)

Op = Aop * FFTop

I also have complex data and complex weights that I can feed to a modified version of the "NormalEquationsInversion",

x = NormalEquationsInversion(Op, data_real, data_imag, Weight_real, Weight_imag)

where

NormalEquationsInversion(Op, data_real, data_imag, Weight_real, Weight_imag, Regs):

    #NOTE: Is there a way to get the real and imag components of this operator as separate operators?
    Op_real = ...
    Op_imag = ...

    OpH_real = Op_real.H
    OpH_imag = Op_imag.H

    Op_normal = (OpH_real * Weight_real * Op_real  + OpH_imag * Weight_imag * Op_imag) + RegH * Reg

    y_normal = OpH_real * Weight_real * data_real + OpH_imag * Weight_imag * data_imag

    xinv, istop = cg(Op_normal, y_normal, **kwargs_cg)
 
    return xinv

Is there a way to fill in the blanks above to make it do what I want?

I am looking into the possibility of using pylops for simultaneous inversion using PP and PS inversion. Although the model parameters are the same between PP and PP+PS joint inversion, the G matrix and the data size would be different when adding the PS data into the inversion. I have not figure out a way to update the G matrix correctly based on the joint inversion scheme. Maybe this issue would help me to sort things out.

Compare performance in FTT operator of performing fft along an axis different than -1 using the following two strategies:

• np.swap+np.fft.fft(..., axis=-1) versus
• np.fft.fft(..., axis=chosen)

Create generic operators in signalprocessing submodule performing Fredholm integrals of first and second kind as in https://en.wikipedia.org/wiki/Fredholm_integral_equation for 3 dimensional operators (note that when the kernel is 2d this is just MatrixMult 1st kind and Identity - MatrixMult 2nd kind).

Such operators have practically being made within MDC and Marchenko but having something more generic will allow to build MDC as a chain of FFT.H Fredholm1 FFT and Marchenko as a Fredholm2 with user defined kernel.

Main advantage: when focusing on speed up this will be limited to the actual integral computation, which is nothing more than a dense matrix-vector or matrix-matrix operation. We can use using anything available open-source for Machine Learning (as fast matrix-matrix is all a neural network needs to be fast). Options:

• pycuda
• scikit-cuda
• tensorflow
• pytorch
• ...

Consider adding
https://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.linalg.lobpcg.html within LinearOperator.eigs method.

Some calls take an N and a dims parameter, e.g., pylops.FirstDerivative. The N is redundant, as it can be derived from dims.

This is there for legacy reasons and should be cleaned-up for v2.0.0 (@mrava87 you can tag that for the milestone/project board).

Also, as it happens in the same instance call, replacing the dir with axis (like used in NumPy/SciPy).

Some of the main linear operators may get a speed-up by simply adding numba <http://numba.pydata.org>_ @jit decorators to _matvec and _rmatvec methods.

Others, such as FirstDerivative may be suited to stencil decorator.

Note: it may be enough to simply add numba to dependencies as it is a fairly common library but to keep dependencies to a minimum it may be better to implement private functions such as _matvec_serial and _matvec_numba and try importingnumba, always falling back to by _matvec_serial if that is not available in the python environment of the user.