Feature_Extraction_and_Assemblage_Clustering.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "import os\n",
    "os.environ['CUDA_VISIBLE_DEVICES'] = '1'\n",
    "import pickle\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import skimage.io as io\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "import warnings\n",
    "warnings.filterwarnings(\"ignore\")\n",
    "\n",
    "from keras.applications.resnet50 import preprocess_input\n",
    "from keras.models import Model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### In this part, we conduct the following procedure to make our data be analytic-ready.\n",
    "**Step 1.** For every species, we select out the **representative images**.\n",
    "\n",
    "**Step 2.** For every species representative image, we calculate its **HSV values with regard of different parts** (body, forewing, hindwing, whole)\n",
    "\n",
    "**Step 3.** For every species representative image, we extract its **2048-dimensional features** from the well-trained neural network model\n",
    "\n",
    "**Step 4.** We cluster species based on either the 2-dimensional t-SNE map or 2048D features into **k assemblages through k-Means Clustering**\n",
    "\n",
    "**Step 5.** We use **t-SNE to compress its 2048-dimensional features** into one dimension as the trait value\n",
    "\n",
    "**Step 6.** We quantify the **assemblage-level color diversity** by calculating the average cosine distance among every pair of species in the same assemblage\n",
    "\n",
    "\n",
    "### output files:\n",
    "1. **all_complete_table.csv**: main result for further analysis where a row implies a **species**\n",
    "2. **trait_analysis.csv**: trait value for T-statistics analysis (T stands for trait), where a row implies an **image**\n",
    "3. **cluster_center.csv**: information about assemblage centers where a row implies an assemblage center\n",
    "4. **in-cluser_pairwise_diversity.csv**: result of pair-wise color distance where a row implies a pair of species"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "model_dirname = '/home/put_data/moth/code/cmchang/regression/fullcrop_dp0_newaug-rmhue+old_species_keras_resnet_fold_20181121_4'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# read testing dataset and set the path to obtain every part's mask\n",
    "Xtest = pd.read_csv(os.path.join(model_dirname, 'test.csv'))\n",
    "Xtest['img_rmbg_path'] = Xtest.Number.apply(lambda x: '/home/put_data/moth/data/whole_crop/'+str(x)+'.png')\n",
    "Xtest['img_keep_body_path'] = Xtest.img_rmbg_path.apply(lambda x: x.replace('whole_crop','KEEP_BODY'))\n",
    "Xtest['img_keep_down_path'] = Xtest.img_rmbg_path.apply(lambda x: x.replace('whole_crop','KEEP_DOWN'))\n",
    "Xtest['img_keep_up_path'] = Xtest.img_rmbg_path.apply(lambda x: x.replace('whole_crop','KEEP_UP'))\n",
    "\n",
    "Xtest = Xtest.reset_index()\n",
    "Xtest.drop(columns='index', inplace=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# get the dictionary to look up the average elevation of a species\n",
    "with open(os.path.join('/home/put_data/moth/metadata/1121_Y_mean_dict.pickle'), 'rb') as handle:\n",
    "    Y_dict = pickle.load(handle)\n",
    "Ytest = np.vstack(Xtest['Species'].apply(lambda x: Y_dict[x]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Species</th>\n",
       "      <th>indices</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>Comibaena cassidara</td>\n",
       "      <td>[2022]</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>Trichoplites albimaculosa</td>\n",
       "      <td>[1597]</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>Comibaena delicatior</td>\n",
       "      <td>[100, 1316, 3361]</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>Ourapteryx venusta</td>\n",
       "      <td>[3266]</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>Euproctis kanshireia</td>\n",
       "      <td>[170, 999, 1024, 1379, 1747, 2597, 2838, 2866,...</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                     Species  \\\n",
       "0        Comibaena cassidara   \n",
       "1  Trichoplites albimaculosa   \n",
       "2       Comibaena delicatior   \n",
       "3         Ourapteryx venusta   \n",
       "4       Euproctis kanshireia   \n",
       "\n",
       "                                             indices  \n",
       "0                                             [2022]  \n",
       "1                                             [1597]  \n",
       "2                                  [100, 1316, 3361]  \n",
       "3                                             [3266]  \n",
       "4  [170, 999, 1024, 1379, 1747, 2597, 2838, 2866,...  "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# aggregate the testing data by Species\n",
    "df_species_group = Xtest.groupby('Species').apply(\n",
    "    lambda g: pd.Series({\n",
    "        'indices': g.index.tolist(),\n",
    "    }))\n",
    "df_species_group = df_species_group.sample(frac=1).reset_index()\n",
    "display(df_species_group.head())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Step 1."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# select out the representative image which is the closest to its average elevation\n",
    "sel = list()\n",
    "for k in range(df_species_group.shape[0]):\n",
    "    row = df_species_group.iloc[k]\n",
    "    i = np.argmin(np.abs(np.array(Xtest.Alt[row['indices']]) - Y_dict[row['Species']]))\n",
    "    sel.append(row['indices'][i])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Number</th>\n",
       "      <th>Family</th>\n",
       "      <th>Genus</th>\n",
       "      <th>Species</th>\n",
       "      <th>Latitude</th>\n",
       "      <th>Longititude</th>\n",
       "      <th>Alt</th>\n",
       "      <th>date</th>\n",
       "      <th>img_path</th>\n",
       "      <th>Filename</th>\n",
       "      <th>...</th>\n",
       "      <th>Alt_Q95</th>\n",
       "      <th>Alt_Q100</th>\n",
       "      <th>Alt_mean</th>\n",
       "      <th>Alt_class</th>\n",
       "      <th>srad</th>\n",
       "      <th>tavg</th>\n",
       "      <th>img_rmbg_path</th>\n",
       "      <th>img_keep_body_path</th>\n",
       "      <th>img_keep_down_path</th>\n",
       "      <th>img_keep_up_path</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>V10-20150116-051</td>\n",
       "      <td>Geometridae</td>\n",
       "      <td>Comibaena</td>\n",
       "      <td>Comibaena cassidara</td>\n",
       "      <td>23.665211</td>\n",
       "      <td>120.633595</td>\n",
       "      <td>296.0</td>\n",
       "      <td>20150116.0</td>\n",
       "      <td>/home/dirl/Desktop/moth/data/flickr/V10-201501...</td>\n",
       "      <td>V10-20150116-051.jpg</td>\n",
       "      <td>...</td>\n",
       "      <td>900.00</td>\n",
       "      <td>1000.0</td>\n",
       "      <td>398.777778</td>\n",
       "      <td>0.0</td>\n",
       "      <td>11926.0</td>\n",
       "      <td>15.500000</td>\n",
       "      <td>/home/put_data/moth/data/whole_crop/V10-201501...</td>\n",
       "      <td>/home/put_data/moth/data/KEEP_BODY/V10-2015011...</td>\n",
       "      <td>/home/put_data/moth/data/KEEP_DOWN/V10-2015011...</td>\n",
       "      <td>/home/put_data/moth/data/KEEP_UP/V10-20150116-...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>V12-20131128-025</td>\n",
       "      <td>Geometridae</td>\n",
       "      <td>Trichoplites</td>\n",
       "      <td>Trichoplites albimaculosa</td>\n",
       "      <td>23.618729</td>\n",
       "      <td>120.928362</td>\n",
       "      <td>2300.0</td>\n",
       "      <td>20131128.0</td>\n",
       "      <td>/home/dirl/Desktop/moth/data/flickr/V12-201311...</td>\n",
       "      <td>V12-20131128-025.jpg</td>\n",
       "      <td>...</td>\n",
       "      <td>3002.00</td>\n",
       "      <td>3002.0</td>\n",
       "      <td>2859.548387</td>\n",
       "      <td>5.0</td>\n",
       "      <td>12090.0</td>\n",
       "      <td>9.000000</td>\n",
       "      <td>/home/put_data/moth/data/whole_crop/V12-201311...</td>\n",
       "      <td>/home/put_data/moth/data/KEEP_BODY/V12-2013112...</td>\n",
       "      <td>/home/put_data/moth/data/KEEP_DOWN/V12-2013112...</td>\n",
       "      <td>/home/put_data/moth/data/KEEP_UP/V12-20131128-...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>A34-20140626-064</td>\n",
       "      <td>Geometridae</td>\n",
       "      <td>Comibaena</td>\n",
       "      <td>Comibaena delicatior</td>\n",
       "      <td>24.815020</td>\n",
       "      <td>121.412850</td>\n",
       "      <td>670.0</td>\n",
       "      <td>20140626.0</td>\n",
       "      <td>/home/dirl/Desktop/moth/data/flickr/A34-201406...</td>\n",
       "      <td>A34-20140626-064.jpg</td>\n",
       "      <td>...</td>\n",
       "      <td>3002.00</td>\n",
       "      <td>3002.0</td>\n",
       "      <td>1920.375000</td>\n",
       "      <td>3.0</td>\n",
       "      <td>18227.0</td>\n",
       "      <td>22.500000</td>\n",
       "      <td>/home/put_data/moth/data/whole_crop/A34-201406...</td>\n",
       "      <td>/home/put_data/moth/data/KEEP_BODY/A34-2014062...</td>\n",
       "      <td>/home/put_data/moth/data/KEEP_DOWN/A34-2014062...</td>\n",
       "      <td>/home/put_data/moth/data/KEEP_UP/A34-20140626-...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>V04-20140721-152</td>\n",
       "      <td>Geometridae</td>\n",
       "      <td>Ourapteryx</td>\n",
       "      <td>Ourapteryx venusta</td>\n",
       "      <td>24.257227</td>\n",
       "      <td>121.008787</td>\n",
       "      <td>2287.0</td>\n",
       "      <td>20140721.0</td>\n",
       "      <td>/home/dirl/Desktop/moth/data/flickr/V04-201407...</td>\n",
       "      <td>V04-20140721-152.jpg</td>\n",
       "      <td>...</td>\n",
       "      <td>3002.00</td>\n",
       "      <td>3002.0</td>\n",
       "      <td>2497.777778</td>\n",
       "      <td>4.0</td>\n",
       "      <td>20143.0</td>\n",
       "      <td>16.200001</td>\n",
       "      <td>/home/put_data/moth/data/whole_crop/V04-201407...</td>\n",
       "      <td>/home/put_data/moth/data/KEEP_BODY/V04-2014072...</td>\n",
       "      <td>/home/put_data/moth/data/KEEP_DOWN/V04-2014072...</td>\n",
       "      <td>/home/put_data/moth/data/KEEP_UP/V04-20140721-...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>V02-20140303-079</td>\n",
       "      <td>Erebidae</td>\n",
       "      <td>Euproctis</td>\n",
       "      <td>Euproctis kanshireia</td>\n",
       "      <td>24.051142</td>\n",
       "      <td>120.780668</td>\n",
       "      <td>312.0</td>\n",
       "      <td>20140303.0</td>\n",
       "      <td>/home/dirl/Desktop/moth/data/flickr/V02-201403...</td>\n",
       "      <td>V02-20140303-079.jpg</td>\n",
       "      <td>...</td>\n",
       "      <td>1072.25</td>\n",
       "      <td>1220.0</td>\n",
       "      <td>347.263158</td>\n",
       "      <td>0.0</td>\n",
       "      <td>14704.0</td>\n",
       "      <td>18.200001</td>\n",
       "      <td>/home/put_data/moth/data/whole_crop/V02-201403...</td>\n",
       "      <td>/home/put_data/moth/data/KEEP_BODY/V02-2014030...</td>\n",
       "      <td>/home/put_data/moth/data/KEEP_DOWN/V02-2014030...</td>\n",
       "      <td>/home/put_data/moth/data/KEEP_UP/V02-20140303-...</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>5 rows × 39 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "             Number       Family         Genus                    Species  \\\n",
       "0  V10-20150116-051  Geometridae     Comibaena        Comibaena cassidara   \n",
       "1  V12-20131128-025  Geometridae  Trichoplites  Trichoplites albimaculosa   \n",
       "2  A34-20140626-064  Geometridae     Comibaena       Comibaena delicatior   \n",
       "3  V04-20140721-152  Geometridae    Ourapteryx         Ourapteryx venusta   \n",
       "4  V02-20140303-079     Erebidae     Euproctis       Euproctis kanshireia   \n",
       "\n",
       "    Latitude  Longititude     Alt        date  \\\n",
       "0  23.665211   120.633595   296.0  20150116.0   \n",
       "1  23.618729   120.928362  2300.0  20131128.0   \n",
       "2  24.815020   121.412850   670.0  20140626.0   \n",
       "3  24.257227   121.008787  2287.0  20140721.0   \n",
       "4  24.051142   120.780668   312.0  20140303.0   \n",
       "\n",
       "                                            img_path              Filename  \\\n",
       "0  /home/dirl/Desktop/moth/data/flickr/V10-201501...  V10-20150116-051.jpg   \n",
       "1  /home/dirl/Desktop/moth/data/flickr/V12-201311...  V12-20131128-025.jpg   \n",
       "2  /home/dirl/Desktop/moth/data/flickr/A34-201406...  A34-20140626-064.jpg   \n",
       "3  /home/dirl/Desktop/moth/data/flickr/V04-201407...  V04-20140721-152.jpg   \n",
       "4  /home/dirl/Desktop/moth/data/flickr/V02-201403...  V02-20140303-079.jpg   \n",
       "\n",
       "                         ...                          Alt_Q95  Alt_Q100  \\\n",
       "0                        ...                           900.00    1000.0   \n",
       "1                        ...                          3002.00    3002.0   \n",
       "2                        ...                          3002.00    3002.0   \n",
       "3                        ...                          3002.00    3002.0   \n",
       "4                        ...                          1072.25    1220.0   \n",
       "\n",
       "      Alt_mean  Alt_class     srad       tavg  \\\n",
       "0   398.777778        0.0  11926.0  15.500000   \n",
       "1  2859.548387        5.0  12090.0   9.000000   \n",
       "2  1920.375000        3.0  18227.0  22.500000   \n",
       "3  2497.777778        4.0  20143.0  16.200001   \n",
       "4   347.263158        0.0  14704.0  18.200001   \n",
       "\n",
       "                                       img_rmbg_path  \\\n",
       "0  /home/put_data/moth/data/whole_crop/V10-201501...   \n",
       "1  /home/put_data/moth/data/whole_crop/V12-201311...   \n",
       "2  /home/put_data/moth/data/whole_crop/A34-201406...   \n",
       "3  /home/put_data/moth/data/whole_crop/V04-201407...   \n",
       "4  /home/put_data/moth/data/whole_crop/V02-201403...   \n",
       "\n",
       "                                  img_keep_body_path  \\\n",
       "0  /home/put_data/moth/data/KEEP_BODY/V10-2015011...   \n",
       "1  /home/put_data/moth/data/KEEP_BODY/V12-2013112...   \n",
       "2  /home/put_data/moth/data/KEEP_BODY/A34-2014062...   \n",
       "3  /home/put_data/moth/data/KEEP_BODY/V04-2014072...   \n",
       "4  /home/put_data/moth/data/KEEP_BODY/V02-2014030...   \n",
       "\n",
       "                                  img_keep_down_path  \\\n",
       "0  /home/put_data/moth/data/KEEP_DOWN/V10-2015011...   \n",
       "1  /home/put_data/moth/data/KEEP_DOWN/V12-2013112...   \n",
       "2  /home/put_data/moth/data/KEEP_DOWN/A34-2014062...   \n",
       "3  /home/put_data/moth/data/KEEP_DOWN/V04-2014072...   \n",
       "4  /home/put_data/moth/data/KEEP_DOWN/V02-2014030...   \n",
       "\n",
       "                                    img_keep_up_path  \n",
       "0  /home/put_data/moth/data/KEEP_UP/V10-20150116-...  \n",
       "1  /home/put_data/moth/data/KEEP_UP/V12-20131128-...  \n",
       "2  /home/put_data/moth/data/KEEP_UP/A34-20140626-...  \n",
       "3  /home/put_data/moth/data/KEEP_UP/V04-20140721-...  \n",
       "4  /home/put_data/moth/data/KEEP_UP/V02-20140303-...  \n",
       "\n",
       "[5 rows x 39 columns]"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Xout: DataFrame only contains representative images\n",
    "Xout = Xtest.iloc[sel]\n",
    "Yout = Ytest[sel]\n",
    "\n",
    "Xout = Xout.reset_index()\n",
    "Xout.drop(columns='index', inplace=True)\n",
    "Xout.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Step 2. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1046\r"
     ]
    }
   ],
   "source": [
    "# extract the HSV features for species representatives\n",
    "import skimage.color as color\n",
    "def img_metrics(img):\n",
    "    hsv = color.rgb2hsv(img)\n",
    "    mask = 1.0 - (np.mean(img, axis=2)==255.0) + 0.0\n",
    "    x,y = np.where(mask)\n",
    "    mean_hsv = np.mean(hsv[x,y], axis=0)\n",
    "    std_hsv = np.std(hsv[x,y], axis=0)\n",
    "    return mean_hsv, std_hsv\n",
    "\n",
    "df_reg_list = list()\n",
    "species_list = list()\n",
    "filename_list = list()\n",
    "for k in range(Xout.shape[0]):\n",
    "    print(k, end='\\r')\n",
    "    species = Xout.iloc[k]['Species']\n",
    "    species_list.append(species)\n",
    "    \n",
    "    body_img = io.imread(Xout.iloc[k]['img_keep_body_path'])\n",
    "    mask = 1.0 - (np.mean(body_img, axis=2)==255.0) + 0.0\n",
    "    body_img[:,:,0] = body_img[:,:,0]*mask\n",
    "    body_img[:,:,1] = body_img[:,:,1]*mask\n",
    "    body_img[:,:,2] = body_img[:,:,2]*mask\n",
    "\n",
    "    img = io.imread(Xout.iloc[k]['img_keep_up_path'])\n",
    "\n",
    "    img += body_img\n",
    "\n",
    "    alt = Y_dict[Xout.iloc[k]['Species']]\n",
    "    \n",
    "    mean_hsv, std_hsv = img_metrics(img)\n",
    "    \n",
    "    whole_img = io.imread(Xout.iloc[k]['img_rmbg_path'])\n",
    "    whole_mean_hsv, whole_std_hsv = img_metrics(whole_img)\n",
    "    \n",
    "    res = np.append(whole_mean_hsv[:3], mean_hsv[:3])\n",
    "    res = np.append(res, [alt])\n",
    "    \n",
    "    df_reg_list.append(res)\n",
    "    \n",
    "df_reg_output = pd.DataFrame(data=df_reg_list,\n",
    "                             columns=['h.whole', 's.whole', 'v.whole', \n",
    "                                      'h.body_fore','s.body_fore', 'v.body_fore','alt'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Step 3. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "# extract 2048-dimensional features\n",
    "from keras.models import load_model\n",
    "model = load_model(os.path.join(model_dirname,'model.h5'))\n",
    "features = model.get_layer('global_average_pooling2d_1')\n",
    "extractor = Model(inputs=model.input, outputs=features.output)\n",
    "\n",
    "TestImg = list()\n",
    "for i in range(Xout.shape[0]):\n",
    "    img = io.imread(list(Xout['img_rmbg_path'])[i])\n",
    "    TestImg.append(img)\n",
    "TestImg = np.stack(TestImg)\n",
    "TestInput = preprocess_input(TestImg.astype(float))\n",
    "\n",
    "Fout = extractor.predict(x=TestInput)\n",
    "Yout = np.array([Y_dict[sp] for sp in Xout.Species])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.save(file='Species_Representative_1047x2048.npy', arr=Fout)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1047, 2048)"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Fout.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Step 4."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "# compress 2048-D features to 2-D map for visualization and clustering\n",
    "from sklearn.manifold import TSNE\n",
    "F_embedded = TSNE(n_components=2, perplexity=120).fit_transform(Fout)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "k-means++\t7.64s\t64836.930\t0.054\n"
     ]
    }
   ],
   "source": [
    "from sklearn.cluster import KMeans\n",
    "from sklearn import metrics\n",
    "from time import time\n",
    "\n",
    "def bench_k_means(estimator, name, data):\n",
    "    t0 = time()\n",
    "    estimator.fit(data)\n",
    "    print('%-9s\\t%.2fs\\t%.3f\\t%.3f'\n",
    "          % (name, (time() - t0), estimator.inertia_,\n",
    "             metrics.silhouette_score(data, estimator.labels_,\n",
    "                                      metric='cosine',\n",
    "                                      sample_size=500)))\n",
    "    return estimator\n",
    "for k in [30]:\n",
    "    km = KMeans(init='k-means++', n_clusters=k, n_init=20)\n",
    "    km = bench_k_means(km, name=\"k-means++\", data=Fout)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Counter({0: 49,\n",
       "         1: 86,\n",
       "         2: 65,\n",
       "         3: 13,\n",
       "         4: 1,\n",
       "         5: 5,\n",
       "         6: 24,\n",
       "         7: 1,\n",
       "         8: 6,\n",
       "         9: 88,\n",
       "         10: 93,\n",
       "         11: 3,\n",
       "         12: 5,\n",
       "         13: 92,\n",
       "         14: 40,\n",
       "         15: 36,\n",
       "         16: 61,\n",
       "         17: 49,\n",
       "         18: 70,\n",
       "         19: 1,\n",
       "         20: 13,\n",
       "         21: 41,\n",
       "         22: 10,\n",
       "         23: 33,\n",
       "         24: 56,\n",
       "         25: 47,\n",
       "         26: 2,\n",
       "         27: 39,\n",
       "         28: 15,\n",
       "         29: 3})"
      ]
     },
     "execution_count": 55,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from collections import Counter\n",
    "Counter(km.labels_)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {},
   "outputs": [],
   "source": [
    "Xout['tsne.0'] = F_embedded[:,0]\n",
    "Xout['tsne.1'] = F_embedded[:,1]\n",
    "Xout['km_label'] = km.labels_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {},
   "outputs": [],
   "source": [
    "# representative image information\n",
    "resout = pd.concat([Xout, df_reg_output], axis=1)\n",
    "resout.to_csv(os.path.join(model_dirname, 'all_complete_table.csv'), index=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### - If clustering based on t-SNE maps"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [],
   "source": [
    "# # cluster information\n",
    "# stat = Xout[['km_label','Alt']].groupby('km_label').apply(np.mean)\n",
    "# stat = stat.sort_values('Alt')\n",
    "# stat.columns = ['km_label', 'class_alt']\n",
    "\n",
    "# # center information\n",
    "# centers = km.cluster_centers_\n",
    "# myk = km.cluster_centers_.shape[0]\n",
    "# centx, centy = list(), list()\n",
    "# for i in range(stat.shape[0]):\n",
    "#     centx.append(centers[int(stat.iloc[i]['km_label']),0])\n",
    "#     centy.append(centers[int(stat.iloc[i]['km_label']),1])\n",
    "\n",
    "# # add center information into clustere information\n",
    "# stat['center_x'] = centx\n",
    "# stat['center_y'] = centy\n",
    "# stat['order'] = np.arange(myk)\n",
    "\n",
    "# # output cluster information\n",
    "# stat.to_csv(os.path.join(model_dirname,'cluster_center.csv'), index=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### - If clustering based on 2048D features"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/usr/local/lib/python3.5/dist-packages/ipykernel_launcher.py:5: FutureWarning: 'km_label' is both an index level and a column label.\n",
      "Defaulting to column, but this will raise an ambiguity error in a future version\n",
      "  \"\"\"\n"
     ]
    }
   ],
   "source": [
    "from sklearn.metrics.pairwise import pairwise_distances\n",
    "\n",
    "# cluster information\n",
    "stat = Xout[['km_label','Alt']].groupby('km_label').apply(np.mean)\n",
    "stat = stat.sort_values('km_label')\n",
    "stat.columns = ['km_label', 'class_alt']\n",
    "\n",
    "# center information\n",
    "centers = km.cluster_centers_\n",
    "myk = km.cluster_centers_.shape[0]\n",
    "centx, centy = list(), list()\n",
    "\n",
    "for i in range(myk):\n",
    "    center = centers[i:(i+1),:]\n",
    "    sel = np.where(km.labels_==i)[0]\n",
    "    nearest_species = np.argmin(pairwise_distances(X=center, Y=Fout2[sel], metric='cosine'))\n",
    "    i_nearest_species = sel[nearest_species]\n",
    "    centx.append(F_embedded[i_nearest_species, 0])\n",
    "    centy.append(F_embedded[i_nearest_species, 1])\n",
    "\n",
    "# add center information into clustere information\n",
    "stat['center_x'] = centx\n",
    "stat['center_y'] = centy\n",
    "\n",
    "stat = stat.sort_values('class_alt')\n",
    "# stat.columns = ['km_label', 'class_alt']\n",
    "stat['order'] = np.arange(myk)\n",
    "\n",
    "# output cluster information\n",
    "stat.to_csv(os.path.join(model_dirname,'cluster_center.csv'), index=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Step 5."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "# compress 2048-D features to 1-D trait for functional trait analysis\n",
    "TestImg = list()\n",
    "for i in range(Xtest.shape[0]):\n",
    "    img = io.imread(list(Xtest['img_rmbg_path'])[i])\n",
    "    TestImg.append(img)\n",
    "TestImg = np.stack(TestImg)\n",
    "TestInput = preprocess_input(TestImg.astype(float))\n",
    "Ftest = extractor.predict(x=TestInput)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.manifold import TSNE\n",
    "F_trait = TSNE(n_components=1, perplexity=100).fit_transform(Ftest)\n",
    "F_trait = F_trait - np.min(F_trait)\n",
    "\n",
    "Xtest['trait'] = F_trait[:,0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.save(file='Species_TestingInstance_4249x2048.npy', arr=Ftest)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {},
   "outputs": [],
   "source": [
    "# image trait information table\n",
    "dtrait = pd.merge(Xtest[['Species', 'trait']], resout[['Species','km_label','alt']], how='left', on='Species')\n",
    "dtrait.to_csv(os.path.join(model_dirname, 'trait_analysis.csv'), index=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Step 6."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {},
   "outputs": [],
   "source": [
    "# calculate in-cluster pairwise distance\n",
    "from sklearn.metrics.pairwise import pairwise_distances\n",
    "\n",
    "# just convert the cluster labels to be ordered for better visualization in the next analysis\n",
    "km_label_to_order = dict()\n",
    "order_to_km_label = dict()\n",
    "for i in range(myk):\n",
    "    km_label_to_order[int(stat.iloc[i]['km_label'])] = i\n",
    "    order_to_km_label[i] = int(stat.iloc[i]['km_label'])\n",
    "\n",
    "pair_diversity = np.array([])\n",
    "order = np.array([])\n",
    "for k in range(myk):\n",
    "    \n",
    "    this_km_label = order_to_km_label[k] \n",
    "    sel = np.where(resout.km_label == this_km_label)[0]\n",
    "    if len(sel) == 1:\n",
    "        t = np.array([[0]])\n",
    "        dist_list = np.array([0])\n",
    "    else:\n",
    "        t = pairwise_distances(Fout[sel, :], metric='cosine')\n",
    "        dist_list = np.array([])\n",
    "        for i in range(t.shape[0]):\n",
    "            dist_list = np.append(dist_list,t[i,(i+1):])\n",
    "    pair_diversity = np.append(pair_diversity, dist_list)\n",
    "    order = np.append(order, np.repeat(k, len(dist_list)))\n",
    "    \n",
    "di = pd.DataFrame({'diversity': pair_diversity,\n",
    "              'order': order})\n",
    "di.to_csv(os.path.join(model_dirname, 'in-cluser_pairwise_diversity.csv'), index=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.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}