KDEExplorersApp_withweight.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Explore Kernel Density Estimate Fitting on WNS 'Endemic Area'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2017-12-12T22:59:37.291663Z",
     "start_time": "2017-12-12T22:59:34.813470Z"
    }
   },
   "outputs": [],
   "source": [
    "import os\n",
    "import ipywidgets as widgets\n",
    "from ipywidgets import interact\n",
    "\n",
    "from osgeo import gdal\n",
    "import geopandas as gpd\n",
    "import numpy as np\n",
    "from scipy import stats\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "from shapely.geometry import mapping, Polygon, MultiPolygon, LineString, Point"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2017-12-12T22:59:37.883904Z",
     "start_time": "2017-12-12T22:59:37.762650Z"
    }
   },
   "outputs": [],
   "source": [
    "wns_fname = os.path.join('data', \"WNS_Status_Provisional_8_30_2018.shp\")\n",
    "wns = gpd.GeoDataFrame.from_file(wns_fname)\n",
    "wns['year'] = wns.WNS_MAP_YR.str[:4].astype(int)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2017-12-12T22:59:37.915971Z",
     "start_time": "2017-12-12T22:59:37.890919Z"
    }
   },
   "outputs": [],
   "source": [
    "target_crs = {'proj': 'aea',\n",
    " 'lat_1': 29.5,\n",
    " 'lat_2': 45.5,\n",
    " 'lat_0': 23,\n",
    " 'lon_0': -96,\n",
    " 'x_0': 0,\n",
    " 'y_0': 0,\n",
    " 'datum': 'NAD83',\n",
    " 'units': 'm',\n",
    " 'no_defs': True}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2017-12-12T22:59:37.986118Z",
     "start_time": "2017-12-12T22:59:37.959062Z"
    }
   },
   "outputs": [],
   "source": [
    "states_fname = os.path.join('data', \"cb_2017_us_state_20m.shp\")\n",
    "states = gpd.read_file(states_fname)\n",
    "states_aea = states.to_crs(target_crs)\n",
    "conus_states = states_aea[~states_aea.NAME.isin(['Alaska', 'Hawaii', 'Puerto Rico'])]\n",
    "\n",
    "bounds = conus_states.bounds\n",
    "x_bounds = [bounds.minx.min(), bounds.maxx.max()]\n",
    "y_bounds = [bounds.miny.min(), bounds.maxy.max()]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2017-12-12T22:59:38.023196Z",
     "start_time": "2017-12-12T22:59:37.996139Z"
    }
   },
   "outputs": [],
   "source": [
    "centroids = wns.centroid\n",
    "centroids = centroids.to_crs(target_crs)\n",
    "wns['x'] = centroids.x\n",
    "wns['y'] = centroids.y\n",
    "\n",
    "wns['area'] = centroids.to_crs(target_crs).geometry.area\n",
    "wns['area_weight'] = wns.area/wns.area.max()\n",
    "# wns.head(5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2017-12-12T22:59:38.059272Z",
     "start_time": "2017-12-12T22:59:38.034219Z"
    }
   },
   "outputs": [],
   "source": [
    "def reclass_as_pcnt(Z):\n",
    "    uniques = np.unique(Z)[::-1]\n",
    "    z_reclass = np.copy(Z)\n",
    "    total_dens = Z.sum()\n",
    "    cum_area = 0.0\n",
    "    for u in uniques:\n",
    "        cum_area += np.count_nonzero(Z==u)*u\n",
    "        z_reclass[Z==u] = cum_area/total_dens\n",
    "    return z_reclass\n",
    "    #area_lookup[unique] = np.count_nonzero(z==unique)\n",
    "\n",
    "def create_kernel(x, y, X, Y, factor=1.2, weights=None):\n",
    "    \n",
    "    positions = np.vstack([X.ravel(), Y.ravel()])\n",
    "    values = np.vstack([x, y])\n",
    "    kernel = gaussian_kde(values, weights=weights)#, bw_method='silverman')#, bw_method=0.15/np.asarray(values).std(ddof=1))\n",
    "#     kernel.silverman_factor()\n",
    "#     kernel.set_bandwidth(bw_method='silverman')\n",
    "    kernel.set_bandwidth( factor)\n",
    "    # kernel.set_bandwidth(bw_method=bw_method)\n",
    "    # kernel.set_bandwidth(bw_method=0.15/np.asarray(values).std(ddof=1))\n",
    "    Z = np.reshape(kernel(positions).T, X.shape)\n",
    "    return reclass_as_pcnt(Z)\n",
    "\n",
    "def create_kernel_contours(x, y, z, levels=[0.5, 0.75, 0.95]):\n",
    "    cset = plt.contour(x, y, z, levels=levels, colors=['red', 'white', 'blue'])\n",
    "    return cset\n",
    "\n",
    "\n",
    "\n",
    "def plot_one(x, y, X, Y, Z, title='', isopleth=0.75):\n",
    "    fig = plt.figure(figsize=(25, 15))\n",
    "    ax = fig.add_subplot(111)\n",
    "\n",
    "    cset = create_kernel_contours(X, Y, Z, levels=[isopleth])   \n",
    "    \n",
    "#     ax.imshow(np.rot90(Z), cmap=plt.cm.Reds_r, \n",
    "#                extent=[X.min(), X.max(), Y.min(), Y.max()], alpha=0.6)\n",
    "\n",
    "    hull = Polygon(zip(x, y)).convex_hull\n",
    "    pnts = np.array(np.asarray(hull.exterior)).transpose()\n",
    "    \n",
    "#     ax.plot(pnts[0],pnts[1])\n",
    "    \n",
    "#     ax.plot(x, y, 'k.', markersize=6, alpha=0.4)\n",
    "#     ax.set_xlim([X.min(), X.max()])\n",
    "#     ax.set_ylim([Y.min(), Y.max()])\n",
    "    ax.set_xlim([-3000000, 3000000])\n",
    "    ax.set_ylim([0, 3700000])\n",
    "#     plt.colorbar()\n",
    "    states_aea.plot(color='None',  edgecolor='black', ax=ax, alpha=0.4)\n",
    "    plt.title(title,fontsize=30)\n",
    "#     plt.show()\n",
    "    ax.set_aspect('equal')\n",
    "    return ax"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# copied from https://gist.github.com/tillahoffmann/f844bce2ec264c1c8cb5\n",
    "\n",
    "import numpy as np\n",
    "from scipy.spatial.distance import cdist\n",
    "\n",
    "class gaussian_kde(object):\n",
    "    \"\"\"Representation of a kernel-density estimate using Gaussian kernels.\n",
    "\n",
    "    Kernel density estimation is a way to estimate the probability density\n",
    "    function (PDF) of a random variable in a non-parametric way.\n",
    "    `gaussian_kde` works for both uni-variate and multi-variate data.   It\n",
    "    includes automatic bandwidth determination.  The estimation works best for\n",
    "    a unimodal distribution; bimodal or multi-modal distributions tend to be\n",
    "    oversmoothed.\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    dataset : array_like\n",
    "        Datapoints to estimate from. In case of univariate data this is a 1-D\n",
    "        array, otherwise a 2-D array with shape (# of dims, # of data).\n",
    "    bw_method : str, scalar or callable, optional\n",
    "        The method used to calculate the estimator bandwidth.  This can be\n",
    "        'scott', 'silverman', a scalar constant or a callable.  If a scalar,\n",
    "        this will be used directly as `kde.factor`.  If a callable, it should\n",
    "        take a `gaussian_kde` instance as only parameter and return a scalar.\n",
    "        If None (default), 'scott' is used.  See Notes for more details.\n",
    "    weights : array_like, shape (n, ), optional, default: None\n",
    "        An array of weights, of the same shape as `x`.  Each value in `x`\n",
    "        only contributes its associated weight towards the bin count\n",
    "        (instead of 1).\n",
    "\n",
    "    Attributes\n",
    "    ----------\n",
    "    dataset : ndarray\n",
    "        The dataset with which `gaussian_kde` was initialized.\n",
    "    d : int\n",
    "        Number of dimensions.\n",
    "    n : int\n",
    "        Number of datapoints.\n",
    "    neff : float\n",
    "        Effective sample size using Kish's approximation.\n",
    "    factor : float\n",
    "        The bandwidth factor, obtained from `kde.covariance_factor`, with which\n",
    "        the covariance matrix is multiplied.\n",
    "    covariance : ndarray\n",
    "        The covariance matrix of `dataset`, scaled by the calculated bandwidth\n",
    "        (`kde.factor`).\n",
    "    inv_cov : ndarray\n",
    "        The inverse of `covariance`.\n",
    "\n",
    "    Methods\n",
    "    -------\n",
    "    kde.evaluate(points) : ndarray\n",
    "        Evaluate the estimated pdf on a provided set of points.\n",
    "    kde(points) : ndarray\n",
    "        Same as kde.evaluate(points)\n",
    "    kde.pdf(points) : ndarray\n",
    "        Alias for ``kde.evaluate(points)``.\n",
    "    kde.set_bandwidth(bw_method='scott') : None\n",
    "        Computes the bandwidth, i.e. the coefficient that multiplies the data\n",
    "        covariance matrix to obtain the kernel covariance matrix.\n",
    "        .. versionadded:: 0.11.0\n",
    "    kde.covariance_factor : float\n",
    "        Computes the coefficient (`kde.factor`) that multiplies the data\n",
    "        covariance matrix to obtain the kernel covariance matrix.\n",
    "        The default is `scotts_factor`.  A subclass can overwrite this method\n",
    "        to provide a different method, or set it through a call to\n",
    "        `kde.set_bandwidth`.\n",
    "\n",
    "    Notes\n",
    "    -----\n",
    "    Bandwidth selection strongly influences the estimate obtained from the KDE\n",
    "    (much more so than the actual shape of the kernel).  Bandwidth selection\n",
    "    can be done by a \"rule of thumb\", by cross-validation, by \"plug-in\n",
    "    methods\" or by other means; see [3]_, [4]_ for reviews.  `gaussian_kde`\n",
    "    uses a rule of thumb, the default is Scott's Rule.\n",
    "\n",
    "    Scott's Rule [1]_, implemented as `scotts_factor`, is::\n",
    "\n",
    "        n**(-1./(d+4)),\n",
    "\n",
    "    with ``n`` the number of data points and ``d`` the number of dimensions.\n",
    "    Silverman's Rule [2]_, implemented as `silverman_factor`, is::\n",
    "\n",
    "        (n * (d + 2) / 4.)**(-1. / (d + 4)).\n",
    "\n",
    "    Good general descriptions of kernel density estimation can be found in [1]_\n",
    "    and [2]_, the mathematics for this multi-dimensional implementation can be\n",
    "    found in [1]_.\n",
    "\n",
    "    References\n",
    "    ----------\n",
    "    .. [1] D.W. Scott, \"Multivariate Density Estimation: Theory, Practice, and\n",
    "           Visualization\", John Wiley & Sons, New York, Chicester, 1992.\n",
    "    .. [2] B.W. Silverman, \"Density Estimation for Statistics and Data\n",
    "           Analysis\", Vol. 26, Monographs on Statistics and Applied Probability,\n",
    "           Chapman and Hall, London, 1986.\n",
    "    .. [3] B.A. Turlach, \"Bandwidth Selection in Kernel Density Estimation: A\n",
    "           Review\", CORE and Institut de Statistique, Vol. 19, pp. 1-33, 1993.\n",
    "    .. [4] D.M. Bashtannyk and R.J. Hyndman, \"Bandwidth selection for kernel\n",
    "           conditional density estimation\", Computational Statistics & Data\n",
    "           Analysis, Vol. 36, pp. 279-298, 2001.\n",
    "\n",
    "    Examples\n",
    "    --------\n",
    "    Generate some random two-dimensional data:\n",
    "\n",
    "    >>> from scipy import stats\n",
    "    >>> def measure(n):\n",
    "    >>>     \"Measurement model, return two coupled measurements.\"\n",
    "    >>>     m1 = np.random.normal(size=n)\n",
    "    >>>     m2 = np.random.normal(scale=0.5, size=n)\n",
    "    >>>     return m1+m2, m1-m2\n",
    "\n",
    "    >>> m1, m2 = measure(2000)\n",
    "    >>> xmin = m1.min()\n",
    "    >>> xmax = m1.max()\n",
    "    >>> ymin = m2.min()\n",
    "    >>> ymax = m2.max()\n",
    "\n",
    "    Perform a kernel density estimate on the data:\n",
    "\n",
    "    >>> X, Y = np.mgrid[xmin:xmax:100j, ymin:ymax:100j]\n",
    "    >>> positions = np.vstack([X.ravel(), Y.ravel()])\n",
    "    >>> values = np.vstack([m1, m2])\n",
    "    >>> kernel = stats.gaussian_kde(values)\n",
    "    >>> Z = np.reshape(kernel(positions).T, X.shape)\n",
    "\n",
    "    Plot the results:\n",
    "\n",
    "    >>> import matplotlib.pyplot as plt\n",
    "    >>> fig = plt.figure()\n",
    "    >>> ax = fig.add_subplot(111)\n",
    "    >>> ax.imshow(np.rot90(Z), cmap=plt.cm.gist_earth_r,\n",
    "    ...           extent=[xmin, xmax, ymin, ymax])\n",
    "    >>> ax.plot(m1, m2, 'k.', markersize=2)\n",
    "    >>> ax.set_xlim([xmin, xmax])\n",
    "    >>> ax.set_ylim([ymin, ymax])\n",
    "    >>> plt.show()\n",
    "\n",
    "    \"\"\"\n",
    "    def __init__(self, dataset, bw_method=None, weights=None):\n",
    "        self.dataset = np.atleast_2d(dataset)\n",
    "        if not self.dataset.size > 1:\n",
    "            raise ValueError(\"`dataset` input should have multiple elements.\")\n",
    "        self.d, self.n = self.dataset.shape\n",
    "            \n",
    "        if weights is not None:\n",
    "            self.weights = weights / np.sum(weights)\n",
    "        else:\n",
    "            self.weights = np.ones(self.n) / self.n\n",
    "            \n",
    "        # Compute the effective sample size \n",
    "        # http://surveyanalysis.org/wiki/Design_Effects_and_Effective_Sample_Size#Kish.27s_approximate_formula_for_computing_effective_sample_size\n",
    "        self.neff = 1.0 / np.sum(self.weights ** 2)\n",
    "\n",
    "        self.set_bandwidth(bw_method=bw_method)\n",
    "\n",
    "    def evaluate(self, points):\n",
    "        \"\"\"Evaluate the estimated pdf on a set of points.\n",
    "\n",
    "        Parameters\n",
    "        ----------\n",
    "        points : (# of dimensions, # of points)-array\n",
    "            Alternatively, a (# of dimensions,) vector can be passed in and\n",
    "            treated as a single point.\n",
    "\n",
    "        Returns\n",
    "        -------\n",
    "        values : (# of points,)-array\n",
    "            The values at each point.\n",
    "\n",
    "        Raises\n",
    "        ------\n",
    "        ValueError : if the dimensionality of the input points is different than\n",
    "                     the dimensionality of the KDE.\n",
    "\n",
    "        \"\"\"\n",
    "        points = np.atleast_2d(points)\n",
    "\n",
    "        d, m = points.shape\n",
    "        if d != self.d:\n",
    "            if d == 1 and m == self.d:\n",
    "                # points was passed in as a row vector\n",
    "                points = np.reshape(points, (self.d, 1))\n",
    "                m = 1\n",
    "            else:\n",
    "                msg = \"points have dimension %s, dataset has dimension %s\" % (d,\n",
    "                    self.d)\n",
    "                raise ValueError(msg)\n",
    "\n",
    "        # compute the normalised residuals\n",
    "        chi2 = cdist(points.T, self.dataset.T, 'mahalanobis', VI=self.inv_cov) ** 2\n",
    "        # compute the pdf\n",
    "        result = np.sum(np.exp(-.5 * chi2) * self.weights, axis=1) / self._norm_factor\n",
    "\n",
    "        return result\n",
    "\n",
    "    __call__ = evaluate\n",
    "\n",
    "    def scotts_factor(self):\n",
    "        return np.power(self.neff, -1./(self.d+4))\n",
    "\n",
    "    def silverman_factor(self):\n",
    "        return np.power(self.neff*(self.d+2.0)/4.0, -1./(self.d+4))\n",
    "\n",
    "    #  Default method to calculate bandwidth, can be overwritten by subclass\n",
    "    covariance_factor = scotts_factor\n",
    "\n",
    "    def set_bandwidth(self, bw_method=None):\n",
    "        \"\"\"Compute the estimator bandwidth with given method.\n",
    "\n",
    "        The new bandwidth calculated after a call to `set_bandwidth` is used\n",
    "        for subsequent evaluations of the estimated density.\n",
    "\n",
    "        Parameters\n",
    "        ----------\n",
    "        bw_method : str, scalar or callable, optional\n",
    "            The method used to calculate the estimator bandwidth.  This can be\n",
    "            'scott', 'silverman', a scalar constant or a callable.  If a\n",
    "            scalar, this will be used directly as `kde.factor`.  If a callable,\n",
    "            it should take a `gaussian_kde` instance as only parameter and\n",
    "            return a scalar.  If None (default), nothing happens; the current\n",
    "            `kde.covariance_factor` method is kept.\n",
    "\n",
    "        Notes\n",
    "        -----\n",
    "        .. versionadded:: 0.11\n",
    "\n",
    "        Examples\n",
    "        --------\n",
    "        >>> x1 = np.array([-7, -5, 1, 4, 5.])\n",
    "        >>> kde = stats.gaussian_kde(x1)\n",
    "        >>> xs = np.linspace(-10, 10, num=50)\n",
    "        >>> y1 = kde(xs)\n",
    "        >>> kde.set_bandwidth(bw_method='silverman')\n",
    "        >>> y2 = kde(xs)\n",
    "        >>> kde.set_bandwidth(bw_method=kde.factor / 3.)\n",
    "        >>> y3 = kde(xs)\n",
    "\n",
    "        >>> fig = plt.figure()\n",
    "        >>> ax = fig.add_subplot(111)\n",
    "        >>> ax.plot(x1, np.ones(x1.shape) / (4. * x1.size), 'bo',\n",
    "        ...         label='Data points (rescaled)')\n",
    "        >>> ax.plot(xs, y1, label='Scott (default)')\n",
    "        >>> ax.plot(xs, y2, label='Silverman')\n",
    "        >>> ax.plot(xs, y3, label='Const (1/3 * Silverman)')\n",
    "        >>> ax.legend()\n",
    "        >>> plt.show()\n",
    "\n",
    "        \"\"\"\n",
    "        if bw_method is None:\n",
    "            pass\n",
    "        elif bw_method == 'scott':\n",
    "            self.covariance_factor = self.scotts_factor\n",
    "        elif bw_method == 'silverman':\n",
    "            self.covariance_factor = self.silverman_factor\n",
    "        elif np.isscalar(bw_method):\n",
    "            self._bw_method = 'use constant'\n",
    "            self.covariance_factor = lambda: self.factor / bw_method\n",
    "        elif callable(bw_method):\n",
    "            self._bw_method = bw_method\n",
    "            self.covariance_factor = lambda: self._bw_method(self)\n",
    "        else:\n",
    "            msg = \"`bw_method` should be 'scott', 'silverman', a scalar \" \\\n",
    "                  \"or a callable.\"\n",
    "            raise ValueError(msg)\n",
    "\n",
    "        self._compute_covariance()\n",
    "\n",
    "    def _compute_covariance(self):\n",
    "        \"\"\"Computes the covariance matrix for each Gaussian kernel using\n",
    "        covariance_factor().\n",
    "        \"\"\"\n",
    "        self.factor = self.covariance_factor()\n",
    "        # Cache covariance and inverse covariance of the data\n",
    "        if not hasattr(self, '_data_inv_cov'):\n",
    "            # Compute the mean and residuals\n",
    "            _mean = np.sum(self.weights * self.dataset, axis=1)\n",
    "            _residual = (self.dataset - _mean[:, None])\n",
    "            # Compute the biased covariance\n",
    "            self._data_covariance = np.atleast_2d(np.dot(_residual * self.weights, _residual.T))\n",
    "            # Correct for bias (http://en.wikipedia.org/wiki/Weighted_arithmetic_mean#Weighted_sample_covariance)\n",
    "            self._data_covariance /= (1 - np.sum(self.weights ** 2))\n",
    "            self._data_inv_cov = np.linalg.inv(self._data_covariance)\n",
    "\n",
    "        self.covariance = self._data_covariance * self.factor**2\n",
    "        self.inv_cov = self._data_inv_cov / self.factor**2\n",
    "        self._norm_factor = np.sqrt(np.linalg.det(2*np.pi*self.covariance)) #* self.n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2017-12-12T22:59:37.951045Z",
     "start_time": "2017-12-12T22:59:37.925992Z"
    }
   },
   "outputs": [],
   "source": [
    "cache = {}\n",
    "\n",
    "def interactive_plot(factor, iso=0.75, year=2017, \n",
    "                     add_counties=False, \n",
    "                     add_centroids=True,\n",
    "                     add_surface=True,\n",
    "                     weight_by_county_area=True,\n",
    "                     weight_by_year=True):\n",
    "    if year in cache:\n",
    "        X, Y, data = cache[year]['data']\n",
    "    else:\n",
    "        cache[year] = {}\n",
    "        cache[year]['kernels'] = {}\n",
    "        data = wns[wns.year <= year]\n",
    "        \n",
    "        pad=500000\n",
    "        X, Y = np.mgrid[x_bounds[0]-pad:x_bounds[1]+pad:100j, y_bounds[0]-pad:y_bounds[1]+pad:100j]\n",
    "        cache[year]['data'] = (X, Y, data)\n",
    "    \n",
    "    \n",
    "    if weight_by_year:\n",
    "        weights = gen_weights(data.years_pres, factor=weight_f.value)\n",
    "    else:\n",
    "        weights = np.ones(data.area_weight.shape)\n",
    "    weights /= np.sum(weights)\n",
    "    \n",
    "    if weight_by_county_area:\n",
    "        county_weights = data.area_weight\n",
    "        county_weights /= np.sum(county_weights)\n",
    "        weights += county_weights\n",
    "        \n",
    "    weights /= np.sum(weights)\n",
    "    weights = list(weights)\n",
    "\n",
    "    Z = create_kernel(data.x, data.y, X, Y, factor, weights=weights)\n",
    "    cache[year]['kernels'][factor] = Z\n",
    "        \n",
    "        \n",
    "        \n",
    "        \n",
    "    ax = plot_one(data.x, data.y, X, Y, Z, isopleth=iso, title=f\"WNS {year}-{year+1}\")\n",
    "    \n",
    "    if add_surface:\n",
    "        ax.imshow(np.rot90(Z), cmap=plt.cm.Reds_r, \n",
    "           extent=[X.min(), X.max(), Y.min(), Y.max()], alpha=0.6)\n",
    "    \n",
    "    if add_counties:\n",
    "        data.to_crs(states_aea.crs).plot(column='year', cmap='plasma', ax=ax, legend=True)\n",
    "\n",
    "    if add_centroids:\n",
    "        data.to_crs(states_aea.crs).centroid.plot(color='grey', ax=ax)\n",
    "\n",
    "        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2017-12-12T22:59:38.187540Z",
     "start_time": "2017-12-12T22:59:38.110379Z"
    }
   },
   "outputs": [],
   "source": [
    "f = widgets.FloatSlider(\n",
    "    value=1.2,\n",
    "    min=0.1,\n",
    "    max=4.0,\n",
    "    step=0.1,\n",
    "    description=\"BW Factor\")\n",
    "\n",
    "f.continuous_update = False"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "y = widgets.IntSlider(\n",
    "    value=2014,\n",
    "    min=wns.year.min()+1,\n",
    "    max=wns.year.max(),\n",
    "    step=1,\n",
    "    description=\"Year\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "i = widgets.FloatSlider(\n",
    "    value=0.75,\n",
    "    min=0.05,\n",
    "    max=.99,\n",
    "    step=0.05,\n",
    "    description=\"Isopleth\")\n",
    "\n",
    "i.continuous_update = False"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "c = widgets.Checkbox(value=False, description='Show Counties')\n",
    "c2 = widgets.Checkbox(value=False, description='Show Centroids')\n",
    "c3 = widgets.Checkbox(value=False, description='Show KDE Surface')\n",
    "c4 = widgets.Checkbox(value=True, description='Weight by county areas')\n",
    "c5 = widgets.Checkbox(value=True, description='Weight by years since detection')\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### In order to weight counties we need to decide on a function to use for weighting.\n",
    "#### Let's start with a simple number of years raised to a power function, as it allows us to explore a range of curves"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*move the slider below to set the shape of this curve*\n",
    "* 0.0 is no weight\n",
    "* 1.0 is linear\n",
    "\n",
    "***Note: Changes to this curve do not update the map until one of the widgets below is triggered***"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'widgets' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-1-271b28361b9e>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m()\u001b[0m\n\u001b[1;32m----> 1\u001b[1;33m weight_f = widgets.FloatSlider(\n\u001b[0m\u001b[0;32m      2\u001b[0m     \u001b[0mvalue\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;36m2\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      3\u001b[0m     \u001b[0mmin\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;36m0.0\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      4\u001b[0m     \u001b[0mmax\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;36m10.0\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      5\u001b[0m     \u001b[0mstep\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;36m0.1\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;31mNameError\u001b[0m: name 'widgets' is not defined"
     ]
    }
   ],
   "source": [
    "weight_f = widgets.FloatSlider(\n",
    "    value=2,\n",
    "    min=0.0,\n",
    "    max=10.0,\n",
    "    step=0.1,\n",
    "    description=\"Weight Exponent\",\n",
    "layout={'width': '500px'})\n",
    "\n",
    "def exponential(x, factor=2.0):\n",
    "    return x**factor\n",
    "\n",
    "wns['years_pres'] = wns.year.max()-wns.year\n",
    "def gen_weights(years_pres, factor=2.0):\n",
    "    weights = [exponential(weight, factor) for weight in years_pres]\n",
    "    \n",
    "#     weights = [w/np.max(weights) for w in weights]\n",
    "    return weights\n",
    "\n",
    "def show_weight_curve(factor=2):\n",
    "    years = list(range(0, wns.years_pres.max()))\n",
    "    weights = gen_weights(years, factor=factor)\n",
    "    weights = weights / np.max(weights)\n",
    "    f, ax1 = plt.subplots(1, 1, figsize=(5,5))\n",
    "    \n",
    "    ax1.plot(years, weights)\n",
    "    ax1.set_title('County Weights by years since detection')\n",
    "    ax1.set_xlabel('Years since Detection (n)')\n",
    "    ax1.set_ylabel('Weight')\n",
    "    \n",
    "    ax1.set_ylim(0,1.1)\n",
    "    \n",
    "    tex = r'$n^{' + str(factor) + '}$'\n",
    "    ax1.text(1, 0.8, tex, fontsize=20, va='bottom')\n",
    "    \n",
    "    plt.tight_layout()\n",
    "\n",
    "_ = interact(show_weight_curve, factor=weight_f)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Once the weight function is finished used the sliders below to see the effect on our KDE"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "c76ac88877014168b6ef8030fb339a81",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(FloatSlider(value=1.2, continuous_update=False, description='BW Factor', max=4.0, min=0.…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "_ = interact(interactive_plot, factor=f, iso=i, year=y, \n",
    "             add_counties=c,\n",
    "             add_centroids=c2,\n",
    "             add_surface=c3,\n",
    "             weight_by_county_area=c4,\n",
    "             weight_by_years=c5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "i.continuous_update = False"
   ]
  },
  {
   "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.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}