deploy: d5aeffb

.buildinfo

@@ -1,4 +1,4 @@
 # Sphinx build info version 1
 # This file hashes the configuration used when building these files. When it is not found, a full rebuild will be done.
-config: b2736e47c27770184fc01163b30d7828
+config: 2b995ed0f2cc84754bcc4807ea617af1
 tags: 645f666f9bcd5a90fca523b33c5a78b7

_downloads/349e24ddf4c2a9f9f1903b6c1d4ed115/plot_post_stack_inversion.py → _downloads/1001629eacd23bb9b9df42b1d6da8075/poststack.py

@@ -187,7 +187,19 @@
                                 sampling=(1, 1, 1), dtype=BDiag.dtype)
 NormEqOp = BDiag.H @ BDiag + epsR * LapOp.H @ LapOp
 dnorm_dist = BDiag.H @ d_dist
-minv3d_reg_dist = pylops_mpi.optimization.basic.cg(NormEqOp, dnorm_dist, x0=mback3d_dist, niter=100, show=True)[0]
+minv3d_ne_dist = pylops_mpi.optimization.basic.cg(NormEqOp, dnorm_dist, x0=mback3d_dist, niter=100, show=True)[0]
+minv3d_ne = minv3d_ne_dist.asarray().reshape((ny, nx, nz))
+
+###############################################################################
+
+# Regularized inversion with regularized equations
+StackOp = pylops_mpi.MPIStackedVStack([BDiag, np.sqrt(epsR) * LapOp])
+d0_dist = pylops_mpi.DistributedArray(global_shape=ny * nx * nz)
+d0_dist[:] = 0.
+dstack_dist = pylops_mpi.StackedDistributedArray([d_dist, d0_dist])
+
+dnorm_dist = BDiag.H @ d_dist
+minv3d_reg_dist = pylops_mpi.optimization.basic.cgls(StackOp, dstack_dist, x0=mback3d_dist, niter=100, show=False)[0]
 minv3d_reg = minv3d_reg_dist.asarray().reshape((ny, nx, nz))
 
 ###############################################################################
@@ -205,7 +217,7 @@
     print('Smooth Distr == Local', np.allclose(d_0, d0_0))
 
     # Visualize
-    fig, axs = plt.subplots(nrows=5, ncols=3, figsize=(9, 14), constrained_layout=True)
+    fig, axs = plt.subplots(nrows=6, ncols=3, figsize=(9, 14), constrained_layout=True)
     axs[0][0].imshow(m3d[5, :, :].T, cmap="gist_rainbow", vmin=m.min(), vmax=m.max())
     axs[0][0].set_title("Model x-z")
     axs[0][0].axis("tight")
@@ -246,12 +258,22 @@
     axs[3][2].set_title('Inverted Model iter x-y')
     axs[3][2].axis('tight')
 
-    axs[4][0].imshow(minv3d_reg[5, :, :].T, cmap="gist_rainbow", vmin=m.min(), vmax=m.max())
-    axs[4][0].set_title("Regularized Inverted Model iter x-z")
+    axs[4][0].imshow(minv3d_ne[5, :, :].T, cmap="gist_rainbow", vmin=m.min(), vmax=m.max())
+    axs[4][0].set_title("Normal Equations Inverted Model iter x-z")
     axs[4][0].axis("tight")
-    axs[4][1].imshow(minv3d_reg[:, 200, :].T, cmap='gist_rainbow', vmin=m.min(), vmax=m.max())
-    axs[4][1].set_title('Regularized Inverted Model iter y-z')
+    axs[4][1].imshow(minv3d_ne[:, 200, :].T, cmap='gist_rainbow', vmin=m.min(), vmax=m.max())
+    axs[4][1].set_title('Normal Equations Inverted Model iter y-z')
     axs[4][1].axis('tight')
-    axs[4][2].imshow(minv3d_reg[:, :, 220].T, cmap='gist_rainbow', vmin=m.min(), vmax=m.max())
-    axs[4][2].set_title('Regularized Inverted Model iter x-y')
+    axs[4][2].imshow(minv3d_ne[:, :, 220].T, cmap='gist_rainbow', vmin=m.min(), vmax=m.max())
+    axs[4][2].set_title('Normal Equations Inverted Model iter x-y')
     axs[4][2].axis('tight')
+
+    axs[5][0].imshow(minv3d_reg[5, :, :].T, cmap="gist_rainbow", vmin=m.min(), vmax=m.max())
+    axs[5][0].set_title("Regularized Inverted Model iter x-z")
+    axs[5][0].axis("tight")
+    axs[5][1].imshow(minv3d_reg[:, 200, :].T, cmap='gist_rainbow', vmin=m.min(), vmax=m.max())
+    axs[5][1].set_title('Regularized Inverted Model iter y-z')
+    axs[5][1].axis('tight')
+    axs[5][2].imshow(minv3d_reg[:, :, 220].T, cmap='gist_rainbow', vmin=m.min(), vmax=m.max())
+    axs[5][2].set_title('Regularized Inverted Model iter x-y')
+    axs[5][2].axis('tight')

_downloads/b1d8138fba244a3a8c36115409e884d3/plot_post_stack_inversion.ipynb → _downloads/25d0d43c288bbcf46dfd5c529fe808c5/poststack.ipynb

@@ -80,7 +80,18 @@
       },
       "outputs": [],
       "source": [
-        "# Regularized inversion with normal equations\nepsR = 1e2\nLapOp = pylops_mpi.MPILaplacian(dims=(ny, nx, nz), axes=(0, 1, 2), weights=(1, 1, 1),\n                                sampling=(1, 1, 1), dtype=BDiag.dtype)\nNormEqOp = BDiag.H @ BDiag + epsR * LapOp.H @ LapOp\ndnorm_dist = BDiag.H @ d_dist\nminv3d_reg_dist = pylops_mpi.optimization.basic.cg(NormEqOp, dnorm_dist, x0=mback3d_dist, niter=100, show=True)[0]\nminv3d_reg = minv3d_reg_dist.asarray().reshape((ny, nx, nz))"
+        "# Regularized inversion with normal equations\nepsR = 1e2\nLapOp = pylops_mpi.MPILaplacian(dims=(ny, nx, nz), axes=(0, 1, 2), weights=(1, 1, 1),\n                                sampling=(1, 1, 1), dtype=BDiag.dtype)\nNormEqOp = BDiag.H @ BDiag + epsR * LapOp.H @ LapOp\ndnorm_dist = BDiag.H @ d_dist\nminv3d_ne_dist = pylops_mpi.optimization.basic.cg(NormEqOp, dnorm_dist, x0=mback3d_dist, niter=100, show=True)[0]\nminv3d_ne = minv3d_ne_dist.asarray().reshape((ny, nx, nz))"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "execution_count": null,
+      "metadata": {
+        "collapsed": false
+      },
+      "outputs": [],
+      "source": [
+        "# Regularized inversion with regularized equations\nStackOp = pylops_mpi.MPIStackedVStack([BDiag, np.sqrt(epsR) * LapOp])\nd0_dist = pylops_mpi.DistributedArray(global_shape=ny * nx * nz)\nd0_dist[:] = 0.\ndstack_dist = pylops_mpi.StackedDistributedArray([d_dist, d0_dist])\n\ndnorm_dist = BDiag.H @ d_dist\nminv3d_reg_dist = pylops_mpi.optimization.basic.cgls(StackOp, dstack_dist, x0=mback3d_dist, niter=100, show=False)[0]\nminv3d_reg = minv3d_reg_dist.asarray().reshape((ny, nx, nz))"
       ]
     },
     {
@@ -98,7 +109,7 @@
       },
       "outputs": [],
       "source": [
-        "if rank == 0:\n    # Check the distributed implementation gives the same result\n    # as the one running only on rank0\n    PPop0 = PoststackLinearModelling(wav, nt0=nz, spatdims=(ny, nx))\n    d0 = (PPop0 @ m3d.transpose(2, 0, 1)).transpose(1, 2, 0)\n    d0_0 = (PPop0 @ m3d.transpose(2, 0, 1)).transpose(1, 2, 0)\n\n    # Check the two distributed implementations give the same modelling results\n    print('Distr == Local', np.allclose(d, d0))\n    print('Smooth Distr == Local', np.allclose(d_0, d0_0))\n\n    # Visualize\n    fig, axs = plt.subplots(nrows=5, ncols=3, figsize=(9, 14), constrained_layout=True)\n    axs[0][0].imshow(m3d[5, :, :].T, cmap=\"gist_rainbow\", vmin=m.min(), vmax=m.max())\n    axs[0][0].set_title(\"Model x-z\")\n    axs[0][0].axis(\"tight\")\n    axs[0][1].imshow(m3d[:, 200, :].T, cmap=\"gist_rainbow\", vmin=m.min(), vmax=m.max())\n    axs[0][1].set_title(\"Model y-z\")\n    axs[0][1].axis(\"tight\")\n    axs[0][2].imshow(m3d[:, :, 220].T, cmap=\"gist_rainbow\", vmin=m.min(), vmax=m.max())\n    axs[0][2].set_title(\"Model y-z\")\n    axs[0][2].axis(\"tight\")\n\n    axs[1][0].imshow(mback3d[5, :, :].T, cmap=\"gist_rainbow\", vmin=m.min(), vmax=m.max())\n    axs[1][0].set_title(\"Smooth Model x-z\")\n    axs[1][0].axis(\"tight\")\n    axs[1][1].imshow(mback3d[:, 200, :].T, cmap=\"gist_rainbow\", vmin=m.min(), vmax=m.max())\n    axs[1][1].set_title(\"Smooth Model y-z\")\n    axs[1][1].axis(\"tight\")\n    axs[1][2].imshow(mback3d[:, :, 220].T, cmap=\"gist_rainbow\", vmin=m.min(), vmax=m.max())\n    axs[1][2].set_title(\"Smooth Model y-z\")\n    axs[1][2].axis(\"tight\")\n\n    axs[2][0].imshow(d[5, :, :].T, cmap=\"gray\", vmin=-1, vmax=1)\n    axs[2][0].set_title(\"Data x-z\")\n    axs[2][0].axis(\"tight\")\n    axs[2][1].imshow(d[:, 200, :].T, cmap='gray', vmin=-1, vmax=1)\n    axs[2][1].set_title('Data y-z')\n    axs[2][1].axis('tight')\n    axs[2][2].imshow(d[:, :, 220].T, cmap='gray', vmin=-1, vmax=1)\n    axs[2][2].set_title('Data x-y')\n    axs[2][2].axis('tight')\n\n    axs[3][0].imshow(minv3d_iter[5, :, :].T, cmap=\"gist_rainbow\", vmin=m.min(), vmax=m.max())\n    axs[3][0].set_title(\"Inverted Model iter x-z\")\n    axs[3][0].axis(\"tight\")\n    axs[3][1].imshow(minv3d_iter[:, 200, :].T, cmap='gist_rainbow', vmin=m.min(), vmax=m.max())\n    axs[3][1].set_title('Inverted Model iter y-z')\n    axs[3][1].axis('tight')\n    axs[3][2].imshow(minv3d_iter[:, :, 220].T, cmap='gist_rainbow', vmin=m.min(), vmax=m.max())\n    axs[3][2].set_title('Inverted Model iter x-y')\n    axs[3][2].axis('tight')\n\n    axs[4][0].imshow(minv3d_reg[5, :, :].T, cmap=\"gist_rainbow\", vmin=m.min(), vmax=m.max())\n    axs[4][0].set_title(\"Regularized Inverted Model iter x-z\")\n    axs[4][0].axis(\"tight\")\n    axs[4][1].imshow(minv3d_reg[:, 200, :].T, cmap='gist_rainbow', vmin=m.min(), vmax=m.max())\n    axs[4][1].set_title('Regularized Inverted Model iter y-z')\n    axs[4][1].axis('tight')\n    axs[4][2].imshow(minv3d_reg[:, :, 220].T, cmap='gist_rainbow', vmin=m.min(), vmax=m.max())\n    axs[4][2].set_title('Regularized Inverted Model iter x-y')\n    axs[4][2].axis('tight')"
+        "if rank == 0:\n    # Check the distributed implementation gives the same result\n    # as the one running only on rank0\n    PPop0 = PoststackLinearModelling(wav, nt0=nz, spatdims=(ny, nx))\n    d0 = (PPop0 @ m3d.transpose(2, 0, 1)).transpose(1, 2, 0)\n    d0_0 = (PPop0 @ m3d.transpose(2, 0, 1)).transpose(1, 2, 0)\n\n    # Check the two distributed implementations give the same modelling results\n    print('Distr == Local', np.allclose(d, d0))\n    print('Smooth Distr == Local', np.allclose(d_0, d0_0))\n\n    # Visualize\n    fig, axs = plt.subplots(nrows=6, ncols=3, figsize=(9, 14), constrained_layout=True)\n    axs[0][0].imshow(m3d[5, :, :].T, cmap=\"gist_rainbow\", vmin=m.min(), vmax=m.max())\n    axs[0][0].set_title(\"Model x-z\")\n    axs[0][0].axis(\"tight\")\n    axs[0][1].imshow(m3d[:, 200, :].T, cmap=\"gist_rainbow\", vmin=m.min(), vmax=m.max())\n    axs[0][1].set_title(\"Model y-z\")\n    axs[0][1].axis(\"tight\")\n    axs[0][2].imshow(m3d[:, :, 220].T, cmap=\"gist_rainbow\", vmin=m.min(), vmax=m.max())\n    axs[0][2].set_title(\"Model y-z\")\n    axs[0][2].axis(\"tight\")\n\n    axs[1][0].imshow(mback3d[5, :, :].T, cmap=\"gist_rainbow\", vmin=m.min(), vmax=m.max())\n    axs[1][0].set_title(\"Smooth Model x-z\")\n    axs[1][0].axis(\"tight\")\n    axs[1][1].imshow(mback3d[:, 200, :].T, cmap=\"gist_rainbow\", vmin=m.min(), vmax=m.max())\n    axs[1][1].set_title(\"Smooth Model y-z\")\n    axs[1][1].axis(\"tight\")\n    axs[1][2].imshow(mback3d[:, :, 220].T, cmap=\"gist_rainbow\", vmin=m.min(), vmax=m.max())\n    axs[1][2].set_title(\"Smooth Model y-z\")\n    axs[1][2].axis(\"tight\")\n\n    axs[2][0].imshow(d[5, :, :].T, cmap=\"gray\", vmin=-1, vmax=1)\n    axs[2][0].set_title(\"Data x-z\")\n    axs[2][0].axis(\"tight\")\n    axs[2][1].imshow(d[:, 200, :].T, cmap='gray', vmin=-1, vmax=1)\n    axs[2][1].set_title('Data y-z')\n    axs[2][1].axis('tight')\n    axs[2][2].imshow(d[:, :, 220].T, cmap='gray', vmin=-1, vmax=1)\n    axs[2][2].set_title('Data x-y')\n    axs[2][2].axis('tight')\n\n    axs[3][0].imshow(minv3d_iter[5, :, :].T, cmap=\"gist_rainbow\", vmin=m.min(), vmax=m.max())\n    axs[3][0].set_title(\"Inverted Model iter x-z\")\n    axs[3][0].axis(\"tight\")\n    axs[3][1].imshow(minv3d_iter[:, 200, :].T, cmap='gist_rainbow', vmin=m.min(), vmax=m.max())\n    axs[3][1].set_title('Inverted Model iter y-z')\n    axs[3][1].axis('tight')\n    axs[3][2].imshow(minv3d_iter[:, :, 220].T, cmap='gist_rainbow', vmin=m.min(), vmax=m.max())\n    axs[3][2].set_title('Inverted Model iter x-y')\n    axs[3][2].axis('tight')\n\n    axs[4][0].imshow(minv3d_ne[5, :, :].T, cmap=\"gist_rainbow\", vmin=m.min(), vmax=m.max())\n    axs[4][0].set_title(\"Normal Equations Inverted Model iter x-z\")\n    axs[4][0].axis(\"tight\")\n    axs[4][1].imshow(minv3d_ne[:, 200, :].T, cmap='gist_rainbow', vmin=m.min(), vmax=m.max())\n    axs[4][1].set_title('Normal Equations Inverted Model iter y-z')\n    axs[4][1].axis('tight')\n    axs[4][2].imshow(minv3d_ne[:, :, 220].T, cmap='gist_rainbow', vmin=m.min(), vmax=m.max())\n    axs[4][2].set_title('Normal Equations Inverted Model iter x-y')\n    axs[4][2].axis('tight')\n\n    axs[5][0].imshow(minv3d_reg[5, :, :].T, cmap=\"gist_rainbow\", vmin=m.min(), vmax=m.max())\n    axs[5][0].set_title(\"Regularized Inverted Model iter x-z\")\n    axs[5][0].axis(\"tight\")\n    axs[5][1].imshow(minv3d_reg[:, 200, :].T, cmap='gist_rainbow', vmin=m.min(), vmax=m.max())\n    axs[5][1].set_title('Regularized Inverted Model iter y-z')\n    axs[5][1].axis('tight')\n    axs[5][2].imshow(minv3d_reg[:, :, 220].T, cmap='gist_rainbow', vmin=m.min(), vmax=m.max())\n    axs[5][2].set_title('Regularized Inverted Model iter x-y')\n    axs[5][2].axis('tight')"
       ]
     }
   ],

_downloads/475cd62e887dac3f4cecf3ef64657324/plot_distributed_array.py

@@ -71,12 +71,25 @@
 pylops_mpi.plot_local_arrays(arr1, "Distributed Array - 1", vmin=0, vmax=1)
 pylops_mpi.plot_local_arrays(arr2, "Distributed Array - 2", vmin=0, vmax=1)
 
+###############################################################################
+# **Scaling** - Each process operates on its local portion of
+# the array and scales the corresponding elements by a given scalar.
+scale_arr = .5 * arr1
+pylops_mpi.plot_local_arrays(scale_arr, "Scaling", vmin=0, vmax=1)
+
 ###############################################################################
 # **Element-wise Addition** - Each process operates on its local portion of
 # the array and adds the corresponding elements together.
 sum_arr = arr1 + arr2
 pylops_mpi.plot_local_arrays(sum_arr, "Addition", vmin=0, vmax=1)
 
+###############################################################################
+# **Element-wise In-place Addition** - Similar to the previous one but the
+# addition is performed directly on one of the addends without creating a new
+# distributed array.
+sum_arr += arr2
+pylops_mpi.plot_local_arrays(sum_arr, "Addition", vmin=0, vmax=1)
+
 ###############################################################################
 # **Element-wise Subtraction** - Each process operates on its local portion
 # of the array and subtracts the corresponding elements together.

_downloads/6b4bdf99ba907d36440c6747c4b6c693/plot_stacked_array.py

@@ -0,0 +1,155 @@
+"""
+Stacked Array
+=========================
+This example shows how to use the :py:class:`pylops_mpi.StackedDistributedArray`.
+This class provides a way to combine and act on multiple :py:class:`pylops_mpi.DistributedArray`
+within the same program. This is very useful in scenarios where an array can be logically
+divided in subarrays and each of them lends naturally to distribution across multiple processes in
+a parallel computing environment.
+"""
+
+from matplotlib import pyplot as plt
+import numpy as np
+from mpi4py import MPI
+
+import pylops
+import pylops_mpi
+
+plt.close("all")
+np.random.seed(42)
+rank = MPI.COMM_WORLD.Get_rank()
+size = MPI.COMM_WORLD.Get_size()
+
+###############################################################################
+# Let's start by defining two distributed array
+subarr1 = pylops_mpi.DistributedArray(global_shape=size * 10,
+                                      partition=pylops_mpi.Partition.SCATTER,
+                                      axis=0)
+subarr2 = pylops_mpi.DistributedArray(global_shape=size * 4,
+                                      partition=pylops_mpi.Partition.SCATTER,
+                                      axis=0)
+# Filling the local arrays
+subarr1[:], subarr2[:] = 1, 2
+
+###############################################################################
+# We combine them into a single
+# :py:class:`pylops_mpi.StackedDistributedArray` object.
+arr1 = pylops_mpi.StackedDistributedArray([subarr1, subarr2])
+if rank == 0:
+    print('Stacked array:', arr1)
+
+# Extract and print full array
+full_arr1 = arr1.asarray()
+if rank == 0:
+    print('Full array:', full_arr1)
+
+# Modify the part of the first array in rank0
+if rank == 0:
+    arr1[0][:] = 10
+full_arr1 = arr1.asarray()
+if rank == 0:
+    print('Modified full array:', full_arr1)
+
+###############################################################################
+# Let's now create a second :py:class:`pylops_mpi.StackedDistributedArray` object
+# and perform different mathematical operations on those two objects.
+subarr1_ = pylops_mpi.DistributedArray(global_shape=size * 10,
+                                       partition=pylops_mpi.Partition.SCATTER,
+                                       axis=0)
+subarr2_ = pylops_mpi.DistributedArray(global_shape=size * 4,
+                                       partition=pylops_mpi.Partition.SCATTER,
+                                       axis=0)
+# Filling the local arrays
+subarr1_[:], subarr2_[:] = 5, 6
+arr2 = pylops_mpi.StackedDistributedArray([subarr1_, subarr2_])
+if rank == 0:
+    print('Stacked array 2:', arr2)
+
+full_arr2 = arr2.asarray()
+if rank == 0:
+    print('Full array2:', full_arr2)
+
+###############################################################################
+# **Negation**
+neg_arr = -arr1
+full_neg_arr = neg_arr.asarray()
+if rank == 0:
+    print('Negated full array:', full_neg_arr)
+
+###############################################################################
+# **Element-wise Addition**
+sum_arr = arr1 + arr2
+full_sum_arr = sum_arr.asarray()
+if rank == 0:
+    print('Summed full array:', full_sum_arr)
+
+###############################################################################
+# **Element-wise Subtraction**
+sub_arr = arr1 - arr2
+full_sub_arr = sub_arr.asarray()
+if rank == 0:
+    print('Subtracted full array:', full_sub_arr)
+
+###############################################################################
+# **Multiplication**
+mult_arr = arr1 * arr2
+full_mult_arr = mult_arr.asarray()
+if rank == 0:
+    print('Multipled full array:', full_mult_arr)
+
+###############################################################################
+# **Dot-product**
+dot_arr = arr1.dot(arr2)
+if rank == 0:
+    print('Dot-product:', dot_arr)
+    print('Dot-product (np):', np.dot(full_arr1, full_arr2))
+
+###############################################################################
+# **Norms**
+l0norm = arr1.norm(0)
+l1norm = arr1.norm(1)
+l2norm = arr1.norm(2)
+linfnorm = arr1.norm(np.inf)
+
+if rank == 0:
+    print('L0 norm', l0norm, np.linalg.norm(full_arr1, 0))
+    print('L1 norm', l1norm, np.linalg.norm(full_arr1, 1))
+    print('L2 norm', l2norm, np.linalg.norm(full_arr1, 2))
+    print('Linf norm', linfnorm, np.linalg.norm(full_arr1, np.inf))
+
+###############################################################################
+# Now that we have a way to stack multiple :py:class:`pylops_mpi.StackedDistributedArray` objects,
+# let's see how we can apply operators to them. More specifically this can be
+# done using the :py:class:`pylops_mpi.MPIStackedVStack` operator that takes multiple
+# :py:class:`pylops_mpi.MPILinearOperator` objects, each acting on one specific
+# distributed array
+x = pylops_mpi.DistributedArray(global_shape=size * 10,
+                                partition=pylops_mpi.Partition.SCATTER,
+                                axis=0)
+# Filling the local arrays
+x[:] = 1.
+
+# Make stacked operator
+mop1 = pylops_mpi.MPIBlockDiag([pylops.MatrixMult(np.ones((5, 10))), ])
+mop2 = pylops_mpi.MPIBlockDiag([pylops.MatrixMult(2 * np.ones((8, 10))), ])
+mop = pylops_mpi.MPIStackedVStack([mop1, mop2])
+
+y = mop.matvec(x)
+y_arr = y.asarray()
+xadj = mop.rmatvec(y)
+xadj_arr = xadj.asarray()
+
+if rank == 0:
+    print('StackedVStack y', y, y_arr, y_arr.shape)
+    print('StackedVStack xadj', xadj, xadj_arr, xadj_arr.shape)
+
+###############################################################################
+# Finally, let's solve now an inverse problem using stacked arrays instead
+# of distributed arrays
+x0 = x.copy()
+x0[:] = 0.
+xinv = pylops_mpi.cgls(mop, y, x0=x0, niter=15, tol=1e-10, show=False)[0]
+xinv_array = xinv.asarray()
+
+if rank == 0:
+    print('xinv_array', xinv_array)