10 SciML - 02 Optimization.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[![Binder](https://mybinder.org/badge_logo.svg)](https://notebooks.gesis.org/binder/v2/gh/jolin-io/fall-in-love-with-julia/main?filepath=10%20SciML%20-%2002%20Optimization.ipynb)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a href=\"https://www.jolin.io\" target=\"_blank\" rel=\"noreferrer noopener\">\n",
    "<img src=\"https://www.jolin.io/assets/Jolin/Jolin-Banner-Website-v1.1-darkmode.webp\">\n",
    "</a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Fall-in-love-with-Julia: Optimization in Julia 101\n",
    "\n",
    "an introduction session\n",
    "\n",
    "\n",
    "I am Stephan Sahm, and today we are going to learn about Optimization.jl\n",
    "\n",
    "1. Common interface to different optimizers\n",
    "2. Derivatives\n",
    "3. Constraints\n",
    "4. Neural networks - special support for minibatches\n",
    "5. Symbolic problem specification"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Optimization.jl: A Unified Optimization Package\n",
    "\n",
    "> Optimization.jl is a package with a scope that is beyond your normal global optimization package. Optimization.jl seeks to bring together all of the optimization packages it can find, local and global, into one unified Julia interface. This means, you learn one package and you learn them all! Optimization.jl adds a few high-level features, such as integrating with automatic differentiation, to make its usage fairly simple for most cases, while allowing all of the options in a single unified interface."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "using CommonSolve: solve  # generic interface\n",
    "import Optimization  # meta package\n",
    "import OptimizationOptimJL, OptimizationBBO, OptimizationMOI  # specific optimizers\n",
    "import ModelingToolkit  # symbolic support\n",
    "import Ipopt  # MOI (MathOptInterface ~ JuMP) optimizers\n",
    "import ForwardDiff, Zygote  # automatic derivatives\n",
    "\n",
    "import Plots  # plotting\n",
    "Plots.plotlyjs();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Common interface to different optimizers"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Rosenbrock function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rosenbrock(u, p) =  (p[1] - u[1])^2 + p[2] * (u[2] - u[1]^2)^2\n",
    "p  = [1.0, 100.0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "![wikipedia](https://upload.wikimedia.org/wikipedia/commons/thumb/6/68/Rosenbrock-contour.svg/450px-Rosenbrock-contour.svg.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The optimum is on (1,1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x, y = -2:0.01:2, 0:0.01:2\n",
    "Plots.surface(x, y, (x,y)->rosenbrock([x,y], p),\n",
    "    linealpha = 0.3, xlabel=\"x (u1)\", ylabel=\"y (u2)\",\n",
    "    zlabel=\"rosenbrock\", zscale=:log, c=:deep)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Specify optimization problem with `Optimization.jl`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "u0 = zeros(2)\n",
    "prob = Optimization.OptimizationProblem(rosenbrock, u0, p)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Solve problem by using specific solvers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "solve(prob, OptimizationOptimJL.ParticleSwarm())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "solve(prob, OptimizationOptimJL.SimulatedAnnealing())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "solve(prob, OptimizationOptimJL.NelderMead())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "overview optimizers, see also http://optimization.sciml.ai/stable\n",
    "<table><tbody><tr><th style=\"text-align: right\">Package</th><th style=\"text-align: center\">Local Gradient-Based</th><th style=\"text-align: center\">Local Hessian-Based</th><th style=\"text-align: center\">Local Derivative-Free</th><th style=\"text-align: center\">Local Constrained</th><th style=\"text-align: center\">Global Unconstrained</th><th style=\"text-align: center\">Global Constrained</th></tr><tr><td style=\"text-align: right\">BlackBoxOptim</td><td style=\"text-align: center\">❌</td><td style=\"text-align: center\">❌</td><td style=\"text-align: center\">❌</td><td style=\"text-align: center\">❌</td><td style=\"text-align: center\">✅</td><td style=\"text-align: center\">❌</td></tr><tr><td style=\"text-align: right\">CMAEvolutionaryStrategy</td><td style=\"text-align: center\">❌</td><td style=\"text-align: center\">❌</td><td style=\"text-align: center\">❌</td><td style=\"text-align: center\">❌</td><td style=\"text-align: center\">✅</td><td style=\"text-align: center\">❌</td></tr><tr><td style=\"text-align: right\">Evolutionary</td><td style=\"text-align: center\">❌</td><td style=\"text-align: center\">❌</td><td style=\"text-align: center\">❌</td><td style=\"text-align: center\">❌</td><td style=\"text-align: center\">✅</td><td style=\"text-align: center\">🟡</td></tr><tr><td style=\"text-align: right\">Flux</td><td style=\"text-align: center\">✅</td><td style=\"text-align: center\">❌</td><td style=\"text-align: center\">❌</td><td style=\"text-align: center\">❌</td><td style=\"text-align: center\">❌</td><td style=\"text-align: center\">❌</td></tr><tr><td style=\"text-align: right\">GCMAES</td><td style=\"text-align: center\">❌</td><td style=\"text-align: center\">❌</td><td style=\"text-align: center\">❌</td><td style=\"text-align: center\">❌</td><td style=\"text-align: center\">✅</td><td style=\"text-align: center\">❌</td></tr><tr><td style=\"text-align: right\">MathOptInterface</td><td style=\"text-align: center\">✅</td><td style=\"text-align: center\">✅</td><td style=\"text-align: center\">✅</td><td style=\"text-align: center\">✅</td><td style=\"text-align: center\">✅</td><td style=\"text-align: center\">🟡</td></tr><tr><td style=\"text-align: right\">MultistartOptimization</td><td style=\"text-align: center\">❌</td><td style=\"text-align: center\">❌</td><td style=\"text-align: center\">❌</td><td style=\"text-align: center\">❌</td><td style=\"text-align: center\">✅</td><td style=\"text-align: center\">❌</td></tr><tr><td style=\"text-align: right\">Metaheuristics</td><td style=\"text-align: center\">❌</td><td style=\"text-align: center\">❌</td><td style=\"text-align: center\">❌</td><td style=\"text-align: center\">❌</td><td style=\"text-align: center\">✅</td><td style=\"text-align: center\">🟡</td></tr><tr><td style=\"text-align: right\">NOMAD</td><td style=\"text-align: center\">❌</td><td style=\"text-align: center\">❌</td><td style=\"text-align: center\">❌</td><td style=\"text-align: center\">❌</td><td style=\"text-align: center\">✅</td><td style=\"text-align: center\">🟡</td></tr><tr><td style=\"text-align: right\">NLopt</td><td style=\"text-align: center\">✅</td><td style=\"text-align: center\">❌</td><td style=\"text-align: center\">✅</td><td style=\"text-align: center\">🟡</td><td style=\"text-align: center\">✅</td><td style=\"text-align: center\">🟡</td></tr><tr><td style=\"text-align: right\">Nonconvex</td><td style=\"text-align: center\">✅</td><td style=\"text-align: center\">✅</td><td style=\"text-align: center\">✅</td><td style=\"text-align: center\">🟡</td><td style=\"text-align: center\">✅</td><td style=\"text-align: center\">🟡</td></tr><tr><td style=\"text-align: right\">Optim</td><td style=\"text-align: center\">✅</td><td style=\"text-align: center\">✅</td><td style=\"text-align: center\">✅</td><td style=\"text-align: center\">✅</td><td style=\"text-align: center\">✅</td><td style=\"text-align: center\">✅</td></tr><tr><td style=\"text-align: right\">QuadDIRECT</td><td style=\"text-align: center\">❌</td><td style=\"text-align: center\">❌</td><td style=\"text-align: center\">❌</td><td style=\"text-align: center\">❌</td><td style=\"text-align: center\">✅</td><td style=\"text-align: center\">❌</td></tr></tbody></table>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sol = solve(prob, OptimizationOptimJL.NelderMead())\n",
    "sol"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fieldnames(typeof(sol))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sol.u"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Array(sol)  # very common"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sol.minimum"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sol.retcode"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sol.original"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "👉 check extra information on ParticleSwarm or SimulatedAnnealing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# your space"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Derivatives\n",
    "\n",
    "You can define derivatives manually, but most often you want to use automatic differentiation.\n",
    "\n",
    "There many kinds of auto differentiations but two are most important:\n",
    "- forward mode\n",
    "    - use for small number of learned parameters\n",
    "    - supports full julia\n",
    "- reverse mode\n",
    "    - use for large number of learned parameters\n",
    "    - supports large subset of julia"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# forward mode\n",
    "optf = Optimization.OptimizationFunction(rosenbrock, Optimization.AutoForwardDiff())\n",
    "prob = Optimization.OptimizationProblem(optf, u0, p)\n",
    "solve(prob, OptimizationOptimJL.BFGS())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# reverse mode\n",
    "optf = Optimization.OptimizationFunction(rosenbrock, Optimization.AutoZygote())\n",
    "prob = Optimization.OptimizationProblem(optf, u0, p)\n",
    "solve(prob, OptimizationOptimJL.BFGS())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "👉 check extra information `original` on those solutions\n",
    "\n",
    "👉 compare `minimum` to our derivative-free NelderMead solution above "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# your space"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. Constraints"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Lower and upper bounds:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Using gradient-free BlackBoxOptimization\n",
    "prob = Optimization.OptimizationProblem(rosenbrock, u0, p, lb = [-1.0,-1.0], ub = [1.0,1.0])\n",
    "solve(prob, OptimizationBBO.BBO_adaptive_de_rand_1_bin_radiuslimited())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# gradient method\n",
    "optf = Optimization.OptimizationFunction(rosenbrock, Optimization.AutoZygote())\n",
    "prob = Optimization.OptimizationProblem(optf, u0, p, lb = [-1.0,-1.0], ub = [1.0,1.0])\n",
    "solve(prob, OptimizationOptimJL.BFGS())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Extra constraint functions with bounds,\n",
    "\n",
    "for example:\n",
    "$$\n",
    "-\\infty <= u_1^2+u_2^2 <= 0.8 \\\\\n",
    "-1.0 <= u_1∗u_2 <= 2.0\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# constraints = (sum of squares, product)\n",
    "function cons(res, x, p)\n",
    "    res .= [x[1]^2+x[2]^2, x[1]*x[2]]\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "optprob = Optimization.OptimizationFunction(rosenbrock, Optimization.AutoForwardDiff(), cons = cons)\n",
    "prob = Optimization.OptimizationProblem(optprob, u0, p, lcons = [-Inf, -1.0], ucons = [0.8, 2.0])\n",
    "sol = solve(prob, OptimizationOptimJL.IPNewton())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sol.original"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To manually inspect our constraints we need to create a little array helper."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "res = zeros(2)\n",
    "cons(res, sol.u, p)\n",
    "res"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can use equality constraints, by having equal lowerbound and upperbounds for your constraints.\n",
    "\n",
    "$$\n",
    "u_1^2​+u_2^2​ = 1.0\\\\\n",
    "u_1​∗u_2​ = 0.5​\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# using symbolic derivatives via ModelingToolkit\n",
    "optprob = Optimization.OptimizationFunction(rosenbrock, Optimization.AutoModelingToolkit(), cons = cons)\n",
    "prob = Optimization.OptimizationProblem(optprob, u0, p, lcons = [1.0, 0.5], ucons = [1.0, 0.5])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# using OptimizationMOI under the hood\n",
    "# here we can directly use the original Optimizer\n",
    "sol = solve(prob, Ipopt.Optimizer())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "res = zeros(2)\n",
    "cons(res, sol.u, p)\n",
    "println(res)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "👉 try out an arbitrary other constraint and see whether you can solve it"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# your space "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4. Neural networks - special support for minibatches\n",
    "\n",
    "⚠️ the mybinder instance with 2GB memory is to small to run everything in one session. It is recommended to restart the kernel and continue here directly."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import Flux, OptimizationOptimisers\n",
    "import IterTools, MLUtils, NNlib\n",
    "import Plots, Optimization\n",
    "using CommonSolve: solve"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "👉 Generate 128 datapoints from the polynomial $y = x² - 2x$ and add some noise.\n",
    "\n",
    "You need `randn`, and you might use `range` for x.\n",
    "\n",
    "Plot it using `Plots.plot(x, y)`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# your space\n",
    "# x = ...\n",
    "# y = ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# making sure x and y are matrix of size (1,128)\n",
    "@assert length(x) == length(y) == 128\n",
    "x = size(x) == (128,) ? collect(x') : x\n",
    "y = size(y) == (128,) ? collect(y') : y"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nn_flux = Flux.Chain(\n",
    "    Flux.Dense(1, 16, NNlib.relu),\n",
    "    Flux.Dense(16, 1),\n",
    ")\n",
    "# we need to make the parameters explicit which we want to optimize\n",
    "# in Flux we can do this via destructure\n",
    "parameters_initial, reconstruct_nn_flux = Flux.destructure(nn_flux)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# this is how to get a prediction\n",
    "y_pred = reconstruct_nn_flux(parameters_initial)(x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "👉 plot both the true solution and our prediction\n",
    "\n",
    "(you may need `Plots.plot!` and `transpose`)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# your space"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# a loss function is a simple function\n",
    "function loss_flux(parameters, x, y)\n",
    "    y_pred = reconstruct_nn_flux(parameters)(x)\n",
    "    sum(abs2, y .- y_pred)\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "👉 calculate our initial loss"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# your space"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Minibatches 🙂"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "k = 10\n",
    "minibatches = MLUtils.DataLoader((x, y), batchsize = k)\n",
    "fieldnames(typeof(minibatches))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x1_batch, y1_batch = first(minibatches)\n",
    "Plots.plot(x1_batch', y1_batch')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "losses = Float64[]\n",
    "function callback(p, l)  # further outputs of the loss function are given as further input arguments to this callback function\n",
    "    push!(losses, l)\n",
    "    if length(losses) % 50 == 0\n",
    "        Plots.plot(losses, show = :inline, yscale = :log10,\n",
    "            label = \"loss\", xlabel = \"#epochs\", ylabel=\"loss (log10 scale)\")\n",
    "    end\n",
    "    return false  # return bool `halt`\n",
    "end\n",
    "\n",
    "_optfun = (θ, _, x_batch, y_batch) -> loss_flux(θ, x_batch, y_batch)\n",
    "optfun = Optimization.OptimizationFunction(_optfun, Optimization.AutoZygote())\n",
    "optprob = Optimization.OptimizationProblem(optfun, parameters_initial)\n",
    "\n",
    "numEpochs = 500\n",
    "sol = solve(\n",
    "    optprob,\n",
    "    OptimizationOptimisers.ADAM(0.01),\n",
    "    IterTools.ncycle(minibatches, numEpochs),\n",
    "    callback = callback,\n",
    ")\n",
    "parameters_learned = sol.minimizer"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "👉 calculate the final loss"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# your space"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "👉 plot our prediction vs the true solution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# more space"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "👉 plot our prediction vs the true solution over a much larger range (remember that x and y need to be row vectors)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# your space"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "⚠️ IMPORTANT: Note that this minibatch data support is not universal, but just supported by a few Optimizer backends. OptimizationOptimisers and OptimizationOptimJL for instance do support this extra data argument. However others like OptimizationMOI does not support it.\n",
    "\n",
    "This is not only a missing feature, but might really be ill-defined, depending on your solver and your problem. Mini-batches introduce a kind of randomness. Some problems (like differential equations) and some solvers loose their mathematical guarantees as soon as you add such stochasticity.\n",
    "\n",
    "⚠️ So be always cautious whether minibatches are really appropriate for your optimization problem!\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 5. Symbolic problem specification"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "using ModelingToolkit: @variables, @parameters, @named\n",
    "using CommonSolve: solve\n",
    "import ModelingToolkit, Optimization, OptimizationOptimJL"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "@variables x y\n",
    "@parameters a b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "loss = (a - x)^2 + b * (y - x^2)^2\n",
    "@named sys = ModelingToolkit.OptimizationSystem(loss, [x, y], [a, b])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Specifying parameters and initial state is a bit more complex: We need to map the symbols to values "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "u0 = [\n",
    "    x => 1.0\n",
    "    y => 2.0\n",
    "]\n",
    "p = [\n",
    "    a => 6.0\n",
    "    b => 7.0\n",
    "];"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we get numerous benefits: Symbolic auto differentiation, auto-parallelism, sparsification, & many more. You can even hierarchically nest systems to have it generate huge optimization problems."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "prob = Optimization.OptimizationProblem(sys, u0, p, grad=true, hess=true)\n",
    "solve(prob, OptimizationOptimJL.Newton())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For more details about symbolic problem descriptions check out the ModelingToolkit.jl [OptimizationSystem documentation](https://mtk.sciml.ai/dev/systems/OptimizationSystem/).\n",
    "\n",
    "Further details about Optimization.jl can be found at its [official documentation](http://optimization.sciml.ai/stable/), which was also the main source for the juypter-notebook at hand."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Thank you for your participation\n",
    "\n",
    "for questions or suggestions please contact me at stephan.sahm@jolin.io\n",
    "\n",
    "\n",
    "#### Sponsored by [Jolin.io](https://www.jolin.io)\n",
    "\n",
    "Jolin.io is an IT-consultancy focussing on Julia\n",
    "\n",
    "We are there to help you, if you want to\n",
    "- try out Julia at your company, or\n",
    "- transition Matlab, Fortran, R, Python, etc. to Julia\n",
    "- or speed up your existing Julia code\n",
    "\n",
    "<a href=\"https://www.jolin.io\" target=\"_blank\" rel=\"noreferrer noopener\">\n",
    "<img src=\"https://www.jolin.io/assets/Jolin/Jolin-Banner-Website-v1.1-darkmode.webp\">\n",
    "</a>"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 1.8.1",
   "language": "julia",
   "name": "julia-1.8"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.8.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}