jakeyeung/sfi_heterogeneous_data_question.ipynb

## sfi_heterogeneous_data_question.ipynb
{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "f48c1e1a-82ec-4855-abcf-03e0ab359fc1",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7f791c31c048>]"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAD4CAYAAAD1jb0+AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAvqklEQVR4nO3dd3gU1RrA4d9JrySEBBJKSELvCKGqSBUQBAS7XpGLYle8NuxdsV+xXmyAoKKIiqJUQUARCFUIvYdAEghJSC977h+zpJAEUjaZ7O73Pg/PzpyZnfmOho/JmVOU1hohhBD2x8XsAIQQQlSNJHAhhLBTksCFEMJOSQIXQgg7JQlcCCHslFtt3iw4OFhHRETU5i2FEMLubdy48aTWOuTc8lpN4BEREcTExNTmLYUQwu4ppQ6XVS5NKEIIYackgQshhJ2SBC6EEHZKErgQQtgpSeBCCGGnJIELIYSdkgQuhBB2yr4T+F/vw4GVZkchhBCmqNWBPDaVnwMxn0PyfojsB5c+BFH9zY5KCCFqjf0+gbt5wl1/wZAX4eQ+mDUGEnaYHZUQQtQa+03gAO5ecPH9cNef4O4NP94NacfNjkoIIWqFfSfws3yCYNircHIPfHwJbPgMZKk4IYSDc4wEDtD9Vvj3YghpAwv/AwdXmR2REELUKMdJ4ABhneGmeeDhB/98a3Y0QghRoxwrgQN4+ECHq2D7D5BzxuxohBCixjheAgfoeiPkZcD+FWZHIoQQNcYxE3hIW0DB1m+M/uJCCOGAHDOB+wTBkOdh90L4drzZ0QghRI24YAJXSrVRSm0p9idNKTVZKRWklFqqlNpr/axfGwFX2MUPQL9HYc9vcPqQ2dEIIYTNXTCBa613a627aq27At2BTOAHYAqwXGvdClhu3a9bLroZXNxgxkjYu8zsaIQQwqYq24QyCNivtT4MjAZmWstnAmNsGJdt1G8O4382uhV+dS0sfhJS48yOSgghbKKyCfx64GvrdiOt9XEA62fDsr6glJqklIpRSsUkJSVVPdKqat4XblsG7UbCuo/h/Z5wdEPtxyGEEDZW4QSulPIARgHfVeYGWuvpWutorXV0SEhIZeOzDU8/uHYW3LcJ8rNh5avmxCGEEDZUmSfw4cAmrXWCdT9BKRUGYP1MtHVwNle/OQQ0hf3L4fBfZkcjhBDVUpkEfgNFzScAC4CzffTGAz/ZKqgadfN8QMHPD0BavNnRCCFElVUogSulfIAhwPxixVOBIUqpvdZjU20fXg0Ibgm3/Gi8zPywN+xbbnZEQghRJRVK4FrrTK11A611arGyU1rrQVrrVtbP5JoL08ai+sOkPyA7FWJ/NDsaIYSoEscciVkRDVqCbwgcWQf7fzc7GiGEqDTnTeAuLjDibcjLhC/Hwt8fmR2REEJUivMmcID2o+DeDdB2BCyaYqzkI4QQdsK5EzgYa2leMxNaDISlz0BeltkRCSFEhbiZHUCd4OpmtIfnphu9Ujx8jRkNw7qYHZkQQpRLEvhZXW4wXmbOvcnY96xnrHYfGG5uXEIIUQ5pQjmrxQB4MBaumw1jPjJWtZ9zDZzab3ZkQghRJkngxbl5QLsrjSXZrpsF6Ynw2RDYPt9I6EIIUYdIAi9Pi4HGLIb1GsO8CdLNUAhR50gCP58GLYwRm62Hw+LHYcUr8iQuhKgzJIFfiIsrjJpmbP/xmozaFELUGZLAK8KvIVz7pbF9bKO5sQghhJUk8IqKugwCwmHzl5CbYXY0QgghCbzCvAJg8LOQcgSOrjc7GiGEkAReKaGdAAU75l/wVCGEqGmSwCsjpA30mAibZsFX18O+ZWZHJIRwYhVdkSdQKTVPKbVLKbVTKdVHKRWklFqqlNpr/axf08HWCUNfhX6PGC8zZ4+DPYvNjkgI4aQq+gT+LrBIa90W6ALsBKYAy7XWrYDl1n3H5+YBA5+CB7dDUAtY9LjRLi6EELXsgglcKVUP6Ad8BqC1ztVapwCjgZnW02YCY2omxDrKzROueB2S98N/O8mcKUKIWleRJ/AoIAn4Qim1WSn1qVLKF2iktT4OYP1sWNaXlVKTlFIxSqmYpKQkmwVeJ7QcbEx+BbDmHXNjEUI4nYokcDegG/CR1voiIINKNJdoradrraO11tEhISFVDLMOa3cl9LnX6B8+80pI2mN2REIIJ1GRBB4HxGmt11n352Ek9ASlVBiA9TOxZkK0AwOehP6Pw/Ft8PElsHGGDPYRQtS4CyZwrfUJ4KhSqo21aBAQCywAxlvLxgM/1UiE9sDDB/pPgXvWQ3Ar+PkB+F8/s6MSQji4iq7Icx8wRynlARwAJmAk/2+VUhOBI8A1NROiHfFvBLevgC+GGd0MtQalzI5KCOGgKpTAtdZbgOgyDg2yaTSOwM0D/MOM7ReDYcIiaNbD3JiEEA5JRmLWhGtnQbfxYMmH1W+ZHY0QwkHJosY1wcUVrnwX8rNhx4+QlwXu3mZHJYRwMPIEXlOUgvZjoCAHXg6Fr28wOyIhhIORBF6T2gw3mlPC+8DuX2H9J2ZHJIRwIJLAa5JS0H40jP8FwrrArw9DZrLZUQkhHIQk8Nrg6gahnY3tE/+YG4sQwmFIAq8tnazd5NMTzI1DCOEwJIHXlqbRENzGGKV55oTZ0QghHIAk8Nri4Qu974S8TPjoYpj3b9j/uzFaUwghqkASeG3qdA30vtvoJ777N/jyKvjlQSjINzsyIYQNHT6Vwa1frOfwqQzOZOfV2H2UrsUnwOjoaB0TE1Nr96vT8nNgxcvw57vQ7RYY9Z7ZEQkhzspOhVVvQI/bwCcYPP0q9fWIKQtL7K9+dADNgnyqHI5SaqPWutR0JjIS0yxunjDkBcg5AzGfg08DGPi08XQuhDBP1ml4LcLY/sv6YHXZFBjweIW+brGUfihetjOBCRdH2ijAItKEYrbLrGtjrHlH+ogLYbasFPh5cunyP6aW+77qr/0niZiykIgpC9kRn8pHf5ReXvH5n2NtG6eVPIGbzb8RtBhovNCcfRVcOQ2adDM7KiGcg8UCR/4ytpv1hteaFx3r9wjZaafw2vI5AAUaXK2zQ2fnFTDzr0O8+tuuEpcbMW1Nubc6mpxZrWaUskgCrwtu+Aa2z4dlz8Gng6DPPTD4eWlOEaKmHd8MM0aUeWhvQF9arbqqcP+n39fg3aglLi6KO77cWOlbebnb/u+zJPC6wM0Tut5gzJ2y7Fmj3S0uBryDoMUA6DgOfILMjlIIxxPWtdxDrX6+qsT+2DUjGZXzItt0ixLlwaQS7bKbHbo5R3VDGpDGKQJKnDOkfSNC/D1tFvZZksDrEu9AYxparwCjdwrA7oXGHCpjPjaSvBDChiq3YtYCz6cZlfMi23UkrVQcX3pMpaFKKff8dZa2XJf7DEtjE8jIycfX07Ypt0LdCJVSh4AzQAGQr7WOVkoFAXOBCOAQcK3W+vT5riPdCCshMxlc3WHuzXBgJaBgwm/QvI/ZkQnhmA6shFmjbX7ZiOw5gGLmv3tyWeuQKl2jvG6ElemFMkBr3bXYRaYAy7XWrYDl1n1hKz5B4OkPt/wEj+wH7/rw7b8grvJtb0KIC8tvbvuFyH8v6MrZp/yuTQNtfv3qdCMcDcy0bs8ExlQ7GlE232C4fg64esJnQ+Dvj8yOSAjHkJkMm76E/Bw6P7+Eh3Lv5Im8iTa7/EDXLdzsupR3rutCgI+7za57VkUTuAaWKKU2KqUmWcsaaa2PA1g/G5b1RaXUJKVUjFIqJikpqfoRO6vmfeHuv6D1MFj0uPGiU+ZREaLq0pPg9UhYcC8H3+xPt/wtLLV056uCqq3VPijnjTLLX3L/gquS/ledSMtV0QR+sda6GzAcuEcpVeHfNbTW07XW0Vrr6JCQqrX/CCuvABj3qdFbZclTRp/V+M1mRyWEfdr5U+FmZHYssz1eZaPnnUx2m1ely+3XTfgo/8qyD/41rUbmPKpQAtdax1s/E4EfgJ5AglIqDMD6mWjz6ERpHj5w/VcQ2c+Yr2HDZzIZlhBV0fXmUkXuqoDJbvOrdLlDXjdyl9vP5R4/cHBfla57PhdM4EopX6WU/9lt4HJgO7AAGG89bTzwU9lXEDanFNw0zxjBuflLeLEBHNtkdlRC2Bd3LyJzvrL5Ze/InVyqLEt7sCc7oPTJ1VSRTomNgB+UUmfP/0prvUgptQH4Vik1ETgCXGPz6ET53DzhhrnwVmtj8p3lz8N1cyo9a5oQzkxraJM9g91et9rkeuNzH+MPSxciskv/w7CrTZmvCatFppO1d5YC+OM1Y+rLJt1h7CcQZPtZz4RwRCmZuXR9YSmt1VGWeD5WrWs9kTex8AVoVIgv39/ZFzdXRdzpLKJCfPF0q/pQ+vL6gUsCdxSxP8G3t0CHsXDNF2ZHI4RdiZiyEA/y2OM1/sInl+O+lst48/pu1UrU5bHFQB5Rl7UfDQ3bw44fYO6/YOkzxlzjQgg4vBZ+uBOSD5bZ/fbgK8PJpRr9tHvewXs396iR5H0+ksAdyYRfoe99cGg1/DkNZo4y2seFcGYJsfDFMNj6NUzrCq82hecCjD/JB+G5ANQL9emsSs/jfa61Hn3Qd6wqfWCjOb/1SgJ3JN714fIX4bFDcMUbEL/JaFbJzzE7MiHM88drJfdz04u2p3Ut3Fzg+XTh9oCct+iQ/VmJr0Vmz6bPE4tQYV2KFmI5y8v2PUwqQmYjdFQ9bzf6if/+Iiz8D1zxJrh7mx2VELWvWS+I/bFSX1nh+RAA7bM/x40C0vBh38tXFJ0w4HG4ZDIk7YLFT8H48vt/1yR5Andklz5ktI1vng0x8mJTOCGtYfWbhbu/BN/G2JznSp7T665yvz7N/T3S8GXbc0Nxcz0nXbp7Q+OLYMJCcDEnlUoCd2RKwTUzwb8xrHpdZjIUzkNr2PEjPB8ImacAeL39j9wbN5BNujXdsj/mtPbj1JUzYPhUiLyszMt0cTHaxVfvOVk7cVeSdCN0Bkl7YM7VkJ4I18+GloPNjkgI2/pnHnxvnUWw4zjY/n2Jw1/3WcDjK9JLfc2NfPZ53VLuZaOzP6J/9w68cXVnrIMZTSHdCJ1ZSGu4bTkERcHscfDTvWZHJITtbPuuKHlDqeQNcMPaUUx0XViqfJzr6jIv+VDunURlz+YkATw9sr2pyft8JIE7C78Q4+nb1cOYP2X/CpmOVjiGLXMqdNrjbl9Tz6tkyvuloHeZ535v6YfFmh73JtTd8RTSC8WZBEXBNTPgmxvhyzHQdiRcN9toKxfCniTuhHkTKWg5mKzcfC40A9D9ufcSPWA0Wwf3IK9A4+FmJOeIKQuJyP4KfzJ5ym02qyydee2ZZ3hg9UE2H02hV2QQnWtgJR1bkTZwZ5SeBGvegb8/gGGvQe87zY5IiMpZ9ITx81tMh+zP2OFV9mo6nbM/YcNL40qNlNx5PI3h7xrNKB/c2I0RncNqJt5qKq8NXJ7AnZFfCAx+FhK2w+8vgU8DiLgE6tXNH14hSjknea8u6EgG3hywhBLlcqLEsczRn7C2/dgyh7m3C6tH7AtD+XFzPJe1sb8FZySBOys3T2PY/ZyrYf5tRtmdayC0k7lxCXGu3Az4631Y+QrUj4AHtqJ9G6IyjDVk3s8fzZv51wFQcE8M1HcBd5/CpkGfC1zex8ONG3uF12AFao68xHRmrYbAv34o2v/4Elj9FqTFmxeTEOeK22Akb4DTh2DRE4XJGyAdX/a+PJxDU0fQqpE/ePg6zXsdSeDOrsVAeC7VWKbNzQuWvwD/7QyLn5Sl2oS5EnfC3qUwa3TJ8nOaT5La3YL7uaMknUSFa62UclVKbVZK/WLdD1JKLVVK7bV+1q+5MEWNazsCnkqAe2OMpL72fWOptrxssyMTzui3x+DD3kYT3wVc3qlpLQRUN1Xmn60HgJ3F9qcAy7XWrYDl1n1h74JbwXVfFu3/fL/MKy5q357F5R460elOIrNn0yV7OhHZX/HNxuO1GFjdUqEErpRqCowAPi1WPBqYad2eCYyxaWTCPG6eMOo9Y3vbXHi7gzFvshC1ZcyHZRb/VHAxvTf0Q+NCqrX39029mtdmZHVKRXuh/Bd4FPAvVtZIa30cQGt9XClV5oqdSqlJwCSA8HD7fNPrlDpfB9lpYMmHZc8av9L2vc/obugkL4iEiVxKr46ztMFNPHLscqKCfTmdmYtSirmTehsvLp3UBRO4UmokkKi13qiU6l/ZG2itpwPTwRjIU9nvC5O4eULfYnOmrHgZ9lp/rX14L/jZfoVtIQo17lpiUqo3867h/WMjAPj94f7mxVXHVKQJ5WJglFLqEPANMFApNRtIUEqFAVg/E8u/hLBrl0yGR4s1oSTsMC0U4QQSdsCLwSUmpZpZMNTEgOquCyZwrfXjWuumWusI4Hrgd631zcAC4OwSzuOBn2osSmE+Tz+4f7Ox/eUYSNxlajjCQZ3YDh/1LdpvM4KMuzZx5oLDcZxTdTpPTgWGKKX2AkOs+8KRBUXBDXONGQ0/7AXrPzE7IuEoUo8Ziwx/fHFR2fVfU3DdHDq8Iw8L5anUUHqt9UpgpXX7FDDI9iGJOq3NMLhnHUy7CH59GIJbQ1TZq5kIUa60eFj6DCTthosfMNaWPFd4b05llFyQ++mR7WspQPvgnMOXRPUERcHd68DD30ji2almRyTqmsxksBSUfcxSAG+3g3++gxPbjMUYVr1RdDyqv7EAiU8QCaklE/jNvaUnW3GSwEXVNGxrzC1+cg9MDZcXm6LImQR4PRI+LWfpvheCzv/9Ue9D02hy8y3c9OnfhcVrHhtQ5oyCzkwSuKi6FgPB19qd8PNhsPBh2DRLVvpxZnuXwVutje34Tcbnjh9hzrXw9Y0ln8oDwuH+LfBkQlHZ3X9DYDMAFmyNJy3bmI/nrykDaVpfXmSeS6aTFVXn4gL/2QnpJ+D722CD9aVm/QiI7GdqaKKWaA3rPgbPenDRTTBnXMnjzwWU3C/+9N12BARFWs8r2QxnsWj2JhZN4fDLtngm9Wthy8gdgiRwUT2ubhDQFP69yJg17sBKOLUPmvUGNw+zoxM1SWv4+FJI+MfY/+nuomMNWkJIW9j1S/nfv+jmMouLr5Jz1oSLI6sbrUOSBC5sp/sEI4H/8qAx9P6O1UZbubB/m2YZk5r1vANy0oyXlPWbFyXvc4V2Mt6RnNhe1DUw4lK44RtjTEEZtNZsOnKaR+dtKyxb+mA/fDzdnHa62AuRBC5sp8MYaJMIP94N2+cZfcUf2Go0qQj7lbQbFtxnbCcfLGoqA7hnA3zQA/zDoHE32L3QKB/wlPEZ2hEmb4f3usGl/yk3eQN0em4J6TlFc9DX93F36nlOKkISuLAtN0+4+jNI3g/xm+HdLtDtFuj3qPGX3FV+5OzOileKtjecM3jrgx7GZ/dbof8UyM00lkDzK7a+ZGAzeDrpvLfIyS8okbzfuLozXZoFVi9uJyB/m0TNuH2F8Wt33AbYPNvYBqjXBIa/brzAklkN7cPY6RD74/nPOdue7eFj/KmkxLSi/t6/PXAp7cLqVfoazkgalkTNUAq6j4fR7xtt4SPfMbodph2DuTcZvVayUsyOUlREdlrJfd9zVm8f95nxIrsaJszYULgtybvi5Alc1LzQjsafDmONaWn3/260kfs1hGGvmh2duJBca3c+rwB4MNZox07cafRCCW4FrqXn7q6Mk+k57EtML9y3WDQuLvLbWUVIAhe1xzsQrngD/plnDJ9uPQzyc6W7YV2UcsRoy/5kEORlGGV5WUUvIRu2s9mtbvlsfYn9jNx8/L2q94+Cs5AELmqfi/XHbtaoorL7txQN6hA17++PYclT0O8RyM8yVlpqfjG4ecHu3+CbG0p/Z/DzNrn1xsPJdGgcQFp2HrP+OsywjqHEHi9qpnnxl1hev7qLTe7l6JSuxWHP0dHROiYmptbuJ+ooiwX+/hCSDxiT9mengLuP0bYaFAWDnoEm3cyO0nFZLPBC/Yqf/+hB8LnA/CUVNOPPgzz3cyz3DmjJ+yv2lXlO60Z+LHlQZrgsTim1UWsdXapcErgwXeoxo6taXiYcXAWZJ6HTtTDybfCUfsA2t/YDWPwENOluvICc1rXs8+5Zb3T99Kr6S8WEtGyCfD04eDKD8CAf2j2z6LxT5bx5TReu7l69F6KOqLwELk0ownwBTWDMB8Z2dhqseQf+/C+kHIZLH4KWg8FFZqGzibwsWPqssR3cxmi2evIEvBxqNKe0GQ6fDIQ7VkFImyrf5tsNR5m6aBfJGbmljnVqEsA/x0pPQfzD3X25KLwSvxkIeQIXddTyF2D1W8a2XyiMeBPaXWnsJx8E/1Bw9zYvPnuTuNPok5+ZDMueLSp/Ih48fG16q+83xvHQd1vLPb7i4f4MeHNl4f7mp4cQ6OOOknEB5aryE7hSygtYBXhaz5+ntX5WKRUEzAUigEPAtVrr07YMWjixSx82Bv14+MIPd8Dcm425NHLOwPEt4NcIut4EPSZWuw9ypVkKjEmaTmw3FrMY8ITRw6au2vqN8d/wXG1HGu8ebCQ+JYvxn69nb7EugS9f1ZHUrDxu6tWc+ZviaNnQj8hgXzY8OZhTGTm0DZU+39VxwSdwZfyz6Ku1TldKuQNrgAeAsUCy1nqqUmoKUF9r/dj5riVP4KJK4rfA9PO81Hom+cJNLNmpsOtXo4ng+FajbfeDHtD8Eji8xjinw1XgXR+GTQVtMXrLuLob/Z3/+Q7ys2Hr3KLzi5u4DJr1qHIVK8RSYPy3WP8/oymk+cWw9WvjHzQw5iKJ3wQ+wTDuE2PgFMBrEZB1uqiOLQZCUAto3temo2Ejpiws3F71yAAaBXjKAgw2YpOXmEopH4wEfhcwC+ivtT6ulAoDVmqtz9toJglcVEtBntFGbsmDk3th5kij/OrPjYQc3sdISCe2gy6AsC5GX+afHzAScEW5eRsJ3JIH/o2N5J9bNDc1fe6FzFPG0l/Fn2wn/QGNu9qipkZTR+JOo/fHie3QvI8xodT+3yt+jctfhiNri6Z0fTIB3L1sE985CiyaFk/8CsBlrUOY+e+eNXIfZ1Wtl5hKKVdgI9AS+EBrvU4p1UhrfRzAmsQblvPdScAkgPBwWc9OVIOrO/g2MLb9Q6H3PfD3BzDv30XnNL7ImEQL4LIpkJtuJG+/RuAdBOjSC+i2Hw1XzzCeXrd9C2fijadY3xBjbnNLnpFQ6zWBy1807n1Wm+HGknJgfL+iCTw/B+JiYNs3xjwxXgHQbTxc9phxbPbYoifrc92/BfYtM57Ie0w0pu7NSIRBzxpP6PNvM85b8iQENAOvQKOZydX2A6aOJmdy+TuryMorWmnnjz3nn7hK2E5ln8ADgR+A+4A1WuvAYsdOa63P+wpZnsCFze1ZbCS/z4cWlbUYZDyJ71tm7If3gQm/FTUXaG00h7h5GXNbe9arelPCtu+KEmaDVnDPOqM5Z8NncGi1cY+2I4pewK54Bf54rWLXDgw3RkQCePgZT/zdJ0CrctaaPKsgH3bMN+bktuGIybMsFk2B1sQcOs0Nn/xd6virYztxQ095WLMlm/UDV0o9C2QAtyNNKKKuSI0zmhoi+xlPmi6uRgJNPmgkTxsNRCkl5wzsWlj2S8Li2o8xVqlZ/WZR2TUzjaf2Bi0htLPRNn9ojRF7q8uhxQA4vg2i6s6gli/XHuLpn8pfwPrRYW24u3/LWozIOVQ5gSulQoA8rXWKUsobWAK8BlwGnCr2EjNIa/3o+a4lCVw4rOzUoqaUs+7fDLELSnbbO+uOVUYbvUm01kQ+/mvh/utXdyY8yIeoYF8a1ivZTp6TX8A9czbTI6I+r/6269xLlfD2tV0Y200G4thaddrAw4CZ1nZwF+BbrfUvSqm1wLdKqYnAEeAam0YshD3xCjB6w+xbBl9da5QFRcElk6FZL/hiWMnzTUreBRbNP8dSSw2kKb6MGcD254dyKj2H1Kw85m2MY9nOBJbtTChxzq19IwA4k51PalYuy3YmEt28hn7TEWWSgTxC2FpBPigXcCljuv3cDEBVadGDqvpgxT6a1vcmyNeD93/fx7qDyYXHvr2jD3GnM/Fwc2FPQjrTlu8t8xqBPu6kZOaVKDs0dUSJ/fwCC26ydmWNkKH0QtSW8y0bZ+NRjxey6chp3li8u8xj397Rh56RQfSMNJ6aLRZNXoGFj1buL3HeV7f3ok9UA7SGzLwCtsWl0K2MIe+SvGufJHAhHNj3G+MKt+8b2JKbejUn0Med+JQsokJKLjDs4qJ4bFhbHhvWtsxrKQV+nm70bRFcozGLipMELoSDOp6axZx1Rwr37xnQEi93Y2Tkuclb2CdJ4EI4qB82HwPghdEd+Ffv5jJZlAOSBC6Eg/niz4M8/3MsAC1CfLmuRzNJ3g5KErgQDiInv4D//XGAt5fuASDYz4PFk/vJy0UHJglcCAeQlp3HFe+uJu50FgAPDWnNvQNbypO3g5MELoSd01oz/vP1hcl789NDqO9r+4mrRN0jCVwIO/X3gVO8+utOtsYVjap8+9oukrydiCRwIezQit2JTPhiQ+H+61d3ZnTXxrKAgpORBC6EHTmanMknqw8wa+1hAJY+2I/wBj6SuJ2UJHAh7ETMoWSu/nht4f7//tWdVo38TYxImE36FwlhJ3bEpxVu3z+oFUM7hJ7nbOEMJIELYQc2Hj7NV8WGxU9bvpfUc2YHFM5HmlCEqMMenLuFXSfOsPN4WonyJoHeeHnI85ezkwQuRB2TmZvPrV9sYH2xebsBHhlqrFjo6qK4o1+UDNIRF07gSqlmwCwgFLAA07XW7yqlgoC5QARwCLhWa3265kIVwjmk5+QXJu+BbRtyXY9mBHq70yuqgcmRibqmIk/g+cBDWutNSil/YKNSailwK7C82JqYU4DHai5UIRzfbTNj+OdYSuF+x8b15GWlKNcFG9G01se11pus22eAnUATYDQw03raTGBMDcUohFPYl5jOsp0JJKTlFJZFhtTuCj7CvlSqDVwpFQFcBKwDGmmtj4OR5JVSDcv5ziRgEkB4eHhZpwghMOY0AXh1bCcGt2tEiL+nyRGJuq7Cr7GVUn7A98BkrXXahc4/S2s9XWsdrbWODgkJqUqMQjiFbdY5Tfw83SR5iwqpUAJXSrljJO85Wuv51uIEpVSY9XgYkFgzIQrh+A6dzOC37ScAmLX2kLnBCLtRkV4oCvgM2Km1frvYoQXAeGCq9fOnGolQCAegtWblniT+PnCKX/85TpNAb3pGNuDu/i344s9DvLZoV+G5V3dvamKkwp6os+1u5Z6g1CXAauAfjG6EAE9gtIN/C4QDR4BrtNbJZV7EKjo6WsfExFQ3ZiHqJK01E2ZsYOXuJPa8NBwPNxd2Hk9j+LurK/T9R4a24e7+LaR/tyhFKbVRax19bvkFn8C11muA8n6iBlU3MCEcwY74VEZMW1O4P/jtP8jJLyjRo2RQ24bEnc6ivq87n9/ag/bPLC48tuHJwdLuLSpNRmIKUU05+QUlkjfAkeTMwu2XxnTkpl7hpZ6sv7uzD9+sP8p9A1tK8hZVIglciGp4/ucdfPHnoRJlX93eiz5RDdiflEF2XgHtw+qV2SzSIyKIHhFBtRSpcESSwIWoggKLpsUTv5YoW/XIABrW88TL3VhcoWVDPzNCE05EErgQlXDoZAY/bD7G/M1xhWVdmgUy49YeshalqHWSwIWooF0n0hj2X6NHSbCfJ2O7NaF9WD0ubx8qyVuYQhK4EBWQV2DhyveKXlS+d8NF9GkhswMKc0kCF+ICLBbNrV+sJ69A0zbUn0WT+5kdkhCAJHAhzmt/UjpvLt7Nn/tOUc/LjbsHtDQ7JCEKSQIXohwbD5/mls/WkZFbwOB2jfjklu4ySlLUKZLAhTjHByv2MX9THAlpOYT4e/L5uM50aRYoyVvUOZLAhSgmO6+ANxbvBox+3LP+3ZPGgd4mRyVE2SSBC2EVn5JF36m/A3B9j2ZMHdfZ5IiEOL8KL+gghD3LyS/g/d/3si8xvdxzYuOL1il5dWyn2ghLiGqRBC4cXn6Bhbtnb+LNJXsY/u4qjpzKLHW8wKLZdsxYEadFiK+0dwu7IE0owuEtjU1g+a5EhncM5bftJxj53mqujW7Gw0PbkJGTT/eXlpU4f0TnxiZFKkTlSAIXDuuXbfHM3XCUvw+cItjPk6njOtOhcT3eXLKHT9ccZHHsCY4mZ5X4ziND23DnZS1MiliIypEELhxC3OlMxn74F12bBdIzMojeUQ2496vNhcf/M6Q1Ad7u3DuwFZd3CGXl7kR+236CdqH1aFTPiwkXRxAVIrMHCvtSkSXVPgdGAola647WsiBgLhABHAKu1VqfvtDNZEk1URPSsvOIfmkZufmWUsceGdqG63o0I9hPFkwQ9qu8JdUq8hJzBjDsnLIpwHKtdStguXVfiFo3f1McvV9ZTm6+BU83Fw68cgV+nsYvlnNu68WkflGSvIXDqsiamKuUUhHnFI8G+lu3ZwIrgcdsGZgQZckvsBB3OouDpzLYdjSVd5btISrYl4mXRtKpSQAuLortzw81O0whakVV28Abaa2PA2itjyulGpZ3olJqEjAJIDw8vIq3EwJOpecwee4WVu89WaL83esvolPTAJOiEsI8Nf4SU2s9HZgORht4Td9P2D+tNTn5FrzcXSmwaCbM2MCqPUmFxz1cXejSLIANh4zXLpK8hbOqagJPUEqFWZ++w4BEWwYlnNey2ARum2W86Pb3dMPT3YWT6bmFx2/sFc4rV8koSSGg6gl8ATAemGr9/MlmEQmn9ds/x7lrzqbC/XHdm5KSmcuQ9qEM7xhKalaeLF0mRDEXTOBKqa8xXlgGK6XigGcxEve3SqmJwBHgmpoMUjiG/UnprNydxIu/xNI21J9Px0ezJ+EMm4+k8N7v+wrPmzGhB5e1Dik1nF2StxAlVaQXyg3lHBpk41iEg0rPyWfa8r1MX3WgsGzXiTNc8tqKUucGeLvTv02578SFEMXISExR415euJOv1x8B4L6BLZl4SSR/7jvFlqOnGdohlE9WH2DxjgT+96/u9GsVYnK0QtgPSeDCpnLzLXi4uZCTX8Da/adYdzC5MHkvmnwpbUPrATCicxgjOocBEB0RZFq8QtgzSeDCJnLzLXy/KY7H5/9D12aBbDmaAoBSMKhtQ96/sRveHq7mBimEg5EELsqVnJHLByv2cfulUYQGeGGxaNYeOMXHf+znuVEdWLErkQBvdxZsjS8xuGbL0RTah9WjS7MAnriiHf5e7ibWQgjHJQlclGlbXAq3frGB5IxcZq09RIifJyfSsrFYh2INeuuPUt/pGRlE+7B63NgrnNaN/Gs5YiGcjyRwAUBegYW1+0/h7eHKN+uPsjT2BFrDqC6NOZmeQ31fD5oH+ZB0JoelOxNIycwr/O7Atg2ZOq4TDf29TKyBEM5HErjghZ9j+fzPgyXKBrdryAODWsswdSHqMEngTmrn8TSWxiZwJjuvRPIe3K4RT49sR/MGviZGJ4SoCEngTuDgyQz8vdyIOZTMgq3xrNpzkvScfAC83V0J8vXgiSvacXX3piZHKoSoDEngDujsKktKKaav2s8rv+4qcdzNRdGmkT9Xd2/K7f2izAhRCGEDksAdTGpmHl1eWAJAl6YBbI1L5bLWIbQJ9ae+jwe39GmOr6f8bxfCEcjfZAfz7ILtJfanDG/LrX0j8HKXQTRCOBpJ4A4iIyef+JQsftwSD8ClrYL5cmIvk6MSQtQkSeB11N6EM+xOOEPPyCAsFvD2cGXVniR2nUijU5MA+rdpiKebC0oplu9M4D/fbiU1y+ibPaJTGB/c1M3kGgghapok8FqSV2BhR3waPh6u5Y5STDyTzYNztxDo48HCbccrdf22of48P6oDJ9NzGNm5sS1CFkLUcZLAa5jWGqUUrZ78rbCsb4sGREcE0TsqiA6NA4iNT2NJ7Am+i4kjPScfX+ukTx2b1GNcN6Nr3+nMPHpGBOHv5cbyXYlMW74XgAFtQmjg58mLozvKZFFCOBl1tstZlb6s1DDgXcAV+FRrPfV850dHR+uYmJgq36+uS0zLJj0nn4Fv/cHQDo1IzsgtXHj3rEtbBbM3IZ0Tadmlvj+4XUPu6t+Cjk0C0JoLvnjMziuQl5NCOAGl1EatdfS55VV+AldKuQIfAEOAOGCDUmqB1jq26mHarz0JZxj1/hqy8ywALN6RUOL4ub1BjqVk8eGKfWyNS+HGns3p3yaExoHelbqnJG8hnFt1mlB6Avu01gcAlFLfAKMBh0/giWnZbDpyGouG1Kw8Hp//D2DMfX1Jy2Aign24/dIoEtJyiArxJSu3gGZBPiWu0STQm5dldXUhRDVUJ4E3AY4W248DSvVbU0pNAiYBhIeHV+N25oo7ncnbS/fw9/5TxKeWbv4AmDGhJ5e1LloSTOYTEULUpOokcFVGWakGda31dGA6GG3g1bhfrbNYNP9bdYDVe5PYcCiZfItmaPtQJl4aRHxKFt2b16ddWD0aB3qRk2+hnixcIISoRdVJ4HFAs2L7TYH46oVjjt0nznAmOw9vD1fyCzRr9p1kyY4T3HlZC15bZMwjctslkfyrT/Nyn6o93aQ9WghRu6qTwDcArZRSkcAx4HrgRptEVYPyCyzEp2Tz0sJYLgqvT30fd578cTsFltK/HNw1ZxMAjw1ry139W9R2qEIIcV5VTuBa63yl1L3AYoxuhJ9rrXfYLLJi9iSc4URqNv2KtS8Xl5iWzYcr95OWlUfTIB82Hk6mga8nmbn5NA70pk2oPydSszmWksX8TccKv7ck1ugpEhbgxeTBrfD1dMPb3ZWkMzm4ubqQlVeAr4crwzuG1US1hBCiWqo1kEdr/Svwq41iKddHK/fzw+ZjjOrSmKdGtCPQx4N5G+OYtnwvGbn5nMnOL3F+h8b1OJCUga+nG38fSCY9Jx9XF4WXmwsAQ9o34pruTenUNIA9Cen0jAiSQTBCCLtjFyMxXx3bieYNfPhwxX4WbC1qZu/SLJAuTQNo3sCXi8IDaRroTUpWXomh6haL5lhKFg38PPDxKF3dsIDK9b0WQoi6wi4SuJe7K5MHt2ZUl8Z8+fdhUjPz6BEZxHXRzXBxKdkZpmG9kgvrurioUn2whRDCEdhFAj8rKsSPZ6/sYHYYQghRJ7iYHYAQQoiqkQQuhBB2ShK4EELYKUngQghhpySBCyGEnZIELoQQdkoSuBBC2ClJ4EIIYaeqtSZmpW+mVBJw2MaXDQZO2viadY3U0TFIHR2DGXVsrrUuNZtfrSbwmqCUiilrsU9HInV0DFJHx1CX6ihNKEIIYackgQshhJ1yhAQ+3ewAaoHU0TFIHR1Dnamj3beBCyGEs3KEJ3AhhHBKksCFEMJO2W0CV0oNU0rtVkrtU0pNMTseW1BKNVNKrVBK7VRK7VBKPWAtD1JKLVVK7bV+1jc71upSSrkqpTYrpX6x7jtUHZVSgUqpeUqpXdb/n30csI4PWn9OtyulvlZKeTlCHZVSnyulEpVS24uVlVsvpdTj1jy0Wyk1tDZjtcsErpRyBT4AhgPtgRuUUu3Njcom8oGHtNbtgN7APdZ6TQGWa61bAcut+/buAWBnsX1Hq+O7wCKtdVugC0ZdHaaOSqkmwP1AtNa6I+AKXI9j1HEGMOycsjLrZf37eT3QwfqdD635qVbYZQIHegL7tNYHtNa5wDfAaJNjqjat9XGt9Sbr9hmMv/RNMOo203raTGCMKQHaiFKqKTAC+LRYscPUUSlVD+gHfAagtc7VWqfgQHW0cgO8lVJugA8QjwPUUWu9Ckg+p7i8eo0GvtFa52itDwL7MPJTrbDXBN4EOFpsP85a5jCUUhHARcA6oJHW+jgYSR5oaGJotvBf4FHAUqzMkeoYBSQBX1ibiT5VSvniQHXUWh8D3gSOAMeBVK31Ehyojucor16m5iJ7TeCqjDKH6Q+plPIDvgcma63TzI7HlpRSI4FErfVGs2OpQW5AN+AjrfVFQAb22ZRQLmsb8GggEmgM+CqlbjY3KlOYmovsNYHHAc2K7TfF+PXN7iml3DGS9xyt9XxrcYJSKsx6PAxINCs+G7gYGKWUOoTR9DVQKTUbx6pjHBCntV5n3Z+HkdAdqY6DgYNa6yStdR4wH+iLY9WxuPLqZWoustcEvgFopZSKVEp5YLxEWGByTNWmlFIY7aY7tdZvFzu0ABhv3R4P/FTbsdmK1vpxrXVTrXUExv+337XWN+NYdTwBHFVKtbEWDQJicaA6YjSd9FZK+Vh/bgdhvLNxpDoWV169FgDXK6U8lVKRQCtgfa1FpbW2yz/AFcAeYD/wpNnx2KhOl2D8+rUN2GL9cwXQAOPN917rZ5DZsdqovv2BX6zbDlVHoCsQY/1/+SNQ3wHr+DywC9gOfAl4OkIdga8x2vXzMJ6wJ56vXsCT1jy0Gxhem7HKUHohhLBT9tqEIoQQTk8SuBBC2ClJ4EIIYackgQshhJ2SBC6EEHZKErgQQtgpSeBCCGGn/g9OURSfyGybYQAAAABJRU5ErkJggg==\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[1::2]\n",
    "tlist2 = tlist[::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",
    "\n",
    "\n",
    "# how to wrangle this data into StochasticTrajectoryData  ?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b2622bd8-2801-4a29-b601-102c5b496edb",
   "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
}
	{
	"cells": [
	{
	"cell_type": "code",
	"execution_count": 1,
	"id": "f48c1e1a-82ec-4855-abcf-03e0ab359fc1",
	"metadata": {},
	"outputs": [
	{
	"data": {
	"text/plain": [
	"[<matplotlib.lines.Line2D at 0x7f791c31c048>]"
	]
	},
	"execution_count": 1,
	"metadata": {},
	"output_type": "execute_result"
	},
	{
	"data": {
	"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAD4CAYAAAD1jb0+AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAvqklEQVR4nO3dd3gU1RrA4d9JrySEBBJKSELvCKGqSBUQBAS7XpGLYle8NuxdsV+xXmyAoKKIiqJUQUARCFUIvYdAEghJSC977h+zpJAEUjaZ7O73Pg/PzpyZnfmOho/JmVOU1hohhBD2x8XsAIQQQlSNJHAhhLBTksCFEMJOSQIXQgg7JQlcCCHslFtt3iw4OFhHRETU5i2FEMLubdy48aTWOuTc8lpN4BEREcTExNTmLYUQwu4ppQ6XVS5NKEIIYackgQshhJ2SBC6EEHZKErgQQtgpSeBCCGGnJIELIYSdkgQuhBB2yr4T+F/vw4GVZkchhBCmqNWBPDaVnwMxn0PyfojsB5c+BFH9zY5KCCFqjf0+gbt5wl1/wZAX4eQ+mDUGEnaYHZUQQtQa+03gAO5ecPH9cNef4O4NP94NacfNjkoIIWqFfSfws3yCYNircHIPfHwJbPgMZKk4IYSDc4wEDtD9Vvj3YghpAwv/AwdXmR2REELUKMdJ4ABhneGmeeDhB/98a3Y0QghRoxwrgQN4+ECHq2D7D5BzxuxohBCixjheAgfoeiPkZcD+FWZHIoQQNcYxE3hIW0DB1m+M/uJCCOGAHDOB+wTBkOdh90L4drzZ0QghRI24YAJXSrVRSm0p9idNKTVZKRWklFqqlNpr/axfGwFX2MUPQL9HYc9vcPqQ2dEIIYTNXTCBa613a627aq27At2BTOAHYAqwXGvdClhu3a9bLroZXNxgxkjYu8zsaIQQwqYq24QyCNivtT4MjAZmWstnAmNsGJdt1G8O4382uhV+dS0sfhJS48yOSgghbKKyCfx64GvrdiOt9XEA62fDsr6glJqklIpRSsUkJSVVPdKqat4XblsG7UbCuo/h/Z5wdEPtxyGEEDZW4QSulPIARgHfVeYGWuvpWutorXV0SEhIZeOzDU8/uHYW3LcJ8rNh5avmxCGEEDZUmSfw4cAmrXWCdT9BKRUGYP1MtHVwNle/OQQ0hf3L4fBfZkcjhBDVUpkEfgNFzScAC4CzffTGAz/ZKqgadfN8QMHPD0BavNnRCCFElVUogSulfIAhwPxixVOBIUqpvdZjU20fXg0Ibgm3/Gi8zPywN+xbbnZEQghRJRVK4FrrTK11A611arGyU1rrQVrrVtbP5JoL08ai+sOkPyA7FWJ/NDsaIYSoEscciVkRDVqCbwgcWQf7fzc7GiGEqDTnTeAuLjDibcjLhC/Hwt8fmR2REEJUivMmcID2o+DeDdB2BCyaYqzkI4QQdsK5EzgYa2leMxNaDISlz0BeltkRCSFEhbiZHUCd4OpmtIfnphu9Ujx8jRkNw7qYHZkQQpRLEvhZXW4wXmbOvcnY96xnrHYfGG5uXEIIUQ5pQjmrxQB4MBaumw1jPjJWtZ9zDZzab3ZkQghRJkngxbl5QLsrjSXZrpsF6Ynw2RDYPt9I6EIIUYdIAi9Pi4HGLIb1GsO8CdLNUAhR50gCP58GLYwRm62Hw+LHYcUr8iQuhKgzJIFfiIsrjJpmbP/xmozaFELUGZLAK8KvIVz7pbF9bKO5sQghhJUk8IqKugwCwmHzl5CbYXY0QgghCbzCvAJg8LOQcgSOrjc7GiGEkAReKaGdAAU75l/wVCGEqGmSwCsjpA30mAibZsFX18O+ZWZHJIRwYhVdkSdQKTVPKbVLKbVTKdVHKRWklFqqlNpr/axf08HWCUNfhX6PGC8zZ4+DPYvNjkgI4aQq+gT+LrBIa90W6ALsBKYAy7XWrYDl1n3H5+YBA5+CB7dDUAtY9LjRLi6EELXsgglcKVUP6Ad8BqC1ztVapwCjgZnW02YCY2omxDrKzROueB2S98N/O8mcKUKIWleRJ/AoIAn4Qim1WSn1qVLKF2iktT4OYP1sWNaXlVKTlFIxSqmYpKQkmwVeJ7QcbEx+BbDmHXNjEUI4nYokcDegG/CR1voiIINKNJdoradrraO11tEhISFVDLMOa3cl9LnX6B8+80pI2mN2REIIJ1GRBB4HxGmt11n352Ek9ASlVBiA9TOxZkK0AwOehP6Pw/Ft8PElsHGGDPYRQtS4CyZwrfUJ4KhSqo21aBAQCywAxlvLxgM/1UiE9sDDB/pPgXvWQ3Ar+PkB+F8/s6MSQji4iq7Icx8wRynlARwAJmAk/2+VUhOBI8A1NROiHfFvBLevgC+GGd0MtQalzI5KCOGgKpTAtdZbgOgyDg2yaTSOwM0D/MOM7ReDYcIiaNbD3JiEEA5JRmLWhGtnQbfxYMmH1W+ZHY0QwkHJosY1wcUVrnwX8rNhx4+QlwXu3mZHJYRwMPIEXlOUgvZjoCAHXg6Fr28wOyIhhIORBF6T2gw3mlPC+8DuX2H9J2ZHJIRwIJLAa5JS0H40jP8FwrrArw9DZrLZUQkhHIQk8Nrg6gahnY3tE/+YG4sQwmFIAq8tnazd5NMTzI1DCOEwJIHXlqbRENzGGKV55oTZ0QghHIAk8Nri4Qu974S8TPjoYpj3b9j/uzFaUwghqkASeG3qdA30vtvoJ777N/jyKvjlQSjINzsyIYQNHT6Vwa1frOfwqQzOZOfV2H2UrsUnwOjoaB0TE1Nr96vT8nNgxcvw57vQ7RYY9Z7ZEQkhzspOhVVvQI/bwCcYPP0q9fWIKQtL7K9+dADNgnyqHI5SaqPWutR0JjIS0yxunjDkBcg5AzGfg08DGPi08XQuhDBP1ml4LcLY/sv6YHXZFBjweIW+brGUfihetjOBCRdH2ijAItKEYrbLrGtjrHlH+ogLYbasFPh5cunyP6aW+77qr/0niZiykIgpC9kRn8pHf5ReXvH5n2NtG6eVPIGbzb8RtBhovNCcfRVcOQ2adDM7KiGcg8UCR/4ytpv1hteaFx3r9wjZaafw2vI5AAUaXK2zQ2fnFTDzr0O8+tuuEpcbMW1Nubc6mpxZrWaUskgCrwtu+Aa2z4dlz8Gng6DPPTD4eWlOEaKmHd8MM0aUeWhvQF9arbqqcP+n39fg3aglLi6KO77cWOlbebnb/u+zJPC6wM0Tut5gzJ2y7Fmj3S0uBryDoMUA6DgOfILMjlIIxxPWtdxDrX6+qsT+2DUjGZXzItt0ixLlwaQS7bKbHbo5R3VDGpDGKQJKnDOkfSNC/D1tFvZZksDrEu9AYxparwCjdwrA7oXGHCpjPjaSvBDChiq3YtYCz6cZlfMi23UkrVQcX3pMpaFKKff8dZa2XJf7DEtjE8jIycfX07Ypt0LdCJVSh4AzQAGQr7WOVkoFAXOBCOAQcK3W+vT5riPdCCshMxlc3WHuzXBgJaBgwm/QvI/ZkQnhmA6shFmjbX7ZiOw5gGLmv3tyWeuQKl2jvG6ElemFMkBr3bXYRaYAy7XWrYDl1n1hKz5B4OkPt/wEj+wH7/rw7b8grvJtb0KIC8tvbvuFyH8v6MrZp/yuTQNtfv3qdCMcDcy0bs8ExlQ7GlE232C4fg64esJnQ+Dvj8yOSAjHkJkMm76E/Bw6P7+Eh3Lv5Im8iTa7/EDXLdzsupR3rutCgI+7za57VkUTuAaWKKU2KqUmWcsaaa2PA1g/G5b1RaXUJKVUjFIqJikpqfoRO6vmfeHuv6D1MFj0uPGiU+ZREaLq0pPg9UhYcC8H3+xPt/wtLLV056uCqq3VPijnjTLLX3L/gquS/ledSMtV0QR+sda6GzAcuEcpVeHfNbTW07XW0Vrr6JCQqrX/CCuvABj3qdFbZclTRp/V+M1mRyWEfdr5U+FmZHYssz1eZaPnnUx2m1ely+3XTfgo/8qyD/41rUbmPKpQAtdax1s/E4EfgJ5AglIqDMD6mWjz6ERpHj5w/VcQ2c+Yr2HDZzIZlhBV0fXmUkXuqoDJbvOrdLlDXjdyl9vP5R4/cHBfla57PhdM4EopX6WU/9lt4HJgO7AAGG89bTzwU9lXEDanFNw0zxjBuflLeLEBHNtkdlRC2Bd3LyJzvrL5Ze/InVyqLEt7sCc7oPTJ1VSRTomNgB+UUmfP/0prvUgptQH4Vik1ETgCXGPz6ET53DzhhrnwVmtj8p3lz8N1cyo9a5oQzkxraJM9g91et9rkeuNzH+MPSxciskv/w7CrTZmvCatFppO1d5YC+OM1Y+rLJt1h7CcQZPtZz4RwRCmZuXR9YSmt1VGWeD5WrWs9kTex8AVoVIgv39/ZFzdXRdzpLKJCfPF0q/pQ+vL6gUsCdxSxP8G3t0CHsXDNF2ZHI4RdiZiyEA/y2OM1/sInl+O+lst48/pu1UrU5bHFQB5Rl7UfDQ3bw44fYO6/YOkzxlzjQgg4vBZ+uBOSD5bZ/fbgK8PJpRr9tHvewXs396iR5H0+ksAdyYRfoe99cGg1/DkNZo4y2seFcGYJsfDFMNj6NUzrCq82hecCjD/JB+G5ANQL9emsSs/jfa61Hn3Qd6wqfWCjOb/1SgJ3JN714fIX4bFDcMUbEL/JaFbJzzE7MiHM88drJfdz04u2p3Ut3Fzg+XTh9oCct+iQ/VmJr0Vmz6bPE4tQYV2KFmI5y8v2PUwqQmYjdFQ9bzf6if/+Iiz8D1zxJrh7mx2VELWvWS+I/bFSX1nh+RAA7bM/x40C0vBh38tXFJ0w4HG4ZDIk7YLFT8H48vt/1yR5Andklz5ktI1vng0x8mJTOCGtYfWbhbu/BN/G2JznSp7T665yvz7N/T3S8GXbc0Nxcz0nXbp7Q+OLYMJCcDEnlUoCd2RKwTUzwb8xrHpdZjIUzkNr2PEjPB8ImacAeL39j9wbN5BNujXdsj/mtPbj1JUzYPhUiLyszMt0cTHaxVfvOVk7cVeSdCN0Bkl7YM7VkJ4I18+GloPNjkgI2/pnHnxvnUWw4zjY/n2Jw1/3WcDjK9JLfc2NfPZ53VLuZaOzP6J/9w68cXVnrIMZTSHdCJ1ZSGu4bTkERcHscfDTvWZHJITtbPuuKHlDqeQNcMPaUUx0XViqfJzr6jIv+VDunURlz+YkATw9sr2pyft8JIE7C78Q4+nb1cOYP2X/CpmOVjiGLXMqdNrjbl9Tz6tkyvuloHeZ535v6YfFmh73JtTd8RTSC8WZBEXBNTPgmxvhyzHQdiRcN9toKxfCniTuhHkTKWg5mKzcfC40A9D9ufcSPWA0Wwf3IK9A4+FmJOeIKQuJyP4KfzJ5ym02qyydee2ZZ3hg9UE2H02hV2QQnWtgJR1bkTZwZ5SeBGvegb8/gGGvQe87zY5IiMpZ9ITx81tMh+zP2OFV9mo6nbM/YcNL40qNlNx5PI3h7xrNKB/c2I0RncNqJt5qKq8NXJ7AnZFfCAx+FhK2w+8vgU8DiLgE6tXNH14hSjknea8u6EgG3hywhBLlcqLEsczRn7C2/dgyh7m3C6tH7AtD+XFzPJe1sb8FZySBOys3T2PY/ZyrYf5tRtmdayC0k7lxCXGu3Az4631Y+QrUj4AHtqJ9G6IyjDVk3s8fzZv51wFQcE8M1HcBd5/CpkGfC1zex8ONG3uF12AFao68xHRmrYbAv34o2v/4Elj9FqTFmxeTEOeK22Akb4DTh2DRE4XJGyAdX/a+PJxDU0fQqpE/ePg6zXsdSeDOrsVAeC7VWKbNzQuWvwD/7QyLn5Sl2oS5EnfC3qUwa3TJ8nOaT5La3YL7uaMknUSFa62UclVKbVZK/WLdD1JKLVVK7bV+1q+5MEWNazsCnkqAe2OMpL72fWOptrxssyMTzui3x+DD3kYT3wVc3qlpLQRUN1Xmn60HgJ3F9qcAy7XWrYDl1n1h74JbwXVfFu3/fL/MKy5q357F5R460elOIrNn0yV7OhHZX/HNxuO1GFjdUqEErpRqCowAPi1WPBqYad2eCYyxaWTCPG6eMOo9Y3vbXHi7gzFvshC1ZcyHZRb/VHAxvTf0Q+NCqrX39029mtdmZHVKRXuh/Bd4FPAvVtZIa30cQGt9XClV5oqdSqlJwCSA8HD7fNPrlDpfB9lpYMmHZc8av9L2vc/obugkL4iEiVxKr46ztMFNPHLscqKCfTmdmYtSirmTehsvLp3UBRO4UmokkKi13qiU6l/ZG2itpwPTwRjIU9nvC5O4eULfYnOmrHgZ9lp/rX14L/jZfoVtIQo17lpiUqo3867h/WMjAPj94f7mxVXHVKQJ5WJglFLqEPANMFApNRtIUEqFAVg/E8u/hLBrl0yGR4s1oSTsMC0U4QQSdsCLwSUmpZpZMNTEgOquCyZwrfXjWuumWusI4Hrgd631zcAC4OwSzuOBn2osSmE+Tz+4f7Ox/eUYSNxlajjCQZ3YDh/1LdpvM4KMuzZx5oLDcZxTdTpPTgWGKKX2AkOs+8KRBUXBDXONGQ0/7AXrPzE7IuEoUo8Ziwx/fHFR2fVfU3DdHDq8Iw8L5anUUHqt9UpgpXX7FDDI9iGJOq3NMLhnHUy7CH59GIJbQ1TZq5kIUa60eFj6DCTthosfMNaWPFd4b05llFyQ++mR7WspQPvgnMOXRPUERcHd68DD30ji2almRyTqmsxksBSUfcxSAG+3g3++gxPbjMUYVr1RdDyqv7EAiU8QCaklE/jNvaUnW3GSwEXVNGxrzC1+cg9MDZcXm6LImQR4PRI+LWfpvheCzv/9Ue9D02hy8y3c9OnfhcVrHhtQ5oyCzkwSuKi6FgPB19qd8PNhsPBh2DRLVvpxZnuXwVutje34Tcbnjh9hzrXw9Y0ln8oDwuH+LfBkQlHZ3X9DYDMAFmyNJy3bmI/nrykDaVpfXmSeS6aTFVXn4gL/2QnpJ+D722CD9aVm/QiI7GdqaKKWaA3rPgbPenDRTTBnXMnjzwWU3C/+9N12BARFWs8r2QxnsWj2JhZN4fDLtngm9Wthy8gdgiRwUT2ubhDQFP69yJg17sBKOLUPmvUGNw+zoxM1SWv4+FJI+MfY/+nuomMNWkJIW9j1S/nfv+jmMouLr5Jz1oSLI6sbrUOSBC5sp/sEI4H/8qAx9P6O1UZbubB/m2YZk5r1vANy0oyXlPWbFyXvc4V2Mt6RnNhe1DUw4lK44RtjTEEZtNZsOnKaR+dtKyxb+mA/fDzdnHa62AuRBC5sp8MYaJMIP94N2+cZfcUf2Go0qQj7lbQbFtxnbCcfLGoqA7hnA3zQA/zDoHE32L3QKB/wlPEZ2hEmb4f3usGl/yk3eQN0em4J6TlFc9DX93F36nlOKkISuLAtN0+4+jNI3g/xm+HdLtDtFuj3qPGX3FV+5OzOileKtjecM3jrgx7GZ/dbof8UyM00lkDzK7a+ZGAzeDrpvLfIyS8okbzfuLozXZoFVi9uJyB/m0TNuH2F8Wt33AbYPNvYBqjXBIa/brzAklkN7cPY6RD74/nPOdue7eFj/KmkxLSi/t6/PXAp7cLqVfoazkgalkTNUAq6j4fR7xtt4SPfMbodph2DuTcZvVayUsyOUlREdlrJfd9zVm8f95nxIrsaJszYULgtybvi5Alc1LzQjsafDmONaWn3/260kfs1hGGvmh2duJBca3c+rwB4MNZox07cafRCCW4FrqXn7q6Mk+k57EtML9y3WDQuLvLbWUVIAhe1xzsQrngD/plnDJ9uPQzyc6W7YV2UcsRoy/5kEORlGGV5WUUvIRu2s9mtbvlsfYn9jNx8/L2q94+Cs5AELmqfi/XHbtaoorL7txQN6hA17++PYclT0O8RyM8yVlpqfjG4ecHu3+CbG0p/Z/DzNrn1xsPJdGgcQFp2HrP+OsywjqHEHi9qpnnxl1hev7qLTe7l6JSuxWHP0dHROiYmptbuJ+ooiwX+/hCSDxiT9mengLuP0bYaFAWDnoEm3cyO0nFZLPBC/Yqf/+hB8LnA/CUVNOPPgzz3cyz3DmjJ+yv2lXlO60Z+LHlQZrgsTim1UWsdXapcErgwXeoxo6taXiYcXAWZJ6HTtTDybfCUfsA2t/YDWPwENOluvICc1rXs8+5Zb3T99Kr6S8WEtGyCfD04eDKD8CAf2j2z6LxT5bx5TReu7l69F6KOqLwELk0ownwBTWDMB8Z2dhqseQf+/C+kHIZLH4KWg8FFZqGzibwsWPqssR3cxmi2evIEvBxqNKe0GQ6fDIQ7VkFImyrf5tsNR5m6aBfJGbmljnVqEsA/x0pPQfzD3X25KLwSvxkIeQIXddTyF2D1W8a2XyiMeBPaXWnsJx8E/1Bw9zYvPnuTuNPok5+ZDMueLSp/Ih48fG16q+83xvHQd1vLPb7i4f4MeHNl4f7mp4cQ6OOOknEB5aryE7hSygtYBXhaz5+ntX5WKRUEzAUigEPAtVrr07YMWjixSx82Bv14+MIPd8Dcm425NHLOwPEt4NcIut4EPSZWuw9ypVkKjEmaTmw3FrMY8ITRw6au2vqN8d/wXG1HGu8ebCQ+JYvxn69nb7EugS9f1ZHUrDxu6tWc+ZviaNnQj8hgXzY8OZhTGTm0DZU+39VxwSdwZfyz6Ku1TldKuQNrgAeAsUCy1nqqUmoKUF9r/dj5riVP4KJK4rfA9PO81Hom+cJNLNmpsOtXo4ng+FajbfeDHtD8Eji8xjinw1XgXR+GTQVtMXrLuLob/Z3/+Q7ys2Hr3KLzi5u4DJr1qHIVK8RSYPy3WP8/oymk+cWw9WvjHzQw5iKJ3wQ+wTDuE2PgFMBrEZB1uqiOLQZCUAto3temo2Ejpiws3F71yAAaBXjKAgw2YpOXmEopH4wEfhcwC+ivtT6ulAoDVmqtz9toJglcVEtBntFGbsmDk3th5kij/OrPjYQc3sdISCe2gy6AsC5GX+afHzAScEW5eRsJ3JIH/o2N5J9bNDc1fe6FzFPG0l/Fn2wn/QGNu9qipkZTR+JOo/fHie3QvI8xodT+3yt+jctfhiNri6Z0fTIB3L1sE985CiyaFk/8CsBlrUOY+e+eNXIfZ1Wtl5hKKVdgI9AS+EBrvU4p1UhrfRzAmsQblvPdScAkgPBwWc9OVIOrO/g2MLb9Q6H3PfD3BzDv30XnNL7ImEQL4LIpkJtuJG+/RuAdBOjSC+i2Hw1XzzCeXrd9C2fijadY3xBjbnNLnpFQ6zWBy1807n1Wm+HGknJgfL+iCTw/B+JiYNs3xjwxXgHQbTxc9phxbPbYoifrc92/BfYtM57Ie0w0pu7NSIRBzxpP6PNvM85b8iQENAOvQKOZydX2A6aOJmdy+TuryMorWmnnjz3nn7hK2E5ln8ADgR+A+4A1WuvAYsdOa63P+wpZnsCFze1ZbCS/z4cWlbUYZDyJ71tm7If3gQm/FTUXaG00h7h5GXNbe9arelPCtu+KEmaDVnDPOqM5Z8NncGi1cY+2I4pewK54Bf54rWLXDgw3RkQCePgZT/zdJ0CrctaaPKsgH3bMN+bktuGIybMsFk2B1sQcOs0Nn/xd6virYztxQ095WLMlm/UDV0o9C2QAtyNNKKKuSI0zmhoi+xlPmi6uRgJNPmgkTxsNRCkl5wzsWlj2S8Li2o8xVqlZ/WZR2TUzjaf2Bi0htLPRNn9ojRF7q8uhxQA4vg2i6s6gli/XHuLpn8pfwPrRYW24u3/LWozIOVQ5gSulQoA8rXWKUsobWAK8BlwGnCr2EjNIa/3o+a4lCVw4rOzUoqaUs+7fDLELSnbbO+uOVUYbvUm01kQ+/mvh/utXdyY8yIeoYF8a1ivZTp6TX8A9czbTI6I+r/6269xLlfD2tV0Y200G4thaddrAw4CZ1nZwF+BbrfUvSqm1wLdKqYnAEeAam0YshD3xCjB6w+xbBl9da5QFRcElk6FZL/hiWMnzTUreBRbNP8dSSw2kKb6MGcD254dyKj2H1Kw85m2MY9nOBJbtTChxzq19IwA4k51PalYuy3YmEt28hn7TEWWSgTxC2FpBPigXcCljuv3cDEBVadGDqvpgxT6a1vcmyNeD93/fx7qDyYXHvr2jD3GnM/Fwc2FPQjrTlu8t8xqBPu6kZOaVKDs0dUSJ/fwCC26ydmWNkKH0QtSW8y0bZ+NRjxey6chp3li8u8xj397Rh56RQfSMNJ6aLRZNXoGFj1buL3HeV7f3ok9UA7SGzLwCtsWl0K2MIe+SvGufJHAhHNj3G+MKt+8b2JKbejUn0Med+JQsokJKLjDs4qJ4bFhbHhvWtsxrKQV+nm70bRFcozGLipMELoSDOp6axZx1Rwr37xnQEi93Y2Tkuclb2CdJ4EI4qB82HwPghdEd+Ffv5jJZlAOSBC6Eg/niz4M8/3MsAC1CfLmuRzNJ3g5KErgQDiInv4D//XGAt5fuASDYz4PFk/vJy0UHJglcCAeQlp3HFe+uJu50FgAPDWnNvQNbypO3g5MELoSd01oz/vP1hcl789NDqO9r+4mrRN0jCVwIO/X3gVO8+utOtsYVjap8+9oukrydiCRwIezQit2JTPhiQ+H+61d3ZnTXxrKAgpORBC6EHTmanMknqw8wa+1hAJY+2I/wBj6SuJ2UJHAh7ETMoWSu/nht4f7//tWdVo38TYxImE36FwlhJ3bEpxVu3z+oFUM7hJ7nbOEMJIELYQc2Hj7NV8WGxU9bvpfUc2YHFM5HmlCEqMMenLuFXSfOsPN4WonyJoHeeHnI85ezkwQuRB2TmZvPrV9sYH2xebsBHhlqrFjo6qK4o1+UDNIRF07gSqlmwCwgFLAA07XW7yqlgoC5QARwCLhWa3265kIVwjmk5+QXJu+BbRtyXY9mBHq70yuqgcmRibqmIk/g+cBDWutNSil/YKNSailwK7C82JqYU4DHai5UIRzfbTNj+OdYSuF+x8b15GWlKNcFG9G01se11pus22eAnUATYDQw03raTGBMDcUohFPYl5jOsp0JJKTlFJZFhtTuCj7CvlSqDVwpFQFcBKwDGmmtj4OR5JVSDcv5ziRgEkB4eHhZpwghMOY0AXh1bCcGt2tEiL+nyRGJuq7Cr7GVUn7A98BkrXXahc4/S2s9XWsdrbWODgkJqUqMQjiFbdY5Tfw83SR5iwqpUAJXSrljJO85Wuv51uIEpVSY9XgYkFgzIQrh+A6dzOC37ScAmLX2kLnBCLtRkV4oCvgM2Km1frvYoQXAeGCq9fOnGolQCAegtWblniT+PnCKX/85TpNAb3pGNuDu/i344s9DvLZoV+G5V3dvamKkwp6os+1u5Z6g1CXAauAfjG6EAE9gtIN/C4QDR4BrtNbJZV7EKjo6WsfExFQ3ZiHqJK01E2ZsYOXuJPa8NBwPNxd2Hk9j+LurK/T9R4a24e7+LaR/tyhFKbVRax19bvkFn8C11muA8n6iBlU3MCEcwY74VEZMW1O4P/jtP8jJLyjRo2RQ24bEnc6ivq87n9/ag/bPLC48tuHJwdLuLSpNRmIKUU05+QUlkjfAkeTMwu2XxnTkpl7hpZ6sv7uzD9+sP8p9A1tK8hZVIglciGp4/ucdfPHnoRJlX93eiz5RDdiflEF2XgHtw+qV2SzSIyKIHhFBtRSpcESSwIWoggKLpsUTv5YoW/XIABrW88TL3VhcoWVDPzNCE05EErgQlXDoZAY/bD7G/M1xhWVdmgUy49YeshalqHWSwIWooF0n0hj2X6NHSbCfJ2O7NaF9WD0ubx8qyVuYQhK4EBWQV2DhyveKXlS+d8NF9GkhswMKc0kCF+ICLBbNrV+sJ69A0zbUn0WT+5kdkhCAJHAhzmt/UjpvLt7Nn/tOUc/LjbsHtDQ7JCEKSQIXohwbD5/mls/WkZFbwOB2jfjklu4ySlLUKZLAhTjHByv2MX9THAlpOYT4e/L5uM50aRYoyVvUOZLAhSgmO6+ANxbvBox+3LP+3ZPGgd4mRyVE2SSBC2EVn5JF36m/A3B9j2ZMHdfZ5IiEOL8KL+gghD3LyS/g/d/3si8xvdxzYuOL1il5dWyn2ghLiGqRBC4cXn6Bhbtnb+LNJXsY/u4qjpzKLHW8wKLZdsxYEadFiK+0dwu7IE0owuEtjU1g+a5EhncM5bftJxj53mqujW7Gw0PbkJGTT/eXlpU4f0TnxiZFKkTlSAIXDuuXbfHM3XCUvw+cItjPk6njOtOhcT3eXLKHT9ccZHHsCY4mZ5X4ziND23DnZS1MiliIypEELhxC3OlMxn74F12bBdIzMojeUQ2496vNhcf/M6Q1Ad7u3DuwFZd3CGXl7kR+236CdqH1aFTPiwkXRxAVIrMHCvtSkSXVPgdGAola647WsiBgLhABHAKu1VqfvtDNZEk1URPSsvOIfmkZufmWUsceGdqG63o0I9hPFkwQ9qu8JdUq8hJzBjDsnLIpwHKtdStguXVfiFo3f1McvV9ZTm6+BU83Fw68cgV+nsYvlnNu68WkflGSvIXDqsiamKuUUhHnFI8G+lu3ZwIrgcdsGZgQZckvsBB3OouDpzLYdjSVd5btISrYl4mXRtKpSQAuLortzw81O0whakVV28Abaa2PA2itjyulGpZ3olJqEjAJIDw8vIq3EwJOpecwee4WVu89WaL83esvolPTAJOiEsI8Nf4SU2s9HZgORht4Td9P2D+tNTn5FrzcXSmwaCbM2MCqPUmFxz1cXejSLIANh4zXLpK8hbOqagJPUEqFWZ++w4BEWwYlnNey2ARum2W86Pb3dMPT3YWT6bmFx2/sFc4rV8koSSGg6gl8ATAemGr9/MlmEQmn9ds/x7lrzqbC/XHdm5KSmcuQ9qEM7xhKalaeLF0mRDEXTOBKqa8xXlgGK6XigGcxEve3SqmJwBHgmpoMUjiG/UnprNydxIu/xNI21J9Px0ezJ+EMm4+k8N7v+wrPmzGhB5e1Dik1nF2StxAlVaQXyg3lHBpk41iEg0rPyWfa8r1MX3WgsGzXiTNc8tqKUucGeLvTv02578SFEMXISExR415euJOv1x8B4L6BLZl4SSR/7jvFlqOnGdohlE9WH2DxjgT+96/u9GsVYnK0QtgPSeDCpnLzLXi4uZCTX8Da/adYdzC5MHkvmnwpbUPrATCicxgjOocBEB0RZFq8QtgzSeDCJnLzLXy/KY7H5/9D12aBbDmaAoBSMKhtQ96/sRveHq7mBimEg5EELsqVnJHLByv2cfulUYQGeGGxaNYeOMXHf+znuVEdWLErkQBvdxZsjS8xuGbL0RTah9WjS7MAnriiHf5e7ibWQgjHJQlclGlbXAq3frGB5IxcZq09RIifJyfSsrFYh2INeuuPUt/pGRlE+7B63NgrnNaN/Gs5YiGcjyRwAUBegYW1+0/h7eHKN+uPsjT2BFrDqC6NOZmeQ31fD5oH+ZB0JoelOxNIycwr/O7Atg2ZOq4TDf29TKyBEM5HErjghZ9j+fzPgyXKBrdryAODWsswdSHqMEngTmrn8TSWxiZwJjuvRPIe3K4RT49sR/MGviZGJ4SoCEngTuDgyQz8vdyIOZTMgq3xrNpzkvScfAC83V0J8vXgiSvacXX3piZHKoSoDEngDujsKktKKaav2s8rv+4qcdzNRdGmkT9Xd2/K7f2izAhRCGEDksAdTGpmHl1eWAJAl6YBbI1L5bLWIbQJ9ae+jwe39GmOr6f8bxfCEcjfZAfz7ILtJfanDG/LrX0j8HKXQTRCOBpJ4A4iIyef+JQsftwSD8ClrYL5cmIvk6MSQtQkSeB11N6EM+xOOEPPyCAsFvD2cGXVniR2nUijU5MA+rdpiKebC0oplu9M4D/fbiU1y+ibPaJTGB/c1M3kGgghapok8FqSV2BhR3waPh6u5Y5STDyTzYNztxDo48HCbccrdf22of48P6oDJ9NzGNm5sS1CFkLUcZLAa5jWGqUUrZ78rbCsb4sGREcE0TsqiA6NA4iNT2NJ7Am+i4kjPScfX+ukTx2b1GNcN6Nr3+nMPHpGBOHv5cbyXYlMW74XgAFtQmjg58mLozvKZFFCOBl1tstZlb6s1DDgXcAV+FRrPfV850dHR+uYmJgq36+uS0zLJj0nn4Fv/cHQDo1IzsgtXHj3rEtbBbM3IZ0Tadmlvj+4XUPu6t+Cjk0C0JoLvnjMziuQl5NCOAGl1EatdfS55VV+AldKuQIfAEOAOGCDUmqB1jq26mHarz0JZxj1/hqy8ywALN6RUOL4ub1BjqVk8eGKfWyNS+HGns3p3yaExoHelbqnJG8hnFt1mlB6Avu01gcAlFLfAKMBh0/giWnZbDpyGouG1Kw8Hp//D2DMfX1Jy2Aign24/dIoEtJyiArxJSu3gGZBPiWu0STQm5dldXUhRDVUJ4E3AY4W248DSvVbU0pNAiYBhIeHV+N25oo7ncnbS/fw9/5TxKeWbv4AmDGhJ5e1LloSTOYTEULUpOokcFVGWakGda31dGA6GG3g1bhfrbNYNP9bdYDVe5PYcCiZfItmaPtQJl4aRHxKFt2b16ddWD0aB3qRk2+hnixcIISoRdVJ4HFAs2L7TYH46oVjjt0nznAmOw9vD1fyCzRr9p1kyY4T3HlZC15bZMwjctslkfyrT/Nyn6o93aQ9WghRu6qTwDcArZRSkcAx4HrgRptEVYPyCyzEp2Tz0sJYLgqvT30fd578cTsFltK/HNw1ZxMAjw1ry139W9R2qEIIcV5VTuBa63yl1L3AYoxuhJ9rrXfYLLJi9iSc4URqNv2KtS8Xl5iWzYcr95OWlUfTIB82Hk6mga8nmbn5NA70pk2oPydSszmWksX8TccKv7ck1ugpEhbgxeTBrfD1dMPb3ZWkMzm4ubqQlVeAr4crwzuG1US1hBCiWqo1kEdr/Svwq41iKddHK/fzw+ZjjOrSmKdGtCPQx4N5G+OYtnwvGbn5nMnOL3F+h8b1OJCUga+nG38fSCY9Jx9XF4WXmwsAQ9o34pruTenUNIA9Cen0jAiSQTBCCLtjFyMxXx3bieYNfPhwxX4WbC1qZu/SLJAuTQNo3sCXi8IDaRroTUpWXomh6haL5lhKFg38PPDxKF3dsIDK9b0WQoi6wi4SuJe7K5MHt2ZUl8Z8+fdhUjPz6BEZxHXRzXBxKdkZpmG9kgvrurioUn2whRDCEdhFAj8rKsSPZ6/sYHYYQghRJ7iYHYAQQoiqkQQuhBB2ShK4EELYKUngQghhpySBCyGEnZIELoQQdkoSuBBC2ClJ4EIIYaeqtSZmpW+mVBJw2MaXDQZO2viadY3U0TFIHR2DGXVsrrUuNZtfrSbwmqCUiilrsU9HInV0DFJHx1CX6ihNKEIIYackgQshhJ1yhAQ+3ewAaoHU0TFIHR1Dnamj3beBCyGEs3KEJ3AhhHBKksCFEMJO2W0CV0oNU0rtVkrtU0pNMTseW1BKNVNKrVBK7VRK7VBKPWAtD1JKLVVK7bV+1jc71upSSrkqpTYrpX6x7jtUHZVSgUqpeUqpXdb/n30csI4PWn9OtyulvlZKeTlCHZVSnyulEpVS24uVlVsvpdTj1jy0Wyk1tDZjtcsErpRyBT4AhgPtgRuUUu3Njcom8oGHtNbtgN7APdZ6TQGWa61bAcut+/buAWBnsX1Hq+O7wCKtdVugC0ZdHaaOSqkmwP1AtNa6I+AKXI9j1HEGMOycsjLrZf37eT3QwfqdD635qVbYZQIHegL7tNYHtNa5wDfAaJNjqjat9XGt9Sbr9hmMv/RNMOo203raTGCMKQHaiFKqKTAC+LRYscPUUSlVD+gHfAagtc7VWqfgQHW0cgO8lVJugA8QjwPUUWu9Ckg+p7i8eo0GvtFa52itDwL7MPJTrbDXBN4EOFpsP85a5jCUUhHARcA6oJHW+jgYSR5oaGJotvBf4FHAUqzMkeoYBSQBX1ibiT5VSvniQHXUWh8D3gSOAMeBVK31Ehyojucor16m5iJ7TeCqjDKH6Q+plPIDvgcma63TzI7HlpRSI4FErfVGs2OpQW5AN+AjrfVFQAb22ZRQLmsb8GggEmgM+CqlbjY3KlOYmovsNYHHAc2K7TfF+PXN7iml3DGS9xyt9XxrcYJSKsx6PAxINCs+G7gYGKWUOoTR9DVQKTUbx6pjHBCntV5n3Z+HkdAdqY6DgYNa6yStdR4wH+iLY9WxuPLqZWoustcEvgFopZSKVEp5YLxEWGByTNWmlFIY7aY7tdZvFzu0ABhv3R4P/FTbsdmK1vpxrXVTrXUExv+337XWN+NYdTwBHFVKtbEWDQJicaA6YjSd9FZK+Vh/bgdhvLNxpDoWV169FgDXK6U8lVKRQCtgfa1FpbW2yz/AFcAeYD/wpNnx2KhOl2D8+rUN2GL9cwXQAOPN917rZ5DZsdqovv2BX6zbDlVHoCsQY/1/+SNQ3wHr+DywC9gOfAl4OkIdga8x2vXzMJ6wJ56vXsCT1jy0Gxhem7HKUHohhLBT9tqEIoQQTk8SuBBC2ClJ4EIIYackgQshhJ2SBC6EEHZKErgQQtgpSeBCCGGn/g9OURSfyGybYQAAAABJRU5ErkJggg==\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[1::2]\n",
	"tlist2 = tlist[::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",
	"\n",
	"\n",
	"# how to wrangle this data into StochasticTrajectoryData ?"
	]
	},
	{
	"cell_type": "code",
	"execution_count": null,
	"id": "b2622bd8-2801-4a29-b601-102c5b496edb",
	"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
	}