simulated_annealing-Solution.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Load data file\n",
    "\n",
    "### Target: Fix comments with `#!`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "using DelimitedFiles, Test, BenchmarkTools, Statistics\n",
    "\n",
    "\"\"\"General Annealing Problem\"\"\"\n",
    "abstract type AnnealingProblem end\n",
    "\n",
    "\"\"\"\n",
    "    SpinAnnealingProblem{T<:Real} <: AnnealingProblem\n",
    "\n",
    "Annealing problem defined by coupling matrix of spins.\n",
    "\"\"\"\n",
    "struct SpinAnnealingProblem{T<:Real} <: AnnealingProblem  # immutable, with type parameter T (a subtype of Real).\n",
    "    num_spin::Int\n",
    "    coupling::Matrix{T}\n",
    "    function SpinAnnealingProblem(coupling::Matrix{T}) where T\n",
    "        size(coupling, 1) == size(coupling, 2) || throw(DimensionMismatch(\"input must be square matrix.\"))\n",
    "        new{T}(size(coupling, 1), coupling)\n",
    "    end\n",
    "end\n",
    "\n",
    "\"\"\"\n",
    "    load_coupling(filename::String) -> SpinAnnealingProblem\n",
    "\n",
    "Load the data file into symmtric coupling matrix.\n",
    "\"\"\"\n",
    "function load_coupling(filename::String)\n",
    "    data = readdlm(filename)\n",
    "    is = Int.(view(data, :, 1)) .+ 1  #! @. means broadcast for the following functions, is here used correctly?\n",
    "    js = Int.(view(data, :, 2)) .+ 1\n",
    "    weights = data[:,3]\n",
    "    num_spin = max(maximum(is), maximum(js))\n",
    "    J = zeros(eltype(weights), num_spin, num_spin)\n",
    "    @inbounds for (i, j, weight) = zip(is, js, weights)\n",
    "        J[i,j] = weight/2\n",
    "        J[j,i] = weight/2\n",
    "    end\n",
    "    SpinAnnealingProblem(J)\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "@testset \"loading\" begin\n",
    "    sap = load_coupling(\"programs/example.txt\")\n",
    "    @test size(sap.coupling) == (300, 300)\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "abstract type AnnealingConfig end\n",
    "\n",
    "struct SpinConfig{Ts, Tf} <: AnnealingConfig\n",
    "    config::Vector{Ts}\n",
    "    field::Vector{Tf}\n",
    "end\n",
    "\n",
    "\"\"\"\n",
    "    random_config(prblm::AnnealingProblem) -> SpinConfig\n",
    "\n",
    "Random spin configuration.\n",
    "\"\"\"\n",
    "function random_config end   # where to put the docstring of a multiple-dispatch function is a problem. Using `abstract function` is proper.\n",
    "\n",
    "function random_config(prblm::SpinAnnealingProblem)\n",
    "    config = rand([-1,1], prblm.num_spin)\n",
    "    SpinConfig(config, prblm.coupling*config)\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "@testset \"random config\" begin\n",
    "    sap = load_coupling(\"programs/example.txt\")\n",
    "    initial_config = random_config(sap)\n",
    "    @test initial_config.config |> length == 300\n",
    "    @test eltype(initial_config.config) == Int\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Main Program for Annealing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "    anneal_singlerun!(config::AnnealingConfig, prblm, tempscales::Vector{Float64}, num_update_each_temp::Int)\n",
    "\n",
    "Perform Simulated Annealing using Metropolis updates for the single run.\n",
    "\n",
    "    * configuration that can be updated.\n",
    "    * prblm: problem with `get_cost`, `flip!` and `random_config` interfaces.\n",
    "    * tempscales: temperature scales, which should be a decreasing array.\n",
    "    * num_update_each_temp: the number of update in each temprature scale.\n",
    "\n",
    "Returns (minimum cost, optimal configuration).\n",
    "\"\"\"\n",
    "function anneal_singlerun!(config, prblm, tempscales::Vector{Float64}, num_update_each_temp::Int)\n",
    "    cost = get_cost(config, prblm)\n",
    "    \n",
    "    opt_config = config\n",
    "    opt_cost = cost\n",
    "    for beta = 1 ./ tempscales\n",
    "        @simd for m = 1:num_update_each_temp  # single instriuction multiple data, see julia performance tips.\n",
    "            proposal, ΔE = propose(config, prblm)\n",
    "            if exp(-beta*ΔE) > rand()  #accept\n",
    "                flip!(config, proposal, prblm)\n",
    "                cost += ΔE\n",
    "                if cost < opt_cost\n",
    "                    opt_cost = cost\n",
    "                    opt_config = config\n",
    "                end\n",
    "            end\n",
    "        end\n",
    "    end\n",
    "    opt_cost, opt_config\n",
    "end\n",
    " \n",
    "\"\"\"\n",
    "    anneal(nrun::Int, prblm, tempscales::Vector{Float64}, num_update_each_temp::Int)\n",
    "\n",
    "Perform Simulated Annealing with multiple runs.\n",
    "\"\"\"\n",
    "function anneal(nrun::Int, prblm, tempscales::Vector{Float64}, num_update_each_temp::Int)\n",
    "    local opt_config, opt_cost\n",
    "    for r = 1:nrun\n",
    "        initial_config = random_config(prblm)\n",
    "        cost, config = anneal_singlerun!(initial_config, prblm, tempscales, num_update_each_temp)\n",
    "        if r == 1 || cost < opt_cost\n",
    "            opt_cost = cost\n",
    "            opt_config = config\n",
    "        end\n",
    "    end\n",
    "    opt_cost, opt_config\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Annealing Problem Interfaces"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "    get_cost(config::AnnealingConfig, ap::AnnealingProblem) -> Real\n",
    "\n",
    "Get the cost of specific configuration.\n",
    "\"\"\"\n",
    "get_cost(config::SpinConfig, sap::SpinAnnealingProblem) = sum(config.config'*sap.coupling*config.config)\n",
    "\n",
    "\"\"\"\n",
    "    propose(config::AnnealingConfig, ap::AnnealingProblem) -> (Proposal, Real)\n",
    "\n",
    "Propose a change, as well as the energy change.\n",
    "\"\"\"\n",
    "@inline function propose(config::SpinConfig, ::SpinAnnealingProblem)  # ommit the name of argument, since not used.\n",
    "    ispin = rand(1:length(config.config))\n",
    "    @inbounds ΔE = -config.field[ispin] * config.config[ispin] * 4 # 2 for spin change, 2 for mutual energy.\n",
    "    ispin, ΔE\n",
    "end\n",
    "\n",
    "\"\"\"\n",
    "    flip!(config::AnnealingConfig, ispin::Proposal, ap::AnnealingProblem) -> SpinConfig\n",
    "\n",
    "Apply the change to the configuration.\n",
    "\"\"\"\n",
    "@inline function flip!(config::SpinConfig, ispin::Int, sap::SpinAnnealingProblem)\n",
    "    @inbounds config.config[ispin] = -config.config[ispin]  # @inbounds can remove boundary check, and improve performance\n",
    "    @simd for i=1:sap.num_spin\n",
    "        @inbounds config.field[i] += 2 * config.config[ispin] * sap.coupling[i,ispin]\n",
    "    end\n",
    "    config\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### **Challege!**\n",
    "Make your program correct and type is stable!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "using Random\n",
    "Random.seed!(2)\n",
    "const tempscales = 10 .- (1:64 .- 1) .* 0.15 |> collect\n",
    "const sap = load_coupling(\"programs/example.txt\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "@testset \"anneal\" begin\n",
    "    opt_cost, opt_config = anneal(30, sap, tempscales, 4000)\n",
    "    @test anneal(30, sap, tempscales, 4000)[1] == -3858\n",
    "    anneal(30, sap, tempscales, 4000)\n",
    "    res = median(@benchmark anneal(30, $sap, $tempscales, 4000))\n",
    "    @test res.time/1e9 < 2\n",
    "    @test res.allocs < 500\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "@benchmark anneal(30, $sap, $tempscales, 4000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Tips for optimization: Find the bottleneck of your program"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "using Profile\n",
    "Profile.clear()\n",
    "@profile anneal(100, sap, tempscales, 4000)\n",
    "Profile.print()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Calling a Fortran program\n",
    "* https://docs.julialang.org/en/v1/manual/calling-c-and-fortran-code/index.html\n",
    "* https://craftofcoding.wordpress.com/2017/02/26/calling-fortran-from-julia-i/\n",
    "* https://craftofcoding.wordpress.com/2017/03/01/calling-fortran-from-julia-ii/"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    ";cd programs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    ";gfortran -shared -fPIC problem.f90 fsa.f90 -o fsa.so"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    ";nm fsa.so"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "@benchmark ccall((:test_, \"fsa.so\"), Int32, ())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# For Python lovers?\n",
    "We can use [PyCall](https://github.com/JuliaPy/PyCall.jl) to call python programs!\n",
    "\n",
    "### **Challenge!**\n",
    "1. use Python package [viznet](https://github.com/GiggleLiu/viznet) and [matplotlib](https://matplotlib.org/) for visualization\n",
    "2. benchmark pure python version of simulated annealing, show the time"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# pip install viznet, matplotlib\n",
    "using PyCall\n",
    "\n",
    "@pyimport viznet\n",
    "@pyimport matplotlib.pyplot as plt\n",
    "brush = viznet.NodeBrush(\"nn.input\")\n",
    "brush >> (0, 0)\n",
    "plt.axis([-1, 1, -1, 1])\n",
    "plt.axis(\"equal\")\n",
    "plt.axis(\"off\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# we benchmark the `test_codec` function in file `testsa.py`\n",
    "pushfirst!(PyVector(pyimport(\"sys\")[\"path\"]), joinpath(@__DIR__, \"programs\"))  # add current folder into path\n",
    "@pyimport testsa\n",
    "@benchmark testsa.test_codec()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 1.8.1",
   "language": "julia",
   "name": "julia-1.8"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.8.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}