{
  "cells": [
    {
      "metadata": {},
      "cell_type": "markdown",
      "source": "# 変数を適切な型で初期化することについて\n\n黒木玄\n\n2017-01-09\n\nJulia言語について「トップレベルにコードを書かずに, 函数の中に書かないと速くならない」という事実はすでに有名だと思う.\n\n**(レベル1) 高速で計算して欲しいコードは必ず函数の中に入れる.**\n\nJulia言語ではこれだけは絶対に守らないとまずい.\n\nその次に注意するべきことは, 変数を適切な型で初期化することです.\n\n例えば, 「計算結果が浮動小数点数になる予定の変数 `s` を `s = 0` と整数型の 0 で初期化する」のはまずいです.  以下は高速に計算したいときに適用できる一般的な原則だと思います.\n\n**一般原則: 計算結果がT型になる予定の変数を別の型で初期化してしまうことは避ける.**\n\nより一般には, \n\n**一般原則: 変数の型が計算の途中で変化してしまうように見える書き方は避ける.**\n\nアルゴリズムの設計がクリアでないと, 変数の型が計算の途中で変化してしまうように見える書き方をしてしまいがちだと思う.  そのようなアルゴリズムの実装の仕方をして計算が速くなるはずがない.  これはJulia言語に限らない一般原則だと思います.\n\nJulia言語を使っても計算があまり速くならない場合には, Julia言語の問題だと思うまえに, 「アルゴリズムの認識が不明瞭なせいで変数の型が途中で変化してしまうような書き方をしていないかどうか」を確認するべきだと思います.\n\n問題はJulia言語におけるその問題の回避法.  最初のうちは, \n\n**(レベル2-1) 計算結果が浮動小数点数になる変数 `s` は `s = 0.0` のように浮動小数点数で初期化する.**\n\n**(レベル2-2) 計算結果が浮動小数点数の長さ2の1次元配列になる予定の変数 `s` は `s = zeros(2)` とか, `s = [1.0, 0.0]` とか, `s = rand(2)` のように浮動小数点数の配列で初期化する.**\n\nというシンプルな処方箋に従うだけでもよいと思います. もちろん, 浮動小数点数型以外についても同様です.\n\nこのようなシンプルな処方箋の利点は **楽なこと** です.  「整数型と浮動小数点型の区別」というようなプログラミングの常識をJulia言語でもきちんと遂行するだけでよい.  **(レベル2)** までなら, Julia言語に特化した特別な頭の使い方は必要ありません.  **クリアに設計されたアルゴリズムをJulia言語に翻訳するだけで高速に計算したい** と思っている人はひとまず **(レベル2)** まで注意しておけば十分な場合が多いと思います.   **(レベル2)** まで注意しておけば, アルゴリズム記述用の擬似コードのようなコードを書くだけで, Julia言語は高速に計算してくれます.\n\nしかし, この処方箋を使い続けて行くと, 次のような不満が出て来ます.\n\n**不満: `s = 0.0` のように浮動小数点数型に決め打ちして変数を初期化してしまうと, 浮動小数点数の計算しか意図した通りにできない函数ができあがってしまう.**\n\n整数型であろうが, 浮動小数点型であろうが, 有理数型であろうが, どの数値型でも通用するような函数を書きたければ, そのように書く必要があります.  Julia言語の世界ではそのように函数を書くことが慣習化しています.\n\n例えば, 変数 `s` の型が引数として与えられる配列 `x` の要素の型と同じになって欲しい場合には `s = zero(eltype(x))` のように変数 `s` を初期化する.  変数 `s` の型が変数 `a` の型と同じになって欲しい場合には `s = zero(a)` とか `s = similar(a)` のように初期化すればよい.  このような書き方はJulia言語特有の書き方であり, Julia言語が強力な言語である理由の一つになっています.\n\n**(レベル3) 変数の初期化を `s = 0.0` のような書き方で浮動小数点型で決め打ちしたりせずに, Julia言語特有のスタイルで, 函数の引数の型から変数の型が適切に決まるような書き方をする.**\n\nこのレベル3まで到達すると, 引数の型が変化しても適切かつ高速に動作する函数を書けるようになります. そのような函数は汎用性が高くなり, 非常に便利です.\n\nこの次のレベルはJulia言語の特徴の一つである **multiple dispatch** を使い始めることだと思う.\n\n**(レベル4) multiple dispatch を使い始める.**\n\nJulia言語では, 同じ名前で引数の型や個数が違う別の函数を書くことができます.\n\n例えば, 特殊函数 $f(x)$ を `x` の型が様々な精度の浮動小数点数と複素数であっても適切に計算してくれるアルゴリズムで実装したとします:\n\n```Julia\nfunction f(x)\n    xがどの型でも適切に計算してくれる一般的な計算法がここに書いてある\nend\n```\n\nしかし一般的な計算法は最適化が十分ではなく, 遅いことが多い.  例えば Float64 型の場合専用に最適化された函数を\n\n```Julia\nfunction f(x::Float64)\n    Float64型専用に最適化された計算法がここに書いてある\nend\n```\n\nと別に同じ名前 `f` で定義すると, 函数 `f` を使っている他の部分で Float64 型の数値 `a` に対する `f(a)` が実行されるときには, Float64 型専用に最適化された `f` の方を使って計算してくれるようになります.  このように, multiple dispatch を使えば, 「特定の場合」に「汎用的な部品」を「最適化された特殊な部品」に自動的に置き換えることによって, 動作を改善することが可能になります.\n\n実際には, その「特定の場合」は「`x` が Float64 型の場合」というような単純な場合だけではなく, 「`x` が `T = Float64, Float32, Float16` に対する `T` 型または `Complex{T}` 型の場合」というような少し込み入った場合も出て来るでしょう. そのような複数の型を持つ引数 `x`  を持つ函数を書くときには **(レベル3)** の処方箋が必要になります.\n\n変数を引数の型が変化しても適切な型で変数を初期化するように書くスタイルと, multiple dispatch でシステム全体を効率的に構成しようとすることは表裏一体だと考えられます.\n\n私が理解しているのは大体以上のようなこと.\n\nまだ全然Julia言語のことは理解していません. \n\nしかし, **(レベル2)** までの使い方ですでにJulia言語は高速でとても便利なプログラミング言語だと私は思っています.  \n\nそして, 単に計算が速いだけではないJulia言語特有の面白さを知るためには **(レベル3)** と **(レベル4)** まで試してみた方がよいと思います."
    },
    {
      "metadata": {
        "trusted": true
      },
      "cell_type": "code",
      "source": "using BenchmarkTools",
      "execution_count": 1,
      "outputs": []
    },
    {
      "metadata": {
        "trusted": true
      },
      "cell_type": "code",
      "source": "function mysum_int(x)\n    s = 0 # sを整数0で初期化\n    for i in 1:endof(x)\n        s += x[i]\n    end\n    s\nend\n\nfunction mysum_float(x)\n    s = 0.0 # sを浮動小数点数0.0で初期化\n    for i in 1:endof(x)\n        s += x[i]\n    end\n    s\nend\n\nfunction mysum_eltype(x)\n    s = zero(eltype(x)) # sをxの要素の型のゼロで初期化\n    for i in 1:endof(x)\n        s += x[i]\n    end\n    s\nend\n\nfunction mysum_inbounds_simd(x)\n    s = zero(eltype(x)) # sをxの要素の型のゼロで初期化\n    @inbounds @simd for i in 1:endof(x)\n        s += x[i]\n    end\n    s\nend",
      "execution_count": 2,
      "outputs": [
        {
          "output_type": "execute_result",
          "execution_count": 2,
          "data": {
            "text/plain": "mysum_inbounds_simd (generic function with 1 method)"
          },
          "metadata": {}
        }
      ]
    },
    {
      "metadata": {
        "trusted": true
      },
      "cell_type": "code",
      "source": "# 非常に遅い\n\nsrand(2018)\nx = randn(10^8)\n@benchmark mysum_int(x)",
      "execution_count": 3,
      "outputs": [
        {
          "output_type": "execute_result",
          "execution_count": 3,
          "data": {
            "text/plain": "BenchmarkTools.Trial: \n  memory estimate:  4.47 GiB\n  allocs estimate:  300000000\n  --------------\n  minimum time:     2.405 s (3.95% GC)\n  median time:      2.473 s (3.97% GC)\n  mean time:        2.456 s (3.96% GC)\n  maximum time:     2.490 s (3.97% GC)\n  --------------\n  samples:          3\n  evals/sample:     1"
          },
          "metadata": {}
        }
      ]
    },
    {
      "metadata": {
        "trusted": true
      },
      "cell_type": "code",
      "source": "# 速い (s=0.0と初期化、Juliaを使い始めた直後の人はこの速さを目標にするとよいと思う)\n\nsrand(2018)\nx = randn(10^8)\n@benchmark mysum_float(x)",
      "execution_count": 4,
      "outputs": [
        {
          "output_type": "execute_result",
          "execution_count": 4,
          "data": {
            "text/plain": "BenchmarkTools.Trial: \n  memory estimate:  16 bytes\n  allocs estimate:  1\n  --------------\n  minimum time:     97.895 ms (0.00% GC)\n  median time:      103.833 ms (0.00% GC)\n  mean time:        103.576 ms (0.00% GC)\n  maximum time:     106.866 ms (0.00% GC)\n  --------------\n  samples:          49\n  evals/sample:     1"
          },
          "metadata": {}
        }
      ]
    },
    {
      "metadata": {
        "trusted": true
      },
      "cell_type": "code",
      "source": "# 速い (Juliaをある程度理解すると s=zero(eltype(x)) のような書き方も覚えるとよい)\n\nsrand(2018)\nx = randn(10^8)\n@benchmark mysum_eltype(x)",
      "execution_count": 5,
      "outputs": [
        {
          "output_type": "execute_result",
          "execution_count": 5,
          "data": {
            "text/plain": "BenchmarkTools.Trial: \n  memory estimate:  16 bytes\n  allocs estimate:  1\n  --------------\n  minimum time:     99.357 ms (0.00% GC)\n  median time:      104.954 ms (0.00% GC)\n  mean time:        104.576 ms (0.00% GC)\n  maximum time:     107.798 ms (0.00% GC)\n  --------------\n  samples:          48\n  evals/sample:     1"
          },
          "metadata": {}
        }
      ]
    },
    {
      "metadata": {
        "trusted": true
      },
      "cell_type": "code",
      "source": "# @inbounds @simd マクロを使うと非常に速い\n\nsrand(2018)\nx = randn(10^8)\n@benchmark mysum_inbounds_simd(x)",
      "execution_count": 6,
      "outputs": [
        {
          "output_type": "execute_result",
          "execution_count": 6,
          "data": {
            "text/plain": "BenchmarkTools.Trial: \n  memory estimate:  16 bytes\n  allocs estimate:  1\n  --------------\n  minimum time:     57.616 ms (0.00% GC)\n  median time:      61.413 ms (0.00% GC)\n  mean time:        61.208 ms (0.00% GC)\n  maximum time:     64.041 ms (0.00% GC)\n  --------------\n  samples:          82\n  evals/sample:     1"
          },
          "metadata": {}
        }
      ]
    },
    {
      "metadata": {
        "trusted": true
      },
      "cell_type": "code",
      "source": "# 組み込み函数の sum は非常に速い\n\nsrand(2018)\nx = randn(10^8)\n@benchmark sum(x)",
      "execution_count": 7,
      "outputs": [
        {
          "output_type": "execute_result",
          "execution_count": 7,
          "data": {
            "text/plain": "BenchmarkTools.Trial: \n  memory estimate:  16 bytes\n  allocs estimate:  1\n  --------------\n  minimum time:     56.367 ms (0.00% GC)\n  median time:      58.750 ms (0.00% GC)\n  mean time:        61.135 ms (0.00% GC)\n  maximum time:     90.089 ms (0.00% GC)\n  --------------\n  samples:          82\n  evals/sample:     1"
          },
          "metadata": {}
        }
      ]
    },
    {
      "metadata": {
        "trusted": true
      },
      "cell_type": "code",
      "source": "# x=randn(10^8) のとき mysum_int(x) が非常に遅い理由\n#\n# 赤字で s::Union{Float64, Int64} の部分が強調されている.\n# 変数 s を整数の 0 で初期化しようとしているせいで, 引数 x に浮動小数点数の配列を与えると, \n# Juliaは s の型を一意に推定し切れずに「整数または浮動小数点数」だとしてしまっている.\n# これが計算が遅くなる理由である.\n#\n# 和 s が浮動小数点数になると期待されるケースでは s=0.0 と初期化すれば速くなる(mysum_float).\n#\n# しかし、それだと引数配列 x の要素が浮動小数点の場合にのみ使用が適切な函数になってしまう.\n# 引数 x の要素の型ごとに適切な方法で和を取るようにするためには s = zero(eltype(x)) のように初期化する.\n# Julia言語では引数の型を指定しない函数を引数の型が変化しても適切に動くように書く習慣がある.\n# s = zero(eltype(x)) のような書き方はそのような習慣の一部分になっている.\n# 複数の引数型に対応している函数を書けることはJulia言語を使う大きな利点の一つである.\n# 単に速さだけを気にするなら s=0.0 という書き方で十分だが, \n# Julia言語の長所を活かし切るためには s=zero(eltype(x)) のような書き方も覚える必要がある.\n#\n# 下の方に s=0.0 と s=zero(eltype(x)) の違いで有理数の和がどのように変わるかの例がある.\n\nsrand(2018)\nx = randn(10^8)\n@code_warntype  mysum_int(x)",
      "execution_count": 8,
      "outputs": [
        {
          "output_type": "stream",
          "text": "Variables:\n  #self# <optimized out>\n  x::Array{Float64,1}\n  i::Int64\n  #temp#@_4::Int64\n  s\u001b[1m\u001b[91m::Union{Float64, Int64}\u001b[39m\u001b[22m\n  #temp#@_6::Core.MethodInstance\n  #temp#@_7::Float64\n\nBody:\n  begin \n      s\u001b[1m\u001b[91m::Union{Float64, Int64}\u001b[39m\u001b[22m = 0 # line 3:\n      $(Expr(:inbounds, false))\n      # meta: location abstractarray.jl endof 134\n      # meta: location abstractarray.jl linearindices 99\n      # meta: location abstractarray.jl indices1 71\n      # meta: location abstractarray.jl indices 64\n      SSAValue(4) = (Base.arraysize)(x::Array{Float64,1}, 1)::Int64\n      # meta: pop location\n      # meta: pop location\n      # meta: pop location\n      # meta: pop location\n      $(Expr(:inbounds, :pop))\n      SSAValue(5) = (Base.select_value)((Base.slt_int)(SSAValue(4), 0)::Bool, 0, SSAValue(4))::Int64\n      SSAValue(6) = (Base.select_value)((Base.sle_int)(1, SSAValue(5))::Bool, SSAValue(5), (Base.sub_int)(1, 1)::Int64)::Int64\n      #temp#@_4::Int64 = 1\n      17: \n      unless (Base.not_int)((#temp#@_4::Int64 === (Base.add_int)(SSAValue(6), 1)::Int64)::Bool)::Bool goto 42\n      SSAValue(7) = #temp#@_4::Int64\n      SSAValue(8) = (Base.add_int)(#temp#@_4::Int64, 1)::Int64\n      i::Int64 = SSAValue(7)\n      #temp#@_4::Int64 = SSAValue(8) # line 4:\n      unless (s\u001b[1m\u001b[91m::Union{Float64, Int64}\u001b[39m\u001b[22m isa Int64)::Bool goto 27\n      #temp#@_6::Core.MethodInstance = MethodInstance for +(::Int64, ::Float64)\n      goto 36\n      27: \n      unless (s\u001b[1m\u001b[91m::Union{Float64, Int64}\u001b[39m\u001b[22m isa Float64)::Bool goto 31\n      #temp#@_6::Core.MethodInstance = MethodInstance for +(::Float64, ::Float64)\n      goto 36\n      31: \n      goto 33\n      33: \n      #temp#@_7::Float64 = (s\u001b[1m\u001b[91m::Union{Float64, Int64}\u001b[39m\u001b[22m + (Base.arrayref)(x::Array{Float64,1}, i::Int64)::Float64)::Float64\n      goto 38\n      36: \n      #temp#@_7::Float64 = $(Expr(:invoke, :(#temp#@_6), :(Main.+), :(s), :((Base.arrayref)(x, i)::Float64)))\n      38: \n      s\u001b[1m\u001b[91m::Union{Float64, Int64}\u001b[39m\u001b[22m = #temp#@_7::Float64\n      40: \n      goto 17\n      42:  # line 6:\n      return s\u001b[1m\u001b[91m::Union{Float64, Int64}\u001b[39m\u001b[22m\n  end\u001b[1m\u001b[91m::Union{Float64, Int64}\u001b[39m\u001b[22m\n",
          "name": "stdout"
        }
      ]
    },
    {
      "metadata": {
        "trusted": true
      },
      "cell_type": "code",
      "source": "# mysum_float (s=0.0と初期化)だと有理数の和が浮動小数点数 2.9289682539682538 になってしまう\n\nx_rational = [1//n for n in 1:10]\n@show x_rational\n@show mysum_float(x_rational);",
      "execution_count": 9,
      "outputs": [
        {
          "output_type": "stream",
          "text": "x_rational = Rational{Int64}[1//1, 1//2, 1//3, 1//4, 1//5, 1//6, 1//7, 1//8, 1//9, 1//10]\nmysum_float(x_rational) = 2.9289682539682538\n",
          "name": "stdout"
        }
      ]
    },
    {
      "metadata": {
        "trusted": true
      },
      "cell_type": "code",
      "source": "# mysum_eltype (s=zero(eltype(x))と初期化)なら有理数の和が有理数になってくれる\n\nx_rational = [1//n for n in 1:10]\n@show x_rational\n@show mysum_eltype(x_rational);",
      "execution_count": 10,
      "outputs": [
        {
          "output_type": "stream",
          "text": "x_rational = Rational{Int64}[1//1, 1//2, 1//3, 1//4, 1//5, 1//6, 1//7, 1//8, 1//9, 1//10]\nmysum_eltype(x_rational) = 7381//2520\n",
          "name": "stdout"
        }
      ]
    },
    {
      "metadata": {
        "trusted": true
      },
      "cell_type": "code",
      "source": "# 組み込み函数の sum でも有理数の和は当然有理数になってくれる\n\nx_rational = [1//n for n in 1:10]\n@show x_rational\n@show sum(x_rational);",
      "execution_count": 11,
      "outputs": [
        {
          "output_type": "stream",
          "text": "x_rational = Rational{Int64}[1//1, 1//2, 1//3, 1//4, 1//5, 1//6, 1//7, 1//8, 1//9, 1//10]\nsum(x_rational) = 7381//2520\n",
          "name": "stdout"
        }
      ]
    },
    {
      "metadata": {
        "trusted": true
      },
      "cell_type": "code",
      "source": "",
      "execution_count": null,
      "outputs": []
    }
  ],
  "metadata": {
    "kernelspec": {
      "name": "julia-0.6",
      "display_name": "Julia 0.6.2",
      "language": "julia"
    },
    "toc": {
      "threshold": 4,
      "number_sections": true,
      "toc_cell": false,
      "toc_window_display": false,
      "toc_section_display": "block",
      "sideBar": true,
      "navigate_menu": true,
      "moveMenuLeft": true,
      "widenNotebook": false,
      "colors": {
        "hover_highlight": "#DAA520",
        "selected_highlight": "#FFD700",
        "running_highlight": "#FF0000",
        "wrapper_background": "#FFFFFF",
        "sidebar_border": "#EEEEEE",
        "navigate_text": "#333333",
        "navigate_num": "#000000"
      },
      "nav_menu": {
        "width": "252px",
        "height": "48px"
      }
    },
    "language_info": {
      "file_extension": ".jl",
      "name": "julia",
      "mimetype": "application/julia",
      "version": "0.6.2"
    },
    "gist": {
      "id": "",
      "data": {
        "description": "変数を適切な型で初期化することについて",
        "public": true
      }
    }
  },
  "nbformat": 4,
  "nbformat_minor": 2
}