s3MG(s).ipynb

{"nbformat":4,"nbformat_minor":0,"metadata":{"colab":{"name":"s3MG(s).ipynb","provenance":[],"collapsed_sections":[],"authorship_tag":"ABX9TyOT6UtutO6vSrbJR6ItL/X2"},"kernelspec":{"name":"python3","display_name":"Python 3"},"accelerator":"GPU"},"cells":[{"cell_type":"code","metadata":{"id":"biGrpEVhCfa0","executionInfo":{"status":"ok","timestamp":1612979361103,"user_tz":-180,"elapsed":2682,"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":"EJHPcclGC4PD"},"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 three layers."]},{"cell_type":"code","metadata":{"id":"Y_sd-pxGCkOS","executionInfo":{"status":"ok","timestamp":1612979363601,"user_tz":-180,"elapsed":1519,"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":2,"outputs":[]},{"cell_type":"markdown","metadata":{"id":"OWQkMqPSAf03"},"source":["# Equations"]},{"cell_type":"code","metadata":{"id":"1try7DlJAfXZ","executionInfo":{"status":"ok","timestamp":1612979370120,"user_tz":-180,"elapsed":6747,"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":3,"outputs":[]},{"cell_type":"markdown","metadata":{"id":"pgSMtUfnDhOD"},"source":["# Training"]},{"cell_type":"code","metadata":{"id":"vWSyjZqWDkqZ","executionInfo":{"status":"ok","timestamp":1612979389388,"user_tz":-180,"elapsed":7098,"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 three layers -- three variables for presmoothing, three variables 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_1 = 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 1\")\n","smoothing_post_var_1 = 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 1\")\n","smoothing_pre_var_2 = 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 2\")\n","smoothing_post_var_2 = 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 2\")\n","smoothing_pre_var_3 = 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 3\")\n","smoothing_post_var_3 = 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 3\")\n","smoothing_pre_var = [smoothing_pre_var_1, smoothing_pre_var_2, smoothing_pre_var_3]\n","smoothing_post_var = [smoothing_post_var_1, smoothing_post_var_2, smoothing_post_var_3]\n","\n","projection_filters = [proj_var for i in range(N_levels)]\n","smoothing_filters_pre = [smoothing_pre_var[i%3] for i in range(N_levels)]\n","smoothing_filters_post = [smoothing_post_var[i%3] 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":5,"outputs":[]},{"cell_type":"code","metadata":{"id":"opSx8L2WDylZ","executionInfo":{"status":"ok","timestamp":1612979472244,"user_tz":-180,"elapsed":81098,"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":6,"outputs":[]},{"cell_type":"code","metadata":{"colab":{"base_uri":"https://localhost:8080/","height":296},"id":"GwLJvaI5D1gJ","executionInfo":{"status":"ok","timestamp":1612979472710,"user_tz":-180,"elapsed":80370,"user":{"displayName":"Vl F","photoUrl":"","userId":"14485856366916106708"}},"outputId":"ebe6e403-1dd6-4c46-956f-c40b2c197f18"},"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":7,"outputs":[{"output_type":"stream","text":["final loss =  0.037872486\n"],"name":"stdout"},{"output_type":"display_data","data":{"image/png":"iVBORw0KGgoAAAANSUhEUgAAAYoAAAEGCAYAAAB7DNKzAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3deXSc9X3v8fd3NKN9Xy1LtiXLO97xgo2BBAyBgHEDbQJJE5LQpMm9SZPe3PaSprknt+empe1pbpuGJpfcAE1CIAuQAIFAgAQHY4z3DRvvi+RFsjZrl0bzu3/MiAjHGmSjmWc083mdo+PRoxnNx49G89Gz/H6POecQEREZic/rACIikthUFCIiEpWKQkREolJRiIhIVCoKERGJyu91gFgoLS11NTU1XscQERk3tmzZctY5V3ahryVlUdTU1LB582avY4iIjBtmdmykr2nXk4iIRKWiEBGRqFQUIiISlYpCRESiUlGIiEhUKgoREYlKRSEiIlGpKIb5wYajrD941usYIiIJRUUR0R8M8fDG49z1wOs8tqXe6zgiIglDRRGR7vfx4z9fwbLaYr700x28vL/J60giIglBRTFMQVaABz6+lLqyHP7m8V109QW9jiQi4jkVxXkyA2nce/t8Gtp6+OaLB7yOIyLiORXFBSytKWbtwok8vPE43f3aqhCR1KaiGMGHl02msy/IL3ee8jqKiIinVBQjWFZbzNSyHB7ddMLrKCIinlJRjMDMuGPpJLYca2X/mQ6v44iIeEZFEcVti6sJpBmPvq6tChFJXSqKKEpzM7h+TgWPb6unLzjodRwREU8kfFGYWY6Z/aeZfdfMPhLv579j6WTaugd4bs+ZeD+1iEhC8KQozOwBM2s0s93nLb/RzN40s4Nmdk9k8W3Az5xznwJujXfWVdNKqS7K4qebtftJRFKTV1sUDwE3Dl9gZmnAfcBNwBzgTjObA1QDQ+/Scd//4/MZN8+vZMOhZs71DsT76UVEPOdJUTjn1gEt5y1eBhx0zh12zvUDjwJrgXrCZQFR8prZp81ss5ltbmoa23maVs+uIBhy/G6/ZpYVkdSTSMcoqvj9lgOEC6IKeBy43cy+DTw10oOdc/c755Y455aUlZWNabBFkwopzA7w4j4dpxCR1OP3OsA7cc51AZ/wMoM/zcd7ZpTx2zebGAw50nzmZRwRkbhKpC2KBmDSsM+rI8sSwrWzK2jp6mf7iVavo4iIxFUiFcUmYLqZ1ZpZOnAH8KTHmd5yzfQy0nzGi3sbvY4iIhJXXp0e+wiwAZhpZvVmdrdzLgh8DngO2Av8xDm3x4t8F1KQHWBpTREv7VNRiEhq8eQYhXPuzhGWPwM8c6nf18zWAGumTZt2qd8iqutmVfD1Z/ZS39pNdVF2TJ5DRCTRJNKup3fNOfeUc+7TBQUFMfn+184uB+A32qoQkRSSVEURa1NLc6gpyeYFHacQkRSiorgIZsZ1syvYcKhZ19MWkZShorhI180qp38wxPqDGqUtIqlBRXGRltQUk5fh19lPIpIykqoozGyNmd3f3t4es+dI9/u4ekYZL+1rJBRyMXseEZFEkVRFEeuznoZcO6ucxo4+dp+MXSGJiCSKpCqKeHnvrHLM0ChtEUkJKopLUJyTzuLJGqUtIqlBRXGJrp1Vzq6Gds6c6/U6iohITKkoLtF1GqUtIilCRXGJZlbkUVWYpVHaIpL0VBSXyMy4fk4F6w400dLV73UcEZGYSaqiiMc4iuE+vHwy/cEQP9504p3vLCIyTiVVUcRrHMWQGRV5XDG1mB++doxBDb4TkSSVVEXhhY+tqKGhrYffvqljFSKSnFQU79L1cyrIy/Tz/J4zXkcREYkJFcW7FEjzcWVdKesONOGcdj+JSPJRUYyBq2eUcaq9l4ONnV5HEREZcyqKMXD1jFIAXt7f5HESEZGxp6IYA9VF2dSV5bDugC5mJCLJR0UxRq6eUcbGw830Dgx6HUVEZEwlVVHEe8DdcFfWldIXDLHteFvcn1tEJJaSqijiPeBuuGVTi/EZbDik3U8iklySqii8lJ8ZYF51Ia8eavY6iojImFJRjKGVdSVsP9FGV1/Q6ygiImNGRTGGVtaVEAw5Nh1t8TqKiMiYUVGMoSVTikn3+3QxIxFJKiqKMZSVnsbq2eX8ctcpgoMhr+OIiIwJFcUYW7uwirOd/bxyUGc/iUhyUFGMsffMLCM/088vtp/0OoqIyJhQUYyxDH8aN8+v5Lk9p+np1yhtERn/kqoovByZPdwt8yfS3T+oSQJFJCkkVVF4OTJ7uOW1xRRlB/jV7lOe5hARGQtJVRSJwp/m4/o5Fby4t5G+oHY/icj4pqKIkZvmVdLRF2S9zn4SkXFORREjV9aVkpfh54W9GnwnIuObiiJG0v0+lk8tZoMmCRSRcU5FEUMr6ko5craLk209XkcREblkKooYWllXAqCtChEZ11QUMTSzIo/inHRdo0JExjUVRQz5fMaKqSVsOHQW55zXcURELomKIsZW1JVwsr2XY83dXkcREbkkKooYGzpOod1PIjJeqShirLY0hwn5mbx6SAPvRGR8SqqiSJRJAYczM1bWlbDhULOOU4jIuJRURZEokwKeb0VdCc1d/ew/0+l1FBGRi5ZURZGoVkSOU2jeJxEZj1QUcVBdlM2Ukmwd0BaRcUlFEScr60rYeLiZ4GDI6ygiIhdFRREnK+pK6egLsufkOa+jiIhcFBVFnKyYqvEUIjI+qSjipCwvgxkVuRpPISLjjooijlbWlbLpaAv9QR2nEJHxQ0URRyvqSugdCLH9RJvXUURERk1FEUdX1JZghnY/ici4oqKIo4LsAHMnFuiAtoiMKyqKOFtZV8K246309A96HUVEZFRUFHG2oq6EgUHH5mMtXkcRERkVFUWcLa0pxu8z1h/U7icRGR9UFHGWk+FnXnUBW7RFISLjRFIVRSJej+JCFk0qYmd9OwOa90lExoGkKopEvR7F+RZOLqQvGOLN0x1eRxEReUdJVRTjxaJJhQBs08A7ERkHVBQeqC7KojQ3ne3HVRQikvhUFB4wMxZOKmTbiVavo4iIvCMVhUcWTS7icFMX7d0DXkcREYlKReGRy6cUAbDhsMZTiEhiU1F45PIpReRl+nlp3xmvo4iIRKWi8Eggzcc1M8p4aV8ToZDzOo6IyIhUFB66bnY5Zzv72NWQ2AMERSS1qSg89J4Z5fgMXtzX6HUUEZERqSg8VJSTzqLJRazb3+R1FBGREY2qKMzsC2aWb2HfM7OtZnZDrMOlgqU1xew52U7vgK5PISKJabRbFJ90zp0DbgCKgI8C98YsVQpZWlPEwKBjh6bzEJEENdqisMi/7wd+4JzbM2yZvAtD4yk2H9MobRFJTKMtii1m9jzhonjOzPIAzZE9Bgqz05lensumo7o+hYgkJv8o73c3sBA47JzrNrNi4BOxi5ValtQU8fTOU4RCDp9PG2oiklhGu0WxAnjTOddmZn8K/C2gk//HyJIpxXT0BtnfqOtTiEjiGW1RfBvoNrMFwJeAQ8D3Y5YqxSyOHKfYpmnHRSQBjbYogs45B6wFvuWcuw/Ii12s1FJTkk1hdkDXpxCRhDTaYxQdZvZlwqfFXmVmPiAQu1ipRdenEJFENtotig8BfYTHU5wGqoF/jlmqFLRwUiEHGjvp6NX1KUQksYyqKCLl8DBQYGa3AL3OOR2jGEOLJhfhHOyq1zkCIpJYRjuFxweB14E/AT4IbDSzP45lsEthZmvM7P729vH3ZruwuhCAbRqhLSIJZrS7nr4CLHXO3eWc+xiwDPhq7GJdGufcU865TxcUFHgd5aIVZAeYWpqjM59EJOGMtih8zrnhc2E3X8RjZZQWTyliy7EWXchIRBLKaN/sf2Vmz5nZx83s48AvgWdiFys1La8tprV7gINNnV5HERF5y6hOj3XO/ZWZ3Q5cGVl0v3PuidjFSk3La0sA2Hi4mRkVGqYiIolhtOMocM49BjwWwywpb1JxFpUFmWw80sJHV9R4HUdEBHiHojCzDuBCO8wNcM65/JikSlFmxrLaYl491IxzDjNNECgi3ot6jMI5l+ecy7/AR55KIjaW15bQ1NHH0eZur6OIiAA6cynhLKstBsLHKUREEoGKIsHUleVQmpvOxiO6kJGIJAYVRYIZOk7xuopCRBKEiiIBLa8toaGthxMtOk4hIt5TUSSgoeMU2qoQkUSgokhAMyvyKMgKsPGIDmiLiPdUFAnI53v7eAoRES+pKBLUe2eWU9/aw/4zmvdJRLylokhQ180uB+CFvWc8TiIiqU5FkaAq8jNZUF3Ar99QUYiIt1QUCez6ORVsP9FG47ler6OISApTUSSw1XMqAPi1dj+JiIdUFAlsZkUetaU5PLPrlNdRRCSFqSgSmJlx87xKNhxq5mxnn9dxRCRFqSgS3M3zKwk5eG7Paa+jiEiKUlEkuFkT8phamsPTO7T7SUS8oaJIcGbGbYur2HC4mV317V7HEZEUpKIYBz62soaCrAD/+sJ+r6OISApSUYwD+ZkBPn31VF7c18j2E21exxGRFKOiGCfuWllDXqafB9cf8TqKiKQYFcU4kZvh57ZFVTy76zQtXf1exxGRFKKiGEc+vHwK/YMhHttS73UUEUkhKopxZOaEPJZMKeLhjccYGAx5HUdEUoSKYpz582vqONrczUPrj3odRURShIpinFk9u5zrZpXzf17Yz8m2Hq/jiEgKUFGMM2bG1269DOfgzu++xv4zHV5HEpEkp6IYhyYVZ/Pwp5bT3T/In3xnA+d6B7yOJCJJLOGLwsymmtn3zOxnXmdJJIsnF/Hdjy2hvWeAJ7Y2eB1HRJJYTIvCzB4ws0Yz233e8hvN7E0zO2hm90T7Hs65w865u2OZc7xaOKmQBdUF/OC1YzjnvI4jIkkq1lsUDwE3Dl9gZmnAfcBNwBzgTjObY2bzzOzp8z7KY5xv3PvTK6ZwsLGT1w63eB1FRJJUTIvCObcOOP8dbBlwMLKl0A88Cqx1zu1yzt1y3kdjLPMlgzULJlKQFeCR1497HUVEkpQXxyiqgBPDPq+PLLsgMysxs+8Ai8zsy1Hu92kz22xmm5uamsYubYLLDKSxZkElz79xmg4d1BaRGEj4g9nOuWbn3Gecc3XOuX+Icr/7nXNLnHNLysrK4hnRcx9YVE3vQIhnd+sqeCIy9rwoigZg0rDPqyPL5BItnlxITUk2j2/VHFAiMva8KIpNwHQzqzWzdOAO4EkPciQNM+P2xdW8driFHbpehYiMsVifHvsIsAGYaWb1Zna3cy4IfA54DtgL/MQ5tyeWOVLBXVfWUJGfwf94bKcmDBSRMRXrs57udM5VOucCzrlq59z3Isufcc7NiBx3+HosM6SK/MwAf7d2LvtOd+jiRiIyphL+YPbFMLM1ZnZ/e3u711E88b7LJnDV9FL+3++O0B/UVoWIjI2kKgrn3FPOuU8XFBR4HcUzn1xVS2NHH8/sOuV1FBFJEklVFALXTC9jalkOD6w/EnVaD+ccO+vb+MX2Bp7acZKjZ7twztHTP8grB87S3q0xGSIS5vc6gIwtn8+4e1UtX3liN3/zxC7+161zSfe//e+Bc70DfPXnu/nF9pNvW56f6cc56OgLMiE/k3+4fR7vmVGGmcXzvyAiCUZFkYTuXDqZU229fOs3B9lz8hx//4F5zK0qYPPRFr7z8iFeOXiW/mCIL66ezs3zKukLhtjd0M7OhnZCIccVU0v495cO8IkHNzGvqoA7lk3ixssmUJKb4fV/TUQ8YMk46+iSJUvc5s2bvY7huWd3neKrv9jN2c5+CrICtPcMUJqbwS3zK7ltcRXzqwtHfGzvwCCPba3nwfVHOdjYSSDNeP+8StYunMiy2hJyM/Q3hkgyMbMtzrklF/xaMhWFma0B1kybNu1TBw4c8DpOQmjvHuBnW+s52NjJlJJsPrZiCtnpo3+Td86x73QHP9l8gp9trqejL0hehp8HP7GUJTXFMUwuIvGUMkUxRFsUsdE7MMiWY6189ee7OXOul79bO5fVcyooyAp4HU1E3qVoRaGznmTUMgNpXDmtlB996goqC7P40k93cMXfv8hzezQZoUgyU1HIRZtQkMlzX7yaxz67gpkT8vjsD7fwv59+g131qTnQUSTZqSjkkqT5jMunFPOjTy3nlvkTefDVo6z51is8sU0z2IokGxWFvCvZ6X6+eecitv7t9SyrLebLj+/SDLYiSUZFIWOiIDvAfR9eTGFWOmvvW8/ab73Cr3afijo6XETGB50ML2OmLC+DX3zuSn6+rYGfbD7BZ364lcqCTKZX5LFmfiVrFkwkM5DmdUwRuUhJdXqsxlEkjuBgiJ9vP8krB5rYUd/OkbNdTCrO4tsfuZy5Vak7aaNIotI4CvGUc451B85yz2M7ae7q589W1fKpq6ZSlJPudTQRidA4CvGUmXHNjDKe/vwq3nfZBL798iGWfv0F7rh/A3tO6pRakUSnLQqJu/1nOnhiWwOPbakn5OA/PrKY9QfPsnp2BfOqtVtKxAva9SQJ6WBjB7f9x6uc6w0CkBVI45/+eD4Lqgvp7AuS7jemled5nFIkNUQrCp31JJ6ZVp7H9+9ezstvNnHtrHLueXwnn39k29vu8/UPzOX6ORU8veMUG480s2paKR9dUeNNYJEUpS0KSRg9/YNsOHyW5s5+cjP8/HRLPS/taySQZgwMOvIz/XT0BfnBJ5czoSCD3IwAEwoyvY4tkhS060nGpb7gIF97cg+ZgTQ+esUUKvIz+aP71nP4bBeDIUd6mo9PXFnD++ZOYO7Egj+4kp+IjJ6KQpLGoaZOvvH8fpbVFrOzvp3HtobnlirPy+Dz107j1gVVFGRr2nORi5UyRaEBd6nnzLlethxr5cH1R9h0tBWfwQ1zJvCNDy0g059GMOS0pSEyCilTFEO0RZF6nHNsO9HGc3tO8911h5lbVcC5ngF6BgZ58OPLmDMx3+uIIglNA+4k6ZkZiycX8eWbZvNvdyxi76lzFGQFMIwP/d8N/GTzCQYGQ17HFBmXtEUhSamrL0h2ehon23v57A+3sLO+ndwMP6W56RRmpzO1NIcvvW8mlfmZ9AVDZKVrskJJbdr1JCnNOcdv3mzk5TebaOkeoK27n81HWzGDdL+Pzt4gn1xVy+evnUZepg6ES2pSUYic50RLN9/49X78PmMw5Hh8WwOF2QFuX1xNZ2+QpbXF3L64CjPzOqpIXKgoRN7Brvp2/vWF/by4r5G8jPDAvtsXV/MX101jSkmO1/FEYk5FITJKA4Mh0sz41xcP8O8vHcA5mFGRy+LJRSyeXMSiyYXUleXi82lLQ5KLikLkEpxs6+HpnSd59VAz24630d4zAEB+pp/3zipn9ewKlk8tpiw3Q7uoZNxTUYi8S6GQ40hzF1uPtfL6kRZe2HuG1u5wcWQGfEwtzeX6ORXctbKGYl2QScahlCkKjcyWeAkOhth98hybj7Zwur2XnfXtbDrWQmV+Jv/ywYUsrSnCn6ZhSjJ+pExRDNEWhXhhd0M7n/nhFupbe8gM+Fg9u4JPXFnD5VOKvY4m8o50PQqROJhbVcAv/+IqXnjjDNtOtPLk9pM8vfMU18+poCg7QHNnP19YPZ351YUAnOsdwGdGboZ+DSWxaYtCJEa6+4M8uP4o//Gbg/jTfPh9Rmt3P++dWU5upp9nd5/GOcfKulL+6n0zmVuly8CKd7TrScRD/cEQaT6jsy/IN188wMv7mzjT3sutCyeSk+Hn8a0NtHb3c/eqWv5y9Qx8PmjtGsBnUJ6vCzNJfKgoRBJYe/cA9/5qL4+8foKi7ABdfYP0RyYwnFuVz53LJnPzvEp+d+AsjR19rFlQSXmeCkTGlopCZBzYeLiZ7284RnVxFlOKc+jsG+DxrQ3sO93xtvv5fcYNl1Vww5wJAEwrz2VGRR6n2nvYfLSV0+d6WTN/IpNLsr34b8g4paIQGaecc+yob+f5PadZWlPMpOJsfrzpOD/dUk9bZBzHhZjBe2eWc/O8Sk629VCWl8Gcifk89OpR2rsHWDmtlD9ZUk2+JkGUCBWFSJLpHRjkeEs3PoOd9e0cbe6mqjCTeVWFFGYHePT14/zo9eOc7ex/2+Oy09OoyM/kyNkuSnLSuaKuhJNtPayeXUF+pp/HtjYwsyKPy6ryaWgNF0xdWS4F2QHyM/3kZwbIzwqQ4fdpNHqSUVGIpKD+YIgDjR1MKcnhcFMnO+rbuWnuBEpzM9hV3869v9rLseZuSnLS2VHfDsCsCXmcaOmmq3+QQJoxMHjh94f0NB/5WeHiuKyqgOtmlVOQFaB/MMTAYIiakhymV+SS4dd1PsYLFYWIRLX/TAddfUEWTiqkLxiitbufirxMmrv6Od7SzbneATp6g5zrGeBc7wDneoJ09A7Q2t3Pa4dbaOnq/4PvmZsRnhNrQXUBA4OOQ02drJpWyrWzyxkIhnh292lOtHQzqzKPyyYWMLU0B3+aj1DIMegcAY1sjysVhYjETHAwxP4znfQGB0lP85HmMw43dbFufxMvvdlIU0cfEJ5M8Vxv8G2P9fuMYCj8HmQGBVkBuvqCBNJ8/NlVU5lfVcDOhnYe31rPZRPzufe2+RTlpOOcoy8YIjOgLZaxkjJFobmeRBLP0NZGYVaAVw81s+dkeDfXldNKmTUhj0NNXbxxqp0jTV20dg+Qm+nnWHMXz+w6DYQLZHltMVuPtZEZ8FGal0FTRx8dvUGW1hQxc0Ie3f2DZKenkZ8ZICfDjxkUZaezoLqQnoEgB850su90B7Mr87jxskoKsgO0dPVz4EwHy2qLY3a8JTgYGjdzfqVMUQzRFoXI+HeoqZPO3iCVhZmU52Wyq76dh149Sm9wkKLsAPmZAV7a18iZc71kp/vp7g/S0Rt8awvlfOl+H/3BEBl+HzdcNoF1+5to7xlgZV0JfcEQb5w8xx8tqqK6KItDjZ0cOttFWW46182uYGppDp19QQ42dhIMOerKwrMFp/mM483drDvQxNXTy946Jdk5x+NbG/jaU3u4blY5//jH8986XuOcY92Bs7R29VNTmsOC6oKEODFARSEiKWFolxTAqfZedta3kZ8VoKYkhynF2ew+2c4jr5/giW31LJpUxNUzyvj2bw9SkpvB/OoCnt19mv5giAn5mdSV53D0bDcNbT0XfK7qoiz8PuNoczcQLqJrZpRxvLmbk+09dPQGmV6ey4HGTi6bmM/VM8rICqSx4VAzGw43v/V9ltUWs7KuhOPN3exqaCfkHFNKcjjV3ovfZyyeXMjlNcVcPqWIqsIsGtp6WLe/iS3HWlk1rZTr51Sw4VAzEwuzmF2Zd8mlo6IQERnGOffWG+rQe6CZ0dEbHpuSFxlf4pzjaHM39a3dZAbSmFGRR4bfxwt7z/CzLfXkZviZMzGfVdNK+d4rR9h2vI3p5blMKs5mTmU+t19ezS93neK+lw5yqCm8NVKam84XrpvOirpSXj10lm++eICznf1MyM9kblU+PjOOt3RTWZBJ70CI7Sfa6BkYBKAwO/DW+Jmc9DS6+gfxGQxtRE0qzuLn/+VKSnIzLnqdqChERDwWjEzLcv4xi+BgiEHnRjyVODgYYu+pDjYfa+GNk+eYOSGPa2aUUVeWy1M7T7Ln5Dmunl5GfWs32463ce/t8y5pq0JFISIiUUUrivFxOF5ERDyjohARkahUFCIiEpWKQkREolJRiIhIVCoKERGJSkUhIiJRqShERCSqpBxwZ2ZNwLFLfHgpcHYM44wV5bp4iZpNuS5OouaCxM12KbmmOOfKLvSFpCyKd8PMNo80OtFLynXxEjWbcl2cRM0FiZttrHNp15OIiESlohARkahUFH/ofq8DjEC5Ll6iZlOui5OouSBxs41pLh2jEBGRqLRFISIiUakoREQkKhVFhJndaGZvmtlBM7vH4yyTzOw3ZvaGme0xsy9Eln/NzBrMbHvk4/0eZDtqZrsiz785sqzYzH5tZgci/xbFOdPMYetku5mdM7MverW+zOwBM2s0s93Dll1wHVnYNyOvu51mtjjOuf7ZzPZFnvsJMyuMLK8xs55h6+47cc414s/OzL4cWV9vmtn74pzrx8MyHTWz7ZHl8VxfI70/xO415pxL+Q8gDTgETAXSgR3AHA/zVAKLI7fzgP3AHOBrwH/3eF0dBUrPW/ZPwD2R2/cA/+jxz/I0MMWr9QVcDSwGdr/TOgLeDzwLGHAFsDHOuW4A/JHb/zgsV83w+3mwvi74s4v8HuwAMoDayO9tWrxynff1fwH+pwfra6T3h5i9xrRFEbYMOOicO+yc6wceBdZ6FcY5d8o5tzVyuwPYC1R5lWcU1gL/Gbn9n8AfeZjlOuCQc+5SR+a/a865dUDLeYtHWkdrge+7sNeAQjOrjFcu59zzzrlg5NPXgOpYPPfF5opiLfCoc67POXcEOEj49zeuuSx8UeoPAo/E4rmjifL+ELPXmIoirAo4MezzehLkjdnMaoBFwMbIos9FNh8fiPcunggHPG9mW8zs05FlFc65U5Hbp4EKD3INuYO3//J6vb6GjLSOEum190nCf3kOqTWzbWb2spld5UGeC/3sEmV9XQWccc4dGLYs7uvrvPeHmL3GVBQJzMxygceALzrnzgHfBuqAhcApwpu+8bbKObcYuAn4r2Z29fAvuvC2rifnXJtZOnAr8NPIokRYX3/Ay3U0EjP7ChAEHo4sOgVMds4tAv4b8CMzy49jpIT82Q1zJ2//gyTu6+sC7w9vGevXmIoirAGYNOzz6sgyz5hZgPCL4GHn3OMAzrkzzrlB51wI+C4x2uSOxjnXEPm3EXgikuHM0KZs5N/GeOeKuAnY6pw7E8no+foaZqR15Plrz8w+DtwCfCTyBkNk105z5PYWwscCZsQrU5SfXSKsLz9wG/DjoWXxXl8Xen8ghq8xFUXYJmC6mdVG/iq9A3jSqzCR/Z/fA/Y6574xbPnw/YofAHaf/9gY58oxs7yh24QPhO4mvK7uitztLuAX8cw1zNv+yvN6fZ1npHX0JPCxyJkpVwDtw3YfxJyZ3Qj8NXCrc6572PIyM0uL3J4KTAcOxzHXSD+7J4E7zCzDzGojuV6PV66I1cA+51z90IJ4rq+R3h+I5WssHpEnPUMAAAK/SURBVEfpx8MH4TMD9hP+S+ArHmdZRXizcSewPfLxfuAHwK7I8ieByjjnmkr4jJMdwJ6h9QSUAC8CB4AXgGIP1lkO0AwUDFvmyfoiXFangAHC+4PvHmkdET4T5b7I624XsCTOuQ4S3n899Dr7TuS+t0d+xtuBrcCaOOca8WcHfCWyvt4Ebopnrsjyh4DPnHffeK6vkd4fYvYa0xQeIiISlXY9iYhIVCoKERGJSkUhIiJRqShERCQqFYWIiESlohCJwsxejfxbY2YfHuPv/TcXei6RRKPTY0VGwczeQ3g201su4jF+9/sJ9y709U7nXO5Y5BOJJW1RiERhZp2Rm/cCV0WuNfCXZpZm4Ws5bIpMXPfnkfu/x8x+Z2ZPAm9Elv08MoninqGJFM3sXiAr8v0eHv5ckRG0/2xmuy187Y8PDfvevzWzn1n4GhIPR0bpisSU3+sAIuPEPQzbooi84bc755aaWQaw3syej9x3MTDXhafBBvikc67FzLKATWb2mHPuHjP7nHNu4QWe6zbCk+EtAEojj1kX+doi4DLgJLAeuBJ4Zez/uyK/py0KkUtzA+H5c7YTnuK5hPD8PgCvDysJgL8wsx2Er/cwadj9RrIKeMSFJ8U7A7wMLB32vetdeLK87YQvmCMSU9qiELk0BnzeOffc2xaGj2V0nff5amCFc67bzH4LZL6L5+0bdnsQ/Q5LHGiLQmR0OghfdnLIc8BnI9M9Y2YzIjPqnq8AaI2UxCzCl6IcMjD0+PP8DvhQ5DhIGeFLcsZ7hlSRt+ivEZHR2QkMRnYhPQT8G+HdPlsjB5SbuPAlYH8FfMbM9hKe7fS1YV+7H9hpZludcx8ZtvwJYAXhWXod8NfOudORohGJO50eKyIiUWnXk4iIRKWiEBGRqFQUIiISlYpCRESiUlGIiEhUKgoREYlKRSEiIlH9f7r2jlUZvPrzAAAAAElFTkSuQmCC\n","text/plain":["<Figure size 432x288 with 1 Axes>"]},"metadata":{"tags":[],"needs_background":"light"}}]},{"cell_type":"markdown","metadata":{"id":"YkybvW-zDkUc"},"source":["# Model evaluation"]},{"cell_type":"code","metadata":{"colab":{"base_uri":"https://localhost:8080/"},"id":"UhxE23eLDma2","executionInfo":{"status":"ok","timestamp":1612979504725,"user_tz":-180,"elapsed":27226,"user":{"displayName":"Vl F","photoUrl":"","userId":"14485856366916106708"}},"outputId":"3384bd0f-2b51-4542-f581-9d7f74532a4e"},"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[i%3] for i in range(N_levels)]\n","  smoothing_filters_post_ = [smoothing_post_var[i%3] 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":8,"outputs":[{"output_type":"stream","text":["J =  3  spectral radius =  0.026539112\n","J =  4  spectral radius =  0.033393394\n","J =  5  spectral radius =  0.03776223\n","J =  6  spectral radius =  0.046081044\n","J =  7  spectral radius =  0.055237\n","J =  8  spectral radius =  0.061990507\n","J =  9  spectral radius =  0.068533845\n","J =  10  spectral radius =  0.07517537\n","J =  11  spectral radius =  0.08136343\n"],"name":"stdout"}]}]}