10_CompareClusteringsPenaltys.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# How to handle Nan-values so that the HA doesn't get marginalized?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It has been a problem this far, that the clustering doesn't work as desired, and the problem is now located in the procedure Nan-values have been handled in clustering. So we need a better way to do that."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    " %matplotlib inline\n",
    "import seaborn\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import os\n",
    "from gradutil import *\n",
    "from pyomo.opt import SolverFactory"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "seedn = 2\n",
    "opt = SolverFactory('glpk')\n",
    "solutions = ideal()\n",
    "revenue, carbon, deadwood, ha = init_boreal()\n",
    "x = np.concatenate((revenue, carbon, deadwood, ha), axis=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Set nan:s to smallest existing option"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's first just try setting all nan:s to the smallest value in the corresponding column"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "norm_data = x.copy()\n",
    "inds = np.where(np.isnan(norm_data))\n",
    "norm_data[inds] = np.take(np.nanmin(norm_data, axis=0),inds[1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then normalize all as before"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "min_norm_x = normalize(norm_data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "%pylab inline\n",
    "pylab.rcParams['figure.figsize'] = (15,12)\n",
    "\n",
    "def hist_plot_norm(data, ax, limits):\n",
    "    ax[0,0].hist(data[:, :7])\n",
    "    ax[0,0].axis(limits)\n",
    "    ax[0,0].set_title('Timber Harvest Revenues')\n",
    "\n",
    "    ax[0,1].hist(data[:, 7:14])\n",
    "    ax[0,1].axis(limits)\n",
    "    ax[0,1].set_title('Carbon storage')\n",
    "\n",
    "    ax[1,0].hist(data[:, 14:21])\n",
    "    ax[1,0].axis(limits)\n",
    "    ax[1,0].set_title('Deadwood')\n",
    "\n",
    "    ax[1,1].hist(data[:, 21:])\n",
    "    ax[1,1].axis(limits)\n",
    "    ax[1,1].set_title('Habitat availability')\n",
    "    return ax\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = min_norm_x\n",
    "fig, ax = plt.subplots(2,2)\n",
    "limits = [.0, 1., 0, 30000]\n",
    "hist_plot_norm(data, ax, limits)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "nclust = 50\n",
    "optim_revenue50, optim_carbon50, optim_deadwood50, optim_ha50 = cNopt(x, min_norm_x, min_norm_x, opt, nclust, seedn)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "print('Relative differences to original values, 50 clusters')\n",
    "print(\"(i) Harvest revenues difference {:.3f}\".format((optim_revenue50-solutions['revenue'])/solutions['revenue']))\n",
    "print(\"(ii) Carbon storage {:.3f}\".format((optim_carbon50-solutions['carbon'])/solutions['carbon']))\n",
    "print(\"(iii) Deadwood index {:.3f}\".format((optim_deadwood50-solutions['deadwood'])/solutions['deadwood']))\n",
    "print(\"(iv) Combined Habitat {:.3f}\".format((optim_ha50-solutions['ha'])/solutions['ha']))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So it looks like this setting is not enough to drive optimization away from these points, and it doesn't tell anything about clustering. We need to adjust values for the optimization part, so we can know how the clustering goes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "no_nan_x = x.copy()\n",
    "inds = np.where(np.isnan(no_nan_x))\n",
    "no_nan_x[inds] = np.take(np.nanmin(no_nan_x, axis=0) - np.nanmax(no_nan_x, axis=0), inds[1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There is then great penalty for choosing the Nan values in optimization."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "nclust = 50\n",
    "penalty_optim_revenue50, penalty_optim_carbon50, penalty_optim_deadwood50, penalty_optim_ha50 = cNopt(x, min_norm_x, no_nan_x, opt, nclust, seedn)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('Relative differences to original values, 50 clusters')\n",
    "print(\"(i) Harvest revenues difference {:.3f}\".format((penalty_optim_revenue50-solutions['revenue'])/solutions['revenue']))\n",
    "print(\"(ii) Carbon storage {:.3f}\".format((penalty_optim_carbon50-solutions['carbon'])/solutions['carbon']))\n",
    "print(\"(iii) Deadwood index {:.3f}\".format((penalty_optim_deadwood50-solutions['deadwood'])/solutions['deadwood']))\n",
    "print(\"(iv) Combined Habitat {:.3f}\".format((penalty_optim_ha50-solutions['ha'])/solutions['ha']))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We have optimization results and it looks like the clustering is not working. We need another paradigma to handle the Nan-values..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Give nan:s some penalty"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "norm_data = x.copy()\n",
    "inds = np.where(np.isnan(norm_data))\n",
    "norm_data[inds] = np.take((np.nanmin(norm_data, axis=0)-np.nanmax(norm_data, axis=0))/2,inds[1])\n",
    "penalty_norm_x = normalize(norm_data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(2,2)\n",
    "limits = [.0, 1., 0, 30000]\n",
    "hist_plot_norm(penalty_norm_x, ax, limits)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "nclust = 50\n",
    "half_optim_revenue50, half_optim_carbon50, half_optim_deadwood50, half_optim_ha50 = cNopt(x, penalty_norm_x, no_nan_x, opt, nclust, seedn)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('Relative differences to original values, 50 clusters')\n",
    "print(\"(i) Harvest revenues difference {:.3f}\".format((half_optim_revenue50-solutions['revenue'])/solutions['revenue']))\n",
    "print(\"(ii) Carbon storage {:.3f}\".format((half_optim_carbon50-solutions['carbon'])/solutions['carbon']))\n",
    "print(\"(iii) Deadwood index {:.3f}\".format((half_optim_deadwood50-solutions['deadwood'])/solutions['deadwood']))\n",
    "print(\"(iv) Combined Habitat {:.3f}\".format((half_optim_ha50-solutions['ha'])/solutions['ha']))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "That is not working either. Need something else."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Give nans ridiculous penalty"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "norm_data = x.copy()\n",
    "inds = np.where(np.isnan(norm_data))\n",
    "norm_data[inds] = np.take((np.nanmin(norm_data, axis=0)-np.nanmax(norm_data, axis=0))*2,inds[1])\n",
    "ridiculous_norm_x = normalize(norm_data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(2,2)\n",
    "limits = [.0, 1., 0, 30000]\n",
    "hist_plot_norm(ridiculous_norm_x, ax, limits)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "nclust = 50\n",
    "ridic_optim_revenue50, ridic_optim_carbon50, ridic_optim_deadwood50, ridic_optim_ha50 = cNopt(x, ridiculous_norm_x, no_nan_x, opt, nclust, seedn)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('Relative differences to original values, 50 clusters')\n",
    "print(\"(i) Harvest revenues difference {:.3f}\".format((ridic_optim_revenue50-solutions['revenue'])/solutions['revenue']))\n",
    "print(\"(ii) Carbon storage {:.3f}\".format((ridic_optim_carbon50-solutions['carbon'])/solutions['carbon']))\n",
    "print(\"(iii) Deadwood index {:.3f}\".format((ridic_optim_deadwood50-solutions['deadwood'])/solutions['deadwood']))\n",
    "print(\"(iv) Combined Habitat {:.3f}\".format((ridic_optim_ha50-solutions['ha'])/solutions['ha']))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Cluster Nan:s separately"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After much effort there rose new hope with an idea to cluster all the Nan values independently! Then there will be one more parameter to decide, but that is small price to pay from success."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "x_nan = x[any(np.isnan(x), axis=1),:]\n",
    "np.shape(x_nan)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "x_num = x[all(~np.isnan(x), axis=1),:]\n",
    "np.shape(x_num)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's cluster all the different \"nan-scenarios\" differently. Details on how the nan:s relate are found in the DataTesting notebook."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "x_nany2y4y6 = x_nan[np.logical_and(np.logical_and(np.isnan(x_nan[:,2]),np.isnan(x_nan[:,4])), np.isnan(x_nan[:,6])),:]\n",
    "x_nany2y4n6 = x_nan[np.logical_and(np.logical_and(np.isnan(x_nan[:,2]),np.isnan(x_nan[:,4])), ~np.isnan(x_nan[:,6])),:]\n",
    "x_nany2n4y6 = x_nan[np.logical_and(np.logical_and(np.isnan(x_nan[:,2]),~np.isnan(x_nan[:,4])), np.isnan(x_nan[:,6])),:]\n",
    "x_nany2n4n6 = x_nan[np.logical_and(np.logical_and(np.isnan(x_nan[:,2]),~np.isnan(x_nan[:,4])), ~np.isnan(x_nan[:,6])),:]\n",
    "x_nann2y4y6 = x_nan[np.logical_and(np.logical_and(~np.isnan(x_nan[:,2]),np.isnan(x_nan[:,4])), np.isnan(x_nan[:,6])),:]\n",
    "x_nann2y4n6 = x_nan[np.logical_and(np.logical_and(~np.isnan(x_nan[:,2]),np.isnan(x_nan[:,4])), ~np.isnan(x_nan[:,6])),:]\n",
    "x_nann2n4y6 = x_nan[np.logical_and(np.logical_and(~np.isnan(x_nan[:,2]),~np.isnan(x_nan[:,4])), np.isnan(x_nan[:,6])),:]\n",
    "x_nann2n4n6 = x_nan[np.logical_and(np.logical_and(~np.isnan(x_nan[:,2]),~np.isnan(x_nan[:,4])), ~np.isnan(x_nan[:,6])),:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.shape(x_nany2y4y6), np.shape(x_nany2y4n6), np.shape(x_nany2n4y6), np.shape(x_nany2n4n6)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.shape(x_nann2y4y6), np.shape(x_nann2y4n6), np.shape(x_nann2n4y6), np.shape(x_nann2n4n6)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The combinations we need are then: x_nany2y4y6, x_nany2y4n6, x_nann2y4n6, x_nann2n4y6."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.shape(x_nany2y4y6)[0]+np.shape(x_nany2y4n6)[0]+np.shape(x_nann2y4n6)[0]+np.shape(x_nann2n4y6)[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Actually last three of these nan-versions are so small, that there is probably no point to cluster them anymore. So we can assume them as single clusters for the optimization part. The biggest one should still be splitted a bit more."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "x_nan1 = x_nany2y4y6\n",
    "x_nan2 = x_nany2y4n6\n",
    "x_nan3 = x_nann2y4n6\n",
    "x_nan4 = x_nann2n4y6"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "clust_x_nan1 = np.concatenate((x_nan1[:,:6],x_nan1[:,7:13],x_nan1[:,14:20], x_nan1[:,21:27]),axis=1)\n",
    "norm_clust_nan1 = normalize(clust_x_nan1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time \n",
    "nclust = 10\n",
    "c, xtoc, dist = cluster(norm_clust_nan1, nclust, seedn, verbose=1)\n",
    "weights = np.array([sum(xtoc == i) for i in range(len(c))])\n",
    "opt_x = np.array([x_nan1[xtoc == i].mean(axis=0)\n",
    "                  for i in range(nclust)])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "c_nan2 = x_nan2.mean(axis=0)\n",
    "c_nan3 = x_nan3.mean(axis=0)\n",
    "c_nan4 = x_nan4.mean(axis=0)\n",
    "\n",
    "w_nan2 = np.shape(x_nan2)[0]\n",
    "w_nan3 = np.shape(x_nan3)[0]\n",
    "w_nan4 = np.shape(x_nan4)[0]\n",
    "\n",
    "combined_data = np.concatenate((opt_x,np.array((c_nan2, c_nan3, c_nan4))), axis=0)\n",
    "combined_weights = np.concatenate((weights, np.array((w_nan2, w_nan3, w_nan4))), axis=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [],
   "source": [
    "res_x = np.concatenate((x_nan1, x_nan2, x_nan3, x_nan4), axis=0)\n",
    "res_xtoc = np.concatenate((xtoc, \n",
    "                           np.ones(np.shape(x_nan2)[0])*(nclust), \n",
    "                           np.ones(np.shape(x_nan3)[0])*(nclust+1), \n",
    "                           np.ones(np.shape(x_nan4)[0])*(nclust+2)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "opt = SolverFactory('glpk')\n",
    "\n",
    "prob1, prob2, prob3, prob4 = optimize_all(normalize(combined_data), combined_weights, opt)\n",
    "\n",
    "val1 = model_to_real_values(res_x[:, :7], prob1[0].model, res_xtoc)\n",
    "val2 = model_to_real_values(res_x[:, 7:14], prob2[0].model, res_xtoc)\n",
    "val3 = model_to_real_values(res_x[:, 14:21], prob3[0].model, res_xtoc)\n",
    "val4 = model_to_real_values(res_x[:, 21:], prob4[0].model, res_xtoc)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [],
   "source": [
    "norm_num_x = normalize(x_num)\n",
    "norm_nan_x = normalize(x_nan)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "opt = SolverFactory('glpk')\n",
    "real_nan_revenue, real_nan_carbon, real_nan_deadwood, real_nan_ha = optimize_all(norm_nan_x, np.ones(len(norm_nan_x)), opt)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [],
   "source": [
    "revenue_list = values_to_list(real_nan_revenue[0].model, x_nan[:,:7])\n",
    "carbon_list = values_to_list(real_nan_carbon[0].model, x_nan[:,7:14])\n",
    "deadwood_list = values_to_list(real_nan_deadwood[0].model, x_nan[:,14:21])\n",
    "ha_list = values_to_list(real_nan_ha[0].model, x_nan[:,21:])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "nclust = 100\n",
    "n_nan_opt_revenue, n_nan_opt_carbon, n_nan_opt_deadwood, n_nan_opt_ha = cNopt(x_nan, norm_nan_x, norm_nan_x, opt, nclust, seedn)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('Relative differences to original values (calculated with Nan:s), 50 clusters')\n",
    "print(\"(i) Harvest revenues difference {:.3f}\".format((n_nan_opt_revenue-sum(revenue_list))/sum(revenue_list)))\n",
    "print(\"(ii) Carbon storage {:.3f}\".format((n_nan_opt_carbon-sum(carbon_list))/sum(carbon_list)))\n",
    "print(\"(iii) Deadwood index {:.3f}\".format((n_nan_opt_deadwood-sum(deadwood_list))/sum(deadwood_list)))\n",
    "print(\"(iv) Combined Habitat {:.3f}\".format((n_nan_opt_ha-sum(ha_list))/sum(ha_list)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "nclust = 25\n",
    "n_num_opt_revenue, n_num_opt_carbon, n_num_opt_deadwood, n_num_opt_ha = cNopt(x_num, norm_num_x, norm_num_x, opt, nclust, seedn)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('Relative differences to original values, 50 clusters')\n",
    "print(\"(i) Harvest revenues difference {:.3f}\".format((n_nan_opt_revenue + n_num_opt_revenue-solutions['revenue'])/solutions['revenue']))\n",
    "print(\"(ii) Carbon storage {:.3f}\".format((n_nan_opt_carbon + n_num_opt_carbon-solutions['carbon'])/solutions['carbon']))\n",
    "print(\"(iii) Deadwood index {:.3f}\".format((n_nan_opt_deadwood + n_num_opt_deadwood-solutions['deadwood'])/solutions['deadwood']))\n",
    "print(\"(iv) Combined Habitat {:.3f}\".format((n_nan_opt_ha + n_num_opt_ha-solutions['ha'])/solutions['ha']))"
   ]
  }
 ],
 "metadata": {
  "git": {
   "suppress_outputs": true
  },
  "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.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}