ksample/ksamp_power_angle.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sys, os\n",
    "import multiprocessing as mp\n",
    "from joblib import Parallel, delayed\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from power_ksamp import power_ksamp_angle\n",
    "from hyppo.independence import CCA, Dcorr, HHG, Hsic, RV, MGC, KMERF\n",
    "from hyppo.tools import *\n",
    "\n",
    "sys.path.append(os.path.realpath('..'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import seaborn as sns\n",
    "sns.set(color_codes=True, style='white', context='talk', font_scale=1.5)\n",
    "PALETTE = sns.color_palette(\"Set1\")\n",
    "sns.set_palette(PALETTE[1:5] + PALETTE[6:])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "These are some constants that are used in this notebook. If running these notebook, please only manipulate these constants if you are not running more tests. They define the sample sizes tested upon and the number of replications. The simulations tested over and the independence tests tested over are defined also."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "from rpy2.robjects import Formula, numpy2ri\n",
    "from rpy2.robjects.packages import importr\n",
    "\n",
    "\n",
    "class Manova:\n",
    "    r\"\"\"\n",
    "    Wrapper of R MANOVA\n",
    "    \"\"\"\n",
    "    def __init__(self):\n",
    "        self.stats = importr('stats')\n",
    "        self.r_base = importr('base')\n",
    "        \n",
    "        numpy2ri.activate()\n",
    "\n",
    "        self.formula = Formula('X ~ Y')\n",
    "        self.env = self.formula.environment\n",
    "\n",
    "    def statistic(self, x, y):\n",
    "        r\"\"\"\n",
    "        Helper function to calculate the test statistic\n",
    "        \"\"\"\n",
    "        self.env['Y'] = y\n",
    "        self.env['X'] = x\n",
    "\n",
    "        stat = self.r_base.summary(self.stats.manova(self.formula), test=\"Pillai\")[3][0, 1]\n",
    "\n",
    "        return stat"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "MAX_ANGLE = 90\n",
    "STEP_SIZE = 5\n",
    "ANGLES = range(1, MAX_ANGLE + STEP_SIZE, STEP_SIZE)\n",
    "POWER_REPS = 5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "SIMULATIONS = {\n",
    "    \"linear\": \"Linear\",\n",
    "    \"exponential\": \"Exponential\",\n",
    "    \"cubic\": \"Cubic\",\n",
    "    \"joint_normal\": \"Joint Normal\",\n",
    "    \"step\": \"Step\",\n",
    "    \"quadratic\": \"Quadratic\",\n",
    "    \"w_shaped\": \"W-Shaped\",\n",
    "    \"spiral\": \"Spiral\",\n",
    "    \"uncorrelated_bernoulli\": \"Bernoulli\",\n",
    "    \"logarithmic\": \"Logarithmic\",\n",
    "    \"fourth_root\": \"Fourth Root\",\n",
    "    \"sin_four_pi\": \"Sine 4\\u03C0\",\n",
    "    \"sin_sixteen_pi\": \"Sine 16\\u03C0\",\n",
    "    \"square\": \"Square\",\n",
    "    \"two_parabolas\": \"Two Parabolas\",\n",
    "    \"circle\": \"Circle\",\n",
    "    \"ellipse\": \"Ellipse\",\n",
    "    \"diamond\": \"Diamond\",\n",
    "    \"multiplicative_noise\": \"Multiplicative\",\n",
    "    \"multimodal_independence\": \"Independence\"\n",
    "}\n",
    "\n",
    "TESTS = [\n",
    "    KMERF,\n",
    "    MGC,\n",
    "    Dcorr,\n",
    "    Hsic,\n",
    "    Manova,\n",
    "    HHG,\n",
    "    CCA,\n",
    "    RV,\n",
    "]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following function calculates the estimated power ``POWER_REPS`` number off times and averages them. It does this iterating over the number of sample sizes.\n",
    "\n",
    "**Note: We only recommend running this code if running the next 2 cells ONCE to generate the csv files used to visualize the plots. This code takes a very long time to run and if running, we recommend using a machine with many cores.**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def estimate_power(sim, test):\n",
    "    est_power = np.array([np.mean([power_ksamp_angle(test, rot_ksamp, sim, angle=i) for _ in range(POWER_REPS)])\n",
    "                          for i in ANGLES])\n",
    "    np.savetxt('../ksample/ksamp_vs_angle/{}_{}.csv'.format(sim, test.__name__),\n",
    "               est_power, delimiter=',')\n",
    "    \n",
    "    return est_power"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[Parallel(n_jobs=-1)]: Using backend LokyBackend with 16 concurrent workers.\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Unable to determine R library path: Command '('/Library/Frameworks/R.framework/Resources/bin/Rscript', '-e', 'cat(Sys.getenv(\"LD_LIBRARY_PATH\"))')' returned non-zero exit status 2.\n",
      "Unable to determine R library path: Command '('/Library/Frameworks/R.framework/Resources/bin/Rscript', '-e', 'cat(Sys.getenv(\"LD_LIBRARY_PATH\"))')' returned non-zero exit status 2.\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[Parallel(n_jobs=-1)]: Done   1 tasks      | elapsed:  2.2min\n",
      "[Parallel(n_jobs=-1)]: Done   2 tasks      | elapsed:  2.2min\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Unable to determine R library path: Command '('/Library/Frameworks/R.framework/Resources/bin/Rscript', '-e', 'cat(Sys.getenv(\"LD_LIBRARY_PATH\"))')' returned non-zero exit status 2.\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[Parallel(n_jobs=-1)]: Done   3 tasks      | elapsed:  2.6min\n",
      "[Parallel(n_jobs=-1)]: Done   4 tasks      | elapsed:  2.7min\n",
      "[Parallel(n_jobs=-1)]: Done   5 tasks      | elapsed:  5.4min\n",
      "[Parallel(n_jobs=-1)]: Done   6 tasks      | elapsed:  7.9min\n",
      "[Parallel(n_jobs=-1)]: Done   7 tasks      | elapsed: 22.0min\n",
      "[Parallel(n_jobs=-1)]: Done   8 tasks      | elapsed: 22.1min\n",
      "[Parallel(n_jobs=-1)]: Done   9 tasks      | elapsed: 25.4min\n",
      "[Parallel(n_jobs=-1)]: Done  10 tasks      | elapsed: 31.1min\n",
      "[Parallel(n_jobs=-1)]: Done  11 tasks      | elapsed: 31.3min\n",
      "[Parallel(n_jobs=-1)]: Done  12 tasks      | elapsed: 31.4min\n",
      "[Parallel(n_jobs=-1)]: Done  13 tasks      | elapsed: 34.8min\n",
      "[Parallel(n_jobs=-1)]: Done  14 tasks      | elapsed: 35.1min\n",
      "[Parallel(n_jobs=-1)]: Done  15 tasks      | elapsed: 49.3min\n",
      "[Parallel(n_jobs=-1)]: Done  16 tasks      | elapsed: 52.3min\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Unable to determine R library path: Command '('/Library/Frameworks/R.framework/Resources/bin/Rscript', '-e', 'cat(Sys.getenv(\"LD_LIBRARY_PATH\"))')' returned non-zero exit status 2.\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[Parallel(n_jobs=-1)]: Done  17 tasks      | elapsed: 52.4min\n",
      "[Parallel(n_jobs=-1)]: Done  18 tasks      | elapsed: 52.6min\n",
      "[Parallel(n_jobs=-1)]: Done  19 tasks      | elapsed: 54.4min\n",
      "[Parallel(n_jobs=-1)]: Done  20 tasks      | elapsed: 55.6min\n",
      "[Parallel(n_jobs=-1)]: Done  21 tasks      | elapsed: 57.0min\n",
      "[Parallel(n_jobs=-1)]: Done  22 tasks      | elapsed: 60.7min\n",
      "[Parallel(n_jobs=-1)]: Done  23 tasks      | elapsed: 70.4min\n",
      "[Parallel(n_jobs=-1)]: Done  24 tasks      | elapsed: 71.9min\n",
      "[Parallel(n_jobs=-1)]: Done  25 tasks      | elapsed: 78.1min\n",
      "[Parallel(n_jobs=-1)]: Done  26 tasks      | elapsed: 80.9min\n",
      "[Parallel(n_jobs=-1)]: Done  27 tasks      | elapsed: 83.3min\n",
      "[Parallel(n_jobs=-1)]: Done  28 tasks      | elapsed: 88.7min\n",
      "[Parallel(n_jobs=-1)]: Done  29 tasks      | elapsed: 89.0min\n",
      "[Parallel(n_jobs=-1)]: Done  30 tasks      | elapsed: 98.5min\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Unable to determine R library path: Command '('/Library/Frameworks/R.framework/Resources/bin/Rscript', '-e', 'cat(Sys.getenv(\"LD_LIBRARY_PATH\"))')' returned non-zero exit status 2.\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[Parallel(n_jobs=-1)]: Done  31 tasks      | elapsed: 106.8min\n",
      "[Parallel(n_jobs=-1)]: Done  32 tasks      | elapsed: 114.8min\n",
      "[Parallel(n_jobs=-1)]: Done  33 tasks      | elapsed: 117.3min\n",
      "[Parallel(n_jobs=-1)]: Done  34 tasks      | elapsed: 119.5min\n",
      "[Parallel(n_jobs=-1)]: Done  35 tasks      | elapsed: 122.7min\n",
      "[Parallel(n_jobs=-1)]: Done  36 tasks      | elapsed: 139.0min\n",
      "[Parallel(n_jobs=-1)]: Done  37 tasks      | elapsed: 171.3min\n",
      "[Parallel(n_jobs=-1)]: Done  38 tasks      | elapsed: 186.9min\n",
      "[Parallel(n_jobs=-1)]: Done  39 tasks      | elapsed: 190.2min\n",
      "[Parallel(n_jobs=-1)]: Done  40 tasks      | elapsed: 190.5min\n",
      "[Parallel(n_jobs=-1)]: Done  41 tasks      | elapsed: 191.2min\n",
      "[Parallel(n_jobs=-1)]: Done  42 tasks      | elapsed: 192.6min\n",
      "[Parallel(n_jobs=-1)]: Done  43 tasks      | elapsed: 192.7min\n",
      "[Parallel(n_jobs=-1)]: Done  44 tasks      | elapsed: 199.6min\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Unable to determine R library path: Command '('/Library/Frameworks/R.framework/Resources/bin/Rscript', '-e', 'cat(Sys.getenv(\"LD_LIBRARY_PATH\"))')' returned non-zero exit status 2.\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[Parallel(n_jobs=-1)]: Done  45 tasks      | elapsed: 209.5min\n",
      "[Parallel(n_jobs=-1)]: Done  46 tasks      | elapsed: 215.7min\n",
      "[Parallel(n_jobs=-1)]: Done  47 tasks      | elapsed: 216.4min\n",
      "[Parallel(n_jobs=-1)]: Done  48 tasks      | elapsed: 218.6min\n",
      "[Parallel(n_jobs=-1)]: Done  49 tasks      | elapsed: 218.9min\n",
      "[Parallel(n_jobs=-1)]: Done  50 tasks      | elapsed: 225.0min\n",
      "[Parallel(n_jobs=-1)]: Done  51 tasks      | elapsed: 233.7min\n",
      "[Parallel(n_jobs=-1)]: Done  52 tasks      | elapsed: 235.4min\n",
      "[Parallel(n_jobs=-1)]: Done  53 tasks      | elapsed: 248.9min\n",
      "[Parallel(n_jobs=-1)]: Done  54 tasks      | elapsed: 251.7min\n",
      "[Parallel(n_jobs=-1)]: Done  55 tasks      | elapsed: 253.8min\n",
      "[Parallel(n_jobs=-1)]: Done  56 tasks      | elapsed: 254.1min\n",
      "[Parallel(n_jobs=-1)]: Done  57 tasks      | elapsed: 256.4min\n",
      "[Parallel(n_jobs=-1)]: Done  58 tasks      | elapsed: 269.7min\n",
      "[Parallel(n_jobs=-1)]: Done  59 tasks      | elapsed: 284.3min\n",
      "[Parallel(n_jobs=-1)]: Done  60 tasks      | elapsed: 287.1min\n",
      "[Parallel(n_jobs=-1)]: Done  61 tasks      | elapsed: 289.5min\n",
      "[Parallel(n_jobs=-1)]: Done  62 tasks      | elapsed: 290.1min\n",
      "[Parallel(n_jobs=-1)]: Done  63 tasks      | elapsed: 291.6min\n",
      "[Parallel(n_jobs=-1)]: Done  64 tasks      | elapsed: 307.0min\n",
      "[Parallel(n_jobs=-1)]: Done  65 tasks      | elapsed: 337.6min\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Unable to determine R library path: Command '('/Library/Frameworks/R.framework/Resources/bin/Rscript', '-e', 'cat(Sys.getenv(\"LD_LIBRARY_PATH\"))')' returned non-zero exit status 2.\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[Parallel(n_jobs=-1)]: Done  66 tasks      | elapsed: 352.2min\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "R[write to console]: There were 50 or more warnings (use warnings() to see the first 50)\n",
      "R[write to console]: \n",
      "\n",
      "/opt/miniconda3/envs/hyppo-dev/lib/python3.10/site-packages/joblib/externals/loky/process_executor.py:702: UserWarning: A worker stopped while some jobs were given to the executor. This can be caused by a too short worker timeout or by a memory leak.\n",
      "  warnings.warn(\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[Parallel(n_jobs=-1)]: Done  67 tasks      | elapsed: 354.9min\n",
      "[Parallel(n_jobs=-1)]: Done  68 tasks      | elapsed: 357.6min\n",
      "[Parallel(n_jobs=-1)]: Done  69 tasks      | elapsed: 360.0min\n",
      "[Parallel(n_jobs=-1)]: Done  70 tasks      | elapsed: 383.6min\n",
      "[Parallel(n_jobs=-1)]: Done  71 tasks      | elapsed: 398.9min\n",
      "[Parallel(n_jobs=-1)]: Done  72 tasks      | elapsed: 412.3min\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Unable to determine R library path: Command '('/Library/Frameworks/R.framework/Resources/bin/Rscript', '-e', 'cat(Sys.getenv(\"LD_LIBRARY_PATH\"))')' returned non-zero exit status 2.\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[Parallel(n_jobs=-1)]: Done  73 tasks      | elapsed: 426.8min\n",
      "[Parallel(n_jobs=-1)]: Done  74 tasks      | elapsed: 429.4min\n",
      "[Parallel(n_jobs=-1)]: Done  75 tasks      | elapsed: 432.1min\n",
      "[Parallel(n_jobs=-1)]: Done  76 tasks      | elapsed: 434.5min\n",
      "[Parallel(n_jobs=-1)]: Done  77 tasks      | elapsed: 450.5min\n",
      "[Parallel(n_jobs=-1)]: Done  78 tasks      | elapsed: 465.4min\n",
      "[Parallel(n_jobs=-1)]: Done  79 tasks      | elapsed: 495.3min\n",
      "[Parallel(n_jobs=-1)]: Done  80 tasks      | elapsed: 509.4min\n",
      "[Parallel(n_jobs=-1)]: Done  81 tasks      | elapsed: 515.8min\n",
      "[Parallel(n_jobs=-1)]: Done  82 tasks      | elapsed: 518.5min\n",
      "[Parallel(n_jobs=-1)]: Done  83 tasks      | elapsed: 520.9min\n",
      "[Parallel(n_jobs=-1)]: Done  84 tasks      | elapsed: 591.7min\n",
      "[Parallel(n_jobs=-1)]: Done  85 tasks      | elapsed: 606.6min\n",
      "[Parallel(n_jobs=-1)]: Done  86 tasks      | elapsed: 658.2min\n",
      "[Parallel(n_jobs=-1)]: Done  87 tasks      | elapsed: 673.4min\n",
      "[Parallel(n_jobs=-1)]: Done  88 tasks      | elapsed: 699.4min\n",
      "[Parallel(n_jobs=-1)]: Done  89 tasks      | elapsed: 701.8min\n",
      "[Parallel(n_jobs=-1)]: Done  90 tasks      | elapsed: 703.8min\n",
      "[Parallel(n_jobs=-1)]: Done  91 tasks      | elapsed: 892.5min\n",
      "[Parallel(n_jobs=-1)]: Done  92 tasks      | elapsed: 909.1min\n",
      "[Parallel(n_jobs=-1)]: Done  93 tasks      | elapsed: 943.3min\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Unable to determine R library path: Command '('/Library/Frameworks/R.framework/Resources/bin/Rscript', '-e', 'cat(Sys.getenv(\"LD_LIBRARY_PATH\"))')' returned non-zero exit status 2.\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[Parallel(n_jobs=-1)]: Done  94 tasks      | elapsed: 959.2min\n",
      "[Parallel(n_jobs=-1)]: Done  95 tasks      | elapsed: 1415.7min\n",
      "[Parallel(n_jobs=-1)]: Done  96 tasks      | elapsed: 1416.7min\n",
      "[Parallel(n_jobs=-1)]: Done  97 tasks      | elapsed: 1418.3min\n",
      "[Parallel(n_jobs=-1)]: Done  98 tasks      | elapsed: 1419.0min\n",
      "[Parallel(n_jobs=-1)]: Done  99 tasks      | elapsed: 1423.3min\n",
      "[Parallel(n_jobs=-1)]: Done 100 tasks      | elapsed: 1433.8min\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Unable to determine R library path: Command '('/Library/Frameworks/R.framework/Resources/bin/Rscript', '-e', 'cat(Sys.getenv(\"LD_LIBRARY_PATH\"))')' returned non-zero exit status 2.\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[Parallel(n_jobs=-1)]: Done 101 tasks      | elapsed: 1444.5min\n",
      "[Parallel(n_jobs=-1)]: Done 102 tasks      | elapsed: 1447.8min\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "R[write to console]: There were 50 or more warnings (use warnings() to see the first 50)\n",
      "R[write to console]: \n",
      "\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[Parallel(n_jobs=-1)]: Done 103 tasks      | elapsed: 1450.5min\n",
      "[Parallel(n_jobs=-1)]: Done 104 tasks      | elapsed: 1452.8min\n",
      "[Parallel(n_jobs=-1)]: Done 105 tasks      | elapsed: 1452.8min\n",
      "[Parallel(n_jobs=-1)]: Done 106 tasks      | elapsed: 1466.6min\n",
      "[Parallel(n_jobs=-1)]: Done 107 tasks      | elapsed: 1467.4min\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Unable to determine R library path: Command '('/Library/Frameworks/R.framework/Resources/bin/Rscript', '-e', 'cat(Sys.getenv(\"LD_LIBRARY_PATH\"))')' returned non-zero exit status 2.\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[Parallel(n_jobs=-1)]: Done 108 tasks      | elapsed: 1481.6min\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "R[write to console]: There were 50 or more warnings (use warnings() to see the first 50)\n",
      "R[write to console]: \n",
      "\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[Parallel(n_jobs=-1)]: Done 109 tasks      | elapsed: 1496.4min\n",
      "[Parallel(n_jobs=-1)]: Done 111 out of 140 | elapsed: 1499.1min remaining: 391.7min\n",
      "[Parallel(n_jobs=-1)]: Done 113 out of 140 | elapsed: 1515.2min remaining: 362.0min\n",
      "[Parallel(n_jobs=-1)]: Done 115 out of 140 | elapsed: 1545.3min remaining: 335.9min\n",
      "[Parallel(n_jobs=-1)]: Done 117 out of 140 | elapsed: 1549.8min remaining: 304.7min\n",
      "[Parallel(n_jobs=-1)]: Done 119 out of 140 | elapsed: 1572.0min remaining: 277.4min\n",
      "[Parallel(n_jobs=-1)]: Done 121 out of 140 | elapsed: 1592.6min remaining: 250.1min\n",
      "[Parallel(n_jobs=-1)]: Done 123 out of 140 | elapsed: 1608.4min remaining: 222.3min\n",
      "[Parallel(n_jobs=-1)]: Done 125 out of 140 | elapsed: 1613.9min remaining: 193.7min\n",
      "[Parallel(n_jobs=-1)]: Done 127 out of 140 | elapsed: 1618.9min remaining: 165.7min\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "R[write to console]: There were 50 or more warnings (use warnings() to see the first 50)\n",
      "R[write to console]: \n",
      "\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[Parallel(n_jobs=-1)]: Done 129 out of 140 | elapsed: 1648.6min remaining: 140.6min\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "R[write to console]: There were 50 or more warnings (use warnings() to see the first 50)\n",
      "R[write to console]: \n",
      "\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[Parallel(n_jobs=-1)]: Done 131 out of 140 | elapsed: 1686.7min remaining: 115.9min\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "R[write to console]: There were 50 or more warnings (use warnings() to see the first 50)\n",
      "R[write to console]: \n",
      "\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[Parallel(n_jobs=-1)]: Done 133 out of 140 | elapsed: 1759.0min remaining: 92.6min\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "R[write to console]: There were 50 or more warnings (use warnings() to see the first 50)\n",
      "R[write to console]: \n",
      "\n",
      "R[write to console]: There were 50 or more warnings (use warnings() to see the first 50)\n",
      "R[write to console]: \n",
      "\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[Parallel(n_jobs=-1)]: Done 135 out of 140 | elapsed: 1868.5min remaining: 69.2min\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "R[write to console]: There were 50 or more warnings (use warnings() to see the first 50)\n",
      "R[write to console]: \n",
      "\n",
      "R[write to console]: There were 50 or more warnings (use warnings() to see the first 50)\n",
      "R[write to console]: \n",
      "\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[Parallel(n_jobs=-1)]: Done 137 out of 140 | elapsed: 2199.8min remaining: 48.2min\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "R[write to console]: There were 50 or more warnings (use warnings() to see the first 50)\n",
      "R[write to console]: \n",
      "\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[Parallel(n_jobs=-1)]: Done 140 out of 140 | elapsed: 2269.4min finished\n"
     ]
    }
   ],
   "source": [
    "outputs = Parallel(n_jobs=-1, verbose=100)(\n",
    "    [delayed(estimate_power)(sim, test) for sim in SIMULATIONS.keys() for test in TESTS]\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following code loops over each saved independence test file and generates absolute power curves for each test and for each simulation modality."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_power():\n",
    "    fig, ax = plt.subplots(nrows=4, ncols=5, figsize=(25,20))\n",
    "    \n",
    "    plt.suptitle(\"Multivariate Three-Sample Testing Increasing Angle\", y=0.93, va='baseline')\n",
    "    \n",
    "    for i, row in enumerate(ax):\n",
    "        for j, col in enumerate(row):\n",
    "            count = 5*i + j\n",
    "            sim = list(SIMULATIONS.keys())[count]\n",
    "            \n",
    "            for test in TESTS:\n",
    "                test_name = test.__name__\n",
    "                power = np.genfromtxt('../ksample/ksamp_vs_angle/{}_{}.csv'.format(sim, test_name),\n",
    "                                      delimiter=',')\n",
    "#                 manova_power = np.genfromtxt('../ksample/ksamp_vs_angle/{}_Manova.csv'.format(sim),\n",
    "#                                           delimiter=',')\n",
    "                \n",
    "                if test_name == \"MGC\":\n",
    "                    col.plot(ANGLES, power, color=\"#e41a1c\", label=test_name, lw=4)\n",
    "                elif test_name == \"KMERF\":\n",
    "                    col.plot(ANGLES, power, color=\"#e41a1c\", label=test_name, lw=4, linestyle='dashed')\n",
    "                elif test_name == \"Manova\":\n",
    "                    col.plot(ANGLES, power, color=\"#000000\", label=test_name, lw=4)\n",
    "                else:\n",
    "                    col.plot(ANGLES, power, label=test_name, lw=2)\n",
    "                col.set_xticks([])\n",
    "                if i == 3:\n",
    "                    col.set_xticks([ANGLES[0] - 1, ANGLES[-1] - 1])\n",
    "                col.set_ylim(-0.05, 1.05)\n",
    "                col.set_yticks([])\n",
    "                if j == 0:\n",
    "                    col.set_yticks([0, 1])\n",
    "                col.set_title(SIMULATIONS[sim])\n",
    "    \n",
    "    fig.text(0.5, 0.08, 'Angle', ha='center')\n",
    "    fig.text(0.08, 0.5, 'Statistical Power Relative to Manova', va='center', rotation='vertical')\n",
    "    leg = plt.legend(bbox_to_anchor=(0.5, 0.07), bbox_transform=plt.gcf().transFigure,\n",
    "                     ncol=len(TESTS), loc='upper center')\n",
    "    leg.get_frame().set_linewidth(0.0)\n",
    "    for legobj in leg.legendHandles:\n",
    "        legobj.set_linewidth(5.0)\n",
    "    plt.subplots_adjust(hspace=.50)\n",
    "    plt.savefig('../ksample/figs/ksamp_power_angle.pdf', transparent=True, bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_power()"
   ]
  },
  {
   "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.9.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}