4_Globally_Diversified_Existing_Pools.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "gross-corps",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "\n",
    "from pymoo.algorithms.soo.nonconvex.ga import GA\n",
    "from pymoo.algorithms.moo.unsga3 import UNSGA3\n",
    "from pymoo.factory import get_crossover, get_mutation, get_sampling, get_reference_directions\n",
    "\n",
    "from pymoo.optimize import minimize\n",
    "from pymoo.core.callback import Callback\n",
    "\n",
    "from pooling_functions import MinConcsProblem, MinNumCntrProblems, calc_pools_conc, calc_pool_conc"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "1809a4c3",
   "metadata": {},
   "source": [
    "This notebook relies on pymoo \"0.5.0\". To install it, one needs to run:\n",
    "\n",
    "> pip install pymoo==0.5.0\n",
    "\n",
    "from the command line."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "f79c1efc",
   "metadata": {},
   "outputs": [],
   "source": [
    "class MyCallback(Callback):\n",
    "    def __init__(self) -> None:\n",
    "        super().__init__()\n",
    "        self.n_evals = []\n",
    "        self.hist_F = []\n",
    "\n",
    "    def notify(self, algorithm):\n",
    "        self.n_evals.append(algorithm.evaluator.n_eval)\n",
    "        feas = np.where(algorithm.opt.get(\"feasible\"))[0]\n",
    "        self.hist_F.append(algorithm.opt.get(\"F\")[feas])\n",
    "        self.last_X = algorithm.opt.get(\"X\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "9fc1dfd9",
   "metadata": {},
   "outputs": [],
   "source": [
    "criff_isos = ['BLZ', 'PAN', 'NIC', 'GTM', 'AIA', 'ATG', 'ABW', 'BHS', 'BRB', 'BMU', \n",
    "                'VGB', 'CYM', 'DMA', 'HTI', 'JAM', 'MSR', 'KNA', 'LCA', 'VCT', 'SXM', 'TTO', 'TCA']\n",
    "\n",
    "pcrafi_isos = ['FJI', 'FSM', 'COK', 'KIR', 'MHL', 'NIU', 'PLW', 'PNG', 'WSM', \n",
    "                'SLB', 'TLS', 'TON', 'TUV', 'VUT'] #'NRU'"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "naked-jacob",
   "metadata": {},
   "source": [
    "Load data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "blank-seven",
   "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>ASM</th>\n",
       "      <th>KHM</th>\n",
       "      <th>COK</th>\n",
       "      <th>FJI</th>\n",
       "      <th>PYF</th>\n",
       "      <th>IDN</th>\n",
       "      <th>KIR</th>\n",
       "      <th>LAO</th>\n",
       "      <th>MYS</th>\n",
       "      <th>MHL</th>\n",
       "      <th>...</th>\n",
       "      <th>MWI</th>\n",
       "      <th>MUS</th>\n",
       "      <th>MOZ</th>\n",
       "      <th>SOM</th>\n",
       "      <th>TZA</th>\n",
       "      <th>ZMB</th>\n",
       "      <th>ZWE</th>\n",
       "      <th>BWA</th>\n",
       "      <th>ZAF</th>\n",
       "      <th>SWZ</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0.000000e+00</td>\n",
       "      <td>1.174686e+05</td>\n",
       "      <td>61567.633482</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>2066.621255</td>\n",
       "      <td>4.932370e+09</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.468178e+05</td>\n",
       "      <td>1.363927e+10</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>...</td>\n",
       "      <td>0.0</td>\n",
       "      <td>4.132596e+08</td>\n",
       "      <td>1.240871e+07</td>\n",
       "      <td>1.518563e+07</td>\n",
       "      <td>8.905124e+09</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.297125e+06</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>0.000000e+00</td>\n",
       "      <td>2.135547e+07</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>1817.411519</td>\n",
       "      <td>0.000000e+00</td>\n",
       "      <td>0.0</td>\n",
       "      <td>2.456694e+06</td>\n",
       "      <td>5.648761e+06</td>\n",
       "      <td>81.644624</td>\n",
       "      <td>...</td>\n",
       "      <td>0.0</td>\n",
       "      <td>3.394774e+05</td>\n",
       "      <td>1.486315e+08</td>\n",
       "      <td>0.000000e+00</td>\n",
       "      <td>7.669974e+06</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.000000e+00</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>2.615973e+07</td>\n",
       "      <td>6.706448e+08</td>\n",
       "      <td>2374.730234</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>42.486543</td>\n",
       "      <td>6.417800e+07</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.000000e+00</td>\n",
       "      <td>3.386322e+05</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>...</td>\n",
       "      <td>0.0</td>\n",
       "      <td>6.093447e+07</td>\n",
       "      <td>1.342145e+07</td>\n",
       "      <td>0.000000e+00</td>\n",
       "      <td>0.000000e+00</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.000000e+00</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>6.236010e+07</td>\n",
       "      <td>3.473528e+07</td>\n",
       "      <td>702363.195433</td>\n",
       "      <td>4825.253640</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>0.000000e+00</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.000000e+00</td>\n",
       "      <td>2.667397e+09</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>...</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.005025e+05</td>\n",
       "      <td>1.356082e+04</td>\n",
       "      <td>8.197128e+05</td>\n",
       "      <td>0.000000e+00</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.000000e+00</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>0.000000e+00</td>\n",
       "      <td>2.488902e+07</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>158.037936</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>3.155947e+09</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.000000e+00</td>\n",
       "      <td>0.000000e+00</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>...</td>\n",
       "      <td>0.0</td>\n",
       "      <td>5.339540e+05</td>\n",
       "      <td>0.000000e+00</td>\n",
       "      <td>0.000000e+00</td>\n",
       "      <td>0.000000e+00</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.000000e+00</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>5 rows × 87 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "            ASM           KHM            COK          FJI          PYF  \\\n",
       "0  0.000000e+00  1.174686e+05   61567.633482     0.000000  2066.621255   \n",
       "1  0.000000e+00  2.135547e+07       0.000000     0.000000  1817.411519   \n",
       "2  2.615973e+07  6.706448e+08    2374.730234     0.000000    42.486543   \n",
       "3  6.236010e+07  3.473528e+07  702363.195433  4825.253640     0.000000   \n",
       "4  0.000000e+00  2.488902e+07       0.000000   158.037936     0.000000   \n",
       "\n",
       "            IDN  KIR           LAO           MYS        MHL  ...  MWI  \\\n",
       "0  4.932370e+09  0.0  1.468178e+05  1.363927e+10   0.000000  ...  0.0   \n",
       "1  0.000000e+00  0.0  2.456694e+06  5.648761e+06  81.644624  ...  0.0   \n",
       "2  6.417800e+07  0.0  0.000000e+00  3.386322e+05   0.000000  ...  0.0   \n",
       "3  0.000000e+00  0.0  0.000000e+00  2.667397e+09   0.000000  ...  0.0   \n",
       "4  3.155947e+09  0.0  0.000000e+00  0.000000e+00   0.000000  ...  0.0   \n",
       "\n",
       "            MUS           MOZ           SOM           TZA  ZMB  ZWE  BWA  \\\n",
       "0  4.132596e+08  1.240871e+07  1.518563e+07  8.905124e+09  0.0  0.0  0.0   \n",
       "1  3.394774e+05  1.486315e+08  0.000000e+00  7.669974e+06  0.0  0.0  0.0   \n",
       "2  6.093447e+07  1.342145e+07  0.000000e+00  0.000000e+00  0.0  0.0  0.0   \n",
       "3  1.005025e+05  1.356082e+04  8.197128e+05  0.000000e+00  0.0  0.0  0.0   \n",
       "4  5.339540e+05  0.000000e+00  0.000000e+00  0.000000e+00  0.0  0.0  0.0   \n",
       "\n",
       "            ZAF  SWZ  \n",
       "0  1.297125e+06  0.0  \n",
       "1  0.000000e+00  0.0  \n",
       "2  0.000000e+00  0.0  \n",
       "3  0.000000e+00  0.0  \n",
       "4  0.000000e+00  0.0  \n",
       "\n",
       "[5 rows x 87 columns]"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "datay = pd.read_csv('data/yearly_damages.csv', index_col=0)\n",
    "datay.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "floppy-first",
   "metadata": {},
   "outputs": [],
   "source": [
    "new_cols_ord = pcrafi_isos + criff_isos + [col for col in datay.columns.values if not col in pcrafi_isos + criff_isos]\n",
    "datay = datay[new_cols_ord]\n",
    "\n",
    "x_pool_int_pcrafi = np.ones(len(pcrafi_isos))\n",
    "x_pool_int_criff = np.ones(len(criff_isos))*2\n",
    "\n",
    "x_fixed_pools = np.hstack([x_pool_int_pcrafi, x_pool_int_criff])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "89a0093b",
   "metadata": {},
   "source": [
    "Define alpha:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "4112230d",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.995"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "RT = 200\n",
    "\n",
    "alpha = 1-1/RT\n",
    "alpha"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "3f581b1c",
   "metadata": {},
   "outputs": [],
   "source": [
    "bools = datay.values >= np.nanquantile(datay.values, alpha, axis=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9e1c8d16",
   "metadata": {},
   "source": [
    "Find optimal global extensions for the existing pools and then find the smallest set of countries in the best extension for each pool. The first optimization should be run using parallel processing, for which one would need to uncomment the commented lines below. Also, nm_iter should be increased (10000 was used)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "tough-medication",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "FIND THE OPTIMAL GLOBAL EXTENSIONS FOR PCRAFI AND CCRIF\n",
      "Step 1 - Find the pools with the minimum concentration\n",
      "\n",
      "Step 2.0 - Find the smallest sub-pool within the best\n",
      "PCRAFI globally extended pool which maintains the same concentration\n",
      "\n",
      "Step 2.1 - Find the smallest sub-pool within the best\n",
      "CRIFF globally extended pool which maintains the same concentration\n",
      "\n"
     ]
    }
   ],
   "source": [
    "#from pymoo.core.problem import starmap_parallelized_eval\n",
    "#from multiprocessing.pool import ThreadPool\n",
    "\n",
    "\n",
    "ex_pools_names = ['PCRAFI', 'CRIFF']\n",
    "\n",
    "N = 2\n",
    "\n",
    "n_countries = datay.shape[1]\n",
    "num_cntries_to_expl = n_countries-x_fixed_pools.shape[0]\n",
    "\n",
    "algorithm = UNSGA3(ref_dirs=get_reference_directions(\"energy\", N, 90),\n",
    "                    n_offsprings=num_cntries_to_expl,\n",
    "                    sampling=get_sampling(\"int_random\"),\n",
    "                    crossover=get_crossover(\"int_sbx\"),\n",
    "                    mutation=get_mutation(\"int_pm\"),\n",
    "                    eliminate_duplicates=True)\n",
    "\n",
    "#n_threads = 4*20\n",
    "#pool = ThreadPool(n_threads)\n",
    "problem = MinConcsProblem(num_cntries_to_expl, datay, bools, alpha, calc_pools_conc, N, x_fixed_pools) #,runner=pool.starmap, func_eval=starmap_parallelized_eval)\n",
    "\n",
    "callback = MyCallback()\n",
    "\n",
    "nm_iter = 200 #10000\n",
    "\n",
    "print(\"FIND THE OPTIMAL GLOBAL EXTENSIONS FOR PCRAFI AND CCRIF\")\n",
    "print(\"Step 1 - Find the pools with the minimum concentration\\n\")\n",
    "\n",
    "res = minimize(problem, algorithm, ('n_gen', nm_iter),\n",
    "                   verbose=False, save_history=False, \n",
    "                   callback=callback)\n",
    "\n",
    "#pool.close()\n",
    "\n",
    "best_X = res.X[res.F.argmin(0)]\n",
    "\n",
    "pool_n = 0\n",
    "for sol in best_X:\n",
    "    print(\"Step 2.{} - Find the smallest sub-pool within the best\\n{} globally extended pool which maintains the same concentration\\n\".format(pool_n, ex_pools_names[pool_n]))\n",
    "\n",
    "    opt_sol_pool = sol == pool_n+1\n",
    "\n",
    "    x_ = np.hstack([x_fixed_pools, opt_sol_pool*(pool_n+1)])\n",
    "    cntries = x_ == pool_n+1\n",
    "\n",
    "    datay_reg = datay.iloc[:, cntries>0]\n",
    "    bools_reg = bools[:, cntries>0]\n",
    "\n",
    "    pool_conc = calc_pool_conc(np.repeat(True, datay_reg.shape[1]), datay_reg.values, bools_reg, alpha)[0]\n",
    "\n",
    "    problem = MinNumCntrProblems(np.sum(opt_sol_pool), datay_reg, bools_reg, alpha,\n",
    "                                     calc_pools_conc, 1, pool_conc, x_fixed_pools[x_fixed_pools == pool_n+1]*0+1)\n",
    "\n",
    "    algorithm = GA(pop_size=np.sum(opt_sol_pool),\n",
    "    sampling=get_sampling(\"bin_random\"),\n",
    "    crossover=get_crossover(\"bin_hux\"),\n",
    "    mutation=get_mutation(\"bin_bitflip\"),\n",
    "    eliminate_duplicates=True)\n",
    "\n",
    "    res2 = minimize(problem,\n",
    "                    algorithm,\n",
    "                    ('n_gen', 20), #750\n",
    "                    verbose=False,\n",
    "                    save_history=True,\n",
    "                    pf=np.array(0.0))\n",
    "\n",
    "    pool_n += 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8c63b003",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "3aca3c628c97e81b205090329ad0a7b5843de12fee5fd7de2a2696721c93c40e"
  },
  "kernelspec": {
   "display_name": "Python 3.9.7 ('pymoo_env')",
   "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.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}