diff --git a/10-12/lesson.ipynb b/10-12/lesson.ipynb
index 232ba27..d247648 100644
--- a/10-12/lesson.ipynb
+++ b/10-12/lesson.ipynb
@@ -1282,7 +1282,7 @@
"include(\"DIS.jl\")\n",
"include(\"UNM.jl\")\n",
"\n",
- "plt = plot(yscale = :log,\n",
+ "plt = plot(yscale = :log10,\n",
" xlims=(0, 35),\n",
" ylims=(1e-15, Inf),\n",
" guidefontsize=16)\n",
@@ -1299,7 +1299,7 @@
],
"metadata": {
"kernelspec": {
- "display_name": "Julia 1.9.3",
+ "display_name": "Julia 1.9.4",
"language": "julia",
"name": "julia-1.9"
},
@@ -1307,7 +1307,7 @@
"file_extension": ".jl",
"mimetype": "application/julia",
"name": "julia",
- "version": "1.9.3"
+ "version": "1.9.4"
}
},
"nbformat": 4,
diff --git a/11-09/NWTN.jl b/11-09/NWTN.jl
index d7b364e..641376d 100644
--- a/11-09/NWTN.jl
+++ b/11-09/NWTN.jl
@@ -12,7 +12,7 @@ function NWTN(f;
MInf::Real=-Inf,
mina::Real=1e-16,
plt::Union{Plots.Plot, Nothing}=nothing,
- plotatend::Bool=true,
+ plotatend::Bool=false,
Plotf::Integer=0,
printing::Bool=true)::Tuple{AbstractArray, String}
@@ -99,7 +99,7 @@ function NWTN(f;
am = 0
a = as
phipm = phip0
- while (feval ≤ MaxFeval ) && ((as - am)) > mina && (phips > 1e-12)
+ while (feval ≤ MaxFeval ) && (as - am) > mina && (phips > 1e-12)
# compute the new value by safeguarded quadratic interpolation
a = (am * phips - as * phipm) / (phips - phipm)
@@ -107,7 +107,7 @@ function NWTN(f;
# compute phi(a)
phia, phip = f2phi(a, true)
-
+
if (phia ≤ phi0 + m1 * a * phip0) && (abs(phip) ≤ -m2 * phip0)
break # Armijo + strong Wolfe satisfied, we are done
end
@@ -372,12 +372,12 @@ function NWTN(f;
end
end
+ if Plotf ≥ 2
+ plot!(plt, gap)
+ elseif Plotf == 1 && n == 2
+ plot!(plt, PXY[1, :], PXY[2, :])
+ end
if plotatend
- if Plotf ≥ 2
- plot!(plt, gap)
- elseif Plotf == 1 && n == 2
- plot!(plt, PXY[1, :], PXY[2, :])
- end
display(plt)
end
diff --git a/11-09/Untitled.ipynb b/11-09/Untitled.ipynb
deleted file mode 100644
index c58f437..0000000
--- a/11-09/Untitled.ipynb
+++ /dev/null
@@ -1,570 +0,0 @@
-{
- "cells": [
- {
- "cell_type": "code",
- "execution_count": 27,
- "id": "cb5170b7-6a4c-415f-8675-8382c5ce837c",
- "metadata": {},
- "outputs": [
- {
- "data": {
- "text/plain": [
- "NWTN (generic function with 1 method)"
- ]
- },
- "execution_count": 27,
- "metadata": {},
- "output_type": "execute_result"
- }
- ],
- "source": [
- "using LinearAlgebra, Printf, Plots\n",
- "\n",
- "function NWTN(f;\n",
- " x::Union{Nothing, Vector}=nothing,\n",
- " eps::Real=1e-6,\n",
- " MaxFeval::Integer=1000,\n",
- " m1::Real=1e-4,\n",
- " m2::Real=0.9,\n",
- " delta::Real=1e-6,\n",
- " tau::Real=0.9,\n",
- " sfgrd::Real=0.2,\n",
- " MInf::Real=-Inf,\n",
- " mina::Real=1e-16,\n",
- " plt::Union{Plots.Plot, Nothing}=nothing,\n",
- " plotatend::Bool=true,\n",
- " Plotf::Integer=0,\n",
- " printing::Bool=true)::Tuple{AbstractArray, String}\n",
- "\n",
- " \n",
- " # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -\n",
- " # inner functions - - - - - - - - - - - - - - - - - - - - - - - - - - - - -\n",
- " # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -\n",
- "\n",
- " function f2phi(alpha, derivate=false)\n",
- " # computes and returns the value of the tomography at alpha\n",
- " #\n",
- " # phi( alpha ) = f( x + alpha * d )\n",
- " #\n",
- " # if Plotf > 2 saves the data in gap() for plotting\n",
- " #\n",
- " # if the second output parameter is required, put there the derivative\n",
- " # of the tomography in alpha\n",
- " #\n",
- " # phi'( alpha ) = < \\nabla f( x + alpha * d ) , d >\n",
- " #\n",
- " # saves the point in lastx, the gradient in lastg, the Hessian in lasth,\n",
- " # and increases feval\n",
- " \n",
- " lastx = x + alpha * d\n",
- " phi, lastg, lastH = f(lastx)\n",
- " \n",
- " if Plotf > 2\n",
- " if fStar > - Inf\n",
- " push!(gap, (phi - fStar) / max(abs(fStar), 1))\n",
- " else\n",
- " push!(gap, phi)\n",
- " end\n",
- " end\n",
- "\n",
- " feval += 1\n",
- "\n",
- " if derivate\n",
- " return (phi, dot(d, lastg))\n",
- " end\n",
- " return (phi, nothing)\n",
- " end\n",
- "\n",
- " # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -\n",
- " # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -\n",
- " \n",
- " function ArmijoWolfeLS(phi0, phip0, as, m1, m2, tau)\n",
- " # performs an Armijo-Wolfe Line Search.\n",
- " #\n",
- " # phi0 = phi( 0 ), phip0 = phi'( 0 ) < 0\n",
- " #\n",
- " # as > 0 is the first value to be tested: if phi'( as ) < 0 then as is\n",
- " # divided by tau < 1 (hence it is increased) until this does not happen\n",
- " # any longer\n",
- " #\n",
- " # m1 and m2 are the standard Armijo-Wolfe parameters; note that the strong\n",
- " # Wolfe condition is used\n",
- " #\n",
- " # returns the optimal step and the optimal f-value\n",
- " \n",
- " lsiter = 1 # count iterations of first phase\n",
- " local phips, phia\n",
- " while feval ≤ MaxFeval \n",
- " phia, phips = f2phi(as, true)\n",
- " \n",
- " if (phia ≤ phi0 + m1 * as * phip0) && (abs(phips) ≤ - m2 * phip0)\n",
- " if printing\n",
- " @printf(\" %2d\", lsiter)\n",
- " end\n",
- " a = as;\n",
- " return (a, phia) # Armijo + strong Wolfe satisfied, we are done\n",
- " \n",
- " end\n",
- " if phips ≥ 0\n",
- " break\n",
- " end\n",
- " as = as / tau\n",
- " lsiter += 1\n",
- " end \n",
- "\n",
- " if printing\n",
- " @printf(\" %2d \", lsiter)\n",
- " end\n",
- " lsiter = 1 # count iterations of second phase\n",
- " \n",
- " am = 0\n",
- " a = as\n",
- " phipm = phip0\n",
- " while (feval ≤ MaxFeval ) && ((as - am)) > mina && (phips > 1e-12)\n",
- " \n",
- " # compute the new value by safeguarded quadratic interpolation\n",
- " a = (am * phips - as * phipm) / (phips - phipm)\n",
- " a = max(am + ( as - am ) * sfgrd, min(as - ( as - am ) * sfgrd, a))\n",
- " \n",
- " # compute phi(a)\n",
- " phia, phip = f2phi(a, true)\n",
- " \n",
- " if (phia ≤ phi0 + m1 * a * phip0) && (abs(phip) ≤ -m2 * phip0)\n",
- " break # Armijo + strong Wolfe satisfied, we are done\n",
- " end\n",
- " \n",
- " # restrict the interval based on sign of the derivative in a\n",
- " if phip < 0\n",
- " am = a\n",
- " phipm = phip\n",
- " else\n",
- " as = a\n",
- " if as ≤ mina\n",
- " break\n",
- " end\n",
- " phips = phip\n",
- " end\n",
- " lsiter += 1\n",
- " end\n",
- "\n",
- " if printing\n",
- " @printf(\"%2d\", lsiter)\n",
- " end\n",
- " return (a, phia)\n",
- " end\n",
- "\n",
- " # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -\n",
- "\n",
- " function BacktrackingLS(phi0, phip0, as, m1, tau)\n",
- " # performs a Backtracking Line Search.\n",
- " #\n",
- " # phi0 = phi( 0 ), phip0 = phi'( 0 ) < 0\n",
- " #\n",
- " # as > 0 is the first value to be tested, which is decreased by\n",
- " # multiplying it by tau < 1 until the Armijo condition with parameter\n",
- " # m1 is satisfied\n",
- " #\n",
- " # returns the optimal step and the optimal f-value\n",
- " \n",
- " lsiter = 1 # count ls iterations\n",
- " while feval ≤ MaxFeval && as > mina\n",
- " phia, _ = f2phi(as)\n",
- " if phia ≤ phi0 + m1 * as * phip0 # Armijo satisfied\n",
- " break # we are done\n",
- " end\n",
- " as *= tau\n",
- " lsiter += 1\n",
- " end\n",
- "\n",
- " if printing\n",
- " @printf(\" %2d\", lsiter)\n",
- " end\n",
- "\n",
- " return (as, phia)\n",
- " end\n",
- " \n",
- " # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -\n",
- " # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -\n",
- "\n",
- " #Plotf = 2\n",
- " # 0 = nothing is plotted\n",
- " # 1 = the level sets of f and the trajectory are plotted (when n = 2)\n",
- " # 2 = the function value / gap are plotted, iteration-wise\n",
- " # 3 = the function value / gap are plotted, function-evaluation-wise\n",
- "\n",
- " Interactive = false\n",
- "\n",
- " PXY = Matrix{Real}(undef, 2, 0)\n",
- "\n",
- " local gap\n",
- " if Plotf > 1\n",
- " if Plotf == 2\n",
- " MaxIter = 50 # expected number of iterations for the gap plot\n",
- " else\n",
- " MaxIter = 70 # expected number of iterations for the gap plot\n",
- " end\n",
- " gap = []\n",
- " end\n",
- "\n",
- " if Plotf == 2 && plt == nothing\n",
- " plt = plot(xlims=(0, MaxIter), ylims=(1e-15, 1e+1))\n",
- " end\n",
- " if Plotf > 1 && plt == nothing\n",
- " plt = plot(xlims=(0, MaxIter))\n",
- " end\n",
- " if plt == nothing\n",
- " plt = plot()\n",
- " end\n",
- "\n",
- " local fStar\n",
- " if isnothing(x)\n",
- " (fStar, x, _) = f(nothing)\n",
- " else\n",
- " (fStar, _, _) = f(nothing)\n",
- " end\n",
- "\n",
- " n = size(x, 1)\n",
- "\n",
- " if m1 ≤ 0 || m1 ≥ 1\n",
- " throw(ArgumentError(\"m1 is not in (0 ,1)\"))\n",
- " end\n",
- "\n",
- " AWLS = (m2 > 0 && m2 < 1)\n",
- "\n",
- " if delta < 0\n",
- " throw(ArgumentError(\"delta must be ≥ 0\"))\n",
- " end\n",
- "\n",
- " if tau ≤ 0 || tau ≥ 1\n",
- " throw(ArgumentError(\"tau is not in (0 ,1)\"))\n",
- " end\n",
- "\n",
- " if sfgrd ≤ 0 || sfgrd ≥ 1\n",
- " throw(ArgumentError(\"sfgrd is not in (0, 1)\"))\n",
- " end\n",
- "\n",
- " if mina < 0\n",
- " throw(ArgumentError(\"mina must be ≥ 0\"))\n",
- " end\n",
- "\n",
- " # \"global\" variables- - - - - - - - - - - - - - - - - - - - - - - - - - - -\n",
- " # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -\n",
- " \n",
- " lastx = zeros(n) # last point visited in the line search\n",
- " lastg = zeros(n) # gradient of lastx\n",
- " lastH = zeros(n, n) # Hessian of lastx\n",
- " d = zeros(n) # Newton's direction\n",
- " feval = 0 # f() evaluations count (\"common\" with LSs)\n",
- "\n",
- " status = \"error\"\n",
- " \n",
- " # initializations - - - - - - - - - - - - - - - - - - - - - - - - - - - - -\n",
- " # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -\n",
- "\n",
- " if fStar > -Inf\n",
- " prevv = Inf\n",
- " end\n",
- "\n",
- " if printing\n",
- " @printf(\"Newton's method\\n\")\n",
- " if fStar > -Inf\n",
- " @printf(\"feval\\trel gap\\t\\t|| g(x) ||\\trate\\t\\tdelta\")\n",
- " else\n",
- " @printf(\"feval\\tf(x)\\t\\t\\t|| g(x) ||\\tdelta\")\n",
- " end\n",
- " @printf(\"\\t\\tls it\\ta*\")\n",
- " @printf(\"\\n\\n\")\n",
- " end\n",
- "\n",
- " \n",
- " v, _ = f2phi(0)\n",
- " ng = norm(lastg)\n",
- " if eps < 0\n",
- " ng0 = -ng # norm of first subgradient: why is there a \"-\"? ;-)\n",
- " else\n",
- " ng0 = 1 # un-scaled stopping criterion\n",
- " end\n",
- "\n",
- " # main loop - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -\n",
- " # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -\n",
- "\n",
- " while true\n",
- "\n",
- " # output statistics - - - - - - - - - - - - - - - - - - - - - - - - - -\n",
- " \n",
- " if fStar > -Inf\n",
- " gapk = ( v - fStar ) / max(abs( fStar ), 1)\n",
- "\n",
- " if printing\n",
- " @printf(\"%4d\\t%1.4e\\t%1.4e\", feval, gapk, ng)\n",
- " if prevv < Inf\n",
- " @printf(\"\\t%1.4e\", ( v - fStar ) / ( prevv - fStar ))\n",
- " else\n",
- " @printf(\"\\t\\t\")\n",
- " end\n",
- " end\n",
- " prevv = v\n",
- " \n",
- " if Plotf > 1\n",
- " if Plotf ≥ 2\n",
- " push!(gap, gapk)\n",
- " end\n",
- " end\n",
- " else\n",
- " if printing\n",
- " @printf(\"%4d\\t%1.8e\\t\\t%1.4e\", feval, v, ng)\n",
- " end\n",
- " \n",
- " if Plotf > 1\n",
- " if Plotf ≥ 2\n",
- " push!(gap, v)\n",
- " end\n",
- " end\n",
- " end\n",
- "\n",
- " # stopping criteria - - - - - - - - - - - - - - - - - - - - - - - - - -\n",
- "\n",
- " if ng ≤ eps * ng0\n",
- " status = \"optimal\"\n",
- " break\n",
- " end\n",
- " \n",
- " if feval > MaxFeval\n",
- " status = \"stopped\"\n",
- " break\n",
- " end\n",
- " \n",
- " # compute Newton's direction- - - - - - - - - - - - - - - - - - - - - -\n",
- "\n",
- " lambdan = eigmin(lastH) # smallest eigenvalue\n",
- " if lambdan < delta\n",
- " if printing\n",
- " @printf(\"\\t%1.4e\", delta - lambdan)\n",
- " end\n",
- " lastH = lastH + (delta - lambdan) * I\n",
- " else\n",
- " if printing\n",
- " @printf(\"\\t0.00e+00\")\n",
- " end\n",
- " end\n",
- " d = -lastH \\ lastg\n",
- " phip0 = lastg' * d\n",
- " \n",
- " # compute step size - - - - - - - - - - - - - - - - - - - - - - - - - -\n",
- " # in Newton's method, the default initial stepsize is 1\n",
- " \n",
- " if AWLS\n",
- " a, v = ArmijoWolfeLS(v, phip0, 1, m1, m2, tau)\n",
- " else\n",
- " a, v = BacktrackingLS(v, phip0, 1, m1, tau)\n",
- " end\n",
- " \n",
- " # output statistics - - - - - - - - - - - - - - - - - - - - - - - - - -\n",
- "\n",
- " if printing\n",
- " @printf(\"\\t%1.2e\", a)\n",
- " @printf(\"\\n\")\n",
- " end\n",
- " \n",
- " if a ≤ mina\n",
- " status = \"error\"\n",
- " break\n",
- " end\n",
- " \n",
- " if v ≤ MInf\n",
- " status = \"unbounded\"\n",
- " break\n",
- " end\n",
- " \n",
- " # compute new point - - - - - - - - - - - - - - - - - - - - - - - - - -\n",
- " \n",
- " # possibly plot the trajectory\n",
- " if n == 2 && Plotf == 1\n",
- " PXY = hcat(PXY, hcat(x, lastx))\n",
- " end\n",
- " \n",
- " x = lastx\n",
- " ng = norm(lastg)\n",
- "\n",
- " # iterate - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -\n",
- " \n",
- " if Interactive\n",
- " readline()\n",
- " end\n",
- " end\n",
- "\n",
- " if plotatend\n",
- " if Plotf ≥ 2\n",
- " plot!(plt, gap)\n",
- " elseif Plotf == 1 && n == 2\n",
- " plot!(plt, PXY[1, :], PXY[2, :])\n",
- " end\n",
- " display(plt)\n",
- " end\n",
- "\n",
- " # end of main loop- - - - - - - - - - - - - - - - - - - - - - - - - - - - -\n",
- " # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -\n",
- " return (x, status)\n",
- "end"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": 10,
- "id": "9fa92a74-19bc-46b3-92de-d39aa74f0075",
- "metadata": {},
- "outputs": [],
- "source": [
- "include(\"./TestFunctions/TestFunctions.jl\")\n",
- "TF = TestFunctions();"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": 29,
- "id": "14e530e5-1b32-4a3e-ac62-b67fd8723f06",
- "metadata": {},
- "outputs": [
- {
- "data": {
- "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlgAAAGQCAIAAAD9V4nPAAAABmJLR0QA/wD/AP+gvaeTAAAgAElEQVR4nO3dd3wU1doH8GfO7KYXkkA6gQQSeglEauiEIgQuCFwpUlTg2t97FRQLoKKUiwUVC02UoggiUlV6EwHTKKGEGkgFQkIqycyc94/lRqQuYXdnd+b3/fjHZp1kn50s88sz58wZgXNOAAAAesXULgAAAEBNCEIAANA1BCEAAOgaghAAAHQNQQgAALqGIAQAAF1DEAIAgK4hCAEAQNcQhAAAoGsIQgAA0DU7CsJXX321oqLCzI1lWcbicCoy/zcF1oD9ryLOuSRJalehX5IkWfzgb0dBuGDBgsLCQjM3Li8vl2XZqvXAPZSVlaldgq5h/6tIURT8IaKi8vJyRVEs+zPtKAgBAABsD0EIAAC6hiAEAABdQxACAICuGczcTlGUkydPHjt2LCoqqlGjRnfcJjc3d/Xq1ZIkDRgwICQkpPL5HTt2HDx4MCoqKj4+njFELwAA2BFzYyk+Pr5z587jxo1buXLlHTfIzMxs1qzZH3/8kZqa2rRp07S0NNPzM2fOHD16dEFBwdSpU5966inLVA0AAGAhgpkXZBQUFHh7ew8bNiwqKmrq1Km3b/DGG2+cOnVqxYoVRPTCCy9UVFR8+eWXxcXFoaGhmzdvjomJuXLlSlhYWHJycmRk5B1fws/PLy0tzdfX15x6SktLjUajwWBuRwuWVVhY6OnpqXYV+oX9ryJZlsvLy11dXdUuRKdKSkqcnZ1FUbTgzzS3I/T29r73Br/++mt8fLzpcd++fX/99VciOnDggJubW0xMDBH5+fm1bdt28+bND1EtAACAhVmso8rMzAwKCjI9Dg4OzszM5Jzf/CQRBQUFZWZm3u0nlJWVTZ48ufLvrObNmw8aNOhuG+/MkII9eD0fXFOvjuvXrzs5OaldhX5h/6vI1BFiuoNaNl+Q6vlSuJe5HaHRaLzvL8tiv0tB+OssK+dcEIRbnrz5+bv9hGo3MRqN93i5pKts5lFLtsYAAGD/JiQ7Xyq7a45UjcU6wqCgoOzsbNPj7OzsoKAgQRBuftL0fLt27e72E5ydnf/zn/+YOUY4Moo3/pmVCsZq+LNYDeXl5c7OzmpXoV/Y/yqSZVkQBOx/Vfyew4kq2gY7GVQZI7yj0tLSEydOmB737Nlz/fr1psfr1q3r2bMnEbVq1aqkpCQhIYGI8vLy9u3bFxcX93AF3+DrxLsECSvPWHjFOQAAsFsLTyhjIrmF+0HzO8LFixdv3Lhx//79KSkpqampTz/9dI8ePZKSktq3b286+fncc8+1bNlyzJgxrq6uP/zww++//05E7u7ur7322mOPPTZixIgNGzYMHjz4blNGq2B0JE0/pIytjzP1AADaVyzRT+eVhD6W73/MDcLmzZu7u7sPHjzY9GWdOnWIqH79+mvWrDE9ExwcfOjQIdMF9SkpKZUX1L/66qtt2rQ5ePDgO++807dvXwuW3jOEnt9Hqfm8YTWL/30AAAD2ZcUZpXMQC3Cx/A34zL2O0AaqcB3h5CShQqH/tsasGVvDdWzqwv5XEa4jVEv7ddLrzcUufqWqXUdon8ZEsSWnlAoMFAIAaNqJAn62kPcMscr5P8cOwkhvIcpb2HABSQgAoGULjitjopjBOpHl2EFIRE9GsUUn7OXsLgAAWJyk0LLTyuhIawWWwwfh4Ai2J0fJLEEWAgBo07p0JcpbiPS21rxIhw9CdwM9VpstSUMQAgBo07dpfEyUFdPK4YOQiMZEsUUnFSQhAID25JfT9izlH7UQhPfULkAwCLQ3G1EIAKA1q88p3UOYtzVX09RCEBLR6Cj29UnMHQUA0JrvTitDI6y7aopGgnBkJFt9TimsULsOAACwnNxSSrjMH61p3ajSSBAGuFKnILbyLJpCAADtWHFGiQ9jrha7T9KdaSQIiWhMlLDoBIIQAEA7vj+jPB5h9ZzSThD2qcnOFtKxfEyZAQDQgvQinlbAu1tnWbWbaScIDYxG1BW+SUNTCACgBd+d5o+FM6P1Y0o7QUhET9Vj35zEGtwAAFrw3WllaB1bhJSmgjDKW4jwEn65iCQEAHBsx/P5pTKKDbDF7WY1FYSENbgBADRh+WllaB2B2eS261oLwiERbEeWklOqdh0AAFBVnOi70/xxm5wXJe0FoaeRBtRmS07h7CgAgKP67SL3MFJMdZv0g9oLQiIaE8UWHMca3AAAjmpuqvJCI9vFkwaDsEOgIAj0Ry6iEADA8aQX8X25triOvpIGg5CIRkViDW4AAIf0+TFlVCRzs/KyajfTZhCOjmKrzipFWIMbAMChXJdp8UllXH2bZpM2gzDQldoHCD+eQ1MIAOBIVpxRWlQXorxtNE3GRJtBSKbb1mMNbgAAhzI3VXmuoWjjF9VsEPYLY2nX+HGswQ0A4CCSrvCcUuoVatN2kDQchAZGw+rggkIAAIfxyVHl+UZMtHUOajcIiWhMFPsmjcvoCQEA7N7V67TuvDImSoVU0nIQNvIRQt3p14tIQgAAezf/hNKvFvNzVuGltRyEZJoygwsKAQDsm8Lpy2PKsw3ViSSNB+HQOmxrhnKpTO06AADg7jZe4P6utltc9BYaD0IvI/WrxZZiygwAgB2bmyo/p1I7SJoPQsIFhQAA9u30NZ54hQ8ORxBaTacgoUymg5cwZQYAwB59fkx5Moq52Poy+r9oPwgFotGYMgMAYJdKJVqSZuvFRW+h/SAkotGRwg9nlBJJ7ToAAODvlp9W2gawcE91psmY6CIIQ9yF1v7CaqzBDQBgZ744pqg4TcZEF0FImDIDAGB/9uXygnLqHqxmO0j6CcL+tdjRfH76GqbMAADYi7mpynMNGVM5B3UThE6MhtZh36ShKQQAsAuXymjjBWVkpPoxpH4FNvN0Pfb1SazBDQBgF+YfVwaFM181Fhe9hY6CsLGPEOhKWzKQhAAAKpM5zT+h/EvVqyYq2UURNoM1uAEA7MG6dCXEjVqotLjoLfQVhMPrst8uKpexBjcAgKpM02TUruIGe6nDNrydqE8YW34aTSEAgGrSCvjhPD6wtr0EkL3UYTNjotgCXFAIAKCez1KVcfWZs3qLi95Cd0HYNVgokSjxMqbMAACooKiClp1SeXHRW9hRKbYhED1RF1NmAADUsfSU0jmYhbrbxTQZE90FIRE9VU/47rRSijW4AQBszh4WF72FfVVjG6HuQkx1Yc15NIUAADa1K5tXKNQ5yI7aQdJnEBLRmCj2Nc6OAgDYlumqCfuKQd0G4cBwlnyFnynElBkAABvJKqEtGcoTdrC46C3sriDbcGL0zwj2LdbgBgCwlXnHlcfrMC+j2nXcRqdBSERj67PFJ7mCnhAAwPokhRbYzeKit7DHmmyjqa/g60zbMpGEAABW99N5pa4XNfG1t/FBIj0HIWENbgAAW7GrxUVvYadl2cbwumzTBeXqdbXrAADQtNR8nlZA/WvZaeLYaVm24etMvWpiDW4AAOv67KgyvgEz2mvg2GtdtvIkzo4CAFjTtQpacUYZW89+48Z+K7ONbsHClTJKvoIpMwAAVrHslNItmAW5qV3H3ek9CJlAoyIFrDIDAGAl848r4xvYddbYdXG2MSaKLT+tXJfVrgMAQHP25/KCcupiZ4uL3gJBSLU9hWa+ws9YgxsAwNLmHVfGN2B2t7ro3yEIibAGNwCAFRSU00/nldH2t7joLey9Ptt4LJwdvMTTizBlBgDAYpaeUuJCmL+r2nXcD4KQiMhFpCER7Js0BCEAgMUsOKGMt8vFRW/hACXaxpP12NcnFazBDQBgEX/k8sIK6hJs38ODRIQgrBRTXfAy0s5sJCEAgAV8dVwZX9/OZ8ncgCD8y+go9vUJTJkBAHhYBeX083lllN1PkzFxjCpt44lIti5dyS9Xuw4AAAe35JTSM9QBpsmYIAj/4udM3UPY91iDGwDg4Sw8oYxzhGkyJg5TqG08iQsKAQAezu85vLCCOtv3ajI3QxD+Tc9QIauEDuVhygwAQBXNO678q4FjTJMxMZi/qaIohw4dYow1adJEEG59j6WlpWVlZTc/U61aNUEQiouLy8tvDLsJglCtWrWHrNiqmEAjI4XFJ5UP24hq1wIA4Hjyy2ltujKrlVHtQh6AuUFYUFDQvXv38vJyWZarVav2yy+/eHh43LzBvHnzPvnkE9PjwsLCkpKSvLw8Jyenp59+esuWLV5eXkTk7e2dmJho2TdgcU/VY61/lqY/IjojCgEAHtC3aUovx5kmY2LuqdFPP/3Uz88vKSkpJSXFYDDMmzfvlg1eeuml0//z2GOPDRo0yMnJyfS/3n77bdPz9p+CRBTuKTTyEdanY6QQAOCBzT/uSNNkTMwtd+XKlaNGjWKMiaI4atSolStX3m3L0tLS77///qmnnqp8pqSk5Ny5c5UnSO0f1uAGAKiCvTm8QqFOjjNNxsTcIExPTw8PDzc9rl279oULF+625apVq2rUqBEbG1v5zOzZs7t37+7j4zN9+vR7vIQkSTt37tzyP8eOHTOzNosbHM725fILxZgyAwDwAOYdV8Y5yGoyNzN3jLC0tLTyVKeLi0txcfHdtly0aNGTTz5ZOZtmzpw5/v7+RJScnNypU6eYmJi4uLg7fuP169c//PBDo/HGEGtsbOzEiRPvUY/RaDQYHmCyzwPpH2pYeFR+paFkpZ/v6IqLi2+fMAU2g/2vIlmWTbMl1C7E7hRUCD+fM77duKLImnfyKSkpqaioEEVzJ3G4uLjcNynMDZLAwMC8vDzT4ytXrgQFBd1xs7Nnz+7du3fZsmWVz5hSkIiaN2/eq1evXbt23S0I3d3df/75Z19fX3PqEUXRqkE4vhEftl2e0soFB5s74pzfMlsKbAn7X0WmIHR1dajZIDax8IjSJ4zX9nOy6qswxpydnc0PQrN+ppnbPfLII3v27DE93r17d6tWre642cKFC3v16hUcHHzH/3vhwgU/P78qVGl7rf0FDyPtxhrcAADmWXBCGd/AwabJmJjbUb300kt9+/atU6eOJElffvnl1q1bTc8HBQX9+OOP7dq1IyJFUb799tvKiyiIqKKiYuzYsd27d/fw8Fi9evXp06eHDh1q8fdgJSMj2dcnlY6BuIoCAOA+9mTzCoU6BDrkSTRzgzA2NnblypWLFy9mjK1duzYmJsb0/OjRowMDA02Ps7KyRowY0adPn8rvEkWxQYMGmzZtqqioqFevXkpKSkBAgGXfgPWMimRRKyuuVYhejnRhKACACuYdV8Y71GoyNxM4t5ezf35+fmlpaWaOEVp7sozJwC1y71BhrKNdE2MDhYWFnp6ealehX9j/KsIY4e3yyyliRcWJwcYaLlZ/rZKSEtXGCPUJa3ADANzX4pNKn5rMBiloJQjCe+kVKpwvotR8e2maAQDs0AKHuunS7Ry4dBswMBoZKeC29QAAd7M7m8ucYh1zmowJgvA+xkSxJaeUCkQhAMCdfHVcGe+Aq8ncDEF4H1HeQqS3sOECkhAA4FZXrtOGdGVkpGNHiWNXbxtPRbFFJzBMCABwq2/TlL5hzNdZ7ToeDoLw/gZHsD05SmYJshAA4C8VCs05orzQyOFzxOHfgA24G2hEXfbpUZwdBQD4y7dpSj1valXDoccHiRCEZprQlM0/rhQ4zB0VAQCsS+Y085DyZrQWFqFEEJqlprvQM5TNx3UUAABERLT8tBLi5qiLi94CQWiuV5uxj48o5YhCANA9hdPMFOUtTbSDhCA0X1NfoYkPLT2FJAQAvVt5VvE0UtdgLbSDhCB8IBObibNSFAWzRwFAxzjR+8nKlBYaaQcJQfhAugQJPs60Lh1NIQDo15pzioFRz1CNtIOEIHxQrzRh7ycjCAFAv95PVqZEO/aaardAED6YAbXZ1XLanY3TowCgRxsu8AqF4mtpKjs09WZsgAn0ShM265CsdiEAACqYniy/qa12kBCEVTAqkiVepiNX0RQCgL5szuB512lgba0Fh9bejw04i/RCI/bfQxgpBAB9mZYkvxnNtNYPIgir5tmGbOMFJb0ITSEA6MXvOTyjhIaEazA1NPiWbMDLSKMj2cdH0BQCgF5MSZTfbM4MWgwNLb4nm/i/xmxxmnLlutp1AABY3/5cfrKAhtXVZmRo813ZQIi7MKAW+yIVTSEAaN87SfIbzZmTRhNDo2/LJiY2Y3NT5VJJ7ToAAKwp6Qo/nEejIjWbF5p9YzZQz1to7c8Wp6EpBAAteydRmdiMOWtnbdFbIQgfyuvN2X8PKRKiEAA06uhVvv+S8lSUlsNCy+/NBlrVEELd6cdzSEIA0KZ3kpSXm4iuBrXrsCYE4cOa2FScnoxbMwGABh3P5zuylPH1NZ4UGn97NtAnTJA5bc1AFAKA1ryXrPxfY9HDqHYdVoYgfFgC0YSmbCaW4QYAbTl9jW+6oDzTQPsxof13aAND67C0Akq4jKYQALTj/WTlhUZiNSe167A+BKEFGBm91BjLcAOAdqQX8bXpyouNdJERuniTNjCuPtuepZy6hqYQALRgeooyrj7zcVa7DptAEFqGu4HG1WcfHkZTCAAOL6uEVp5R/q+xdi+h/zsEocW80FD8/oySXap2HQAAD2dGijwmitVwUbsOW0EQWoy/Kz0eweamYvooADiwnFJaekr5dxMdpYOO3qoNTGjKvjqmFFaoXQcAQFXNPiSPqMuC3TR3H/q7QxBaUrin0DWYLTiBkUIAcEhXrtOik8rLemoHCUFocZOasw8PK+WIQgBwQB8elodEsDAPHbWDhCC0uGa+Qv1q9P1pJCEAOJiCcpp/XJnYVHe5oLs3bAOvNhVnpChYhxsAHMvHR5T4MBbuqa92kBCE1tA9RHAz0MYLSEIAcBhFFfTFMfnVZnoMBT2+ZxuY0JTNwjLcAOA4PktV4kJYlLfu2kFCEFrJoHCWVUK/56ApBAAHUCLRx0d02g4SgtBKRIH+04TNwjLcAOAIvjymdApijX302A4SgtB6xkSxg5f40atoCgHArpXJ9OER5Y3m+o0D/b5za3MR6ZkGWIYbAOzdghPKI9WFpr46bQcJQWhVzzVka84rF4rRFAKAnapQaPYh5XUdt4OEILQqH2caFck+OYKmEADs1NcnlUY+9EgN/baDhCC0tn83YV+fVPLL1a4DAOA2kkIzUpQ3muvlvoN3gyC0rpruwj9qsw8O45pCALA7S08pEZ7ULkDX7SAhCG3grWj2RaqSixv2AoA9kTnNSFHeitZ7O0gIQhuo5SE8Xof9FwvNAIA9WXFG8XOhTkF6bwcJQWgbbzQXF51UMjB9FADsAyeanqxMbYF2kAhBaBtBbjQmis1IwfRRALALP55V3I0UF4J2kAhBaDOvNRO/O62cLURTCAAq44TRwb9BENpIdRf6VwP2fjKaQgBQ2brziszp0ZpoB29AENrOK03Fn88rJwrQFAKAmqYlK5OjGWKwEoLQdqo50YuNxHeT0BQCgGp+vchLJfpHbRz8/4J9YVMvNWZbMhTckgIAVMGJ3vxTntIC7eDfIAhtytNIrzQVpyaiKQQAFXx3WmECPRaOI//fYHfY2vMN2R+5/OAlNIUAYFPlCk1OUGa0EtEO3gJBaGsuIr3ajL2diIVmAMCm5qYqjX2ELlhK5jYIQhWMq89S82l3NppCALCR/HKamSK/F4Nj/h1gp6jAidHrzdEUAoDtzEiR+9dijXzQDt4BglAdY6LYxWLanoWmEACsLqOYLziuvBWNA/6dYb+oQxTozWj25p9oCgHA6t5MUJ5ryELd0Q7eGYJQNcPqsKIK2nQBTSEAWNHhPP7LBeXlplhZ9K4QhKphAk2OZm8lyEhCALCeiQfkN6NFL6PaddgxBKGaBoYzJtCac7i+HgCsYkcWP1lAY+vjUH8v2DtqEoimtBDfSlAUdIUAYGmc6JX98qxWzAlH+nvC7lFZn5qCrzN9fwZNIQBY2PJTiijQQCyodj/YQeqb0kKcmqhIiEIAsJxyhaYkKh+0xoJq92duEObm5vbu3dvDwyMkJGT58uW3b7Bo0aI6Nzl37pzp+ePHj7du3drV1bVevXo7d+60VN1a0i1YCHWjb08hCQHAYj47qjTxEWIDkYP3ZzBzu1deeaV69ep5eXkJCQk9evTo0KFDzZo1b94gPz8/Ojp61qxZpi9DQkJMD0aNGvXoo4/u27fvhx9+GDx4cHp6uouLiwXfgDZMixGHbZeH12HOmOEMAA8tv5xmHpK39zH3CK9zZnWExcXFK1eufOONN5ycnNq2bdutW7elS5fevpmnp2fE/xiNRiJKTU09fPjwhAkTGGOPP/64n5/fhg0bLPwONKFdgNDQhxaeQFMIABYwPVkeUIs1rIZ20CxmBWF6erokSfXq1TN92ahRo7S0tNs3W7t2bWBgYIsWLb766ivTMydPnoyIiHBzc6v8xpMnT97tVTjn+fn5V/9HkqQHeysO7t2W4nvJSom+3jQAWF5GMf/6pDKlBc4vmcusxvnq1avu7u6CcOOPCy8vryNHjtyyTVxcXLdu3UJCQvbv3z9q1Ch3d/cRI0bk5+e7u7tXbuPl5ZWXl3e3VyksLGzRokXlq/Ts2XP+/Pl327i0tNRoNBoM2mn8o5wpxtdpTnLF8/UcIAyLiorULkHXsP9VJMtyeXm5Pf+lPmG/8em63EMuKyxUuxQrKC0tLS8vF0VzY97FxcV0hvIezAqS6tWrFxUVKYrCGCOi/Px8f3//W7Zp0qSJ6UGfPn1efPHFVatWjRgxws/P79q1a5Xb5Ofn169f/26v4uXllZaW5uvra05JBoNBY0FIRO+35l03Ss83c/V0hDUgPD091S5B17D/1WIKQldXV7ULubNDeXx7jnSik9EhDiNVIIqis7Oz+UFoDrNOjdasWdPFxaWyC0xJSak8TXpHlV1d/fr1z549W5mFKSkp9whCaOQjdA9mnxzFSCEAVNHEA/JbWFDtAZkVhK6ursOHD58yZUpeXt769ev37NkzYsQIIjp+/Hj//v1N2yxduvTYsWOXL1/euHHjJ598MmjQICKKjIxs06bNlClTrl27Nnfu3OvXr/fu3dt6b0YDprZgHx+Rr15Xuw4AcEDbs/jZQiyo9sDMPbU4a9asF154oWHDhkFBQStWrAgMDCQiSZIuX75s2iAlJeXdd9+9evVq7dq1p0+fbkpKIlqyZMkzzzwTGRkZFRW1du3a+56r1blIb6FfGPvoiPxOSwx0A8AD4EQT9svTH2FG5OADEji3l2Uu/fz8zB8j1N5kmUrni3iLn6Rjg4z+djoGQURUWFiIMSoVYf+ryG7HCJeeUj5LVfb1M2j7momSkhJ1xgjBlmp5CI/XYbMP4569AGCucoWmJiqzW2FBtapAENqjN5uLC08oGcX20qwDgJ379KjS1BcLqlURgtAeBbnR6Cg28xCmjwLA/eWX06xD8rQYHM+rCDvOTk1qJi4/pZwtRFMIAPfxfrI8sDYWVKs6BKGdqu5C/2rA3k9GUwgA93KxmC8+qUyOxjzzqkMQ2q9Xmoo/n1dOFKApBIC7ev2g8mxDFuSmdh2ODEFov6o50YuNxHeT0BQCwJ0dzuNbMpVXmqAdfCgIQrv2f43Z1gzlUB6aQgC4gymJysSmogfWKXk4CEK75mGk/zRBUwgAd5B0he/P5eOxoNpDwx60dy80Yvty+cFLaAoB4G+mJCiTmjNXDa6vZWsIQnvnItJrzdg7SVhoBgD+knCZJ17hT0XhGG4B2IkOYFx9dvQq7ctFUwgAN7yVIL+JdtBCEIQOwInR683ZW3+iKQQAIqI/L/PUq/RkPRzALQP70TGMjmTpxbQtE00hANAbB+U3mjMnHL8tBDvSMRgYTY9h//eHLGECKYC+7cnmp67RaIwOWg52pcN4LJyFutNnqUhCAF2bkii/GY2771oS9qUjmdNGfD9ZzizBCVIAndqTzdOL6Im6OHRbEvamI4n0Fp6ux149gKYQQKfeTJAnt2AGHLktCrvTwbwZLe7J4duz0BQC6M6WDJ5ZQkMjcNy2MOxQB+NmoP+2Yi/8LlegLQTQmXeS5LfRDloB9qjjGRTOwjwwawZAX365yK+U0T/RDloB9qlD+qStOB2zZgD0ZGqiPLUlY7gLvRUgCB1SXS/h6XpsImbNAOjD+nReXEGP1cYR2yqwWx3Vm9HiXsyaAdCHd5Lkd9AOWg2C0FG5GWg2Zs0A6MCa84qk0D/QDloN9qwDeyyc1fKgT48iCQE0ixO9nai83RLdoBUhCB3bnLbi9BQ5oxgnSAG06ceziihQ3zAcq60IO9ex1fUSxtXHrBkAbVI4vZukvNtSRDtoVQhCh/dmc3FfLmbNAGjQD2cUdwP1rokctC4EocNzNdDs1uz5vZg1A6ApMqd3kpS3W4pqF6J9CEItGFib1fakTzBrBkBDvjut+DpTXAjaQatDEGrEnLbiDMyaAdAKmdO0JOXdGLSDtoAg1Ii6XsJ4zJoB0IolaUoNV+oShHbQFhCE2vFGc3FfLt+WiaYQwLHJnKanKNMwOmgrCELtcDXQB62x1gyAw1t8UqnlQZ3QDtoKglBTBmDWDICDq1Do/WTlrWi0g7aDINSaOW3FmZg1A+CwFp1UorypQyDaQdtBEGpNXS9hfAM2AbNmABxQuUIzU5SpLdAO2hSCUINebyb+gVkzAA5o/nGlkQ+19kc7aFMIQg3CrBkAR1Qm04wUZTJGB20OQahNA2qzcE+ag1kzAI7jq2NKTHXhkRpoB20NQahZc9qKszBrBsBBlMk0+7AypQWOySrATtesOl7C+AbsFcyaAXAEc1OVNv5Ccz+0gypAEGoZZs0AOITLZTT7kIx2UC3Y71rmaqCP2mDWDIBd40RP7ZbHRLHGPmgH1YEg1Lh/1GLhnvTxESQhgJ2akaLklnLcd1BFCELtm9NWnHUIs2YA7NEfuXzOEXlFV9GIg7F6sO+1r46X8EwD9vJ+NIUA9uVSGQ3ZKi/uZAjzwElRNSEIdWFSM/HAJf7rRTSFAPZC4ZxBYCwAACAASURBVDRyhzQyUugVihRUGYJQF0yzZl7cJ1+X1S4FAIiIaHqKUiIRlhW1BwhCvehfi0V6Ya0ZALuwK5vPTZWXdxENOAbbAfwSdMR0h6bzRThBCqCmS2U0fLv8dUdDiDtOitoFBKGO1PESnm/IXsVaMwDqUTg9sUMaHSX0xNCg3UAQ6strzcSDl/gvmDUDoJL3kpUyGUOD9gVBqC+YNQOgop1Z/Itj8vIuoohu0J4gCHWnXy1WzxuzZgBsLbeURuyQF3cyBLshBu0LglCPPm6DWTMANmUaGnwySugRghS0OwhCParjJbzQiE3ErBkAW5mWrFxXaDKGBu0SglCnXm0q/olZMwA2sTOLf4mhQTuGINQpzJoBsA3T0OA3GBq0YwhC/TLNmsEdmgCsR+E0Yof0VD0hDkODdgxBqGumOzRh1gyAlbyTJFco9FY0hgbtGoJQ1yI8hRcasQm4QxOAFezI4vOP8+VdDBgatHMIQr2b1ExMzsOsGQALy7lx1aAY5KZ2KXA/CEK9cxbpk7YiZs0AWJDpqsFx9RmGBh0CghCoV6hQ31v4CLNmACzk7URZUuiN5jjAOgb8noCI6OO27IPDcjpmzQA8tC0ZfOFJ/l1XDA06DAQhEBFFeAoTmopDt8sVaAsBHkJWCY3aKS/tLAa4ql0KmA1BCDdMaMoCXYWX92OoEKCKFE4jd0rPNGCdg9AMOhIEIdwgEC3uJG7O4N+moSsEqIopibKs0CQMDToa/MLgL55G+qGr+PJ+OekKBgsBHsz2LL7oBF+OoUEHhCCEv2niK3zSVvznNrmgXO1SABxHTik9sUP+trMYiKFBB4QghFsNrcPiQoSRO2V0hQDmUDgN3y4924B1C0Yz6JAM5m/6/fffL1q0iDE2fvz4AQMG3PJ/jx49umjRopSUFBcXl/j4+LFjxzLGiGjBggUHDx40bePm5vbRRx9ZqnSwno/biF03SrNSlFeb4U8lgPt4K0HmRPjH4rjMDcLNmzc///zzS5YskSRp1KhRNWrUiI2NvXmDDRs2eHp6Tpo0qaio6KWXXsrLy5s0aRIRbd261cPDo0ePHkTk5ORk8TcA1mBk9ENXwyM/S9HVcUNtgHvZlskXn+QJAzA06MDMDcJPP/305Zdf7t27NxE9//zzc+fOvSUIJ06cWPk4MzNzyZIlpiAkoujo6MGDB1uoYLCRIDda0lkcvl060N8Q6o5/4gB3kF1KT+yQl2Bo0MGZ28snJye3bt3a9LhNmzZJSUn32Dg1NTUiIqLyy+XLlw8ePHjSpEmZmZlVLhRsr0uQ8GIjcdBWLEMKcAeSQkO2Ss83Yl0xNOjgzO0Ic3JyfHx8TI99fX2zs7PvtuXOnTu//fbbhIQE05fdu3c3GAxeXl6rV6+Ojo4+fPiwv7//Hb+xqKioc+fOBsONktq1azdjxoy7vUppaanRaKzcGKzkuTr0R7bxxT3SBy2lm58vLi4WBPzjVw32v4pkWS4vL5dleWqK6MrYcxFlRUVq16QnJSUlFRUVomjuLR5dXFzumxTmBomnp2dpaanpcXFxsbe39x03+/PPP4cMGfLDDz/UrVvX9MxTTz1lejBgwIC2bdsuW7bs3//+9x2/183Nbc6cOV5eXqYva9as6eHhcbd6RFFEENrGkq7U+mfpxyynUZF/nT/gnN/jtwPWhv2vIlMQ7rjs8kO6nDjA4OXirHZF+sIYc3Z2Nj8IzWFukNSqVev06dPt2rUjotOnT4eFhd2+TXJycnx8/Pz583v27Hm3H5KXl3e3l2CMNWvWzNfX18ySwDY8jfRDN7HLBqmprxDthy4EgDJK6Knd0oquhuouapcClmDuGOHQoUMXLFggSVJ5efmiRYuGDh1qen727NkXLlwgouPHj/ft2/ejjz7q169f5XcpipKammp6nJKSsmnTpo4dO1q0frCFxj7Cp23FIVvlfFxlD7onKTR6r/hSI7FDIP4u1Ahzg/DZZ591dXUNDw8PDw8PDAysPOH5+uuvnzlzhog++OCDjIyMoUOHCoIgCELt2rWJSJbljh07hoSEREVFdezYcdKkSXFxcdZ5I2Bdj9dhvUKFkTtwlT3o3eQk7m7gE5riqkHtEDh/gCNbRkYGYywoKMj8b+GcZ2RkSJIUGhp67yE9Pz+/tLQ0M0+NYrKM7VUo1G2j1Lsmm9SMFRYWenp6ql2RfmH/q2XTBf6vvfLeXlJoNZwVVUdJSYlqY4QmISEhD/oCgiCEhoY+6HeBHTIy+qGb4ZE1Ugs/od2dJ0sBaNnFYv7kLmlFF+bnjDMjmoLuHh5AoCst6SyO3CmlF2N0BPRFUmjodvk/TcT2AWqXApaGIIQH0zlIeLmJOOp3I66yB115/U/Zy0ivYGhQi/BLhQc2oSmr6U7//gNJCHqx8QJfcYZ/08mAMyGahCCEByYQfdGqYkcW//ok7mUP2nehmD+1S1reRcRVg1qFIISqcDfwn+LE1w7KiZcxawC0TFJo6Db55aZi+wB0g5qFIIQqquctfNJWfGyrfOW62qUAWM1rB+VqTvRyExwqtQy/Xai6f0aw/rWE0TslBW0haNGGC3zlWf5NZwwNahyCEB7K7FZiQTlNT8FgIWjNhWL+9C5peRfRD6tqax2CEB6KgdEP3QxfHlN+uYiuELSjQqHHt8kTMDSoDwhCeFiBrrS0szhmp5RehCwELVA4Pb1b9ncR/o2hQX3ArxksoFOQ8EpTsf9muVS6/8YAdm7CAfnUNb60i4hmUCcQhGAZ/2nCIr0EXGUPju6NP+VtmXxDT4M7lvTXDQQhWIZAtKijuDubL8JV9uCwPj6irDrLf+llqOakdilgQwhCsBgPI62OEycdlBNwlT04oK9PKp+lKtv7iAGuapcCtoUgBEuq5y18FSsOwlX24GiWnVKmJCi/9RaD3TAyqDsIQrCwf9RiA2oJo3bgKntwGGvPKxMOyL/0FiM8kYJ6hCAEy5vVSiysoPeSMVgIDmBbJh+7R17bw9CwGlJQpxCEYHmmq+znHVc2XUBXCHZtfy4ful1a2c0QUx0pqF8IQrCKAFf6oZs4Zpd0thBZCHbqcB7vv1n6uqOhYyBSUNcQhGAtbf2FV5uJA7fgKnuwR6eu8d6/ynPaio/WRArqHYIQrOjfjVmUt/ASrrIHO3OxmPfYJE9twf4ZgWMgIAjByhZ2EH/P4QtPYOIM2ItLZRS3SX6uIXu6Hg6AQIQgBGvzMNLq7uLrf+Iqe7ALBeXU6xdpeB2Ge+1CJXwUwOqivIV5seKgrfLlMrVLAX0rkajvb1JsgPBmNA598Bd8GsAW+tdij9UWhm2XZLSFoJJyhQZukep6CR+3FdWuBewLghBsZMYjYrlC05IwWAgqkDkN3y67G4QFHXBzJbgVghBsxMBoRVfDwhO4yh5sjRON3S1fK+fLu4iIQbgdghBsJ8CVVnQTx+ySzuAqe7AVTvTsXvlcIV8TZ3DGOVG4EwQh2FRbf2FSM3Eg7mUPtjLpoJx4mf/cw+CKG+3CXSAIwdZeasya+grj9+Aqe7C695KVDel8Yy+Dp1HtUsCOIQhBBV/EiolX+AJcZQ/W9Hmqsvik8ltvg5+z2qWAfUMQggrcDfRTd/GNP+U/cZU9WMfSU8qMFOW33mKQm9qlgN1DEII6Ir2Fz9uJg3Eve7CCL44prx1UNj8qhuNGu2AGBCGo5rFwNjhc+OdWXGUPFsOJpibKHxxWtj8q1vNGCoJZEISgpumPiEygdxIxcQYsoFyhJ3bIv17k+/oZIpGCYDYEIahJFGhpZ8PXJ/ma85g4Aw+lsIL6/SaVSLTtUUMNF7WrAYeCIASV+bvSD93Ef+2RT1/DGVKooswS3nG9VM9bWNVNxPWC8KAQhKC+Nv7CG83FgVvkElxlDw/uyFXedq3cv5Ywp63IcEIUHhyCEOzCC41YtB+usocHtjWTd9sozXyETW2B9dOgihCEYC++aC8evcrnHcdgIZjr2zRl+HZpZTfD43VwKIOqw9l0sBeuBlrRVYxdL0X7CY/UwBkuuI+ZKcpXx5UdfQz1q+HTAg8Ff0aBHYn0FubHioO3ypdwL3u4O5nTM3vl5aeV3X1FpCA8PAQh2Jd+tdjjdYSh23CVPdxZsUT9f5POXON74g0h7khBsAAEIdid92NEUaCpuMoebpNdSp3WSwGuwvqeuKEEWAyCEOwOE2hJZ8M3J/lP5zBxBv5yLJ+3XSv1DRMWdhSNOHSB5eDTBPboxlX2e3GVPdywK5t32SC90xKXSYDlIQjBTrXxFyZHi49twb3sgXZl88FbpWVdDE/UxSELLA+fKrBfzzVk0dWFsbjKXt+Sr/AhW6VvOxm6BWNqDFgFghDs2uftxNSr/MtjGCzUqSNXee9fpPkdxJ6hSEGwFgQh2DVXA/3QTZycIO/NwWCh7py6xnv/In/URowPw5EKrAgfL7B3db2EhR3FYdtxlb2+XCzmPTbJk6MZlk8Da8MnDBxAfBgbXkd4HFfZ68alMorbJD/bkI2tj2MUWB0+ZOAYpsWIRkaTEzBxRvvyy6nnJml4HfZKExygwBbwOQPHwARa3sXw3Wm+GlfZa1qxRH1/lToFCW9G4+gENoKPGjgMX2f6vqs4fo98PB9nSLWpVKK+v0r1vIUP2+CqebAd3IYJHEmrGsI7LcUh2+Q/+hnc8OHVlgqFBm+VQt2F+R1EXCphKR999NGyZcvUrsIyhg0b9p///McaPxnHEnAwzzRgf+Ty8XvkBR1EZ7QNWiFzGr5ddhKFrzuKDDFoOcnJyf369evTp4/ahTysjRs3JicnW+mHIwjB8XzRXhyyVQpYVtE9hMWHCX1qsuouatcED4ETPb1bLijna3sYDBiusbTatWu3bNlS7SoeVmpqalpampV+OIIQHI+bgdb3NORdp62Zyrrz/N9/VER4Cn3DhPgw1rI6ugnH88p++UQ+/623AS0+qAJBCI7K15kGh7PB4SQp4h+5fOVZZeAW2SCQKRE7BQm4U49DeP2gvC2Tb3vU4IH7C4JKEITg8AyMYgOF2EBxTls6epWvT+dTE+Xj+bxrMOsbJvSvxbyd1C4R7uLjI8rqc3xnX4OPs9qlgI4hCEFTGvkIjXyEV5uxS2W06YKyPp2/uK+isY8QH8YG1BaivHHi1I58nqrMTVV29hUDXNUuBfQNQQjaVMOFRkaykZFUKolbMpX16bzLBqWaE8XXEvrWZO0CBExNVNeSU8qMFGVnXzHYDb8JUBmCEDTO1UDxYSw+jBROSVf4unTl//6Q04t4r1AWX0voHcowNGV7a84rrx6Qtz5qCPdECoL6EISgF0ygltWFltXFqS3obCHfnMG/TVPG7pYfqSH0rckGRwhoTWxjcwYfv0fe1NPQoBp2ONygKMr333+fkJBw4cKF6dOn16lTx5avjnl1oEfhnsK4+mxdD8PZfxrH1WcJl3mTH6WYNdLURDnhMtZvs6Lfc/iIHdKP3QwtcKEL3ESSpKVLl7q5ua1bty4vL8/Gr44gBF3zcabB4ezbzmL2cOPMVuLV6zRoq9zqZ2l9OuLQ8g5e4gO2SN93NcQGIgX16/333z9//nzll++99156erqTk9PGjRvfffddo1GFsQoEIQARkZFRt2BhTlvxzD8Nb0WzqYlys9XSyrMK8tBSjlzl/X6TFnYwdAlCCupaZmbmF198YXq8f//+uXPnBgUFqVsSxggB/kYgig9jfcPY+nTl7URlWpLyZjQbFI5Jpg/l1DXe+xf5ozZi3zDsSJXNOaJ8ctRG9zITBJr5CHss/G8d13PPPdepU6e3337b2dl5/vz5Tz/9tCpd4M0QhAB3cHMcvpOovJukvIU4rKqLxbzHJnlyNHu8Dk5Bqe/Jeiy+lu0+yLfPQWvQoEGDBg1+/vnn3r17r1q1KiUlxWbF3A2CEOCuTHEYH8a2ZPBJB2XEYRXkllLcJvnZhmxsfaSgXfA0kqdR5Y/wM888M2/evCtXrsTGxtaqVUvdYghjhADm6B4iHPiHYVpLNiNFafGT9NM5jB2aJbOEd98kDa/DXmmCQw38ZeDAgUePHp05c+a4cePUroXogTrCbdu2bd261d/ff8yYMV5eXrdvcPbs2WXLlsmyPGTIkAYNGpie5JyvXr364MGDdevWHTVqlOrnggGqRiDqV4vF17oxdvh2ojK5BRtQG93hrU4U8H05fE8O35vDM4v5v5uwN6ORgvA3Tk5Oo0ePXrp06c03SoyMjDx16hQRtWrViojOnz8fFhZmm3rM/YAuXrx4+PDhPj4+u3fv7tChgyRJt2xw9uzZli1b5ufny7Lcpk2bQ4cOmZ5/4403Jk+eHBAQsGTJkiFDhliydgCbM50s/fMfhtmtxZkpSpMfpW/TFFnf7aGkUMJlPueIMmSrHLCsossGaeMF3sRHmB8rXn7COLUFbq0Ed5Cenj527FhR/OvjkZaWxm9isxQkMztCzvn777//+eefDxgwQFGUxo0br127duDAgTdv89lnn/Xr12/27NlEVFZW9uGHHy5evLigoODTTz/dv39/w4YNx44dGxQUdOTIkcaNG1vlrQDYUPcQoXuIYcMF/naiPC1ZGVRbGBjOWlYXdNIg5l2n33P43hxlbw5PusIjvYTYQGFgbeHDNoZQd53sA6ii5OTkJUuWbN68ec6cOWrXcoNZQXjx4sW0tLSePXsSEWOse/fuO3bsuCUIt2/fPnHiRNPjnj17Pvnkk0T0559/+vj4NGzYkIg8PDzatWu3fft2BCFoRp+aQp+ahoOX+E/nlCd2yCUS/aOWMKA26xAoiJqLg7QC/nsu35vD92bzi8W8tb/QPoBNjmat/QVPjHjAgwgLC9u7d2/16tXVLuQGs4IwKyvLw8PDzc3N9GVAQEBSUtIt22RnZ/v7+1dukJ2dzTnPysqqfNL0fFZW1t1epbS09KWXXnJ2vnFfspiYmNGjR99t47KyMlmWDQbMelVHWVkZhnsrNfGkJk1ochM6W0QbLgivHxBOFFCvEOpTk/cM5u5W+JDaZv9LCp28RvsuCb/nCrtzhQqFt/CjdjXo8VY8pjp3qhxXkalMtnYtdkSW5fLycsFBmn9ZtrvfTfPmzZs3b16Fb5RluaysrKysjHN+8znVezMajffd2Kx/owaDQVH+ugBTkqTb/xGKoli5xyVJEkVREIRbvlGW5Xv86xVFsVmzZu7u7qYvw8PD71G9+D/m1A8Wh51/R3W96SVveqkxXSimXzNo+VnhmX3UIZAeq0X9wsjLcsllvf1fWEEHLtHvl+j3XNp/iWq6UfsA6h5CE5tSw7/WyHaMDLAeB/r8M6admUqMMfEmZn6XOX+ymBWEwcHBJSUlV69e9fHxIaLMzMzbV8QJDg7OzMw0Pc7MzAwODiaioKCgzMxMzrmplIyMjA4dOtztVZycnJ588klfX19zSjKFMTpCtRiNRnSE9xBRjZ6pRs80oivXaUO6svKs8n/7eSt/oW9N9s86LPCh70Nr2f2fWcL35vA92XxvDj9ZwJv6CrGBwkuNhdgAhhvH344xxjl3lM+/o3Su5hAEwfg/lv1DxKw/FgIDA2NiYn788UciKikp2bRpU9++fYmosLBw7969pm369OmzatUq0+Mff/zRtEGbNm1kWd6zZw8RZWZmHjhwoHfv3hasHsDO+TnTyEi2rocha7jxxUYs4TJvuKoidp00M0U5dU21yaaSQkev8nnHlZE75NrfSzFrpG/TlGA34eM24uUnjHviDTMeEePDkIKgF+Z2VNOmTRs2bNj+/fuTkpKaNm3apUsXIjp8+HBsbCznnIieffbZNm3axMfHu7m57dmzZ9++fUTk7Oz87rvvDhkypH///lu3bh0/frw9LCIAYHtu/7s/cJksbs5QfjrH266Vwz2FuBAh2E0IdiN/VyHAlYLcBGuMKRJRQTn9nsN/z1X2ZPOEyzzcU4gNFHqGCu/GsFoe2mkaAKpAMMWYOc6dO7dz587g4OCuXbua2tLCwsJDhw61b9/etEFRUdGvv/4qSVLPnj2rVatW+Y2HDx9OSEiIiopq167dPX6+n59fWlqamadGS0tLcWpURYWFhZ6enmpX4dhkTjuz+N4cnlPKM0sot5TnlFJWCedEQW5CoCvVcBFC3KmGixDoRoGu5O8qBLlSgKvgajB3/58t5Htzbvx3rpDHVBdiA4V2AaxdgGDBAUu9MU2WcXV96BPcNjFu3Liff/7Zw8ND7UIeVlFRUb9+/ebPn19SUuLs7GzZU6MPEITWhiB0IAhC6ymWKLuEZ5fSpTKeUUyXynh2CWWV0qVSnlVKOaXcIFCgCw90Z4GuQqAb+bsIQW4U4Cr4u1KwG2WX0t7sG+HHibcPYLEBQrsAIdpPMGhn2oSaHCsIS0pKsrOz1a7CMkJCQpydna0RhAgSAPvibqA6XkIdL7rb5MxrFXTqUlExc88u5dkllFvGf8+h3DIlp5SySqiaE8UGCvG1hJmtWLgnznnqnZubW0REhNpV2DsEIYCD8TJSpCf39BRwGQOAReBcCQAA6BqCEAAAdA1BCAAAuuaoQbhx48bDhw+rXYV+zZ07t7i4WO0qdKq0tPTTTz9Vuwr9Sk1NXbdundpV6NeaNWuOHz9u2Z/pwEG4f/9+tavQr3nz5mlmTrbDyc3N/fLLL9WuQr8OHDiwYcMGtavQr3Xr1iUkJFj2ZzpqEAIAAFgEghAAAHQNQQgAALpmR0usubq6BgYGmnn3rMuXLzs7O2OVL7VcvHgxMDAQS9ypQpblzMzMmjVrql2IThUVFZWWltaoUUPtQnTq0qVLbm5ulXeuva9hw4a9++67997Gjg5kp06dun79upkbV1RUiKKopXtOOpbr1687O+MmParB/leRoij3vsc4WNWDHvxvv3vu7eyoIwQAALA9dFQAAKBrCEIAANA1BCEAAOgaghAAAHTNjmaNmqOsrOzQoUNHjhzx9/fv27fvHbe5evXqggULsrKy4uLievfubeMKNW/Pnj1r1qzx8fEZM2ZMcHDwLf/39OnTW7durfwyPj7enClbcA+bNm3avHlzUFDQ2LFjq1WrdvsGJ06cWLJkiSzLw4cPb9y4se0r1DDO+fLlyxMSEiIiIp5++mkXF5dbNtizZ09qamrll+PGjbNtgVpWVlaWnJycmpoaFBR0tyP5lStXFixYkJOT06tXrx49elT5tRysI/zvf/87YsSIjz/++KOPPrrjBpIkdezYMTExMSIiYvz48fPmzbNxhdq2YcOG/v37h4SEZGRktG7dOj8//5YNDh48OGPGjDP/U1ZWpkqdmjFv3rzx48dHREQkJiZ26NBBkqRbNkhLS2vdurUgCB4eHu3bt09JSVGlTq167bXXZs6cGRkZuXbt2kGDBt2+wfLly5ctW1b5gbd9hRr23nvvjRo16sMPP7zbEvPXr19v37794cOHw8PDx4wZ880331T9xbhDkWWZc/7ZZ5917dr1jhusWrWqfv36ps3Wr18fHh5uegwW0b59+y+++ML0uFu3bh9//PEtG3z33Xc9evSweV3aJMtyeHj4unXrTI8bNGjw448/3rLNCy+8MHbsWNPjiRMnjhw50tZVald+fr67u3tqairnvLi42MvLKykp6ZZtnnnmmWnTpqlRnfaZDt0ffvhh796977jB8uXLmzRpoigK53z16tX16tUzPa4CB+sI73sR5a5du7p162barHv37ufPn09PT7dJadpXUVGxb9++uLg405dxcXE7d+68fbPMzMzZs2cvXLgwJyfHtgVqTXp6+rlz57p3705EjLFu3brdvsN37dp1398IVE1CQoKPj0+DBg2IyM3NrX379rt27brjZrNmzVqxYoX564GAOe57tN+5c2dcXJwgCEQUFxd34sSJKt8Sx8GC8L6ysrIqlz5ydnb29vbOyspStyTNyMnJURTF39/f9GVAQMDt+9bT07NZs2b5+fnr1q1r0KBBYmKizcvUjqysLC8vr8pxqYCAgMzMzNu3qfzA+/v7Z2dncyyRYSHZ2dk3r6N2x/0fGhoaGBiYn5//wQcftGjR4tq1a7atUddu/vB7eHi4ublV+Whvd0E4a9Ysw218fHzM/HaDwSDLcuWXFRUVTk5O1qlUm/71r3/dvv9jYmKIyLSmVOUw1R33bZ8+fZYuXTpt2rQ1a9aMHDly8uTJNq5fS4xG482DghUVFbcvq2YwGCq3kSTJYDCY/kCGh3f7weT2/f/6669//vnn77///r59+1xdXb/66ivb1qhrt/wDkSSpykd7uwvCiRMnSre5evWqmd9umsdhelxQUFBUVHT7zEa4hy+//PL2/f/nn38SUY0aNYxGY+XuzcjIuPeM0Hbt2mH6wMMIDg4uLi6unJF0xx0eEhJS2aZkZGSEhITYtERNCw4OzszMrOyw7/2BF0WxTZs2+MDb0s0f/tzc3PLy8iof7e0uCKtm+/btBQUFRBQfH//LL78UFRUR0apVq2JiYjB931IYY3369Fm5ciURSZK0Zs2afv36EdH169e3bdtWXl5ORJXTRDnn69evx2z+hxEcHBwTE7Nq1SoiKioq2rRpk2mHFxQUbN++3bRNfHy86TdCRCtXroyPj1erWu1p1aqVKIqmXZ2RkbF///5HH32UiLKzs//44w/TNqWlpaYHJSUlW7duxQfe2jjn27ZtKywsJKL4+PgNGzaUlJQQ0apVq9q3b+/r61v1n+tAtm7d2rJly7CwME9Pz5YtW77++uum5xlju3fvNj0eOHBg06ZNR44cWb169c2bN6tXrAYlJSVVr1592LBhbdu2jY2NLSsr45yfO3eOiEx/O/fv379Lly4jRoxo3rx5RETEqVOn1C7Zsf32229+fn4jR45s2rTpwIEDTU/u3r2bMWZ6fOnSpaioqJ49e/bv379WrVoXL15Ur1gNWrhwob+//5gxYyIiIiZMmGB6cvHixZGRkabHQUFBffv2HT58eGhoaFxcnOlfBFjEpk2bvHv3jgAAAWZJREFUWrZsGRoa6u3t3bJlyylTpnDOKyoqiOjAgQOcc0VR+vbt27x585EjR/r5+e3YsaPKr+Vgd5+4evXqzScffH19w8PDiejgwYMNGjTw8PAgIs75zp07s7OzY2NjQ0NDVatVoy5durR9+/Zq1ap16dLFNGpYXl6enJwcHR1tNBqvXr164MCBvLy84ODgtm3bYoD24V28eHHPnj2BgYGdOnUyjf8VFRUdO3bskUceMW1QUlKyZcsWWZa7d++OO3Ra3IkTJxITE+vUqdOqVSvTM5cvX87KymrSpAkRpaenJyUllZWVRUZGtmjRQtVKtSYvL+/s2bOVX/r5+dWuXZuIDhw40KhRI9P9CBVF2bFjR25ubseOHR9mFMzBghAAAMCyNDJGCAAAUDUIQgAA0DUEIQAA6BqCEAAAdA1BCAAAuoYgBAAAXUMQAgCAriEIAQBA1xCEAACgawhCAADQNQQhAADo2v8DJhMMqV8JuGEAAAAASUVORK5CYII=",
- "image/svg+xml": [
- "\n",
- "\n"
- ],
- "text/html": [
- "\n",
- "\n"
- ]
- },
- "metadata": {},
- "output_type": "display_data"
- },
- {
- "data": {
- "text/plain": [
- "([0.9999999959286234, 0.9999999907475938], \"optimal\")"
- ]
- },
- "execution_count": 29,
- "metadata": {},
- "output_type": "execute_result"
- }
- ],
- "source": [
- "NWTN(TF[6], printing=false, plotatend=true, Plotf=1)"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "id": "ae61290b-22fa-493a-9b56-3a7458e3b3ed",
- "metadata": {},
- "outputs": [],
- "source": []
- }
- ],
- "metadata": {
- "kernelspec": {
- "display_name": "Julia 1.9.4",
- "language": "julia",
- "name": "julia-1.9"
- },
- "language_info": {
- "file_extension": ".jl",
- "mimetype": "application/julia",
- "name": "julia",
- "version": "1.9.4"
- }
- },
- "nbformat": 4,
- "nbformat_minor": 5
-}
diff --git a/project/thinQR.jl b/project/thinQR.jl
index 31be4dd..cf4d561 100644
--- a/project/thinQR.jl
+++ b/project/thinQR.jl
@@ -1,7 +1,7 @@
module thinQR
import Base: size, show, getproperty, getfield, propertynames, \, *
-using LinearAlgebra: norm, I, triu, diagm, ldiv!
+using LinearAlgebra: norm, I, triu, diagm, ldiv!, dot
export QRthin, qrhous!, qyhous, qyhoust
@@ -18,16 +18,22 @@ function qrhous!(A::AbstractMatrix{T})::QRthin{T} where T
m, n = size(A)
d = zeros(n)
- for j ∈ 1:n
- s = norm(A[j:m, j])
- # iszero(s) && throw(ArgumentError("The matrix A is singular"))
- d[j]= copysign(s, -A[j,j])
- fak = sqrt(s * (s + abs(A[j,j])))
- A[j,j] -= d[j]
- A[j:m, j] ./= fak
+ @inbounds begin
+ for j ∈ 1:n
+ s = norm(A[j:m, j])
+ # iszero(s) && throw(ArgumentError("The matrix A is singular"))
+ d[j]= copysign(s, -A[j,j])
+ fak = sqrt(s * (s + abs(A[j,j])))
+ A[j,j] -= d[j]
+
+ A[j:m, j] ./= fak
- if j < n
- A[j:m, j+1:n] -= A[j:m, j] * (A[j:m, j]' * A[j:m, j+1:n])
+ if j < n
+ A[j:m, j+1:n] -= A[j:m, j] * (A[j:m, j]' * A[j:m, j+1:n])
+ #for k ∈ j+1:n
+ # A[j:m, k] -= A[j:m, j] .* dot(A[j:m, j], A[j:m, k])
+ #end
+ end
end
end
@@ -47,7 +53,8 @@ function qyhoust(A::QRthin{T}, y::AbstractArray{T}) where T
m, n = size(A.A)
z = deepcopy(y)
for j ∈ 1:n
- z[j:m] = z[j:m] - A.A[j:m, j] * (A.A[j:m, j]' * z[j:m])
+ # z[j:m] = z[j:m] - A.A[j:m, j] * (A.A[j:m, j]' * z[j:m])
+ z[j:m] -= A.A[j:m, j] .* dot(A.A[j:m, j], z[j:m])
end
return z
end
@@ -101,7 +108,7 @@ function (\)(A::QRthin{T}, b::AbstractVector{T}) where T
v = qyhoust(A, b)
x = zeros(m)
for j ∈ m:-1:1
- x[j] = (v[j] - x[j+1:m]' * A.A[j, j+1:m]) * A.d[j]^-1
+ x[j] = (v[j] - dot(x[j+1:m], A.A[j, j+1:m])) * A.d[j]^-1
end
return x
end