ipython/bundle_adjustment.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Large-scale bundle adjustment in scipy\n",
    "\n",
    "TAGS: Optimization and fitting\n",
    "\n",
    "AUTHORS: Nikolay Mayorov"
   ]
  },
  {
   "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": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from __future__ import print_function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import urllib\n",
    "import bz2\n",
    "import os\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First download the data file:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "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": true
   },
   "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": true
   },
   "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": {},
   "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",
    "    \"\"\"Compute 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": true
   },
   "outputs": [],
   "source": [
    "x0 = np.hstack((camera_params.ravel(), points_3d.ravel()))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": true
   },
   "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 0x1067696d8>]"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYYAAAEACAYAAAC3adEgAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztnXm0FcW1xr/NJaiAIAYFBVEUUHBCFCQS9SqOmDjGIdFE\nxfiCs1GJoK4HJnlOcTaOMRo1cVYUjZFBvMZZCDgxieKAVwERQVQkgPX+qNOe7r49z33O91vrrNOn\nurtqd5/u2lW7du0SpRQIIYQQg1Z5C0AIIaRYUDEQQgixQMVACCHEAhUDIYQQC1QMhBBCLFAxEEII\nsZCIYhCRjiLykIjMFpGZIrKriHQSkYkiMldEJohIxyTKIoQQki5J9RiuA/CUUqovgB0BzAEwCsBk\npdTWAKYAGJ1QWYQQQlJE4k5wE5EOAGYopbaypc8BsKdSapGIdAXQpJTaJlZhhBBCUieJHkNPAEtE\n5E4RmS4it4lIWwBdlFKLAEAptRDAxgmURQghJGWSUAytAQwAcKNSagCAr6HNSPauCGNvEEJICWid\nQB4fA1iglJpW+f0ItGJYJCJdTKakxU4niwgVBiGEREApJWnkG7vHUDEXLRCRPpWkoQBmAhgP4IRK\n2vEAHvfIo7SfMWPG5C4D5c9fjnqUv8yy14L8aZJEjwEAzgTwDxH5AYD5AE4E0ADgQREZDuBDAEcl\nVBYhhJAUSUQxKKXeADDQYdc+SeRPCCEkOzjzOSaNjY15ixALyp8vZZa/zLID5Zc/TWLPY4gtgIjK\nWwZCCCkbIgJV1MFnQgghtQUVAyGEEAtUDIQQQixQMRBCCLFAxUAIIcQCFQMhhBALVAyEEEIsUDEQ\nQgixQMVACCHEAhUDIYQQC1QMhBBCLFAxEEIIsUDFQAghxAIVAyGEEAtUDIQQQixQMRBCCLFAxUAI\nIcQCFQMhhBALVAyEEEIsUDEQQgixQMVACsfRRwNff523FITUL6KUylcAEZW3DKRYiABvvglsv33e\nkhBSXEQESilJI2/2GAghhFioKcWwdi0wdWreUpA4GJ1HdiIJyY+aUgzjxwODBuUtBYnD4sV5S0AI\nqSnF8N//5i0ByYrvvgNefjlvKQipTRJTDCLSSkSmi8j4yu9OIjJRROaKyAQR6ZhUWYRMngzstlve\nUhBSmyTZYzgLwCzT71EAJiultgYwBcDoBMsiJWXffavjB7NmaQ8kM0HHGFavTl42QogmEcUgIt0B\nDANwuyn5EAB3VbbvAnBoEmV5y5F2CSQukydXt+fNy08OQog7SfUYrgEwEoC5nddFKbUIAJRSCwFs\nnFBZpIYJqtzZCCAkPWIrBhE5CMAipdTrALxe10gOiCLAkiWRRPPlq69YweSBl5mI7qqE5E/rBPIY\nAuBgERkGYD0A64vIPQAWikgXpdQiEekKwNURcezYsd9vNzY2orGx0bJ/xQqgc2d/QcJW8itWhDue\nJMvEiXlLEIyvvwbatctbClLvNDU1oampKZvClFKJfQDsCWB8ZfsKAOdXts8HcJnLOcoLQKn58z0P\n+Z4HH9THB6W5OdzxSXLKKUrtsks+ZefF2rX6fq9dq3/369fy/n/6qU6bMaOaNn++UldcYT3uqaey\n+e/efrtazlFHKfXCC+mXSUgQKnVnonW48UlzHsNlAPYVkbkAhlZ+R6IWzQoTJgDTpuUtRbY895z+\nNv7PhoaWxxi9PvN/ftttwO9+l65sbixdWt1+8EH9IaTWScKU9D1KqecAPFfZXgpgnyTz94PjBcXm\n9tutv50UQy02AggpG6WY+ZxWZcFKKFvuvdf626wY3norXF5ZNQLc5lkQUsuUQjGQ2sLJlLTDDsCC\nBdZjlAI++IChTgjJmppSDDQllYvWNkPmmjXW3888A/TsaVUYBnn91+wxkHqgFIqhll/G+fOLf31H\nHQXssUf1tzGIHBXjelu1ck43uKsyb948AGznm2/iyeIHGxukHimFYkiLIlTIW20FjBuXtxTeTJgA\nPP+89sgRAWzTTALhpEy8Kt2VK4HmZv/j2rWrxk16+WXgT38KLxshxEpNKYaytu7KMtHOHOb6hBPC\nnfvBB/7HXHop8PDDeru5ufp/muMriQAXX2z9r7/7Tn//8Y/Ju7UWafB5zRquhU2yoaYUQ1kpQs8l\nCOZK0jDzRMG4Xnule/vtwBln+J9vmiifOf/+d35ljxwJtG+fX/mkfqgJxbB4cbTwCmWpkLPk22+B\nd9/V2++8A4wYkU9PTKng5ab5P9plCOtWmyRz5uRXNqkvakIxXHghsP/+1Zf4yCPzlScsRVJQl1wC\n9O6ttx96CLj11uq+vE11kyZVt82yrLde9rIQUsuUQjGErTgNOzUJz/Ll6ZfhZkqyY98/fHi445Mg\nTJ5NTUDXrsnLoBTwz38WqwFBaptEQ2KkyddfAwsXai+eWqOIL/zNN1cHxY3KMe8eQ9F58UVg0aLk\n8zXcejfZJPm8CXGiND2G884DevVy3s+KKx7vvQd8+qk17dRTdSs1L5yUpTmtHv9r+39ESFqUQjEA\nwLJlyedZlJZ63nL06qXXYrZjXzQnTmVsPjfI9YYZfHYqIynC5FmPyorUJqVRDF7whYyPk398kgor\nSl5e/2tW/znjNJF6pBSKIWjrsawKIu8egxv2nkKW97co9+Sll/KWgJDsKYVi8KOsCqFI+Nn0gWxN\nSXHLy4OyyUuIG6VRDGm8dEVplZLwOD0PeY8xEFIrlEYxeFF2r6QiKCine2fMtE3ivpqvMergs/m8\nefPiy5Q0ZX3+CLFTGsXAly5diqCczPjJM2KE9bd9LYekWbu2ul20e0VI0pRCMfi9iFGVRlFe8KLI\nkSZRTD9PPx2ujFdfDXd8EAwZb7yxmmZfR8J+bNqsWhU+ui0hYSiFYgBqzyvJLGtZgqMlpcCCmpLC\n8uWX4c8JShozmqPS3Kyj2x58sHvI9hkzgA8/zFYuUjuURjF4USaFYGCu+KZNy0+OWiHKhLggFHGC\nmzHZ84kn9BoWTgwYABx0UDbykNqjFIohrZeetMQIuZ0FSf+naT4jRXr+dt65uv3ZZ+7HGQsYkeiI\nOEcFqHVKoRj8iOqVVBTbflHkAICnnmqZlnSlGOR6nSo1r/PqqcdgJq14ViLAK6+kk3fZMK8gWC+U\nRjGEfemmTElHDie6dcu2vDQoknKKSr30GMx49Qri/qd5LkpE8qUUiiGoV5J5sG3oUH3ee+8lL4+I\nddDvk090yOWoKKU9TdIcPE2COBVNEhPSvI6vV3Nj0oPi48YxPhQpiWIAgr30zz9v/f3AA+6huuMS\ndFH2jTYK5h0yfDjQsWM8mYpMlAluQdLMOD0jK1cC33zjX55fnmX0iouiyA8/HNh66+RlIeWiNIrB\ny3PH7YUcM8Y5XSngtNPiyxSEJUuAWbO8j1GqOuj7xRf6nNat9aStZcuAzz9PX04vglR4Cxbk2+Nx\n6zEMHQpsv3328pSZDz7IWwKSN6VRDPPnhz/HrcW0Zg1w003xTCNhznU61q2y3XVXoG9fPdN21Sqg\nUyegc+dq937oUGDu3PDyJoGXgujRI/ykqyy8kt56K9qzY8/TPPMZ0Lb9jTe2pi1cGL2cNAj7fD/6\naLzza4U33gC++ipvKfIltmIQke4iMkVEZorIWyJyZiW9k4hMFJG5IjJBRGIZSuwvplUG67eB24Od\nxAMfVzG47Z83T/cY7Olnn62/p0wBnn1Wh4OeOTO4DEngdx2G3E5EWagnTPlOPYbPP6++4C+84F+m\nF9dcY/29dm1LV9FbbolXRtIEfUYPPVTfuyOOSFeestC/P3D88da0elOSSfQY1gA4Rym1LYAfAThN\nRLYBMArAZKXU1gCmABgdtYBly6INiPn9matWWX8vWBC+jKgECWltPsa8rONZZwFDhgDDhqUjm5cc\nSRwXFb8Wuf0emiu6++6LVqaR58qVwY8tG48/Hv3cFSuynS+xzTbASSelX86jj1qvyx6bq9aJrRiU\nUguVUq9Xtr8CMBtAdwCHALircthdAA6NWoZfgDS3F9KvxzB4sDW9R4/g5oCkewxOx5iv27ydtdfI\n0qXJ5peGAnHqMSxfnl6ZZWhBxpXxiy/8j+nQoWVvyokVK3RDLIiCdeL004HXX9dm1DvuiJaHH8OH\nAwccUP39wAPV7dtuS6fMopLoGIOIbAGgP4BXAHRRSi0CtPIAsLH7mUmVb/3t92I4DZYGrXSzqBjM\nLZYnn2x5fR99VCwX17jmmrjY74/5dxaylbXH4Ma//hXsuCCD1R06AOuuC7RtG37S4PXX6yCGaQ+K\njxsHTJhQ/e22zvzq1bX3X9tpnVRGItIewMMAzlJKfSUi9qrTtSodO3bs99uNjY1obGy0nhixEn7/\n/ZZpX32l/1g35szRPYck8ZN/+XL94tix2zmd6NgxuzGTsOWMHg2ceCLQp0868tiP93pZnZ6FIJS5\nAghyD4P0CvJmxgz97efdB+h3qaEBaN8+fDl207Ib334bPu8kaGpqQlNTUyZlJaIYRKQ1tFK4Ryll\nWCwXiUgXpdQiEekKYLHb+WbFEARjDkG7dkb5wc/t00fPVHbjb38D9tvPPx+3l66pCdhzz3AyzZ4N\nDBzYMj1ouIM1a7R7a1z23jt+HmYuu0y/pH/8ozU9K1OS+beX80LU8moBL4eBvPj2W/0e2u36ft5l\nt98OnHyyjiUVNDDlnDnaC1Cp6GaurLA3mi+++OLUykrKlHQHgFlKqetMaeMBnFDZPh5AjCEuKwMH\nAoMGtUz3q4yV0oO4s2d7H2NwyinuJgh7xWCsKLbXXi3zD+KuOnWqu0x+JDXm8OyzyeRjJog3WVJ4\nKYaoC/kUPVaSF3kqLxHvnuKbb7rve/FF/e6F5eST9XcYd+4oKwEWXYEkQRLuqkMAHAtgbxGZISLT\nReQAAJcD2FdE5gIYCuCyqGXYH/DZs63dymuv1d9+5oJLLvEvy2zXv+UW4M47g8k4fbpzHkDttC6j\nENVjJWlTUtI9hjIQ5B56mUXefjtYOX/+s3O6Uek6yRFlUaWgijfMf+1lPjIvzmQmqlmyTMQ2QCil\nXgTQ4LJ7nyh5LlyobYVhp+a//LJzulFpXHRR9bcbQSuyZcuAu+8G3nlH/545sxqSQylg5MiqnfPQ\nQ/XCKkcfDayzjo6tlGQMp6wUT5RynFrqaclrrzjML31UBVXmHkMQjOfXiTRn3Ce5FKtS1p79ypW6\n8u7Z0//cI49035f1PKEiUciZz4ccov2Vk8I+i9Erdo69ArEPVBuum6tXA3/5C3DPPdV9xx5b3b7y\nSsA8dHL88dorA7C6wSVBkXskV1/tvs+oSJMaPLdXzEFbvLVKkPv6ySfpy+GEl6K2/49/+5v+dmv4\nvf8+sMce1jT7bxKOQioGe8WdZcVnf2Dtspi7wPYH2Jgg5yfvOedEk82NrO5PnBbx73/vvu8//2mZ\nVmRlBzhPeCpjjyEpd+elS91b2FECIjrR3Oyc7jTGFmTczXztYdyZy/g/h6VQimHpUuCRR1re+EmT\nrL/D/jHffBP8QfzuO+3PfN113sd52bSLXqm99BIwfny2ZZpXhrPfn7CL8jgxZ046rV+3/9hoxRaZ\nLJ/Dk08GtttOL+5jd4F1WvAnimxO/8U332jvNzuLTT6Qbi65/ftXt5OexFl2EpvHkASDBmnb+w47\nWNPt69oaD9XkycG6jF27atNOEJqbdehhA+NhbG62+kYfcYR2xXQia8XgVN7ixdrGvtlm+ve77+oK\nuE8f4Oc/15PjJk4Mt2xh0OsaP15PZEo6Xy/sAeBuvTV+nmExV1yjR+sWqdsAZi1gXmLU8NT50Y9a\nehQ5ebt5/eduytgp/dRT9fidFxtuqFcmPPBAa3o9zm0JSmF6DAsWVAdkvVzZzOy7b7BJYABw3nne\n+w37v92sYZiSune3xt756KPieCc4vWR77VWdqHfooUDv3nowX0TLDgBXXZWOPIcc4h2QzS5vEqaG\nyy+vbq9a5f9/ByVqJXDDDTqCb56kMWnRjNkTz2vOiHlRqyTKNRN0HsZi11lUGvYYrBRCMRxxBLDj\njtHOTSoEtZvbmnkN5KCrZRWhx2B253ULkmY+L0gAwTAVjZ/tes89qybCpO/XHXfkEzbZXDkWoVUZ\n5L4mJac5H7PZEABapVjLuPXaAR1bKShJj/uVnUIohkcfjT4135gunwVBvVyKPsYQFTevEDNPPhks\nr3//u7rtFOogzj187bXo59qJWnGmWRkmSe/eyeRjvk92Tz6/yMF2wsxD8Hr/zY0dv+cpjPtsEZR+\n2hT68Y2zJGMY7PZpO889V90u4kORlSIK4ukRpEtul9dpJmmca0pyYNhvvsnIkTrO1YknVtM++aR6\nTVmGcrcT5B4aJtSw2F2uvRoEYZWkMe5ln5Xs9GwFvb9O5iwzYZ63JN2si0qhFcO992ZTjt8CJeaZ\nnUEGse1eVGmz8cb+A3BFpkgvmH19bq95GIB+Hlas0MrIqDC6dau2ms2myFrimGPc99nH6cI0psxj\nBlGCL5oxP1d+Pdl6X7HNTqEVQ5YLgATlmWf8jxk1Kn05zKxZ09JzKyu+/VZ7hsShKIphyRJgiy2S\nzTOjYJiOBLmvW22VfLn2nn4YU1Jaz0KS+Rp5FeW5TQMqhhph7lxgyy11b+X224OdE9csJgKstx5w\n883hzrPboJ2U2t13R5crLP366XAmZrluvNE5FLoXTvfz/vvjyZY2YaLyZhGBNElTrXk9BfYIwlFo\nxVCUwGcPP5xcXmmOUbz/PvCLXwQbJAb0i243nWSBfUKS06p5Qa8hCWbPbtmyP/10f7t00QnSog1j\ngnzkkXDzU7zkcJMtyffD7MqexvPEHkNOxFmLtl5ZsiT40oezZoVrSZn91t145BH/Y8zrVxeJuJVS\n0dY2CFJxhTFBGo0Ic0ywqGShGNImT8eCtCm0Yih7i63orFyZ/Iv4xBMt0+zKp0wvP2lJWKeQvHoM\naRMkemtZKbRiqOWuWhFIaoEfP+wun0UxEZoRKVelFISivD9OcriNH5bhPyjKfU2TQisGp8U8/OYc\nkOBktai5vYwg5qY8qIcXPg5Rn5WkoqsWBb8lRmuBQisGJ/zmHJDg5KUYikQtL9OYZUiMsORtSgoT\nZtuO1xyOWqF0ioGEw2sVrqxMOkVWDLfdpr9r0ZQUhCyuOcqs4rT5+ONsyikrVAw1jl+IinrsMSxd\nWu0pOCnHJJedzJMi9BjOPz+459PKlbVz78tOodZjIMnj9+JnUWlvt136ZYThhz/U61QMHw6sv75O\nM/cYatkN0U6rVun2HK+4Ivixm2wCDBmSniwkOOwx1Dh+FX+9ttAWLAAuvtialtaAaF49piBBKMPI\n9tln0WWx43Svly9Pd53uiRN1uc3N0WZCT5sG/OxnyctVRNhjqHH8Xvw0X0QSnzj/T5DKL4himDVL\nV6ZOS2iWif3317Psu3ePdv7AgcnKU2SoGGqcIpiSiozRci3q4PP226ebf5Br3nbb6hKxSZGVu6pT\nuBXiD01JdU4RK8MscaqgJk7MXg4nmpvzlqBK0NULi0aScc7qCSqGGsev4q+ngVY/5szR37fckq8c\nBgcfHOw4r7GEW2/VoWVeeMF5GVy3JW0NDH//pGfJZ9VjKMtKekWDpqQax08x1OpCMm7Yo6gadngR\nvQ41UIxgeJ9/HixoIQC0a2etaM3BJ0eMANq3B447Tv+2V8h+FfQbbwSTISzmcv/732DrnETBa01o\n4g71aY3j9+L7tRhrjb32sv6+8caWx2Qx+em///UeHI4iQ8+ewLhxLVcrM6+oZszfCFoRP/10eDmC\ncvXVWrbHHgOGDUunDPYYosHbVuOccor3/noaY3C6VmNmeNb34eST9RwKYxnYadOqrVsR4L77wuf5\nwQfA4Ye3XKjJPD6w4YZ6OdJ99gmW55Qp4eUIwooVwLnnak+n//3favpHHwU7f/Vq50i+dsaMsf7u\n2jW4jEH49ttk8ysMSqlUPwAOADAHwDsAznfYr3S7lh9+8vuce276Zfz2t9XtDTZwPuYvf9Hf++8f\nLu9ly8LL8/LL+d/3qJ/ddstfBuO+54WuvtOpt0Xnnw4i0qqiEIYC+ATAVADHKKXmmI5RQHoyEEJI\nmowfD/z0p9mXKyJQSqXS103blDQIwDyl1IdKqdUA7gdwSMplEkJIZgT1HisTaSuGbgDMDpEfV9II\nIYQUlIK4q441bTdWPoQQQgyamprQZPe3Tom0xxgGAxirlDqg8nsU9IDJ5aZjOMZACCk1KVajrpR5\njGEqgF4isrmItAFwDIDxXie0LkgfhhBCgjBgQN4SJE+qikEptRbA6QAmApgJ4H6l1GyvcxYvrm6f\ncUaa0hGSL9tuq7+5jnl5ee016wTCWiH1CW5KqaeVUlsrpXorpTwD9yoFdOqUtkSEZM9dd7VMe+st\n4OWXgcMOqwZ7mzcvWv4nneS+79xzo+VJ/OnXL28J0qHQhps87HaEJM1HH+mw1fPnA19+qUNBGAwe\nrL+POELvW3993Qrt1Ano3Tt4GddfD/z1r877+B6lx7rr5i1BOhRaMQShVSvgu+/yloIQd9q1099j\nx3ofZywzOnBg+NXS2rbVYTR+/vPQ4pEY1GqQvsLESurQoWVakJZOWePEk/ohqzhMxxzjvb9/f+Ci\ni8KP3V1/fXSZAGDo0Hjnk+wphGIYM8a5gvdTDEWJm0+IF3kHKlQK2GQTYNAg4A9/0EH0Hnoo3Plx\nOPDAeOeT7CmEKal162i2ut/8JtkFyglJgyiKIU5l/PjjwCG2wDMLFlRDULdpE25R+7im2jZt4p1P\nsqcQPQY3OGhGaoEsewyrVzvH7mloiC6H8R4aCxmFZZ11op2XNEcfDVxySd5SlINCKoZttkk2v/PO\nSzY/QsKQpWJwmiAat4FlnD9kSLx88ub++4HRo6Of/6tfJSdL0SmEYhgxwjnd6YHu2TNc3u+9pwfd\nDLp3D3c+IXFp2zZvCeKxySb6O++xkrzZaKO8JciOQiiGzp2d050Uw847+x9jZsstq9tHHQVcfHE4\n2QiJS94ujXF7DBtuSLNuvVEIxWBn222d3VcB4O67o+fbvn1x7J1R2W679PI2T7wixCApU1TWuA16\nb7pptPzMPaaklwgtGoVUDPfeC3z6qfMDtd56wfM5/njrbxHgJz+JJ1vepGnnzbtlS4qJ8R66VfA/\n/nF2soRhs82c05ubs5WjjBRSMbRpk4xd1v5gtGsHdOwYP988ufnmvCWoH3bdNW8JkiGpFr9bPkYw\nQDd69IhXflTat8+n3FqgkIrBIMgD7XWM3UPjN7+JJ08RSHsAsFevdPPPCre4QWHwq/AMdtghfll2\nitSA8XsPjz3We3+YmE9JcvbZyeZndmKpdUqpGD78UAcmA7zNH0bsGYMwZqh6RCngzTfzliIZknAt\nDKok+/aNX5adIq1L4tdjKKq30gknJJvfscfWzyB8oRWDGz16VM1Ebh5NQEtF0MrjaoO2DovI2rXJ\n5RV1nOHpp53T87qvSVSsbg4QdooexDHtweO4iuHvf3dOf/ZZ4NVX4+V97bXxzq9XCq0Ygj7QBxzg\nnB7GA8mIgFkk/vSnYMd5KTwvvvsOOOWUaOfa2X9/4NJL/Y/zMzskwcSJyee5erX7vvXXB8Z7rkuY\nLh9/nG7+xnuY1nwMN7NZY6OONPv889HzjqK0goyJ1HrPodCKIShO7mejRgHHHed/brdu+ruI3eEw\nM7Y//hh44QUdJC0oIsANN1hbbHHuw0EHOZdhJkyMnqi4xd366U+j5+nVA1l/fWDHHaPn7USYisd4\nhtPCqLjd7qvfM+N1LeefD+y+e/X3wIEt847j9bTxxuHPmT49enm1QqEVQxyt/OMfu/sxmx9EA/PD\nfdVVLR9QN7bZBvjBD8LL54e9leQU/8ZMt27alfWii4Dly52PcbofDQ3A3nvr7bjK0Zgha+ZHP/I+\nZ9w45xZhHI8go6doN4vtu2+w88OOGShlNWkOGxbu/LQw5rzEeY8++cQ/RpL5uTnnHOs+vzUoLrvM\n2mNIep6R3/PnxA9/mKwMZaRmFYNXt/fxx1umeT3cXiiVThd7ww2tv2+6Kfi5brbxKVOcB+SMCj3s\n/W7XzloJdu5szWOzzYBbb7WeY1c+hx7q7O1x0UUtTQxBzX1BlbobYZ0U7M+AU8MjLEUxVTgpezsi\n+j9saNCNKjNjxiRXjhlzRIM+fcKdG4Si3P+8KLRiiMrcudo+6UaQdaWNB9rPHKVU9SEyWt5xmTrV\n2mrdZ594D2qvXsDIkbpHceed0fKwL+6ydKnumXh5fvzhD869ELtN3OmYhobwlYVXfl4ccUS0cgyM\n/8YwaY4alV3Fsv32/sckJYtXPq++CixbpreDNJTuvFM/1wZPPVXdPvJI//Pfew/4z3/09muvubsn\nB3UgIFYKrRiMBzFIRW6mTx9r5bDTTv7n2GdEd+miv0WqD7wbRlfb64VwmkJ/4onOx+6yS8vKLezL\nPW2a7hktXKgXmL/iiuo+L1NH69YtTQfTp+tVvPbbr5rWqZOuvI880t0rx8323a0b8OST3vI7XW/Y\nCt/eenVyYVUKePhh93KCBE4zrv9//qflvpde8j/fiaAOBXG8vt57L/q5ZkS0mdKYUBbEI+xnPwMm\nTar+7t9f9x6PO875f77oopZpAwYA77+ve5ZbbeVcTqdO0RsY9UyhFYPB0qXxzu/XT3saeLm2ej08\nRqvjppv0gu5mlAIeeEBvmweLV6zwr8zvuMN7fxx23lmPSxgKzsxf/+ruBihSvR4Dvxfd/iIbx5tb\nhPZjDzoovLILe7wxLnPQQXqhmiiTxn79a/9jtthCf0f1DnMiqNtwnIifbiEj3LDPCzJIynFj3Dj3\nSaiDBjmnG/fei1tucXaMIO4UWjEk2R3/8ENnG7XxUB92mHO68X3ggdq1s2dPYPJkq4zrrae/zS1t\no/VkmG6MFs0FF7RUULvsUt12c/ncdNNg7qBB6NrV/UUDqhVc1PkRjz2mP0Fxq1ji/v89e2qz1ZNP\nOodbf/ZZf3n8Kr333wfOPdd9f9oT1ezRhp1wu49hnSbsHj5uc1eCmCvDKhM/x4HBg4H/+z/nfQcf\n7N9DJVZwZ2bJAAAPKElEQVQKNL+yJeYHeqedgBkzksu7sRFoatJeER9/DGywgW59GWG5jRdaRH/M\nNtChQ/VLsnhxy5fu7rutnhDHHacr4t13177wG2ygbf5LllSPOfNMva9zZ3eXyoYGbbvOgo020p5C\nrVoBf/5z+IWTvFpncVqXUc61m7NGjKiuFe41DhWUDh3cewpTpwaruOMQZUncqBjvhPHs778/8Oij\nLa/x8MPjlXPllS3Hrvyuc511dMPtwgvjlU00he4xmAkSbyWMD78xqHzssdqXGtAPuzHha4cd9GpP\nbnkuWqS/7Yrhl7+0hlJo3VpPwGvXTisFQI8tjBxpPW/4cN2yKcp8CsN3/LTTqi1Lo9UdtiV/zDHB\nF4Q33BWj9BaCxMa5+WbdyrdjHhAP8x94Hes0VhQFr3EEt4bEK6+E67UF4Sc/0XNlevSojqcddpi/\n2cvN/u/GZpsFf16i4mVWJgVXDIMHh3MdDPMA7rGH98QZEb0+rJcddu+9W5qgkuacc4Czzkq3jKDc\neCPw2Wfhz7vvPmuPy4uhQ933eVWyn37qPgPezhZbtFQ83bpVB4+Ncnr2bFlmljPkL7hAfzuNExm4\n9VZ23bX6PiRlkm1o0J5tzzxTjVUWBCcZi9IAikoagROLRKEVwymnAN98k07evXvHm2oP6BfkyiuT\nkceNq64qzhoS665b3JaWV+UZlrPO0qY9J++tk05Krhw/DJv573/vfkwe/vYdOvhPAmtuBlauzEae\nrLn5Zm0yruW5DoVWDCRbsgiP7NRSdGs9prlanRe/+IV1jMlMkp5HQTjjjGDu1kVj002zHf+w47e2\ne5z/ccQIruDmiYhcISKzReR1EXlERDqY9o0WkXmV/ft55UPyR6nihHIwMDxR9tgjm/LyWlDGi+uv\njz+zvl+/4McOHhyvrLJQdlNW2sRt/0wEsK1Sqj+AeQBGA4CI9ANwFIC+AA4EcJMI/4qkcBo8LQtB\nnwJzNz2Ir3oST9eoUeHnzBT9qe7UyXk+iYF9ZvBGG2mz3D//GT/ktRdR7tsNNwQ/NuwcDWIllmJQ\nSk1WShnzXl8BYHTgDgZwv1JqjVLqA2il4eE5T8IQpKKsJZzGD5ziXcWloaHlLPuiVPzvvBMtyujS\npd7uxrvvrmcQG7Rrp2fLDxvmPdclKk6zw4Ny+un6220lNXMAviATE53g0rmaJC2mwwEYvifdACww\n7WuupBEHonq6FKXSSgvDRn344cVy5V24sLqdlUy9e+uZwbNm6VX29tpLp8cd8xg/XscayoohQ+Kd\nf8EFVW8tO2aTm9//4qYsa3lAOQy+E9xEZBIAc5tNACgAFyqlnqgccyGA1Uqp+6IIMdYUm7exsRGN\nIWcemSMtlo25c6Ovs7zTTtVAYmXB74Xt06fq2nrssTpI3IABuoewySbAl18GyycpnJ4tNw+oX/0q\n3dXcOneueoVtvbWeuR03TLVdsaQ1uP7JJ7rS7dpVD6hHDVXvNrs5LI2NwHPPJZNXVjQ1NaGpqSmT\nsnwVg1LKczK6iJwAYBgAc2zRZgBmK1/3SpojY/2CtvswapSePVzGae9RQwaffbY1lEaZMSr5r77S\nFd211+oKsHVrq5nj1Vf1/5yGGcmNnj2d0zfYQAdXNCuozTcPHma6aMybp2fjR22k+GGORea2XkiS\nGJNJ3XBrWBSlV+qEvdF8sRGmIQVihcQQkQMAjASwh1JqlWnXeAD/EJFroE1IvQCk1mFtaHAP8FWr\nXHNN3hJEw+nFW3ddbR4xm9Sc/s8ePaqTtvr31942s2alI6cfX3xR7EokLL16pacUsuaMM3SkVi+C\n/HfrrAOsqtRqX38dX64yETdW0g0A2gCYVHE6ekUpdapSapaIPAhgFoDVAE5VKp71buBAbWMl5cYt\nim2QdQUAvXbC/Pna62TmzNqqnEl8pkzRs77jLDdqsMsuwIsv6u16e85iKQallOuUKKXUpQASigeq\nQ1qHWQOZFI81a4KHk3Zjt93cGwhZhqswqLcKo+gYg/J+cJDZm0JHVyW1RVyl4IYxHjdkCPDuu+mU\n4cSOO4ZfBjQpzjyz3E4XRSVMyPVahoqBlJopU6rrYIiEj+QZh9dfz64sO337Wpd/JeEIW+nXm5Jg\nrCRSaoKaDggJQtzwI7UCewyEEALg7bfDL0pVq1AxEBKCLE1VJFu8FkSqN6gYCAnIN9/kG0qaJEfY\nKLL1NsZQU4ohLa8XUr/svDPw29/q7bw8kEjyuM2nIZqaGnw+/PCq6yIpN1Fj6SRN+/bA1VfnLQVJ\nGs5j8KamFEPr1lXXRVJeZs6MHkOKkKQwjznUmylJYkaqiC+ASNxoGYQQEoo33mi5roO9Glqzptpz\nXbUKaNMmG9mCIiJQSqWismqqx0AIIUnRuqZGYMNBxUAIqTvqzTQUFioGQkjdEdZ6XW+KhIqBEEKI\nBSoGQgghFqgYCCF1zxtveO+nKYkQQmoce5SEHXbIR46iwnkMhJC6QymgVSvrbyeMnsLq1cVzX+U8\nBkIISRAu1OMNFQMhhBALVAyEEEIsUDEQQogPNCURQgix0KrOaso6u1xCCCF+UDEQQgixQMVACCHE\nAhUDIYQQCwWby0cIIcXhvPNarvRWDyQSEkNEzgXwJwCdlVJLK2mjAQwHsAbAWUqpiS7nMiQGISRz\nzC6oZayC0gyJEbvHICLdAewL4ENTWl8ARwHoC6A7gMki0psagBBCik8SYwzXABhpSzsEwP1KqTVK\nqQ8AzAMwKIGyCCGEpEwsxSAiBwNYoJR6y7arG4AFpt/NlTRCCCEFx9eUJCKTAHQxJwFQAC4CcAG0\nGYkQQkiN4KsYlFKOFb+IbAdgCwBviIhAjyVMF5FB0D2EHqbDu1fSHBk7duz3242NjWhsbPSXnBBC\n6oimpiY0NTVlUlZiC/WIyPsABiilvhCRfgD+AWBXaBPSJACOg8/0SiKE5IHhlbTffsCECfnKEoVC\neyWZUNBmJiilZonIgwBmAVgN4FTW/oSQotGhQzmVQtpwaU9CSF0iohXD8uV5SxINLu1JCCEkM6gY\nCCGEWKBiIIQQYoGKgRBCiAUqBkIIIRaoGAghhFigYiCEEGKBioEQQogFKgZCCCEWqBgIIYRYoGIg\nhBBigYqBEEKIBSoGQgghFqgYCCGEWKBiIIQQYoGKgRBCiIUkV3AjhJDScOmleqEe0hKu4EYIISWE\nK7gRQgjJDCoGQgghFqgYCCGEWKBiIIQQYoGKgRBCiAUqBkIIIRaoGAghhFigYiCEEGKBioEQQogF\nKgZCCCEWYisGETlDRGaLyFsicpkpfbSIzKvs2y9uOYQQQrIhlmIQkUYAPwWwvVJqewBXVtL7AjgK\nQF8ABwK4SURSiemRN01NTXmLEAvKny9llr/MsgPllz9N4vYYTgFwmVJqDQAopZZU0g8BcL9Sao1S\n6gMA8wAMillWISn7w0X586XM8pdZdqD88qdJXMXQB8AeIvKKiDwrIjtX0rsBWGA6rrmSRgghpOD4\nrscgIpMAdDEnAVAALqqc30kpNVhEBgJ4CMCWaQhKCCEkG2KtxyAiTwG4XCn1XOX3PACDAZwMAEqp\nyyrpTwMYo5R61SEPLsZACCERSGs9hrgruD0GYG8Az4lIHwBtlFKfi8h4AP8QkauhTUi9ALzmlEFa\nF0YIISQacRXDnQDuEJG3AKwC8CsAUErNEpEHAcwCsBrAqVymjRBCykHuS3sSQggpFrnOfBaRA0Rk\njoi8IyLn5yjHX0VkkYi8aUrrJCITRWSuiEwQkY6mfY6T90RkgIi8Wbmea03pbUTk/so5L4tIj4Tl\n7y4iU0RkZmWi4ZllugYRWUdEXhWRGRX5x5RJ/kr+rURkesWMWjbZPxCRNyr3/7USyt9RRB6qyDNT\nRHYti/wi0qdy36dXvpeLyJm5y6+UyuUDrZTeBbA5gB8AeB3ANjnJ8mMA/QG8aUq7HMDvKtvnQ8/X\nAIB+AGZAm+G2qFyD0fN6FcDAyvZTAPavbJ8C4KbK9tHQczySlL8rgP6V7fYA5gLYpmTX0Lby3QDg\nFeh5L2WS/7cA/g5gfAmfn/nQ3oXmtDLJ/zcAJ1a2WwPoWCb5TdfRCsAnADbLW/7ELy7ETRgM4F+m\n36MAnJ+jPJvDqhjmAOhS2e4KYI6TnAD+BWDXyjGzTOnHALi5sv00gF0r2w0APkv5Wh4DsE8ZrwFA\nWwDTAAwsi/wAugOYBKARVcVQCtkreb4P4Ie2tFLID6ADgPcc0kshv03m/QA8XwT58zQl2SfBfYxi\nTYLbWCm1CACUUgsBbFxJd5u81w36GgzM1/P9OUqptQCWiciGaQgtIltA935egX6wSnENFVPMDAAL\nAUxSSk0tkfzXABgJPb/HoCyyoyL3JBGZKiK/Lpn8PQEsEZE7K+aY20SkbYnkN3M0gHsr27nKz+iq\nwUlylD4VF10RaQ/gYQBnKaW+QkuZC3sNSqnvlFI7Qbe+B4nItiiB/CJyEIBFSqnXffIsnOwmhiil\nBgAYBuA0EdkdJbj3FVoDGADgxso1fA3dqi6L/DpDkR8AOBh6kjCQs/x5KoZmAOZBkO6VtKKwSES6\nAICIdAWwuJLeDG0DNDDkdku3nCMiDQA6KKWWJimsiLSGVgr3KKUeL+M1AIBS6ksATQAOKIn8QwAc\nLCLzAdwHYG8RuQfAwhLIDgBQSn1a+f4M2gw5COW494BuGS9QSk2r/H4EWlGURX6DAwH8R1XjzeUq\nf56KYSqAXiKyuYi0gbaJjc9RHoFVk44HcEJl+3gAj5vSj6mM9PdEZfJepbu3XEQGiYhAz+kwn3N8\nZftIAFNSkP8OaBvjdWW7BhHpbHhdiMh6APYFMLsM8iulLlBK9VBKbQn9DE9RSv0SwBNFlx0ARKRt\npacJEWkHbed+CyW49wBQMbcsED3BFgCGAphZFvlN/By6YWGQr/xpDKKEGGw5ANqDZh6AUTnKcS+0\nN8AqAB8BOBFAJwCTK/JNBLCB6fjR0N4AswHsZ0rfGfqlmgfgOlP6OgAerKS/AmCLhOUfAmAttGfX\nDADTK/d2wzJcA4DtKzK/DuBNABdW0kshv6mMPVEdfC6F7NA2euO5ect4D8sifyX/HaEbmq8DeBTa\nK6lM8rcF8BmA9U1pucrPCW6EEEIscPCZEEKIBSoGQgghFqgYCCGEWKBiIIQQYoGKgRBCiAUqBkII\nIRaoGAghhFigYiCEEGLh/wHjTnp/SG4rlgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x105203cc0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(f0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": true
   },
   "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.18e+02       2.73e+04    \n",
      "       7             10         1.3422e+04      1.37e+01       6.84e+01       2.20e+03    \n",
      "       8             11         1.3418e+04      3.70e+00       1.28e+02       7.88e+03    \n",
      "       9             12         1.3414e+04      4.19e+00       2.64e+01       6.24e+02    \n",
      "      10             13         1.3412e+04      1.88e+00       7.55e+01       2.61e+03    \n",
      "      11             14         1.3410e+04      2.09e+00       1.77e+01       4.97e+02    \n",
      "      12             15         1.3409e+04      1.04e+00       4.02e+01       1.32e+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.32e+03.\n"
     ]
    }
   ],
   "source": [
    "t0 = time.time()\n",
    "res = least_squares(fun, x0, jac_sparsity=A, verbose=2, x_scale='jac', ftol=1e-4, method='trf',\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 33 seconds\n"
     ]
    }
   ],
   "source": [
    "print(\"Optimization took {0:.0f} seconds\".format(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 0x10de73438>]"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYYAAAEACAYAAAC3adEgAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xm81NTZB/DfAygIKKIIVGVxF6yKqIjF6ogbVkWktiLi\ngopYF6yv1qX0Uy6uILXa+tatovXFFbUqWKXicrUuVRSpKIggZREFUSxW3Ljc5/3jJE6SSWaSmWSS\nuff3/XzmM5lMlmdmMnmSc05ORFVBRERka5F2AERElC1MDERE5MLEQERELkwMRETkwsRAREQuTAxE\nRORScWIQkW1F5DkReVdE5orIGGt8RxF5WkQWiMjfRaRD5eESEVHSpNLrGESkK4CuqjpHRNoDeBPA\nsQBGAvhMVa8TkUsBdFTVyyqOmIiIElXxGYOqrlTVOdbwlwDmA9gWJjncbU12N4Ahla6LiIiSV/EZ\ng2thIj0B1AP4IYDlqtrR8d4aVd0itpUREVEiYqt8toqRHgZwgXXm4M047HuDiKgGtIpjISLSCiYp\nTFHVx63Rq0Ski6qusuohPgmYlwmDiKgMqipJLDeuM4Y7AcxT1T84xk0DcJo1fCqAx70z2VS1Zh/j\nxo1LPQbGn34czTH+Wo69KcSfpIrPGERkAICTAMwVkbdgiox+DWAigKkicjqApQB+Xum6iIgoeRUn\nBlV9GUDLgLcPrXT5RERUXbzyuUK5XC7tECrC+NNVy/HXcuxA7cefpFibq5YVgIimHQMRUa0REWjG\nK5+JiKiJYGIgIiIXJgYiInJhYiAiIhcmBiIicmFiICIiFyYGIiJyYWIgIiIXJgYiInJhYiAiIhcm\nBiIicmFiIKJma+VK4Mgj044ie5gYiKjZmjULmDEj7Siyh4mBiIhcmBiIiMiFiYGIiFyYGIiIyCWW\nxCAik0VklYi87Rg3TkQ+FJHZ1mNQHOsiIqJkxXXGcBeAI3zG/15V+1oP1v0TEdWAWBKDqr4E4HOf\ntxK5HykRESUn6TqG80RkjojcISIdEl4XERHFoFWCy74ZwBWqqiJyFYDfAzjDb8K6urrvh3O5HHK5\nXIJhERHVnvr6etTX11dlXaKq8SxIpAeA6aq6R8T3NK4YiIiimD4dGDwYqMVdkIhAVRMpro+zKEng\nqFMQka6O94YCeCfGdRERUUJiKUoSkfsA5ABsKSLLAIwDcLCI9AHQCGAJgNFxrIuIiJIVS2JQ1eE+\no++KY9lERFRdvPKZiIhcmBiIiMiFiYGIiFyYGIiIyIWJgYiIXJgYiIjIhYmBiIhcmBiIiMiFiYGI\niFyYGIiIyIWJgYiaLeGtxHwxMRARkQsTAxERuTAxEBGRCxMDERG5MDEQEZELEwMREbkwMRARAVi+\nHFi/Pu0osoGJgYianfXrgQcfBI45Jj+ue3fg+uvDzT9vHtCrVzKxZUEsiUFEJovIKhF52zGuo4g8\nLSILROTvItIhjnUREVXqxReBYcMKx3/2Wbj5X30VeO+9eGPKkrjOGO4CcIRn3GUAnlHVXQA8B+Dy\nmNZFREQJiiUxqOpLAD73jD4WwN3W8N0AhsSxLiKipKimHUE2JFnH0FlVVwGAqq4E0DnBdRERVU1T\n72OpVRXXFZiL6+rqvh/O5XLI5XJVCIeIqHbU19ejvr6+KutKMjGsEpEuqrpKRLoC+CRoQmdiICJK\nS9iipDTOGLwHzePHj09sXXEWJYn1sE0DcJo1fCqAx2NcFxERJSSu5qr3AXgFwM4iskxERgKYAOAw\nEVkA4BDrNRFRzWMdQwiqOjzgrUPjWD4RUTWwVZLBK5+JiMiFiYGIyPLRR+a5oaH4dE29KImJgYjI\n8uCDwPTpwEYbBU9z8cXAkiVVCykV1byOgYgoE7xH/Lfdlh9evLj4vNdfD2y1VfwxZQnPGIio2bvo\nomjTr16dTBxZwcRAREQuTAxEROTCxEBERC5MDERE5MLEQERELkwMRNTsrVuXdgTZwsRAREQuTAxE\nROTCxEBERC5MDBSrDz5IOwKi0op1gseut5kYKEZffgnsuGPaURBRpZgYKDYbNqQdAVE4a9ZEn+eF\nF4B9940/lixiYiCiZmf9+ujzzJgBvPFG/LFkERMDESXmv/9NOwJ/U6emHUG2JX4/BhFZAmAtgEYA\n61W1X9LrJKJs2Gwz4Lvvit/4Jg1ZTVhZUY0b9TQCyKnq51VYF6Woqd/usByLFwM77NC8W7o0NqYd\nQSFuq8VVoyhJqrQeosxZujTtCIiiq8YOWwHMFJFZIjKqCusjIqIKVKMoaYCqfiwiW8EkiPmq+pJz\ngjFj6rDFFmY4l8shl8tVISx/330HtG0LNDSkFoKv4cOBL74Anngi7UiIKA319fWor6+vyroSTwyq\n+rH1vFpEHgXQD4ArMRxySB2OPTb8Mj/6CHjkEeD884OnWbbMXHDVu3e0eNety2Z7/EcfBb75Ju0o\niCgt3oPm8ePHJ7auRIuSRKStiLS3htsBOBzAO97polbMTZ4MjBkDvPRS8DQDBwK77RZtuUB2K6Wy\nGhfVtmOPBY45pvLlXH89cNJJ/u9lcdttwVrPopI+Y+gC4FERUWtd96rq096JoiYGe0P78Y+D5y3n\nApYsy+KfK0veecccCPB7imbatHiWc8cdwHvvAffeG8/yksbtpLhE86aq/ltV+6jqXqq6u6pO8Jtu\n6FDgyitLL++VV5Jv+pbVDcYZ1+WXAw89lF4sQdL67pYuBXbfHXjuuXTWXwvWrmVRZFjNuWmxLTMn\nVDNnlp5mwADz5y+2A9qwAXjqqfLjSDsxPP+8fx2HM64JE4DrrqteTHF4663klt2zp3n++uvk1pFV\nN9wA9O9ferpOnYATT0w+nrjcdBPw7bdpR9F8ZSYxhM3Shx1mKmKDvPAC8JOflB9H1MSwYUN+hzRj\nRuWJZeBA4OmCwrboy129urI4bEuWVH4EtWIF0LdvefPuvTcwaZJ73Cef+PdZE+VscskSYNGi8mKq\nliuvBA4/vPg006YBr71WelkNDcDChfHEVQ1jxpjfeOxYoK4O+Oqr8pazaBGw006F49M+AMy6zCSG\nKH/qN98Mfi/MTuyWW4Bddw2/vmIuuMA0bwVMGWuQ1q3DX+zk911E2ZA//RTo3Dn89MVstx3w8sul\np2vfHvjwQ//3Kmn6O3s28OST7nGjR/v3chklge2zj/8OI0umTs2fSd98M/C//1vZ8uIqhr3ySuDU\nU+NZVjGqwDXXAFdfDbRrB5x3nmlOHsWbb8ZzADBjhjlTby4ykxjiKNdTBQ491AwvW1b4vgiwcqU5\nIl+woPL1AcC8ee71B/nuO+DMM6Nv2F7eYqYHHig8mqp0HV5hbpS+bl08R6RffgmMHFk4funS/OcM\nak4cppmxiGmYUGvl7eeeW7x5dhjl/sfWrjVFnGPGmNe33w783/+Fn7/co3NvvH/6U3yV5VE9+2w6\n601Lk0oMYaxYUfz9MBux81qHKBv9M88A//M/4af3i+v2293jTzwRePxxc1RuT5PV0+S//KX0znv+\nfDOdV8+ewIUXmmH7882aBUycmJ8mzJkNYGLI6nfk552CBt7lKfc/tvnmpojzppsqW05UfuuJa91R\nf//m1rw1Mx/XPs2dO9dUqEW1ZEm+ErISYTaY9u1N2WdDg7slTJiNttjRd7FyVDuuV18tfO/66907\n3Lh3evbyhgyprBnwyJFAqzIaSNsXe9o3V7Hjue464LLL8tOFLSppbAz+oy9cCOy5Z/QYk3TBBfEs\nJ66daqmDq7jW+/775S3r669NEWScaulAIg6ZSQyvv26er7suf1Q9eHDpDeG3vwVOP92UhfsVH0UV\ndgNYuLDwD3LRRaXnK/Z57KKYo48OjsvvlPbNN93FR1GPbl591Xznp51WfLrHHzfdcgSxP1slO6Cz\nzio9TdDnC5sYVIN/59deA95+O9xywq6rmssQyf+XkoilGs47zzyP8ulZLcxn+N3vTKOFIFtvXVhv\nVUpzSwzV6CupLA0NwPTp5ki42FFmmOsfvMLcCNy581ixwhzNF6us/PjjcOsu989pxxI0v/NMJMpG\nrAr86Ef518OHm8q6c87xX16xsvmg2BYv9h//zjtAhw5At275cXPmBC+/VFFZ2O+2sTF4GZWUh993\nn6m/uvvu8pYRFMvzz+fHffmlOWMtZv58oJ/PXU+c38+6daZCtxoaGoCNNw4//Z/+FPxemORfqv4o\n7H/ViUVJGWFvHI2NZsdSbsVmuc02nX+igw8Gdt65+PR//nP05UbhTQzFbjFYauemai4qPPHEwngu\nushUdAYJ833aSdVupTRwoP90u+8ODBqUf71qVellA8klhunTgeXLwy1j0SL3Mj7+GBgxorBSNqhF\nzIoV7ms7VE0z3FKOOir4PecOcfXqws9o71Q3bCidXML4xz+C3xs1CvjZz8xwu3amf7M4pHXW09zO\nGDKVGKZPzw//8pfm+bHHTBv4UjvmILmc2ZiiHiW8/z5w5JFm2E5Kw4fnyy7/+ld3nYa3UrWx0f/I\nJepFWN98Y7p6+Owz87rUH+Pll0vvYEeMMNeCPPRQ4fKcFxXZ8Yu4j1qPO85/B+pd1mefme/J67vv\n8p/HWQQWtjVV1MTQ0GCafto7xsZG/3UNHmyuKvfGCpg7kU2das4cL7648EAlaN1BBzQnnOC+tuPB\nB4EuXcywCHD//f7z2cWlDQ3Ar37lfu/TT/PDK1cWzrt4sWlSHaWTyKA4gOJFj/ffDzz8cP71Bx+E\nX2cxlSaGMHVkcVxHVOsylRgGDy4c9/rrlV0BuWaNaYO89dbmddgziJkzzXxO998fXHbpPPoTAa69\nFthkk8LpHn4YuPPO8OXhRxwRvkksABxwAHD22cWnefFF89zYWLg8547MrkNRdZ+hPPZYfhlOfnUM\nP/1p4XStW5srcQH3EXWp03UR4J573Dscv/XbScf28stmR2wXt0Vpz9+6tTm7+u9/TV3M3Lmmsr+Y\n3r1NUn3lleDfyxvDHXe4X8+ZU7yhwurVpiw9ql69ou1cvWempSp17TP9pIpeovx2fkVS3sTvZ9as\nwnEsSkqZd6O9/vrK25w7Lzz79NPg7H/++cAZZ5S3Dm85unNn7nXGGUDLlu5xIu5uQeyeY71NFcP8\nqf2OFm0NDfn3Sy3LvtbjzDOBW28tvV5b2OIYmwhw0EHh/nwnnxz8nv15OnVyt/DyLrdTp3DXZtiV\n0PaV9sUqrZ3mzzcHJAMGBDehtWOyzxL8zvKCtvt//9s/Dvt7f+SR4kfGpQ60pkwJ/pzOAyO/mO2K\n46Df0u/sff36cMVo9vylXHONefa7EDbMRaZ+6+AZQ8peeCHe5a1cWXjtQFB33bfdZk7pgeinrO++\n6359332l55k7112kccst+eEf/9g82000bc4znqAYi13pOWmS+0rkYp/T3jEuWxZcgRyXF190f/5y\nNDa6y9FtYRKOX4s2+4IuJ3sHUaoF2m23FcbhZMdk7xD9DgD8yuWXLAG2395/p2c3Ipg+3RR7Au7i\nJZtffdi6dfkkdsop7jiCrFtnKsP9Eubatf7z2Gfvzu3+kkvyCdKPM4bx48NfXOcXu/ds0o/3Pwcw\nMaQuqFuFONl/RhHzePVVcxTl3JA+/zyeda1fH7yjHjgw350GEH3n+8or+WG/jXnKFPfrL75w1xUA\nwX/8jz6KfvrsLUqKmlzLaWHmXb+9Q/r44/y2FOZP3aNH4TjvdSWLFuWXP3++efZ+xza7bN6vWAII\n/m7t4stSF7XZsQVdwW+P97tS2K+eq317UwzZtWvx9Xp17154lv3HPxZO9957puze7sds8uT8e85h\nP87kunhx8HUdW20FbLNN/rVqYV1SsebWxQRtQ3H3MpAZqprqA4Can7A6j3vu8R//hz/4j//rXytb\n30knlZ5G1X/8qlXh1zNhQvCybX7vf/ON//iLLlKdNKlw/EYbub/HKVPy711xhXn+4APz/Oyz5X9v\nxWK2H8cf73595pmqa9aY4Q4dzPOFF6r+7W9m+Isvwq+noSF8rIsWqX74YbjPYzv4YP/f/+KLzfP2\n20f7vt5/P/i9o45yv7Z/p7Vriy/zkEPMd1pq3TvsoNqrV7R4r7mm8Lvfc89w/xf7O7v/fjP80EOq\n06cXTrPffoXf/957R98Gv/hCdbPN/Kf55pvC37ZazO47of1yUgsOHQCqmxiCHhMn+o8fMCD5dTt3\nruU+Bg3yH6/WrxwlydiP1q2D37v11uD3pk41z5UkhilTVL/6Kto8Z5yh2qlT4fj+/c3zU0/5fz/2\nd1Tuo64uWmJobAxODM3lceWVlX33qqqHHhptelXVffaJPs/kycHTMDEkFQCykRia6kOtXzkocST5\n+OMfy593yJDo83TtWvz9004rHNeunerMmZV9zrPOUn3ttdK/g23sWPd7v/lN+ttJtR+HH+7cwUV/\nqJb+vf2+/6jrePDB4tOsXRtiD54QJgY+yn6MGJF+DFl5tGrlP95bJJXEw/2Hdj+iFhs1hUfPnsHf\nR9jvM8r0b76pet550dcxbFjp6d54I+IePSZJJgYxy0+PuR90ujEQJU3VXJC2++6F9zLo2rV4E+Om\nqE0b0yrvs8/C3YHO64kn/PsUS8OeexbvyiUpIgJVTaS9VOKJQUQGAbgRpgXUZFWd6HmfiYGavEce\n8b/Yj5qGNI6vazYxiEgLAO8DOATARwBmARimqu85pmFiIKKa1tQSQ9LXMfQDsFBVl6rqegAPADg2\n4XUSEVEFkk4M2wBwdpDwoTWOiIgyKiP3Y6hzDOesBxER2err61Fv384wYUnXMfQHUKeqg6zXl8E0\nsZromIZ1DERU01jHEM0sADuKSA8R2RjAMAA+vbcQEVFWJFqUpKobROQ8AE8j31x1fpLrJCKiyiTe\nu6qqzlDVXVR1J1WdkPT6iCj7fvc7YNNNgWOr1EbxBz9Ibtl2F/lNSea63SZqbrbcMu0Iqq9PH9MF\n9mOPVWd9Ye6P4jViBLDxxqWn897vpSlgYqDEtGmTdgTlK3av43KsWWPukXzSSYXvtcpI28Bq2nzz\nyua3b6gVVrdu4W/wY5sypfStQBcuBIYMibbcWsDEUMIJJ+TvT1zL7DuKRTFuXGXrtG+JmRVRihOG\nDSscN2ZM+UUSHTsCd91l7lntvbHNYYeVt8w47btvddfnTIb9+kWf/+c/Bzp3zr8eOrT0PCeeGH09\nQXfgs0W9sVGtaPKJYezYyuZ/4IHkb+s3bJj/XbXiFOaP453+gANKTzduHPDWW+5xffuaZ+99raO6\n7DLz/Le/lb+MqVNLT2PfpzjIQQeZ5w4d8ve+7tXLPc2ee4aPyXtD+ylTgH32CT+/7cMP/c9Awnr1\n1fzw66/7N7lMqn8n53+q3P+X3aS/ZUugrq709OWcmXl/K6+mesvPJp8YrrqqcNwvflE47q674ltn\n9+7u1w88UHz6li2TKXbZdFPzPHasOeuJso5DDsnvKLbaKni6ujpTXuxk34+4kj/NwQcD115rYthl\nl/DzzZzpfv2zn+WHg+IZONB//HbbmecBAwrf896O0t5+/G6l6TV8uInLeU/nQw4pPZ/Trbea21gG\ntZ+3k7Mf+2Bpn32AjTYKnu6bb4B27QrH+43bfvvg5Tj5nbk6b6u6//6ll7HDDubZTs5//rPptbZH\nD/d9w/1uMRr17ORXvzLb1ODB0eardU0+MQBA69b5Yb8/0y23mPLfrbd2j3fejzmKJ54wG5Qt6M/X\nvr15jqOM2W+nb9+f2N4hdu1qijSKGTTIPKvmj5ac3UTbp87duwOnn54f7/3MxTjvVW275BL3a+dO\n3N4RlDJxYr4LZ7u3fKD4zk8VOO44//euuMJ/esB9JDlyJLDFFmb4mGPMs3dbcjr0UHMm84Mf5Jdn\n7+TOOit4Puf67WIU+/Ubb+SneeEF4O9/z7++9NL88PLl+QMXv0S5fj0wwWo76PzfONmJ4aWX8gnN\nWazT0OBfadu7d36dznXb94EOisnLPuABzOc++WQzvGQJcPbZ+ftpn39+fjr7M3foUHr5Tptvbn6v\nnj393/dLkk1Bs0gMTt4biT/+uNmYgOBEEPXIt0UL958qaMdk/3mC3v/qq8L6jW7d/Kf1O9Ly/gln\nzTI3ZXey/1Q2u9JVNb9Dvu66/A7I7nd+6FD3UfNRRxUWawR9rt69C8dNmFD6tL2UvfbyP4KePt08\nn3JK6WU4jwxHjDBHuOecY147twNnrJttBmyyiXs55TZhdB5R77pr8HT2b2PHsffe+ff239+93TiL\ng9q0ycfqt123amUOnmzFihMHDDD/pwsucH8ffkWI06YBb7/tvxzn9hC0A7YdeKC7FGDvvQsPrH75\ny/zwDTeYZ3tbbKo78rhlPjFUUr5sc+4sRPKvlyzJH+F5p/N7HZa3+WGxI1bA7NAA91EXYP7A3mTl\n/MMvW5YfLlZJZu8AOnUqXIe3pYbdWkQV2HFH8+y3A/GWsQNmB+kU1NTPeyR6zDFmHWES8KRJxd8v\n9l1vu23xeYcOLUxaZ51ldpTdupmdvb1NOCsdW7Qwv7mzPuuoo4qvqxS/xLJqlXlWBfbYIz9sO/RQ\n82zvmJ99tnAZLVqYoiwg+Pt27uRHjSoe5zHHADfemJ/noYf8p2vVysRlb7/O4tZ99zWfY80ac7Dh\n3UadLr649Hd75JH54RNPdB+w2Gd2YeokmrNMJYbnny8cZ/8Bwnr44cJxQUeiPXr4/zl++EPz7E0M\nv/lNuBi6ds3Pe9JJ7nLNNm1MpWFDA/DrX5siFLvOw+80Nyj23r3dO1g7MSxYUFhU4/2M9g7miCOC\nP4PzdN1PC58tZ7vtzB25bEEtNtq0Md/PVVeZCkS/cvmLLnK/7t/f7HQuvrhw2t12cy/722/d7++z\nj2nFUsqWWwYfDCxblt/xAsBOO+WH+/Y1Oz3nkezuu5den9MJJ+QPgkaMAM49t3Aavx2mM95p08y1\nAfZv41d30qKFifW774ITQzkHRPY8xx/v/74d05AhZvvz29Y7djS/n7PupRzOOosuXUxLMK9SB2vN\nXaYSg98pqPcoz9t++d5786eSnTub02bnhgG4N/Sddw5evz3d3Lnu12efbcqRnTsXbxND59G70z33\nmKMUu1KsRQtzBNqypdn5TZzo/we1i3S8RUfOsmu7XuHoo/M7om23LWwj7t2hdO5sjticZ0tOixYV\nFjF5Be1U7CMy1XzlbZCxY/Mtfryc5c6AaUHj3enYV82+8457vPdMZcstzXaTyxW/jWTUtvXvv28S\n94gRhe+JRGup1KZN/jNPmWISRRjOA4dNNgmf0IvtGPv39z8jDBsHYOpPnPVZ9npFip8RxGGnndLp\n1K4pyVRicJ4d7Lef/6mwt5hm+PB8qxj7SPhHP3JP49xI2rYN3mgmTDDl6YDZ2drlzePHA3feaXa+\n9rzeSlzvDty7DjtZFSsucc5jt6OfMcM9jXNHahfd9OgB3HSTGW7Rwr2D++gjYPTownUtXep/VAqY\n8utSTU39zhjiEKbCv0uX8pa9227uJppOL74Yrmmz8zfaaafg31PE1MeUapFWqag7wKAKZadddgHm\nzYu2XG9i8NZnVdJCLc4WexMn+pdMOJVzlXRTk6lrLjt0AK65xhSxbLqp+1S4a1dT7OLd6QPmSNLZ\nKsP+s9gbt7eOIYizaaNdYennX/8yR+bF2sifcopp4WGz4y5WmbflluZo3clbbu81b545+rd35CLu\n1i5hL8iaMiVak90022/vtZdJmOXG8PDDhWcfcfZ306ZN/ky30it8SwlTYW9vG336FCaGuH7Hc891\n/wedTZwnTTKVxlnQqZM5c7QPIv0kddBTSzL3FRS7BP3CC92tP8aPN8/nn+8uV/c2AYz7tHKPPUyR\nybJlwZ2A7bJL/uwDMElv1ari7dyDjtLXrAmep1cvd0uLcv/ogwf7n6EFCbseZ5NBwL9tudM11/hf\nexK0bG9dRBg//Wm4PnD82EVlxXz9dbjp4lBq237llXzjhjiMHOk/ftSo4KvrTz452vddjWKgYsV8\nfusv56rpWpaZxODcid5+u9lB2G6+2X+H8tvfmue2bcNdGAPEe6TbrRtw9dXhp+/cufgfJKiytmPH\n8Ecx1TqSD3vthfN369GjMFF4XXBB+KvVBw40vXTGqdRO6aCDgJUrwy+v0qu/Syl1xrD//vFuE3fe\nGX2eqOtv2dLd/5Cz2W5cHQ5GTT79+5d3EFKrMpMYnJW5o0a5+275xS/cxTxZNXt2vqy/HN4K13KU\nuxOIMt8LL4Rr5WOzmxd6K4nj5r2WoJhKLiqMUsdx8MGly7SD2AcRxc6gaqGSNeo2KeLuZ+ull8zz\nDTf4FyXHrRa+06RlJjHEWa7n98M6L/aK84ffYYd85e5ee5Xud6cY+0roYkqdGVWjfPTAA6M197Pr\na5K8uGjBAv+uK4L861+FF/sBpq6h3MptPy1bmjLtcjz5pEmmYRssVEu/fsVb93lVesZibzelrtqn\n+GQmMcR5yu39szz6aHL9vrdpk+9YrVLFemocNMi0oy/VsV0WO/WqRkw77xxtPTvu6N8H09FHRysq\nStI225iWVMV2/pVeLR7GW2+5Lxp79NFoZ3+V/P7ldC4YRrHvlGcMGWqVFKUYICpneWXcZwxxCmrT\nD4S/AjyLiYEqk/ZOrE8f9/8zajFcudvkZpsVdtBI1ZGJxPCTn4TvnREwV/mm/WfJGpHSfcenKUw/\nReQv7TOGSpWbGFavdpckRO0Aj8qXWGIQkXEARgH4xBr1a1Wd4TdtqStkvUqVoxf7I4nU9p3FALOT\n9etWOcvtr+++O+0IateBBxZezW8bPbp0H1CAqcgu1ilfkspNDM4WfEuXBncgWQ4eWBaX9BnD71X1\n96UmcvYMGYdSP+zVVxfeRauWxL2TPfro8rsYp+QddJDp8NHP0KHhbsK0Zk1l/QNdcYXZTsoRR/Gm\n9x4naWhOxbRJJ4ZQX2XQRTPlKnXGkOWyyz32CO6eOCnFrvKmpqHSFmG77eburDCKWtuh8owh+VZJ\n54nIHBG5Q0RYQhhCLVyvkbaJE1k0VUtqLTFQhWcMIjITgLPVtwBQAGMB3AzgClVVEbkKwO8B+Bbg\n1Dk6R8+VZADzAAAKf0lEQVTlcsiV2/DbUuqMgWrbHntE746d0jFjRukeX7OkWCebaauvr0e9faPr\nhFWUGFT1sNJTAQD+DCCwwKIu5rtm7L9/cL8/xa4VyIJq9bFDVA3F7vlB0XgPmsfbncUlILGiJBFx\n7oKHAki4Q4S8Sy4xNyzxUnXfYCWLRo8Ormgkitvzz5suToickqx8vk5E+gBoBLAEgM9dAcirZcvg\npolEcatG30NZxOaqxSWWGFSVlzQ1QWzWStT0ZeLKZ6oNPJKi5oDbeYY60SMiqhbu/ItjYiAicgiT\nNMLeMrdWMTEQUbPXuXPaEWQLEwMRNXt2P1JZ7pa/mpgYiKjZY48IbkwMRETkwsRARM0OL3ArjomB\niIhcmBiIqNk7zNEdKM8YmBiIiHDnndGmb+rJg4mBiJod1jEUx8RAREQuTAxEROTCxEBEzU6xvo5Y\nlMTEQETNUOfO0RNA9+7JxJJFTAxERA5BCePcc4H//Ke6saSFiYGImrWwPau2aAF06JBsLFnBxEBE\nzVrv3u7XrGOoMDGIyPEi8o6IbBCRvp73LheRhSIyX0QOryxMIiKqlkrv+TwXwHEAbnOOFJFeAH4O\noBeAbQE8IyI7qTIXExFlXUVnDKq6QFUXAvD2Zn4sgAdUtUFVlwBYCKBfJesiIqoGHr4mV8ewDYDl\njtcrrHFERJnSrVvaEWRPyaIkEZkJoItzFAAFMFZVp8cRRF1d3ffDuVwOuVwujsUSERX16adAu3bu\ncVk9Y6ivr0d9fX1V1iVxFPuLyPMALlLV2dbrywCoqk60Xs8AME5VX/OZl1UPRJQ6EZMkJk0Czjkn\nOEGIAF27AitXpptERASqmshNSeMsSnIGOA3AMBHZWES2A7AjgNdjXBcRUSJ4nFp5c9UhIrIcQH8A\nT4jIUwCgqvMATAUwD8CTAM7haQERUW2oqLmqqj4G4LGA964FcG0lyyciqjYewvLKZyIi8mBiICKy\nbLRR2hFkQ6VXPhMRNRlnngkcdRQwe3bakaSLZwxERJaNNwZ69gQmT047knQxMRARkQsTAxFRRE29\n5RITAxERuTAxEBGRCxMDERG5MDEQEZELEwMREbkwMRARkQsTAxERuTAxEBGRC/tKIiICcOONwJAh\naUeRDbHc2rOiAHhrTyKqISJAly7AqlW8tScRETUTTAxERORS6T2fjxeRd0Rkg4j0dYzvISJfichs\n63Fz5aESEVE1VFr5PBfAcQBu83lvkar29RlPREQZVlFiUNUFACAifhUgiVSKEBFRspKsY+hpFSM9\nLyIHJLgeIiKKUckzBhGZCaCLcxQABTBWVacHzPYRgO6q+rlV9/CYiPRW1S8rjpiIiBJVMjGo6mFR\nF6qq6wF8bg3PFpEPAOwMwPcW23V1dd8P53I55HK5qKskImrS6uvrUV9fX5V1xXKBm4g8D+BiVX3T\net0JwBpVbRSR7QG8AGB3Vf2Pz7y8wI2IagYvcCtBRIaIyHIA/QE8ISJPWW8dCOBtEZkNYCqA0X5J\ngYiIsoddYhARRcAzBiIianaYGIiIyIWJgYiIXJgYiIjIhYmBiIhcmBiIiCKYMgUYPTrtKJLFxEBE\nFMGIEUCnTmlHkSwmBiIicmFiICIiFyYGIiJyYWIgIoqoffu0I0gW+0oiIopowwZg0SJgl13SiyHJ\nvpKYGIiIahA70SMioqphYiAiIhcmBiIicmFiICIiFyYGIiJyYWIgIiKXihKDiFwnIvNFZI6IPCIi\nmzneu1xEFlrvH155qEREVA2VnjE8DWA3Ve0DYCGAywFARHoD+DmAXgCOBHCziCTS3jZt9fX1aYdQ\nEcafrlqOv5ZjB2o//iRVlBhU9RlVbbRe/hPAttbwYAAPqGqDqi6BSRr9KllXVtX6xsX401XL8ddy\n7EDtx5+kOOsYTgfwpDW8DYDljvdWWOOIiCjjWpWaQERmAujiHAVAAYxV1enWNGMBrFfV+xOJkoiI\nqqbivpJE5DQAowAMVNVvrXGXAVBVnWi9ngFgnKq+5jM/O0oiIipDJjvRE5FBAK4HcKCqfuYY3xvA\nvQD2gylCmglgJ/aWR0SUfSWLkkq4CcDGAGZajY7+qarnqOo8EZkKYB6A9QDOYVIgIqoNqXe7TURE\n2ZLqlc8iMkhE3hOR90Xk0hTjmCwiq0Tkbce4jiLytIgsEJG/i0gHx3u+F++JSF8Redv6PDc6xm8s\nIg9Y87wqIt1jjn9bEXlORN4VkbkiMqaWPoOItBaR10TkLSv+cbUUv7X8FiIyW0Sm1WDsS0TkX9b3\n/3oNxt9BRB6y4nlXRParlfhFZGfre59tPa8VkTGpx6+qqTxgktIiAD0AbARgDoBdU4rlAAB9ALzt\nGDcRwCXW8KUAJljDvQG8BVMM19P6DPaZ12sA9rWGnwRwhDX8CwA3W8MnwFzjEWf8XQH0sYbbA1gA\nYNca+wxtreeWMNfE9Kux+C8EcA+AaTW4/SwG0NEzrpbi/wuAkdZwKwAdail+x+doAeAjAN3Sjj/2\nDxfhS+gP4CnH68sAXJpiPD3gTgzvAehiDXcF8J5fnACegqlk7wpgnmP8MAC3WMMzAOxnDbcEsDrh\nz/IYgENr8TMAaAvgDQD71kr8MBd2zgSQQz4x1ETs1jL/DWBLz7iaiB/AZgA+8BlfE/F7Yj4cwD+y\nEH+aRUnei+A+RLYuguusqqsAQFVXAuhsjQ+6eG8bmM9gc36e7+dR1Q0A/iMiWyQRtIj0hDn7+SfM\nhlUTn8EqinkLwEoAM1V1Vg3FfwOAX8Fc32OrldhhxT1TRGaJyJk1Fv92AD4Vkbus4pjbRaRtDcXv\ndAKA+6zhVONn76rhxVlLn0zbY5H2AB4GcIGqfonCmDP7GVS1UVX3gjn67iciu6EG4heRowCsUtU5\nJZaZudgdBqhqXwA/AXCuiPwYNfDdW1oB6AvgT9ZnWAdzVF0r8ZsFimwE05XQQ9aoVONPMzGsAOCs\nBNnWGpcVq0SkCwCISFcAn1jjV8CUAdrsuIPGu+YRkZYANlPVNXEGKyKtYJLCFFV9vBY/AwCo6hcA\n6gEMqpH4BwAYLCKLAdwPYKCITAGwsgZiBwCo6sfW82qYYsh+qI3vHjBHxstV9Q3r9SMwiaJW4rcd\nCeBNVf3Uep1q/GkmhlkAdhSRHiKyMUyZ2LQU4xG4M+k0AKdZw6cCeNwxfphV078dgB0BvG6d7q0V\nkX4iIgBO8cxzqjX8MwDPJRD/nTBljH+otc8gIp3sVhcisgmAwwDMr4X4VfXXqtpdVbeH2YafU9WT\nAUzPeuwAICJtrTNNiEg7mHLuuaiB7x4ArOKW5SKyszXqEADv1kr8DifCHFjY0o0/iUqUCJUtg2Ba\n0CwEcFmKcdwH0xrgWwDLAIwE0BHAM1Z8TwPY3DH95TCtAeYDONwxfm+YP9VCAH9wjG8NYKo1/p8A\nesYc/wAAG2Badr0FYLb13W5RC58BwO5WzHMAvA3TDxdqJX7HOg5CvvK5JmKHKaO3t5u59v+wVuK3\nlr8nzIHmHAB/hWmVVEvxtwWwGsCmjnGpxs8L3IiIyIWVz0RE5MLEQERELkwMRETkwsRAREQuTAxE\nROTCxEBERC5MDERE5MLEQERELv8PE6d3QGfYlB4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1068b1198>"
      ]
     },
     "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
}