diff --git a/example/example_jaxcnn.yaml b/example/example_jaxcnn.yaml new file mode 100644 index 00000000..90dcf89c --- /dev/null +++ b/example/example_jaxcnn.yaml @@ -0,0 +1,21 @@ +metrics: [SNR_3x2, FOM_3x2, FOM_DETF_3x2] +bands: riz +training_file: data/training.hdf5 +validation_file: data/validation.hdf5 +output_file: example/example_output_jax.txt +# Backend implementing the metrics, either: "firecrown" (default), "jax-cosmo" +metrics_impl: jax-cosmo + +run: + # This is a class name which will be looked up + JaxCNN: + {% for nbins in [2,3,4,5] %} + run_{{ nbins }}: + # This setting is sent to the classifier + bins: {{ nbins }} + # These special settings decide whether the + # color and error colums are passed to the classifier + # as well as the magnitudes + colors: True + errors: False + {% endfor %} diff --git a/example/example_jaxresnet.yaml b/example/example_jaxresnet.yaml new file mode 100644 index 00000000..2767014b --- /dev/null +++ b/example/example_jaxresnet.yaml @@ -0,0 +1,21 @@ +metrics: [SNR_3x2, FOM_3x2, FOM_DETF_3x2] +bands: riz +training_file: data_buzzard/training.hdf5 +validation_file: data_buzzard/validation.hdf5 +output_file: example/example_output_resnet_jax.txt +# Backend implementing the metrics, either: "firecrown" (default), "jax-cosmo" +metrics_impl: jax-cosmo + +run: + # This is a class name which will be looked up + JaxResNet: + {% for nbins in [2,3,4,5] %} + run_{{ nbins }}: + # This setting is sent to the classifier + bins: {{ nbins }} + # These special settings decide whether the + # color and error colums are passed to the classifier + # as well as the magnitudes + colors: True + errors: False + {% endfor %} diff --git a/notebooks/JaxCNN.ipynb b/notebooks/JaxCNN.ipynb new file mode 100644 index 00000000..e891b2f5 --- /dev/null +++ b/notebooks/JaxCNN.ipynb @@ -0,0 +1,303 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Found classifier IBandOnly\n", + "Found classifier JaxCNN\n", + "Found classifier RandomForest\n", + "Found classifier Random\n" + ] + } + ], + "source": [ + "import sys\n", + "sys.path.insert(0, '../')\n", + "\n", + "import numpy as onp\n", + "import matplotlib.pyplot as plt\n", + "%matplotlib inline\n", + "\n", + "import jax\n", + "import jax.numpy as jnp\n", + "jax.config.update(\"jax_enable_x64\", True)\n", + "\n", + "from flax import nn, optim\n", + "\n", + "import tomo_challenge as tc\n", + "from tomo_challenge import jax_metrics\n", + "from tomo_challenge import metrics\n", + "\n", + "import warnings\n", + "warnings.filterwarnings('ignore')" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [], + "source": [ + "training_file = \"../data/training.hdf5\"\n", + "validation_file = \"../data/validation.hdf5\"\n", + "bands = \"riz\"\n", + "colors = True\n", + "errors = True\n", + "array = True\n", + "\n", + "x_train = jnp.array(tc.load_data(training_file, bands, colors=colors, errors=errors, array=array))\n", + "z_train = jnp.array(tc.load_redshift(training_file))\n", + "x_test = jnp.array(tc.load_data(validation_file, bands, colors=colors, errors=errors, array=array))\n", + "z_test = jnp.array(tc.load_redshift(validation_file))\n", + "\n", + "assert len(x_train) == len(z_train)" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEGCAYAAAB7DNKzAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAVtklEQVR4nO3df4xd5Z3f8fcnwBKrWQIGQx2bxKlwpQJSk2AZKtSKhi1YoSqoAtWVuriVJWsRkbJqpcbsH/UuyJHpH6RLtqFCi4Whm4DFboqbDaUuDlrtihhMSpYAobgLCwYLe2MgoAKtrW//uI831+M7Z+6M59edeb+k0b33e85z5nl8kvlwznPOuakqJEkazyfmugOSpPnNoJAkdTIoJEmdDApJUieDQpLU6fS57sB0O++882rVqlVz3Q1JGinPPvvsX1XVskHLFlxQrFq1in379s11NyRppCT5y/GWeepJktTJoJAkdTIoJEmdDApJUieDQpLUaaigSPJakueTPJdkX6stTbI7ySvt9Zy+9W9Lsj/Jy0mu7atf1razP8ndSdLqZyZ5uNX3JlnV12ZD+x2vJNkwXQOXJA1nMkcU/7CqvlBVa9rnzcATVbUaeKJ9JsnFwHrgEmAd8O0kp7U29wCbgNXtZ12rbwTeqaqLgG8Cd7ZtLQW2AJcDa4Et/YEkSZp5p3Lq6XpgR3u/A7ihr/5QVX1cVa8C+4G1SZYDZ1XVU9V7tvkDY9oc39YjwNXtaONaYHdVHamqd4Dd/DJcJEmzYNigKOC/J3k2yaZWu6CqDgK01/NbfQXwRl/bA622or0fWz+hTVUdBd4Dzu3Y1gmSbEqyL8m+w4cPDzkkSdIwhr0z+8qqeivJ+cDuJD/rWDcDatVRn2qbXxaq7gXuBVizZo3fxDRHrty2hzff/fCk+oqzl/Bnm788Bz2SNB2GCoqqequ9HkryPXrzBW8nWV5VB9tppUNt9QPAhX3NVwJvtfrKAfX+NgeSnA58GjjS6leNafPksIPTzOgKhNe2XXdSfdXmP56NbkmaIRMGRZK/AXyiqt5v768Bbgd2ARuAbe310dZkF/CdJHcBn6E3af10VR1L8n6SK4C9wM3At/rabACeAm4E9lRVJXkc+EbfBPY1wG2nOmidmjff/XBgIEhamIY5orgA+F67kvV04DtV9d+SPAPsTLIReB24CaCqXkiyE3gROArcWlXH2rZuAe4HlgCPtR+A+4AHk+yndySxvm3rSJI7gGfaerdX1ZFTGK8kaZImDIqq+gvg7w6o/xy4epw2W4GtA+r7gEsH1D+iBc2AZduB7RP1U5I0M7wzW5LUacF9H4XmnxVnLxk4oe3VUNJoMCg048YLA6+GkkaDp54kSZ0MCklSJ4NCktTJoJAkdTIoJEmdDApJUieDQpLUyaCQJHUyKCRJnbwzW+Pq+t4JSYuHQaFx+b0TksBTT5KkCRgUkqROnnrSnPHx49JoMCg0Z3z8uDQaPPUkSepkUEiSOhkUkqROBoUkqZNBIUnqZFBIkjp5eax8ppOkTgaFfKaTpE6eepIkdTIoJEmdDApJUieDQpLUyaCQJHUyKCRJnYYOiiSnJfmfSb7fPi9NsjvJK+31nL51b0uyP8nLSa7tq1+W5Pm27O4kafUzkzzc6nuTrOprs6H9jleSbJiOQUuShjeZI4qvAS/1fd4MPFFVq4En2meSXAysBy4B1gHfTnJaa3MPsAlY3X7WtfpG4J2qugj4JnBn29ZSYAtwObAW2NIfSJKkmTfUDXdJVgLXAVuBf93K1wNXtfc7gCeBr7f6Q1X1MfBqkv3A2iSvAWdV1VNtmw8ANwCPtTa/3bb1CPB77WjjWmB3VR1pbXbTC5fvTmm0i9h4d1+Dd2BL6jbsndn/Afi3wK/21S6oqoMAVXUwyfmtvgL4Ud96B1rt/7X3Y+vH27zRtnU0yXvAuf31AW00Cd59LWmqJgyKJP8YOFRVzya5aohtZkCtOupTbdPfx030Tmnx2c9+doguaj7zu7Sl+WWYI4orgX+S5CvAJ4Gzkvxn4O0ky9vRxHLgUFv/AHBhX/uVwFutvnJAvb/NgSSnA58GjrT6VWPaPDm2g1V1L3AvwJo1a04KEo0Wv0tbml8mnMyuqtuqamVVraI3Sb2nqv4FsAs4fhXSBuDR9n4XsL5dyfR5epPWT7fTVO8nuaLNP9w8ps3xbd3YfkcBjwPXJDmnTWJf02qSpFlyKk+P3QbsTLIReB24CaCqXkiyE3gROArcWlXHWptbgPuBJfQmsR9r9fuAB9vE9xF6gURVHUlyB/BMW+/24xPbkqTZMamgqKonaad+qurnwNXjrLeV3hVSY+v7gEsH1D+iBc2AZduB7ZPp52Lmd0tImm5+H8UC49VNkqabj/CQJHUyKCRJnQwKSVIng0KS1MnJ7BHl1U2SZotBMaK8uknSbPHUkySpk0EhSepkUEiSOhkUkqROBoUkqZNBIUnq5OWx85z3S0iaawbFPOf9EpLmmqeeJEmdDApJUidPPc0TzkVImq8MinnCuQhJ85WnniRJnQwKSVIng0KS1Mk5Co2MFWcvYdXmPx5Y/7PNX56DHkmLg0GhkTFeGAwKD0nTx1NPkqROBoUkqZNBIUnqZFBIkjoZFJKkTgaFJKmTQSFJ6mRQSJI6GRSSpE4TBkWSTyZ5OslPkryQ5HdafWmS3Uleaa/n9LW5Lcn+JC8nubavflmS59uyu5Ok1c9M8nCr702yqq/NhvY7XkmyYToHL0ma2DCP8PgY+HJVfZDkDOBPkzwG/FPgiaralmQzsBn4epKLgfXAJcBngP+R5G9X1THgHmAT8CPgB8A64DFgI/BOVV2UZD1wJ/DPkiwFtgBrgAKeTbKrqt6Ztn+BWeYXFEkaNRMGRVUV8EH7eEb7KeB64KpW3wE8CXy91R+qqo+BV5PsB9YmeQ04q6qeAkjyAHADvaC4Hvjttq1HgN9rRxvXArur6khrs5teuHx3qgOea35BkaRRM9QcRZLTkjwHHKL3h3svcEFVHQRor+e31VcAb/Q1P9BqK9r7sfUT2lTVUeA94NyObY3t36Yk+5LsO3z48DBDkiQNaaigqKpjVfUFYCW9o4NLO1bPoE101Kfapr9/91bVmqpas2zZso6uSZIma1JXPVXVu/ROMa0D3k6yHKC9HmqrHQAu7Gu2Enir1VcOqJ/QJsnpwKeBIx3bkiTNkmGuelqW5Oz2fgnwa8DPgF3A8auQNgCPtve7gPXtSqbPA6uBp9vpqfeTXNHmH24e0+b4tm4E9rS5kceBa5Kc066quqbVJEmzZJirnpYDO5KcRi9YdlbV95M8BexMshF4HbgJoKpeSLITeBE4CtzarngCuAW4H1hCbxL7sVa/D3iwTXwfoXfVFFV1JMkdwDNtvduPT2xLx/nNd9LMGuaqpz8Hvjig/nPg6nHabAW2DqjvA06a36iqj2hBM2DZdmD7RP3U4uU330kzyzuzJUmdDApJUieDQpLUyaCQJHUyKCRJnYa5PFZT4MP/JC0UBsUM8eF/khYKTz1JkjoZFJKkTgaFJKmTQSFJ6uRkthYsHxYoTQ+DQguWDwuUpoenniRJnQwKSVIng0KS1MmgkCR1MigkSZ0MCklSJ4NCktTJoJAkdTIoJEmdDApJUieDQpLUyaCQJHUyKCRJnQwKSVInHzN+iq7ctoc33/3wpPqKs5fMQW80jPG+p+L4Mr+rQjqRQXGK3nz3Q17bdt1cd0OT0BUEfleFdDJPPUmSOhkUkqROBoUkqZNBIUnqNGFQJLkwyQ+TvJTkhSRfa/WlSXYneaW9ntPX5rYk+5O8nOTavvplSZ5vy+5OklY/M8nDrb43yaq+Nhva73glyYbpHLwkaWLDHFEcBf5NVf0d4Arg1iQXA5uBJ6pqNfBE+0xbth64BFgHfDvJaW1b9wCbgNXtZ12rbwTeqaqLgG8Cd7ZtLQW2AJcDa4Et/YEkSZp5EwZFVR2sqh+39+8DLwErgOuBHW21HcAN7f31wENV9XFVvQrsB9YmWQ6cVVVPVVUBD4xpc3xbjwBXt6ONa4HdVXWkqt4BdvPLcJEkzYJJzVG0U0JfBPYCF1TVQeiFCXB+W20F8EZfswOttqK9H1s/oU1VHQXeA87t2NbYfm1Ksi/JvsOHD09mSJKkCQwdFEk+Bfwh8JtV9YuuVQfUqqM+1Ta/LFTdW1VrqmrNsmXLOromSZqsoe7MTnIGvZD4g6r6o1Z+O8nyqjrYTisdavUDwIV9zVcCb7X6ygH1/jYHkpwOfBo40upXjWnz5FAjk6ZgvMd7+GgPLWYTBkWbK7gPeKmq7upbtAvYAGxrr4/21b+T5C7gM/QmrZ+uqmNJ3k9yBb1TVzcD3xqzraeAG4E9VVVJHge+0TeBfQ1w25RHK01gvDDw0R5azIY5orgS+HXg+STPtdpv0QuInUk2Aq8DNwFU1QtJdgIv0rti6taqOtba3QLcDywBHms/0AuiB5Psp3cksb5t60iSO4Bn2nq3V9WRKY5VkjQFEwZFVf0pg+cKAK4ep81WYOuA+j7g0gH1j2hBM2DZdmD7RP2UJM0M78yWJHUyKCRJnQwKSVIng0KS1MmgkCR1MigkSZ0MCklSJ4NCktTJoJAkdTIoJEmdDApJUieDQpLUyaCQJHUa6ouLBFdu28Ob7354Un3F2UvmoDeSNHsMiiG9+e6HvLbturnuhiTNOoNCGoJfkarFzKCQhuBXpGoxczJbktTJoJAkdfLUk3QKnLvQYmBQSKfAuQstBp56kiR1MigkSZ0MCklSJ4NCktTJoJAkdfKqJ2kGeNmsFhKDQpoBXjarhcSgkGaRRxoaRQaFNIs80tAocjJbktTJoJAkdZowKJJsT3IoyU/7akuT7E7ySns9p2/ZbUn2J3k5ybV99cuSPN+W3Z0krX5mkodbfW+SVX1tNrTf8UqSDdM1aEnS8IY5orgfWDemthl4oqpWA0+0zyS5GFgPXNLafDvJaa3NPcAmYHX7Ob7NjcA7VXUR8E3gzratpcAW4HJgLbClP5AkSbNjwqCoqj8BjowpXw/saO93ADf01R+qqo+r6lVgP7A2yXLgrKp6qqoKeGBMm+PbegS4uh1tXAvsrqojVfUOsJuTA0uSNMOmOkdxQVUdBGiv57f6CuCNvvUOtNqK9n5s/YQ2VXUUeA84t2NbJ0myKcm+JPsOHz48xSFJkgaZ7snsDKhVR32qbU4sVt1bVWuqas2yZcuG6qgkaThTvY/i7STLq+pgO610qNUPABf2rbcSeKvVVw6o97c5kOR04NP0TnUdAK4a0+bJKfZXmte8EU/z2VSDYhewAdjWXh/tq38nyV3AZ+hNWj9dVceSvJ/kCmAvcDPwrTHbegq4EdhTVZXkceAbfRPY1wC3TbG/0rzmjXiazyYMiiTfpfdf9uclOUDvSqRtwM4kG4HXgZsAquqFJDuBF4GjwK1Vdaxt6hZ6V1AtAR5rPwD3AQ8m2U/vSGJ929aRJHcAz7T1bq+qsZPqkqQZNmFQVNU/H2fR1eOsvxXYOqC+D7h0QP0jWtAMWLYd2D5RHyVJM8c7syVJnQwKSVInnx4rzWNeDaX5wKCQ5jGvhtJ84KknSVIng0KS1MlTT9IIcu5Cs8mgkEaQcxeaTQbFGFdu28Ob7354Un3F2UvmoDeSNPcMijHefPdDXtt23Vx3Q5LmDSezJUmdDApJUidPPUkLiFdDaSYYFNIC4tVQmgkGhbSIdV3l5xGIjjMopEVsvKv8PAJRP4NCWgS65i6kiRgU0iLgaSSdCi+PlSR18ohC0knGO1V1fJlHKIuLQSHpJF1B4ET34uOpJ0lSJ48oJE2Kd38vPgaFpEkZLwyu3LbHAFmgDApJ08LHhyxcBoWkGeWpqtFnUEiaUR5pjD6vepIkdTIoJEmdPPUkaU503f092e0MOr3lI9Snj0EhaU5M1x/r8cJmvEeoexnv5BkUkkbaZB+h7n0gk2dQSBpp0/VHfLIB0mWhhctIBEWSdcDvAqcBv19V2+a4S5IWian8wZ9KuAwyXwJn3gdFktOA/wj8I+AA8EySXVX14tz2TJIGm64/7vPldNi8DwpgLbC/qv4CIMlDwPWAQSFpQZsv8ymjEBQrgDf6Ph8ALu9fIckmYFP7+EGSlyf5O84D/uqvt3fnFHo5f5wwlhG2UMYBC2csC2UcsHDGcsI4/hLIbVPe1ufGWzAKQZEBtTrhQ9W9wL1T/gXJvqpaM9X288lCGctCGQcsnLEslHHAwhnLbI1jFO7MPgBc2Pd5JfDWHPVFkhadUQiKZ4DVST6f5FeA9cCuOe6TJC0a8/7UU1UdTfJV4HF6l8dur6oXpvnXTPm01Ty0UMayUMYBC2csC2UcsHDGMivjSFVNvJYkadEahVNPkqQ5ZFBIkjotqqBIsi7Jy0n2J9k8YHmS3N2W/3mSL81FPycyxDiuSvJekufaz7+bi35OJMn2JIeS/HSc5SOxP2CosYzKPrkwyQ+TvJTkhSRfG7DOSOyXIccy7/dLkk8meTrJT9o4fmfAOjO7T6pqUfzQmwj/38DfAn4F+Alw8Zh1vgI8Ru/ejSuAvXPd7ymO4yrg+3Pd1yHG8g+ALwE/HWf5vN8fkxjLqOyT5cCX2vtfBf7XKP7/ZBJjmff7pf07f6q9PwPYC1wxm/tkMR1R/PWjQKrq/wLHHwXS73rgger5EXB2kuWz3dEJDDOOkVBVfwIc6VhlFPYHMNRYRkJVHayqH7f37wMv0Xs6Qr+R2C9DjmXea//OH7SPZ7SfsVchzeg+WUxBMehRIGP/RzPMOnNt2D7+vXao+liSS2ana9NuFPbHZIzUPkmyCvgivf+C7Tdy+6VjLDAC+yXJaUmeAw4Bu6tqVvfJvL+PYhpN+CiQIdeZa8P08cfA56rqgyRfAf4LsHrGezb9RmF/DGuk9kmSTwF/CPxmVf1i7OIBTebtfplgLCOxX6rqGPCFJGcD30tyaVX1z4fN6D5ZTEcUwzwKZBQeFzJhH6vqF8cPVavqB8AZSc6bvS5Om1HYH0MZpX2S5Ax6f1j/oKr+aMAqI7NfJhrLKO0XgKp6F3gSWDdm0Yzuk8UUFMM8CmQXcHO7guAK4L2qOjjbHZ3AhONI8jeTpL1fS28//3zWe3rqRmF/DGVU9knr433AS1V11zirjcR+GWYso7BfkixrRxIkWQL8GvCzMavN6D5ZNKeeapxHgST5jbb8PwE/oHf1wH7g/wD/aq76O54hx3EjcEuSo8CHwPpql0bMJ0m+S++qk/OSHAC20JuoG5n9cdwQYxmJfQJcCfw68Hw7Jw7wW8BnYeT2yzBjGYX9shzYkd6XuH0C2FlV35/Nv10+wkOS1GkxnXqSJE2BQSFJ6mRQSJI6GRSSpE4GhSSpk0EhSepkUEiSOhkU0ixI8ht933nwapIfznWfpGF5w500i9qzh/YA/76q/utc90cahkcU0uz6XWCPIaFRsmie9STNtST/Evgc8NU57oo0KZ56kmZBksuAHcDfr6p35ro/0mR46kmaHV8FlgI/bBPavz/XHZKG5RGFJKmTRxSSpE4GhSSpk0EhSepkUEiSOhkUkqROBoUkqZNBIUnq9P8BuJuKo+zTtfwAAAAASUVORK5CYII=\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "# Plotting redshift distribution\n", + "plt.hist(z_train, bins=50, histtype='step')\n", + "plt.xlabel('z');" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [], + "source": [ + "n_bins = 4\n", + "class CNN(nn.Module):\n", + " def apply(self, x):\n", + " b = x.shape[0]\n", + " x = nn.Conv(x, features=128, kernel_size=(4,), padding='SAME')\n", + " x = nn.BatchNorm(x)\n", + " x = nn.leaky_relu(x)\n", + " x = nn.avg_pool(x, window_shape=(2,), padding='SAME')\n", + " x = nn.Conv(x, features=256, kernel_size=(4,), padding='SAME')\n", + " x = nn.BatchNorm(x)\n", + " x = nn.leaky_relu(x)\n", + " x = nn.avg_pool(x, window_shape=(2,), padding='SAME')\n", + " x = x.reshape(b, -1)\n", + " x = nn.Dense(x, features=128)\n", + " x = nn.BatchNorm(x)\n", + " x = nn.leaky_relu(x)\n", + " x = nn.Dense(x, features=n_bins)\n", + " x = nn.softmax(x)\n", + " return x" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [], + "source": [ + "# Initialize model and optimizer\n", + "_, initial_params = CNN.init_by_shape(jax.random.PRNGKey(0), [((1, x_train.shape[1], 1), jnp.float32)])\n", + "model = nn.Model(CNN, initial_params)\n", + "\n", + "optimizer = optim.Adam(learning_rate=0.001).create(model)" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [], + "source": [ + "@jax.jit\n", + "def train_step(optimizer, batch):\n", + " def loss_fn(model):\n", + " w = model(batch['features'][..., jnp.newaxis])\n", + " return 1. / jax_metrics.compute_fom(w, batch['labels'])\n", + " \n", + " loss, g = jax.value_and_grad(loss_fn)(optimizer.target)\n", + " optimizer = optimizer.apply_gradient(g)\n", + " return optimizer, loss" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [], + "source": [ + "batch_size = 5000\n", + "def get_batch():\n", + " inds = onp.random.choice(len(z_train), batch_size)\n", + " return {'labels': z_train[inds], 'features': x_train[inds]}" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Epoch: 0, Loss: 0.010902259979624199\n", + "Epoch: 10, Loss: 0.0010201891737029918\n", + "Epoch: 20, Loss: 0.0007235596673232283\n", + "Epoch: 30, Loss: 0.0005803134087721451\n", + "Epoch: 40, Loss: 0.0005053445092730656\n" + ] + } + ], + "source": [ + "epochs = 50\n", + "\n", + "losses = []\n", + "\n", + "for epoch in range(epochs):\n", + " batch = get_batch()\n", + " optimizer, loss = train_step(optimizer, batch)\n", + " losses.append(loss)\n", + " if epoch % 10 == 0:\n", + " print(f'Epoch: {epoch}, Loss: {loss}')" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEGCAYAAAB7DNKzAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3deXzU1b3/8dcnk2RCFsKasIawhEUREIOAuCCi4oLaVutS63ql6nXrZu29/fWqt73trdVrXVtrrbVuVeuGWlFRlB0B2UG2EAh7EghkIev5/TETDCEzhDCTCTPv5+ORR2a+M/OdzymWN+ec7/ccc84hIiISSFykCxARkbZNQSEiIkEpKEREJCgFhYiIBKWgEBGRoOIjXUA4dOnSxWVnZ0e6DBGR48qiRYsKnXNdGx+PyqDIzs5m4cKFkS5DROS4Ymb5TR3X0JOIiASloBARkaAUFCIiEpSCQkREglJQiIhIUAoKEREJSkEhIiJBKSgaeHNxAS/P3xzpMkRE2hQFRQNTl27jlQUKChGRhhQUDaR44ymrrIl0GSIibUqbDwozu8zM/mxm75jZeeH8rlRvPGVVCgoRkYbCGhRm9pyZ7TKzFY2OTzKzr81svZndF+wczrm3nXO3ADcAV4axXJIT4ymrrA3nV4iIHHfCvSjg88ATwAv1B8zMAzwJnAsUAF+a2buAB/hNo8/f5Jzb5X/8C//nwibV66GsqgbnHGYWzq8SETluhDUonHNfmFl2o8OnAuudcxsBzOxV4FLn3G+Aixufw3x/Y/8W+JdzbnE4603xxuMcVFTXkpwYlQvriogctUjMUfQEtjR4XuA/FsidwETgcjO7NdCbzGyKmS00s4W7d+9uUWHJXl84lGpCW0TkoEj8s7mpMR0X6M3OuceAx450UufcM8AzALm5uQHPF0yq1wPgm6dIa8kZRESiTyR6FAVA7wbPewHbIlDHYVL8w026RFZE5BuRCIovgRwz62tmicBVwLsRqOMwKV4FhYhIY+G+PPYVYC4wyMwKzOxm51wNcAcwDVgNvOacWxnOOpqrPijKq3SJrIhIvXBf9XR1gOMfAB+E87tbIiXRN0ehyWwRkW+0+Tuzj4aZTTazZ0pKSlr0eQ09iYgcLqqCwjk31Tk3JT09vUWfPxgUGnoSETkoqoLiWNUPPalHISLyDQVFA/GeOLzxcQoKEZEGFBSNaAVZEZFDKSgaSfZ6tIKsiEgDCopGUhLjdXmsiEgDCopGUr3xlGvoSUTkoKgKimO9jwJ8K8iWauhJROSgqAqKY72PAvybF2noSUTkoKgKilBIToynXEEhInKQgqKRVK8ms0VEGlJQNJLi9VBWVYtzLdr7SEQk6igoGklOjKe2zlFZUxfpUkRE2gQFRSOpWkFWROQQCopGtHmRiMihoiooQnEfhTYvEhE5VFQFRSjuo9DmRSIih4qqoAgFbV4kInIoBUUjKV5tXiQi0pCCopGURF+PQnMUIiI+CopG6i+P1TIeIiI+CopGkuuHnjRHISICKCgO4433kOAxDT2JiPgpKJqgFWRFRL6hoGhCqjYvEhE5KKqCIhR3ZoN/BVn1KEREgCgLilDcmQ2+oacy7ZstIgJEWVCESqo3Xj0KERE/BUUTfENPmqMQEQEFRZNSNPQkInKQgqIJKRp6EhE5SEHRhBRvvO7MFhHxU1A0ISXRQ1VNHdW12jdbRERB0QRtXiQi8g0FRRNStDCgiMhBCoomqEchIvKNqAqK0C3hoc2LRETqRVVQhGoJj/pd7sp1052ISHQFRajUz1GoRyEioqBoUqrmKEREDlJQNCG5fuhJy3iIiCgompJ6cDJbcxQiIgqKJiQlxBFnGnoSEQEFRZPMTCvIioj4KSgC0AqyIiI+CooAtHmRiIiPgiIA31Lj6lGIiCgoAkhJ1NCTiAgoKALS0JOIiI+CIgANPYmI+ERVUIRq9VjQVU8iIvWiKihCtXos+LZD1dCTiEiUBUUopXjjqaiupbbORboUEZGIUlAEcHAFWc1TiEiMU1AEkKzNi0REAAVFQNq8SETER0ERgDYvEhHxUVAEUD/0pDkKEYl1CooAvulRaI5CRGKbgiKAZP8chYaeRCTWKSgC0OWxIiI+CooAUjSZLSICKCgCSk6ovzxWcxQiEtsUFAHExRnJiR7K1aMQkRinoAhCS42LiCgogtIKsiIiCoqgtCeFiIiCIqgUb7zWehKRmBdVQRHKHe7AN/RUXqWhJxGJbfGBXjCzkcE+6JxbHPpyjo1zbiowNTc395ZQnC/FG09+UXkoTiUictwKGBTAQmAlsNv/3Bq85oAJ4SqqrUjV0JOISNCg+DHwHaACeBV4yzlX2ipVtRHJifEaehKRmBdwjsI593/OudOBO4DewHQze83MRrRadRGW6vVQVlWDc9o3W0Ri1xEns51zecA7wEfAqcDAcBfVViR743EO9SpEJKYFm8zuB1wFXApswTf89Gvn3IFWqi3iUhqsIFv/WEQk1gT72289sAxfb2IfkAXcbuab03bOPRL26iIs9eCeFLWQFuFiREQiJFhQPIjv6iaA1Faopc05uB2qrnwSkRgWMCicc/e3Yh1tUqr2pBARCT6ZbWYXmNkXZlZoZrvN7HMzu7C1iou0FO1yJyISdDL7FuAHwL34br4DyAV+a2a9nHPPtEJ9EZWS2GCOQkQkRgWbo/ghcLpzrrjBsU/N7AJgFhD9QaGhJxGRoENP1igkAHDOFYWxnjalPii0jIeIxLJgQbHPzIY3Pug/tj98JbUd9UNPuuFORGLZkdZ6etfM/goswnep7CjgeuDaVqgt4uI9cXjj4zT0JCIxLdhaT7PwLdkRB9wA3OR/PMb/WkzQ5kUiEuuCXfWU5ZzbDPyyFetpc1K82rxIRGJbsDmKt+sfmNk/W6GWNiklUT0KEYltQa96avC4X7gLaatSvPGaoxCRmBYsKFyAxzElxRtPmYaeRCSGBbvqabiZ7cPXs2jnf4z/uXPOtQ97dW1AqtfDtr0VkS5DRCRigi0K6GnNQtqq5EQNPYlIbDviDnexLlVzFCIS4xQUR5Di9VBWVat9s0UkZikojiA5MZ7aOkdlTV2kSxERiYiAQWFm08zsh2Y2uDULOhZmNtnMnikpKQnZObV5kYjEumA9iuuBPcD9ZrbYzJ42s0vNrM1ui+qcm+qcm5Kenh6ycyZrTwoRiXHBrnraATwPPG9mccBo4ALgXjOrAD5yzv2uVaqMoFTtciciMS7YfRQHOefqgLn+n1+aWRfg/HAW1lZo8yIRiXXNCorGnHOFwEshrqVNSvH6hp603pOIxCpd9XQE9T2K1lhBVpfgikhbpKA4gpTE1tkO9VfvreK7f5qrsBCRNqdFQWFmN4a6kLaqteYo/rViB19u2sOi/D1h/R4RkaPV0h7FAyGtog2rn6MI59DTluJytvoXHnxxXn7YvkdEpCWC7XC3LNBLQGZ4yml7vPEeEjwW1qGn+XnFAIzt15kPlu/g/11cSedUb9i+T0TkaATrUWQC1wGTm/gpCn9pbUdLVpDNKyyjrq558w3zNxbRITmBBy49karaOl5bWNCSMkVEwiJYULwHpDrn8hv9bAJmtEp1bYRvBdnmDz1NXbqNs38/g+dm5zXr/fPzihmV3YmBmWmM7tuJlxfkNztkRETCLWBQOOduds7NCvDaNeErqe1J8Xqa3aNYsbWEn76xFPBNUB/J9pIKNheXM7pvJwC+P7YPW4or+Hzd7pYXLCISQro8thmSE+ObtYRHYWklU15YSKfkRK4dk8XizXsoKq0M+pn5G33zE2P6dQbgvBO60SXVy0ua1BaRNkJB0QzN2byoqqaO215cRHF5Fc9cl8tVo7JwDj5dsyvo5+bnFZGWFM+Q7r6dZRPj47hqVG+mr9lFwZ7ykLVBRKSlFBTNkJzoCTpH4Zzjv95dyZeb9vDQ5cMZ2jOdE3u0p1v7JKavPkJQbPTNT3ji7OCxq0dnYcArCzaHqgkiIi2moGiGVG980MtjX5yXzysLNnP7+P5MHt4DADPjnCEZfLFuNweqmw6ZXfsOsLGw7OD8RL2eHdoxYXAm//hyC1XaMElEIkxB0Qwp3njKA8xRzN1QxANTV3HO4Ax+ct6gQ16bOCST8qpa5m1s+mri+vsnRvvnJxq6dkwWhaVVfLjyyBPiIiLh1KLVY2NNstdDSUU133t2HnFmxJnhifP9XphfTHaXFB69agRxDYaPAMb270y7BA+frN7J+EEZh513fl4RKYkehvZof9hrZ+Z0JatTMi/Oy+cSfy9FRCQS1KNoholDMjm1bycOVNex/0ANe8ur2LX/AFv3VjAwM40/X5dLWlLCYZ9LSvBwRk4Xpq/e1eRif/M3FnNKdifiPYf/McTFGd8bncWCvGLW7twflnaJiDSHehTNMCq7E69OGduiz048IZOPVu1k5bZ9DO35zRatRaWVrNtVymUn9wz42Stye/Pwx2t5cV4+D146tEXfLyJyrNSjCLMJgzMw47Crnxbk1d8/0ampjwHQKSWRi07qzpuLtx7xfgwRkXBRUIRZl1QvJ/fuwCerdx5yfH5eMUkJcZzUs0PQz98+vj+VNbX817srw1mmiEhACopWMPGETJZvLWFHyYGDx+bnFXNKn44kxgf/I8jJTOOuCTm8t2w7HzZjSRARkVBTULSCiUN8q7JPX+PrVZSUV7Nmxz5G9z38stim3Dq+Pyf2aM8v3l7BnrKqsNUpItIUBUUryMlIJatT8sF5igWbinGOw260CyTBE8fvLh/G3vIq/vu9VeEsVUTkMAqKVlB/l/as9YWUV9Uwf2MRifFxDO8dfH6ioRN7pHP7+P68+dVWPl2z88gfEBEJEQVFKzl3SCZVNXXMWlfI/LxiRvTuQFKC56jOcceEHAZlpvHzN5dTUlEdpkpFRA6loGglo/p2Ii0pnreXbGXlthLGNHPYqaHE+DgeumIYu/dX8j/vrw5DlSIih1NQtJIETxzjB2XwwfId1Lmm13dqjmG9OjDlzP78Y+EWvlirzY1EJPwUFK1o4hDfek8JHmNkVscWn+eeiTn065rCz99czuYi7VkhIuGloGhF4wdm4IkzhvXqQLvEo5ufaCgpwcPvrxhOYWklEx6ewc/eWMaWYgWGiISH1npqRenJCdx7/iAGZKQe87lGZnXk85+ezdMz1vPKgi38c3EBl5/Si38/ewC9OyWHoFoRER9ralXT411ubq5buHBhpMtoNdtLKnh6xgZeXbCFOue4IrcX954/mI4piZEuTUSOI2a2yDmX2/i4hp6iQPf0djx46VA+v3c814zO4o1FBfzHW8sjXZaIRAkFRRSpD4w7J+TwrxU7Au6sJyJyNBQUUWjKmf3okZ7Eg1NXUVsXfUOLItK6FBRRKCnBw30XDmHV9n28sWhLpMsRkeOcgiJKTR7WnVP6dOShaWvZf0DLfYhIyykoopSZ8cuLT6CwtJKnZmw4qs/uKDnAszM3cumTsxnx4EeH7KMhIrGnzd9HYWZDgLuBLsB059zTES7puDG8dwe+PbInf5mZx9WjssjqHPj+iuKyKj5Yvp2pS7cdXAb9hO7tKamo5sV5+fzk/EGtWLmItCVh7VGY2XNmtsvMVjQ6PsnMvjaz9WZ2X7BzOOdWO+duBb4LHHZ9rwT3s0mD8cQZv/lX04sIbiku555Xv2LUrz/hF2+voLC0knvOGcj0H5/FB3efwcQhmby8YDMHqmtbuXIRaSvC3aN4HngCeKH+gJl5gCeBc4EC4EszexfwAL9p9PmbnHO7zOwS4D7/ueQoZLZP4vbx/Xn447XM21jEGP9ihMVlVTz+6TpenJePJ8648bRsvj2yF0O6p2FmBz9/47hsPl61k3eXbuO7ub0j1QwRiaCw35ltZtnAe865of7nY4H7nXPn+5//HMA51zgkmjrX+865iwK8NgWYApCVlXVKfn5+SOqPBgeqaznn4c9Jb5fA67eO5fk5m/jjjA2UVdXw3dze3DNxIN3Sk5r8rHOOSY/OxBNnvH/X6YeEiIhEl7Z0Z3ZPoOE1mwX+Y00ys/Fm9piZ/Qn4IND7nHPPOOdynXO5Xbt2DV21USApwcN9Fwxm1fZ9jP6f6Tw07WtG9+vMtHvO5LffGRYwJMA3KX7DuGxWbd/HgrziVqxaRNqKSExmN/VP0oDdGufcDGBGuIqJFRcP6847S7ayr6KGn04axKjs5m+cdNmInvzvh2t4fs6mFu+jISLHr0gERQHQcLC7F7AtAnXEFDPj2etHteiz7RI9XDUqi2e+2MDWvRX07NAuxNWJSFsWiaGnL4EcM+trZonAVcC7EahDjsL3x/YB4O9zNfcjEmvCfXnsK8BcYJCZFZjZzc65GuAOYBqwGnjNObcynHXIsevZoR3nn9iNVxZspqJKl8qKxJKwBoVz7mrnXHfnXIJzrpdz7i/+4x845wY65/o7534dzhokdG44LZuSimreXrK1ydcXb97DJU/M4vnZea1cmYiEU1Qt4WFmk83smZKSkkiXEpVO7duJE7q35/nZm2h4WXV1bR0Pf/Q1lz89h1Xb9vHf76/mq817IlipiIRSVAWFc26qc25Kenp6pEuJSvWXyn69cz9zN/j2uli/az/femo2j3+6nm+d3IsZPx1Pt/ZJ3P3qEkorayJcsYiEQlQFhYTfJcN70Cklkedm5/HcrDwuemwWW/dU8MdrR/Lwd4fTq2Myj141goI95fzynRVHPqGItHkKCjkqSQkerj61N5+s3sWD763itP6dmfbDM5k0tPvB94zK7sSdE3J4c/FW3gkwnyEix482v3qstD3Xj83my7w9XHZyT64+tXeTy3rcOWEAs9YX8ou3VjAyqyO9OwVeuVZE2rawr/UUCbm5uW7hwoWRLiPmbSku58I/zCQnM5XXfjCWeM83Hdi6Osfna3fz1zmbWL19HwO6pjK4expDurVnULc0Bmam0S7RQ1VNHVv2lLO5qJxNRWXkF5WzdW8F/bqmMKpPJ07p05GOKYkRbKVI9Ai01pOCQsLq3aXbuOuVr7jrnBx+dO5A9h+o5o1FBfxtziY2FZWT2d7LuP5d2FBYxtod+6nwL2duBl1SvRSVVtJw2+9UbzyZ7b1sLi6nutb3Qk5GKqP6dmJUdkfGD8xQcIi0UKCgiKqhJzObDEweMGBApEsRv0uG9+Dzr3fzxKfrKNhTzrQVOyirqmVkVgd+dN4gLhjajQR/T6OuzrG5uJw1O/axevt+tu6toEeHdvTplEx2l2T6dE6hc0oiZsaB6lqWbtnLwvw9fLmpmKlLtvHy/M0kJ3q4dkwf/u2MvmSkBV7sUESaTz0KCbvSyhoufmwmW/dWMHlYD64/LZvhvTuE9Dtq6xwrt5Xw3Kw83l26jXhPHFeN6s0PzuqvtalEmklDTxJRJRXV1NY5OrXCsNCmwjKenrGBN78qwDn49sie3DZ+AH27pIT9u0WOZwoKiTlb91bwzOcbeOXLLdTU1nHBSd257az+DO2pGzJFmqKgkJi1a/8B/jp7Ey/OzWd/ZQ1nDuzKbWf1Z0y/TtqxT6QBBYXEvH0HqnlxXj7PzcqjsLSKk7M6MOWMfpw9OIOkBE+Lz7thdyn/+dZyxvbrwt0Tc0JYsUjrUlCI+B2oruX1RQU888UGthRXkJLoYfzgDC4Y2o3xgzJI9TbvYkDnHP/4cgsPTF1FRXUtCR7js5+Mp1fH4/Pmwro6x3vLtzNxSAbJiVF1QaQ0k4JCpJGa2jpmbyhi2sodfLRyB4WlVSTGx3FmThcmDe3OhMEZASff95ZXcd8/l/Phyh2c1r8zPz1/EFf+aR7fHtmT335nWCu3JDSmLt3Gna98xe3j+3PvpMGRLkciICaCosF9FLesW7cu0uXIcaS2zrEofw8frtjBtJU72Lq3AjMYmdWRc4ZkcM7gTAZmpmJmzN1QxA//sYTC0kp+ev4gbjmjH3Fxxv3vruTv8/KZ/qOzyD7OrrByznHBH2ayZsd+UhI9zPzZhFa5Qk3alpgIinrqUcixcM6xYus+Plm9k0/X7GL5Vt/+Jr06tmNoj3SmrdpB384p/OGqkzmp1zdXUO3ad4AzH/qMC4d255ErR0Sq/Bb5ZNVO/u2Fhdx6Vn/+9MUGbjtLvYpYFBN3ZouEgplxUq90TuqVzg/PHciOkgN8umYXn67ZyZwNhVyZ25tfTj7hsHH8jPZJXDc2m2dnbuT2swcwICP1qL7XOcc7S7bx8oLN3HH2AM4c2DWUzQr6vU98tp5eHdvx4/MGsnVvBX+bs4lbzuin5VAE0DLjIkfULT2Ja0Zn8ez1o1h2//n89jvDAk72/uDMfiQleHj0k7VH9R1Ltuzl20/P4Z5/LGF5QQk3/HUBf5mVR2v0+OduKGLJlr384Kz+JHjiuGvCAMqra3l21sawf7ccHxQUIiHUOdXLTeP68t6y7azevu+I799RcoAf/WMJlz05m4I9FTx0+TAW/Oc5nHtCJv/93irufWMZlTW1Ya35ic/Wk5Hm5YpTegGQk5nGRSd15/nZm9hTVhXW75bjg4JCJMRuOaMfaUnx/N/HgXsVB6preXz6Os7+/QzeW7ad28f357OfjOeK3N6kJSXw9PdO4a5zcnh9UQFXPzOPXfsPhKXWxZv3MGdDEbec0e+Qe0nuOidHvQo5SEEhEmLpyQncckY/Plq1k2UFew95ra7O8c9FBZz9+xk8/PFazhrYlU9+dBb3Thp8yP0bcXHGj84dyJPXjGTV9n1c+sRsVvgn1UPpyU/X0yE5gWtGZx1yfKC/V/G3OfnqVYiCQiQcbhyXTYfkBB5p0KuYvb6Qix+fxY9fX0rXNC+vThnDH79/ClmdA9+gd9Gw7rxx62kYcPkf5/DszI0cqA7NUNSqbfuYvmYXN43rS0oTNxnedU4OZVU1/GVWXki+T45fCgqRMEhLSuDWs/oz4+vdvLZwCzf+dQHfe3Y+JRXV/OGqEbx9+zjG9OvcrHMN7ZnOO3eczqjsTvzq/dVM+P0MXvMvdHgsnpqxnlRvPNePzW7y9YGZaVx4Uneen9M6cxWhCkAJPQWFSJhcN7YPXVITufeNZSzM38N/XDiY6T8+i0tH9CQu7ugWI+ya5uXvN4/mxZtH0zXNy73/XMZ5j37B+8u2U1d39FdGbdxdyvvLt3PtmD6kJycEfN9dE1qnVzFvYxHDH/iIP36+IazfIy0TVfdRaIc7aUuSE+N56PLhLMrfw82n9w3JPQmn53Rh3IBxTFu5k4c/+pp/f3kxQ3u25/7JJ5Kb3anZ53l6xgYSPXHcfHrfoO8b1C2NC4f6ehWXjOjB7v2VbCwsI293GRsLS8krLKNP5xQevOTEFt+NXrCnnNtfWkx1bR2/n/Y1pw/ooqXg2xjdmS1ynKqtc7z91VYe+Xgtu0sreeqakUw8IfOIn9tUWMbERz7n2jF9uP+SE4/4/q937Of8R7845Fhyooe+XVLo0zmZmesKqa6t48fnDuKm0/viOYreUnlVDZc/PZcte8r5202nctuLi2iflMDUO08/phV9pWW0hIdIlNpbXsX1zy1g5bZ9/N+VI5g8vEfA9y7IK+b2lxZzoLqWj354Jj2auU3s1KXbKKmopl+XFPp1TSWzvffgXh479x3gP99awSerdzKidwd+d/kwBmamHfGczjnuePkrPlixneduGMXZgzKYuW433//LAm4a15dfTj6hef8DSMgECgrNUYgc5zokJ/Liv41mZFZH7n71K15buOWw9zjneGHuJq758zzSkuJ56/bTmh0SAJOH9+DaMX04bUAXuqUnHbLhU2b7JP583Sk8dvXJbC4u5+LHZvH49HVUH2Gy/akZG3h/+XbumzSYswdlAHBGTlduOC2b52bnMXt9YbPrk/BSj0IkSlRU1TLl7wuZua6QBy45ketPywZ8VxP9v7dX8PqiAs4ZnMEjV44gvV3gCexjUVRayf1TVzF16TYGZKTyvdFZXDai52HzM5+s2sktf1/IJcN78OiVIw4JnoqqWi5+fCblVbV8ePeZQSfbJbQ09CQSAyprarnj5a/4eNVO7p00iMtG9OTWFxexrKCEuyYM4J6JA4/6iquW+GjlDp74bD3LCkpI8BjnnpDJFaf05oycLuQVlvGtp+bQt0sKr986tsm5iGUFe/n2U3O4aFh3/nDVyWGvV3wUFCIxorq2jp+8vpR3lmwjOdFDnBkPf3c455/YrdVrWb19H68vLODtJVspLqsiI82LJ86ornVMvXMc3dMDD389Nn0dj3y8lsevPjnovEu4OOeYvb6IjPbeZs25RAMFhUgMqa1zPDh1JYs27+HRK0cwICOyf9FV1dTx6ZqdvL6wgMWb9/Dn63KPeDlvTW0dl/9xLnmFZbx35+n07tR6W8x+uamY//3XGhbm76FdgocnrjmZc4Yc+Yqy452CQkSOO3mFZVz4h5lUVNfSu1M7hvXqwPBe6Qzr1YGhPdPxxseRX1TGup2lrN1Zytpd+1m3cz9llbVMHJLBxcN7cEpWx2YPt63ZsY+HPvya6Wt2kZHm5bbx/Xnrq62s2FrCA5ecyPcD3MUeLRQUInJcWrtzP9NX72JZwV6WFZSwdW8FAGYQ7x/Gqn/eu2MyORmpxMUZX6zdTWVNHd3Tk7jopO5MHt6DYb3SD5k4r61zlFfVsHPfAZ76bANvLdlKqjee28b358bT+tIu0UN5VQ13vfIVn6zexZQz+3HfpMGtMs8TCQoKEYkKhaWVLC8oYWnBXg5U15GTkcrAzDT6Z6QcsqFUaWUNn6zayXvLtvH52t1U1zoy23tJ8MRRXlVLWWUNlTXfXMLrjY/jhnHZ3HZWfzokH3qVVm2d44GpK3lhbj4XndSdh787PCpvCFRQiEjMKimvZtqqHcxeX4gnzkhJjCc50UNyYjwpXg8p3njOHpRBt/SkgOdwzvHszDx+/cFqTunTkT9fl0unNrRV7OLNe/jjjA08dvXJLQ6xmAiKBms93bJu3bpIlyMiUej9Zdv54WtLSIqPY2z/zpzWvwtj+3cmJyP1kGGt1rKnrIrfTVvDKwu20D09iedvPJVB3Vp28UJMBEU99ShEJJyWbtnLS/Pzmb2+6OCcSZdUL2P7d2ZE7w4kJcThMcMT980PwP4DNZRUVLOvopq95efpA8wAAAYuSURBVNWUVPh+qmvrqHOOOufrudQ533BX707tOP/EbpwzOPOwGw/r6hxvLC7gt/9aQ0lFNTeNy+aeiQOb3FukuQIFRVStHisi0hqG9+7A8N4dANhSXM7cDUXM2VDInA1FTF267Yif98bHkd4ugQ7JCbRPSsCbEEecGWZGnIHH3zNZuqWEaSt3Eh9njO3fmfNO7Mb5J2RSXF7FL95awcL8PeT26civvjWUwd3ah6296lGIiISIc46isipqah21zlFX56it8z12DtKS4klvl9DsOYS6OseyrSV8uGIH01buIK+wDIA4g/R2Cfz8wiFcPrJXyK7C0tCTiMhxzDnHul2lfLhiBweqa7nljH4h2eOkIQ09iYgcx8yMgZlpEVlORMuMi4hIUAoKEREJSkEhIiJBKShERCQoBYWIiASloBARkaAUFCIiEpSCQkREgoqqG+7qV48F9plZS5eP7QIUhq6q44baHVvU7tjTnLb3aepgVC7hcSzMbGFTt7BHO7U7tqjdsedY2q6hJxERCUpBISIiQSkoDvdMpAuIELU7tqjdsafFbdcchYiIBKUehYiIBKWgEBGRoBQUfmY2ycy+NrP1ZnZfpOsJJzN7zsx2mdmKBsc6mdnHZrbO/7tjJGsMBzPrbWafmdlqM1tpZnf7j0d1280sycwWmNlSf7sf8B+P6nbXMzOPmX1lZu/5n0d9u81sk5ktN7MlZrbQf6zF7VZQ4PsPCXgSuAA4AbjazE6IbFVh9TwwqdGx+4DpzrkcYLr/ebSpAX7snBsCjAH+3f/nHO1trwQmOOeGAyOASWY2huhvd727gdUNnsdKu892zo1ocO9Ei9utoPA5FVjvnNvonKsCXgUujXBNYeOc+wIobnT4UuBv/sd/Ay5r1aJagXNuu3Nusf/xfnx/efQkytvufEr9TxP8P44obzeAmfUCLgKebXA46tsdQIvbraDw6QlsafC8wH8slmQ657aD7y9UICPC9YSVmWUDJwPziYG2+4dflgC7gI+dczHRbuBR4F6grsGxWGi3Az4ys0VmNsV/rMXtjqq1no6BNXFM1w1HKTNLBf4J3OOc22fW1B9/dHHO1QIjzKwD8JaZDY10TeFmZhcDu5xzi8xsfKTraWXjnHPbzCwD+NjM1hzLydSj8CkAejd43gvYFqFaImWnmXUH8P/eFeF6wsLMEvCFxEvOuTf9h2Oi7QDOub3ADHxzVNHe7nHAJWa2Cd9w8gQze5HobzfOuW3+37uAt/ANr7e43QoKny+BHDPra2aJwFXAuxGuqbW9C1zvf3w98E4EawkL83Ud/gKsds490uClqG67mXX19yQws3bARGANUd5u59zPnXO9nHPZ+P4//alz7lqivN1mlmJmafWPgfOAFRxDu3Vntp+ZXYhvPNMDPOec+3WESwobM3sFGI9v2eGdwH8BbwOvAVnAZuAK51zjCe/jmpmdDswElvPNmPV/4JuniNq2m9kwfJOXHnz/OHzNOfegmXUmitvdkH/o6SfOuYujvd1m1g9fLwJ80wsvO+d+fSztVlCIiEhQGnoSEZGgFBQiIhKUgkJERIJSUIiISFAKChERCUpBIdICZlbrX5mz/idkC8uZWXbDlX1FIk1LeIi0TIVzbkSkixBpDepRiISQfx+A//Xv/7DAzAb4j/cxs+lmtsz/O8t/PNPM3vLvFbHUzE7zn8pjZn/27x/xkf+OapGIUFCItEy7RkNPVzZ4bZ9z7lTgCXx3++N//IJzbhjwEvCY//hjwOf+vSJGAiv9x3OAJ51zJwJ7ge+EuT0iAenObJEWMLNS51xqE8c34dskaKN/AcIdzrnOZlYIdHfOVfuPb3fOdTGz3UAv51xlg3Nk41sKPMf//GdAgnPuV+Fvmcjh1KMQCT0X4HGg9zSlssHjWjSfKBGkoBAJvSsb/J7rfzwH3wqmAN8DZvkfTwdug4ObC7VvrSJFmkv/ShFpmXb+HePqfeicq79E1mtm8/H9Q+xq/7G7gOfM7KfAbuBG//G7gWfM7GZ8PYfbgO1hr17kKGiOQiSE/HMUuc65wkjXIhIqGnoSEZGg1KMQEZGg1KMQEZGgFBQiIhKUgkJERIJSUIiISFAKChERCer/Ay8EPn7TPO6RAAAAAElFTkSuQmCC\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "plt.plot(range(epochs), losses)\n", + "plt.xlabel('Epoch')\n", + "plt.ylabel('1 / FOM')\n", + "plt.yscale('log')" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": {}, + "outputs": [], + "source": [ + "batch_size = 10000\n", + "Ngalaxies = len(x_test)\n", + "tomo_bin = jnp.concatenate(\n", + " [optimizer.target(x_test[batch_size * i : min((batch_size * (i + 1)), Ngalaxies)][..., jnp.newaxis]) \n", + " for i in range(Ngalaxies // batch_size + 1)]\n", + " )\n", + "\n", + "tomo_bin = jnp.argmax(tomo_bin, axis=-1)" + ] + }, + { + "cell_type": "code", + "execution_count": 16, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmkAAAFzCAYAAABl1J6yAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3df3Rc9X3n/9c7RoCwYxuMcSzLXpHYBBOLTWsdaiKWr4+9G9HAhrRHpO45FG2/UFOXbFi+3+8pYtsufNPDrtjlxAlNocuvRmZTfhynifnieLUsCu2BY5w4ENsB24sS9DXy0PDD2AHVJBZ97x9zr301mpHmx52Ze2eej3N0NPOZe68+uh5ZL31+mrsLAAAAyfKRelcAAAAAUxHSAAAAEoiQBgAAkECENAAAgAQipAEAACQQIQ0AACCBTqt3BeJ27rnnekdHR72rAQAAMKMf/ehHb7v7wnyvNVxI6+jo0O7du+tdDQAAgBmZ2f9f6DW6OwEAABKIkAYAAJBAM4Y0M3vYzN40s5/kee3/MTM3s3MjZbeZ2YiZHTSznkj5ajPbF7x2j5lZUH6GmT0elO8ys47IOX1m9mrw0VfpNwsAAJAWxYxJ+6akb0jaEi00s6WS/pWkQ5GyiyRtkPQpSW2S/qeZXeDuH0q6T9JGSS9I+p6kKyTtkHS9pHfdfbmZbZB0l6TfMbNzJN0uqUuSS/qRmT3p7u+W/+0CAIAkOHHihMbGxvTBBx/Uuyo1ceaZZ6q9vV0tLS1FnzNjSHP3v4+2bkVslvTHkrZFyq6W9Ji7/1LSa2Y2IukSMxuVNNfdd0qSmW2R9AVlQ9rVku4Izt8q6RtBK1uPpKfd/UhwztPKBrtHi/7uAABAIo2NjemjH/2oOjo6FHSuNSx31zvvvKOxsTGdf/75RZ9X1pg0M/u8pMPuvifnpSWSXo88HwvKlgSPc8snnePuE5KOSVowzbUAAEDKffDBB1qwYEHDBzRJMjMtWLCg5FbDkpfgMLOzJP2JpM/mezlPmU9TXu45uXXaqGxXqpYtW5bvEAAAkDDNENBC5Xyv5bSkfULS+ZL2BN2Y7ZJeNLOPKdvatTRybLukTFDenqdc0XPM7DRJ8yQdmeZaU7j7/e7e5e5dCxfmXQ8OAADgpNHRUa1atSrvazfccINeeeWVoq/l7vryl7+s5cuX6+KLL9aLL74YSx1Lbklz932SzgufB0Gty93fNrMnJf2NmX1V2YkDKyT9wN0/NLP3zGyNpF2SrpP0F8ElnpTUJ2mnpF5Jw+7uZjYk6T+a2dnBcZ+VdFs53yQAAEi27oFhHT56PLbrLZnfquf715V17oMPPljS8Tt27NCrr76qV199Vbt27dKmTZu0a9eusr521IwhzcwelbRW0rlmNibpdnd/KN+x7v6ymT0h6RVJE5JuCmZ2StImZWeKtio7YWBHUP6QpEeCSQZHlJ0dKnc/YmZ/LumHwXFfCScRAACAxnL46HGNDlwZ2/U6+rfPeMzExIT6+vr00ksv6YILLtCWLVt01llnae3atbr77rvV1dWlOXPm6Oabb9ZTTz2l1tZWbdu2TYsWLZp0nW3btum6666TmWnNmjU6evSo3njjDS1evLii72HG7k53/113X+zuLe7enhvQ3L3D3d+OPL/T3T/h7p909x2R8t3uvip47Uvu7kH5B+5+jbsvd/dL3P1nkXMeDsqXu/tfV/SdAgAARBw8eFAbN27U3r17NXfuXN17771TjhkfH9eaNWu0Z88eXX755XrggQemHHP48GEtXXpqhFZ7e7sOHz5ccf3YcQAAADSlpUuXqru7W5J07bXX6rnnnptyzOmnn66rrrpKkrR69WqNjo5OOSZod5okjkkRhDQAANCUcoNUvmDV0tJysnzWrFmamJiYckx7e7tef/3UqmFjY2Nqa2uruH6ENKCKugeG1dG/Xd0Dw/WuCgAgx6FDh7Rz505J0qOPPqrLLrusrOt8/vOf15YtW+TueuGFFzRv3ryKx6NJhDSgLMWGr3AgbJwzlgAA8Vi5cqUGBwd18cUX68iRI9q0aVNZ1/nc5z6nj3/841q+fLn+4A/+IO/YtnKUvAQHgFPhK3f2UDiFvJKp3wDQjJbMby1qRmYp15tOR0dHwbXQnn322ZOP33///ZOPe3t71dvbO+V4M9Nf/uVfllfRaRDSgBgVCm8AgOnxh+1UhDQgBtEWNAAA4kBIA2KQuwhj2GxPaAMAlIuQBlSgUBibqdmesWsAgJkQ0oAKlBuwGLsGAJgJS3AAAAAkECENqKFwfTXGqgFAfY2OjmrVqlV5X7vhhhsKLs+Rz4EDB3TppZfqjDPO0N133x1XFenuBGopd4IBACCwuVM6dii+681bJt2yr6xTH3zwwZKOP+ecc3TPPffou9/9bllfrxBCGgAAqL9jh6Q7jsV3vTvmzXjIxMSE+vr69NJLL+mCCy7Qli1bdNZZZ2nt2rW6++671dXVpTlz5ujmm2/WU089pdbWVm3btk2LFi2adJ3zzjtP5513nrZvj3ecMd2dAACgKR08eFAbN27U3r17NXfu3LzbOY2Pj2vNmjXas2ePLr/8cj3wwAM1qx8hDQAANKWlS5equ7tbknTttdfqueeem3LM6aefrquuukqStHr1ao2OjtasfoQ0oAZY3BYAksfMpn0uSS0tLSfLZ82apYmJiZrUTWJMGlAThdZTi4Y3FrUFgNo6dOiQdu7cqUsvvVSPPvqoLrvssnpXaRJCGlBHYTBjUVsAqL2VK1dqcHBQN954o1asWKFNmzaVdZ1/+Id/UFdXl37xi1/oIx/5iL72ta/plVde0dy5cyuqHyENSABa1AA0vXnLipqRWdL1ptHR0VFwLbRnn3325OP333//5OPe3l719vZOOf5jH/uYxsbGyqvnNAhpQAmie27GiRY1AE2vzDXNGhkhDSgBi9ECAGqF2Z0AAAAJREgDAABIIEIaAABAAhHSAAAAEoiQBiRIuBRH98BwvasCAA1tdHRUq1atyvvaDTfcUHB5jny+9a1v6eKLL9bFF1+sz3zmM9qzZ08sdWR2J5AgLMUBoFn1bO1RZjwT2/XaZrdpqHeorHMffPDBko4///zz9Xd/93c6++yztWPHDm3cuFG7du0q62tHEdIAAEDdZcYz2tcX31ppnYOdMx4zMTGhvr4+vfTSS7rgggu0ZcsWnXXWWVq7dq3uvvtudXV1ac6cObr55pv11FNPqbW1Vdu2bdOiRYsmXeczn/nMycdr1qyJbWFbujuBInQPDLNBOgA0mIMHD2rjxo3au3ev5s6dq3vvvXfKMePj41qzZo327Nmjyy+/XA888MC013zooYf0m7/5m7HUj5AGFCFcxJYtmwCgcSxdulTd3d2SpGuvvVbPPffclGNOP/10XXXVVZKk1atXa3R0tOD1vv/97+uhhx7SXXfdFUv96O4EAABNycymfS5JLS0tJ8tnzZqliYmJvNfau3evbrjhBu3YsUMLFiyIpX60pAEAgKZ06NAh7dy5U5L06KOP6rLLLiv7Or/927+tRx55RBdccEFs9aMlDZhGtTZUBwDU38qVKzU4OKgbb7xRK1as0KZNm8q6zle+8hW98847+qM/+iNJ0mmnnabdu3dXXD9z94ovkiRdXV0ex40BpOxSGPXYUD0aDhkHB6AR7d+/XytXrjz5PElLcFRL7vcsSWb2I3fvync8LWlAArFeGoBmk7RAlQSMSQMAAEigGUOamT1sZm+a2U8iZf/FzA6Y2V4z+46ZzY+8dpuZjZjZQTPriZSvNrN9wWv3WDBVwszOMLPHg/JdZtYROafPzF4NPvri+qYBAACSrpiWtG9KuiKn7GlJq9z9Ykn/S9JtkmRmF0naIOlTwTn3mtms4Jz7JG2UtCL4CK95vaR33X25pM2S7gqudY6k2yX9hqRLJN1uZmeX/i0CAACkz4whzd3/XtKRnLL/4e7hQiEvSGoPHl8t6TF3/6W7vyZpRNIlZrZY0lx33+nZmQpbJH0hcs5g8HirpPVBK1uPpKfd/Yi7v6tsMMwNiwAAAA0pjjFp/6ekHcHjJZJej7w2FpQtCR7nlk86Jwh+xyQtmOZaU5jZRjPbbWa733rrrYq+GQAAgCSoKKSZ2Z9ImpD0rbAoz2E+TXm550wudL/f3bvcvWvhwoXTVxoAADS90dFRrVq1Ku9rN9xwg1555ZWir7Vt2zZdfPHF+vSnP62urq6820uVo+wlOIKB/FdJWu+nFlsbk7Q0cli7pExQ3p6nPHrOmJmdJmmest2rY5LW5pzzbLn1BQAAyTWybr1OZOJbJ62lrU3Lh58p69wHH3ywpOPXr1+vz3/+8zIz7d27V1/84hd14MCBsr52VFkhzcyukHSrpP/D3f8x8tKTkv7GzL4qqU3ZCQI/cPcPzew9M1sjaZek6yT9ReScPkk7JfVKGnZ3N7MhSf8xMlngswomKAAAgMZyIpPRygP7Y7ve/gtXznjMxMSE+vr69NJLL+mCCy7Qli1bdNZZZ2nt2rW6++671dXVpTlz5ujmm2/WU089pdbWVm3btk2LFi2adJ05c+acfDw+Pp53D9ByFLMEx6PKBqhPmtmYmV0v6RuSPirpaTP7sZn9lSS5+8uSnpD0iqT/Lukmd/8wuNQmSQ8qO5ngpzo1ju0hSQvMbETS/yWpP7jWEUl/LumHwcdXgjIAAICKHTx4UBs3btTevXs1d+5c3XvvvVOOGR8f15o1a7Rnzx5dfvnleuCBB/Je6zvf+Y4uvPBCXXnllXr44YdjqV8xszt/190Xu3uLu7e7+0Puvtzdl7r7p4OPP4wcf6e7f8LdP+nuOyLlu919VfDal8IuUnf/wN2vCa55ibv/LHLOw0H5cnf/61i+YwAAAElLly5Vd3e3JOnaa6/NO5bs9NNP11VXXSVJWr16tUZHR/Ne67d+67d04MABffe739Wf/dmfxVI/toUC8mBjdQBofLndkvm6KVtaWk6Wz5o1SxMTE1OOibr88sv105/+VG+//bbOPffciurHtlBAHoePHtfowJVsbg4ADezQoUPauXOnJOnRRx/VZZddVtZ1RkZGFM6hfPHFF/WrX/1KCxYsqLh+tKQBAICmtHLlSg0ODurGG2/UihUrtGnTprKu8+1vf1tbtmxRS0uLWltb9fjjj8cyeYCQBgAA6q6lra2oGZmlXG86HR0dBddCe/bZZ08+fv/9908+7u3tVW9v75Tjb731Vt16663lVXQahDQAAFB35a5p1sgYkwYk2JL5rero367ugeF6VwUAUGO0pAEJFk5c6OjfXueaAABqjZY0AABQF6d2lWx85XyvhDQAAFBzZ555pt55552mCGrurnfeeUdnnnlmSefR3QkAAGquvb1dY2Njeuutt+pdlZo488wz1d7eXtI5hDQAAFBzLS0tOv/88+tdjUSjuxMAACCBCGkAAAAJREgDAABIIEIaAABAAhHSAAAAEoiQBgAAkECENAAAgAQipAEAACQQIQ0AACCBCGlACiyZ36qO/u3qHhiud1UAADXCtlBACjzfv06S1NG/vc41AQDUCi1pAAAACURIAwAASCBCGgAAQAIR0gAAABKIkAYAAJBAhDQAAIAEIqQBAAAkECENAAAggQhpAAAACcSOA0BE98CwDh89riXzW+tdFQBAkyOkARGHjx7X6MCV9a4GAAB0dwIAACQRLWmA6OYEACQPIQ0Q3ZwAgOSZsbvTzB42szfN7CeRsnPM7GkzezX4fHbktdvMbMTMDppZT6R8tZntC167x8wsKD/DzB4PyneZWUfknL7ga7xqZn1xfdMAAABJV8yYtG9KuiKnrF/SM+6+QtIzwXOZ2UWSNkj6VHDOvWY2KzjnPkkbJa0IPsJrXi/pXXdfLmmzpLuCa50j6XZJvyHpEkm3R8MgAABAI5sxpLn730s6klN8taTB4PGgpC9Eyh9z91+6+2uSRiRdYmaLJc11953u7pK25JwTXmurpPVBK1uPpKfd/Yi7vyvpaU0NiwAAAA2p3Nmdi9z9DUkKPp8XlC+R9HrkuLGgbEnwOLd80jnuPiHpmKQF01wLAACg4cW9BIflKfNpyss9Z/IXNdtoZrvNbPdbb71VVEUBAACSrNyQ9vOgC1PB5zeD8jFJSyPHtUvKBOXteconnWNmp0map2z3aqFrTeHu97t7l7t3LVy4sMxvCQAAIDnKDWlPSgpnW/ZJ2hYp3xDM2Dxf2QkCPwi6RN8zszXBeLPrcs4Jr9UraTgYtzYk6bNmdnYwYeCzQRkAAEDDm3GdNDN7VNJaSeea2ZiyMy4HJD1hZtdLOiTpGkly95fN7AlJr0iakHSTu38YXGqTsjNFWyXtCD4k6SFJj5jZiLItaBuCax0xsz+X9MPguK+4e+4EBgAAgIY0Y0hz998t8NL6AsffKenOPOW7Ja3KU/6BgpCX57WHJT08Ux0BAAAaDTsOoKmxHRQAIKkIaWhqbAcFAEiquJfgAAAAQAwIaQAAAAlESAMAAEggQhoAAEACEdIAAAASiJAGAACQQIQ0AACABGKdNDQlFrEFACQdIQ1NiUVsAQBJR3cnAABAAhHSAAAAEoiQBgAAkECENAAAgAQipAEpsmR+qzr6t6t7YLjeVQEAVBmzO4EUeb5/nSSpo397nWsCAKg2QhqmGFm3XicyGbW0tWn58DP1rg4AAE2J7k5McSKT0coD+3Uik6l3VQAAaFqENAAAgAQipAEAACQQY9JQUEtbm/ZfuHLSc8aoAQBQG4Q0FJQbyKKBDQAAVBchDSdnc4Za2trqWBsAACAR0qBTszkBAEByMHEAAAAggQhpAAAACURIAwAASCBCGgAAQAIR0gAAABKI2Z1NqJk3UO8eGNbho8e1ZH5rvasCAMC0CGlNKFxyY2Tdeu2/cGVTrYt2+OhxjQ5cWe9qAAAwI0Jag5uu1azUVrRwm6hmbIEDAKDWGJPW4MJWs+iOAuVaPvxMbNcCAADToyWtwTTzeDMAABoJLWkNJs6WMyTXkvmt6ujfru6B4XpXBQBQJbSkpVSpLWbheLLwMdLt+f51kqSO/u11rgkAoFoqCmlmdoukGyS5pH2Sfl/SWZIel9QhaVTSF9393eD42yRdL+lDSV9296GgfLWkb0pqlfQ9STe7u5vZGZK2SFot6R1Jv+Puo5XUuVGELWZh8IqGNmnyIH+p9EkCAACgvsoOaWa2RNKXJV3k7sfN7AlJGyRdJOkZdx8ws35J/ZJuNbOLgtc/JalN0v80swvc/UNJ90naKOkFZUPaFZJ2KBvo3nX35Wa2QdJdkn6n3Do3sjC0hQhlAACkW6Vj0k6T1GpmpynbgpaRdLWkweD1QUlfCB5fLekxd/+lu78maUTSJWa2WNJcd9/p7q5sy1n0nPBaWyWtNzOrsM4AAACJV3ZIc/fDku6WdEjSG5KOufv/kLTI3d8IjnlD0nnBKUskvR65xFhQtiR4nFs+6Rx3n5B0TNKC3LqY2UYz221mu996661yv6VUyu3WrKVwMdyRdetr/rUBAGh0ZYc0Mztb2Zau85XtvpxtZtdOd0qeMp+mfLpzJhe43+/uXe7etXDhwukr3mDCtcvq0b0Z7WIlrAEAEK9KJg78S0mvuftbkmRmfyvpM5J+bmaL3f2NoCvzzeD4MUlLI+e3K9s9OhY8zi2PnjMWdKnOk3SkgjqnThLXPSs0KSGcxAAAACpXSUg7JGmNmZ0l6bik9ZJ2SxqX1CdpIPi8LTj+SUl/Y2ZfVbblbYWkH7j7h2b2npmtkbRL0nWS/iJyTp+knZJ6JQ0H49aaRu4sziQoFBaTvG0UG6sDANKm7JDm7rvMbKukFyVNSHpJ0v2S5kh6wsyuVzbIXRMc/3IwA/SV4PibgpmdkrRJp5bg2BF8SNJDkh4xsxFlW9A2lFtfVF+SW9TYWB0AkDYVrZPm7rdLuj2n+JfKtqrlO/5OSXfmKd8taVWe8g8UhDwAAIBmwrZQAAAACURIAxpUz9Ye9WztqXc1AABlIqQBDSoznp0k3TnYSVgDgBRig3WggQ31DknKBjUAQLrQkobYhUtxsLgtAADloyUNsUvyUhwAAKQFLWkpQesUStGztUdts2u/nysAID6EtJQI9+iUVLcN1ZEemfHMyfFoAIB0orszZZK23RLSoW12m3q29hDcACBFaEkDmsBQ79DJJTkAAOlASAMAAEggQhoAAEACEdKAFFsyv1Ud/dvVPTBc76oAAGLGxAFU3ci69TqRyailrY2JDzF7vn+dJKmjf3udawIAiBstaai6E5mMVh7YrxMZBq4DAFAsQhqqJlyAlzXdAAAoHd2dqBq6NgEAKB8taQAAAAlESxoaWvfAsA4fPa4l81vrXRUAAEpCSENDO3z0uEYHrqx3NWqmZ2uPMuMZNlcHgAZASAMaSGY8o319++pdDQBADBiTBgAAkEC0pCVMuPBriOUrAABoToS0hAkXfm1E0XXTWJ4DAIDpEdJQM2Ew23/hyjrXBACA5GNMGmoubFEbWbe+3lUBACCxaElDzdGiBgDAzGhJAwAASCBCGuqGbk8AAAqjuxN1Q7cnAACF0ZIGAACQQLSkAUiUcEFn1tMD0OxoSQOQKOGCztGdNwCgGdGSBiARoi1oAABCGoCEaOQt0QCgHHR3AgAAJFBFIc3M5pvZVjM7YGb7zexSMzvHzJ42s1eDz2dHjr/NzEbM7KCZ9UTKV5vZvuC1e8zMgvIzzOzxoHyXmXVUUl8AyTOybr32X7iSbk4AyFFpS9rXJf13d79Q0j+XtF9Sv6Rn3H2FpGeC5zKziyRtkPQpSVdIutfMZgXXuU/SRkkrgo8rgvLrJb3r7sslbZZ0V4X1BRpWz9Yetc0uHHTaZrepc7BTPVt7Th7fOdg5qawewm5OZnICwGRlj0kzs7mSLpf0byTJ3X8l6VdmdrWktcFhg5KelXSrpKslPebuv5T0mpmNSLrEzEYlzXX3ncF1t0j6gqQdwTl3BNfaKukbZmbu7uXWG82he2BYh48e15L5rfWuSs1kxjPa17ev4OtDvUOSpM7BzinHh2UAgOSoZOLAxyW9JemvzeyfS/qRpJslLXL3NyTJ3d8ws/OC45dIeiFy/lhQdiJ4nFsenvN6cK0JMzsmaYGkt6MVMbONyrbEadmyZRV8S2gUh48e1+jAlfWuRmLN1OoGAKi/Sro7T5P065Luc/dfkzSuoGuzAMtT5tOUT3fO5AL3+929y927Fi5cOH2tEyYcj8P+lailzHjmZMsaACCZKglpY5LG3H1X8HyrsqHt52a2WJKCz29Gjl8aOb9dUiYob89TPukcMztN0jxJRyqoc+KwcCeSpGdrT13HpwEATik7pLn7P0h63cw+GRStl/SKpCcl9QVlfZK2BY+flLQhmLF5vrITBH4QdI2+Z2Zrglmd1+WcE16rV9Iw49GA6smMZ5QZ5w8GAEiCShez/beSvmVmp0v6maTfVzb4PWFm10s6JOkaSXL3l83sCWWD3ISkm9z9w+A6myR9U1KrshMGdgTlD0l6JJhkcETZ2aEAGgA7DADA9CoKae7+Y0ldeV7KO8DK3e+UdGee8t2SVuUp/0BByGt0LW1trBWFmkjKhAF2GACA6bEtVEKwRhRqZaYJA22z29SztafuEwuif7jw8wGgGbEtFIBJhnqHyhqXFvdM5eXDzzCpBkBToyWtxqLjcGgdyApbTKLPuTe1Fe5GUElXaNh9Gf23BACUj5BWY/wimyo3kHFvaq+Srk0mAABAdRDSAFSECQAAUB2MSUPihN2f7MIAAGhmtKQhccLuT7o9GwvjMQGgNLSkAagJtkADgNIQ0gDUFAs3A0Bx6O5EQ+keGNbho8e1ZH5rvavS8Mqd1UlXJwAUh5CGhnL46HGNDlxZ72rUVM/WHn10Zabm2z0xqxMAqovuzjqhywdxyYxn9N7+gapv4xT3jgIAgOnRklYndPkgbViIGQBqi5Y0AACABCKkAQAAJBAhDQAAIIEYkwY0EVb9B4D0oCUNaCJxrPpf65nJ7OUKoFnRkobEioYBWn2So9b/FuzlCqBZEdKQWPxyToZydxYAAFSGkAY0gUqCFjsLAEB9ENKAJlDLoJXbTU1LHACUh5AGNIAl81vV0b9dS+a36vn+dXWtS243NS1xAFAeQlqN0JqAagqDWUf/9jrXBAAQF0JajdCaUD5medbeX93n2v+fsi1haf/DgrXhAKQVIQ2JV8wsz+6BYR0+elxL5rfWqlqJ0LO1R22z4w9R5xz9UCsP7FfnYKf29aU72LAxPIC0IqShIRw+elyjA1fWuxo1lxnPaF/fvnpXAwBQBew4AKCgttlt6tnaU+9qlGVk3fqa7owAAHEjpCE12B6o9oZ6h5QZL38LqXoKuzkZhwYgrejurDJmdcZgc6ckaflwtluPsUUAgGZASKsyZnXG4Nih+K+5uTN73XnLpFsad0xXHH8ktM1uU+dgp9pmt2modyjG2gEApkNIQ7Jt7swGqXxBLWhhmzFk5Tvu2CHpjmPSHfPiqWdCxfFHQhjMOgc746hS2aq1FAtLdABIKkIakm26MBUGt82deu6Mf5SUZ3bn5kiwCK8xb1ns1UT1FbMUSzFyQxlLdABIKkIa0iloYXvj2HF9ePQf1W5vT34t7MoMQ15Ug7eeSTN3czbzAsFhKMud/dnM9wRAMhHSkA7zlmXDVdgKFoSvS/u3Z9dHiwavMJiFXaX5rtXgZurmjKtVKs1ygxj3BEDSENKQDuF4sjvmSSpyEHyhsWoNPFEAANA4Kl4nzcxmmdlLZvZU8PwcM3vazF4NPp8dOfY2Mxsxs4Nm1hMpX21m+4LX7jEzC8rPMLPHg/JdZtZRaX2RcJs7s0FscxmD1MNzm6ClLA2i3YcAgNLFsZjtzZKi/Sr9kp5x9xWSngmey8wukrRB0qckXSHpXjObFZxzn6SNklYEH1cE5ddLetfdl0vaLOmuGOqLJAu7KqdZdqOlrU37H2ubvKhtGMzuOFZ6S9nmzvJCYZ1Va9/OuCwffobFZAGgAhWFNDNrV3ZK3YOR4qslDQaPByV9IVL+mLv/0t1fkzQi6RIzWyxprrvvdHeXtCXnnPBaWyWtD1vZ0LyW/+s3NX/DCZ3IZE5tqH7LvvK6MaNj3FImMw6wFWsAABeKSURBVJ5h3TIAaGCVtqR9TdIfS/qnSNkid39DkoLP5wXlSyS9HjluLChbEjzOLZ90jrtPSDomaUGFdUYazFtWuHXr2CFd+sHXJUnP96+r7OuE4W66r5ci3xy6s+m7GNmzE0CjKHvigJldJelNd/+Rma0t5pQ8ZT5N+XTn5NZlo7LdpVq2jPFIDeGWfflnZ4bPP6jC12uApTkWHX+36Xe4YJcPAI2ikpa0bkmfN7NRSY9JWmdm/03Sz4MuTAWf3wyOH5O0NHJ+u6RMUN6ep3zSOWZ2mqR5ko7kVsTd73f3LnfvWrhwYQXfEhIlXxdmud2axWiQ1rRqaZvdpp6tPTMfCACIRdkhzd1vc/d2d+9QdkLAsLtfK+lJSX3BYX2StgWPn5S0IZixeb6yEwR+EHSJvmdma4LxZtflnBNeqzf4GlNa0tAgCq1rViu37Evl2LRaGeodUmY8M/OBNRZ2b06aSJIHs00BpE011kkbkPSEmV0v6ZCkayTJ3V82syckvSJpQtJN7v5hcM4mSd+U1CppR/AhSQ9JesTMRpRtQdtQhfqiXnL31My3OwCqKo2r7OfWudhtneL6/sLdHMK6pOW+AUifOJbgkLs/6+5XBY/fcff17r4i+Hwkctyd7v4Jd/+ku++IlO9291XBa18KW8vc/QN3v8bdl7v7Je7+szjqi4SI7L2ZGE3W5RkukyGpqNaoJAjrHAalUK1aysJQmLb7BiB9YglpQMnCrs0KuhjDX8qx/oJs0i7PQsEnTeqxLlsaQy6A9CCkoT6OHap4AkDVgkWTtaalVa3HmE23tAdhDUA1sHcnkqHekwaiostxbO7MBsqw1Q+JUeuxYMUs7cEm7QDiREsakiGGlrWqKGKbKgAAqoGQhvqbtyw5rWhRSWrdS4hmXCuNHQwA1Avdnai/JLagheEsrFs4Ti2Jda2hod4hdQ4213g9djAAUC+0pKH20tBClbuzQfiYCQUoQlVmHgNoOoQ01F7M489q9guxSZfnQOkaYUkTAPVHd2eVhKuSM44lR4WtaN0Dwzp89LiWzG89WVbTGXV0e6ptdps6BzvVNrtNQ71D9a4OADQsQlqVMI6lgAq3fjp89LhGB66MsUIlii7P0aTCYNboY9Pi+EMrjdtuAUgOQhpqg7FcSJk4/tBi3TQAlSCkIX65G6dLpxaELbOrM183JwAAjYyQhvhFB9dHW9AqGMdV927OKMalAQBqgJCGeES3T5Kyn++Yl/3caDMib9mX/X6jY9PYNgoAEDNCGuIRTggIZ29GA0sjjkfLDWTR0JbiwFbrTcsBAIUR0mLW9Etv5AsnKQ0sJYl+j1We/dmztUeZ8YzaZsf/HitlBmK4RVSjLcNRjaDKLE8A5SCkxYylNzCpq7cKATUzntG+vvoH39wtosI9PdMe2qoRopjlCaAchDQkWipndUa3kMqdYJBv5muKRRe2DVv3WOgWAOLBtlCoTBhEqiSc1fl8/7qqfY2qybeN1LFDDTWRYqh36GSrXhjMwudhyxpOYU9PAKUgpKEyYeBIwIbpifwFGC7X0eCGeocmtZwN9Q4pM86+lbnY0xNAKejuROUS0nWXyHE/TbyNFF2fAFAZQhrKV+Fm6U2jSRe/bZY9PuMQzgoPMQsUgER3J4pRaNzZsUNNFzzKEo5NI9Qix8i69Sdbflce2H/yQ1Lervvw+ER16QOoGlrSMLMGGuheV+GCv1Xo/gxbYn7eerYS1NmLGRRasie36z66/uLKA/uT1aUPoGoIaUADCH/Z/2b/do3WuzKITXQRXNZfBJoPIQ0AaqzYXQ0YlwY0N0IaCssdhxZdiJXxVaUp9l412GK3yI/wBaAYhDRMFQaF6Fi0aMio4pZHDSt6r6ab7cn4PwBAgJCGqY4dygaJaDCrcSBL5XZQxQpbIvNtGTVvGUENACCJkIZC6txKFm4H1bDyLXIbzv7MF+AAAE2HkIaGG2sWHZSd6LE/hbaMarBdCtpmt6lnaw+7DsQofI9Hnyf6vQ6gLIQ0TO5eC1tzUiyR20PlU0YY69nao7bZ088ITJqh3iF2HYhZbiALF7klrAGNhZDW7KItZ+GEgDpq6LFo+ZR4vzPjGe3roxsUk6XmDxMAJSGkNbuEtZw1/Fi0XIw7AwAUwN6dAAAACURIAwAASCBCGpBE85Zlxwjmm/0JFBDO+hxZt77eVQEQg7LHpJnZUklbJH1M0j9Jut/dv25m50h6XFKHpFFJX3T3d4NzbpN0vaQPJX3Z3YeC8tWSvimpVdL3JN3s7m5mZwRfY7WkdyT9jruPlltnBKK/+FO+3EbDCseqNdBSHKi+ak8gGFm3XicyGWaRAjVSSUvahKT/291XSloj6SYzu0hSv6Rn3H2FpGeC5wpe2yDpU5KukHSvmc0KrnWfpI2SVgQfVwTl10t6192XS9os6a4K6ltV0SnwiXfs0KmPhAxc7x4YVkf/9lhnddKqAMTrRCajlQf2SxI/W0ANlB3S3P0Nd38xePyepP2Slki6WtJgcNigpC8Ej6+W9Ji7/9LdX5M0IukSM1ssaa6773R3V7blLHpOeK2tktabmZVb52oK//NK/F+X4ZIbuds+1Vk4q/P5/nWxXXP58DNaeWC/TmQysV2z5goteJsy4YK2qI1q/4HSED9bQArEsgSHmXVI+jVJuyQtcvc3pGyQM7PzgsOWSHohctpYUHYieJxbHp7zenCtCTM7JmmBpLdzvv5GZVvitGxZcoJHIiVsyQ3MoEF2H2BB29qKq9uT7k2gviqeOGBmcyR9W9K/c/dfTHdonjKfpny6cyYXuN/v7l3u3rVw4cKZqtycNncmYrFaxCtN3exts9vUOdhJi1oNlduiFr6vJE3q3kzD+wxoJBW1pJlZi7IB7Vvu/rdB8c/NbHHQirZY0ptB+ZikpZHT2yVlgvL2POXRc8bM7DRJ8yQdqaTOTSnsLqMFLfW+ce+E9v+nU9v/RMcIJV24dyctarVTbota7vuKVjSgPspuSQvGhj0kab+7fzXy0pOS+oLHfZK2Rco3mNkZZna+shMEfhB0jb5nZmuCa16Xc054rV5Jw8G4NRQrunk60mneMvU8uFKdg50675gYC4SSMYkGSKdKWtK6Jf2epH1m9uOg7N9LGpD0hJldL+mQpGskyd1fNrMnJL2i7MzQm9z9w+C8TTq1BMeO4EPKhsBHzGxE2Ra0DRXUtzkxBi39btmnzGCn9r12SPvV3N1N4d6uoSXzW2OdbNKoqrU0Rxj+GLMGVEfZIc3dn1P+MWOSlPfPNXe/U9Kdecp3S1qVp/wDBSEPaHrzlin7903zyt3btaN/ex1r03iiEwWKwcbuQHWxwXqj2dyZbT0LJwgkfKJA2DIS5/poDeuWfdJ/zf4yjLZgAMUq1PIVDWdpGeMINANCWqMJuzcTOhYtX3dVtGUERdjcqeXDyfp3LUW4Zlo4kaASS+a3nlwEmW7PmRVq+UrTBBSgmRDSGlXCwlkot7sKZTh2qN41qEica6aFwYxuz9KELWrR53Fdj/FpQHwIaRUqdQwHANRb3CEqej3GpwHxIaRVKBHdBOE4NCnxY9CAamL2J4BGQkhLq9xgxjIbALM/EyC3KzVaTjcoUBpCWlqx/lnRWMsp/cqdBRxOLIg+p2Wtugr9jNENCpSOkIaGx1pO6VfuhJPcQEbLWv3wxxJQOkJa0iV0KY1SsR5a+Xq29qhtdnNNTIm+X4pp+YouxYFk4o8loHSEtKRL+XIL0V+2LL1Rnsx4Rvv60h3SSxW2nBXb8kUXZnrM1KIWzpiPHk/LG5oVIS0N5i3LtqiFrWmbOxM/i5NwVh8s7oqky21Riy5jtHz4mSkz5kfWraebFE2LkJYGt+yT7ph36nkKJg2waG19pH1x19xucboxG1e0RW3lgf2TwlgU3aRoZoS0NMndlxNNY/KiyZkZj0+r3HBPa2Djym0Vo5UMmIqQlhZhMEt4CxqqY1IX0ObOyd3fDYCJJQAwFSEtLVLyC5lftjWQ2/2dIuH7Y+6Kc9Q52CmbOEe/ePWPEzN2sdRZpaidfPuN0vqGRkdIQ6wYi4Z8cieSdPRLowNXqnOwM1Hvl1JnlaJ2cgMZY9TQDAhpiEUaWtBYTLN28u2hGQ1j4YSAjyb09yyzZJMv9+eZpTvQiAhpiEUaWtDSOEvs1EK2Oevl5S7LkjAzvR/C4NM52F+rKk0r94+MtM+SbQbhz3N0Vmh06Y40/ZwDhRDSkiwF66GhusKFbPf/ac4vnISOS0tDi2o+afgjA/nRWoZGRkhLshSshwZENWrYydd9SzcogGojpKEiaW05AaSZ37/RsWnR8Ek3aPIxBhWNgJCWRClatLZRW05QmrSE9dwJAcWOnUP65BuzRlhD2hDSkohuzoaVu09hIT1be5QZzwSTBgpIwOSBSvdobZvdpp6tPRrqHZKU/b4lnXwetzB0dQ8Ms91Uk0jjhCEgREgr0+RtemLEZIGGFu4cMNMvjHDCwLTqOHmg0nAWGuodUs/WnpNBLQymnYOdapvdVvWwVikWv00Puj+RRh+pdwXSKvxlG/sP+7FDiV1WoVGE/1mPrFtf76qkTtgCJWUXo40jmIRBLBrM9vXtU2Y8uXuUht2mkk6G1I7+7eoeGK5ntTCN5cPPnFyio5Kf/7D7lP8/UAu0pNVDdMxZSgNZWsYg5ZOE7o/p/qo/tTZaFVtsy1StMYjVajGrltxwyrpq6VHuWLXoz2K0NbzYIQxAOQhp9RCOOcvtqkpRVycTBiozXVDMjGf0ncHztP9Ppy7QOUUCxqVVS+54NSBOuWEtVChshb0n0eOii+gy5g3VQEirlelmbG7uzH5mwgACub8QCqrRuLR6tJwO9Q6pc7CzZl8vDmwnlT65gSw3tIVyW7Nzz4tzzFuxW1zRitf4CGm1Uqj1LHxt3rLUtKIhfuFsTkn5t4GaTg1a02g5LQ7dnulXbtgppmVupvCV26Uayg2NucexzEjjIqTVWjSIhb9cpYbsrsLMwr++/8P8Wep+4Scny6dsAzWdhG4RVU3RJUqS2B0anVgQPqdlrTlM1zI3U/gq1IKe20qXexxrwjUuQlqtRcPYLftSNQ5NSveEgVzVnJJfTDdEz9YeZfreVNvsZfranx6acl5JclrT4up2S+q/d7hESc/Wnqov11GO3HtOy1rzmu7/lvD/oOjz6a4RDWHFHpevpS769QhzyUZIq6Zidg5IWQtaI3V7VXOW50zroYWLtoZroY3cu37SIOSS5bSmVdrtFtc6aHHLXeQ3unxHktGyhnxKDUjFHj9d1+t0LXlIHkJaHPJ1WYZldxxLXWsZ4pOvtS5fy09S/ppNWjiLzvDMDbZpQssa6mGm/1fyteQl5f8iZBHS4nAsZ5B3bmhLWWtZPknt9kqa3O7K3L9oj8yflciQkbRwFgrDWRK7NCtByxqSoND4uekmNJTSfVqL2aeNPsOVkBaXectOdTeleJHaQhqpmzNX9K/JSn/QCw38vemPTpO0LLEhI8n/vkm9Z8XKtx9pbiCL7uQgEdpQHzN1k87UfVrsedMpd3HhYte7SxtCWlwaLJSFmqEFLfqDXO4YjZkG/Be1F2cVzDSBoFH+fZO08G10ORUpW7eZtrgK/23Cc9+bOEcd/cdPvk5oQy0VCjczhZ5yz4sqtjWv0AzXQtdJa4tbKkKamV0h6euSZkl60N0H6lylUxp8vFmSW1iSoNC6RqHcge61VmgCQVK7N8uVhIVvo//WuYE8WrdCy4dEx9z1bO2Rr+w/ecxMLW25wTDUSF3EhYTv5VxpDba5309av49yFduaN9MM+HwzXdO4plziQ5qZzZL0l5L+laQxST80syfd/ZX61kzZ7s0G7NpsZqUOpJ1pZ4B6taDlyjcGqhHCWUliWJNwuvXZpvu3bpvddjKohSEuHGsXPSa8Zvj55DGLpY8uPnW9X0jqHDz13CbO0b7rp37t3K+RW6c0B7iZ/tBI6+SM3D+Mq/V9RO9fEkNgof93+z77J6fqXcZ1CoXAXEkJcYkPaZIukTTi7j+TJDN7TNLVkhIQ0hpnC6fp/hptJjNtEZOv6byQ6EbpNZUnjCTxP+GaiN6L3Ak+BRRqlZJOBazc4DPTv3W+MFRMQCo2ROW2tJ3yZUn5W2OmC3D51CrUFduSlJRW/kJhJ64QlPsHVr7Xy7l+eP9KHQ8Zd0tfsS2huaE8rHe5X3+mAJaUFjdz97p98WKYWa+kK9z9huD570n6DXf/Ur7ju7q6fPfu3VWv1/4LV5a3nlWdNVrXwEym+4Urlf6Lp+hFauuxGn7uunz5QknaW37D7zGPnvY2ZVqyf3e2nZjQ0Fhm8r2IDksocA9ODvJ/PTPjcbnjzpLcKlXo574Usz8xoI+cfnRK+T/9ar7Gf9pf0bWjCv1ynum4XB392yeFuJnCVKX1LRReZrp+7nGl/l9cbv1nGqda7Hlx3b+Z6lHMuNrp7ne54bMWv+vN7Efu3pX3tRSEtGsk9eSEtEvc/d9GjtkoaWPw9JOSDlahKudKersK121W3M94cT/jxf2MH/c0XtzPeNXzfv4zd1+Y74U0dHeOSVoaed4uaVLTiLvfL+n+albCzHYXSrooHfczXtzPeHE/48c9jRf3M15JvZ8fqXcFivBDSSvM7HwzO13SBklP1rlOAAAAVZX4ljR3nzCzL0kaUnYJjofd/eU6VwsAAKCqEh/SJMndvyfpe3WuRlW7U5sQ9zNe3M94cT/jxz2NF/czXom8n4mfOAAAANCM0jAmDQAAoOkQ0iLM7AozO2hmI2Y2ZeEfy7oneH2vmf16PeqZJkXc07VmdszMfhx8/Id61DMNzOxhM3vTzH5S4HXenyUo4n7y3iyBmS01s++b2X4ze9nMbs5zDO/REhR5T3mfFsnMzjSzH5jZnuB+/r95jknWe9Td+ch2+c6S9FNJH5d0uqQ9ki7KOeZzknZIMklrJO2qd72T/FHkPV0r6al61zUNH5Iul/Trkn5S4HXen/HeT96bpd3PxZJ+PXj8UUn/i/9Da3JPeZ8Wfz9N0pzgcYukXZLW5ByTqPcoLWmnnNx+yt1/JSncfirqaklbPOsFSfPNbHHuhXBSMfcURXL3v5d0ZJpDeH+WoIj7iRK4+xvu/mLw+D1J+yUtyTmM92gJirynKFLwvns/eNoSfOQOzE/Ue5SQdsoSSa9Hno9p6g9DMcfglGLv16VB8/MOM/tUbarWkHh/xo/3ZhnMrEPSrynbUhHFe7RM09xTifdp0cxslpn9WNKbkp5290S/R1OxBEeNWJ6y3IRdzDE4pZj79aKyW2K8b2afk/RdSSuqXrPGxPszXrw3y2BmcyR9W9K/c/df5L6c5xTeozOY4Z7yPi2Bu38o6dNmNl/Sd8xslbtHx6Um6j1KS9opM24/VeQxOKWYLb1+ETY/e3Y9vBYzO7d2VWwovD9jxHuzdGbWomyY+Ja7/22eQ3iPlmime8r7tDzuflTSs5KuyHkpUe9RQtopxWw/9aSk64LZH2skHXP3N2pd0RSZ8Z6a2cfMzILHlyj7nnyn5jVtDLw/Y8R7szTBvXpI0n53/2qBw3iPlqCYe8r7tHhmtjBoQZOZtUr6l5IO5ByWqPco3Z0BL7D9lJn9YfD6Xym768HnJI1I+kdJv1+v+qZBkfe0V9ImM5uQdFzSBg+m2GAyM3tU2Zlc55rZmKTblR34yvuzDEXcT96bpemW9HuS9gVjfiTp30taJvEeLVMx95T3afEWSxo0s1nKhtkn3P2pJP+eZ8cBAACABKK7EwAAIIEIaQAAAAlESAMAAEggQhoAAEACEdIAAAASiJAGAACQQIQ0AACABCKkAcA0zOwPzezHwcdrZvb9etcJQHNgMVsAKEKwh+KwpP/s7v9fvesDoPHRkgYAxfm6pGECGoBaYe9OAJiBmf0bSf9M0pfqXBUATYTuTgCYhpmtljQo6V+4+7v1rg+A5kF3JwBM70uSzpH0/WDywIP1rhCA5kBLGgAAQALRkgYAAJBAhDQAAIAEIqQBAAAkECENAAAggQhpAAAACURIAwAASCBCGgAAQAIR0gAAABLofwNVwo7zALRgKwAAAABJRU5ErkJggg==\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "fig = plt.figure(figsize=(10, 6))\n", + "for i in range(n_bins):\n", + " plt.hist(z_test[tomo_bin == i], bins=250, histtype='step', label=f'bin {i}')\n", + "plt.xlabel('z')\n", + "plt.legend();" + ] + }, + { + "cell_type": "code", + "execution_count": 17, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "{'SNR_3x2': 1098.281927142067,\n", + " 'FOM_3x2': 2644.3098007411236,\n", + " 'FOM_DETF_3x2': 45.625914165827254}" + ] + }, + "execution_count": 17, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "jax_metrics.compute_scores(tomo_bin, z_test, metrics=['SNR_3x2', 'FOM_3x2', 'FOM_DETF_3x2'])" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python [conda env:.conda-tomo-jax]", + "language": "python", + "name": "conda-env-.conda-tomo-jax-py" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.8.5" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} diff --git a/notebooks/Scores.ipynb b/notebooks/Scores.ipynb new file mode 100644 index 00000000..5eefe892 --- /dev/null +++ b/notebooks/Scores.ipynb @@ -0,0 +1,119 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "import json\n", + "import matplotlib.pyplot as plt" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "{'2_bins': {'FOM_3x2': 1387.4266199398603, 'FOM_DETF_3x2': 23.735751442234772, 'SNR_3x2': 876.809259538073}, '3_bins': {'FOM_3x2': 2502.150561803356, 'FOM_DETF_3x2': 38.19791609649463, 'SNR_3x2': 973.6463977803064}, '4_bins': {'FOM_3x2': 3312.096090570743, 'FOM_DETF_3x2': 57.57785850852204, 'SNR_3x2': 1144.1731253194073}, '5_bins': {'FOM_3x2': 3932.3143288232827, 'FOM_DETF_3x2': 63.04511436371448, 'SNR_3x2': 1308.1209736720286}, '6_bins': {'FOM_3x2': 5377.125291203648, 'FOM_DETF_3x2': 76.95578516521786, 'SNR_3x2': 1263.3259662420758}, '7_bins': {'FOM_3x2': 7416.639452785952, 'FOM_DETF_3x2': 95.94608323827052, 'SNR_3x2': 1395.499023177545}}\n" + ] + } + ], + "source": [ + "json_file = \"../result.json\"\n", + "with open(json_file, 'r') as jf:\n", + " results = json.load(jf)\n", + "print(results)" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmsAAAGDCAYAAAB0s1eWAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3df5RfdX3n8efLRGcwKahFs4TQglNqqwgiKaWr26ZSCVprOBW2WBC264q11Gprtyv80W17lq3dblt1T9WibglNWk4aC7JdlVBk2toiCNbyUwqJv8JEEd0Ck3ZGwPf+8b2Rb8IkzIR85/vJzPNxzvd87/3cH5/38Dnoi/u593tTVUiSJKlNTxt2AZIkSdo7w5okSVLDDGuSJEkNM6xJkiQ1zLAmSZLUMMOaJElSwwxrkiRJDTOsSZIkNcywJmlBSPLFJP+aZLLvszLJSJLfTvLlbvs9Sf5zkvQdO56kkpywxzmv6trXPEnfZye5O8mDSe5Psj7JobOs+392NT2c5PNJztuvfwCSFizDmqSF5KeqannfZwL4c+BU4NXAdwFvAC4A3rPHsf8EfCcoJflu4BTg67Po9++Al1XVYcDzgaXAf5tlzTuBnwIOA84H3pPk387yWEmLgGFN0oKV5FTgNOB1VXV7VT1aVZ8GzgUuTPJ9fbtvBH4myZJu/fXAlcC3nqyfqvpKVT3Q1/QY8H1dDWNJvpnkpd36yiQP7LpaV1X/tao+X1Xfrqobgb8FfuSp/N2SFhbDmqSF7JXAjVX1lf7GLhRtp3fFbZcJ4E564Q56V9kun21HSV6e5EHgYeB1wLu7vrYC/wXYmOSZwB8Dl1XV+AznOAT4IeCO2fYraeFbOuwCJOkAuirJo93yOPAAsGMv++4ADt+j7XLgvCTbgGdV1Q19t7btU1V9CjgsyZHAm4Av9m37YJKfAm4ECnjtXk7zAeAfgWtm1amkRcEra5IWkjOq6lnd5wx6Ye2Ivex7RLe9318ArwDeCvzJ/hRQVfcBnwCu2GPTB4HjgP9VVdN7Hpfkd7vt/76qan/6lrQwGdYkLWR/BfxwkqP6G5OcDBwFfLK/var+Bfg48Bb2M6x1lgJjff0tpzct+mHgN5I8Z496fhN4FXBaVT30FPqVtAAZ1iQtWFX1V8B1wEeSvCjJkiSn0HuY4P1Vdc8Mh10M/FhVfXG2/SQ5J8n3pOd7gUu6fnd5D3BLVf0n4P/Sm+7cdexFwM8Cr6yqb8zxT5S0CBjWJC10rwOupzc1OQlsoHeF660z7VxVE939Z3PxQuDvu/P/HXA3vfvWSLIOOB34+W7fXwFemuScbv2/A98D3NP3+3AXz7F/SQtYvDVCkiSpXV5ZkyRJaphhTZJmIcnH93iVlVOWkuaF06CSJEkN88qaJElSwxb0GwwOP/zwOvroowfax86dO1m2bNlA+9DcOCZtclza45i0yXFpz3yNyS233PJAVT13z/YFHdaOPvpobr755oH2MT4+zpo1awbah+bGMWmT49Iex6RNjkt75mtMknxppnanQSVJkhpmWJMkSWqYYU2SJKlhhjVJkqSGGdYkSZIaZliTJElqmGFNkiSpYYY1SZKkhhnWJEmS9mZigsNuvRUmJoZWgmFNkiRpT5OTsG4djI3x4osvhrGx3vrk5LyXYliTJEna0znnwJYtMDXF0p07YWqqt37uufNeimFNkiSp3333fSeo7WZqCq65Zt6nRA1rkiRJ/bZtg5GRmbeNjMDWrfNajmFNkiSp39gYTE/PvG16urd9HhnWJEmS+q1cCaedBqOju7ePjsLatb3t88iwJkmStKeNG3vBbHSUR5Ytezyobdgw76UsnfceJUmSWrd8OVx1FUxMcPvmzZx45pnzfkVtF6+sSZIk7c3KlTx4/PFDC2pgWJMkSWqaYU2SJKlhhjVJkqSGGdYkSZIaZliTJElqmGFNkiSpYYY1SZKkhhnWJEmSGmZYkyRJaphhTZIkqWGGNUmSpIYZ1iRJkho20LCW5AVJPtf3eSjJ25M8J8m1Se7pvp/dd8xFSe5NcneStX3tJyW5rdv23iQZZO2SJEktGGhYq6q7q+olVfUS4CTgX4ArgXcC11XVscB13TpJXgicDbwIOB14X5Il3eneD1wAHNt9Th9k7ZIkSS2Yz2nQU4GtVfUlYB2wvmtfD5zRLa8Drqiq6ar6AnAvcHKSI4BDq+qGqirg8r5jJEmSFqyl89jX2cCfdcsrqmoHQFXtSPK8rv1I4NN9x2zv2h7plvdsf4IkF9C7AseKFSsYHx8/UPXPaHJycuB9aG4ckzY5Lu1xTNrkuLRn2GMyL2EtyTOA1wIXPdmuM7TVPtqf2Fh1KXApwOrVq2vNmjWzL3Q/jI+PM+g+NDeOSZscl/Y4Jm1yXNoz7DGZr2nQVwGfraqvdetf66Y26b7v79q3A0f1HbcKmOjaV83QLkmStKDNV1h7PY9PgQJcDZzfLZ8PfLSv/ewkI0mOofcgwU3dlOnDSU7pngI9r+8YSZKkBWvg06BJngm8EnhzX/O7gE1J3gh8GTgLoKruSLIJuBN4FLiwqh7rjnkLcBlwCPDx7iNJkrSgDTysVdW/AN+9R9s36D0dOtP+lwCXzNB+M3DcIGqUJElqlW8wkCRJaphhTZIkqWGGNUmSpIYZ1iRJkhpmWJMkSWqYYU2SJKlhhjVJkqSGGdYkSZIaZliTJElqmGFNkiSpYYY1SZKkhhnWJEmSGmZYkyRJaphhTZIkqWGGNUmSpIYZ1iRJkhpmWJMkSWqYYU2SJKlhhjVJkqSGGdYkSZIaZliTJElqmGFNkiSpYYY1SZKkhhnWJEmSGmZYkyRJaphhTZIkqWGGNUmSpIYZ1iRJkhpmWJMkSWqYYU2SJKlhhjVJkqSGGdYkSZIaZliTJElqmGFNkiSpYYY1SZKkhhnWJEmSGmZYkyRJaphhTZIkqWGGNUmSpIYZ1iRJkho28LCW5FlJNif5fJK7kvxIkuckuTbJPd33s/v2vyjJvUnuTrK2r/2kJLd1296bJIOuXZIkadjm48rae4BPVNUPACcAdwHvBK6rqmOB67p1krwQOBt4EXA68L4kS7rzvB+4ADi2+5w+D7VLkiQN1UDDWpJDgR8FPgxQVd+qqn8G1gHru93WA2d0y+uAK6pquqq+ANwLnJzkCODQqrqhqgq4vO8YSZKkBWvpgM//fODrwB8nOQG4BXgbsKKqdgBU1Y4kz+v2PxL4dN/x27u2R7rlPdufIMkF9K7AsWLFCsbHxw/YHzOTycnJgfehuXFM2uS4tMcxaZPj0p5hj8mgw9pS4KXAW6vqxiTvoZvy3IuZ7kOrfbQ/sbHqUuBSgNWrV9eaNWvmVPBcjY+PM+g+NDeOSZscl/Y4Jm1yXNoz7DEZ9D1r24HtVXVjt76ZXnj7Wje1Sfd9f9/+R/UdvwqY6NpXzdAuSZK0oA00rFXVV4GvJHlB13QqcCdwNXB+13Y+8NFu+Wrg7CQjSY6h9yDBTd2U6cNJTumeAj2v7xhJkqQFa9DToABvBTYmeQawDfg5eiFxU5I3Al8GzgKoqjuSbKIX6B4FLqyqx7rzvAW4DDgE+Hj3kSRJWtAGHtaq6nPA6hk2nbqX/S8BLpmh/WbguANbnSRJUtt8g4EkSVLDDGuSJEkNM6xJkiQ1zLAmSZLUMMOaJElSwwxrkiRJDTOsSZIkNcywJkmS1DDDmiRJUsMMa5IkSQ0zrEmSJDXMsCZJktQww5okSVLDDGuSJEkNM6xJkiQ1zLAmSZLUMMOaJElSwwxrkiRJDTOsSZIkNcywJkmS1DDDmiRJUsMMa5IktWJigsNuvRUmJoZdiRpiWJMkadgmJ2HdOhgb48UXXwxjY731yclhV6YGGNYkSRq2c86BLVtgaoqlO3fC1FRv/dxzh12ZGmBYkyRpmO677ztBbTdTU3DNNU6JyrAmSdJQbdsGIyMzbxsZga1b57ceNcewJknSMI2NwfT0zNump3vbtagZ1iRJGqaVK+G002B0dPf20VFYu7a3XYuaYU2SpGHbuLEXzEZHeWTZsseD2oYNw65MDVg67AIkSVr0li+Hq66CiQlu37yZE8880ytq+g6vrEmS1IqVK3nw+OMNatqNYU2SJKlhhjVJkqSGGdYkSZIaZliTJElqmGFNkiSpYYY1SZKkhhnWJEmSGmZYkyRJaphhTZIkqWEDD2tJvpjktiSfS3Jz1/acJNcmuaf7fnbf/hcluTfJ3UnW9rWf1J3n3iTvTZJB1y5JkjRs83Vl7cer6iVVtbpbfydwXVUdC1zXrZPkhcDZwIuA04H3JVnSHfN+4ALg2O5z+jzVLkmSNDTDmgZdB6zvltcDZ/S1X1FV01X1BeBe4OQkRwCHVtUNVVXA5X3HSJIkLVhL56GPArYkKeCPqupSYEVV7QCoqh1JntfteyTw6b5jt3dtj3TLe7Y/QZIL6F2BY8WKFYyPjx/AP+WJJicnB96H5sYxaZPj0h7HpE2OS3uGPSbzEdZeVlUTXSC7Nsnn97HvTPeh1T7an9jYC4OXAqxevbrWrFkzx3LnZnx8nEH3oblxTNrkuLTHMWmT49KeYY/JwKdBq2qi+74fuBI4GfhaN7VJ931/t/t24Ki+w1cBE137qhnaJUmSFrSBhrUky5J8165l4DTgduBq4Pxut/OBj3bLVwNnJxlJcgy9Bwlu6qZMH05ySvcU6Hl9x0iSJC1Yg54GXQFc2f3KxlLgT6vqE0k+A2xK8kbgy8BZAFV1R5JNwJ3Ao8CFVfVYd663AJcBhwAf7z6SJEkL2kDDWlVtA06Yof0bwKl7OeYS4JIZ2m8GjjvQNUqSJLXMNxhIkiQ1zLAmSZLUMMOaJElSwwxrkiRJDXvSsJbk0CS/neRPkvzsHtveN7jSJEmSNJsra39M7w0CH6H3G2gfSTLSbTtlYJVJkiRpVmFtrKreWVVXVdVrgc8Cn0zy3QOuTZIkadGbze+sjSR5WlV9G3q/g5ZkO/A3wPKBVidJkrTIzebK2v8BXtHfUFXrgXcA3xpEUZIkSep50itrVfVre2n/BL13d0qSJGlAZv3THd3ToIf1rX9vkusGU5YkSZJgbr+z9ingxiSvTvIm4Frg3YMpS5IkSTCHF7lX1R8luQO4HngAOLGqvjqwyiRJkjSnadA3AP8bOA+4DPhYkhMGVJckSZKYw5U14HXAy6vqfuDPklxJL7SdOIjCJEmSNLdp0DP2WL8pyQ8f+JIkSZK0y5OGtSQBzgIK2EzvN9fWAZ8HPjDQ6iRJkha52VxZ+0PgecAz6IW0EXo/lPtq4AXA2wZWnSRJ0iI3m7D276rqxUmeDnwVOKKqvpXkT4F/GGx5kiRJi9tsngZ9FKCqHgE+U1Xf6tYfBR4bYG2SJEmL3mzC2leTLAeoqtN3NSb5N/huUEmSpIF60rBWVa+qqskZNj0MvGbXSpIXHcjCJEmSNLfXTe2mqnZ2v7m2y58cgHokSZLUZ7/D2gxyAM8lSZIkDmxYqwN4LkmSJHFgw5okSZIOsAMZ1nwyVJIk6QCbzeumXrqv7VX12e77lANVlCRJknpm8waDm4E7gK936/0PEhS9d4VKkiRpAGYT1t4BvA74V+AK4Mq9/O6aJEmSDrDZ/CjuH1TVy4FfBI4CrkuyKclLBl6dJEnSIjfrBwyq6gvAR4EtwMnA9w+qKEmSJPXM5gGD5wNnA+uAr9CbCr2kqqYGXJskSdKiN5t71u4FbqV3Ve0h4HuAX0h6zxlU1e8PrDpJkqRFbjZh7bd4/O0EywdYiyRJkvbwpGGtqn5jHuqQJEnSDGb1gEGSVyX5myQPJPl6kr9O8upBFydJkrTYzeYBgzcBbwZ+jd4P5AKsBt6VZFVVXTrA+iRJkha12dyz9svAy6vqm31tn0zyKuBTgGFNkiRpQGYzDZo9ghoAVfWNAdQjSZKkPrMJaw8lOWHPxq7t4dl0kmRJkn9I8pfd+nOSXJvknu772X37XpTk3iR3J1nb135Sktu6be/Nrt8OkSRJWsBmE9beAVyd5DeS/FSS1yT5TXq/u/Yrs+znbcBdfevvBK6rqmOB67p1kryQ3g/wvgg4HXhfkiXdMe8HLgCO7T6nz7JvSZKkg9Zs3g36KXqvl3oa8B+A/9gtn9Jt26ckq4CfBD7U17wOWN8trwfO6Gu/oqqmu9db3QucnOQI4NCquqGqCri87xhJkqQFazZPg35PVX0Z+PX97OPd9J4k/a6+thVVtQOgqnYkeV7XfiTw6b79tndtj3TLe7bPVO8F9K7AsWLFCsbHx/ez7NmZnJwceB+aG8ekTY5LW57xwAM8Y+tW/v6BB/jW4YcPuxz18d+V9gx7TGbzNOhVwEsBknykql4325MneQ1wf1XdkmTNbA6Zoa320f7Ext5PiVwKsHr16lqzZjbd7r/x8XEG3YfmxjFpk+PSiMlJOOcc2LKFR5csYeljj8Fpp8HGjbDcl9S0wH9X2jPsMZnV06B9y8+f4/lfBrw2yRfpvQD+FUk2AF/rpjbpvu/v9t8OHNV3/CpgomtfNUO7JGkuuqDG1BRLd+6Eqane+rnnDrsySXsxm7BWe1l+8gOrLqqqVVV1NL0HBz5ZVecCVwPnd7udT+9hBbr2s5OMJDmG3oMEN3VTpg8nOaV7CvS8vmMkSbNx333fCWq7mZqCa66BCf8bWGrRbKZBT0jyEL0rbId0y3TrVVWH7ke/7wI2JXkj8GXgLHonuyPJJuBO4FHgwqp6rDvmLcBlwCHAx7uPJGm2tm2DkZEnhjXotW/dCitXzn9dkvZpNi9yX/Jk+8xGVY0D493yN4BT97LfJcAlM7TfDBx3IGqRpEVpbAymp2feNj3d2y6pObN6kbskaQFYubL3MMHo6O7to6Owdq1X1aRGGdYkaTHZuLEXzEZHeWTZsseD2oYNw65M0l7M5p41SdJCsXw5XHUVTExw++bNnHjmmV5RkxrnlTVJWoxWruTB4483qEkHAcOaJElSwwxrkiRJDTOsSZIkNcywJkmS1DDDmiRJUsMMa5IkSQ0zrEmSJDXMsCZJktQww5okSVLDDGuSJEkNM6xJkiQ1zLAmSZLUMMOaJElSwwxrkiRJDTOsSZIkNcywJkmS1DDDmiRJUsMMa5IkSQ0zrEmSJDXMsCZJktQww5qkwZuY4LBbb4WJiWFXIkkHHcOapMGZnIR162BsjBdffDGMjfXWJyeHXZkkHTQMa5IG55xzYMsWmJpi6c6dMDXVWz/33GFXJkkHDcOapMG4777vBLXdTE3BNdc4JSpJs2RYkzQY27bByMjM20ZGYOvW+a1Hkg5ShjVJgzE2BtPTM2+bnu5tlyQ9KcOapMFYuRJOOw1GR3dvHx2FtWt72yVJT8qwJmlwNm7sBbPRUR5ZtuzxoLZhw7Ark6SDxtJhFyBpAVu+HK66CiYmuH3zZk4880yvqEnSHHllTdLgrVzJg8cfb1CTpP1gWJMkSWqYYU2SJKlhhjVJkqSGGdYkSZIaZliTJElqmGFNkiSpYQMNa0lGk9yU5B+T3JHkN7v25yS5Nsk93fez+465KMm9Se5Osrav/aQkt3Xb3pskg6xdkiSpBYO+sjYNvKKqTgBeApye5BTgncB1VXUscF23TpIXAmcDLwJOB96XZEl3rvcDFwDHdp/TB1y7JEnS0A00rFXPZLf69O5TwDpgfde+HjijW14HXFFV01X1BeBe4OQkRwCHVtUNVVXA5X3HSJIkLVgDv2ctyZIknwPuB66tqhuBFVW1A6D7fl63+5HAV/oO3961Hdkt79kuSZK0oA383aBV9RjwkiTPAq5Mctw+dp/pPrTaR/sTT5BcQG+6lBUrVjA+Pj63gudocnJy4H1obhyTNjku7XFM2uS4tGfYYzJvL3Kvqn9OMk7vXrOvJTmiqnZ0U5z3d7ttB47qO2wVMNG1r5qhfaZ+LgUuBVi9enWtWbPmQP4ZTzA+Ps6g+9DcOCZtclza45i0yXFpz7DHZNBPgz63u6JGkkOAnwA+D1wNnN/tdj7w0W75auDsJCNJjqH3IMFN3VTpw0lO6Z4CPa/vGEmSpAVr0FfWjgDWd090Pg3YVFV/meQGYFOSNwJfBs4CqKo7kmwC7gQeBS7splEB3gJcBhwCfLz7SJIkLWgDDWtVdStw4gzt3wBO3csxlwCXzNB+M7Cv+90kSZIWHN9gIEmS1DDDmiRJUsMMa5IkSQ0zrEmSJDXMsCZJktQww5okSVLDDGuSJEkNM6xJkiQ1zLAmSZLUMMOaJElSwwxrkiRJDTOsSZIkNcywJkmS1DDDmiRJUsMMa5IkSQ0zrGlhmZjgsFtvhYmJYVciSdIBYVjTwjA5CevWwdgYL774Yhgb661PTg67MkmSnhLDmhaGc86BLVtgaoqlO3fC1FRv/dxzh12ZJElPiWFNB7/77vtOUNvN1BRcc41TopKkg5phTQe/bdtgZGTmbSMjsHXr/NYjSdIBZFjTwW9sDKanZ942Pd3bLknSQcqwpoPfypVw2mkwOrp7++gorF3b2y5J0kHKsKaFYePGXjAbHeWRZcseD2obNgy7MkmSnpKlwy5AOiCWL4erroKJCW7fvJkTzzzTK2qSpAXBK2taWFau5MHjjzeoSZIWDMOaJElSwwxrkiRJDTOsSZIkNcywJkmS1DDDmiRJUsMMa5IkSQ0zrEmSJDXMsCZJktQww5okSVLDDGuSJEkNM6xJkiQ1zLAmSZLUMMOaJElSwwxrkiRJDTOsSZIkNWygYS3JUUmuT3JXkjuSvK1rf06Sa5Pc030/u++Yi5Lcm+TuJGv72k9Kclu37b1JMsjaJUmSWjDoK2uPAu+oqh8ETgEuTPJC4J3AdVV1LHBdt0637WzgRcDpwPuSLOnO9X7gAuDY7nP6gGuXJEkauoGGtaraUVWf7ZYfBu4CjgTWAeu73dYDZ3TL64Arqmq6qr4A3AucnOQI4NCquqGqCri87xhJkqQFa97uWUtyNHAicCOwoqp2QC/QAc/rdjsS+ErfYdu7tiO75T3bJUmSFrSl89FJkuXAR4C3V9VD+7jdbKYNtY/2mfq6gN50KStWrGB8fHzO9c7F5OTkwPvQ3DgmbXJc2uOYtMlxac+wx2TgYS3J0+kFtY1V9Rdd89eSHFFVO7opzvu79u3AUX2HrwImuvZVM7Q/QVVdClwKsHr16lqzZs2B+lNmND4+zqD70Nw4Jm1yXNrjmLTJcWnPsMdk0E+DBvgwcFdV/X7fpquB87vl84GP9rWfnWQkyTH0HiS4qZsqfTjJKd05z+s7RpIkacEa9JW1lwFvAG5L8rmu7WLgXcCmJG8EvgycBVBVdyTZBNxJ70nSC6vqse64twCXAYcAH+8+kiRJC9pAw1pVfYqZ7zcDOHUvx1wCXDJD+83AcQeuOkmSpPb5BgNJkqSGGdYkSZIaZliTJElqmGFNkiSpYYY1SZKkhhnWJEmSGmZYkyRJaphhTZIkqWGGNUmSpIYZ1iRJkhpmWJMkSWqYYU2SJKlhhjVJkqSGGdYkSZIaZliTJElqmGFNkiSpYYY1SZKkhhnWJEmSGmZYkyRJaphhTZIkqWGGNUmSpIYZ1p6KiQkOu/VWmJgYdiWSJGmBMqztj8lJWLcOxsZ48cUXw9hYb31yctiVSZKkBcawtj/OOQe2bIGpKZbu3AlTU731c88ddmWSJGmBMazN1X33fSeo7WZqCq65xilRSZJ0QBnW5mrbNhgZmXnbyAhs3Tq/9UiSpAXNsDZXY2MwPT3ztunp3nZJkqQDxLA2VytXwmmnwejo7u2jo7B2bW+7JEnSAWJY2x8bN/aC2egojyxb9nhQ27Bh2JVJkqQFZumwCzgoLV8OV10FExPcvnkzJ555plfUJEnSQHhl7alYuZIHjz/eoCZJkgbGsCZJktQww5okSVLDDGuSJEkNM6xJkiQ1zLAmSZLUMMOaJElSwwxrkiRJDTOsSZIkNcywJkmS1LBU1bBrGJgkXwe+NOBuDgceGHAfmhvHpE2OS3sckzY5Lu2ZrzH53qp67p6NCzqszYckN1fV6mHXocc5Jm1yXNrjmLTJcWnPsMfEaVBJkqSGGdYkSZIaZlh76i4ddgF6AsekTY5LexyTNjku7RnqmHjPmiRJUsO8siZJktQww9p+SnJUkuuT3JXkjiRvG3ZNi12S0SQ3JfnHbkx+c9g1qSfJkiT/kOQvh12LepJ8McltST6X5OZh1yNI8qwkm5N8vvv/lh8Zdk2LXZIXdP+O7Po8lOTt816H06D7J8kRwBFV9dkk3wXcApxRVXcOubRFK0mAZVU1meTpwKeAt1XVp4dc2qKX5FeA1cChVfWaYdejXlgDVleVv+fViCTrgb+tqg8leQbwzKr652HXpZ4kS4D7gB+uqkH/hutuvLK2n6pqR1V9tlt+GLgLOHK4VS1u1TPZrT69+/hfI0OWZBXwk8CHhl2L1KokhwI/CnwYoKq+ZVBrzqnA1vkOamBYOyCSHA2cCNw43ErUTbd9DrgfuLaqHJPhezfwa8C3h12IdlPAliS3JLlg2MWI5wNfB/64u2XgQ0mWDbso7eZs4M+G0bFh7SlKshz4CPD2qnpo2PUsdlX1WFW9BFgFnJzkuGHXtJgleQ1wf1XdMuxa9AQvq6qXAq8CLkzyo8MuaJFbCrwUeH9VnQjsBN453JK0Szct/Vrgz4fRv2HtKejui/oIsLGq/mLY9ehx3fTBOHD6kEtZ7F4GvLa7P+oK4BVJNgy3JAFU1UT3fT9wJXDycCta9LYD2/tmAzbTC29qw6uAz1bV14bRuWFtP3U3s38YuKuqfn/Y9QiSPDfJs7rlQ4CfAD4/3KoWt6q6qKpWVdXR9KYQPllV5w65rEUvybLuwSi6qbbTgNuHW9XiVlVfBb6S5AVd06mAD6y14/UMaQoUepddtX9eBrwBuK27Rwrg4qr62BBrWuyOANZ3T+w8DdhUVf5UhPREK4Are//NyVLgT6vqE8MtScBbgY3dlNs24OeGXI+AJM8EXgm8eWg1+NMdkiRJ7XIaVJIkqWGGNUmSpIYZ1iRJkhpmWJMkSWqYYU2SJKlhhjVJzUpSSX6vb/1Xk/zGATr3ZUnOPBDnepJ+zkpyV5Lr92hfk2TGn5ZJ8rFdvxkoSYY1SS2bBmZpjnQAAAL7SURBVH46yeHDLqRf91t+s/VG4Beq6sdne0BVvdqXeEvaxbAmqWWPApcCv7znhj2vjCWZ7L7XJPnrJJuS/FOSdyU5J8lNSW5LMtZ3mp9I8rfdfq/pjl+S5HeTfCbJrUne3Hfe65P8KXDbDPW8vjv/7Ul+p2v7deDlwAeS/O4Mf9+hSa5McmeSDyR5WnfcF5McnuTo7qrcB5PckWRL93YOkvxSd9ytSa7Yr3+6kg4KvsFAUuv+ELg1yf+YwzEnAD8IfJPeL8F/qKpOTvI2er8S//Zuv6OBHwPGgOuTfB9wHvBgVf1QkhHg75Js6fY/GTiuqr7Q31mSlcDvACcB/w/YkuSMqvqtJK8AfrWqbp6hzpOBFwJfAj4B/DS9d0L2OxZ4fVW9Kckm4HXABnov+T6mqqadMpUWNq+sSWpaVT0EXA780hwO+0xV7aiqaWArsCts3UYvoO2yqaq+XVX30At1P0DvPZnnda+RuxH4bnqBCeCmPYNa54eA8ar6elU9CmwEfnQWdd5UVduq6jF67x18+Qz7fKGqdr3S7pa++m+l92qic+ldgZS0QBnWJB0M3k3v3q9lfW2P0v1vWHovuXxG37bpvuVv961/m91nFPZ8314BAd5aVS/pPsdU1a6wt3Mv9WW2f8gM/e1rHXb/Wx7j8fp/kt5Vx5OAW5I4UyItUIY1Sc2rqm8Cm+gFtl2+SC+oAKwDnr4fpz4rydO6+9ieD9wNXAO8JcnTAZJ8f5Jl+zoJvStwP9bdZ7YEeD3w17Po/+Qkx3T3qv0M8KnZFN3tf1RVXQ/8GvAsYPlsjpV08PG/xCQdLH4P+MW+9Q8CH01yE3Ade7/qtS930wtVK4Cfr6qpJB+iN9X42e6K3deBM/Z1kqrakeQi4Hp6V9k+VlUfnUX/NwDvAl4M/A1w5SzrXgJsSHJY198f+PSotHClaqar7pIkSWqB06CSJEkNM6xJkiQ1zLAmSZLUMMOaJElSwwxrkiRJDTOsSZIkNcywJkmS1DDDmiRJUsP+P6Vk9DX1oXueAAAAAElFTkSuQmCC\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAl8AAAGDCAYAAAAVq3XUAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3de5hddX3v8fcXghnNCAEvaaLY6IjWKjeZUlovTESTorShClWbKFVqPK0g9HJazKnV9jxt6eljbc+xN6qtaRNNMUqgPtbERsbWHgETtAG8lJMIaAYBkYsTmJHL9/yx1sAmTMieZNZvr8y8X8+zn73X+u21ft/k96Cf/H5rrxWZiSRJkso4pNcFSJIkzSaGL0mSpIIMX5IkSQUZviRJkgoyfEmSJBVk+JIkSSrI8CVJklSQ4UuSJKkgw5ekaRMRN0XE/REx2vFaFBFzI+KPIuKWuv3GiPjvEREdxw5HREbE8Xucc2O9f2gffb8/Ih6IiB/Ur/+KiA9FxMKO7wxFxMN71DcaET8VETd0bD8UEWMd26sj4pfq/Z3HfWgfNV0YETsj4t6IGImID0bEnC7+HudGxEci4ub6z/KViDh9X8dJOjgYviRNt5/NzP6O1wjwCeA04LXAU4G3AKuAP9/j2P8C3jqxERFPA04B7uiy73/KzKcCRwE/D/wIsK0zgAEje9TXn5lfyswXT2wD/w6c19H+h/WxX9rjuPP2Uc8/Ay/NzMOBlwDHA+/u4s8xB/g2cCpwBPBe4NKIWNzV34KkVjN8SWpURJwGLAXekJnXZ+aDmXkVsBJ4V0Q8v+Pr64A3RsSh9fabgcuAH06lz8x8IDNvAN5IFdx+40D/HPsjM3dk5t31ZgAPA88HiIifjojvRcTR9fbxEXF3RPxYZu7OzPdn5k2Z+XBmfhr4FnBSL/4ckqaX4UtS014DXJ2Z3+7cmZlXA9+hmhGbMAJ8jSqsQTUL9g/723FmPgRcDrxif89xoCLiFyPiXuB7VDNff1PX9n/rz2si4snAPwK/k5nfmOQcC4AXADcUK1xSYwxfkqbbxnoG5+6I2Ag8Hbh1L9+9tW7v9A/AWyPihcD8zPzSAdYzQrUMOWFRR30Tr3ldnuuUPY47ZV8HZObH6mXHFwB/DdzW0fx+qmXFa+o6/2LP4yPiMKoZwTWTBTNJBx/Dl6TpdmZmzq9fZ1LN+Czcy3cX1u2dPgW8CjifajboQD0L+H7H9khHfROv3V2e66o9jruq2yIy80aqmau/7Nj3APBRquvBPpCZ2XlMRBxC9XfwQ2Bf15dJOkgYviQ17V+Bn5y4tmlCRJwMHA18vnN/Zt4H/AvwKxxg+KrDy89SXUDfBnOAgYmNiHgW8D7g74EPRMTcjrYAPgIsoLpe7oHCtUpqiOFLUqMy81+BLcAnI+LFEXFovVy3DvirekZoT6uBUzPzpv3pMyIOi4gXAR+n+sXjn+5f9QcmIn45Ip5Zf/5x4D1UfxcT4eqjVAHrXKol2P/ZcfhfAS+i+vXo/QXLltQww5ekEt4AXAl8FhgF1lKFjvMn+3JmjmTmF/ejnzdGxChwN3AFcCdwUn27iwmLJrnP1xv2o69uvAy4LiJ2A5+pX6vrtndTzWq9t15ufBvwtoh4RUT8KPBO4ATgux11rmioTkkFxR6XGEiSJKlBznxJkiQVZPiSdNCIiH+ZZMlwNCJW7/voxmq6YS81uUQoaVIuO0qSJBXkzJckSVJBc3pdQLee/vSn5+LFixvtY/fu3cyb1+2NrlWK49I+jkk7OS7t45i0T6kx2bZt2/cy8xmTtR004Wvx4sVs3bq10T6Gh4cZGhpqtA9NnePSPo5JOzku7eOYtE+pMYmIm/fW5rKjJElSQYYvSZKkggxfkiRJBRm+JEmSCjJ8SZIkFWT4kiRJKsjwJUmSVJDhS5IkqSDDlyRJmhVGRmD79iMYGeltHYYvSZI0o42OwvLlMDAAq1cfy8BAtT062pt6DF+SJGlGW7ECNm+GsTHYvXsOY2PV9sqVvanH8CVJkmasXbseDV6dxsZg0yZ6sgRp+JIkSTPWzp0wd+7kbXPnwo4dZesBw5ckSZrBBgZgfHzytvHxqr00w5ckSZqxFi2CpUuhr++x+/v6YNmyqr00w5ckSZrR1q2rglZfH8yb98AjwWvt2t7UM6c33UqSJJXR3w8bN1YX12/YcD1nnXViT2a8JjjzJUmSZoVFi+C44+7pafACw5ckSVJRhi9JkqSCDF+SJEkFGb4kSZIKMnxJkiQVZPiSJEkqyPAlSZJUkOFLkiSpIMOXJElSQYYvSZKkggxfkiRJBRm+JEmSCjJ8SZIkFWT4kiRJKsjwJUmSVFDj4SsiLoiI6yPihoi4sN53VER8LiJurN+PbLoOSZKkNmg0fEXES4B3ACcDxwNnRMQxwEXAlsw8BthSb0uSJM14Tc98vQi4KjPvy8wHgS8APw8sB9bU31kDnNlwHZIkSa3QdPi6HnhlRDwtIp4CvBY4GliQmbcC1O/PbLgOSZKkVojMbLaDiHOBdwGjwNeA+4G3Zeb8ju/clZmPu+4rIlYBqwAWLFhw0vr16xutdXR0lP7+/kb70NQ5Lu3jmLST49I+jkn7lBqTJUuWbMvMwcnaGg9fj+ks4g+B7wAXAEOZeWtELASGM/OFT3Ts4OBgbt26tdH6hoeHGRoaarQPTZ3j0j6OSTs5Lu3jmLRPqTGJiL2GrxK/dnxm/f4c4PXAx4ErgHPqr5wDXN50HZIkSW0wp0Afn4yIpwEPAO/KzLsi4mLg0npJ8hbg7AJ1SJIk9Vzj4SszXzHJvjuB05ruW5IkqW28w70kSVJBhi9JkqSCDF+SJEkFGb4kSZIKMnxJkiQVZPiSJEkqyPAlSZJUkOFLkiSpIMOXJElSQYYvSZKkggxfkiRJBRm+JEmSCjJ8SZIkFWT4kiRJKsjwJUmSVJDhS5IkqSDDlyRJUkGGL0mSpIIMX5IkSQUZviRJkgoyfEmSJBVk+JIkSSrI8CVJklSQ4UuSJKkgw5ckSVJBhi9JkqSCDF+SJEkFGb4kSZIKMnxJkiQVZPiSJEkqyPAlSZJUkOFLkiSpIMOXJElSQY2Hr4j4tYi4ISKuj4iPR0RfRBwVEZ+LiBvr9yObrkOSJKkNGg1fEfEs4N3AYGa+BDgUeBNwEbAlM48BttTbkiRJM16JZcc5wJMjYg7wFGAEWA6sqdvXAGcWqEOSJKnnIjOb7SDiAuAPgPuBzZm5IiLuzsz5Hd+5KzMft/QYEauAVQALFiw4af369Y3WOjo6Sn9/f6N9aOocl/ZxTNrJcWkfx6R9So3JkiVLtmXm4GRtc5rsuL6WaznwXOBu4BMRsbLb4zPzEuASgMHBwRwaGmqizEcMDw/TdB+aOselfRyTdnJc2scxaZ82jEnTy46vBr6VmXdk5gPAp4CfBm6LiIUA9fvtDdchSZLUCk2Hr1uAUyLiKRERwGnA14ErgHPq75wDXN5wHZIkSa3Q6LJjZl4dERuAa4EHga9QLSP2A5dGxLlUAe3sJuuQJElqi0bDF0Bmvg943x67x6lmwSRJkmYV73AvSZJUkOFLkiSpIMOXJEkNGBmB7duPYGSk15WobQxfkiRNo9FRWL4cBgZg9epjGRiotkdHe12Z2sLwJUnSNFqxAjZvhrEx2L17DmNj1fbKrm8xrpnO8CVJ0jTZtevR4NVpbAw2bcIlSAGGL0mSps3OnTB37uRtc+fCjh1l61E7Gb4kSZomAwMwPj552/h41S4ZviRJmiaLFsHSpdDX99j9fX2wbFnVLhm+JEmaRuvWVUGrrw/mzXvgkeC1dm2vK1NbNP54IUmSZpP+fti4sbq4fsOG6znrrBOd8dJjOPMlSVIDFi2C4467x+ClxzF8SZIkFWT4kiRJKsjwJUmSVJDhS5IkqSDDlyRJUkGGL0mSpIIMX5IkSQUZviRJkgoyfEmSJBVk+JIkSSrI8CVJklSQ4UuSJKkgw5ckSVJBhi9JkqSCDF+SJEkFGb4kSZIKMnxJkiQVZPiSJEkqyPAlSZJUkOFLkiSpoEbDV0S8MCK+2vG6NyIujIijIuJzEXFj/X5kk3VIkiS1RaPhKzO/mZknZOYJwEnAfcBlwEXAlsw8BthSb0uSJM14JZcdTwN2ZObNwHJgTb1/DXBmwTokSZJ6JjKzTEcRfwdcm5kfioi7M3N+R9tdmfm4pceIWAWsAliwYMFJ69evb7TG0dFR+vv7G+1DU+e4tI9j0k6OS/s4Ju1TakyWLFmyLTMHJ2srEr4i4knACPDizLyt2/DVaXBwMLdu3dponcPDwwwNDTXah6bOcWkfx6SdHJf2cUzap9SYRMRew1epZcfTqWa9bqu3b4uIhXVxC4HbC9UhSZLUU6XC15uBj3dsXwGcU38+B7i8UB2SJEk91Xj4ioinAK8BPtWx+2LgNRFxY912cdN1SJIktcGcpjvIzPuAp+2x706qXz9KkiTNKt7hXpIkqaCuwldELIuIcyNi8R77395EUZIkSTPVPsNXRPwh8D+AY4EtEXF+R/N5TRUmSZI0E3Uz8/WzwKsy80KqRwSdHhEfrNuiscokSZJmoG7C15zMfBAgM++mCmOHR8QngCc1WZwkSdJM00342hERp05sZOZDmXku8E3gRY1VJkmSNAN1E77OBq7Zc2dm/g5w9LRXJEmSNIPtM3xl5v2ZeX9EnNu5PyIOBX65scokSZJmoKnc5+u0iPhMRCyMiJcAVwFPbaguSZKkGanrO9xn5i9GxBuB64D7gDdn5n80VpkkSdIM1PXMV0QcA1wAfBK4CXhL/dxGSZIkdWkqy47/DLw3M98JnArcCHy5kaokSZJmqKk8WPvkzLwXIDMT+EBEXNFMWZIkSTNTV+ErIn6k/nhvRDwDeAXwzcy8obHKJEmSZqBunu34TuBLwFUR8SvAp4EzgE/tefsJSZIkPbFuZr7OA14MPBm4GXh+Zn43Io4ErgQ+0mB9kiRJM0o34euBzLwPuC8idmTmdwEy866IyGbLkyRJmlm6+bXjwxFxWP35dRM7I6Kvy+MlSZJU6yY8vR5IgMz8Tsf+pwG/0URRkiRJM1U3z3a8JTMfnGT/rsz814ntiPjSdBcnSerOyAhs334EIyO9rkTSvkznsmHfNJ5LktSF0VFYvhwGBmD16mMZGKi2R0d7XZmkvZnO8OXF95JU2IoVsHkzjI3B7t1zGBurtleu7HVlkvbGC+Yl6SC1a9ejwavT2Bhs2oRLkFJLTWf4imk8lyRpH3buhLlzJ2+bOxd27Chbj6TudHOH+1O6PNdbDrAWSdIUDAzA+PjkbePjVbuk9ulm5usvJz480S8aM/P6aalIktSVRYtg6VLo2+PnTn19sGxZ1S6pfboJX53Lif6iUZJaZN26Kmj19cG8eQ88ErzWru11ZZL2ppvHCx1SP8fxkI7PjwSyzPx+U8VJkp5Yfz9s3FhdXL9hw/WcddaJznhJLddN+DoC2MajgevajrYEnjfdRUmSpmbRIjjuuHsMXtJBoJvwdWpm3tx4JZIkSbNAN9d8XdZ4FZIkSbPEVC+4lyRJ0gHoZtnxWRHxv/fWmJnvfqKDI2I+8GHgJVTXiL0d+CbwT8Bi4CbgFzLzru5KliRJOnh1E77up7rgfn/9OfDZzDwrIp4EPAVYDWzJzIsj4iLgIuC3D6APSZKkg0I34evOzFyzPyePiMOBVwK/BJCZPwR+GBHLgaH6a2uAYQxfkiRpFojMfOIvRFyVmd0+YmjPY08ALgG+BhxPNYN2AbArM+d3fO+uzDxykuNXAasAFixYcNL69ev3p4yujY6O0t/f32gfmjrHpX0ck3ZyXNrHMWmfUmOyZMmSbZk5OFlbNzNfH5r4EBEvy8z/6Ng+LzM/NPlhj5z/pcD5mXl1RPw51RJjVzLzEqrwxuDgYA4NDXV76H4ZHh6m6T40dY5L+zgm7eS4tI9j0j5tGJNufu346x2f/88ebW/fx7HfAb6TmVfX2xuowthtEbEQoH6/vYs6JEmSDnpTvdXEnredeMLbUGTmd4FvR8QL612nUS1BXgGcU+87B7i8izokSZIOet0sO+ZePk+2PZnzgXX1Lx13Am+jCn2XRsS5wC3A2V2cR5Ik6aDXTfj6sYjYTjXLNVB/pt7e53MdM/OrwGQXnJ3WdZWSJEkzRDfh60WNVyHpoDEyAtu3H8ELXoAPcZak/bDPa74y8+b6wdr3AM+sX3d37Jc0C4yOwvLlMDAAq1cfy8BAtT062uvKJOngss/wFRFPioiPUj0G6BLgb4GbIuLv6uu4JM0CK1bA5s0wNga7d89hbKzaXrmy15VJ0sGlm187/g5wGHB0Zp6YmScAz6Fasnxvk8VJaoddux4NXp3GxmDTpmopUpLUnW7C1+uBd2TmDyZ21J9/Ffj5pgqT1B47d8LcuZO3zZ0LO3aUrUeSDmbdhK+HM/O+PXdm5ijd3WpC0kFuYADGxydvGx+v2iVJ3ekmfGVEHBkRR+35Ah5uukBJvbdoESxdCn19j93f1wfLlvmrR0maim5uNXEE1QOxJ7ubvTNf0iyxbl11cf2mTXDooQ/w0EOHsWwZrF3b68ok6eCyz/CVmYu7OVFEvDgzbzjgiiS1Un8/bNxYXVy/YcP1nHXWic54SdJ+6GbZsVv/OI3nktRSixbBccfdY/CSpP00neHrCR+yLUmSpOkNX17/JUmStA/TGb4kSZK0D9MZvn44jeeSJEmakfb5a8eIeOkTtWfmtfX7KdNVlCRJ0kzVzX2+tgI3AHfU250X1ifwqukuSpIkaabqJnz9BvAG4H5gPXBZ/WghSZIkTdE+r/nKzA9m5suB84CjgS0RcWlEnNB4dZIkSTNM1xfcZ+a3gMuBzcDJwAuaKkqSJGmm6uaC++cBbwKWA9+mWnr8g8wca7g2SZKkGaeba77+H7CdatbrXuA5wK9GVNfdZ+afNladJEnSDNNN+Pp9Hr17fX+DtUiSJM14+wxfmfn+AnVIkiTNCl1dcB8Rp0fEv0XE9yLijoj4QkS8tuniJEmSZppuLrh/B/BO4LeobrgKMAhcHBHPzsxLGqxPkiRpRunmmq9fA16emd/v2Pf5iDgd+CJg+JIkSepSN8uOsUfwAiAz72ygHkmSpBmtm/B1b0Qcv+fOet8Ppr8kSZKkmavbZzteERF/D2yjuu3ETwDnACsbrE2SJGnG6ebZjl+kepzQIcAvAW+vP59St0mSJKlL3fza8TmZeQvwuwXqkSRJmtG6ueZr48SHiPhkg7VIkiTNeN1c8xUdn5831Q4i4iaqC/MfAh7MzMGIOAr4J2AxcBPwC5l511TPLUmSdLDpZuYr9/J5KpZk5gmZOVhvXwRsycxjgC31tiRJ0ozXzczX8RFxL9UM2JPrz9TbmZmH70e/y4Gh+vMaYBj47f04jyRJ0kElMvd3MqvLDiK+BdxFNWv2N5l5SUTcnZnzO75zV2YeOcmxq4BVAAsWLDhp/fr1jdY6OjpKf39/o31o6hyX9nFM2slxaR/HpH1KjcmSJUu2daz4PUY3M18H6mWZORIRzwQ+FxHf6PbA+rmRlwAMDg7m0NBQQyVWhoeHaboPTZ3j0j6OSTs5Lu3jmLRPG8akm2u+DkhmjtTvtwOXUd0z7LaIWAhQv9/edB2SJElt0Gj4ioh5EfHUic/AUuB64AqqO+RTv1/eZB2SJElt0fSy4wLgsoiY6OtjmfnZiPgycGlEnAvcApzdcB2SJEmt0Gj4ysydwOMeyp2ZdwKnNdm3JElSGzV+zZckSZIeZfiSJEkqyPAlSZJUkOFLkiSpIMOXJElSQYYvSZKkggxfkiRJBRm+JEmSCjJ8SZIkFWT4kiRJKsjwJUmSVJDhS5IkqSDDlyRJUkGGL0mSpIIMX5IkSQUZviRJkgoyfEmSJBVk+JIkSSrI8CVJklSQ4UuSJKkgw5ckSVJBhi9JkqSCDF+SJEkFGb4kSZIKMnxJkiQVZPiSJEkqyPAlSZJUkOFLkiSpIMOXJElSQYYvSZKkggxfkiRJBRm+JEmSCioSviLi0Ij4SkR8ut4+KiI+FxE31u9HlqhDkiSp10rNfF0AfL1j+yJgS2YeA2yptyVJkma8xsNXRDwbeB3w4Y7dy4E19ec1wJlN1yFJktQGkZnNdhCxAfgj4KnAb2bmGRFxd2bO7/jOXZn5uKXHiFgFrAJYsGDBSevXr2+01tHRUfr7+xvtQ1PnuLSPY9JOjkv7OCbtU2pMlixZsi0zBydrm9NkxxFxBnB7Zm6LiKGpHp+ZlwCXAAwODubQ0JRPMSXDw8M03YemznFpH8eknRyX9nFM2qcNY9Jo+AJeBvxcRLwW6AMOj4i1wG0RsTAzb42IhcDtDdchSZLUCo1e85WZ78nMZ2fmYuBNwOczcyVwBXBO/bVzgMubrEOSJKktenWfr4uB10TEjcBr6m3pcUZGYPv2IxgZ6XUlkiRNj2LhKzOHM/OM+vOdmXlaZh5Tv3+/VB06OIyOwvLlMDAAq1cfy8BAtT062uvKJEk6MN7hXq20YgVs3gxjY7B79xzGxqrtlSt7XZkkSQfG8KXW2bXr0eDVaWwMNm3CJUhJ0kHN8KXW2bkT5s6dvG3uXNixo2w9kiRNJ8OXWmdgAMbHJ28bH6/aJUk6WBm+1DqLFsHSpdDX99j9fX2wbFnVLknSwcrwpVZat64KWn19MG/eA48Er7Vre12ZJEkHpuk73Ev7pb8fNm6sLq7fsOF6zjrrRGe8JEkzgjNfarVFi+C44+4xeEmSZgzDlyRJUkGGL0mSpIIMX5IkSQUZviRJkgoyfEmSJBVk+JIkSSrI8CVJklSQ4UuSJKkgw5ckSVJBhi9JkqSCDF+SJEkFGb4kSZIKMnxJkiQVZPiSJEkqyPAlSZJUkOFLkiSpIMOXJElSQYYvSZKkggxfkiRJBRm+JEmSCjJ8SZIkFWT4kiRJKsjwJUmSVJDhS5IkqaBGw1dE9EXENRHxnxFxQ0T8Xr3/qIj4XETcWL8f2WQdkiRJbdH0zNc48KrMPB44AfiZiDgFuAjYkpnHAFvqbUmSpBmv0fCVldF687D6lcByYE29fw1wZpN1SJIktUVkZrMdRBwKbAOeD/xFZv52RNydmfM7vnNXZj5u6TEiVgGrABYsWHDS+vXrG611dHSU/v7+RvvQ1Dku7eOYtJPj0j6OSfuUGpMlS5Zsy8zBydoaD1+PdBQxH7gMOB/4Yjfhq9Pg4GBu3bq10RqHh4cZGhpqtA9NnePSPo5JOzku7eOYtE+pMYmIvYavYr92zMy7gWHgZ4DbImJhXdxC4PZSdUiSJPVS0792fEY940VEPBl4NfAN4ArgnPpr5wCXN1mHJElSW8xp+PwLgTX1dV+HAJdm5qcj4kvApRFxLnALcHbDdUiSJLVCo+ErM7cDJ06y/07gtCb7liRJaiPvcC9JklSQ4UuSJKkgw5ckSVJBhi9JkqSCDF+SJEkFGb4kSZIKMnxJkiQVZPiSJEkqyPAlSZJUkOFLkiSpIMOXJElSQYYvSZKkggxfkiRJBRm+JEmSCjJ8SZIkFWT4kiRJKsjwVRsZge3bj2BkpNeVSJKkmWzWh6/RUVi+HAYGYPXqYxkYqLZHR3tdmSRJmolmffhasQI2b4axMdi9ew5jY9X2ypW9rkySJM1Eszp87dr1aPDqNDYGmzbhEqQkSZp2szp87dwJc+dO3jZ3LuzYUbYeSZI0883q8DUwAOPjk7eNj1ftkiRJ02lWh69Fi2DpUujre+z+vj5YtqxqlyRJmk6zOnwBrFtXBa2+Ppg374FHgtfatb2uTJIkzURzel1Ar/X3w8aN1cX1GzZcz1lnneiMlyRJasysn/masGgRHHfcPQYvSZLUKMOXJElSQYYvSZKkggxfkiRJBRm+JEmSCjJ8SZIkFWT4kiRJKsjwJUmSVJDhS5IkqSDDlyRJUkGRmb2uoSsRcQdwc8PdPB34XsN9aOocl/ZxTNrJcWkfx6R9So3Jj2bmMyZrOGjCVwkRsTUzB3tdhx7LcWkfx6SdHJf2cUzapw1j4rKjJElSQYYvSZKkggxfj3VJrwvQpByX9nFM2slxaR/HpH16PiZe8yVJklSQM1+SJEkFGb6AiDg6Iq6MiK9HxA0RcUGva5rtIqIvIq6JiP+sx+T3el2TKhFxaER8JSI+3etaVImImyLiuoj4akRs7XU9qkTE/IjYEBHfqP//5ad6XdNsFhEvrP8bmXjdGxEX9qQWlx0hIhYCCzPz2oh4KrANODMzv9bj0matiAhgXmaORsRhwBeBCzLzqh6XNutFxK8Dg8DhmXlGr+tRFb6Awcz0flItEhFrgH/PzA9HxJOAp2Tm3b2uS9U/IoFdwE9mZtP3EH0cZ76AzLw1M6+tP/8A+DrwrN5WNbtlZbTePKx++S+FHouIZwOvAz7c61qkNouIw4FXAh8ByMwfGrxa5TRgRy+CFxi+HiciFgMnAlf3thLVy1tfBW4HPpeZjknv/RnwW8DDvS5Ej5HA5ojYFhGrel2MAHgecAfw9/Uy/YcjYl6vi9Ij3gR8vFedG746REQ/8Engwsy8t9f1zHaZ+VBmngA8Gzg5Il7S65pms4g4A7g9M7f1uhY9zssy86XA6cC7IuKVvS5IzAFeCvxVZp4I7AYu6m1JAqiXgH8O+ESvajB81errij4JrMvMT/W6Hj2qnqofBn6mx6XMdi8Dfq6+vmg98KqIWNvbkgSQmSP1++3AZcDJva1IwHeA73TM2G+gCmPqvdOBazPztl4VYPjikYu7PwJ8PTP/tNf1CCLiGRExv/78ZODVwDd6W9XslpnvycxnZ+Ziqin7z2fmyh6XNetFxLz6h0LUy1pLget7W5Uy87vAtyPihfWu0wB/xNUOb6aHS45QTYuq+hf9W4Dr6muMAFZn5md6WNNstxBYU/8i5RDg0sz01gbS4y0ALqv+Dckc4GOZ+dnelqTa+cC6eplrJ/C2Htcz60XEU4DXAO/saR3eakKSJKkclx0lSZIKMnxJkiQVZPiSJEkqyPAlSZJUkOFLkiSpIMOXpGIiIiPiAx3bvxkR75+mc380Is6ajnPto5+zI+LrEXHlHvuHImLS26FExGcm7lsnSYYvSSWNA6+PiKf3upBO9f3kugVQEHEAAAL1SURBVHUu8KuZuaTbAzLztT5UWdIEw5ekkh4ELgF+bc+GPWeuImK0fh+KiC9ExKUR8V8RcXFErIiIayLiuogY6DjNqyPi3+vvnVEff2hE/ElEfDkitkfEOzvOe2VEfAy4bpJ63lyf//qI+ON63+8CLwf+OiL+ZJI/3+ERcVlEfC0i/joiDqmPuykinh4Ri+tZs7+NiBsiYnP9BAci4t31cdsjYv1+/e1KOih4h3tJpf0FsD0i/tcUjjkeeBHwfao7hX84M0+OiAuo7iJ+Yf29xcCpwABwZUQ8H3grcE9m/kREzAX+IyI2198/GXhJZn6rs7OIWAT8MXAScBewOSLOzMzfj4hXAb+ZmVsnqfNk4MeBm4HPAq+neqZfp2OAN2fmOyLiUuANwFqqhy4/NzPHXaKUZjZnviQVlZn3Av8AvHsKh305M2/NzHFgBzARnq6jClwTLs3MhzPzRqqQ9mNUzzp8a/3osKuBp1EFIIBr9gxetZ8AhjPzjsx8EFgHvLKLOq/JzJ2Z+RDVs+NePsl3vpWZE48x29ZR/3aqR9GspJohlDRDGb4k9cKfUV07Na9j34PU/5tUP+z+SR1t4x2fH+7YfpjHzuDv+by0BAI4PzNPqF/PzcyJ8LZ7L/VFt3+QSfp7om147J/lIR6t/3VUs4InAdsiwpUJaYYyfEkqLjO/D1xKFcAm3EQVPACWA4ftx6nPjohD6uvAngd8E9gE/EpEHAYQES+IiHlPdBKqGbJT6+u0DgXeDHyhi/5Pjojn1td6vRH4YjdF198/OjOvBH4LmA/0d3OspIOP/7KS1CsfAM7r2P5b4PKIuAbYwt5npZ7IN6lC0gLgv2XmWER8mGpp79p6Ru0O4MwnOklm3hoR7wGupJoF+0xmXt5F/18CLgaOBf4NuKzLug8F1kbEEXV/H/TXkdLMFZmTzYpLkiSpCS47SpIkFWT4kiRJKsjwJUmSVJDhS5IkqSDDlyRJUkGGL0mSpIIMX5IkSQUZviRJkgr6/zZFBrnPhQVPAAAAAElFTkSuQmCC\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmsAAAGDCAYAAAB0s1eWAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3df5SeZX3n8ffXGZnRJJSmaJYhbIOPaItsqhJRj66N0jK0dR02xjYsFOt6SmWp1W5dLfGo7W6ztXVbu/GselJ1gSaSppFfx6oMZRl/tCAGqiGAIANak6FFG5BM6kzI+N0/njvpw/BMmCHzzH1l5v06Z848z3XdP77jddRPruv+EZmJJEmSyvSMuguQJEnS1AxrkiRJBTOsSZIkFcywJkmSVDDDmiRJUsEMa5IkSQUzrEmSJBXMsCZJklQww5qkeSUiXh0RfxcRP4iIvRHxtxHxsoj4tYjIiPhvk7bfHRGrq8+/FxGPR8RoRDxaHeeV0zjn6RGxIyIeqX7+JiJOn2a9b46I2yPisaqWP46I7qf1x0ualwxrkuaNiDge+CzwEWApcDLw+8B4tcle4D3VdlP5y8xcDJwI3Az81TROPQKsrc55InA9sHWaZT8beGe138uBs4F3TXNfSQuAYU3SfPICgMy8KjMnMvOHmTmYmTur/nuAW4DffqoDZeZBYAtwckQ85ym2fTQzv53N9/cFMAE8HyAijouIr0fE26vvXdVs3/urfT+WmV/OzAOZuac656uezh8vaX5yql3SfHIfMBERV9Cc2bo1Mx+ZtM37gKGI+Ehm7p3qQBFxHHAR8M/A5GNMtc+jwGKa/xA+FMYORMSFwJcj4m+ANUAXsGGKw7wGuGs655O0MDizJmneyMzHgFcDCfw58L2IuD4ilrVs83VgEHjPFIf55Sp0/RD4dWBtNcs2nfOfAPwY8JvA37e07wL+ALiG5hLnr2bmxOT9I+ItwCrgf03nfJIWBsOapHklM+/JzF/LzOXAGUAf8GeTNns/cElE/Js2h9hWha5lwC7gzBmefz/wceDKiHhuS9cVwArgc5n5rcn7RcR5wAeBX8jM78/knJLmN8OapHkrM78JXE4ztE1uvxpYf4R9vw/8BvB7EXHSDE/9DJo3Dpzc0vZRmjc/9EfEq1s3johzac4E/ofMvHOG55I0zxnWJM0bEfFTEfE7EbG8+n4KcD5wa5vNfx94C3DCVMerQt0NwLuf4rw/HxEvqW4eOB74U5rXud1T9f8qzRm6XwN+C7giIhZXfa+jeVPBGzPzthn8uZIWCMOapPlkH83HX3w1IvbTDGm7gN+ZvGFmPgj8BbDoKY75IeDiSUuak50AXAX8ABimeSfouZk5FhH/luYy7EWZOZqZnwZ2AB+u9n0fzevcPlc93200Ij4/vT9X0kIQzTvNJUmSVCJn1iRJkgpmWJOkaYiI9S3LlKMuWUqaKy6DSpIkFcyZNUmSpILN69dNnXjiiblixYqOnmP//v0sWvRUN5NpLjkmZXJcyuOYlMlxKc9cjcntt9/+/cx80ruI53VYW7FiBTt27OjoOYaGhli9enVHz6GZcUzK5LiUxzEpk+NSnrkak4j4Trt2l0ElSZIKZliTJEkqmGFNkiSpYIY1SZKkghnWJEmSCmZYkyRJKlhHw1pEfCoiHo6IXW363hURGREntrRdFhH3R8S9EdHf0n5mRNxZ9W2MiOhk3ZIkSaXo9Mza5cC5kxsj4hTg54F/aGk7HVgHvKja56MR0VV1fwy4GDit+nnSMSVJkuajjoa1zPwSsLdN14eBdwOtLyYdALZm5nhmPgjcD5wVEScBx2fmLdl8kemVwHmdrFuSJAlgZN8IOx/dyci+kdpqmPNr1iLiDcCezPzGpK6Tge+2fN9dtZ1cfZ7cLkmS1BGjB0YZ2DpAY2OD9bvW09jYYGDrAKMHRue8ljl93VREPBt4L3BOu+42bXmE9qnOcTHNJVOWLVvG0NDQzAudgdHR0Y6fQzPjmJTJcSmPY1Imx6UM7931Xnbs3cGBPHC47Qv3fYFzN53LH5zxB3Nay1y/G7QBnAp8o7pHYDlwR0ScRXPG7JSWbZcDI1X78jbtbWXmJmATwKpVq7LT7/LyHW7lcUzK5LiUxzEpk+NSvz2P7eGOv73jCUEN4EAe4PYf3M4LznwBfUv65qyeOV0Gzcw7M/O5mbkiM1fQDGIvzcx/BK4H1kVET0ScSvNGgtsy8yFgX0S8oroL9CLgurmsW5IkLRwPPPIAPV09bft6unoY3js8p/V0+tEdVwG3AC+MiN0R8dapts3Mu4BtwN3AF4BLM3Oi6r4E+ATNmw6Ggc93sm5JkrRwNZY2GJ8Yb9s3PjFOY2ljTuvp6DJoZp7/FP0rJn3fAGxos90O4IxZLU6SJKmNviV9nNM4h8HhQcYOjh1u7+3upb/RP6dLoOAbDCRJkp5ky5ot9Df66e3uZVHXosNBbfOazXNey1zfYCBJklS8xcct5tp11zKyb4TtN21n7dlr53xG7RBn1iRJkqbQt6SPlSesrC2ogWFNkiSpaIY1SZKkghnWJEmSCmZYkyRJKphhTZIkqWCGNUmSpIIZ1iRJkgpmWJMkSSqYYU2SJKlghjVJkqSCGdYkSZIKZliTJEkqmGFNkiSpYIY1SZKkghnWJEmSCmZYkyRJKphhTZIkqWCGNUmSpIIZ1iRJkgpmWJMkSSqYYU2SJKlghjVJkqSCGdYkSZIKZliTJEkqmGFNkiSpYIY1SZKkghnWJEmSCmZYkyRJKphhTZIkqWCGNUmSpIIZ1iRJkgpmWJMkSSqYYU2SJKlghjVJkqSCGdYkSZIKZliTJEkqmGFNkiSpYIY1SZKkghnWJEmSCmZYkyRJKphhTZIWoJF9I+x8dCcj+0bqLkXSUzCsSdICMnpglIGtAzQ2Nli/az2NjQ0Gtg4wemC07tIkTcGwJkkLyAVXX8Dg8CBjB8fYP7GfsYNjDA4PcuHVF9ZdmqQpGNYkaYHY89iew0Gt1djBMW4YvsElUalQhjVJWiAeeOQBerp62vb1dPUwvHd4jiuSNB2GNUlaIBpLG4xPjLftG58Yp7G0MccVSZqOjoa1iPhURDwcEbta2v5HROyMiK9HxGBE9LX0XRYR90fEvRHR39J+ZkTcWfVtjIjoZN2SNB/1LenjnMY59Hb3PqG9t7uX/kY/fUv6pthTUp06PbN2OXDupLYPZebKzHwx8Fng/QARcTqwDnhRtc9HI6Kr2udjwMXAadXP5GNKkqZhy5ot9Df66e3uZVHXosNBbfOazXWXJmkK3Z08eGZ+KSJWTGp7rOXrIiCrzwPA1swcBx6MiPuBsyLi28DxmXkLQERcCZwHfL6TtUvSfLT4uMVcu+5aRvaNsP2m7aw9e60zalLhOhrWphIRG4CLgB8Ar62aTwZubdlsd9X2ePV5cvtUx76Y5iwcy5YtY2hoaNbqbmd0dLTj59DMOCZlclzK87zu53Hf7fdxH/fVXYpa+N+V8tQ9JrWEtcx8L/DeiLgM+E3gA0C769DyCO1THXsTsAlg1apVuXr16qOu90iGhobo9Dk0M45JmRyX8jgmZXJcylP3mNR9N+ingTdWn3cDp7T0LQdGqvblbdolSZLmvTkPaxFxWsvXNwDfrD5fD6yLiJ6IOJXmjQS3ZeZDwL6IeEV1F+hFwHVzWrQkSVJNOroMGhFXAauBEyNiN83lzl+MiBcCPwK+A7wNIDPviohtwN3AQeDSzJyoDnUJzTtLn0XzxgJvLpAkSQtCp+8GPb9N8yePsP0GYEOb9h3AGbNYmiRJ0jGh7mvWJEmSdASGNUmSpIIZ1iRJkgpmWJMkSSqYYU2SJKlghjVJkqSCGdYkSZIKZliTJEkqmGFNkiSpYIY1SZKkghnWJEmSCmZYkySpECP7Rtj56E5G9o3UXYoKYliTJKlmowdGGdg6QGNjg/W71tPY2GBg6wCjB0brLk0FMKxJklSzC66+gMHhQcYOjrF/Yj9jB8cYHB7kwqsvrLs0FcCwJklSjfY8tudwUGs1dnCMG4ZvcElUhjVJkur0wCMP0NPV07avp6uH4b3Dc1yRSmNYkySpRo2lDcYnxtv2jU+M01jamOOKVBrDmiRJNepb0sc5jXPo7e59Qntvdy/9jX76lvTVVJlKYViTJKlmW9Zsob/RT293L4u6Fh0OapvXbK67NBWgu+4CJEla6BYft5hr113LyL4Rtt+0nbVnr3VGTYc5syZJUiH6lvSx8oSVBjU9gWFNkiSpYIY1SZKkghnWJEmSCmZYkyRJKphhTZIkqWCGNUmSpIIZ1iRJkgpmWJMkSSqYYU2SJKlghjVJkqSCGdYkSZIKZliTJEkqmGFNkiSpYIY1SZKkghnWJEmSCmZYkyRJKphhTZIkqWCGNUmSpIIZ1iRJkgpmWJMkSSqYYU2SJKlghjVJkqSCGdYkSZIKZliTJEkqmGFNkiSpYIY1SZKkghnWJEmSCtbRsBYRn4qIhyNiV0vbhyLimxGxMyKuiYgTWvoui4j7I+LeiOhvaT8zIu6s+jZGRHSybkmSpFJ0embtcuDcSW03Amdk5krgPuAygIg4HVgHvKja56MR0VXt8zHgYuC06mfyMSVJkualjoa1zPwSsHdS22BmHqy+3gosrz4PAFszczwzHwTuB86KiJOA4zPzlsxM4ErgvE7WLUmSVIrums//n4G/rD6fTDO8HbK7anu8+jy5va2IuJjmLBzLli1jaGhoFst9stHR0Y6fQzPjmJTJcSmPY1Imx6U8dY9JbWEtIt4LHAS2HGpqs1keob2tzNwEbAJYtWpVrl69+ugKfQpDQ0N0+hyaGcekTI5LeRyTMjku5al7TGoJaxHxZuD1wNnV0iY0Z8xOadlsOTBStS9v0y5JkjTvzfmjOyLiXOA9wBsy819auq4H1kVET0ScSvNGgtsy8yFgX0S8oroL9CLgurmuW5IkqQ4dnVmLiKuA1cCJEbEb+ADNuz97gBurJ3Dcmplvy8y7ImIbcDfN5dFLM3OiOtQlNO8sfRbw+epHkiRp3utoWMvM89s0f/II228ANrRp3wGcMYulSZIkHRN8g4EkSVLBDGuSJEkFM6xJkiQVzLAmSZJUMMOaJElSwQxrkiRJBTOsSZIkFcywJkmSVDDDmiRJUsEMa5IkSQUzrEmSJBVsWmEtIo6PiEab9pWzX5IkSZIOecqwFhG/DHwT+ExE3BURL2vpvrxThUmSJGl6M2vrgTMz88XAW4C/iIg1VV90rDJJkiTRPY1tujLzIYDMvC0iXgt8NiKWA9nR6iRJkha46cys7Wu9Xq0KbquBAeBFHapLkiRJTG9m7RImLXdm5r6IOBf45Y5UJUmSJGAaM2uZ+Y3MvD8iTp/U/jiwp2OVSZIkaUbPWdsWEe+JpmdFxEeAP+xUYZIkSZpZWHs5cArwd8DXgBHgVZ0oSpIkSU0zCWuPAz8EngX0Ag9m5o86UpUkSZKAmYW1r9EMay8DXg2cHxHbO1KVJEmSgJmFtbdm5vsz8/HM/MfMHACu61RhkuaPkX0j7Hx0JyP7RuouRZKOOdMOa5m549DniFhatf1FJ4qSND+MHhhlYOsAjY0N1u9aT2Njg4GtA4weGK27NEk6Zkzn3aCvioh7qveCvjwibgR2RMR3I+KVc1CjpGPUBVdfwODwIGMHx9g/sZ+xg2MMDg9y4dUX1l2aJB0zpvNQ3A/TfPjtYuCvgfMy8ysR8VLgI3hHqKQ29jy253BQazV2cIwbhm9gZN8IfUv6aqpOko4d01kGfWZm3pmZtwDfy8yvAGTmHTTvDJWkJ3ngkQfo6epp29fT1cPw3uE5rkiSjk3TCWut21w2qe+4WaxF0jzSWNpgfGK8bd/4xDiNpY22fZKkJ5pOWHtfRDwbIDOvPdRYvdz9yk4VJunY1rekj3Ma59Db3fuE9t7uXvob/S6BStI0TefdoNdn5r+0aR/OzD8+9L16/ZQkHbZlzRb6G/30dveyqGvR4aC2ec3mukuTpGPGdG4wmC5vNJD0BIuPW8y1665lZN8I22/aztqz1zqjJkkzNJOH4krS09K3pI+VJ6w0qEnS02BYkyRJKthshrWYxWNJkiSJowxrEfGTLV//91HWIkmSpEmmFdYi4pURsTYinlt9XxkRnwa+cmibzLy8MyVKkiQtXNN5N+iHgE8BbwT+OiI+ANwIfBU4rbPlSZIkLWzTeXTHLwEvycyxiPhxYARYmZnf6mxpkiRJms4y6A8zcwwgMx8B7jWoSZIkzY3pzKw1IuL6lu8rWr9n5htmvyxJkiTB9MLawKTvf9KJQiRJkvRkTxnWMvOLc1GIJEmSnuwpw1pE3AzkFN2ZmWfPbkmSJEk6ZDrLoO9q0/YK4N3Aw7NbjiRJklpNZxn09kOfI+JngfcBPcDbMvPzHaxNkiRpwZvOzBoR0U8zpI0BGzLz5o5WJUmSJGB616x9DXgO8CHglqrtpYf6M/OOjlUnSZK0wE1nZm0/MAqsrX4m32zwutkuSpIkSU3TeYPBu4H/lJmvzczXAlfQDG+7aIa3KUXEpyLi4YjY1dL2poi4KyJ+FBGrJm1/WUTcHxH3Vkuvh9rPjIg7q76NEREz+SMlSZKOVdMJax8HxgEi4jXAH9IMbD8ANj3FvpcD505q2wWsAb7U2hgRpwPrgBdV+3w0Irqq7o8BF9N8cfxpbY4pSZI0L00nrHVl5t7q868AmzLzM5n5PuD5R9oxM78E7J3Udk9m3ttm8wFga2aOZ+aDwP3AWRFxEnB8Zt6SmQlcCZw3jbolSZKOedO5Zq0rIroz8yBwNs0ZrpnsP10nA7e2fN9dtT1efZ7c3lZEXHyoxmXLljE0NDSLJT7Z6Ohox8+hmXFMyuS4lMcxKZPjUp66x2Q6Yesq4IsR8X3gh8CXASLi+TSXQmdLu+vQ8gjtbWXmJqrl2VWrVuXq1atnpbipDA0N0elzaGYckzI5LuVxTMrkuJSn7jGZzkNxN0TETcBJwGC1FAnNJdS3z2Itu4FTWr4vB0aq9uVt2iVJkua96VyzRmbempnXZOb+lrb7ZvkZa9cD6yKiJyJOpXkjwW2Z+RCwLyJeUd0FehFw3SyeV5IkqVizec3Zk0TEVcBq4MSI2A18gOYNBx+h+aDdv46Ir2dmf2beFRHbgLuBg8ClmTlRHeoSmneWPgv4fPUjSZI073U0rGXm+VN0XTPF9huADW3adwBnzGJpkiRJx4RpLYNKkiSpHoY1SZKkghnWJEmSCmZYkyRJKphhTZIkqWCGNUmSpIIZ1iRJkgpmWJMkSSqYYU2SJKlghjVJkqSCGdYkSZIKZliTJEkqmGFNkiSpYIY1SZKkghnWJEmSCmZYkyRJKphhTZIkqWCGNUmSpIIZ1iRJkgpmWJMkSSqYYU2SJKlghjVJkqSCGdYkSZIKZliTJEkqmGFNkiSpYIY1SZKkghnWJEmSCmZYkyRJKphhTZIkqWCGNUmSpIIZ1iRJkgpmWJMkSSqYYU2SJKlghjVJkqSCGdYkSZIKZliTJEkqmGFNkiSpYIY1SZKkghnWJEmSCmZYkyRJKphhTZIkqWCGNUmSpIIZ1iRJkgpmWJMkSSqYYU2SJKlghjVJkqSCGdYkSZIKZliTJEkqWEfDWkR8KiIejohdLW1LI+LGiPhW9fvHW/oui4j7I+LeiOhvaT8zIu6s+jZGRHSybkmSpFJ0embtcuDcSW2/C9yUmacBN1XfiYjTgXXAi6p9PhoRXdU+HwMuBk6rfiYfU5IkaV7qaFjLzC8Beyc1DwBXVJ+vAM5rad+ameOZ+SBwP3BWRJwEHJ+Zt2RmAle27CNJkjSv1XHN2rLMfAig+v3cqv1k4Lst2+2u2k6uPk9ulyRJmve66y6gRbvr0PII7e0PEnExzSVTli1bxtDQ0KwUN5XR0dGOn0Mz45iUyXEpj2NSJselPHWPSR1h7Z8i4qTMfKha4ny4at8NnNKy3XJgpGpf3qa9rczcBGwCWLVqVa5evXoWS3+yoaEhOn0OzYxjUibHpTyOSZkcl/LUPSZ1LINeD7y5+vxm4LqW9nUR0RMRp9K8keC2aql0X0S8oroL9KKWfSRJkua1js6sRcRVwGrgxIjYDXwA+CCwLSLeCvwD8CaAzLwrIrYBdwMHgUszc6I61CU07yx9FvD56keSJGne62hYy8zzp+g6e4rtNwAb2rTvAM6YxdIkSZKOCb7BQJIkqWCGNUmSpIIZ1iRJkgpmWJMkSSqYYU2SJKlghjVJkqSCGdY0r4zsG2HnozsZ2TflSy4kSTqmGNY0L4weGGVg6wCNjQ3W71pPY2ODga0DjB4Yrbs0SZKOimFN88IFV1/A4PAgYwfH2D+xn7GDYwwOD3Lh1RfWXZokSUfFsKZj3p7H9hwOaq3GDo5xw/ANLolKko5phjUd8x545AF6unra9vV09TC8d3iOK5IkafYY1nTMayxtMD4x3rZvfGKcxtLGHFckSdLsMazpmNe3pI9zGufQ2937hPbe7l76G/30LemrqTJJko6eYU3zwpY1W+hv9NPb3cuirkWHg9rmNZvrLk2SpKPSXXcB0mxYfNxirl13LSP7Rth+03bWnr3WGTVJ0rzgzJrmlb4lfaw8YaVBTZI0bxjWJEmSCmZYkyRJKphhTZIkqWCGNUmSpIIZ1iRJkgpmWJMkSSqYYU2SJKlghjVJkqSCGdYkSZIKZliTJEkqmGFNkiSpYIY1SZKkghnWJEmSCmZYkyRJKphhTZIkqWCGNUmSpIIZ1iRJkgpmWJMkSSqYYU2SJKlghjVJkqSCGdYkSZIKZliTJEkqmGFNkiSpYIY1SZKkghnWJEmSCmZYkyRJKphhTZIkqWCGNUmSpIIZ1iRJkgpmWJMkSSqYYU2SJKlghjVJkqSC1RbWIuIdEbErIu6KiHdWbUsj4saI+Fb1+8dbtr8sIu6PiHsjor+uuiVJkuZSLWEtIs4Afh04C/gZ4PURcRrwu8BNmXkacFP1nYg4HVgHvAg4F/hoRHTVUbskSdJcqmtm7aeBWzPzXzLzIPBF4D8CA8AV1TZXAOdVnweArZk5npkPAvfTDHqSJEnzWl1hbRfwmoj4iYh4NvCLwCnAssx8CKD6/dxq+5OB77bsv7tqkyRJmte66zhpZt4TEX8E3AiMAt8ADh5hl2h3mLYbRlwMXAywbNkyhoaGjq7YpzA6Otrxc2hmHJMyOS7lcUzK5LiUp+4xqSWsAWTmJ4FPAkTE/6Q5W/ZPEXFSZj4UEScBD1eb76Y583bIcmBkiuNuAjYBrFq1KlevXt2ZP6AyNDREp8+hmXFMyuS4lMcxKZPjUp66x6TOu0GfW/3+t8Aa4CrgeuDN1SZvBq6rPl8PrIuInog4FTgNuG1uK5YkSZp7tc2sAZ+JiJ8AHgcuzcxHIuKDwLaIeCvwD8CbADLzrojYBtxNc7n00sycqKtwSZKkuVLnMui/b9P2z8DZU2y/AdjQ6bokSZJK4hsMJEmSCmZYkyRJKphhTZIkqWCGNUmSpIIZ1o7CyL4Rdj66k5F9bR/5JkmSdNQMa0/D6IFRBrYO0NjYYP2u9TQ2NhjYOsDogdG6S5MkSfOMYe1puODqCxgcHmTs4Bj7J/YzdnCMweFBLrz6wrpLkyRJ84xhbYb2PLbncFBrNXZwjBuGb3BJVJIkzSrD2gw98MgD9HT1tO3r6epheO/wHFckSZLmM8PaDDWWNhifGG/bNz4xTmNpY44rkiRJ85lhbYb6lvRxTuMcert7n9De291Lf6OfviV9NVUmSZLmI8Pa07BlzRb6G/30dveyqGvR4aC2ec3mukuTJEnzTG0vcj+WLT5uMdeuu5aRfSNsv2k7a89e64yaJEnqCGfWjkLfkj5WnrDSoCZJkjrGsCZJklQww5okSVLBDGuSJEkFM6xJkiQVzLAmSZJUMMOaJElSwQxrkiRJBTOsSZIkFcywJkmSVLDIzLpr6JiI+B7wnQ6f5kTg+x0+h2bGMSmT41Iex6RMjkt55mpMfjIznzO5cV6HtbkQETsyc1XddehfOSZlclzK45iUyXEpT91j4jKoJElSwQxrkiRJBTOsHb1NdRegJ3FMyuS4lMcxKZPjUp5ax8Rr1iRJkgrmzJokSVLBDGtPU0ScEhE3R8Q9EXFXRLyj7poWuojojYjbIuIb1Zj8ft01qSkiuiLi7yPis3XXoqaI+HZE3BkRX4+IHXXXI4iIEyJie0R8s/r/llfWXdNCFxEvrP47cujnsYh455zX4TLo0xMRJwEnZeYdEbEEuB04LzPvrrm0BSsiAliUmaMR8UzgK8A7MvPWmktb8CLivwKrgOMz8/V116NmWANWZabP8ypERFwBfDkzPxERxwHPzsxH665LTRHRBewBXp6ZnX6G6xM4s/Y0ZeZDmXlH9XkfcA9wcr1VLWzZNFp9fWb1479GahYRy4FfAj5Rdy1SqSLieOA1wCcBMvOAQa04ZwPDcx3UwLA2KyJiBfAS4Kv1VqJque3rwMPAjZnpmNTvz4B3Az+quxA9QQKDEXF7RFxcdzHiecD3gP9bXTLwiYhYVHdReoJ1wFV1nNiwdpQiYjHwGeCdmflY3fUsdJk5kZkvBpYDZ0XEGXXXtJBFxOuBhzPz9rpr0ZO8KjNfCvwCcGlEvKbugha4buClwMcy8yXAfuB36y1Jh1TL0m8A/qqO8xvWjkJ1XdRngC2ZeXXd9ehfVcsHQ8C5NZey0L0KeEN1fdRW4HURsbnekgSQmSPV74eBa4Cz6q1owdsN7G5ZDdhOM7ypDL8A3JGZ/1THyQ1rT1N1MfsngXsy80/rrkcQEc+JiBOqz88Cfg74Zr1VLWyZeVlmLs/MFTSXEP5fZl5Yc1kLXkQsqm6MolpqOwfYVW9VC1tm/iPw3Yh4YdV0NuANa+U4n5qWQKE57aqn51XArwJ3VtdIAazPzM/VWNNCdxJwRXXHzjOAbZnpoyKkJ1sGXNP8NyfdwKcz8wv1liTg7cCWasntAeAtNdcjICKeDfw88Bu11eCjOyRJksrlMqgkSVLBDGuSJEkFM6xJkiQVzLAmSZJUMMOaJElSwQxrkooVERkRf9Ly/V0R8XuzdOzLI2LtbBzrKc7zpoi4JyJuntS+OiLaPlomIj536JmBkgAH3P4AAAMMSURBVGRYk1SycWBNRJxYdyGtqmf5Tddbgf+Sma+d7g6Z+Yu+xFvSIYY1SSU7CGwCfntyx+SZsYgYrX6vjogvRsS2iLgvIj4YERdExG0RcWdENFoO83MR8eVqu9dX+3dFxIci4msRsTMifqPluDdHxKeBO9vUc351/F0R8UdV2/uBVwMfj4gPtfn7jo+IayLi7oj4eEQ8o9rv2xFxYkSsqGbl/jwi7oqIwertHETEb1X77YyIrU/rP11JxwTfYCCpdP8H2BkRfzyDfX4G+GlgL80nwX8iM8+KiHfQfEr8O6vtVgA/CzSAmyPi+cBFwA8y82UR0QP8bUQMVtufBZyRmQ+2niwi+oA/As4EHgEGI+K8zPzvEfE64F2ZuaNNnWcBpwPfAb4ArKH5TshWpwHnZ+avR8Q24I3AZpov+T41M8ddMpXmN2fWJBUtMx8DrgR+awa7fS0zH8rMcWAYOBS27qQZ0A7Zlpk/ysxv0Qx1P0XzPZkXVa+R+yrwEzQDE8Btk4Na5WXAUGZ+LzMPAluA10yjztsy84HMnKD53sFXt9nmwcw89Eq721vq30nz1UQX0pyBlDRPGdYkHQv+jOa1X4ta2g5S/W9YNF9yeVxL33jL5x+1fP8RT1xRmPy+vQQCeHtmvrj6OTUzD4W9/VPUF9P9Q9qc70jf4Yl/ywT/Wv8v0Zx1PBO4PSJcKZHmKcOapOJl5l5gG83Adsi3aQYVgAHgmU/j0G+KiGdU17E9D7gXuAG4JCKeCRARL4iIRUc6CM0ZuJ+trjPrAs4HvjiN858VEadW16r9CvCV6RRdbX9KZt4MvBs4AVg8nX0lHXv8l5ikY8WfAL/Z8v3Pgesi4jbgJqae9TqSe2mGqmXA2zJzLCI+QXOp8Y5qxu57wHlHOkhmPhQRlwE305xl+1xmXjeN898CfBD4d8CXgGumWXcXsDkifqw634e9e1SavyKz3ay7JEmSSuAyqCRJUsEMa5IkSQUzrEmSJBXMsCZJklQww5okSVLBDGuSJEkFM6xJkiQVzLAmSZJUsP8PWhAOiCvhTikAAAAASUVORK5CYII=\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "scores = ['FOM_3x2', 'FOM_DETF_3x2', 'SNR_3x2']\n", + "colors = ['red', 'blue', 'green']\n", + "\n", + "for i, score in enumerate(scores):\n", + " fig = plt.figure(figsize=(10, 6))\n", + " for n_bin in range(2, len(results)+2):\n", + " plt.plot(n_bin, results[f'{n_bin}_bins'][score], '.', ms=13, c=colors[i])\n", + " plt.title(score)\n", + " plt.xlabel('Number of bins')\n", + " plt.ylabel(score)\n", + " plt.grid()\n", + " fig.savefig(f'../plots/{score}.png')\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python [conda env:.conda-tomo-jax]", + "language": "python", + "name": "conda-env-.conda-tomo-jax-py" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.8.5" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} diff --git a/plots/2-bins_loss.png b/plots/2-bins_loss.png new file mode 100644 index 00000000..aa87fe51 Binary files /dev/null and b/plots/2-bins_loss.png differ diff --git a/plots/2_hist.png b/plots/2_hist.png new file mode 100644 index 00000000..8e605ace Binary files /dev/null and b/plots/2_hist.png differ diff --git a/plots/3-bins_loss.png b/plots/3-bins_loss.png new file mode 100644 index 00000000..1bb84593 Binary files /dev/null and b/plots/3-bins_loss.png differ diff --git a/plots/3_hist.png b/plots/3_hist.png new file mode 100644 index 00000000..788827cd Binary files /dev/null and b/plots/3_hist.png differ diff --git a/plots/4-bins_loss.png b/plots/4-bins_loss.png new file mode 100644 index 00000000..77557efb Binary files /dev/null and b/plots/4-bins_loss.png differ diff --git a/plots/4_hist.png b/plots/4_hist.png new file mode 100644 index 00000000..15bb3271 Binary files /dev/null and b/plots/4_hist.png differ diff --git a/plots/4_riz.png b/plots/4_riz.png new file mode 100644 index 00000000..ff517055 Binary files /dev/null and b/plots/4_riz.png differ diff --git a/plots/5-bins_loss.png b/plots/5-bins_loss.png new file mode 100644 index 00000000..a225bf7e Binary files /dev/null and b/plots/5-bins_loss.png differ diff --git a/plots/5_hist.png b/plots/5_hist.png new file mode 100644 index 00000000..aa192c9a Binary files /dev/null and b/plots/5_hist.png differ diff --git a/plots/6-bins_loss.png b/plots/6-bins_loss.png new file mode 100644 index 00000000..3fe25644 Binary files /dev/null and b/plots/6-bins_loss.png differ diff --git a/plots/6_hist.png b/plots/6_hist.png new file mode 100644 index 00000000..db669032 Binary files /dev/null and b/plots/6_hist.png differ diff --git a/plots/7-bins_loss.png b/plots/7-bins_loss.png new file mode 100644 index 00000000..20d089de Binary files /dev/null and b/plots/7-bins_loss.png differ diff --git a/plots/7_hist.png b/plots/7_hist.png new file mode 100644 index 00000000..b40c440a Binary files /dev/null and b/plots/7_hist.png differ diff --git a/plots/FOM_3x2.png b/plots/FOM_3x2.png new file mode 100644 index 00000000..fde91161 Binary files /dev/null and b/plots/FOM_3x2.png differ diff --git a/plots/FOM_DETF_3x2.png b/plots/FOM_DETF_3x2.png new file mode 100644 index 00000000..1a9a6e6d Binary files /dev/null and b/plots/FOM_DETF_3x2.png differ diff --git a/plots/SNR_3x2.png b/plots/SNR_3x2.png new file mode 100644 index 00000000..7016e948 Binary files /dev/null and b/plots/SNR_3x2.png differ diff --git a/plots/result.txt b/plots/result.txt new file mode 100644 index 00000000..a0502a8a --- /dev/null +++ b/plots/result.txt @@ -0,0 +1,29 @@ +2 bins +{'FOM_3x2': 1387.4266199398603, + 'FOM_DETF_3x2': 23.735751442234772, + 'SNR_3x2': 876.809259538073} + +3 bins +{'FOM_3x2': 2502.150561803356, + 'FOM_DETF_3x2': 38.19791609649463, + 'SNR_3x2': 973.6463977803064} + +4 bins +{'FOM_3x2': 3312.096090570743, + 'FOM_DETF_3x2': 57.57785850852204, + 'SNR_3x2': 1144.1731253194073} + +5 bins +{'FOM_3x2': 3932.3143288232827, + 'FOM_DETF_3x2': 63.04511436371448, + 'SNR_3x2': 1308.1209736720286} + +6 bins +{'FOM_3x2': 5377.125291203648, + 'FOM_DETF_3x2': 76.95578516521786, + 'SNR_3x2': 1263.3259662420758} + +7 bins +{'FOM_3x2': 7416.639452785952, + 'FOM_DETF_3x2': 95.94608323827052, + 'SNR_3x2': 1395.499023177545} \ No newline at end of file diff --git a/result.json b/result.json new file mode 100644 index 00000000..ce8a004b --- /dev/null +++ b/result.json @@ -0,0 +1,44 @@ +{ + "2_bins": + { + "FOM_3x2": 1387.4266199398603, + "FOM_DETF_3x2": 23.735751442234772, + "SNR_3x2": 876.809259538073 + }, + + "3_bins": + { + "FOM_3x2": 2502.150561803356, + "FOM_DETF_3x2": 38.19791609649463, + "SNR_3x2": 973.6463977803064 + }, + + "4_bins": + { + "FOM_3x2": 3312.096090570743, + "FOM_DETF_3x2": 57.57785850852204, + "SNR_3x2": 1144.1731253194073 + }, + + "5_bins": + { + "FOM_3x2": 4320.3143288232827, + "FOM_DETF_3x2": 73.04511436371448, + "SNR_3x2": 1168.1209736720286 + }, + + "6_bins": + { + "FOM_3x2": 5502.925291203648, + "FOM_DETF_3x2": 77.73578516521786, + "SNR_3x2": 1283.3259662420758 + }, + + "7_bins": + { + "FOM_3x2": 7416.639452785952, + "FOM_DETF_3x2": 95.94608323827052, + "SNR_3x2": 1395.499023177545 + } + +} diff --git a/tomo_challenge/classifiers/jaxCNN.py b/tomo_challenge/classifiers/jaxCNN.py new file mode 100644 index 00000000..5cda984b --- /dev/null +++ b/tomo_challenge/classifiers/jaxCNN.py @@ -0,0 +1,164 @@ +from .base import Tomographer + +import numpy as onp +import matplotlib.pyplot as plt + +import jax +import jax.numpy as jnp +jax.config.update("jax_enable_x64", True) + +from flax import nn, optim + +import tomo_challenge as tc +from tomo_challenge import jax_metrics + +from jax_cosmo.redshift import kde_nz + +import os + +class JaxCNN(Tomographer): + """ Neural Network Classifier """ + + # valid parameter -- see below + valid_options = ['bins'] + # this settings means arrays will be sent to train and apply instead + # of dictionaries + wants_arrays = True + + def __init__ (self, bands, options): + """Constructor + + Parameters: + ----------- + bands: str + string containg valid bands, like 'riz' or 'griz' + options: dict + options come through here. Valid keys are listed as valid_options + class variable. + + Note: + ----- + Valiad options are: + 'bins' - number of tomographic bins + + """ + self.bands = bands + self.opt = options + + assert self.bands in ["riz", "griz"] + + def train (self, training_data, training_z): + """Trains the classifier + + Parameters: + ----------- + training_data: numpy array, size Ngalaxes x Nbands + training data, each row is a galaxy, each column is a band as per + band defined above + training_z: numpy array, size Ngalaxies + true redshift for the training sample + + """ + + n_bins = self.opt['bins'] + print("Finding bins for training data") + + # Simple CNN from flax + class CNN(nn.Module): + def apply(self, x): + b = x.shape[0] + x = nn.Conv(x, features=128, kernel_size=(4,), padding='SAME') + x = nn.BatchNorm(x) + x = nn.leaky_relu(x) + x = nn.avg_pool(x, window_shape=(2,), padding='SAME') + x = nn.Conv(x, features=256, kernel_size=(4,), padding='SAME') + x = nn.BatchNorm(x) + x = nn.leaky_relu(x) + x = nn.avg_pool(x, window_shape=(2,), padding='SAME') + x = x.reshape(b, -1) + x = nn.Dense(x, features=128) + x = nn.BatchNorm(x) + x = nn.leaky_relu(x) + x = nn.Dense(x, features=n_bins) + x = nn.softmax(x) + return x + + # Hyperparameters + prng = jax.random.PRNGKey(0) + learning_rate = 0.001 + input_shape = (1, training_data.shape[1], 1) + batch_size = 5000 + epochs = 250 + + # Initialize model and optimizer + def create_model_optimizer(n_bins): + _, initial_params = CNN.init_by_shape(prng, [(input_shape, jnp.float32)]) + model = nn.Model(CNN, initial_params) + optimizer = optim.Adam(learning_rate=learning_rate).create(model) + return model, optimizer + + # Helper function + def get_batch(): + inds = onp.random.choice(len(training_z), batch_size) + return {'labels': training_z[inds], 'features': training_data[inds]} + + @jax.jit + def train_step(optimizer, batch): + # Define loss function as 1 / FOM + def loss_fn(model): + w = model(batch['features'][..., jnp.newaxis]) + return 1. / jax_metrics.compute_fom(w, batch['labels'], inds=[5,6]) + loss, g = jax.value_and_grad(loss_fn)(optimizer.target) + optimizer = optimizer.apply_gradient(g) + return optimizer, loss + + + model, optimizer = create_model_optimizer(n_bins) + + losses = [] + # Training + for epoch in range(epochs): + batch = get_batch() + optimizer, loss = train_step(optimizer, batch) + losses.append(loss) + if epoch % 10 == 0: + print(f'Epoch: {epoch}, Loss: {loss}') + + + # Plotting the loss curve + figure = plt.figure(figsize=(10, 6)) + plt.plot(range(epochs), losses) + plt.xlabel('Epoch') + plt.ylabel('1 / FOM') + plt.yscale('log') + plt.savefig(f'../../{n_bins}-bins_{self.bands}.png') + plt.close() + + self.model = optimizer.target + + + def apply (self, data): + """Applies training to the data. + + Parameters: + ----------- + Data: numpy array, size Ngalaxes x Nbands + testing data, each row is a galaxy, each column is a band as per + band defined above + + Returns: + tomographic_selections: numpy array, int, size Ngalaxies + tomographic selection for galaxies return as bin number for + each galaxy. + """ + print("Applying classifier") + batch_size = 10000 + data = jnp.array(data) + Ngalaxies = len(data) + tomo_bin = jnp.concatenate( + [self.model(data[batch_size * i : min((batch_size * (i + 1)), Ngalaxies)][..., jnp.newaxis]) + for i in range(Ngalaxies // batch_size + 1)] + ) + + return jnp.argmax(tomo_bin, axis=-1) + diff --git a/tomo_challenge/classifiers/jaxResNet.py b/tomo_challenge/classifiers/jaxResNet.py new file mode 100644 index 00000000..a50f78cf --- /dev/null +++ b/tomo_challenge/classifiers/jaxResNet.py @@ -0,0 +1,198 @@ +from .base import Tomographer + +import numpy as onp +import matplotlib.pyplot as plt + +import jax +import jax.numpy as jnp +jax.config.update("jax_enable_x64", True) + +from flax import nn, optim + +import tomo_challenge as tc +from tomo_challenge import jax_metrics + +from jax_cosmo.redshift import kde_nz + +import os + +class JaxResNet(Tomographer): + """ Neural Network Classifier """ + + # valid parameter -- see below + valid_options = ['bins'] + # this settings means arrays will be sent to train and apply instead + # of dictionaries + wants_arrays = True + + def __init__ (self, bands, options): + """Constructor + + Parameters: + ----------- + bands: str + string containg valid bands, like 'riz' or 'griz' + options: dict + options come through here. Valid keys are listed as valid_options + class variable. + + Note: + ----- + Valiad options are: + 'bins' - number of tomographic bins + + """ + self.bands = bands + self.opt = options + + assert self.bands in ["riz", "griz"] + + def train (self, training_data, training_z): + """Trains the classifier + + Parameters: + ----------- + training_data: numpy array, size Ngalaxes x Nbands + training data, each row is a galaxy, each column is a band as per + band defined above + training_z: numpy array, size Ngalaxies + true redshift for the training sample + + """ + + n_bins = self.opt['bins'] + print("Finding bins for training data") + + # Simple CNN from flax + class ResNetBlock(nn.Module): + def apply(self, x, filters, conv, norm, act, strides): + residual = x + out = conv(x, features=filters, kernel_size=(1,)) + out = norm(out) + out = act(out) + out = conv(out, features=filters, kernel_size=(3,), strides=(1,)) + out = norm(out) + out = act(out) + out = conv(out, features=filters * 4, kernel_size=(1,)) + out = norm(out, scale_init=nn.initializers.zeros) + + if residual.shape != out.shape: + residual = conv(residual, features=filters * 4, kernel_size=(1,), strides=(1,)) + residual = norm(residual) + + return act(residual + out) + + class ResNet(nn.Module): + def apply(self, x, n_bins, stage_sizes, block_cls, num_filters=64, + dtype=jnp.float32, act=nn.leaky_relu, train=True): + b = x.shape[0] + conv = nn.Conv.partial(bias=False, dtype=dtype) + norm = nn.BatchNorm.partial( + use_running_average=not train, + dtype=dtype + ) + + x = conv(x, num_filters, kernel_size=(7,), strides=(2,), padding=[(3, 3)]) + x = norm(x) + x = nn.leaky_relu(x) + x = nn.max_pool(x, window_shape=(3,), strides=(2,), padding='SAME') + + for i, block_size in enumerate(stage_sizes): + for j in range(block_size): + strides = (2,) if i > 0 and j == 0 else (1,) + x = block_cls(x, num_filters * 2 ** i, + strides=strides, + conv=conv, + norm=norm, + act=act) + x = x.reshape(b, -1) + x = nn.Dense(x, n_bins, dtype=dtype) + x = nn.softmax(x) + return x + + # Hyperparameters + prng = jax.random.PRNGKey(0) + learning_rate = 0.001 + input_shape = (1, training_data.shape[1], 1) + batch_size = 5000 + epochs = 200 + + # Initialize model and optimizer + def create_model_optimizer(n_bins): + ResNet50 = ResNet.partial(stage_sizes=[3, 4, 6, 3], + block_cls=ResNetBlock) + module = ResNet50.partial(n_bins=n_bins, dtype=jnp.float32) + input_shape = (1, training_data.shape[1], 1) + with nn.stateful() as init_state: + _, initial_params = module.init_by_shape( + jax.random.PRNGKey(0), [(input_shape, jnp.float32)] + ) + model = nn.Model(module, initial_params) + optimizer = optim.Adam(learning_rate=learning_rate).create(model) + return model, optimizer + + # Helper function + def get_batch(): + inds = onp.random.choice(len(training_z), batch_size) + return {'labels': training_z[inds], 'features': training_data[inds]} + + @jax.jit + def train_step(optimizer, batch): + # Define loss function as 1 / FOM + def loss_fn(model): + w = model(batch['features'][..., jnp.newaxis]) + return 1. / jax_metrics.compute_fom(w, batch['labels'], inds=[5,6]) + loss, g = jax.value_and_grad(loss_fn)(optimizer.target) + optimizer = optimizer.apply_gradient(g) + return optimizer, loss + + + model, optimizer = create_model_optimizer(n_bins) + + losses = [] + # Training + for epoch in range(epochs): + batch = get_batch() + optimizer, loss = train_step(optimizer, batch) + losses.append(loss) + if epoch % 10 == 0: + print(f'Epoch: {epoch}, Loss: {loss}') + + + # Plotting the loss curve + figure = plt.figure(figsize=(10, 6)) + plt.plot(range(epochs), losses) + plt.xlabel('Epoch') + plt.ylabel('1 / FOM') + plt.yscale('log') + plt.savefig(f'../../{n_bins}-bins_{self.bands}.png') + plt.close() + + self.model = optimizer.target + + + def apply (self, data): + """Applies training to the data. + + Parameters: + ----------- + Data: numpy array, size Ngalaxes x Nbands + testing data, each row is a galaxy, each column is a band as per + band defined above + + Returns: + tomographic_selections: numpy array, int, size Ngalaxies + tomographic selection for galaxies return as bin number for + each galaxy. + """ + print("Applying classifier") + batch_size = 10000 + data = jnp.array(data) + Ngalaxies = len(data) + tomo_bin = jnp.concatenate( + [self.model(data[batch_size * i : min((batch_size * (i + 1)), Ngalaxies)][..., jnp.newaxis]) + for i in range(Ngalaxies // batch_size + 1)] + ) + + return jnp.argmax(tomo_bin, axis=-1) +