jaimeide/tutorial_dpcca_computation.ipynb

## tutorial_dpcca_computation.ipynb
{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Detrended Partial Cross Correlation in Python\n",
    "**Author: Jaime Ide,\n",
    "Date: Feb 15, 2018**\n",
    "\n",
    "Simple code to compute the detrended partial cross correlation analysis (DPCCA) coefficient. It is simple, but since I didn't find any code in Python, I decided to post it here. In addition to DPCCA, I also compute the DCCA, correlation and partial correlation matrixes as well so that you can compare them. Credits to Pitithat Puranachot who posted a version of the DCCA in R (CrossValidated). I have the Matlab code version too. Just send me a message if you need it.\n",
    "\n",
    "- **References:** \n",
    "\n",
    "[1] Ladislav Kristoufek. _Measuring correlations between non-stationary series with DCCA coefficient_, Physica A: Statistical Mechanics and its Applications, Volume 402, Pages 291-298, ISSN 0378-4371, 2014.\n",
    "https://doi.org/10.1016/j.physa.2014.01.058.\n",
    "\n",
    "[2] Yuan, N. et al. _Detrended Partial-Cross-Correlation Analysis: A New Method for Analyzing Correlations in Complex System_, Scientific Reports, 5, Issue 8143, 2015.\n",
    "\n",
    "- **If you want to give credit for the code or know more how DPCCA can improve supervised learning, please reference my work at NIPS2017:**\n",
    "\n",
    "[3] Ide,JS; Cappabianco,FA; Faria,FA; Li,CSR. _Detrended Partial Cross Correlation for Brain Connectivity Analysis_. Advances in Neural Information Processing Systems 30, pages: 889--897, 2017 http://papers.nips.cc/paper/6690-detrended-partial-cross-correlation-for-brain-connectivity-analysis.pdf\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 1. Computing DPCCA (step-by-step)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7f63453cfa20>"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAAEICAYAAABLdt/UAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzsnXV8VMf6h5/ZuAfiIUGCuwUtxaHFKVahSg1+dbuVW6G9\nLXU3bm+9pYUWd4eWAgkEglsgEBKIEvfs7vz+OBsIkITIZjcb5mn3s3LOmXk3Id+dfecVIaVEoVAo\nFA0HnbUNUCgUCoV5UcKuUCgUDQwl7AqFQtHAUMKuUCgUDQwl7AqFQtHAUMKuUCgUDQwl7AqbQwgx\nWwjxq5nGaiqEyBVC2JljPEthsjnM2nYo6idK2BVVRggxQAixQwiRJYRIF0JsF0L0srZd1UEIcUYI\nMbz0uZTyrJTSXUppsKZd1cVkc6y17VDUT+ytbYDCNhBCeAIrgVnAH4AjcCNQZE27rjeEEPZSSr21\n7VDUb9SKXVFV2gBIKX+XUhqklAVSyvVSygMAQoiWQojNQogLQog0IcQ8IYR36cWmlfJzQogDQog8\nIcR3QogAIcQaIUSOEGKjEKKR6dzmQggphHhICHFeCJEohHimIsOEEH1N3yQyhRD7hRCDKzjvF6Ap\nsMLkyvhXmbnsTedsFUK8aRovVwixQgjhY3o/2UKI3UKI5mXGbCeE2GD6BnNcCDGtEjvvFULEmt7v\naSHE9DLHZgghjgohMoQQ64QQzcock0KIR4QQMUBMmddamR47CSE+EEKcFUIkCyHmCiFcTMd8hRAr\nTT+bdCHENiGE+rtv6Egp1U3drnkDPIELwE/AKKDRFcdbASMAJ8AP+Bv4pMzxM0AEEAA0AVKAvUB3\n0zWbgddM5zYHJPA74AZ0BlKB4abjs4FfTY+bmOwajbZQGWF67lfB+zhTOs4Vc9mbnm8FTgItAS/g\nCHACGI72Dfdn4AfTuW5APHCf6VgPIA3oWM68bkA20Nb0PKj0PGCiac72pnFeBnaUuVYCG4DGgEuZ\n11qZHn8CLDcd9wBWAG+bjr0NzAUcTLcbAWHtf0/qVrc39cmtqBJSymxgAJqg/A9IFUIsF0IEmI6f\nlFJukFIWSSlTgY+AQVcM87mUMllKeQ7YBkRKKaOllEXAEjSRL8vrUso8KeVB4Afg9nJMuxNYLaVc\nLaU0Sik3AFFoQl9TfpBSnpJSZgFrgFNSyo1Sc4H8WcbOscAZKeUPUkq9lHIvsAiYUsG4RqCTEMJF\nSpkopTxsev1hNCE+appjDtCt7KrddDxdSllQdkAhhAAeBJ4yHc8xXX+b6ZQStA+RZlLKEinlNiml\nKhDVwFHCrqgyJuG5V0oZAnQCgtFWiwgh/IUQ84UQ54QQ2cCvgO8VQySXeVxQznP3K86PL/M4zjTf\nlTQDpppcDZlCiEy0D6Cgar69mtjZDOhzxdzTgcArB5RS5gG3AjOBRCHEKiFEuzLjfFpmjHRAoH0b\nKSWe8vEDXIE9Za5fa3od4H20bwPrTW6gF6rw/hU2jhJ2RY2QUh4DfkQTeNC+8kugi5TSE20lLWo5\nTWiZx02B8+WcEw/8IqX0LnNzk1K+U5HptbTpyrn/umJudynlrHInlnKdlHIE2ofOMbRvPqXjPHzF\nOC5Syh1VsDsN7cOmY5lrvaSU7qY5c6SUz0gpw4BxwNNCiGG1f+uK+owSdkWVMG0SPiOECDE9D0Vz\njUSYTvEAcoFMIUQT4DkzTPuKEMJVCNERzY+9oJxzfgXGCSFuEkLYCSGchRCDS+0sh2TAXPHfK4E2\nQoi7hBAOplsvIUT7K080bRSPF0K4oUUS5QKlIZZzgRdN7xMhhJcQYmpVDJBSGtE+ID4WQvibrm8i\nhLjJ9HisEKKVyWWTbZrTpkI7FdVHCbuiquQAfYBIIUQemqAfAkqjVV5H2zzMAlYBi80w519oboRN\nwAdSyvVXniCljAcmAC+hbbDGo32oVPRv+23gZZPb4tnaGGfyZ49E82efB5KAd9E2g69Eh/azOo/m\nahkE/J9pnCWm6+ab3FiH0Daoq8rzaD+nCNP1G4G2pmOtTc9zgZ3AV1LKrdUYW2GDCLWPoqhvmMIJ\nTwMOUsVsKxTVRq3YFQqFooGhhF2hUCgaGMoVo1AoFA0MtWJXKBSKBoZVioD5+vrK5s2bW2NqhUKh\nsFn27NmTJqX0u9Z5VhH25s2bExUVZY2pFQqFwmYRQsRV5TzlilEoFIoGhhJ2hUKhaGAoYVcoFIoG\nhuqgpFAoGiQlJSUkJCRQWFhobVOqjbOzMyEhITg4ONToeiXsCoWiQZKQkICHhwfNmzdHq4FmG0gp\nuXDhAgkJCbRo0aJGYyhXjEKhaJAUFhbi4+NjU6IOIITAx8enVt80lLArFIoGi62Jeim1tfu6EPZ9\nKfs4kHrA2mYoFAqFRTCLsAshvIUQC4UQx0yd1vuZY1xzoDfqeWbrM7y7611rm6JQKBQWwVybp58C\na6WUU4QQjmg9GOsF289tJ6UgBSNGa5uiUCgUFqHWK3YhhCcwEPgOQEpZLKXMrO245mJRzCIA0grS\nKDYUW9kahUJxvTFkyBA2bNgAwMsvv8zjjz9e53OaY8UehtaS7AchRFdgD/CEqSv7RYQQDwEPATRt\n2tQM016b1PxU/k74myC3IBLzEknKS6Kpp2XmVigU9YfXVxzmyPlss47ZIdiT18Z1vPbcr7/Oq6++\nSkpKCtHR0SxfvtysdpSHOXzs9mi9Lr+WUnYH8oAXrjxJSvmNlDJcShnu53fN4mRmYdmpZRikgQe7\nPAhAYl6iReZVKBSKUgYOHIiUko8++oj58+djZ2dHbGws999/P1OmTKmTOc2xYk8AEqSUkabnCylH\n2C2NURpZHLOY8IBw+gb1BZSwKxTXK1VZWdcVBw8eJDExEV9fXzw8PAAICwvju+++qzNhr/WKXUqZ\nBMQLIUq7og8DjtR23NoSlRRFfE48k1pPIsA1AIFQwq5QKCxKYmIi06dPZ9myZbi5ubFu3TqLzGuu\nOPbHgHlCiANAN2COmcatMYtiFuHh4MGIZiNwtHPE18WXpLwka5ulUCiuE/Lz85k0aRIffvgh7du3\n55VXXmH27NkWmdsswi6l3Gfyn3eRUk6UUmaYY9yaklWUxca4jYwJG4OzvTMAQW5BnM89b02zFArF\ndYSrqys7d+5kxIgRgOZr37lzJwAXLlxg5syZREdH8/bbb5t97gZZBGxl7EqKjcVMaXPJfxXoFsiJ\njBNWtEqhUCg0fHx8mDt3bp2N3+BKCkgpWRSziI4+HWnbuO3F14Pdg0nMS0RKaUXrFAqFou5pcMJ+\nKO0QMRkxTGo96bLXA90CKTIUkVFkVS+RQqFQ1DkNTtgXxSzCxd6F0S1GX/Z6kFsQAIm5KjJGoVA0\nbBqUsOeX5LPm9BpGNhuJu6P7ZccuCrsKeVQoFA2cBiXs686sI1+ff9mmaSlK2BUKxfVCgxL2RTGL\nCPMKo6tf16uOeTl54WLvooRdoVA0eBqMsJ/MOMn+1P1Maj2p3O4jQgiC3IJUkpJCoWjwNBhhXxSz\nCHudPeNajqvwHJWkpFAorgcahLAXG4pZEbuCoaFDaezcuMLzAt0ClStGoVA0eBqEsG8+u5msoiwm\nt55c6XlBbkGkF6ZTqK9592+FQqGoDrbaaMPqLIpZRLBbMH2D+1Z6XrB7MADJ+ck082xmCdMUCkV9\nYM0LkHTQvGMGdoZR71zzNFtttGFVEnISiEiM4JbWt6ATlb+dQLdAQIU8KhQKy1Feo42lS5fy4IMP\nMmHCBNavX2/2OW1+xb44ZjE6oWNiq4nXPFdlnyoU1ylVWFnXFeU12pg4cSITJ04kIyODZ599lpEj\nR5p1TpteseuNepadXMYNwTdcXI1Xhmq4oVAoLMm1Gm28+eabPPLII2af16aFffu57aQUpFxz07QU\nBzsH/Fz9lLArFIo6p7JGG1JKnn/+eUaNGkWPHj3MPrdNu2IWxSzCx9mHgaEDq3xNkFuQEnaFQlHn\nlDbaKKVso43PP/+cjRs3kpWVxcmTJ5k5c6ZZ57ZZYU/NT+XvhL+5p+M9OOgcqnxdkFsQR9OP1qFl\nliUhJ4Ffj/7K9HbTCfUMtbY5CoWiCjz++ON1GvZos66YZaeWYZCGq+quX4sgtyAScxMxSmMdWWZZ\n5h+bz7yj85i4bCJf7ftKxegrFArbFHajNLI4ZjHhAeHVjkcPdAuk2FhMemF6HVlnWSISI+jo05Fh\nTYfx9f6vmbhsIn/F/2VtsxQKhRWxSWGPSooiPie+2qt1uJSk1BCKgaUXpnM84zhDmw7lvUHv8e3I\nb3G0c+TRzY/y2KbHSMhJsLaJCoXCCtiksC+KWYSHgwcjmo2o9rUNqS77rsRdAPQN0jJu+wT1YdG4\nRTzV8ykikyKZuGwic/fPpchQZE0zFQqFhbE5Yc8qymJj3EbGhI3B2d652teXxrs3hCqPEYkRuDu4\n08Gnw8XXHOwcmNFpBssnLmdQyCC+3Pclk5ZN4p9z/1jRUoVCYUlsTthXxq6k2FhcbpekquDp6Imr\nvWuDcMVEJEbQK7AX9rqrg5sC3QL5cPCH/HfEf9EJHbM2zuLJLU82iA80hUJROTYl7FJKFsUsoqNP\nR9o2blujMUobbti6KyY+J55zuefoE9Sn0vP6B/dn0fhFPNHjCbaf286EpRP434H/UWwotpClCoXC\n0tiUsB9KO0RMRkyNNk3LEuRu+8IemRgJQL+gftc819HOkQc6P8DyicsZ0GQAn0V/xuTlk9mdtLuu\nzVQoFFbApoR9UcwiXOxdGN1idK3GaQgt8iITI/Fz8aOFV4sqXxPkHsTHQz7m6+FfozfqeXzz45QY\nS+rQSoVCYQ1sKvP03o730je4L+6O7rUap7ThRoG+ABd7FzNZZzmM0khkYiQDmgwot7/rtRjQZABP\n9nySZ/96liMXjpTb/FuhaEi8u+tdjqUfM+uY7Rq34/nez1/zvCFDhvDSSy8xYsQIXn75ZbKzs/ns\ns8/MasuV2JSwN/dqTnOv5rUepzQyJikvqVor3vpCTEYMGUUZ1/SvV0bPgJ4A7Eneo4RdoahDrNFo\nw6aE3VyUjWW3RWGPSIwAqJWw+7r40sKrBVFJUczoNMNcpikU9ZKqrKzrirKNNrZu3YqdnR1Hjx7l\n008/JS0tjWHDhjFr1iyzzmlTPnZzYevZpxGJETT3bF6lGvSVER4QTnRKNAajwUyWKRSKKylttOHk\n5HSx0Ub79u2ZO3cuf/zxB1FRUWaf87oUdj9XP3RCZ5ORMSWGEvYk77mYbVobwgPCyS3J5ViGeX2P\nCoVCo7JGG8uXL2fAgAEMGzbM7PNel8LuoHPAz8XPJpN1DqQdoEBfYBZhL/WzRyWZf8WgUFzvVNZo\nA2D8+PHs2LGDefPmmX3u69LHDrYb8hiRGIFO6AgPDK/1WAFuAYR6hBKVHMU9He8xg3UKhaKUyhpt\nbN26lcWLF1NUVMTo0bUL3y6P61rYD104ZG0zqk1kYiQdGnfAy8nLLOOFB4Sz6ewmjNKITlyXX+AU\nCoszePBgBg8eXGfjX7d/yUHu2ordlhpu5JXkcTD1IH2Da++GKSU8MJzs4mxiMmLMNqZCobAuZhN2\nIYSdECJaCLHSXGPWJUFuQZQYS2yq4cae5D3opb5WYY5XEh6guXSikpWfXaFoKJhzxf4EYDPNREtj\n2W1pAzUiMQInOye6+3c325jB7sEEuwWzJ3mP2cZUKBTWxSzCLoQIAcYA35pjPEtQGgNuSyGPEYkR\ndPPvhpOdk1nHDQ8MZ0/yHqSUZh1XobA2tvpvurZ2m2vF/gnwL6BCh7UQ4iEhRJQQIio1NdVM09Yc\nW0tSSitIIyYjxixhjlcSHhBOemE6sVmxZh9bobAWzs7OXLhwwebEXUrJhQsXcHaufiOhUmodFSOE\nGAukSCn3CCEGV3SelPIb4BuA8PBwq/+kPRw9cHdwt5kV+5Vt8MzJRT97UhQtvVuafXyFwhqEhISQ\nkJBAfVhIVhdnZ2dCQkJqfL05wh1vAMYLIUYDzoCnEOJXKeWdZhi7Tgl0C7QZH3tkUiQejh60b9ze\n7GOHeITg7+pPVHIUt7a71ezjKxTWwMHBgRYtbK8WlDmotStGSvmilDJEStkcuA3YbAuiDraTpCSl\nJOJ8BL0De2OnszP7+EIIegb0JCo5yua+tioUiqu5buPYAZtpkZeQk8D5vPN14oYpJTwgnLSCNOKy\n4+psDoVCYRnMKuxSyq1SyrHmHLMuCXIPIrMok/ySfGubUik7E7U0ZHPGr19JaYkCFc+uUNg+1/2K\nHSApv367YyITIwlwDaC5Z/M6m6OFZwt8nH2UsCsUDQAl7EBibv11xxilkV1Ju+gT1KdGbfCqykU/\ne5LysysUto4Sdup3ktLx9ONkFmXWqX+9lPDAcJLzkzmXe67O51IoFHXHdS3sfq5+2Am7ei3s5miD\nV1VU3RiFomFwXQu7vc4ef1f/eh3yGJkYSZhXGP6u/nU+V0vvlng7eavGGwqFjXNdCzvU75DHYkOx\n2drgVQWd0NHDv4dasSsUNs51L+z1Oft0f+p+Cg2FFnHDlBIeGM653HP1+luMQqGonOte2IPcgkjO\nT8ZgNFjblKsobYPXK7CXxeYs9bPvTtptsTkVCoV5UcLuFoTeqOdC4QVrm3IVkYmRdPLphIejh8Xm\nbNOoDR4OHqo+u0Jhwyhhd6+fIY+5xbkcSjtkUTcMgJ3Ojh4Bys+uUNgyStjraSx7VHIUBmmw2MZp\nWcIDwonLjiM13/bKnSoUCiXs9Tb7NCIxAmc7Z7r6d7X43KV1Y5Q7RqGwTa57YXd3dMfDwaPerdgj\nEyPp7t/d7G3wqkK7xu1wc3BT7hiFwka57oUdINA9sF4Je2p+KiczT9I32PJuGNASt7r5d1OJSgqF\njaKEHQh2C65XcduRSZGAZcoIVER4QDinsk6RXphuNRsUCkXNUMJO/UtSijgfgZeTF+0atbOaDaXx\n7MrPrlDYHkrY0TZQs4uzySvJs7YpSCmJTIqsszZ4VaWjT0ec7ZyVO0ahsEGUsFOm4UY9cMfEZceR\nlJdklTDHsjjYOdDVv6vaQFUobBAl7NSvJKXIROv710sJDwgnJiOGrKIsa5uiUCiqgRJ26leSUkRi\nBEFuQTT1aGptUwgPCEci2Zu819qmKBSKaqCEHfBzMTXcsHKSksFosEgbvKrS2a8zjjpH5Y5RKGwM\nJexo9VECXAOsvmI/ln6M7OJsq/vXS3Gyc6KLXxcl7AqFjaGE3USgm/WTlCzZBq+q9AzoybH0Y+QU\n51jbFIVCUUWUsJsIdrd+klJEYgStvFvh6+JrVTvKEh4YjlEaiU6JtrYpCoWiiihhNxHkFkRynvUa\nbhQZiohOia43bphSuvp1xV5nr9wxCoUNoYTdRKBbIHqpJ7XAOqVqD6QeoMhQVK/cMAAu9i508unE\nniSVgapQ2ApK2E1YO0lpT/IeBILu/t2tMn9lhAeGc/jCYfJL8q1tikKhqAJK2E1YO5Y9OiWaVo1a\n4eXkZZX5KyM8IByDNLAvZZ+1TVEoFFVACbsJa2af6o169qXso4d/D4vPXRW6+XfDTtgpP7tCYSMo\nYTfh5uCGp6OnVZKUYjJiyNfn10s3DGg/mw4+HVSlR4XCRlDCXoYgtyCrrNj3pmgp+z0Delp87qoS\nHhDOwbSDFOoLrW2KQqG4BkrYy2A1YU/eS5BbEIFugRafu6r0DOhJibGEA6kHrG2KQqG4BkrYy2CN\n7FMpJdEp0fXWDVNK94DuCITysysUNoAS9jIEuweTU5xDbnGuxeZMyE0gtSC13m6cluLp6Em7xu2U\nsCsUNoAS9jJYI+SxNFW/e0D9XrGD5o45kHqAYkOxtU1RKBSVoIS9DKU+bksK+97kvXg4etDKu5XF\n5qwp4YHhFBmKOJh20NqmKBSKSqi1sAshQoUQW4QQR4UQh4UQT5jDMGtgjezTvSl76e7fHZ2o/5+x\nPf21qB3VB1WhqN+YQ030wDNSyvZAX+ARIUQHM4xrcXxdfLEX9hZbsacXpnM663S93zgtxdvZm9aN\nWrMraZe1TVEoFJVQa2GXUiZKKfeaHucAR4EmtR3XGtjp7Ahws1zDjdIU/fq+cVqWIaFDiEqOIq0g\nzdqmKBSKCjDr938hRHOgOxBpznEtSZBbkMWyT6NTonHQOdDRt6NF5jMHY1qMwSiNrD291tqmKBSK\nCjCbsAsh3IFFwJNSyuxyjj8khIgSQkSlplqnNG5VsGSS0t6UvXT27YyTnZNF5jMHYd5htG/cntWn\nV1vbFIVCUQFmEXYhhAOaqM+TUi4u7xwp5TdSynApZbifn585pq0TAt0CSclPQW/U1+k8BfoCjqQd\nqb1/PSsBdn4Jx9dAdiJIaR4DK2F0i9EcTDtIXHZcnc+lUCiqj31tBxBCCOA74KiU8qPam2RdgtyD\nMEgDaQVpdZrifyjtEHqpp0dALfzrifth3lTITb70mps/BHXVbsHdtHuvUBCi9kabuLnFzXy05yNW\nn17NrK6zzDauQqEwD7UWduAG4C7goBCitGD3S1JKm/yuHuwWDGix7HUp7HuTtcJfXf261myAkxvh\nj3vApRE8uAUMxZrQJ+6H8/vg1GaQpjZ/Lo1MYt/tkug3agG6mn1hC3QLJDwwnNWxq5nZZSbCjB8a\nCoWi9tRa2KWU/wAN5i+7NJb9fO75Og1DjE6JppV3DRtr7P0FVjwBAR3gjj/BU7OZpmX6pZYUQPIR\nSNxnuu3XXDbGEu24kye0HArjPtGEv5qMaTGG2TtncyT9CB19bGfzV6G4HjDHir1BYYnsU4PRwL7U\nfYwNG1u9C6WEv96FrW9rojztZ3DyKP9cBxcI6andStEXQ+pRTeTP7YXoXyHlCNzxBzRuUS1Thjcb\nzluRb7EqdpUSdoWinlH/0x0tjKuDK15OXnWafXoi4wR5JXnV+0ZgKIHlj2qi3m26JsYViXpF2Dtq\nbpged2sr9buXQm4KfDsc4quXdOTl5MWNTW5k7em1GIyG6tmhUCjqFCXs5VDXIY+ljTWqnJhUlAO/\n3aqtsAc9DxO+BDuH2hvSfAA8sAmcPeHHsXBoUbUuHx02mtSCVHYn7669LQqFwmwoYS+Huhb26JRo\nAt0CL/ZZrZScJPhhNMRuhfGfw5CXzBrhgm8ruH8jNOkBC2fA3+9XOWRyUMgg3BzcWB1rk/vkCkWD\nRQl7OdRl9qmUkujkKjbWSD0O346AC6fgjgWaC6UucPOBu5dB52mw+U1Y+n+aP/4aONs7M6zpMDbE\nbaDIUFQ3tikUimqjhL0cgtyCyC3JJac4x+xjn8s9R0pByrXdMHE74LuRoC+E+1ZB6xFmt+Uy7J1g\n0jcw+EXY/xv8Ogny06952ZiwMeSW5LItYVvd2qdQKKqMiooph0D3S5ExHo7V3KC8Bhf965UlJh1e\nAosfAu9mcOdCaNS8xvOVGIwkZxdyPrOQxKyCMvcFBHg6858JndDpTK4dIWDwC9A4DJY9At+NgOl/\nas8roHdgb3ycfVgVu4rhzYbX2E6FQmE+lLCXQ2mSUlJeEm0atTHr2HuT9+LhUEljjZ1fwrqXILQv\n3P47uDaudDwpJYfPZ5OQkc+5zEISMwtIzCrkvEm8U3KKrnKZe7k44OPuyMajKbQJ8OCe/s0vP6HL\nNPAKgfnT4X/DNDvKxsiXwV5nz80tbubP43+SXZyNp6NnFX8SCoWirlDCXg5lk5TMTXRKNN38u13d\nWMNo1AQ98mtoPx4m/Q8cnCsdq1hv5Kk/9rHqwKX9AGcHHcFeLgR7uzCwtR9B3i408XYmyMuFYNO9\nm5M9Ukru+3E376w5xuC2fjTzcbt88Gb94YGNWsmCn8bBhK+gy9Ry7RjTYgzzjs5jU9wmbml9S41+\nLgqFwnwoYS8HHxcf7HXmb7iRUZhBbFYs41qOu/rg3h81Ue8zC26ac810/8ISA7N+3cOW46k8Obw1\nw9sH0MTbBW9Xhyql+AsheHtSZ0Z+/DfPLTzA/Af7XnLJlOLTUhP3BXfC4gcgPRYG/euqqJxOvp0I\n9Qhl1elVStgVinqA2jwtB53QEegaaHZhL22scVVETEEGbPoPNBsAN799TVHPKSzhnu93sfVEKnNu\n6cyTw9vQqYkXjdwcq1W3JcjLhVfGdmDX6XR+2nmm/JNcG8NdS6DLbbB1DiyZCfrLI2CEEIwJG8Ou\nxF2k5KdUeX6FQlE3KGGvgCD3ILNnn+5N2YuDzoFOvp0uP7DlbSjMhFHvXDNGPSOvmDu/jWRPXAaf\n3NqNO/o0rZVNU3uGMKStH++uPcaZtLzyT7J3glvmwpCX4cB8bWP1Cka3GI1EqgYcDQwpJWezz2KU\nRmuboqgGStgroC6SlPam7KWTb6fLG2skH4bd30L4DAjsXOn1KdmF3PZNBEeTcph7Z08mdKt9B0LN\nJdMFBzsd/1p4AKOxguQkIWDQc9D3/+DQYsi7vDVeC68WdPDpwKrTq2ptk6J+YDAaeHvX24xZMoZp\nK6ax+exmpAXq/StqjxL2CghyCyIlP4WS0mqItaRAX8CRC1c01pAS1jyv1XwZ8u9Kr49Pz2fqf3cS\nn5HPj/f2YniHALPYBRDo5cyrYzuw60w6P+44U/nJ3e/UygEfXnLVodEtRnPkwhFOZ502m20K65Bf\nks+TW5/k92O/MyZsDAX6Ap7Y8gS3rryVv+L/UgJfz1HCXgFBbkEYpZHUfPO08TuUdgi9UX95YtLR\n5XBmGwx9udKwxlOpuUz7704y8or59YE+9G/laxabyjKlZwhD2/nz3rpjnK7IJQMQ0BH8O8KBBVcd\nGtViFAKh2ubZOGkFacxYN4O/E/7mpT4v8c6N77Bs4jLevOFNcopzeHTzo9yx6g7+OfePEvh6ihL2\nCigNeTSXOyY6JRqAbv7dtBeK82HdvyGgE/S8r8LrDp/PYtrcnZQYjCx4uB89mla/dnpVKI2ScbTT\n8a+F+yt2yYAW556wWyt1UAZ/V396B/Zmdexq9Qdvo8RmxjJ91XRis2L5dMin3N7udkDLV5jQagLL\nb1nO6/1fJ70wnVkbZ3HXmrvYcX6H+n3XM5SwV0DZ7FNzsDdl7+WNNXZ8BlnxMOpdsCs/6nRPXDq3\nfROBk73PNTyXAAAgAElEQVSOPx7uR/uguk3+CfB05rVxHdl9JoMfKnPJdJ4CCDi48KpDY8LGcDbn\nLIcvHK4zOxV1w+6k3dy55k6KDEX8cNMPDA4dfNU5DjoHJrWexMpbVvJK31dIzk/m4Q0Pc+/ae9mV\nWL3Sz4q6Qwl7BZSu2M0RGWMwGtifsv+SGybzLPzzMXS8RSudWw7/xKRx57e78HV34s9Z/Qnzc6+1\nHVVhUo8mDGvnz3trjxGbmlv+SV4hmt0HFlxVCXJYs2E46BxYFas2UW2JFadW8NCGh/Bz8WPemHl0\n9K28eYqDnQPT2k5j1S2r+Heff5OQk8D96+9nxroZRCVFWchqRUUoYa8AF3sX/Fz82H5ue62/ZsZk\nxpBbkkv3ANPG6fpXAAEj/lPu+esPJzHjx90083Hlj4f70cTbpVbzVwchBHMmdcbJXsdzCw9gqMgl\n02UapJ/SOjGVwdPRk0Ehg1hzeo1qwGEDSCmZu38uL/3zEj38e/DL6F9o4l71aCtHO0dua3cbqyev\n5oXeL3A66zT3rbuPB9Y/cDFvQ2F5lLBXwsyuM4lKjmJhzNUuh+pQ2ri6h38POP03HFkKNz4N3qFX\nnbskOoFZ8/bSIdiT+Q/1xc/D6apz6poAT2dmj+/InrgMftheQYRL+/Fg5wQH/7jq0Oiw0VwovEBk\nUmQdW6qoDSXGEl7b8Rpf7vuScWHjmDt8bo1r/TjZOTG9/XTWTFrDc+HPEZMRw91r7mbT2U1mtlpR\nFZSwV8LUNlPpE9iHD6M+rFV99uiUaAJcAwhy9tPCG72bQv/Hrjrv14g4nv5jP72bN+bXB/rg7epY\nG/NrxS3dmzC8vT/vrzvOqfJcMi7e0OYmreuSQX/ZoYEhA3F3cFcNOOoxOcU5PLLxEZacXMLMrjN5\na8BbOJihK5ezvTN3d7ybNZPW0Mm3Ey9ue5GjF46awWJFdVDCXglCCGb3n41RGpm9c3aNXDJSSvYm\n76VHQA/E3h+15tEj39KaTZdhzcFEXl56iKFt/fnhvl64O1m3jI8Qgjm3dMbZwY7n/txfvkumy62Q\nl6p1dyqDk50Tw5sNZ+PZjRTqCy1jsKLKJOUlcc/ae9idtJs3+r/BI90eqbgUhdEIMRsg5ShUw7Xm\n6uDKZ0M/w8vJi0c3P6pKTVgYJezXIMQjhKd6PsWO8ztYcvLqpJxrcbGxhndbrTtRi0HQ/vIiYMnZ\nhby45CBdQrz4+s6eODvYmcv8WuHv6czs8R3YezaT7/8pxyXTegQ4e5cb0z66xWjySvL4O+FvC1h6\nfSClZH/qft6OfJvZO2bz1b6vWHRiEdsStnEi4wRZRVnXXHwcSz/G9FXTScxN5KvhX127aFvEVzBv\nCnzVF95pCj+M0faIDi+BjLhK2yj6uvjyxdAvyC3O5bHNj1GgL6jJ21bUAFXdsQrc2vZW1p9Zz/u7\n36d/cH8C3QKrfG1p/Hr32J1aU+pR715WD0ZKyXMLD1BYYuDjW7vhaF+/PmsndmvCqgNJfLD+OEPa\n+dPKv0x0jr0TdJwIB/6AolxwunSsd2BvfF18WRW7ipHNR1rB8oZDcl4yK2JXsOzkMs5kn8HJzgl3\nB3fSC9ORXC6sznbO+Lv6E+AWgL+rv/bYNYAA1wAK9AW8GfEmnk6e/DTqp2v3GrhwCjb/B1qPhI6T\n4PxeOLcHIueCwdQ60dVX65fbpCcE99Aeu11KoGvbuC3vDXyPxzY/xkvbXuLDwR9eXbJaYXaUsFcB\nndDxRv83mLR8Eq/vfJ2vhn1V5SqKe1P24mHvSqt9C6HPw+Df/rLjv0TE8feJVP4zoSMtLRTSWB20\nKJlOpvK++1k4sz92Zcv7drkV9vwIx1drkTIm7HR23Nz8ZhYcX0BWUdal+H1FlSjUF7L57GaWnVrG\nzvM7kUh6+Pfgvk73MbLZSNwd3SkxlJBakEpKfgpJ+Umk5KVwPjeJhJwkEnOTOJWxh8yiNIxc2gPx\n1DWjr/OLbD5gx16Xs3i7OtLI1YFGbo54uzjg7eqoLS6MRlj2qPbhPf5z8AiEblqyEvpiSD5kEvpo\nTexjNkDph4x3U03km/WHbtMZFDqIZ8Of5f2o9/k8+nOe6PGE5X+g1xlK2KtIqGcoT/R4gnd3v8vy\nU8uZ0GpCla6LTo6ma4kRO9fGWtu5MpxMyeWtVUcZ1MaPO/s2qwuzzYK/hzOvj+/IE/P38e22WB4e\n1PLSwdC+4NVUc8eUEXaAsWFj+fXor2yM28jkNpMtbLXtIaVkX+o+lp1cxroz68gtySXILYiHujzE\nuLBxOIsAzqTlsfpABslZiaTlFpGWW0xqbhFpuc6k5QSQXegDlI1BNyLs8nFzzcPDrYji/GYsOJZN\nsSGzQjvcHO14wHEDT+l38N9Gz3JkVSLtg/LoEuJF5yZeeDg7mlbpPaCX6aKiHEjcr4W/lq7sjyyF\nbR/C4Be5q9udnM4+zbcHv6W5Z/Mq//0oaoYS9mpwR/s72BC3gXd3v0u/4H74u/pXen5mYSansk4x\nNiMThr4BLpfKAZQYjDy1YB+ujna8P6VLteqoW4PxXYNZdSCRDzecYFh7f1r5m3rB6nRaJur2TyA3\nBdwv/Uw6+HSgmWczVp9erYS9EhJzEy+6Ws7mnMVR50wb9/74udxAQXZzVmwt5MvFRykoOXTZdZ7O\n9vh6OOHr7kT7QE98Wzni6+508TVfd9NzdydcHC/t20gpKSgxkJFfQkZeMZn5JWQWFJORX0JmXjFk\nnuGhQ79ywKU3a+2HkHImg2X7tG5iQkArP3e6hnrTNcSLrqHetAv0xNHJQ0taK5twF79L88evfBIR\n8TUvDXuV+MB4Zu+cTYhHCD0Delrk53s9IqxR4yE8PFxGRdlmdlpcdhyTl0+mb1BfPh/6eaWCvCV2\nDY9v+xc/FrnT84F/QHfpj+vD9cf5fPNJ5t7Zg5s7BVnC9FqTklPIyI//xsPZnvcmd6VfSx/TgWPw\nVR+4+V3oO/Oya77e9zVf7/+aDVM2EOBmvoqUDYHv9qxj3rEfSS05DEJizA+jKLMH+pzOYHTCyV5H\nMx9XmjZ2o7mPK8183WjW2JVmPq4EejnjZF8Hm+xGI/w8Xlt9/18EeGnJShl5xexPyORAQhb74zPZ\nn5BJWq7mZ3e009E+2JNuJqHvEuJNmK+b1pFLSji2Cja+BhdOktWsH3e6G8g0FPLb6N8I9bw6l0NR\nMUKIPVLK8Guep4S9+vx0+Cc+iPqAOQPmlN/mzsRHi6bwa84xdg76GqcWN158fU9cBlPn7uCW7iF8\nOK2rJUw2G3viMnhqwT7OpudzZ9+mvDCqvRaaOfdG0NnDQ1suOz8uO46xS8bybPiz3NPxHitZXf94\nZdN3LIn/FPTeeBv608ZtCO18m9Hcx42mPq4093HD38Pp6naFdc3u72DV0zDuM+hZ8e9LSsn5rEJN\n5OMz2RefyaFzWeQVayGRHk729GrRmOdvbkfbQA8wlMDen2Hr25wtyuCO0FAauwfx67gFqgF6NVDC\nXocYjAbuWXsPp7NOs2ziMnxdyimjm36aOxfejHBpzC937bz4cl6RntGfbUNvkKx98kY8nGufFGJp\n8ov1fLj+BN9vP02wlwtvT+rMwLT5sP5leHQP+La67PzbV96OQRr4Y9zVWarXG1JKnlz7EZtTfsTN\n0JHl077B372eCFvmWfiqH4T00tohVtM9aDBKTqXmsi8+kwMJmaw+mEROYQmPD23NzMEtcbDTab74\nHV+wO+orHvLzopejH1+Om4+Dh/o2VxWqKuwq7qgG2OnseOOGNyjUF/Kfnf8pN3a4cN1LHHZ0oHvL\nUZe9/p+VRzibns/Ht3azSVEHcHW055WxHVg4sz/ODjru/n4Xb8Z1QCIqLDFwNP0osZmxVrC2/mCU\nRu5b/gqbU36kkbEP62//sWaiHvsXZJwxr3FSwvLHtcfjP6u2qAPY6QRtAjyYFh7KmxM7s+GpgdzU\nMZAPN5xgwhfbOXw+y9RU5kV6zdzNa+4d2FmSxju/3Ij8+0MoqXmcu5RSlQ4ugxL2GhLmFcYj3R9h\nc/xm1p65os/nqc0cituEXgh6hg66+PKGI8nM3x3PwwNb0rtFxY01bIWezRqx6vEbeWRIS344WMQu\n0Zn8qN+uSlq5ufnN6ISOlbErrWSp9SkxljB14RPsyVxGIMNYd+fXeLo4V2+Q7ESYP13zgf93EJzZ\nbj4Do3+B2C0w4nUtXNEM+Lg78cUdPZh7Z09ScoqY8MV2Plp/nGK9ETwCmTjlD+4Pm8gfbk7Mi/oY\nPu8J0fOqnOGaX5LP5rObeXX7qwz+YzBD/xzKZ3s/IyEnwSz22zLKFVML9EY9d6+5m/iceJZOWIqP\ni49WN+Xr/nxjX8jnLpJ/bvsHLycvUnOKuPmTv/H3dGbpI/3rZuPLihw6l8W6Xz/kmYJP+ajpl9x7\n6zQau12qdfPYpsfYn7qf9VPW42xfTUGzcfJL8pm8aBYJRXtpYTeZRbe9gkN1fv9GI+z9CTa8qiUG\n3fCkKfPzNEz4CrpMrZ2BWee0zNKgrnD3ci3Sycxk5hfzxoojLI4+R9sAD96f2oUuId4YpZGntz7N\nlrOb+bzEg4EJhyC0D9y7CsqpXZOan8pfCX+xNX4rEYkRFBmK8HDwYECTAeTr89l2bhtGaaRfUD+m\ntJnCkNAhZqmBYzZSjoF/uxpfrnzsFuJU5immrpjK4NDBfDT4IzixDn6bxsyuQ0nWwZIJS5BS8uDP\nUfwdk8bKxwbQJsDD2mbXCcV5meg+bMPvJYP4xPEhXp/QkTGdgxBCsDtpNzPWzeCVvq8wre20aw/W\nQMgqzGLCohmklcTQ2flefp321OUJXtciLQZWPAFx26H5jTDuU/BpCQUZsOAurbXikJdh4LM1cp8g\nJcybqo0/awc0blH9MarBpqPJvLTkIGm5xTw0MIwnhrXGSBH3rr2XuOw4fmk2mTYb34IRb8ANTyCl\n5GTmSbbEb2Fr/FYOph0EoIl7E4aEDmFw6GB6BPTAQaeJd1JeEktOLmFxzGKS8pJo7NyYCa0mMLn1\nZJp5WjlXZO8vsPwxmPYTdKhZHL8Sdgvy7cFv+XTvp3ww6ANuivgZw9kIBoT4MarFaF7t9yq/7zrL\ni4sP8srYDtw/oG7/cKzOn/ehP7WVqW4/EH0uj5s6BvCfiZ3wc3fitlW3kV+Sz7KJy6qeVm40XBYm\nakucz01iypL7yNYn0s/zMb6ZNKPq+QqGEtj+Kfz1Hjg4a4Xjut95uXjrizShOLBAOzb2k3JXuZWy\n7zdYOqvcUNW6IqughDmrjrIgKp5W/u68N6ULTXyKuWPVHdjp7Pi52JO487vY2udutqREcS73HACd\nfTszOHQwg0MH09q7daU/S4PRwI7zO1h4YiF/JfyFQRroHdibKW2mMKzpMBztLFw5NXoeLHsEWg6B\n237Xfqc1QAm7BdEb9UxfPZ2k3PMsOXmM1G63MuXCX8wZMIfOXkMZ/dk2ujf15pcZfSwfvmZpjq+B\n32/DcNt8/pfcho82nMDFwY5Xx3bAudF+Xtj2Ap8P/bzctmtXEb8bfp0EjZpBl9ug81SwkeiJUxmx\n3Lb8fgoM2Yz0fYEPx02puqif26NtZCYfgg4TYdR7Fb9vKWHr2/DXuxA2GKb9DM5VLN+QnajlH/h3\ngHtX14kLpjL+PpHKi4sPcj6rgPtvaMGYcCMPb5pBob4QicQJQd+QgQwOHcygkEH4ufrVaJ7U/FSW\nnlzKophFnMs9h7eTN+Nbjmdym8mEeYWZ+V2Vw77fYOn/ab+f23+/qrJrdVDCbmFOZJzg1uVTGZGb\nQ/fejzPnyHesmriGJ+bFcSoll3VPDSTIy3KdkKyGvhg+bKutTKZ8z6nUXP618AB74jKYNbgZG3Of\nIsQ9hB9u/qHycZIOwY+jteqRrj5amrrQQcuh0PV2aDsaHF0t856qyb7kg9y39mGK9UamhrzO7Jtu\nqtqFxXmw+S2I/BrcA2DMh9BuTNWujZ4HKx4H3zZwxx/lNnG5DClh/h1warPmgvFpWfn5dUROYQnv\nrDnGvMizNPdx5e6hxSQURzIgM42+u37G9Y6F0Hq4WeYySiMRiREsPLGQLWe3oJd6evj3oINPBwAk\nWmRNaWG1Kx9fPAeJDh29gnoxLHRY5T78/fNhyUwIGwS3z6+VqIOFhV0IcTPwKWAHfCulfKey8xui\nsCMlc//Xgy+d9AS7BWOQBsY1+oqPN8bw6W3dmNCt6u3GbJ5Vz2hC81wMOHlgMEqeWrCPtYeTeGzi\nOb45/Cnzx8yvuK/mhVPw/c2aW2HGWi1KI/UEHJgP+xdAdgI4emh+yq63QrMBl602z2akMmvtbDxd\n7GjeuBGNXdxxsXe5/Obggqu961Wvezh64OFY8z2QbfE7eXTTY+hLXJjR8m2eGVp+T9urOLkRVj6l\nxZKHz4Dhs6u+8i4ldqvmd3dwhTsWQHC3is898CcsfgBGvllu0xdLs+NkGs8vPkBCRgHTeoYyuasv\nvdaORxj1WgZsDV0XFZFWkMbyU8tZenIpqfmpAAgE2v/i4rcrUfpfmW9bAkGxoZickhwaOzfmlla3\nMLnNZEI9rvgw3b8AljwMLQZqom6GhYjFhF0IYQecAEYACcBu4HYp5ZGKrmmQwn42kpLvR3JH+3CO\nFabQ138Ym7eNZEznID67vbu1rbMsZyPh+5Ewce7FioDx6fkM/XArE3v68E/hk9zY5EbeG/Te1ddm\nJWiiXpIP960FvytKyxqN2kbfgflweBkU54BniFaArOtt4NeWifOf5mThJmSJN+iKsbMrQYqiKpvv\n4+xDS++WtPBqQZhXGC29WxLmFYavi2+l7pQVJ9fy739eRF/sw5Md3+OhG3pce7L8dFj7ovZ+fFpr\nMeTN+lfZ1qtIPgK/TdPGnfqD1uXqSnJT4Mve4NMKZqyrN3sYeUV63l93nN92naVYb2SU63G+Nr5O\nXOfHCZ74upbgVE8wGA3sTNzJH8f/4O+EvzFIA/2D+zO1zVQGhQ7C4dASTdSbD4DbF5jt26Ulhb0f\nMFtKeZPp+YsAUsq3K7qmQQr7skfg8FKOz1jJ9I0P4pg5BZnTi7VPDMTLtR6FW1kCKeHTrtA4DO5e\nevHlV5Ye4vddZ7lj1H6Wn17AmklrCHIvUycnNxV+GAW5yXDPispXnADF+Vq54AML4OQmkAaOBHTi\nNpdsWtoPYs6wOWw+lsLGo8nsi89AUoK/l6BfK0/CW7jRPtgJI8UU6Asu3jKKMjiddZrYrFhiM2PJ\nLbnUFtDD0eMyoQ/zCiPMO4wgtyB+ObSAD/a8jaGwKa+Ef8Dt4ZWEtBn0kHQAzvyjbZAWZsKAp+DG\nZ82zMs1J0sQ96SCMfh96PXD58QV3adFbM7eBX9vaz2dmcov0bDmWwtpDSYw+8W+Gs5vJuo9o16Er\nozoFMqC1b52FC+cV6UnNKSI1t0i7zykiJafw4uPS1zPySujUxJOxXYLp09qOv5NWsejEIpLzk/Gz\nd2di2nmmeLUn+I7FZnUZWlLYpwA3SykfMD2/C+gjpXz0ivMeAh4CaNq0ac+4uLhazVuvKMqBD9pC\np0kw4QteXLqL3yNSmPdAX25oVU65geuBzW/Btg/g6aNaLW+0TlED39vC0E4O7Cx+luntp/Ncr+e0\n8wsy4aexkHYS7lpc/VVrbgocXMi9B+ZywsHAqvMpNOr1EAx6Hpw9ScstYuvxVDYdTebvE6nkFRtw\nstdxQytfhrbzZ1h7/6v2QKSUpBakcirz1EWhj83SbumF6RfPc9Q5U2wsxJDbjndufJ/xXZpfbltR\nLpyLgrMRcHantilckqcdC+kN4z6BgArcUjWlKBcWzoCYdZqrZfgbmrvq8BL4814Y9prWUL2eU5ie\ngP1XvYlx6si0vGfIKTTg7mTP0Hb+jOoUyKC2frg6XrtIrdEoSc0t4nxmAUlZhZzPKiQxs4DErMLL\nhLu01k1ZdAJ83Z3w83DC30O793B2YPvJNI4l5SAE9GremDFdAvAvnMfKYz+xzcUFhGBAkwFMbTOV\nG0NuxF5X+2K6lhT2qcBNVwh7byllhY67Brdi3/OTtnF1/wa+j/PjjZVHuH9AC14Z28HallmPtBj4\nIhxumgP9Hrn48pzVR/nftlhGDdnE3rQdbJiyAQ/s4JdbtGiQ2+fXeLPs533reH//swx2HsvnLrla\n3LCbn5ZN2eW2i374Ir2B3acz2Hg0mU3HkolP11LZOwR50rtFY/RGI7mFenKL9OSY7nOL9OQW6skp\n0muZk3Z52DmmonNKRueYgg5Xvhj1LMPaB2vfPM7uvCTkiftBGgABgZ2gaT9o2le79wyu9Y+6Qgx6\nWPsC7P4ftB8PN7+tZax6h8L9G8HORqp27/wK1r1IyZSf+cehH2sPJrH+SBIZ+SU4O+gY3MafUZ0D\naeHrRmIZwT6fVUhSVgHnMwtJzi5Ef0XfXid7HcHeLgR4OuHn4YzfFeJdemvk6lhh7sHJlFxWHjjP\nygOJtEvbwKcOX3DCqRPbb3ibTJe9rD6zjNSCVPxd/ZncejKTWk+qVge2K1GuGEvy7XAozObHbvOZ\nvfIoN3cM5PM7utcrn6BV+GYwSCM8fKnvaXpeMQPf20L3VrnsM87m2R5Pcs++VVo6+5TvoeM1enBW\nQImhhH4/j6XIUMxft6+isZur1vRhzb8gYbdW2GrUe1pziDJIqRWu2ng0hc1HUzhwLhNXR3vcnUw3\nZ3s8TPeXPXeyx93ZAXcnexrLTNrk7MQ7dY8m5BdOaoPbOUFIuEnE+0Nor+pviNYWKbW+pev+rUVk\nGEq030eADS06DHr4ZpD2re6RSHByR28wsut0OmsOJbHucBIpOZfvoTja6Qj0cibIy5lgbxeCvJwJ\n8nYhyNOZIG9ngr1c8HZ1MFsfBHloCSy6nwT3zjxoeJ5j6RJ7neCG1t60ah5PXPEmdiVFIITg48Ef\nM7Tp0BrNY0lht0fbPB0GnEPbPL1DSnm4omsalLCbapFHtX2aKfvDGdkhgC+n91CiDhDxtbZifGTX\nZb7cjzac4LNNMfTu9zsZWUdYHXsKh/FfQI+7ajzVG399y59nPmVC0Eu8OfL2SweMRm1jcsNrkJeq\nzTHstcv6ctYIo0ELFdzzI5xYC0a9FppZdjUe3E1rLVcfOLJc28wb9C/Nn29rlG7I3/CElpVaBqNR\nEh2fQVpuMcFeLgR5O+Pj5mi55jWHl2pur9DeMH0h0tGNQ+eyL67kz2UW4Gino3cb8Pbfy78HzCLI\no2a1oiwd7jga+AQt3PF7KeVblZ3foIR93b8xRsylV8HndG/fhq+m96h3DamtRk4yfNQOBjwNw165\n+HJ2YQkD39nEg40+Za53Au8GDmP0TZ/UeJrMwmwG/X4TOn0gkfctxLG8jbXCLC2LM3IuOLjBkJe0\nTcXquiMyz0L0r9ot+5wWY9/1di0ix7+jxZN8qoW+GOwtnHFpTpY9osWFz9xeq3orZuXIMvjzPu0b\n4Z0LteqVZZBSEh2fycr9iaw6eJ7k7CK+uasnIzvWzB2jEpQsgb6YwvfasKWgFX+2fJuv7+zR4Ip7\n1ZpfJsGFGHjiwKV0eCnZ/93/0TnhN8a07ICHdygLxi6o8Qrr4ZVvsD1tIU+0/4IH+wys/OTU47Dm\nec31499Bc8+UaYJSLvpiLfpm78/aKh20RKked2uJUrYslrZEXppWATKgE9y7sma1cczJkeWw8D5o\n0hPuXHSVqF+J0SjZfSadrqHeODvUTCdUPXYLsH31LzgXZ3AsaIIS9Yrocqu2yo2PvPTaX+/RNeE3\nFojRGIomcDT9KFHJNfugP5l+lh1pi/HU9+GB3tcQaNBcQnctgVt/heJcLRLnz3u1+PkrST2hNQ/5\nqD38eQ+kHtNcGU/s1yJ3Ok5Uom5J3Hy1xK24f+Dgn9azo6QADi7URD24O0y/eqVeHjqdoE+YT41F\nvTrYyLZ4/ePPqHh8d/9EuqMvs+5/WIl6RbQbo2VCHlig+Z4jvoatc6DrHRT5PcuJVQcJ6uTNT4d/\noldgr2uPdwXPbpqDlDpm3/hs1Vf8QkD7cdBquBZH/s/HcHwt3PgM9Lpfi/He+5O2EaqzhzY3Q497\noNWwepPMc93S4x6tdvy6f0PrkeDibf45CrMgM15bkGSVvY/X7vO0TFWahGsrded60gGrDErYa8Ci\nPQl8tGgL250OYuz7FPZOatVWIU7umrgfXgKBnbXN1HZjYfzn3C4F32w7g13+DfyVsIrYrNhqFWX6\nK24Xp/K3E6qbwMi2NUi0cXCBwS9oPvL1/4Ytb2o3gMYtYfjr2jEbKTx2XaDTaTV0/jcUtsyB0eVk\nL1eVkkI4tAgS910S7cx4KMq6/Dx7Z/AK1cJEAztr997Noe0o7d93PUQJezVZGn2OZxfu5z2/KHTZ\nRnQ9ax7Jcd3QeZr21XnlUxCmFQfDzh4n4PFhrXlxWRrebdfzy5FfeK3fa1Ua0iiNvLJtDsYST94f\n/UTt7GvUTHPNnNqi1WxpOwqa3WB9H66ifIK7axvfu/8H3e64dobylRRkaE27I/8LeSng5KnVI/IK\n1RLjSkXcq6l27+Znc/8WlLBXg2X7zvH0H/vo17wRkwu2asV96rgxQYOg5RDtj8QzCG6bd1kI4OSe\nIcz9y5+Cwt6sOLWCR7s9qnWiugY/H1hKhuEUPTxm0im4ZuVcy7Wz5RDzjKWoW4b8WwszXPUM3L+h\natFImfFaTP+en7TM31YjtPDJ5gNsTrivhdo8rSIr9p/nqQX76NW8Md8PLkCXGQfd77a2WbaBnQPM\n2g73rQFHt8sOOdjpeGpEG1IT+lBkKGLB8QXXHK5QX8gX+z5FFjXh/VH31ZXVivqMi7dWmfJcFET/\nXPm5SYdg8UNa/aJd32j7KzO3a+GJLW5scKIOStirxKoDiTy5YB/hzRrzw329cD74m5ZB2H6stU2z\nHZw9K9x4HNclmNaNWuJY1Infj82nUF9Y6VDv7vgvRaQzOvhhAj3rZ012hQXoMk0r2bxxNuRduPyY\nlADBS+YAAAs7SURBVBD7lxZuO/cGOLYK+s7SIpom/Vcr7dCAUcJ+DdYeSuTx+dF0D/Xmh/t64WrI\ngaMrNL9xLYvmKzR0OsHTI9uQmdSPzKIMVsSuqPDc1PxUFsX+gi6/M6+OGGdBKxX1DiFgzAdaEb6N\npr0Zg17bEP1mMPw8XqtwOexVeOoQ3PQWeIVY1WRL0eB97AXFBr75O5bD57NwcrDD0U6Ho70OJ3vt\nvvR56WMnh0uvXcgtZs7qo3QN8eLHGb1xc7KHyD/BUFSr9HfF1YzsEEDHLd2JKwnhp8M/M7n15HL7\nor645X2M6Hmo46O4OzX4f76Ka+HfHvr+H+z4TItzP7QYMuO0WvPjPtPyKMzcpMMWaNB/GZuOJvPa\n8sMkZBTQyt8do1FSpDdSbDBSVGKg2GCkWG/EWEnybbdQb36a0fuSiET/DEFdtZvCbAgheHZkO2Ys\nHECcw3y2JWxjUOigy845nHaUyNS1uBYOZtYNfa1kqaLeMeh5bZX+z8daGeSb5mgZwfW5vEMd0yCF\nPSEjnzdWHGH9kWRa+7sz/6G+9A2rONJCbzBeFPlivZEi001vNNLSz/1SQa/z+0zNCz6w0Du5vrix\ntS/dfQZyTL+WHw79eJmwSyl5fstbSIMLL/Z7TBVZU1zCyV1rzFKQCSE9rW1NvaBBCXux3sh3/5zm\ns00xALwwqh33D2hxTRGwt9Nhb6fD9Vp5RtG/aMkKnaeayWJFWYQQPDeyI9P/7M8e+9UcvnCYjj5a\nA4oNZ7YQl7+fAOOtTOzSysqWKuodVmrGXV9pMMuenacuMPqzbby79hgD2/iy8ZlBzBzU0nwru5IC\nrQFw+/F1k8asAKB3i8b09h0FRie+O/AjACXGEt7Y8R6GIj/mDH/YcuVYFQobxeZX7Kk5RcxZfZQl\n0ecIbezC9/eGM7RdHaSAH12hpRqrTdM65/mR3Zj6Ry826taTmPs0y0+uJ0t/js4uT9KnhZmSkRSK\nBozNCrvBKPktMo731h2nsMTA/7d390FW1XUcx99fV5aVpwB3MQJ041EImIXWeNgsIgUkDJmwsqbQ\npmSQSCpnMp9qmtF/yAJHxiRyBEXNyUByrAGpCQKiXZYFFiFBXAFZeUweJILd/fbHudQO+wR7H86e\ncz+vmZ1777l793zmN2c/c+Z37/2dOeP7c8+4/lyVm6ZFmsqXQrfC4HOzklbDe3dlTP5tlNVtYP7m\nhax+dw21H/bnselfDjuaSCREsti37v+Ah1ZUsv29E5T0v5qfTR1Kv4I0LsZzfC9UrYPxD2X1O+2Z\n9MDEsXzx5WG8XvUq7sbNH/02/Xq0vDSqiESs2E+cOc+8VbtYtmkf+Z3a88QdI7h1eM/0z7lueR7s\nCij6enr3I/8z8JrOjL16GqXnt8KpG3h4eusucC2SjSJV7I+srOQPWw9y59hCvn/zQLrktUv/Tmtr\noOKFYO3udF5RXhp4ZMIkvrDofWaNGUd+pzZy7VCRCIhUsd83YRDfubEvQ3tl8Ervb6+BU9UweV7m\n9ikAFOZ3pPS+mXTIjdRhKhK6SP3H9OnegT6Z3mn50mA95oGTMr1nAZW6SCvoncDmnD4Mb/0puAJ9\nTgamfUREUkDF3pytL0FdjdZdF5FIUbE3pa4umIbpMwoKBoadRkTkkqnYm7J7FRzbHVxbUUQkQlTs\nTVm/ILio7SemhZ1EROSyqNgbs78U9m2AMbP1pqmIRI6KvTEbFkBeVxihBb9EJHpU7Bc7ugd2vhbM\nrbdP4/ozIiJpomK/2MYnIScXRs0MO4mISKuo2Os7fThYF6boa9CpR9hpRERaRcVe36anofYcjJ0T\ndhIRkVZTsV/wn9NQuhgGT9H1E0Uk0lTsF2x5Ds5+ACVzw04iIpIUFTtA7XnYuBCuK4HexWGnERFJ\nioodYMdyOLEfxn4v7CQiIklTsbsHywcUXA8DJoSdRkQkaUkVu5nNM7NdZrbNzJabWddUBcuYt/8M\nhyqDs3VdqFpEYiDZJlsNDHX34cBbwI+Tj5Rh6xdA554w7Pawk4iIpERSxe7uq9y9JvHw70Dv5CNl\n0MEt8M5fYfQsuDI37DQiIimRyrmHbwF/bOpJM7vbzMrMrOzIkSMp3G0S1j8B7bvAJ+8MO4mISMq0\nWOxm9oaZVTbyM7Xe7zwI1ADLmvo77r7I3YvdvbigoCA16ZPxryp4cwUU3wV5Hwk7jYhIyrR4CXh3\nv6m5581sBjAF+Ly7e6qCpd3GhWA5MGpW2ElERFKqxWJvjplNAn4EfNbdz6QmUgZ8eAzKn4PhX4Eu\nPcNOIyKSUsnOsT8JdAZWm1mFmf0qBZnSr3Qx1Pxbi32JSCwldcbu7v1TFSRjzp2BfzwNA2+BHteH\nnUZEJOWy7xs5FcvgzDEo0fIBIhJP0Sr2o7vh0Jutf31dbXCFpN43wLVjUpdLRKQNiVaxr/05PDUG\nFt8E5UuDNdQvx86VwcccS+4Fs7REFBEJW7SKfeKjMPExOHsSVs6BxwcFtwfKgsW8muMOf5sP3fvB\noMmZySsiEoKk3jzNuI75MGY2jL4HDpRC+RLY/rvg7L3HEBj5zeAjjB26N3xt1TqoroAp8+GKnMxn\nFxHJEAvjO0XFxcVeVlaWmj929iRUvhKU+8FyyMmFwbfCyBlQeOP/V2x8/ktQvRXmVkK7vNTsW0Qk\ng8xss7u3eDWgaJ2xNyavS7AsQPFd8P724ItH234blH23QhjxDfhYEex5A8Y/rFIXkdiL/hl7Y86f\nhV2vBVM176wNtrXrCD/YAVd1S99+RUTSKHvO2BvTLg+GTQ9+ju+Fihchf4BKXUSyQjyLvb7ufWH8\ng2GnEBHJmGh93FFERFqkYhcRiRkVu4hIzKjYRURiRsUuIhIzKnYRkZhRsYuIxIyKXUQkZkJZUsDM\njgDvtvLl+cDRFMaJA41J4zQuDWlMGorSmFzn7gUt/VIoxZ4MMyu7lLUSsonGpHEal4Y0Jg3FcUw0\nFSMiEjMqdhGRmIlisS8KO0AbpDFpnMalIY1JQ7Ebk8jNsYuISPOieMYuIiLNULGLiMRMpIrdzCaZ\n2T/NbI+Z3R92nrbAzKrMbLuZVZhZGq832HaZ2TNmdtjMKutt625mq81sd+I2qy6f1cSY/NTM3ksc\nKxVmNjnMjJlmZn3M7C9mttPMdpjZvYntsTtWIlPsZpYDLARuAYYAd5jZkHBTtRmfc/eiuH0W9zI8\nC0y6aNv9wBp3HwCsSTzOJs/ScEwAfpk4Vorc/fUMZwpbDfBDdx8MjAZmJzokdsdKZIod+BSwx933\nuvs54CVgasiZpA1w97XA8Ys2TwWWJO4vAW7LaKiQNTEmWc3dq929PHH/FLAT6EUMj5UoFXsvYH+9\nxwcS27KdA6vMbLOZ3R12mDbkGnevhuAfGugRcp624rtmti0xVRP5KYfWMrNCYASwiRgeK1Eqdmtk\nmz6rCSXuPpJgimq2mX0m7EDSZj0F9AOKgGrg8XDjhMPMOgGvAHPd/WTYedIhSsV+AOhT73Fv4GBI\nWdoMdz+YuD0MLCeYshI4ZGY9ARK3h0POEzp3P+Tute5eB/yaLDxWzKwdQakvc/ffJzbH7liJUrGX\nAgPM7ONmlgt8FVgZcqZQmVlHM+t84T4wAahs/lVZYyUwI3F/BvBqiFnahAvllTCNLDtWzMyA3wA7\n3f0X9Z6K3bESqW+eJj6eNR/IAZ5x90dDjhQqM+tLcJYOcCXwQjaOiZm9CIwjWH71EPATYAXwMnAt\nsA+43d2z5s3EJsZkHME0jANVwMwLc8vZwMw+DawDtgN1ic0PEMyzx+pYiVSxi4hIy6I0FSMiIpdA\nxS4iEjMqdhGRmFGxi4jEjIpdRCRmVOwiIjGjYhcRiZn/An1sPrqvpGCNAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f6390c40fd0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from numpy.matlib import repmat\n",
    "\n",
    "%matplotlib inline\n",
    "\n",
    "# Simple time-series (from CrossValidated so that you can compare the results)\n",
    "x1 = [-1.042061,-0.669056,-0.685977,-0.067925,0.808380,1.385235,1.455245,0.540762 ,0.139570,-1.038133,0.080121,-0.102159,-0.068675,0.515445,0.600459,0.655325,0.610604,0.482337,0.079108,-0.118951,-0.050178,0.007500,-0.200622]\n",
    "x2 = [-2.368030,-2.607095,-1.277660,0.301499,1.346982,1.885968,1.765950,1.242890,-0.464786,0.186658,-0.036450,-0.396513,-0.157115,-0.012962,0.378752,-0.151658,0.774253,0.646541,0.311877,-0.694177,-0.412918,-0.338630,0.276635]\n",
    "x3 = np.array(x1)+np.array(x2)**2 # ad hoc\n",
    "\n",
    "# Plot\n",
    "cdata = np.array([x1,x2,x3]).T\n",
    "plt.plot(cdata)\n",
    "plt.title('Sample time series')\n",
    "plt.legend(['$x_1$','$x_2$','$x_3$'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(23, 3)"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAAEICAYAAABLdt/UAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xd4FNX+x/H3Se+QnpAQEnqXXqRIEaWJ2AAVxIr8xGvB\n3rDrVa8oil4reEUQEEXpAqEqNfTeAgmBENJI7zm/P2ZBREpCNpns5vt6njybnczOfHez+WT2zJlz\nlNYaIYQQ9sPB7AKEEEJYlwS7EELYGQl2IYSwMxLsQghhZyTYhRDCzkiwCyGEnZFgFzWOUmqVUurB\nq3xshFIqWynlaO26hLAWCXZx1ZRSdymlYixBl6iUWqyU6m52XdaklDqmlLr+7H2tdbzW2ktrXWJm\nXUJcjgS7uCpKqfHAx8A7QDAQAXwO3GxmXUIICXZxFZRStYA3gHFa61+01jla6yKt9Xyt9TOWdb5T\nSr113mN6KaUSzrt/TCn1jFJqp1IqRyn1rVIq2HLUn6WUWq6U8rWsG6mU0kqp+5RSx5VS6UqpsUqp\njpbHn1FKTT5v268ppX447/7Zxztd5Lk0UEqtUEqlKqVSlFLTlVK1LT+bhvEPa77lU8mz529LKTVC\nKRVzwfaeVErNs3zvqpT6j1IqXimVpJT6QinlfonXtKFSarVSKsNSx6xL1X5+U5JS6l6l1J9KqY8s\nr0OsUupay/LjSqnTSqnRZf/tCnsgwS6uRlfADZhbwe3cBvQDGgM3AYuBF4EAjPfmYxes3xloBAzH\n+LTwEnA90AIYppS67ipqUMC7QB2gGVAXeA1Aaz0KiAdusjS/vH/BY+cBTZRSjc5bdhcww/L9e5bn\n1gZoCIQBEy5Rx5vAUsAXCAc+Lcdz6AzsBPwt+54JdLTscyQwWSnlVY7tCRsnwS6uhj+QorUuruB2\nPtVaJ2mtTwBrgY1a621a6wKMfxptL1j/Ta11vtZ6KZAD/Ki1Pn3e4y9c/4q01oe11su01gVa62Rg\nIlCmfxBa61zgN+BOAEvANwXmKaUU8BDwpNY6TWudhdFsNeISmysC6gF1LM/xj3I8jaNa66mWdv9Z\nGP+c3rA8p6VAIUbIixpCgl1cjVQg4GJNG+WUdN73eRe5f+FRZnnXvyKlVJBSaqZS6oRSKhP4AeMT\nQ1nNwBLsGEfrv1oCPxDwALZYmkjOAEssyy/mWYxPD5uUUnuUUveXo4YLXwe01hV+bYTtkmAXV2M9\nkA8Mvcw6ORjBdlZIpVZ09ft+F9BAa621D0bThTrv51ca/nQpxj+5NhgBf7YZJgUjUFtorWtbvmpp\nrS8asFrrU1rrh7TWdYCHgc+VUg0tz4VyPB8hJNhF+WmtMzDaij9TSg1VSnkopZyVUgOUUmfbobcD\nA5VSfkqpEOCJKixxO9DT0ue8FvDCZdb1BrKBM0qpMOCZC36eBNS/1IMtzVFzgA8AP2CZZXkp8DXw\nkVIqCEApFaaUuvFi21FK3aGUCrfcTcf4h1JiaR46AYxUSjlajuQbXOb5CCHBLq6O1noiMB54GUgG\njgOPAr9aVpkG7ACOYRzVzqrC2pZZ9rcT2AIsuMzqrwPtgAxgIfDLBT9/F3jZ0pzy9CW2MQPjJO5P\nF5x3eA44DGywNPMsB5pcYhsdgY1KqWyMk7KPa62PWn72EMY/nFSME8XrLvN8hEDJRBtCCGFf5Ihd\nCCHsjAS7EELYGQl2IYSwMxLsQghhZyp6gclVCQgI0JGRkWbsWgghbNaWLVtStNaXusjtHFOCPTIy\nkpiYmCuvKIQQ4hylVFxZ1pOmGCGEsDMS7EIIYWck2IUQws6Y0sYuhBCVraioiISEBPLz880updzc\n3NwIDw/H2dn5qh4vwS6EsEsJCQl4e3sTGRmJMTy+bdBak5qaSkJCAlFRUVe1DWmKEULYpfz8fPz9\n/W0q1AGUUvj7+1fok4YEuxDCbtlaqJ9V0bol2EWZlepSVsSvYE3CGmRUUCGqL2ljF1ektWbF8RV8\ntv0zDqUfAqB7WHde6vwS4d7hV3i0EKKqlfmIXSk1RSl1Wim1+7xlr1nmitxu+RpYOWUKM2itWZOw\nhhELR/DEyicoKinivR7v8VzH59iatJVbfruFKbunUFRaZHapQojzlKcp5jug/0WWf6S1bmP5WmSd\nsoSZtNasP7meUYtHMS56HBkFGbzV7S3m3jyXgfUHMrL5SH4b+htd63Tloy0fMWLBCHYk7zC7bCGq\npd69e7Ns2TIAXn75ZR577LFK32eZm2K01muUUpGVV4qoDrYkbWHytsnEJMUQ7BHMhK4TGNpwKM4O\nf+9PG+IZwid9PiE6Ppp3N77LqEWjGNZkGI+3exxvF2+Tqhfi4l6fv4e9JzOtus3mdXx49aYWV973\n668zYcIETp8+zbZt25g3b55V67gYa7SxP6qUugeIAZ7SWqdfbCWl1BhgDEBERIQVdiusaWfyTiZv\nm8z6xPUEuAfwfKfnub3x7bg6ul72cX0j+tIltAuTt01mxv4ZRMdH83yn57mh3g022yNBCGvq2bMn\nWmsmTpzIqlWrcHR0JDY2lrfffpuMjAzmzJlj9X2Wa85TyxH7Aq11S8v9YCAFY0b1N4FQrfX9V9pO\nhw4dtIzuWD3sS93HZ9s/Y3XCanxdfXmg1QMMazIMdyf3cm9rT+oeXl/3OvvS9tEjrAcvdXmJMK+w\nSqhaiCvbt28fzZo1M7sMdu3axW233UZAQADr1v19HvLbb7/9ksF+sfqVUlu01h2utM8KdXfUWidp\nrUu01qXA10CnimxPVJ3E7ETGrxrPsAXD2Hp6K4+1fYzFty1mdIvRVxXqAC38WzBj0Aye6/gcW5K2\nMPTXoUzdPVVOrooaKzExkbvvvpvffvsNT09Pfv/99yrZb4WCXSkVet7dW4Ddl1pXVC8T1k3gjxN/\nMPaasSy5bQkPtX4IT2fPCm/XycHpbydXJ26ZKCdXRY2Um5vLrbfeyocffkizZs145ZVXeO2116pk\n3+Xp7vgjsB5oopRKUEo9ALyvlNqllNoJ9AaerKQ6hRXtSd3DhsQNjL1mLOPajMPHxcfq+zh7cvXj\n3h+TUZDBqEWj+GHvD1bfjxDVlYeHB+vXr6dfv36A0da+fv16AFJTUxk7dizbtm3j3Xfftfq+y9Mr\n5s6LLP7WirWIKjJ191S8nL24o/Edlb6vsydXn1/zPB9u+ZDOoZ1p5Nuo0vcrRHXm7+/PF198UWnb\nlyEFapjjmcdZFreMYU2GVVm3RE9nT17v9jrezt5M+HMCxaXFVbJfIWoqCfYa5rs93+GoHBnZbGSV\n7tfPzY8XO7/I7tTdTNs7rUr3LURNI8Feg6TkpfDr4V8Z0mAIgR5XnOjc6m6MvJG+EX2ZvG0yRzOO\nVvn+hagpJNhrkBn7ZlBUWsS9Le41Zf9KKV7u8jJuTm5M+HMCJaUlptQhhL2TYK8hcopymHlgJtfX\nu57IWpGm1XH2qtbtydv5cf+PptUhhD2TYK8h5hycQ1ZhFve3vOKFwZVucP3B9AjrwaStkzieedzs\ncoSwOxLsNUBRSRHf7/2eTiGdaBnQ0uxyUEoxoesEnByceHX9q5TqUrNLEsKuSLDXAAtiF3A693S1\nOFo/K8QzhKc7PM3mU5uZc9D6gyAJUZNJsNu5Ul3KlN1TaOrXlGvrXGt2OX9za6Nb6RLahQ9jPuRk\n9kmzyxHCbkiw27mVx1dyLPMY97W4r9oNo6uU4rVrX0OjeX396zKPqrBL1XqiDWF7tNZM2T2FMK8w\nboi8wexyLirMK4zx7cfz9sa3+fXwr9zS6BazSxL2aPHzcGqXdbcZ0goG/PuKq5kx0YYcsduxLUlb\n2Jm8k9EtRuPkUH3/hw9rMowOwR34YPMHJOUkmV2OEFZ1/kQbM2fOxNHRkV9//ZWHHnqIm2++maVL\nl1p9n9X3r11U2JTdU/Bz82Now6Fml3JZDsqB1699ndvm3cabG97k0z6fVrtmI2HjynBkXVl27dpF\nYmIiAQEBeHsb4zMNHTqUoUOHkp6eztNPP80NN1j3E7Ucsdupg+kHWXtiLXc1veuqJ86oShE+Efyr\n7b9YnbCahUcXml2OEFZxpYk23nrrLcaNG2f1/Uqw26mpu6fi7uTOiKYjzC6lzO5udjfXBF7Dvzf9\nm5S8FLPLEaJCLjfRhtaa5557jgEDBtCuXTur71uC3Q6dzD7J4qOLub3x7dRyrWV2OWXm6ODIG93e\nIK8oj3c2vmN2OUJUyOUm2vj0009Zvnw5c+bMqZRx2aWN3Q59v/d7FIp7mt9jdinlVr9WfR5p8wgf\nb/2Y34/9zo2RN5pdkhBW99hjj1Vqt0c5YrczZ/LP8MuhXxhYfyAhniFml3NVRrcYTQv/Fryz8R3S\n89PNLkcImyPBbmd+3P8jecV51Wr4gPJycnDijW5vkFmYybubrD8fpBD2ToLdjuQW5TJj/wx6hfei\nQe0GZpdTIY19GzOm9RgWH13MrP2zzC5HCJsiwW5H5h6ey5mCM9zfynaP1s/3YKsH6RbWjbc2vsVb\nG96iqKTI7JKEsAkS7HaiqLSI/+35H22D2tI2qK3Z5ViFs4Mzn/X5jPta3sesA7N4cOmD0g1SiDKQ\nYLcTS44uITEn0abb1i/G0cGR8e3H837P99mbupcRC0awJ2WP2WUJUa1JsNsBrTVT90ylQa0G9Azv\naXY5lWJA1ACmDZyGo3LknsX3MP/IfLNLEqLakmC3A2tPrOVQ+iHua3kfDsp+f6VN/Zoyc/BM2gS1\n4cU/XuT9ze9TXFpsdllCVDv2mwI1yJTdUwj2CGZg1ECzS6l0vm6+fNHvC0Y2G8m0vdMYu2ys9HUX\n4gJy5amN25G8gy1JW3imwzM4OzqbXU6VcHZw5rlOz9HErwlvrn+TOxfeyaTek2ji18Ts0kQ19d6m\n99iftt+q22zq15TnOj13xfV69+7Niy++SL9+/Xj55ZfJzMzkk08+sWotF5Ijdhs3fd90vJ29ub3x\n7WaXUuWGNhzKd/2/o6i0iJGLRrLk6BKzSxLiH15//XXefvttpk+fzrZt2/joo48qfZ9yxG7DUvJS\nWBa3jBFNRuDh7GF2OaZoFdiKWYNnMX7VeJ5Z8wz70vbxWNvHcHRwNLs0UY2U5ci6spw/0caqVatw\ndHRk3759TJo0iZSUFPr27cv//d//WXWfcsRuw34++DPFpcUMbzLc7FJMFeAewLc3fMsdje9gyu4p\njFsxjoyCDLPLEgL4a6INV1fXcxNtNGvWjC+++ILZs2cTExNj9X1KsNuo4tJifjr4E11DuxJZK9Ls\nckzn7OjMhK4TeKXLK2xM3MhdC+8i9kys2WWJGu5yE23MmzeP7t2707dvX6vvV4LdRq06voqk3CSb\nmkijKgxrMowpN04hpyiHuxbdxZqENWaXJGqoy020ATBkyBDWrVvH9OnTrb5vaWO3UTP3zyTUM5Tr\nwq8zu5Rqp21QW2YOnsljKx7j0ehHebL9k9zb4l6ZR1VUqbMTbZx1/kQbq1at4pdffqGgoICBA63f\nTVmC3QbFnoll46mNPN7ucTlJeAkhniF81/87XvnzFSZumcih9EO8eu2ruDq6ml2aEPTq1YtevXpV\n2valKcYGzTwwE2cHZ25tdKvZpVRrHs4e/Oe6/zCuzTjmx87n/iX3k5ybbHZZQlQ6CXYbk1OUw7wj\n87gx8kb83PzMLqfaU0ox9pqxfNTrIw6dOcSIhTKImLB/Euw2ZsGRBeQU5chJ03K6vt71TBswDSfl\nxOglo1l8dLHZJQlRacoc7EqpKUqp00qp3ect81NKLVNKHbLc+lZOmQKMURxnHphJM79mtA5obXY5\nNqeJXxNmDJpBC/8WPLvmWT7Z+gmlutTsskQl0lqbXcJVqWjd5Tli/w7of8Gy54ForXUjINpyX1SS\nmKQYDp85zJ1N75QeHlfJ392fb274htsa3cbXu77miZVPkFOUY3ZZohK4ubmRmppqc+GutSY1NRU3\nN7er3kaZe8VordcopSIvWHwz0Mvy/f+AVYB51+7auZn7Z+Lj4kP/qAv/v4rycHZ05tWur9LItxEf\nbP6AkYtG8mmfTwn3Dje7NGFF4eHhJCQkkJxseyfM3dzcCA+/+vdjRbs7BmutEwG01olKqaBLraiU\nGgOMAYiIiKjgbmuepJwkVsSv4O5md+Pu5G52OTZPKcXdze6mfq36PL36ae5ceCcTe02kY0hHs0sT\nVuLs7ExUVJTZZZiiyk6eaq2/0lp30Fp3CAwMrKrd2o05h+ZQoktq/Lgw1ta1TldmDJqBr5svY5aO\nYfaB2WaXJESFVTTYk5RSoQCW29MVL0lcqKikiDkH59AtrBt1feqaXY7dqedTj+kDp9O1Tlfe3PAm\nE2MmyklVYdMqGuzzgNGW70cDv1Vwe+IiouOjSclL4c6md5pdit3ydvHm0z6fMrzJcKbumcoLa1+g\nqKTI7LKEuCplbmNXSv2IcaI0QCmVALwK/BuYrZR6AIgH7qiMImu6H/f/SJhXGN3qdDO7FLvm6ODI\nS51fIsQzhElbJ5Gan8pHvT7C28Xb7NKEKJfy9Iq51OGi9cecFOccTD/I1tNbGd9+vIwLUwWUUjzY\n6kGCPYKZ8OcE7l1yL5/3/Zxgz2CzSxOizOTK02pu5v6ZuDq6ckvDW8wupUa5qcFNfHb9ZyRkJTBy\n8UiOnDlidklClJkEezWWVZjFgtgFDIgaQG232maXU+NcW+davuv/HcWlxYxaPIotSVvMLkmIMpFg\nr8bmHZlHXnGejAtjomb+zfhh4A8EuAcwZukYlh5banZJQlyRBHs1pbVm5v6ZtA5oTQv/FmaXU6OF\neYUxbcA0WgS04OnVT/PD3h/MLkmIy5Jgr6Y2JG7gWOYxOVqvJmq51uKrfl/RJ6IP721+jw9jPpS+\n7qLakmCvpmbun4mvqy83RN5gdinCws3JjQ+v+5ARTUbw3Z7veH7N8xSWFJpdlhD/IFPjVUOJ2Yms\nSljFfS3uk6ncqhlHB0de7PwiIZ4hfLz1Y6Ove++P8HHxMbs0Ic6RI/Zq6KeDPwEwrMkwkysRF6OU\n4oFWD/Buj3fZenoroxeP5lTOKbPLEuIcCfZqprCkkJ8P/UzP8J7U8apjdjniMgbXH8znfT8nMSeR\n4QuG88PeHygoKTC7LCEk2KubpXFLSctP484mMi6MLehapyvTBkyjYe2GvLf5PQbPHczPB3+mqFTG\nmRHmqRHBvilxE9/v+Z5dybuq/R/czP0zqedTjy51uphdiiijRr6N+PbGb/n6hq8Jcg/itfWvMfTX\noSyMXSg9Z4Qp7P7kqdaaV/58hZM5JwFwd3KndWBr2gW1o11wO1oHtMbD2cPkKg37UvexI3kHz3Z8\nFgdVI/7n2pUuoV3oPLAzqxNW8+m2T3l+7fN8s+sbHm37KH3q9pHpDEWVsftg35+2n5M5J3my/ZPU\n8arDtqRtbD29lS92fIFG46gcaebXjLbBbWkf1J62wW3xc/MzpdaZB2bi7uTOzQ1vNmX/ouKUUvSq\n24ue4T35/djvfL79c55Y+QQt/Vvyr7b/omudrhLwotLZfbBHx0fjoBwY2nAofm5+9I805gvNKsxi\nR/IOtiZtZevprczaP4tpe6cBEOkTSfvg9rQNakvviN5V0pUtJS+FRbGLGFR/kHSdswMOyoEBUQPo\nV68f84/M5787/svDyx+mfXB7Hmv7GO2C25ldorBjyowZvDt06KBjYmKqZF+3/HYLtV1rM7X/1Muu\nV1hSyN7UvWxJ2sK208ZRfVZhFuFe4XzW9zPq165faTUezTjKI8sfITkvmZmDZtLQt2Gl7UuYo7Ck\nkDkH5/D1rq9JyUuhe1h3/tX2XzT3b252acKGKKW2aK07XHE9ew72+Mx4Bs0dxLMdn2VU81Hlemyp\nLiXmVAzPrnmWwpJC/tPrP1xb51qr1xhzKobHVz6Ok4MTn/b5lNaBra2+D1F95BXn8eP+H5myewoZ\nBRn0COvB0IZDua7udXIxmriisga7XZ+hWxG/AoA+EX3K/VgH5UCn0E7MGDSDEK8QHln+iNUnOl4Y\nu5Axy8bg5+bHDwN/kFCvAdyd3Lm/5f0svnUxj1zzCAfSDvDU6qfoPbs3b65/kx3JOzDjYEvYF7s+\nYh+1aBQFJQXMvqligZxTlMOza55lTcIaRjYbydMdnq7QbEZaa77a+RWTt0+mQ3AHPu79MbVca1Wo\nRmGbSkpL2Ji4kXmx84iOiya/JJ9In0huanATN9W/iVCvULNLFNVIjW+KSclLoc/sPjzS5hHGXjO2\nwtsrKS3hwy0fMm3vNHqE9eD9nu/j5eJV7u0UlRbx5vo3mXt4LoPqD+KNa9/AxdGlwvUJ25ddmM2y\nuGXMOzKPmKQYFIqOIR0Z0mAI/er1qzbdcoV5anyw/3TwJ95Y/wY/D/mZxr6Nrbbd2Qdm887Gd4iq\nFcVnfT8r12X/WYVZjF81ng2JG3i49cOMazNOur6Ji0rISmB+7HzmH5nP8azjuDu5069eP25qcBOd\nQjrJdQ41VI0P9rHLxxKfGc/CWxZaPTzXn1zPU6uewtnRmU/6fMI1gddc8TGJ2Yk8Ev0IxzKOMaHr\nBG5pJHOYiivTWrM9eTvzjszj96O/k1WURYhnCEMbDuW2RrcR4hlidomiCtXoYM8qzKLnrJ6MbDaS\npzo8VSn7iM2I5dHoR0nKSeKt7m8xIGrAJdfdk7qHR6MfJb84n496f0SXUBkuQJRffnE+qxJW8dvh\n3/jzxJ84KAd61e3F8CbD6RzaWY7ia4CyBrtdXqD0x4k/KC4tvqreMGVVv1Z9pg+czhMrn+DZNc9y\nLOMYY68Z+49PB6uPr+aZNc9Q27U2Xw34ika+jSqtJmHf3Jzc6B/Zn/6R/UnISmDOwTn8cugXouOj\nqedTj2GNh3Fzw5vlRLywzyP2p1c/TcypGFYMW1HpRzGFJYW8vv515h2Zx4CoAbzZ7c1z/ZF/3P8j\n/970b5r6NWVyn8kEegRWai2i5iksKWRp3FJm7Z/F9uTtuDm60T+qPyOajKBFgMyVa29q7BF7QUkB\naxPWMrD+wCr5aOri6MJb3d6ifq36fLz1Y05kn+DjXh8zdc9Upu2dRq/wXrzX8z3p0SAqhYujC4Pr\nD2Zw/cEcSDvArAOzWBC7gF8P/0pL/5YMbzqc/pH9cXNyM7tUUYXs7oh9TcIaxkWP4/O+n9MjvEel\n7ONSlsct54W1L6DRFJQUcFfTu3i247MV6vMuRHllFWYx/8h8Zh+YzZGMI/i4+DC04VCGNRlGPZ96\nZpcnKqDGnjx9bd1rLDm2hDXD15jSP3xP6h5eX/c6QxoMYWTzkVW+fyHO0loTkxTDrAOziI6LplgX\n4+HkgbeLN94u3vi4+ODj4nPu/qWW13atTYhniHTNrQZqZFNMSWkJK4+vpGdYT9Mu+mnh36LCV7oK\nYQ1KGRc4dQzpSHJuMouOLiIpN4nMgkyyCrPIKsriVO4pDp05RGZhJtmF2WgufqAX5BFEl9AudK3T\nlS6hXQhwD6jiZyPKw66CfUfyDtLy0+hTr/J6wwhhiwI9AhndYvRl1ynVpWQXZRuhb/nKLMwkOTeZ\nzac2szphNfOOzAOgYe2G50K+Q3AHOYdUzdhVsEfHR+Ps4Ez3Ot3NLkUIm+OgHM41w1xoRNMRlOpS\n9qftZ/3J9WxI3HBuDgMn5UTrwNbngr5lQEucHOwqWmyO3bSxa60Z8MsA6teqz+fXf27VbQsh/im/\nOJ9tp7exIXEDGxI3sC91HxqNl7MXHUM60jO8J0MbDpWQt6Ia18Z+MP0gJ7JP8FCrh8wuRYgawc3J\nja51utK1TlcA0vPT2XRqExsSN7D+5HpWHl/JgtgFvNfjPYI9g02utmaxm2uQV8SvQKG4ru51Zpci\nRI3k6+bLjZE38mrXV1ly2xLe6f4Oe1P3csf8O1ibsNbs8moUuwn26Pho2ga1lbP1QlQTNzW4iVmD\nZxHoEcgj0Y8wcctEikqLzC6rRrCLYE/ISuBA+oFKHRtGCFF+UbWimD5wOsMaD2Pq7qncu+ReTmaf\nNLssu2cXwV6RKfCEEJXLzcmNV7q+wgfXfcCRM0e4Y/4d5/5mReWwSrArpY4ppXYppbYrpSp/zrsL\nRMdH09i3MXW961b1roUQZdQ/sj8/Df6JcO9wHl/5OO9teo/CkkKzy7JL1jxi7621blOWrjjWlJqX\nyrbT2+gb0bcqdyuEuAp1feoybcA0RjYbyQ/7fmDU4lEczzxudll2x+abYlYdX4VGS7ALYSNcHF14\nrtNzfNz7Y45nHeeOBXew5NgSs8uyK9YKdg0sVUptUUqNudgKSqkxSqkYpVRMcnKylXYLK46vIMwr\nzKrzmgohKl/fiL7MuWkODWo34JnVz/Dm+jfJL843uyy7YK1g76a1bgcMAMYppXpeuILW+iutdQet\ndYfAQOtMOJFTlMP6k+vpE9FHRp4TwgbV8arDd/2/476W9zH74GzuXnQ3sRmxZpdl86wS7Frrk5bb\n08BcoJM1tnsla0+spai0SJphhLBhzg7OjG8/ns/7fk5ybjLD5g/jm13fUFQifd6vVoWDXSnlqZTy\nPvs9cAOwu6LbLYsV8Svwc/OjTWCbqtidEKIS9QjvwZwhc+gR1oNJWycxbMEwtiZtNbssm2SNI/Zg\n4A+l1A5gE7BQa13pZ0KKSopYm7CWXnV7yQxFQtiJII8gPur9EZP7TCanKIfRS0bz6rpXOZN/xuzS\nbEqFBwHTWscC11ihlnLZeGoj2UXZ0gwjhB26ru51dAzpyBc7v+D7Pd+zMn4lT3V4iiENhsj5tDKw\n2e6O0fHReDh50Dm0s9mlCCEqgYezB+Pbj2f2TbOJ8Ing5T9f5oGlD8jJ1TKwyWAv1aWsjF9J97Du\nuDq6ml2OEKISNfZtzPcDvmdC1wnsT9vPbfNuY/K2ydI18jJsMth3Ju8kNT9VmmGEqCEclAN3NL6D\neUPn0T+yP1/u/JJb593KupPrzC6tWrLJYI+Oj8bJwYke4T3MLkUIUYUC3AN4t8e7fNXvKxyUAw8v\ne5hn1zxLSl6K2aVVKzYX7FprouOj6RzSGW8Xb7PLEUKYoGudrvw85Gf+75r/Y3nccobMHcLcQ3Mx\nY6rP6siG0OaQAAAgAElEQVTmgv3wmcMczzouQ/QKUcO5OrrySJtH+HnIzzT1b8qEdROYsG6CtL1j\ng8EeHR+NQkmwCyEAYzKPr/t9zdhrxvLr4V+5Z/E9JGQlmF2WqWwu2FfEr+CawGtkCjwhxDmODo6M\nazOOyX0mk5CdwPAFw1mTsMbsskxjU8F+Mvsk+9L2ydG6EOKirqt7HbMGzyLUM5RHox/l8+2fU6pL\nzS6rytlUsJ+dTku6OQohLqWud12mDZzGTQ1u4r87/su46HFkFGSYXVaVsqlgzyrMonVAayJ8Iswu\nRQhRjbk7ufNWt7d4pcsrbEjcwPAFw9mbutfssqqMMqN7UIcOHXRMzNVNjaq1lrEihBBltjN5J+NX\njSc9P52Xu7zMLY1uMbukq6aU2lKW6Udt6ogdkFAXQpRL68DWzL5pNm2D2zJh3QReW/caBSUFZpdV\nqWwu2IUQorz83Pz48vovebDVg/x86GdGLx7NyeyTZpdVaSTYhRA1gqODI4+3e5xJvScRlxnHsAXD\nWHfCPseakWAXQtQofSL6MHPwTALdAxm7fCxf7fzK7rpESrALIWqcej71mD5wOgOiBvDptk95fMXj\nZBZmml2W1UiwCyFqJA9nD/7d49+80OkF/jjxByMWjOBA2gGzy7IKCXYhRI2llOKuZncxtf9UCooL\nGLloJPOOzDO7rAqTYBdC1Hhtgtow66ZZtApsxUt/vMRbG96isKTQ7LKumgS7EEJgTOLxVb+vuK/F\nfcw6MIt7l9zLqZxTZpd1VSTYhRDCwsnBifEdxjOx10RiM2IZNn8YGxI3mF1WuUmwCyHEBfrV68eP\ng37Ez82Ph5c9zDe7vrGp2Zkk2IUQ4iKiakUxY9AMbqh3A5O2TuLxlY+TVZhldlllIsEuhBCX4OHs\nwfs93+e5js+xNmEtIxaM4GD6QbPLuiIJdiGEuAylFCObj+TbG78lrziPkYtGsiB2gdllXZYEuxBC\nlEG74HbMvmk2zf2b88LaF5i0dVK1bXeXYBdCiDIKcA/g6xu+5rZGt/HNrm94e+Pb1XKcGSezCxBC\n1GCFuZC0G05uh8Ttxq0ugY4PQpu7wcXD7Ar/wdnBmVe7voqPiw9T90wlpyiHN7q9gbODs9mlnSPB\nLoSwupJSTXJWAS5ODni4OOLq5IAqyoNTu/4K8MTtkLwfzh7xegRAnTaQlw6LnoZV70LnsUbIe/iZ\n+4QuoJTiyfZP4u3izSfbPiGnKIcPrvsAV0dXs0sDbHBqPCFE9XMmt5Bt8WfYGp/O1vh0dh1Pp2Hh\nflo7xNLK4Sgt1VEaqhM4KiNv0lVtYp0bEu/amBPuTTjt1YxCjxDcXZ3wc3emk+MBWh37Do+45eDs\nAe1GQ9dHoHb1m+94xr4ZvLvpXTqHduaT3p/g4Vx5nzLKOjWeBLsQolxKSjWHTmexNe6vII9NzgEg\nQGXySO313FL8O75FxuX4uc7+nPJsSoJ7Y465NOKwU0MSS/3ILy4lt7CE3MIS8gqLLbclZBUUn9tX\nK+cTPOGxhF4Fq1AKkusNwrH7E/g3aFetpsmcd2Qer/z5Ci0DWvJ538+p5VqrUvYjwS6EsIqM3CK2\nHk9nW1w6W+PPsOP4mXPh6+vhTLsIX/r7nuC6jF8JjFuEKimAyB7Q/l6ody14h0I5Qji3sJhDSdkc\nOJXFgaQsDiZlkZ54lJvzf+NOxxV4qXz+oA3RvndSFHEtTUJr0STYmwaBnvh5upgW+NFx0Tyz5hmi\nakXxZb8vCXAPsPo+JNhFpUjLKSQpMx8w/lYV6tzfrOLs3++FyxQOCsJ9PXB0qD5HWeLSjiRns2xv\nEsv2JrE1Ph2twUFBkxAf2kXUpl2EL+3D3KmXuAS1+Rs4uRVcvOCaO4028aCmVq8pLaeQI/EJOG6Z\nQuNjP+BVnM4u3YDPiwbze2lHSjHa88N93anr60FdPw/Cfd0J9/Wgrp87df088HGr3BOc606u44mV\nTxDkEcTX/b4m1CvUqtuXYBdWkV1QzKajqfx5OJU/D6ew/9TVX1LdMsyH/97dnrp+1a+nQ01XUqrZ\nGp/OckuYx6YYTSst6vjQt1kwXaL8aF23Nl6uTpAeBzFTYOv3kJcGAU2g00PQeji4+VRNwUX5sGMG\net2nqLRYcr3qsTv0NmKc27EtL4Tj6XkkpOeRfV6zDkAtd+fzgt8I+2AfNwK8XAn0ciXA2wUPl4r1\nKdl2ehvjlo/D08WTr/t9TWStyApt73wS7OKqFBSXsDXuDOuOpPDn4RR2JGRQUqpxcXKgQz1fujUM\noH6AJwBn3zlag0Zz9q2k4W8XbmhtnFz7cNlBHB0UHw9vQ68mQVX7xMQ/5BYWs/ZQCsv2JrFi/2nS\ncgpxdlR0qe9Pv+bB9G0WTFhtd2Pl0lKIXQmbv4GDS4xlTQdBx4cgqme5mlqsqrQE9i+APyfBiS3G\nMp8waNAb3aAvmaHdiM9z43h6LgnpuRxPy7N8n0dCei75Rf/sg+7h4kiAlysBXi7GrberJfj/uh/o\n5Updv0t/At2Xuo+xy8cC8FW/r2ji18QqT1eCXZRJSalm94kM/jySwvojqWw+lkZ+USkOClqH16Zb\nQ3+6NQigXT1f3JwdK7SvYyk5jP1hCweSsni8byMe69MIB2maqVKns/JZse80y/Ym8cfhFAqKS/Fx\nc6J30yCubxbMdU0C/95cUZgDW6fB5q8h9TB4Bho9VDrcB7XCzXsiF3PmOByJhsPRELsaCjJAOUCd\ndtCwLzToC2HtwdE4Itdak5xdwOnMApKzC0jJKiAlu5CU7IK/vrKM+2m5hVwYlZ4ujrSN8KVdPV/a\n1/OlTd3a1HL/67WLzYhlzNIx5Bbn8nnfz2kT1KbCT7FKg10p1R+YBDgC32it/3259SXYzbfpaBrf\nrI1lQ2wqmfnGx9Umwd5c29CfaxsE0Lm+3z/bI7OSYO4YSDsKjfsbR2z1rgXHsrdb5hWW8NKvu/hl\n6wmuaxzIpBFtqO3hYs2nJi5w8kwei3YlsmhXItuOn0FrCPd1p1/zYPo1C6ZjlB/OjhdchF6Yaxyd\n/zkJclMgvJPR3NL8ZnCqHn21L6uk2DiCPxv0J7YAGtxqQdR1fwV97bpl2lxxSSlpuYXngv5UZj67\nT2SwJS6dfYmZlGrjQ0vjIO9zQd++ni8urmcYs2wMyXnJTOo9ia51ulboaVVZsCulHIGDQD8gAdgM\n3Km13nupx0iwm2vOlgSe/3kn/l4u9G4SRNcGRpgHel/mD/bEFpg5EvLPGGF+7A8ozge32paQH2j8\nobh6XXH/Wmumb4znjfl7CfJx5YuR7WkZVjndw2qqUxn5LNqVyMJdiWyJSweM9vL+LUK4vnkwTUO8\nL957pDDXaD//82PISYb6vaHXCxDRuYqfgZXlpkHsKkvQr4Csk8bygMZGU1LdzsZX7YhyNytlFxSz\n8/gZtsSlExNndP/Mshws+Xm60DLCgTjniWQWJ/Lute/QP+r6c58ayqsqg70r8JrW+kbL/RcAtNbv\nXuoxEuzm0Frz0bKDfLLiMN0bBvDZ3e3+9tHxknbMhHmPgXcwjJgBIa2Mj+hHVsD+hXBgsRH4jq7Q\noLdxJN94AHgFXnaz24+f4ZEftpCSU8ibN7dgeMfzLj7Jz4SMBMg8CWjjKNHJDRxdjFsny62j63k/\nq9kXUidl5rPYEuabjxlh3jzUh0GtQxnUKpRIy7mRiyrKg5ipRqBnJxlHtb1fhIguVVR9FdLauOL1\ncLQR9Mc3QWG28TPvUKjbyRL0XSC0dbk+kQKUlmoOn85i74H9pBzZSumpPXgWHGROnZMcddU8FjKa\n+/s/c1WlV2Ww3w7011o/aLk/CuistX70gvXGAGMAIiIi2sfFxVVov6J8CopLeHbOTn7bfpJhHcJ5\n+5ZW//z4faGSYlj+KqyfbPRLvuN/4Ol/8fXi1xshv38hZMQDyvjjaDrI+PJvcN76RZCVCBkJZJ8+\nxsK1mylMi6eDbw5N3DNwyDhhtI+Wl3L4K/xdvSG8g3HEWb8X+NYr//ZswOmsfJbsPsWCnYlsPpaG\n1tA0xJtBrUIZ2DqUBoFX+ARVlA9bvoM/PoLsU8bvufeLxqeymqKkGE7vheMbja/4jZb3MODkDmHt\n/jqir9vpn8Mb5GfA6X2QtMfYTtJeOL3HWH52F951OOFZn/84F/N4j9dp0Pjq/mFWZbDfAdx4QbB3\n0lr/61KPkSP2qpWeU8jD07aw6Vgaz9zYhEd6NbjyRRy5aTDnfqMnRKeH4ca3y3bkorUxHsiBRUZv\nhVO7jOWBTY32zYwEI9QvGBEvz8mH2EJfslxDaNGsBd7BkcbJOZ8wI7CLC4yvkgKjCejs/XPLLlie\nlwZx64x9AfhGGQHfoLcRXtVs7JHySM4qYMmeUyzceZKNR40wbxzsxaBWdRjUOoSGQd5X3khRvtFd\n8Y+JxmtUrzv0fgEiu1f+E7AFmSctQb8J4jfAqZ1Qauk6GdAYwjoY77GkvX/9EwBw9YGgZhDcAoKa\nW26bgbuvVcqSphgBGD1R7v9uMwnpeXxwR2tubhN25Qed3gc/3mmE8OCJ0O6eqy8gPc5oqjm42Oia\nViv8gq+6Rni7erF8bxJPzt6Og1J8PKINvSvaJVJrSD5gtK3GrjLOCxRmAcoYbKp+L+OIvm5ncHar\n2L4qkdaaQ0fjOByzlOKjf1IvZydOlJDlHIhHQF3qRNQnIDQKfELBu45x61b74m3FxQVGoK+daLQz\nR1xrBHpUz6p/YrakMNe4COts2J/YCh7+RnAHN4cgy22tupXa9bMqg90J4+RpX+AExsnTu7TWey71\nGAn2qrElLo2Hvt+C1pqv7ulAx8gyHKXuXwi/jDEGXhr+Q5WfNDu/S+RjfRrxeF8rdoksKTL+IGNX\nGkGfsNk4CnNyg4iuRtBHdofa9Yw/WgfzpisoSE/gyOZlZB9cTWBqDFH6uLEcF07XaoVvrVp4FpxG\nZSVCbuo/N+Dk/veg9w4FF0+j62JmgtF+3PsFoy29Go25Ii6vqrs7DgQ+xujuOEVr/fbl1pdgr3zz\nd5zkqZ92EFbbnSn3diTqcifOwLgAZc0HsOodo9/v8B+gVhmO7ivBhV0ib20XhoeLEx4ujrg5O+Lh\nYny5uzji4eKEu7Pj1Q1VUJBlNNfEroIjKyF5318/c3QxwtAnDHws4Xjue8utVzA4VKxvP2C5giuO\nrAOrSd6zEq9TmwgqOgFAtnbnqEdLSuteS912/fBr2Nk4cXy+4gKjOSUz0TgKz0y03D/599uSQqPb\nYu8XjE8qEug2Ry5QqqG01vx39RHeX3KAjpG+fDWqA76eV+gnXpANv46FffOh9Qi4aZLpTRPnd4ks\nLLnyDDVnx/32cDYC38vViXr+njQK8qJhkBeNgr2o5+95+RPGWacgIcYIwswTltvzvi8p+Pv6yhG8\nQ4yQ9/AHByfLl+Nf3yuHSy7TyoHMU0dwiF+Hd0ESAOnai50OzckO7URQyz60at8dN1cr9BvXGgoy\njTZgCXSbJcFeAxWVlPLKr7uZufk4Q66pw/u3t77y1aJpR2HmXUb3rxvegi6PVKs//IzcIpKzC8gr\nLCG3sJjcohLyLUO95hadN9xrUYllHeM2M7+I2OQcTpzJO7ctJwdFZIAnDQONoG9oCf0GgV5Xfp20\nNk4onwv8C4I/L8341FNabHzpEuOcQmnxuVutS9AlxZSWGMscKCFF12JTaVMSfNri1fg62rTvQouw\n2tVqSFpRfZQ12Gt2x187kplfxLjpW1l7KIVHezdkfL/GV26bjl0FP91rhNbIn6FBn6ootVxqeThT\ny+PqR+TLKSgmNjmHw8lZHErK5vDpbA4mZbF07ylKLcc0SkFdX49zQR/p70k9fw/q+XsQWsvdaOZR\nyujq6elv9G0ug9JSzf5TWWw8msqmo2lsOppGak4hAEHernSu78+1Dfzp0zSIwT7V9+StsD0S7Hbg\nxJk87p+6mSPJ2bx/e2uGdSjDZdLbpsO8f0FAI+Oio/P7mdsRT1cnWoXXolX4369sLSgu4VhKLodO\nZ3H4dDaHTmdz5HQ2fxxK+VvTj4ujA3X93C1h/1fgR/p7Eubr/remneKSUvaczPxbkJ8driGstjvX\nNQmkc5QfnaP8qefvIUflotLYVLCfPJNHcYkmwl+GfT1rx/EzPPh9DPmFJfzv/k50a1iGwf0PLzdC\nPaqHcZLUtQz9nu2Mq5MjTUK8aRLy9+deUqo5lZlPXEoOcWm5HEvNIS7FuF13JJW8opJz6zo6KMJ9\n3YmwDEO8NS6dnELj51EBngxsFUqnKD86RfkR7ivvWVF1bCrY//P7AebvPMldnSJ4tE+jy49tUgPM\n3ZbAcz/vItDLlemPdKZxcBkC+tRumH2vcfFEDQ31y3F0UITVdiestjsXXnuptTFBc1xaLsdScohL\ntQR/ai5FJaXc0i6MzlH+dI7yI0iaVoSJbCrYnxvQFDcXR37YGM9PWxJ4sEd9HuoRhXclz4pS3ZSU\nat5fsp8v18TSOcqP/45sj9+Ver6AcaJvxjAjzO+eLaFeTkopgnzcCPJxK9s1AUKYxCZ7xcQmZ/Ph\n0oMs3JWIn6cL43o3ZGSXCFydrNCnuJrLyCvi8ZnbWHUgmVFd6jHhpuZXHvMFjD7bUwcYvWDuX2IM\n5CWEsCk1orvjzoQzvLdkP38eTiWstjvj+zVmaNswu51X80hyNg/9L4b4tFzeuLkld3WOuPKDwBjk\naOadxmh2d82GRtdXbqFCiEpR1mA375ppK2gdXpvpD3Zh2gOd8PV05qmfdjDok7VE70vCjH9YlWnl\ngdMM/exPMvKKmPFQl7KHutaw+Fk4tBQG/UdCXYgawKaD/awejQKZN647k+9qS35RCQ/8L4ZhX65n\nS1ya2aVVmNaaL1cf4f7vNlPX14PfHu1Gp6hytO+unwwx30K3x6HD/ZVXqBCi2rDpppiLKSopZdbm\n40yKPkRyVgHXNwvm2f5NytZjpJrJLyrhhV92MXfbCQa1CuWDO1qXbwb1vb/B7NHGdGa3TzV1UCsh\nRMXViDb2y8ktLGbKH0f5cnUs2YXF3NYunKduaExoLfdK3a+1nMrIZ8y0GHYmZPD0DY0Z17th+S5o\nOb4Z/jcYQlrD6HngbBvPWwhxaTU+2M9Kzynks5WH+X59HErB/d2j+L9eDf45UXM1sjU+nYenbSG3\noJiPR7SlX/Pg8m0g7Sh8c70x/+iD0eBZhouWhBDVngT7BY6n5TJx2UHmbjuBr4cz/+rTiLurYRfJ\nn2KO89Lc3YTWduPrezqUvwkpNw2+vcGYiPjB5caQAUIIu2CfvWLOxMPuX67qoXX9PPhoeBsW/Ks7\nzev48MaCvVw/cTXzdpyktNT8HjS5hcW8MX8vz8zZSacoP34b1638oV5cALNGwZk4Y/wXCXUhaiTb\nOmKf+3+w+2d4dBP4Rl71/rXWrDmUwruL9rH/VBatw2vxwoBmdG1wkYmaK9mJM3l8v/4YP26MJzO/\nmPu6RfLSwGY4leWio/NpDXMfhp2z4NZvoPUdlVKvEMI89tkUk3ECJneAhn2NcU4qqKRUM3fbCSYu\nPcDJjHz6NA3iuf5N/zEwlLVprdkan86UP46xZM8ptNYMaBnK/d0jaV/vKi9VX/kOrH4P+rwMPZ+x\nbsFCiGrBPoMdjOnbVrwF98yD+tdZpZ78ohK+W3eMz1YeJqegmNvbhzO+XxNCall3IKfC4lIW705k\nyh9H2ZGQgY+bE3d2imBU13oVG/1v+wz49f+g7UgYMrlaTZQhhLAe+w32onz4rCO4eMHDa8HReuOY\npecUMnnlYaatj8PBAR7oHsVt7cIJ9/XAxenqT0ek5RTy46Z4vl9/jKTMAuoHenJftyhus8zlWSFH\n18K0ocYkzHfPAcfq29tHCFEx9hvsAHvnwexRMPA/0Okh6xVmcTwtl/8sPcBv208C4KAgzNeYbCHC\nz+PcDDuRAcb9S02rduBUFlP/PMrcbScoKC6lR6MA7u8exXWNAq88u1FZpMfBV72M7owPLge3Wld8\niBDCdtl3sGsN3w+BxJ3w2DbwqJwhVA8lZbEjIYM4y5jbcak5HE3JOTcrzlkhPm7nZtWJ8PcgwMuF\n+TsS+eNwCq5ODtzaLpz7ukVa9+rXwhyjW2PGcXhopd3OgCSE+It9z3mqFPR/D77oDivfhkEfVspu\nGgV70+giYXwmt/BvkyycDf3o/adJyTZmsg/xceOZG5twV6cIfMsyVnp5aA2/PgKn98JdP0moCyH+\nxjaDHSC4OXR8ADZ/YwxuFdyiynZd28OF2h4uXFO39j9+ll1QzKmMPOr5e5ZtnPSr8cdE2PsrXP+6\njNYohPgH27pA6UK9XjDalRc/ZxzFVgNerk40DPKuvFA/+DtEvwktbzdGbBRCiAvYdrB7+EHvl+DY\nWtg3z+xqKl/yQfj5QWP2oyGfSrdGIcRF2XawA7S/D4JawO8vQ1Ge2dVUnvwMmHkXOLoYwwW4yKz3\nQoiLs/1gd3SCAf+GjHhYN9nsaipHaQn8/BCkH4Vh30PtumZXJISoxmw/2AGiekKzIcZJxYwTZldj\nfSvfhkO/w4D3ILKb2dUIIao5+wh2gBveAl0KyyaYXYl17f4F1n4I7UZDhwfMrkYIYQPsJ9h968G1\nj8HuORC33uxqrOPULvhtHNTtbFxlKydLhRBlYD/BDtD9CfAJg8XPGu3StiwnFX68C9xqw7Bp4GTl\ni5yEEHbLvoLdxRP6vQGndsK2ig/ra5qSIvhpNGQnwYjp4F3OqfGEEDWafQU7QMvbIKIrRL8BeWfM\nrubq/G7pmz/kEwhrZ3Y1QggbY3/BrpTReyQ3FVa/b3Y15bftB9j0JXQZB9eMMLsaIYQNsr9gBwi9\nBtrdYwRk8kGzqym745thwZNQv5fRpCSEEFfBPoMdoM8r4OwJv79QbcaRuay0ozBrJPjUgdunWnUC\nESFEzVKhYFdKvaaUOqGU2m75GmitwirMKxB6PQeHlxsDZ1VnJ7bAt/2gpABG/Fhp48sLIWoGaxyx\nf6S1bmP5WmSF7VlPpzEQ0Ng4ai8uMLuaizuwBL4bDM4e8MAyYzhiIYSoAPttigFj/s/+70JaLPz5\nidnV/NPmb2HmnRDYxJjaLqCR2RUJIeyANYL9UaXUTqXUFKWU76VWUkqNUUrFKKVikpOTrbDbMmp4\nPTS/GVa+BSvegtLSqtv3pWgNy1+HheOhYT+4dyF4BZldlRDCTlxxzlOl1HIg5CI/egnYAKQAGngT\nCNVa33+lnVZ4ztPyKi6EhU8aXQlb3AJD/wvO7lW3/wtr+W0c7JptDDk88D9yolQIUSZWm/NUa12m\nudeUUl8DC8qybpVzcoEhk4329mWvwpnjxpjmVX1FZ36G0fPl6BroOwG6j5fxX4QQVlfRXjGh5929\nBdhdsXIqkVLGVHLDfzAmgf6mLyTtqbr9ZyTAlP4Qtw5u+RJ6PCWhLoSoFBVtY39fKbVLKbUT6A08\naYWaKlezwXDfYigthm9vgINLK3+fp3bDN/2McB/5s1xRKoSoVBUKdq31KK11K611a631EK11orUK\nq1R12sBDK8CvPvw4HDZ8UXkXMR1ZaRypA9y/xLiqVAghKpF9d3e8HJ86RtA2GQhLnoNFT0NJsXX3\nsWMmTL8dakcY3RmDW1h3+0IIcRE1N9jBGOZ32DSj7X3zNzDjDuMEZ0VpDWs+gLkPQ71r4f7FUCus\n4tsVQogykH52Dg7GgFv+DY0BuL69Ae6aBb6R5d9W1ilI3GnM4rRzFrQebvTGkUkyhBBVSIL9rHb3\nGGE+axR83dfoDhnR+eLrag3px4wJPRJ3GGF+aqcxMQYACno8DX1elp4vQogqJ8F+vqie8GC00STz\nv8Fw82fQ4lZIOWgJcUuQn9oFBZYmG+UIQc2gQV8IbW0MGRzcEtx8zH0uQoga64pXnlaGKr/ytLxy\n04wj97g/wMkNivON5U7uxgnQswEe0hqCmoOzm7n1CiFqBKtdeVojefjBqLnwx0dQkGkEeGhr8G8k\nl/8LIao9SalLcXIxxnMXQggbU7O7OwohhB2SYBdCCDsjwS6EEHZGgl0IIeyMBLsQQtgZCXYhhLAz\nEuxCCGFnJNiFEMLOmDKkgFIqGYi7yocHYEygLf4ir8nFyevyT/Ka/JMtvSb1tNaBV1rJlGCvCKVU\nTFnGSqhJ5DW5OHld/klek3+yx9dEmmKEEMLOSLALIYSdscVg/8rsAqoheU0uTl6Xf5LX5J/s7jWx\nuTZ2IYQQl2eLR+xCCCEuQ4JdCCHsjE0Fu1Kqv1LqgFLqsFLqebPrqQ6UUseUUruUUtuVUtV4vsHK\no5SaopQ6rZTafd4yP6XUMqXUIcutr5k1VrVLvCavKaVOWN4r25VSA82ssaoppeoqpVYqpfYppfYo\npR63LLe794rNBLtSyhH4DBgANAfuVEo1N7eqaqO31rqNvfXFLYfvgP4XLHseiNZaNwKiLfdrku/4\n52sC8JHlvdJGa72oimsyWzHwlNa6GdAFGGfJELt7r9hMsAOdgMNa61itdSEwE7jZ5JpENaC1XgOk\nXbD4ZuB/lu//Bwyt0qJMdonXpEbTWidqrbdavs8C9gFh2OF7xZaCPQw4ft79BMuymk4DS5VSW5RS\nY8wuphoJ1longvEHDQSZXE918ahSaqelqcbmmxyullIqEmgLbMQO3yu2FOzqIsukryZ001q3w2ii\nGqeU6ml2QaLa+i/QAGgDJAIfmluOOZRSXsDPwBNa60yz66kMthTsCUDd8+6HAydNqqXa0FqftNye\nBuZiNFkJSFJKhQJYbk+bXI/ptNZJWusSrXUp8DU18L2ilHLGCPXpWutfLIvt7r1iS8G+GWiklIpS\nSrkAI4B5JtdkKqWUp1LK++z3wA3A7ss/qsaYB4y2fD8a+M3EWqqFs+FlcQs17L2ilFLAt8A+rfXE\n835kd+8Vm7ry1NI962PAEZiitX7b5JJMpZSqj3GUDuAEzKiJr4lS6kegF8bwq0nAq8CvwGwgAogH\n7onQZwAAAABjSURBVNBa15iTiZd4TXphNMNo4Bjw8Nm25ZpAKdUdWAvsAkoti1/EaGe3q/eKTQW7\nEEKIK7OlphghhBBlIMEuhBB2RoJdCCHsjAS7EELYGQl2IYSwMxLsQghhZyTYhRDCzvw/8BY7vAZ6\nYu8AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f6345307278>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Define\n",
    "nsamples,nvars = cdata.shape\n",
    "\n",
    "# Cummulative sum after removing mean\n",
    "cdata = cdata-cdata.mean(axis=0)\n",
    "xx = np.cumsum(cdata,axis=0)\n",
    "plt.plot(xx)\n",
    "plt.title('Cummulative sum')\n",
    "plt.legend(['$x_1$','$x_2$','$x_3$'])\n",
    "xx.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Return sliding windows\n",
    "def sliding_window(xx,k):\n",
    "    # Function to generate boxes given dataset(xx) and box size (k)\n",
    "    import numpy as np\n",
    "    # generate indexes. O(1) way of doing it :)\n",
    "    idx = np.arange(k)[None, :]+np.arange(len(xx)-k+1)[:, None]\n",
    "    return xx[idx],idx"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Compute Eq.4\n",
    "k = 4\n",
    "F2_dfa_x = np.zeros(nvars)\n",
    "allxdif = []\n",
    "for ivar in range(nvars): # do for all vars\n",
    "    xx_swin , idx = sliding_window(xx[:,ivar],k)\n",
    "    nwin = xx_swin.shape[0]\n",
    "    b1, b0 = np.polyfit(np.arange(k),xx_swin.T,deg=1) # linear fit\n",
    "    #x_hat = [[b1[i]*j+b0[i] for j in range(k)] for i in range(nwin)] # slow version\n",
    "    x_hatx = repmat(b1,k,1).T*repmat(range(k),nwin,1) + repmat(b0,k,1).T\n",
    "    # Store differences to the linear fit\n",
    "    xdif = xx_swin-x_hatx\n",
    "    allxdif.append(xdif)\n",
    "    # Eq.4\n",
    "    F2_dfa_x[ivar] = (xdif**2).mean()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 1.        ,  0.62664048,  0.62095623],\n",
       "       [ 0.62664048,  1.        ,  0.1783183 ],\n",
       "       [ 0.62095623,  0.1783183 ,  1.        ]])"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Get the DCCA matrix\n",
    "dcca = np.zeros([nvars,nvars])\n",
    "for i in range(nvars): # do for all vars\n",
    "    for j in range(nvars): # do for all vars\n",
    "        # Eq.5 and 6\n",
    "        F2_dcca = (allxdif[i]*allxdif[j]).mean()\n",
    "        # Eq.1: DCCA\n",
    "        dcca[i,j] = F2_dcca / np.sqrt(F2_dfa_x[i] * F2_dfa_x[j])\n",
    "\n",
    "dcca"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-1.        ,  0.66890233,  0.66406175],\n",
       "       [ 0.66890233, -1.        , -0.34508556],\n",
       "       [ 0.66406175, -0.34508556, -1.        ]])"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Get DPCCA\n",
    "C = np.linalg.inv(dcca)\n",
    "dpcca = np.zeros([nvars,nvars])\n",
    "for i in range(nvars):\n",
    "    for j in range(nvars):\n",
    "        dpcca[i,j] = -C[i,j]/np.sqrt(C[i,i]*C[j,j]);\n",
    "dpcca"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 1.        ,  0.84980211,  0.49768411],\n",
       "       [ 0.84980211,  1.        , -0.45946314],\n",
       "       [ 0.49768411, -0.45946314,  1.        ]])"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Compute Corr (for the parCorr)\n",
    "corr = np.corrcoef(cdata.T)\n",
    "cov = np.cov(cdata.T)\n",
    "# Get parCorr\n",
    "C0 = np.linalg.inv(cov)\n",
    "mydiag = np.sqrt(np.abs(np.diag(C0)))\n",
    "# Efficient computation\n",
    "(-C0/repmat(mydiag,nvars,1).T)/repmat(mydiag,nvars,1)+2*np.eye(nvars)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "# 2. Function to compute DPCCA"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Return sliding windows\n",
    "def sliding_window(xx,k):\n",
    "    # Function to generate boxes given dataset(xx) and box size (k)\n",
    "    import numpy as np\n",
    "\n",
    "    # generate indexes! O(1) way of doing it :)\n",
    "    idx = np.arange(k)[None, :]+np.arange(len(xx)-k+1)[:, None]\n",
    "    return xx[idx],idx\n",
    "\n",
    "def compute_dpcca_others(cdata,k):\n",
    "    # Input: cdata(nsamples,nvars), k: time scale for dpcca\n",
    "    # Output: dcca, dpcca, corr, partialCorr\n",
    "    #\n",
    "    # Date(last modification): 02/15/2018\n",
    "    # Author: Jaime Ide (jaime.ide@yale.edu)\n",
    "    \n",
    "    # Code distributed \"as is\", in the hope that it will be useful, but WITHOUT ANY WARRANTY;\n",
    "    # without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. \n",
    "    # See the GNU General Public License for more details.\n",
    "    \n",
    "    import numpy as np\n",
    "    from numpy.matlib import repmat\n",
    "    \n",
    "    # Define\n",
    "    nsamples,nvars = cdata.shape\n",
    "\n",
    "    # Cummulative sum after removing mean\n",
    "    #cdata = signal.detrend(cdata,axis=0) # different from only removing the mean...\n",
    "    cdata = cdata-cdata.mean(axis=0)\n",
    "    xx = np.cumsum(cdata,axis=0)\n",
    "    \n",
    "    F2_dfa_x = np.zeros(nvars)\n",
    "    allxdif = []\n",
    "    # Get alldif and F2_dfa\n",
    "    for ivar in range(nvars): # do for all vars\n",
    "        xx_swin , idx = sliding_window(xx[:,ivar],k)\n",
    "        nwin = xx_swin.shape[0]\n",
    "        b1, b0 = np.polyfit(np.arange(k),xx_swin.T,deg=1) # linear fit (UPDATE if needed)\n",
    "        \n",
    "        #x_hat = [[b1[i]*j+b0[i] for j in range(k)] for i in range(nwin)] # Slower version\n",
    "        x_hatx = repmat(b1,k,1).T*repmat(range(k),nwin,1) + repmat(b0,k,1).T\n",
    "    \n",
    "        # Store differences to the linear fit\n",
    "        xdif = xx_swin-x_hatx\n",
    "        allxdif.append(xdif)\n",
    "        # Eq.4\n",
    "        F2_dfa_x[ivar] = (xdif**2).mean()\n",
    "    # Get the DCCA matrix\n",
    "    dcca = np.zeros([nvars,nvars])\n",
    "    for i in range(nvars): # do for all vars\n",
    "        for j in range(nvars): # do for all vars\n",
    "            # Eq.5 and 6\n",
    "            F2_dcca = (allxdif[i]*allxdif[j]).mean()\n",
    "            # Eq.1: DCCA\n",
    "            dcca[i,j] = F2_dcca / np.sqrt(F2_dfa_x[i] * F2_dfa_x[j])   \n",
    "    \n",
    "    # Get DPCCA\n",
    "    C = np.linalg.inv(dcca)\n",
    "    \n",
    "    # (Clear but slow version)\n",
    "    #dpcca = np.zeros([nvars,nvars])\n",
    "    #for i in range(nvars):\n",
    "    #    for j in range(nvars):\n",
    "    #        dpcca[i,j] = -C[i,j]/np.sqrt(C[i,i]*C[j,j])\n",
    "    \n",
    "    # DPCCA (oneliner version)\n",
    "    mydiag = np.sqrt(np.abs(np.diag(C)))\n",
    "    dpcca = (-C/repmat(mydiag,nvars,1).T)/repmat(mydiag,nvars,1)+2*np.eye(nvars)\n",
    "    \n",
    "    # Include correlation and partial corr just for comparison ;)\n",
    "    # Compute Corr\n",
    "    corr = np.corrcoef(cdata.T)\n",
    "    # Get parCorr\n",
    "    cov = np.cov(cdata.T)\n",
    "    C0 = np.linalg.inv(cov)\n",
    "    mydiag = np.sqrt(np.abs(np.diag(C0)))\n",
    "    parCorr = (-C0/repmat(mydiag,nvars,1).T)/repmat(mydiag,nvars,1)+2*np.eye(nvars)\n",
    "    \n",
    "    return corr,parCorr,dcca,dpcca"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Pearson:\n",
      " [[ 1.          0.80626524  0.22904353]\n",
      " [ 0.80626524  1.         -0.0799021 ]\n",
      " [ 0.22904353 -0.0799021   1.        ]]\n",
      "PartialCorr:\n",
      " [[ 1.          0.84980211  0.49768411]\n",
      " [ 0.84980211  1.         -0.45946314]\n",
      " [ 0.49768411 -0.45946314  1.        ]]\n",
      "DCCA(k=6):\n",
      "[[ 1.          0.8198706   0.82176979]\n",
      " [ 0.8198706   1.          0.60633106]\n",
      " [ 0.82176979  0.60633106  1.        ]]\n",
      "DPCCA(k=6):\n",
      "[[ 1.          0.70974722  0.71306451]\n",
      " [ 0.70974722  1.         -0.20663271]\n",
      " [ 0.71306451 -0.20663271  1.        ]]\n"
     ]
    }
   ],
   "source": [
    "# Example using the function\n",
    "k = 6\n",
    "corr,parCorr,dcca, dpcca = compute_dpcca_others(cdata,k)\n",
    "print('Pearson:\\n',corr)\n",
    "print('PartialCorr:\\n',parCorr)\n",
    "print('DCCA(k={}):\\n{}'.format(k,dcca))\n",
    "print('DPCCA(k={}):\\n{}'.format(k,dpcca))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}