2020/Week16/OrderReductionMethods.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Method 1 (Combastel 2003)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "using Revise, ReachabilityAnalysis, Test, BenchmarkTools, Plots\n",
    "using StaticArrays"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/svg+xml": [
       "<?xml version=\"1.0\" encoding=\"utf-8\"?>\n",
       "<svg xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\" width=\"600\" height=\"400\" viewBox=\"0 0 2400 1600\">\n",
       "<defs>\n",
       "  <clipPath id=\"clip3600\">\n",
       "    <rect x=\"0\" y=\"0\" width=\"2400\" height=\"1600\"/>\n",
       "  </clipPath>\n",
       "</defs>\n",
       "<path clip-path=\"url(#clip3600)\" d=\"\n",
       "M0 1600 L2400 1600 L2400 0 L0 0  Z\n",
       "  \" fill=\"#ffffff\" fill-rule=\"evenodd\" fill-opacity=\"1\"/>\n",
       "<defs>\n",
       "  <clipPath id=\"clip3601\">\n",
       "    <rect x=\"480\" y=\"0\" width=\"1681\" height=\"1600\"/>\n",
       "  </clipPath>\n",
       "</defs>\n",
       "<path clip-path=\"url(#clip3600)\" d=\"\n",
       "M141.865 1487.47 L2352.76 1487.47 L2352.76 47.2441 L141.865 47.2441  Z\n",
       "  \" fill=\"#ffffff\" fill-rule=\"evenodd\" fill-opacity=\"1\"/>\n",
       "<defs>\n",
       "  <clipPath id=\"clip3602\">\n",
       "    <rect x=\"141\" y=\"47\" width=\"2212\" height=\"1441\"/>\n",
       "  </clipPath>\n",
       "</defs>\n",
       "<polyline clip-path=\"url(#clip3602)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  234.389,1487.47 234.389,47.2441 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip3602)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  552.458,1487.47 552.458,47.2441 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip3602)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  870.527,1487.47 870.527,47.2441 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip3602)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  1188.6,1487.47 1188.6,47.2441 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip3602)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  1506.67,1487.47 1506.67,47.2441 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip3602)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  1824.73,1487.47 1824.73,47.2441 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip3602)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  2142.8,1487.47 2142.8,47.2441 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip3602)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  141.865,1306.82 2352.76,1306.82 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip3602)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  141.865,1101.31 2352.76,1101.31 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip3602)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  141.865,895.79 2352.76,895.79 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip3602)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  141.865,690.274 2352.76,690.274 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip3602)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  141.865,484.758 2352.76,484.758 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip3602)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  141.865,279.242 2352.76,279.242 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip3602)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  141.865,73.7256 2352.76,73.7256 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip3600)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  141.865,1487.47 2352.76,1487.47 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip3600)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  141.865,1487.47 141.865,47.2441 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip3600)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  234.389,1487.47 234.389,1470.19 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip3600)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  552.458,1487.47 552.458,1470.19 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip3600)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  870.527,1487.47 870.527,1470.19 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip3600)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  1188.6,1487.47 1188.6,1470.19 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip3600)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  1506.67,1487.47 1506.67,1470.19 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip3600)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  1824.73,1487.47 1824.73,1470.19 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip3600)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  2142.8,1487.47 2142.8,1470.19 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip3600)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  141.865,1306.82 168.395,1306.82 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip3600)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  141.865,1101.31 168.395,1101.31 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip3600)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  141.865,895.79 168.395,895.79 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip3600)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  141.865,690.274 168.395,690.274 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip3600)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  141.865,484.758 168.395,484.758 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip3600)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  141.865,279.242 168.395,279.242 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip3600)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  141.865,73.7256 168.395,73.7256 \n",
       "  \"/>\n",
       "<g clip-path=\"url(#clip3600)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 234.389, 1541.47)\" x=\"234.389\" y=\"1541.47\">-3</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip3600)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 552.458, 1541.47)\" x=\"552.458\" y=\"1541.47\">-2</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip3600)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 870.527, 1541.47)\" x=\"870.527\" y=\"1541.47\">-1</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip3600)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 1188.6, 1541.47)\" x=\"1188.6\" y=\"1541.47\">0</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip3600)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 1506.67, 1541.47)\" x=\"1506.67\" y=\"1541.47\">1</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip3600)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 1824.73, 1541.47)\" x=\"1824.73\" y=\"1541.47\">2</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip3600)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 2142.8, 1541.47)\" x=\"2142.8\" y=\"1541.47\">3</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip3600)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:end;\" transform=\"rotate(0, 117.865, 1324.32)\" x=\"117.865\" y=\"1324.32\">-2</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip3600)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:end;\" transform=\"rotate(0, 117.865, 1118.81)\" x=\"117.865\" y=\"1118.81\">-1</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip3600)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:end;\" transform=\"rotate(0, 117.865, 913.29)\" x=\"117.865\" y=\"913.29\">0</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip3600)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:end;\" transform=\"rotate(0, 117.865, 707.774)\" x=\"117.865\" y=\"707.774\">1</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip3600)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:end;\" transform=\"rotate(0, 117.865, 502.258)\" x=\"117.865\" y=\"502.258\">2</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip3600)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:end;\" transform=\"rotate(0, 117.865, 296.742)\" x=\"117.865\" y=\"296.742\">3</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip3600)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:end;\" transform=\"rotate(0, 117.865, 91.2256)\" x=\"117.865\" y=\"91.2256\">4</text>\n",
       "</g>\n",
       "<path clip-path=\"url(#clip3602)\" d=\"\n",
       "M204.437 1446.71 L678.833 1410 L1270.03 1164.01 L1619.62 1003.95 L1952.89 737.483 L2192.37 492.242 L2290.18 88.0053 L1815.79 124.721 L1224.59 370.71 L875 530.769 \n",
       "  L541.732 797.236 L302.248 1042.48 L204.437 1446.71 L204.437 1446.71  Z\n",
       "  \" fill=\"#009af9\" fill-rule=\"evenodd\" fill-opacity=\"0.5\"/>\n",
       "<polyline clip-path=\"url(#clip3602)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:0.5; fill:none\" points=\"\n",
       "  204.437,1446.71 678.833,1410 1270.03,1164.01 1619.62,1003.95 1952.89,737.483 2192.37,492.242 2290.18,88.0053 1815.79,124.721 1224.59,370.71 875,530.769 \n",
       "  541.732,797.236 302.248,1042.48 204.437,1446.71 \n",
       "  \"/>\n",
       "</svg>\n"
      ]
     },
     "execution_count": 47,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Z = Zonotope(rand(2), rand(2, 6)); # order 3\n",
    "plot(Z)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2×10 Array{Float64,2}:\n",
       " 0.569054  0.122998  0.295936  0.0510516  …  0.82904   0.65224  0.576509\n",
       " 0.187848  0.36725   0.783097  0.174121      0.238775  0.58694  0.404778"
      ]
     },
     "execution_count": 44,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# compute the interval hull\n",
    "M = Z.generators"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2×1 Array{Float64,2}:\n",
       " 2.2718354808591217\n",
       " 2.6949894545555884"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sum(abs.(M), dims=2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "interval_hull_v0 (generic function with 1 method)"
      ]
     },
     "execution_count": 37,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# interval hull\n",
    "function interval_hull_v0(G::Matrix{N}) where {N}\n",
    "    n = size(G, 1)\n",
    "    Gsum = sum(abs.(G), dims=2)\n",
    "    Gbox = zeros(N, n, n)\n",
    "    @inbounds for i in 1:n\n",
    "        Gbox[i, i] = Gsum[i]\n",
    "    end\n",
    "    return Gbox\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  535.608 ns (9 allocations: 576 bytes)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "2×2 Array{Float64,2}:\n",
       " 3.08868  0.0\n",
       " 0.0      3.00576"
      ]
     },
     "execution_count": 42,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "@btime interval_hull_v0($M) # order 3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  552.107 ns (9 allocations: 640 bytes)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "2×2 Array{Float64,2}:\n",
       " 4.31308  0.0\n",
       " 0.0      3.86407"
      ]
     },
     "execution_count": 45,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "@btime interval_hull_v0($M) # order 10"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "interval_hull (generic function with 1 method)"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# interval hull\n",
    "function interval_hull(G::Matrix{N}) where {N}\n",
    "    n, p = size(G)\n",
    "    Gbox = zeros(N, n, n)\n",
    "    @inbounds for i in 1:n\n",
    "        for j in 1:p\n",
    "            Gbox[i, i] += abs(G[i, j])\n",
    "        end\n",
    "    end\n",
    "    return Gbox\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  38.217 ns (1 allocation: 112 bytes)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "2×2 Array{Float64,2}:\n",
       " 6.57925  0.0\n",
       " 0.0      4.75817"
      ]
     },
     "execution_count": 35,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "@btime interval_hull($M) # order 10"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  38.298 ns (1 allocation: 112 bytes)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "2×2 Array{Float64,2}:\n",
       " 4.31308  0.0\n",
       " 0.0      3.86407"
      ]
     },
     "execution_count": 46,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "@btime interval_hull($M) # order 3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "interval_hull (generic function with 2 methods)"
      ]
     },
     "execution_count": 50,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function interval_hull(Z::Zonotope)\n",
    "    Zonotope(Z.center, interval_hull(Z.generators))\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 82,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Zonotope{Float64,Array{Float64,1},Array{Float64,2}}([0.5680414549351669, 0.6271652184216445], [0.28161605196515205 -0.3209667418511202 … -0.2670767663612974 0.2444286112146754; -0.4478900355327835 -0.37906366892683785 … -0.2530432395046349 0.14618542154448555])"
      ]
     },
     "execution_count": 82,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Z"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 89,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/svg+xml": [
       "<?xml version=\"1.0\" encoding=\"utf-8\"?>\n",
       "<svg xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\" width=\"600\" height=\"400\" viewBox=\"0 0 2400 1600\">\n",
       "<defs>\n",
       "  <clipPath id=\"clip4000\">\n",
       "    <rect x=\"0\" y=\"0\" width=\"2400\" height=\"1600\"/>\n",
       "  </clipPath>\n",
       "</defs>\n",
       "<path clip-path=\"url(#clip4000)\" d=\"\n",
       "M0 1600 L2400 1600 L2400 0 L0 0  Z\n",
       "  \" fill=\"#ffffff\" fill-rule=\"evenodd\" fill-opacity=\"1\"/>\n",
       "<defs>\n",
       "  <clipPath id=\"clip4001\">\n",
       "    <rect x=\"480\" y=\"0\" width=\"1681\" height=\"1600\"/>\n",
       "  </clipPath>\n",
       "</defs>\n",
       "<path clip-path=\"url(#clip4000)\" d=\"\n",
       "M141.865 1487.47 L2352.76 1487.47 L2352.76 47.2441 L141.865 47.2441  Z\n",
       "  \" fill=\"#ffffff\" fill-rule=\"evenodd\" fill-opacity=\"1\"/>\n",
       "<defs>\n",
       "  <clipPath id=\"clip4002\">\n",
       "    <rect x=\"141\" y=\"47\" width=\"2212\" height=\"1441\"/>\n",
       "  </clipPath>\n",
       "</defs>\n",
       "<polyline clip-path=\"url(#clip4002)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  224.59,1487.47 224.59,47.2441 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip4002)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  876.818,1487.47 876.818,47.2441 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip4002)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  1529.05,1487.47 1529.05,47.2441 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip4002)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  2181.27,1487.47 2181.27,47.2441 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip4002)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  141.865,1302.6 2352.76,1302.6 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip4002)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  141.865,973.658 2352.76,973.658 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip4002)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  141.865,644.72 2352.76,644.72 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip4002)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  141.865,315.781 2352.76,315.781 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip4000)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  141.865,1487.47 2352.76,1487.47 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip4000)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  141.865,1487.47 141.865,47.2441 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip4000)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  224.59,1487.47 224.59,1470.19 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip4000)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  876.818,1487.47 876.818,1470.19 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip4000)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  1529.05,1487.47 1529.05,1470.19 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip4000)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  2181.27,1487.47 2181.27,1470.19 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip4000)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  141.865,1302.6 168.395,1302.6 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip4000)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  141.865,973.658 168.395,973.658 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip4000)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  141.865,644.72 168.395,644.72 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip4000)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  141.865,315.781 168.395,315.781 \n",
       "  \"/>\n",
       "<g clip-path=\"url(#clip4000)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 224.59, 1541.47)\" x=\"224.59\" y=\"1541.47\">-1</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip4000)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 876.818, 1541.47)\" x=\"876.818\" y=\"1541.47\">0</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip4000)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 1529.05, 1541.47)\" x=\"1529.05\" y=\"1541.47\">1</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip4000)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 2181.27, 1541.47)\" x=\"2181.27\" y=\"1541.47\">2</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip4000)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:end;\" transform=\"rotate(0, 117.865, 1320.1)\" x=\"117.865\" y=\"1320.1\">-1</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip4000)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:end;\" transform=\"rotate(0, 117.865, 991.158)\" x=\"117.865\" y=\"991.158\">0</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip4000)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:end;\" transform=\"rotate(0, 117.865, 662.22)\" x=\"117.865\" y=\"662.22\">1</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip4000)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:end;\" transform=\"rotate(0, 117.865, 333.281)\" x=\"117.865\" y=\"333.281\">2</text>\n",
       "</g>\n",
       "<path clip-path=\"url(#clip4002)\" d=\"\n",
       "M204.437 888.425 L571.793 1183.08 L1158.09 1446.71 L1476.93 1350.54 L1825.32 1184.07 L2244.01 934.693 L2290.18 646.294 L1922.83 351.638 L1336.53 88.0053 L1017.69 184.177 \n",
       "  L669.296 350.649 L250.609 600.026 L204.437 888.425 L204.437 888.425  Z\n",
       "  \" fill=\"#009af9\" fill-rule=\"evenodd\" fill-opacity=\"0.5\"/>\n",
       "<polyline clip-path=\"url(#clip4002)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:0.5; fill:none\" points=\"\n",
       "  204.437,888.425 571.793,1183.08 1158.09,1446.71 1476.93,1350.54 1825.32,1184.07 2244.01,934.693 2290.18,646.294 1922.83,351.638 1336.53,88.0053 1017.69,184.177 \n",
       "  669.296,350.649 250.609,600.026 204.437,888.425 \n",
       "  \"/>\n",
       "<path clip-path=\"url(#clip4002)\" d=\"\n",
       "M2290.18 88.0053 L204.437 88.0053 L204.437 1446.71 L2290.18 1446.71 L2290.18 88.0053 L2290.18 88.0053  Z\n",
       "  \" fill=\"#e26f46\" fill-rule=\"evenodd\" fill-opacity=\"0.5\"/>\n",
       "<polyline clip-path=\"url(#clip4002)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:0.5; fill:none\" points=\"\n",
       "  2290.18,88.0053 204.437,88.0053 204.437,1446.71 2290.18,1446.71 2290.18,88.0053 \n",
       "  \"/>\n",
       "</svg>\n"
      ]
     },
     "execution_count": 89,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#Z = rand(Zonotope)\n",
    "plot(Z)\n",
    "plot!(interval_hull(Z))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Chaning the order"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 308,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "┌ Info: Precompiling LoopVectorization [bdcacae8-1622-11e9-2a5c-532679323890]\n",
      "└ @ Base loading.jl:1260\n"
     ]
    }
   ],
   "source": [
    "using LoopVectorization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 319,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "interval_hull_avx (generic function with 1 method)"
      ]
     },
     "execution_count": 319,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# interval hull\n",
    "function interval_hull_avx(G::Matrix{N}) where {N}\n",
    "    n, p = size(G)\n",
    "    Gbox = zeros(N, n, n)\n",
    "    @avx for j in 1:p\n",
    "        for i in 1:n\n",
    "            Gbox[i, i] += abs(G[i, j])\n",
    "        end\n",
    "    end\n",
    "    return Gbox\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 314,
   "metadata": {},
   "outputs": [],
   "source": [
    "G = rand(3, 20) .- 0.5;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 320,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  47.666 ns (1 allocation: 160 bytes)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "3×3 Array{Float64,2}:\n",
       " 4.66442  0.0      0.0\n",
       " 0.0      4.55711  0.0\n",
       " 0.0      0.0      3.54845"
      ]
     },
     "execution_count": 320,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "@btime interval_hull_avx($G)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 318,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  62.164 ns (1 allocation: 160 bytes)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "3×3 Array{Float64,2}:\n",
       " 4.66442  0.0      0.0\n",
       " 0.0      4.55711  0.0\n",
       " 0.0      0.0      3.54845"
      ]
     },
     "execution_count": 318,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "@btime interval_hull($G)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Sorting**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 80,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2×6 Array{Float64,2}:\n",
       "  0.281616  -0.320967  -0.449456  -0.0353958  -0.267077  0.244429\n",
       " -0.44789   -0.379064   0.400732  -0.438377   -0.253043  0.146185"
      ]
     },
     "execution_count": 80,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Z = Zonotope(rand(2), rand(2, 6) .- 0.5)\n",
    "M = Z.generators"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 138,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "sort_gens (generic function with 1 method)"
      ]
     },
     "execution_count": 138,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function sort_gens(G::Matrix{N}) where {N}\n",
    "    p = size(G, 2)\n",
    "    weights = zeros(N, p)\n",
    "    @inbounds for i in 1:p\n",
    "        weights[i] = norm(view(G, :, i), 2)\n",
    "    end\n",
    "    indices = sortperm(weights, rev=false) # ordena de menor a mayor\n",
    "\n",
    "    weights, indices\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 139,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "weights = [0.5290681285562354, 0.496698011345726, 0.6021596560622862, 0.43980407256177734, 0.3679142293918612, 0.28480787112099604]\n",
      "indices = [6, 5, 4, 2, 1, 3]\n",
      "  255.749 ns (12 allocations: 624 bytes)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "([0.5290681285562354, 0.496698011345726, 0.6021596560622862, 0.43980407256177734, 0.3679142293918612, 0.28480787112099604], [6, 5, 4, 2, 1, 3])"
      ]
     },
     "execution_count": 139,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "weights, indices = sort_gens(M)\n",
    "@show weights\n",
    "@show indices;\n",
    "@btime sort_gens($M)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 204,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "sort_gens2 (generic function with 1 method)"
      ]
     },
     "execution_count": 204,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function sort_gens2(G::Matrix{N}) where {N}\n",
    "    p = size(G, 2)\n",
    "    weights = zeros(N, p)\n",
    "    @inbounds for i in 1:p\n",
    "        weights[i] = norm(view(G, :, i), 2)\n",
    "    end\n",
    "    indices = Vector{Int}(undef, p)\n",
    "    sortperm!(indices, weights, rev=true, initialized=false)\n",
    "\n",
    "    weights, indices\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 137,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "weights = [0.5290681285562354, 0.496698011345726, 0.6021596560622862, 0.43980407256177734, 0.3679142293918612, 0.28480787112099604]\n",
      "indices = [6, 5, 4, 2, 1, 3]\n",
      "  255.234 ns (12 allocations: 624 bytes)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "([0.5290681285562354, 0.496698011345726, 0.6021596560622862, 0.43980407256177734, 0.3679142293918612, 0.28480787112099604], [6, 5, 4, 2, 1, 3])"
      ]
     },
     "execution_count": 137,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "weights, indices = sort_gens2(M)\n",
    "@show weights\n",
    "@show indices;\n",
    "@btime sort_gens2($M)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 221,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "sort_gens_partial (generic function with 1 method)"
      ]
     },
     "execution_count": 221,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function sort_gens_partial(G::Matrix{N}, m::Int) where {N}\n",
    "    p = size(G, 2)\n",
    "    weights = zeros(N, p)\n",
    "    @inbounds for i in 1:p\n",
    "        weights[i] = norm(view(G, :, i), 2)\n",
    "    end\n",
    "    indices = [1:p;]\n",
    "    partialsortperm!(indices, weights, 1:m, rev=true, initialized=true)\n",
    "\n",
    "    weights, indices\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 352,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "ordenar_parcial2 (generic function with 1 method)"
      ]
     },
     "execution_count": 352,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function ordenar(weights)\n",
    "    p = length(weights)\n",
    "    indices = Vector{Int}(undef, p)\n",
    "    sortperm!(indices, weights, rev=false, initialized=false)\n",
    "    return indices\n",
    "end\n",
    "\n",
    "function ordenar_parcial(weights, m)\n",
    "    p = length(weights)\n",
    "    partialsortperm(weights, m, rev=false, alg=MergeSort)\n",
    "end\n",
    "\n",
    "function ordenar_parcial2(weights, m)\n",
    "    sortperm(weights; alg=PartialQuickSort(m))\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 353,
   "metadata": {},
   "outputs": [],
   "source": [
    "weights = rand(1000) .- 0.5;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 354,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ordenar(weights) = [414, 836, 422, 613, 838, 116, 473, 823, 911, 220, 709, 272, 167, 423, 596, 343, 32, 553, 498, 786, 995, 852, 165, 700, 839, 15, 741, 628, 107, 627, 545, 534, 584, 48, 742, 2, 252, 810, 97, 617, 146, 999, 233, 89, 219, 952, 389, 788, 996, 803, 143, 392, 757, 351, 530, 383, 14, 693, 703, 245, 63, 202, 344, 589, 106, 948, 691, 196, 960, 438, 234, 884, 330, 640, 848, 427, 629, 725, 13, 795, 470, 608, 487, 845, 739, 454, 745, 962, 391, 754, 226, 466, 949, 214, 173, 352, 972, 50, 254, 273, 285, 893, 804, 318, 633, 859, 336, 561, 400, 376, 357, 105, 432, 868, 112, 380, 632, 684, 369, 60, 894, 382, 609, 335, 819, 809, 736, 76, 523, 482, 396, 918, 247, 953, 653, 144, 337, 440, 184, 671, 88, 157, 210, 387, 978, 618, 892, 869, 286, 820, 259, 443, 408, 639, 201, 908, 966, 998, 563, 462, 682, 968, 390, 359, 923, 217, 576, 669, 95, 940, 783, 664, 56, 195, 573, 955, 890, 502, 244, 109, 275, 843, 572, 29, 826, 346, 153, 853, 651, 969, 986, 546, 197, 605, 215, 883, 5, 85, 769, 483, 491, 751, 930, 266, 24, 680, 72, 284, 624, 388, 832, 660, 965, 588, 881, 264, 270, 79, 985, 578, 652, 687, 780, 127, 910, 43, 516, 547, 140, 138, 793, 536, 189, 262, 937, 11, 474, 512, 433, 614, 730, 69, 540, 729, 324, 424, 650, 55, 301, 110, 903, 973, 456, 513, 620, 696, 775, 974, 265, 641, 527, 204, 537, 552, 874, 590, 977, 670, 909, 111, 807, 569, 661, 964, 406, 151, 122, 674, 665, 114, 615, 844, 950, 381, 84, 791, 420, 160, 878, 331, 762, 255, 698, 731, 360, 945, 341, 967, 772, 183, 885, 750, 119, 415, 942, 261, 113, 12, 279, 6, 493, 126, 726, 299, 571, 566, 480, 448, 499, 417, 159, 46, 476, 936, 251, 231, 248, 789, 761, 549, 525, 913, 306, 587, 211, 27, 118, 139, 963, 905, 51, 323, 827, 924, 654, 601, 83, 808, 734, 178, 798, 835, 681, 38, 176, 776, 900, 855, 431, 521, 288, 805, 446, 418, 866, 857, 901, 860, 634, 271, 984, 62, 441, 373, 260, 556, 689, 326, 870, 871, 567, 307, 797, 293, 858, 902, 102, 281, 386, 555, 850, 439, 212, 954, 147, 662, 484, 329, 604, 496, 685, 117, 311, 888, 91, 934, 416, 834, 941, 928, 459, 846, 23, 190, 830, 86, 428, 287, 577, 238, 101, 350, 164, 635, 637, 70, 994, 579, 485, 465, 481, 385, 339, 430, 322, 593, 507, 560, 185, 477, 216, 145, 586, 951, 611, 28, 997, 929, 535, 333, 946, 120, 108, 779, 699, 349, 802, 463, 206, 679, 849, 374, 645, 705, 161, 25, 289, 243, 291, 767, 296, 305, 583, 297, 692, 753, 701, 599, 714, 551, 9, 814, 784, 90, 519, 35, 544, 782, 933, 690, 277, 818, 327, 777, 898, 880, 274, 568, 437, 623, 1, 712, 847, 676, 622, 492, 18, 104, 425, 228, 636, 720, 149, 694, 156, 961, 897, 794, 815, 558, 461, 338, 873, 550, 959, 489, 221, 64, 971, 768, 192, 241, 724, 447, 925, 854, 230, 796, 678, 472, 518, 125, 659, 657, 865, 171, 99, 919, 154, 225, 59, 907, 328, 514, 53, 175, 16, 715, 94, 115, 648, 378, 309, 812, 186, 208, 98, 575, 479, 510, 524, 533, 401, 141, 862, 445, 500, 96, 723, 181, 831, 177, 879, 926, 708, 379, 856, 891, 876, 863, 743, 364, 792, 402, 697, 367, 686, 128, 531, 763, 778, 607, 813, 68, 253, 47, 137, 36, 582, 727, 649, 198, 295, 395, 626, 77, 58, 71, 943, 460, 103, 864, 434, 630, 377, 78, 449, 182, 348, 517, 188, 142, 585, 134, 405, 707, 721, 538, 765, 842, 511, 411, 728, 515, 276, 237, 292, 80, 824, 735, 54, 548, 886, 702, 508, 124, 61, 591, 312, 542, 932, 976, 1000, 368, 594, 41, 365, 610, 828, 668, 133, 562, 362, 310, 989, 67, 529, 435, 471, 993, 65, 939, 283, 37, 816, 602, 917, 872, 737, 722, 811, 467, 771, 366, 716, 595, 766, 980, 129, 677, 746, 625, 442, 981, 162, 904, 758, 656, 419, 17, 163, 75, 494, 452, 574, 592, 148, 501, 912, 229, 631, 600, 44, 667, 970, 817, 315, 749, 421, 20, 495, 250, 168, 57, 242, 800, 412, 557, 713, 526, 543, 603, 294, 916, 224, 235, 947, 642, 982, 541, 718, 409, 19, 944, 205, 458, 207, 82, 232, 313, 958, 647, 179, 30, 169, 282, 191, 833, 673, 520, 877, 308, 130, 213, 121, 781, 223, 821, 303, 394, 987, 249, 355, 66, 257, 353, 851, 506, 509, 40, 194, 706, 619, 302, 74, 429, 683, 621, 384, 451, 486, 347, 829, 688, 464, 957, 8, 861, 321, 158, 174, 135, 644, 956, 371, 638, 393, 290, 150, 166, 704, 218, 504, 906, 616, 200, 258, 475, 172, 532, 938, 332, 300, 92, 711, 646, 426, 488, 825, 31, 334, 921, 655, 979, 240, 278, 490, 990, 497, 42, 733, 867, 397, 280, 436, 612, 7, 505, 227, 132, 222, 444, 87, 927, 203, 774, 342, 152, 889, 170, 920, 773, 155, 187, 263, 21, 372, 785, 457, 822, 193, 770, 468, 875, 732, 717, 131, 695, 554, 3, 759, 975, 580, 316, 413, 317, 666, 81, 559, 320, 136, 469, 361, 887, 26, 298, 73, 404, 896, 675, 564, 528, 895, 522, 356, 806, 606, 450, 340, 915, 39, 93, 988, 992, 33, 246, 740, 398, 267, 268, 239, 256, 49, 455, 100, 799, 738, 354, 922, 539, 752, 403, 790, 363, 840, 478, 764, 399, 34, 899, 672, 841, 787, 747, 304, 269, 744, 345, 914, 123, 236, 180, 663, 375, 581, 931, 983, 453, 325, 209, 801, 658, 837, 52, 597, 719, 407, 991, 565, 755, 4, 756, 503, 760, 570, 598, 10, 319, 410, 45, 199, 748, 358, 935, 643, 710, 314, 22, 882, 370]\n"
     ]
    },
    {
     "ename": "MethodError",
     "evalue": "MethodError: no method matching partialsortperm!(::Array{Int64,1}, ::Array{Float64,1}, ::Int64; rev=false, alg=Base.Sort.MergeSortAlg(), initialized=false)\nClosest candidates are:\n  partialsortperm!(::AbstractArray{#s659,1} where #s659<:Integer, ::AbstractArray{T,1} where T, ::Union{Integer, OrdinalRange}; lt, by, rev, order, initialized) at sort.jl:845 got unsupported keyword argument \"alg\"",
     "output_type": "error",
     "traceback": [
      "MethodError: no method matching partialsortperm!(::Array{Int64,1}, ::Array{Float64,1}, ::Int64; rev=false, alg=Base.Sort.MergeSortAlg(), initialized=false)\nClosest candidates are:\n  partialsortperm!(::AbstractArray{#s659,1} where #s659<:Integer, ::AbstractArray{T,1} where T, ::Union{Integer, OrdinalRange}; lt, by, rev, order, initialized) at sort.jl:845 got unsupported keyword argument \"alg\"",
      "",
      "Stacktrace:",
      " [1] kwerr(::NamedTuple{(:rev, :alg, :initialized),Tuple{Bool,Base.Sort.MergeSortAlg,Bool}}, ::Function, ::Array{Int64,1}, ::Array{Float64,1}, ::Int64) at ./error.jl:157",
      " [2] partialsortperm(::Array{Float64,1}, ::Int64; kwargs::Base.Iterators.Pairs{Symbol,Any,Tuple{Symbol,Symbol},NamedTuple{(:rev, :alg),Tuple{Bool,Base.Sort.MergeSortAlg}}}) at ./sort.jl:790",
      " [3] ordenar_parcial(::Array{Float64,1}, ::Int64) at ./In[352]:10",
      " [4] top-level scope at show.jl:613",
      " [5] top-level scope at In[354]:2"
     ]
    }
   ],
   "source": [
    "@show ordenar(weights);\n",
    "@show ordenar_parcial(weights, 4);\n",
    "@show ordenar_parcial2(weights, 4);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 357,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "search: \u001b[0m\u001b[1ms\u001b[22m\u001b[0m\u001b[1mo\u001b[22m\u001b[0m\u001b[1mr\u001b[22m\u001b[0m\u001b[1mt\u001b[22m \u001b[0m\u001b[1ms\u001b[22m\u001b[0m\u001b[1mo\u001b[22m\u001b[0m\u001b[1mr\u001b[22m\u001b[0m\u001b[1mt\u001b[22m! \u001b[0m\u001b[1ms\u001b[22m\u001b[0m\u001b[1mo\u001b[22m\u001b[0m\u001b[1mr\u001b[22m\u001b[0m\u001b[1mt\u001b[22mperm \u001b[0m\u001b[1ms\u001b[22m\u001b[0m\u001b[1mo\u001b[22m\u001b[0m\u001b[1mr\u001b[22m\u001b[0m\u001b[1mt\u001b[22mperm! \u001b[0m\u001b[1ms\u001b[22m\u001b[0m\u001b[1mo\u001b[22m\u001b[0m\u001b[1mr\u001b[22m\u001b[0m\u001b[1mt\u001b[22m_gens \u001b[0m\u001b[1ms\u001b[22m\u001b[0m\u001b[1mo\u001b[22m\u001b[0m\u001b[1mr\u001b[22m\u001b[0m\u001b[1mt\u001b[22mslices \u001b[0m\u001b[1ms\u001b[22m\u001b[0m\u001b[1mo\u001b[22m\u001b[0m\u001b[1mr\u001b[22m\u001b[0m\u001b[1mt\u001b[22m_gens2\n",
      "\n"
     ]
    },
    {
     "data": {
      "text/latex": [
       "\\begin{verbatim}\n",
       "sort(v; alg::Algorithm=defalg(v), lt=isless, by=identity, rev::Bool=false, order::Ordering=Forward)\n",
       "\\end{verbatim}\n",
       "Variant of \\href{@ref}{\\texttt{sort!}} that returns a sorted copy of \\texttt{v} leaving \\texttt{v} itself unmodified.\n",
       "\n",
       "\\section{Examples}\n",
       "\\begin{verbatim}\n",
       "julia> v = [3, 1, 2];\n",
       "\n",
       "julia> sort(v)\n",
       "3-element Array{Int64,1}:\n",
       " 1\n",
       " 2\n",
       " 3\n",
       "\n",
       "julia> v\n",
       "3-element Array{Int64,1}:\n",
       " 3\n",
       " 1\n",
       " 2\n",
       "\\end{verbatim}\n",
       "\\rule{\\textwidth}{1pt}\n",
       "\\begin{verbatim}\n",
       "sort(A; dims::Integer, alg::Algorithm=DEFAULT_UNSTABLE, lt=isless, by=identity, rev::Bool=false, order::Ordering=Forward)\n",
       "\\end{verbatim}\n",
       "Sort a multidimensional array \\texttt{A} along the given dimension. See \\href{@ref}{\\texttt{sort!}} for a description of possible keyword arguments.\n",
       "\n",
       "To sort slices of an array, refer to \\href{@ref}{\\texttt{sortslices}}.\n",
       "\n",
       "\\section{Examples}\n",
       "\\begin{verbatim}\n",
       "julia> A = [4 3; 1 2]\n",
       "2×2 Array{Int64,2}:\n",
       " 4  3\n",
       " 1  2\n",
       "\n",
       "julia> sort(A, dims = 1)\n",
       "2×2 Array{Int64,2}:\n",
       " 1  2\n",
       " 4  3\n",
       "\n",
       "julia> sort(A, dims = 2)\n",
       "2×2 Array{Int64,2}:\n",
       " 3  4\n",
       " 1  2\n",
       "\\end{verbatim}\n",
       "\\rule{\\textwidth}{1pt}\n",
       "\\begin{verbatim}\n",
       "sort(parameters::Vector, order::Tuple)\n",
       "\\end{verbatim}\n",
       "Filter and sort the vector \\texttt{parameters} according to \\texttt{order}.\n",
       "\n",
       "\\subsubsection{Input}\n",
       "\\begin{itemize}\n",
       "\\item \\texttt{parameters} – vector of tuples\n",
       "\n",
       "\n",
       "\\item \\texttt{order}      – tuple of symbols\n",
       "\n",
       "\\end{itemize}\n",
       "\\subsubsection{Output}\n",
       "A new vector of tuples corresponding to \\texttt{parameters} filtered and sorted according to \\texttt{order}.\n",
       "\n",
       "\\subsubsection{Examples}\n",
       "\\begin{verbatim}\n",
       "julia> parameters= [(:U1, :U), (:X1, :X), (:W1, :W)];\n",
       "\n",
       "julia> sort(parameters, (:X, :U, :W))\n",
       "3-element Array{Tuple{Any,Symbol},1}:\n",
       " (:X1, :X)\n",
       " (:U1, :U)\n",
       " (:W1, :W)\n",
       "\n",
       "julia>  parameters = [(:const, :c), (:A, :A)];\n",
       "\n",
       "julia> sort(parameters, (:A, :B, :c, :D))\n",
       "2-element Array{Tuple{Any,Symbol},1}:\n",
       " (:A, :A)\n",
       " (:const, :c)\n",
       "\\end{verbatim}\n",
       "\\subsubsection{Notes}\n",
       "\\texttt{parameters} is a vector that contains tuples whose second element is considered for the sorting according to \\texttt{order}.\n",
       "\n",
       "If a value of \\texttt{order} is not contained in \\texttt{parameters}, the corresponding entry of \\texttt{order} will be omitted.\n",
       "\n"
      ],
      "text/markdown": [
       "```\n",
       "sort(v; alg::Algorithm=defalg(v), lt=isless, by=identity, rev::Bool=false, order::Ordering=Forward)\n",
       "```\n",
       "\n",
       "Variant of [`sort!`](@ref) that returns a sorted copy of `v` leaving `v` itself unmodified.\n",
       "\n",
       "# Examples\n",
       "\n",
       "```jldoctest\n",
       "julia> v = [3, 1, 2];\n",
       "\n",
       "julia> sort(v)\n",
       "3-element Array{Int64,1}:\n",
       " 1\n",
       " 2\n",
       " 3\n",
       "\n",
       "julia> v\n",
       "3-element Array{Int64,1}:\n",
       " 3\n",
       " 1\n",
       " 2\n",
       "```\n",
       "\n",
       "---\n",
       "\n",
       "```\n",
       "sort(A; dims::Integer, alg::Algorithm=DEFAULT_UNSTABLE, lt=isless, by=identity, rev::Bool=false, order::Ordering=Forward)\n",
       "```\n",
       "\n",
       "Sort a multidimensional array `A` along the given dimension. See [`sort!`](@ref) for a description of possible keyword arguments.\n",
       "\n",
       "To sort slices of an array, refer to [`sortslices`](@ref).\n",
       "\n",
       "# Examples\n",
       "\n",
       "```jldoctest\n",
       "julia> A = [4 3; 1 2]\n",
       "2×2 Array{Int64,2}:\n",
       " 4  3\n",
       " 1  2\n",
       "\n",
       "julia> sort(A, dims = 1)\n",
       "2×2 Array{Int64,2}:\n",
       " 1  2\n",
       " 4  3\n",
       "\n",
       "julia> sort(A, dims = 2)\n",
       "2×2 Array{Int64,2}:\n",
       " 3  4\n",
       " 1  2\n",
       "```\n",
       "\n",
       "---\n",
       "\n",
       "```\n",
       "sort(parameters::Vector, order::Tuple)\n",
       "```\n",
       "\n",
       "Filter and sort the vector `parameters` according to `order`.\n",
       "\n",
       "### Input\n",
       "\n",
       "  * `parameters` – vector of tuples\n",
       "  * `order`      – tuple of symbols\n",
       "\n",
       "### Output\n",
       "\n",
       "A new vector of tuples corresponding to `parameters` filtered and sorted according to `order`.\n",
       "\n",
       "### Examples\n",
       "\n",
       "```jldoctest\n",
       "julia> parameters= [(:U1, :U), (:X1, :X), (:W1, :W)];\n",
       "\n",
       "julia> sort(parameters, (:X, :U, :W))\n",
       "3-element Array{Tuple{Any,Symbol},1}:\n",
       " (:X1, :X)\n",
       " (:U1, :U)\n",
       " (:W1, :W)\n",
       "\n",
       "julia>  parameters = [(:const, :c), (:A, :A)];\n",
       "\n",
       "julia> sort(parameters, (:A, :B, :c, :D))\n",
       "2-element Array{Tuple{Any,Symbol},1}:\n",
       " (:A, :A)\n",
       " (:const, :c)\n",
       "```\n",
       "\n",
       "### Notes\n",
       "\n",
       "`parameters` is a vector that contains tuples whose second element is considered for the sorting according to `order`.\n",
       "\n",
       "If a value of `order` is not contained in `parameters`, the corresponding entry of `order` will be omitted.\n"
      ],
      "text/plain": [
       "\u001b[36m  sort(v; alg::Algorithm=defalg(v), lt=isless, by=identity, rev::Bool=false, order::Ordering=Forward)\u001b[39m\n",
       "\n",
       "  Variant of \u001b[36msort!\u001b[39m that returns a sorted copy of \u001b[36mv\u001b[39m leaving \u001b[36mv\u001b[39m itself\n",
       "  unmodified.\n",
       "\n",
       "\u001b[1m  Examples\u001b[22m\n",
       "\u001b[1m  ≡≡≡≡≡≡≡≡≡≡\u001b[22m\n",
       "\n",
       "\u001b[36m  julia> v = [3, 1, 2];\u001b[39m\n",
       "\u001b[36m  \u001b[39m\n",
       "\u001b[36m  julia> sort(v)\u001b[39m\n",
       "\u001b[36m  3-element Array{Int64,1}:\u001b[39m\n",
       "\u001b[36m   1\u001b[39m\n",
       "\u001b[36m   2\u001b[39m\n",
       "\u001b[36m   3\u001b[39m\n",
       "\u001b[36m  \u001b[39m\n",
       "\u001b[36m  julia> v\u001b[39m\n",
       "\u001b[36m  3-element Array{Int64,1}:\u001b[39m\n",
       "\u001b[36m   3\u001b[39m\n",
       "\u001b[36m   1\u001b[39m\n",
       "\u001b[36m   2\u001b[39m\n",
       "\n",
       "  ────────────────────────────────────────────────────────────────────────────\n",
       "\n",
       "\u001b[36m  sort(A; dims::Integer, alg::Algorithm=DEFAULT_UNSTABLE, lt=isless, by=identity, rev::Bool=false, order::Ordering=Forward)\u001b[39m\n",
       "\n",
       "  Sort a multidimensional array \u001b[36mA\u001b[39m along the given dimension. See \u001b[36msort!\u001b[39m for a\n",
       "  description of possible keyword arguments.\n",
       "\n",
       "  To sort slices of an array, refer to \u001b[36msortslices\u001b[39m.\n",
       "\n",
       "\u001b[1m  Examples\u001b[22m\n",
       "\u001b[1m  ≡≡≡≡≡≡≡≡≡≡\u001b[22m\n",
       "\n",
       "\u001b[36m  julia> A = [4 3; 1 2]\u001b[39m\n",
       "\u001b[36m  2×2 Array{Int64,2}:\u001b[39m\n",
       "\u001b[36m   4  3\u001b[39m\n",
       "\u001b[36m   1  2\u001b[39m\n",
       "\u001b[36m  \u001b[39m\n",
       "\u001b[36m  julia> sort(A, dims = 1)\u001b[39m\n",
       "\u001b[36m  2×2 Array{Int64,2}:\u001b[39m\n",
       "\u001b[36m   1  2\u001b[39m\n",
       "\u001b[36m   4  3\u001b[39m\n",
       "\u001b[36m  \u001b[39m\n",
       "\u001b[36m  julia> sort(A, dims = 2)\u001b[39m\n",
       "\u001b[36m  2×2 Array{Int64,2}:\u001b[39m\n",
       "\u001b[36m   3  4\u001b[39m\n",
       "\u001b[36m   1  2\u001b[39m\n",
       "\n",
       "  ────────────────────────────────────────────────────────────────────────────\n",
       "\n",
       "\u001b[36m  sort(parameters::Vector, order::Tuple)\u001b[39m\n",
       "\n",
       "  Filter and sort the vector \u001b[36mparameters\u001b[39m according to \u001b[36morder\u001b[39m.\n",
       "\n",
       "\u001b[1m  Input\u001b[22m\n",
       "\u001b[1m  –––––––\u001b[22m\n",
       "\n",
       "    •    \u001b[36mparameters\u001b[39m – vector of tuples\n",
       "\n",
       "    •    \u001b[36morder\u001b[39m – tuple of symbols\n",
       "\n",
       "\u001b[1m  Output\u001b[22m\n",
       "\u001b[1m  ––––––––\u001b[22m\n",
       "\n",
       "  A new vector of tuples corresponding to \u001b[36mparameters\u001b[39m filtered and sorted\n",
       "  according to \u001b[36morder\u001b[39m.\n",
       "\n",
       "\u001b[1m  Examples\u001b[22m\n",
       "\u001b[1m  ––––––––––\u001b[22m\n",
       "\n",
       "\u001b[36m  julia> parameters= [(:U1, :U), (:X1, :X), (:W1, :W)];\u001b[39m\n",
       "\u001b[36m  \u001b[39m\n",
       "\u001b[36m  julia> sort(parameters, (:X, :U, :W))\u001b[39m\n",
       "\u001b[36m  3-element Array{Tuple{Any,Symbol},1}:\u001b[39m\n",
       "\u001b[36m   (:X1, :X)\u001b[39m\n",
       "\u001b[36m   (:U1, :U)\u001b[39m\n",
       "\u001b[36m   (:W1, :W)\u001b[39m\n",
       "\u001b[36m  \u001b[39m\n",
       "\u001b[36m  julia>  parameters = [(:const, :c), (:A, :A)];\u001b[39m\n",
       "\u001b[36m  \u001b[39m\n",
       "\u001b[36m  julia> sort(parameters, (:A, :B, :c, :D))\u001b[39m\n",
       "\u001b[36m  2-element Array{Tuple{Any,Symbol},1}:\u001b[39m\n",
       "\u001b[36m   (:A, :A)\u001b[39m\n",
       "\u001b[36m   (:const, :c)\u001b[39m\n",
       "\n",
       "\u001b[1m  Notes\u001b[22m\n",
       "\u001b[1m  –––––––\u001b[22m\n",
       "\n",
       "  \u001b[36mparameters\u001b[39m is a vector that contains tuples whose second element is\n",
       "  considered for the sorting according to \u001b[36morder\u001b[39m.\n",
       "\n",
       "  If a value of \u001b[36morder\u001b[39m is not contained in \u001b[36mparameters\u001b[39m, the corresponding entry\n",
       "  of \u001b[36morder\u001b[39m will be omitted."
      ]
     },
     "execution_count": 357,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "?sort"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 351,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  21.414 μs (4 allocations: 7.98 KiB)\n",
      "  11.735 μs (3 allocations: 7.97 KiB)\n",
      "  11.869 μs (3 allocations: 7.97 KiB)\n"
     ]
    }
   ],
   "source": [
    "@btime ordenar($weights);\n",
    "@btime ordenar_parcial($weights, 4);\n",
    "@btime ordenar_parcial2($weights, 4);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 222,
   "metadata": {},
   "outputs": [],
   "source": [
    "M = rand(3, 20) .- 0.5;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 223,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  866.949 ns (26 allocations: 1.53 KiB)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "([0.5418556576014054, 0.6376114530568366, 0.6517988869384466, 0.549504845032359, 0.6920728992613562, 0.6289673404321229, 0.3681561542158848, 0.6085334801996619, 0.24222904667797338, 0.4743478377653587, 0.49198424779384486, 0.5081295421282085, 0.2905420486410917, 0.430616149318117, 0.6544704420391821, 0.6569162461795516, 0.6077802596790824, 0.5714812817550053, 0.3941155390512953, 0.5520421444336712], [5, 16, 15, 3, 2, 6, 8, 17, 18, 20, 4, 1, 12, 11, 10, 14, 19, 7, 13, 9])"
      ]
     },
     "execution_count": 223,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "weights, indices = sort_gens_partial(M, 3)\n",
    "@btime sort_gens_partial($M, 3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 224,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  606.802 ns (26 allocations: 1.48 KiB)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "([0.5418556576014054, 0.6376114530568366, 0.6517988869384466, 0.549504845032359, 0.6920728992613562, 0.6289673404321229, 0.3681561542158848, 0.6085334801996619, 0.24222904667797338, 0.4743478377653587, 0.49198424779384486, 0.5081295421282085, 0.2905420486410917, 0.430616149318117, 0.6544704420391821, 0.6569162461795516, 0.6077802596790824, 0.5714812817550053, 0.3941155390512953, 0.5520421444336712], [5, 16, 15, 3, 2, 6, 8, 17, 18, 20, 4, 1, 12, 11, 10, 14, 19, 7, 13, 9])"
      ]
     },
     "execution_count": 224,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "weights, indices = sort_gens2(M)\n",
    "@btime sort_gens2($M)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 178,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "search: \u001b[0m\u001b[1mp\u001b[22m\u001b[0m\u001b[1ma\u001b[22m\u001b[0m\u001b[1mr\u001b[22m\u001b[0m\u001b[1mt\u001b[22m\u001b[0m\u001b[1mi\u001b[22m\u001b[0m\u001b[1ma\u001b[22m\u001b[0m\u001b[1ml\u001b[22m\u001b[0m\u001b[1ms\u001b[22m\u001b[0m\u001b[1mo\u001b[22m\u001b[0m\u001b[1mr\u001b[22m\u001b[0m\u001b[1mt\u001b[22m\u001b[0m\u001b[1mp\u001b[22m\u001b[0m\u001b[1me\u001b[22m\u001b[0m\u001b[1mr\u001b[22m\u001b[0m\u001b[1mm\u001b[22m\u001b[0m\u001b[1m!\u001b[22m \u001b[0m\u001b[1mp\u001b[22m\u001b[0m\u001b[1ma\u001b[22m\u001b[0m\u001b[1mr\u001b[22m\u001b[0m\u001b[1mt\u001b[22m\u001b[0m\u001b[1mi\u001b[22m\u001b[0m\u001b[1ma\u001b[22m\u001b[0m\u001b[1ml\u001b[22m\u001b[0m\u001b[1ms\u001b[22m\u001b[0m\u001b[1mo\u001b[22m\u001b[0m\u001b[1mr\u001b[22m\u001b[0m\u001b[1mt\u001b[22m\u001b[0m\u001b[1mp\u001b[22m\u001b[0m\u001b[1me\u001b[22m\u001b[0m\u001b[1mr\u001b[22m\u001b[0m\u001b[1mm\u001b[22m \u001b[0m\u001b[1mp\u001b[22m\u001b[0m\u001b[1ma\u001b[22m\u001b[0m\u001b[1mr\u001b[22m\u001b[0m\u001b[1mt\u001b[22m\u001b[0m\u001b[1mi\u001b[22m\u001b[0m\u001b[1ma\u001b[22m\u001b[0m\u001b[1ml\u001b[22m\u001b[0m\u001b[1ms\u001b[22m\u001b[0m\u001b[1mo\u001b[22m\u001b[0m\u001b[1mr\u001b[22m\u001b[0m\u001b[1mt\u001b[22m! \u001b[0m\u001b[1mp\u001b[22m\u001b[0m\u001b[1ma\u001b[22m\u001b[0m\u001b[1mr\u001b[22m\u001b[0m\u001b[1mt\u001b[22m\u001b[0m\u001b[1mi\u001b[22m\u001b[0m\u001b[1ma\u001b[22m\u001b[0m\u001b[1ml\u001b[22m\u001b[0m\u001b[1ms\u001b[22m\u001b[0m\u001b[1mo\u001b[22m\u001b[0m\u001b[1mr\u001b[22m\u001b[0m\u001b[1mt\u001b[22m\n",
      "\n"
     ]
    },
    {
     "data": {
      "text/latex": [
       "\\begin{verbatim}\n",
       "partialsortperm!(ix, v, k; by=<transform>, lt=<comparison>, rev=false, initialized=false)\n",
       "\\end{verbatim}\n",
       "Like \\href{@ref}{\\texttt{partialsortperm}}, but accepts a preallocated index vector \\texttt{ix} the same size as \\texttt{v}, which is used to store (a permutation of) the indices of \\texttt{v}.\n",
       "\n",
       "If the index vector \\texttt{ix} is initialized with the indices of \\texttt{v} (or a permutation thereof), \\texttt{initialized} should be set to \\texttt{true}.\n",
       "\n",
       "If \\texttt{initialized} is \\texttt{false} (the default), then \\texttt{ix} is initialized to contain the indices of \\texttt{v}.\n",
       "\n",
       "If \\texttt{initialized} is \\texttt{true}, but \\texttt{ix} does not contain (a permutation of) the indices of \\texttt{v}, the behavior of \\texttt{partialsortperm!} is undefined.\n",
       "\n",
       "(Typically, the indices of \\texttt{v} will be \\texttt{1:length(v)}, although if \\texttt{v} has an alternative array type with non-one-based indices, such as an \\texttt{OffsetArray}, \\texttt{ix} must also be an \\texttt{OffsetArray} with the same indices, and must contain as values (a permutation of) these same indices.)\n",
       "\n",
       "Upon return, \\texttt{ix} is guaranteed to have the indices \\texttt{k} in their sorted positions, such that\n",
       "\n",
       "\\begin{verbatim}\n",
       "partialsortperm!(ix, v, k);\n",
       "v[ix[k]] == partialsort(v, k)\n",
       "\\end{verbatim}\n",
       "The return value is the \\texttt{k}th element of \\texttt{ix} if \\texttt{k} is an integer, or view into \\texttt{ix} if \\texttt{k} is a range.\n",
       "\n",
       "\\section{Examples}\n",
       "\\begin{verbatim}\n",
       "julia> v = [3, 1, 2, 1];\n",
       "\n",
       "julia> ix = Vector{Int}(undef, 4);\n",
       "\n",
       "julia> partialsortperm!(ix, v, 1)\n",
       "2\n",
       "\n",
       "julia> ix = [1:4;];\n",
       "\n",
       "julia> partialsortperm!(ix, v, 2:3, initialized=true)\n",
       "2-element view(::Array{Int64,1}, 2:3) with eltype Int64:\n",
       " 4\n",
       " 3\n",
       "\\end{verbatim}\n",
       "\\begin{verbatim}\n",
       "\n",
       "\\end{verbatim}\n"
      ],
      "text/markdown": [
       "```\n",
       "partialsortperm!(ix, v, k; by=<transform>, lt=<comparison>, rev=false, initialized=false)\n",
       "```\n",
       "\n",
       "Like [`partialsortperm`](@ref), but accepts a preallocated index vector `ix` the same size as `v`, which is used to store (a permutation of) the indices of `v`.\n",
       "\n",
       "If the index vector `ix` is initialized with the indices of `v` (or a permutation thereof), `initialized` should be set to `true`.\n",
       "\n",
       "If `initialized` is `false` (the default), then `ix` is initialized to contain the indices of `v`.\n",
       "\n",
       "If `initialized` is `true`, but `ix` does not contain (a permutation of) the indices of `v`, the behavior of `partialsortperm!` is undefined.\n",
       "\n",
       "(Typically, the indices of `v` will be `1:length(v)`, although if `v` has an alternative array type with non-one-based indices, such as an `OffsetArray`, `ix` must also be an `OffsetArray` with the same indices, and must contain as values (a permutation of) these same indices.)\n",
       "\n",
       "Upon return, `ix` is guaranteed to have the indices `k` in their sorted positions, such that\n",
       "\n",
       "```julia\n",
       "partialsortperm!(ix, v, k);\n",
       "v[ix[k]] == partialsort(v, k)\n",
       "```\n",
       "\n",
       "The return value is the `k`th element of `ix` if `k` is an integer, or view into `ix` if `k` is a range.\n",
       "\n",
       "# Examples\n",
       "\n",
       "```jldoctest\n",
       "julia> v = [3, 1, 2, 1];\n",
       "\n",
       "julia> ix = Vector{Int}(undef, 4);\n",
       "\n",
       "julia> partialsortperm!(ix, v, 1)\n",
       "2\n",
       "\n",
       "julia> ix = [1:4;];\n",
       "\n",
       "julia> partialsortperm!(ix, v, 2:3, initialized=true)\n",
       "2-element view(::Array{Int64,1}, 2:3) with eltype Int64:\n",
       " 4\n",
       " 3\n",
       "```\n",
       "\n",
       "```\n",
       "\n",
       "```\n"
      ],
      "text/plain": [
       "\u001b[36m  partialsortperm!(ix, v, k; by=<transform>, lt=<comparison>, rev=false, initialized=false)\u001b[39m\n",
       "\n",
       "  Like \u001b[36mpartialsortperm\u001b[39m, but accepts a preallocated index vector \u001b[36mix\u001b[39m the same\n",
       "  size as \u001b[36mv\u001b[39m, which is used to store (a permutation of) the indices of \u001b[36mv\u001b[39m.\n",
       "\n",
       "  If the index vector \u001b[36mix\u001b[39m is initialized with the indices of \u001b[36mv\u001b[39m (or a\n",
       "  permutation thereof), \u001b[36minitialized\u001b[39m should be set to \u001b[36mtrue\u001b[39m.\n",
       "\n",
       "  If \u001b[36minitialized\u001b[39m is \u001b[36mfalse\u001b[39m (the default), then \u001b[36mix\u001b[39m is initialized to contain the\n",
       "  indices of \u001b[36mv\u001b[39m.\n",
       "\n",
       "  If \u001b[36minitialized\u001b[39m is \u001b[36mtrue\u001b[39m, but \u001b[36mix\u001b[39m does not contain (a permutation of) the\n",
       "  indices of \u001b[36mv\u001b[39m, the behavior of \u001b[36mpartialsortperm!\u001b[39m is undefined.\n",
       "\n",
       "  (Typically, the indices of \u001b[36mv\u001b[39m will be \u001b[36m1:length(v)\u001b[39m, although if \u001b[36mv\u001b[39m has an\n",
       "  alternative array type with non-one-based indices, such as an \u001b[36mOffsetArray\u001b[39m,\n",
       "  \u001b[36mix\u001b[39m must also be an \u001b[36mOffsetArray\u001b[39m with the same indices, and must contain as\n",
       "  values (a permutation of) these same indices.)\n",
       "\n",
       "  Upon return, \u001b[36mix\u001b[39m is guaranteed to have the indices \u001b[36mk\u001b[39m in their sorted\n",
       "  positions, such that\n",
       "\n",
       "\u001b[36m  partialsortperm!(ix, v, k);\u001b[39m\n",
       "\u001b[36m  v[ix[k]] == partialsort(v, k)\u001b[39m\n",
       "\n",
       "  The return value is the \u001b[36mk\u001b[39mth element of \u001b[36mix\u001b[39m if \u001b[36mk\u001b[39m is an integer, or view into\n",
       "  \u001b[36mix\u001b[39m if \u001b[36mk\u001b[39m is a range.\n",
       "\n",
       "\u001b[1m  Examples\u001b[22m\n",
       "\u001b[1m  ≡≡≡≡≡≡≡≡≡≡\u001b[22m\n",
       "\n",
       "\u001b[36m  julia> v = [3, 1, 2, 1];\u001b[39m\n",
       "\u001b[36m  \u001b[39m\n",
       "\u001b[36m  julia> ix = Vector{Int}(undef, 4);\u001b[39m\n",
       "\u001b[36m  \u001b[39m\n",
       "\u001b[36m  julia> partialsortperm!(ix, v, 1)\u001b[39m\n",
       "\u001b[36m  2\u001b[39m\n",
       "\u001b[36m  \u001b[39m\n",
       "\u001b[36m  julia> ix = [1:4;];\u001b[39m\n",
       "\u001b[36m  \u001b[39m\n",
       "\u001b[36m  julia> partialsortperm!(ix, v, 2:3, initialized=true)\u001b[39m\n",
       "\u001b[36m  2-element view(::Array{Int64,1}, 2:3) with eltype Int64:\u001b[39m\n",
       "\u001b[36m   4\u001b[39m\n",
       "\u001b[36m   3\u001b[39m\n",
       "\n",
       "\u001b[36m  \u001b[39m"
      ]
     },
     "execution_count": 178,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "?partialsortperm!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 96,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  272.756 ns (7 allocations: 704 bytes)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "6-element Array{Float64,1}:\n",
       " 0.5290681285562354\n",
       " 0.496698011345726\n",
       " 0.6021596560622862\n",
       " 0.43980407256177734\n",
       " 0.3679142293918612\n",
       " 0.28480787112099604"
      ]
     },
     "execution_count": 96,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "@btime sort_gens($M)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 94,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.5290681285562354"
      ]
     },
     "execution_count": 94,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "M[1, 1]^2 + M[2, 1]^2 |> sqrt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 92,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "6-element Array{Float64,1}:\n",
       " 0.5290681285562354\n",
       " 0.496698011345726\n",
       " 0.6021596560622862\n",
       " 0.43980407256177734\n",
       " 0.3679142293918612\n",
       " 0.28480787112099604"
      ]
     },
     "execution_count": 92,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "@btime sort_gens($M)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 282,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "_reduce_order_combastel (generic function with 2 methods)"
      ]
     },
     "execution_count": 282,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# return a zonotope whose order is r (if possible), otherwise the order is the\n",
    "# closest integer (floor)\n",
    "# assume that r >= 1\n",
    "using LazySets: center\n",
    "\n",
    "function _reduce_order_combastel(Z::Zonotope{N, MN}, r::Number) where {N, MN}\n",
    "    r >= 1 || throw(ArgumentError(\"the target order should be at least 1, but it is $r\"))\n",
    "    c = Z.center\n",
    "    G = Z.generators\n",
    "    n, p = size(G)\n",
    "    \n",
    "    (r * n >= p) && return Z  # do not reduce\n",
    "\n",
    "    # compute interval hull of L\n",
    "    Lred = zeros(N, n, n)\n",
    "    @inbounds for i in 1:n\n",
    "        for j in 1:p\n",
    "            Lred[i, i] += abs(G[i, j])\n",
    "        end\n",
    "    end\n",
    "\n",
    "    if r == 1\n",
    "        return Zonotope(c, Lred)\n",
    "    end\n",
    "\n",
    "    # split Z into K and L\n",
    "    # sort generators by decreasing 2-norm\n",
    "    weights = zeros(N, p)\n",
    "    @inbounds for i in 1:p\n",
    "        weights[i] = norm(view(G, :, i), 2)\n",
    "    end\n",
    "    indices = Vector{Int}(undef, p)\n",
    "    sortperm!(indices, weights, rev=true, initialized=false)\n",
    "\n",
    "    m = floor(Int, n * (r - 1))\n",
    "    Gred = hcat(G[:, indices[1:m]], Lred)\n",
    "    return Zonotope(c, Gred)\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 262,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "10.0"
      ]
     },
     "execution_count": 262,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Z = Zonotope(rand(2), rand(2, 20) .- 0.5);\n",
    "order(Z) * 1.0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Reduction to order 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 264,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/svg+xml": [
       "<?xml version=\"1.0\" encoding=\"utf-8\"?>\n",
       "<svg xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\" width=\"600\" height=\"400\" viewBox=\"0 0 2400 1600\">\n",
       "<defs>\n",
       "  <clipPath id=\"clip4800\">\n",
       "    <rect x=\"0\" y=\"0\" width=\"2400\" height=\"1600\"/>\n",
       "  </clipPath>\n",
       "</defs>\n",
       "<path clip-path=\"url(#clip4800)\" d=\"\n",
       "M0 1600 L2400 1600 L2400 0 L0 0  Z\n",
       "  \" fill=\"#ffffff\" fill-rule=\"evenodd\" fill-opacity=\"1\"/>\n",
       "<defs>\n",
       "  <clipPath id=\"clip4801\">\n",
       "    <rect x=\"480\" y=\"0\" width=\"1681\" height=\"1600\"/>\n",
       "  </clipPath>\n",
       "</defs>\n",
       "<path clip-path=\"url(#clip4800)\" d=\"\n",
       "M141.865 1487.47 L2352.76 1487.47 L2352.76 47.2441 L141.865 47.2441  Z\n",
       "  \" fill=\"#ffffff\" fill-rule=\"evenodd\" fill-opacity=\"1\"/>\n",
       "<defs>\n",
       "  <clipPath id=\"clip4802\">\n",
       "    <rect x=\"141\" y=\"47\" width=\"2212\" height=\"1441\"/>\n",
       "  </clipPath>\n",
       "</defs>\n",
       "<polyline clip-path=\"url(#clip4802)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  327.844,1487.47 327.844,47.2441 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip4802)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  774.781,1487.47 774.781,47.2441 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip4802)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  1221.72,1487.47 1221.72,47.2441 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip4802)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  1668.66,1487.47 1668.66,47.2441 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip4802)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  2115.59,1487.47 2115.59,47.2441 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip4802)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  141.865,1449.05 2352.76,1449.05 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip4802)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  141.865,1167.3 2352.76,1167.3 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip4802)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  141.865,885.552 2352.76,885.552 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip4802)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  141.865,603.803 2352.76,603.803 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip4802)\" style=\"stroke:#000000; stroke-width:2; stroke-opacity:0.1; fill:none\" points=\"\n",
       "  141.865,322.054 2352.76,322.054 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip4800)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  141.865,1487.47 2352.76,1487.47 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip4800)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  141.865,1487.47 141.865,47.2441 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip4800)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  327.844,1487.47 327.844,1470.19 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip4800)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  774.781,1487.47 774.781,1470.19 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip4800)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  1221.72,1487.47 1221.72,1470.19 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip4800)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  1668.66,1487.47 1668.66,1470.19 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip4800)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  2115.59,1487.47 2115.59,1470.19 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip4800)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  141.865,1449.05 168.395,1449.05 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip4800)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  141.865,1167.3 168.395,1167.3 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip4800)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  141.865,885.552 168.395,885.552 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip4800)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  141.865,603.803 168.395,603.803 \n",
       "  \"/>\n",
       "<polyline clip-path=\"url(#clip4800)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:1; fill:none\" points=\"\n",
       "  141.865,322.054 168.395,322.054 \n",
       "  \"/>\n",
       "<g clip-path=\"url(#clip4800)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 327.844, 1541.47)\" x=\"327.844\" y=\"1541.47\">-4</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip4800)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 774.781, 1541.47)\" x=\"774.781\" y=\"1541.47\">-2</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip4800)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 1221.72, 1541.47)\" x=\"1221.72\" y=\"1541.47\">0</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip4800)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 1668.66, 1541.47)\" x=\"1668.66\" y=\"1541.47\">2</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip4800)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:middle;\" transform=\"rotate(0, 2115.59, 1541.47)\" x=\"2115.59\" y=\"1541.47\">4</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip4800)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:end;\" transform=\"rotate(0, 117.865, 1466.55)\" x=\"117.865\" y=\"1466.55\">-4</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip4800)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:end;\" transform=\"rotate(0, 117.865, 1184.8)\" x=\"117.865\" y=\"1184.8\">-2</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip4800)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:end;\" transform=\"rotate(0, 117.865, 903.052)\" x=\"117.865\" y=\"903.052\">0</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip4800)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:end;\" transform=\"rotate(0, 117.865, 621.303)\" x=\"117.865\" y=\"621.303\">2</text>\n",
       "</g>\n",
       "<g clip-path=\"url(#clip4800)\">\n",
       "<text style=\"fill:#000000; fill-opacity:1; font-family:Arial,Helvetica Neue,Helvetica,sans-serif; font-size:48px; text-anchor:end;\" transform=\"rotate(0, 117.865, 339.554)\" x=\"117.865\" y=\"339.554\">4</text>\n",
       "</g>\n",
       "<path clip-path=\"url(#clip4802)\" d=\"\n",
       "M204.437 625.425 L229.183 698.221 L271.497 762.022 L339.047 843.113 L450.734 965.162 L579.045 1077.07 L758.866 1196.46 L918.626 1283.29 L983.81 1315.97 L1080.01 1360.1 \n",
       "  L1234.8 1412.58 L1296.67 1422.37 L1455.36 1433.03 L1544.37 1438.86 L1740.34 1446.71 L1945.76 1337.86 L2007.1 1298.65 L2134.81 1198.1 L2278.45 1075.65 L2286.5 967.741 \n",
       "  L2290.18 909.294 L2265.44 836.498 L2223.12 772.697 L2155.57 691.605 L2043.89 569.557 L1915.58 457.652 L1735.75 338.256 L1575.99 251.429 L1510.81 218.745 L1414.61 174.617 \n",
       "  L1259.82 122.143 L1197.95 112.348 L1039.26 101.693 L950.254 95.8538 L754.276 88.0053 L548.861 196.858 L487.524 236.073 L359.806 336.615 L216.169 459.067 L208.117 566.978 \n",
       "  L204.437 625.425 L204.437 625.425  Z\n",
       "  \" fill=\"#009af9\" fill-rule=\"evenodd\" fill-opacity=\"0.5\"/>\n",
       "<polyline clip-path=\"url(#clip4802)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:0.5; fill:none\" points=\"\n",
       "  204.437,625.425 229.183,698.221 271.497,762.022 339.047,843.113 450.734,965.162 579.045,1077.07 758.866,1196.46 918.626,1283.29 983.81,1315.97 1080.01,1360.1 \n",
       "  1234.8,1412.58 1296.67,1422.37 1455.36,1433.03 1544.37,1438.86 1740.34,1446.71 1945.76,1337.86 2007.1,1298.65 2134.81,1198.1 2278.45,1075.65 2286.5,967.741 \n",
       "  2290.18,909.294 2265.44,836.498 2223.12,772.697 2155.57,691.605 2043.89,569.557 1915.58,457.652 1735.75,338.256 1575.99,251.429 1510.81,218.745 1414.61,174.617 \n",
       "  1259.82,122.143 1197.95,112.348 1039.26,101.693 950.254,95.8538 754.276,88.0053 548.861,196.858 487.524,236.073 359.806,336.615 216.169,459.067 208.117,566.978 \n",
       "  204.437,625.425 \n",
       "  \"/>\n",
       "<path clip-path=\"url(#clip4802)\" d=\"\n",
       "M2290.18 88.0053 L204.437 88.0053 L204.437 1446.71 L2290.18 1446.71 L2290.18 88.0053 L2290.18 88.0053  Z\n",
       "  \" fill=\"#e26f46\" fill-rule=\"evenodd\" fill-opacity=\"0.5\"/>\n",
       "<polyline clip-path=\"url(#clip4802)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:0.5; fill:none\" points=\"\n",
       "  2290.18,88.0053 204.437,88.0053 204.437,1446.71 2290.18,1446.71 2290.18,88.0053 \n",
       "  \"/>\n",
       "<path clip-path=\"url(#clip4802)\" d=\"\n",
       "M2290.18 88.0053 L204.437 88.0053 L204.437 1446.71 L2290.18 1446.71 L2290.18 88.0053 L2290.18 88.0053  Z\n",
       "  \" fill=\"#3da44d\" fill-rule=\"evenodd\" fill-opacity=\"0.5\"/>\n",
       "<polyline clip-path=\"url(#clip4802)\" style=\"stroke:#000000; stroke-width:4; stroke-opacity:0.5; fill:none\" points=\"\n",
       "  2290.18,88.0053 204.437,88.0053 204.437,1446.71 2290.18,1446.71 2290.18,88.0053 \n",
       "  \"/>\n",
       "</svg>\n"
      ]
     },
     "execution_count": 264,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "plot(Z)\n",
    "plot!(_reduce_order_combastel(Z, 1))\n",
    "plot!(LazySets.reduce_order(Z, 1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 265,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  62.446 ns (2 allocations: 144 bytes)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "Zonotope{Float64,Array{Float64,1},Array{Float64,2}}([0.11452054536894174, 0.8389944773210667], [4.666755122634283 0.0; 0.0 4.822403971457481])"
      ]
     },
     "execution_count": 265,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "@btime _reduce_order_combastel($Z, 1) # reduccion a orden 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 266,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  2.500 μs (60 allocations: 5.77 KiB)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "Zonotope{Float64,Array{Float64,1},LinearAlgebra.Diagonal{Float64,Array{Float64,1}}}([0.11452054536894174, 0.8389944773210667], [4.666755122634284 0.0; 0.0 4.822403971457481])"
      ]
     },
     "execution_count": 266,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "@btime LazySets.reduce_order($Z, 1) # reduccion a orden 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 267,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "40.322580645161295"
      ]
     },
     "execution_count": 267,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "2.5e-6 / 62e-9 # 40x speedup !"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Reduction to other orders"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 268,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "10//1"
      ]
     },
     "execution_count": 268,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "order(Z)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 272,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(2, 20)"
      ]
     },
     "execution_count": 272,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "size(Z.generators)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "10 = 20 / 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Quiero 6 generadores (o menos) => orden 3 (o menos)\n",
    "- Lred tiene 2 columnas\n",
    "- Entonces, K tiene que tener m = 4 columnas (o menos)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 289,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4//1"
      ]
     },
     "execution_count": 289,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "_reduce_order_combastel(Z, 4) |> order"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 290,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  829.360 ns (33 allocations: 2.19 KiB)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "Zonotope{Float64,Array{Float64,1},Array{Float64,2}}([0.11452054536894174, 0.8389944773210667], [0.4596068963239488 0.40234086310604833 … 4.666755122634283 0.0; 0.386344382438087 -0.42376795343057627 … 0.0 4.822403971457481])"
      ]
     },
     "execution_count": 290,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "@btime _reduce_order_combastel($Z, 4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 292,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3.9179734620024127"
      ]
     },
     "execution_count": 292,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "3.248e-6 / 829e-9"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 291,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  3.248 μs (94 allocations: 8.83 KiB)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "Zonotope{Float64,Array{Float64,1},SparseArrays.SparseMatrixCSC{Float64,Int64}}([0.11452054536894174, 0.8389944773210667], \n",
       "  [1, 1]  =  0.285763\n",
       "  [2, 1]  =  0.356851\n",
       "  [1, 2]  =  -0.287089\n",
       "  [2, 2]  =  0.397177\n",
       "  [1, 3]  =  -0.357456\n",
       "  [2, 3]  =  0.308171\n",
       "  [1, 4]  =  0.321379\n",
       "  [2, 4]  =  0.434615\n",
       "  [1, 5]  =  0.459607\n",
       "  [2, 5]  =  0.386344\n",
       "  [1, 6]  =  0.402341\n",
       "  [2, 6]  =  -0.423768\n",
       "  [1, 7]  =  2.55312\n",
       "  [2, 8]  =  2.51548)"
      ]
     },
     "execution_count": 291,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "@btime LazySets.reduce_order($Z, 4);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Using LoopVectorization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 330,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "_reduce_order_combastel_avx (generic function with 1 method)"
      ]
     },
     "execution_count": 330,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "using LoopVectorization\n",
    "\n",
    "function _reduce_order_combastel_avx(Z::Zonotope{N, MN}, r::Number) where {N, MN}\n",
    "    r >= 1 || throw(ArgumentError(\"the target order should be at least 1, but it is $r\"))\n",
    "    c = Z.center\n",
    "    G = Z.generators\n",
    "    n, p = size(G)\n",
    "    \n",
    "    (r * n >= p) && return Z  # do not reduce\n",
    "\n",
    "    # compute interval hull of L\n",
    "    Lred = zeros(N, n, n)\n",
    "    @avx for i in 1:n\n",
    "        for j in 1:p\n",
    "            Lred[i, i] += abs(G[i, j])\n",
    "        end\n",
    "    end\n",
    "\n",
    "    if r == 1\n",
    "        return Zonotope(c, Lred)\n",
    "    end\n",
    "\n",
    "    # split Z into K and L\n",
    "    # sort generators by decreasing 2-norm\n",
    "    weights = zeros(N, p)\n",
    "    @avx for i in 1:p\n",
    "        weights[i] = norm(G[:, i], 2)\n",
    "    end\n",
    "    indices = Vector{Int}(undef, p)\n",
    "    sortperm!(indices, weights, rev=true, initialized=false)\n",
    "\n",
    "    m = floor(Int, n * (r - 1))\n",
    "    Gred = hcat(G[:, indices[1:m]], Lred)\n",
    "    return Zonotope(c, Gred)\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 322,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "10.0"
      ]
     },
     "execution_count": 322,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Z = Zonotope(rand(2), rand(2, 20) .- 0.5);\n",
    "order(Z) * 1.0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 324,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  61.238 ns (2 allocations: 144 bytes)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "Zonotope{Float64,Array{Float64,1},Array{Float64,2}}([0.7185079816644413, 0.6034170985427747], [5.22933020283825 0.0; 0.0 4.279036914578358])"
      ]
     },
     "execution_count": 324,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "@btime _reduce_order_combastel($Z, 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 325,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  55.817 ns (2 allocations: 144 bytes)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "Zonotope{Float64,Array{Float64,1},Array{Float64,2}}([0.7185079816644413, 0.6034170985427747], [5.22933020283825 0.0; 0.0 4.279036914578358])"
      ]
     },
     "execution_count": 325,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "@btime _reduce_order_combastel_avx($Z, 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 326,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  784.286 ns (33 allocations: 2.03 KiB)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "Zonotope{Float64,Array{Float64,1},Array{Float64,2}}([0.7185079816644413, 0.6034170985427747], [-0.2926313230153139 -0.3210852113390412 5.22933020283825 0.0; -0.48850751712190177 -0.4319938110264734 0.0 4.279036914578358])"
      ]
     },
     "execution_count": 326,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "@btime _reduce_order_combastel($Z, 2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 332,
   "metadata": {},
   "outputs": [
    {
     "ename": "MethodError",
     "evalue": "MethodError: no method matching -(::Colon, ::Int64)\nClosest candidates are:\n  -(!Matched::Complex{Bool}, ::Real) at complex.jl:307\n  -(!Matched::Missing, ::Number) at missing.jl:115\n  -(!Matched::Base.CoreLogging.LogLevel, ::Integer) at logging.jl:108\n  ...",
     "output_type": "error",
     "traceback": [
      "MethodError: no method matching -(::Colon, ::Int64)\nClosest candidates are:\n  -(!Matched::Complex{Bool}, ::Real) at complex.jl:307\n  -(!Matched::Missing, ::Number) at missing.jl:115\n  -(!Matched::Base.CoreLogging.LogLevel, ::Integer) at logging.jl:108\n  ...",
      "",
      "Stacktrace:",
      " [1] _reduce_order_combastel_avx(::Zonotope{Float64,Array{Float64,1},Array{Float64,2}}, ::Int64) at ./In[330]:26",
      " [2] top-level scope at In[332]:1"
     ]
    }
   ],
   "source": [
    "_reduce_order_combastel_avx(Z, 2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 335,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "search: \u001b[0m\u001b[1mM\u001b[22m\u001b[0m\u001b[1me\u001b[22m\u001b[0m\u001b[1mr\u001b[22m\u001b[0m\u001b[1mg\u001b[22m\u001b[0m\u001b[1me\u001b[22m\u001b[0m\u001b[1mS\u001b[22m\u001b[0m\u001b[1mo\u001b[22m\u001b[0m\u001b[1mr\u001b[22m\u001b[0m\u001b[1mt\u001b[22m\n",
      "\n"
     ]
    },
    {
     "data": {
      "text/latex": [
       "\\begin{verbatim}\n",
       "MergeSort\n",
       "\\end{verbatim}\n",
       "Indicate that a sorting function should use the merge sort algorithm. Merge sort divides the collection into subcollections and repeatedly merges them, sorting each subcollection at each step, until the entire collection has been recombined in sorted form.\n",
       "\n",
       "Characteristics:\n",
       "\n",
       "\\begin{itemize}\n",
       "\\item \\emph{stable}: preserves the ordering of elements which compare equal (e.g. \"a\" and \"A\" in a sort of letters which ignores case).\n",
       "\n",
       "\n",
       "\\item \\emph{not in-place} in memory.\n",
       "\n",
       "\n",
       "\\item \\emph{divide-and-conquer} sort strategy.\n",
       "\n",
       "\\end{itemize}\n"
      ],
      "text/markdown": [
       "```\n",
       "MergeSort\n",
       "```\n",
       "\n",
       "Indicate that a sorting function should use the merge sort algorithm. Merge sort divides the collection into subcollections and repeatedly merges them, sorting each subcollection at each step, until the entire collection has been recombined in sorted form.\n",
       "\n",
       "Characteristics:\n",
       "\n",
       "  * *stable*: preserves the ordering of elements which compare equal (e.g. \"a\" and \"A\" in a sort of letters which ignores case).\n",
       "  * *not in-place* in memory.\n",
       "  * *divide-and-conquer* sort strategy.\n"
      ],
      "text/plain": [
       "\u001b[36m  MergeSort\u001b[39m\n",
       "\n",
       "  Indicate that a sorting function should use the merge sort algorithm. Merge\n",
       "  sort divides the collection into subcollections and repeatedly merges them,\n",
       "  sorting each subcollection at each step, until the entire collection has\n",
       "  been recombined in sorted form.\n",
       "\n",
       "  Characteristics:\n",
       "\n",
       "    •    \u001b[4mstable\u001b[24m: preserves the ordering of elements which compare equal\n",
       "        (e.g. \"a\" and \"A\" in a sort of letters which ignores case).\n",
       "\n",
       "    •    \u001b[4mnot in-place\u001b[24m in memory.\n",
       "\n",
       "    •    \u001b[4mdivide-and-conquer\u001b[24m sort strategy."
      ]
     },
     "execution_count": 335,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "?MergeSort"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Method 2 (Girard 2005)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Reference implementation (LazySets)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "function reduce_order_LazySets(Z::Zonotope{N}, ::Type{<:Zonotope}, r::Union{Integer, Rational}) where {N<:Real}\n",
    "    c, G = Z.center, Z.generators\n",
    "    d, p = dim(Z), ngens(Z)\n",
    "\n",
    "    if r * d >= p || r < 1\n",
    "        # do not reduce\n",
    "        return Z\n",
    "    end\n",
    "\n",
    "    h = zeros(N, p)\n",
    "    for i in 1:p\n",
    "        h[i] = norm(G[:, i], 1) - norm(G[:, i], Inf)\n",
    "    end\n",
    "    ind = sortperm(h)\n",
    "\n",
    "    m = p - floor(Int, d * (r - 1)) # subset of ngens that are reduced\n",
    "    rg = G[:, ind[1:m]] # reduced generators\n",
    "\n",
    "    # interval hull computation of reduced generators\n",
    "    Gbox = Diagonal(sum(abs.(rg), dims=2)[:])\n",
    "    if m < p\n",
    "        Gnotred = G[:, ind[m+1:end]]\n",
    "        Gred = [Gnotred Gbox]\n",
    "    else\n",
    "        Gred = Gbox\n",
    "    end\n",
    "    return Zonotope(c, Gred)\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Reference implementation (ReachabilityAnalysis)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "function reduce_order_RA(Z::Zonotope{N, VN, MN}, r::Union{Integer, Rational}) where {N<:Real, VN, MN}\n",
    "    c = Z.center\n",
    "    G = Z.generators\n",
    "    d, p = size(G)\n",
    "\n",
    "    if (r * d >= p) || r < 1\n",
    "        # do not reduce\n",
    "        return Z\n",
    "    end\n",
    "\n",
    "    # subset of ngens that are reduced\n",
    "    m = p - floor(Int, d * (r - 1))\n",
    "\n",
    "    h = zeros(N, p)\n",
    "    @inbounds for i in 1:p\n",
    "        Gi = view(G, :, i)\n",
    "        h[i] = norm(Gi, 1) - norm(Gi, Inf)\n",
    "    end\n",
    "    ind = sortperm(h)\n",
    "\n",
    "    Gbox = zeros(N, d, p)\n",
    "    for j in ind[1:m]\n",
    "        for i in 1:d\n",
    "            @inbounds Gbox[i, j] += abs(G[i, j])\n",
    "        end\n",
    "    end\n",
    "\n",
    "    if m < p\n",
    "        Gnotred = G[:, ind[m+1:end]]\n",
    "        Gred = hcat(Gnotred, Gbox)\n",
    "    else\n",
    "        Gred = Gbox\n",
    "    end\n",
    "    return Zonotope(c, Gred)\n",
    "end"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 1.4.0",
   "language": "julia",
   "name": "julia-1.4"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.4.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}