surrogate-optimization.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Surrogate Optimization\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": {},
   "outputs": [],
   "source": [
    "include(\"support_code.jl\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "include(\"gp.jl\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "p = let\n",
    "\tN = 4\n",
    "\ty_dom = (-2,2)\n",
    "\tG = GroupPlot(N, 1, groupStyle=\"horizontal sep=0.25cm, ylabels at=edge left\",\n",
    "\t                    style=\"width=4.5cm, xlabel=\\$x\\$, ylabel=\\$y\\$, xtick=\\\\empty, ytick=\\\\empty, xmin=0, xmax=8, ymin=$(y_dom[1]), ymax=$(y_dom[2])\")\n",
    "\n",
    "\tf_true = x -> sin(x) + 0.25*sin(x/2) - 0.3*sin(2x)\n",
    "\tx_train = [3.5, 5.5]\n",
    "\tx_arr = unique(sort!(append!(collect(range(0, stop=8, length=201)), x_train)))\n",
    "\tX = [[x] for x in collect(x_arr)]\n",
    "\ty = (x -> f_true(x[1])).(X)\n",
    "\n",
    "\tp_true = Plots.Linear(x_arr, y, style=\"solid, black, mark=none\", legendentry=\"true\")\n",
    "\n",
    "\t# create training set\n",
    "\tGP = GaussianProcess()\n",
    "\tfor x in x_train\n",
    "\t\tpush!(GP, [x], f_true(x))\n",
    "\tend\n",
    "\n",
    "\tfor n in 1 : N\n",
    "\t    p_train = Plots.Scatter([x[1] for x in GP.X], GP.y, style=\"only marks, mark=*, mark size=1, mark options={draw=black, fill=black}\", legendentry=\"fit points\")\n",
    "\n",
    "\t    # predict\n",
    "\t    μₚ, νₚ = predict(GP, X)\n",
    "\t    σₚ = sqrt.(νₚ)\n",
    "\t    p_pred_μ = Plots.Linear(x_arr, μₚ, style=\"solid, thick, pastelBlue, mark=none\", legendentry=\"predicted\")\n",
    "\t    upperConfidence = μₚ + 1.96*σₚ\n",
    "\t    lowerConfidence = μₚ - 1.96*σₚ\n",
    "\n",
    "\t    # prediction-based optimization\n",
    "\t    pred_μ = x -> predict(GP, Vector{Float64}[Float64[x]])[1][1]\n",
    "\t    x_slice = myopt(pred_μ)\n",
    "\n",
    "\t    # slice\n",
    "\t    y_dom = (-2,2)\n",
    "\t    p_sample = Plots.Scatter([x_slice], [f_true(x_slice)], style=\"only marks, mark=*, mark size=1, mark options={draw=pastelRed, fill=pastelRed}\", legendentry=\"next sample\")\n",
    "\t    Plots.Linear([x_slice, x_slice], [y_dom[1], y_dom[2]], style=\"pastelBlue!40, mark=none\")\n",
    "\n",
    "\t    ax = Axis([p_true, p_pred_μ, p_train,  p_sample])\n",
    "\t    if n == N\n",
    "\t        ax.legendPos=\"outer north east\"\n",
    "\t    else\n",
    "\t        push!(ax.plots, Plots.Command(\"\\\\legend{}\"))\n",
    "\t    end\n",
    "\t    push!(G, ax)\n",
    "\n",
    "\t    push!(GP, [x_slice], f_true(x_slice))\n",
    "\t    push!(x_train, x_slice)\n",
    "\t    x_arr = unique(sort!(append!(collect(range(0, stop=8, length=201)), x_train)))\n",
    "\t    X = [[x] for x in collect(x_arr)]\n",
    "\tend\n",
    "\n",
    "\tG\n",
    "end\n",
    "plot(p)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "using Optim\n",
    "p = let\n",
    "\tN = 4\n",
    "\tM = 2\n",
    "\ty_dom = (-2,2)\n",
    "\tG = GroupPlot(N, M, groupStyle=\"horizontal sep=0.25cm, vertical sep=0.25cm, ylabels at=edge left, xlabels at=edge bottom\",\n",
    "\t                    style=\"axis on top=true, clip marker paths=true, width=4.5cm, xlabel=\\$x\\$, ylabel=\\$y\\$, xtick=\\\\empty, ytick=\\\\empty, xmin=0, xmax=8, ymin=$(y_dom[1]), ymax=$(y_dom[2])\")\n",
    "\n",
    "\tf_true = x -> sin(x) + 0.25*sin(x/2) - 0.3*sin(2x)\n",
    "\tx_train = [3.5, 5.5]\n",
    "\tx_arr = unique(sort!(append!(collect(range(0, stop=8, length=101)), x_train)))\n",
    "\tX = [[x] for x in collect(x_arr)]\n",
    "\ty = (x -> f_true(x[1])).(X)\n",
    "\n",
    "\tp_true = Plots.Linear(x_arr, y, style=\"solid, black, mark=none\", legendentry=\"true\")\n",
    "\n",
    "\t# create training set\n",
    "\tGP = GaussianProcess()\n",
    "\tfor x in x_train\n",
    "\t\tpush!(GP, [x], f_true(x))\n",
    "\tend\n",
    "\n",
    "\tfor n in 1 : N*M\n",
    "\t    p_train = Plots.Scatter([x[1] for x in GP.X], GP.y, style=\"only marks, solid, mark=*, mark size=1, mark options={draw=black, fill=black}\", legendentry=\"fit points\")\n",
    "\n",
    "\t    # predict\n",
    "\t    μₚ, νₚ = predict(GP, X)\n",
    "\t    σₚ = sqrt.(νₚ)\n",
    "\t    p_pred_μ = Plots.Linear(x_arr, μₚ, style=\"solid, thick, pastelBlue, mark=none\", legendentry=\"predicted\")\n",
    "\t    upperConfidence = μₚ + 1.96*σₚ\n",
    "\t    lowerConfidence = μₚ - 1.96*σₚ\n",
    "\t    p_pred_hi = Plots.Linear(x_arr, upperConfidence, style=\"draw=none, mark=none, name path=A, forget plot\")\n",
    "\t\tp_pred_lo = Plots.Linear(x_arr, lowerConfidence, style=\"draw=none, mark=none, name path=B, forget plot\")\n",
    "\t\tp_pred = Plots.Command(\"\\\\addplot[pastelBlue!40] fill between[of=A and B];\\\\addlegendentry{confidence region}\")\n",
    "\n",
    "\t    # prediction-based optimization\n",
    "\t    pred_ν = x->predict(GP, Vector{Float64}[Float64[x]])[2][1]\n",
    "\t    x_slice = myopt(x->-pred_ν(x))\n",
    "\n",
    "\t    # slice\n",
    "\t    y_dom = (-2,2)\n",
    "\t    p_sample = Plots.Scatter([x_slice], [f_true(x_slice)], style=\"only marks, mark=*, mark size=1, mark options={draw=pastelRed, fill=pastelRed}\", legendentry=\"sampled\")\n",
    "\t    Plots.Linear([x_slice, x_slice], [y_dom[1], y_dom[2]], style=\"pastelBlue!40, mark=none\")\n",
    "\n",
    "\t    ax = Axis([p_pred_hi, p_pred_lo, p_pred, p_true, p_pred_μ, p_train, p_sample])\n",
    "\t    if n == N\n",
    "\t        ax.legendPos=\"outer north east\"\n",
    "\t    else\n",
    "\t        push!(ax.plots, Plots.Command(\"\\\\legend{}\"))\n",
    "\t    end\n",
    "\t    push!(G, ax)\n",
    "\n",
    "\t    push!(GP, [x_slice], f_true(x_slice))\n",
    "\t    push!(x_train, x_slice)\n",
    "\n",
    "\t    i = findfirst(x_arr .> x_slice)\n",
    "\t    if min(abs(x_arr[i] - x_slice), abs(x_arr[i-1] - x_slice)) > 1e-4\n",
    "\t    \tinsert!(x_arr, i, x_slice)\n",
    "\t    end\n",
    "\t    X = [[x] for x in collect(x_arr)]\n",
    "\tend\n",
    "\n",
    "\tG\n",
    "end\n",
    "plot(p)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "using Optim\n",
    "p = let\n",
    "\tN = 4\n",
    "\tM = 2\n",
    "\ty_dom = (-2,2)\n",
    "\tG = GroupPlot(N, M, groupStyle=\"horizontal sep=0.25cm, vertical sep=0.25cm, ylabels at=edge left, xlabels at=edge bottom\",\n",
    "\t                    style=\"axis on top=true, clip marker paths=true, width=4.5cm, xlabel=\\$x\\$, ylabel=\\$y\\$, xtick=\\\\empty, ytick=\\\\empty, xmin=0, xmax=8, ymin=$(y_dom[1]), ymax=$(y_dom[2])\")\n",
    "\n",
    "\n",
    "\tf_true = x -> sin(x) + 0.25*sin(x/2) - 0.3*sin(2x)\n",
    "\tx_train = [3.5, 5.5]\n",
    "\tx_arr = unique(sort!(append!(collect(range(0, stop=8, length=201)), x_train)))\n",
    "\tX = [[x] for x in collect(x_arr)]\n",
    "\ty = (x -> f_true(x[1])).(X)\n",
    "\n",
    "\tp_true = Plots.Linear(x_arr, y, style=\"solid, black, mark=none\", legendentry=\"true\")\n",
    "\n",
    "\t# create training set\n",
    "\tGP = GaussianProcess()\n",
    "\tfor x in x_train\n",
    "\t\tpush!(GP, [x], f_true(x))\n",
    "\tend\n",
    "\n",
    "\tfor n in 1 : N*M\n",
    "\t    p_train = Plots.Scatter([x[1] for x in GP.X], GP.y, style=\"only marks, mark=*, mark size=1, mark options={draw=black, fill=black}\", legendentry=\"fit points\")\n",
    "\n",
    "\t    # predict\n",
    "\t    μₚ, νₚ = predict(GP, X)\n",
    "\t    σₚ = sqrt.(νₚ)\n",
    "\t    p_pred_μ = Plots.Linear(x_arr, μₚ, style=\"solid, thick, pastelBlue, mark=none\", legendentry=\"predicted\")\n",
    "\t    upperConfidence = μₚ + 1.96*σₚ\n",
    "\t    lowerConfidence = μₚ - 1.96*σₚ\n",
    "\t    p_pred_hi = Plots.Linear(x_arr, upperConfidence, style=\"draw=none, mark=none, name path=A, forget plot\")\n",
    "\t\tp_pred_lo = Plots.Linear(x_arr, lowerConfidence, style=\"draw=none, mark=none, name path=B, forget plot\")\n",
    "\t\tp_pred = Plots.Command(\"\\\\addplot[pastelBlue!40] fill between[of=A and B];\\\\addlegendentry{confidence region}\")\n",
    "\n",
    "\t    # prediction-based optimization\n",
    "\t    pred_μ = x->predict(GP, Vector{Float64}[Float64[x]])[1][1]\n",
    "\t    pred_ν = x->predict(GP, Vector{Float64}[Float64[x]])[2][1]\n",
    "\t    f_opt = x->pred_μ(x) - sqrt(pred_ν(x))\n",
    "\t    x_slice = myopt(f_opt)\n",
    "\t    p_opt = Plots.Linear(x_arr, f_opt.(x_arr), style=\"solid, thick, pastelRed, mark=none\", legendentry=\"lower bound\")\n",
    "\t    p_opt_vert = Plots.Linear([x_slice, x_slice], [f_opt(x_slice), f_true(x_slice)], style=\"solid, pastelRed!40, mark=none, forget plot\")\n",
    "\n",
    "\t    # slice\n",
    "\t    y_dom = (-2,2)\n",
    "\t    p_sample = Plots.Scatter([x_slice], [f_true(x_slice)], style=\"only marks, mark=*, mark size=1, mark options={draw=pastelRed, fill=pastelRed}\", legendentry=\"sampled\")\n",
    "\t    Plots.Linear([x_slice, x_slice], [y_dom[1], y_dom[2]], style=\"pastelBlue!40, mark=none\")\n",
    "\n",
    "\t    ax = Axis([p_pred_hi, p_pred_lo, p_pred, p_true, p_opt_vert, p_opt, p_pred_μ, p_train, p_sample])\n",
    "\t    if n == N\n",
    "\t        ax.legendPos=\"outer north east\"\n",
    "\t    else\n",
    "\t        push!(ax.plots, Plots.Command(\"\\\\legend{}\"))\n",
    "\t    end\n",
    "\t    push!(G, ax)\n",
    "\n",
    "\t    push!(GP, [x_slice], f_true(x_slice))\n",
    "\t    push!(x_train, x_slice)\n",
    "\t    x_arr = unique(sort!(append!(collect(range(0, stop=8, length=201)), x_train)))\n",
    "\t    X = [[x] for x in collect(x_arr)]\n",
    "\tend\n",
    "\n",
    "\tG\n",
    "end\n",
    "plot(p)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "using Random\n",
    "using Vec\n",
    "\n",
    "p = let\n",
    "\n",
    "\tGP = GaussianProcess()\n",
    "\n",
    "\t# create random \"true\" function\n",
    "\tRandom.seed!(3)\n",
    "\tx_arr = range(0, stop=8, length=101)\n",
    "\tX = [[x] for x in collect(x_arr)]\n",
    "\ty = rand(GP, X)\n",
    "\n",
    "\tp_true = Plots.Linear(x_arr, y, style=\"solid, black, mark=none\", legendentry=\"true\")\n",
    "\n",
    "\t# create training set\n",
    "\tfor i in [5, 20, 35, 50]\n",
    "\t\tpush!(GP, X[i], y[i])\n",
    "\tend\n",
    "\tp_train = Plots.Scatter([x[1] for x in GP.X], GP.y, style=\"only marks, mark=*, mark size=1, mark options={draw=black, fill=black}\", legendentry=\"fit points\")\n",
    "\n",
    "\t# predict\n",
    "\tμₚ, νₚ = predict(GP, X)\n",
    "\tσₚ = sqrt.(νₚ)\n",
    "\tp_pred_μ = Plots.Linear(x_arr, μₚ, style=\"solid, thick, pastelBlue, mark=none\", legendentry=\"predicted\")\n",
    "\tupperConfidence = μₚ + 1.96*σₚ\n",
    "\tlowerConfidence = μₚ - 1.96*σₚ\n",
    "\n",
    "\t# slice\n",
    "\tslice_index = 70\n",
    "\ty_dom = (-3,3)\n",
    "\tx_slice = X[slice_index][1]\n",
    "\tμ_slice = μₚ[slice_index]\n",
    "\tσ_slice = σₚ[slice_index]\n",
    "\tN = Normal(μ_slice, σ_slice)\n",
    "\ty_arr = collect(range(y_dom[1], stop=y_dom[2], length=101))\n",
    "\tx_arr = [x_slice + pdf(N, y) for y in y_arr]\n",
    "\tp_slice = Plots.Linear(x_arr, y_arr, style=\"pastelBlue!40, mark=none\")\n",
    "\tp_slice_vert = Plots.Linear([x_slice, x_slice], [y_arr[1], y_arr[end]], style=\"black, mark=none\")\n",
    "\n",
    "\ty_min = minimum(GP.y)\n",
    "\tx_min = GP.X[argmin(GP.y)][1]\n",
    "\tp_fill_lo = Plots.Linear([x_slice, x_slice], [y_arr[1], y_min], style=\"draw=none, mark=none, name path=C, forget plot\")\n",
    "\ty_arr_hi = collect(range(y_dom[1], stop=y_min, length=51))\n",
    "\tx_arr_hi = [x_slice + pdf(N, y) for y in y_arr_hi]\n",
    "\tp_fill_hi = Plots.Linear(x_arr_hi, y_arr_hi, style=\"draw=none, mark=none, name path=D, forget plot\")\n",
    "\tp_fill = Plots.Command(\"\\\\addplot[pastelBlue!40, forget plot] fill between[of=C and D];\")\n",
    "\tp_flat = Plots.Linear([x_min, x_slice], [y_min, y_min], style=\"solid, pastelBlue!40, mark=none, forget plot\")\n",
    "\n",
    "\tarow_best_so_far = Plots.Command(\"\\\\draw [->, thick, shorten >=5pt] (axis description cs:0.3,0.75) -- (axis cs:$(x_min),$(y_min));\")\n",
    "\tnode_best_so_far = Plots.Node(\"\\\\small best so far\", 0.3, 0.75, style=\"anchor=south east\", axis=\"axis description cs\")\n",
    "\n",
    "\tA = VecE2(x_slice-1, y_min-1)\n",
    "\tarow_prob_of_imp = Plots.Command(\"\\\\draw [->, thick] (axis cs:$(A.x),$(A.y)) -- (axis cs:$(x_slice+0.08),$(A.y+0.3));\")\n",
    "\tnode_prob_of_imp = Plots.Node(\"\\\\small probability of improvement\", A.x, A.y, style=\"anchor=east\", axis=\"axis cs\")\n",
    "\n",
    "\tB = A + VecE2(0, -0.75)\n",
    "\tarow_query_point = Plots.Command(\"\\\\draw [->, thick, shorten >=5pt] (axis cs:$(B.x),$(B.y)) -- (axis cs:$(x_slice),$(B.y));\")\n",
    "\tnode_query_point = Plots.Node(\"\\\\small query point\", B.x, B.y, style=\"anchor=east\", axis=\"axis cs\")\n",
    "\n",
    "\tAxis([p_true, p_pred_μ, p_train, p_slice, p_slice_vert, p_flat, p_fill_lo, p_fill_hi, p_fill,\n",
    "\t    node_best_so_far, arow_best_so_far, node_prob_of_imp, arow_prob_of_imp, node_query_point, arow_query_point],\n",
    "\t    xmin=0, xmax=8, style=\"ytick=\\\\empty, xtick=\\\\empty\",\n",
    "\t    ymin=y_dom[1], ymax=y_dom[2],\n",
    "\t    xlabel=L\"x\", ylabel=L\"y\", width=\"9cm\", legendPos=\"outer north east\",\n",
    "\t)\n",
    "end\n",
    "\n",
    "plot(p)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "prob_of_improvement(y_min, μ, ν) = isapprox(ν, 0, atol=1e-4) ? 0.0 : cdf(Normal(μ, sqrt(ν)), y_min)\n",
    "\n",
    "p = let\n",
    "\n",
    "\tN = 4\n",
    "\ty_dom = (-2,2)\n",
    "\tG = GroupPlot(N, 2, groupStyle=\"horizontal sep=0.25cm, ylabels at=edge left, vertical sep=0.25cm, xlabels at=edge bottom\",\n",
    "\t                    style=\"axis on top=true, clip marker paths=true, width=4.5cm, xlabel=\\$x\\$, xtick=\\\\empty, ytick=\\\\empty, xmin=0, xmax=8, ymin=$(y_dom[1]), ymax=$(y_dom[2])\")\n",
    "\n",
    "\tf_true = x -> sin(x) + 0.25*sin(x/2) - 0.3*sin(2x)\n",
    "\tx_train = [3.5, 5.5]\n",
    "\tx_arr = unique(sort!(append!(collect(range(0, stop=8, length=201)), x_train)))\n",
    "\tX = [[x] for x in collect(x_arr)]\n",
    "\ty = (x -> f_true(x[1])).(X)\n",
    "\n",
    "\tp_true = Plots.Linear(x_arr, y, style=\"solid, black, mark=none\", legendentry=\"true\")\n",
    "\n",
    "\t# create training set\n",
    "\tGP = GaussianProcess()\n",
    "\tfor x in x_train\n",
    "\t\tpush!(GP, [x], f_true(x[1]))\n",
    "\tend\n",
    "\n",
    "\taxes_lo = Axis[]\n",
    "\n",
    "\tfor n in 1 : N\n",
    "\t    p_train = Plots.Scatter([x[1] for x in GP.X], GP.y, style=\"only marks, mark=*, mark size=1, mark options={draw=black, fill=black}\", legendentry=\"fit points\")\n",
    "\n",
    "\t    # predict\n",
    "\t    μₚ, νₚ = predict(GP, X)\n",
    "\t    σₚ = sqrt.(νₚ)\n",
    "\t    p_pred_μ = Plots.Linear(x_arr, μₚ, style=\"solid, thick, pastelBlue, mark=none\", legendentry=\"predicted\")\n",
    "\t    upperConfidence = μₚ + 1.96*σₚ\n",
    "\t    lowerConfidence = μₚ - 1.96*σₚ\n",
    "\t    p_pred_hi = Plots.Linear(x_arr, upperConfidence, style=\"draw=none, mark=none, name path=A, forget plot\")\n",
    "\t\tp_pred_lo = Plots.Linear(x_arr, lowerConfidence, style=\"draw=none, mark=none, name path=B, forget plot\")\n",
    "\t\tp_pred = Plots.Command(\"\\\\addplot[pastelBlue!40] fill between[of=A and B];\\\\addlegendentry{confidence region}\")\n",
    "\n",
    "\t    # prediction-based optimization\n",
    "\t    pred_μ = x->predict(GP, Vector{Float64}[Float64[x]])[1][1]\n",
    "\t    pred_ν = x->predict(GP, Vector{Float64}[Float64[x]])[2][1]\n",
    "\t    y_min = minimum(GP.y)\n",
    "\t    f_opt = x->prob_of_improvement(y_min, pred_μ(x), pred_ν(x))\n",
    "\t    x_slice = myopt(x->-f_opt(x))\n",
    "\t    p_opt = Plots.Linear(x_arr, f_opt.(x_arr), style=\"solid, pastelRed, mark=none\")\n",
    "\t    p_opt_vert = Plots.Linear([x_slice, x_slice], [f_opt(x_slice), f_true(x_slice)], style=\"solid, pastelRed!40, mark=none\")\n",
    "\n",
    "\t    # slice\n",
    "\t    y_dom = (-2,2)\n",
    "\t    p_sample = Plots.Scatter([x_slice], [f_true(x_slice)], style=\"only marks, mark=*, mark size=1, mark options={draw=pastelRed, fill=pastelRed}\", legendentry=\"sampled\")\n",
    "\t    Plots.Linear([x_slice, x_slice], [y_dom[1], y_dom[2]], style=\"pastelBlue!40, mark=none\")\n",
    "\n",
    "\t    push!(axes_lo, Axis([p_opt,\n",
    "\t                Plots.Scatter([x_slice], [f_opt(x_slice)], style=\"axis on top, only marks, mark=*, mark size=1, mark options={draw=pastelRed, fill=pastelRed}\")\n",
    "\t                ], ymin=0, ymax=1))\n",
    "\n",
    "\t    ax = Axis([p_pred_hi, p_pred_lo, p_pred, p_true, p_pred_μ, p_train, p_sample])\n",
    "\t    if n == N\n",
    "\t        ax.legendPos=\"outer north east\"\n",
    "\t    else\n",
    "\t        push!(ax.plots, Plots.Command(\"\\\\legend{}\"))\n",
    "\t    end\n",
    "\t    if n == 1\n",
    "\t        ax.ylabel=L\"y\"\n",
    "\t        axes_lo[1].ylabel=L\"P[f(\\vect x) < y_\\text{min}]\"\n",
    "\t        axes_lo[1].style = \"ytick={0,1}\"\n",
    "\t    end\n",
    "\t    push!(G, ax)\n",
    "\t    push!(GP, [x_slice], f_true(x_slice))\n",
    "\t    push!(x_train, x_slice)\n",
    "\t    x_arr = unique(sort!(append!(collect(range(0, stop=8, length=201)), x_train)))\n",
    "\t    X = [[x] for x in collect(x_arr)]\n",
    "\tend\n",
    "\tfor ax in axes_lo\n",
    "\t    push!(G, ax)\n",
    "\tend\n",
    "\n",
    "\tG\n",
    "end\n",
    "plot(p)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "function expected_improvement(y_min, μ, ν)\n",
    "\tif isapprox(ν, 0.0, atol=1e-4)\n",
    "\t\treturn 0.0\n",
    "\tend\n",
    "\tσ = sqrt(ν)\n",
    "    tmp = (y_min - μ)/σ\n",
    "    p_imp = prob_of_improvement(y_min, μ, ν)\n",
    "    p_ymin = pdf(Normal(μ, σ), y_min)\n",
    "    return (y_min - μ)*p_imp + σ^2*p_ymin\n",
    "end\n",
    "\n",
    "p = let\n",
    "\n",
    "\tN = 4\n",
    "\ty_dom = (-2,2)\n",
    "\tG = GroupPlot(N, 2, groupStyle=\"horizontal sep=0.25cm, ylabels at=edge left, vertical sep=0.25cm, xlabels at=edge bottom\",\n",
    "\t                    style=\"axis on top, clip marker paths=true, width=4.5cm, xlabel=\\$x\\$, xtick=\\\\empty, ytick=\\\\empty, xmin=0, xmax=8, ymin=$(y_dom[1]), ymax=$(y_dom[2])\")\n",
    "\n",
    "\tf_true = x -> sin(x) + 0.25*sin(x/2) - 0.3*sin(2x)\n",
    "\tx_train = [3.5, 5.5]\n",
    "\tx_arr = unique(sort!(append!(collect(range(0, stop=8, length=201)), x_train)))\n",
    "\tX = [[x] for x in collect(x_arr)]\n",
    "\ty = (x -> f_true(x[1])).(X)\n",
    "\n",
    "\tp_true = Plots.Linear(x_arr, y, style=\"solid, black, mark=none\", legendentry=\"true\")\n",
    "\n",
    "\t# create training set\n",
    "\tGP = GaussianProcess()\n",
    "\tfor x in x_train\n",
    "\t\tpush!(GP, [x], f_true(x))\n",
    "\tend\n",
    "\n",
    "\taxes_lo = Axis[]\n",
    "\n",
    "\tfor n in 1 : N\n",
    "\t    p_train = Plots.Scatter([x[1] for x in GP.X], GP.y, style=\"only marks, mark=*, mark size=1, mark options={draw=black, fill=black}\", legendentry=\"fit points\")\n",
    "\n",
    "\t    # predict\n",
    "\t    μₚ, νₚ = predict(GP, X)\n",
    "\t    σₚ = sqrt.(νₚ)\n",
    "\t    p_pred_μ = Plots.Linear(x_arr, μₚ, style=\"solid, thick, pastelBlue, mark=none\", legendentry=\"predicted\")\n",
    "\t    upperConfidence = μₚ + 1.96*σₚ\n",
    "\t    lowerConfidence = μₚ - 1.96*σₚ\n",
    "\t    p_pred_hi = Plots.Linear(x_arr, upperConfidence, style=\"draw=none, mark=none, name path=A, forget plot\")\n",
    "\t\tp_pred_lo = Plots.Linear(x_arr, lowerConfidence, style=\"draw=none, mark=none, name path=B, forget plot\")\n",
    "\t\tp_pred = Plots.Command(\"\\\\addplot[pastelBlue!40] fill between[of=A and B];\\\\addlegendentry{confidence region}\")\n",
    "\n",
    "\t    # prediction-based optimization\n",
    "\t    pred_μ = x->predict(GP, Vector{Float64}[Float64[x]])[1][1]\n",
    "\t    pred_ν = x->predict(GP, Vector{Float64}[Float64[x]])[2][1]\n",
    "\t    y_min = minimum(GP.y)\n",
    "\t    f_opt = x->expected_improvement(y_min, pred_μ(x), pred_ν(x) + 1e-8)\n",
    "\t    x_slice = myopt(x->-f_opt(x))\n",
    "\t    p_opt = Plots.Linear(x_arr, f_opt.(x_arr), style=\"solid, pastelRed, mark=none\")\n",
    "\t    p_opt_vert = Plots.Linear([x_slice, x_slice], [f_opt(x_slice), f_true(x_slice)], style=\"solid, pastelRed!40, mark=none\")\n",
    "\n",
    "\t    # slice\n",
    "\t    y_dom = (-2,2)\n",
    "\t    p_sample = Plots.Scatter([x_slice], [f_true(x_slice)], style=\"only marks, mark=*, mark size=1, mark options={draw=pastelRed, fill=pastelRed}\", legendentry=\"sampled\")\n",
    "\t    Plots.Linear([x_slice, x_slice], [y_dom[1], y_dom[2]], style=\"pastelBlue!40, mark=none\")\n",
    "\n",
    "\t    push!(axes_lo, Axis([p_opt,\n",
    "\t                Plots.Scatter([x_slice], [f_opt(x_slice)], style=\"axis on top, only marks, mark=*, mark size=1, mark options={draw=pastelRed, fill=pastelRed}\")\n",
    "\t                ], ymin=0, ymax=1))\n",
    "\n",
    "\t    ax = Axis([p_pred_hi, p_pred_lo, p_pred, p_true, p_pred_μ, p_train, p_sample])\n",
    "\t    if n == N\n",
    "\t        ax.legendPos=\"outer north east\"\n",
    "\t    else\n",
    "\t        push!(ax.plots, Plots.Command(\"\\\\legend{}\"))\n",
    "\t    end\n",
    "\t    if n == 1\n",
    "\t        ax.ylabel=L\"y\"\n",
    "\t        axes_lo[1].ylabel = L\"E[I(y)]\"\n",
    "\t        axes_lo[1].style = \"ytick={0,1}\"\n",
    "\t    end\n",
    "\t    push!(G, ax)\n",
    "\t    push!(GP, [x_slice], f_true(x_slice))\n",
    "\t    push!(x_train, x_slice)\n",
    "\t    x_arr = unique(sort!(append!(collect(range(0, stop=8, length=201)), x_train)))\n",
    "\t    X = [[x] for x in collect(x_arr)]\n",
    "\tend\n",
    "\tfor ax in axes_lo\n",
    "\t    push!(G, ax)\n",
    "\tend\n",
    "\n",
    "\tG\n",
    "end\n",
    "plot(p)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "if !@isdefined(GaussianProcess)\n",
    "\tinclude(\"../chapter/gp.jl\")\n",
    "end\n",
    "\n",
    "using Distributions\n",
    "p = let\n",
    "\tZ = Normal(0.0, 1.0)\n",
    "\tdom=(-2,6)\n",
    "\n",
    "\tfunction f_(x::Float64, z::Float64=0.0)\n",
    "\t    A = Normal(1.0, 1.0)\n",
    "\t    B = Normal(3.0, 0.5)\n",
    "\t    C = Normal(5.0, 0.25)\n",
    "\t    return 1 - (2pdf(A, x+z) + pdf(B, x+z) + pdf(C, x+z))\n",
    "\tend\n",
    "\n",
    "\ty_max = 0.8\n",
    "\tlo, hi = -0.8, 1.0\n",
    "\n",
    "\tp = Plots.Plot[]\n",
    "\n",
    "\t# append!(p, plot_transparent_interval((-0.6635, 3.8425), (lo,hi), style=\"green!50\"))\n",
    "\t# append!(p, plot_transparent_interval((4.4838, 5.5095), (lo,hi), style=\"green!50\"))\n",
    "\tappend!(p, plot_transparent_intervals([(-0.6635, 3.8425), (4.4838, 5.5095)], (lo,hi), \"green\", 0.5))\n",
    "\n",
    "\tpush!(p, Plots.Linear(f_, dom, xbins=201, style=\"solid, black, mark=none\"))\n",
    "\tpush!(p, Plots.Linear([dom[1], dom[2]], [y_max, y_max], style=\"solid, gray, mark=none\"))\n",
    "\tpush!(p, Plots.Scatter([5.0], [f_(5.0)], style=\"only marks, black, mark=*, mark size=1, mark options={draw=black, fill=black}\"))\n",
    "\tpush!(p, Plots.Node(L\"x^*\", 5.0, f_(5.0), style=\"right\"))\n",
    "\n",
    "\tAxis(p, width=\"9cm\", xlabel=L\"x\", ylabel=L\"y\", ymin=lo, ymax=hi, style=\"enlarge x limits=0, xtick=\\\\empty, ytick={$(y_max)}, yticklabel={\\$y_\\\\text{max}\\$}, axis on top\",\n",
    "\t)\n",
    "\n",
    "end\n",
    "plot(p)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "using Random\n",
    "\n",
    "p = let\n",
    "\tP_safe = 0.9\n",
    "\tσ_safe = optimize(x->(cdf(Normal(0.0,1.0), x) - P_safe)^2, -5, 5).minimizer\n",
    "\tβ = σ_safe*σ_safe\n",
    "\ty_max = 0.8\n",
    "\n",
    "\tfunction f_(x::Float64, z::Float64=0.0)\n",
    "\t    A = Normal(1.0, 1.0)\n",
    "\t    B = Normal(3.0, 0.5)\n",
    "\t    C = Normal(5.0, 0.25)\n",
    "\t    return 1 - (2pdf(A, x+z) + pdf(B, x+z) + pdf(C, x+z))\n",
    "\tend\n",
    "\tf_true = x -> f_(x,0.0)\n",
    "\tdom = (-2,6)\n",
    "\tGP = GaussianProcess(ν = 0.01)\n",
    "\n",
    "\tx_train = [0.5,1.0,2.0,3.0,3.5]\n",
    "\tx_arr = unique(sort!(append!(collect(range(dom[1], stop=dom[2], length=201)), x_train)))\n",
    "\tX = [[x] for x in collect(x_arr)]\n",
    "\ty = (x -> f_true(x[1])).(X)\n",
    "\tp_true = Plots.Linear(x_arr, y, style=\"solid, black, mark=none\", legendentry=\"true\")\n",
    "\n",
    "\t# create training set\n",
    "\tRandom.seed!(0)\n",
    "\tfor x in x_train\n",
    "\t    push!(GP, [x], f_true(x) + randn()*sqrt(GP.ν))\n",
    "\tend\n",
    "\tp_train = plot_GP_data(GP, legendentry=\"fit points\")\n",
    "\n",
    "\t# predict\n",
    "\t(μₚ, νₚ) = predict(GP, X)\n",
    "\tσₚ = sqrt.(νₚ)\n",
    "\tp_pred_μ = Plots.Linear(x_arr, μₚ, style=\"solid, thick, pastelBlue, mark=none\", legendentry=\"predicted\")\n",
    "\tupperConfidence = μₚ + σ_safe*σₚ\n",
    "\tlowerConfidence = μₚ - σ_safe*σₚ\n",
    "\tp_pred_hi = Plots.Linear(x_arr, upperConfidence, style=\"draw=none, mark=none, name path=A, forget plot\")\n",
    "\tp_pred_lo = Plots.Linear(x_arr, lowerConfidence, style=\"draw=none, mark=none, name path=B, forget plot\")\n",
    "\tp_pred = Plots.Command(\"\\\\addplot[pastelBlue!40] fill between[of=A and B];\\\\addlegendentry{confidence region}\")\n",
    "\tp_lowerbound = Plots.Linear(x_arr, lowerConfidence, style=\"solid, thick, pastelGreen, mark=none\", legendentry=L\"\\ell(x)\")\n",
    "\tp_upperbound = Plots.Linear(x_arr, upperConfidence, style=\"solid, thick, pastelRed, mark=none\", legendentry=L\"u(x)\")\n",
    "\n",
    "\tAxis([p_pred_hi, p_pred_lo, p_pred, p_true, p_pred_μ, p_lowerbound, p_upperbound, p_train],\n",
    "\t\t\t\tstyle=\"legend pos=outer north east, axis on top=true, clip marker paths=true, width=9cm, height=4cm, xlabel=\\$x\\$, xtick=\\\\empty, ytick=\\\\empty, enlarge x limits=0, legend style={draw=none}\")\n",
    "end\n",
    "plot(p)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "using Random\n",
    "\n",
    "p = let\n",
    "\tlo, hi = -1.0, 1.2\n",
    "\n",
    "\tP_safe = 0.9\n",
    "\tσ_safe = optimize(x->(cdf(Normal(0.0,1.0), x) - P_safe)^2, -5, 5).minimizer\n",
    "\tβ = σ_safe*σ_safe\n",
    "\ty_max = 0.8\n",
    "\n",
    "\tfunction f_(x::Float64, z::Float64=0.0)\n",
    "\t    A = Normal(1.0, 1.0)\n",
    "\t    B = Normal(3.0, 0.5)\n",
    "\t    C = Normal(5.0, 0.25)\n",
    "\t    return 1 - (2pdf(A, x+z) + pdf(B, x+z) + pdf(C, x+z))\n",
    "\tend\n",
    "\tf_true = x -> f_(x,0.0)\n",
    "\tdom = (-2,6)\n",
    "\tGP = GaussianProcess(ν = 0.01)\n",
    "\n",
    "\tx_train = [0.5,1.0,2.0,3.0,3.5]\n",
    "\tx_arr = unique(sort!(append!(collect(range(dom[1], stop=dom[2], length=201)), x_train)))\n",
    "\tX = [[x] for x in collect(x_arr)]\n",
    "\ty = (x -> f_true(x[1])).(X)\n",
    "\n",
    "\tp_true = Plots.Linear(x_arr, y, style=\"solid, black, mark=none\", legendentry=L\"\\text{objective function } f\")\n",
    "\n",
    "\t# create training set\n",
    "\tRandom.seed!(0)\n",
    "\tfor x in x_train\n",
    "\t    push!(GP, [x], f_true(x) + randn()*sqrt(GP.ν))\n",
    "\tend\n",
    "\tp_train = plot_GP_data(GP, legendentry=\"fit points\")\n",
    "\n",
    "\t# predict\n",
    "\t(μₚ, νₚ) = predict(GP, X)\n",
    "\tσₚ = sqrt.(νₚ)\n",
    "\tp_pred_μ = Plots.Linear(x_arr, μₚ, style=\"solid, thick, pastelBlue, mark=none\", legendentry=\"predicted\")\n",
    "\tupperConfidence = μₚ + σ_safe*σₚ\n",
    "\tlowerConfidence = μₚ - σ_safe*σₚ\n",
    "\tp_pred_hi = Plots.Linear(x_arr, upperConfidence, style=\"draw=none, mark=none, name path=A, forget plot\")\n",
    "\tp_pred_lo = Plots.Linear(x_arr, lowerConfidence, style=\"draw=none, mark=none, name path=B, forget plot\")\n",
    "\tp_pred = Plots.Command(\"\\\\addplot[pastelBlue!40] fill between[of=A and B];\\\\addlegendentry{confidence interval}\")\n",
    "\tp_lowerbound = Plots.Linear(x_arr, lowerConfidence, style=\"solid, thick, pastelGreen, mark=none\", legendentry=L\"\\ell(x)\")\n",
    "\tp_upperbound = Plots.Linear(x_arr, upperConfidence, style=\"solid, thick, pastelRed, mark=none\", legendentry=L\"u(x)\")\n",
    "\n",
    "\tp = Plots.Plot[]\n",
    "\tpush!(p, p_pred_hi)\n",
    "\tpush!(p, p_pred_lo)\n",
    "\tpush!(p, p_pred)\n",
    "\tpush!(p, p_true)\n",
    "\tpush!(p, p_train)\n",
    "\tpush!(p, Plots.Linear([dom[1], dom[2]], [y_max, y_max], style=\"solid, thick, gray, mark=none\", legendentry=\"safety threshold\"))\n",
    "\n",
    "\tsafe_regions = get_safe_regions(GP, β, dom..., y_max)\n",
    "\tappend!(p, plot_transparent_intervals(safe_regions, (lo,hi), \"green\", 0.5, \"estimated safe region \\$\\\\mathcal{S}\\$\"))\n",
    "\n",
    "\tG = GroupPlot(1, 2, groupStyle=\"vertical sep=0.5cm, xlabels at=edge bottom\",\n",
    "\t    style=\"axis on top=true, width=9cm, height=4cm, xlabel=\\$x\\$, xtick=\\\\empty, enlarge x limits=0\")\n",
    "\tpush!(G, Axis(p, xlabel=L\"x\", ylabel=L\"y\", ymin=lo, ymax=hi, style=\"legend pos=outer north east, ytick=\\\\empty, legend style={draw=none}\"))\n",
    "\n",
    "\tpush!(G, Axis(\n",
    "\t    [\n",
    "\t    Plots.Linear(x->prob_is_safe(predict(GP, [x]), y_max), dom, xbins=300, style=\"solid, green, thick, mark=none\"),\n",
    "\t    Plots.Linear([dom[1], dom[2]], [P_safe, P_safe], style=\"solid, gray, mark=none\"),\n",
    "\t    Plots.Node(L\"P_\\text{safe}\", 3, P_safe, style=\"below\")\n",
    "\t    ], ymin=0, ymax=1, ylabel=\"safety probability\", style=\"enlarge x limits=0, ytick={0,1} \"\n",
    "\t))\n",
    "end\n",
    "plot(p)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "using Random\n",
    "\n",
    "p = let\n",
    "\tlo, hi = -1.0, 1.2\n",
    "\n",
    "\tP_safe = 0.9\n",
    "\tσ_safe = optimize(x->(cdf(Normal(0.0,1.0), x) - P_safe)^2, -5, 5).minimizer\n",
    "\tβ = σ_safe*σ_safe\n",
    "\ty_max = 0.8\n",
    "\n",
    "\tfunction f_(x::Float64, z::Float64=0.0)\n",
    "\t    A = Normal(1.0, 1.0)\n",
    "\t    B = Normal(3.0, 0.5)\n",
    "\t    C = Normal(5.0, 0.25)\n",
    "\t    return 1 - (2pdf(A, x+z) + pdf(B, x+z) + pdf(C, x+z))\n",
    "\tend\n",
    "\tf_true = x -> f_(x,0.0)\n",
    "\tdom = (-2,6)\n",
    "\tGP = GaussianProcess(ν = 0.01)\n",
    "\n",
    "\tx_train = [0.5,1.0,2.0,3.0,3.5]\n",
    "\tx_arr = unique(sort!(append!(collect(range(dom[1], stop=dom[2], length=201)), x_train)))\n",
    "\tX = [[x] for x in collect(x_arr)]\n",
    "\ty = (x -> f_true(x[1])).(X)\n",
    "\n",
    "\tp_true = Plots.Linear(x_arr, y, style=\"solid, black, mark=none\", legendentry=\"objective function\")\n",
    "\n",
    "\t# create training set\n",
    "\tRandom.seed!(0)\n",
    "\tfor x in x_train\n",
    "\t    push!(GP, [x], f_true(x) + randn()*sqrt(GP.ν))\n",
    "\tend\n",
    "\tp_train = plot_GP_data(GP, legendentry=\"fit points\")\n",
    "\n",
    "\t# predict\n",
    "\t(μₚ, νₚ) = predict(GP, X)\n",
    "\tσₚ = sqrt.(νₚ)\n",
    "\tp_pred_μ = Plots.Linear(x_arr, μₚ, style=\"solid, thick, pastelBlue, mark=none\", legendentry=\"predicted\")\n",
    "\tupperConfidence = μₚ + σ_safe*σₚ\n",
    "\tlowerConfidence = μₚ - σ_safe*σₚ\n",
    "\tp_pred_hi = Plots.Linear(x_arr, upperConfidence, style=\"draw=none, mark=none, name path=A, forget plot\")\n",
    "\tp_pred_lo = Plots.Linear(x_arr, lowerConfidence, style=\"draw=none, mark=none, name path=B, forget plot\")\n",
    "\tp_pred = Plots.Command(\"\\\\addplot[pastelBlue!40] fill between[of=A and B];\\\\addlegendentry{confidence region}\")\n",
    "\n",
    "\tp = Plots.Plot[]\n",
    "\tpush!(p, p_pred_hi)\n",
    "\tpush!(p, p_pred_lo)\n",
    "\tpush!(p, p_pred)\n",
    "\tpush!(p, p_true)\n",
    "\tpush!(p, p_train)\n",
    "\tpush!(p, Plots.Linear([dom[1], dom[2]], [y_max, y_max], style=\"solid, thick, gray, mark=none\", legendentry=\"safety threshold\"))\n",
    "\n",
    "\t(best_hi, i) = findmin(upperConfidence)\n",
    "\tpush!(p, Plots.Linear([dom[1], dom[2]], [best_hi, best_hi], style=\"solid, thick, pastelRed, mark=none\", legendentry=\"best upper bound\"))\n",
    "\n",
    "\tax3 = Axis(p, xlabel=L\"x\", ylabel=L\"y\", ymin=lo, ymax=hi, style=\"legend pos=outer north east, ytick=\\\\empty, legend style={draw=none}, axis on top=true, width=9.5cm, height=4cm, xlabel=\\$x\\$, xtick=\\\\empty, enlarge x limits=0\")\n",
    "\n",
    "\tsafe_regions = get_safe_regions(GP, β, dom..., y_max)\n",
    "\tax1 = Axis(plot_transparent_intervals(safe_regions, (lo,hi), \"green\", 0.5, \"estimated safe region \\$\\\\mathcal{S}\\$\"),\n",
    "\t        xmin=dom[1], xmax=dom[2], ymin=lo, ymax=hi, style=\"ytick=\\\\empty, xtick=\\\\empty, legend pos=outer north east, axis lines=none, legend style={draw=none}, width=9.5cm, height=2.05cm\")\n",
    "\n",
    "\tM_regions = get_potential_maximizers(GP, β, safe_regions)\n",
    "\tax2 = Axis(plot_transparent_intervals(M_regions, (lo,hi), \"magenta\", 0.5, \"potential minimizers \\$\\\\mathcal{M}\\$\"),\n",
    "\t        xmin=dom[1], xmax=dom[2], ymin=lo, ymax=hi, style=\"ytick=\\\\empty, xtick=\\\\empty, legend pos=outer north east, axis lines=none, legend style={draw=none}, width=9.5cm, height=2.05cm\")\n",
    "\n",
    "\tG = GroupPlot(1, 3, groupStyle=\"vertical sep=0cm\")\n",
    "\tpush!(G, ax1)\n",
    "\tpush!(G, ax2)\n",
    "\tpush!(G, ax3)\n",
    "\tG\n",
    "end\n",
    "plot(p)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "using Random\n",
    "\n",
    "p = let\n",
    "    lo, hi = -1.0, 1.2\n",
    "\n",
    "    P_safe = 0.9\n",
    "    σ_safe = optimize(x->(cdf(Normal(0.0,1.0), x) - P_safe)^2, -5, 5).minimizer\n",
    "    β = σ_safe*σ_safe\n",
    "    y_max = 0.8\n",
    "\n",
    "    function f_(x::Float64, z::Float64=0.0)\n",
    "        A = Normal(1.0, 1.0)\n",
    "        B = Normal(3.0, 0.5)\n",
    "        C = Normal(5.0, 0.25)\n",
    "        return 1 - (2pdf(A, x+z) + pdf(B, x+z) + pdf(C, x+z))\n",
    "    end\n",
    "    f_true = x -> f_(x,0.0)\n",
    "    dom = (-2,6)\n",
    "    GP = GaussianProcess(ν = 0.01)\n",
    "\n",
    "    lo, hi = -1.0, 1.2\n",
    "\n",
    "    x_train = [0.5,1.0,2.0,3.0,3.5]\n",
    "    x_arr = unique(sort!(append!(collect(range(dom[1], stop=dom[2], length=201)), x_train)))\n",
    "    X = [[x] for x in collect(x_arr)]\n",
    "    y = (x -> f_true(x[1])).(X)\n",
    "\n",
    "    p_true = Plots.Linear(x_arr, y, style=\"solid, black, mark=none\", legendentry=\"objective function\")\n",
    "\n",
    "    # create training set\n",
    "    Random.seed!(0)\n",
    "    GP = GaussianProcess(ν = 0.01)\n",
    "    for x in x_train\n",
    "        push!(GP, [x], f_true(x) + randn()*sqrt(GP.ν))\n",
    "    end\n",
    "    p_train = plot_GP_data(GP, legendentry=\"fit points\")\n",
    "\n",
    "    # predict\n",
    "    (μₚ, νₚ) = predict(GP, X)\n",
    "    σₚ = sqrt.(νₚ)\n",
    "    p_pred_μ = Plots.Linear(x_arr, μₚ, style=\"solid, thick, pastelBlue, mark=none\", legendentry=\"predicted\")\n",
    "    upperConfidence = μₚ + σ_safe*σₚ\n",
    "    lowerConfidence = μₚ - σ_safe*σₚ\n",
    "    p_pred_hi = Plots.Linear(x_arr, upperConfidence, style=\"draw=none, mark=none, name path=A, forget plot\")\n",
    "    p_pred_lo = Plots.Linear(x_arr, lowerConfidence, style=\"draw=none, mark=none, name path=B, forget plot\")\n",
    "    p_pred = Plots.Command(\"\\\\addplot[pastelBlue!40] fill between[of=A and B];\\\\addlegendentry{confidence region}\")\n",
    "\n",
    "    p = Plots.Plot[]\n",
    "    push!(p, p_pred_hi)\n",
    "    push!(p, p_pred_lo)\n",
    "    push!(p, p_pred)\n",
    "    push!(p, p_true)\n",
    "    push!(p, p_train)\n",
    "    push!(p, Plots.Linear([dom[1], dom[2]], [y_max, y_max], style=\"solid, thick, gray, mark=none\", legendentry=\"safety threshold\"))\n",
    "\n",
    "    ax3 = Axis(p, xlabel=L\"x\", ylabel=L\"y\", ymin=lo, ymax=hi, style=\"legend pos=outer north east, ytick=\\\\empty, legend style={draw=none}, axis on top=true, width=9.5cm, height=4cm, xlabel=\\$x\\$, xtick=\\\\empty, enlarge x limits=0\")\n",
    "\n",
    "    safe_regions = get_safe_regions(GP, β, dom..., y_max)\n",
    "    ax1 = Axis(plot_transparent_intervals(safe_regions, (lo,hi), \"green\", 0.5, \"estimated safe region \\$\\\\mathcal{S}\\$\"),\n",
    "            xmin=dom[1], xmax=dom[2], ymin=lo, ymax=hi, style=\"ytick=\\\\empty, xtick=\\\\empty, legend pos=outer north east, axis lines=none, legend style={draw=none}, width=9.5cm, height=2.05cm\")\n",
    "\n",
    "    d = (x,x′) -> norm(x-x′,2)\n",
    "    L = 4.0\n",
    "\n",
    "    E_regions = get_potential_expanders(GP, β, safe_regions, L, d, y_max)\n",
    "    ax2 = Axis(plot_transparent_intervals(E_regions, (lo,hi), \"orange\", 0.5, \"expected expanders \\$\\\\mathcal{E}\\$\"),\n",
    "            xmin=dom[1], xmax=dom[2], ymin=lo, ymax=hi, style=\"ytick=\\\\empty, xtick=\\\\empty, legend pos=outer north east, axis lines=none, legend style={draw=none}, width=9.5cm, height=2.05cm\")\n",
    "\n",
    "    G = GroupPlot(1, 3, groupStyle=\"vertical sep=0cm\")\n",
    "    push!(G, ax1)\n",
    "    push!(G, ax2)\n",
    "    push!(G, ax3)\n",
    "    G\n",
    "end\n",
    "plot(p)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "using Random\n",
    "\n",
    "p = let\n",
    "\tM = 9\n",
    "\tG = GroupPlot(1, 4M, groupStyle=\"vertical sep=0cm, xlabels at=edge bottom, xticklabels at=edge bottom\", style=\"xlabel=\\$x\\$, width=10cm\")\n",
    "\tregion_height = \"1.75cm\"\n",
    "\n",
    "\tP_safe = 0.9\n",
    "\tσ_safe = optimize(x->(cdf(Normal(0.0,1.0), x) - P_safe)^2, -5, 5).minimizer\n",
    "\tβ = σ_safe*σ_safe\n",
    "\ty_max = 0.8\n",
    "\n",
    "\tfunction f_(x::Float64, z::Float64=0.0)\n",
    "\t    A = Normal(1.0, 1.0)\n",
    "\t    B = Normal(3.0, 0.5)\n",
    "\t    C = Normal(5.0, 0.25)\n",
    "\t    return 1 - (2pdf(A, x+z) + pdf(B, x+z) + pdf(C, x+z))\n",
    "\tend\n",
    "\tf_true = x -> f_(x,0.0)\n",
    "\tdom = (-2,6)\n",
    "\tGP = GaussianProcess(ν = 0.01)\n",
    "\n",
    "\tlo, hi = -1.0, 1.2\n",
    "\n",
    "\td = (x,x′) -> norm(x-x′,2)\n",
    "\tL = 4.0\n",
    "\n",
    "\tx_train = [2.5]\n",
    "\tx_arr = unique(sort!(append!(collect(range(dom[1], stop=dom[2], length=201)), x_train)))\n",
    "\tX = [[x] for x in collect(x_arr)]\n",
    "\ty = (x -> f_true(x[1])).(X)\n",
    "\n",
    "\tp_true = Plots.Linear(x_arr, y, style=\"solid, black, mark=none\")\n",
    "\n",
    "\t# create training set\n",
    "\tRandom.seed!(0)\n",
    "\tfor x in x_train\n",
    "\t    push!(GP, [x], f_true(x) + randn()*sqrt(GP.ν))\n",
    "\tend\n",
    "\n",
    "\tfor i in 1 : M\n",
    "\t    p_train = plot_GP_data(GP)\n",
    "\n",
    "\t    # predict\n",
    "\t    (μₚ, νₚ) = predict(GP, X)\n",
    "\t    σₚ = sqrt.(νₚ)\n",
    "\t    p_pred_μ = Plots.Linear(x_arr, μₚ, style=\"solid, thick, pastelBlue, mark=none\")\n",
    "\t    upperConfidence = μₚ + σ_safe*σₚ\n",
    "\t    lowerConfidence = μₚ - σ_safe*σₚ\n",
    "\t    p_pred_hi = Plots.Linear(x_arr, upperConfidence, style=\"draw=none, mark=none, name path=A, forget plot\")\n",
    "\t    p_pred_lo = Plots.Linear(x_arr, lowerConfidence, style=\"draw=none, mark=none, name path=B, forget plot\")\n",
    "\t    p_pred = Plots.Command(\"\\\\addplot[pastelBlue!40, forget plot] fill between[of=A and B];\")\n",
    "\n",
    "\t    p = Plots.Plot[]\n",
    "\t    push!(p, p_pred_hi)\n",
    "\t    push!(p, p_pred_lo)\n",
    "\t    push!(p, p_pred)\n",
    "\t    push!(p, p_true)\n",
    "\t    push!(p, p_train)\n",
    "\t    push!(p, Plots.Linear([dom[1], dom[2]], [y_max, y_max], style=\"solid, gray, mark=none\"))\n",
    "\n",
    "\t    if i == 1\n",
    "\t        safe_regions = get_safe_regions(GP, β, dom..., y_max)\n",
    "\t        push!(G, Axis(plot_transparent_intervals(safe_regions, (lo,hi), \"green\", 0.5, \"estimated safe region\"),\n",
    "\t                      xmin=dom[1], xmax=dom[2], ymin=lo, ymax=hi,\n",
    "\t                      style=\"ytick=\\\\empty, xtick=\\\\empty, axis lines=none, height=$region_height, legend style={draw=none, at={(1.01,1)}, anchor=north west}\"))\n",
    "\n",
    "\t        M_regions = get_potential_maximizers(GP, β, safe_regions)\n",
    "\t        push!(G, Axis(plot_transparent_intervals(M_regions, (lo,hi), \"magenta\", 0.5, \"potential minimizers\"),\n",
    "\t                      xmin=dom[1], xmax=dom[2], ymin=lo, ymax=hi,\n",
    "\t                      style=\"ytick=\\\\empty, xtick=\\\\empty, axis lines=none, height=$region_height, legend style={draw=none, at={(1.01,0.5)}, anchor=north west}\"))\n",
    "\n",
    "\t        E_regions = get_potential_expanders(GP, β, safe_regions, L, d, y_max)\n",
    "\t        push!(G, Axis(plot_transparent_intervals(E_regions, (lo,hi), \"orange\", 0.5, \"potential expanders\"),\n",
    "\t                      xmin=dom[1], xmax=dom[2], ymin=lo, ymax=hi,\n",
    "\t                      style=\"ytick=\\\\empty, xtick=\\\\empty, axis lines=none, height=$region_height, legend style={draw=none, at={(1.01,0)}, anchor=north west}\"))\n",
    "\t    else\n",
    "\t        safe_regions = get_safe_regions(GP, β, dom..., y_max)\n",
    "\t        push!(G, Axis(plot_transparent_intervals(safe_regions, (lo,hi), \"green\", 0.5),\n",
    "\t                      xmin=dom[1], xmax=dom[2], ymin=lo, ymax=hi,\n",
    "\t                      style=\"ytick=\\\\empty, xtick=\\\\empty, axis lines=none, height=$region_height\"))\n",
    "\n",
    "\t        M_regions = get_potential_maximizers(GP, β, safe_regions)\n",
    "\t        push!(G, Axis(plot_transparent_intervals(M_regions, (lo,hi), \"magenta\", 0.5),\n",
    "\t                      xmin=dom[1], xmax=dom[2], ymin=lo, ymax=hi,\n",
    "\t                      style=\"ytick=\\\\empty, xtick=\\\\empty, axis lines=none, height=$region_height\"))\n",
    "\n",
    "\t        E_regions = get_potential_expanders(GP, β, safe_regions, L, d, y_max)\n",
    "\t        push!(G, Axis(plot_transparent_intervals(E_regions, (lo,hi), \"orange\", 0.5),\n",
    "\t                      xmin=dom[1], xmax=dom[2], ymin=lo, ymax=hi,\n",
    "\t                      style=\"ytick=\\\\empty, xtick=\\\\empty, axis lines=none, height=$region_height\"))\n",
    "\t    end\n",
    "\n",
    "\t    x_sample = [NaN]\n",
    "\t    w_best = -Inf\n",
    "\t    Δx = 0.01\n",
    "\t    for (a,b) in vcat(M_regions, E_regions)\n",
    "\t        m = ceil(Int, (b-a)/Δx)\n",
    "\t        for x in range(a, stop=b, length=m)\n",
    "\t            w = width(GP, [x], β)\n",
    "\t            if w > w_best\n",
    "\t                w_best = w\n",
    "\t                x_sample[1] = x\n",
    "\t            end\n",
    "\t        end\n",
    "\t    end\n",
    "\t    push!(p, Plots.Scatter(x_sample, [f_true(x_sample[1])], style=\"only marks, black, mark=o, mark size=1.5, mark options={draw=black}\"))\n",
    "\n",
    "\t    push!(GP, x_sample, f_true(x_sample[1]) + randn()*sqrt(GP.ν))\n",
    "\n",
    "\t    ax = Axis(p, ylabel=L\"y\", ymin=lo, ymax=hi, style=\"ytick=\\\\empty, axis on top=true, height=2.5cm, xtick=\\\\empty, enlarge x limits=0\")\n",
    "\t    push!(G, ax)\n",
    "\tend\n",
    "\tG\n",
    "end\n",
    "plot(p)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "using Colors\n",
    "using Random\n",
    "\n",
    "p = let\n",
    "\tG = GroupPlot(4, 4, groupStyle=\"horizontal sep=0.25cm, vertical sep=0.25cm, ylabels at=edge left, xlabels at=edge bottom, xticklabels at=edge bottom, yticklabels at=edge left\",\n",
    "\t    style=\"xlabel=\\$x_1\\$, ylabel=\\$x_2\\$, width=4cm, height=4cm, contour/labels=false, view={0}{90}\"\n",
    "\t)\n",
    "\n",
    "\tfunction flower(x; a=1, b=1, c=4)::Float64\n",
    "\t    if isapprox(norm(x), 0.0)\n",
    "\t        return 0.0\n",
    "\t    end\n",
    "\t\treturn a*norm(x) + b*sin(c*atan(x[2], x[1]))\n",
    "\tend\n",
    "\tf_(x::Vector{Float64}) = flower(x)\n",
    "\tf_(x1::Float64, x2::Float64) = f_([x1,x2])\n",
    "\n",
    "\txdomain = ( -3, 3)\n",
    "\tydomain = ( -3, 3)\n",
    "\n",
    "\ty_max = 2.0\n",
    "\n",
    "\tfunction plot_GP_mean(GP, β)\n",
    "\t    m, k, ν = GP.m, GP.k, GP.ν\n",
    "\t    iK = inv(K(GP.X, GP.X, k) + ν*I)\n",
    "\n",
    "\t    u_pred = (x,y) -> begin\n",
    "\t                X_pred = Vector{Float64}[[x,y]]\n",
    "\t                tmp = K(X_pred, GP.X, k) * iK\n",
    "\t                μₚ = μ(X_pred, m) + tmp*(GP.y - μ(GP.X, m))\n",
    "\t                νₚ = diag(K(X_pred, X_pred, k) + ν*I - tmp*K(GP.X, X_pred, k))\n",
    "\t                return μₚ[1] + sqrt(β*νₚ[1])\n",
    "\t            end\n",
    "\n",
    "\t    plots = Plots.Plot[]\n",
    "\t    push!(plots, Plots.Image(u_pred, xdomain, ydomain, xbins=600, ybins=600, colormap=pasteljet, colorbar = false))\n",
    "\t    push!(plots, Plots.Contour(u_pred, xdomain, ydomain, levels=[y_max], contour_style=\"draw color=white\", style=\"white\", xbins=101, ybins=101))\n",
    "\t    push!(plots, Plots.Scatter([x[1] for x in GP.X], [x[2] for x in GP.X], style=\"only marks, white, mark=*, mark size=0.75, mark options={draw=white, fill=white}\"))\n",
    "\n",
    "\t    Axis(plots, xmin=xdomain[1], xmax=xdomain[2], ymin=ydomain[1], ymax=ydomain[2])\n",
    "\tend\n",
    "\n",
    "\tGP = GaussianProcess(m=x->y_max + 0.5, ν = 0.01)\n",
    "\tβ = 10.0\n",
    "\tn = 51\n",
    "\tX = Array{Vector{Float64}}(undef, n*n)\n",
    "\ti = 0\n",
    "\tfor x1 in range(xdomain[1], stop=xdomain[2], length=n)\n",
    "\t    for x2 in range(ydomain[1], stop=ydomain[2], length=n)\n",
    "\t        X[i+=1] = [x1,x2]\n",
    "\t    end\n",
    "\tend\n",
    "\ti = argmin([norm([-2,1]-x,2) for x in X])\n",
    "\n",
    "\tm = length(X)\n",
    "\tu, ℓ = fill(Inf, m), fill(-Inf, m)\n",
    "\tS, M, E = falses(m), falses(m), falses(m)\n",
    "\n",
    "\tRandom.seed!(0)\n",
    "\tpush!(GP, X[i], f_(X[i]) + randn()*GP.ν)\n",
    "\tupdate_confidence_intervals!(GP, X, u, ℓ, β);\n",
    "\tS[:] = u .≤ y_max\n",
    "\tbest_val, i_best = findmin(u[S])\n",
    "\ti_best = something(findfirst(isequal(i_best), cumsum(S)), 0)\n",
    "\n",
    "\tax = plot_GP_mean(GP, β)\n",
    "\tpush!(ax, Plots.Scatter([X[i_best][1]], [X[i_best][2]], style=\"only marks, white, mark=+, mark size=3, mark options={draw=white, fill=white}\"))\n",
    "\tpush!(G, ax)\n",
    "\n",
    "\tfor k in 2 : prod(G.dimensions)\n",
    "\t    compute_sets!(GP, S, M, E, X, u, ℓ, y_max, β)\n",
    "\t    i = get_new_query_point(M, E, u, ℓ)\n",
    "\n",
    "\t    push!(GP, X[i], f_(X[i]))\n",
    "\t    update_confidence_intervals!(GP, X, u, ℓ, β);\n",
    "\t    S[:] = u .≤ y_max\n",
    "\t    best_val, i_best = findmin(u[S])\n",
    "\t    i_best = something(findfirst(isequal(i_best), cumsum(S)), 0)\n",
    "\n",
    "\n",
    "\t    ax = plot_GP_mean(GP, β)\n",
    "\t    push!(ax, Plots.Scatter([X[i_best][1]], [X[i_best][2]], style=\"only marks, white, mark=+, mark size=3, mark options={draw=white, fill=white}\"))\n",
    "\t    push!(G, ax)\n",
    "\tend\n",
    "\n",
    "\tG\n",
    "end\n",
    "plot(p)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p = let\n",
    "function flower(x; a=1, b=1, c=4)::Float64\n",
    "\tif isapprox(norm(x), 0.0)\n",
    "\t\treturn 0.0\n",
    "\tend\n",
    "\treturn a*norm(x) + b*sin(c*atan(x[2], x[1]))\n",
    "end\n",
    "\n",
    "y_max = 2.0\n",
    "xdomain = ( -3, 3)\n",
    "ydomain = ( -3, 3)\n",
    "\n",
    "plots = Plots.Plot[]\n",
    "push!(plots, Plots.Image((x,y)->flower([x,y]), xdomain, ydomain, xbins=600, ybins=600, colormap=pasteljet, colorbar = false))\n",
    "push!(plots, Plots.Contour((x,y)->flower([x,y]), xdomain, ydomain, levels=[y_max], contour_style=\"draw color=white\", style=\"white\", xbins=101, ybins=101))\n",
    "Axis(plots, width=\"4cm\", height=\"4cm\", xlabel=L\"x_1\", ylabel=L\"x_2\", style=\"contour/labels=false, axis equal, view={0}{90}\")\n",
    "end\n",
    "plot(p)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "@webio": {
   "lastCommId": null,
   "lastKernelId": null
  },
  "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
}