{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Optimizing in a unbounded space\n",
    "\n",
    "Sometimes we have the inclusion bounds, like $[0, \\inf)$, or $[0, 1]$. An option is to lean on `optimize.minimize` to enforce these bounds, but alternatively we can express the problem a domain that has no constraint, and then project back to the original domain. We'll use an example to make this concrete.\n",
    "\n",
    "1. To make a parameter cover $[0, \\inf)$\n",
    "  1. Use `bounds = (0, None)`. \n",
    "  2. Use a parameter over $(-\\inf, \\inf)$, and put it into the exponential function, which has range $(0, \\inf)$. \n",
    "2. To make a parameter cover $[0, 1]$:\n",
    "  1. Use `bounds = (0, None)`. \n",
    "  2. Use a parameter over $(-\\inf, \\inf)$, and put it into the expit function, which has range $(0, 1)$. \n",
    "  \n",
    "$$ \\text{expit}(x) = \\frac{1}{1+\\exp(-x)} $$\n",
    "\n",
    "The B. cases typically have better convergence, and can take advantage of a larger family of optimization algorithms. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from autograd import numpy as np\n",
    "from autograd import elementwise_grad, value_and_grad, hessian\n",
    "from scipy.optimize import minimize\n",
    "\n",
    "df = pd.read_csv(\"../churn_data.csv\")\n",
    "T = df['T'].values\n",
    "E = df['E'].values\n",
    "\n",
    "breakpoints = np.array([28,  33,  58,  63,  88,  93,  117, 122, 148, 153])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.00103402 0.0332657  0.00104655 0.01516591 0.0005534  0.01317368\n",
      " 0.00099481 0.01216889 0.00125312 0.0068996 ]\n"
     ]
    }
   ],
   "source": [
    "# same model as last time, and note we have the `bounds` argument inactive.\n",
    "\n",
    "def cumulative_hazard(log_params, times):\n",
    "    # this is NumPy black magic to get piecewise hazards, let's chat after. \n",
    "    times = np.atleast_1d(times)\n",
    "    n = times.shape[0]\n",
    "    times = times.reshape((n, 1))\n",
    "    M = np.minimum(np.tile(breakpoints, (n, 1)), times)\n",
    "    M = np.hstack([M[:, tuple([0])], np.diff(M, axis=1)])\n",
    "    return np.dot(M, np.exp(log_params)) # diff here, use to be np.dot(M, params)\n",
    "\n",
    "hazard = elementwise_grad(cumulative_hazard, argnum=1)\n",
    "\n",
    "def survival_function(params, t):\n",
    "    return np.exp(-cumulative_hazard(params, t))\n",
    "\n",
    "def log_hazard(params, t):\n",
    "    return np.log(np.clip(hazard(params, t), 1e-25, np.inf))\n",
    "\n",
    "def log_likelihood(params, t, e):\n",
    "    return np.sum(e * log_hazard(params, t)) - np.sum(cumulative_hazard(params, t))\n",
    "\n",
    "def negative_log_likelihood(params, t, e):\n",
    "    return -log_likelihood(params, t, e)\n",
    "\n",
    "from autograd import value_and_grad\n",
    "\n",
    "results = minimize(\n",
    "        value_and_grad(negative_log_likelihood), \n",
    "        x0 = np.zeros(len(breakpoints)),\n",
    "        method=None, \n",
    "        args=(T, E),\n",
    "        jac=True,\n",
    "        bounds=None\n",
    ")\n",
    "\n",
    "log_estimates_ = results.x\n",
    "estimates_ = np.exp(log_estimates_)\n",
    "print(estimates_)\n",
    "# see previous Part 7 to confirm these are \"really close enough\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Congratulations! You've sucessfully implemented a pretty good approximations to the Python package [lifelines](https://lifelines.readthedocs.io/)! \n",
    "\n",
    "Let's move onto Part 9 (slides), which is off this track: interpreting output. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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
}