Q1.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "bizarre-water",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\u001b[32m\u001b[1m  Activating\u001b[22m\u001b[39m environment at `~/devel/HW0_S23/Project.toml`\n"
     ]
    }
   ],
   "source": [
    "# here is how we activate an environment in our current directory\n",
    "import Pkg; Pkg.activate(@__DIR__)\n",
    "\n",
    "# instantate this environment (download packages if you haven't)\n",
    "Pkg.instantiate();\n",
    "\n",
    "# let's load LinearAlgebra in\n",
    "using LinearAlgebra\n",
    "using Test"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "logical-footage",
   "metadata": {
    "tags": []
   },
   "source": [
    "# Question 1: Differentiation in Julia (10 pts)\n",
    "Julia has a fast and easy to use forward-mode automatic differentiation package called [ForwardDiff.jl](https://github.com/JuliaDiff/ForwardDiff.jl) that we will make use of throughout this course. In general it is easy to use and very fast, but there are a few quirks that are detailed below. This notebook will start by walking through general usage for the following cases:\n",
    "- functions with a single input \n",
    "- functions with multiple inputs\n",
    "- composite functions\n",
    "\n",
    "as well as a guide on how to avoid the most common ForwardDiff.jl error caused by creating arrays inside the function being differentiated. First, let's look at the ForwardDiff.jl functions that we are going to use:\n",
    "- `FD.derivative(f,x)` derivative of scalar or vector valued f wrt scalar x \n",
    "- `FD.jacobian(f,x)` jacobian of vector valued f wrt vector x\n",
    "- `FD.gradient(f,x)` gradient of scalar valued f wrt vector x \n",
    "- `FD.hessian(f,x)` hessian of scalar valued f wrt vector x "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "confused-wheel",
   "metadata": {},
   "source": [
    "### Note on gradients:\n",
    "For an arbitrary function $f(x):\\mathbb{R}^N \\rightarrow \\mathbb{R}^M$, the jacobian is the following:\n",
    "\n",
    "\n",
    "$$\\frac{\\partial f(x)}{\\partial x} = \\left[\\begin{array}{ccc}\n",
    "\\frac{\\partial f_{1}}{\\partial x_{1}} & \\cdots & \\frac{\\partial f_{1}}{\\partial x_{n}} \\\\\n",
    "\\vdots & \\ddots & \\vdots \\\\\n",
    "\\frac{\\partial f_{m}}{\\partial x_{1}} & \\cdots & \\frac{\\partial f_{m}}{\\partial x_{n}}\n",
    "\\end{array}\\right] $$\n",
    "\n",
    "\n",
    "Now if we have a scalar valued function (like a cost function) $f(x):\\mathbb{R}^N \\rightarrow \\mathbb{R}$, the jacobian is the following row vector:\n",
    "\n",
    "$$\\frac{\\partial f(x)}{\\partial x} = \\left[\\begin{array}{ccc}\n",
    "\\frac{\\partial f_{1}}{\\partial x_{1}} & \\cdots & \\frac{\\partial f_{1}}{\\partial x_{n}}\n",
    "\\end{array}\\right] $$\n",
    "\n",
    "The transpose of this jacobian for scalar valued functions is called the gradient:\n",
    "\n",
    "$$ \\nabla f(x) = \\bigg[\\frac{\\partial f(x)}{\\partial x}\\bigg]^T $$\n",
    "\n",
    "TLDR:\n",
    "- the jacobian of a scalar value function is a row vector \n",
    "- the gradient is the transpose of this jacobian, making the gradient a column vector \n",
    "- ForwardDiff.jl will give you an error if you try to take a jacobian of a scalar valued function, use the gradient function instead"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "nearby-blond",
   "metadata": {},
   "source": [
    "## Part (a): General usage (2 pts)\n",
    "The API for functions with one input is detailed below:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "challenging-punishment",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "∂foo1_∂x = FD.derivative(foo1, x1) = -0.5995527297708665\n",
      "∇foo2 = FD.gradient(foo2, x2) = [-0.20062732761046076, -0.20258465256181787]\n",
      "∇²foo2 = FD.hessian(foo2, x2) = [0.9796676351783216 0.0; 0.0 -0.9792647540611311]\n",
      "∂foo3_∂x = FD.jacobian(foo3, x2) = [-0.04092728957750974 -0.9796676351783216; 0.9792647540611311 -0.913735943846448]\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "2×2 Matrix{Float64}:\n",
       " -0.0409273  -0.979668\n",
       "  0.979265   -0.913736"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# NOTE: this block is a tutorial, you do not have to fill anything out.\n",
    "\n",
    "#---------load the package-----------\n",
    "# using ForwardDiff # this puts all exported functions into our namespace \n",
    "# import ForwardDiff # this means we have to use ForwardDiff.<function name>\n",
    "import ForwardDiff as FD # this let's us do FD.<function name> \n",
    "\n",
    "function foo1(x)\n",
    "    #scalar input, scalar output\n",
    "    return sin(x)*cos(x)^2\n",
    "end\n",
    "\n",
    "function foo2(x)\n",
    "    # vector input, scalar output\n",
    "    return sin(x[1]) + cos(x[2])\n",
    "end\n",
    "function foo3(x)\n",
    "    # vector input, vector output\n",
    "    return [sin(x[1])*x[2];cos(x[2])*x[1]]\n",
    "end\n",
    "\n",
    "\n",
    "let # we just use this to avoid creating global variables\n",
    "    \n",
    "    # evaluate the derivative of foo1 at x1\n",
    "    x1 = 5*randn();\n",
    "    @show ∂foo1_∂x = FD.derivative(foo1, x1);\n",
    "    \n",
    "    # evaluate the gradient and hessian of foo2 at x2\n",
    "    x2 = 5*randn(2);\n",
    "    @show ∇foo2 = FD.gradient(foo2, x2);\n",
    "    @show ∇²foo2 = FD.hessian(foo2, x2);\n",
    "    \n",
    "    # evluate the jacobian of foo3 at x2\n",
    "    @show ∂foo3_∂x = FD.jacobian(foo3,x2);\n",
    "    \n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "maritime-acceptance",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "foo4_expansion (generic function with 1 method)"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# here is our function of interest\n",
    "function foo4(x)\n",
    "    Q = diagm([1;2;3.0]) # this creates a diagonal matrix from a vector\n",
    "    return 0.5*x'*Q*x/x[1] - log(x[1])*exp(x[2])^x[3] \n",
    "end\n",
    "\n",
    "function foo4_expansion(x)\n",
    "    # TODO: this function should output the hessian H and gradient g of the function foo4\n",
    "    \n",
    "    # TODO: calculate the gradient of foo4 evaluated at x\n",
    "    g = zeros(length(x))\n",
    "    \n",
    "    # TODO: calculate the hessian of foo4 evaluated at x\n",
    "    H = zeros(length(x),length(x))\n",
    "    \n",
    "    return g, H\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "weighted-sucking",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\u001b[0m\u001b[1mTest Summary: | \u001b[22m\u001b[32m\u001b[1mPass  \u001b[22m\u001b[39m\u001b[36m\u001b[1mTotal\u001b[22m\u001b[39m\n",
      "1a            | \u001b[32m   2  \u001b[39m\u001b[36m    2\u001b[39m\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "Test.DefaultTestSet(\"1a\", Any[], 2, false, false)"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "@testset \"1a\" begin                        \n",
    "    x = [.2;.4;.5]\n",
    "    g,H = foo4_expansion(x)\n",
    "    @test isapprox(g,[-18.98201379080085, 4.982885952667278, 8.286308762133823],atol = 1e-8)        \n",
    "    @test norm(H -[164.2850689540042 -23.053506895400425 -39.942805516320334;\n",
    "                             -23.053506895400425 10.491442976333639 2.3589262864014673;\n",
    "                             -39.94280551632034 2.3589262864014673 15.314523504853529]) < 1e-8 \n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "attractive-shade",
   "metadata": {},
   "source": [
    "## Part (b): Derivatives for functions with multiple input arguments (2 pts)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "removed-sponsorship",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\u001b[32m\u001b[1mTest Passed\u001b[22m\u001b[39m"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# NOTE: this block is a tutorial, you do not have to fill anything out. \n",
    "\n",
    "# calculate derivatives for functions with multiple inputs \n",
    "function dynamics(x,a,b,c)\n",
    "    return [x[1]*a; b*c*x[2]*x[1]]\n",
    "end\n",
    "\n",
    "let \n",
    "    x1 = randn(2)\n",
    "    a = randn()\n",
    "    b = randn()\n",
    "    c = randn()\n",
    "    \n",
    "    # this evaluates the jacobian with respect to x, given a, b, and c\n",
    "    A1 = FD.jacobian(dx -> dynamics(dx, a, b, c), x1)\n",
    "    \n",
    "    # it doesn't matter what we call the new variable\n",
    "    A2 = FD.jacobian(_x -> dynamics(_x, a, b, c), x1)\n",
    "    \n",
    "    # alternatively we can do it like this using a closure\n",
    "    dynamics_just_x(_x) = dynamics(_x, a, b, c) \n",
    "    A3 = FD.jacobian(dynamics_just_x, x1)\n",
    "    \n",
    "    @test norm(A1 - A2) < 1e-13 \n",
    "    @test norm(A1 - A3) < 1e-13\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "superb-aspect",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "eulers_jacobians (generic function with 1 method)"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function eulers(x,u,J)\n",
    "    # dynamics when x is angular velocity and u is an input torque\n",
    "    ẋ = J\\(u - cross(x,J*x))\n",
    "    return ẋ\n",
    "end\n",
    "\n",
    "function eulers_jacobians(x,u,J)\n",
    "    # given x, u, and J, calculate the following two jacobians \n",
    "    \n",
    "    # TODO: fill in the following two jacobians\n",
    "    \n",
    "    # ∂ẋ/∂x\n",
    "    A = zeros(3,3)\n",
    "    \n",
    "    # ∂ẋ/∂u\n",
    "    B = zeros(3,3)\n",
    "    \n",
    "    return A, B\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "expensive-spencer",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\u001b[0m\u001b[1mTest Summary: | \u001b[22m\u001b[32m\u001b[1mPass  \u001b[22m\u001b[39m\u001b[36m\u001b[1mTotal\u001b[22m\u001b[39m\n",
      "1b            | \u001b[32m   2  \u001b[39m\u001b[36m    2\u001b[39m\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "Test.DefaultTestSet(\"1b\", Any[], 2, false, false)"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "@testset \"1b\" begin                                                \n",
    "    \n",
    "    x = [.2;-7;.2]\n",
    "    u = [.1;-.2;.343]\n",
    "    J = diagm([1.03;4;3.45])\n",
    "    \n",
    "    A,B = eulers_jacobians(x,u,J)\n",
    "\n",
    "    skew(v) = [0 -v[3] v[2]; v[3] 0 -v[1]; -v[2] v[1] 0]\n",
    "    @test isapprox(A,-J\\(skew(x)*J - skew(J*x)), atol = 1e-8)  \n",
    "\n",
    "    @test norm(B - inv(J)) < 1e-8                \n",
    "\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "public-technique",
   "metadata": {},
   "source": [
    "## Part (c): Derivatives of composite functions (1 pts)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "dress-listening",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "norm(∇f_1 - ∇f_2) = 0.0\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "0.0"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# NOTE: this block is a tutorial, you do not have to fill anything out. \n",
    "function f(x)\n",
    "    return x[1]*x[2]\n",
    "end\n",
    "function g(x)\n",
    "    return [x[1]^2; x[2]^3]\n",
    "end\n",
    "\n",
    "let \n",
    "    x1 = 2*randn(2)\n",
    "    \n",
    "    # using gradient of the composite function\n",
    "    ∇f_1 = FD.gradient(dx -> f(g(dx)), x1)\n",
    "    \n",
    "    # using the chain rule \n",
    "    J = FD.jacobian(g, x1)\n",
    "    ∇f_2 = J'*FD.gradient(f, g(x1))\n",
    "    \n",
    "    @show norm(∇f_1 - ∇f_2)\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "inappropriate-sweet",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "composite_derivs (generic function with 1 method)"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function f2(x)\n",
    "    return x*sin(x)/2\n",
    "end\n",
    "function g2(x)\n",
    "    return cos(x)^2 - tan(x)^3\n",
    "end\n",
    "\n",
    "function composite_derivs(x)\n",
    "    \n",
    "    # TODO: return ∂y/∂x where y = g2(f2(x)) \n",
    "    # (hint: this is 1D input and 1D output, so it's ForwardDiff.derivative)\n",
    "    return 0.0\n",
    "end    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "analyzed-press",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\u001b[0m\u001b[1mTest Summary: | \u001b[22m\u001b[32m\u001b[1mPass  \u001b[22m\u001b[39m\u001b[36m\u001b[1mTotal\u001b[22m\u001b[39m\n",
      "1c            | \u001b[32m   1  \u001b[39m\u001b[36m    1\u001b[39m\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "Test.DefaultTestSet(\"1c\", Any[], 1, false, false)"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "@testset \"1c\" begin                                           \n",
    "    x = 1.34 \n",
    "    deriv = composite_derivs(x)\n",
    "\n",
    "    @test isapprox(deriv,-2.390628273373545,atol = 1e-8)  \n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "lightweight-injection",
   "metadata": {},
   "source": [
    "## Part (d): Fixing the most common ForwardDiff error (2 pt)\n",
    "First we will show an example of this error:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "attended-master",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\u001b[36m\u001b[1m[ \u001b[22m\u001b[39m\u001b[36m\u001b[1mInfo: \u001b[22m\u001b[39mthis error is expected:\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-------types of input x---------\n",
      "typeof(x) = Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(f_zero_1), Float64}, Float64, 2}}\n",
      "eltype(x) = ForwardDiff.Dual{ForwardDiff.Tag{typeof(f_zero_1), Float64}, Float64, 2}\n",
      "-------types of output xdot---------\n",
      "typeof(xdot) = Vector{Float64}\n",
      "eltype(xdot) = Float64\n",
      "MethodError: no method matching Float64(::ForwardDiff.Dual{ForwardDiff.Tag{typeof(f_zero_1), Float64}, Float64, 2})\n",
      "\u001b[0mClosest candidates are:\n",
      "\u001b[0m  (::Type{T})(::Real, \u001b[91m::RoundingMode\u001b[39m) where T<:AbstractFloat at rounding.jl:200\n",
      "\u001b[0m  (::Type{T})(::T) where T<:Number at boot.jl:760\n",
      "\u001b[0m  (::Type{T})(\u001b[91m::AbstractChar\u001b[39m) where T<:Union{AbstractChar, Number} at char.jl:50\n",
      "\u001b[0m  ..."
     ]
    }
   ],
   "source": [
    "# NOTE: this block is a tutorial, you do not have to fill anything out. \n",
    "function f_zero_1(x)\n",
    "    println(\"-------types of input x---------\")\n",
    "    @show typeof(x) # print out type of x \n",
    "    @show eltype(x) # print out the element type of x \n",
    "    \n",
    "    \n",
    "    xdot = zeros(length(x)) # this default creates zeros of type Float64\n",
    "    println(\"-------types of output xdot---------\")\n",
    "    @show typeof(xdot)\n",
    "    @show eltype(xdot)\n",
    "    \n",
    "    # these lines will error because i'm trying to put a ForwardDiff.dual \n",
    "    # inside of a Vector{Float64}\n",
    "    xdot[1] = x[1]*x[2]\n",
    "    xdot[2] = x[2]^2\n",
    "    \n",
    "    return xdot \n",
    "end\n",
    "\n",
    "let \n",
    "    # try and calculate the jacobian of f_zero_1 on x1\n",
    "    x1 = randn(2)\n",
    "    @info \"this error is expected:\"\n",
    "    try \n",
    "        FD.jacobian(f_zero_1,x1)\n",
    "    catch e \n",
    "        buf = IOBuffer()\n",
    "        showerror(buf,e)\n",
    "        message = String(take!(buf))\n",
    "        Base.showerror(stdout,e)\n",
    "    end\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "realistic-afghanistan",
   "metadata": {},
   "source": [
    "This is the most common ForwardDiff error that you will encounter. ForwardDiff works by pushing `ForwardDiff.Dual` variables through the function being differentiated. Normally this works without issue, but if you create a vector of `Float64` (like you would with `xdot = zeros(5)`, it is unable to fit the `ForwardDiff.Dual`'s in with the `Float64`'s. To get around this, you have two options:"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "disturbed-oregon",
   "metadata": {},
   "source": [
    "### Option 1 \n",
    "Our first option is just creating xdot directly, without creating an array of zeros to index into. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "preceding-destiny",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-------types of input x---------\n",
      "typeof(x) = Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(f_zero_1), Float64}, Float64, 2}}\n",
      "eltype(x) = ForwardDiff.Dual{ForwardDiff.Tag{typeof(f_zero_1), Float64}, Float64, 2}\n",
      "-------types of output xdot---------\n",
      "typeof(xdot) = Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(f_zero_1), Float64}, Float64, 2}}\n",
      "eltype(xdot) = ForwardDiff.Dual{ForwardDiff.Tag{typeof(f_zero_1), Float64}, Float64, 2}\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "2×2 Matrix{Float64}:\n",
       " -0.322914  -0.376712\n",
       " -0.0       -0.645828"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# NOTE: this block is a tutorial, you do not have to fill anything out. \n",
    "function f_zero_1(x)\n",
    "    \n",
    "    # let's create xdot directly, without first making a vector of zeros \n",
    "    xdot = [x[1]*x[2], x[2]^2]\n",
    "    \n",
    "    # NOTE: the compiler figures out which type to make xdot, so when you call the function normally\n",
    "    # it's a Float64, and when it's being diffed, it's automatically promoted to a ForwardDiff.Dual type\n",
    "    \n",
    "    println(\"-------types of input x---------\")\n",
    "    @show typeof(x) # print out type of x \n",
    "    @show eltype(x) # print out the element type of x \n",
    "    \n",
    "    println(\"-------types of output xdot---------\")\n",
    "    @show typeof(xdot)\n",
    "    @show eltype(xdot)\n",
    "    \n",
    "    return xdot \n",
    "end\n",
    "\n",
    "let \n",
    "    # try and calculate the jacobian of f_zero_1 on x1\n",
    "    x1 = randn(2)\n",
    "    FD.jacobian(f_zero_1,x1) # this will work\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "mathematical-cause",
   "metadata": {},
   "source": [
    "### Option 2\n",
    "The second option is to create the array of zeros in a way that accounts for the input type. This can be done by replacing `zeros(length(x))` with `zeros(eltype(x),length(x))`. The first argument `eltype(x)` simply creates a vector of zeros that is the same type as the element type in vector x. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "hungarian-calcium",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-------types of input x---------\n",
      "typeof(x) = Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(f_zero_1), Float64}, Float64, 2}}\n",
      "eltype(x) = ForwardDiff.Dual{ForwardDiff.Tag{typeof(f_zero_1), Float64}, Float64, 2}\n",
      "-------types of output xdot---------\n",
      "typeof(xdot) = Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(f_zero_1), Float64}, Float64, 2}}\n",
      "eltype(xdot) = ForwardDiff.Dual{ForwardDiff.Tag{typeof(f_zero_1), Float64}, Float64, 2}\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "2×2 Matrix{Float64}:\n",
       " 1.01609  -0.635973\n",
       " 0.0       2.03218"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# NOTE: this block is a tutorial, you do not have to fill anything out. \n",
    "function f_zero_1(x)\n",
    "    \n",
    "    xdot = zeros(eltype(x), length(x))\n",
    "    \n",
    "    xdot[1] = x[1]*x[2]\n",
    "    xdot[2] = x[2]^2\n",
    "    \n",
    "    println(\"-------types of input x---------\")\n",
    "    @show typeof(x) # print out type of x \n",
    "    @show eltype(x) # print out the element type of x \n",
    "    \n",
    "    println(\"-------types of output xdot---------\")\n",
    "    @show typeof(xdot)\n",
    "    @show eltype(xdot)\n",
    "    \n",
    "    return xdot \n",
    "end\n",
    "\n",
    "let \n",
    "    # try and calculate the jacobian of f_zero_1 on x1\n",
    "    x1 = randn(2)\n",
    "    FD.jacobian(f_zero_1,x1) # this will fail! \n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "creative-transport",
   "metadata": {},
   "source": [
    "Now you can show that you understand these two options by fixing two broken functions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "sixth-education",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "dynamics (generic function with 2 methods)"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# TODO: fix this error when trying to diff through this function\n",
    "# hint: you can use promote_type(eltype(x),eltype(u)) to return the correct type if either x or u is a ForwardDiff.Dual (option 1)\n",
    "\n",
    "function dynamics(x,u)\n",
    "    xdot = zeros(length(x))\n",
    "    xdot[1] = x[1]*sin(u[1])\n",
    "    xdot[2] = x[2]*cos(u[2])\n",
    "    return xdot\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "suspected-lancaster",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\u001b[0m\u001b[1mTest Summary: | \u001b[22m\u001b[32m\u001b[1mPass  \u001b[22m\u001b[39m\u001b[36m\u001b[1mTotal\u001b[22m\u001b[39m\n",
      "1d            | \u001b[32m   2  \u001b[39m\u001b[36m    2\u001b[39m\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "Test.DefaultTestSet(\"1d\", Any[], 2, false, false)"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "@testset \"1d\" begin                                     \n",
    "    x = [.1;.4]\n",
    "    u = [.2;-.3]\n",
    "    A = FD.jacobian(_x -> dynamics(_x,u),x) \n",
    "    B = FD.jacobian(_u -> dynamics(x,_u),u) \n",
    "    @test typeof(A) == Matrix{Float64}                  \n",
    "    @test typeof(B) == Matrix{Float64}                  \n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e1149952",
   "metadata": {},
   "source": [
    "## Finite Difference Derivatives \n",
    "If you ever have trouble working through a ForwardDiff error, you should always feel free to use the [FiniteDiff.jl](https://github.com/JuliaDiff/FiniteDiff.jl) FiniteDiff.jl package instead. This computes derivatives through a [finite difference method](https://en.wikipedia.org/wiki/Finite_difference_method). This is slower and less accurate than ForwardDiff, but it will always work so long as the function works.  \n",
    "\n",
    "Before with ForwardDiff we had this:\n",
    "\n",
    "- `FD.derivative(f,x)` derivative of scalar or vector valued f wrt scalar x \n",
    "- `FD.jacobian(f,x)` jacobian of vector valued f wrt vector x\n",
    "- `FD.gradient(f,x)` gradient of scalar valued f wrt vector x \n",
    "- `FD.hessian(f,x)` hessian of scalar valued f wrt vector x \n",
    "\n",
    "Now with FiniteDiff we have this:\n",
    "\n",
    "- `FD2.finite_difference_derivative(f,x)` derivative of scalar or vector valued f wrt scalar x \n",
    "- `FD2.finite_difference_jacobian(f,x)` jacobian of vector valued f wrt vector x\n",
    "- `FD2.finite_difference_gradient(f,x)` gradient of scalar valued f wrt vector x \n",
    "- `FD2.finite_difference_hessian(f,x)` hessian of scalar valued f wrt vector x "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "id": "67e88724",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "∂foo1_∂x = FD2.finite_difference_derivative(foo1, x1) = 0.3341157713242085\n",
      "∇foo2 = FD2.finite_difference_gradient(foo2, x2) = [0.885051153162438, 0.08516940844437484]\n",
      "∇²foo2 = FD2.finite_difference_hessian(foo2, x2) = [-0.4654937833547592 0.0; 0.0 -0.996366485953331]\n",
      "∂foo3_∂x = FD2.finite_difference_jacobian(foo3, x2) = [-0.07547071296721697 0.4654937768355012; 0.996366485953331 0.04123838245868683]\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "\u001b[32m\u001b[1mTest Passed\u001b[22m\u001b[39m"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# NOTE: this block is a tutorial, you do not have to fill anything out.\n",
    "\n",
    "# load the package \n",
    "import FiniteDiff as FD2 \n",
    "\n",
    "function foo1(x)\n",
    "    #scalar input, scalar output\n",
    "    return sin(x)*cos(x)^2\n",
    "end\n",
    "\n",
    "function foo2(x)\n",
    "    # vector input, scalar output\n",
    "    return sin(x[1]) + cos(x[2])\n",
    "end\n",
    "function foo3(x)\n",
    "    # vector input, vector output\n",
    "    return [sin(x[1])*x[2];cos(x[2])*x[1]]\n",
    "end\n",
    "\n",
    "\n",
    "let # we just use this to avoid creating global variables\n",
    "    \n",
    "    # evaluate the derivative of foo1 at x1\n",
    "    x1 = 5*randn();\n",
    "    @show ∂foo1_∂x = FD2.finite_difference_derivative(foo1, x1);\n",
    "    \n",
    "    # evaluate the gradient and hessian of foo2 at x2\n",
    "    x2 = 5*randn(2);\n",
    "    @show ∇foo2 = FD2.finite_difference_gradient(foo2, x2);\n",
    "    @show ∇²foo2 = FD2.finite_difference_hessian(foo2, x2);\n",
    "    \n",
    "    # evluate the jacobian of foo3 at x2\n",
    "    @show ∂foo3_∂x = FD2.finite_difference_jacobian(foo3,x2);\n",
    "    \n",
    "    @test norm(∂foo1_∂x - FD.derivative(foo1, x1)) < 1e-4 \n",
    "    @test norm(∇foo2 - FD.gradient(foo2, x2)) < 1e-4 \n",
    "    @test norm(∇²foo2 - FD.hessian(foo2, x2)) < 1e-4 \n",
    "    @test norm(∂foo3_∂x - FD.jacobian(foo3, x2)) < 1e-4 \n",
    "    \n",
    "    \n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fc24d707",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 1.6.7",
   "language": "julia",
   "name": "julia-1.6"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.6.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}