{ "cells": [ { "cell_type": "markdown", "id": "4fad6f64-ece7-46a4-bc17-0b928b19053a", "metadata": {}, "source": [ "# Solving the heat equation\n", "\n", "In this notebook we will learn how to use tensorized linear solvers for solving a\n", "partial differential equation (PDE). The PDE in question will be the heat equation\n", "$\\Delta x = b$. $\\Delta$ is the laplace operator and $b$ some boundary condition.\n", "We will solve the equation in question with finite difference on a two dimensional\n", "grid." ] }, { "cell_type": "code", "execution_count": 1, "id": "4564d77b-314a-4765-89a8-1e820b42f5a4", "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "from trainsum.numpy import trainsum as ts\n", "from trainsum.typing import UniformGrid\n", "import matplotlib.pyplot as plt\n", "\n", "def laplace(grid: UniformGrid, idx: int) -> ts.TensorTrain[NDArray]:\n", " \"\"\"\n", " Get the tensor train which represents the finite difference\n", " laplace operator with the dimension defined by idx.\n", " \"\"\"\n", " train = None\n", " with ts.exact():\n", " for i, dim in enumerate(grid.dims):\n", " if i == idx:\n", " tmp = -2.0 * ts.shift(dim, 0)\n", " tmp += 1.0 * ts.shift(dim, -1)\n", " tmp += 1.0 * ts.shift(dim, 1)\n", " tmp *= -1/grid.spacings[0]**2\n", " else:\n", " tmp = ts.shift(dim, 0)\n", " if train is None:\n", " train = tmp\n", " else:\n", " train.extend(tmp)\n", " return train" ] }, { "cell_type": "markdown", "id": "c2096489-7164-46d7-a6fa-62bf87ad0fba", "metadata": {}, "source": [ "As a first step we define the laplace operator as a tridiagonal toeplitz\n", "matrix with -2 on the main diagonal and 1 on the diagonals above and below." ] }, { "cell_type": "code", "execution_count": 2, "id": "25c045bd-5eb8-448d-9658-da69352eab76", "metadata": {}, "outputs": [], "source": [ "# define the domain on which the problem is solved\n", "shape = ts.trainshape(32, 32, mode=\"block\")\n", "domains = ts.domain(0.0, 1.0), ts.domain(0.0, 1.0)\n", "grid = ts.uniform_grid(shape.dims, domains)\n", "\n", "xshape = ts.trainshape(shape.dims[0])\n", "yshape = ts.trainshape(shape.dims[1])" ] }, { "cell_type": "markdown", "id": "72ade7b4-7f69-42e3-97da-f79992ba099e", "metadata": {}, "source": [ "Next we define solution space with an uniformly spaced grid. For doing so we will \n", "need some dimensions which are associated to some intervals (domains)." ] }, { "cell_type": "code", "execution_count": 3, "id": "de7869f1-a42f-4ca8-8f0c-60e025865f63", "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaAAAAGdCAYAAABU0qcqAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjgsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvwVt1zgAAAAlwSFlzAAAPYQAAD2EBqD+naQAAHKdJREFUeJzt3Xts1fUd//HX6e0A0p5aSnvaUVgBBRXpsk7qicpQOkqXmCKY4GVZcQQDK2bQOZXF67akDhOvQfhjv8lMRByLQDQ/cVpsiVtho5MgOvujrBsl0KIsPacUeyjt5/fH4tmOUOhpT3lzyvORfBN6zud8z/vLl/bp6bnocc45AQBwkSVZDwAAuDwRIACACQIEADBBgAAAJggQAMAEAQIAmCBAAAATBAgAYCLFeoCv6+vr09GjR5Weni6Px2M9DgAgRs45dXZ2Kj8/X0lJ/T/OueQCdPToURUUFFiPAQAYotbWVk2YMKHf64ctQOvWrdMzzzyjtrY2FRUV6aWXXtKsWbMueLv09HRJ0pz8pUpJShvQfblR3gHP1Zc+asBrJelM+sBmkKSwLzWmfYczB/4b0O6s2B4NhrMG/glLZ8b1xLRvX3ZXTOunZH0x4LXXpR+Nad8zRh0Z8Nqr0wY+hyT5k5MHvHZsUmz/rmJ1sq97wGvbentj2vf/O5094LUHuvv/YXIun3TmD3jtoX8PfA5JCn5xxYDXppyI7XvT++/Yvt9G/Xvg32/ejr7YZgkO/PszpfN0TPtO6hz4vytPd3jAa8/0nVbd0f8T+Xnen2EJ0BtvvKHq6mpt2LBBJSUlev7551VWVqampibl5OSc97Zf/dotJSlNKUkDC4tLjiFAMaz9zyADX9+bGts/8jNpAw9Qsje2b4ikUQP/hkgaPfAftJKUPOZMTOtTrxh4xEeNje3vcEwMs4+N4e9bkjKSB75+7Hl+zRAPSX0D3//J3tg+3nHM6YH/HY5Kie38pPYN/Nwnd8f2vZk0euDRTxoV29yxfr8lpw387zwlNbYApaQM/PykJMf4cyJ54HN7BvFP/EJPowzLd82zzz6rZcuW6b777tO1116rDRs2aMyYMfrtb387HHcHAEhAcQ/Q6dOn1djYqNLS0v/eSVKSSktL1dDQcNb6cDisUCgUtQEARr64B+iLL75Qb2+vcnNzoy7Pzc1VW1vbWetramrk8/kiGy9AAIDLg/n7gNasWaNgMBjZWltbrUcCAFwEcX8RQnZ2tpKTk9Xe3h51eXt7u/x+/1nrvV6vvN4YXxgAAEh4cX8ElJaWpuLiYtXW1kYu6+vrU21trQKBQLzvDgCQoIblZdjV1dWqrKzUd77zHc2aNUvPP/+8urq6dN999w3H3QEAEtCwBGjx4sX6/PPP9fjjj6utrU3f+ta3tGPHjrNemAAAuHx5nHOxvWttmIVCIfl8Ps1RhVI8sb15DABg74zrUZ22KxgMKiMjo9915q+CAwBcnggQAMAEAQIAmCBAAAATBAgAYIIAAQBMECAAgAkCBAAwQYAAACYIEADABAECAJggQAAAEwQIAGCCAAEATBAgAIAJAgQAMEGAAAAmCBAAwAQBAgCYIEAAABMECABgggABAEwQIACACQIEADBBgAAAJggQAMAEAQIAmCBAAAATBAgAYIIAAQBMECAAgAkCBAAwQYAAACYIEADABAECAJggQAAAEwQIAGCCAAEATBAgAIAJAgQAMEGAAAAmCBAAwAQBAgCYIEAAABMECABgggABAEwQIACAibgH6Mknn5TH44napk+fHu+7AQAkuJTh2Ol1112n999//793kjIsdwMASGDDUoaUlBT5/f7h2DUAYIQYlueADh48qPz8fE2ePFn33nuvDh8+3O/acDisUCgUtQEARr64B6ikpEQbN27Ujh07tH79erW0tOiWW25RZ2fnOdfX1NTI5/NFtoKCgniPBAC4BHmcc24476Cjo0OTJk3Ss88+q6VLl551fTgcVjgcjnwdCoVUUFCgOapQiid1OEcDAAyDM65HddquYDCojIyMftcN+6sDMjMzdfXVV6u5ufmc13u9Xnm93uEeAwBwiRn29wGdPHlShw4dUl5e3nDfFQAggcQ9QA8++KDq6+v1z3/+U3/+8591xx13KDk5WXfffXe87woAkMDi/iu4I0eO6O6779aJEyc0fvx43Xzzzdq9e7fGjx8f77sCACSwuAdo8+bN8d4lAGAE4rPgAAAmCBAAwAQBAgCYIEAAABMECABgggABAEwQIACACQIEADBBgAAAJggQAMAEAQIAmCBAAAATBAgAYIIAAQBMECAAgAkCBAAwQYAAACYIEADABAECAJggQAAAEwQIAGCCAAEATBAgAIAJAgQAMEGAAAAmCBAAwAQBAgCYIEAAABMECABgggABAEwQIACACQIEADBBgAAAJggQAMAEAQIAmCBAAAATBAgAYIIAAQBMECAAgAkCBAAwQYAAACYIEADABAECAJggQAAAEwQIAGCCAAEATBAgAICJmAO0a9cu3X777crPz5fH49G2bduirnfO6fHHH1deXp5Gjx6t0tJSHTx4MF7zAgBGiJgD1NXVpaKiIq1bt+6c169du1YvvviiNmzYoD179uiKK65QWVmZuru7hzwsAGDkSIn1BuXl5SovLz/ndc45Pf/883r00UdVUVEhSXr11VeVm5urbdu26a677hratACAESOuzwG1tLSora1NpaWlkct8Pp9KSkrU0NBwztuEw2GFQqGoDQAw8sU1QG1tbZKk3NzcqMtzc3Mj131dTU2NfD5fZCsoKIjnSACAS5T5q+DWrFmjYDAY2VpbW61HAgBcBHENkN/vlyS1t7dHXd7e3h657uu8Xq8yMjKiNgDAyBfXABUWFsrv96u2tjZyWSgU0p49exQIBOJ5VwCABBfzq+BOnjyp5ubmyNctLS3at2+fsrKyNHHiRK1atUq/+tWvdNVVV6mwsFCPPfaY8vPztWDBgnjODQBIcDEHaO/evbr11lsjX1dXV0uSKisrtXHjRj300EPq6urS/fffr46ODt18883asWOHRo0aFb+pAQAJz+Occ9ZD/K9QKCSfz6c5qlCKJ9V6HABAjM64HtVpu4LB4Hmf1zd/FRwA4PJEgAAAJggQAMAEAQIAmCBAAAATBAgAYIIAAQBMECAAgAkCBAAwQYAAACYIEADABAECAJggQAAAEwQIAGCCAAEATBAgAIAJAgQAMEGAAAAmCBAAwAQBAgCYIEAAABMECABgggABAEwQIACACQIEADBBgAAAJggQAMAEAQIAmCBAAAATBAgAYIIAAQBMECAAgAkCBAAwQYAAACYIEADABAECAJggQAAAEwQIAGCCAAEATBAgAIAJAgQAMEGAAAAmCBAAwAQBAgCYIEAAABMECABgggABAEzEHKBdu3bp9ttvV35+vjwej7Zt2xZ1/ZIlS+TxeKK2+fPnx2teAMAIEXOAurq6VFRUpHXr1vW7Zv78+Tp27Fhke/3114c0JABg5EmJ9Qbl5eUqLy8/7xqv1yu/3z/ooQAAI9+wPAdUV1ennJwcTZs2TStWrNCJEyf6XRsOhxUKhaI2AMDIF/cAzZ8/X6+++qpqa2v161//WvX19SovL1dvb+8519fU1Mjn80W2goKCeI8EALgEeZxzbtA39ni0detWLViwoN81//jHPzRlyhS9//77mjt37lnXh8NhhcPhyNehUEgFBQWaowqleFIHOxoAwMgZ16M6bVcwGFRGRka/64b9ZdiTJ09Wdna2mpubz3m91+tVRkZG1AYAGPmGPUBHjhzRiRMnlJeXN9x3BQBIIDG/Cu7kyZNRj2ZaWlq0b98+ZWVlKSsrS0899ZQWLVokv9+vQ4cO6aGHHtLUqVNVVlYW18EBAIkt5gDt3btXt956a+Tr6upqSVJlZaXWr1+v/fv363e/+506OjqUn5+vefPm6Ze//KW8Xm/8pgYAJLyYAzRnzhyd73UL77777pAGAgBcHvgsOACACQIEADBBgAAAJggQAMAEAQIAmCBAAAATBAgAYIIAAQBMECAAgAkCBAAwQYAAACYIEADABAECAJggQAAAEwQIAGCCAAEATBAgAIAJAgQAMEGAAAAmCBAAwAQBAgCYIEAAABMECABgggABAEwQIACACQIEADBBgAAAJggQAMAEAQIAmCBAAAATBAgAYIIAAQBMECAAgAkCBAAwQYAAACYIEADABAECAJggQAAAEwQIAGCCAAEATBAgAIAJAgQAMEGAAAAmCBAAwAQBAgCYIEAAABMxBaimpkY33HCD0tPTlZOTowULFqipqSlqTXd3t6qqqjRu3DiNHTtWixYtUnt7e1yHBgAkvpgCVF9fr6qqKu3evVvvvfeeenp6NG/ePHV1dUXWrF69Wm+99Za2bNmi+vp6HT16VAsXLoz74ACAxOZxzrnB3vjzzz9XTk6O6uvrNXv2bAWDQY0fP16bNm3SnXfeKUn67LPPdM0116ihoUE33njjBfcZCoXk8/k0RxVK8aQOdjQAgJEzrkd12q5gMKiMjIx+1w3pOaBgMChJysrKkiQ1Njaqp6dHpaWlkTXTp0/XxIkT1dDQcM59hMNhhUKhqA0AMPINOkB9fX1atWqVbrrpJs2YMUOS1NbWprS0NGVmZkatzc3NVVtb2zn3U1NTI5/PF9kKCgoGOxIAIIEMOkBVVVU6cOCANm/ePKQB1qxZo2AwGNlaW1uHtD8AQGJIGcyNVq5cqbffflu7du3ShAkTIpf7/X6dPn1aHR0dUY+C2tvb5ff7z7kvr9crr9c7mDEAAAkspkdAzjmtXLlSW7du1c6dO1VYWBh1fXFxsVJTU1VbWxu5rKmpSYcPH1YgEIjPxACAESGmR0BVVVXatGmTtm/frvT09MjzOj6fT6NHj5bP59PSpUtVXV2trKwsZWRk6IEHHlAgEBjQK+AAAJePmAK0fv16SdKcOXOiLn/llVe0ZMkSSdJzzz2npKQkLVq0SOFwWGVlZXr55ZfjMiwAYOQY0vuAhgPvAwKAxHZR3gcEAMBgESAAgAkCBAAwQYAAACYIEADABAECAJggQAAAEwQIAGCCAAEATBAgAIAJAgQAMEGAAAAmCBAAwAQBAgCYIEAAABMECABgggABAEwQIACACQIEADBBgAAAJggQAMAEAQIAmCBAAAATBAgAYIIAAQBMECAAgAkCBAAwQYAAACYIEADABAECAJggQAAAEwQIAGCCAAEATBAgAIAJAgQAMEGAAAAmCBAAwAQBAgCYIEAAABMECABgggABAEwQIACACQIEADBBgAAAJggQAMAEAQIAmCBAAAATMQWopqZGN9xwg9LT05WTk6MFCxaoqakpas2cOXPk8XiituXLl8d1aABA4ospQPX19aqqqtLu3bv13nvvqaenR/PmzVNXV1fUumXLlunYsWORbe3atXEdGgCQ+FJiWbxjx46orzdu3KicnBw1NjZq9uzZkcvHjBkjv98fnwkBACPSkJ4DCgaDkqSsrKyoy1977TVlZ2drxowZWrNmjU6dOtXvPsLhsEKhUNQGABj5YnoE9L/6+vq0atUq3XTTTZoxY0bk8nvuuUeTJk1Sfn6+9u/fr4cfflhNTU168803z7mfmpoaPfXUU4MdAwCQoDzOOTeYG65YsULvvPOOPvzwQ02YMKHfdTt37tTcuXPV3NysKVOmnHV9OBxWOByOfB0KhVRQUKA5qlCKJ3UwowEADJ1xParTdgWDQWVkZPS7blCPgFauXKm3335bu3btOm98JKmkpESS+g2Q1+uV1+sdzBgAgAQWU4Ccc3rggQe0detW1dXVqbCw8IK32bdvnyQpLy9vUAMCAEammAJUVVWlTZs2afv27UpPT1dbW5skyefzafTo0Tp06JA2bdqk73//+xo3bpz279+v1atXa/bs2Zo5c+awHAAAIDHF9ByQx+M55+WvvPKKlixZotbWVv3gBz/QgQMH1NXVpYKCAt1xxx169NFHz/t7wP8VCoXk8/l4DggAEtSwPAd0oVYVFBSovr4+ll0CAC5TfBYcAMAEAQIAmCBAAAATBAgAYIIAAQBMECAAgAkCBAAwQYAAACYIEADABAECAJggQAAAEwQIAGCCAAEATBAgAIAJAgQAMEGAAAAmCBAAwAQBAgCYIEAAABMECABgggABAEwQIACACQIEADBBgAAAJggQAMAEAQIAmCBAAAATBAgAYIIAAQBMECAAgAkCBAAwQYAAACYIEADABAECAJggQAAAEwQIAGCCAAEATBAgAIAJAgQAMEGAAAAmCBAAwAQBAgCYIEAAABMECABgggABAEwQIACACQIEADARU4DWr1+vmTNnKiMjQxkZGQoEAnrnnXci13d3d6uqqkrjxo3T2LFjtWjRIrW3t8d9aABA4ospQBMmTNDTTz+txsZG7d27V7fddpsqKir0ySefSJJWr16tt956S1u2bFF9fb2OHj2qhQsXDsvgAIDE5nHOuaHsICsrS88884zuvPNOjR8/Xps2bdKdd94pSfrss890zTXXqKGhQTfeeOOA9hcKheTz+TRHFUrxpA5lNACAgTOuR3XarmAwqIyMjH7XDfo5oN7eXm3evFldXV0KBAJqbGxUT0+PSktLI2umT5+uiRMnqqGhod/9hMNhhUKhqA0AMPLFHKCPP/5YY8eOldfr1fLly7V161Zde+21amtrU1pamjIzM6PW5+bmqq2trd/91dTUyOfzRbaCgoKYDwIAkHhiDtC0adO0b98+7dmzRytWrFBlZaU+/fTTQQ+wZs0aBYPByNba2jrofQEAEkdKrDdIS0vT1KlTJUnFxcX661//qhdeeEGLFy/W6dOn1dHREfUoqL29XX6/v9/9eb1eeb3e2CcHACS0Ib8PqK+vT+FwWMXFxUpNTVVtbW3kuqamJh0+fFiBQGCodwMAGGFiegS0Zs0alZeXa+LEiers7NSmTZtUV1end999Vz6fT0uXLlV1dbWysrKUkZGhBx54QIFAYMCvgAMAXD5iCtDx48f1wx/+UMeOHZPP59PMmTP17rvv6nvf+54k6bnnnlNSUpIWLVqkcDissrIyvfzyy8MyOAAgsQ35fUDxxvuAACCxDfv7gAAAGAoCBAAwQYAAACYIEADABAECAJggQAAAEwQIAGCCAAEATBAgAICJmD8Ne7h99cEMZ9QjXVKf0QAAGIgz6pH035/n/bnkAtTZ2SlJ+lD/13gSAMBQdHZ2yufz9Xv9JfdZcH19fTp69KjS09Pl8Xgil4dCIRUUFKi1tfW8ny2U6DjOkeNyOEaJ4xxp4nGczjl1dnYqPz9fSUn9P9NzyT0CSkpK0oQJE/q9PiMjY0Sf/K9wnCPH5XCMEsc50gz1OM/3yOcrvAgBAGCCAAEATCRMgLxer5544gl5vV7rUYYVxzlyXA7HKHGcI83FPM5L7kUIAIDLQ8I8AgIAjCwECABgggABAEwQIACAiYQJ0Lp16/TNb35To0aNUklJif7yl79YjxRXTz75pDweT9Q2ffp067GGZNeuXbr99tuVn58vj8ejbdu2RV3vnNPjjz+uvLw8jR49WqWlpTp48KDNsENwoeNcsmTJWed2/vz5NsMOUk1NjW644Qalp6crJydHCxYsUFNTU9Sa7u5uVVVVady4cRo7dqwWLVqk9vZ2o4kHZyDHOWfOnLPO5/Lly40mHpz169dr5syZkTebBgIBvfPOO5HrL9a5TIgAvfHGG6qurtYTTzyhv/3tbyoqKlJZWZmOHz9uPVpcXXfddTp27Fhk+/DDD61HGpKuri4VFRVp3bp157x+7dq1evHFF7Vhwwbt2bNHV1xxhcrKytTd3X2RJx2aCx2nJM2fPz/q3L7++usXccKhq6+vV1VVlXbv3q333ntPPT09mjdvnrq6uiJrVq9erbfeektbtmxRfX29jh49qoULFxpOHbuBHKckLVu2LOp8rl271mjiwZkwYYKefvppNTY2au/evbrttttUUVGhTz75RNJFPJcuAcyaNctVVVVFvu7t7XX5+fmupqbGcKr4euKJJ1xRUZH1GMNGktu6dWvk676+Puf3+90zzzwTuayjo8N5vV73+uuvG0wYH18/Tuecq6ysdBUVFSbzDJfjx487Sa6+vt45959zl5qa6rZs2RJZ8/e//91Jcg0NDVZjDtnXj9M557773e+6n/zkJ3ZDDZMrr7zS/eY3v7mo5/KSfwR0+vRpNTY2qrS0NHJZUlKSSktL1dDQYDhZ/B08eFD5+fmaPHmy7r33Xh0+fNh6pGHT0tKitra2qPPq8/lUUlIy4s6rJNXV1SknJ0fTpk3TihUrdOLECeuRhiQYDEqSsrKyJEmNjY3q6emJOp/Tp0/XxIkTE/p8fv04v/Laa68pOztbM2bM0Jo1a3Tq1CmL8eKit7dXmzdvVldXlwKBwEU9l5fch5F+3RdffKHe3l7l5uZGXZ6bm6vPPvvMaKr4Kykp0caNGzVt2jQdO3ZMTz31lG655RYdOHBA6enp1uPFXVtbmySd87x+dd1IMX/+fC1cuFCFhYU6dOiQfv7zn6u8vFwNDQ1KTk62Hi9mfX19WrVqlW666SbNmDFD0n/OZ1pamjIzM6PWJvL5PNdxStI999yjSZMmKT8/X/v379fDDz+spqYmvfnmm4bTxu7jjz9WIBBQd3e3xo4dq61bt+raa6/Vvn37Ltq5vOQDdLkoLy+P/HnmzJkqKSnRpEmT9Pvf/15Lly41nAxDddddd0X+fP3112vmzJmaMmWK6urqNHfuXMPJBqeqqkoHDhxI+OcoL6S/47z//vsjf77++uuVl5enuXPn6tChQ5oyZcrFHnPQpk2bpn379ikYDOoPf/iDKisrVV9ff1FnuOR/BZedna3k5OSzXoHR3t4uv99vNNXwy8zM1NVXX63m5mbrUYbFV+fucjuvkjR58mRlZ2cn5LlduXKl3n77bX3wwQdR/9sUv9+v06dPq6OjI2p9op7P/o7zXEpKSiQp4c5nWlqapk6dquLiYtXU1KioqEgvvPDCRT2Xl3yA0tLSVFxcrNra2shlfX19qq2tVSAQMJxseJ08eVKHDh1SXl6e9SjDorCwUH6/P+q8hkIh7dmzZ0SfV0k6cuSITpw4kVDn1jmnlStXauvWrdq5c6cKCwujri8uLlZqamrU+WxqatLhw4cT6nxe6DjPZd++fZKUUOfzXPr6+hQOhy/uuYzrSxqGyebNm53X63UbN250n376qbv//vtdZmama2trsx4tbn7605+6uro619LS4v70pz+50tJSl52d7Y4fP2492qB1dna6jz76yH300UdOknv22WfdRx995P71r38555x7+umnXWZmptu+fbvbv3+/q6iocIWFhe7LL780njw25zvOzs5O9+CDD7qGhgbX0tLi3n//ffftb3/bXXXVVa67u9t69AFbsWKF8/l8rq6uzh07diyynTp1KrJm+fLlbuLEiW7nzp1u7969LhAIuEAgYDh17C50nM3Nze4Xv/iF27t3r2tpaXHbt293kydPdrNnzzaePDaPPPKIq6+vdy0tLW7//v3ukUcecR6Px/3xj390zl28c5kQAXLOuZdeeslNnDjRpaWluVmzZrndu3dbjxRXixcvdnl5eS4tLc194xvfcIsXL3bNzc3WYw3JBx984CSdtVVWVjrn/vNS7Mcee8zl5uY6r9fr5s6d65qammyHHoTzHeepU6fcvHnz3Pjx411qaqqbNGmSW7ZsWcL9x9O5jk+Se+WVVyJrvvzyS/fjH//YXXnllW7MmDHujjvucMeOHbMbehAudJyHDx92s2fPdllZWc7r9bqpU6e6n/3sZy4YDNoOHqMf/ehHbtKkSS4tLc2NHz/ezZ07NxIf5y7eueR/xwAAMHHJPwcEABiZCBAAwAQBAgCYIEAAABMECABgggABAEwQIACACQIEADBBgAAAJggQAMAEAQIAmCBAAAAT/x9jSFla1Y648gAAAABJRU5ErkJggg==", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# create a tensor train, which has ones at one side else zeros\n", "with ts.variational(max_rank=1):\n", " data = np.zeros(shape.dims[0].size())\n", " data[0] = 1.0\n", " rhs = ts.tensortrain(xshape, data)\n", " rhs.extend(ts.full(yshape, 1))\n", "\n", "# multiply the train along the y-axis with a gauss-shapes function\n", "with ts.variational(max_rank=10):\n", " y = np.linspace(domains[1].lower, domains[1].upper, shape.dims[1].size())\n", " data = np.exp(-10*(y-0.5)**2) # change this part to get another boundary condition\n", " #data = (y-0.5)**2\n", " tmp = ts.tensortrain(yshape, data)\n", " rhs = ts.einsum(\"ab,b->ab\", rhs, tmp)\n", "\n", "# plot the resulting matrix\n", "plt.figure()\n", "plt.imshow(rhs.to_tensor())\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "bdd5649f-1bc0-45ea-94ae-3672da171ef4", "metadata": {}, "source": [ "In the next cell we define the boundary condition $b$. It is chosen so one site\n", "of the square represents a gaussian shaped function." ] }, { "cell_type": "code", "execution_count": 4, "id": "21899e26-7cc2-4b68-b24c-f76ff7d2928a", "metadata": {}, "outputs": [], "source": [ "# define the laplace operators\n", "shape.ranks = 15\n", "laplace_ops = [laplace(grid, i) for i in range(len(shape.dims))]\n", "laplace_maps = [ts.linear_map(\"imjn,mn->ij\", op, shape) for op in laplace_ops]" ] }, { "cell_type": "markdown", "id": "65a72be9-3158-4615-98e5-62a2a15c157a", "metadata": {}, "source": [ "After defining the right hand side of the equation we will see how to define the linear\n", "operators used to define the linear equation system. Since $\\Delta$ is the sum of its $x$\n", "and $y$ component, we can define it as two seperate operators $\\Delta_x$ and $\\Delta_y$." ] }, { "cell_type": "code", "execution_count": 5, "id": "7343afe8-953b-46d2-aa81-f7d20352655c", "metadata": {}, "outputs": [], "source": [ "# define the strategy and local solvers\n", "strat = ts.sweeping_strategy(ncores=2, nsweeps=5)\n", "decomp = ts.svdecomposition(max_rank=10, cutoff=1e-6)\n", "solver = ts.gmres(nsteps=25, subspace=25, eps=1e-6)\n", "\n", "# define the tensorized linear solver\n", "lin_solver = ts.linsolver(\n", " rhs,\n", " *laplace_maps,\n", " method=\"dmrg\",\n", " strategy=strat,\n", " decomposition=decomp,\n", " solver=solver)" ] }, { "cell_type": "markdown", "id": "e3ac63dc-ae16-4cd7-add1-7398c4fff75e", "metadata": {}, "source": [ "At this step all components of the equation are defined and we can instantiate the solver.\n", "For doing so we need a sweeping strategy, a decomposition (for strategies with ncores >= 2)\n", "and a non tensorized solver for linear equation systems. Since it is very convenient to\n", "implicitly define the local linear maps the library implements a GMRES solver." ] }, { "cell_type": "code", "execution_count": 6, "id": "a8c27a32-0500-45a9-a237-3e6e2677b215", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "7.08044032E-07" ] } ], "source": [ "# callback for printing and retrieval of local results\n", "vals = []\n", "def callback(lrange, res):\n", " print(f\"{res.residuals[-1]:.8E}\", end=\"\\r\", flush=True)\n", " vals.append(res.residuals[-1])\n", " return False\n", "\n", "# create some start guess and solve the system\n", "guess = ts.full(shape, 1)\n", "res = lin_solver(guess, callback=callback)" ] }, { "cell_type": "markdown", "id": "489d037b-2445-404f-824d-892e4aa3d6d5", "metadata": {}, "source": [ "With the solver we can now start the process of solving the equation." ] }, { "cell_type": "code", "execution_count": 7, "id": "311fe2e7-3f51-429a-ac14-1070ee3cc8df", "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaAAAAGdCAYAAABU0qcqAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjgsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvwVt1zgAAAAlwSFlzAAAPYQAAD2EBqD+naQAAKMFJREFUeJzt3X9wVOd97/HP2ZW0/JJWFiAJFUHAJGAbQ6fUVjROKDEqPzrjgYDn2knmFqcee0yFpzZNk9BJ7LjtjFx7xnGSIfiPTk1yJ5jEHQNjz7VdW0Ri3ApaqBliJ9YYqhY8IGHTshICraTd5/7B9aaykX2+0h6eXfF+zewMaB8dPWfPnvPR6pz9bOCccwIA4CqL+Z4AAODaRAABALwggAAAXhBAAAAvCCAAgBcEEADACwIIAOAFAQQA8KLE9wQ+KpvN6vTp0yovL1cQBL6nAwAwcs6pr69PdXV1isVGf51TcAF0+vRp1dfX+54GAGCcTp06pdmzZ496f2QBtH37dj355JPq7u7W0qVL9aMf/Ui33nrrp35feXm5JGnFzE0qiZWF+lmusjz0vIaTk0OPlaTBZLg5SNJgZdy07IHrwv8FNH2dadEarMqGH1yVNi27enqvafyC5Aehx9407Yxp2Tck3gs9dn7pedOya+Lht/2UkM/VsbqYHQw9ticTfqwk/ftQZeixv0mPfjC5krcvzAo99nhqhmnZZ89VhB/8XwnTssv+y3Z2IvHf4cdO+m/Dvimp7Hwm/NiUbduXpC6FHhuc7ws9djg7qLb3f5I7no/680Mv0eDnP/+5tm7dqmeeeUYNDQ16+umntXr1anV2dqq6uvoTv/fDP7uVxMrCB1Dc8OQqmRR+rKRsafgDS6bUFkDxsvBPcssqSlJskuFJPsX2p86SqbbAKpsW/jGcNM32lJw6KfxjPq3UdlCpiIcfP+UT/syQDyXZ8Mvvz9jmMnUo/GM4qdS2fcoUftuXDNue5LFLhn35om3Z8YTtMTT8rqJ4qS2ASkrDB1BJiW3eJfHwcwlitv1e0qeeRolkr3nqqad033336etf/7puvPFGPfPMM5oyZYr+/u//PoofBwAoQnkPoMHBQR05ckRNTU2//SGxmJqamtTR0fGx8el0Wr29vSNuAICJL+8B9MEHHyiTyaimpmbE12tqatTd3f2x8S0tLUomk7kbFyAAwLXB+/uAtm3bplQqlbudOnXK95QAAFdB3i9CmDFjhuLxuHp6ekZ8vaenR7W1tR8bn0gklEgYz7ADAIpe3l8BlZWVadmyZWptbc19LZvNqrW1VY2Njfn+cQCAIhXJZdhbt27Vpk2b9Pu///u69dZb9fTTT6u/v19f//rXo/hxAIAiFEkA3XXXXXr//ff1yCOPqLu7W7/7u7+rV1555WMXJgAArl2RNSFs2bJFW7ZsGfP3u+FhubBv7ssY3kxlGCtJQdYZxpoWrSBjWPaw7c2isXT48cMXbU+D/+qbahrfFQv/wGScbT1TmfDNFv9ubEKYWRL+nd9Tx/AmPYv+bPh3/b8/HL4ZRJLOGJoQui5ONy37ZF9V6LHW51XW8LwtMewPkhQMm4bb9mXrccJyDDIe3yzHTjcc/kFx2XBjvV8FBwC4NhFAAAAvCCAAgBcEEADACwIIAOAFAQQA8IIAAgB4QQABALwggAAAXhBAAAAvIqviGTfnLt/CyBqqeMK3Wlweb6jYiBnGSlJsKHw9SHzQtGjFB8IvO3sxblr2QGmZafzZYFrosUNZ2+9EvYPhq3jeS1xnWnZ56UDosYmYsbvFKJ0Nv6v2DU0yLftcOnwFzvuXbHU5/31hSuixAxdsz6uY4Xlr2R8k+/4WGzKMNR4nTDU/xuOb5dgZ+nhsGMsrIACAFwQQAMALAggA4AUBBADwggACAHhBAAEAvCCAAABeEEAAAC8IIACAFwQQAMALAggA4EVhd8Flw/UJBSHHSbJ1H0kKDMODjGnRihv6o+Jp47LD15gpW2bryXJx29PmkmFsJmP7nejSYGnosR+U2XrMJpeG30ClMePGNxrKhu89uzQU/jGRpIuGx7D/UsK07MGL4ZcdXLA9r+KXDF2Khv1BGsP+ZtiXrccJyzHIfnwLf+x0luMsXXAAgEJGAAEAvCCAAABeEEAAAC8IIACAFwQQAMALAggA4AUBBADwggACAHhBAAEAvCjgKp6spJC1EiFrHyQpyBjqJCTFMuGrLWJDxmUbxscGTYtWPB2+psRdNFbxBLbfWzKGp9ngsG3Zw4PhK2r6y2wdKCUl4cfHYrZtb5XNht9Gw8PhHxPJ9hhmB2yHjGAg/PaM99u2fYnheWut1okNRrgvW48ThmOQ9fhmOXZePibndyyvgAAAXhBAAAAvCCAAgBcEEADACwIIAOAFAQQA8IIAAgB4QQABALwggAAAXhBAAAAvCCAAgBeF2wVnkTV0FA0bxkoKhg0dT4axkhQ3dEJZut0kKWvYss5WHSYX2OYSuPC/52SGjes5FH78YInt6T5YYugYjLgLzhm64GR8DC3jY2nb76zxAUNf2yXbvOMD0YyV7N1xln3ZepywHIOsxzfTsTMCvAICAHiR9wD63ve+pyAIRtwWLVqU7x8DAChykfwJ7qabbtLrr7/+2x9i/NMHAGDiiyQZSkpKVFtbG8WiAQATRCTngN59913V1dVp/vz5+trXvqaTJ0+OOjadTqu3t3fEDQAw8eU9gBoaGrRz50698sor2rFjh7q6uvTFL35RfX19Vxzf0tKiZDKZu9XX1+d7SgCAAhQ4Z/lMVrvz589r7ty5euqpp3Tvvfd+7P50Oq10+rfXPPb29qq+vl4rK/+3SoKyUD8jSFaEnk+2fGrosZKUSU4KPXaovNS07KFp4fN/cKrtd4WhqYaPcLY9JBqebBufmRz+KZaZZHs6ZhOGy0gNl1Vbx3MZ9pVFeRl2ySXD2H7TolXab9ueZf3hn4elF2yXPpf2DYUeG0/ZrjeP9YV/YFwq/F+nht2gWs//H6VSKVVUjH58jvzqgMrKSn3uc5/T8ePHr3h/IpFQIpGIehoAgAIT+fuALly4oBMnTmjWrFlR/ygAQBHJewB94xvfUHt7u/7jP/5D//zP/6wvf/nLisfj+spXvpLvHwUAKGJ5/xPce++9p6985Ss6d+6cZs6cqS984Qs6ePCgZs6caVtQ1klByL/DZsL/TTUwVk8EQ5nQY2NDtk6bmKFGJj5oPDdiquIxnjMwCgznL4LwD7ckKTsY/jF3xnNAlooiF/E5oFiEj2FgOQcU/nSEJOM5IGv9jaWKJ22syTLubzFLFc9QdMcg6/HNcuxU1vCYhLy0IO8BtHv37nwvEgAwAdEFBwDwggACAHhBAAEAvCCAAABeEEAAAC8IIACAFwQQAMALAggA4AUBBADwggACAHgR+ccxjJVzTk7h+oQCy0caWbqPJAWG8bFh27Jjg+HzPx5hj5n115DA2brjAsPDEmRsy84aPoLJGZ/tpi64IOI+PcPmL6guuEHDWGsXnKHfrSTqLjjDeOtxwnIMsh7fwna2XR6a/7G8AgIAeEEAAQC8IIAAAF4QQAAALwggAIAXBBAAwAsCCADgBQEEAPCCAAIAeEEAAQC8KNgqHmWz4TtcMobuEctYSYGhNiMYstVgxAfDj3cltqoXFw8/3tisY6rvkCRlw/8Aa41M1lANY34MLcOjbeJRyFYqSbbaHkkKhsOPjRnGSlLMUsVjrL+xVPdYKoHGNBfDvmw9TliOQdbjm2l81jAPF24sr4AAAF4QQAAALwggAIAXBBAAwAsCCADgBQEEAPCCAAIAeEEAAQC8IIAAAF4QQAAALwggAIAXhdsFZ+AM3WRBxtbDJEMPU8zY8eRKDcsetJWNxeOW0bbfQ8xdY4aHJWvtgjP0uznTYyI5y8NSSF1wxqe4pX8vNmzb+DFDV5+9f80y1vagxKzjIzxOmPoojcc3y7EzCrwCAgB4QQABALwggAAAXhBAAAAvCCAAgBcEEADACwIIAOAFAQQA8IIAAgB4QQABALwggAAAXhRuF5xzCl2AlTX0GWWtXXCGoixDZ5MkBYZOqFjc9ruCi1k6noyPibMVn1m6yTIZ27Jj8fDrGWUXnPEhMbP070XZBWcZK0nxofATjxnGmpedtnbY2cZb9mXrccJ0DLIe3yzHTktvXMixvAICAHhhDqADBw7ojjvuUF1dnYIg0N69e0fc75zTI488olmzZmny5MlqamrSu+++m6/5AgAmCHMA9ff3a+nSpdq+ffsV73/iiSf0wx/+UM8884wOHTqkqVOnavXq1RoYGBj3ZAEAE4f5HNDatWu1du3aK97nnNPTTz+t73znO1q3bp0k6ac//alqamq0d+9e3X333eObLQBgwsjrOaCuri51d3erqakp97VkMqmGhgZ1dHRc8XvS6bR6e3tH3AAAE19eA6i7u1uSVFNTM+LrNTU1ufs+qqWlRclkMnerr6/P55QAAAXK+1Vw27ZtUyqVyt1OnTrle0oAgKsgrwFUW1srSerp6Rnx9Z6entx9H5VIJFRRUTHiBgCY+PIaQPPmzVNtba1aW1tzX+vt7dWhQ4fU2NiYzx8FAChy5qvgLly4oOPHj+f+39XVpaNHj6qqqkpz5szRQw89pL/5m7/RZz/7Wc2bN0/f/e53VVdXp/Xr1+dz3gCAImcOoMOHD+tLX/pS7v9bt26VJG3atEk7d+7UN7/5TfX39+v+++/X+fPn9YUvfEGvvPKKJk2aZPtB2Wz4TpGsoaoiY62qCD8+GLL1lMRKwr8AdTFjB4qpGsb4QtjWUmKqhrFWvVjqdbJxY1+OYXghVfFYt08sE/4brNsnNmyoyzGMlaTYYPjx8UHbfh8bNO7Lhioe63HCVK9jPr5FVPPjwo0NnLMU/ESvt7dXyWRSt0/6XyoJykJ9TzB1SujlB5Mnm+bjJifCj02Em+9vl10aemwmYSsyyxrGZ8psAZQtsx1tsyXhx1vGSgTQFRFAH2MOoLRtReOG8cGlIdOyg/SgYdlp07LdpUvhx/ZfDD122A1q/8AvlEqlPvG8vver4AAA1yYCCADgBQEEAPCCAAIAeEEAAQC8IIAAAF4QQAAALwggAIAXBBAAwAsCCADghbkLriBlDRUeGVvFRjBsGF9i7CkZCp//MXONTHTdMIGz/d5iqW+x1MJItnqdWMy2bBeLuF8nIoFlf5Ctq8+6fYJhw7KtVTyG/rWYtYpn2Dbe0u8WRHkMMi7bdOyMAK+AAABeEEAAAC8IIACAFwQQAMALAggA4AUBBADwggACAHhBAAEAvCCAAABeEEAAAC8KtorHOckpXE1E4Ay1Gc5YPZE1LDtjq++Qoe7DUvUhSYGhisf8W4jxIQyy4eeSzdjqb4J4+LHOMPbywsOvqIuw+kiSAsvz1rp9DE+twFjFY6nuCSKs4gkirNaRZNqXzccJyzHIenwzHDudYdlhh/IKCADgBQEEAPCCAAIAeEEAAQC8IIAAAF4QQAAALwggAIAXBBAAwAsCCADgBQEEAPCCAAIAeFGwXXAmWUNHkbGHKTCMD4aN/VFxQ/4P2RYdM3STWTqe/v832IZnw69nELd1qjnDeBczLts03PgYRshQYXd5vGH/sXbBWcZb+9pili7FiLvggkz48ebjhOEYZD2+WY6dUeAVEADACwIIAOAFAQQA8IIAAgB4QQABALwggAAAXhBAAAAvCCAAgBcEEADACwIIAOBF4VbxZF3oThFLlUzgjFUVhooNZeO2ZRsqOQJDtc7lbzDM21itEzPWdzhLVZKlnki2Kh4ZH0JnecyNyzYzPOSBuVrJsGxzFU+EdTlRLts43rIvKxvhMch4fDPVcFn2+5DL5RUQAMALAggA4IU5gA4cOKA77rhDdXV1CoJAe/fuHXH/PffcoyAIRtzWrFmTr/kCACYIcwD19/dr6dKl2r59+6hj1qxZozNnzuRuzz333LgmCQCYeMwXIaxdu1Zr1679xDGJREK1tbVjnhQAYOKL5BxQW1ubqqurtXDhQm3evFnnzp0bdWw6nVZvb++IGwBg4st7AK1Zs0Y//elP1draqr/9279Ve3u71q5dq8wolxK2tLQomUzmbvX19fmeEgCgAOX9fUB333137t8333yzlixZouuvv15tbW1auXLlx8Zv27ZNW7duzf2/t7eXEAKAa0Dkl2HPnz9fM2bM0PHjx694fyKRUEVFxYgbAGDiizyA3nvvPZ07d06zZs2K+kcBAIqI+U9wFy5cGPFqpqurS0ePHlVVVZWqqqr02GOPaePGjaqtrdWJEyf0zW9+UwsWLNDq1avzOnEAQHEzB9Dhw4f1pS99Kff/D8/fbNq0STt27NCxY8f0k5/8ROfPn1ddXZ1WrVqlv/7rv1YikcjfrD/K0q1k6I+SJBe39E0ZOptk63cz9ZIp4moyY5VVyEo/SbbeOElyGUtfm+1RcVH3u0XE8nhLMnUBBsbtY+qOs/a1Gfb7YMi2b5q63WTc943LdhEe38y9dHlmDqAVK1Z8YoHdq6++Oq4JAQCuDXTBAQC8IIAAAF4QQAAALwggAIAXBBAAwAsCCADgBQEEAPCCAAIAeEEAAQC8IIAAAF7k/fOA8sZlFbp0zNBldXm5Bpbuq1E+dG9U8fD5b+2ZszR2mSvPSoxlY5btkzV23ln63aLsgotFXBxneB5G2QVnGitjF5yxlyywdMdF2e0mGfsorcuO8PgW1bEz5FheAQEAvCCAAABeEEAAAC8IIACAFwQQAMALAggA4AUBBADwggACAHhBAAEAvCCAAABeFG4Vj4WlBsNSayFJ2fC1Gc5aI5Mx1mZYlm0Ya21uCYx1LIpFV8VjqtcxVvGYxxcI8/Yp0ioeGfYfc7WOtS7HMBdnXU/DMch+fIvuGBQGr4AAAF4QQAAALwggAIAXBBAAwAsCCADgBQEEAPCCAAIAeEEAAQC8IIAAAF4QQAAALwggAIAXBdsF57JOLgjXa2TqPbN2WVm6law9TMPDhsHGTWXoMQss05D9MbR0wQV0wY1bIXXBmfYJYy+ZqUsxwp45SbZ92dzXFn68ed80cBHMg1dAAAAvCCAAgBcEEADACwIIAOAFAQQA8IIAAgB4QQABALwggAAAXhBAAAAvCCAAgBcFW8VjYamfCMyVHJnwY2O26hZLaUYQGOZhFTfWE1nrPoLwj7mLG38nirCKpziLeMYgwioeU12OuebHsGxrtY5lv5fkTLVAxn3ZMhfj8S3K6p4weAUEAPDCFEAtLS265ZZbVF5erurqaq1fv16dnZ0jxgwMDKi5uVnTp0/XtGnTtHHjRvX09OR10gCA4mcKoPb2djU3N+vgwYN67bXXNDQ0pFWrVqm/vz835uGHH9aLL76o559/Xu3t7Tp9+rQ2bNiQ94kDAIpb4MbxR8D3339f1dXVam9v1/Lly5VKpTRz5kzt2rVLd955pyTpnXfe0Q033KCOjg59/vOf/9Rl9vb2KplMakVsg0qC0nArURr+VFZQYjvtFcTj4Qcb5iFJioVfdlBimIckWeZtPe8Si+48TSGdA7pmcA7oCuON54CGLedpjOeAhsJ/1IMzz9uwbMM8ht2Q2rIvKJVKqaKiYtRx4zoHlEqlJElVVVWSpCNHjmhoaEhNTU25MYsWLdKcOXPU0dFxxWWk02n19vaOuAEAJr4xB1A2m9VDDz2k2267TYsXL5YkdXd3q6ysTJWVlSPG1tTUqLu7+4rLaWlpUTKZzN3q6+vHOiUAQBEZcwA1Nzfrrbfe0u7du8c1gW3btimVSuVup06dGtfyAADFYUzvA9qyZYteeuklHThwQLNnz859vba2VoODgzp//vyIV0E9PT2qra294rISiYQSicRYpgEAKGKmV0DOOW3ZskV79uzR/v37NW/evBH3L1u2TKWlpWptbc19rbOzUydPnlRjY2N+ZgwAmBBMr4Cam5u1a9cu7du3T+Xl5bnzOslkUpMnT1YymdS9996rrVu3qqqqShUVFXrwwQfV2NgY6go4AMC1wxRAO3bskCStWLFixNefffZZ3XPPPZKk73//+4rFYtq4caPS6bRWr16tH//4x3mZLABg4hjX+4CiMKb3ARne82J5z5Ak0/tpAut7WCzvSbK8r0dSYHmvjvV9PebxhvffRPgeIytn7PYrFEE2wl06yvfqWOdtWraxIy3KzkjDe28kY8+c9X1AEb3H6Kq8DwgAgLEigAAAXhBAAAAvCCAAgBcEEADACwIIAOAFAQQA8IIAAgB4QQABALwggAAAXozp4xiKmrViI8qPfLZ+TLCBpdQkiLqNyRkeF2sdS4R1OUG2OKt4zHU5FtbtE+HHfVtqZ8yNY8ZKG9O+XECVQ77xCggA4AUBBADwggACAHhBAAEAvCCAAABeEEAAAC8IIACAFwQQAMALAggA4AUBBADwggACAHhRuF1wLispZK+RoWvMWXrJZOtJs/ZNBda+qYi4wPh7iHU9o+zTi0X4O5R1LoUi0i44Y9dYhPuPqVPNGedt7WmMspcuysfQ8rhEMJZXQAAALwggAIAXBBAAwAsCCADgBQEEAPCCAAIAeEEAAQC8IIAAAF4QQAAALwggAIAXhVvFExVLfYckBcZKDgMXj4efRpS1PTFrfYetosZU9RMz1t8YqmFMlUDXEHN9i0WUdTmWeVv3e+P+ZnoMrfuypf7Iup6e8QoIAOAFAQQA8IIAAgB4QQABALwggAAAXhBAAAAvCCAAgBcEEADACwIIAOAFAQQA8IIAAgB4MSG64Jyh/ygIX792edmGjqfA0tkkSYZuMhez/a5g6o4zdrvJ0u0m2brmrNV7lsfQuGjzehYKa6eaadkF0hsn2dbTuGxzP56pr822faLs6rMcO6NQpHsYAKDYmQKopaVFt9xyi8rLy1VdXa3169ers7NzxJgVK1YoCIIRtwceeCCvkwYAFD9TALW3t6u5uVkHDx7Ua6+9pqGhIa1atUr9/f0jxt133306c+ZM7vbEE0/kddIAgOJnOgf0yiuvjPj/zp07VV1drSNHjmj58uW5r0+ZMkW1tbX5mSEAYEIa1zmgVColSaqqqhrx9Z/97GeaMWOGFi9erG3btunixYujLiOdTqu3t3fEDQAw8Y35KrhsNquHHnpIt912mxYvXpz7+le/+lXNnTtXdXV1OnbsmL71rW+ps7NTL7zwwhWX09LSoscee2ys0wAAFKnAjfEav82bN+vll1/WG2+8odmzZ486bv/+/Vq5cqWOHz+u66+//mP3p9NppdPp3P97e3tVX1+vFcF6lQSl4SZjuFw2MHwMtiTTR0SbP/LZMhfrZdiWuVg/Btt8GXaEH4Ud5cdscxn2FZbNZdhXXr5hLlF+3Ld1PU1v1wi/jsNuSG1ur1KplCoqKkYdN6ZXQFu2bNFLL72kAwcOfGL4SFJDQ4MkjRpAiURCiURiLNMAABQxUwA55/Tggw9qz549amtr07x58z71e44ePSpJmjVr1pgmCACYmEwB1NzcrF27dmnfvn0qLy9Xd3e3JCmZTGry5Mk6ceKEdu3apT/6oz/S9OnTdezYMT388MNavny5lixZEskKAACKkymAduzYIenym03/p2effVb33HOPysrK9Prrr+vpp59Wf3+/6uvrtXHjRn3nO9/J24QBABOD+U9wn6S+vl7t7e3jmlDkrCdos+FPRDvjOWtbX5vx5KLl5Lxx4kFgPEFrOeca9QURJhGezC9W5v0nwh4z08l547wjvAjBfoGDYXyUF6BEoEgv8wEAFDsCCADgBQEEAPCCAAIAeEEAAQC8IIAAAF4QQAAALwggAIAXBBAAwAsCCADgxZg/kC5yzkkKW0FhqMEwVOtIUhCzfOaIbdmWBpzAWiVi+fwg6+eTRPkZPNY+I8O2N39e0zVijB8JFo71eWtRQPOO8jN7LPU6LsrPVLKsY8ixvAICAHhBAAEAvCCAAABeEEAAAC8IIACAFwQQAMALAggA4AUBBADwggACAHhBAAEAvCCAAABeFG4XXIGwdCuZeuMkU3ecC92L9+E3GOYddUdahL10FubmsGLtjouyIy1KEfbGmfvurJ1qFpb+NY2h362I8AoIAOAFAQQA8IIAAgB4QQABALwggAAAXhBAAAAvCCAAgBcEEADACwIIAOAFAQQA8GJiVPGYajaMdR+BoS7HWJlhru6xiLLmx6pYq2FwdRVS5YyxLse26MKp+fG9b/IKCADgBQEEAPCCAAIAeEEAAQC8IIAAAF4QQAAALwggAIAXBBAAwAsCCADgBQEEAPCCAAIAeDExuuAszN1HUXZCRZn/mciWHMQC2zdENxVcwyLtVItShD1zvrvdrHgFBADwwhRAO3bs0JIlS1RRUaGKigo1Njbq5Zdfzt0/MDCg5uZmTZ8+XdOmTdPGjRvV09OT90kDAIqfKYBmz56txx9/XEeOHNHhw4d1++23a926dXr77bclSQ8//LBefPFFPf/882pvb9fp06e1YcOGSCYOAChugXPj+6NhVVWVnnzySd15552aOXOmdu3apTvvvFOS9M477+iGG25QR0eHPv/5z4daXm9vr5LJpFZonUqC0vFMLT8C4/kO07KL8y+g5nNAQAQ4B3SlZRfGYzLshtSmfUqlUqqoqBh13JiPgJlMRrt371Z/f78aGxt15MgRDQ0NqampKTdm0aJFmjNnjjo6OkZdTjqdVm9v74gbAGDiMwfQr371K02bNk2JREIPPPCA9uzZoxtvvFHd3d0qKytTZWXliPE1NTXq7u4edXktLS1KJpO5W319vXklAADFxxxACxcu1NGjR3Xo0CFt3rxZmzZt0q9//esxT2Dbtm1KpVK526lTp8a8LABA8TC/D6isrEwLFiyQJC1btkz/+q//qh/84Ae66667NDg4qPPnz494FdTT06Pa2tpRl5dIJJRIJOwzBwAUtXGfBc9ms0qn01q2bJlKS0vV2tqau6+zs1MnT55UY2PjeH8MAGCCMb0C2rZtm9auXas5c+aor69Pu3btUltbm1599VUlk0nde++92rp1q6qqqlRRUaEHH3xQjY2Noa+AAwBcO0wBdPbsWf3xH/+xzpw5o2QyqSVLlujVV1/VH/7hH0qSvv/97ysWi2njxo1Kp9NavXq1fvzjH0cy8asmyssanaGjJsrLwY0s0wauCQVy+XOxGff7gPKt4N4HVCgKKIAAfERhHUa9i/x9QAAAjAcBBADwggACAHhBAAEAvCCAAABeEEAAAC8IIACAFwQQAMALAggA4IW5DTtqHxYzDGtI4s3F/wNNCEDBoglhhGENSfrt8Xw0BRdAfX19kqQ39H89z6TA8PwGUGT6+vqUTCZHvb/guuCy2axOnz6t8vJyBf+j/6y3t1f19fU6derUJ3YLFTvWc+K4FtZRYj0nmnysp3NOfX19qqurUyw2+pmegnsFFIvFNHv27FHvr6iomNAb/0Os58RxLayjxHpONONdz0965fMhLkIAAHhBAAEAvCiaAEokEnr00UeVSCR8TyVSrOfEcS2so8R6TjRXcz0L7iIEAMC1oWheAQEAJhYCCADgBQEEAPCCAAIAeFE0AbR9+3Z95jOf0aRJk9TQ0KB/+Zd/8T2lvPre976nIAhG3BYtWuR7WuNy4MAB3XHHHaqrq1MQBNq7d++I+51zeuSRRzRr1ixNnjxZTU1Nevfdd/1Mdhw+bT3vueeej23bNWvW+JnsGLW0tOiWW25ReXm5qqurtX79enV2do4YMzAwoObmZk2fPl3Tpk3Txo0b1dPT42nGYxNmPVesWPGx7fnAAw94mvHY7NixQ0uWLMm92bSxsVEvv/xy7v6rtS2LIoB+/vOfa+vWrXr00Uf1b//2b1q6dKlWr16ts2fP+p5aXt100006c+ZM7vbGG2/4ntK49Pf3a+nSpdq+ffsV73/iiSf0wx/+UM8884wOHTqkqVOnavXq1RoYGLjKMx2fT1tPSVqzZs2Ibfvcc89dxRmOX3t7u5qbm3Xw4EG99tprGhoa0qpVq9Tf358b8/DDD+vFF1/U888/r/b2dp0+fVobNmzwOGu7MOspSffdd9+I7fnEE094mvHYzJ49W48//riOHDmiw4cP6/bbb9e6dev09ttvS7qK29IVgVtvvdU1Nzfn/p/JZFxdXZ1raWnxOKv8evTRR93SpUt9TyMyktyePXty/89ms662ttY9+eSTua+dP3/eJRIJ99xzz3mYYX58dD2dc27Tpk1u3bp1XuYTlbNnzzpJrr293Tl3eduVlpa6559/PjfmN7/5jZPkOjo6fE1z3D66ns459wd/8Afuz/7sz/xNKiLXXXed+7u/+7urui0L/hXQ4OCgjhw5oqamptzXYrGYmpqa1NHR4XFm+ffuu++qrq5O8+fP19e+9jWdPHnS95Qi09XVpe7u7hHbNZlMqqGhYcJtV0lqa2tTdXW1Fi5cqM2bN+vcuXO+pzQuqVRKklRVVSVJOnLkiIaGhkZsz0WLFmnOnDlFvT0/up4f+tnPfqYZM2Zo8eLF2rZtmy5evOhjenmRyWS0e/du9ff3q7Gx8apuy4IrI/2oDz74QJlMRjU1NSO+XlNTo3feecfTrPKvoaFBO3fu1MKFC3XmzBk99thj+uIXv6i33npL5eXlvqeXd93d3ZJ0xe364X0TxZo1a7RhwwbNmzdPJ06c0F/+5V9q7dq16ujoUDwe9z09s2w2q4ceeki33XabFi9eLOny9iwrK1NlZeWIscW8Pa+0npL01a9+VXPnzlVdXZ2OHTumb33rW+rs7NQLL7zgcbZ2v/rVr9TY2KiBgQFNmzZNe/bs0Y033qijR49etW1Z8AF0rVi7dm3u30uWLFFDQ4Pmzp2rX/ziF7r33ns9zgzjdffdd+f+ffPNN2vJkiW6/vrr1dbWppUrV3qc2dg0NzfrrbfeKvpzlJ9mtPW8//77c/+++eabNWvWLK1cuVInTpzQ9ddff7WnOWYLFy7U0aNHlUql9A//8A/atGmT2tvbr+ocCv5PcDNmzFA8Hv/YFRg9PT2qra31NKvoVVZW6nOf+5yOHz/ueyqR+HDbXWvbVZLmz5+vGTNmFOW23bJli1566SX98pe/HPGxKbW1tRocHNT58+dHjC/W7Tnael5JQ0ODJBXd9iwrK9OCBQu0bNkytbS0aOnSpfrBD35wVbdlwQdQWVmZli1bptbW1tzXstmsWltb1djY6HFm0bpw4YJOnDihWbNm+Z5KJObNm6fa2toR27W3t1eHDh2a0NtVkt577z2dO3euqLatc05btmzRnj17tH//fs2bN2/E/cuWLVNpaemI7dnZ2amTJ08W1fb8tPW8kqNHj0pSUW3PK8lms0qn01d3W+b1koaI7N692yUSCbdz507361//2t1///2usrLSdXd3+55a3vz5n/+5a2trc11dXe6f/umfXFNTk5sxY4Y7e/as76mNWV9fn3vzzTfdm2++6SS5p556yr355pvuP//zP51zzj3++OOusrLS7du3zx07dsytW7fOzZs3z126dMnzzG0+aT37+vrcN77xDdfR0eG6urrc66+/7n7v937Pffazn3UDAwO+px7a5s2bXTKZdG1tbe7MmTO528WLF3NjHnjgATdnzhy3f/9+d/jwYdfY2OgaGxs9ztru09bz+PHj7q/+6q/c4cOHXVdXl9u3b5+bP3++W758ueeZ23z729927e3trquryx07dsx9+9vfdkEQuH/8x390zl29bVkUAeSccz/60Y/cnDlzXFlZmbv11lvdwYMHfU8pr+666y43a9YsV1ZW5n7nd37H3XXXXe748eO+pzUuv/zlL52kj902bdrknLt8KfZ3v/tdV1NT4xKJhFu5cqXr7Oz0O+kx+KT1vHjxolu1apWbOXOmKy0tdXPnznX33Xdf0f3ydKX1k+SeffbZ3JhLly65P/3TP3XXXXedmzJlivvyl7/szpw542/SY/Bp63ny5Em3fPlyV1VV5RKJhFuwYIH7i7/4C5dKpfxO3OhP/uRP3Ny5c11ZWZmbOXOmW7lyZS58nLt625KPYwAAeFHw54AAABMTAQQA8IIAAgB4QQABALwggAAAXhBAAAAvCCAAgBcEEADACwIIAOAFAQQA8IIAAgB4QQABALz4f6WB0szg9ZrxAAAAAElFTkSuQmCC", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# plot the heat map\n", "plt.figure()\n", "plt.imshow(res.to_tensor())\n", "plt.show()" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "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.14.2" } }, "nbformat": 4, "nbformat_minor": 5 }