Pre.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# data analysis and wrangling\n",
    "# import pandas as pd\n",
    "import numpy as np\n",
    "import random as rnd\n",
    "import math\n",
    "\n",
    "# visualization\n",
    "# import seaborn as sns\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "from sklearn.linear_model import LinearRegression\n",
    "from sklearn.preprocessing import PolynomialFeatures"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "def get_test_track(s):\n",
    "    i,moves,des = s.split(' ')\n",
    "    moves=moves.split(';')[:-1]\n",
    "    moves=[move.split(',') for move in moves]\n",
    "    moves=[[int(ll) for ll in l] for l in moves]\n",
    "    des=[float(p) for p in des.split(',')]\n",
    "#     flag=int(flag)\n",
    "#     print(i,moves[0],'->',des,flag)\n",
    "    m = np.matrix(moves)\n",
    "    m = m[np.lexsort(m.T)][0]\n",
    "#     print(m)\n",
    "    return (i,(m[0,0],m[0,1]), m, des)\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "def sort_track(track):\n",
    "    #此函数用于将矩阵沿-1属性排序\n",
    "    return track[np.lexsort(track.T)][0]\n",
    "    \n",
    "def get_track(s):\n",
    "    i,moves,des,flag = s.split(' ')\n",
    "    moves=moves.split(';')[:-1]\n",
    "    moves=[move.split(',') for move in moves]\n",
    "    moves=[[int(ll) for ll in l] for l in moves]\n",
    "    des=[float(p) for p in des.split(',')]\n",
    "    flag=int(flag)\n",
    "#     print(i,moves[0],'->',des,flag)\n",
    "    m = sort_track(np.matrix(moves))\n",
    "#     print(m)\n",
    "    return (i,(m[0,0],m[0,1]), m, des, flag)\n",
    "def to_percent(mat):\n",
    "    return np.nan_to_num(mat/mat.sum(axis=1).reshape(-1,1))\n",
    "def plot(data, rows=3, cols=3):\n",
    "    for i,s in enumerate(data):\n",
    "        if i >= rows*cols:\n",
    "            break\n",
    "        start, track, des, flag = get_track(s)\n",
    "        canvas = plts[i//cols][i%cols]\n",
    "        canvas.plot(*start, '.b')\n",
    "        canvas.plot(track[:,0],track[:,1],'-g')\n",
    "        canvas.plot(*des,'.r')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "\n",
    "#按首列分布聚合\n",
    "def get_distr(mat, part=4, min=None, max=None, percent=True):\n",
    "    assert part!=0\n",
    "    assert mat.shape[1] >= 2\n",
    "    min = mat[:,0].min() if min is None else min\n",
    "    max = mat[:,0].max() if max is None else max\n",
    "    parts = []\n",
    "    space = (max-min)/part\n",
    "    for i in range(part):\n",
    "        RangeBool = np.all([mat[:,0]>=space*i+min,mat[:,0]<space*(i+1)+min],axis=0)\n",
    "        parts.append(RangeBool.flat)\n",
    "    distributions = np.array([mat[np.ix_(p, [0,1])].sum(axis=0).flat[1] for p in parts])\n",
    "    distributions = to_percent(distributions.reshape(1, -1)) if percent else distributions\n",
    "    return distributions\n",
    "    \n",
    "def theta_distr(theta, part=4):\n",
    "    return get_distr(theta, part, min=0, max=math.pi)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "human = []\n",
    "machine = []\n",
    "test = []\n",
    "with open('dsjtzs_txfz_training.txt') as f:\n",
    "    people = list(filter(lambda l:l.strip().endswith('1'),f.readlines()))\n",
    "    f.seek(0)\n",
    "    machine = list(filter(lambda l:l.strip().endswith('0'),f.readlines()))\n",
    "with open('dsjtzs_txfz_test1.txt') as f:\n",
    "    test = f.readlines()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "matrix([[ 1.,  1.]])"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def get_acceleration(track,bias = 1):\n",
    "    '''\n",
    "    此函数用于获得加速度，至少需要三个点，具体公式是解三个点用加速度表示时的方程组\n",
    "    这里将时间差强行+1，是为了防止同一时间出现了两个点造成的inf\n",
    "    固有偏差影响，视平均时间差不同，加速度总体上大约缩小了0.25%\n",
    "    '''\n",
    "    assert track.shape[0] >= 3\n",
    "    B_sub_A = track[1:-1] - track[:-2]\n",
    "    C_sub_A = track[2:] - track[:-2]\n",
    "    C_sub_B = track[2:] - track[1:-1]\n",
    "    a = (- B_sub_A[:,:-1] / (B_sub_A[:,-1].reshape((-1,1))+bias) + \\\n",
    "        C_sub_A[:,:-1] / (C_sub_A[:,-1].reshape((-1,1)) + bias)) / \\\n",
    "        (C_sub_B[:,-1].reshape((-1,1))+bias) * 2\n",
    "    return a\n",
    "\n",
    "index, start, track, des, flag = get_track(people[3])\n",
    "# get_acceleration(track)\n",
    "# print(track)\n",
    "get_acceleration(np.matrix([[0,0,0],[1,1,1],[3,3,2]]), bias=0) #测试数据，函数正确的话加速度应为1\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "matrix([[  5.02491694e-01,   1.20930233e-01],\n",
       "        [  9.01034363e-01,   1.12808053e-01],\n",
       "        [  1.29769231e+00,   0.00000000e+00],\n",
       "        [  1.58739939e+00,   0.00000000e+00],\n",
       "        [  1.79058893e+00,   0.00000000e+00],\n",
       "        [  1.50737855e+00,   0.00000000e+00],\n",
       "        [  1.01941748e+00,   0.00000000e+00],\n",
       "        [  9.73406593e-01,   0.00000000e+00],\n",
       "        [  7.09560440e-01,   1.22500000e-01],\n",
       "        [  6.11650485e-01,   1.26213592e-01],\n",
       "        [  8.15533981e-01,   0.00000000e+00],\n",
       "        [  9.51456311e-01,   0.00000000e+00],\n",
       "        [  1.09928673e+00,   0.00000000e+00],\n",
       "        [  1.36913510e+00,   0.00000000e+00],\n",
       "        [  2.44493882e+00,   0.00000000e+00],\n",
       "        [  2.21776463e+00,   0.00000000e+00],\n",
       "        [  7.19203040e-01,   0.00000000e+00],\n",
       "        [  4.17971905e-01,   0.00000000e+00],\n",
       "        [  3.32707911e-01,   0.00000000e+00],\n",
       "        [  2.78320099e-01,   0.00000000e+00],\n",
       "        [  4.46171618e-01,   1.71421721e-01],\n",
       "        [  1.15294118e+00,   3.05882353e-01],\n",
       "        [  1.17533685e+00,   2.81605641e-01],\n",
       "        [  3.61270217e-01,   2.22462203e-01],\n",
       "        [  6.17716870e-02,   0.00000000e+00],\n",
       "        [  6.06020282e-02,   0.00000000e+00],\n",
       "        [  5.08017797e-01,   0.00000000e+00],\n",
       "        [  4.93424608e-01,   0.00000000e+00],\n",
       "        [  1.52565881e-01,  -1.12444837e-01],\n",
       "        [  7.06473870e-03,  -1.87647467e-01],\n",
       "        [  3.80395867e-02,   0.00000000e+00],\n",
       "        [  2.16608392e-01,   0.00000000e+00],\n",
       "        [  3.07788102e-01,   0.00000000e+00],\n",
       "        [  3.12523191e-01,   0.00000000e+00],\n",
       "        [  2.80494505e-01,   0.00000000e+00],\n",
       "        [  2.71844660e-01,   0.00000000e+00],\n",
       "        [  2.06208791e-01,   0.00000000e+00],\n",
       "        [  1.40247253e-01,   0.00000000e+00],\n",
       "        [  1.29206921e-01,   0.00000000e+00],\n",
       "        [  1.07895024e-01,   0.00000000e+00],\n",
       "        [  2.44058321e-01,   0.00000000e+00],\n",
       "        [  4.80888580e-01,   0.00000000e+00],\n",
       "        [  5.06484441e-01,   0.00000000e+00],\n",
       "        [  2.73762376e-01,  -2.57425743e-02],\n",
       "        [  2.43953287e-01,  -2.17059891e-02],\n",
       "        [  3.86981267e-01,  -1.00100435e-01],\n",
       "        [  3.86981267e-01,  -1.00100435e-01],\n",
       "        [  2.25936759e-01,   0.00000000e+00],\n",
       "        [  1.21916239e-01,   0.00000000e+00],\n",
       "        [ -1.29688651e-01,   0.00000000e+00],\n",
       "        [ -1.30308012e-01,  -2.95974312e-04],\n",
       "        [ -1.35860403e-01,  -5.36401271e-04],\n",
       "        [ -4.01519262e-02,  -2.56248730e-01],\n",
       "        [ -5.98650148e-04,  -3.32640333e-01]])"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def get_speed(track, bias=1):\n",
    "    '''\n",
    "    此函数用于获得速度，至少需要三个点，具体公式是解三个点用速度表示时的方程组\n",
    "    这里默认时间偏差+1，是为了防止同一时间出现了两个点造成的inf\n",
    "    固有偏差影响，视平均时间差不同，速度总体上大约缩小了0.25%\n",
    "    '''\n",
    "    assert track.shape[0] >= 3\n",
    "    B_sub_A = track[1:-1] - track[:-2]\n",
    "    C_sub_A = track[2:] - track[:-2]\n",
    "    C_sub_B = track[2:] - track[1:-1]\n",
    "    sp = np.multiply((C_sub_B[:,:-1] / np.square(C_sub_B[:,-1]+bias).reshape((-1,1)) + \\\n",
    "        B_sub_A[:,:-1] / np.square(B_sub_A[:,-1]+bias).reshape((-1,1))) ,\n",
    "        (B_sub_A[:,-1].reshape((-1,1))+bias))\n",
    "    sp = np.multiply(sp, (C_sub_B[:,-1].reshape((-1,1))+bias))\n",
    "    sp = sp/(C_sub_A[:,-1].reshape((-1,1))+bias)\n",
    "    return sp\n",
    "\n",
    "index, start, track, des, flag = get_track(people[0])\n",
    "get_speed(track)\n",
    "# print(track)\n",
    "# get_speed(np.matrix([[0,0,0],[1,1,1],[3,3,2]]), bias=0) #测试数据，函数正确的话加速度应为1.5\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "matrix([[  353. ,  2607. ,   349. ],\n",
       "        [  360. ,  2607. ,   362.5],\n",
       "        [  367. ,  2607. ,   376. ],\n",
       "        [  377.5,  2613.5,   397. ],\n",
       "        [  388. ,  2620. ,   418. ],\n",
       "        [  402. ,  2620. ,   430. ],\n",
       "        [  416. ,  2620. ,   442. ],\n",
       "        [  458. ,  2620. ,   467.5],\n",
       "        [  500. ,  2620. ,   493. ],\n",
       "        [  542. ,  2620. ,   520. ],\n",
       "        [  584. ,  2620. ,   547. ],\n",
       "        [  629.5,  2620. ,   569.5],\n",
       "        [  675. ,  2620. ,   592. ],\n",
       "        [  699.5,  2620. ,   617.5],\n",
       "        [  724. ,  2620. ,   643. ],\n",
       "        [  752. ,  2620. ,   668.5],\n",
       "        [  780. ,  2620. ,   694. ],\n",
       "        [  801. ,  2620. ,   718. ],\n",
       "        [  822. ,  2620. ,   742. ],\n",
       "        [  836. ,  2626.5,   767.5],\n",
       "        [  850. ,  2633. ,   793. ],\n",
       "        [  867.5,  2633. ,   818.5],\n",
       "        [  885. ,  2633. ,   844. ],\n",
       "        [  909.5,  2633. ,   869.5],\n",
       "        [  934. ,  2633. ,   895. ],\n",
       "        [  958.5,  2633. ,   920.5],\n",
       "        [  983. ,  2633. ,   946. ],\n",
       "        [ 1021.5,  2633. ,   976. ],\n",
       "        [ 1060. ,  2633. ,  1006. ],\n",
       "        [ 1102. ,  2633. ,  1034.5],\n",
       "        [ 1144. ,  2633. ,  1063. ],\n",
       "        [ 1189.5,  2633. ,  1078. ],\n",
       "        [ 1235. ,  2633. ,  1093. ],\n",
       "        [ 1259.5,  2633. ,  1118.5],\n",
       "        [ 1284. ,  2633. ,  1144. ],\n",
       "        [ 1298. ,  2633. ,  1177. ],\n",
       "        [ 1312. ,  2633. ,  1210. ],\n",
       "        [ 1319. ,  2633. ,  1226.5],\n",
       "        [ 1326. ,  2633. ,  1243. ],\n",
       "        [ 1329.5,  2633. ,  1298.5],\n",
       "        [ 1333. ,  2633. ,  1354. ],\n",
       "        [ 1343.5,  2633. ,  1381. ],\n",
       "        [ 1354. ,  2633. ,  1408. ],\n",
       "        [ 1364.5,  2639.5,  1429. ],\n",
       "        [ 1375. ,  2646. ,  1450. ],\n",
       "        [ 1413.5,  2652.5,  1471. ],\n",
       "        [ 1452. ,  2659. ,  1492. ],\n",
       "        [ 1462.5,  2665.5,  1517.5],\n",
       "        [ 1473. ,  2672. ,  1543. ],\n",
       "        [ 1476.5,  2672. ,  1748.5],\n",
       "        [ 1480. ,  2672. ,  1954. ],\n",
       "        [ 1483.5,  2672. ,  2002. ],\n",
       "        [ 1487. ,  2672. ,  2050. ],\n",
       "        [ 1490.5,  2672. ,  2141.5],\n",
       "        [ 1494. ,  2672. ,  2233. ],\n",
       "        [ 1497.5,  2672. ,  2239. ],\n",
       "        [ 1501. ,  2672. ,  2245. ],\n",
       "        [ 1508. ,  2672. ,  2269. ],\n",
       "        [ 1515. ,  2672. ,  2293. ],\n",
       "        [ 1515. ,  2665.5,  2320. ],\n",
       "        [ 1515. ,  2659. ,  2347. ],\n",
       "        [ 1518.5,  2659. ,  2450.5],\n",
       "        [ 1522. ,  2659. ,  2554. ],\n",
       "        [ 1525.5,  2659. ,  2638. ],\n",
       "        [ 1529. ,  2659. ,  2722. ],\n",
       "        [ 1536. ,  2659. ,  2747.5],\n",
       "        [ 1543. ,  2659. ,  2773. ],\n",
       "        [ 1546.5,  2659. ,  2783.5],\n",
       "        [ 1550. ,  2659. ,  2794. ],\n",
       "        [ 1557. ,  2659. ,  2818. ],\n",
       "        [ 1564. ,  2659. ,  2842. ],\n",
       "        [ 1571. ,  2659. ,  2867.5],\n",
       "        [ 1578. ,  2659. ,  2893. ],\n",
       "        [ 1585. ,  2659. ,  2918.5],\n",
       "        [ 1592. ,  2659. ,  2944. ],\n",
       "        [ 1595.5,  2659. ,  2968. ],\n",
       "        [ 1599. ,  2659. ,  2992. ],\n",
       "        [ 1602.5,  2659. ,  3017.5],\n",
       "        [ 1606. ,  2659. ,  3043. ],\n",
       "        [ 1609.5,  2659. ,  3071.5],\n",
       "        [ 1613. ,  2659. ,  3100. ],\n",
       "        [ 1616.5,  2659. ,  3139. ],\n",
       "        [ 1620. ,  2659. ,  3178. ],\n",
       "        [ 1627. ,  2659. ,  3199. ],\n",
       "        [ 1634. ,  2659. ,  3220. ],\n",
       "        [ 1641. ,  2659. ,  3232. ],\n",
       "        [ 1648. ,  2659. ,  3244. ],\n",
       "        [ 1658.5,  2659. ,  3272.5],\n",
       "        [ 1669. ,  2659. ,  3301. ],\n",
       "        [ 1672.5,  2652.5,  3373. ],\n",
       "        [ 1676. ,  2646. ,  3445. ],\n",
       "        [ 1683. ,  2646. ,  3467.5],\n",
       "        [ 1690. ,  2646. ,  3490. ],\n",
       "        [ 1704. ,  2639.5,  3518.5],\n",
       "        [ 1718. ,  2633. ,  3547. ],\n",
       "        [ 1725. ,  2633. ,  3569.5],\n",
       "        [ 1732. ,  2633. ,  3592. ],\n",
       "        [ 1735.5,  2633. ,  3619. ],\n",
       "        [ 1739. ,  2633. ,  3646. ],\n",
       "        [ 1735.5,  2633. ,  4288. ],\n",
       "        [ 1732. ,  2633. ,  4930. ],\n",
       "        [ 1728.5,  2633. ,  4955.5],\n",
       "        [ 1725. ,  2633. ,  4981. ],\n",
       "        [ 1721.5,  2626.5,  5723.5],\n",
       "        [ 1718. ,  2620. ,  6466. ],\n",
       "        [ 1711. ,  2620. ,  6514. ],\n",
       "        [ 1704. ,  2620. ,  6562. ],\n",
       "        [ 1704. ,  2613.5,  6580. ],\n",
       "        [ 1704. ,  2607. ,  6598. ],\n",
       "        [ 1700.5,  2607. ,  6917.5],\n",
       "        [ 1697. ,  2607. ,  7237. ]])"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def expand_track(track, max_sample=200):\n",
    "    '''\n",
    "    此函数用来扩充取样点，具体扩充方法是每两个点取似然中间点，符合x和y轴上速度相似的设定\n",
    "    为了扩充至少需要三个点，来计算加速度。\n",
    "    我认为，如果只有两个点，那无论什么模型都很难判断轨迹是否是robot，所以不如当作脏数据清洗掉\n",
    "    输入为track的matrix,返回扩充后的matrix，shape'=(shape[0]*2-1, shape[1])\n",
    "    '''\n",
    "    max_exp = max_sample - track.shape[0]#最大扩充点数，使得扩充后的样点数量不大于max_sample\n",
    "    v = get_speed(track)\n",
    "    a = get_acceleration(track)\n",
    "    t = track[1:,-1] - track[:-1,-1]\n",
    "    #加速度相等的推测点算法，结果轨迹不正常\n",
    "#     Xb = track[1:-1,:-1] - np.multiply(v,t[:-1]/2) - np.multiply(a, np.square(t[:-1]/2))*0.5\n",
    "#     Xb_e = track[-2,:-1] + np.multiply(v[-1],t[-1]/2) + np.multiply(a[-1], np.square(t[-1]/2))*0.5\n",
    "#     to_exp = np.vstack([Xb, Xb_e])\n",
    "    #速度相等的推测点算法，结果能接受\n",
    "#     Xb = track[1:-1,:-1] + np.multiply(v,t[1:]/2)\n",
    "#     Xb_s = track[1,:-1] - np.multiply(v[0],t[0]/2)\n",
    "#     to_exp = np.vstack([Xb_s, Xb])\n",
    "    #直接取中间点，对轨迹影响最小\n",
    "    Xb = (track[:-1] + track[1:])/2\n",
    "    to_exp = Xb\n",
    "    expanded =  np.vstack([track, to_exp]) if to_exp.shape[0] <= max_exp else np.vstack([track, to_exp[:max_exp]])\n",
    "    return sort_track(expanded)\n",
    "\n",
    "index, start, track, des, flag = get_track(people[0])\n",
    "# print(track)\n",
    "expand_track(track)\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[False]\n",
      " [False]\n",
      " [ True]\n",
      " [ True]]\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "matrix([[ 6,  7],\n",
       "        [ 9, 10]])"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "a = np.arange(12).reshape((-1,3))\n",
    "a = np.matrix(a)\n",
    "print(np.argsort(a[:,2],axis=0)>1)\n",
    "a[np.ix_((np.argsort(a[:,2],axis=0)>1).flat,[0,1])]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def drop_track(track, max_sample=100):\n",
    "    '''\n",
    "    此函数用于将轨迹中超过数量的点抛弃\n",
    "    抛弃规则是时间距离越近越容易被抛弃\n",
    "    '''\n",
    "    t = track[1:,-1] - track[:-1, -1]\n",
    "    to_drop = track.shape[0] - max_sample\n",
    "    assert to_drop >= 0\n",
    "    return track[np.ix_((np.argsort(track[:,-1],axis=0)>to_drop-1).flat,[0,1])]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 269 2503]\n",
      " [ 444 2594]\n",
      " [ 696 2594]\n",
      " [ 773 2490]\n",
      " [ 850 2594]\n",
      " [ 885 2594]\n",
      " [ 899 2581]\n",
      " [ 906 2594]\n",
      " [ 976 2594]]\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "0.67096774193548392"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def find_stab(track, stride=100, height=30, max_t=50):\n",
    "    pace = stride / 2\n",
    "    stabs = []\n",
    "    trackT = track.copy()\n",
    "    trackT[1:,2] = trackT[1:,0]>=trackT[:-1,0]\n",
    "    trackT[0,2] = 1\n",
    "#     print(track)\n",
    "    for i in range(1, track.shape[0]):\n",
    "        trend = track[i,0] >= track[i-1,0]\n",
    "        Range = np.all([track[:,0]>=track[i,0]-pace,\n",
    "                        track[:,0]<=track[i,0]+pace,\n",
    "                        trackT[:,2]==trend ,\n",
    "                        np.abs(track[:,2] - track[i,2]) < max_t,\n",
    "                       ],axis=0)\n",
    "        r = track[np.ix_(Range.flat, [0,1])]\n",
    "        s = np.abs(track[i,1] - (r[0,1]+r[-1,1])/2)\n",
    "        if s > height:\n",
    "            stabs.append(np.hstack([track[i,0:2], np.matrix([s])]))\n",
    "#         print(r[0,1], r[-1,1], (r[0,1]+r[-1,1])/2, track[i,1] - (r[0,1]+r[-1,1])/2)\n",
    "    if len(stabs) > 0:\n",
    "        return np.vstack([*stabs])\n",
    "    else:\n",
    "        return np.zeros((1,3))#固有\n",
    "    \n",
    "def find_stab_2(track,height=0):\n",
    "    stabs = []\n",
    "    track[1:,2][track[1:,1]>track[:-1,1]] = 1\n",
    "    track[1:,2][track[1:,1]<track[:-1,1]] = -1\n",
    "    track[1:,2][track[1:,1]==track[:-1,1]] = 0\n",
    "    track[0,2] = track[1,2]\n",
    "    track[:-1,2] = np.abs(track[1:,2] - track[:-1,2])\n",
    "    track[-1,2] = 0\n",
    "#     print(track)\n",
    "    stabs = track[np.ix_((track[:,2]!=0).flat, [0,1])]\n",
    "    if stabs.shape[0] == 0:\n",
    "        return None\n",
    "    condx = np.all([np.abs(stabs[1:,1]-stabs[:-1,1])>=height,], axis=0)\n",
    "    return np.vstack([track[0,0:2],stabs[0],stabs[1:][np.ix_(condx.flat, [0,1])], track[-1,0:2]])\n",
    "    \n",
    "index, start, track, des, flag = get_track(machine[315])\n",
    "stabs = find_stab_2(track)\n",
    "print(stabs)\n",
    "np.abs((stabs[1:-1,1] - (stabs[2:,1]+stabs[:-2,1])/2)/(stabs[2:,0]-stabs[:-2,0]+1)).max()\n",
    "\n",
    "# track[1:,2][track[1:,1]>track[:-1,1]] = 1\n",
    "# track[1:,2][track[1:,1]<track[:-1,1]] = -1\n",
    "# track[1:,2][track[1:,1]==track[:-1,1]] = 0\n",
    "# track[0,2] = track[1,2]\n",
    "# track[:-1,2] = np.abs(track[1:,2] - track[:-1,2])\n",
    "# print(track)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "124 2121\n",
      "160 2541\n",
      "366 2177\n",
      "387 2198\n",
      "392 2058\n",
      "397 2086\n",
      "404 2205\n",
      "465 2086\n",
      "466 2128\n",
      "496 2919\n",
      "516 2100\n",
      "553 2471\n",
      "566 2128\n",
      "576 2002\n",
      "602 2478\n",
      "677 2023\n",
      "704 3241\n",
      "883 6258\n",
      "889 2002\n",
      "897 2345\n",
      "924 2191\n",
      "967 2303\n",
      "972 2870\n",
      "984 2128\n",
      "1078 2401\n",
      "1136 3010\n",
      "1157 2373\n",
      "1158 2037\n",
      "1216 2352\n",
      "1223 2205\n",
      "1375 2002\n",
      "1428 2121\n",
      "1444 4459\n",
      "1470 2002\n",
      "1488 2009\n",
      "1501 3248\n",
      "1525 2954\n",
      "1581 2030\n",
      "1627 2128\n",
      "1683 2023\n",
      "1716 2086\n",
      "1723 4655\n",
      "1726 2303\n",
      "1794 2058\n",
      "1813 2324\n",
      "1832 2002\n",
      "1903 2632\n",
      "1917 2205\n",
      "1976 2534\n",
      "2018 2569\n",
      "2047 2303\n",
      "2074 2058\n",
      "2081 2933\n",
      "2129 2044\n",
      "2151 2037\n",
      "2198 2520\n",
      "2204 3136\n",
      "2209 2254\n",
      "2236 2121\n",
      "2264 2527\n",
      "2330 2170\n",
      "2332 2205\n",
      "2373 2226\n",
      "2380 2093\n",
      "2450 2051\n",
      "2459 2135\n",
      "2510 4305\n",
      "2530 3178\n",
      "2539 2002\n",
      "2543 2107\n",
      "2545 2415\n"
     ]
    }
   ],
   "source": [
    "import random\n",
    "# random.shuffle(people)\n",
    "# random.shuffle(machine)\n",
    "def show_data_attr(data, limit=100):\n",
    "    for i,s in enumerate(data):\n",
    "        if i >= limit:\n",
    "            break\n",
    "        index, start, track, des, flag = get_track(s)\n",
    "        track = track - track[0]\n",
    "        if  track[:,0].max()>2000:\n",
    "            print(index, track[:,0].max())\n",
    "#         track = np.vstack([track[0], track, track[-1]])#添加固定偏差来保证数据有效\n",
    "#         delta = track[1:] - track[:-1]\n",
    "#         delta[:,2]+=1#固定偏差用来消除0\n",
    "\n",
    "#         angular = np.arccos(delta[:,0]/np.sqrt(np.square(delta[:,0])+np.square(delta[:,1])))\n",
    "#         angular = np.nan_to_num(angular)\n",
    "#         angularDeltaTime = np.hstack((angular, delta[:,2]))#x1角度,x2为时间差\n",
    "#         timeAngular = np.hstack((track[1:,2],angular))#x1为时间，x2为角度\n",
    "\n",
    "#         angularVelocity = (angular[1:] - angular[:-1])/delta[1:,2]#角速度\n",
    "#         angularVelocityDeltaTime = np.hstack((np.nan_to_num(angularVelocity), delta[1:,2]))#x1角度,x2为时间差\n",
    "#         timeAngularVelocity = np.hstack((track[2:,2],np.nan_to_num(angularVelocity)))#x1为时间，x2为角度\n",
    "\n",
    "#         speed = np.sqrt(np.square(delta[:,0])+np.square(delta[:,1]))/delta[:,2]\n",
    "#         timeSpeed = np.hstack((track[1:,2],np.nan_to_num(speed)))#x1为时间，x2为速度\n",
    "#         speedDeltaTime = np.hstack((np.nan_to_num(speed), delta[:,2]))#x1为速度，x2为时间差\n",
    "\n",
    "#         acceleration = (speed[1:,0] - speed[:-1,0])/delta[1:,2]#加速度\n",
    "#         accelerationDeltaTime = np.hstack((np.nan_to_num(acceleration), delta[1:,2]))\n",
    "#         timeAcceleration = np.hstack((track[2:,2], np.nan_to_num(acceleration)))\n",
    "\n",
    "#         trackGrad = (track[1:,1] - track[:-1,1])/(track[1:,0] - track[:-1,0])\n",
    "        \n",
    "        \n",
    "    #     trackHessen = (trackGrad[1:,1] - trackGrad[:-1,1])/(trackGrad[1:,0] - trackGrad[:-1,0])\n",
    "    #     print('speed-max:%f'%speed[:,0].max()\n",
    "    #           ,\"mean:%f\"%speed[:,0].mean()\n",
    "    #           ,\"std:%f\"%speed[:,0].std()\n",
    "    #           ,\"%%25time speed:%f\"%(get_distr(timeSpeed,part=4)[0][0])\n",
    "    #           ,\"%%25speed time:%f\"%(get_distr(speedDeltaTime,part=4)[0][0])\n",
    "    #          )    \n",
    "\n",
    "    #     print('acceleration-max:%f'%acceleration[:,0].max()\n",
    "    #           ,\"mean:%f\"%acceleration[:,0].mean()\n",
    "    #           ,\"std:%f\"%acceleration[:,0].std()\n",
    "    #           ,\"%%25time acceleration:%f\"%(get_distr(timeAcceleration,part=4)[0][0])\n",
    "    #           ,\"%%25acceleration time:%f\"%(get_distr(accelerationDeltaTime,part=4)[0][0])\n",
    "    #          )\n",
    "\n",
    "    #     print('angular-max:%f'%angular[:,0].max()\n",
    "    #       ,\"mean:%f\"%angular[:,0].mean()\n",
    "    #       ,\"std:%f\"%angular[:,0].std()\n",
    "    #       ,\"%%25time angular:%f\"%(get_distr(timeAngular,part=4)[0][0])\n",
    "    #       ,\"longest angular time:%f\"%(get_distr(angularDeltaTime,part=16)[0].max())\n",
    "    #      )\n",
    "\n",
    "    #     print('angularVelocity-max:%f'%angularVelocity[:,0].max()\n",
    "    #       ,\"mean:%f\"%np.abs(angularVelocity[:,0]).mean()\n",
    "    #       ,\"std:%f\"%angularVelocity[:,0].std()\n",
    "    #       ,\"%%25 time angularVelocity:%f\"%(get_distr(np.abs(timeAngularVelocity),part=4)[0][0])\n",
    "    #       ,\"%%25 angularVelocity time:%f\"%(get_distr(np.abs(angularVelocityDeltaTime),part=8)[0][0])\n",
    "    #      )\n",
    "\n",
    "    #     print('track-range-x:%f -y:%f'%((track[:,0].max() - track[:,0].min()), (track[:,1].max() - track[:,1].min()))\n",
    "    #       ,\"mean-x:%f -y:%f\"%((track[:,0].mean()), (track[:,1].mean()))\n",
    "    #       ,\"std-x:%f -y%f\"%((track[:,0].std()), (track[:,1].std()))\n",
    "    #       ,\"hessen:%f\"%(1)\n",
    "    # #     ,\"%%25 time track:%f\"%(get_distr(np.abs(track[:,0:]),part=4)[0][0])\n",
    "    # #     ,\"%%25 track time:%f\"%(get_distr(np.abs(angularVelocityDeltaTime),part=8)[0][0])\n",
    "    #      )\n",
    "\n",
    "    #     print(trackHessen)\n",
    "#         stabs = find_stab_2(track,height=0)\n",
    "#         if stabs is None:\n",
    "# #             print(index,0)\n",
    "#             pass\n",
    "#         elif stabs.shape[0] >= 3:\n",
    "#             winX = track[:,0].max() - track[:,0].min()\n",
    "#             winY = track[:,1].max() - track[:,1].min()\n",
    "#             stabs = np.min(\n",
    "#                 np.abs(\n",
    "#                 np.hstack([stabs[1:-1,1]-stabs[2:,1],stabs[1:-1,1]-stabs[:-2,1]])\n",
    "#                 ),axis=1)/winY*winX/(stabs[2:,0]-stabs[:-2,0]+10)\n",
    "#             t = stabs.max()\n",
    "#             if stabs.max() > 1.5:\n",
    "#                 print(index,stabs.max(), np.count_nonzero(stabs>1.5))\n",
    "#         else:\n",
    "#             print(index,0)\n",
    "    #     elif stabs.shape[0] == 1:\n",
    "    #         print(1, 0, 0)\n",
    "    #     else:\n",
    "    #         print(stabs.shape[0], np.abs(stabs[1:,1]-stabs[:-1,1]).mean(), np.abs(stabs[1:,1]-stabs[:-1,1]).max())\n",
    "    #         print(0)\n",
    "    #         if stabs.shape[0] >= 3:\n",
    "    #             print(((stabs[1:-1,1] - (stabs[2:,1]+stabs[:-2,1])/2)/(stabs[2:,0]-stabs[:-2,0])).max())\n",
    "    #     print('-'*6)\n",
    "    #     print(get_distr(speed_grad))\n",
    "\n",
    "show_data_attr(people[:],limit=10000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/miller/anaconda3/lib/python3.6/site-packages/ipykernel/__main__.py:14: RuntimeWarning: invalid value encountered in true_divide\n",
      "/home/miller/anaconda3/lib/python3.6/site-packages/ipykernel/__main__.py:31: RuntimeWarning: invalid value encountered in true_divide\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYEAAAD8CAYAAACRkhiPAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzsnXmYFNXV/z+HHWTfcQYYxkFkRwQEVNSIgkTBFUET4otL\nNKgx8dWf0ddofDUxebOLSojGiFFQowRcwAUENxABQWAEBhyWGVlEdtmH+/vjdNE1Pd0z1dPV+/08\nTz/VXV3LmTvV9a1z7rnnijEGi8VisWQnNZJtgMVisViShxUBi8ViyWKsCFgsFksWY0XAYrFYshgr\nAhaLxZLFWBGwWCyWLMaKgMVisWQxVgQsFosli7EiYLFYLFlMrWQbUBUtW7Y0eXl5yTYjrixZsmSH\nMaaV1+2zoU0gunaxbRKedG6XPXv2sHnzZgBatmxJ27ZtK2yzb98+1q5dWwasAXYYY86t6rjp3CbR\n4PVaSXkRyMvLY/Hixck2I66IyMZots+GNoHo2sW2SXjStV3Kyso49dRTKSwsJDc3l/79+zNlyhS6\ndet2Ypvdu3czePBggEJjTC8Rae3l2OnaJtHi9Vqx4SCLxZJyLFq0iIKCAvLz86lTpw5jxoxhxowZ\n5bZ58cUXueKKKwCOABhjtifB1LTHioDFYkk5SktLad++/YnPubm5lJaWlttm7dq17Nq1C6CLiCwR\nkXGRjiciN4vIYhFZ/M0338TL7LQk5cNB6cDX+77m8U8fZ9qqaRw+drjC97Vr1qZ3m94Mbj+Yri27\nUkNUe/Ob5dO9dXff7DAGRo6E8ePh8st9O6wlDsyfD7/5DbzxBtRK4K9w7174/vfhb38DV2QlLTl2\n7BhLliwBKAKGAQtEZKExZm3otsaYycBkgH79+lUonXzjjXDGGXDrrXE2OgWxIhADx81x7px9J5MW\nT6LMlDGi8wjaNWxXYbsDRw+w+OvFvL729XLrfz7w5/xh2B98s+fQIb2pHDhgRSDVef99ePtt2LoV\ncnMTd97CQvjoI1iwILVFICcn50SnMEBJSQk5OTnltsnNzaVFixZ8+umnx40xO0TkA6A3UEEEquLV\nV+G99+CWW0AkVuvTCysCMXD/nPt5fNHj3NT3Ju49+17ym+VXuv2OAzvYuDvYV9P6JE/9WJ7Zu1eX\n8+fDrl3QrJmvh7f4yLff6jLRIrBtmy737EncOatD//79KSoqori4mJycHKZNm8aLL75YbptRo0Zx\n2223ASAiDYAzgT9Fe66yMm2P3bth8WLo39+PvyB9sCJQTZ5b9hyPffwYPz7jxzz1/acQD48PLRu0\npGWDlnGzyRGBsjKYPRvGjo3bqSwx4oiAc1NOFM75du9O7HmjpVatWkycOJFhw4ZRVlbG+PHj6d69\nO5MmTQLglltuoWvXrgwfPpy5c+d2BxYBTxtjVkZ7rj17NJQK8PLL2ScCtmO4Gny48UNuev0mLuh0\nAY9f/LgnAUgE7qe7mTOTZ0c0jB8/ntatW9OjR48T6x566CHH9e8mIstEZETSDIwTyRKBrVt1meoi\nADBixAjWrl3L+vXruf/++wG9+d9yyy0ntrn77rsBVhljehhj/lyd82jfsvbNvPJKUBCyBSsCUVJ2\nvIxrX7uWTs068crVr1C7Zu1km3QCxxPIz4dZs+Do0eTa44Xrr7+e2bNnV1j/s5/9DDT/u48x5q2E\nGxZnku0JpHo4KJHs3KnLUaNg40b47LPk2pNorAhEyUebPqJkbwmPnP8IzeqnVtDdEYEf/EB/5B9+\nmFx7vDBkyBCaN2+ebDMSzo4durSeQPJxPIHrr4fatTUklE1YEYiSVwpfoX6t+lzc+eJkm1IB5+nu\nqqugbt30CQmF4/HHHwcNB/1DRFJLbX3A3TGcSNKlTyCROJ5Afj5cdFH2hYSsCETBcXOcV798lYs7\nX0zDOg2TbU4FHE+gXTsYOlRFIB0v5ltvvZWvvvoKoBDYAoTNo03XAUBHjsD+/freegLJx/EEmjeH\nq6+GTZtg0aLk2pRIrAhEwcebPmbr/q1c3e3qZJsSFkcEGjfWQWPFxbBqVXJtqg5t2rShZs2azse/\nAwPCbWeMmWyM6WeM6deqleeaaknH8QLA9gmkAo4n0KyZ9gvUrq3eQLZgRSAKnl32LHVr1uX7nb+f\nbFPCsmcP1KsHderAJZfoutdfr3yfVGTLli3uj5cDUaf9pTKOCLRqlVgR+O479UBErCfgZtcuaNBA\nQ6hNm8KwYdovkI5edHWwIuCRJz97kmeXPcuPz/gxjeo2SrY5Ydm7V70AgJNP1nznVO8XGDt2LIMG\nDWLNmjXk5ubyzDPPcM8999CzZ0+AbsD5wM+Sa6W/OCLQrZs+hR45kpjzOoLToYNeK8ePJ+a8qU7o\nwMqrr4bNm+HTT8tvt38/PPlksFM/U7Ai4IHpX07ntrdu45JTL/G1zIPfuEUA4NJL9UJ2dz4ePQo3\n3ADnnw9ff514G0OZOnUqW7Zs4ejRo5SUlHDDDTfw/PPPs2LFCtAU0ZHGmC1VHSdabroJ7rnH76N6\nwy0CANsTVPvSuQ66dNGnXCd8mO3s3Kn9AQ4jR6o37Q4Jbd4MZ58NEyZA7946Kj9TsCJQBSu3r+Ta\n165lQM4Apl05jVo1UneQ9Z490KRJ8PPIkfpjf/NN/Xz4MFxzDfzjH7BwIQwYAEuXJsfWZHLsGLz4\nIjz9tI6uTjShIpCokJBzni5ddGn7BZRQT8AJCb3yinpLixfDmWdqH9sTT0DDhvC978GvfpWc68dv\nrAhUwZOfPQnAzLEzOanOSUm2pnJCPYFevdT1nzlTi8tdcQVMnw5//auKQI0acM458NprybM5GaxY\noUX2du1KzsAgRwS6BwrIJkoE3J4A2H4Bh1BPAIIhoQcegCFDtL/gk0/gJz9RUbjuOnjoIbjgAgip\ncJ12WBGohMPHDvPSqpe47LTLfC/2Fg9CRUBEvYF339WO4lmztITw7berS7tokQrFlVdqWeNs6Qhb\nsCD4Psxg5bjz7bfage/McJhIT0AEOnfWz1YElHDFFp2Q0K9/DX36aFjVEe1GjWDKFPjnP/Uholev\n1O97qwwrApXwVtFb7Dy4k3G9Is5VkVKEhoNAL+aDB7V08bPPws03B79r21bXjx0L992nIyYPV5wO\nIeNYuBDatIGBA5MnAi1aqA2QWE+gZUt9gQ0HOezcWVEEmjSBu+/W0tJz50LrMM+AP/qRhlM7dtTU\n0ttvV4873bAiUAnPf/E8bU5qw4WnXJhsUzwR6gkAnHsuXHstvPSSXrSh1KsHL7wADz+sTzdDh0Ia\njbuqFgsWwKBBcPHF6g258/YTgSMCDRpofDlRo4a3bVPhadpUP1tPQB96DhyoGA4CeOQReOop/Y1E\noksXvZ5+9jOYOBEGD04/IbAiEIFvD3zLG2vf4Nqe16Z0Z7CDk+0RKgJ16uhN/qqrIu8rorHPl14K\ndoK55vPIKHbsgHXrVASGD9d2e/fdxNrgiADoTTmRnkDbtkFv0YpAcLRwLHNv1K0Lf/wjTJoEn3+e\nfskWnkVARGqKyOci8kbgc3MReVdEigLLZq5tfyEi60RkjYgMc60/Q0RWBL77q6RKDeYwvLzqZY4e\nP8q43rGFgmbPnk2XLl0oKCjgsccei7idiPQXkWMiUsntOjIHDmimQmg4KBpGj4Z58zQL4oUXqn+c\nVGbhQl0OHKjTCbZooX0licQtAm3bJrZPoE0bKwJu3CUjYmVY4E6XbqP0o/EEfgp86fp8LzDHGNMZ\nmBP4jIh0A8YA3YHhwJMi4tQAeAq4CegceA2Pyfo4MuWLKfRs3ZPebXpX+xhlZWVMmDCBWbNmUVhY\nyNSpUyksLIy0+W+Bd6p7LnfJiFg480yNf65fH9txUpUFC6BmTejXT5cXXaTTPPo5cKoq4d+xw7B6\n9R0UFBSwYkUvNmxYWuW+sc6xYEzQE6hVS8NQqd4n4PUBCmhQ3Qcod8mIWOnQQdt1ZZqNb/ckAiKS\nC3wfeNq1ehTwXOD9c8BlrvXTjDGHjTHFwDpggIi0AxobYxYaYwwwxbVPSvHM0mdYWLKQH/b6YUwT\nxixatIiCggLy8/OpU6cOY8aMYcaMGeE2bQ28ClR72JBfIgBQUKAhk0xk4ULNjGrQQD8PH65PyMuX\n+3P8qoTfGNi5cxYHDhRRVFTEBRdMZvPmWz3tG8scC/v3a4KA0xndtGlqewJeH6DKNFE/l2o+QPnp\nCdSooRlEGSkCwJ+BewD381Ib10jOrUDg8iIHcEeUSwLrcgLvQ9dXIFnVIcuOl3H3O3dz4+s3ctEp\nF3Fr/1tjOl5paSnt27c/8Tk3N5fSkKTiwOdmqJdUbZynuljCQQ6ZKgJlZdoRPGhQcN1FF+nSryyh\nqoR/zx44fnwGgwaNQ0To3Xsgx47tZtOmLdE8NESNE3Jq21aXTZqktgh4bYtAyfFdVPMByo8+ATc9\nemSgCIjIJcB2Y8ySSNsEnux9yzJPRnXIw8cOc9UrV/H7Bb9nQv8JvHntmwkpF33nnXcClBhjKg1I\nVCWMfnsCJSX65JhJrFqlT8QDBwbXtW0Lp5/unwhUJfyaiVRKx466jT6Z57JyZWmV+8Yyx4KTgZQu\nnoDXB6jp06cDVPtJ0QkH+TWvUY8eWgYkUaVA/MCLJ3AWMFJENgDTgO+JyL+AbYEQD4Gl82eXAu1d\n++cG1pUG3oeuTzrHjh9j7Ktj+c/q//CX4X9h4oiJvmQE5eTksNmVZlNSUuLEdU+wePFigPxA+16F\n9qFUCJNVJYx+iwCAlvTPHJxBYm5PADQk9MkniYmRO+mozv/JuSk7N6NIeJ1jAcI/MIR6Ak2bpn6f\nQFXceeed/Pa3v/W0baSHKMcT8MODBhUBSK/O4SpFwBjzC2NMrjEmD+3wnWuM+QEwE3Ayz38EOL7a\nTGCMiNQVkU5oB/CiQOhor4gMDGQFjXPtkzSOm+OMnzGe6aun8+dhf+aOM+/w7dj9+/enqKiI4uJi\njhw5wrRp0xg5cmS5bYqLiwFWBNr338BPjDH/ifZcfoeDIPM6hxcu1IFS+fnl1w8frvWE5s6N/RxV\nCb+KQA6HDuk2elMuQSSn0n29zrEA4R8YQj2BVA8HeX2AGjNmDEBPKnmAgsgPUTt3qiAGmzY2HBFI\np5BQLOMEHgMuFJEiYGjgM8aYVcDL6BPLbGCCMcYps/QTtHN5HbAeSHByXnmMMdz+1u08/8XzPHL+\nI/x04E99PX6tWrWYOHEiw4YNo2vXrowePZru3bszadIkJk2a5Ou54uEJZFq/gDNILLSvf9AgLQXg\nR0ioKuFXERjJnDlTMMbw9dcLgSYcPdqu0n1jnWNh2zbtuHRGC6d6OMjrA9SGDRsAVlDNB6hwJSNi\noU0bTf9NJxGIKuZhjJkHzAu8/xa4IMJ2jwKPhlm/GOgRrZHx4tcf/ponFz/J3YPv5r5z7ovLOUaM\nGMGIEeWz+W655Zaw2xpjrq/ueRwRaOTDVAfNmmmMNJNEYOdOWLMGxoUZ9lG7to6UfvNNrS5qjKaM\n5uRotchocAt/WVkZx4+PZ/bs7nz4oYr+4cO3ACPo3PktCgoKqFevAfAs27ZV3Hf8+PF0DxSsueee\ne1i2bBkE51j4cTR2bd2qk9g4T7xOOMiYiqKYCkRqC+fhKdJvKFrClYyIBZE07Bw2xqT064wzzjDx\n4NnPnzU8hPnhaz80x48fj8s5vAIsNjG2yV13GdOggX82DRhgzIUX+ne86hBNu1R1nbz1ljFgzNy5\n4b9/7jn9PvRVWFh9+/ft02OMHh1c98ADxogYc+xYcN1JJxlz553ejlnda2XkSGN69w4e53e/U9v2\n76/uX5c6xPL7GTTImKFD/bVnwgRjGjc2Jsm3Fc/tkpVlI+YWz+XGmTdyYf6FPD3y6ZjGAqQK4UpG\nxEKmpYkuWKDhkP79w3//wx9CURGsXq0ew9KlWnLjySerf841a3TpTub69lt98nTHoBNROmLr1mB/\nANhRww5+ewKgnsDevZphlw5knQh8d+Q7bph5AwXNC3h19KvUqVkn2Sb5gt8icMopsHFj4qY+jDcL\nF0LPnjqiMxwiKnxdusCpp2ra6DXXaLng6s7AtXq1LkNFwCkZ4ZCI0hHbtgUzg8AWkXPYtcu/9FCH\ndOsczjoR+NX8X7Fh9wYmXzo5ZecKrg7hykjHQkGBxsW13y29OX5c68GHpoZWxW236biC55+v3nm9\nikC8PQGnZITbE3BEIN3TRGNBR2/77wk48w5YEUhBlm1dxh8X/JEbT7+RIR2HJNscX4lHOAjiHxIa\nP348rVu3pkePYL7Azp07ufDCCwF6hBYnrA5ffqnt4x4k5oUBAzR8NHFi9SbccURgx45gbaJkiMDe\nvVoyOZwIZLMn8N13mhrstyfQrJkmFVgRSDHKjpdx0+s30bJBS3534e+SbY7vxEsE4j1W4Prrr2d2\nSG7mY489xgUXXACaBnmiOGF1iTRIzAu33aY38+qMIXBEoKwsOCgpkgjs2AFHj0Z/Di84YwTc4SDb\nJ+B/yQg30WYIrVkD+/b5b4cXskYE/rnsnyz+ejF/Hv5nmtWPw389yfgdDmrVStNN4+0JDBkyhOYh\nj2IzZszgR8EZcNzFCavFggX6tOdMqxgNo0drbv3EidHtV1YGa9fqEyEEQ0KRRMC9jd84Xob1BMrj\nd8kINz16QGGhtz61gwe1rPlDD/lvhxeyRgReXPkiXVp04Zru1yTblLjgtyfgdJQmI0No27ZttGvX\nzvnoLk5YDq+FBlev1nlgq5MEVq8e3HSTziG7caP3/TZs0BvAOefo52++0ZDMd99FFoF4hYQq8wSy\nuU8gnp7AWWfpDGOffVb1tgsW6HUxf77/dnghK0Tg2wPfMn/DfC4/7fKMSAcNJdKsYrGSCmmigXzn\nsBF547HQ4JYtwSfy6uCMS/rb37zv44SC3CLg1A0Klx0E8ROBcJ5AvXo6I5b1BOIjAueeqw8dXsKI\nzjbLlmkiQqLJChF4Y+0blJkyruh6RbJNiQv796sQ+BkOAhWB4mLtPPNCdVMpQ2nTps2JMgkhxQmj\nxhj4+ms4+eTq29Ohgz7ZRfOkFk4EduzQ98nwBGrWrHjeVC8dEW/8nEsglObNNc3Yiwi8/74KclmZ\nN8/Bb7JCBKavnk5u41z6ndwv2abEBT/rBrkpKNDOyk2bqt72q680du5HEbaRI0fy3HPOfEXlihNG\nza5dGoYJRpeqR/fummXkNUto9WrtVzn1VP1cmScQbxHYtk1ni6sR8mvPhEqisRBPTwC05Mgnn2jM\nf8cOnec7dO7u/ft1jov/+i/9/Mkn8bGlMjJeBL478h1vr3+by7pclpGhIIifCDj5zl5m3fr4YxWM\njz6K7hxjx45l0KBBrFmzhtzcXJ555hnuvfde3tXZ33vgKk5YHZy6a7F4AgDduqmgeL1Rr14Np52m\nT3hNmmh9+Ugi0LChznTmxO795gc/gHCzM1pPIDjVZjz43ve0X+iTT+CJJ+DVV7Vvyc1HH6mnfcUV\neo19/HF8bKmM2Ivmpzhvr3+bQ8cOcXnXy5NtStzws4y0m169NIywZAlcXkXzLQ1Mk7tiRXTnmDp1\natj1c+bMQURWGmOGRnfE8nz9tS5j9QS6dtVlYWH5DtZIrF4dbLNWrSr3BCC+YwXOOy/8+lQvJx1v\nnIFi8Xo2PPtsFZlZs+Bf/9J1n39efpv339cChmedpa9XXtExJaFeWzzJeE9g+urpNK/fPOMGh7mJ\nlydQv756A84NvjKcizvVBsj46QmAhoSqYscOfZ12mn5OtghEIlmegNNfUl3++ld/QjjxKBnhplEj\nHXD4+OP6v23RouJv6f33dRBjgwYweLD+P2Jtn2jJaBE4WnaUN9a+waWnXurLTGGpSrxEAKBvX/UE\nKouFHz+uIlCjhhZhO3TIfzuqi1+eQLt22r5h5jqvgFM4LpwI1K+vr3DHd2xNFMnoE1iyRL0qZwBf\ndfjpT/0Rr3iUjAjFCQl17Qo33KAPSc7YgT17tD3OP18/Dx6sy0SHhDJaBD7Y+AG7D+3m8tMyNxQE\n8QsHgQ5i2b4dSkv1Zv/00xXT2IqLVYguuEAzHLw8LSeKr7/Wm/dJJ8V2HBH1Brz8bc6TXDgRCOcF\nAOTl6TiE6pSnqC7J8AQcgQztIK0OsbZVvD0BAK1+osLVt6/2mzkPEh98oL8pRwQ6d9bkCqdzePt2\nnejomWfiN5ocMlwE3v3qXWrVqMUF+WHnvskY4ukJnHGGLpcs0Qvyppvg5ZfLb+OEgpwJW1IpJLRl\nS+yhIIeuXb15AqtXa4dwx476uVWrYIgokgh06gQHDsRv1HA4mjRRry2RnptTXtkP8XHqMVUXv2cV\nC8c55+jN/qabNGUUgiGhuXN1vIZT00pEvYF//xtyczVEePHFcOONKgTxIqNFYE7xHAbmDqRhnTh1\n/6cIfs4qFkrv3hrmWbJEL07QkI+bpUu1A+zyy7UGf7Sdw/Hk669jDwU5dOumsd2qJoVfvVpTQ505\nA1q10ie5r76q3BMA9aoSRTIqiToegBcRmD17Nl26dKGgoIDHwqQ3/etfL9CrVy+AbiLyiYj0jsaW\nnTvj7wmIqBDUqKEp1w0bBh+a3n9fb/r16gW3HzdOtzv/fPjjH2HePBWJxx6LnzeQsSKw6+Aulny9\nhAs6ZbYXAPojbtjQv8my3TRooE/An34K/wnM3hoqAp9/rjfIk07SbVNNBPz0BKDqkJCTHurQurUu\n162r3BOAzBcBr55AWVkZEyZMYNasWRQWFjJ16lQKQ9ywjh07MV9H8BUC/wtM9mpHWZn+3fH2BNzU\nqKEPVZ9/rqHB5cuDoSCHK6/U759/Hn72Mx15/D//o6HCF1+Mk13xOWzymbdhHgaTFSIQj5IRbs44\nA955R93nxo3Li4Ax6gn07aufe/ZMnXCQMf6Gg5wMocpCQocP6xO/WwScihZHj1btCSRy/oZkFJFz\nRMAZrRuJRYsWUVBQQH5+PnXq1GHMmDHMmFF+zODAgYNpFryLLwRyvdrhzK+cSBEA/Z0sWxYcVOll\nDusRI6BPH/jNb1S8/CZjRWBO8Rwa1G7AmblnJtuUuBNvEXBu8I0a6cCjdeuCnXJbtmgHlhPv7NlT\nf+hV/cgTgV+jhR06dtTMnso8gXXrNFYdTgQgsgg0bKidgon0BJJRTtprOKi0tJT27duf+Jybm0tp\naWm5bUL6BG4AZnm1I54lIyrj9NO1WNzf/66ec6TpTt2IwH33aaf6q6/6b1NGi8CQjkMyZvrIyvC7\njHQoTufwyJFaIvfAgWA6o9PJ5QhFKk2t59cYAYcaNfTmXpknEJoZBN5EADQklAxPIFHhoCNHgmMh\nfO4YboSKwP+LtG1oxdl4l4yIhPOw9O67Opisdm1v+11xhU59+uij/meQZaQIlO4tZfWO1VkRCoLE\nhIOGDYPbbw/W5HdCQp9/rk8qvQNdcj176nLFCo1j/uMfOon7n/4UP/si4dcYATddu1buCTgi4NQM\nAu8ikJeXnD6BRHkC7nEQVZ0zJyeHza480pKSEnJCSsEePw5ffPEFQEdglDHm20jHC604myxPoFu3\n4I0/tD+gMmrWVG/giy/gjTf8tSkjRWBusQbcrAj4Q/36mh565pnhRaCgIJiZlJurXsnPf643tRtu\ngGnT4P77Ez9zknPT8csTAP0Rb9oUueTv6tXQvn35ejT16gU/V+UJbNwYe+qjVxIdDnLu6V7KVfTv\n35+ioiKKi4s5cuQI06ZNY+TIkeW22bhxE1dccQVAsTFmbTS2JMsTqFMn+KDkpT/Azdix+pvy2xvI\nSBGYUzyHFvVb0LttVBljacvevfENB7lp315z4B0RcHcKg3oFd9wBw4fDX/6iYaG5c7WSopNdlCic\ncJDfngBEHtofmhnk4HgDVXkCR44E7Y43J52kT5iJEgGnU7hHj6r7jGrVqsXEiRMZNmwYXbt2ZfTo\n0XTv3p1JkyYBkwD43e8e5lutxdFRRJaJyGKvtvTtq9dnrueuZP8YOFCvAyc05JXateGuuzRTb9Uq\n/+zJOBEwxjCneA7ndzqfGpL8P6+qXOcXXngBNM95RXVynUFjuvH0BNzUqAGnnKIisHOnPrmGXswP\nP6w3/Dvu0NpDZ52lNziniFai8Gu0sBsnQ0ijEOUxJrIIOGmiVXkCkLh+AZHElo5wRKBnT2/CM2LE\nCNauXcv69eu5//77AbjlllsAneXnT396ml2qJoXGmD7GGM+14k89Va/PRP1u3Pz61zpvQK1qVLL5\n/vd1+cEH/tmT/Lukz6zftZ6SvSV8Ly9KXysOeMl17qS//DXGmJ5EmesMGjrYty+xF3NBgYqAM+il\nqieaGjXguuvgvffiVy45HH6OEXAoKID8fPjDHyoO3vn6aw0TxeIJQOL7BRIZDmrcWCfpOXw49pHK\niQqb+U2TJkHBj5a8PPVe/JyKMuNE4P3i9wE4v1MUvS5xwkuu82CtGuVk/0aV6wzB2HSiwkGg/QLr\n1+soYvDm1l53nf5op02Lr21u/Bwj4FCrlnZyFxZqjXg34TKDHFq1UjF0OmPDkQwRSGQ56ZISvYHF\n0iHtjoWnqwjEgogOIPvgA//6BTJPBDa8T9uGbenSokuyTfGU6xxCVLnOEHTlE+kJdO6sT3Gvv64/\n6kqm9z1B164ah01kSMjPkhFuLr1U+zwefLB8+ecPPtAfqdPx52bYMBg9uvI68fXqqb2JThNNJxFw\ne1/ZKAIAQ4aoRx06cr+6ZJQIGGOYt2Ee5+Wdl3aziInI+VSS6xya5+zg1A1KtCcAOiuSu1O4Kn7w\nA/UeElEv3e/Rwm5EtFPx4EH41a+C53vpJX1KCyeK11wDEebPKUcy0kQT1SewebMmFsQiAgcPBt8n\nsuJqKnHuubr0q18go0Rg7bdr2bJ/C+fnJT8UBN5ynQPUB56mklzn0Dxnh3hWEI2EIwIQXYbDmDH6\nJKx94fHF79HCoZx6qt7YX3pJn06XL9cRnWPGxHbcZAwYS4Qn4AwUc3sC1RlV7haBbPUETj1VK4z6\n1S+QUSIwb8M8gJQRAS+5zpt0FvdTgB9Gm+sMyQkH5eQEKx9GIwLt2umcAy+84Pkprmcgayqq9D+I\nzxiBUK5o43zuAAAgAElEQVS+WjOk5s5VMahZUwuAxUKnTjoOwZl4JN7Eo0/g668rCv2WLfo/z80N\n5ubH6glkqwiIaEjIegJheH/D+5zc6GQKmhck2xSg8lxnzXeGhx9+GHSu5yerc7NLhifglMWF6MJB\noCGh4uKoZpY6P9r0P/C/ZEQ4LrpIB8m9/LJ2eA8dqvV/YuG007RI2Fdf+WNjVTRtqskFx475d8zJ\nk/X/7C657aSH+hkOylYRAA0Jbdrkj9eYMSLg9Aecn3d+SvUHRMp11nxnePrppwGWBW50Ud/snB9S\nZVkn8aBLF73hRTvY5vLLdbDZ9OnxscshHiUjQqlXT+spPf+8/hivuSb2YzqZRYmanc25bpyHCYfC\nQs1nr07c3YmAOjd+97pYO4atCCjnnKNLP6aizBgRWL1jNdu+25YyoaBEEc+pJSvj17+G115T1zQa\nGjXSEaPLl3ve5T0RWSIiN0dznkSIAGhI6OhRHc15uQ+zmHYJJLUlarLxSKUjnnhCS31UZ5S3c/N3\ni4DzPjdXxbNuXSsCsdCli/72/MgQyhgReH+Djg84L++85BqSYHbv1li0n6NivXDqqcGnkWjp2dPz\nxDOrjTF9gIuBCSIyxP1lpIwp0HCQ36OFwzFsmJ7n4ouDT7hVjRI3xnDHHXdQUFBAr169WOqUYgU+\n/ng2tWp14Te/Kb/vzp07uVAnrO0hIu+KiC9VbyI9lS9apMv/+Z/oa9g7N3x3NnRJiT4AOKITrkP6\n0KHyRebCYUVAqVtXQ2vr18d+rCpFQETqicgiEVkuIqtE5FeB9c0DF2NR6EUpIr8QkXUiskZEhrnW\nnxHo6FsnIn8VH+M28zbMo33j9uQ3y/frkGmBU0Y6hSJgVdKzp+Y5e5hP9yiAMWY7MB0Y4P4yUsYU\nxGe0cDjq1dMsjUAXj6dR4rNmzaKoqIiioiImT57MrbfeWm7fAQNm0aVL+X0fe+wxLrjgAoCVwBzg\nXj/sD1dO+vBh9dS6d9ewULTZXOE8gc2by4cOmzatmB304IN6zgMHIh/bikCQ/PwEiQBwGPieMaY3\n0AcYLiID0YtwjjGmM66LUkS6AWOA7sBwtMPTmfjwKeAmoHPgNTz2P0GfrOZvnJ+W4wNiZc+exPcH\nxIpOC1u5N/Ddd99B4PoUkZOAi9AboCcSJQKgsz45YScvo8RnzJjBuHHjEBEGDhzI7t272bJly4l9\n+/bNZ82aOlxzTXDfGTNm8KMf/cg5xHPAZX7YHs4TWL5cQ1wPPqgd/w8+6D1bae/eYP9CaDgoVARC\nPYFZs3Td7NmRj+8WgWXLvNmUqZxyij8JBFWKgFGcwrm1Ay8DjEIvRih/UY4CphljDhtjioF1wAAR\naQc0NsYsNMYYYAo+XcjFu4vZ/t12zmp/lh+HSyt27058f0CsuOcciMQ2HYp7mogsBxYBbxpjKrk9\nlGfLlvj3B4TDyyjxSNs46087TetBNWwY3Hfbtm20C/5BW4E2kWyoLEwWSrg+gc8+0+WZZ2rfz4YN\nmvHjhdAQkPu960+mWbPy59y+PXg9VDZ7lruEdxT9ShnJKafo2ItIZc294qlPQERqisgyYDvwrjHm\nU6CNMcYpeuu+KHOAza7dSwLrcgLvQ9fHzKJSDWAOyBlQxZaZR7xnFYsHbdroqNrQSpyffhqccyA/\nPx+0OmRvY0x3Y8yjXo9vTGI9Ab9xylVHKikdeIiKmLdTWZgslHCewGefadXT9u01Dfbcc+GRR3Ra\nxKpwRKBly6AIHD2qf0tlnsC8ebrs00fLkUQqLufWtGwdMexwyim6jNUb8CQCxpiyQAddLvpU3yPk\n+0ovymiJ5kkGVATq1apHj9Y9qtw209i9O/3CQVCxc3jfPq2zPnp07MfetUvDF8nwBLyMEo+0jbPe\nSRMtLAzu26ZNG7YEVCHgVW/3w97GjbU/yd0nsGgRDBig60XUG9i2Df76V33qXLUK3nwz/BSizo1/\n4MDge/dAMYdQEZgzR2353//Va+Hdd8Pb674d2D4BXSZEBByMMbuB99FY/rbAxRh6UZYCLseP3MC6\nUspXyHTWhzuP5ycZgE9LP+WMdmdQu6bHCTsziHT0BED7BVatCmaeOE+ZTnnqWNgeuBKdGv6JxMso\n8ZEjRzJlyhSMMSxcuJAmTZrQrl27E/seOlRMw4ZH+Pjj4L4jR47kueec6Cs/Asp3NFSTGjX05uvc\nkPft0/RU9wTogwfDJZdoyqiT4nvJJcHa9m6cG/+AAcH+AfdAMQenY9h5mp87Vz2Oiy7S7/797/D2\nWhEI4ngCsXYOe8kOaiUiTQPv6wMXAquBmejFCOUvypnAGBGpKyKd0A7gRYHQ0V4RGRjIChqHDxfy\n0bKjLN2yNCtDQZDensCBA8GnGKfj0evE25Xh3Ci8VDf1Gy+jxEeMGEF+fj4FBQXcdNNNPPnkk+X2\nHT58GGVlXalRQ/cFuPfee3lXH497AEOBirmn1cRdOmLJEr0xu0UAtGDe7bfDb36jhfAmTNARq9+G\nVLoqKdF2d25QpaXlB4o5NG2qo5QPHNDjrFunJUXq1IFRo2DmTHj99Yqptjt2BI9RVqaptmja7Bci\nEuX49fSmWTN9xSoCXua2aQc8F8jwqQG8bIx5Q0QWAC+LyA3ARmA0gDFmlYi8DBQCx4AJxhgn0/gn\nwD/RgmmziLJscjhWbl/JoWOHslIEnAll0tETcHcOd+4czPpIdxEAvcmPGDGi3DpnhDiAiPBE6GQE\nIfv++c/ws59piY1OnaBFixbMmTMHEVlpjBnqp73u0IwzPiBUBPLzVQgcmjXTAWUrVsB55wXXO1lA\nzg2/pCSyJwB63vd1iM+JOXevvBKee66Mm2+ewMcfv0tubi79+/dn5MiRfPNNtxPHKC6exZEjRaBZ\nY7eh2YdnVqsR0hQ/0kS9ZAd9YYw53RjTyxjTwxjzcGD9t8aYC4wxnY0xQ40xO137PGqMOcUY08UY\nM8u1fnHgGKcYY24L9CXExKelnwJwZk5W/e8BdbWNSU9PoHt3jTc7ncOOCNSpE/uxky0CfuBoyKyY\nH5Oqxl1O+rPPVHSqqoHkpPmGdu6HikBpqa5r2LB8fSt3Ebm5c/V/FXB6uPBCaNBgETVrVky1dYeD\n1q/XVFsAY8xCoKkTos4W/EgTTfsRw4tKF9GyQUvymuYl25SEk6ySEX7QoIEWoXM6F50BQpngCfhB\n5876A3/rrfifyx0OcjqFq6JtWxWKSCLgZGaVlAQHirmH8LjLSc+dC+efH5xwp1496Nu3lB072p8o\nbOek0X7zjY6WBdi/v3yqLT5mHKYLp5yi3uLNN2tShZYii46MEIEBOQOybpAYJK94nF/k5QXjxX56\nAjt2aAemc7NIR0S0JMX8+fHvAHXCQdu3a3w+NBQUyb5evcqLwMGD2kfg1Adq1SoYDip/rw5es599\npt87oSCHgQN15LK7Zr7Th+B0+EcTR4g24zBdGDZM52z++9/hlVfgk0+iP0Zai8Dew3sp/KaQASdn\nX38ApLcnAJrC6eTC+90nkM5egEPfvpqSGe+y0k44yBkk5kUEQEVg5cpghpczRsAJBeXkBEUgtNqs\nIwLOwLBQERgxIocaNTafyBIqKSmhcWMnXVbXNWhQPtUWHzMO04Vzzy1fbLA6D1FpLQJLvl6CwXBm\nbvb1B0BQBNLVE2jXTmsIGeN/OCgTfue9e+sy3iNjHRH49FMNyXidI6J3bxVvp2PSXSnUWW7YUHGg\nmHNO0CfX3Nzg/BQO55zTn/r1i5gypZg1azTVtl8/TZd1RKBDB021BQiUstnjGsCaNdSpE7z5Z50I\nOCOF+5/s8dElw3DCQensCRw9qiEEvzuGM0EEunfXm3Jo3N1vmjTRkNPcudCtm3bieiG0czicCHz5\npR47UjjIGE0NDY3mOumyBw8Oo0+frlx11WgaN+4OTGLnzkmB449wRpb3AP6OZh9mJY0a6bI6D1Fp\nLQIfb/6YguYFtGjQItmmJIVMCAeBegPWE6hI/fpasjveIuDckBcu9B4KAhUMt0itXas3c2eAdG5u\nsD8j1BOoU0eTA6BiKMjh+utH8Oyzazl0aD11694f6PC/hXPP1XRbY06k2q40xvQ0xkQ1K18m4WRe\nZZUnUHa8jA83fch5Hc9LtilJI909gbZtdblli3+egDGZIwKgN9TtvhSIiIwjAmVl3jKDHOrV08lN\nvvhC2/3f/9Y5Jpw5HNw3/nAz0DnnPb+SeaDGjYMrrtB5DZxOTycclO0jht1kpSfwxbYv2H1oN+fm\nnZtsU5LGnj36tOhHCCUZOJ7Ali3+eQL79uno40wRgUaNKk796DfuPqVoPAEIZgitWKGhnzFjgt+5\nb/yh4SDnvJ07h//OQQSefFIzvf72N70+nDEGVgSC1K+vy6zyBOZv1Nyxcztmrwika8kIB7cIOJ5A\nrJm+mTBGwE3jxsHKqvHC8STr1g2O5PZKr16apz55ss5wd+WVwe8cETjppPDe6h13wEMPVX2ONm20\nbpExOjahZmB2EisCQerV02V1RMBL2YiUZP7G+eQ3y6d9k0oeIzKcdC0e59Cwob62bAkOFDp6NLZj\nZpoINGoUfxFwHiT69In+JuJ0Dk+erB287qJ97r6BcOL+4x97P8+dd6on0Lx58FqxIhDEGRNTHU86\nLUXguDnOBxs/YFSXUck2JamkuycAwbECjovvjBCtLpkoAk55kHiNh3SuoWhDQRAUgaNHYezY8t81\nbKgPKZWFe7xSrx688456jM5snVYEgjgikDWewMrtK9l5cGfWTSofyp49wZtnuuKMFXDcWSsC5Wnc\nWNvk8OFgG/lNixZaFfSGG6Lft317vdEfPAiXhZkn8NJLow8xRaJzZ106g6OsCARxro2s8QTmb7D9\nAaCeQKdOybYiNtq1g6VLg2EEKwLlcXL29+2LnwiIwMSJ1d93+HDtmAznlT7/fGy2hcOGgyqSdeGg\n+Rvn07FJRzo27ZhsU5JKuvcJQDAcdOqp+tkPEWjQIJimmO7UCvxCU/mGN21aYs/nhMVSuU0SjSMC\n1elTS7vsIGMM8zfOT5tQ0OzZFSfGcBOopt1eRNZFOzFGpojA/v3BXHhHBGbPng06Wcg6EbnX6/Ey\naYyAJTyOJ5Dtcwy7cbzEw4ej3zetRKDseBlPLX6KHQd2pEUoqKysjAkTJjBr1iwKCwuZOnUqhU6v\nVoBZWjC+HjoD283oxBhVcviwTsad7h3DzoAxp0jasWPBdgPWAt2AsSLSLfwRymNFIPOx4aCKOCJw\n6FD0+6aNCMz5ag5nTD6DCW9NYFDuIC7venmyTaqSRYsWUVBQcWIMN4HP3xrF88QY6V4ywsEZK+BM\nU3jsWLDdgCPGmCPANMBTKpgVgczHikBFnHBQxorA7z/5PUOfH8ruQ7uZduU0Ph7/MU3rpf4jcGlp\n+UkvnIkxQrcBjrhWeZoYI93nEnAILSdw7FjFdiOKyUKsCGQ+LQKlwjpmd5dgOfLydFnVjHDhSIuO\n4au6XUXZ8TJ+OvCn1KsVpxSJFEdEbkbDRXTo0IHmzeH3v4d+/ZJsWIycdho8+6xOBGNM+cFGVRHa\nJgD33ZdZN4dBg+D//s97Zc9sYPBg+M9/dEIVi3LDDfpA6B6x7ZW0EIG8pnn8v7P/X7LNiJqcnPKT\nXpSUlJCTk1NhG8A9xCPsxBjGmMnAZIB+/fqZli3hrrviYHSCEYHrry+/bsECb5OFhLYJ6DR7mUSf\nPvqylGdUdo8TrUCNGnD11dXc119TLG769+9PUVERxcXFHDmiE2OMHDmy3DaBzy1EydqJMdw47QbU\nEZE6wBhgZnKtslgyEysCccSZGGPYsGF07dqV0aNH0717dyZNmsSkSToxxogRIwAOA+vI8okxHJx2\nA04FvgReNsasSq5VFktmIibFk21F5BtgY4yHaQns8MGceJ23ozHGc3emiOwD1lTbqsQRa7t7bpfA\ndfJdjOdLFLG0S7TXih+/Hzfx+C0l7DqBtGkTP47rqV1SXgT8QEQWG2MS3oUar/Mm6++JlkTbadsl\n/sTD9nRuD0j/37kNB1ksFksWY0XAYrFYsphsEYHJGXbeZP090ZJoO227xJ942J7O7QFp/jvPij4B\ni8VisYQnWzwBi8VisYTBioDFYrFkMRkhAiKyQURWiMgyEVkcWNdcRN4VkaLAsplr+18E6tSvERHP\nFUhE5B8isl1EVrrWRX0eETkjYO86EfmriPfZY0VkeOB4UdXZ9wu/2jqWNghjU1LbJGBDyrWLR7vb\ni8j7IlIoIqtE5KeB9Q+JSGng71kmIiOitH1vOrZHnNvkryl5nRhj0v4FbABahqz7HXBv4P29wG8D\n77sBy4G6QCdgPVDT43mGAH2BlbGcB1gEDAQEmAVc7PH8NQPHyUfrDS0HuqVjW1e3DVKxTVKxXaKw\nux3QN/C+EcE5HB4C/jvM9l5tPwBck27tEec2mQVsS7XrJCM8gQiMAp4LvH8OuMy1fpox5rAxphgt\n1zDAywGNMR8AO2M5j+hcAY2NMQuN/jenuPapigHAOmPMVybKOvtxJpFtEEqqtgkkt108YYzZYoxZ\nGni/Dy3TUVnZbq+2fweMCLNvSrcHxLVNpgANIuyftHbJFBEwwHsiskS0vDBAGxMsxLYVaBN4nwO4\nS1R6rlUfgWjPkxN4X53z+217dfCjrWNpg1BSoU0g9dolakQkDzgd+DSw6nbRKU//4QpReLX9KHBZ\nOrcH+N4mJajnmlLXSaaIwNnGmD7AxcAEERni/jKglnHPhU3UeZJMSrR1CpLW7SIiDYFXgTuNMXvR\naU7zgT7AFuAPUR5yAvAhadoeEJc2Afgk1a6TjBABY0xpYLkdmI6GCLYF3CYCy8BU5pQC7mmrwtaq\nj4Joz1MaeF+d8/tte9T41NaxtEEoSW8TSMl28YyI1EZvdi8YY14DMMZsM8aUGWOOo9VtnZCpV9vr\nAaXp2B4B++LRJrloXD+lrpO0FwEROUlEGjnvgYuAlWj9+R8FNvsR4EzuOxMYIyJ1RaQTOsH7ohhM\niOo8Abdvr4gMDPToj3PtUxWfAZ1FpJMkoc6+X20dYxuEktQ2gZRtF6+2C/AM8KUx5o+u9e55ri8P\n/D1ebT8J+C9gRrq1B8StTQRtk3cCx0qddqluj3KqvFD3bHngtQq4P7C+BTAHKALeA5q79rkfVeQ1\nRNGrDkxF3cCjaBzuhuqcB+gX+OevByYSGLnt0YYRaLbCeudvTce2jqUNUqlNUrldPNp+Nhp++AJY\nFniNAJ4HVgTWzwTaRWH7RrQEctq1RxzbZD3agZty14ktG2GxWCxZTNqHgywWi8VSfawIWCwWSxZj\nRcBisViymFrJNqAqWrZsafLy8pJtRlxZsmTJDhPFHKnZ0CYQXbvYNglPNrSLbZPweG2XKkVARNqj\nvdpt0B7zycaYv7i+vwv4PdDKGLNDRC4EHkPruBwB7jbGzA1sOw+ty3EwsPtFRvNlI5KXl8fixYur\nMjOtEZGoJr3OhjaB6NrFtkl4sqFdRGSXiKxBR+M+bYx5rLLts6FNwPu14sUTOAbcZYxZGsiFXiIi\n7xpjCgMCcRGwybX9DuBSY8zXItIDeJvyQ5qvM8Zk/n/AYokPjaO54WU6ZWVlAB3QYmslwGciMtMY\nU5hUw9KIKkXA6MCELYH3+0TEKaZUCPwJuAfXQAVjzOeu3VcB9UWkrjHmcHWN3LB7A6u2r6ru7gxu\nP5hm9ZtVvaElKRw5AvPnw/HjMGQI1K8fv3Nt2AAtW0LDhvE7R7zIxhvetm1gDLRtG/77RYsWARw2\nxnwFICJOAcFK26S0FJYt89fWRCKi7QKwfDmUlECNGvDLX0Lr1tEdK6o+AXcxJREZhQ4LX15JKesr\ngaUhAvCciBxFh2Q/YsIMVAgUVroZoEOHDry+5nXumH1HNKaW47b+t/H4iMervb8lvkybBj8KjJf8\ny1/gjur/q6tk8GC46Sb41a/id454Ud0bXnEx7NgRf/viwYBAYYZFgTH933yjf0+7dnDFFVBaWgoa\ndnYoAc4MPU7oPWX+fLjuurianhSeeCIoDl7xLALuYkpoiOg+NBQUafvuwG9DtrnOGFMaCCu9CvwQ\n7W8ohzFmMoFJlvv162dGdx/NwNyBXk0tx4CnB7D/6P5q7ZsIDh06BNBVRJaj/49/G2MeTK5VieXA\ngeD7/XH8VxkDW7fqKx3xesML5eGH4Z//jJNRCWJAmGLvH37off/Qe8qwYUFhSTectvjJT+D664Mi\nX6OGfo4WTyIQWkxJRHqikxw4XkAusFREBhhjtopILlocaZwxZr1zHBMssrVPRF5EiydVEIFQ2jRs\nQ5uGbaraLCwdmnSo1n6Jom7dugBrjDF9A+38kYjMMsYsTLJpSSGeA9gPHtTjf/dd/M6RCoQ+9d55\nJ1x1VZKNqia33aYhvLvvhnPPhUsu0fV//7t6dZ9+mgOahOLgqZhaixb6SkcaNYJ9+1QM+vfXVyx4\nyQ6qUEzJGLMCaO3aZgPQL5Ad1BR4E50p52PXNrWApoFtagOXoHUy4k4ql8YIiOjxwMfagVfqGhxn\n4vmvcm7+6SoCOTnebnihT729e0Pv3gkx0XcGD4aZM2HcOI2D5+TA8OFw4436fX+9A9YLFFgrRQsI\nXps0gxOIXxNtehksdhYatvmehJlbMwy3AQXAL13bt0anSHtbRJyiTKVoOda4IggmDe6pIrIMLSH7\nrjHm06q2z1SsCETGfcNLVsXURNOsmfYXRbrh1apVCzQ78W10BrCXjTHVzyJJI/wSAS/ZQR+h81hW\ntk2e6/0jwCMRNj0jGuP8oJJO65TCGNMn4EVNF5EexpiV7u9DXfxMJREi4O6DSCdCbng1gX9kyw2v\nCvYYY/ol24hEk0hPIO1J5XCQG2PMbuB9YHiY7yYbY/oZY/q1auV5cGTa4de/av58OHas/Drn5p+u\nnkCAPcaYU40xpxhjHk22MZbkYUXAI6keDvrmm29An+oQkfrAhcDqZNqUTPwQga++gvPOg7feKr8+\n3cNBFoubhIWD0h0RSWlPYMuWLQBdAn0lNdCY5hvJtSp5+PGvctJMd+8uv96KgCWTsCLgEam8OyPp\n9OrVC6AwG2Oa4fBDBJww0JEj5ddbEbBkEjYcFAWpHA7KNsaPH0/r1q3p0aNH2O+Li+fRpEkT+vTp\nA9BNRH4Z7TkcETgcUqjEioAlk7Ai4JFUDwdlG9dffz2zZ8+udJtzzjmHZVrYpdAY83C05zh6VJeR\nRODYsYpegsWSblgR8Eiqh4OyjSFDhtC8efOI3yciHBT63mJJR6wIRIENB6UPxsAnn3zi9JV0DtSg\nCouI3Cwii0VkcSDLCqg6HBT63mJJR6wIeMSGg9KLNm36smnTJr744gvQEdT/ibRtpLETVgQs2YAV\nAY+k+jgBS3nq1GlMw2Cx/z1AbRFpGc0xbDjIkg1YEfBIupSNsCj79m11e24N0Gv022iOYT0BSzZg\nxwlEgQ0HpQ5jx45l3rx57Nixg9zcXH71q1/xwQeBdB5uYc2af9Ojx1NOnZwO6DzUUf0DrQhYsgEr\nAh6x4aDUYurUqRXWHT0KzurTT7+N9967DQARWW2M+STac1QWDqpfX+cVsCJgSXdsOMgjNhyUXviZ\nIhrOE3DmX7UiYEl3rAhEgQ0HpQ9WBCwWb1gR8IgNB6UX8R4sZkXAkilYEfCIHSeQXvjxr6qsbIQz\nnCBdJ5axWBysCHjElo1IL+IdDmrcGOrWtZ6AJf2xIhAFNhyUPsQ7HHTSSfqyImBJd6wIeMSGg9KL\neHkCR47oeisClkzBioBHbDgovYiXCDh9AOkqAq+88grdu3cHOENE7AREFisC0WDDQelDvMJBzk0/\nXUWgR48evPbaawD7k22LJTWwI4Y9YsNB6UW8PAHnpt+gQXqKQNeuXZNtgiXFsJ6AR1J9nMDmzZsB\nThWRQhFZJSI/TbZNySTeIpCunoDFEor1BDyS6mUjAoXSSowx3USkEbBERN41xhQm2bSkkAgRaNAg\nNUVg6NChbN26tcL6Rx99lFGjRnk+jojcDNwM0KFDB9/ss6QWCRMBEWkPTAHaAAaYbIz5i+v7u4Df\nA62MMTtE5ELgMaAOcAS42xgzN7DtGcA/gfrAW8BPo60QWR1SORzUrl07gAMAxph9IvIlkANYEagm\n6don8N577/lyHGPMZGAyQL9+/VL34rfERCLDQceAu4wx3YCBwAQR6aZGSHvgImCTa/sdwKXGmJ7A\nj4DnXd89BdwEdA68hsf8F1RBqoeD3IhIHnA68GlyLUkefo8Ydo7nFoFGjWDTJl02agQPPKDfLVsG\np58Oe/dWfY5Dh2DfvthttViqS8JEwBizxRizNPB+H+A8qQL8CbgHgndZY8znxpivAx9XAfVFpK6I\ntAMaG2MWBp7+pwCX+fNnRCbVw0EOItIQeBW40xhT4TYUaT7dTMNPT8CY4PsdO3TZogXceivcdRfc\nfLPWEnr7bf1uwQIVgqKiqs9xzz0wbFjstnpl+vTp5ObmApwEvCkibyfu7JZUJCl9Au4nVREZBZQa\nY5ZXcqO9ElhqjDksIjlAieu7EoJiEldSORwUQFABeMEY81q4DbLBxRfxVwRAQ0K1a8O2bfq5dWvI\nzYXf/14/790LM2fq+507yy8rY8MGWLs2dlu9cvnll3P55ZcjIkuNMXacgCXx2UHuJ1U0RHQf8MtK\ntu8O/Bb4cbRG+fnUm+rhoIBAdQS+NMb8McnmxJ3x48fTunVrevToUeE7ETh+3HDHHXdQUFAA0E1E\n+kZ7DrcIOJ3DW7dCs2ZQp075bTt2hO3bdaKZbwOTWHoRgX37dLuysmits1j8IaEiICK1Kf+kegrQ\nCVguIhuAXGCpiLQNbJ8LTAfGGWPWBw5TGtjOITewrgLGmMnGmH7GmH6tnLKP1STVxwl8/PHHAC2A\n74nIssBrRJLNihvXX389s2fPDvudCGzePIuioiKKNCazEe1HiopwIrBtG7RtW3Hbjh11uWlTUAS+\n9ZF0roUAABN9SURBVDCj8f796rXs3h2tdRaLPyRMBERjPc/gelI1xqwwxrQ2xuQZY/LQ0E5fY8xW\nEWkKvAnca4z52DmOMWYLsFdEBgaOOQ6Y4c+fUYn9KV424uyzzwZYYozpZYzpE3i9lWy74sWQIUNo\n3rx52O9EYOPGGYwbN87py/kOaBroT/JMaDgIVATatKm4rSMCGzdGFw7aHxi36/Q1WCyJJpGewFnA\nD/H+pHobUAD80rV9YCoPfgI8DawD1gOzYrDdM6kcDrIEEYEDB0pp3769e3XUfUeRPIGqRCBaT8Dr\nthZLPEhYx7Ax5iOo/HE64A047x8BHomw3WKgYjA4jqR6OMgSJNqLOtKgqGhEICcHata0noAl/bBl\nIzyS6uEgSxARqF8/xyml4RB131FoOOjgQc0CCicCtWppttCGDcGn+q1bYerU8sdZuhS+/BJWrdI0\n0kiewJ49MCsh/q0l27EiEAU2HJQeiEBu7kimTJnieG8nAXsC/UmeCfUEnPTQcCIA0KFDeU/gnXfg\n2mth0qTgNj/4AfzkJ9Cjhw4oCx1/4HD11TBiBGTwUA5LimBFwCM2HJRajB07lkGDBrFmzRpyc3N5\n5pln+OCDScAkRODkk0eQn5/vpIh2RPuRoiKSCITLDgJo1w7WrIHjx8uvf+QRHWl88KB+v3x5xX1D\nPYHPP9elveQs8cYWkPNIqo8TyDamTp1aYd3Roxp+0YtaeOKJJwAQkcJAP1JUHD0KNWroTf3IkWB5\nh0iewMknB5/cnbpCjRqpeDzxBAwdqsfatavivqGewMGDBGyP1mqLJTqsJ+CRdCkbYfF3xPBJJ+n7\nw4c1xg+RRaCdKwFVHRAtCXH22fCvf8GKFZHPFeoJOCJQ2SCydeugpCTy9xaLF6wIRIENB6UH8RIB\nd8mIcLhF4OSTddmxI1x5pQrAf/6jpSfCEeoJOCGlykSgc2conwVrqQ7ZPuWmFQGP2HBQ+hAPEfjZ\nz+Cjj6BpU6hbN/z2bhFwBpedfDJcFihv+J//QM+emk7qpnVrzRi64w54+uny34X2L1j8J9un3LQi\n4BEbDkof/BSBDh3gtNPUC3jnncihIAg+/QPcf7/e7MeMgbw8uOIKrTl0xRVw3XXQz/W8ecklerN/\n/HGtSOomkidgnVL/6Nq1K126dEm2GUnDdgxHgQ0HpQd+ikDTpjB3rubtP/po5eEXtydwzjnl4/Wv\nvlp+28WLoX9/fX/vvfDMM5pF9MAD5SexieQJpOJkNpb0xIqAR2w4KH3wUwRqBa7sJk3gd7+rfHsn\nVFSvXnC/SLRoEXzfsGH55X5XUCKSJ2BHGEdH6JSb27apMF96qZ1y04qAR0SE4zZAmxbEQwS8nrdd\nO00rrYpwItCokS7dIhDpkrMiEB2hU27m5sLw4RDF/R/IzPk4bJ+AR2zZiPQhWSIA2g/QsmXV2zVq\nFDx2gwa6rI4nUKMGrF8fnArTYokWKwJRYMNB6UEyReC3v4X/+7+qtxOB5s1VAGrW1HWOCLjnHI4k\nAs64guPHdUzCXXdFZ6clSLZPuWnDQR6xZSPSh2SKwFlned82dDqEWMJBc+d6P6+lPNk+5ab1BDxi\nw0Hpg18icPRo9CIQDS1aBG/8EPQE3CWoI3kCTukJp//B8Saq4u677+a0004DnXJzemDyJksWY0Ug\nCmw4KD1IpicQDXl55dNKHRFwJbFE9ASc4nbO9146owEuvPBCVq5cCVAIrAV+4dlgS0Ziw0EeseGg\n9CFdRODxx8uPCXC8ArcIOJ7AsWP6tO/8YEM9BK+ewEUXXeT+uBC4KgqTLRmI9QQ8YsNB6UO6iECz\nZuVHIDuegFOjCIJP+rVrw+23B9eHikA17RxPgqZmtaQezm/EikAUpEE4KE9EtovIymQbkig2bYrf\nsY8di1zwLR44dYpCPQEn/TNQGfvEejduT2Do0KH06NGjwmvGjBnuXdoCx4AXItkjIjeLyGIRWfyN\nnd0mY7HhII+kSThoB3AtMCXZhiSCDz+EIUPguedg3Ljgej88AWPi7wmEUqOGCoHbEygrC5aVdhMq\nAu4+gdCBUaH885//BGgKXGcquagzcWCUpSLWE/BImpSN2A94mN48M1gZ8HcWLCi/3g8RcMIwiRQB\n0HEDi13T32zfDu6pkvfs0aVbKMB7n8Ds2bP5nda/WGeMORCLrZbMwIqAR2wV0dRj1arZQBdefLGA\nxx577MR6RwTmzZtHkyZNQNMhl4nIL70e28m+SbQINGtW/vPo0TofsUPTpioS06aV386rnbfddhv7\ndDTaqYE2mVTVPpbMxopAFKRBOKhKMiXOW1ZWxksvTQBmcc01hUydOpUtWwqB8p7AOeecA1BojOlj\njHnY6/GTJQLlw/bhWRxmoky3J/DWW/DVV+H3XbduHZvVtXDa5Jbq2GnJHKwIeCRNwkFVYoyZbIzp\nZ4zp16pVq2SbU20WLVpEq1YFQD41a9ZhzJgxLF+ud1A/wkHJEgEdx1U5hw5VXOfuE/j+96FbN/9s\nsmQ2VgQ8YsNBqUVpaSnNmgWL++fm5rJ7dylQXgQ++eQT0HDQLBHp7vX4TkZOokUAYNCgyr8PJwKh\nfQKHD/tnjyWzsSIQBWkQDuoELAC6iEiJiNyQbIOSgSMCffv2ZZPmkBYCjwP/Cb99xRBZsjwBgLff\n1oFkkQgnAo6dlc1JbLGEI2EiICLtReR9ESkUkVUi8tOQ7+8SESMiLQOfWwS23y8iE0O2nSciawId\nW8tEJMLU3/6RJuGgYmNMO2NMbWNMrjHmmWQbFC9ycnLYtSuYNlNSUkLTpjp5ryMCjRs3pmFgBJYx\n5i2gtnN9uQkXIkumCDRqBGefHfn7A2FyehxPIJxAWCyVkUhP4BhwlzGmGzAQmCAi3dQIaQ9cBLiH\n/hwCHgD+O8Lxrgt0bPUxxmyvvuneSJNxAllD//792b69CCimrOwI06ZNo1evkUBQBLZu3XrifyYi\nA9Dr9Fsvx3dEIJGDxdxUlvLpLjXt4PQJWBGwREvCRMAYs8UYszTwfh/wJZAT+PpPwD0QfNQ2xnxn\njPkIFYOkY8tGpBa1atXimmsmAsN46aWujB49mpNP7g5MYt++SRgD//73v+mh+ZXdgL8CYyobHOUm\nmZ4AVF4Qbu/eyNuHG1hmsVRGUvoERCQPOB34VERGAaXGmOVRnvO5QCjoAUlQr20ahIOyiu7dRwBr\nufba9dx///2BtbfQpMktGKM58atWrQJNhxxojPnE67GTLQKVeQLhRMDpC7CegCVaEi4CItIQeBW4\nEw0R3Qd4HsQT4DpjTHfgnMDrhxHO5VtOvA0HpS6hF3E6p4g6VOYJhAsHOSKwaFF87LFkLgkVARGp\njQrAC8aY14BT0IyW5SKyAcgFlopI28qOY4wpDSz3AS8CAyJs51tOvA0HpQ+ZIALRegKOvdddFx97\nLJlLwgrIBUI2zwBfGmP+CGCMWQG0dm2zAehnjNkR9iC6TS2gqTFmR0BULgEqr5jlEzYclJqE3vD9\nEAGnzn+yOoaj7RM4dqz83AQWi1cSWUX0LDRss0JElgXW3RdI3QtLQBQaA3VE5DI0g2gj8HZAAGqi\nAvD3GGz3hA0HpQ9+iIDTwdqgQez2VIfKPIH16yuuO3YsfJjIYqmKhIlAINOn0tMZY/Iq++ziDK+G\n+UWajBPIKiJdvH6IgJOLX79+bMepLtEO+iorKz9BvcXiFTti2CO2bET64KcIJMsTiLbsw/vv65zF\nbm6+WZfz52ubLF3qi2mWDMOKQBTYcFBqEenfkQki0KUL3H8//Nd/Vf8Yfw8ESc87T5eXXhqzWZYM\nxIqAR2w4KHWJR4poskVABB55BP7xD2jRQtfNCswG/MYb3o/jTI4T+t5icbDTS3rEhoNSj0g3tUzo\nGHazYYN2/DZtCiUlkJNT+fbdukGhTq1QbvCYLS5nCYf1BKLAhoNSi0g3tUzwBNw0bKgCABUF4PPP\ng2EfB/d8y+4yEra8tCUc1hPwiA0HpR7xFoEaNZI3TsArffqUt7FmTZ2s3sFdcdSWlLCEw3oCHrHj\nBFIPRwTiMVjswAH1AtIhChioln0Cd1rrhg3B93YwmSUcVgQ8YstGpB7x7BNwRCAdqFu3/Hu33UOG\nJN4eS3phRSAKbDgotYh3OCiVReCpp4JZQm3awD33wGWXwUcfVR7CeuCBB+jVqxfolJvviMjJibDX\nkrpYEfCIDQelHo4IxCNF9ODB1BaBW27RCeVB/97f/hamT4fTT4etWyPvd/fdd/PFF1+ATrn5BtFX\n8LVkGFYEPGLDQalHIvoE0pGRIyN/17hxY/fHk8C6txZ/yPjsILDhoFTDEYHQvgGR2AdGHTiQvLpB\nsdKhA2zfDq0jzLwdmICnF1APOD/ScUTkZuBmPWYH3+20ZBaZ7wmkRziosYisEZF1InJvso2JN2vX\nzga6MHVqAY899tiJ9Y4nYIzhjjvuAOghIl+ISF+vx05nT2Do0KGcf34PIPQ1A4BHH30U4AvgBeC2\nSMfxcz4OS+aT8Z5Aqo8TKNPH4g7ofLolwGciMtMYU5hUw+JEWVkZ77wzAXiXyy7LZerU/lx22Uig\n2wkRmDVrFkVFRQAr0ZvdU8CZXo5/4EDVI3NTlffe0+k1wsV69+2DRo1OfHwBeAt4MDGWWTKZzBcB\nEb7a9RV/WvCnZJtSgb7t+lJnSx2Aw8aYrwBEZBowCu0AjMjOnfDcc/G30W82bFjE3r0FQD7PPw8j\nRozhySdnAN1YtAh69oQZM2Ywbtw4Zs+ejTFmoYg0FZF2xpgtVR0/nT0Bh9xcLTPhZsuWIho16ux8\nHAWsTrBZlgwl40Ugr0keAD9/5+fJNSQMPx/4cwbtGwTgHg5UQpin3tA47/bt8PPU+5M8UAq0P/Hp\nrbdygU9PfM7Lg9LSUtq3b+/eqQTIAcqJQLjYd6pnB3lh6VKYORNuvFGrko4bB7/4xb2sWbMG1GO8\nCLgluVZaksUvfqGVauvV8+d4GS8CD5//MP89+L+TbUZY6tSsw5sz3vS0rTFmMjAZoF+/fqZzZ9i9\nO57WxYcZM+C997RuzqFDMG2aTrL+4IMa7mjcGEaN8nas0DYBKCqKPcMo2bRqBddfD1ddpZ3cdeqA\nTvENIlJojMna4tKFhcmbPzpVuO8+fflFxjeniNCkXpNkmxGRHA1g13GtykUflyulZk1okrp/VkQ6\nd87hhRc2U7eujpLdubOEgoIc3EksOTk5bN682b2bpzaB8vV30pl0/f/Gm/KZshY/yPjsoFSnf//+\nAPVEpJOI1AHGADOTa1X86N+/P0VFRRQXF3PkyBGmTZvGyJAE+ZEjRzJlyhQARGQgsMdLf4DFYoke\nKwJJppb6tpuAt4EvgZeNMauSalQcqVWrFhMnTmTYsGF07dqV0aNH0717dyZNmsSkSZMAGDFiBPn5\n+aD5kX8HfpJEky2WjEZSPYdeRL4BNvp82JbADp+PGcuxOxpjPCd0x9gm8fzb/T6H53axbRKeGNrF\ntklF0qlNwGO7pLwIxAMRWWyM6Zdux/aDRNiX6m0Qim2Titg2qUimtokNB1ksFksWY0XAYrFYsphs\nFYHJaXpsP0iEfaneBqHYNqmIbZOKZGSbZGWfgMVisViUbPUELBaLxUKGiICItBeR90WkUERWichP\nA+sfEpFSEVkWeI1w7fOLQOnmNSIyzLX+DBFZEfjur6JsCKxbJiKLA9s1F5F3RaQosGxWnWPHoS18\nsTXkmP8Qke0istK1LiX//kj43S7p3ib2OglPVraLMSbtX0A7oG/gfSNgLVpo6yHgv8Ns3w1YDtQF\nOgHrgZqB7xYBAwEBZgEXAxuAliHH+B1wb+D9vcBvq3PsOLSFL7aG7D8E6AusTPW/P1Htku5tYq8T\n2y4nbElU4ybyhc7CcSGRReAXwC9cn98GBqFistq1fizwtwgXxhqgXeB9O2BNdY6doIs4KlsjHDcv\n5CJOyb8/ke2Szm1irxPbLs4rI8JBbkQkDzidYH3i20Vnp/qHy+XKAdwVypxSxTmB96HrDfCeiCwR\nLV8M0MYE69lsBdpU89h+44etXkjVvz8SiWiXdGoTe52EJ+vaJaOqiIpIQ7Tm7p3GmL0i8hTwv+g/\n9v+3d/8qDQRBHMe/U4hFsFEQUir4BiIWeQI7OysD+hgprbWytRKxU7AS1Bew8i8iYpkiKW1Fx2I3\nmMgF4nkel+zvAwfhuJvbnSwMd5vc7gC7wFaO0A13b5vZPHBhZgMLeri7m1lVfmZVelsr1v9hSs3L\nGORE4yRbcnmZmDsBM5siFIAjdz8BcPeOu3+4+yfhRWQr8fDBlU2+X1Xcjp8H9rt7O8brAqcxTsfM\n6vHadaCbJ/Zf+pyloLaOopL9H6akvIxNTjROsqWYl4koAnGW/AB4cve9vv31vsPWCWvWQnhV84aZ\nTZvZArAEXMfbszczW40xN4FzM5uJ8WqEVZ0eYoxmjNektxr472L3zikqD7Ui2jri5SrX/2FKzMtY\n5ETjJFuyeSlrwuU/N6BBeORzB9zEbQ04BO7j/jPiREw8p0WYeX+mb5YdWCZ88a/APrBImK2/BR6B\nVjxuDrgCXoBLYDZHbCs4D4W19UfcY8LSju+EZ5HbVex/mXkZ55xonCgv/Zv+MSwikrCJeBwkIiL5\nqAiIiCRMRUBEJGEqAiIiCVMREBFJmIqAiEjCVARERBKmIiAikrAvE6gOvKAq8QkAAAAASUVORK5C\nYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f820646f748>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "saves = 10000\n",
    "rows,cols=2,5\n",
    "# data = machine\n",
    "data = people\n",
    "fig,plts = plt.subplots(rows,cols)\n",
    "for i,s in enumerate(data[882:]):\n",
    "    if i >= rows:\n",
    "        break\n",
    "    index, start, track, des, flag = get_track(s)\n",
    "    track = np.vstack([track[0], track, track[-1]])#添加固定偏差来保证数据有效\n",
    "    delta = track[1:] - track[:-1]\n",
    "    delta[:,2]+=1#固定偏差用来消除0\n",
    "    \n",
    "    angular = np.arccos(delta[:,0]/np.sqrt(np.square(delta[:,0])+np.square(delta[:,1])))\n",
    "    angular = np.nan_to_num(angular)\n",
    "    angularDeltaTime = np.hstack((angular, delta[:,2]))#x1角度,x2为时间差\n",
    "    timeAngular = np.hstack((track[1:,2],angular))#x1为时间，x2为角度\n",
    "    \n",
    "    angularVelocity = (angular[1:] - angular[:-1])/delta[1:,2]#角速度\n",
    "    angularVelocityDeltaTime = np.hstack((np.nan_to_num(angularVelocity), delta[1:,2]))#x1角度,x2为时间差\n",
    "    timeAngularVelocity = np.hstack((track[2:,2],np.nan_to_num(angularVelocity)))#x1为时间，x2为角度\n",
    "    \n",
    "    speed = np.sqrt(np.square(delta[:,0])+np.square(delta[:,1]))/delta[:,2]\n",
    "    timeSpeed = np.hstack((track[1:,2],np.nan_to_num(speed)))#x1为时间，x2为速度\n",
    "    speedDeltaTime = np.hstack((np.nan_to_num(speed), delta[:,2]))#x1为速度，x2为时间差\n",
    "    \n",
    "    acceleration = (speed[1:,0] - speed[:-1,0])/delta[1:,2]#加速度\n",
    "    accelerationDeltaTime = np.hstack((np.nan_to_num(acceleration), delta[1:,2]))\n",
    "    timeAcceleration = np.hstack((track[2:,2], np.nan_to_num(acceleration)))\n",
    "    \n",
    "    trackGrad = (track[1:,1] - track[:-1,1])/(track[1:,0] - track[:-1,0])\n",
    "    speedY = delta[:,1]/delta[:,2]\n",
    "#     trackHessen = (trackGrad[1:] - trackGrad[:-1])/(trackGrad[1:,0] - trackGrad[:-1,0])\n",
    "    \n",
    "    \n",
    "\n",
    "    canvas = plts[i][0]#轨迹\n",
    "    canvas.plot(track[:,2],track[:,1],'-g')\n",
    "#     canvas.plot(*start, '.r')\n",
    "#     stabs = find_stab(track, stride=100, height=25, max_t=100)\n",
    "    stabs = find_stab_2(track, height=0)\n",
    "    deltax = track[:,0].max() - track[:,0].min()\n",
    "    deltay = track[:,1].max() - track[:,1].min()\n",
    "#     stabs_t = np.abs((stabs[1:-1,1] - (stabs[2:,1]+stabs[:-2,1])/2)/deltay*deltax/(stabs[2:,0]-stabs[:-2,0]+1))\n",
    "#     t = stabs_t.max()\n",
    "#     print(stabs_t)\n",
    "#     print(index,stabs_t.max(), np.count_nonzero(stabs_t>1.5))\n",
    "#     canvas.plot(stabs[:,0], stabs[:,1], '.r') if stabs is not None else False\n",
    "#     canvas.plot(*des,'.r')\n",
    "    canvas = plts[i][1]#角度\n",
    "    canvas.plot(track[1:,0],angular,'-b')\n",
    "    canvas = plts[i][2]#移动速度\n",
    "    canvas.plot(track[1:,0],speed,'-b')\n",
    "    canvas = plts[i][3]#角速度\n",
    "    canvas.plot(track[2:,0],angularVelocity,'-b')\n",
    "#     canvas = plts[i][4]#加速度\n",
    "#     canvas.plot(track[2:,0],acceleration,'-b')\n",
    "    canvas = plts[i][4]#轨迹Hessen\n",
    "    canvas.plot(track[1:,0],trackGrad,'-b')\n",
    "#     print(track, np.nan_to_num(trackGrad))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "metadata": {
    "collapsed": true
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/miller/anaconda3/lib/python3.6/site-packages/ipykernel/__main__.py:14: RuntimeWarning: invalid value encountered in true_divide\n",
      "/home/miller/anaconda3/lib/python3.6/site-packages/ipykernel/__main__.py:15: RuntimeWarning: invalid value encountered in true_divide\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Done\n"
     ]
    }
   ],
   "source": [
    "data = people\n",
    "saves = len(data)\n",
    "paramNums = 21\n",
    "parameters = np.empty((saves, paramNums))\n",
    "for i,s in enumerate(data):\n",
    "    if i >= saves:\n",
    "        break\n",
    "    index, start, track, des, flag = get_track(s)\n",
    "#     index, start, track, des = get_test_track(s)\n",
    "    track = np.vstack([track[0], track, track[-1]])#添加固定偏差来保证数据有效\n",
    "    delta = track[1:] - track[:-1]\n",
    "    delta[:,2]+=1#固定偏差用来消除0\n",
    "    \n",
    "    angular = np.arccos(delta[:,0]/np.sqrt(np.square(delta[:,0])+np.square(delta[:,1])))\n",
    "    angular = np.nan_to_num(angular)\n",
    "    angularDeltaTime = np.hstack((np.nan_to_num(angular), delta[:,2]))#x1角度,x2为时间差\n",
    "    timeAngular = np.hstack((track[1:,2],np.nan_to_num(angular)))#x1为时间，x2为角度\n",
    "    \n",
    "    angularVelocity = (angular[1:] - angular[:-1])/delta[1:,2]#角速度\n",
    "    angularVelocityDeltaTime = np.hstack((np.nan_to_num(angularVelocity), delta[1:,2]))#x1角度,x2为时间差\n",
    "    timeAngularVelocity = np.hstack((track[2:,2],np.nan_to_num(angularVelocity)))#x1为时间，x2为角度\n",
    "    \n",
    "    speed = np.sqrt(np.square(delta[:,0])+np.square(delta[:,1]))/delta[:,2]\n",
    "    timeSpeed = np.hstack((track[1:,2],np.nan_to_num(speed)))#x1为时间，x2为速度\n",
    "    speedDeltaTime = np.hstack((np.nan_to_num(speed), delta[:,2]))#x1为速度，x2为时间差\n",
    "    \n",
    "    acceleration = (speed[1:,0] - speed[:-1,0])/delta[1:,2]#加速度\n",
    "    accelerationDeltaTime = np.hstack((np.nan_to_num(acceleration), delta[1:,2]))\n",
    "    timeAcceleration = np.hstack((track[2:,2], np.nan_to_num(acceleration)))\n",
    "    \n",
    "    stabs = find_stab_2(track,height=0)\n",
    "    if stabs is not None and stabs.shape[0] >= 3:\n",
    "        stabMax = np.abs((stabs[1:-1,1] - (stabs[2:,1]+stabs[:-2,1])/2)/(stabs[2:,0]-stabs[:-2,0]+1)).max()\n",
    "    else:\n",
    "        stabMax = 0\n",
    "#     indexs[i] = index\n",
    "    parameters[i] = [speed[:,0].max(),\n",
    "                     speed[:,0].mean(), \n",
    "                     speed[:,0].std(), \n",
    "                     get_distr(timeSpeed,part=4)[0][0],\n",
    "                     get_distr(speedDeltaTime,part=4)[0][0],\n",
    "                     \n",
    "                     acceleration[:,0].max(),\n",
    "                     acceleration[:,0].mean(),\n",
    "                     acceleration[:,0].std(),\n",
    "                     get_distr(timeAcceleration,part=4)[0][0],\n",
    "                     get_distr(accelerationDeltaTime,part=4)[0][0],\n",
    "                     \n",
    "                     angular[:,0].max(),\n",
    "                     angular[:,0].mean(),\n",
    "                     angular[:,0].std(),\n",
    "                     get_distr(timeAngular,part=4)[0][0],\n",
    "                     get_distr(angularDeltaTime,part=16)[0].max(),\n",
    "                     \n",
    "                     angularVelocity[:,0].max(),\n",
    "                     np.abs(angularVelocity[:,0]).mean(),\n",
    "                     angularVelocity[:,0].std(),\n",
    "                     get_distr(np.abs(timeAngularVelocity),part=4)[0][0],\n",
    "                     get_distr(np.abs(angularVelocityDeltaTime),part=8)[0].max(),\n",
    "                     \n",
    "                     stabMax,\n",
    "                    ]\n",
    "    if np.any(np.isnan(parameters[i])) or np.any(parameters[i]>10000):\n",
    "        print(index)\n",
    "        print(parameters[i])\n",
    "        break\n",
    "    \n",
    "root = 'processed'\n",
    "\n",
    "with open('%s/human_%d'%(root, paramNums), 'w') as f:\n",
    "    parameters.tofile(f)\n",
    "    print(\"Done\")\n",
    "# with open('test_index', 'w') as f:\n",
    "#     indexs.tofile(f)\n",
    "#     print(\"Done\")\n",
    "#     print('speed-max:%f'%speed[:,0].max()\n",
    "#           ,\"mean:%f\"%speed[:,0].mean()\n",
    "#           ,\"std:%f\"%speed[:,0].std()\n",
    "#           ,\"%%25time speed:%f\"%(get_distr(timeSpeed,part=4)[0][0])\n",
    "#           ,\"%%25speed time:%f\"%(get_distr(speedDeltaTime,part=4)[0][0])\n",
    "#          )    \n",
    "\n",
    "#     print('acceleration-max:%f'%acceleration[:,0].max()\n",
    "#           ,\"mean:%f\"%acceleration[:,0].mean()\n",
    "#           ,\"std:%f\"%acceleration[:,0].std()\n",
    "#           ,\"%%25time acceleration:%f\"%(get_distr(timeAcceleration,part=4)[0][0])\n",
    "#           ,\"%%25acceleration time:%f\"%(get_distr(accelerationDeltaTime,part=4)[0][0])\n",
    "#          )\n",
    "\n",
    "#     print('angular-max:%f'%angular[:,0].max()\n",
    "#       ,\"mean:%f\"%angular[:,0].mean()\n",
    "#       ,\"std:%f\"%angular[:,0].std()\n",
    "#       ,\"%%25time angular:%f\"%(get_distr(timeAngular,part=4)[0][0])\n",
    "#       ,\"longest angular time:%f\"%(get_distr(angularDeltaTime,part=16)[0].max())\n",
    "#      )\n",
    "    \n",
    "#     print('angularVelocity-max:%f'%angularVelocity[:,0].max()\n",
    "#       ,\"mean:%f\"%np.abs(angularVelocity[:,0]).mean()\n",
    "#       ,\"std:%f\"%angularVelocity[:,0].std()\n",
    "#       ,\"time angularVelocity:%f\"%(get_distr(np.abs(timeAngularVelocity),part=4)[0][0])\n",
    "#       ,\"angularVelocity time:%f\"%(get_distr(np.abs(angularVelocityDeltaTime),part=8)[0][0])\n",
    "#      )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": true
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/miller/anaconda3/lib/python3.6/site-packages/ipykernel/__main__.py:18: RuntimeWarning: invalid value encountered in true_divide\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD3CAYAAAAT+Z8iAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAGuZJREFUeJzt3X+MG+d95/H3StQPq7sWiIY6t6gdO7H01WYBn2vVNmJb\n8iJxosgpasDAXYM0SIpLFV+gxkgbwEEiFYULJYbTwknVIu5BaeH40kOBc1ygUaPorhbkSmoA6RQb\nsGXtVz+cixCgbrbJWlrC+mHK2z9muJqluFxydkgOn/m8/lny4XD4fIezXz7zPM/MDM3MzCAiIuFa\n0u8KiIhIdynRi4gEToleRCRwSvQiIoFTohcRCVyp3xVoZnJyOjdTgcrlVUxNvdXvanSVYgxHEeJU\njPOrVEaGmpWrRb+AUmlpv6vQdYoxHEWIUzF2ToleRCRwSvQiIoFTohcRCZwSvYhI4JToRUQC19b0\nSjO7G3jS3ccTZR8HPufu74+fbwUeAWrATnffY2bXAd8B1gDTwKfcfTLbEEREpJUFW/Rm9hjwLWBl\nouzXgU8DQ/HzG4BHgXuBzcATZrYC+CzwirtvBJ4FdmQdQN9Uq5SOHYVqtd81EZmf9lOhva6bM8DD\n9Sdm9svAV4HPJ5a5Czjs7pfc/RxwGrgNuA/4QbzMXuCBLCrdd9Uq5c3jlLd8kPLmcf0TST5pP5XY\ngl037v5dM7sZwMyWAn8N/CFwIbHY9cC5xPNpYHVDeb1sQeXyqlydFFGpjMwteP01OHUSgNKpk1R+\ndhZuuTu7D6xW4fhxGBuD4eGerPOaGANUhBghEWe399M+KsJ3mWWMnV4CYQOwFniaqCvnfWb2DWA/\nkKzVCPAmcD5RXi9bUJ5Ob65URpicnJ5buOYmymvXUTp1ktradUytuQkal0krboXNrnvfgcUn+wXW\n2TTGwBQhRmiIs5v7aR8V4btMG+N8Pw4dJXp3PwKMAcSt/L9z98/HffRfMbOVwApgFHgVOAw8CBwB\ntgAHO655VqpVSn6Cmo0uPnEODzO170B260so+QlKiVZYyU9Q23Bn7tYpA6CL+6kMlkymV7r7G8Au\nokS+H9ju7heJWv5jZnYI+AzweBaf17Fu9FUOD0fJMuN/npqNUlu7Lnq8dl30D5rDdcqA6NJ+KoNl\nKI/3jM366pWlY0cpb/ng7POpvS+03aLty2FilkcfbaxTh8LhKEKcirHl+4p79cqBa9F2oxWmlp1I\nYeXyevSZU1+liBRYIVr0gFq0Mlh0opNkqDiJvpH+kSSvdKJT2PqQe4qZ6PWPJDnWbDqsBKJPuaeQ\niV7/SJJnAzd5QNrWr9xTyESvfyTJtXjywNTeF7I5M1pyo1+5pxizbhppFo7kXX3ygISlT7mnmIke\n9I8kIv3Rh9xTyK4bEZEiUaIXEQmcEr1Ip3QOhgwYJXqRTugcDBlASvQiHdA5GDKIlOhFOqBzMGQQ\nFXd6pUgaOgdDBpASvUindA6GDBh13YiIBE6JXqRdmlYpA0qJPm+UTPJJ0yplgCnR5ymxKpnkVqbT\nKvO0z0khFDvR5yyxao52fmU2rTJn+5wUQ6ETfd4Sq+Zo51hG14jP2z4nxVDo6ZX1xFo6dTIfiVVz\ntPMtg2mVudvnpBAKnehzmVg1Rzt4008+BUDt9jvysc9J8NpK9GZ2N/Cku4+b2e3AXwBXgEvAJ939\n38xsK/AIUAN2uvseM7sO+A6wBpgGPuXuk90IJDUlVumVuH++3pqf2neg3zWSgliwj97MHgO+BayM\ni/4c+Jy7jwPPA180sxuAR4F7gc3AE2a2Avgs8Iq7bwSeBXZkHoHIgOhq/7xm8kgL7QzGngEeTjz/\nmLu/HD8uAReBu4DD7n7J3c8Bp4HbgPuAH8TL7gUeyKTWIgOoa4PtiZk83Hmnkv0g6dEP9IJdN+7+\nXTO7OfH8XwHM7B7g94FNRK34c4m3TQOrgesT5fWyBZXLqyiVlrazaE9UKiP9rkLXKcZeVGAEfnQM\njh+nNDZGJav++ddfg/hIgYkJKj87C7fcnc26c6rv32UWqlXY9AGYmID16+Ho0TljNlnGmGow1sx+\nG9gOfNTdJ83sPJCs1QjwJpAsr5ctaGrqrTTV6opKZYTJyel+V6OrFGOPved9cGEGLmRUnzU3UY5n\n8rB+PZNrboK8xNoFufouF6F07CjliYnoycQEU4eOzI4Xpo1xvh+HjhO9mX2CaNB13N1/ERcfAb5i\nZiuBFcAo8CpwGHgwfn0LcLDjmotIa4nZY+X77op+RCT3ejnVtqNEb2ZLgV3AWeB5MwN40d3/2Mx2\nESXyJcB2d79oZk8D3zazQ8Bl4OOZ1l5EIvXZY8PD2R0pSHf1cHr30MxM/n79Jyenc1OpUA4TW1GM\n4ShCnIqx5fuGmpUX+hIIIiJFoEQvIhI4JXoRkcAp0YuIBE6JXkQkcEr0IdJ1T0QkQYk+NLqDkYg0\nUKLPgwxb4LqDkYg0UqLvt4xb4LodoYg0UqLvs8xb4Bnd21RyRuMusghK9H3WlRZ48ronMvg07iKL\nVOx7xuZBHu9bK7nS7KhPt7+UTqhFnwdqgUsLGneRxVKLflBVqzoKKAod9ckiKdEPorjPtn7DAg26\nFkD9qE8kBXXdDCDNlReRToSb6AOejqY+WxHpRJhdN6F3bajPVkQ6EGSLvhBdG5qpI9JawEf1nQoy\n0atrQ6TgdJLZHGF23ahrQ6TQdJLZXEG26AF1bYgUmI7q5wqzRS8ixaaj+jmU6EUkTDrJbFa4XTci\nIgK02aI3s7uBJ9193MxuBZ4BZoBXgW3u/o6ZbQUeAWrATnffY2bXAd8B1gDTwKfcfbILcYiIyDwW\nbNGb2WPAt4CVcdFTwA533wgMAQ+Z2Q3Ao8C9wGbgCTNbAXwWeCVe9llgR/YhiIhIK+206M8ADwP/\nM36+AXgxfrwX+DBwBTjs7peAS2Z2GrgNuA/4WmLZP2qnUuXyKkqlpW0F0AuVyki/q9B1ijEcRYhT\nMXZmwUTv7t81s5sTRUPuPhM/ngZWA9cD5xLLNCuvly1oauqtdhbriUplhMnJ6X5Xo6sUYziKEKdi\nbP2+ZtIMxr6TeDwCvAmcjx+3Kq+XiYhID6VJ9C+Z2Xj8eAtwEDgCbDSzlWa2GhglGqg9DDzYsKyI\niPRQmkT/BeBxM/shsBx4zt3fAHYRJfL9wHZ3vwg8DYyZ2SHgM8Dj2VRbRETaNTQzM7PwUj02OTmd\nm0qpPzAMRYgRihGnYmz5vqFm5TphSkQkcEr0IiKBU6IXEQmcEr2ISOCU6EVEAqdELyISOCV6EZHA\nKdGLiAROib6ZapXSsaOFv3O8iIRBtxJsVK1S3jxO6dTJ6ObCPzrW7xqJiCyKWvQNSn6C0qmT0eNT\nJ+H48T7XSERkcZToG9RsNGrJQ/R3bKzPNRIRWRx13TQaHmZq3wFKfoKajVIZHoYLYV9ASUTCpkTf\nzPAwtQ139rsWIiKZUNeNiEjglOhFRAKnRC8iEjglehGRwCnRi4gEToleRCRwSvQiIoFTohcRCZwS\nvYhI4JToRUQCl+oSCGa2DPg2cDNwBdgK1IBngBngVWCbu79jZluBR+LXd7r7nsVXW0RE2pW2Rf8g\nUHL3e4A/Ab4CPAXscPeNwBDwkJndADwK3AtsBp4wsxWLr7aIiLQrbaI/CZTMbAlwPfA2sAF4MX59\nL/AAcBdw2N0vufs54DRw2+KqLCIinUh79coqUbfNBPAu4DeBTe4+E78+Dawm+hE4l3hfvbylcnkV\npdLSlFVrVttqdAORsTEYHu747ZXKSHZ1ySnFGI4ixKkYO5M20f8BsM/dv2RmNwL7geWJ10eAN4Hz\n8ePG8pampt5KWa0mGm4NOLXvQEfJvlIZYXIy7OvRK8ZwFCFOxdj6fc2k7bqZ4mpL/RfAMuAlMxuP\ny7YAB4EjwEYzW2lmq4FRooHanmm8NWDJT/Ty40VE+i5tov86cIeZHSRqzX8Z2AY8bmY/JGrdP+fu\nbwC7iJL+fmC7u19cfLXb13hrwJqN9vLjRUT6LlXXjbtXgf/a5KX7myy7G9id5nMy0XBrwDR99CIi\ng6wYtxLUrQFFpMB0ZqyISOCU6EVEAqdELyISOCV6EZHAKdGLiAROiV46V61SOnY0urSEDB59f4Wj\nRC+diS8pUd7yQcqbx5UsBo2+v0JSopeO6JISg03fX4708MhKiV46oktKDDZ9fznR4yOrYpwZK9nR\nJSUGm76/XGh2ZNXNs/fVopfO1S8poSQxmPT99V2vj6zUohcR6bUeH1kp0YuI9EMPL7YYVtdNEecH\nFzFmEelIOIm+iPODixiziHQsmERfxPnBRYxZJFhdPDoPJtEXcX5wEWMWCVKXj87DGYwt4vzgIsYs\nEqDGo3OOH4f3vC+z9QfTogeKOT+4iDGLBKbx6JyxsUzXH06LXkRkUDUcnVeGh+HCdGarV6IXEcmD\nLs6rD6vrRkRErqFELyISOCV6EZHApe6jN7MvAb8FLAe+CbwIPAPMAK8C29z9HTPbCjwC1ICd7r5n\nsZUW6ZpqVdNVJTipWvRmNg7cA9wL3A/cCDwF7HD3jcAQ8JCZ3QA8Gi+3GXjCzFZkUG+R7OmSEhKo\ntF03m4FXgL8HvgfsATYQteoB9gIPAHcBh939krufA04Dty2qxiJdoktKSKjSdt28C3g38JvALcA/\nAEvcfSZ+fRpYDVwPnEu8r17eUrm8ilJpacqqZa9SGel3FbpOMQL33QXr18PEBKxfT/m+uway+0bf\nZRiyjDFtov85MOHulwE3s4tE3Td1I8CbwPn4cWN5S1NTb6WsVvYqlREmJ7M7cSGPFGPC9/df7aO/\nMJPpSSu9oO8yDGljnO/HIW3XzSHgI2Y2ZGa/CvwS8ELcdw+wBTgIHAE2mtlKM1sNjBIN1Irkky4p\nIQFK1aJ39z1mtokokS8BtgE/Bnab2XLgBPCcu18xs11ESX8JsN3dL2ZTdRERaUfq6ZXu/liT4vub\nLLcb2J32c0REZHF0wpSISOCU6EVEAqdELyISOCV6EZHAKdGLiAROiV5EJHBK9CIigVOil/ZUq5SO\nHQ37io5FiFEKSYleFlaEy/cWIUYpLCV6WVARLt9bhBiluJToZUE1G6W2dl30eO266MqOgSlCjFJc\nqa91IwUyPMzUvgNh32KvCDFKYSnRS3vql+8NWRFilEJS142ISOCU6EVEAqdELyISOCV6EZHAKdGL\n9JrOwJUeU6IX6SWdgSt9oEQv0kM6A1f6IaxEr0NiyTmdgSv9EM4JU/EhcenUSWpr1zG174DObpT8\n0Rm40gfBtOh1SCwDo34GrpK89EgwiV6HxCIizYXTdaNDYhGRphaV6M1sDXAM+BBQA54BZoBXgW3u\n/o6ZbQUeiV/f6e57FlXjVnRRKhGRa6TuujGzZcD/AC7ERU8BO9x9IzAEPGRmNwCPAvcCm4EnzGzF\n4qosIiKdWEyL/s+AvwK+FD/fALwYP94LfBi4Ahx290vAJTM7DdwGHG214nJ5FaXS0kVULVuVyki/\nq9B1ijEcwcdZrVJ5/TUYGwu6izbL7zFVojez3wUm3X2fmdUT/ZC7z8SPp4HVwPXAucRb6+UtTU29\nlaZaXVGpjDA5Od3vanSVYgxH8HFWq1Qe/ABMTAQ9jTrt9zjfj0Parpv/BnzIzA4AtwPPAmsSr48A\nbwLn48eN5SIiHSv5CZiYiB5rGnXbUiV6d9/k7ve7+zjwMvBJYK+ZjceLbAEOAkeAjWa20sxWA6NE\nA7WSJ0U+o7jIsQ+gmo3C+vXRY02jbluW0yu/AOw2s+XACeA5d79iZruIkv4SYLu7X8zwM2WxqlXY\n9AHKgR8KN6WzqQfP8DAcPcrUoSOaRt2BRSf6uFVfd3+T13cDuxf7OX0Tt/hC3amaHQoXZYpqs7Op\nixL7QNM06o4Fc2ZsV1SrcOedQV9StsiHwjqbWooinDNju6AQrd0iHwrrbGopCLXoWyhMa7fIF9kq\ncuyS3oAN4qtF30qRW7si0twADuKrRb+Q+Vp8A/aLPjC0XaWZHO0Xg3hJdCX6NHTfz+7QdpVmcrZf\nDOIgvhJ9CoP4iz4ItF2lmcb9guPH+1uheBB/au8LA9FtA0r0qQziL/og0HaVZhr3C8bG+lwjBm4Q\nX4OxaWhaXndou0ozDftFZXgYLgR84bYuUIs+rQH7RR8Yjds1R4Nw0kdZ/78VbL9Sopf8ytkgnASi\ngPtVeIm+YL/UIdPgrHRDEfersBJ9AX+pQ6bBWemGIu5XQQ3G6mqEgdHgrHRDAferoFr0RfylDp4G\nvaUb+rlfNete7nKXc1At+iL+UovIAGl2nRy4tizjG7yHlehBNyUQkdyabyD4mrJbfiXTzw2q60ZE\nJM+adS/3oss5vBa9iMxVrao7My/m6V7udpezEr1IyAbw2unBa9a93OUuZ3XdiASsiCcHybWU6EUC\nlqspxzprvW/bQF03IlnKW394XqYcqwupr9tALXqRrOT1Ehw5OOlMXUj93QapWvRmtgz4G+BmYAWw\nE3gNeAaYAV4Ftrn7O2a2FXgEqAE73X3P4qstkj+6BMf86l1I9dZsEc9a7+c2SNui/wTwc3ffCHwE\n+EvgKWBHXDYEPGRmNwCPAvcCm4EnzGzF4qvdhnb7wtRveFVI26IPseSqPzxvBvD2e5nr4zZI20f/\nv4Hn4sdDRK31DcCLcdle4MPAFeCwu18CLpnZaeA24GjqGrej3b4w9RteFdK2aBZLxqeUN9Vuf3je\n+vF7RWet920bpEr07l4FMLMRooS/A/gzd5+JF5kGVgPXA+cSb62Xt1Qur6JUWpqmapFX/x8kDqEr\n/38CVq2Cd78bTsT9YqOjsP/7bS1X+clPorLk37GxeX88OH68+XLVKhw9ChcuwHXXwZ3Z3jFn9nPr\nde9g/ZWfnZ27LX52Fm65e+HYICprtj3qyy+0rRrX02qdreKu1+n1ZrH8CpV2k339e4Jrt+FCMVVG\notPXW8W26QMwMRHV+YUX4N//vb1Y2zQb50J1bVfjerq13g60/V0mP2u+77Tf5tkOHcfYQupZN2Z2\nI/D3wDfd/X+Z2dcSL48AbwLn48eN5S1NTb2VtlpRa+7Tvzcb2MzSElc+/XuUfvw6M8uWMfT221F5\naRlDtbeZWbacobcvU3vvrbD1M5TOnJ6zHMuXw+XLs8vNLt+s1ZtoSV6zHFD+0CZKZ07PLl57761M\n/d9/XvxON+dzr9a93fVXKiNMrrmJcqL/cGrNTTCZuC9ns9jeeysApTOnr90eCx0hJF9PrqfVOlvG\nndjez//jNbFUgMnJNu4zWq3O+Z7mbMM0R4oN8Uw/+RTliYlouZ/8hBlbz9CVWmZHUZXKSBRnVkdo\njet5/h8pP/zR7NfbwXpmY+zks+b7Tvttnu3QcYyx+X4cUvXRm9l/Av4P8EV3/5u4+CUzG48fbwEO\nAkeAjWa20sxWA6NEA7VdU/ITlH78+uzzoSu12eezyRsYqr0dl13m/Nf/kuk//cbsjpBcjsuXZ5dL\n/m02ap4cjGtcruQn5iR5iP75sxh5n/u5V+ve0foX6D9sGtuZ07MxNW6PhWYYzHk9uZ4W62wdd2J7\n//Rs6r7Qxu8puQ3bnTUxb2xx2ZUbb5pdduhKra1YO5XVDI/G9az4p31dWW83Z6C0+k77rVfbIe1g\n7JeBMvBHZnbAzA4Qdd88bmY/BJYDz7n7G8AuoqS/H9ju7hcXX+351Wx0thUFULvlPbPPZ5Ytmy2f\nKUWPa2vXcemhh6ndfsfsQFpyOZYvj8vm/m022JYcjGtcrrFeELUsshiwm/u5V+ve8fpbTMNrGtt7\nb52NqXF7LDQwOef15HparLOtOtXfk3JK4TX7T2IbtjvYOm9sa9dRu/0OfvH9F6jdeGPzemckq4Hh\nxvVcemBzV9bbzYHrVt9pv/VqOwzNzMwsvFSPTU5OL65S1Sqll38EQO32O4Dol7P2azdROuVR+Vqj\n9NOzcwfE6oNkieXK79/A1MuvRWU/PXv173wDacl1NFv/yz+a7aOv3X5Htne1b4yxzfW3fZjYLDaY\nf2BxoUHH5OvJ9bRaZzt1avKejg6FG/efhi6ntuo2X2zN9rVW+1OH5sSZ1aBv43q6td42perWaPWd\n9luT7bCIrpuhZuVhJvoMpd3gg0QxhqMIcSrGlu9rmuh1ZqyISOCU6EVEAqdELyISOCV6EZHAKdGL\niAROiV5EJHC5nF4pIiLZUYteRCRwSvQiIoFTohcRCZwSvYhI4JToRUQCp0QvIhI4JXoRkcClvpVg\nSMzsS8BvEd0w5ZtENzl/BpghuiPWNnd/x8y2Ao8Q3Qx9p7vv6U+NO2Nmy4BvAzcT3bB9K1EMzxBA\njGZ2N/Cku4+b2a20GZeZXQd8B1hDdD/jT7n7ZF+CWEBDjLcDf0H0XV4CPunu/xZSjImyjwOfc/f3\nx8+DidHM1gC7iW7itJToezzTjRgL36KPb394D3AvcD9wI/AUsMPdNwJDwENmdgPwaLzcZuAJM1vR\nl0p37kGg5O73AH8CfIVAYjSzx4BvASvjok7i+izwSrzss0R3ScudJjH+OVHyGweeB74YYIyY2a8D\nnyb6Hgkwxq8Bf+vum4jqvL5bMRY+0RNtzFeIbnT+PWAPsIGoVQ+wF3gAuAs47O6X3P0ccBq4rffV\nTeUkUDKzJcD1wNuEE+MZ4OHE807iug/4QcOyedQY48fc/eX4cQm4SGAxmtkvA18FPp9YJqgYiZL5\nr5nZPwG/AxygSzEq0cO7gN8A/gvw34G/BZa4e/3aENPAaqIEeS7xvnr5IKgSddtMEB0q7gKGQojR\n3b9L9MNV10lcyfLcxtoYo7v/K4CZ3QP8PvB1AorRzJYCfw38IVGd64KJMXYzMOXuDwBngS/SpRiV\n6OHnwD53v+zuTtQ6Sm7EEeBN4Hz8uLF8EPwBUYzrgP9M1F+/PPF6CDHWvZN4vFBcyfKBitXMfhv4\nK+CjcV9tSDFuANYCTwN/B7zPzL5BWDFClHv+IX78PaIGZ1diVKKHQ8BHzGzIzH4V+CXghbjvHmAL\ncBA4Amw0s5VmthoYJRrsGwRTXG0N/AJYBrwUWIx1ncR1mGj8Irls7pnZJ4ha8uPu/npcHEyM7n7E\n3cfiMYiPAa+5++cJKMbYIa7WexNwnC7FWPhEH88qeYloA38P2AZ8AXjczH5I1PJ9zt3fIOryOAjs\nB7a7+8X+1LpjXwfuMLN63b9MFGdIMdZ18t09DYyZ2SHgM8Djfapz2+JujV1ErbrnzeyAmT0eUozz\nCTDGLwCfNLN/AT4CfLVbMeoyxSIigSt8i15EJHRK9CIigVOiFxEJnBK9iEjglOhFRAKnRC8iEjgl\nehGRwP0HCmdqXByfRqQAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f1fc307d4a8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "saves = 400\n",
    "# fig,plts = plt.subplots(2)\n",
    "# for i,s in enumerate(people):\n",
    "#     if i >= 1000:\n",
    "#         break\n",
    "#     index, start, track, des, flag = get_track(s)\n",
    "# #     track = np.vstack([track[0], track, track[-1]])#添加固定偏差来保证数据有效\n",
    "# #     des = np.array(des) + np.array(start)\n",
    "#     plt.plot(des[0], des[1],'.b')\n",
    "for i,s in enumerate(machine[-100:]):\n",
    "    if i >= saves:\n",
    "        break\n",
    "    index, start, track, des, flag = get_track(s)\n",
    "    track = np.vstack([track[0], track, track[-1]])#添加固定偏差来保证数据有效\n",
    "    delta = track[1:] - track[:-1]\n",
    "    delta[:,2]+=1#固定偏差用来消除0\n",
    "    \n",
    "    angular = np.arccos(delta[:,0]/np.sqrt(np.square(delta[:,0])+np.square(delta[:,1])))\n",
    "    angular = np.nan_to_num(angular)\n",
    "#     des = np.array(des) + np.array(start)\n",
    "#     track = np.vstack([track[0], track, track[-1]])#添加固定偏差来保证数据有效\n",
    "    plt.plot(des[0], des[1],'.r')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 105,
   "metadata": {
    "collapsed": true
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/miller/anaconda3/lib/python3.6/site-packages/ipykernel/__main__.py:19: RuntimeWarning: invalid value encountered in true_divide\n"
     ]
    }
   ],
   "source": [
    "from random import shuffle\n",
    "import os\n",
    "# shuffle(machine)\n",
    "saves = 10000\n",
    "\n",
    "dirs = ['data/people/speed_grad', 'data/people/angle_grad', 'data/people/speed', 'data/people/angle',\n",
    "       'data/machine/speed_grad', 'data/machine/angle_grad', 'data/machine/speed', 'data/machine/angle',\n",
    "       'data/people/y_track', 'data/machine/y_track']\n",
    "for dir_ in dirs:\n",
    "    if not os.path.isdir(dir_) and not os.path.isfile(dir_):\n",
    "        os.makedirs(dir_)\n",
    "        print(\"mkdir %s\"%dir_)\n",
    "    \n",
    "for i,s in enumerate(machine):\n",
    "    if i >= saves:\n",
    "        break\n",
    "    index, start, track, des, flag = get_track(s)\n",
    "    delta = track[1:] - track[:-1]\n",
    "    theta = np.arccos(delta[:,0]/np.sqrt(np.square(delta[:,0])+np.square(delta[:,1])))\n",
    "    theta_speed = (theta[1:] - theta[:-1])/delta[1:,2]\n",
    "    speed = np.sqrt(np.square(delta[:,0])+np.square(delta[:,1]))/delta[:,2]\n",
    "    speed_grad = (speed[1:] - speed[:-1])/delta[1:,2]\n",
    "#     print(speed)\n",
    "#     canvas = plt\n",
    "#     canvas.plot(*start, '.b')\n",
    "#     plt.figure(figsize=(5,5),dpi=20,edgecolor='k')\n",
    "#     plt.plot(track[1:,2],theta,'-b',linewidth=10)\n",
    "#     plt.axis('off')\n",
    "#     plt.savefig('data/people/angle/%s.png'%index)\n",
    "#     plt.close()\n",
    "#     plt.plot(track[1:,2],speed,'-b',linewidth=10)\n",
    "#     plt.axis('off')\n",
    "#     plt.savefig('data/people/speed/%s.png'%index)\n",
    "#     plt.close()\n",
    "#     plt.plot(track[2:,2],theta_speed,'-b',linewidth=10)\n",
    "#     plt.axis('off')\n",
    "#     plt.savefig('data/people/angle_grad/%s.png'%index)\n",
    "#     plt.close()\n",
    "#     plt.plot(track[2:,2],speed_grad,'-b',linewidth=10)\n",
    "#     plt.axis('off')\n",
    "#     plt.savefig('data/people/speed_grad/%s.png'%index)\n",
    "#     plt.close()\n",
    "\n",
    "    plt.figure(figsize=(5,5),dpi=20,edgecolor='k')\n",
    "    plt.plot(track[:,2],track[:,1],'-b',linewidth=10)\n",
    "    plt.axis('off')\n",
    "    plt.savefig('data/people/y_track/%s.png'%index)\n",
    "    plt.close()\n",
    "    \n",
    "\n",
    "#     misc.imsave('data/machine/lena_new_sz.png', track)"
   ]
  }
 ],
 "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.6.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}