A Monte Carlo simulation to estimate Pi with Clojure and Incanter.

core.clj

(ns montecarlopi.core
  (require [incanter.core   :refer [sqrt sq view]]
           [incanter.charts :refer [function-plot
                                    add-points
                                    add-text
                                    set-x-label
                                    set-y-label
                                    set-y-range
                                    xy-plot]]))

(defn in-circle? [point]
  (let [[x y] point]
    (if (<= (+ (sq x)
               (sq y))
            1.0)
      true
      false)))

(defn outside-circle? [point]
  (not (in-circle? point)))

(defn generate-random-point []
  [(rand) (rand)])

(defn generate-random-points [n]
  (take n (repeatedly generate-random-point)))

(defn points-in-circle [points]
  (filter in-circle? points))

(defn sort-points [points]
  {:inside  (filter in-circle? points)
   :outside (filter outside-circle? points)})

(defn point-ratio [points]
  (float (/ (count (points-in-circle points))
            (count points))))

(defn estimate-pi [points]
  (* 4 (point-ratio points)))

(defn circle [x] (sqrt (- 1.0 (sq x))))

(defn draw-circle []
  (function-plot circle 0 1)) 

(defn plot-points [chart points label]
  (let [xs (map first  points)
        ys (map second points)]
    (add-points chart xs ys :series-label label)))

(defn generate-and-sort-points [n]
  (sort-points (generate-random-points n)))

(defn label [point]
  (if (in-circle? point)
    "inside"
    "outside"))

(defn make-plot [n]
  (let [points (generate-random-points n)
        sorted (sort-points points)]
    (doto (draw-circle)
      (set-y-range -0.25 1)
      (plot-points (:inside  sorted) "inside")
      (plot-points (:outside sorted) "outside")
      (set-x-label "")
      (set-y-label "")
      (add-text 0.10 -0.10 (str "Total: "
                                (count points)))
      (add-text 0.10 -0.05 (str "Inside: "
                                (count (:inside sorted))))
      (add-text 0.10 -0.15 (str "Ratio: "
                                (point-ratio points)))
      (add-text 0.10 -0.20 (str "Pi: "
                                (format "%4f" (estimate-pi points))))
      (view :width 500 :height 600))))

(defn -main [& args]
 (make-plot 100000))

core_spec.clj

(ns montecarlopi.core-spec
  (:require [speclj.core :refer :all]
            [montecarlopi.core :refer :all]))


(def inside  [[0.0 0.0]
              [0.2 0.3]
              [0.1 0.8]
              [0.0 1.0]
              [1.0 0.0]]) 

(def outside [[1.0 0.2]
              [0.8 0.8]
              [0.5 0.9]
              [0.8 0.7]
              [0.9 0.9]])

(def points (concat inside outside))

(describe "in-circle?"
  (it "should return true for points inside the circle."
      (should (in-circle? [0 0.1])))

  (it "should return false for points outside the circle."
      (should-not (in-circle? [0.8 0.9]))))

(describe "generate-random-point"
  (it "should return a point with x and y values between 0 and 1."
      (let [[x y] (generate-random-point)]
        (should (<= x 1))
        (should (>= x 0))

        (should (<= y 1))
        (should (>= y 0)))))

(describe "generate-random-points"
  (it "should return the specified number of random points."
    (should= 100 (count (generate-random-points 100)))))

(describe "points-in-circle"
  (it "should return only points inside the circle."
    (should= 5 (count (points-in-circle points)))
    (should= inside (points-in-circle points))
    (should-not= outside (points-in-circle points))))

(describe "point-ratio"
  (it "should return the correct ratio of inside to outside points"
      (should= 0.5 (point-ratio points))))

(describe "outside-circle?"
  (it "should return true for points outside the circle."
      (should (outside-circle? [0.9 0.7]))))

(describe "sort-points"
  (it "should return a map of correctly sorted points."
      (should= {:inside  inside  
                :outside outside}
               (sort-points points))))

(describe "circle"
  (it "should equal 1 when x is 0."
    (should= 1.0 (circle 0)))
  (it "should equal 0 when x is 1."
      (should= 0.0 (circle 1)))
  (it "should equal 0.866 when x is 0.5"
      (should= "0.866" (format "%.3f" (circle 0.5)))))

(describe "draw-circle"
  (it "should be a JFreeChart object."
      (should= "class org.jfree.chart.JFreeChart"
               (str (.getClass (draw-circle))))))

(describe "generate-and-sort-points"
  (it "should return the correct number of points, sorted."
      (let [points (generate-and-sort-points 10)]
        (should= 2 (count points))
        (should= 10 (+ (count (:inside points))
                       (count (:outside points))))
        (should (not (nil? (:inside points))))
        (should (not (nil? (:outside points)))))))

(describe "label"
  (it "should return the correct label for a given point"
      (should= "inside" (label [0.1 0.3]))
      (should= "outside" (label [0.9 0.8]))))


(describe "estimate-pi"
 (it "should return four times the ratio of inside to outside points."
      (should= 2.0 (estimate-pi points)))

  (it "should return something close to pi when given a lot of points."
      (let  [estimate (estimate-pi (generate-random-points 70000))]
        (should (> estimate 3.13))
        (should (< estimate 3.15)))))

(run-specs)

project.clj

(defproject montecarlopi "0.1.0-SNAPSHOT"
  :dependencies [[org.clojure/clojure "1.5.0"]
                 [incanter "1.5.0-SNAPSHOT"]]
  :main montecarlopi.core
  :profiles {:dev {:dependencies [[speclj "2.5.0"]]}}
  :plugins [[speclj "2.5.0"]]
  :test-paths ["spec/"])