01_Make_models.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "successful-responsibility",
   "metadata": {},
   "source": [
    "## Code accompanyment to\"Machine learning applied to a modern-Pleistocene petrographic dataset: The global prediction of sand modal composition (GloPrSM) model\"\n",
    "### J. Isaac Johnson, Glenn R. Sharman, Eugene Szymanski, and Xiao Huang\n",
    "### University of Arkansas, Department of Geosciences\n",
    "### Please direct questions and correspondance to gsharman@uark.edu"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "stone-headset",
   "metadata": {},
   "source": [
    "## Step 1: Load sand modal composition data and make random forests models"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "differential-exposure",
   "metadata": {},
   "source": [
    "### Import required modules"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "tender-heather",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "pd.set_option('display.max_columns', None)\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "from sklearn.ensemble import RandomForestRegressor\n",
    "from sklearn.model_selection import KFold, ShuffleSplit, train_test_split\n",
    "from sklearn.metrics import *\n",
    "import time\n",
    "import pickle\n",
    "import pathlib\n",
    "import os"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "hollow-contribution",
   "metadata": {},
   "source": [
    "### Load dependent variable data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "inclusive-moral",
   "metadata": {},
   "outputs": [],
   "source": [
    "rf_data = pd.read_excel('GloPrSM_Input_v1.0.xlsx', engine='openpyxl')\n",
    "rf_data.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "prescription-furniture",
   "metadata": {},
   "source": [
    "### Feature selection"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "pointed-blink",
   "metadata": {},
   "outputs": [],
   "source": [
    "feature_list = ['Relief_AVG', 'Area_sq_km', 'Pre_mm_AVG', 'Tmp_dc_AVG', 'Slope_AVG', 'Lith_PYVAVI',\n",
    "                'Lith_PAPI', 'Lith_VB', 'Lith_PB', 'Lith_EVSC', 'Lith_SMSSSU', 'Lith_MT']\n",
    "features = rf_data[feature_list]\n",
    "features.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "suffering-black",
   "metadata": {},
   "source": [
    "### Feature correlation (optional)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "voluntary-accuracy",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "feature_labels =['Relief','Catchment\\n Area','Precipitation','Temperature','Slope','PY+VA+VI','PA+PI','VB','PB','EV+SC','SM+SS+SU','MT']\n",
    "corr = features.corr()\n",
    "top_corr_features = corr.index\n",
    "fig, ax = plt.subplots(figsize=(12,10))\n",
    "ax.axis('on'), ax.patch.set_edgecolor('black')\n",
    "g=sns.heatmap(features[top_corr_features].corr(), vmin=-1., vmax=1., xticklabels=feature_labels, yticklabels=feature_labels, annot=True, cmap=\"bwr\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "selective-major",
   "metadata": {},
   "source": [
    "### Model training\n",
    "Note: this code can take many hours to run, depending on the number of splits chosen. A large amount of hard drive space is required (~2.5 GB per split used)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "spectacular-formation",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Specify output folder (current directory used by default)\n",
    "base_path = os.getcwd() + '/v1.0'\n",
    "\n",
    "model_path = base_path + '/' + 'models'\n",
    "val_path = base_path + '/' + 'validation'\n",
    "test_path = base_path + '/' + 'test_labels'\n",
    "\n",
    "# Recursively creates the directory and does not raise an exception if the directory already exists\n",
    "pathlib.Path(model_path).mkdir(parents=True, exist_ok=True)\n",
    "pathlib.Path(val_path).mkdir(parents=True, exist_ok=True)\n",
    "pathlib.Path(test_path).mkdir(parents=True, exist_ok=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "technological-villa",
   "metadata": {},
   "outputs": [],
   "source": [
    "labels = ['FQ_QFL_IJ', 'LQ_QFL_IJ', 'QmQch_QmQpQch_IJ', 'QpQch_QmQpQch_IJ', 'FkFp_FpFk_IJ', 'LsLv_LvLsLm_IJ', 'LmLv_LvLsLm_IJ']\n",
    "\n",
    "splits = 1 # Note, 100 splits are used in the article\n",
    "rs = ShuffleSplit(n_splits=splits, test_size=.2, random_state=0)\n",
    "stats = np.zeros(shape=(splits,len(labels)))\n",
    "\n",
    "start = time.time()\n",
    "\n",
    "for j, label in enumerate(labels):\n",
    "#for j, label in enumerate(labels[0:2]): # Just QFL\n",
    "    \n",
    "    val_df = pd.DataFrame()\n",
    "    tst_df = pd.DataFrame()\n",
    "\n",
    "    model = RandomForestRegressor(n_estimators=len(rf_data),random_state=0,max_features='auto')\n",
    "    data = rf_data[rf_data[label].isnull()==False].reset_index()\n",
    "    print(label,len(data))\n",
    "    \n",
    "    i = 0\n",
    "    for train_index, test_index in rs.split(data):\n",
    "        start = time.time()\n",
    "        \n",
    "        train_features = data.iloc[train_index].loc[:, feature_list]\n",
    "        test_features  = data.iloc[test_index].loc[:, feature_list]\n",
    "        train_labels = data.iloc[train_index].loc[:,label]\n",
    "        test_labels  = data.iloc[test_index].loc[:,label]\n",
    "        \n",
    "        model.fit(train_features, train_labels)\n",
    "        prediction = model.predict(test_features)\n",
    "        r2 = r2_score(test_labels, prediction)\n",
    "\n",
    "        stats[i,j] = r2\n",
    "        val_df.loc[:,'{}_valid_{}'.format(label, i)] = prediction\n",
    "        tst_df.loc[:,'{}_label_{}'.format(label, i)] = test_labels\n",
    "\n",
    "        # Export the validation results\n",
    "        val_df.to_csv(val_path+'/'+'{}_validation_rlf.csv'.format(label),index=False)\n",
    "        tst_df.to_csv(test_path+'/'+'{}_label_rlf.csv'.format(label),index=False)\n",
    "        \n",
    "        # Save the model\n",
    "        model_filename = 'model_'+str(i)+'.sav'\n",
    "        model_filepath = model_path+'/'+str(label)\n",
    "        pathlib.Path(model_filepath).mkdir(parents=True, exist_ok=True) # Recursively creates the directory and does not raise an exception if the directory already exists\n",
    "        pickle.dump(model, open(model_filepath+'/'+model_filename, 'wb'))\n",
    "            \n",
    "        print(i, 'R2: {}, {} sec'.format(round(r2,6), round(time.time()-start,1)))    \n",
    "        i += 1\n",
    "    print()\n",
    "\n",
    "r2_df = pd.DataFrame(stats, columns=[x+'_R2' for x in labels])\n",
    "r2_df.to_csv(base_path+'/'+'R2_stats_rlf.csv',index=False)\n",
    "\n",
    "# Save the list of features used in the models, so you know what is going on\n",
    "features = pd.DataFrame()\n",
    "features['Inputs'] = feature_list\n",
    "features.to_csv(base_path+'/'+'feature_list.csv', index=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "referenced-friend",
   "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.9.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}