{ "cells": [ { "cell_type": "code", "execution_count": 491, "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": 492, "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": 493, "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": 494, "metadata": {}, "outputs": [], "source": [ "SECOND = 1\n", "MINUTE = 60*SECOND\n", "HOUR = 60*MINUTE\n", "DAY = 24*HOUR\n", "DURATION_OF_A_STEP = 1*DAY\n", "\n", "experiment_steps = 1000\n", "experiment_duration = experiment_steps * DURATION_OF_A_STEP\n", "time_array = np.arange(0,experiment_steps*DURATION_OF_A_STEP,DURATION_OF_A_STEP)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# MECHANISMS (dimension)" ] }, { "cell_type": "code", "execution_count": 495, "metadata": {}, "outputs": [], "source": [ "mechanisms_names = ['mech_one', 'mech_two', 'mech_three']" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# STATES (dimension)" ] }, { "cell_type": "code", "execution_count": 496, "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)\n", "# states_data is a 3-dimensional array - montecarlo, time, states\n", "# montecarlo[time[states]]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Initial Conditions" ] }, { "cell_type": "code", "execution_count": 497, "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": 498, "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": 499, "metadata": {}, "outputs": [], "source": [ "mech_one_coef_A = 0.05" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# MECHANISMS EQUATIONS (func)" ] }, { "cell_type": "code", "execution_count": 500, "metadata": {}, "outputs": [], "source": [ "# state/mechanism matrix\n", "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": 501, "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": 502, "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": 503, "metadata": {}, "outputs": [], "source": [ "env_proc_names = ['proc_one']" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Stochastic Process Coef (params)" ] }, { "cell_type": "code", "execution_count": 504, "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": 505, "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": 506, "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) #behaviors have access to exogenous data @ step-1, not @ step" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# DATA COLLECTION" ] }, { "cell_type": "code", "execution_count": 507, "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
86400-70293
172800-200280
259200-860214
345600-671234
432000-671233
518400-911209
604800-991201
691200-931207
777600-1041196
864000-1011199
950400-861214
1036800-921208
1123200-921210
1209600-991203
1296000-1221180
1382400-1351167
1468800-1611141
1555200-1611141
1641600-1871115
1728000-1951107
1814400-2011101
1900800-1891113
1987200-1891112
2073600-1891111
2160000-1911109
2246400-2001100
2332800-206194
2419200-206194
2505600-218182
...............
9983808000-38600
83894400-38600
83980800-38600
84067200-38600
84153600-38600
84240000-38600
84326400-38600
84412800-38600
84499200-38600
84585600-38600
84672000-38600
84758400-38600
84844800-38600
84931200-38600
85017600-38600
85104000-38600
85190400-38600
85276800-38600
85363200-38600
85449600-38600
85536000-38600
85622400-38600
85708800-38600
85795200-38600
85881600-38600
85968000-38600
86054400-38600
86140800-38600
86227200-38600
86313600-38600
\n", "

100000 rows × 3 columns

\n", "
" ], "text/plain": [ " a b c\n", "0 0 0 0 300\n", " 86400 -7 0 293\n", " 172800 -20 0 280\n", " 259200 -86 0 214\n", " 345600 -67 1 234\n", " 432000 -67 1 233\n", " 518400 -91 1 209\n", " 604800 -99 1 201\n", " 691200 -93 1 207\n", " 777600 -104 1 196\n", " 864000 -101 1 199\n", " 950400 -86 1 214\n", " 1036800 -92 1 208\n", " 1123200 -92 1 210\n", " 1209600 -99 1 203\n", " 1296000 -122 1 180\n", " 1382400 -135 1 167\n", " 1468800 -161 1 141\n", " 1555200 -161 1 141\n", " 1641600 -187 1 115\n", " 1728000 -195 1 107\n", " 1814400 -201 1 101\n", " 1900800 -189 1 113\n", " 1987200 -189 1 112\n", " 2073600 -189 1 111\n", " 2160000 -191 1 109\n", " 2246400 -200 1 100\n", " 2332800 -206 1 94\n", " 2419200 -206 1 94\n", " 2505600 -218 1 82\n", "... ... .. ...\n", "99 83808000 -386 0 0\n", " 83894400 -386 0 0\n", " 83980800 -386 0 0\n", " 84067200 -386 0 0\n", " 84153600 -386 0 0\n", " 84240000 -386 0 0\n", " 84326400 -386 0 0\n", " 84412800 -386 0 0\n", " 84499200 -386 0 0\n", " 84585600 -386 0 0\n", " 84672000 -386 0 0\n", " 84758400 -386 0 0\n", " 84844800 -386 0 0\n", " 84931200 -386 0 0\n", " 85017600 -386 0 0\n", " 85104000 -386 0 0\n", " 85190400 -386 0 0\n", " 85276800 -386 0 0\n", " 85363200 -386 0 0\n", " 85449600 -386 0 0\n", " 85536000 -386 0 0\n", " 85622400 -386 0 0\n", " 85708800 -386 0 0\n", " 85795200 -386 0 0\n", " 85881600 -386 0 0\n", " 85968000 -386 0 0\n", " 86054400 -386 0 0\n", " 86140800 -386 0 0\n", " 86227200 -386 0 0\n", " 86313600 -386 0 0\n", "\n", "[100000 rows x 3 columns]" ] }, "execution_count": 507, "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 }