Chapter03.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 3 Nonlinear Equations and Complementarity Problems "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 194,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "\n",
    "import numpy as np\n",
    "from numpy import append, array, diagonal, tril, triu\n",
    "from numpy.linalg import inv\n",
    "from scipy.linalg import lu\n",
    "#from scipy.linalg import solve\n",
    "from pprint import pprint\n",
    "from numpy import array, zeros, diag, diagflat, dot\n",
    "\n",
    "from sympy import *\n",
    "import sympy as sym\n",
    "init_printing()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One of the most basic numerical operations encountered in computational economics\n",
    "is to find the solution of a system of nonlinear equations. Nonlinear equations generally\n",
    "arise in one of two forms. In the nonlinear *rootfinding problem*, a function f\n",
    "mapping $R^n$ to $R^n$ is given and one must compute an n-vector $x$, called a *root* of $f$,\n",
    "that satisfies\n",
    "\n",
    "$$f(x) = 0$$\n",
    "\n",
    "\n",
    "In the nonlinear fixed-point problem, a function $g$ from $R^n$ to $R^n$ is given and one\n",
    "must compute an n-vector x called a fixed-point of $g$, that satisfies\n",
    "\n",
    "$$x = g(x)$$\n",
    "\n",
    "\n",
    "The two forms are equivalent. The rootfinding problem may be recast as a fixed-point\n",
    "problem by letting $g(x) = x - f(x)$; conversely, the fixed-point problem may be recast\n",
    "as a rootfinding problem by letting $f(x) = x - g(x)$.\n",
    "\n",
    "\n",
    "In the related complementarity problem, two n-vectors $a$ and $b$, with $a < b$, and\n",
    "a function f from $R^n$ to $R^n$ are given, and one must compute an n-vector $x \\in [a; b]$,\n",
    "that satisfies\n",
    "\n",
    "\n",
    "$$x_i > a_i \\rightarrow f_i(x) \\forall i = 1,...,n$$\n",
    "\n",
    "$$x_i < b_i \\rightarrow f_i(x) \\forall i = 1,...,n$$\n",
    "\n",
    "\n",
    "The rootfinding problem is a special case of complementarity problem in which $a_i =\n",
    "-\\inf$ and $b_i = +\\inf$ for all i. However, the complementarity problem is not simply to\n",
    "find a root that lies within specified bounds. An element $f_i(x)$ may be nonzero at a\n",
    "solution of the complementarity problem, provided that $x_i$ equals one of the bounds\n",
    "$a_i$ or $b_i$.\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "https://github.com/QuantEcon/QuantEcon.lectures.code\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3.1 Bisection Method"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The bisection method is perhaps the simplest and most robust method for computing\n",
    "the root of a continuous real-valued function defined on a bounded interval of the real line.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The bisection method is an iterative procedure. Each iteration begins with an\n",
    "interval known to contain or to bracket a root of f, meaning the function has diffierent\n",
    "signs at the interval endpoints. The interval is bisected into two subintervals of equal\n",
    "length."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The bisection method's greatest strength is its robustness. In contrast to other\n",
    "rootfinding methods, the bisection method is guaranteed to compute a root to a\n",
    "prescribed tolerance in a known number of iterations, provided valid data are input.\n",
    "\n",
    "\n",
    "https://lectures.quantecon.org/py/scipy.html#roots-and-fixed-points"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "https://en.wikipedia.org/wiki/Bisection_method\n",
    "\n",
    "![](https://upload.wikimedia.org/wikipedia/commons/thumb/8/8c/Bisection_method.svg/375px-Bisection_method.svg.png)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 195,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#%pylab inline\n",
    "%matplotlib notebook\n",
    "# pylab Populating the interactive namespace from numpy and matplotlib\n",
    "# numpy for numerical computation\n",
    "# matplotlib for ploting\n",
    "from matplotlib import pyplot as plt\n",
    "from mpl_toolkits.mplot3d import Axes3D     # 3d"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 196,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "application/javascript": [
       "/* Put everything inside the global mpl namespace */\n",
       "window.mpl = {};\n",
       "\n",
       "mpl.get_websocket_type = function() {\n",
       "    if (typeof(WebSocket) !== 'undefined') {\n",
       "        return WebSocket;\n",
       "    } else if (typeof(MozWebSocket) !== 'undefined') {\n",
       "        return MozWebSocket;\n",
       "    } else {\n",
       "        alert('Your browser does not have WebSocket support.' +\n",
       "              'Please try Chrome, Safari or Firefox ≥ 6. ' +\n",
       "              'Firefox 4 and 5 are also supported but you ' +\n",
       "              'have to enable WebSockets in about:config.');\n",
       "    };\n",
       "}\n",
       "\n",
       "mpl.figure = function(figure_id, websocket, ondownload, parent_element) {\n",
       "    this.id = figure_id;\n",
       "\n",
       "    this.ws = websocket;\n",
       "\n",
       "    this.supports_binary = (this.ws.binaryType != undefined);\n",
       "\n",
       "    if (!this.supports_binary) {\n",
       "        var warnings = document.getElementById(\"mpl-warnings\");\n",
       "        if (warnings) {\n",
       "            warnings.style.display = 'block';\n",
       "            warnings.textContent = (\n",
       "                \"This browser does not support binary websocket messages. \" +\n",
       "                    \"Performance may be slow.\");\n",
       "        }\n",
       "    }\n",
       "\n",
       "    this.imageObj = new Image();\n",
       "\n",
       "    this.context = undefined;\n",
       "    this.message = undefined;\n",
       "    this.canvas = undefined;\n",
       "    this.rubberband_canvas = undefined;\n",
       "    this.rubberband_context = undefined;\n",
       "    this.format_dropdown = undefined;\n",
       "\n",
       "    this.image_mode = 'full';\n",
       "\n",
       "    this.root = $('<div/>');\n",
       "    this._root_extra_style(this.root)\n",
       "    this.root.attr('style', 'display: inline-block');\n",
       "\n",
       "    $(parent_element).append(this.root);\n",
       "\n",
       "    this._init_header(this);\n",
       "    this._init_canvas(this);\n",
       "    this._init_toolbar(this);\n",
       "\n",
       "    var fig = this;\n",
       "\n",
       "    this.waiting = false;\n",
       "\n",
       "    this.ws.onopen =  function () {\n",
       "            fig.send_message(\"supports_binary\", {value: fig.supports_binary});\n",
       "            fig.send_message(\"send_image_mode\", {});\n",
       "            fig.send_message(\"refresh\", {});\n",
       "        }\n",
       "\n",
       "    this.imageObj.onload = function() {\n",
       "            if (fig.image_mode == 'full') {\n",
       "                // Full images could contain transparency (where diff images\n",
       "                // almost always do), so we need to clear the canvas so that\n",
       "                // there is no ghosting.\n",
       "                fig.context.clearRect(0, 0, fig.canvas.width, fig.canvas.height);\n",
       "            }\n",
       "            fig.context.drawImage(fig.imageObj, 0, 0);\n",
       "        };\n",
       "\n",
       "    this.imageObj.onunload = function() {\n",
       "        this.ws.close();\n",
       "    }\n",
       "\n",
       "    this.ws.onmessage = this._make_on_message_function(this);\n",
       "\n",
       "    this.ondownload = ondownload;\n",
       "}\n",
       "\n",
       "mpl.figure.prototype._init_header = function() {\n",
       "    var titlebar = $(\n",
       "        '<div class=\"ui-dialog-titlebar ui-widget-header ui-corner-all ' +\n",
       "        'ui-helper-clearfix\"/>');\n",
       "    var titletext = $(\n",
       "        '<div class=\"ui-dialog-title\" style=\"width: 100%; ' +\n",
       "        'text-align: center; padding: 3px;\"/>');\n",
       "    titlebar.append(titletext)\n",
       "    this.root.append(titlebar);\n",
       "    this.header = titletext[0];\n",
       "}\n",
       "\n",
       "\n",
       "\n",
       "mpl.figure.prototype._canvas_extra_style = function(canvas_div) {\n",
       "\n",
       "}\n",
       "\n",
       "\n",
       "mpl.figure.prototype._root_extra_style = function(canvas_div) {\n",
       "\n",
       "}\n",
       "\n",
       "mpl.figure.prototype._init_canvas = function() {\n",
       "    var fig = this;\n",
       "\n",
       "    var canvas_div = $('<div/>');\n",
       "\n",
       "    canvas_div.attr('style', 'position: relative; clear: both; outline: 0');\n",
       "\n",
       "    function canvas_keyboard_event(event) {\n",
       "        return fig.key_event(event, event['data']);\n",
       "    }\n",
       "\n",
       "    canvas_div.keydown('key_press', canvas_keyboard_event);\n",
       "    canvas_div.keyup('key_release', canvas_keyboard_event);\n",
       "    this.canvas_div = canvas_div\n",
       "    this._canvas_extra_style(canvas_div)\n",
       "    this.root.append(canvas_div);\n",
       "\n",
       "    var canvas = $('<canvas/>');\n",
       "    canvas.addClass('mpl-canvas');\n",
       "    canvas.attr('style', \"left: 0; top: 0; z-index: 0; outline: 0\")\n",
       "\n",
       "    this.canvas = canvas[0];\n",
       "    this.context = canvas[0].getContext(\"2d\");\n",
       "\n",
       "    var rubberband = $('<canvas/>');\n",
       "    rubberband.attr('style', \"position: absolute; left: 0; top: 0; z-index: 1;\")\n",
       "\n",
       "    var pass_mouse_events = true;\n",
       "\n",
       "    canvas_div.resizable({\n",
       "        start: function(event, ui) {\n",
       "            pass_mouse_events = false;\n",
       "        },\n",
       "        resize: function(event, ui) {\n",
       "            fig.request_resize(ui.size.width, ui.size.height);\n",
       "        },\n",
       "        stop: function(event, ui) {\n",
       "            pass_mouse_events = true;\n",
       "            fig.request_resize(ui.size.width, ui.size.height);\n",
       "        },\n",
       "    });\n",
       "\n",
       "    function mouse_event_fn(event) {\n",
       "        if (pass_mouse_events)\n",
       "            return fig.mouse_event(event, event['data']);\n",
       "    }\n",
       "\n",
       "    rubberband.mousedown('button_press', mouse_event_fn);\n",
       "    rubberband.mouseup('button_release', mouse_event_fn);\n",
       "    // Throttle sequential mouse events to 1 every 20ms.\n",
       "    rubberband.mousemove('motion_notify', mouse_event_fn);\n",
       "\n",
       "    rubberband.mouseenter('figure_enter', mouse_event_fn);\n",
       "    rubberband.mouseleave('figure_leave', mouse_event_fn);\n",
       "\n",
       "    canvas_div.on(\"wheel\", function (event) {\n",
       "        event = event.originalEvent;\n",
       "        event['data'] = 'scroll'\n",
       "        if (event.deltaY < 0) {\n",
       "            event.step = 1;\n",
       "        } else {\n",
       "            event.step = -1;\n",
       "        }\n",
       "        mouse_event_fn(event);\n",
       "    });\n",
       "\n",
       "    canvas_div.append(canvas);\n",
       "    canvas_div.append(rubberband);\n",
       "\n",
       "    this.rubberband = rubberband;\n",
       "    this.rubberband_canvas = rubberband[0];\n",
       "    this.rubberband_context = rubberband[0].getContext(\"2d\");\n",
       "    this.rubberband_context.strokeStyle = \"#000000\";\n",
       "\n",
       "    this._resize_canvas = function(width, height) {\n",
       "        // Keep the size of the canvas, canvas container, and rubber band\n",
       "        // canvas in synch.\n",
       "        canvas_div.css('width', width)\n",
       "        canvas_div.css('height', height)\n",
       "\n",
       "        canvas.attr('width', width);\n",
       "        canvas.attr('height', height);\n",
       "\n",
       "        rubberband.attr('width', width);\n",
       "        rubberband.attr('height', height);\n",
       "    }\n",
       "\n",
       "    // Set the figure to an initial 600x600px, this will subsequently be updated\n",
       "    // upon first draw.\n",
       "    this._resize_canvas(600, 600);\n",
       "\n",
       "    // Disable right mouse context menu.\n",
       "    $(this.rubberband_canvas).bind(\"contextmenu\",function(e){\n",
       "        return false;\n",
       "    });\n",
       "\n",
       "    function set_focus () {\n",
       "        canvas.focus();\n",
       "        canvas_div.focus();\n",
       "    }\n",
       "\n",
       "    window.setTimeout(set_focus, 100);\n",
       "}\n",
       "\n",
       "mpl.figure.prototype._init_toolbar = function() {\n",
       "    var fig = this;\n",
       "\n",
       "    var nav_element = $('<div/>')\n",
       "    nav_element.attr('style', 'width: 100%');\n",
       "    this.root.append(nav_element);\n",
       "\n",
       "    // Define a callback function for later on.\n",
       "    function toolbar_event(event) {\n",
       "        return fig.toolbar_button_onclick(event['data']);\n",
       "    }\n",
       "    function toolbar_mouse_event(event) {\n",
       "        return fig.toolbar_button_onmouseover(event['data']);\n",
       "    }\n",
       "\n",
       "    for(var toolbar_ind in mpl.toolbar_items) {\n",
       "        var name = mpl.toolbar_items[toolbar_ind][0];\n",
       "        var tooltip = mpl.toolbar_items[toolbar_ind][1];\n",
       "        var image = mpl.toolbar_items[toolbar_ind][2];\n",
       "        var method_name = mpl.toolbar_items[toolbar_ind][3];\n",
       "\n",
       "        if (!name) {\n",
       "            // put a spacer in here.\n",
       "            continue;\n",
       "        }\n",
       "        var button = $('<button/>');\n",
       "        button.addClass('ui-button ui-widget ui-state-default ui-corner-all ' +\n",
       "                        'ui-button-icon-only');\n",
       "        button.attr('role', 'button');\n",
       "        button.attr('aria-disabled', 'false');\n",
       "        button.click(method_name, toolbar_event);\n",
       "        button.mouseover(tooltip, toolbar_mouse_event);\n",
       "\n",
       "        var icon_img = $('<span/>');\n",
       "        icon_img.addClass('ui-button-icon-primary ui-icon');\n",
       "        icon_img.addClass(image);\n",
       "        icon_img.addClass('ui-corner-all');\n",
       "\n",
       "        var tooltip_span = $('<span/>');\n",
       "        tooltip_span.addClass('ui-button-text');\n",
       "        tooltip_span.html(tooltip);\n",
       "\n",
       "        button.append(icon_img);\n",
       "        button.append(tooltip_span);\n",
       "\n",
       "        nav_element.append(button);\n",
       "    }\n",
       "\n",
       "    var fmt_picker_span = $('<span/>');\n",
       "\n",
       "    var fmt_picker = $('<select/>');\n",
       "    fmt_picker.addClass('mpl-toolbar-option ui-widget ui-widget-content');\n",
       "    fmt_picker_span.append(fmt_picker);\n",
       "    nav_element.append(fmt_picker_span);\n",
       "    this.format_dropdown = fmt_picker[0];\n",
       "\n",
       "    for (var ind in mpl.extensions) {\n",
       "        var fmt = mpl.extensions[ind];\n",
       "        var option = $(\n",
       "            '<option/>', {selected: fmt === mpl.default_extension}).html(fmt);\n",
       "        fmt_picker.append(option)\n",
       "    }\n",
       "\n",
       "    // Add hover states to the ui-buttons\n",
       "    $( \".ui-button\" ).hover(\n",
       "        function() { $(this).addClass(\"ui-state-hover\");},\n",
       "        function() { $(this).removeClass(\"ui-state-hover\");}\n",
       "    );\n",
       "\n",
       "    var status_bar = $('<span class=\"mpl-message\"/>');\n",
       "    nav_element.append(status_bar);\n",
       "    this.message = status_bar[0];\n",
       "}\n",
       "\n",
       "mpl.figure.prototype.request_resize = function(x_pixels, y_pixels) {\n",
       "    // Request matplotlib to resize the figure. Matplotlib will then trigger a resize in the client,\n",
       "    // which will in turn request a refresh of the image.\n",
       "    this.send_message('resize', {'width': x_pixels, 'height': y_pixels});\n",
       "}\n",
       "\n",
       "mpl.figure.prototype.send_message = function(type, properties) {\n",
       "    properties['type'] = type;\n",
       "    properties['figure_id'] = this.id;\n",
       "    this.ws.send(JSON.stringify(properties));\n",
       "}\n",
       "\n",
       "mpl.figure.prototype.send_draw_message = function() {\n",
       "    if (!this.waiting) {\n",
       "        this.waiting = true;\n",
       "        this.ws.send(JSON.stringify({type: \"draw\", figure_id: this.id}));\n",
       "    }\n",
       "}\n",
       "\n",
       "\n",
       "mpl.figure.prototype.handle_save = function(fig, msg) {\n",
       "    var format_dropdown = fig.format_dropdown;\n",
       "    var format = format_dropdown.options[format_dropdown.selectedIndex].value;\n",
       "    fig.ondownload(fig, format);\n",
       "}\n",
       "\n",
       "\n",
       "mpl.figure.prototype.handle_resize = function(fig, msg) {\n",
       "    var size = msg['size'];\n",
       "    if (size[0] != fig.canvas.width || size[1] != fig.canvas.height) {\n",
       "        fig._resize_canvas(size[0], size[1]);\n",
       "        fig.send_message(\"refresh\", {});\n",
       "    };\n",
       "}\n",
       "\n",
       "mpl.figure.prototype.handle_rubberband = function(fig, msg) {\n",
       "    var x0 = msg['x0'];\n",
       "    var y0 = fig.canvas.height - msg['y0'];\n",
       "    var x1 = msg['x1'];\n",
       "    var y1 = fig.canvas.height - msg['y1'];\n",
       "    x0 = Math.floor(x0) + 0.5;\n",
       "    y0 = Math.floor(y0) + 0.5;\n",
       "    x1 = Math.floor(x1) + 0.5;\n",
       "    y1 = Math.floor(y1) + 0.5;\n",
       "    var min_x = Math.min(x0, x1);\n",
       "    var min_y = Math.min(y0, y1);\n",
       "    var width = Math.abs(x1 - x0);\n",
       "    var height = Math.abs(y1 - y0);\n",
       "\n",
       "    fig.rubberband_context.clearRect(\n",
       "        0, 0, fig.canvas.width, fig.canvas.height);\n",
       "\n",
       "    fig.rubberband_context.strokeRect(min_x, min_y, width, height);\n",
       "}\n",
       "\n",
       "mpl.figure.prototype.handle_figure_label = function(fig, msg) {\n",
       "    // Updates the figure title.\n",
       "    fig.header.textContent = msg['label'];\n",
       "}\n",
       "\n",
       "mpl.figure.prototype.handle_cursor = function(fig, msg) {\n",
       "    var cursor = msg['cursor'];\n",
       "    switch(cursor)\n",
       "    {\n",
       "    case 0:\n",
       "        cursor = 'pointer';\n",
       "        break;\n",
       "    case 1:\n",
       "        cursor = 'default';\n",
       "        break;\n",
       "    case 2:\n",
       "        cursor = 'crosshair';\n",
       "        break;\n",
       "    case 3:\n",
       "        cursor = 'move';\n",
       "        break;\n",
       "    }\n",
       "    fig.rubberband_canvas.style.cursor = cursor;\n",
       "}\n",
       "\n",
       "mpl.figure.prototype.handle_message = function(fig, msg) {\n",
       "    fig.message.textContent = msg['message'];\n",
       "}\n",
       "\n",
       "mpl.figure.prototype.handle_draw = function(fig, msg) {\n",
       "    // Request the server to send over a new figure.\n",
       "    fig.send_draw_message();\n",
       "}\n",
       "\n",
       "mpl.figure.prototype.handle_image_mode = function(fig, msg) {\n",
       "    fig.image_mode = msg['mode'];\n",
       "}\n",
       "\n",
       "mpl.figure.prototype.updated_canvas_event = function() {\n",
       "    // Called whenever the canvas gets updated.\n",
       "    this.send_message(\"ack\", {});\n",
       "}\n",
       "\n",
       "// A function to construct a web socket function for onmessage handling.\n",
       "// Called in the figure constructor.\n",
       "mpl.figure.prototype._make_on_message_function = function(fig) {\n",
       "    return function socket_on_message(evt) {\n",
       "        if (evt.data instanceof Blob) {\n",
       "            /* FIXME: We get \"Resource interpreted as Image but\n",
       "             * transferred with MIME type text/plain:\" errors on\n",
       "             * Chrome.  But how to set the MIME type?  It doesn't seem\n",
       "             * to be part of the websocket stream */\n",
       "            evt.data.type = \"image/png\";\n",
       "\n",
       "            /* Free the memory for the previous frames */\n",
       "            if (fig.imageObj.src) {\n",
       "                (window.URL || window.webkitURL).revokeObjectURL(\n",
       "                    fig.imageObj.src);\n",
       "            }\n",
       "\n",
       "            fig.imageObj.src = (window.URL || window.webkitURL).createObjectURL(\n",
       "                evt.data);\n",
       "            fig.updated_canvas_event();\n",
       "            fig.waiting = false;\n",
       "            return;\n",
       "        }\n",
       "        else if (typeof evt.data === 'string' && evt.data.slice(0, 21) == \"data:image/png;base64\") {\n",
       "            fig.imageObj.src = evt.data;\n",
       "            fig.updated_canvas_event();\n",
       "            fig.waiting = false;\n",
       "            return;\n",
       "        }\n",
       "\n",
       "        var msg = JSON.parse(evt.data);\n",
       "        var msg_type = msg['type'];\n",
       "\n",
       "        // Call the  \"handle_{type}\" callback, which takes\n",
       "        // the figure and JSON message as its only arguments.\n",
       "        try {\n",
       "            var callback = fig[\"handle_\" + msg_type];\n",
       "        } catch (e) {\n",
       "            console.log(\"No handler for the '\" + msg_type + \"' message type: \", msg);\n",
       "            return;\n",
       "        }\n",
       "\n",
       "        if (callback) {\n",
       "            try {\n",
       "                // console.log(\"Handling '\" + msg_type + \"' message: \", msg);\n",
       "                callback(fig, msg);\n",
       "            } catch (e) {\n",
       "                console.log(\"Exception inside the 'handler_\" + msg_type + \"' callback:\", e, e.stack, msg);\n",
       "            }\n",
       "        }\n",
       "    };\n",
       "}\n",
       "\n",
       "// from http://stackoverflow.com/questions/1114465/getting-mouse-location-in-canvas\n",
       "mpl.findpos = function(e) {\n",
       "    //this section is from http://www.quirksmode.org/js/events_properties.html\n",
       "    var targ;\n",
       "    if (!e)\n",
       "        e = window.event;\n",
       "    if (e.target)\n",
       "        targ = e.target;\n",
       "    else if (e.srcElement)\n",
       "        targ = e.srcElement;\n",
       "    if (targ.nodeType == 3) // defeat Safari bug\n",
       "        targ = targ.parentNode;\n",
       "\n",
       "    // jQuery normalizes the pageX and pageY\n",
       "    // pageX,Y are the mouse positions relative to the document\n",
       "    // offset() returns the position of the element relative to the document\n",
       "    var x = e.pageX - $(targ).offset().left;\n",
       "    var y = e.pageY - $(targ).offset().top;\n",
       "\n",
       "    return {\"x\": x, \"y\": y};\n",
       "};\n",
       "\n",
       "/*\n",
       " * return a copy of an object with only non-object keys\n",
       " * we need this to avoid circular references\n",
       " * http://stackoverflow.com/a/24161582/3208463\n",
       " */\n",
       "function simpleKeys (original) {\n",
       "  return Object.keys(original).reduce(function (obj, key) {\n",
       "    if (typeof original[key] !== 'object')\n",
       "        obj[key] = original[key]\n",
       "    return obj;\n",
       "  }, {});\n",
       "}\n",
       "\n",
       "mpl.figure.prototype.mouse_event = function(event, name) {\n",
       "    var canvas_pos = mpl.findpos(event)\n",
       "\n",
       "    if (name === 'button_press')\n",
       "    {\n",
       "        this.canvas.focus();\n",
       "        this.canvas_div.focus();\n",
       "    }\n",
       "\n",
       "    var x = canvas_pos.x;\n",
       "    var y = canvas_pos.y;\n",
       "\n",
       "    this.send_message(name, {x: x, y: y, button: event.button,\n",
       "                             step: event.step,\n",
       "                             guiEvent: simpleKeys(event)});\n",
       "\n",
       "    /* This prevents the web browser from automatically changing to\n",
       "     * the text insertion cursor when the button is pressed.  We want\n",
       "     * to control all of the cursor setting manually through the\n",
       "     * 'cursor' event from matplotlib */\n",
       "    event.preventDefault();\n",
       "    return false;\n",
       "}\n",
       "\n",
       "mpl.figure.prototype._key_event_extra = function(event, name) {\n",
       "    // Handle any extra behaviour associated with a key event\n",
       "}\n",
       "\n",
       "mpl.figure.prototype.key_event = function(event, name) {\n",
       "\n",
       "    // Prevent repeat events\n",
       "    if (name == 'key_press')\n",
       "    {\n",
       "        if (event.which === this._key)\n",
       "            return;\n",
       "        else\n",
       "            this._key = event.which;\n",
       "    }\n",
       "    if (name == 'key_release')\n",
       "        this._key = null;\n",
       "\n",
       "    var value = '';\n",
       "    if (event.ctrlKey && event.which != 17)\n",
       "        value += \"ctrl+\";\n",
       "    if (event.altKey && event.which != 18)\n",
       "        value += \"alt+\";\n",
       "    if (event.shiftKey && event.which != 16)\n",
       "        value += \"shift+\";\n",
       "\n",
       "    value += 'k';\n",
       "    value += event.which.toString();\n",
       "\n",
       "    this._key_event_extra(event, name);\n",
       "\n",
       "    this.send_message(name, {key: value,\n",
       "                             guiEvent: simpleKeys(event)});\n",
       "    return false;\n",
       "}\n",
       "\n",
       "mpl.figure.prototype.toolbar_button_onclick = function(name) {\n",
       "    if (name == 'download') {\n",
       "        this.handle_save(this, null);\n",
       "    } else {\n",
       "        this.send_message(\"toolbar_button\", {name: name});\n",
       "    }\n",
       "};\n",
       "\n",
       "mpl.figure.prototype.toolbar_button_onmouseover = function(tooltip) {\n",
       "    this.message.textContent = tooltip;\n",
       "};\n",
       "mpl.toolbar_items = [[\"Home\", \"Reset original view\", \"fa fa-home icon-home\", \"home\"], [\"Back\", \"Back to  previous view\", \"fa fa-arrow-left icon-arrow-left\", \"back\"], [\"Forward\", \"Forward to next view\", \"fa fa-arrow-right icon-arrow-right\", \"forward\"], [\"\", \"\", \"\", \"\"], [\"Pan\", \"Pan axes with left mouse, zoom with right\", \"fa fa-arrows icon-move\", \"pan\"], [\"Zoom\", \"Zoom to rectangle\", \"fa fa-square-o icon-check-empty\", \"zoom\"], [\"\", \"\", \"\", \"\"], [\"Download\", \"Download plot\", \"fa fa-floppy-o icon-save\", \"download\"]];\n",
       "\n",
       "mpl.extensions = [\"eps\", \"jpeg\", \"pdf\", \"png\", \"ps\", \"raw\", \"svg\", \"tif\"];\n",
       "\n",
       "mpl.default_extension = \"png\";var comm_websocket_adapter = function(comm) {\n",
       "    // Create a \"websocket\"-like object which calls the given IPython comm\n",
       "    // object with the appropriate methods. Currently this is a non binary\n",
       "    // socket, so there is still some room for performance tuning.\n",
       "    var ws = {};\n",
       "\n",
       "    ws.close = function() {\n",
       "        comm.close()\n",
       "    };\n",
       "    ws.send = function(m) {\n",
       "        //console.log('sending', m);\n",
       "        comm.send(m);\n",
       "    };\n",
       "    // Register the callback with on_msg.\n",
       "    comm.on_msg(function(msg) {\n",
       "        //console.log('receiving', msg['content']['data'], msg);\n",
       "        // Pass the mpl event to the overriden (by mpl) onmessage function.\n",
       "        ws.onmessage(msg['content']['data'])\n",
       "    });\n",
       "    return ws;\n",
       "}\n",
       "\n",
       "mpl.mpl_figure_comm = function(comm, msg) {\n",
       "    // This is the function which gets called when the mpl process\n",
       "    // starts-up an IPython Comm through the \"matplotlib\" channel.\n",
       "\n",
       "    var id = msg.content.data.id;\n",
       "    // Get hold of the div created by the display call when the Comm\n",
       "    // socket was opened in Python.\n",
       "    var element = $(\"#\" + id);\n",
       "    var ws_proxy = comm_websocket_adapter(comm)\n",
       "\n",
       "    function ondownload(figure, format) {\n",
       "        window.open(figure.imageObj.src);\n",
       "    }\n",
       "\n",
       "    var fig = new mpl.figure(id, ws_proxy,\n",
       "                           ondownload,\n",
       "                           element.get(0));\n",
       "\n",
       "    // Call onopen now - mpl needs it, as it is assuming we've passed it a real\n",
       "    // web socket which is closed, not our websocket->open comm proxy.\n",
       "    ws_proxy.onopen();\n",
       "\n",
       "    fig.parent_element = element.get(0);\n",
       "    fig.cell_info = mpl.find_output_cell(\"<div id='\" + id + \"'></div>\");\n",
       "    if (!fig.cell_info) {\n",
       "        console.error(\"Failed to find cell for figure\", id, fig);\n",
       "        return;\n",
       "    }\n",
       "\n",
       "    var output_index = fig.cell_info[2]\n",
       "    var cell = fig.cell_info[0];\n",
       "\n",
       "};\n",
       "\n",
       "mpl.figure.prototype.handle_close = function(fig, msg) {\n",
       "    fig.root.unbind('remove')\n",
       "\n",
       "    // Update the output cell to use the data from the current canvas.\n",
       "    fig.push_to_output();\n",
       "    var dataURL = fig.canvas.toDataURL();\n",
       "    // Re-enable the keyboard manager in IPython - without this line, in FF,\n",
       "    // the notebook keyboard shortcuts fail.\n",
       "    IPython.keyboard_manager.enable()\n",
       "    $(fig.parent_element).html('<img src=\"' + dataURL + '\">');\n",
       "    fig.close_ws(fig, msg);\n",
       "}\n",
       "\n",
       "mpl.figure.prototype.close_ws = function(fig, msg){\n",
       "    fig.send_message('closing', msg);\n",
       "    // fig.ws.close()\n",
       "}\n",
       "\n",
       "mpl.figure.prototype.push_to_output = function(remove_interactive) {\n",
       "    // Turn the data on the canvas into data in the output cell.\n",
       "    var dataURL = this.canvas.toDataURL();\n",
       "    this.cell_info[1]['text/html'] = '<img src=\"' + dataURL + '\">';\n",
       "}\n",
       "\n",
       "mpl.figure.prototype.updated_canvas_event = function() {\n",
       "    // Tell IPython that the notebook contents must change.\n",
       "    IPython.notebook.set_dirty(true);\n",
       "    this.send_message(\"ack\", {});\n",
       "    var fig = this;\n",
       "    // Wait a second, then push the new image to the DOM so\n",
       "    // that it is saved nicely (might be nice to debounce this).\n",
       "    setTimeout(function () { fig.push_to_output() }, 1000);\n",
       "}\n",
       "\n",
       "mpl.figure.prototype._init_toolbar = function() {\n",
       "    var fig = this;\n",
       "\n",
       "    var nav_element = $('<div/>')\n",
       "    nav_element.attr('style', 'width: 100%');\n",
       "    this.root.append(nav_element);\n",
       "\n",
       "    // Define a callback function for later on.\n",
       "    function toolbar_event(event) {\n",
       "        return fig.toolbar_button_onclick(event['data']);\n",
       "    }\n",
       "    function toolbar_mouse_event(event) {\n",
       "        return fig.toolbar_button_onmouseover(event['data']);\n",
       "    }\n",
       "\n",
       "    for(var toolbar_ind in mpl.toolbar_items){\n",
       "        var name = mpl.toolbar_items[toolbar_ind][0];\n",
       "        var tooltip = mpl.toolbar_items[toolbar_ind][1];\n",
       "        var image = mpl.toolbar_items[toolbar_ind][2];\n",
       "        var method_name = mpl.toolbar_items[toolbar_ind][3];\n",
       "\n",
       "        if (!name) { continue; };\n",
       "\n",
       "        var button = $('<button class=\"btn btn-default\" href=\"#\" title=\"' + name + '\"><i class=\"fa ' + image + ' fa-lg\"></i></button>');\n",
       "        button.click(method_name, toolbar_event);\n",
       "        button.mouseover(tooltip, toolbar_mouse_event);\n",
       "        nav_element.append(button);\n",
       "    }\n",
       "\n",
       "    // Add the status bar.\n",
       "    var status_bar = $('<span class=\"mpl-message\" style=\"text-align:right; float: right;\"/>');\n",
       "    nav_element.append(status_bar);\n",
       "    this.message = status_bar[0];\n",
       "\n",
       "    // Add the close button to the window.\n",
       "    var buttongrp = $('<div class=\"btn-group inline pull-right\"></div>');\n",
       "    var button = $('<button class=\"btn btn-mini btn-primary\" href=\"#\" title=\"Stop Interaction\"><i class=\"fa fa-power-off icon-remove icon-large\"></i></button>');\n",
       "    button.click(function (evt) { fig.handle_close(fig, {}); } );\n",
       "    button.mouseover('Stop Interaction', toolbar_mouse_event);\n",
       "    buttongrp.append(button);\n",
       "    var titlebar = this.root.find($('.ui-dialog-titlebar'));\n",
       "    titlebar.prepend(buttongrp);\n",
       "}\n",
       "\n",
       "mpl.figure.prototype._root_extra_style = function(el){\n",
       "    var fig = this\n",
       "    el.on(\"remove\", function(){\n",
       "\tfig.close_ws(fig, {});\n",
       "    });\n",
       "}\n",
       "\n",
       "mpl.figure.prototype._canvas_extra_style = function(el){\n",
       "    // this is important to make the div 'focusable\n",
       "    el.attr('tabindex', 0)\n",
       "    // reach out to IPython and tell the keyboard manager to turn it's self\n",
       "    // off when our div gets focus\n",
       "\n",
       "    // location in version 3\n",
       "    if (IPython.notebook.keyboard_manager) {\n",
       "        IPython.notebook.keyboard_manager.register_events(el);\n",
       "    }\n",
       "    else {\n",
       "        // location in version 2\n",
       "        IPython.keyboard_manager.register_events(el);\n",
       "    }\n",
       "\n",
       "}\n",
       "\n",
       "mpl.figure.prototype._key_event_extra = function(event, name) {\n",
       "    var manager = IPython.notebook.keyboard_manager;\n",
       "    if (!manager)\n",
       "        manager = IPython.keyboard_manager;\n",
       "\n",
       "    // Check for shift+enter\n",
       "    if (event.shiftKey && event.which == 13) {\n",
       "        this.canvas_div.blur();\n",
       "        // select the cell after this one\n",
       "        var index = IPython.notebook.find_cell_index(this.cell_info[0]);\n",
       "        IPython.notebook.select(index + 1);\n",
       "    }\n",
       "}\n",
       "\n",
       "mpl.figure.prototype.handle_save = function(fig, msg) {\n",
       "    fig.ondownload(fig, null);\n",
       "}\n",
       "\n",
       "\n",
       "mpl.find_output_cell = function(html_output) {\n",
       "    // Return the cell and output element which can be found *uniquely* in the notebook.\n",
       "    // Note - this is a bit hacky, but it is done because the \"notebook_saving.Notebook\"\n",
       "    // IPython event is triggered only after the cells have been serialised, which for\n",
       "    // our purposes (turning an active figure into a static one), is too late.\n",
       "    var cells = IPython.notebook.get_cells();\n",
       "    var ncells = cells.length;\n",
       "    for (var i=0; i<ncells; i++) {\n",
       "        var cell = cells[i];\n",
       "        if (cell.cell_type === 'code'){\n",
       "            for (var j=0; j<cell.output_area.outputs.length; j++) {\n",
       "                var data = cell.output_area.outputs[j];\n",
       "                if (data.data) {\n",
       "                    // IPython >= 3 moved mimebundle to data attribute of output\n",
       "                    data = data.data;\n",
       "                }\n",
       "                if (data['text/html'] == html_output) {\n",
       "                    return [cell, data, j];\n",
       "                }\n",
       "            }\n",
       "        }\n",
       "    }\n",
       "}\n",
       "\n",
       "// Register the function which deals with the matplotlib target/channel.\n",
       "// The kernel may be null if the page has been refreshed.\n",
       "if (IPython.notebook.kernel != null) {\n",
       "    IPython.notebook.kernel.comm_manager.register_target('matplotlib', mpl.mpl_figure_comm);\n",
       "}\n"
      ],
      "text/plain": [
       "<IPython.core.display.Javascript object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<img src=\"data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAoAAAAHgCAYAAAA10dzkAAAgAElEQVR4Xu3dCbhVZdn/8S+IICriazmPZTkrDpAaKCkYmiO+aqm8ouXrhDmgpmk5S+ZAloKaplmaaQ6ZU+RYhn8VxQElcS4py/Q1wBRB4H89rH3yeDznsPe597y++7r2xRHWvdZ6Puve8HPttZ7VDV8KKKCAAgoooIACuRLolqvROlgFFFBAAQUUUEABDIA2gQIKKKCAAgookDMBA2DODrjDVUABBRRQQAEFDID2gAIKKKCAAgookDMBA2DODrjDVUABBRRQQAEFDID2gAIKKKCAAgookDMBA2DODrjDVUABBRRQQAEFDID2gAIKKKCAAgookDMBA2DODrjDVUABBRRQQAEFDID2gAIKKKCAAgookDMBA2DODrjDVUABBRRQQAEFDID2gAIKKKCAAgookDMBA2DODrjDVUABBRRQQAEFDID2gAIKKKCAAgookDMBA2DODrjDVUABBRRQQAEFDID2gAIKKKCAAgookDMBA2DODrjDVUABBRRQQAEFDID2gAIKKKCAAgookDMBA2DODrjDVUABBRRQQAEFDID2gAIKKKCAAgookDMBA2DODrjDVUABBRRQQAEFDID2gAIKKKCAAgookDMBA2DODrjDVUABBRRQQAEFDID2gAIKKKCAAgookDMBA2DODrjDVUABBRRQQAEFDID2gAIKKKCAAgookDMBA2DODrjDVUABBRRQQAEFDID2gAIKKKCAAgookDMBA2DODrjDVUABBRRQQAEFDID2gAIKKKCAAgookDMBA2DODrjDVUABBRRQQAEFDID2gAIKKKCAAgookDMBA2DODrjDVUABBRRQQAEFDID2gAIKKKCAAgookDMBA2DODrjDVUABBRRQQAEFDID2gAIKKKCAAgookDMBA2DODrjDVUABBRRQQAEFDID2gAIKKKCAAgookDMBA2DODrjDVUABBRRQQAEFDID2gAIKKKCAAgookDMBA2DODrjDVUABBRRQQAEFDID2gAIKKKCAAgookDOBPAbA7wE7A2sC7wK/B74FTG917OcDs4EPgWS0ANgaeC5n/eFwFVBAAQUUUKAJBfIYAM8BbgKmAEsClwIbAJu1CYBDgAea8Jg7JAUUUEABBRTIuUAeA2DbQ94PmAwsB8wo/GE6AzgUuD/n/eHwFVBAAQUUUKAJBQyA2de/hwGfbXMG8O/A4sCfgcuAK5vw+DskBRRQQAEFFMihQN4DYDrLdyuwJ3BPq+O/HfAwMA/YAbgO+DZweZseSX6rALNy2DsOWQEFFFBAgUYW6AP8rXCdfyOPo0v7nucAuAvwc2Ak8JtF6J0KfBkY1Ga5VdvcPNKlg2CRAgoooIACCtREYDXgrzXZco03mtcAuD9wCbA3cG8RxyAFwGHAwDbLLpOuG3z99ddZZpn0o6+IwMknn8yYMWMiq7AW0LF8baClluUTKM+a7Mm44003wamnzuSvf109rawvMDO+1sZbQx4D4JHAmcCuwMR2Dlm6Gzi5pLuE080g6W7g64EUAse1FwBnzJhhACxD748ePZqxY8eWYU35XoWO5Tv+WmpZPoHyrMmejDnOnw/9+sHXvz6T0aNT9jMAxkQbqzqFurnAB4Xdbpnnb6dCIExfDZ8HpNPCaR7AdBPIeOCKdoa58AygAbA8DeBfbDqWR6B8a7EntSyfQHnWZE/GHG+7DQ4/HJ56aiYrrmgAjGnmu9oAWMbjP2HCBIYNS9+0+4oI6BjR+3itllqWT6A8a7Inu+64YAF84Quw337wjW/MpG9fA2DXNa00ANoDCiiggAIKNIDAhAkwYgS89hrMm2cAzOM1gOVsUwNgOTVdlwIKKKCAAhUS2HZb2Gkn+Pa3YeZMA6ABMNZoBsCYn9UKKKCAAgpUXOAPf4Ddd4c//xnSpB0GwOxuV19dFzAAdt3OSgUUUEABBaoi8OUvw1ZbwZlpDpA074tnAA2Awc4zAAYBLVdAAQUUUKCSAo89BkOGZNf+fepTBsAWa88AxrrOABjzs1oBBRRQQIGKCuy2G6y7Lpx//keb8QygXwFHm84AGBW0XgEFFFBAgQoJPPkkDBoEr7wCK65oAGzN7BnAWNMZAGN+ViuggAIKKFAxgeHD4TOfgbYPmfIMoGcAo01nAIwKWq+AAgoooEAFBJ5+GrbeGl5+GVZe+eMbMAAaAKMtZwCMClqvgAIKKKBABQT23jsLfj/60SdXbgA0AEZbzgAYFbReAQUUUECBMgs8+ywMGAAvvQSrrmoAbI/XawBjTWcAjPlZrYACCiigQNkFvva1bMqXcePaX7VnAD0DGG06A2BU0HoFFFBAAQXKKPCnP8Fmm8ELL8AaaxgAO6L1DGCs6QyAMT+rFVBAAQUUKKvAiBGw1FJw+eUdr9YzgJ4BjDadATAqaL0CCiiggAJlEkhn/TbZBJ5/HtZaywDYGatnAGNNZwCM+VmtgAIKKKBA2QQOOAAWXxx+8pPOV+kZQM8ARpvOABgVtF4BBRRQQIEyCEybBv36QboGME3+3NnLAGgAjLacATAqaL0CCiiggAJlEEjX/vXuDVdcseiVGQANgIvuks6XMABGBa1XQAEFFFAgKJCu+dt000Vf+9eyGQOgATDYchgAo4LWK6CAAgooEBTYbz9Yemn48Y+LW5EB0ABYXKd0vJQBMCpovQIKKKCAAgGBqVNh880hXQO45prFrcgAaAAsrlMMgFEn6xVQQAEFFKiIQHrqx7LLwmWXFb96A6ABsPhuaX9JzwBGBa1XQAEFFFCgiwLPPQf9+2dn/zp66kd7qzYAGgC72HL/KTMARgWtV0ABBRRQoIsC++wDn/40jB9f2goMgAbA0jrmk0sbAKOC1iuggAIKKNAFgSlT4AtfgBdfhNVWK20FBkADYGkdYwCMelmvgAIKKKBAWQT22gtWWgkuuaT01RkADYCld83HKzwDGBW0XgEFFFBAgRIFnnwSBg6El16CVVYpsRgwABoAS+8aA2DUzHoFFFBAAQVCArvuCuusAxde2LXVGAANgF3rnI+qPAMYFbReAQUUUECBEgQeeQSGDoVXXoEVViihsNWiBkADYNc6xwAYdbNeAQUUUECBLgnssAMMGABjxnSpfGGRAdAA2PXuySo9AxgVtF4BBRRQQIEiBX7/e9htN3j1VVhuuSKL2lnMAGgA7Hr3GACjdtYroIACCihQtMCCBTB4MAwZAqedVnRZuwsaAA2AsQ7yDGDUz3oFFFBAAQWKErjnHkiPfUvX/vXtW1RJhwsZAA2AsQ4yAEb9rFdAAQUUUGCRAuns31ZbwfDhcNJJi1x8kQsYAA2Ai2ySRSzgNYBRQesVUEABBRRYhMDtt8PBB2dn/5ZaKs5lADQARrvIABgVtF4BBRRQQIFOBObPh803h5Ej4dhjy0NlADQARjvJABgVtF4BBRRQQIFOBH75SzjhhOyZv0ssUR4qA6ABMNpJBsCooPUKKKCAAgp0IDB3LmywAZx4YvYVcLleBkADYLSXDIBRQesVUEABBRToQOCKK+D882HqVOjRo3xMBkADYLSbDIBRQesVUEABBRRoR2D2bPjc57Ln/X71q+UlMgAaAKMdZQCMClqvgAIKKKBAOwJjx8LPfw5PPAHdu5eXyABoAIx2lAEwKmi9AgoooIACbQRmzoTPfhZ+9jP4ylfKz2MANABGu8oAGBW0XgEFFFBAgTYCZ54Jv/sdPPQQdOtWfh4DoAEw2lUGwKig9QoooIACCrQSeOut7OzfnXfCNttUhsYAaACMdpYBMCpovQIKKKCAAq0Ejj8ennsO7r67ciwGQANgtLsMgFFB6xVQQAEFFCgI/OUvsN56MHEibLZZ5VgMgAbAaHcZAKOC1iuggAIKKFAQ+PrXIU3/8otfVJbEAGgAjHaYATAqaL0CCiiggAJkX/v27w/PPgtrr11ZEgOgATDaYQbAqKD1CiiggAIKALvvDquvDpdcUnkOA6ABMNplBsCooPUKKKCAArkXSNf87bgjvPQSrLhi5TkMgAbAaJcZAKOC1iuggAIK5FpgwQLYdlvYfns444zqUBgADYDRTjMARgWtV0ABBRTItcDtt0O6+ePll2GZ9K9qFV4GQANgtM0MgFFB6xVQQAEFciswbx706wf/+79w9NHVYzAAGgCj3WYAjApar4ACCiiQW4FrroHTT4fnn4devarHYADMZwD8HrAzsCbwLvB74FvA9FattwlwMbAF8C/gCqC9KxMMgNX7vLolBRRQQIEmEnj/fVh3XRgzBkaMqO7ADID5DIDnADcBU4AlgUuBDYCWOceXBl4ArgLOBNYB0gNpLgB+2KZFDYDV/cy6NQUUUECBJhE47zz45S/h8cehe/fqDsoAmM8A2LbL+gGTgeWAGcBI4PvAKsD8wsJHAd8EPm8ArO6H1K0poIACCjSfwNtvZ5M933wzDBlS/fEZAA2AqevS17+HAZ8ttOBYYH1gp1YtuTXwR6Bv4Wvjlj/yDGD1P7duUQEFFFCgwQWOPRamTYO77qrNQAyABsChwK3AnsA9hTa8ElgK2LdVW66XnlIDrA78rdXvGwBr89l1qwoooIACDSrwyiuw4Ybw2GOw8ca1GYQBMN8BcBfg54WvfH/TqgVLPgM4atQoevbsuXAVw4YNW/j2pYACCiiggAKfFNh3X+jdG65KV9pX8TVhwgTSO73mzJnDuHHj0o/pm72ZVdyNutlUt7rZk+ruyP5Aetrg3sC9bTZ9AHBem2sA0+xER3oNYHUPkltTQAEFFGgugUmT4Etfyr7+XW212o3NM4D5PAOYgly6u3dXYGI77ZfuAp5WuAs43TGcbvy4E7jQu4Br92F1ywoooIACjS2QHvm23XYwcCCck/51reHLAJjPAJju7J0LfFDovXQWdEHhpo+WQLgRML4wD2C6MzhNFXNWO73qNYA1/AC7aQUUUECBxhG44w446KDqPvKtIx0DYD4DYDk/LQbAcmq6LgUUUECBphSYOxc22QRGjYIj0/dwNX4ZAA2A0RY0AEYFrVdAAQUUaHqB8ePhRz+CKVNg8cVrP1wDoAEw2oUGwKig9QoooIACTS0wYwZ87nPZXb+7pqvv6+BlADQARtvQABgVtF4BBRRQoKkFTjopm/PvvvugW53MPWIANABGP3QGwKig9QoooIACTSvw2muwwQYwcSJstln9DNMAaACMdqMBMCpovQIKKKBA0wqkSZ979YKf/rS+hmgANABGO9IAGBW0XgEFFFCgKQUefRSGDMkmfV511foaogHQABjtSANgVNB6BRRQQIGmE0iTPg8aBEOHwhln1N/wDIAGwGhXGgCjgtYroIACCjSdwE03wTe/CS++CEun52vV2csAaACMtqQBMCpovQIKKKBAUwnMng3rrw/f/S58/ev1OTQDoAEw2pkGwKig9QoooIACTSVw7rlw440waRIstlh9Ds0AaACMdqYBMCpovQIKKKBA0wi88Qassw7ceSdsu239DssAaACMdqcBMCpovQIKKKBA0wh84xswcyb86lf1PSQDoAEw2qEGwKig9QoooIACTSEweXJ25+/UqbDWWvU9JAOgATDaoQbAqKD1CiiggAINL5CmfRk8GLbZBs45p/6HYwA0AEa71AAYFbReAQUUUKDhBdJXvkcdBS+8AH361P9wDIAGwGiXGgCjgtYroIACCjS0wPvvZ9O+nH46HHhgYwzFAGgAjHaqATAqaL0CCiigQEMLjBkDt9wCjz0G3bs3xlAMgAbAaKcaAKOC1iuggAIKNKzA9Omw3nrw299mN4A0yssAaACM9qoBMCpovQIKKKBAwwrsvz906wbXXttYQzAAGgCjHWsAjApar4ACCijQkAJ//CPstBM8/zysumpjDcEAaACMdqwBMCpovQIKKKBAwwnMmwf9+8M++8C3v91wu48B0AAY7VoDYFTQegUUUECBhhO4/HI4/3x49llYYomG230DIAbAaNcaAKOC1iuggAIKNJTAO+/A5z8PV18Nu+7aULv+n531DKABMNq5BsCooPUKKKCAAg0lkCZ8fvFFuOuu7AaQRnwZAA2A0b41AEYFrVdAAQUUaBiBKVPgC1+AJ5/Mpn9p1JcB0AAY7V0DYFTQegUUUECBhhBIz/vdfnvYfHO48MKG2OUOd9IAaACMdrABMCpovQIKKKBAQwhcfz2MHg3TpsEy6V+/Bn4ZAA2A0fY1AEYFrVdAAQUUqHuBWbOyr3zTnb/77Vf3u7vIHTQAGgAX2SSLWMAAGBW0XgEFFFCg7gVOOAEmTYIHHmjcGz9aIxsADYDRD50BMCpovQIKKKBAXQtMnQpbbJEFwI02qutdLXrnDIAGwKKbpYMFDYBRQesVUEABBepWIN34MWQIbLopjB1bt7tZ8o4ZAA2AJTdNmwIDYFTQegUUUECBuhW44QY45pjmuPHDr4A/3mYNOoVj3XxWDIB1cyjcEQUUUECBcgq03Phx3nmw//7lXHPt1+UZQM8ARrvQABgVtF4BBRRQoC4FvvUtePRRePDB5rjxwzOAngEs5wfNAFhOTdelgAIKKFAXAs8+mz3xIwXAjTeui10q6054BtAzgNGGMgBGBa1XQAEFFKgrgXTjx+DBsOWW2bx/zfgyABoAo31tAIwKWq+AAgooUFcCP/0pnHoqpOlfll66rnatbDtjADQARpvJABgVtF4BBRRQoG4E3n47e+LHj38Mw4fXzW6VfUcMgAbAaFMZAKOC1iuggAIK1I3AoYfC9Olwxx3Nd+NHa2QDoAEw+qEzAEYFrVdAAQUUqAuBRx7JJn2eMgU++9m62KWK7YQB0AAYbS4DYFTQegUUUECBmgt8+CH07w977w2nnFLz3an4DhgADYDRJjMARgWtV0ABBRSoucBFF8Fll8HTT0OvXjXfnYrvgAHQABhtMgNgVNB6BRRQQIGaCrz+OmywAdx2G2y/fU13pWobNwAaAKPNZgCMClqvgAIKKFBTgT32gL594ZprarobVd24AdAAGG04A2BU0HoFFFBAgZoJ/PrXcPDB8Kc/wfLL12w3qr5hA6ABMNp0BsCooPUKKKCAAjURmDUL1l8fzj4bDjywJrtQs40aAA2A0eYzAEYFrVdAAQUUqInA0UfDM8/A/fc395x/7eEaAA2A0Q+dATAqaL0CCiigQNUFJk3Knvf75JOw7rpV33zNN2gANABGm9AAGBW0XgEFFFCgqgJpzr8BAyDd/HHaaVXddN1szABoAIw2owEwKmi9AgoooEBVBcaOzZ71m5c5//wKuP326lbVrmu+jRkAm++YOiIFFFCgaQVefRU23hjuvDP7CjivL88AegYw2vsGwKig9QoooIACVRFYsAB23BHWWgsuv7wqm6zbjRgADYDR5jQARgWtV0ABBRSoisDPfgYnnQRTp8Kyy1Zlk3W7EQOgATDanAbAqKD1CiiggAIVF3jzzWzOvyuvhOHDK765ut+AATCfAfCrwCigH7A0sDgwv1W3pp9nAx+S+SwAtgaea6ejDYB1/zF3BxVQQAEF9tsPPvgAbr5ZiyRgAMxnANwBWA5YEriygwA4BHigiI+JAbAIJBdRQAEFFKidQLrhY//9s8e9rbxy7fajnrZsAMxnAGzpwXT/0/0dBMChhT9bVL8aABcl5J8roIACCtRMID3ubcMN4dRTs2f++soEDIAGwI4C4N8LwfDPwGWFM4XtfW4MgP5tooACCihQtwJHHgnPPZfPx711dlAMgAbA9gLgdsDDwDwgfV18HfBtoL2b5g2AdfvXnjumgAIK5FvgoYeyaV+eego+//l8W7QdvQHQANheAGzbJ6cCXwYGtfPxWRgAR40aRc+ePRf+8bBhwxa+fSmggAIKKFArgffeg3794PDDYfToWu1FfW13woQJpHd6zZkzh3HjxqUf+6ZvhOtrT6uzN3l+EkhH1wC2FwBTohvYUQCcMWMGyyyTsqAvBRRQQAEFai9wwgmQzgBOnAiLLVb7/am3PfAMYD7PAHYvXN+XAuDdQJ/C171zgE0LU79MKUwNk+4Gvh5IZwEX/q9Cm5dfAdfbp9r9UUABBXIu8OijsN128PjjsMEGOcfoYPgGwHwGwJHA1YX5/VJrtMz1l679S4HuPGC1wjyA6SaQ8cAVHfSQAdC/WxRQQAEF6kYgzfW3+ebZtC8nn1w3u1V3O2IAzGcALGcjGgDLqem6FFBAAQVCAt/9LqR5/9JZwMXTYw58tStgADQARj8aBsCooPUKKKCAAmURePJJGDgQHn4YNk0XNPnqUMAAaACMfjwMgFFB6xVQQAEFwgJz5sCWW8Iuu8BZZ4VX1/QrMAAaAKNNbgCMClqvgAIKKBAWOP10uPVWmDQJCrOShdfZzCswABoAo/1tAIwKWq+AAgooEBKYPBkGDcqmfNlss9CqclNsADQARpvdABgVtF4BBRRQoMsC6a7f/v1hr73gtNO6vJrcFRoADYDRpjcARgWtV0ABBRTossApp8Ddd3vXb6mABkADYKk903Z5A2BU0HoFFFBAgS4JpOv9Bg/Owt/GG3dpFbktMgAaAKPNbwCMClqvgAIKKFCywOzZ2YTPI0Y44XPJeOnhvzNn0rdvegywzwLuip812ZNDZvgsYFtBAQUUUKCaAt/6Fjz4YDbnX48e1dxyc2zLAOgZwGgnGwCjgtYroIACCpQk8Mc/wrBh2bN+11+/pFIXLggYAA2A0Q+DATAqaL0CCiigQNECs2ZBv35w1FFwzDFFl7lgGwEDoAEw+qEwAEYFrVdAAQUUKFrgkEPgpZfg3nuhe/eiy1zQAPiJHuhmV4QEDIAhPosVUEABBYoVuPNO2G8/mDIF1lij2CqXa0/AM4CeAYx+MgyAUUHrFVBAAQUWKfDWW7DRRvD978PIkYtc3AUWIWAANABGPyQGwKig9QoooIACnQosWAB77w3z5sEtt0A3v7sLd4wB0AAYbSIDYFTQegUUUECBTgWuuw5Gj4Znn4XllxerHAIGQANgtI8MgFFB6xVQQAEFOhT485+zu36vuQZ2312ocgkYAA2A0V4yAEYFrVdAAQUUaFcgfeW7/faw7rrw4x+LVE4BA6ABMNpPBsCooPUKKKCAAu0KnHsuXHUVPPkkLLWUSOUUMAAaAKP9ZACMClqvgAIKKPAJgSeegG22gd//HgYMEKjcAgZAA2C0pwyAUUHrFVBAAQU+JvDee7D55jBiBHznO+JUQsAAaACM9pUBMCpovQIKKKDAxwSOOAKeeSY7+7fYYuJUQsAAaACM9pUBMCpovQIKKKDAfwTuuAP23x+eego+8xlhKiVgADQARnvLABgVtF4BBRRQYKHAG29kU75ccAEccIAolRQwABoAo/1lAIwKWq+AAgoowPz5sOOO2UTP117r0z4q3RIGQANgtMcMgFFB6xVQQAEFFp71Gz8+++p3mfQvi6+KChgADYDRBjMARgWtV0ABBXIu8PjjsO22cP/9sNVWOceo0vANgAbAaKsZAKOC1iuggAI5Fpg1K5vy5aCD4OSTcwxR5aEbAA2A0ZYzAEYFrVdAAQVyLHDggZCe93vvvU75Us02MAAaAKP9ZgCMClqvgAIK5FTgF7+Ao46Cp5+GVVfNKUKNhm0ANABGW88AGBW0XgEFFMihwEsvwRZbwM9+BrvvnkOAGg/ZAGgAjLagATAqaL0CCiiQM4EPPoAvfjF71u9FF+Vs8HUyXAOgATDaigbAqKD1CiigQM4Ejj4aJk7M3r165WzwdTJcA6ABMNqKBsCooPUKKKBAjgR+/WsYORImT4a1187RwOtsqAZAA2C0JQ2AUUHrFVBAgZwIpLt9N9sMLrsM9tknJ4Ou02EaAA2A0dY0AEYFrVdAAQVyIDB3LgweDJtskgVAX7UVMAAaAKMdaACMClqvgAIK5EDgpJPgrrvg0Uehd+8cDLjOh2gANABGW9QAGBW0XgEFFGhygTvugP32g8ceg/XWa/LBNsjwDIAGwGirGgCjgtYroIACTSzw2mvZo97GjYN9923igTbY0AyABsBoyxoAo4LWK6CAAk0qkOb7S3P99e8P48c36SAbdFgGQANgtHUNgFFB6xVQQIEmFfjmN+H//T/n+6vHw2sANABG+9IAGBW0XgEFFGhCgRtugMMOy+b7+8xnmnCADT4kA6ABMNrCBsCooPUKKKBAkwlMmwYDBsC118JuuzXZ4JpkOAZAA2C0lQ2AUUHrFVBAgSYS+Pe/YautYMcd4fzzm2hgTTYUA6ABMNrSBsCooPUKKKBAkwgsWAAjRsD06XDffdCjR5MMrAmHYQA0AEbb2gAYFbReAQUUaBKBSy6BMWOy6/5WWqlJBtWkwzAAGgCjrW0AjApar4ACCjSBwMMPww47wO9+BwMHNsGAmnwIBkADYLTFDYBRQesVUECBBhf4xz+yyZ7T497S1C++6l/AAGgAjHapATAqaL0CCijQwAIffghDh8Kqq2Z3/Xbr1sCDydGuGwANgNF2NwBGBa1XQAEFGljg+ONhwgR45BFYaqkGHkjOdt0AaACMtrwBMCpovQIKKNCgAtdfD4cfDo89Buus06CDyOluGwANgNHWNwBGBa1XQAEFGlDgqadg0CC48Ub4ylcacAA532UDoAEw+hEwAEYFrVdAAQUaTODtt6F/fzj4YDjllAbbeXd3oYAB0AAY/SgYAKOC1iuggAINJJBu+khP+ejbF266yZs+GujQfWxXDYD5DIBfBUYB/YClgcWB+a06YxPgYmAL4F/AFcAZHTS5AbBRP/3utwIKKNAFgeOOg9/+Nrvpo0+fLqzAkroQMADmMwDuACwHLAlc2SYApkD4AnAVcCaQLuu9G7gA+GE7XWsArIuPsjuhgAIKVF7gF7+AUaOymz4+//nKb88tVE7AAJjPANjSUYOB+9sEwJHA94FVWp0VPApIU3u293E3AFbu8+maFVBAgboRePxx+NKXvOmjbg5IcEcMgAbAtgFwLLA+sFOr3toa+CPQF3i3Tc8ZAIMfQssVUJ+YFCMAACAASURBVECBehd44w0YMACOPhpOOKHe99b9K0bAAGgAbBsA01fCaSrPfVs10HrAc8DqwN8MgMV8tFxGAQUUaA6B2bOzM39pnr9rrvGmj+Y4qt4FnI5jnh9a095XwF06Azhq1Ch69uy58HMxbNiwhW9fCiiggAKNLbBgAYwcCS+8AA8+CEss0djjyfveT5gwgfROrzlz5jBu3Lj0Y/p2b2YebQyAH78L+ADgvDbXAB4NHOk1gHn8eDhmBRTIs8AFF8BFF8GkSbDyynmWaL6x+xVwPs8Adi/c+JHOAKY7fNON/PPS/xAUvv6dVrgL+JxC6LsTuNC7gJvvLwBHpIACCnQkcNddsM8+2Zm/NOmzr+YSMADmMwCmO32vBhYU2jmdBU0/bwf8AdgIGF+YB3AGcClwVget700gzfV3gqNRQAEFePZZGDgQLr8cvvY1QZpRwACYzwBYzl42AJZT03UpoIACNRZ4803Ycsvs2r/TT6/xzrj5igkYAA2A0eYyAEYFrVdAAQXqRCDd8bv99rDmmpAmfe6W56vk6+SYVGo3DIAGwGhvGQCjgtYroIACdSCQ7vgdMQJeeQUeeMA7fuvgkFR0FwyABsBogxkAo4LWK6CAAnUgcNZZcOWV2WPeVlyxDnbIXaiogAHQABhtMANgVNB6BRRQoMYCN9wAhxwCf/wjbLxxjXfGzVdFwABoAIw2mgEwKmi9AgooUEOBiRPTBP4+47eGh6AmmzYAGgCjjWcAjApar4ACCtRI4MUXYeut4cwz4YgjarQTbrYmAgZAA2C08QyAUUHrFVBAgRoIvPVWFv722APOP78GO+AmaypgADQARhvQABgVtF4BBRSoskCa7mXIkOzxbjfeCN3T86F85UrAAGgAjDa8ATAqaL0CCihQRYH582HffeEvf4H774fevau4cTdVNwIGQANgtBkNgFFB6xVQQIEqCpx4Itx0EzzyCCy/fBU37KbqSsAAaACMNqQBMCpovQIKKFAlgUsuyR7v9vDDsM46Vdqom6lLAQOgATDamAbAqKD1CiigQBUEbr45e77vvffCVltVYYNuoq4FDIAGwGiDGgCjgtYroIACFRZ46CHYaafs+b677Vbhjbn6hhAwABoAo41qAIwKWq+AAgpUUGDqVBg0CM49N3vahy8FkoAB0AAY/SQYAKOC1iuggAIVEvjrX7O5/g46CM44o0IbcbUNKWAANABGG9cAGBW0XgEFFKiAwDvvwODBMGAAXHkldOtWgY24yoYVMAAaAKPNawCMClqvgAIKlFngvffgy1+G5ZaDW26BHj3KvAFX1/ACBkADYLSJDYBRQesVUECBMgrMnQvDh6drvGDCBCd6LiNtU63KAGgAjDa0ATAqaL0CCihQJoH0lI8DD4RnnoEHH4Rlly3Til1N0wkYAA2A0aY2AEYFrVdAAQXKILBgARx3HNx2G0ycCCutVIaVuoqmFTAAGgCjzW0AjApar4ACCpRB4Hvfgx/+MAt/a69dhhW6iqYWMAAaAKMNbgCMClqvgAIKBAXGj4dTToEHHoBNNw2uzPJcCBgADYDRRjcARgWtV0ABBQICP/85HHFEdsPHF78YWJGluRIwABoAow1vAIwKWq+AAgp0UeDWW2HEiOy6v6FDu7gSy3IpYAA0AEYb3wAYFbReAQUU6ILAPffAHnvAdddlv/pSoBQBA6ABsJR+aW9ZA2BU0HoFFFCgRIF0o8eOO8Kll2ZnAH0pUKqAAdAAWGrPtF3eABgVtF4BBRQoQeCJJ2DIEEh3/R5+eAmFLqpAKwEDoAEw+oEwAEYFrVdAAQWKFHj6adh+++yO39GjiyxyMQXaETAAGgCjHwwDYFTQegUUUKAIgalTYfDgLPh9+9tFFLiIAp0IGAANgNEPiAEwKmi9AgoosAiBF17Iwt8hh8AZZ8ilQFzAAGgAjHaRATAqaL0CCijQicArr8C222Y3e6Tr/rp1k0uBuIAB0AAY7SIDYFTQegUUUKADgddegy99CYYPh7FjDX82SvkEDIAGwGg3GQCjgtYroIAC7Qi0hL9dd4Uf/cjwZ5OUV8AAaACMdpQBMCpovQIKKNBGwPBnS1RawABoAIz2mAEwKmi9Agoo0ErA8Gc7VEPAAGgAjPaZATAqaL0CCihQEDD82QrVEjAAGgCjvWYAjApar4ACCgDpbt80ybPX/NkO1RAwABoAo31mAIwKWq+AArkXSPP8pce7/fd/ww9+4A0fuW+IKgAYAA2A0TYzAEYFrVdAgVwLpCd8pPA3cqTz/OW6Eao8eAOgATDacgbAqKD1CiiQW4FnnoGhQ+Hww+H00z3zl9tGqMHADYAGwGjbGQCjgtYroEAuBSZPhh12gOOOg5NPziWBg66hgAHQABhtPwNgVNB6BRTIncDDD8POO8Mpp8Dxx+du+A64DgQMgAbAaBsaAKOC1iugQK4E7rkH9twTzjsv++rXlwK1EDAAGgCjfWcAjApar4ACuRH49a9hxAi4/HLYf//cDNuB1qGAAdAAGG1LA2BU0HoFFMiFwLXXwmGHwXXXwe6752LIDrKOBQyABsBoexoAo4LWK6BA0wtceimceCLcems25YsvBWotYAA0AEZ70AAYFbReAQWaVmDBAjj77Gxy5zvvhK23btqhOrAGEzAAGgCjLWsAjApar4ACTSkwfz4cfTTccgtMmAAbbdSUw3RQDSpgADQARlvXABgVtF4BBZpOYM6c7MkeTzwBv/sdrLVW0w3RATW4gAHQABhtYQNgVNB6BRRoKoF3382e6fv223DXXbDCCk01PAfTJAIGQANgtJUNgFFB6xVQoGkE/vlP2GUX6NMnu+Ej/epLgXoUMAAaAKN9aQCMClqvgAJNIfDyy7DjjrDFFnDNNdCrV1MMy0E0qYAB0AAYbW0DYFTQegUUaHiBSZOyR7sdcED2hI/u3Rt+SA6gyQUMgAbAaIsbAKOC1iugQEMLpOldvvY1OOssOOaYhh6KO58jAQOgAbC9dj8N+C7wHpnPAuB2oL0HFxkAc/QXhkNVQIGPC/zkJ9lUL1dfDXvvrY4CjSNgADQAdhQA01z12xbRygbAIpBcRAEFmksgzfH33e/C+PGQnu87eHBzjc/RNL+AAdAAaABs/s+5I1RAgTIKzJ4NBx4Ijz2WTfOy3nplXLmrUqBKAgZAA2BHAfD4wlfA6Wvgh4FTgNfaWdgzgFX6sLoZBRSovUCa5mWPPSA94u2222D55Wu/T+6BAl0RMAAaANvrmw2AWcDrwMrA+cBWwCaFUNi6xgDYlU+eNQoo0HAC06bBV77y0TQvvXs33BDcYQX+I2AANAAW83HoCcwAdgXubVOwMACOGjWKnj3TYjBs2LCFb18KKKBAswjcd192k8ehh8I55zjNS7Mc17yNY8KECaR3es2ZM4dx48alH/sCM/Nmkcab7nL11blASwDcDbinvQA4Y8YMllkmZUFfCiigQHMJXHopHH88pH8r07V/vhRoBgHPABoA2+vjNJnB/cDbwIqFr4AHARsD/zYANsNH3zEooMCiBD78MJvX74Ybsse6DUp/C/pSoEkEDIAGwPZa+bbCNX9LAe8AfyjMC/hKOwt7DWCT/GXgMBRQ4COBd96BffaBv/8dbr8d1lpLHQWaS8AAaACMdrQBMCpovQIK1JXA88/D7rvDuuvCdddBnz51tXvujAJlETAAGgCjjWQAjApar4ACdSOQzvaNGAFHHAFnnw2LLVY3u+aOKFBWAQOgATDaUAbAqKD1CihQc4H0ZI8xY+D734f0eLf09a8vBZpZwABoAIz2twEwKmi9AgrUVODdd2HkSJg8OXusW79+Nd0dN65AVQQMgAbAaKMZAKOC1iugQM0EXnwR9twTVlgBbrwRPvWpmu2KG1agqgIGQANgtOEMgFFB6xVQoCYC6VFu6czfwQfDuedCjx412Q03qkBNBAyABsBo4xkAo4LWK6BAVQXmzYNTT4WLL4arroK99qrq5t2YAnUhYAA0AEYb0QAYFbReAQWqJvDWW7DvvjB9OtxyC6y/ftU27YYUqCsBA6ABMNqQBsCooPUKKFAVgUceye7u3XLL7Myf8/tVhd2N1KmAAdAAGG1NA2BU0HoFFKiowIIFcNFF8J3vZHP7pce7dfMp8BU1d+X1L2AANABGu9QAGBW0XgEFKiaQHun29a9nU7ykZ/putVXFNuWKFWgoAQOgATDasAbAqKD1CihQEYHHH4e994aNNoJrroHllqvIZlypAg0pYAA0AEYb1wAYFbReAQXKKtD6K9/TToPjj4fu3cu6CVemQMMLGAANgNEmNgBGBa1XQIGyCbz5Jhx0EDz3HFx/PWy9ddlW7YoUaCoBA6ABMNrQBsCooPUKKFAWgfvugxEjYJtt4Mc/hmWXLctqXYkCTSlgADQARhvbABgVtF4BBUICc+bA6afDj34EP/hB9mQP7/INkVqcAwEDoAEw2uYGwKig9Qoo0GWBadNg//1h7tzsK98NNujyqixUIFcCBkADYLThDYBRQesVUKBkgXSjx+WXZzd4HHYYnHMO9OpV8mosUCC3AgZAA2C0+Q2AUUHrFVCgJIF0o8c3vgFPPplN7zJkSEnlLqyAAoAB0AAY/SAYAKOC1iugQNECt90GhxwCgwfDZZc5t1/RcC6oQBsBA6ABMPqhMABGBa1XQIFFCvzrX3D00fCb38DFF2fX/XmjxyLZXECBDgUMgAbA6MfDABgVtF4BBToVuOee7HFuG24IV14Jq60mmAIKRAUMgAbAaA8ZAKOC1iugQLsC774LJ54IP/sZXHBB9tWvZ/1sFgXKI2AANABGO8kAGBW0XgEFPiFw773ZfH5rrglXXQVrry2SAgqUU8AAaACM9pMBMCpovQIK/EdgxoxsapcbboDvfx8OPdTn+NoeClRCwABoAIz2lQEwKmi9AgosFLjzzizwpcmcr7giO/vnSwEFKiNgADQARjvLABgVtF6BnAv84x9wzDFw990wdiwcdJDX+uW8JRx+FQQMgAbAaJsZAKOC1iuQU4H0NI+rr86+8h06FH74Q1h55ZxiOGwFqixgADQARlvOABgVtF6BHAq88EL2de9LL8H48bDrrjlEcMgK1FDAAGgAjLafATAqaL0CORKYPRvOPRfOOy+7yzc9w7dPnxwBOFQF6kTAAGgAjLaiATAqaL0CORH47W/hyCNh2WXh0kthwICcDNxhKlCHAgZAA2C0LQ2AUUHrFWhygenT4dhjIT3RY8yY7KvfxRZr8kE7PAXqXMAAaACMtqgBMCpovQJNKvDBB3DRRXD22TB8OJx/Pqy4YpMO1mEp0GACBkADYLRlDYBRQesVaEKBu+7Kpnbp3Rsuvhi23bYJB+mQFGhgAQOgATDavgbAqKD1CjSRQLqrNwW/hx/Ozvyl5/f26NFEA3QoCjSJgAHQABhtZQNgVNB6BZpAID3CLQW+ceNg5Mjs5099qgkG5hAUaFIBA6ABMNraBsCooPUKNLDAhx9mj2079VTYZJPsSR79+jXwgNx1BXIiYAA0AEZb3QAYFbRegQYVSNO6HHcczJ0LF1yQTebcrVuDDsbdViBnAgZAA2C05Q2AUUHrFWgwgcmT4cQT4Ykn4LTT4PDDoWfPBhuEu6tAzgUMgAbA6EfAABgVtF6BBhF49VX4znfg1lvhm9+Ek06C//qvBtl5d1MBBT4mYAA0AEY/EgbAqKD1CtS5wD//mU3gfNll8LWvwRlnwBpr1PlOu3sKKNCpgAHQABj9iBgAo4LWK1CnAunO3nRTR3qnefy+973sRg9fCijQ+AIGQANgtIsNgFFB6xWoM4H33sumczn3XNhoo+zs38CBdbaT7o4CCoQEDIAGwFADAQbAqKD1CtSJwOzZ2ZQu6UzfKqtkwW+HHbyzt04Oj7uhQFkFDIAGwGhDGQCjgtYrUGOB1sFvhRWyO3v32MPgV+PD4uYVqKiAAdAAGG0wA2BU0HoFaiTw/vvwk59kZ/w+/Wk4/XTYfXfo3r1GO+RmFVCgagIGQANgtNkMgFFB6xWossCsWdkdvRdeCCuu+NEZP4NflQ+Em1OghgIGQANgtP0MgFFB6xWoksA778DFF8NFF8G668Ipp8DOO/tVb5X43YwCdSVgADQARhvSABgVtF6BCgtMnw4//CFcfjn0759N5rzddga/CrO7egXqWsAAaACMNqgBMCpovQIVEnjuuewZvddfn53pO+EE2GqrCm3M1SqgQEMJGAANgNGGNQBGBa1XoIwCCxbAgw9mkzffey/8z//AccdlX/n6UkABBVoEDIAGwOinwQAYFbRegTIIfPAB/PKX8IMfwOuvw2GHZc/rXWmlMqzcVSigQNMJGAANgNGmNgBGBa1XICDw5pvZtX3jx0PfvnDssdlZvyWXDKzUUgUUaHoBA6ABMNrkBsCooPUKlCiQvuZ97DG45BL41a+y5/QecwzsuKNz+JVI6eIK5FbAAGgAjDa/ATAqaL0CRQqkiZtvvDELftOmwYEHwhFHwHrrFbkCF1NAAQUKAgZAA2BnH4YzgIMLz/t9AhgFPNemwADoXycKVFjg+eezr3mvuSa7pu/II7Ovefv0qfCGXb0CCjStgAHQANhRc58AHAnsBLwMnAYcAKwDvNeqyADYtH89OLBaCqTn8956axb8HnkE9toLDj0UBg1y/r5aHhe3rUCzCBgADYAd9fIrwFjgksICiwF/A0YD1xkAm+WvAMdRbwJPPQVXXQXXXgvLLw+HHAIjR2bP6vWlgAIKlEvAAGgAbK+X0lm9fwFbA4+2WmACMAU43gBYro+g61EA/u//silcfvITSF/37rMPfOMbMHCgZ/vsDwUUqIyAAdAA2F5nrQb8BVgfmNZqgV8CM4FD2gbA1998nWWWSbnRlwIKFCMwZw7cc0/2lI7f/hY23TS7rm/PPb22rxg/l1FAgZhACoCrr7B6Wknfwr/tsRU2YHW3BtznSu9yyWcAGQCkL4nT63OFd6X30vUroIACCiigQPECLwHpnV7zgEkLfzIAFi+YiyXbuwbwDeDY9q4B9AxgLnrCQXZRYOrUbL6+m2+Gf/0LdtsN9t0Xtt7aefu6SGqZAgoEBTwD6FfAHbVQus4v3QW8M5DC4KnACCA9UdS7gIMfPMubXyDN05dC3w03wMsvfxT60mTNvXo1//gdoQIK1LeA1wAaADvr0NOBQ4E029jjzgNY3x9m9672Ai2hLwW/9POwYbD33rD77l7XV/uj4x4ooEBrAQOgATD6iXAewKig9Q0rkB7JNnlyNl9feqczfekMXwp9u+4K3hfVsIfWHVeg6QUMgAbAaJMbAKOC1jeUQLp79w9/gN/8Bn796+yavp13huHDYaedPNPXUAfTnVUgxwIGQANgtP0NgFFB6+te4K234O674fbbYcIEWHpp2GUX2GMP2H57r+mr+wPoDiqgwCcEDIAGwOjHwgAYFbS+7gTmz4cnnshCX3pPmpTN05dCX/pqd/PNnaC57g6aO6SAAiUJGAANgCU1TDsLGwCjgtbXhcAbb2QTM//ud9lZvrlzYYcdsq9103V9q6xSF7vpTiiggAJlETAAGgCjjWQAjApaXxOBWbPgoYey0Jfef/oTbLEFDB2ahb40R1+PHjXZNTeqgAIKVFzAAGgAjDaZATAqaH1VBN57DyZOhAceyN6PPw5rrJEFvnSmL13Lt9xyVdkVN6KAAgrUXMAAaACMNqEBMCpofUUE0t25KfCls3zpna7jW3ll2G677P2lL8Gaa1Zk065UAQUUqHsBA6ABMNqkBsCooPVhgTQf36uvwsMPZ+8U/KZMgc99DrbZJntvuy185jPevBHGdgUKKNAUAgZAA2C0kQ2AUUHrSxaYOTO7S/exx+CRR7LQ98472d25X/xi9h40CFZaqeRVW6CAAgrkQsAAaACMNroBMCpofacC6dq9p5/OnriRrttLoS/dsLHaavCFL8CWW2aBL4W/3r3FVEABBRQoRsAAaAAspk86W8YAGBW0/j8C6SxeCntPPZW901m+FPb+67+yO3TTOwW+AQOy6/l8KaCAAgp0TcAAaADsWud8VGUAjArmsD7NsffCC/Dss9m1eumdAt9f/pLdmduvXzbxcgp86cxeOtvXrVsOoRyyAgooUCEBA6ABMNpaBsCoYBPXf/BBFvTSWbyW99Sp8Pzz0LMnbLghbLRR9k6BLwU/p2Jp4oZwaAooUDcCBkADYLQZDYBRwQavnzcP/vpXePHFLOxNm5b9mt7pztwll4T11oMNNoD118/eG28Ma60F3bs3+ODdfQUUUKBBBQyABsBo6xoAo4INUP/uu1mYa/1+6SV4+WV45RX48MNsTr1114V11vn4r3592wAH2F1UQIHcCRgADYDRpjcARgVrXJ/C29//np3Fe/317Dq81u/XXoO334Y+fbJ59NKZu/RrmmNv7bWzX1P4S1/p+lJAAQUUaAwBA6ABMNqpBsCoYIXq040Wb74J//hH9v7b3+CNNz7+6/TpWfhLEymvsEJ2s0UKc+lGjJZf088p8KVr87wRo0IHy9UqoIACVRYwABoAoy1nAIwKFlE/fz6kyY//7/+ys3Gt3//8J6T3W29lv7aEvjSlSnp96lOw4orZtCmrrJL92vJzCnyrrpr9t2fwijgQLqKAAgo0iYAB0AAYbWUDYCeC6cxamsj43/+GdB1des+a9dE7hbr03zNmfPROv5eeY5sCXHqnn9M7hcAePbIzcSnUtfy6/PLQ9p0CX3qn31988eghtl4BBRRQoNkEDIAGwGhPVyUApiCVAlB6t/45/Xe6C7Xlz9LP7b3TdW6t3+nr0fbec+ZA63eaxqT1e/bs7L/ffx/Sz+nX1u8U9lreKfSld8trscVg6aWza+mWWSb7teXdty+0fafJj9N72WWzX9M71fk1bLRlrVdAAQUUMAAaAKOfgoUBsHfvGXTrln7s+JWCW3q1/Nryc9vfbwl46dfWy5a6oykopdCV3uksWDp71vqdfq/l91t+7tUr+yo0vVt+Tr+2vJdYIvs5PXIs/Zx+bfl5qaWyKU9av9PvpdCXfk3rNLyVehRdXgEFFFCgEgIGQANgtK8WBsBJk2aw9NIfD4DthZ2W32v9Z+nntr+f5odr+f2WX1OQSz+3/Fn6Nf1e+rX1O4W8lmWjg7NeAQUUUECBZhQwABoAo31dla+AoztpvQIKKKCAAgp8JGAANABGPw8GwKig9QoooIACClRZwABoAIy2nAEwKmi9AgoooIACVRYwABoAoy1nAIwKWq+AAgoooECVBQyABsBoyxkAo4LWK6CAAgooUGUBA6ABMNpyBsCooPUKKKCAAgpUWcAAaACMtpwBMCpovQIKKKCAAlUWMAAaAKMtZwCMClqvgAIKKKBAlQUMgAbAaMsZAKOC1iuggAIKKFBlAQOgATDacgbAqKD1CiiggAIKVFnAAGgAjLacATAqaL0CCiiggAJVFjAAGgCjLWcAjApar4ACCiigQJUFDIAGwGjLGQCjgtYroIACCihQZQEDoAEw2nIGwKig9QoooIACClRZwABoAIy2nAEwKmi9AgoooIACVRYwABoAoy1nAIwKWq+AAgoooECVBQyABsBoyxkAo4LWK6CAAgooUGUBA6ABMNpyBsCooPUKKKCAAgpUWcAAaACMtpwBMCpovQIKKKCAAlUWMAAaAKMtZwCMClqvgAIKKKBAlQUMgAbAaMsZAKOC1iuggAIKKFBlAQOgATDacgbAqKD1CiiggAIKVFnAAGgAjLacATAqaL0CCiiggAJVFjAAGgCjLWcAjApar4ACCiigQJUFDIAGwGjLGQCjgtYroIACCihQZQEDoAEw2nIGwKig9QoooIACClRZwABoAIy2nAEwKmi9AgoooIACVRYwABoAoy1nAIwKWq+AAgoooECVBQyABsBoyxkAo4LWK6CAAgooUGUBA6ABMNpyBsCooPUKKKCAAgpUWcAAaACMtpwBMCpovQIKKKCAAlUWMAAaAKMtZwCMClqvgAIKKKBAlQUMgAbAti03GHgAeLfwB92Ad4A1OuhNA2CVP7RuTgEFFFBAgaiAAdAA2F4AvB/oASwoosEMgEUgFbvIhAkTGDZsWLGLu1wHAjqWrzW01LJ8AuVZkz1ZHkcDoAGwowDYE5hXRJsZAItAKnaR0aNHM3bs2GIXd7kOBHQsX2toqWX5BMqzJnuyPI4GQANgRwFwOtALmAKcBfyhg5YzAJbns7hwLf7FVh5MHcvjaE+Wz1HL8ln6+S6PpQEwPwHwamBk4WvddF1f29eDwPbAisAKwHNAb+Aw4GxgS+CZduoWBsDXX3+dZZZJP/qKCJx88smMGTMmsgprAR3L1wZaalk+gfKsyZ4sj2MKgKuvvnpaWV9gZnnW2lhraS8MNdYIitvbJYElOll0LjCrgz9P1wQ+DHynnT9fFUhnC30poIACCiigQOMJrAb8tfF2O77HeQmAEan7gEeAU9pZSfJbpZPwGNmutQoooIACCihQOYE+wN+KvOmzcntRozUbAD8O/2XgBeDPhTOGhwLfAwYCk2t0jNysAgoooIACCihQVgED4Mc509e8/wssB7xfuAnkTOD3ZVV3ZQoooIACCiigQA0FDIA1xHfTCiiggAIKKKBALQQMgItWPwM4GEi3+T4BjCrcJdxe5bLAOGBnYD5wJ3BkulN40ZvJxRKlWKY7s7cGPiC7Wz1NzP0t4LJcSHU8yK8WerAfsDSweKHXOqqwJ8tnaU+2b5kuk0l/561ZeIpS+sYkfVY7u0HOviyfpX35SctTgQOATwNzCv92nwQ83cm/H7nrSQNg52nihEKA2wl4GTit0FTrAO+1U5oCX/oHOf0jnWxvAP4N7JHz0JKGX6pleiRfmn8xmfv6SGCHwiUK6c72K4sIgPZkx91TqqU92b7lOcBNhUtmUl9eCmwAbNbJB9e+LJ+lfflJy88DbxZOvqQnex1V+J+SlTu54SN3MTHvYQAABYtJREFUPWkA7DxavAKkR1NcUlhsscIdQ6OB69qUpucFvwZsAjxb+LP081OFZwnnfbqYUiwTX/pL7SEg/Z+cr08KpOdWpymKOjsDaE8W1znFWNqTxVmmpdLZ6XTTXLqWur1vP+zL8lnal4u2TA91OBy4sDDP79vtlOSyJw2AHTdP+sr3X4WvIR9ttdiEwv/pHt+mdDfgl0D6P+DWr9nAXsAdi+7Tpl2iVMuWv9Q2AroD/wBuK0zKnc6o+oJiQos9WVynFGNpTxZnmZZKX/+mSfQ/20GJfVk+S/uyY8uvFE7UpIme0yVZPyh8E9VeRS570gDYcfOkySH/AqwPTGu1WAp5adbwQ9qUjgDOB9Ip5tavv6ennAG/KP4z33RLlmqZALYCni+E8I2BawrHYd+m0+nagIoJLfZkcbbFWNqTxVkOBW4F9gTu6aDEviyfpX25aMt0bV96Elj6Fu5me/IjAQNgx81T6lmrXP4fxKI/ewuXKNWyvdVuC9wLpIk7040heX8VE1rsyeK6pBhLe3LRlrsAPy/8Y/ubTha3L8tnaV8u2jItkbLOO8A2hW/w2lblsicNgJ03T3vXrb0BHNvBNYCvFq5/abkGsOVamHR3nNcAfvJ6yo4sO/tLLYXJ9LV63l/FhJZ0XYs9uehOKcbSnuzccf/CtdJ7F/5HrbOl7cvyWdqXi/58pyXSjSDpetT/AW5ppySXPWkA7Lx50nV+aRqXNMVBCoPphoT09cW6HdwFfHvhovz0l2Gyvb4wLcLw4nq0qZcqxXKFwh2E6SaQdLf1hsBPCzfZpH9g8vxK10SmGz9SaLm7cEZ0XmGqgzRVTtuXPdlxt5RiaU927Jj+jkwT5u8KTCzyw2lftg9VqqV92b5juus3Xa6V7gReHkh3qqdr8dMlXema8vZeuetJA+Ci/7Y6HUiPhEtfPT7eah7A1YGpwI6t/tJL1xqkO4bTVyHpH+PUUOkDna4Z9AXFWqb/G/sVkKbbSXdep+so07UbZxem1cmzZbqW5epWUxm0zJG4XeFsnz1ZfHeUYmlPduyaLrCf2+rSjJaeTNNnpUDo35XF92SplvZlx2Guf2Gu1PTv76TC/6S0PNLVniycpSq+NV1SAQUUUEABBRRQoOEFPAPY8IfQASiggAIKKKCAAqUJGABL83JpBRRQQAEFFFCg4QUMgA1/CB2AAgoooIACCihQmoABsDQvl1ZAAQUUUEABBRpewADY8IfQASiggAIKKKCAAqUJGABL83JpBRRQQAEFFFCg4QUMgA1/CB2AAgoooIACCihQmoABsDQvl1ZAAQUUUEABBRpewADY8IfQASiggAIKKKCAAqUJGABL83JpBRRQQAEFFFCg4QUMgA1/CB2AAgoooIACCihQmoABsDQvl1ZAAQUUUEABBRpewADY8IfQASiggAIKKKCAAqUJGABL83JpBRRQQAEFFFCg4QUMgA1/CB2AAgoooIACCihQmoABsDQvl1ZAAQUUUEABBRpewADY8IfQASiggAIKKKCAAqUJGABL83JpBRRQQAEFFFCg4QUMgA1/CB2AAgoooIACCihQmoABsDQvl1ZAAQUUUEABBRpewADY8IfQASiggAIKKKCAAqUJGABL83JpBRRQQAEFFFCg4QUMgA1/CB2AAgoooIACCihQmoABsDQvl1ZAAQUUUEABBRpewADY8IfQASiggAIKKKCAAqUJGABL83JpBRRQQAEFFFCg4QUMgA1/CB2AAgoooIACCihQmoABsDQvl1ZAAQUUUEABBRpewADY8IfQASiggAIKKKCAAqUJGABL83JpBRRQQAEFFFCg4QUMgA1/CB2AAgoooIACCihQmoABsDQvl1ZAAQUUUEABBRpewADY8IfQASiggAIKKKCAAqUJ/H/A0lsr75/UNgAAAABJRU5ErkJggg==\">"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7f3d90f10668>]"
      ]
     },
     "execution_count": 196,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x = np.linspace(0,3,100)\n",
    "plt.plot(x,x**3-2)\n",
    "plt.plot(x, x-x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 197,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def bisect(f, a, b, tol=10e-5):\n",
    "    \"\"\"\n",
    "    Implements the bisection root finding algorithm, assuming that f is a\n",
    "    real-valued function on [a, b] satisfying f(a) < 0 < f(b).\n",
    "    \"\"\"\n",
    "    lower, upper = a, b\n",
    "\n",
    "    while upper - lower > tol:\n",
    "        middle = 0.5 * (upper + lower)\n",
    "        # === if root is between lower and middle === #\n",
    "        if f(middle) > 0:  \n",
    "            lower, upper = lower, middle\n",
    "        # === if root is between middle and upper  === #\n",
    "        else:              \n",
    "            lower, upper = middle, upper\n",
    "\n",
    "    return 0.5 * (upper + lower)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In fact SciPy provides it’s own bisection function,"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 198,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAALQAAAASCAYAAADyiPTBAAAABHNCSVQICAgIfAhkiAAABfhJREFU\naIHt2nuoZ1UVB/DPdaZ5OM6Yj/ERE+qMkwlWjmkPKSmIoAcSQ5HRJP5nEkGFkkqZaWMToT0EI4s6\nOYMpEzWlA2mkYWShYYUFmfQYm572srKUcG5/rH26+57f/v1++5w7/hO/L/zgnrXX+a6111n7tfZl\nhhn+jzA3Rn4ovofnL4H7RbgUq7EB9+EK/KajtwNfxU/xJE7CudiJ/ZneRnwAT+AprMF78fsOX61e\nrX8tpsXkFHwY+zCPo3EJ/jBGv8XFeEZ6d0g/NuCyZHNl8nMHHhzI9xJciH+L2KzG9gJfrd0W0+JX\n618t3/9wFu5Pjg7FGbgTz0zPh+Ee/BEndnTnO7+n8P6Ozkn4E7Zmsm34EZYP0OvjH9NjcrgYfNsy\n2eX4MVaMeQdOwOO4siOv7cdR2IPjO+8+hM0D+Lbga1iVyW7A33H6ALstpsWv1r9aPnAq9qIRmb+U\nhN6LkzuyLYnzlo58H24UAboWLyjw7cFfcEgmWy1G8wUD9Gr9q43JdjEY8uAfif/gojHvEP2eN5rQ\ntf24BO8s8F6Ojwzg+1jy502Z7PVJ9skBdmvjV+vf4Bxt+igX8E88gmM68r+KkZjjW1O4VojEuK/Q\n9hBu76nX178WjfEx+RluK8gfxF1j3tmK84wmdJ9+fApfKei9B9cN4Dsfj+HVmewtyccdPe120SjH\nr49/U/kOGdU7KPgljhX7oBxPipHXB0eKme+JQttjYi/cR+9g+7dWLLOPFNp+ixcW5IfhtUZXK/r1\n44d4A27GEUm2UizXzQC+m8T26c5MdobYBn6xp91a9PFvMBpLm6HX4LiO7FmJ8+6O/F5xuLgWH8WX\n8ZysfZnYZ5ZG8P7EubyHXl//WjTKMTktybuHOqIv7aEpx3ZsSn93Z+g+/VhpYS/5O5FQu/CagXxd\nbMKjRrdNNXa7aJTjN9S/cXy9jC8F14iRfnZH/nM8O3veJioDecJ9TmwF8qrMcTiQ/Dymp14f/1o0\nyjE5O8k/WGjbmdqOzWSniwHcorSH7tOPtfi6hUP17R17ffmI2fd6sZW6TLkaVmM3R2N8Tg35bpP4\nlqZcgY34Bz5UaFtWeP6XxYeQ9SK4F6bn5YnrgeTnUT31+vjXolGOyYuVk5LYUsyL2Z/Y4n3B4spH\n6d0+/bgUn8Hr8AsLs+bzBvLlWCFWrO8UdGrs5miMz6kh/k3iW5ryFKwUJ9Jxh4US9qVfjiNEnfI6\n0dkT8H1RL50boNfXv0Y5JhuNT+jbUtva9HwRXtnRGfduTT/eLU79LdaISsUBUfLqy1fCK5KPuzNZ\nH7stGpNzqq9/0/iGK0/BTlw1pu0efLsg3698SOji1/jmEvUm+ZejUY7JGvEhP1Fou0tUTogl9PqC\nzriELiHvx1ziPq2g9/bEe2oPPniu0YuKdYnrgDjMDrXb6J9Tk75bke/pqnK0uFIsJVdksvOzv7dY\nmL1yHG10hu5ivbit2r0EvWn+1eBx/MDic0CLk0VFAF4lEmZP9mtnufPS89YuQYZuP9aLi6GHC7qf\nFpchqwpt4/jWpX48YOHASpwriERedhDs1qL2+1ahMXk0bTa9vLVNeea7Mfv7SyJ5c7QXHPlt4bvE\nxcWGTHaxmMlXDdCr9S9HY3xMrhJX5vnSuCnpv2PMO8StZGmGrunHnLgWflmBd23SbasrNXyrRD34\nYVFKa3Fm8vH+AXZzNMbHr893q+EbQXuYObTQdo5Yfu6Y8P45otyzq/O7xeJ65lliyW8DMCdqofda\nHJT34VcWKh9bRCWkW5Go1av1L8ekmByPv+Ftmezj+InJV9+bE+f2gf14c7JxYiZbJ2a1Nw7gu0Yk\nVz4wd4qLqDMH2M0xKX61/tXyIUojd4j/P2hLMY+KU+5bM71TkrEbJhj7c8bR/V3d0X25CMRNYtm9\n2uioXC1q1J8Xs/o38NKC3Vq9Wv9qY0KU4/aKw9FnRQ26tA0hPv7d4uJlXpwXvivKZX36QQzOPbhV\nfOTd4iCXow/fBYlnV9K71eJ7gT52a+NX61+f7zHDDDPMMMMMM8zwdOG/b1tkJJL6cWsAAAAASUVO\nRK5CYII=\n",
      "text/latex": [
       "$$1.2599210498938191$$"
      ],
      "text/plain": [
       "1.2599210498938191"
      ]
     },
     "execution_count": 198,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from scipy.optimize import bisect\n",
    "f =lambda x: x**3-2\n",
    "#f = lambda x: np.sin(4 * (x - 0.25)) + x + x**20 - 1\n",
    "bisect(f, 0, 2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3.2 Function Iteration\n",
    "Function iteration is a relatively simple technique that may be used to compute a\n",
    "fixed-point, $x = g(x)$, of a function from $R^n$ to $R^n$. The technique is also applicable\n",
    "to a rootfinding problem $f(x) = 0$, by recasting it as the equivalent fixed-point\n",
    "problem $x = x - f(x)$.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "https://en.wikipedia.org/wiki/Fixed-point_iteration\n",
    "\n",
    "\n",
    "![](https://upload.wikimedia.org/wikipedia/commons/thumb/8/8b/Sine_fixed_point.svg/375px-Sine_fixed_point.svg.png)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 199,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def myfixpoint(f, x0, maxit = 1000,tol=10e-5):\n",
    "    # \n",
    "    x = x0\n",
    "    for it in range(maxit):\n",
    "        xold = x\n",
    "        x = f(x) \n",
    "        diff = np.linalg.norm(x - xold)\n",
    "\n",
    "        if diff < tol:\n",
    "            \n",
    "            return x\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 200,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def g(x):\n",
    "    return x**0.5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 201,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "g = lambda x: x**0.5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 202,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAALUAAAASCAYAAAAdSp//AAAABHNCSVQICAgIfAhkiAAABl1JREFU\naIHt2nmMX1UVB/DPlK50QYcpTJGG1lorKkoV0BRsNOJutCEuNKlKXFJUUFQaKW7FiBYTFTdwixmp\nC3GtikZN3KMQUFFxo+AC1I02VdBCW2jxj3PfzP29ee/93u83CTMx800mM/e87zvnnnOXd869wzSm\n8X+GgVJ7Fd6FW3AfhrAJ/+xD94PxNuzDQczHG/GPKcIr43zMEv73w2sTu2fgI/gV7sZ+HMqeX4vL\n0t9b8VX8IfGW4znYhp0l2w/Flkzn3an9n4r+H45r8Kga/9raPRabk69zkt6tuKGk7/HYmPo0L/1c\nXMFr60cv+hyROr0hk12I32B21QsNWI7dOCOTbRCDOXMK8Mo4DntFAJtQx2sbuzeISVD1cwhPy7jl\n5wfxloo+PRy3Yk1qD+PPNb6cjOuSvjq0sXsktmNJJluOG7Eyk63G1zA3k12GO3FiH370og8x22/X\nOfiDuAevrHqhAduxBzMy2Tyxe541BXhlfEwM4JYGThOvbewux1Kx0+d9PBXvL+m8Jdnbjvfg0RX9\nmSkm0nmZbCl24bWZ7Hh8AyNil26a1G3sbsK5FfILcUnWfl+y9fxM9uwk+0AffrTVN4od+HqF/AZ8\nr+qFGswWg3ltxbMbcdUk88o4A2fqPqmbeG1j96EKzgJ8Uyy+HD9o6EuBl+IAHtCCW2BE86RuY/dy\nfKVC/nq8N2u/GHfgqZlsfbK/NZO19aOVvmK3WCg+G7dWKPobHtvFWI5BsfL2VTy7A6dMMi/HAjwT\nV1Y8a8vrJXbnVHDejbeKHLFXrMfN+Hcf704Ev8Q6fBYPTLI5Iv0ayXhXiNTsO5nsMSKl+Vwma+tH\nK33F5/K49PvOCkV7sSh1en8Xo8Qn4y6deU+BY7A42Z0s3r2ZfLPuhWE33kRid6rYWH5W8Wx2sjsk\n8u0VuEB8FYgi/zTxZVordq8FWIaLcH2DP03oZpeYuC8Xk/FJIh15Ot6EXzfoXiHSwHNFnTNRP6r0\njWKN2MIvqnhxW3p2dIPyMj4pCrb8dGVYBOk+HDXJPKKw2Jy169KPbryJxO7nYvCq8EeRVxbYIE5S\nhlN7KOn+Pc7OeE8UpwWPqNE7ojn96Ga3wEJ8y1hBeZV6P9fhg2JhbNY5Pv340aRvFI9TP6hXpmfH\n1HS4CouTwY2pPRPvwC+SriMnmTcDn9J5MlHlfxtev7F7Mn5XIS9wWEX7LmMF0dFJ9z7j8/Gd6muI\nEc2TupvdAhfg43gW/pR0/h0nNOieje/jJ8bGol8/6vSN5tS7Gl6cn35XnXnWYZcY7GFROGwRATgk\nOr9nknkbxeAe6OJHG16/sXuVOImow8GK9i48N7ULX242Ph/fiaeItKdXdLMLr8MT8ApxqnICLhUT\n9NMNug+IL9oacV7PxPyo0jeK+WLgy8dKRPX+r4aO9oLb8N1J5g2LT1cZ5d22La+f2M3Cf8VRYBV+\nhB9XyHfqLIRvT9wyfpj6uaTi2Yj6nbqN3QHh0yMreGcn3cen9sOMv+RZZOxcfkGStfWjlb6iUNwr\nEvKlxuMhotqdKBaLW6i6gby/eKeL4GzPOLPS7zNFDn2FuCFrw/uy3mN3ilgMu2v6vFrktmUMiXPk\nAtfovOwoUBSmTV+Rfu0uFkdvN1XwPirOqeeKyXa9iNmqTG/xJRgwluq08aMXfaN4O/6qM+leIVbB\nq0vclcbnPznOE6vv2Ex2vljxc6cAr4xl2l2+1PF6iR28JD2ru9T6ophIOVand/LbvfUi383HYkAc\njX2+RveI+p26jd0B8S8Hp1W8v1DEeo6I9z1i8g9mnJOSvut69KMXfaNYkpS8KJNdit/qLJTWiq3+\n21VKEt6MvxirmFeLCnrNFOGVsVIEptuuX8drG7sCm5Kel9XYOVmcnBS55ID4KvxUZ345A1eLhVvg\nBWJnW1ajuyheD5+A3RcK33Ibi/AFPC+TvVNsNPli3yZSr5P68KOVvvJRyIliwHaIVTcorilvyzir\nRP7zJVHsVGGe2L2Gkp4jxAXD1VOEV2CR+OedVWJi7hefuEt0ph1teG1iV2Cd2DFPV31GTRRirxHF\n0yJxO3mx8ZdLg6Ioni92sgGxuPM04igx+A8ydkS2W/xvyifwmT7srhU3iPvFIjkMHzb+RvIscYZ9\nrygk94hdf0eJ18aPXvRNYxrTmMY0pjGNadxf+B89BH6Fs8LFAQAAAABJRU5ErkJggg==\n",
      "text/latex": [
       "$$0.9999440756158363$$"
      ],
      "text/plain": [
       "0.9999440756158363"
      ]
     },
     "execution_count": 202,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "myfixpoint(g, x0=0.4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "SciPy has a function for finding (scalar) fixed points too"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 203,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array(1.0)"
      ]
     },
     "execution_count": 203,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from scipy.optimize import fixed_point\n",
    "\n",
    "fixed_point(lambda x: x**0.5, 0.4)  # 0.4 is an initial guess"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3.3 Newton's Method\n",
    "\n",
    "\n",
    "In practice, most nonlinear rootfinding problems are solved using *Newton's method*\n",
    "or one of its variants. Newton's method is based on the principle of *successive linearization*. Successive linearization calls for a hard nonlinear problem to be replaced\n",
    "with a sequence of simpler linear problems whose solutions converge to the solution\n",
    "of the nonlinear problem. Newton's method is typically formulated as a rootfinding\n",
    "technique, but may be used to solve a fixed-point problem $x = g(x)$ by recasting it\n",
    "as the rootfinding problem $f(x) = x - g(x) = 0$.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$x:f(x)=0\\,.$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The method starts with a function f defined over the real numbers x, the function's derivative f ′, and an initial guess x0 for a root of the function f. If the function satisfies the assumptions made in the derivation of the formula and the initial guess is close, then a better approximation x1 is\n",
    "\n",
    "$$ x_{1}=x_{0}-{\\frac {f(x_{0})}{f'(x_{0})}}\\,.$$\n",
    "\n",
    "\n",
    "Geometrically, $(x1, 0)$ is the intersection of the x-axis and the tangent of the graph of f at $(x0, f (x0))$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The process is repeated as\n",
    "\n",
    "$$ x_{n+1}=x_{n}-{\\frac {f(x_{n})}{f'(x_{n})}}\\,$$\n",
    "until a sufficiently accurate value is reached."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Derivation\n",
    "\n",
    "Suppose $f : [a, b] → ℝ$ is a differentiable function defined on the interval $[a, b]$ with values in the real numbers $ℝ$. The formula for converging on the root can be easily derived. Suppose we have some current approximation xn. Then we can derive the formula for a better approximation, $x_{n + 1}$ by referring to the diagram on the right. The equation of the tangent line to the curve $y = f (x)$ at the point $x = x_n$ is\n",
    "\n",
    "\n",
    "\n",
    "$$ y=f'(x_{n})\\,(x-x_{n})+f(x_{n}),$$\n",
    "\n",
    "\n",
    "where $f′$ denotes the derivative of the function $f$.\n",
    "\n",
    "\n",
    "The x-intercept of this line (the value of $x$ such that $y = 0$) is then used as the next approximation to the root, $x_{n+1}$. In other words, setting $y$ to zero and $x$ to $x_{n+1}$ gives\n",
    "\n",
    "$$ 0=f'(x_{n})\\,(x_{n+1}-x_{n})+f(x_{n}).$$\n",
    "Solving for $x_{n+1}$ gives\n",
    "\n",
    "$$ {\\displaystyle x_{n+1}=x_{n}-{\\frac {f(x_{n})}{f'(x_{n})}}.}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "https://en.wikipedia.org/wiki/Newton%27s_method\n",
    "\n",
    "![](https://upload.wikimedia.org/wikipedia/commons/thumb/e/e0/NewtonIteration_Ani.gif/450px-NewtonIteration_Ani.gif)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In SciPy this algorithm is implemented by scipy.optimize.newton\n",
    "\n",
    "Unlike bisection, the Newton-Raphson method uses local slope information\n",
    "\n",
    "This is a double-edged sword:\n",
    "\n",
    "- When the function is well-behaved, the Newton-Raphson method is faster than bisection\n",
    "- When the function is less well-behaved, the Newton-Raphson might fail\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 204,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAALUAAAASCAYAAAAdSp//AAAABHNCSVQICAgIfAhkiAAABmFJREFU\naIHt2WusXUUVB/BfS+mD3hYBKUhKkFZEElSKb6pEY6LxEVObmGhskA8mQBQVhQgkagWrRUMNYloF\nH0dKFIKPKtQojWgwPlLEaCAGND5aKyr4qoi0Cr1+WLPvnbvv7N7Z51Y/mPNPds7Za9asx8zaM7PW\nMMII/2eY00E/Aj/EM2Yh+7m4FIuwHDvxXvyuxbcRX8V92I+T8RpsxZ6MbwXeh314HIvxbvyhJa+W\nr9a+BjONyan4EHZhHE/EJfhjB3+Di3F46juMH8txWdK5INm5EfcMKa/Wvhq9r8An8FM8Kub3QNa+\nE5vT/+fjvMS3KD0bCn70nTfwHNyVjB0WZ+J2PCG9j+FOPIgnt3jHW8/jeE+L52T8CWsz2joxWPOG\n4OtjHzOPyZHiA1yX0S7HvZjf0QdOwiNY36LX+nEMtuFJrb7345Qh5NXaV6v3XabPb/McwMsT3yp8\nDQuzvpvxd5yR0frOm9OwHQOxIs0mqLfjKS3aqiTzphZ9F64Tg3Q1nlmQtw1/wdyMtkisOucOwVdr\nX+2YbBADmwfI0fg3LujoQ/g9bnrQ1PpxCS4syL0cVw0hr9a+Wr1bcKJY6XPdq3FN9v7RpOd1Ge3V\nifaxjNYnrqZhYHZB/Q/sxrIW/a9ixcjxnRlkzRfBsbPQdj9u68nX174GA91j8nPcWqDfgzs6+qzF\n600Pmj5+bMFXCnzvxKYh5NXYV6sXPl7gGcPXxUfV4Bzsxcsy2huS7o0ZrWre5vrv4Nc4Tpzbcuw3\n1ZkaHC1WwH2Ftr3ijNWH71Dbt0RsubsLbQ/gWQX6GF6pvLr08eMnWIPP46hEWyCOFoMh5NXYV6sX\n3lro+2FxBn40o90gjnC3Z7QzxVH0CxltVvM2MLuVejGOb9FOSDK/3aJ/XyQcV+Mj+DKemrUfJs51\npZVmT5I5rwdfX/saDJTH5PREbydShC9NIpVjA1am/+2VsI8fC0ye9X8vgupGkaANI6/Gvlq9JawW\nieNMWImHTD+6DTNvExiYXVCX8EHx5Z3Vov9SnLsarBMVg9z4z4jtJa/WHC+SjXGT21EtXx/7GgyU\nx+SsRH9/oW1rajsuo50hPuIGpaDp48cSfMNkAnZbS19feTX21ept424dCV3CGlwrjnOX6a7O5Zhp\n3iYwcGiDegUexgcKbYcV3v9paoJwrHD0vPQ+L8n6cbLzmJ58fexrMFAek+fpnvibUtsJ6X0uPmdq\nRaTUt48fl+J6vAq/Mrl6Pn0IebX21erN8VL8rKOtjfli5f2e7jmjbt4mMHDognqBqBxsmokxw670\n5DhK1Fk3CSdOwo/E2WzOEHx97Rsoj8kK3RN/a2pbkt4vwEtaPF19a/y4SFQEGiwWlYQDolzXV16t\nfX30NviS2DFq8eKk+5aO9t5xNXDognorruhouxPfLdD3KCc2bfwW35ol38HsyzFQHpPFYjKvKbTd\nITJzYru/tsDTFdQl5H7MSbJPL/Cdn+Se1kNerX3D6D1cVC42dNjxNNMvtZaarGePFfrUztsEBg5N\nUK83/SLlnOz/wyKTbmOfKDcdDMcKG8+fBd9M9uUY6B6Tu0VS2MZukwnMOuwQNePm2Z5k3pfe1xZk\nNGj7sUw5CSUCb6+o4dbKq7VvGL2rU5+LCn2Wit3iMZPJKbFYNOf1I1t91quftwkMHDyoTzFzCWWd\n8pd0Xfb/i+I6OUdTTM+Nfoe43Fie0S4WK/rCIfhq7csx0D0mV4hr2vx4szLxv6WjD5E0lVbqGj/m\niCvuFxbkLkm8TeD1GZeZ7Oujt8GbkpzSRdRCUUf/hSg/Nnh26nNXi3/Geeu6Im0cPUIkbTnOFhcm\nO0xec7ZxtjhjfVOUenJ9eWBclfjeLGqNc8TX/ANRz2wwlux4LL2vErdarzX1mFLLV2tfjoONyRa8\nTQz41kS7UCRG13fII7bl/LePH+N4Oz4pkrXfJPpSfFoE8v4e8mrt66O3QVNd+VdBxz5Ryn3Q5FFN\n0vGIqR9C73lblpjvNbnsPyS2zzdmfKeKktvmtoAMf85ktJ8rW7wvEsnADWKLu9L01WORcPyzYnXf\ngRcU9Nby1dpXOyZEKWy7GPRPiePIicpYmmQ8kGTuEx/ymp5+EBO9DTeLasstIsnK0UdejX21ehus\nwd/E6tuFc5OcG5N9N5t6X0G/uBphhBFGGGGEEUb4X+E/2El+R8U7ZLYAAAAASUVORK5CYII=\n",
      "text/latex": [
       "$$1.2599210498948732$$"
      ],
      "text/plain": [
       "1.2599210498948732"
      ]
     },
     "execution_count": 204,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from scipy.optimize import newton\n",
    "\n",
    "f = lambda x: x**3 -2\n",
    "\n",
    "newton(f, 0.2)   # Start the search at initial condition x = 0.2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 205,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "newton?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Newton \n",
    "\n",
    "Signature: newton(func, x0, fprime=None, args=(), tol=1.48e-08, maxiter=50, fprime2=None)\n",
    "Docstring:\n",
    "Find a zero using the Newton-Raphson or secant method.\n",
    "\n",
    "Find a zero of the function `func` given a nearby starting point `x0`.\n",
    "The Newton-Raphson method is used if the derivative `fprime` of `func`\n",
    "is provided, otherwise the secant method is used.  If the second order\n",
    "derivate `fprime2` of `func` is provided, parabolic Halley's method\n",
    "is used.\n",
    "\n",
    "### Parameters\n",
    "\n",
    "\n",
    "func : function\n",
    "\n",
    "    The function whose zero is wanted. It must be a function of a\n",
    "    **single variable\n",
    "    ** of the form f(x,a,b,c...), where a,b,c... are extra\n",
    "    arguments that can be passed in the `args` parameter.\n",
    "    \n",
    "x0 : float\n",
    "\n",
    "    An initial estimate of the zero that should be somewhere near the\n",
    "    actual zero.\n",
    "    \n",
    "fprime : function, optional\n",
    "\n",
    "    The derivative of the function when available and convenient. If it\n",
    "    is None (default), then the secant method is used.\n",
    "\n",
    "\n",
    "\n",
    "Since `newton` only take one variable function. We use `fsolve` in scipy.optimize"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example: Cournot Duopoly\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To illustrate the use of this function, consider a simple Cournot duopoly model, in which the inverse demand for a good is\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 206,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "q = np.array([1, 1])\n",
    "c = np.array([0.6, 0.8])\n",
    "eta = 1.6\n",
    "e = -1/eta    \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 207,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.68765683,  0.55048639])"
      ]
     },
     "execution_count": 207,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fval = np.sum(q)**e + e * np.sum(q)**(e-1)*(q) - np.diag(c).dot(q)\n",
    "fjac = e*np.sum(q)**(e-1)*(np.ones([2,2]))+e*np.sum(q)**(e-1)*(np.eye(2))+(e-1)*e*np.sum(q)**(e-2)*(q).dot(np.array([1, 1]))-np.diag(c)\n",
    "q = q - np.linalg.inv(fjac).dot(fval)\n",
    "q"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 208,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-0.1542114, -0.3542114])"
      ]
     },
     "execution_count": 208,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fval"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 209,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-0.67598669,  0.12664449],\n",
       "       [ 0.12664449, -0.87598669]])"
      ]
     },
     "execution_count": 209,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fjac"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Newton method"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 210,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.83961268,  0.68883816])"
      ]
     },
     "execution_count": 210,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "maxit = 1000\n",
    "tol = 1/10000\n",
    "# start from [0.2 0.2] will break \n",
    "#q = np.array([1, 1])\n",
    "#c = np.array([0.6, 0.8])\n",
    "#eta = 1.6\n",
    "#e = -1/eta\n",
    "\n",
    "for it in np.arange(maxit):\n",
    "\n",
    "    fval = np.sum(q)**e + e * np.sum(q)**(e-1)*(q) - np.diag(c).dot(q)\n",
    "    fjac = e*np.sum(q)**(e-1)*(np.ones([2,2]))+e*np.sum(q)**(e-1)*(np.eye(2))+(e-1)*e*np.sum(q)**(e-2)*(q).dot(np.array([1, 1]))-np.diag(c)\n",
    "    q = q - np.linalg.inv(fjac).dot(fval)\n",
    "    if np.linalg.norm(fval)<tol:\n",
    "        break\n",
    "        \n",
    "q"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 211,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([  5.34511334e-05,   6.51508959e-05])"
      ]
     },
     "execution_count": 211,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fval"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 212,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def cournot(q):\n",
    "    c = np.array([0.6, 0.8])\n",
    "    eta = 1.6\n",
    "    e = -1 / eta\n",
    "    fval = np.sum(q)**e + e * np.sum(q)**(e-1)*(q) - np.diag(c).dot(q)\n",
    "    fjac = e*np.sum(q)**(e-1)*(np.ones([2,2]))+e*np.sum(q)**(e-1)*(np.eye(2))+\\\n",
    "            (e-1)*e*np.sum(q)**(e-2)*(q).dot(np.array([1, 1]))-np.diag(c)\n",
    "    return fval, fjac"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 213,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def mynewton(f, x0, maxit=1000, tol=1/1000 ):\n",
    "    x = x0\n",
    "    for it in np.arange(maxit):\n",
    "        fval, fjac = f(x)\n",
    "        x = x - np.linalg.inv(fjac).dot(fval)\n",
    "        if np.linalg.norm(fval)<tol:\n",
    "            break\n",
    "    return x    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 214,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.84001971,  0.689215  ])"
      ]
     },
     "execution_count": 214,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mynewton(cournot, x0 = np.array([1, 1])) # converge "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### without backstepping, it does not converge with x0= np.array([0.2, 0.2]). It is overstepping too far.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 215,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ nan,  nan])"
      ]
     },
     "execution_count": 215,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mynewton(cournot, x0 = np.array([0.2, 0.2])) # break ? without "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "More robust implementation can be found in\n",
    "https://github.com/randall-romero/CompEcon-python/blob/master/compecon/nonlinear.py"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### fsolve function\n",
    "\n",
    "Since scipy.optimize.newton only take one variable function. We use fsolve in scipy.optimize"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 216,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from scipy.optimize import fsolve"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 217,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "\n",
    "def cournot(q):\n",
    "    c = np.array([0.6, 0.8])\n",
    "    eta = 1.6\n",
    "    e = -1/eta\n",
    "    fval = np.sum(q)**e + e * np.sum(q)**(e-1)*(q) - np.diag(c).dot(q)\n",
    "    fjac = e*np.sum(q)**(e-1)*(np.ones([2,2]))+e*np.sum(q)**(e-1)*(np.eye(2))+(e-1)*e*np.sum(q)**(e-2)*(q).dot(np.array([1, 1]))-np.diag(c)\n",
    "    return fval\n",
    "\n",
    "def cournotjac(q):\n",
    "    c = np.array([0.6, 0.8])\n",
    "    eta = 1.6\n",
    "    e = -1/eta\n",
    "    fval = np.sum(q)**e + e * np.sum(q)**(e-1)*(q) - np.diag(c).dot(q)\n",
    "    fjac = e*np.sum(q)**(e-1)*(np.ones([2,2]))+e*np.sum(q)**(e-1)*(np.eye(2))+(e-1)*e*np.sum(q)**(e-2)*(q).dot(np.array([1, 1]))-np.diag(c)\n",
    "    return fjac    \n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 218,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "q = fsolve(func = cournot, x0= np.array([0.2, 0.2]))  # without jacobian"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 219,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.8395676 ,  0.68879643])"
      ]
     },
     "execution_count": 219,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "q"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 220,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "q = fsolve(func = cournot,x0= np.array([0.2, 0.2]), fprime= cournotjac)  # with jacobian"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 221,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.8395676 ,  0.68879643])"
      ]
     },
     "execution_count": 221,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "q"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3.4 Quasi-Newton Methods\n",
    "\n",
    "Quasi-Newton methods offer an alternative to Newton's method for solving rootfinding\n",
    "problems. Quasi-Newton methods are based on the same successive linearization\n",
    "principle as Newton's method, except that they replace the Jacobian $f'$ with an\n",
    "estimate that is easier to compute. Quasi-Newton methods are easier to implement\n",
    "and less likely to fail due to programming errors than Newton's method because the\n",
    "analyst need not explicitly code the derivative expressions. Quasi-Newton methods,\n",
    "however, often converge more slowly than Newton's method and additionally require\n",
    "the analyst to supply an initial estimate of the function's Jacobian.\n",
    "\n",
    "\n",
    "\n",
    "The **secant method** is the most widely used univariate quasi-Newton method. The\n",
    "secant method is identical to the univariate Newton method, except that it replaces\n",
    "the derivative of f with a finite-difference approximation constructed from the function\n",
    "values at the two previous iterates:\n",
    "\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "https://en.wikipedia.org/wiki/Secant_method\n",
    "\n",
    "![](https://upload.wikimedia.org/wikipedia/commons/thumb/9/92/Secant_method.svg/450px-Secant_method.svg.png)\n",
    "\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Newton's method to find zeroes of a function ${\\displaystyle g}$ of multiple variables is given by: $ x_{n+1}=x_{n}-[J_{g}(x_{n})]^{-1}g(x_{n})\\,\\!$ where ${\\displaystyle [J_{g}(x_{n})]^{-1}}$ is the left inverse of the Jacobian matrix $ J_{g}(x_{n})$ of ${\\displaystyle g}$ evaluated for ${\\displaystyle x_{n}}$.\n",
    "\n",
    "Strictly speaking, any method that replaces the exact Jacobian ${\\displaystyle J_{g}(x_{n})}$ with an approximation is a **quasi-Newton method.** "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$f'(x_{n-1}) ~= {\\frac {f(x_{n-1})-f(x_{n-2})} {x_{n-1}-x_{n-2}}   }$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "\n",
    "$${\\displaystyle x_{n}=x_{n-1}-f(x_{n-1}){\\frac {x_{n-1}-x_{n-2}}{f(x_{n-1})-f(x_{n-2})}}={\\frac {x_{n-2}f(x_{n-1})-x_{n-1}f(x_{n-2})}{f(x_{n-1})-f(x_{n-2})}}}  $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "#### *Broyden's method* \n",
    "\n",
    "is the most popular multivariate generalization of the univariate\n",
    "secant method. Broyden's method generates a sequence of vectors $x^{(k)}$ and matrices\n",
    "$A^{(k)}$ that approximate the root of $f$ and the Jacobian $f'$ at the root, respectively.\n",
    "\n",
    "\n",
    "\n",
    "Broyden's method begins with the analyst supplying a guess $x^{(0)}$ for the root of the\n",
    "function and a guess $A^{(0)}$ for the Jacobian of the function at the root. Often, $A^{(0)}$\n",
    "is set equal to the numerical Jacobian of f at $x^{(0)}$. \n",
    "\n",
    "\n",
    "Alternatively, some analysts use\n",
    "a rescaled identity matrix for $A^{(0)}$, though this typically will require more iterations\n",
    "to obtain a solution than if a numerical Jacobian is computed at the outset. Given\n",
    "$x^{(k)}$ and $A^{(k)}$, one updates the root approximation by solving the linear rootfinding\n",
    "problem obtained by replacing f with its first-order Taylor approximation about $x^{(k)}$:\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "https://en.wikipedia.org/wiki/Broyden%27s_method\n",
    "\n",
    "\n",
    "\n",
    "To solve a system of $k$ nonlinear equations\n",
    "\n",
    "$${\\displaystyle \\mathbf {f} (\\mathbf {x} )=\\mathbf {0} ,}$$ \n",
    "\n",
    "\n",
    "where $f$ is a vector-valued function of vector $x$:\n",
    "\n",
    "\n",
    "$${\\displaystyle \\mathbf {x} =(x_{1},x_{2},x_{3},\\dotsc ,x_{k})}$$\n",
    "\n",
    "$${\\displaystyle \\mathbf {f} (\\mathbf {x} )=(f_{1}(x_{1},x_{2},\\dotsc ,x_{k}),f_{2}(x_{1},x_{2},\\dotsc ,x_{k}),\\dotsc ,f_{k}(x_{1},x_{2},\\dotsc ,x_{k}))} $$\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "For such problems, Broyden gives a generalization of the one-dimensional Newton's method, replacing the derivative with the Jacobian J. The Jacobian matrix is determined iteratively based on the secant equation in the finite difference approximation:\n",
    "\n",
    "$${\\displaystyle \\mathbf {J} _{n}(\\mathbf {x} _{n}-\\mathbf {x} _{n-1})\\simeq \\mathbf {f} (\\mathbf {x} _{n})-\\mathbf {f} (\\mathbf {x} _{n-1}),} $$\n",
    "\n",
    "where n is the iteration index. For clarity, let us define:\n",
    "\n",
    "$${\\displaystyle \\mathbf {f} _{n}=\\mathbf {f} (\\mathbf {x} _{n}),} $$\n",
    "$${\\displaystyle \\Delta \\mathbf {x} _{n}=\\mathbf {x} _{n}-\\mathbf {x} _{n-1},} $$\n",
    "$${\\displaystyle \\Delta \\mathbf {f} _{n}=\\mathbf {f} _{n}-\\mathbf {f} _{n-1},}$$\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "so the above may be rewritten as:\n",
    "\n",
    "$${\\displaystyle \\mathbf {J} _{n}\\Delta \\mathbf {x} _{n}\\simeq \\Delta \\mathbf {f} _{n}.} $$\n",
    "\n",
    "\n",
    "The above equation is underdetermined when $k$ is greater than one. \n",
    "\n",
    "Broyden suggests using the current estimate of the Jacobian matrix $\\mathbf {J} _{n-1}$ and improving upon it by taking the solution to the secant equation that is a minimal modification to $\\mathbf {J} _{n-1}$:\n",
    "\n",
    "$${\\displaystyle \\mathbf {J} _{n}=\\mathbf {J} _{n-1}+{\\frac {\\Delta \\mathbf {f} _{n}-\\mathbf {J} _{n-1}\\Delta \\mathbf {x} _{n}}{\\|\\Delta \\mathbf {x} _{n}\\|^{2}}}\\Delta \\mathbf {x} _{n}^{\\mathrm {T} }} $$\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This minimizes the following Frobenius norm:\n",
    "\n",
    "$${\\displaystyle \\|\\mathbf {J} _{n}-\\mathbf {J} _{n-1}\\|_{\\mathrm {f} }.} $$\n",
    "We may then proceed in the Newton direction:\n",
    "\n",
    "$${\\displaystyle \\mathbf {x} _{n+1}=\\mathbf {x} _{n}-\\mathbf {J} _{n}^{-1}\\mathbf {f} (\\mathbf {x} _{n}).}$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Broyden also suggested using the Sherman-Morrison formula to update directly the inverse of the Jacobian matrix:\n",
    "\n",
    "\n",
    "$${\\displaystyle \\mathbf {J} _{n}^{-1}=\\mathbf {J} _{n-1}^{-1}+{\\frac {\\Delta \\mathbf {x} _{n}-\\mathbf {J} _{n-1}^{-1}\\Delta \\mathbf {f} _{n}}{\\Delta \\mathbf {x} _{n}^{\\mathrm {T} }\\mathbf {J} _{n-1}^{-1}\\Delta \\mathbf {f} _{n}}}\\Delta \\mathbf {x} _{n}^{\\mathrm {T} }\\mathbf {J} _{n-1}^{-1}} $$\n",
    "\n",
    "\n",
    "This first method is commonly known as the **\"good Broyden's method\".**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In theory, Broyden's method converges if $f$ is continuously differentiable, if $x_0$\n",
    "is \"sufficiently\" close to a root of $f$ at which $f_0$ is invertible, and if $J_0^{-1}$ are\n",
    "\"sufficiently\" close to the Jacobian or inverse Jacobian of $f$ at that root.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Like Newton's method, the robustness of Broyden's method depends on the regularity of\n",
    "f and its derivatives. Broyden's method may also have diÆculty computing a precise\n",
    "root estimate if f0 is ill-conditioned near the root."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 222,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#https://github.com/randall-romero/CompEcon-python/blob/master/compecon/tools.py\n",
    "# The script also computes an initial guess for the inverse Jacobian by inverting the finite difierence derivative computed\n",
    "# using the toolbox function fdjac, which is discussed in Chapter 5 (page 107).\n",
    "def jacobian(func, x, *args, **kwargs):\n",
    "\n",
    "    # if type(func(x, *args, **kwargs)) is tuple:\n",
    "    #     F = lambda x: func(x, *args, **kwargs)[0]\n",
    "    # else:\n",
    "    #     F = lambda x: func(x, *args, **kwargs)\n",
    "    F = lambda z: func(z, *args, **kwargs)\n",
    "\n",
    "    x = x.flatten()\n",
    "    dx = x.size\n",
    "    f = F(x)\n",
    "    df = f.size\n",
    "    x = x.astype(float)\n",
    "\n",
    "    ''' Compute Jacobian'''\n",
    "    tol = np.spacing(1) ** (1/3)\n",
    "\n",
    "    h = tol * np.maximum(abs(x), 1)\n",
    "    x_minus_h = x - h\n",
    "    x_plus_h = x + h\n",
    "    deltaX = x_plus_h - x_minus_h\n",
    "    fx = np.zeros((dx, df))\n",
    "\n",
    "    for k in range(dx):\n",
    "        xx = x.copy()\n",
    "        xx[k] = x_plus_h[k]\n",
    "        fplus = F(xx)\n",
    "\n",
    "        xx[k] = x_minus_h[k]\n",
    "        fminus = F(xx)\n",
    "\n",
    "        fx[k] = np.squeeze((fplus - fminus) / deltaX[k])  # fixme doing this to deal with broadcasting\n",
    "\n",
    "    return fx.T\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 229,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def cournot(q):\n",
    "    c = np.array([0.6, 0.8])\n",
    "    eta = 1.6\n",
    "    e = -1 / eta\n",
    "    fval = np.sum(q)**e + e * np.sum(q)**(e-1)*(q) - np.diag(c).dot(q)\n",
    "    fjac = e*np.sum(q)**(e-1)*(np.ones([2,2]))+e*np.sum(q)**(e-1)*(np.eye(2))+\\\n",
    "            (e-1)*e*np.sum(q)**(e-2)*(q).dot(np.array([1, 1]))-np.diag(c)\n",
    "    return fval"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 230,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "f = cournot\n",
    "x = np.array([0.2, 0.2])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 231,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 1.09894808,  1.05894808])"
      ]
     },
     "execution_count": 231,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fval = f(x)\n",
    "fval"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 232,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-3.88977465, -0.5194381 ],\n",
       "       [-0.5194381 , -4.08977465]])"
      ]
     },
     "execution_count": 232,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# using function to calculate Jacobian\n",
    "fjac = jacobian(cournot,np.array([0.2, 0.2]) )\n",
    "fjac"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 233,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-0.26151986,  0.03321537],\n",
       "       [ 0.03321537, -0.2487309 ]])"
      ]
     },
     "execution_count": 233,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# using numpy to calculate inversed Jacobian\n",
    "fjacinv = np.linalg.pinv(np.atleast_2d(fjac))\n",
    "fjacinv "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 235,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.83480764,  0.68819305])"
      ]
     },
     "execution_count": 235,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "maxit = 100\n",
    "tol = 1/100\n",
    "\n",
    "for it in range(maxit):\n",
    "    fnorm = np.linalg.norm(fval)\n",
    "    if fnorm<tol: \n",
    "        break \n",
    "    d = -(fjacinv.dot(fval))\n",
    "    x = x+d\n",
    "    fold = fval\n",
    "    fval = f(x)\n",
    "    u = fjacinv.dot((fval-fold))\n",
    "    fjacinv = fjacinv + np.outer((d-u), np.dot(d.T, fjacinv))/np.dot(d.T,u) ## ? np.outer!!!!! Key\n",
    "x   "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In linear algebra, an outer product is the tensor product of two coordinate vectors, a special case of the Kronecker product of matrices. The outer product of two coordinate vectors ${\\displaystyle \\mathbf {u} } $  and ${\\displaystyle \\mathbf {v} } $ , denoted ${\\displaystyle \\mathbf {u} \\otimes \\mathbf {v} } $, is a matrix ${\\displaystyle \\mathbf {w} } $  such that the coordinates satisfy ${\\displaystyle w_{ij}=u_{i}v_{j}} $. The outer product for general tensors is also called the tensor product.\n",
    "\n",
    "The outer product contrasts with the dot product, which takes as input a pair of coordinate vectors and produces a scalar.\n",
    "\n",
    "(source: https://en.wikipedia.org/wiki/Outer_product)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Without analytic Jacobian, we use Broyden method to find the root."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 246,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def mybroyden(f, x0, maxit = 1000, tol = 1/1000):\n",
    "    \n",
    "    x=x0\n",
    "    A = f(x)\n",
    "    _is_there_jacobian = (type(A) is tuple) and (len(A) == 2)\n",
    "    \n",
    "    if _is_there_jacobian:\n",
    "        print('Jacobian was  provided by user!')\n",
    "        fval,fjac = f(x)\n",
    "    else:    \n",
    "        print('Jacobian was not provided by user!')\n",
    "        fval = f(x)\n",
    "        try:\n",
    "            fjac = jacobian(f,x)\n",
    "        except NameError:\n",
    "            print(\"jacobian function Not in scope!\\n Using identity matrix as jacobian matrix\")\n",
    "            fjac = np.identity(x.size)\n",
    "        else:\n",
    "            print(\"jacobian function In scope!\")\n",
    "            #fjac = jacobian(f,x)  \n",
    "    # using numpy to calculate inversed Jacobian, only for initialize. later, we update it without calculation of new inversed Jacobian.\n",
    "    fjacinv = np.linalg.pinv(np.atleast_2d(fjac))\n",
    "    # or we can use identity matrix as an initiative matrix\n",
    "    #fjacinv = - np.identity(x.size)\n",
    "    for it in range(maxit):\n",
    "        fnorm = np.linalg.norm(fval)\n",
    "        if fnorm<tol: \n",
    "            break \n",
    "        d = -(fjacinv.dot(fval))\n",
    "        # update x\n",
    "        x = x+d\n",
    "        fold = fval\n",
    "        fval = f(x)[0] if _is_there_jacobian else f(x)# two outputs\n",
    "        # update jacobian\n",
    "        u = fjacinv.dot((fval-fold))\n",
    "        fjacinv = fjacinv + np.outer((d-u), np.dot(d.T, fjacinv))/np.dot(d.T,u) ## ? np.outer !!! key\n",
    "    return x   \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 247,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Jacobian was not provided by user!\n",
      "jacobian function In scope!\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([ 0.83945538,  0.6889152 ])"
      ]
     },
     "execution_count": 247,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mybroyden(cournot, np.array([0.2, 0.2]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Scipy.optimize.broyden1\n",
    "\n",
    "https://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.optimize.broyden1.html\n",
    "\n",
    "scipy.optimize.broyden1\n",
    "\n",
    "\n",
    "    scipy.optimize.broyden1(F, xin, iter=None, alpha=None, reduction_method='restart', max_rank=None, verbose=False, maxiter=None, f_tol=None, f_rtol=None, x_tol=None, x_rtol=None, tol_norm=None, line_search='armijo', callback=None, **kw)\n",
    "\n",
    "Find a root of a function, using Broyden’s first Jacobian approximation.\n",
    "\n",
    "This method is also known as “Broyden’s good method”.\n",
    "\n",
    "This algorithm implements the inverse Jacobian Quasi-Newton update\n",
    "\n",
    "![](https://docs.scipy.org/doc/scipy-0.14.0/reference/_images/math/55c80f4be075d374bb3a6683ef676f504969d3a9.png)\n",
    "\n",
    "which corresponds to Broyden’s first Jacobian update\n",
    "\n",
    "![](https://docs.scipy.org/doc/scipy-0.14.0/reference/_images/math/4b7a3d3206120c879b3b7b94a95c2262ffb7b62a.png)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.8395676 ,  0.68879643])"
      ]
     },
     "execution_count": 42,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from scipy.optimize import broyden1\n",
    "broyden1(cournot,np.array([0.2, 0.2]) )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3.5 Problems With Newton Methods\n",
    "\n",
    "\n",
    "\n",
    "Several difficulties commonly arise in the application of Newton and quasi-Newton\n",
    "methods to solving multivariate non-linear equations. The most common cause of\n",
    "failure of Newton-type methods is coding errors committed by the analyst. The next\n",
    "most common cause of failure is the specification of a starting point that is not sufficiently\n",
    "close to a root. And yet another common cause of failure is an ill-conditioned\n",
    "Jacobian at the root. These problems can often be mitigated by appropriate action,\n",
    "though they cannot always be eliminated altogether."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "The *first cause of failure*, coding error, may seem obvious and not specific to\n",
    "rootfinding problems. It must be emphasized, however, that with Newton's method,\n",
    "the likelihood of committing an error in coding the analytic Jacobian of the function\n",
    "is often high. A careful analyst can avoid Jacobian coding errors in two ways. First,\n",
    "the analyst could use Broyden's method instead of Newton's method to solve the\n",
    "rootfinding problem. Broyden's method is derivative-free and does not require the\n",
    "explicit coding of the function's analytic Jacobian. Second, the analyst can perform\n",
    "a simple, but highly effiective check of his code by comparing the values computed\n",
    "by his analytic derivatives to those computed using finite diffierence methods. Such a\n",
    "check will almost always detect an error in either the code that returns the function's\n",
    "value or the code that returns its Jacobian."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "The *second problem, a poor starting value*, can be partially addressed by **\"backstepping\"**.\n",
    "If taking a full Newton (or quasi-Newton) step $x+d$ does not offer an improvement\n",
    "over the current iterate $x$, then one \"backsteps\" toward the current iterate\n",
    "$x$ by repeatedly cutting $d$ in half until $x+d$ does offer an improvement. Whether a step\n",
    "$d$ offers an improvement is measured by the Euclidean norm $||f(x)|| = \\frac{1}{2} f(x)^{T}f(x)$.\n",
    "Clearly, $||f(x)||$ is precisely zero at a root of f, and is positive elsewhere. Thus,\n",
    "one may view an iterate as yielding an improvement over the previous iterate if it\n",
    "reduces the function norm, that is, if $||f(x)|| > ||f(x+ d)||$. Backstepping prevents\n",
    "Newton and quasi-Newton methods from taking a large step in the wrong direction,\n",
    "substantially improving their robustness."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "necessarily prevent Newton type methods\n",
    "from getting stuck at a local minimum of $||f(x)||$. If $||f(x)||$ must decrease with\n",
    "each step, it may be difficult to find a step length that moves away from the current\n",
    "value of x. Most good root-finding algorithms employ so mechanism for getting\n",
    "unstuck. We use a very simple one in which the backsteps continue until either\n",
    "$||f(x)|| > ||f(x+ d)||$ or $ ||f(x+ d/2)|| > ||f(x+ d)||$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def cournot(q):\n",
    "    c = np.array([0.6, 0.8])\n",
    "    eta = 1.6\n",
    "    e = -1 / eta\n",
    "    fval = np.sum(q)**e + e * np.sum(q)**(e-1)*(q) - np.diag(c).dot(q)\n",
    "    fjac = e*np.sum(q)**(e-1)*(np.ones([2,2]))+e*np.sum(q)**(e-1)*(np.eye(2))+\\\n",
    "            (e-1)*e*np.sum(q)**(e-2)*(q).dot(np.array([1, 1]))-np.diag(c)\n",
    "    return fval, fjac"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### mynewton_backstep"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.83901993,  0.68828942])"
      ]
     },
     "execution_count": 44,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x0= np.array([0.2, 0.2])\n",
    "f = cournot\n",
    "def mynewton_backstep(f, x0, maxit=1000, tol=1/1000, maxsteps = 1000 ):\n",
    "    x = x0\n",
    "    for it in np.arange(maxit):\n",
    "        fval, fjac = f(x)\n",
    "        fnorm = np.linalg.norm(fval)\n",
    "        if np.linalg.norm(fval)<tol:\n",
    "            break\n",
    "        d =   - np.linalg.inv(fjac).dot(fval)  \n",
    "        fnormold = np.inf # initial fnormold\n",
    "        for backstep in np.arange(maxsteps):\n",
    "            fvalnew = f(x+d)[0] # first output/outcome\n",
    "            fnormnew = np.linalg.norm(fvalnew)\n",
    "            if fnormnew < fnorm: \n",
    "                break\n",
    "            if fnormold < fnormnew: \n",
    "                d=2*d\n",
    "                break\n",
    "            fnormold = fnormnew;\n",
    "            d = d/2;\n",
    "        x = x + d\n",
    "\n",
    "    return x    \n",
    "mynewton_backstep(f,x0)  # even start from [0.2,0.2], still converge with backstepping."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### even start from [0.2,0.2], still converge with backstepping."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Note: Backtracking\n",
    "Newton’s method may not converge for a variety of reasons. One potential problem occurs when the step from $x_k$ to $x_{k+1}$ is so large that the root is stepped over completely. Backtracking is a strategy that combats the problem of overstepping by moving only a fraction of the full step from $x_k$ to $x_{k+1}$.\n",
    "\n",
    "https://github.com/OpenSourceMacro/BootCamp2017/blob/master/Computation/Wk4_DifIntOpt/ACME_Newtons.pdf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Object `scipy.optimize.newton` not found.\n"
     ]
    }
   ],
   "source": [
    "scipy.optimize.newton?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### mybroyden_backstep"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def cournot(q):\n",
    "    c = np.array([0.6, 0.8])\n",
    "    eta = 1.6\n",
    "    e = -1 / eta\n",
    "    fval = np.sum(q)**e + e * np.sum(q)**(e-1)*(q) - np.diag(c).dot(q)\n",
    "    fjac = e*np.sum(q)**(e-1)*(np.ones([2,2]))+e*np.sum(q)**(e-1)*(np.eye(2))+\\\n",
    "            (e-1)*e*np.sum(q)**(e-2)*(q).dot(np.array([1, 1]))-np.diag(c)\n",
    "    return fval # for Broyden, no need for analytic jacobian\n",
    "\n",
    "x0= np.array([0.2, 0.2])\n",
    "f = cournot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.8395465 ,  0.68872313])"
      ]
     },
     "execution_count": 47,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def mybroyden_backstep(f, x0, maxit = 1000, tol = 1/1000, maxsteps = 100):\n",
    "    x=x0\n",
    "    fval = f(x)\n",
    "    # using function to calculate Jacobian\n",
    "    fjac = jacobian(f,x)\n",
    "    # using numpy to calculate inversed Jacobian\n",
    "    fjacinv = np.linalg.pinv(np.atleast_2d(fjac))\n",
    "    for it in range(maxit):\n",
    "        fnorm = np.linalg.norm(fval)\n",
    "        if fnorm<tol: \n",
    "            break \n",
    "        d = -(fjacinv.dot(fval))\n",
    "        fnormold = np.inf # initial fnormold\n",
    "        for backstep in np.arange(maxsteps):\n",
    "            fvalnew = f(x+d) # first output/outcome\n",
    "            fnormnew = np.linalg.norm(fvalnew)\n",
    "            if fnormnew < fnorm: \n",
    "                break\n",
    "            if fnormold < fnormnew: \n",
    "                d=2*d\n",
    "                break\n",
    "            fnormold = fnormnew;\n",
    "            d = d/2;           \n",
    "        # update x\n",
    "        x = x+d\n",
    "        fold = fval\n",
    "        fval = f(x)\n",
    "        # update jacobian\n",
    "        u = fjacinv.dot((fval-fold))\n",
    "        fjacinv = fjacinv + ((d-u).dot(d.T)*(fjacinv))/(d.T.dot(u))\n",
    "    return x  \n",
    "\n",
    "mybroyden_backstep(f,x0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The **third problem, an ill-conditioned Jacobian at the root**, occurs less often, but\n",
    "should not be ignored. An ill-conditioned Jacobian can render inaccurately computed\n",
    "Newton step dx, creating severe difficulties for the convergence of Newton and\n",
    "Newton-type methods. In some cases, ill-conditioning is a structural feature of the\n",
    "underlying model and cannot be eliminated. However, in many cases, ill-conditioning\n",
    "is inadvertently and unnecessarily introduced by the analyst. A common source of\n",
    "avoidable ill-conditioning arises when the natural units of measurements for model\n",
    "variables yield values that vary vastly in order of magnitude. When this occurs, the\n",
    "analyst should consider rescaling the variables so that their values have comparable\n",
    "orders of magnitude, preferably close to unity. Rescaling will generally lead to faster\n",
    "execution time and more accurate results.\n",
    "\n",
    "#### ill-conditioned\n",
    "https://github.com/OpenSourceMacro/BootCamp2017/blob/master/Computation/Wk4_DifIntOpt/ACME_CondStab.pdf\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3.6 Choosing a Solution Method\n",
    "\n",
    "\n",
    "Numerical analysts have special terms that they use to classify the rates at which\n",
    "iterative routines converge. Specifically, a sequence of iterates x(k) is said to converge\n",
    "to $x^{k}$ at a rate of order p if there is constant $C > 0$ such that\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "$$ ||x^{(k+1)} - x^*|| <=  C ||x^{(k)} - x^*||^p $$\n",
    "\n",
    "\n",
    "for sufficiently large k. In particular, the rate of convergence is said to be linear if\n",
    "$C < 1$ and $p = 1$, superlinear if $1 < p < 2$, and quadratic if $p = 2$.\n",
    "\n",
    "\n",
    "\n",
    "The asymptotic rates of convergence of the nonlinear equation solution methods\n",
    "discussed earlier are well known. The bisection method converges at a linear rate\n",
    "with $C = 1/2$.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3.7 Complementarity Problems\n",
    "\n",
    "\n",
    "Many economic models naturally take the form of a complementary problem rather\n",
    "than a rootfinding or fixed point problem. In the complementarity problem, two n-\n",
    "vectors a and b, with a < b, and a function f from $R^n$ to $R^n$  are given, and one must\n",
    "find an n-vector $x \\in [a; b]$, that satisfies\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "$$x_i > a_i \\rightarrow f_i(x) \\forall i = 1,...,n$$\n",
    "\n",
    "$$x_i < b_i \\rightarrow f_i(x) \\forall i = 1,...,n$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "The complementarity conditions require that $f_i(x) = 0$ whenever $a_i < x_i < b_i$. The\n",
    "complementarity problem thus includes the rootfinding problem as a special case in\n",
    "which $a_i = - \\inf$ and $bi = + \\inf$ for all $i$. The complementarity problem, however, is not\n",
    "to find a root that lies within specified bounds. An element $f_i(x)$ may be nonzero at\n",
    "a solution of a complementarity problem, though only if $x_i$ equals one of its bounds.\n",
    "For the sake of brevity, we denote the complementarity problem $CP(f; a; b)$.\n",
    "\n",
    "\n",
    "Complementarity problems arise naturally in economic equilibrium models. In\n",
    "this context, $x$ is an n-vector that represents the levels of certain economic activities.\n",
    "For each $i = 1; 2; : : : ; n$, $a_i$ denotes a lower bound on activity $i$, $b_i$ denotes an upper\n",
    "bound on activity $i$, and $f_i(x)$ denotes the marginal arbitrage profit associated with\n",
    "activity $i$. Disequilibrium arbitrage profit opportunities exist if either $x_i < b_i$ and\n",
    "$f_i(x) > 0$, in which case an incentive exists to increase $x_i$, or $x_i > a_i$ and $f_i(x) < 0$, in\n",
    "which case an incentive exists to decrease $x_i$. An arbitrage-free economic equilibrium\n",
    "obtains if and only if $x$ solves the complementarity problem $CP(f; a; b)$.\n",
    "\n",
    "\n",
    "Complementarity problems also arise naturally in economic optimization models."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3.8 Complementarity Methods\n",
    "\n",
    "\n",
    "\n",
    "Although the complementarity problem appears quite diffierent from the ordinary\n",
    "rootfinding problem, it actually can be reformulated as one. In particular, x solves the\n",
    "complementarity problem $CP(f; a; b)$ if and only if it solves the rootfinding problem\n",
    "\n",
    "\n",
    "$$ \\hat{f}(x) =min(max(f(x), a - x), b - x) = 0$$\n",
    "\n",
    "\n",
    "(source, http://www.karenkopecky.net/Teaching/eco613614/Notes_ComplementarityMethods.pdf)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "\n",
    "\n",
    "The equivalence, however, can easily be demonstrated graphically for the univariate complementarity problem.\n",
    "\n",
    "Figure 3.7 illustrates minmax rootfinding formulation of the same four univariate complementarity problems examined in Figure 3.6. In all four plots, the curves y = a — x and y = b — x are drawn with narrow dashed lines, the curve y = /(x) is drawn with a narrow solid line, and the curve y = f(x) is drawn with a thick solid line; clearly, in all four figures, / lies between the lines y = x — a and y = x — b and coincides with f inside the lines. In Figure 3.7a, $f(a) < 0$ and the unique solution to the complementarity problem is $x^* = a$, which coincides with the unique root of f; in Figure 3.7b,$ f (b) > 0$ and the unique solution to the complementarity problem is $x^* = b$, which coincides with the unique root of f ; in Figure 3.7c, $f (a) > 0 > f (b)$ and the unique solution to the complementarity problem lies between a and b and coincides with the unique root of f (and f). In Figure 3.7d, f is upwardly sloped and possesses multiple roots, all of which, again, coincide with roots of f .\n",
    "\n",
    "![](https://www.ajjacobson.us/bellman-equation/images/2348_33_7.jpg)\n",
    "\n",
    " f'<0, f(a)>0>f(b)\n",
    " \n",
    " \n",
    "![](https://www.ajjacobson.us/bellman-equation/images/2348_33_8.jpg)\n",
    "\n",
    "f'<0, f(b)>0\n",
    "\n",
    "![](https://www.ajjacobson.us/bellman-equation/images/2348_33_9.jpg)\n",
    "\n",
    "b) f'<0, f(b)>0\n",
    "\n",
    "![](https://www.ajjacobson.us/bellman-equation/images/2348_33_10.jpg)\n",
    "\n",
    "d) f'>0\n",
    "\n",
    "The reformulation of the complementarity problem as a rootfinding problem suggests that it may be solved using standard rootfinding algorithms, such as Newton's method. To implement Newton's method for the minmax rootfinding formulation requires computation of the Jacobian J of f . The ith row of J may be derived directly from the Jacobian J of f:\n",
    "\n",
    "$$\\hat J_i (x) = \\begin{cases}\n",
    " & J_i (x),  \\text{ for } a_i - x_i < f_i(x)<b_i -x_i \\\\ \n",
    " & - I_i, \\,\\;\\; \\:\\:  \\text{otherwise} \n",
    "\\end{cases}$$\n",
    "\n",
    "Here, Ii. is the ith row of the identity matrix.\n",
    "\n",
    "\n",
    "\n",
    "(source https://www.ajjacobson.us/bellman-equation/complementarity-methods.html\n",
    "\n",
    "http://www.codecogs.com/latex/eqneditor.php)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# def cournot(q):\n",
    "#     c = np.array([0.6, 0.8])\n",
    "#     eta = 1.6\n",
    "#     e = -1 / eta\n",
    "#     fval = np.sum(q)**e + e * np.sum(q)**(e-1)*(q) - np.diag(c).dot(q)\n",
    "#     fjac = e*np.sum(q)**(e-1)*(np.ones([2,2]))+e*np.sum(q)**(e-1)*(np.eye(2))+\\\n",
    "#             (e-1)*e*np.sum(q)**(e-2)*(q).dot(np.array([1, 1]))-np.diag(c)\n",
    "#     return fval, fjac\n",
    "# x0= np.array([0.2, 0.2])\n",
    "# f = cournot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# user-supplied \n",
    "x0, maxit, tol , maxsteps =np.array([2]), 10000,  1/1000, 100 # x0 = 0 converge to 0, x0 = 2 converge to 2.005"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "f =lambda x: (1.01 - (1 - x) ** 2,  2 * (1 - x))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#%pylab inline\n",
    "#%pylab notebook\n",
    "# pylab Populating the interactive namespace from numpy and matplotlib\n",
    "# numpy for numerical computation\n",
    "# matplotlib for ploting\n",
    "from matplotlib import pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "\n",
    "a = 0\n",
    "b = np.inf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.LineCollection at 0x7f3d91469a58>"
      ]
     },
     "execution_count": 55,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x = np.linspace(-0.5,2.5,100)\n",
    "plt.plot(x,np.fmin(np.fmax(f(x)[0], a-x),b-x))\n",
    "plt.hlines(y = 0, xmin=-0.5, xmax=2.5, colors='r', linestyles='--')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 2.00498756])"
      ]
     },
     "execution_count": 57,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x = x0 # x0 = 0 converge to 0, x0 = 2 converge to 2.005\n",
    "for it in np.arange(maxit):\n",
    "    fval,fjac = f(x)\n",
    "    fhatval = np.fmin(np.fmax(fval, a-x),b-x)\n",
    "    fhatjac = -np.identity(x.size)\n",
    "    i = (fval> (a-x)) & (fval<(b-x))\n",
    "    if np.any(i):\n",
    "        fhatjac[i] = fjac[i]\n",
    "    x = x - np.linalg.inv(fhatjac).dot(fhatval)\n",
    "    if np.linalg.norm(fhatval)<tol: \n",
    "        break\n",
    "x        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 91,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def myncpsolve(f, a,b,x0,maxit = maxit, tol = tol):\n",
    "    x = x0\n",
    "    for it in np.arange(maxit):\n",
    "        fval,fjac = f(x)\n",
    "        fhatval = np.fmin(np.fmax(fval, a-x),b-x)\n",
    "        fhatjac = -np.identity(x.size)\n",
    "        i = (fval> (a-x)) & (fval<(b-x))\n",
    "        if np.any(i):\n",
    "            fhatjac[i] = fjac[i]\n",
    "        x = x - np.linalg.inv(fhatjac).dot(fhatval)\n",
    "        if np.linalg.norm(fhatval)<tol: \n",
    "            break\n",
    "    return x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 165,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 2.00498756])"
      ]
     },
     "execution_count": 165,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x0=np.array([1.40])\n",
    "myncpsolve(f, a,b,x0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 93,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 2.00498756])"
      ]
     },
     "execution_count": 93,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x0=np.array([2.5])\n",
    "myncpsolve(f, a,b,x0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 166,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-0.005])"
      ]
     },
     "execution_count": 166,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x0=np.array([15])\n",
    "myncpsolve(f, a,b,x0)\n",
    "# overshooting"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "Using Newton's method to find a root of f will often work well. However, in many cases, the **nondifferentiable kinks** in f create difficulties for Newton's method, undermining its ability to converge rapidly and possibly even causing it to cycle. One way to deal with the kinks is to replace f with a function that has the same roots, but is smoother and therefore less prone to numerical difficulties. One function that has proven very effective for solving the complementarity problem in practical applications is Fischer's function \n",
    "\n",
    "$$\\tilde{f}(x)  =  \\phi^- (\\phi^+(f(x),a - x),b - x)$$\n",
    "\n",
    ", where $$\\phi^{\\pm}(u,v) = u_i + v_i \\pm \\sqrt{u_i^2 + v_i^2}$$\n",
    "\n",
    "In Figures 3.8a and 3.8b, the functions $\\hat{f}$ and $\\tilde{f}$ , respectively, are drawn as thick solid lines for a representative complementarity problem. Clearly, $\\hat{f}$ and $\\tilde{f}$ can differ substantially. What is important for solving the complementarity problem, however, is that $\\hat{f}$ and $\\tilde{f}$ possess the same signs and roots and that  $\\tilde{f}$ is smoother than $\\hat{f}$.\n",
    "\n",
    "\n",
    "a) Minimax Formulation \n",
    "\n",
    "![](https://www.ajjacobson.us/bellman-equation/images/2348_33_11.jpg)\n",
    "\n",
    "b) Semismooth Formulation\n",
    "\n",
    "![](https://www.ajjacobson.us/bellman-equation/images/2348_33_12.jpg)\n",
    "\n",
    "(source https://www.ajjacobson.us/bellman-equation/complementarity-methods.html)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Figure 3.9a displays $\\hat f$ (dashed) and $\\tilde f$ (solid) for the complementarity problem and Figure 3.9b magnifies the plot near the origin, making it clear why the problem is hard. Newton's method starting at any value slightly less than 1 will tend to move toward 0. In order to avoid convergence to this false root, Newton's method must take a sufficiently large step to exit the **region of attraction**. This will not happen with $\\hat f$ because 0 poses an upper bound on the positive Newton step. With $\\tilde f$ , however, the function is smooth at its local maximum near the origin, meaning that the Newton step can be very large.\n",
    "\n",
    "![](https://www.ajjacobson.us/bellman-equation/images/2348_33_13.png)\n",
    "\n",
    "(source https://www.ajjacobson.us/bellman-equation/complementarity-methods.html)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 184,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#https://github.com/randall-romero/CompEcon-python/blob/master/compecon/nonlinear.py\n",
    "\n",
    "def fischer(u, v, plus=True):\n",
    "    \"\"\" Computes Fischer's function\n",
    "    phi^±(u, v) = u + v ± sqrt(u^2 + v^2)\n",
    "    In turn, it is assumed that u and v are functions of x. If the Jacobian of Fischer's\n",
    "    function wrt x is required, then partial derivatives du and dv are required.\n",
    "    Fischer's function is useful to transform a complementarity problem into\n",
    "    a nonlinear root-finding problem.\n",
    "    Args:\n",
    "        u:    first term\n",
    "        v:    second term\n",
    "        plus: if True (default), compute  u + v + sqrt(u^2+v^2), else u + v - sqrt(u^2+v^2)\n",
    "    Returns:\n",
    "        phi:                  if either du or dv is None\n",
    "\n",
    "    References:\n",
    "        Miranda and Fackler 2002 Applied Computational Economics and Finance, pp. 49-50\n",
    "    \"\"\"\n",
    "    s = 1 if plus else -1\n",
    "    sq = np.sqrt(u * u + v * v)\n",
    "    ftildeval = u + v + s * sq\n",
    "    return ftildeval\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 185,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "f = lambda x: (1.01 - (1 - x) ** 2, 2 * (1 - x))\n",
    "a = 0\n",
    "b = np.inf\n",
    "x0=np.array([2.5])\n",
    "x = x0\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 186,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "\n",
    "for it in np.arange(maxit):\n",
    "    fval,fjac = f(x)\n",
    "    \n",
    "    ftildeval = fval\n",
    "    da = a-x\n",
    "    db = b-x\n",
    "    hasLowerBound = np.isfinite(a)\n",
    "    hasUpperBound = np.isfinite(b)\n",
    "    if np.any(hasLowerBound):  # apply the Fischer + transform\n",
    "        ftildeval= fischer(ftildeval, da)\n",
    "    if np.any(hasUpperBound):  # apply the Fischer - transform\n",
    "        ftildeval = fischer(ftildeval, db, plus=False)\n",
    "    \n",
    "    ftildejac = -np.identity(x.size)\n",
    "    i = (fval> (a-x)) & (fval<(b-x))\n",
    "    if np.any(i):\n",
    "        ftildejac[i] = fjac[i]\n",
    "    x = x - np.linalg.inv(ftildejac).dot(ftildeval)\n",
    "    if np.linalg.norm(ftildeval)<tol: \n",
    "        break"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 187,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 2.00498756])"
      ]
     },
     "execution_count": 187,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Performance does not improve than CP ？？"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 188,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def myssmooth(f, a,b,x0,maxit = maxit, tol = tol):\n",
    "    \"\"\"f gives two outputs, value; jacobian\"\"\"\n",
    "    x = x0\n",
    "    for it in np.arange(maxit):\n",
    "        fval,fjac = f(x)\n",
    "\n",
    "        ftildeval = fval\n",
    "        da = a-x\n",
    "        db = b-x\n",
    "        hasLowerBound = np.isfinite(a)\n",
    "        hasUpperBound = np.isfinite(b) # only apply Fisher if there is a finite bound  ??\n",
    "        if np.any(hasLowerBound):  # apply the Fischer + transform\n",
    "            ftildeval= fischer(ftildeval, da)\n",
    "        if np.any(hasUpperBound):  # apply the Fischer - transform\n",
    "            ftildeval = fischer(ftildeval, db, plus=False)\n",
    "\n",
    "        ftildejac = -np.identity(x.size)\n",
    "        i = (fval> (a-x)) & (fval<(b-x))\n",
    "        if np.any(i):\n",
    "            ftildejac[i] = fjac[i]\n",
    "        x = x - np.linalg.inv(ftildejac).dot(ftildeval)\n",
    "        if np.linalg.norm(ftildeval)<tol: \n",
    "            break\n",
    "    return x    \n",
    "\n",
    "# source\n",
    "#https://github.com/randall-romero/CompEcon-python/blob/master/compecon/nonlinear.py\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 189,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "f = lambda x: (1.01 - (1 - x) ** 2, 2 * (1 - x))\n",
    "a = 0\n",
    "b = 10\n",
    "tol = 1/10000\n",
    "maxit = 10000"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 190,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.00448019])"
      ]
     },
     "execution_count": 190,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x0=np.array([1.1]) ##??? it does not improve from fhat\n",
    "myssmooth(f, a,b,x0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 191,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 2.00498756])"
      ]
     },
     "execution_count": 191,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x0=np.array([2.5])\n",
    "myssmooth(f, a,b,x0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 192,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.00448019])"
      ]
     },
     "execution_count": 192,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x0=np.array([15])\n",
    "myssmooth(f, a,b,x0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 88,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Reference\n",
    "\n",
    "\n",
    "scipy.optimize does not have ncpsolve, but dolo python package provides a ncpsolve.\n",
    "\n",
    "import scipy.optimize as opt\n",
    "\n",
    "https://github.com/EconForge/econforge/wiki/CompEcon\n",
    "\n",
    "\n",
    "https://github.com/EconForge/dolo/blob/master/trash/dolo/numeric/solver.py\n",
    "\n",
    "\n",
    "https://github.com/EconForge/dolo/blob/master/dolo/numeric/optimize/ncpsolve.py\n",
    "\n",
    "\n",
    "\n",
    "http://nbviewer.jupyter.org/github/sbustamante/ComputationalMethods/blob/master/material/one-variable-equations.ipynb\n",
    "\n",
    "\n",
    " - Bisection Method\n",
    "\n",
    " - Fixed-point Iteration\n",
    "\n",
    " - Newton-Raphson Method\n",
    "\n",
    " - Secant Method\n",
    " \n",
    " \n",
    "https://github.com/dingliumath/economics/blob/master/3.2%20-%20Finding%20Roots%2C%20Newton's%20Method.ipynb \n",
    " \n",
    "https://github.com/dingliumath/economics/blob/master/3.3%20-%20Ramsey%20via%20Line%20Search%2C%20Shooting%20Method.ipynb\n",
    "\n",
    "https://github.com/dingliumath/economics/blob/master/3.2%20-%20Finding%20Roots%2C%20Built-in%20Methods.ipynb \n",
    "\n",
    "\n",
    "https://github.com/EconForge/econforge/wiki/CompEcon\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}