Chia-data-analysis.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import csv\n",
    "import datetime\n",
    "import math\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.patches as mpatches\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import subprocess                               # ssh command\n",
    "import time\n",
    "\n",
    "from colour import Color\n",
    "from matplotlib.lines import Line2D\n",
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "from scipy.optimize import curve_fit\n",
    "from scipy.stats import binom\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Check notebook server type"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# % indicates magic funciton -> runs command. Run %magic in a cell for more info\n",
    "# 'notebook' is used for jupyter notebook\n",
    "# 'widget'   is used for jupyter lab\n",
    "\n",
    "def is_lab_notebook():\n",
    "        import re\n",
    "        import psutil\n",
    "        \n",
    "        return any(re.search('jupyter-lab', x)\n",
    "                   for x in psutil.Process().parent().cmdline())\n",
    "    \n",
    "\n",
    "if is_lab_notebook():\n",
    "    %matplotlib widget\n",
    "else:\n",
    "    %matplotlib notebook \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Sync latest file from server"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# TODO: move these definitions to a config file. This would allow\n",
    "#       more portability between types of jupyter nb instances (notebook/lab/hub)\n",
    "\n",
    "SSH_PATH='ethanblagg@10.0.0.9:~/Documents/scripts/'\n",
    "SSH_IDENTITY_FILE='~/.ssh/id_rsa_ubuntu_server'\n",
    "LOCAL_PATH='/home/ethanblagg/jupyter-lab-files/chia-analytics/data/'\n",
    "F_LOG='netspace-size.log'\n",
    "\n",
    "\n",
    "try:\n",
    "    print(subprocess.check_output(['scp', '-i', SSH_IDENTITY_FILE,\n",
    "                                   SSH_PATH + F_LOG,\n",
    "                                   LOCAL_PATH\n",
    "                                  ]))\n",
    "        \n",
    "except subprocess.CalledProcessError as err:\n",
    "    print(err)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Get data from CSV"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "file_name = LOCAL_PATH + F_LOG\n",
    "\n",
    "epoch_time_arr = []\n",
    "chia_farmed_arr = []\n",
    "num_plots_arr = []\n",
    "plots_size_PiB_arr = []\n",
    "net_size_PiB_arr = []\n",
    "win_prob_arr = []\n",
    "headers = []\n",
    "\n",
    "with open(file_name, newline='') as csvfile:\n",
    "    data = csv.reader(csvfile, delimiter=',', quotechar='|')\n",
    "    \n",
    "    headers = next(data)\n",
    "    \n",
    "    for i, row in enumerate(data):\n",
    "        #print(i)\n",
    "        epoch_time_arr.append(float(row[0]))\n",
    "        chia_farmed_arr.append(float(row[1]))\n",
    "        num_plots_arr.append(float(row[2]))\n",
    "        plots_size_PiB_arr.append(float(row[3]))\n",
    "        net_size_PiB_arr.append(float(row[4]))\n",
    "        win_prob_arr.append(float(row[5]))\n",
    "\n",
    "        \n",
    "epoch_times    = np.array(epoch_time_arr)\n",
    "chia_farmed    = np.array(chia_farmed_arr)\n",
    "num_plots      = np.array(num_plots_arr)\n",
    "plots_size_PiB = np.array(plots_size_PiB_arr)\n",
    "net_size_PiB   = np.array(net_size_PiB_arr)\n",
    "win_prob       = np.array(win_prob_arr)\n",
    "\n",
    "\n",
    "# create date-time (number) array for plotting\n",
    "tts = np.copy(epoch_times)                         # seconds since epoch\n",
    "ttd = np.copy(epoch_times)                         # days since epoch\n",
    "\n",
    "\n",
    "\n",
    "tts_zeroed = np.copy(tts) - tts[0]\n",
    "\n",
    "for i, t in enumerate(ttd):\n",
    "    ttd[i] = matplotlib.dates.epoch2num(ttd[i])\n",
    "\n",
    "# Convert times to local time (Denver, -6:00 UTC)\n",
    "print(ttd[0])\n",
    "UTC_CORRECTION = -6 / 24\n",
    "ttd = ttd + UTC_CORRECTION\n",
    "print(ttd[0])\n",
    "    \n",
    "ttd_zeroed = np.copy(ttd) - ttd[0]\n",
    "\n",
    "first_sample_index = 219\n",
    "\n",
    "# Deprecated....\n",
    "xx = np.copy(epoch_times)\n",
    "\n",
    "\n",
    "print(\"epoch_time_arr[-1]: {}\".format(epoch_time_arr[-1]))\n",
    "print(\"tts[0]: {} \\t\\t tts[-1]: {}\".format(tts[0], tts[-1]))\n",
    "print(\"ttd[0]: {} \\t ttd[-1]: {}\".format(ttd[0], ttd[-1]))\n",
    "print(\"ttdz[0]: {} \\t\\t\\t ttdz[-1]: {}\".format(ttd_zeroed[0], ttd_zeroed[-1]))\n",
    "print(\"ttsz[0]: {} \\t\\t\\t ttsz[-1]: {}\".format(tts_zeroed[0], tts_zeroed[-1]))\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Netspace size prediction\n",
    "##### Data for when updates are released"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "updates = [\n",
    "    [np.datetime64('2021-04-14T13:00'), '1.0.5'],\n",
    "    [np.datetime64('2021-04-21T13:00'), '1.1.0'],\n",
    "    [np.datetime64('2021-04-22T13:00'), '1.1.1'],\n",
    "    [np.datetime64('2021-04-25T13:00'), '1.1.2'],\n",
    "    [np.datetime64('2021-05-01T13:00'), '1.1.3'],\n",
    "    [np.datetime64('2021-05-04T12:00'), '1.1.4'],\n",
    "    [np.datetime64('2021-05-09T13:00'), '1.1.5'],\n",
    "    [np.datetime64('2021-05-20T18:30'), '1.1.6'],\n",
    "    [np.datetime64('2021-06-05T18:30'), '1.1.7']\n",
    "          ]\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "PLOT_EXP = False\n",
    "PLOT_SIG = True\n",
    "PLOT_UPDATES = True\n",
    "\n",
    "prediction_interval = 50\n",
    "\n",
    "sigmoid_offset = int(net_size_PiB[0]) # correction for y offset for sigmoid (doesn't want to work in func_sigmoid)\n",
    "\n",
    "# Curve fitting functions\n",
    "def func_exp(x, a, b, c):\n",
    "    return (a * (x**b)) + c\n",
    "\n",
    "def func_sigmoid(x, a, b, c, d):\n",
    "    return (d / (a + b*math.e**(-c*x))) + sigmoid_offset\n",
    "\n",
    "\n",
    "# Create arrays for time (x-axis) that will be used for predicting the future\n",
    "# We will predict with seconds, but plot with days (matplot lib uses days since epoch - ?)\n",
    "last_day = math.floor(tts_zeroed[-1])\n",
    "tt_future_tail = np.arange(last_day, last_day+365*24*3600, 24*3600)\n",
    "tt_future_s = np.concatenate((tts_zeroed, tt_future_tail)).astype(int)\n",
    "tt_future_d_z = (tt_future_s / 3600 / 24)#.astype(int)\n",
    "tt_future_d = tt_future_d_z + ttd[0]\n",
    "\n",
    "\n",
    "# Fit curves\n",
    "popt_exp, pcov_exp = curve_fit(func_exp, ttd_zeroed, net_size_PiB, maxfev=1000)\n",
    "a_exp, b_exp, c_exp = popt_exp\n",
    "\n",
    "popt_sigmoid, pcov_sigmoid = curve_fit(func_sigmoid, ttd_zeroed, net_size_PiB, maxfev=1000)\n",
    "a_sigmoid, b_sigmoid, c_sigmoid, d_sigmoid = popt_sigmoid\n",
    "\n",
    "\n",
    "# Create prediction y-values\n",
    "yy_exp     = func_exp(tt_future_d_z, a_exp, b_exp, c_exp)\n",
    "yy_sigmoid = func_sigmoid(tt_future_d_z, a_sigmoid, b_sigmoid, c_sigmoid, d_sigmoid)\n",
    "\n",
    "\n",
    "# Get 365 day netsize prediction for every new point (or sparse set of time points?)\n",
    "ttd_yy_preds = np.zeros( int(len(ttd) / prediction_interval) + 1 )\n",
    "\n",
    "yy_preds = np.zeros(int(len(ttd) / prediction_interval) + 1)\n",
    "count = 0\n",
    "\n",
    "for i, t in enumerate(ttd):\n",
    "    if (i % prediction_interval == 0 and i > 50):\n",
    "        count = count + 1\n",
    "        try:\n",
    "            popt_sig, pcov_sig = curve_fit(func_sigmoid, ttd_zeroed[0:i], net_size_PiB[0:i], maxfev=1000)\n",
    "            a_sigmoid, b_sigmoid, c_sigmoid, d_sigmoid = popt_sig\n",
    "        \n",
    "            ttd_yy_preds[count] = ttd[i]\n",
    "            yy_preds[count] = func_sigmoid(tt_future_d_z[-1], a_sigmoid, b_sigmoid, c_sigmoid, d_sigmoid)\n",
    "        except:\n",
    "            print('All frigged up')\n",
    "            ttd_yy_preds[count] = ttd[i]\n",
    "            yy_preds[count] = yy_preds[count-1]\n",
    "            \n",
    "        #print('i: {}'.format(i))\n",
    "\n",
    "# use np.concatenate and np.full by chrisaycock\n",
    "def shift(arr, num, fill_value=np.nan):\n",
    "    if num >= 0:\n",
    "        return np.concatenate((np.full(num, fill_value), arr[:-num]))\n",
    "    else:\n",
    "        return np.concatenate((arr[-num:], np.full(-num, fill_value)))\n",
    "\n",
    "ttd_yy_preds[ttd_yy_preds == 0] = ttd_yy_preds[shift(ttd_yy_preds, 2) == 0]\n",
    "        \n",
    "# Set y limit to be a bit more than the sigmoid prediction (b/c exponential is ridiculous)\n",
    "y_max = max(yy_sigmoid) * 1.1\n",
    "\n",
    "# Print-outs\n",
    "print(\"tt_future_s[0]: {} \\t\\t tt_future_s[-1]: {}\".format(tt_future_s[0], tt_future_s[-1]))\n",
    "print(\"tt_future_d_z[0]: {} \\t\\t tt_future_d_z[-1]: {}\".format(tt_future_d_z[0], tt_future_d_z[-1]))\n",
    "print(\"tt_future_d[0]: {} \\t\\t tt_future_d[-1]: {}\".format(tt_future_d[0], tt_future_d[-1]))\n",
    "print()\n",
    "#print(\"yy_sigmoid[50:64]: {}\".format(yy_sigmoid[50:64]))\n",
    "#print(\"tt_future_d_z[50:64]: {}\".format(tt_future_d_z[50:64]))\n",
    "print()\n",
    "print(\"End of data at {} days\".format(ttd[-1]))\n",
    "print(\"Prediction ends at {} days\".format(tt_future_d[-1]))\n",
    "print(\"Value at {} days: {}\".format(tt_future_d[-1], \n",
    "                                    func_sigmoid(tt_future_d[-1], \n",
    "                                    a_sigmoid, b_sigmoid, c_sigmoid, d_sigmoid)))\n",
    "print()\n",
    "print(\"EXP:  {} * t^{} + {}\".format(a_exp, b_exp, c_exp))\n",
    "print(\"SIG:  ({} / ({}*e^-{}t)) + {}\".format(d_sigmoid, a_sigmoid, b_sigmoid, c_sigmoid, sigmoid_offset))\n",
    "print(sigmoid_offset)\n",
    "print()\n",
    "print('len(yy_preds): {}'.format(len(yy_preds)))\n",
    "#print('ttd_yy_preds:\\n {}'.format(ttd_yy_preds))\n",
    "#print('yy_preds:\\n {}'.format(yy_preds))\n",
    "print()\n",
    "print('ttd[500]: {}'.format(ttd[500]))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Prepare the figure\n",
    "fig_growth = plt.figure(num=None, figsize=(15, 30), dpi=80, facecolor='w', edgecolor='k')\n",
    "ax_now  = fig_growth.add_subplot(411)\n",
    "ax_zoom = fig_growth.add_subplot(412)\n",
    "ax_full = fig_growth.add_subplot(413)\n",
    "ax_log1 = fig_growth.add_subplot(414)\n",
    "\n",
    "\n",
    "# Common plotting/formatting\n",
    "t = time.strftime('%Y-%m-%d %H:%M:%S', time.localtime(epoch_times[0]))\n",
    "plt.suptitle(\"Netspace history & predictions. Day 0 is {}\".format(t))\n",
    "\n",
    "for ax in [ax_now, ax_zoom, ax_full, ax_log1]:\n",
    "    if PLOT_EXP: \n",
    "        ax.plot(tt_future_d, yy_exp, \n",
    "                color='blue', \n",
    "                linewidth=1,\n",
    "                label='Prediction curve [exponential]') \n",
    "        \n",
    "    if PLOT_SIG: \n",
    "        ax.plot(tt_future_d, yy_sigmoid, \n",
    "                color='green', \n",
    "                linewidth=2,\n",
    "                label='Prediction curve [sigmoid]')\n",
    "    \n",
    "    if PLOT_UPDATES:\n",
    "        for u in updates:\n",
    "            x = matplotlib.dates.date2num(u[0])\n",
    "            ax.axvline(x=x, linewidth=1, color='red', linestyle=':')\n",
    "            ax.text(x, 0, u[1])\n",
    "            \n",
    "    ax.plot(ttd, \n",
    "            net_size_PiB, \n",
    "            color='red', \n",
    "            linewidth=1,\n",
    "            label='Actual Netspace Size')\n",
    "    \n",
    "    ax.plot(ttd_yy_preds[:-2], \n",
    "            yy_preds[:-2], \n",
    "            color='orange', \n",
    "            linewidth=1,\n",
    "            linestyle='-',\n",
    "            label='Netspace Size in ~1yr [sigmoid]')\n",
    "    \n",
    "    ax.set_xlabel(\"Days\")\n",
    "    ax.set_ylabel(\"Netspace size [PiB]\")\n",
    "    ax.grid(which='major', axis='both', linestyle='-' , color='k')\n",
    "    ax.grid(which='minor', axis='both', linestyle='--')\n",
    "    ax.minorticks_on()\n",
    "    ax.xaxis.axis_date()\n",
    "    ax.legend()\n",
    "\n",
    "    \n",
    "# Unique plotting/formatting\n",
    "ax_full.text(tt_future_d[-1], yy_sigmoid[-1]*1.01, str(yy_sigmoid[-1]) + ' PiB', \n",
    "             fontweight='bold', \n",
    "             horizontalalignment='right')\n",
    "\n",
    "now_y_lo = min(net_size_PiB[-500:-1])\n",
    "now_y_hi = max(net_size_PiB[-500:-1]) + 100\n",
    "\n",
    "ax_now.set_xlim([ttd[-500],ttd[-1]+.5])\n",
    "ax_now.set_ylim([now_y_lo, now_y_hi])\n",
    "ax_now.set_title(\"Netsize: last 500 samples (500*5min = 2500min = 1day18hr)\")\n",
    "\n",
    "ax_zoom.set_xlim([ttd[0],ttd[-1]+5])\n",
    "ax_zoom.set_ylim([0,yy_preds[-2]*1.1])#net_size_PiB[-1]*1.1])\n",
    "ax_zoom.set_title(\"Netsize: Full History\")\n",
    "\n",
    "ax_full.set_ylim([0, y_max])\n",
    "ax_full.set_title(\"Netsize: 1 Year prediction\")\n",
    "\n",
    "ax_log1.semilogy()\n",
    "ax_log1.set_xlim([ttd[0],ttd[-1]+5])\n",
    "ax_log1.set_ylim([1 , net_size_PiB[-1]*1.1]) #net_size_PiB[0]\n",
    "ax_log1.set_title(\"Log scale base 10 means nothing\")\n",
    "ax_log1.legend()\n",
    "\n",
    "display(fig_growth)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Growth Rate\n",
    "\n",
    "\"Convolved\" == \"filtered\"\n",
    "\n",
    "Many things can be plotted here. What I have found useful to watch is listed in this chart. Enable plots for the others.\n",
    "\n",
    "|                          | netsize convolved [N]  |netsize convolved [Y]| netsize sigmoid prediction |\n",
    "| ---                      | :---                   | :---                | :---                       |\n",
    "| derivative convolved [N] | ---                    | cyan                | green                      |\n",
    "| derivative convolved [Y] | ~red~                  | blue                | ---                        |\n",
    "\n",
    "\n",
    "Change kernel size for more or less smoothing. Convolution delays the response by `kernel_size` data samples on the t-axis. This is why the convolved data trails when plotted. This is _sort of halfway_ corrected here.\n",
    "\n",
    "Here, I slide convolved plots for \"clock correction\". Currently, this is just the blue plot, where the netsize data was convolved. Then a derivative was taken. Then the derivative was convolved. So this is a smooth plot, that follows the cyan data, where the derivative was not \"re-convolved\". The cyan and blue should actually be slid to the left as well.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Convolution parameters\n",
    "days_to_smooth = 4\n",
    "samples_per_day = 24 * 60 / 5\n",
    "kernel_size = int(days_to_smooth * samples_per_day)\n",
    "kernel = np.ones(kernel_size) / kernel_size\n",
    "conv_mode = 'full'\n",
    "\n",
    "# Create derivative dealio\n",
    "net_size_PiB_dydt = np.zeros(len(net_size_PiB))\n",
    "sigmoid_dydt = np.zeros(len(yy_sigmoid))\n",
    "\n",
    "\n",
    "# Calculate derivative for real data\n",
    "for i, y in enumerate(net_size_PiB[:-2]):\n",
    "    net_size_PiB_dydt[i] = (net_size_PiB[i+1] - net_size_PiB[i]) / (ttd[i+1] - ttd[i])\n",
    "\n",
    "# Calculate derivative for sigmoid prediction\n",
    "for i, y in enumerate(yy_sigmoid[:-2]):\n",
    "    sigmoid_dydt[i] = (yy_sigmoid[i+1] - yy_sigmoid[i]) / (tt_future_d[i+1] - tt_future_d[i])\n",
    "\n",
    "    \n",
    "    \n",
    "# Convolve read data, then calculate derivative of it\n",
    "net_size_PiB_conv = np.convolve(net_size_PiB, kernel, mode=conv_mode)\n",
    "net_size_PiB_conv_dydt = np.zeros(len(net_size_PiB_conv))\n",
    "\n",
    "for i, y in enumerate(net_size_PiB_conv[:-2-kernel_size]):\n",
    "    net_size_PiB_conv_dydt[i] = (net_size_PiB_conv[i+1] - net_size_PiB_conv[i]) / (ttd[i+1] - ttd[i])\n",
    "\n",
    "\n",
    "    \n",
    "# Convolve the dydt of actual plot data\n",
    "net_size_PiB_dydt_convolved = np.convolve(net_size_PiB_dydt, kernel, mode=conv_mode)\n",
    "\n",
    "# Convolve the dydt of convolved actual data\n",
    "net_size_PiB_conv_dydt_conv = np.convolve(net_size_PiB_conv_dydt, kernel, mode=conv_mode)\n",
    "\n",
    "\n",
    "\n",
    "# Time arrays for plotting convolutions\n",
    "n_extra_points = len(net_size_PiB_conv_dydt_conv) - len(ttd)\n",
    "t_per_point = 1/24/12 # 5 minutes = 1d / 24hr / (60/5)\n",
    "\n",
    "ttd_conv_tail = np.arange(0, n_extra_points*t_per_point, t_per_point) + ttd[-1]\n",
    "ttd_conv = np.concatenate((ttd, ttd_conv_tail))\n",
    "ttd_conv_shft = ttd_conv - (kernel_size/2 * t_per_point)\n",
    "\n",
    "t_valid0_i = int(-(kernel_size/2))\n",
    "t_valid0 = ttd[t_valid0_i]\n",
    "t_valid1 = ttd[-1]\n",
    "t_valid_conv_shift_index = int(-n_extra_points-t_valid0_i-(kernel_size/2))\n",
    "\n",
    "print(kernel_size)\n",
    "print()\n",
    "print(\"n_extra_points: {}\".format(n_extra_points))\n",
    "print(\"t_per_point:    {}\".format(t_per_point))\n",
    "print()\n",
    "print(\"len(ttd):                         {}\".format(len(ttd)))\n",
    "print(\"ttd[0]: {:.2f} \\t\\t  ttd[-1]: {:.2f}\".format(ttd[0], ttd[-1]))\n",
    "print()\n",
    "print(\"len(tt_future_d):                 {}\".format(len(tt_future_d)))\n",
    "print(\"tt_future_d[0]: {:.2f} \\t  tt_future_d[-1]: {:.2f}\".format(tt_future_d[0], tt_future_d[-1]))\n",
    "print()\n",
    "print(\"len(ttd_conv_tail)                {}\".format(len(ttd_conv_tail)))\n",
    "print(\"ttd_conv_tail[0]: {:.2f} \\t  ttd_conv_tail[-1]: {:.2f}\".format(ttd_conv_tail[0], ttd_conv_tail[-1]))\n",
    "print()\n",
    "print(\"len(ttd_conv):                    {}\".format(len(ttd_conv)))\n",
    "print(\"ttd_conv[0]: {:.2f} \\t\\t  ttd_conv[-1]: {:.2f}\".format(ttd_conv[0], ttd_conv[-1]))\n",
    "print()\n",
    "print(\"len(ttd_conv_shft):               {}\".format(len(ttd_conv_shft)))\n",
    "print(\"ttd_conv_shft[0]: {:.2f} \\t  ttd_conv_shft[-1]: {:.2f}\".format(ttd_conv_shft[0], ttd_conv_shft[-1]))\n",
    "print()\n",
    "\n",
    "print(\"len(net_size_PiB_dydt_convolved): {}\".format(len(net_size_PiB_dydt_convolved)))\n",
    "print(\"len(net_size_PiB_conv_dydt_conv): {}\".format(len(net_size_PiB_conv_dydt_conv)))\n",
    "\n",
    "\n",
    "#print(\"net_size_PiB[0:4]: {}\".format(net_size_PiB[0:4]))\n",
    "#print(\"net_size_PiB_dydt[0:4]: {}\".format(net_size_PiB_dydt[0:4]))\n",
    "#print(\"yy_sigmoid[50:64]: {}\".format(yy_sigmoid[50:64]))\n",
    "#print(\"sigmoid_dydt[50:64]: {}\".format(sigmoid_dydt[50:64]))\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create the figure and subplot\n",
    "fig_growth_dydt = plt.figure(num=None, figsize=(12, 15), dpi=80, facecolor='w', edgecolor='k')\n",
    "ax_dydt = fig_growth_dydt.add_subplot(211)\n",
    "ax_shft = fig_growth_dydt.add_subplot(212)\n",
    "ax_netsize = ax_dydt.twinx()\n",
    "\n",
    "\n",
    "# Plotting\n",
    "for ax in [ax_dydt, ax_shft]:\n",
    "    \n",
    "    #ax.plot(ttd[:],\n",
    "    #        net_size_PiB[0:len(ttd)],\n",
    "    #        color='magenta',\n",
    "    #        linewidth=1,\n",
    "    #        label=\"Actual Netsize (Convolved)\")\n",
    "    \n",
    "    ax.plot(tt_future_d, sigmoid_dydt,\n",
    "            color='green',\n",
    "            linewidth=1, \n",
    "            label=\"Sigmoid Growth Rate Prediction\")\n",
    "\n",
    "    ax.plot(ttd[:],\n",
    "            net_size_PiB_conv_dydt[0:len(ttd)],\n",
    "            color='cyan',\n",
    "            linewidth=.2, \n",
    "            label=\"Netsize Growth Rate\")\n",
    "    \n",
    "    ax.plot(ttd_conv_shft[:t_valid_conv_shift_index],\n",
    "            net_size_PiB_conv_dydt_conv[:t_valid_conv_shift_index], \n",
    "            color='blue',  \n",
    "            linewidth=1,\n",
    "            label=\"Netsize Growth Rate (Convolved)\")\n",
    "    \n",
    "    ax.plot(ttd_conv_shft[t_valid_conv_shift_index:],\n",
    "            net_size_PiB_conv_dydt_conv[t_valid_conv_shift_index:], \n",
    "            color='blue',  \n",
    "            linewidth=1,\n",
    "            linestyle=':',\n",
    "            label=\"Netsize Growth Rate (Convolved, invalid)\")\n",
    "    \n",
    "    # Common formatting\n",
    "    ax.xaxis.axis_date()\n",
    "    ax.grid()\n",
    "    ax.set_xlabel(\"Date\")\n",
    "    ax.set_ylabel(\"PiB\")\n",
    "\n",
    "\n",
    "    \n",
    "ax_netsize.plot(ttd[:],\n",
    "                net_size_PiB[0:len(ttd)],\n",
    "                color='grey',\n",
    "                linewidth=.5,\n",
    "                label=\"Actual Netsize (Convolved)\")\n",
    "    \n",
    "# Unique formatting\n",
    "ax_dydt.set_title(\"Netspace growth [PiB/day]\")\n",
    "ax_dydt.set_ylim([min(net_size_PiB_conv_dydt[600:])*1.1, max(net_size_PiB_conv_dydt[600:])*1.1])\n",
    "ax_dydt.set_xlim([ttd[0], ttd[-1]+10])\n",
    "ax_dydt.legend()\n",
    "ax_dydt.set_ylabel(\"PiB (Growth Rate)\")\n",
    "ax_netsize.set_ylabel(\"PiB (Netsize)\")\n",
    "\n",
    "ax_shft.set_title(\"Netspace growth [PiB/day], with clock correction\")\n",
    "ax_shft.set_xlim([ttd[250], ttd_conv_shft[-1]])\n",
    "ax_shft.set_ylim([-10, 1000])\n",
    "ax_shft.legend()\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "# Update the legend for both axes in subplot 1\n",
    "handles, labels = ax_shft.get_legend_handles_labels()\n",
    "patch = Line2D([0], [0], label = 'Actual Netsize',color='grey')\n",
    "handles.append(patch)\n",
    "plt.legend(handles=handles)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Insert: Acceleration of growth"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sigmoid_accel = np.zeros(len(yy_sigmoid))\n",
    "netsize_accel = np.zeros(len(net_size_PiB_conv_dydt_conv))\n",
    "\n",
    "\n",
    "# Calc d2y/dt2 for sigmoid\n",
    "for i, y in enumerate(sigmoid_dydt[:-2]):\n",
    "    sigmoid_accel[i] = (sigmoid_dydt[i+1] - sigmoid_dydt[i]) / (tt_future_d[i+1] - tt_future_d[i])\n",
    "\n",
    "\n",
    "for i, y in enumerate(net_size_PiB_conv_dydt_conv[:-2]):\n",
    "    netsize_accel[i] = (net_size_PiB_conv_dydt_conv[i+1] - net_size_PiB_conv_dydt_conv[i]) / \\\n",
    "                       (ttd_conv[i+1] - ttd_conv[i])\n",
    "\n",
    "# filter for invalid values, then convolve\n",
    "netsize_accel[netsize_accel < -10000] = 0\n",
    "netsize_accel_conv = np.convolve(netsize_accel, kernel, mode=conv_mode)\n",
    "    \n",
    "# Time arrays for plotting convolutions\n",
    "n_extra_points = len(netsize_accel_conv) - len(ttd)\n",
    "t_per_point = 1/24/12 # 5 minutes = 1d / 24hr / (60/5)\n",
    "\n",
    "ttd_conv_tail2 = np.arange(0, n_extra_points*t_per_point, t_per_point) + ttd[-1]\n",
    "ttd_conv2 = np.concatenate((ttd, ttd_conv_tail2))\n",
    "ttd_conv_shft2 = ttd_conv2 - (kernel_size/2 * t_per_point)\n",
    "\n",
    "t_valid0_i = int(-(kernel_size/2))\n",
    "t_valid_conv_shift_index = int(-n_extra_points-t_valid0_i-(kernel_size/2))\n",
    "\n",
    "print(\"len(ttd): {}\".format(len(ttd)))\n",
    "print(\"len(sigmoid_dydt):    {}\".format(len(sigmoid_dydt)))\n",
    "print(\"len(sigmoid_accel):   {}\".format(len(sigmoid_accel)))\n",
    "print(\"len(netsize_accel):   {}\".format(len(netsize_accel)))\n",
    "print(\"len(net_size_PiB_conv_dydt_conv): {}\".format(len(net_size_PiB_conv_dydt_conv)))\n",
    "print(\"len(netsize_accel_conv): {}\".format(len(netsize_accel_conv)))\n",
    "print(\"len(ttd_conv_shft2): {}\".format(len(ttd_conv_shft2)))\n",
    "\n",
    "print(\"t_valid_conv_shift_index: {}\".format(t_valid_conv_shift_index))\n",
    "print()\n",
    "print(\"ttd_conv_shft2[t_valid_conv_shift_index-5:t_valid_conv_shift_index+5]: \\n\\t{}\".format(ttd_conv_shft2[t_valid_conv_shift_index-5:t_valid_conv_shift_index+5]))\n",
    "print()\n",
    "print(\"netsize_accel_conv[t_valid_conv_shift_index-5:t_valid_conv_shift_index+5]: \\n\\t{}\".format(netsize_accel_conv[t_valid_conv_shift_index-5:t_valid_conv_shift_index+5]))\n",
    "print()\n",
    "print(\"netsize_accel_conv[t_valid_conv_shift_index-5+kernel_size:t_valid_conv_shift_index+5+kernel_size]: \\n\\t{}\".format(netsize_accel_conv[t_valid_conv_shift_index-5+kernel_size:t_valid_conv_shift_index+5+kernel_size]))\n",
    "print()\n",
    "print(\"netsize_accel[t_valid_conv_shift_index-5+kernel_size:t_valid_conv_shift_index+5+kernel_size]: \\n\\t{}\".format(netsize_accel[t_valid_conv_shift_index-5+kernel_size:t_valid_conv_shift_index+5+kernel_size]))\n",
    "# ttd_conv_shft2[:t_valid_conv_shift_index]\n",
    "# netsize_accel_conv[:t_valid_conv_shift_index]\n",
    "#netsize_accel\n",
    "\n",
    "print(netsize_accel[(netsize_accel < -10000)])\n",
    "print(netsize_accel_conv[(netsize_accel_conv < -10000)])\n",
    "print(netsize_accel_conv[(netsize_accel_conv > 10000)])\n",
    "\n",
    "print(type(ttd_conv_shft))\n",
    "print(type(ttd_conv_shft[0]))\n",
    "for i, a in enumerate(ttd_conv_shft2[:-2]):\n",
    "    #if ttd_conv_shft2[i+1] - ttd_conv_shft2[i] > 0.03:\n",
    "    #    print(i)\n",
    "    if not isinstance(a, np.float64):\n",
    "        print(i)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create the figure and subplot\n",
    "fig_growth_accel = plt.figure(num=None, figsize=(12, 8), dpi=80, facecolor='w', edgecolor='k')\n",
    "ax_accel = fig_growth_accel.add_subplot(111)\n",
    "ax_netsize2 = ax_accel.twinx()\n",
    "\n",
    "\n",
    "# Plot\n",
    "ax_accel.plot(tt_future_d, \n",
    "              sigmoid_accel,\n",
    "              color='green',\n",
    "              linewidth=1, \n",
    "              label=\"Sigmoid Acceleration Prediction\")\n",
    "\n",
    "ax_accel.plot(ttd_conv[0:len(netsize_accel)], \n",
    "              netsize_accel,\n",
    "              color='cyan',\n",
    "              linewidth=.2,\n",
    "              label=\"Netsize Acceleration\")\n",
    "\n",
    "ax_accel.plot(ttd_conv_shft2[0:int(t_valid_conv_shift_index/2 - kernel_size/2)], \n",
    "              netsize_accel_conv[0:int(t_valid_conv_shift_index/2 - kernel_size/2)],\n",
    "              color='blue',\n",
    "              linewidth=1,\n",
    "              label=\"Netsize Acceleration - Smoothed\")\n",
    "\n",
    "ax_accel.plot(ttd_conv_shft2[t_valid_conv_shift_index+kernel_size:], \n",
    "              netsize_accel_conv[t_valid_conv_shift_index+kernel_size:],\n",
    "              color='blue',\n",
    "              linestyle=':',\n",
    "              linewidth=1,\n",
    "              label=\"Netsize Acceleration - Smoothed (invalid)\")\n",
    "\n",
    "ax_netsize2.plot(ttd[:],\n",
    "                net_size_PiB[0:len(ttd)],\n",
    "                color='grey',\n",
    "                linewidth=.5,\n",
    "                label=\"Actual Netsize (Convolved)\")\n",
    "\n",
    "# Formatting\n",
    "ax_accel.xaxis.axis_date()\n",
    "ax_accel.grid()\n",
    "ax_accel.set_xlim([ttd[250], ttd_conv_shft[-1]])\n",
    "ax_accel.set_ylim([-500, 500])\n",
    "ax_accel.legend()\n",
    "\n",
    "ax_netsize2.set_ylabel(\"PiB (Netsize)\")\n",
    "\n",
    "# Update the legend for both axes in subplot 1\n",
    "handles, labels = ax_accel.get_legend_handles_labels()\n",
    "patch = Line2D([0], [0], label = 'Actual Netsize',color='grey')\n",
    "handles.append(patch)\n",
    "plt.legend(handles=handles)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### PMF - Create a matrix with win percentages from binomial pmf blah blah stats bullshit\n",
    "\n",
    "Change 'r' variable to see probabilities for winning r times\n",
    "e.g. to see a distribution for likelihood of winning 6 times, set r = 6"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "# n = number of trials\n",
    "# p = probability of winning\n",
    "# r = # of wins \n",
    "\n",
    "pmf_dist_days = 365\n",
    "challenges_per_day = 4608\n",
    "num_data_points= len(ttd)\n",
    "r = 1\n",
    "win_percents = np.zeros((num_data_points, pmf_dist_days))\n",
    "challenges_range = np.arange(challenges_per_day, pmf_dist_days*challenges_per_day, challenges_per_day)\n",
    "\n",
    "for i in range(0, num_data_points):\n",
    "    p = win_prob[i]\n",
    "    \n",
    "    for j, n in enumerate(challenges_range):\n",
    "        win_percents[i][j] = binom.pmf(r, n, p)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plot PMFs across real time"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"num_data_points: {}\".format(num_data_points))\n",
    "print(\"days: {}\".format(pmf_dist_days))\n",
    "\n",
    "# Get range of max probabilty of win indices, then set color range up to that value\n",
    "max_day_index = 0\n",
    "for sample in range(0, num_data_points):    \n",
    "    m = win_percents[sample].argmax()\n",
    "    if (m > max_day_index):\n",
    "        max_day_index = m\n",
    "    \n",
    "print(\"max_day_index: {}\".format(max_day_index))\n",
    "\n",
    "green = Color(\"green\")\n",
    "colors_year = list(green.range_to(Color(\"red\"), max_day_index+1))\n",
    "\n",
    "\n",
    "fig_pmf = plt.figure(num=None, figsize=(12, 12), dpi=80, facecolor='w', edgecolor='k')\n",
    "ax_pmf = fig_pmf.add_subplot(111, projection='3d')\n",
    "\n",
    "x_array = np.ones(pmf_dist_days)\n",
    "max_prob_indices = np.zeros(num_data_points)\n",
    "\n",
    "print(win_percents.shape)\n",
    "\n",
    "line_width = 100/len(ttd) #0.1\n",
    "print(\"line_width: {}\".format(line_width))\n",
    "\n",
    "print(\"ttd[0]: {}\".format(ttd[0]))\n",
    "\n",
    "for sample in range(0, num_data_points):    \n",
    "    max_prob_indices[sample] = win_percents[sample].argmax()\n",
    "    \n",
    "    ax_pmf.plot(ttd[sample]*x_array, np.arange(0, pmf_dist_days), win_percents[sample][:],\n",
    "                color=colors_year[int(max_prob_indices[sample])].hex_l,\n",
    "                linewidth= line_width\n",
    "               )\n",
    "    ax_pmf.plot([ttd[sample]], max_prob_indices[sample], 0,\n",
    "                marker='o',\n",
    "                markersize=2,\n",
    "                color=colors_year[int(max_prob_indices[sample])].hex_l\n",
    "    )\n",
    "    \n",
    "elev = 19 #12       \n",
    "angle = -1 #-12\n",
    "ax_pmf.view_init(elev, angle)   \n",
    "ax_pmf.set_title(\"\")\n",
    "ax_pmf.set_ylabel(\"Days til win\")\n",
    "ax_pmf.set_xlabel(\"Plotting Date\\n\\n\\n\")\n",
    "ax_pmf.set_zlabel(\"Chance of win\")\n",
    "ax_pmf.xaxis.axis_date()\n",
    "ax_pmf.grid(which='major', axis='y', linestyle='-')\n",
    "ax_pmf.grid(which='minor', axis='y', linestyle='--', bottom='off')\n",
    "#ax_pmf.yaxis.minorticks_on()\n",
    "\n",
    "#ax_pmf.set_xlim([58, ttd[-1]])\n",
    "ax_pmf.set_xlim([ttd[first_sample_index], ttd[-1]])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Insert: PMF plot with meshgrid style"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"coming soon!\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plot fraction of owned netspace"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "netfrac_fig = plt.figure(num=None, figsize=(12, 6), dpi=80, facecolor='w', edgecolor='k')\n",
    "netfrac_ax = netfrac_fig.add_subplot(111)\n",
    "\n",
    "#netfrac_ax.plot(xx[219:], win_prob[219:])\n",
    "\n",
    "\n",
    "netfrac_ax.plot(ttd[219:], win_prob[219:])\n",
    "netfrac_ax.xaxis.axis_date()\n",
    "\n",
    "netfrac_ax.set_title(\"Fraction of Netspace Owned\")\n",
    "netfrac_ax.set_xlabel(\"Days\")\n",
    "netfrac_ax.set_ylabel(\"Fraction Owned\")\n",
    "\n",
    "plt.grid(which='major', axis='both', linestyle='-' , color='k')\n",
    "plt.grid(which='minor', axis='both', linestyle='--')\n",
    "plt.minorticks_on()\n",
    "\n",
    "#netfrac_ax.set_xlim([58, xx[-1]])\n",
    "\n",
    "#for i,n in enumerate(xx):\n",
    "#    if xx[i+1] < xx[i]:\n",
    "#        print(i+1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "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.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}