{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "b0e50c8f",
   "metadata": {},
   "source": [
    "# Deconvolution Problems \n",
    "\n",
    "We consider the two-dimensional deconvolution problems to find a non-negative function f given data \n",
    "$$\n",
    "    d \\sim \\mathrm{Pois}(h*f)\n",
    "$$\n",
    "with a non-negative convolution kernel $h$, and $\\mathrm{Pois}$ denotes the element-wise Poisson distribution.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9b627178",
   "metadata": {},
   "outputs": [],
   "source": [
    "import sys\n",
    "from pathlib import Path\n",
    "\n",
    "import numpy as np\n",
    "\n",
    "from regpy.operators import GaussianBlur\n",
    "from regpy.solvers import Setting\n",
    "from regpy.solvers.linear import TikhonovCG\n",
    "from regpy.functionals import *\n",
    "from regpy.hilbert import L2\n",
    "\n",
    "\n",
    "example_dir = \"../../../../examples/deconvolution\"\n",
    "Path(example_dir).resolve()\n",
    "sys.path.insert(0, str(example_dir))\n",
    "from plots import comparison_plot, convergence_plot, plot_recos\n",
    "from test_images import mixed\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d1dc8a11",
   "metadata": {},
   "source": [
    "## Generate synthetic data with impulsive noise"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "485d1481",
   "metadata": {},
   "outputs": [],
   "source": [
    "grid, exact_sol = mixed(M=256,N=256,c_bubbles=1.,fac=50.)\n",
    "r\"\"\"grid is the underlying UniformGridFcts vector space, and exact_sol the exact solution.\"\"\"\n",
    "a=0.05\n",
    "conv =  GaussianBlur(grid,a,pad_amount=16)\n",
    "r\"\"\"Convolution operator \\(f\\mapsto h*f\\) for the convolution kernel \\(h(x)=\\exp(-|x|_2^2/a^2)\\).\"\"\"\n",
    "blur = conv(exact_sol)\n",
    "blur[blur<0] = 0.\n",
    "data = blur.copy()\n",
    "n,m=grid.shape\n",
    "for j in range(64**2):\n",
    "    nn = np.random.randint(n)\n",
    "    mm = np.random.randint(m)\n",
    "    data[nn,mm] =  np.random.randint(2000)-1000\n",
    "\n",
    "comparison_plot(grid,exact_sol,data,title_left='noisy measurement data')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "28688e65",
   "metadata": {},
   "source": [
    "## Reconstruction with a quadratic data fidelity term\n",
    "... fails completely for this data!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3b58f161",
   "metadata": {},
   "outputs": [],
   "source": [
    "alpha = 1e-4\n",
    "setting = Setting(conv,L2,L2)\n",
    "solver = TikhonovCG(setting=setting, data = data, regpar = alpha,reltolx=1e-6)\n",
    "fal,_ = solver.run()\n",
    "\n",
    "comparison_plot(grid,exact_sol,fal,title_left='quadratic Tikhonov regularization')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d54633ef",
   "metadata": {},
   "source": [
    "## setting with Huber data fidelity and quadratic penalty term\n",
    "We now set up a generalized Tikhonov regularization setting\n",
    "$$\n",
    "\\hat{f} \\in \\mathrm{argmin}_f\\left[\\frac{1}{\\alpha}H_{\\sigma}(Tf-g^{\\mathrm{obs}})+\\|f\\|^2\\right]    \n",
    "$$\n",
    "where the Huber data fidelity term $H_{\\sigma}$ is quadratic on coordinates in $[-\\sigma,\\sigma]$ and $|\\cdot|-\\sigma/2$ for large coordinates.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b831ccb9",
   "metadata": {},
   "outputs": [],
   "source": [
    "from regpy.functionals import Huber,LppPower\n",
    "sigma = 0.1\n",
    "Sdat = (1./sigma)*Huber(grid,sigma=sigma,eps =1e-10)\n",
    "L2pen = LppPower(grid,p=2)\n",
    "huber_setting = Setting(conv,L2pen,Sdat,alpha,data=data)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fb4ec86a",
   "metadata": {},
   "source": [
    "display the available methods"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cf0f8924",
   "metadata": {},
   "outputs": [],
   "source": [
    "huber_setting.display_all_methods(full_names=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bc717299",
   "metadata": {},
   "outputs": [],
   "source": [
    "methods =  ['FISTA','dual_FISTA','PDHG','dual_PDHG','dual_FB']\n",
    "\n",
    "recos = {}\n",
    "for method in methods:\n",
    "    #huber_setting.set_stopping_rule(method, OptimalityCondStopping(huber_setting,tol=0.01,max_iter=500,logging_level=logging.INFO))\n",
    "    recos[method]= huber_setting.run(method)[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "58bd3402",
   "metadata": {},
   "source": [
    "plot the convergence histories "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b93001bd",
   "metadata": {},
   "outputs": [],
   "source": [
    "convergence_plot(huber_setting,methods,'Huber data fid., quadratic penalty')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a74baafe",
   "metadata": {},
   "source": [
    "visualize the reconstructions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "52d58a3f",
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_recos(huber_setting,methods,recos,truth=exact_sol)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "RegPy_testing",
   "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.12.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}