cd_tests.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import os \n",
    "os.environ[\"ISISROOT\"] = \"/usgs/cpkgs/anaconda3_linux/envs/isis3.9.0\"\n",
    "\n",
    "from pysis import isis\n",
    "from pysis.exceptions import ProcessError\n",
    "\n",
    "import matplotlib\n",
    "from matplotlib import pyplot as plt\n",
    "\n",
    "from plio.io.io_gdal import GeoDataset \n",
    "import numpy as np\n",
    "\n",
    "from autocnet.examples import get_path\n",
    "from autocnet.cg import change_detection as cd\n",
    "from autocnet.utils.utils import bytescale\n",
    "\n",
    "import matplotlib"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def cartesian_product(*arrays):\n",
    "    la = len(arrays)\n",
    "    dtype = \"object\"\n",
    "    arr = np.empty([len(a) for a in arrays] + [la], dtype=dtype)\n",
    "    for i, a in enumerate(np.ix_(*arrays)):\n",
    "        arr[...,i] = a\n",
    "    return arr.reshape(-1, la)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Set up images for now change \n",
    "nochange_before_geo = GeoDataset(get_path(\"no_change.before.cub\"))\n",
    "nochange_after_geo = GeoDataset(get_path(\"no_change.after.cub\"))\n",
    "\n",
    "nochange_before_arr = nochange_before_geo.read_array()\n",
    "nochange_after_arr = nochange_after_geo.read_array()\n",
    "nochange_before_arr[nochange_before_arr==nochange_before_arr.min()] = np.nan\n",
    "nochange_after_arr[nochange_after_arr==nochange_after_arr.min()] = np.nan\n",
    "\n",
    "\n",
    "boulder_shift_before_geo = GeoDataset(get_path(\"boulder_shift.before.cub\"))\n",
    "boulder_shift_after_geo = GeoDataset(get_path(\"boulder_shift.after.cub\"))\n",
    "\n",
    "boulder_shift_before_arr = boulder_shift_before_geo.read_array()\n",
    "boulder_shift_after_arr = boulder_shift_after_geo.read_array()\n",
    "boulder_shift_before_arr[boulder_shift_before_arr==boulder_shift_before_arr.min()] = np.nan\n",
    "boulder_shift_after_arr[boulder_shift_after_arr==boulder_shift_after_arr.min()] = np.nan"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "9X3 Matrix\n"
     ]
    }
   ],
   "source": [
    "#{\"nfeatures\": 2000, \"scaleFactor\": 1.1, \"nlevels\": 1}\n",
    "# cluster_params={\"min_samples\": 10, \"max_eps\": 10, \"eps\": None, \"xi\":.5}\n",
    "\n",
    "nfeatures = [500, 1000, 2000]\n",
    "scalefactor = [1.1]\n",
    "nlevels = [1] \n",
    "extractor_params = [dict(zip([\"nfeatures\", \"scaleFactor\", \"nlevels\"], ps)) for ps in cartesian_product(nfeatures, scalefactor, nlevels)] \n",
    "\n",
    "\n",
    "min_samples = [10, 15, 20]\n",
    "max_eps = [10, 20, 30]\n",
    "xi = [.5]\n",
    "eps = [.5]\n",
    "cluster_params =  [dict(zip([\"min_samples\", \"max_eps\", \"eps\", \"xi\"], ps)) for ps in cartesian_product(min_samples, max_eps, eps, xi)] \n",
    "\n",
    "print(f\"{len(cluster_params)}X{len(extractor_params)} Matrix\")\n",
    "parameters = cartesian_product(extractor_params, cluster_params)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/krodriguez/.conda/envs/autocnetdev/lib/python3.6/site-packages/sklearn/cluster/_optics.py:502: UserWarning: All reachability values are inf. Set a larger max_eps or all data will be considered outliers.\n",
      "  UserWarning)\n",
      "/home/krodriguez/.conda/envs/autocnetdev/lib/python3.6/site-packages/sklearn/cluster/_optics.py:502: UserWarning: All reachability values are inf. Set a larger max_eps or all data will be considered outliers.\n",
      "  UserWarning)\n",
      "/home/krodriguez/.conda/envs/autocnetdev/lib/python3.6/site-packages/sklearn/cluster/_optics.py:502: UserWarning: All reachability values are inf. Set a larger max_eps or all data will be considered outliers.\n",
      "  UserWarning)\n"
     ]
    }
   ],
   "source": [
    "boulder_shift_results = [cd.okubogar_detector(boulder_shift_before_geo, boulder_shift_after_geo, extractor_kwargs=params[0], cluster_params=params[1]) for params in parameters]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 180,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "nx = len(extractor_params)\n",
    "ny = len(cluster_params)\n",
    "\n",
    "fig, axs = plt.subplots(ny, nx, sharex='col', sharey='row')\n",
    "fig.set_figheight(50)\n",
    "fig.set_figwidth(20)\n",
    "for ax in axs.flat:\n",
    "    ax.set(xlabel='x-label', ylabel='y-label')\n",
    "    \n",
    "xs, ys = np.mgrid[:ny, :nx]\n",
    "cmap = matplotlib.cm.get_cmap('plasma')\n",
    "\n",
    "for i, coord in enumerate(zip(xs.flatten(),ys.flatten())):\n",
    "    x,y = coord\n",
    "    polys, weights, diff = boulder_shift_results[i]\n",
    "    \n",
    "    indices = [i for (i,w) in enumerate(weights) if w > 0]\n",
    "    filtered_polys = []\n",
    "    filtered_polys = [polys[i] for i in indices]\n",
    "    \n",
    "    axs[x,y].set(xlabel=f\"ex {parameters[i][1]}\", ylabel=f\"cl {parameters[i][0]}\")\n",
    "    axs[x,y].imshow(diff, cmap=\"Greys_r\")\n",
    "    axs[x,y].imshow(boulder_shift_after_arr, cmap=\"Greys_r\", alpha=.5)\n",
    "\n",
    "    if weights:\n",
    "        colors = bytescale(np.asarray(weights))\n",
    "        ax = [axs[x,y].fill(*p.exterior.xy, c=\"r\", alpha=.6) for i,p in enumerate(filtered_polys)]\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/krodriguez/.conda/envs/autocnetdev/lib/python3.6/site-packages/sklearn/cluster/_optics.py:502: UserWarning: All reachability values are inf. Set a larger max_eps or all data will be considered outliers.\n",
      "  UserWarning)\n",
      "/home/krodriguez/.conda/envs/autocnetdev/lib/python3.6/site-packages/sklearn/cluster/_optics.py:502: UserWarning: All reachability values are inf. Set a larger max_eps or all data will be considered outliers.\n",
      "  UserWarning)\n",
      "/home/krodriguez/.conda/envs/autocnetdev/lib/python3.6/site-packages/sklearn/cluster/_optics.py:502: UserWarning: All reachability values are inf. Set a larger max_eps or all data will be considered outliers.\n",
      "  UserWarning)\n"
     ]
    }
   ],
   "source": [
    "no_change_results = [cd.okubogar_detector(nochange_before_geo, nochange_after_geo, extractor_kwargs=params[0], cluster_params=params[1]) for params in parameters]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 181,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "nx = len(extractor_params)\n",
    "ny = len(cluster_params)\n",
    "\n",
    "fig, axs = plt.subplots(ny, nx, sharex='col', sharey='row')\n",
    "fig.set_figheight(50)\n",
    "fig.set_figwidth(20)\n",
    "for ax in axs.flat:\n",
    "    ax.set(xlabel='x-label', ylabel='y-label')\n",
    "    \n",
    "xs, ys = np.mgrid[:ny, :nx]\n",
    "cmap = matplotlib.cm.get_cmap('plasma')\n",
    "\n",
    "for i, coord in enumerate(zip(xs.flatten(),ys.flatten())):\n",
    "    x,y = coord\n",
    "    polys, weights, diff = no_change_results[i]\n",
    "    \n",
    "    indices = [i for (i,w) in enumerate(weights) if w > 6]\n",
    "    filtered_polys = []\n",
    "    filtered_polys = [polys[i] for i in indices]\n",
    "    \n",
    "    axs[x,y].set(xlabel=f\"ex {parameters[i][1]}\", ylabel=f\"cl {parameters[i][0]}\")\n",
    "    axs[x,y].imshow(diff, cmap=\"Greys_r\")\n",
    "    axs[x,y].imshow(nochange_after_arr, cmap=\"Greys_r\", alpha=.5)\n",
    "\n",
    "    if weights:\n",
    "        colors = bytescale(np.asarray(weights))\n",
    "        ax = [axs[x,y].fill(*p.exterior.xy, c=\"r\", alpha=.6) for i,p in enumerate(filtered_polys)]\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 80,
   "metadata": {},
   "outputs": [],
   "source": [
    "nc_weights = []\n",
    "for r in no_change_results:\n",
    "    nc_weights.extend(r[1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 176,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(-0.2109375, 4.7621055)"
      ]
     },
     "execution_count": 176,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "min(nc_weights), max(nc_weights), mean()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 81,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "bs_weights = []\n",
    "for r in boulder_shift_results:\n",
    "    bs_weights.extend(r[1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 182,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(1, figsize=(10,5))\n",
    "_ = plt.hist(nc_weights, color=\"orange\", alpha=.5, label=\"no change\")\n",
    "plt.legend(loc='upper right')\n",
    "plt.figure(2, figsize=(10,5))\n",
    "_ = plt.hist(bs_weights, alpha=.5, label=\"change\")\n",
    "plt.legend(loc='upper right')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 183,
   "metadata": {},
   "outputs": [],
   "source": [
    "# plt.figure(1, figsize=(10,5))\n",
    "# _ = plt.hist(boulder_shift_results[0][2].flatten(), bins=60, alpha=.6, label=\"change\")\n",
    "# _ = plt.hist(no_change_results[0][2].flatten(), bins=40, alpha=.6, label=\"no change\")\n",
    "# plt.legend(loc='upper right')\n",
    "\n",
    "plt.figure(1, figsize=(10,5))\n",
    "_ = plt.hist(all_weights, bins=20, alpha=.6, label=\"change\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 166,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(-43.0, 67.0, 9.827137)"
      ]
     },
     "execution_count": 166,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.min(boulder_shift_results[0][2]), np.max(boulder_shift_results[0][2]), np.mean(boulder_shift_results[0][2])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 167,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(-16.0, 25.0, 2.4594834)"
      ]
     },
     "execution_count": 167,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.min(no_change_results[0][2]), np.max(no_change_results[0][2]), np.mean(no_change_results[0][2])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "autocnetdev",
   "language": "python",
   "name": "autocnetdev"
  },
  "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.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}