Common Lisp vs. Python 3.11 vs. Cython vs. C++ Performance for Simulations

common-lisp-vs-python-3.11-vs-cython-vs-c++-performance-for-simulations.org

Common Lisp vs. Python 3.11 vs. Cython vs. C++ Performance for Simulations

Introduction

We compare the median run times of the agent-based simulation code discussed in the article titled The Bitter Truth: Python 3.11 vs Cython vs C++ Performance for Simulations published on Dec 19, 2022 against a corresponding Common Lisp implementation.

The Common Lisp implementation is (emphatically) written in non-idiomatic Lisp closely following the reference implementations, with some type declarations and compiled using SBCL.

The fastest implementation turns out to be Lisp, with a slight advantage over C++ (see table below).

Language	Factor by which Lisp is faster
Python 3.11	7.19
Cython	1.36
C++	1.06
Lisp	1

These results show that Common Lisp can be used for both prototyping and production implementations without re-implementing the hot spots in different languages (be it Cython or C++).

Methodology and results

We carried out one hundred repetitions for each method (multitime -n 100).

example_02.cpp (GCC 11.3 with -O2)

1: ./example_02
            Mean        Std.Dev.    Min         Median      Max
real        0.494       0.181       0.173       0.509       0.886       
user        0.490       0.182       0.173       0.509       0.886       
sys         0.003       0.003       0.000       0.004       0.012  

example_02.py (Python 3.11)

1: python example_02.py
            Mean        Std.Dev.    Min         Median      Max
real        3.951       2.535       1.430       3.453       16.334      
user        3.945       2.535       1.409       3.446       16.329      
sys         0.006       0.006       0.000       0.004       0.052    

example_02_cython.py (Python 3.11 + Cython 0.29.32)

1: python example_02_cython.py
            Mean        Std.Dev.    Min         Median      Max
real        1.020       0.796       0.579       0.655       2.913       
user        1.014       0.795       0.574       0.650       2.913       
sys         0.005       0.004       0.000       0.004       0.020

simulation.lisp (SBCL 2.1.11)

Executable built with:

(declaim (optimize speed))
(compile-file "simulation.lisp")
(load "simulation.fasl")
(save-lisp-and-die "simulation" :toplevel #'main :executable t)

The results:

1: ./simulation
            Mean        Std.Dev.    Min         Median      Max
real        0.528       0.239       0.422       0.480       2.188       
user        0.476       0.230       0.362       0.432       2.051       
sys         0.052       0.017       0.020       0.050       0.136

simulation.lisp

;;; Translation of example_02 in
;;;   https://github.com/multi-agent-ai/examples
;;; to non-idiomatic Common Lisp code.

(defconstant +world-width+ 2560d0)
(defconstant +world-height+ 1440d0)

(defclass agent ()
  ((vmax :initform 0.0d0 :type double-float)
   (x :initarg :x :type double-float)
   (y :initarg :y :type double-float)
   (dx :initform 0.0d0 :type double-float)
   (dy :initform 0.0d0 :type double-float)
   (alive :initform t :type boolean :reader alive-p)
   (target :initform nil)
   (age :initform 0 :type fixnum)
   (energy :initform 0 :type fixnum)))

(defmethod initialize-instance :after ((agent agent) &rest initargs)
  (declare (ignore initargs))
  (setf (slot-value agent 'x) (random +world-width+)
	(slot-value agent 'y) (random +world-height+)))

(defun update (agent food)
  (with-slots (vmax x y dx dy target age energy) agent
    (declare (type double-float vmax x y dx dy)
	     (type fixnum age energy))

    (incf age)
    
    ;; We can't move
    (when (zerop vmax)
      (return-from update nil))

    ;; Target is dead, don't chase it further
    (when (and target (not (alive-p target)))
      (setf target nil))

    ;; Eat the target if close enough
    (when target
      (with-slots ((target-x x) (target-y y)) target
	(declare (type double-float target-x target-y))
	(let ((squared-dist (+ (expt (- x target-x) 2) (expt (- y target-y) 2))))
	  (when (< squared-dist 400)
	    (setf (slot-value target 'alive) nil
		  energy (1+ energy))))))

    ;; Agent doesn't have a target, find a new one
    (when (null target)
      (loop :with min-dist :of-type double-float := 9999999d0
	    :with min-agent := nil
	    :for a :in food
	    :when (alive-p a) :do
	      (with-slots ((a-x x) (a-y y)) a
		(declare (type double-float a-x a-y))
		(let ((squared-dist (+ (expt (- x a-x) 2) (expt (- y a-y) 2))))
		  (when (< squared-dist min-dist)
		    (setf min-dist squared-dist
			  min-agent a))))
	    :finally (when (< min-dist 100000d0)
		       (setf target min-agent))))

    ;; Initialize forces to zero
    (let ((fx 0d0)
	  (fy 0d0))
      (declare (dynamic-extent fx fy))

      ;; Move in the direction of the target, if any
      (when target
	(with-slots ((target-x x) (target-y y)) target
	  (declare (type double-float target-x target-y))
	  (incf fx (* 1d-1 (- target-x x)))
	  (incf fy (* 1d-1 (- target-y y)))))

      ;; Update our direction based on the 'force'
      (incf dx (* 5d-2 fx))
      (incf dy (* 5d-2 fy))

      ;; Slow down agent if it moves faster than its max velocity
      (let ((velocity (sqrt (+ (expt dx 2) (expt dy 2)))))
	(unless (< velocity vmax)
	  (setf dx (* (/ dx velocity) vmax)
		dy (* (/ dy velocity) vmax))))

      (incf x dx)
      (incf y dy)

      (setf x (min (max x 0d0) +world-width+)
	    y (min (max y 0d0) +world-height+)))))

(defclass predator (agent) ())

(defmethod initialize-instance :after ((predator predator) &rest initargs)
  (declare (ignore initargs))
  (setf (slot-value predator 'vmax) 2.5d0))

(defclass prey (agent) ())

(defmethod initialize-instance :after ((prey prey) &rest initargs)
  (declare (ignore initargs))
  (setf (slot-value prey 'vmax) 2d0))

(defclass plant (agent) ())

(defmethod initialize-instance :after ((plant plant) &rest initargs)
  (declare (ignore initargs))
  (setf (slot-value plant 'vmax) 0d0))

(defun main ()
  (setf *random-state* (make-random-state t))

  (with-open-file (stream "output.csv" :direction :output :if-exists :supersede)
    (format stream "0, Title, Predator Prey Relationship / Example 02 / Lisp~%")

    (loop :with predators := (loop :repeat 10 :collect (make-instance 'predator))
	  :and preys := (loop :repeat 10 :collect (make-instance 'prey))
	  :and plants := (loop :repeat 100 :collect (make-instance 'plant))

	  :for timestep :below 10000 :do
	    ;; Update all agents
	    (loop :for p :in predators :do (update p preys))
	    (loop :for p :in preys :do (update p plants))

	    ;; Handle eaten and create new plants
	    (setf plants (remove-if-not #'alive-p plants))
	    (loop :repeat 2 :do (push (make-instance 'plant) plants))

	    ;; Handle eaten and create new preys
	    (setf preys (remove-if-not #'alive-p preys))
	    (loop :for p :in preys
		  :when (< 5 (slot-value p 'energy)) :do
		    (with-slots (energy x y) p
		      (declare (type double-float x y))
		      (setf energy 0)
		      (push (make-instance 'prey
					   :x (+ x (- (random 4d1) 2d1))
					   :y (+ y (- (random 4d1) 2d1)))
			    preys)))

	    ;; Handle old and create new predators
	    (setf predators (remove-if (lambda (a) (< 2000 (slot-value a 'age))) predators))
	    (loop :for p :in predators
		  :when (< 10 (slot-value p 'energy)) :do
		    (with-slots (energy x y) p
		      (declare (type double-float x y))
		      (setf energy 0)
		      (push (make-instance 'predator
					   :x (+ x (- (random 4d1) 2d1))
					   :y (+ y (- (random 4d1) 2d1)))
			    predators)))
	  :finally (format t "~D ~D ~D~%"
			   (length predators) (length preys) (length plants)))))