Skip to content

Commit

Permalink
feat: added dottest
Browse files Browse the repository at this point in the history
  • Loading branch information
mrava87 committed Oct 30, 2024
1 parent 57b793e commit df966d1
Show file tree
Hide file tree
Showing 10 changed files with 198 additions and 6 deletions.
19 changes: 19 additions & 0 deletions docs/source/api/index.rst
Original file line number Diff line number Diff line change
Expand Up @@ -104,3 +104,22 @@ Basic

cg
cgls


Utils
-----

.. currentmodule:: pylops_mpi.DistributedArray

.. autosummary::
:toctree: generated/

local_split


.. currentmodule:: pylops_mpi.utils.dottest

.. autosummary::
:toctree: generated/

dottest
2 changes: 2 additions & 0 deletions pylops_mpi/utils/__init__.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,2 @@
# isort: skip_file
from .dottest import *
107 changes: 107 additions & 0 deletions pylops_mpi/utils/dottest.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,107 @@
__all__ = ["dottest"]

from typing import Optional

import numpy as np

from pylops_mpi.DistributedArray import DistributedArray
from pylops.utils.backend import to_numpy


def dottest(
Op,
u: DistributedArray,
v: DistributedArray,
nr: Optional[int] = None,
nc: Optional[int] = None,
rtol: float = 1e-6,
atol: float = 1e-21,
raiseerror: bool = True,
verb: bool = False,
) -> bool:
r"""Dot test.
Perform dot-test to verify the validity of forward and adjoint
operators using user-provided random vectors :math:`\mathbf{u}`
and :math:`\mathbf{v}` (whose Partition must be consistent with
the operator being tested). This test can help to detect errors
in the operator ximplementation.
Parameters
----------
Op : :obj:`pylops_mpi.LinearOperator`
Linear operator to test.
u : :obj:`pylops_mpi.DistributedArray`
Distributed array of size equal to the number of columns of operator
v : :obj:`pylops_mpi.DistributedArray`
Distributed array of size equal to the number of rows of operator
nr : :obj:`int`
Number of rows of operator (i.e., elements in data)
nc : :obj:`int`
Number of columns of operator (i.e., elements in model)
rtol : :obj:`float`, optional
Relative dottest tolerance
atol : :obj:`float`, optional
Absolute dottest tolerance
.. versionadded:: 2.0.0
raiseerror : :obj:`bool`, optional
Raise error or simply return ``False`` when dottest fails
verb : :obj:`bool`, optional
Verbosity
Returns
-------
passed : :obj:`bool`
Passed flag.
Raises
------
AssertionError
If dot-test is not verified within chosen tolerances.
Notes
-----
A dot-test is mathematical tool used in the development of numerical
linear operators.
More specifically, a correct implementation of forward and adjoint for
a linear operator should verify the following *equality*
within a numerical tolerance:
.. math::
(\mathbf{Op}\,\mathbf{u})^H\mathbf{v} =
\mathbf{u}^H(\mathbf{Op}^H\mathbf{v})
"""
if nr is None:
nr = Op.shape[0]
if nc is None:
nc = Op.shape[1]

if (nr, nc) != Op.shape:
raise AssertionError("Provided nr and nc do not match operator shape")

y = Op.matvec(u) # Op * u
x = Op.rmatvec(v) # Op'* v

yy = np.vdot(y.asarray(), v.asarray()) # (Op * u)' * v
xx = np.vdot(u.asarray(), x.asarray()) # u' * (Op' * v)

# convert back to numpy (in case cupy arrays were used), make into a numpy
# array and extract the first element. This is ugly but allows to handle
# complex numbers in subsequent prints also when using cupy arrays.
xx, yy = np.array([to_numpy(xx)])[0], np.array([to_numpy(yy)])[0]

# evaluate if dot test passed
passed = np.isclose(xx, yy, rtol, atol)

# verbosity or error raising
if (not passed and raiseerror) or verb:
passed_status = "passed" if passed else "failed"
msg = f"Dot test {passed_status}, v^H(Opu)={yy} - u^H(Op^Hv)={xx}"
if not passed and raiseerror:
raise AssertionError(msg)
else:
print(msg)

return passed
17 changes: 15 additions & 2 deletions tests/test_blockdiag.py
Original file line number Diff line number Diff line change
@@ -1,18 +1,23 @@
"""Test the MPIBlockDiag and MPIStackedBlockDiag classes
Designed to run with n processes
$ mpiexec -n 10 pytest test_blockdiag.py --with-mpi
"""
from mpi4py import MPI
import numpy as np
from numpy.testing import assert_allclose
np.random.seed(42)

import pytest

import pylops
import pylops_mpi
from pylops_mpi.utils.dottest import dottest

par1 = {'ny': 101, 'nx': 101, 'dtype': np.float64}
par1j = {'ny': 101, 'nx': 101, 'dtype': np.complex128}
par2 = {'ny': 301, 'nx': 101, 'dtype': np.float64}
par2j = {'ny': 301, 'nx': 101, 'dtype': np.complex128}

np.random.seed(42)


@pytest.mark.mpi(min_size=2)
@pytest.mark.parametrize("par", [(par1), (par1j), (par2), (par2j)])
Expand All @@ -30,10 +35,14 @@ def test_blockdiag(par):
y[:] = np.ones(shape=par['ny'], dtype=par['dtype'])
y_global = y.asarray()

# Forward
x_mat = BDiag_MPI @ x
# Adjoint
y_rmat = BDiag_MPI.H @ y
assert isinstance(x_mat, pylops_mpi.DistributedArray)
assert isinstance(y_rmat, pylops_mpi.DistributedArray)
# Dot test
dottest(BDiag_MPI, x, y, size * par['ny'], size * par['nx'])

x_mat_mpi = x_mat.asarray()
y_rmat_mpi = y_rmat.asarray()
Expand Down Expand Up @@ -73,10 +82,14 @@ def test_stacked_blockdiag(par):
y = pylops_mpi.StackedDistributedArray(distarrays=[dist1, dist2])
y_global = y.asarray()

# Forward
x_mat = StackedBDiag_MPI @ x
# Adjoint
y_rmat = StackedBDiag_MPI.H @ y
assert isinstance(x_mat, pylops_mpi.StackedDistributedArray)
assert isinstance(y_rmat, pylops_mpi.StackedDistributedArray)
# Dot test
dottest(StackedBDiag_MPI, x, y, size * par['ny'] + par['nx'] * par['ny'], size * par['nx'] + par['nx'] * par['ny'])

x_mat_mpi = x_mat.asarray()
y_rmat_mpi = y_rmat.asarray()
Expand Down
27 changes: 24 additions & 3 deletions tests/test_derivative.py
Original file line number Diff line number Diff line change
@@ -1,10 +1,15 @@
"""Test the derivative classes
Designed to run with n processes
$ mpiexec -n 10 pytest test_derivative.py --with-mpi
"""
import numpy as np
from mpi4py import MPI
from numpy.testing import assert_allclose
import pytest

import pylops
import pylops_mpi
from pylops_mpi.utils.dottest import dottest

np.random.seed(42)
rank = MPI.COMM_WORLD.Get_rank()
Expand Down Expand Up @@ -194,6 +199,8 @@ def test_first_derivative_forward(par):
# Adjoint
y_adj_dist = Fop_MPI.H @ x
y_adj = y_adj_dist.asarray()
# Dot test
dottest(Fop_MPI, x, y_dist, np.prod(par['nz']), np.prod(par['nz']))

if rank == 0:
Fop = pylops.FirstDerivative(dims=par['nz'], axis=0,
Expand Down Expand Up @@ -226,6 +233,8 @@ def test_first_derivative_backward(par):
# Adjoint
y_adj_dist = Fop_MPI.H @ x
y_adj = y_adj_dist.asarray()
# Dot test
dottest(Fop_MPI, x, y_dist, np.prod(par['nz']), np.prod(par['nz']))

if rank == 0:
Fop = pylops.FirstDerivative(dims=par['nz'], axis=0,
Expand Down Expand Up @@ -259,6 +268,9 @@ def test_first_derivative_centered(par):
# Adjoint
y_adj_dist = Fop_MPI.H @ x
y_adj = y_adj_dist.asarray()
# Dot test
dottest(Fop_MPI, x, y_dist, np.prod(par['nz']), np.prod(par['nz']))

if rank == 0:
Fop = pylops.FirstDerivative(dims=par['nz'], axis=0,
sampling=par['dz'],
Expand Down Expand Up @@ -290,6 +302,8 @@ def test_second_derivative_forward(par):
# Adjoint
y_adj_dist = Sop_MPI.H @ x
y_adj = y_adj_dist.asarray()
# Dot test
dottest(Sop_MPI, x, y_dist, np.prod(par['nz']), np.prod(par['nz']))

if rank == 0:
Sop = pylops.SecondDerivative(dims=par['nz'], axis=0,
Expand Down Expand Up @@ -322,6 +336,8 @@ def test_second_derivative_backward(par):
# Adjoint
y_adj_dist = Sop_MPI.H @ x
y_adj = y_adj_dist.asarray()
# Dot test
dottest(Sop_MPI, x, y_dist, np.prod(par['nz']), np.prod(par['nz']))

if rank == 0:
Sop = pylops.SecondDerivative(dims=par['nz'], axis=0,
Expand Down Expand Up @@ -354,6 +370,8 @@ def test_second_derivative_centered(par):
# Adjoint
y_adj_dist = Sop_MPI.H @ x
y_adj = y_adj_dist.asarray()
# Dot test
dottest(Sop_MPI, x, y_dist, np.prod(par['nz']), np.prod(par['nz']))

if rank == 0:
Sop = pylops.SecondDerivative(dims=par['nz'], axis=0,
Expand Down Expand Up @@ -385,6 +403,8 @@ def test_laplacian(par):
# Adjoint
y_adj_dist = Lop_MPI.H @ x
y_adj = y_adj_dist.asarray()
# Dot test
dottest(Lop_MPI, x, y_dist, np.prod(par['n']), np.prod(par['n']))

if rank == 0:
Lop = pylops.Laplacian(dims=par['n'], axes=par['axes'],
Expand All @@ -409,6 +429,7 @@ def test_gradient(par):
x_fwd = pylops_mpi.DistributedArray(global_shape=np.prod(par['n']), dtype=par['dtype'])
x_fwd[:] = np.random.normal(rank, 10, x_fwd.local_shape)
x_global = x_fwd.asarray()

# Forward
y_dist = Gop_MPI @ x_fwd
assert isinstance(y_dist, pylops_mpi.StackedDistributedArray)
Expand All @@ -421,15 +442,15 @@ def test_gradient(par):
x_adj_dist2[:] = np.random.normal(rank, 20, x_adj_dist2.local_shape)
x_adj_dist3 = pylops_mpi.DistributedArray(global_shape=int(np.prod(par['n'])), dtype=par['dtype'])
x_adj_dist3[:] = np.random.normal(rank, 30, x_adj_dist3.local_shape)

x_adj = pylops_mpi.StackedDistributedArray(distarrays=[x_adj_dist1, x_adj_dist2, x_adj_dist3])

x_adj_global = x_adj.asarray()
y_adj_dist = Gop_MPI.H @ x_adj
assert isinstance(y_adj_dist, pylops_mpi.DistributedArray)

y_adj = y_adj_dist.asarray()

# Dot test
dottest(Gop_MPI, x_fwd, y_dist, len(par['n']) * np.prod(par['n']), np.prod(par['n']))

if rank == 0:
Gop = pylops.Gradient(dims=par['n'], sampling=par['sampling'],
kind=kind, edge=par['edge'],
Expand Down
8 changes: 7 additions & 1 deletion tests/test_fredholm.py
Original file line number Diff line number Diff line change
@@ -1,3 +1,7 @@
"""Test the MPIFredholm1 class
Designed to run with n processes
$ mpiexec -n 10 pytest test_fredholm.py --with-mpi
"""
import numpy as np
from numpy.testing import assert_allclose
from mpi4py import MPI
Expand All @@ -9,12 +13,12 @@
from pylops_mpi import DistributedArray
from pylops_mpi.DistributedArray import local_split, Partition
from pylops_mpi.signalprocessing import MPIFredholm1
from pylops_mpi.utils.dottest import dottest

np.random.seed(42)
rank = MPI.COMM_WORLD.Get_rank()
size = MPI.COMM_WORLD.Get_size()


par1 = {
"nsl": 6,
"ny": 6,
Expand Down Expand Up @@ -130,6 +134,8 @@ def test_Fredholm1(par):
# Adjoint
y_adj_dist = Fop_MPI.H @ y_dist
y_adj = y_adj_dist.asarray()
# Dot test
dottest(Fop_MPI, x, y_dist, par["nsl"] * par["nx"] * par["nz"],par["nsl"] * par["ny"] * par["nz"])

if rank == 0:
Fop = pylops.signalprocessing.Fredholm1(
Expand Down
7 changes: 7 additions & 0 deletions tests/test_linearop.py
Original file line number Diff line number Diff line change
@@ -1,3 +1,7 @@
"""Test the MPILinearOperator class
Designed to run with n processes
$ mpiexec -n 10 pytest test_linearop.py --with-mpi
"""
import numpy as np
from numpy.testing import assert_allclose
from mpi4py import MPI
Expand Down Expand Up @@ -58,6 +62,7 @@ def test_transpose(par):
"""Test the TransposeLinearOperator"""
Op = pylops.MatrixMult(A=((rank + 1) * np.ones(shape=(par['ny'], par['nx']))).astype(par['dtype']))
BDiag_MPI = MPIBlockDiag(ops=[Op, ])

# Tranposed Op
Top_MPI = BDiag_MPI.T

Expand All @@ -76,6 +81,7 @@ def test_transpose(par):
Top_y = Top_MPI.H @ y
assert isinstance(Top_y, DistributedArray)
Top_y_np = Top_y.asarray()

if rank == 0:
ops = [pylops.MatrixMult((i + 1) * np.ones(shape=(par['ny'], par['nx'])).astype(par['dtype'])) for i in
range(size)]
Expand Down Expand Up @@ -109,6 +115,7 @@ def test_scaled(par):
Sop_y = Sop_MPI.H @ y
assert isinstance(Sop_y, DistributedArray)
Sop_y_np = Sop_y.asarray()

if rank == 0:
ops = [pylops.MatrixMult((i + 1) * np.ones(shape=(par['ny'], par['nx'])).astype(par['dtype'])) for i in
range(size)]
Expand Down
4 changes: 4 additions & 0 deletions tests/test_solver.py
Original file line number Diff line number Diff line change
@@ -1,3 +1,7 @@
"""Test solvers
Designed to run with n processes
$ mpiexec -n 10 pytest test_solver.py --with-mpi
"""
import numpy as np
from numpy.testing import assert_allclose
from mpi4py import MPI
Expand Down
Loading

0 comments on commit df966d1

Please sign in to comment.