nmayorov/bundle_adjustment_scipy.ipynb

## bundle_adjustment_scipy.ipynb
{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Large-scale bundle adjustment in scipy"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A bundle adjusmtent problem arises in 3-D reconstruction and it can be formulated as follows (taken from https://en.wikipedia.org/wiki/Bundle_adjustment):\n",
    "\n",
    "> Given a set of images depicting a number of 3D points from different viewpoints, bundle adjustment can be defined as the problem of simultaneously refining the 3D coordinates describing the scene geometry as well as the parameters of the relative motion and the optical characteristics of the camera(s) employed to acquire the images, according to an optimality criterion involving the corresponding image projections of all points."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "More precisely. We have a set of points in real world defined by their coordinates $(X, Y, Z)$ in some apriori chosen \"world coordinate frame\". We photograph these points by different cameras, which are characterized by their orientation and translation relative to the world coordinate frame and also by focal length and two radial distortion parameters (9 parameters in total). Then we precicely measure 2-D coordinates $(x, y)$ of the points projected by the cameras on images. Our task is to refine 3-D coordinates of original points as well as camera parameters, by minimizing the sum of squares of reprojecting errors."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let $\\pmb{P} = (X, Y, Z)^T$ - a radius-vector of a point, $\\pmb{R}$ - a rotation matrix of a camera, $\\pmb{t}$ - a translation vector of a camera, $f$ - its focal distance, $k_1, k_2$ - its distortion parameters. Then the reprojecting is done as follows:\n",
    "\n",
    "\\begin{align}\n",
    "\\pmb{Q} = \\pmb{R} \\pmb{P} + \\pmb{t} \\\\\n",
    "\\pmb{q} = -\\begin{pmatrix} Q_x / Q_z \\\\ Q_y / Q_z \\end{pmatrix} \\\\\n",
    "\\pmb{p} = f (1 + k_1 \\lVert \\pmb{q} \\rVert^2 + k_2 \\lVert \\pmb{q} \\rVert^4) \\pmb{q}\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The resulting vector $\\pmb{p}=(x, y)^T$ contains image coordinates of the original point. This model is called \"pinhole camera model\", a very good notes about this subject I found here http://www.comp.nus.edu.sg/~cs4243/lecture/camera.pdf"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "--------------"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's start solving some real bundle adjusment problem. We'll take a problem from http://grail.cs.washington.edu/projects/bal/."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import urllib\n",
    "import bz2\n",
    "import os\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "First download the data file:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "BASE_URL = \"http://grail.cs.washington.edu/projects/bal/data/ladybug/\"\n",
    "FILE_NAME = \"problem-49-7776-pre.txt.bz2\"\n",
    "URL = BASE_URL + FILE_NAME"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "if not os.path.isfile(FILE_NAME):\n",
    "    urllib.request.urlretrieve(URL, FILE_NAME)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now read the data from the file:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def read_bal_data(file_name):\n",
    "    with bz2.open(file_name, \"rt\") as file:\n",
    "        n_cameras, n_points, n_observations = map(\n",
    "            int, file.readline().split())\n",
    "\n",
    "        camera_indices = np.empty(n_observations, dtype=int)\n",
    "        point_indices = np.empty(n_observations, dtype=int)\n",
    "        points_2d = np.empty((n_observations, 2))\n",
    "\n",
    "        for i in range(n_observations):\n",
    "            camera_index, point_index, x, y = file.readline().split()\n",
    "            camera_indices[i] = int(camera_index)\n",
    "            point_indices[i] = int(point_index)\n",
    "            points_2d[i] = [float(x), float(y)]\n",
    "\n",
    "        camera_params = np.empty(n_cameras * 9)\n",
    "        for i in range(n_cameras * 9):\n",
    "            camera_params[i] = float(file.readline())\n",
    "        camera_params = camera_params.reshape((n_cameras, -1))\n",
    "\n",
    "        points_3d = np.empty(n_points * 3)\n",
    "        for i in range(n_points * 3):\n",
    "            points_3d[i] = float(file.readline())\n",
    "        points_3d = points_3d.reshape((n_points, -1))\n",
    "\n",
    "    return camera_params, points_3d, camera_indices, point_indices, points_2d"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "camera_params, points_3d, camera_indices, point_indices, points_2d = read_bal_data(FILE_NAME)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we have numpy arrays: \n",
    "\n",
    "1. `camera_params` with shape `(n_cameras, 9)` contains initial estimates of parameters for all cameras. First 3 components in each row form a rotation vector (https://en.wikipedia.org/wiki/Rodrigues%27_rotation_formula), next 3 components form a translation vector, then a focal distance and two distortion parameters.\n",
    "2. `points_3d` with shape `(n_points, 3)` contains initial estimates of point coordinates in the world frame.\n",
    "3. `camera_ind` with shape `(n_observations,)` contains indices of cameras (from 0 to `n_cameras - 1`) involved in each observation.\n",
    "4. `point_ind` with shape `(n_observations,)` contatins indices of points (from 0 to `n_points - 1`) involved in each observation.\n",
    "5. `points_2d` with shape `(n_observations, 2)` contains measured 2-D coordinates of points projected on images in each observations.\n",
    "\n",
    "And the numbers are:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "n_cameras: 49\n",
      "n_points: 7776\n",
      "Total number of parameters: 23769\n",
      "Total number of residuals: 63686\n"
     ]
    }
   ],
   "source": [
    "n_cameras = camera_params.shape[0]\n",
    "n_points = points_3d.shape[0]\n",
    "\n",
    "n = 9 * n_cameras + 3 * n_points\n",
    "m = 2 * points_2d.shape[0]\n",
    "\n",
    "print(\"n_cameras: {}\".format(n_cameras))\n",
    "print(\"n_points: {}\".format(n_points))\n",
    "print(\"Total number of parameters: {}\".format(n))\n",
    "print(\"Total number of residuals: {}\".format(m))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We chose a relatively small problem to reduce computation time, but scipy's algorithm is capable of solving much larger problems, although required time will grow proportionally."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "Now define the function which returns a vector of residuals. We use numpy vectorized computations:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def rotate(points, rot_vecs):\n",
    "    \"\"\"Rotate points by given rotation vectors.\n",
    "    \n",
    "    Rodrigues' rotation formula is used.\n",
    "    \"\"\"\n",
    "    theta = np.linalg.norm(rot_vecs, axis=1)[:, np.newaxis]\n",
    "    with np.errstate(invalid='ignore'):\n",
    "        v = rot_vecs / theta\n",
    "        v = np.nan_to_num(v)\n",
    "    dot = np.sum(points * v, axis=1)[:, np.newaxis]\n",
    "    cos_theta = np.cos(theta)\n",
    "    sin_theta = np.sin(theta)\n",
    "\n",
    "    return cos_theta * points + sin_theta * np.cross(v, points) + dot * (1 - cos_theta) * v"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def project(points, camera_params):\n",
    "    \"\"\"Convert 3-D points to 2-D by projecting onto images.\"\"\"\n",
    "    points_proj = rotate(points, camera_params[:, :3])\n",
    "    points_proj += camera_params[:, 3:6]\n",
    "    points_proj = -points_proj[:, :2] / points_proj[:, 2, np.newaxis]\n",
    "    f = camera_params[:, 6]\n",
    "    k1 = camera_params[:, 7]\n",
    "    k2 = camera_params[:, 8]\n",
    "    n = np.sum(points_proj**2, axis=1)\n",
    "    r = 1 + k1 * n + k2 * n**2\n",
    "    points_proj *= (r * f)[:, np.newaxis]\n",
    "    return points_proj"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def fun(params, n_cameras, n_points, camera_indices, point_indices, points_2d):\n",
    "    \"\"\"Function which computes residuals.\n",
    "    \n",
    "    `params` contains camera parameters and 3-D coordinates.\n",
    "    \"\"\"\n",
    "    camera_params = params[:n_cameras * 9].reshape((n_cameras, 9))\n",
    "    points_3d = params[n_cameras * 9:].reshape((n_points, 3))\n",
    "    points_proj = project(points_3d[point_indices], camera_params[camera_indices])\n",
    "    return (points_proj - points_2d).ravel()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can see that computing Jacobian of `fun` is cumbersome, thus we will rely on the finite difference approximation. To make this process time feasible we provide Jacobian sparsity structure (i. e. mark elements which are known to be non-zero):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from scipy.sparse import lil_matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def bundle_adjustment_sparsity(n_cameras, n_points, camera_indices, point_indices):\n",
    "    m = camera_indices.size * 2\n",
    "    n = n_cameras * 9 + n_points * 3\n",
    "    A = lil_matrix((m, n), dtype=int)\n",
    "\n",
    "    i = np.arange(camera_indices.size)\n",
    "    for s in range(9):\n",
    "        A[2 * i, camera_indices * 9 + s] = 1\n",
    "        A[2 * i + 1, camera_indices * 9 + s] = 1\n",
    "\n",
    "    for s in range(3):\n",
    "        A[2 * i, n_cameras * 9 + point_indices * 3 + s] = 1\n",
    "        A[2 * i + 1, n_cameras * 9 + point_indices * 3 + s] = 1\n",
    "\n",
    "    return A"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we are ready to run optimization. Let's visualize residuals evaluated with the initial parameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "x0 = np.hstack((camera_params.ravel(), points_3d.ravel()))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "f0 = fun(x0, n_cameras, n_points, camera_indices, point_indices, points_2d)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x107584828>]"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYMAAAEACAYAAABRQBpkAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztnXm4JVV1t98f8ygIKI3QsVFBAUWFMIlII0PQKJgYRSMG\no5JPUSGJgoD5pI2JiiYqmE8TEyVEBUWiDI60aKNEQEBAoMEGBaVRGgwK4siwvj+qzr116tY8nKo6\nZ73Pc59btWvX3uucU7V/e1xbZobjOI4z26zVtQGO4zhO97gYOI7jOC4GjuM4jouB4ziOg4uB4ziO\ng4uB4ziOQwNiIGlzSedKuknSSkl7SdpC0nJJqyRdJGnzJox1HMdx2qGJlsFpwJfMbCdgV+Bm4ERg\nuZntCFwcnjuO4zg9RXUWnUnaDLjGzJ4QC78Z2N/M1khaBKwws6fUM9VxHMdpi7otg+2BeySdIem7\nkv5d0sbA1ma2JoyzBti6Zj6O4zhOi9QVg3WA3YAPm9luwK+IdQlZ0PRwnxeO4zg9Zp2a968GVpvZ\nleH5ucBJwF2SFpnZXZK2Ae6O3yjJBcJxHKcCZqam06wlBmFhf4ekHc1sFXAQcGP4dxRwavj/vJT7\nG/9Ak0LSMjNb1rUdVXH7u8Xt744h2w7tVaTrtgwA3gR8StJ6wA+AvwTWBs6R9BrgduClDeTjOI7j\ntERtMTCz64A9Ei4dVDdtx3EcZzL4CuTqrOjagJqs6NqAmqzo2oCarOjagJqs6NqAGqzo2oA+Umud\nQa2MJRvymIHjOE4XtFV2esvAcRzHcTFwHMdxXAwcx3EcXAwcx3EcXAwcx3EcXAwcx3EcXAwcx3Ec\nXAwcx3EcXAwcx3EcXAwcx3EcXAwcx3EcXAwcx3EcXAwcx3EcXAwcx3EcXAwcx3EcXAwcx3EcXAwc\nx3EcXAwcx3EcXAwcx3EcXAwcx3EcXAycHiFxusQTurbDcWYRFwOnT7wJeHHXRjjOLOJi4DiO40yP\nGEjs1bUNTiOs3bUBjjOLTIUYSGwKXN61HU4jPL1rAxxnFpkKMWB6PodTAIlndW2D40wbjRSiktaW\ndI2kC8PzLSQtl7RK0kWSNm8iH8eRWBv4n67tcJxpo6ka9XHASsDC8xOB5Wa2I3BxeO7MMBKfltgq\nPF5HmntWHMfpAbXFQNJ2wPOB/wAUBh8GnBkenwm8qG4+OXjB0n+OAPYIj7MGiZVxzXGclmiiZfAB\n4HjgkUjY1ma2JjxeA2zdQD6OAy78jtMKtcRA0guAu83sGlJqdGZmVHiBJUziwDr2FUjfRWqyeK3f\ncXrKOjXvfxZwmKTnAxsAj5L0CWCNpEVmdpekbYC7k26WtCxyusLMVsSijMYc8qhaW1xE0HJxJsvi\nuglIyKzdVoLExmb8qs08HCcPSUuBpW3nU0sMzOxk4GQASfsDbzGzV0p6L3AUcGr4/7yU+5fVyX+I\nSLwKOMNspmvJ22dcm/teJDYB3m7GCe2blMgDEhsC/xdYZTY3DuY4EyOsJK8YnUs6pY18mp6fP6qp\nvQc4WNIq4LnheRXaLjC76H/eIz/K1DKqfDycESf6m+9OMB7VJWsTVHhO7tgOx2mVut1Ec5jZJcAl\n4fG9wEFNpV0k+wnm5ZREWuB8LksMCifLZH73WW7BOTNE31fu+os4HZwbO58TA4mdw4VkI/oq7P4s\nOlNN38XAmS62CP9HWwY3AkdGzgUgsYRQGCSSVrB74ew4DTItYtDX2qQzzsbh/4di4RtFjiWxPnBb\nJGz3Vq3KRrH/jjOV9F0MpnEAGQCJbSU26Cr/Iki8WeJ3kfPnSI38Jo/EzqNpGrBd7HrSepCNEsIc\nx6lI38VgmlkNvLdrI3LYE1hP4smhL6FLgF3LJCDxqITgLBH+A+ZbEKN4SQJ0v8Q+YR6LJU4rY1cJ\nvEXgzATTIgZVa/hddy89puP8i7IkcnxlyXs3jByPFayRVsaREq8Pj/di/nf5XCTubyS+EUv7seH/\nFwHHlrSrLJ2JgoRSRNVxGmNaxGCoDKXWGRXNdWukk/Z59wE+nBC+VeS+DQhWYU76Oxvl98QJ5xvl\nxcB9HebvzACDFgOJDSRe0rUd00BY+9wlPN5U4nyabzkVLciz8o0u2huKmNZl264NcKafQYsB8MfA\nOcxPQby6W3NK06fCbH/ghvB4RwI35CPa6k4r+vmPihxH1yR8vkFb0ujTb+Q4rdF3MSj7Iu7WihWz\nwYb5UWpTtWUQXc2elEanYz8Si9rarEfiIOj3rDNnOmjMHUVbSKwF7GTGjS0k3/UAcu9qnRIvZX5R\nmMX+O8m00o0jcQXBjK7VbaTvOFH6LgYiGDw7h+waoRdWFQj3c3iUGbdEgj/DfHdRa1nH/veZLm3c\nM/wfX3fhOI3T924imO7FRV0Xhl8CViWExwvrOmIbvbeJz9u7bqIe5O84tRmCGGRR9yWc9Zd4k5Tw\nJkVKKcdZ9Ol36VqwHWci9F0MJuWmuCv6WtDEWwRtTTHt6+ePsn7XBjjOJOi7GOQxzULRJfFCuqlu\noiboatFZFv4cOoNnCGLgL9r00FZB7s+I49RkCGKQxdDHDPraTTKya+3MWOXSih9XZZexxNWIjXXp\n+jlynNoMQQz8RZs8o+eibbGqkn7ch9HBTRiSwZyNDbnvdpxe0ncxyHv5hi4UXRcued9fE89Hld9o\ncYm4c2tlWmolPCNyHN+HYcTEnkOJM/q+D4YzTPouBjD8Aj9O9PO8qDMrijGpFkKcqu6a27CzN66j\nw5bJq4AvSOyUEucJUrk9JxwH+r8COY8hCkXTfeht0nRlIf550z5/me+lD9/npJ7DkduLA4GVJH/e\nbxBsENT3Z8vpGX1vGbS9zmCIYtIWfxA5Hm1c3/Tz0fb3Pe0F4B0F4vT9nR4EEte05Xywrwz9wZmp\nH6tl/jVyvGX4v+nCtejzVuZ3naWWQRFq2SJxk8TrmjJmwDwjP8p0MQQxKPVwS5zQVtoJeX1F4rg6\nafScpgvXot1ETaXfBH0q6CfBU2h/hpbTQ/ouBkVnE60XCTsVQGp+m0IJk9grEvRHwEtrprmWNFcT\nnzRFZxPVKWQt5bgp+tAy6BOVvmOJAzt8Dp0e0HcxgGIP92bRE4mnAre2Y07yLI44EpdKHFIg6tHA\nz+qZVJm8wrOJ52OShXV0TcA60vhz0SLT0Hr4Gt09h71FYseubZgUQxCDrDnVaS/hf6fdIHGqxKb1\nTCrEvgTbcubxOACJDSW2k/ihxBPD/Z273vs2t/CW2CrcF6FKek27tI4evwv4RQPpD41pEKY+sXHX\nBkyKIYjBevlRFpBVyJzA/KYhVag67XFE2sv6QYLZItsTDF79BlgtsRWAxEclnl/G0AIU7SYaiyeN\nbZH5HeD7TRpVg+j3vX1rmYgbJBZFggbt2VTiufGgTgzpmLAStig/5nRSSwwkLZb0DUk3SrpB0rFh\n+BaSlktaJekiSZtXzYLsAivtWpt+8+uKQRo7p9x3T/j/aOCvJJ4isV+JdOuQNmYQtXUbyOyOmfiY\nQbhC96nh8VE100+qOOzCuI+kZ9bMo0kKf8cSrwmnT17coj1D4lbgzlhYmdXwg6Zuy+BB4G/MbBdg\nb+ANknYCTgSWm9mOBA/aiRXTfz7Fulri5BXCY5u/SyyeoN+ZNPfQz86IM+Jw4Cbgm00blUIRtwd1\nBLWJRWdJ9x1PMCsG4O0V0xpRpEuxT10zZWz5j6qZhN2aVVrtVfM7Q2plH/Q48TLx/Ank2QtqiYGZ\n3WVm14bHDxAUVNsChwFnhtHOpLrbhUcDj80yISU8rzC5MHb/j4E/K2hT0y2DpDhPL5FHm+yTEh61\nuUzh00a3ZNKYwUYp151iFBk0vYNgv+xMJDaVWLvqOJ3EOyUeQ+CGY+ec6JWQOEaaHzyXgnG8WaMx\ndxSSlhA0l68AtjazNeGlNVB4gDGJNrtyovEe3XDadXj8WIYJKyEltjTjfydgSx5l9qie1Gyi6PdV\nd+wg7flru/urKk3Yskt+FLYk7IrL4f7RgQRmxZ4BiSsJ9uf+c9qbGThiPxibVps6BiTxC2AnM37a\nsk0TpxExkLQJwQye48zsl9L8721mJinxAZW0LHK6wsxWxKNQ7eFeUABIrE92t8eBwEcr5JVF1V2y\njixw388Kpt8EpfKReANws9lE+qJ9ncF08ofAduFxrjfacFLDhmbc26pVwfjYE2ByYiBpKbC07Xxq\ni4GkdQmE4BNmdl4YvEbSIjO7S9I2wN1J95rZsgJZzLkNDvsoNzKbmzJYRijOBI6ImxA5fmnC9SQS\nCxyJfYBrzPhtCZtqIbGjGatqJGFhOp9syKQR/wJcTfBCZ9HW1NI+1dQnTa8/u4TM0m2U+GvgtFhw\nZveixAuBjxB0URdteawNPBS2VMrO0ptopSOsJK+Yy1w6pY186s4mEvAxYKWZfTBy6QKYm8VxFHBe\n/N4SRB+c04GfV0ynSNMXAInjJV6Zdjl2/qzw/7eB/5MTFxa+rE8ralcCRZrpWYzse0XF+7LoYgey\nSYpBX7uJOkXishwHb6/PSeIDlK+kXgCl1+REy76ybsoHPZU4jboDevsSdGkcIOma8O9Q4D3AwZJW\nAc8Nz6sgxqctxrt/lob/M3+ccB51UsGZ9tC+l3SbswrCdUvEHbFFfpRBMqk1LN5NNM4khelJCWF7\n59zzrjYMqUB0JlT8O3sW2dQZA+0ttbqJzOxS0l/6g8qmJ7Ex8BQzro4ER1U7/qONZj0kFqgS65vx\nO+Cfy9qSwQsk1gBzTbWI2wNJvBrm3FC8TuI3wEmhHQBvaNCWuoVf0YKjysrhpGmHk3JUN+nFlH1q\nGfTJliQ2zI9S/DmJr7mROMyMCwrk8UDk2GJ57FHg/qmjbyuQTwKuioVVWXQ24rDwf5o72qz7x6aX\nRdYhbA48L5Zm1Fnd0YyPPRwL/DacXrdbSRvy6HNN+CkJYXF7y8xESkORdEf/j20gXacdynQfjlYD\n/01GnPiam1LrAqQFrXlgqj0Rp9I3MYh394hg5H5UGDdd66lamMbt2DeSXryZPHr414KxFk8TTEoM\n4p+3zO+QtQK46Zdu9H205k9GSpxb36fa+ERtkVgv5sk3jyplTtoWn5V6NmJ7SJ9Bv36/zuiFGEis\nK/EqFv4o0b1ck8Qg70fcoMwDI7GbxPvzoiXk+0jkWm8Jfa8c3URSJeK+L+O+pFpZ2fQ3JJi9VPa+\nqtRyWT6FvAa4XGJnaWy3vGhreiy4QJpFHCSuDdlegSUelfL+RzdyyhsfSEy6wj29pxdiALycQKE3\nyYm3CYDEXqHv9Twx+C8yPJgm3P9/iTVJJbaIOa/aN+G+OmKQN/0yiwX5SWwszXfRSDxWYvfw9ATg\noxUEoeh0vT0zZmElpdfE8xf9/iSVH6uqSOJsIolXSFwyIRuSaL2WK42txxkJ+o0E+y9PiiOBL+bE\nuQ/4+4TwxyeElcHFoA0kHs2864qsaWeCucGiywkWXRV58A/LuihxReQ0yW3GxcBPYmFxO18TsbEs\nddxpJ+X3bgK3IEi8j2AF+FXhdL+/CuPUXVyX9jlPIxDguumU4Y2R461YOL130vwJ8JyObWibP4kc\nR9/BIr6smqLoe5O3Av3xeDcR0KA7iipIXAX8oGj0Fkww8t1Zl9kLddLimvSdRNdTvGVihiQPxC2M\nFnTDHR6ep91TderemXTjQbRPaw4mvcYiOiDclE+fIu/6w6k3i9eaFXbCt1ZWWrNEp2IA7E69AnT0\nUPbFB/kkNs3JI+6bflIUnT8e7YZL23h915TwPNoQgrRuvLRC95GU8EnRpRjEKevGpExZkFV2RZ+D\nIgs/y35n3k3UEjvEAyS2SYiX9aMWqZUmcWDWxQq+8IvMoW6SrnwTJeW7XUJYHhNzgVyWyGBomsC8\nRsGe2BdG7tmC0OGhVGlgsveE8/qjnoT/qcHkRzX0pPc/TlEPpksSwg6InbsY0A8xSBo0/lBCWJld\nw4qS13f+L5HjIgOjVV11V+UsiQ/mR2ucssv3RxR9iRZUENpGWjCo+KOcW0bPwwuYfw6vBA4Oj49v\nyLSytN0ymMReGkUG4LMmm0R7CopMe83zYDCVhX+cPohBEkUXI728VSvGOTM/yvjUuglRZ65+1YJj\nA4AWhaiLAdjbJZ5cM43o7z/pisE0UWR3sawKSdYeKFVwMeiQJNUv4xNoppC4UOKvpdJ95lvmR0nk\nfIkfsFCI8sTl/pzrE0Xis9LY7KMNJPbPcbSWxCh+12Nw0FLLQKrUBdpmIZq0wn1ElS7LLFwMOiSp\nZZD1g5zdkh156x7K8PsG04rzAgJvj2W7V+o43HpCQljazmhFr0+aPwP+InIuIq6C69Ln7gVpbJvV\nvLiLgV+3aE4VslrhSxrOK15O9vZ3rUNfxWD3/ChOArnbEFbkqwXivKxAnKLTiIdGmvvkrNprWxRt\nGZTpYn0UgMS3y5szFcQL/yLdWIOjr2KQxFSq8TQjLXBx3PUc/DyqPmOnNmpFPYp+x1V+i7617CZF\nvJx8ZydWtIyLgdMmL46d93Eqad8Fqi0m8bnbfGcn+bvNRNkzJDGIzw12+k/8hZ3YiugZpmghWWZx\n3KwK5ggBSI2tsO4lQxKDL3RtgFOa3hYiExrc7e3np9+29Y3Rs3Jnp1a0zJDEwBkefS5wRmsK+mxj\nFdocM5hVvJvIcWrSqwJHYpPIFqVJz37T+1HP7PoXiT8FHiwYd72K6xicBnExcNqkSb81TbAC+IXE\nx5iv7UUFq8j02DJ08X61MWZQhax9ROJ8Dri1ZPo+gNwwLgbOLDFalPfqCeV3bUHX3k2yWX4UoFxh\n+pgqhpTg6TTn/noBEgeErY8tJbaqkMRjJK5s3LCe0Yfl844zzVSqcEmJK7yLkLeZy4hcMQg9t/6e\nBldlp2XVcvpfJ1hp/g6q7Y89rYslx3AxcGaV0YK4Xo1rRGi7ACryuX8E/LBlO1pBWrAp1VoEPotm\ndhwnD+8mcmaVpK6PLtxHLECqv0mSVMvvVJSqLZS2yWtN/ENCWGPC32e/U1VxMXBmicTnPbKYqLDz\ntpY5uUikHNF4qcQSiZ0l9k64nrVDGRJvKmJDh+QVxkkD5E22AqdODLybyJkJJJ7KuBfaUcGwFvDx\nyVuUyYkF491PpFAKvYuOOJ1gY5elwLYSa5uNFZB/m5P2IQVtqExYuz6speST9jV2McjAWwZTSmQ+\nvRMQL/BHfpM2oXn/91HGCiCJtSQ2bzIDiS9LnATsErv0CuCh8Ph5YdwDpOxWQUh857fGkDhC4jAC\nb6jnUa1gzSvY224ZTF1F2sVgevlF1sVZWuQTblazRyz4eeH/dVlYiLbJq4GfSxwR2rbZaDMdiRsk\nPlAhzUOBdwFfTrg2EoMvSLyFYGbNfgXSfGIFO4ryaeBTMLfHeOa2k3Ek/oT8/cb3jZ3/Z5k8CvDH\nDafXOa2JgaRDJd0s6RZJb20rH6cyfduspCuaXnUc53cSb5CwsND/xzD80+H5SBROJxClA8skXmCw\nOVqovy/8X2Qz+aJbz1ZlE+C08LjwjnsSBxEsUssjaQC9yQrQjg2m1Qtk1vzMOklrA98HDiJw7nQl\n8HIzuykSx/o7q89xHCeXz5rx0klnKsnMrPExi7ZaBnsCt5rZ7Wb2IEGz8PCW8nIcx+mCl3RtQJO0\nJQbbAndEzleTvjWg4ziO0zFtjYgX7P9ZFjleGv45juM4IyQtZQKFY1ticCfjm0YvJmgdxFjWUvaO\n4zjTgZmtIOIfStIpbeTTVjfRVcAOkpZIWo9gxsQFLeXlOI7j1KQVMTCzh4A3Al8FVgKfic4kSiHa\nkpi6ObyOE2G06Yvv6z1czmPK1mm19mHM7Mtm9mQze5KZvTs7LjJL6kZynMGzDFgVC9sOeJYZK5if\nD194rn2MI3Ou310xXSebT5pN19z4virb1Pn9cGaSE8x4B/B24KthpUdm3G3GZQBm3A08yox7gZ0Y\nbyEX4Wzgkkatdorwra4NaJq+ikERbujaAMfJ4X4AMz5jxqFpkcz4Zfj/5rIt5ND53GtrWemUJhTx\nqaKvYpDbMjDjaZMwxHFq8PtJZGJWaP/gk4Anl0z61ArmRCm7r7HTIX0Qg00SwvLE4M1tGOI4DfPL\nrg0ALgWuN+M9wC3An5e4t2537Ydq3u9MkM7FwIxfVbjn/W3Y4jgN88CE8/vnhLCXQLAFpBlmxtkl\n0qsrBitr3u9MkM7FIMZl4X8fQJ4OJtJN0mMmOdtkHeD4BQYYj8Q2tSnDqHzYteL911e8r2k+ApzQ\ntRF9p29icEXD6f19w+k55ZiqqXcVmNjnN+PhFqY6jiplv6t4fy9+fzOOMZtz3126JwJ4eZP29JWu\nxSC+AcWIpJbB+xLCslgCfDhy/o8p8RynLfowZlCH0Xt4X8X7eyEGDfCzrg2YBJ2KgRnfTrmUJAY/\nLpn2j5jf5QngHWXudxph1rv7runagJrIDFF97KOPYtBHm3pB1y2DOD/KuPaRmmkP/SFoc2HRDfi6\njcYx47dd21CTkZhXfXf69s5dBHymwn3Rz3FXQ7b0jr6JwenAZiTUKM14uE7CZmOthCHSdmHdxmBv\nqdac0zvyWna35FzvlRiY8UfAR7u2o6/0SgzCmQ/3d21HTzmz5fTb6NLpomb8+Q7ynFbynonv5lzv\nSgx+0VG+g6ZXYhAh7SGsWjv+n6qG9Agj2Eu6DUTgdnwa+JuG0vm3htIZMnli8IOc612JwUkNp3dn\nw+n1kqGJwf7AzgXu/37s/Mv1zOkFbb5YAl7XQroTLwzCiQN1uY+edXF0RJ4Y5H1HXX2H3yH9eS5t\nU+gzauonQ/RVDBIx414z8vZFgIVikDU1rO9NyuXhf2Phy/n4pjKpsTAp1QHbgIk+Y3n94tNM22KQ\n5pTv9TUL37tJ71Z1kU+hr2JQ9EF4fUr4bwrc+4Xwfx/nEL8kcjwqpBc8xGbVBmjDF22PKvcmpPVV\n0h2gRW2+LCVOk5we/q87SBi1e7ecuDvWzKsOL8mPUou6teG8Ckaau4qvhv/3q5Jp6Pk1rdDPEoP/\nCv9nTTaZ2hZCX8WgKP+bEPYt4JgC965JSePCWhY1gBnnRk7zpvdtCTwN+H/Mv0RF8riKeUGo9YCb\nLdi8Ze5S5Ph7dfIoyGh16a9j4Xl923GMYjXI75t12nL4Wsvp/zTnevQ7ujjnepy/Bd6Ucu3HAGZc\nmpN/FlVaummVy5mgr2JQp3C6LtwoJIkFvlsYf2BXA/9eIq8ycaNkTeP8p9j5eeH/xBcr7Dq7wYw3\nkr5daKLAhYIA7ayUzfNL81bg2Q3nOZqJFv9+492GeRRdcRt/To8tmU8b/KShdB4LZO5QyPgz+cnY\ntQ9lOaE04wNplYi608jDNB5Mu5Rxz68J3G7XEaHBMo1ikCYEAKclhEUfjrKzEOI10KJk1VrGun7M\n5hbb5dZUw5coqYvs9WS3GrK+szT+IZZ3/Dd7euzcGBeIL5olzvJaBbymgj0AH6h4X5yihdHoM49m\nME3aS2dagVcbM+4psDbHCL9zM/4zdn+aMMYFIvqbXVTAtLOYn7Lc+DoWM3YgGICeOfoqBmUZzSDZ\nkQwfRGFtIa8WvDznepQiolWmpvYM5qc0xgWj6MBXtFZ8PrDCjDuzdtrK4UuMT83dCtgQ+HjGPeen\nOE0r2hIoO5Vv9B2XXTj39YxrRb7vMd89Zlw8yVknVdy/VyTrM70FeHR4/NkCab2P8YHjaEs4d3Ww\nGa8wY8Pw9BnhXxmK/K4zOcjcVzEYPXxFagoQ9vubcYvZWIFQ5GVZEzm+LeJCwIBNc+69rkD6SZ9h\nWVJEM66L2T92uUBeEIwDHA7saMaLzDggcu24nHvj3V4nm/HHwH9Hwh4w47dm3Eb68zP6XsZsDhcU\njnaoS5vFVaUwvTl2fkbsfEHfdLgX8YEZaZ6bcW3E6PMnrX/ps8fc9RtKx8KFoqPfssgzeqPZ/D7P\nZvyEoAL3oZT790pJZ3szfm6W+Q7mdXM5EfoqBsDc8vG0wcki9/+WoFYb9esTf+Cihe+YYzGzOQdd\nHwX2SUh/VDs+PRL8mLzaoVml7QQLiUEoiBekDGz+B3BIxr1/FQvK7CZIqP2PZmi9Myl6eM8NYUHc\n2kIeszlxOBfYz4wfVkimiC+oQfpzyqhwlE6qkUSMv8voVkp09mfG7QWSPpOFW2/OZK2/CH0Vg2hh\nGh+Yyos/hhnPjhQOY5fC/2+LxE0aA/gUQQ35cmBpQvoym69xm81NVR01nUddGF8h+bOMZmxkdU+9\nlRqiGLHt12aFusGeWTGLtxPMER+JSNUXr4mulu3MeEnKjJR3RY6PTrk/z/ajyZ7aWWR6cx5XZ1yr\n6la6Cf4iJTw++SGJss9EEZ9if52YkfH9cAzAKcA6XRuQQrQwuKdA/DKbb7wB+ARBV8LDZtwq8WUW\n9u0/DGDGkaMAMy6RuJTkvu9DiAwom/ELiRcSDNyeZsbdEvsz7g74LILayxPJcH9gxnsBNJne6P0I\nunmOAr6YcD11cNWMaxivyVnKcdMssCmh5fFaghbcRjC2bWrVQdhfZ9Swd6UB8SZbFCfpw2v0230K\neIUZn5BYA1w+Fsm4suYz+lli4z5mmMRZZOzdbMZpEh8smEeZsaCZoq8tgyhFCvpXlUjvGwBmfNss\nWLJuxvPNeG0kzltIWZRlxn4kzFs3Y3l8dowZXzDjQTPuDs8vMRufy2zGRWZ8pMYK4EYx41IL9sr9\nL7O5NRijMRkVmGES5bOUX19Q5UXMnQpoxsfM2Dg8jT5Tn66QH2RPUbzerPLuYDD/nWV9FyemhD+X\nhhYUjghbzPsTWXcQPreNCpIZD5jxqSbTTMqm5fQHS1/FIN7Pl8fdJeLeS/a+CZjxz2ZcmxHl84Si\n0iJvIn/Ad1J8Eti27E1mvMdswRTTNLLEMOsFfiXj3T55NikyFkRYaP9pJMrDwLUps6EmxWi6ZdbU\n5cRrYddI404Hzfgm8HdUeA7iSTVgTttk2VhlGvYg6GU3kRlX0lJTzYzfEGyJWSeNpMVrjWI2tmVn\np4StlqYWM1XhV8AmKdf+14yHGupCW0awfiIqTKOUbwO2bySXHMz4T4kzCMaZymwN2yqhcOY9B1tS\nfJ1GX0ljj2xVAAAMdklEQVT7bl8M3Go2kdX0E6evLQNncjxAe66xo5SpEcZrXxMRIjPeYekby9fp\n9qnC2wnGtgaFBSviuxzczqNyd6wZn5tWIYAaYiDpfZJuknSdpM9J2ixy7SRJt0i6WVLqVEanNFX3\nok3FjE3NeHvT6dZktPbjyMxYAd9sIL8kH1edYsY7rb477hszrsW7VvtcgDfpL8wXnaVQp2VwEbCL\nmT2dYObESQCSdgaOINh34FDgw5K8BdIMbXupbJOiL1i0iT5qEWQN1tZehRv2hz+mgD1jt9XNt2W+\nDZkivyIhbDVwEPDCNgwKqfK9TcPmVL2nciFtZsvNbNTkugLYLjw+HDjbzB40s9sJBoP3rGWlA4yt\nYZg1ok3zplbPjtHj73Y7CKYWU2IzdjP2NeNzGVHi3R1mxmILXGp8IfGOeqzJj5KMGXdQzJsAVN/2\n9PyK900NTdXYX03gwwbgcYz7HllN/RkI00qVmQl9LbSaZDQAuZLAuV1fauGLIscTsSlcL3EKQYVq\nL+ZnTpWZQZfEu4B1a6ZRhqQ1K2U4jmKOJPPWjaR1hynleGbInE0kaTnjL8CIk83swjDO24Dfm9lZ\nGUklvjiSlkVOV5jZikxrF1L3heiSJdSwP8/lRQ/JKzyvJpwfb8Zqib3MWAPsKs1tOlIknaaZG3A0\nY03KrKUbCfaTaIXQrcqVAFLgUTO+pqVCmgZjs7DamgG0iOA7/CWwNynuJfIw42MN2TO4qaGSlpLg\n/aBpMsXAzA7Oui7pVcDzYczh150w74iKoJmb6IfGzJYVMTLdPi6R2JwBzoqqMTh4CoGv+WliKcFq\n1t8S9HVjNuZG+DiCFsJ7F9zZPvFC8gbgqUREKVyU9cZJGtUgTwI2YN4tdKOEgj5ilzbyiFFmQ54i\n4Z0TVpJXjM4lndJGPpXXGUg6lGCzmP3NLPogXQCcJen9BN1DrfoHN+M+ac4P0NTTp/UHJUl72XYF\nbgjdDkCCJ1gzfi7NLQJ8mKCpv1k83iQw42lSfwuOspiV3gWuz9xA4MerLlGXGEXdXAyeOovOPgSs\nByxX8BZfZmbHmNlKSecQ9Pc+BBxjZnVenrKrkZ1+krjxi1nujmgjrifYxe4RYPNpKpCdRvhD4Cc1\n3YCMiKYxM89ZZTEws1RvgGb2Lkq4CMjOh28xowM6U8Q61FjsA2DGXeRvZCIm9/LOTCExBMwyPbxG\nKVuW/LysLUOll+4onOnCGtjTNoV/iRw/aYLO/r4KXDahvOJ8k3Gvq04zpHnYHeSeFVVwMXCGyv5E\nFiNNsu/bqm8h2kTePwfe3FX+U0DZFt3MtABdDJxBEq4adpym6LM7jongYuA4zqyzAyVWd08rLgaO\nk889kLm/hTNgzHzGIrgYOE4RFjN8H/1OQNpWpWn4mMHAmJkfzJko/wb8ri9bkjr1MePhCe0lPjim\nQgzCzecPzI/pDIDG92yoymiPbMeZBaZCDADM+HrXNji1eTLww66NcGaeH0eOZ6bXYXAO3pzpxYxV\nZjzUtR3OzPPJrg3oAhcDx3GcCCl7YE89LgaO4ziOi4HjOE4GM9NKcDFwHMdxXAwcx3EcFwPHcWab\n/XKuezeR4zjOtGPGpV3b0BdcDBzHcRwXA8dxnAy8m8hxHMeZHVwMHMdxHBcDx3GcDGZm4xuZddMl\nJsnMzD2LO44zUaT5cQAzEsugUZy0613SVtnpLQPHcRzHxcBxHMdxMXAcx3FwMXAcx3GYom0vHcdx\nCnI5sHdOnNcBiyZgS2+o3TKQ9GZJj0jaIhJ2kqRbJN0s6ZC6eTiO40wSM/7NjHd0bcckqdUykLQY\nOBj4USRsZ+AIYGdgW+BrknY0s0fq5OU4juO0R92WwfuBE2JhhwNnm9mDZnY7waKNPWvm4ziO47RI\nZTGQdDiw2sy+F7v0OGB15Hw1QQvBcRzH6SmZ3USSlpM8iPI24CQgOh6QtSJuZjz/OY7jDJFMMTCz\ng5PCJT0V2B64ThLAdsDVkvYC7gQWR6JvF4YlpbMscrrCzFYUNdxxHGcWkLQUWNp6Pk34JpJ0G7C7\nmd0bDiCfRTBOsC3wNeBJFsvIfRM5jtMFEpcRTC39oRlP7NqesrRVdja1ziDi+MlWSjoHWAk8BBwT\nFwLHcZyuGaIQtIl7LXUcZ6YYtQz66JG0CO611HEcx2kNFwPHcRzHxcBxHMdxMXAcx3FwMXAcx3Fw\nMXAcx3FwMXAcx3FwMXAcx3FwMXAcx3FwMXAcx3FwMXAcx3FwMXAcx3FwMXAcx3FwMXAcx3FwMXAc\nx3FwMXAcx3Fobqczx3GcofB+YJ+ujegbvtOZ4zjOgPCdzhzHcZzWcDFwHMdxXAwcx3EcFwPHcRwH\nFwPHcRwHFwPHcRwHFwPHcRwHFwPHcRwHFwPHcRwHFwPHcRyHmmIg6U2SbpJ0g6RTI+EnSbpF0s2S\nDqlvpuM4jtMmlcVA0gHAYcCuZvZU4J/C8J2BI4CdgUOBD0uauhaIpKVd21AHt79b3P7uGLLtbVKn\nkH498G4zexDAzO4Jww8HzjazB83sduBWYM9aVvaTpV0bUJOlXRtQk6VdG1CTpV0bUJOlXRtQg6Vd\nG9BH6ojBDsBzJF0uaYWkPwzDHwesjsRbDWxbIx/HcRynZTL3M5C0HFiUcOlt4b2PNrO9Je0BnAM8\nISWpbvxkO47jOIWovJ+BpC8D7zGzS8LzW4G9gdcCmNl7wvCvAKeY2RWx+10gHMdxKtDGfgZ1djo7\nD3gucImkHYH1zOxnki4AzpL0foLuoR2A78Rv9o1tHMdx+kMdMfg48HFJ1wO/B/4CwMxWSjoHWAk8\nBBxjXW2n5jiO4xSis20vHcdxnP7Qyfx/SYeGC9JukfTWLmwI7fi4pDVh62YUtoWk5ZJWSbpI0uaR\na4mL6STtLun68NppkfD1JX0mDL9c0uMbtn+xpG9IujFc+HfskD6DpA0kXSHpWkkrJb17SPaH6a8t\n6RpJFw7Q9tslfS+0/zsDtH9zSecqWPi6UtJeQ7Ff0pPD7330d5+kYzu138wm+gesTbD2YAmwLnAt\nsNOk7Qht2Q94JnB9JOy9wAnh8VsJBskhWER3bWjzkvAzjFpW3wH2DI+/BBwaHh8DfDg8PgL4dMP2\nLwKeER5vAnwf2Glgn2Gj8P86wOXAswdm/98CnwIuGODzcxuwRSxsSPafCbw68vxsNiT7I59jLeCn\nwOIu7W/8gxX44PsAX4mcnwicOGk7IvkvYVwMbga2Do8XATeHxycBb43E+wrB7KltgJsi4S8D/jUS\nZ6/Iw3pPy5/lPOCgIX4GYCPgSmCXodgPbAd8DTgAuHBozw+BGGwZCxuE/QQF/w8Twgdhf8zmQ4Bv\ndW1/F91E2wJ3RM77tihtazNbEx6vAbYOj9MW08XD72T+88x9VjN7CLhP0hZtGC1pCUEr5woG9Bkk\nrSXp2tDOb5jZjQOy/wPA8cAjkbCh2A7B+p+vSbpK0tEDs3974B5JZ0j6rqR/l7TxgOyP8jLg7PC4\nM/u7EIPBjFhbIKm9t1fSJsB/A8eZ2S+j1/r+GczsETN7BkEt+zkKfF5Fr/fSfkkvAO42s2uAxGnS\nfbU9wr5m9kzgecAbJO0Xvdhz+9cBdiPoBtkN+BVBL8McPbcfAEnrAS8EPhu/Nmn7uxCDOwn6xkYs\nZlzZumaNpEUAkrYB7g7D43ZvR2D3neFxPHx0zx+Eaa0DbGZm9zZprKR1CYTgE2Z23hA/A4CZ3Qd8\nEdh9IPY/CzhM0m0EtbrnSvrEQGwHwMx+Gv6/B/g8gQ+xodi/GlhtZleG5+cSiMNdA7F/xPOAq23e\nt1tn338XYnAVsIOkJaEqHgFc0IEdaVwAHBUeH0XQDz8Kf5mk9SRtT7iYzszuAu4PZzIIeCVwfkJa\nfwZc3KShYX4fA1aa2QeH9hkkbTWaLSFpQ+Bg4Joh2G9mJ5vZYjPbnqCZ/3Uze+UQbAeQtJGkTcPj\njQn6ra8fiv1hvncoWPAKwVjZjcCFQ7A/wsuZ7yKK5zlZ+9sYECkwYPI8gpkvtwIndWFDaMfZwE8I\nFs3dAfwlsAXBoOAq4CJg80j8k0Obbwb+KBK+O8GLdCtweiR8fQKfTbcQzJRZ0rD9zybor76WoBC9\nhsBt+CA+A/A04Luh/d8Djg/DB2F/JI/9mZ9NNAjbCfrcrw3/bhi9h0OxP0z/6QSTDq4DPkcwqDwk\n+zcGfgZsGgnrzH5fdOY4juP4tpeO4ziOi4HjOI6Di4HjOI6Di4HjOI6Di4HjOI6Di4HjOI6Di4Hj\nOI6Di4HjOI4D/H+RMcklx9PqsgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x105b34e80>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(f0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "A = bundle_adjustment_sparsity(n_cameras, n_points, camera_indices, point_indices)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import time\n",
    "from scipy.optimize import least_squares"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "   Iteration     Total nfev        Cost      Cost reduction    Step norm     Optimality   \n",
      "       0              1         8.5091e+05                                    8.57e+06    \n",
      "       1              3         5.0985e+04      8.00e+05       1.46e+02       1.15e+06    \n",
      "       2              4         1.6077e+04      3.49e+04       2.59e+01       2.43e+05    \n",
      "       3              5         1.4163e+04      1.91e+03       2.86e+02       1.21e+05    \n",
      "       4              7         1.3695e+04      4.67e+02       1.32e+02       2.51e+04    \n",
      "       5              8         1.3481e+04      2.14e+02       2.24e+02       1.54e+04    \n",
      "       6              9         1.3436e+04      4.55e+01       3.17e+02       2.72e+04    \n",
      "       7             10         1.3422e+04      1.36e+01       6.82e+01       2.19e+03    \n",
      "       8             11         1.3418e+04      3.67e+00       1.25e+02       7.77e+03    \n",
      "       9             12         1.3414e+04      4.11e+00       2.64e+01       6.29e+02    \n",
      "      10             13         1.3412e+04      1.88e+00       7.44e+01       2.59e+03    \n",
      "      11             14         1.3410e+04      2.07e+00       1.77e+01       4.94e+02    \n",
      "      12             15         1.3409e+04      1.04e+00       3.98e+01       1.33e+03    \n",
      "`ftol` termination condition is satisfied.\n",
      "Function evaluations: 15, initial cost: 8.5091e+05, final cost 1.3409e+04, first-order optimality 1.33e+03.\n"
     ]
    }
   ],
   "source": [
    "t0 = time.time()\n",
    "res = least_squares(fun, x0, jac_sparsity=A, verbose=2, scaling='jac', ftol=1e-4,\n",
    "                    args=(n_cameras, n_points, camera_indices, point_indices, points_2d))\n",
    "t1 = time.time()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Optimization took 45 seconds\n"
     ]
    }
   ],
   "source": [
    "print(\"Optimization took {} seconds\".format(int(t1 - t0)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Setting `scaling='jac'` was done to automatically scale the variables and equalize their influence on the cost function (clearly the camera parameters and coordinates of the points are very different entities). This option turned out to be crucial for successfull bundle adjustment."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's plot residuals at the found solution:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x10b47a2b0>]"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYMAAAEACAYAAABRQBpkAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xm4LFV97vHvyywIMhjPYfSgggISRSIaJ44KCoag5CZB\nEmVQ0cSgJibKpOGIGtFcFePVJDeicomiqGEKgxwM+4ZERWQ8gMcjCobxABEZDITplz9q1d61a1d1\nVw+1u3vv9/M8++nqqlVVq3dX16/WqrVWKSIwM7PFbZ1RZ8DMzEbPwcDMzBwMzMzMwcDMzHAwMDMz\nHAzMzIwBg4Gk7SVdIul6SddJeneav6WklZLWSLpI0ubDya6ZmbVBg/QzkLQUWBoRV0t6MnAF8Abg\nCOCeiPiEpKOBLSLimKHk2MzMhm6gkkFE3BkRV6fpB4EfAdsCBwKnpmSnkgUIMzMbU0O7ZyBpGbAH\ncBmwJCLWpkVrgSXD2o+ZmQ3fUIJBqiL6FvCeiHiguCyyeiiPeWFmNsbWG3QDktYnCwSnRcRZafZa\nSUsj4k5JWwN3VaznAGFm1oeI0LC3OVAwkCTgFOCGiDi5sOgc4DDg4+n1rIrVW/lA80XSiohYMep8\n9Mv5Hy3nf3QmOe/Q3oX0oCWDlwJvAq6VdFWadyxwEnCGpLcCNwO/P+B+zMysRQMFg4j4N+rvO+wz\nyLbNzGz+uAdy/6ZGnYEBTY06AwOaGnUGBjQ16gwMaGrUGRjA1KgzMI4G6nQ20I6lmOR7BmZmo9DW\nudMlAzMzczAwMzMHAzMzw8HAzMxwMDAzMxwMzMwMBwMzM8PBwMzMcDAwMzMcDMzMDAcDMzPDwcDM\nFiGJdSSuGHU+xokHqjOzRUdifeCRCCbuHOSB6szMrDUOBmZm5mBgZmYOBmZmxhCCgaQvSloraVVh\n3gpJt0q6Kv3tN+h+zMysPcMoGXwJKJ/sA/hUROyR/i4cwn7MzKwlAweDiLgUuLdi0cQ12TIzW6za\nvGfwLknXSDpF0uYt7sfMzAa0Xkvb/VvgxDT9YeCTwFvLiSStKLydioiplvJjZjaRJC0Hlre+n2H0\nQJa0DDg3InZvusw9kM1sVNwDea5WqokkbV14exCwqi6tmZmN3sDVRJJOB/YGnirpFuAEYLmk55O1\nKroJeMeg+zEzs/Z4oDozW3RcTTSXeyCbmZmDgZmZORiYmRkOBmZmhoOBmZnhYGBmZjgYmJkZDgZm\nZoaDgZmZ4WBgZmY4GJjZ4jSacXjGmIOBmZk5GJiZmYOBmZnhYGBmZjgYmJkZDgZmZoaDgZktchJP\nk9h01PkYNQcDM1tUJP4QuLgway3wzYbrbizxcCsZG7GBg4GkL0paK2lVYd6WklZKWiPpIkmbD7of\nM7Mh+Udg79K8bRquuwWw4XCzMx6GUTL4ErBfad4xwMqI2Bn4TnpvZmZjauBgEBGXAveWZh8InJqm\nTwXeMOh+zMxapFFnYNTaumewJCLWpum1wJKW9mNmNp8W7JhG67W9g4gISZX/QEkrCm+nImKq7fyY\nmU0SScuB5W3vp61gsFbS0oi4U9LWwF1ViSJiRUv7NzPrRdNqonkvGaSL5Kn8vaQT2thPW9VE5wCH\npenDgLNa2o+ZmQ3BMJqWng58F3i2pFskHQGcBOwraQ3wqvTezGzS+Z5BnYg4pGbRPoNu28ysTRKb\n5ZMjzcgYcA9kM1vMthh1BsaFg4GZGewCIHWtLXE1kZnZQiaxNXA7NVVGEn8K3DKvmZpHDgZmtpid\nUZjerDZV5tPAbS3mZaRcTWRmi9lePabftpVcjAEHAzMzczAwMzMHAzMzw8HAzMxwMDAzMxwMzMwM\nBwMzM8PBwMzMcDAwMzMcDGxIJJ456jyYDWhRD2PtYGDDcqPE80edCTPrj4OBDdOTRp0Bs2GT2Ebi\nl6POR9scDGyYFnUx2xasnYCnjDoTbXMwMLOhk9h01HmoIrF01HkYV60+z0DSzcD9wOPAoxHR63Cx\nNlkW7FOgrGf3S7w0gu+OOiMl6486A+Oq7YfbBLA8In7R8n5sPLiaqEAigPUjeGzUeRmRXxt1Bqy5\n+agm8gnCFrN1R50BsybaDgYBXCzph5KObHlfZmbWp7ariV4aEXdI+jVgpaTVEXFpvlB6yVfgez9J\nb6ciYqrl/HQkcQ+wawR3jTIfRRLvBY6OYMmo82Jm80/ScmB52/tpNRhExB3p9W5JZ5I9b/TSmRTf\nvT+CFU23J7ER8P4ITuyQZnPgORF8v48sbwU8A8YnGAAvB5426kyY2QwJRcxPg4l0kTw1s2+d0MZ+\nWqsmkrSxpE3T9CbAa4BV5WQ9bnZ34EMSr+uQ5qPA93rc7jhzCx0bOonjJdYOYTv7SVxds9jH7gRp\ns2SwBDhTUr6fr0TERaU0vQaD/OA6r8O6g36mcbvh7R9UBxI7AbdE8PCo8zJhXsZwSpyvBZ43hO3M\nF/+earQWDCLiJug6Vs3bJZZE8IZOiST2ANYMLXOdjW0wkDgS2CqCk0aYn07m9X+XOhCtAf4KOH4+\n9z0pJJ5E1rz1/lHnZQL0cvyKBRZYxqEH8usbpLkSOIEu/3yJ3xpCfkYSDCT2ktisYlHxM58EfGye\nsjQULQ9ed0d6XfBDBZRJvFLioQZJzwNubjk7QyPxJsmNJUZhHIJBU+8DntUlzT8zTz0MJSSxSZpe\nkjoYDeIy4EMV83varjScjj4S20lD+V9elapyet3/FyS+Wpq3kcSrKpI3Po4ltkolzbElcaDUtRHD\nXsBGDTa3M7DF4LmaN6cBR6Wg8LX0nfd8npJYX+K/W8jf9C5a3PZITFIwAAYbM1/iN4dw0s4dDDyY\npmt/bBLXSOzXcJtV1Xa95veuIT1b4Bbg/d0SSVwh8eouyfqpjnwzcEhp3uHAd6qy0cN2P0tW0hxn\nLyf13k03aD81wLaGVpUhcZDEFcPaXqddAUeQ/cYeAi6ReGqP29gY2GDgjIinU2jJs5BNWjCoPbCl\n6R/4ERXLbpdYDvx6g300PbHs0HCdXwe+PugAWXkppPB+P4ltapJvUjO/V01+gC8gayk2MImvSrMD\nRyp1bZXe1h2vXQdFk4hUZTXwCWKevQf4swHW7zsYSGwisZPEOWnW68i+77b3Xf49vQJGdp+srjbC\nJYMx1qnovzXw0obbqf2SU5F1w/S2lwN9M2bqt3uV7+eLpfkXkO4fpBNd8aAdt+/1A1KjYRkOYW79\n/xqY7jMSMD2+/GcLaf6wYT6WMUE3/SSeBo1LlXUG+bz3k/3/fzvP0oB5aapqP6M6+T4xov3Ou7E5\naUj8msQnuiSbc2BLrCsN9QvrdNBdwUw1xQcartOLqjr6/DP/fsWyQwvT27eQn/zk+w8SWw+wnT+A\nngdry6/gN4PpUlX+v9gfOKqPfKxDzY9bYkOJe/vYZpsOG8I2BgkG5fPDW3tcv9/jcBtgx163le7j\nNb3oa2piLh4GNTbBgKyq4X0AEl+uqSMsNrP8pMQbyU4y83XVsCszJZBiy5/rh7T9d1TMa3owFm8c\nN/p/SLxY4ohCFVudt0Hljdue9tclL3/RINmgQb9Tc8BNgc0H3P7ISNws8faKRRNzMpP40zR5GH0E\nA+DFwL912P6ZQLmvUzd1/z9XE82Tw8i+2LJiVch76f3KqdeWOZtJnetIpdY73DTN8zMK002/1++R\nVT/tIbGLNKfZanHfyxpuc1rdfRKJHSR2K83+6w6bitJrv9bpsI1B6tb3lTh/iNvN19m9sI+tatLm\nng68smL+rABavvfUsiatnYo+3WFZk2O6W0OF1wG7NM8O4Gqi0SmcXNeReKrEnoXFe5aSd7qhvKzf\nLBSmPwZdW0/8bZ/7aar4GbdsuE7lVUsqRh8n8ZOKxW8BjumwzX2b7Dfdvzggva+7T7ISuK7B9soG\nDQYbkX7cxRvVEi+hewfJPK0qWqRdRFZ1NStdzfpPlti7NK+uXX3x+76nQfY2S9sLzTxprJzXBwft\n+yHxex2W/bHET9Pbr6eGG8MwjCvxQQJzmUsGbZB4T+FtPs7JicDngR92WHX/DstuStvuta57E4lr\n0/S2aRsfkfjdNG9jiR8V0peLs3mvz4FIrCdxOs1vju4m8Zx89Zo0f0I2dlNVC4mqQLNR8UQlcbLE\nyxrmZe/yzPSZNiNr+96Pnn/MEocUTswHFLZR/Bz/DlxcWm/9dOKfkrKWLOmm9aDPJ3gvswYdYzvg\nzjR9bZNmpKkEV3XV/fTCdF6Nmd/3eW5hWbdSBoVjqcoZHZa9mtml1CYXEU0MdPJNx8CGXZJVNaWe\nmGq2QY1FMABOrpj3PObeUO2ps1A68dye3i5vuNpSZornee/o44FvFNI8p5Q+319IvBn4r5r8nKDm\nHWhWAG9smBbgTTAdxOp+OJ2qtN5SMW8HZo9f8x6q72vk+yv+2KYq0q0B7uuQhzqRgkvH+wo1VSlf\nheng/EKaVzk9QlaK2JuZ4+AoOvR2llgp8Uxlndo27rL9XPHG/+4wq8/GM6h2DFSerKs+U14CKg4S\n2aTqo5gvGlSX5iMJt1Wt0kvHwtMqZh/YYNXfqZjnaqIRqDqBDdpTtHhy2IOaE4DEkRKXD7iv3Es6\nLFsBPF7sH5ACyCcL7/PSzpwSRwN58Kz7XrdrshFpuiXPAcBnetj/oV2Wz/pMEj+tqHKp85dUnwCL\n7lF1D+X82NqRmWaStSrudahmumwfsjbxV1LfYS+vpso/d6er1Vn12xLbF6q3VJif5/e5FdVTywrp\n8mV1v4PdOnwfV5TSlvfzwfRae/Isl9JTyWvbuvTl1RukyUvvVZ0gm/SZ6TYcTK/5mSjjFAyWVcx7\nesW8XlxSel938+i3gd9I04N+ya9tkObt0qxi/nsL0/mNyD8YIA9zPkOqF2/aZr1YzVV1U3JYngEg\ncewQt1lVLVg8zmdd2as0fEc6yR3eYfuHFdJVyTs21p3kyyfLZ6Xt5XkUWau1Kv/BTAe04v6L92b+\nJL1WXcnnw4IUW+VJYp/0tngfp9vvYN/0PyifLCs7dualdGlWdeSrgVu77Cd3iNT12SfvznfXcJtN\nOBiMQKf6/36Vr/D+DKavxkPiucqasBa/2F5bQJTtmPaxvlQbfE6ARoOMdaXqzlzHltKsR8NOd92q\nA1ryVw3SNP3x7S3x7NK8Tsd5+Z6UKtJvUyjN/Xl6rXtmRl6NdnTN8nIwyING8YElnQZczE+mddUe\neWe890lzWtfkx0rxBPcZYKXU8T5AlW9TalKdqkjLx/yu6eo/rx4sdp7cgd5UPtRF4uxSiUYMb4yy\numDQ5kXSSIxTMOhUvdKv4gm36ktdBXyJ2UMU5DcLa5+m1tChwA39rCjNrq+tWF4cu7+q+HtA6f25\n0LVDX+7/1cyfviEs8VhqNRKlZW1dLW1G8x/fkcBqifMLgbJTvsonpPWAd5bmbQrclqbzkseLqm6S\nw5zGA+V911WjNL0Yylt8bQQg1Q6X8HLg0ZplT2LmXsK70ut0CyGJU2FOQK1SPvFXHTtvYPbVf3GU\n4iel/f1fac6Dr2pJ7JEu5j6rbKTicmBcAtzddHsd9iNqAhDVVUoTbZyCQRuKnYjq6gzzp7Dlnpxe\nP1iRthdf6JZAqh2D/z+6rFqsgqgdyTX9YPaBOW36ywGlaDc636jdguwK8/OFeS/slpce/O/CPYs6\nc/63Et8uzdqfmY5yje6VJP8Njeux+xkvpy4w5f/D3WuW16lqJlznsvR6IUxXDVU5lJkg0c1v9rD/\nsrz566thVmunbvLv9SiykYqb6KdV0NOg9qmKriZagEZZ3PvIELbRrSnip6GypNHpxuUHOizrVIWR\nnxia9oeocjj0NfRw1WB5eW/TqoYIU9J0/Xq/dinfc+hE2eNaTyq8/zYtP4e8pOuAfgVNhx2v6hza\nVF7KqWs1Vaef51fMuSDqRNmw63d2SNK0tdjEmM8D0dpRWVyVpluz9HLFNSxdh76eZ3UtswYZDRSy\nk1K35w4UnVd6/xqyPg6LVb8PJeraIiynbMiapzC7kUYT3RqvnCLxWERtterEUcRo+lRIikXUn8MW\nr43I+oA8EypbTX2I+nrpherXyer0+xnJ97n013u9FRHzX10kKSJi6PttLRhI2o+sM9m6wBci4uOl\n5Q4Gthj8Ocz0I7GFxcGg20aldYEfk92kug24HDgkIn5USONgYGYTbSEFg7ZuIO8F3BgRN0fEo8DX\naPbgezMzG4G2gsG2ZM/Qzd1K8+Z6ZmY2z9pqTdSw/mdFYXo5zceSMzNbHCQtZx5Ojm0Fg9uY3bZ9\neyrHIFnR0u7NzBaGiJhi1rDnaqX1WVvVRD8EdpK0TNIGwMHAOS3ty8zMBtRKySAiHpN0FNlgVusC\npxRbEpmZ2XhxpzMzm2+3kI1U2/YjYwG+BfyvtjbupqVmZn2KYIcI/m6edndZ9yRzfBRY3SDdfASz\nebPQg8F3R52BeXBb9yQ962W8nSpjM1xA8vnuSaZVDUs9iGeTPS702m4JF4mruyfp6LXA53pI/1/0\nXh3+YeCILmneTfORXSfCQg8G/Yx9kvsv4B+HlZEWreyy/B/62Gbdk7aaumDA9YftX5smjKhMex3V\n4wo1cV8EfxPB8yg8MyD5BfDWPrc7TA/M476+VpjOR7ltvP8ILiIb3SBX+bzx0jqP97KPpOqhUUV3\npO0uGAs5GFzIYGOOb0r7NzV+OYQ6xx93WX4qcE+P2/xFgzQXAzuX5p2dXgd9iPiH0+sLO6bqrMmI\npE0fPvQQM0OFl7e7osu6xRNG+fe2cw95KK/35j7Wy72PVH0SgSIqR7590wDbr5IfE8P4TZ1emC4H\n2Dq9Xth1OzcuuBueCzYYRLA/cz9f1ZAY5Xn3pvX7OaGV1zm9MtWMYVxZdMrnFRH8O/BvPWzv3IhG\nB/q7IuY8WOXtDfLUVQR/mU5S5UdSdvKi0jZOLr6tWaeuuXP+VLgrK/J2cun9h9Jk3UNWiv+Ly9M2\n86vOX9F7ldwF6f9edxHxsQ7r7pleV1fst/hApadS/R1eUZjO/6dNn6BXNYz4g+m1pwuiiOmLm69E\ncD7Z91j83FXPtujlNxBkjzU9iLlPvcs91sP2JsKCDQZJ+SArv/9ABOcw+7msg0T8Y8ie2pR7sCbd\nNek1P6CuqEnXRNUPKX/ecf5Z7m2wnfwZBMXtnVqT9hTg5jS9BzNPRsvXrQsGVQ8S+oMGeevmXyP4\nQZp+RaG0dXLdCinY1FX9fLNiXtX/uVhF8ZmabU3/LyK4KYI9I3gi7f9hsoAAcFZhnTkXCYXPlA/z\nkr8vPrP6UzAdnCB73Gnum8yUEKuO8WWkYBHBf9Z8lkfS61HMHC/F/8F6pX0WRekVBns05S6QPZwo\ngtdHcBwzFyOXFtLlz3bu6eQdwWMRnMXMcV7W9AlrE2OhB4NOvhTBR9N08YCOmukmnmD2U7rqnkH7\ncGn57aXl6wNf72GfZeUf3rvIfuxFf116n5/Miie9P2LmGDl8euPB29KJjAiujpgeojlPWxkMImZ6\nURbmnc7gx2GxBDH9nUVMV+k06eOSP0ieCC4HfofsecqztllSHH8rr1orf2+NSkkRHFSVl5KHgO+n\n6XXSelcVlk9FzDr+LmemBPAwHapqUmmwWDpYU7H/W1PazwHnk534p7+7VIde3vbnoPYxpnkgXE3l\nCAWz3A0zT6aLYHUE95XSXJhenyA9PS1iTsDpeo+BBr/7hXa/AMY7GJTro/tRVzLYluxEV5VukJLB\n9aX3dcEgv5k1lV6LeSGCx2DOgV60bs309CaKrxH8KoKfl/bx/tL7/KpP6b0ieLiiyujMDvnKP2/T\n/+EVaV9N0tedIMuq7o881GC9bxXfRHBmBFeSnfS+RnXJID+5nsDMd39TKU0/j/AsW5bytHEEX0rz\nivnJhyf4aXqtei7xzyKmT7h1DSuKJ/aq0upDheU3RHBgYZ23V6QHeDyCR5m5AJpuXBDBg6nEs5ys\n0UL5SXBF34jo3CosYjo4P0EW2C4sLk6v3Z5gtmiNSzA4CaaL+QBU1EdXKdYT7l+xvPLzRXB74eRX\ndDWzr/Zyz2uQFyJm3bT+BFl1StF2ZEXpz5DdbD00za96OH1V3vMrs+LVZh4MjgS2TtP5s1vLJ9km\nz+vtVnSvu9LdqlCXW9uMMv34DwUOjeA3KpIcVXr/98BHInh3RdqzS+83jJjTPvziuryUVNZbR/Bb\nEXy6Zp2LU5oTC1fk5SBedZwV3UVWCgH4ILOfevZo2v7PyyuV8vtXZP//G1L671ekz0tHG5RKE3Xb\nbCq/eKhrtfZEWn43sCRibrPjdLHyQAQH9LH/yn1G8Hi6b1g2lJLBQjQuweAJut9M/Wrp/REU2tin\nE/FfltIUD+5vUS8/oPcA7i/M/2eyKojyFX/R71TNjODoCFZRGLAvgtvSQfrPEbyh5or4j9Nr0/rU\n/ETxnxFzgkC5NHAPWd108R5J0XbQ9SHxlT+UiKyaJJUoOv2vieC0CE4rzb4lLftcKe0fRfDBUto/\nTsveUEo758Qbwb5kwbFTM+Ne+wAsIws8VQEqmF3/31EEEZGVtiL4SAQnFhaX/0dFxav4x/L/fwf5\nSbmutApZIPtVh+Ud85GU81G8ZzJo/5Wu0vHXqWpuUZ7omxiXYPCPzD5xVz3wela9egRfJl0BF26u\nla+gitu8mforny8z0zHpX0hF/Qi+EcEL0wlcZAd2+WZsucVLeR9NOoVNrxMzPTNXlNJMVay3GvhK\nmi7+AB4kaxFyeHmFCA6KqG7bnoJVt5PBoM1G61zeIE1fj49M1RHb1Cx+PjS6Ii1+Rz+vKVlmi7P6\n/y16zGbHffa4rErXknaqFnxyj9stn0OOYnbLoXE5+Z5N9wdsvS69jkue59VYBIMIfgQzRfGI6Tv1\n+RXb4VS3DCnXz+Y/kL1K76HDFxzB5RHTLRNOjMhuPlXYiaz1TCfX0nub5jk/7HSyKZ4gyyfhpwN7\nF0oX+f0BRfCf6a9JXfnxwNoe8jrKH0pehddvHl5FqQdsBNcU6poHdR0paEfwyyFts06T3+70/YqK\nqqFhfY/nAf+/sJ9fRUy3wDkT+HjVSvMtgkciaznYqVVRr6WiBWUsggFABGdQ3wrl1IhZV9h5XfxV\nzG7KmVf3XF58P8Q8/izV3+5AffXQfREzHYLSyXpzsk5svXolNe2jI/iPUrG73x/3jREs7SF90/0c\nV3r/gS7pT6G+KWvZSuCfGqadFsElzG4x04uuJaIIdo+Yfe9rQIOUDPYE/maIefl21cwILo2offDK\n2TG3Nc9IdakmqxpGpe7+yoIzDsGgWLXzephVD3wcM23Yp0XwtvQaEfxLh20XP9/QrmgjuCWv522Y\n/r6I2j4HUB8Ef8Xsexgdd9M0PwNq1Domgo9RKHHETDPeuvTnR8yt1qpJ+7OIoY9E2fHkmm4Q9xI0\nG223i06Bq+NvN4IrIyobJvQlgv36Wa2Pdb5cs17VTfS2FJsn3wnsO4/7HplxCAbToxemG6tnF96f\nF9FfPXEyX8PLPpPuV76dDKPNcr/BoJf19qTZMA+5vAR3SA/r9KOXKpm+/9cRPVWnAby6y83MOnlV\nRqcb0fM+dHIfej4mIziCmftyeU/iC+jzfpE1Nw7BYJg3JDv1OI6K5YP6ArA2Xal2vPIdgm90Wd7W\njd1p6WqzSW/mPP3xabJTn4lB7Uhv9dJ/RGnoiqR8/2lgXUqtnXwYeAnDvYE8DL0+oGrQ0upUer0v\n5reTl28gj8gwv+TyD+R4ZppqDl0ER3ZoVdKLTvWYlwEPRPfx3xflARzBzb18BxHcVVWvn7YzFlfb\nEdwfwffofMKfj9/uzjAd0CHrO7NRD+sPdEx2qd+3IRuHYNBkhMzcv/ey4Qh+UDiJtlEyGIqI+o4w\nqXVT1aiSZeUhLWzyjbRkEFnHz7sL75+I6KlHdb/B4DJ6L4UMTTTrEb/gtPIM5B41bUEC8HI653ks\nT/bzYN0+66bnw2fJRoC03o26ZDCofk+qBzL7s99Zl9CGp5VgIGkF8DZmriqOjYgLK5Le2stJLEXs\nTkXHTj+ebutOgv8Dc/tAjHEgoKaXrjXTqTXRxcwenbPOKuA5w8lOz/oKBhGz+gIspbcGAoNYrBeT\nQHslgwA+FRGf6pLuoiHvt9uXeR7wsiHvc95EDL0Dz5fp4SlgNr9S66W6MZNWAa9osJnfovtTuzo5\nE9iqz3UHrm7powVXGxZFtVGbRc0mUXbYN3c7lgzSsBI93XeYRyd2TzJcERwxJj82a0kE/93pnlSD\n9e+J4KR+V+93vyOyqEsGbQaDd0m6RtIpkjZvcT8LRd1DQWzG12D20Ns21sa2+tLm6ruaSNJKqntk\nHg/8LTNXuh8m6zBSMTja+h+QHssPmKmImOo3P3m2OiybtKsUK0mlmPJDeWw87cPkVUGOZclA0nKo\nHfJjaPoOBhHRqIu2pC9Qc9Ub8Wh5yOlBVT2dKbdqyPsatr6L8mbjJoLvjDoPC0W6SJ7K30s6oTbx\nANpqTbR1ROTjxx/EPJ2I00BkVSOAjmXEL4rgBolnjjoftmj8Bb0167YFrq3WRB+X9HyyqpmbgHe0\ntJ8FJYKfjToPtjikDmVNnia4mIz9RWObWgkGEXFo91Q2YX486gyYWXvGoQeyjblJqGYzG4JFfZxP\nQpd2MzNrmYOBmVmmScngju5JJpODgZmZORiYmSW+Z2BmZoubg4GZmTkYmJklriYyM7PFzcHAzBaz\n6wrTLhmYmS1GEeze4yoLNmA4GJiZZRbsib4JBwMzM3MwMDMzBwMzs5yriczMFos+hmT/RSsZGTMO\nBmZmmcogEcFVwBbznJd552BgZtZFBL8cdR7a5mBgZovRoxXzfM+gH5J+T9L1kh6X9ILSsmMl/UTS\nakmvGTybZmZDFaPOwLgZ5BnIq4CDgL8vzpS0K3AwsCuwLXCxpJ0j4okB9mVmZi3qu2QQEasjYk3F\notcDp0fEoxFxM3AjsFe/+zEzmyeuJhqybYBbC+9vJSshmJmNizXAz0ediXHSsZpI0kpgacWi4yLi\n3B72U1k/J2lF4e1UREz1sE0zs369hLklgbEsGUhaDixvez8dg0FE7NvHNm8Dti+83y7Nq9r+ij62\nb2Y2kAhVlRsZAAAH1ElEQVQeGHUemkoXyVP5e0kntLGfYVUTFSPqOcAbJW0gaUdgJ+AHQ9qPmVlb\nxrJkMF8GaVp6kKRbgBcD50m6ACAibgDOAG4ALgDeGRFuxmVmNsY0qvO0pIiIRR2JzWz0pOl7ms8C\nbqwbuyiluxNY2sf4RkPT1rnTPZDNzMzBwMzMHAzMzABuJmv1eN6I8zEygwxHYWa2UPwigoeBA0ad\nkVFxycDMzBwMzMx6sGBbQDoYmJmZg4GZmTkYmJkZDgZmZoaDgZmZ4WBgZmY4GJiZGQ4GZmaGh6Mw\nMzuYbGyiRc3PMzAzayA9z2AtsMTPMzAzswXJwcDMzAZ6BvLvSbpe0uOSXlCYv0zSQ5KuSn+fH05W\nzcysLYPcQF4FHAT8fcWyGyNijwG2bWZm86jvYBARqwEk3wM2M5t0bd0z2DFVEU1JellL+zAzsyHp\nWDKQtBJYWrHouIg4t2a124HtI+LedC/hLEm7RcQDA+bVzMxa0jEYRMS+vW4wIh4BHknTV0r6KbAT\ncGU5raQVhbdTETHV6/7MzBYyScuB5a3vZ9BOZ5IuAf4iIq5I758K3BsRj0t6BvCvwHMj4pel9dzp\nzMwmhjud1ZB0kKRbgBcD50m6IC3aG7hG0lXAN4B3lAOBmZmNFw9HYWbWgEsGZma24DkYmJmZg4GZ\nmTkYmJkZDgZmZoafdGZm1tT+wPrAOaPOSBvctNTMrCGJnYA1blpqZmYLkoOBmZk5GJiZmYOBmVkv\nHhl1BtriYGBm1lAEPwd2HXU+2uDWRGZmE8SticzMrDUOBmZm5mBgZmYOBmZmhoOBmZnhYGBmZgwQ\nDCT9taQfSbpG0j9Jekph2bGSfiJptaTXDCerZmbWlkFKBhcBu0XE84A1wLEAknYFDibrmLEf8HlJ\nC64EImn5qPMwCOd/tJz/0ZnkvLep75N0RKyMiCfS28uA7dL064HTI+LRiLgZuBHYa6Bcjqflo87A\ngJaPOgMDWj7qDAxo+agzMKDlo87AAJaPOgPjaFhX7G8Bzk/T2wC3FpbdCmw7pP2YmVkLOj7pTNJK\nYGnFouMi4tyU5njgkYj4aodNjWbMCzMza2SgsYkkHQ4cCbw6Ih5O844BiIiT0vsLgRMi4rLSug4Q\nZmZ9aGNsor6DgaT9gE8Ce0fEPYX5uwJfJbtPsC1wMfCsGNWIeGZm1lXHaqIuPgtsAKyUBPC9iHhn\nRNwg6QzgBuAx4J0OBGZm421kQ1ibmdn4GEn7f0n7pQ5pP5F09CjykPLxRUlrJa0qzNtS0kpJayRd\nJGnzwrLKznSS9pS0Ki37TGH+hpK+nuZ/X9LTh5z/7SVdIul6SddJevckfQZJG0m6TNLVkm6Q9LFJ\nyn/a/rqSrpKUN6iYpLzfLOnalP8fTGD+N5f0TWWdX2+Q9KJJyb+kZ6f/e/53n6R3jzT/ETGvf8C6\nZH0PlgHrA1cDu8x3PlJeXg7sAawqzPsE8P40fTRwUpreNeV1/ZT3G5kpWf0A2CtNnw/sl6bfCXw+\nTR8MfG3I+V8KPD9NPxn4MbDLhH2GjdPresD3gZdNWP7fC3wFOGcCj5+bgC1L8yYp/6cCbykcP0+Z\npPwXPsc6wB3A9qPM/9A/WIMP/pvAhYX3xwDHzHc+CvtfxuxgsBpYkqaXAqvT9LHA0YV0FwIvBrYG\nflSY/0bg7wppXlQ4WO9u+bOcBewziZ8B2Bi4HNhtUvJP1tHyYuCVwLmTdvyQBYOtSvMmIv9kJ/6f\nVcyfiPyX8vwa4NJR538U1UTbArcU3o9bp7QlEbE2Ta8FlqTpus505fm3MfN5pj9rRDwG3CdpyzYy\nLWkZWSnnMiboM0haR9LVKZ+XRMT1E5T/TwPvA54ozJuUvEPW/+diST+UdOSE5X9H4G5JX5J0paR/\nkLTJBOW/6I3A6Wl6ZPkfRTCYmDvWkYXUsc+vpCcD3wLeExEPFJeN+2eIiCci4vlkV9mvkPTK0vKx\nzL+kA4C7IuIqoLLN97jmveClEbEHsD/wJ5JeXlw45vlfD3gBWTXIC4BfkdUyTBvz/AMgaQPgt4Fv\nlJfNd/5HEQxuI6sby23P7Mg2amslLQWQtDVwV5pfzvd2ZPm+jZlxmYrz83V2SNtaD3hKRPximJmV\ntD5ZIDgtIs6axM8AEBH3AecBe05I/l8CHCjpJrKruldJOm1C8g5ARNyRXu8GziTrGzQp+b8VuDUi\nLk/vv0kWHO6ckPzn9geuSN8BjPD/P4pg8ENgJ0nLUlQ8GDhnBPmocw5wWJo+jKwePp//RkkbSNoR\n2An4QUTcCdyfWjIIeDNwdsW2fhf4zjAzmvZ3CnBDRJw8aZ9B0lPz1hKSngTsC1w1CfmPiOMiYvuI\n2JGsmP8vEfHmScg7gKSNJW2apjchq7deNSn5T/u9RdLOadY+wPXAuZOQ/4JDmKkiKu9zfvPfxg2R\nBjdM9idr+XIjcOwo8pDycTpwO/AIWd3aEcCWZDcF15AN0715If1xKc+rgdcW5u9J9kO6EfibwvwN\ngTOAn5C1lFk25Py/jKy++mqyk+hVZMOGT8RnAHYHrkz5vxZ4X5o/Efkv7GNvZloTTUTeyercr05/\n1+W/w0nJf9r+88gaHVwD/BPZTeVJyv8mwD3ApoV5I8u/O52ZmZkfe2lmZg4GZmaGg4GZmeFgYGZm\nOBiYmRkOBmZmhoOBmZnhYGBmZsD/AK67rdaovoU2AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10703b1d0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(res.fun)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see much better picture of residuals now, with the mean being very close to zero. There are some spikes left. It can be explained by outliers in the data, or, possibly, the algorithm found a local minimum (very good one though) or didn't converged enough. Note that the algorithm worked with Jacobian finite difference aproximate, which can potentially block the progress near the minimum because of insufficient accuracy (but again, computing exact Jacobian for this problem is quite difficult)."
   ]
  }
 ],
 "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.4.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}