figure.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "119ab91e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Standard and GIS Modules\n",
    "import os\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import pandas as pd\n",
    "import time\n",
    "import scipy\n",
    "\n",
    "import warnings\n",
    "import seaborn as sns\n",
    "warnings.filterwarnings(\"ignore\")\n",
    "plt.rcParams.update({\"font.size\":20})"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3712c602",
   "metadata": {},
   "source": [
    "### Real data"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5b5f28a9",
   "metadata": {},
   "source": [
    "##### Compute LOSH"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6eba7823",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "from esda.losh import LOSH\n",
    "import libpysal\n",
    "# ls = LOSH(connectivity=w, inference=\"chi-square\").fit(boston_ds['NOX'])\n",
    "\n",
    "dataset_target = {\n",
    "    \"plants\": \"richness_species_vascular\",\n",
    "    \"meuse\": \"zinc\",\n",
    "    \"atlantic\": \"Rate\",\n",
    "    \"deforestation\": \"deforestation_quantile\",\n",
    "    \"california_housing\": \"median_house_value\",\n",
    "}\n",
    "\n",
    "dataset_losh = {}\n",
    "for d in os.listdir(\"data\"):\n",
    "#     if \"meuse\" in d: #  or \"deforestation\" in d:\n",
    "#         continue\n",
    "#     print(\"------------- \", d)\n",
    "    f = pd.read_csv(os.path.join(\"data\", d))\n",
    "    coords = f[[\"x\", \"y\"]]\n",
    "    \n",
    "    # with KNN:\n",
    "    w = libpysal.weights.KNN(coords, k=20)\n",
    "    \n",
    "    # extract target var\n",
    "    target_var = dataset_target[d.split(\".\")[0]]\n",
    "    \n",
    "    ls = LOSH(connectivity=w, inference=\"chi-square\").fit(f[target_var].values)\n",
    "    \n",
    "    dataset_losh[(d.split(\".\")[0]).replace(\"_\", \" \")] = np.mean(ls.Hi)\n",
    "#     f.drop([\"x\", \"y\", ]], axis=1, inplace=True)\n",
    "#     w_cutoff = (np.max(coords, axis=0) - np.min(coords, axis=0)).sum() / 10\n",
    "#     print(w_cutoff)\n",
    "#     w = get_weights_as_array(np.array(coords), w_cutoff)\n",
    "dataset_losh"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "116756e9",
   "metadata": {},
   "outputs": [],
   "source": [
    "# 100\n",
    "{'plants': 0.9909525898042943,\n",
    " 'california_housing': 0.9199100144573773,\n",
    " 'deforestation': 0.9968442529301278,\n",
    " 'atlantic': 0.9620820900732183}\n",
    "# 10\n",
    "{'plants': 1.0165645863370958,\n",
    " 'california_housing': 0.8750507024465561,\n",
    " 'deforestation': 0.9969507834341067,\n",
    " 'atlantic': 1.0031870804023555}\n",
    "# 20\n",
    "{'plants': 1.0609760745783343,\n",
    " 'california_housing': 0.8773202806724744,\n",
    " 'deforestation': 0.9978335532280674,\n",
    " 'atlantic': 1.004762023375391}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5109aacb",
   "metadata": {},
   "source": [
    "##### Compute sample size"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c75396c3",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "sample_dict = {}\n",
    "feature_dict = {}\n",
    "for out_file in os.listdir(\"data\"):\n",
    "    if \"folds\" in out_file or \"synthetic\" in out_file or out_file[0] == \".\" or \"deprecated\" in out_file:\n",
    "        continue\n",
    "    dataset_name = out_file[:-4].replace(\"_\", \" \")\n",
    "    data = pd.read_csv(os.path.join(\"data\", out_file))\n",
    "    sample_dict[dataset_name] = len(data)\n",
    "    feature_dict[dataset_name] = data.shape[1] - 2\n",
    "sample_dict"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e7de9f24",
   "metadata": {},
   "source": [
    "#### Model recommentation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e38ab3ec",
   "metadata": {},
   "outputs": [],
   "source": [
    "recommendation_dict = {\"california housing\": \"RF\", \"meuse\": \"GWR\", \"plants\": \"GWR / Kriging\", \"atlantic\": \"RF\", \"deforestation\": \"RF\"}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "980e47e2",
   "metadata": {},
   "source": [
    "#### Read results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "54c21e36",
   "metadata": {},
   "outputs": [],
   "source": [
    "res_path = \"outputs/real_feb_23\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a1de3c1d",
   "metadata": {},
   "outputs": [],
   "source": [
    "all_res = []\n",
    "for out_file in os.listdir(res_path):\n",
    "    if \"folds\" in out_file or \"synthetic\" in out_file or out_file[0] == \".\" or \"deprecated\" in out_file:\n",
    "        continue\n",
    "    dataset_name = \" \".join(out_file.split(\"_\")[1:])[:-4]\n",
    "    res = pd.read_csv(os.path.join(res_path, out_file))\n",
    "    res[\"Dataset\"] = dataset_name\n",
    "    raw_data = pd.read_csv(os.path.join(\"data\", out_file[8:]))\n",
    "#     res[\"Samples\"] = len(raw_data)\n",
    "#     res[\"LOSH\"] = round(dataset_losh[out_file[8:-4]], 2)\n",
    "#     res.sort_values(\"Method\", ascending=False, inplace=True)\n",
    "    all_res.append(res)\n",
    "all_res = pd.concat(all_res)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fed1e3b4",
   "metadata": {},
   "outputs": [],
   "source": [
    "all_res.loc[all_res[\"Method\"] == \"linear regression\", \"Method\"] = \"OLS\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1a1f56cb",
   "metadata": {},
   "outputs": [],
   "source": [
    "# remove SAR\n",
    "all_res = all_res[all_res[\"Method\"] != \"SAR\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fc2625c7",
   "metadata": {},
   "outputs": [],
   "source": [
    "# check if all the metrics align in terms of ranking\n",
    "all(all_res.sort_values([\"Dataset\", \"R-Squared\"]).reset_index() == all_res.sort_values([\"Dataset\", \"MAE\"]).reset_index())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ff6f1f14",
   "metadata": {},
   "outputs": [],
   "source": [
    "pivoted = all_res.pivot(index='Dataset', columns='Method', values='RMSE')\n",
    "pivoted.loc[\"california housing\"] = pivoted.loc[\"california housing\"].round()\n",
    "# sorted_columns\n",
    "pivoted[\"LOSH\"] = pd.Series(dataset_losh)\n",
    "pivoted[\"Samples\"] = pd.Series(sample_dict)\n",
    "pivoted[\"k\"] = pd.Series(feature_dict)\n",
    "pivoted[\"Recommended model\"] = pd.Series(recommendation_dict)\n",
    "col_order = ['Dataset', \"Samples\", \"k\", \"LOSH\", \"Recommended model\", 'OLS', 'SLX', 'GWR', 'RF', 'RF (coordinates)', 'spatial RF', 'Kriging']\n",
    "\n",
    "final = pivoted.reset_index().reset_index(drop=True)[col_order] #.drop_index(\"Method\", axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b1e9c0c9",
   "metadata": {},
   "outputs": [],
   "source": [
    "final"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1db0cddf",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(final.to_latex(index=False, float_format=\"%.2f\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "75e44ce3",
   "metadata": {},
   "source": [
    "#### Make table"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "587f6a77",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(all_res.groupby([\"Dataset\", \"Samples\", \"LOSH\", \"Method\"]).mean().round(5).to_latex(float_format=\"%.2f\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ae6a910f",
   "metadata": {},
   "source": [
    "#### Make plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a936cde1",
   "metadata": {},
   "outputs": [],
   "source": [
    "# how big is each real dataset:\n",
    "data_path = \"data/\"\n",
    "data_len = {}\n",
    "for dataset in os.listdir(data_path):\n",
    "    test = pd.read_csv(data_path+dataset)\n",
    "#     print(dataset, len(test))\n",
    "    dataset_name = \" \".join(dataset[:-4].split(\"_\"))\n",
    "    data_len[dataset_name] = f\"{dataset_name}\\n({len(test)})\"\n",
    "all_res[\"Dataset\"] = all_res[\"Dataset\"].map(data_len)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c02ffaf7",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.rcParams.update({\"font.size\":20})\n",
    "plt.figure(figsize=(12,6))\n",
    "sns.barplot(data=all_res.reset_index(), x=\"Dataset\", y=\"R-Squared\", hue=\"Method\")\n",
    "# plt.xlabel(\"Number of samples\")\n",
    "plt.legend(ncol=4, fontsize=15.5)\n",
    "plt.ylim(0., 1.35)\n",
    "plt.tight_layout()\n",
    "plt.xlabel(\"Dataset\", weight=\"bold\")\n",
    "plt.ylabel(\"R-Squared score\", weight=\"bold\")\n",
    "plt.savefig(\"outputs/real_dataset_barplot.pdf\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d312e4c2",
   "metadata": {},
   "outputs": [],
   "source": [
    "results.columns"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "907e595b",
   "metadata": {},
   "source": [
    "## Explanatory plots"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9f20540b",
   "metadata": {},
   "outputs": [],
   "source": [
    "weights = np.array([-0.95, 0.38, 0.66, -0.43, 0.22])\n",
    "\n",
    "nr_data = 10000\n",
    "nr_feats = 5\n",
    "feat_cols = [\"feat_\" + str(i) for i in range(nr_feats)]\n",
    "\n",
    "coords = np.array([[i, j] for i in range(int(np.sqrt(nr_data))) for j in range(int(np.sqrt(nr_data)))])\n",
    "coords = coords / np.max(coords) * 2 - 1\n",
    "\n",
    "spatial_variation = np.zeros((nr_data, nr_feats))\n",
    "for i in range(nr_feats):\n",
    "    spatial_variation[:, i] = 0.5 * (\n",
    "        np.sin(coords[:, 0] * np.pi * 2 + i)\n",
    "        + np.cos(coords[:, 1] * np.pi * 2 + i)\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8d003bf7",
   "metadata": {},
   "outputs": [],
   "source": [
    "cmap = plt.cm.Spectral\n",
    "fig = plt.figure(figsize=(20,4)) # TODO\n",
    "\n",
    "########## Number 1\n",
    "for i in range(5):\n",
    "    ax = fig.add_subplot(1, 5, i+1)\n",
    "#     print(weights[i] + spatial_variation[:, i])\n",
    "    cols = [cmap((val+1) / 2) for val in (weights[i] + 0.5 * spatial_variation[:, i])]\n",
    "    im = ax.scatter(coords[:, 0], coords[:, 1], c=cols)\n",
    "    ax.set_xlabel(rf\"Coefficient $\\beta_{i+1}$\")\n",
    "    ax.set_xticks([])\n",
    "    ax.set_yticks([])\n",
    "    ax.spines['bottom'].set_color('white')\n",
    "    ax.spines['top'].set_color('white') \n",
    "    ax.spines['right'].set_color('white')\n",
    "    ax.spines['left'].set_color('white')\n",
    "\n",
    "import matplotlib as mpl\n",
    "cmap = mpl.cm.viridis\n",
    "bounds = [-1, 2, 5, 7, 12, 15]\n",
    "cmap = mpl.cm.Spectral\n",
    "norm = mpl.colors.Normalize(vmin=-1, vmax=1)\n",
    "\n",
    "fig.subplots_adjust(right=0.88)\n",
    "cbar_ax = fig.add_axes([0.9, 0.15, 0.02, 0.7])\n",
    "cbar = fig.colorbar(mpl.cm.ScalarMappable(norm=norm, cmap=cmap),\n",
    "             cax=cbar_ax)\n",
    "\n",
    "plt.savefig(\"outputs/coefficient_figure_5.png\")\n",
    "plt.show()\n",
    "\n",
    "\n",
    "# ########## Number 2\n",
    "# for i in range(5):\n",
    "#     ax = fig.add_subplot(2, 5, i+1 + 5)\n",
    "# #     print(weights[i] + spatial_variation[:, i])\n",
    "#     cols = [cmap((val+1) / 2) for val in (weights[i] + 0.4 * spatial_variation[:, i])]\n",
    "#     im = ax.scatter(coords[:, 0], coords[:, 1], c=cols)\n",
    "#     ax.set_xlabel(\"Coefficient \"+str(i+1))\n",
    "#     ax.set_xticks([])\n",
    "#     ax.set_yticks([])\n",
    "#     ax.spines['bottom'].set_color('white')\n",
    "#     ax.spines['top'].set_color('white') \n",
    "#     ax.spines['right'].set_color('white')\n",
    "#     ax.spines['left'].set_color('white')\n",
    "\n",
    "# import matplotlib as mpl\n",
    "# cmap = mpl.cm.viridis\n",
    "# bounds = [-1, 2, 5, 7, 12, 15]\n",
    "# cmap = mpl.cm.Spectral\n",
    "# norm = mpl.colors.Normalize(vmin=-1, vmax=1)\n",
    "\n",
    "# fig.subplots_adjust(right=0.825)\n",
    "# cbar_ax = fig.add_axes([0.85, 0.15, 0.05, 0.7])\n",
    "\n",
    "# cbar = fig.colorbar(mpl.cm.ScalarMappable(norm=norm, cmap=cmap),\n",
    "#              cax=cbar_ax)\n",
    "             \n",
    "# plt.tight_layout(w_pad=5, h_pad = 5)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5bffe3f3",
   "metadata": {},
   "source": [
    "### Noise figure\n",
    "\n",
    "To generate the noise plots, run the above code with the three different types of noise, and execute the plot below with each type"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7b4531e5",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# CODE FOR ALL THREE NOISE TYPES\n",
    "noise_level = 0.5\n",
    "\n",
    "nr_data = 10000\n",
    "nr_feats = 5\n",
    "feat_cols = [\"feat_\" + str(i) for i in range(nr_feats)]\n",
    "\n",
    "coords = np.array([[i, j] for i in range(int(np.sqrt(nr_data))) for j in range(int(np.sqrt(nr_data)))])\n",
    "coords = coords / np.max(coords) * 2 - 1\n",
    "synthetic_data = pd.DataFrame(coords, columns=[\"x_coord\", \"y_coord\"])\n",
    "\n",
    "spatial_variation = np.zeros((nr_data, nr_feats))\n",
    "for i in range(nr_feats):\n",
    "    spatial_variation[:, i] = 0.5 * (\n",
    "        np.sin(coords[:, 0] * np.pi * 2 + i)\n",
    "        + np.cos(coords[:, 1] * np.pi * 2 + i)\n",
    "    )\n",
    "\n",
    "    \n",
    "for noise_type in [\"constant\", \"heterogeneous - same\", \"heterogeneous - different\"]:\n",
    "\n",
    "    if noise_type == \"constant\":\n",
    "        noise = np.random.normal(0, noise_level, nr_data)\n",
    "    elif noise_type == \"heterogeneous - different\":\n",
    "        spatial_variation_different = noise_level * (\n",
    "            0.5\n",
    "            * (\n",
    "                synthetic_data[\"x_coord\"].values\n",
    "                + synthetic_data[\"y_coord\"].values\n",
    "            )\n",
    "            + 1\n",
    "        )\n",
    "        noise = np.random.normal(\n",
    "            0,\n",
    "            spatial_variation_different,\n",
    "            len(spatial_variation_different),\n",
    "        )\n",
    "    elif noise_type == \"heterogeneous - same\":\n",
    "        # e.g. high noise level (0.5), spatial variation is from\n",
    "        # sin and cos so it's between -1 and 1, so we make + 1\n",
    "        # so on average we multiply by 1, but varying variance\n",
    "        # between 0.5 * 0 and 0.5 * 2\n",
    "        spatially_dependent_noise = noise_level * (\n",
    "            spatial_variation[:, 0] + 1  # without locality level!\n",
    "        )\n",
    "        noise = np.random.normal(\n",
    "            0, spatially_dependent_noise, nr_data\n",
    "        )\n",
    "    else:\n",
    "        raise RuntimeError(\"Noise must be one of above\")\n",
    "\n",
    "    plt.figure(figsize=(6,4))\n",
    "    plt.scatter(coords[:, 0], coords[:, 1], c=noise, vmin=-3, vmax=3)\n",
    "    plt.colorbar()\n",
    "    # plt.title(f\"Distribution $\\epsilon$ ({noise_type})\", fontsize=15)\n",
    "    plt.axis(\"off\")\n",
    "    plt.tight_layout()\n",
    "    plt.savefig(f\"outputs/noise_{noise_type}.png\")\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bbd0f97b",
   "metadata": {},
   "source": [
    "# Synthetic experiment - load results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f6db9914",
   "metadata": {},
   "outputs": [],
   "source": [
    "def add_discrete_noise(results):\n",
    "    results[\"noise_discrete\"] = pd.NA\n",
    "    results[\"locality_discrete\"] = pd.NA\n",
    "    results.loc[results[\"noise\"] < 0.3,  \"noise_discrete\"] = \"low\"\n",
    "    results.loc[results[\"noise\"] >= 0.3,  \"noise_discrete\"] = \"high\"\n",
    "    results.loc[results[\"locality\"] < 0.3,  \"locality_discrete\"] = \"low\"\n",
    "    results.loc[results[\"locality\"] >= 0.3,  \"locality_discrete\"] = \"high\"\n",
    "    return results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "68ae2751",
   "metadata": {},
   "outputs": [],
   "source": [
    "path = \"outputs/syn_feb_23\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d905af50",
   "metadata": {},
   "outputs": [],
   "source": [
    "# # merge the three files\n",
    "# results = []\n",
    "# for noise in [\"uniformly_distributed\", \"heterogeneous_-_same\", \"heterogeneous_-_different\"]:\n",
    "#     results.append(pd.read_csv(os.path.join(path, \"synthetic_data_results_\"+ noise+\".csv\")))\n",
    "# results = pd.concat(results)\n",
    "# results.loc[results[\"model\"] == \"linear regression\", \"model\"] = \"OLS\"\n",
    "# results = add_discrete_noise(results)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e87187e0",
   "metadata": {},
   "outputs": [],
   "source": [
    "noise_level_range = [0, 0.1, 0.2, 0.3, 0.4, 0.5]\n",
    "locality_range = [0, 0.1, 0.2, 0.3, 0.4, 0.5]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5d6923fb",
   "metadata": {},
   "outputs": [],
   "source": [
    "use_function = \"non-linear 2\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "79662156",
   "metadata": {},
   "outputs": [],
   "source": [
    "# take files from the testing-directory and merge them with the linear results\n",
    "non_linear_dir = \"new_non_linear\"\n",
    "test_dir = \"multiple_nonlinear_tests\"\n",
    "\n",
    "# merge the three files\n",
    "results = []\n",
    "for noise in [\"uniformly_distributed\", \"heterogeneous_-_same\", \"heterogeneous_-_different\", \"train\"]:\n",
    "    results2 = pd.read_csv(os.path.join(path, non_linear_dir, test_dir, \"synthetic_data_results_\"+ noise+\".csv\"))\n",
    "    results1 = pd.read_csv(os.path.join(path, \"synthetic_data_results_\"+ noise+\".csv\"))\n",
    "    # only use non-linear results from file 1 and linear ones from file 2\n",
    "    results2 = results2[results2[\"data mode\"] == use_function]\n",
    "    results1 = results1[results1[\"data mode\"] == \"linear\"]\n",
    "    # concat\n",
    "    results_one = pd.concat((results1, results2))\n",
    "    # post proressing\n",
    "    results_one.loc[results_one[\"model\"] == \"linear regression\", \"model\"] = \"OLS\"\n",
    "    results_one.loc[results_one[\"data mode\"] == use_function, \"data mode\"] = \"non-linear\"\n",
    "    results_one.to_csv(os.path.join(path, non_linear_dir, \"synthetic_data_results_\"+ noise+\".csv\"))\n",
    "    if noise != \"train\":\n",
    "        results.append(results_one)\n",
    "\n",
    "results = pd.concat(results)\n",
    "# general preprocesing steps\n",
    "results = add_discrete_noise(results)\n",
    "print(len(results)//2)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f24e27eb",
   "metadata": {},
   "source": [
    "### Main plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1ea6447f",
   "metadata": {},
   "outputs": [],
   "source": [
    "include_models = ['OLS', 'SLX', 'GWR', 'RF', 'RF (coordinates)', 'spatial RF', 'Kriging']\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "05a3385e",
   "metadata": {},
   "outputs": [],
   "source": [
    "def main_plot(results, nr_data=500, noise_type=\"uniformly distributed\", save_path=\"outputs/main_plot.pdf\", score_col=\"RMSE\"):\n",
    "    include_models = [m for m in include_models if m in results[\"model\"].unique()]\n",
    "    include_function = [\"linear\", \"non-linear\"]\n",
    "    # [model for model in results[\"model\"].unique() if \"geo\" not in model]\n",
    "    nr_models = len(include_models)\n",
    "    fig = plt.figure(figsize=(15, 4))\n",
    "    #fig = plt.figure(figsize=(16, 6.5))\n",
    "    for mode_ind, mode in enumerate(include_function):\n",
    "    #     print(\"----------------\")\n",
    "        for model_ind, model in enumerate(include_models):\n",
    "    #         print(mode, \"data, --> model:\", model)\n",
    "            results_filter = results[\n",
    "                (results[\"data mode\"] == mode) & \n",
    "                (results[\"model\"] == model) & \n",
    "                (results[\"nr_data\"] == nr_data) & \n",
    "                (results[\"noise_type\"] == noise_type)\n",
    "            ]\n",
    "            results_filter.set_index([\"noise\", \"locality\"], inplace=True)\n",
    "            visualize_scores = np.zeros((len(noise_level_range), len(locality_range)))\n",
    "            for i, noise in enumerate(noise_level_range):\n",
    "                for j, locality in enumerate(locality_range):\n",
    "                    score = results_filter.loc[noise, locality][score_col].mean()\n",
    "                    visualize_scores[i, j] = score\n",
    "\n",
    "            ax1 = fig.add_subplot(len(include_function), nr_models+1, ((nr_models+1) * mode_ind) + model_ind+1)\n",
    "            imshow_plot = ax1.imshow(visualize_scores, vmin=0, vmax=0.9)\n",
    "    #         plt.axis(\"off\")\n",
    "#             if model_ind==0:\n",
    "#                 ax1.set_ylabel(\"$\\longleftarrow$ Increasing \\n noise\", fontsize=15)\n",
    "#             ax1.yaxis.set_label_position(\"right\")\n",
    "#             ax1.yaxis.tick_right()\n",
    "            plt.xticks([])\n",
    "            plt.yticks([])\n",
    "#             ax1.set_xlabel(\"$\\longrightarrow$ decreasing \\n stationarity\", fontsize=10)\n",
    "            if model_ind == 0:\n",
    "    #             ax2 = ax1.twinx()\n",
    "    #             ax2.set_ylabel(mode)\n",
    "    #             ax2.yaxis.set_label_position(\"right\")\n",
    "#                 pad = 2\n",
    "                mode_new = \"non-linear\\n(simple)  \" if mode == \"non-linear (simple)\" else mode\n",
    "                ax1.annotate(mode_new, xy=(0, 0.5), xytext=(-50, 0), # ax1.yaxis.labelpad - pad\n",
    "                    xycoords=ax1.yaxis.label, textcoords='offset points',\n",
    "                    size=18, ha='right', va='center', rotation=90, weight=\"bold\")\n",
    "            if mode_ind == 0:\n",
    "                ax1.set_title(model, weight=\"bold\", fontsize=15)\n",
    "    \n",
    "    fig.text(0.5, 0.0, \"$\\longrightarrow$ decreasing stationarity\", ha='center')\n",
    "#     fig.text(0.5, 0.36, \"$\\longrightarrow$ decreasing stationarity\", ha='center')\n",
    "#     fig.text(0.5, 0.7, \"$\\longrightarrow$ decreasing stationarity\", ha='center')\n",
    "    \n",
    "    fig.text(0.06, 0.45, \"$\\longleftarrow$ Increasing noise\", va='center', rotation='vertical')\n",
    "    # make colorbar\n",
    "    # fig.subplots_adjust(right=0.95)\n",
    "    cbar_ax = fig.add_axes([0.88, 0.05, 0.02, 0.9])\n",
    "    fig.colorbar(imshow_plot, cax=cbar_ax, label=score_col)\n",
    "    plt.tight_layout()\n",
    "    if save_path is not None:\n",
    "        plt.savefig(save_path)\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f45f6842",
   "metadata": {},
   "outputs": [],
   "source": [
    "main_plot(results, nr_data=1000, save_path =\"outputs/main_plot.pdf\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "39700ebb",
   "metadata": {},
   "source": [
    "### Barplot - low noise, 0.3 non-stationarity, over samples\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5e76ad28",
   "metadata": {},
   "outputs": [],
   "source": [
    "score_col = \"RMSE\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3ecbf08c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# only look at local models\n",
    "subset = results[\n",
    "    (results[\"model\"] != \"SAR\")\n",
    "     # .isin([\"GWR\", \"RF\", \"RF (coordinates)\", \"spatial RF\", \"Kriging\"])) # \"spatial RF\",\n",
    "]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ac28b3eb",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.rcParams.update({\"font.size\":22})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "441d1439",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# subset.groupby([\"nr_data\", \"data mode\", \"model\", \"noise_discrete\", \"locality_discrete\"]).agg({\"R2 score\": \"mean\"})\n",
    "plt.figure(figsize=(18,6))\n",
    "counter = 1\n",
    "modes = [\"linear\", \"non-linear\"]\n",
    "for mode, save_name in zip(modes, [\"linear\", \"non_linear\"]):\n",
    "    plt.subplot(1, len(modes), counter)\n",
    "    counter += 1\n",
    "    subset_2 = subset[\n",
    "        (subset[\"data mode\"] == mode) &\n",
    "        (subset[\"noise_discrete\"] == \"low\") & \n",
    "        (subset[\"locality_discrete\"] == \"high\") & \n",
    "#         (subset[\"noise\"] == 0.1) & (subset[\"locality\"] == 0.4)&\n",
    "        (subset[\"noise_type\"] == \"uniformly distributed\")\n",
    "    ]\n",
    "    subset_2 = subset_2.groupby([\"nr_data\", \"model\"]).agg({score_col: \"mean\"})\n",
    "\n",
    "    ax = sns.barplot(data=subset_2.reset_index().set_index(\"model\").loc[include_models].reset_index(), x=\"nr_data\", y=score_col, hue=\"model\")\n",
    "    plt.ylim(0, 0.7)\n",
    "    plt.xlabel(\"Number of samples\")\n",
    "#     if mode == \"non-linear (simple)\":\n",
    "#         plt.legend(title=\"Model\", loc=\"lower right\", framealpha=1, ncol=2)\n",
    "#     else:\n",
    "    plt.legend([],[], frameon=False)\n",
    "    plt.title(mode+\" DGP\")\n",
    "    \n",
    "handles, labels = ax.get_legend_handles_labels()\n",
    "plt.tight_layout()\n",
    "plt.figlegend(handles, labels, loc = 'upper center', ncol=7, labelspacing=0., bbox_to_anchor=(0.5,1.09))\n",
    "plt.savefig(f\"outputs/barplot_main.pdf\", bbox_inches=\"tight\")\n",
    "plt.show()\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f89fd4ac",
   "metadata": {},
   "source": [
    "### Noise type analysis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "28fd4877",
   "metadata": {},
   "outputs": [],
   "source": [
    "fontsize=18\n",
    "plt.rcParams.update({\"font.size\":fontsize})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "467e6faa",
   "metadata": {},
   "outputs": [],
   "source": [
    "subset = results[\n",
    "    results[\"model\"].isin([\"OLS\", \"GWR\", \"RF (coordinates)\", \"Kriging\"]) # \"spatial RF\",\n",
    "]\n",
    "subset[\"noise_type\"] = subset[\"noise_type\"].map({\n",
    "    'uniformly distributed':'uniformly distributed noise', 'heterogeneous - same': 'heterogeneous (trigonometric)',\n",
    "       'heterogeneous - different': \"heterogeneous (linear)\" \n",
    "})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1d974285",
   "metadata": {},
   "outputs": [],
   "source": [
    "# fig = plt.figure(figsize=(12, 9))\n",
    "fig = plt.figure(figsize=(13, 4))\n",
    "counter = 1\n",
    "for i, mode in enumerate([\"linear\", \"non-linear\"]): #  \"non-linear (simple)\" # linear\", \n",
    "    for j, model in enumerate([\"GWR\", \"Kriging\"]):\n",
    "        subset2 = subset[\n",
    "                (subset[\"model\"] == model) &\n",
    "                (subset[\"data mode\"] == mode) &\n",
    "        #         (subset[\"noise_discrete\"] == \"low\") & \n",
    "                (subset[\"locality_discrete\"] == \"high\") & \n",
    "        #         (subset[\"noise\"] == 0.3) & \n",
    "#                 (subset[\"locality\"] == 0.4) & \n",
    "        #         (subset[\"noise_type\"] == \"constant\") &\n",
    "                (subset[\"nr_data\"] == 500)\n",
    "        ]\n",
    "#         subset2[\"noise_type\"] = subset2[\"noise_type\"] + \" noise\"\n",
    "        ax = fig.add_subplot(1, 4, counter)\n",
    "        sns.lineplot(ax=ax, data=subset2.reset_index(), x =\"noise\", y=\"RMSE\", hue=\"noise_type\")\n",
    "#         if counter == 1:\n",
    "#             plt.legend(title=\"Noise (spatial distribution)\")# , loc=(1, 1))\n",
    "#         else:\n",
    "        plt.legend([], [], frameon=False)\n",
    "        plt.xlabel(r\"Noise level $\\sigma$\")\n",
    "    \n",
    "        plt.title(f\"{model}\\n{mode} DGP\", fontsize=fontsize)\n",
    "        if counter > 1:\n",
    "            plt.ylabel(\"\")\n",
    "            plt.yticks([], [])\n",
    "#         if j == 0:\n",
    "#             ax.annotate(mode, xy=(0, 0.5), xytext=(-20, 0), # ax1.yaxis.labelpad - pad\n",
    "#                     xycoords=ax.yaxis.label, textcoords='offset points',\n",
    "#                     size='large', ha='right', va='center', rotation=90, weight=\"bold\")\n",
    "#         ax.set_xlabel(\"Noise level ($\\sigma$)\")\n",
    "        ax.set_ylim(0, 0.8)\n",
    "#         if counter in [2, 4]:\n",
    "#             plt.ylabel(\"\")\n",
    "#         if counter in [1, 2]:\n",
    "#             plt.title(model, weight=\"bold\", fontsize=16)\n",
    "        counter += 1\n",
    "\n",
    "handles, labels = ax.get_legend_handles_labels()\n",
    "plt.tight_layout()\n",
    "plt.figlegend(handles, labels, loc = 'upper center', ncol=3, labelspacing=0., bbox_to_anchor=(0.5,1.09))\n",
    "plt.savefig(\"outputs/noise_analysis.pdf\", bbox_inches=\"tight\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5e970481",
   "metadata": {},
   "source": [
    "### Runtime"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "acd774aa",
   "metadata": {},
   "outputs": [],
   "source": [
    "time_plot_data = results.groupby([\"model\", \"nr_data\"])[\"time\"].mean()\n",
    "time_plot_data = time_plot_data.loc[['OLS', 'SLX', 'GWR', 'RF', 'RF (coordinates)', 'spatial RF', 'Kriging']]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bd508470",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(10,4))\n",
    "sns.barplot(data=time_plot_data.reset_index(), x=\"nr_data\", y=\"time\", hue=\"model\")\n",
    "plt.yscale(\"log\")\n",
    "plt.ylabel(\"Runtime (seconds)\")\n",
    "plt.xlabel(\"Number of samples\")\n",
    "plt.legend(ncol=3, loc=\"upper left\")\n",
    "plt.ylim(0, 800)\n",
    "plt.tight_layout()\n",
    "plt.savefig(\"outputs/runtime_plot.pdf\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "de3f6fc3",
   "metadata": {},
   "source": [
    "### Train vs test score"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "831a9c34",
   "metadata": {},
   "outputs": [],
   "source": [
    "train_results = pd.read_csv(os.path.join(path, 'new_non_linear', \"synthetic_data_results_train.csv\"))\n",
    "train_results.loc[train_results[\"model\"] == \"linear regression\", \"model\"] = \"OLS\"\n",
    "test_results = results[results[\"noise_type\"] == \"uniformly distributed\"]# .drop(\"Unnamed: 0\", axis=1)\n",
    "train_results[\"evaluation\"] = \"training error\"\n",
    "test_results[\"evaluation\"] = \"testing error\"\n",
    "print(len(train_results), len(test_results))\n",
    "traintest = pd.concat((train_results, test_results))\n",
    "# traintest[\"R2 score\"] = traintest[\"R2 score\"].clip(-5, 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "05d12297",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(15, 4)) # 10\n",
    "subset = traintest[\n",
    "    (traintest[\"model\"].isin([\"SLX\", \"GWR\"])) # & # \"spatial RF\",\n",
    "    & (traintest[\"nr_data\"] == 1000)\n",
    "].sort_values(\"model\", ascending=False)\n",
    "counter = 1\n",
    "modes_considered = [\"linear\"] # [\"linear\", \"non-linear\"]\n",
    "signal_to_noise_modes = [\n",
    "        \"Weak non-stationarity ($\\lambda$) \\n and weak noise ($\\sigma$)\",\n",
    "        \"Strong non-stationarity ($\\lambda$) \\n and weak noise ($\\sigma$)\",\n",
    "        \"Weak non-stationarity ($\\lambda$) \\n and strong noise ($\\sigma$)\",\n",
    "        \"Strong non-stationarity ($\\lambda$) \\n and strong noise ($\\sigma$)\",\n",
    "]\n",
    "for i, mode in enumerate(modes_considered): #  \"non-linear (simple)\" # linear\", \n",
    "    for j, greatersmaller in enumerate(signal_to_noise_modes):\n",
    "        #\"Low $\\lambda$ and low $\\sigma$\",\n",
    "         #                               \"High $\\lambda$ and low $\\sigma$\", \n",
    "          #                              \"Low $\\lambda$ and high $\\sigma$\"]):\n",
    "        if greatersmaller == \"Weak non-stationarity ($\\lambda$) \\n and weak noise ($\\sigma$)\":\n",
    "            subset2 = subset[\n",
    "                (subset[\"noise\"] < 0.4) & \n",
    "                (subset[\"locality\"] < 0.4)\n",
    "            ]\n",
    "        elif greatersmaller == \"Strong non-stationarity ($\\lambda$) \\n and weak noise ($\\sigma$)\":\n",
    "            subset2 = subset[\n",
    "                (subset[\"noise\"] < 0.4) & \n",
    "                (subset[\"locality\"] >= 0.4)\n",
    "            ]\n",
    "        elif greatersmaller == \"Strong non-stationarity ($\\lambda$) \\n and strong noise ($\\sigma$)\":\n",
    "            subset2 = subset[\n",
    "                (subset[\"noise\"] >= 0.4) & \n",
    "                (subset[\"locality\"] >= 0.4)\n",
    "            ]\n",
    "        elif greatersmaller == \"Weak non-stationarity ($\\lambda$) \\n and strong noise ($\\sigma$)\":\n",
    "            subset2 = subset[\n",
    "                (subset[\"noise\"] >= 0.4) & \n",
    "                (subset[\"locality\"] < 0.4)\n",
    "            ]\n",
    "        subset2 = subset2[subset2[\"data mode\"] == mode]\n",
    "        \n",
    "        if i==0 and j==1: \n",
    "            print(subset2.groupby([\"model\", \"evaluation\"]).agg({\"RMSE\": \"mean\"}))\n",
    "        ax = fig.add_subplot(len(modes_considered), len(signal_to_noise_modes), counter)\n",
    "        sns.barplot(ax=ax, data=subset2, x=\"model\", y=\"RMSE\", hue=\"evaluation\")\n",
    "#         if counter == 2:\n",
    "#             plt.legend(title=\"Evaluation data\", ncol=1, framealpha=1, loc=\"lower center\")\n",
    "#         else:\n",
    "        ymax = 0.59\n",
    "        plt.legend([], [], frameon=False)\n",
    "        if j == 0:\n",
    "            ax.annotate(mode, xy=(0, 0.5), xytext=(-20, 0), # ax1.yaxis.labelpad - pad\n",
    "                    xycoords=ax.yaxis.label, textcoords='offset points',\n",
    "                    size=19, ha='right', va='center', rotation=90, weight=\"bold\")\n",
    "        if i == len(modes_considered)-2:\n",
    "            plt.xticks([], [])\n",
    "        if j > 0:\n",
    "            plt.yticks([], [])\n",
    "        else:\n",
    "            plt.yticks(np.arange(0, ymax, 0.1), np.around(np.arange(0, ymax, 0.1), 1))\n",
    "#         ax.set_xlabel(\"Noise level ($\\sigma$)\")\n",
    "        ax.set_ylim(0, ymax)\n",
    "        ax.set_xlabel(\"\")\n",
    "        if counter in [2, 3, 4]:\n",
    "            plt.ylabel(\"\")\n",
    "        if counter <= len(signal_to_noise_modes):\n",
    "            plt.title(greatersmaller, weight=\"bold\", fontsize=15)\n",
    "        counter += 1\n",
    "plt.tight_layout()\n",
    "\n",
    "handles, labels = ax.get_legend_handles_labels()\n",
    "plt.tight_layout()\n",
    "plt.figlegend(handles, labels, loc = 'upper center', ncol=2, labelspacing=0., bbox_to_anchor=(0.5,1.1))\n",
    "plt.savefig(\"outputs/train_analysis.pdf\", bbox_inches=\"tight\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "592a99ca",
   "metadata": {},
   "outputs": [],
   "source": [
    " 0.198287 / 0.153864, 0.22/0.189239, 0.195645 / 0.149066"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a027c8f8",
   "metadata": {},
   "source": [
    "### Geo rf vs spatial rf (with old results)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3bfd2ea4",
   "metadata": {},
   "outputs": [],
   "source": [
    "res_old = pd.read_csv(\"outputs/synthetic_results_nov_22/synthetic_results.csv\")\n",
    "res_old = res_old[res_old[\"model\"].isin(['geographical RF', 'spatial RF']) & (res_old[\"nr_data\"] < 1000)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cb1177f1",
   "metadata": {},
   "outputs": [],
   "source": [
    "spatial_rf_better = res_old.set_index(\"model\")\n",
    "print(\"Spatial RF has better R2 score in percent cases:\",\n",
    "    sum(spatial_rf_better.loc[\"spatial RF\"][\"R2 score\"].values >= spatial_rf_better.loc[\"geographical RF\"][\"R2 score\"].values) / (len(spatial_rf_better) / 2)\n",
    "     )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8120f832",
   "metadata": {},
   "outputs": [],
   "source": [
    "# # First try: model on x axis and function mode on y axis\n",
    "# plt.figure(figsize=(10, 4))\n",
    "# res_old.loc[res_old[\"R2 score\"] < 0, \"R2 score\"] = 0\n",
    "# sns.barplot(data=res_old, x=\"model\", y=\"R2 score\", hue=\"data mode\")\n",
    "# plt.legend(ncol=3, fontsize=17, loc=\"upper center\")\n",
    "# plt.ylim(0, 1)\n",
    "\n",
    "plt.rcParams.update({\"font.size\":18})\n",
    "plt.figure(figsize=(7, 4.4))\n",
    "res_old.loc[res_old[\"R2 score\"] < 0, \"R2 score\"] = 0\n",
    "sns.barplot(data=res_old, x=\"data mode\", y=\"R2 score\", hue=\"model\")\n",
    "plt.legend(ncol=2, fontsize=18, loc=\"upper center\") # , title=\"Model\")\n",
    "plt.xlabel(\"Function mode\")\n",
    "plt.ylim(0, 1)\n",
    "plt.tight_layout()\n",
    "plt.savefig(\"outputs/geo_vs_spatial_rf.pdf\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8259dff5",
   "metadata": {},
   "source": [
    "## GWR coefficient analysis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "248253fc",
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import pandas as pd\n",
    "import time\n",
    "import scipy\n",
    "import warnings\n",
    "\n",
    "from sklearn.metrics import r2_score\n",
    "from mgwr.gwr import GWR\n",
    "from mgwr.sel_bw import Sel_BW\n",
    "\n",
    "np.random.seed(42)\n",
    "\n",
    "\n",
    "\n",
    "def non_linear_function_simple(feat_arr, weights):\n",
    "    if len(weights.shape) == 1:\n",
    "        weights = np.expand_dims(weights, 0)\n",
    "    function_zoo = [\n",
    "        np.sin,\n",
    "        np.exp,\n",
    "        lambda x: x ** 2,\n",
    "        lambda x: x,\n",
    "        np.cos,\n",
    "        lambda x: np.log(x ** 2),\n",
    "    ]\n",
    "    feature_transformed = np.zeros(feat_arr.shape)\n",
    "    for i in range(feat_arr.shape[1]):\n",
    "        feature_transformed[:, i] = (\n",
    "            function_zoo[i](feat_arr[:, i]) * weights[:, i]\n",
    "        )\n",
    "    return np.sum(feature_transformed, axis=1)\n",
    "\n",
    "nr_feats = 5\n",
    "max_depth = 30\n",
    "noise_type = \"constant\"\n",
    "\n",
    "locality = 0.3\n",
    "\n",
    "weights = np.array([-0.95, 0.38, 0.66, -0.43, 0.22])\n",
    "\n",
    "nr_data = 1000\n",
    "\n",
    "# MAKE MAIN DATA\n",
    "train_cutoff = int(nr_data * 0.9)\n",
    "feat_cols = [\"feat_\" + str(i) for i in range(nr_feats)]\n",
    "synthetic_data = pd.DataFrame(\n",
    "    np.random.rand(nr_data, 2 + nr_feats) * 2 - 1,\n",
    "    columns=[\"x_coord\", \"y_coord\"] + feat_cols,\n",
    ")\n",
    "\n",
    "# simulate spatial variation of features (varying per weight)\n",
    "spatial_variation = np.zeros((nr_data, nr_feats))\n",
    "for i in range(nr_feats):\n",
    "    spatial_variation[:, i] = 0.5 * (\n",
    "        np.sin(synthetic_data[\"x_coord\"].values * np.pi * 2 + i)\n",
    "        + np.cos(synthetic_data[\"y_coord\"].values * np.pi * 2 + i)\n",
    "    )\n",
    "                \n",
    "spatially_dependent_weights = weights + locality * spatial_variation\n",
    "\n",
    "synthetic_data[\"label\"] = non_linear_function_simple(\n",
    "                        synthetic_data[feat_cols].values,\n",
    "                        spatially_dependent_weights,\n",
    "                    )\n",
    "\n",
    "param_arr = [spatially_dependent_weights[:train_cutoff, 0]]\n",
    "for noise_level in [0, 0.3, 0.5]:\n",
    "    noise = np.random.normal(0, noise_level, nr_data)\n",
    "    synthetic_data[\"label\"] = synthetic_data[\"label\"] + noise\n",
    "\n",
    "    train_data, test_data = (\n",
    "                        synthetic_data[:train_cutoff],\n",
    "                        synthetic_data[train_cutoff:],\n",
    "                    )\n",
    "\n",
    "    train_coords = np.array(train_data[[\"x_coord\", \"y_coord\"]])\n",
    "    train_y = np.expand_dims(train_data[\"label\"].values, 1)\n",
    "    train_x = np.array(train_data[feat_cols])\n",
    "    # bandwidth selection\n",
    "    # import pickle\n",
    "    gwr_selector = Sel_BW(\n",
    "        train_coords, train_y, train_x, fixed=True, kernel=\"exponential\"\n",
    "    )\n",
    "    gwr_bw = gwr_selector.search(criterion=\"AICc\")\n",
    "    # create and train model\n",
    "    gwr_model = GWR(\n",
    "        train_coords,\n",
    "        train_y,\n",
    "        train_x,\n",
    "        gwr_bw,\n",
    "        kernel=\"exponential\",\n",
    "        fixed=True,\n",
    "    )\n",
    "    gwr_results = gwr_model.fit()\n",
    "\n",
    "    test_coords = np.array(test_data[[\"x_coord\", \"y_coord\"]])\n",
    "    test_x = np.array(test_data[feat_cols])\n",
    "    # predict\n",
    "    test_pred = gwr_model.predict(\n",
    "        test_coords, test_x, gwr_results.scale, gwr_results.resid_response\n",
    "    ).predictions\n",
    "\n",
    "\n",
    "    score = r2_score(test_pred, test_data[\"label\"])\n",
    "    print(score)\n",
    "    param_arr.append(gwr_results.params[:, 1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "583a60ec",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(16,4))\n",
    "names = [r\"Real $\\beta_1$\", r\"GWR $\\beta_1$ ($\\sigma=0$)\", r\"GWR $\\beta_1$ ($\\sigma=0.3$)\",r\"GWR $\\beta_1$ ($\\sigma=0.5$)\"]\n",
    "for i in range(4):\n",
    "    plt.subplot(1,4,i+1)\n",
    "    plt.scatter(train_data[\"x_coord\"], train_data[\"y_coord\"], c=param_arr[i])\n",
    "    plt.title(names[i], fontsize=18)\n",
    "    plt.axis(\"off\")\n",
    "plt.tight_layout()\n",
    "plt.savefig(\"outputs/gwr_beta_comparison.png\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "45fb53a2",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.9.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}