U-Net.ipynb

{"nbformat":4,"nbformat_minor":0,"metadata":{"colab":{"name":"U-Net.ipynb","provenance":[],"collapsed_sections":[],"authorship_tag":"ABX9TyNPbPrl2hULhdY9284EcU8n"},"kernelspec":{"name":"python3","display_name":"Python 3"},"accelerator":"GPU"},"cells":[{"cell_type":"code","metadata":{"id":"gNuX205rFSU1","executionInfo":{"status":"ok","timestamp":1612980622742,"user_tz":-180,"elapsed":2006,"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":"K2hkkbQJF5v4"},"source":["# The model\n","\n","+ Modified Richardson iteration with U-Net preconditioner\n","+ U-Net coincides with filtering preconditioner\n","+ No serialization is performed\n","+ U-Net architecture:\n","\n","U-Net = skip connection starts here $\\rightarrow$ convolution with stride $1$ $\\rightarrow$ convolution with stride $>1$ $\\rightarrow$ U-Net $\\rightarrow$ transposed convolution with stride $>1$ $\\rightarrow$ transposed convolution with stride $1$ $\\rightarrow$ skip connection ends here"]},{"cell_type":"code","metadata":{"id":"tJU5DGXkF5Fq","executionInfo":{"status":"ok","timestamp":1612980622744,"user_tz":-180,"elapsed":1314,"user":{"displayName":"Vl F","photoUrl":"","userId":"14485856366916106708"}}},"source":["# U-Net\n","def u_net(x, projection_filter_pre, projection_filter_post, pre_smoothing_filters, post_smoothing_filters, strides, N_layers, k=0):\n","  input_shape = tf.shape(x)\n","  if k == N_layers:\n","    y = x\n","    x = tf.nn.convolution(x, pre_smoothing_filters[k], strides=[1, 1], padding=\"SAME\")\n","    x = tf.nn.conv_transpose(x, post_smoothing_filters[k], x.shape, strides=[1, 1], padding=\"SAME\")\n","    return x + y\n","  y = x\n","  x = tf.nn.convolution(x, pre_smoothing_filters[k], strides=[1, 1], padding=\"SAME\")\n","  x = tf.nn.convolution(x, projection_filter_pre[k], strides=strides)\n","  x = u_net(x, projection_filter_pre, projection_filter_post, pre_smoothing_filters, post_smoothing_filters, strides, N_layers, k=k+1)\n","  x = tf.nn.conv_transpose(x, projection_filter_post[k], input_shape, strides=strides, padding=\"VALID\")\n","  x = tf.nn.conv_transpose(x, post_smoothing_filters[k], x.shape, strides=[1, 1], padding=\"SAME\")\n","  return x + y\n","\n","# modified Richardson solver\n","def u_net_solver(x, b, A, projection_filter_pre, projection_filter_post, pre_smoothing_filters, post_smoothing_filters, strides, N_levels):\n","  # diagonal is always 1\n","  r = b - A(x)\n","  x = x + u_net(r, projection_filter_pre, projection_filter_post, pre_smoothing_filters, post_smoothing_filters, strides, N_levels)\n","  return x\n","\n","# approximate spectral radius\n","def stochastic_trace_u_net_solver(A, J, N_levels, projection_filter_pre, projection_filter_post, pre_smoothing_filters, post_smoothing_filters, strides, 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_filter_pre[0].dtype)\n","  for _ in range(N_sweeps):\n","    x = u_net_solver(x, x*0, A, projection_filter_pre, projection_filter_post, pre_smoothing_filters, post_smoothing_filters, strides, N_levels)\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":1612980624246,"user_tz":-180,"elapsed":1598,"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":"ELQslfHcG92J"},"source":["# Training"]},{"cell_type":"code","metadata":{"id":"mSjiThilG_px","executionInfo":{"status":"ok","timestamp":1612980624248,"user_tz":-180,"elapsed":393,"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","u_net_strides = [2, 2]\n","\n","# training all present convolutions. note, that N_levels = J - 1 but additional pre- and postsmoothing filters are present!\n","pre_smoothing_filters = [tf.Variable(np.outer([1/2, 1, 1/2], [[1/2], [1], [1/2]]).reshape((3, 3, 1, 1))/2, dtype=tf.float32, trainable=True, name=\"pre_smoothing_filter_\"+str(i)) for i in range(N_levels+1)]\n","post_smoothing_filters = [tf.Variable(np.outer([1/2, 1, 1/2], [[1/2], [1], [1/2]]).reshape((3, 3, 1, 1))/2, dtype=tf.float32, trainable=True, name=\"post_smoothing_filter_\"+str(i)) for i in range(N_levels+1)]\n","projection_filter_pre = [tf.Variable(np.outer([1/2, 1, 1/2], [[1/2], [1], [1/2]]).reshape((3, 3, 1, 1))/2, dtype=tf.float32, trainable=True, name=\"pre_projection_filter_\"+str(i)) for i in range(N_levels)]\n","projection_filter_post = [tf.Variable(np.outer([1/2, 1, 1/2], [[1/2], [1], [1/2]]).reshape((3, 3, 1, 1))/2, dtype=tf.float32, trainable=True, name=\"post_projection_filter_\"+str(i)) for i in range(N_levels)]\n","\n","variables = pre_smoothing_filters + post_smoothing_filters + projection_filter_pre + projection_filter_post"],"execution_count":4,"outputs":[]},{"cell_type":"code","metadata":{"id":"uNYHT57UHa16","executionInfo":{"status":"ok","timestamp":1612980814458,"user_tz":-180,"elapsed":189813,"user":{"displayName":"Vl F","photoUrl":"","userId":"14485856366916106708"}}},"source":["opt = tf.keras.optimizers.Adam(learning_rate=0.01)\n","N_epoch = 2000\n","\n","history = np.zeros(N_epoch)\n","for i in range(N_epoch):\n","  with tf.GradientTape() as tape:\n","    loss = stochastic_trace_u_net_solver(A, J, N_levels, projection_filter_pre, projection_filter_post, pre_smoothing_filters, post_smoothing_filters, u_net_strides, N_sweeps=10, batch_size=10)\n","  opt.apply_gradients(zip(tape.gradient(loss, variables), variables))\n","  history[i] = loss.numpy()"],"execution_count":5,"outputs":[]},{"cell_type":"code","metadata":{"colab":{"base_uri":"https://localhost:8080/","height":296},"id":"l2V4LZ8DHcem","executionInfo":{"status":"ok","timestamp":1612980815054,"user_tz":-180,"elapsed":190076,"user":{"displayName":"Vl F","photoUrl":"","userId":"14485856366916106708"}},"outputId":"dd2a630c-f70f-4ab7-cfde-1bc028fb1649"},"source":["# training history\n","final_loss = stochastic_trace_u_net_solver(A, J, N_levels, projection_filter_pre, projection_filter_post, pre_smoothing_filters, post_smoothing_filters, u_net_strides, 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":6,"outputs":[{"output_type":"stream","text":["final loss =  0.64342713\n"],"name":"stdout"},{"output_type":"display_data","data":{"image/png":"iVBORw0KGgoAAAANSUhEUgAAAZsAAAEGCAYAAACzYDhlAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3dd5wU5f0H8M/3OvXoiIiCqEEsgCLR2AWVpmiMURNjjUrKLyYxKtagokGNsUVJsHeN2BtV6dJ7vwNP+nHAcXfccW33+f2xs3uzszO7M3c7s+U+79eLF7uzszPPze4+33m6KKVARETkpoxEJ4CIiNIfgw0REbmOwYaIiFzHYENERK5jsCEiItdlJToByapTp06qZ8+eiU4GEVFKWbp06V6lVGfjdgYbCz179sSSJUsSnQwiopQiIj+abWc1GhERuY7BhoiIXMdgQ0RErmOwISIi1zHYEBGR6xhsiIjIdQw2RETkOgYbD323cQ+2l1YlOhlERJ5jsPHQja8txsVPz050MoiIPMdg47HKWl+ik0BE5DkGGyIich2DDRERuY7BhoiIXMdgQ0RErmOwISIi1zHYEBGR6xhsiIjIdQw2RETkOgYbIiJyHYMNERG5jsGGiIhcx2BDRESuY7AhIiLXMdgQEZHrGGyIiMh1DDYeUUolOglERAnDYOMRxhoias4YbIiIyHUMNh5hwYaImjMGG4+wzYaImjMGG48w1BBRc8ZgQ0RErmOw8Qhr0YioOWOw8YhiRRoRNWMMNh5hyYaImjMGGyIich2DDRERuY7BxiOsRiOi5ozBxiPsIEBEzRmDDRERuY7BxiOsRiOi5ozBxiOMNUTUnDHYEBGR6xhsPMJZn4moOWOw8QhDDRE1Zww2HmHBhoiaMwYbIiJyHYONV1iyIaJmjMHGZaWVtVi+tRQHDtWGts0t2IvqOl8CU0VE5K2sRCcg3fzqpQXo3CYXF/TpguVbD+D1+UUR+1z7ykIc3akVtuytxLS/nIPl2w7gkpMPR4ucTO8TTETkAWGXXHMDBw5US5Yscfy+0x+bgd3l1Y7f1yonE9/+7Ty0ys1C61zeAxBRahKRpUqpgcbtzNXirDGBBgAqa3346WMzAABF40fEM0lERAnHNps4O6J9i7DnY4b1cXwMv5+lTSJKLyzZxNn0v56L9bvKUV5dj65tc3Fslzbolp+H/ZW1eOiLdbaOUVpVi46tc11OKRGRdxhs4iwvOxMDjmwftm1U/+74aOn20PNF9w1GlzZ5AID1u8qRnSkY8q/Zode/Xr0LvzmjpyfpJSLyAqvRPJKT1XCp9R0Aju/WFsd0aYOF9w5G17aB0syTUzZ6nj4iIjcx2HgkO7PhUudlRXZx7to2D7ecfTSA8MBERJQOmKt5JFcXQDIyxHSf63/WEwAwuE9XL5JEROQZBhuP5NoorWRnZuC4rq3DZhsgIkoHDDYesVs11io3C1W1nMqGiNILg41HjuzY0tZ+LbIzcYjBhojSDIONR4JdnWNpmZPJkg0RpR0GmyTTIicLhzgjNBGlGQYbj3Vv1yLq6y2zM7HvYA1Wby/zKEVERO5jsPHQpnHDMPPO86Lu0yInE+XV9bjk33NxoIq90ogoPXC6Gg/Z6ZGmX9Omus7vZnKIiDzDkk2SaZnNBdSIKP0w2CQZrtZJROmIwSbJ6IONmM9qQ0SUchhskkxLlmyIKA0x2CSZFtnss0FE6YfBJsnoSzaKq0MTUZpgsEky7CBAROmIwSbJmC2sRkSU6hhskgxX6SSidMScLclkZzb0d2bXZyJKF82i65OItALwIoBaADOVUu8kOEmW9CUbdhAgonSRsiUbEXlVRPaIyBrD9qEislFECkVkjLb55wAmKaVuAXCp54l1ICdTF2zAaENE6SFlgw2A1wEM1W8QkUwALwAYBqAvgGtEpC+AIwBs03ZL6sViWLIhonSUssFGKTUbwH7D5kEACpVSW5RStQDeBzAKwHYEAg6Q5H9zdljJhogoPSR1xtsI3dFQggECQaY7gI8BXCEiEwB8YfVmEblVRJaIyJKSkhJ3U2ohvGTDcENE6aFZdBBQSlUCuNHGfhMBTASAgQMHJiSnz8po6ILGWENE6SLdSjY7APTQPT9C25YyhP2diSgNpVuwWQzgWBHpJSI5AK4G8HmC09RoLNkQUbpI2WAjIu8B+B7AT0Rku4jcrJSqB/BHAFMArAfwP6XU2kSmsymKK6oTnQQiorhI2TYbpdQ1Ftu/BvC1x8lxxZX/+R5L7h+CTq1zE50UIqImSdmSTXNRUHww0UkgImoyBpskV3aoLtFJICJqMgabJHewpj7RSSAiajIGmyR3sJolGyJKfQw2SY4lGyJKBww2Sa6CwYaI0gCDTZI7WM1gQ0Spj8EmybEajYjSAYNNkmPJhojSAYNNkmPJhojSAYNNkquuCywsWlFdx/VtiChlMdgkuapaH2rqfThp7FTc+8kaBhwiSkkMNkmuqtaHfQdrAQDvLdqKM8d/m+AUERE5x2CThH710yNDj2t9/rD50XaWcdkBIko9DDZJKCez4WOpqfPB52fVGRGlNlvBRkRuF5G2EvCKiCwTkYvcTlxz1b9HOwBAvx7tUOvzM9gQUcqzW7K5SSlVDuAiAO0B/AbAeNdS1cyN6n845o25AGcf0wm19X742CmAiFKc3WAj2v/DAbylLbUsUfanJhARdG/XAjlZGfAroLben+gkERE1id1gs1REpiIQbKaISBsAaZkDisglIjKxrKws0UlBblbg4wmOtSEiSlV2g83NAMYAOE0pVQUgG8CNrqUqgZRSXyilbs3Pz090UpDDYENEacJusDkDwEal1AERuRbA/QASf+uf5oLBpqqWwYaIUpvdYDMBQJWI9ANwB4DNAN50LVUEoKEL9CGtZDPwqPYAgB0HDiUsTUREjWE32NSrwDwpowD8Wyn1AoA27iWLACA3OxMAUF0XaB4LlnQ4iwARpZosm/tViMg9CHR5PltEMhBotyEX5WnB5ZEv1wEAurTJTWRyiIgazW7J5ioANQiMt9kN4AgAT7qWKgIAdGwdHlx6dGiZoJQQETWNrWCjBZh3AOSLyEgA1Uopttm4rF3L8MJjmzy7BVEiouRid7qaXwJYBOBKAL8EsFBEfuFmwghoobXZBLVrmRN6/NmKHV4nh4io0ezeKt+HwBibPQAgIp0BTAcwya2EEZBnCDb9e7RDblYGaur9uP39FRjVv3uCUkZE5IzdNpuMYKDR7HPwXmqkvOzwS5yZIfjdeb0TlBoiosazW7KZLCJTALynPb8KwNfuJImC8rLCSzaZIsgQTklHRKnHVrBRSt0pIlcAOFPbNFEp9Yl7ySIAyMgQ5GRlhCbizBCBnzNAE1EKst29SSn1EYCPXEwLmcjTB5sMwM+1bYgoBUUNNiJSAcAsdxMASinV1pVUUUhedibKq+sBBNpsuLYNEaWiqMFGKcUpaRKsTV4W9lTUAAi02bBgQ0SpiD3KktwxXVqHHgvbbIgoRTHYJLl2LRoGcmZmCBhriCgVMdgkuda6KWoyReAzqUcrr67zMklERI4x2CS5r1btCj3OyEBEsCnaW4mTx07Fuwu3ep00IiLbGGyS3O7y6tDjDBEoXT3agapalBwMdB74cOk2z9NGRGQXg02SG3fZiaHHmRmCM3p3Cj2vqK5nGw4RpQQGmyR3qrYUNACIAENPPAx3D+0DIFClVlPvS1TSiIhsY7BJctmZDR9RcF60w9vlAQD8SoVmF1AqMLvA09M2Ya9WtUZElCwYbJJcji7YBAOPaEHHr4Cq2oaSzeKi/Xh2RgHumrTK20QSEcXAYJPksjIjZ3nO0DbV1PtQrHUg6NgqJxSEyg+xKzQRJReuM5zk9NVoQWVaMLn349VokRNYhiAvJxPZWmCq8/m9SyARkQ0s2RiIyCUiMrGsrCzRSQEQXo0WFBxrs3J7GRZs2Q8gMB6nzhfYHvyfiChZMNgYKKW+UErdmp+fn+ikADCvRhOLBdT+9uFKAEBFDavRiCi5MNgkObNqNKvBNVv3VwEAisvZG42IkguDTZLLNinZxKwkYy0aESUZBpskZ1ZlFmvWAMVoQ0RJhsEmBSnOUdNo09cVY3tpVaKTQdTssOtzChhyfFf0Oaxh0dRYoYaxyNpv31yCNnlZWD324kQnhahZYbBJAS9fPzDseWOCyZodZTi6cyu0zEnMR15VW4/szAzzDg8eq6iuT3QSiJqdxP/yybGYJRvD8/LqOox8fi5uf3+FW0mKqe+DU3DDa4sSdn4iiu31eT/gulfd+Z0y2KSgWG02xtdr6gIzCizfWtrkc5cdqsPgp2Ziw+5yx++dV7ivyedPFQu27MPot5bCb7KyKlGyGvvFOszeVOLKsRls0oTFOM+4m1uwF5tLKvHcjAJvTpiibnljCSav3Y2Dtc2zyu7FmYU49r6v437cTcUV2FRcEffjkvvYZpOC8rIzI7ZlZ2Y0LDcQ5b31Pj+mrSuGiOC4rq1xdOfWLqWSklHZoTq0zs1CZoa7dydPTN7oynEveno2AKBo/AhXjk/uYckmBV11Wo+Ibfo51KLVsr3w3Wb87p1lGP32Ulzw1Cw3kkdJyudX6PfQVNz78epEJ4WaIQabFGTWoysny95HuavsULyTQymi3h8o+X6wZJtn5+SYsNTkxufGYJMmzKa1IdJLRL7PWJOa6l3o2MJgk6L+fknfsOfJMH6Fkps/ATl/Is5JTVfvwjIlzKFSVCvD4EyzdW+I9BLRC9utc7JLubvq/PFfgJE5VIrKMPQmYsmGYkmnkk0tV6N1Rfd2LZCTmYEWJj1em4o5VIoyxpbsrIbg069HO8v32RmPU7S3Ej9/cR7Kq7kImx0Lt+xDcXl1QtNQWlmLLSUHo+6jEpA/uxXfGGzc0bVtLgb16uDKzSuDTYrKEOuSTe9OrZp07Kenb8KyrQcwY32x5T6sim9w1cQF+OljMxKahnOf/C5mV/Z0Ktnw++cOBfcGiHNQZ4rKygi/T9AHG/0PfNamEmS5PICPAnx+5fpgSSvlNiYXTa9gw2jjBr+KvJGNF5ZsUtRRHVuGPQ8b1Knbfv2ri/Drlxc26hwC6y+dV9PjJLt6i+qc4GeQTHmiz4XEbN1XhctfnIeyQ+ZVrm614yfTdU0nSinXftsMNinqxO75Yc/142ysfuDx/IHyxx6gbzsI3m2XHarDwRqtpJFE18mNz+y5bwuwfOsBTFm72+Kc7lwAdql2h2LJhmLRV6PV1vsSmJLk5UbGF5yPDmiIKyUViessEO1vTEw1mjvHZahxh18puFUTzGCTJvTT1UxZW4ydByKnpYnnDUsqVqO5kdeGBZvQ8RsujvI4W4w28jsx42xYskklge8ISzYUhXFQ5w97KxOUkuZF3w4SDCz6QOx1nlgXpUtwIgZCsjdaalEs2ZCZGXecG3ps7Bc/cfYWi3elYJEkTtzIn5JtIHtVrXUVajrNjcZg4w622ZCp3rq1aPSDOoFAl+eV2w5EvGfdzjLX05Wsyi16TDWFvo0k+FD/SXidJ/7u7aWWr7la9WRxaFajpRa/UshwKSow2KQJsxG/o16YF/ZcKWDldvvBxuv2Brfd/sGKuB/TLM8T3Z2h1+NBFheV4pBF6YYdBCgWv1JRhzw0BYNNmojnRJxN/aq9NHsLlm0tjUta4unHffFvx9Ln38lys71uV7npdleDjcWXxq12Ik7E6Q43ZxBgsDEQkUtEZGJZWWpVN8VrLiO/X2F/VdOqmx79ej1+/uL80PP/LfZusa5o3Mhr/WYdBPTn1P6fvGYXqmpjj/J3UyLyZ7bZpBa22XhIKfWFUurW/Pz82Dsnkbzs2B9lrDvbbfurcPS9X2P2ppJ4JQsA8PHy7XE9XmO5US3oN2mz0Z9FKWDNjjKMfnsZHvh0bdzPb87870yr6WpYkeYKP2cQoFjybEwJ7jO5tZ28pmHk96biCtvnS8U7SzfSrEweGzPYCm3esm2lVfFPgAP6JUri3pbkeQcBVw7b7LFkQ5b+cH5v3DOsj61gYzZZ4+govZfSjSvBJqxkoyK36XLhRHc6NyuFNVWsvykYFPx+hb0Ha+JzUnAiTrewZEOW7ry4D247tzc6t8kN2946190Jve18IV+eYzXWJ32EdRDQ/vcbijteV/lY5cNmaW3yuWKmJbDHP6duxMBx01FSEZ+Aw5KNO5SKPgFvUzDYpIkj2rcIe94ix/lKe05uFu3s++z0AsdpcJMbd8P6TC94+IiqoyTJGPXpinv1llVvNO0007W1kfZVNj7YmJUiKb44gwDFZFzG1U6HgaCPlpo34Df195xs2YE7MwhEHjWsbUS3PdHzyenTWrgn+qqe8T5nw4DXxl8Evwsls0Sq8/lx96RV2J7gtjw9rmdDMbXJyw577mQN8Ts+XIkNu8sd/YDtfB+jzdOVLsJijUnJRv+6W9UTRlafjT5dw56dE5dzTbK4UTE7J9C0gKvv4JIOMwgs2LIPHyzZhrs/WpXopISwzYZiMrbZOAk2ALB9f+Qs0dHY+a07yQ+umDAfA8dNc5QGp7waZxPeNuJ9pmj1d7razhGjnSgep3baweEP7y7DZyt2xOHM7vLqJsQOvwqfASOeGGzSlJ3eaXpufL/q/fZLNkt/LMXeg7XxT4SO2xm/VZtNstyDJ8Osz035mjltc/pq1S7c/n78pyhKb2yzIYdWbo+chDOaDBFHja52ljAI5m3JUuPhRl4bXrIx2aYa1phJfJtN4s4Z/G7FqxotWb5T6YZtNuRYdZ2z9hKzAZ9rdpjPsQUAG3bbHwCaLNypRtMfX6tG071eVVsPn4MSXjxU1JhPi5OIHlxl2kzb8RhtpL+M5dXxn8E7UZJpNgS22ZDrfCryK//qvB8SkpZUokxKNvptt761FPW+wPNMt+onDG58bbHpdldLNhZ/2vWvLgo8MJx7wZZ9OGgRFAHgmembsGZH+PyE+hLjda8salQyk1HR3uTpjcYZBMix9i2zY++k4/cr16omFv6w350DO+bVOJuGbVtKKkOlRrcaXu1KxHo2RkP+NQsrth3A1RMX4M/vL7fc75npBRj5/NywbfpVUaMtf51qdpgs4Z4oLNmQY89dM8DR/la/3bkFe6O+b09FNWZu3BN6nsyD7dxJWmRvNGNDfJ32PDPhbTaJ+2z0Z172Y2D5ifW7zKtirb5D6dDdOdlxBgFyrFt+i9g76Vj9kJ+csiHs+Xe6wFJQXIGfvzgfN+iqbeKdHxyq9VkuBuaUO4M6I09gDNzBNpvvNsZ3Nm2nEhpsdOeOdedsdePjcdOX67Y5HG7gBc4gQI5lOfzGBDIik1+5IWfQD+K78OnZ2F4a+MEEM5OPl8d3XMPxD07GgEemxuVYdkpdu8uq8db3RbZLaGbdiY3vDbbZJJqbmfVdMQYmml0Bq2tstd2XZiWbez9ZnegkRPArIMOlaMNgk6acNkZb3fWu3HYAB6pij38Jvv21OHYqKNMWcXPas86KnWr+OyetxAOfrcWZ47+1dUyTgo1JySY5MsnElmwaHtudKTpie5Jcx3QWWBbaHQw2acppxvKXD1Zi0lLzUsktby6xfT5jiaq6rvFVYEOentXo95qxU1oJTrGzs6za1jHNRrUbr31dkmSSxmQUOFi/qKm27m/ocRXsKGF1Vay6Aqdzm01NfXyqiptKgTMIkAOz7zzf8XQ1QMPMvEYFNiZtDGYDK7eHd1cNVrM1Rrymo3fCaYmwtLJhvEeog4AhU6xpQsCNJ2Ow3VzizWSctfXhJVPj5JxGVtuTpYTohgc+XZPoJABgmw3Z9NbNg/DGTYNwZMeW6NI2D1P/cg5OObKd6b79ephvbyyru06PhpbYYierWuSwm/Z+XRVjaB4ww4nMFq0zs7nkIFZuczbzgxORebU3H86/vw1faiI0q4DFJxKPud2SuVekmWQZHuDmDALurrBFnjr72M5hz4/r2gYFxeZ3r23z4vvRW/227Xxxd5U1lH6Wby3FgCPbxytZ4WzkP3UOG/NNB3UaTmSnzQsABj8VqDacNPoMDOzZwVE67Ii8IfAmQ15cVBr2PDg5pnVQaXo1WorFmqSZipPjbKjxLL44j152ku1D2PnhWmUExq3BRn+9LSUN86zNWL8n4nUAcen+7ErXZ93t9rqdgel9dpeFV/8dMPmbo/nFf75vcrqumDA/Ylui2jyMmdcqQ1WrkT6V+nYlYzVaZZQZCFKtfSce7STl1XV4cWZhkzpSKMU2G2okqy7Q+Q5mGLCzLo3V9/uRL9eFbzBJThtdKcuq2u2TOHSpNq70+NmKHVi/q9z0dbv0f3ewI4WxS2u0TNGKk1HlxeWRnRmW/lgasc2YAa/eET3TjxervMvqauvTeeHTs023A4Epb6ykVqgJTGxb0cT53sZ+vhZPTN6IWZsaN54r+P1nmw01ilWjqpOMtabeTrAxP963GwwlFZPdsjIavoZWd1XxuNnSn9rnV7j9/RVhi4i9Mb/I8TGNf7e+SjCouhE9jR79KjxI1/n8GPXCPNzy5pKIVTbtNpwbx9m88N1mx+nygmX1miH90UqMqVayAYCbXjef086uCq1t0M7v1Uzwa8QZBKhRjI3TwTnTnJS07WRmyub322xgnj5jaMoa9bHoT/2/JZErTC4qct5Ia/xzzAZwNmac0Nerd4c9Ly6vxsptBzBtXTEuf3Fe2Gt2B/DGOwO2e8NilXlZ90azN6hzToH1HXwKxppQ21adz4+npm50XNJpaohgyYbiZsE9g0MN9nHPeGxWXBgD176DNXh2RkNvpbcXbI1ruvT0aVy2NbKaybjgnJ3M1M6yx/Fob8rJavipVhh7t9nMHOKdAdu9YbEulTrtjRb+wqcrdtpLQIr5fMVOPP9tIf45ZaOj9zX14w1+npxBgBqlv9bFeeWDF+Gw/DzcM/x4ZAjQNi87rt2f7WY8xgzj7o9WYdo68/E98aY/9SaTAY3GYBNtPZ8g499tllGu07UL2Wn/MmU4rr66zu40NPG+wXCyEqsZq+Ss321+3Z00fOv/VuM4n2QQ7ExiZv7mQFvUoUaOz2pslbPbVY8MNmnuw9FnYPXYi0IdAn5x6hHY8o8RyMnKwHu3/DRu57H7RTWWbEpt9tSavamkyWMn9O826xGVafiV2vmxG0t0saa+n1cYfRZtPX1PrDe//zHstd26GQ4sewIatsd7TKTdWGPVDrev0rxL+OvziszP52icjf59yVen9s+p5qWWPeXV+GhZoIrX62S7fT4GmzSXnZmBNnnmPc9a5sRvrE1jg43dO/1v1uzGu4vcq2IzY6cXmfHPLoqxXLaTUfD6YPfv7wrDXtNP/2IdbMKf2/2M9h6swf8Wb4u5n92JMTdalFQA81KH1VGdXDv935osAyb1rGqqBj02I/TY6xAZbAN7bkZBjD0bh8GmmSsaPwJF40eYvnb74GPDnkcrWdi9K7r8xflh80CZDaK0CkBzNtkvFZiySGOdz48lRfvx1oLw0sNDX6yNOWWOsWrnxhg9ipzcnQcb1leblMJuf3+FLg3m7zcGA7slw9+/swx3fbQKW/dFX0HSZ/jsrObBKy63vob6KsYgp4M9zej3bEwvQ7dNtxhPpue0pDG7kV2eg4Ltgsd1bdOk41hhsCEAwGOXRw7ybJkT3oYRLaN0ctf8mq6apN4ksDz29XrT905eu9t0u11WnRgqqusx16R6q2hfVcw5q5yOVXGSYQZrn/YejBHwrHpvGT4wuyWDYICti1FPZgxmwbE9Tqo7b3vLbJLX+M4gkGxzqm3cbW8C1K9X73LUuaSxXZ6DWucGajruuOi4Jh3HSrMKNiJytIi8IiKTEp2WZDOoV+QUMfqpZuYW7MX2Uus7XSej5Et107eYtXGsjdJ42hRWec6q7QfwzHTzqoNYY2SmOuzc4CQjXqx1xY5V1WiVCRszH7O/f02UYBmrndlYBRbsJeekq7dZqcfJRJylFu0++uvc1I4M8VZZa2+Q76E6Hx78zLsJOoO/xZxMd8KCq8FGRNqJyCQR2SAi60XkjEYe51UR2SMiEVdeRIaKyEYRKRSRMdGOo5TaopS6uTFpSHdtW4S36ww94TB0zc8LPX930Y8498mZlu930k1TH8TMMtLsKOsnN6WTgNV79SuNGjm9K7aa+LThePaP9dAX62ylwerl0W8tNewXuePz3za+fv6JyeGruH6/OVA6rLKZmdoVLP2apf+6VxeZvke/a3WdP6nWwjF2RImmaF/0NsB4Cl6jVO36/CyAyUqpPgD6AQirHxGRLiLSxrDtGJPjvA5gqHGjiGQCeAHAMAB9AVwjIn1F5CQR+dLwr0t8/qT01KlVbtjzF399Ckae1C303DjIEACe/MXJocczNuyJ+EH/rHdH03NNmNkwct1sEOS8QutpSF73uP7d6SqbsXZfuf0Ayh0O1ovVw804o0DQ91v2oXBPQ5WNWaxtygy/xkGwb2g95qritIx30H9nbwFg3ja1ekeZaVuFPjAt/bEU93ycPKtiOlnGwunEsABw0OYs40bB75nTVX7tci3YiEg+gHMAvAIASqlapZRx/vRzAXwqIrnae24B8LzxWEqp2QDMupQMAlColVhqAbwPYJRSarVSaqThX+wWuUAaLhGRiWVl3swblSyMdzMZGYKMDEH3di0s39O+ZU7Y8xsMjeO3nH205Xu/Xr0LALDbZF6vaB76Yh0ONmKuMaBxXX+dnqv8UPRAMnH2Ftz25tKo++gppUznPtO7L8rywvpSm1nJwFj96fMr/KD1qIt1vazGrzgdH7LbsFCd8bRPaqVmq95vZqUb454fLIndu84rTgK83ZK1/objjg9Xht1k2D6XSt2STS8AJQBeE5HlIvKyiLTS76CU+hDAFAAfiMivAdwE4EoH5+gOQP8t2q5tMyUiHUXkPwAGiMg9Zvsopb5QSt2an5/vIBnpYdad5+G3Z/XCZ384M7Rt+l/PNd33+jOOwhmGkovxDvP8PtaFyd+/s8y0l1VQtLaEq/77PbZ4tPDXfos2ATM9O7YMZdRBZneJ35tMIPm2oSdc0N8+XGUZwIIlyWjZkX7xumC+dcPPeoalRT+h5YSZDV2sYwVasyrQep/fcclmvaFHWpnJ3+vzKzz0+Vrbx0zGsTVBTozr/AIAABJcSURBVNK2ekcZ5hbsjToIFIi8hiu2Ob9Z9qdqyQaBtXJOATBBKTUAQCWAiDYVpdQTAKoBTABwqVLKtVxEKbVPKTVaKdVbKfUPt86Tqo7q2Ar3j+wbNrNAi5xMFD46DIfr2m+evqofxl56AlrlZll2m7ay6N7BocfvLjLPYAFg5PNzw57rqx7W7izHBU/Nwp/eW266ZEE8BE9nd/blj3//Mww9sVvE9r9f0td0/55jvgpL+/0Wvd4+WrYdz31baPraQ18EMl99W9T/botsFl27M5DxHKqthwjwwMi+eGBkQ7qunrgg9Fg/2PWyF+ahaG+l5YJuZlWMK7eXOW6zufH1xWFB3WzG6uLy6tBS3Sd2bxvxurHbdXVtZCB067vilNNAeO0rCzH8uTnYE6WEaywt/e3DlTF7MRoFq9HcWjzNzWCzHcB2pdRC7fkkBIJPGBE5G8CJAD4B8HeH59gBoIfu+RHaNoqjrMwMPHLZiTimS2ssuncwLh9wRNio8OUPXIinruwX9p5zjgss5HbucQ0LuvXUVhANem9ReNXGkOOtS0Lf3H52xLbPV+5Ev4enYuPuCiz9cT/mFkR2X67z+aGUctyxYN3DDU2EPcd8hZ5jvsLAcdNxw2uLQnf0+jaqU45sj+8MM1zfPvhYZOl69uQbOmH0e3hqkxrT3/j+R5RU1IRmYbjujKMwqFcHLLp3ML78v7NC+414bi4mzNyMXWXV6NgqF5kZgpvP6oUhx3cN7TP6raX4dkNxRO+68/45E6NemGd6/Sp0JZ9u2s3IFRPmh612asy39J/x787rHXp8yiPTcMFTM01LNUB4d/iHLj0Bc+46P+z1P7yzDNe/uigUWM16fPV7eGqjq2DjKVYbnBX9gE8js/gwcNx0R8cPfp+dLo1ul2vBRim1G8A2EfmJtmkwgLB500VkAICJAEYBuBFARxEZ5+A0iwEcKyK9RCQHwNUAPm9y4inC4OO7Yvpfzw0LFkHtW+XgilOPwMvXDQxte/aq/gCAidedin5HBKokgz+x284xb8t5+frT8MM/hkds/+b2s3Fc1zZh1T96Fz8zG1dM+B7XvrIQN762CAPHTcdJY6fgrkkrcex93+CcJ7/Db14Jr9d/4oqTTY8FAP/+1YCIedKAwHiXmRtLMPzZOXh17g/YZKgX79MtfDDcXy48Dr06BWqOxwzrg7l3n49v7wivlpy+fg+27Y/sUv7H8836yQCX9jscyx64MPT8tEcbMpRgFtGlbR5O7B5eDfz45A34cOl2dGnT0BHk0ctPDD2evHY3bnrdbMxLwIGqOkxbV4yeY77Cdxsjmz9f0n32wW7kb940CBseCe/XM+HaUzHzb+ehaPwI3D20T1gpZUtJJYY9Mzts/4VaSfjLVbtC26rr/OjRoSUW3zcklDHO2LAHszaVYMRzgRKx1ewPMzfuwaFaH+p8fpRU1ODtBT+Gd5P2+fHG/KKwoOT3K9zz8eqYVVl2GdthrhrYw2LPSD3HfGVaxdzYAGZ2DLeq0cTNtbpFpD+AlwHkANgC4EalVKnu9TMBlCulVmvPswHcoJR6yXCc9wCcB6ATgGIAf1dKvaK9NhzAMwAyAbyqlHo0HmkfOHCgWrLE+sdH5uYV7sXx3dqiQ6uGzgM19T4MfmoWHh51Ai7oE7ib/vXLC8J6nZ3YvS2+/L9A6aVwz0EM+VdgieRnruqPywY0NMP1HPNVXNK5adww5GRl4MWZhTiyQ0t8sHgb5hTsxas3DAylcUvJQTwzvQCfr4w9u3DR+BGorvPhzPHfhub8ClYxrttZjuO7tQmVBqtq69H3wSmWx2qZk4l1Dw/F45M3hPXcA4A/DzkWfx5yHGZtKsH1hobxuXefjyPat2z4G4srcNHT4Zl3325t8bWulFhZU48T/m6dFru2PDYcR9/7tWl67v90dWg2b7Nq12ifadH4ERGvBz87AKiorsNJY6eGvd6vRzus21kW6sl1QZ8ukesqGfxlyHG4fEB3XPPSglDV6bATD8M3a3ajU+sc7D1Yi8Pz8zD/nsFh75u6djeenLIR/7yyH3p1boVnphVAQeGmM3uh5GANTjw8H6VVtWjXMhu5WYEbmBHPzQkbS7b5seHorbt2OZkZ+NPgY/Dpip2WPQ3fvGkQDsvPQ9u8bNTU+3DbW0uxwWSw6J+HHBsK/mOG9cHoc3ujtLIWT0zZiPtHHI+87EyM/2Y9fnN6TyzfVorb31+BGXeci96dW0e9XtGIyFKl1MCI7W4Gm1TGYOO+dTvL8Yd3l+GaQT3w27OOjugFs21/FXp0aBm2rbKmHhkiOP7ByU06tzHT8/sVKmvrTeeR8/tVREaq17NjS8y8s6Fa552FP+Lw/BZRO0i8Mb8Ifzc0eP+kaxuUHarDV386Cx1bB0ogz04vwNPTNwEAfntWL9w59CehTGvguOmhevlL+h2O568ZEHGe6jof+jzQcK3a5GVh9diLY6YFCASmESd3C/UGs/Lbs3rh/pF9I4LCsgcuRIdWOVBK4ckpGzHy5MPR9/DI9pbhz84xnbbm2av7Y1T/7vjtG4vDpncxfnZzCkoiSq56H/3uZyipqMHot+33AvTKigcvRLuWOSipqIFSCoMem4HObXKx+L4hAIAPl2zDnZNW2T7eAyP7orSyNmIuPTOd2+Ti6E6tIuaO+/6eC9At37oXaiwMNg4x2CS3PRXVyMnMQLuWOVr34Boclp+HOp8f475cFxrzYcVpx4Z5hXvx+OQNOFBVFzYJJgAUPjosrG3GrkU/7Mcv//t96PnsO8/HkR1bRnlHuDqtymf4Sd1weJQu6tV1Pjz85Tq8u3ArHhjZFzef1SvsdaUUnplegGdnFODq03rg2tOPwjFdWoeqEsd/swH/nb3ZcmT/E1ecjF+e1gM19T68OrcIj2uDPX/4x3Bb69nX1PtwqNaH/g9PAwC0ysnEWl2b2YGqWlwxYT42l1Ri47ihoWCrt21/FdbuLA8LKN3y8/DSdQNDVYrPzSjAv6ZtipkeLxm/hxNmbsaFfbvimC6BksW+gzU4ddx03HxWL7wy94eYx9vwyFDkZWdixbYDuOyFeTH3N7N67EWWk/fawWDjEINNaltctB+vzfsBz19zCjIzBAXFFWiVmxU1U7ajcE8Fxn21Hr8/7xg8PnkD3r75p2iRE5n52bVuZzk27C7HhX27NukHbofPr6I2/s7cuAdnHtMJ2RaBs6C4Aut2laNT61zM3LgH63dV4MqBR+DSfofbCiqxjP18LV6fXxQqETVWZU09lm89gH498iOuadmhOtz0+mJ0aJUTWkfpVz89Eu8ubJhR/NzjOmNWEye1jMVYGrZjTkEJ8ltk49qXF0aswAsEFkc8TNdrdMGWfWE9De3a8tjwJo21YbBxiMGGKL3NL9yL3OxMnHpU+7BOAiKC3WXVuPWtJXj6qv44PL8FtpVW4cuVO3Fqzw6oqqnHpuKDuH1IYFZ0n19hV9khvLtwK0af1xsLNu/DrW8txeUDuuOpK/vh+W8L4VMKp/fqgNZ5WZi6thi/P793k5b42F9Zi/Yts7F6Rxnu+N9KPHr5SRjUq0PEfvU+P5b+WIptpYfws94d8Z9Zm9G1bR4uPuGwULvoQ5eeEKpGfe+W0yPGzznFYOMQgw0RNdb20ip0aZMXtpx3svlq1S7kZGXgwr5dY+/sgFWwid/qWUREBABhvQKT1YiTIwchuyl5wy4REaUNBhsiInIdgw0REbmOwYaIiFzHYENERK5jsCEiItcx2BARkesYbIiIyHWcQcCCiJQAiD6bo7VOACJX8ko8pssZpssZpsuZdE3XUUqpzsaNDDYuEJElZtM1JBrT5QzT5QzT5UxzSxer0YiIyHUMNkRE5DoGG3dMTHQCLDBdzjBdzjBdzjSrdLHNhoiIXMeSDRERuY7BhoiIXMdgE0ciMlRENopIoYiM8fjcPUTkOxFZJyJrReR2bftYEdkhIiu0f8N177lHS+tGEbnYxbQVichq7fxLtG0dRGSaiBRo/7fXtouIPKela5WInOJSmn6iuyYrRKRcRP6cqOslIq+KyB4RWaPb5vgaicj12v4FInK9S+l6UkQ2aOf+RETaadt7isgh3bX7j+49p2rfgUIt7Y1f5N46XY4/u3j/Zi3S9YEuTUUiskLb7uX1ssofvPuOKaX4Lw7/AGQC2AzgaAA5AFYC6Ovh+bsBOEV73AbAJgB9AYwF8DeT/ftqacwF0EtLe6ZLaSsC0Mmw7QkAY7THYwA8rj0eDuAbAALgdAALPfrsdgM4KlHXC8A5AE4BsKax1whABwBbtP/ba4/bu5CuiwBkaY8f16Wrp34/w3EWaWkVLe3DXEiXo8/Ojd+sWboMrz8F4MEEXC+r/MGz7xhLNvEzCEChUmqLUqoWwPsARnl1cqXULqXUMu1xBYD1ALpHecsoAO8rpWqUUj8AKETgb/DKKABvaI/fAHCZbvubKmABgHYi4vb6tYMBbFZKRZsxwtXrpZSaDWC/yTmdXKOLAUxTSu1XSpUCmAZgaLzTpZSaqpSq154uAHBEtGNoaWurlFqgAjnWm7q/JW7pisLqs4v7bzZaurTSyS8BvBftGC5dL6v8wbPvGINN/HQHsE33fDuiZ/auEZGeAAYAWKht+qNWFH41WEyGt+lVAKaKyFIRuVXb1lUptUt7vBtA1wSkK+hqhGcAib5eQU6vUSLSeBMCd8BBvURkuYjMEpGztW3dtbR4kS4nn53X1+tsAMVKqQLdNs+vlyF/8Ow7xmCTZkSkNYCPAPxZKVUOYAKA3gD6A9iFQDHea2cppU4BMAzAH0TkHP2L2t1bQvrgi0gOgEsBfKhtSobrFSGR18iKiNwHoB7AO9qmXQCOVEoNAPBXAO+KSFsPk5SUn53ONQi/qfH8epnkDyFuf8cYbOJnB4AeuudHaNs8IyLZCHyR3lFKfQwASqlipZRPKeUH8BIaqn48S69Saof2/x4An2hpKA5Wj2n/7/E6XZphAJYppYq1NCb8euk4vUaepVFEbgAwEsCvtUwKWjXVPu3xUgTaQ47T0qCvanMlXY347Ly8XlkAfg7gA116Pb1eZvkDPPyOMdjEz2IAx4pIL+1u+WoAn3t1cq0++BUA65VS/9Jt17d3XA4g2EvmcwBXi0iuiPQCcCwCjZLxTlcrEWkTfIxA4/Ia7fzBnizXA/hMl67rtN4wpwMo0xXz3RB2t5no62Xg9BpNAXCRiLTXqpAu0rbFlYgMBXAXgEuVUlW67Z1FJFN7fDQC12iLlrZyETld+55ep/tb4pkup5+dl7/ZIQA2KKVC1WNeXi+r/AFefsea0sOB/yJ6fAxHoJfHZgD3eXzusxAoAq8CsEL7NxzAWwBWa9s/B9BN9577tLRuRBN7u0RJ19EI9PJZCWBt8LoA6AhgBoACANMBdNC2C4AXtHStBjDQxWvWCsA+APm6bQm5XggEvF0A6hCoB7+5MdcIgTaUQu3fjS6lqxCBevvg9+w/2r5XaJ/xCgDLAFyiO85ABDL/zQD+DW32kjiny/FnF+/frFm6tO2vAxht2NfL62WVP3j2HeN0NURE5DpWoxERkesYbIiIyHUMNkRE5DoGGyIich2DDRERuY7BhshlIjJf+7+niPwqzse+1+xcRMmGXZ+JPCIi5yEwK/FIB+/JUg2TXpq9flAp1Toe6SNyE0s2RC4TkYPaw/EAzpbA2iV/EZFMCawNs1ibPPI2bf/zRGSOiHwOYJ227VNtItO1wclMRWQ8gBba8d7Rn0sb+f2kiKyRwLooV+mOPVNEJklgTZp3tNHlRK7KSnQCiJqRMdCVbLSgUaaUOk1EcgHME5Gp2r6nADhRBabEB4CblFL7RaQFgMUi8pFSaoyI/FEp1d/kXD9HYELKfgA6ae+Zrb02AMAJAHYCmAfgTABz4//nEjVgyYYocS5CYP6pFQhM994RgfmxAGCRLtAAwJ9EZCUC68f00O1n5SwA76nAxJTFAGYBOE137O0qMGHlCgQW8SJyFUs2RIkjAP5PKRU2kaHWtlNpeD4EwBlKqSoRmQkgrwnnrdE99oH5AHmAJRsi71QgsCRv0BQAv9OmfoeIHKfNjG2UD6BUCzR9EFimN6gu+H6DOQCu0tqFOiOwXLHbs1QTWeIdDZF3VgHwadVhrwN4FoEqrGVaI30JzJf/nQxgtIisR2DW4gW61yYCWCUiy5RSv9Zt/wTAGQjMtq0A3KWU2q0FKyLPseszERG5jtVoRETkOgYbIiJyHYMNERG5jsGGiIhcx2BDRESuY7AhIiLXMdgQEZHr/h/rXtZNspU0lAAAAABJRU5ErkJggg==\n","text/plain":["<Figure size 432x288 with 1 Axes>"]},"metadata":{"tags":[],"needs_background":"light"}}]},{"cell_type":"markdown","metadata":{"id":"NJOqjNGXH7HU"},"source":["# Model evaluation"]},{"cell_type":"code","metadata":{"colab":{"base_uri":"https://localhost:8080/"},"id":"phUd-6PfH87I","executionInfo":{"status":"ok","timestamp":1612980826195,"user_tz":-180,"elapsed":5516,"user":{"displayName":"Vl F","photoUrl":"","userId":"14485856366916106708"}},"outputId":"610ce970-1848-4750-eee7-553a6b4b0073"},"source":["for J in [3, 4, 5, 6, 7, 8, 9, 10, 11]:\n","  N_levels_ = min(N_levels, J-1)\n","\n","  spectral_radius = stochastic_trace_u_net_solver(A, J, N_levels_, projection_filter_pre, projection_filter_post, pre_smoothing_filters, post_smoothing_filters, u_net_strides, N_sweeps=10, batch_size=10).numpy()\n","  print(\"J = \", J, \" spectral radius = \", spectral_radius)"],"execution_count":7,"outputs":[{"output_type":"stream","text":["J =  3  spectral radius =  0.49641305\n","J =  4  spectral radius =  0.54704124\n","J =  5  spectral radius =  0.6452189\n","J =  6  spectral radius =  0.85153264\n","J =  7  spectral radius =  0.9834595\n","J =  8  spectral radius =  1.1335756\n","J =  9  spectral radius =  1.2149988\n","J =  10  spectral radius =  1.3108406\n","J =  11  spectral radius =  1.4078285\n"],"name":"stdout"}]}]}