Clinic_Notes_and_Tutorial.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Uncomment and run these commands if you're on the CUAHSI JupyterHub to prepare the execution environment.\n",
    "# This may take 10 minutes\n",
    "# ! conda uninstall -y ncl esmf\n",
    "# ! pip install statsmodels -U\n",
    "# ! sed '0,/dakota_clinic/s/dakota_clinic/base/' environment_everything_but_dakota.yml > environment_cuahsijh.yml\n",
    "# ! conda env update --file environment_cuahsijh.yml\n",
    "# ! make install-hydroshare"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import time\n",
    "import glob\n",
    "import subprocess\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "from pandas.api.types import CategoricalDtype\n",
    "from plotnine import *\n",
    "from sklearn.linear_model import LinearRegression\n",
    "\n",
    "TOP_DIR = os.getcwd()\n",
    "os.environ[\"TOP_DIR\"] = TOP_DIR"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Goals of this Clinic\n",
    "\n",
    "The scope of both calibration and the mechanics of using [Dakota](https://dakota.sandia.gov) are large. This two-hour clinic was designed to deliver the following information. \n",
    "- Conceptual introduction to calibration and a discussion of some of the major options in calibration methods. \n",
    "- Introduction to the [Dakota](https://dakota.sandia.gov) package. \n",
    "- Background for a simple example using data from [Clow (2014)](https://www.earth-syst-sci-data.net/6/201/2014/essd-6-201-2014.pdf) and a simple model of 1D heat diffusion. \n",
    "- Experience with a simple example in Dakota --- including understanding how you must set up your model so Dakota can run it, and what files Dakota uses and makes. \n",
    "- Knowledge of what User Guides and resources are available to you for future work. \n",
    "\n",
    "# Step 1: Introduction to calibration\n",
    "\n",
    "## What is calibration?\n",
    "\n",
    "Calibration is the act of finding the best available input values to a model in order to minimize (or sometimes maximize) one or more outputs.\n",
    "\n",
    "Often the outputs are are a comparison between an observation and its **simulated equivalent**. \n",
    "\n",
    "Sometimes also called optimization, regression, or parameter estimation. \n",
    "\n",
    "If we could write an analytical expression for our model, we could approach calibration by finding finding places in the input parameter space where the derivative of the function was zero and the curvature is convex. \n",
    "\n",
    "Many models (in general and for Earth's surface) do not have an analytical expression that relates inputs and outputs. Thus, we rely on numerical methods for calibration. \n",
    "\n",
    "These calibration methods think of our models as a black box that transform our inputs into outputs. They evaluate one or more sets of input values and then use the output values to determine a new set of inputs. Sometimes they require the gradient (or Hessian) of the outputs with respect to the inputs. \n",
    "\n",
    "We call the quantity the calibration method is minimizing the **objective function**. If we want to minimize more than one output (say, the difference between observed and simulated hydrographs at two locations within a region), then we either need to indicate the correct way to combine (weight) these two components of our objective function, or we need to use a approach called **multi-objective optimization**. These are both interesting topics, but we won't discuss either today. \n",
    "\n",
    "Packages like [Dakota](https://dakota.sandia.gov) or methods within the [scipy.optimize](https://docs.scipy.org/doc/scipy-0.14.0/reference/optimize.html) package are set up to take a black box model, and iteratively modify the input parameters in order to determine the value that minimizes the objective function.  \n",
    "\n",
    "<img src=\"figures/black_box.png\" alt=\"blackbox\" width=\"400\"/>\n",
    "\n",
    "I will call each time a model is used to calculate output values for a given set of input values *function evaluation* and *model evaluation* interchangeably. \n",
    "\n",
    "Depending on the method used for calibration, the output of the calibration method might be a single set of parameter values at which the best model output occurred. Other methods may additionally provide an uncertainty estimate of our parameter values. For example, a gradient-based method (defined in a moment) typically provides parameter estimates based on linear assumptions. Bayesian calibration methods provide an empirical posterior parameter distribution. \n",
    "\n",
    "\n",
    "## Why does this even matter?\n",
    "\n",
    "If you have a model that runs in a millisecond, and only a few model input parameters. Then it probably doesn't. You can do a grid search (we will do this below) and brute force your way to a good model data fit with little computational expenditure. You will know you have a global optimum. The number of model evaluations you need to explore a parameter space expands geometrically.\n",
    "\n",
    "If you have a model that runs in a millisecond and many parameters, maybe it matters. \n",
    "\n",
    "But if your model takes more than a few seconds (or minutes, or hours) to run, and you have more parameters, then you probably need a more computationally frugal way to determine the best parameter set. \n",
    "\n",
    "\n",
    "## Comments on objective functions, methods,  and parameters. \n",
    "\n",
    "Calibration/optimization methods are not one size fits all. Some of the major choices and trade offs are listed below. \n",
    "\n",
    "- Gradient vs derivative-free optimization methods\n",
    "    * [Gradient-based methods](https://en.wikipedia.org/wiki/Gradient_descent) follow the gradient of the objective function. They often take steps (in parameter space) proportional to the magnitude of the gradient. If the objective function is smooth and no local minima exist, then they are a great, computationally frugal approach. If there are local minima, then their results may be dependent on the starting point. Gradient-based methods often converge quickly. If there is substantial computational cost in numerically calculating gradients and/or Hessians then these methods may not be appropriate. \n",
    "    * [Derivative-free methods](https://en.wikipedia.org/wiki/Derivative-free_optimization) are methods designed for situations when information about the derivative of the objetive function is not available or impractical. These methods often require more function evaluations than a gradient-based methods. However, they may be more likely to find values at or near a global minimum. \n",
    "- Complex model vs statistical surrogate\n",
    "    * Some methods toggle back and forth between making function evaluations of the \"real\" model (termed the **complex** model and constructing/using a statistical surrogate. A method we will use below (Efficient Global Optimization) does this using a Gaussian process surrogate. \n",
    "- Constrained vs. Unconstrained Parameters\n",
    "    * If parameter values can be constrained (e.g. by what is physically reasonable) this may be helpful. \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Step 2: Problem statement and dataset\n",
    "\n",
    "Today's clinic will use a model of 1D diffusion of heat in the Earth's crust and data from [Clow (2014)](https://www.earth-syst-sci-data.net/6/201/2014/essd-6-201-2014.pdf).\n",
    "\n",
    "In Earth's crust, the equation for temperature, $T$, as a function of depth, $z$, is given as\n",
    "\n",
    "$\\frac{\\partial T}{\\partial t} = \\kappa \\frac{\\partial T}{\\partial z}$\n",
    "\n",
    "where $\\kappa$ [m2/s] is the thermal diffusivity. \n",
    "\n",
    "Given a heat flux from below of $Q_m$ [J/(m2 s)] and a surface temperature of $T_0$ [K], the steady solution for $T(z)$ at steady state is given as\n",
    "\n",
    "$T(z) = \\frac{Q_m}{k} z + T_0$\n",
    "\n",
    "where $k$ [J / (m s K)] is the thermal conductivity. \n",
    "\n",
    "Here, the relationship between $k$ and $\\kappa$ depends on the heat capacity, $c$ [J/(kg K)], and the density, $\\rho$ [kg/m3].\n",
    "\n",
    "$ \\kappa = \\frac{k}{\\rho c}$\n",
    "\n",
    "Typical values:\n",
    "- $\\kappa$ = 3E-7 m2/s.\n",
    "- $Q_m$ = 0.05 W/m2\n",
    "- $k$ = 2 J / (m2 s)\n",
    "- $dT/dz$ = 0.03 K/m (30 degrees per kilometer)\n",
    "\n",
    "On the North Slope of Alaska, repeat observations of temperature in boreholes show two main features [(Clow, 2014)](https://www.earth-syst-sci-data.net/6/201/2014/essd-6-201-2014.pdf):\n",
    " \n",
    "  1. Steady state profiles at depth and warming at the surface. This makes sense given modern anthropogenic climate change. \n",
    "  2.  Cooling through time at depth. This makes sense given the thermal disturbance of drilling. \n",
    "\n",
    "This paper and the associated data files can be found in the folder `resources`. \n",
    "\n",
    "The next cell plots the temperatures as a function of time and depth for a well at Awuna (AWU). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "files = glob.glob(os.path.join(TOP_DIR, \"resources\", \"clow_2014\", \"G10015\", \"AWU*.txt\"))\n",
    "dfs = []\n",
    "for path in files:\n",
    "    date = os.path.split(path)[-1].split(\".\")[0].split(\"_\")[-1]\n",
    "    site = os.path.split(path)[-1].split(\".\")[0].split(\"_\")[0]\n",
    "    year = int(date[:2])\n",
    "    if year < 19:\n",
    "        year += 2000\n",
    "    else:\n",
    "        year += 1900\n",
    "    tdf = pd.read_csv(path, header=22, skip_blank_lines=False, delim_whitespace=True)\n",
    "    tdf[\"Date\"] = date\n",
    "    tdf[\"Year\"] = year\n",
    "    tdf[\"Site\"] = site\n",
    "    dfs.append(tdf)\n",
    "df = pd.concat(dfs).sort_values([\"Site\", \"Date\", \"Depth\"]).reset_index(drop=True)\n",
    "\n",
    "(\n",
    "    ggplot(df, aes(x=\"Temperature\", y=\"Depth\", color=\"factor(Year)\"))\n",
    "    + geom_path()\n",
    "    + facet_wrap(\"~Site\")\n",
    "    + scale_y_reverse()\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From this plot, we can see that the temperature disturbance (kink in the profile) reaches further and further down in Earth's surface with increasing year. \n",
    "\n",
    "We also see that (once the thermal disturbance of drilling has relaxed) the surface temperature is increasing through time. But how much did the surface temperature change? What pattern of surface temperature change through time is most consistent with these observations?\n",
    "\n",
    "We will use a model for heat diffusion address the following question:\n",
    "\n",
    "**What surface temperature change history is most consistent with the bore hole measurements at AWU?**\n",
    "\n",
    "To begin we will simulate surface temperature as a simple step change of size $\\Delta T$ at some number of years $Y$ before present. We will use only the 2012 bore hole. We will determine the initial temperature profile by making a linear fit to the bottom portion of the the profile. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_2012 = df[df[\"Year\"] == 2012]\n",
    "deep_portion = df_2012[\"Depth\"] > 100\n",
    "fit = LinearRegression().fit(\n",
    "    df_2012.Depth[deep_portion].values.reshape(-1, 1),\n",
    "    df_2012.Temperature[deep_portion].values.reshape(-1, 1),\n",
    ")\n",
    "\n",
    "df_2012 = df_2012.assign(\n",
    "    run_initial=fit.intercept_[0] + fit.coef_[0][0] * df_2012[\"Depth\"]\n",
    ")\n",
    "\n",
    "(\n",
    "    ggplot(\n",
    "        df_2012.melt(id_vars=[\"Depth\"], value_vars=[\"Temperature\", \"run_initial\"]),\n",
    "        aes(x=\"value\", y=\"Depth\", color=\"variable\"),\n",
    "    )\n",
    "    + geom_path()\n",
    "    + scale_y_reverse()\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Our parameters are:\n",
    "- `delta_T`, the temperature change, $\\Delta T$, between 0 and 20 degrees K\n",
    "- `duration_years`, the time before present, $Y$, when the temperature change occurred, between 15 and 120 years.\n",
    "\n",
    "Our objective function is the **root mean squared error** between observed and modeled temperature at the end of the model run. \n",
    "\n",
    "$ RMSE = \\left( \\frac{1}{N} \\sum^{N}_{i=1} \\left(o_i - s_i \\right)^2 \\right)^{0.5}$\n",
    "\n",
    "where there are $N$ observations ($o_i$) and simulated equivalents ($s_i$). \n",
    "\n",
    "*We are going to skip over the details of how the forward difference solution is implemented. We will use a dx of 10 m and a dt of 1 year. If you want to know more, look at the folder `heat` which contains the source code.** \n",
    "\n",
    "## Sneak Preview\n",
    "When we get to the end of the prepared demo, you will be set loose with a few options for exploration. I'll encourage you to either:\n",
    "\n",
    "  - Try out a different optimization method supported by Dakota.\n",
    "  - Change the objective function.\n",
    "  - Change the surface temperature model.\n",
    "  - Change something else about the model."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Step 3: Introduction to Dakota\n",
    "\n",
    "[Dakota](https://dakota.sandia.gov) is a software package for model analysis, calibration, and prediction under uncertainty developed by Sandia National Labs. Here are a couple of useful links:\n",
    "\n",
    "- [Online Documentation](https://dakota.sandia.gov/content/69-reference-manual)\n",
    "- [PDF manuals to download](https://dakota.sandia.gov/content/manuals)\n",
    "- [Downloads binaries or source code](https://dakota.sandia.gov/download.html)\n",
    "\n",
    "\n",
    "At most basic, Dakota provides A LOT of different methods for running your model using the same interface. As you will see below, we will only need to change a few lines in the Dakota input file to change optimization methods. \n",
    "\n",
    "In addition, Dakota has well thought out utilities for restarting partially completed computation (this gets you past issues with HPCC wall times) and writes output that documents exactly what/how/when it did what it did. \n",
    "\n",
    "It has a GUI if you are into that (I've never used it, but a recent update overhauled it and I've heard its great). \n",
    "\n",
    "\n",
    "## Why should I think about using Dakota?\n",
    "\n",
    "Dakota is a mature software package developed by a national lab. My experience is that it is well designed and well documented. It is extensive, which means that sometimes it may take a while to find what you are looking for in the documentation, but it is there and well written. \n",
    "\n",
    "Its been designed to deal with most any implementation of a black box model. \n",
    "\n",
    "Some reasons to not use Dakota:\n",
    "  - You really want to use a method it doesn't support\n",
    "  - You only need to do a few one-off things and would rather use matlab/python builtin-ish tools. \n",
    " "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## What is Dakota?\n",
    "[From Dakota's website](https://dakota.sandia.gov/content/about)\n",
    "### Summary\n",
    "The Dakota project delivers both state-of-the-art research and robust, usable software for optimization and UQ. Broadly, the Dakota software's advanced parametric analyses enable design exploration, model calibration, risk analysis, and quantification of margins and uncertainty with computational models. The Dakota toolkit provides a flexible, extensible interface between such simulation codes and its iterative systems analysis methods, which include:\n",
    "\n",
    "- optimization with gradient and nongradient-based methods;\n",
    "- uncertainty quantification with sampling, reliability, stochastic expansion, and epistemic methods;\n",
    "- parameter estimation using nonlinear least squares (deterministic) or Bayesian inference (stochastic); and\n",
    "- sensitivity/variance analysis with design of experiments and parameter study methods.\n",
    "\n",
    "These capabilities may be used on their own or as components within advanced strategies such as hybrid optimization, surrogate-based optimization, mixed integer nonlinear programming, or optimization under uncertainty.\n",
    "\n",
    "### Intended Use and Impact\n",
    "Computational methods developed in structural mechanics, heat transfer, fluid mechanics, shock physics, and many other fields of engineering can be an enormous aid to understanding the complex physical systems they simulate. Often, it is desired to use these simulations as virtual prototypes to obtain an optimized design for a particular system, or to develop confidence in performing predictions for systems that cannot be observed or tested directly. This effort seeks to enhance the utility of these computational simulations by enabling their use for a variety of iterative analyses, so that simulations may be used not just for single-point solutions, but also achieve broader impact in the areas of credible prediction and optimal design.\n",
    "\n",
    "This allows analysts to address the fundamental engineering questions of foremost importance to our programs, such as\n",
    "\n",
    "- \"What is the best design?\"\n",
    "- \"How safe is it?\"\n",
    "- \"How much confidence do I have in my answer?\"\n",
    "with respect to system performance objectives such as:\n",
    "\n",
    "- minimizing weight, cost, or defects;\n",
    "- limiting a critical temperature, stress, or vibration response;\n",
    "- maximizing performance, reliability, throughput, reconfigurability, agility, or design robustness\n",
    "\n",
    "A systematic approach for addressing these questions will lead to higher confidence in our computational simulations, in turn leading to better designs and improved system performance. Moreover, a reduced dependence on physical testing can shorten the design cycle and reduce development costs.\n",
    "\n",
    "### C++ Software Toolkit\n",
    "Written in C++, the Dakota toolkit provides a flexible, extensible interface between simulation codes and a variety of iterative systems analysis methods, including optimization, uncertainty quantification, deterministic/stochastic calibration, and parametric/sensitivity/variance analysis. These capabilities may be used on their own or as foundational components within advanced solution strategies that coordinate multiple computational models and iteration methods.  Initiated in 1994, it strives to span the research to production spectrum, providing both a mature tool for production use as well as a foundation for new algorithm research.\n",
    "\n",
    "Dakota is open source under GNU LGPL, with applications spanning defense programs for DOE and DOD, climate modeling, computational materials, nuclear power, renewable energy, and many others."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Step 4: The nitty gritty of how Dakota works\n",
    "\n",
    "Dakota is specified by an input file. We are now going to look at one, and talk about each part. \n",
    "\n",
    "A note on commands: Here when we use the `!` we are calling a system command (like you might do in a terminal) from a Jupyter notebook cell. Elsewhere in this notebook we will also use [``subprocess.call``](https://docs.python.org/2/library/subprocess.html#subprocess.call) to execute system commands from with a Python command. \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "! cat $TOP_DIR/analysis/dakota_01_grid.in"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are six blocks that may be used to specify a Dakota run. The order they appear does not matter. \n",
    "\n",
    "variables, interface, responses, model, method, and environment\n",
    "\n",
    "Parts of the input file\n",
    "- variables: input parameters\n",
    "- interface: how to call the model that transforms inputs into outputs\n",
    "- responses: output values\n",
    "- model: (optional in simple cases) a statement of which variables, interfaces, and responses blocks go together\n",
    "- method: which iterative method is used\n",
    "- environment: general settings like how Dakota creates output. \n",
    "\n",
    "**protip**: A single Dakota input file can specify running a very complex set of evaluations. You can chain together multiple methods (e.g. derivative free optimization to find the approximate location of the global minimum, followed by gradient-based methods). This is why multiple \"models\" can be defined for Dakota. \n",
    "\n",
    "We saw from the .in file that the *analysis_driver* is called to run the model. \n",
    "\n",
    "It has three parts:\n",
    " 1. Transform Dakotas params.in into something usable for our model.\n",
    " 2. Runs our thermal model.\n",
    " 3. Puts output into results.out \n",
    " \n",
    "Note that we specified the names `results.out` and `params.in`\n",
    " \n",
    "We also ask Dakota to keep track of all the files for us (useful for debugging, can be turned off if you don't need them or if you don't want to store them). This can be changed or removed. \n",
    "\n",
    "Lets look at it in a small amount of detail. Here I've turned line numbering on so that it is easier to discuss. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "! cat -n $TOP_DIR/analysis/template_dir/driver.py"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Also in `template_dir` is the input template. Its just a yaml file. \n",
    "\n",
    "Dakota doen't care what this is. From its perspective it is our job to ensure that when the `analysis_driver` is run, we transform input variables to output responses. \n",
    "\n",
    "**protip**: Dakota can run multiple analysis drivers per interface, it is also set up to run an `input_filter` and and and `output_filter` for each model evaluation. This can be a nice way to split up preprocessing, model evaluation, and post processing. \n",
    "\n",
    "**protip**: Dakota has lots of parallelization options. asynchronous evaluation_concurrency >1 is just the start. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "! cat $TOP_DIR/analysis/template_dir/input_template.yml"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "These curly braces are important. Dakota's dprepro tool (which we called from the driver script) replaces {deltaT} with the correct value for the input variable deltaT for that model iteration. \n",
    "\n",
    "Let's see what happens when we call Dakota from the command line. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "! dakota"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "How about the shell script that I've written to start Dakota. \n",
    "\n",
    "Note that this is where we specify how the log file, the restart file (.rst) and the .out file are named. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "! cat $TOP_DIR/analysis/start_01_grid.sh"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now lets call Dakota to run our project. This may take 5 minutes or more to execute as we will make Dakota execute all 100 necessary model evaluations. \n",
    "\n",
    "While this is happening, click\n",
    "\n",
    "File > Open \n",
    "\n",
    "and then explore the ``MULTIDIM_PARAM`` file structure within the ``analysis`` folder. \n",
    "\n",
    "This looks a little awkward because I've had to make it work for the following situations:\n",
    "1) Running on your computer\n",
    "2) Running on Binder\n",
    "3) Running on Travis for testing\n",
    "4) Running on Hydroshare\n",
    "\n",
    "Some of them don't wait until the command completes. This means that the remainder of the notebook will break. To get around this I wrote a silly function that checks how many files exist every 10 seconds. When that number doesn't change then its safe to move on. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "os.chdir(os.path.join(TOP_DIR, \"analysis\"))\n",
    "\n",
    "subprocess.call(\"./start_01_grid.sh\")\n",
    "\n",
    "def wait_until_done():\n",
    "    path = os.path.join(\"*\", \"*\", \"driver.py\")\n",
    "    nowfiles = len(glob.glob(path))\n",
    "    time.sleep(10)\n",
    "    while nowfiles != len(glob.glob(path)):\n",
    "        time.sleep(10)\n",
    "        nowfiles = len(glob.glob(path))\n",
    "\n",
    "wait_until_done()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next lets change the directory to the one we told Dakota to make (MULTIDIM_PARAM) and look at it. \n",
    "\n",
    "I've made all this work within the Jupyter Notebook, but if you want to see the file structure, you should be able to click File>Open and end up at the Jupyter filetree view. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "os.chdir(os.path.join(TOP_DIR, \"analysis\", \"MULTIDIM_PARAM\"))\n",
    "! ls"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Dakota has made a folder for each of the 100 model evaluations we have requested as part of our multi-dimensional parameter study. \n",
    "\n",
    "Lets look at what is in run.1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "os.chdir(os.path.join(TOP_DIR, \"analysis\", \"MULTIDIM_PARAM\", \"run.1\"))\n",
    "! ls"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It has the copied files, less the input_template.yml that we deleted. \n",
    "\n",
    "It also has two other files, the params.in and result.out files that Dakota uses to indicate the input variable values and to read the result values. Lets look at each as well as what inputs.yml now looks like."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "! cat $TOP_DIR/analysis/MULTIDIM_PARAM/run.1/params.in"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "! cat $TOP_DIR/analysis/MULTIDIM_PARAM/run.1/inputs.yml"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "! cat $TOP_DIR/analysis/MULTIDIM_PARAM/run.1/results.out"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, lets look at the Dakota output file. It starts by stating the time and Dakota version. \n",
    "\n",
    "It repeats the input file, then writes all of the output. You can make it more or less verbose depending on how the method blocks are specified. \n",
    "\n",
    "It ends with a summary of information relevant given the method and a summary of computational time used. \n",
    "\n",
    "**protip**: Dakota uses something called an Active Set Vector which helps it keep track of what combination of value, gradient, and Hessian it expects. \n",
    "\n",
    "If we scroll down to the end we can see that Dakota has summarized which function evaluation gave us the best output. This includes the parameter values and the residuals. It also gives us information about correlation matrices. \n",
    "\n",
    "The best value is something like a temperature change of 8.3 degrees that occurred 64 years ago. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "! cat -n $TOP_DIR/analysis/dakota_01_grid.out"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can plot the RMSE as a function of parameter space. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = pd.read_csv(os.path.join(TOP_DIR, \"analysis\", \"dakota_01_grid.dat\"), engine=\"python\", delim_whitespace=True)\n",
    "df = df.set_index([\"deltaT\", \"duration_years\"]).drop(columns=[\"interface\"])\n",
    "\n",
    "(\n",
    "    ggplot(df.reset_index(), aes(x=\"deltaT\", y=\"duration_years\", color=\"rmse\"))\n",
    "    + geom_point()\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Other methods\n",
    "Now, a grid search is really a lazy (or anti-lazy) approach to optimization. So next we will use Dakota to implement two additional methods. Next we will do a gradient based method (NL2SOL) and a derivative-free method that uses a surrogate (EGO). I'm not going to say anything more about how the methods work, but you can look at their references in the Dakota documentation. We can talk about them in discussion. \n",
    "\n",
    "First, look at what is different about their input files. This is summarized below, but its in white text, so you'll have to mouse over to reveal. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "! cat $TOP_DIR/analysis/dakota_02_nl2sol.in"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "! cat dakota_03_ego.in"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Some major differences:\n",
    "\n",
    "  - <span style=\"color:white\">They have different names for the tabular data file and the work directory</span> \n",
    "  - <span style=\"color:white\">The method block has a different keyword (e.g. efficient_global or nl2sol)</span> \n",
    "  - <span style=\"color:white\">NL2SOL has an initial point specified for input variables. This is where the gradient descent starts. </span> \n",
    "  - <span style=\"color:white\">Gradients are necessary for NL2SOL, so we have specified how they are calculated in the responses block. </span> \n",
    "  - <span style=\"color:white\">Seed is necessary for reproducible EGO (a default is used otherwise)</span> \n",
    "\n",
    "Lets run each of these new .in files. Again, this may take a few minutes. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "os.chdir(os.path.join(TOP_DIR, \"analysis\"))\n",
    "\n",
    "subprocess.call(\"./start_02_nl2sol.sh\")\n",
    "wait_until_done()\n",
    "subprocess.call(\"./start_03_ego.sh\")\n",
    "wait_until_done()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we will compare the results. First: how many function evaluations are there? And where did they occur?\n",
    "\n",
    "If the plotting commands look strange to you, its because I'm not using straight [matplotlib](http://matplotlib.org) but a package called [plotnine](https://plotnine.readthedocs.io/en/stable/index.html) which is the best implementation I'm aware of for R/Tidyverse's [ggplot2](http://ggplot2.tidyverse.org). If you know of a better one, let me know.  \n",
    "\n",
    "I find using [grammar of graphics](https://vita.had.co.nz/papers/layered-grammar.html) implementations to be the most efficient when plotting structured and/or self-described data (e.g., things stored in [pandas DataFrames](https://pandas.pydata.org))."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "files = glob.glob(os.path.join(TOP_DIR, \"analysis\", \"*.dat\"))\n",
    "dfs = []\n",
    "for file in files:\n",
    "    df = pd.read_csv(file, engine=\"python\", delim_whitespace=True)\n",
    "    df[\"method\"] = file.split(\".\")[0].split(\"_\")[-1]\n",
    "    dfs.append(df)\n",
    "\n",
    "df = pd.concat(dfs, ignore_index=True)\n",
    "method_cats = CategoricalDtype(categories=[\"grid\", \"nl2sol\", \"ego\"], ordered=True)\n",
    "df[\"method\"] = df[\"method\"].astype(method_cats)\n",
    "df = df.set_index([\"method\", \"deltaT\", \"duration_years\"]).drop(columns=[\"interface\"])\n",
    "\n",
    "(\n",
    "    ggplot(df.reset_index(), aes(x=\"deltaT\", y=\"duration_years\", color=\"%eval_id\"))\n",
    "    + geom_point()\n",
    "    + scale_color_cmap(name=\"jet\")\n",
    "    + facet_wrap(\"~method\")\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can see that the three methods make function calls in very different parts of parameter space. Now lets color the dots by the RMSE."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "(\n",
    "    ggplot(df.reset_index(), aes(x=\"deltaT\", y=\"duration_years\", color=\"rmse\"))\n",
    "    + geom_point()\n",
    "    + facet_wrap(\"~method\")\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "How are the results and number of evaluations influenced by method?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sum_df = df.drop(columns=[\"%eval_id\"]).groupby(\"method\").agg([np.count_nonzero, np.min])\n",
    "sum_df.columns = sum_df.columns.map(\"|\".join).str.strip(\"|\")\n",
    "\n",
    "(\n",
    "    ggplot(\n",
    "        sum_df.reset_index(),\n",
    "        (aes(x=\"rmse|count_nonzero\", y=\"rmse|amin\", color=\"method\")),\n",
    "    )\n",
    "    + geom_point()\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And what are the estimated best values for `deltaT` and `duration_years`?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "best_df = df[df.rmse.isin(sum_df[\"rmse|amin\"].values)].reset_index()\n",
    "print(best_df)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Discussion and Questions?\n",
    "\n",
    "\n",
    "# Exploration if time:\n",
    "* Explore other optimization methods. Start by going to the [Online Reference Manual](https://dakota.sandia.gov/content/69-reference-manual) and selecting Topics Area > Methods > Optimization and Calibration\n",
    "* Change or add a second component of the objective function.\n",
    "* Make the model (of surface temperature history) more complex.\n",
    "* Change something about the thermal model in the driver (or if you really want, in the `heat` module).\n"
   ]
  }
 ],
 "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.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}