1_lamps.ipynb

{
 "metadata": {
  "name": "",
  "signature": "sha256:0aa34b53130266f7974004049239e72e9be0771ace00c9bd878a312b11d57d56"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import numpy as np\n",
      "import cvxopt as cvx\n",
      "import matplotlib\n",
      "import matplotlib.cm\n",
      "import matplotlib.pyplot as plt\n",
      "from numpy import log,exp\n",
      "%matplotlib inline \n",
      "matplotlib.rcParams['figure.figsize'] = (10, 6)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 1
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "Let us generate some data:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "np.random.seed(420) #makes data reproducible\n",
      "\n",
      "Ides = 0.1\n",
      "pmax = 2.\n",
      "N = 30\n",
      "M = 10\n",
      "slope = 0.6\n",
      "\n",
      "xxx = [-0.5]\n",
      "yyy = [0.]\n",
      "x = []\n",
      "y = []\n",
      "for i in range(N):\n",
      "    s = np.random.randn()\n",
      "    yyy.append(yyy[-1]+slope * s)\n",
      "    xxx.append(xxx[-1]+1.)\n",
      "    x.append(0.5*(xxx[-1]+xxx[-2]))\n",
      "    y.append(0.5*(yyy[-1]+yyy[-2]))\n",
      "    \n",
      "ymax = max(yyy)    \n",
      "\n",
      "ylamps = [ymax +2.] * M\n",
      "xlamps = sorted(np.random.rand(M)*N)\n",
      "\n",
      "plt.plot(xxx,yyy)\n",
      "plt.scatter(x,y,color = 'red')\n",
      "plt.scatter(xlamps,ylamps, color = 'green')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 2,
       "text": [
        "<matplotlib.collections.PathCollection at 0x7f7cf725b810>"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAlUAAAFwCAYAAACRj46qAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XmczWX/x/HXZRmGLMVEoZn2u52yRNuMyFaJNu6S3ELr\n3b7XL3fLfbcX7UpllFEpUlkipgVFllAqLUNkOSTCbLh+f1xHBrOf7znfs7yfj4dHM2fO+X4/feeY\nebuu6/u5jLUWEREREQlNFb8LEBEREYkHClUiIiIiHlCoEhEREfGAQpWIiIiIBxSqRERERDygUCUi\nIiLigWqhHsAYUwP4DEgKHm+MtfY/oR5XREREJJYYL/pUGWNqWWu3GmOqAjOAf1trZ4d8YBEREZEY\n4cn0n7V2a/DDGrjRKnUUFRERkYTiSagyxlQxxswHVgNTrLVzvDiuiIiISKzwaqRqh7W2BdAUaGOM\nOdqL44qIiIjEipAXqhdlrd1kjJkOdAa+K/o1Y4ymBEVERCRmWGtNRZ4f8kiVMaahMaZe8ONkoCPw\nfQnF6c8ef+677z7fa4i2P7omui66Lrouuia6Ln7/qQwvRqoOAEYYY6rgQtpb1toJHhxXREREJGaE\nHKqstYuAEz2oRURERCRmqaO6z9LT0/0uIeromhRP16V4ui7F03XZm65J8XRdvONJ889yncgYG6lz\niYiIiITCGION9EJ1EREREVGoEhEREfGEQpWIiIiIBxSqRERERDygUCUiIiLiAYUqEREREQ8oVImI\niIh4QKFKRERExAMKVSIiIiIeUKgSERER8YBClYiIiIgHFKpEREREPKBQJSIiIuIBhSoRERERDyhU\niYiIiHhAoUpERETEAwpVIiIiIh5QqBIRERHxgEKViIiIiAcUqkREREQ8oFAlIiIi4gGFKhEREREP\nKFSJiIiIeEChSkRERMQDClUiIiIiHlCoEhEREfGAQpWIiIiIBxSqRERERDygUCUiIiLiAYUqERER\nEQ8oVImIiIh4QKFKKiWwJcCclXMIbAlE9LXi6BqGTzivbXmOre9t5a9BqNfOr/NK/FCokgrLWpRF\n6tOpdBzZkdSnU8lanBWR14qjaxg+4by25Tm2vreVvwahXju/zivxxVhrI3MiY2ykziXhE9gSIPXp\nVHK35f79WHK1ZJbdsIyU2ilhe604uobhE85rW55j63tb+WsQ6rXz67wS3YwxWGtNRV6jkSqpkJw/\nc0iqmrTbY9WrVifnz5ywvlYcXcPwCee1Lc+x9b2t/DUI9dr5dV6JPwpVUiFp9dMo2F6w22OF2wtJ\nq58W1teKo2sYPuG8tuU5tr63lb8GoV47v84r8UehSiokpXYKw7sPJ7laMnVr1CW5WjLDuw8v11B3\nKK8VR9cwfMJ5bctzbH1vK38NQr12fp1X4k/Ia6qMMU2BTKARsAN42Vo7tJjnaU1VHAlsCZDzZw5p\n9dMq/AMklNeKo2sYPuG8tuU5tr63lb8GoV47v84r0akya6q8CFWNgcbW2gXGmH2AuUB3a+33ezxP\noUpERERigi8L1a21q621C4IfbwaWAE1CPa6IiIhILPF0TZUxJg1oDnzl5XFFREREop1noSo49TcG\nuD44YiUiIiKSMKp5cRBjTDVcoBpprX2/pOcNHjz474/T09NJT0/34vQiIiIiIcnOziY7OzukY3jS\nUd0Ykwmss9beVMpztFBdREREYoJfd/+dAnwGLAJs8M9d1tpJezxPoUpERERigi+hqtwnUqgSERGR\nGKG9/0RERER8olAlIiIi4gGFKhEREREPKFSJiIiIeEChSkRERMQDClUiIiIiHlCoEhEREfGAQpWI\niIiIBxSqRERERDygUCUiIiLiAYUqEREREQ8oVImIiIh4QKFKRERExAMKVSIiIiIeUKgSERER8YBC\nlYiIiIgHFKpEREREPKBQJSIiIuIBhSoRERERDyhUiYiIiHhAoUpERETEAwpVIiIiIh5QqBIRERHx\ngEKViIiIiAcUqkREREQ8oFAlIiIi4gGFKhEREREPKFSJiIiIeEChSkRERMQDClUiIiIiHlCoEhER\nEfGAQpWIiIiIBxSqRERERDygUCUiIiLiAYUqEREREQ8oVImIiIh4QKFKRERExAMKVSIiIiIeUKgS\nERER8YAnocoYM9wYs8YYs9CL44mIiIjEGq9Gql4DOnl0LBEREZGY40mostZ+AWzw4lgiIiIisUhr\nqkREREQ8oFAlIiIi4oFqkTzZ4MGD//44PT2d9PT0SJ5eREREpFjZ2dlkZ2eHdAxjrfWkGGNMGvCB\ntfa4Er5uvTqXiIiISDgZY7DWmoq8xquWCqOAmcARxpjlxph+XhxXREREJFZ4NlJV5ok0UiUiIiIx\nwreRKhEREZFEp1AlIiIi4gGFKhEREREPKFSJiIiIeEChSkRERMQDClUiIiIiHlCoEhEREfGAQpWI\niIiIBxSqRERERDygUCUiIiLiAYUqEREREQ8oVImIiIh4QKFKRERExAMKVSIiIiIeUKgSERER8YBC\nlYiIiIgHFKpEREREPKBQJSIiIuIBhSoRERERDyhUiYiIiHhAoUpERETEAwpVIiIiIh5QqBIRERHx\ngEKViIiIiAcUqkREREQ8oFAlIiIi4gGFKhEREREPKFSJiIiIeEChSkRERMQDClUiIiIiHlCoEhER\nEfGAQpWIiIiIBxSqRERERDygUCUiIiLiAYUqEREREQ8oVIlEiW+/hVtv9bsKERGpLIUqkSjx5pvw\n+OMwdqzflYiISGUYa21kTmSMjdS5RGJR8+bQqxc8+6wbtapXz++KREQSlzEGa62p0GsUqkT8t3Il\nHH88rF0LV10F1avDc8/5XZWISOKqTKjyZPrPGNPZGPO9MeZHY8ztXhxTJJFMnAidOkHVqvDII24K\ncOZMv6sSEZGKCDlUGWOqAM8CnYBjgN7GmH+EelyRRDJhAnTt6j7ed194+mkYOBAKCvytS0REys+L\nkarWwFJr7TJrbSEwGujuwXFFEkJBAUyb5kaqdrrwQkhLg8ce860sERGpIC9CVRPgtyKfrwg+JsBf\nf/ldgUS7L76Af/wDUlJ2PWaMW1P11FPw44/+1SYiIuVXLZInGzx48N8fp6enk56eHsnT+6JnT+jc\nGW6+2e9KJFoVnforKjUV7rkHBg1yI1mmQsslRUSkIrKzs8nOzg7pGCHf/WeMORkYbK3tHPz8DsBa\nax/Z43kJefffihWQkQFXXqlgJcU7+mjIzISWLff+2vbtcPLJcHWfTfRr+4ObEyw6pCUiImHh191/\nc4DDjDGpxpgkoBcw3oPjxoWmTWH6dHjxRdfYUaSoX3+F9evhxBOL/3rVqjCs5yTuuD6XtWf2dsNX\nWVmRLVJERMrFkz5VxpjOwBBcSBturX24mOck5EjVTitWQHq6G7G65Ra/q5Fo8dxzMGcOvP56CU8I\nBCA1ldtyB7OSJrzJpZCcDMuWacRKRCSMKjNS5cmaKmvtJOBIL44Vr5o2hexsNxUIwWAVCEBOjqZ0\nEtiECdC3bylPyMmBpCTuy/0Px7GISXSic/VZ7nG9Z0REooo6qkfYzjVWA1vN59Zxp0BSkrunfvhw\n6N3b7/IkgnJzYf/9Yfly15uqWMGRKnJzmcxZXMmLLK7ZitrLlyhUiYiEkW8d1aX8mjaF6WPWM2x0\nHR7LvQY2bnS/Xfv3d79AJWFkZ0OLFqUEKnDBafhwSE6mU90vOaXqlww+83MFKhGRKKRQ5YOmBb+Q\nvc85vMwAHiO4wKp6dTelIwmjpFYKe+nd262hmjqVJxd2JHPOUcybF/byRESkgiLap0qC0tJosm0Z\n08kgnWwshtsKn3VrqyQhWOtC1dix5XxBSgqkpLA/bm/AgQPhyy+hmv4Gi4hEDY1U+SE4pdMkeYMb\nsTIDefS8mZrSSSA//gj5+XDccRV/bd++ULcuPPOM93WJiEjlaaG6n4J3/62scQjp5zdgwAC47Ta/\ni5JIeOopWLIEhg2r3OuXLoW2beHrrzXAKSISDlqoHmtSUqBVK5oc34DsbHjlFTe1I/Gv3OupSnD4\n4XDTTXD11W4qUURE/KeRqiiycqVrt9C/P9x+u9/VSLhs3gwHHAC//w516lT+OAUFcNJJbn/Aiy/2\nrj4REdFIVcxr0sRtaTN8uEas4tknn0CbNqEFKnAtzoYNgxtvhA0bvKlNREQqT6EqyuwMVs88A3Pn\n+l2NeC4QYMKItXQ9fbMnh2vbFnr2hFtv9eRwIiISAoWqKNSkiftFOWWK35WIp7KysAelMmFcAV3/\ne6pnGyP/978webJrJioiIv5RqIpS7du7ESuJE4EA9O/P4rxDqW4LODL/G8+66NetC88+63pX5eZ6\nUKuIiFSKQlWUOuMMmDnTLUaWOBDcGHkCXenKBAx42kW/e3c4/nh48EFPDiciIpWgUBWl9t3X3TY/\ne7bflYgn0tKgoIAJdKUbH7nHCgs9bTL1zDNu4frChZ4dUkREKkChKoppCjCOpKTw59BM5tOC9Drz\nIDnZ3ebpYRf9Aw5w66sGDIDt2z07rIiIlJNCVRTLyFCoiidT6l3AaR1qkPzJh26D5N69PT9H//5Q\nsyY895znhxYRkTKo+WcU27QJDjwQ1q1zvygltvXrBy1bwjXXhPc8P/wAp57qWnIcdFB4zyUiEq/U\n/DPO1K0Lxx4Ls2b5XYmEascOmDgRunQJ/7mOPBKuv15b2IiIRJpCVZTTFGB8mD/f3XxwyCGROd9t\nt7kZxrdf2QRz5njSukFEREqnUBXlMjJg2jS/q5BQhbqBckUlJcHLF37MjQO38MeZF0JqqmfNRkVE\npHhaUxXltmyBRo1gzRqoXdvvaqSy2rZ1PaTOPDNCJwwEIDWV63IfYSu1GM4V7o7DZcs8veNQRCRe\naU1VHKpdG1q0gBkz/K5EKmvdOvjuO7d4PGKCzUb/y11MoSPTSfe02aiIiOxNoSoGaAowtk2e7HqO\n1agRwZMGm43WYTPPcQ0DGUZuQVVPm42KiMjuFKpigJqAxrZIr6cC3BTf8OGQnMw5dT+jRdWFPHDW\n55r6ExEJI62pigF5edCwIfz+u2uzILFj+3a3Jm7BAmja1IcCAgHIyWF18sEc374hU6bACSf4UIeI\nSIzRmqo4VbMmtG4Nn3/udyVSUbNnQ5MmPgUqcCNTrVrR+NiG2sJGRCTMFKpiRPv2WlcVi3yZ+itB\n//5Qq5bbeFlERLynUBUj1AQ0NkVTqDIGhg1zrR2WLfO7GhGR+KM1VTGioMCtq8rJgf3287saKY9V\nq+CYY2DtWqhWze9qdnnoIdei46OPXNASEZG9aU1VHEtKgnbt4NNP/a5EymvSJOjYMboCFcCtt8Jv\nv8Ho0X5XIiISXxSqYoimAGNLNE39FZWUBK+8AjfdBOvX+12NRNLGjXDuubB0qd+ViMQnhaoYolAV\nOwoLYepU6NzZ70qK16YNXHQR9OsH+fl+VyORkJ8PPXrA3Lnw8st+VyMSnxSqYsiJJ7ppm7Vr/a5E\nyjJzJhx2mOtRFa0efdRNTXbvDlu3+l2NhNOOHS5A16/vOvyPHAnbtvldlUj8UaiKIdWqwWmnQXa2\n35VIWaJ16q+oGjXg7bdd8OvUyU0NSXy6/XZYvhzefBOOPRZSU124EhFvKVTFmIwMmD4hF+bMcd2y\nJSrFQqgCF9Rfew2aN3e90PSWij9DhsCHH8L770NysnusXz/3fRcRbylUxZiM/ElMy/zN3VaWmgpZ\nWX6XJHtYvhxWr4aWLf2upHyqVIGhQ6FLFzj9dFi5EpeuFNxjV/D7N+bVjTz6KEycCA0a7Ppyr15u\nzd+6df6VKBKPFKpiSSDACfefzzrbgN831oLcXNcmW7/4osro0XDeeVC1qt+VlJ8xrinov/4Fp7bY\nzE/NMhTcY1VWFqSm8lnGfVzdP5+Prp1IWtruT6lXD7p1g1GjfKlQJG4pVMWSnByq1KjOGXzKdDLc\nY9Wru46gEhWshREjoG9fvyupnFsvD3DHn3dyRv5kFm9squAeawIB6N+fb3MP5sItrzGKf9L8gfOL\n/f5dfrmmAEW8plAVS9LSoKCADKYzjfbuscJC9vpnqPhm7lzIy4NTTvG7kkrKyWFQrZE8zi10YCqz\naaXgHktyclhZLZWuTOAJbqYDn5T4/Wvf3vUpW7Ag8mWKxKuQQpUx5gJjzGJjzHZjzIleFSUlSEmB\n4cNpX2Mm082ZbtXp8OHucYkKI0bAZZfF8PYvweDem9G8zAC68RHZeScruMeIjfsdTJfNb3M1z3Mp\nb7oHS/iHV9WqbkT19dcjWqJIXAtp7z9jzJHADuAl4BZr7bxSnqu9/zxi1wZofFR9vvp4E2knNSj7\nBRIRBQXQpAnMng0HH+x3NSHIynJTftWrMz2vLRfXHMerb9bk7LP9LkxKk5/vbjY4pvqPDP2sOSap\nugtUw4dD797Fvubnn6FtW1ixwnXaF5FdIr73n7X2B2vtUiBW/10ek8z+KaR3qM70hQpU0WTCBDj6\n6BgPVOB+AS9bBlOnkrFiJB98XJMrrtBegdGsaHPPpyccgVnuvn8sW1ZioAI49FA46ijXckFEQqc1\nVTGqfXttWRNtYnmB+l5SUqBVK0hJoU0b9/v55pth2DC/C5Pi3HHHruaeVauy2/evLOpZJeKdamU9\nwRgzBSi62YYBLHC3tfaDipxs8ODBf3+cnp5Oenp6RV4uRWRkuFvgrY3h9TtxZN06F3JHjPC7kvA4\n9lj49FPXZWHTJrjlFr8rkp2GDoUPPoAZM3Y196yICy6AG25wvdUaN/a+PpFYkZ2dTXaIW5aEtKbq\n74MYMx24WWuqIsdaaNrUbVlz+OF+VyPPPANffulGCuLZihUuWF1/PVx5pd/VyJgxLhDNmOFailVW\nv35wzDEKyyJFRXxN1Z7n9/BYUgZjNAUYTTIz42jqrxRNm8LTT2t9VTT4/HO4+mq3HiqUQAW7pgD1\n716R0ITaUuE8Y8xvwMnAh8aYid6UJeWRkaFQFQ2++w5WrYIzz/S7ksho08b14yos9LuSxGWt636/\nc9/GUJ12mrt7cM7HG7Q9kUgIQr37b5y1tpm1Ntlae4C1totXhUnZdoYq/evSXyNGwKWXxta2NKGo\nX9+1PVq40O9KEte337pQ69Wm3cbA5Scu5LVuY7Q9kUgIdPdfDDv4YLcwdckSvytJXNu3wxtvuIaf\niaRtW5g50+8qEtfYsW5/Sc9uUgkE6Dv+fN7e3pPcjfnankikkhSqYpymAP31ySdw4IGuP1UiaddO\nocpP48a5UOWZnBya1QxwEnMZR/DA2p5IpMIUqmKcQpW/4qo3VQW0awezZvldRWJavtz9OfVUDw8a\n3J7ocl7nNfq5x7SvqEiFKVTFuIwM11Zhxw6/Kwm/P/6AQYPgr7/8rsTZtAk++gh69fK7ksg7/HDY\nvBlWrvS7ksQzbhycfTZUK7PLYAUE9xXtUXMSX9OK32ocpn1FRSpBoSrGNW0K++0Hixb5XUn4Pf44\nvPee23Vj2zafiwkEeOfRX8lol0/Dhj7X4gNj3LoqjVZF3rhx0KNHGA7cuzfJy3/goh6FZN60oNTt\nbUSkeApVcSARpgDXrIGXXoKvvnK3ft90k4/FZGVBaiojHl5F308uS9i7pLSuKvLWr3ftLDp2DNMJ\nUlLod/v+vP5Obd1VLFIJClVxoH17mDbN7yrC6+GH4ZJL4JBD4J133ALxZ57xoZBAAPr355fcxny/\n/TC6FoxN2LuktK4q8j780PVDq8x2NOXVurVbo/7FF+E7h0i8UqiKA+nprrvy9u1+VxIeK1a4BeF3\n3eU+r1/f/XL53//cmqaIysmBpCRG0odejCaJwoS9S6pVK9erKi/P70oSx9ixYZr6K8IYbbIsUlkK\nVXGgUSN3W//8+X5XEh4PPggDBuy+2evBB7v1VZdfDgsWRLCYtDRsfgGZXEZfgrsnJ+hdUrVqwVFH\nuekoCb+tW92IdLdu4T/XpZe6v1+bN4f/XCLxRKEqTsTrFOAvv7hNY2+7be+vnXwyPP88nHsurFy4\nPjLba6Sk8MVt46lp8jmxzk9uHiaB75LSFGDkTJ7spub22y/85zrgANey4d13w38ukXiiUBUn4nWx\n+uDBcN110KBB8V+/8EK4qt03nNP8NzZ3OC8i22uMWNmBvvcchPlkKixbltB3SWmxeuR43vCzDJoC\nFKk4YyN0i4cxxkbqXIlo/Xo3JbZ+vVviEw+++86tF/vpJ6hbt4QnBQLYg1K5Iu8ZAqQwlh5UTa7h\nwk4YRo9yc6FJE1i82E25Jrrly93oyapVHm6ZInspLHTT3wsWQLNmkTlnQYF7r3/1lbtBRCTRGGOw\n1lboJ5tGquJEgwZw6KFuBixe3Hcf3HJLKYEKICcHUyOJF7iKzezDLTwe1oXj48a5EKFA5TRr5jaS\n/vVXvyuJb59/7oJNpAIVQFIS/POf8PrrkTunSKxTqIoj8TQFOH8+zJgB11xTxhOD22skUci7nM8k\nOvN8br+wLRxP1G1pSmKM1lVFQiTu+itOv37uPZ8IOzaIeEGhKo7ERagKBGDOHO69PZ8774Tatct4\nfnB7DZKT2bfuDj6qcT4P1P4fE7/2furv999h9mzo3t3zQ8c0rasKL2sjv55qp+bN3cL4eLwJRiQc\nFKriyOmnw5dfuluvY1KwU/msjLtYNHUNA+u9Vb7X9e7t1lBNncohv33Kux8l07ev66HkpTfegJ49\nXSsB2UWhKrzmzt3VvsIPl1+uKUCR8tJC9TjTp48LVW+95fGGq+EWCLg793JzOZOp9CaLK5JHVXrB\n+ejRcPvtLmQecEDo5VkLxx4LL74Ip50W+vHiSUEB7LsvrF4Nder4XU38uecet1D9kUf8Of+6dXDY\nYe6vYr16/tQg4gctVBdeeQX++guuvJLY2rsr2Kl8Ghks5yDXWDOEBee9ermGoeecA1u2hF7evHmu\nc/ipp4Z+rHiTlAQtWsTXTRLRxK/1VDs1bOi2xnmrnAPHIolMoSrO1KjhOiEvWgR33OF3NRUQ7FR+\nDw8ymMFUZ1vIncrvvtuNLl16aehb+IwYAZddprYBJdEUYHj8+OUfbFhbQOuD/d1bUj2rRMpHoSoO\n7bMPTJjg9sd79FG/qymnlBQmXDeRTaYevepM8KRTuTEwbBhs2ADXXgt//lm54xQUuOVeffpUupS4\np1AVBllZjDvtCbpvGkmVg8Pf1LY0nTu7thk//+xbCSIxQaEqTjVoAB9/DC+84LJJtNuxA+6degYP\nvNqUqp987Fmn8qQkN3L3559u0Ovyy90v/7+nRoN3G5a2vc2ECW6RsBoglqxtW7d+TbfeeyQQgP79\nGbetG+cVvO26zvbvH/5tmEpQrRp06QKTJvlyepGYoVAVx5o0cfuF3XuvCxblCRB+ee89qFIFzutb\nD1q18rQb+n77uX/kL13qpgMvvxyOOw6G9p3LhoNOgI4dS93eJjNTvanK0qiRW6z+ww9+VxIncnJY\nVa0ZSziKDIJ9UsLY1LY8OnVSqBIpi+7+SwDz50On9DxG5/WgffIsN581fHjU7Fm3fbsLOU8+6aYZ\nws1a+PT9Pxl2wWQmbO/EuYxnIMM4peY8zPLd7zZcv951ql++vIzO7kKfPm5bof79/a4kDgQCvNTk\nfj4tbMsoLnGPJSeHbful8li/3o3Wrl3r1m6KxDvd/SfFatE0wDt559KrYARzNh7u+1TCnkaNcqNJ\nnTpF5nzGQHqTpYzaZxA/cRjNWcAVvMKxBXMZ8nAuf/yx67lZWdCtmwJVebRtq3VVnklJYexRd9Ej\n6SP35vNgjWGoGjRw0+AzZvhWgkjUU6hKBDk5nJE8m1e4gnP4gO850vephJ0KC2HwYHjooQjfWRfc\n3qYh67mJp1jCUbxQ7d/MWbY/hxziRl0+/3zXXX9SNm1X452NG2HmrwfQ+fshMHWqZ2sMQ6UpQJHS\nKVQlgmCAOJcPeJTb6MRkluc3Ctv+eBXx2mtueu2MMyJ84iLb21C3LiY5mdNf/xdvjKnJzz/DSSfB\nwIGuoWWHDhGuLUYdeyys+G0Hf0ydFzWjoLFqwgS3Q0Kdgxt6vsYwFJ07u3WaIlI8ralKFFlZbsqv\nenWeyh3ESw3v4fNv6vr6szovDw4/HN59F1q39qmIQMCN2KWl7fWLy1rYvFldwsstK4sOlzbippov\n0NV+FFXr9mLNxRe7+yeuuMLvSna3bRvsvz8sXgwHHuh3NSLhVZk1VQpViaRIgLhnSAqTJrkNmP0K\nDUOGuI1a33/fn/OLh4LbDN2beycWw4Pc6/vC6liVlweNG7s7KRs18ruavV10EXTt6u6iFYlnWqgu\npUtJ+Xsq4YEHoGVL6N7d/RCPtC1b4OGH4YEHIn9uCYPgNkPtmMks2rrHomTdXqyZNs3dDRuNgQrc\nFKDWVYkUT6EqQRkDzz3nclbv3m5YP5Keecatozr++MieV8IkuG7vZL5kDq3YRtWQtxlKVGPHwnnn\n+V1Fyc46C6ZMCX3rJ5F4pFCVwKpWhZEjYetWtyg7UsFq+nR44gn4z38icz6JgODC/32T82lWZSWL\narTyvQVALNq+HcaPj+5Q1bSpW0/19dd+VyISfRSqEtzObVx+/x1OPhm++SZ859q6Fa6/3rUryMyE\nI48M37nEB717w7JltD27ITMHf6xF6pUwa5ab9jv0UL8rKZ1aK4gUT6FKqF0bJk6Ea65xdxzdfbdH\n66yKbIszaxa0aOG6Mi9c6PYRkziUkkK77inMWuz/LZPjx7v3WqTcfjs8+2yRfSUrYdw46NHDu5rC\nRa0VRIqnUCWAW2PVr5/7JfTDD9C8uWt+CVRuz8CsLEhNJb9DN+488HV6dMnloYfgjTdc93SJX+3a\n+dtZvbAQbrgBBg2Cnj1dW4xwGzPGrYXKzISzz3ZbuVSUtdG/nmqnU091bRWK7j4gIgpVsofGjd0v\niP/9D3r1gqs7LmXTQceWuenwbgIB6N+f+blH0mrTVL7fdhjf5B/FBWeoIWQiOOII1xF89erIn3vV\nKsjIgJ9+gu++cw00b745vOdcvRquvdatT5wxA044wf2jpKLTY4sXuzVVzZuHp04v1awJp53mmr2L\nyC4KVVIc/1GkAAAajklEQVSsHj3g20/XUTj9C47Nm8OHG08t956B237K4YEdd9OJydzKY7xHTxol\nbdDt9QmiShW3Pi8iW9YUGUX97DPXJqRTJzf1t+++8PTT8PHH8OGH4Tm9te4mjyuugDZtXBeJ//7X\n7Wc5cKAbMcvLo1yjvTtHqSK6XVMINAUosjeFKilR/Q2/8vI+N/I6l3MDT9ObUaytekCp4WjJEmh3\nTXO+KGzNPE6kD29gQLfXJ5iITAEGp5hth4481eQxLjwnj1dfhXvvdcEO3F7EI0a4gBOOnXNeew1+\n+w3+7/92fzw9HRYsgBUroM2RG/i2WecyR3tjZT3VTjsXq6uns8guClVSsmDvofZMZyHH05QVHLd5\nJiNnH7nXD9Lt2+HJJ910S/9B1Zk0ch1Nk/9wv9WSk3V7fYIJe6gKTjFvzq1Cr00v8UbhxXyV35xO\nJ+6dnE4/3d1xOnCgtwEgJ8ctTh850t1Fu6f99oN3ng/w71V3kZ4/iec3/hNbwmhvTo4LZ6ec4l19\n4Xb44VCjBnz7rd+ViESPkEKVMeZRY8wSY8wCY8y7xpi6XhUmUaDIpsO16lbnseT7+OjBBTw+rC5d\nurgdSAB++cWtYxk3Dr76yi0QNv90t9czdar7r26vTyitWrn2HPn5YTpBTg7fVz2G1symDn8xg1NI\nq7GqxFHU++9379MRI7w5/Y4dbpuW225zG0mXxCzLoX+tLGZwCq/yL7rzPoGqjfeq8/334ZxzoFo1\nb+qLBGPUXV1kT6GOVH0MHGOtbQ4sBe4MvSSJKr13D0ct7+7E11+7f/2fdJJrw9CmjVsLMn06HHJI\nkdcW2RZHEss++7gF6/Pnh+f47y46gtM2T+AmnuQVBlCT/FKnmGvUcHee3nor/Ppr6OcfMsSNzt50\nUxlPDI72HsFSZtKOo1hC881fMCXn8N2eFit3/e1J/apEdufZhsrGmPOA8621fUr4ujZUjjM//OC2\nurnqKjjqKL+rkWhz7bUuZJcZPCpg2za480545x0YM2AyLR/q4VaHFxa6UdUyRkQff9wtYp/+doCq\nv+W40FOR0B8I8N201ZxxzTF8+VWV8jXpzMpyU37BOqfd9CF9R7Tnoovcova//nLNPlevdjPlseSv\nv1x39dWrXb87kXhSmQ2VvQxV44HR1tpRJXxdoUokgYwa5br1jxnjzfHWrHFtPpKS4M03oWFD3Nqk\nnJxyh6Pt2+HM49bQdelQbqv9HBQUlCuMAZCVReG/BtG28FMGVHmVQSPalX9ae48616+HAQPcqFmX\nLu4fKO++W75DRZuMDLhlwEa6Hf5jxUOqSBQLS6gyxkwBiu6XbgAL3G2t/SD4nLuBE62155dyHIUq\nkQSSk+MWrK9cGXqbgFmz4KKL3DqmwYPdvpWVEgiwrNmptMz/gql04AQWuuGhZctKDwOBAKSmMjj3\nNr6iDRPoiinP60phLbz8smu78PLLcMkllfx/8tnDvRfw+zszGbrPXRULqSJRrjKhqsxlkdbajmWc\n9HKgK9C+rGMNHjz474/T09NJT08v6yUiEqNSU11wWL7cfVxZI0a4tVCvvuq6lYckJ4fUmmt4Iv9m\n+jCS2bSmZvXqLgGWFo5ycvi6Smte4Crm08K1CSnP60phjLsjsWfPGN5lIBCg89gruWh7Jmy8xj3W\nvz906KARK4k52dnZZGdnh3SMkKb/jDGdgSeA062168t4rkaqRBLM+efDBRdUfuDirbfcmqxp0zza\ngDs44mRzc7mAMRzCLzyWfF+ZI065ywOcmLae++xgevGWezDEkaq4MGcOOzqcxYGbljCTdhzCr66N\nytSp7iYVkRhWmZGqUO/+ewbYB5hijJlnjHk+xOOJSBwJpV/Vhx/C9de7u8s8CVTwd5sQk5zMS/vc\nwigu4dNbxpcZjO5+OoXj2yTTK3m8eq8VlZZGlcJ8OjGZyXRyj6nRryQwzxaql3kijVSJJJyZM+G6\n62Du3Iq9bvp0uPhiF6xatw5DYcGF4xN+OoKr76zHN99AvXrFPzU72613WrgQGuyo2ML4hJCVRVbf\nSYy2F/F+9Qu1pkrihq93/5V5IoUqkYSTlwcNGsDateW/5f6rr1wjzLffdtu9hNtVV7ltLV9/fe+v\nbdrkNkh+7jno2jX8tcSqdd+v49CW9Qn8sIGkJgqbEh/8mP4TESlRzZoulMyZU8YTgxsOL8z+g3PP\ndXvqReo+lscfhxkzim9pcOONbss+BarSNfxHQ448uhozlypQSWJTqBKRsCpzXVVwY+Sl7QfRJSOP\noRfPoFu3iJVH7dpu/75rroFVq3Y9/sEHbhryiSciV0ssU3d1EYUqEQmztm1dn6liBTdGXp7bkI6b\n3+N+7uXiVzruteFwuJ18smtvcMUVrg1EIOD2sBwxAurUiWgpMatzZ5g82e8qRPylUCUiYdW2rRup\nKnZJZU4Oa6o1oQNTuYGn6c+ru/o/Rdi997qu7S+95NZZXXIJnHZaxMuIWW3auG/b6tV+VyLiH4Uq\nEQmrAw90XQh+/HHvr22ofzBnbX6XS3iTGxjiHvTplvzq1d2my7ffDt9/Dw88EPESYlq1anDmmfDx\nx35XIuIfhSoRCbvi1lVt3gxdL2tIhy5J/F/Nx6Ki/9M//uGWeI0Z4xbZS8V07qx1VZLY1FJBRMLu\n2Wddn6dhw9zneXnQrRsccoh7zKxT/6d48NtvcOKJbgqw0vszikQJtVQQkahUdKSqsNBtjpySAi++\nGNxsOSXFbWuiQBXTmjWD/feHefP8rkTEHwpVIhJ2xx/vtslbvx769nWL1keO1GhGPNIUoCQyhSoR\nCbtq1dxAVJcubmro7bfdwnCJP506qbWCJC6FKhGJiPR0N9X3/vtuPbrEp9NPd+vn/vzT70pEIk8L\n1UUkIgoLXaiqVs3vSiTcunRxjVTPP9/vSqJTfr4Lnq1a+V2JlEYL1UUkalWvrkCVKLSuqnRPPAFn\nnAEbNkTmfDt2wPbtkTlXolOoEhERT+0MVZqc2NuqVS5UtW7ttkGKhI8+cnfcSvgpVImIiKeOOMLd\n2blkid+VRJ+773ZTow89BM8/70aRwm3oUOjRI/znEdBgvIiIeMqYXaNVRx/tdzXRY+5cmDgRfvjB\nbdSdnAzTpkGHDuE753ffweLFcOGF4TuH7KKRKhER8ZxaK+zOWrjxRrj/frcjkzFu4+4XXgjveYcO\nhSuvhBo1wnsecXT3n4iIeG7jRmjaFNasgVq1/K7Gf2PGuE26583b1fT2r78gNdXdCdi0qffn3LDB\nbQW1ZAk0buz98eOd7v4TEZGoUK+e2wfw00/9rsR/eXlw663w9NO77yJQpw7885/w8svhOe8rr8A5\n5yhQRZJClYiIhIWmAJ2nn4bmzSEjY++vXXWVC1WFhd6ec9s2eO45uO46b48rpVOoEhGRsFC/Krct\n0+OPw2OPFf/1Y46Bww93Ow146YMP4MAD1WA00hSqREQkLJo3dyMwkerHFFUCAZgzh3tuyaVfPzjs\nsJKfevXVrr2Cl4YMgeuv9/aYUjaFKhERCYsqVVzjyTvugHff9buaCMrKgtRU5re/mQ/f3MQ9/xhT\n6tN79HCtD7zq6/XNN/DTT9CzpzfHk/LT3X8iIhJW8+e7qcDMTLfOKq4FApCais3NJYPp9CaLQckj\nYdkySEkp8WX33OPuBhwyJPQS+veHQw+Fu+4K/ViJTHf/iYhI1GnRAsaOhT594PPP/a4mzHJyICmJ\n9+jJH+xHf4a7jS9zckp92cCB8MYbsGVLaKdftw7eew8GDAjtOFI5ClUiIhJ27drBqFFw/vmus3g4\njB4NJ5wA+fnhOX65pKWRl2+4lcd4ihupxna3sCwtrdSXHXQQnHqqu0ahePllN51YyqCYhJFClYiI\nRESHDu6Xfrdubg3RzsXcBAIhHTc/3y32vuce97Gvo2EpKQw5bzrHVfmOM+t+7faiGT68XCnn6qvh\n+aGF2NmVuyaFhW7B+7//XZnCxQsKVSIiEjHdu7sWA2edtpVfmp0BHTu6tuJZWZU63q+/wimnuNYF\nc+e6ZpoTJ3pcdAWsXg2PTWnO4zPbwdSpbi1V797lem3HdVn89e1yvmp/Z8WvSSDA2Md+4pBmBTRv\nXsniJWRaqC4iIpEVCPBikwd4tPAGPuc0mvC7G9EpYzF30deTk8P474/gipvrceedcMMNbj+92bOh\nXz/49tvw/28UZ8AA103+8ccr+MLgAvfHc69mIceTSd/yX5OsLOjfn1MKpnFz1SH0fP3ccgc5KZkW\nqouISPTLyeHKWplcyYt0ZAoBGpZrMTcAWVkUHnQot506k2v7buL9az7mxhtdoAJo2RLWroXly8P6\nf1CsBQtc08177qnEi4ML3PvxGh9wDutoUL5rEghA//58nXs0K7YfwLkF77jb/0KcUpXKUagSEZHI\nSkuDggJu4zF68h6dmMzGguQyF3MTCLDyX/fSPu8jFhYcyTzbgraPnLdbgKhSxbVtiPQUoLVutGzw\nYKhfvxIHCF6TBvxBd97nNfqVa4H7zjD2DNdxDc+5hfHlDajiOYUqERGJrJQUt3g7OZkH6jzGqVW/\npNtBi9hSq/Rprqlv/0HL/C/oxGQm0JWGrC82QHTpEvntccaNg/Xr4YorKnmAItfkqlqZvGiuYsfL\n5VjgnpbGmvz6jOdcruAV91h5wpiEhdZUiYiIP4Jro3YclEb/O1L4/XcYPx5q1Nj9adu3w4MPwksv\nbOeNDWfTvqBIYipm3VEg4PbTW7sWkpLC/7+Rnw9HHw0vvujW3YckEMD+mkPLAc156JHqdO5c9kvu\nv3AhK8fO4aXaN7lANXy41lR5oDJrqhSqRETEd9u2Qa9esGMHvP02VKvmHg8E4NJLIS/Prcc+8FO3\nKJvq1UsNEK1bwyOPQEZG+Gt/7DHXxmH8eO+OOXy422S5rGMWFLgbBae+/QfH1PzZjVCpSZUnFKpE\nRCRm5efDeee5TPD66zBrlgtal14KDzywK2jtHOEqLUDcd58LYo88Et6a16yBY46BmTPhiCO8O+7W\nrdCsmWsTUdpM3ptvwmuvue4N4i2FKhERiWlbt7p9AqtWdQ1Chw+Hs8+u+HFmzYJBg2DhQu9rLGrQ\nIKhdG5580vtj33AD1KoF//1v8V+3Ftq0cXcbnnuu9+dPdApVIiIS8zZtcnfR/fvflV9vvX077L8/\nfPMNNG3qZXW7/PwznHwy/Pgj7Luv98f/4Qc44wy3ZGzPdWYAX34Jl1zizl+1qvfnT3TqUyUiIjGv\nbl038hPKDWxVq8JZZ8HkyZ6VtZfMTBdqwhGoAI48Eo491m2QXJyhQ+HaaxWooklIocoYc78x5htj\nzHxjzCRjTGOvChMREQlF587h61dlLYwcCX36hOf4O111ldvPb08rV7q2Ef/6V3jPLxUT6kjVo9ba\nE6y1LYCPgPs8qElERCRknTvDJ5+4mwS9NmOG6+Zw4oneH7uoc8+FX36BRYt2f/zFF90+h/Xqhff8\nUjEhhSpr7eYin9YGdoRWjoiIiDcaNYJDDnFrj7yWmQmXXbZre5xwqV4dBg6EF17Y9VheHgwbBtdd\nF95zS8WFvFDdGPMgcBnwJ5BhrV1fwvO0UF1ERCLq7rvdVF1Jd9BVRl4eHHigu7MwXIvgi1q50q2t\nWr4c6tRx7SbeeivyW/EkmrAsVDfGTDHGLCzyZ1Hwv+cAWGvvsdYeBLwJKDeLiEjU6NLF+/DxwQdu\n2i8SgQqgSRM480x44w0XEIcOdXdGSvSpVtYTrLXlbbo/CpgADC7pCYMH7/pSeno66enp5Ty0iIhI\nxZ18susTumoVHHCAN8fcOfUXSVdfDddf70astmxxm0aLt7Kzs8nOzg7pGCFN/xljDrPW/hT8+Drg\nNGvtRSU8V9N/IiIScRdeCN26weWXh36stWtd5/QVK2CffUI/XnlZC0cd5T6+9lr3R8LLjz5VDwen\nAhcAHYDrQzyeiIiIp7ycAhw9Gs45J7KBCtyC+KuuciNufftG9txSfuqoLiIice3339202dq1RfYP\nrKSWLeF//4OO5V0Y46G8PFi82NUg4aeO6iIiIns48EBodsA2Zo9Y4jZjrqTvvnMjRe3be1hcBdSs\nqUAV7RSqREQkvmVl0eXHp5l4zYeQmgpZWZU6zMiRblsabQsjJdH0n4iIxK9AAFJT+TS3FbfwOHNo\n7VqhL1sGKSnlPsz27W4vwokT3VSixD9N/4mIiBSVkwNJSbRjJks5nLWkuDblOTkVOkx2tstgClRS\nGoUqERGJX2lpUFBAdbbRnmlMppPbDDAtrUKHycwM/+bJEvsUqkREJH6lpMDw4ZCcTJea2Uyserb7\nvAJTf1u2wPvvQ+/eYaxT4oLWVImISPwLBPjtq99p0fc41qytUv7F5oEAb7zwF1mfN+GjKTXCWqJE\nF62pEhERKU5KCs3OPoHGB1Th66/L+ZqsLEhNJfOBHPp8OqDSdw1K4tBIlYiIJIxbb4XataHIVrTF\nC941uDJ3X45jEStpQnKyqfBdgxK7NFIlIiJSii5dYNKkcjwxeNfgKP5JT94jmbxK3TUoiUWhSkRE\nEsYpp8CSJbBuXRlPTEvD5heQyWX0YaR7rBJ3DUpiUagSEZGEUaMGpKfDxx+X8cSUFL75v3f5y9Tl\ntDrfuIahFbxrUBKPQpWIiCSU8k4BZga60OfGhlT5ZIpbS6WeClIGLVQXEZGEkpMDbdq4zZGrlDC0\nsG0bNG0Kn30GRxwR0fIkSmihuoiISBnS0mC//WDevJKfM2WKe54ClVSEQpWIiCScsqYAMzPhsssi\nV4/EB4UqERFJOJ07w8SJxX9t40aYMAEuvjiyNUnsU6gSEZGEc/rpsGgR/PHH3l97911o3x4aNIh8\nXRLbFKpERCTh1KzpgtXUqXt/TVN/UlkKVSIikpCKmwLMyYHFi6FrV19KkhinUCUiIglp52L1HTt2\nPfbGG3DRRa5JqEhFKVSJiEhCOvRQqFMHFi50n1sLI0dq6k8qT6FKREQSVtEpwNmzXbBq08bfmiR2\nKVSJiEjC6tJlV6jKzIQ+fcBUqIe2yC7apkZERBJWbi40agQ//wxHH+1Gqw4+2O+qJBpUZpsahSoR\nEUloXbpAw4awfDl8+qnf1Ui00N5/IiIiFdS5s7vrr08fvyuRWKeRKhERSWhLl8JJJ7mRqvr1/a5G\nooWm/0RERCph82bYZx+/q5BoolAlIiIi4gGtqRIRERHxiUKViIiIiAcUqkREREQ8oFAlIiIi4gGF\nKhEREREPKFSJiIiIeEChSkRERMQDClUiIiIiHlCoEhEREfGAJ6HKGHOzMWaHMWY/L44nIiIiEmtC\nDlXGmKZAR2BZ6OUknuzsbL9LiDq6JsXTdSmerkvxdF32pmtSPF0X73gxUvUUcKsHx0lIejPvTdek\neLouxdN1KZ6uy950TYqn6+KdkEKVMeZc4Ddr7SKP6hERERGJSdXKeoIxZgrQqOhDgAXuAe7CTf0V\n/ZqIiIhIwjHW2sq90JhjganAVlyYagqsBFpba9cW8/zKnUhERETEB9baCg0WVTpU7XUgY34FTrTW\nbvDkgCIiIiIxxMs+VRZN/4mIiEiC8mykSkRERCSRRbSjujHmPmPMCmPMvOCfzpE8fzQxxnQ2xnxv\njPnRGHO73/VEC2NMjjHmG2PMfGPMbL/r8YsxZrgxZo0xZmGRx/Y1xnxsjPnBGDPZGFPPzxr9UMJ1\nSeifK8aYpsaYacaYb40xi4wx/w4+ntDvl2Kuy3XBxxP9/VLDGPNV8GfsImPMfcHH04wxXwZ/J2UZ\nY8q8kS1elHJNXjPG/BJ8fJ4x5vgyjxXJkapgoX9Za5+M2EmjkDGmCvAjcCbwOzAH6GWt/d7XwqKA\nMeYX4KREX5tnjDkV2AxkWmuPDz72CLDeWvtoMIjva629w886I62E65LQP1eMMY2BxtbaBcaYfYC5\nQHegHwn8finlulxMAr9fAIwxtay1W40xVYEZwPXATcAYa+07xpgXgAXW2pd8LTSCSrgmVwIfWGvf\nK+9x/Nj7T+uuoDWw1Fq7zFpbCIzG/WUX9/5I+D0prbVfAHsGy+7AiODHI4DzIlpUFCjhukAC/1yx\n1q621i4IfrwZWIK7Gzuh3y8lXJcmwS8n7PsFwFq7NfhhDVxrJQtkAO8GHx8B9PChNN8Uc012BD+v\n0HvFj19e1xhjFhhjXkm04egimgC/Ffl8Bbv+sic6C0w2xswxxgzwu5gos7+1dg24XxjA/j7XE030\ncwU3hQM0B74EGun94hS5Ll8FH0ro94sxpooxZj6wGpgC/Az8aa3dGSRWAAf6VZ8f9rwm1to5wS89\nGHyvPGGMqV7WcTwPVcaYKcaYhUX+LAr+9xzgeeBQa23zYOEJO/wqJTrFWtsS6Ir7wXeq3wVFMd1l\n4ujnChCc4hoDXB8cmdnz/ZGQ75dirkvCv1+stTustS1wI5qtgX/4XJLv9rwmxpijgTustUcBrYAG\nQJnrnz1fiGat7Vj2swB4GfjA6/PHiJXAQUU+39k4NeFZa1cF/xswxozF/YX/wt+qosYaY0wja+2a\n4HqRvZrsJiJrbaDIpwn5cyW4qHgMMNJa+37w4YR/vxR3XfR+2cVau8kYkw20BeobY6oER6sS9ndS\nkWvSeee6O2ttoTHmNeDmsl4f6bv/Ghf5tCewOJLnjyJzgMOMManGmCSgFzDe55p8Z4ypFfxXJcaY\n2sBZJO57BNxcftH5/PHA5cGP+wLv7/mCBLHbddHPFQBeBb6z1g4p8pjeL8Vcl0R/vxhjGu6c8jTG\nJOO2mvsOmA5cGHxaQr1fSrgm3+98rxhjDG5NYpnvlUjf/ZeJm9feAeQAg3bO+Sea4G28Q3DBdri1\n9mGfS/KdMeZgYCxumqIa8GaiXhdjzCggHTfkvAa4DxgHvAM0A5YBF1lr//SrRj+UcF0ySOCfK8aY\nU4DPgEW4vzsWty/rbOBtEvT9Usp1+SeJ/X45DrcQvUrwz1vW2oeCP39HA/sC84FLgzdSxb1Srskn\nQEPcP+IWAFcWWdBe/LHU/FNEREQkdAl/67qIiIiIFxSqRERERDygUCUiIiLiAYUqEREREQ8oVImI\niIh4QKFKRERExAMKVSIiIiIeUKgSERER8cD/A2e/YnBuy4hZAAAAAElFTkSuQmCC\n",
       "text": [
        "<matplotlib.figure.Figure at 0x7f7cf732c6d0>"
       ]
      }
     ],
     "prompt_number": 2
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The blue profile is a collection of $N=30$ patches. The green dots above are lamps. The goal of this exercize is to set the power of each lamp, in order to obtain an intensity *as uniform as possible* on all patches. The quantity of interest is in fact $\\log I$. We propose to minimize\n",
      "$$\\max_{k=1,\\ldots,N}\\ | \\log I_k - \\log I_{des} | ,$$\n",
      "subject to the constraint that the power of each lamp must be in the interval $[0,p_{max}]$.\n",
      "\n",
      "Physics tell us that the intesity on the $i$th patch is given by\n",
      "$ I_i = \\sum_{j=1}^M A_{i,j}\\, p_j$,\n",
      "where $p_j$ is the power of the $j$th lamp and\n",
      "$$A_{i,j} = \\left\\{\\begin{array}{cl}\n",
      "                    \\frac{1}{(r_i,j)^2} \\cos \\theta_{i,j} & \\textrm{if patch $i$ receives light from lamp $j$}\\\\\n",
      "                    0 & \\textrm{otherwise.}\n",
      "                    \\end{array} \\right.                    \n",
      "$$\n",
      "In this formula, $r_{i,j}$ denotes the distance from patch $i$ to lamp $j$, and $\\theta_{i,j}$ is the illumination angle."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "--"
     ]
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "Computation of Matrix $A$"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "#r[i,j] is the distance from lamp j to patch i\n",
      "r = np.zeros((N,M))\n",
      "for i in range(N):\n",
      "    for j in range(M):\n",
      "        r[i,j] = np.linalg.norm([x[i]-xlamps[j],y[i]-ylamps[j]])\n",
      "\n",
      "\n",
      "#theta[i,j] is illumination angle from lamp j to patch i\n",
      "costheta = np.zeros((N,M))\n",
      "for i in range(N):\n",
      "    for j in range(M):\n",
      "        #index of xxx_k's between xi and xlamp_j\n",
      "        index_x_between = [k for k in range(N+1) if np.sign((xxx[k]-x[i])*(xxx[k]-xlamps[j])) == -1]\n",
      "        y_patch = [yyy[k] for k in index_x_between]\n",
      "        y_line_from_lamp = [ymax +2. + (xxx[k]-xlamps[j]) * (y[i]-ymax-2.)/(x[i]-xlamps[j]) for k in index_x_between]\n",
      "\n",
      "        if any([y2<=y1 for y1,y2 in zip(y_patch,y_line_from_lamp)]):\n",
      "            #this patch gets no line from this lamp (it is hidden by some other patch)\n",
      "            costheta[i,j] = 0.\n",
      "        else:\n",
      "            dy = yyy[i+1]-yyy[i]\n",
      "            normal_dir = np.array([-dy,1.])\n",
      "            normal_dir /= np.linalg.norm(normal_dir)\n",
      "            lamp_dir = np.array([xlamps[j]-x[i],ylamps[j]-y[i]])\n",
      "            costheta[i,j] = (normal_dir.dot(lamp_dir))/np.linalg.norm(lamp_dir)\n",
      "            \n",
      "assert np.min(costheta)>=0.\n",
      "assert np.min(costheta)<=1.\n",
      "\n",
      "A = cvx.matrix(0.,(N,M))\n",
      "for i in range(N):\n",
      "    for j in range(M):\n",
      "        A[i,j] = 1./r[i,j]**2 * costheta[i,j]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 3
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "short function to visualize a solution"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def visualize_sol(pp,normalized=False):\n",
      "    colors = (np.log(A*pp)-np.log(Ides)).ravel()\n",
      "    if normalized:\n",
      "        mm,MM = -1.2,1.2\n",
      "    else:\n",
      "        mm,MM = min(colors),max(colors)\n",
      "    logdifs = np.array(colors)\n",
      "    colors -= mm\n",
      "    colors /= (MM-mm)\n",
      "    \n",
      "    for i in range(N):\n",
      "        plt.plot([xxx[i],xxx[i+1]],[yyy[i],yyy[i+1]],color = matplotlib.cm.jet(colors[i]),linewidth=2)\n",
      "    fg = plt.scatter(x,y,c=logdifs)\n",
      "    plt.scatter(xlamps,ylamps, s=20*pp)\n",
      "    cbar = plt.colorbar(fg)\n",
      "    yt = cbar.ax.get_yticks()\n",
      "    cbar.ax.set_yticklabels(['Ides * {0:.2f}'.format(yi) if yi>=1 else 'Ides / {0:.2f}'.format(1./yi) for yi in exp((MM-mm)*yt+mm)])\n",
      "    \n",
      "    II = A*pp\n",
      "    print 'des. Intensity: '+str(Ides)\n",
      "    print 'min  Intensity: '+str(min(II))\n",
      "    print 'max  Intensity: '+str(max(II))\n",
      "    print 'ratio Imax/Ides:'+str(max(II)/Ides)\n",
      "    print 'ratio Ides/Imin:'+str(Ides/min(II))\n",
      "    "
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 4
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Let us now visualize the solution with all lamps on, with a power p=1:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "p_unit = cvx.matrix([1.]*M)\n",
      "print p_unit\n",
      "visualize_sol(p_unit)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "[ 1.00e+00]\n",
        "[ 1.00e+00]\n",
        "[ 1.00e+00]\n",
        "[ 1.00e+00]\n",
        "[ 1.00e+00]\n",
        "[ 1.00e+00]\n",
        "[ 1.00e+00]\n",
        "[ 1.00e+00]\n",
        "[ 1.00e+00]\n",
        "[ 1.00e+00]\n",
        "\n",
        "des. Intensity: 0.1"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "min  Intensity: 0.0545211480435\n",
        "max  Intensity: 0.2531736417\n",
        "ratio Imax/Ides:2.531736417\n",
        "ratio Ides/Imin:1.83415066609\n"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAkYAAAFwCAYAAABdFrRBAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3XeYVOXZx/HvvX2X3kWKgCCCSLNhw7WiL4iKBnsLUWOI\nJhpjI0ZsidHErlETxW5i7KJYcW0IooJ0BBRQQBGQurBl5n7/mNllxa2cKbvD73Ndc+3smXOe555h\n2XPvU83dERERERFIS3YAIiIiIvWFEiMRERGRKCVGIiIiIlFKjERERESilBiJiIiIRCkxEhEREYnK\nCFqAmWUD7wNZ0fKedffrgpYrIiIikmgWi3WMzCzP3QvNLB34CLjY3T8JXLCIiIhIAsWkK83dC6NP\ns4m0GmnVSBEREWlwYpIYmVmamU0DvgPecvepsShXREREJJFi1WIUdvcBQEdgPzPrHYtyRURERBIp\n8ODritx9vZm9CxwNzKn4mpmpe01ERHYo7m6Jqqu5ma8LXswSd+8SvJiGKxaz0loDJe6+zsxygSOB\nmys7N5U3rB07dixjx45Ndhhxk8rvL5XfG+j9NXR6fw2XWcJyIgDWAWMDljEWdgkeScMWixaj9sCj\nZpZGpGvuv+7+WgzKFRERkTqIaTfQDirwZ+juM4GBMYhFREREAshMdgApQMlljOTn5yc7hLhK5feX\nyu8N9P4aOr0/qQvd1IOLyQKPtarIzFN5jJGIiEhFZpbQwddm5v8IWMYfSOyA8fpIyaWIiEiKUFda\ncEqMREREUoRu6sHpMxQREUkRajEKLiYrX4uIiIikArUYiYiIpAjd1IPTZygiIpIi1JUWnBIjERGR\nFKGbenAaYyQiIiISpeRSREQkRagrLTglRiIiIilCiVFwSoxERERShG7qwekzFBERSRFqMQpOg69F\nREREotRiJCIikiJ0Uw9On6GIiEiKUFdacEqMREREUoRu6sFpjJGIiIhIlBIjERGRFJEZ8FEZM9tQ\nxfFxZjYiFnGb2bVVHO9nZpPMbKaZTTezkVWcd4mZzY6e85aZdYoe72xmn5nZ59EyLqgpFrW6iYiI\npIg43dQ9PsWCmd0ETAFamdkdwMPuPqPCKZuAM919kZm1Bz4zs9fdff02RX0O7OXuW8zs18CtwCnA\nCmCQu5eYWR4w28xecvfvqopJiZGIiEiKiPfgazO7Bzgc+AYoqXB8IHAb0AhYBZzj7t+b2cXABdFz\n57j7aRXLc/cxZnYfcDKRBGbBNq8vrPB8hZmtBNoA67c5770K304GTo8eL6lwPBewmt6jEiMRERGp\nUbTbrIe794q23swBHjKzDOBuYLi7r452d/0FGAVcAXSJttg0raTMG4AJQCkw2swecveZVdS/L5Dp\n7otqCHVUtMyy6zoCrwK7An+srrUIlBiJiIikjDjf1A8Gnoby1puJ0eM9gT7AW2ZmRMYvL4++9gXw\nlJm9CLy4bYHufg2AmQ1w9+urqjiaiD0GnFldgGZ2BrAXcEiFOr4F+pnZTsBLZvasu/9QVRlKjERE\nRFJEXbvSJkcfARkwy90PrOS1ocBgYDgwxsz6uHt425NqSIqaAOOBq9x9ajXnHQFcBQzepgutrI7v\nzGwWkQTv+arK0aw0ERGRFFHXWWgHA3+s8KhC2bic94GTzSwt2oJzaPT4fKCNmQ0CMLMMM+sdfa1z\ndPzPlUBToHFd3o+ZZRJpaXrU3V+o5rwBwP1Eu/MqHO9gZjnR5y2Ag6LxVkktRiIiIlIdB3D3F8zs\nMGA2sBSYFD1eYmYnAXebWTMgHbjDzL4EnoiOLTLgzkpmk9VkJJFkpoWZnRuN5Rx3n2Fm1wFT3X08\ncAuRgd//i3bnLXH344FewD/MLByN4RZ3n11dheYet1l4P63IzBNVl4iISLKZGe5e4yyoGNZX/aji\nWtgJEhpzfaQWIxERkRSRGfSuXhqTMBo0JUYiIiIpIkOJUWAafC0iIiISpRYjERGRFJGZnuwIGj4l\nRiIiIikicFeaqCstFS1btozzz7+IYcNO5V//eojazAZ88smnOe640zn33AtZtKim1dZTz8qVK/nN\nby5h2LBTufvu+wiHf7b+2A5p+fLlXHDBxQwbdioPPPCvWv0sVWbz5s2MGTOWoUNP4ZprrmfLli3l\nr7k7jzzyGMceexqjRo1myZIlsQo/ad544w1GjDiL0077FdOmTav23GeffY7jjz+Ds866gPnzq11e\npdxzzz1ffs28efOqPC8cDnPXXfcybNip/OY3l7By5co6vQ9peDIzgj2EyC+lRDwiVUm8/fDDD96m\nTWdPT7/c4TFv1Kifjxkzttprbr/9Ls/L283hEU9Lu9abNdvJly5dmqCIk2/dunXeoUMPz8z8ncPj\nnpe3n48efWmyw0q61atXe9u2XTwj4zKHxzwvb4Bffvmf6lxOKBTyQYMO95yckxye8Jyc433w4GM8\nHA67u/uNN/7N8/J6Ozzq6eljvGXLDr5ixYpYv52EeeGFFzwvb2eHBxz+7o0atfZp06ZVeu6DD/7b\n8/K6OjzsZjd6kyZtfdGiRdWW/+9/P+x5eV2i19zkTZq08QULFlR67ujRl3pe3n4Oj3tm5sXeoUMP\nX7duXeD3KLUTve8l9j7bNtgj0THXx0di/8Ek7h544AHPzT3ZwaOPxZ6b27zaa9q27eYwrfyazMwL\n/eabb05QxMn39NNPe+PGx1T4zFZ5Rka2l5SUJDu0pHrooYc8L+/ECp/Lt56d3bg8oamtGTNmeKNG\n3RxKo+UUe15eJ583b567uzdr1t5hXnk92dnn+J133hmPt5QQAwbkO7xY4XP7q5977oWVntu5cx+H\nj8rPTUv7g48Zc0215Xfp0tfhgwrX/NGvvHLMz84rKSnxjIxsh1Xl5zZufIw//fTTMXmfUjMlRg3z\noYazFFNaWop7boUjuYTD1c+/DIVKgZzy78PhHEpKdpw5mz//zLIr/qLZYVX+sxTarnLMstjac5+O\nWSalpZGfsW1//txzyl9riCKx//RzKy6u/P38/P9ebo3/9yLlb3tN0c/O2/ozXPGzzW3Qn63Ugu7q\ngWmMUYoZPnw4mZmvYXYXMJG8vFM466xzq73m/PPPJS/vLOAt4EFych7nF784KRHh1gtDhgwhO3sK\naWm3ABPJzT2ZE088lczMum7HmFqGDRtGVtbbmN1O5GdpJGeccQ6R1fZrr0+fPnTs2IjMzIuBd8nK\n+g1du7Zh9913B2DUqHPJyzsdeBu4j6ys5zjhhBNi/XYS5uKLzyUv77fAq8DT5Ob+hQsuqHxD8Asv\nPJe8vF8CbwDjyMu7n9NPP6Xa8kePPpe8vFHRax4hL+8+zjzz1J+dl5mZyYknnkpu7khgImlpt5Cd\nPYUhQ4YEe4NSv2UEfIi60lLRrFmzfMiQE71fv8H+pz9dX2OXUCgU8ptv/ocPGJDvhx463KdOnZqg\nSOuPBQsW+NChI71v34P9ssvGeFFRUbJDqhfmzJlT/rN09dVjvbi4eLvKWbVqlZ9++nm+554H+Vln\nXeBr1qwpf620tNSvv/6v3r//IX7EEcf79OnTYxV+0jz00Djfe+/D/cADj/E33nijyvPC4bDfdttd\nPnDgoX7IIcP8o48+qrHscDjsd955jw8ceKgPHjzMP/zwwyrPLSoq8ssuG+N9+x7sQ4eOrHIsksQH\nyehK6xbskeiY6+Mj8F5pZtYReAxoB4SBf7n7XZWc50HrEhERaSiSsVeadwtYxlfaKy0WDWelwKXu\nPt3MGgOfmdmb7l71HFIRERGJPS3wGFjgxMjdvwO+iz7faGZzgQ6AEiMREZFE0jihwGL6EZpZF6A/\nMCWW5YqIiEgtKDEKLGaz0qLdaM8Cv3P3jbEqV0RERCRRYpJbmlkGkaTocXd/qarzxo4dW/48Pz+f\n/Pz8WFQvIiKSdAUFBRQUFCQ3CI0xCizwrDQAM3sMWOXul1ZzjmaliYjIDiMps9L2DljGp5qVFrjF\nyMwOBE4HZprZNMCBq9399aBli4iISB1ojFFgsZiV9hFqvBMREZEUoNxSREQkVaiZIjAlRiIiIqlC\nd/XA9BGKiIikCt3VA9NHKCIikip0Vw8sZgs8ioiIiDR0yi1FRERShQZfB6bESEREJFXorh6YPkIR\nEZFUobt6YBpjJCIiIhKl3FJERCRVaIxRYEqMREREUoXu6oHpIxQREUkVuqsHpjFGIiIiIlHKLUVE\nRFKFxhgFpsRIREQkVeiuHpi60kRERFJFRsBHJcxsQxXHx5nZiFiEbWbXVvPaBDP70cxeruacTmY2\n0cw+N7PpZnZM9PgRZvapmX1hZlPN7NCaYlFiJCIiItXxeBVsZjeZ2XCglZndYWZ9KzntFuCMGor6\nE/Bfdx8InArcFz3+AzDM3fsB5wCP1xSTEiMREZFUEYcWo4rM7B4zm2tmbwJtKxwfaGYF0VaZCWbW\nLnr8YjObHW3FeWrb8tx9DHA0cDpwr7vPqOScd4GNNYQWBppGnzcHlkWv/cLdv4s+nw3kmFlmdQWp\nN1JERCRVxHHwdbTbrIe79zKz9sAc4CEzywDuBoa7+2ozGwn8BRgFXAF0cfcSM2taSZk3ABOAUmC0\nmT3k7jO3I7zrgDfN7GIgDziikrpOAj5395LqClJiJCIikirie1c/GHgawN1XmNnE6PGeQB/gLTMz\nIr1Ry6OvfQE8ZWYvAi9uW6C7XwNgZgPc/foAsZ0KjHP3281sEPAEsEfZi2a2B/BX4MiaClJiJCIi\nkirqeFcvWBp5BGTALHc/sJLXhgKDgeHAGDPr4+7hbU8KmBRBpHVqSLSsyWaWY2at3X2VmXUEngfO\ndPfFNRWkMUYiIiI7qPzOMPagrY8qWPTr+8DJZpYW7Uorm+E1H2gTbanBzDLMrHf0tc7u/h5wJZEx\nQI23M1SrEEdllhDtPjOzXkB2NClqDowHrnD3ybWpSC1GIiIiqSI+Y4wcwN1fMLPDgNnAUmBS9HhJ\ndPzO3WbWLBrFHWb2JfBEdGyRAXe6+/q6Vm5m7xPprmtsZkuBUe7+lpldB0x19/HAZcC/zOwSIgOx\nz45ePhrYFfhzdEkAB45y91VV1ucet1l4P63IzBNVl4iISLKZGe5eXStHrOtz/3PAMq4noTHXR2ox\nEhERSRW6qwemMUYiIiIiUcotRUREUoXu6oHpIxQREUkVcVzgcUehxEhERCRV6K4emMYYiYiIiEQp\ntxQREUkVuqsHpo9QREQkVWiMUWBKjERERFKF7uqBaYyRiIiISJRySxERkVShu3pg+ghFRERShcYY\nBabESEREJFXorh6YPkIREZFUobt6YBp8LSIiIhIVk8TIzB4ys+/NbEYsyhMREZHtkBHwITFrMRoH\nDIlRWSIiIrI90gM+JDb5obt/aGa7xKIsERER2U5q9QlMY4xEREREopRbioiIpArd1QNL6Ec4duzY\n8uf5+fnk5+cnsnoREZG4KSgooKCgILlBaJxQYObusSnIrAvwirvvWcXrHqu6RERE6jszw90tgfW5\nvxOwjMNJaMz1Uaym6z8FTAJ2M7OlZnZuLMoVERERSaRYzUo7LRbliIiISAAaYxSYPkIREZFUoTFG\ngSkxEhERSRW6qwemdYxEREREopRbioiIpArd1QPTRygiIpIqdFcPTB+hiIhIqtDg68CUGImIiKQK\n3dUD0+BrERERkSjlliIiIqlCd/XA9BGKiIikCo0xCkyJkYiISKrQXT0wjTESERGRKpnZhiqOjzOz\nETGq49pqXptgZj+a2cu1KOdEMwub2cDo95lm9rCZzTCzaWZ2SE1lKLcUERFJFfG5q3tcSgXM7CZg\nCtDKzO4AHnb3GducdguQB1xQQ1mNgYuByRUOnwe4u/c1szbABGDv6spRi5GIiEiqyAj4qIGZ3WNm\nc83sTaBtheMDzazAzKZGW3jaRY9fbGazzWy6mT21bXnuPgY4GjgduLeSpAh3fxfYWIt3fwNwM1BU\n4VhvYGK0nB+AtWamxEhERGRH4OnBHtWJdpv1cPdewNnAAdHjGcDdwInuvg8wDvhL9LIrgP7u3h/4\ndSVl3kCkFedJYLSZ7bk979vMBgAd3X3CNi99AQw3s3Qz6wrsBXSqrix1pYmIiOygCj6Ego9qffrB\nwNMA7r7CzCZGj/cE+gBvmZkRaXRZHn3tC+ApM3sReHHbAt39GogkNu5+/fa8h2idtxFJ1soPR78+\nDPQCpgJLgI+AUHXlKTESERFJEaE63tUPzo88ylx/y3ZVa8Asdz+wkteGAoOB4cAYM+vj7uFtT9re\npCiqCbAHUBBNknYCXjKz4e7+OXBpeaBmHwFfVleYutJERERSRCgj2KMKZa0v7wMnm1mambUHDo0e\nnw+0MbNBEOlaM7Pe0dc6u/t7wJVAU6Dxdr41qxDHT7j7endv6+7d3L0rkcHXx7r752aWa2Z50biO\nBErcfV51FanFSEREJEWUpgdt7/hZYw5EZ6W5+wtmdhgwG1gKTIoeLzGzk4C7zawZkWUm7zCzL4En\nzKwpkaTmTndfX9eIzOx9It11jc1sKTDK3d8ys+uAqe4+vpJ4y5KotsAbZhYClgFn1life9xm4f20\nIjNPVF0iIiLJZma4e6WtHHGqzzduCZYYNc4JJzTm+kgtRiIiIikilBH0tl4ckzgaMiVGIiIiKSKU\nrs3SglJiJCIikiJC2kU2MCVGIiIiKaJUiVFgmq4vIiIiEqUWIxERkRQR0m09MH2CIiIiKUJjjIJT\nYiQiIpIilBgFpzFGIiIiIlFqMRIREUkRajEKTomRiIhIitB0/eCUGImIiKQIzUoLTmOMRERERKKU\nWoqIiKQIjTEKTomR7BCKPMSbbzxPz6Zt2O2A/GSHIyISF0qMglNiJDuEWY89RNtzLmDGvn3pMXk6\nZpbskEREYk6Dr4PTGCPZIXxx0hDWtWlJp09msPa115IdjoiI1FNKjGSH8FkjY/wVvwRg2bXX4u5J\njkhEJPZCZAR6iBIj2QE4zmds5u0LTyatXTs2ffYZa195JdlhiYjEXIj0QA+JUWJkZkeb2Twz+9LM\nrohFmSKxspQSVlFK47zGdLzqSgC+VauRiKQgJUbBBU6MzCwNuAcYAuwBnGpmuwctVyRWPqMQgIHk\n0e78C8hs357C6dP58cUXkxyZiEhsKTEKLhYtRvsCC9x9ibuXAP8BjotBuSIxUZYY7UUeabm57Hz1\n1UC01SgcTmZoIiJSz8QiMeoAfFPh+2+jx1LOG8vggo9BPTANS8XECKDtr35FZocObJ45kzXPP5/M\n0EREYqqU9EAPSfA6RmPHji1/np+fT35+fiKrD2RdMZzyAawthm5N4Io+yY5IamMtpSykiCyMPcgB\nIC0nhw5jxrD4N79h2dixtBwxAkvTPAQRCaagoICCgoKkxqCZZcFZ0AGoZjYIGOvuR0e/vxJwd//b\nNud5Qx/s+vI3cNy7YMCEw2FISraLpZaJbOB8lrI3efyHruXHw0VFTOncmbSVK7nWjE0DB/L0iy/S\nsWPHJEYrIqnEzHD3hK0ma2Y+3g8PVMYweyehMddHsfgzeSrQ3cx2MbMs4BTg5RiUW+8M7wTX9gUn\n0nq0aEOyI5KaVBx4XdGS5cu5Z+1aAK51p9e0aRw/ZEjC4xMRkfolcGLk7iHgt8CbwGzgP+4+N2i5\n9dWf+8HwjpEutePfhY0lyY5IqvP5NuOLykyaNIkfMzOByH+CX4bDzP7ySzZsULYrIg2XZqUFF5PO\nSHd/HegZi7LquzSDxw+CfV+DWWvh9IJifp8xiVatWrLnnntqD656pJgwM9gMwEByf/JamzZt+MqM\nEJAOlALpZuTl5f2sHBGRhkIDqIPTiNPt0DQLXjwUGqWFeHlFFsc8MY/9jxzGmb88X4sG1iOz2UIR\nzq5k02KbvwGOOOIIeh5wAL9t1Ii1QCYw7rTTSE/XLxURabi0JUhwSoy20+7NoPH/LgGg6JBfU3jT\nl7z0wae8qEUD642t0/Rzf/ZaWloaz7/2Glc/8girTz0VgB4ffoiXliY0RhERqV+UGAWw5q2HoGRL\n5JusHLYM/AULFy5MblBSbtv1i7aVnp7OSSedxKmPPUZ29+4ULVrEqscfT2SIIiIxpTFGwSkxCqBn\n775Q8GD596F9TqJHnwFJjEjKRDaOrT4xKmMZGXS89loAlt1wA+ESjagXkYZJiVFwSowCeP7px9j5\nw7tpNLYffDcfb7cbT+YdoZWx64ElFLOGEK1IZxeyajy/1amnktOzJ0Vff82qRx9NQIQiIrGnla+D\nU2IUQI8ePVj85Rw+f/1/TPpFG5pmwrNL4OZZyY5MKrYWGTXPFLT0dDpUbDUqLo5rfCIiUj8pMQoo\nMzOT3Xbbjf27tuTJgyOrYo+ZBhOWJTuyHdunVSzsWJ1WI0eS06sXxUuX8sO4cfEKTUQkbuIxK83M\nKl3gzczGmdmIWMRtZtdW89oEM/vRzKpcPNrMLjGz2WY23czeMrNO0eP9zGySmc2MvjaypliUGMXQ\nsI5wXf/IytinfQAL1yc7oh1XVQs7VsfS0+kY3c9v+Y03Ei4qikdoIiJxE6cxRnEbIGJmN5nZcKCV\nmd1hZn0rOe0W4Iwaivoc2Mvd+wPPAbdGjxcCZ7r7nsAxwB1m1rS6gpQYxdiYPeH4TnBwpxd4Y/0f\nWV2yONkh7XB+pJRFFJNdYePY2mp50knk9ulDeu63rH3p0jhFKCISH/EefG1m95jZXDN7E2hb4fhA\nMysws6nRFp520eMXV2jJeWrb8tx9DHA0cDpwr7vPqOScd4GN1cXl7u+5e3SaOJOBDtHjC9x9UfT5\nCmAl0Ka6spQYxViawWMHwbCu73JAx1d5acUXyQ5phzJ+/Hj2+d15AOQtXAbFdVuXyNLS2OXW09nz\nLWi86/2EN6+NR5giInERz8Qo2m3Ww917AWcDB0SPZwB3Aye6+z7AOOAv0cuuAPpHW3J+XUmZNwAT\ngCeB0Wa2Zww+hlHRMreta18gsyxRqooSozhokgkHtOwDwBGdNBI7UaZOncopo86l5HenALD80xlc\n9Mc/1LmcpkdexpYleWS1DeMb7o91mCIiDdXBwNNQ3voyMXq8J9AHeMvMpgFjgJ2jr30BPGVmpwOh\nbQt092vc/RVglbv/3t1nBgnQzM4A9mJrV1rZ8fbAY8A5NZWh9b/jpHNOHxYBW2x2skPZYbw8/hWy\nzxtBZreOAGQeMpDnL7+dB+68u07lWHoGOQOfgB9HkF5yB4QvhjTtoSYi9V9dp9wvKFjOwoLlQas1\nYJa7H1jJa0OBwcBwYIyZ9XH38LYnufv1gYMwOwK4Chjs7iUVjjcBxgNXufvUmspRi1Gc5NEbgELm\n4mibiURo2rgJfLe6/PuSr74lr3Hj7SrLGh0PmXtD+HvYeF+sQhQRiau6zkLrlt+Zo8YOKn9UoWzN\nk/eBk80sLdoCc2j0+HygjZkNgkjXmpn1jr7W2d3fA64EmgLb90s5EkOVa6+Y2QDgfmC4u6+ucDwT\neBF41N1fqE1FSoziJIPmZNEJp4jNaJuQRPjlL39J3rpNAIRWrmH9L/7Irddt5x8hZtAseu2Gv0G4\n2nF/IiL1QjxnpUUTi4XAbOARYFL0eAlwEvA3M5sOTAP2j449esLMvgA+A+509zrP1zaz94H/AoeZ\n2VIzOzJ6/DozGxY97RagEfA/M5tmZmUbl44EDgLOiR7/vIqZb1vrS9Ru8GbmO9rO81/zB9Yygc7c\nQCtOTHY4O4Q7Ni3hnkYb6TLlS8aWtOGggw7a/sLcYeUBUDwZmv0Fml4Vu0BFJOWZGe5e8wqzsavP\nb/GLApVxud2d0JjrI7UYxVEeewBQiAZgJ8q8RpH/z6P3OzRYUgTRVqMbIs833AphLUwlIvWb9koL\nTolRHOURmZlWiAZgJ0LFjWPrsuJ1tbIPh+yDIfwjbLgzNmWKiMSJ9koLTolRHEUGYBubmUcY7b0V\nb19TzI+EaE0GncmMTaFm0LSs1egfhEpWxaZcEZE4iMeWIDsaJUZxlE5jsumKU8oWvkx2OClv68ax\nubXaOLbWcg5hffE+4Ov4y1Vt6LlbR6ZOrXHGp4iINEBKjOJM44wS57Pt2B+tNjZt2sSZo+cB8KeL\n4ZbLlzH82CNZv15jjkSkftEYo+CUGMWZxhklTrwSo/nz5/P1N5EWKDM4bgjs3M6ZO3duTOsREQlK\niVFw6lCMs7LEaJNajOJqNaV8TTE5GL3JjWnZbdq0YdmKYtaug+bNIrP4CwuLaNu2bc0Xi4gkkJKb\n4NRiFGd57A6ksYWFhNmc7HBS1ufR1qK+5JIZy/FFQKdOnfjtb3/HXsc24ou5aZjBg7f1oWvXrjGt\nR0REkk+JUZylkUsO3YEQm5mf7HBSVllitHeMu9HKXHf9zTz8yKt8+cMVABw8YDaEVsSlLhGR7aXp\n+sGpKy0B8ujDFr6kkFk0on+yw0lJ8RpfVNEhhxwCHAKr5sHmF2D9zdBCaxuJSP2hKffBqcUoATTO\nKL6KCDOTLRgwII6JUbmmYyNfNz4Apd/Gvz4RkVrS4OvglBglQNmU/c1KjOJiJpspwelBNk0T8R87\nqy/kjgSKYP1f4l+fiEgtKTEKTolRAuTSEyODLXxNiE3JDiflJKIb7WeaXQsYbPo3lC5JXL0iIhJX\nSowSII0scugJOJvR2jex9ll0tl9CE6PM3pB3GlAC629KXL0iItXQ4OvglBglSFl32iZmJjmS1BLG\nmZaMFiOApn8G0mDTOCj9KrF1S5WKVyxk7h//j5KFnyQ7FJGE015pwSkxShCtgB0fX0U3jm1LBh1j\ntXFsbWXuBnlnsTyzFdO33JLYuqVK315/Ar22TGDhDacnOxSRhNMYo+CUGCWIBmDHR9n4ooHkxXbj\n2Fpa3+wSHm97Mq82ymNh6dsJr19+au4LL/DBK5H/Yz06h5McjYg0REqMEiSX7hjZFLGUUtYlO5yU\n8XmyutGimmb0ZVBRGm5pPJ/+Dt+zPClxCGxYvpxXfvUrlq6A0rQcMtZ+BT8sTHZYIgmlFqPglBgl\niJFJLrsD6k6Lpc/ivOJ1bQzOvY3e9KfYSvgvD7OR9UmLZUfl4TAvnn02m9esoeuRQ0jvf1zkhZmv\nJDcwkQQXXfH9AAAgAElEQVRTYhScEqMEKSkpYep7KwG4+saR/OM2jUkJalm4lMUUk4uxOzlJi8Mw\nhnMyHdmF9azlGcZRQnHS4tkRTb7jDr56+23yWrfmuHHjsD2HR16Y9XJyAxNJMM1KC05D0BPk2uvG\nUNhkPnsc0pRTfteIs/b+Cx077sLJI09OdmgN1iObCqEJZBbmkZmX+PFFFWWQyS84h4e5i+V8w79X\n3sUPD2xi5/YdOPPMM8nOzk5qfKmopKSEmTNnsmHBAj646ioAhj/8ME3at4dmx0BaOiz6AAp/hLwW\nSY5WRBoKtRglyBtvvsJ+QyKzpho1SefUS3J54039NRvEpC2lhEvT6FKSm+xQAGhEE05hFBQZq9t+\nx5KD5/CPZ27hiGMOp6SkJNnhpZTVq1ezV9++nDB4MM+feiqh4mL6n3cePY89NnJCXgvY9WAIh2DO\nhOQGK5JAmq4fnBKjBGnZsjXzphWVf//9N2FatWqXxIgaNneYvqQlX3+yO5elt052OOVahdvy1C9e\nwt3pmt+JU18/luWF3zJhgm7OsfSHiy+m8aJFnLtpEy3d2WDGp82b//SkPtHuNI0zkh2IxhgFp8Qo\nQf564+3ce/lGli6ItBysXpbFHy65PMlRNVzTN8HyYtg5y9i/Uf35z1xcXMzCN77G3QFISzN6DO3C\n2rVrkxxZapkzaxb9SkrKF2hY4c6s2dtMatgz2no0dwKE1GInOwYlRsEFSozM7CQzm2VmITMbGKug\nUtHee+/NlMnTsPWRj+kf9/2enXbaKclRNVyv/hj5OrQFWHKHF/1ETk4O+x00iDcv+YDiTZGb8V6/\n7UO/Q/dMcmSppf/uu1M2aigMzMjNpd9ee/30pDbdYadesHkdLPwg0SGKSAMVtMVoJnAC8F4MYkl5\nu+66KwftdVbkmzxtIRHE+GhiNKwejql94b8v0GxpW+7p9jjfFnxHboscPuz0BoXaQDgm3J0D1qwh\nG/g2LY17c3PZad99ufLqq39+cp9oq9EsdafJjkGz0oILlBi5+3x3XwBJWHK4gdq6Z5pWwN5e3xfD\nJxsh2+DwZsmO5udat27Nqy+8yprv13BX/kO0Y2d+ZBXP8iillCY7vAZv6r33suTtt8lp0YJfv/02\nH02bxpvvvktOTiVLNlScth/t3hRJZRp8HZzGGCVYNp1Jpwml/EAJK5MdToM0YS04cFgzqEfDiyqV\nRTYn80ua0JSlfMVrPIujG/T2Wjl7Nm9edhkAx/7rX+x16KH07NkTq6o/tcsgaNwaVn0F381NYKQi\nyaExRsHVmB6a2VtAxelTRuS+NMbd69Q+PXbs2PLn+fn55Ofn1+XylGCkkcsebGQyhcyiGYclO6QG\nZ/yayNf62I1WmaY0ZyS/5DHuZSozWE0LTuMIsvVLqE5Kt2zh+dNOI1RUxIBRo+h94ok1X5SWDr2H\nwiePwsyXoX3v+AcqO6yCggIKCgqSGoOSm+DMY9C8bGbvAn9w98+rOcdjUVcqWMZtrOTf7MSFtOei\nZIfToBSHofVU2BCCxQNhl+QteF1n85nFVcxgE435Pe05hn2SHVKD8sallzL59ttp2b07F0ybRlbj\nxrW7cPpz8NBJ0PUAuPSj+AYpUoGZ4e4JG2piZp7vwZYGKbBjEhpzfRTLDsUd+oOsi63jjGYmOZKG\n54P1kaSoT17DSooAetKHXfmGGZRQSKNkh9OgLLr990y+/U7SMtIZ8dRTtU+KAHY/CjKyYPHHsGEl\nNGkbv0BFkkwDqIMLOl3/eDP7BhgEjDczrWJXC3n0AWAzszXepI7q82y02tiH7gB8RyjJkTQc7s67\ndz8KQP75I+iwTx1b2nKasLH9vuDOw78/ljfffDMOUYrUDxp8HVzQWWkvunsnd8919/bufkysAktl\nWexMBi0o5UdKWJ7scBoMd3ilgSdG3YhMo/ua9UmOpOGwcIgzji0lfxAceP3f63z9vHnzuOG5qQDk\nl3zCuaccz3PPPRfrMEXqBQ2+Dk6z0pLAMHKj3WmFzK7hbCnz5RZYtAVaZsCgJsmOZvt0oykQSYzU\nWlhL335Ojm3kkGE9SGvVuc6XP/jPe2iXFdmOp1szeOigzdx207WxjlJEUoQSoyQp607Teka1VzYb\n7f9aQHoDHdHWgmyakcVGSviBLckOp2FYWBD52j1/uy4vKS4mFN76fcfGka1bRFKRWoyCU2KUJI2i\niVFhA02MljCZiUxNaJ3jK2wD0lAZRtdoq9FXrEtyNA1EwMTo1DPP4e+z85i9OvL9rLUZnPmrC2MS\nmkh9E4/EyMw2VHF8nJmNiEXcZlZlM66ZTTCzH83s5WrOyTKz/5jZAjP72Mw6V3jtqujxuWZ2VE2x\nKDFKkoY8AHsdW7iCBdzFImbwYULqXFsKH26AdGBI8xpPr9cqdqdJDUIl8FV0n7PtTIwOOOAAHn/m\nBe5bFVnDaGjPxlx08e9iFKBI/RKnLUHidpMys5vMbDjQyszuMLO+lZx2C3BGDUWNAta4ew/gjug1\nmFlvYCTQCzgGuM+qXBE2QolRkmTSlgzaEGIDxSxNdjh1NIPdmU2YdP7FJkriOMPK3bn51r/T+bTR\nlDp0WLOY5ukNK5HcVlclRrX3zedQtBHa7AbNdt7uYo466ijufX0mNG1Pk9K12LIvYhikyI7DzO6J\ntry8CbStcHygmRWY2dRoC0+76PGLzWy2mU03s6e2Lc/dxwBHA6cD97r7jErOeRfYWENoxwGPRp8/\nC+WrJw8H/uPupe6+GFgA7FtdQUqMkqihjjP6nmc4jIm0IcQS1vNUHAeQP/LoY9zwwDg2jLoRgBXv\nvMhtd94Vt/oSQTPT6mDhu5GvPQ4NXlZaGuwZ3VR2ZpUt8iINWjyn60e7zXq4ey/gbOCA6PEM4G7g\nRHffBxgH/CV62RVAf3fvD/y6kjJvACYATwKjzWzP7XzrHYBvANw9BKwzs5YVj0ctix6rkhKjJGqI\n44yKWckaCsjG+R17kQa8wHzmsiou9f3n5fEUnnkNNIkMLCrpPoBnXn41LnUlSmcak4bxLRsp1npG\n1Qs4vuhn+kQTo1l12s1IpMGI8+Drg4GnAdx9BTAxerwn0Ad4y8ymAWOAsibeL4CnzOx0+PkvPHe/\nJrq92Cp3/727x2rl4+2eoqPVnJKobMr+5gY0Zf97ngdCtOBIdqMHJ7CF55jH7XzCnRxFbox/pFo3\nb4aVbNnawb1oJi2bN4tpHYmWRTodacxSNrCEDfSggQ+aipdQCXwVHcO26yGxKXO3wyEzF775DNYu\ng+bV/uEo0uDUdWbZ5oJP2FwQeCKNAbPc/cBKXhsKDCbSpTXGzPq4e3jbk9z9+oAxfAt0ApabWTrQ\n1N3XmNmy6PEyHYm0GlVJLUZJlFe+ltEcvAG0HIQpYSWRhfHaMRKA0+hNF5rxHZt4hJ91DQc29qor\nyCyMzN6yae/R+LHr+dvYa2JeT6J1K5+Zpu60Ki39FIo3QdvdoVn72JSZlQu7Hxl5Pmt8bMoUacBy\n8/el5djR5Y8qlLW+vA+cbGZpZtYeKOvjng+0MbNBEOlaiw56Bujs7u8BVwJNgTrs5/OzGKprBXqF\nSPcewC/Y2pr1MnBKdNZaV6A78El1FSkxSqJMWpFJe8IUsoWvkx1OjX6kgBJWkUtXmrIXAJmkcwn7\nku7GBBax90UjOHzY0cycGZvW0B49erDrqZGp1eekL2f6lI/p27eySQsNiwZg10Ksu9HK7Dk88lXd\naZKC4jkrzd1fABYCs4FHgEnR4yXAScDfzGw6MA3YPzr26Akz+wL4DLjT3ev8S8/M3gf+CxxmZkvN\n7Mjo8evMbFj0tIeA1ma2APg9kUQMd58DPAPMAV4DflPTjvbqSkuyRuzJWlZQyGxyo/to1Vff8wwQ\naS2yCol7V5qT/vxcQifuzk5/P4llT3zA4CMOZeZn0+nYsWOgOpcWwdySLBqnwT/PP5XsFEnllRjV\nQiwHXle0x9DI1/lvQ9EmyNaGvpI64rHfmbs3rfD8oirOmQFU1ud9cAzqH1zF8WsrPC+CaFfGz8/7\nK/DX2taXIreZhmvrOKP6PQC7kEWs51PSyKU1Q3/yWjgc5pWz/0q4NERadiadRh1GqyF9mTAh+J7C\nr0YXdTyqOSmTFMFPu9Ia2jpWCVFaDF9/FHnePUbji8o03Qm67AelRZHkSCSFaOXr4FLoVtMw5TWQ\nPdO+538AtGYoGfx0ozIzI83SKF61tfUjr3dHsrKyAtdbtg1IQ900tiqtyaExmaynmDUUJTuc+ueb\nT6G4ENr1gibtYl9+H03bl9SkxCg4JUZJtjUxmotTkuRoKheikFVEBqq24xc/e93M+OMfL2P6UX9j\nzUfzAOh4wWEcesIxgeotDMHEaK71fymWGBmmAdjVWRCnbrQyZeOMZo+H8M8myIjIDkyJUZJl0Ixs\nOrMZ2MyiZIfzE+7O+PHjefy1iwixicben0bsVum5111zLbf+8Vo6PvQluQvWkdGiEY81XRCom2ji\nOtgShn0aQ7vgjU/1ztZxRtoz7WfiNfC6TPs+0HIX2LASllQ7QUWkQQmF0wM9RIlRvbCMXnzEgczn\nvWSH8hNXXn0Zv7v8LFoPmALA8/d/T1WD+c2Ms888i/88/Dj39DiNRmTyCcuZyJLtrr9s09hU60Yr\nowHYVSgt2jq+KFbrF23LTLPTJCWVlqYHeogSo3ohnXY4abzPXNbzY7LDAWDlypX885/38cDkDrRu\nn4m789BtHzFjRs1rFbUhj/PoD8C/mMYPFNa5fnclRjuspVOhZDPstAc0aVvz+dtL44wkBYVKMwI9\nRIlRvXAw59GFXSkmzHgeJVwPFntcu3YtzVrm0Lhp5D+KmdGqXS5r166t1fWHsguD6EAhpdzJVMJ1\n7FL7ohCWFcPOWTAgRWdTd6EJBizV1iDltmzZwooPno58E69utDLdD4GcJrBiFqyq/+uIiUhiKDGq\nBzJpxVB+SWOasYyv+Ijg09yD6tq1K23bNSMUctydic+uZcXXJfTr169W1xvGaPaiGdl8Fl7D8163\n8VNls9GGtoj0eqSiHDLoQCNCON/UuHF06lu4cCF79OjKktcfAOBf78ypsus2JjKy8H0OJDwYwsuu\ni189IgkUKk0P9BAlRvVGHk0YylkYxmTeYjHzkxpPZmYmj7z6W9LTjckTNvG/m3N5/bWJNG9e+329\nmpHNaPZiwabunLWqJX/etIFNtbzRlXWjDU3RbrQy6k7batQZJ/P77t8xqF2k9ezpN6bw7LPPxrVO\n77sr4SHgWdoeRFKDEqPglBjVI53pwf4cDTiv8hirN69k48bktCQ4IUpavgPAuf/3KNM+nVPr1qKK\n9vOdsZJ2FIWzuWljGt1XFXLrhiIuve56OvTcne4DBvDcc8/95JqVxfDJRsg2OLxh7xdbIyVGW82Z\n9yWnVZj0eECbzcyeFd+FT63jeQB4+y1Usq+lSINTWpIe6CFKjOqd/RlCJ+9OIRu4+qMLaNOuFceO\nGEZhYd0HMAfxIx9QzHdk04lmDNrucsyMaS2a81DzEPtkpLEyDFcWlnLnKb9kxTPj+epPN3HWby+i\noKCg/JoJayMb8xzaDBqn+P9TJUZb7b7briyMrlxQEoLXl+fRq3fv6i8KqlFfsJ0hbROEk9tKKyL1\ngxKjeiaNNH54NMzmH4vofkQnrl13HssyvuLKP12R0DjKVrpux0lYwB+TxpbFOdlN+bhlDq80zyZz\n7ixo3gLbqT120CFsvu9hnhy/tSsj1WejVdRNiVG5fz/+X0JpmQD8+v0sBhx+PCNHVrr1UcyYGZZ+\nAAAe+jiudYkkQjiUEeghSozqpQ/fmcQ3n34PQHpGGoeM6cdHkz9MWP1bWMo6JmFk04bhMSvXzPi/\n7Ay6Xfk7/LOpW48fcDBPj76cWzcVs6bUeSM68S3VxxcBtCOPPDJYQxFrd/CtQXp224X9d478Srrq\nsfd4cNzjWAJG3lv6/gB4aFLc6xKJu9L0YA9RYlQf7dJxF2Y/83X5jJyWXZrRrWeXhNX/PZEBr60Z\nQia1H2xdW7decw0555yM3/43bNyDpM+Yxua8Rly5sYRdv93ChhD0znW65MS86nrHsPLutEW+g6+A\nvWQKVloEO/ele79BCUmKoGJipBYjSQFKjAJTYlQPXX3FGNZ+Usy4w17jh3lryWmWxfH35CdkF/Yw\nW1jJSwC0Iz7dGEOHDuWdl17k96HNXL15LYt6dWFC82z2z0xj7aZIU25Wox1jXZ8NGzYw5/XIKs+n\njPktd9x9V5IjSqKy/dHivX7RttIHAlkQnoN77dbpEpHUpQ7Feqh58+Z8Nvlz3n77bQoXr+e73b5g\neaNFTON9BhKnLRKiVvEmIdbTiD1oHN3gNh72339/9t9///LvOwFHZqUztjjMvaEQV7XdMXL28y76\nDSsOacVO9KLzlSO4rv8f6Nm9B8ccE2wD3gapfH+0OG0cWwWzbEjfC0If46EpWMaQhNYvElOlKbrw\nWwLtGHefBig3N5djjz2Wk48+naPTTgOggBf5nm/iWu/3PAPEr7WoOmbGdTuls6pvOiOb7xg/mgUF\nBTQ/si8AGU3zaD3qUN4pmJjkqJKgeDMs/jiymueugxNefXl3WqnGGUkDVxrwIUqMGoKe9Kc/BxEi\nxMs8QjFb4lLPRmazidmk05TWHBWXOuSn2rRtw7pPFwLg7hTN/Iad27VPclRJsGQyhIqhfV9o1DLh\n1ZclRmickTR0SowCU2LUQBzKCbShA2v5gdd5Mi7jjRbzIABtOY40doCRz/XAP2+7i6/Ou5/ilesw\nM1qEsjj//POTHVbilXWj9UhsN1qZrVP2p+C+Y4xvkxSlxCgwJUYNRAaZDOccskhjBe/zNv+OaXI0\nn/Gs5wMAGpEfs3KlegcddBDTp3xKl8JsAG558gEaN26c5KiSIFkDr6MsbWewzsB6CM9JSgwiUj8o\nMWpAWtKO3vQiDWcxH/ABDxOiJFCZTpjPeZqpjGMdTYEetKR/bAKWWunWrRuHdomMM1qWFZ9u0m19\nt3EaocIVsS94yo0w5f66XVNcCEumJG18URlN25eUUBLwIUqMGpojuYCeDCGdTObzHuP5K4Wsxd0p\nKandT/X8+fN5++23+XbFEt7jDmbyImGy6MYY9uO/pOnHIuESuTXIii/vpuXD+7Hkg9MJEcNuo2Vf\nwM1/husvhDf/UvvrFn8cGV/UoT/kJW9VT62ALSkhFPAhugM2RIM5nWP5E41oyUoW8sSmP7DnEZ3J\na5RH/pGDWbVqVZXXXnfjWA44ZF+uue9CHl5+PkuYQia5HM6V9ORIDE31TIayrUG+inNiNJ8veabN\nQtLDIbrOLsB+mBm7wjOfgkEe+eV62xh47h+1u65sfNGu+bGLZTtYhlqMJAVojFFgSowaqDZ05QSu\nJ3fdTlijEi5+K597Np9MZq81nDXqjEqv+eKLL7jnn3dw/exD+fXzA9llr+asXlLIkSXX0oF+CX4H\nUlF7GpFNOqvYwnqK41LHHObyJE+zskUzFvcbiuGkvR+jPfhKlsKGO2EEcM5FkWP/ugwe/AOEw4RC\n1fwpmuSB1+XS+gE5EJ6Ph1cnNxYRSRolRg1YLk355v7mLProByzNSEtP4/i/9WXSpMrXYvnqq684\n5LxuNG6VXX7sjmM+pnSVZqAlW1qFrUFi1Z0WCoV4//33mTBhApM3TuFp/kuIEAcwiG6DHobsZrDk\nTfj69eCVrbkGvAianAyn3AWXPwHpGfD8bbyR34LGWRl0aNOS11/fpq6K44u6HRw8jgDMsiB9bwA8\nNDmpsYhsN7UYBabEqIFr12YnXhs7l3A4MkMtKzeDi149lM2V3FxbHbCRodd0L/9++vhlFK8z2rZt\nm7B4pWqxTIyKioo47OjDOG306dw27XZeynmFMGEO5iD+j2Ow3Naw358iJ79/GYQD/EYsmg4bHgcy\noVV0bNFhp+PXv8qmsDEkbz2FR8AzXX/kzJEnsnjx4q3Xfj0JQiXQYQDkxX5fvrrSOCNp8JQYBabE\nqIE744wzaFa6M3cc/B4Trp/L+u+20HnfZrzItaxiMQBhwnzCI3zdbjxp6caEv33J1bu+w1Pnz+PF\n514hPV0bB9YHsUyM/nn/P/kh+wfO/OJs9r16EGkZaSx9YjFDKo4j6/9baNoFVs+G2Y9sf2WrLgcc\nmo2GzG7lh9d02YtjpkV+tszgwJYwvG0aU6dO3Xrtwug0/WR3o0VpoUdp8JQYBRZorzQzuwU4FigC\nFgHnunv8p9VIuaysLCa++R6vvPIKP/74I/sVDmARr7CSRbzMjRzAGSxjCt/yOWlkcAAXcMKF/fl+\nxPd06tSJnBx1o9UXsRyA/dXir+h30QDS0iJ/+2z6YRPvjS3AzqgwuD4jBw6+GV49BSZdAz1Pgaw6\nrqFU+CZsfgvSmkHLP/3kpaZNm/L5xjS+LoSueZFj13faxPKMoq0nle+Pll+3euNk65T9KbiXYqbt\nJEV2NEFbjN4E9nD3/sAC4KrgIUldZWZmMmLECEaNGkWfbgMZxtXsxmBCFPEx9/Mtn5NNY47iT+zK\nYJo2bUqPHj2UFNUzXWgCwGLWEwq4eGf3U3ej25BdAfCwM/nGSey9994/P3G3kbDTfrDpO/j01rpV\n4qFoaxHQ4mpIb/WTlzMzM7njzrs48Is8fr8ohxmb0uiQ7ez91O9gzsdQtAmWfAKWlvTxRWUsrR1Y\nN2AThGclOxyRulOLUWCBEiN3f9vdw9FvJwMdg4ckQaWTyWBG0YcjMBwjkyFcSzt6JTs0qUYTsmhL\nLsWEWc6m7S7nIybx/b4rAXj38onc2+4uQp+W8uA9D/78ZDM45O8AFE26iXOGD+Lvt9xc7Syyb775\nhhHHH8XVF7WF4i8IWQdodlGl5/7q/At4ZeL77P6H21n+h2fwfYdiG9bAVYfDa/+IjG3qOBBym233\n+401TduXBk2JUWDmHpttJczsZeA/7v5UFa97rOqS2pvBS3RgAK3onOxQpBb+zBSm8D1j2JvB7Fzn\n69/nA97gLQCOZRi91vWksLCQnXbaCbPK16has2YNU8a055gexSxeA2f+L48Bh5/OXff9PJHavHkz\n/fvtxtnHLueyX4XJyoRr7t6Ja//xDRkZteh2CpXCnRfAmw9DG6AVcOhlcFwdW6viKFx8L+Etv2Xl\n6iPJafYfWrZM/Ka2khrMDHdP2OJwZub8J+B99pTExlwf1dhiZGZvmdmMCo+Z0a/HVjhnDFBSVVIk\nydOX45QUNSDbOwDbHe5dVcAbvIVhHM9xDGJfmjVrRvv27atMigBeffVVXpkTGSTdpSWMP6uQBx8a\nV2mr0fTp02mUvZ6rL4wkRe7w9MvrWbBgQe0CTc+AS/4Np10D0XFHbK4/v4Pdnb/+7Q0ASoreoXfv\nrnzyySdJjkpEEqnGP/Hc/cjqXjezc4D/Aw6rqayxY8eWP8/Pzyc/P7+mS0R2KF23cwD2zcveY2PH\nibgbJ9kJDKjDfnfuzvcbt/6NlJcFVDHGKTc3F6MU90gvXHEJrN8QJjc3t/bBmsGZ18GaT2Dxh3D0\nJbW/Ns5effVV/ve/d7j8Yuiwc5j771jPOef8gjlzliQ7NGkACgoKKCgoSG4Q2tYjsEBdaWZ2NPAP\nYLC7V7tUrLrSRGq2lA2cx7u0I5fHqPZvknJTl8PIiYsZcdKT9Fx9LOd36FunOlevXs2APXdn9F5r\nuHRwmMx0uPeboxh9+xs/OzccDjN+3C4MP+RbliyDs6/Mo2PXo3n8yWerbZWqUigE9Wi5iDvuuIMF\nc6/kzlsjM+eKiqBFh3SKizX4QuouKV1pjwa8z56trrSgs9LuBhoDb5nZ52Z2XwxiEtlhdaARmaTx\nPZvZVIutrn/cDCNfgMVLu1A86ZI6J0UArVq14oOPP2V+81/w6Je7APCbfssqXfQxrfQrjj3kO9yN\nJ94+jpFn3MKjjz+zfUkR1KukCKBv3768/lY6mwoj38+ZB3337F79RSL1SRwGX5vZhiqOjzOzEbEI\n28w+NbPMbY6NNrMFZhYys0oH+5lZPzObFB3mM93MRlZ47Qkzmxcd/vNvM6vVL5xAi3S4e48g14vI\nT6WTRheasIB1fM16+tCqynPd4dxXYfE62Ls9/H1wXpXn1mSXXXbh4cf/A6Vb4JFe2Jrooo97/uqn\nJ64Zg1EKTc9hzA3jtru++uqwww7j7HMu5bzRN/PUuFJKSjJ54skXkx2WSO3Fp3Ezrt09ZtYF+Nbd\nt/1r8EPgFaCgmss3AWe6+yIzaw98ZmavR9dUfMLdz4jW8RTwK+CBmuLRytci9UxtB2Df/gm89CU0\ny4ZnToDsWKxFWLboI0QWfSzeuPW1LVNg4zNgOdDy+hhUVj/9+c83cN/9MwHYZ+8Mevb8//buPD6q\n+tzj+OeZ7EjYBEFWZREXbFFUJIBSEcRLxQ3trdar3i7aWu1ildZ61bbWq9jrUu2iXVxqrVVcq3Wr\nGlEoCIqIIAgIGHYFNARZsjz3j3MCY8gyyZlkJpPv+/WaF5MzZ37nOR4n8+T3e36/MyDFEYmkDzO7\n08zeM7MXgP3ith9pZsVmNsfMnjWz7uH2y8xsYdibU9cErQnAXjdtdPf57v4hUGeXtLsvc/fl4fN1\nwEaCOa+4e3ybb5DgkkJKjETSTH0F2CUlJdx222384Ja/M+Xl4I+4e78MBybzNmO1LfroDh9fETzv\n+H3I6ZPEA6afrt0OhtggjO1Q9U6qwxFJXDOuYxQOmw1y90OA84GicHs2QWnNme5+NHAPEN44kSnA\n0HAh6IvraLrWxKixzOwYIKc6UYrbng2cl+gxlBiJpJn+dfQYLV68mCFHHM2UR1dy2+YTqXDjwoM+\n4bTBSQ4gbtFH5t4MW9fAtqdgx2sQ2xc6/zjJB0xPe24PMjPFkYg0QvMu8Dga+Bvs7p15Odw+GBhC\nUG88D/gp7F6IbT7woJmdSy1z5sK6ol7uvrKxp1qjnf2B+4ELann5t8Cr7j4jkbZ0IyCRNHNAXGJU\nheK/9TgAABwySURBVBMLe5F/fM0vKD36chgd9tyUrqH0iWth8h+TH0SvUTDwDFj2GMy8Gg4OV4Hu\ncg1kpc8q1c3Jsorw8vvxin9Dbu0re4ukncbWGC0rhuXFUY9qwLvuPrKW1yYCxwGTgJ+a2ZC4O2ZA\nkGy93kD79dY4mVkh8DTwE3efU+O1a4Cu7v6tBo6xm3qMRNJMJ/LoQh47qGQ9n+3e/tGmLTB44p4d\n181n08cfNV8go2+CWDaU3QvlSyBnAHSsqyc88+zpMdKtQSSDDRwDJ12351G76hqf6cBXzCwW9tB8\nKdy+BOhmZsdCMHRlZoeGr/V191eBHwMdCGayx5sAPNtAlEYddUZhj9MTwH3u/niN174BnAR8tYH2\nP0eJkUga6k/QKxM/nHbEyRdC14OCH0rX0e7VqznzlAnNF0TngTD0m3Bw+HOXG8Bym+946SZ2GFAI\nvhKvWpfqaEQSUx7xUTsHCBOPZcBC4F5gZri9HJgM3GRmbwPzgBFhbc8DZjYfeBO4PZwtFm8M8Gpt\nBzWzS82sBOgFzDezu8Ptw6qfA2cDo4ALzGxeuHRQ9bolvyMoEJ8Vbr+6zjOMP25LLbqoBR5FEvdH\nFvEIyziPwXyNwXy0DYb+yVlbZuS+cQd503/BpZd8h+t/fm3T1xBKxKfPw5qTYUceHFwC7bo237HS\nUOW2E/HKl4gVPEosJynLtUgbkpIFHn8Z8Xv2py0Xs5n1Au5294kN7tyCVGMkkob6756Z9ilVDl97\nCtaWGSN7wytTLiUnq4VqXjqeBDunQ8eukNe2kiIAsoqg8qVgOE2JkbQGrWiRdndfQ1CDlFaUGImk\nofgp+zfMgBdWQNcCeOg0yGnpxaL3G9XCB0wfljUCR3VGIm2JEiORNNSb9mRjLFxZwF9ecwzjgVOh\nd4dUR9a2WPaxwZPKubjvwtpSjZW0Tq2oxyhdqfhaJA3lEKNrWRfeeXIYVW5cVQQn9U91VG2PWWeI\nHQLshMp5qQ5HpGHNu45Rm6DESCRNzXp6CLu25XN43x1cd1yqo2m72tK0/aodW9i+veGbF0saa55Z\naW2KEiORNPXdUTvos/82fnXaNrL1SU2ZtpIY/f2eZ+nb+3pu+uE1qQ5FJKX061YkTV3auzurLtiH\n8e33TXUobVpbSYw6lT3Jmk0deOSF5JeezuJZ3uctKtQl0fwqIz5EiZFIOmvOJYokQbFDgI7gJXjV\n6lRH0zy2r+aE/e+hc/vtLPogm0WLkrei+jZK+TfP8E/upZxdSWtX6qAao8iUGImI1MMshmUFs9My\nttdo2e3kZO3itDE7AJg2bVHSml7KPBznAA6hgH2S1q7UQYlRZEqMREQaYNkZPJxW/imsvAuAs84b\nDyQ3MVrCmwAMZljS2pR6qPg6MiVGIiINySoCMjQxWnE3VGyFrmMYe9pJdOqUz4IFG1my5OPITZey\nmTUsJ5scBvCFht8gkgaUGImINMCyhgMGlW/iviPV4SRP1S5YfnvwfNAV5OZmceqpg4Hk9Bq9z1sA\n9GcIueRHbk8SoOLryJQYiYg0wKwDxA4DyqHyrVSHkzyrH4Ida6DwUOh+MgCTJx8KwCOPRE+M9gyj\nHRW5LUmQaowiU2IkIpKA3dP2K2amOJIkcYelvwqeD/rR7imQ48b1p0OHPObP38DSpZua3PwWNrKB\nD8klnwM5LBkRSyKUGEWmxEhEJAFmRxFbDax/INWhJMfG56F0AeTvD73P2b05Ly+bSZOiD6dV9xYN\n5ItkkxMtVpEWpMRIRCQBVmbEtoCtnY+Xzk51OE1WVlbGtGnTWP/aj4INA74HWXmf2+ess6INpznO\nYuYCmo3W4jQrLTIlRiIiidj3G/h+IzAH3v8KlG9OdUSNtmnTJo4ZNoSX7/0vemQtpGwHLK86ca/9\nxo8fQGFhLvPmrWf58saf58esZTPryWcf+nJwMkKXRKn4OjIlRiIiCTAzbEAxtD8a27kKll0Y1Om0\nIlNvvJ7jB6zjt1/fDsDi1XDllKv22i8/P5tTTmn6cNqSsLfoII4gi6wIEUujqcYoMiVGIiKJiuXC\n4IchqxNsfgrW3hK5SXdn2up3uWLRK1RWVSUhyLqtW7OKEw7dc1uOyipYt7b225xUD6dNm/Zeo47h\nuBZ1lFZNiZGISGPkHwCD7guer5wCpdFmqdmqF7jysxx+1amIhR+8GD2+eow8bhxVFtwktqIS/vfJ\nAkYdN7bWfU86aQDt2+cyd+5aVqzYkvAx1rOKT9nEPnSkFwOTErc0gnqMIlNiJCLSWPtOgp6XA5Ws\nn348A/u1Y+yY4axatapx7Sx7Ap6cRNHHwYraM3qckPxY43zzgrM4fXjw/NifxCjoPY6f/3JqrfsW\nFORw/PE9ABg79ofceedv8ASGDhd7ddH1kcT0FdPyVHwdmf6vFRFpgs0drmDu4ix6dK5g4X3bGTt4\nDhNPHkNlZYIVrO/9Ff4xGSp3UVQQrAo9s5m/mGKr7iY3q4Kqricy493P+NvDT5KfX/uK1CUlJbz8\n8u8BWLGiO1OmTOWaa66rs+3i4mL6DerDq+ueBqBwbc+kxy8JUPF1ZEqMRESa4M233mHq39sBkJcD\nPznP+WTzRkpKSmrd/+OPP+arkydx0AH7c9v5/fFnzwOvhOFXU3TYWQDM3NWMNUaVO2D5rwGIDb6S\nvLy8enefNm0aFRXVXxEFfPbZ6dxxx29q3XfNmjWcftZpnPbAcDr03Iftn+zknPEXJNTDJJJuslMd\ngIhIa9S5c2feXFzFrnLIDdcvHHLgLjp06LDXvu7Oqf9xIke1W8RrF5bTvTDYvv2oaykYeR1D3Glv\n8EElrK90emRZ8gMueQB2boCOQ6Hb3lP0a4vZLD5Ry6tzEt6cOXM4YHgP+g3vDkB+x1zWr1/Pxo0b\n6d69exKCl4SpTigy9RiJiDTBsGHDOObYsRz/vX145a3gjhoPX59Hl/Z7fzOtW7eOpUsXc+uX9yRF\nv56Zz/TtwW1Gss0YnhMkQ/9ujl4jr6r19h/1mTx5Mvn5y4HFAMRiZVxyybdr3bdbt258vOKT3T1E\npes+Y9f22pNEaWYqvo5MiZGISBOYGX996HG+e+VdzPj0WjaVH0aH/G3w/rnBEFmcgvx8rhpTTiz8\njVtVBX+Yk027du1271OUG7w4c1czDD+tfxrKlkBBH+h1dkJv6du3L7Nnv84ppzixWCVVVftx8cWX\n17pvUVERJ58/FjNj67rP+O3wp7jxxhspKChI5llIIlR8HZm11BiwmbnGm0UkY+1cC/OHQvlH0Odn\n0PeaYLtXwcvfhfm/o6ISnl0Md80pYGeXoTz30mtkZQULID6/o4oJm8spyjFmdMtNbmzTR8Om1+Hw\nW2DgDxr99smTH+bRR9/j1ltP4vvfP7bWfZ71+3jP3iD27305xidQVFQUNepWz8xw92YYF63zeM4R\nEb9n57VszOlIPUYiIsmQ1xMO+itgUHIdfPISVFXA8xfC/N/hWXlM7/wDXmr3bcb99w08/fwru5Mi\ngOG5hgFzy50dyfwjcvOsICnK6Qj9vtGkJvYs9lj7KtgVlLPc3gHg/BGXKCmSVk3F1yIiydJpHPT5\nHyj5OSz+Kmw4BpY/Azn7YKc+yQl9x1LXSkWdYsZh2ca7Fc5b5U5RbpL+aF96c/Dvgd+GnMImNTFx\n4kHk52czY0YJa9aU0qvX52uHVrCQXexgP/rQGRVbp5TqhCJTj5GISDL1uQbaHw8LPwqSoryOcOYL\n0Lf2FabjVSdDyZq2v73sn/jax3HLgQGXNbmd9u1zOfnkYBXrRx/d+xYh1bcAOZijmnwMSRIVX0em\nxEhEJJksC3pOhc9ikG3w5T9Cz8SGlpJdgF217BcYzvY+X4D8/SO1Vddw2i528AELADiIIyMdQ5JA\nxdeRaShNRCTZuh0Dk+6HnO7Qs+E1g6rF9xgF6whFG077qF+M8l3tyRlU+2yyxvjylw8iLy+L11//\nkHXrtrL//sGw3HIWUEE5PelPB7pEPo5IqkXqMTKzn5vZfDObZ2bPmVmPZAUmItKq9Tu3UUkRwMAs\no1sMNlTBioi3Z9jFako7f8zqYwZR0OHMaI0BhYV5TJgwEHd47LE9w2kaRkszuiVIZFGH0qa6+xfd\n/QjgGeDaJMQkItImmdnu4bQZEeuMSnkVgEJGESM50/8nTw6G0x55JBhO2842VrIIwxjEEUk5hkSk\nGqPIIiVG7l4W9+M+QDPe6EdEJPMV5SSnAHvr7sTo+MgxVTvllIPIzc1i+vRVbNhQxjLmU0UlfRjM\nPmiV67SgxCiyyMXXZna9mX0InANcEz0kEZG2KxkF2FVsp4xZABRyXFLiAujYMZ/x4wfsHk6rHkYb\nzLCkHUPSj5ltrWP7PWZ2RpKOMdfMcmpsu8TMlppZpZnVW8BmZoVmVmJmv47b9pWw3GeBmf1vorE0\nmBiZ2Ytm9k7cY0H47ykA7n61u/cF/gpcmuiBRURkb8NyjRxgQYVTWtW05KiM2Tg7KWAIOXRLanzV\ns9Oe+tfblLCEGFkMYmhSjyERNM+stGa9bYWZHQCsdveaEbwOjAVWJdDMLyDsJg3a7AJMBb7k7ocD\nPczsS4nE0+CsNHcfl0hDwIPAP4Hr6trhuuv2vDRmzBjGjBmTYNMiIm1DgRlH5hizy53Zu5xx+Y2f\nmbaVYgAKGZPc4IBJkwaTkxOjrOdyHOdADiWfdg2/sQ0oLi6muLg4tUE0cwG1md1JkKyUEJdKmdmR\nwC0EZTUfAxe4+wYzuwy4KNx3kbufU0uzE4Dnam509/lh2/V+CMxsGLBf2Eb1LID+wPvuvjn8+SXg\nTOCVhs4x0nR9Mxvo7svCH08D9l75K058YiQiIrUbmRtjdnklM8urGJffuIoHx3cXXndIYn1RtU6d\n8hk3bgD9v/IGoGG0eDX/4P/Zz37W8kE0Y51QOGw2yN0PMbP9gUXAn8wsG7gDmOTum8zsbOAG4OvA\nFOAAdy83s7oK0SYA329iTAb8CjgXiO/IWQYMNrO+wFqCHCVn7xb2FnUdoxvN7CCCoutVwMUR2xMR\nafOKco1btjVtZtpOllLOWrLpSgFDmiE6OOOCvnw6agaVO2IMyD+8WY4haWk08DcAd19nZi+H2wcD\nQ4AXw0QlRpCMAMwHHjSzJ4AnajYY1hX1cveVTYzpO8Az7r427FiyML5PzOzbwMME/WgzgQGJNBgp\nMXL3yVHeLyIiexsRFmDP2uVUupPViIUeS3cPox2HNdPNDYpO6sJflmZRua4TWSNzIavh90gLaWyP\nUWUxVBVHPaoB77r7yFpemwgcB0wCfmpmQ9w9PuMfTVBLVJ/6apxGAKPM7DtAIZBjZlvd/Sp3f4Zg\nKSHM7JskONBonsy7ONd3IDNvqWOJiLR23VdsYWNeOy586PfcdN45dOuWWBH1cs5lG3Ppy+10YkKz\nxbe1bCd57Z1c8pvtGK2dmeHuSbobcELHc7Ijfs9W7B1zmGgUmtnpwLcIkp3uwELgG8A/wuf/5e6z\nwqG1g9x9kZn1c/dVYc/QCuBQdy+Na3sq8C93f6Ge81oBHOXum+oL3czOB4a5+2Xhz93c/SMz6wy8\nDJwVV/5TJ90rTUQkzdx6++1sejUYpbgvvzNfPHYEW7ZsafB9FXzKNuYB2RRS2x/vyVPYPk9JUTpq\nnpWvHcDdHyeo3VkI3EswPEU4m2wycJOZvQ3MA0aECdIDZjYfeBO4PT4pCo0hbjZZPDO71MxKgF7A\nfDO7O9w+rPp5A243s4XAa8ANiSRFoB4jEZG007F7D0qfnw37BXdZKrj4XG6dOI6LLrqo3vd9wjN8\nyA/Zh+EM4P6WCFXqkZIeo8gz61suZjPrBdzt7hNb4niJ0k1kRUTSTPnOHZCz5zYelZ27sH379gbf\nV11f1KEZpumLJJu7ryEYlksrGkoTEUkzk88+m/zvfx0qg7GN7H+/xsSJ9X9/OJVsZTqQ3NuAiLQ1\nSoxERNLMH+64g28OPYz8d4JbbvzPI48xaNCget+zoXQGlXxCLn3Io39LhCmSkZQYiYikmby8PH59\n881cMXoEAJsOqDsp2rp1KydPPIGpd50NwJvFVXgTbyUiIkqMRETS1p4byta90OMVU75HdpcFnP+j\nYFHh555YyF13/75F4hPJREqMRETS1LG5weSgueXOzjpm9c5+Ywbn/qCA6ttJDTgcZr/xWovFKOmm\nee4i25YoMRIRSVOdYsZh2cYu4K3y2hOjfn0PZMvG4DV3583iKvr1HdiCUUp6qYj4ECVGIiJpbGTY\na1TXcNr/3fwbYhbcG/Pem7axbkkPfnT5lS0Wn6Qb9RhFpcRIRCSNVdcZzdhVe4/RgQP6UDQ+qC86\n+dhfMeO1uRQWFrZYfCKZRgs8ioiksaK4HiN3311LVG0bc3DbTj6DGT/mq6kIUdKKhsOiUmIkIpLG\nBmYZXWOwoQpWVEL/Gr+1t4arXRdqtWsBNBwWnYbSRETSmJnVO22/NLz/Zgetdi2AaoyiU2IkIpLm\ninJqL8DeyQp2sYosOtGOoakITSTjaChNRCTNjcyNAZXMrFGAXd1bVMgojKwURCbpRzVGUSkxEhFJ\nc8NyjRxgQYVTWuV0iAU9SKovkr1pOCwqDaWJiKS5AjOOzDGqgNlhr1ElZWxjLhCjkNEpjU/SiRZ4\njEqJkYhIK7C7ALs8qDMqYyZOOe0YSjadUhmapBUVX0elxEhEpBUoqrECdmk4jNZBw2giSaUaIxGR\nVqC6x2jWLqfCK9lq0wEo1DR9+RwNh0WlxEhEpBXomWX0y4JVlbCoYhHkfEQOPchncKpDk7Si4bCo\nNJQmItJKjAx7jUq8GAh6iwyr5x3S9qj4OiolRiIirUT1cFp2rHq16zEpjEYkM2koTUSklSjKNbrE\nNtE9612MXNpzbKpDkrSjobSolBiJiLQSh2cbJ+S/TsycHB9OzNqlOiRJOxoOi0qJkYhIK5FtxikF\nrwGwtnw0B+emOCBJQ+oxikqJkYhIK9KOIykp30BHH5PqUEQykrl7w3sl40Bm3lLHEhERSTUzw91b\nbNqgmTk8FbGVSS0aczpSj5GIiEjG0FBaVEqMREREMoYSo6i0jpGIiIhISD1GIiIiGUPT9aNSYiQi\nIpIxNJQWlRIjERGRjKEeo6iUGImIiGQM9RhFpeJrERERkZASIxERkYxREfGxNzPbWsf2e8zsjGRE\nbWZzzSynxrZLzGypmVWaWZd63nuTmb1rZgvN7LZwW4GZPW1m75nZAjO7IdFYkpIYmdnlZlZVX+Ai\nIiLS3MojPmrVrLetMLMDgNXuXjOA14GxwKp63jsCKHL3IcAQ4BgzOy58+WZ3PwQ4AhhlZiclEk/k\nxMjMegPjqCfwtqC4uDjVITSrTD6/TD430Pm1djo/aZzk9xjFM7M7w16YF4D94rYfaWbFZjbHzJ41\ns+7h9svCnpy3zezBOpqdADxXc6O7z3f3D4H6blHiQL6Z5QMFBLXTG9x9u7u/GrZTAbwF9G7wBElO\nj9GtwBVJaKdVy/QPdyafXyafG+j8Wjudn6SLcNhsUNgLcz5QFG7PBu4AznT3o4F7gOqhqynAUHcf\nClxcR9O1JkaJcPdZQDGwDlgDPO/uS2rE3Qk4BXgpkTYjzUozs0lAibsvMGvT95wTERFJA42dlfYB\nsCLRnUcDfwNw93Vm9nK4fTDBMNaLFiQDMWBt+Np84EEzewJ4omaDYV1RL3df2cjAq98/ADgY6EnQ\ns/QvM3vO3WeEr2cBDwK3JXqMBhMjM3sR6B6/iaDr6mrgKoJhtPjXREREJCUau45R3/BR7ZWmHNSA\nd919ZC2vTQSOAyYBPzWzIe5eFff6aIJaovrUV+N0OjDL3bcDmNmzwAhgRvj63cASd7+j4dOoPpp7\nkx4E2eF69qSb5cBKYL869nc99NBDDz30aEuPpn7HNvF7eWUSYl5ZS7tbw39PB54l6BHaH9gMnAHk\nAO8Dx4b7ZQOHhs/7hf/mAKuBDjXangqMb+C8VgD71vHa2cALQFZ4jH8BE8PXrgceaex/RwvfHJmZ\nrQCOdPctSWlQREREUs7MSt29Q/j8DuBE4EOCDpE/u/tjZvYFgjqjjgRJym3AfQRdUB0IepX+4u43\n12j7DWC0u++s5biXAlcSjFptBP7p7t8ys2HAReHzGPBbgl6pKuBZd7/CzHoBJcB7wC6CpO9Od/9z\ng+ebxMToA+Aod9+clAZFREQkY4XJy93uPjHVscRLWmIkIiIi0tq16MrXZnatma02s7fCx4SWPH5z\nMLMJZrbYzN43sympjifZzGylmc03s3lhl2erZmZ/MrMNZvZO3LbOZvaCmS0xs+fNrGMqY4yijvPL\niM+dmfU2s5fDNVEWmNll4faMuH61nN+l4fZMuX55ZjY7/F2ywMyuDbcfYGazwt+hfwunfrc69Zzf\nPWb2Qbj9rXDISdJYi/YYhf+jbHX3W1rsoM0oHNt8n2BlzrXAHOA/3X1xSgNLonCIdFim1I6Z2Sig\nDLjf3b8QbrsJ2OTuU8PktrO7/ziVcTZVHeeXEZ87M+sB9HD3t82sPfAmcCpwIRlw/eo5v6+QAdcP\nwMzauftn4RTqGcD3gB8C09z9ETP7HfC2u9+V0kCbqI7zuxj4h7s/ltroJFGpuFdaJk3pPwZY6u6r\nPFjK/CGCX2SZpHpNiozg7q8DNZO8UwmKBAn/Pa1Fg0qiOs4PMuBz5+7r3f3t8HkZQVFlbzLk+tVx\nfr3Cl1v99QNw98/Cp3kEM5cc+BLwaLj9PoKZT61SLedXPS09I65fW5GKL7xLwqXB/9hau7zjVFe9\nV1vNnl9kmcKB58Nl3r+Z6mCayX7uvgGCLyfilrnPIJn0uau+t9JQYBbQPdOuX9z5zQ43ZcT1M7OY\nmc0jWOrlRWA58InvWddmNcFCfa1SzfNz9znhS9eH1+//rMaNUiX9JD0xMrMXzeyduMeC8N9TCKbU\nDfBgafD1QKvvGm4DRrr7UcB/EPxyHpXqgFpAps1IyKjPXTjMNA34XtizUvN6terrV8v5Zcz1c/cq\ndz+CoKfvGIIVizNGzfMzs0OBH3twC42jgX0JbpEhaSzpRW7uPq7hvQD4A/CPZB+/ha3h80uG9g63\nZQx3Xxf++5GZPU7wy6yhVUpbmw1m1t3dN4R1HhtTHVAyuftHcT+26s9dWJg7jWA9lCfDzRlz/Wo7\nv0y6ftXcvdTMiglWKO5kZrGw1ygjfofGnd+E6towdy83s3uAy1ManDSopWel9Yj78Qzg3ZY8fjOY\nAww0s35mlgv8J/BUimNKGjNrF/71ipntA4yn9V8zCMb748f8nwIuCJ+fDzxZ8w2tzOfOL8M+d38G\nFrn77XHbMun67XV+mXL9zKxr9TCgmRUQ3E5qEcECgGeFu7Xa61fH+S2uvn5mZgT1b63y+rUlLT0r\n7X6CcfMqgqXLL6quDWitwqmztxMkmX9y9xtTHFLSmNmBwOMEQxPZwF9b+/mZ2YPAGIIu7Q3AtQQ3\nNnwE6AOsAs52909SFWMUdZzfl8iAz52ZjQSmAwvYc/uCq4A3gIdp5devnvM7h8y4focTFFfHwsff\n3f2X4e+Zh4DOwDzga+FkllalnvN7CehK8MfK28DFcUXakoa0wKOIiIhIKGOmYYuIiIhEpcRIRERE\nJKTESERERCSkxEhEREQkpMRIREREJKTESERERCSkxEhEREQkpMRIREREJPT/Z3yhvmcYzIkAAAAA\nSUVORK5CYII=\n",
       "text": [
        "<matplotlib.figure.Figure at 0x7f7cf731fad0>"
       ]
      }
     ],
     "prompt_number": 5
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "Rounded Least-square approximation"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import picos as pic\n",
      "P = pic.Problem()\n",
      "p = P.add_variable('p',M)\n",
      "P.minimize(pic.norm(A*p-Ides)**2,verbose=False)\n",
      "p_sol = p.value\n",
      "p_rlsq = cvx.matrix([min(max(pi,0),pmax) for pi in p_sol])"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stderr",
       "text": [
        "Freeing MIP data.\n"
       ]
      }
     ],
     "prompt_number": 6
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now, let's visualize the solution !"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "visualize_sol(p_rlsq)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 0
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "Constrained Least-square approximation (QP)"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import picos as pic\n",
      "P = pic.Problem()\n",
      "p = P.add_variable('p',M)\n",
      "#TODO\n",
      "p_lsq = p.value"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 0
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "visualize_sol(p_lsq)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 0
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "L1-approximation (LP)"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import picos as pic\n",
      "P = pic.Problem()\n",
      "p = P.add_variable('p',M)\n",
      "#TODO ...\n",
      "p_l1 = p.value"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 0
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "visualize_sol(p_l1)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 0
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "Convex-Programming solution"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import picos as pic\n",
      "P = pic.Problem()\n",
      "p = P.add_variable('p',M)\n",
      "I = A*p\n",
      "#TODO\n",
      "p_cvx = p.value"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 0
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "visualize_sol(p_cvx)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 0
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "LP solution + bisection"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import picos as pic\n",
      "\n",
      "\n",
      "\n",
      "alpha_lo = 1\n",
      "alpha_up = 10\n",
      "\n",
      "alpha = (alpha_lo+alpha_up)/2.\n",
      "\n",
      "while alpha_up-alpha_lo > 0.0001:\n",
      "    P = pic.Problem()\n",
      "    p = P.add_variable('p',M)\n",
      "\n",
      "    #TODO\n",
      "    if P.status == 'optimal':\n",
      "        alpha_up = alpha\n",
      "        p_lpb = p.value\n",
      "    else:\n",
      "        alpha_lo = alpha\n",
      "        \n",
      "    alpha = (alpha_lo+alpha_up)/2.\n",
      "    print alpha_lo,alpha_up\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 0
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "visualize_sol(p_lpb)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 0
    }
   ],
   "metadata": {}
  }
 ]
}