Analysis.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%reload_ext autoreload\n",
    "%autoreload 2\n",
    "\n",
    "import matplotlib.pylab as plt\n",
    "import numpy as np\n",
    "import seaborn as sns\n",
    "\n",
    "from angle_set import AngleSet, combine_fullrank_matrix\n",
    "from helpers import savefig"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Constraints analysis\n",
    "\n",
    "Visualize the constraints matrices and compare number of obtained constraints with degrees of freedom.\n",
    "\n",
    "To reproduce Figure 2 of the paper, you can jump to the end of the file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "d = 2  # do not change.\n",
    "N = 5\n",
    "\n",
    "np.random.seed(51)\n",
    "\n",
    "angle_set = AngleSet(N=N, d=d)\n",
    "angle_set.set_points(mode='random')\n",
    "angle_set.plot_all()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Apoly, bpoly = angle_set.get_polygon_constraints([3, 4], print_out=True)\n",
    "assert np.allclose(Apoly.dot(angle_set.theta), bpoly)\n",
    "plt.matshow(np.c_[Apoly, bpoly.reshape((-1, 1))])\n",
    "plt.title('polygon constraints')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Arays, brays = angle_set.get_ray_constraints(verbose=False)\n",
    "assert np.allclose(Arays.dot(angle_set.theta), brays)\n",
    "plt.matshow(np.c_[Arays, brays.reshape((-1, 1))])\n",
    "plt.title('single constraints')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Afull, bfull = combine_fullrank_matrix(Arays, Apoly, brays, bpoly, False)\n",
    "print('kept rows: {} of {}'.format(Afull.shape[0] - Arays.shape[0], Apoly.shape[0]))\n",
    "plt.matshow(np.c_[Afull, bfull])\n",
    "plt.title('full rank matrix')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_rays, n_poly = angle_set.get_n_rays(), angle_set.get_n_poly()\n",
    "\n",
    "assert n_rays == Arays.shape[0]\n",
    "print('number of ray constraints found:', n_rays)\n",
    "\n",
    "assert n_poly == Afull.shape[0] - Arays.shape[0]\n",
    "print('number of poly constraints found:', n_poly)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Create plot for paper\n",
    "\n",
    "### Generate data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "d = 2  # do not change.\n",
    "\n",
    "print('N \\t M \\t DOF \\t necessary \\t linear \\t nonlinear')\n",
    "\n",
    "Ns = np.arange(4, 8) # to make it notebook run faster. \n",
    "#Ns = np.arange(4, 20) # uncomment for paper results\n",
    "\n",
    "Ms = np.empty((len(Ns), ))\n",
    "DOFs = np.empty((len(Ns), ))\n",
    "n_linear_list = np.empty((len(Ns), ))\n",
    "n_nonlinear_list = np.empty((len(Ns), ))\n",
    "n_poly_list = np.empty((len(Ns), ))\n",
    "n_rays_list = np.empty((len(Ns), ))\n",
    "\n",
    "for i, N in enumerate(Ns):\n",
    "    angle_set = AngleSet(N=N, d=d)\n",
    "    angle_set.set_points(mode='random')\n",
    "    n_rays = angle_set.get_n_rays()\n",
    "    n_poly = angle_set.get_n_poly()\n",
    "\n",
    "    Apoly, bpoly = angle_set.get_polygon_constraints([3])\n",
    "    Aray, bray = angle_set.get_ray_constraints()\n",
    "\n",
    "    # Sanity check only. is taking too long for N >= 11\n",
    "    if N < 11:\n",
    "        Afull, bfull = combine_fullrank_matrix(Aray, Apoly, bray.flatten(), bpoly, False)\n",
    "        assert n_rays == Aray.shape[0]\n",
    "        assert n_poly == Afull.shape[0] - Aray.shape[0]\n",
    "\n",
    "    Afull = np.vstack([Aray, Apoly[:n_poly]])\n",
    "    bfull = np.hstack([bray, bpoly[:n_poly]])\n",
    "\n",
    "    #DOF = angle_set.N * angle_set.d - 2*angle_set.d - 1\n",
    "    DOF = angle_set.get_DOF()\n",
    "    necessary = angle_set.num_angles - DOF\n",
    "\n",
    "    nonlinear = sum([1] + [i for i in range(2, angle_set.N - 2)])\n",
    "    assert nonlinear == angle_set.get_n_sine()\n",
    "\n",
    "    assert nonlinear + n_rays + n_poly == necessary\n",
    "    print('{} \\t {} \\t {} \\t {} \\t {} \\t {}'.format(angle_set.N, angle_set.num_angles, DOF, necessary, n_rays + n_poly,\n",
    "                                                    nonlinear))\n",
    "\n",
    "    n_poly_list[i] = n_poly\n",
    "    n_rays_list[i] = n_rays\n",
    "    n_linear_list[i] = n_rays + n_poly\n",
    "    n_nonlinear_list[i] = nonlinear\n",
    "    Ms[i] = angle_set.num_angles\n",
    "    DOFs[i] = DOF"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Create plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_here(ax, range_):\n",
    "    colors = sns.color_palette('tab20')\n",
    "    ax.stackplot(Ns[range_], [n_rays_list[range_], n_poly_list[range_], \n",
    "                              n_nonlinear_list[range_], DOFs[range_]],\n",
    "                 labels=['single', 'triangle', 'non-linear', 'DOF'],\n",
    "                 colors=colors)\n",
    "    ax.plot(Ns[range_], Ms[range_], color='black', label='no. angles $M$', ls=':')\n",
    "    ax.set_xlim(Ns[range_][0], Ns[range_][-1])\n",
    "\n",
    "\n",
    "def plot_inset(ax, loc=[0.5, 0.5, 0.47, 0.47], \n",
    "               range_full=[0, 10], range_inset=[14, 15, 1000, 1500]):\n",
    "    \"\"\" Create zoom inside the plot. \"\"\"\n",
    "    axins = ax.inset_axes(loc)  # location (x,y), size\n",
    "    plot_here(axins, range_full)\n",
    "    # sub region of the original image\n",
    "    x1, x2, y1, y2 = range_inset\n",
    "    axins.set_xlim(x1, x2)\n",
    "    axins.set_ylim(y1, y2)\n",
    "    axins.set_xticklabels('')\n",
    "    axins.set_yticklabels('')\n",
    "    axins.set_xticks([])\n",
    "    axins.set_yticks([])\n",
    "    axins.set_facecolor((1, 1, 1, 0.5))\n",
    "    r_, lines_ = ax.indicate_inset_zoom(axins, label=None)\n",
    "    [l.set_color('black') for l in lines_]\n",
    "    [l.set_alpha(1.0) for l in lines_]\n",
    "    r_.set_alpha(1.0)\n",
    "    r_.set_edgecolor('black')\n",
    "\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "fig.set_size_inches(6, 2.5)\n",
    "\n",
    "range_ = np.arange(len(Ns))\n",
    "#range_ = np.arange(14) # uncomment for paper results\n",
    "\n",
    "plot_here(ax, range_)\n",
    "plot_inset(ax, loc=[0.02, 0.25, 0.5, 0.35], range_full=range_, range_inset=[4, 8, 0, 200])\n",
    "h, l = ax.get_legend_handles_labels()\n",
    "ax.legend(h[::-1], l[::-1], loc='upper left', ncol=3)\n",
    "ax.set_ylabel('count', fontsize=12)\n",
    "ax.set_xlabel('number of points $N$', fontsize=12)\n",
    "plt.tight_layout()\n",
    "ax.grid()\n",
    "\n",
    "#savefig(fname='plots/number_constraints.pdf')"
   ]
  },
  {
   "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.6.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}