06-emulation_GP.ipynb

{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div style=\"float:right\">\n",
    "    <table>\n",
    "    <tr>\n",
    "        <td> <img src=\"figs/kcl_logo.png\" alt=\"Drawing\" style=\"width:54px;height:40px;\"/> \n",
    "        <td> <img src=\"figs/nceo_logo.png\" alt=\"Drawing\" style=\"width:200px;height:40px;\"/> \n",
    "        <td> <img src=\"figs/multiply_logo.png\" alt=\"Drawing\" style=\"width:40px;height:40px;\"/>\n",
    "    </tr>\n",
    "    </table>\n",
    "</div>\n",
    "\n",
    "# Emulating a typical radiative transfer model\n",
    "\n",
    "**Author:** Jose Gómez-Dans (NCEO & UCL)  `jose.gomez-dans@kcl.ac.uk`\n",
    "\n",
    "\n",
    "## The `gp_emulator` package\n",
    "\n",
    "This package was designed to emulate radiatie transfer models. The process entails the following steps:\n",
    "\n",
    "* Decide on what input parameters are required\n",
    "* Decide their ranges\n",
    "* Generate a training input parameter set\n",
    "* Run the model for each element of the input training set and store the outputs\n",
    "* Pass the input and output training pairs to the library, and let it fit the hyperparameters\n",
    "\n",
    "We can show how this works with an example of the PROSPECT+SAIL model.\n",
    "\n",
    "\n",
    "## Setting the input parameter ranges\n",
    "\n",
    "We can set the parameter names and their ranges simply by having lists with minimum and maximum values. This assumes a uniformly-distributed parameter distribution between those two boundaries, but other distributions are possible (we never had any reason to try them though!). We additionally set up the SRFs and other variables that need to be defined here… We train the model on 250 samples and test on (say) 100. 100 validation samples is probably too few, but for the sake of not waiting too much… ;-)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from functools import partial\n",
    "\n",
    "import numpy as np\n",
    "\n",
    "import gp_emulator\n",
    "\n",
    "import prosail\n",
    "\n",
    "# Spectral band definition. Just a top hat, with start and\n",
    "# end wavlengths as an example\n",
    "b_min = np.array( [ 620., 841, 459, 545, 1230, 1628, 2105] )\n",
    "b_max = np.array( [ 670., 876, 479, 565, 1250, 1652, 2155] )\n",
    "wv = np.arange ( 400, 2501 )\n",
    "passband = []\n",
    "\n",
    "# Number of training and validation samples\n",
    "n_train = 250\n",
    "n_validate = 100\n",
    "# Validation number is small, increase to a more realistic value\n",
    "# if you want\n",
    "\n",
    "# Define the parameter names and their ranges\n",
    "# Note that we are working here in transformed coordinates...\n",
    "\n",
    "# Define geometry. Each emulator is for one geometry\n",
    "sza = 30.\n",
    "vza = 0.\n",
    "raa = 0. # in degrees\n",
    "\n",
    "parameters = [ 'n', 'cab', 'car', 'cbrown', 'cw', 'cm', 'lai', 'ala', 'bsoil', 'psoil']\n",
    "min_vals = [ 0.8       ,  0.46301307,  0.95122942,  0.        ,  0.02829699,\n",
    "        0.03651617,  0.04978707,  0.44444444,  0.        ,  0.]\n",
    "max_vals = [ 2.5       ,  0.998002  ,  1.        ,  1.        ,  0.80654144,\n",
    "        0.84366482,  0.99501248,  0.55555556,  2.   , 1     ]\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We then require a function for calling the RT model. In the case of PROSAIL, we can do that easily from Python, in other models available in e.g. Fortran, you could have a function that calls the external model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def inverse_transform ( x ):\n",
    "    \"\"\"Inverse transform the PROSAIL parameters\"\"\"\n",
    "    x_out = x*1.\n",
    "    # Cab, posn 1\n",
    "    x_out[1] = -100.*np.log ( x[1] )\n",
    "    # Cab, posn 2\n",
    "    x_out[2] = -100.*np.log ( x[2] )\n",
    "    # Cw, posn 4\n",
    "    x_out[4] = (-1./50.)*np.log ( x[4] )\n",
    "    #Cm, posn 5\n",
    "    x_out[5] = (-1./100.)*np.log ( x[5] )\n",
    "    # LAI, posn 6\n",
    "    x_out[6] = -2.*np.log ( x[6] )\n",
    "    # ALA, posn 7\n",
    "    x_out[7] = 90.*x[7]\n",
    "    return x_out\n",
    "\n",
    "\n",
    "\n",
    "def rt_model ( x, passband=None, do_trans=True ):\n",
    "    \"\"\"A coupled land surface/atmospheric model, predicting refl from\n",
    "    land surface parameters. Thisfunction provides estimates of refl for\n",
    "    a particular illumination geometry.\n",
    "\n",
    "    The underlying land surface reflectance spectra is simulated using\n",
    "    PROSAIL. The input parameter ``x`` is a vector with the following components:\n",
    "\n",
    "        * ``n``\n",
    "        * ``cab``\n",
    "        * ``car``\n",
    "        * ``cbrown``\n",
    "        * ``cw``\n",
    "        * ``cm``\n",
    "        * ``lai``\n",
    "        * ``ala``\n",
    "        * ``bsoil``\n",
    "        * ``psoil``\n",
    "\n",
    "    \"\"\"\n",
    "    x, sza, vza, raa = x\n",
    "\n",
    "    # Invert parameter LAI\n",
    "    if do_trans:\n",
    "        x = inverse_transform ( x )\n",
    "    ################# surface refl with prosail #####################\n",
    "\n",
    "    surf_refl = prosail.run_prosail(x[0], x[1], x[2], x[3], x[4], x[5], x[6], x[7], 0.01, sza, vza, raa,\n",
    "                    rsoil=x[8], psoil=x[9])\n",
    "    if passband is None:\n",
    "        return surf_refl\n",
    "    else:\n",
    "        return surf_refl[passband].mean()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we loop over all the bands, and prepare the emulators. we do this by using the create_emulator_validation function, that does everything you’d want to do… We just stuff the emulator, training and validation sets in one list for convenience."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Doing band 1\n"
     ]
    },
    {
     "ename": "KeyboardInterrupt",
     "evalue": "",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mKeyboardInterrupt\u001b[0m                         Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-3-b8ba7696c5e3>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[1;32m     11\u001b[0m     x = gp_emulator.create_emulator_validation (simulator, parameters, min_vals, max_vals,\n\u001b[1;32m     12\u001b[0m                                 \u001b[0mn_train\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mn_validate\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdo_gradient\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mTrue\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 13\u001b[0;31m                                 n_tries=15, args=(30, 0, 0) )\n\u001b[0m\u001b[1;32m     14\u001b[0m     \u001b[0mretval\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mappend\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0mx\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/.local/lib/python3.6/site-packages/gp_emulator-1.6.10-py3.6.egg/gp_emulator/emulation_helpers.py\u001b[0m in \u001b[0;36mcreate_emulator_validation\u001b[0;34m(f_simulator, parameters, minvals, maxvals, n_train, n_validate, do_gradient, fix_params, thresh, n_tries, args, n_procs)\u001b[0m\n\u001b[1;32m    327\u001b[0m     \u001b[0;32mif\u001b[0m \u001b[0mtraining_set\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mndim\u001b[0m \u001b[0;34m==\u001b[0m \u001b[0;36m1\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    328\u001b[0m         \u001b[0mgp\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mGaussianProcess\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0msamples\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mtraining_set\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 329\u001b[0;31m         \u001b[0mgp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mlearn_hyperparameters\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mn_tries\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mn_tries\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    330\u001b[0m     \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    331\u001b[0m         gp = MultivariateEmulator(\n",
      "\u001b[0;32m~/.local/lib/python3.6/site-packages/gp_emulator-1.6.10-py3.6.egg/gp_emulator/GaussianProcess.py\u001b[0m in \u001b[0;36mlearn_hyperparameters\u001b[0;34m(self, n_tries, verbose, x0)\u001b[0m\n\u001b[1;32m    255\u001b[0m                 \u001b[0mfirst\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;32mFalse\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    256\u001b[0m                 \u001b[0mtheta\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mx0\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 257\u001b[0;31m             \u001b[0mT\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_learn\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mtheta\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mverbose\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    258\u001b[0m             \u001b[0mlog_like\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mappend\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mT\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    259\u001b[0m             \u001b[0mparams\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mappend\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mT\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/.local/lib/python3.6/site-packages/gp_emulator-1.6.10-py3.6.egg/gp_emulator/GaussianProcess.py\u001b[0m in \u001b[0;36m_learn\u001b[0;34m(self, theta0, verbose)\u001b[0m\n\u001b[1;32m    216\u001b[0m                 \u001b[0mfactr\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m0.1\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    217\u001b[0m                 \u001b[0mpgtol\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m1e-20\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 218\u001b[0;31m                 \u001b[0miprint\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0miprint\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    219\u001b[0m             )\n\u001b[1;32m    220\u001b[0m         \u001b[0;32mexcept\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mlinalg\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mLinAlgError\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/.local/lib/python3.6/site-packages/scipy/optimize/lbfgsb.py\u001b[0m in \u001b[0;36mfmin_l_bfgs_b\u001b[0;34m(func, x0, fprime, args, approx_grad, bounds, m, factr, pgtol, epsilon, iprint, maxfun, maxiter, disp, callback, maxls)\u001b[0m\n\u001b[1;32m    197\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    198\u001b[0m     res = _minimize_lbfgsb(fun, x0, args=args, jac=jac, bounds=bounds,\n\u001b[0;32m--> 199\u001b[0;31m                            **opts)\n\u001b[0m\u001b[1;32m    200\u001b[0m     d = {'grad': res['jac'],\n\u001b[1;32m    201\u001b[0m          \u001b[0;34m'task'\u001b[0m\u001b[0;34m:\u001b[0m \u001b[0mres\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'message'\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/.local/lib/python3.6/site-packages/scipy/optimize/lbfgsb.py\u001b[0m in \u001b[0;36m_minimize_lbfgsb\u001b[0;34m(fun, x0, args, jac, bounds, disp, maxcor, ftol, gtol, eps, maxfun, maxiter, iprint, callback, maxls, **unknown_options)\u001b[0m\n\u001b[1;32m    333\u001b[0m             \u001b[0;31m# until the completion of the current minimization iteration.\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    334\u001b[0m             \u001b[0;31m# Overwrite f and g:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 335\u001b[0;31m             \u001b[0mf\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mg\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mfunc_and_grad\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mx\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    336\u001b[0m         \u001b[0;32melif\u001b[0m \u001b[0mtask_str\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mstartswith\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34mb'NEW_X'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    337\u001b[0m             \u001b[0;31m# new iteration\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/.local/lib/python3.6/site-packages/scipy/optimize/lbfgsb.py\u001b[0m in \u001b[0;36mfunc_and_grad\u001b[0;34m(x)\u001b[0m\n\u001b[1;32m    283\u001b[0m     \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    284\u001b[0m         \u001b[0;32mdef\u001b[0m \u001b[0mfunc_and_grad\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mx\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 285\u001b[0;31m             \u001b[0mf\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mfun\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mx\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    286\u001b[0m             \u001b[0mg\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mjac\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mx\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    287\u001b[0m             \u001b[0;32mreturn\u001b[0m \u001b[0mf\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mg\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/.local/lib/python3.6/site-packages/scipy/optimize/optimize.py\u001b[0m in \u001b[0;36mfunction_wrapper\u001b[0;34m(*wrapper_args)\u001b[0m\n\u001b[1;32m    291\u001b[0m     \u001b[0;32mdef\u001b[0m \u001b[0mfunction_wrapper\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mwrapper_args\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    292\u001b[0m         \u001b[0mncalls\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m+=\u001b[0m \u001b[0;36m1\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 293\u001b[0;31m         \u001b[0;32mreturn\u001b[0m \u001b[0mfunction\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mwrapper_args\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0margs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    294\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    295\u001b[0m     \u001b[0;32mreturn\u001b[0m \u001b[0mncalls\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mfunction_wrapper\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/.local/lib/python3.6/site-packages/gp_emulator-1.6.10-py3.6.egg/gp_emulator/GaussianProcess.py\u001b[0m in \u001b[0;36mloglikelihood\u001b[0;34m(self, theta)\u001b[0m\n\u001b[1;32m    117\u001b[0m             \u001b[0mHyperparameters\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    118\u001b[0m         \"\"\"\n\u001b[0;32m--> 119\u001b[0;31m         \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_set_params\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mtheta\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    120\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    121\u001b[0m         loglikelihood = (\n",
      "\u001b[0;32m~/.local/lib/python3.6/site-packages/gp_emulator-1.6.10-py3.6.egg/gp_emulator/GaussianProcess.py\u001b[0m in \u001b[0;36m_set_params\u001b[0;34m(self, theta)\u001b[0m\n\u001b[1;32m    177\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    178\u001b[0m         \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mtheta\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mtheta\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 179\u001b[0;31m         \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_prepare_likelihood\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    180\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    181\u001b[0m     \u001b[0;32mdef\u001b[0m \u001b[0m_learn\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mtheta0\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mverbose\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/.local/lib/python3.6/site-packages/gp_emulator-1.6.10-py3.6.egg/gp_emulator/GaussianProcess.py\u001b[0m in \u001b[0;36m_prepare_likelihood\u001b[0;34m(self)\u001b[0m\n\u001b[1;32m    101\u001b[0m         \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mQ\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mZ\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0mexp_theta\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mD\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m*\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0meye\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mn\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    102\u001b[0m         \u001b[0mL\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mlinalg\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcholesky\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mQ\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 103\u001b[0;31m         \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0minvQ\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mlinalg\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0minv\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mL\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mT\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdot\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mlinalg\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0minv\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mL\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    104\u001b[0m         \u001b[0;31m# self.invQ = np.linalg.inv ( self.Q )\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    105\u001b[0m         \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0minvQt\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdot\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0minvQ\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mtargets\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m<__array_function__ internals>\u001b[0m in \u001b[0;36minv\u001b[0;34m(*args, **kwargs)\u001b[0m\n",
      "\u001b[0;32m~/anaconda/lib/python3.6/site-packages/numpy/linalg/linalg.py\u001b[0m in \u001b[0;36minv\u001b[0;34m(a)\u001b[0m\n\u001b[1;32m    545\u001b[0m     \u001b[0msignature\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m'D->D'\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0misComplexType\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mt\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32melse\u001b[0m \u001b[0;34m'd->d'\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    546\u001b[0m     \u001b[0mextobj\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mget_linalg_error_extobj\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0m_raise_linalgerror_singular\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 547\u001b[0;31m     \u001b[0mainv\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0m_umath_linalg\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0minv\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0ma\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0msignature\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0msignature\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mextobj\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mextobj\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    548\u001b[0m     \u001b[0;32mreturn\u001b[0m \u001b[0mwrap\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mainv\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mastype\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mresult_t\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcopy\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mFalse\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    549\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mKeyboardInterrupt\u001b[0m: "
     ]
    }
   ],
   "source": [
    "retval = []\n",
    "for iband,bmin in enumerate ( b_min ):\n",
    "    # Looping over the bands....\n",
    "    print(\"Doing band %d\" % (iband+1))\n",
    "    passband = np.nonzero( np.logical_and ( wv >= bmin, wv <= b_max[iband] ) )\n",
    "    # Define the SRF for the current band\n",
    "    # Define the simulator for convenience\n",
    "    simulator = partial (rt_model, passband=passband)\n",
    "    # Actually create the training and validation parameter sets, train the emulators\n",
    "    # and return all that\n",
    "    x = gp_emulator.create_emulator_validation (simulator, parameters, min_vals, max_vals,\n",
    "                                n_train, n_validate, do_gradient=True,\n",
    "                                n_tries=15, args=(30, 0, 0) )\n",
    "    retval.append (x)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a rel=\"license\" href=\"http://creativecommons.org/licenses/by-nc/4.0/\"><img alt=\"Creative Commons Licence\" style=\"border-width:0\" src=\"https://i.creativecommons.org/l/by-nc/4.0/88x31.png\" /></a><br />This work is licensed under a <a rel=\"license\" href=\"http://creativecommons.org/licenses/by-nc/4.0/\">Creative Commons Attribution-NonCommercial 4.0 International License</a>."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "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.10.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}