{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Juliaで学ぶ古典モンテカルロシミュレーション\n",
    "## おまけ：書き方による速度の違い\n",
    "せっかくなので、一番コンパクトな書き方をしてみよう。そして、少しづつ書き方を変えてみる。速度差は出るだろうか。  \n",
    "まず、そのまま書き下してみる。\n",
    "計算環境はJulia Boxを使用  \n",
    "https://www.juliabox.com  \n",
    "し、計算機環境は\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Julia Version 0.6.0\n",
      "Commit 9036443 (2017-06-19 13:05 UTC)\n",
      "Platform Info:\n",
      "  OS: Linux (x86_64-pc-linux-gnu)\n",
      "  CPU: Intel(R) Xeon(R) CPU E5-2673 v3 @ 2.40GHz\n",
      "  WORD_SIZE: 64\n",
      "  BLAS: libopenblas (USE64BITINT DYNAMIC_ARCH NO_AFFINITY Haswell)\n",
      "  LAPACK: libopenblas64_\n",
      "  LIBM: libopenlibm\n",
      "  LLVM: libLLVM-3.9.1 (ORCJIT, haswell)\n"
     ]
    }
   ],
   "source": [
    "versioninfo()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "となっている。  \n",
    "二次元イジング模型のメトロポリス法によるモンテカルロシミュレーションのコードは、"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "モンテカルロ法! (generic function with 1 method)"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function スイープ!(配置,Lx,Ly,w,mz)\n",
    "    for ix in 1:Lx\n",
    "        for iy in 1:Ly  \n",
    "            σi = 配置[ix,iy]\n",
    "            ΔE =  2σi*(配置[ifelse(ix==Lx,1,ix+1),iy]+配置[ifelse(ix==1,Lx,ix-1),iy]+配置[ix,ifelse(iy==Ly,1,iy+1)]+配置[ix,ifelse(iy==1,Ly,iy-1)])\n",
    "            配置[ix,iy] = ifelse(rand() <= w[ΔE+9],-σi,σi)\n",
    "            mz += ifelse(配置[ix,iy] == -σi,-2*σi,0)\n",
    "        end\n",
    "    end\n",
    "    return 配置,mz\n",
    "end\n",
    "\n",
    "function モンテカルロ法!(配置,Lx,Ly,T,最大ステップ,熱化ステップ)\n",
    "    w = exp.(-[i for i in -8:8]/T) #ボルツマン重み e^[-ΔE/T]\n",
    "    mz = sum(配置)\n",
    "    Mz = 0\n",
    "    for i=1:熱化ステップ\n",
    "        配置,mz=スイープ!(配置,Lx,Ly,w,mz)        \n",
    "    end    \n",
    "    for i=1:最大ステップ\n",
    "        配置,mz=スイープ!(配置,Lx,Ly,w,mz)\n",
    "        Mz += abs(mz)\n",
    "    end\n",
    "    return Mz/(Lx*Ly*最大ステップ)\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "となり、このくらい短く書くことができる。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "106.624404 seconds (1.11 M allocations: 36.195 MiB, 0.01% gc time)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "0.9920440608"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Lx = 100\n",
    "Ly = 100\n",
    "T = 1.0\n",
    "最大ステップ = 1000000\n",
    "熱化ステップ = 200\n",
    "srand(123)\n",
    "配置 = rand([-1,1],Lx,Ly)\n",
    "T = 1.0\n",
    "@time Mz = モンテカルロ法!(配置,Lx,Ly,T,最大ステップ,熱化ステップ)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "次に、乱数を少し先に確保してみる。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "スイープ! (generic function with 1 method)"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function スイープ!(配置,Lx,Ly,w,mz)\n",
    "    r = rand(Lx,Ly) #乱数確保\n",
    "    for ix in 1:Lx\n",
    "        for iy in 1:Ly  \n",
    "            σi = 配置[ix,iy]\n",
    "            ΔE =  2σi*(配置[ifelse(ix==Lx,1,ix+1),iy]+配置[ifelse(ix==1,Lx,ix-1),iy]+配置[ix,ifelse(iy==Ly,1,iy+1)]+配置[ix,ifelse(iy==1,Ly,iy-1)])\n",
    "            配置[ix,iy] = ifelse(r[ix,iy] <= w[ΔE+9],-σi,σi)\n",
    "            mz += ifelse(配置[ix,iy] == -σi,-2*σi,0)\n",
    "        end\n",
    "    end\n",
    "    return 配置,mz\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "114.048818 seconds (3.03 M allocations: 74.626 GiB, 2.45% gc time)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "0.988445197"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Lx = 100\n",
    "Ly = 100\n",
    "T = 1.0\n",
    "最大ステップ = 1000000\n",
    "熱化ステップ = 200\n",
    "srand(123)\n",
    "配置 = rand([-1,1],Lx,Ly)\n",
    "T = 1.0\n",
    "@time Mz = モンテカルロ法!(配置,Lx,Ly,T,最大ステップ,熱化ステップ)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "少し遅くなった？　念のためもう一回実行する。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "115.313045 seconds (3.00 M allocations: 74.625 GiB, 2.53% gc time)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "0.988445197"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Lx = 100\n",
    "Ly = 100\n",
    "T = 1.0\n",
    "最大ステップ = 1000000\n",
    "熱化ステップ = 200\n",
    "srand(123)\n",
    "配置 = rand([-1,1],Lx,Ly)\n",
    "T = 1.0\n",
    "@time Mz = モンテカルロ法!(配置,Lx,Ly,T,最大ステップ,熱化ステップ)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "やはり遅くなっている。配列を確保したのが原因だろうか？  \n",
    "\n",
    "乱数に関しては元に戻して、次に、ifelse関数を三項演算子「a ? b : c 」に変えてみる"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "スイープ! (generic function with 1 method)"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function スイープ!(配置,Lx,Ly,w,mz)\n",
    "    for ix in 1:Lx\n",
    "        for iy in 1:Ly  \n",
    "            σi = 配置[ix,iy]\n",
    "            ΔE =  2σi*(配置[ix==Lx?1:ix+1,iy]+配置[ix==1?Lx:ix-1,iy]+配置[ix,iy==Ly?1:iy+1]+配置[ix,iy==1?Ly:iy-1])\n",
    "            配置[ix,iy] = ifelse(rand() <= w[ΔE+9],-σi,σi)\n",
    "            mz += ifelse(配置[ix,iy] == -σi,-2*σi,0)\n",
    "        end\n",
    "    end\n",
    "    return 配置,mz\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "118.716406 seconds (1.01 M allocations: 31.153 MiB, 0.01% gc time)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "0.9920440608"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Lx = 100\n",
    "Ly = 100\n",
    "T = 1.0\n",
    "最大ステップ = 1000000\n",
    "熱化ステップ = 200\n",
    "srand(123)\n",
    "配置 = rand([-1,1],Lx,Ly)\n",
    "T = 1.0\n",
    "@time Mz = モンテカルロ法!(配置,Lx,Ly,T,最大ステップ,熱化ステップ)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "遅い気がする。念のためもう一度。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "121.250223 seconds (1.00 M allocations: 30.525 MiB, 0.01% gc time)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "0.9920440608"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Lx = 100\n",
    "Ly = 100\n",
    "T = 1.0\n",
    "最大ステップ = 1000000\n",
    "熱化ステップ = 200\n",
    "srand(123)\n",
    "配置 = rand([-1,1],Lx,Ly)\n",
    "T = 1.0\n",
    "@time Mz = モンテカルロ法!(配置,Lx,Ly,T,最大ステップ,熱化ステップ)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "うーん。これは遅い。つまり、現時点では"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "スイープ! (generic function with 1 method)"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function スイープ!(配置,Lx,Ly,w,mz)\n",
    "    for ix in 1:Lx\n",
    "        for iy in 1:Ly  \n",
    "            σi = 配置[ix,iy]\n",
    "            ΔE =  2σi*(配置[ifelse(ix==Lx,1,ix+1),iy]+配置[ifelse(ix==1,Lx,ix-1),iy]+配置[ix,ifelse(iy==Ly,1,iy+1)]+配置[ix,ifelse(iy==1,Ly,iy-1)])\n",
    "            配置[ix,iy] = ifelse(rand() <= w[ΔE+9],-σi,σi)\n",
    "            mz += ifelse(配置[ix,iy] == -σi,-2*σi,0)\n",
    "        end\n",
    "    end\n",
    "    return 配置,mz\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "104.852433 seconds (1.01 M allocations: 31.152 MiB, 0.01% gc time)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "0.9920440608"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Lx = 100\n",
    "Ly = 100\n",
    "T = 1.0\n",
    "最大ステップ = 1000000\n",
    "熱化ステップ = 200\n",
    "srand(123)\n",
    "配置 = rand([-1,1],Lx,Ly)\n",
    "T = 1.0\n",
    "@time Mz = モンテカルロ法!(配置,Lx,Ly,T,最大ステップ,熱化ステップ)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "が一番速い。  \n",
    "次に、ifelse関数を余りを計算するmod1(ix-1,Lx)で置き換えてみよう。この関数は"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "100\n",
      "1\n"
     ]
    }
   ],
   "source": [
    "ix = 1\n",
    "println(mod1(ix-1,Lx))\n",
    "ix = Lx\n",
    "println(mod1(ix+1,Lx))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "と周期境界条件を表現できる。これを使ってみる。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "スイープ! (generic function with 1 method)"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function スイープ!(配置,Lx,Ly,w,mz)\n",
    "    for ix in 1:Lx\n",
    "        for iy in 1:Ly  \n",
    "            σi = 配置[ix,iy]\n",
    "            ΔE =  2σi*(配置[mod1(ix+1,Lx),iy]+配置[mod1(ix-1,Lx),iy]+配置[ix,mod1(iy+1,Ly)]+配置[ix,mod1(iy-1,Ly)])\n",
    "            配置[ix,iy] = ifelse(rand() <= w[ΔE+9],-σi,σi)\n",
    "            mz += ifelse(配置[ix,iy] == -σi,-2*σi,0)\n",
    "        end\n",
    "    end\n",
    "    return 配置,mz\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "628.313394 seconds (1.01 M allocations: 31.204 MiB, 0.00% gc time)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "0.9920440608"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Lx = 100\n",
    "Ly = 100\n",
    "T = 1.0\n",
    "最大ステップ = 1000000\n",
    "熱化ステップ = 200\n",
    "srand(123)\n",
    "配置 = rand([-1,1],Lx,Ly)\n",
    "T = 1.0\n",
    "@time Mz = モンテカルロ法!(配置,Lx,Ly,T,最大ステップ,熱化ステップ)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "とても遅くなった。ifelse関数が一番速い。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "次に、コードが少し汚くなるが、ifelse関数も除去してみよう。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "スイープ! (generic function with 1 method)"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function スイープ!(配置,Lx,Ly,w,mz)\n",
    "    ix = 1\n",
    "    iy = 1\n",
    "    σi = 配置[ix,iy]\n",
    "    ΔE =  2σi*(配置[ix+1,iy]+配置[Lx,iy]+配置[ix,iy+1]+配置[ix,Ly])\n",
    "    配置[ix,iy] = ifelse(rand() <= w[ΔE+9],-σi,σi)\n",
    "    mz += ifelse(配置[ix,iy] == -σi,-2*σi,0)  \n",
    "    \n",
    "    for iy in 2:Ly-1\n",
    "        σi = 配置[ix,iy]\n",
    "        ΔE =  2σi*(配置[ix+1,iy]+配置[Lx,iy]+配置[ix,iy+1]+配置[ix,iy-1])\n",
    "        配置[ix,iy] = ifelse(rand() <= w[ΔE+9],-σi,σi)\n",
    "        mz += ifelse(配置[ix,iy] == -σi,-2*σi,0)        \n",
    "    end\n",
    "    iy = Ly\n",
    "    σi = 配置[ix,iy]\n",
    "    ΔE =  2σi*(配置[ix+1,iy]+配置[Lx,iy]+配置[ix,1]+配置[ix,iy-1])\n",
    "    配置[ix,iy] = ifelse(rand() <= w[ΔE+9],-σi,σi)\n",
    "    mz += ifelse(配置[ix,iy] == -σi,-2*σi,0)      \n",
    "        \n",
    "    for ix in 2:Lx-1\n",
    "        iy = 1\n",
    "        σi = 配置[ix,iy]\n",
    "        ΔE =  2σi*(配置[ix+1,iy]+配置[ix-1,iy]+配置[ix,iy+1]+配置[ix,Ly])\n",
    "        配置[ix,iy] = ifelse(rand() <= w[ΔE+9],-σi,σi)\n",
    "        mz += ifelse(配置[ix,iy] == -σi,-2*σi,0)  \n",
    "        \n",
    "        for iy in 2:Ly-1  \n",
    "            σi = 配置[ix,iy]\n",
    "            ΔE =  2σi*(配置[ix+1,iy]+配置[iy-1,iy]+配置[ix,iy+1]+配置[ix,iy-1])\n",
    "            配置[ix,iy] = ifelse(rand() <= w[ΔE+9],-σi,σi)\n",
    "            mz += ifelse(配置[ix,iy] == -σi,-2*σi,0)\n",
    "        end\n",
    "        \n",
    "        iy = Ly\n",
    "        σi = 配置[ix,iy]\n",
    "        ΔE =  2σi*(配置[ix+1,iy]+配置[ix-1,iy]+配置[ix,1]+配置[ix,iy-1])\n",
    "        配置[ix,iy] = ifelse(rand() <= w[ΔE+9],-σi,σi)\n",
    "        mz += ifelse(配置[ix,iy] == -σi,-2*σi,0)             \n",
    "    end\n",
    "\n",
    "    ix = Lx\n",
    "    iy = 1\n",
    "    σi = 配置[ix,iy]\n",
    "    ΔE =  2σi*(配置[1,iy]+配置[ix-1,iy]+配置[ix,iy+1]+配置[ix,Ly])\n",
    "    配置[ix,iy] = ifelse(rand() <= w[ΔE+9],-σi,σi)\n",
    "    mz += ifelse(配置[ix,iy] == -σi,-2*σi,0)  \n",
    "    \n",
    "    for iy in 2:Ly-1\n",
    "        σi = 配置[ix,iy]\n",
    "        ΔE =  2σi*(配置[1,iy]+配置[ix-1,iy]+配置[ix,iy+1]+配置[ix,iy-1])\n",
    "        配置[ix,iy] = ifelse(rand() <= w[ΔE+9],-σi,σi)\n",
    "        mz += ifelse(配置[ix,iy] == -σi,-2*σi,0)        \n",
    "    end\n",
    "    iy = Ly\n",
    "    σi = 配置[ix,iy]\n",
    "    ΔE =  2σi*(配置[1,iy]+配置[ix-1,iy]+配置[ix,1]+配置[ix,iy-1])\n",
    "    配置[ix,iy] = ifelse(rand() <= w[ΔE+9],-σi,σi)\n",
    "    mz += ifelse(配置[ix,iy] == -σi,-2*σi,0)  \n",
    "    \n",
    "    \n",
    "    return 配置,mz\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "101.889886 seconds (1.04 M allocations: 32.397 MiB, 0.01% gc time)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "0.9990793904"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Lx = 100\n",
    "Ly = 100\n",
    "T = 1.0\n",
    "最大ステップ = 1000000\n",
    "熱化ステップ = 200\n",
    "srand(123)\n",
    "配置 = rand([-1,1],Lx,Ly)\n",
    "T = 1.0\n",
    "@time Mz = モンテカルロ法!(配置,Lx,Ly,T,最大ステップ,熱化ステップ)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "すごい汚いコードになったが速くはなった。どうせならさらにいってみよう。mzの部分はもう少し高速化できる。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "スイープ! (generic function with 1 method)"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function スイープ!(配置,Lx,Ly,w,mz)\n",
    "    ix = 1\n",
    "    iy = 1\n",
    "    σi = 配置[ix,iy]\n",
    "    ΔE =  2σi*(配置[ix+1,iy]+配置[Lx,iy]+配置[ix,iy+1]+配置[ix,Ly])\n",
    "    配置[ix,iy] = ifelse(rand() <= w[ΔE+9],-σi,σi)\n",
    "    mz += 配置[ix,iy]-σi\n",
    "    \n",
    "    for iy in 2:Ly-1\n",
    "        σi = 配置[ix,iy]\n",
    "        ΔE =  2σi*(配置[ix+1,iy]+配置[Lx,iy]+配置[ix,iy+1]+配置[ix,iy-1])\n",
    "        配置[ix,iy] = ifelse(rand() <= w[ΔE+9],-σi,σi)\n",
    "        mz += 配置[ix,iy]-σi       \n",
    "    end\n",
    "    iy = Ly\n",
    "    σi = 配置[ix,iy]\n",
    "    ΔE =  2σi*(配置[ix+1,iy]+配置[Lx,iy]+配置[ix,1]+配置[ix,iy-1])\n",
    "    配置[ix,iy] = ifelse(rand() <= w[ΔE+9],-σi,σi)\n",
    "    mz += 配置[ix,iy]-σi     \n",
    "        \n",
    "    for ix in 2:Lx-1\n",
    "        iy = 1\n",
    "        σi = 配置[ix,iy]\n",
    "        ΔE =  2σi*(配置[ix+1,iy]+配置[ix-1,iy]+配置[ix,iy+1]+配置[ix,Ly])\n",
    "        配置[ix,iy] = ifelse(rand() <= w[ΔE+9],-σi,σi)\n",
    "        mz += 配置[ix,iy]-σi\n",
    "        \n",
    "        for iy in 2:Ly-1  \n",
    "            σi = 配置[ix,iy]\n",
    "            ΔE =  2σi*(配置[ix+1,iy]+配置[iy-1,iy]+配置[ix,iy+1]+配置[ix,iy-1])\n",
    "            配置[ix,iy] = ifelse(rand() <= w[ΔE+9],-σi,σi)\n",
    "            mz += 配置[ix,iy]-σi\n",
    "        end\n",
    "        \n",
    "        iy = Ly\n",
    "        σi = 配置[ix,iy]\n",
    "        ΔE =  2σi*(配置[ix+1,iy]+配置[ix-1,iy]+配置[ix,1]+配置[ix,iy-1])\n",
    "        配置[ix,iy] = ifelse(rand() <= w[ΔE+9],-σi,σi)\n",
    "        mz += 配置[ix,iy]-σi          \n",
    "    end\n",
    "\n",
    "    ix = Lx\n",
    "    iy = 1\n",
    "    σi = 配置[ix,iy]\n",
    "    ΔE =  2σi*(配置[1,iy]+配置[ix-1,iy]+配置[ix,iy+1]+配置[ix,Ly])\n",
    "    配置[ix,iy] = ifelse(rand() <= w[ΔE+9],-σi,σi)\n",
    "    mz += 配置[ix,iy]-σi \n",
    "    \n",
    "    for iy in 2:Ly-1\n",
    "        σi = 配置[ix,iy]\n",
    "        ΔE =  2σi*(配置[1,iy]+配置[ix-1,iy]+配置[ix,iy+1]+配置[ix,iy-1])\n",
    "        配置[ix,iy] = ifelse(rand() <= w[ΔE+9],-σi,σi)\n",
    "        mz += 配置[ix,iy]-σi      \n",
    "    end\n",
    "    iy = Ly\n",
    "    σi = 配置[ix,iy]\n",
    "    ΔE =  2σi*(配置[1,iy]+配置[ix-1,iy]+配置[ix,1]+配置[ix,iy-1])\n",
    "    配置[ix,iy] = ifelse(rand() <= w[ΔE+9],-σi,σi)\n",
    "    mz += 配置[ix,iy]-σi\n",
    "    \n",
    "    \n",
    "    return 配置,mz\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "101.138176 seconds (1.04 M allocations: 32.277 MiB, 0.01% gc time)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "0.9990793904"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Lx = 100\n",
    "Ly = 100\n",
    "T = 1.0\n",
    "最大ステップ = 1000000\n",
    "熱化ステップ = 200\n",
    "srand(123)\n",
    "配置 = rand([-1,1],Lx,Ly)\n",
    "T = 1.0\n",
    "@time Mz = モンテカルロ法!(配置,Lx,Ly,T,最大ステップ,熱化ステップ)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "これは余り速くならなかった。しかし、あえてifelse関数を使う必要はないので、この変更は採用しよう。  \n",
    "結局、見やすさと速さを両立しているのは、"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "モンテカルロ法! (generic function with 1 method)"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function スイープ!(配置,Lx,Ly,w,mz)\n",
    "    for ix in 1:Lx\n",
    "        for iy in 1:Ly  \n",
    "            σi = 配置[ix,iy]\n",
    "            ΔE =  2σi*(配置[ifelse(ix==Lx,1,ix+1),iy]+配置[ifelse(ix==1,Lx,ix-1),iy]+配置[ix,ifelse(iy==Ly,1,iy+1)]+配置[ix,ifelse(iy==1,Ly,iy-1)])\n",
    "            配置[ix,iy] = ifelse(rand() <= w[ΔE+9],-σi,σi)\n",
    "            mz += 配置[ix,iy]-σi\n",
    "        end\n",
    "    end\n",
    "    return 配置,mz\n",
    "end\n",
    "\n",
    "\n",
    "function モンテカルロ法!(配置,Lx,Ly,T,最大ステップ,熱化ステップ)\n",
    "    w = exp.(-[i for i in -8:8]/T) #ボルツマン重み e^[-ΔE/T]\n",
    "    mz = sum(配置)\n",
    "    Mz = 0\n",
    "    for i=1:熱化ステップ\n",
    "        配置,mz=スイープ!(配置,Lx,Ly,w,mz)        \n",
    "    end    \n",
    "    for i=1:最大ステップ\n",
    "        配置,mz=スイープ!(配置,Lx,Ly,w,mz)\n",
    "        Mz += abs(mz)\n",
    "    end\n",
    "    return Mz/(Lx*Ly*最大ステップ)\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "101.150472 seconds (1.02 M allocations: 31.624 MiB, 0.01% gc time)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "0.9920440608"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Lx = 100\n",
    "Ly = 100\n",
    "T = 1.0\n",
    "最大ステップ = 1000000\n",
    "熱化ステップ = 200\n",
    "srand(123)\n",
    "配置 = rand([-1,1],Lx,Ly)\n",
    "T = 1.0\n",
    "@time Mz = モンテカルロ法!(配置,Lx,Ly,T,最大ステップ,熱化ステップ)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "この形が良さそうだ。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Plots.PyPlotBackend()"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "using Plots\n",
    "pyplot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Lx = 100\n",
    "Ly = 100\n",
    "最大ステップ = 100000\n",
    "熱化ステップ = 200\n",
    "srand(123)\n",
    "配置 = rand([-1,1],Lx,Ly)\n",
    "nt = 20\n",
    "time = [(nt+1-i)*0.2 for i in 1:nt]\n",
    "Mztdep = zeros(Float64,nt)\n",
    "i = 0\n",
    "@time for T in time\n",
    "    i += 1\n",
    "    Mztdep[i]  = モンテカルロ法!(配置,Lx,Ly,T,最大ステップ,熱化ステップ)\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot(time,Mztdep)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 0.6.0",
   "language": "julia",
   "name": "julia-0.6"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "0.6.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}