Monte Carlo Pi (Clojure)

monte_carlo_pi_fast.md

Monte Carlo estimation for the value of π

Iteration 1

(defn pair
  []
  [(rand) (rand)])

(defn pairs
  [n]
  (let [pairs (take n (repeatedly pair))]
    pairs))

(defn within
  [pair]
  (let [fst (first pair)
        snd (second pair)
         r (+ (* fst fst) (* snd snd)) ]
    (if (<= r 1.0)
      :in
      :out)))

(defn pi
  [size]
  (double (* 4 (/ (count (filter #(= :in %) (map within (pairs size)))) size))))

Iteration 2

This is a really bad idea. Not random and much worse than the random versions.

(defn range-float-0-1
  "Return a range of floats between 0 and 1 with 1/resolution steps
  (range-float 10)
  (0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0)
  This is inclusive on both ends"
  [resolution] 
  (map #(float (/ % resolution)) (range 0 (inc resolution))))

(defn coordinates
  "Returns coordinates on a 2D space between [0 0] and [1 1] using 1/resolution as step"
  [resolution]
  (let [r (range-float-0-1 resolution)]
    (for [x r y r] 
      [x y])))

(defn within?
  "Is the dot within a 1 unit distance from the origin"
  [v]
  (if (<= 
        (+ 
          (Math/pow (first v) 2)
          (Math/pow (second v) 2))
        1) true false))
  
(defn calc-pi 
  [rez]
  (let [
         coordz (coordinates rez)
         freq (frequencies (pmap within? coordz))
         wthn (get freq 'true)
         total (count coordz)
         ]
  (float (* 4 (/ wthn total)))))

Iteration 3

(def pi-seq 
  (pmap first 
    (iterate 
      (fn [[pi within total]] 
        [ (double (* (/ within total) 4)) 
          (if (<= (+ (Math/pow (rand) 2) (Math/pow (rand) 2)) 1) (inc within) within)
          (inc total) ] ) 
      [0.0 0N 1N])))

(nth pi-seq 1000000)
3.140492

(def pi-seq-with-error 
  (pmap 
    (fn [[pi _ _ err]] [pi err])
      (iterate 
        (fn [[pi within total error]] 
          [ (double (* (/ within total) 4)) 
            (if (<= (+ (Math/pow (rand) 2) (Math/pow (rand) 2)) 1) (inc within) within)
            (inc total) 
            (* (Math/abs (/ (- pi 3.141592653589793238462643383280) 3.141592653589793238462643383280)) 100) ] ) 
        [0.0 0N 1N 0.0])))

(doseq [i (range 10)] (println (map #(format "%8f" %) (nth pi-seq-with-error (* 1000000 i)))))
(0.000000 -1.000000)
(3.138925 -0.000849)
(3.139586 -0.000638)
(3.140864 -0.000232)
(3.140490 -0.000351)
(3.140713 -0.000280)
(3.140729 -0.000275)
(3.140905 -0.000219)
(3.141455 -0.000044)
(3.141520 -0.000023)

Benchmark

=> (criterium/bench (nth pi-seq 10000))
Evaluation count : 271680 in 60 samples of 4528 calls.
             Execution time mean : 217.588823 µs
    Execution time std-deviation : 5.384881 µs
   Execution time lower quantile : 213.853079 µs ( 2.5%)
   Execution time upper quantile : 228.646360 µs (97.5%)
                   Overhead used : 8.250169 ns

Found 5 outliers in 60 samples (8.3333 %)
        low-severe       3 (5.0000 %)
        low-mild         2 (3.3333 %)
 Variance from outliers : 12.5774 % Variance is moderately inflated by outliers

Iteration 4

This is clearly the winner

Defs

(defn pi-error
  [pi-estimated]
  (let [pi-known 3.141592653589793238462643383280]
    (* (Math/abs (/ (- pi-estimated pi-known) pi-known)) 100)))

(defn experiment 
  [] 
  (if (<= (+ (Math/pow (rand) 2) (Math/pow (rand) 2)) 1) 1 0))

(defn pi-estimate 
  [n] 
  (* 4 (float (/ (reduce + (take n (repeatedly experiment))) n))))  

(doseq [i (range 1000000 10000000 1000000)] 
   (let [pi-estimate (pi-estimate i) pi-err (pi-error pi-estimate) ] 
      (println (format "%d  %10f  %10f" i pi-estimate pi-err))))
1000000    3.142592    0.031809
2000000    3.143582    0.063326
3000000    3.142743    0.036605
4000000    3.143556    0.062499
5000000    3.141274    0.010129
6000000    3.141701    0.003456
7000000    3.140828    0.024343
8000000    3.143265    0.053248
9000000    3.141982    0.012403      

Benchmarks

=> (criterium/bench (pi-estimate 10000))
Evaluation count : 40980 in 60 samples of 683 calls.
             Execution time mean : 1.469525 ms
    Execution time std-deviation : 14.329410 µs
   Execution time lower quantile : 1.456673 ms ( 2.5%)
   Execution time upper quantile : 1.507841 ms (97.5%)
                   Overhead used : 8.250169 ns

Found 4 outliers in 60 samples (6.6667 %)
        low-severe       2 (3.3333 %)
        low-mild         2 (3.3333 %)
 Variance from outliers : 1.6389 % Variance is slightly inflated by outliers

Iteration 5

Make is a little faster (not actually faster).

Defs

(defn pi-error
  ^Double [^Double pi-estimated]
  (let [pi-known 3.141592653589793238462643383280]
    (* (Math/abs (/ (- pi-estimated pi-known) pi-known)) 100)))

(defn experiment 
  ^Long [] 
  (if (<= (+ (Math/pow (rand) 2) (Math/pow (rand) 2)) 1) 1 0))

(defn pi-estimate 
  ^Double [n] 
  (* 4 (double (/ (reduce + (take n (repeatedly experiment))) n))))

Benchmarks

=> (criterium/bench (pi-estimate 10000))
Evaluation count : 39900 in 60 samples of 665 calls.
             Execution time mean : 1.513959 ms
    Execution time std-deviation : 7.673476 µs
   Execution time lower quantile : 1.504963 ms ( 2.5%)
   Execution time upper quantile : 1.530864 ms (97.5%)
                   Overhead used : 8.250169 ns

Found 2 outliers in 60 samples (3.3333 %)
        low-severe       1 (1.6667 %)
        low-mild         1 (1.6667 %)
 Variance from outliers : 1.6389 % Variance is slightly inflated by outliers