using Plots

, LinearAlgebra

, ForwardDiff

, Roots

, PlutoUI

1.3 s

using BenchmarkTools

; TableOfContents(title="Table of Contents", indent=true)

266 ms

Plots.PlotlyBackendPlotlyBackend

plotly()

141 ms

Code Improvement

md"

# Code Improvement

6.4 ms

Tips for speeding it up

focus on kernel routines: e.g. mvup, mvdown, findneighbors
there are only finite number of results for findneighbors, so a look-up table for staple / Esite or dedicated fast function will be efficient
simplicity is king!
specify the types of inputs
in-line, avoid nesting, check closure, type stability
fixed / preallocated memory runs faster
- e.g. hard code the dim = 4 case would give us faster code
tools: @code_warntype, ProfileView, etc.

last but not least...

DO NOT OVERDO IT

md"

# Tips for speeding it up

* focus on kernel routines: e.g. mvup, mvdown, findneighbors

* there are only finite number of results for findneighbors, so a look-up table for staple / Esite or dedicated fast function will be efficient

* simplicity is king!

* specify the types of inputs

* in-line, avoid nesting, check closure, type stability

* fixed / preallocated memory runs faster

- e.g. hard code the dim = 4 case would give us faster code

* tools: @code_warntype, ProfileView, etc.

last but not least...

* DO NOT OVERDO IT

25.0 ms

old (generic function with 1 method)

function old(Tnow; Ns = 8, dim = 3)

# MC implementation of ND Ising Model

jcoup = 1.0

mu = 1.0

Bnow = 0.

beta = 1.0 / Tnow

latt = ones(Int, ntuple(i1 -> Ns, dim))

sites = CartesianIndices(latt)

# handle periodicity

function Ip(i1::Int)::Int

ret = (i1 + Ns) % Ns

# map 0 to Ns

ret += (ret == 0) * Ns

return ret

end

function Ipc(x::CartesianIndex)::CartesianIndex

tmp = ntuple(i1 -> Ip(x[i1]), dim)

return CartesianIndex(tmp)

end

function findneighbors(latt, site::CartesianIndex; bothway = true)

# find all the neighbors given a site (index)

14.9 ms

new (generic function with 3 methods)

function new(Tnow, Ns=8, dim=3)

# simplify mvup, mvdown

# dim must be known ntuple in the same scope

# remove garbage in findneighbors: e.g. bothway: i<j = (i != j) / 2

# tabulate / make fast function for staple to find Esite

# MC implementation of ND Ising Model

jcoup = 1.0

mu = 1.0

Bnow = 0.01

beta = 1.0 / Tnow

latt = ones(Int64, ntuple(i1 -> Ns, dim))

sites = CartesianIndices(latt)

function mvup(x::CartesianIndex, d::Int64)

local dim::Int64 = length(x)

i_tmp::Int64 = x[d] == Ns ? 1 : x[d] + 1

return CartesianIndex(ntuple(i1-> i1!=d ? x[i1] : i_tmp, dim))

11.6 ms

Enter cell code...

55.6 μs

BenchmarkTools.Trial

BenchmarkTools.Trial: 4 samples with 1 evaluation.
 Range (min … max):  1.654 s …  1.675 s  ┊ GC (min … max): 1.42% … 1.46%
 Time  (median):     1.666 s             ┊ GC (median):    1.50%
 Time  (mean ± σ):   1.666 s ± 9.644 ms  ┊ GC (mean ± σ):  1.58% ± 0.22%

  █                █                             █       █  
  █▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁█▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁█▁▁▁▁▁▁▁█ ▁
  1.65 s        Histogram: frequency by time        1.68 s <

 Memory estimate: 826.21 MiB, allocs estimate: 30466372.

BenchmarkTools.Trial: 777 samples with 1 evaluation.
 Range (min … max):  6.357 ms …  6.957 ms  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     6.408 ms              ┊ GC (median):    0.00%
 Time  (mean ± σ):   6.437 ms ± 72.294 μs  ┊ GC (mean ± σ):  0.00% ± 0.00%

   ▅█▆▄▂                 ▁▃▁                                  
  ▅█████▇▇▇▄▄▄▄▄▄▄▃▄▃▄▃▄▄███▇▇▆▄▄▃▂▃▂▃▃▃▂▁▂▃▂▃▂▁▂▁▁▂▁▂▁▂▂▁▁▃ ▃
  6.36 ms        Histogram: frequency by time        6.67 ms <

 Memory estimate: 26.81 KiB, allocs estimate: 1328.

let

p1 = @benchmark $old(2.0)

p2 = @benchmark $new(2.0)

[p1, p2]

end

31.5 s

Code Analysis

md"

# Code Analysis

136 μs

FASTEST ISING CODE ON EARTH

(constantly update :D) DO NOT GO CRAZY!!

md"

# FASTEST ISING CODE ON EARTH

(constantly update :D)

DO NOT GO CRAZY!!

164 μs

pokemon1 (generic function with 1 method)

function pokemon1(Tnow; Ns=12, dim=3)

# code improvement due to Biplab and Sasha!

# MC implementation of ND Ising Model

jcoup = 1.0

mu = 1.0

Bnow = 0.01

beta = 1.0 / Tnow

latt = ones(Int64, ntuple(i1 -> Ns, dim))

sites = CartesianIndices(latt)

function mvup(x::CartesianIndex, d::Int64)

local dim::Int64 = length(x)

i_tmp::Int64 = x[d] == Ns ? 1 : x[d] + 1

return CartesianIndex(ntuple(i1-> i1!=d ? x[i1] : i_tmp, dim))

end

function mvdown(x::CartesianIndex, d::Int64)

local dim::Int64 = length(x)

i_tmp::Int64 = x[d] == 1 ? Ns : x[d] - 1

return CartesianIndex(ntuple(i1-> i1!=d ? x[i1] : i_tmp, dim))

12.0 ms

pokemon2 (generic function with 1 method)

function pokemon2(Tnow; Ns=12, dim=3)

# MC implementation of ND Ising Model

# using iterators and tuples

jcoup = 1.0

mu = 1.0

Bnow = 0.01

beta = 1.0 / Tnow

latt = ones(Int64, ntuple(i1 -> Ns, dim))

# sites = Iterators.product(ntuple(i1->1:Ns, dim)...)

sites = Tuple.(CartesianIndices(latt))

function mvup(x, d::Int64)

local dim::Int64 = length(x)

i_tmp::Int64 = x[d] == Ns ? 1 : x[d] + 1

return ntuple(i1-> i1!=d ? x[i1] : i_tmp, dim)

end

function mvdown(x, d::Int64)

local dim::Int64 = length(x)

i_tmp::Int64 = x[d] == 1 ? Ns : x[d] - 1

return ntuple(i1-> i1!=d ? x[i1] : i_tmp, dim)

12.1 ms

sasha1 (generic function with 1 method)

function sasha1(T::Float64; B=0.01, Nd = 3, Ns = 12, Nstir = 150, Nm = 100, dNm = 10)

J = 1.0

μ = 1.0

β = 1 / T

# Spins

latt = ones(Int, ntuple(i -> Ns, Nd))

# Coordinate Grid

sites = CartesianIndices(latt)

# Neighbours (directional vectors)

neigsv = [CartesianIndex(ntuple(i -> Int(x == i), Nd)) for x in 1:Nd]

append!(neigsv, [CartesianIndex(ntuple(i -> -Int(x == i), Nd)) for x in 1:Nd])

# Periodic Index

PIndex(ind::CartesianIndex)::CartesianIndex = CartesianIndex(ntuple(i -> mod1(ind[i], Ns), Nd))

# Neighbours

neigs = [[PIndex(ind + n) for n in neigsv] for ind in sites]

# Energy of one spin

12.4 ms

sasha2 (generic function with 1 method)

function sasha2(T::Float64, ; B = 0.01, Nd = 3, Ns = 12, Nstir = 150, Nm = 100, dNm = 10)

J = 1.0

μ = 1.0

β = 1 / T

# Spins

latt = ones(Int, ntuple(i -> Ns, Nd))

# Coordinate Grid

sites = eachindex(latt)

# Strides and increments

str = strides(latt)

steps = ntuple(i -> Ns * str[i], Nd)

# Neighbours (with periodicity)

nh(x::Int64, s::Int64, a::Int64)::Int64 = x + s - a * (cld(x + s, a) - cld(x, a))

nl(x::Int64, s::Int64, a::Int64)::Int64 = x - s + a * (cld(x, a) - cld(x - s, a))

nf = [ntuple(i -> nh(site, str[i], steps[i]), Nd) for site in sites]

nb = [ntuple(i -> nl(site, str[i], steps[i]), Nd) for site in sites]

#nf = [[nh(site, str[i], steps[i]) for i in 1:Nd] for site in sites]

14.7 ms

biplab1 (generic function with 1 method)

function biplab1(T;B=0.01,Ns=12,dim=3)

## Define Parameters

J = 1.0

μ = 1.0

β = 1/T

## Create Lattice

latt = ones(Int64,ntuple(i->Ns,dim))

sites = CartesianIndices(latt)

## Unit Vectors

unitvec = [CartesianIndex(ntuple(i->(x==i),dim)) for x in 1:dim]

unitvec2 = .- unitvec

append!(unitvec,unitvec2)

## Periodic Boundary Condition

PBC(site) = CartesianIndex(mod1.(Tuple(site),Ns))

## List to save the neighbors data

neighbors = [[PBC(site+eáµ¢) for eáµ¢ in unitvec] for site in sites]

## List of possible energy for lookup

list_energy = [-J*neig - μ*B for neig in -2dim:2dim]

12.6 ms

biplab2 (generic function with 1 method)

function biplab2(T::Float64;B=0.01,Ns=12,dim=3)

## Define Parameters

J = 1.0

μ = 1.0

β = 1/T

## Create Lattice

latt = ones(Int64,ntuple(i->Ns,dim))

sites = CartesianIndices(latt)

## Unit Vectors

unitvec = [CartesianIndex(ntuple(i->(x==i),dim)) for x in 1:dim]

unitvec2 = .- unitvec

append!(unitvec,unitvec2)

## Periodic Boundary Condition

PBC(site::CartesianIndex) = CartesianIndex(mod1.(Tuple(site),Ns))

neighbors = [[PBC(site+dx) for dx in unitvec] for site in sites]

list_energy = [-J*neig - μ*B for neig in -2dim:2dim]

Esite(site::CartesianIndex) = list_energy[sum(latt[neighbor] for neighbor in neighbors[site]) + 2dim + 1]

14.8 ms

BenchmarkTools.Trial

BenchmarkTools.Trial: 127 samples with 1 evaluation.
 Range (min … max):  38.551 ms …  41.251 ms  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     39.276 ms               ┊ GC (median):    0.00%
 Time  (mean ± σ):   39.444 ms ± 667.939 μs  ┊ GC (mean ± σ):  0.00% ± 0.00%

      █▃█      ▂▂ ▆                                             
  ▅▄▅████▄▇▅▇▇███▇█▅█▇▄▅█▅▄▄▁▄▅▄▇▄▄█▄▄▅▇▁▇▁▁▅▁▄▁▁▄▄▅▄▄▁▄▄▄▁▁▁▄ ▄
  38.6 ms         Histogram: frequency by time         41.2 ms <

 Memory estimate: 1.19 MiB, allocs estimate: 27273.

BenchmarkTools.Trial: 181 samples with 1 evaluation.
 Range (min … max):  27.083 ms …  29.023 ms  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     27.649 ms               ┊ GC (median):    0.00%
 Time  (mean ± σ):   27.661 ms ± 376.113 μs  ┊ GC (mean ± σ):  0.00% ± 0.00%

    ▁           ▁ ▁ █▄▆▁▂ ▁                                     
  ▄▆██▅▆▄▇█▅█▇▆▄█▅███████▇██▆▄▃▆▁▃▁▁▁▃▃▁▁▁▅▁▁▃▁▄▁▃▃▃▄▃▃▁▁▁▁▁▁▃ ▃
  27.1 ms         Histogram: frequency by time         28.9 ms <

 Memory estimate: 257.58 KiB, allocs estimate: 5006.

BenchmarkTools.Trial: 168 samples with 1 evaluation.
 Range (min … max):  28.893 ms … 39.723 ms  ┊ GC (min … max): 0.00% … 25.38%
 Time  (median):     29.576 ms              ┊ GC (median):    0.00%
 Time  (mean ± σ):   29.770 ms ±  1.152 ms  ┊ GC (mean ± σ):  0.37% ±  2.38%

     ▁▇█▄█▃                                                    
  ▃▆████████▅▅▄▃▃▁▁▃▃▁▁▁▁▁▁▁▁▁▁▁▁▁▁▃▁▁▁▁▁▁▁▁▁▁▃▁▁▁▁▁▃▁▁▁▁▁▁▁▃ ▃
  28.9 ms         Histogram: frequency by time        35.1 ms <

 Memory estimate: 3.04 MiB, allocs estimate: 2158.

BenchmarkTools.Trial: 178 samples with 1 evaluation.
 Range (min … max):  27.148 ms …  36.560 ms  ┊ GC (min … max): 0.00% … 24.76%
 Time  (median):     27.985 ms               ┊ GC (median):    0.00%
 Time  (mean ± σ):   28.104 ms ± 834.281 μs  ┊ GC (mean ± σ):  0.35% ±  2.30%

        ▁▃▇█▂▁▅▃▁                                               
  ▃▁▃▁▄▇█████████▆▆▃▃▃▃▃▁▁▁▁▁▃▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▃ ▃
  27.1 ms         Histogram: frequency by time         31.9 ms <

 Memory estimate: 3.04 MiB, allocs estimate: 2158.

BenchmarkTools.Trial: 238 samples with 1 evaluation.
 Range (min … max):  20.838 ms …  21.462 ms  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     21.084 ms               ┊ GC (median):    0.00%
 Time  (mean ± σ):   21.086 ms ± 125.525 μs  ┊ GC (mean ± σ):  0.00% ± 0.00%

          ▃      █▄  ▂▅ ▃    ▅▇▃ ▃▂ ▂     ▂ ▂                   
  ▅▁▃▃▃▃█▆█▇█▆██▃██▆▅██▇██▇▆▅███▇██▇█▅▅████▆██▇▅▇▇█▃▃▅▃█▃▃▅▁▅▃ ▅
  20.8 ms         Histogram: frequency by time         21.4 ms <

 Memory estimate: 36.31 KiB, allocs estimate: 1328.

BenchmarkTools.Trial: 247 samples with 1 evaluation.
 Range (min … max):  19.938 ms …  20.688 ms  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     20.242 ms               ┊ GC (median):    0.00%
 Time  (mean ± σ):   20.241 ms ± 142.488 μs  ┊ GC (mean ± σ):  0.00% ± 0.00%

         ▂ ▂█▄ ▅▁   ▂▁▄▁▁▂ ▅ ▂▇▂ ▂ ▂▅ ▇▅▇▄▄▄▁▄                  
  ▃▁▁▃▁▅▅█▆███▆███▆█████████████▆█▅███████████▆▅▆▃▆▅▁▆▃▁▁▁▃▁▁▅ ▅
  19.9 ms         Histogram: frequency by time         20.6 ms <

 Memory estimate: 76.86 KiB, allocs estimate: 1331.

let

p1 = @benchmark sasha1(1.4)

p2 = @benchmark sasha2(1.4)

p3 = @benchmark biplab1(1.4)

p4 = @benchmark biplab2(1.4)

p5 = @benchmark pokemon1(1.4)

p6 = @benchmark pokemon2(1.4)

[p1, p2, p3, p4, p5, p6]

end

66.2 s

Tuple{Float64, Float64, Float64, Vararg{Float64}}

0.99478

0.00516408

-2.97916

0.995093

0.00480562

-2.99528

0.994757

0.00441708

-2.97405

0.0729466

let

Tnow = 5 * rand()

[sasha1(Tnow), biplab1(Tnow), pokemon1(Tnow)]

end

108 ms