Comparing DS with RA

compare with recurrent analysis.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Comparison of DynamicalSystems.jl and RecurrenceAnalysis.jl"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "using DynamicalSystems, RecurrenceAnalysis, BenchmarkTools"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the following RA stands for RecurrenceAnalysis while DS stands for DynamicalSystems (including ChaosTools).\n",
    "\n",
    "# 1. Delay Embedding\n",
    "\n",
    "## 1.1. Methods give same result\n",
    "\n",
    "First confirm that both methods give the same result"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "towel = Systems.towel();\n",
    "N = 1000\n",
    "x = trajectory(towel, N)[:, 1]\n",
    "const τ = 2\n",
    "D = 1 # for DS this is the amount of temporal neighbors, symbol subject to change"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2-dimensional Dataset{Float64} with 999 points\n",
       " 0.085     0.76827 \n",
       " 0.285813  0.681871\n",
       " 0.76827   0.837347\n",
       " 0.681871  0.51969 \n",
       " 0.837347  0.966676\n",
       " 0.51969   0.112748\n",
       " 0.966676  0.386547\n",
       " 0.112748  0.910741\n",
       " 0.386547  0.306095\n",
       " 0.910741  0.824263\n",
       " 0.306095  0.545332\n",
       " 0.824263  0.954994\n",
       " 0.545332  0.165792\n",
       " ⋮                 \n",
       " 0.446527  0.221204\n",
       " 0.941787  0.649266\n",
       " 0.221204  0.877588\n",
       " 0.649266  0.400512\n",
       " 0.877588  0.924354\n",
       " 0.400512  0.269831\n",
       " 0.924354  0.764945\n",
       " 0.269831  0.681131\n",
       " 0.764945  0.841502\n",
       " 0.681131  0.499625\n",
       " 0.841502  0.962963\n",
       " 0.499625  0.137614\n"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "r1 = reconstruct(x, D, τ)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "999×2 Array{Float64,2}:\n",
       " 0.085     0.76827 \n",
       " 0.285813  0.681871\n",
       " 0.76827   0.837347\n",
       " 0.681871  0.51969 \n",
       " 0.837347  0.966676\n",
       " 0.51969   0.112748\n",
       " 0.966676  0.386547\n",
       " 0.112748  0.910741\n",
       " 0.386547  0.306095\n",
       " 0.910741  0.824263\n",
       " 0.306095  0.545332\n",
       " 0.824263  0.954994\n",
       " 0.545332  0.165792\n",
       " ⋮                 \n",
       " 0.446527  0.221204\n",
       " 0.941787  0.649266\n",
       " 0.221204  0.877588\n",
       " 0.649266  0.400512\n",
       " 0.877588  0.924354\n",
       " 0.400512  0.269831\n",
       " 0.924354  0.764945\n",
       " 0.269831  0.681131\n",
       " 0.764945  0.841502\n",
       " 0.681131  0.499625\n",
       " 0.841502  0.962963\n",
       " 0.499625  0.137614"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "r2 = RecurrenceAnalysis.embed(x, D+1, τ)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "true"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Matrix(r1) == r2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "scrolled": true
   },
   "source": [
    "## 1.2. Return type\n",
    "`reconstruct` returns `Dataset` which is a Vector of `SVectors`. I have found this to be more performant than a matrix for operations that truly understand the dataset as a vector of points. For example entropy computation.\n",
    "\n",
    "RA returns `Matrix`. More comparisons will follow later in this notebook.\n",
    "\n",
    "## 1.3. Benchmarks\n",
    "\n",
    "Now benchmark the functions for various dimensions and lengths"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "N = 1000, D = 1\n",
      "  4.206 μs (6 allocations: 15.91 KiB)\n",
      "  24.320 μs (3 allocations: 31.52 KiB)\n",
      "N = 1000, D = 10\n",
      "  16.213 μs (8 allocations: 84.94 KiB)\n",
      "  131.413 μs (5 allocations: 168.80 KiB)\n",
      "N = 100000, D = 1\n",
      "  436.480 μs (7 allocations: 1.53 MiB)\n",
      "  3.331 ms (5 allocations: 3.05 MiB)\n",
      "N = 100000, D = 10\n",
      "  2.660 ms (8 allocations: 8.39 MiB)\n",
      "  18.878 ms (5 allocations: 16.78 MiB)\n"
     ]
    }
   ],
   "source": [
    "for N in [1000, 100000]\n",
    "    x = trajectory(towel, N)[:, 1]\n",
    "    for D in [1, 10]\n",
    "        println(\"N = $(N), D = $(D)\")\n",
    "        @btime reconstruct($x, $D, $τ);\n",
    "        @btime RecurrenceAnalysis.embed($x, $(D+1), $τ);\n",
    "    end\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "RA is much slower than the method in DS. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1.4. Feature comparison\n",
    "\n",
    "DS also supports:\n",
    "\n",
    "- multiple timeseries embedding\n",
    "- embeddings with different delay times (aka delay vectors)\n",
    "- provides an embedding struct that can embed on demand any entry of the \"full embedding\"\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 2. Distance matrix\n",
    "\n",
    "## 2.1. Implementation\n",
    "\n",
    "The implementation in RA uses `pairwise` from `Distances` _and in addition_ requires all inpute to be transposed. This introduces at least two problems:\n",
    "1. Transposition, that may or may not be costly\n",
    "2. Using `views` several times over, that even though optimized, do have some cost.\n",
    "\n",
    "In addition, the entire RA package uses some string to identify the function that does the distance computation. I am sure that this leads to type instability. Although it becomes clear later that this is not a performance hit, it is still not Julia syntax, but Pythonic. I stand that users should directly use `Euclidean()` etc. as per multiple dispatch. It is \"bad education\" to use strings, I believe.\n",
    "\n",
    "---\n",
    "\n",
    "Let's try to create the same pairwise-like function that takes as an input a vector of vectors, something much more suited for this kind of operation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "pairwiseDS (generic function with 1 method)"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "VV = Vector{<:SVector}\n",
    "using Distances\n",
    "\n",
    "# main API\n",
    "distancematrixDS(x::Dataset, metric::Metric = Chebyshev()) = distancematrixDS(x, x, metric)\n",
    "function distancematrixDS(x::Dataset, y::Dataset, metric::Metric = Chebyshev())\n",
    "    pairwiseDS(x.data, y.data, metric)\n",
    "end\n",
    "# Core function: pairwise of vectors of svectors\n",
    "function pairwiseDS(x::VV, y::VV, metric::Metric)\n",
    "    d = zeros(eltype(eltype(x)), length(x), length(y))\n",
    "    for j in 1:length(y)\n",
    "        for i in 1:length(x)\n",
    "            @inbounds d[i,j] = evaluate(metric, x[i], y[j])\n",
    "        end\n",
    "    end\n",
    "    return d\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "999×999 Array{Float64,2}:\n",
       " 0.0       0.21861    0.686753    …  0.653867   0.781154    0.754746 \n",
       " 0.21861   0.0        0.506891       0.435305   0.622739    0.584749 \n",
       " 0.686753  0.506891   0.0            0.348783   0.145404    0.749531 \n",
       " 0.646566  0.427978   0.329197       0.0200783  0.471141    0.423315 \n",
       " 0.778068  0.620728   0.14662        0.492483   0.00557238  0.895209 \n",
       " 0.786553  0.615305   0.766052    …  0.41921    0.909082    0.0319517\n",
       " 0.960762  0.742153   0.492529       0.307119   0.58985     0.529249 \n",
       " 0.145147  0.286937   0.659619       0.701481   0.730623    0.864523 \n",
       " 0.551848  0.389043   0.654172       0.352467   0.799036    0.20291  \n",
       " 0.827637  0.640946   0.14307        0.397632   0.155021    0.800315 \n",
       " 0.313982  0.138038   0.546698    …  0.377811   0.679026    0.451318 \n",
       " 0.76248   0.603759   0.130292       0.477334   0.0189918   0.879488 \n",
       " 0.758212  0.577657   0.707592       0.360397   0.85041     0.0536949\n",
       " ⋮                                ⋱                                  \n",
       " 0.655731  0.487897   0.695091       0.364084   0.840364    0.0990287\n",
       " 0.865012  0.656784   0.255895       0.300556   0.329337    0.676237 \n",
       " 0.174648  0.206106   0.548545       0.595306   0.626146    0.79062  \n",
       " 0.67353   0.459632   0.452754    …  0.104109   0.594395    0.302503 \n",
       " 0.807811  0.639528   0.139716       0.467964   0.0528473   0.872821 \n",
       " 0.589906  0.427706   0.676254       0.362701   0.821525    0.165242 \n",
       " 0.839361  0.643923   0.172059       0.359933   0.214652    0.757588 \n",
       " 0.204343  0.0159982  0.522345       0.449569   0.637367    0.590098 \n",
       " 0.683877  0.505025   0.00532237  …  0.352001   0.143575    0.752232 \n",
       " 0.653867  0.435305   0.348783       0.0        0.490307    0.404965 \n",
       " 0.781154  0.622739   0.145404       0.490307   0.0         0.893354 \n",
       " 0.754746  0.584749   0.749531       0.404965   0.893354    0.0      "
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "distancematrix(r2, \"euclidean\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "999×999 Array{Float64,2}:\n",
       " 0.0       0.21861    0.686753    …  0.653867   0.781154    0.754746 \n",
       " 0.21861   0.0        0.506891       0.435305   0.622739    0.584749 \n",
       " 0.686753  0.506891   0.0            0.348783   0.145404    0.749531 \n",
       " 0.646566  0.427978   0.329197       0.0200783  0.471141    0.423315 \n",
       " 0.778068  0.620728   0.14662        0.492483   0.00557238  0.895209 \n",
       " 0.786553  0.615305   0.766052    …  0.41921    0.909082    0.0319517\n",
       " 0.960762  0.742153   0.492529       0.307119   0.58985     0.529249 \n",
       " 0.145147  0.286937   0.659619       0.701481   0.730623    0.864523 \n",
       " 0.551848  0.389043   0.654172       0.352467   0.799036    0.20291  \n",
       " 0.827637  0.640946   0.14307        0.397632   0.155021    0.800315 \n",
       " 0.313982  0.138038   0.546698    …  0.377811   0.679026    0.451318 \n",
       " 0.76248   0.603759   0.130292       0.477334   0.0189918   0.879488 \n",
       " 0.758212  0.577657   0.707592       0.360397   0.85041     0.0536949\n",
       " ⋮                                ⋱                                  \n",
       " 0.655731  0.487897   0.695091       0.364084   0.840364    0.0990287\n",
       " 0.865012  0.656784   0.255895       0.300556   0.329337    0.676237 \n",
       " 0.174648  0.206106   0.548545       0.595306   0.626146    0.79062  \n",
       " 0.67353   0.459632   0.452754    …  0.104109   0.594395    0.302503 \n",
       " 0.807811  0.639528   0.139716       0.467964   0.0528473   0.872821 \n",
       " 0.589906  0.427706   0.676254       0.362701   0.821525    0.165242 \n",
       " 0.839361  0.643923   0.172059       0.359933   0.214652    0.757588 \n",
       " 0.204343  0.0159982  0.522345       0.449569   0.637367    0.590098 \n",
       " 0.683877  0.505025   0.00532237  …  0.352001   0.143575    0.752232 \n",
       " 0.653867  0.435305   0.348783       0.0        0.490307    0.404965 \n",
       " 0.781154  0.622739   0.145404       0.490307   0.0         0.893354 \n",
       " 0.754746  0.584749   0.749531       0.404965   0.893354    0.0      "
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "distancematrixDS(r1, Euclidean())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2.2. Confirmation that it works"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "true"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "distancematrixDS(r1) == distancematrix(r2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2.3. Benchmarks"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "N = 1000, D = 1\n",
      "  4.694 ms (2 allocations: 7.61 MiB)\n",
      "  6.247 ms (28 allocations: 7.62 MiB)\n",
      "N = 1000, D = 5\n",
      "  14.265 ms (2 allocations: 7.49 MiB)\n",
      "  11.362 ms (28 allocations: 7.49 MiB)\n",
      "N = 1000, D = 10\n",
      "  34.079 ms (2 allocations: 7.34 MiB)\n",
      "  18.724 ms (28 allocations: 7.34 MiB)\n"
     ]
    }
   ],
   "source": [
    "for N in [1000]\n",
    "    x = trajectory(towel, N)[:, 1]\n",
    "    for D in [1, 5, 10]\n",
    "        println(\"N = $(N), D = $(D)\")\n",
    "        r1 = reconstruct(x, D, τ);\n",
    "        r2 = RecurrenceAnalysis.embed(x, (D+1), τ);\n",
    "        @btime distancematrixDS($r1)\n",
    "        @btime RecurrenceAnalysis.distancematrix($r2)\n",
    "    end\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "See also the discussion here: https://discourse.julialang.org/t/pairwise-distance-using-matrix-or-vector-svector-the-first-is-faster/17312/4 and there is also some kind of weird optimization happening with the tranpose... I am not sure though. \n",
    "\n",
    "For very small dimensionalitities SVector is faster, but then Matrix wins. It is very weird though, because the actual distance evaluation is faster for SVector:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "m = Chebyshev()\n",
    "v = rand(10)\n",
    "sv = SVector{10}(v);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  30.293 ns (0 allocations: 0 bytes)\n",
      "  25.173 ns (0 allocations: 0 bytes)\n"
     ]
    }
   ],
   "source": [
    "@btime evaluate($m, $v, $v)\n",
    "@btime evaluate($m, $sv, $sv);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The benchmarks are better if one uses Euclidean distance though!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "N = 1000, D = 1\n",
      "  4.498 ms (2 allocations: 7.61 MiB)\n",
      "  6.294 ms (2027 allocations: 7.69 MiB)\n",
      "N = 1000, D = 10\n",
      "  7.055 ms (2 allocations: 7.34 MiB)\n",
      "  6.087 ms (1991 allocations: 7.41 MiB)\n",
      "N = 1000, D = 100\n",
      "  61.214 ms (2 allocations: 4.90 MiB)\n",
      "  5.571 ms (1631 allocations: 4.95 MiB)\n"
     ]
    }
   ],
   "source": [
    "for N in [1000]\n",
    "    x = trajectory(towel, N)[:, 1]\n",
    "    for D in [1, 10, 100]\n",
    "        println(\"N = $(N), D = $(D)\")\n",
    "        r1 = reconstruct(x, D, τ);\n",
    "        r2 = RecurrenceAnalysis.embed(x, (D+1), τ);\n",
    "        @btime distancematrixDS($r1, Euclidean())\n",
    "        @btime RecurrenceAnalysis.distancematrix($r2, \"euclidean\")\n",
    "    end\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As you can see the version with `D = 100` is _much_ slower with static vectors. That is expected as these are the limits where static vectors are usable. We can have an if clause in the top level and if D near 100 the data is converted to matrix and then the `Distances.pairwise` is used."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 3. Recurrence Matrix\n",
    "\n",
    "The implementation of the Recurrence matrix is as fast as it can get. BUT. It's return value is not very reasonable. It returns a sparse matrix with values the booleans... My question is, why? We know what the booleans mean.\n",
    "\n",
    "Should it return just a vector of indices, while only keeping the indices that correspond to `true` ? Something like:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "recDS (generic function with 2 methods)"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function recDS(x, radius, metric = \"max\")\n",
    "    r = RecurrenceAnalysis.distancematrix(x, metric)\n",
    "    out = Tuple{Int, Int}[]\n",
    "    s = size(x)[1]\n",
    "    @inbounds for j in 1:s\n",
    "        for i in 1:s\n",
    "            # i == j && continue  <- This is worth considering...\n",
    "            r[i, j] <= radius && push!(out, (i, j))\n",
    "        end\n",
    "    end\n",
    "    return out\n",
    "end        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = trajectory(towel, 10000)[:, 1]\n",
    "data = Matrix(reconstruct(x, 1, 1));"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "10000×10000 SparseArrays.SparseMatrixCSC{Bool,Int64} with 33506 stored entries:\n",
       "  [1    ,     1]  =  true\n",
       "  [2    ,     2]  =  true\n",
       "  [3    ,     3]  =  true\n",
       "  [142  ,     3]  =  true\n",
       "  [4142 ,     3]  =  true\n",
       "  [4    ,     4]  =  true\n",
       "  [1273 ,     4]  =  true\n",
       "  [4143 ,     4]  =  true\n",
       "  [7390 ,     4]  =  true\n",
       "  [9999 ,     4]  =  true\n",
       "  [5    ,     5]  =  true\n",
       "  [1956 ,     5]  =  true\n",
       "  ⋮\n",
       "  [5328 ,  9997]  =  true\n",
       "  [5994 ,  9997]  =  true\n",
       "  [6071 ,  9997]  =  true\n",
       "  [9997 ,  9997]  =  true\n",
       "  [2068 ,  9998]  =  true\n",
       "  [9998 ,  9998]  =  true\n",
       "  [4    ,  9999]  =  true\n",
       "  [2069 ,  9999]  =  true\n",
       "  [7390 ,  9999]  =  true\n",
       "  [9999 ,  9999]  =  true\n",
       "  [4725 , 10000]  =  true\n",
       "  [10000, 10000]  =  true"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "recurrencematrix(data, 0.001)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "33506-element Array{Tuple{Int64,Int64},1}:\n",
       " (1, 1)        \n",
       " (2, 2)        \n",
       " (3, 3)        \n",
       " (142, 3)      \n",
       " (4142, 3)     \n",
       " (4, 4)        \n",
       " (1273, 4)     \n",
       " (4143, 4)     \n",
       " (7390, 4)     \n",
       " (9999, 4)     \n",
       " (5, 5)        \n",
       " (1956, 5)     \n",
       " (6, 6)        \n",
       " ⋮             \n",
       " (5328, 9997)  \n",
       " (5994, 9997)  \n",
       " (6071, 9997)  \n",
       " (9997, 9997)  \n",
       " (2068, 9998)  \n",
       " (9998, 9998)  \n",
       " (4, 9999)     \n",
       " (2069, 9999)  \n",
       " (7390, 9999)  \n",
       " (9999, 9999)  \n",
       " (4725, 10000) \n",
       " (10000, 10000)"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "recDS(data, 0.001)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "But once again my method is actually muuuch slower:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  1.956113 seconds (75 allocations: 777.565 MiB, 1.23% gc time)\n",
      "  9.819708 seconds (389.78 M allocations: 6.555 GiB, 7.44% gc time)\n"
     ]
    }
   ],
   "source": [
    "@time recurrencematrix(data, 0.001);\n",
    "@time recDS(data, 0.001);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "I guess I am a complete noob in LinearAlgebra optimizations.... :(\n",
    "\n",
    "Anyway, same discussion as above goes to `crossrecurrence` and `joinrecurrence`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 4. Estimating embedding dimension\n",
    "## 4.1. Features\n",
    "\n",
    "DS has only Cao's method while RA has also the FNN and the even more recent False First Nearest Neighbors. Very cool!\n",
    "\n",
    "## 4.2. Same output & Benchmark\n",
    "I'll compare Cao's method as its the only one we got"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  4.403417 seconds (9.36 M allocations: 468.482 MiB, 4.21% gc time)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "5-element Array{Float64,1}:\n",
       " 0.4225089613465132\n",
       " 0.8484447761147434\n",
       " 0.9113140842398368\n",
       " 0.924133227174861 \n",
       " 0.9594589837340649"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "using ChaosTools, PyPlot\n",
    "\n",
    "x = trajectory(towel, 10000)[:, 1]\n",
    "Ds = 1:5\n",
    "\n",
    "@time Es = estimate_dimension(x, 1, Ds) # compute E1 for Dimensions Ds .+ 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "([0.433097, 0.849603], [1.2272, 1.32898])"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ERA = afnn(x, Ds .+ 1, 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Alright the return value of the above is not really useful... We need E1 in many various dimensions to actually estimate when it saturates. Just two values cannot let us deduce this. Lets try again:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " 54.894057 seconds (593.86 k allocations: 56.627 GiB, 18.10% gc time)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "4-element Array{Float64,1}:\n",
       " 0.43309679708094556\n",
       " 0.8496026271703787 \n",
       " 0.9199791355071634 \n",
       " 0.9204688035749184 "
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "@time ERA = [afnn(x, (Ds[i], Ds[i+1]) .+ 1, 1)[1][1] for i in 1:4]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Okay, `afnn` is much, _much_ slower than our implementation. Several orders of magnitude."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiMAAAGUCAYAAADnDSLhAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvqOYd8AAAIABJREFUeJzs3Xl4VFW6/v1vZR4rgSRkYggzBAgEFIUEogwyqo0gtIjYaL8OnNOnvWht+teelrZf2z7n9fg7dKsgCiIKioCAgMEgqIDKjChDAAkzhCFA5rGq9vtHQSBWgFQg2Rnuz3VxWVVr1+aJgdTN2ms/y2IYhoGIiIiISTzMLkBEREQaN4URERERMZXCiIiIiJhKYURERERMpTAiIiIiplIYEREREVMpjIiIiIipFEZERETEVAojIiIiYiqFERERETGVwoiIiIiYysvsAmqCw+Hg9OnTBAcHY7FYzC5HRESkUTAMg7y8PGJiYvDwqPp8R4MMI6dPn6ZFixZmlyEiItIonThxgubNm1f5+AYZRoKDgwHn/wyr1WpyNSIiIo1Dbm4uLVq0KP8crqoGGUauXJqxWq0KIyIiIrXM3SUSWsAqIiIiplIYEREREVMpjIiIiIipFEZERETEVAojIiIiYiqFERERETGVwoiIiIiYSmFERERETNUgm55Vl8PhICMjg8OHD2MYhva1qSeufK/atGlD27Zt3doPQUREzKcwcplhGGzcuJGwsDAGDBiAt7e32SWJG8rKyjhw4AAbN26kf//+CpIiIvWI/gl52YkTJwgICKBr164KIvWQt7c3Xbt2JSAggJMnT5pdjohInXepoJQTFwvNLgPQzEi5EydO0LVrV7PLkFvUoUMH9uzZo12bRUQuMwyDk5eK2Hs6l32nc9iXmcve07lk5hRzX3wk70y8w+wSFUauyMvL06Z6DYDVaiUvL8/sMkRETFFmd3DoXD77TjsDx97L4SOv2Fbp8fsyc2u5wsopjFxD6wzqP30PRaSxKCixkX55lmPf6Vz2ZuZw8Ew+pXbHTd8b7OdFfLSVLjEh2B0Gnh7m/uxUGLlMH2INh76XItLQnMsrLp/t2JfpDB9HLxRgGDd/b5TVjy4xVrrEWImPcQaQ5k3869TPSoURERGROsLhMDh2sfBy8Li6vuN8XslN3+thgTYRQZdnPJzBIz7aSliQby1UfmsURkRERExQYrPz89l8Z+i4POuRnplLQan9pu/19fKgU7S1QvDoFBVMgE/9/Fivn1WLiIjUIzlFZRXXd5zO4dC5fGyOm19nCQ3wvnyZJaQ8fLQOD8TLs+F051AYERERuU0Mw+BMbrHL3SwnLhZV6f3Nm/iXLyyNv7zOIzrEr06t76gJCiPilg0bNpCWlsYbb7xR4Rba0NBQ/Pz8sNlsxMbGkpiYyMiRI/nVr36Fp6eny3nKysp47bXX+Pjjj8v/krVq1YoHH3yQ0tJSrFYrEyZMqLWvS0TEXXaHwZGs/GtmO5yLSy8WlN70vZ4eFto3c67viL+ysDQ6hJCAxtl0U2FE3NK/f3/69++Pr68v06ZN484772T16tWEhYUBcPHiRZYtW8bf/vY33n//fbp168aiRYvo1KlT+TlsNhsjRowgOzub1NTU8gZlGzduZPLkyezZs4cPP/zQlK9PRKQyRaV2DpzNq7C+Y/+ZXIrLbn4bbYCPJ52vrO24POvRPjIIP2/Xf6g1VgojUi3NmzcHnE3GrgQRgKZNm/Lkk08yevRo7r//fr799ltSUlLYsmULcXFxAMyZM4cvv/zSpVNqv379WL9+PYmJibX6tYiIXOtSQenlu1hyymc9Ms7nU4XlHYQH+RAfE3JN8LDSKizQ9D4edZ3CiFTLzXbGDQ0NZdmyZSQkJJCZmcmTTz7JunXrAFi9ejUAMTExLu9r2rQpf/jDH25/wSIiv1ChTXrm5Vbpp3M5nVNcpffHhQWU9+24EjyaWf1quOqGSWFEakx4eDhTp07lueee46uvvuK7774jKSkJh8M5rTlt2jT+9a9/ubzvvvvuY/v27bVdrog0YGV2Bxnn89l7Krd81mPf6Vxyr9Mm/VrenhY6RAZfne2IDaFTVDDBfo1zfUdNUBiRGjV+/Hiee+45ABYvXkxSUhIDBw5k5cqVvPHGG5w4cYK3336byMjI8vd06tSpwhoTERF3FJTY2H/m8t0sl8PHgbN5lNqq0Cbd14vOMRXXd7RrFoSPV8O5jbYuUhipovvf+LZKHfDqqohgX1b+Lrn2f9+ICOLi4jh69Cg7d+4E4Nlnn+WTTz5h06ZNLF++nK+//pr//M//5He/+x2+vnW/U6CI1B3n80oqdCpNP53LETfapMf/Ini0aFq32qQ3FgojVXQ+r4QzuVW7jigVRUREcPToUc6ePQuAj48PaWlpPPnkkyxevJicnBxeeOEFZsyYweuvv86oUaNMrlhE6hqHw+D4xcLL6zuuLiw9V4V/JFos0CY8sELvjs7RVsLrQZv0xkJhpIoiguv3H1oz67+y2NXb++r11eDgYBYtWsSKFSt47rnnOHLkCEeOHOGhhx5i3LhxzJs3T7MkIo3UlTbp1+7Pkp6ZR37Jzdd3+Hp50CkqmPhrgkd9bpPeWOi7U0VmXOJoKM6dOwdAVFSUy9gDDzzAkCFDmD59On//+9/Jy8vjk08+wdvbW71GRBqB3OIy9v2iadihc3mU2W9+nSXE39tlN9o2DaxNemOhMCI16sKFCxw9ehSAvn37VnqMr68vU6dOZcKECYwZM4bNmzczf/58XnrpJdq3b1+L1YpITTEMg7O5JRWahu3NzKlym/TYUP/ymY4rl1tiGkGb9MZCYURq1KJFizAMA4vFwvjx4wH485//zKuvvupybGxsLKtWraJDhw5cvHiRH374QWFEpB5ytkkvKA8eVxaXVrVNeruIoPLZjvjLi0tDA3xqoXIxi8KI1JicnBz++7//G4AJEyaU36577tw5tm/fzh133OHynrCwMAYMGMCSJUsIDQ2t1XpFpGpsdgc5RWVkF5WRU1RGTmEZmTnF5es79mfmUVRmv+l5/L096RwdXGFhaYfIYLVJb4QURqRabLYbLyQrKiri17/+NceOHaNr16688cYbFcanTZvGypUrXTq5GoZBeno60dHR9O/f/7bXLSJOhmFQUGp3horCUnIKy8oDRvblxzlFpZfHr32trEoLSX/pSpv0K51K42OsxKlNulymMCLVcvz4cQBKS0ux2+3lO/PabDbWrFnDn/70J3bv3s2oUaOYO3cuISEhFd6fmprKyJEjee211+jSpQsAhYWFTJ06lcOHD7Ny5Ur8/NRWWUx2pVmFYQC3+Li88cX1HlOFY1wflzkc5BWVkVtURm5xGfnFpeQWlZFX7Hwtv7iMvGIbuUWl5BeXkVtiI//y63bDwILz1xVXHlsqPDawAMGAFQOL5eoxXHNMAf4cMaJpFRZQoXdHfIyVZsG+Wt8h16UwIm7ZsGEDq1ev5p///Cfg3Gk3KCiIVq1a4efnR35+PqGhoQwePJg5c+Zw5513upwjLCyM1atXk56ezvPPP8+5c+ew2+1kZ2fTp08ftm/fTnx8fG1/aWK2siLY/zns+gjO/ATG5W6Z1/2Q5zqv36bH9YQ30PTyL7fUwBKMnOhkLI8vx6o26eImi2FUpU9d/ZKbm0tISAg5OTlYrdYqvSctLY0hQ4bUcGVSG/S9rEcMA07tgB/mw56lUJJjdkVyK9rcAxM/M7sKMVF1Pn9BMyMiYoa8M/DjQucsSNYB1/HACPAJ5MqFACyW2/SYmx5jWCzYDbDbDewG2BwGdoeBzeF8bPvF4zL71ed2x9V5lWsvfhjlFzquPr/RMcY1x1DJMR4W8PL0wMvTEy9PD3w8PfDycj738fTA29MDby9PvD0tl//riY+X8/Wr67Qsl7/uyh5ThWMqeRzeAZHqUBgRkdphK4GDX8APC+DQWjB+cbeFdwDE/woSH4WWfcHj1hpXldjs5Xd65FyzKDO7qIycwtJKFmte/WV31M6EcbCfFyH+3oQGeBPq70OIvzchAd6E+nuXvx7i702Iv0/549AAb/y9PbX+QhoUhRERqVmZPzoDyO7FUHTRdbxlX2cAiX8QfIMrDDkcBnklNnIrhInSCgEiu7C0Qti48rgqt5beDj6eHpUECJ+rIaM8UFwZ9yHU35tgPy91ChW5rNph5NChQ7z00kt8/fXX2Gw2BgwYwKuvvkrbtm2rdb7CwkJee+01PvnkE06dOoW/vz9Dhgzh5ZdfJi4urrpliogZCrLgp0XOyzBnd7uOW5tDj0eg+yMQ5vyZUVRq55vdmazancm+07nlIaOWJimw+nkRGnA1RFj9nQGjfEbC38f5WkDF1/y8PTRLIXKLqhVG1qxZw+jRoxk3bhwHDhzA09OT5557jp49e7JmzRruuusut853/vx57rnnHqxWK5999hlt27Zl/fr1TJw4kWXLlpGamkpysvaGEanT7GXOyy8/zIeDaeAoqzju5QedRjpnQVqngIcnxWV21u89w+c/ZbI2/SyFpbc2m+Hj5VEhQFS4vHFtyLgSOi6/FuznrX4XIiZy+26ajIwMevToQbt27dixY0f5YiibzUaXLl24dOkS6enphIWFVel8hmFwzz33sGfPHg4ePFjhfdu2baNPnz40adKE3bt3V7rRWmV0N03jpu9lLTuX7gwgPy2CgnOu47F3OANIl4fAP5QSm52NB7P4fHcmX+47W2kDrQAfTyKCfX9xeaPiuoqrYeLqbIY6d4qYq9buppkyZQr5+flMnjy5QvdMLy8vnnnmGaZMmcLUqVOZPXt2lc736aefsmHDBp566imXAHPnnXcyZMgQUlNTefHFF5kzZ4675YpITSi6BLuXOC/DnN7pOh4UCd1/Dd3HQ7NOlNocfJeRxaoff2TNvjPkFbsGkBB/b4Z2iWJEQjR92obhrfUUIo2GW2HkyJEjrFixAoBBgwa5jF/51+j8+fP57//+7yrNjrz//vsAdOvWrdLxsWPHkpqayvz583n99de1X4mIWRx2OPy1czHq/s/BXlJx3MMbOg6DxAnQdiBleLAp4wKfb/iJL/aeIaeozOWUwX5e3Bcfxcju0SS1DcfHSwFEpDFyK4ykpqYCEBgYSOvWrV3GO3bsiJ+fH8XFxSxfvpwnn3zypufctm0bAP7+/pWO33333YCz7fjGjRu5//773SlZRG5V1iHYtcDZFyTvtOt4VIIzgHR7GJtvKFuOXGTVZ+l8sSeTS4WuASTI14vB8ZGM6BZNvw7h+Hrp0opIY+dWGElLSwOgRYsWlY57enoSGxtLRkYGW7durVIYuXjReatfTk7lnRebN29e/vjgwYPulCsi1VWcC3uXOS/DnNjsOh4QBgnjoMd47M26su3oRVatOcUXe3aQle+6TXyAjyeDOkcyIiGalA4RWtshIhW4FUYOHz4MVAwIvxQWFkZGRgb79++v0jlDQkK4cOHCdY/39r66x0F2dnalx2RmZpKZmVn+PD8/H4Bdu3YRFBRU/np0dDTR0dFVqkuk0XE44OhGZwDZ9xnYiiqOWzyh/X2Q+CiOdvex41QBn2/NJHX3Os7llbiczs/bg4GdIhmZEM09HZvh76MAIiKVcyuMnD9/HoDg4ODrHuPr6wvApUuXqnTO3r17s3r1alavXo3D4XDZUv7a8wQGBlZ6jlmzZvHyyy+7vJ6SklLh+bRp0/jrX/9apbpEGo1LR2HXx/DjR5B93HU8ojMkPorRbSw7L/rw+U+ZpC77ljO5xS6H+np5cG/HZozsHs2ATs0I8FFfRRG5Obd+Uly5pBIQEHDdYxwO506bxcWuP6gqM3nyZFavXs3Jkyd57733+O1vf1th/MiRI+WPr9f87Omnn+aBBx4of56fn09KSgrr1693mRkREaC0APatcK4FObrRddwvBLo9jNHjUX6yt2bV7kxS39rLqewil0N9PD1I6RjByIRoBnaOJMhXAURE3OPWTw0fHx9sNtdb8q5VVuZcsNakSZMqnXPkyJH8/ve/55///Ce///3vCQwM5OGHHyY3N5eFCxfyt7/9rfzY6zVT++Xll9zcXAB69Ojh1n3OIg2aYcDxzc4Asnc5lOZVHLd4QNsBGN3Hs8+azMp9l/h8wWlOXPze5VTenhb6t49gREI0g+IjtWW8iNwSt8JIWFgYhYWFFBW5/uvoiisLUcPDw6t83unTp9OrVy/efPNNnnnmGf74xz/Sq1cvxo0bR2JiIl988QWJiYmV3sEjIjeRcwp+/Ni5FuRihut407YYPR7l5+iRfHbY4PMvMjl6YbvLYV4eFpLbhzOiWzT3xUcREqAAIiK3h1thpHPnzpw4cYKzZ89e95gLFy4AuB0cHnvsMR577LEKr50/f56JEycC8B//8R9unU+kUSsrhv2rnLMgGV9zdWP7y3yCoesojrcYxZJzMazadobD5w+5nMbTw0LftmGMTHAGkCaBPrVTv4g0Km6FkeTkZNasWVNhHce1CgoKyMrKAmDw4MG3XNxf//pXbDYb3bt3Z8KECbd8PpEGzTCc3VB/WAB7lkBxJbfLt+7P2TajWVyYyGd7s/n5+3yg4myJhwXubhPGiIRohnaJIizIt3bqF5FGy60wMnr0aF566SVOnz7NmTNnXPaK2b3buTunj48PAwYMuKXCvv/+e95++238/f15//338fLSori6YMOGDaSlpfHGG2+Ql+dcc+Dt7U3Tpk0pKioiJiaGpKQkXnjhBTp27Fjl8z777LNMnDiRPn361FTpDVfeWfjpE+dlmPPpruOhLbnU4WFWkMLHBy3s/zwPOFnhEIsFesc1ZWRCNEO6RtEs2K92ahcRwc0wEh8fz7Bhw8pvxZ00aVKF8XXr1gEwceLEG97+ezPHjx9n7NixeHh48OGHH9KjR49qn0tur/79+9O/f3/8/Px46aWXSE5O5uuvv8bLy4tLly7xt7/9jenTp/PRRx+xatWqKoXS7OxsPvzwQ/Ly8hRGqspWCge/cF6G+flLMH6x2613APltR7DWdyCzj8ewZ0M+kO9ymjtaNWFkQjTDukUTaVUAERFzuL0RxPTp0/H39+edd96p8HphYSFz5swhLCyMV155xeV9EyZMwGq1MmPGjBuef+fOnfTr14+cnByWLl3K6NGj3S1RakFsbCzg7CtzZdaqSZMm/O///i9jx46lqKiIxx9/vEq3eM+ZM4eCggIWL17MuXOV7PoqV2X+BKv/BP+3Eyx6zBlIrgkiJTG9+abjS4wLnkfXXaN4bouVPZkVQ0hiy1D+c0RnNv2fASx5ti+/SWqtICIipnI7jHTo0IG5c+eyfft2pk6dSklJCZmZmYwdO5bs7GxWrFhBZGRkhfdkZWWxYMEC8vLymDlzpss5c3NzWbduHZMmTeKuu+6iVatW7Ny5U/vQ1GG/bE53rSvbAJw8eZLNmytpJX4Nh8PBjBkzCAoKorS0lHffffe21tkgFFyAzW/D28kwqx9smQmFF8qH7UHR7Ip7kmebvEvHw8/xmx87seV0xT1hujcP4c/DO/Ht1HtZNjmJ3/ZrQ3RI5ftBiYjUtmotxBg3bhzNmzdn2rRpxMTEEBwczMiRI5k9e7bLOhJw3uY7YcIEli9fzrPPPlthbNKkSXz00UfExMTQt29fVq5cydChQ6v31Uid0LJly/LHVxY0X8+qVavw8vLiL3/5C1OnTmXWrFn86U9/wtOzkbcOt9vg0FrYNR8OfAGOiuHC8PTlSPi9zC9J5v0zcTiyXMNh11grI7rFMDIhmhZNr9+oUETEbNVeFZqUlMTatWurfPyHH35Y6etz585l7ty51S1D6qBDh67eItq5c+cbHvuvf/2LZ599lkcffZS//OUvnDhxghUrVjBq1KiaLrNuOrffuQ7kp08g3/UW+vMhXVlu3MOb5xPIORbkMt4pKpj7u8cwvFs0rcMr3z5BRKSu0S0qVTUrBfLr8XqGoGbw9Poa/23Kysp49dVXARg6dChdunS57rH79u1j8+bNLFmyhNDQUB566CEWLlzIm2++2bjCSFE27PnUGUJO7XAd9gljrfe9vHmpNwfOum5S2SEyiJEJzgDSrplrQBERqesURqoq/xzknTa7ijrLbrezZcsWXnzxRTZt2kSfPn2YN2/eDd/zr3/9i1//+teEhoYC8NRTT7Fw4UK++uor0tPTbzqrUq857HD4G2cASV8F9oq73totXuzwu5t3cvvwTXE3bL/4q9omIpCRCc5LMB0iq3/nmohIXaAwUlVBzcyu4NbUUP1bt26ld+/eZGRkcPHiRZo3b05aWhqDBw/GYrFc933Z2dnMnz+fDRs2lL9277330qFDBw4ePMiMGTN44403aqRmU13IcPYD+fFjyD3lMnzMuw0fFCWz1NaXS0UV91WKCwtgZEIMIxKi6RQVfMP/vyIi9YnCSFXVwiWO+ighIYFvv/2WrVu3kpyczMmTJ8nOzr7pB+WcOXPo2rUrPXv2rPD6U089xfPPP8+8efP4xz/+UWHX5XqrJM+5Md2uBXB8k8twnoeVpWV9+cTWn33FcRXGWjT1dwaQbtF0ibEqgIhIg6QwIrdF7969+fvf/84f//hHfvvb35KYmEj79u0rPdbhcPDWW2/h6enJ3XffXWGstLQULy8v8vLy+OCDD5g8eXJtlH/7ORxw7DtnANn3GZQVVhi248F6Rw8+sfXnK0dPyq75qxgb6s+IhGhGJkTTLTZEAUREGjyFEbltnn/+eb766iu++OILxowZw+bNm/H3d+1lsXLlShwOBwcPHqy0zf9TTz3Fu+++y1tvvVX/wsilY/DjQmcIyT7mMnzIiGWRrT/L7Mmcp0n561FWv/IA0qNFqAKIiDQqCiNy21gsFj744AO6d+/OTz/9xL/927/x3nvvuRz3xhtvMGXKlOvuNzRlyhRmz57Nvn37+Prrr7n33ntruvRbU1oI6SudPUGObHAZzjUCWGHvw2J7Cj8abQFn0GgW7Mvwbs4A0rNlEzw8FEBEpHFSGJFqudLmvaSk4l0gERERzJ8/n8GDBzN37lwSEhJ47rnnysd37NjBtm3bWL58+XXP3alTJwYOHMjatWt5/fXX62YYMQw4sdUZQPYsg9K8CsMOw8K3jq4stqewxnEHJfgAEB7ky/BuUYzoFs2dcU0VQEREUBiRajAMg6+++gqAgwcPcuHCBcLCwsrHBwwYwJ///GdeeeUVpkyZwrlz5/jjH/+It7c3//7v/05MTEyll2+u1a5dO9auXcvnn3/Oxx9/zCOPPFKjX1OV5Z6+fBnmI7jws8vwEUckS+wpLLX3IxPn/5OmgT6M7hrFyIRo7modhqcCiIhIBRbDMAyzi7jdcnNzCQkJIScnB6vVevM3AGlpaQwZMqSGK6v/XnvtNaZPn87p01d7rgQHB9OlSxc2bbp6p4jdbufee+9l48aNAHh5eREYGEhOTg4AzZo1Y9q0aS5rQvLz8+nZsyc//1zxg759+/YcPHiwSjXe9u9lWTEcSHWuA8n4CgxHheF8w4/P7Xez2N6f7UZHwEJogDdDu0QxIiGaPm3C8PJ0exsoEZF6pzqfv6CZEXHTCy+8wAsvvHDT4zw9PSv0EKmqoKCgKoeOGmUYcPoHZwDZvQSKs10O2WSPZ7G9P6sdvSnCD6ufFw9fDiBJ7cLxVgAREakShRGRa+Wfg58WOUPIuX0uwyeNcD6192eJvR8njEiCfb0YFh/JyO7RJLeLwMdLAURExF0KIyK2Uvh5DexagPHzGiwOW4XhIsOH1Y7eLLansNnRmQAfbwZ1jeSlhBj6tQ/Hz7uR7zAsInKLFEak8TqzxxlAflqEpTALuHLTrdN2RwcW21NItd+FzTuYgV2bMTMhmns6NlMAERG5jRRGpHEpvAi7l2Dsmo8l80egYgA5YzRhqb0fS+z9Oe3VnAHxzfivbjEM6NQMfx8FEBGRmqAwIg2f3QYZX2Hsmo+xfzUejtIKAaTE8OJLxx0stqewxaM7/TtE8lz3GAZ2akagr/6KiIjUNP2klQbN2DYH29f/hXfhOSxUnAX50dGGxfYU0kiie4c4RiXE8GbnZgT7eZtVrohIo6QwIg1WyYG1+H4+hWujxXnDynJ7MkuNe4hql8iIhBheiI8kxF8BRETELAoj0mBlffk6sZcff2XvwULHQMraDGJYQgs+7hJJaICPqfWJiIiTwshlDbARbaNlGAaOM3uJzfoegOOOCI4PmcN/JbakaaACiIhIXaMOTdew2+1mlyC3yG63Y7FYOLPm/5a/9lXoaH6T3E5BRESkjlIYuSwmJoaTJ0+aXYbcopMnTxIdGkjE4c8AyDUCaDHg/zG5KhERuRGFkcs6duzInj17OHbsmGZI6iG73c6xY8fYs2cPEZlr8KYMgFXeQ7gnoa3J1YmIyI1ozchlvr6+DBo0iAMHDrB//34Mw8Bi0Vbv9cGV71V0dDSD+idhm/4sAGWGJx53P4Onh76PIiJ1mcLINXx9fUlISCAhIcHsUqSaCr57l0B7DgBp9GFEci+TKxIRkZvRZRppOBwOSr99o/zp6fjfqoGZiEg9oDAiDYbtwBc0KToGwCZHPEMHDTG5IhERqQqFEWkwLq2bXv54e/QjtAwLMLEaERGpKoURaRgyfyQiawsAGY5o7hj8iMkFiYhIVSmMSINwYe3/lj9ODRrF3W3DTaxGRETcoTAi9V/OKUIzVgBw0QgiNmWSbssWEalHFEak3svbOANPnI3qlnoMYXhPNTkTEalPFEakfivJx/uHec6Hhhf2O36Ln7enyUWJiIg7FEakXivd/gF+9jwAVhrJjOrf0+SKRETEXQojUn857JR8+1b50yPtHqdZsJ+JBYmISHUojEi95UhfRXCRc6flDfZuDBs40OSKRESkOhRGpN7K+/pqk7NvI8bRNTbExGpERKS6FEakfjqxjZCsnQAccDSn14AxJhckIiLVVe0wcujQIcaPH090dDQRERGMGzeOjIyMahdy5MgRnnrqKVq1akVQUBBNmzZl8ODBrFixotrnlIYr75ursyJL/R5kUHyUidWIiMitqFYYWbNmDYmJiQQEBHDgwAGOHj2K1WqlZ8+ebNmyxe3zbd68mR49erBnzx6WLVtGbm6DvcNuAAAgAElEQVQuR44cYeDAgfzqV79i6tSp1SlTGqpLxwjMSAXgvBFCVPLjeHqoyZmISH1lMQzDcOcNGRkZ9OjRg3bt2rFjxw48PJx5xmaz0aVLFy5dukR6ejphYWFVOl9paSkdOnSgqKiIAwcOEBoaWmH8t7/9LXPmzGHdunUMGDCgSufMzc0lJCSEnJwcrFarO1+e1APFq6bit/1tAN5wjOXxP8/A6udtclUiIlLdz1+3Z0amTJlCfn4+kydPLg8iAF5eXjzzzDOcP3/erZmMjRs3cuzYMZKTk12CCMCQIc5t4D/77DN3S5WGqDgHjx8+dD40vCnq/riCiIhIPedWGDly5Ej5Go5Bgwa5jF8JDvPnz+fChQtVOmdWVhYAp06dqnTcz8/ZNyIwMNCdUqWBsm+fh4+9AIBPHf0Zm5JockUiInKr3AojqanO6/SBgYG0bt3aZbxjx474+flRUlLC8uXLq3TOzp07A7Blyxa++uorl/HNmzfj4eHBo48+6k6p0hDZyyj9bkb50/SWE4gLV0gVEanv3AojaWlpALRo0aLScU9PT2JjYwHYunVrlc6ZkJBAUlISAOPGjWPnzp3lYydOnGDmzJn8/e9/p0uXLu6UKg3Rvs/wL8oEYK09keED+ptckIiI3A5uhZHDhw8D0Lx58+sec2Xh6v79+6t83g8++ICoqCiysrLo168fS5cu5fz584waNYrXXnuNP/3pT+6UKQ2RYVCw/p/lT9dYx9CnTdUWSYuISN3m5c7B58+fByA4OPi6x/j6+gJw6dKlKp+3TZs2fPPNN9x3330cP36cMWPGEBsby/z580lJSbnp+zMzM8nMzCx/np+fD8CuXbsICgoqfz06Opro6Ogq1yV1yPFNBGb9BMAeRxx33PMAFotu5xURaQjcmhm5ePEiAAEBAdc9xuFwAFBcXOxWIR07dmT+/Pm0b9+ewMBATp48yaOPPsqOHTtu+t5Zs2bRq1ev8l9XAkxKSkqF12fNmuVWTVJ3FG+4Oiuy0PMBHugRa2I1IiJyO7k1M+Lj44PNZrvhMWVlZQA0adLErULS0tJ4/vnn2bBhA+fOnWPkyJGcOHGC/v37s3z5cgYPHnzd9z799NM88MAD5c/z8/NJSUlh/fr1LjMjUg9dyMA3w7leKdNoSnifX+Pn7WlyUSIicru4FUbCwsIoLCykqKjousfk5OQAEB4eXuXzfv311zzwwAMsWbKEqKgooqKi+O677xg8eDAHDhxg1KhRbNu2rfzOm1/65eWX3NxcAHr06KGmZw2A7bu38MLZm2++YwiP921rckUiInI7uXWZ5koYOHv27HWPudJfpLJbfytTWlrK448/TmRkJCNHjix/vUWLFqxbt44WLVpQUFDA888/706p0lAUXoRdCwAoMHzJjn+UZsF+JhclIiK3k1thJDk5GXA2P6tMQUFBeROzG11WudaqVas4ceIEiYmJLgsSY2NjWbDA+UG0Zs2aG87ISMNkbH8PL4dz/dEi+z080j/B5IpEROR2cyuMjB49GoDTp09z5swZl/Hdu3cDzrUlVd1H5uDBg4CznXxl+vXrR/v27bHZbOWXgKSRsJVQ+r1zDxqHYWFH9K/pGhticlEiInK7uRVG4uPjGTZsGACrV692GV+3bh0AEydOvOHtv9dq1aoVAOnp6dc9xs/PjyZNmhAREeFOuVLf7VmKb7HzdvI0xx2MTOljckEiIlIT3N4ob/r06fj7+/POO+9UeL2wsJA5c+YQFhbGK6+84vK+CRMmYLVamTFjRoXXH3zwQWJjY0lPT+fzzz93eV96ejrp6ek888wzeHrqDopGwzAo2Xj1dt4V/qMYHB9lYkEiIlJT3A4jHTp0YO7cuWzfvp2pU6dSUlJCZmYmY8eOJTs7mxUrVhAZGVnhPVlZWSxYsIC8vDxmzpxZYSwgIIBPP/2UJk2a8Jvf/IZly5Zhs9koKSlhzZo1jBw5kqFDh/Lyyy/f2lcq9cuR9fhecM6W/eBoR6/koXh6qMmZiEhD5HYYAeceMt988w07duwgJiaGPn36EBcXx759++jbt6/L8eHh4UyYMIGgoCCeffZZl/G77rqLXbt28fDDD/OHP/yB0NBQoqOj+dvf/saf//xnPvvsM7y9tU18Y1L27b/KH3/ISMb2bmliNSIiUpMshmEYZhdxu+Xm5hISEkJOTo76jNRH5/bDjLsAOGmE817Ppbz0YHeTixIRkZup7udvtWZGRGqSY9Nb5Y/n2ofyWFI7E6sREZGapjAidUv+eYwfFwKQa/hztu3DtA4PNLkoERGpSQojUrdsm42noxSAhfYBjO/X1eSCRESkpimMSN1RVkTZlncBsBkebGz6EH3ahplclIiI1DSFEak7fvoE72Ln3kapjru4v99dLlsEiIhIw6MwInWDw4Htu6sLVz/xfpAHesSYWJCIiNQWhRGpGw6txeuic5+iLY5O9Lp7AH7e6rgrItIYKIxInWD//s3yx+87hjPh7lYmViMiIrVJYUTMd2Y3nkfXA3DEEYl/15E0s/qZXJSIiNQWhRExnbHp6qzIe/ZhTEpWkzMRkcZEYUTMlZuJsftTALKNQI7EPki35iEmFyUiIrVJYUTMtfUdPBxlACywD+TRfp1NLkhERGqbwoiYp7QA+7b3nA8NT9YEPsjg+EiTixIRkdqmMCLm2fURniXZAKx09GVkUk+8PPVHUkSksdFPfjGHw479+6tNzuZbRjL2zhYmFiQiImZRGBFzHFiNZ/YRAL61dyGhVzIh/t4mFyUiImZQGBFTXHs772z7CH6T1NrEakRExEwKI1L7Tu3AcnwTAD87YvFqP4jW4YEmFyUiImZRGJHat+nqWpHZ9uFM6tfWxGJERMRsCiNSu7JPYOxdDkCWYWVf2BD6tg0zuSgRETGTwojUri1vYzHsAMy3D2JCv45YLBaTixIRETN5mV2ANCLFuTh2zMMDKDG8WeE9nNQesWZXJSIiJtPMiNSeHz7EozQPgKX2ZIbfnYCft6fJRYmIiNkURqR22G0Ym2eWP53nGM5jfVqZWJCIiNQVCiNSO9JXYMk5AcDX9u50SriTSKufyUWJiEhdoDAiNc8wKjQ5e9c+gieS1eRMREScFEak5p3YguXUDgDSHS0pbZ5MQvNQk4sSEZG6QmFEat61syK24TzRr42JxYiISF2jMCI16+JhjPRVAJw1QtkRPID74iNNLkpEROoShRGpWZvfxoIBwDzbEB5NaoeXp/7YiYjIVWp6JjWn6BLGD/OxAIWGL8s87+OLO1qaXZWIiNQx+ieq1Jwd72MpKwBgsb0/g3t1IiTA2+SiRESkrlEYkZphK8XYMgsAh2Fhrn0ov+kbZ25NIiJSJymMSM3YuwxLXiYAax09adOxO20igkwuSkRE6iKFEbn9DOMXt/OO4IkkNTkTEZHKKYzI7Xd0I5z5CYAfHW3IiehFUrswk4sSEZG6SmFEbr9Nb5U/nG0bzqTkNlgsFhMLEhGRukxhRG6v8wfh4BcAnDLC2OybxKjEWJOLEhGRukxhRG6vzTPKH861DWXs3W3w8/Y0sSAREanrqh1GDh06xPjx44mOjiYiIoJx48aRkZFRrXP9/ve/x2Kx3PTXa6+9Vt1ypTYUZGH8+DEAeYY/S4wBPHZ3nLk1iYhInVetMLJmzRoSExMJCAjgwIEDHD16FKvVSs+ePdmyZYtb5yoqKuKDDz6o0rEjRoyoTrlSW7a/h8VWDMAn9ntISWhLVIifyUWJiEhd53Y7+IyMDEaPHk27du1455138PBw5pmZM2eyYcMG7r//ftLT0wkLq9rdEwsXLsTPz4+ZM2eSmJhISEiIyzEPP/wwdrud+Ph4d8uV2lJWjLH1HSyA3bDwvn0ob+p2XhERqQK3w8iUKVPIz89n8uTJ5UEEwMvLi2eeeYYpU6YwdepUZs+eXaXzrVq1ip07dxIdHV3p+OnTp9m7dy9/+ctf3C1VatPuxVgKzgOw2nEXkS070KNFqMlFiYhIfeDWZZojR46wYsUKAAYNGuQyPmTIEADmz5/PhQsXbnq+srIyXnzxxesGEYClS5diGAYPP/ywO6VKbTIMl9t51eRMRESqyq0wkpqaCkBgYCCtW7t+2HTs2BE/Pz9KSkpYvnz5Tc/n7e1Nz549b3jMkiVL6NSpE127dnWnVKlNGevgfDoA2xwdOGftypAukSYXJSIi9YVbYSQtLQ2AFi1aVDru6elJbKyzp8TWrVtvsTQ4d+4c3377LWPGjLnlc0kN+v5q6/fZtuE83jcOL0/dNS4iIlXj1pqRw4cPA9C8efPrHhMWFkZGRgb79++/tcqAZcuWYbfbb3qJJjMzk8zMzPLn+fn5AOzatYugoKubs0VHR9/wkpBUw9m9cPhrAI45mvGt5138f3e2NLkoERGpT9wKI+fPOxcoBgcHX/cYX19fAC5dunQLZTktWbKEDh06kJCQcMPjZs2axcsvv+zyekpKSoXn06ZN469//est1yXX2HS1ydl79mE8dGdLQgK8TSxIRETqG7fCyMWLFwEICAi47jEOhwOA4uLiWyjL+Xt98803TJ069abHPv300zzwwAPlz/Pz80lJSWH9+vUuMyNyG+Wdxdi9CAuQYwSw2J7CyqQ4s6sSEZF6xq0w4uPjg81mu+ExZWVlADRp0qT6VQHLly/HZrNVab3ILy+/5ObmAtCjRw+sVust1SE3sO1dLPZSAD6yD6R3xxa0jQi6yZtEREQqcmuV4ZVGZkVFRdc9JicnB4Dw8PBbKAs+/fRT2rVrR48ePW7pPFJDSgth2xwAygxP5tnu0+28IiJSLW6Fkc6dOwNw9uzZ6x5zpb9IZbf+VlVOTg5r165Vb5G67MePoch52W6low9BzVrRr/2tBVAREWmc3AojycnJgLP5WWUKCgrIysoCYPDgwdUuasWKFZSWliqM1FUOR4XdeedcbnJmsVhMLEpEROort8LI6NGjAWeL9jNnzriM7969G3CuLRkwYEC1i/r0009p06YNiYmJ1T6H1KCf0+DCIQC+t8dzyr89oxJjTS5KRETqK7fCSHx8PMOGDQNg9erVLuPr1q0DYOLEiTe8/fdG8vPzSUtL06xIXXZtkzP7cMb3bom/j6eJBYmISH3mdpvM6dOn4+/vzzvvvFPh9cLCQubMmUNYWBivvPKKy/smTJiA1WplxowZLmPX+vzzzykuLlYYqatO/wDHvgUgwxHNRhJ5rE8rk4sSEZH6zO0w0qFDB+bOncv27duZOnUqJSUlZGZmMnbsWLKzs1mxYgWRkRX3JcnKymLBggXk5eUxc+bMG55/yZIltG7dml69erlbmtSGazbEm2MfzrBusUSH+JtYkIiI1HfV2kBk3LhxfPPNN+zYsYOYmBj69OlDXFwc+/bto2/fvi7Hh4eHM2HCBIKCgnj22Weve96ioiJWr16tvWjqqpyTGHuXAXDRCOJTez+eSNbtvCIicmvcanp2raSkJNauXVvl4z/88MObHuPv71++r4zUQVtmYXE4m959aB9Ml5bN6NEi1OSiRESkvtPWqlI1JXmwY57zoeHFfNtgzYqIiMhtoTAiVfPDAihxdtddbk/GOySKoV2iTC5KREQaAoURuTmHvWKTM/swHusTh5en/viIiMit06eJ3Nz+VZB9DID19gSOe7Xikd4tTC5KREQaCoURublfNDkb3bM5oQE+JhYkIiINSbXvppFG4sRWOLkVgP2OFmx0dGNtUpy5NYmISIOimRG5sU1XZ0Xm2IeR0qEZ7ZpVr9W/iIhIZRRG5PouHYX0lQCcN0L4zJ6k23lFROS2UxiR69v8NhgOAObZ7qNlsyb0bx9uclEiItLQKIxI5Yqy4Qdn19wiw4cF9oFMSorDYrGYXJiIiDQ0CiNSuZ3zoNTZmv9Tez8c/mE8lNjc5KJERKQhUhgRV/Yy2DKr/Okc+3DG39USfx9PE4sSEZGGSmFEXO37DHJPAfClvSfHLTFM7NPK5KJERKShUhiRigwDvn+j/Okc+3CGd4smOsTfxKJERKQhU9MzqejY95C5C4Ddjjg2OzqzTE3ORESkBmlmRCq6psnZbNtwerRoQmLLJiYWJCIiDZ3CiFyVdQgOrAYg02jK54671eRMRERqnMKIXLV5BmAA8L5tCOHWIIZ1jTK3JhERafAURsSp8CLs+giAAsOXj+0DmNi3Fd6e+iMiIiI1S5804rR9DtiKAFhkv4dS72AeubOlyUWJiEhjoDAiYCuBre8CYDcsvGcfykM9m9Mk0MfkwkREpDFQGBHYvQTyzwKQ5riTE0Ykk/rGmVuTiIg0GgojjZ1hwKa3yp/Otg2nf4cI2kcGm1iUiIg0Jmp61tgd/hrO7QVgp6MdO40OvK8mZyIiUos0M9LY/WJWpG1EIP3bR5hYkIiINDYKI43ZuXQ4tBaAE44I0hx3MimpNR4eFpMLExGRxkRhpDG7ZlZkrn0oQf5+PNQz1sSCRESkMVIYaazyz8FPiwDINfz5xH4Pj/RuSYCPlhGJiEjtUhhprLbNBnsJAB/bB1DsEcDEPq1MLkpERBojhZHGqKzIGUYAm+HB+7ahDO0aRUyov8mFiYhIY6Qw0hj9uBAKLwDwueNuMgnjiSTtzisiIuZQGGlsHI7Lu/M6zbYNp3uLUHq2DDWxKBERacwURhqbQ19C1kEAtjg6sdtowxNJcVgsup1XRETMoTDS2Gx6s/zhu7YRRFp9Gd4t2sSCRESksVMYaUwyf4IjGwA47IhinSORiX3i8PbUHwMRETGPPoUak2uanL1nH4aPlxfje7c0sSARERGFkcYj9zTsWQLAJSOIT+39eKhnc5oE+phcmIiINHYKI43F1nfAYQNggX0gRfjxhHbnFRGROkBhpDEoyYft7wFQangyz3Yf/dqH0z4y2OTCREREbiGMHDp0iPHjxxMdHU1ERATjxo0jIyPjthWWn5/Pu+++y6hRoxgyZAhPPPEE27Ztu23nb1R2fQTFOQCscCRxniY8kawmZyIiUjdUK4ysWbOGxMREAgICOHDgAEePHsVqtdKzZ0+2bNlyy0V9/PHHtG3blsWLF/Pyyy+TlpbGe++9x5133nnL5250HHaXJmdtIgJJaR9hYlEiIiJXuR1GMjIyGD16NO3ateOdd97BarUSGBjIzJkziYqK4v777+fChQvVKsbhcPDv//7vjB8/nsmTJ5OWlkZCQkK1ziWXHUiFS0cA2Gjvyn6jJZOSWuPhoSZnIiJSN7gdRqZMmUJ+fj6TJ0/Gw+Pq2728vHjmmWc4f/48U6dOrVYxv/vd73jrrbeYNm0a06ZNU1fQ2+Ga23ln20dg9fNidM9YEwsSERGpyK0wcuTIEVasWAHAoEGDXMaHDBkCwPz5892eHfmf//kfZsyYwcCBA5k2bZpb75XrOLkDjm8C4KAjlvWOBB7p3ZIAHy+TCxMREbnKrTCSmpoKQGBgIK1buy6A7NixI35+fpSUlLB8+fIqn/fQoUO8+OKLWCwW3nzzTc2I3C7XtH6fbR+Op4cHE/vGmVePiIhIJdwKI2lpaQC0aNGi0nFPT09iY52XALZu3Vrl8z7//POUlpZy//3306lTJ3dKkuvJPg77PgPgvGHlM3sSQ7tEERvqb3JhIiIiFbk1X3/48GEAmjdvft1jwsLCyMjIYP/+/VU658mTJ8sv/QwePJhXX32V1NRU0tPTsVgs9O3bl5dffpnExMTrniMzM5PMzMzy5/n5+QDs2rWLoKCg8tejo6OJjm4km8JtmQWGHYAPbfdRgg9PJMeZW5OIiEgl3JoZOX/+PADBwddvluXr6wvApUuXqnTOpUuXYhgGAIsWLaJXr158+eWX7Nu3j4cffpiVK1fSp08fVq9efd1zzJo1i169epX/SklJASAlJaXC67NmzapSTfVecQ7smOd8aHgz3z6I7s1D6NmyicmFiYiIuHJrZuTixYsABAQEXPcYh8MBQHFxcZXOuXHjRgDatm3L+vXry9eL+Pv7M3PmTI4fP05qaiqPPvoohw4domnTpi7nePrpp3nggQfKn+fn55OSksL69etdZkYahZ0fQmkeAEvt/biIlWnJrbUWR0RE6iS3woiPjw82m+2Gx5SVlQHQpEnV/hV+6tQpwLkOpbIPyz/84Q+kpqZy6dIlPvjgA5577jmXY355+SU3NxeAHj16YLVaq1RHg2G3wZa3y5/OsQ8j0urLsK6NJIiJiEi949ZlmrCwMACKioque0xOjrPteHh4eJXOmZWVBVw/vNx7773lQWPTpk1VrrXRSv8Mck4A8JW9BxlGLBP7xOHjpW2IRESkbnLrE6pz584AnD179rrHXOkvUtmtv5UJDQ0Frl7e+SWLxUK7du2Aqq9DabQMA76/ejvvu/YR+Hp58EjvliYWJSIicmNuhZHk5GTA2fysMgUFBeUzHYMHD67SOdu0aQNcvVxTmZiYGODqzIxcx/HNcHonAPscrdjkiOehnrE0DfQxuTAREZHrcyuMjB49GoDTp09z5swZl/Hdu3cDzrUlAwYMqNI5hw4dCsCPP/5IYWFhpcdcmRG5++673Sm38bmmydm7tuGAhUlJ2p1XRETqNrfCSHx8PMOGDQOo9FbbdevWATBx4sQb3v57rTFjxhAZGUlZWRmffvqpy7hhGOzfvx9/f3/Gjx/vTrmNy4UM2P85AGeMJqxy9KFf+3A6RFbt+yAiImIWt1c1Tp8+HX9/f955550KrxcWFjJnzhzCwsJ45ZVXXN43YcIErFYrM2bMqPB6UFAQr732GgCvvvoqBQUFFcYXLlzI8ePH+etf/0pEhLa9v64tbwPOfi3zbEMow4tJSXGmliQiIlIVboeRDh06MHfuXLZv387UqVMpKSkhMzOTsWPHkp2dzYoVK4iMjKzwnqysLBYsWEBeXh4zZ850Oedjjz3Giy++yP79+xkzZgwnT57EMAxWrVrF5MmTefrpp/njH/9Y/a+yoSu8CD/Mdz40fPnIPoDW4YHc06GZyYWJiIjcXLXu9xw3bhzffPMNO3bsICYmhj59+hAXF8e+ffvo27evy/Hh4eFMmDCBoKAgnn322UrP+corr7Bw4UIuXbpE+/btiYqK4vXXX2f27Nm8/fbblb5HLtvxPpQ519sssqeQQxCTkuLw8FCTMxERqfssxpVe7A1Ibm4uISEh5OTkNPymZ7ZS+GcC5GXiMCzcW/o6F32bs/n/DCTQ162ediIiIrekup+/6oRV3+1dCnnOTQLXOO7gmBHFI71bKoiIiEi9oTBSn/2iydls2zA8LDCxTysTixIREXGPwkh9dmQDnHX2dtnlaMt2oyNDu0bRvMn1NzIUERGpaxRG6rNN186KOJucPaEmZyIiUs8ojNRX5w/Az2sAOGmEs9rRm4TmIfRqVbXdkkVEROoKhZH6avPV5nFzbUOw48kTSa2xWHQ7r4iI1C8KI/VRQRb8uBCAPMOfRfZ7aRbsy/Bu0SYXJiIi4j6Fkfpo2xywFQOw0H4veQQwsU8rfLz07RQRkfpHn171TVkxbHsXAJvhwfu2Ifh6efBI75YmFyYiIlI9CiP1ze5FUHAegNWO3pwigl/1iCUsyNfkwkRERKpHYaQ+MQzY9Fb5U+ftvDApOc6kgkRERG6dwkh9cmgdnN8PwFZHR3402pHULoxOUQ18/x0REWnQFEbqk01vlD+cc3lWRE3ORESkvlMYqS/O7IHD3wBw1BHJl45exIUFcG/HZubWJSIicosURuqLa9aKvGcfigMPJiW1xsNDTc5ERKR+UxipD/LOwO7FAGQbgSy2pxDs58WYXs1NLkxEROTWKYzUB1vfAUcZAB/ZB1KEH7++swWBvl4mFyYiInLrFEbqutIC2P4eADY8mWe7Dw8LTOwTZ25dIiIit4nCSF3348dQdAmAz+x9OEtThnSJokXTAJMLExERuT0URuoyhwM2Xd2dt/x23mTdzisiIg2HwkhddvALuJgBwHf2Luwz4ugWG8IdrZqYXJiIiMjtozBSl216s/zhbPuVWZE4LBbdzisiIg2HwkhddWonHPsOgEOOGL5xdCci2JcR3WJMLkxEROT2Uhipq65pcjbHPgwDDx67uxU+XvqWiYhIw6JPtroo5yTsXQbARYJZau+Hj5cH4+9qaXJhIiIit5/CSF205W0w7AB8aBtMCT78qkcM4UG+JhcmIiJy+ymM1DUlebBjHgClePOhbTAAk7Q7r4iINFAKI3XNzg+hJBeApbYksgihb9swOkdbTS5MRESkZiiM1CV2G2yZWf50zpXbeTUrIiIiDZjCSF2yfyVkHwfgG3t3fjaa0yosgAGdmplcmIiISM1RGKlLrrmd90qTs0l94/DwUJMzERFpuBRG6ooTW+HkNgD2Gy351tGVYF8vxtzRwuTCREREapbCSF3x/RvlD2fbhgEWxt3ZgiBfL/NqEhERqQUKI3XBxSOwfxUAWYSywt4XDws83jfO3LpERERqgcJIXbDlbTAcAMwtu49SvLkvPooWTQNMLkxERKTmKYyYrSjb2VsEKMGXBfaBAExKijOxKBERkdqjMGK2He9DWQEAi2z9yCaYLjFWerduam5dIiIitURhxEz2MtgyCwAHFt6zDwOcTc4sFt3OKyIijUO1w8ihQ4cYP3480dHRREREMG7cODIyMm6pmNLSUlq0aIHFYnH5NXr06Fs6d520dxnknQZgnb0nR4xowoN8Gdk92uTCREREak+1wsiaNWtITEwkICCAAwcOcPToUaxWKz179mTLli3VLmbevHmcPHmy0rEXX3yx2uetkwwDNr1Z/vRdm7PJ2WN3t8LXy9OsqkRERGqdxTAMw503ZGRk0KNHD9q1a8eOHTvw8HDmGZvNRpcuXbh06RLp6emEhYW5VYjdbqdTp0689dZbtGzZssKYl5cX7dq1q/K5cnNzCQkJIScnB6u1jm4wd2QjzBsJwF7aMKL4/8XH05Pv/88AwoN8TS5ORFULY8gAABdJSURBVETEfdX9/HW7o9aUKVPIz89n8uTJ5UEEnIHhmWeeYcqUKUydOpXZs2e7dd6FCxfSoUMH7rvvPndLqp+uaf0+q3Q4YOHBHjEKIiIi0ui4dZnmyJEjrFixAoBBgwa5jA8ZMgSA+fPnc+HChSqf1zAM/vGPfzBmzBjcnKipn7J+hoOrAThnCSPV0RuASdqdV0REGiG3wkhqaioAgYGBtG7t+sHZsWNH/Pz8KCkpYfny5VU+7/Lly///9u48Oqry/uP4ezJkgSwIEwiJSEKOQoOCWKIsglgQkFooyw9QjKgoa/WUH3iwVvsjKCi2tCK1YCIQKlpRgUaqpIQjgnJkE4JoY0ADKGpIswHZyDb398eYkZiF3BByZ+DzOuceZ+7mR54j+ea593ke/vOf/zB16lRCQ0OZNm0aaWlpZqJ5lz0r3B9XlQ+nklb0j3bQI8JDHymJiIhcQqaKka1btwJwzTV1L95mt9u5+uqrAdi3b1+j7/vss8+6P+fn57Nq1SpiY2OZOXMmZWVlZiJ6vuI8OPQGAKW21qyvGgLA1IHqFRERkSuTqXdGjh07BkDnzp3rPcfhcJCZmUlGRkaj7llSUsKjjz5KXl4eGRkZbN26la+//hqn00lCQgLp6emkpqYSEBBgJqrn+mQNVJYC8EbFYM4SSKSjDUN+1tHiYCIiItYwVYzk5OQAEBwcXO85/v6uFzALCgoadc82bdowZcoU93fDMFi1ahVPPPEEeXl5fPTRR8yZM4eXX3653ntkZWWRlZXl/l5UVATAoUOHCAoKcu8PDw8nPNzCOTwqy2BfIgBOfFhTdScADwyIwu6jSc5EROTKZOoxTX5+PuAqIOrjdLoWfDt37lyTAtlsNqZNm8aBAwfcj3xeeeUVvvrqq3qvSUhIoE+fPu5t8ODBAAwePLjG/oSEhCZlajafvQ3F/wUg1Xkz3xodCfZvxYTYuh97iYiIXAlM9Yz4+flRWVnZ4DkVFRUAtGvXrumpgMjISFJSUoiNjaW8vJzNmzczd+7cOs+dMWMGo0ePdn8vKipi8ODB7Ny5s1bPiGUMo+Zw3grXJGcTb76GIH/TI6xFREQuG6Z+CjocDkpKSigtLa33nDNnzgAQGhp6ccmAnj17MnXqVF5++eUGp5r/6eOXs2fPAtC7d2/PmfQsczv8Nx2Aw7bupBnXYbPB/f2jrM0lIiJiMVOPaWJiYgDIzs6u95zq+UXqGvrbFOPGjQPA19e3We5nmfN6RVaUuRbEGxYTRhdH/Y+8RERErgSmipGBAwcCrsnP6lJcXExubi4Aw4YNu8hoLpGRkYBrDhOvlZ0Ome8DcMonjFRnLKDhvCIiImCyGKleOff777/n1KlTtY5/9tlngOvdkiFDhjRDPNeoHF9fX8aMGdMs97PEnh97RV4uG4ETH3qEh9C3a3sLQ4mIiHgGU8VIjx49GDnS9YghJSWl1vH333f99j9lypQGh/+asWXLFh555BFrXz69GIXZcPgtAEp8Anm7yjXSZ+rArthsGs4rIiJiqhgBWLZsGa1btyYxMbHG/pKSElavXo3D4WDRokW1rouLiyMkJIQVK1bU2P/NN9+wadMmSkpKal2TkZHB3r17a8zQ6nX2r4KqcgDWlf+CYloTGuTHqBu9tLgSERFpZqaLkW7dupGUlMQnn3zC448/TllZGVlZWUycOJHTp0+zefNmwsLCalyTm5vL66+/TmFhIStXrqxx7KGHHmL8+PH06tWL5ORkysvLKS8vZ/369axYsYKNGzd67+yr5SWuYgSostlJqnQtJBjXLxL/VnYrk4mIiHgM08UIwKRJk9ixYwcHDhwgIiKC/v37ExUVRXp6OgMGDKh1fmhoKHFxcQQFBTFr1qwax/785z8zfPhw8vPzufvuu+nZsyezZs0iIiKC5cuXExgY2LT/Mk9weD2UuiaK+7fRn1M48LP7cG/fSIuDiYiIeA6bYRiG1SGa29mzZ2nbti1nzpyxbp4RpxP+dgvkfQnAr8oW8bkRzf/06czSCTdak0lEROQSaurP3yb1jEgjfJnqLkQO2W/gcyMagAdvjbIwlIiIiOdRMXKp7H7J/fGvpa53RfpFt+f6iLZWJRIREfFIKkYuhaxP4cRHro+tOrPdeRMAU2/VJGciIiI/pWLkUjhv6veXSodj4EOX9m0YGhPWwEUiIiJXJhUjze3Md/D5RgCK7W3ZWDUIgPsHRGH30SRnIiIiP6VipLntSwRnJQB/rxjCOfwJ8m/FxNjOFgcTERHxTCpGmlNZERxIAqDK5ktSuWuxwAmxnQkO8PJVh0VERC4RFSPN6dDrcO4MAFtsg8jhKmw2eGBAlLW5REREPJiKkebirII9P66781LpcADuiAkj0uHFs8iKiIhcYipGmkvGe1BwAoA035s4YnQBNJxXRETkQlSMNJfzJjl7odjVKxITHkK/6PZWJRIREfEKKkaaw8n9cHIvAN/7RfGhsxcAU2+NwmbTcF4REZGGqBhpDnt+nOTsxZLhgI3QID9G3RhhXSYREREvoWLkYhV8DenvAFDcqh3JlQMAuLdvJAG+diuTiYiIeAUVIxdrbwIYTgDWVg6nDD/87D7c26+LxcFERES8g4qRi3HuDBx8FYBKH39Wn/sFAKNujKBjcICVyURERLyGipGLcfBVKC8EIMV+O/mEAPDgrVEWhhIREfEuKkaaqqoC9rzs/rqs6A4AbunanhuubmtVKhEREa+jYqSp0t+Bs98CkBbQj0zjakCTnImIiJilYqQpDKPGJGd/PDsUgM7tWjOsR5hVqURERLySipGm+GY3fJ8GwHetu7Hb2QNwLYhn99EkZyIiImaoGGmKj3/sFVlWPAywEehnZ+LN11iXSURExEupGDErLxOObAGgyK8DyeV9AZgQew0hAb5WJhMREfFKKkbM2rMSMAD4e9UIKmiFzeZ6RCMiIiLmqRgxoyQfDr0OQKW9DQnFtwEw9GdhRIUGWplMRETEa6kYMeNAElSUAPBvvzs4SxAAUwdGWRhKRETEu6kYaazKctibCICBjT+edk39/rNOwfSPdliZTERExKupGGmszzdC0SkAPg0axDeGaz6RqQO7YrNpOK+IiEhTqRhprMpSaN0OgGcLhgDgCPRj9I0RVqYSERHxeipGGit2KvxvOu9cu4h9VdcBcG+/SAJ87RYHExER8W4qRkwoxZ//y+wO2PC124jr18XqSCIiIl5PxYgJm9K+5UxpBQCjekXQMTjA4kQiIiLeT8VIIzmdBmt2HXd/f1Cr84qIiDQLFSON9NFXuWTmFANwS1R7enZua3EiERGRy4OKkUbafzzf/VmTnImIiDSfVlYH8BaPjejOqBsj2HTwW4b16GR1HBERkcuGihETuncK5olfxlgdQ0RE5LKixzQiIiJiKRUjIiIiYqkmFyNfffUVkydPJjw8nA4dOjBp0iQyMzObMxszZszAZrOxY8eOZr2viIiIeI4mFSOpqancdNNNtGnThiNHjnDixAlCQkL4+c9/zt69e5sl2HvvvUdiYmKz3EtEREQ8l+liJDMzk/Hjx3PttdeSmJhISEgIgYGBrFy5kk6dOjFq1Cjy8vIuKlRubi7Tpk27qHuIiIiIdzBdjMydO5eioiJmz56Nj8+Pl7dq1YqZM2eSk5PD448/flGhpk+fzqRJky7qHiIiIuIdTBUjx48fZ/PmzQDccccdtY6PGDECgNdee63JvSNJSUkcPXqU5557rknXi4iIiHcxVYxs2bIFgMDAQLp2rb02S/fu3QkICKCsrIzk5GTTYU6cOMH8+fN57bXXCAjQInQiIiJXAlPFyNatWwG45ppr6jxut9u5+uqrAdi3b5+pIE6nk/vvv5/58+fTu3dvU9e2lKysLOLj48nKyrI6ijSS2sz7qM28j9rM+3ham5kqRo4dOwZA586d6z3H4XAAkJGRYSrI0qVLsdlszJs3z9R14PpDPXjwoHs7dOgQAIcOHaqx/2L/0LOysli4cKHHNJ5cmNrM+6jNvI/azPt4WpuZmg4+JycHgODg4HrP8ff3B6CgoKDR9z18+DBLly7lk08+qfFSbGMlJCSwcOHCWvsHDx5c4/uCBQuIj483fX8RERG5dEwVI/n5rpVr27RpU+85TqcTgHPnzjXqnmVlZcTFxfGXv/yFLl26mInjNmPGDEaPHu3+XlRUxODBg9m5cydBQUHu/eHh4U26v4iIiFw6pooRPz8/KisrGzynoqICgHbt2jXqnk8++SQxMTHExcWZiVJDeHh4jULjzJkzAERHRxMSElLj3LNnzzb531NUVOT+58XcR1qO2sz7qM28j9rM+1yqNqu+l2EYpq4zVYw4HA5KSkooLS2t95zqQiA0NPSC99u5cydvvfUWn376qZkYF1RYWAjU/6Ltxfrp4x/xfGoz76M28z5qM+9zqdqssLCQtm3bNvp8U8VITEwMJ0+eJDs7u95zqucXqWvo708tWrSIrKysCxYNI0eOxG63M2jQIFJSUi5434iICE6ePElwcDA2m+2C5zfWoUOH3I9/PHXEj9SkNvM+ajPvozbzPpeqzQzDoLCwkIiICFPXmSpGBg4cSGpqKsePH6/zeHFxMbm5uQAMGzbsgverqKigsrLygo9+qt8/aahH5nw+Pj4Njvhpqur3T4KCgmo9/hHPpDbzPmoz76M28z6Xss3M9IhUMzV0Zfz48QB8//33nDp1qtbxzz77DHC9WzJkyJAL3m/Hjh0YhlHvVu2DDz7AMAyt3isiInIZMlWM9OjRg5EjRwLU+bjk/fffB2DKlCkNDv8VERERqWZ6Uo9ly5bRunVrEhMTa+wvKSlh9erVOBwOFi1aVOu6uLg4QkJCWLFiRdPTWiw8PJwFCxZoiLAXUZt5H7WZ91GbeR9PazObYXb8DfDmm28SFxfH3Llzefrpp8nPz2fatGl8/PHHvPvuuwwYMKDG+bm5uXTo0AGAG264wf0454Lhfnj59IMPPuD22283G1NERES8gPnpToFJkyaxY8cODhw4QEREBP379ycqKor09PRahQi4hvnGxcURFBTErFmzLjq0iIiIXD6a1DMiIiIi0lya1DMiIiIi0lxUjIiIiIilVIw0wueff87dd99NWFgY/v7+REdH8+ijj3rM0stS0/Hjx4mLi6Njx474+/vTrVs34uPjG714o3iO5557DpvNxtq1a62OIg3o27cvNput1tanTx+ro8kFlJeX849//IN77rmHIUOG8MADD7B169aWD2JIg1JSUoyAgAADqLV16NDBOHjwoNUR5TwZGRmGw+EwfH19DYfDUaO9Ro4caXU8MSEtLc3w8/MzACMpKcnqOFKP1NTUOv9+BIxNmzZZHU8akJqaanTt2tWIjY01PvzwQ0uzqGekAXl5eUyePJmbb76Z5ORkjhw5wo4dO7jrrrsAyMnJ4de//jUlJSUWJxVwLS8wYcIEFi5cSGFhIbm5uRw9etQ9wislJYV//vOfFqeUxjh37hxxcXHuVcDFcy1evJi1a9fyxRdf1NrGjBljdTypx5IlSxgxYgRDhw7l448/ZtCgQdYGsrQU8nBLliwx7rvvPsPpdNbY73Q6jcmTJ7ur/1WrVlmUUM73t7/9zXj33Xdr7c/OzjauuuoqAzB++9vfWpBMzJozZ44xZswYIzIyUj0jHmzXrl3G9ddfX+vvSPFszz//vAEY999/v9VR3NQz0oD9+/ezcuXKWiv/2mw2XnjhBVq1cq0zePDgQSviyU889NBD7l6r83Xs2JF+/foBcNVVV7V0LDFp+/btvPnmm7zyyitWR5ELWLx4MRMmTKixlph4trfffpvf/e53xMTEkJCQYHUcNxUjDViyZAmBgYF1HuvYsSPXX389AAEBAS0ZS+rh7+9f77GAgADsdjuTJ09uwURi1unTp3nwwQdZtWoVoaGhVseRBqSlpZGSkkJ8fDxt27blnnvuYefOnVbHkgYUFBQwa9YsDMNg2bJlDf6d2dJUjDTg2muvbfB4dUPecMMNLRFHmqiyspI9e/YQHx9Pt27drI4jDfjNb37Dr371K375y19aHUUuYPHixe7PRUVFrF+/nttvv53x48dz+vRpC5NJfRYuXEheXh69evVi+PDhVsepQcVIExmGQWZmJv7+/npJy8MtWrSICRMm8NRTT1kdRRrw9ttvc+DAAf70pz9ZHUUuwDAMxo4dy/Lly3nkkUeIiYlxH9u0aRMDBgwgNzfXwoTyU+fOnWPNmjUAjBo1iuXLlzN06FDCwsJo3749Q4cOZfv27dYFtPaVFe+1e/duAzDmzJljdRSpR1ZWlvHwww8bgNGzZ09j27ZtVkeSenz33XdGWFiYsX///hr79QKr99i4caPRpUsX94v9I0aMsDqSnOedd95xt02vXr2MDRs2GGfOnDEKCgqMBQsWGDabzbDZbJYNyFAx0kSTJk0ywsPDjYKCAqujSB1eeOEFo2fPnjXmPGjVqpWxYcMGq6NJHUaMGGE888wztfarGPEu+fn5Nf6/2759u9WR5AePPfaYARj+/v5GeXl5reOzZ882AMPPz8/IyMho8Xx6TNMEu3fvZsOGDaxdu1ajMzzUnDlzOHz4MEePHnW/tFpZWcnMmTMpLi62OJ2c76WXXqKwsJAnnnjC6ihykdq1a8fWrVtp3749gOb18SDfffcdAOHh4fj6+tY6Pm/ePGw2G+Xl5axcubKl4+mdEbOKi4t5+OGHefrppz3uBSCp7brrruP1119nwYIFAOTm5pKSkmJxKql25MgRnnnmGdatW4fdbrc6jjSD8PBw5s+fD0BmZqbFaaRa9Ts87dq1q/N4dHQ0t9xyC+D6hbulqRgxafr06fTv35/f//73VkcRE5588kk6deoEwLFjxyxOI9WWLl3qfrs/KCio1vbNN98AMHPmTIKCgtzD6cWzjRs3DqDO38DFGtW9+E6ns95zqkcbFhQUtEim86kYMSE+Pp7S0lKPmihGGsfX15eRI0cC0Lp1a4vTSLWKigqqqqooLi6uczN+mEyrrKzMvU88X2RkJADdu3e3OIlUi46OBn58XFOXiIgIABwOR4tkOp+KkUZKTEzko48+4o033lB3speq7hmp7ooU661duxbD9SJ9nVv1D7WkpCQMw+DEiRPWBpZGqf7NeuLEiRYnkWp33nkn4Hpc8+WXX9Z5TnW7Vc9Y3ZJUjDTC2rVrWbNmDcnJyXXOWFdZWcmGDRssSCZmHD16lN69e9O3b1+ro4hc1rZs2cLYsWPp06eP1VHkB7fddhu9evUCYP369XWek56ejs1mY+rUqS0ZDVAxckGvvvoqzz//POvWraOsrIzc3Fxyc3PJycnh5MmTbNu2jTvvvJOysjKro17xzp49W+8Lc59++inbtm1j9erVLZxK5PKTk5PDxo0byc/Pr3UsOzubdevWaW0hD+Pj48OLL76I3W7nr3/9K6dOnapxfPfu3ezatYuZM2fSs2fPlg/Y4oOJvciKFSsMm81WY66KuragoCCjqKjI6rhXvOr5Dfr27Wu89957Rnl5uVFVVWVs3rzZuOmmm4zdu3dbHVFM0jwjnql6MsGwsDAjKSnJKCkpMSorK40tW7YY06ZNM7Kzs62OKPVISEgwbDabERsba3zxxReGYbhWX+7SpYtx11131TkHSUtQz0g9kpOTmT17dqNWoxw7dmy9C+pJy3nsscfo3r07aWlpjBkzhq5duzJ69GiOHz/Orl27LHkOKnI5+sMf/sC4ceOoqqpixowZ9OjRg/vuu4+qqioSExPp2LGj1RGlHtOnTyc1NZXg4GBiY2Pp0KED8+bN46mnnuJf//qXZSOgbEZjftqKiIiIXCLqGRERERFLqRgRERERS6kYEREREUupGBERERFLqRgRERERS6kYEREREUupGBERERFLqRgRERERS6kYEREREUupGBERERFLqRgRERERS6kYEREREUupGBERERFLqRgRERERS/0/rEWqK6GuAa8AAAAASUVORK5CYII=",
      "text/plain": [
       "Figure(PyObject <Figure size 600x400 with 1 Axes>)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "rc(\"figure\", figsize = (6,4))\n",
    "plot( Ds .+ 1, Es, label = \"DS\");\n",
    "plot( (Ds .+ 1)[1:end-1], ERA, label = \"RA\");\n",
    "legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "At least our methods do yield the same results (so both are correct).\n",
    "\n",
    "In DS we have a _dedicated and different_ function for the E2 quantity of Cao. Most of the time you do not need both quantities at the same time, so there is no reason to compute them both..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Conclusions\n",
    "\n",
    "In this notebook I have only reviewed common features and the most central function `recurrencematrix`. I have also read the README and have an overview of all features. I think RA is a great package and it has a lot of functionality and also very useful functionality!\n",
    "\n",
    "Here is what I believe are the best steps forward:\n",
    "\n",
    "0. Move `RecurrenceAnalysis` to JuliaDynamics org. (When you do I will give you owner priviliges for it)\n",
    "1. Put `DynamicalSystemsBase` to the dependencies of RA.\n",
    "1. Use `reconstruct` from DynamicalSystemsBase exclusively. It is too many more features and it is much faster! This also means that all `embed` code is discarded.\n",
    "2. As far as `distancematrix` is concerned, we could mainly use the method I wrote in this notebook, which uses a `Dataset`. I suggest to write a higher level if-clause that if the dimensionality of the dataset is more than let's say 10, it converts it to a matrix and uses your method. This means that both mine and your methods are written and we choose which one to use!\n",
    "3. I think `recurrencematrix` could be made more intuitive. But on the other hand your method is just so darn fast!!! I don't know why, but I suggest that we keep your method for now and see how it goes.\n",
    "4. Move all delay embedding dimension estimation functionality `estimate_dimension` to DynamicalSystemsBase. At the moment it lives in `ChaosTools`. (I'll do that)\n",
    "5. Use `estimate_dimension` from DynSysBase exclusively. Contribute your extra two methods FNN and FFNN to DynSysBase. I will be overviewing the contribution PRs and be sure we also get better performance (they should also be reworked to \"just return the values at the given dimensions\" instead of only first and last).\n",
    "6. Compare the mutual information methods between the one in ChaosTools (in `estimate_delay`) and the one in RA. Use the best one exclusively in both ChaosTools and RA. It is already clear that the method in RA **has more features** than then one in ChaosTools, the only thing left to check is the performance.\n",
    "6. [Optional/Suggestion] : Rework the handling of the metrics: instead of asking for a string ask for the `Metric` instance. E.g. `\"euclidean\" => Euclidean()`. This is important for education as it teaches the Julia way. But this is just my humble opinion.\n",
    "8. Add citations to documentation strings.\n",
    "\n",
    "Once these are done, and if you agree of course, we make RA an official part of DS by \n",
    "\n",
    "1. adding it to its documentation page as a set of dedicated pages (much like `ChaosTools`). For this to happen we will also have to write a proper documentation (i.e. separate / expand the current README to documenter acceptable files that expand the docstrings. I will take care of that).\n",
    "2. re-export it. \n",
    "3. Make sure the source code is clear. I will read the source and wherever there is something unclear I'll either clarify myself or ask for your help in claryfying.\n",
    "\n",
    "You might think number 3 is irrelevant or useless or unecessary but the philoshopy of JuliaDynamics is that the source code is clear and understandable. "
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 1.0.0",
   "language": "julia",
   "name": "julia-1.0"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.0.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}