main_loop.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "cd29baab-2e1d-47fc-9567-819accb25ec8",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "import os                                                 # to set current working directory \n",
    "import numpy as np                                        # arrays and matrix math\n",
    "import pandas as pd                                       # DataFrames\n",
    "import matplotlib.pyplot as plt                           # plotting\n",
    "cmap = plt.cm.inferno                                     # color map\n",
    "import geostatspy.geostats as geostats\n",
    "import geostatspy.GSLIB as GSLIB\n",
    "import math # For n_effective calculations\n",
    "import scipy # For n_effective calculations\n",
    "import time\n",
    "import json\n",
    "\n",
    "from statistics import mean\n",
    "from sklearn.model_selection import train_test_split\n",
    "from sklearn.ensemble import BaggingRegressor\n",
    "from sklearn.tree import DecisionTreeRegressor\n",
    "from sklearn.metrics import mean_squared_error\n",
    "from sklearn.model_selection import GridSearchCV\n",
    "from sklearn.preprocessing import StandardScaler\n",
    "from IPython.display import clear_output\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "2b35cea8-1f70-4ec0-9b55-96c6b713e0dd",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1/4\n"
     ]
    }
   ],
   "source": [
    "with open('.workflow_ignore.txt') as f:\n",
    "    ignore_list = f.read().splitlines()\n",
    "\n",
    "filename_list = [f.path.split('/')[-1] for f in os.scandir('dataset/') if(f.is_dir() and f.path.split('/')[-1] not in ignore_list)] # Getting files\n",
    "\n",
    "for filenum, filename in enumerate(filename_list):\n",
    "    #### READ DATA AND PARSE\n",
    "    df = pd.read_csv('dataset/' + filename + '/' + filename + '.csv')\n",
    "    \n",
    "    # We are sampling in every 500 meter for a realistic scenario\n",
    "    # The data where every 100x100 isotropic cell is stored is going\n",
    "    # to be our truth model for the results.\n",
    "    df = df.iloc[::5] \n",
    "\n",
    "    predictors = df.iloc[:,2:-1].columns.values\n",
    "    response = df.columns[-1]\n",
    "    predictors_and_response = df.iloc[:,2:].columns.values\n",
    "\n",
    "    num_predictors = predictors.shape[0] # NO NEED, but defined just in case\n",
    "    num_response = 1 # NO NEED, but defined just in case\n",
    "    num_predictors_and_response = predictors_and_response.shape[0] # Used for loops\n",
    "\n",
    "    df_train, df_test = train_test_split(df, test_size=0.5, random_state=73073)\n",
    "    X_train = df_train.loc[:, predictors].values\n",
    "    y_train = df_train.loc[:, response].values\n",
    "    X_test = df_test.loc[:, predictors].values\n",
    "    y_test = df_test.loc[:, response].values\n",
    "    #### READ REQUIRED PARAMS\n",
    "    with open('dataset/' + filename + '/' + filename + '_var_params_and_n_eff.json', 'r') as f:\n",
    "        data = json.load(f)\n",
    "\n",
    "    #### ASSIGN VARIABLES FROM PARAMS\n",
    "    n_eff_calculated = data['n_eff_calculated']\n",
    "    n_standard = y_train.shape[0]\n",
    "\n",
    "    #### HYPERPARAMETER TUNING:\n",
    "\n",
    "    L = 20\n",
    "    max_num_leaf_nodes = [i for i in range(2,n_standard//3,20)]\n",
    "    vals = list(range(y_train.shape[0])) # We will draw from the indexes and then we will slice the dataset with these indexes. We do this because data is multivariate\n",
    "    n_eff_list = sorted(set([n_eff_calculated] + [i for i in range(10, df_train.shape[0]+1, 10)] + [df_train.shape[0]]))\n",
    "    n_eff_params = {}\n",
    "    n_eff_train_params = {}\n",
    "    hyperparameter_error_spatial = []\n",
    "    hyperparameter_error_standard = []\n",
    "\n",
    "    hyperparameter_mse_test_dict = {}\n",
    "    hyperparameter_mse_train_dict = {}\n",
    "\n",
    "    for n_eff in n_eff_list:\n",
    "        hyperparameter_mse_test_dict[n_eff] = {}\n",
    "        hyperparameter_mse_train_dict[n_eff] = {}\n",
    "        mse_init = 9999\n",
    "        mse_train_init = 9999\n",
    "        node_init = 0\n",
    "        node_train_init = 0\n",
    "        for node in max_num_leaf_nodes:\n",
    "            regressor = DecisionTreeRegressor(max_leaf_nodes = node,random_state = 73073) # One decision tree\n",
    "            y_pred = np.zeros(y_test.shape)\n",
    "            y_pred_train = np.zeros(y_train.shape)\n",
    "\n",
    "            for i in range(L):\n",
    "                instance = np.random.choice(vals, size=n_eff, replace=True) # Size is equal to n_effective\n",
    "                # Fitting standard bootstrap and storing the predictions for one realization\n",
    "                dt = regressor.fit(X_train[instance], y_train[instance])\n",
    "                y_pred = y_pred + dt.predict(X_test)\n",
    "                y_pred_train = y_pred_train + dt.predict(X_train)\n",
    "\n",
    "            y_pred /= (L+1)\n",
    "            y_pred_train /= (L+1)\n",
    "            mse = mean_squared_error(y_test, y_pred)\n",
    "            mse_train = mean_squared_error(y_train, y_pred_train)\n",
    "            if(mse<mse_init):\n",
    "                mse_init = mse\n",
    "                node_init = node\n",
    "            hyperparameter_mse_test_dict[n_eff][node] = mse\n",
    "\n",
    "            if(mse_train<mse_train_init):\n",
    "                mse_train_init = mse_train\n",
    "                node_train_init = node\n",
    "            hyperparameter_mse_train_dict[n_eff][node] = mse_train\n",
    "\n",
    "        n_eff_params[n_eff] = node_init\n",
    "        n_eff_train_params[n_eff] = node_train_init\n",
    "\n",
    "    # Obtain hyperparameter tuning with training data\n",
    "    json_string = json.dumps(hyperparameter_mse_train_dict)\n",
    "    with open('dataset/' + filename + '/' + filename + '_hyperparameter_mse_train_dict.json', 'w') as f:\n",
    "        f.write(json_string)\n",
    "\n",
    "    # Obtain hyperparameter tuning with testing data\n",
    "    json_string = json.dumps(hyperparameter_mse_test_dict)\n",
    "    with open('dataset/' + filename + '/' + filename + '_hyperparameter_mse_test_dict.json', 'w') as f:\n",
    "        f.write(json_string)\n",
    "\n",
    "    # Obtain tuned hyperparameter for each n\n",
    "    json_string = json.dumps(n_eff_params)\n",
    "    with open('dataset/' + filename + '/' + filename + '_hyperparameter_for_each_n.json', 'w') as f:\n",
    "        f.write(json_string)\n",
    "\n",
    "    vals = list(range(y_train.shape[0])) # We will draw from the indexes and then we will slice the dataset with these indexes. We do this because data is multivariate\n",
    "    len_vals = len(vals) # Total size of train data for easiness of syntax\n",
    "    L = 100 # Number of realizations for each n_effective\n",
    "    L_2 = 10 # Bootstrap of Bootstrap number\n",
    "    n_bins = 50 # Number of bins for histograms\n",
    "\n",
    "    n_eff_broadcast = np.array([n_eff_list for i in range(L_2)])\n",
    "\n",
    "    # Creating empty lists to store bootstrapped samples\n",
    "\n",
    "    boots_arr = [[[None for i in range(L)] for j in range(len(n_eff_list))] for k in range(L_2)]\n",
    "\n",
    "    # Creating empty lists to store Decision Tree regression results\n",
    "\n",
    "    reals = [[[None for i in range(L)] for j in range(len(n_eff_list))] for k in range(L_2)]\n",
    "\n",
    "    # Creating empty lists to store Decision Tree regression Training Fits\n",
    "\n",
    "    trains = [[[None for i in range(L)] for j in range(len(n_eff_list))] for k in range(L_2)]\n",
    "\n",
    "    # Creating empty lists to store the aggregate (bagged) of L realizations for each n_effective\n",
    "\n",
    "    bags = [[None for i in range(len(n_eff_list))] for j in range(L_2)]\n",
    "\n",
    "    # Creating empty lists to store the aggregate (bagged) of L realizations for each n_effective (TRAINING DATA)\n",
    "\n",
    "    bags_trains = [[None for i in range(len(n_eff_list))] for j in range(L_2)]\n",
    "\n",
    "    # Creating empty lists to store MSE for each bagged result.\n",
    "\n",
    "    mse_test = [[None for i in range(len(n_eff_list))] for j in range(L_2)]\n",
    "\n",
    "    # Creating empty lists to store MSE for each bagged result. (TRAINING DATA)\n",
    "\n",
    "    mse_train = [[None for i in range(len(n_eff_list))] for j in range(L_2)]\n",
    "\n",
    "    # Main Loop\n",
    "\n",
    "    # For each n_effective (number of draws or n_draw here) in n_effective list\n",
    "    for n_idx, n_draw in enumerate(n_eff_list):\n",
    "        node = n_eff_params[n_draw]\n",
    "        regressor = DecisionTreeRegressor(max_leaf_nodes=node, random_state = 73073) # One decision tree\n",
    "        for j in range(L_2):\n",
    "            # We will have L (100) realizations. For each of them, do the following:\n",
    "            for i in range(L):\n",
    "                # Sample with replacement (bootstrap). Note that we sample INDEXES, not the values. Then, using these indexes, we will draw from X_train\n",
    "                spt_idxs = np.random.choice(vals, size=n_draw, replace=True) # Size is equal to n_effective\n",
    "\n",
    "                #Storing the Indexes as bootstrap realizations so that we can access them later if we need\n",
    "                boots_arr[j][n_idx][i] = spt_idxs\n",
    "\n",
    "                # Fitting spatial bootstrap and storing the predictions for one realization\n",
    "                dt = regressor.fit(X_train[spt_idxs], y_train[spt_idxs])\n",
    "                trains[j][n_idx][i] = dt.predict(X_train)\n",
    "                reals[j][n_idx][i] = dt.predict(X_test)\n",
    "\n",
    "            # After L realizations made, we aggregate the predictions (bagging) and store them in respective lists    \n",
    "            bags[j][n_idx] = np.mean(reals[j][n_idx], axis=0)\n",
    "            bags_trains[j][n_idx] = np.mean(trains[j][n_idx], axis=0)\n",
    "\n",
    "            # Last, we calculate MSE for each aggregation.\n",
    "            mse_test[j][n_idx] = mean_squared_error(y_test, bags[j][n_idx])\n",
    "            mse_train[j][n_idx] = mean_squared_error(y_train, bags_trains[j][n_idx])\n",
    "\n",
    "\n",
    "    boots_arr = np.asarray(boots_arr, dtype=object)\n",
    "    reals = np.asarray(reals)\n",
    "    bags = np.asarray(bags)\n",
    "    bags_trains = np.asarray(bags_trains)\n",
    "    mse_test = np.asarray(mse_test)\n",
    "    mse_train = np.asarray(mse_train)\n",
    "\n",
    "    mse_train_dict = {}\n",
    "    mse_test_dict = {}\n",
    "\n",
    "    for idx in range(len(n_eff_list)):\n",
    "        mse_train_dict.update({n_eff_list[idx]: mse_train[:,idx].tolist()})\n",
    "        mse_test_dict.update({n_eff_list[idx]: mse_test[:,idx].tolist()})\n",
    "\n",
    "    # Obtain mse for training\n",
    "    json_string = json.dumps(mse_train_dict)\n",
    "    with open('dataset/' + filename + '/' + filename + '_mse_train.json', 'w') as f:\n",
    "        f.write(json_string)\n",
    "\n",
    "    # Obtain mse for training\n",
    "    json_string = json.dumps(mse_test_dict)\n",
    "    with open('dataset/' + filename + '/' + filename + '_mse_test.json', 'w') as f:\n",
    "        f.write(json_string)\n",
    "\n",
    "    mse_ratio_dict = {filename + '/' + filename: mean(mse_test_dict[n_standard]) / mean(mse_test_dict[n_eff_calculated])}\n",
    "\n",
    "    # Obtain mse ratios\n",
    "    json_string = json.dumps(mse_ratio_dict)\n",
    "    with open('dataset/' + filename + '/' + filename + '_mse_ratio.json', 'w') as f:\n",
    "        f.write(json_string)\n",
    "    print(str(filenum+1) + '/' + str(len(filename_list)))\n",
    "    clear_output(wait=True)\n",
    "    \n",
    "    break"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.10.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}