Gauss-Newton.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Implementação do método Gauss-Newton com NLPModels e aplicação nos dados covid.csv"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "using NLPModels"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "ADNLSModel(Minimization problem Generic\n",
       "nvar = 2, ncon = 0 (0 linear)\n",
       ", NLSMeta(2, 2, [0.0, 0.0], 4, 3), NLSCounters(Counters(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), 0, 0, 0, 0, 0, 0, 0), var\"#3#4\"(), NLPModels.var\"#6#8\"())"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "nlp = ADNLSModel(\n",
    "    x -> [x[1] * x[2] - 1.0; x[1]^2 + x[2]^2 - 4.0],\n",
    "    [0.5; 0.7], 2\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2-element Array{Float64,1}:\n",
       " -0.65\n",
       " -3.26"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "residual(nlp, nlp.meta.x0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2×2 Array{Float64,2}:\n",
       " 0.7  0.5\n",
       " 1.0  1.4"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "jac_residual(nlp, nlp.meta.x0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "J (generic function with 1 method)"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "F(x) = residual(nlp, x)\n",
    "J(x) = jac_residual(nlp, x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "([0.5, 0.7], 3.324169069105842, 6.140321327096815)"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "using LinearAlgebra\n",
    "\n",
    "x = copy(nlp.meta.x0)\n",
    "x, norm(F(x)), norm(J(x)' * F(x))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "([0.5176380902050415, 1.9318516525781366], 1.1102230246251565e-16, 2.220446049250313e-16)"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "d = -(J(x)' * J(x)) \\ (J(x)' * F(x))\n",
    "x = x + d\n",
    "x, norm(F(x)), norm(J(x)' * F(x))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2-element Array{Float64,1}:\n",
       " 0.5\n",
       " 0.7"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x = copy(nlp.meta.x0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "([0.5176380902050414, 1.9318516525781368], 8.95090418262362e-16, 3.4470592372182254e-15)"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "using Krylov\n",
    "\n",
    "d, krylov_stats = Krylov.lsmr(J(x), -F(x))\n",
    "x = x + d\n",
    "x, norm(F(x)), norm(J(x)' * F(x))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 385,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "gauss_newton (generic function with 1 method)"
      ]
     },
     "execution_count": 385,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "using Krylov, NLPModels, SolverTools # JuliaSmoothOptimizers\n",
    "\n",
    "function gauss_newton(nls :: AbstractNLSModel;\n",
    "        α = 0.1,\n",
    "        atol = 1e-6, rtol = 1e-6, max_eval = 10_000, max_time = 60.0)\n",
    "    x = copy(nls.meta.x0)\n",
    "    F(x) = residual(nls, x)\n",
    "    J(x) = jac_residual(nls, x)\n",
    "    \n",
    "    Fx = F(x)\n",
    "    Jx = J(x)\n",
    "    JᵀF = Jx' * Fx\n",
    "    \n",
    "    ϵ = atol + rtol * norm(JᵀF)\n",
    "    start_time = time()\n",
    "    Δt = 0.0\n",
    "    \n",
    "    solved = norm(JᵀF) < ϵ\n",
    "    tired = neval_residual(nls) > max_eval || Δt > max_time\n",
    "    \n",
    "    while !(solved || tired)\n",
    "        d, krylov_stats = lsmr(Jx, -Fx)\n",
    "        x += α * d\n",
    "        \n",
    "        Fx = F(x)\n",
    "        Jx = J(x)\n",
    "        JᵀF = Jx' * Fx\n",
    "        \n",
    "        Δt = time() - start_time\n",
    "        solved = norm(JᵀF) < ϵ\n",
    "        tired = neval_residual(nls) > max_eval || Δt > max_time\n",
    "    end\n",
    "    \n",
    "    status = :unknown\n",
    "    if solved\n",
    "        status = :first_order\n",
    "    elseif tired\n",
    "        if Δt > max_time\n",
    "            status = :max_time\n",
    "        else\n",
    "            status = :max_eval\n",
    "        end\n",
    "    end\n",
    "    \n",
    "    return GenericExecutionStats(status, nls, solution=x, dual_feas=norm(JᵀF),\n",
    "                                 elapsed_time=Δt, objective=dot(Fx, Fx) / 2)\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 384,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "STATUSES:\n",
      "  :acceptable     => solved to within acceptable tolerances\n",
      "  :exception      => unhandled exception\n",
      "  :first_order    => first-order stationary\n",
      "  :infeasible     => problem may be infeasible\n",
      "  :max_eval       => maximum number of function evaluations\n",
      "  :max_iter       => maximum iteration\n",
      "  :max_time       => maximum elapsed time\n",
      "  :neg_pred       => negative predicted reduction\n",
      "  :not_desc       => not a descent direction\n",
      "  :small_residual => small residual\n",
      "  :small_step     => step too small\n",
      "  :stalled        => stalled\n",
      "  :unbounded      => objective function may be unbounded from below\n",
      "  :unknown        => unknown\n",
      "  :user           => user-requested stop\n"
     ]
    }
   ],
   "source": [
    "SolverTools.show_statuses()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 386,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\"Execution stats: first-order stationary\""
      ]
     },
     "execution_count": 386,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "nls = ADNLSModel(\n",
    "    x -> [x[1] * x[2] - 1.0; x[1]^2 + x[2]^2 - 4.0],\n",
    "    [0.5; 0.7], 2\n",
    ")\n",
    "\n",
    "output = gauss_newton(nls)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 387,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Generic Execution stats\n",
      "  status: \"first-order stationary\"\n",
      "  objective value: 1.3742013128069374e-12\n",
      "  primal feasibility: 0.0\n",
      "  dual feasibility: 6.864157090622998e-6\n",
      "  solution: [0.5176380043262384  1.9318512573626214]\n",
      "  iterations: -1\n",
      "  elapsed time: 0.0006451606750488281\n"
     ]
    }
   ],
   "source": [
    "print(output)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 399,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Generic Execution stats\n",
      "  status: \"first-order stationary\"\n",
      "  objective value: 3.4910339246502494e-25\n",
      "  primal feasibility: 0.0\n",
      "  dual feasibility: 1.8699247906322293e-11\n",
      "  solution: [0.9999999999999658  1.000000000000015]\n",
      "  iterations: -1\n",
      "  elapsed time: 5.602836608886719e-5\n"
     ]
    }
   ],
   "source": [
    "nls = ADNLSModel(\n",
    "    x -> [x[1] - 1.0; 10 * (x[2] - x[1]^2)],\n",
    "    [-1.2; 1.0], 2\n",
    ")\n",
    "\n",
    "output = gauss_newton(nls, α = 1.0)\n",
    "print(output)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 400,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3"
      ]
     },
     "execution_count": 400,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "neval_residual(nls)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 395,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Generic Execution stats\n",
      "  status: \"first-order stationary\"\n",
      "  objective value: 0.4489276974679337\n",
      "  primal feasibility: 0.0\n",
      "  dual feasibility: 1.9008628166288058e-13\n",
      "  solution: [-0.07005677942035861  0.9114550457707855  -0.047297441804754106]\n",
      "  iterations: -1\n",
      "  elapsed time: 7.295608520507812e-5\n"
     ]
    }
   ],
   "source": [
    "A = rand(10, 3)\n",
    "b = rand(10)\n",
    "nls = ADNLSModel(x -> A * x - b, zeros(3), 10)\n",
    "\n",
    "output = gauss_newton(nls, α = 1.0)\n",
    "print(output)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 396,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2"
      ]
     },
     "execution_count": 396,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "neval_residual(nls)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Dados do COVID"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 401,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "ADNLSModel(Minimization problem Generic\n",
       "nvar = 3, ncon = 0 (0 linear)\n",
       ", NLSMeta(25, 3, [0.0, 0.0, 0.0], 75, 6), NLSCounters(Counters(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), 0, 0, 0, 0, 0, 0, 0), var\"#54#56\"(), NLPModels.var\"#6#8\"())"
      ]
     },
     "execution_count": 401,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "using CSV, DataFrames\n",
    "\n",
    "data = CSV.read(\"covid.csv\")\n",
    "rename!(data, Dict(\n",
    "        Symbol(\"Countries and territories\") => \"CountryExp\",\n",
    "        :Cases => \"NewConfCases\"\n",
    "    ))\n",
    "df = data[findall(data[!,:CountryExp] .== \"Italy\"),:]\n",
    "y = reverse(df.NewConfCases)\n",
    "y = cumsum(y[findlast(y .== 0)+1:end])\n",
    "\n",
    "n = length(y)\n",
    "x = range(0, 1, length=n)\n",
    "\n",
    "h(β, x) = β[3] / (1 + exp(-β[1] - β[2] * x))\n",
    "nls = ADNLSModel(\n",
    "    β -> [y[i] - h(β, x[i]) for i = 1:n],\n",
    "    [-5.53; 5.76; 4e4], n\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 402,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "25-element Array{Float64,1}:\n",
       " -144.01288618111252\n",
       " -124.65873714096207\n",
       " -125.74062027262167\n",
       "  -97.28034501301562\n",
       "  -91.07089093487843\n",
       " -122.85666231413438\n",
       "  -11.54579054625151\n",
       "   51.54621925278343\n",
       "   71.4272149886649 \n",
       "  356.1476464210357 \n",
       "  156.53571768829215\n",
       "  392.995262587885  \n",
       "  446.4396203512497 \n",
       "  558.4672748834896 \n",
       "  533.9186622361076 \n",
       "  809.9767951143713 \n",
       " 1140.9663950501063 \n",
       " 1568.9373593874816 \n",
       "  957.9579618640746 \n",
       " 1463.7799115929283 \n",
       " 2102.2109032925455 \n",
       " 2461.257293723609  \n",
       "  234.0600354319322 \n",
       " 4076.9991666750066 \n",
       " 5687.085816057781  "
      ]
     },
     "execution_count": 402,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "residual(nls, nls.meta.x0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 403,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3-element Array{Float64,1}:\n",
       "    -5.53\n",
       "     5.76\n",
       " 40000.0 "
      ]
     },
     "execution_count": 403,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "F(β) = residual(nls, β)\n",
    "J(β) = jac_residual(nls, β)\n",
    "β = copy(nls.meta.x0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 408,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Generic Execution stats\n",
      "  status: \"maximum number of function evaluations\"\n",
      "  objective value: 4.4688635649275e6\n",
      "  primal feasibility: 0.0\n",
      "  dual feasibility: 1.0873342719422732\n",
      "  solution: [-5.73117604392885  5.057879880801879  80835.9720485137]\n",
      "  iterations: -1\n",
      "  elapsed time: 0.19988489151000977\n"
     ]
    }
   ],
   "source": [
    "reset!(nls)\n",
    "output = gauss_newton(nls, atol=1e-12, rtol=1e-12)\n",
    "β = output.solution\n",
    "print(output)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 409,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "25-element Array{Float64,1}:\n",
       "   261.3245933691231 \n",
       "   322.3861880646415 \n",
       "   397.6451364527162 \n",
       "   490.36579268103856\n",
       "   604.5440416377118 \n",
       "   745.0613643564338 \n",
       "   917.8662093968711 \n",
       "  1130.184753411804  \n",
       "  1390.7616176283746 \n",
       "  1710.1283976379402 \n",
       "  2100.893411964594  \n",
       "  2578.0391849980288 \n",
       "  3159.2040500598127 \n",
       "  3864.9101078317285 \n",
       "  4718.681129008411  \n",
       "  5746.971207154713  \n",
       "  6978.800096130133  \n",
       "  8444.969125400341  \n",
       " 10176.721463476424  \n",
       " 12203.72631150056   \n",
       " 14551.326246363134  \n",
       " 17237.107448626346  \n",
       " 20267.040224236443  \n",
       " 23631.673118209077  \n",
       " 27303.09067090625   "
      ]
     },
     "execution_count": 409,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "y_pred = [h(β,xi) for xi in x]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 410,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Plots.GRBackend()"
      ]
     },
     "execution_count": 410,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "using Plots\n",
    "gr()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 411,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/svg+xml": [
       "<?xml version=\"1.0\" encoding=\"utf-8\"?>\n",
       "<svg xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\" width=\"600\" height=\"400\" viewBox=\"0 0 2400 1600\">\n",
       "<defs>\n",
       "  <clipPath id=\"clip2400\">\n",
       "    <rect x=\"0\" y=\"0\" width=\"2400\" height=\"1600\"/>\n",
       "  </clipPath>\n",
       "</defs>\n",
       "<path clip-path=\"url(#clip2400)\" d=\"\n",
       "M0 1600 L2400 1600 L2400 0 L0 0  Z\n",
       "  \" fill=\"#ffffff\" fill-rule=\"evenodd\" fill-opacity=\"1\"/>\n",
       "<defs>\n",
       "  <clipPath id=\"clip2401\">\n",
       "    <rect x=\"480\" y=\"0\" width=\"1681\" height=\"1600\"/>\n",
       "  </clipPath>\n",
       "</defs>\n",
       "<path clip-path=\"url(#clip2400)\" d=\"\n",
       "M287.95 1487.47 L2352.76 1487.47 L2352.76 47.2441 L287.95 47.2441  Z\n",
       "  \" fill=\"#ffffff\" fill-rule=\"evenodd\" fill-opacity=\"1\"/>\n",
       "<defs>\n",
       "  <clipPath id=\"clip2402\">\n",
       "    <rect x=\"287\" y=\"47\" width=\"2066\" height=\"1441\"/>\n",
       "  </clipPath>\n",
       "</defs>\n",
       "<polyline clip-path=\"url(#clip2402)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  346.388,1487.47 346.388,47.2441 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip2402)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  833.371,1487.47 833.371,47.2441 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip2402)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  1320.35,1487.47 1320.35,47.2441 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip2402)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  1807.34,1487.47 1807.34,47.2441 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip2402)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  2294.32,1487.47 2294.32,47.2441 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip2402)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  287.95,1447.39 2352.76,1447.39 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip2402)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  287.95,1204.45 2352.76,1204.45 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip2402)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  287.95,961.499 2352.76,961.499 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip2402)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  287.95,718.551 2352.76,718.551 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip2402)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  287.95,475.604 2352.76,475.604 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip2402)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  287.95,232.656 2352.76,232.656 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip2400)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  287.95,1487.47 2352.76,1487.47 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip2400)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  287.95,1487.47 287.95,47.2441 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip2400)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  346.388,1487.47 346.388,1470.19 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip2400)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  833.371,1487.47 833.371,1470.19 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip2400)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  1320.35,1487.47 1320.35,1470.19 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip2400)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  1807.34,1487.47 1807.34,1470.19 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip2400)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  2294.32,1487.47 2294.32,1470.19 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip2400)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  287.95,1447.39 312.728,1447.39 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip2400)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  287.95,1204.45 312.728,1204.45 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip2400)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  287.95,961.499 312.728,961.499 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip2400)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  287.95,718.551 312.728,718.551 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip2400)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  287.95,475.604 312.728,475.604 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip2400)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  287.95,232.656 312.728,232.656 \n",
       "  \"/>\n",
       "<g clip-path=\"url(#clip2400)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 346.388, 1541.47)\" x=\"346.388\" y=\"1541.47\">0.00</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip2400)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 833.371, 1541.47)\" x=\"833.371\" y=\"1541.47\">0.25</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip2400)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 1320.35, 1541.47)\" x=\"1320.35\" y=\"1541.47\">0.50</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip2400)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 1807.34, 1541.47)\" x=\"1807.34\" y=\"1541.47\">0.75</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip2400)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 2294.32, 1541.47)\" x=\"2294.32\" y=\"1541.47\">1.00</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip2400)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:end;\" transform=\"rotate(0, 263.95, 1464.89)\" x=\"263.95\" y=\"1464.89\">0</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip2400)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:start;\" transform=\"rotate(0, 93.6653, 1228.17)\" x=\"93.6653\" y=\"1228.17\">5.0×10</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip2400)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:38px; text-anchor:start;\" transform=\"rotate(0, 242.206, 1200.76)\" x=\"242.206\" y=\"1200.76\">3</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip2400)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:start;\" transform=\"rotate(0, 93.6653, 985.226)\" x=\"93.6653\" y=\"985.226\">1.0×10</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip2400)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:38px; text-anchor:start;\" transform=\"rotate(0, 242.206, 957.816)\" x=\"242.206\" y=\"957.816\">4</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip2400)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:start;\" transform=\"rotate(0, 93.6653, 742.279)\" x=\"93.6653\" y=\"742.279\">1.5×10</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip2400)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:38px; text-anchor:start;\" transform=\"rotate(0, 242.206, 714.868)\" x=\"242.206\" y=\"714.868\">4</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip2400)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:start;\" transform=\"rotate(0, 93.6653, 499.331)\" x=\"93.6653\" y=\"499.331\">2.0×10</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip2400)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:38px; text-anchor:start;\" transform=\"rotate(0, 242.206, 471.921)\" x=\"242.206\" y=\"471.921\">4</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip2400)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:start;\" transform=\"rotate(0, 93.6653, 256.384)\" x=\"93.6653\" y=\"256.384\">2.5×10</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip2400)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:38px; text-anchor:start;\" transform=\"rotate(0, 242.206, 228.973)\" x=\"242.206\" y=\"228.973\">4</text>\n",
       "</g>\n",
       "<circle clip-path=\"url(#clip2402)\" style=\"fill:#000000; stroke:none; fill-opacity:1\" cx=\"346.388\" cy=\"1446.71\" r=\"18\"/>\n",
       "<circle clip-path=\"url(#clip2402)\" style=\"fill:#009af9; stroke:none; fill-opacity:1\" cx=\"346.388\" cy=\"1446.71\" r=\"14\"/>\n",
       "<circle clip-path=\"url(#clip2402)\" style=\"fill:#000000; stroke:none; fill-opacity:1\" cx=\"427.552\" cy=\"1443.7\" r=\"18\"/>\n",
       "<circle clip-path=\"url(#clip2402)\" style=\"fill:#009af9; stroke:none; fill-opacity:1\" cx=\"427.552\" cy=\"1443.7\" r=\"14\"/>\n",
       "<circle clip-path=\"url(#clip2402)\" style=\"fill:#000000; stroke:none; fill-opacity:1\" cx=\"508.716\" cy=\"1441.13\" r=\"18\"/>\n",
       "<circle clip-path=\"url(#clip2402)\" style=\"fill:#009af9; stroke:none; fill-opacity:1\" cx=\"508.716\" cy=\"1441.13\" r=\"14\"/>\n",
       "<circle clip-path=\"url(#clip2402)\" style=\"fill:#000000; stroke:none; fill-opacity:1\" cx=\"589.879\" cy=\"1436.41\" r=\"18\"/>\n",
       "<circle clip-path=\"url(#clip2402)\" style=\"fill:#009af9; stroke:none; fill-opacity:1\" cx=\"589.879\" cy=\"1436.41\" r=\"14\"/>\n",
       "<circle clip-path=\"url(#clip2402)\" style=\"fill:#000000; stroke:none; fill-opacity:1\" cx=\"671.043\" cy=\"1431.89\" r=\"18\"/>\n",
       "<circle clip-path=\"url(#clip2402)\" style=\"fill:#009af9; stroke:none; fill-opacity:1\" cx=\"671.043\" cy=\"1431.89\" r=\"14\"/>\n",
       "<circle clip-path=\"url(#clip2402)\" style=\"fill:#000000; stroke:none; fill-opacity:1\" cx=\"752.207\" cy=\"1428.1\" r=\"18\"/>\n",
       "<circle clip-path=\"url(#clip2402)\" style=\"fill:#009af9; stroke:none; fill-opacity:1\" cx=\"752.207\" cy=\"1428.1\" r=\"14\"/>\n",
       "<circle clip-path=\"url(#clip2402)\" style=\"fill:#000000; stroke:none; fill-opacity:1\" cx=\"833.371\" cy=\"1415.96\" r=\"18\"/>\n",
       "<circle clip-path=\"url(#clip2402)\" style=\"fill:#009af9; stroke:none; fill-opacity:1\" cx=\"833.371\" cy=\"1415.96\" r=\"14\"/>\n",
       "<circle clip-path=\"url(#clip2402)\" style=\"fill:#000000; stroke:none; fill-opacity:1\" cx=\"914.534\" cy=\"1404.39\" r=\"18\"/>\n",
       "<circle clip-path=\"url(#clip2402)\" style=\"fill:#009af9; stroke:none; fill-opacity:1\" cx=\"914.534\" cy=\"1404.39\" r=\"14\"/>\n",
       "<circle clip-path=\"url(#clip2402)\" style=\"fill:#000000; stroke:none; fill-opacity:1\" cx=\"995.698\" cy=\"1392.73\" r=\"18\"/>\n",
       "<circle clip-path=\"url(#clip2402)\" style=\"fill:#009af9; stroke:none; fill-opacity:1\" cx=\"995.698\" cy=\"1392.73\" r=\"14\"/>\n",
       "<circle clip-path=\"url(#clip2402)\" style=\"fill:#000000; stroke:none; fill-opacity:1\" cx=\"1076.86\" cy=\"1365.47\" r=\"18\"/>\n",
       "<circle clip-path=\"url(#clip2402)\" style=\"fill:#009af9; stroke:none; fill-opacity:1\" cx=\"1076.86\" cy=\"1365.47\" r=\"14\"/>\n",
       "<circle clip-path=\"url(#clip2402)\" style=\"fill:#000000; stroke:none; fill-opacity:1\" cx=\"1158.03\" cy=\"1358.38\" r=\"18\"/>\n",
       "<circle clip-path=\"url(#clip2402)\" style=\"fill:#009af9; stroke:none; fill-opacity:1\" cx=\"1158.03\" cy=\"1358.38\" r=\"14\"/>\n",
       "<circle clip-path=\"url(#clip2402)\" style=\"fill:#000000; stroke:none; fill-opacity:1\" cx=\"1239.19\" cy=\"1325.97\" r=\"18\"/>\n",
       "<circle clip-path=\"url(#clip2402)\" style=\"fill:#009af9; stroke:none; fill-opacity:1\" cx=\"1239.19\" cy=\"1325.97\" r=\"14\"/>\n",
       "<circle clip-path=\"url(#clip2402)\" style=\"fill:#000000; stroke:none; fill-opacity:1\" cx=\"1320.35\" cy=\"1297.45\" r=\"18\"/>\n",
       "<circle clip-path=\"url(#clip2402)\" style=\"fill:#009af9; stroke:none; fill-opacity:1\" cx=\"1320.35\" cy=\"1297.45\" r=\"14\"/>\n",
       "<circle clip-path=\"url(#clip2402)\" style=\"fill:#000000; stroke:none; fill-opacity:1\" cx=\"1401.52\" cy=\"1260.08\" r=\"18\"/>\n",
       "<circle clip-path=\"url(#clip2402)\" style=\"fill:#009af9; stroke:none; fill-opacity:1\" cx=\"1401.52\" cy=\"1260.08\" r=\"14\"/>\n",
       "<circle clip-path=\"url(#clip2402)\" style=\"fill:#000000; stroke:none; fill-opacity:1\" cx=\"1482.68\" cy=\"1222.28\" r=\"18\"/>\n",
       "<circle clip-path=\"url(#clip2402)\" style=\"fill:#009af9; stroke:none; fill-opacity:1\" cx=\"1482.68\" cy=\"1222.28\" r=\"14\"/>\n",
       "<circle clip-path=\"url(#clip2402)\" style=\"fill:#000000; stroke:none; fill-opacity:1\" cx=\"1563.84\" cy=\"1161.69\" r=\"18\"/>\n",
       "<circle clip-path=\"url(#clip2402)\" style=\"fill:#009af9; stroke:none; fill-opacity:1\" cx=\"1563.84\" cy=\"1161.69\" r=\"14\"/>\n",
       "<circle clip-path=\"url(#clip2402)\" style=\"fill:#000000; stroke:none; fill-opacity:1\" cx=\"1645.01\" cy=\"1089.19\" r=\"18\"/>\n",
       "<circle clip-path=\"url(#clip2402)\" style=\"fill:#009af9; stroke:none; fill-opacity:1\" cx=\"1645.01\" cy=\"1089.19\" r=\"14\"/>\n",
       "<circle clip-path=\"url(#clip2402)\" style=\"fill:#000000; stroke:none; fill-opacity:1\" cx=\"1726.17\" cy=\"1001.88\" r=\"18\"/>\n",
       "<circle clip-path=\"url(#clip2402)\" style=\"fill:#009af9; stroke:none; fill-opacity:1\" cx=\"1726.17\" cy=\"1001.88\" r=\"14\"/>\n",
       "<circle clip-path=\"url(#clip2402)\" style=\"fill:#000000; stroke:none; fill-opacity:1\" cx=\"1807.34\" cy=\"954.405\" r=\"18\"/>\n",
       "<circle clip-path=\"url(#clip2402)\" style=\"fill:#009af9; stroke:none; fill-opacity:1\" cx=\"1807.34\" cy=\"954.405\" r=\"14\"/>\n",
       "<circle clip-path=\"url(#clip2402)\" style=\"fill:#000000; stroke:none; fill-opacity:1\" cx=\"1888.5\" cy=\"842.017\" r=\"18\"/>\n",
       "<circle clip-path=\"url(#clip2402)\" style=\"fill:#009af9; stroke:none; fill-opacity:1\" cx=\"1888.5\" cy=\"842.017\" r=\"14\"/>\n",
       "<circle clip-path=\"url(#clip2402)\" style=\"fill:#000000; stroke:none; fill-opacity:1\" cx=\"1969.66\" cy=\"713.206\" r=\"18\"/>\n",
       "<circle clip-path=\"url(#clip2402)\" style=\"fill:#009af9; stroke:none; fill-opacity:1\" cx=\"1969.66\" cy=\"713.206\" r=\"14\"/>\n",
       "<circle clip-path=\"url(#clip2402)\" style=\"fill:#000000; stroke:none; fill-opacity:1\" cx=\"2050.83\" cy=\"589.449\" r=\"18\"/>\n",
       "<circle clip-path=\"url(#clip2402)\" style=\"fill:#009af9; stroke:none; fill-opacity:1\" cx=\"2050.83\" cy=\"589.449\" r=\"14\"/>\n",
       "<circle clip-path=\"url(#clip2402)\" style=\"fill:#000000; stroke:none; fill-opacity:1\" cx=\"2131.99\" cy=\"585.076\" r=\"18\"/>\n",
       "<circle clip-path=\"url(#clip2402)\" style=\"fill:#009af9; stroke:none; fill-opacity:1\" cx=\"2131.99\" cy=\"585.076\" r=\"14\"/>\n",
       "<circle clip-path=\"url(#clip2402)\" style=\"fill:#000000; stroke:none; fill-opacity:1\" cx=\"2213.15\" cy=\"282.363\" r=\"18\"/>\n",
       "<circle clip-path=\"url(#clip2402)\" style=\"fill:#009af9; stroke:none; fill-opacity:1\" cx=\"2213.15\" cy=\"282.363\" r=\"14\"/>\n",
       "<circle clip-path=\"url(#clip2402)\" style=\"fill:#000000; stroke:none; fill-opacity:1\" cx=\"2294.32\" cy=\"88.0053\" r=\"18\"/>\n",
       "<circle clip-path=\"url(#clip2402)\" style=\"fill:#009af9; stroke:none; fill-opacity:1\" cx=\"2294.32\" cy=\"88.0053\" r=\"14\"/>\n",
       "<polyline clip-path=\"url(#clip2402)\" style=\"stroke:#ff0000; stroke-width:8; stroke-opacity:1; fill:none\" points=\"\n",
       "  346.388,1434.7 366.064,1434.03 385.74,1433.33 405.416,1432.6 425.092,1431.83 444.769,1431.02 464.445,1430.16 484.121,1429.26 503.797,1428.32 523.473,1427.32 \n",
       "  543.149,1426.28 562.825,1425.17 582.501,1424.02 602.177,1422.8 621.853,1421.52 641.529,1420.17 661.205,1418.75 680.881,1417.27 700.557,1415.7 720.233,1414.05 \n",
       "  739.909,1412.32 759.585,1410.5 779.262,1408.58 798.938,1406.57 818.614,1404.45 838.29,1402.23 857.966,1399.89 877.642,1397.43 897.318,1394.85 916.994,1392.13 \n",
       "  936.67,1389.28 956.346,1386.28 976.022,1383.13 995.698,1379.82 1015.37,1376.34 1035.05,1372.69 1054.73,1368.85 1074.4,1364.82 1094.08,1360.58 1113.75,1356.14 \n",
       "  1133.43,1351.47 1153.11,1346.57 1172.78,1341.43 1192.46,1336.04 1212.13,1330.37 1231.81,1324.43 1251.49,1318.2 1271.16,1311.66 1290.84,1304.8 1310.52,1297.61 \n",
       "  1330.19,1290.08 1349.87,1282.18 1369.54,1273.9 1389.22,1265.23 1408.9,1256.14 1428.57,1246.63 1448.25,1236.67 1467.92,1226.25 1487.6,1215.34 1507.28,1203.93 \n",
       "  1526.95,1192 1546.63,1179.53 1566.3,1166.49 1585.98,1152.87 1605.66,1138.65 1625.33,1123.8 1645.01,1108.3 1664.68,1092.13 1684.36,1075.27 1704.04,1057.7 \n",
       "  1723.71,1039.4 1743.39,1020.34 1763.06,1000.5 1782.74,979.862 1802.42,958.406 1822.09,936.111 1841.77,912.958 1861.44,888.927 1881.12,864.002 1900.8,838.165 \n",
       "  1920.47,811.402 1940.15,783.698 1959.83,755.041 1979.5,725.421 1999.18,694.828 2018.85,663.256 2038.53,630.7 2058.21,597.157 2077.88,562.629 2097.56,527.119 \n",
       "  2117.23,490.631 2136.91,453.174 2156.59,414.761 2176.26,375.406 2195.94,335.127 2215.61,293.946 2235.29,251.888 2254.97,208.981 2274.64,165.257 2294.32,120.75 \n",
       "  \n",
       "  \"/>\n",
       "</svg>\n"
      ]
     },
     "execution_count": 411,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "scatter(x, y, leg=false)\n",
    "\n",
    "xg = range(0, 1.0, length=100)\n",
    "hx = [h(β, xi) for xi in xg]\n",
    "\n",
    "plot!(xg, hx, c=:red, lw=2)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 1.3.0",
   "language": "julia",
   "name": "julia-1.3"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.3.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}