compressive_imaging.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Compressive Imaging with Deep Untrained Decoder Prior"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "num GPUs 1\n"
     ]
    }
   ],
   "source": [
    "from __future__ import print_function\n",
    "import matplotlib.pyplot as plt\n",
    "import os\n",
    "import warnings\n",
    "warnings.filterwarnings('ignore')\n",
    "from include import *\n",
    "from PIL import Image\n",
    "import PIL\n",
    "import pywt\n",
    "import numpy as np\n",
    "import torch\n",
    "import torch.optim\n",
    "from torch.autograd import Variable\n",
    "from sklearn import linear_model\n",
    "\n",
    "GPU = True\n",
    "if GPU == True:\n",
    "    torch.backends.cudnn.enabled = True\n",
    "    torch.backends.cudnn.benchmark = True\n",
    "    dtype = torch.cuda.FloatTensor\n",
    "    os.environ['CUDA_VISIBLE_DEVICES'] = '0'\n",
    "    print(\"num GPUs\",torch.cuda.device_count())\n",
    "    device = 'cuda'\n",
    "    if torch.cuda.device_count()==0:\n",
    "        dtype = torch.FloatTensor\n",
    "        device = 'cpu'\n",
    "else:\n",
    "    dtype = torch.FloatTensor\n",
    "    device = 'cpu'\n",
    "from scipy.fftpack import dct\n",
    "from scipy.fftpack import idct\n",
    "from scipy import io as sio \n",
    "import time"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Load image and pre-process"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Dimensions of input image: (1, 28, 28)\n"
     ]
    }
   ],
   "source": [
    "#dataset = 'mnist' # 'mnist' or 'celeba'\n",
    "dataset = 'mnist'\n",
    "path = './test_data/' + dataset + '/' \n",
    "img_name = dataset + '1' # 1-5 (for celeba), 1-6 (for mnist)\n",
    "img_path = path + img_name + \".jpg\"\n",
    "img_pil = Image.open(img_path)\n",
    "if dataset == 'celeba':\n",
    "    img_pil = img_pil.crop((60,80+20,60+64,80+84)) #crop to 3 x 64 x 64\n",
    "img_np = pil_to_np(img_pil)\n",
    "print('Dimensions of input image:', img_np.shape)\n",
    "img_np = img_np / np.max(img_np)\n",
    "img_np_orig = 1*img_np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Display image x and convert to pytorch variable"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPsAAAD4CAYAAAAq5pAIAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAACrdJREFUeJzt3ctLVmsfxvH1mJqHMstTJ9rSIEIoCoKGHSCoCKOIiAYRDSSIoPoDiogGTZpUo5w0iUYRBB0oEkQqGqTQJCgqTMKOmloeMn1H74YXWtfvzZWPz97X9zPs4nY9PXbtBfu37nvlJicnEwD/fkUz/QEA5AdlB0xQdsAEZQdMUHbARHE+L5bL5fhf/79QVKT/mzsxMSHzsrKy1GxkZGRKn+lPqaysTM2+ffuWx0/ye+bOnSvzwcHBPH2S3zc5OZn71Z9zZwdMUHbABGUHTFB2wARlB0xQdsAEZQdM5PK56811zq7m4Ekys7Pw8vJymY+Pj8v8x48ff/Lj/I+KigqZj46Oylx979GMf9asWTL/+fOnzGcSc3bAHGUHTFB2wARlB0xQdsAEZQdMUHbABHP2AjBnzhyZR7PssbGx1Cz6/UbXHhoaknkWxcX6OIVc7pfj4r9Fs251DkBDQ4Nc+/79e5kX8hyeOTtgjrIDJig7YIKyAyYoO2CCsgMm8nqUNH4t2qoZjd6yHCUdjbciJSUlMlefPTquua+vT+bR9tzh4eHULBqtVVdXy7y/v1/mhYg7O2CCsgMmKDtggrIDJig7YIKyAyYoO2CCOXsBKC0tlXk0Z1fz5qdPn8q1jx8/lnm0xXXt2rUyP3nyZGrW1tYm186ePVvmao6eJEkyf/781Cya4UfPJxTyFtc03NkBE5QdMEHZAROUHTBB2QETlB0wQdkBExwlnQdZ90Y3NjbKvL29PTWrqamRa6PXIkfUMdZJkiSvXr1KzQ4ePCjXdnV1yTw6B6CoKP1eVllZKdcODg7KPHoGIPps04mjpAFzlB0wQdkBE5QdMEHZAROUHTBB2QETzNkLQDSzvXnzpsy3bNmSmmV95XL0WuXu7m6Zr1ixIjV7+/atXLt9+3aZv379WuZqT3q037yurk7mHz9+lPlMYs4OmKPsgAnKDpig7IAJyg6YoOyACUZveRCNr86cOSPz48ePy1yNgVpaWuTa5uZmmV+6dEnm0Rbaa9eupWYLFy6Ua6PxVn19vcyrqqpSs4GBAblWbY9NkvhV1zN5lDSjN8AcZQdMUHbABGUHTFB2wARlB0xQdsDEv2bOXlJSIvPotccR9YreaKa6dOlSmb9582YqH+lvra2tqdnhw4flWvW65ySJX4sczaOXLVuWmnV2dsq10RHc9+/fl/np06dTs46ODrk2mqNHz05k/feWBXN2wBxlB0xQdsAEZQdMUHbABGUHTFB2wERe5+yzZs2SF5uYmJDr1Uw365w9urb6+dHPvnHjhsw3btwo8+j1wn/99Vdq1tvbK9eq5weSJP67RfNm9TvbvHmzXHv79m2ZR6+6Vr/TNWvWyLXRXnp1TPVMY84OmKPsgAnKDpig7IAJyg6YoOyACcoOmNBD0j8smmVnMTY2JvPoeYJoX/b4+Phvf6b/WrlypczV+eZJEu+9fvfuXWoW/b2iOXppaanMo7386vp37tyRa0+cOCHz8+fPy1zZv3+/zC9evCjz6PmCLP9epgt3dsAEZQdMUHbABGUHTFB2wARlB0xQdsDEP+rceDXbzDrXzDI3jfbS3717V+abNm2S+eXLl2V+7Nix1Oz79+9ybVbRHD/Kleh3Gs34BwcHU7PoXPi9e/fK/N69ezKfzmdKIuxnB8xRdsAEZQdMUHbABGUHTFB2wERet7hGojGNGnFFY5royORofKbGPLW1tXLt6tWrZT40NCTzaJupGq9F3+m8efNk3tfXJ/PIdG71PHTokMzVq6yjUeuuXbtkHo1TCxF3dsAEZQdMUHbABGUHTFB2wARlB0xQdsBEXufs0bbCaFtglm2D0VbeaLukWl9fXy/X1tTUyDzS1dU15bXRd/b169cp/+wkib/XLK+6jp4BuHLlisyPHDmSmkXPPkTfW3l5ucyHh4dlPhO4swMmKDtggrIDJig7YIKyAyYoO2CCsgMmCmo/eySahSvR3DSa+Sp1dXUyHxgYkHlZWZnMr1+/LvPq6urUrL+/X66NvpeKigqZR0dVZ/leo9dwR2cQqGcIZs+eLdd+/vxZ5qOjozIvRNzZAROUHTBB2QETlB0wQdkBE5QdMEHZARN5nbNHc9Forprl9b/RufFZXl09f/58mVdVVcn80aNHMu/p6ZG5mpVH33kk6yuf1RkG0Xce7QmPnk9Q5whEz2wsX75c5jP5Suap4s4OmKDsgAnKDpig7IAJyg6YoOyACcoOmMjrnD2ao0d7jLPsIY7mzSMjI1P+2Tt27JB59Pdev369zNetWydzda58dO25c+fKPNqPHs261fvZs57l39TUJPPa2trULHruoru7W+bROxCyPLcxXbizAyYoO2CCsgMmKDtggrIDJig7YKKgjpKezm2DWY40TpIkKS5O/6q+fPki10Zbc6N80aJFMn/y5ElqFo2Yom2k0Ygpy8gyEm0NvnXrlswbGhpSs5cvX8q1N2/elHkhjtYi3NkBE5QdMEHZAROUHTBB2QETlB0wQdkBE3mds0czX7UdMkmyHUscbZeM5slKZ2enzNWMPkmSpLe3V+Zbt26VeUdHR2oWvXo4qyxbXEtLS+Xac+fOyXzBggUyV1pbW2X+8OHDKf/sQsWdHTBB2QETlB0wQdkBE5QdMEHZAROUHTCR1zl7NOuOqKOmo73w0X72aE6v5sXt7e1y7cDAgMwXLlw45WsnyfTO0hsbG2X+5s2bKa+/cOGCXBsd0R0dLX7q1KnULJrhR8+ERM8IROcEzATu7IAJyg6YoOyACcoOmKDsgAnKDpig7ICJXD7Pvy4qKpIXiz6L2jsdzaKjfDpF55tv27ZN5s+ePZP5nj17UrOenh659vv37zKPHD16VOZnz55NzQYHB+XaxYsXy/zBgwcyb25uTs2i5zIKcU7+/5qcnPzl4Qzc2QETlB0wQdkBE5QdMEHZAROUHTCR1y2u0ZHK0TZUNS6Z7tFaeXl5ahaNafr7+2U+NDQk81WrVslcHWV9/fp1ufbFixcy3717t8yXLFkic7UVNBqtXb16VeYHDhyQuRJtt66vr5f5hw8fpnztmcKdHTBB2QETlB0wQdkBE5QdMEHZAROUHTCR1y2uxcXF8mLR7FPN6bO87jlJ4u21c+bMSc2iOXl07ba2Nplv2LBB5p8+fUrNamtr5dpI9L329fXJ/MuXL6nZvn375Nrnz5/LPNqmOjY2JnOloqJC5lm3Bk8ntrgC5ig7YIKyAyYoO2CCsgMmKDtggrIDJvI6Z8/lcpkupvZGRzPV6BW70fqSkpLULNqHH2lqapL5zp07Zd7S0pKaRa9cjo5zjs4gUNdOkiS5du1aaha9Fjl6JXOkuro6NYvOGIgUFen7ZPQMwHRizg6Yo+yACcoOmKDsgAnKDpig7IAJyg6Y+EfN2QHEmLMD5ig7YIKyAyYoO2CCsgMmKDtggrIDJig7YIKyAyYoO2CCsgMmKDtggrIDJig7YIKyAyYoO2CCsgMmKDtggrIDJig7YIKyAyYoO2CCsgMmKDtggrIDJig7YIKyAyYoO2CCsgMmKDtgIq+vbAYwc7izAyYoO2CCsgMmKDtggrIDJig7YIKyAyYoO2CCsgMmKDtggrIDJig7YIKyAyYoO2CCsgMmKDtggrIDJig7YIKyAyYoO2CCsgMmKDtggrIDJv4DVngh/+1u1yoAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "if dataset == 'celeba':\n",
    "    plt.imshow(img_np.transpose(1,2,0))\n",
    "else:\n",
    "    plt.imshow(img_np[0,:,:])\n",
    "    plt.gray()\n",
    "plt.axis('off')\n",
    "img_var = np_to_var(img_np).type(dtype)\n",
    "d = img_np.shape[1]\n",
    "out_ch = img_np.shape[0]\n",
    "d_image = img_np.size"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Select and set up model to run - denoise, CS, PR "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Compression/denoising mode\n"
     ]
    }
   ],
   "source": [
    "#choose mode, 1 (denoising) , 2 (compressed sensing), 3 (phase retrieval)\n",
    "mode = 1\n",
    "if mode==1:\n",
    "    f = 1 #(default)\n",
    "    print('Compression/denoising mode')\n",
    "    Ameas_var = 1\n",
    "    img_var_meas = img_var\n",
    "elif mode==2:\n",
    "    print('Compressed sensing mode')\n",
    "    f = 0.2 #compression rate\n",
    "    print('Compression rate is ', f)\n",
    "    m_image = int(f*d_image)\n",
    "    print('Number of measurements is ',m_image, ' for signal of length ', d_image)\n",
    "    # random Gaussian measurement matrix : A\n",
    "    Ameas = np.random.randn(m_image,d_image).astype(float)/np.sqrt(m_image)\n",
    "    Ameas_var = torch.from_numpy(Ameas).float().to(device)\n",
    "    # measurements : y = A*x\n",
    "    img_var_meas = torch.matmul(Ameas_var,img_var.to(device).reshape(d_image,1))\n",
    "if mode==3:\n",
    "    print('Compressed phase retrieval mode')    \n",
    "    f = 0.6 #compression rate\n",
    "    print('Compression rate is ', f)\n",
    "    m_image = int(f*d_image)\n",
    "    print('Number of measurements is ',m_image, ' for signal of length ', d_image)\n",
    "    # random Gaussian measurement matrix : A    \n",
    "    Ameas = np.random.randn(m_image,d_image).astype(float)/np.sqrt(m_image)\n",
    "    Ameas_var = torch.from_numpy(Ameas).float().to(device)\n",
    "    # full measurements : A*x\n",
    "    img_var_meas = torch.matmul(Ameas_var,img_var.to(device).reshape(d_image,1))      \n",
    "    # absolute valued measurements : y = |A*x|\n",
    "    img_var_meas = torch.abs(img_var_meas)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Define network"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {},
   "outputs": [],
   "source": [
    "#use decoder architecture or DC GAN architecture\n",
    "decodetype = 'upsample' # transposeconv / upsample"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "n_scales= 3 num_channels_up= [25, 15, 10]\n",
      "number of parameters:  1800\n",
      "Sequential(\n",
      "  (dconv0): Sequential(\n",
      "    (0): ReflectionPad2d((0, 0, 0, 0))\n",
      "    (1): Conv2d(25, 15, kernel_size=(1, 1), stride=(1, 1), bias=False)\n",
      "  )\n",
      "  (drelu0): ReLU()\n",
      "  (dbn0): BatchNorm2d(15, eps=1e-05, momentum=0.1, affine=True, track_running_stats=True)\n",
      "  (dups0): Upsample(scale_factor=2, mode=bilinear)\n",
      "  (dconv1): Sequential(\n",
      "    (0): ReflectionPad2d((0, 0, 0, 0))\n",
      "    (1): Conv2d(15, 10, kernel_size=(1, 1), stride=(1, 1), bias=False)\n",
      "  )\n",
      "  (drelu1): ReLU()\n",
      "  (dbn1): BatchNorm2d(10, eps=1e-05, momentum=0.1, affine=True, track_running_stats=True)\n",
      "  (dups1): Upsample(scale_factor=2, mode=bilinear)\n",
      "  (dconv2): Sequential(\n",
      "    (0): ReflectionPad2d((0, 0, 0, 0))\n",
      "    (1): Conv2d(10, 1, kernel_size=(1, 1), stride=(1, 1), bias=False)\n",
      "  )\n",
      "  (sig): Sigmoid()\n",
      ")\n"
     ]
    }
   ],
   "source": [
    "if dataset == 'mnist':\n",
    "    num_channels = [25,15,10] \n",
    "elif dataset== 'celeba':    \n",
    "    num_channels = [120,25,15,10] \n",
    "else:\n",
    "    num_channels = [512,256,128]\n",
    "output_depth = img_np.shape[0] # number of output channels\n",
    "net = autoencodernet(num_output_channels=output_depth,num_channels_up=num_channels,need_sigmoid=True, Ameas=Ameas_var,\n",
    "                        decodetype=decodetype\n",
    "                        ).type(dtype)\n",
    "\n",
    "print(\"number of parameters: \", num_param(net))\n",
    "if decodetype == 'upsample':\n",
    "    print(net.decoder)\n",
    "elif decodetype == 'transposeconv':\n",
    "    print(net.convdecoder)\n",
    "net_in = copy.deepcopy(net)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Invert Image from Measurements with Deep Network Prior"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "shape of latent code B1:  [1, 25, 7, 7]\n",
      "initializing latent code B1...\n",
      "not optimizing over latent code Z1\n",
      "optimize decoder with adam 0.0001\n",
      "optimizing with gradient descent...\n",
      "\n",
      "LR is set to 0.0001\n",
      "\n",
      "\n",
      "Iteration 02999   Train loss 0.020388  \n",
      "LR is set to 7e-05\n",
      "\n",
      "\n",
      "Iteration 05999   Train loss 0.005642  \n",
      "LR is set to 4.9e-05\n",
      "\n",
      "\n",
      "Iteration 08999   Train loss 0.003089  \n",
      "LR is set to 3.4299999999999993e-05\n",
      "\n",
      "\n",
      "Iteration 09999   Train loss 0.002713  \n",
      "time elapsed: 37.747268199920654\n"
     ]
    }
   ],
   "source": [
    "#pick optimization procedure\n",
    "optim = 'gd'             #'pgd' (projected gradient descent), 'gd' (gradient descent)\n",
    "if optim == 'pgd':\n",
    "    optimizer2='SGD'      #outer loop optimizer - 'SGD' (or try 'adam')\n",
    "    numit = 1000          #number of outer iterations of LS\n",
    "    LR_LS = 10            #typically 5-10 ; required for outer loop of LS\n",
    "\n",
    "    OPTIMIZER='SGD'       #inner loop optimizer - SGD or adam\n",
    "    numit_inner = 10      #number of inner loop iterations for projection\n",
    "    LR = 0.5             #typically 0.02-0.05 for pgd/inner loop of projection, higher for more complex structures\n",
    "\n",
    "    lr_decay_epoch = 500  #decay learning rates of both inner and outer optimizers\n",
    "    \n",
    "elif optim == 'gd':\n",
    "    OPTIMIZER='adam'       #optimizer - SGD or adam \n",
    "    numit = 10000         #number of iterations for SGD\n",
    "    LR = 0.0001              #typically 0.02-0.5 for gd , higher for more complex structures\n",
    "\n",
    "    optimizer2 = None                                    \n",
    "    numit_inner = None\n",
    "    LR_LS = None\n",
    "\n",
    "    lr_decay_epoch = 3000\n",
    "    \n",
    "t0 = time.time()\n",
    "mse_t, ni, net, ni_mod, in_np_img = fit( \n",
    "                            net=net,\n",
    "                            num_channels=num_channels,\n",
    "                            num_iter=numit,\n",
    "                            numit_inner = numit_inner,\n",
    "                            LR=LR,\n",
    "                            LR_LS = LR_LS,\n",
    "                            OPTIMIZER = OPTIMIZER,                          \n",
    "                            optimizer2 = optimizer2,             \n",
    "                            lr_decay_epoch = lr_decay_epoch,             \n",
    "                            img_clean_var=img_var_meas,\n",
    "                            find_best=True,\n",
    "                            Ameas = Ameas_var,\n",
    "                            model = mode,\n",
    "                            code='uniform',\n",
    "                            decodetype=decodetype,\n",
    "                            optim=optim,\n",
    "                            out_channels=out_ch        \n",
    "                            )\n",
    "t1 = time.time()\n",
    "print('\\ntime elapsed:',t1-t0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Plot loss"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7efe514728d0>]"
      ]
     },
     "execution_count": 60,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAEOCAYAAACetPCkAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAH7pJREFUeJzt3Xl4VfW97/H3d2eHDCSEISEEEhIggKCiSFBoHai1OOJQBetQH4cDV1s9R2tra2/vc+552nM6eE/vc6xa63Ss1euEreKIlmrFmYCMIoKMQSDMhISM+3f/2Bu6jQGyIDtr77U/r+dZz87+7d9a65sfIZ+s2ZxziIiIdFbI7wJERCS1KDhERMQTBYeIiHii4BAREU8UHCIi4omCQ0REPFFwiIiIJwoOERHxRMEhIiKeKDhERMSTsN8FJEJhYaGrqKjwuwwRkZQxf/78bc65os70DWRwVFRUUF1d7XcZIiIpw8zWdbavdlWJiIgnCg4REfFEwSEiIp4oOERExBMFh4iIeKLgEBERTwJ5Om59cyvz1u74Upt10M86amzXs6M+HS/LOtGno2Udfn0dad8vZBabop9Z3PuQWVzbP95/aR4MC/GVeUIWrTB+GSKS3gIZHKu31jP1/vf9LiOQ/hE4sbCxDsImFB82RjhkZISMzIzoazgUIpwRbQ9nhGJt0a/DB742MkIhMmPzRvuH4pYTOtAvMyNEVjhEj3CIHhmx17ivs8IZsdcO+sTeZ4VDCkWRTgpkcAwp7MlDN5xy4L3DfaWP+2rTV3q5Djp1MNtXGju9vq/M11Gfw9cQ7eKIOIi46KtzDhf3PuIc7sBnX35/VPMQmyfSwTw42iLRz1oiEdoijtY2R2vs65Y2F22LRGhujdDQ3EZrJELrgXZ34H1rxMXm+fJyIh3+gxyZ/YGSnZlBbo/olLP/NTP8lbbcHmFy9vfNClOQk0nvnEz65Pagd89M8rPCCiMJpEAGR15WmFOHF/pdhnSDSMTR3BaJTq1xU1uEppYIzW1tNLVr/2qfyIE+Ta1tNLa00dAcnfY1t7GvpY0d9fvY19x6oK2hpY22w6RWOGT0zs2kd24P+uRmUpDTg749MynKz6J/fnbsNYui2JTbI5D/HSWA9JMqKS0UMrJDGWRnZnTrep2LBta+5v0h08rufS3srG9hZ0MzuxqirzsbWtjV0MzOhmZqdjawuKaZ7fXNHYZOXlb4QIiUFGRT2ieH0j65B14H9s4mK9y936dIRxQcIkfAzMgKZ5AVzqB3rrd5IxHHjoZmttY1UVvXxNbYVFvXeKDt4/W7eHnxJlrjAsYMivOjgTKsKI/hxXmMKM5nRHE+xb2ytFtMuo2CQ6SbhUJGYV4WhXlZjCo5eL+2iGPznkZqdjRQs3MfNTv3sWFnA+t3NPDX5Vt4unrDgb752WFGFOczvH80TI4bVMCxA3vRM0v/xaXrWUcHX1OVmU0BplRWVk5fuXKl3+WIJNT2vU18tmUvK2vrWLllL59tqWNl7V521DcD0S2UyqI8xpT2ZkxpAceXFjC6pFe379aT1GBm851zVZ3qG6Tg2K+qqsrptuqSrmrrGlm6cTeLa3azpGY3i2p2s21vExA9YH9iWW8mDuvHxGH9OGlwHwWJAAoOBYdIHOeiu7wW1+xmwfqdfLB6B0tqdhFxkBUOcWplId8aXcw3RxVTlJ/ld7niEy/BoR2gIgFnZpQU5FBSkMPZxw4AYE9jC/PW7GDuym288ckW5nxai9kSqsr7cMnYUs4fU0JBTqbPlUuy0haHSJpzzrF8Ux1vfLKFlxZ/wcravfQIh5g8upjvjB/M1yv76YytNKBdVQoOkSPinGPJxt08N7+GFxZ9wa6GFkYW53P9qRVcdOIgHQ8JMAWHgkPkqDW1tvHiok08/M4alm/aQ1F+FtNPG8JVp5TrNN8AUnAoOES6jHOO9z/fzn1vfc47q7bROzeT6acN5ZqJ5eRn6zhIUCg4FBwiCbFg/U5+N2clb67YSu/cTGacPpTrvjaEnB7ahZXqFBwKDpGEWrRhF/81ZyV/+7SWAb2yue1bw7n0pFLCGXo2XKryEhz6VxYRz04o680j147n6RkTGFCQzY+fW8K5/zWXv36ypcNHAUiwKDhE5IidMrQff/ne1/j9VSfRGnH802PVXP6HD1iwfqffpUkCKThE5KiYGeceX8Lrt53Ozy8+jtXb6vn2fe9x61Mfs2VPo9/lSQIoOESkS2RmhPjuhHL+/qNJ3HJmJa8s2cyZ/+ctHnx7NS1tEb/Lky6k4BCRLtUzK8ztk0fy+m2nc/KQvvz7K8u56J53WbG5zu/SpIsoOEQkISoKe/LIteO5/+pxbNnTyJR73uGhuauJdOWD4sUXCg4RSRgz45zjBjD7ttM5fXgRv3h5Odc9Oo89jS1+lyZHQcEhIglXmJfFg9eM4+cXH8e7q7Zxyb3vsmZbvd9lyRFScIhItzAzvjuhnD/dcAo76pu55L53mb9Op+2mIgWHiHSricP68fz3v07vnEyueugD5izf4ndJ4pGCQ0S6XXm/nsy86WuMKM5nxp/m88LCjX6XJB4oOETEF4V5WTw5fQLjK/pw29MLFR4pRMEhIr7pmRXmkWvHc/KQvgqPFKLgEBFf5faIhsf4imh4vL5ss98lyWEoOETEd7k9wvz3deM5flABtz69kGVf7Pa7JDkEBYeIJIXcHmEevKaKgpxMpv+xmto63SAxWSk4RCRp9O+VzYPXVLGzoYUZj82nsaXN75KkA0kfHGY21MweNrOZftciIol33KAC/u/lJ7Jwwy7umLlYD4ZKQgkNDjN7xMxqzWxpu/ZzzGyFma0ys58cahnOudXOuRsSWaeIJJdzjhvAj84eyaxFX3DP31b5XY60E07w8h8F7gEe299gZhnAvcC3gBpgnpnNAjKAX7ab/3rnXG2CaxSRJPS9ScNYVbuX/3zjM4b1z+O840v8LkliEhoczrm3zayiXfPJwCrn3GoAM3sKuMg590vggkTWIyKpw8z45bePZ932en7wzELK++Vy7MACv8sS/DnGMQjYEPe+JtbWITPrZ2b3A2PN7M5D9JthZtVmVr1169auq1ZEfJOdmcEfvhs90+rWpxbqYHmSSPqD48657c65G51zw2JbJQfr94Bzrso5V1VUVNSdJYpIAhXlZ/Gby05gZe1efvPaCr/LEfwJjo1AWdz70libiEiHzhhRxDUTy3nk3TW8t2qb3+WkPT+CYx4w3MyGmFkP4DvALB/qEJEUcue5oxha2JMfPruI3fv0BEE/Jfp03CeB94GRZlZjZjc451qBm4HZwHLgGefcskTWISKpL6dHBr+9/ES21DXxb7P0K8NPiT6r6oqDtL8CvNLV6zOzKcCUysrKrl60iCSBE8t68/1vVHL3nJWcNbpYp+j6JOkPjnvhnHvROTejoECn7IkE1S1nVjKmtICf/mUJtXt0Pys/BCo4RCT4MjNC/HbaiexrbuOO53RLEj8oOEQk5VT2z+POc4/hrRVb+X8frfe7nLSj4BCRlHTNxApOrSzkFy8tZ+22er/LSSuBCg4zm2JmD+zerYfAiARdKGTcNXUMmRnG7c8uIhLRLqvuEqjg0MFxkfRSUpDDv045lvnrdvLkPO2y6i6BCg4RST/fPmkQE4f249evfsrWuia/y0kLCg4RSWlmxi8uOY7Glgi/ePkTv8tJCwoOEUl5w4ryuGnSMF5Y+AVzV+ru2Imm4BCRQLhp0jCGFPbkfz2/lObWiN/lBFqggkNnVYmkr+zMDP51ymjWbm/gTx+s87ucQAtUcOisKpH0Nmlkf04bXsjdc1ayu0F30E2UQAWHiMhPzxvFnsYW7nlzpd+lBJaCQ0QCZVRJL6aOK+WP761j/fYGv8sJJAWHiATO7ZNHkhEyfjP7U79LCSQFh4gETnGvbGacPpSXFm9iwfqdfpcTOAoOEQmkGacPpSg/i/94ebluvd7FAhUcOh1XRPbrmRXm9m+NoHrdTl5butnvcgIlUMGh03FFJN7UqjJGFufzq9c+1UWBXShQwSEiEi8jZNx53jGs297AU7p7bpdRcIhIoJ0xoojxFX24783PaWpt87ucQFBwiEigmRm3njWCzXsaeWbeBr/LCQQFh4gE3teG9WNceR/u//tqWtp0rONoKThEJPDMjJu/UcnGXft4YeEXfpeT8hQcIpIWJo0sYlRJL+57a5WeT36UAhUcuo5DRA7GzLhp0jBWb61nzqe1fpeT0gIVHLqOQ0QO5bzjBjCodw4Pvr3a71JSWqCCQ0TkUMIZIa4/dQgfrd3Bx7qH1RFTcIhIWrl8fBn52WEemrvG71JSloJDRNJKXlaYq04p59Wlm/S8jiOk4BCRtHPd1ysImfHY+2v9LiUlKThEJO0U98rm7OMG8Ez1BvY16zYkXik4RCQtXTOhnD2Nrby4SBcEeqXgEJG0dPKQvowszuexD9bqQU8eKThEJC2ZGVdPLGfpxj18vGGX3+WklEAFh64cFxEvLhk7iLysMI+/v87vUlJKoIJDV46LiBd5WWG+fdIgXlq8ie17m/wuJ2UEKjhERLz67oRymtsiPFNd43cpKUPBISJpbXhxPhOG9uWJD9fprrmdpOAQkbR3xcmDqdm5jw/WbPe7lJSg4BCRtHf2sQPIzw7zrHZXdYqCQ0TSXnZmBlNOGMirSzexp7HF73KSnoJDRASYVlVGY0uElxdv8ruUpKfgEBEBTigtYHj/PJ6t3uB3KUlPwSEiQvRK8mlVZSxYv4tVtXV+l5PUFBwiIjEXjx1ERsh4dr4Okh+KgkNEJKYoP4tvjOzPnxdspLUt4nc5SUvBISISZ2pVKVvrmnj3c13TcTCBCg7d5FBEjtakkUX0yg7zwscb/S4laXUqOMzsX8ysl0U9bGYLzGxyoovzSjc5FJGjlRXO4LzjS5i9bLOeDngQnd3iuN45tweYDPQBvgv8KmFViYj46KITB1Hf3MZfl2/xu5Sk1NngsNjrecCfnHPL4tpERALllCF9GdArmxcWandVRzobHPPN7HWiwTHbzPIBnXIgIoEUChkXnjiQt1ZsZWd9s9/lJJ3OBscNwE+A8c65BiATuC5hVYmI+OzCEwbSGnHMXrbZ71KSTmeDYyKwwjm3y8yuBn4G6NQlEQmsYwf2orxfLi8v0b2r2utscPweaDCzE4Dbgc+BxxJWlYiIz8yMC8aU8N7n2/VY2XY6GxytzjkHXATc45y7F8hPXFkiIv67YMxA2iKOV5dqd1W8zgZHnZndSfQ03JfNLET0OIeISGAdMyCfoUU9eWnxF36XklQ6GxyXA01Er+fYDJQCdyWsKhGRJBDdXTWQD9fsoLau0e9ykkangiMWFk8ABWZ2AdDonNMxDhEJvAvGlOAcvLpEu6v26+wtR6YBHwFTgWnAh2Z2WSILExFJBiOK8xlRnKfdVXHCnez3P4lew1ELYGZFwF+BmYkqTEQkWVwwZiC/feMzNu9uZEBBtt/l+K6zxzhC+0MjZruHeUVEUtr5Y0oAdE1HTGd/+b9mZrPN7FozuxZ4GXglcWWJiCSPYUV5jCrppd1VMZ09OP4j4AFgTGx6wDn340QWJiKSTC48YSAfr9/Fmm31fpfiu07vbnLOPeec+0Fs+ksiixIRSTaXjB1EyODpeRv8LsV3hwwOM6szsz0dTHVmtqe7iuwsPQFQRBJlQEE2k0cP4MmP1tPQ3Op3Ob46ZHA45/Kdc706mPKdc726q8jO0hMARSSRpp8+hN37Wni2usbvUnylM6NERDppXHlfxg7uzcPvrKEt4vwuxzcKDhERD6afNpT1Oxp445P0faysgkNExIOzjx1AWd8cHpq72u9SfKPgEBHxICNkXPe1IVSv28nH63f6XY4vFBwiIh5NG19GfnaYh+au8bsUXyg4REQ8yssKc+Upg3l16SZqdjb4XU63U3CIiByBayZWAPDkR+v9LcQHCg4RkSMwqHcOZx7Tn6fnbaC5NeJ3Od1KwSEicoSumlDOtr3NzF6WXg95UnCIiByhM4YXUdY3h8c/WOd3Kd1KwSEicoRCIePKk8v5cM0OVm6p87ucbqPgEBE5CtOqSumREeKJD9PnILmCQ0TkKPTLy+Lc4wfw3PyatLlrroJDROQoXT2hnLqmVmYtTI8nBCo4RESOUlV5H0YW5/P4h+twLvh3zVVwiIgcJTPj6gmDWbpxD4tqgv8gOQWHiEgXuHjsIHJ7ZKTFqbkKDhGRLpCfnclFJw7ixUVfsKuh2e9yEkrBISLSRa6eMJim1gjPLdjodykJpeAQEekixw4sYOzg3jwR8IPkCg4RkS509SnlrN5az/ufb/e7lIRRcIiIdKHzx5TQOzeTxz8M7kFyBYeISBfKzsxg6rhSXl+2hdq6Rr/LSQgFh4hIF7t8fBmtEcfzHwfzILmCQ0Ski1X2z+ekwb15promkAfJkz44zOxiM3vQzJ42s8l+1yMi0hnTqspYVbuXjzfs8ruULpfQ4DCzR8ys1syWtms/x8xWmNkqM/vJoZbhnHveOTcduBG4PJH1ioh0lfPHlJCTmcGz1Rv8LqXLJXqL41HgnPgGM8sA7gXOBUYDV5jZaDM73sxeajf1j5v1Z7H5RESSXn52JuePKeHFRZsCd7v1hAaHc+5tYEe75pOBVc651c65ZuAp4CLn3BLn3AXtplqL+jXwqnNuQSLrFRHpStOqytjb1MqrS4L1THI/jnEMAuK33WpibQdzC3AWcJmZ3XiwTmY2w8yqzax669atXVOpiMhRGF/Rh4p+uTwTsN1VSX9w3Dl3t3NunHPuRufc/Yfo94Bzrso5V1VUVNSdJYqIdMjMmFpVxodrdrB2W73f5XQZP4JjI1AW97401iYiEjiXnlRKyGDm/Bq/S+kyfgTHPGC4mQ0xsx7Ad4BZPtQhIpJwAwqyOWNEETPn19AWCcY1HYk+HfdJ4H1gpJnVmNkNzrlW4GZgNrAceMY5tyyRdYiI+GlaVRmb9zQyd2Uwjr+GE7lw59wVB2l/BXilq9dnZlOAKZWVlV29aBGRI/bNUcX07dmDZ6trmDSy/+FnSHJJf3DcC+fci865GQUFBX6XIiJyQI9wiItPHMTrn2xmR33qPx0wUMEhIpKspo0vpaUtGDc+VHCIiHSDYwb0YkxpAc9Ub0j5Gx8qOEREusnUqjI+3VzH0o17/C7lqAQqOMxsipk9sHv3br9LERH5igtPGEhWOJTyV5IHKjh0cFxEkllBTibnHjeAFxZupLGlze9yjliggkNEJNlNqypjT2Mrs5el7o0PFRwiIt1owtB+lPbJ4dnq1L0FiYJDRKQbhULG1HFlvPv5NjbsaPC7nCOi4BAR6WaXjos+SeK5Bam51RGo4NBZVSKSCkr75HJqZSHPVtcQScEbHwYqOHRWlYikiqlVZWzctY/3V2/3uxTPAhUcIiKpYvLoYgpyMlPymg4Fh4iID7IzM7jwhIG8tnQzexpb/C7HEwWHiIhPLhtXSlNrhJcWbfK7FE8UHCIiPhlTWsDw/nnMnJ9au6sUHCIiPjEzplaVsmD9LlbV7vW7nE4LVHDodFwRSTUXjx1ERshS6pqOQAWHTscVkVTTPz+bSSOK+POCGtpS5JqOQAWHiEgqumxcKVv2NDFn+Ra/S+kUBYeIiM/OGl1Meb9cfvvGZylxJbmCQ0TEZ5kZIW6fPJJPN9cxa9EXfpdzWAoOEZEkcMHxJYwu6cV/vrGC5taI3+UckoJDRCQJhELGHeeMZMOOfTz50Xq/yzkkBYeISJI4Y0QRE4b25Xd/W0l9U6vf5RxUoIJD13GISCozM+445xi27W3mkXfW+F3OQQUqOHQdh4ikupMG9+HsY4v5w9ur2VHf7Hc5HQpUcIiIBMEPJ4+kobmV+95c5XcpHVJwiIgkmeHF+VwytpQ/fbCOLXsa/S7nKxQcIiJJ6F++OZy2iOPeJNzqUHCIiCShwf1ymTa+jCc/Wk/Nzga/y/kSBYeISJK65cxKzIy756z0u5QvUXCIiCSpkoIcrjx5MM8t2MiGHcmz1aHgEBFJYjdNGkZGyJLqWIeCQ0QkiRX3yuaK8WXMnF+TNMc6AhUcunJcRILoxknDCJlx31uf+10KELDg0JXjIhJEJQU5TK0q5dnqDWzctc/vcoIVHCIiQfW9b1QCcH8SbHUoOEREUsCg3jlcNq6Up+dtYNNuf7c6FBwiIinie5MqaXOOh+f6e+dcBYeISIoo65vLBWNKePKj9exuaPGtDgWHiEgKufGMYdQ3t/H4h+t8q0HBISKSQkaV9OK04YU8+t5amlrbfKlBwSEikmKmnzaUrXVNzFr4hS/rV3CIiKSY04YXcsyAfB5+Zw3OuW5fv4JDRCTFmBn/dNpQPt1cx9yV27p9/QoOEZEUdOEJAynulcVds1fQ0hbp1nUrOEREUlCPcIh/nXIsSzbu5uF3uve6jkAFh25yKCLp5LzjSzhrVDG/m7OS2m58NnmggkM3ORSRdPOz80fR0ub4zewV3bbOQAWHiEi6qSjsyfWnDmHm/BoWbtjVLetUcIiIpLibz6ykKD+L/z1rGZFI4k/PVXCIiKS4vKwwd5w9kr1NrWzb25Tw9YUTvgYREUm4S08q5eKxg8jMSPz2gIJDRCQAQiEjhHXPurplLSIiEhgKDhER8UTBISIinig4RETEEwWHiIh4ouAQERFPFBwiIuKJ+fH0qEQzs63A/ie5FwDtb5cb39b+80IgUU9G6aiWrprnUP0O9tnhxuZgbfHvNV4aL42Xt37JOl7lzrmiTvV0zgV6Ah44VFv7z4Hq7qylq+Y5VL+DfXa4sTnEGMWPn8ZL46XxSoPxip/SYVfVi4dp6+jzRDmSdXV2nkP1O9hnhxubg7V115hpvLzReHmj8TpCgdxVdTTMrNo5V+V3HalC4+WNxssbjZc33TVe6bDF4dUDfheQYjRe3mi8vNF4edMt46UtDhER8URbHCIi4omCQ0REPFFwiIiIJwqOQzCznmb2RzN70Myu8rueVGBmQ83sYTOb6XctqcDMLo79fD1tZpP9rifZmdkoM7vfzGaa2U1+15MKYr/Hqs3sgq5aZtoFh5k9Yma1Zra0Xfs5ZrbCzFaZ2U9izd8GZjrnpgMXdnuxScLLmDnnVjvnbvCn0uTgcbyej/183Qhc7ke9fvM4XsudczcC04Cv+1Gv3zz+DgP4MfBMV9aQdsEBPAqcE99gZhnAvcC5wGjgCjMbDZQCG2Ld2rqxxmTzKJ0fMzmy8fpZ7PN09CgexsvMLgReBl7p3jKTxqN0crzM7FvAJ0BtVxaQdsHhnHsb2NGu+WRgVeyv5WbgKeAioIZoeEAajtV+Hscs7XkZL4v6NfCqc25Bd9eaDLz+fDnnZjnnzgXScvexx/GaBEwArgSmm1mX/B4Ld8VCAmAQ/9iygGhgnALcDdxjZufjw2X9Sa7DMTOzfsC/A2PN7E7n3C99qS75HOxn7BbgLKDAzCqdc/f7UVwSOtjP1ySiu5CzSN8tjo50OF7OuZsBzOxaYJtzLtIVK1NwHIJzrh64zu86UolzbjvR/fXSCc65u4n+gSKd4Jx7C3jL5zJSjnPu0a5cXtrufmlnI1AW97401iYHpzHzRuPljcbLm24dLwVH1DxguJkNMbMewHeAWT7XlOw0Zt5ovLzReHnTreOVdsFhZk8C7wMjzazGzG5wzrUCNwOzgeXAM865ZX7WmUw0Zt5ovLzReHmTDOOlmxyKiIgnabfFISIiR0fBISIinig4RETEEwWHiIh4ouAQERFPFBwiIuKJgkPShpndama5ce9fMbPeHua/sN3tqruqrvdirxVmdmUXL/unHa1L5GjoOg5JG2a2Fqhyzm3zsYZw7GKtjj6bBPzQOdfpB+4canmxz/c65/K8VypycNrikJRlZj8ws6Wx6dZYW4WZfWpmT5jZ8tiT4nLN7J+BgcCbZvZmrO9aMyuMm+dRM/ssNu9ZZvauma00s5Nj/a81s3tiXy+Mm/aZ2RkWfdLaI2b2kZl9bGYXxc03y8z+Bszp4PvYG/vyV8BpsWXeZmYZZnaXmc0zs8Vm9j9i/SeZ2Vwzm0X0WQuY2fNmNt/MlpnZjFjbr4Cc2PKeiF9X7Hbud8XGbomZXR637Ldi47Z/HC0B/3ySypxzmjSl3ASMA5YAPYE8YBkwFqgAHPD1WL9HiP4VD7AWKIxbxlqgMDZPK3A80T+m5sfmM6LPNHg+1v9a4J52dUwB5gKZwH8AV8faewOfxeq7luhtrvse5HvZG3udBLwU1z4D+Fns6yygGhgS61cPDInr2zf2mgMsBfrFL7uDdV0KvAFkAMXAeqAktuzdRG+SFyJ6a4tT/f731pRck7Y4JFWdCvzFOVfvnNsL/Bk4LfbZBufcu7GvH4/1PZw1zrklLvq8gmXAHOecIxpOFR3NYGbDgbuAac65FmAy8BMzW0j01t/ZwOBY9zecc+0fvnM4k4FrYsv7EOgHDI999pFzbk1c3382s0XAB0TvkjqcQzsVeNI51+ac2wL8HRgft+ya2Fgs5CDfv6QvPY9Dgqj9gbvOHMhrivs6Evc+Qgf/T8wsj+hznKc75zbtbwYudc6taNf3FKJbCF4ZcItzbna75U2KX17s/VnAROdcg5m9RTS0jlT8WLSh3xPSjrY4JFXNBS6OHb/oCVwSawMYbGYTY19fCbwT+7oOyO+i9T8C/Ldzbm5c22zglv3HBMxsrMdltq9vNnCTmWXGljci9r22VwDsjIXGMUQfFbpfy/7525kLXB47jlIEnA585LFeSVMKDklJLvp87keJ/rL7EHjIOfdx7OMVwPfNbDnQB/h9rP0B4LX9B8ePlJmVA5cB18cdIK8Cfk70WMdiM1sWe+/FYqDNzBaZ2W3AQ0QPfi8ws6XAH+j4r//XgHDs+/0V0d1V+z0Qq+eJdvP8Jba+RcDfgDucc5s91itpSqfjSqCYWQXRA8zH+VyKSGBpi0NERDzRFoeIiHiiLQ4REfFEwSEiIp4oOERExBMFh4iIeKLgEBERTxQcIiLiyf8HqOaDZLUzRiAAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.xlabel('optimizer iteration')\n",
    "plt.ylabel('loss')\n",
    "plt.loglog(mse_t)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Compute initialization error"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.5498978\n",
      "0.75879157\n",
      "0.7395467\n"
     ]
    }
   ],
   "source": [
    "lvls = len(num_channels)\n",
    "if decodetype == 'upsample':\n",
    "    nettype = net.decoder\n",
    "    netintype = net_in.decoder\n",
    "elif decodetype == 'transposeconv':\n",
    "    nettype = net.convdecoder\n",
    "    netintype = net_in.convdecoder\n",
    "ComputeInitErr(nettype,netintype,lvls)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Display reconstructed image and statistics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {},
   "outputs": [],
   "source": [
    "out_img_np = net( ni.type(dtype) ).data.cpu().numpy()[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Image inversion with Deep-Decoder, SNR: 25.62715799918684\n",
      "MSE: 0.022399362\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPsAAAD4CAYAAAAq5pAIAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAACsJJREFUeJzt3T1sjf8fxvG7HnpaDxUUTSRVQVMRicRAhBiIgRSLxGqyGLpUiAmDwcDQwcNCDUJMQqykAymhUotEBImH0Iq0jT7zn//ivj7/f7+5z6le79d65XvO7Zxevzv5fc73e1f9/v07AzDzzar0BQAoD8oOmKDsgAnKDpig7ICJOeV8s+rq6n/2f/1XVVXlZrNm6f9mqrVZlmXRRGRycnLK+a9fv+TaVNG/TeWzZ8+Wa6M8eu8U0XdSdJ6ydnx8/K8fDHd2wARlB0xQdsAEZQdMUHbABGUHTFB2wERZ5+wzVerOwSJ3HqbO+It+fyW6tiLn7DMRd3bABGUHTFB2wARlB0xQdsAEZQdMUHbABHP2/1GR+4+L3PscSZ1VV3LWzRz+/8OdHTBB2QETlB0wQdkBE5QdMEHZARM2o7fUMYwa8xQ9Oouuvbq6OjebM0d/xdFxzdEx2dFR1RMTE1PKsiw+QjuiPreUI7ArbarXxp0dMEHZAROUHTBB2QETlB0wQdkBE5QdMFHWOXuR8+ai56Lq9aNZc5RH115bWyvzjRs35mbt7e1ybUtLi8xHR0dl3tvbK/POzs7crLu7W64dHByUecoW1+gzT30Md0StL2pLM3d2wARlB0xQdsAEZQdMUHbABGUHTFB2wMSM2c/+Lx8rHO0p37Bhg8xPnTqVm23fvl2urampkfnY2JjMV6xYIfOmpqbc7Pz583LtgwcPZB7th085gyCSOqevhOl3RQAKQdkBE5QdMEHZAROUHTBB2QETlB0wUVXk44D/VF1dXb43+0PRc9WUtcuXL5f57du3Za72s0eiWfX3799lHs3p6+vrc7MfP37Itfv27ZP5y5cvZR7921JEc/TotxMp+9mjfGxs7K8vzp0dMEHZAROUHTBB2QETlB0wQdkBEzZbXKPjnFO2yEZjmLq6Opm3trbKvLm5Webq33bnzh259tatWzL/9OmTzKMttOfOncvNFixYINeeOXNG5qdPn5Z5T09PbjY+Pi7Xpo6/itxyzSObAUiUHTBB2QETlB0wQdkBE5QdMEHZARNl3eI6d+7cpDdT88XJyUm5Npqzp4jm7Pv375f52bNnZd7Q0CDzmzdv5mZtbW1ybernEs18N2/enJtduXJFrl23bp3MBwYGZH7kyJHcrKurS66NHlUd9SZ1C2yK0dFRtrgCzig7YIKyAyYoO2CCsgMmKDtggrIDJsq6n306Pza5SIcOHZJ5Y2OjzKMjlzs7O3OzIn9fkGXxvPnFixe5WbRf/fr16zJfunSpzNvb23Ozr1+/yrWvXr2SefS7joj63IrqCXd2wARlB0xQdsAEZQdMUHbABGUHTFB2wERZ97PX1tYmvZm61ujxvKnnxiulUknmatacZVnW1NQk8+7ubpkfPnw4N/vy5YtcW0nR59bR0SHzgwcPylx9pw8fPpRrjx49KvNoL30lf1MyPj7OfnbAGWUHTFB2wARlB0xQdsAEZQdMUHbAxIzZzx6d013k7wnmzZsn8/r6eplH1/7+/XuZR88an66i67548aLM16xZI/Nt27blZrt27ZJrm5ubZf7y5UuZR/vd1d9jUX+r3NkBE5QdMEHZAROUHTBB2QETlB0wUdbRWzRSiEZzakQVvXa0xTVlLLho0SKZV1dXT/m1syzLvn//LvOampopv3c0IooeLZz6uStv376V+dWrV2W+Y8eO3Cz6zvbu3SvzN2/eyHxoaEjmakt29JlOdTTHnR0wQdkBE5QdMEHZAROUHTBB2QETlB0wUdY5e+qsW818o7Wp22vVjH/JkiVJ7x3Nunt6emQ+PDycm82dO1eujeboUZ5yhHc0Lx4dHZX5o0ePZK6uLfr9wapVq2Q+f/58mUdzdra4AigMZQdMUHbABGUHTFB2wARlB0xQdsDEPzVnr+RjcNWcPTpKOvp3qzl5lmXZ48ePZT44OCjzFKmPwo5+Q5Cir69P5q9fv87NVq9eLdf29/fLvJJHl7OfHYBE2QETlB0wQdkBE5QdMEHZAROUHTAxrc6NL3q9Es3w1f7naO9ztHc62rcd5Sl7xiNFfuapUvbiR/+uaD969LjpSv6t5uHODpig7IAJyg6YoOyACcoOmKDsgAnKDpgo65w92gMcKfKs7eja1DnhW7dulWtLpZLM58zRX8PKlStl/vHjx9ws2m9eyTl6NC+OPpctW7bIvLm5OTeLvu/oc/v586fMUz7X1J7kvm4hrwpg2qHsgAnKDpig7IAJyg6YoOyAibKO3qJRSzSuSBm9pebqSOVoTJM6Som20D59+jTp9YukvvPocdLr16+X+cmTJ6f83iMjI3JtdHx3tD6irq2oI9W5swMmKDtggrIDJig7YIKyAyYoO2CCsgMmyjpnj6TO4VNEjxYeGBjIzaKZ7LFjx2QezeGjLa7qyOTokcuVVFdXJ/Nr167JfO3atTL/9OlTbnb9+nW59smTJzIfGxuTecqsvKhHk3NnB0xQdsAEZQdMUHbABGUHTFB2wARlB0xMq6OkU+boqfvVI2pe3dvbO+W1WRYfmbxx40aZL1q0KDfr6+uTa6O9+JHo2hsbG3OzEydOyLXRHD26djWnv3z5slw7PDws8+jvqag96Sm4swMmKDtggrIDJig7YIKyAyYoO2CCsgMmyjpnj2ay0dxU7Tkves6u1n/+/Fmu7e/vl/nSpUtlvmnTJpnv3r07N7t3755cG+3LXrx4scyPHz8uc7WXX+3Dz7IsGxoakvndu3dlfuXKldws+k6iv5foNyMpc/aizm3gzg6YoOyACcoOmKDsgAnKDpig7ICJqiKPZ/7TwoUL5ZtFoze1VTQ6Cjp67ZQti6VSSa7t6OiQeWtrq8yj8di3b99ys2fPnsm1P378kLnaopplWbZz506ZNzQ05GbRZ37//n2Zt7W1yfzDhw+5WfT3Eo3OUkdvSmonx8bG/vrm3NkBE5QdMEHZAROUHTBB2QETlB0wQdkBE2Wds9fW1so3S9mmGs3RU+fsSjRzbWlpkXm0TXTPnj0yr6+vz82ibcWR0dFRmX/58kXmz58/z80uXbok13Z1dck8urYURc/ZVZ56vDdzdsAcZQdMUHbABGUHTFB2wARlB0xQdsBEWY+STp0fpszCUx+Rq947uq53797J/MKFCzIfGRmR+YEDB3KzZcuWybXRtUdHLt+4cUPmapb+9etXuTZ61HWKoh+pXOQjmafaA+7sgAnKDpig7IAJyg6YoOyACcoOmKDsgImy7mcvlUpJ+9lTpD7Sucj97qm5evRxtDYSna+ecl5/9Jmm/i4jZdadup89ytW/PfpMo89tYmKC/eyAM8oOmKDsgAnKDpig7IAJyg6YoOyAibLuZy/nTP9PKXPPKE+d0Uf7tov8DUCqlH3hqbPsIhW5H71SuLMDJig7YIKyAyYoO2CCsgMmKDtgoqyjtyKlHg0c5SnbLVNHY5UcrUWKHL2lvjf+G3d2wARlB0xQdsAEZQdMUHbABGUHTFB2wMQ/NWcv8mjgotcr03mOXqSi5+Qprz8TvxPu7IAJyg6YoOyACcoOmKDsgAnKDpig7ICJsj6yGUDlcGcHTFB2wARlB0xQdsAEZQdMUHbABGUHTFB2wARlB0xQdsAEZQdMUHbABGUHTFB2wARlB0xQdsAEZQdMUHbABGUHTFB2wARlB0xQdsAEZQdM/AcXDXZWCMLhAwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "maxv = np.max(img_np) \n",
    "print(\"Image inversion with Deep-Decoder, SNR: \" + str(psnr(img_np_orig,out_img_np,maxv)))  \n",
    "reconstruction_err = mse(img_np_orig,out_img_np,maxv)\n",
    "print('MSE:',reconstruction_err)\n",
    "if dataset == 'celeba':\n",
    "    plt.imshow(out_img_np.transpose(1,2,0))\n",
    "else:\n",
    "    plt.imshow(out_img_np[0,:,:])\n",
    "plt.axis('off')\n",
    "plt.show()\n",
    "save_path= optim+'_'+img_name+str(int(10*f))+'.png'\n",
    "savefig=False\n",
    "if savefig:\n",
    "    plt.savefig(save_path,bbox_inches='tight') "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Display initialization and statistics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Image at random initialization of Deep-Decoder, SNR: 6.450037630631322\n",
      "MSE: 0.89944655\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPsAAAD4CAYAAAAq5pAIAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAC5NJREFUeJzt3UluFFsehfGgx6Y3FrItISQQW2DADtgRK2E7rIE5AumJHmRjGtNDjWpGnM/KW0QVdb7fkL/CGRmZ56X0TtwbR379+jVJ+v939L99ApKWYdilEoZdKmHYpRKGXSpxfMkXu3fvXvxf/0eP5v/2HDt2bHb25cuXeOze3l6cP3v2LM4fPnw4O3v//n08dn19Pc7Pnj0b5+fPn4/ztbW12dmJEyfisSdPnozz06dPx/mlS5fifHt7e3a2sbERjz137lyc7+7uxvmDBw9mZ48ePYrHvnnzJs6pxdrc3Izz9N7pmtJncvfu3SO/+3d/2aUShl0qYdilEoZdKmHYpRKGXSph2KUSi/bsP3/+jPMjR35bDx5qTh39yN+eptzxp9l/Al2379+//7HXPn48f0Xo3NLx1PGP3HcxTdN06tSp2Rnd+0Bzet90f0N6b9Th02vPvuZKR0n66xh2qYRhl0oYdqmEYZdKGHaphGGXSizas1N/SPPUTVLnSn0x9aJpPnLsYXz79u2PHfvjx4+V//Y0TdPFixfjPK21py6bPlO6rmk9PK2lJ3RvA723dO7Uo6/6ffCXXSph2KUShl0qYdilEoZdKmHYpRKLVm/7+/txTlXL58+fZ2dUIR0cHMQ5HZ+WS9KxVBHR+6a/n6qYT58+xWOp7qRlqHTuN27cmJ1RHZpqu8PY2tqandHyWarO0ndxFH0mVm+SIsMulTDsUgnDLpUw7FIJwy6VMOxSiUV79qdPn8Y5bec8YqSrnqb8mNzRLZFHlzS+fft2dkbX/N27d3FOPnz4EOfXr1+fnV29ejUeS4+qpuW1qcenx2Rfvnw5zum60WO8030f1OGvunW4v+xSCcMulTDsUgnDLpUw7FIJwy6VMOxSiUV7duoPR7aapo6eum6anzlzZuVj6dyoR//48WOcp+tCffDr16/jnO5PuHDhQpzv7u7OzmitPaH7G9IeBGl2mL9Na+3pO5H+Pt274Hp2SZFhl0oYdqmEYZdKGHaphGGXShh2qcSiPTutIaaendZ9J7RHOe1/nuajHf7Xr1/jnK5L6oxpz3o6N+rZaW11ukdgdE34ly9f4jx1/HR/QTp2mvgzoUdCp7X49IwD+r7M8ZddKmHYpRKGXSph2KUShl0qYdilEotWb+nxvdPE1Vqa07FUlVDFlGoeWsKatqGeJj53Oj5VMS9fvozH0nJJqnlG6lRaykn1F9WpqdqjZcO0/JZem+ZpyTQtr6Xv6hx/2aUShl0qYdilEoZdKmHYpRKGXSph2KUSi/bst2/fjnPqm9NySlruSL3q3t5enD9//nyl85omfrQwbWtMfz9ta0x9MfXk1MNvbm7GebpHgJZyvnjxIs5p+W66R4Duu6C/TT06fWZpTj07Lcee4y+7VMKwSyUMu1TCsEslDLtUwrBLJQy7VGLRnn1nZ+eP/W3qk2ltNHX8qYenLnp9fT3Oqese2ap6f38/HksdPz1mmx5dnK4rnRutpadzT1346P4H1KOPPHaZril1/HP8ZZdKGHaphGGXShh2qYRhl0oYdqmEYZdKLNqzU7dJfXLqVakXHX1kczqe9o2n16a107S++cKFC7MzureB3vfbt2/jnO4xSI9dpv3Pab98un8h7c1O3zXaH4HOne77SK9P50bft9m/u9JRkv46hl0qYdilEoZdKmHYpRKGXSph2KUSi/bsjx8/jnPqfFNvSr0nrZ2mXjWdG/WihO4RoF419fSja+mpR6f9+NN6ePrbtK88febp+0Jr4ema07nRcwjS8bRWnu5XmeMvu1TCsEslDLtUwrBLJQy7VMKwSyUWrd7u378f57QUNG2xS1UKLSMdQedNFRFtmUzVXqqw6G/TudFrU12aqjfapppee2QJLD1qmurQV69exfk///wT50+ePJmdUZ1Jn+kcf9mlEoZdKmHYpRKGXSph2KUShl0qYdilEov27E+fPo1z6mzTlsr0mFta6kmdbXrt0Q6flizSPPWu1KPTUk7axpqWqSbUZf9J9L6p46dzp648LblO229PEy/HnuMvu1TCsEslDLtUwrBLJQy7VMKwSyUMu1Ri0Z49be17GKn7pD6Zek/qLlPfTF00rbWn9cnUZacenvpgureBzp22PU73IFDXPXpvBH3mCZ0bPcp61S58mvi+DR/ZLCky7FIJwy6VMOxSCcMulTDsUgnDLpVYtGcfXVud5rTmm/42nVvqTem1R/depx4+vTfqbOmRzXRd6NxTz09/m3p0ugcg7UtPn9novRPXr1+P87Rv/adPn+KxdG/DHH/ZpRKGXSph2KUShl0qYdilEoZdKmHYpRKL9uxpr+xp4uecp951dO3zyDPQCb0vmo/09PS+R+8/oD0Kzp8/PzujvpjW4tPxqz7HfJr4fW1vb8f5rVu34jzd/7C3txePpefaz/GXXSph2KUShl0qYdilEoZdKmHYpRKLVm9UGdCywlRX0JJFWopJ9VaqcUYfuUxzqpjSudP7Gn00MdWGqdobrTvp+HRuVNVubGzE+dWrV+P85s2bcZ5en7aptnqTFBl2qYRhl0oYdqmEYZdKGHaphGGXSizas9NyypE+mo6lvpjmaaknLQOlrYFHH6tM7z2hz4QePUxdeLo/gY6l+wvovoy0vHZraysee+3atTinJa60RDad+8WLF+Oxqy639pddKmHYpRKGXSph2KUShl0qYdilEoZdKrFoz05riKlPTnPqmlddA/xv1JUn1KNTlz2CrindI0CdLr23tAcBrYWnx03TY5PX1tZWPpbOjd73x48f4zx9H+n7QJ/Z3GO4/WWXShh2qYRhl0oYdqmEYZdKGHaphGGXSizas1+5ciXOR/ZPp7XP7969i/ORddmja8IJPXo43QMw0vdOE38mtG57Z2dndkZrwjc3N4deO30uBwcH8dgXL17E+YcPH+Kc7ilJ1/X9+/fxWPo+3blz57f/7i+7VMKwSyUMu1TCsEslDLtUwrBLJRat3ka3yE11ycix08RLWEeqN1pOSfUW1TypqqHH/9J1o6WeVI+l90bXJW0FPU1cb6XrQvUWXXM6dzq3kdemKnaOv+xSCcMulTDsUgnDLpUw7FIJwy6VMOxSiUV7dlqSOPp44IS6Seo2R7aiPn36dJzTtsQjnTDdP0D3CNBndu7cuZXndCz17EeP5t+q9H2h+w/o86Ytuum9pUc2jz5+fI6/7FIJwy6VMOxSCcMulTDsUgnDLpUw7FKJRXt26nSpN01z+tv0mFvqNmme0P0BNB/ZDpre9+i6bOqT5x4ffJjXJvSZpO8EfV9GH7NN1z19l0e+a/E1/8hflfQ/x7BLJQy7VMKwSyUMu1TCsEslDLtUYtGendaU03ykf6QOn/ZHX1tbm52NdvjUs9PjqNPr07rr9L6maZo2NjaG5iP7p9NjtukzSz0+nTfdP0D3Puzu7sZ52qOA1tLT92GOv+xSCcMulTDsUgnDLpUw7FIJwy6VMOxSiUV7dtr/nNYQpzn1ydSznzhxIs5Tp0v3B1BvSsdTrzpyXWhPe9o3nnr0dN3o/gTa85569vSZ0vum+w9G753Y399faTZNPp9dEjDsUgnDLpUw7FIJwy6VMOxSiUWrt729vTgfqYlGaphp4qokbT1MNQxVTFStjVRvVDmmRwdPE2/3TNc9XZvRLbbptdN1ofdNc7qu9JmnOpaWzx4cHMT5HH/ZpRKGXSph2KUShl0qYdilEoZdKmHYpRKL9uy0ZJE63dSzjzzu+TDz1LPTsdTDU48+soU23bswOqf3nrpu6qJpyTOheyuSke/DNI3dWzG65focf9mlEoZdKmHYpRKGXSph2KUShl0qYdilEkdGu0xJfwd/2aUShl0qYdilEoZdKmHYpRKGXSph2KUShl0qYdilEoZdKmHYpRKGXSph2KUShl0qYdilEoZdKmHYpRKGXSph2KUShl0qYdilEoZdKmHYpRL/Aq/BZ0A9/9sJAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Initialization error: 1.3114233\n"
     ]
    }
   ],
   "source": [
    "in_img = in_np_img.data.cpu().numpy()\n",
    "maxv = np.max(img_np) \n",
    "print(\"Image at random initialization of Deep-Decoder, SNR: \" + str(psnr(img_np_orig,in_img,maxv)))  \n",
    "reconstruction_err = mse(img_np_orig,in_img,maxv)\n",
    "print('MSE:',reconstruction_err)\n",
    "if dataset == 'celeba':\n",
    "    plt.imshow(in_img[0,:,:,:].transpose(1,2,0))\n",
    "else:\n",
    "    plt.imshow(in_img[0,0,:,:])\n",
    "plt.axis('off')\n",
    "plt.show()\n",
    "init_err = np.linalg.norm(out_img_np-in_img)/np.linalg.norm(out_img_np)\n",
    "print('Initialization error:', init_err)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Image inversion with sparsity priors"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Pick sparsifying basis - Compute wavelet or cosine transforms operators"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "metadata": {},
   "outputs": [],
   "source": [
    "basis = 'DCT'              # spatial, DWT, DCT\n",
    "if basis == 'DWT': #wavelet transform\n",
    "    Winv = construct_Wminv(d=img_np.shape[1],wave_name='db1') #imported from wavelet_DCT_basis.py\n",
    "    d_wav_image = int(np.sqrt(Winv.shape[1]))\n",
    "    print('Size of image:', img_np.size, ', size of wavelet transform matrix:', Winv.shape)\n",
    "elif basis == 'DCT': #DCT transform    \n",
    "    Dinv = construct_IDCT2mat(d=d)                            #imported from wavelet_DCT_basis.py"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Modify measurement operator to incorporate transform "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 78,
   "metadata": {},
   "outputs": [],
   "source": [
    "dsize = img_np.shape[1]*img_np.shape[2]\n",
    "if not (mode==1):\n",
    "    A = Ameas\n",
    "    if basis=='DWT': #wavelet\n",
    "        dwsize = Winv.shape[1]\n",
    "        Aeff = np.zeros((Ameas.shape[0],dwsize*out_ch))\n",
    "        for i in range(img_np.shape[0]):   \n",
    "            Anew = np.dot(A[:,i*dsize:(i+1)*dsize],Winv)\n",
    "            Aeff[:,i*dwsize:(i+1)*dwsize] = Anew\n",
    "    elif basis=='DCT': #DCT\n",
    "        Aeff = np.zeros((Ameas.shape[0],dsize*out_ch))\n",
    "        for i in range(img_np.shape[0]):   \n",
    "            Aeff[:,i*dsize:(i+1)*dsize] = np.dot(A[:,i*dsize:(i+1)*dsize],Dinv)    \n",
    "    else:    \n",
    "        Aeff = A\n",
    "    y = img_var_meas.detach().cpu().numpy()\n",
    "elif mode==1:\n",
    "    Aeff = np.identity(np.ravel(img_var.cpu()).shape[0])\n",
    "    y = np.ravel(img_var.cpu())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Store measurements and operator for MATLAB codes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 80,
   "metadata": {},
   "outputs": [],
   "source": [
    "#used for comparisons - TVAL3, SPARTA (supplementary code not provided in current version, will be updated)\n",
    "if not (mode==1):\n",
    "    save_mat = True\n",
    "    if save_mat:\n",
    "        sio.savemat('A.mat', {'A':Aeff})\n",
    "        sio.savemat('y.mat', {'y':y})\n",
    "        sio.savemat('xtrue.mat',{'x_':np.ravel(img_np)})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 81,
   "metadata": {},
   "outputs": [],
   "source": [
    "if not (mode==1):\n",
    "    img_vec = sio.loadmat('xtrue.mat')['x_']\n",
    "    #img_vec = -img_vec"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 84,
   "metadata": {},
   "outputs": [],
   "source": [
    "if not (mode==1):\n",
    "    if basis=='DWT': #wavelet\n",
    "        img_vec = np.reshape(img_vec,[img_np.shape[0],dwsize])\n",
    "        img_rec = np.zeros(img_np.shape)\n",
    "        for i in range(img_np.shape[0]):    \n",
    "            img_vec2 = np.dot(Winv,img_vec[i,:].T)\n",
    "            img_rec[i,:,:] = np.reshape(img_vec2,[img_np.shape[1],img_np.shape[2]])    \n",
    "    elif basis=='DCT': #DCT\n",
    "        img_vec = np.reshape(img_vec,[img_np.shape[0],dsize])\n",
    "        img_rec = np.zeros(img_np.shape)\n",
    "        for i in range(img_np.shape[0]):    \n",
    "            img_vec2 = np.dot(Dinv,img_vec[i,:].T)\n",
    "            img_rec[i,:,:] = np.reshape(img_vec2,[img_np.shape[1],img_np.shape[2]])                \n",
    "    else:\n",
    "        img_rec = np.reshape(img_vec,img_np.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 107,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(784,)"
      ]
     },
     "execution_count": 107,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "img_vec.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 103,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(784, 784)\n"
     ]
    }
   ],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Compressed sensing and denoising with Lasso"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 108,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Time taken for Lasso: 0.3528320789337158\n"
     ]
    }
   ],
   "source": [
    "if mode==1 or mode==2:\n",
    "    t0 = time.time()\n",
    "    reg = linear_model.Lasso(alpha=1e-6) #(max_iter=10000, fit_intercept=False, tol=0.00001)\n",
    "    reg.fit(Aeff, y)  \n",
    "    img_vec = reg.coef_\n",
    "    if mode==2:\n",
    "        if basis=='DWT':\n",
    "            img_vec = np.reshape(img_vec,[img_np.shape[0],dwsize])\n",
    "            img_rec = np.zeros(img_np.shape)\n",
    "            for i in range(img_np.shape[0]):    \n",
    "                img_vec2 = np.dot(Winv,img_vec[i,:].T)\n",
    "                img_rec[i,:,:] = np.reshape(img_vec2,[img_np.shape[1],img_np.shape[2]])    \n",
    "        elif basis=='DCT':\n",
    "            img_vec = np.reshape(img_vec,[img_np.shape[0],dsize])\n",
    "            img_rec = np.zeros(img_np.shape)\n",
    "            for i in range(img_np.shape[0]):    \n",
    "                img_vec2 = np.dot(Dinv,img_vec[i,:].T)\n",
    "                img_rec[i,:,:] = np.reshape(img_vec2,[img_np.shape[1],img_np.shape[2]])                \n",
    "        else:\n",
    "            img_rec = np.reshape(img_vec,img_np.shape)\n",
    "    else:\n",
    "        img_rec = np.reshape(img_vec,img_np.shape)\n",
    "    t1 = time.time()\n",
    "    print('Time taken for Lasso:',t1-t0)            "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Display reconstructed image and statistics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 109,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Image inversion with Lasso, SNR: 61.025452383696795\n",
      "MSE: 6.249879684849042e-06\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPsAAAD4CAYAAAAq5pAIAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAACnRJREFUeJzt3c1vjPsfxvF7OqotrXqoajRRsZBaIGJhiYiEBEktRCIisWgkIsIfQCQsWNhgxcaGWBAS8RRBkOhCUo2FBSEtovHYVj20RX+r31n1vj7O3GY6zvV+Lc+V78y0euVOzme+329udHQ0AfDfVzHeHwBAaVB2wARlB0xQdsAEZQdMTCjlm+VyOf7X/xhyuZzMo4lJVVVVajY0NFTQZ/pTJk2alJp9/fq1hJ/k36mtrZX54OBgiT7Jvzc6OjrmHxRPdsAEZQdMUHbABGUHTFB2wARlB0xQdsBErpS73lzn7GoOniTjOwuvrq6W+Y8fPzLlWdTU1Mg8+r2pny2a8VdU6Ofgr1+/ZD6emLMD5ig7YIKyAyYoO2CCsgMmKDtggrIDJpizl4HJkyfLfGRkRObDw8NFe+8vX74U/NqRfD4v82if/8+fP2Wu/rYbGxvl2rdv38q8nOfwzNkBc5QdMEHZAROUHTBB2QETlB0wUdKjpDG2aKtmtI00y1HS0XgrMmGC/hNSn72urk6u7evrk3m0Pff79++pWTRaq6+vl3l/f7/MyxFPdsAEZQdMUHbABGUHTFB2wARlB0xQdsAEc/YyMHHiRJlHc3Z15HJnZ6dc29HRIfPoauIlS5bIfN++fanZ7du35dro96Lm6EmSJFOnTk3Nohl+9NrlvMU1DU92wARlB0xQdsAEZQdMUHbABGUHTFB2wARHSZdA1r3Rc+fOlfndu3dTsxkzZsi1kyZNknkkOsb6+fPnqdm2bdvk2kePHmV6b7VXPzpCO/p+QfQdgCzHe2fFUdKAOcoOmKDsgAnKDpig7IAJyg6YoOyACebsZSCa2V6+fFnmq1evTs2ieXEkOhe+p6dH5vPnz0/NXr58KdeuXbtW5i9evJC52pMe7TdvaGiQ+fv372U+npizA+YoO2CCsgMmKDtggrIDJig7YILRWwnk83mZHzx4UOZ79uyR+bt371Kz9vZ2uXbDhg0yP3HihMynT58u83PnzqVmTU1Ncq36uZIkSRobG2WuroT+/PmzXBtdZR3l43mUNKM3wBxlB0xQdsAEZQdMUHbABGUHTFB2wMR/Zs4ebcWMrj2OqCt6o5lqc3OzzLu7uwv6TP936tSp1GzHjh1ybXV1tcyjq4ujeXNLS0tqFl0nra5cTpIkuXnzpswPHDiQmt2/f1+ujRT77y0L5uyAOcoOmKDsgAnKDpig7IAJyg6YoOyAiZLO2SsqKuSbRZ9FzXQrKyvl2pGREZlH763mqtFM9eLFizJfsWKFzKPrhefMmZOa9fb2yrXRXvvoZ4vWq+8nrFq1Sq69evWqzPv6+mSuvv+wePFiuTbaSz80NCTz8cScHTBH2QETlB0wQdkBE5QdMEHZAROUHTChN+X+YcWc6Q8PD2daH+3LzrI/ubW1VeZTpkyRebT3+s2bN6lZ1p8r+v7Cz58/Za7m7NeuXZNr9+7dK/OjR4/KXNmyZYvMjx07JvPo+wXR72U88GQHTFB2wARlB0xQdsAEZQdMUHbABGUHTPxV58ar2WbWuWaWuWl0hviNGzdkvnLlSpmfPHlS5rt3707Nvn37JtdmFc3x1Zw9Ev2bRrm6gz363Js2bZJ59G9ayl6N8d7sZwecUXbABGUHTFB2wARlB0xQdsBESbe4RqJxiNpumWWrZfTa0es3NDTItYsWLZL54OCgzCdOnChzNV6Lfqf19fUyj45rjhRzq+f27dtlrq6yjsalbW1tMr9+/brMyxFPdsAEZQdMUHbABGUHTFB2wARlB0xQdsBEWc3Zo22B6grerK+dZR7c2Ngo8xkzZhT82kmSJI8ePSp4bfRz9/f3F/zav/P6Wa66jo7YPn36tMx37tyZmkXffYj+1qqrq2X+/ft3mY8HnuyACcoOmKDsgAnKDpig7IAJyg6YoOyAibKas0eyzMKjefDIyEjBrx3N2QcGBmQezWzPnz8vc7UnPZqjR7+XmpoamUdHVWe56jq6hjvak65+9qqqKrn2w4cPMh8aGpJ5OeLJDpig7IAJyg6YoOyACcoOmKDsgAnKDpgo6Zw9Ops9mnWrs9+znhsfzZtVPm3aNLk22pf94MEDmb969Urm6rNFs+hIsa98VqI94dGsXJ0jEP29zJs3T+bjeSVzoXiyAyYoO2CCsgMmKDtggrIDJig7YIKyAyZKOmeP5ujRPeTR/mYlmvFn2Z+8bt06mUefe9myZTJfunSpzLu6ulKz6HdeW1sr8+ju+GjWrfazZ70nYMGCBTJvaGhIzfL5vFzb09Mj878RT3bABGUHTFB2wARlB0xQdsAEZQdMlNVR0lmuZI5kOSo6SfSo5uPHj3JttL02ypubm2X+8OHDgl876xbWYh6pXFdXJ/Nr167JfNasWanZs2fP5NpLly7J/G/Ekx0wQdkBE5QdMEHZAROUHTBB2QETlB0wUdI5ezTzzXK9b6SYM/zOzk6ZR8c59/b2ynzNmjUyv3fvXmoWfQcgqyxbXKMtzUeOHJH59OnTZa6cOnVK5h0dHQW/drniyQ6YoOyACcoOmKDsgAnKDpig7IAJyg6YKOmcPeusW81lo9fOOsNXV/zevXtXrh0YGJB5U1OTzKPPXsxZektLi8y7u7sLXn/8+HG5NjqiO9pLv3///tTs8OHDcm30nZDoOwLRddPjgSc7YIKyAyYoO2CCsgMmKDtggrIDJig7YCIXXZv7R98sl8v0ZmrvdDSLVnPyYrty5YrM165dK/PHjx/LfOPGjanZ69ev5dqs58bv2rVL5ocOHUrNPn/+LNfOnj1b5rdu3ZL5+vXrU7PoexnlOCf/XaOjo7mx/jtPdsAEZQdMUHbABGUHTFB2wARlB0yUdItrdKRyND5T45Jij9aqq6tTs2hM09fXJ/PBwUGZL1y4UOZdXV2p2YULF+Tap0+fylyN9ZIkvk5abQWNRmtnzpyR+datW2WuRKO3mTNnyvzdu3cFv/d44ckOmKDsgAnKDpig7IAJyg6YoOyACcoOmCjpFtd8Pi/fLJp95vP51KzYc/bJkyenZl++fMn02nfu3JH58uXLZf7+/fvUrKGhoZCP9I/ouw+fPn2SuTrmevPmzXLtkydPZB79vYyMjMhcqampkXnWrcHFxBZXwBxlB0xQdsAEZQdMUHbABGUHTFB2wMRfdZR0ZWVlahbNVNXa31mv9uJnvQ66tbVV5m1tbTJvb29PzebOnSvXRsc5R2cQqPdOkiQ5e/Zsaqa+N5EkSTI8PCzzSH19fWrW39+f6bVzuTFH2f8oZa/GeG/m7IAzyg6YoOyACcoOmKDsgAnKDpig7ICJv2rODiDGnB0wR9kBE5QdMEHZAROUHTBB2QETlB0wQdkBE5QdMEHZAROUHTBB2QETlB0wQdkBE5QdMEHZAROUHTBB2QETlB0wQdkBE5QdMEHZAROUHTBB2QETlB0wQdkBE5QdMEHZAROUHTBB2QETJb2yGcD44ckOmKDsgAnKDpig7IAJyg6YoOyACcoOmKDsgAnKDpig7IAJyg6YoOyACcoOmKDsgAnKDpig7IAJyg6YoOyACcoOmKDsgAnKDpig7IAJyg6Y+B/yhAom+zqcOwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "if mode==1 or mode==2:\n",
    "    maxi = np.max(img_rec) \n",
    "    print(\"Image inversion with Lasso, SNR: \" + str(psnr(img_np_orig,img_rec,maxi))) \n",
    "    reconstruction_err = mse(img_np_orig,img_rec,maxi)\n",
    "    print('MSE:',reconstruction_err)\n",
    "    if out_ch==3:\n",
    "        plt.imshow(img_rec.transpose(1,2,0))\n",
    "    else: \n",
    "        plt.imshow(np.clip(img_rec[0,:,:],0,1))\n",
    "        plt.gray()\n",
    "    plt.axis('off') \n",
    "    plt.show()  \n",
    "    save_path= 'lasso_'+img_name+str(int(10*f))+'.png'\n",
    "    savefig=False\n",
    "    if savefig:\n",
    "        plt.savefig(save_path,bbox_inches='tight') "
   ]
  }
 ],
 "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.6.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}