s1MG(s).ipynb

{"nbformat":4,"nbformat_minor":0,"metadata":{"colab":{"name":"s1MG(s).ipynb","provenance":[],"collapsed_sections":[],"authorship_tag":"ABX9TyO7Xgr7F7EigKAOIfwJyyIy"},"kernelspec":{"name":"python3","display_name":"Python 3"},"accelerator":"GPU"},"cells":[{"cell_type":"code","metadata":{"id":"FKdbL4Dx_YHJ","executionInfo":{"status":"ok","timestamp":1612978491713,"user_tz":-180,"elapsed":2307,"user":{"displayName":"Vl F","photoUrl":"","userId":"14485856366916106708"}}},"source":["import tensorflow as tf\n","import numpy as np\n","import matplotlib.pyplot as plt\n","import scipy"],"execution_count":1,"outputs":[]},{"cell_type":"markdown","metadata":{"id":"DNZugWEQ_66B"},"source":["# The model\n","\n","+ Restriction operator is fixed.\n","+ Petrov-Galerkin coarse-grid correction.\n","+ Both smoothing operators are adjastable $\\text{conv}D(A)^{-1}\\left(b_k - A_k x_k\\right)$.\n","+ Serialization of a single layer."]},{"cell_type":"code","metadata":{"id":"VWiS_5ht_7zT","executionInfo":{"status":"ok","timestamp":1612978606155,"user_tz":-180,"elapsed":731,"user":{"displayName":"Vl F","photoUrl":"","userId":"14485856366916106708"}}},"source":["# computes A_k*x_k\n","def apply_A(x, A, projection_filters, shapes, strides, k):\n","  if k == 0:\n","    return A(x)\n","  else:\n","    for i in range(k-1, -1, -1):\n","      x = tf.nn.conv_transpose(x, projection_filters[i], shapes[i], strides[i], padding=\"VALID\")\n","    x = A(x)\n","    for i in range(k):\n","      x = tf.nn.convolution(x, projection_filters[i], strides[i])\n","    return x\n","\n","# find diagonal of an emplicitly defined matrix -- works for constant diagonals only\n","def extract_diagonal(A, projection_filters, shapes, strides, k):\n","  y = np.zeros(shapes[k].numpy())\n","  m = shapes[k].numpy()[1]//2 + 1\n","  y[:, m, m, :] = 1\n","  y = tf.Variable(y, dtype=projection_filters[0].dtype)\n","  D = apply_A(y, A, projection_filters, shapes, strides, k)[0, m, m, 0]\n","  return D\n","\n","# slightly generalized Jacobi smoother\n","def smooth_r(x, b, A, D, projection_filters, shapes, strides, smoothing_filter, k):\n","  r = b - apply_A(x, A, projection_filters, shapes, strides, k)\n","  x = x + tf.nn.convolution(r/D, smoothing_filter, padding=\"SAME\")\n","  return x\n","\n","# for some reason tf.nn.conv_transpose requires shape of an output -- this function compute required shapes\n","def compute_shapes(J, projection_filters, strides, batch_size=10):\n","  shapes = []\n","  x = tf.Variable(np.random.randn(batch_size, 2**J-1, 2**J-1, 1), dtype=projection_filters[0].dtype)\n","  shapes.append(tf.shape(x))\n","  for filter, stride in zip(projection_filters, strides):\n","    x = tf.nn.convolution(x, filter, stride)\n","    shapes.append(tf.shape(x))\n","  return shapes\n","\n","# multigrid\n","def multigrid_r(x, b, A, projection_filters, shapes, strides, smoothing_filters_pre, smoothing_filters_post, k=0):\n","  if k == len(projection_filters):\n","    return b/apply_A(tf.ones_like(b), A, projection_filters, shapes, strides, k)\n","  else:\n","    D = extract_diagonal(A, projection_filters, shapes, strides, k)\n","    x = smooth_r(x, b, A, D, projection_filters, shapes, strides, smoothing_filters_pre[k], k)\n","    r = b - apply_A(x, A, projection_filters, shapes, strides, k)\n","    r = tf.nn.convolution(r, projection_filters[k], strides=strides[k])\n","    e = multigrid_r(tf.zeros_like(r), r, A, projection_filters, shapes, strides, smoothing_filters_pre, smoothing_filters_post, k=k+1)\n","    e = tf.nn.conv_transpose(e, projection_filters[k], shapes[k], strides=strides[k], padding=\"VALID\")\n","    x = x + e\n","    x = smooth_r(x, b, A, D, projection_filters, shapes, strides, smoothing_filters_post[k], k)\n","    return x\n","\n","# approximate spectral radius\n","def stochastic_trace_r(A, J, projection_filters, shapes, strides, smoothing_filters_pre, smoothing_filters_post, N_sweeps=10, batch_size=10):\n","  x = tf.Variable(np.random.choice([-1, 1], (batch_size, 2**J-1, 2**J-1, 1)), dtype=projection_filters[0].dtype)\n","  for _ in range(N_sweeps):\n","    x = multigrid_r(x, x*0, A, projection_filters, shapes, strides, smoothing_filters_pre, smoothing_filters_post)\n","  return ((tf.norm(x)**2)/batch_size)**(1/(2*N_sweeps))"],"execution_count":5,"outputs":[]},{"cell_type":"markdown","metadata":{"id":"OWQkMqPSAf03"},"source":["# Equations"]},{"cell_type":"code","metadata":{"id":"1try7DlJAfXZ","executionInfo":{"status":"ok","timestamp":1612978608471,"user_tz":-180,"elapsed":756,"user":{"displayName":"Vl F","photoUrl":"","userId":"14485856366916106708"}}},"source":["# equations\n","\n","poisson_filter = tf.Variable(np.array([[0, -1/4, 0], [-1/4, 1, -1/4], [0, -1/4, 0]]).reshape(3, 3, 1, 1), dtype=tf.float32, trainable=False, name=\"Poisson\")\n","A_poisson = lambda x: tf.nn.convolution(x, poisson_filter, strides=[1, 1], padding=\"SAME\")\n","\n","poisson_large_filter = tf.Variable(np.array([[0, 0, 1/60, 0, 0], [0, 0, -16/60, 0, 0], [1/60, -16/60, 1, -16/60, 1/60], [0, 0, -16/60, 0, 0], [0, 0, 1/60, 0, 0]]).reshape(5, 5, 1, 1), dtype=tf.float32, trainable=False, name=\"Poisson large\")\n","A_poisson_large = lambda x: tf.nn.convolution(x, poisson_large_filter, strides=[1, 1], padding=\"SAME\")\n","\n","mehrstellen_filter = tf.Variable(np.array([[-1/20, -1/5, -1/20], [-1/5, 1, -1/5], [-1/20, -1/5, -1/20]]).reshape(3, 3, 1, 1), dtype=tf.float32, trainable=False, name=\"Mehrstellen\")\n","A_mehrstellen = lambda x: tf.nn.convolution(x, mehrstellen_filter, strides=[1, 1], padding=\"SAME\")\n","\n","epsilon = 10\n","poisson_anisotropic_filter = tf.Variable(np.array([[0, -1/(2+2*epsilon), 0], [-epsilon/(2+2*epsilon), 1, -epsilon/(2+2*epsilon)], [0, -1/(2+2*epsilon), 0]]).reshape(3, 3, 1, 1), dtype=tf.float32, trainable=False, name=\"Anisotropic Poisson\")\n","A_anisotropic = lambda x: tf.nn.convolution(x, poisson_anisotropic_filter, strides=[1, 1], padding=\"SAME\")\n","\n","tau = 3/4\n","mixed_filter = tf.Variable(np.array([[-tau/8, -1/4, tau/8], [-1/4, 1, -1/4], [tau/8, -1/4, -tau/8]]).reshape(3, 3, 1, 1), dtype=tf.float32, trainable=False, name=\"Mixed\")\n","A_mixed = lambda x: tf.nn.convolution(x, mixed_filter, strides=[1, 1], padding=\"SAME\")"],"execution_count":6,"outputs":[]},{"cell_type":"markdown","metadata":{"id":"Nm_rV-IMAyYP"},"source":["# Training"]},{"cell_type":"code","metadata":{"id":"4WRaF3JnAppc","executionInfo":{"status":"ok","timestamp":1612978616314,"user_tz":-180,"elapsed":7290,"user":{"displayName":"Vl F","photoUrl":"","userId":"14485856366916106708"}}},"source":["A = A_poisson # change this line to try different linear equation\n","N_batch = 10\n","\n","J = 5\n","N_levels = J - 1\n","\n","# we are training a single layer -- one variable for presmoothing, one variable for postsmoothing, restriction is fixed\n","proj_var = tf.Variable(np.outer([1/2, 1, 1/2], [[1/2], [1], [1/2]]).reshape((3, 3, 1, 1))/2, dtype=tf.float32, trainable=False, name=\"projection_var\")\n","smoothing_pre_var = tf.Variable(np.outer([1/2, 1, 1/2], [[1/2], [1], [1/2]]).reshape((3, 3, 1, 1)), dtype=tf.float32, trainable=True, name=\"presmoothing_var\")\n","smoothing_post_var = tf.Variable(np.outer([1/2, 1, 1/2], [[1/2], [1], [1/2]]).reshape((3, 3, 1, 1)), dtype=tf.float32, trainable=True, name=\"postsmoothing_var\")\n","\n","projection_filters = [proj_var for i in range(N_levels)]\n","smoothing_filters_pre = [smoothing_pre_var for i in range(N_levels)]\n","smoothing_filters_post = [smoothing_post_var for i in range(N_levels)]\n","variables = [smoothing_pre_var, smoothing_post_var]\n","\n","strides = [[2, 2]]*N_levels\n","shapes = compute_shapes(J, projection_filters, strides, batch_size=N_batch)"],"execution_count":7,"outputs":[]},{"cell_type":"code","metadata":{"id":"XXM3R6sEBGeT","executionInfo":{"status":"ok","timestamp":1612978721185,"user_tz":-180,"elapsed":86697,"user":{"displayName":"Vl F","photoUrl":"","userId":"14485856366916106708"}}},"source":["opt = tf.keras.optimizers.Adam(learning_rate=0.01)\n","N_epoch = 200\n","\n","history = np.zeros(N_epoch)\n","for i in range(N_epoch):\n","  with tf.GradientTape() as tape:\n","    loss = stochastic_trace_r(A, J, projection_filters, shapes, strides, smoothing_filters_pre, smoothing_filters_post)\n","  opt.apply_gradients(zip(tape.gradient(loss, variables), variables))\n","  history[i] = loss.numpy()"],"execution_count":8,"outputs":[]},{"cell_type":"code","metadata":{"colab":{"base_uri":"https://localhost:8080/","height":296},"id":"vZxFjeUJBAC7","executionInfo":{"status":"ok","timestamp":1612978742946,"user_tz":-180,"elapsed":1321,"user":{"displayName":"Vl F","photoUrl":"","userId":"14485856366916106708"}},"outputId":"9ba26840-c009-40a1-f4f4-2d6984bf996a"},"source":["# training history\n","final_loss = stochastic_trace_r(A, J, projection_filters, shapes, strides, smoothing_filters_pre, smoothing_filters_post, N_sweeps=10, batch_size=10).numpy()\n","fig, ax = plt.subplots()\n","ax.set_yscale(\"log\")\n","ax.plot(history)\n","ax.set_xlabel(\"iteration\")\n","ax.set_ylabel(\"loss\")\n","print(\"final loss = \", final_loss);"],"execution_count":9,"outputs":[{"output_type":"stream","text":["final loss =  0.05901901\n"],"name":"stdout"},{"output_type":"display_data","data":{"image/png":"iVBORw0KGgoAAAANSUhEUgAAAYoAAAEGCAYAAAB7DNKzAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3deXjc1X3v8fd3Rvu+WrY25EUyGLxim80GUpYCgZCQDZJmaRaatMkNt+1NSblNl6e3y22ThrZpCGl4SAIJSSFcSEIKhITVsfG+b5It25Jlydr3bebcP2akCMcSsq2Z34zm83oePR79NJr5+mg0H53f+Z1zzDmHiIjIZHxeFyAiIrFNQSEiIlNSUIiIyJQUFCIiMiUFhYiITCnJ6wIioaioyFVVVXldhohI3Ni6dWurc674bF+blUFRVVXFli1bvC5DRCRumNmxyb6mU08iIjIlBYWIiExJQSEiIlNSUIiIyJQUFCIiMiUFhYiITElBISIiU1JQTPDYxmO8frjV6zJERGKKgiJsaDTAYxuP8anvbmbTkTavyxERiRkKirDUJD+PfeoKyvMz+MSjm6lv7fO6JBGRmKCgmKAoK5XvfXItAH/33H6PqxERiQ0KijPMy03nD9+xiBf2NbOhVuMVIiIKirP45Lr5lOWl87WXDntdioiI5xQUZ5GW7OeDayp482g7TV0DXpcjIuIpBcUk7lheCsBPdzZ5XImIiLcUFJOYX5TJ0rJcnt150utSREQ8paCYwruWl7K7sYujulRWRBKYgmIKty+fB8BP1KsQkQSmoJjCvNx01lYV8OzOkzjnvC5HRMQTCoq3cceKUmpbejlwqsfrUkREPBHzQWFmmWb2HTP7lpl9ONrPf9tlc/H7TIPaIpKwPAkKM3vEzFrMbM8Zx28xs4NmVmtm94cP3wU86Zz7NPCuaNdamJXK1QsLeX7vqWg/tYhITPCqR/EocMvEA2bmB74O3AosAe4xsyVAOXAifLdAFGscd11NMUdO93GyU5PvRCTxeBIUzrlXgfYzDq8Fap1zR5xzw8ATwJ1AA6GwgCnqNbN7zWyLmW05ffr0jNZ7zaIiAN7Q2k8ikoBiaYyijN/0HCAUEGXAj4H3mtk3gJ9M9s3OuYedc6udc6uLi4tntLCL52ZTlJXC6woKEUlASV4X8Hacc33A73tZg5lxzaIi3qhtxTmHmXlZjohIVMVSj6IRqJjweXn4WEy4ZlERrb3DHGzWZbIiklhiKSg2A9VmNt/MUoC7gWc9rmnc2DiF9tQWkUTj1eWxPwB+DSw2swYz+6RzbhT4HPA8sB/4kXNurxf1nU1ZXjoLijI1oC0iCceTMQrn3D2THH8OeC7K5UzbuuointzawPBokJSkWOqMiYhEzqx6tzOzO8zs4a6urog8/jWLiugfDrD9eEdEHl9EJBbNqqBwzv3EOXdvbm5uRB7/ygWF+EzzKUQkscyqoIi03PRklpXnaT6FiCQUBcU5WreoiJ0NXXQPjnhdiohIVCgoztG66iICQcfGujavSxERiQoFxTlaWZlHerJf4xQikjAUFOcoNcnP2vkFGqcQkYQxq4Ii0pfHjlm3qIi60300dWnZcRGZ/WZVUET68tgx66q1nIeIJI5ZFRTRsrgktOy4xilEJBEoKM6Dz2dcvbCI12vbcM55XY6ISEQpKM7TukVFtPYOcai51+tSREQiSkFxnq4ZG6fQ6ScRmeUUFOdpbNnx1w/P7P7cIiKxRkFxAdZXF/HrI230DY16XYqISMQoKC7ArUvnMTgS5KUDLV6XIiISMbMqKKI14W7MmqoC5mSn8tOdJ6PyfCIiXphVQRGtCXdj/D7jncvm8fKh0/RoNVkRmaVmVVB44fZlpQyPBvnF/mavSxERiQgFxQVaWZFHUVYqvzqgq59EZHZSUFwgn8+4tqaI1w6fJhDULG0RmX0UFDPguppiOvpH2N0YnUF0EZFoUlDMgPXVxZjBKwd1+klEZh8FxQwoyExhWXkerxzSfAoRmX0UFDNk/aIidpzo1CxtEZl1FBQz5PKL8gk62NWgcQoRmV1mVVBEe2b2RCsq8gDYfqIj6s8tIhJJsyoooj0ze6L8zBQWFGWy/Xhn1J9bRCSSZlVQeG1FZR7bj3dq1zsRmVUUFDNoZWU+rb1DNHQMeF2KiMiMUVDMoFWVoXGKbcc1TiEis4eCYgYtLskmI8XPlnoFhYjMHgqKGZTk97G6qoBNR9u8LkVEZMYoKGbYFfMLONTcS3vfsNeliIjMCAXFDLtyQQEAb6pXISKzhIJihi0tyyMt2cfGI+1elyIiMiMUFDMsJcnH5Rfls/GIehQiMjsoKCLgyvmFHDjVQ1vvkNeliIhcsFkVFF6u9TTRdYuLAXj1sPanEJH4N6uCwsu1nia6rDSXoqwUXtZGRiIyC8yqoIgVPp9xbXUxrx7SPtoiEv8UFBFy3eLQPto7G7SarIjENwVFhFxbXYzP0OknEYl7CooIyc9MYUVFHq8c1D7aIhLfFBQRdP3iOexq7KJVl8mKSBxTUETQ9YuLcQ5ePaTTTyISvxQUEaTLZEVkNlBQRND4ZbKHdZmsiMQvBUWEXbe4mM7+EfY0ejtbXETkfCkoIuzqhUUAvFHX6nElIiLnR0ERYcXZqVw8N5s3ahUUIhKfFBRRcM2iIjbXdzA4EvC6FBGRczargiJWVo890zWLChkeDbL1WIfXpYiInLNZFRSxsnrsmdbOLyTJZ7yu008iEodmVVDEqqzUJC4tzWH7cfUoRCT+KCiiZFl5HnsauwlqPoWIxBkFRZQsK8+ld2iUI629XpciInJOFBRRsrwiD4CdJ2JroF1E5O0oKKJkYXEWGSl+bWQkInFHQRElfp+xtCyXnQ3qUYhIfFFQRNHyijz2n+xmeDTodSkiItOmoIiiZeW5DAeCHDzV43UpIiLTpqCIouXl4QFtjVOISBxRUERReX46+RnJ7FJQiEgcUVBEkZmxrDyPXRrQFpE4oqCIsuXluRxq7qF/eNTrUkREpkVBEWXLyvMIOtjT2O11KSIi06KgiLJlFaGVbTVOISLxQkERZXOy0yjNTdPEOxGJGwoKD4QGtNWjEJH4oKDwwLKKXI619dPZP+x1KSIib2taQWFmXzCzHAv5tpltM7ObI13cuYrVrVDPNDbxTpfJikg8mG6P4hPOuW7gZiAf+AjwDxGr6jzF6laoZ7qsTAPaIhI/phsUFv73NuB7zrm9E47JOcpNT2ZBUSY7tDeFiMSB6QbFVjN7gVBQPG9m2YCWQL0Ay8pz1aMQkbgw3aD4JHA/sMY51w8kA78fsaoSwPKKPFp6hjjVNeh1KSIiU5puUFwFHHTOdZrZ7wH/G9B5kwuwTCvJikicmG5QfAPoN7PlwJ8AdcB3I1ZVAri0NIckn+n0k4jEvOkGxahzzgF3Av/unPs6kB25sma/tGQ/NSXZukRWRGLedIOix8y+ROiy2J+ZmY/QOIVcgOUVuew80Ukog0VEYtN0g+KDwBCh+RSngHLgnyJWVYJYVp5H9+Ao9W39XpciIjKpaQVFOBweB3LN7HZg0DmnMYoLtKxcE+9EJPZNdwmPDwBvAu8HPgBsMrP3RbKwRFBTkk1qko+dmngnIjEsaZr3e4DQHIoWADMrBn4BPBmpwhJBst/H0rJcth5r97oUEZFJTXeMwjcWEmFt5/C9MoWrFxWxu7GLrv4Rr0sRETmr6b7Z/7eZPW9mHzezjwM/A56LXFmJY311EUEHG+pavS5FROSspjuY/b+Ah4Fl4Y+HnXN/FsnCEsWKijyyUpN4rVZBISKxabpjFDjnngKeimAtCSnZ7+PKBQW8flhBISKxacoehZn1mFn3WT56zKw7WkXOduurizne3s9xzacQkRg0ZY/COadlOqJg7fwCALaf6KCyMMPjakRE3kpXLsWAhcVZpPh97GtSJ01EYo+CIgakJPlYNCeL/U09XpciIvJbFBQxYklpDvtOqkchIrFHQREjLpmXQ2vvEC092vFORGKLgiJGLJmXA6DTTyIScxQUMeI3QaHTTyISWxQUMSI3I5myvHSNU4hIzJlVQWFmd5jZw11d8bls9yXzsjlwSkEhIrFlVgWFc+4nzrl7c3NzvS7lvFSXZHO0tY+RQNDrUkRExs2qoIh31XOyGAk4jrX1eV2KiMg4BUUMqSkJrZhyqLnX40pERH5DQRFDFhZnYQaHmnWJrIjEDgVFDElP8VORn8HhFvUoRCR2KChiTE1JFofVoxCRGKKgiDG68klEYo2CIsboyicRiTUKihizeG7oyqe9mqEtIjFCQRFjFpdkk5HiZ9uxDq9LEREBFBQxJ8nvY0VFHlsUFCISIxQUMWj1Rfnsb+qmd2jU61JERBQUsWh1VQFBB9uPq1chIt5TUMSglZV5+Ay21CsoRMR7CooYlJ2WzOK5OWzVOIWIxAAFRYy6Yn4BW491MDQa8LoUEUlwCooYdfXCQgZGAmw/3ul1KSKS4BQUMeqKBYX4DDbUtnpdiogkOAVFjMpNT2ZZeR6vKyhExGMKihh2zaJCdjZ00TM44nUpIpLAFBQx7JpFRQSCjk1H2r0uRUQSmIIihq2qzCfJZ2zTxDsR8ZCCIoalJfu5eF42O07oyicR8Y6CIsatqMhjV0MXwaDzuhQRSVAKihi3vDyP3qFR6k5rH20R8YaCIsatrMwDYLtOP4mIRxQUMW5BURbZaUkapxARzygoYpzPZywvz2OHlvIQEY8oKOLAyso8DpzSRkYi4g0FRRxYo42MRMRDCoo4MLaR0WZtZCQiHlBQxIHstGQumZfDlnot5SEi0aegiBNrqgrYfryTkUDQ61JEJMEoKOLEmqoCBkYC7DvZ7XUpIpJgFBRxYk1VPgCbjrZ5XImIJBoFRZyYk5NG9ZwsXjusjYxEJLoUFHFkfXUxbx5tZ3Ak4HUpIpJAFBRxZH1NEUOjQTbr6icRiSIFRRy5Yn4BKX6fTj+JSFQpKOJIRkoSq6vyefXQaa9LEZEEEvNBYWYLzOzbZvak17XEgutqijlwqofGzgGvSxGRBBHRoDCzR8ysxcz2nHH8FjM7aGa1Znb/VI/hnDvinPtkJOuMJzcuKQHgpf3NHlciIoki0j2KR4FbJh4wMz/wdeBWYAlwj5ktMbOlZvbTMz7mRLi+uLOwOIsFRZm8uE9BISLRkRTJB3fOvWpmVWccXgvUOueOAJjZE8Cdzrm/B24/3+cys3uBewEqKyvP92Hiwk1LSnjkjaP0DI6QnZbsdTkiMst5MUZRBpyY8HlD+NhZmVmhmT0ErDSzL012P+fcw8651c651cXFxTNXbQy6cUkJIwHHywc1qC0ikRfRHsVMcM61AZ/xuo5Ysqoyn7k5afzgzePcsbzU63JEZJbzokfRCFRM+Lw8fEymye8zfv+aKjbUtbH3ZJfX5YjILOdFUGwGqs1svpmlAHcDz3pQR1y7e20lmSl+vvXqkRl9XOccLd2DBIJuRh9XROJXRE89mdkPgOuBIjNrAP7SOfdtM/sc8DzgBx5xzu2NZB2zUW56Mh+6opJvvXaU5u4h/vrOS6kpyT7nxwkGHXWne9l0tJ03wx+nugcpzU3j5kvnkp+RQmaqn3m56fzOxXNIT/ETCDr8PovA/0pEYpE5N/v+cly9erXbsmWL12VE3NBogMc2Huc/flVLeoqfn31+PbkZb38VVH1rHz/ccoJf7GvmWHs/w6OhzZDmZKdyxYJClpbl8NrhVrbUdzAwYQHCzBQ/SX4fvUOj1JRkc8/aCj5y5UWYKTRE4p2ZbXXOrT7r12ZTUJjZHcAdixYt+vThw4e9Lidqth/v4P0P/Zp3XDyHb/7e5fgm+Wt/cCTAf7xcx0Mv1xFwjqsXFrJkXg4Li7NYO7+AiwozfutNfzQQpG84wP6mbp7deRIjtDXrm0fb2Ha8k7VVBRTnpJLkM+bmpvGH1y2aVliJSGxJmKAYkyg9iokeef0of/PTfbz/8nI+elUVx9v7uWlJCSlJoWGojUfa+NKPd3O0tY87V5TywG2XMCcn7byfzznHoxvqeXRDPX6fEQg6GjsGWFCcyZ/dcjFtvcM0dQ2SmernygWFzMlJpSAjhSR/zK8aI5KQFBQJ4l9ePMSDL/2mJ1U9J4tPrZ9PY8cA//arWioLMvjbd1/G+urIzDPZUNvKvd/bSu/Q6Fm/Xp6fzoN3r2TniU7eqG1lOBBkbVUB71lVRnl+BgCtvUP8fM8prl5YyMLiLCDUqzna2seiOVk6zSUSIQqKBPLz3U30DQfITPHztz/bP7544J0rSvm79ywlMzWyU2dOdQ1yrK2P0rx0SnLSaOsbYkt9B+19w3zzlTpOdg0CsGhOFil+H/uaQnuAr51fwEggyO6GLkaDjoqCdH76+fV0D4xw3w93sPVYB1fML+D25aUUZaZw05ISkvw+OvqGyc9MeUsNwaBjYCQw/n/tHRolM8WvkBGZgoIiQY0GgjR2DtAzOMqlpTmev1G29Q7xjZfruH7xHNZVFwFwor2fp7c38tzuJnLSkrm8Kp+L52bzJz/aSWVBBg0dA6Qm+bjnikp+vK2B1t5hAFZV5pGfkcJLB1q478ZqvnBDNXsau/nKiwd5o7aVkYDjncvmUT0ni3//ZS3vuHgOX/nAclL8PlKTfONtcbpniMbOAZaW5eL3GUOjAVKT/IwGgpzuHWJuThqtvcNsrm/nupriiAetiFcUFBJ3vrOhnn9+4SB3rSzj3usWUpaXzkggSEffMBvq2viLZ/bgHCwrz2VDXRvF2amc7hkiLyOZ960qDz3Gr+sZCTiuXFDAlvoOgs4RdDAvN43qkmxOdg5Q29ILQGluGmZGY+cARVmp9A6NMDgSpCgrle6BEYYDQeYXZXLrZXPZcaITM1hZkc99N1aT5PfR0j3IV188hJlRUZBOQUYK77h4DiU5aTjnaOgYYCQQZEH4dJpIrFFQyKzT3jeM34zstCS+/qta9p7s5qqFhbxnVRk54YUSa1t6Odk5wLU1xWw73sGL+5rJSPZz4FQPx9r7mJuTzsrKPOblpvHc7iaS/T5qSrJp6hogMzWJyoIMdp7oJC8jhdVV+fzdz/ZzqnuQS0tz8fmMnSc6ubammCvmF/Dohnq6B0bISPHT0T8CQGqSj3WLith2vGP82BduqOa+G6vpHhzlKy8cpCAzhc9ev5CmzkGy0pIoykr1rE0lsSkoRGbA8GiQodHA+Iq93990nL94Zg+BoKOmJIt/u2cVi+dm0z88SkPHAN94uY4tx9pZU1XA5Rfls/VYBz/e1khpbhrDAUd73xBBB9lpSfQMjpKR4uePb6phbm4aG+ra2HSkjY9dXcWH1layq7GLlu5BTrQPcLilhxsvKeGysly++UodVy8q4ncvnetx60i8S5igSNR5FOKd7sER/GbTGrtwzvHk1gZeO9xKz+AI991YQ0f/MP+1tYFVlfm8fLBlfD/0tGQfC4qy2NfUTUaKn/7ht0587BsOkOQzRsNLrbxz6TxuuGQOrxw6zb6T3Tz2qSsouYDLnyXxJExQjFGPQuKRc469J7tJ9vsoz08nPdnPd35dz4GmHq5bXExVYSZzc9PITkvi8Y3HqDvdxyfWzefp7Y088vrR8au7RoKO62uK+eZHLudwSy9ffHIXo8EgH7uqiveuKh+fkDkwHMDhyEjRAL0oKERmvdFAkAOneqjIz+CJzcf5+58fYE1VPrsbu8hKTaIwM5WDzT3ctaqMj19dxRObT/DM9kbyM1P40R9cRWleOoGgY3djF4tLsklP8Xv9X5IoU1CIJJDRQJDPPLaN5u5BLi3N4Y9vqqE4O5V/+2UtX33xEBAaaL/1srm8tL+FgqwUrq0u5vXaVo629jEvN43bls6jsWOAd60o5bal837r8TXDfvZRUIgIAM/saKSzf4R3rygjNyOZLfXt3PfDHfQNjVJVlMl7Vpbx5NYG9p3sJj8zhdM9Q7x7RSmfvX4RO0508NTWRraf6OD6xXN48O4VOm01iygoROScBIKOoHP860uH+earR8ZXGK4pyWJ5eR5PbWsYX0xyX1M3dS29fOb6hXx6/QKS1duISwoKETlvp3uGeGZHIzUl2ayvLsLM+MW+Zh586TDH2vqoKMigMCuVVw+dprIgg49edRG56cmMBBwd/cP88kALlQUZfOX9yydd2Vi8p6AQkYj75YFmHnyplp0nOt9yvKowg/q2fu6/9WI+c91Cj6qTt5MwQaF5FCLeGluuBCDZ7yM92U9OehJ/9P1tPL+3mfL8dMrz07lnbSWleemUhRePFO8lTFCMUY9CJLZ0D47wLy8eor1vmK3HOsbDJC3Zxz+/fzlHTvdR29LLp9cvYGl5rsfVJiYFhYjEjEDQsaW+nb7hUR78xWF2NnQBkJWaRO/QKCsr81h9UT45acncfOlcFs89973g5dwpKEQkJg0MB3jolTquXFDIZWU5fPfXx3hxXzP7TnYzHAiS7Df+8PpFfOiKSkpy0ugfHuX7m45TmJXCncvL8IV3Vzza2seCokwNll8ABYWIxJ32vmG+/MwefrqrCbPQDon9QwHa+kJ7kiyvyOOy0hw21LWN74D4obWVrK7KZ1l5nsfVxx8FhYjErdqWHv57zylqW3oZDgT5xDXzqTvdy3++dpTTvUNUFYb2CXl6eyMHTvUA8MBtl/DpaxcA8NArdbx5tJ3P/84iTvcMcby9n6VlueRmJFOclUqhlnYHFBQikgCcc5zqHuQvn9nLLw+08O8fWsXuxk6+/qs6UpJ845MGJ0pP9vPUZ6+mJCeVH21poLall5uWlHDLZXM51NxDZmoSZXnpjAZC3zubly5RUIhIwugaGOG2B18b3y/+rlVlfPn2JTy9vZGqwkwuLcth78lu+oZG+Zuf7CMzNYmRQJCGjgGyUpPoHx7l7rWV/GjzCVKSfHzkqot4ZvtJMlL8fO3uFSwty/V8W+FIUFCISEJp6R5kc30HpXlpLC/Pm3SQe0NtKx/+9iYKM1P4z4+toaYki4898iab6zu4tqaY4dEAG4+0s7wij5buQZq6BgGYX5TJB9dU8P7Ly8dPXZ3sHKCpa4BVlfnjQeKcY9PRdo639fO+y8vfUkd7X2jW+p0rSmNi2RMFhYjIJDbXt1ORn8Hc3NDEv96hUTbUtnLDJSUAHDzVwyXzsukeGOWJzcfpGxpl49F23jzaTorfxzsuLmZBcRbf2VBP/3CAS0tzuOGSEpxzPLe7ibrTfQD87qUlfO2DK0lP8XOyc4Df+/YmjpzuY92iIr54y2IOnOrhkdeP4jPjrlVlHGruoTQvnS/cUP2WHszgSIAj4cc0g8LMFOaE92Y/1tZPVVHmebVDwgSFZmaLSLQcbu7h8U3HeWHvKU52DXJtTTE3LynhsY3HONQcGlRfU1XAey8vp3tghP/z3H7K8tJ559J5/GjLCUYDjo9efRHffOXI+E6FF8/Nxjk42NxDWrKPwZEgf/SOhfzpzYtxDt6oa+WBp/dwvL1/vA6z0PM0dgzQ2jvEtr+4aVo7Lp4pYYJijHoUIhItzjna+4YpyEwZ/8t/cCTA0GiQ3PTk8fttPNLGXz6zl4PNPVxXU8yf33YJi+dmc+BUN/Wt/czNTWN5eFb60dbQYotffmYPP3jzBFmpSZhBz+AoFxVm8IUbqslI8Y+Hys92NVGen85tS+dxx/JS0pLPfeMpBYWISAwYDQRp7x9mTvb01rcKBB1Pb29kV0MngaBjVWU+ty2dF5EdCKcKCu06IiISJUl+37RDAsDvM953eTnvu7w8glW9Pe+H2kVEJKYpKEREZEoKChERmZKCQkREpqSgEBGRKSkoRERkSgoKERGZkoJCRESmNCtnZpvZaeDYeX57EdA6g+XMFNV17mK1NtV1bmK1Lojd2s6nroucc8Vn+8KsDIoLYWZbJpvG7iXVde5itTbVdW5itS6I3dpmui6dehIRkSkpKEREZEoKit/2sNcFTEJ1nbtYrU11nZtYrQtit7YZrUtjFCIiMiX1KEREZEoKChERmZKCIszMbjGzg2ZWa2b3e1xLhZn9ysz2mdleM/tC+PhfmVmjme0If9zmQW31ZrY7/PxbwscKzOxFMzsc/jc/yjUtntAmO8ys28zu86q9zOwRM2sxsz0Tjp21jSzkX8Ovu11mtirKdf2TmR0IP/fTZpYXPl5lZgMT2u6hKNc16c/OzL4Ubq+DZva7Ua7rhxNqqjezHeHj0Wyvyd4fIvcac84l/AfgB+qABUAKsBNY4mE984BV4dvZwCFgCfBXwJ963Fb1QNEZx/4vcH/49v3AP3r8szwFXORVewHXAquAPW/XRsBtwM8BA64ENkW5rpuBpPDtf5xQV9XE+3nQXmf92YV/D3YCqcD88O+tP1p1nfH1rwBf9qC9Jnt/iNhrTD2KkLVArXPuiHNuGHgCuNOrYpxzTc65beHbPcB+oMyreqbhTuA74dvfAd7tYS03AHXOufOdmX/BnHOvAu1nHJ6sje4EvutCNgJ5ZjYvWnU5515wzo2GP90IRH3PzUnaazJ3Ak8454acc0eBWkK/v1Gty8wM+ADwg0g891SmeH+I2GtMQRFSBpyY8HkDMfLGbGZVwEpgU/jQ58Ldx0eifYonzAEvmNlWM7s3fKzEOdcUvn0KKPGgrjF389ZfXq/ba8xkbRRLr71PEPrLc8x8M9tuZq+Y2XoP6jnbzy5W2ms90OycOzzhWNTb64z3h4i9xhQUMczMsoCngPucc93AN4CFwAqgiVDXN9rWOedWAbcCf2Rm1078ogv1dT255trMUoB3Af8VPhQL7fVbvGyjyZjZA8Ao8Hj4UBNQ6ZxbCfwx8H0zy4liSTH5s5vgHt76B0nU2+ss7w/jZvo1pqAIaQQqJnxeHj7mGTNLJvQieNw592MA51yzcy7gnAsC3yJCXe6pOOcaw/+2AE+Ha2ge68qG/22Jdl1htwLbnHPN4Ro9b68JJmsjz197ZvZx4Hbgw+E3GMKndtrCt7cSGguoiVZNU/zsYqG9koC7gB+OHYt2e53t/YEIvsYUFCGbgWozmx/+q/Ru4Fmvigmf//w2sN8599UJxyeeV3wPsOfM741wXZlmlj12m9BA6B5CbfWx8N0+BjwTzbomeMtfeV631xkma6NngTe4yJgAAAMJSURBVI+Gr0y5EuiacPog4szsFuCLwLucc/0TjhebmT98ewFQDRyJYl2T/eyeBe42s1Qzmx+u681o1RV2I3DAOdcwdiCa7TXZ+wORfI1FY5Q+Hj4IXRlwiNBfAg94XMs6Qt3GXcCO8MdtwPeA3eHjzwLzolzXAkJXnOwE9o61E1AIvAQcBn4BFHjQZplAG5A74Zgn7UUorJqAEULngz85WRsRuhLl6+HX3W5gdZTrqiV0/nrsdfZQ+L7vDf+MdwDbgDuiXNekPzvggXB7HQRujWZd4eOPAp85477RbK/J3h8i9hrTEh4iIjIlnXoSEZEpKShERGRKCgoREZmSgkJERKakoBARkSkpKESmYGYbwv9WmdmHZvix//xszyUSa3R5rMg0mNn1hFYzvf0cvifJ/WbBvbN9vdc5lzUT9YlEknoUIlMws97wzX8A1of3GvifZua30F4Om8ML1/1B+P7Xm9lrZvYssC987P+FF1HcO7aQopn9A5AefrzHJz5XeAbtP5nZHgvt/fHBCY/9spk9aaE9JB4Pz9IViagkrwsQiRP3M6FHEX7D73LOrTGzVOANM3shfN9VwGUutAw2wCecc+1mlg5sNrOnnHP3m9nnnHMrzvJcdxFaDG85UBT+nlfDX1sJXAqcBN4ArgFen/n/rshvqEchcn5uJrR+zg5CSzwXElrfB+DNCSEB8D/MbCeh/R4qJtxvMuuAH7jQonjNwCvAmgmP3eBCi+XtILRhjkhEqUchcn4M+Lxz7vm3HAyNZfSd8fmNwFXOuX4zexlIu4DnHZpwO4B+hyUK1KMQmZ4eQttOjnke+Gx4uWfMrCa8ou6ZcoGOcEhcTGgryjEjY99/hteAD4bHQYoJbckZ7RVSRcbprxGR6dkFBMKnkB4FHiR02mdbeED5NGffAva/gc+Y2X5Cq51unPC1h4FdZrbNOffhCcefBq4itEqvA77onDsVDhqRqNPlsSIiMiWdehIRkSkpKEREZEoKChERmZKCQkREpqSgEBGRKSkoRERkSgoKERGZ0v8HHdJTsLhrTN8AAAAASUVORK5CYII=\n","text/plain":["<Figure size 432x288 with 1 Axes>"]},"metadata":{"tags":[],"needs_background":"light"}}]},{"cell_type":"markdown","metadata":{"id":"U3S1NOmjBRUB"},"source":["# Model evaluation"]},{"cell_type":"code","metadata":{"colab":{"base_uri":"https://localhost:8080/"},"id":"APwvXRj-Be1Z","executionInfo":{"status":"ok","timestamp":1612978792581,"user_tz":-180,"elapsed":26182,"user":{"displayName":"Vl F","photoUrl":"","userId":"14485856366916106708"}},"outputId":"6bfd87dd-f5d3-4fcc-cfd3-3fce06f97c91"},"source":["N_batch = 10\n","\n","for J in [3, 4, 5, 6, 7, 8, 9, 10, 11]:\n","  N_levels = J - 1\n","\n","  # serialization\n","  projection_filters_ = [proj_var for i in range(N_levels)]\n","  smoothing_filters_pre_ = [smoothing_pre_var for i in range(N_levels)]\n","  smoothing_filters_post_ = [smoothing_post_var for i in range(N_levels)]\n","\n","  strides_ = [[2, 2]]*N_levels\n","  shapes_ = compute_shapes(J, projection_filters_, strides_, batch_size=N_batch)\n","  spectral_radius = stochastic_trace_r(A, J, projection_filters_, shapes_, strides_, smoothing_filters_pre_, smoothing_filters_post_, N_sweeps=10, batch_size=10).numpy()\n","  print(\"J = \", J, \" spectral radius = \", spectral_radius)"],"execution_count":10,"outputs":[{"output_type":"stream","text":["J =  3  spectral radius =  0.045686502\n","J =  4  spectral radius =  0.051948443\n","J =  5  spectral radius =  0.05819692\n","J =  6  spectral radius =  0.066308185\n","J =  7  spectral radius =  0.073735125\n","J =  8  spectral radius =  0.08033249\n","J =  9  spectral radius =  0.08683526\n","J =  10  spectral radius =  0.09427497\n","J =  11  spectral radius =  0.10188268\n"],"name":"stdout"}]}]}