03-example-cython-vs-cpp-vs-julia.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[![Binder](https://mybinder.org/badge_logo.svg)](https://mybinder.org/v2/gh/jolin-io/workshop-accelerate-Python-with-Julia/main?filepath=03-example-cython-vs-cpp-vs-julia.ipynb)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a href=\"https://www.jolin.io\" target=\"_blank\" rel=\"noreferrer noopener\">\n",
    "<img src=\"https://www.jolin.io/assets/Jolin/Jolin-Banner-Website-v1.3-darkmode.webp\">\n",
    "</a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# **Simulation example:** Python vs Cython vs C++ vs Julia\n",
    "\n",
    "The code is adapted from the [blog post by The Multi-Agent AI Guy](https://medium.com/agents-and-robots/the-bitter-truth-python-3-11-vs-cython-vs-c-performance-for-simulations-babc85cdfef5). Big thank you Multi-Agent AI Guy!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%html\n",
    "<iframe width=\"560\" height=\"315\" src=\"https://www.youtube.com/embed/hbFrPjeBpqg\" title=\"YouTube video player\" frameborder=\"0\" allow=\"accelerometer; autoplay; clipboard-write; encrypted-media; gyroscope; picture-in-picture; web-share\" allowfullscreen>\n",
    "</iframe>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Again let's activate julia"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from juliacall import Main as jl\n",
    "%load_ext juliacall.ipython\n",
    "# JuliaCall comes with its own Julia dependency file juliapkg.json\n",
    "# however for binder it is much simpler to just reuse binder's installation mechanism\n",
    "%julia Pkg.activate(Base.current_project())\n",
    "%julia using PythonCall\n",
    "%julia set_var(k, v) = @eval $(Symbol(k)) = $v"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import math\n",
    "import random\n",
    "\n",
    "from datetime import datetime\n",
    "random.seed(datetime.now().timestamp())\n",
    "\n",
    "WORLD_WIDTH = 2560\n",
    "WORLD_HEIGHT = 1440\n",
    "TIMESTEPS = 2000\n",
    "\n",
    "jl.set_var(\"WORLD_WIDTH\", WORLD_WIDTH)\n",
    "jl.set_var(\"WORLD_HEIGHT\", WORLD_HEIGHT)\n",
    "jl.set_var(\"TIMESTEPS\", TIMESTEPS)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Python implementation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Python definition of simulated Agent"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class PythonAgent():\n",
    "    def __init__(self, x=None, y=None, world_width=0, world_height=0):\n",
    "        super().__init__()\n",
    "\n",
    "        # default values\n",
    "        self.vmax = 2.0\n",
    "\n",
    "        # initial position\n",
    "        self.world_width = world_width\n",
    "        self.world_height = world_height\n",
    "        self.x = x if x else random.randint(0, self.world_width)\n",
    "        self.y = y if y else random.randint(0, self.world_height)\n",
    "\n",
    "        # initial velocity\n",
    "        self.dx = 0\n",
    "        self.dy = 0\n",
    "\n",
    "        # inital values\n",
    "        self.is_alive = True\n",
    "        self.target = None\n",
    "        self.age = 0\n",
    "        self.energy = 0\n",
    "\n",
    "    def update(self, food=()):\n",
    "        self.age = self.age + 1\n",
    "\n",
    "        # we can't move\n",
    "        if self.vmax == 0:\n",
    "            return\n",
    "\n",
    "        # target is dead, don't chase it further\n",
    "        if self.target and not self.target.is_alive:\n",
    "            self.target = None\n",
    "\n",
    "        # eat the target if close enough\n",
    "        if self.target:\n",
    "            squared_dist = (self.x - self.target.x) ** 2 + (self.y - self.target.y) ** 2\n",
    "            if squared_dist < 400:\n",
    "                self.target.is_alive = False\n",
    "                self.energy = self.energy + 1\n",
    "\n",
    "        # agent doesn't have a target, find a new one\n",
    "        if not self.target:\n",
    "            min_dist = 9999999\n",
    "            min_agent = None\n",
    "            for a in food:\n",
    "                if a is not self and a.is_alive:\n",
    "                    sq_dist = (self.x - a.x) ** 2 + (self.y - a.y) ** 2\n",
    "                    if sq_dist < min_dist:\n",
    "                        min_dist = sq_dist\n",
    "                        min_agent = a\n",
    "            if min_dist < 100000:\n",
    "                self.target = min_agent\n",
    "\n",
    "        # initalize 'forces' to zero\n",
    "        fx = 0\n",
    "        fy = 0\n",
    "\n",
    "        # move in the direction of the target, if any\n",
    "        if self.target:\n",
    "            fx += 0.1*(self.target.x - self.x)\n",
    "            fy += 0.1*(self.target.y - self.y)\n",
    "\n",
    "        # update our direction based on the 'force'\n",
    "        self.dx = self.dx + 0.05 * fx\n",
    "        self.dy = self.dy + 0.05 * fy\n",
    "\n",
    "        # slow down agent if it moves faster than it max velocity\n",
    "        velocity = math.sqrt(self.dx ** 2 + self.dy ** 2)\n",
    "        if velocity > self.vmax:\n",
    "            self.dx = (self.dx / velocity) * (self.vmax)\n",
    "            self.dy = (self.dy / velocity) * (self.vmax)\n",
    "\n",
    "        # update position based on delta x/y\n",
    "        self.x = self.x + self.dx\n",
    "        self.y = self.y + self.dy\n",
    "\n",
    "        # ensure it stays within the world boundaries\n",
    "        self.x = max(self.x, 0)\n",
    "        self.x = min(self.x, self.world_width)\n",
    "        self.y = max(self.y, 0)\n",
    "        self.y = min(self.y, self.world_height)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def main(Agent):\n",
    "\n",
    "    class Predator(Agent):\n",
    "        def __init__(self, **kwargs):\n",
    "            super().__init__(**kwargs)\n",
    "            self.vmax = 2.5\n",
    "\n",
    "    class Prey(Agent):\n",
    "        def __init__(self, **kwargs):\n",
    "            super().__init__(**kwargs)\n",
    "            self.vmax = 2.0\n",
    "\n",
    "    class Plant(Agent):\n",
    "        def __init__(self, **kwargs):\n",
    "            super().__init__(**kwargs)\n",
    "            self.vmax = 0.0\n",
    "\n",
    "            \n",
    "    # open the ouput file\n",
    "    # f = open('output.csv', 'w')\n",
    "    # print(0, ',', 'Title', ',', 'Predator Prey Relationship / Example 02 / Cython', file=f)\n",
    "\n",
    "    # create initial agents\n",
    "    kwargs = dict(\n",
    "        world_width=WORLD_WIDTH,\n",
    "        world_height=WORLD_HEIGHT,\n",
    "    )\n",
    "    preys = [Prey(**kwargs) for i in range(10)]\n",
    "    predators = [Predator(**kwargs) for i in range(10)]\n",
    "    plants = [Plant(**kwargs) for i in range(100)]\n",
    "\n",
    "    timestep = 0\n",
    "    while timestep < TIMESTEPS:\n",
    "        # update all agents\n",
    "        # no need to update the plants; they do not move\n",
    "        for a in preys:\n",
    "            a.update(plants)\n",
    "        \n",
    "        for a in predators:\n",
    "            a.update(preys)\n",
    "\n",
    "        # handle eaten and create new plant\n",
    "        plants = [p for p in plants if p.is_alive is True]\n",
    "        plants = plants + [Plant(**kwargs) for i in range(2)]\n",
    "\n",
    "        # handle eaten and create new preys\n",
    "        preys = [p for p in preys if p.is_alive is True]\n",
    "\n",
    "        for p in preys[:]:\n",
    "            if p.energy > 5:\n",
    "                p.energy = 0\n",
    "                preys.append(Prey(x = p.x + random.randint(-20, 20), y = p.y + random.randint(-20, 20), **kwargs))\n",
    "\n",
    "        # handle old and create new predators\n",
    "        predators = [p for p in predators if p.age < 2000]\n",
    "\n",
    "        for p in predators[:]:\n",
    "            if p.energy > 10:\n",
    "                p.energy = 0\n",
    "                predators.append(Predator(x = p.x + random.randint(-20, 20), y = p.y + random.randint(-20, 20), **kwargs))\n",
    "\n",
    "        # write data to output file\n",
    "        #[print(timestep, ',', 'Position', ',', 'Predator', ',', a.x, ',', a.y, file=f) for a in predators]\n",
    "        #[print(timestep, ',', 'Position',  ',', 'Prey', ',', a.x, ',', a.y, file=f) for a in preys]\n",
    "        #[print(timestep, ',', 'Position',  ',', 'Plant', ',', a.x, ',', a.y, file=f) for a in plants]\n",
    "\n",
    "        timestep = timestep + 1\n",
    "\n",
    "    print(len(predators), len(preys), len(plants))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Cython definition of Agent\n",
    "\n",
    "Cython has the unique advantage that you can easily reuse classes within Python\n",
    "\n",
    "The Cython package was already compiled and is available as `agent` module"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from agent import Agent as CythonAgent"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Benchmark Python against Cython class"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%timeit -r 3 main(PythonAgent)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%timeit main(CythonAgent)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Full language switch"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Cython\n",
    "\n",
    "The Cython module `agent` also includes a redefinition of the `main` routine, which gives a little extra speedup."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import agent\n",
    "%timeit agent.main(WORLD_WIDTH=WORLD_WIDTH, WORLD_HEIGHT=WORLD_HEIGHT, TIMESTEPS=TIMESTEPS)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## C++\n",
    "\n",
    "The blog post also came with a C++ implementation. Only a small bug was fixed, so that C++ and Python version have identical simulations.\n",
    "\n",
    "C++ needs be compiled first"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!g++ agent.cpp -o agent -std=c++11 -O3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "then we can time it from the commandline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!time ./agent $TIMESTEPS"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "and from python"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import subprocess\n",
    "%timeit subprocess.run([\"./agent\", str(TIMESTEPS)])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Julia\n",
    "\n",
    "Julia is not Object Oriented like Python, hence we need to restructure the class hierarchy.\n",
    "\n",
    "Luckily it is super simple, as Predator, Prey and Plant don't really make use of Object Orientation inheritance anyway."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "--------\n",
    "#### 💻 your space\n",
    "- 👉 DELETE ALL the following julia code\n",
    "- 👉 build the julia version yourself (🙂 you can always ask for help 🙂)\n",
    "- start untyped and finally add types (only needed for the `struct` fields)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%julia\n",
    "# your space"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%julia\n",
    "# Base.@kwdef is a nice little helper which enables default arguments and construction by keywords\n",
    "Base.@kwdef mutable struct Agent\n",
    "    vmax::Float64 = 2.5\n",
    "    world_width::Int = 10\n",
    "    world_height::Int = 10\n",
    "    x::Int = rand(0:world_width)\n",
    "    y::Int = rand(0:world_height)\n",
    "    \n",
    "    # initial velocity\n",
    "    dx::Float64 = 0.0\n",
    "    dy::Float64 = 0.0\n",
    "\n",
    "    # inital values\n",
    "    is_alive::Bool = true\n",
    "    target::Union{Nothing, Agent} = nothing\n",
    "    age::Int = 0\n",
    "    energy::Int = 0\n",
    "end\n",
    "\n",
    "Predator(; kwargs...) = Agent(; vmax = 2.5, kwargs...)\n",
    "Prey(; kwargs...) = Agent(; vmax = 2.0, kwargs...)\n",
    "Plant(; kwargs...) = Agent(; vmax = 0.0, kwargs...)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%julia\n",
    "function update!(self, food)\n",
    "    self.age = self.age + 1\n",
    "\n",
    "    # we can't move\n",
    "    if self.vmax == 0.0\n",
    "        return\n",
    "    end\n",
    "\n",
    "    # target is dead, don't chase it further\n",
    "    if self.target !== nothing && !self.target.is_alive\n",
    "        self.target = nothing\n",
    "    end\n",
    "    \n",
    "    # eat the target if close enough\n",
    "    if self.target !== nothing\n",
    "        squared_dist = (self.x - self.target.x) ^ 2 + (self.y - self.target.y) ^ 2\n",
    "        if squared_dist < 400\n",
    "            self.target.is_alive = false\n",
    "            self.energy = self.energy + 1\n",
    "        end\n",
    "    # agent doesn't have a target, find a new one\n",
    "    else\n",
    "        min_dist = 9999999\n",
    "        min_agent = nothing\n",
    "        for a in food\n",
    "            if a !== self && a.is_alive\n",
    "                sq_dist = (self.x - a.x) ^ 2 + (self.y - a.y) ^ 2\n",
    "                if sq_dist < min_dist\n",
    "                    min_dist = sq_dist\n",
    "                    min_agent = a\n",
    "                end\n",
    "            end\n",
    "        end\n",
    "        if min_dist < 100000\n",
    "            self.target = min_agent\n",
    "        end\n",
    "    end\n",
    "\n",
    "    # initalize 'forces' to zero\n",
    "    fx = 0.0\n",
    "    fy = 0.0\n",
    "\n",
    "    # move in the direction of the target, if any\n",
    "    if self.target !== nothing\n",
    "        fx += 0.1 * (self.target.x - self.x)\n",
    "        fy += 0.1 * (self.target.y - self.y)\n",
    "    end\n",
    "\n",
    "    # update our direction based on the 'force'\n",
    "    self.dx = self.dx + 0.05 * fx\n",
    "    self.dy = self.dy + 0.05 * fy\n",
    "\n",
    "    # slow down agent if it moves faster than it max velocity\n",
    "    velocity = sqrt(self.dx ^ 2 + self.dy ^ 2)\n",
    "    if velocity > self.vmax\n",
    "        self.dx = (self.dx / velocity) * (self.vmax)\n",
    "        self.dy = (self.dy / velocity) * (self.vmax)\n",
    "    end\n",
    "\n",
    "    # update position based on delta x/y\n",
    "    self.x = self.x + Int(round(self.dx))\n",
    "    self.y = self.y + Int(round(self.dy))\n",
    "\n",
    "    # ensure it stays within the world boundaries\n",
    "    self.x = max(self.x, 0)\n",
    "    self.x = min(self.x, self.world_width)\n",
    "    self.y = max(self.y, 0)\n",
    "    self.y = min(self.y, self.world_height)\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%julia\n",
    "function main(; show=true, world_width=WORLD_WIDTH, world_height=WORLD_HEIGHT, timesteps=TIMESTEPS)\n",
    "    kwargs = (; world_width, world_height)\n",
    "    preys = [Prey(; kwargs...) for i in 1:10]\n",
    "    predators = [Predator(; kwargs...) for i in 1:10]\n",
    "    plants = [Plant(; kwargs...) for i in 1:100]\n",
    "\n",
    "    timestep = 0\n",
    "    while timestep < timesteps\n",
    "        # update all agents\n",
    "        #[f.update([]) for f in plants]  # no need to update the plants; they do not move\n",
    "        for a in preys\n",
    "            update!(a, plants)\n",
    "        end\n",
    "        for a in predators\n",
    "            update!(a, preys)\n",
    "        end\n",
    "\n",
    "        # handle eaten and create new plant\n",
    "        # plants = [p for p in plants if p.is_alive]\n",
    "        filter!(p -> p.is_alive, plants)\n",
    "        append!(plants, [Plant(; kwargs...) for i in 1:2])\n",
    "\n",
    "        # handle eaten and create new preys\n",
    "        # preys = [p for p in preys if p.is_alive]\n",
    "        filter!(p -> p.is_alive, preys)\n",
    "        for p in preys\n",
    "            if p.energy > 5\n",
    "                p.energy = 0\n",
    "                push!(preys, Prey(; x = p.x + rand(-20:20), y = p.y + rand(-20:20), kwargs...))\n",
    "            end\n",
    "        end\n",
    "\n",
    "        # handle old and create new predators\n",
    "        # predators = [p for p in predators if p.age < 2000]\n",
    "        filter!(p -> p.age < 2000, predators)\n",
    "        for p in predators\n",
    "            if p.energy > 10\n",
    "                p.energy = 0\n",
    "                push!(predators, Predator(; x = p.x + rand(-20:20), y = p.y + rand(-20:20), kwargs...))\n",
    "            end\n",
    "        end\n",
    "\n",
    "        # write data to output file\n",
    "        #[print(timestep, ',', 'Position', ',', 'Predator', ',', a.x, ',', a.y, file=f) for a in predators]\n",
    "        #[print(timestep, ',', 'Position',  ',', 'Prey', ',', a.x, ',', a.y, file=f) for a in preys]\n",
    "        #[print(timestep, ',', 'Position',  ',', 'Plant', ',', a.x, ',', a.y, file=f) for a in plants]\n",
    "\n",
    "        timestep += 1\n",
    "    end\n",
    "    show && println(\"$(length(predators)), $(length(preys)), $(length(plants))\")\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "👉 delete up to this line\n",
    "\n",
    "----------"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Benchmark the performance of you julia code\n",
    "\n",
    "Try both `BenchmarkTools.@btime` and `BenchmarkTools.@benchmark`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%julia using BenchmarkTools\n",
    "%julia @btime main()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In order to compare C++ and Julia in more detail, it is informative to see the amount of garbage collection in Julia. (Because C++ does not have a garbage collector)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%julia @benchmark main(show=false)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Further resources\n",
    "\n",
    "- [learning Julia](https://julialang.org/learning/)\n",
    "- [juliacall / PythonCall.jl](https://cjdoris.github.io/PythonCall.jl/stable/) - my recommended interface to access Julia from Python\n",
    "- [pyjulia / PyCall.jl](https://github.com/JuliaPy/PyCall.jl) - older interface, but from julia it has the lovely `py\"\"` string to create python objects\n",
    "- [julia-numpy / TyPython.jl](https://github.com/Suzhou-Tongyuan/jnumpy) - implement a python package in Julia (think of this more like Cython pyx files)\n",
    "- [PackageCompiler.jl](https://julialang.github.io/PackageCompiler.jl/dev/) - create C libraries which you then can load from python\n",
    "\n",
    "ways to add Julia itself as a dependency:\n",
    "- `juliacall` uses [`juliapkg`](https://github.com/cjdoris/pyjuliapkg/) internally which also installs julia if it is not available\n",
    "- alternatively, the python package [`jill`](https://pypi.org/project/jill/) is widely used to install julia from python\n",
    "- the conda-forge [`julia`](https://github.com/conda-forge/julia-feedstock) package is well maintained\n",
    "- when compiling julia to a C library with PackageCompiler.jl, you do not need any extra Julia runnable - the C library is self-contained.\n",
    "\n",
    "# Thank you for participating ❤️\n",
    "\n",
    "You can always reach me at stephan.sahm@jolin.io or <a style=\"justify-content: center; padding: 7px; text-align: center; outline: none; text-decoration: none !important; color: #ffffff !important; width: 200px; height: 32px;border-radius: 16px; background-color: #0A66C2;\" href=\"https://www.linkedin.com/company/jolin-io/\" target=\"_blank\">Follow on LinkedIn</a>\n",
    "<!-- alternative follow https://www.linkedin.com/comm/mynetwork/discovery-see-all?usecase=PEOPLE_FOLLOWS&followMember=stephan-sahm-918656b7 -->"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a href=\"https://www.jolin.io\" target=\"_blank\" rel=\"noreferrer noopener\">\n",
    "<img src=\"https://www.jolin.io/assets/Jolin/Jolin-Banner-Website-v1.3-darkmode.webp\">\n",
    "</a>"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "workshop-accelerate-Python-with-Julia",
   "language": "python",
   "name": "workshop-accelerate-python-with-julia"
  },
  "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.10.9"
  },
  "vscode": {
   "interpreter": {
    "hash": "520e995520d0f28b9f1e7cacfd9ba1493aa60b57e5f0cc1543205df7dd9220a2"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}