{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# L1-L1 IRLS\n\nThis example shows how to use the :py:class:`pylops.optimization.sparsity.irls` solver to\nsolve problems in the form:\n\n    .. math::\n        J = \\left\\| \\mathbf{y}-\\mathbf{Ax}\\right\\|_{1} + \\epsilon \\left\\|\\mathbf{x}\\right\\|_{1}\n\nThis can be easily achieved by recasting the problem into this equivalent formulation:\n\n    .. math::\n        J = \\left\\|\\left[\\begin{array}{c}\n        \\mathbf{A} \\\\\n        \\epsilon \\mathbf{I}\n        \\end{array}\\right] \\mathbf{x}-\\left[\\begin{array}{l}\n        \\mathbf{y} \\\\\n        \\mathbf{0}\n        \\end{array}\\right]\\right\\|_{1}\n\nand solving it using the classical version of the IRLS solver with L1 norm on the data term. In PyLops,\nthe creation of the augmented system happens under the hood when users provide the following optional\nparameter (``kind=\"datamodel\"``) to the solver.\n\nWe will now consider a 1D deconvolution problem where the signal is contaminated with Laplace noise.\nWe will compare the classical L2-L1 IRLS solver that works optimally under the condition of Gaussian\nnoise with the above descrived L1-L1 IRLS solver that is best suited to the case of Laplace noise.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import random\n\nimport matplotlib.pyplot as plt\nimport numpy as np\n\nimport pylops\n\nplt.close(\"all\")\nnp.random.seed(10)\nrandom.seed(0)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Let's start by creating a spiky input signal and convolving it with a Ricker\nwavelet.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "dt = 0.004\nnt = 201\nt = np.arange(nt) * dt\n\nnspikes = 5\nx = np.zeros(nt)\nx[random.sample(range(0, nt - 1), nspikes)] = -1 + 2 * np.random.rand(nspikes)\n\nh, th, hcenter = pylops.utils.wavelets.ricker(t[:101], f0=20)\nCop = pylops.signalprocessing.Convolve1D(nt, h=h, offset=hcenter)\n\ny = Cop @ x"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We add now a realization of Laplace-distributed noise to our signal and\nperform a standard spiky deconvolution\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "yn = y + np.random.laplace(loc=0.0, scale=0.05, size=y.shape)\n\nxl2l1 = pylops.optimization.sparsity.irls(\n    Cop,\n    yn,\n    threshR=True,\n    kind=\"model\",\n    nouter=100,\n    epsR=1e-4,\n    epsI=1.0,\n    warm=True,\n    **dict(iter_lim=100),\n)[0]\n\nxl1l1 = pylops.optimization.sparsity.irls(\n    Cop,\n    yn,\n    threshR=True,\n    kind=\"datamodel\",\n    nouter=100,\n    epsR=1e-4,\n    epsI=1.0,\n    warm=True,\n    **dict(iter_lim=100),\n)[0]\n\nfig, axs = plt.subplots(2, 1, sharex=True, figsize=(12, 5))\naxs[0].plot(t, y, \"k\", lw=4, label=\"Clean\")\naxs[0].plot(t, yn, \"r\", lw=2, label=\"Noisy\")\naxs[0].legend()\naxs[0].set_title(\"Data\")\naxs[1].plot(t, x, \"k\", lw=4, label=\"L2-L1\")\naxs[1].plot(\n    t,\n    xl2l1,\n    \"r\",\n    lw=2,\n    label=f\"L2-L1 (NMSE={(np.linalg.norm(xl2l1 - x)/np.linalg.norm(x)):.2f})\",\n)\naxs[1].plot(\n    t,\n    xl1l1,\n    \"c\",\n    lw=2,\n    label=f\"L1-L1 (NMSE={(np.linalg.norm(xl1l1 - x)/np.linalg.norm(x)):.2f})\",\n)\naxs[1].legend()\naxs[1].set_xlabel(\"t\")\nplt.tight_layout()"
      ]
    }
  ],
  "metadata": {
    "kernelspec": {
      "display_name": "Python 3",
      "language": "python",
      "name": "python3"
    },
    "language_info": {
      "codemirror_mode": {
        "name": "ipython",
        "version": 3
      },
      "file_extension": ".py",
      "mimetype": "text/x-python",
      "name": "python",
      "nbconvert_exporter": "python",
      "pygments_lexer": "ipython3",
      "version": "3.9.15"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}