bracketing.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Bracketing\n",
    "\n",
    "This notebook was automatically generated from the Algorithms for Optimization source code. Each cell generates a figure from the original text. While this code is not optimized for use in lectures, we provide it here to be adapted for such projects. We hope you find it useful."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "include(\"support_code.jl\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "p = let\n",
    "\tN = 4\n",
    "\tg = GroupPlot(N,1,groupStyle=\"horizontal sep=0.25cm, ylabels at=edge left\",\n",
    "\t                  style=\"width=5cm, xlabel=\\$x\\$, ylabel=\\$y\\$, ytick=\\\\empty\")\n",
    "\n",
    "\tf = x->sin(x)\n",
    "\n",
    "\ts = 0.5\n",
    "\tk = 2.0\n",
    "\n",
    "\ta = 7.0\n",
    "\tb = a + s\n",
    "\tya = f(a)\n",
    "\tyb = f(b)\n",
    "\n",
    "\tx_arr = collect(range(3π/2-2, stop=3π/2+5.5, length=101))\n",
    "\n",
    "\tfor i in 1 : N\n",
    "\n",
    "\t    ax = Axis([Plots.Linear(x_arr, f.(x_arr), style=\"solid, black, mark=none\"),\n",
    "\t          Plots.Scatter([a,b], [f(a), f(b)], style=\"black, mark=*, mark options={solid, draw=black, fill=black}\"),\n",
    "\t            ], style=\"xtick={$a,$b}, xticklabels={\\$a\\$,\\$b\\$}, xticklabel style={text height=2ex}\")\n",
    "\n",
    "\t    push!(g, ax)\n",
    "\n",
    "\t    if yb > ya\n",
    "\t        a, b = b, a\n",
    "\t        ya, yb = yb, ya\n",
    "\t        s = -s\n",
    "\t    end\n",
    "\n",
    "\t    c, yc = b + s, f(b + s)\n",
    "\t    if yc > yb\n",
    "\t        a, b = a, c\n",
    "\t    else\n",
    "\t        a, ya, b, yb = b, yb, c, yc\n",
    "\t        s *= k\n",
    "\t    end\n",
    "\tend\n",
    "\n",
    "\tg\n",
    "end\n",
    "\n",
    "plot(p)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "φ = MathConstants.φ\n",
    "p = let\n",
    "\n",
    "\tf = x-> -exp(-x^2)\n",
    "\ta₀, b₀ = -2.0, 2.0\n",
    "\tlo, hi = -1.0, 0.25\n",
    "\n",
    "\tN = 5\n",
    "\tg = GroupPlot(N,2,groupStyle=\"horizontal sep=0.25cm, vertical sep=1.25cm, ylabels at=edge left\",\n",
    "\t                  style=\"clip marker paths, width=4.25cm, ymin=$lo, ymax=$hi, xmin=$a₀, xmax=$b₀, xlabel=\\$x\\$, ylabel=\\$y\\$, ytick=\\\\empty,  xtick=\\\\empty, axis on top=true\")\n",
    "\n",
    "\ta, b = a₀, b₀\n",
    "\tn = N+1\n",
    "\ts = (1-√5)/(1+√5)\n",
    "    ρ = 1 / (φ*(1-s^(n+1))/(1-s^n))\n",
    "    d = ρ*b + (1 - ρ)*a\n",
    "    fd = f(d)\n",
    "    for i in 1 : n - 1\n",
    "        if i == n - 1\n",
    "            c = 0.01*a + 0.99*d\n",
    "        else\n",
    "            c = ρ*a + (1 - ρ)*b\n",
    "        end\n",
    "        fc = f(c)\n",
    "        if fc < fd\n",
    "            b, d, fd = d, c, fc\n",
    "        else\n",
    "            a, b, fb = b, c, fc\n",
    "        end\n",
    "        ρ = 1 / (φ*(1-s^(n-i+1))/(1-s^(n-i)))\n",
    "\n",
    "        p = Plots.Plot[]\n",
    "        push!(p, Plots.Linear(f, (a₀, b₀), style=\"solid, thick, black, mark=none\"))\n",
    "        push!(p, Plots.Linear([a,a], [lo,hi], style=\"name path=A, draw=none, mark=none\"))\n",
    "        push!(p, Plots.Linear([b,b], [lo,hi], style=\"name path=B, draw=none, mark=none\"))\n",
    "        push!(p, Plots.Command(\"\\\\addplot[pastelBlue!50] fill between[of=A and B]\"))\n",
    "        push!(p, Plots.Scatter([a,b], [f(a), f(b)], style=\"mark=*, mark size=1, mark options={solid, draw=black, fill=black}\"))\n",
    "        push!(g, Axis(p))\n",
    "    end\n",
    "\n",
    "\t###\n",
    "\n",
    "\tfunction golden_section_search(f, a, b, n)\n",
    "\t    ρ = φ-1\n",
    "\t    d = ρ * b + (1 - ρ)*a\n",
    "\t    fd = f(d)\n",
    "\t    for i = 1 : n-1\n",
    "\t        c = ρ*a + (1 - ρ)*b\n",
    "\t        fc = f(c)\n",
    "\t        if fc < fd\n",
    "\t            b, d, fd = d, c, fc\n",
    "\t        else\n",
    "\t            a, b, fb = b, c, fc\n",
    "\t        end\n",
    "\t    end\n",
    "\t    return a < b ? (a, b) : (b, a)\n",
    "\tend\n",
    "\n",
    "\n",
    "\tfor n in 2 : N+1\n",
    "\t    a, b = golden_section_search(f, a₀, b₀, n)\n",
    "\n",
    "\t    p = Plots.Plot[]\n",
    "\t    push!(p, Plots.Linear(f, (a₀, b₀), style=\"solid, thick, black, mark=none\"))\n",
    "\t    push!(p, Plots.Linear([a,a], [lo,hi], style=\"name path=A, draw=none, mark=none\"))\n",
    "\t    push!(p, Plots.Linear([b,b], [lo,hi], style=\"name path=B, draw=none, mark=none\"))\n",
    "\t    push!(p, Plots.Command(\"\\\\addplot[pastelBlue!50] fill between[of=A and B]\"))\n",
    "\t    push!(p, Plots.Scatter([a,b], [f(a), f(b)], style=\"mark=*, mark size=1, mark options={solid, draw=black, fill=black}\"))\n",
    "\t    push!(g, Axis(p))\n",
    "\tend\n",
    "\n",
    "\tg.axes[3].title = \"Fibonacci Search\"\n",
    "\tg.axes[8].title = \"Golden Section Search\"\n",
    "\n",
    "\tg\n",
    "end\n",
    "\n",
    "plot(p)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "\n",
    "p = let\n",
    "\tf = x-> (sin(x) + sin(x/2))/4\n",
    "\ta₀, b₀ = -5.0, 9.0\n",
    "\tlo, hi = -0.6, 0.6\n",
    "\n",
    "\tN = 5\n",
    "\tg = GroupPlot(N,2,groupStyle=\"horizontal sep=0.25cm, vertical sep=1.25cm, ylabels at=edge left\",\n",
    "\t                  style=\"clip marker paths, width=4.25cm, ymin=$lo, ymax=$hi, xmin=$a₀, xmax=$b₀, xlabel=\\$x\\$, ylabel=\\$y\\$, ytick=\\\\empty,  xtick=\\\\empty, axis on top=true\")\n",
    "\n",
    "\ta, b = a₀, b₀\n",
    "\tn = N+1\n",
    "\ts = (1-√5)/(1+√5)\n",
    "    ρ = 1 / (φ*(1-s^(n+1))/(1-s^n))\n",
    "    d = ρ*b + (1 - ρ)*a\n",
    "    fd = f(d)\n",
    "    for i in 1 : n - 1\n",
    "        if i == n - 1\n",
    "            c = 0.01*a + 0.99*d\n",
    "        else\n",
    "            c = ρ*a + (1 - ρ)*b\n",
    "        end\n",
    "        fc = f(c)\n",
    "        if fc < fd\n",
    "            b, d, fd = d, c, fc\n",
    "        else\n",
    "            a, b, fb = b, c, fc\n",
    "        end\n",
    "        ρ = 1 / (φ*(1-s^(n-i+1))/(1-s^(n-i)))\n",
    "\n",
    "        p = Plots.Plot[]\n",
    "        push!(p, Plots.Linear(f, (a₀, b₀), style=\"solid, thick, black, mark=none\"))\n",
    "        push!(p, Plots.Linear([a,a], [lo,hi], style=\"name path=A, draw=none, mark=none\"))\n",
    "        push!(p, Plots.Linear([b,b], [lo,hi], style=\"name path=B, draw=none, mark=none\"))\n",
    "        push!(p, Plots.Command(\"\\\\addplot[pastelBlue!50] fill between[of=A and B]\"))\n",
    "        push!(p, Plots.Scatter([a,b], [f(a), f(b)], style=\"mark=*, mark size=1, mark options={solid, draw=black, fill=black}\"))\n",
    "        push!(g, Axis(p))\n",
    "    end\n",
    "\n",
    "\t###\n",
    "\n",
    "\tfunction golden_section_search(f, a, b, n)\n",
    "\t    ρ = φ-1\n",
    "\t    d = ρ * b + (1 - ρ)*a\n",
    "\t    fd = f(d)\n",
    "\t    for i = 1 : n-1\n",
    "\t        c = ρ*a + (1 - ρ)*b\n",
    "\t        fc = f(c)\n",
    "\t        if fc < fd\n",
    "\t            b, d, fd = d, c, fc\n",
    "\t        else\n",
    "\t            a, b, fb = b, c, fc\n",
    "\t        end\n",
    "\t    end\n",
    "\t    return a < b ? (a, b) : (b, a)\n",
    "\tend\n",
    "\n",
    "\n",
    "\tfor n in 2 : N+1\n",
    "\t    a, b = golden_section_search(f, a₀, b₀, n)\n",
    "\n",
    "\t    p = Plots.Plot[]\n",
    "\t    push!(p, Plots.Linear(f, (a₀, b₀), style=\"solid, thick, black, mark=none\"))\n",
    "\t    push!(p, Plots.Linear([a,a], [lo,hi], style=\"name path=A, draw=none, mark=none\"))\n",
    "\t    push!(p, Plots.Linear([b,b], [lo,hi], style=\"name path=B, draw=none, mark=none\"))\n",
    "\t    push!(p, Plots.Command(\"\\\\addplot[pastelBlue!50] fill between[of=A and B]\"))\n",
    "\t    push!(p, Plots.Scatter([a,b], [f(a), f(b)], style=\"mark=*, mark size=1, mark options={solid, draw=black, fill=black}\"))\n",
    "\t    push!(g, Axis(p))\n",
    "\tend\n",
    "\n",
    "\tg.axes[3].title = \"Fibonacci Search\"\n",
    "\tg.axes[8].title = \"Golden Section Search\"\n",
    "\n",
    "\tg\n",
    "end\n",
    "\n",
    "plot(p)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "\n",
    "   p = let\n",
    "\n",
    "   \t\tf = x->(x-1.25)^2 + 0.25\n",
    "   \t\txpts = [0.3,1.0,1.75]\n",
    "\n",
    "    \tAxis([\n",
    "    \t\t\tPlots.Linear(f, (0.25,2), style=\"solid, black, mark=none\"),\n",
    "    \t\t\tPlots.Scatter(xpts,f.(xpts), style=\"black, mark=*, mark options={solid, draw=black, fill=black}\"),\n",
    "    \t\t\tPlots.Scatter([1.25], [f(1.25)], style=\"pastelBlue, mark=*, mark options={solid, draw=pastelBlue, fill=pastelBlue}\"),\n",
    "    \t\t],\n",
    "    \t\twidth=\"1.2*9cm\", xlabel=L\"x\", ylabel=L\"y\", style=\"axis lines=middle, xtick=\\\\empty, ytick=\\\\empty\",\n",
    "    \t\txmin=0, xmax=2.25, ymin=0,\n",
    "    \t)\n",
    "   end\n",
    "\n",
    "   plot(p)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "p = let\n",
    "\tf = x->x - 2x^2 + 0.2x^3\n",
    "\n",
    "\tfunction get_q(a,b,c)\n",
    "\t    fa = f(a)\n",
    "\t    fb = f(b)\n",
    "\t    fc = f(c)\n",
    "\t    return x-> fa*(x-b)*(x-c)/((a-b)*(a-c)) +\n",
    "\t               fb*(x-a)*(x-c)/((b-a)*(b-c)) +\n",
    "\t               fc*(x-a)*(x-b)/((c-a)*(c-b))\n",
    "\tend\n",
    "\n",
    "\tfunction get_axis(a,b,c)\n",
    "\n",
    "\t    x = 0.5*(ya*(b^2-c^2)+yb*(c^2-a^2)+yc*(a^2-b^2)) /\n",
    "\t            (ya*(b-c)    +yb*(c-a)    +yc*(a-b))\n",
    "\t    yx = f(x)\n",
    "\n",
    "\t    Axis([Plots.Linear(f, (0,10), style=\"black, solid, mark=none\"),\n",
    "\t      Plots.Scatter([a,b,c], f.([a,b,c]), style=\"mark=*, mark size=1.2, mark options={solid, draw=black, fill=black}\"),\n",
    "\t    Plots.Scatter([x], [f(x)], style=\"mark=*, mark size=1.2, mark options={solid, draw=pastelBlue, fill=pastelBlue}\"),\n",
    "\t      Plots.Linear(get_q(a,b,c), (0,10), style=\"pastelBlue, solid, mark=none\"),\n",
    "\t      ], style=\"xtick={$a,$b,$c}, xticklabels={\\$a\\$,\\$b\\$,\\$c\\$}\") #$\n",
    "\tend\n",
    "\n",
    "\tg = GroupPlot(4,1,groupStyle=\"horizontal sep=0.25cm, ylabels at=edge left\",\n",
    "\t\t\t                    style=\"width=5cm, xlabel=\\$x\\$, ylabel=\\$y\\$, ytick=\\\\empty, ymin=-25, ymax=30, xmin=-1, xmax=11, xticklabel style={text height=2ex}\")\n",
    "\n",
    "\ta = 0.5\n",
    "\tb = 4\n",
    "\tc = 9.75\n",
    "\n",
    "\tya, yb, yc = f(a), f(b), f(c)\n",
    "\n",
    "\tfor N in 1 : 4\n",
    "\t    push!(g, get_axis(a,b,c))\n",
    "\n",
    "\t    x = 0.5*(ya*(b^2-c^2)+yb*(c^2-a^2)+yc*(a^2-b^2)) /\n",
    "\t            (ya*(b-c)    +yb*(c-a)    +yc*(a-b))\n",
    "\t    yx = f(x)\n",
    "\n",
    "\t    a, b, c,ya, yb, yc = x > b ?\n",
    "\t\t\t\t(b, x, c, yb, f(x), yc) :\n",
    "\t\t\t\t(a, x, b, ya, f(x), yb)\n",
    "\n",
    "\t    if x > b\n",
    "\t        if yx > yb\n",
    "\t            c, yc = x, yx\n",
    "\t        else\n",
    "\t            a, ya, b, yb = b, yb, x, yx\n",
    "\t        end\n",
    "\t    elseif x < b\n",
    "\t        if yx > yb\n",
    "\t            a, ya = x, yx\n",
    "\t        else\n",
    "\t            c, yc, b, yv = b, yb, x, yx\n",
    "\t        end\n",
    "\t    end\n",
    "\tend\n",
    "\n",
    "\tg.axes[end].style *= \", xticklabels={{\\\\hspace{-4pt}\\$a\\$},{\\\\hspace{4pt}\\$b\\$},\\$c\\$}\"\n",
    "\n",
    "\tg\n",
    "end\n",
    "\n",
    "plot(p)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "using Vec\n",
    "\n",
    "p = let\n",
    "\tf = x-> 1.75 - (exp(-((x-2.5)^2)) + 1.25*exp(-0.5*(x-6)^2))\n",
    "\n",
    "\ta = 0.5\n",
    "\tb = 6.5\n",
    "\n",
    "\tfunction _get_shubert_piyavskii_intersection(A::VecE2, B::VecE2, L::Float64)\n",
    "\t    t = ((A.y - B.y) - L*(A.x - B.x)) / (2L)\n",
    "\t    A + VecE2(t, -t*L)\n",
    "\tend\n",
    "\n",
    "\n",
    "\tN = 1\n",
    "\tplots = Plots.Plot[]\n",
    "\tpush!(plots, Plots.Linear(f, (a,b), xbins=151, style=\"solid, thick, black, mark=none\", legendentry=L\"f(x)\"))\n",
    "\n",
    "\tL = 1.0\n",
    "\tx1 = (a+b)/2\n",
    "\tP1 = VecE2(x1, f(x1))\n",
    "\tδ1 = (b-a)/2\n",
    "\tpts = [VecE2(a, P1.y-L*δ1), P1, VecE2(b, P1.y-L*δ1)]\n",
    "\tpush!(plots, Plots.Linear([P.x for P in pts], [P.y for P in pts], style=\"solid, thick, gray, mark=none\", legendentry=\"lower bound\"))\n",
    "\tpush!(plots, Plots.Scatter([P1.x], [P1.y], style=\"mark=*, mark size=1.5, mark options={solid, draw=pastelBlue, fill=pastelBlue}\"))\n",
    "\n",
    "\tA = P1 + VecE2(-1.0, -L*1.0)\n",
    "\th = 1.0\n",
    "\tpush!(plots, Plots.Linear([A.x, A.x, A.x-h/L], [A.y, A.y-h, A.y-h], style=\"solid, gray, mark=none\"))\n",
    "\tpush!(plots, Plots.Node(L\"\\ell\",A.x,A.y-h/2, style=\"right\"))\n",
    "\tpush!(plots, Plots.Node(L\"1\",A.x-h/(2L),A.y-h, style=\"below\"))\n",
    "\n",
    "\n",
    "\tAxis(plots, style=\"width=10cm, height=6cm, ymin=-2, ymax=2, xmin=$(a-0.5), xmax=$(b+0.5), xlabel=\\$x\\$, ylabel=\\$y\\$, xtick={$a, $((b+a)/2), $b}, xticklabels={\\$a\\$, \\$x^{(1)}\\$, \\$b\\$}, ytick=\\\\empty, legend pos=outer north east, xticklabel style={text height=2ex}\")\n",
    "end\n",
    "\n",
    "plot(p)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "using Vec\n",
    "\n",
    "p = let\n",
    "\tf = x-> 1.75 - (exp(-((x-2.5)^2)) + 1.25*exp(-0.5*(x-6)^2))\n",
    "\n",
    "\ta = 0.5\n",
    "\tb = 6.5\n",
    "\n",
    "\tfunction _get_shubert_piyavskii_intersection(A::VecE2, B::VecE2, L::Float64)\n",
    "\t    t = ((A.y - B.y) - L*(A.x - B.x)) / (2L)\n",
    "\t    A + VecE2(t, -t*L)\n",
    "\tend\n",
    "\n",
    "\n",
    "\tN = 1\n",
    "\tplots = Plots.Plot[]\n",
    "\tpush!(plots, Plots.Linear(f, (a,b), xbins=151, style=\"solid, thick, black, mark=none\", legendentry=L\"f(x)\"))\n",
    "\n",
    "\tL = 1.0\n",
    "\tx1 = (a+b)/2\n",
    "\tP1 = VecE2(x1, f(x1))\n",
    "\tδ1 = (b-a)/2\n",
    "\tpts = [VecE2(a, P1.y-L*δ1), P1, VecE2(b, P1.y-L*δ1)]\n",
    "\n",
    "\n",
    "\tN = 3\n",
    "\tfor k in 2 : N\n",
    "\t    i = argmin([P.y for P in pts])\n",
    "\t    x = pts[i].x\n",
    "\t    P = VecE2(x, f(x))\n",
    "\n",
    "\t    if i == 1\n",
    "\t        pts[1] = P\n",
    "\t        insert!(pts, 2, _get_shubert_piyavskii_intersection(P, pts[2], L))\n",
    "\t    elseif i == length(pts)\n",
    "\t        pts[end] = P\n",
    "\t        insert!(pts, length(pts), _get_shubert_piyavskii_intersection(pts[end-1], P, L))\n",
    "\t    else\n",
    "\t        P_prev = _get_shubert_piyavskii_intersection(pts[i-1], P, L)\n",
    "\t        P_next = _get_shubert_piyavskii_intersection(P, pts[i+1], L)\n",
    "\t        deleteat!(pts, i)\n",
    "\t        insert!(pts, i, P_next)\n",
    "\t        insert!(pts, i, P)\n",
    "\t        insert!(pts, i, P_prev)\n",
    "\t    end\n",
    "\tend\n",
    "\n",
    "\tpush!(plots, Plots.Linear([P.x for P in pts], [P.y for P in pts], style=\"solid, thick, gray, mark=none\", legendentry=\"lower bound\"))\n",
    "\tpush!(plots, Plots.Scatter([P.x for P in pts[1:2:end]], [P.y for P in pts[1:2:end]], style=\"mark=*, mark size=1.5, mark options={solid, draw=pastelBlue, fill=pastelBlue}\"))\n",
    "\n",
    "\ti = argmin([P.y for P in pts])\n",
    "\tx = pts[i].x\n",
    "\tP = VecE2(x, f(x))\n",
    "\tP_prev = _get_shubert_piyavskii_intersection(pts[i-1], P, L)\n",
    "\tP_next = _get_shubert_piyavskii_intersection(P, pts[i+1], L)\n",
    "\n",
    "\tpush!(plots, Plots.Linear([P_prev.x, P.x, P_next.x], [P_prev.y, P.y, P_prev.y], style=\"solid, thick, pastelBlue, mark=none\"))\n",
    "\tpush!(plots, Plots.Scatter([P.x], [P.y], style=\"mark=*, mark size=1.5, mark options={solid, draw=pastelBlue, fill=pastelBlue}\"))\n",
    "\tpush!(plots, Plots.Command(\"\\\\draw[->, ] ($x,$(pts[i].y + 0.2)) -- ($x,$(P.y - 0.2))\"))\n",
    "\n",
    "\txticks = \"$a, $(pts[3].x), $(P.x), $b\"\n",
    "\txticklabels = \"\\$a = x^{(2)}\\$, \\$x^{(1)}\\$, \\$x^{(4)}\\$, \\$b = x^{(3)}\\$\"\n",
    "\n",
    "\tAxis(plots, style=\"width=10cm, height=6cm, ymin=-0.75, ymax=2, xmin=$(a-0.5), xmax=$(b+0.5), xlabel=\\$x\\$, ylabel=\\$y\\$, xtick={$xticks}, xticklabels={$xticklabels}, ytick=\\\\empty, legend pos=outer north east, xticklabel style={text height=2ex}\")\n",
    "end\n",
    "\n",
    "plot(p)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "using Vec\n",
    "\n",
    "p = let\n",
    "\n",
    "\tf = x-> 1.75 - (exp(-((x-2.5)^2)) + 1.25*exp(-0.5*(x-6)^2))\n",
    "\n",
    "\ta = 0.5\n",
    "\tb = 6.5\n",
    "\n",
    "\tfunction _get_shubert_piyavskii_intersection(A::VecE2, B::VecE2, L::Float64)\n",
    "\t    t = ((A.y - B.y) - L*(A.x - B.x)) / (2L)\n",
    "\t    A + VecE2(t, -t*L)\n",
    "\tend\n",
    "\tfunction get_shubert_piyavskii_axis(N)\n",
    "\n",
    "\t    plots = Plots.Plot[]\n",
    "\n",
    "\t    push!(plots, Plots.Linear(f, (a,b), xbins=151, style=\"solid, thick, black, mark=none\"))\n",
    "\n",
    "\t    L = 1.0\n",
    "\t    x1 = (a+b)/2\n",
    "\t    P1 = VecE2(x1, f(x1))\n",
    "\t    δ1 = (b-a)/2\n",
    "\t    pts = [VecE2(a, P1.y-L*δ1), P1, VecE2(b, P1.y-L*δ1)]\n",
    "\n",
    "\t    for k in 2 : N\n",
    "\t        i = argmin([P.y for P in pts])\n",
    "\t        x = pts[i].x\n",
    "\t        P = VecE2(x, f(x))\n",
    "\n",
    "\t        if i == 1\n",
    "\t        \tpts[1] = P\n",
    "\t        \tinsert!(pts, 2, _get_shubert_piyavskii_intersection(P, pts[2], L))\n",
    "\t        elseif i == length(pts)\n",
    "\t        \tpts[end] = P\n",
    "\t        \tinsert!(pts, length(pts), _get_shubert_piyavskii_intersection(pts[end-1], P, L))\n",
    "\t        else\n",
    "\t        \tP_prev = _get_shubert_piyavskii_intersection(pts[i-1], P, L)\n",
    "\t        \tP_next = _get_shubert_piyavskii_intersection(P, pts[i+1], L)\n",
    "\t        \tdeleteat!(pts, i)\n",
    "\t        \tinsert!(pts, i, P_next)\n",
    "\t        \tinsert!(pts, i, P)\n",
    "\t        \tinsert!(pts, i, P_prev)\n",
    "\t       \tend\n",
    "\t    end\n",
    "\n",
    "\t    push!(plots, Plots.Linear([P.x for P in pts], [P.y for P in pts], style=\"solid, gray, mark=none\"))\n",
    "\n",
    "\t    # now plot the intersections\n",
    "\t    offset = pts[2].y < pts[1].y ? 1 : 2\n",
    "\n",
    "\t    i = offset + 2*(argmin([P.y for P in pts[offset:2:end]]) - 1)\n",
    "\t    xmin = pts[i].x\n",
    "\t    ymin = pts[i].y\n",
    "\n",
    "\t    for j in (offset == 1 ? 2 : 1) : 2 : length(pts)\n",
    "\t        if pts[j].y < ymin\n",
    "\t            xj = pts[j].x\n",
    "\t            delta = ymin - pts[j].y\n",
    "\t            x_lo = max(a, xj-delta/L)\n",
    "\t            x_hi = min(b, xj+delta/L)\n",
    "\t            push!(plots, Plots.Linear([x_lo,x_hi],[ymin,ymin], style=\"pastelBlue, ultra thick, solid, mark=none\"))\n",
    "\t        end\n",
    "\t    end\n",
    "\n",
    "\t    Axis(plots)\n",
    "\tend\n",
    "\n",
    "\tg = GroupPlot(3,3,groupStyle=\"horizontal sep=0.25cm, vertical sep=0.25cm, xlabels at=edge bottom, ylabels at=edge left\",\n",
    "\t                    style=\"width=5cm, ymin=0, ymax=2, xmin=$(a-0.5), xmax=$(b+0.5), xlabel=\\$x\\$, ylabel=\\$y\\$, xtick=\\\\empty, ytick=\\\\empty\")\n",
    "\n",
    "\tfor N in 1 : 9\n",
    "\t    push!(g, get_shubert_piyavskii_axis(N))\n",
    "\tend\n",
    "\n",
    "\tg\n",
    "end\n",
    "\n",
    "plot(p)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "p = let\n",
    "\tN = 4\n",
    "\n",
    "\tf = x-> (sin(x) + sin(x/2))/4\n",
    "\ta₀, b₀ = -5.0, 9.0\n",
    "\ta, b = a₀, b₀\n",
    "\tya = f(a)\n",
    "\tyb = f(b)\n",
    "\n",
    "\tg = GroupPlot(N,1,groupStyle=\"horizontal sep=0.25cm, ylabels at=edge left\",\n",
    "\t                  style=\"clip marker paths, width=4.25cm, ymin=-0.6, ymax=0.6, xmin=$a₀, xmax=$b₀, xlabel=\\$x\\$, ylabel=\\$f'\\$, ytick={0}, yticklabels=\\\\empty, xtick=\\\\empty, axis on top=true\")\n",
    "\n",
    "\tfor i in 1 : N\n",
    "\n",
    "\t    x = (a+b)/2\n",
    "\t    y = f(x)\n",
    "\t    if y == 0\n",
    "\t        a, b = x, x\n",
    "\t    elseif sign(y) == sign(ya)\n",
    "\t        a = x\n",
    "\t    else\n",
    "\t        b = x\n",
    "\t    end\n",
    "\n",
    "\t    p = Plots.Plot[]\n",
    "\t    push!(p, Plots.Linear([a₀, b₀], [0,0], style=\"solid, black!70, mark=none\"))\n",
    "\t    push!(p, Plots.Linear(f, (a₀, b₀), style=\"solid, thick, black, mark=none\"))\n",
    "\t    push!(p, Plots.Linear([a,a], [-0.6,0.6], style=\"name path=A, draw=none, mark=none\"))\n",
    "\t    push!(p, Plots.Linear([b,b], [-0.6,0.6], style=\"name path=B, draw=none, mark=none\"))\n",
    "\t    push!(p, Plots.Command(\"\\\\addplot[pastelBlue!50] fill between[of=A and B]\"))\n",
    "\t    push!(p, Plots.Scatter([a,b], [f(a), f(b)], style=\"mark=*, mark size=1, mark options={solid, draw=black, fill=black}\"))\n",
    "\n",
    "\t    push!(g, Axis(p)) # style=\"xtick={$a,$b}, xticklabels={\\$a\\$,\\$b\\$}\"\n",
    "\tend\n",
    "\n",
    "\tg\n",
    "end\n",
    "\n",
    "plot(p)"
   ]
  }
 ],
 "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
}