diff --git a/__init__.py b/__init__.py index 54f4ba8..e18f7f4 100644 --- a/__init__.py +++ b/__init__.py @@ -1,4 +1,5 @@ from .mvcomp import mysort, feature_gen, norm_covar_inv, mah_dist_feat_mat, mah_dist_mat_2_roi, subject_list, feature_list, compute_average, model_comp, spatial_mvcomp, dist_plot, correlation_fig +from .mvcomp import model_comp_simplified, compute_average_simplified from .version import __version__ # initial hack to not import optional plotting functions if necessary packages do not exist diff --git a/examples/Example4_model_comp_simplified.ipynb b/examples/Example4_model_comp_simplified.ipynb new file mode 100644 index 0000000..378d27d --- /dev/null +++ b/examples/Example4_model_comp_simplified.ipynb @@ -0,0 +1,263 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [], + "source": [ + "from glob import glob\n", + "import sys\n", + "sys.path.append('../')\n", + "import mvcomp as mvc\n", + "import nibabel as nb\n", + "import numpy as np\n", + "from matplotlib import pyplot as plt\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Generate false data\n", + "- pull data for three contrasts from one individual from nitrc\n", + " - https://www.nitrc.org/frs/?group_id=1205#\n", + "- use this individual's three metrics as a toy example, modifying each metric by noise\n", + "\n", + "```text\n", + " sub001_sess1_INV2.nii.gz\n", + "22.95 MB\t660\tAny\t.gz\t10.1038/sdata.2014.54\n", + " \t\n", + " sub001_sess1_T1map.nii.gz\n", + "40.05 MB\t684\tAny\t.gz\t10.1038/sdata.2014.54\n", + " \t\n", + " sub001_sess1_T1w.nii.gz\n", + "43.06 MB\t677\tAny\t.gz\t10.1038/sdata.2014.54\n", + "\n", + "...\n", + "```" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "(240, 320, 320)\n", + "./sub002_sess1_INV2.nii.gz\n", + "(240, 320, 320)\n", + "./sub002_sess1_T1w.nii.gz\n", + "(240, 320, 320)\n", + "./sub002_sess1_T1map.nii.gz\n", + "(240, 320, 320)\n", + "./sub003_sess1_INV2.nii.gz\n", + "(240, 320, 320)\n", + "./sub003_sess1_T1w.nii.gz\n", + "(240, 320, 320)\n", + "./sub003_sess1_T1map.nii.gz\n", + "(240, 320, 320)\n", + "new_ID: sub004 should be an outlier\n", + "./sub004_sess1_INV2.nii.gz\n", + "(240, 320, 320)\n", + "new_ID: sub004 should be an outlier\n", + "./sub004_sess1_T1w.nii.gz\n", + "(240, 320, 320)\n", + "new_ID: sub004 should be an outlier\n", + "./sub004_sess1_T1map.nii.gz\n" + ] + } + ], + "source": [ + "# simulate data by simply multiplying the one subject by noise\n", + "# here we make the final subject MUCH more different (by scaling the noise by 10) to\n", + "# ensure that we can distinguish them properly by their D2\n", + "\n", + "ID='sub001'\n", + "s1 = glob('./sub001*')\n", + "idx=2\n", + "for rep in range(3):\n", + " idx2 = rep+idx\n", + " for s in s1:\n", + " new_ID = ID.replace(\"1\",str(idx2))\n", + " img = nb.load(s)\n", + " d=img.get_fdata()\n", + " print(img.shape)\n", + " d[np.isnan(d)]=0\n", + " d[np.isinf(d)]=0\n", + " if rep == 2:\n", + " print(f\"new_ID: {new_ID} should be an outlier\")\n", + " scale=10 #create an outlier\n", + " else:\n", + " scale=1 \n", + " d=d*scale*np.random.random(d.shape) #scale up the data\n", + " out_fname = s.replace(ID,new_ID)\n", + " print(out_fname)\n", + " out_img = nb.Nifti1Image(d,affine=img.affine,header=img.header)\n", + " out_img.update_header()\n", + " out_img.to_filename(out_fname)\n" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "[['./sub001_sess1_INV2.nii.gz', './sub001_sess1_T1w.nii.gz', './sub001_sess1_T1map.nii.gz'], ['./sub002_sess1_T1map.nii.gz', './sub002_sess1_T1w.nii.gz', './sub002_sess1_INV2.nii.gz'], ['./sub003_sess1_INV2.nii.gz', './sub003_sess1_T1map.nii.gz', './sub003_sess1_T1w.nii.gz'], ['./sub004_sess1_INV2.nii.gz', './sub004_sess1_T1map.nii.gz', './sub004_sess1_T1w.nii.gz']]\n", + "No 'model_feature_image_list' was provided\n", + "\t- Model features will be iteratively computed as the mean of all other subjects (leave one out)\n", + "subject None feature matrix creation in 5.66 s\n", + "subject None feature matrix creation in 5.9 s\n", + "subject None feature matrix creation in 5.9 s\n", + "subject None feature matrix creation in 6.03 s\n", + "Total time for mahalanobis distance calculation on 4 subjects with 4702024 voxels: 7.03s\n" + ] + } + ], + "source": [ + "s1 = glob('./sub001*')\n", + "s2 = glob('./sub002*')\n", + "s3 = glob('./sub003*')\n", + "s4 = glob('./sub004*')\n", + "all_Ss = [s1,s2,s3,s4]\n", + "print(all_Ss)\n", + "\n", + "img = nb.load(s1[0])\n", + "_mask_d = (img.get_fdata()>100).astype(int)\n", + "mask_img = nb.Nifti1Image(_mask_d,affine=img.affine,header=img.header)\n", + "res = mvc.model_comp_simplified(all_Ss,mask=mask_img,verbosity=2)" + ] + }, + { + "cell_type": "code", + "execution_count": 15, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 15, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmkAAAGbCAYAAACfwwddAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAdJUlEQVR4nO3df6zd5X0f8PcnmAAq4GACEeGSmQykBEpL4RaYEqFsWYCwKGQNqVx1w1VcWcqoRNVNG1n/oEnUKpmUJouyZGNNFMI6fixtBYoEqQuNKlUpxhRaAwmzG1i44AYHE4dIgdjOsz/O18nBXF9f+177PPf69ZKOzvd8zvd5zvc891i8eb6/qrUWAAD68ppJbwAAAK8mpAEAdEhIAwDokJAGANAhIQ0AoEMrJr0Bi+31r399W7169aQ3AwDggB566KHvtdZOm+29ZRfSVq9enU2bNk16MwAADqiq/t/+3rO7EwCgQ0IaAECHhDQAgA4tu2PSAICjy65duzIzM5OXXnpp0puyX8cff3ympqZy7LHHzruNkAYALGkzMzM56aSTsnr16lTVpDfnVVpref755zMzM5Ozzz573u3s7gQAlrSXXnopp556apcBLUmqKqeeeupBz/QJaQDAktdrQNvrULZPSAMA6JBj0gCAZWXzzM5F7e+CqZUHXOell17K5Zdfnpdffjm7d+/Otddem4985CML+lwhDQBggY477rjcf//9OfHEE7Nr1668/e1vz7vf/e5cdtllh9yn3Z0AAAtUVTnxxBOTjC4JsmvXrgUfJyekAQAsgj179uTCCy/M6aefnne961259NJLF9SfkAYAsAiOOeaYPPLII5mZmcnGjRvz6KOPLqg/IQ0AYBG97nWvyzve8Y7ce++9C+pHSAMAWKDt27fn+9//fpLkRz/6Uf7iL/4ib3nLWxbUp7M7AYBlZT6XzFhs27Zty9q1a7Nnz5785Cc/ya/+6q/mPe95z4L6FNIAABboF37hF/Lwww8vap92dwIAdEhIAwDokJAGANAhIQ0AoENCGgBAh4Q0AIAOuQQHALC8PLu4l8LIG3/pgKs8/fTTue666/KP//iPec1rXpP169fnhhtuWNDHCmkAAAu0YsWKfPKTn8xFF12UF198MRdffHHe9a535bzzzjvkPu3uBABYoDPOOCMXXXRRkuSkk07KW9/61jzzzDML6lNIAwBYRE899VQefvjhXHrppQvqR0gDAFgkP/zhD/P+978/n/70p3PyyScvqK95hbSqeqqqNlfVI1W1aaitqqoNVbVleD5lbP0PV9XWqnqiqq4cq1889LO1qj5TVTXUj6uqO4b6A1W1eqzN2uEztlTV2gV9WwCAw2TXrl15//vfn1//9V/Pr/zKryy4v4OZSfvnrbULW2vTw+sbk9zXWjs3yX3D61TVeUnWJDk/yVVJPldVxwxtPp9kfZJzh8dVQ31dkhdaa+ck+VSSTwx9rUpyU5JLk1yS5KbxMAgA0IPWWtatW5e3vvWt+Z3f+Z1F6XMhZ3dek+Qdw/ItSb6e5D8N9dtbay8nebKqtia5pKqeSnJya+0bSVJVX07yviT3DG1+b+jrK0k+O8yyXZlkQ2ttx9BmQ0bB7rYFbDcAsJzN45IZi+2v//qvc+utt+aCCy7IhRdemCT5gz/4g1x99dWH3Od8Q1pL8udV1ZL8j9bazUne0FrbliSttW1Vdfqw7plJ/mas7cxQ2zUs71vf2+bpoa/dVbUzyanj9Vna/FRVrc9ohi5vetOb5vmVAAAWx9vf/va01ha1z/mGtLe11p4dgtiGqvrWHOvWLLU2R/1Q2/ysMAqNNyfJ9PT04o4QAMAEzOuYtNbas8Pzc0n+LKPjw75bVWckyfD83LD6TJKzxppPJXl2qE/NUn9Fm6pakWRlkh1z9AUAsKwdMKRV1c9V1Ul7l5NckeTRJHcn2Xu25dokdw3LdydZM5yxeXZGJwhsHHaNvlhVlw3Hm123T5u9fV2b5P42mjP8WpIrquqU4YSBK4YaAMCyNp/dnW9I8mfD1TJWJPnfrbV7q+rBJHdW1bok30nygSRprT1WVXcmeTzJ7iTXt9b2DH19KMmXkpyQ0QkD9wz1LyS5dTjJYEdGZ4emtbajqj6W5MFhvY/uPYkAAGA5O2BIa619O8kvzlJ/Psk799Pm95P8/iz1TUl+fpb6SxlC3izvfTHJFw+0nQAAy4k7DgAAdGgh10kDAOjOY88/tqj9nX/q+Qdc54Mf/GC++tWv5vTTT8+jjz66KJ9rJg0AYIF+4zd+I/fee++i9imkAQAs0OWXX55Vq1Ytap9CGgBAh4Q0AIAOCWkAAB0S0gAAOuQSHIdo88zOV9UumFo5gS0BAMbN55IZi+3Xfu3X8vWvfz3f+973MjU1lY985CNZt27dgvoU0gAAFui2225b9D7t7gQA6JCQBgDQISENAKBDQhoAQIeENACADglpAAAdcgkOAGBZ+dGjjy1qfyf8/Pyuu3bvvffmhhtuyJ49e/Kbv/mbufHGGxf0uWbSAAAWaM+ePbn++utzzz335PHHH89tt92Wxx9/fEF9CmkAAAu0cePGnHPOOXnzm9+c1772tVmzZk3uuuuuBfUppAEALNAzzzyTs84666evp6am8swzzyyoTyENAGCBWmuvqlXVgvoU0gAAFmhqaipPP/30T1/PzMzkjW9844L6FNIAABbol3/5l7Nly5Y8+eST+fGPf5zbb789733vexfUp0twAADLynwvmbGYVqxYkc9+9rO58sors2fPnnzwgx/M+ecvbDuENACARXD11Vfn6quvXrT+7O4EAOiQkAYA0CEhDQBY8ma7BEZPDmX7hDQAYEk7/vjj8/zzz3cb1Fpref7553P88ccfVDsnDgAAS9rU1FRmZmayffv2SW/Kfh1//PGZmpo6qDZCGgCwpB177LE5++yzJ70Zi87uTgCADglpAAAdEtIAADokpAEAdEhIAwDokJAGANAhIQ0AoENCGgBAh4Q0AIAOCWkAAB0S0gAAOiSkAQB0SEgDAOiQkAYA0CEhDQCgQ0IaAECHhDQAgA4JaQAAHRLSAAA6JKQBAHRISAMA6JCQBgDQISENAKBDQhoAQIeENACADglpAAAdEtIAADokpAEAdEhIAwDokJAGANCheYe0qjqmqh6uqq8Or1dV1Yaq2jI8nzK27oeramtVPVFVV47VL66qzcN7n6mqGurHVdUdQ/2Bqlo91mbt8BlbqmrtonxrAIDOHcxM2g1Jvjn2+sYk97XWzk1y3/A6VXVekjVJzk9yVZLPVdUxQ5vPJ1mf5NzhcdVQX5fkhdbaOUk+leQTQ1+rktyU5NIklyS5aTwMAgAsV/MKaVU1leRfJfmjsfI1SW4Zlm9J8r6x+u2ttZdba08m2Zrkkqo6I8nJrbVvtNZaki/v02ZvX19J8s5hlu3KJBtaaztaay8k2ZCfBTsAgGVrvjNpn07yH5P8ZKz2htbatiQZnk8f6mcmeXpsvZmhduawvG/9FW1aa7uT7Exy6hx9vUJVra+qTVW1afv27fP8SgAA/TpgSKuq9yR5rrX20Dz7rFlqbY76obb5WaG1m1tr06216dNOO22emwkA0K/5zKS9Lcl7q+qpJLcn+RdV9b+SfHfYhZnh+blh/ZkkZ421n0ry7FCfmqX+ijZVtSLJyiQ75ugLAGBZO2BIa619uLU21VpbndEJAfe31v5NkruT7D3bcm2Su4blu5OsGc7YPDujEwQ2DrtEX6yqy4bjza7bp83evq4dPqMl+VqSK6rqlOGEgSuGGgDAsrZiAW0/nuTOqlqX5DtJPpAkrbXHqurOJI8n2Z3k+tbanqHNh5J8KckJSe4ZHknyhSS3VtXWjGbQ1gx97aiqjyV5cFjvo621HQvYZgCAJaFGE1bLx/T0dNu0adNh/5zNMztfVbtgauVh/1wAYPmoqodaa9OzveeOAwAAHRLSAAA6JKQBAHRISAMA6JCQBgDQISENAKBDQhoAQIeENACADglpAAAdEtIAADokpAEAdEhIAwDokJAGANAhIQ0AoENCGgBAh4Q0AIAOCWkAAB0S0gAAOiSkAQB0SEgDAOiQkAYA0CEhDQCgQ0IaAECHVkx6A5aTzTM7X1W7YGrlBLYEAFjqzKQBAHRISAMA6JCQBgDQISENAKBDQhoAQIeENACADglpAAAdEtIAADokpAEAdEhIAwDokJAGANAhIQ0AoENCGgBAh4Q0AIAOCWkAAB0S0gAAOiSkAQB0SEgDAOiQkAYA0CEhDQCgQ0IaAECHhDQAgA4JaQAAHRLSAAA6JKQBAHRISAMA6JCQBgDQISENAKBDKya9Acvd5pmdr6pdMLVyAlsCACwlZtIAADokpAEAdEhIAwDokJAGANAhIQ0AoENCGgBAh4Q0AIAOHTCkVdXxVbWxqv6uqh6rqo8M9VVVtaGqtgzPp4y1+XBVba2qJ6rqyrH6xVW1eXjvM1VVQ/24qrpjqD9QVavH2qwdPmNLVa1d1G8PANCp+cykvZzkX7TWfjHJhUmuqqrLktyY5L7W2rlJ7htep6rOS7ImyflJrkryuao6Zujr80nWJzl3eFw11NcleaG1dk6STyX5xNDXqiQ3Jbk0ySVJbhoPgwAAy9UBQ1ob+eHw8tjh0ZJck+SWoX5LkvcNy9ckub219nJr7ckkW5NcUlVnJDm5tfaN1lpL8uV92uzt6ytJ3jnMsl2ZZENrbUdr7YUkG/KzYAcAsGzN65i0qjqmqh5J8lxGoemBJG9orW1LkuH59GH1M5M8PdZ8ZqidOSzvW39Fm9ba7iQ7k5w6R1/7bt/6qtpUVZu2b98+n68EANC1eYW01tqe1tqFSaYymhX7+TlWr9m6mKN+qG3Gt+/m1tp0a236tNNOm2PTAACWhoM6u7O19v0kX89ol+N3h12YGZ6fG1abSXLWWLOpJM8O9alZ6q9oU1UrkqxMsmOOvgAAlrX5nN15WlW9blg+Icm/TPKtJHcn2Xu25dokdw3LdydZM5yxeXZGJwhsHHaJvlhVlw3Hm123T5u9fV2b5P7huLWvJbmiqk4ZThi4YqgBACxrK+axzhlJbhnO0HxNkjtba1+tqm8kubOq1iX5TpIPJElr7bGqujPJ40l2J7m+tbZn6OtDSb6U5IQk9wyPJPlCkluramtGM2hrhr52VNXHkjw4rPfR1tqOhXxhAICloEYTVsvH9PR027Rp02H/nM0zOw+57QVTKxdxSwCApaqqHmqtTc/2njsOAAB0SEgDAOiQkAYA0CEhDQCgQ0IaAECHhDQAgA4JaQAAHZrPxWxZZLNdY8210wCAcWbSAAA6JKQBAHRISAMA6JCQBgDQISENAKBDQhoAQIeENACADglpAAAdEtIAADokpAEAdEhIAwDokJAGANAhIQ0AoENCGgBAh4Q0AIAOCWkAAB0S0gAAOiSkAQB0SEgDAOjQiklvACObZ3a+qnbB1MoJbAkA0AMzaQAAHRLSAAA6JKQBAHRISAMA6JCQBgDQISENAKBDQhoAQIeENACADglpAAAdEtIAADokpAEAdEhIAwDokJAGANAhIQ0AoEMrJr0B7N/mmZ2vql0wtXICWwIAHGlm0gAAOiSkAQB0SEgDAOiQkAYA0CEhDQCgQ0IaAECHhDQAgA4JaQAAHRLSAAA6JKQBAHRISAMA6JCQBgDQISENAKBDQhoAQIdWTHoDODibZ3a+qnbB1MoJbAkAcDiZSQMA6JCQBgDQISENAKBDBwxpVXVWVf1lVX2zqh6rqhuG+qqq2lBVW4bnU8bafLiqtlbVE1V15Vj94qraPLz3maqqoX5cVd0x1B+oqtVjbdYOn7GlqtYu6rcHAOjUfGbSdif59621tya5LMn1VXVekhuT3NdaOzfJfcPrDO+tSXJ+kquSfK6qjhn6+nyS9UnOHR5XDfV1SV5orZ2T5FNJPjH0tSrJTUkuTXJJkpvGwyAAwHJ1wJDWWtvWWvvbYfnFJN9McmaSa5LcMqx2S5L3DcvXJLm9tfZya+3JJFuTXFJVZyQ5ubX2jdZaS/Llfdrs7esrSd45zLJdmWRDa21Ha+2FJBvys2AHALBsHdQxacNuyF9K8kCSN7TWtiWjIJfk9GG1M5M8PdZsZqidOSzvW39Fm9ba7iQ7k5w6R1/7btf6qtpUVZu2b99+MF8JAKBL8w5pVXVikj9J8tuttR/MteostTZH/VDb/KzQ2s2ttenW2vRpp502x6YBACwN8wppVXVsRgHtj1trfzqUvzvswszw/NxQn0ly1ljzqSTPDvWpWeqvaFNVK5KsTLJjjr4AAJa1+ZzdWUm+kOSbrbU/HHvr7iR7z7Zcm+Susfqa4YzNszM6QWDjsEv0xaq6bOjzun3a7O3r2iT3D8etfS3JFVV1ynDCwBVDDQBgWZvPbaHeluTfJtlcVY8Mtf+c5ONJ7qyqdUm+k+QDSdJae6yq7kzyeEZnhl7fWtsztPtQki8lOSHJPcMjGYXAW6tqa0YzaGuGvnZU1ceSPDis99HW2o5D+6oAAEtHjSaslo/p6em2adOmw/45s91Dsyfu5wkA/auqh1pr07O9544DAAAdEtIAADokpAEAdEhIAwDokJAGANAhIQ0AoENCGgBAh4Q0AIAOCWkAAB0S0gAAOiSkAQB0aD43WGcJmu3eou7nCQBLh5k0AIAOCWkAAB0S0gAAOiSkAQB0SEgDAOiQkAYA0CEhDQCgQ0IaAECHhDQAgA4JaQAAHRLSAAA65N6dR5HZ7ueZuKcnAPTITBoAQIeENACADglpAAAdEtIAADokpAEAdEhIAwDokJAGANAhIQ0AoENCGgBAh9xxgFnvROAuBAAwWWbSAAA6JKQBAHRISAMA6JCQBgDQISENAKBDQhoAQIeENACADrlOGrNy7TQAmCwzaQAAHRLSAAA6JKQBAHRISAMA6JCQBgDQISENAKBDQhoAQIdcJ415c+00ADhyzKQBAHRISAMA6JCQBgDQISENAKBDQhoAQIeENACADrkEBwvishwAcHiYSQMA6JCQBgDQISENAKBDBwxpVfXFqnquqh4dq62qqg1VtWV4PmXsvQ9X1daqeqKqrhyrX1xVm4f3PlNVNdSPq6o7hvoDVbV6rM3a4TO2VNXaRfvWAACdm89M2peSXLVP7cYk97XWzk1y3/A6VXVekjVJzh/afK6qjhnafD7J+iTnDo+9fa5L8kJr7Zwkn0ryiaGvVUluSnJpkkuS3DQeBgEAlrMDhrTW2l8l2bFP+ZoktwzLtyR531j99tbay621J5NsTXJJVZ2R5OTW2jdaay3Jl/dps7evryR55zDLdmWSDa21Ha21F5JsyKvDIgDAsnSol+B4Q2ttW5K01rZV1elD/cwkfzO23sxQ2zUs71vf2+bpoa/dVbUzyanj9VnavEJVrc9oli5vetObDvErsVhclgMAFm6xTxyoWWptjvqhtnllsbWbW2vTrbXp0047bV4bCgDQs0MNad8ddmFmeH5uqM8kOWtsvakkzw71qVnqr2hTVSuSrMxo9+r++gIAWPYONaTdnWTv2ZZrk9w1Vl8znLF5dkYnCGwcdo2+WFWXDcebXbdPm719XZvk/uG4ta8luaKqThlOGLhiqAEALHsHPCatqm5L8o4kr6+qmYzOuPx4kjural2S7yT5QJK01h6rqjuTPJ5kd5LrW2t7hq4+lNGZoickuWd4JMkXktxaVVszmkFbM/S1o6o+luTBYb2Pttb2PYEBAGBZqtGk1fIxPT3dNm3adNg/Z7aD49k/Jw4AwKtV1UOttenZ3nPHAQCADh3qJTjgoLgsBwAcHDNpAAAdEtIAADokpAEAdEhIAwDokJAGANAhZ3cyMc74BID9M5MGANAhIQ0AoENCGgBAh4Q0AIAOOXGArjiZAABGzKQBAHRISAMA6JCQBgDQIcek0T3HqQFwNDKTBgDQISENAKBDQhoAQIcck8aS5Dg1AJY7M2mH6Ns/+NakNwEAWMaENACADtndybJhFygAy4mZNACADglpAAAdsruTZc0uUACWKjNpAAAdMpPGUcfsGgBLgZk0AIAOmUmDzD67lphhA2ByzKQBAHTITBrMwfFrAEyKmTQAgA6ZSYODZHYNgCPBTBoAQIfMpMEiMLsGwGIzkwYA0CEzaXCYmF0DYCHMpAEAdMhMGhxBZtcAmC8hDSZMcANgNkIadEhwA0BIgyVCcAM4ughpsIQJbgDLl5AGy4zgBrA8CGlwFBDcAJYeIQ2OUrMFt/0R6ACOPCENOCAzcQBHnpAGHJL5zsQJcwCHRkgDDithDuDQCGlAF4Q5gFcS0oAlRZgDjhZCGrAsHczZq/sS8IAeCGkA+3B5EqAHQhrAApixAw4XIQ1gQgQ8YC5CGsASJODB8iekARxlnCELS4OQdoiO+/7WfDvJm09+y6Q3BeCwWMhsXSLkwUIJaQAcFmbsYGGWREirqquS/NckxyT5o9baxye8SQAsEmEOZtd9SKuqY5L8tyTvSjKT5MGquru19vgkt+u1T23Lyxeek2//4Ft2eQIcAU6W4GjTfUhLckmSra21bydJVd2e5JokEw1p4wQ1gL4t9Pi6xSY0Mh9LIaSdmeTpsdczSS4dX6Gq1idZP7z8YVU9cQS26/VJvncEPodXM/aTZfwnx9hPjrGfrOU8/v9kf28shZBWs9TaK160dnOSm4/M5oxU1abW2vSR/ExGjP1kGf/JMfaTY+wn62gd/9dMegPmYSbJWWOvp5I8O6FtAQA4IpZCSHswyblVdXZVvTbJmiR3T3ibAAAOq+53d7bWdlfVbyX5WkaX4Phia+2xCW9WcoR3r/IKxn6yjP/kGPvJMfaTdVSOf7XWDrwWAABH1FLY3QkAcNQR0gAAOiSkHaSquqqqnqiqrVV146S3Z7moqqeqanNVPVJVm4baqqraUFVbhudTxtb/8PA3eKKqrhyrXzz0s7WqPlNVs13C5ahXVV+squeq6tGx2qKNd1UdV1V3DPUHqmr1Ef2CHdvP2P9eVT0z/P4fqaqrx94z9oukqs6qqr+sqm9W1WNVdcNQ99s/zOYYe7/9ubTWPOb5yOjEhX9I8uYkr03yd0nOm/R2LYdHkqeSvH6f2n9JcuOwfGOSTwzL5w1jf1ySs4e/yTHDexuT/LOMrq93T5J3T/q79fhIcnmSi5I8ejjGO8m/S/Lfh+U1Se6Y9Hfu5bGfsf+9JP9hlnWN/eKO/RlJLhqWT0ryf4cx9tuf3Nj77c/xMJN2cH56i6rW2o+T7L1FFYfHNUluGZZvSfK+sfrtrbWXW2tPJtma5JKqOiPJya21b7TRv9Ivj7VhTGvtr5Ls2Ke8mOM93tdXkrzTrObIfsZ+f4z9ImqtbWut/e2w/GKSb2Z0Vxu//cNsjrHfH2MfuzsP1my3qJrrR8b8tSR/XlUP1eg2X0nyhtbatmT0DzzJ6UN9f3+HM4flfevMz2KO90/btNZ2J9mZ5NTDtuXLw29V1d8Pu0P37m4z9ofJsCvsl5I8EL/9I2qfsU/89vdLSDs4B7xFFYfsba21i5K8O8n1VXX5HOvu7+/g73N4HMp4+1scnM8n+adJLkyyLcknh7qxPwyq6sQkf5Lkt1trP5hr1Vlqxn8BZhl7v/05CGkHxy2qDpPW2rPD83NJ/iyjXcvfHaa2Mzw/N6y+v7/DzLC8b535Wczx/mmbqlqRZGXmv4vvqNNa+25rbU9r7SdJ/mdGv//E2C+6qjo2o5Dwx621Px3KfvtHwGxj77c/NyHt4LhF1WFQVT9XVSftXU5yRZJHMxrbtcNqa5PcNSzfnWTNcCbP2UnOTbJx2E3xYlVdNhyHcN1YGw5sMcd7vK9rk9w/HD/CLPYGhMG/zuj3nxj7RTWM1ReSfLO19odjb/ntH2b7G3u//QOY9JkLS+2R5OqMzkr5hyS/O+ntWQ6PjM6W/bvh8djecc3oWIL7kmwZnleNtfnd4W/wRMbO4EwyndE/8n9I8tkMd9XweNWY35bRroVdGf3f57rFHO8kxyf5Pxkd7LsxyZsn/Z17eexn7G9NsjnJ32f0H5ozjP1hGfu3Z7T76++TPDI8rvbbn+jY++3P8XBbKACADtndCQDQISENAKBDQhoAQIeENACADglpAAAdEtIAADokpAEAdOj/AxgiJ8yHJqbRAAAAAElFTkSuQmCC", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAnMAAAGbCAYAAACvTxZ8AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAesElEQVR4nO3dfYyd1X0n8O8PjA0qGDApiDBUBoGWlyQlgQJSoigtDRAahaghlaO0uIorpIiVqFptF/oPSyKq5I822SibaFGDQmgLQWkjULSQutCo2ioLOIUuL3HKrPGGATckGGOyCgacs3/MYzKYsT3Y47lz7nw+0tV97u8+5zxnOEL+6jwvt1prAQCgT4eMegAAAOw/YQ4AoGPCHABAx4Q5AICOCXMAAB1bNuoBzLe3vOUtbfXq1aMeBgDAPn3ve9/7SWvtlw+kj7ELc6tXr86GDRtGPQwAgH2qqv97oH04zQoA0DFhDgCgY8IcAEDHxu6aOQBgaXnllVcyNTWVl156adRD2aPDDz88ExMTOeyww+a9b2EOAOja1NRUjjrqqKxevTpVNerhvEFrLc8991ympqZyyimnzHv/TrMCAF176aWXctxxxy3KIJckVZXjjjvuoK0cCnMAQPcWa5Db5WCOT5gDAOiYa+YAgLHyyNQL89rf2yeO3uc+L730Ut773vdmx44defXVV3PFFVfkhhtumNdx7IkwBwBwgFasWJH77rsvRx55ZF555ZW85z3vyQc+8IFceOGFB/3YTrMCABygqsqRRx6ZZPpRKa+88sqCXccnzAEAzIOdO3fmnHPOyfHHH5/3v//9ueCCCxbkuMIcAMA8OPTQQ/Pwww9namoqDzzwQB599NEFOa4wBwAwj4455pi8733vyz333LMgxxPmAAAO0I9//ONs27YtSfKzn/0s//AP/5AzzjhjQY7tblYAYKzM5VEi823Lli1Zu3Ztdu7cmZ///Of5nd/5nXzwgx9ckGMvmTD32HOPvbZ99nFnj3AkAMC4ecc73pGHHnpoJMd2mhUAoGPCHABAx4Q5AICOCXMAAB0T5gAAOibMAQB0bMk8mgQAWCKemedHhLz1nfvc5amnnsqVV16Zf//3f88hhxySq666Ktdcc838jmMPlmSY88w5AGA+LVu2LH/+53+ed73rXXnxxRdz7rnn5v3vf3/OOuusg35sp1kBAA7QiSeemHe9611JkqOOOipnnnlmnn766QU5tjAHADCPNm/enIceeigXXHDBghxPmAMAmCc//elP85GPfCSf//zns3LlygU55pzCXFVtrqpHqurhqtow1FZV1fqqemJ4P3bG/tdV1WRV/aCqLplRP3foZ7KqvlBVNdRXVNXXh/r9VbV6Rpu1wzGeqKq18/aXAwDMo1deeSUf+chH8vGPfzy//du/vWDHfTMrc7/eWjuntXbe8PnaJPe21k5Pcu/wOVV1VpI1Sc5OcmmSL1XVoUObLye5Ksnpw+vSob4uyfOttdOSfC7JZ4e+ViW5PskFSc5Pcv3M0AgAsBi01rJu3bqceeaZ+aM/+qMFPfaB3M16eZL3Ddu3JPlOkv881G9vre1I8mRVTSY5v6o2J1nZWvtuklTV15J8OMndQ5v/MvT1jSRfHFbtLkmyvrW2dWizPtMB8LYDGDcAMM7m8CiR+fbP//zPufXWW/P2t78955xzTpLkz/7sz3LZZZcd9GPPNcy1JH9fVS3Jf2+t3ZTkhNbaliRprW2pquOHfU9K8r9mtJ0aaq8M27vXd7V5aujr1ap6IclxM+uztHlNVV2V6RW//Mqv/Moc/yQAgPnxnve8J621kRx7rmHu3a21Z4bAtr6qNu5l35ql1vZS3982vyhMh8ubkuS8884bzX9JAIARmNM1c621Z4b3Z5N8M9PXr/2oqk5MkuH92WH3qSQnz2g+keSZoT4xS/11bapqWZKjk2zdS18AAGQOYa6qfqmqjtq1neTiJI8muSvJrrtL1ya5c9i+K8ma4Q7VUzJ9o8MDwynZF6vqwuF6uCt3a7OrryuS3Nem1yq/neTiqjp2uPHh4qEGAEDmdpr1hCTfHJ4isizJ37TW7qmqB5PcUVXrkvwwyUeTpLX2WFXdkeTxJK8mubq1tnPo65NJvprkiEzf+HD3UP9KkluHmyW2Zvpu2LTWtlbVp5M8OOz3qV03QwAAMIcw11rblORXZ6k/l+SiPbS5McmNs9Q3JHnbLPWXMoTBWb67OcnN+xonAMBS5BcgAAA6diDPmQMAWHQee+6xee3v7OPO3uc+n/jEJ/Ktb30rxx9/fB599NF5Pf6+WJkDADhAv//7v5977rlnJMcW5gAADtB73/verFq1aiTHFuYAADomzAEAdEyYAwDomDAHANCxsXs0yc9e3plHpl54Q33T9p/Ouv/Pf/bGfQEOlrdPHD3qIcDYm8ujRObbxz72sXznO9/JT37yk0xMTOSGG27IunXrFuTYYxfmAAAW2m233TayYzvNCgDQsbFemdu0feOohwAAcFBZmQMA6JgwBwDQMWEOAKBjY33N3FzMvK7u1JVnjHAkAABv3pIPcwDAePnZo4/Na39HvG1uz6275557cs0112Tnzp35gz/4g1x77bXzOo49cZoVAOAA7dy5M1dffXXuvvvuPP7447ntttvy+OOPL8ixhTkAgAP0wAMP5LTTTsupp56a5cuXZ82aNbnzzjsX5NjCHADAAXr66adz8sknv/Z5YmIiTz/99IIcW5gDADhArbU31KpqQY4tzAEAHKCJiYk89dRTr32emprKW9/61gU5tjAHAHCAfu3Xfi1PPPFEnnzyybz88su5/fbb86EPfWhBju3RJADAWJnro0Tm07Jly/LFL34xl1xySXbu3JlPfOITOfvshRmHMAcAMA8uu+yyXHbZZQt+XKdZAQA6JswBAHRMmAMAujfbo0EWk4M5PmEOAOja4Ycfnueee27RBrrWWp577rkcfvjhB6V/N0AAAF2bmJjI1NRUfvzjH496KHt0+OGHZ2Ji4qD0LcwBAF077LDDcsopp4x6GCPjNCsAQMeEOQCAjglzAAAdE+YAADomzAEAdEyYAwDomDAHANAxYQ4AoGPCHABAx4Q5AICOCXMAAB0T5gAAOibMAQB0TJgDAOiYMAcA0DFhDgCgY8IcAEDHhDkAgI4JcwAAHVs26gEsJpu2b3xt+9SVZ4xwJAAAc2NlDgCgY8IcAEDHhDkAgI4JcwAAHRPmAAA6JswBAHRMmAMA6JgwBwDQMWEOAKBjwhwAQMfmHOaq6tCqeqiqvjV8XlVV66vqieH92Bn7XldVk1X1g6q6ZEb93Kp6ZPjuC1VVQ31FVX19qN9fVatntFk7HOOJqlo7L381AMCYeDMrc9ck+f6Mz9cmube1dnqSe4fPqaqzkqxJcnaSS5N8qaoOHdp8OclVSU4fXpcO9XVJnm+tnZbkc0k+O/S1Ksn1SS5Icn6S62eGRgCApW5OYa6qJpL8VpK/nFG+PMktw/YtST48o357a21Ha+3JJJNJzq+qE5OsbK19t7XWknxttza7+vpGkouGVbtLkqxvrW1trT2fZH1+EQABAJa8ua7MfT7JnyT5+YzaCa21LUkyvB8/1E9K8tSM/aaG2knD9u7117Vprb2a5IUkx+2lr9epqquqakNVbXh+63Nz/JMAAPq3zzBXVR9M8mxr7Xtz7LNmqbW91Pe3zS8Krd3UWjuvtXbesauOm+MwAQD6N5eVuXcn+VBVbU5ye5LfqKq/SvKj4dRphvdnh/2nkpw8o/1EkmeG+sQs9de1qaplSY5OsnUvfQEAkDmEudbada21idba6kzf2HBfa+13k9yVZNfdpWuT3Dls35VkzXCH6imZvtHhgeFU7ItVdeFwPdyVu7XZ1dcVwzFakm8nubiqjh1ufLh4qAEAkGTZAbT9TJI7qmpdkh8m+WiStNYeq6o7kjye5NUkV7fWdg5tPpnkq0mOSHL38EqSryS5taomM70it2boa2tVfTrJg8N+n2qtbT2AMQMAjJWaXgAbH2e/453t9v/xnSTJpu0b97ufU1eeMU8jAviFt08cPeohAItIVX2vtXbegfThFyAAADomzAEAdGxJhbkV2yazYtvkqIcBADBvllSYAwAYN0syzFmdAwDGxZIMcwAA40KYAwDomDAHANAxYQ4AoGPCHABAx4Q5AICOCXMAAB0T5gAAOrZkwpwHBQMA42jJhDkAgHEkzAEAdEyYAwDomDAHANCxJRvmVmybdFMEANC9JRvmAADGgTAHANAxYQ4AoGNjHeZW/J/N+9120/aN2bR94/wNBgDgIFg26gEALCWPTL0w6iEAY2asV+YAAMbd+Ia5f3OKFAAYf+Mb5gAAloCxD3MHchMEAMBiN/ZhDgBgnAlzAAAdE+YAADomzAEAdEyYAwDomDAHANAxYQ4AoGPCHABAx5Z8mFuxbXLUQwAA2G9LIswJbADAuFoSYW755i2jHgIAwEGxJMIcAMC4EuYAADomzAEAdEyYAwDomDAHANAxYQ4AoGPCHABAx4Q5AICOCXMAAB0bzzD3bxtHPQIAgAUxnmEOAGCJEOYAADomzAEAdEyYAwDo2JIJc8s3bxn1EAAA5t2SCXMAAONImAMA6Jgwl2TFtslRDwEAYL8IcwAAHRPmAAA6JswBAHRsn2Guqg6vqgeq6l+r6rGqumGor6qq9VX1xPB+7Iw211XVZFX9oKoumVE/t6oeGb77QlXVUF9RVV8f6vdX1eoZbdYOx3iiqtbO618PANC5uazM7UjyG621X01yTpJLq+rCJNcmube1dnqSe4fPqaqzkqxJcnaSS5N8qaoOHfr6cpKrkpw+vC4d6uuSPN9aOy3J55J8duhrVZLrk1yQ5Pwk188MjQAAS90+w1yb9tPh42HDqyW5PMktQ/2WJB8eti9PcntrbUdr7ckkk0nOr6oTk6xsrX23tdaSfG23Nrv6+kaSi4ZVu0uSrG+tbW2tPZ9kfX4RAAEAlrw5XTNXVYdW1cNJns10uLo/yQmttS1JMrwfP+x+UpKnZjSfGmonDdu711/XprX2apIXkhy3l752H99VVbWhqjY8v/W5ufxJAABjYU5hrrW2s7V2TpKJTK+yvW0vu9dsXeylvr9tZo7vptbaea21845dddxehgYAMF7e1N2srbVtSb6T6VOdPxpOnWZ4f3bYbSrJyTOaTSR5ZqhPzFJ/XZuqWpbk6CRb99IXAACZ292sv1xVxwzbRyT5zSQbk9yVZNfdpWuT3Dls35VkzXCH6imZvtHhgeFU7ItVdeFwPdyVu7XZ1dcVSe4brqv7dpKLq+rY4caHi4fagtm0feNrLwCAxWbZHPY5Mcktwx2phyS5o7X2rar6bpI7qmpdkh8m+WiStNYeq6o7kjye5NUkV7fWdg59fTLJV5MckeTu4ZUkX0lya1VNZnpFbs3Q19aq+nSSB4f9PtVa23ogfzAAwDjZZ5hrrf3vJO+cpf5ckov20ObGJDfOUt+Q5A3X27XWXsoQBmf57uYkN+9rnHOxfPOWvLz6xPnoCgBgUfALEAAAHRPmAAA6JswBAHRMmAMA6JgwBwDQMWEOAKBjwhwAQMeEOQCAjglzAAAdE+YAADomzA1WbJvMim2Tox4GAMCbss/fZu3Z8hefGvUQAAAOKitzAAAdE+YAADomzAEAdEyYAwDo2JILc8s3bxn1EAAA5s3Y3c264+cv5en/tzlJsny0QwEAOOiW3MocAMA4EeYAADomzAEAdEyYAwDomDAHANAxYQ4AoGPCHABAx4Q5AICOCXMAAB0T5gAAOibMAQB0TJgDAOiYMAcA0DFhDgCgY8IcAEDHlmSYW755y6iHAAAwL5ZkmNubFdsmRz0EAIA5E+YAADomzAEAdEyYAwDomDAHANAxYQ4AoGPCHABAx5aNegA92bR942vbp648Y4QjAQCYZmUOAKBjwhwAQMeEOQCAjglzAAAdE+YAADomzAEAdEyYAwDo2NiFuUN2vDyn/ZZv3nKQRwIAcPCNXZgDAFhKhDkAgI4JcwAAHRPmAAA6JswBAHRMmAMA6JgwBwDQMWFuFiu2TWbFtslRDwMAYJ+EOQCAjglzAAAdE+YAADq2zzBXVSdX1T9W1fer6rGqumaor6qq9VX1xPB+7Iw211XVZFX9oKoumVE/t6oeGb77QlXVUF9RVV8f6vdX1eoZbdYOx3iiqtbO618PANC5uazMvZrkj1trZya5MMnVVXVWkmuT3NtaOz3JvcPnDN+tSXJ2kkuTfKmqDh36+nKSq5KcPrwuHerrkjzfWjstyeeSfHboa1WS65NckOT8JNfPDI0AAEvdPsNca21La+1fhu0Xk3w/yUlJLk9yy7DbLUk+PGxfnuT21tqO1tqTSSaTnF9VJyZZ2Vr7bmutJfnabm129fWNJBcNq3aXJFnfWtvaWns+yfr8IgACACx5b+qaueH05zuT3J/khNbalmQ68CU5ftjtpCRPzWg2NdROGrZ3r7+uTWvt1SQvJDluL33tPq6rqmpDVW144YXtb+ZPAgDo2pzDXFUdmeRvk/xha21vialmqbW91Pe3zS8Krd3UWjuvtXbe0Uev3MvQAADGy5zCXFUdlukg99ettb8byj8aTp1meH92qE8lOXlG84kkzwz1iVnqr2tTVcuSHJ1k6176AgAgc7ubtZJ8Jcn3W2t/MeOru5Lsurt0bZI7Z9TXDHeonpLpGx0eGE7FvlhVFw59Xrlbm119XZHkvuG6um8nubiqjh1ufLh4qAEAkGTZHPZ5d5LfS/JIVT081P40yWeS3FFV65L8MMlHk6S19lhV3ZHk8UzfCXt1a23n0O6TSb6a5Igkdw+vZDos3lpVk5lekVsz9LW1qj6d5MFhv0+11rbu358KADB+9hnmWmv/M7Nfu5YkF+2hzY1JbpylviHJ22apv5QhDM7y3c1Jbt7XOAEAliK/AAEA0LGxDXPLX3xq3zsBAHRubMMcAMBSIMwBAHRMmAMA6JgwBwDQMWEOAKBjwhwAQMeEOQCAjs3l57yWrBXbJrPjmNNm/W7T9o2vbZ+68oyFGhIAwOtYmQMA6NiSDnPLN28Z9RAAAA7Ikg5zAAC9E+YAADomzAEAdEyYAwDomDAHANAxYQ4AoGPCHABAx4Q5AICOCXMAAB0T5gAAOibMAQB0TJgDAOiYMAcA0DFhDgCgY0s+zC3fvGXUQwAA2G9LPswBAPRMmAMA6Jgwtw8rtk2OeggAAHskzAEAdEyYAwDomDAHANAxYQ4AoGPCHABAx4Q5AICOCXMAAB1bNuoBjINN2ze+tn3qyjNGOBIAYKmxMgcA0DFhDgCgY8JckuWbt4x6CAAA+0WYAwDomDAHANAxYQ4AoGPCHABAx4Q5AICOCXMAAB0T5gAAOibMzcGKbZNZsW1y1MMAAHgDYQ4AoGPCHABAx4Q5AICOCXMAAB0T5gbLN28Z9RAAAN40YQ4AoGPCHABAx4Q5AICOCXMAAB1bNuoBjJtN2ze+tn3qyjNGOBIAYCmwMgcA0LF9hrmqurmqnq2qR2fUVlXV+qp6Yng/dsZ311XVZFX9oKoumVE/t6oeGb77QlXVUF9RVV8f6vdX1eoZbdYOx3iiqtbO218NADAm5rIy99Ukl+5WuzbJva2105PcO3xOVZ2VZE2Ss4c2X6qqQ4c2X05yVZLTh9euPtcleb61dlqSzyX57NDXqiTXJ7kgyflJrp8ZGgEAmEOYa639U5Ktu5UvT3LLsH1Lkg/PqN/eWtvRWnsyyWSS86vqxCQrW2vfba21JF/brc2uvr6R5KJh1e6SJOtba1tba88nWZ83hkoAgCVtf6+ZO6G1tiVJhvfjh/pJSZ6asd/UUDtp2N69/ro2rbVXk7yQ5Li99PUGVXVVVW2oqg0vvLB9P/8kAID+zPcNEDVLre2lvr9tXl9s7abW2nmttfOOPnrlnAa6P1ZsmzxofQMA7I/9DXM/Gk6dZnh/dqhPJTl5xn4TSZ4Z6hOz1F/XpqqWJTk606d199TXQeP3WQGA3uxvmLsrya67S9cmuXNGfc1wh+opmb7R4YHhVOyLVXXhcD3clbu12dXXFUnuG66r+3aSi6vq2OHGh4uHGgAAg30+NLiqbkvyviRvqaqpTN9h+pkkd1TVuiQ/TPLRJGmtPVZVdyR5PMmrSa5ure0cuvpkpu+MPSLJ3cMrSb6S5Naqmsz0ityaoa+tVfXpJA8O+32qtbb7jRgAAEvaPsNca+1je/jqoj3sf2OSG2epb0jytlnqL2UIg7N8d3OSm/c1RgCApcovQAAAdEyYAwDomDAHANCxfV4zx/7btH3ja9unrjxjhCMBAMaVlTkAgI4JcwAAHRPmduNXIACAnghzAAAdE+YAADomzL1JK7ZNZsW2yVEPAwAgiTAHANA1YQ4AoGPCHABAx4S5BbJp+8bX/SIEAMB8EOYAADo2dmGudr6c5S8+NephAAAsiLELc/PBr0AAAL0Q5gAAOibMAQB0TJgDAOiYMLef/KQXALAYLBv1AJaamc+aO3XlGSMcCQAwDqzMAQB0TJgDAOiYMLcHnjUHAPRAmAMA6JgwBwDQMWEOAKBjHk0yQh5TAgAcKCtzB8CDgwGAURPmAAA6JszthceTAACLnTAHANAxN0AsEm6GAAD2h5U5AICOCXMAAB0T5gAAOuaauX1YvnlLXl594h6/3/WsuR3HnDZvx3T9HAAwV1bmAAA6JswBAHRMmAMA6Jhr5uZgX9fNHUyunwMA9sbKHABAx4S5ebLrrlYAgIXkNGtHnHIFAHZnZQ4AoGNW5jpllQ4ASKzMzasV2yZdOwcALChhbo6Wb94y6iHs0abtG1+3UgcALB3CHABAx1wzN0ZcRwcAS4+VuTdhrqdaXTcHACwUK3NjyiodACwNVubepB5X53bdIOEmCQAYP1bm9sPyzVvy8uoTRz2M/WLFDgDGizC3hAl2ANA/YW4/zWV1btep1h3HnLYQQzoggh0A9EmYOwBzPd26YttkF4Ful9murRPwAGBxEuYO0K4bIsZplW42e7p5QsgDgNES5ubJmwl1vQa62czlDlmBDwAOHmFuns18dMmegt24Bbp92VfgE/YAYP91Eeaq6tIk/zXJoUn+srX2mREPaU5meybdroA38zl0SynYzeZgPv9OUARg3C36MFdVhyb5b0nen2QqyYNVdVdr7fHRjmz/7B7wXl594hseMLzUw9186ulByYInAPtj0Ye5JOcnmWytbUqSqro9yeVJugxzu5tt9W559ryiNxvhbzz0FDzng/AKMD96CHMnJXlqxuepJBfM3KGqrkpy1fBxx0Uf/0+PLtDYmH9vSfKTUQ+C/Wb++mXu+mb++vUfDrSDHsJczVJrr/vQ2k1JbkqSqtrQWjtvIQbG/DN/fTN//TJ3fTN//aqqDQfaxyHzMZCDbCrJyTM+TyR5ZkRjAQBYVHoIcw8mOb2qTqmq5UnWJLlrxGMCAFgUFv1p1tbaq1X1H5N8O9OPJrm5tfbYXprctDAj4yAxf30zf/0yd30zf/064Lmr1tq+9wIAYFHq4TQrAAB7IMwBAHRsrMJcVV1aVT+oqsmqunbU4+GNqurmqnq2qh6dUVtVVeur6onh/dgZ3103zOcPquqS0YyaJKmqk6vqH6vq+1X1WFVdM9TN3yJXVYdX1QNV9a/D3N0w1M1dR6rq0Kp6qKq+NXw2f52oqs1V9UhVPbzrUSTzOX9jE+Zm/OzXB5KcleRjVXXWaEfFLL6a5NLdatcmube1dnqSe4fPGeZvTZKzhzZfGuaZ0Xg1yR+31s5McmGSq4c5Mn+L344kv9Fa+9Uk5yS5tKoujLnrzTVJvj/js/nry6+31s6Z8TzAeZu/sQlzmfGzX621l5Ps+tkvFpHW2j8l2bpb+fIktwzbtyT58Iz67a21Ha21J5NMZnqeGYHW2pbW2r8M2y9m+h+Vk2L+Fr027afDx8OGV4u560ZVTST5rSR/OaNs/vo2b/M3TmFutp/9OmlEY+HNOaG1tiWZDgxJjh/q5nSRqqrVSd6Z5P6Yvy4Mp+geTvJskvWtNXPXl88n+ZMkP59RM3/9aEn+vqq+N/wEaTKP87fonzP3JuzzZ7/ojjldhKrqyCR/m+QPW2vbq2abpuldZ6mZvxFpre1Mck5VHZPkm1X1tr3sbu4Wkar6YJJnW2vfq6r3zaXJLDXzN1rvbq09U1XHJ1lfVRv3su+bnr9xWpnzs1/9+lFVnZgkw/uzQ92cLjJVdVimg9xft9b+biibv4601rYl+U6mr8Uxd314d5IPVdXmTF9C9BtV9Vcxf91orT0zvD+b5JuZPm06b/M3TmHOz371664ka4fttUnunFFfU1UrquqUJKcneWAE4yNJTS/BfSXJ91trfzHjK/O3yFXVLw8rcqmqI5L8ZpKNMXddaK1d11qbaK2tzvS/bfe11n435q8LVfVLVXXUru0kFyd5NPM4f2NzmnU/fvaLEaiq25K8L8lbqmoqyfVJPpPkjqpal+SHST6aJK21x6rqjiSPZ/pOyquHU0WMxruT/F6SR4Zrr5LkT2P+enBikluGO+IOSXJHa+1bVfXdmLue+X+vDydk+tKGZDp3/U1r7Z6qejDzNH9+zgsAoGPjdJoVAGDJEeYAADomzAEAdEyYAwDomDAHANAxYQ4AoGPCHABAx/4/ipt8+MgBG/oAAAAASUVORK5CYII=", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "plt.figure(figsize=(10,7))\n", + "for idx in range(res['all_dist'].shape[-1])[::-1]: #reverse this, so that the one with more distance is plotted on the backgroun\n", + " plt.hist(res['all_dist'][:,idx],bins=100,alpha=0.2, label=idx)\n", + "plt.legend()\n", + "\n", + "plt.figure(figsize=(10,7))\n", + "for idx in range(res['all_dist'].shape[-1])[::-1]: #reverse this, so that the one with more distance is plotted on the backgroun\n", + " plt.hist(res['all_dist'][:,idx],bins=100,alpha=0.2, label=idx)\n", + "\n", + "plt.xlim(0,500)\n", + "plt.legend()" + ] + }, + { + "cell_type": "code", + "execution_count": 16, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "4259246379.799023\n", + "17036985519.196093\n" + ] + } + ], + "source": [ + "D2_mean = res['all_dist'].mean(axis=-1)\n", + "print(D2_mean.sum())\n", + "out_d = np.zeros(img.shape)\n", + "out_d[_mask_d.astype(bool)] = D2_mean\n", + " \n", + "out_img = nb.Nifti1Image(out_d,affine=img.affine,header=img.header)\n", + "out_img.update_header()\n", + "out_img.to_filename(f'./all_mean_D2.nii.gz')\n", + "\n", + "\n", + "D2_sum = res['all_dist'].sum(axis=-1)\n", + "print(D2_sum.sum())\n", + "out_d = np.zeros(img.shape)\n", + "out_d[_mask_d.astype(bool)] = D2_sum\n", + " \n", + "out_img = nb.Nifti1Image(out_d,affine=img.affine,header=img.header)\n", + "out_img.update_header()\n", + "out_img.to_filename(f'./all_sum_D2.nii.gz')\n", + "\n", + "for idx in range(res['all_dist'].shape[-1]):\n", + " out_d = np.zeros(img.shape)\n", + " out_d[_mask_d.astype(bool)] = res['all_dist'][:,idx]\n", + " out_img = nb.Nifti1Image(out_d,affine=img.affine,header=img.header)\n", + " out_img.update_header()\n", + " out_img.to_filename(f'./{str(idx+1).zfill(3)}_D2.nii.gz')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "py3p9", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.9.7" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +}