{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "title: Nonlinear Modelling\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Originally Contributed by**: Arpit Bhatia"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This tutorial provides a breif introduction to nonlinear modelling in JuMP.\n",
    "For more details and specifics, visit the [JuMP docs](http://www.juliaopt.org/JuMP.jl/stable/nlp/)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Nonlinear Programs\n",
    "While we have already seen examples of linear, quadratic and conic programs,\n",
    "JuMP also supports other general smooth nonlinear (convex and nonconvex) optimization problems."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A JuMP model object can contain a mix of linear, quadratic, and nonlinear contraints or objective functions.\n",
    "Thus, a model object for a nonlinear program is constructed in the same way as before."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "using JuMP, Ipopt\n",
    "model = Model(optimizer_with_attributes(Ipopt.Optimizer, \"print_level\" => 0));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Variables\n",
    "Variables are modelled using the `@variable` macro as usual and\n",
    "a starting point may be provided by using the `start` keyword argument"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "@variable(model, x, start = 4)\n",
    "@variable(model, y, start = -9.66);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Parameters\n",
    "Only in the case of nonlinear models, JuMP offers a syntax for \"parameter\" objects\n",
    "which can refer to a numerical value."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "@NLparameter(model, p == 0.003); # Providing a starting value is necessary for parameters\n",
    "@NLparameter(model, l[i = 1:10] == 4 - i); # A collection of parameters"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `value` and `set_value` functions are used to query and update the value of a parameter respectively."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3.0"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "value(l[1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-4.0"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "set_value(l[1], -4)\n",
    "value(l[1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Parameters are useful since it's faster to modify a model in-place by changing the value of the parameter\n",
    "compared to creating an entirely new model object."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Expressions\n",
    "JuMP also supports the creation of arithmetic expressions which can then be inserted into\n",
    "constraints, the objective and other expressions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "@NLexpression(model, expr_1, sin(x))\n",
    "@NLexpression(model, expr_2, asin(expr_1)); # Inserting one expression into another"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are some [syntax rules](https://pkg.julialang.org/docs/JuMP/DmXqY/0.19.2/nlp/#Syntax-notes-1)\n",
    "which must be followed while writing a nonlinear expression."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that JuMP also supports linear and quadratic expression.\n",
    "You can find out more about this functionality in the [docs](https://pkg.julialang.org/docs/JuMP/DmXqY/0.19.2/expressions/)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Nonlinear Objectives and Constraints\n",
    "Nonlinear objectives and constraints are specified by using the `@NLobjective` and `@NLconstraint` macros."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "@NLconstraint(model, exp(x) + y^4 <= 0)\n",
    "@NLobjective(model, Min, tan(x) + log(y))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### User-defined Functions\n",
    "In addition to supporting a library of built-in functions,\n",
    "JuMP supports the creation of user-defined nonlinear functions to use within nonlinear expressions.\n",
    "The `register` function is used to enable this functionality."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "my_function(a,b) = (a * b)^-6 + (b / a)^3\n",
    "register(model, :my_function, 2, my_function, autodiff = true)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The arguements to the function are:\n",
    "- model for which the function is being registered\n",
    "- Julia symbol object corresponding to the name of the function\n",
    "- Number of arguments the function takes\n",
    "- name of the Julia method\n",
    "- instruction for JuMP to compute exact gradients automatically"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## MLE using JuMP"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since we already have a bit of JuMP experience at this point,\n",
    "let's try a modelling example and apply what we have learnt.\n",
    "In this example, we compute the maximum likelihood estimate (MLE) of\n",
    "the parameters of a Gaussian distribution i.e. the sample mean and variance."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If $X_{1}, \\ldots, X_{n}$ are an id sample from a population with pdf or pmf\n",
    "$f\\left(x | \\theta_{1}, \\ldots, \\theta_{k}\\right),$ the likelihood function is defined by"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "L(\\theta | \\mathbf{x})=L\\left(\\theta_{1}, \\ldots, \\theta_{k} | x_{1}, \\ldots, x_{n}\\right)=\\prod_{i=1}^{n} f\\left(x_{i} | \\theta_{1}, \\ldots, \\theta_{k}\\right)\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For each sample point $\\mathbf{x}$, let $\\hat{\\theta}(\\mathbf{x})$ be a parameter value\n",
    "at which $L(\\theta | \\mathbf{x})$ attains its maximum as a function of $\\theta,$ with $\\mathbf{x}$ held fixed.\n",
    "A maximum likelihood estimator (MLE) of the parameter $\\theta$ based on a sample $\\mathbf{X}$ is\n",
    "$\\hat{\\theta}(\\mathbf{X})$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Gaussian likelihood is -"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "L(\\theta | \\mathbf{x})=\\prod_{i=1}^{n} \\frac{1}{(2 \\pi)^{1 / 2}} e^{-(1 / 2)\\left(x_{i}-\\theta\\right)^{2}}=\\frac{1}{(2 \\pi)^{n / 2}} e^{(-1 / 2) \\Sigma_{i=1}^{n}\\left(x_{i}-\\theta\\right)^{2}}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In most cases, the natural logarithm of\n",
    "$L(\\theta | \\mathbf{x}), \\log L(\\theta | \\mathbf{x})$ (known as the log likelihood),\n",
    "is used rather than $L(\\theta | \\mathbf{x})$ directly.\n",
    "The reason is that the log likelihood is easier to differentiate.\n",
    "This substituion is possible because the log function is strictly increasing on $(0, \\infty)$,\n",
    "which implies that the extrema of $L(\\theta | \\mathbf{x})$ and $\\log L(\\theta | \\mathbf{x})$ coincide."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "┌ Info: Starting guess, mean: -1.0091959881909511, std: 1.1820542647634424\n",
      "└ @ Main In[9]:12\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "******************************************************************************\n",
      "This program contains Ipopt, a library for large-scale nonlinear optimization.\n",
      " Ipopt is released as open source code under the Eclipse Public License (EPL).\n",
      "         For more information visit http://projects.coin-or.org/Ipopt\n",
      "******************************************************************************\n",
      "\n",
      "μ = -0.02294386457202795\n",
      "mean(data) = -0.022943864572027916\n",
      "σ^2 = 1.0096978289614724\n",
      "var(data) = 1.0107085374810936\n",
      "MLE value: 0.0\n"
     ]
    }
   ],
   "source": [
    "using Random, Statistics\n",
    "\n",
    "Random.seed!(1234)\n",
    "\n",
    "n = 1_000\n",
    "data = randn(n)\n",
    "\n",
    "mle = Model(optimizer_with_attributes(Ipopt.Optimizer, \"print_level\" => 0))\n",
    "@NLparameter(mle, problem_data[i = 1:n] == data[i])\n",
    "μ0 = randn()\n",
    "σ0 = rand() + 1\n",
    "@info \"Starting guess, mean: $μ0, std: $σ0\"\n",
    "@variable(mle, μ, start = μ0)\n",
    "@variable(mle, σ >= 0.0, start = σ0)\n",
    "@NLexpression(mle, loglikelihood,\n",
    "    -(n / 2) * (log(2π) + 2 * log(σ)) - inv(2 * σ^2) * sum((xi - μ)^2 for xi in problem_data)\n",
    ")\n",
    "\n",
    "@NLobjective(mle, Max, loglikelihood)\n",
    "\n",
    "optimize!(mle)\n",
    "\n",
    "println(\"μ = \", value(μ))\n",
    "println(\"mean(data) = \", mean(data))\n",
    "println(\"σ^2 = \", value(σ)^2)\n",
    "println(\"var(data) = \", var(data))\n",
    "println(\"MLE value: \", exp(objective_value(mle)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "μ = -0.02294386457202795\n",
      "mean(data) = -0.010731153129307803\n",
      "σ^2 = 1.0096978289614724\n",
      "var(data) = 1.0173868093390144\n",
      "MLE objective: -1423.76408661786\n"
     ]
    }
   ],
   "source": [
    "# Changing the data\n",
    "\n",
    "data = randn(n)\n",
    "optimize!(mle)\n",
    "\n",
    "println(\"μ = \", value(μ))\n",
    "println(\"mean(data) = \", mean(data))\n",
    "println(\"σ^2 = \", value(σ)^2)\n",
    "println(\"var(data) = \", var(data))\n",
    "println(\"MLE objective: \", objective_value(mle))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Writing Convex Models"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Nonlinear solvers like Ipopt are usually local solvers.\n",
    "For convex problems, the local optima is also the global optima,\n",
    "and thus Ipopt is able to provide us with the correct solution.\n",
    "However, in case a problem is not written in the convex form, then we may be unable to solve it."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A tool that helps us in dealing with this issue is [Disciplined Convex Programming](https://dcp.stanford.edu) (DCP).\n",
    "DCP is a system for constructing mathematical expressions with known curvature from a given library of base functions.\n",
    "Specifically, it helps us to construct convex optimization models when possible,\n",
    "i.e. minimize convex function or maximize concave function and use constraints that are convex $f$ <= concave $g$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 1.5.1",
   "language": "julia",
   "name": "julia-1.5"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.5.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}