Q_PDE_Poisson_outsample10D.ipynb

{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "1kQQuxm6USWe"
      },
      "source": [
        "# BM, N = 10\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 1,
      "metadata": {
        "id": "ri5r6bp3USWf"
      },
      "outputs": [],
      "source": [
        "import numpy as np\n",
        "import torch\n",
        "import torch.nn as nn\n",
        "import torch.nn.functional as F\n",
        "import torch.optim as optim\n",
        "from torch.optim.lr_scheduler import LambdaLR\n",
        "from tqdm import tqdm, trange\n",
        "from math import exp, sqrt, log\n",
        "import matplotlib.pyplot as plt\n",
        "from mpl_toolkits.mplot3d import Axes3D\n",
        "from scipy.special import i0, i1, iv\n",
        "from numpy import random\n",
        "from torch.nn.functional import normalize\n",
        "import matplotlib.pyplot as plt\n",
        "%matplotlib qt"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 2,
      "metadata": {
        "id": "dtgNI6YnUSWg"
      },
      "outputs": [],
      "source": [
        "#from torch.utils.tensorboard import SummaryWriter\n",
        "#writer = SummaryWriter()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 3,
      "metadata": {
        "id": "K4unI4tqUSWg"
      },
      "outputs": [],
      "source": [
        "#device = torch.device(\"cuda:0\" if torch.cuda.is_available() else \"cpu\")\n",
        "#torch.set_default_tensor_type('torch.cuda.FloatTensor')\n",
        "device = torch.device(\"cpu\")"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 4,
      "metadata": {
        "id": "j6OZ0n4ZUSWh"
      },
      "outputs": [],
      "source": [
        "class PDEQnet(nn.Module):\n",
        "  def __init__(self, dim, width, beta):\n",
        "    super(PDEQnet, self).__init__()\n",
        "    self.dim = dim\n",
        "    self.width = width\n",
        "    self.innerwidth = int(width/2)\n",
        "    self.beta = beta\n",
        "\n",
        "    self.wb = nn.Linear(self.dim, self.width).to(device)\n",
        "    self.wb2 = nn.Linear(self.width, self.width).to(device)\n",
        "    self.c = nn.Linear(self.width, 1, bias=False).to(device)\n",
        "\n",
        "  def forward(self, x):\n",
        "    x = self.wb(x)\n",
        "    x = torch.sigmoid(x)\n",
        "    #x = 1/(1+torch.pow(x,2))\n",
        "    x = self.wb2(x)\n",
        "    x = torch.sigmoid(x)\n",
        "    #x = 1/(1+torch.pow(x,2))\n",
        "    x = self.c(x)\n",
        "    return x\n",
        "\n",
        "  def assign_value(self):\n",
        "    self.c.weight.data = torch.as_tensor(np.random.uniform(-1, 1, size=self.c.weight.shape), dtype=torch.float32).to(device)\n",
        "    self.wb.weight.data = torch.as_tensor(np.random.normal(0, 1, size=self.wb.weight.shape),  dtype=torch.float32).to(device)\n",
        "    self.wb.bias.data = torch.as_tensor(np.random.normal(0, 1, size=self.wb.bias.shape) ,dtype=torch.float32).to(device)\n",
        "    self.wb2.weight.data = torch.as_tensor(np.random.normal(0, 1, size=self.wb2.weight.shape),  dtype=torch.float32).to(device)\n",
        "    self.wb2.bias.data = torch.as_tensor(np.random.normal(0, 1, size=self.wb2.bias.shape) ,dtype=torch.float32).to(device)\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 5,
      "metadata": {
        "id": "QKXV8Kj3USWi"
      },
      "outputs": [],
      "source": [
        "# Parameters\n",
        "dim = 10\n",
        "gamma = 0.2\n",
        "\n",
        "#d = [1 for i in range(dim)]\n",
        "\n",
        "# Hyper parameters\n",
        "N = 2**6\n",
        "beta = 0.5+0.01\n",
        "# Learning rate\n",
        "initial_lr = 0.02 #* N**(2*beta-1)\n",
        "\n",
        "# Auxiliary functions\n",
        "def eta(grid):\n",
        "  res = 1.0\n",
        "  for i in range(dim):\n",
        "    res = res -  grid[:,i]*grid[:,i]\n",
        "  return res\n",
        "\n",
        "def eta_diff(grid, count):\n",
        "  return -2*grid[:,count]\n",
        "\n",
        "def eta_diff_2(grid, count):\n",
        "  return -2\n",
        "\n",
        "\n",
        "# Monte Carlo\n",
        "Nmc_max = 2000\n",
        "Nmc_initial = 500\n",
        "Nbasepoints = int(1e7)\n",
        "# Default type\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 6,
      "metadata": {
        "id": "pAobmnpKUSWi"
      },
      "outputs": [],
      "source": [
        "# Q fit, fixed grid\n",
        "qnet = PDEQnet(dim, N, beta).to(device)\n",
        "# initialization of PDEQnet paramters\n",
        "qnet.assign_value()\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 7,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "w4FeS4VbUSWj",
        "outputId": "078962cc-a791-4aa7-95a1-4baa24cb3dc0"
      },
      "outputs": [],
      "source": [
        "#qnet.load_state_dict(torch.load('/content/6d, Q-learning, 30k epochs.pkl'))"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 8,
      "metadata": {},
      "outputs": [],
      "source": [
        "# Create reservoir of points, so we don't have to sample every loop\n",
        "#this is maybe the one thing that's worth doing on a GPU\n",
        "\n",
        "source = torch.randn(size=(Nbasepoints, dim))\n",
        "source = normalize(source, p=2.0) #Normalize to sphere\n",
        "radius = torch.rand(size = (Nbasepoints,1))\n",
        "radius = torch.pow(torch.rand(size = (Nbasepoints,1)),1/dim)\n",
        "source = radius*source #renormalize\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 9,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "dk8QAmWgUSWk",
        "outputId": "d23ac1fa-7eb8-4254-b3ba-f252211ad845"
      },
      "outputs": [
        {
          "name": "stderr",
          "output_type": "stream",
          "text": [
            " 39%|███▉      | 11659/30000 [3:51:13<6:03:44,  1.19s/it]      \n"
          ]
        },
        {
          "ename": "KeyboardInterrupt",
          "evalue": "",
          "output_type": "error",
          "traceback": [
            "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
            "\u001b[1;31mKeyboardInterrupt\u001b[0m                         Traceback (most recent call last)",
            "\u001b[1;32m<ipython-input-9-23647ca948b2>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m\u001b[0m\n\u001b[0;32m     36\u001b[0m   \u001b[0ml2\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;36m0.0\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     37\u001b[0m   \u001b[1;32mfor\u001b[0m \u001b[0mcount\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mentry\u001b[0m \u001b[1;32min\u001b[0m \u001b[0menumerate\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mgrid_list\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m---> 38\u001b[1;33m     \u001b[0mdf_dx\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mtorch\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mautograd\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mgrad\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mout\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mentry\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mgrad_outputs\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mtorch\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mones\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mentry\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0msize\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mcreate_graph\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;32mTrue\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m     39\u001b[0m     \u001b[0md2f_dx2\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mtorch\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mautograd\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mgrad\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mdf_dx\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mentry\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mgrad_outputs\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mtorch\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mones\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mentry\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0msize\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mcreate_graph\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;32mTrue\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     40\u001b[0m     \u001b[0mdf_dx\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mtorch\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mreshape\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mdf_dx\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;33m(\u001b[0m\u001b[1;33m-\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
            "\u001b[1;32m~\\AppData\\Local\\Programs\\Python\\Python38\\lib\\site-packages\\torch\\autograd\\__init__.py\u001b[0m in \u001b[0;36mgrad\u001b[1;34m(outputs, inputs, grad_outputs, retain_graph, create_graph, only_inputs, allow_unused, is_grads_batched)\u001b[0m\n\u001b[0;32m    273\u001b[0m         \u001b[1;32mreturn\u001b[0m \u001b[0m_vmap_internals\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m_vmap\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mvjp\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;36m0\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;36m0\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mallow_none_pass_through\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;32mTrue\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mgrad_outputs\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    274\u001b[0m     \u001b[1;32melse\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 275\u001b[1;33m         return Variable._execution_engine.run_backward(  # Calls into the C++ engine to run the backward pass\n\u001b[0m\u001b[0;32m    276\u001b[0m             \u001b[0moutputs\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mgrad_outputs_\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mretain_graph\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mcreate_graph\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0minputs\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    277\u001b[0m             allow_unused, accumulate_grad=False)  # Calls into the C++ engine to run the backward pass\n",
            "\u001b[1;31mKeyboardInterrupt\u001b[0m: "
          ]
        }
      ],
      "source": [
        "#Training algorithm, main\n",
        "\n",
        "# Num of epoch\n",
        "Num_epo = 30000\n",
        "\n",
        "# Loss level\n",
        "loss_list = []\n",
        "\n",
        "# Optmizer, scheduler\n",
        "Qoptimizer = optim.RMSprop(qnet.parameters(), lr = initial_lr, alpha = 0.99, eps = 1e-08)\n",
        "Qscheduler = LambdaLR(Qoptimizer, lr_lambda= lambda epoch: initial_lr/(1+(epoch//1000)))\n",
        "\n",
        "#fig, axs = plt.subplots(2,1)\n",
        "#fig,ax =plt.subplots()\n",
        "# Training\n",
        "for count in trange(Num_epo):\n",
        "\n",
        "  Nmc = int(Nmc_initial+(Nmc_max-Nmc_initial)*(1+count)/(1+Num_epo))\n",
        "\n",
        "  epoch_sample = torch.randint(0,len(source), (Nmc,1))\n",
        "\n",
        "  grid_list = []\n",
        "  for i in range(dim):\n",
        "    ent = [[u] for u in source[epoch_sample,i]] #grab sample from reservoir\n",
        "    ent = torch.tensor(ent, requires_grad = True)\n",
        "    grid_list.append(ent)\n",
        "  grid_tuple = tuple(grid_list)\n",
        "  grid = torch.cat(grid_tuple, 1)\n",
        "  \n",
        "  # Net output\n",
        "  out = qnet(grid)\n",
        "\n",
        "  # Tensor reshape\n",
        "  out_r =  torch.reshape(out, (-1,))\n",
        "  # Compute partial derivatives and the Laplacian\n",
        "  l2 = 0.0\n",
        "  for count, entry in enumerate(grid_list):\n",
        "    df_dx = torch.autograd.grad(out, entry, grad_outputs = torch.ones(entry.size()), create_graph = True)[0]\n",
        "    d2f_dx2 = torch.autograd.grad(df_dx, entry, grad_outputs = torch.ones(entry.size()), create_graph = True)[0]\n",
        "    df_dx = torch.reshape(df_dx, (-1,))\n",
        "    d2f_dx2 = torch.reshape(d2f_dx2, (-1,))\n",
        "    l2 += d2f_dx2*eta(grid) + 2*df_dx*eta_diff(grid, count) + out_r* eta_diff_2(grid, count)\n",
        "  \n",
        "  #evaluate function\n",
        "  l = out_r*eta(grid)\n",
        "  #evaluate PDE operator\n",
        "  lq = 1 - gamma*l + 0.5*l2\n",
        "  LQ = lq.clone().detach()\n",
        "\n",
        "  # Q-learning\n",
        "  loss_to_min =-1* torch.dot(LQ, l)\n",
        "  #with torch.cuda.stream(s):\n",
        "  Qoptimizer.zero_grad()\n",
        "  loss_to_min.backward()\n",
        "  Qoptimizer.step()\n",
        "  Qscheduler.step()\n",
        "\n",
        "  loss = float(torch.dot(LQ, LQ))\n",
        "  loss/= Nmc\n",
        "  loss_list.append(loss)\n",
        "  \n",
        "  if count%10 == 0:\n",
        "    print(loss)\n",
        "  #  ax.clear()\n",
        "  #  axs[1].clear()\n",
        "  #  axs[0].scatter(grid.detach().cpu(),l.detach().cpu())\n",
        "  #  ax.scatter(grid[:,0].detach().cpu(), LQ.detach().cpu())\n",
        "  #  plt.pause(0.01)\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "lT6HSkdEUSWm"
      },
      "outputs": [],
      "source": [
        "model_save_name = 'BM, 15d, 30k epochs.pkl'\n",
        "path = F\"/content/{model_save_name}\" \n",
        "torch.save(qnet.state_dict(), path)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 307
        },
        "id": "w2IaSkcpUSWm",
        "outputId": "45aee8d8-4ee4-486a-edcc-192a4d2b18e0"
      },
      "outputs": [],
      "source": [
        "# Loss level\n",
        "plt.figure(figsize=(7,4.5))\n",
        "#ax = fig.add_subplot(1, 2, 1)\n",
        "axis=[i for i in range(Num_epo)]\n",
        "axis = np.cumsum([1/(1+(i//500)) for i in range(Num_epo)])\n",
        "plt.xlabel('Cumulative training time')\n",
        "plt.ylabel('Loss level')\n",
        "plt.yscale('log') \n",
        "fig1 = plt.plot(axis,loss_list,'blue')"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 19,
      "metadata": {},
      "outputs": [],
      "source": [
        "# Loss level\n",
        "plt.figure(figsize=(7,4.5))\n",
        "#ax = fig.add_subplot(1, 2, 1)\n",
        "axis=[i for i in range(len(loss_list))]\n",
        "plt.xlabel('Training epoch')\n",
        "plt.ylabel('Loss level')\n",
        "plt.yscale('log') \n",
        "fig1 = plt.plot(axis,loss_list,'blue')"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 307
        },
        "id": "D_Xm8eAbUSWm",
        "outputId": "49cc20b7-58aa-426d-b3d1-99f9e419e639"
      },
      "outputs": [],
      "source": [
        "log_loss_list = [log(x) for x in loss_list]\n",
        "plt.figure(figsize=(7,4.5))\n",
        "#ax = fig.add_subplot(1, 2, 1)\n",
        "axis=[i for i in range(Num_epo)]\n",
        "plt.xlabel('Number of epochs')\n",
        "plt.ylabel('Loss level')\n",
        "fig1 = plt.plot(axis[400:600], log_loss_list[400:600],'blue')"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 16,
      "metadata": {
        "id": "Kseai8SHUSWq"
      },
      "outputs": [],
      "source": [
        "# Plotting part 1: # Check if the approximator is symmetric\n",
        "Nmc = 400\n",
        "mesh = 1000\n",
        "\n",
        "test_axis = [i/mesh for i in range(1,mesh+1)]\n",
        "relative_err_list = []\n",
        "approx_list = []\n",
        "sq_list = []\n",
        "exact_list = []\n",
        "\n",
        "radius = 1.0\n",
        "for x in test_axis:\n",
        "  exact = (1 - iv(9, sqrt(2*gamma)*x)/(x**9*iv(9, sqrt(2*gamma))))/gamma\n",
        "  exact_list.append(exact)\n",
        "  test_source = random.normal(0, 1, size=(Nmc//10, dim))\n",
        "  test_source = torch.Tensor(test_source)\n",
        "  test_source = normalize(test_source, p=2.0)\n",
        "  test_source = x * test_source\n",
        "\n",
        "  test_grid_list = []\n",
        "  for i in range(dim):\n",
        "    ent = [[u] for u in test_source[:,i]]\n",
        "    ent = torch.tensor(ent, requires_grad = True).to(device)\n",
        "    test_grid_list.append(ent)\n",
        "  test_grid_tuple = tuple(test_grid_list)\n",
        "  test_grid = torch.cat(test_grid_tuple, 1).to(device)\n",
        "\n",
        "  test_out = qnet(test_grid).to(device)\n",
        "\n",
        "# Tensor reshape\n",
        "  test_out_r =  torch.reshape(test_out, (-1,)).to(device)\n",
        "  test_l = test_out_r*eta(test_grid)\n",
        "\n",
        "  sq_list.append(float(torch.dot(test_l-exact, test_l-exact)))\n",
        "  relative_err_list.append(float(torch.max(abs(exact - test_l)/exact)))\n",
        "  #relative_err_list.append(100*float(torch.mean(abs(exact - test_l))))\n",
        "  approx_list.append(float(torch.mean(test_l)))"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 17,
      "metadata": {
        "id": "goz-5zTXUSWr"
      },
      "outputs": [],
      "source": [
        "# Plotting part 2\n",
        "# Scatter plot #Test: the solution vs the approximator\n",
        "scatter_axis = test_axis * (Nmc//10)\n",
        "relative_err_list = []\n",
        "approx_list = []\n",
        "std_list = []\n",
        "\n",
        "\n",
        "for x in scatter_axis:\n",
        "  test_source = random.normal(0, 1, size=(1, dim))\n",
        "  test_source = torch.Tensor(test_source)\n",
        "  test_source = normalize(test_source, p=2.0)\n",
        "  test_source = x * test_source\n",
        "\n",
        "  test_grid_list = []\n",
        "  for i in range(dim):\n",
        "    ent = [[u] for u in test_source[:,i]]\n",
        "    ent = torch.tensor(ent, requires_grad = True).to(device)\n",
        "    test_grid_list.append(ent)\n",
        "  test_grid_tuple = tuple(test_grid_list)\n",
        "  test_grid = torch.cat(test_grid_tuple, 1).to(device)\n",
        "\n",
        "  test_out = qnet(test_grid).to(device)\n",
        "\n",
        "# Tensor reshape\n",
        "  test_out_r =  torch.reshape(test_out, (-1,)).to(device)\n",
        "  test_l = test_out_r*eta(test_grid)\n",
        "  approx_list.append(float(torch.mean(test_l)))\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 18,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 386
        },
        "id": "xpG0aEc_USWr",
        "outputId": "bed5dbd3-676d-47bd-907a-0ffbf1dd8c1d"
      },
      "outputs": [
        {
          "data": {
            "text/plain": [
              "(<matplotlib.patches.Rectangle at 0x164b2f59c40>,\n",
              " (<matplotlib.patches.ConnectionPatch at 0x164b2f07ee0>,\n",
              "  <matplotlib.patches.ConnectionPatch at 0x164b2f6edf0>,\n",
              "  <matplotlib.patches.ConnectionPatch at 0x164b2f59040>,\n",
              "  <matplotlib.patches.ConnectionPatch at 0x164b2f20550>))"
            ]
          },
          "execution_count": 18,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "# Exact vs Approximate\n",
        "fig = plt.figure()\n",
        "ax = fig.add_subplot()\n",
        "ax.set_xlabel('Distance from the Origin')\n",
        "ax.set_ylabel('Value function')\n",
        "ax.set_title('Solution of the PDE')\n",
        "ax.scatter(scatter_axis, approx_list, s= .2, label = 'Approximated solution')\n",
        "ax.plot(test_axis, exact_list,'red', label = 'Exact solution')\n",
        "ax.legend()\n",
        "axins = ax.inset_axes([0.05,0.05, 0.5, 0.5])\n",
        "scatter_data_zoom = [approx_list[i] for i in range(len(scatter_axis)) if (0.48 < scatter_axis[i]  < 0.51)]\n",
        "scatter_axis_zoom = [item for item in scatter_axis if (0.48 < item < 0.51)]\n",
        "line_data_zoom = [exact_list[i] for i in range(len(test_axis)) if (0.48 < test_axis[i] < 0.51)]\n",
        "line_axis_zoom = [item for item in test_axis if (0.48 < item < 0.51)]\n",
        "axins.scatter(scatter_axis_zoom, scatter_data_zoom, s=.1)\n",
        "axins.plot(line_axis_zoom, line_data_zoom, 'red')\n",
        "axins.set_xticks([])\n",
        "axins.set_yticks([])\n",
        "ax.indicate_inset_zoom(axins)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 313
        },
        "id": "xW8Sgs6-USWs",
        "outputId": "3d81ed9d-19c1-48d4-9bb4-d407d4ccf9c6"
      },
      "outputs": [],
      "source": [
        "# Mean square error\n",
        "plt.xlabel('Distance from the Origin')\n",
        "plt.ylabel('Mean square error')\n",
        "plt.scatter(test_axis, sq_list,color = 'blue', s= .2)\n",
        "plt.title('Mean square error of the approximator')"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 313
        },
        "id": "1V_Z4lt8USWt",
        "outputId": "1583f032-024e-4b28-ad6b-1e8ade85453f"
      },
      "outputs": [],
      "source": [
        "# Max relative error\n",
        "plt.xlabel('Distance from the Origin')\n",
        "plt.ylabel('Percentage')\n",
        "plt.scatter(test_axis, relative_err_list, s=0.2, color='gray',)\n",
        "plt.gca().set_yticklabels([f'{x:.2%}' for x in plt.gca().get_yticks()]) \n",
        "plt.title('MaxRelative error in percentage')"
      ]
    }
  ],
  "metadata": {
    "accelerator": "GPU",
    "colab": {
      "collapsed_sections": [],
      "name": "Q-PDE playground.ipynb",
      "provenance": []
    },
    "kernelspec": {
      "display_name": "Python 3",
      "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.8.0"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}