From 528b427c432571b10b46b6039ff239fd78aaa413 Mon Sep 17 00:00:00 2001
From: Dimitri Lozeve I needed a small project to try APL while I was learning. Something array-based, obviously. Since I already implemented a Metropolis-Hastings simulation of the Ising model, which is based on a regular lattice, I decided to reimplement it in Dyalog APL. It is only a few lines long, but I will try to explain what it does step by step. The first function simply generates a random lattice filled by elements of . The first function simply generates a random lattice filled by elements of \(\{-1,+1\}\). Let’s deconstruct what is done here: Sample output: Sample output, for site in a random lattice: Sample output, for site \((3, 3)\) in a random \(5\times 5\) lattice: We can now bring everything together for display: Final output, with a random lattice, after 50000 update steps: Final output, with a \(80\times 80\) random lattice, after 50000 update steps: The Ising model is a model used to represent magnetic dipole moments in statistical physics. Physical details are on the Wikipedia page, but what is interesting is that it follows a complex probability distribution on a lattice, where each site can take the value +1 or -1. We have a lattice consisting of sites . For each site, there is a moment . is called the configuration of the lattice. The total energy of the configuration is given by the Hamiltonian where denotes neighbours, and is the interaction matrix. The configuration probability is given by: The configuration probability is given by: \[
\pi_\beta(\sigma) = \frac{e^{-\beta H(\sigma)}}{Z_\beta}
- where is the inverse temperature, and the normalisation constant. For our simulation, we will use a constant interaction term . If , the probability will be proportional to , otherwise it would be . Thus, adjacent spins will try to align themselves.The Ising model in APL
L←{(2×?⍵ ⍵⍴2)-3}
⍴
function: ⍵ ⍵⍴2
?
draws a random number between 1 and its argument. We give it our matrix to generate a random matrix of 1 and 2.L
.N
is the size of the lattice.xn
and yn
are respectively the vertical and lateral neighbours of the site. N|
takes the coordinates modulo N
(so the lattice is actually a torus). (Note: we used ⎕IO←0
to use 0-based array indexing.)+/
sums over all neighbours of the site, and then we multiply by the value of the site itself to get .+/
sums over all neighbours of the site, and then we multiply by the value of the site itself to get \(\Delta E\).
@@ -88,21 +90,21 @@
}
3 3ising.∆E ising.L 5
¯4
-
?
function.new
is the lattice but with the site flipped.*
function (exponential) and our previous ∆E
function.?0
returns a uniform random number in . Based on this value, we decide whether to update the lattice, and we return it.?
function.new
is the lattice but with the \((x,y)\) site flipped.*
function (exponential) and our previous ∆E
function.?0
returns a uniform random number in \([0,1)\). Based on this value, we decide whether to update the lattice, and we return it.Ising←{' ⌹'[1+1=({10 U ⍵}⍣⍵)L ⍺]}
-L ⍺
.⍣
function, which applies a function times.⍣
function, which applies a function \(n\) times. 80ising.Ising 50000
⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹ ⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹
⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹
diff --git a/_site/posts/ising-model.html b/_site/posts/ising-model.html
index e521987..b7b73a6 100644
--- a/_site/posts/ising-model.html
+++ b/_site/posts/ising-model.html
@@ -6,6 +6,8 @@
Mathematical definition
-
For our simulation, we will use a constant interaction term \(J > 0\). If \(\sigma_i = \sigma_j\), the probability will be proportional to \(\exp(\beta J)\), otherwise it would be \(\exp(\beta J)\). Thus, adjacent spins will try to align themselves.
The Ising model is generally simulated using Markov Chain Monte Carlo (MCMC), with the Metropolis-Hastings algorithm.
The algorithm starts from a random configuration and runs as follows:
The simulation is in Clojure, using the Quil library (a Processing library for Clojure) to display the state of the system.
@@ -57,7 +59,7 @@ H(\sigma) = -\sum_{i\sim j} J_{ij}\, \sigma_i\, \sigma_j, (:require [quil.core :as q] [quil.middleware :as m]))The application works with Quil’s functional mode, with each function taking a state and returning an updated state at each time step.
-The setup
function generates the initial state, with random initial spins. It also sets the frame rate. The matrix is a single vector in row-major mode. The state also holds relevant parameters for the simulation: , , and the iteration step.
The setup
function generates the initial state, with random initial spins. It also sets the frame rate. The matrix is a single vector in row-major mode. The state also holds relevant parameters for the simulation: \(\beta\), \(J\), and the iteration step.
(defn setup [size]
"Setup the display parameters and the initial state"
(q/frame-rate 300)
@@ -68,13 +70,13 @@ H(\sigma) = -\sum_{i\sim j} J_{ij}\, \sigma_i\, \sigma_j,
:beta 10
:intensity 10
:iteration 0}))
Given a site , we reverse its spin to generate a new configuration state.
+Given a site \(i\), we reverse its spin to generate a new configuration state.
(defn toggle-state [state i]
"Compute the new state when we toggle a cell's value"
(let [matrix (:matrix state)]
(assoc state :matrix (assoc matrix i (* -1 (matrix i))))))
In order to decide whether to accept this new state, we compute the difference in energy introduced by reversing site :
+In order to decide whether to accept this new state, we compute the difference in energy introduced by reversing site \(i\): \[ \Delta E = +J\sigma_i \sum_{j\sim i} \sigma_j. \]
The filter some?
is required to eliminate sites outside of the boundaries of the lattice.
(defn get-neighbours [state idx]
"Return the values of a cell's neighbours"
diff --git a/_site/posts/lsystems.html b/_site/posts/lsystems.html
index d7ebb83..e07aa59 100644
--- a/_site/posts/lsystems.html
+++ b/_site/posts/lsystems.html
@@ -6,6 +6,8 @@
Dimitri Lozeve - Generating and representing L-systems
+
+
@@ -40,29 +42,29 @@
Definition
An L-system is a set of rewriting rules generating sequences of symbols. Formally, an L-system is a triplet of:
-- an alphabet (an arbitrary set of symbols)
-- an axiom , which is a non-empty word of the alphabet ()
-- a set of rewriting rules (or productions) , each mapping a symbol to a word: . Symbols not present in are assumed to be mapped to themselves.
+- an alphabet \(V\) (an arbitrary set of symbols)
+- an axiom \(\omega\), which is a non-empty word of the alphabet (\(\omega \in V^+\))
+- a set of rewriting rules (or productions) \(P\), each mapping a symbol to a word: \(P \subset V \times V^*\). Symbols not present in \(P\) are assumed to be mapped to themselves.
-During an iteration, the algorithm takes each symbol in the current word and replaces it by the value in its rewriting rule. Not that the output of the rewriting rule can be absolutely anything in , including the empty word! (So yes, you can generate symbols just to delete them afterwards.)
+During an iteration, the algorithm takes each symbol in the current word and replaces it by the value in its rewriting rule. Not that the output of the rewriting rule can be absolutely anything in \(V^*\), including the empty word! (So yes, you can generate symbols just to delete them afterwards.)
At this point, an L-system is nothing more than a way to generate very long strings of characters. In order to get something useful out of this, we have to give them meaning.
Drawing instructions and representation
-Our objective is to draw the output of the L-system in order to visually inspect the output. The most common way is to interpret the output as a sequence of instruction for a LOGO-like drawing turtle. For instance, a simple alphabet consisting only in the symbols , , and could represent the instructions “move forward”, “turn right by 90°”, and “turn left by 90°” respectively.
+Our objective is to draw the output of the L-system in order to visually inspect the output. The most common way is to interpret the output as a sequence of instruction for a LOGO-like drawing turtle. For instance, a simple alphabet consisting only in the symbols \(F\), \(+\), and \(-\) could represent the instructions “move forward”, “turn right by 90°”, and “turn left by 90°” respectively.
Thus, we add new components to our definition of L-systems:
-- a set of instructions, . These are limited by the capabilities of our imagined turtle, so we can assume that they are the same for every L-system we will consider:
+
- a set of instructions, \(I\). These are limited by the capabilities of our imagined turtle, so we can assume that they are the same for every L-system we will consider:
Forward
makes the turtle draw a straight segment.
TurnLeft
and TurnRight
makes the turtle turn on itself by a given angle.
Push
and Pop
allow the turtle to store and retrieve its position on a stack. This will allow for branching in the turtle’s path.
Stay
, which orders the turtle to do nothing.
-- a distance , i.e. how long should each forward segment should be.
-- an angle used for rotation.
-- a set of representation rules . As before, they will match a symbol to an instruction. Symbols not matched by any rule will be associated to
Stay
.
+- a distance \(d \in \mathbb{R_+}\), i.e. how long should each forward segment should be.
+- an angle \(\theta\) used for rotation.
+- a set of representation rules \(R \subset V \times I\). As before, they will match a symbol to an instruction. Symbols not matched by any rule will be associated to
Stay
.
-Finally, our complete L-system, representable by a turtle with capabilities , can be defined as
+Finally, our complete L-system, representable by a turtle with capabilities \(I\), can be defined as \[ L = (V, \omega, P, d, \theta,
+R). \]
One could argue that the representation is not part of the L-system, and that the same L-system could be represented differently by changing the representation rules. However, in our setting, we won’t observe the L-system other than by displaying it, so we might as well consider that two systems differing only by their representation rules are different systems altogether.
Implementation details
The LSystem
data type
diff --git a/_site/skills.html b/_site/skills.html
index 30f1774..be07964 100644
--- a/_site/skills.html
+++ b/_site/skills.html
@@ -6,6 +6,8 @@
Dimitri Lozeve - Skills in Statistics, Data Science and Machine Learning
+
+
diff --git a/site.hs b/site.hs
index 1240cc4..84407fc 100644
--- a/site.hs
+++ b/site.hs
@@ -78,6 +78,6 @@ customPandocCompiler =
newExtensions = defaultExtensions `mappend` customExtensions
writerOptions = defaultHakyllWriterOptions
{ writerExtensions = newExtensions
- , writerHTMLMathMethod = MathML
+ , writerHTMLMathMethod = MathJax ""
}
in pandocCompilerWith defaultHakyllReaderOptions writerOptions
diff --git a/templates/default.html b/templates/default.html
index 2606a35..638121a 100644
--- a/templates/default.html
+++ b/templates/default.html
@@ -7,6 +7,7 @@
Dimitri Lozeve - $title$
+