sampling-plans.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Sampling Plans\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": [
    "using Random\n",
    "p = let\n",
    "\tL = 6π\n",
    "\tf = x -> (cos(x) - 1)*min(abs(L - x), abs(-x - L))\n",
    "\tX = collect(-L:2π:L)\n",
    "\tRandom.seed!(4)\n",
    "\tX2 = X + rand(length(X))*2π .- π\n",
    "\tplots = Plots.Plot[]\n",
    "\tpush!(plots, Plots.Linear(f, (-L, L), xbins=151, style=\"solid, black, mark=none\", legendentry=L\"f(x)\"))\n",
    "\tpush!(plots, Plots.Linear(X, f.(X), style=\"only marks, mark=*, mark size=1, mark options={draw=pastelBlue, fill=pastelBlue}\", legendentry=\"sampling on grid\"))\n",
    "\tpush!(plots, Plots.Linear(X2, f.(X2), style=\"only marks, mark=*, mark size=1, mark options={draw=pastelRed, fill=pastelRed}\", legendentry=\"stratified sampling\"))\n",
    "\tAxis(plots, width=\"9cm\", xlabel=L\"x\", ylabel=L\"y\", style=\"xtick=\\\\empty, ytick=\\\\empty, legend cell align=left, legend style={draw=none, at={(0.5,-0.15)}, anchor=north, legend columns=1}\")\n",
    "end\n",
    "plot(p)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "using Random\n",
    "using LinearAlgebra\n",
    "\n",
    "function mutate!(Xs)\n",
    "    m, n = length(Xs), length(Xs[1])\n",
    "    j = rand(1:n)\n",
    "    i = randperm(m)[1:2]\n",
    "    Xs[i[1]][j], Xs[i[2]][j] = Xs[i[2]][j], Xs[i[1]][j]\n",
    "    return Xs\n",
    "end\n",
    "function uniform_projection_plan(n)\n",
    "\tp = randperm(n)\n",
    "\t[[u,v] for (u,v) in zip(p,1:n)]\n",
    "end\n",
    "function pairwise_distances(Xs, p=2)\n",
    "    n = length(Xs)\n",
    "    [norm(Xs[i]-Xs[j], p) for i in 1:(n-1) for j in (i+1):n]\n",
    "end\n",
    "phiq(Xs, q, p) = sum(d^(-q) for d in pairwise_distances(Xs, p))^(1/q)\n",
    "\n",
    "p = let\n",
    "\n",
    "\tN = 20\n",
    "\tcapacity = 6*3\n",
    "\n",
    "\tRandom.seed!(0)\n",
    "\tX = Vector{Vector{Float64}}[]\n",
    "\tpush!(X, [Float64[i,i] for i in 1:N])\n",
    "\tfor i in 2 : capacity\n",
    "\t\tpush!(X, mutate!(mutate!(deepcopy(X[end]))))\n",
    "\tend\n",
    "\tΦ = [phiq(x, 1, 2) for x in X]\n",
    "\n",
    "\tfor i in 1 : 50\n",
    "\t\t# remove the index with minimum distance to neighbors\n",
    "\t\ti_min = 1\n",
    "\t\td_min = Inf\n",
    "\t\tp = sortperm(Φ)\n",
    "\t\tX[:] = X[p]\n",
    "\t\tΦ[:] = Φ[p]\n",
    "\t\tfor i in 2 : length(X)-1\n",
    "\t\t\td = Φ[i+1] - Φ[i-1]\n",
    "\t\t\t@assert d >= 0\n",
    "\t\t\tif d < d_min\n",
    "\t\t\t\td_min = d\n",
    "\t\t\t\ti_min = 1\n",
    "\t\t\tend\n",
    "\t\tend\n",
    "\t\tX[i_min] = mutate!(X[i_min])\n",
    "\t\tΦ[i_min] = phiq(X[i_min], 1, 2)\n",
    "\tend\n",
    "\n",
    "\tsort!(X, by=x->phiq(x, 1, 2))\n",
    "\n",
    "    g = GroupPlot(6,3,groupStyle=\"horizontal sep=0.5cm, vertical sep=0.75cm\")\n",
    "\n",
    "   \tfor Y in X\n",
    "   \t\tpush!(g, Axis(\n",
    "   \t\t\tPlots.Scatter([x[1] for x in Y], [x[2] for x in Y], style=\"mark=*, mark size=0.7, mark options={draw=black, fill=black}\"),\n",
    "   \t\t\twidth=\"4cm\", height=\"4cm\", style=\"ytick=\\\\empty, xtick=\\\\empty, title style={yshift=-0.5em}\",\n",
    "   \t\t\ttitle = \"\\$\\\\Phi_1 = $(round(phiq(Y, 1, 2), digits=1))\\$\", # $\n",
    "   \t\t\t)\n",
    "\t\t)\n",
    "   \tend\n",
    "\n",
    "    g\n",
    "end\n",
    "\n",
    "plot(p)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "using Random\n",
    "\n",
    "min_dist(a, B, p) = minimum(norm(a-b, p) for b in B)\n",
    "d_max(A, B, p=2) = maximum(min_dist(a, B, p) for a in A)\n",
    "function greedy_local_search(X, m, d=d_max)\n",
    "\tS = [X[rand(1:m)]]\n",
    "\tfor i in 2 : m\n",
    "\t\tj = argmin([x ∈ S ? Inf : d(X, push!(copy(S), x))\n",
    "\t\t           for x in X])\n",
    "\t\tpush!(S, X[j])\n",
    "\tend\n",
    "\treturn S\n",
    "end\n",
    "function exchange_algorithm(X, m, d=d_max)\n",
    "\tS = X[randperm(length(X))[1:m]]\n",
    "\tδ, done = d(X, S), false\n",
    "\twhile !done\n",
    "\t\tbest_pair = (0,0)\n",
    "\t\tfor i in 1 : m\n",
    "\t\t\ts = S[i]\n",
    "\t\t\tfor (j,x) in enumerate(X)\n",
    "\t\t\t\tif !in(x, S)\n",
    "\t\t\t\t\tS[i] = x\n",
    "\t\t\t\t\tδ′ = d(X, S)\n",
    "\t\t\t\t\tif δ′ < δ\n",
    "\t\t\t\t\t\tδ = δ′\n",
    "\t\t\t\t\t\tbest_pair = (i,j)\n",
    "\t\t\t\t\tend\n",
    "\t\t\t\tend\n",
    "\t\t\tend\n",
    "\t\t\tS[i] = s\n",
    "\t\tend\n",
    "\t\tdone = best_pair == (0,0)\n",
    "\t\tif !done\n",
    "\t\t\ti,j = best_pair\n",
    "\t\t\tS[i] = X[j]\n",
    "\t\tend\n",
    "\tend\n",
    "\treturn S\n",
    "end\n",
    "function multistart_local_search(X, m, alg, k_max, d=d_max)\n",
    "\tsets = [alg(X, m, d) for i in 1 : k_max]\n",
    "\treturn sets[argmin([d(X, S) for S in sets])]\n",
    "end\n",
    "\n",
    "p = let\n",
    "\n",
    "\tRandom.seed!(0)\n",
    "\tX = Vector{Float64}[]\n",
    "\tfor i in 1 : 100\n",
    "\t\tpush!(X, rand(2))\n",
    "\tend\n",
    "\n",
    "\tXs_gls = multistart_local_search(X, 8, greedy_local_search, 10)\n",
    "\tXs_ea = multistart_local_search(X, 8, exchange_algorithm, 10)\n",
    "\n",
    "\tpX = Plots.Scatter([x[1] for x in X], [x[2] for x in X], style=\"mark=*, mark size=1, mark options={draw=black, fill=black}\")\n",
    "\n",
    "    g = GroupPlot(2,1,groupStyle=\"horizontal sep=0.5cm\")\n",
    "\n",
    "    push!(g, Axis(\n",
    "    \t[pX,\n",
    "    \t Plots.Scatter([x[1] for x in Xs_gls], [x[2] for x in Xs_gls], style=\"mark=*, mark size=1.5, mark options={draw=pastelRed, fill=pastelRed}\")\n",
    "    \t],\n",
    "    \twidth=\"6.25cm\", xlabel=L\"x_1\", ylabel=L\"x_2\", title=\"greedy local search\",\n",
    "    \tstyle=\"ytick=\\\\empty, xtick=\\\\empty\",\n",
    "    ))\n",
    "\n",
    "    push!(g, Axis(\n",
    "    \t[pX,\n",
    "    \t Plots.Scatter([x[1] for x in Xs_ea], [x[2] for x in Xs_ea], style=\"mark=*, mark size=1.5, mark options={draw=pastelRed, fill=pastelRed}\")\n",
    "    \t],\n",
    "    \twidth=\"6.25cm\", xlabel=L\"x_1\", title=\"exchange algorithm\",\n",
    "    \tstyle=\"ytick=\\\\empty, xtick=\\\\empty\",\n",
    "    ))\n",
    "\n",
    "\n",
    "    g\n",
    "end\n",
    "\n",
    "plot(p)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "using Sobol\n",
    "import OnlineStats: Mean, Series, value\n",
    "using Random\n",
    "\n",
    "p = let\n",
    "\n",
    "\tf = x -> sin(10x)\n",
    "\n",
    "\n",
    "\tRandom.seed!(0)\n",
    "\tn_tries_per = 4\n",
    "\tf_true = sin(5)^2 / 5\n",
    "\n",
    "\tarr_n = [round(Int, n) for n in 10 .^ range(0, stop=4, length=101)]\n",
    "\tarr_err_rand = map(arr_n) do n\n",
    "\t    o = Series(Mean())\n",
    "\t    for i in 1 : n_tries_per\n",
    "\t        v = Series(Mean())\n",
    "\t        for j in 1 : n\n",
    "\t            fit!(v, f(rand()))\n",
    "\t        end\n",
    "\t        fit!(o, value(v)[1])\n",
    "\t    end\n",
    "\t    abs(value(o)[1] - f_true) / f_true\n",
    "\tend\n",
    "\tarr_err_sobol = map(arr_n) do n\n",
    "\t    o = Series(Mean())\n",
    "\t    for x in Iterators.take(SobolSeq(1), n)\n",
    "\t        fit!(o, f(x[1]))\n",
    "\t    end\n",
    "\t    abs(value(o)[1] - f_true) / f_true\n",
    "\tend\n",
    "\n",
    "\tAxis([\n",
    "\t    Plots.Linear(arr_n, arr_err_rand, style=\"mark=none, pastelRed\", legendentry=\"rand\"),\n",
    "\t    Plots.Linear(arr_n, arr_err_sobol, style=\"mark=none, pastelBlue\", legendentry=\"Sobol\"),\n",
    "\t    ], xlabel=\"number of samples\", ylabel=\"relative error\", xmode=\"log\", ymode=\"log\", width=\"9cm\",\n",
    "\t       style=\"trim axis left, trim axis right, legend cell align=left, legend style={draw=none, at={(0.5,-0.4)}, anchor=north, legend columns=2}\")\n",
    "end\n",
    "\n",
    "plot(p)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "using Primes\n",
    "function halton(i, b)\n",
    "    result, f = 0.0, 1.0\n",
    "    while i > 0\n",
    "        f = f / b;\n",
    "        result = result + f * mod(i, b)\n",
    "        i = floor(Int, i / b)\n",
    "    end\n",
    "    return result\n",
    "end\n",
    "get_filling_set_halton(m; b=2) = [halton(i,b) for i in 1: m]\n",
    "function get_filling_set_halton(m, n)\n",
    "    bs = primes(max(ceil(Int, n*(log(n) + log(log(n)))), 6))\n",
    "    seqs = [get_filling_set_halton(m, b=b) for b in bs[1:n]]\n",
    "    return [collect(x) for x in zip(seqs...)]\n",
    "end\n",
    "\n",
    "p = let\n",
    "\n",
    "\tn = 19\n",
    "\tseqs = [get_filling_set_halton(n, b=b) for b in [19,23]]\n",
    "\tX = [collect(x) for x in zip(seqs...)]\n",
    "\n",
    "\tAxis(\n",
    "\t    Plots.Scatter([x[1] for x in X], [x[2] for x in X],\n",
    "\t                      style=\"mark=*, mark size=1, mark options={draw=black, fill=black}\"),\n",
    "\t    width=\"5cm\", height=\"5cm\",\n",
    "\tstyle=\"xtick=\\\\empty, ytick=\\\\empty, clip marker paths=true, axis on top=true\"\n",
    "\t    )\n",
    "end\n",
    "\n",
    "plot(p)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "using Primes\n",
    "using Random\n",
    "import Printf: @sprintf\n",
    "\n",
    "function get_filling_set_random_walk(n)\n",
    "    Xs = [rand()]\n",
    "    for i in 2 : n\n",
    "        push!(Xs, mod(Xs[end] + 0.5 + rand(), 1))\n",
    "    end\n",
    "    return Xs\n",
    "end\n",
    "\n",
    "function get_filling_set_additive_recurrence(n; c=φ-1)\n",
    "    Xs = [rand()]\n",
    "    for i in 2 : n\n",
    "        push!(Xs, mod(Xs[end] + c, 1))\n",
    "    end\n",
    "    return Xs\n",
    "end\n",
    "function get_filling_set_additive_recurrence(n, k)\n",
    "    ps = primes(max(ceil(Int, k*(log(k) + log(log(k)))), 6))\n",
    "    seqs = [get_filling_set_additive_recurrence(n, c=sqrt(p))\n",
    "            for p in ps[1:k]]\n",
    "    return [collect(x) for x in zip(seqs...)]\n",
    "end\n",
    "\n",
    "p = let\n",
    "\n",
    "\tg = GroupPlot(5,3, groupStyle=\"horizontal sep=0.25cm, vertical sep=0.25cm\")\n",
    "\n",
    "\tns = [10,100,1000]\n",
    "\n",
    "\tfunction add_plot(Xs,mark_size)\n",
    "\t\tpush!(g, Axis(\n",
    "\t\t    Plots.Scatter([x[1] for x in Xs], [x[2] for x in Xs], [i/length(Xs) for i in 1 : length(Xs)],\n",
    "\t\t                      style=\"mark=*, mark size=$(mark_size), scatter/use mapped color={draw=mapped color, fill=mapped color}, colormap name=pasteljet\"),\n",
    "\t\t    width=\"4.5cm\", height=\"4.5cm\", enlargelimits=0,\n",
    "\t\tstyle=\"xtick=\\\\empty, ytick=\\\\empty, clip marker paths=true, axis on top=true\"\n",
    "\t\t    ))\n",
    "\tend\n",
    "\n",
    "\tRandom.seed!(0)\n",
    "\tmark_sizes = [1,0.75,0.5]\n",
    "\tfor (j,n) in enumerate(ns)\n",
    "\t\tmark_size = mark_sizes[j]\n",
    "\t    add_plot([rand(2) for i in 1 : n], mark_size)\n",
    "\t    g.axes[end].ylabel=@sprintf(\"\\\\small \\$m = %d\\$\", n)\n",
    "\n",
    "\t    add_plot(get_filling_set_additive_recurrence(n, 2), mark_size)\n",
    "\t    add_plot(get_filling_set_halton(n, 2), mark_size)\n",
    "\n",
    "\t    Xs = Array{Vector{Float64}}(undef, n)\n",
    "\t    for (i,x) in enumerate(Iterators.take(SobolSeq(2), n))\n",
    "\t        Xs[i] = x\n",
    "\t    end\n",
    "\t    add_plot(Xs, mark_size)\n",
    "\n",
    "\t    add_plot(uniform_projection_plan(n)./n, mark_size)\n",
    "\tend\n",
    "\n",
    "\tg.axes[end-4].xlabel=\"\\\\small Random\"\n",
    "\tg.axes[end-3].xlabel=\"\\\\small Additive Recurrence\"\n",
    "\tg.axes[end-2].xlabel=\"\\\\small Halton\"\n",
    "\tg.axes[end-1].xlabel=\"\\\\small Sobol\"\n",
    "\tg.axes[end].xlabel=\"\\\\small Uniform Projection\"\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
}