{
"cells": [
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Convex Optimization in Julia\n",
"\n",
"## Madeleine Udell | ISMP 2015"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"## Convex.jl team\n",
"\n",
"* [Convex.jl](https://github.com/cvxgrp/Convex.jl): Madeleine Udell, Karanveer Mohan, David Zeng, Jenny Hong"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"## Collaborators/Inspiration:\n",
"\n",
"* [CVX](http://www.cvxr.com): Michael Grant, Stephen Boyd\n",
"* [CVXPY](https://github.com/cvxgrp/cvxpy): Steven Diamond, Eric Chu, Stephen Boyd\n",
"* [JuliaOpt](https://github.com/JuliaOpt): Miles Lubin, Iain Dunning, Joey Huchette"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"# initial package installation\n",
"Pkg.add(\"Convex\")\n",
"Pkg.add(\"SCS\")\n",
"Pkg.add(\"Gadfly\")\n",
"Pkg.add(\"Interact\")"
]
},
{
"cell_type": "code",
"execution_count": 28,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"problem status is Optimal\n",
"optimal value is 43.82521706423979\n"
]
}
],
"source": [
"# Make the Convex.jl module available\n",
"using Convex\n",
"using SCS # first order splitting conic solver [O'Donoghue et al., 2014]\n",
"set_default_solver(SCSSolver(verbose=0)) # could also use Gurobi, Mosek, CPLEX, ...\n",
"\n",
"# Generate random problem data\n",
"m = 50; n = 100\n",
"A = randn(m, n)\n",
"x♮ = sprand(n, 1, .5) # true (sparse nonnegative) parameter vector\n",
"noise = .1*randn(m) # gaussian noise\n",
"b = A*x♮ + noise # noisy linear observations\n",
"\n",
"# Create a (column vector) variable of size n.\n",
"x = Variable(n)\n",
"\n",
"# nonnegative elastic net with regularization\n",
"λ = 1\n",
"μ = 1\n",
"problem = minimize(norm(A * x - b)^2 + λ*norm(x)^2 + μ*norm(x, 1), \n",
" x >= 0)\n",
"\n",
"# Solve the problem by calling solve!\n",
"solve!(problem)\n",
"\n",
"println(\"problem status is \", problem.status) # :Optimal, :Infeasible, :Unbounded etc.\n",
"println(\"optimal value is \", problem.optval)"
]
},
{
"cell_type": "code",
"execution_count": 29,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/html": [],
"text/plain": [
"Slider{Float64}([Input{Float64}] 2.5,\"λ\",2.5,0.0:0.1:5.0)"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/html": [],
"text/plain": [
"Slider{Float64}([Input{Float64}] 2.5,\"mu\",2.5,0.0:0.1:5.0)"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"image/png": [
"iVBORw0KGgoAAAANSUhEUgAAAhcAAAF6CAYAAACqW3pRAAAABmJLR0QA/wD/AP+gvaeTAAAgAElEQVR4nOzdeYAU1bk+/vc9Vd2zwAgIiCBMdwPigkxVD6OIBtSgJkYlGjNuUYNEcWGMMRrJdhOTm83ExCxgjCaKmkQjWV1iFiPuiM50VfWIosB0dYOgIrLO2lXn/f3RjHK58Re4HPpM1/f9/AXDTM3T3fVUH6pO1wFgjDHGGGOMMcYYY4wxxhhjjDHGGGOMMcYYY4wxxhhjjDHGGGOMMcYYY4wxxhhjjDHGGGOMMcYYY4wxxhhjjDHGGGOMMcYYY4wxxhhjjDHGGGOMMcYYY4wxxhhjjDHG/q0VK3IHExHqzsH2XS6Xq85m88N052BqZDKFMbozMDVaW9fXtrauGaI7B1OjnN0U5fpFqvX14V3Z7Fu1unOwfbd5M0yTkq7TnYOpgUj3687A1IjFgpNjMeNy3TmYGuXsZsUOLhBhdWdnUerOwfadYeB2AFqvOwdThV7XnYCpgRhsIcI3dedgqpSvm2W7rNDy6ZbjAQAW3rPwOQCAG+fOresMa05DxL5YZ+zvty65tbtcWRhjjDG2/5TlzMXnLpo3GoH+DAAnAgDMmzcv1iVrngGAZgA4pVjb+9DebtPzCjOWLiVTbVKmQyazaqTr5ifrzsHUcJzcibozMDUymcIYz+s4THcOpkY5u1mOwQUGZvxOQvoXIhEAQLwYPw8Q8ovuWdS8cPHC+YA4ouXiltTebFRK+aUDD3yrav9EZuVlHglA5+lOwVQR39CdgKlhGNAEIM7UnYOpUr5u7vf/+V8zZ/4NQHKpFBhH2nkZhmAyAHrvfRPhq2jKqQCQ29PtItIS09xUVB6YlZ1hYIFIPqs7B1MDUf5WdwamRhDIVbGY3KI7B1OjnN3cr4OLlrktTUTypIWLbzt9/qXzv/jePxCNQoT33kwQaQMAHLw327bt1N3qkjKdLCuZg70YWLKBzbZTv9CdgakxdWryVQB4VXcOpkY5u7l/z1xIuAkB4/PnzL8PCSYTAM6/dP4bRJAFoEPe/zY40Ajhkd1/3HH8vwoBIy0reXR7e8eoMBSPAkDWtpNzXde/BQBOAqBf23bqVtfNXQeAn5ISv9jYmHjcdfN3A9AUw5CnT5ky/i3P81+SEjam08mPZbNrp0gZ3g1A/7Lt1ALPy11AhNcD4I9tO/Fr181/F4BOkVJ+prFxvOe6/sMAMNqyEse0tW040DT7/kZEK9Lp1Kdd1z8JAH6ACA9YVvIWz8t/loguAaCv2Hbq747j34kIaSnF7MbG+vWu678AQNtsO3VqW5t/hGHAfQDwpG0nb3Dd/LkAdCMRLUynU4td1/9vADgNEeZZVjLjefk/E9HYWKzzuDAcPEhK+icArLTt5EXZbMdMKcWPAOj3tp36nuP48xHhUiL4ejqdfNR187cDUJNhBGdPmTJxrev6zwFAj20nZ7numkkAxm8R8RnLSlznef45RPAlRPi5ZSV/5Ti5mxDxDCJxVTpd/5Lr+n8AgERXlzFDiJ5YdXXsCQBYZdvJCxwnfzwi/YQI/pROJ7/turkrAfAyRPFNy6p/yPP8RUQwDRGaLSuZc93c04gYJ5JfMwxzTRjK3wHAMttOXuM4ubMQ8asAeIdtJ+5wXf+rAHCWlHBNY2Nymef5DxLB+Fis9sTu7kCaZt/TRJBLp5PNnpc7lggXAuBDtp34puf5lxHBlQDwbdtO/sl1/Z8CwHFhaJw/deq41a7rPwEAMdtOznCcXBIRf08EL6bTyas9L38mEX0dkX5lWamfe17uS0R4DgB+zrYTz7qufz8AHBoE4awhQ2RPZ2fsOQAo2HbyE45TOBpR/hyAHrXt1NddN3cpAM5HxO9ZVuL3npe/lYhmIMpPWdb41zzP/ycADLKs5HGet24sUfBnAGqz7dQV2ax/mpTw3wCw2LaTCx0nfyMinSslXd/YmHrKdf37AOCIWMw49e23x24bNiz/AgCtt+3UbNf10wBwJyL8zbKSX/W8/CVE9FlEusWyUg84jv8DRDhJSryksTHxiuvm/wYgD7Tt1DGtrf5o04SHAcC17eRlrttxKoD4DhH9Op1O/dh1/esB4AJEuNGykk94nv88EcQRg481NEzY6Hn5lxDhLctKnt7e3tEQhuIuAPynbSe+5Dj+pxDhOiK4NZ1O/sZ1c98DwJOFMC5taBjX7nn+o0QwyrISR2eza0YSmX8FwHbbTlzqeflZRHQzANxv28kfOk7uc4h4ERF+KZ1O/NN1/V8BgBUEcGZTU3KD6+ZeBBDv2nbio6U5PnQPESxNp5Nf8Lzc+UR4AyL+1LIS9zqO/21E+AgAXG7bScd1cw8B4JjNmxPHHnTQugOKxfAfAPCqbScvdpzciYh4CxE+mE4nvu+6fgsAzBEC/quhIfmY5/l3EEEjonmWZY1d53n+8wDQaVnJU7LZ3OFS4q8B4GnbTn7e8/KfJKIvAtAi207d7bq5bwDg6VKKKxsb61td1/8TAIwbNKh4/Natoto0jX8RwevpdPJCzyvMIJK3ItIfLCv1Xc/LXUWEn0HEb1hW4mHXzf8cgI6WUp7T2Dg+77r+MwBQtO3khx1nzaGIxv0A9Jxtp651Xf9sAPgKAPxFCPmUlMaJADQbkeZbVmq54/hLECEVBPGZNTWmKBa7ngSgNbadOi+T8acLAT8DgD/bdvJbrpufB0DziOhb6XTqz67r/wwApguB5zY0JDo8z3+SiIRtp2Z6np8igiUA8IJtJ1s8rzCbSH4NgH5p26nbHcf/CiKcTQSfTaeTz7tu/gEAmtjTU/ywlNXF2trwGQDwbTv5SdftOAZA3EZEj6TTqZs8z/8MEVyFCN+1rOQfHCf3Y0T8EFF4QTo9YZXr+v8CgGrbTh7f2rqm3jSNPwLgS7aduMpx/NMR4RtEcHc6nVzkurkvAuAnEcV1llX/jOf5vyGCw4TAUwxjR2exOOh5RFxnWYmzstn8VCnpFwDwmG0n/8txcnMQsQUAv2/biQdd1/8hAJwQhnDx1KnJV1039w8APMC2k8dmMoUxQsiHiMBJp5OXe17uo0T4LUS817ISP/U8/wYiOB8AvmDbyaWu698LAEcGQfyjU6eOftfz8i8CwAbbTp7pODkbAP4CIH6TTie+7Lr5iwDoc4j0Q8tK3e95/vfDsPiDxsZDN+7jOz8A7OdPi3xuzpXJAMy60i+iOQQQA4HfgxBSgPSthYsXnXTVhVcNE3GRAQlNi+5dtGlPt10aeFQ3W9bBnfvvEbByyGRyJwiBs2w7+TXdWdi+c5z8U+l04gTdOdi+87zCbAA5ybKSt+jOwvZdObu5X89c/Hjx7X7/n6+ZM/9NBIj97K6F65ubm98aNWjkW9fMmf8iAdQjwXd+thcDCwAAIcR33313VK/y0EyD4BWA+Du6UzBV5Nd1J2BqhCG0GoZ8TXcOpsr/I91smdsy5sa5c+t052CMMcYYA8fJfyeXy1XrzsH2XXt7R4Pj5ObozsHUcF3/R7ozMDUyGX+64+SbdedgapSzmxV8Eyqyt22rMXSnYPuuWMRhQuB43TmYGkQ4VXcGpoZhiJEAMqE7B1OjnN2s2MFFPE5zjzxyVJfuHGzfDRsGy7dvx6zuHEwNIrxAdwamRrFoPg7Q/ZTuHEwN7iZjjDHGWLk5jv9Xz3tzkO4cbN9lMrkTXNf/pu4cTA3HyfP/dCPC8wqzPc+/QXcOpkY5u1mxS64LATuKRUm6c7B9hygCROjRnYOpgUjbdWdgaiDKIhF3Myq4m4wxxhhj5dbWtnYiEVXsmRf2vpUrV9ZlMoUxunMwNXiJ7uhwnNzQ9vaOUbpzMDXK2c2KfXMWIvxpNvtWje4cbN91dVU1CiGv1J2DqSGlcYfuDEwNIYyZUoqLdedgapSzmxX7UVQh4Nna2u2B7hxMBbERQK7QnYKpgSif1J2BqRGGsN4wkOe2RQR3c08spRZYSrwUMGOMMTbAVOxlkelm11FAlXvmhb2vtdUfnckUmnTnYGp4Xv5M3RmYGplMR6K9vaNBdw6mRjm7WbGDi6PMYAYC7dcl41l5CEGThJCzdedgakgJfF+EiDAM05JSnKo7B1OjnN2s2P/5vx7ElwPiDN052L6LxcIOKeNSdw6mCv1CdwKmRhDgK4Yh39Sdg6nC3fzPnqDPwlJ6XXcMxhhjjP1PFXtZ5OzqzlMQgC+LREAmkz/SdXPn6c7B1HDd3Dd0Z2BqZDKFJp5DEx3l7GbFDi5GYDgeiOdcRIMcCYCTdadgahCJE3VnYGoYBowBIL4pWkSUs5sVO+fi+WLVo4R4lu4cbN+Zpsgi4jrdOZgaRHi97gxMjTDsXQZgeLpzMDW4m3uC51wwxhhjA1LFXha5uHrbhTznIhpct+MY181dpzsHU8N1/ft1Z2BqeF5+luf5l+nOwdQoZzcrdnBRjVDHcy6iQUqsAcBhunMwNYiQF6GLDBwEAEN1p2BqcDf3BF8WYYwxxgakij1zUSdkle4MTI2lS8l8/vm1vMJtRKxcubJOdwamRmtrayyXy1XrzsHUKGc3K3ZwcX7VjjlIVLH52fuGDPGPr60Nv6Q7B1Oju7vmEd0ZmBqx2EGnbduGLbpzMDXK2c2KfXPuJHwHkJcCjgLDwO0AtF53DqYKX66MCsRgCxHy7b8jg7v5n/GcC8YYY2xAqtgzF2mjdwIQfxQ1CjKZVSNdN8936IwIx8mdqDsDUyOTKYzxvA6+Q2dElLObFTu4aIr3nsxLrkeFeSQA8doikSF4bZGIMAxoAhC8tkhklK+bFTu42CDNV3nORTQIITcgilbdOZgaQtDDujMwNcIwyEuJfPvviOBu7gmec8EYY4wNSBV75uLUeOc0vv13NHien3LdjlN152BquG7uCt0ZmBptbf4R2WzHTN05mBrl7GbFDi4SRmjz7b+jIQypHkB8SHcOpgaRuFB3BqaGaYpDicQxunMwNcrZzYpdcr21r+pxAvyY7hxMheAVgPg7ulMwVeTXdSdgaoQhtBqGfE13DqYKd/M/4zkXjDHG2IBUsZdFmqt2nMFzLqKhvb2jwXFyc3TnYGq4rv8j3RmYGpmMP91x8s26czA1ytnNir0sMlTIQ3jORTQUizhMCByvOwdTgwin6s7A1DAMMRJAJnTnYGqUs5sVe+bioZ7a+wlR6s7B9t2wYbBcCLxVdw6mBhFeoDsDU6NYNB8vFsM7dedganA39wTPuWCMMcYGpIo9c3F5zbYreMn1aMhkcie4rv9N3TmYGo6Tf0p3BqaG5xVme55/g+4cTI1ydrNi35wDoj6+/Xc0IIoAEXp052BqINJ23RmYGoiySMTdjAru5p7gyyKMMcbYgFSxZy4SIhihOwNTY+XKlXWZTGGM7hxMDV6iOzocJze0vb1jlO4cTI1ydrNiBxenVnWdw3MuoqGrq6pRCHml7hxMDSmNO3RnYGoIYcyUUlysOwdTo5zdrNg35y1SvMFzLqIhFqPNRNShOwdTA5HadGdgaoSh3Cgl5nXnYGpwN/cEz7lgjDHGBqSKPXMxPdZ1JBDf/jsKWlv90ZlMoUl3DqaG5+XP1J2BqZHJdCTa2zsadOdgapSzmxU7uDjKDGYg8O2/o0AImiSEnK07B1NDSuD7IkSEYZiWlOJU3TmYGuXsZsWuLfJ6EF8OiDN052D7LhYLO6SM863cI4N+oTsBUyMI8BXDkG/qzsFU4W7+ZzzngjHGGBuQKvayyBlVXTN5yfVocN01kzyvwJdFIoJvFx0dmUyH5Tj5U3TnYGqUs5sVO7gYLYIjeMn1aJBSjCaSPKEzIqREntAZEYZhJoQgS3cOpkY5u1mxcy6eL1Y9Sohn6c7B9p1piiwirtOdg6lBhNfrzsDUCMPeZQCGpzsHU4O7uSd4zgVjjDE2IFXsZZGLq7ddyHMuosF1O45x3dx1unMwNVzXv193BqaG5+VneZ5/me4cTI1ydrNiBxfVCHU85yIapMQaABymOwdTgwh5EbrIwEEAMFR3CqZGObtZsYOLB3oHLyZEvjdCBGzdmnyuq8v4ru4cTI2amu4zdGdgahSLbz92wAG0UHcOpkY5u1mxgwvGGGOMDUwVO7g4v2rHHF5yPRqGDPGPr60Nv6Q7B1Oju7vmEd0ZmBqx2EGnbduGLbpzMDXK2c2KfXPuJHyHl1yPBsPA7QC0XncOpgp/iisqEIMtRMi3/46M8nVzv0+IvGruVeNNibNAGLnh24c/fdOSm/oAAG6cO7euM6w5DRH7Yp2xv9+65NbuvdrwE/RZQGiBk3DSfgnOGGOMsf+T/Xrm4urPXD3RkMazRKKBZHj+O4Pe+QcAwLx582JdsuYZAGgGgFOKtb0P7e22J5u9Y1XnZXpks/lh7e2FCbpzMDUymQLfbTUiMplVIzOZjoTuHEyNcnZzvw4uRChOR6CFC+9ZeM3CxbddBkBHXHvZtaPixfh5gJBfdM+i5oWLF84HxBEtF7ek9mbbx8V6T+c5F9EQBLIhDOWndedgaiDSD3VnYGoYRtV0wxDNunMwNcrZzf16+++Fixf+BADgmrnXWCTpbEBa85Nf/uStljktkwHw/VvKEr6KppwKALk93fYGab4KiHzP+wgQQm5AjLXqzsHUEIIe1p2BqRGGQR7R2LtL1mzAKmc3y7K2CIU0DZBmAFFw1YVXDQOiUYjwbP+/I9IGADh4959rbV1TH4+bZkNDomPpUjKHDy/UBwH2NDbWr1/aV/uCCTT12eXrhk+bNnbT8uXrhtfUhEOIqt6yrIM7M5nCGNOk6k2b6gsnnYRBNpsf39cXBE1NEwqrVq2q6u6OH4IYdk6ZMv6t1tY1Q+Jxc3hfX7CpqWnCVs9bfRBibHBdnVyfSqV6ds3x4INkHH54ISGl6LXtcW943puDEHtHdXcbW6dNG7tpxYq1B4ahHGoYNW9PnnzQjt1zeJ6fCkMpGxvH53fP8cILqw6orY2PCEP5bjqd2rJ7jvb21eOIYrEpU+pzS5aA2DVHa+v62ni8eHBfX2xbU9OYd3bP0drqj47Hsaavb+Papqam4q45VqxYEQ/DwWMNQ3ZNnpx6c/ccmcyqkaYZr9uxQ2w47rhx3bvmAABsby8kiYw+yxq7bvcc2Wx+GAAMi8e7Nx5++OHbd8/hOLkkEfTZdv1Du+dYuXJlXV9fzUgA2NzQkNi8ew7PWzcWMYxPmVLvAwDsmuP559fWDB4sRwdB3/bGxkM39ufo6up759hjD922YkXu4DAUtYaxY93kyZP7HCeXFALRspK51tbWWDw+clxfH3U3NSU39OcwDLFl8uRx77a2rh8RjxcP6OuLvdnUNKarP8fKlfX55maQ7e2FFGKxOGXKxLX9OYiKOyxr4tuOkxtqGOLA3XPU1PS9ceihh/ZmMh0JwxBi1xz9+/yKFW8PDsPugz4oh+uuPUQIWbVrjv59PpfLVW/fLsbsnqN/n29v7xhFZAzaNYdpGsa/615/jv59fvfuERn3Z7P58f+ue/05Pqh7/Tn+Xff6c3xQ9/pzfFD3stn8+CAIw3/Xvf4cH9S9/hz/rnv9OT6oe/05/l33+nN8UPf6c/y77vXn+KDu9ef4oO55np+SkiidTvkfdAwAgEJDQ8L7/zkG/K/u9efYvXu759i1ewAAu+bYvXv9x4AP6t7OHPTvutef44O615/j33WvP8f/pXsqjgEf1L3+HLt3b9cc556L4f/uHv7R81YfZFkT3/533Xv11dQ7556LoYr3/f16WaFlbkvTvHnzRiy8Z+EdCxcvmgUEvWYMzySELAEd0v99EuBADGHF7j9vmmK+lHAjAMDw4YU6KeUCRLoIAOC0qu2zDkAaFo+HHwYAiMfDD5f+vfdwAAAh6GIp5YLhwwt1AABSwo2mKeYDAHR1xcZIKReEYWnhs1jMbJJSLjAM85jS98Y+LqVcsGOHcQgAgGGYV0opFxCRGD++Y3Dp9wSXAAAQ9UySUi6oqiqeDABQLIYnSikX9PX1HAEAgCg/JaVcUFe3YejOHF9AFC0AANu3GwdLKRdIiZ8AAKiujk2VUi4QAo4tbds8Q0q5YOtWWV/6WfMKKeWCJ58Eo6HhtdpSpnBO6fH2HFp6DH2nAgD09QUzpZQLgqDzqNJzCRdIKRfEYkOHlbYN1wshri1tt2qUlHJBsYifBACorTVsKeUCADyu9BjM06WUCwYPDhIAAGFoXi6lXLB69ep4W9uG6p055paey2C8lHKBaRY/Wvo9OENKuaC7u6ah9FzSeaVtjxhR2jZch4hfc92OU3t760aUMotzAQB6emqs0rbhQ6XHGP+olHJBdbXceQkt/IyUcsGyZeuqfN+Pl57L4mWl51KmSs9H/LRSZjpOSrmgpqbKKr1Oorn0mA8YWdq2uJYIrgMAqKkZPrz0s3ReaX+pPqr0mgYzSo+h7yNSygXxeDih9BiDS6WUCw477K3qtrY2s7RvmZcDAAwaVKwvPQbj9NLP4vTSz5rp0uuE50gpF/T04EE7c1wDADeU/jxkaGl/kBcAAARB5+RS5vCE0nPde3Lp33sOLeUI50gpF9TXrx60ZAmI0mMwrgAA2LJFjis97+aZO5/3aaX9wZha2gfwE6V93hxdyinmSwlfKO0Pbwwp7fPyU6XMPUeUXofwJACAqqpwVmnbPYftfD4WSikXDB3qDyYiLO0P4moAgB07jENKz4/x8dJjNI8uPQbz6NLrZHxcSrlg+3YxZucx4Oqd3cOhQ/3Bpd8TXFx6/XsOK3UvnFV6zcOTpJQLent7Di89RrpISrmgtvaNIf3dM4zSMWDHDnP0zmPA2aXn0phaeowwrbRtc7aUcsG775pjS8+HcYWUcsGSJSDq61cPKn1v+OnSdrsPLT2Xvf3HgBN2dm9y6THKC0qPcUj/nS5v2Pk6w9atYlRp3yodA+JxM1163XB66bk0TpdSLhg0qFi/s3vzpJQL2trazMMOe2tn94I5pZ8NJ+w8Bnyk9DoFM6SUC7q6qnceA+j80vN14IGlbcN1/ceAYvGAkaV9q3QJpKamyio9BnleNtsx0zDip+18rpOlfdH4jJRyge/78WXL1lWVXpfwM6XfU0yVHm985zEAPlTax2us0nON55Zep7oRO5+PzyGWutfdPWh46bnEc0t/r2koPUacsXPbHy39ezB+5z4/t/R8bKhevXp1fOe+ddnO/TJZep3M00u/Bo8r7Q+GXXrM+MnSYxzS371rieD60s8eeGBpv4X+7vUfA2aW9oe+U0vPR8/E0rbDS6WUCxoaXqt98kkwSsci84pSf/rGlR6DeUbp98Cxpf7EppZep1L3tm4Vo0r7rWjp715d3YahpddUXljKvOPInceAE0s5iyeXflf3oaWfDS6RUi4YP75jMBHtPAaYVwIAvPuuOZZILgYwZ5ceg3lM6bk0m0r7Fp516KF+HSiyXz8tcs2c+V+WBOaiexZ9c+ffXwQhLqeQBgPStxYuXnTSVRdeNUzERQYkNC26d9GmPd325cu2rvhld10VfVhM3H+PgJVDJpM7QQicZdvJr+nOwvad4+SfSqcTJ+jOwfad5xVmA8hJlpW8RXcWtu/K2c39e1nEFPdiIH/ZMqelFVAKInp64V0/85qbm41Rg0a+dc2c+S8SQD0SfOdnezGwAAB4OTCfIcBZ+ys6Kx8p8XUA0ak7B1NDCOA3oogIw8CLxcDXnYOpEbluXnvZtaNuar4pvvvXW+a2jLlx7tz/22kYXnKdMcYYYyo1P7f9H7iUVunOwfZde3tHg+Pk5ujOwdRwXf9HujMwNTIZf7rj5PmjqBFRzm6W5dMi+8NQIQ/hJdejoVjEYULgeN05mBpEOFV3BqaGYYiRAJJvohUR5exmxd6E6qGe2vt5yfVoGDYMlguBt+rOwdQgwgt0Z2BqFIvm48VieKfuHEwN7uae4DkXjDHG2IBUsWcuLq3eeimWYeE1tv85Tv54z/O/rDsHU8N1fV5yPSKyWf801/V5yfWIKGc3K3ZwwRhjjDGmFl8WYYwxxgakij1zkRDBiP/8XawSrFy5si6TKYzRnYOp4Xkdh+nOwNRwnNzQ9vaOUbpzMDXK2c2KHVycWtV1Di+5Hg1dXVWNQsgrdedgakhp3KE7A1NDCGOmlOJi3TmYGuXsZsW+OW+R4g1AJN052L6LxWgzEXXozsHUQKQ23RmYGmEoN0qJed05mBrczT3Bcy4YY4yxAaliz1xMj3UdCcQfRY2C1lZ/dCZTaNKdg6nhefkzdWdgamQyHYn29o4G3TmYGuXsZsUOLo4ygxkIfPvvKBCCJgkhZ+vOwdSQEm7QnYGpYRimJaU4VXcOpkY5u1mxa4u8HsSXA+IM3TnYvovFwg4p43wr98igX+hOwNQIAnzFMOSbunMwVbib/xnPuWCMMcYGpIq9LHJGVddMvv13NLjumkmeV+DLIhHheT5fFomITKbDcpz8KbpzMDXK2c2KHVyMFsERvOR6NEgpRhNJntAZEVIiT+iMCMMwE0KQpTsHU6Oc3azYORfPF6seJcSzdOdg+840RRYR1+nOwdQgwut1Z2BqhGHvMgDD052DqcHd3BM854IxxhgbkCr2ssiF1dvP4TkX0eB5fqPr5vj23xHhOP6dujMwNbLZjpmum79Idw6mRjm7WbGDi0FII3jORTSEIdUBIC9cFhk4SXcCpgaRORSRDtadg6lSvm5W7ODigd7BiwmR740QAVu3Jp/r6jK+qzsHU6OmpvsM3RmYGsXi248dcAAt1J2DqVHOblbs4IIxxhhjA1PFDi7Or9oxh5dcj4YhQ/zja2vDL+nOwdTo7q55RHcGpkYsdtBp27Zhi+4cTI1ydrNi35x7CLbzkuvRIAR1A9Bm3TmYGoi0XncGpgp1AsAW3SmYGtzNPcEfRWWMMcYGpIo9czHZ7B2rOwNTI5vND2tvL0zQnYOpkckU+G6rEZHJrBqZyXQkdOdgapSzmxU7uDgu1ns6z7mIhiCQDWEoP607B1MDkX6oOwNTwzCqphuGaNadg6lRzm5W7JvzBmm+ynMuokEIuQFRtI362qoAACAASURBVOrOwdQQgh7WnYGpEYZBXkrk239HBHdzT/CcC8YYY2xAqtgzFyfEemy+/Xc0tLevHud5hRm6czA1HMe/UHcGpkZb29qJrttxjO4cTI1ydrNiBxeTzL5pfPvvaCgWjfFE8hTdOZgqeIXuBEwN06QjEcVM3TmYKuXrZsUuuf5yYD5DgLN052D7Tkp8HUB06s7B1BACbtGdgakRhoEXi4GvOwdTg7u5J3jOBWOMMTYgVexlkeaqHWfwnItoaG/vaHCc3BzdOZgaruv/SHcGpkYm4093nDx/FDUiytnNir0sMlTIQ3jORTQUizhMCByvOwdTgwin6s7A1DAMMRJA8k20IqKc3azYMxf/6K39Ay+5Hg21tb0ZKcXtunMwNYQI5+nOwNSQMnxaCHmf7hxMDe7mnuA5F4wxxtiAVLFnLi6t3nopz7mIBsfJH+95/pd152BquK7PS65HRDbrn+a6Pi+5HhHl7GbFDi4YY4wxxtTiyyKMMcbYgFSxZy5GYXCA7gxMjVwuV53N5ofpzsHUyGQKY3RnYGq0tq6vbW1dM0R3DqZGObtZsYOL2dVdF/CS69GweTNMk5Ku052DqYFI9+vOwNSIxYKTYzHjct05mBrl7GbFvjlvkeINXnI9GmIx2kxEHbpzMDUQqU13BqZGGMqNUmJedw6mBndzT/CcC8YYY2xAqtgzF9NjXUcC8UdRo6C11R+dyRSadOdganhe/kzdGZgamUxHor29o0F3DqZGObtZsYOLo8xgBgLf/jsKhKBJQsjZunMwNaSEG3RnYGoYhmlJKU7VnYOpUc5uVuzgIh8aLs+5iAbDwAKAfFZ3DqYGovyt7gxMjSCQqxDli7pzMDW4m3uC51wwxhhjA1LFnrk4o6prJt/+Oxpcd80kzyvwZZGI8DyfL4tERCbTYTlO/hTdOZga5exmxQ4uRovgCF5yPRqkFKOJJE/ojAgpkSd0RoRhmAkhyNKdg6lRzm6a5fpFqj1frHqUEM/SnYPtO9MUWURcpzsHU4MIr9edgakRhr3LAAxPdw6mBndzT/CcC8YYY2xAqtjLIhdWbz+H51xEg+f5ja6bu1J3DqaG4/h36s7A1MhmO2a6bv4i3TmYGuXsZsUOLgYhjeA5F9EQhlQHgLzYVWTgJN0JmBpE5lBEOlh3DqZK+bpZsYOLB3oHLyZEqTsH23dbtyaf6+oyvqs7B1Ojpqb7DN0ZmBrF4tuPHXAALdSdg6lRzm5W7OCCMcYYYwNTxQ4uzq/aMYeXXI+GIUP842trwy/pzsHU6O6ueUR3BqZGLHbQadu2YYvuHEyNcnazYt+cewi28+2/o0EI6gagzbpzMDUQab3uDEwV6gSALbpTMDXK2c39PiHy2suuHSWD4mkScXMxVvzrHXfcUQQAuHHu3LrOsOY0ROyLdcb+fuuSW7v3asNP0GcBoQVO4sljjDHG2ECyX89czL9k/vAwCJYTwEcBYEa8L7b6xrlz6+bNmxfrkjXPAEAzAJxSrO19aG+3PdnsHas8MNMim80Pa28vTNCdg6mRyRT4bqsRkcmsGpnJdCR052BqlLOb+3VwIQw6iQD+tXDxbecvunvRDQDwUndYfXa8GD8PEPKL7lnUvHDxwvmAOKLl4pbU3mz7uFjv6TznIhqCQDaEofy07hxMDUT6oe4MTA3DqJpuGKJZdw6mRjm7uV/fnIuBfNo0zS8DAFxzzTVVBHBUQOQCwWQgfP+WsoSvoimn7s22N0jzVZ5zEQ1CyA2IolV3DqaGEPSw7gxMjTAM8lIi3/47IsrZzf26tsjt993+NgDANXOuOpZ2hLchiAd/fu/Psy2fnv85RHi2//sQaQMA/K8btTiO/ylEHGTbiTtWrHh7cLHYfSGRfCudTv3l5TC+YjDKGU9l/OmNjcllmYw/XQicIgQ83tCQ6HCc3McRxahYrOa3kycftMN18/OIqDOdTv7G81YfRBQ7iwhy6XTin9ns2ilSyumGEb4wZcr4rOflZxHBhFhMPjR5cupNx/EvRMTBllV/Z1vbhhrTLF5EJN9Op1N/9jw/RYSnENHL6XTyec/LTSMSlmHgv6ZMqV/jeYXZRHRwdXX3/Ycffvh2z/MvI8Je207cl8msGilE/GyA0Lft8f9w3fxkADieSL6YTqdc1/VPAsBDg4AebmpKbvC83PlE4oDXX6//1dFH+7GtW8UlAPSObSf/mMl0JIQwPgIAr9h24lnHKRyNSGmiYGk6PWGV6+bPAIAxQRD8rqlpwlbX9eciYmBZiXuXL183vKoqPAdRFiwr9be2Nv8Iw8AZAPSSbScdx8mdiCgmAYhHbXvcG66bPxcAhm7eXH/X2LGrjc7O+KcR4V3LSvy+vX31uDCMnYaIr1pW/TOZTKFJCGpEDJ+yrPGvOY5/OiIeQiQfTKdTW1w3dykRkG3XL16xYu2BxaL8pBC0tqEh+Vg2mztcSjFTCGhraEi0ZbMdM6U0Dkc0/mpZY9c5Tr4ZEYbFYjsWd3d3k2mOvJQINqfTiSWet24sUfgxKeVrjY2ppzzPbyTCpjCkZ6ZOTb6azfqnSYnjYjHx+8mTx73rOLlPI6Jh28m7HCc3FFGcS0RvpNPJR113zSQA80QidNLp+pc8rzCDiI4IguBvTU0TCp7nn0OEw4cMkff6fjIYNqwwF0Bute3U7zKZwhgh6AwAWmXbyaWu66cB8Ggp4dnGxsQrrpv7CIBIBEHsj01NY97xvPwlABSzrOSvXnhh1QHV1fHzEWGDZSUedpw1hyKaJwGErm2Pf9Fx8scjwmQpw783No7Pu67/CQAc0dUl7ps+fWyv5xUuQ5TbLSt1f2urP9o08UxEWm1ZyScymQ5LCGMaIj5vWfUvu27HqQBGUsq+PzU2HrrRdfMXEUFNOp24s797iPimZdU/1N5emBCGNAtRZi0r9UJ/9xDpn5aVzElJq103Pw+x6jeWdXDnzu7tSKeTv21v7xgVhsbHpYSOxsbE4/3dE0Isa2gY157J5E8WAsYbRviXKVPGv9XfPdtO3OF5bw4i6v3U7t2TktobG5PLPC93LJFo2L17/ccAx8lfjgjdtp349e7HAM8rHEVEx0kZLm9sHO/tfgzwvNwFRKLOsup/uWzZuqraWnkxAG207eSfHCeXRBSnEsGKdDrxnOt2HANg2P3d87z8mUQwuqen74Fjjz10m+f5nwHAomUl7m1tXT/CNIufAJB52079PZPJHykEfKi/e/3HACnxkcbG+vWumzsPQAzZvLn+rmTSN7duFZcg0ibLSv6htXVNvWmaH+3vXv8xACB40rYnvP6/u+fPJaIwnU7ds3v3+o8BiNRqWYl/ZjK5E4QQh+1+DAiCjXfX1NRgsTh4zu7HACHClQ0N45/OZvNTpYSpQsinGxpSKx3H/xgijhUCljQ0JDY7Tm4OIqBtp+7OZvPDpIRmIlqXTif/6nkdhxEZJ0iJmcbG+tb+7hlG8bEpUyau9bz8J4ngwEGD+u5Zt25iWOoebLHtxIOuu/YQAHk6kXw9nU492d+9/mOA5+U+SiTqe3uNP0ybNnaT5+UvISLTtpN3tbauGWKa5nkAsN62E4/0d6//GOC6+Q8BwJG7d2/IEHnvSy8li5MmFT6DKLdZVuqB3Y8BjpOzEcUxAPCcbSdW7H4McN38xYhUZVnJX65cubKup6fmgv5jQFvb2omGIT+MKD3LSi13HP84RDyKSP4jnU75ruufDYAjgyD266lTR3d7XuHy/u6tWJE7uFjEbZ6Xn2VZiX+1t3c0hKFxbH/3HCd/ShjGnmtqGtOl4v1/v19WaJlz9dUE4seAYt7CxQu/BgBACFkCOqT/eyTAgRjCir3Z7mSj70jVWZkeiPJAzyvM0J2DKTNTdwCmhpRiXGmwxKKBovFazp8z/5T5c+Y/19zcbOz69ZZPtxzfMmf+UgCAqy68atj8OfNz8y+ZP3xvtn35sq0r8Am5WmVepkcmkzvBdf1v6s7B1HCc/FO6MzA1PK8w2/P8G3TnYGqUs5v7d8l1hI8gwRGjBh2Ua5kzv/QlpG+/ueOtX44aNPKta+bMf5EA6pHgOz+7d9Gmvdn0y4H5DAHO2i+5WVlJia8DiE7dOZgaQsAtujMwNcIw8GIx8HXnYGr8P9PNlrktY26cO7fu//TDvOQ6Y4wxxlQ6+/kdD+NSWqU7B9t3mUz+yNJENRYFrpv7hu4MTI1MptDkefkzdedgapSzmxV7n4gRGI7nJdejQo4EwMm6UzA1iMSJujMwNQwDxgDQYbpzMDXK2c2KHVz8o7f2D7zkejTU1vZmpBS3687B1BAinKc7A1NDyvBpIeR9unMwNbibe4LnXDDGGGMDUsWeubi0euulWIaF19j+5zj54z3P/7LuHEwN1/V5yfWIyGb901zX5yXXI6Kc3azYwYWJGOc5F9FAJE0iqNadg6lBhP+3T4CxAYdIxBC5m1HB3dwTfFmEMcYYG5Aq9szFKAwO0J2BqZHL5aqz2fww3TmYGplMYYzuDEyN1tb1ta2ta4bozsHUKGc3K3ZwMbu66wJecj0aNm+GaVLSdbpzMDUQ6X7dGZgasVhwcixmXK47B1OjnN2s2DfnLVK8wUuuR0MsRpuJqEN3DqYGIrXpzsDUCEO5UUrM687B1OBu7gmec8EYY4wNSBV75iJt9E4A4o+iRkEms2qk6+b5Dp0R4Ti5E3VnYGpkMoUxntfBd+iMiHJ2s2IHF03x3pMR+KOo0WAeCUC8tkhkCF5bJCIMA5oABK8tEhnl62bFDi7yoeHynItoMAwsAMhndedgaiDK3+rOwNQIArkKUb6oOwdTg7u5J3jOBWOMMTYgVeyZizOqumby7b+jwXXXTPK8wmzdOZganuffoDsDUyOT6bAcJ3+K7hxMjXJ2s2IHF6NFcATf/jsapBSjiWST7hxMDSmRr9FHhGGYCSHI0p2DqVHObprl+kWqtfZVPU6AH9Odg6kQvAIQf0d3CqaK/LruBEyNMIRWw5Cv6c7BVOFu/mc854IxxhgbkCr2ssiF1dvP4TkX0eB5fqPr5q7UnYOp4Tj+nbozMDWy2Y6Zrpu/SHcOpkY5u1mxg4tBSCN4zkU0hCHVASAvdhUZOEl3AqYGkTkUkQ7WnYOpUr5uVuzg4oHewYsJUerOwfbd1q3J57q6jO/qzsHUqKnpPkN3BqZGsfj2YwccQAt152BqcDf3BM+5YIwxxgakij1zcXnNtit4yfVoyGRyJ7iu/03dOZgajpN/SncGpobnFWbzfUuio5zdrNg35x6C7Xz772gQgroBaLPuHEwNRFqvOwNThToBYIvuFEwN7uae4MsijDHG2IBUsWcuJpu9Y3VnYGpks/lh7e2FCbpzMDUymQLfbTUiMplVIzOZjoTuHEyNcnazYgcXx8V6T+c5F9EQBLIhDOWndedgaiDSD3VnYGoYRtV0wxDNunMwNcrZzYp9c36HjA6ecxEVYiMArdCdgqmBKJ/UnYGpEYawHgD59t8Rwd3cEzzngjHGGBuQKvbMxQmxHptv/x0N7e2rx3leYYbuHEwNx/Ev1J2BqdHWtnai63YcozsHU6Oc3azYwcUks28a3/47GopFYzyRPEV3DqYKXqE7AVPDNOlIRDFTdw6mSvm6WbFLrr8cmM8Q4CzdOdi+kxJfBxCdunMwNYSAW3RnYGqEYeDFYuDrzsHU4G7uCZ5zwRhjjA1IFXtZ5OzqzlN4zkU0ZDL5I103d57uHEwN1819Q3cGpkYmU2jyvPyZunMwNcrZzYodXIzAcDzPuYgKORIAJ+tOwdQgEifqzsDUMAwYA0CH6c7B1ChnNyt2cPGP3to/8JLr0VBb25uRUtyuOwdTQ4hwnu4MTA0pw6eFkPfpzsHU4G7uCZ5zwRhjjA1IFXvm4tLqrZfynItocJz88Z7nf1l3DqaG6/qP6M7A1Mhm/dNc12/RnYOpUc5uVuzgwkSM85yLaCCSJhFU687B1CDCOt0ZmBpEIobI3YwK7uae4MsijDHG2IBUsWcuRmFwgO4MTI1cLledzeaH6c7B1MhkCmN0Z2BqtLaur21tXTNEdw6mRjm7WbGDi9nVXRfwkuvRsHkzTJOSrtOdg6mBSPfrzsDUiMWCk2Mx43LdOZga5exmxb45b5HiDV5yPRpiMdpMRB26czA1EKlNdwamRhjKjVJiXncOpgZ3c0/wnAvGGGNsQKrYMxdpo3cCEH8UNQoymVUjXTfPd+iMCMfJnag7A1MjkymM8bwOvkNnRJSzmxU7uGiK956MwB9FjQbzSADitUUiQ/DaIhFhGNAEIHhtkcgoXzcrdnCRDw2X51xEg2FgAUA+qzsHUwNR/lZ3BqZGEMhViPJF3TmYGtzNPcFzLhhjjLEBqWLPXJwa75zGt/+OBs/zU67bcaruHEwN181doTsDU6OtzT8im+2YqTsHU6Oc3azYwUXCCG2+/Xc0hCHVA4gP6c7B1CASF+rOwNQwTXEokThGdw6mRjm7aZbrF6nW2lf1OAF+THcOpkLwCkD8Hd0pmCry67oTMDXCEFoNQ76mOwdThbv5n/GcC8YYY2xAqtjLIhdWbz+H51xEg+f5ja6bu1J3DqaG4/h36s7A1MhmO2a6bv4i3TmYGuXsZsUOLgYhjeA5F9EQhlQHgLzYVWTgJN0JmBpE5lBEOlh3DqZK+bpZsYOLh3pq7ydEqTsH23fDhsFyIfBW3TmYGkR4ge4MTI1i0Xy8WAz5TFREcDf3BM+5YIwxxgakij1zcXnNtit4yfVoyGRyJ7iu/03dOZgajpN/SncGpobnFWZ7nn+D7hxMjXJ2sywfRb2u+bqa4qDiGQsXL1zS/7Ub586t6wxrTkPEvlhn7O+3Lrm1e2+22UOwHRDj6tOychOCugFgs+4cTA1EWq87A1OFOgFgi+4UTI1ydnO/T4i8+tKrxwmC6wHxmIV3LzoOAGDevHmxeF/sJSJYhYhvA8lJC++57ZS92vAT9FlAaIGTePIYY4wxNpDs98sKQuJtAGLarl+LF+PnAUJ+0T2LmhcuXjgfEEe0XNyS2pvtJkQwQm1SpsvKlSvrMpkCf1okIniJ7uhwnNzQ9vaOUbpzMDXK2c39PrhYeM+iM0nS5//HFwkmA6H3/t/xVTTl1L3Z7qlVXefwnIto6OqqahRC8n0uIkJK4w7dGZgaQhgzpRQX687B1ChnN/W8ORONQiC//6+ItAEA9uqz1O+Q0cFLrkeF2AhAK3SnYGogyid1Z2BqhCGsB0C+/XdElLObWtYWIYQsAB3S/3cJcKARwiO7f5/r5m4GEENsO3FlNpsfJqX8HpFYk04nvr8qiK08EMMZjzv55nQ6scRx8s2I8mQhxB0NDYk2180vAJDjhRBfbGhIbHbd/O0AcqttpxZ4np8ioi8CQMa2U7/wvPwsInkukfh9Op34p+PkL0eUTYZhfH/KlPo1jpP/DqIcblnJq9raOupMU3wfEXzLSn3Xdf00AF0JAE/Ydup3rut/AoA+QmT8Mp2uf8lx/C8g0sQgqPpKU9OYdxzHvw0AOtPp5BcymY6EEPhlRHAtK/Vz1/VPAqDzEeFPlpX6m+v6cwFoGoD8oW1PeN3z/G8R0cjNm5PzR49+raanp+oWIiyk08lvZzIdlhB4NRE+lU4nf+s4ubMQ4TREuNuyUi+4rv95ADoMMfwvy5r4tuv6CwGgz7aTn29vXz0uDI2vAmC7bScXZrMdM6XETxHhX9Lp5F8dJzcHEaYLAbc2NKRWlj7ZQaMGDQo+u3XrIMM0e28FEG/YduKbnlc4iii8BkA8Y9uJX3teYTZReDqRuDedTjznOLnPIcIRQYA3NTUlN7hu7idEkmw79blMpjBGiPDriGKFZSV+6rr5DwHIixHFI5aVeNh18xcDyA9JKX7S2Jh4xfPyXyeSY7q6zM+NHh3Q1q3wEyLYkE6nbspk8kcKIa9FFM9ZVuJex/FPR6TZiMavLav+GdfNXwMgjwIwv2nb497wvPytRGTYdvKzK1bkDi4W4RtE8Go6nfpxJuNPF4LmEMFf0+nUXxzH/xQizRTCXNjQMK7ddf3/AqCxiDWfLxbX9Znm8IUA+JZtJ7/meR2HEeHnAeAF207dnc36p0lJZ0kJv21sTD3lOP58RGowjPBbU6ZMXOu6/i2IUGNZyfmZzKqRQpjfAsDXbTv5Q8/LTSOCuQD4N9tO/snzcucTwUlS0m2NjeM9x/G/gkj1PT3BFwqFiZ2TJvm3EeE76XTyK46z5lBEcQMivmhZyV+5bu4jAPAJRPydZSWf8LzcVURgE8F30+mU73n+94mwzrYTVy1fvm54VVXxO0S4Op1O/sBxCkcjhpchin9aVuL3rps/F0DOQsRfWFYyAwDdrpv7BREssO3kVs/zbweAzbad+mJ7e2FCGIY3Aog2207ckcnkTxZCNkspljQ2Jh533fw8ADlVCHFzQ0Oiw3Vz3wOAYZaVvNJ1/SGIcDMA5Gw79T3P8xuJ6AoA8S/bTjzoeflPEslTpDTubGysb3Wc/I2IckJvb+zL06aN3eS6+Z8j0nbLSt7oOLkkInwJABzbTt3uef6Hieg8APijbaf+7nn+ZUR0dBiaP5g6ddxqx/G/jUgjXn89eXV9/epB1dXmDxAxb1nJ7zhOzkaEqxBhqWWlHnBd/2wA+igi3GVZqeWu618PQJOkDL7a2HjoRs/zFxFBt20nb2htXVNvmuIrROil08nbMpncCULAhULgnxsako+5bu5SADgWkX5kWeNfc13/vwHooCDY1BKLjY0Tdf8IANfadvJb2ezaKVIGLUT4dDqd/I3j5D6OCB+TEhc3NiaXuW7uOgA43DDoa1OmjH/Ldf2fImJoWYnrXHftIQDB1wDEy7ad+JnnFWYQhRcRyYcsK/mw5+UvIZLHhyH+eOrU5KuOk7sJEUYPGQLXbthgYm1t8GNEsd6yEt9w3fxkAPlZAPGsbSfu87z8mUTyDABxn20nnvW8/LVE8kgpjW80Ntavd5zcjxEBbTt1bWurP9o06SZE8YplJX7iOPnjEeUliMajllX/UOluoXIGovEzy6p/2XXzXwOQhwRB1XVDhnSGnZ3mTwHgTdtOfb2tzT/CMOhzRLAsnU4tdhz/Y4j0cSHoNw0N4592Xb8FgKYgxv7bssauc13/RwAQt+1kS3t7x6gwxG8C4Gu2nfyR5+WOJYJLieCxdDr1Z8fxL0SkEwyDFk2ZMj7ruv5XAWhcdXXvDRs2HNY9bJi/CBE3Wlbyq667ZhKAuB4Al9t28i7Py32UCM4GwAdsO7nUcfyrEckKAvntpqYJBcfxfwAAg9Lp5NWtretHmGbvtxFxlWUlb3HdjmMA8DMA+HfbTv7RdXPnAcCHieDn6XTK9bzcl4ggGQTyxqlTx2/3PP/nRGJTOp34clvb2okAwcGe519mWclfOk7+FET5SUTxoGUl/uW6uSuI4HfpdErJBF4tgwskfAkQvgUA/33VhVcNA4ATQ8D/9XGnIJCL4nFhAgBs2lS/ffjwws1BgD0AAAcbctRrQWxLX5/xBABAX5/xRE0NtBFVvQUAICXeZ5qietOm+u0AAELA9/v6ZAAAUFtbXN/dHb8ZMewEACgWg9Z43Mz19QWbSt9b/Ati7F+DB4frAQDCMLg9HjdNRJQPPkg7Dj+8cLOUohcAALH6dcTem7u7ja0AALGY8WQYStcwqt8GACASvzEMqt6+ffSWnTl+EIZSAgDU1YVv7pqjp6fYVlsbz4ehfLe07eARxNiTdXWwvvSzwS+IYrETT4RwyZLDunbNIWX1KtMs3tzXF9sGABCPm0+HocwaRs3bpecS7o/HRU1f36bNpW3DD/tzCNH7FtHgmw1DdgEAdHWFbm1tfF1/DqLgUcOIP71jh9gAAGAYwZ1EsdjEiRP7AADb2ws3Exl9pefS7IjH6b0ciPQMong5Hu/eWHou8XfxONYAbHyntG241TBwtOcVZpjm9uVh+H6O6upur6+vZj3s/CSJlH1/M834c11dpRwAxq+EgPj06WN7AQB2zdHTI3KDB8PNQdC3vZQZnwfAV7u6et8pvU5ySRiKRw1j28bStuVPhEAEAOju3rQpHh95c18fdZf2l56X+/pqbjYMsaX0GOJ/j8eLy/r6jDdLj9G8GzGMr1w5qqe5eZRsby/cjFgsAgB0dsYKgwfLm4mKO0o/S8sMQ7zW09P3Tul1oj+EoXisujp4e2eOnxmGEKU/b91imiPf2+dNc9CKMOx+L0exWPV4PF58MQhiO3MYi4WQVYXCxM7mZpDt7eLmvr4gAAAYOlSs3b5dvJeDCJYbhlj1/j5PfyQy/l5T07ehlFMuMk3DKO0Ph2ytqXm/e/F49ath2P3ePt/ba/yrpgZa+7tHZGw3DLp506b6HYhI2Wz+vRyDB4dvbN8u3tvnpQxeMk2zIwhKOQwj/AuR8XhdnVy/8xhw287u0dKltGP48Pf3eYDq14R4v3s9PcbSmhrIxGL93cNfG4ao7uo6ZGt/94JAhqUcwYb/eQwI2+Jx0+/f5wGCh4SIPTFkSLDzGBD+Ih43zeZmkEuWTOzctXtC1KxC7H1vn4/FjKfCUHr93ZNS3G+aVB0EW/sP2rdIWerekCHyre5u8719vq8vcGpr42vfPwaEjyLGntqle3cQxWJTp04NliwBefjh7+/zfX3Gmnhc7noMeCYMZXt1dal7QYAPxONY09u7cee24VYpJZUyb9u4a/e6u3u92tr4G8UiDHPdjmPCsO8x04w/29890wx/RRSLJZP1fcnk/+xeEMRy8Xjxve4hwrOIYkX/MaBYpAfjcVETi23vX6zwx7TzHHRNTeemMBz8XvdqarqzfX01NwPQ5p3b/ls8Xny+r8/s3+fvQoT41KmjewCAdu1eb6/hDx4s38sBQM8LIVZ2dfXtPAbQ78NQ/NUwtvZ37yf93evtfffdXY8Bpjno5V27F4bxf8TjxRd6e0vdAzDuFkJWZbOHdfV3rz9HV1d87a7HtzeoCAAAGiBJREFUACnhBcMQr/fniMfpj2Eo/jZkSLCzP3Jhf47t20dv2fV9LxYb/MquOXp7Y4/X1IQvvd89814hZFVHR3JHUxPKXbt34IHBui1b4FXE4KHSYwhejMfNNf3HAMOgP7/6amrnc7XvynL77PmXzJ+OBvyw/9Mizc3NxqhBI3+DgOMJoB4JvvOzexb9dG+2efmyrSt+2V1XRR8WE/dPalYupf+t4SzbTn5Ndxa27xwn/1Q6nThBdw627zyvMBtATrKs5C26s7B9V85uluXMxaJ7Fy0DgOP6/75kyZIQAM5vmdsypha6tn//rrv2erT0ehBfDogzVOZkesRiYYeUcb6Ve2TQL3QnYGoEAb5iGPJN3TmYKtzN/4xv/80YY4wNSBX7Uc6zqztP4SXXoyGTyR+5c2ISiwDXzX1DdwamRiZTaPK8/Jm6czA1ytnNih1cjMBwPC+5HhVyJABO1p2CqUEkTtSdgalhGDAGgPimaBFRzm5W7ODiH721f+Al16OhtrY3I6W4XXcOpoYQ4TzdGZgaUoZPCyHv052DqcHd3BM854IxxhgbkCr2zMXF1dsu5DkX0eC6HcfsvMEPiwDX9e/XnYGp4Xn5WZ7nX6Y7B1OjnN2s2MFFNUIdz7mIBimxBgCH6c7B1CBCXoQuMnAQAAzVnYKpwd3cE3xZ5P9r787jpKjOvYH/zqnumWFwVDZRok4PavK6ZLp6QIlREa8RRcWrUXDXERXFGUJc4k30atDsm8vr4IJeBfVKriS57rsO7hqhq2rmoig63T0aFIGg4AzMdNd57h9GE/PJ+8pS9KErv+9f9DJVv+HzeaZPnzr1HCIioq1Sxc5cDFelbW1noGjkcrmajo4CZy5iIpvt5rejmFi4cFntwoXvbGc7B0WjnLVZsYOLY2p6T+aW6/GwejXGGCNccxETSgnXXMREMln6VjLpnGs7B0WjnLVZsR/OPaJWcsv1eHActRaQZbZzUFR4uTIulCp9JKLY/js2WJtfjmsuiIiItkoVO3ORcfp2g/BW1DjIZpcO8/0CO3TGhOflxtnOQNHIZrtHBEEXO3TGRDlrs2IHF6Or+r6lwFtR4yGxFyDcWyQ2NPcWiQnHwWhAc2+R2ChfbVbs4KIQOj7XXMSD46huwLxgOwdFQylzj+0MFI1SySxVyvzRdg6KBmtzQ3DNBRER0VapYmcuxlf1jGH773gIgnyD73eNt52DouH7ufNsZ6BoLFqU37Ojo2us7RwUjXLWZsUOLuqd0GX773gIQ9kV0AfazkHRENGn2M5A0Ugk9B4iej/bOSga5azNRLlOFLWF/dVPCdSRtnNQFEqvA1UrbaegqJgf2k5A0QhDLHQc86btHBQV1uaX45oLIiKirVLFXhY5pWbt8VxzEQ9BkG/y/dz5tnNQNDwvf6vtDBSNjo6usb5fOM12DopGOWuzYgcXA5UM5ZqLeAhDqQO4FXB8qK/aTkDREElsr5TsaDsHRaV8tVmxg4sH1tfOE6WM7Ry0+QYNwqtaq2tt56BoiKiTbWegaBSLiaeKxZAzUTHB2twQXHNBRES0VarYmYtzB6w5j1uux0M2mzvY9/NX285B0fC8wrO2M1A0gqD7mCDIX2I7B0WjnLVZsR/O6wVr2f47HrSWdYCstp2DoqGULLOdgaIiPQA+sp2CosHa3BC8LEJERLRVqtiZi3pdGmo7A0VjyZIlddlsN+8WiQlu0R0fnpfbvrOza7jtHBSNctZmxQ4uxlf3Hs81F/HQ21vdpLVhn4uYMMaZbTsDRUNrZ6wx+nTbOSga5azNiv1wXilOF9dcxIVeAchi2ykoGkqZBbYzUDTCEMsAxfbfMcHa3BBcc0FERLRVqtiZi4OT6122/46Hzs63dwmC7oNs56BoeF6eu6LGxKJF7+7u+13cFTUmylmbFTu4+Gqifwzbf8dDseiMFDGH2c5BUVHn2U5A0UgkZC+l9FjbOSgq5avNit1y/a1S1atQit92YyCZDLuMqWIr99iQW2wnoGiUSup1xzEf2M5BUWFtfjmuuSAiItoqVexlkeNqeg7jmot4yGYLe/l+7kTbOSgavp+7ynYGikY22z06CAoTbeegaJSzNit2cDFUhSO55iIuzDBA7W07BUVDRI+znYGi4TgYAQibosVEOWuzYgcXLxWrH458y/UFci6ekW9HeswoPCNHoV1abMfYUhIJ3eE4eq7tHBQNEXWx7QwUjTDsezkMzXzbOSgarM0NsSXWXCyQh9EuN0R6zCi0y8/RLgtsxyAiItoQFTtzcXrNmlO45iIefL9rP9/PXWg7B0XD9/PzbGegaARB4dAgyJ9jOwdFo5y1WbGDixqFOq65iAdj1ABADbKdg6IhorgJXWyogQC2t52CosHa3BC8LEJERLRVqtiZizptqm1noGi0t0vipZfeHWA7B0VjyZIldbYzUDQWLlyYzOVyNbZzUDTKWZsVO7g4qfqTZm65Hg/bbZc/oLY2/IHtHBSNdesGPGQ7A0Ujmdxhwpo1qtV2DopGOWuzYj+ce0St5Jbr/8ACuQLt8iPbMTaG46i1gCyznYOiws65caFU6SMRxfbfsVG+2qzYvUXuWV/3eyhwRP33BLsBqKhpzHQ6lQWQtZ2DopHJpM61nYGi0dg48jnbGSg65azNip25yDh9u0F4K2ocZLNLh/l+gR06Y8LzcuNsZ6BoZLPdI4Kgix06Y6KctVmxg4vRVX3fUuCtqPGQ2AsQ7i0SG5p7i8SE42A0oLm3SGyUrzYrdnDxvkm8wTUX8aC1eV8pvdB2DoqG1vKg7QwUjTAsFYxRge0cFI1y1mbFrrl4qK/2OSg02s5Bm891d3sLABcBxkQ6nfq17QwUjaamkRxYxEg5a7NiZy7GV/WM+UL7703p1skOn1uFIMg3+H7XeNs5KBq+nzvPdgaKxqJF+T07OrrG2s5B0ShnbVbs4KLeCV2I1GGBfIx2+T9YAIN2SaFd+vGMHPOlB3hMBmMBDJ6W/Tb65O3yHtply/Tbb5cS2uXof/jaM3IK2uXPG3GsU7FAVv1/Xn8M7XLdRmeMWBjKroA+0HYOioaIPsV2BopGIqH3ENEb/zeStkrlrM2KvSyysL/6KYGabDsHRaH0OlC10nYKior5oe0EFI0wxELHMW/azkFRKV9tVuzgwgur37GdgaLR1LTHCgArbOegaGQyDQtsZ6BoNDXtyuZ2MVLO2qzYyyKTqj85WnHFRCx0dnY1el6u2XYOiobv56+xnYGikc3m9/e8wiTbOSga5azNih1cbK/NV7ggMx6KRTVIKTXSdg6KhogaZTsDRcNx9DCtpd52DopGOWuzYgcXD6yvnSdgn4s4GDQIr2qtrrWdg6Ihok62nYGiUSwmnioWw1tt56BolLM2K3bNxXJJrAEACBQMMtAABKP+cnPqbnhSxiEBQKOAfoRwUId/UYvxiGyLKoyC+pv9N56RUSji3S+c4DnZCYIdcbDyvvB8u3wDnw3KnpU9AaxHP0IkMBLAR1iHLlRhFDT6sRKvYQeMgoMl6MMuUFiFQ9Wf/uZYQ2EwEgPQiXXIYBX+2kjqKWmEg8EA8n/3qys8JYdAIUQIBcFCHK56Ps9jsA2AEp6TnVBCfdkapLdLAgb7AaiBRhaHqI/+8vzu0FA4WC0FADwp20FjbxyqXvrsRxsaGtYDWL8R5xoNBwWMVdGt01goSXyMfdELHxNVb2TH/SfE6/TxMXr0CNZCjJSzNit25uLcAWvOUxAFBQca9wAAFH6HT3tfXI4E5gC4GQYXw8H5ULgJADAATXDwDDQe+fxgCvehGv/6hROEmAzBb//BqZ8HUAUAMPgNQlyIBKYBuA3Ab1GLUXDwDBRewCDsCMHzKOEAJHAHNJq/cCSFw6DxONZjJDRexHAM//w1B3OgcAsULvy78zt/yf8skliAKuz5+SsG10BhLwBACScDaNmQ/8tI9GMwNF6ExtNQ2P/z5wVXIsRfW84msR80XsS94nz2VDabO9j381dv8LkUHkaIIyNK/qleDIXGixiIPSI97j8hzys8azsDRSMIuo8JgvwltnNQNMpZm9ZmLi6dMqWuJxwwQSnVn+xJPn7t/GvXbczPl0T62f47HpTSJaVkw2cuaKumlKy1nYGioZQpGrMRs4q0VStnbVoZXEydOjXZ2598HsBSAB8Wa/taABy2Mce4Y/12dwD4AfdFrXyZTP2LAF60nYOi4bqpf9wEjipOY2PqUdsZKDrlrE0rl0WqilUnQqEwa+6sSW1z2lqg1NDW01sbNuYY9bo0dEvlo/JasmRJXTbbPcJ2DooGt+iOD8/Lbd/Z2TX8y99JlaCctWlnzYVgb8jf7LQn6g2VMBt1i8z46t7jueV6PPT2Vjdpbc63nYOiYYwz23YGiobWzlhj9Om2c1A0ylmblgYXMlxB8p89VEreB7DjxhxipThdXHMRF3oFIIttp6BoKGUW2M5A0QhDLAMU23/HRDlr08qaC1HoAOQrnz02wGAnxEN//z7fz88DZIjrNozPZpcOc5zkPSLyuus2zKhV4RMTkj3p9yTxWkep+o2mquLXdlD9e7xeql7UbRLd4xK9B22nZVBgql7Ll5IfXlC9+ripfn6euwr/PiwZ/nbf5PoDVhnn/VdLNe+P1utfOLvm43Pv66tb9Hix9mnPy1/+fGnFsdet2/7hrk9z3K0Uhjc21h9e96y56aDkumMPq+o9/KJPhv3hiJqeYWMTvZP/0LdN98LigMfatlmx//19tX96xyRf7kLVmkOTvd6UAWsuu7Bnh/s+LOlXfD//qFJYm06nJh9btXZdElh9Ss0nZx23Zqdrfle77MC7++ve2zfZ97XLe4fe9c3E+uOOqPrk4Gf6am5dEG5TuKlu+cn3923z8ZPFgfOSxqhx1euPOKfm4++cAJyRzeb3n9u/uv6V4oAPXi3VPHTRwNX7vFWqLnWVEo++/unv8GulkBZxml13lz/5fv6h2es/HnZj73Y3dXZ272aMuVkEr7hu6grfLxytlMwwBnMymdR/el7+e1pjvAgudd2U5/uFO5SSnbfdVib29g5MlEq9/73avPvuuI93uWaEKjX8epsPr9jTy30tk2m47oSaNT0jndKBpwa5I9Lphsduq31/8v3Fuu7v775i5wOAgu/nHlAKJp1uONbzcimt1a3GqNcymfrLOjryE0RwEaDuSqfr7/T9/EVKYcK3Py7e/06YXOz7+duUQj1Qc2xPT9EMHBg+IIKC66bO6egojBKRn4vgMddN/cb3C6cpJWcag+symdTDnlf4qdayL4Cp6XQq99y2783+Rc/2XQ+H23zY2fn2LsYkbheRrOs2/Jvvd41XSn/PGPnPTKZhThAUZgByNCD/nk43vOr7uVuUUiPXrSseP2QI+np7kw+JqPdct/4sz8u5WqtfGaOezGTqfxkEuZMBNQVQ/zedrn/Q9/M/UgrfKJWcaaNG7fK27+d/J4LaTCZ1ZDbbPcJxzFwRBK6buiQICocC8n0Av02nU//h+/lWpfCvxuCHmUzqJc/L36g19lBKTV61ate1gwcXHhXBB66bOr2j492vi4TXAPJMOt3wM98vTFZKzjVGZmUyDfd5Xm6m1uoAwLSm0yPf9P3cfwF6e9etP7yzs2u4MfpuQP1POl1/oeflxmmtLhdR8123fnYQ5KYB6tuAvjqd3vV5QA0OgvyTxWLVyaNG7bSqo6PwhDFYkcmkTslmC3s5jlwvgmddN/XjIMgfD+B8ADen06nf+37+CqUwNgzVjKam+tc9L3+P1hjW2Fg/ftGi94ckk/3zRPCG66a+EwTdBwHmSkD+kE433OT7halKySRj5CeZTMMCz8tdp7XaW2tz2te/PnK57xceB8xHrttw4qfTw7oNUC+k0/VXeV7uWK1Vi4i61XXr7w2C/GUADlHKuaixcZdO38/fpRR2/POf6ycMGdJdJyL3AngrnU61eF7+m1rjKhHc77qptiDInw3gJED9PJ2uf9r3879RCo1hqM9satp1meflH1EKva6bOsHz3tlDa+dGEbzsuqkrg6AwEZDvAHJ7Ot0wz/MKl2oth4WhuaSpaWQQBIU5gHyltrZ49KpVqB4wIPl7EfWO69afHwS5bwDqR4B6KJ2uv97zcs1aq1NFzK9cd+QTQZD/JYAMkDgrnd75Pd/PPwig5Lqp4zo6CiNF5BZj8MdMJnW55+WP0hrfFZG56XT93b6fv1gpHAHg39LpVNb38/+hFHbt6XGOGTgwqYH19xmDfCaTOtfzuvfV2vxUBI+6buqaICicAcjpSuGaxsbUo75f+JlSMjoMzTlNTSMLQVC4T8Ro1204JpvtqnccfZuILHLdhu8HQe4IQF1sjNydyTTM9bzcd7VWRxmjL8tkdn0tCPKzATQkErXH1db2lNasUQ+K4F3XTU0JgnwTgF8YgycymdSvPC9/qtZoFlHXu279Q0GQ/zGAMUqp8xob67t8P/97ANWumzo6CN7bGSjdYQz8TCb1Pc8rHKa1XCqCea6buj0ICt8BZCIgV6TTDa/4fuEmpWT3YjE8AVjdm0wOfcQYWZbJNJyZzXalHUf/WkSedt2GnwdB7iRAnW2MtGUyDff7fu4qpdQ3jQkvyGR2WxoE+XtFUOe6qQmLF+d2LJXUXSLodN3URUGQ/xcAPzBG3ZvJ1N/qefkWrXGsMWpmJlP/ou/n25TC1xIJ58S99tr5o46OwuMiWO66qdOCoHsfwHzT8/KXZzKpn3heYZLWMlUEN7lu6g++X7iyr8+ZNWbMzv/vzS63dq1nth7Q2tzSDgDTTpk2qKW5JddyRsuQjTlGR0d+Qnu7VGyfDvqrhQvzO2Wz3aNt56BofPoBSXGQzXbVd3Z2NdrOQdEoZ21auSyyvHf5K4Asn97c8kenSr+hBdfOunPWRo2WwhDTBw9eXr2lMlL5aC1f1docYzsHRcMYsC9CTDhOIm2MHm87B0WjnLVp5Zv//PnzQwAntU5pHVGL3rW/vP32Tbj3Vt1RLL7XH3k4KrtkMuwypsrYzkFRkVtsJ6BolErqdccxH9jOQVFhbRIRERGVVxDkL1u6dCkvi8RANlvYy/dzJ9rOQdHw/dxVX/4uqgTZbPdorqGJj3LWZsXuLWIMDuztreOCzlgwwwC1t+0UFA0RPc52BoqG42AEIGyKFhPlrM2K/XBWCpe/+eZw9ryPgURCdyil3rOdg6Ihoi62nYGiEYZ9LwNO8OXvpErA2iQiIiIiIiIiIiIiIiIiIiLaulXkrqKXTplS1xMOmKCU6k/2JB+/dv6162xnok134aQLBxQHFo9um9M233YW2jwzzpkx3JSKE4xSq4vJ4iOzZ88u2s5Em2balGkjE0YdCu3khqwd8tzM+TPZtDAGWs9sPQAA2ua2vbglz1Nxt6JOnTo12WsGPA9gEoDDirV9D9jORJvugrMu2KU4sO9nUHKh7Sy0eVrOaBkSlkqvCnAEgIOq+pNvXzplSp3tXLTxLjj7gt0d47wgohvFhCetHLjyCduZaPN997SpOynIfQDGbelzVdzgoqpYdSIUCrPmzprUNqetBUoNbT29tcF2Lto02qgbAT3Gdg7afNqRQwR4um3OjSfNumPWJQBeWxfWHGc7F208HeqjFKStbW7b9LY5N54DyJ4zzpkx3HYu2iyqlKi6VZQ8rZTIlj5Z5fW5EOwNqL/edy3qDZUwowDk7IWiTdU2d9bEljNa9lcOfmM7C22eYsk8V11d/TwATJ8+vdqsNfuEIlfbzkUbr21O2/UAMH3K9LQYOQ5K3rn+tuuX285Fm256c8slENNutKpSsuWXRFTczAVEhitI/rOHSsn7AHa0F4iIAODmu27+8Prbrl8+vXnaN+ST8GUFde9Nd97UYTsXbToJZQwgB0GkNO2UaYNs56FN0zqldbQoOeSGuTdeU65zVtzgQhQ6BPKVzx4bYLAKsdhmJiL6VGvzBRcI9HVQemrbnLYrbeehTdM6pXX01KlTh7bNbZvdNmfWoRD0JZKKe4xUKoOZSlRVS3PLXUowWYDJLWe1nLklT1lxl0WUqNeg8GMAP/rLSHpcCFW2PeqJ6B9raW45TIBTP+z58ID58+eHtvPQplNGxif7kgkAVwOAUmo7aM024BUqgVJrCYk6BUBBmgVIKqWe3LLnrDDLe5e/MnzgsOXTm1v+KMCuSvDTG+6ctcp2LqJ/egqHK8GewwfukGttbvn0KSU/ueGOG2+xnIw2VkLfqUrmttbm1oVQRovIc22338DBRYW6bs7N+c/+Pb255QMFJG+4vW2ZxUhbr9YprSN4mxsR0ZYz45wZw2dOmlllOwcRERERERERERERERERERERERERERERERERERERERERERERERERERERERERERERERERERERERERERERERERERERERERERERERERERERERERbVHKdgAiqnyTj5p8JsQc6TjJlnkPzlsJACceOelnEPT/16Pzf2g7HxGVl7YdgIgqX7Gq+BCgxpqwdD0AnHjUCScJ8D1AP2o7GxGVH2cuiCgSJx55wtEC9SBEnQ0lv1CCGzlrQfTPybEdgIjiYfHS19/aZ4+9d4HCTACLP1y34sx8Pm9s5yKi8uNlESKKjChZDgBQWLFgwYKS5ThEZAlnLogoEicdefwYgb4dSn4CUWftvfs+7y9+e3HWdi4iKj/OXBDRZps4cWKtgb4TwH33Pvy7KwH1S2i5ZtLESQ22sxFR+XFwQUSbrdbU/ALADqUwnA4Aa2XtVRC8ixBzZ86cyb8zREREREREREREREREREREREREREREREREREREREREREREREREREREREREREREREREREREREREREREREREREREREREREREREREREREREREREREREREREREREREREREREREREREREREREREREREREREREREREREREREREREREREREREREREREREREREcfG/ZETv/yFWrG0AAAAASUVORK5CYII="
],
"image/svg+xml": [
"\n",
"\n"
],
"text/html": [
"\n",
"\n"
],
"text/plain": [
"Plot(...)"
]
},
"execution_count": 29,
"metadata": {
"comm_id": "c84ce65f-154b-4ea3-8847-174338405f68",
"reactive": true
},
"output_type": "execute_result"
}
],
"source": [
"using Gadfly, Interact\n",
"@manipulate for λ=0:.1:5, mu=0:.1:5\n",
" problem = minimize(norm(A * x - b)^2 + λ*norm(x)^2 + μ*norm(x, 1), \n",
" x >= 0)\n",
" solve!(problem)\n",
" plot(x=x.value, Geom.histogram(minbincount = 20), \n",
" Scale.x_continuous(minvalue=0, maxvalue=3.5))#, Scale.y_continuous(minvalue=0, maxvalue=6))\n",
"end"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Quick convex prototyping"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"## Variables"
]
},
{
"cell_type": "code",
"execution_count": 30,
"metadata": {
"collapsed": false,
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [
{
"data": {
"text/plain": [
"Variable of\n",
"size: (1, 1)\n",
"sign: NoSign()\n",
"vexity: AffineVexity()"
]
},
"execution_count": 30,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Scalar variable\n",
"x = Variable()"
]
},
{
"cell_type": "code",
"execution_count": 31,
"metadata": {
"collapsed": false,
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [
{
"data": {
"text/plain": [
"Variable of\n",
"size: (4, 1)\n",
"sign: NoSign()\n",
"vexity: AffineVexity()"
]
},
"execution_count": 31,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# (Column) vector variable\n",
"y = Variable(4)"
]
},
{
"cell_type": "code",
"execution_count": 32,
"metadata": {
"collapsed": false,
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [
{
"data": {
"text/plain": [
"Variable of\n",
"size: (4, 4)\n",
"sign: NoSign()\n",
"vexity: AffineVexity()"
]
},
"execution_count": 32,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Matrix variable\n",
"Z = Variable(4, 4)"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"# Expressions"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"Convex.jl allows you to use a [wide variety of functions](http://convexjl.readthedocs.org/en/latest/operations.html) on variables and on expressions to form new expressions."
]
},
{
"cell_type": "code",
"execution_count": 34,
"metadata": {
"collapsed": false,
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [
{
"data": {
"text/plain": [
"AbstractExpr with\n",
"head: +\n",
"size: (1, 1)\n",
"sign: NoSign()\n",
"vexity: AffineVexity()\n"
]
},
"execution_count": 34,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"x + 2x"
]
},
{
"cell_type": "code",
"execution_count": 35,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"AbstractExpr with\n",
"head: +\n",
"size: (1, 1)\n",
"sign: NoSign()\n",
"vexity: ConcaveVexity()\n"
]
},
"execution_count": 35,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"e = y[1] + logdet(Z) + sqrt(x) + minimum(y)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Examine the expression tree"
]
},
{
"cell_type": "code",
"execution_count": 38,
"metadata": {
"collapsed": false,
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [
{
"data": {
"text/plain": [
"AbstractExpr with\n",
"head: logdet\n",
"size: (1, 1)\n",
"sign: NoSign()\n",
"vexity: ConcaveVexity()\n"
]
},
"execution_count": 38,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"e.children[2]"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"# Constraints\n",
"\n",
"A constraint is convex if convex combinations of feasible points are also feasible. Equivalently, feasible sets are convex sets.\n",
"\n",
"In other words, convex constraints are of the form\n",
"\n",
"* `convexExpr <= 0`\n",
"* `concaveExpr >= 0`\n",
"* `affineExpr == 0`"
]
},
{
"cell_type": "code",
"execution_count": 39,
"metadata": {
"collapsed": false,
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [
{
"data": {
"text/plain": [
"Constraint:\n",
"<= constraint\n",
"lhs: Variable of\n",
"size: (1, 1)\n",
"sign: NoSign()\n",
"vexity: AffineVexity()\n",
"rhs: 0\n",
"vexity: AffineVexity()"
]
},
"execution_count": 39,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"x <= 0"
]
},
{
"cell_type": "code",
"execution_count": 40,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"Constraint:\n",
"<= constraint\n",
"lhs: AbstractExpr with\n",
"head: qol_elem\n",
"size: (1, 1)\n",
"sign: Positive()\n",
"vexity: ConvexVexity()\n",
"\n",
"rhs: AbstractExpr with\n",
"head: sum\n",
"size: (1, 1)\n",
"sign: NoSign()\n",
"vexity: AffineVexity()\n",
"\n",
"vexity: ConvexVexity()"
]
},
"execution_count": 40,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"x^2 <= sum(y)"
]
},
{
"cell_type": "code",
"execution_count": 42,
"metadata": {
"collapsed": false,
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [
{
"data": {
"text/plain": [
"Constraint:\n",
"sdp constraint\n",
"expression: AbstractExpr with\n",
"head: +\n",
"size: (4, 4)\n",
"sign: NoSign()\n",
"vexity: AffineVexity()\n",
"\n"
]
},
"execution_count": 42,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"M = Z \n",
"for i = 1:length(y)\n",
" M += rand(size(Z))*y[i]\n",
"end\n",
"M ⪰ 0"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"# Problems"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {
"collapsed": false,
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [
{
"data": {
"text/plain": [
"Problem:\n",
"minimize AbstractExpr with\n",
"head: +\n",
"size: (1, 1)\n",
"sign: NoSign()\n",
"vexity: ConvexVexity()\n",
"\n",
"subject to\n",
"Constraint:\n",
">= constraint\n",
"lhs: Variable of\n",
"size: (1, 1)\n",
"sign: NoSign()\n",
"vexity: AffineVexity()\n",
"rhs: AbstractExpr with\n",
"head: maximum\n",
"size: (1, 1)\n",
"sign: NoSign()\n",
"vexity: ConvexVexity()\n",
"\n",
"vexity: ConvexVexity()\n",
"current status: not yet solved"
]
},
"execution_count": 12,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"x = Variable()\n",
"y = Variable(4)\n",
"objective = 2*x + 1 - sqrt(sum(y))\n",
"constraint = x >= maximum(y)\n",
"p = minimize(objective, constraint)"
]
},
{
"cell_type": "code",
"execution_count": 43,
"metadata": {
"collapsed": false,
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [
{
"data": {
"text/plain": [
":Optimal"
]
},
"execution_count": 43,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# solve the problem\n",
"solve!(p)\n",
"p.status"
]
},
{
"cell_type": "code",
"execution_count": 44,
"metadata": {
"collapsed": false,
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [],
"source": [
"x.value"
]
},
{
"cell_type": "code",
"execution_count": 45,
"metadata": {
"collapsed": false,
"slideshow": {
"slide_type": "subslide"
}
},
"outputs": [
{
"data": {
"text/plain": [
"1x1 Array{Float64,2}:\n",
" 0.500004"
]
},
"execution_count": 45,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# can evaluate expressions directly\n",
"evaluate(objective)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Pass to solver\n",
"\n",
"call a `MathProgBase` solver suited for your problem class\n",
"\n",
"* see the [list of Convex.jl operations](http://convexjl.readthedocs.org/en/latest/operations.html) to find which cones you're using\n",
"* see the [list of solvers](http://www.juliaopt.org/) for an up-to-date list of solvers and which cones they support"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"to solve problem using a different solver, just import the solver package and pass the solver to the `solve!` method: eg\n",
"\n",
" using Mosek\n",
" solve!(p, MosekSolver())"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Warmstart"
]
},
{
"cell_type": "code",
"execution_count": 23,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"elapsed time: 0.021169205 seconds (3092632 bytes allocated)\n",
"elapsed time: 0.0075224 seconds (3081120 bytes allocated)\n"
]
}
],
"source": [
"# Generate random problem data\n",
"m = 50; n = 100\n",
"A = randn(m, n)\n",
"x♮ = sprand(n, 1, .5) # true (sparse nonnegative) parameter vector\n",
"noise = .1*randn(m) # gaussian noise\n",
"b = A*x♮ + noise # noisy linear observations\n",
"\n",
"# Create a (column vector) variable of size n.\n",
"x = Variable(n)\n",
"\n",
"# nonnegative elastic net with regularization\n",
"λ = 1\n",
"μ = 1\n",
"problem = minimize(norm(A * x - b)^2 + λ*norm(x)^2 + μ*norm(x, 1), \n",
" x >= 0)\n",
"@time solve!(problem)\n",
"λ = 1.5\n",
"@time solve!(problem, warmstart = true)"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# DCP examples"
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {
"collapsed": false,
"slideshow": {
"slide_type": "subslide"
}
},
"outputs": [
{
"data": {
"text/plain": [
"AbstractExpr with\n",
"head: +\n",
"size: (1, 1)\n",
"sign: NoSign()\n",
"vexity: AffineVexity()\n"
]
},
"execution_count": 17,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# affine\n",
"x = Variable(4)\n",
"y = Variable (2)\n",
"sum(x) + y[2]"
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
"AbstractExpr with\n",
"head: +\n",
"size: (1, 1)\n",
"sign: NoSign()\n",
"vexity: ConvexVexity()\n"
]
},
"execution_count": 18,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"2*maximum(x) + 4*sum(y) - sqrt(y[1] + x[1]) - 7 * minimum(x[2:4])"
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {
"collapsed": false,
"slideshow": {
"slide_type": "subslide"
}
},
"outputs": [
{
"data": {
"text/plain": [
"AbstractExpr with\n",
"head: +\n",
"size: (4, 1)\n",
"sign: NoSign()\n",
"vexity: NotDcp()\n"
]
},
"execution_count": 19,
"metadata": {},
"output_type": "execute_result"
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"WARNING: Expression not DCP compliant. Trying to solve non-DCP compliant problems can lead to unexpected behavior.\n"
]
}
],
"source": [
"# not dcp compliant\n",
"log(x) + x^2"
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {
"collapsed": false,
"slideshow": {
"slide_type": "subslide"
}
},
"outputs": [
{
"data": {
"text/plain": [
"AbstractExpr with\n",
"head: qol_elem\n",
"size: (4, 1)\n",
"sign: Positive()\n",
"vexity: ConvexVexity()\n"
]
},
"execution_count": 20,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# $f$ is convex increasing and $g$ is convex\n",
"square(pos(x))"
]
},
{
"cell_type": "code",
"execution_count": 21,
"metadata": {
"collapsed": false,
"slideshow": {
"slide_type": "subslide"
}
},
"outputs": [
{
"data": {
"text/plain": [
"AbstractExpr with\n",
"head: qol_elem\n",
"size: (4, 1)\n",
"sign: Positive()\n",
"vexity: ConvexVexity()\n"
]
},
"execution_count": 21,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# $f$ is convex decreasing and $g$ is concave \n",
"invpos(sqrt(x))"
]
},
{
"cell_type": "code",
"execution_count": 22,
"metadata": {
"collapsed": false,
"slideshow": {
"slide_type": "subslide"
}
},
"outputs": [
{
"data": {
"text/plain": [
"AbstractExpr with\n",
"head: geomean\n",
"size: (4, 1)\n",
"sign: Positive()\n",
"vexity: ConcaveVexity()\n"
]
},
"execution_count": 22,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# $f$ is concave increasing and $g$ is concave \n",
"sqrt(sqrt(x))"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Julia 0.3.8",
"language": "julia",
"name": "julia-0.3"
},
"language_info": {
"name": "julia",
"version": "0.3.8"
}
},
"nbformat": 4,
"nbformat_minor": 0
}