{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Rcpp入門 その3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "「オリジナルのR」、「RcppArmadillo」、「RcppEigen」で処理速度の比較をしてみたいと思います。\n",
    "\n",
    "今回はソートで比較してみます。\n",
    "\n",
    "「オリジナルのR」は`order()`、「Armadillo」は`sort_index()`、「Eigen」は`std::sort()`を使います。\n",
    "\n",
    "※ バージョンによっては正常に動作しない場合があります。私が試した環境は、\n",
    "\n",
    "* R version 3.4.1\n",
    "* Rcpp 0.12.12\n",
    "* RcppArmadillo 0.7.960.1.2\n",
    "* RcppEigen 0.3.3.3.0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##  Armadilloは、sort_index()で実装"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[1] '0.12.12'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "[1] '0.7.960.1.2'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# パッケージの読込み\n",
    "library(Rcpp)\n",
    "library(RcppArmadillo)\n",
    "\n",
    "# パッケージのバージョンの確認\n",
    "packageVersion(\"Rcpp\")\n",
    "packageVersion(\"RcppArmadillo\")\n",
    "\n",
    "# C+11にてコンパイルします\n",
    "Sys.setenv(\"PKG_CXXFLAGS\"=\"-std=c++11\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# C++のコード\n",
    "code <- '\n",
    "arma::uvec IndexSortArma(arma::vec src)\n",
    "{\n",
    "    // ソート\n",
    "    arma::uvec index = sort_index(src);\n",
    "\n",
    "    return index + 1; //Rのインデックスは1スタートなので、1追加\n",
    "} \n",
    "'\n",
    "\n",
    "# コンパイルとリンクを行います（少し時間がかかります）\n",
    "cppIndexSortArma <- cppFunction(code, depends=\"RcppArmadillo\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Eigenは、std::sort()で実装"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Eigenにはソート関数が無いようなので、std::sort()を使います。\n",
    "\n",
    "[std::vectorをインデックスをそのままにラムダ式でsortする](http://qiita.com/Nasupl_r/items/35a903a4a69932c26c83)を参考にしました。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\n",
      "Attaching package: 'RcppEigen'\n",
      "\n",
      "The following objects are masked from 'package:RcppArmadillo':\n",
      "\n",
      "    fastLm, fastLmPure\n",
      "\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "[1] '0.12.12'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "[1] '0.3.3.3.0'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# パッケージの読込み\n",
    "library(Rcpp)\n",
    "library(RcppEigen)\n",
    "\n",
    "# パッケージのバージョンの確認\n",
    "packageVersion(\"Rcpp\")\n",
    "packageVersion(\"RcppEigen\")\n",
    "\n",
    "# C+11にてコンパイルします\n",
    "Sys.setenv(\"PKG_CXXFLAGS\"=\"-std=c++11\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# C++のコード\n",
    "code <- '\n",
    "Eigen::VectorXi IndexSortEigen(Eigen::VectorXd src)\n",
    "{\n",
    "    // 出力するインデックス\n",
    "    int num = src.size();\n",
    "    Eigen::VectorXi index(num);\n",
    "    std::iota(index.data(), index.data() + num, 0);\n",
    "\n",
    "    // ソート\n",
    "    std::sort(\n",
    "        index.data(),\n",
    "        index.data() + num,\n",
    "        [&src](double x, double y) -> bool { return src[x] < src[y]; }\n",
    "    );\n",
    "\n",
    "    return index.array() + 1; //Rのインデックスは1スタートなので、1追加\n",
    "}\n",
    "'\n",
    "\n",
    "# コンパイルとリンクを行います（少し時間がかかります）\n",
    "cppIndexSortEigen <- cppFunction(code, depends=\"RcppEigen\", includes=c('#include <RcppEigen.h>'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 結果"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "元データを作成します。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "src <- rnorm(10000000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**「オリジナルのR」の`order()`**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "   user  system elapsed \n",
      "   0.57    0.05    0.62 \n"
     ]
    }
   ],
   "source": [
    "t <- proc.time()\n",
    "\n",
    "dst_r <- order(src)\n",
    "\n",
    "print(proc.time() - t)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**「Armadillo」の`sort_index()`**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "   user  system elapsed \n",
      "   1.03    0.05    1.08 \n"
     ]
    }
   ],
   "source": [
    "t <- proc.time()\n",
    "\n",
    "dst_arma <- cppIndexSortArma(src)\n",
    "\n",
    "print(proc.time() - t)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**「Eigen」については、`std::sort()`**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "   user  system elapsed \n",
      "   2.15    0.03    2.18 \n"
     ]
    }
   ],
   "source": [
    "t <- proc.time()\n",
    "\n",
    "dst_eigen <- cppIndexSortEigen(src)\n",
    "\n",
    "print(proc.time() - t)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**結果のチェック**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1] \"index_sort()の結果チェック\"\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "TRUE"
      ],
      "text/latex": [
       "TRUE"
      ],
      "text/markdown": [
       "TRUE"
      ],
      "text/plain": [
       "[1] TRUE"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1] \"std::sort()の結果チェック\"\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "TRUE"
      ],
      "text/latex": [
       "TRUE"
      ],
      "text/markdown": [
       "TRUE"
      ],
      "text/plain": [
       "[1] TRUE"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "print(\"index_sort()の結果チェック\")\n",
    "all(dst_r == dst_arma)\n",
    "print(\"std::sort()の結果チェック\")\n",
    "all(dst_r == dst_eigen)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "今回の結果としては、以下の順位になりました。\n",
    "\n",
    "1. 「オリジナルのR」の`order()` \n",
    "1. 「Armadillo」の`sort_index()` \n",
    "1. 「Eigen」については、`std::sort()`\n",
    "\n",
    "「Armadillo」や「Eigen」については、この結果でライブラリの性能についての良し悪しは言えないと思います。「Eigen」にいたっては、`std::sort()`を使っていますし。\n",
    "\n",
    "ただし、以下のことは言えると思っています。\n",
    "\n",
    "* 「オリジナルのR」の標準の関数は結構速い\n",
    "* 「Armadillo」や「Eigen」を使いさえすれば高速化できるというわけではない"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "----\n",
    "[ホームへ](https://waku-take-a.github.io/index.html)  \n",
    "[↑Rの記事Topへ](https://waku-take-a.github.io/R.html)  \n",
    "\n",
    "※  \n",
    "上記のipynbファイル等は、[こちらに](https://github.com/WAKU-TAKE-A/RTips)あります。  \n",
    "解凍後、`The jupyter notebook`にアップロードすれば、ローカルな環境で実行することも、編集することもできます。\n",
    "\n",
    "<(_ _)>"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "3.4.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}