dimitri/blog
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

Add posts

Browse source
This commit is contained in:
Dimitri Lozeve 2018-11-13 21:16:39 +01:00
parent 24e8c3b734
commit 3d56dad938
53 changed files with 1248 additions and 1 deletions
Download patch file Download diff file Expand all files Collapse all files

272
posts/ising-apl.org Normal file
View file

@ -0,0 +1,272 @@
---
title: Ising model simulation in APL
date: 2018-03-05
tags: ising simulation montecarlo apl
---
* The APL family of languages
** Why APL?
I recently got interested in [[https://en.wikipedia.org/wiki/APL_(programming_language)][APL]], an /array-based/ programming
language. In APL (and derivatives), we try to reason about programs as
series of transformations of multi-dimensional arrays. This is exactly
the kind of style I like in Haskell and other functional languages,
where I also try to use higher-order functions (map, fold, etc) on
lists or arrays. A developer only needs to understand these
abstractions once, instead of deconstructing each loop or each
recursive function encountered in a program.
APL also tries to be a really simple and /terse/ language. This
combined with strange Unicode characters for primitive functions and
operators, gives it a reputation of unreadability. However, there is
only a small number of functions to learn, and you get used really
quickly to read them and understand what they do. Some combinations
also occur so frequently that you can recognize them instantly (APL
programmers call them /idioms/).
** Implementations
APL is actually a family of languages. The classic APL, as created by
Ken Iverson, with strange symbols, has many implementations. I
initially tried [[https://www.gnu.org/software/apl/][GNU APL]], but due to the lack of documentation and
proper tooling, I went to [[https://www.dyalog.com/][Dyalog APL]] (which is proprietary, but free
for personal use). There are also APL derivatives, that often use
ASCII symbols: [[http://www.jsoftware.com/][J]] (free) and [[https://code.kx.com/q/][Q/kdb+]] (proprietary, but free for personal
use).
The advantage of Dyalog is that it comes with good tooling (which is
necessary for inserting all the symbols!), a large ecosystem, and
pretty good [[http://docs.dyalog.com/][documentation]]. If you want to start, look at [[http://www.dyalog.com/mastering-dyalog-apl.htm][/Mastering
Dyalog APL/]] by Bernard Legrand, freely available online.
* The Ising model in APL
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.html][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 $\{-1,+1\}$.
#+BEGIN_SRC apl
L←{(2×?⍵ ⍵⍴2)-3}
#+END_SRC
Let's deconstruct what is done here:
- ⍵ is the argument of our function.
- We generate a ⍵×⍵ matrix filled with 2, using the ~⍴~ 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.
- We multiply everything by 2 and subtract 3, so that the result is in $\{-1,+1\}$.
- Finally, we assign the result to the name ~L~.
Sample output:
#+BEGIN_SRC apl
ising.L 5
1 ¯1 1 ¯1 1
1 1 1 ¯1 ¯1
1 ¯1 ¯1 ¯1 ¯1
1 1 1 ¯1 ¯1
¯1 ¯1 1 1 1
#+END_SRC
Next, we compute the energy variation (for details on the Ising model,
see [[./ising-model.html][my previous post]]).
#+BEGIN_SRC apl
∆E←{
⎕IO←0
(x y)←⍺
N←⊃⍴⍵
xn←N|((x-1)y)((x+1)y)
yn←N|(x(y-1))(x(y+1))
⍵[x;y]×+/⍵[xn,yn]
}
#+END_SRC
- ⍺ is the left argument (coordinates of the site), ⍵ is the right argument (lattice).
- We extract the x and y coordinates of the site.
- ~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 $\Delta E$.
Sample output, for site $(3, 3)$ in a random $5\times 5$ lattice:
#+BEGIN_SRC apl
3 3ising.∆E ising.L 5
¯4
#+END_SRC
Then comes the actual Metropolis-Hastings part:
#+BEGIN_SRC apl
U←{
⎕IO←0
N←⊃⍴⍵
(x y)←?N N
new←⍵
new[x;y]×←(2×(?0)>*-⍺×x y ∆E ⍵)-1
new
}
#+END_SRC
- ⍺ is the $\beta$ parameter of the Ising model, ⍵ is the lattice.
- We draw a random site $(x,y)$ with the ~?~ function.
- ~new~ is the lattice but with the $(x,y)$ site flipped.
- We compute the probability $\alpha = \exp(-\beta\Delta E)$ using the ~*~ 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.
We can now bring everything together for display:
#+BEGIN_SRC apl
Ising←{' ⌹'[1+1=({10 U ⍵}⍣⍵)L ⍺]}
#+END_SRC
- We draw a random lattice of size ⍺ with ~L ⍺~.
- We apply to it our update function, with $\beta$=10, ⍵ times (using the ~⍣~ function, which applies a function $n$ times.
- Finally, we display -1 as a space and 1 as a domino ⌹.
Final output, with a $80\times 80$ random lattice, after 50000 update
steps:
#+BEGIN_SRC apl
80ising.Ising 50000
⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹ ⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹
⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹
⌹⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹ ⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹ ⌹
⌹⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹ ⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹ ⌹
⌹⌹⌹⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹ ⌹
⌹⌹⌹⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹ ⌹⌹⌹ ⌹⌹⌹⌹⌹⌹ ⌹
⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹ ⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹ ⌹⌹ ⌹⌹⌹⌹⌹ ⌹
⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹ ⌹ ⌹⌹⌹⌹ ⌹⌹⌹
⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹ ⌹⌹⌹ ⌹⌹⌹
⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹ ⌹⌹⌹ ⌹⌹⌹⌹⌹⌹
⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹
⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹
⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹
⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹ ⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹
⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹ ⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹
⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹
⌹⌹⌹⌹ ⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹ ⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹
⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹ ⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹
⌹ ⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹ ⌹⌹⌹⌹ ⌹⌹⌹ ⌹⌹⌹⌹⌹
⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹ ⌹⌹ ⌹
⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹⌹
⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹⌹⌹
⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹⌹ ⌹
⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹ ⌹
⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹ ⌹⌹⌹ ⌹⌹⌹⌹⌹⌹ ⌹⌹ ⌹⌹⌹
⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹ ⌹ ⌹⌹⌹ ⌹⌹⌹
⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹ ⌹⌹⌹
⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹
⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹
⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹
⌹⌹⌹⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹⌹⌹
⌹⌹⌹⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹⌹
⌹⌹⌹⌹⌹⌹⌹⌹⌹ ⌹ ⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹
⌹⌹⌹⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹
⌹⌹⌹⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹
⌹⌹⌹⌹⌹⌹⌹⌹⌹ ⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹
⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹ ⌹⌹ ⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹
⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹
⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹ ⌹⌹ ⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹
⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹
⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹
⌹⌹⌹⌹ ⌹ ⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹
⌹⌹⌹⌹⌹ ⌹ ⌹ ⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹
⌹⌹⌹⌹⌹ ⌹ ⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹
⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹
⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹
⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹
⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹⌹
⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹
⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹
⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹⌹⌹ ⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹
⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹⌹ ⌹ ⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹
⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹
⌹⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹ ⌹⌹ ⌹⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹
⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹
⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹
⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹⌹⌹
⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹
⌹⌹⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹
⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹⌹⌹
⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹ ⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹⌹⌹
⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹ ⌹⌹⌹ ⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹ ⌹⌹⌹
⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹ ⌹⌹⌹ ⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹ ⌹⌹
⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹ ⌹⌹⌹ ⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹ ⌹⌹
⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹ ⌹⌹⌹ ⌹⌹⌹ ⌹⌹ ⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹ ⌹⌹
⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹ ⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹
⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹ ⌹⌹ ⌹⌹⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹
⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹⌹
⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹
⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹ ⌹⌹⌹
⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹⌹ ⌹⌹⌹
⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹⌹ ⌹⌹⌹
⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹
⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹⌹ ⌹ ⌹⌹⌹⌹⌹ ⌹⌹⌹
⌹⌹⌹⌹⌹⌹⌹⌹⌹ ⌹ ⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹ ⌹⌹⌹ ⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹
⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹
⌹⌹⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹⌹ ⌹⌹
⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹⌹
⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹
⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹ ⌹⌹⌹⌹ ⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹⌹
#+END_SRC
Complete code, with the namespace:
#+BEGIN_SRC apl
:Namespace ising
L←{(2×?⍵ ⍵⍴2)-3}
∆E←{
⎕IO←0
(x y)←⍺
N←⊃⍴⍵
xn←N|((x-1)y)((x+1)y)
yn←N|(x(y-1))(x(y+1))
⍵[x;y]×+/⍵[xn,yn]
}
U←{
⎕IO←0
N←⊃⍴⍵
(x y)←?N N
new←⍵
new[x;y]×←(2×(?0)>*-⍺×x y ∆E ⍵)-1
new
}
Ising←{' ⌹'[1+1=({10 U ⍵}⍣⍵)L ⍺]}
:EndNamespace
#+END_SRC
* Conclusion
The algorithm is very fast (I think it can be optimized by the
interpreter because there is no branching), and is easy to reason
about. The whole program fits in a few lines, and you clearly see what
each function and each line does. It could probably be optimized
further (I don't know every APL function yet...), and also could
probably be golfed to a few lines (at the cost of readability?).
It took me some time to write this, but Dyalog's tools make it really
easy to insert symbols and to look up what they do. Next time, I will
look into some ASCII-based APL descendants. J seems to have a [[http://code.jsoftware.com/wiki/NuVoc][good
documentation]] and a tradition of /tacit definitions/, similar to the
point-free style in Haskell. Overall, J seems well-suited to modern
functional programming, while APL is still under the influence of its
early days when it was more procedural. Another interesting area is K,
Q, and their database engine kdb+, which seems to be extremely
performant and actually used in production.
Still, Unicode symbols make the code much more readable, mainly
because there is a one-to-one link between symbols and functions,
which cannot be maintained with only a few ASCII characters.

197
posts/ising-model.org Normal file
View file

@ -0,0 +1,197 @@
---
title: Ising model simulation
author: Dimitri Lozeve
date: 2018-02-05
tags: ising visualization simulation montecarlo
---
The [[https://en.wikipedia.org/wiki/Ising_model][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.
[[../images/ising.gif]]
* Mathematical definition
We have a lattice $\Lambda$ consisting of sites $k$. For each site,
there is a moment $\sigma_k \in \{ -1, +1 \}$. $\sigma =
(\sigma_k)_{k\in\Lambda}$ is called the /configuration/ of the
lattice.
The total energy of the configuration is given by the /Hamiltonian/
\[
H(\sigma) = -\sum_{i\sim j} J_{ij}\, \sigma_i\, \sigma_j,
\]
where $i\sim j$ denotes /neighbours/, and $J$ is the
/interaction matrix/.
The /configuration probability/ is given by:
\[
\pi_\beta(\sigma) = \frac{e^{-\beta H(\sigma)}}{Z_\beta}
\]
where $\beta = (k_B T)^{-1}$ is the inverse temperature,
and $Z_\beta$ the normalisation constant.
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.
* Simulation
The Ising model is generally simulated using Markov Chain Monte Carlo
(MCMC), with the
[[https://en.wikipedia.org/wiki/Metropolis%E2%80%93Hastings_algorithm][Metropolis-Hastings]]
algorithm.
The algorithm starts from a random configuration and runs as follows:
1. Select a site $i$ at random and reverse its spin: $\sigma'_i = -\sigma_i$
2. Compute the variation in energy (hamiltonian) $\Delta E = H(\sigma') - H(\sigma)$
3. If the energy is lower, accept the new configuration
4. Otherwise, draw a uniform random number $u \in ]0,1[$ and accept the new configuration if $u < \min(1, e^{-\beta \Delta E})$.
* Implementation
The simulation is in Clojure, using the [[http://quil.info/][Quil
library]] (a [[https://processing.org/][Processing]] library for
Clojure) to display the state of the system.
This post is "literate Clojure", and contains
[[https://github.com/dlozeve/ising-model/blob/master/src/ising_model/core.clj][=core.clj=]]. The
complete project can be found on
[[https://github.com/dlozeve/ising-model][GitHub]].
#+BEGIN_SRC clojure
(ns ising-model.core
(:require [quil.core :as q]
[quil.middleware :as m]))
#+END_SRC
The application works with Quil's
[[https://github.com/quil/quil/wiki/Functional-mode-(fun-mode)][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: $\beta$, $J$, and the iteration step.
#+BEGIN_SRC clojure
(defn setup [size]
"Setup the display parameters and the initial state"
(q/frame-rate 300)
(q/color-mode :hsb)
(let [matrix (vec (repeatedly (* size size) #(- (* 2 (rand-int 2)) 1)))]
{:grid-size size
:matrix matrix
:beta 10
:intensity 10
:iteration 0}))
#+END_SRC
Given a site $i$, we reverse its spin to generate a new configuration
state.
#+BEGIN_SRC clojure
(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))))))
#+END_SRC
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.
#+BEGIN_SRC clojure
(defn get-neighbours [state idx]
"Return the values of a cell's neighbours"
[(get (:matrix state) (- idx (:grid-size state)))
(get (:matrix state) (dec idx))
(get (:matrix state) (inc idx))
(get (:matrix state) (+ (:grid-size state) idx))])
(defn delta-e [state i]
"Compute the energy difference introduced by a particular cell"
(* (:intensity state) ((:matrix state) i)
(reduce + (filter some? (get-neighbours state i)))))
#+END_SRC
We also add a function to compute directly the hamiltonian for the
entire configuration state. We can use it later to log its values
across iterations.
#+BEGIN_SRC clojure
(defn hamiltonian [state]
"Compute the Hamiltonian of a configuration state"
(- (reduce + (for [i (range (count (:matrix state)))
j (filter some? (get-neighbours state i))]
(* (:intensity state) ((:matrix state) i) j)))))
#+END_SRC
Finally, we put everything together in the ~update-state~ function,
which will decide whether to accept or reject the new configuration.
#+BEGIN_SRC clojure
(defn update-state [state]
"Accept or reject a new state based on energy
difference (Metropolis-Hastings)"
(let [i (rand-int (count (:matrix state)))
new-state (toggle-state state i)
alpha (q/exp (- (* (:beta state) (delta-e state i))))]
;;(println (hamiltonian new-state))
(update (if (< (rand) alpha) new-state state)
:iteration inc)))
#+END_SRC
The last thing to do is to draw the new configuration:
#+BEGIN_SRC clojure
(defn draw-state [state]
"Draw a configuration state as a grid"
(q/background 255)
(let [cell-size (quot (q/width) (:grid-size state))]
(doseq [[i v] (map-indexed vector (:matrix state))]
(let [x (* cell-size (rem i (:grid-size state)))
y (* cell-size (quot i (:grid-size state)))]
(q/no-stroke)
(q/fill
(if (= 1 v) 0 255))
(q/rect x y cell-size cell-size))))
;;(when (zero? (mod (:iteration state) 50)) (q/save-frame "img/ising-######.jpg"))
)
#+END_SRC
And to reset the simulation when the user clicks anywhere on the screen:
#+BEGIN_SRC clojure
(defn mouse-clicked [state event]
"When the mouse is clicked, reset the configuration to a random one"
(setup 100))
#+END_SRC
#+BEGIN_SRC clojure
(q/defsketch ising-model
:title "Ising model"
:size [300 300]
:setup #(setup 100)
:update update-state
:draw draw-state
:mouse-clicked mouse-clicked
:features [:keep-on-top :no-bind-output]
:middleware [m/fun-mode])
#+END_SRC
* Conclusion
The Ising model is a really easy (and common) example use of MCMC and
Metropolis-Hastings. It allows to easily and intuitively understand
how the algorithm works, and to make nice visualizations!

221
posts/lsystems.org Normal file
View file

@ -0,0 +1,221 @@
---
title: Generating and representing L-systems
author: Dimitri Lozeve
date: 2018-01-18
tags: lsystems visualization algorithms haskell
---
L-systems are a formal way to make interesting visualisations. You can
use them to model a wide variety of objects: space-filling curves,
fractals, biological systems, tilings, etc.
See the Github repo: [[https://github.com/dlozeve/lsystems]]
* What is an L-system?
** A few examples to get started
[[../images/lsystems/dragon.png]]
[[../images/lsystems/gosper.png]]
[[../images/lsystems/plant.png]]
[[../images/lsystems/penroseP3.png]]
** Definition
An [[https://en.wikipedia.org/wiki/L-system][L-system]] is a set of
rewriting rules generating sequences of symbols. Formally, an L-system
is a triplet of:
+ 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 $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
$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/, $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/ $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 $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
The mathematical definition above translate almost immediately in a
Haskell data type:
#+BEGIN_SRC haskell
-- | L-system data type
data LSystem a = LSystem
{ name :: String
, alphabet :: [a] -- ^ variables and constants used by the system
, axiom :: [a] -- ^ initial state of the system
, rules :: [(a, [a])] -- ^ production rules defining how each
-- variable can be replaced by a sequence of
-- variables and constants
, angle :: Float -- ^ angle used for the representation
, distance :: Float -- ^ distance of each segment in the representation
, representation :: [(a, Instruction)] -- ^ representation rules
-- defining how each variable
-- and constant should be
-- represented
} deriving (Eq, Show, Generic)
#+END_SRC
Here, ~a~ is the type of the literal in the alphabet. For all
practical purposes, it will almost always be ~Char~.
~Instruction~ is just a sum type over all possible instructions listed
above.
** Iterating and representing
From here, generating L-systems and iterating is straightforward. We
iterate recursively by looking up each symbol in ~rules~ and replacing
it by its expansion. We then transform the result to a list of ~Instruction~.
** Drawing
The only remaining thing is to implement the virtual turtle which will
actually execute the instructions. It goes through the list of
instructions, building a sequence of points and maintaining an
internal state (position, angle, stack). The stack is used when ~Push~
and ~Pop~ operations are met. In this case, the turtle builds a
separate line starting from its current position.
The final output is a set of lines, each being a simple sequence of
points. All relevant data types are provided by the
[[https://hackage.haskell.org/package/gloss][Gloss]] library, along
with the function that can display the resulting ~Picture~.
* Common file format for L-systems
In order to define new L-systems quickly and easily, it is necessary
to encode them in some form. We chose to represent them as JSON
values.
Here is an example for the [[https://en.wikipedia.org/wiki/Gosper_curve][Gosper curve]]:
#+BEGIN_SRC json
{
"name": "gosper",
"alphabet": "AB+-",
"axiom": "A",
"rules": [
["A", "A-B--B+A++AA+B-"],
["B", "+A-BB--B-A++A+B"]
],
"angle": 60.0,
"distance": 10.0,
"representation": [
["A", "Forward"],
["B", "Forward"],
["+", "TurnRight"],
["-", "TurnLeft"]
]
}
#+END_SRC
Using this format, it is easy to define new L-systems (along with how
they should be represented). This is translated nearly automatically
to the ~LSystem~ data type using
[[https://hackage.haskell.org/package/aeson][Aeson]].
* Variations on L-systems
We can widen the possibilities of L-systems in various ways. L-systems
are in effect deterministic context-free grammars.
By allowing multiple rewriting rules for each symbol with
probabilities, we can extend the model to
[[https://en.wikipedia.org/wiki/Probabilistic_context-free_grammar][probabilistic
context-free grammars]].
We can also have replacement rules not for a single symbol, but for a
subsequence of them, thus effectively taking into account their
neighbours (context-sensitive grammars). This seems very close to 1D
cellular automata.
Finally, L-systems could also have a 3D representation (for instance
space-filling curves in 3 dimensions).
* Usage notes
1. Clone the repository: =git clone [[https://github.com/dlozeve/lsystems]]=
2. Build: =stack build=
3. Execute =stack exec lsystems-exe -- examples/penroseP3.json= to see the list of options
4. (Optional) Run tests and build documentation: =stack test --haddock=
Usage: =stack exec lsystems-exe -- --help=
#+BEGIN_SRC
lsystems -- Generate L-systems
Usage: lsystems-exe FILENAME [-n|--iterations N] [-c|--color R,G,B]
[-w|--white-background]
Generate and draw an L-system
Available options:
FILENAME JSON file specifying an L-system
-n,--iterations N Number of iterations (default: 5)
-c,--color R,G,B Foreground color RGBA
(0-255) (default: RGBA 1.0 1.0 1.0 1.0)
-w,--white-background Use a white background
-h,--help Show this help text
#+END_SRC
Apart from the selection of the input JSON file, you can adjust the
number of iterations and the colors.
=stack exec lsystems-exe -- examples/levyC.json -n 12 -c 0,255,255=
[[../images/lsystems/levyC.png]]
* References
1. Prusinkiewicz, Przemyslaw; Lindenmayer, Aristid (1990). /The Algorithmic Beauty of Plants./ Springer-Verlag. ISBN 978-0-387-97297-8. [[http://algorithmicbotany.org/papers/#abop]]
2. Weisstein, Eric W. "Lindenmayer System." From MathWorld--A Wolfram Web Resource. [[http://mathworld.wolfram.com/LindenmayerSystem.html]]
3. Corte, Leo. "L-systems and Penrose P3 in Inkscape." /The Brick in the Sky./ [[https://thebrickinthesky.wordpress.com/2013/03/17/l-systems-and-penrose-p3-in-inkscape/]]