diff --git a/examples/examples-sim-gLV/examples-sim-gLV.ipynb b/examples/examples-sim-gLV/examples-sim-gLV.ipynb index b412553..c59464a 100644 --- a/examples/examples-sim-gLV/examples-sim-gLV.ipynb +++ b/examples/examples-sim-gLV/examples-sim-gLV.ipynb @@ -2,7 +2,7 @@ "cells": [ { "cell_type": "code", - "execution_count": null, + "execution_count": 2, "id": "f07fa1f2-187e-4ce0-af95-31d6120977fe", "metadata": { "pycharm": { @@ -62,7 +62,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 3, "id": "3e0845a5", "metadata": { "collapsed": false, @@ -73,7 +73,38 @@ "name": "#%%\n" } }, - "outputs": [], + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "number of species: 5\n", + "specific growth rates: [1.27853844 0.55683415 2.06752757 0.86387608 0.70448068]\n", + "interaction matrix: \n", + "[[-0.05 0. -0.025 0. 0. ]\n", + " [ 0. -0.1 0. 0.05 0. ]\n", + " [ 0. 0. -0.15 0. 0. ]\n", + " [ 0. 0. 0. -0.01 0. ]\n", + " [ 0.02 0. 0. 0. -0.2 ]]\n", + "metabolite production: \n", + "None\n", + "perturbation matrix: \n", + "[]\n" + ] + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEGCAYAAABiq/5QAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAA+c0lEQVR4nO3dd3hUVfrA8e/JpPceQgmhFxEQIiqgFBGxgg3rihV111133V1X17qWXbuy/twVVhFsKIpYsCEgUlR6JwLSAiGk18n0Ob8/7oyZQIQEM5kk836e5zz3zmTmzjuTmfeee+455yqtNUIIIYJHSKADEEII0bIk8QshRJCRxC+EEEFGEr8QQgQZSfxCCBFkQgMdQGOkpqbq7OzsQIchhBBtyrp160q01mlH3u/XxK+U+hNwC6CBLcCNQDTwHpAN7AMma63Lj7Wd7Oxs1q5d689QhRCi3VFK7W/ofr819SilOgF/AHK01gMAE3AVcC+wWGvdC1jsuS2EEKKF+LuNPxSIUkqFYtT0DwETgdmev88GJvk5BiGEED78lvi11vnAs0AeUABUaq0XAhla6wLPYwqA9Iaer5SaqpRaq5RaW1xc7K8whRAi6PizqScJo3bfDegIxCilrmvs87XWM7TWOVrrnLS0o85NCCGEOEH+bOoZB+zVWhdrrR3Ah8BwoFAplQngWRb5MQYhhBBH8GfizwNOV0pFK6UUcDaQC3wCTPE8ZgrwsR9jEEIIcQS/defUWq9SSn0ArAecwAZgBhALzFVK3Yyxc7jCXzEIIYQ4ml/78WutHwYePuJuG0btXwghRAO01uw64GDJWjNXnhNPUpypWbffJkbuCiFEMDhQ6GDJ2loWrzFzsMiJKQQG9oxg+MDoZn0dSfxCCBFAhWVOlq6rZclaM7sOOFAKBvWMYPK4eM4cHEVCbPPW9kESvxBCtLjDpU6Wbajl2/W15O6zA9Cnazh3XJbI6KHRpCX6NzVL4hdCiBaQX+xg5SZLvWTfq0sYt0xMYNQp0XRKD2uxWCTxCyGEH9jsbjbusrFmm4XV260cLHICPsl+SDSd0lou2fuSxC+EEM3EYnWzeG0tKzbWsnGXDbtDEx6mGNw7gkmj4jj95Cg6pgY+7QY+AiGEaOPyCh188m01X/1gxmzVdE4P5cKRsZx2UiQDe0YQEd66rnkliV8IIU6Ay635YauFj5bWsO5HK6EmGD0kmkmj4+iXHY4xYUHrJIlfCCEayeXSbNltY8UmCys21lJU7iI10cRNFyVw/ohYkuObv+ulP0jiF0KIY7Da3azLtbJik4Xvt1ioMrsJC4WhfSP57eVJjBgYhcnUemv3DZHEL4QQDdh7yM78pTUsWm3GatfERilOPzmKkYOiObVfJFGRravdvikk8QshhIfLrflhi4UPl1azYYeN8DDFuFOjGZsTw8BeEYS2sZr9L5HEL4QIejUWN198V8NHS6spKHWRnmTi1kmJnD88xi9TJgSaJH4hRNAqqXAyb0k1n66oodaqOblnBFMvSWLkoLbXbt8UkviFEEHnQKGD9xZV8fUqMy4XjBoSzZXnxNM7KzzQobUISfxCiKCgtWbbHjvvL65ixSYLYaGK886I5YpxcQGbOiFQJPELIdq16lo3X68ys2BFDfsKHMRGKa4ZH88lY+LaTL/75iaJXwjR7mit2brbxoIVNXy7wYLdoenTNZw/X5vM2KHRbborZnPwW+JXSvUB3vO5qzvwEPCG5/5sYB8wWWtd7q84hBDBw+3WrNhk4Y3PK9mT7yA6UjHh9BguGBlLry7B0X7fGP682PoOYDCAUsoE5APzgXuBxVrrJ5VS93pu/81fcQgh2j+tNd9ttjDrs0p2H3TQJSOUv1ybzBip3TeopZp6zgZ2a633K6UmAqM9988GliKJXwhxArTWrNpq5fUFFew64KBTWij3TUlh7KnRmELab3fMX6ulEv9VwBzPeobWugBAa12glEpv6AlKqanAVICsrKwWCVII0XbszLMz7d0ycvfZyUwxcc9vkjlnWEy77n/fXPye+JVS4cDFwH1NeZ7WegYwAyAnJ0f7ITQhRBtksbmZtaCSeUuqSYwL4c/XJnPu6THtZjqFltASNf7zgPVa60LP7UKlVKantp8JFLVADEKIdmDNdgsvzCnjcKmLi0bGcuukRGKjpQ2/qVoi8V9NXTMPwCfAFOBJz/LjFohBCNGGlVe7+M8H5SxeU0tWRigv3p3OwJ6RgQ6rzfJr4ldKRQPnALf53P0kMFcpdTOQB1zhzxiEEG2X3aH5ZHk1b31RRa3VzZQLErh6fDzhYdKs82v4NfFrrWuBlCPuK8Xo5SOEEA1yODWff1fD219WUVLhYkifCH4/OZmumcE1tYK/yMhdIUSr4XJpFq428+bnlRwudTGgRwR/vyGFwb2lWac5SeIXQrQKKzfXMv3DCg4WOemTFc6frk4mp19kq75oeVsliV8IEVDFFU5eeq+cFZssZGeG8dhtqQwfGCUJ348k8QshAsLt1ny6vIb/fVyB0wVTJyVy+dlx0h+/BUjiF0K0uL2H7Dz3dhnb99oZ2jeSP16dFHRz4geSJH4hRIuxWN289WUVcxdVERsdwn1TUhg3LFqadVqYJH4hhN9prVmytpZXPqygtNLFuafHcPulie3yQuZtgSR+IYRf7Tpg56W55WzdbaN3VjiP3JrKSd0jAh1WUJPEL4Twi4pqFzM/reSzlTUkxITwl2uTmXBGDCEyXXLASeIXQjQrt1vz+Xdm/vdRBWarm8vGxHH9+QkymVorIolfCNFsdh+088Ico7fOoF4R3HVVMtkyzUKrI4lfCPGrWaxuZn1WybxvqomPDuHeKSmcI711Wi1J/EKIX2X5xlr+b245xRXGHPk3T0wgPkZ667RmkviFECfEbHHz7/fK+Hp1Ld07hfHQLdJbp62QxC+EaLLcvTYef72UwlInUy5I4LoJ8XKt2zZEEr8QotFcbs27C6t4fUElaYkmXrw7gwE9pJbf1kjiF0I0SlGZk3/NLmXTLhtjcqL501XJ0kWzjZLEL4Q4rqXrzLwwpxyHS/O365MZf1qM9NhpwyTxCyF+UWWNixffLefb9bX0zQ7n7zek0Dld+uW3df6+2Hoi8CowANDATcAO4D0gG9gHTNZal/szDiFE063YWMsLc8qornVz88UJXHWOnMBtL/zdQDcN+FJr3RcYBOQC9wKLtda9gMWe20KIVqK61s2/ZpXw0IwSUhJM/PdvHbh2QoIk/XbEbzV+pVQ8cBZwA4DW2g7YlVITgdGeh80GlgJ/81ccQojG+36LheffKaO82sX158dz7YQEwkIl4bc3/mzq6Q4UA68rpQYB64C7gAytdQGA1rpAKZXe0JOVUlOBqQBZWVl+DFMIcbjUycvvl7Nys3Hd2yfuSKN3VnigwxJ+4s/EHwoMAX6vtV6llJpGE5p1tNYzgBkAOTk52j8hChHcHE7NB0uqefPzSsC47u1lY+Oklt/O+TPxHwQOaq1XeW5/gJH4C5VSmZ7afiZQ5McYhBC/YONOK9PeLWP/YScjB0XxuyuSyEiWjn7BwG//Za31YaXUAaVUH631DuBsYLunTAGe9Cw/9lcMQoij2exuXppbzuffmemQYuKJO9I44+SoQIclWpC/d++/B95WSoUDe4AbMXoSzVVK3QzkAVf4OQYhhEdxuZOHppewI8/O1ePj+c358USGy+jbYOPXxK+13gjkNPCns/35ukKIo23dbePh/xVjs2seuz2VEQOjAx2SCBBp0BMiCCxYUcO/3ysjIzmU5+5Kk6tiBTlJ/EK0Y06X5uX3y/l4WQ2n9o/kgZtSiZOJ1YKeJH4h2qmSCiePzyxl8082rhwXxy2TEjGFSDdNIYlfiHZpzXYL/5pVitWuuf/GFM4+NSbQIYlWRBK/EO2Iy6WZtaCSdxZWkZ0ZxkM3p9JV2vPFESTxC9FOFJc7efz1Urb8ZOP8ETHceUWSdNUUDZLEL0Q7sGqb0bTjcErTjjg+SfxCtGEOp+a1TyqYu6ia7p3CeOiWVLIypGlHHJskfiHaqPwiB4/PLGVHnp2Lz4zljssSiZCmHdEIkviFaIMWrTbzwpwyTCHwyK2pnHWKjMIVjSeJX4g2xGJ1M+29chauMnNyjwj+fmOKzKgpmky+MUK0ET8dsPPoayUcKnZy/fnx/OY8uRyiODGS+IVoAxb+UMPzc8qJiw7hubvSGdQ7MtAhiTZMEr8QrZjDqfnPB8ZcO4N7RfDAzakkx5sCHZZo4yTxC9FKFVc4+cf/Sti+184VZ8cxdVKiNO2IZiGJX4hWaNNOK4/OLMFi0zx0Syqjh0ivHdF8JPEL0YporZm/tIb/zCunY6rMnS/8QxK/EK2Ew6mZ9m4Zn39nZsTAKO6dkkJMlAzIEs3Pr4lfKbUPqAZcgFNrnaOUSgbeA7KBfcBkrXW5P+MQorUrr3bxyIwStuy2cd158dxwQQIhMne+8JOWqE6M0VoP1lp7r717L7BYa90LWOy5LUTQ2n3Qzm+fOsyOPDsP3JTCTRclStIXfhWI48iJwGzP+mxgUgBiEKJVWL6xlt8/V4jTBdPuTmdsjsyqKfzP3238GliolNLAdK31DCBDa10AoLUuUEqlN/REpdRUYCpAVlaWn8MUomVprXn7yypmflpJ3+xwHrstjZQE6Z8vWoa/E/8IrfUhT3L/Win1Y2Of6NlJzADIycnR/gpQiJZmd2iefauURWtqGXdqNH+5LoXwMGnaES3Hr4lfa33IsyxSSs0HhgGFSqlMT20/EyjyZwxCtCbl1S4eml7Mtj12broogWsnxKOUJH3RsvzWxq+UilFKxXnXgfHAVuATYIrnYVOAj/0VgxCtyd5Ddn739GF2HXDw8C2pXHdegiR9ERDHrPErpe5uxDbMWuvpDdyfAcz3fLFDgXe01l8qpdYAc5VSNwN5wBVNjFmINmf1NguPvlZCZLjixT+l0zc7ItAhiSB2vKaevwL/BY5VLbkdOCrxa633AIMauL8UOLsJMQrRZv08EveDcrp3CuPx29NIl/nzRYAd7xv4ptb60WM9wNOMI4Q4QmGZkxfeKWP1divDB0Zx/w0pREXKSFwReMdM/Frre463gcY8Rohg4nZrPl1Rw4z5FWjgziuSmDQqVgZliVajUcecSqm7gNcxpl94FTgFuFdrvdCPsQnR5hwscvDsW2Vs/snG0L6R/PnaZDqkSNOOaF0a+428SWs9TSl1LpAG3IixI5DELwTgcms+WFzN6wsqCQuFv16XzIQzYqTXjmiVGpv4vd/e84HXtdablHyjhQDAanfzxOulrNxkYcTAKO66KonURKnli9arsd/OdUqphUA34D5P/3y3/8ISom0oq3LxwCvF7Nhv53eXJ3LpmDip5YtWr7GJ/2ZgMLBHa12rlErBaO4RImjlHXZw38tFlFW5eXRqKiMGyVWyRNvQ2L5lGugP/MFzOwaI9EtEQrQBm3ZZ+f2zhVjtmuf/lC5JX7QpjU38/wHOAK723K4GXvZLREK0covXmLnnpSKS4kL4v792oJ+MwhVtTGObek7TWg9RSm0A0FqXK6XC/RiXEK2O1po3v6hi1oJKBvWK4NHb0oiLlgFZou1pbOJ3KKVMGE0+KKXSkJO7IohY7W6eebOMb9bVcs6waP58rUylLNquxib+fwPzgXSl1BPA5cADfotKiFakuMLJg6+UsOuAnamTErnyHOm5I5qH0+1k0+FNbCnaQoW1gkprpbG0VRrFWslT455iaMehzfq6jUr8Wuu3lVLrMCZXU8AkrXVus0YiRCv04z4bD04vodbq5rHbUhk+UE7iihNXYa3g+wPf892B71h5YCWr8ldR66it95jY8FgSIhJIiEwgMTIRu8ve7HEcb1rmeK11lVIqGeOCKXN8/pastS5r9oiEaCUWrzHzzFtlJMeH8H9/zaBbRzmtJRrP7rKzuXAzq/NXszp/NWsOrSG3OBeNJkSFMLjDYG4afBMjskYwNHMoKdEpxEfEExri/8F/x3uFd4ALgXV42vc9lOd2dz/FJUTAuN2aWZ9V8tYXVQzsGcEjt6aSGCfXwxXHVlhTyPK85Szfv5wf8n9g4+GNP9fW06LTGNZpGFeedCUjuozgtM6nERseG7BYjzc754WeZbeWCUeIwLLY3Dw5u5TlGy2cPzyGu65KJixU2vPF0faW72XZ/mUsz1vOsv3L2FW2C4DosGhO7Xgqd512F8M6DePUjqeSlZDVqs4LNXZ2zkuAJVrrSs/tRGC01voj/4UmRMsqKnPywPRi9hx08NvLE7lMpl8QPiqtlSzZu4Sv93zNwt0L2V2+G4CkyCTO7HomU4dO5ayuZ3FKh1MIM4UFONpja2xj0sNa6/neG1rrCqXUw8BHfolKiBaWu9fGg9OLsdo1T/w2jdNOigp0SCKA7C47P5X9xI8lP7Lp8CYW7V3EqoOrcGkXseGxjMkew12n3cWYbmPon9afENW2xnM0NvE39K4ae7RgAtYC+VrrCz0nit8DsoF9wGStdXkj4xCi2S1eY+bpN0tJTTDxzB/S5SRukLG77CzZu4Qle5fwY8mP/FjyI3vK9+DSLgAUipyOOdw78l7G9xjP6Z1PJ9zUtr8jjU38a5VSz2NM06CB32Oc8G2Mu4BcIN5z+15gsdb6SaXUvZ7bf2t8yEI0D5dbM2tBJW9/aZzE/cfUVBJi5SRuMLA6rSzcvZB5ufP4ZMcnVFgrCDeF0zulN4M6DOLKk66kb2pf+qb2pXdKb+Ii4gIdcrNqbOL/PfAgRk0djAuwHHcAl1KqM3AB8ARwt+fuicBoz/psYCmS+EULq6xx8c9ZpazZbuX8ETHcdaWcxG3vrE4rn+/6nPe3v8+CnQuosdeQFJnEpL6TuKzfZZzT/RwiQoNj3qXGDuAyA/cqpWK11jVN2P6LwD2A7+4yQ2td4NlugVIqvaEnKqWmAlMBsrKymvCSQhzbzjw7D88opqzKxd3XJHPBCLlSVnvlcDlYvHcxc7bOYX7ufKrt1aRFp3H1gKu5vP/ljMke0+pPxPpDY9vph2NcazcWyFJKDQJu01r/9hjPuRAo0lqvU0qNbmpgWusZwAyAnJwcfZyHC9Eon6+sYdp7ZSTFmZh2dwZ9ZWbNdsfpdrIybyXvbXuP97e/T0ltCQkRCVze/3KuHnA1Y7qNaZFBUq1ZY9/9C8C5wCcAnksvnnWc54wALlZKnY8xd3+8UuotoFAplemp7WdijAgWwq/sDs2/55bx+UozQ/tG8sBNKdKe345U26r5avdXfLLjEz7b9RllljKiQqO4uM/FXD3gaib0nBA0zTiN0ejdntb6wBGHw67jPP4+4D4AT43/L1rr65RSzwBTgCc9y4+bFrIQTXO41Mkj/ythZ56dayfEc8OFCZhCpGmnrTtcc5iPfvyIj3d8zJK9S7C77CRHJXNBrwu4uM/FTOg5IaCjY1uzxib+A57mHu2Zh/8PGD11TsSTwFyl1M1AHnDFCW5HiONam2vh8ZmluFyax26TyyO2dflV+XyY+yEf5H7A8v3L0Wh6JPXgzlPvZGLfiQzvMjzom3Eao7Gf0O3ANKATkA98BfyusS+itV6K0XsHrXUpxiyfQviN262Zs7CKmZ9Wkp0Zxj+mptI5PfhO4rUHeZV5RrLf/gErD6wE4KS0k3h41MNc1v8yTko7SU7ON1Fje/WUANf6ORYhmkWNxc1Ts0tZudnC2Jxo/nxtMlERbWtkZbDbUbKDD3M/5MMfP2TtobUADMoYxGNjHuOyfpfRL61fgCNs2xrbq6c7Ro3/dIwBXN8Df9Ja7/FjbEI02d5Ddh6aXsLhUid3XpHEJaNjpTbYRmwv3s67W99lXu48thdvB+C0Tqfx1LinuKTvJfRK6RXgCNuPxjb1vIMxavcSz+2rMObmP80fQQlxIr5dX8tTb5QSHal4/o/pnNwzMtAhieMot5Tz7tZ3mbVpFqvzVxOiQhjVdRR35NzBpL6T6BzfOdAhtkuNTfxKa/2mz+23lFJ3+iMgIZpKa81bX1Tx+oJK+ncL5x9T00hJkK6arZXL7WLRnkW8vvF1PvrxI2wuGyenn8zz45/n2oHXkh7T4JhO0Ywam/i/8cyr8y5GU8+VwGeeCdeQK3GJQLHZ3Tz9VhnfrJWLoLd2B6sOMnPDTF5d/yoHqg6QHJXMrUNu5cZTbuSUDqdIk1wLamziv9KzvO2I+29CrsQlAqS00sUDrxSzM8/OLRMTuHp8vCSPVsbldvHV7q+Yvm46C3YuwK3djO8xnufGP8fFfS6WQVUB0thePXIFLtGq7Dpg5/7/FlNjcfOPW1MZOVj657cm+VX5Ru1+w6vkVeaRHpPOPcPv4daht9I9SeqJgdbYXj1XAF9qrauVUg8AQ4DHtNYb/BqdEA1YtqGWJ2eXEh8Twkt/zqBH57Y9N3p74XK7+PKnL5mxfgaf7fwMl3Yxrvu4n2v3bX0O+/aksU09D2qt31dKjcSYs+dZ4BWkV49oQVpr3vmqitc+MU7iPjo1jWQ5iRtwByoPMHPDTF7b8BoHqg6QEZPBX4f/lVuG3EKP5B6BDk80oLGJ3zsvzwXAf7XWHyulHvFPSEIcze7QPPd2KV+vruXsU6P563VyEjdQyi3lfLv/W77Z+w3f7PuGLUVbUCjG9xjPixNe5KLeFwXlVMdtSWMTf75SajowDnhKKRVBw5djFKLZlVe7eGh6Mdv22LnxogSumyAncVvartJdzFg3gyX7lrChYAMaTVRoFCOyRnD1gKu5asBVdEuSU4FtRWMT/2RgAvCs50LrmcBf/ReWEIa9h4yTuOVVbh66JZXRQ+Qkbks6XHOYR799lBnrZmAKMTG8y3AeGf0IY7LHMKzTMOmV00YdM/ErpdZrrYdorWuBD733e66gVeD7GP+GKYLRD1ssPP56CVERIbx4dzp9ukqSaSlVtiqe/e5Znvv+OewuO7fn3M6DZz1IRmxGoEMTzeB4Nf5+SqnNx/i7AhKaMR4hcLo0Mz+t5N2FVfTsEsYTt6eRliRT7bYEm9PG9HXTeWzZY5TUlnDlSVfy+NjH6ZncM9ChiWZ0vF9T30Zs45gXZBGiKQrLnDz2Wgnb99q5aGQsv708kYhwOZ3kb1sKtzBzw0ze3PwmpZZSxnYby1PjniKnY06gQxN+cMzEr7Xe31KBCLFyUy1Pv1mGy6158KYUxuTEBDqkdq3CWsGcLXOYuXEmaw+tJSwkjEl9J3Hb0NsY222snEBvx+T4WQScw6mZPr+CD7+ppndWOA/enEKnNOkO6C+7y3bzxPInmLN1DlanlYEZA5k2YRrXnHwNqdGpgQ5PtABJ/CKgDhQ6eOL1Unbm2bl0TBxTJyVK/3w/2Vexj8eXPc6sjbMIM4Vxw6AbuGXILQzJHCK1+yDjt8SvlIoElgERntf5QGv9sGdGz/eAbGAfMFlrXe6vOETrpLXms5Vm/vNBOeFhikenynw7/pJXmccTy55g5saZmJSJ3536O+4deS+ZcZmBDk0EiD9r/DZgrNa6RikVBqxQSn0BXAos1lo/6Znq+V7gb36MQ7QylTUunn2rjJWbLQzpE8HfpqSQligHn83tUPUhnlj2BP9b/z8Abht6G/eNvI9O8Z0CHJkINL/92rTWGqjx3AzzFA1MBEZ77p+NcRF2SfxBYs12C0+9UUp1rZvbL03k8rFxhIRIM0NzKjYX89TKp3h5zcs43U5uGnwT9591P1kJWYEOTbQSfq1mKaVMwDqgJ/Cy1nqVUirDMwAMrXWBUqrBy+0opaYCUwGysuQL29bZHZr/fVzBvCXVdO0QypO/S6dnF5mtsTlVWCt47rvneHHVi9Q6avnNwN/w0KiHZBpkcRS/Jn6ttQsYrJRKBOYrpQY04bkzgBkAOTk52j8RipaQV+jg8ddK+Omgg0mjYrntEumb35wKqguYuWEmz37/LBXWCiafNJlHRj1Cv7R+gQ5NtFIt0rDqmd9nKcZ8P4VKqUxPbT8TKGqJGETL01qzcJWZae+VEx6qePz2VIYPlBO4zaGwppB5ufOYu20uy/YvQ6O5sPeFPDbmMQZ3GBzo8EQr589ePWmAw5P0o/DM7Al8AkwBnvQsP/ZXDCJwaq1uXpxTxqI1tQzqFcHfb5QTuL9WpbWSd7e+y9ztc1m6bylu7aZfaj8eHvUwk0+aLDV80Wj+/CVmArM97fwhwFyt9QKl1PfAXKXUzUAecIUfYxABsGO/jcdmlnK4xMkNFyZw7YR4THIC94QVmYuY9sM0Xl7zMpW2Snqn9Ob+M+9n8kmTGZDe6NZTIX7mz149m4FTGri/FDjbX68rAsfl1sxdVM3rn1aQFGfi+T+lM7BnZKDDarPyKvN49rtneXX9q1idVi7rfxn3DL+HnI45MuBK/Cpy7C2aRX6RgyffKGXbHjtnDo7iz9cmEx8jl0U8EbtKd/GvFf/izc1vAvCbgb/hnhH30De1MXMmCnF8kvjFr6K15pNlNUyfX0GoCf5+QwpnnxotNdITkFeZx6PfPsqsjbMIN4VzR84d/GX4X6T/vWh2kvjFCSsud/LMW2WszbWS0y+Sv16XLPPmn4DDNYf55/J/Mn3ddADuHHYn9428Ty56IvxGfqWiybTWLFpdy7/nluFywR+vSuKiM2Ollt9EZZYynln5DP9e/W9sThs3Dr6RB0c9KDV84XeS+EWTVFS7ePHdMpZtsDCgRwR/uz5ZplBuIqvTykurXuKJ5U9QZavi6pOv5pFRj9ArpVegQxNBQhK/aLSVm2t5/u0yaixupk5K5IpxcdJNswnc2s3bm9/mgW8eIK8yj/N7nc+/zv4XAzMGBjo0EWQk8YvjMlvcvPxBOV9+b6ZH5zCe+UM63TvJPDtNsWjPIu75+h42HN7AkMwhvD7xdcZ2GxvosESQksQvjmnDDitPv1lKcbmLa8+N5/oLEggLlVp+Y9hddj7b+Rn/Xftfvt7zNV0TuvL2pW9z1YCrCFEyV5EIHEn8okEWm5v/fVTBR9/W0Dk9lGl/zuCk7hGBDqvV01qzvmA9szbOYs7WOZRaSukQ24FnznmGO4fdSWSoDGgTgSeJXxxl0y4rT79ZRkGJk0tHx3LzxESiIqSGeiwV1gpeW/8aszbNYmvRViJMEUzqO4kpg6ZwTo9zCA2Rn5poPeTbKH5mtbt59eNK5i+tpkNKKC/8MZ1BvaWGeixlljJe/OFFpq2aRpWtijM6n8ErF7zC5JMmkxSVFOjwhGiQJH4BwNbdNp56o5T8YicTR8UydWIiUZFSy/8lJbUlPP/987y0+iVq7DVc3v9yHjjzAQZ1GBTo0IQ4Lkn8Qc5idfPaJxXM/7aGjGQTz92Vzil9pJb/Sw7XHOb575/nP2v+Q62jlisHXMn9Z94vs2SKNkUSfxBbm2vh+XfKOFzqYuKoWG6dmEi01PIbtKNkB899/xyzN83G6XZyzcnX8PeRf5c58EWbJIk/CFXXuvnvPKNffuf0UF68W6ZP/iU/HPyBp1c+zUc/fkS4KZybBt/E3WfcLaNsRZsmiT/IrNhYy4vvllFR4+bq8fFMuSCB8DDpl+/L4XKwYOcCXvjhBZbnLScpMon7z7yfO4fdKROniXZBEn+Q2F/g4L8flrN6m5WencP452/T6Z0lo299/VjyIzM3zGT2ptkUmYvoEt+FF859gVuG3EJseGygwxOi2Ujib+cqa1y88XklHy+rISpCcfuliVw6Jo5Qk9TyAWrsNby/7X1e2/AaKw+sxKRMXNj7Qm4+5WbO63We9L8X7ZJ8q9spp8u4QMrszyoxW9xcODKWGy5MIDFOrooFYHFYeGn1Szy54knKreX0TunNU+Oe4vpB19MhtkOgwxPCr/yW+JVSXYA3gA6AG5ihtZ6mlEoG3gOygX3AZK11ub/iCEZrtlv4v/fLOVDoZGjfSH57eSLdOkqzDoDT7WTWxlk8svQR8qvzOa/nefz9zL8zossIuZ6ACBr+rPE7gT9rrdcrpeKAdUqpr4EbgMVa6yeVUvcC9wJ/82McQaO43MnLH5SzbIOFzumhPHFHGqcPiJSEhjGHzoe5H3L/kvvZUbqD0zufztuXvs2o7FGBDk2IFue3xK+1LgAKPOvVSqlcoBMwERjtedhsYCmS+H8Vp0szb0k1sz+vxO2GGy9K4Mpx8dJbB6iyVfH+tvd5Zd0rrD20ln6p/Zh/5Xwm9pkoO0QRtFqkjV8plQ2cAqwCMjw7BbTWBUqp9F94zlRgKkBWllyK7pds2mVl2rvl7CtwcPqASH4/OZnM1OA+deNyu1i8dzGzN81mfu58LE4LfVL68NrFr3H9oOvlhK0Ien7/BSilYoF5wB+11lWNrWVprWcAMwBycnK0/yJsm/IKHcxeUMk362rJSDbx2O2pjBgYHeiwAkZrzZaiLbyz5R3e2vwW+dX5JEYmcsPgG5gyaArDOg2TGr4QHn5N/EqpMIyk/7bW+kPP3YVKqUxPbT8TKPJnDO3NoRInb35eyderzISHKa47L55rzo0nMjz4plrQWrO1aCtzt81l7va57CzdiUmZmNBzAi+c+wIX9blI5r8PdlqDzQZWK1gsxlIpiI2FmBiIjDRun8h2q6uNbbpcRnE669aVgtDQo4vDAWVlUF5uLL2lshJqa41isdRf/vOfMHRos34s/uzVo4DXgFyt9fM+f/oEmAI86Vl+7K8Y2pOiMidvfVnFF9/VYDIpLhsbx1Xj40kKwu6Zu0p38ebmN5m7bS47SncQokIYnT2au0+/m0v6XUJ6TIOth+JEaW0kNYCwsKY9126HoiIoLKxfLBYwmepKaKixVMpIjkcWp9P4W0hIXfEm7KqqoxNpeTnU1BiJ/lhCQoydQGwsxMVBcnJdSUoyljExxnvIz4dDh4xlfr6RlJtTRARER0NUVP2l3d68rwMorf3TiqKUGgksB7ZgdOcE+DtGO/9cIAvIA67QWpcda1s5OTl67dq1fomztbNY3by+oJKPl1WjNVwwMpZrz40nNTG42qmdbief7vj058sYepP9Ff2v4NJ+l7bPZO92GwmspASKi41SXQ3h4XUlIsJYhoRAaamRVH0TbVGRkfx8a6RHLhsqTqdRHA4jDq+oKEhMNEpSkrGMizNeo7raSLbV1XWlsvLXfQYhIcbOJjTU2AFpbcTjLVpDfPzRyTopyYgrMtKI2buMijKeZzYbsfqWhnYgVVVGHOHh0LEjdOpUVzp2NBKz747LW6DuM/QtoaFHx5mcDAkJdc9rRkqpdVrrnCPv92evnhXALx1Dne2v121P1my38Nw7ZRSXuzjvjBiuOy+BDinBlfALqgt4df2rTF83nfzqfDrHd+axMY9x8yk3kxmX2fIBaW3U9MzmulJTYywtll+ulR4+DPv3Q15e3TIvz0iYYWF1JTTUWNrtRiL3TbpNkZQEGRmQlmase5OTb4I6Mln5Fm8s3ni8ibeyEioqjFJebuxcdu0yEmpsrPFaWVlG0o2Lg5QUI44jS3R0wzscret/HiEBbsJ0Oo3/bXz8iTUJtVLBlUXaiMoaF/+dV8HCVWayMkKZdncGA3oEz/VurU4rn+/6nLc2v8WnOz/F6XYyvsd4Xj7/ZS7ofUHTe+U4HMaP126vX5P1Ls1mo2ZXWWksvetlZXU17eLiupr3rzn0TkszEmO/fnDuuUYC9G3O8K6HhhqPTUuD1NS6ZXy88Ti73Wi7ttvr3pc3yaalGTXU1s67Y2nNQkON2ng708o/9eCitebbDRZeeq+MKrOb35wXz7UTgmP2TJfbxdJ9S3lnyzvMy51Hpa2SjJgM7jrtLm4belv9aZC1Nmqc+/cfXfLz65ocvMVmO7Gg4uPrkm9WFgwZYqwnJ9edHPQtUVFGrdC3KcLbHJGRAV26GIleiACTxN8KuNyaTTttzPummu+3WOiTFc7Tv0+mR+c2UGv7FSwOC8v2L+OLn77g/e3vc6j6EHGhMVyadhbXRJ/GWFtHQtcWw4KXoKCgfrFY6m8sMhK6doXOnY22V98Tdt4kHR5evznFW+OMjTWSvLckJBjP80ObqxCtgST+ANpf4GDhKjOLVpsprnARE6W47ZJELh8bh6m9zJ5ptxu188pKdGkpubt/4KuDS/mqeiPfmg5gDXET4YQJu+GazXDRDjNRzi+AL+q2kZAAmZlGOe00Y9mxo5HovSUtrV21wQrhT5L4W5jZ4mbhKjMLV5nZsd9OSAgM6x/JHZclMnxgdOtv1nE4YONGWL4cdu6s63tsNtet19TUnQS0WLCEwvQcePF02J9obKZvleK2kkTOdWYzKv5kort1hlNTj27TTk83mlCEEM1GEn8LKa10Me+baj5dVo3ZqunZOYw7Lkvk7FNjSI5vxU0KZjOsWmUk+uXL4YcfjPvASMzeZpToaKOkpUF2NiQmUpsYw/S4nTzNSg7rKkbFDuD+HpMYP+gyumYPkhq6EAEiid/P8godzP26iq9Xm3G54KxTopk8Lo6+2a2wl05FhVGbX7++ruzYYZygVAoGDoQbb4QzzzRKZsPdKWsdtbyy9hWeXvk0heZCxnYby3ujHuasrme16NsRQjRMEr+f5O6z8c5XVXy32UJYqOL84bFcfnYcndKaOPLx16qtrT/asKDA6B/uO0ilrMzoqnjgQN3zOnUyerFMngzDhsHw4cZgHQ+tNYU1h9lfsZ+8yryfy/7K/Xx34DuKa4sZ130cH4z6gJFZI1v2PQshjkkSfzPbvtfGG59Vsnq7lfiYEH5zXjyTRsX558pXZrMxeMab1A8erL/Mzzdq8Ufyjh70lk6d4OSToW9fI9mfcorRtu5Da82Pxbks3beUb/Z9w9J9SymuLa73mLjwOLomdmVU9ij+eNofGZE1ovnfsxCtnHcyhNbckimJv5ls22Pjjc8rWeNJ+LdOSmTiWbFERzbjyMPiYli50mhrX7HCaIrxzqECxjetQwejS2OvXjB69NFDzDt2bNQoxGJzMduKt7GlcAvfHfyOpfuWcrjmMACd4jpxbs9zOa3TaWQnZpOVkEXXhK4kRB490MU7yl7rukGtjWWzGQcsbrcxqNN3qVT9kfj+HOCpdf1pY7yzE1RXG2O9fJcNlZoaI07vYNa4OONf4O0xeuQsCk6n0dvUOzNCQkLdemxswzMX+A7w9Z2FxTvNjnecl28JDa3r7eqNKzbWiKmhGQ1qa+vi9P4vvINtIyPrTvUcuYyONt7PkV85h8M4APUekB46ZPQJiIysK97/cUSE8Xxvgbql1Vp/ILW3n0FkpDGmLSXFOB3lXQ8Pr5tZwvf9eQdRe2fD8K57Yz/ytbU2Yt671yj79tUtoW4ISHp63XpiYv2Byd7X8Mbq7dOQmurfIR9+m6unObXWuXq01mzcaWPOwirW5lpJiA3hynHxTDwrlqgTTfhOp1Fj936DvN+m1auN9nYwfgXDhsHIkUbtvEsXI7F36ND0SbSAwppCvvjpCzYUbGBzQS5bdpdSWhAFlVlQ1YW4sER6d+zIwC7Z5HTrTa/MDiQkKOz2ui/9nj11y4ICIyE0JDa27geYnGwsk5KMbvm+g2NLSowfZmOFhxtJIiKiLin6JketjUTkTbreEhdnfOTeAbu+xWIxEuQvvZdj8Q4NiIszXtdmq7+zOJFttmUmU92cY9HRRmIuLq6/k2pOoaH160T+lpgI3boZJTvbqIgUFdUf+F1UdPTwk2OJijJ2ALNnw5gxJxZXi8/V054VVzj56nszX/5g5lCxk8TYEKZeksjEM5uQ8LU25mrZuLGubNpk3OebFZQykvqgQcaJ1ZEjISfHyHA+myoqgv0b6qaCOXTImErFO6WKd1ldbfwIQ0IdODBj0ZVY3JUQMhBVew66ugPo+s1S1cA6T3n9F95Op07QvTuMHWush4bWryV5a0gVFcYphdJSo+TlGbejoupqO7161fXojInxxBtSf+l2Gz8i70y73nWbrf40Od51MJKNb2Lfu9dYhoXV7Qi6datbj4qqX+vzloiIuqTuW3v3lpiYYx+BaF131OB21582x7t0OOqmxPEtNTX1pwLyfa++NWrfdd9arG9t1uk8el61mhrj/iPHv8XGGp+H7/Q+3v+FUnW1bt+evb49fH1nGTabjRqu74Godz0x8ehZlL3Fe6Tj/Qy9696jDd8SFlY3A3JpqVGR8C6dzvrvy7seEVF3VOc9KvKu/9JrZ2b+3ImtUbyzcni37V23WOri81Z8vOWIVtdmITX+RnI4Nd9vsfDFdzWs2W7FrWFwrwgmDI/lrFOijj0fvttttMWvW2eU9euNRO9pf9coynueSl63URSm9KcoOpvi8I4UudMossVTVGLCYmm4yaOiom6uL19RUXWT/yUmasJja3GFl1BNPnvK91NRYwZXOInhGXSIyiI9siPdOsWT1UWRlWUcRHiXYWFHN2tUVRk//O7djfFTkTLtvRCtzi/V+CXxH4PWmh/32fl6tZlv1tVSWeMmJcHEhNNjmHBGDJ3SG2hW8SR595p1HFr2ExXr91Cx/RCVljAqSKQyNJXSDv05ENefPNWVvNpU8oqjMJuPbnMPD69rH/TWIr01Le96bGz9AaydOruojtrKlqrlbCnazObCzWwt2orZYfS9NykTo7JHManPJC7uczFdE7v6+2MUQgSINPU0QX6Rg0Vralm02kx+sZOwUBg+MJpzT4vh1P6RddMpHFGT37fiIF9vSmeRbSSLGU8p1xy9cSdw0EjoWVnQry+cm2WsZ2UZzfTp6UZpzEywWmtyS3JZsncJs/YuYenXSym3lgOQHJXMoIxB3DLkFgZmDGRgxkD6p/UnOkwmChMimEni99Bas2qrlbe/qmTbHjtKweDeEVwzIZ4zB0cTG+rCvimXopdyKV2zh5LNhyjdVUapLYb1DGGR+i27dQ8AOiaauWCkjTMmuEhJN/3cK8O3d0ZEI8dvaa0pNBeSV5nHoepDR5WNhzdSaC4EoFtiNy7tdyljssdwVtez6BzfWa4zK4Q4StAnfq01G3bYmPlpBdv32slMDeXmC+LpWFHG/m/yWPi6nWf2xrOlqislDAQGHrWN2GgXY8Yo/jAezjkH+vaNQamYJsVhtpspqClgR8kOcktyyS3ONZYluVRYK+o91qRMdIjtQMe4jozrPo6x3cYyJnsM3ZK6/YpPQggRLII68W/dbeOV90rYftBFpMtJwg4zm9fGMNechpNEoDtR1DIgbj+TBuyma78CUvpnkNI3jdT0kJ+7JWZkmBrsRenWbkpqSyioLqCgpqDe8rD5MIdrDlNYU8jhmsNU2+v3XUyPSadfaj+uOukq+qX1o3tSdzrGdaRjXEfSotMwhbTi+X2EEK2aPy+2PhO4ECjSWg/w3JcMvAdkA/uAyVrrcn/FcBSzmTXvbOGdL2CzTkFnRmCrDWH/+iQKtsfS0Z3PwIQ8Lhr6E4OGRTDo/E70Gt8NU3i/epuxOW3sq9jHnvI9LD+0m0M7DlFaW0qJpYSS2hJjvbaEUkspTvfRnYkTIhLIjMukQ2wHhnYcSoeYDnSINUrP5J70S+tHclRyS30qQogg48+LrZ8F1ABv+CT+p4EyrfWTSql7gSSt9d+Ot60T7dWTf2AbO5YsZdE3UWwo6U1VfGciEo1ul06rxpVfQGLYVtL7FpE8sJrQTBcOtxOnT3G5XTi1E4fLwaHqQ+wp38PBqoNo6j43kzKRGp1KSnSKsYwylqnRqWTGZtIxriOZcZlkxhrJPipMphkWQvhfIC62vkwplX3E3ROB0Z712cBS4LiJ/0Q98vhSklLi2RAzHFtoODWV1RwoW0N5xlvUpC+BrnXj3EP2hRCWF0aYKYywkDBMISZCQ0IJDQnFpIz1DrEdGJ09mh5JPeie1J3uSd3pkdyDjJgMOYkqhGgzWrqNP0NrXQCgtS5QSv3imDSl1FRgKkBWVtYJvdjN/b/k9LQFmF3xlKechT3rXFwZY4iOGEdkaCThpvB6iV4IIYKBXwdweWr8C3yaeiq01ok+fy/XWicdbzsnPIDLZYWChZD3AeR/DI4qCEuAzhMhazJkjoeQFp4mWQghWkhrGcBVqJTK9NT2M4Eiv76aKRI6X2wUlw0OL4K89+Hgx7D3DYhIha5XQbfrITmndc+jKoQQzcSPk9k26BNgimd9CvBxi72yKQI6XQBnzIJLC+GsTyBjDPz0P/hqGHzWD7Y+Aeb9LRaSEEIEgj979czBOJGbChQCDwMfAXOBLCAPuEJrXXa8bfl1rh57hdEUtPcNKF5u3BeWCDFZEN2lbhmdBYknQcJJ0jwkhGgTZJK2xqjZCwc/gurdUJsHtQfAnAd2n31TSAQkDTKahpJzICUH4vtBSFCPhRNCtEKtpY2/dYvtBn3/dPT9TrPRBFS+GcrWGmXvm7DrP8bfQ8Ihrhck9Dd2AvH9jPW4HmCKlnMHQohWRRJ/Y4TGGIk8oT9kX2Xcp91QvQtK10DFFqjKhbINcGCe8TcvFQrhCUZvorAECIuHiBRIHAhJQyB5CER1lJ2DEKLFSOI/USoE4vsYxZfLClU7oXI7mPcZXUgdlUaxV4KzCiq3wYH54B39G5lRtxOIyYKIdIhMN+6PTIfQWNkxCCGajST+5maKhKSBRjkWRw1UbIKy9VC2DsrXw/aFoBu4GKspEiI7QHQniPKU6I7GMizO2Nn4FrfVOOoIT4LwZOMIw3dpipQdiRBBTBJ/oITFQtoIo3i57GArAqtPsRWBtRAsBVCbD+UbIH8BuGpP/LVDwuqanbxNUOEJxhFGVGdjBxPd2bOD6Qyh0cYOyVvczrodVEiYUZR36dmhaA1uG7gsPjsku3FyPDTGOIqRHZAQASGJvzUxhRuJNrrzsR+ntdF0VJtv7ABMUUYSNUVCiGcJYC83eiTZSsFeCrYyY+moMpqdfm6GqoKaPVD8HdiKf917UKFGM5jb3pgHe3YCMRCeCOEp9Y9OIlKMHYrL6tmJ+BzRoOqaw3ybxsKTjMc6a43Pxrt0WYwdq3aA26dop7EDjMo0SmQmRKYZ7wGMIydrMVgLoPYQWA4Zn6tSRgwqxGeJz7bt9Yv3/aJ8nquM9+OoNv4Hzmpj3VltvK53B/nzMvbo+8JiwRRjLEPjPDvzuLp1U6TxGrZio1h9ltpl7NRNUUbxrqPqN1F6vy8us8978HnfymR8dtFZRlNlTJbxOXqnQXFaoPag0Uuu9oCxrt0Q1cE4kvUuIzOM3nG2MqOyYz0MlsPG0lbq+QhNnhJSt3Q76n83vOvK5PldRNVfhkTUVVh+rrSEGt9d45/uWXiWKsT4LMMTja7e4UnGuinCU1krPqKiVmxsLzKjfglPqvuf20p9fpclxv/95/d0REkdDlEZTf8tHoMk/rZIKU+iTDz248JiIaZL07btshnJrfagsWOxHKz7ESmTJ7F71uHoROq2A+66HVC9H1y4sS1njdFTyltcZiOZ2kqNxFC+0Vj3PapRIT7biTSSlq244aaxX0uZjB8qykhA+uiptRu9nZBwz2emPInEU7zrpkhPkvZJ1tGdjOd6Px9rITj3GJ+bo8ZYNjYmFVK/s8GJ8u5owLM9Xbd02404j3zvUR2N/6E3aTcqXlPD/1NvUtYuoIEu6CFhnu9cRF1y125jR+C01O0QmlNIeCMrOD4xqrCmH62P/gKiJjTtOcchiV/UZ4owurXGtoKreTktxg/dFNnwOAntNnYY1kJPKTIG5JkiPbXX6Po12ZDwo5umQkKN51gK6orVs0QbySuqo+eIwLMM914rwZvA3XWJPCTMk+zD6mq8/uCyGztMh3cn6nO04KjyOYqoMRJ2ZJoxRUlEmlEiUz1JyOI5MrLUHRlpd13zX1g8hMYf/704qsB8wBj/Ys7zLA8YO4yYLkYToncZ3dnYIVmL6tfqrYXGEdqRRwJRHYydom8zonZ7dgLuxn/W2rOTctuOOOrzrjt9mh59ltrlOeopB0eF50i6wjgaCos3jjYj0jzLdOOzdjt8vpc+xWU3Pnvv0W1EirEeFo/xfXIdXWK7n8g35JhkAJcQQrRTvzSAq6Xn6hFCCBFgkviFECLISOIXQoggI4lfCCGCjCR+IYQIMpL4hRAiyEjiF0KIICOJXwghgkybGMCllCoGTvRiuKlASTOG01bI+w4+wfre5X3/sq5a67Qj72wTif/XUEqtbWjkWnsn7zv4BOt7l/fddNLUI4QQQUYSvxBCBJlgSPwzAh1AgMj7Dj7B+t7lfTdRu2/jF0IIUV8w1PiFEEL4kMQvhBBBpl0nfqXUBKXUDqXUT0qpewMdT0tQSs1UShUppbYGOpaWpJTqopT6RimVq5TappS6K9AxtQSlVKRSarVSapPnff8j0DG1JKWUSSm1QSm1INCxtBSl1D6l1Bal1Eal1AldoardtvErpUzATuAc4CCwBrhaa709oIH5mVLqLKAGeENrPSDQ8bQUpVQmkKm1Xq+UigPWAZOC4P+tgBitdY1SKgxYAdyltf4hwKG1CKXU3UAOEK+1vjDQ8bQEpdQ+IEdrfcKD1tpzjX8Y8JPWeo/W2g68C0wMcEx+p7VeBpQFOo6WprUu0Fqv96xXA7lAp8BG5X/aUOO5GeYp7bM2dwSlVGfgAuDVQMfS1rTnxN8JOOBz+yBBkAgEKKWygVOAVQEOpUV4mjs2AkXA11rroHjfwIvAPYA7wHG0NA0sVEqtU0pNPZENtOfErxq4LyhqQsFMKRULzAP+qLWuCnQ8LUFr7dJaDwY6A8OUUu2+iU8pdSFQpLVeF+hYAmCE1noIcB7wO0/zbpO058R/EOjic7szcChAsYgW4Gnjnge8rbX+MNDxtDStdQWwFJgQ2EhaxAjgYk9797vAWKXUW4ENqWVorQ95lkXAfIxm7SZpz4l/DdBLKdVNKRUOXAV8EuCYhJ94TnK+BuRqrZ8PdDwtRSmVppRK9KxHAeOAHwMaVAvQWt+nte6stc7G+G0v0VpfF+Cw/E4pFePpvIBSKgYYDzS5B1+7TfxaaydwJ/AVxom+uVrrbYGNyv+UUnOA74E+SqmDSqmbAx1TCxkB/Aaj5rfRU84PdFAtIBP4Rim1GaOy87XWOmi6NgahDGCFUmoTsBr4TGv9ZVM30m67cwohhGhYu63xCyGEaJgkfiGECDKS+IUQIshI4hdCiCAjiV8IIYKMJH4hjqCUSlRK/daz3lEp9UGgYxKiOUl3TiGO4JnrZ0EwzW4qgktooAMQohV6EujhmfhsF9BPaz1AKXUDMAkwAQOA54BwjIFjNuB8rXWZUqoH8DKQBtQCt2qt2/1oWtF2SFOPEEe7F9jtmfjsr0f8bQBwDcb8KE8AtVrrUzBGS1/vecwM4Pda66HAX4D/tETQQjSW1PiFaJpvPPP9VyulKoFPPfdvAQZ6ZgcdDrxvTB8EQETLhynEL5PEL0TT2HzW3T633Ri/pxCgwnO0IESrJE09QhytGog7kSd6rgGwVyl1BRizhiqlBjVncEL8WpL4hTiC1roUWOm5YP0zJ7CJa4GbPTMobiMILvkp2hbpzimEEEFGavxCCBFkJPELIUSQkcQvhBBBRhK/EEIEGUn8QggRZCTxCyFEkJHEL4QQQeb/ARQ43rwLL3gcAAAAAElFTkSuQmCC", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], "source": [ "# In this example n >> p and it it is basically same as standard regression\n", "# We have to be careful as most of these gLV models are very weakly identifiable\n", @@ -133,7 +164,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 5, "id": "cbab2390", "metadata": { "collapsed": false, @@ -144,7 +175,42 @@ "name": "#%%\n" } }, - "outputs": [], + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "number of species: 5\n", + "specific growth rates: [1.27853844 0.55683415 2.06752757 0.86387608 0.70448068]\n", + "interaction matrix: \n", + "[[-0.05 0. -0.025 0. 0. ]\n", + " [ 0. -0.1 0. 0.05 0. ]\n", + " [ 0. 0. -0.15 0. 0. ]\n", + " [ 0. 0. 0. -0.01 0. ]\n", + " [ 0.02 0. 0. 0. -0.2 ]]\n", + "metabolite production: \n", + "None\n", + "perturbation matrix: \n", + "[[ 0.]\n", + " [-1.]\n", + " [ 0.]\n", + " [-1.]\n", + " [ 0.]]\n" + ] + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEGCAYAAABiq/5QAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAA/ZklEQVR4nO3dd3hUVfrA8e/JzKQ3UgklhA4CChhRgVUULCB2UVRcdFXsbW3g7qqr6GLDVX8qouKCooANFBUFBBWkhSYdKaGEENLbZDLt/P64MyShBshkksz7eZ7z3DtJ5t53JjPvPffcc85VWmuEEEIEjiB/ByCEEKJ+SeIXQogAI4lfCCECjCR+IYQIMJL4hRAiwJj9HUBtJCQk6LS0NH+HIYQQjcrKlSvztNaJh/68UST+tLQ0MjIy/B2GEEI0KkqpXUf6uTT1CCFEgJHEL4QQAUYSvxBCBBhJ/EIIEWAk8QshRICRxC+EEAFGEr8QQgQYSfxCCNEAZec5eeeLQkqt7jrfdqMYwCWEEIFic2YlM+aV8utqK0FB0LNTCH1PD6/TfUjiF0IIP3O7Ncs22Jgxr4S1f1YSEaq4flAUV18QRWJs3adpSfxCCOEnlXY381ZY+eLnUnZlO0hqZuKea2MZ0jeSiDDftcRL4hdCiHpWUOLim19LmfVrGcVlbjq0svDUrfEMODMcs0n5fP+S+IUQop7syLLzxc+lzF9RjtMF53QPY9jAKM7oGIJSvk/4XpL4hRDCh+wOza+rrXy3uIy1f1YSGqwY0i+Say+IolWSxS8xSeIXQggf2LnPzneLy5m3vJyScjcpCWbuvCqWy/pFEB1h8mtskviFEKKOuNyahSutfL2wlI077ZhN8Jee4VzWL5KenUIICqq/5pxjkcQvhBCnSGvN8o023p9ZxI4sB62Tzdx9TSwXnx1BbJR/a/dHIolfCCFOwebMSibOLGLN1kpSEsz862/xnN87vMHU7o9EEr8QQpyErAMOPvymmIWrrMREBnH/sGZc/pdILOaGm/C9fJb4lVKdgenVftQOeBqY4vl5GpAJXK+1LvRVHEIIUVe01vyxrZKvF5ayaE0FwRbFLYOjuX5QtE8HXNU1nyV+rfUWoCeAUsoEZAFfA6OB+VrrcUqp0Z7HT/oqDiGEOFV2h2Z+RjlfLyhl214HUeFBDBsYxXUDo4mPaXht+MdTX009A4HtWutdSqkrgQGen08GFiKJXwjRABWXufhyQSnf/maMsE1LsfD3m+IY1Cec0ODGU8M/VH0l/uHAZ571ZK11NoDWOlsplXSkJyilRgGjAFJTU+slSCGEAKOGP/OXUj75oZhym+bcHmFcMyCKXp3rd4Str/g88SulgoErgDEn8jyt9URgIkB6err2QWhCCFGD1ppfVll5f2YR2fku+pwWyqirY2nXMtjfodWp+qjxDwZWaa1zPI9zlFIpntp+CnCgHmIQQohjWr+9kglfFbJxp512LS28dH8cZ50W5u+wfKI+Ev+NVDXzAHwDjATGeZaz6iEGIYQ4otxCJxO+LmJBhpX4GBOP3xLHxWdHYGrA/fBPlU8Tv1IqHLgIuKvaj8cBM5RStwO7gWG+jEEIIY7E7tB8uaCUj38oxu2Gvw6J5oaLogkLabwXbWvLp4lfa20F4g/5WT5GLx8hhPCLFRsreGtGIXsPOOl3Rhj3XtuMlITAGc8aOK9UCBHwsvOcvPNlIYvXVtAqycy4+xLp061ptuMfiyR+IUSTV1Ds4uM5xXy3qAyzSXHHlTFcd2E0wZam245/LJL4hRBNVqnVzbSfSvhqQSlOl2ZI30hGDIn2yQ3MG5PAfvVCiCapwubmq4WlTJtbgtWmuTA9nFsvi6Gln+541dBI4hdCNCm/rbHy5vRC8otd9D09jNuGxtC+VdMagHWqJPELIZqEvCInb04vZNHaCtq3svDMHQl0bx/i77AaJEn8QohGze3WfLe4jIlfF+FwwZ1XxTJsYBRmU2BeuK0NSfxCiEZr934Hr31awLptlfTqHMLfb4yTdvxakMQvhGh0tNZ8t7ict2YUEBocxOO3xHHpORFNYubM+iCJXwjRqNjsbt6YVsiPS8tJ7xrK6JHxxEU3vpuh+JMkfiFEo7H3gINn389j5z4Hfx0SzS1DYpr0ZGq+IolfCNEo/LbGystT8jGZFP+5NzCnWqgrkviFEA2a06X5YFYRM+aV0rlNMM/ckUDzeEldp0LePSFEg7VlVyWvTi1g+14HV54XyT3XNgvY+XXqkiR+IUSDY7O7+d/sYr6YX0qzaBPP35VAvzPC/R1WkyGJXwjRoKzaYuO1qQVk5zkZ2j+SUVfFEhne9G+OUp8k8QshGoQyq5sJXxXy/e/ltEw0M/7hJHp2CvV3WE2SJH4hhN9t3W3n3+/nklPoYvhFUYy8LIaQYKnl+4okfiGE32it+fa3Mt7+opBmUSbe+Hsy3drJxGq+5uubrccCHwDdAQ38DdgCTAfSgEzgeq11oS/jEEI0PBU2N69/VsC8FVbOOi2Up26NJyZSRuDWB1+fS70BzNFadwHOADYBo4H5WuuOwHzPYyFEAMnMdnDvy/v5OcPKbZfH8J97EyXp1yOf1fiVUtHAecCtAFprO2BXSl0JDPD82WRgIfCkr+IQQjQs85aXM/6zAsKCFS8/kETvLnIBt775sqmnHZALfKSUOgNYCTwEJGutswG01tlKqaQjPVkpNQoYBZCamurDMIUQ9aGsws2b04ymnR4dQvjX3+JJCPB73/qLL991M9AbeEBrvUwp9QYn0KyjtZ4ITARIT0/XvglRCFEf1m+v5IWP8sgtcnHb0BhuuiQak9woxW98mfj3Anu11ss8j7/ASPw5SqkUT20/BTjgwxiEEH7kcmk+/qGYT34oITnezJuPJnNaW+m1428+S/xa6/1KqT1Kqc5a6y3AQGCjp4wExnmWs3wVgxDCf/blOXnxozw27rRzyTkR3D+sGRFh0je/IfB1A9sDwFSlVDCwA7gNoyfRDKXU7cBuYJiPYxBC1LP5K8p5/bMClIJ//S2eC9Ij/B2SqManiV9rvQZIP8KvBvpyv0II/7DZ3fzfDGPahW7tgvnHbTKFckMk/xEhRJ3Yuc/Ocx/ms3u/g5sviWbk0BjMcgG3QZLEL4Q4Jd4bn//f54VEhCleuj+R9K5yd6yGTBK/EOKkVdjcvPppAQsyrJzZJZQxI+OJi5ERuA2dJH4hxEk5UODkHxNy2Znl4PYrYrjx4miC5MbnjYIkfiHECducWck/J+Ris2teuDeRs+XG542KJH4hxAlZuLKccVMKiIsO4pUHk2jbItjfIYkTJIlfCFErWms+mVPCR98W0719CM+NSiA2StrzGyNJ/EKI47I7NK9+ks+8FVYu6hPOozfHE2yR9vzGShK/EOKYcgqcPDsxjy277dx+hTHBmlKS9BszSfxCiKNatdnG85PycDo1z9+VQL8zwv0dkqgDkviFEIfRWjN9XikfzCyidXMLz41KoHWyxd9hiToiiV8IUYPV5uaVTwr4ZZWV83uH88SIOMJCZVbNpkQSvxDioN05Dp6ZmMee/Q7uujqW6wdFSXt+EySJXwhxcL6dd74oJCRY8fKDSfTuLPfCbaok8QsR4IrLXLw2tYBFayvo3TmE0SPlXrhNnfx3hQhgKzfbGDc5n+IyF3dfE8t1F0bJfDsBQBK/EAHI7tBM+raIGfNKSU028+K9zenYWqZeCBSS+IUIMEWlLp56N5fNmXau+Eskd18bS2iw9NoJJJL4hQgg+/OdPPHWAQ4Uunj2zgTO6yUDsgKRJH4hAsTOfXaeeCuXSrubVx5IpEcH6bUTqHya+JVSmUAp4AKcWut0pVQcMB1IAzKB67XWhb6MQ4hAt357JU+9c4CQ4CD++/dk2rWU9vxAVh8NexdorXtqrdM9j0cD87XWHYH5nsdCCB9Zuq6Cx988QGyUibcek6Qv6ifxH+pKYLJnfTJwlR9iECIg/LS0jH++l0ubFAtvPJpM83hp3RW+T/wa+EkptVIpNcrzs2StdTaAZ5l0pCcqpUYppTKUUhm5ubk+DlOIpsXt1vxvdhHjphTQs2MI4x9OopncNEV4+Prw309rvU8plQTMVUptru0TtdYTgYkA6enp2lcBCtHU2Oxuxk3O59fVFVx6bgSP3BiHxSyDskQVnyZ+rfU+z/KAUuproA+Qo5RK0VpnK6VSgAO+jEGIQJJb5ORfE/L4c4+du6+JZdhAmWRNHM5nTT1KqQilVJR3HbgYWA98A4z0/NlIYJavYhAikGzZVcm9L+WwJ8fB2LsTuX6Q3ClLHJkva/zJwNeeD54Z+FRrPUcptQKYoZS6HdgNDPNhDEIEhIWrrLw0OZ/YqCDpuSOOy2eJX2u9AzjjCD/PBwb6ar9CBBK3WzPl+2KmfF9Ct3bBPHdXolzEFcclfbuEaKTKK9z8Z3I+v/9RwSXnGBdxgy3StCOOTxK/EI3QnhwHT7+Xy54DTu4f1oyrB0RKe76oNUn8QjQyyzZUMHZSHqYgxSsPJNFL7pQlTtAxE79S6u+12Ea51vq9OopHCHEUWmumzS3lg1lFtGtp4fm7EmUkrjgpx/vUPA68CxzrHPJuQBK/ED5kd2henZrPvOVWLjgznMdGxBEWInPoi5NzvMT/sdb6uWP9gaePvhDCR4rLXDz9Xh7rtldy2+UxjLhU+ueLU3PMxK+1fuJ4G6jN3wghTs7uHAdPvZNLbqGTf/0tngvSpZ4lTl2tzhWVUg8ppaKV4UOl1Cql1MW+Dk6IQLZmq40HXsnBWuFm/MPJkvRFnaltI+HftNYlGNMuJAK3AeN8FpUQAW7OkjKeeOsAcdFBvP1Ec7q1C/F3SKIJqW2XAG+D4hDgI631WiWNjELUOa01H80u5pMfSjizSyjP3JFAZLhcxBV1q7aJf6VS6iegLTDGM/ma23dhCRF4nC7Na1ML+HFpOUP6RfDw8DjMJqlfibpX28R/O9AT2KG1tiql4jGae4QQdcBqc/PMxDxWbrZx29AYRgyWnjvCd2p7DqmB04AHPY8jABkuKEQdyC928fD4HFZvtfH4LXHcMiRGkr7wqdom/neAc4EbPY9Lgbd9EpEQAWRXtoP7X9nP3lwnL96TyOBzI/0dkggAtW3qOVtr3VsptRpAa12olJIJv4U4Beu22fjnhDzMZvjvI8l0SpWvlKgftU38DqWUCaPJB6VUInJxV4iT9utqKy98lEdynJmX7k8iJUHm3BH1p7aftjeBr4EkpdQLwHXAP30WlRBN2MxfSnlrRiFd04J54Z5EYiLlxiniyBwuB0opzEF1WzGo1da01lOVUisx7pylgKu01pvqNBIhmjitNR9+U8ynP5bQ9/Qw/vm3eEKDpY9+oCuoKGDp3qXsKtrFruJd7C7eze7i3ewq3sW+0n3Mu2UeF7S9oE73ebxpmaO11iVKqTjgAPBZtd/Faa0L6jQaIZoop0vz6icF/LSsnKH9I3nohmaYpI9+QNJaszlvM7O3zmb2n7NZvHsxLu0CwBJkoXVMa1JjUhnYdiCpMam0jG5Z5zEcr8b/KTAUWImnfd9DeR63O94OPNcGMoAsrfVQz0FkOpAGZALXa60LTzhyIRoJq83Ns+/nkbFJ+ugHKpvTxq+7fuX7P7/n263fsqNwBwA9m/dkTP8xXNT+IjrEdaB5ZHOClO/PAo83O+dQz7LtKezjIWATEO15PBqYr7Uep5Qa7Xn85ClsX4gGa0+Og+cn5bEjy8HjI+IY3Fe6awYCrTVb87cyZ9scftz+IwszF1LhrCDUHMrAtgN5vO/jXNbxMlrHtPZLfLVq41dKXQ38rLUu9jyOBQZorWce53mtgMuAFwDv3byuBAZ41icDC5HEL5oYt1sz85cy3p9ZRLBFMfbuRM7pHubvsIQP2V12FmYuZObmmfyw7QcyizIB6BzfmTt738mlHS7l/LTzCbeE+zdQat+r5xmt9dfeB1rrIqXUM8DM4zzvv8ATQFS1nyVrrbM928lWSiUd6YlKqVHAKIDU1NRahimE/+3Pd/LKJ/ms3lJJn26hPHZzHAmx0l2zKSqzl/HDnz8wc8tMvtv6HcWVxYRbwrmo3UU82e9JLml/CW2bnUqDiW/U9tN4pEan410YHgoc0FqvVEoNOMG40FpPBCYCpKen6+P8uRB+p7VmzpJy3v6iEK3h0ZvjGNI3Qtrzm5gyexmzNs9i2oZpzN0+l0pXJQnhCVzb9Vqu6nIVg9oNIszSsM/uapv4M5RS4zGmadDAAxgXfI+lH3CFUmoIxrw+0UqpT4AcpVSKp7afgtFbSIhGrbjMxUtT8lm63kbPjiE88dd4uRF6E2J32ZmzbQ6frvuUb7Z8Q4WzgtbRrbn3rHu5qstV9GvdD1NQ4xmPUdtP5gPAvzB64wD8xHEGcGmtxwBjADw1/se01iOUUq8AIzFu5DISmHXCUQvRgOzPd/Lk/x1gf76T+66L5eoBUQQFSS2/sdtdvJvFuxczf+d8vtr0FYW2QhLCE7it523c2ONG+rbuWy89cHyhtgO4yoHRSqlIrXXZKe5zHDBDKXU7sBsYdorbE8Jvtu2xM+adXCrtbl55MInTO8iktY2RW7tZnb2axXsW8/ue31m8ZzF7S/YCEBUcxVVdruKmHjcxsO1ALCaLn6M9dbXt1dMX+ACIBFKVUmcAd2mt763N87XWCzF676C1zscYASxEo7Zqi42n38slIjSINx5Npm0LmWStscktz+WjNR/x3sr3Dvatbx3dmv6p/enXuh/9WvejR3KPOp8ywd9q+2peBy4BvgHw3HrxPJ9FJUQD93NGOeMm59M62cK4+xJJbNa0EkNTprVm0e5FvJvxLl9u+hK7y875bc7nmfOf4YK0C/zWt74+1frTqrXec0jvBFfdhyNEw/f5/BLe/bKI0zuEMPbuRLknbiNRWlnKlLVTeDfjXTbkbiAmJIZ70u/hrjPvomtiV3+HV69qm/j3eJp7tGce/gcxRuMK0WDkFTn5bU0Fg/tG+GTyM5dL8+6XhXy1sIzze4czZmQ8wRa5iNvQbS/Yztsr3ubD1R9SUllCeot0Jl0xiRu639AgBlP5Q20T/93AG0BLIAv4EbjPV0EJcaJ2ZNkZ83YuuUUuflxazvN3J5BYh4OmyircjJ2Ux/INNq67MIq7r4mVnjsNmNaan3f+zBvL3mD21tmYgkwMO20YD539EGe3Otvf4fmd0rrhj41KT0/XGRkZ/g5DNFArN9t4dmIuoSFBDL8oiknfFhMeGsRzdyXQNS3klLefnefkqXdz2Zvj4KHhcQztL/PtNFS55bl8/MfHfLDqAzblbSIxPJG70+/m7vS7aRHVwt/h1Tul1EqtdfqhP69tr552GDX+czAGcC0BHtFa76jTKIU4QT8sKWP81ALaNLfw4r2JJMWZ6d0llH+8m8sjrx/g8RFxDDwr4qS3v26bjacn5uF2w0sPJNG7s3TXbGhcbhdzd8zlw9UfMmvzLBxuB+e0OoePrvyI4d2HE2qW/9mhansu/CnGqN2rPY+HY8zNL+dMwi+01vxvdjEf/1DCmV1CeebOBCLDjHb9ti2CeeeJ5jz7fh4vfJRPZraD24bGnHDTzE/Lynltaj7JcWZeuCeR1smNv/92U7K3ZC8frPqASasnsadkDwnhCTzQ5wH+1utvdEvq5u/wGrTaJn6ltf642uNPlFL3+yIgIY7H4dS8+kk+c5dbGXxuBI/cFIf5kJuaxEaZeOXBJN6YXsDUOSVkZjt4amQ8YaHHv+jrdmsmfWvcKatX5xCeuSOB6IjGMxy/KXNrN/N2zOOdFe/w7dZv0VpzcfuLGX/JeK7ofAXBJhlLURu1TfwLPHPnT8No6rkB+M5zUxXkTlyivpRZ3Tzzfi6rt1Ry2+UxjLj06Dc1sZgVj94UR7sWwbzzRSH3vryfZ+9MpE3K0WvuFTY3L07OZ/HaCob2j+TBG5oddlAR9S/fms9Haz5iQsYEthduJzE8kSf6PsGoM0c1yNkvG7paXdxVSu08xq+11vq4d+I6FXJxVwDkFDgZ83Yuew84eHxEPBedXfu2+9VbbIydlIfNrnl8RBwDzjz8uTkFTv45IZedWQ7uva4ZVw+IlJk1/UhrzeI9i3lv5Xt8vuFzKl2V9E/tzz3p93Bt12sJMZ/6hfum7mgXd6VXj2gUvHPi2OxunhuVSK+TuMiaW+Tk3+/nsXGnnWsvjOKuq2MP1uY37Kjk6fdysTs1T9+ewFmnNexpdZuyIlsRH6/9mPdWvseG3A1EBUcx4vQR3JN+Dz2Se/g7vEblVHv1DAPmaK1LlVL/BHoDz2utV9dxnEIcZsXGCp59P4+o8CDePIU5cRJjzbz+SDLvfV3Elz+XsmWXnX/dHs+aLZW8OjWfxGZmXr8nkdTmchHXH1Zlr+Kt5W8xff10KpwVnNXiLD64/AOGdx9ORPDJ98wSh6ttU88fWuvTlVL9gf8ArwJPaa3rpVeP1PgD1/eLyxj/WQHtWhjdNevqTlY/Z5Tz6tQCTArKbZqeHUN45s4EYiLlIm590lrz0/afeOX3V5i/cz6RwZHc3ONm7jrzLnql9PJ3eI3eKdX4qZqX5zLgXa31LKXUs3UVnBCHcrk1//u2mKk/lnDWaaE8c0cC4bXokVNbF6ZH0K5lMC9NyadrWjD3XicXceuTw+Vg+obpvPL7K/yR8wctolrwykWvMOrMUUSHRPs7vCavtok/Syn1HjAIeEkpFcKRb8coxCkrLnMxdlI+KzfbuKxfBA8NP7y7Zl1IS7Hw7pPN63y74ujK7eW8v+p9xi8Zz56SPXRL7Mb/rvwfN/a4Ubpi1qPaJv7rgUuBVz03Wk8BHvddWCJQbcqs5N/v51FY6uKxm+MY0k+mR2gKim3FvL3ibV5f+jp51jzOa3MeE4ZOYHCHwdJzyg+Od8P0VVrr3lprK/CV9+da62wgu/rf+DZM0dRprfnmtzLe/ryQhFgTbz3WnE6pUgNs7PKsebyx9A3eWv4WxZXFDOk4hH/85R/0bd3X36EFtOPV+Lsqpf44xu8VEFOH8YgAZLO7ef3TAuYut3J2t1DG3BovI2UbuaySLMYvGc+ElROocFRwTddreOovT9E7ReqIDcHxEn+XWmxDbsgiTtqubAfPfZh3cD6dmy+NlumOG7FtBdt4efHLTF47GZfbxU09bmJM/zEBd6OThu6YiV9rvetkN6yUCgV+BUI8+/lCa/2MZ5qH6UAakAlcr7UuPNn9iMZJa82PS8t5c3ohIcGKcfclyqCpRmzN/jWMWzSOzzd+jiXIwh297uCxvo/JdAoNlC9vFFoJXKi1LlNKWYBFSqkfgGuA+VrrcZ75f0YDT/owDtHAWG1u/jutgHnLrfTsGMJTt8XXWf98UX+KbEXM3jqbqeumMmfbHKKCo3i87+M8fM7DNI+U3lINmc++bdoYGVbmeWjxFA1cCQzw/HwysBBJ/AFj2x47z32Yx75cJ7d6mnZM0rTTaOSU5TBryyy+2vQV83fOx+l20iKqBWMvGMt9fe4jNjTW3yGKWvBpNUspZQJWAh2At7XWy5RSyZ5eQWits5VSSUd57ihgFEBqaqovwxT1QGvNN7+W8c6XhURHmHjtoSTO6CQ3yGgs5mybw4u/vcii3YvQaNo3a88j5zzCNV2voU/LPgQpGdbTmPg08WutXUBPpVQs8LVSqvsJPHciMBGMKRt8E6GoD+UVbl6dWsAvq6z06RbK6L/GExslvXYagz3Fe3j4x4f5atNXtG/WnmfOf4Zrul5D96Tu0v++EauXhlXPoK+FGIPAcpRSKZ7afgpwoD5iEP6xI8vOs+/nsS/PyZ1XxXLDoCjptdMIOFwO/rv0v/z7l3/j1m5evPBFHu37qIyubSJ8lviVUomAw5P0w/BM9wB8A4wExnmWs3wVg/CvOUvKeGNaIRFhymja6ShNO43Br7t+5d7v7mVD7gau6HwFb1z6Bmmxaf4OS9QhX9b4U4DJnnb+IGCG1nq2UmoJMEMpdTuwGxjmwxiEH1Ta3bw5vZAflpTTs1MI/7wtgbgYadppyLTWLNq9iPFLxzNz80zaxLRh1vBZXNH5Cn+HJnzAl716/gAOm1dVa50PDPTVfoV/7ckxBmRt3+vg5kujuXVojPTaacDsLjszNszg9aWvsyp7FXFhcTx93tM82f9Jwi3h/g5P+Ih0nhZ1QmvNd4vLeeeLQoItiv/cl8jZ3WRAVkOVb81nQsYE3l7xNtll2XRJ6MKEyyZwyxm3SMIPAJL4xSkrLnPx2tQCFq2toHfnEEaPlAFZDZXNaeO/S//Li7+9SKm9lIvbX8ykKydxcfuLpUtmAJFvpzglqzbb+M/kfIrLXNx1dSzDBkqvnYZIa8209dMYM38Mu4p3cXmny3lx4It0T6p1D2vRhEjiFyfF4dRM+qaIGfNLaZVk5oV7ZBrlhur3Pb/z9x//zrKsZfRs3pNJV07iwrYX+jss4UeS+MUJ273fwdiP8ti2x8Hl/SO557pYQoOlmaCh2VG4g9HzRvP5xs9JiUzhoys/4pbTb8EUJD2sAp0kflFr3gu4b39uzKj53KgE+veUC4ENTZGtiBd/e5E3lr2BOcjMs+c/y2N9HyMiOMLfoYkGQhK/qBW5gNvwOVwOJq6cyLO/PEu+NZ9be97K2AvH0iKqhb9DEw2MfHPFca3aYmPc5HyKSuUCbkOktea7P7/j8bmPszlvMxekXcD4S8bTs3lPf4cmGihJ/OKoHE7NR7OLmT63hFZJZsbeLRdwG5I/8//k4z8+5pM/PmFn0U46xXfim+HfMLTTUJlATRyTJH5xRDv32fnP//LZttfBZf0iuPe6ZoSFyAVcf8uz5jF9/XQ+/uNjlmUtI0gFMbDtQJ674Dlu6HYDFpPF3yGKRkASv6jB7dZ8vbCUiTOLiAgN4vm7E+h3ulzA9bc/8//khd9eYOq6qTjdTk5PPp1XLnqFm3rcJG344oRJ4hcH5RY5eXlKASs32zineyiPjYgnLlq6/vnTlrwtjP1tLJ+u+5QQUwj3pt/L7b1v5/Tk0/0dmmjEJPELABauLOf1zwpxODV/vymOy/pFSDuxH23K3cTY38Yybf00QkwhPHLOIzze93GSI5P9HZpoAiTxB7jiMhdvzihkQYaVrmnBjLk1nlZJ0k7sLzsKd/D0gqf5dN2nhFnCePTcR3ms72MkRRzxDqVCnBRJ/AHs9z+sjP+0gJJyN3+7PIYbL47GZJJavj/klOUw9texvLfyPcxBZh7v+ziP9X2MxIhEf4cmmiBJ/AGorMLN258X8uPSctq1tPDS/Um0byXdNP2h2FbMa0teY/yS8dicNu7sfSf/Ov9fcsFW+JQk/gCzYmMFr35SQH6JixGDo7llcAwWs9Ty61uls5J3VrzDC7+9QH5FPjd0u4HnL3iejvEd/R2aCACS+ANEfrGL974uZN5yK22am/n3qGS6pIX4O6yA49Zupq+fzlM/P0VmUSYXt7+Y/wz8D71Tevs7NBFAJPE3cU6X5qsFpUz+rhinS3PzpdGMuDSaEJlNs979kvkLj819jIx9GfRs3pO5t8xlULtB/g5LBCBJ/E3Yqs023ppRwK79Ts7uFsr9w5rRUnrs1LuNuRt5ct6TzN46m9bRrZl81WRGnD5C7ngl/MZniV8p1RqYAjQH3MBErfUbSqk4YDqQBmQC12utC30VRyDKLXTyzpdF/LLKSkq8ibF3J9BXRt/WqzxrHp9v+Jyp66ayeM9iokOiGTdwHA+e/SBhFrkXsfAvX9b4ncCjWutVSqkoYKVSai5wKzBfaz1OKTUaGA086cM4AobLrfn2tzI+mFWE0wW3Do3hhkFR0qxTT6wOK99s+Yap66YyZ9scnG4n3RK78eKFL3LnmXeSEJ7g7xCFAHyY+LXW2UC2Z71UKbUJaAlcCQzw/NlkYCGS+E/Zjiw74z8tYONOO2d2CeXhG5vRMlGaderD3pK9jF8ynvdXvU+ZvYyWUS155JxHuLnHzZyefLqMgBYNTr208Sul0oBewDIg2XNQQGudrZQ64pBEpdQoYBRAampqfYTZKNkdmo9/KGbaTyVEhgcxZmQ8g/qES7KpB5vzNvPy4pf55I9PcGs3w7sP5/Zet3Nem/Pk9oaiQfN54ldKRQJfAg9rrUtqm5C01hOBiQDp6enadxE2Xis2VvDWjEL2HnByyTkR3H1NLDGRgZ1wtNY+P+ityFrBuMXj+HrT14SYQ7jrzLt4tO+jpMWm+XS/QtQVnyZ+pZQFI+lP1Vp/5flxjlIqxVPbTwEO+DKGpmjLrkren1nEqi2VtEg08+qDSfTuEurvsPxmf9l+ZmyYwafrPmVtzlpuOf0Wnuz3JO3j2tfZPpxuJ7M2z+LN5W/y665fiQmJ4am/PMWDZz8o8+iIU+NyQUkJWK1GqaioWrdaoU8fSKrbz5gve/Uo4ENgk9Z6fLVffQOMBMZ5lrN8FUNTsyfHwaRvi/lllZWYyCDuuy6Wy/8SRbAl8Jp1CisK+WrTV3y2/jMWZC7Ard2ckXwG13a9lilrp/Dh6g8Z3n04Y/qPoXtS95PeT741nw9WfcDbK95mT8ke2sS04ZWLXmHUmaOIDomuw1cUQNxucDjAZDLKkc7QtDb+xlucTggKMopSVetgJM38fKMUFFQty8qMJHpoUQoiIyEqquYyOhri4iA+vubSYoGiIti7t2bZt8/YnstlxFd9qRSYzUaxWKrWHY7D4ywqMl7v0Xz/PQweXKf/AqWPtcNT2bBS/YHfgHUY3TkBnsJo558BpAK7gWFa64JjbSs9PV1nZGT4JM7GIL/YxZTvivnu9zKCLYrrB0YxbGA0EWGB1VtHa83SvUt5c/mbfLXpK+wuOx3iOnBj9xu5sfuNdE3sCkB2aTavLXmNCRkTKHeUc1WXq3iq/1Oc1fKsWu8nY18G7618j6nrpmJz2riw7YU82OdBhnYa2vDb77WG4mLYvx8OHKiZnKonKO+6t1R/7E22TmfVOkBMDDRrVrNER4PNZiTg0tKqUlICublGHNVLTo6xDy9vkjR53leHo+bvT5ZSEBoKYWFVJTTUeH/KyqpKZeWxt2OxGDEduu2kJIiIqIrdu/S+jurvnbeYTIcfWOLijPcxIgLCw6tKWJix7NTJeN9P6i1QK7XW6Yf93FeJvy4FauLXWvP97+W8+2UhdodmaP9IRgyOCbibo1Q6K5mxYQZvLn+TjH0ZxITEMPKMkYw4fQTpLdKP2qafb83nzWVv8ubyNymyFdE1oStDOw3l8k6Xc27rczEH1TzhXX9gPdPWT2Pa+mlsL9xOmDmMW06/hQfOfuCUzhpqKCuDXbsgM9NY2mxGYgkONpbeYrcbSTs3t+aypMT42+BgCAmpWgYFQV6ekVT37zeeX5csFiNheg8AtWU2Q3IyNG9es4SH1zzoeA883n0dWsxmY/9ud9XSW2JiDk+m8fHGPmpzvcfhMP4vxcU1a+Les4iyMuM1tGpVVVJSjLgaOEn8jUxukZPXphawfIONnh1DePTmuIAbdbuvdB8TV05kQsYEcspz6JLQhQf7PMgtZ9xCZHBkrbdTUlnC5DWTmbVlFr/s+gWn20lcWByDOwxmSMchbC/YzrQN09iYu5EgFcSFbS9keLfhXNPlapq5g42EW714k0F5uVG86xUVxg69zRDeZgmA7Gwj2eflndibEBxs1CwTE41ldLSRICsrjeTuXTqdRrKrnlyTk40SHl5VI61eKz3auslUM+F6X4fWRptzYWHNUlJi1E6jooz4oqKqSmxsVZOMqHeS+BsJrTVzl1v5vxkFOJxw51WxXHV+JEFBgdGOX+GoYNaWWUxeO5mftv+EW7u5rONlPHj2gwxqN+jY0xxYrZCVZSSisjKjuaH60m6n2FHGT87NfKu38L3aRr4ykvVfrIkMz0nk2swwknMrqmp/3mR+KKWMU/OICKN9OCLCSH5K1ayNemuozZtDWppR2rSpWkZEVLVj2+1V6xaLkeijompXaxXiCI6W+GWungakoMTF658WsPiPCrq1C+bJvwbG3bC01izavYgpa6cwY+MMSipLaB3dmtH9RnNbr9voENfBSKC5ecYFtX37YM8eowa9c6exzMw0mjmOIwYYBgwzmXAFm1nVMozmOpLWpmYQE2E0G3RNq7rQl5R0eDmRZgQhGiBJ/A3A9r12fvi9jJ+WlVPp0Nx9TSzXXhiFqQnX8iscFSzY+TPfb5jJ7D+/Y1dFNhEqhOtMPfirvQMDNoQTtHA97L/ZSPT79x/evmw2V9WeL7/cWLZubSRvb1ODt9dGRITRFu5tvlAKE1C7y71CNC2S+P2k1Orm5xXl/LCknK277VjM0P+McP46JIY2KY20lq+10eZcVmZ0USssNJae9cz87XyXv4TvXZv5OSoPm0kTboeBO+G5jXDNpkoi7RlgXlPVpp2UBKedZlxMa9GiZmnZsqoHhRCi1iTx17Otu+18Pr+E39ZUYHdo2reycP+wZgzqE050RANPYm43bNgACxfCL7/A1q1VFzjLy402drf7sKdtToBnBsCM7kAUtC8LZlR+W4aE9eD8lHMJ7dUe7kupSvYxMdKMIoQPSeKvJ39sszF1TgkrNtqICFMM7hvB4HMj6ZTagO91a7fD+vWwaJGR7H/91ejRApCaCj17Gs0o4eFVFzq9pVkzdobZ+HfBV3yc/RNh5lCeSr+Pkel30Cm+kz9flRABTxK/D2mtydhk45M5JazbVklsZBB3XBnDledFNbzBV04nbNoEGRlVZe3aqsEt3nb088+HAQOMx0eRVZLF2F/H8sHqDzApEw+f8zCj+48mMSKxPl6JEOI4JPH7gNaaJesqmPJ9CVt320mMNXH/sGYM6RdBaH3Oje90Gv3HDx1mnpd35CHu3tGSUVFw5pnwwAOQng7nnGNcRK2m0lnJ1vytZBZlsqt4F5lFmQfX1+Wsw6Vd3Nn7Tv7xl3/QMrpl/b1mIcRxSeKvY6s22/jwmyI2ZdpJSTDz6M1xXNQnwjfz6RQVwcaNRtfGvXsPX+7ff3ibe2io0ZbuHd3Yo0fVaMfTTjMSfceOhw26cbldrNm/hvk75zNvxzwW7V5EhbOqj3uoOZQ2MW1Ii03jnvR7ePDsB2nbrG3dv2YhGjin07hE1ZD7HUjiryPrt1cy6dsi1mytJDHWxN9viuPScyMwm+ow4eflGe3sv/xiLNeurTm5U0SE0Z2xVSu4+OKq9VatqtZjY2t14TTPmse6nHWszVnLb7t/Y8HOBRTajDtknpZ4Gnf0voO+rfvSNrYtabFpJEUkHXM6ZG+YJ3LN1uEwrhe7XMbxq/pSqZrTr/jyWrDWVdPX2O3GLAveqWgOLd6fV1+WlRk9SaOjqwa2etdNpiNPn2M2G/+q6lPixMYazz10xgJvqR5v9XVv3N7iHexrNtccZBsdbfR+NZmqBiRXH/9WXn749D4ul7GP0NCal3i8Y9q8l3/MR8g03hPSrCyjnpKVZdRlvNPrVF96/8eHFjDG2FUfQO1dhoVBQsLhJTi4ajqh6q/PZquaOcM7K4b38aH79A5kzsqCHTtg+3ZjuWOHMROH1jU7piUnG8vY2Jr78C5DQ43YEhOrSqgPJ9yVkbunQGvN5kw7U74vZtkGG82igrjpkuhTmzFTa6Ppxfsp2rnT+FQtXWr0qAHjE33uuUZ7+5lnGhdavf3XTyIDFtuKmbNtDiuzV7J2/x+s3bmHnKwQKE6F4lRigxPo1ro1vdp0oG+H00hrHk90tJE8tm8/vOzbVzVg1Vu8YmKMD3h8fNUXMS7OSPCHTktTVFT71+A9CISEHHlKF62NBBQdbcQQE1OVfB0OI0kXF9dcWq1Vg2lP5GtiMtWcvSAy0ki21Q8Ox5sXzJ+8Sa0uWSxVc49FRBjv7ZFOSOtKaKjxHtdXektIgHbtjNK2rXHCfOhMHwcOGP/72oqMNA4AkyYZl9VOhozcrUNWm5t5y8uZvaiMbXsdRIUbF22vHhBFWMgJtOGXlxu19lWrjLJmDWzbdvinIzERevWCm24ykv1ZZxnVhGqcTsjaXTWINSvLaLYvLKxaeqdVMZnAZHFS4S6mxJlHiTMPrZLBehtBJW1wV9a8MXsRsNhT/u8oLyUhAdq3h379jBMLzxipg9O8eGcyKCoyjmt5ecYXf8MG43F4eNWUNL16Va17a6BBQTWXbnfNmXa905hXVh599t7y8qqkXlho1MyKi42k5D0YJCcbkyFGRxsxVa/1eZfVa+/VizfZe2duOJbKSuPf7HbXnEbHO2WOw1H1P/MOiSgsNGqm1acCqv5aq++z+nr1Gqx3XjeLxfjMVJ9M01uczppj37zr4eGHT+njnVXZZqs5fVH1Hr7epXfdWxuvPudZy5ZVJ6SVlcb2KipqLqtXJKpXKMLCqs4wvHGaTMbZSGGh8VnLyzMqFLm5xntb/UzH+zpDQ2vOnFH9LKn6/qqvt2hhJProWs7Q7Z2Ruvr2HQ7jPakeY/WS4INbNUuN/wRs3W3n299KmZ9hxVZp9MEf2j+SQWdFHL+XjtUKq1fDihVGj5lVq2Dz5qpPUGIixT36s7/lmeTEduJAWBo55pbkuBI4UBR8sIv8oc0ehYVGot+z5/CZbMPDq5oJYmLdmCNKqAjaz66ivRwoKUQ7LYQFxZAY0pJmwcm0bRlFWpsg2rQxruWmpholOLiqtupNnN4DSLt2RsKv7QdfCFF/pMZ/korLXCxYaWWOZ4RtiEVxQXo4l/ePpEta8JHbtR0Ooyq7fDnu5RmULNtE0cZ9FLmjKKQZ+XGd2J3yKJlndWcXbcgsjWfXPjPFPx++LaWMI35ERFVtt3rNNzraqGV75/9KS4PEFuXsU8vZWLiKNTlrWLN/DUtzN+HSxpGhc3xnbu96Ldeedi29mveq1a0KT3I6cCFEAyQ1/iNwuTQrNtqYs7SMJesqcDihQysLg/tGclGfCCLDq9XutTaaZ1asQC9bztpfivhxQyvmOAeyml6UEI3myGcDUVE1J2tMTTVOHb0Xg5KTjbbwI10Yq87pdpKxL4O52+cyd8dcluxdgtNtzGvTMqolPZv3PFh6Ne9Fu2bt5GbsQgQAqfHXQkGJiy/ml/DjsnIKS9zERAZxxXlRXHpOBO1bBRtJPjMTVq7EvWIlRUs3k7t6LytLO/Ijl/AjY8ghGYAz2hRyc39IaKeI9fTI8JZmzYxkX8sONoBxIbmksoS9JXvJKs0iqySLrNIsVmavZMHOBRRXFqNQ9E7pzaPnPsoFaRfQO6W3DJoSQhxGEj9QXuHm8/klzJhfit2hObdHGJeeE0E7Uznrvt7KF+8Ws2ZjMFl5weQ5Y8njPPK5GjdVHXXjY51cfGkQl1xq9KRMSWl2UrE4XA7+LPiTDQc2sDF3IxtyjWVmUSbljvLD/r5tbFuu73Y9g9oN4sK2F5IQ7oMrQUKIJiWgE7/doZm9qIyPfyimuMxNhxA7MTsLWfN1EB/ltCDLlQIYNeY2lizSEsro1sJMQltI6KhJaG60v3fuDL17m497o6Fyezn7Svexr3Qf2WXZxrI0m+yy7IOPtxVsO9hMo1C0j2vPaYmncXH7i2kZ1ZKW0S0PLltEtSDU7MPOvkKIJslniV8pNQkYChzQWnf3/CwOmA6kAZnA9VrrQl/FcERak7V0Nx9MKWBJZTT2UAule01sWZrMgtwQgnDRxbyNAS3+pNfpm+k1KJ6e17QjLvXI0w6U2cvYV7qPJXtzybPmkWfNI9dqrOeU5xxM9PtK91FSWXLY80NMIaREpZASmcJpiadxVeer6JbUjW6J3eiS0IUwS5iv3xEhRIDx2cVdpdR5QBkwpVrifxko0FqPU0qNBppprZ883rZO9eKu3rGTxZ9tYNqSSNarFujmYZgsmvJcM6wppnNQFt26VdC5r5n2g+KwtIrB6XbWKA6Xg6zSLLYVbGN74Xa2FWxjW8E29pftP+I+wy3hJIYnHqyZt4hsYSyjWpASlWIsI1OIDY2VC61CCJ/wyz13lVJpwOxqiX8LMEBrna2USgEWaq07H287J5v473noZTL2pxMS0Y7gUKMdxlrqotCxhfzYL8hPnAxm53G2crgWUS3oENeBDs060CGuA61jWpMYnkhiRCIJ4QkkhCcQbgk//oaEEMKHGkqvnmStdTaAJ/knHe0PlVKjgFEAqampJ7WzPlHbeaTrFH7YOZgFhd0w9dhL9w4HCLeEE2ZuTYj5eSxBFoJNwVhMnmWQBYvJgjnIjDnIjEmZDq6nRKXQrlk7SepCiEatwV7c1VpPBCaCUeM/mW1cOvQvxOX+wUPRr/IQQPw50GY4pA6D8BZ1GK0QQjQe9X03kBxPEw+e5QFf7izlnBGEXL4ErtgBPceB2warHoaZrWD+hbBjMjgP7yIphBBNWX0n/m+AkZ71kcCsetlrZFs47UkYvBou2wQ9noHyPbD0VviqOSy7A3IX199UfkII4Ue+7NXzGTAASABygGeAmcAMIBXYDQzTWhccb1s+mbJBa8hdBDs+gt0zjJp/VCdoNxKanQkRqRDRBszSni+EaJz80qunrvh8rh5HGez+3DgI5P5W83chCcYBIKINxHSH+D5GCZWpEIQQDVtD6dXTMFkiof1tRqnIhtLtUL4LrLuMZfluKN4Ae2eC9tw5IiKt6iAQ1xtiukHoUTspCSFEgyGJ/1BhKUah/+G/c5RB4SrIX+4py4xmIq+QROMAENMNYrtBVEcIbgaWGLDEQnAMBFnq65UIIcQRSeI/EZZISDrPKF62A1C41jgjKN4AReth5xRwHuUea6ZwCImH2DMg/iyISzeW0nQkhKgnkvhPVWgSpFxkFC+twboXynaAoxjsRcbSu27bDwWrYN93gOcaS3gqxKdDeBsIS4bQ5hCaXK0kQZD8u4QQp04yiS8oBRGtjXIsjhIoWA0FKyA/w2hG2jcHXNYjbDMIQlMgvFXNYo4CVwW4bJ7iWccNwXHG2UVwvHGROiTe8zgWzNEQZDp8P0KIJk8Svz9ZoiH5fKNU5ygDW061sh+s+6Bir3EmUbwBsuccefBZUAiYPDN6Ooo5eEZxJOZI4yBgiTGWoc0hvLXnoNK6at0cAdpVrTiNJYCyQFCwce0iyLOugmoehFwVnmIHU6jRZGaK8CzDjL8XQtQbSfwNkSXSKFHtj/43WhuJ3VluJE9TGJhCaiZRtwvshVCZB/Z8qMw31h3FYC+uan7yNkEd64DiS+YI4yL4wTMT79lJgnEwOfRsxmXznAF5msHCvM1izY3tuCvBaTVeh8tqrLus4KoE7QB3taIdxoEvrEVVCU2uOhtyu4yDb0UWVOwDaxbYCwDluX1aUNUSjO257OA+pEC156iqdZfNOPNzlBjXhRwl4Cg1DqzmSKNYIqvWzZFgiTry7yzRhxdTuPE+VOYa16O8S9sBzz7CPZ+dcM+6Z9zKoZ8N72cNPJ8x5VkGgTIZ75u323NkGoS1rGqadJYbPePKd4N1N1j3GL3jvB0pQj3LsOZGxaEyz+hd5y22bOOzC8a+lMnYtzIZ+3fbjffRfcjnRJmqfTdCa35PgiyeSku1orzp0FNZ8nZ1V0HGexnczDiLDm5mFHOYsR/bgUMqageM1xHavOZnMyTO2JazwniN1YujuOo11Shm45piWPO6/crV6dZE/VHKqKUHxx79b4JMEJpglNryHlCse4yzC+se44OqTMb2lLnqQwlVydNlr0qq2nXIF82zHhRsfFGc5eAs85RyI9E5CsHmOUAVrjK+6PZCjC+hMp5vDoMgz7a00/iSuSpO/j08GhVkfFFRxtmW9+zmRBw8+/EmFeVJJJ7iXTeFGknFHG0k9Mj2RvNdkMnz3njeJ+teY+korXr/jnU2V/Viavl3x3ktlhjjAA2e2N1VS7e9KjEf3K3nYOAs9xwoq//Oe5B0HyFc05Hf74P7dhnP8559gvGZPPg585SgEOP3Nc44K6odhOtAkMX4vNeWMnsqMif4mR3wA4RdemLPOQ5J/KKm6geU2B7+jcXt+XJ7E+ehtDYSoC0HKvYbS3uB8cU3R1TVYs0RngNPSLXanbemZzbOfir2VRVrllHDRxs11/CWnrMBzzLEeyB1VyVwbxLzNnv5+h4L2m0kEO+BwVnqOYCWHF4sUUbngJBEY+ldD7JUnQ0dXFYYr8cSY3Q/tsR4muOO83pcNk+t3jv2JdN4bIk0Oi5EpFYtw1oAyqjp2rJr1u5dFVU15epnBOaj3JBI6xN7r7XnQHXoWZ/bAW5ntW15lkoZn0NHsVERqV4cRcZB+tBOGKFJRsWkYr9Rcai+dNurndVWK5Zo4333NqO6qzWtRpzc7MTHIiN3hRCiiTrayF25qiaEEAFGEr8QQgQYSfxCCBFgJPELIUSAkcQvhBABRhK/EEIEGEn8QggRYCTxCyFEgGkUA7iUUrnArpN8egKQV4fhNBbyugNPoL52ed1H10ZrfdjNPhpF4j8VSqmMI41ca+rkdQeeQH3t8rpPnDT1CCFEgJHEL4QQASYQEv9EfwfgJ/K6A0+gvnZ53SeoybfxCyGEqCkQavxCCCGqkcQvhBABpkknfqXUpUqpLUqpbUqp0f6Opz4opSYppQ4opdb7O5b6pJRqrZRaoJTapJTaoJR6yN8x1QelVKhSarlSaq3ndf/b3zHVJ6WUSSm1Wik129+x1BelVKZSap1Sao1S6qTuUNVk2/iVUiZgK3ARsBdYAdyotd7o18B8TCl1HlAGTNFad/d3PPVFKZUCpGitVymlooCVwFUB8P9WQITWukwpZQEWAQ9prZf6ObR6oZT6O5AORGuth/o7nvqglMoE0rXWJz1orSnX+PsA27TWO7TWdmAacKWfY/I5rfWvQMFx/7CJ0Vpna61XedZLgU1AS/9G5XvaUOZ5aPGUplmbO4RSqhVwGfCBv2NpbJpy4m8J7Kn2eC8BkAgEKKXSgF7AMj+HUi88zR1rgAPAXK11QLxu4L/AE4Dbz3HUNw38pJRaqZQadTIbaMqJXx3hZwFREwpkSqlI4EvgYa11ib/jqQ9aa5fWuifQCuijlGryTXxKqaHAAa31Sn/H4gf9tNa9gcHAfZ7m3RPSlBP/XqB1tcetgH1+ikXUA08b95fAVK31V/6Op75prYuAhcCl/o2kXvQDrvC0d08DLlRKfeLfkOqH1nqfZ3kA+BqjWfuENOXEvwLoqJRqq5QKBoYD3/g5JuEjnoucHwKbtNbj/R1PfVFKJSqlYj3rYcAgYLNfg6oHWusxWutWWus0jO/2z1rrEX4Oy+eUUhGezgsopSKAi4ET7sHXZBO/1toJ3A/8iHGhb4bWeoN/o/I9pdRnwBKgs1Jqr1Lqdn/HVE/6Abdg1PzWeMoQfwdVD1KABUqpPzAqO3O11gHTtTEAJQOLlFJrgeXAd1rrOSe6kSbbnVMIIcSRNdkavxBCiCOTxC+EEAFGEr8QQgQYSfxCCBFgJPELIUSAkcQvxCGUUrFKqXs96y2UUl/4OyYh6pJ05xTiEJ65fmYH0uymIrCY/R2AEA3QOKC9Z+KzP4GuWuvuSqlbgasAE9AdeA0Ixhg4VgkM0VoXKKXaA28DiYAVuFNr3eRH04rGQ5p6hDjcaGC7Z+Kzxw/5XXfgJoz5UV4ArFrrXhijpf/q+ZuJwANa6zOBx4B36iNoIWpLavxCnJgFnvn+S5VSxcC3np+vA073zA7aF/jcmD4IgJD6D1OIo5PEL8SJqay27q722I3xfQoCijxnC0I0SNLUI8ThSoGok3mi5x4AO5VSw8CYNVQpdUZdBifEqZLEL8QhtNb5wGLPDetfOYlN3Azc7plBcQMBcMtP0bhId04hhAgwUuMXQogAI4lfCCECjCR+IYQIMJL4hRAiwEjiF0KIACOJXwghAowkfiGECDD/D2gQM5yxmbhVAAAAAElFTkSuQmCC", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], "source": [ "set_all_seeds(1234)\n", "\n", @@ -165,11 +231,14 @@ "mu = np.random.lognormal(0.01, 0.5, num_species)\n", "\n", "# construct perturbation matrix\n", - "epsilon = np.array([0, -1, 0, -1, 0])\n", + "npert = 1\n", + "epsilon = np.zeros([num_species,npert])\n", + "epsilon[:,0] = [0, -1, 0, -1, 0]\n", "\n", "# instantiate simulator\n", "simulator = gMLV_sim(num_species=num_species,\n", " num_metabolites=num_metabolites,\n", + " num_perturbations=npert,\n", " M=M,\n", " mu=mu,\n", " epsilon=epsilon)\n", @@ -180,13 +249,18 @@ "init_species = 10 * np.ones(num_species)\n", "init_metabolites = 10 * np.ones(num_metabolites)\n", "\n", - "# perturbation\n", - "tp = 2\n", + "# perturbation information encoded in a function\n", + "def pert_fn(t):\n", + " if 2.0 <= t < 2.2:\n", + " return np.array([1])\n", + " else: \n", + " return np.array([0])\n", + "\n", "\n", "times = np.arange(0, 5, 0.1)\n", "yobs, sobs, sy0, mu, M, _ = simulator.simulate(times=times, \n", " sy0=np.hstack((init_species, init_metabolites)),\n", - " tp=tp)\n", + " u=pert_fn)\n", "\n", "\n", "# add some gaussian noise\n", @@ -214,7 +288,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.11.7" + "version": "3.8.8" } }, "nbformat": 4, diff --git a/examples/run_gMLV_sims.py b/examples/run_gMLV_sims.py index f343ced..c6902ae 100644 --- a/examples/run_gMLV_sims.py +++ b/examples/run_gMLV_sims.py @@ -1,6 +1,6 @@ ''' Simulation script for gMLV model with perturbations and metabolites. -Usage: python run_gMLV_sims.py +Usage: python run_gMLV_sims.py Example: python run_gMLV_sims.py outputs/ 100 The script will create a folder named outputs/ and save the results there. @@ -8,26 +8,21 @@ a different number as the second argument. ''' - import logging -from time import time -import math -from gMLV import * -from numpy import linalg as la -import copy -from scipy.integrate import odeint -import matplotlib.pyplot as plt import random import os import sys import numpy as np -import matplotlib as mpl -mpl.use('tkagg') +import matplotlib +import matplotlib.pyplot as plt +matplotlib.use('tkagg') -sys.path.append('../') -# testing the the linter g;l +import sys +sys.path.append("../") +sys.path.append("../gMLV") +from gMLV import * # work around for retracing warning logging.disable(logging.WARNING) @@ -50,9 +45,245 @@ def set_all_seeds(seed): np.random.seed(seed) random.seed(seed) -# some plotting functions +def generate_params(num_species, num_pert, zero_prop=0, hetergeneous=False): + ''' + generates parameters for GLV simulation according to Cao et al 2017 + (Inferring human microbial dynamics from temporal metagenomics data: Pitfalls and lessons) + Method in the supplimentary + num_species: number of microbial strains + num_perterbations: number of perterbations + zero_prop: proportion of the interaction matrix that should be zeros + ''' + + N = numpy.random.normal(0, 1, (num_species, num_species)) + + if hetergeneous: + y = 1.2 + u = numpy.random.uniform(0, 1, size=(num_species)) + H = (1-u)**(1/(1-y)) + H = numpy.diag(H) + s = numpy.sum(H) + else: + H = numpy.eye(num_species) + # s = 3 #from the paper + s = numpy.sum(H) # this seems to prevent instability when more species + + a = numpy.random.binomial(1, 1-zero_prop, size=(num_species, num_species)) + # the interaction matrix + A = 1/s*N@H*a + + # set all diagonal elements to -1 to ensure stability + numpy.fill_diagonal(A, -1) + # generate feasible growth rate + r = numpy.random.uniform(0.00001, 1, size=(num_species)) + ss = -numpy.linalg.inv(A)@r + + while not numpy.all(ss >= 0): + + # changed max from 1 to 0.5 for stability of binary perts with few species + r = numpy.random.uniform(0.00001, 1., size=(num_species)) + ss = -numpy.linalg.inv(A) @ r + + C = numpy.random.uniform(-3, 3, size=(num_species, num_pert)) * 1/s + + # for the binary pert scheme choose ICs to be close to the ss + ICs = ss # this can be changed to start slightly away from ss + return r, A, C, ICs + +def generate_data_perts(simulator, tmax, sampling_time, dt, num_timecourses, ICs, num_pert, species_prob=1, num_metabolites=0, noise_std=0): + '''' + Generates data with external perturbations e.g. antibiotics or food. + + simulator: simulator object of the gMLV_sim class above + tmax: max time (days) + sampling_time: time between different perturbations + dt: time between different simulated points + num_timecourses:number of time courses to simulate + ICs: intial conditions + num_pert: number of different perturbations + species_prob: probability of each species appearing in each timecourse + num_metabolites: number of metabolites + noise_std: standard dev of measruement noise + ''' + + ryobs = [] # species + rsobs = [] # metabolites + rysim = [] + rssim = [] + ry0 = [] + rs0 = [] + all_perts = [] + + times = numpy.arange(0, tmax, dt) + + num_species = simulator.nsp + + for timecourse_idx in range(num_timecourses): + + pert_matrix = numpy.random.binomial(1, 0.5, size=(tmax//sampling_time, num_pert)) + + #print( "perturbations: ") + #print(pert_matrix ) + + all_perts.append(pert_matrix) + + # initial conditions + init_species = numpy.random.uniform(low=0, high=2, size=( + num_species,)) * ICs * numpy.random.binomial(1, species_prob, size=(num_species,)) + init_metabolites = numpy.random.uniform( + low=10, high=50, size=num_metabolites) + + ysim, ssim, sy0, mu, M, _ = simulator.simulate(times=times, sy0=numpy.hstack((init_species, init_metabolites)), + u=lambda t: binary_step_pert(t, pert_matrix, sampling_time)) + if numpy.sum(ysim > 10) < 0: # instability + print('unstable') + else: + yobs = ysim[0:-1:int(sampling_time // dt)] + sobs = ssim[0:-1:int(sampling_time // dt)] + # add some gaussian noise + yobs = yobs + \ + numpy.random.normal(loc=0, scale=noise_std, size=yobs.shape) + sobs = sobs + \ + numpy.random.normal(loc=0, scale=noise_std, size=sobs.shape) + + # append results + ryobs.append(yobs) + rsobs.append(sobs) + rysim.append(ysim) + rssim.append(rssim) + + ry0.append(init_species) + rs0.append(init_metabolites) + # Xs, Fs = linearize_time_course_16S(yobs,times) + # X = numpy.vstack([X, Xs]) + # F = numpy.vstack([F, Fs]) + + ryobs = numpy.array(ryobs) + rysim = numpy.array(rysim) + all_perts = numpy.array(all_perts) + + return ryobs, rysim, all_perts + + +def generate_data_transplant(simulator, tmax, sampling_time, dt, num_timecourses, ICs, species_prob=1, num_metabolites=0, noise_std=0): + '''' + Generates data with transplant perturbations + + simulator: simulator object of the gMLV_sim class above + tmax: max time (days) + sampling_time: time between different perturbations + dt: time between different simulated points + num_timecourses:number of time courses to simulate + ICs: intial conditions + species_prob: probability of each species appearing in each timecourse + num_metabolites: number of metabolites + noise_std: standard dev of measruement noise + ''' + + ryobs = [] # species + rsobs = [] # metabolites + rysim = [] + rssim = [] + ry0 = [] + rs0 = [] + all_perts = [] + + times = numpy.arange(0, sampling_time, dt) + + num_species = simulator.nsp + + for timecourse_idx in range(num_timecourses): + + # initial conditions + init_species = numpy.random.uniform(low=0, high=2, size=( + 1, num_species)) * ICs * numpy.random.binomial(1, species_prob, size=(1, num_species)) + init_metabolites = numpy.random.uniform( + low=10, high=50, size=(1, num_metabolites)) + + ysim = [] + ssim = [] + + p_matrix = [] + ys = init_species + ss = init_metabolites + yobs = [ + ys[0] + numpy.random.normal(loc=0, scale=noise_std, size=ys[0].shape)] + sobs = [ + ss[0] + numpy.random.normal(loc=0, scale=noise_std, size=ss[0].shape)] + + p = numpy.zeros((num_species,)) + + perturbed = False + for i in range(int(tmax//sampling_time)): + + # print(yo.shape, ss.shape) + + ys, ss, sy0, mu, M, _ = simulator.simulate( + times=times, sy0=numpy.hstack((ys[-1, :], ss[-1, :]))) + + if numpy.random.uniform() < 0.1 and not perturbed and i < int(tmax//sampling_time)-1: + perturbed = True + + p_rem = numpy.random.uniform(low=0, high=1, size=(num_species,)) * numpy.random.binomial(1, species_prob, + size=( + num_species,)) + + p_add = numpy.random.uniform(low=0, high=1, size=(num_species,)) * numpy.random.binomial(1, species_prob, + size=( + num_species,)) + p = p_add - 2*p_rem + else: + p = numpy.zeros((num_species,)) + p_matrix.append(p) + + ys[-1, :] += p + ys[ys < 0] = 0 + + # print(yo.shape, ss.shape) + yo = ys[-1] + so = ss[-1] + # add some gaussian noise + + yo = yo + numpy.random.normal(loc=0, scale=noise_std, size=yo.shape) + so = so + numpy.random.normal(loc=0, scale=noise_std, size=so.shape) + + ysim.extend(ys) + ssim.extend(ss) + + if i < int(tmax//sampling_time)-1: + + yobs.append(yo) + sobs.append(so) + all_perts.append(p_matrix) + # append results + ryobs.append(yobs) + rsobs.append(sobs) + rysim.append(ysim) + rssim.append(rssim) + + ry0.append(init_species) + rs0.append(init_metabolites) + # Xs, Fs = linearize_time_course_16S(yobs,times) + # X = numpy.vstack([X, Xs]) + # F = numpy.vstack([F, Fs]) + + ryobs = numpy.array(ryobs) + rysim = numpy.array(rysim) + all_perts = numpy.array(all_perts) + + return ryobs, rysim, all_perts + +def binary_step_pert(t, pert_matrix, dt): + # solver sometimes goes slightly past end of time interval + i = min(int(t//dt), len(pert_matrix)-1) + + p = pert_matrix[i] + return p + + +# some plotting functions def plot_fit_gMLV_pert(yobs, yobs_h, perts, sobs, sobs_h, sampling_times, ysim, times): # plot the fit fig, axs = plt.subplots(1, 2, figsize=(16., 6.)) @@ -106,7 +337,7 @@ def plot_fit_gMLV_pert(yobs, yobs_h, perts, sobs, sobs_h, sampling_times, ysim, if __name__ == '__main__': - if len(sys.argv) == 3: + if len(sys.argv) == 4: # check if the third argument is a number, and if the second one is a path if not sys.argv[2].isdigit(): print("Please enter a valid number of simulations") @@ -117,18 +348,19 @@ def plot_fit_gMLV_pert(yobs, yobs_h, perts, sobs, sobs_h, sampling_times, ysim, num_sims = int(sys.argv[2]) save_path = sys.argv[1] + '/' + mode = int(sys.argv[3]) os.makedirs(save_path, exist_ok=True) else: print("Using default values for number of simulations and save path") - print("Usage: python run_gMLV_sims.py ") + print("Usage: python run_gMLV_sims.py ") num_sims = 100 save_path = 'outputs/' + mode = 0 os.makedirs(save_path, exist_ok=True) # set_all_seeds(0) # total number of time courses will be num_sims x num_timecourses (per timecourse) - num_timecourses = 1 # 9*100 num_species = 3 @@ -156,22 +388,24 @@ def plot_fit_gMLV_pert(yobs, yobs_h, perts, sobs, sobs_h, sampling_times, ysim, all_ryobs = np.zeros([num_sims, sampling_times.shape[0], num_species]) all_rysim = np.zeros([num_sims, times.shape[0], num_species]) - all_perts = np.zeros([num_sims, sampling_times.shape[0], num_pert]) - all_parms = np.zeros( - [num_sims, num_species + num_species*num_species + num_species]) + all_parms = np.zeros([num_sims, num_species + num_species*num_species + num_species]) + + if mode == 0: + # This is parameter perturbations + all_perts = np.zeros([num_sims, sampling_times.shape[0], num_pert]) + else: + # This is transplant perturbations + all_perts = np.zeros([num_sims, sampling_times.shape[0], num_species]) for nsim in range(num_sims): # print("nsim",nsim) - # QUESTION: what is the purpose of this if loop? - if nsim % 100 == 0: + if nsim % 10 == 0: print('percent data generated:', nsim/num_sims * 100) # generate params according to paper approach - # C is perturbation interaction vector/m - # TODO: #25 generate_params is not defined anywhere - mu, M, C, ss = generate_params( - num_species, num_pert, zero_prop=zero_prop, hetergeneous=False) + # C is perturbation interaction vector/m (also called epsilon) + mu, M, C, ss = generate_params(num_species, num_pert, zero_prop=zero_prop, hetergeneous=False) # print("mu: ", mu) # print("M: ", M) @@ -186,11 +420,13 @@ def plot_fit_gMLV_pert(yobs, yobs_h, perts, sobs, sobs_h, sampling_times, ysim, num_metabolites=num_metabolites, M=M, mu=mu, - C=C) + epsilon=C) + + if mode == 0: + ryobs, rysim, perts = generate_data_perts(simulator, tmax, sampling_time, dt, num_timecourses, ss, num_pert, species_prob=species_prob, noise_std=0.00) - # FIXME: #26 generate_data_perts is not defined anywhere - ryobs, rysim, perts = generate_data_perts( - simulator, tmax, sampling_time, dt, num_timecourses, ss, num_pert, species_prob=species_prob, noise_std=0.00) + else: + ryobs, rysim, perts = generate_data_transplant(simulator, tmax, sampling_time, dt, num_timecourses, ss, species_prob=1, num_metabolites=0, noise_std=0.00) # print(ryobs.shape, rysim.shape, all_perts.shape) diff --git a/gMLV/gMLV_sim.py b/gMLV/gMLV_sim.py index 363b387..596198c 100644 --- a/gMLV/gMLV_sim.py +++ b/gMLV/gMLV_sim.py @@ -61,8 +61,8 @@ def __init__(self, num_species=2, num_metabolites=0, num_perturbations=0, mu=Non else: self.epsilon = epsilon - def simulate(self, times, sy0, tp=None): - syobs = odeint(gMLV, sy0, times, args=(self.nsp, self.np, self.mu, self.M, self.beta, self.epsilon, tp)) + def simulate(self, times, sy0, u=None): + syobs = odeint(gMLV, sy0, times, args=(self.nsp, self.np, self.mu, self.M, self.beta, self.epsilon, u)) yobs = syobs[:, 0:self.nsp] sobs = syobs[:, self.nsp:] return yobs, sobs, sy0, self.mu, self.M, self.beta @@ -75,7 +75,7 @@ def print(self): print(f'perturbation matrix: \n{self.epsilon}') -def gMLV(sy, t, nsp, np, mu, M, beta, epsilon, tp): +def gMLV(sy, t, nsp, np, mu, M, beta, epsilon, u): """ generalised Lotka Volterra with metabolite production @@ -85,7 +85,7 @@ def gMLV(sy, t, nsp, np, mu, M, beta, epsilon, tp): :param mu: specific growth rates vector :param M: interaction matrix :param beta: metabolite production rate matrix - :param p: a tuple containing time-dependent perturbation and perturbation matrix + :param u: a function that returns the perturbation signal at time t :return: change in species + metabolites vector """ @@ -93,25 +93,20 @@ def gMLV(sy, t, nsp, np, mu, M, beta, epsilon, tp): y = sy[0:nsp] s = sy[nsp:] - if np > 0: - for p_i in range(np): - if tp[p_i][0] <= t < tp[p_i][1]: - instantaneous_growth = mu + M @ y + epsilon[:, p_i] - else: - instantaneous_growth = mu + M @ y - else: + #if np > 0: + # for p_i in range(np): + # if tp[p_i][0] <= t < tp[p_i][1]: + # instantaneous_growth = mu + M @ y + epsilon[:, p_i] + # else: + # instantaneous_growth = mu + M @ y + #else: + # instantaneous_growth = mu + M @ y + + if u is None: instantaneous_growth = mu + M @ y + else: + instantaneous_growth = mu + M @ y + epsilon @ u(t) - # if p[0] is None: - # instantaneous_growth = mu + M @ y - # # dN = np.multiply(mu, y) + np.multiply(y, M @ y) - # else: - # if p[0] <= t < (p[0] + 1): - # instantaneous_growth = mu + M @ y + p[1] - # # dN = np.multiply(mu, y) + np.multiply(y, M @ y) + np.multiply(y, p[1]) - # else: - # instantaneous_growth = mu + M @ y - # # dN = np.multiply(mu, y) + np.multiply(y, M @ y) dN = numpy.multiply(y, instantaneous_growth) if beta is None: @@ -129,3 +124,4 @@ def gMLV(sy, t, nsp, np, mu, M, beta, epsilon, tp): dS = q @ y return numpy.hstack((dN, dS)) + diff --git a/gMLV/gMLV_sim_merged.py b/gMLV/gMLV_sim_merged.py deleted file mode 100644 index a2ab6ba..0000000 --- a/gMLV/gMLV_sim_merged.py +++ /dev/null @@ -1,384 +0,0 @@ -import random -from scipy import stats -import numpy as np -from scipy.integrate import odeint - - -class gMLV_sim: - def __init__(self, num_species=2, num_metabolites=0, num_perturbations=0, mu=None, M=None, beta=None, epsilon=None, C=None): - self.nsp = num_species - self.nm = num_metabolites - self.np = num_perturbations - - if mu is None: - self.mu = np.random.lognormal(0.01, 0.5, self.nsp) - else: - self.mu = mu - - if M is None: - self.M = np.zeros((self.nsp, self.nsp)) - # add self repression on the diagonal - for species_idx in range(self.nsp): - self.M[species_idx, species_idx] = random.uniform(-0.5, 0.0) - - # add random interactions - for i in range(self.nsp): - for j in range(self.nsp): - if i == j: - continue - else: - tau = stats.halfcauchy.rvs(loc=0, scale=0.001) - lam = stats.halfcauchy.rvs(loc=0, scale=1) - M = stats.norm.rvs(loc=0, scale=tau*lam) - # if i == j: - # self.M[i, j] = -abs(M) - self.M[i, j] = M - # i = random.randint(0, self.nsp-1) - # j = random.randint(0, self.nsp-1) - # self.M[i, j] = random.normalvariate(mu=0, sigma=0.1) - else: - self.M = M - - if beta is None and self.nm > 0: - self.beta = np.zeros((self.nm, self.nsp)) - for _ in range(self.nm): - i = random.randint(0, self.nm-1) - j = random.randint(0, self.nsp-1) - self.beta[i, j] = random.uniform(a=0, b=1) - else: - self.beta = beta - - if epsilon is None: - self.epsilon = np.zeros((self.nsp, self.np)) - - # add random interactions - for i in range(self.nsp): - for j in range(self.np): - tau = stats.halfcauchy.rvs(loc=0, scale=1) - lam = stats.halfcauchy.rvs(loc=0, scale=1) - epsilon = stats.norm.rvs(loc=0, scale=tau * lam) - self.epsilon[i, j] = -abs(epsilon) - else: - self.epsilon = epsilon - self.C = C - - def simulate(self, times, sy0, tp=None): - syobs = odeint(gMLV, sy0, times, args=(self.nsp, self.np, - self.mu, self.M, self.beta, self.C, self.epsilon, tp)) - yobs = syobs[:, 0:self.nsp] - sobs = syobs[:, self.nsp:] - return yobs, sobs, sy0, self.mu, self.M, self.beta - - def print(self): - print(f'number of species: {self.nsp}') - print(f'specific growth rates: {self.mu}') - print(f'interaction matrix: \n{self.M}') - print(f'metabolite production: \n{self.beta}') - print(f'perturbation matrix: \n{self.epsilon}') - -# FIXME: This merge brought C and epsilon. Decide if both are necessary or not -def gMLV(sy, t, nsp, np, mu, M, beta, C, epsilon, tp): - """ - generalised Lotka Volterra with metabolite production - - :param sy: species + metabolites vector - :param t: time - :param nsp: number of species - :param mu: specific growth rates vector - :param M: interaction matrix - :param beta: metabolite production rate matrix - :param p: perturbation function that returns the perturbation vector as a function of time - :return: change in species + metabolites vector - """ - - # separate species and metabolites - - sy[sy < 0] = 0 - y = sy[0:nsp] - s = sy[nsp:] - - if np > 0: - for p_i in range(np): - if tp[p_i][0] <= t < tp[p_i][1]: - instantaneous_growth = mu + M @ y + epsilon[:, p_i] - else: - instantaneous_growth = mu + M @ y - else: - instantaneous_growth = mu + M @ y - - # if p[0] is None: - # instantaneous_growth = mu + M @ y - # # dN = np.multiply(mu, y) + np.multiply(y, M @ y) - # else: - # if p[0] <= t < (p[0] + 1): - # instantaneous_growth = mu + M @ y + p[1] - # # dN = np.multiply(mu, y) + np.multiply(y, M @ y) + np.multiply(y, p[1]) - # else: - # instantaneous_growth = mu + M @ y - # # dN = np.multiply(mu, y) + np.multiply(y, M @ y) - dN = np.multiply(y, instantaneous_growth) - - if beta is None: - dS = [] - else: - # this is simple production - # dS = beta @ y - - # metabolite production as in Clark et al., 2021: eqs(4 & 5) - if len(beta.shape) == 3: - rho = np.dot(beta, y) # eq(6) - else: - rho = beta - q = np.multiply(rho, instantaneous_growth) - dS = q @ y - - return np.hstack((dN, dS)) - - -def generate_params(num_species, num_pert, zero_prop=0, hetergeneous=False): - ''' - generates parameters for GLV simulation according to Cao et al 2017 - (Inferring human microbial dynamics from temporal metagenomics data: Pitfalls and lessons) - Method in the supplimentary - num_species: number of microbial strains - num_perterbations: number of perterbations - zero_prop: proportion of the interaction matrix that should be zeros - ''' - - N = np.random.normal(0, 1, (num_species, num_species)) - - if hetergeneous: - y = 1.2 - u = np.random.uniform(0, 1, size=(num_species)) - H = (1-u)**(1/(1-y)) - H = np.diag(H) - s = np.sum(H) - else: - H = np.eye(num_species) - # s = 3 #from the paper - s = np.sum(H) # this seems to prevent instability when more species - - a = np.random.binomial(1, 1-zero_prop, size=(num_species, num_species)) - # the interaction matrix - A = 1/s*N@H*a - - # set all diagonal elements to -1 to ensure stability - np.fill_diagonal(A, -1) - # generate feasible growth rate - r = np.random.uniform(0.00001, 1, size=(num_species)) - ss = -np.linalg.inv(A)@r - - while not np.all(ss >= 0): - - # changed max from 1 to 0.5 for stability of binary perts with few species - r = np.random.uniform(0.00001, 1., size=(num_species)) - ss = -np.linalg.inv(A) @ r - - C = np.random.uniform(-3, 3, size=(num_species, num_pert)) * 1/s - - # for the binary pert scheme choose ICs to be close to the ss - ICs = ss # this can be changed to start slightly away from ss - return r, A, C, ICs - - -def binary_step_pert(t, pert_matrix, dt): - # solver sometimes goes slightly past end of time interval - i = min(int(t//dt), len(pert_matrix)-1) - - p = pert_matrix[i] - return p - - -def generate_data_perts(simulator, tmax, sampling_time, dt, num_timecourses, ICs, num_pert, species_prob=1, num_metabolites=0, noise_std=0): - '''' - Generates data with external perturbations e.g. antibiotics or food. - - simulator: simulator object of the gMLV_sim class above - tmax: max time (days) - sampling_time: time between different perturbations - dt: time between different simulated points - num_timecourses:number of time courses to simulate - ICs: intial conditions - num_pert: number of different perturbations - species_prob: probability of each species appearing in each timecourse - num_metabolites: number of metabolites - noise_std: standard dev of measruement noise - ''' - - ryobs = [] # species - rsobs = [] # metabolites - rysim = [] - rssim = [] - ry0 = [] - rs0 = [] - all_perts = [] - - times = np.arange(0, tmax, dt) - - num_species = simulator.nsp - - for timecourse_idx in range(num_timecourses): - if timecourse_idx % 100 == 0: - print('percent data generated:', - timecourse_idx/num_timecourses * 100) - - pert_matrix = np.random.binomial( - 1, 0.5, size=(tmax//sampling_time, num_pert)) - - all_perts.append(pert_matrix) - - # initial conditions - init_species = np.random.uniform(low=0, high=2, size=( - num_species,)) * ICs * np.random.binomial(1, species_prob, size=(num_species,)) - init_metabolites = np.random.uniform( - low=10, high=50, size=num_metabolites) - - ysim, ssim, sy0, mu, M, _ = simulator.simulate(times=times, sy0=np.hstack((init_species, init_metabolites)), - p=lambda t: binary_step_pert(t, pert_matrix, sampling_time)) - if np.sum(ysim > 10) < 0: # instability - print('unstable') - else: - yobs = ysim[0:-1:int(sampling_time // dt)] - sobs = ssim[0:-1:int(sampling_time // dt)] - # add some gaussian noise - yobs = yobs + \ - np.random.normal(loc=0, scale=noise_std, size=yobs.shape) - sobs = sobs + \ - np.random.normal(loc=0, scale=noise_std, size=sobs.shape) - - # append results - ryobs.append(yobs) - rsobs.append(sobs) - rysim.append(ysim) - rssim.append(rssim) - - ry0.append(init_species) - rs0.append(init_metabolites) - # Xs, Fs = linearize_time_course_16S(yobs,times) - # X = np.vstack([X, Xs]) - # F = np.vstack([F, Fs]) - - ryobs = np.array(ryobs) - rysim = np.array(rysim) - all_perts = np.array(all_perts) - - return ryobs, rysim, all_perts - - -def generate_data_transplant(simulator, tmax, sampling_time, dt, num_timecourses, ICs, species_prob=1, num_metabolites=0, noise_std=0): - '''' - Generates data with transplant perturbations - - simulator: simulator object of the gMLV_sim class above - tmax: max time (days) - sampling_time: time between different perturbations - dt: time between different simulated points - num_timecourses:number of time courses to simulate - ICs: intial conditions - species_prob: probability of each species appearing in each timecourse - num_metabolites: number of metabolites - noise_std: standard dev of measruement noise - ''' - - ryobs = [] # species - rsobs = [] # metabolites - rysim = [] - rssim = [] - ry0 = [] - rs0 = [] - all_perts = [] - - times = np.arange(0, sampling_time, dt) - - num_species = simulator.nsp - - for timecourse_idx in range(num_timecourses): - - if timecourse_idx % 100 == 0: - print('percent data generated:', - timecourse_idx/num_timecourses * 100) - - # initial conditions - init_species = np.random.uniform(low=0, high=2, size=( - 1, num_species)) * ICs * np.random.binomial(1, species_prob, size=(1, num_species)) - init_metabolites = np.random.uniform( - low=10, high=50, size=(1, num_metabolites)) - - ysim = [] - ssim = [] - - p_matrix = [] - ys = init_species - ss = init_metabolites - yobs = [ - ys[0] + np.random.normal(loc=0, scale=noise_std, size=ys[0].shape)] - sobs = [ - ss[0] + np.random.normal(loc=0, scale=noise_std, size=ss[0].shape)] - - p = np.zeros((num_species,)) - - perturbed = False - for i in range(int(tmax//sampling_time)): - - # print(yo.shape, ss.shape) - - ys, ss, sy0, mu, M, _ = simulator.simulate( - times=times, sy0=np.hstack((ys[-1, :], ss[-1, :]))) - - if np.random.uniform() < 0.1 and not perturbed and i < int(tmax//sampling_time)-1: - perturbed = True - - p_rem = np.random.uniform(low=0, high=1, size=(num_species,)) * np.random.binomial(1, species_prob, - size=( - num_species,)) - - p_add = np.random.uniform(low=0, high=1, size=(num_species,)) * np.random.binomial(1, species_prob, - size=( - num_species,)) - p = p_add - 2*p_rem - else: - p = np.zeros((num_species,)) - p_matrix.append(p) - - ys[-1, :] += p - ys[ys < 0] = 0 - - # print(yo.shape, ss.shape) - yo = ys[-1] - so = ss[-1] - # add some gaussian noise - - yo = yo + np.random.normal(loc=0, scale=noise_std, size=yo.shape) - so = so + np.random.normal(loc=0, scale=noise_std, size=so.shape) - - ysim.extend(ys) - ssim.extend(ss) - - if i < int(tmax//sampling_time)-1: - - yobs.append(yo) - sobs.append(so) - - all_perts.append(p_matrix) - # append results - ryobs.append(yobs) - rsobs.append(sobs) - rysim.append(ysim) - rssim.append(rssim) - - ry0.append(init_species) - rs0.append(init_metabolites) - # Xs, Fs = linearize_time_course_16S(yobs,times) - # X = np.vstack([X, Xs]) - # F = np.vstack([F, Fs]) - - ryobs = np.array(ryobs) - rysim = np.array(rysim) - all_perts = np.array(all_perts) - - return ryobs, rysim, all_perts - - -def set_all_seeds(seed): - np.random.seed(seed) - random.seed(seed) diff --git a/gMLV/run_gMLV_sims.py b/gMLV/run_gMLV_sims.py deleted file mode 100644 index f343ced..0000000 --- a/gMLV/run_gMLV_sims.py +++ /dev/null @@ -1,225 +0,0 @@ -''' -Simulation script for gMLV model with perturbations and metabolites. -Usage: python run_gMLV_sims.py -Example: python run_gMLV_sims.py outputs/ 100 - -The script will create a folder named outputs/ and save the results there. -The number of simulations is 100 by default, but can be changed by passing -a different number as the second argument. -''' - - -import logging -from time import time -import math -from gMLV import * -from numpy import linalg as la -import copy -from scipy.integrate import odeint -import matplotlib.pyplot as plt -import random -import os -import sys - -import numpy as np -import matplotlib as mpl -mpl.use('tkagg') - - -sys.path.append('../') -# testing the the linter g;l - -# work around for retracing warning -logging.disable(logging.WARNING) -os.environ["TF_CPP_MIN_LOG_LEVEL"] = "3" - -SMALL_SIZE = 13 -MEDIUM_SIZE = 17 -BIGGER_SIZE = 20 - -plt.rc('font', size=SMALL_SIZE) # controls default text sizes -plt.rc('axes', titlesize=SMALL_SIZE) # fontsize of the axes title -plt.rc('axes', labelsize=MEDIUM_SIZE) # fontsize of the x and y labels -plt.rc('xtick', labelsize=SMALL_SIZE) # fontsize of the tick labels -plt.rc('ytick', labelsize=SMALL_SIZE) # fontsize of the tick labels -plt.rc('legend', fontsize=SMALL_SIZE) # legend fontsize -plt.rc('figure', titlesize=BIGGER_SIZE) # fontsize of the figure title - - -def set_all_seeds(seed): - np.random.seed(seed) - random.seed(seed) - -# some plotting functions - - -def plot_fit_gMLV_pert(yobs, yobs_h, perts, sobs, sobs_h, sampling_times, ysim, times): - # plot the fit - fig, axs = plt.subplots(1, 2, figsize=(16., 6.)) - - for species_idx in range(yobs.shape[1]): - axs[0].plot(times, ysim[:, species_idx], '--', label='simulation') - - axs[0].set_prop_cycle(None) - - for species_idx in range(yobs.shape[1]): - axs[0].scatter(sampling_times, yobs[:, species_idx], - s=100, marker='x', label='observed') - - axs[0].set_prop_cycle(None) - - # for species_idx in range(yobs.shape[1]): - # axs[0].scatter(sampling_times, yobs_h[:, species_idx], s= 100,marker ='x', label = 'prediction') - - axs[0].set_xlabel('time (days)') - axs[0].set_ylabel('[species]') - - handles, labels = axs[0].get_legend_handles_labels() - newLabels, newHandles = [], [] - for handle, label in zip(handles, labels): - if label not in newLabels: - newLabels.append(label) - newHandles.append(handle) - - axs[0].legend(newHandles, newLabels) - - axs[1].set_prop_cycle(None) - # perts = np.vstack((perts[0], perts[0], perts)) - # sampling_times = np.append(sampling_times, 100) - - for pert_idx in range(perts.shape[1]): - axs[1].scatter(sampling_times[1:], - perts[:, pert_idx], marker='o', s=100) - axs[1].set_xlim(left=0, right=100) - - axs[1].set_ylabel('transplant perturbation') - axs[1].set_xlabel('time') - - # for metabolite_idx in range(sobs.shape[1]): - # axs[1].plot(timepoints, sobs[:, metabolite_idx], color=cols[metabolite_idx]) - # axs[1].plot(timepoints, sobs_h[:, metabolite_idx], '--', color=cols[metabolite_idx]) - # axs[1].set_xlabel('time') - # axs[1].set_ylabel('[metabolite]'); - -# set_all_seeds(1234) - - -if __name__ == '__main__': - - if len(sys.argv) == 3: - # check if the third argument is a number, and if the second one is a path - if not sys.argv[2].isdigit(): - print("Please enter a valid number of simulations") - sys.exit(1) - if not os.path.isdir(sys.argv[1]): - print("Please enter a valid path to save the outputs") - sys.exit(1) - - num_sims = int(sys.argv[2]) - save_path = sys.argv[1] + '/' - os.makedirs(save_path, exist_ok=True) - else: - print("Using default values for number of simulations and save path") - print("Usage: python run_gMLV_sims.py ") - num_sims = 100 - save_path = 'outputs/' - os.makedirs(save_path, exist_ok=True) - - # set_all_seeds(0) - - # total number of time courses will be num_sims x num_timecourses (per timecourse) - - num_timecourses = 1 # 9*100 - - num_species = 3 - - # controls probability of dropout - species_prob = 1.0 - - # npert is number of independent perturbations - # FIXME: change num_pert to 0 and see fix issue with input array sizes or shapes - num_pert = 1 - num_metabolites = 0 - - # construct interaction matrix - zero_prop = 0. # the proportion of zeros in the interaction matrix - - tmax = 100 - sampling_time = 10 - dt = 1 - - times = np.arange(0, tmax, dt) - sampling_times = np.arange(0, tmax, sampling_time) - - # print("npert", num_pert) - # print("nsims", num_timecourses) - - all_ryobs = np.zeros([num_sims, sampling_times.shape[0], num_species]) - all_rysim = np.zeros([num_sims, times.shape[0], num_species]) - all_perts = np.zeros([num_sims, sampling_times.shape[0], num_pert]) - all_parms = np.zeros( - [num_sims, num_species + num_species*num_species + num_species]) - - for nsim in range(num_sims): - # print("nsim",nsim) - - # QUESTION: what is the purpose of this if loop? - if nsim % 100 == 0: - print('percent data generated:', nsim/num_sims * 100) - - # generate params according to paper approach - # C is perturbation interaction vector/m - # TODO: #25 generate_params is not defined anywhere - mu, M, C, ss = generate_params( - num_species, num_pert, zero_prop=zero_prop, hetergeneous=False) - - # print("mu: ", mu) - # print("M: ", M) - # print("C: ", C) - - all_parms[nsim, :] = np.concatenate( - (mu.flatten(), M.flatten(), C.flatten()), axis=None) - # print("p:", all_parms[nsim,:] ) - - # instantiate simulator - simulator = gMLV_sim(num_species=num_species, - num_metabolites=num_metabolites, - M=M, - mu=mu, - C=C) - - # FIXME: #26 generate_data_perts is not defined anywhere - ryobs, rysim, perts = generate_data_perts( - simulator, tmax, sampling_time, dt, num_timecourses, ss, num_pert, species_prob=species_prob, noise_std=0.00) - - # print(ryobs.shape, rysim.shape, all_perts.shape) - - # species levels and perturbations for each time point - all_ryobs[nsim, :, :] = ryobs.astype(np.float32) - all_rysim[nsim, :, :] = rysim.astype(np.float32) - # export each simulation as csv - # create a numpy array concatenating the time points and the simulated data - data_export = np.concatenate( - (times.reshape(-1, 1), rysim[0, :, :]), axis=1) - np.savetxt(save_path + '/simulations' + str(nsim) + - '.csv', data_export, delimiter=',') - - # np.savetxt(save_path + '/simulations' + str(nsim) + '.csv', rysim, delimiter=',') - # np.savetxt(save_path + '/simulations.csv', rysim[0,:,:], delimiter=',') - all_perts[nsim, :, :] = perts.astype(np.float32) - - np.save(save_path + '/abundances_sampled.npy', all_ryobs) - np.save(save_path + '/abundances.npy', all_rysim) - np.save(save_path + '/perts.npy', all_perts) - np.save(save_path + '/parms.npy', all_parms) - - # plot some of the results - for i in range(10): - plot_fit_gMLV_pert(all_ryobs[i], 0, # pred[-i-1, :, :], - all_perts[i, 0:-1, :], None, None, sampling_times, all_rysim[i], times) - - # print("new timecourse") - # print( all_ryobs[i] ) - # print( all_perts[i, 0:-1, :] ) - - plt.savefig(save_path + '/test_plot_' + str(i) + '.pdf')