Ising model のPython実装

pyising.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Ising Python version "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "簡単な解説　\n",
    "Ising model - http://enakai00.hatenablog.com/entry/20150106/1420538321"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Pure Python"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Reference: http://nbviewer.jupyter.org/gist/genkuroki/4fa46c68c56ee0f3b1a6fc8ec628b9d7\n",
    "# Author MathSorcerer\n",
    "\n",
    "from math import log, sqrt, exp\n",
    "from random import choice, random\n",
    "from copy import deepcopy\n",
    "import time\n",
    "from matplotlib import pyplot as plt\n",
    "%matplotlib inline\n",
    "\n",
    "def ising2d_sum_of_adjacent_spins(s, m, n, i, j):\n",
    "    i_bottom = i+1 if i+1 < m else 0\n",
    "    i_top = i-1 if i-1 >= 0 else m-1\n",
    "    j_right = j+1 if j+1 < n else 0\n",
    "    j_left = j-1 if j-1 >= 0 else n-1\n",
    "    return s[i_bottom][j]+s[i_top][j]+s[i][j_right]+s[i][j_left]\n",
    "\n",
    "def ising2d_sweep(s, beta, niters):\n",
    "    m, n = len(s),len(s[0])\n",
    "    prob = [exp(-2*beta*k) for k in [-4, -3, -2, -1, 0, 1, 2, 3, 4]]\n",
    "    for _ in range(int(niters/(m*n))):\n",
    "        for i in range(m):\n",
    "            for j in range(n):\n",
    "                s1 = s[i][j]\n",
    "                k = s1*ising2d_sum_of_adjacent_spins(s, m, n, i, j)\n",
    "                s[i][j] = -s1 if random() < prob[k+4] else s1\n",
    "    return s\n",
    "\n",
    "\n",
    "def plot_result(s_begin,s_end):\n",
    "    fig = plt.figure()\n",
    "    ax1 = fig.add_subplot(121)\n",
    "    ax2 = fig.add_subplot(122)\n",
    "    ax1.imshow(s_begin, cmap='gray')\n",
    "    ax2.imshow(s_end, cmap='gray')\n",
    "    plt.show()\n",
    "\n",
    "n = 100\n",
    "beta_critical = log(1+sqrt(2))/2\n",
    "\n",
    "\n",
    "def main():\n",
    "    rand_ising2d = [[choice([-1, 1]) for j in range(n)] for i in range(n)]\n",
    "    s_begin = deepcopy(rand_ising2d)\n",
    "    begin = time.time()\n",
    "    s_end = ising2d_sweep(rand_ising2d, beta_critical, 1e6)\n",
    "    end = time.time()\n",
    "    print(\"Elapsed=\", end-begin)\n",
    "\n",
    "    #plot_result(s_begin,s_end)\n",
    "\n",
    "if __name__ == '__main__':\n",
    "    main()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Python with Numba(Macで動いた）"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Reference: http://nbviewer.jupyter.org/gist/genkuroki/4fa46c68c56ee0f3b1a6fc8ec628b9d7\n",
    "# Author MathSorcerer\n",
    "\n",
    "from math import log, sqrt, exp\n",
    "from random import choice, random\n",
    "from copy import deepcopy\n",
    "import time\n",
    "from array import array\n",
    "from matplotlib import pyplot as plt\n",
    "from numba import jit \n",
    "import numpy as np \n",
    "\n",
    "@jit('i8(i8[:,:],i8,i8,i8,i8)')\n",
    "def ising2d_sum_of_adjacent_spins(s, m, n, i, j):\n",
    "    i_bottom = i+1 if i+1 < m else 0\n",
    "    i_top = i-1 if i-1 >= 0 else m-1\n",
    "    j_right = j+1 if j+1 < n else 0\n",
    "    j_left = j-1 if j-1 >= 0 else n-1\n",
    "    return s[i_bottom][j]+s[i_top][j]+s[i][j_right]+s[i][j_left]\n",
    "\n",
    "@jit('i8[:,:](i8[:,:],f8,i8)')\n",
    "def ising2d_sweep(s, beta, niters):\n",
    "    m, n = s.shape\n",
    "    prob = [exp(-2*beta*k) for k in [-4, -3, -2, -1, 0, 1, 2, 3, 4]]\n",
    "    for _ in range(int(niters/(m*n))):\n",
    "        for i in range(m):\n",
    "            for j in range(n):\n",
    "                s1 = s[i][j]\n",
    "                k = s1*ising2d_sum_of_adjacent_spins(s, m, n, i, j)\n",
    "                s[i][j] = -s1 if random() < prob[k+4] else s1\n",
    "    return s\n",
    "\n",
    "\n",
    "\n",
    "def plot_result(s_begin,s_end):\n",
    "    fig = plt.figure()\n",
    "    ax1 = fig.add_subplot(121)\n",
    "    ax2 = fig.add_subplot(122)\n",
    "    ax1.imshow(s_begin, cmap='gray')\n",
    "    ax2.imshow(s_end, cmap='gray')\n",
    "    plt.show()\n",
    "\n",
    "n = 100\n",
    "beta_critical = log(1+sqrt(2))/2\n",
    "\n",
    "\n",
    "def main():\n",
    "    rand_ising2d = np.array([[choice([-1, 1]) for j in range(n)] for i in range(n)])\n",
    "    s_begin = deepcopy(rand_ising2d)\n",
    "    begin = time.time()\n",
    "    s_end = ising2d_sweep(rand_ising2d, beta_critical, 1e9)\n",
    "    end = time.time()\n",
    "    print(\"Elapsed=\", end-begin)\n",
    "\n",
    "    #plot_result(s_begin,s_end)\n",
    "\n",
    "\n",
    "\n",
    "if __name__ == '__main__':\n",
    "    main()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Python with Numba(Windowsで動いた）\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Reference: http://nbviewer.jupyter.org/gist/genkuroki/4fa46c68c56ee0f3b1a6fc8ec628b9d7\n",
    "# Author MathSorcerer\n",
    "\n",
    "from math import log, sqrt, exp\n",
    "from random import choice, random\n",
    "from copy import deepcopy\n",
    "import time\n",
    "from array import array\n",
    "from matplotlib import pyplot as plt\n",
    "from numba import jit, prange\n",
    "import numpy as np\n",
    "\n",
    "\n",
    "@jit('i4(i4[:,:],i4,i4,i4,i4)')\n",
    "def ising2d_sum_of_adjacent_spins(s, m, n, i, j):\n",
    "    i_bottom = i+1 if i+1 < m else 0\n",
    "    i_top = i-1 if i-1 >= 0 else m-1\n",
    "    j_right = j+1 if j+1 < n else 0\n",
    "    j_left = j-1 if j-1 >= 0 else n-1\n",
    "    return s[i_bottom, j]+s[i_top, j]+s[i, j_right]+s[i, j_left]\n",
    "\n",
    "\n",
    "@jit('i4[:,:](i4[:,:],f4,i4)', parallel=True, nopython=True)\n",
    "def ising2d_sweep(s, beta, niters):\n",
    "    m, n = s.shape\n",
    "    prob = np.exp(-2*beta*np.array([-4, -3, -2, -1, 0, 1, 2, 3, 4]))\n",
    "    for _ in range(int(niters/(m*n))):\n",
    "        for i in range(m):\n",
    "            for j in range(n):\n",
    "                s1 = s[i, j]\n",
    "                k = s1*ising2d_sum_of_adjacent_spins(s, m, n, i, j)\n",
    "                s[i, j] = -s1 if np.random.random() < prob[k+4] else s1\n",
    "                # s[i, j] = -s1 if random() < prob[k+4] else s1 <- very slow\n",
    "    return s\n",
    "\n",
    "\n",
    "def plot_result(s_begin, s_end):\n",
    "    fig = plt.figure()\n",
    "    ax1 = fig.add_subplot(121)\n",
    "    ax2 = fig.add_subplot(122)\n",
    "    ax1.imshow(s_begin, cmap='gray')\n",
    "    ax2.imshow(s_end, cmap='gray')\n",
    "    plt.show()\n",
    "\n",
    "n = 100\n",
    "beta_critical = log(1+sqrt(2))/2\n",
    "\n",
    "\n",
    "def main():\n",
    "    rand_ising2d = np.array([[choice([-1, 1])\n",
    "                              for j in range(n)] for i in range(n)])\n",
    "    s_begin = deepcopy(rand_ising2d)\n",
    "    begin = time.time()\n",
    "    s_end = ising2d_sweep(rand_ising2d, beta_critical, 1e9)\n",
    "    end = time.time()\n",
    "    print(\"Elapsed=\", end-begin)\n",
    "\n",
    "    # plot_result(s_begin,s_end)\n",
    "\n",
    "\n",
    "if __name__ == '__main__':\n",
    "    main()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Python with Cython"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%load_ext Cython"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%%cython --annotate\n",
    "\n",
    "# Reference: http://nbviewer.jupyter.org/gist/genkuroki/4fa46c68c56ee0f3b1a6fc8ec628b9d7\n",
    "# Author MathSorcerer\n",
    "\n",
    "from math import exp\n",
    "import array\n",
    "from cpython cimport array as carray\n",
    "from random import choice, random\n",
    "from libc.stdlib cimport rand, RAND_MAX\n",
    "import time\n",
    "cimport cython\n",
    "import numpy as np\n",
    "\n",
    "@cython.boundscheck(False)\n",
    "cdef int ising2d_sum_of_adjacent_spins(int[:,:] s, int m, int n,int  i, int j):\n",
    "    cdef:\n",
    "        int i_bottom = i+1 if i+1 < m else 0\n",
    "        int i_top = i-1 if i-1 >= 0 else m-1\n",
    "        int j_right = j+1 if j+1 < n else 0\n",
    "        int j_left = j-1 if j-1 >= 0 else n-1\n",
    "    return s[i_bottom,j]+s[i_top,j]+s[i,j_right]+s[i,j_left]\n",
    "\n",
    "@cython.boundscheck(False)\n",
    "@cython.wraparound(False)\n",
    "cpdef int[:,:] ising2d_sweep(int[:,:] s, double beta, int niters):\n",
    "    cdef int m, n, s1, loop,num_iteration,k\n",
    "    m,n= s.shape[0],s.shape[1]\n",
    "    num_iteration=int(niters/(m*n))\n",
    "    cdef carray.array prob = array.array('d',[exp(-2*beta*k) for k in [-4, -3, -2, -1, 0, 1, 2, 3, 4]])\n",
    "    for loop in range(num_iteration):\n",
    "        for i in range(m):\n",
    "            for j in range(n):\n",
    "                s1 = s[i,j]\n",
    "                k = s1*ising2d_sum_of_adjacent_spins(s, m, n, i, j)\n",
    "                s[i,j] = -s1 if rand()/RAND_MAX <  prob[k+4] else s1\n",
    "    return s"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "\n",
    "from copy import deepcopy\n",
    "from math import log,sqrt \n",
    "import numpy as np\n",
    "n=100\n",
    "beta_critical = log(1+sqrt(2))/2\n",
    "\n",
    "def main():\n",
    "    rand_ising2d = np.array([[choice([-1, 1]) for j in range(n)] for i in range(n)]).astype(np.int32)\n",
    "    s_begin = deepcopy(rand_ising2d)\n",
    "    begin = time.time()\n",
    "    s_end = ising2d_sweep(rand_ising2d, beta_critical, 1e8)\n",
    "    end = time.time()\n",
    "    print(\"Elapsed=\", end-begin)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "main()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}