std_dev_vs_num_runs.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "# Imports\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import os\n",
    "from datetime import datetime\n",
    "from matplotlib.ticker import MultipleLocator\n",
    "from scipy.stats import t as t_dist\n",
    "from code.Utils import ESFitness\n",
    "\n",
    "# Default variable definitions\n",
    "num_samples = 100\n",
    "available_files = [ \n",
    "    # (2, 3), \n",
    "    (5, 1),  (5, 3),  (5, 6),  (5, 9),  (5, 10),  (5, 14),  (5, 17),  (5, 20),  (5, 23),\n",
    "    (10, 1), (10, 3), (10, 6), (10, 9), (10, 10), (10, 14), (10, 17), (10, 20), (10, 23),\n",
    "    (20, 1), (20, 3), (20, 6), (20, 9), (20, 10), (20, 14), (20, 17), (20, 20), (20, 23),\n",
    "]\n",
    "# available_files = [ \n",
    "#     (5, 3),  (5, 9),  (5, 14),  (5, 17),  (5, 23),\n",
    "#     (10, 3), (10, 9), (10, 14), (10, 17), (10, 23),\n",
    "#     (20, 3), (20, 9), (20, 14), (20, 17), (20, 23),\n",
    "# ]\n",
    "\n",
    "folder_name = 'C:\\\\src\\\\master-thesis\\\\experiments\\\\num_runs_vs_std_dev'\n",
    "data_file_name = 'raw_data\\\\GA_results_{ndim}dim_f{fid}.tdat'\n",
    "save_file_name = 'processed_data\\\\samples_data_{ndim}dim_f{fid}.npz'\n",
    "other_save_file_name = 'processed_data\\\\normalized_means_and_spread.npz'\n",
    "distances_save_file = 'processed_data\\\\distances.npz'\n",
    "plot_file_prefix = 'plots\\\\'\n",
    "\n",
    "plot_file_extension = '.pdf'\n",
    "fig_size = (8,6)\n",
    "\n",
    "os.chdir(folder_name)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "samples_data_5dim_f1.npz already exists, skipping...\nsamples_data_5dim_f3.npz already exists, skipping...\nsamples_data_5dim_f6.npz already exists, skipping...\nsamples_data_5dim_f9.npz already exists, skipping...\nsamples_data_5dim_f10.npz already exists, skipping...\nsamples_data_5dim_f14.npz already exists, skipping...\nsamples_data_5dim_f17.npz already exists, skipping...\nsamples_data_5dim_f20.npz already exists, skipping...\nsamples_data_5dim_f23.npz already exists, skipping...\nsamples_data_10dim_f1.npz already exists, skipping...\nsamples_data_10dim_f3.npz already exists, skipping...\nsamples_data_10dim_f6.npz already exists, skipping...\nsamples_data_10dim_f9.npz already exists, skipping...\nsamples_data_10dim_f10.npz already exists, skipping...\nsamples_data_10dim_f14.npz already exists, skipping...\nsamples_data_10dim_f17.npz already exists, skipping...\nsamples_data_10dim_f20.npz already exists, skipping...\nsamples_data_10dim_f23.npz already exists, skipping...\nsamples_data_20dim_f1.npz already exists, skipping...\nsamples_data_20dim_f3.npz already exists, skipping...\nsamples_data_20dim_f6.npz already exists, skipping...\nsamples_data_20dim_f9.npz already exists, skipping...\nsamples_data_20dim_f10.npz already exists, skipping...\nsamples_data_20dim_f14.npz already exists, skipping...\nsamples_data_20dim_f17.npz already exists, skipping...\nsamples_data_20dim_f20.npz already exists, skipping...\nsamples_data_20dim_f23.npz already exists, skipping...\n0:00:00.032000\n"
     ]
    }
   ],
   "source": [
    "objects = {}\n",
    "\n",
    "start = datetime.now()\n",
    "\n",
    "# Create ESFitness objects from data files\n",
    "for ndim, fid in available_files:\n",
    "    if save_file_name.format(ndim=ndim, fid=fid)[15:] in os.listdir('processed_data'):\n",
    "        print(\"{} already exists, skipping...\".format(save_file_name.format(ndim=ndim, fid=fid)[15:]))\n",
    "        continue\n",
    "    with open(data_file_name.format(ndim=ndim, fid=fid), 'r') as f:\n",
    "        lines = [line for line in f]\n",
    "    objects[(ndim, fid)] = [eval(line) for line in lines]\n",
    "\n",
    "    num_ESs = len(objects[(ndim, fid)])\n",
    "    num_runs = len(objects[(ndim, fid)][0].min_fitnesses)\n",
    "    \n",
    "    means = np.zeros((num_runs, num_ESs, num_samples))\n",
    "    medians = np.zeros((num_runs, num_ESs, num_samples))\n",
    "    std_devs = np.zeros((num_runs, num_ESs, num_samples))\n",
    "    \n",
    "    for obj in objects[(ndim, fid)]:\n",
    "        obj.min_fitnesses = np.array(obj.min_fitnesses)\n",
    "    \n",
    "    for sample_size in range(2, num_runs):\n",
    "        samples = np.zeros((num_ESs, num_samples, sample_size))\n",
    "        for sample_num in range(num_samples):\n",
    "            sample_indices = np.random.choice(num_runs, sample_size, replace=False)\n",
    "            for ES_num in range(num_ESs):\n",
    "                obj = objects[(ndim, fid)][ES_num]\n",
    "                samples[ES_num,sample_num,:] = obj.min_fitnesses[sample_indices]\n",
    "    \n",
    "        means[sample_size, :, :] = np.mean(samples, axis=2)\n",
    "        # medians[sample_size, :, :] = np.median(samples, axis=2)\n",
    "        std_devs[sample_size, :, :] = np.std(samples, axis=2)\n",
    "\n",
    "    save_file = save_file_name.format(ndim=ndim, fid=fid)\n",
    "    # np.savez(save_file, means=means, medians=medians, std_devs=std_devs)\n",
    "    np.savez(save_file, means=means, std_devs=std_devs)\n",
    "\n",
    "stop = datetime.now()\n",
    "print(stop-start)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [],
   "source": [
    "for ndim, fid in available_files:\n",
    "    \n",
    "    with open(save_file_name.format(ndim=ndim, fid=fid), 'rb') as save_file:\n",
    "        npzfile = np.load(save_file)\n",
    "        means = npzfile['means']\n",
    "        num_runs, num_ESs, num_samples = means.shape\n",
    "    \n",
    "        ESs_to_be_plotted = range(0, num_ESs, 30)\n",
    "        plot_length = num_runs//4\n",
    "            \n",
    "        for ES in ESs_to_be_plotted:\n",
    "            x_data = range(2, plot_length)\n",
    "            y_data = np.mean(means[2:plot_length,ES,:], axis=1)\n",
    "            y_error = np.std(means[2:plot_length,ES,:], axis=1)\n",
    "            data_mean = np.mean(y_data)\n",
    "            y_data = y_data / data_mean\n",
    "            y_error = y_error / data_mean\n",
    "            \n",
    "            plot_file_name = \"mean_std_dev_errorbar_{}dim_f{}_ES{}_normalized\".format(ndim, fid, ES)\n",
    "            \n",
    "            plt.figure(figsize=fig_size)\n",
    "            plt.axhline(y=1, color='k')\n",
    "            plt.errorbar(x=x_data, y=y_data, yerr=y_error, linestyle='None', marker='o')\n",
    "            plt.title(\"Normalized Mean and Standard Deviation for ES {}/{} in {}dim F{}\".format(ES, num_ESs, \n",
    "                                                                                                ndim, fid))\n",
    "            plt.xlabel(\"Number of runs\")\n",
    "            plt.tight_layout()\n",
    "            plt.savefig(plot_file_prefix + plot_file_name + plot_file_extension)\n",
    "            plt.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "all_data = []\n",
    "all_errors = []\n",
    "total_ESs = 0\n",
    "\n",
    "for ndim, fid in available_files:\n",
    "    \n",
    "    with open(save_file_name.format(ndim=ndim, fid=fid), 'rb') as save_file:\n",
    "        npzfile = np.load(save_file)\n",
    "        means = npzfile['means']\n",
    "        num_runs, num_ESs, num_samples = means.shape\n",
    "        plot_length = num_runs\n",
    "        total_ESs += num_ESs\n",
    "            \n",
    "        for ES in range(num_ESs):\n",
    "            y_data = np.mean(means[2:plot_length,ES,:], axis=1)\n",
    "            y_error = np.std(means[2:plot_length,ES,:], axis=1)\n",
    "            data_mean = np.mean(y_data)\n",
    "            y_data = (y_data / data_mean)\n",
    "            y_error = y_error / data_mean\n",
    "            \n",
    "            all_data.append(y_data)\n",
    "            all_errors.append(y_error)\n",
    "\n",
    "all_data = np.mean(np.array(all_data), axis=0)\n",
    "all_errors = np.mean(np.array(all_errors), axis=0)\n",
    "x_data = range(2, plot_length)\n",
    "np.savez(other_save_file_name, means=all_data, std_devs=all_errors, x_range=x_data)\n",
    "\n",
    "plot_file_name = \"std_err_plot_normalized_aggregated\"\n",
    "\n",
    "_, ax = plt.subplots(figsize=fig_size)\n",
    "plt.grid(True, which='both', axis='y', linestyle=':', color='k', alpha=0.75)\n",
    "plt.plot(x_data, all_errors, 'b-')\n",
    "plt.title(\"Standard Error, Aggregated over {} different ESs\".format(total_ESs))\n",
    "plt.xlabel(\"Number of runs\")\n",
    "plt.ylabel(\"Relative standard error\")\n",
    "plt.xlim(xmax=x_data[-1])\n",
    "plt.minorticks_on()\n",
    "plt.tight_layout()\n",
    "plt.savefig(plot_file_prefix + plot_file_name + plot_file_extension)\n",
    "plt.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "271694\n[0.0, 0.029405521658226862, 0.073823005157169691, 0.13007347685890025, 0.20347927289976137, 0.3005954465849387, 0.43402301913608449, 0.63528328460738781, 0.95456442504846895, 1.4722950644871693, 2.3872675169699238, 4.2920310674249098, 8.6178515443898984, 22.392458282697451, 73.941823770233682, 416.35203358269729, 5169.766250033229, 140494.2572353656, 9165357.8847258277, 457467627.84533823, 51590306299660.875]\n"
     ]
    }
   ],
   "source": [
    "def plot_cdf(data, num_cases, num_bins=200, ndim=None, fid=None):\n",
    "    start_index = 0\n",
    "    while data[start_index] == 0:\n",
    "        start_index += 1\n",
    "    \n",
    "    plot_data = data[start_index:]\n",
    "    # Create bins so each contains an equal number of points, but remove 0-bins if many values are the same \n",
    "    bins = sorted(set(plot_data[::len(plot_data)//num_bins]))\n",
    "    if plot_data[-1] > bins[-1]:\n",
    "        bins.append(plot_data[-1])\n",
    "    \n",
    "    if ndim is not None and fid is not None:\n",
    "        title = \"Distance histogram for {} distances between {} ESs\".format(len(data), num_cases)\n",
    "        title += \" in {}dim F{}\".format(ndim, fid)\n",
    "        plot_file_name = \"distances_hist_{}dim_F{}\".format(ndim, fid)\n",
    "    else:\n",
    "        title = \"Aggregate distance histogram for {} distances between {} ESs\".format(len(data), num_cases)\n",
    "        plot_file_name = \"distances_hist_aggregate\"\n",
    "    plt.figure(figsize=fig_size)\n",
    "    plt.title(title)\n",
    "    _, _, patches = plt.hist(plot_data, bins=bins, cumulative=True, histtype='step', normed=1)\n",
    "    patches[0].set_xy(patches[0].get_xy()[:-1])  # Remove the downward line at the end \n",
    "    plt.grid(True, axis='both', linestyle=':', color='k', alpha=0.75)\n",
    "    plt.xscale('log')\n",
    "    plt.ylim(ymax=1)\n",
    "    plt.yticks(np.arange(0, 1.1, 0.1))\n",
    "    plt.tight_layout()\n",
    "    plt.savefig(plot_file_prefix + plot_file_name + plot_file_extension)\n",
    "    plt.close()\n",
    "\n",
    "all_distances = []\n",
    "total_ESs = 0\n",
    "\n",
    "for ndim, fid in available_files:\n",
    "    with open(data_file_name.format(ndim=ndim, fid=fid), 'r') as f:\n",
    "        lines = [line for line in f]\n",
    "    objects = [eval(line) for line in lines]\n",
    "    FCEs = [obj.FCE for obj in objects]\n",
    "    ERTs = [obj.ERT for obj in objects]\n",
    "    \n",
    "    num_ESs = len(FCEs)\n",
    "    total_ESs += num_ESs\n",
    "    distances = []\n",
    "    for ES in range(num_ESs-1):\n",
    "        for other_ES in range(ES+1, num_ESs):\n",
    "            this = FCEs[ES]\n",
    "            other = FCEs[other_ES]\n",
    "            dist = np.abs(this - other)\n",
    "            if dist != 0:\n",
    "                distances.append(dist / min(this, other))\n",
    "            else:\n",
    "                this = ERTs[ES]\n",
    "                other = ERTs[other_ES]\n",
    "                if this is None or other is None:\n",
    "                    distances.append(0)  # Apparently, distance here is actually 0... :/\n",
    "                else:\n",
    "                    distances.append(np.abs(this - other) / min(this, other))\n",
    "    \n",
    "    distances.sort()\n",
    "    # plot_cdf(distances, num_ESs, ndim=ndim, fid=fid)\n",
    "    \n",
    "    all_distances.extend(distances)\n",
    "\n",
    "all_distances.sort()\n",
    "np.savez(distances_save_file, distances=all_distances)\n",
    "print(len(all_distances))\n",
    "print(all_distances[::len(all_distances)//20])\n",
    "\n",
    "plot_cdf(all_distances, total_ESs, num_bins=1000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "distances = np.load(distances_save_file)['distances']\n",
    "means = np.load(other_save_file_name)['means']\n",
    "std_devs = np.load(other_save_file_name)['std_devs']\n",
    "x_range = np.load(other_save_file_name)['x_range']\n",
    "five_percent_points = distances[len(distances)//20::len(distances)//20]\n",
    "\n",
    "max_num_runs = 128\n",
    "probabilities = {}\n",
    "mean_A = 1\n",
    "for i in range(20):\n",
    "    probabilities[i] = []\n",
    "    mean_B = 1 + five_percent_points[i]\n",
    "    \n",
    "    for j, n in enumerate(x_range):\n",
    "        std_dev_A = std_devs[j]\n",
    "        std_dev_B = std_dev_A * mean_B\n",
    "        std_error = np.sqrt((std_dev_A**2 + std_dev_B**2) / n)\n",
    "        t = (mean_B - mean_A) / std_error\n",
    "        \n",
    "        df = 2*n - 2\n",
    "        p_value = 2 * t_dist.sf(t, df=df)\n",
    "        probabilities[i].append(p_value)\n",
    "\n",
    "markers = ['bs', 'bo', 'b^', 'bv', 'b*', \n",
    "           'gs', 'go', 'g^', 'gv', 'g*', \n",
    "           'rs', 'ro', 'r^', 'rv', 'r*', \n",
    "           'ys', 'yo', 'y^', 'yv', 'y*', ]\n",
    "\n",
    "_, ax = plt.subplots(figsize=fig_size)\n",
    "ax.yaxis.set_major_locator(MultipleLocator(0.2))\n",
    "ax.yaxis.set_minor_locator(MultipleLocator(0.05))\n",
    "ax.xaxis.set_major_locator(MultipleLocator(20))\n",
    "ax.xaxis.set_minor_locator(MultipleLocator(10))\n",
    "plt.title('P-value of ES comparison at relative distance VS number of runs')\n",
    "plt.xlabel('number of runs')\n",
    "plt.ylabel('p-value')\n",
    "\n",
    "for i in range(10):\n",
    "    p = (i+1) * 5\n",
    "    dist = five_percent_points[i]\n",
    "    plt.plot(x_range, probabilities[i], markers[i], label='{}th %-ile: {:.2}'.format(p, dist))\n",
    "\n",
    "plt.xlim(xmax=max_num_runs)\n",
    "plt.minorticks_on()\n",
    "plt.grid(True, axis='both', which='both', linestyle=':', color='k', alpha=0.75)    \n",
    "plt.legend(numpoints=1)\n",
    "\n",
    "plot_file_name = 'certainty_at_distance_vs_num_runs'\n",
    "plt.tight_layout()\n",
    "plt.savefig(plot_file_prefix + plot_file_name + plot_file_extension)\n",
    "plt.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_histogram(min_fitnesses, nbins=20, ndim=None, fid=None):\n",
    "    if ndim is not None and fid is not None:\n",
    "        plot_file_name = \"fitness_hist_{}dim_F{}\".format(ndim, fid)\n",
    "        title = \"Distribution of normalized fitness minimum values for {}dim F{}\".format(ndim, fid)\n",
    "    else:\n",
    "        plot_file_name = \"fitness_hist_aggregate\"\n",
    "        title = \"Aggregate distribution of {} normalized fitness minimum values\".format(len(min_fitnesses))\n",
    "\n",
    "    plt.figure(figsize=fig_size)\n",
    "    plt.title(title)\n",
    "    plt.hist(min_fitnesses, bins=nbins, histtype='step')\n",
    "    plt.yscale('log')\n",
    "    plt.tight_layout()\n",
    "    plt.savefig(plot_file_prefix + plot_file_name + plot_file_extension)\n",
    "    plt.close()\n",
    "\n",
    "all_FCEs = []\n",
    "all_ERTs = []\n",
    "all_min_fitnesses = []\n",
    "\n",
    "for ndim, fid in available_files:\n",
    "    with open(data_file_name.format(ndim=ndim, fid=fid), 'r') as f:\n",
    "        min_fitnesses = []\n",
    "        for line in f:\n",
    "            obj = eval(line)\n",
    "            min_fits = np.array(obj.min_fitnesses)\n",
    "            min_fits = (min_fits / np.mean(obj.min_fitnesses)).tolist()\n",
    "            min_fitnesses.extend(min_fits)\n",
    "\n",
    "        all_min_fitnesses.extend(min_fitnesses)\n",
    "\n",
    "plot_histogram(all_min_fitnesses, nbins=25)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    ""
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2.0
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}