{ "cells": [ { "cell_type": "code", "execution_count": 455, "metadata": {}, "outputs": [], "source": [ "import pandas as pd\n", "from scipy.stats import poisson\n", "import numpy as np\n", "import math\n", "import seaborn as sns\n", "import matplotlib as mpl\n", "import matplotlib.pyplot as plt" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# UTILS" ] }, { "cell_type": "code", "execution_count": 456, "metadata": {}, "outputs": [], "source": [ "def bound_norm_random(low, high):\n", " res = np.random.normal((high+low)/2,(high-low)/6)\n", " if (reshigh):\n", " res = bound_norm_random(low, high)\n", " return res" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# TYPE/PHASE" ] }, { "cell_type": "code", "execution_count": 457, "metadata": {}, "outputs": [], "source": [ "EXPERIMENT_TYPES = ['1 off run', 'Monte Carlo', 'Monte Carlo Parameter Sweep', 'Monte Carlo Pairwise']\n", "\n", "experiment_type = EXPERIMENT_TYPES[1]\n", "monte_carlo_runs = 100\n", "\n", "#correct number of runs if inconsistent with experiment type\n", "if (experiment_type == EXPERIMENT_TYPES[0]):\n", " monte_carlo_runs = 1" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# TIMESCALE" ] }, { "cell_type": "code", "execution_count": 458, "metadata": {}, "outputs": [], "source": [ "SECOND = 1\n", "MINUTE = 60*SECOND\n", "HOUR = 60*MINUTE\n", "DAY = 24*HOUR\n", "DURATION_OF_UNIT_OF_TIME = 1*DAY\n", "\n", "experiment_steps = 1000\n", "experiment_duration = experiment_steps * DURATION_OF_UNIT_OF_TIME\n", "time_array = np.arange(0,experiment_steps*DURATION_OF_UNIT_OF_TIME,DURATION_OF_UNIT_OF_TIME)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# MECHANISMS (dimension)" ] }, { "cell_type": "code", "execution_count": 459, "metadata": {}, "outputs": [], "source": [ "mechanisms_names = ['mech_one', 'mech_two', 'mech_three']" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# STATES (dimension)" ] }, { "cell_type": "code", "execution_count": 460, "metadata": {}, "outputs": [], "source": [ "states_names = ['a', 'b', 'c']\n", "states_data = [[np.zeros(experiment_steps)]*len(states_names)]*monte_carlo_runs\n", "states_data = np.zeros((monte_carlo_runs, experiment_steps, len(states_names)), dtype=int)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Initial Conditions" ] }, { "cell_type": "code", "execution_count": 461, "metadata": {}, "outputs": [], "source": [ "states_0 = {\n", " 'a': 0,\n", " 'b': 0,\n", " 'c': 300\n", "}\n", "# an initial condition must be set for every state\n", "assert np.array([k in states_0 for k in states_names]).all(), 'Error: The initial condition of one or more states is unkonwn'\n", "\n", "# copy initial condition to the states dataset\n", "for i in range(len(states_names)):\n", " states_data[:,0,i] = states_0[states_names[i]]" ] }, { "cell_type": "code", "execution_count": 462, "metadata": {}, "outputs": [], "source": [ "T_0 = 0\n", "time_array = T_0 + time_array" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Mechanisms Coef (params)" ] }, { "cell_type": "code", "execution_count": 463, "metadata": {}, "outputs": [], "source": [ "mech_one_coef_A = 0.05" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# MECHANISMS EQUATIONS (func)" ] }, { "cell_type": "code", "execution_count": 464, "metadata": {}, "outputs": [], "source": [ "def mech_one(_states_data, _time_array, _run, _step, args):\n", "# print('mech 1')\n", " _states_data[_run, _step, states_names.index('a')] += (1-mech_one_coef_A)*args[0]\n", " _states_data[_run, _step, states_names.index('b')] += mech_one_coef_A*args[0]\n", " return _states_data\n", "\n", "def mech_two(_states_data, _time_array, _run, _step, args):\n", "# print('mech 2')\n", " _states_data[_run, _step, states_names.index('a')] -= args[0]\n", " return _states_data\n", "\n", "def mech_three(_states_data, _time_array, _run, _step, args):\n", "# print('mech 3')\n", " _states_data[_run, _step, states_names.index('b')] -= args[0]\n", " return _states_data\n", "\n", "def mech_four(_states_data, _time_array, _run, _step):\n", "# print('mech 4')\n", " _states_data[_run, _step, states_names.index('a')] = _states_data[_run, _step-1, states_names.index('a')]\n", " _states_data[_run, _step, states_names.index('b')] = _states_data[_run, _step-1, states_names.index('b')] \n", " return _states_data\n", "\n", "mechanisms = [eval(m) for m in mechanisms_names]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Behavioral Model Coef (params) " ] }, { "cell_type": "code", "execution_count": 465, "metadata": {}, "outputs": [], "source": [ "behavior_one_coef_A = 0.01\n", "behavior_one_coef_B = -0.01\n", "\n", "behavior_two_coef_A = 1" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# BEHAVIORAL MODEL (func)" ] }, { "cell_type": "code", "execution_count": 466, "metadata": {}, "outputs": [], "source": [ "behaviors_names = ['behavior_one', 'behavior_two']\n", "def behavior_one(_states_data, _time_array, _run, _step): \n", " c_var = ( _states_data[_run, _step, states_names.index('c')]\n", " - _states_data[_run, _step-1, states_names.index('c')] )\n", " c_var_perc = c_var / _states_data[_run, _step-1, states_names.index('c')]\n", " \n", " if (c_var_perc > behavior_one_coef_A):\n", " return mech_one(_states_data, _time_array, _run, _step, [c_var])\n", " elif (c_var_perc < behavior_one_coef_B):\n", " return mech_two(_states_data, _time_array, _run, _step, [-c_var])\n", " return _states_data\n", "\n", "def behavior_two(_states_data, _time_array, _run, _step):\n", " b_balance = _states_data[_run, _step-1, states_names.index('b')]\n", " if (b_balance > behavior_two_coef_A):\n", " return mech_three(_states_data, _time_array, _run, _step, [b_balance])\n", " return _states_data\n", "\n", "behaviors = [eval(b) for b in behaviors_names] " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# ENVIRONMENTAL PROCESS (dimension)" ] }, { "cell_type": "code", "execution_count": 467, "metadata": {}, "outputs": [], "source": [ "env_proc_names = ['proc_one']" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Stochastic Process Coef (params)" ] }, { "cell_type": "code", "execution_count": 468, "metadata": {}, "outputs": [], "source": [ "proc_one_coef_A = 0.7\n", "proc_one_coef_B = 1.3" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# ENVIRONMENTAL PROCESS (func)" ] }, { "cell_type": "code", "execution_count": 469, "metadata": {}, "outputs": [], "source": [ "def proc_one(_states_data, _time_array, _run, _step):\n", " _states_data[_run, _step, states_names.index('a')] = _states_data[_run, _step-1, states_names.index('a')]\n", " _states_data[_run, _step, states_names.index('b')] = _states_data[_run, _step-1, states_names.index('b')] \n", " _states_data[_run, _step, states_names.index('c')] = ( _states_data[_run, _step-1, states_names.index('c')]\n", " * bound_norm_random(proc_one_coef_A, proc_one_coef_B) )\n", " return _states_data\n", "\n", "env_proc = [eval(p) for p in env_proc_names]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# ENGINE" ] }, { "cell_type": "code", "execution_count": 470, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/Users/markusbkoch/.local/share/virtualenvs/DiffyQ-SimCAD-4_qpgnP9/lib/python3.6/site-packages/ipykernel_launcher.py:5: RuntimeWarning: invalid value encountered in long_scalars\n", " \"\"\"\n" ] } ], "source": [ "for i in range(monte_carlo_runs):\n", " for t in range(1,experiment_steps):\n", " for p in env_proc:\n", " states_data = p(_states_data=states_data,\n", " _time_array=time_array, \n", " _run=i, \n", " _step=t)\n", " for b in behaviors:\n", " states_data = b(_states_data=states_data,\n", " _time_array=time_array, \n", " _run=i, \n", " _step=t)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# DATA COLLECTION" ] }, { "cell_type": "code", "execution_count": 471, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
abc
0000300
86400150316
172800-290272
259200283332
345600190323
432000260331
518400561363
604800201327
691200402349
777600220331
864000-230286
950400-520257
1036800-330276
1123200-330274
1209600-210286
1296000-150292
1382400-200287
1468800-10306
155520090317
164160020310
172800060315
1814400-660243
1900800-1060203
1987200-1190190
2073600-1170192
2160000-1450164
2246400-1520157
2332800-1640145
2419200-1660143
2505600-1680141
...............
9983808000-30200
83894400-30200
83980800-30200
84067200-30200
84153600-30200
84240000-30200
84326400-30200
84412800-30200
84499200-30200
84585600-30200
84672000-30200
84758400-30200
84844800-30200
84931200-30200
85017600-30200
85104000-30200
85190400-30200
85276800-30200
85363200-30200
85449600-30200
85536000-30200
85622400-30200
85708800-30200
85795200-30200
85881600-30200
85968000-30200
86054400-30200
86140800-30200
86227200-30200
86313600-30200
\n", "

100000 rows × 3 columns

\n", "
" ], "text/plain": [ " a b c\n", "0 0 0 0 300\n", " 86400 15 0 316\n", " 172800 -29 0 272\n", " 259200 28 3 332\n", " 345600 19 0 323\n", " 432000 26 0 331\n", " 518400 56 1 363\n", " 604800 20 1 327\n", " 691200 40 2 349\n", " 777600 22 0 331\n", " 864000 -23 0 286\n", " 950400 -52 0 257\n", " 1036800 -33 0 276\n", " 1123200 -33 0 274\n", " 1209600 -21 0 286\n", " 1296000 -15 0 292\n", " 1382400 -20 0 287\n", " 1468800 -1 0 306\n", " 1555200 9 0 317\n", " 1641600 2 0 310\n", " 1728000 6 0 315\n", " 1814400 -66 0 243\n", " 1900800 -106 0 203\n", " 1987200 -119 0 190\n", " 2073600 -117 0 192\n", " 2160000 -145 0 164\n", " 2246400 -152 0 157\n", " 2332800 -164 0 145\n", " 2419200 -166 0 143\n", " 2505600 -168 0 141\n", "... ... .. ...\n", "99 83808000 -302 0 0\n", " 83894400 -302 0 0\n", " 83980800 -302 0 0\n", " 84067200 -302 0 0\n", " 84153600 -302 0 0\n", " 84240000 -302 0 0\n", " 84326400 -302 0 0\n", " 84412800 -302 0 0\n", " 84499200 -302 0 0\n", " 84585600 -302 0 0\n", " 84672000 -302 0 0\n", " 84758400 -302 0 0\n", " 84844800 -302 0 0\n", " 84931200 -302 0 0\n", " 85017600 -302 0 0\n", " 85104000 -302 0 0\n", " 85190400 -302 0 0\n", " 85276800 -302 0 0\n", " 85363200 -302 0 0\n", " 85449600 -302 0 0\n", " 85536000 -302 0 0\n", " 85622400 -302 0 0\n", " 85708800 -302 0 0\n", " 85795200 -302 0 0\n", " 85881600 -302 0 0\n", " 85968000 -302 0 0\n", " 86054400 -302 0 0\n", " 86140800 -302 0 0\n", " 86227200 -302 0 0\n", " 86313600 -302 0 0\n", "\n", "[100000 rows x 3 columns]" ] }, "execution_count": 471, "metadata": {}, "output_type": "execute_result" } ], "source": [ "data = pd.DataFrame(states_data[0], \n", " index=[[0]*experiment_steps, time_array], \n", " columns=states_names)\n", "for i in range(1,monte_carlo_runs):\n", " b = pd.DataFrame(states_data[i],\n", " index=[[i]*experiment_steps, time_array], \n", " columns=states_names)\n", " data = data.append(b)\n", "data" ] } ], "metadata": { "kernelspec": { "display_name": "DiffyQ-SimCAD Env", "language": "python", "name": "diffyq-simcad" }, "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.5" } }, "nbformat": 4, "nbformat_minor": 2 }