diff --git a/Sim_CAD_Engine.ipynb b/Sim_CAD_Engine.ipynb new file mode 100644 index 0000000..93e7aad --- /dev/null +++ b/Sim_CAD_Engine.ipynb @@ -0,0 +1,874 @@ +{ + "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 +}