examples/understanding_array_thresholding.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "a3b87237-8575-4798-a392-278aebf77801",
   "metadata": {},
   "source": [
    "# Deep dive into thresholding continuous fields to make categorical fields\n",
    "\n",
    "In this notebook we will explore the categorization process that is in pyEMU.  Note this process was inspired by Todaro and other (2023) \"Experimental sandbox tracer tests to characterize a two-facies aquifer via an ensemble smoother\" [https://doi.org/10.1007/s10040-023-02662-1](https://doi.org/10.1007/s10040-023-02662-1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "37e0fb1f-51b0-442c-88e1-2812b54ad902",
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import shutil\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import pyemu"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "62977378-37d4-4a2a-8801-b0ade16e7440",
   "metadata": {},
   "source": [
    "First, lets generate just a single multivariate gaussian field.  We will use this as the \"original array\" thru out the rest of this notebook:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "72b9a89c-84d8-4109-8dc7-d18714817191",
   "metadata": {},
   "outputs": [],
   "source": [
    "nrow = ncol = 50\n",
    "delx = np.ones(ncol)\n",
    "dely = np.ones(nrow)\n",
    "v = pyemu.geostats.ExpVario(contribution=1.0,a=500)\n",
    "gs = pyemu.geostats.GeoStruct(variograms=v)\n",
    "ss = pyemu.geostats.SpecSim2d(delx=delx,dely=dely,geostruct=gs)\n",
    "np.random.seed(122341)\n",
    "org_arr = ss.draw_arrays(1,mean_value=10)[0,:,:]\n",
    "assert org_arr.min() > 0.0\n",
    "cb = plt.imshow(org_arr)\n",
    "_ = plt.colorbar(cb)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f7094e59-f763-41eb-988d-da9193249ebe",
   "metadata": {},
   "source": [
    "Now let's setup a workspace"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "322ca92f-06ce-4492-a43a-c864e75b6df7",
   "metadata": {},
   "outputs": [],
   "source": [
    "ws = \"temp_thresh\"\n",
    "if os.path.exists(ws):\n",
    "    shutil.rmtree(ws)\n",
    "os.makedirs(ws)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "faabba44-0fa6-49d2-9213-5b5c261c02dc",
   "metadata": {},
   "source": [
    "And save the original array in that workspace"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "209f5fbc-4f2d-48c5-b3d6-7495c5cac86c",
   "metadata": {},
   "outputs": [],
   "source": [
    "orgarr_file = os.path.join(ws,\"orgarr.dat\")\n",
    "np.savetxt(orgarr_file,org_arr)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "14a06478-8c58-4fe0-9ce1-b3a2808315f2",
   "metadata": {},
   "source": [
    "The categorization process in pyEMU currently only supports 2 categories/facies.  So we need to define a `dict` that contains the proportions of each category and the initial value used to fill each category.  For example purposes, we will use the extrema of the original array as the fill values, so that we will end up with a categorical array that has only two unique values: the min and max of the original array:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a86bcae0-e99b-4abd-8a70-3574f6a8e3a9",
   "metadata": {},
   "outputs": [],
   "source": [
    "cat_dict = {1:[0.95,org_arr.min()],2:[0.05,org_arr.max()]}\n",
    "thresharr_file,threshcsv_file = pyemu.helpers.setup_threshold_pars(orgarr_file,cat_dict=cat_dict,\n",
    "                                                         testing_workspace=ws)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0cf681cc-cd3f-48ec-9b01-a0c78f34ac77",
   "metadata": {},
   "source": [
    "Now lets see what was created as part of the setup process:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "55aefaae-cad3-4b3c-a026-271de70d3737",
   "metadata": {},
   "outputs": [],
   "source": [
    "os.listdir(ws)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8d46d36f-b377-4c53-b1b2-82eff652a2d0",
   "metadata": {},
   "source": [
    "Notice the what we have are files with the original array name and some suffix.  Let's check them out:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "47047de9-27f0-4fad-aeb4-895129c7acca",
   "metadata": {},
   "outputs": [],
   "source": [
    "df = pd.read_csv(threshcsv_file)\n",
    "df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "323471e6-728d-42e6-8b4c-49e6b7e30e72",
   "metadata": {},
   "outputs": [],
   "source": [
    "catarr_file = orgarr_file+\".threshcat.dat\"\n",
    "print(catarr_file)\n",
    "assert os.path.exists(catarr_file)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1defbda2-e13b-4a8a-98bd-a038f4718cf9",
   "metadata": {},
   "outputs": [],
   "source": [
    "thresharr_file = orgarr_file+\".thresharr.dat\"\n",
    "print(thresharr_file)\n",
    "assert os.path.exists(thresharr_file)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ffc78187-afee-4665-94da-895da6b1f6cc",
   "metadata": {},
   "outputs": [],
   "source": [
    "threshcsv_file_results = orgarr_file+'.threshprops_results.csv'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9e80d025-891f-4ca7-943b-ddb708286707",
   "metadata": {},
   "outputs": [],
   "source": [
    "def load_and_plot(save_name=None):\n",
    "    cat_arr = np.loadtxt(catarr_file)\n",
    "    new_arr = np.loadtxt(orgarr_file)\n",
    "    thresh_arr = np.loadtxt(thresharr_file)\n",
    "    thresh_arr = (thresh_arr-thresh_arr.min())/thresh_arr.max()\n",
    "    ddf = pd.read_csv(threshcsv_file_results)\n",
    "    cat1_prop = ddf.loc[0,\"proportion\"]/ ddf.loc[:,\"proportion\"].sum()\n",
    "    cat2_prop = ddf.loc[1,\"proportion\"]/ ddf.loc[:,\"proportion\"].sum()\n",
    "    cat1_thresh = ddf.loc[0,\"threshold\"]\n",
    "    cat2_thresh = ddf.loc[1,\"threshold\"]\n",
    "    fig,axes = plt.subplots(1,3,figsize=(10,2.5))\n",
    "    cb = axes[0].imshow(org_arr)\n",
    "    plt.colorbar(cb,ax=axes[0])\n",
    "    #cb = axes[2].imshow(cat_arr)#,vmin=org_arr.min(),vmax=org_arr.max())\n",
    "    #plt.colorbar(cb,ax=axes[2])\n",
    "    cb = axes[2].imshow(new_arr,vmin=org_arr.min(),vmax=org_arr.max())\n",
    "    plt.colorbar(cb,ax=axes[2])\n",
    "    cb = axes[1].imshow(thresh_arr)#,vmin=org_arr.min(),vmax=org_arr.max())\n",
    "    plt.colorbar(cb,ax=axes[1])\n",
    "    axes[1].contour(thresh_arr,levels=[cat1_thresh,cat2_thresh],colors=['w',\"w\"])\n",
    "    axes[0].set_title(\"original array\\n\",loc=\"left\",fontsize=8)\n",
    "    axes[1].set_title(\"thresholding array\\ncat 1 threshold:{0:3.4f}\".\\\n",
    "                      format(cat1_thresh),loc=\"left\",fontsize=8)\n",
    "    #axes[2].set_title(\"categorized array\\ncat 1 proportion: {0:3.4f}\".\\\n",
    "    #                  format(cat1_prop)\\\n",
    "    #                  ,loc=\"left\",fontsize=10)\n",
    "    axes[2].set_title(\"new array\\n{1:2.2f}% max, {0:2.2f}% min\".format(cat1_prop*100,cat2_prop*100)\n",
    "                      ,loc=\"left\",fontsize=8)\n",
    "    plt.tight_layout()\n",
    "    for ax in axes:\n",
    "        ax.set_yticks([])\n",
    "        ax.set_xticks([])\n",
    "    plt.tight_layout()\n",
    "    if save_name is not None:\n",
    "        plt.savefig(save_name)\n",
    "    else:    \n",
    "        plt.show()\n",
    "    plt.close(fig)\n",
    "    return cat_arr,new_arr,ddf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d3e6858a-9892-4afc-a5d9-0840ca9ad99b",
   "metadata": {},
   "outputs": [],
   "source": [
    "cat_arr,new_arr,newnew_df = load_and_plot()\n",
    "newnew_df"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "968ee1a0-6d37-4f14-bcd3-c00a4c9e87f1",
   "metadata": {},
   "source": [
    "So there it is! The original array for reference, the \"thresholding array\" (which is just a scaled and normed version of the original array) and the resulting \"new array\".\n",
    "\n",
    "\n",
    "Now let's experiment - feel free to change the quantities in `new_df`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ee5c2cc6-d9f1-4546-a55d-24d812a53b2e",
   "metadata": {},
   "outputs": [],
   "source": [
    "new_df = df.copy()\n",
    "new_df.loc[0,\"threshproportion\"] = .25\n",
    "new_df.to_csv(threshcsv_file)\n",
    "pyemu.helpers.apply_threshold_pars(threshcsv_file)\n",
    "_,_,newnew_df = load_and_plot()\n",
    "newnew_df.iloc[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1c451fea-bd3a-46b2-b7b8-dddd30947fb8",
   "metadata": {},
   "source": [
    "Now lets sweep over a range of category 1 proportions and make some figs:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ffdadd80-4b39-4b03-b91a-280f4273aaff",
   "metadata": {},
   "outputs": [],
   "source": [
    "cat1_props = np.linspace(0.01,0.99,100)\n",
    "cat1_props"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "330c370a-e286-4cda-aeb3-ec35aa7c1efe",
   "metadata": {},
   "outputs": [],
   "source": [
    "for i,prop in enumerate(cat1_props):\n",
    "    new_df = df.copy()\n",
    "    new_df.loc[0,\"threshproportion\"] = prop\n",
    "    new_df.to_csv(threshcsv_file)\n",
    "    pyemu.helpers.apply_threshold_pars(threshcsv_file)\n",
    "    save_name = os.path.join(ws,\"fig_{0:04d}.png\".format(i))\n",
    "    _,_,newnew_df = load_and_plot(save_name=save_name)\n",
    "    print(i,\" \",end='')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0b48cf4e-7f48-4b47-9fe5-9940f855cb2f",
   "metadata": {},
   "source": [
    "And if you have `ffmpeg` installed, we can make an sweet-as animated gif:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0c1144ba-7410-43df-896d-2525c093c1e2",
   "metadata": {},
   "outputs": [],
   "source": [
    "fps = 15\n",
    "pyemu.os_utils.run(\"ffmpeg -i fig_{0:04d}.png -vf palettegen=256 palette.png\".format(int(len(cat1_props)/2)),cwd=ws)\n",
    "pyemu.os_utils.run(\"ffmpeg -r {0} -y -s 1920X1080 -i fig_%04d.png -i palette.png -filter_complex \\\"scale=720:-1:flags=lanczos[x];[x][1:v]paletteuse\\\" fancy.gif\".format(fps),\n",
    "        cwd=ws)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7bfa7e9d-c49b-4ad2-b223-7d161598de5e",
   "metadata": {},
   "source": [
    "![SegmentLocal](temp_thresh/fancy.gif \"segment1\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "59eb4b05-9d3f-4baa-9193-169d7385b355",
   "metadata": {},
   "source": [
    "So how does this work within the PEST world?  Well we can treat the thresholding array as an array we want to parameterize (maybe with pilot points?) as well as parameterizing the fill values and proportions in the \"threshprops.csv\" file.  This will let us manipulate the shape of the resulting categorical array that forward model will use as an input. In turn, this yield variability in the simulated response of the system.  And away we go!"
   ]
  }
 ],
 "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.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}