Skip to content

Instantly share code, notes, and snippets.

@jakeyeung

jakeyeung/sfi_playground_multi_traj_example_cat_inputs.ipynb

Created October 24, 2023 08:10
Show Gist options
  • Save jakeyeung/ecb5c05d405f109be198857a1c54c025 to your computer and use it in GitHub Desktop.
Save jakeyeung/ecb5c05d405f109be198857a1c54c025 to your computer and use it in GitHub Desktop.
Download ZIP
Concatenating input trajectories gives error
Display the source blob
Display the rendered blob
Raw
sfi_playground_multi_traj_example_cat_inputs.ipynb
{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"id": "68b3068f-3edf-4e90-9a82-8495da628743",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"[<matplotlib.lines.Line2D at 0x7f1c38185080>]"
]
},
"execution_count": 1,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAD4CAYAAAD1jb0+AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAuzUlEQVR4nO3dd3gU1frA8e+bDiEBQiCEUAJIkV4igig2EAQEVOwFvCp2xQ7qvRcLiuX606tYUFDsYgFRiiJciopgkN47JIQEAgkhkLrn98cs6QkJyWZ2N+/nefLMzJnZ3XcgeXNy5hQxxqCUUsrz+NgdgFJKqTOjCVwppTyUJnCllPJQmsCVUspDaQJXSikP5VedHxYeHm6io6Or8yOVUsrjrVq16rAxpmHR8mpN4NHR0cTGxlbnRyqllMcTkb0llWsTilJKeShN4Eop5aE0gSullIfSBK6UUh5KE7hSSnkoTeBKKeWhNIErpZSH8uwEvuRV2DoPcrPtjkQppapdtQ7kqVJZ6fDXB3A8EUIi4b4VEFTX7qiUUqraeG4NPCAYHt4I/Z+FtAR4qyfsXmZ3VEopVW08N4ED+PpD34fg2k8hMw0WT7I7IqWUqjaencABRKDDMOh1J+z/02paUUqpGsDzE/gp4W3BkQMr3rM7EqWUqhbek8A7XwutL4Gl/4Eju+yORimlXM57Erh/EFz8NBgHTO4NPz0CGal2R6WUUi7jPQkcoGkMPLAK2l4GsVPhq5vsjkgppVzGuxI4QN0ouO4za3/PMojTBSSUUt7J+xL4KXf/bm0/HgLG2BuLUkq5wGkTuIi0E5E1Bb6OichYEQkTkQUist25rV8dAZdb407Q5XrIyYCZd0FOlt0RKaVUlTptAjfGbDXGdDPGdAN6AieAmcA4YKExpg2w0HnsXga/Auc/DOu+hh8fgsM77I5IKaWqTEWbUC4Fdhpj9gLDgenO8unAiCqMq2oE1YVz77b2134B718Ah7baG5NSSlWRiibw64EvnfsRxpgEAOe2UUkvEJExIhIrIrGHDh0680jPVEhjuOk76HErZJ+AqZfB+m+rPw6llKpiYsr5gE9EAoADQEdjTKKIpBhj6hU4f9QYU2Y7eExMjImNtbFXyNE9MGMUJKyBAc9Z86gopZSbE5FVxpiYouUVqYFfDvxtjEl0HieKSKTzzSOBpMqH6WL1o+HWWdb+gn/ZGYlSSlVaRRL4DeQ3nwDMBkY590cBP1RVUC5Vqz60GWjtnzhibyxKKVUJ5UrgIlIbGAB8X6B4EjBARLY7z3nOXK4dr7S22xdoH3GllMcq14o8xpgTQIMiZclYvVI8T5sBUCcCZo6BXYvhynftjkgppSrMe0diliU4HO5bCefcaXUvXPqq3REppVSF1cwEDlCrHgx41tpf/LL2D1dKeZyam8DBWlfz9gUgPvDBpbBiCmSdsDsqpZQql5qdwAGa9YJ7l1vzic97HP7+xO6IlFKqXDSBAzRoDZe/Yu3Pf1K7FyqlPIIm8FM6XgltL7f2t83X7oVKKbenCfwUkfzuhLPu0cWRlVJuTxN4QbXqw13LrP357jc7rlJKFaQJvKjILtYgHwCHw95YlFKqDJrAS3L2MGv71weQm21vLEopVQpN4CVpO8jaznsCJp8LR/faG49SSpVAE3hJ2vSHu5Zac4anH4aPh0JOpt1RKaVUIZrASxPZ1VrwofNISN0HC/5td0RKKVWIJvDTaX2xtV3xLmyZa28sSilVgCbw0zn7CqtrYURn+OoGeLMrfPsPyEq3OzKlVA2nCbw8IrvAnQshqqe1ruaG7+C1djpaUyllq/KuyFNPRL4VkS0isllE+ohImIgsEJHtzm2ZCxp7PL9Aa+bCcfshoA5kpcGf79gdlVKqBitvDfxNYL4xpj3QFdgMjAMWGmPaAAudx97NxxeCQuGx7dCsN/z8lC4GoZSyzWkTuIiEAv2AqQDGmCxjTAowHJjuvGw6MMI1IbqhgNoweg50GAGLXoCvboLsDLujUkrVMOWpgbcCDgEfichqEflQRIKBCGNMAoBz26ikF4vIGBGJFZHYQ4cOVVngtvP1g6s+sPa3/ARLX7E3HqVUjVOeBO4H9ADeNcZ0B9KpQHOJMWaKMSbGGBPTsGHDMwzTTfkFwAN/W/vL/gPHEuyNRylVo5QngccBccaYFc7jb7ESeqKIRAI4t0muCdHNNWgNFz9j7c++X7sXKqWqzWkTuDHmILBfRNo5iy4FNgGzgVHOslHADy6J0BNc+Dj0HQs7foVFEyH7pN0RKaVqgPL2QnkA+FxE1gHdgBeBScAAEdkODHAe11y974XQpvDnZHj3PMg8bndESikvV64EboxZ42zH7mKMGWGMOWqMSTbGXGqMaePc1uyFJEMi4N4/oNcYOLIL5jxqd0RKKS+nIzGrUlBdGPwq9LgV1n1lNakopZSLaAJ3hf7PgvjA4km6qo9SymU0gbtC7TC45J8Q9xcs+Kfd0SilvJSf3QF4rfMfhpR9sPxtCI2CPvfaHZFSystoDdxVROCy5639n8fD/pX2xqOU8jqawF0pMASueBNqh8OMW7V/uFKqSmkTiqv1HA21G8DXN8NHl8MtM6GWd8+8q5TX+OE+WP2ZtX/RePANsJZbbH2J9Vd2SY4dgNkPQlQP6HIdhDYB/1ouCU9MNS5KEBMTY2JjY6vt89yGIxeeC8s/npBqXyxKqfKbULf0c2OWQJNuxctfawvHEwuX9Z9gPRc7QyKyyhgTU7Rcm1Cqg48vPLgm//iH++BEzR73pJTHW/KytS3aVbho8gb4dYJLVvDSBF5dwlpaMxe2ON/6k+zV1nBoq91RKaVKkrK/eO3bLwiu+6xw2YS68Fx9azuhLr//7yfS+k0o+T2Td1Z5mJrAq1OD1nDbHBg+GYwDPhmh62oq5Y7e6FS87MHVkHUi/zjm9mKX9F1yE3GLpxUuvOwFGB8H4WdVcZCawO3R5Tprm3YAFj5nbyxKqfLZvxJWf5p//PsbJV52ts++wgV97rd6pLmAJnA7+PrDiHet/d9eh3UztCauVHVy5ELCWlj8srUc4kdD4I+385pC6Hxt8dd8Mwr2LMs/LrhflhzXLbeovVDslLIfpg2CY3HQuDPc/ZvdESlVM5TVu6QKjM26lz2mMbP+dWuVdBvWXijuqF4zuHORtX9kj62hKFVjVEOz5SxHX9aYs1w+5kMTuN1CIiCkCWSlwd4/7I5GKe+37D/V8CHC1hcGufxTNIG7g1O18I8u166FSrlauyHV8jEz/453+WeUK4GLyB4RWS8ia0Qk1lkWJiILRGS7c6vjw89UaCSM/MjaX/Bve2NRypvNexK2znHZ2/fKmEx0xhcAjPt+vcs+55SK1MAvNsZ0K9CQPg5YaIxpAyx0Hqsz1ekqaDMQts2DjTPtjkYp77TivUKHuaHNSrws3jSo8FtHZ3xBEvn12EA/1zdwVOYThgPTnfvTgRGVjqamO+cOa/vNaFj1sdW9SSnlMr7H9hcrSzahDM58ib4Zb7Igtwd9M9487fucumbzc4NoHlabpvVrseV592kDN8AvIrJKRMY4yyKMMQkAzm2jkl4oImNEJFZEYg8dOlT5iL1Z28tg1I/W/o8PWX3ElVJV57HtZZ6eaobRM/M9UqlDPA25M/sx4mnI89k3513TLuNjbs4an3f8XPYtxNOQsf3bUCvAl6VPXMxvT16ClDZbYRUq73SyfY0xB0SkEbBARLaU9wOMMVOAKWD1Az+DGGuWlv3g6qnw3e1wuOxvNqVUBQU3LPN0tGM/fuSQUyQ19vNZl7e/NWh0oXOzcvsC8Mav2/H39eHamGbUCfTDxwcC/XyrJu5SlKsGbow54NwmATOBXkCiiEQCOLdJrgqyxonqYc0hvvH7/LmIlVKVZ8peZPxS39U0keS84xBO8Ib/21zou67U1xwhNG//1Z+3cs7EXzn7X/N5/JvSX1NVTpvARSRYREJO7QOXARuA2cAo52WjgB9cFWSNE9YKxm6ARh2sqWePHbA7IqW8Q8F5+UvxoF9+J4L1QXcwwrfs8Rk3+y4osXz2Wtf/3JanBh4B/CYia4GVwBxjzHxgEjBARLYDA5zHqqoE1Ia2A6396VdAenLZ1yulqsRI36VYj/3Kp7fPphLLp/+jVxVFVLrTtoEbY3YBXUsoTwYudUVQyumip6w/+f58D15vD1E9rZU9mve2OzKlPNO9K+Cdc097mT+59PNZW663bCjWCltv39idwZ0i2ZqYRqOQQBrUCaxUqOWhIzHdmV8ADHgO7lhgzT28bzlMGwhHdtsdmVKeqVF7uPL9Mi9ZktuFALKZGlC+IffnNq3NnklDGNqlCT4+wtmRodWSvEEXNfYMkV2hcRdw5MBfH8DnI+GBVXZHpZTnObgBZt5V6umeGe+STF26yY7yvV94W7hjYRUFV3GawD2FCAx5zUrgoU3sjkYpz5RcdmJeFXRPxd7vrmXgY19DhjaheJrQprB7KXzYHxJc301JKY/icMCLTeHZ+pDiXBnn/X6w5FX4T3trUYYCFvX9nI/qPUBqx5tLeLMCnijSbHnjNzAhFfyDqjD4itMauKe57lNY/jZs+A6mDoALn4Seo6H26btHKeX1EjdYUzMDvNE5vzyh5AeSl7RvBANeAOCu1eG8H/BG8YvuWGj9fE1ItX5BHN1trW/rBrQG7mmiesDIaVYNwL82LHzWqmEc2mZ3ZErZKzcHIjpW7DVTB8DxQ2yIT+VnRy+iMz4vfk3TAgvh+Pi4TfIGTeCeq+1l8MQuay7xrONW75QTR+yOSinX2PEr/F8naym0GYWbQfjfS1b58w1g0QvwwN8Ve+/gcIa+dWo5Q+G6zH/mnxv8WqXCdjVdE9MbxK+CDy6x9oe+ATG32RqOUlVq+6/w+dWFy0bPgejzrf1KrG/puGcF135/hNi9RwuVfzWmN71bVXxKWVcpbU1MbQP3BlE94ZJnrNrHT2OtofitLrQ7KqUq79DW4skb4OPKr6qzatAPXP1/O0s8507JuyzahOIt+j0OD2+COhHwyXCrVvLZ1ZCRandkSp0ZY2BygeHo4+PK/dJpOYPINGXXT6+elX6mkbkNTeDepG4U3LM8f6j9jl9hUnNdZ1N5DmPg709g6mXwbL388pEfQWAITEgl66J/lfkWv+d25LmcW2mX+QmH716P6TG60Pl3coblLXsGcGX3KLY8P4izI0O5q18r3r+lZxXekGtpE4q3CW4A/5hv7Z9qGzx2ABq2sy8mpcorZR/MfqBwWbebrSUHndrOb4/wGbuD8vtuF0zIBYU3bg7D3rS+nOqt2AczrfUqt74wKG/O7nkPXVBVd1FtNIF7s4Evws9PwewH4ba5UK/k9f+UchtFxzM8kwR+1rwiDofhvi+sHiYGH1pnfMptvvP5MvcS/jW0A3/uSmb1/hT+ero/SWkZNAgueT6SG89tTu9WYWw8cMzlCy64mvZC8WaO3MLzH1/+itXVsNsNUD/atrCUKtOJI/BKy2LFCfftpM9/VhQrX/3PAdQPDqiOyGxTWi8UTeA1wR9vwS/P5B/XbwniA0ecT+AveAw6XQ0RHeyJT6mifn7aGnFc0BX/5e3U83jtl220Cg9mwSMX4uvj+nUn3YF2I6zJ+txvTYCVk2k9JJo/HjIL9E5Z9pr1BfDYDqhT9rqBSrmcr3/xsh8fJGr4RgD6d4ioMcm7LOWugYuILxALxBtjhopIGPA1EA3sAa41xhwt/R20Bu42MtPgeJLVX/zobvh+DMT9lX8++gK48WsICLYvRqVOKfAXpDXUXZh1X1+6Natna1jVqbQaeEW6ET4EbC5wPA5YaIxpAyx0HitPEBhizecgYiXxO36FW2bln9+zDF5sAn9/Ckf32hamUnxxfaHmvz1BNwEwb32CXRG5lXI1oYhIU2AIMBF4xFk8HLjIuT8dWAw8WbXhqWrT+mJrDomTR2HTLKvWM/t+61xoUzgWB+IL3W6EvX9A485w7t3Qoo+tYSvvltFqAEHb5hUrv7pnUxuicT/lakIRkW+Bl4AQ4DFnE0qKMaZegWuOGmPql/DaMcAYgObNm/fcu1drdB7h5FHYsRCOJ1rbnaWsOnL/Kgg/q3pjU14vIfUkWxLS+GLlPhZsSuSShmlMS7NW0km+ZyMNImpWAj/jh5giMhRIMsasEpGLKvrBxpgpwBSw2sAr+nplk1r1ofNIa7/Pfdb26F4IiXQm9wXww30wfag1hN/GVUmU9/nnrA38ujmJALIBf6ZdUR+cY3UahEfYGps7KU8TSl9gmIgMBoKAUBH5DEgUkUhjTIKIRAJJrgxUuYH6LaxtSAR0vxkWTYS0A7D6E2tRCaWqSMusbewJuje/oOBAS+Oo9njc1WmrTcaY8caYpsaYaOB6YJEx5mZgNnBqYt5RwA8ui1K5p2s+sra7FtsahvJg236G/3aH5AKzAm5fwNPx95Z8/b0r8kZmqsr1A58EzBCR24F9wDVVE5LyGH++Y203zoSrPii5765SBWWdgLQE8PGFNV/Akpet8rd6lP26O/9nrUalCqlQAjfGLMbqbYIxJhm4tOpDUh5j8H+smlPiBusH8J7lEFjH7qiUuzIGXows9+VtMj4hGz9eHdmFa6J0Hp+S6JMndebqNIS7lln7Kfvgl6fhwBrIOGZrWMpNJW44/TU3fsP2od8RnfEF2c765YjuUS4OzHNpAleV4+MDdzvXE1z1MUy5ECY1gxVTbA1LuaGITsXLzhpgrSYFcNFTZLXqz1U/5uad/nP8pfj7apoqjc6FoiqvcWcY9pY1j7N/MGSnw7zHoUl3aHaO3dEpu/35nvW8pFMJS6Pd/K217fc4ALNXxZGWmQPA8vGX0LhuUHVF6ZH0V5uqGj1uhX+nwNMH4LrPrLKp/a0FaVXNlZsN85+ElL3w2+tW2aBJMG4fTCi83J/DYdielJZ3/POGg9UZqUfSGriqOuKcHe7sKyCyGySsyV+Q9o6F1hzkweE2Bads8dPD+fsjP7KmMe44othlO5KO0//1JYXKbu7dwsXBeT5N4Mo1xiyGRS/kT1P7YZEOS0/sLr76inIfxuT/Qi7JsQPwwaXWL+vBr5R8zdE9sPpTa3/gi4WWRcv/GMOG+GM8+s2avLJfHu5HcKAfftr2fVqawJVriMCl/7RmPZx1j1XW9UZY6xxS90pL6DUGDm8HR441Z3mT7tYQ/iO7IDfLalsvK4ko11g3w3qekZMB130Omcfy/w9vmWVNRzzjFut45fvWaNxTzWanJO8s3Lc7quSFgvtOWsSB1Iy845BAP9pGhFThzXg3XZFHVb8VU6wBHCcOAwKU8j0ovtDucqvmVtIDMFX1krbAO+dW/HWdroaR0/KPs04U7vP9wN/WL/MCjDG0HD8373jSVZ3p3rw+7RprAi9KV+RR7uPcMdDrTqvvuMmFlR/Cn5Pzz4e3tWrhjhzY8pP1tfknaHkBdL0B/GvZF7u32/tb8TL/2pB9onBZi/Phtjnwegc4Fg8bvoOrp0L6Ydj1P2v1p4KKJG+A487eJgDf3t2HmGhtUqsorYEr95F9EvyCrGYTR641lW3SJvisSO379gXQrFfJ7+FwQFYaBNV1fbzeJDsDfn8TFr9oHYe1gvtjrSHvp5w4YrVpt7wQmnSzyiaU49/5tvklzht/x/RYft2cCMDulwYj2lxWKq2BK/dXsGbt42ut4xnaxEokaQkw/Qrr3NQB+V3Q4v+Gw9ug/RBrpaEXI62220uegQ4jILxNtd+GW3I4rF+MRZOkMfBsvcJlV0/Nn0q4oNph0PehwmXXfprfHl5UCc0mp6SezM5L3gC5DoOfrybwitIauPIsp2p8rS+BnYtOf337odYDtppYuyv6IBHg7t+hwVngH2T10X6+SLfOh9Za3T0rInGj9X8R1hrqRsGR3dDmMgioXepLRk1byZJth/KO/3q6Pw1DdJbB0lTFmphK2e+ORVCnMSTvKFwe1ip/v+WFMOJda3/LT/BydLWF5zYcuSXP8PdeX5gYYdXIf3uj+Pmt8yv+WREd4bwHoP1giOxq9fMukrzXx6VyLCOb45k5vDx/Cxe2bVjo/LM/bqz45yptQlEepmlPeGyrte/IhdS4/IUmDqyGhLX5i0uEtYJpAyEjBd49H9oNgvMftnq3+DuHaP86AX77P2t/+GRroQpPlXYQPhtpDZYqOGCqaS+rSWTeE/llzxVZ/TCyK/S6C7pcV+VhrYtLYdjbv9P/7Eb8urnkdV/+2Jlc5Z9bE2gTivJu5XnIVtBTCWX+6V/tDm2FuY/D7iXWL6QHV5d8XU4mvNCoePn4OOvZAEDsR/DT2OLXNDsXbv+lykI+kp5FkL8Ph9OyqBfsz40f/MmG+NJnqHzpqs7c0Kt5lX2+N9KHmKpmOvWwc8G/rF4Wp/NuH6sd2G6JG61h5+/0zi87sqv4dVnpEBAMk0volXPvivzkDRBzm7VN3gEdr4IGrayBU1Vkzf4UHp2xhp2H0oud69K0LuviUouVfzT6HC5uX8IvHlUumsBVzTDgOeh8jVWLPbzdWpbrZIq1mtDAifD1LbBtnjX8+7/doV5za5j4OXdUb5y52VZT0NQBZV+35BX438SSz10zvcT5RoD8JF7Fdh46zojJv5d6/rpzmhVK4Esfv5hmYbW062AlnbYJRUSCgKVAIFbC/9YY828RCQO+BqKBPcC1xpijZb2XNqEot5VxzJrHvCQPb7IW0v3xQauHxuWvWMvHZRyzkuiWuXDfiqppeplyMRz4u/Tz//jZatcvzei5EN238nGU0/HMHP5vwTam/rY7r2zKLT35Y2cyj1zWll82JlLL35fLOzUmPSuH3YfT6dK0XrXF5y1Ka0IpTwIXINgYc1xE/IHfgIeAq4AjxphJIjIOqG+MebKs99IErtza8SR4rYR+44Nfs5pfUvfnl3W7CdZ8nn98x0JoGgPHEqz+0gUX3jXGmk/k+CEIP6vkz969NL+fe3k17QU3fm21k6fuh04jrQU2qtHV7/7Bqr1Wve3VkV0Y2qUJtQJ8T/MqVVFnnMCLvEltrAR+D/AJcJExJkFEIoHFxph2Zb1eE7jyKK+1tUaDnon7V1nNMC80LH7utvkQ2QVebAIN2kDy9uLXXDbRmjqgUUdrZOmHA/KvqxVmzfZ4qveNjaLHzcnb19GUrlOpBC4ivsAq4CxgsjHmSRFJMcbUK3DNUWNMsSciIjIGGAPQvHnznnv37j3zu1CqOm2ZC1/dULjsmUOw6Qf43tk23v0Wa16X9/tV/vMueBQufsZqrvF138dTuQ7DiMm/sz6+8EPJuQ9eQIcmoTZF5d2qqgZeD5gJPAD8Vp4EXpDWwJXHSU+2prYNCrV6e5ySGm+N7gxtYh07cq35z0+tOlNQ9AUw+idrP2EdvH9B/jn/YOg5Cuq1gHPvcusRo1k5Dg6mZvCv2RtYvPVQoXMXtAln2uhzdP1KF6mSboTGmBQRWQwMAhJFJLJAE0rJPfSV8mTBDUour1tkpXQfX+j/byvZb5oNY9dZ5cYUbpeO7GLV4tfPgA7DC3fzc1PJxzPp+ULZS+N9evsZTEGrKu20CVxEGgLZzuRdC+gPvAzMBkYBk5zbH1wZqFIeYeBE6+uUkmrUfgG2jfiM3XOEke8tB6B94xBeGdmFHIehe7N6xdqv31+yk7VxKZzbspRfYsp25amBRwLTne3gPsAMY8xPIrIcmCEitwP7gGtcGKdSqhIOH89kR9Jxrp/yZ17ZloNpDHs7v+/2Ozf1oHNUXXIchnVxKbw0bwsAc9fnLy58Q6/mZOc6OJ6RQ5dmdXll/lb+c03X6rsRVYgOpVfKC62PS2XuhgSGdW3C+rhUnvhuXd6562KaMbhLJCknsgj08+Huz0rud960fi3ijp4sVLZn0pBCxzm5Dl27shroUHqlapAr3rZW1nl38c5i514e2aXQ8Xf39GHke8spWJdrFR7MvLEX4O9sv1+9/yhR9YoPVNLkbS9N4Ep5mcyc3ELHCx7uR0TdIJKOZdKiQfEk3LNFGLtfGlKsvOg1yv1oAlfKy7y+YFve/sQrO9HGucp7aJC/XSEpF9EErpQXSTqWwZKth6gT6MdfT/fXYe1eThO4Ul7A4TC0emouAL4+wrPDOmryrgE0gSvl4bYnpvGP6X/lHX80+hz6tS1hDhbldTSBK+XB/tyVnNe3289HWPbkxUTWrWVzVKq6aAJXykMdz8zJS94DO0bw9o09dC6SGkYTuFIeJCM7l3kbEnj468LLvk28srMm7xpIE7hSHuSfszbwzaq4vOOZ955H9+ZVt66l8iyawJVyc8YY1sen5s1b0io8mJn39iW0lp8uoFDDaQJXys29OHczHyzLX3Pym7v7ULe2DspR1uyCSik3tmhL/lT739zdhwZ1Asu4WtUkmsCVcmOf/rmXnYfS845n/LW/jKtVTaNNKEq5mYzsXIa9/RtB/r6siyu87mRHXXNSFaAJXCk3YS24sJzDx7MKlX9wawyLtiQxqFNjLtQRlqqA8iyp1gz4BGgMOIApxpg3RSQM+BqIBvYA1xpjjrouVKW8W9KxjLzkfWX3KIZ2iaRtRAjNwmozoEOEzdEpd1SeGngO8Kgx5m8RCQFWicgCYDSw0BgzSUTGAeOAJ10XqlLeKfVENkPeWlZoQYXbz29Jp6i69gWlPMJpE7gxJgFIcO6nichmIAoYDlzkvGw6sBhN4EpV2Bcr9xVbuizIX2cSVKdXoTZwEYkGugMrgAhncscYkyAijUp5zRhgDEDz5s0rFaxS3qhRiNUt8KcHzieybpB2E1TlVu5uhCJSB/gOGGuMOVbe1xljphhjYowxMQ0b6gMYpYpaH2/1NGlQJ0CTt6qQciVwEfHHSt6fG2O+dxYnikik83wkkFTa65VSxRlj2HgglU//3AvAnHUJNkekPE15eqEIMBXYbIx5vcCp2cAoYJJz+4NLIlTKgzkchl82HWTp9sMs35lMo5BAhnVrwuBOkTwzawNz1ucnbZ2USlWUmIKPvku6QOR8YBmwHqsbIcBTWO3gM4DmwD7gGmPMkbLeKyYmxsTGxlY2ZqXcijGGW6au5Lcdh9kx8XL8fH1Yuz+F4ZN/L9frP/lHL11BR5VJRFYZY2KKlpenF8pvQGlTnl1a2cCU8mSbDhxj8H+X5R1f/uYyDqScJD0rN7+sU2PW7k+hT+twJgzrQOcJv+Sd2/L8IO1xos6YjsRU6gxl5uQWSt4Auw6nk+uw/qp98/puDO8WVex1b9/Ynb92H+GBS9to8laVoglcqTMwad4W3luys1DZ12N6c050GLsOHyczx0GHyJLnLRnapQlDuzSpjjCVl9MErlQFGGNoOX5uobLFj11E47pBebXpsxqF2BGaqoE0gStVDgdSTvLdqjhmronPK2sbUYevx/ShfnCAjZGpmkwTuFKnsTc5nQtfXQxAgJ8PV3WP4uzIUPp3iNDkrWylCVypMhhjuPb95XnH79zYg/46M6ByE5rAlSqFMYZnf9xE4rFM6gT6seHZgXaHpFQhmsCVKsGBlJO8tWg7X67cT4CvD+MHt7c7JKWK0QSuVBF7Dqdz7fvLSUrLpEvTunx/z3n4+eryscr9aAJXymnh5kSe+HYdPj5CTq6DaaNj6NWygSZv5bY0gSuFNenU7dOteXqC/H349u7zdEUc5fY0gasaz+EwtHrKGpzTITKUuQ9dYHNESpWP/m2ovN53q+JYtCWx1PMJxzLy9j+5vVd1hKRUldAauPJqH/++mwk/bgKKTy7lcBgMsGZfSl5ZuK6IozyIJnDltdIzc3hhzmZiWtRn56HjjP16Db9sSmTcoPY0DAlk4BtL2Zt8Iu/6gR11gI7yLJrAldfZkZTGmwt3sHznYXIchrH925Kd6+Cxb9YyZ10CS7ceIi0zp9BrhnSO5PXrutoUsVJnRhO48mjGGB7+eg1bDqZx/lnhjIxpyqA38ufovqR9I85vEw7A/LH9+HVzIvM2HMRXoHagH7edF01MdJhd4StVKeVZUm0aMBRIMsZ0cpaFAV8D0cAe4FpjzNHTfZguqaaqUq7DsGhLEnd+Uvx76pzo+rx5fXcahQRqP27l8UpbUq0839kfA4OKlI0DFhpj2gALncdKVZsDKScZMfn3vOT91ZjevHaN1QTy6IC2fHjrOTSpV0uTt/Jq5VkTc6mIRBcpHg5c5NyfDiwGnqzKwJQqKiktg4SUDDJzHHkzBD4z5GzCggM4t2UYvVsJI3s2tTlKparPmbaBRxhjEgCMMQki0qi0C0VkDDAGoHnz5mf4caqm+25VHI9+s7ZQWS1/X+64oJVNESllP5c/xDTGTAGmgNUG7urPU54tIzs3b2myNftTGDH597xztfx9qR3gS3J6FgDzx+qISVWznWkCTxSRSGftOxJIqsqgVM1zJD2L95fs5P2luwBoFBJIUlpmoWvWTbgMf23TVirPmSbw2cAoYJJz+0OVRaRqnKwcB/d8tooVu48AVk374naNqBPkx139WhEc6Ievj2jyVqqI0yZwEfkS64FluIjEAf/GStwzROR2YB9wjSuDVJ5v4eZEflqXwMzV8bxzUw+iGwSzPSmNR2asJddhtay1bhjMzPv6Ehrkb3O0SnmG8vRCuaGUU5dWcSzKC8WnnGT4279x+HhWXtm9n/9d4rUTr+ysyVupCtCRmMplcnId9J20KO/4vZt7ct5ZDXh/yU4iQoPo2CSUq99dTo/m9XhmaAd6NK9vY7RKeR5N4KrKZOU4CPDz4fDxTJbvTOaHNfGA1TQy58EL8nqXPD4wf33JPZOG2BKrUt5AE7iqtNST2Tw7eyPfr45nWNcmzF57AIBAPx/GX96euy5sbXOESnknTeCqRPuPnGD6H3t4sH8bQgL9MAY+X7mP1XuPct8lZ7FwcyK+Pj78vPEgK529RwBmrz3AuS3DuKxjY27t00J7jijlQprAVTHLth/ilqkrAZgRux9/Xx+OnMji1Lxn36+OL3R964bBRNWvTbuIOjzUvy11AvXbSqnqoD9piozsXFbuPoKPCN+u2s/PGxNpHBpE12Z1ST2ZTVS92jQKDWR7YhpLtx0mK9eR99q7L2zNIwPaEuCnNW2lqpsm8Bpu/Pfr+XLlvkJlQzpH8shlbWndsI5NUSmlykMTeA20OeEYCzYlknoyu1DyHtolkmeGdKBx3SAbo1NKlZcmcC+3NzmdOoF+xO49yo9rD/C/LUmkZ+UC1pD1sOAAnhvekaFdmtgcqVKqojSBeyFjDCLCVyv3Me779YXO+Qi0iwjhmpimOhWrUh5OE7gXyczJpefzv3I8M4eYFvWJ3XuUXi3D6N6sHiFBftx6XrQOVVfKi2gC9yKf/LGX487V1o9n5vDIgLbc1jeaEE3aSnklTeBeIDMnl/1HTjBx7mYAwoIDmD+2n81RKaVcTRO4GzqansUvmw5ycbtGZOU6CK8TyPKdyWxKOEagnw83925BoJ8PIsLa/Snc+/nfxKecBKBbs3rMuq+vzXeglKoOmsCryf4jJ0hIzaBL07p5kzoVZIzh37M3svtwOsu2Hy7zvV6Ys7nQsY/AKyO7cOxkNhe3L3V5UqWUl9EE7mLGGGbE7ufJ7/J7g9xzUWsaBAcw6rxoth5MY/W+o0xZtov9R04S4Jw7xM9HeGrw2fj6CIePZ9KxSV0ahgSyLTGN8c6eJR0iQ2kfGcLD/dvSLKy2LfenlLKPGHPm6wyLyCDgTcAX+NAYM6ms62NiYkxsbOwZf547y8l1sO/ICWaujuetRTu48dzmzN9wkCPp+QsZ1K/tT9P6tVkfn1rs9S3Dg7mrXyuujWlGelbOaR88Flz8Vynl3URklTEmpmj5GdfARcQXmAwMAOKAv0RktjFm05mH6bmemrmeGbFxecdfrMgf4Xhh24Y8ellbujStB1gPHb9csY856xNo3bAO/zi/JS3Dg/Nm7itPrxFN3kqpyjSh9AJ2GGN2AYjIV8BwwKsTeE6ugy0H01gfn0pEaCC/bEzkq7/2A9CqYTAhQf5MuKIDmTkOWjSozZH0LNo3DsXXR/LeI9DPl9F9WzK6b0u7bkMp5QUqk8CjgP0FjuOAc4teJCJjgDEAzZs3r8TH2ccYw88bE5n2+27Wx6VyMju32DVR9Wox76ELCPQrXDOOrFurusJUStUwlUngUkJZsQZ1Y8wUYApYbeCV+LxqdyDlJG8t2k7snqNsTzpOeJ1ArjunGS3Dgzl4LIOL2jakQ5NQcnINdWv54+NT0j+JUkq5RmUSeBzQrMBxU+BA5cKpfg6HYcXuI9QP9udEVi4Bvj5MnLOZni3qszs5nTnrEujVMoyJV3bi6h5Nte1ZKeU2KpPA/wLaiEhLIB64HrixSqJyocwca/GCL1bs46J2DdmckMbHf+wpdt3yXcl5+/+9vrtOsaqUcjtnnMCNMTkicj/wM1Y3wmnGmI1VFlkB2xPTSEjNoF/bhiXFQXzKSd78dTuB/j4E+fmy5WAa4XUCSMvIoXWjOrQMD2Zv8gn+3nuUlXvy12+ct+EgAL1bhXFtTDNCg/zJcTg4fDyLAF8fsh0OWjYI1uStlHJLlRrIY4yZC8ytolhK9c7incxcHc/wbk14ZkgHQoL8+HLlPt5ZvJOMrFzSnBM4gTUqsXNUXbYlplG3lj/LdhwmK8eBv6+QnWs1wd91YSsGdmxMnUA/jqZn0atlGCLafq2U8iweMRLzpas60yysNu8u3sEPa/Kb2XtFh9E+MoSzGtUhpkUYIUF+5DoM0eHBeddk5zqIP3qSqPq1dIV0pZRX8YgEHuTvyyMD2jKsaySf/bmPYyezOb9NOFd2jzptzdnf16dQQldKKW/hEQn8lLMahTBhWEe7w1BKKbegbQpKKeWhNIErpZSH0gSulFIeShO4Ukp5KE3gSinloTSBK6WUh9IErpRSHkoTuFJKeahKrYlZ4Q8TOQTsreK3DQfKXsbds3n7/YH336Pen+ez+x5bGGOKzeZXrQncFUQktqTFPr2Ft98feP896v15Pne9R21CUUopD6UJXCmlPJQ3JPApdgfgYt5+f+D996j35/nc8h49vg1cKaVqKm+ogSulVI2kCVwppTyUxyZwERkkIltFZIeIjLM7nqogIs1E5H8isllENorIQ87yMBFZICLbndv6dsdaGSLiKyKrReQn57HX3J+I1BORb0Vki/P/sY+X3d/Dzu/NDSLypYgEefr9icg0EUkSkQ0Fykq9JxEZ78w7W0VkoD1RWzwygYuILzAZuBzoANwgIh3sjapK5ACPGmPOBnoD9znvaxyw0BjTBljoPPZkDwGbCxx70/29Ccw3xrQHumLdp1fcn4hEAQ8CMcaYToAvcD2ef38fA4OKlJV4T86fx+uBjs7XvOPMR7bwyAQO9AJ2GGN2GWOygK+A4TbHVGnGmARjzN/O/TSsH/4orHub7rxsOjDClgCrgIg0BYYAHxYo9or7E5FQoB8wFcAYk2WMScFL7s/JD6glIn5AbeAAHn5/xpilwJEixaXd03DgK2NMpjFmN7ADKx/ZwlMTeBSwv8BxnLPMa4hINNAdWAFEGGMSwEryQCMbQ6usN4AnAEeBMm+5v1bAIeAjZxPRhyISjJfcnzEmHngN2AckAKnGmF/wkvsrorR7cqvc46kJvKSl6L2mP6SI1AG+A8YaY47ZHU9VEZGhQJIxZpXdsbiIH9ADeNcY0x1Ix/OaE0rlbAceDrQEmgDBInKzvVFVO7fKPZ6awOOAZgWOm2L9KefxRMQfK3l/boz53lmcKCKRzvORQJJd8VVSX2CYiOzBava6REQ+w3vuLw6IM8ascB5/i5XQveX++gO7jTGHjDHZwPfAeXjP/RVU2j25Ve7x1AT+F9BGRFqKSADWQ4XZNsdUaSIiWO2nm40xrxc4NRsY5dwfBfxQ3bFVBWPMeGNMU2NMNNb/2SJjzM14z/0dBPaLSDtn0aXAJrzk/rCaTnqLSG3n9+qlWM9pvOX+CirtnmYD14tIoIi0BNoAK22Iz2KM8cgvYDCwDdgJPG13PFV0T+dj/Tm2Dljj/BoMNMB6Er7duQ2zO9YquNeLgJ+c+15zf0A3INb5fzgLqO9l9/cssAXYAHwKBHr6/QFfYrXpZ2PVsG8v656Ap515ZytwuZ2x61B6pZTyUJ7ahKKUUjWeJnCllPJQmsCVUspDaQJXSikPpQlcKaU8lCZwpZTyUJrAlVLKQ/0/neXGqxHYTNkAAAAASUVORK5CYII=\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"import sys\n",
"sys.path.append('/nfs/scistore12/hpcgrp/jyeung/projects/StochasticForceInference')\n",
"# Import the package:\n",
"from StochasticForceInference import *\n",
"\n",
"# run SFI\n",
"# Diffusion parameters: a linear diffusion gradient (multiplicative noise)\n",
"\n",
"dim=2\n",
"diffusion_coeff = 0.5\n",
"D = diffusion_coeff * np.identity(dim) \n",
"\n",
"# Force field parameters (stochastic Lorenz process)\n",
"a, b, g = 1., 2., 3.\n",
"a = 10\n",
"b = 0.1\n",
"g = 0.2\n",
"\n",
"force = lambda X : np.array([[ a - b * x[0],\n",
" b * x[0] - g * x[1]] for x in X ])\n",
"\n",
"# Simulation parameters\n",
"initial_position = np.array([[-2 for i in range(dim)]]) \n",
"initial_position2 = np.array([[0, 75]])\n",
"\n",
"dt = 0.01\n",
"oversampling = 4\n",
"prerun = 0\n",
"Npts = 10000\n",
"tau = dt * Npts\n",
"tlist = np.linspace(0.,tau,Npts)\n",
"\n",
"# heterogeneous case: time list is staggered so each time point only measures one of the two particles\n",
"tlist1 = tlist[::2]\n",
"tlist2 = tlist\n",
"\n",
"# tlist2 = tlist[1::2]\n",
"\n",
"\n",
"# Run the simulation using our OverdampedLangevinProcess class\n",
"np.random.seed(1)\n",
"X1 = OverdampedLangevinProcess(force,D,tlist1,initial_position=initial_position,oversampling=oversampling,prerun=prerun )\n",
"X2 = OverdampedLangevinProcess(force,D,tlist2,initial_position=initial_position2,oversampling=oversampling,prerun=prerun )\n",
"\n",
"plt.plot([x[0][0] for x in X1.data], [x[0][1] for x in X1.data])\n",
"plt.plot([x[0][0] for x in X2.data], [x[0][1] for x in X2.data])\n"
]
},
{
"cell_type": "code",
"execution_count": 2,
"id": "597cca4b-9b5b-4bcc-9291-89318df3cc47",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[0. 0.020002 0.040004 0.060006 0.080008 0.10001 0.120012 0.140014\n",
" 0.160016 0.180018]\n",
"[0. 0.010001 0.020002 0.030003 0.040004 0.050005 0.060006 0.070007\n",
" 0.080008 0.090009]\n"
]
}
],
"source": [
"print(tlist1[:10])\n",
"print(tlist2[:10])"
]
},
{
"cell_type": "code",
"execution_count": 3,
"id": "986e3cf9-b0cc-470e-98e8-b1961af326ff",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"(5000, 1, 2)\n",
"(10000, 1, 2)\n"
]
}
],
"source": [
"print(X1.data.shape)\n",
"print(X2.data.shape)"
]
},
{
"cell_type": "code",
"execution_count": 4,
"id": "4b9ee065-13fe-4d90-b9e4-dc67b7e3625d",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"[<matplotlib.lines.Line2D at 0x7f1c28afe0b8>]"
]
},
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAD4CAYAAAD1jb0+AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAuzUlEQVR4nO3dd3gU1frA8e+bDiEBQiCEUAJIkV4igig2EAQEVOwFvCp2xQ7qvRcLiuX606tYUFDsYgFRiiJciopgkN47JIQEAgkhkLrn98cs6QkJyWZ2N+/nefLMzJnZ3XcgeXNy5hQxxqCUUsrz+NgdgFJKqTOjCVwppTyUJnCllPJQmsCVUspDaQJXSikP5VedHxYeHm6io6Or8yOVUsrjrVq16rAxpmHR8mpN4NHR0cTGxlbnRyqllMcTkb0llWsTilJKeShN4Eop5aE0gSullIfSBK6UUh5KE7hSSnkoTeBKKeWhNIErpZSH8uwEvuRV2DoPcrPtjkQppapdtQ7kqVJZ6fDXB3A8EUIi4b4VEFTX7qiUUqraeG4NPCAYHt4I/Z+FtAR4qyfsXmZ3VEopVW08N4ED+PpD34fg2k8hMw0WT7I7IqWUqjaencABRKDDMOh1J+z/02paUUqpGsDzE/gp4W3BkQMr3rM7EqWUqhbek8A7XwutL4Gl/4Eju+yORimlXM57Erh/EFz8NBgHTO4NPz0CGal2R6WUUi7jPQkcoGkMPLAK2l4GsVPhq5vsjkgppVzGuxI4QN0ouO4za3/PMojTBSSUUt7J+xL4KXf/bm0/HgLG2BuLUkq5wGkTuIi0E5E1Bb6OichYEQkTkQUist25rV8dAZdb407Q5XrIyYCZd0FOlt0RKaVUlTptAjfGbDXGdDPGdAN6AieAmcA4YKExpg2w0HnsXga/Auc/DOu+hh8fgsM77I5IKaWqTEWbUC4Fdhpj9gLDgenO8unAiCqMq2oE1YVz77b2134B718Ah7baG5NSSlWRiibw64EvnfsRxpgEAOe2UUkvEJExIhIrIrGHDh0680jPVEhjuOk76HErZJ+AqZfB+m+rPw6llKpiYsr5gE9EAoADQEdjTKKIpBhj6hU4f9QYU2Y7eExMjImNtbFXyNE9MGMUJKyBAc9Z86gopZSbE5FVxpiYouUVqYFfDvxtjEl0HieKSKTzzSOBpMqH6WL1o+HWWdb+gn/ZGYlSSlVaRRL4DeQ3nwDMBkY590cBP1RVUC5Vqz60GWjtnzhibyxKKVUJ5UrgIlIbGAB8X6B4EjBARLY7z3nOXK4dr7S22xdoH3GllMcq14o8xpgTQIMiZclYvVI8T5sBUCcCZo6BXYvhynftjkgppSrMe0diliU4HO5bCefcaXUvXPqq3REppVSF1cwEDlCrHgx41tpf/LL2D1dKeZyam8DBWlfz9gUgPvDBpbBiCmSdsDsqpZQql5qdwAGa9YJ7l1vzic97HP7+xO6IlFKqXDSBAzRoDZe/Yu3Pf1K7FyqlPIIm8FM6XgltL7f2t83X7oVKKbenCfwUkfzuhLPu0cWRlVJuTxN4QbXqw13LrP357jc7rlJKFaQJvKjILtYgHwCHw95YlFKqDJrAS3L2MGv71weQm21vLEopVQpN4CVpO8jaznsCJp8LR/faG49SSpVAE3hJ2vSHu5Zac4anH4aPh0JOpt1RKaVUIZrASxPZ1VrwofNISN0HC/5td0RKKVWIJvDTaX2xtV3xLmyZa28sSilVgCbw0zn7CqtrYURn+OoGeLMrfPsPyEq3OzKlVA2nCbw8IrvAnQshqqe1ruaG7+C1djpaUyllq/KuyFNPRL4VkS0isllE+ohImIgsEJHtzm2ZCxp7PL9Aa+bCcfshoA5kpcGf79gdlVKqBitvDfxNYL4xpj3QFdgMjAMWGmPaAAudx97NxxeCQuGx7dCsN/z8lC4GoZSyzWkTuIiEAv2AqQDGmCxjTAowHJjuvGw6MMI1IbqhgNoweg50GAGLXoCvboLsDLujUkrVMOWpgbcCDgEfichqEflQRIKBCGNMAoBz26ikF4vIGBGJFZHYQ4cOVVngtvP1g6s+sPa3/ARLX7E3HqVUjVOeBO4H9ADeNcZ0B9KpQHOJMWaKMSbGGBPTsGHDMwzTTfkFwAN/W/vL/gPHEuyNRylVo5QngccBccaYFc7jb7ESeqKIRAI4t0muCdHNNWgNFz9j7c++X7sXKqWqzWkTuDHmILBfRNo5iy4FNgGzgVHOslHADy6J0BNc+Dj0HQs7foVFEyH7pN0RKaVqgPL2QnkA+FxE1gHdgBeBScAAEdkODHAe11y974XQpvDnZHj3PMg8bndESikvV64EboxZ42zH7mKMGWGMOWqMSTbGXGqMaePc1uyFJEMi4N4/oNcYOLIL5jxqd0RKKS+nIzGrUlBdGPwq9LgV1n1lNakopZSLaAJ3hf7PgvjA4km6qo9SymU0gbtC7TC45J8Q9xcs+Kfd0SilvJSf3QF4rfMfhpR9sPxtCI2CPvfaHZFSystoDdxVROCy5639n8fD/pX2xqOU8jqawF0pMASueBNqh8OMW7V/uFKqSmkTiqv1HA21G8DXN8NHl8MtM6GWd8+8q5TX+OE+WP2ZtX/RePANsJZbbH2J9Vd2SY4dgNkPQlQP6HIdhDYB/1ouCU9MNS5KEBMTY2JjY6vt89yGIxeeC8s/npBqXyxKqfKbULf0c2OWQJNuxctfawvHEwuX9Z9gPRc7QyKyyhgTU7Rcm1Cqg48vPLgm//iH++BEzR73pJTHW/KytS3aVbho8gb4dYJLVvDSBF5dwlpaMxe2ON/6k+zV1nBoq91RKaVKkrK/eO3bLwiu+6xw2YS68Fx9azuhLr//7yfS+k0o+T2Td1Z5mJrAq1OD1nDbHBg+GYwDPhmh62oq5Y7e6FS87MHVkHUi/zjm9mKX9F1yE3GLpxUuvOwFGB8H4WdVcZCawO3R5Tprm3YAFj5nbyxKqfLZvxJWf5p//PsbJV52ts++wgV97rd6pLmAJnA7+PrDiHet/d9eh3UztCauVHVy5ELCWlj8srUc4kdD4I+385pC6Hxt8dd8Mwr2LMs/LrhflhzXLbeovVDslLIfpg2CY3HQuDPc/ZvdESlVM5TVu6QKjM26lz2mMbP+dWuVdBvWXijuqF4zuHORtX9kj62hKFVjVEOz5SxHX9aYs1w+5kMTuN1CIiCkCWSlwd4/7I5GKe+37D/V8CHC1hcGufxTNIG7g1O18I8u166FSrlauyHV8jEz/453+WeUK4GLyB4RWS8ia0Qk1lkWJiILRGS7c6vjw89UaCSM/MjaX/Bve2NRypvNexK2znHZ2/fKmEx0xhcAjPt+vcs+55SK1MAvNsZ0K9CQPg5YaIxpAyx0Hqsz1ekqaDMQts2DjTPtjkYp77TivUKHuaHNSrws3jSo8FtHZ3xBEvn12EA/1zdwVOYThgPTnfvTgRGVjqamO+cOa/vNaFj1sdW9SSnlMr7H9hcrSzahDM58ib4Zb7Igtwd9M9487fucumbzc4NoHlabpvVrseV592kDN8AvIrJKRMY4yyKMMQkAzm2jkl4oImNEJFZEYg8dOlT5iL1Z28tg1I/W/o8PWX3ElVJV57HtZZ6eaobRM/M9UqlDPA25M/sx4mnI89k3513TLuNjbs4an3f8XPYtxNOQsf3bUCvAl6VPXMxvT16ClDZbYRUq73SyfY0xB0SkEbBARLaU9wOMMVOAKWD1Az+DGGuWlv3g6qnw3e1wuOxvNqVUBQU3LPN0tGM/fuSQUyQ19vNZl7e/NWh0oXOzcvsC8Mav2/H39eHamGbUCfTDxwcC/XyrJu5SlKsGbow54NwmATOBXkCiiEQCOLdJrgqyxonqYc0hvvH7/LmIlVKVZ8peZPxS39U0keS84xBO8Ib/21zou67U1xwhNG//1Z+3cs7EXzn7X/N5/JvSX1NVTpvARSRYREJO7QOXARuA2cAo52WjgB9cFWSNE9YKxm6ARh2sqWePHbA7IqW8Q8F5+UvxoF9+J4L1QXcwwrfs8Rk3+y4osXz2Wtf/3JanBh4B/CYia4GVwBxjzHxgEjBARLYDA5zHqqoE1Ia2A6396VdAenLZ1yulqsRI36VYj/3Kp7fPphLLp/+jVxVFVLrTtoEbY3YBXUsoTwYudUVQyumip6w/+f58D15vD1E9rZU9mve2OzKlPNO9K+Cdc097mT+59PNZW663bCjWCltv39idwZ0i2ZqYRqOQQBrUCaxUqOWhIzHdmV8ADHgO7lhgzT28bzlMGwhHdtsdmVKeqVF7uPL9Mi9ZktuFALKZGlC+IffnNq3NnklDGNqlCT4+wtmRodWSvEEXNfYMkV2hcRdw5MBfH8DnI+GBVXZHpZTnObgBZt5V6umeGe+STF26yY7yvV94W7hjYRUFV3GawD2FCAx5zUrgoU3sjkYpz5RcdmJeFXRPxd7vrmXgY19DhjaheJrQprB7KXzYHxJc301JKY/icMCLTeHZ+pDiXBnn/X6w5FX4T3trUYYCFvX9nI/qPUBqx5tLeLMCnijSbHnjNzAhFfyDqjD4itMauKe57lNY/jZs+A6mDoALn4Seo6H26btHKeX1EjdYUzMDvNE5vzyh5AeSl7RvBANeAOCu1eG8H/BG8YvuWGj9fE1ItX5BHN1trW/rBrQG7mmiesDIaVYNwL82LHzWqmEc2mZ3ZErZKzcHIjpW7DVTB8DxQ2yIT+VnRy+iMz4vfk3TAgvh+Pi4TfIGTeCeq+1l8MQuay7xrONW75QTR+yOSinX2PEr/F8naym0GYWbQfjfS1b58w1g0QvwwN8Ve+/gcIa+dWo5Q+G6zH/mnxv8WqXCdjVdE9MbxK+CDy6x9oe+ATG32RqOUlVq+6/w+dWFy0bPgejzrf1KrG/puGcF135/hNi9RwuVfzWmN71bVXxKWVcpbU1MbQP3BlE94ZJnrNrHT2OtofitLrQ7KqUq79DW4skb4OPKr6qzatAPXP1/O0s8507JuyzahOIt+j0OD2+COhHwyXCrVvLZ1ZCRandkSp0ZY2BygeHo4+PK/dJpOYPINGXXT6+elX6mkbkNTeDepG4U3LM8f6j9jl9hUnNdZ1N5DmPg709g6mXwbL388pEfQWAITEgl66J/lfkWv+d25LmcW2mX+QmH716P6TG60Pl3coblLXsGcGX3KLY8P4izI0O5q18r3r+lZxXekGtpE4q3CW4A/5hv7Z9qGzx2ABq2sy8mpcorZR/MfqBwWbebrSUHndrOb4/wGbuD8vtuF0zIBYU3bg7D3rS+nOqt2AczrfUqt74wKG/O7nkPXVBVd1FtNIF7s4Evws9PwewH4ba5UK/k9f+UchtFxzM8kwR+1rwiDofhvi+sHiYGH1pnfMptvvP5MvcS/jW0A3/uSmb1/hT+ero/SWkZNAgueT6SG89tTu9WYWw8cMzlCy64mvZC8WaO3MLzH1/+itXVsNsNUD/atrCUKtOJI/BKy2LFCfftpM9/VhQrX/3PAdQPDqiOyGxTWi8UTeA1wR9vwS/P5B/XbwniA0ecT+AveAw6XQ0RHeyJT6mifn7aGnFc0BX/5e3U83jtl220Cg9mwSMX4uvj+nUn3YF2I6zJ+txvTYCVk2k9JJo/HjIL9E5Z9pr1BfDYDqhT9rqBSrmcr3/xsh8fJGr4RgD6d4ioMcm7LOWugYuILxALxBtjhopIGPA1EA3sAa41xhwt/R20Bu42MtPgeJLVX/zobvh+DMT9lX8++gK48WsICLYvRqVOKfAXpDXUXZh1X1+6Natna1jVqbQaeEW6ET4EbC5wPA5YaIxpAyx0HitPEBhizecgYiXxO36FW2bln9+zDF5sAn9/Ckf32hamUnxxfaHmvz1BNwEwb32CXRG5lXI1oYhIU2AIMBF4xFk8HLjIuT8dWAw8WbXhqWrT+mJrDomTR2HTLKvWM/t+61xoUzgWB+IL3W6EvX9A485w7t3Qoo+tYSvvltFqAEHb5hUrv7pnUxuicT/lakIRkW+Bl4AQ4DFnE0qKMaZegWuOGmPql/DaMcAYgObNm/fcu1drdB7h5FHYsRCOJ1rbnaWsOnL/Kgg/q3pjU14vIfUkWxLS+GLlPhZsSuSShmlMS7NW0km+ZyMNImpWAj/jh5giMhRIMsasEpGLKvrBxpgpwBSw2sAr+nplk1r1ofNIa7/Pfdb26F4IiXQm9wXww30wfag1hN/GVUmU9/nnrA38ujmJALIBf6ZdUR+cY3UahEfYGps7KU8TSl9gmIgMBoKAUBH5DEgUkUhjTIKIRAJJrgxUuYH6LaxtSAR0vxkWTYS0A7D6E2tRCaWqSMusbewJuje/oOBAS+Oo9njc1WmrTcaY8caYpsaYaOB6YJEx5mZgNnBqYt5RwA8ui1K5p2s+sra7FtsahvJg236G/3aH5AKzAm5fwNPx95Z8/b0r8kZmqsr1A58EzBCR24F9wDVVE5LyGH++Y203zoSrPii5765SBWWdgLQE8PGFNV/Akpet8rd6lP26O/9nrUalCqlQAjfGLMbqbYIxJhm4tOpDUh5j8H+smlPiBusH8J7lEFjH7qiUuzIGXows9+VtMj4hGz9eHdmFa6J0Hp+S6JMndebqNIS7lln7Kfvgl6fhwBrIOGZrWMpNJW44/TU3fsP2od8RnfEF2c765YjuUS4OzHNpAleV4+MDdzvXE1z1MUy5ECY1gxVTbA1LuaGITsXLzhpgrSYFcNFTZLXqz1U/5uad/nP8pfj7apoqjc6FoiqvcWcY9pY1j7N/MGSnw7zHoUl3aHaO3dEpu/35nvW8pFMJS6Pd/K217fc4ALNXxZGWmQPA8vGX0LhuUHVF6ZH0V5uqGj1uhX+nwNMH4LrPrLKp/a0FaVXNlZsN85+ElL3w2+tW2aBJMG4fTCi83J/DYdielJZ3/POGg9UZqUfSGriqOuKcHe7sKyCyGySsyV+Q9o6F1hzkweE2Bads8dPD+fsjP7KmMe44othlO5KO0//1JYXKbu7dwsXBeT5N4Mo1xiyGRS/kT1P7YZEOS0/sLr76inIfxuT/Qi7JsQPwwaXWL+vBr5R8zdE9sPpTa3/gi4WWRcv/GMOG+GM8+s2avLJfHu5HcKAfftr2fVqawJVriMCl/7RmPZx1j1XW9UZY6xxS90pL6DUGDm8HR441Z3mT7tYQ/iO7IDfLalsvK4ko11g3w3qekZMB130Omcfy/w9vmWVNRzzjFut45fvWaNxTzWanJO8s3Lc7quSFgvtOWsSB1Iy845BAP9pGhFThzXg3XZFHVb8VU6wBHCcOAwKU8j0ovtDucqvmVtIDMFX1krbAO+dW/HWdroaR0/KPs04U7vP9wN/WL/MCjDG0HD8373jSVZ3p3rw+7RprAi9KV+RR7uPcMdDrTqvvuMmFlR/Cn5Pzz4e3tWrhjhzY8pP1tfknaHkBdL0B/GvZF7u32/tb8TL/2pB9onBZi/Phtjnwegc4Fg8bvoOrp0L6Ydj1P2v1p4KKJG+A487eJgDf3t2HmGhtUqsorYEr95F9EvyCrGYTR641lW3SJvisSO379gXQrFfJ7+FwQFYaBNV1fbzeJDsDfn8TFr9oHYe1gvtjrSHvp5w4YrVpt7wQmnSzyiaU49/5tvklzht/x/RYft2cCMDulwYj2lxWKq2BK/dXsGbt42ut4xnaxEokaQkw/Qrr3NQB+V3Q4v+Gw9ug/RBrpaEXI62220uegQ4jILxNtd+GW3I4rF+MRZOkMfBsvcJlV0/Nn0q4oNph0PehwmXXfprfHl5UCc0mp6SezM5L3gC5DoOfrybwitIauPIsp2p8rS+BnYtOf337odYDtppYuyv6IBHg7t+hwVngH2T10X6+SLfOh9Za3T0rInGj9X8R1hrqRsGR3dDmMgioXepLRk1byZJth/KO/3q6Pw1DdJbB0lTFmphK2e+ORVCnMSTvKFwe1ip/v+WFMOJda3/LT/BydLWF5zYcuSXP8PdeX5gYYdXIf3uj+Pmt8yv+WREd4bwHoP1giOxq9fMukrzXx6VyLCOb45k5vDx/Cxe2bVjo/LM/bqz45yptQlEepmlPeGyrte/IhdS4/IUmDqyGhLX5i0uEtYJpAyEjBd49H9oNgvMftnq3+DuHaP86AX77P2t/+GRroQpPlXYQPhtpDZYqOGCqaS+rSWTeE/llzxVZ/TCyK/S6C7pcV+VhrYtLYdjbv9P/7Eb8urnkdV/+2Jlc5Z9bE2gTivJu5XnIVtBTCWX+6V/tDm2FuY/D7iXWL6QHV5d8XU4mvNCoePn4OOvZAEDsR/DT2OLXNDsXbv+lykI+kp5FkL8Ph9OyqBfsz40f/MmG+NJnqHzpqs7c0Kt5lX2+N9KHmKpmOvWwc8G/rF4Wp/NuH6sd2G6JG61h5+/0zi87sqv4dVnpEBAMk0volXPvivzkDRBzm7VN3gEdr4IGrayBU1Vkzf4UHp2xhp2H0oud69K0LuviUouVfzT6HC5uX8IvHlUumsBVzTDgOeh8jVWLPbzdWpbrZIq1mtDAifD1LbBtnjX8+7/doV5za5j4OXdUb5y52VZT0NQBZV+35BX438SSz10zvcT5RoD8JF7Fdh46zojJv5d6/rpzmhVK4Esfv5hmYbW062AlnbYJRUSCgKVAIFbC/9YY828RCQO+BqKBPcC1xpijZb2XNqEot5VxzJrHvCQPb7IW0v3xQauHxuWvWMvHZRyzkuiWuXDfiqppeplyMRz4u/Tz//jZatcvzei5EN238nGU0/HMHP5vwTam/rY7r2zKLT35Y2cyj1zWll82JlLL35fLOzUmPSuH3YfT6dK0XrXF5y1Ka0IpTwIXINgYc1xE/IHfgIeAq4AjxphJIjIOqG+MebKs99IErtza8SR4rYR+44Nfs5pfUvfnl3W7CdZ8nn98x0JoGgPHEqz+0gUX3jXGmk/k+CEIP6vkz969NL+fe3k17QU3fm21k6fuh04jrQU2qtHV7/7Bqr1Wve3VkV0Y2qUJtQJ8T/MqVVFnnMCLvEltrAR+D/AJcJExJkFEIoHFxph2Zb1eE7jyKK+1tUaDnon7V1nNMC80LH7utvkQ2QVebAIN2kDy9uLXXDbRmjqgUUdrZOmHA/KvqxVmzfZ4qveNjaLHzcnb19GUrlOpBC4ivsAq4CxgsjHmSRFJMcbUK3DNUWNMsSciIjIGGAPQvHnznnv37j3zu1CqOm2ZC1/dULjsmUOw6Qf43tk23v0Wa16X9/tV/vMueBQufsZqrvF138dTuQ7DiMm/sz6+8EPJuQ9eQIcmoTZF5d2qqgZeD5gJPAD8Vp4EXpDWwJXHSU+2prYNCrV6e5ySGm+N7gxtYh07cq35z0+tOlNQ9AUw+idrP2EdvH9B/jn/YOg5Cuq1gHPvcusRo1k5Dg6mZvCv2RtYvPVQoXMXtAln2uhzdP1KF6mSboTGmBQRWQwMAhJFJLJAE0rJPfSV8mTBDUour1tkpXQfX+j/byvZb5oNY9dZ5cYUbpeO7GLV4tfPgA7DC3fzc1PJxzPp+ULZS+N9evsZTEGrKu20CVxEGgLZzuRdC+gPvAzMBkYBk5zbH1wZqFIeYeBE6+uUkmrUfgG2jfiM3XOEke8tB6B94xBeGdmFHIehe7N6xdqv31+yk7VxKZzbspRfYsp25amBRwLTne3gPsAMY8xPIrIcmCEitwP7gGtcGKdSqhIOH89kR9Jxrp/yZ17ZloNpDHs7v+/2Ozf1oHNUXXIchnVxKbw0bwsAc9fnLy58Q6/mZOc6OJ6RQ5dmdXll/lb+c03X6rsRVYgOpVfKC62PS2XuhgSGdW3C+rhUnvhuXd6562KaMbhLJCknsgj08+Huz0rud960fi3ijp4sVLZn0pBCxzm5Dl27shroUHqlapAr3rZW1nl38c5i514e2aXQ8Xf39GHke8spWJdrFR7MvLEX4O9sv1+9/yhR9YoPVNLkbS9N4Ep5mcyc3ELHCx7uR0TdIJKOZdKiQfEk3LNFGLtfGlKsvOg1yv1oAlfKy7y+YFve/sQrO9HGucp7aJC/XSEpF9EErpQXSTqWwZKth6gT6MdfT/fXYe1eThO4Ul7A4TC0emouAL4+wrPDOmryrgE0gSvl4bYnpvGP6X/lHX80+hz6tS1hDhbldTSBK+XB/tyVnNe3289HWPbkxUTWrWVzVKq6aAJXykMdz8zJS94DO0bw9o09dC6SGkYTuFIeJCM7l3kbEnj468LLvk28srMm7xpIE7hSHuSfszbwzaq4vOOZ955H9+ZVt66l8iyawJVyc8YY1sen5s1b0io8mJn39iW0lp8uoFDDaQJXys29OHczHyzLX3Pym7v7ULe2DspR1uyCSik3tmhL/lT739zdhwZ1Asu4WtUkmsCVcmOf/rmXnYfS845n/LW/jKtVTaNNKEq5mYzsXIa9/RtB/r6siyu87mRHXXNSFaAJXCk3YS24sJzDx7MKlX9wawyLtiQxqFNjLtQRlqqA8iyp1gz4BGgMOIApxpg3RSQM+BqIBvYA1xpjjrouVKW8W9KxjLzkfWX3KIZ2iaRtRAjNwmozoEOEzdEpd1SeGngO8Kgx5m8RCQFWicgCYDSw0BgzSUTGAeOAJ10XqlLeKfVENkPeWlZoQYXbz29Jp6i69gWlPMJpE7gxJgFIcO6nichmIAoYDlzkvGw6sBhN4EpV2Bcr9xVbuizIX2cSVKdXoTZwEYkGugMrgAhncscYkyAijUp5zRhgDEDz5s0rFaxS3qhRiNUt8KcHzieybpB2E1TlVu5uhCJSB/gOGGuMOVbe1xljphhjYowxMQ0b6gMYpYpaH2/1NGlQJ0CTt6qQciVwEfHHSt6fG2O+dxYnikik83wkkFTa65VSxRlj2HgglU//3AvAnHUJNkekPE15eqEIMBXYbIx5vcCp2cAoYJJz+4NLIlTKgzkchl82HWTp9sMs35lMo5BAhnVrwuBOkTwzawNz1ucnbZ2USlWUmIKPvku6QOR8YBmwHqsbIcBTWO3gM4DmwD7gGmPMkbLeKyYmxsTGxlY2ZqXcijGGW6au5Lcdh9kx8XL8fH1Yuz+F4ZN/L9frP/lHL11BR5VJRFYZY2KKlpenF8pvQGlTnl1a2cCU8mSbDhxj8H+X5R1f/uYyDqScJD0rN7+sU2PW7k+hT+twJgzrQOcJv+Sd2/L8IO1xos6YjsRU6gxl5uQWSt4Auw6nk+uw/qp98/puDO8WVex1b9/Ynb92H+GBS9to8laVoglcqTMwad4W3luys1DZ12N6c050GLsOHyczx0GHyJLnLRnapQlDuzSpjjCVl9MErlQFGGNoOX5uobLFj11E47pBebXpsxqF2BGaqoE0gStVDgdSTvLdqjhmronPK2sbUYevx/ShfnCAjZGpmkwTuFKnsTc5nQtfXQxAgJ8PV3WP4uzIUPp3iNDkrWylCVypMhhjuPb95XnH79zYg/46M6ByE5rAlSqFMYZnf9xE4rFM6gT6seHZgXaHpFQhmsCVKsGBlJO8tWg7X67cT4CvD+MHt7c7JKWK0QSuVBF7Dqdz7fvLSUrLpEvTunx/z3n4+eryscr9aAJXymnh5kSe+HYdPj5CTq6DaaNj6NWygSZv5bY0gSuFNenU7dOteXqC/H349u7zdEUc5fY0gasaz+EwtHrKGpzTITKUuQ9dYHNESpWP/m2ovN53q+JYtCWx1PMJxzLy9j+5vVd1hKRUldAauPJqH/++mwk/bgKKTy7lcBgMsGZfSl5ZuK6IozyIJnDltdIzc3hhzmZiWtRn56HjjP16Db9sSmTcoPY0DAlk4BtL2Zt8Iu/6gR11gI7yLJrAldfZkZTGmwt3sHznYXIchrH925Kd6+Cxb9YyZ10CS7ceIi0zp9BrhnSO5PXrutoUsVJnRhO48mjGGB7+eg1bDqZx/lnhjIxpyqA38ufovqR9I85vEw7A/LH9+HVzIvM2HMRXoHagH7edF01MdJhd4StVKeVZUm0aMBRIMsZ0cpaFAV8D0cAe4FpjzNHTfZguqaaqUq7DsGhLEnd+Uvx76pzo+rx5fXcahQRqP27l8UpbUq0839kfA4OKlI0DFhpj2gALncdKVZsDKScZMfn3vOT91ZjevHaN1QTy6IC2fHjrOTSpV0uTt/Jq5VkTc6mIRBcpHg5c5NyfDiwGnqzKwJQqKiktg4SUDDJzHHkzBD4z5GzCggM4t2UYvVsJI3s2tTlKparPmbaBRxhjEgCMMQki0qi0C0VkDDAGoHnz5mf4caqm+25VHI9+s7ZQWS1/X+64oJVNESllP5c/xDTGTAGmgNUG7urPU54tIzs3b2myNftTGDH597xztfx9qR3gS3J6FgDzx+qISVWznWkCTxSRSGftOxJIqsqgVM1zJD2L95fs5P2luwBoFBJIUlpmoWvWTbgMf23TVirPmSbw2cAoYJJz+0OVRaRqnKwcB/d8tooVu48AVk374naNqBPkx139WhEc6Ievj2jyVqqI0yZwEfkS64FluIjEAf/GStwzROR2YB9wjSuDVJ5v4eZEflqXwMzV8bxzUw+iGwSzPSmNR2asJddhtay1bhjMzPv6Ehrkb3O0SnmG8vRCuaGUU5dWcSzKC8WnnGT4279x+HhWXtm9n/9d4rUTr+ysyVupCtCRmMplcnId9J20KO/4vZt7ct5ZDXh/yU4iQoPo2CSUq99dTo/m9XhmaAd6NK9vY7RKeR5N4KrKZOU4CPDz4fDxTJbvTOaHNfGA1TQy58EL8nqXPD4wf33JPZOG2BKrUt5AE7iqtNST2Tw7eyPfr45nWNcmzF57AIBAPx/GX96euy5sbXOESnknTeCqRPuPnGD6H3t4sH8bQgL9MAY+X7mP1XuPct8lZ7FwcyK+Pj78vPEgK529RwBmrz3AuS3DuKxjY27t00J7jijlQprAVTHLth/ilqkrAZgRux9/Xx+OnMji1Lxn36+OL3R964bBRNWvTbuIOjzUvy11AvXbSqnqoD9piozsXFbuPoKPCN+u2s/PGxNpHBpE12Z1ST2ZTVS92jQKDWR7YhpLtx0mK9eR99q7L2zNIwPaEuCnNW2lqpsm8Bpu/Pfr+XLlvkJlQzpH8shlbWndsI5NUSmlykMTeA20OeEYCzYlknoyu1DyHtolkmeGdKBx3SAbo1NKlZcmcC+3NzmdOoF+xO49yo9rD/C/LUmkZ+UC1pD1sOAAnhvekaFdmtgcqVKqojSBeyFjDCLCVyv3Me779YXO+Qi0iwjhmpimOhWrUh5OE7gXyczJpefzv3I8M4eYFvWJ3XuUXi3D6N6sHiFBftx6XrQOVVfKi2gC9yKf/LGX487V1o9n5vDIgLbc1jeaEE3aSnklTeBeIDMnl/1HTjBx7mYAwoIDmD+2n81RKaVcTRO4GzqansUvmw5ycbtGZOU6CK8TyPKdyWxKOEagnw83925BoJ8PIsLa/Snc+/nfxKecBKBbs3rMuq+vzXeglKoOmsCryf4jJ0hIzaBL07p5kzoVZIzh37M3svtwOsu2Hy7zvV6Ys7nQsY/AKyO7cOxkNhe3L3V5UqWUl9EE7mLGGGbE7ufJ7/J7g9xzUWsaBAcw6rxoth5MY/W+o0xZtov9R04S4Jw7xM9HeGrw2fj6CIePZ9KxSV0ahgSyLTGN8c6eJR0iQ2kfGcLD/dvSLKy2LfenlLKPGHPm6wyLyCDgTcAX+NAYM6ms62NiYkxsbOwZf547y8l1sO/ICWaujuetRTu48dzmzN9wkCPp+QsZ1K/tT9P6tVkfn1rs9S3Dg7mrXyuujWlGelbOaR88Flz8Vynl3URklTEmpmj5GdfARcQXmAwMAOKAv0RktjFm05mH6bmemrmeGbFxecdfrMgf4Xhh24Y8ellbujStB1gPHb9csY856xNo3bAO/zi/JS3Dg/Nm7itPrxFN3kqpyjSh9AJ2GGN2AYjIV8BwwKsTeE6ugy0H01gfn0pEaCC/bEzkq7/2A9CqYTAhQf5MuKIDmTkOWjSozZH0LNo3DsXXR/LeI9DPl9F9WzK6b0u7bkMp5QUqk8CjgP0FjuOAc4teJCJjgDEAzZs3r8TH2ccYw88bE5n2+27Wx6VyMju32DVR9Wox76ELCPQrXDOOrFurusJUStUwlUngUkJZsQZ1Y8wUYApYbeCV+LxqdyDlJG8t2k7snqNsTzpOeJ1ArjunGS3Dgzl4LIOL2jakQ5NQcnINdWv54+NT0j+JUkq5RmUSeBzQrMBxU+BA5cKpfg6HYcXuI9QP9udEVi4Bvj5MnLOZni3qszs5nTnrEujVMoyJV3bi6h5Nte1ZKeU2KpPA/wLaiEhLIB64HrixSqJyocwca/GCL1bs46J2DdmckMbHf+wpdt3yXcl5+/+9vrtOsaqUcjtnnMCNMTkicj/wM1Y3wmnGmI1VFlkB2xPTSEjNoF/bhiXFQXzKSd78dTuB/j4E+fmy5WAa4XUCSMvIoXWjOrQMD2Zv8gn+3nuUlXvy12+ct+EgAL1bhXFtTDNCg/zJcTg4fDyLAF8fsh0OWjYI1uStlHJLlRrIY4yZC8ytolhK9c7incxcHc/wbk14ZkgHQoL8+HLlPt5ZvJOMrFzSnBM4gTUqsXNUXbYlplG3lj/LdhwmK8eBv6+QnWs1wd91YSsGdmxMnUA/jqZn0atlGCLafq2U8iweMRLzpas60yysNu8u3sEPa/Kb2XtFh9E+MoSzGtUhpkUYIUF+5DoM0eHBeddk5zqIP3qSqPq1dIV0pZRX8YgEHuTvyyMD2jKsaySf/bmPYyezOb9NOFd2jzptzdnf16dQQldKKW/hEQn8lLMahTBhWEe7w1BKKbegbQpKKeWhNIErpZSH0gSulFIeShO4Ukp5KE3gSinloTSBK6WUh9IErpRSHkoTuFJKeahKrYlZ4Q8TOQTsreK3DQfKXsbds3n7/YH336Pen+ez+x5bGGOKzeZXrQncFUQktqTFPr2Ft98feP896v15Pne9R21CUUopD6UJXCmlPJQ3JPApdgfgYt5+f+D996j35/nc8h49vg1cKaVqKm+ogSulVI2kCVwppTyUxyZwERkkIltFZIeIjLM7nqogIs1E5H8isllENorIQ87yMBFZICLbndv6dsdaGSLiKyKrReQn57HX3J+I1BORb0Vki/P/sY+X3d/Dzu/NDSLypYgEefr9icg0EUkSkQ0Fykq9JxEZ78w7W0VkoD1RWzwygYuILzAZuBzoANwgIh3sjapK5ACPGmPOBnoD9znvaxyw0BjTBljoPPZkDwGbCxx70/29Ccw3xrQHumLdp1fcn4hEAQ8CMcaYToAvcD2ef38fA4OKlJV4T86fx+uBjs7XvOPMR7bwyAQO9AJ2GGN2GWOygK+A4TbHVGnGmARjzN/O/TSsH/4orHub7rxsOjDClgCrgIg0BYYAHxYo9or7E5FQoB8wFcAYk2WMScFL7s/JD6glIn5AbeAAHn5/xpilwJEixaXd03DgK2NMpjFmN7ADKx/ZwlMTeBSwv8BxnLPMa4hINNAdWAFEGGMSwEryQCMbQ6usN4AnAEeBMm+5v1bAIeAjZxPRhyISjJfcnzEmHngN2AckAKnGmF/wkvsrorR7cqvc46kJvKSl6L2mP6SI1AG+A8YaY47ZHU9VEZGhQJIxZpXdsbiIH9ADeNcY0x1Ix/OaE0rlbAceDrQEmgDBInKzvVFVO7fKPZ6awOOAZgWOm2L9KefxRMQfK3l/boz53lmcKCKRzvORQJJd8VVSX2CYiOzBava6REQ+w3vuLw6IM8ascB5/i5XQveX++gO7jTGHjDHZwPfAeXjP/RVU2j25Ve7x1AT+F9BGRFqKSADWQ4XZNsdUaSIiWO2nm40xrxc4NRsY5dwfBfxQ3bFVBWPMeGNMU2NMNNb/2SJjzM14z/0dBPaLSDtn0aXAJrzk/rCaTnqLSG3n9+qlWM9pvOX+CirtnmYD14tIoIi0BNoAK22Iz2KM8cgvYDCwDdgJPG13PFV0T+dj/Tm2Dljj/BoMNMB6Er7duQ2zO9YquNeLgJ+c+15zf0A3INb5fzgLqO9l9/cssAXYAHwKBHr6/QFfYrXpZ2PVsG8v656Ap515ZytwuZ2x61B6pZTyUJ7ahKKUUjWeJnCllPJQmsCVUspDaQJXSikPpQlcKaU8lCZwpZTyUJrAlVLKQ/0/neXGqxHYTNkAAAAASUVORK5CYII=\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"\n",
"plt.plot([x[0][0] for x in X1.data], [x[0][1] for x in X1.data])\n",
"plt.plot([x[0][0] for x in X2.data], [x[0][1] for x in X2.data])\n",
"# plt.plot([x[0][0] for x in X2_hetero.data], [x[0][1] for x in X2_hetero.data])\n"
]
},
{
"cell_type": "code",
"execution_count": 5,
"id": "923fc327-418b-44e2-888f-b821d548682b",
"metadata": {},
"outputs": [],
"source": [
"# try running it twice, then set up as dictionary? \n",
"homo1 = StochasticTrajectoryData(X1.data, tlist1, data_type = \"homogeneous\")\n",
"homo2 = StochasticTrajectoryData(X2.data, tlist2, data_type = \"homogeneous\")\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 7,
"id": "a3d6b6df-12b0-45bc-a00d-21d6afaa17f4",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"(4998, 1, 2)\n"
]
}
],
"source": [
"print(homo1.dX_plus.shape)"
]
},
{
"cell_type": "code",
"execution_count": 15,
"id": "51ab5ec0-63b1-478a-bfa3-b2e418194518",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"(5000,)\n",
"(10000,)\n",
"(15000,)\n",
"check size for X_ito\n",
"(4998, 1, 2)\n",
"(9998, 1, 2)\n",
"(14996, 1, 2)\n",
"check size for dX_plus\n",
"(4998, 1, 2)\n",
"(9998, 1, 2)\n",
"(14996, 1, 2)\n",
"check size for dX_minus\n",
"(4998, 1, 2)\n",
"(9998, 1, 2)\n",
"(14996, 1, 2)\n",
"check size for dt\n",
"(4998,)\n",
"(9998,)\n",
"(14996,)\n"
]
}
],
"source": [
"# try just concatenate them as list\n",
"# tcat = [tlist1, tlist2]\n",
"\n",
"tcat = np.concatenate((tlist1, tlist2), axis = 0)\n",
"\n",
"print(tlist1.shape)\n",
"print(tlist2.shape)\n",
"print(tcat.shape)\n",
"\n",
"data = {}\n",
"data['X_ito'] = np.concatenate((homo1.X_ito, homo2.X_ito), axis = 0)\n",
"data['dX_plus'] = np.concatenate((homo1.dX_plus, homo2.dX_plus), axis = 0)\n",
"data['dX_minus'] = np.concatenate((homo1.dX_minus, homo2.dX_minus), axis = 0)\n",
"data['dt'] = np.concatenate((homo1.dt, homo2.dt), axis = 0)\n",
"\n",
"# print(homo1.X_ito.shape)\n",
"# print(homo2.X_ito.shape)\n",
"# print(data['X_ito'].shape)\n",
"\n",
"attr_lst = ['X_ito', 'dX_plus', 'dX_minus', 'dt']\n",
"for attr in attr_lst:\n",
" print('check size for ' + attr)\n",
" print(getattr(homo1, attr).shape)\n",
" print(getattr(homo2, attr).shape)\n",
" print(data[attr].shape)\n",
"\n",
"# data['X_ito'] = [homo1.X_ito, homo2.X_ito]\n",
"# data['dX_plus'] = [homo1.dX_plus, homo2.dX_plus]\n",
"# data['dX_minus'] = [homo1.dX_minus, homo2.dX_minus]\n",
"# data['dt'] = [homo1.dt, homo2.dt]\n",
"\n",
"hetero = StochasticTrajectoryData(data, tcat, data_type = \"heterogeneous\")"
]
},
{
"cell_type": "code",
"execution_count": 18,
"id": "58c48e4d-4c94-4b3e-8d93-b6ec6ffc8aae",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"2\n"
]
}
],
"source": [
"print(hetero.d)"
]
},
{
"cell_type": "code",
"execution_count": 19,
"id": "38714df3-0cb4-4a56-92c9-67bb32b15ec3",
"metadata": {},
"outputs": [
{
"ename": "IndexError",
"evalue": "index 14996 is out of bounds for axis 0 with size 14996",
"output_type": "error",
"traceback": [
"\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
"\u001b[0;31mIndexError\u001b[0m Traceback (most recent call last)",
"\u001b[0;32m<ipython-input-19-d7a48b9d9945>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m 23\u001b[0m \u001b[0mmethod\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m'MSD'\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 24\u001b[0m \u001b[0;31m#method='WeakNoise',\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 25\u001b[0;31m \u001b[0mbasis\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m{\u001b[0m \u001b[0;34m'type'\u001b[0m \u001b[0;34m:\u001b[0m \u001b[0;34m'polynomial'\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m'order'\u001b[0m \u001b[0;34m:\u001b[0m \u001b[0;36m1\u001b[0m\u001b[0;34m}\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 26\u001b[0m ) \n\u001b[1;32m 27\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;32m~/projects/StochasticForceInference/SFI_inference.py\u001b[0m in \u001b[0;36mcompute_diffusion\u001b[0;34m(self, basis, method, space_dependent_error, regularize)\u001b[0m\n\u001b[1;32m 212\u001b[0m \u001b[0;32mraise\u001b[0m \u001b[0mKeyError\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\"Wrong diffusion_method argument:\"\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mdiffusion_method\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 213\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 214\u001b[0;31m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mD_average\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0meinsum\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m't,tmn->mn'\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdata\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdt\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0marray\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m[\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0meinsum\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'imn->mn'\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mD_local\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mt\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m)\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0mt\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mrange\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mlen\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdata\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mt\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m/\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdata\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mtauN\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 215\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mD_average_inv\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mlinalg\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0minv\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mD_average\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 216\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;32m~/projects/StochasticForceInference/SFI_inference.py\u001b[0m in \u001b[0;36m<listcomp>\u001b[0;34m(.0)\u001b[0m\n\u001b[1;32m 212\u001b[0m \u001b[0;32mraise\u001b[0m \u001b[0mKeyError\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\"Wrong diffusion_method argument:\"\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mdiffusion_method\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 213\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 214\u001b[0;31m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mD_average\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0meinsum\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m't,tmn->mn'\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdata\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdt\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0marray\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m[\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0meinsum\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'imn->mn'\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mD_local\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mt\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m)\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0mt\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mrange\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mlen\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdata\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mt\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m/\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdata\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mtauN\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 215\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mD_average_inv\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mlinalg\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0minv\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mD_average\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 216\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;32m~/projects/StochasticForceInference/SFI_inference.py\u001b[0m in \u001b[0;36m__D_MSD__\u001b[0;34m(self, t)\u001b[0m\n\u001b[1;32m 514\u001b[0m \u001b[0;31m# adapted to the problem at hand.\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 515\u001b[0m \u001b[0;32mdef\u001b[0m \u001b[0m__D_MSD__\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mt\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 516\u001b[0;31m \u001b[0;32mreturn\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0meinsum\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'im,in->imn'\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdata\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdX_plus\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mt\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdata\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdX_plus\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mt\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m/\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdata\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdt\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mt\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 517\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 518\u001b[0m \u001b[0;32mdef\u001b[0m \u001b[0m__D_Vestergaard__\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mt\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;31mIndexError\u001b[0m: index 14996 is out of bounds for axis 0 with size 14996"
]
}
],
"source": [
"# do inference\n",
"freq = 1\n",
"\n",
"# data = StochasticTrajectoryData(Xcat[::freq],tlist[::freq]) \n",
"\n",
"# center = data.X_ito.mean(axis=(0,1)) \n",
"# width = 2.1 * abs(data.X_ito-center).max(axis=(0,1)) \n",
"\n",
"S = StochasticForceInference(hetero) \n",
"\n",
"S.compute_drift(basis = { 'type' : 'polynomial', 'order' : 2} ,\n",
" #basis = { 'type' : 'Fourier', 'order' : 3, 'center' : center, 'width' : width, } ,\n",
" #diffusion_mode = 'WeakNoise', # Best for space-dependent noise with large dt\n",
" diffusion_mode = 'MSD', # Best for space-dependent noise with short trajectories\n",
" #diffusion_mode = 'constant', \n",
" #diffusion_mode = 'Vestergaard', # Best for space-dependent noise with large measurement error \n",
" #mode='Ito'\n",
") \n",
"\n",
"\n",
"S.compute_diffusion(\n",
" #method='Vestergaard',\n",
" method='MSD',\n",
" #method='WeakNoise',\n",
" basis = { 'type' : 'polynomial', 'order' : 1}\n",
") \n",
"\n",
"S.compute_force()\n",
"S.compute_drift_error() \n",
"S.compute_diffusion_error()\n",
"S.compute_entropy()\n",
"\n",
"S.print_report()\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "2482d1e9-adb7-4d4d-92a3-339cbfc25de8",
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "sfi",
"language": "python",
"name": "sfi"
},
"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.6.13"
}
},
"nbformat": 4,
"nbformat_minor": 5
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment