optimizing computer simulations of the Ising model - Gambit-list

25 Fév 2009

      Hi Gambitters,
I want to use Scheme to implement a computer simulation of the Ising model of
spins. Background: http://en.wikipedia.org/wiki/Ising_model
Basically all the algorithm does is allocate a vector of spins and flip them
randomly in such a way that a representative sampling at a specified temperature
is produced.
Efficiency of the algorithm is very important as it would allow me to study
larger samples.
Currently on my own system after compilation: with a sample size of 40x40 I can
do 1000 flips/site in a little over 3 seconds. Currently the simulation starts
from an all up state and displays some variables and the end configuration of
the sample. It looks like this:
$ time gsi ising -e "(main 1000)"
D: 2, N: 1600, L: 40, BETA*J: .3, RANDOM-STATE:
10624525222961816100342112271285465503733219408383542344109
E: 72, M: -40
+-++++---++---+-+-+-++++--+----+++-+-+--
+++-+---+-+------++---+-+--+--+++++++---
----+-+--+--+-+--+-+-+++-++---------++-+
---+---+--++-++++-++++-+--+++-++---+--+-
++++--++-+----++++----+-++---+----+-+++-
-----+++--+---+-+++++-++++-+---+++-+++-+
++-+--+-+---+--+++++-+-+---+++++-++---+-
-+---+--++--++--+---++++--+-----+-+-+--+
----+-+--+---++--+-++-+--+----+++---++--
--++-++---+-++++--+++-++++----+-+++++---
++-+-+-+-+-+++-+-++++++-+---+---++-+-+--
+--+---++-++++----+-++-+-+-++-++-------+
+--------+++-+-+++---+-+--+-+---++--+-++
+---++-+-+++++++--+-++--+-+++-++++---++-
+--+----+++--+--+-+---+++-+-+---+--+++++
++-+++-+++--+++--+-++++--++-+---+-+-++++
+-+---+++--+-+--+++-+--+---+-++++++-+++-
++---++---+-+--++-+++----+-+++--+++-+---
++-++-+--+---+++-----+-+----+--+-+-+-+-+
++--++--+--+-+++--++++--+++++--++-+--+-+
-+---++--+-++--++-+--++-+-++--+-++++-+-+
+----+-+--+++--+---+--+-+++--++--++-+++-
++++--+-++--++++----+-++-+--------+-++++
+-----+++++++-+-+--++-++-+--+++--++-++++
+--++-+-++-------+-+++-----++-+-++-+-++-
+--+++--+-+++--+---+---+-----++-+-++-+++
--+-+-+---+-+++----+-+--------+----+-+-+
-+--+--+-+++-+-----+-+----++++-++-++-++-
-+-++-+-+--+--+----+++--+-+-+-++--+++-+-
+-+++-+-++-+-+++++-+----+-+-+--+-++----+
+------+-+---+-++-+-++++-++-++-++---++-+
++----++-+-+--+-+-+++-+-+---+++-----++-+
++++-++++-++-+-+++--++++-++-+--+++++-++-
++-----+++++--++--+-+-+-++-+-+++------+-
++++-+++---+---+-+-++-+---+-+----++++--+
+---++--++---++-+----------+-----+---++-
---++-++--+--++--+--++-+----+--+-++----+
++--++--+--++++-+---++-++------++-++----
-++-++-++-++---+----+-++---++++++--++-++
+-----+++++--++-+-+---+-++---+-++--++-+-
real    0m3.507s
user    0m3.198s
sys     0m0.007s
to compile the program I use:
gsc -cc-options "-O2 -march=native" ising
I do not know enough about Gambit to know what's fast and what's not and I have
no experience with profiling, so I would appreciate some pointers on how to go
about optimizing the program.
As a side note I was also hoping to compare the performance of various Schemes
for this workload, but it seems the ones I tried so far (bigloo and guile) do
not support srfi-27 (Sources of Random Bits) much. (It worked in guile until
just now when I added a few bits to be able to change the pseudorandom number
generator seed.) If anyone knows other implementations where it would run that
would be cool.
Finally, the source is attached. Any comments appreciated.
Thanks in advance,
Marijn
-- 
Sarcasm puts the iron in irony, cynicism the steel.

Marijn Schouten (hkBst), Gentoo Lisp project, Gentoo ML
http://www.gentoo.org/proj/en/lisp/, #gentoo-{lisp,ml} on FreeNode

;;;calling

; $ gsc -cc-options "-O2 -march=native" ising 
; $ time gsi ising -e "(main 1000)"

;;;end calling

;;;compatibility

;bigloo compat
;(define random-integer random)

;bigloo has no random real generator!!! (define (random-real n) (random ?))

;guile compat
;guile does not support srfi-27 enough (define random-integer random) (define random-real random:uniform)

;;;end compat 

(define (SPIN+) +1)
(define (SPIN-) -1)

;;; coupling constant
(define J 1) ;ferromagnet
;(define J -1) ;antiferromagnet

;;; inverse temperature times coupling constant
(define (BETA*J) 0.3)
;(define (BETA*J) 0.4)
;(define (BETA*J) 0.5)

(define (display-spin spin)
  (cond 
   ((= spin (SPIN+)) (display "+"))
   ((= spin (SPIN-)) (display "-"))
   (#t display "?")))

;;; sample length in one dimension
(define L 40)

;;; sample dimension
(define D 2)

;;; total number of spins
;;; if this is a power of 2
;;; then modulo can be fast
(define N (expt L D))

;round N down to power of 2
;(set! N (expt 2 (inexact->exact (floor (/ (log N) (log 2))))))

(define GRID (make-vector N (SPIN+)))

(define (display-spatial-sample grid)
  (let loop ((site 0))
    (cond 
     ((< site N)
      (display-spin (vector-ref grid site))
      (if (= (- L 1) (modulo site L)) (newline))
      (loop (+ 1 site)))))
  (if (not (= (expt L D) N)) (newline))
  (newline))

(define (display-sample grid)
  (for-each display-spin (vector->list grid))
  (newline))

(define (dimension-loop site kons)
  (let loop ((d 0) (ret '()))
    (cond
     ((< d D)
      (loop (+ 1 d) (kons d site ret)))
     (#t ret))))

;;;helical boundary conditions
(define (helical-neighbors site)
  (dimension-loop site (lambda (d site ret)
    		 (cons (modulo (+ (expt L d) site) N)
    		       (cons (modulo (- site (expt L d)) N)
    			     ret)))))

(define (helical-plus-neighbors site)
  (dimension-loop site (lambda (d site ret)
    		 (cons (modulo (+ (expt L d) site) N)
    		       ret))))

;;; Metropolis => equal probabilities for each transition (single-spin-flip)
(define (random-site)
  (random-integer N))

(define (flip! grid site)
  (vector-set! grid site (* -1 (vector-ref grid site))))

(define (energy-jump/J grid site)
  (* 2 (vector-ref grid site)
     (apply + (map (lambda (n) (vector-ref grid n)) (helical-neighbors site)))))

;;; Metropolis => always accept energy-lowering transitions
(define (accept-flip grid site)
  (if (> 0 (energy-jump/J grid site))
      #t
      (< (random-real) (exp (* (BETA*J) (energy-jump/J grid site))))))

(define (single-spin-flip! grid)
  (let ((site (random-site)))
    (if (accept-flip grid site)
    (flip! grid site))))

(define (magnetization grid)
  (apply + (vector->list grid)))

(define (energy/J grid)
  (define (-energy/J-plus site)
    (* (vector-ref grid site)
       (apply + (map (lambda (n) (vector-ref grid n))
    	     (helical-plus-neighbors site)))))
  (let loop ((s 0) (E 0))
    (cond
     ((< s N)
      (loop (+ s 1) (- E (-energy/J-plus s))))
     (#t
      E))))

(define (energy grid)
  (* J (energy/J grid)))

(define (time-step)
  (let loop ((n 0))
    (cond
     ((< n N)
      (single-spin-flip! GRID)
      (loop (+ n 1))))))

(define (display-params . p)
  (let loop ((p p))
    (cond
     ((not (null? p))
      (display (car p))(display ": ")(display (cadr p))
      (cond
       ((not (null? (cddr p)))		; not last elt?
    (display ", ")
    (loop (cddr p)))))))
  (newline))

(define (display-config)
  (display-params "D" D "N" N "L" L "BETA*J" (BETA*J) "RANDOM-STATE" (random-source-state-ref default-random-source)))

(define TIMES #f)

(define (main args)
  (random-source-pseudo-randomize! default-random-source 0 0)
;  (write args)(newline)
  (set! TIMES args)
  (display-config)
;  (display-spatial-sample GRID)
  (let loop ((t 0))
    (cond
     ((< t TIMES)
      (time-step)
      (loop (+ t 1)))))
  (display-params "E" (energy GRID) "M" (magnetization GRID))
  (display-spatial-sample GRID)
)