{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import scipy\n",
    "from glob import glob\n",
    "from notebooks.lib import read_xyz_alt, read_load, linear, mean_across"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "outputs": [],
   "source": [
    "c_direct_neighbours = read_load(\"../data/c-3d-direct-neighbours\", read_xyz_alt)\n",
    "c_off_axis_neighbours = read_load(\"../data/c-3d-off-axis-neighbours\", read_xyz_alt)"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "outputs": [],
   "source": [
    "direct_meaned_by_N = mean_across(c_direct_neighbours)\n",
    "offaxis_meaned_by_N = mean_across(c_off_axis_neighbours)"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "outputs": [
    {
     "data": {
      "text/plain": "<Figure size 640x480 with 1 Axes>",
      "image/png": "[encoded data, 116160 characters]"
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": "<Figure size 640x480 with 1 Axes>",
      "image/png": "[encoded data, 144708 characters]"
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "for _, c in c_direct_neighbours.groupby('run'):\n",
    "    plt.plot(c.fd)\n",
    "\n",
    "for _, c in c_off_axis_neighbours.groupby('run'):\n",
    "    plt.plot(c.fd)\n",
    "\n",
    "plt.show()\n",
    "\n",
    "for _, run in c_direct_neighbours.groupby('run'):\n",
    "    plt.plot(run.N, run.cr)"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "outputs": [
    {
     "data": {
      "text/plain": "<Figure size 640x480 with 1 Axes>",
      "image/png": "[encoded data, 40160 characters]"
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": "(array([1.89170474, 0.58907793]),\n array([[ 5.09493737e-07, -1.86623148e-06],\n        [-1.86623148e-06,  6.96759885e-06]]))"
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "without_prefix = c_direct_neighbours\n",
    "\n",
    "p, pcov = scipy.optimize.curve_fit(linear, np.log(without_prefix.cr), np.log(without_prefix.N))\n",
    "linear_extent = np.linspace(0, np.max(np.log(c_direct_neighbours.cr)))\n",
    "\n",
    "plt.scatter(np.log(c_direct_neighbours.cr), np.log(c_direct_neighbours.N), s=1, marker='.', color=\"tab:blue\")\n",
    "plt.plot(linear_extent, linear(linear_extent, *p), color=\"tab:red\")\n",
    "\n",
    "plt.xlabel(\"$\\\\log r_{max}$\")\n",
    "plt.ylabel(\"$\\\\log N$\")\n",
    "plt.show()\n",
    "p, pcov"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/var/folders/yz/hjlw4lyn3rq0s5y_r82b55cw0000gn/T/ipykernel_38665/3346041710.py:1: FutureWarning: The default value of numeric_only in DataFrameGroupBy.mean is deprecated. In a future version, numeric_only will default to False. Either specify numeric_only or select only columns which should be valid for the function.\n",
      "  mn = c_direct_neighbours.groupby(\"N\").mean().reset_index()\n"
     ]
    },
    {
     "data": {
      "text/plain": "<Figure size 640x480 with 1 Axes>",
      "image/png": "[encoded data, 24380 characters]"
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": "(array([1.89170474, 0.58907793]),\n array([[ 5.09493737e-07, -1.86623148e-06],\n        [-1.86623148e-06,  6.96759885e-06]]))"
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mn = c_direct_neighbours.groupby(\"N\").mean().reset_index()\n",
    "mn_p, mn_pcov = scipy.optimize.curve_fit(linear, np.log(mn.cr), np.log(mn.N))\n",
    "linear_extent = np.linspace(0, np.max(np.log(mn.cr)))\n",
    "\n",
    "plt.scatter(np.log(mn.cr), np.log(mn.N), s=1, marker='.', color=\"tab:blue\")\n",
    "plt.plot(linear_extent, linear(linear_extent, *mn_p), color=\"tab:red\")\n",
    "\n",
    "plt.xlabel(\"$\\\\log r_{max}$\")\n",
    "plt.ylabel(\"$\\\\log N$\")\n",
    "plt.show()\n",
    "p, pcov"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "outputs": [],
   "source": [
    "rust_direct = read_load(\"../data/rust-3d\", read_xyz_alt)\n",
    "rust_direct_meaned_by_N = mean_across(rust_direct)"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "outputs": [],
   "source": [
    "rust_offaxis = read_load(\"../data/rust-3d-offaxis\", read_xyz_alt)\n",
    "rust_offaxis_meaned_by_N = mean_across(rust_direct)"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "outputs": [
    {
     "data": {
      "text/plain": "<Figure size 600x600 with 1 Axes>",
      "image/png": "[encoded data, 57200 characters]"
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(figsize=(6, 6))\n",
    "\n",
    "\n",
    "def nc_fd(df, label, color=None, ignore_prefix=0):\n",
    "    plt.fill_between(\n",
    "        df.N,\n",
    "        # TODO Check error math here\n",
    "        (df[ignore_prefix:]['fd']['mean'] - df[ignore_prefix:]['fd']['stderr']),\n",
    "        (df[ignore_prefix:]['fd']['mean'] + df[ignore_prefix:]['fd']['stderr']),\n",
    "        color=color,\n",
    "        alpha=0.2,\n",
    "        label=f\"{label}, standard error band\"\n",
    "    )\n",
    "\n",
    "    plt.plot(\n",
    "        df.N,\n",
    "        df['fd']['mean'],\n",
    "        color=color,\n",
    "        label=f\"{label}, fd mean\"\n",
    "    )\n",
    "\n",
    "\n",
    "nc_fd(direct_meaned_by_N, label=\"IPC + 3D Direct\", color=\"tab:blue\")\n",
    "# nc_fd(offaxis_meaned_by_N, label=\"IPC + 3D Off-axis\", color=\"tab:green\")\n",
    "nc_fd(rust_direct_meaned_by_N, label=\"NF 3D Direct\", color=\"tab:orange\")\n",
    "nc_fd(rust_offaxis_meaned_by_N, label=\"NF 3D Off-axis\", color=\"tab:purple\")\n",
    "\n",
    "plt.plot([np.min(without_prefix.N), np.max(without_prefix.N)], [2.5, 2.5], color='tab:red', label='Theory')\n",
    "plt.fill_between(without_prefix.N, 2.5 - 0.01, 2.5 + 0.01, alpha=0.2, color='tab:red', label='Theory error band')\n",
    "\n",
    "plt.xlabel(\"$N_C$\")\n",
    "plt.ylabel(\"$fd$ (instantaneous)\")\n",
    "plt.legend()\n",
    "\n",
    "plt.savefig('../figures/3d-nc-fd-convergence.svg')\n",
    "plt.savefig('../figures/3d-nc-fd-convergence.png')"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "outputs": [],
   "source": [
    "from notebooks.graphs import cr_n"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "outputs": [
    {
     "data": {
      "text/plain": "(array([1.89170474, 0.58907793]),\n array([[ 5.09493737e-07, -1.86623148e-06],\n        [-1.86623148e-06,  6.96759885e-06]]))"
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/plain": "<Figure size 640x480 with 1 Axes>",
      "image/png": "[encoded data, 40160 characters]"
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "cr_n(c_direct_neighbours)"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "outputs": [
    {
     "data": {
      "text/plain": "(array([1.90130968, 0.47834345]),\n array([[ 5.61079868e-07, -2.07748262e-06],\n        [-2.07748262e-06,  7.83527791e-06]]))"
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/plain": "<Figure size 640x480 with 1 Axes>",
      "image/png": "[encoded data, 40180 characters]"
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "cr_n(rust_direct)"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "outputs": [
    {
     "data": {
      "text/plain": "(array([1.89843243, 0.71870004]),\n array([[ 5.17525369e-07, -1.85359693e-06],\n        [-1.85359693e-06,  6.77175199e-06]]))"
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/plain": "<Figure size 640x480 with 1 Axes>",
      "image/png": "[encoded data, 38660 characters]"
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "cr_n(rust_offaxis)"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [],
   "metadata": {
    "collapsed": false
   }
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}
