{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "O5RXyRXm4tIf"
   },
   "source": [
    "# A Cooling Coffee Cup with Polynomial Chaos Expansion\n",
    "\n",
    "In this tutorial we will perform a Polynomial Chaos Expansion for a\n",
    "model of a cooling coffee cup. The model uses Newton's law of cooling to\n",
    "evolve the temperature, $T$, over time ($t$) in an environment at\n",
    "$T_{env}$ :\n",
    "\n",
    "$$\\frac{dT(t)}{dt} = -\\kappa (T(t) -T_{env})$$\n",
    "\n",
    "The constant $\\kappa$ characterizes the rate at which the coffee cup\n",
    "transfers heat to the environment. In this example we will analyze this\n",
    "model using the polynomial chaos expansion (PCE) UQ algorithm. We will\n",
    "use a constant initial temperature $T_0 = 95 ^\\circ\\text{C}$, and vary\n",
    "$\\kappa$ and $T_{env}$ using a given probability distribution."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "ugDD5kSbAtw_"
   },
   "source": [
    "## Installing EasyVVUQ\n",
    "\n",
    "This example is based on https://github.com/vecma-project/VECMA-tutorials/blob/master/EasyVVUQ/Cooling_Cup/EasyVVUQ_Cooling_Cup.ipynb and was updated based on the latest EasyVVUQ interface available at https://github.com/UCL-CCS/EasyVVUQ/blob/dev/tutorials/easyvvuq_pce_tutorial.py"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "O5RXyRXm4tIf"
   },
   "source": [
    "## EasyVVUQ Script Overview\n",
    "\n",
    "We illustrate the intended workflow using the following basic example\n",
    "script, a python implementation of the cooling coffee cup model used in\n",
    "the textit{uncertainpy} documentation (code for which is in the\n",
    "tests/cooling/ subdirectory of the EasyVVUQ distribution directory). The\n",
    "code takes a small key/value pair input and outputs a comma separated\n",
    "value CSV) file.\n",
    "\n",
    "The input files for this tutorial are the *cooling\\_model* application\n",
    "(`cooling_model.py`), and an input template (`cooling.template`)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "execution": {
     "iopub.execute_input": "2025-07-21T13:16:37.840673Z",
     "iopub.status.busy": "2025-07-21T13:16:37.840566Z",
     "iopub.status.idle": "2025-07-21T13:21:12.056193Z",
     "shell.execute_reply": "2025-07-21T13:21:12.039824Z",
     "shell.execute_reply.started": "2025-07-21T13:16:37.840662Z"
    },
    "id": "Yd2XHwIE4tIg"
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import easyvvuq as uq\n",
    "import chaospy as cp\n",
    "import numpy as np\n",
    "import os\n",
    "from shutil import rmtree"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "uhwgG1Wc4tIj"
   },
   "source": [
    "## Create a new Campaign\n",
    "\n",
    "We start by creating an\n",
    "EasyVVUQ Campaign. The Campaign functions as as state machine for the VVUQ workflows. It uses a\n",
    "database (CampaignDB) to store information on both the target application\n",
    "and the VVUQ algorithms being employed. It also collects data from the simulations\n",
    "and can be used to store and resume your state."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2025-07-21T13:21:12.121017Z",
     "iopub.status.busy": "2025-07-21T13:21:12.120815Z",
     "iopub.status.idle": "2025-07-21T13:21:12.160008Z",
     "shell.execute_reply": "2025-07-21T13:21:12.159149Z",
     "shell.execute_reply.started": "2025-07-21T13:21:12.120996Z"
    }
   },
   "outputs": [],
   "source": [
    "work_dir = os.getcwd()\n",
    "campaign_work_dir = os.path.join(work_dir, \"easyvvuq_fd\")\n",
    "# clear the target campaign dir\n",
    "if os.path.exists(campaign_work_dir):\n",
    "    rmtree(campaign_work_dir)\n",
    "os.makedirs(campaign_work_dir)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "execution": {
     "iopub.execute_input": "2025-07-21T13:21:12.168965Z",
     "iopub.status.busy": "2025-07-21T13:21:12.168845Z",
     "iopub.status.idle": "2025-07-21T13:21:12.465169Z",
     "shell.execute_reply": "2025-07-21T13:21:12.461694Z",
     "shell.execute_reply.started": "2025-07-21T13:21:12.168956Z"
    },
    "id": "Cf5cFQXr4tIj"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "db_location = sqlite:////Volumes/UserData/dpc/GIT/EasyVVUQ/tutorials/correlated/easyvvuq_fd/campaign.db\n",
      "active_sampler_id = None\n",
      "campaign_name = coffee_pce\n",
      "campaign_dir = /Volumes/UserData/dpc/GIT/EasyVVUQ/tutorials/correlated/easyvvuq_fd/coffee_pcen17odu1s\n",
      "campaign_id = 1\n",
      "\n"
     ]
    }
   ],
   "source": [
    "db_location = \"sqlite:///\" + campaign_work_dir + \"/campaign.db\"\n",
    "my_campaign = uq.Campaign(name='coffee_pce', db_location=db_location, work_dir=campaign_work_dir)\n",
    "print(my_campaign)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "DgnFaMSP4tIn"
   },
   "source": [
    "## Parameter space definition\n",
    "\n",
    "The parameter space is defined using a dictionary. Each entry in the\n",
    "dictionary follows the format:\n",
    "\n",
    "`\"parameter_name\": {\"type\" : \"<value>\", \"min\": <value>, \"max\": <value>, \"default\": <value>}`\n",
    "\n",
    "With a defined type, minimum and maximum value and default. If the\n",
    "parameter is not selected to vary in the Sampler (see below) then the\n",
    "default value is used for every run. In this example, our full parameter\n",
    "space looks like the following: :"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "execution": {
     "iopub.execute_input": "2025-07-21T13:21:12.470896Z",
     "iopub.status.busy": "2025-07-21T13:21:12.470782Z",
     "iopub.status.idle": "2025-07-21T13:21:12.475200Z",
     "shell.execute_reply": "2025-07-21T13:21:12.474385Z",
     "shell.execute_reply.started": "2025-07-21T13:21:12.470887Z"
    },
    "id": "bHHJZWbf4tIo"
   },
   "outputs": [],
   "source": [
    "params = {\n",
    "    \"temp_init\": {\"type\": \"float\", \"min\": 0.0, \"max\": 100.0, \"default\": 95.0},\n",
    "    \"kappa\": {\"type\": \"float\", \"min\": 0.0, \"max\": 0.1, \"default\": 0.025},\n",
    "    \"t_env\": {\"type\": \"float\", \"min\": 0.0, \"max\": 40.0, \"default\": 15.0},\n",
    "    \"out_file\": {\"type\": \"string\", \"default\": \"output.csv\"}\n",
    "}"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "sW-haeK24tIq"
   },
   "source": [
    "## App Creation\n",
    "\n",
    "GenericEncoder is used for substituting values into application template input, which are then used to run the application.\n",
    "\n",
    "On the other hand, the application output file (in CSV format) is parsed and converted to the EasyVVUQ internal dictionary. The file is parsed in such a way that each column will appear as a vector QoI in the output dictionary.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "execution": {
     "iopub.execute_input": "2025-07-21T13:21:12.478701Z",
     "iopub.status.busy": "2025-07-21T13:21:12.478515Z",
     "iopub.status.idle": "2025-07-21T13:21:12.481838Z",
     "shell.execute_reply": "2025-07-21T13:21:12.481524Z",
     "shell.execute_reply.started": "2025-07-21T13:21:12.478691Z"
    },
    "id": "nHx45Wks4tIq"
   },
   "outputs": [],
   "source": [
    "from pathlib import Path\n",
    "\n",
    "encoder = uq.encoders.GenericEncoder(\n",
    "    template_fname=\"cooling.template\",\n",
    "    delimiter=\"$\",\n",
    "    target_filename=\"cooling_in.json\"\n",
    ")\n",
    "\n",
    "decoder = uq.decoders.SimpleCSV(\n",
    "    target_filename=\"output.csv\",\n",
    "    output_columns=[\"te\"]\n",
    ")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Define actions\n",
    "\n",
    "Actions are applied to each relevant run in the database."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2025-07-21T13:21:12.483494Z",
     "iopub.status.busy": "2025-07-21T13:21:12.483414Z",
     "iopub.status.idle": "2025-07-21T13:21:12.494323Z",
     "shell.execute_reply": "2025-07-21T13:21:12.493786Z",
     "shell.execute_reply.started": "2025-07-21T13:21:12.483486Z"
    }
   },
   "outputs": [],
   "source": [
    "execute = uq.actions.ExecuteLocal(\n",
    "    \"python3 {}/cooling_model.py cooling_in.json\".format(work_dir)\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2025-07-21T13:21:12.496237Z",
     "iopub.status.busy": "2025-07-21T13:21:12.496043Z",
     "iopub.status.idle": "2025-07-21T13:21:12.499241Z",
     "shell.execute_reply": "2025-07-21T13:21:12.498843Z",
     "shell.execute_reply.started": "2025-07-21T13:21:12.496228Z"
    }
   },
   "outputs": [],
   "source": [
    "actions = uq.actions.Actions(\n",
    "    uq.actions.CreateRunDirectory(root=campaign_work_dir, flatten=True),\n",
    "    uq.actions.Encode(encoder),\n",
    "    execute,\n",
    "    uq.actions.Decode(decoder)\n",
    ")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "-4HKkpZP4tIv"
   },
   "source": [
    "GenericEncoder performs simple text substitution into a supplied\n",
    "template, using a specified delimiter to identify where parameters\n",
    "should be placed. The template is shown below (\\$ is used as the\n",
    "delimiter). The template substitution approach is likely to suit most\n",
    "simple applications but in practice many large applications have more\n",
    "complex requirements, for example the multiple input files or the\n",
    "creation of a directory hierarchy. In such cases, users may write their\n",
    "own encoders by extending the BaseEncoder class.\n",
    "\n",
    "As can be inferred from its name SimpleCSV reads CVS files produced by the cooling model code. Again many applications output results in different formats, potentially requiring bespoke Decoders."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "FXcdBc0B4tIv"
   },
   "source": [
    "```\n",
    "{\n",
    "   \"T0\":\"$temp_init\",\n",
    "   \"kappa\":\"$kappa\",\n",
    "   \"t_env\":\"$t_env\",\n",
    "   \"out_file\":\"$out_file\"\n",
    "}\n",
    "```"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "uEig_muB4tIx"
   },
   "source": [
    "## Add the app to the campaign\n",
    "\n",
    "Having created an encoder, decoder and parameter space definition for our\n",
    "$cooling$ app, we can add it to our campaign. :"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "execution": {
     "iopub.execute_input": "2025-07-21T13:21:12.500734Z",
     "iopub.status.busy": "2025-07-21T13:21:12.500634Z",
     "iopub.status.idle": "2025-07-21T13:21:12.516356Z",
     "shell.execute_reply": "2025-07-21T13:21:12.516022Z",
     "shell.execute_reply.started": "2025-07-21T13:21:12.500724Z"
    },
    "id": "Ftd-tq6a4tIx"
   },
   "outputs": [],
   "source": [
    "my_campaign.add_app(\n",
    "    name=\"cooling\",\n",
    "    params=params,\n",
    "    actions=actions\n",
    ")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "EBUKUTyL4tIz"
   },
   "source": [
    "## The Sampler\n",
    "\n",
    "The user specified which parameters will vary and their corresponding\n",
    "distributions. In this case the `kappa` and `t_env` parameters are varied,\n",
    "both according to a specific distribution:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2025-07-21T13:21:12.518243Z",
     "iopub.status.busy": "2025-07-21T13:21:12.518071Z",
     "iopub.status.idle": "2025-07-21T13:21:12.536997Z",
     "shell.execute_reply": "2025-07-21T13:21:12.536682Z",
     "shell.execute_reply.started": "2025-07-21T13:21:12.518228Z"
    }
   },
   "outputs": [],
   "source": [
    "# Select distribution of the parameters defined next\n",
    "params_dist = 'normal_dep' # 'uniform' or 'normal' or 'normal_dep'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "execution": {
     "iopub.execute_input": "2025-07-21T13:21:12.550457Z",
     "iopub.status.busy": "2025-07-21T13:21:12.550239Z",
     "iopub.status.idle": "2025-07-21T13:21:12.580057Z",
     "shell.execute_reply": "2025-07-21T13:21:12.579721Z",
     "shell.execute_reply.started": "2025-07-21T13:21:12.550445Z"
    },
    "id": "skA-fE-c4tI0"
   },
   "outputs": [],
   "source": [
    "# [kappa, t_env]\n",
    "mu = [0.05, 20]\n",
    "stddev = [0.008, 1.5]\n",
    "\n",
    "\n",
    "    \n",
    "nvar = len(mu)\n",
    "cov_xy = 0.0096 #0.005 #[0, 0.0024, 0.0048, 0.0072, 0.0096, 0.012] #use the values in order to get corr = range(0,0.2,1)\n",
    "cov = np.array([[stddev[0]**2,cov_xy], [cov_xy,stddev[1]**2]])\n",
    "corr = np.array([ [cov[i][j]/(stddev[i]*stddev[j])  for j in range(nvar)]  for i in range(nvar)])\n",
    "\n",
    "joint = None\n",
    "if params_dist == 'uniform':\n",
    "    vary = {\n",
    "        \"kappa\": cp.Uniform(0.5*mu[0], 1.5*mu[0]),\n",
    "        \"t_env\": cp.Uniform(0.75*mu[1], 1.25*mu[1])\n",
    "    }\n",
    "elif params_dist == 'normal':\n",
    "    vary = {\n",
    "        \"kappa\": cp.Normal(mu[0], stddev[0]),\n",
    "        \"t_env\": cp.Normal(mu[1], stddev[1])\n",
    "    }\n",
    "    #cov = np.array([[stddev[0]**2,0], [0,stddev[1]**2]])\n",
    "    #joint = cp.MvNormal(mu, cov) # will use Rosenblatt\n",
    "elif params_dist == 'normal_dep':\n",
    "    vary = {\n",
    "        \"kappa\": cp.Normal(mu[0], stddev[0]),\n",
    "        \"t_env\": cp.Normal(mu[1], stddev[1])\n",
    "    }\n",
    "    joint = cp.MvNormal(mu, cov) # will use Rosenblatt\n",
    "    #joint = corr # will use Cholesky transform"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "execution": {
     "iopub.execute_input": "2025-07-21T13:21:12.581782Z",
     "iopub.status.busy": "2025-07-21T13:21:12.581694Z",
     "iopub.status.idle": "2025-07-21T13:21:12.607164Z",
     "shell.execute_reply": "2025-07-21T13:21:12.606867Z",
     "shell.execute_reply.started": "2025-07-21T13:21:12.581773Z"
    },
    "id": "1eDGvW3J4tI2"
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "2025-07-21 15:21:12,592:easyvvuq.sampling.fd:INFO:Performing relative analysis: True\n",
      "2025-07-21 15:21:12,593:easyvvuq.sampling.fd:INFO:Performing perturbation of the nodes, base value = mean, with delta = 0.03\n",
      "2025-07-21 15:21:12,595:easyvvuq.sampling.fd:INFO:Using user provided joint distribution with Rosenblatt transformation\n",
      "2025-07-21 15:21:12,595:easyvvuq.sampling.fd:INFO:Generated 5/5 samples for the FD scheme\n",
      "2025-07-21 15:21:12,596:easyvvuq.sampling.fd:INFO:Performing Rosenblatt transformation\n",
      "2025-07-21 15:21:12,598:easyvvuq.sampling.fd:DEBUG:The independent distribution consists of: J(Normal(mu=0.05, sigma=0.008), Normal(mu=20, sigma=1.5))\n",
      "2025-07-21 15:21:12,598:easyvvuq.sampling.fd:DEBUG:Using parameter permutation: ['kappa', 't_env']\n"
     ]
    }
   ],
   "source": [
    "perturbation = 0.03\n",
    "relative_analysis = True\n",
    "\n",
    "if params_dist == 'normal_dep':\n",
    "    my_sampler = uq.sampling.FDSampler(vary=vary, distribution=joint, perturbation=perturbation, relative_analysis=relative_analysis)\n",
    "else:\n",
    "    my_sampler = uq.sampling.FDSampler(vary=vary, perturbation=perturbation, relative_analysis=relative_analysis)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "snnsrp174tI4"
   },
   "source": [
    "Finally we set the campaign to use this sampler. :"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "execution": {
     "iopub.execute_input": "2025-07-21T13:21:12.609133Z",
     "iopub.status.busy": "2025-07-21T13:21:12.608971Z",
     "iopub.status.idle": "2025-07-21T13:21:12.627292Z",
     "shell.execute_reply": "2025-07-21T13:21:12.626872Z",
     "shell.execute_reply.started": "2025-07-21T13:21:12.609121Z"
    },
    "id": "E-Ruh53E4tI5"
   },
   "outputs": [],
   "source": [
    "my_campaign.set_sampler(my_sampler)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "IIDKFbhQ4tI7"
   },
   "source": [
    "Calling the campaign's draw\\_samples() method will cause the specified\n",
    "number of samples to be added as runs to the campaign database, awaiting\n",
    "encoding and execution. If no arguments are passed to draw\\_samples()\n",
    "then all samples will be drawn, unless the sampler is not finite. In\n",
    "this case PCESampler is finite (produces a finite number of samples) and\n",
    "we elect to draw all of them at once: :"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "execution": {
     "iopub.execute_input": "2025-07-21T13:21:12.629847Z",
     "iopub.status.busy": "2025-07-21T13:21:12.629547Z",
     "iopub.status.idle": "2025-07-21T13:21:12.651040Z",
     "shell.execute_reply": "2025-07-21T13:21:12.650615Z",
     "shell.execute_reply.started": "2025-07-21T13:21:12.629834Z"
    },
    "id": "RD5cupoa4tI7"
   },
   "outputs": [],
   "source": [
    "samples = my_campaign.draw_samples()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2025-07-21T13:21:12.653057Z",
     "iopub.status.busy": "2025-07-21T13:21:12.652873Z",
     "iopub.status.idle": "2025-07-21T13:21:12.655872Z",
     "shell.execute_reply": "2025-07-21T13:21:12.655598Z",
     "shell.execute_reply.started": "2025-07-21T13:21:12.653048Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "temperature,kappa\n",
      "20.000000,0.050000\n",
      "20.225000,0.051500\n",
      "19.775000,0.048500\n",
      "20.360000,0.050000\n",
      "19.640000,0.050000\n"
     ]
    }
   ],
   "source": [
    "kappas = [s['kappa'] for s in samples]\n",
    "temperatures = [s['t_env'] for s in samples]\n",
    "print(\"temperature,kappa\")\n",
    "for s in zip(temperatures, kappas):\n",
    "    print(\"%f,%f\" % (s[0], s[1]))"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "oUMfmupj4tI-",
    "tags": []
   },
   "source": [
    "## Execute Runs and Collation\n",
    "\n",
    "Create a directory hierarchy containing the encoded input files for every run that has not yet been\n",
    "completed.\n",
    "Calling `collate()` at any stage causes the campaign to\n",
    "aggregate decoded simulation output for all runs which have not yet been\n",
    "collated."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2025-07-21T13:21:12.657558Z",
     "iopub.status.busy": "2025-07-21T13:21:12.657446Z"
    }
   },
   "outputs": [],
   "source": [
    "my_campaign.execute().collate()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "nhtEXW0m4tJB"
   },
   "source": [
    "## Analysis\n",
    "\n",
    "This collated data is stored in the campaign database. An analysis\n",
    "element, here PCEAnalysis, can then be applied to the campaign's\n",
    "collation result. :"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Post-processing analysis\n",
    "output_columns = [\"te\"]\n",
    "\n",
    "my_analysis = uq.analysis.FDAnalysis(sampler=my_sampler, qoi_cols=output_columns)\n",
    "my_campaign.apply_analysis(my_analysis)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "rGgiN7cA4tJF"
   },
   "source": [
    "The output of the analysis can be next queried in order to get detailed statistics:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "-TvW8CIj4tJG"
   },
   "outputs": [],
   "source": [
    "# Get Descriptive Statistics\n",
    "results = my_campaign.get_last_analysis()\n",
    "mean = {}\n",
    "std = {}\n",
    "p10 = {}\n",
    "p90 = {}\n",
    "\n",
    "sobols_first = {}\n",
    "derivatives_first = {}\n",
    "sobols_first_conf = {}\n",
    "sobols_total = {}\n",
    "\n",
    "t = np.linspace(0, 200, 150)\n",
    "\n",
    "for output_column in output_columns:\n",
    "    mean[output_column] = results.describe(\"te\", \"mean\")\n",
    "    std[output_column] = results.describe(\"te\", \"std\")\n",
    "    p10[output_column] = results.describe(\"te\", \"10%\")\n",
    "    p90[output_column] = results.describe(\"te\", \"90%\")\n",
    "    sobols_first[output_column] = {}\n",
    "    derivatives_first[output_column] = {}\n",
    "    sobols_total[output_column] = {}    \n",
    "    for param in vary.keys():\n",
    "        sobols_first[output_column][param] = results._get_sobols_first(output_column, param)\n",
    "        #s1_kappa = results._get_sobols_first('te', 'kappa')\n",
    "        #s1_t_env = results._get_sobols_first('te', 't_env')\n",
    "        \n",
    "        #a = np.array(sobols_first[output_column][param])\n",
    "        #np.savetxt(f'000_{param}.csv', a, delimiter=',')\n",
    "        \n",
    "        sobols_total[output_column][param] = results._get_sobols_total(output_column, param)\n",
    "        #st_kappa = results._get_sobols_total('te', 'kappa')\n",
    "        #st_t_env = results._get_sobols_total('te', 't_env')\n",
    "        \n",
    "        derivatives_first[output_column][param] = results._get_derivatives_first(output_column, param)\n",
    "\n",
    "\n",
    "#print(f'Mean:\\n{mean}')\n",
    "#print(f'Std:\\n{std}')\n",
    "#print(f'p10:\\n{p10}')\n",
    "#print(f'p90:\\n{p90}')"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Data serialization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from serializer import serialize_yaml\n",
    "from serializer import deserialize_yaml\n",
    "from joblib import dump, load\n",
    "\n",
    "data = serialize_yaml(my_campaign=my_campaign,\n",
    "        output_name= \"sobols.yml\", output_columns=output_columns,\n",
    "        params_dist=params_dist, vary=vary,\n",
    "        mu=mu, stddev=stddev, cov=cov,\n",
    "        my_sampler=my_sampler, polynomial_order=0, regression=0,\n",
    "        time=t, mean=mean, variance=std, p10=p10, p90=p90,\n",
    "        derivatives_first=derivatives_first,\n",
    "        sobols_first=sobols_first, sobols_total=sobols_total)\n",
    "\n",
    "# Serialize the surrogate models into a file\n",
    "filename = campaign_work_dir + \"/surrogates.joblib\"\n",
    "dump(results.surrogate(), filename)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "XxFsGUPdQma0"
   },
   "source": [
    "## Visualization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import matplotlib\n",
    "from matplotlib import cm\n",
    "from matplotlib.ticker import LinearLocator\n",
    "import numpy as np"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Sensitivity indices"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "font = {'weight' : 'normal',\n",
    "        'size'   : 15}\n",
    "\n",
    "matplotlib.rc('font', **font)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_data(data, method, params, output_column):\n",
    "    fig, (ax1) = plt.subplots(nrows=1, ncols=1, figsize=(16,6))\n",
    "    \n",
    "    if output_column == \"te\":\n",
    "        output_label = \"Temperature $T(t)$ °C\"\n",
    "        time_label = ('Time t (min)')\n",
    "        param_labels = {\"kappa\": \"$\\kappa$\", \"t_env\":\"$T_{env}$\"}\n",
    "    else:\n",
    "        output_label = output_column\n",
    "        time_label = ('Time t (hr)')\n",
    "\n",
    "    # empty for e.g. SCsampler since this statistic is not supported\n",
    "    if np.array(data[output_column][\"model\"][\"p10\"]).size == 0:\n",
    "        data[output_column][\"model\"][\"p10\"] = np.zeros(len(data[output_column][\"model\"][\"time\"]))        \n",
    "    if np.array(data[output_column][\"model\"][\"p90\"]).size == 0:\n",
    "        data[output_column][\"model\"][\"p90\"] = np.zeros(len(data[output_column][\"model\"][\"time\"]))\n",
    "    \n",
    "    ax1.plot(data[output_column][\"model\"][\"time\"], data[output_column][\"model\"][\"mean\"], 'g-', alpha=0.75, label='Mean')\n",
    "    ax1.plot(data[output_column][\"model\"][\"time\"], data[output_column][\"model\"][\"p10\"], 'b-', alpha=0.25)\n",
    "    ax1.plot(data[output_column][\"model\"][\"time\"], data[output_column][\"model\"][\"p90\"], 'b-', alpha=0.25)\n",
    "    ax1.fill_between(\n",
    "      data[output_column][\"model\"][\"time\"],\n",
    "      data[output_column][\"model\"][\"p10\"],\n",
    "      data[output_column][\"model\"][\"p90\"],\n",
    "      alpha=0.25,\n",
    "      label='90% prediction interval')\n",
    "    ax1.set_xlabel(time_label)\n",
    "    ax1.set_ylabel(output_label, color='b')\n",
    "    ax1.tick_params('y', colors='b')\n",
    "    ax1.legend()\n",
    "\n",
    "    ax1t = ax1.twinx()\n",
    "    ax1t.plot(data[output_column][\"model\"][\"time\"], data[output_column][\"model\"][\"variance\"], 'r-', alpha=0.5, label='Variance')\n",
    "    ax1t.set_ylabel('Variance', color='r')\n",
    "    ax1t.tick_params('y', colors='r')\n",
    "    ax1t.legend(frameon=False)\n",
    "\n",
    "    ax1.grid()\n",
    "    ax1.set_title(f'Statistical Moments {method}')\n",
    "    \n",
    "    plt.subplots_adjust(wspace=0.35)\n",
    "    plt.savefig(f'SAmodel_{method}.pdf', bbox_inches='tight')\n",
    "    plt.show()\n",
    "    \n",
    "        \n",
    "    #################################################################\n",
    "    has_derivative_analysis = \"derivatives_first\" in data[output_column][params[0]]\n",
    "    nrows = 1\n",
    "    ncols = 2 if has_derivative_analysis else 1\n",
    "    figure, axes = plt.subplots(nrows=nrows, ncols=ncols, figsize=(15,5))\n",
    "\n",
    "    idx = 1\n",
    "    plt.subplot(nrows, ncols, idx)\n",
    "    \n",
    "    sum_s1 = np.array(data[output_column][params[0]][\"sobols_first\"])\n",
    "    sum_st = np.array(data[output_column][params[0]][\"sobols_total\"])\n",
    "    for param in params[1:]:\n",
    "        sum_s1 = sum_s1 + np.array(data[output_column][param][\"sobols_first\"])\n",
    "        sum_st = sum_st + np.array(data[output_column][param][\"sobols_total\"])\n",
    "    \n",
    "    # empty for e.g. SCsampler since this statistic is not supported\n",
    "    if sum_st.size == 0:\n",
    "        sum_st = np.zeros(len(data[output_column][\"model\"][\"time\"]))\n",
    "    \n",
    "    for param in params:\n",
    "        if output_column == \"te\":\n",
    "            param_label = param_labels[param]\n",
    "        else:\n",
    "            param_label = param\n",
    "        plt.plot(data[output_column][\"model\"][\"time\"], data[output_column][param][\"sobols_first\"], label=param_label)\n",
    "        \n",
    "    plt.plot(data[output_column][\"model\"][\"time\"], sum_s1, '--', label='Sum')\n",
    "    #plt.plot(data[output_column][\"model\"][\"time\"], sum_st, '--', label='Sobol Total Sum')\n",
    "\n",
    "    plt.xlabel(time_label)\n",
    "    plt.ylabel('First-order Sobol indices')\n",
    "    plt.title(f'Variance-based {method}')\n",
    "    plt.grid(which='both')\n",
    "    plt.ylim([-0.1, 1.2])\n",
    "    plt.legend()\n",
    "    \n",
    "    ####################################################################\n",
    "    if has_derivative_analysis:\n",
    "        idx = idx + 1\n",
    "        plt.subplot(nrows, ncols, idx)\n",
    "        for param in params:\n",
    "            if output_column == \"te\":\n",
    "                param_label = param_labels[param]\n",
    "            else:\n",
    "                param_label = param\n",
    "            plt.plot(data[output_column][\"model\"][\"time\"], data[output_column][param][\"derivatives_first\"], label=param_label)\n",
    "            #plt.semilogy(data[output_column][\"model\"][\"time\"], np.abs(data[output_column][param][\"derivatives_first\"]), label=param_label)\n",
    "\n",
    "        plt.legend()\n",
    "        plt.title(f\"Derivative-based {method}\")\n",
    "        plt.xlabel(time_label)\n",
    "        plt.ylabel(\"Sensitivity index\")\n",
    "        plt.grid(visible=True, which='both')    \n",
    "    ####################################################################\n",
    "\n",
    "    figure.tight_layout(pad=3.0)\n",
    "    plt.savefig(f'SAerror_SobolSum_{method}.pdf', bbox_inches='tight')\n",
    "    plt.show()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The FDSampler approximates the model derivatives using finite-differences, it is the only output this analysis can provide. The Statistical Moments and Variance-based indices are undefined."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "name = \"\"\n",
    "params = list(vary.keys())\n",
    "outputs = output_columns[0]\n",
    "plot_data(data, name, params , outputs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_derivatives(data, method, params, output_column):\n",
    "    \n",
    "    if output_column == \"te\":\n",
    "        output_label = \"Temperature $T(t)$ °C\"\n",
    "        time_label = ('Time t (min)')\n",
    "        param_labels = {\"kappa\": \"$\\kappa$\", \"t_env\":\"$T_{env}$\"}\n",
    "    else:\n",
    "        output_label = output_column\n",
    "        time_label = ('Time t (hr)')\n",
    "    \n",
    "\n",
    "    has_derivative_analysis = \"derivatives_first\" in data[output_column][params[0]]\n",
    "    nrows = 1\n",
    "    ncols = 2 if has_derivative_analysis else 1\n",
    "    figure, axes = plt.subplots(nrows=nrows, ncols=ncols, figsize=(15,5))\n",
    "\n",
    "    \n",
    "    if has_derivative_analysis:\n",
    "        idx = 1\n",
    "        plt.subplot(nrows, ncols, idx)\n",
    "    \n",
    "        for idd, param in enumerate(params):\n",
    "            if output_column == \"te\":\n",
    "                param_label = param_labels[param]\n",
    "            else:\n",
    "                param_label = param\n",
    "            plt.plot(data[output_column][\"model\"][\"time\"], data[output_column][param][\"derivatives_first\"], label=param_label)\n",
    "\n",
    "        plt.legend()\n",
    "        plt.title(f\"Derivative-based {method} (linear scale)\")\n",
    "        plt.xlabel(time_label)\n",
    "        plt.ylabel(\"Sensitivity index\")\n",
    "        plt.grid(visible=True, which='both') \n",
    "    \n",
    "    ####################################################################\n",
    "    if has_derivative_analysis:\n",
    "        idx = idx + 1\n",
    "        plt.subplot(nrows, ncols, idx)\n",
    "        for idd, param in enumerate(params):\n",
    "            if output_column == \"te\":\n",
    "                param_label = param_labels[param]\n",
    "            else:\n",
    "                param_label = param\n",
    "            plt.semilogy(data[output_column][\"model\"][\"time\"], np.abs(data[output_column][param][\"derivatives_first\"]), label=param_label)\n",
    "\n",
    "        plt.legend()\n",
    "        plt.title(f\"Derivative-based {method} (logarithmic scale)\")\n",
    "        plt.xlabel(time_label)\n",
    "        plt.ylabel(\"Sensitivity index\")\n",
    "        plt.grid(visible=True, which='both')    \n",
    "    ####################################################################\n",
    "\n",
    "    figure.tight_layout(pad=3.0)\n",
    "    plt.savefig(f'SAerror_derivativesLog_{method}.pdf', bbox_inches='tight')\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_derivatives(data, name, params , outputs)"
   ]
  }
 ],
 "metadata": {
  "colab": {
   "collapsed_sections": [],
   "name": "EasyVVUQ - Cooling Coffee Cup Example",
   "provenance": []
  },
  "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.12.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}