Merge pull request #5 from skleinbo/makie

Makie

Project.toml

@@ -0,0 +1,18 @@
+name = "Ising"
+authors = ["Stephan Kleinboelting <[email protected]>"]
+version = "0.1.0"
+
+[deps]
+AbstractPlotting = "537997a7-5e4e-5d89-9595-2241ea00577e"
+CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
+DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
+DataFramesMeta = "1313f7d8-7da2-5740-9ea0-a2ca25f37964"
+Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
+GLM = "38e38edf-8417-5370-95a0-9cbb8c7f171a"
+GeometryTypes = "4d00f742-c7ba-57c2-abde-4428a4b178cb"
+LaTeXStrings = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a"
+Observables = "510215fc-4207-5dde-b226-833fc4488ee2"
+Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
+StatsPlots = "f3b207a7-027a-5e70-b257-86293d7955fd"

README.md

@@ -2,35 +2,72 @@
 
 This package implements Markov-chain-Monte-Carlo methods (MCMC) to study the equilibrium thermodynamic behavior of the Ising model on a square lattice. Currently it implements the classical Metropolis algorithm as well as Wolff's cluster algorithm.
 
-You'll find a source file with all the simulation routines and a notebook that ties them together into parametric studies, produces nice plots and so on.  
-Also included is a notebook for visualizing the dynamics See __"Live visualization"__ below.
+You'll find a source file with all the simulation routines and a notebook that ties them together into parametric studies, produces nice plots and so on.
+To watch the algorithms in action, a visualization app based on [Makie](http://makie.juliaplots.org/stable/index.html) is included. ![visualization](img/window.png)
+
 
 It is mainly intended for students of statistical physics learning about phase transitions and critical phenomena.
 
-For now the simulation is restricted to a square lattice with constant couplings and external fields. It should however be quite easy to extend the code to cover more general systems.
+For now the simulation is restricted to a square lattice with constant couplings and external fields. It should however be easy to extend the code to cover more general systems.
 
-__If you spot an error or are missing a feature, please feel free to open an issue or a pull request!__
+__If you spot an error or miss a feature, please feel free to open an issue or a pull request!__
 
-## Installation
+## Compatibility and Installation
+Julia >= 1.0 is required.
 
-### Compatibility
-Julia >= 0.6 is required for the simulation and data analysis portion of the notebook.
+The simulation routines assembled in `src/mcmc.jl` use base Julia only.
 
-At the moment the live-visualization requires Julia 0.6 and is not (yet) compatible with Julia 0.7.
+The analysis notebook utilizes various packages. Install them once with
+`]add Plots LaTeXStrings StatsPlots GLM DataFrames DataFramesMeta CSV Distributions`
 
-### Setup
+The visualization uses Makie and requires
+`]add Makie AbstractPlotting Colors Observables`
 
-Clone the files to your computer
-```bash
-  git clone https://github.com/skleinbo/Ising.git
-```
+Clone the repository to your computer
+`git clone https://github.com/skleinbo/Ising.git`
 
-The simulation routines assembled in `src/mcmc.jl` use base Julia only.
+## Usage
+Load the MCMC methods with `include(src/mcmc.jl)`.
 
-If you want to follow along the analysis in the accompanying notebook, a few packages need to be installed.  
-Run `src/build.jl` once to install any missing packages.
+Help is available
+```julia
+help?> run_metropolis
+search: run_metropolis _run_metropolis!
 
-## Usage
+  run_metropolis(L, beta, h;Tmax=1, sweep=0, sample_interval=1)
+
+  Sets up a random state and runs the Metropolis algorithm for a given set of parameters. Samples in defined intervals along the Markov-Chain.
+
+  An initial thermal sweep to go to equilibrium may be specified.
+
+  Returns an array of averaged observables: [E, E^2, m, m^2, m^4] with total energy E and magnetisation per spin m.
+
+  Arguments
+  ≡≡≡≡≡≡≡≡≡≡≡
+
+    •    L::Integer: Linear system size
+
+    •    beta::Float: Inverse temperature
+
+    •    h::Float: External field
+
+    •    Tmax::Integer: Number of steps
+
+    •    sweep::Integer: Length of the initial sweep
+
+    •    sample_interval::Integer: sample interval
+
+  Example
+  ≡≡≡≡≡≡≡≡≡
+
+  julia> run_metropolis(50, 0., 0.;Tmax=50*10^3*50^2,sample_interval=10*50^2,sweep=10^3*50^2)
+  5-element Array{Float64,1}:
+    -39.9664      # <E>
+   1753.5168      # <E^2>
+      0.010382079999999908 # <m>
+      0.00011572633600000021 # <m^2>
+      1.6924325969920386e-8 # <m^4>
+```
 
 Look around the source code to see which methods are available. Here's a minimal example to get you going:
 
@@ -51,11 +88,17 @@ julia> run_metropolis(L, 1/2.26, 0.; Tmax=25*10^3*L^2,sample_interval=10*L^2,swe
      1.0010325963705214e-10
 ```
 
-Admittedly the notebook is not that simple to understand if you hadn't had much exposure to Julia and DataFrames. Write your own routines!
+Admittedly the notebook may not be simple to understand if you hadn't had much exposure to Julia and DataFrames. And even if, it's not well documentes. __Write your own routines!__
 
 ### Live Visualization
-<center><img src="img/window.png" width=300></img></center>
-
-A few more packages are need to give us interactivity and render a scene. They will be installed automatically once you load `src/visualization.jl`.
-
-If you execute the cells at the end of the notebook in order, two things should happen. You should be presented with three sliders inside the notebook that allow you to control temperature, field strength, and simulation speed. A window should open in which the system is presented as a checker board.
+`include("Visual_Ising.jl")` opens a window with sliders and buttons to adjust
+parameters on the left, and a depiction of the current state on the right. Zoom in
+and out with the mouse wheel. "Speed"
+determines how many steps the algorithm takes per frame, i.e a value of `-2.0` means
+`1/100*L^2` updates per frame (default: 60fps; change with `BASE_FPS[]=$fps`).
+
+* _To change the system size (default 128x128):_ `g_L[]=_linear system size_`
+* Press <kbd>p</kbd> to run/pause the simulation
+* Reset the state by pressing <kbd>r</kbd>
+* Speed up/slow down with <kbd>=</kbd> / <kbd>-</kbd> (QWERTY layout)
+* Reset the view with <kbd>c</kbd>

Visual_Ising.jl

@@ -0,0 +1,121 @@
+using Makie
+using GeometryTypes, AbstractPlotting, Colors
+import Observables
+using Observables: AbstractObservable, on, off, async_latest
+
+include("src/mcmc.jl");
+include("src/visualization.jl")
+
+### GUI
+color1 = colorant"black"
+color2 = colorant"cornflowerblue"
+# s1, g_L = AbstractPlotting.textslider(16:16:512, "L")
+s2, g_mult = AbstractPlotting.textslider(-3.0:0.1:+2.0, start=-2.0, "speed [log]")
+
+s3, g_h = AbstractPlotting.textslider(-5f0:0.1f0:5f0, start=0f0,"field")
+s4, g_T = AbstractPlotting.textslider(0f0:0.005f0:5f0, "temperature")
+
+b1 = AbstractPlotting.button(Theme(raw = true, camera = campixel!),"Pause/Play" )
+b2 = AbstractPlotting.button(Theme(raw = true, camera = campixel!),"Metropolis" )
+
+b1_click = on(b1[end][:clicks]) do c
+    run_signal[] = !run_signal[]
+end
+
+b2_click = on(b2[end][:clicks]) do c
+    if algorithm[] == :metropolis
+        algorithm[] = :wolff
+        b2[end].input_args[end][] = "Wolff"
+    else
+        algorithm[] = :metropolis
+        b2[end].input_args[end][] = "Metropolis"
+    end
+end
+
+### Logic
+
+algorithm = Node(:metropolis)
+
+# Timey-Wimey
+frame_node = Node(0)
+run_signal = Node(false)
+BASE_FPS = Node(60)
+
+
+g_L = Node(128)
+# L = async_latest(g_L,1) #
+
+
+positions, square = get_primitives(g_L[])
+config0 = frustratedConfiguration(g_L[]);
+cluster = zeros(Bool, g_L[],g_L[]);
+
+color_node = Node([RGB(0.,0.,0.) for j in 1:g_L[]^2])
+
+render_map = lift(run_signal) do rs
+    global simtask = @async begin
+        if !rs
+            @info "Render task ended."
+            return nothing
+        end
+        while run_signal[]
+            t1 = time()
+            if algorithm[] == :metropolis
+                metropolis_sweep!(config0, round(Int,10^g_mult[] * g_L[]^2), 1/g_T.val,g_h.val)
+            elseif algorithm[] == :wolff
+                wolff_sweep!(round(Int,10^g_mult[] * g_L[]^2), config0, cluster, 1/g_T[], g_h[])
+            end
+            color_node[] = reshape(color_gen(config0, color1, color2),g_L[]^2)
+            t2 = time()
+            if (t2-t1) < 1/BASE_FPS[]
+                sleep(1/BASE_FPS[] - (t2-t1))
+            end
+            nothing
+        end
+    end
+end
+
+init_map = lift(g_L) do L
+    run_signal[] = false
+
+    global positions,square = get_primitives(L)
+    global config0 = frustratedConfiguration(L);
+    global cluster = zeros(Bool, L,L);
+
+    empty!(color_node.listeners)
+    color_node[] = reshape(color_gen(config0,color1, color2),L^2)
+
+    global state_plot = Makie.meshscatter(reshape(positions,L^2); color=color_node,
+        markersize=1,marker=GLNormalMesh(square), axis_type=axis2d!, camera=cam2d!,
+        raw=true, shading=false
+        )
+    # state_plot.attributes[:padding][] = [0f0, 0f0, 0f0]
+    # cam = cameracontrols(state_plot)
+    # cam.area[] = HyperRectangle(-1.05f0, -1.05f0, 2.1f0,2.1f0)
+    state_plot.theme.attributes[:backgroundcolor][] = :gray
+    update_cam!(state_plot, FRect(-1,-1,2,2))
+    state_plot[end][:light] = Vec{3,Float32}[[1.0, 1.0, 1.0], [0.1, 0.1, 0.1], [0.9, 0.9, 0.9], [0.0, 0.0, -20.0]]
+
+    global scene = vbox(hbox(s2,s3,s4,b2,b1),state_plot, parent=Scene(windowtitle="Ising"))
+    on(scene.events.keyboardbuttons) do buttons
+        if ispressed(scene, Keyboard.p)
+            run_signal[] = !run_signal[]
+        elseif ispressed(scene, Keyboard.c)
+            update_cam!(state_plot, FRect(-1,-1,2,2))
+        elseif ispressed(scene, Keyboard.r)
+            config0 = frustratedConfiguration(g_L[]);
+            cluster = zeros(Bool, g_L[],g_L[]);
+        elseif ispressed(scene, Keyboard.equal)
+            g_mult[] = min(g_mult[]+0.1, 2.0)
+        elseif ispressed(scene, Keyboard.minus)
+            g_mult[] = max(g_mult[]-0.1, -3.0)
+        end
+    end
+    on(scene.events.window_area) do wa
+        update_cam!(state_plot, FRect(-1,-1,2,2))
+    end
+    display(scene);
+    # push!(run_signal, true)
+    nothing
+end
+#

src/mcmc.jl

@@ -71,11 +71,11 @@ function metropolis_step!(state::Matrix{Int8},beta,h)
 end
 
 """
-    sweep!(state,n,beta,h)
+    metropolis_sweep!(state,n,beta,h)
 
 Perform n Metropolis steps.
 """
-function sweep!(state,n,beta,h)
+function metropolis_sweep!(state,n,beta,h)
     for _ in 1:n
         metropolis_step!(state,beta,h)
     end
@@ -90,7 +90,7 @@ Setup a configuration with `LxL` spins and perform an intial thermal sweep of
 function init(L,beta,h,sweep)
     state = randomConfiguration(L)
     # Initial sweep to get into the steady state
-    sweep!(state,sweep,beta,h)
+    metropolis_sweep!(state,sweep,beta,h)
     return state
 end
 
@@ -153,7 +153,7 @@ function _run_metropolis!(state::Matrix{Int8},beta,h;Tmax::Int=1,sample_interval
         while(t<Tmax)
             ## Take the defined no. of steps before
             ## recording a measurement
-            sweep!(state,sample_interval,beta,h)
+            metropolis_sweep!(state,sample_interval,beta,h)
 
             ## Record observables
             e = H(state,1.,h)
@@ -182,12 +182,12 @@ magnetetisation every `sample_interval` until `Tmax`.
 """
 function metropolis_timeseries(L::Int, beta, Tmax; sample_interval=L^2, sweep=1000)
     state = init(L,beta,0.,sweep*L^2)
-    ts = Vector{Float64}(div(Tmax,sample_interval))
+    ts = Vector{Float64}(undef, div(Tmax,sample_interval))
     t=0
     k=1
     mag0 = m(state)
     while(t<Tmax)
-        sweep!(state,sample_interval,beta,0.)
+        metropolis_sweep!(state,sample_interval,beta,0.)
         ts[k] = m(state)
         k+=1
         t+=sample_interval
@@ -212,7 +212,7 @@ function cluster!(cluster_state, state, i,j, p)
         cluster_state[i,j] = true
         for neighbor in [(i%L+1,j),(i,j%L+1),(i==1 ? L : i-1,j),(i,j==1 ? L : j-1)]
             if state[neighbor...] == s && !cluster_state[neighbor...] && rand()<p
-                cluster_state[neighbor...]=true
+                cluster_state[neighbor...] = true
                 cluster!(cluster_state, state, neighbor..., p)
             end
         end
@@ -224,7 +224,7 @@ end
 
 Build a cluster and flip it.
 """
-function wolff_step!(state, cluster_state, beta)
+function wolff_step!(state, cluster_state, beta, h)
     @inbounds begin
         L = size(state,1)
         i,j = rand(1:L,2)
@@ -235,6 +235,12 @@ function wolff_step!(state, cluster_state, beta)
     return nothing
 end
 
+function wolff_sweep!(n, args...)
+    for _ in 1:n
+        wolff_step!(args...)
+    end
+end
+
 """
     _run_wolff!(state,cluster,beta,h;Tmax=1,sample_interval=1)
 
@@ -264,7 +270,7 @@ function _run_wolff!(state::Matrix{Int8},cluster,beta,h;Tmax::Int=1,sample_inter
             ## Do the defined no. of steps before
             ## taking a measurement
             for _ in 1:sample_interval
-                wolff_step!(state, cluster, beta)
+                wolff_step!(state, cluster, beta, h)
             end
             ## Record observables
             e = H(state,1.,0.)
@@ -296,7 +302,7 @@ function run_wolff(L::Int, beta,h;Tmax::Int=1,sweep::Int=0,sample_interval::Int=
         cluster = zeros(Bool,L,L)
         # Initial sweep to get into the steady state
         for _ in 1:sweep
-            wolff_step!(state, cluster, beta)
+            wolff_step!(state, cluster, beta, h)
         end
         return state,cluster
     end