Finalized first version of notebooks

DIG-Kaust · Mar 31, 2024 · 442f183 · 442f183
1 parent fe29e7c
commit 442f183
Show file tree

Hide file tree

Showing 7 changed files with 634 additions and 867 deletions.
diff --git a/README.md b/README.md
@@ -1,22 +1,25 @@
 ![LOGO](https://github.com/DIG-Kaust/Project_Template/blob/master/logo.png)
 
-Reproducible material for **XXX -
-Ravasi M., Author M., Author C.** submitted to XXX.
+This library contains all the required building blocks to perform Full Waveform Inversion (FWI) with Devito (for modelling). 
+The vision is to seamlessly integrate this with PyLops and PyProximal to benefit from their modular handling of linear 
+and proximal operators.
 
 
 ## Project structure
 This repository is organized as follows:
 
-* :open_file_folder: **package**: python library containing routines for ....;
-* :open_file_folder: **data**: folder containing data (or instructions on how to retrieve the data
-* :open_file_folder: **notebooks**: set of jupyter notebooks reproducing the experiments in the paper (see below for more details);
-* :open_file_folder: **scripts**: set of python scripts used to run multiple experiments ...
+* :open_file_folder: **devitofwi**: python library containing routines to perform FWI with Devito modelling engines;
+* :open_file_folder: **data**: folder containing sample data to run notebooks;
+* :open_file_folder: **notebooks**: set of jupyter notebooks used to showcase different ways to run FWI with ``devitofwi``;
 
 ## Notebooks
 The following notebooks are provided:
 
-- :orange_book: ``X1.ipynb``: notebook performing ...;
-- :orange_book: ``X2.ipynb``: notebook performing ...
+- :orange_book: ``Modelling_filtering.ipynb``: notebook comparing modelling with a filtered wavelet and filtering of original data computed with unfiltered wavelet;
+- :orange_book: ``AcousticVel_L2_1stage.ipynb``: notebook performing acoustic FWI parametrized in velocity with entire data;
+- :orange_book: ``AcousticVel_L2refr_1stage.ipynb``: notebook performing acoustic FWI parametrized in velocity with only refracted waves;
+- :orange_book: ``AcousticVel_L2_Nstages.ipynb``: notebook performing acoustic FWI parametrized in velocity with entire data in N frequency stages;
+- :orange_book: ``AcousticVel_L2refr_Nstages.ipynb``: notebook performing acoustic FWI parametrized in velocity with only refracted waves in N frequency stages.
 
 
 ## Getting started :space_invader: :robot:
@@ -37,13 +40,15 @@ pip install -e .
 
 Remember to always activate the environment by typing:
 ```
-conda activate my_env
+conda activate devitofwi
 ```
 
-Finally, to run tests simply type:
-```
-pytest
-```
+## To do list :memo:
+
+The following list is intended to define some of the improvements that should be implemented moving forward:
 
-**Disclaimer:** All experiments have been carried on a Intel(R) Xeon(R) CPU @ 2.10GHz equipped with a single NVIDIA GEForce RTX 3090 GPU. Different environment 
-configurations may be required for different combinations of workstation and GPU.
+- [ ] Optimize time step in N stages FWI to avoid using a too large time step when working at low-frequency 
+      (not chosen from the velocity model and frequency of the observed data  but from that of each stage)
+- [ ] Create ``acousticfwi`` wrapper that can implement multi-stage acoustic FWI with time-space masking
+      (and can specialize to any alternative case)
+- [ ] Improve handling of user-defined source wavelets
diff --git a/devitofwi/preproc/masking.py b/devitofwi/preproc/masking.py
@@ -1,18 +1,23 @@
 import numpy as np
 import matplotlib.pyplot as plt
-from scipy.ndimage import gaussian_filter
+from scipy.signal import filtfilt
 
 
-def timespace_mask(arrival, nt, dt, toff=0., smooth=5):
-    """Design time-space mask based on arrival time
+def timespace_mask(arrival, nt, dt, toff=0., nsmooth=None):
+    """Design time-space mask based on arrival times
 
     """
     ns, nr = arrival.shape
 
+    # create mask based on arrival times
     mask = np.zeros((ns, nr, nt))
     for isrc in range(ns):
         for irec in range(nr):
             it = np.round((arrival[isrc, irec] + toff) / dt).astype('int')
             mask[isrc, irec, it:] = 1.
 
+    # smooth mask along time axis
+    if nsmooth is not None:
+        mask = filtfilt(np.ones(nsmooth) / nsmooth, 1, mask, axis=-1)
+
     return mask
diff --git a/notebooks/acoustic/AcousticVel_L2_1stage.ipynb b/notebooks/acoustic/AcousticVel_L2_1stage.ipynb
diff --git a/notebooks/acoustic/AcousticVel_L2_Nstages.ipynb b/notebooks/acoustic/AcousticVel_L2_Nstages.ipynb
@@ -7,7 +7,11 @@
     "id": "_vT1QXwvePfu"
    },
    "source": [
-    "# Acoustic Multi-stage FWI parametrized in velocity using full data"
+    "# Multi-stage Acoustic FWI(VP) with entire data\n",
+    "\n",
+    "## Author: M. Ravasi\n",
+    "\n",
+    "This notebook performs acoustic FWI parametrized in velocity using the entire dataset and multiple stages where lower frequency components are inverted first and higher frequency components are included gradually in the process. "
    ]
   },
   {
@@ -53,8 +57,7 @@
     "from devitofwi.waveengine.acoustic import AcousticWave2D\n",
     "from devitofwi.loss.l2 import L2\n",
     "\n",
-    "configuration['log-level'] = 'ERROR'\n",
-    "#configuration['log-level'] = 'WARNING'"
+    "configuration['log-level'] = 'ERROR'"
    ]
   },
   {
@@ -67,6 +70,18 @@
     "clear_devito_cache()"
    ]
   },
+  {
+   "cell_type": "code",
+   "execution_count": 13,
+   "id": "3ba40d47",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Callback to track model error\n",
+    "def fwi_callback(xk, vp, vp_error):\n",
+    "    vp_error.append(np.linalg.norm((xk - vp.reshape(-1))/vp.reshape(-1)))"
+   ]
+  },
   {
    "cell_type": "markdown",
    "id": "BgkMPqmyePfx",
@@ -115,7 +130,7 @@
     "\n",
     "# Velocity model\n",
     "path = '../../data/'\n",
-    "velocity_file = path + 'Marm.bin' # true model "
+    "velocity_file = path + 'Marm.bin'"
    ]
   },
   {
@@ -362,7 +377,7 @@
   {
    "cell_type": "code",
    "execution_count": 9,
-   "id": "373e9ab9",
+   "id": "d76d067a",
    "metadata": {},
    "outputs": [
     {
@@ -401,7 +416,7 @@
   {
    "cell_type": "code",
    "execution_count": 10,
-   "id": "a5d59942",
+   "id": "941c87f4",
    "metadata": {},
    "outputs": [],
    "source": [
@@ -425,7 +440,7 @@
   {
    "cell_type": "code",
    "execution_count": 11,
-   "id": "4b5586eb",
+   "id": "58dee018",
    "metadata": {},
    "outputs": [
     {
@@ -459,7 +474,7 @@
   {
    "cell_type": "code",
    "execution_count": 12,
-   "id": "e93c3a43",
+   "id": "e30900d6",
    "metadata": {},
    "outputs": [
     {
@@ -485,18 +500,6 @@
     "              vmin=-d_vmax, vmax=d_vmax)"
    ]
   },
-  {
-   "cell_type": "code",
-   "execution_count": 13,
-   "id": "fdb5c93d",
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "# Callback to track model error\n",
-    "def fwi_callback(xk, vp, vp_error):\n",
-    "    vp_error.append(np.linalg.norm((xk - vp.reshape(-1))/vp.reshape(-1)))"
-   ]
-  },
   {
    "cell_type": "code",
    "execution_count": 14,
@@ -4395,7 +4398,7 @@
   {
    "cell_type": "code",
    "execution_count": 15,
-   "id": "96bca601",
+   "id": "b23ae1ff",
    "metadata": {},
    "outputs": [
     {
@@ -4436,7 +4439,7 @@
   {
    "cell_type": "code",
    "execution_count": 16,
-   "id": "518ff9c5",
+   "id": "b0d555b9",
    "metadata": {},
    "outputs": [
     {
@@ -4469,7 +4472,7 @@
   {
    "cell_type": "code",
    "execution_count": 17,
-   "id": "2f8debf9",
+   "id": "9bbea043",
    "metadata": {},
    "outputs": [
     {

diff --git a/notebooks/acoustic/AcousticVel_L2refr_1stage.ipynb b/notebooks/acoustic/AcousticVel_L2refr_1stage.ipynb
diff --git a/notebooks/acoustic/AcousticVel_L2refr_Nstages.ipynb b/notebooks/acoustic/AcousticVel_L2refr_Nstages.ipynb
@@ -7,7 +7,13 @@
     "id": "_vT1QXwvePfu"
    },
    "source": [
-    "# Acoustic Multi-stage FWI parametrized in velocity using only refractions"
+    "# Multi-stage Acoustic FWI(VP)  using first refractions and then entire data\n",
+    "\n",
+    "## Author: M. Ravasi\n",
+    "\n",
+    "This notebook performs acoustic FWI parametrized in velocity using the entire dataset and multiple stages where lower frequency components are inverted first and higher frequency components are included gradually in the process. \n",
+    "\n",
+    "Here, only refractions are inverted for the low-frequency components, whilst the entire data is inverted in the last stage with higher frequency components. This notebook showcases how to combine the filtering capabilities of ``devitofwi`` provided by ``devitofwi.preproc.filtering.Filter``  with the masking capabilities of ``devitofwi.preproc.masking.timespace_mask``."
    ]
   },
   {
@@ -394,7 +400,7 @@
   {
    "cell_type": "code",
    "execution_count": 10,
-   "id": "21e516c2",
+   "id": "b57fc4f7",
    "metadata": {
     "scrolled": false
    },
@@ -443,7 +449,7 @@
     "distance = np.sqrt((x_s[:, 0][:, None]-x_r[:, 0][None, :])**2 + (x_s[:, 1][:, None]-x_r[:, 1][None, :])**2)\n",
     "directtime = distance / vwater\n",
     "\n",
-    "tsmask = 1 - timespace_mask(directtime, amod.geometry.nt, amod.geometry.dt * 1e-3, toff, vwater)\n",
+    "tsmask = 1 - timespace_mask(directtime, amod.geometry.nt, amod.geometry.dt * 1e-3, toff, nsmooth=None)\n",
     "\n",
     "fig, axs = plt.subplots(1, 3, sharey=True, figsize=(14, 9))\n",
     "for ax, ishot in zip(axs, [0, par['ns']//2, par['ns']-1]):\n",
@@ -477,7 +483,7 @@
   {
    "cell_type": "code",
    "execution_count": 11,
-   "id": "3915f422",
+   "id": "a1c17de7",
    "metadata": {},
    "outputs": [],
    "source": [
@@ -492,7 +498,7 @@
   {
    "cell_type": "code",
    "execution_count": 12,
-   "id": "2801f374",
+   "id": "3d1aec79",
    "metadata": {
     "scrolled": false
    },
@@ -530,7 +536,7 @@
   {
    "cell_type": "code",
    "execution_count": 13,
-   "id": "8f0ed489",
+   "id": "d7036b79",
    "metadata": {},
    "outputs": [
     {
@@ -558,7 +564,7 @@
   },
   {
    "cell_type": "raw",
-   "id": "85a9201a",
+   "id": "b6bbbef4",
    "metadata": {},
    "source": [
     "# Testing first gradient\n",
@@ -590,7 +596,7 @@
     "              vmin=-d_vmax, vmax=d_vmax)\n",
     "\n",
     "# Mask observed data to have direct arrival and refractions\n",
-    "tsmask = 1 - timespace_mask(directtime, amod.geometry.nt, amod.geometry.dt * 1e-3, toff[ifreq], vwater)\n",
+    "tsmask = 1 - timespace_mask(directtime, amod.geometry.nt, amod.geometry.dt * 1e-3, toff[ifreq], nsmooth=None)\n",
     "TSmaskop = [Diagonal(tsmask[isrc].T.ravel()) for isrc in range(par['ns'])]\n",
     "dobs_masked = np.hstack([(TSmaskop[isrc] @ dobsfilt[isrc].ravel())[None, :] for isrc in range(par['ns'])])\n",
     "dobs_masked = dobs_masked.reshape(par['ns'], amod.geometry.nt, par['nr'])\n",
@@ -630,7 +636,7 @@
   },
   {
    "cell_type": "raw",
-   "id": "c9ed32c3",
+   "id": "131aaa05",
    "metadata": {},
    "source": [
     "# Compare with 7Hz ricker\n",
@@ -661,7 +667,7 @@
     "                      space_order=space_order, nbl=nbl)\n",
     "\n",
     "\n",
-    "tsmask = 1 - timespace_mask(directtime, amod1.geometry.nt, amod1.geometry.dt * 1e-3, toff, vwater)\n",
+    "tsmask = 1 - timespace_mask(directtime, amod1.geometry.nt, amod1.geometry.dt * 1e-3, toff, nsmooth=None)\n",
     "\n",
     "# Mask observed data to have direct arrival and refractions\n",
     "TSmaskop = [Diagonal(tsmask[isrc].T.ravel()) for isrc in range(par['ns'])]\n",
@@ -4548,7 +4554,7 @@
     "    # Mask observed data to have direct arrival and refractions\n",
     "    if applymask[ifreq]:\n",
     "        tsmask = 1 - timespace_mask(directtime, amod.geometry.nt, amod.geometry.dt * 1e-3, \n",
-    "                                    wavpad * amod.geometry.dt * 1e-3 + toff[ifreq], vwater)\n",
+    "                                    wavpad * amod.geometry.dt * 1e-3 + toff[ifreq], nsmooth=None)\n",
     "        TSmaskop = [Diagonal(tsmask[isrc].T.ravel()) for isrc in range(par['ns'])]\n",
     "        dobs_masked = np.hstack([(TSmaskop[isrc] @ dobs_filt[isrc].ravel())[None, :] for isrc in range(par['ns'])])\n",
     "        dobs_masked = dobs_masked.reshape(par['ns'], amod.geometry.nt, par['nr'])\n",

diff --git a/notebooks/acoustic/Modelling_filtering.ipynb b/notebooks/acoustic/Modelling_filtering.ipynb