{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Cython\n", "\n", "* Example from [here](https://github.com/shekhar270779/Learn_Pandas/blob/master/Cython%20in%20Jupyter-notebook.ipynb) and [this tutorial](https://www.youtube.com/watch?v=CC0IkiNByV4)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Cython** is a superset of python programming language that gives access into c/c++ constructs. \n", "- Cython acts as a bridge between Python and C/C++." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Load cython notebook extension (and i have already installed cython )" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Requirement already satisfied: Cython in /srv/conda/envs/notebook/lib/python3.7/site-packages (0.29.24)\n", "Note: you may need to restart the kernel to use updated packages.\n" ] } ], "source": [ "pip install Cython" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "%load_ext cython" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "%%cython\n", "'''Inside this cell is Cython Code only\n", "'''\n", "\n", "def first_cython_function(int i):\n", " return i * 3.1415926;\n" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "31.415926" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# use cython function in a separate (python) cell\n", "first_cython_function(10)" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "\u001b[0;31mDocstring:\u001b[0m\n", "::\n", "\n", " %cython [-a] [-+] [-3] [-2] [-f] [-c COMPILE_ARGS]\n", " [--link-args LINK_ARGS] [-l LIB] [-n NAME] [-L dir] [-I INCLUDE]\n", " [-S SRC] [--pgo] [--verbose]\n", "\n", "Compile and import everything from a Cython code cell.\n", "\n", "The contents of the cell are written to a `.pyx` file in the\n", "directory `IPYTHONDIR/cython` using a filename with the hash of the\n", "code. This file is then cythonized and compiled. The resulting module\n", "is imported and all of its symbols are injected into the user's\n", "namespace. The usage is similar to that of `%%cython_pyximport` but\n", "you don't have to pass a module name::\n", "\n", " %%cython\n", " def f(x):\n", " return 2.0*x\n", "\n", "To compile OpenMP codes, pass the required `--compile-args`\n", "and `--link-args`. For example with gcc::\n", "\n", " %%cython --compile-args=-fopenmp --link-args=-fopenmp\n", " ...\n", "\n", "To enable profile guided optimisation, pass the ``--pgo`` option.\n", "Note that the cell itself needs to take care of establishing a suitable\n", "profile when executed. This can be done by implementing the functions to\n", "optimise, and then calling them directly in the same cell on some realistic\n", "training data like this::\n", "\n", " %%cython --pgo\n", " def critical_function(data):\n", " for item in data:\n", " ...\n", "\n", " # execute function several times to build profile\n", " from somewhere import some_typical_data\n", " for _ in range(100):\n", " critical_function(some_typical_data)\n", "\n", "In Python 3.5 and later, you can distinguish between the profile and\n", "non-profile runs as follows::\n", "\n", " if \"_pgo_\" in __name__:\n", " ... # execute critical code here\n", "\n", "optional arguments:\n", " -a, --annotate Produce a colorized HTML version of the source.\n", " -+, --cplus Output a C++ rather than C file.\n", " -3 Select Python 3 syntax.\n", " -2 Select Python 2 syntax.\n", " -f, --force Force the compilation of a new module, even if the\n", " source has been previously compiled.\n", " -c COMPILE_ARGS, --compile-args COMPILE_ARGS\n", " Extra flags to pass to compiler via the\n", " `extra_compile_args` Extension flag (can be specified\n", " multiple times).\n", " --link-args LINK_ARGS\n", " Extra flags to pass to linker via the\n", " `extra_link_args` Extension flag (can be specified\n", " multiple times).\n", " -l LIB, --lib LIB Add a library to link the extension against (can be\n", " specified multiple times).\n", " -n NAME, --name NAME Specify a name for the Cython module.\n", " -L dir Add a path to the list of library directories (can be\n", " specified multiple times).\n", " -I INCLUDE, --include INCLUDE\n", " Add a path to the list of include directories (can be\n", " specified multiple times).\n", " -S SRC, --src SRC Add a path to the list of src files (can be specified\n", " multiple times).\n", " --pgo Enable profile guided optimisation in the C compiler.\n", " Compiles the cell twice and executes it in between to\n", " generate a runtime profile.\n", " --verbose Print debug information like generated .c/.cpp file\n", " location and exact gcc/g++ command invoked.\n", "\u001b[0;31mFile:\u001b[0m /srv/conda/envs/notebook/lib/python3.7/site-packages/Cython/Build/IpythonMagic.py\n" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "%%cython?" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "def foo(a, b):\n", " return a + b" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " 2 0 LOAD_FAST 0 (a)\n", " 2 LOAD_FAST 1 (b)\n", " 4 BINARY_ADD\n", " 6 RETURN_VALUE\n" ] } ], "source": [ "from dis import dis \n", "# disassebles compiled python objects\n", "dis(foo)" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [], "source": [ "%%cython\n", "def cyfoo(a, b):\n", " return a + b" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "120 ns ± 3.57 ns per loop (mean ± std. dev. of 7 runs, 10000000 loops each)\n" ] } ], "source": [ "%timeit foo(1000000, 2000000)" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "110 ns ± 2.29 ns per loop (mean ± std. dev. of 7 runs, 10000000 loops each)\n" ] } ], "source": [ "%timeit cyfoo(1000000, 2000000)" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "139 ns ± 3.65 ns per loop (mean ± std. dev. of 7 runs, 10000000 loops each)\n" ] } ], "source": [ "%timeit foo('x', 'y')" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "138 ns ± 5.23 ns per loop (mean ± std. dev. of 7 runs, 10000000 loops each)\n" ] } ], "source": [ "%timeit cyfoo('x', 'y')" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [], "source": [ "def pyfac(n):\n", " if n <= 1:\n", " return 1\n", " return n * pyfac(n -1)" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "3.51 µs ± 151 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)\n" ] } ], "source": [ "%timeit pyfac(20.0)" ] }, { "cell_type": "code", "execution_count": 44, "metadata": {}, "outputs": [], "source": [ "%%cython\n", "\n", "def cyfac(n):\n", " if n <= 1:\n", " return 1\n", " return n * cyfac(n-1)\n", "\n", "def cyfac_double(double n):\n", " if n <= 1:\n", " return 1.0\n", " return n * cyfac_double(n-1)" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "1.46 µs ± 47.7 ns per loop (mean ± std. dev. of 7 runs, 1000000 loops each)\n" ] } ], "source": [ "%timeit cyfac(20.0)" ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "914 ns ± 71 ns per loop (mean ± std. dev. of 7 runs, 1000000 loops each)\n" ] } ], "source": [ "%timeit cyfac_double(20.0)" ] }, { "cell_type": "code", "execution_count": 43, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", "\n", "\n", "\n", " \n", " Cython: _cython_magic_88dfdb3a94d86c859cfd2127a1590a68.pyx\n", " \n", "\n", "\n", "

Generated by Cython 0.29.24

\n", "

\n", " Yellow lines hint at Python interaction.
\n", " Click on a line that starts with a \"+\" to see the C code that Cython generated for it.\n", "

\n", "
 1:
\n", "
+2: cpdef double cyfac_double_fast(double n):
\n", "
static PyObject *__pyx_pw_46_cython_magic_88dfdb3a94d86c859cfd2127a1590a68_1cyfac_double_fast(PyObject *__pyx_self, PyObject *__pyx_arg_n); /*proto*/\n",
       "static double __pyx_f_46_cython_magic_88dfdb3a94d86c859cfd2127a1590a68_cyfac_double_fast(double __pyx_v_n, CYTHON_UNUSED int __pyx_skip_dispatch) {\n",
       "  double __pyx_r;\n",
       "  __Pyx_RefNannyDeclarations\n",
       "  __Pyx_RefNannySetupContext(\"cyfac_double_fast\", 0);\n",
       "/* … */\n",
       "  /* function exit code */\n",
       "  __pyx_L0:;\n",
       "  __Pyx_RefNannyFinishContext();\n",
       "  return __pyx_r;\n",
       "}\n",
       "\n",
       "/* Python wrapper */\n",
       "static PyObject *__pyx_pw_46_cython_magic_88dfdb3a94d86c859cfd2127a1590a68_1cyfac_double_fast(PyObject *__pyx_self, PyObject *__pyx_arg_n); /*proto*/\n",
       "static PyObject *__pyx_pw_46_cython_magic_88dfdb3a94d86c859cfd2127a1590a68_1cyfac_double_fast(PyObject *__pyx_self, PyObject *__pyx_arg_n) {\n",
       "  double __pyx_v_n;\n",
       "  PyObject *__pyx_r = 0;\n",
       "  __Pyx_RefNannyDeclarations\n",
       "  __Pyx_RefNannySetupContext(\"cyfac_double_fast (wrapper)\", 0);\n",
       "  assert(__pyx_arg_n); {\n",
       "    __pyx_v_n = __pyx_PyFloat_AsDouble(__pyx_arg_n); if (unlikely((__pyx_v_n == (double)-1) && PyErr_Occurred())) __PYX_ERR(0, 2, __pyx_L3_error)\n",
       "  }\n",
       "  goto __pyx_L4_argument_unpacking_done;\n",
       "  __pyx_L3_error:;\n",
       "  __Pyx_AddTraceback(\"_cython_magic_88dfdb3a94d86c859cfd2127a1590a68.cyfac_double_fast\", __pyx_clineno, __pyx_lineno, __pyx_filename);\n",
       "  __Pyx_RefNannyFinishContext();\n",
       "  return NULL;\n",
       "  __pyx_L4_argument_unpacking_done:;\n",
       "  __pyx_r = __pyx_pf_46_cython_magic_88dfdb3a94d86c859cfd2127a1590a68_cyfac_double_fast(__pyx_self, ((double)__pyx_v_n));\n",
       "  int __pyx_lineno = 0;\n",
       "  const char *__pyx_filename = NULL;\n",
       "  int __pyx_clineno = 0;\n",
       "\n",
       "  /* function exit code */\n",
       "  __Pyx_RefNannyFinishContext();\n",
       "  return __pyx_r;\n",
       "}\n",
       "\n",
       "static PyObject *__pyx_pf_46_cython_magic_88dfdb3a94d86c859cfd2127a1590a68_cyfac_double_fast(CYTHON_UNUSED PyObject *__pyx_self, double __pyx_v_n) {\n",
       "  PyObject *__pyx_r = NULL;\n",
       "  __Pyx_RefNannyDeclarations\n",
       "  __Pyx_RefNannySetupContext(\"cyfac_double_fast\", 0);\n",
       "  __Pyx_XDECREF(__pyx_r);\n",
       "  __pyx_t_1 = PyFloat_FromDouble(__pyx_f_46_cython_magic_88dfdb3a94d86c859cfd2127a1590a68_cyfac_double_fast(__pyx_v_n, 0)); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 2, __pyx_L1_error)\n",
       "  __Pyx_GOTREF(__pyx_t_1);\n",
       "  __pyx_r = __pyx_t_1;\n",
       "  __pyx_t_1 = 0;\n",
       "  goto __pyx_L0;\n",
       "\n",
       "  /* function exit code */\n",
       "  __pyx_L1_error:;\n",
       "  __Pyx_XDECREF(__pyx_t_1);\n",
       "  __Pyx_AddTraceback(\"_cython_magic_88dfdb3a94d86c859cfd2127a1590a68.cyfac_double_fast\", __pyx_clineno, __pyx_lineno, __pyx_filename);\n",
       "  __pyx_r = NULL;\n",
       "  __pyx_L0:;\n",
       "  __Pyx_XGIVEREF(__pyx_r);\n",
       "  __Pyx_RefNannyFinishContext();\n",
       "  return __pyx_r;\n",
       "}\n",
       "
+3:     if n <= 1:
\n", "
  __pyx_t_1 = ((__pyx_v_n <= 1.0) != 0);\n",
       "  if (__pyx_t_1) {\n",
       "/* … */\n",
       "  }\n",
       "
+4:         return 1.0
\n", "
    __pyx_r = 1.0;\n",
       "    goto __pyx_L0;\n",
       "
+5:     return n * cyfac_double_fast(n - 1)
\n", "
  __pyx_r = (__pyx_v_n * __pyx_f_46_cython_magic_88dfdb3a94d86c859cfd2127a1590a68_cyfac_double_fast((__pyx_v_n - 1.0), 0));\n",
       "  goto __pyx_L0;\n",
       "
" ], "text/plain": [ "" ] }, "execution_count": 43, "metadata": {}, "output_type": "execute_result" } ], "source": [ "%%cython -a\n", "\n", "cpdef double cyfac_double_fast(double n):\n", " if n <= 1:\n", " return 1.0\n", " return n * cyfac_double_fast(n - 1)" ] }, { "cell_type": "code", "execution_count": 19, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "127 ns ± 2.17 ns per loop (mean ± std. dev. of 7 runs, 10000000 loops each)\n" ] } ], "source": [ "%timeit cyfac_double_fast(20.0)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Integral Types" ] }, { "cell_type": "code", "execution_count": 20, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "0 1 -3 1 True 1.0 -2.0 100000.0 abc\n" ] } ], "source": [ "%%cython\n", "\n", "# cdef is an directive , telling objects are c objects\n", "cdef:\n", " int i = 0\n", " unsigned long j = 1\n", " signed short k = -3\n", " bint flag = True\n", " long long ll = 1LL\n", " float a = 1.0\n", " double b = -2.0\n", " long double c= 1e5\n", " str s = \"abc\"\n", " \n", " \n", "print(i, j, k, ll, flag, a, b, c, s) " ] }, { "cell_type": "code", "execution_count": 21, "metadata": {}, "outputs": [], "source": [ "%%cython \n", "\n", "import datetime\n", "cimport cpython.datetime \n", "\n", "import array\n", "cimport cpython.array" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Example spin field" ] }, { "cell_type": "code", "execution_count": 22, "metadata": {}, "outputs": [], "source": [ "import numpy as np" ] }, { "cell_type": "code", "execution_count": 23, "metadata": {}, "outputs": [], "source": [ "def random_spin_field(N, M):\n", " return np.random.choice([-1, 1], size = (N, M))" ] }, { "cell_type": "code", "execution_count": 24, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([[-1, -1, -1, -1, -1, 1, -1, -1, 1, -1],\n", " [-1, 1, -1, 1, -1, 1, 1, -1, 1, -1],\n", " [-1, -1, 1, -1, -1, 1, -1, 1, -1, 1],\n", " [ 1, -1, -1, -1, -1, -1, 1, 1, 1, -1],\n", " [-1, -1, -1, 1, -1, 1, 1, 1, 1, 1],\n", " [-1, 1, -1, 1, -1, -1, 1, -1, 1, -1],\n", " [ 1, 1, -1, 1, 1, 1, 1, -1, 1, 1],\n", " [-1, -1, 1, 1, -1, 1, 1, -1, 1, -1],\n", " [ 1, -1, 1, -1, -1, -1, 1, -1, -1, 1],\n", " [ 1, 1, -1, 1, -1, -1, 1, -1, 1, -1]])" ] }, "execution_count": 24, "metadata": {}, "output_type": "execute_result" } ], "source": [ "field = random_spin_field(10, 10)\n", "field" ] }, { "cell_type": "code", "execution_count": 25, "metadata": {}, "outputs": [], "source": [ "# pip install pillow\n", "from PIL import Image" ] }, { "cell_type": "code", "execution_count": 26, "metadata": {}, "outputs": [], "source": [ "def display_spin_field(field):\n", " return Image.fromarray(np.uint8((field +1 )* 0.5 * 255)) # expects nu. between 0 and 255" ] }, { "cell_type": "code", "execution_count": 27, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAMgAAADICAAAAACIM/FCAAAi3ElEQVR4nGVd25bbOAwDfeb/fxl9IHFR2m1nktiWKBIEQcXJAgRJECAAkgQJgPcEIPecfbrnYs/dh7pcr+/DHenGup/0HHsiPKuu02s2SVdpZJuoUz0NNZAfMab/WBibtBQv++z2ujUHen01Bu0p2isazhe/Ts3i9E9+ADkccDgAhxgQQ+zveG306I7NmbYPORwCs2NB/7BXEsO5UXeI2TGJAbiD3y/stGcPYk4NcsNCp7yHb/RhrJQbh8+1u547YKN2WbwjshKZJM9y8XpgoL9eqq3Ercjjr6vKgrsEwy9PZwxbwwRDzDwZAQzIGYym3cNzts9gLvQzg3MwZeAAnAGGG03wCR8wLBiM8oOARgZmhDgMZzDgjJGhJcejNDaQuM4Lvvas8fH4VP6Ua+tMO/cCdHCgoEHFZJGjFU2lgQ174W6YGio1OgblWuFd6yu3ThnjCX+wE2+1yZzEB1Ox1CR7bizLcN+S5ZlDMmavDaLHgWxRXNYfXsfMLJT2xVmAhjZGnprF4QDgxVpgA3GPMXu5kI2JA2QF537vA8dPjwF7rj3qP5VlDvmDlwRKkFGkWcOPAMMMoBUrepXOttJkZwQCxHwyh4LPJQ8mtk3XLSeYIEtwHX+hAMAbYL3voXf2lzvZ2T0jh40OzuhVHiPxfs2ElQZfXC3j9jzFaM0YXMy5RuIYaPbMIcfoHM6MQ7p+uVUr6ppqwTgEJ7XxGMd0xUPnqAIWaNZbPFReYS6BgKrXXraOu+bfISdSDuhMx0wLoWIoC1uVOORU8e4JJZe8vrPwTg8AVa2UoOE1FTTXxJFnjZmXocKHd2nGFK09lO/UTQqYtClyE2+ETcsA4itMLQOICy5uwtUhfz0+SR1Bk0L5zTxx2Hr3QLWmLS7v+BQINgfm8I+j1AvMuM522UqtgjRiKTbQpVzhxRPcKNx7nGcSngInxCTCTHR1gcYnGax+WoiixWljEeTfBm2UT0do8ngXLGwALsJT06QGwjS7NaLUVEqCiUE4vagUuIK9PR/ia2MpFdkIDvxAZYvgbJymjIrDoyGtY3S9MyijalHSkGdGZY4CXcejmiwm6ipXGwyHf4hH/1NYmaJUcUtkjoufJViwy6sscqyMUB2UH0UXxmzpLufGyFNOUY59jSrq9kWyz/NJ8tUMFaQfd4vdTPa+FnKpCVKu7RMVGFGdRcKwA9eCEcP5Q4Vh5NU5xPHQFw0sm+SRCAmFkPJVaPjh61zjpEFWluzhzCbHBXqsBC5N4k/wyY/SshFXfJynBLGDDBevLa/oSQa0pIvur8JRussa9xYVIMBAf/pOzDe8Ndk5p0+OaymdWqrOwKRLkBQArj+zVOaJpZvTa6wYnLZ2aiU15oznlbh9aVJavI6q3kmw4L/Lp3whdBbCjTw711LrQXtFxt0CPLJJQBGwji4tXiMAwci4sKggqWxZSMG/XcQuVC5pcBks42ovJQdd77QR4rBrpmzF3Ogq0hkMMUMjErhkdw9waJltnCDUDbPuYz0hyDCg6lpVOic/zTTuvO3SQ9hI2AytmM3Ce3EI7EacPmT59ian2daDmmMzCT13F1SUicGQrXiQY4kAq0m/qGKXWcf+0Fiy+4N7pDZDiNsBqB72PEqlXOOcK/ZGBEBAmx20nxXJwEDMOkObcPoSdEu1aXiKUn67ZAKAmb+1MKlvNu1qhIYVsU12gTEYs6NEU6PL4qXQvEhBtQvxz47gNkARvuIG17ZEQGak4KpFMOFoOtqmrMDxNftPBV9VJ7Iyr2umygTnkShQ/ol0Em27+J2PvlH27jIPDRdEUvhajfSzjqMBy96cweu0R2XIYhMpDvDUEE8du2zJWUSq754ZUwuy2bPASoqL611+lelJuaAGylljCU91sRRgP1BpNjNRABHeEhRDxwrCMMZ7tgrsZ9gcxHhbBM6nZNWlkyr1GrzQ1vYHfaIwydq0cBpd9deQMt8PMJYhzUSUL25gsttRaGch9cZlLqWy+jzXxd5KgBvDO6o2sXYmCgzIcGeCD9Y2g2xQMaQK4U9hJIBvl3nNs9s1EWWxJV3WlCOqiHb0kRkdu21ISMxcczKp9vIaiBw0N7o0ncw1SudM1QoGN2VVNaGPSjshwZekVJU+jnYSAeLhapNMaO7nQk1Yg1TedRqHc8W0prY08sJLtp0ccGPw3WzISxX/RZ8RZogU2JIvgphhDOY6nZ9j1mF+kTrHGM1UXprWwswncP0sXtbGdpadehITZKOHVxpm6Zm6Nk4Qa0rV5vza+vFpzCgtQ+M/O4zPtFoA4LnkpmIVJXVzhqxlHXasZUNWb98QBD5Kv4BU62ZI0/SK29Xe1NQ6wp473m0Zk3MprCQjzNFa6m0Qn1aikpDKw5mRezf5R8LJTRq8py29qdxRAurifZSkv6dRNanaudgarFLSI9SQtXtQeqQlgE7WSXcqy04RFD7wVo/ZrvdUBW+pcyQjgTBHrCWg7DPGA+hTs2sP12Y4olrvTYJRAMVOK6jn4n+9inXETfoHn4GTm8evp6IWMffTI8sMDXdoHBEwdMX4wnPkJDJTq6abESAnnQ7K+k+Y/oQyOVek2DTH5GJxhBKwNl6Vxel7fSUrbcUsfpQBi2jNr+jZvN4eIZzFr7yiJpclCGcFcYTvuP66uh/s6B32lQbaELyMtKWD7cEg5XyybSCWgBKGThSTrAQzg9EzzlUlZS6BKXLvyhNvyrl9mchUNJ7IxNkhYxeXZuEqcsJGaLlOzwUSqkJ5MJkNu+SBqdGa3/I9ZFcbOId6D1P/5OlIKvk9VPezQ1AHUATXHUC7xT6rkVXrf/3BHPcApSIcqU4bC2P2FJVyysAutJEbTp2EhFe5UhrKOftyueHcq2fnum7l1T6I6BPVas3Ea5mv68yol5oYEYu6t1I0xV749s0wxVPnEyog41AyAfafK/tmHMU+ubpyYc25OqzyTZHw0cxVAZ5XQN8sQKiyJ1vujWuONlEFYQFcwXD65CBe2V8ZUanmjhbOJl8kGlPhxZsgw4Y8Ez8F3xokZV+/1BFnMI1YqailVEySrJqTmUJwFdp+tIsknX84yq6v5b9As5yVRI/H0IiRIWaVsidXepzaVlARjrHwBFNMZsqr7aOnlys/wDnBclTJgcrmOFa+QM0JEWUP6pMSCsTE/3RoWF1Q6T2aKQ+W1WWGFiVV24kQZ+/jH6FvjrufqTvyhFQHi5+ZQ2/h02h+VdeKleuhqq5qcPG5ubuaJ2TsvyS75EHIpMC6MsCstIG0SqDYxoEs9amNvHG6ZVPL5VHA0bOIVo8a2o/+VOgHf6kQDr0B369D0lUrWo0U3uokqJjeWefC0eRarpEo4Xspc9IhlVgy/Lq7Spf1oVPApKnk1I/WIzlaRRTlhKp4crm90QzevOBKqoKarC6SS2Tm8brG+mQ/MnEsaDKjOlcm2fXGK247a67/oAg1gAMwEXfQZr4tZqYy50FoHN0QtcVWulopevemdGPRnciPunHeX+bn7/5gJ7yGg3ggFf3OSPOCtDp+GmFmZuEz1zMzAHyAh77ujAUDeOD9fQBe6J+MWO9JbqFjEV1o7jD7+V0cCriFBtO9YgG/SZSi2RC3vJPbw+k3VNXJH16yWIrO8BMowVL4PEOrn5aHzj6Z9OTG+aMSCFosh8CnY3aI1j5Y3LPAeu95DARXRYSgN27VXg22u3YV5VnB0o7jthGQGA+OhtmElhlL5hBl7ox+68TlyHg+/zlpqj+oElUgQl2q7PCOG3uOlN94jx7cU3vI1FLnGz2LsIskj+pmyqpPrwmf55kOZVW/qjZLPvb1cpyfFI14vbIow7ptM1fcrw/E3buyrFbbf3DWeL/BBV7+EWn8eOhgQt2nRr0DNxyo2s3VVCJZ6qRQ3t3u5BowVQIsvpZ5HuclfAZRwKTjFcQcMNPKs69vMyhY/yHaqdhXVztidbCPu1B0JBN2+/g3K0q3JRV8jqYVZHv1yjJk3IKpYJmTHwPjD1lSoKw0+yOut11S5tzekm9Inlcs3fkaYMy4EZ9mbTPtvib6auXGwd2Z4XpUT3h9uPe+zkiMAa6q8HkXDQTZN6LqfmmKeem2WhvCuJKVXN+PhejEgj8xIw9etQiIMMw7gNHIiNjxfsLodChZVhp9C4eR9rxgWF2mMXzZ69h9Lr66t1HaXSlwrreQVdM6KjVnye5gMOs6o4SMup9r7vZG3FsdZ44KiausnlqUxhd0VQ1A3BgKJj6hOjyjJXOlrmcXosSHwBXVHZBKorOG6m5frzyyPO2ImyhflFktcqT+jUkLFffXvW0zbt2rmXp0ffnLEknNlSFHAPPteSExSJYDCofESLGvYQWBgtwVklq5Abb3E5ER7/epmQRzK8SqxPMbLzs1XuHaoLG1/0vBJ9PUAcm52Te44NPEc+Gwy71eqnVLtHtr4d1AaflkUxUHYzU608Wb822qOXtdwk3ft9HCeT7mtvt7ohDpXPhcB/AeqtsSQnga/TJv7uAknBGy/pTJqNILCghJ/MG+q1ZazD1iqrHLqDWr0U7ZwPGy2zu39C4BOw+XuEwj8s9o8CIFJ9STcZW7d9WfEFZQm6uK2QKIIi+eDxrdkFDlq2AycgCFQqkooWNS4ZPJ2i6gbdpkfnLebxtfNll8WlfsUamOlC7CYsGCFvB1MCB1+B6lDHm2d8rnLJe7s5KyBmWfW2AbITKS9pPGg2du9d0jtjqSPEqHHZ1djurpTIWyNgPf9FGBuaY1ZwQgkPXEtMr2e1waMFqulhILogBqLqIvjTR5O6VaLGNq/YhubdscssebDm6NnlChLoavMA4ds45fnv3K4AJvIiVLDL+Y9rzwQo74rkqMGAWQznRJkSxMrK/8wXqVtZex2guYWHzPLitHiX9prMadkEI+qaDYQkV1BxgXBjvhiLw29nRd6RCXLtTO4zzSLGxThS2V0rbxymhV24feevPxLnqKq9jV7UR0WsU7iEn6/SD3B62NmyfyBfEkUSWJswdGUPCTJHWuOiUC6Ic2lAVJTY1dl/fSkradL0mSTl28FoYYehVPttlZFEXUAsn3l2c+P7kOb/hKE7X+pFEbXEXDQq/AKG2FZEUr1Aovkd80PiIjYUkHVdXId4szwSx/HvpQrSsyZJNcPfiPWu3OgKTiZ7L2CDpfgwTkDu0zq3+FvIjt2Wt5d4kflhMoH90nsOxdpyNHotCd/O0VMjpavvbxzldrZx6HSijbGitKmckIQOS32yVFvHR6tLNWNHbWaWxflI5ex8N897oN1PwvZ/0v+rV+E1ohD5g/qUeAbtyD2mYCOfHsjJRT95APTon0PbH0sMzbl5RwV8r2ZJ4lTMir/WCwcI3HXfetqCdW7K6odqZyrhRJx6vmLTStxf0Z4cN4GttKwP0AQW9Kc5wjat955lyBVk7kRn8oFYjrOs0yXe4AB7Iw5ubML8McV4hL8+0LH45xCjA1Tcsp6pl0C+0Od/Oy0eZQiuQl5UjAFDzUC/ndtbCrAlIYUrQcLldM/s6YANwAsin0ZXqEaFWDiq+Kis3uAkDarnrdu0MdlQTTfWERf7LYrjQKqj1zSVKcOc8oRQauMg+0VQGghTbmQ9iRqfEskTDal+wjDriLRTSM/IvMH4mQqsNnBHRM4nMDoaPteuZgZm4YlMhiclnm42N5l0thqWaVOWWTxioD6QvL7Xvf73IFoLeX0w071bwPgQp8k8Jxzs+OGMc7VQeg3T+hAzC6McyOyc5uLS/bowgob/sF3r+Z0MECz8aXdioa6loj/EYdhP2s5Go3NLtgPdG0g5q3VBr/P6yUOdpLSlgSPdiIfEKHMukSnIddjDbjsWYKmpJ8yZVOD6corK34DEYPgU+xUbX28sMn2So72NExzWjz1rqLMrfY3dsunLvFPB98m+7EhpIXxobNPpK6sCgIVKXnF/DNHhmvGLuzdx++pptfeuGa/B5RnbIL82Cq1AP3dqViv8C1pj+Rw801eLt406pWqMZ3dBc/S2spf/wdGZc/12Wftjkz4WTb/I3cq4yg3vU5XljpXIOB71uvKOG5GoznoR2Tugy+P19pwr+zT98zcb5WJp4eEN1UXSqpsE9Gq16fRWuLTh0DSXVtJvPeZyL2jV5QykiYpS/z7BJhN/hUE1/R7k7RGWNKhOW1FkoRjOeW6CWyVitwGykkKPyl6rTEsB21BFsWKaC8SF3S37eVM2WodqU+WXvddfJ+VwJfBV8sImpCqsSOiNKYRaAiKyWk+fE9KFOq+kYAPmbj56KsIUW36rGs1IjQKkqh2BItyFRr6WN/5DkJfiftMctY3SrZXRjHNx/MMZmj2UyG7LBzn1PnJC7gdRtcvLDvbXgpXW8sm2CkFAgRmToVBwYRRoJR8HHO8TkKNss9j1bqUuqTK4zMOZkzSPbCWD9ejAOxiTEuJz3rMF6N4TemflpT2A45tqYxnPSsbPSB2G7o54p40OOdjZ+RPleEFMZ7Df6aBKOGUZRECGoSHBWYxcB4xgHumylVTqWhqON5xxz3LZJ356Re5o5wFHhfWzn4EAu5Bd+1Wrg/j6pu0SuMPTp7fMMazwe7ypQUdUTI13GoDtZ3zex9SueRS2BNhNmqzcn04qRAAwFAgFty0jAQEwl5vi4w9s+fxHrSoaGJeugkKXTZB2+ykAA/bfEISdSnwyiQCXWrfM6m/by36jdPFq3zR2oB5jt7U21PcCfQsE7MN+pYT94mTByyQExpPuqFegyXU2oQ9/VCXIq2DEwDDxV5kfBJo6fdUcR/TqL6ED+XDdoZYC3NbfvW+M8+kgelGyWYV31aPCQMM6L5M/uSM5pjFKog7zxGC/iLp1eUxOHoZlC1b0qTLCHc0Ippq0zehtEZClNtZDi6Cp87twTBPr2dEUYfrZ3eiVdwLZ9hzWVBbMYsflCgAeCTh/pFIXaex7MfZxqTukjp3mzzl6oaNUFssDYqYppsl+0CfcG7R6P7LsxfJjRvc27a/pWDTXeBM0Y06RSTp4fyGyGRX1deHZkK+k1VLDyX5Sc1qOcVZu8c0Ml6L0omqVA9xKaC9tRVigp38a6rJtSjYYslV6Didtcq6JhY2/XvIexXorQmKVa+yVI7THaeOwrmxsx4JUTaKlcL1mvhez3161mwjKzXNEZAVrUFeWSdROfkxXtjH/YUrYgKdsiiDTGierNwAutpEe3PDhDDvz5N74McUJPVT1mArxkOPt1aRE6iZnaZi0ZqhLx1775D2wxQo1haRm+YTKU9qm+4DFyfXm1QQortj4CpNamF3pky+4f4xbVXLucWC+mzxZuu9PbdFX6Mp1VZuoVbL3mcEhNXpwDtSc6V9hHSh7hdHF4tvx0mLQec3TKdwsjGJutataIwYljQeZLkFEQ+dAPdz1uUFsZwYzXlOesiFVjq9jYqst4WOMmqe3rxjW5eW+NeyEJ3+N4XdzuSR7q4INwS8JNBB4FdztW1sRcgMPgrz0bl+9CrWu8Pw1wVgusBvUQTf+WT3tso0EnsXLYoOsIG6lGydPKvKMNhNdQfQjAvOK+aSnpXxcgcELqn0YNKvXq/XB95DtIBuYjwjqP6pZEq/08myDnewxLzi1mYFNFr40DzvtyOLgweUpTjSosnDiWE6Qj4kX6O8HEBLXbSCxphSkCJsYzF2qlU2ByrqvJ45FzPldA3Tzu0Fs/PLp+UiCBi4MAqxNa+xSkNQLDY1SoCUUkktOolbToqAof4kqDnS5vo1sBLiJccBUtareFZjP1pWyzwC97edKRnSKETEvxaTVY+gu77VSmIRY6uWdmeBizDFrdbG4cliA/ivSEuchyVR+NFkw6ue7lsOcI84elqNT7xah2AoWQ8bCjtjv0ivqkPGXlQExGcs8cZt3ngWmYubqFeumrg710RpR9pgGogqGFxZTHvAlzd5Hy+9Jy0UT8kbuXSBr+dtqob+l7hms4K4kqhlqDlHvGpGh3eQAfxqPHeErljukHbZqtKJ+KoupJtxjvZCtN0R/2LnnYwW4XnV3b6PHhUuoPpw0KET4JDFwkeI6gr8xFXCBrDUUfNrHPbQm+n2a3+HzVAbX1Sauu1C5d+7VsbdBGPRrhXxE53iX86uQTAGTfQxFfNl5TO6gbz835lpfAfbXt/58jQZeWyDzfWZYD9LzVkMZEcpWsbJLDW88bc5Iuyd1+DVwqp79fy5iVNkPUO3Tls/M2pkO2nhtTvXBZO/HcekFeQUsHsRO5sR89UjQyVBqY70ZQvMu6V8nKmLy/WUT1JpVJcUrR6pAgtRKCpfEE1o+tu2D2S/6eyQYVQhcxj/1Bv+mYoRZKEJs1OSLeeBTgWKG4QhZqwY5XKburd9asJr2bcjbevySxoQ+tlaFwPI8rq3QVmNU1kWVr9CdOIZ/uxXdLMidiYOZNn9lT4TpxiIVOkJH3haFvcEHlSD5+Qt4ygaag47PDJzNbyU0N09z55sWGNAfCX8w/jJye8FP0cu2acZ4/zie4ylmk9BrLO5JOvmFgkJlZMbBXHW7NalP23RBTBWKq9FScez1hqq9VS3JztKjVnUpWFFsQXMpeC1qalR5XgpvbL95pW177J4RRRyguoUCYoG5ywSc67IBTg8qtEEF8Y20/OVeJmdnrQkJB/KW+KNWCzfSwbbHmI/ONjfE2hB/U3djMneNJyR7OlH5Q6yXSPwfgA3PdMumqr2liSbjyzh1//o4U9/hSBTTKufg5v4t6QOxTGGF3NzQMXORLKihRo0oJZ2R2fBkQJtyFVgRSifIWRhYevHbCX7RvB9j9ydmLk15qRNa9/Fj5RpfMpxSJMpbWz10+7JZWPz5Fqv53wTtBizPPlb4utvdKcEf04opKid0Q1KFEtkOPs5HUSqlIEcZhei8cLxEG/HZ6LWBMGI0n4fpJBHzbZNfztolyYMJ5SklPcLxjvEKL2rTHn9rv/6mLVUXBBfQqt+F3FUNUB1eCWIXA8NDIviOSHmdLDl6D3SOE3DMXlp4EtLzeZZTJ8Mm8bUxWP3i9bRkgrHpIShrG57msv/Y+KSp9P2+hwV6YHPwZZwy4AMQkYY4ZdoJHczNnBoqbrTA6LZBiPwB/TLK8mSVcaGqnlT4FWlucKv1hx+d27kVJRb6DSbAUjCWUVoHFMEtYS2RgS39jLDowW7GXLyz7sSv4rye39yuCXKl8EKC4Od+Ziz2Js6LpCyD2wuKr2phLLTSCVLuJkU+lzWbX0HaQOqlUUKtClwEyn9r387s4u241m6xkH6nQhBJP8srgT7mS1TSk+MDjKxK4uNb76APnx6SUkSdNSCo92aTiVmM9y9XhB3tYWvW82rkq1cuNOv02vvWJj7rvYB9CtKS1vdy55dOSq0ds3oioIUjH9aPneyDEh7Vtvu4Fw9Lv/F4RrgW8vwzpAJHvpSu945J2yU1r6tnAQ3j7294Q/Up7uqgiyIOwFEKoRd80txrCZmfnSuSlBXXNKZszpvPFNK7hdrt1yRLa13PxsPEVPtwEpUTjwLcrqi7D3nmzlNcwvlPdVfo8eXDwuEq9zLS5VnqOxq8x/aTDND41c92nvuC2YRAdpmlXMzRI1a8Bokn6GMTmX7DFFhnK1+JawbpVC0mLEFKkeTqdKxYZXRardM3VtFe2bmj17KwGC4J9T/7SUQilpJQcRuvcv/r8rlHUXIzPwqbKVYBRGYdLibQVUf62OiL0JESZuPhZGL2G+sZujnESjy1m3/pH+9L4Cj9Nokh65X+i43WSKaJKHUo3GymTXX7cN+D2N4nvaLaIJUWlQpEwqlqongh8NAKPgjXhK9g/hZR4f018D5cEybUSpDCdCTY79fq1dqPRkYL22jd7p2bH1lpNYcFwazaZXPOjbPM7zKvVR0TsjriI5LXHF5NZGYdVktSdOD/RwXNdab8MApfdEySJhBdsH5cBmvLg1CkWsqGldO8anpnAjGiKi3Uog1cJx/qEHEUdA9WJMR409t9UFu2Cw18PyLko84j3XkzLPg2BjDAD/QLRk6z+8UnWpPWHDq85e5Z5oPBy3mfLWfwr26ASFLgM3/WXkw3Ee+Ps4Lfw4cc/jxnJDqBpKHJaTFZ9zbTkqCZ8DFbya1jP4nB9NX4ShxoEy556Eod/SHbXcLzh4kf/zQlfYtyBWBVeeJRNYAUASBsm26NtIYGUTbJ9n+NkgyMUZvmW925PSJ5IPTzbCSsLxd5rDjkvy6kREaKeviAh6lyO/7f/r7d5Rf7Jcp/o95v274psgchMB/Y1VuDp/FXTRN65V8Lpui4O4IV0BfOG5ac84nR3/OEXXBZ60SCikkFQx+bAKluWRrzWYGSQnxxg3OeN6ML9On+Y62dkHmxiyBa8fgUGt2+fULskZBvAS0H2gw9UP51R7RA0flTY6Py0pgNwtPBLUUsE4WWMuXk1FOkG8A01g5u/Amzp3YycpOKYEXxcMCU/LqsyUMBpGcRohR7BlJtAAcz4R8RhgQKruiTWeIlOaKcKHokfiJLRmRj+DsPLQckgVaCiiTtMpJUA0ef2+a0PJhmglnaWaay9iUH49CVJ5wLKpnEAbH9vbQT+XeHR5MkmpZE3mxqfwUy8HWbPcm2Tr4ygbPX7cjKzVtc4XaNj2VU/tYKJO8jqbTFJcUQPti5+gLX+vxqSo1ZWk+wPZMbfhMdXMeMcwHwmqCxbORZ2AE/vUJ+8zK9yxvG0moj4VrcGBrbDf4/rB+Ct2M3uY8qlKt/ATP5eFlCodVbztEC3eud7JHgTGracs307xYjTKBBmvTkOcpe+DiaRg3Jk2H54vcr0EdGmHMYGo3JeQKXXg6IAXwR+4KwDg0wi3bLK8kKaQXwPgQ6YLTYpX92NOjefFAm/AE74z/XmwtwT6VB2ks4u1Tl+XDIG2g7wVozXneXYwo9Tg1uv8ZLWpoFQj9vgM3u750eGCOmjQuGN6x9dpLe+SfigrBRUIGz9Cj56b0sRaq+GkHLF7CjSng68o8jcdPVmdUOrZElWjH+o+CF3WqjfQEYYKuDqXoT8JzSvPvJTGfb5ogHwB9Pnm8plSvP4CiXENwbGCYDQdJSdS1fyxDMbQN2fKlougGGEhdpLAUvjoiGZSKZFC8yliPb8PHuhdR1GkxjAaRK+efHAfv5xk8cq9u3VxR/tr918m5Zmb2BAYbgt6jhvlB6c69I0URs02y2Lpc4kyQ9Ql5ZWbateet9kxcM5P4m35N3Y2SyXo6TwjbY584Me5PIXnre+AFRlElVJH4c4c0/zmRWsbGW79jy9sqguGV9ckpOe2IJSAFQUvg2Eapt5hS6Q/+mAmcrncAU6X6y7XjoLm3aI/XhiVSCPWchPqooiwsdNcPYS576145gkrjapjXend8YsTnbgPZ6oqquYm+R/CrB8arQrwAfO1QekiVlQ4wg+pGx6qcF37BCk/UTznTbMsYrClgwvUBdcPHwckVaswP0I1mWQN7JRSvhCl1TO0fvmCbm2dFFXR+WRLn+QZ6Rf2t4HxboK41suHqUsmV0XiKcyBzwNQqMcNl3m7wvBO2H25EKeRWlxUoAH+Hayxyb10vem6gaitGeHQog7MsPrg3n3RjyrXaD0XoSj4GjA5F9COLRHSROau+EqgIQZ/JXGS0FOOUr0wL/ngXsItDrzPkHiWpy5RwS4V7WZBa/OKw1SKGl2/KlHyv4fRasrplfA/PFIX+ZVQQ7JyzC2V4G4Bwkgi2CLDRwuWRjNV6FpjNGTyHTW/oNysmiRyFRJnsjh9/IlIun6IFVYHwQUKicK5dgTEu1dI3RP6syH6rOSimXtw39nRmyROv1f3lv9FpqJqJ+sTqsrJTnThXwEoVI6Hn1yfimOo2wMMYOxH/ZoHVPgpgoKDTS7WKFkvlp9nMV1XTIslxc1lykRD8o2qqQHlUqWE3cr5sMpvoE8bbRQZ4bhAUAgXJcpjC0Qvhnr/9NYzuDsg8gkLCVIVKVPuCcGwhMaZqGbTWhYIEPgH71oP6dx25sIAAAAASUVORK5CYII=\n", "text/plain": [ "" ] }, "execution_count": 27, "metadata": {}, "output_type": "execute_result" } ], "source": [ "display_spin_field(random_spin_field(200,200))" ] }, { "cell_type": "raw", "metadata": {}, "source": [ "### Ising model logic\n", "For every point in grid:\n", " energy = my spin * sum of all spins (+1 or -1) of neighbouring points\n", " if energy is improved by switching:\n", " switch\n", " else if we are particularly unlucky:\n", " switch anyway\n", " " ] }, { "cell_type": "code", "execution_count": 28, "metadata": {}, "outputs": [], "source": [ "def ising_step(field, beta=0.5):\n", " N, M = field.shape\n", " for n_offset in range(2):\n", " for m_offset in range(2):\n", " for n in range(n_offset, N, 2):\n", " for m in range(m_offset, M, 2):\n", " _ising_update(field, n, m, beta)\n", " return field\n", "\n", "\n", "def _ising_update(field, n, m, beta):\n", " total = 0\n", " N, M = field.shape\n", " for i in range(n-1, n+2):\n", " for j in range(m-1, m+2):\n", " if i == n and j == m:\n", " continue\n", " total += field[i%N, j%M]\n", " dE = 2 * field[n, m] * total\n", " if dE < 0:\n", " field[n, m] *= -1 # switch\n", " elif np.exp(-dE * beta) > np.random.rand():\n", " field[n, m] *= -1\n", " " ] }, { "cell_type": "code", "execution_count": 29, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAMgAAADICAAAAACIM/FCAAAauklEQVR4nI1dWZbsOgqEPHf/W6Y/xBAR4Hqd3a/SaUmIMdBkX7f8xPtyi/xjZt638eNh5t3Ewzxe3bzxipKQhZl51C+hV+WvQTLybkLlrhVyXY3CzP4NzemDCagc2VvJ6kelMEv5srhEHQrA1QhgWA3puvQR5uGRUruH+c8iIliVee2mH3dPNRD3zjV2q+Cf0K6MygKEmSVTU5f04C13a72qu9xo50GHWNoCF3tmciomTl4BkHCsCgpYjf0sYf38VLAUvth/+nN//4mSXrtoCuVO7fBi2OUfI86KHzCcB3HHDPh0he1boT6KzlsUkdNo/CIW096eM5bbHqsmaQO/xhWT3G3V6OD/YQk5OMmB+MAEwcuT3XbuZ1pnAS/8oNL2A0ezVEFS3rHooLujO4/lGiDReKk0RaQcCY+62JmbBUod0LCij2kU525mP4sHNW5uEe3zECbT3YOrMLNwd0SSDI02WBCDvkAodg1iNExp9aW3hgJbjMrPVNVOilU20osHCtiJRhyVzebJXKdZIyERo4dpusUPReYocTPzp/WAGPBSfXIzjSJQaUOHHRDhJUyzxRauqUSD4Gg5/4vYFvGJBIlSG5Xg7wVkHC/g8V8AZ8M+5xYUo/9O+8BU8kMZBtLrIvJ/BjVGnyhHaOeIDGHm5uE+VoGKEUuOCUGUjDAzop3BHUMgRuLWHnokEHGIBQ5SZbLaOLWr25C1KBJRXd0SiFNw5NWPC7K40+XTicdU8kc1IlQGsecYCSAIwzgjLZNvDrC0UlE8Mgf2moJYgNZR9gIuTgv91bnLUDeAtR3+jc8eqjY0U/rdQhyFeWrvKUi7X7FSEEWtIcUgkUrma0yA2YzZDAsU1CdqApwqBpK4/XxDUjTjwcVqAmrootiVAsyFpWHqJpxEgMo1jBvGII9oms4kE7+HJKREgKAt25ZjgBVGohlE6CY4EVkwUdfl40M9zDyivNtjqSwsHiAiuosqVJoqQ5m7sYR+eOgdaHWYhPI3mVJmXgcWZLCT01XtwyrYNCU/CnPyQpy4u794ACSLNCGx3FPUvzRpE0ESOF8NtCNQC9lwQi3vj/qDNAl2xF9QB1Y+8KdJY0qgbvajOPD5QxBIeSWJsTlcLnRKrdVQRR9YWB1CpzgsNBv/o5wdoE9OwNMvmnsNUYnP8DX2k4iYdK6WWplPFTbTPLJIA3+sZmY0hnpSRpLYeIQgfPk13XbMBJnaV94H7umqU9sb7wUHu0W88NHB/GbFkf8aWfiVM8/msZkL6cmB6YPm8tjfrlUGD+0R+8UbE1SVjqivNvReEiCKPen5L6sSC54DBYyRIzQAz5cA0CjaeW1dLa6aU7+nHx/RhVwcXvcvCyidBmqlkHVAsZisomh31xWTk5eo9eIQEO8Re4/vcaRQ4w9y+7nT6nNwo//yc/vSTHDJKEYoeuAYTzU+gwRRCo9BaBJjv6HDntd/axJGRFc+HV6Cf9qKVjdagDHlV8YLOvWEnnz6CgyBwG8C048Ev6CD3W/ZF2toVpGOfEiEa3blfYz3+eWIxWHEO+YgNM6/0/nGu2QkLAiaa+Gw2Q1VM1MxktPHKs4erKpdDjxlt95sVHeMyPsntucLqoc33ajzO/ML56/ND0IKF8DRiZ+KJOFbWC1RYjVAP8c2E/cQbMCLUy1bE7I2yVgJiLj7b7obfFeqdihS5brq1b0Yf0knm5w7uE3cP0g/8gX2ibL8sMAXOp3cVdSAQ1iYJO5kcjhGSgvsws39sAJn44gIWHUOhwmOI5uaU6HbI2O4LcaowQ6Ac2yfOklliHIm2ihNk8reL4yRzLPlCTUU5ATi2Bq0UYu+26Dj19WW3DvpLFvB2GxAeVOP7Ps3qcTry8VQyBRfOvAH5avP09AjR969cm7UslgQeo0Qr27Evxku8NIQWS+xKGycQHhh/671G+TqLYfMSIimwL2mhtkUJrJ4FdMWGHZ0THY9UgyUhmqYROhGb0VX1qt4rb9ErAChik2tmd4XQOIHWl0MOXrpFahmmmAYuY5VqL8/wXLQDG88IejmlWwFFm/Fy1xOql6JvlUIuLO94EwaFBfRDXcmqzzSgHJl1JhLYIKjAnWZMYDLH4tBF/8Aq3s2AdNAplwXr03BbxtNrODAQLETbyG6ahQix7v/VwStD6cx1f8H4pr5Vvivx0zT/suvPVHbdDlgtItxmqu/Mw4uKo3zI8mUyaETBK0iNhXe53nB2nzawyoUkQW4UjW7DAUjkir1N2q6JrQm4wWBKVjIBNPCaKwFBCACY76O8hjvl1XInAd+YlZ0s+HGD03WHYzQnGqRFjvY/XEVAZYGaTBWxilgc7RFiZEk7tNcgwz/FzSP33ZPsVvPHiJAZk8P1omgwpa+75Qq+FTCu7gwbkRCVrWsjpHlkQuCXIy8VjJvo/ZAQhPwY2tYcCC0T+94144aP3ylDI3PY5R7RVwp4Xn9v9cPbWZQe9gfMFrvZYXwjubCCw31KmrDoVawkK4hty+lpGttw5dXkRz74/QFlBGIzeweydynNWSpgYqGG1mlJNTi6fTwmJlj1st1iR4pYq8TKQFV5shC9wvmZGSZeXjyVImGYNTFfmJbWjplvtEdFDOdisff44o+Ufr4/QSv1CDysxf542rs3746ImMkiuaWq0pLWN0391iupk6BLdGxwZKL7Px1dRLVU9LH/nHCT+UzWCLDhxtNS/ww4rjCETd6K7aC4ycV8iLM4ICCgk643M2tXLZhFKHOlZW+vQdlzhRJOmaMToLsur4sYrfuFxRWA3B5hS7wYnV3TV37UCVUJCcJup4WdMr06JE+iSZeACequXAWDLNvE23voCvtRy45QMQhoYD/SNVH4hXU0hMP3FQOOvwRcx0zalsIK3IMOUVxDCMwRlZIuHNPdBjUawgwAHPOLsnfenGHQIfZeq0OcBs5WB6zPNMIkm+yU7S24Acvmwz5G5zVpTYSLF03Wku+3UV/qfXf2fgjKr4yb7spD0mhY9zhguTQT0Q40dnDkGgQw/zvM+gcytj2FwY4HQO5apI1vbUeLSgvlAJTzA+qXBt8pSTuLuboHTwkAu1+Ngwh7EAgtTsTYrBLv4qVpt09SkHhQ/IQI9qElfbFGI4OUazAb08E+re9DXT0+heNozPO1rCDG7WEvU6BqBPYmNL598Q42Sivpuc63o3fLUc8tQL5eBecQKJSLXuTyS/e/xJ2wVBhdBX7uqgSLI1/E0ZMuo5JJL4NY4JBHH2VLkCc8LdqMCOt0Ykk49BUA64vmb4+vc9O5xFHjlpRQDm+ZiP4tIdP0vEJuufPqF3sm8gGpm3KfffQw/sQxUA5J/sr9QPz8IPTUDLH0RSrRZC6NZNwx22KuWhqv7rp+2g1k5JlyoiKVY5ywqbQYQtjn3XQNabKqS0zMLRm+7d8lT/SIrhMTut4c3sYmSc8CHOkt9darBhTRgaVVhUJLhakttCCSSIxfabK9Yyy4x4OHMkLHAkRL1TY1DLTAJNi3wr7qAgTyQWIU30DUQvUOVwOTW8nssBrKj+MfQoC4nLHi5Vf7u0ZQNQhh6z4wzM8kwSCShXGikE+LVEdOVCs/jJ7VlpjOXj0N4aDDpE+iWLkAf1TMCdYaIR+MzHkhpVnO9hU1DQyXjotwrzzSIA2qMfKmg4FjtiEYe9mBgo/Dooud4XDJTWxmVNYwYDY8OazauDT98iJrOmU4MwogiNobQkpFvgDsbCfvZZAzbaFf6XumFPPPRw166eYNrwD1yjPDLL7NtEznCcevLY3fey5fX3cYOzDSp3fSwRAl9hcfNgrumFmKoBY7QGN4Ncwb/H7o/umeirt7F40IWjrXPyAc0Q+F7jZEekEggPlAjjTvM/dCV5zcm1O+OiqYog0Se2C80+BDlpXZBB1rDu3zxFEOW0IYz5JvQMPWhcuhm9eesdyDxC7h0NxCtVSl+P684nc/VKA1bmDWUQxxExITmK+iYvVGJ7kqdLbdliD1jwhRjxdVQagmTWiqeKToislTzXs8/xMDDUmOpeeckwNTwaecuMyrEwjpBNJGWYVGJhz/DIJf4rOis4zw5AMUCE7+lUOvUcP54T6/UATOC5ptdUccBVP+lZD59tAPdjwZseBVlm3gncqXIkhEF/Z2rjLwRsj11l/sg6OhAk0PtLWkfClpke+hYPBZKnJHRNCPS7TAsJxe80HSVRSuQItc+BLDRVEcI4lsO31ugoiq9Yxo8EX1AG4lCE/klayXVH3Fr9Tvaa/aodPvSUKnJJpqFS+ra0MHzki8/SZMev3mEjhKqlYwhJPYCLvos69GFgE+iqWJtZHwGr8f7FmzE6L0EOvIpKSzV7VYWD6GWbmP9UcY8FueoxKvfEKPDf/vpCTXDf0eeNrGnYoGNqtqS4wvEMTBWlMG71vl4FlCECTbORm16zRoAo6Ir5SRafbJ5d5e55DjLnSl5hYiN2cNCYJTCILdK8461MEbHrtsKipNOnUVD1akDlGLmkIyUNrNpdZdApjlqZtkPj0NgKp22c2pGBTLqd6eoRnrIopaIxY6yzYZmRGlpBVvviF2+2euAQM6zNeaQKPLUXEuxs+B2htHnZe/XGeGZ9udUDVGBsZ4hJDlkkVyi2XvFWLT7l3/xUjz+jISK1rKKKUuA72aYbW0x8KR/kb366m5WHb4iQcPj0tblaj307HASQ3aO6OCvQdsKv0v6INE87K0AyOh0BdZToGyr/yNwIS33Ks9rsf59+Y5PbBhGf7jAfoCiDhSA4KPmXnXnyofDOjjsW64w9i8jIddhoVpyhuYeRxrLWrbczUzwQf77N7cxblhoEPDbrlmdpOPWJrOiuCu+T/wdDOEitjf/kCJM+JPF6kFWTYwdw9CKaE4NkW4IOr0qM+dsLg2abHY+Q/AYZCrL7u6JmVSKQcaEMg0vn82wYSUqnOp89CzKh3BVgEblE96XBNfk7kQhZANxbVSwDwaLmvuNaFGhvZFmkbhp7JBZ75yOzPcG5QlAbwOypX43ebo34ROkSoe31mxdNwPBMCxb/fJUeOPoa9tKH1DKL8E/hZT3QQi3dmgXqMY00Y7U96xcPSuFZsJsvFQ7Edik+ZRvAwCZLquzkhdA1FRwCnLhEZmWTAyPolsz6cHY6Dy7hmHHjW/5+klvnwuiFlvMyfBLXMYrL2zvQ83kYzX/DH6Oej970E+yqiLn9KF3oqx+WMRWMw74O84CGDvnrSYXcFj2pAADjQC/QODkDRQTWAIWBN9DDDdBVaxCLkWA65OuIje16m3Fw7EiQm6DN1xiI+naVjBAcCnDPwD26L2aTlK7fT+MlNXJK1/8eH+HKzsPh3bgeB8/BbITwMhh5VpGdHBpHatdbD1LkAHzPGGHGh/+NcdmEB3/mRBcxbadVzRXkakVMDL/hlAiqlOvgqCm1meRybuetcNQcWndLkoxuyhpUxj2dMAJUkY7AycC6oZ7XhlPuwgFkTPX/kvCJgJV/zVXEYqGE8b2+MGEhsDVBDmBQRHP4/3DRRSnitCVVhrJ4Jq9o1jv1sOeTS96RH+aHPd4fqvU0RbDDo5om2wyKk3i6pd2LvBV9YfIdIlR1zzlAQP2vbw6Uq9sKrmzDE0x0QdhjJ3g6JKyu0IgnD1V573WnvuY5J2qBMEwXdimff4zBeVPJFYTVJn2ew5q0z3EZha47lIjOV/jGchV7LQImlaAyEAb2lvpjj+Wt7WlCE4mFkclb6kYxK/3vdw/knd2rOXYMd7txen98eVY5F+6QwiwpPdPRpUGzdPtT9h23dECC28uBMxAfTujQdnQGHfStX3y/EakYC4Cc1LtvnBK5BeEb8SbRSSDSNvf1F4Tc03zuxA/K4mLEt4jOYoMQEi/GFMhEkx5aBlnYttyP3IxH7lULDgGpDMvOKi0Wo2nOKyWYrQxyUOdKAO2buo2xHSt74La5XMjzZq4jlpaI45DDVEKGKVAuqBYsmX6k9nSjgqX5wfbwAppMbtStWlYQAMp2ZggX7749mEe8vfS1ahyemPeABMuceMJqR9cRNPyT5DzH+OJJWzvxY+nfXgkH3pD49uqW5ssu+4f6vJ0Tu8/BrtC86qAU8eVZ3PLQJvd9liUzPpNneD+d9QnofBO06YhyM7hghTe9D+Ee9aj6PYDqVWtvJDEzWxScidVVkUBMJmvMYD8CtKxFqnKGXQ9Wasx/5/5PgmGFKWsWeGivfHBfBJ6iwmWhIoPb1hQ7Lup5lBDnyF1CpiarSyHgouqaQrQD0cp6ZbCeAW8Hdb1hdeaQcbmC9N7TG8Ozw8xkenbs14y3iqz3gIhzZcxutryf3oKVsxdb1RpcZUpXaJlv4GpsZLAZAYvpCM5kzsxE5rLAE2I6fYMIhQH1hXw9868U9hR9ZU9kdBx3ArhCqBxehtubwG9DRrerJlO7ItZapZl0qGt3WzqDhymRZeMcmiHEFqx16X3nk/qRj8En9Kcvbbojsmw3lxvpoF8x76uBorbX8lfpF4R8osgQ5DGxwJqFXWhjRSt8HrIJi8a0P6igYoN1eo/LXZYUUEa0NObbTRGKmaXWcFFeGKwApBRIfSL+qTI1oEvXk77t1Qm/fuy3XCgqrQ+al9n0iRRBeaH8gFZUd0XlF43Cx6f12P4his1FUFAJ03PBugaPFBrrJ/gNSKvxHCDxlhWOT7v860yoCdrfl+aw4rGTKxDXPv8zB54G0kPlw+Ko2L0Ic+DPLs/GtOcrYeKd2uAV3M5HoaarhiaPGNGe/MDtealf0CyXh8EXJwfmF9nAZ6r8ypH0pco20rjQmb9kah2Cr931dJV+GrOIfVOJucQw7j0I/2nD8QP2pR6q+CXYodadoqhcy8hbBdaSef02Wwlaar6FH+kW5npSqy7/3CABuEUIXmfUaG0x3SrbW2fCYE0ALfbT7ojhKuap9wO6DHTpIJydWDihzTjlSZPYGjVU8CexzSJG/PdUnxg69o13m14FlYaRAH3a0+QfhZ5HwahcLmnQf1vydynG4pVz1ZDD6d3TkxltRxhMw62R7d31/Qr57AFIeTrzlLKMqYEjgSlDT5GdJtRl11USUM2JAWs6tU75jpbFcYJ0rgrXy6BvQGUTusAJT9dduHmW4ng3W6xGgD54cjklATjRwJLuQog3Tf6CA6LSNWboqu/LBedchy5/Z7rRIUahxFtRgPIPHnG8ixla6arW79xykYjUAgfBYdpalSWeccIwNwgzeX2IcG82WlCFzn84vVFybxNQ5LINJZ3EsM8SWtQY62OwA+PFB9PflJnpumvpr/7fKhI1iMMNxfIK7YQpcv//o5woKHFao231gJQGxYB892mlQSdDyMGxgtk+LTFIYQbN8PyYYeKWlZL4ZqDHZDiztjZRsY62MlU5cOYxzqOVm++VuIyt4LFuw6h7xRVK5/DbIVhhjSb4GcYH1kC/pht7f8kOHbsNwLgfO6r1dPuZvhSO6UbrLXpM9FyByQCavpxvFUjH8jd4D77x/3I+VFIdFqpJilwHG81rKQL9hbeBOBDVje6m5aw8ZW0DIwQm6dTynfsq5Mta+0x98BDTPFhbqywzSo6t3ILTofG7HzPg1ZH3dWWhvqCN7j3IsqgpWQffHLjSr1g5g7EBaQeQCkadNcTaeYNbD+C0HmIk4kZnE3Cc58G/XhMCrofY0pdVfE6jwatMWKCNnj76fkZMuT32ueTQi304nIX34x5ABeF8+8PGUOlT62eCGNDWrlb5GhveVJjn/+RfFrkcOX+zHNkbTba06DveCC5Yo+sKpBbmt0/H4/Q/knThvfg27+iakkNWrUvpw5uH5h8lgzm3xR/zX5gShbC1tOa7F9PRHGiR899p4VavL+DDTLLnIUcNox1l7IKI7tUmMz9GNbRLfvz6zPZXFfMPyfJb/O053RmaPCbTTtChPoUf8gbnFSnSkrcMAkoePECRycwHZMsAMobWD9Ag4EsDgt/aZV0zIun3vVPVJszbxD7WG/0YviKofCxhjvGhQgmWc9HUgk+kbNAP28mFoBdjrIczwBDv8p/A2dhHPtLmfZeySmss5iewAONCLY+jonWWjMN2mXv5BtDpNhcs5Ip9yDNcRyoOQnTzIG1PCtrL6Mf5Nk+oZNYRF7ERgG2eZZynxYaNZBFJYJgUlU1fYh+YQKOgvsOFk1v2+6mr0ywqqt3RuJE9s/l1wn+W1nIyI/HXYxxowMWLS9DURQhkzq/hUbkcXaJJ4YBWonUHHf8v5kUHa1w49efc0iFX5LmvQ6FfAF2WFAviCHAtrUqc7TZupjhU1m3196D0Ei3fOPaDxlSBITQjoUDOWSqi350NdkXBvDQ6wmSOfZK6AYM8nwHOXvV4Y4HPa9zVXwNx0H5Y93fc5CHgMfLZ/ixYmWElg6z2bId8ty8qQ57iWwR4UheFwLVqOaPdH5yrU9NmKHkcJI0+B6v84fnyCft6gi8nGmnioHD18cp5B42g+1mhUZORVt9d4VBStLnKTsvb54JvBMMhJBZwbGfaWoiQpxHIUWRGY+Ay6Ix8OFbf8N0M/ViAaY/gYiLyRLqA+pQ/FodkoMTNYiHx/3Pp0EWw5eafWW445Q3TFhAoDFnEjFce6RafpwJ0Z99or66C0ROFUMztMokY1e2dR5PlPAHs4JDp/YRPD8Z8KaLL9zwjEIBerRjyM5ajDQOLRZuaNYvrUKntj0eJo/FIslmsOUJucFIF3jSeNFGaWElhbZHHW70dq/QxLfqzPt0tfufadmJYVaU3qVpEWmRAw7IRg/1GVjgLRZ4+4UdeXUuHfON00DWxAYUENEAa0V0wu+MnTQV3rxK/3JhRHnYgIYLbotJIgFIV9DHeWQ4mGJT00XNd6QnL3fnJ1LzRk3VhNMJsO0YpUzFIosfCDWxtwEOXLfifOrl3dPMPLSIao5ThJvxOQyQBtTD0bwX18CK3cm2eIcvhmh0ON52J0tP/VMnzBDtmJJr1JnhM0JRtdzg9uMY/dUMYub3cWYjlMMolvYIahbVR+qBXz87V70fiRbJQ6+ewPGVSAwcuZKkOb1MEsFNWGeVDhSCcXlmpGUE8fkgg7aDfeHuOB2l9ql7ykbuQWPjHCukjPZeHD/ogPN9PpBOj1CtDv5A499MHik27SDot/d2EC5ng3poATNboZOJTewwP8rfSXc8ZFh3pu3ZFpFKmb1v8AFpLs7giBfDcAAAAASUVORK5CYII=\n", "text/plain": [ "" ] }, "execution_count": 29, "metadata": {}, "output_type": "execute_result" } ], "source": [ "display_spin_field(ising_step(random_spin_field(200, 200)))" ] }, { "cell_type": "code", "execution_count": 30, "metadata": {}, "outputs": [], "source": [ "# Animate using ipywidgets" ] }, { "cell_type": "code", "execution_count": 31, "metadata": {}, "outputs": [], "source": [ "from ipywidgets import interact" ] }, { "cell_type": "code", "execution_count": 32, "metadata": {}, "outputs": [], "source": [ "def display_ising_sequence(images):\n", " def _show(frame=(0, len(images) -1)):\n", " return display_spin_field(images[frame])\n", " return interact(_show)" ] }, { "cell_type": "code", "execution_count": 33, "metadata": {}, "outputs": [ { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "08161af3fb2547b89e91e5acb9f80ccc", "version_major": 2, "version_minor": 0 }, "text/plain": [ "interactive(children=(IntSlider(value=25, description='frame', max=50), Output()), _dom_classes=('widget-inter…" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/plain": [ "._show(frame=(0, 50))>" ] }, "execution_count": 33, "metadata": {}, "output_type": "execute_result" } ], "source": [ "images = [random_spin_field(200, 200)]\n", "for i in range(50):\n", " images.append(ising_step(images[-1].copy()))\n", "display_ising_sequence(images)" ] }, { "cell_type": "code", "execution_count": 34, "metadata": {}, "outputs": [], "source": [ "%%cython\n", "import numpy as np\n", "cimport numpy as np # gives you access to Numpy C API\n", "def cy_ising_step(field, beta=0.5):\n", " N, M = field.shape\n", " for n_offset in range(2):\n", " for m_offset in range(2):\n", " for n in range(n_offset, N, 2):\n", " for m in range(m_offset, M, 2):\n", " _cy_ising_update(field, n, m, beta)\n", " return field\n", "\n", "\n", "def _cy_ising_update(field, n, m, beta):\n", " total = 0\n", " N, M = field.shape\n", " for i in range(n-1, n+2):\n", " for j in range(m-1, m+2):\n", " if i == n and j == m:\n", " continue\n", " total += field[i%N, j%M]\n", " dE = 2 * field[n, m] * total\n", " if dE < 0:\n", " field[n, m] *= -1 # switch\n", " elif np.exp(-dE * beta) > np.random.rand():\n", " field[n, m] *= -1" ] }, { "cell_type": "code", "execution_count": 35, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "417 ms ± 21.2 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n", "393 ms ± 33.6 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n" ] } ], "source": [ "field = random_spin_field(200, 200)\n", "\n", "%timeit ising_step(field)\n", "\n", "%timeit cy_ising_step(field)" ] }, { "cell_type": "code", "execution_count": 36, "metadata": {}, "outputs": [], "source": [ "%%cython\n", "import numpy as np # gives access to python functions\n", "cimport numpy as np # gives you access to Numpy C API\n", "\n", "def cy_ising_step_fast(long[:, :]field, float beta=0.5):\n", " cdef int N = field.shape[0]\n", " cdef int M = field.shape[1]\n", " cdef int n_offset, m_offset, n, m\n", " \n", " for n_offset in range(2):\n", " for m_offset in range(2):\n", " for n in range(n_offset, N, 2):\n", " for m in range(m_offset, M, 2):\n", " _cy_ising_update_fast(field, n, m, beta)\n", " return field\n", "\n", "\n", "cdef _cy_ising_update_fast(long[:,:]field, int n, int m, float beta):\n", " cdef int total = 0\n", " cdef int N = field.shape[0]\n", " cdef int M = field.shape[1]\n", " cdef int i, j\n", " for i in range(n-1, n+2):\n", " for j in range(m-1, m+2):\n", " if i == n and j == m:\n", " continue\n", " total += field[i%N, j%M]\n", " cdef float dE = 2 * field[n, m] * total\n", " if dE < 0:\n", " field[n, m] *= -1 # switch\n", " elif np.exp(-dE * beta) > np.random.rand():\n", " field[n, m] *= -1" ] }, { "cell_type": "code", "execution_count": 37, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "414 ms ± 10.7 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n", "382 ms ± 15.1 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n", "64.5 ms ± 2.41 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)\n" ] } ], "source": [ "field = random_spin_field(200, 200)\n", "\n", "%timeit ising_step(field)\n", "\n", "%timeit cy_ising_step(field)\n", "\n", "%timeit cy_ising_step_fast(field)" ] }, { "cell_type": "code", "execution_count": 38, "metadata": {}, "outputs": [], "source": [ "%%cython \n", "import numpy as np # gives access to python functions\n", "cimport numpy as np # gives you access to Numpy C API\n", "\n", "from libc.math cimport exp\n", "from libc.stdlib cimport rand \n", "cdef extern from \"limits.h\":\n", " int RAND_MAX\n", "\n", "def cy_ising_step_faster(long[:, :]field, float beta=0.5):\n", " cdef int N = field.shape[0]\n", " cdef int M = field.shape[1]\n", " cdef int n_offset, m_offset, n, m\n", " \n", " for n_offset in range(2):\n", " for m_offset in range(2):\n", " for n in range(n_offset, N, 2):\n", " for m in range(m_offset, M, 2):\n", " _cy_ising_update_faster(field, n, m, beta)\n", " return field\n", "\n", "\n", "cdef _cy_ising_update_faster(long[:,:]field, int n, int m, float beta):\n", " cdef int total = 0\n", " cdef int N = field.shape[0]\n", " cdef int M = field.shape[1]\n", " cdef int i, j\n", " for i in range(n-1, n+2):\n", " for j in range(m-1, m+2):\n", " if i == n and j == m:\n", " continue\n", " total += field[i%N, j%M]\n", " cdef float dE = 2 * field[n, m] * total\n", " if dE < 0:\n", " field[n, m] *= -1 # switch\n", " elif exp(-dE * beta) > rand()/RAND_MAX:\n", " field[n, m] *= -1" ] }, { "cell_type": "code", "execution_count": 39, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "418 ms ± 13.4 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n", "386 ms ± 18.9 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n", "63.6 ms ± 1.61 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)\n", "5.54 ms ± 259 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)\n" ] } ], "source": [ "field = random_spin_field(200, 200)\n", "\n", "%timeit ising_step(field)\n", "\n", "%timeit cy_ising_step(field)\n", "\n", "%timeit cy_ising_step_fast(field)\n", "\n", "%timeit cy_ising_step_faster(field)" ] }, { "cell_type": "code", "execution_count": 40, "metadata": {}, "outputs": [], "source": [ "%%cython \n", "import numpy as np # gives access to python functions\n", "cimport numpy as np # gives you access to Numpy C API\n", "cimport cython\n", "\n", "from libc.math cimport exp\n", "from libc.stdlib cimport rand \n", "cdef extern from \"limits.h\":\n", " int RAND_MAX\n", "\n", "@cython.boundscheck(False) # dont check index boundaries\n", "@cython.wraparound(False) # can not use negative index\n", "def cy_ising_step_fastest(long[:, :]field, float beta=0.5):\n", " cdef int N = field.shape[0]\n", " cdef int M = field.shape[1]\n", " cdef int n_offset, m_offset, n, m\n", " \n", " for n_offset in range(2):\n", " for m_offset in range(2):\n", " for n in range(n_offset, N, 2):\n", " for m in range(m_offset, M, 2):\n", " _cy_ising_update_fastest(field, n, m, beta)\n", " return field\n", "\n", "\n", "cdef _cy_ising_update_fastest(long[:,:]field, int n, int m, float beta):\n", " cdef int total = 0\n", " cdef int N = field.shape[0]\n", " cdef int M = field.shape[1]\n", " cdef int i, j\n", " for i in range(n-1, n+2):\n", " for j in range(m-1, m+2):\n", " if i == n and j == m:\n", " continue\n", " total += field[i%N, j%M]\n", " cdef float dE = 2 * field[n, m] * total\n", " if dE < 0:\n", " field[n, m] *= -1 # switch\n", " elif exp(-dE * beta) > rand()/RAND_MAX:\n", " field[n, m] *= -1" ] }, { "cell_type": "code", "execution_count": 41, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "417 ms ± 22.6 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n", "377 ms ± 19.4 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n", "63.9 ms ± 1.95 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)\n", "5.48 ms ± 134 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)\n", "5.47 ms ± 120 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)\n" ] } ], "source": [ "field = random_spin_field(200, 200)\n", "\n", "%timeit ising_step(field)\n", "\n", "%timeit cy_ising_step(field)\n", "\n", "%timeit cy_ising_step_fast(field)\n", "\n", "%timeit cy_ising_step_faster(field)\n", "\n", "%timeit cy_ising_step_fastest(field)" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "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.7.10" }, "toc": { "base_numbering": 1, "nav_menu": {}, "number_sections": true, "sideBar": true, "skip_h1_title": false, "title_cell": "Table of Contents", "title_sidebar": "Contents", "toc_cell": false, "toc_position": {}, "toc_section_display": true, "toc_window_display": false } }, "nbformat": 4, "nbformat_minor": 4 }