montecarlo.md

package org.entish.montecarlo
import scala.util.Random
import System.{ currentTimeMillis ⇒ _time }
import com.nicta.rng._, scalaz._, Scalaz._

object Main extends App {

  def profile[R](code: ⇒ R, t: Long = _time, msg: String = "Profile") = (msg, code, _time - t)

  def row[R](i: Int, p: (String, R, Long)) = s"${p._1} |  $i | ${p._2} | ${p._3}"

  def circleTest() = {
    val (x, y) = (Random.nextDouble, Random.nextDouble)
    x * x + y * y <= 1
  }

  def monteCarloRecurse(trials: Int, test: () ⇒ Boolean) = {
    def bool2double(b: Boolean) = if (b) 1.0d else 0.0d
    @scala.annotation.tailrec
    def recurse(n: Int, sum: Double): Double =
      if (n <= 0) sum / trials
      else recurse(n - 1, sum + bool2double(test()))
    recurse(trials, 0.0d)
  }

  def monteCarloMap(trials: Int, test: () ⇒ Boolean) =
    (1 to trials).map(_ ⇒ if (test()) 1.0d else 0.0d).sum / trials

  def monteCarloFold(trials: Int, test: () ⇒ Boolean) =
    (1 to trials).foldLeft(0.0d)((s, i) ⇒ s + (if (test()) 1.0d else 0.0d)) / trials

  def monteCarloIterator(trials: Int, test: () ⇒ Boolean) =
    (1 to trials).iterator.foldLeft(0.0d)((s, i) ⇒ s + (if (test()) 1.0d else 0.0d)) / trials

  def monteCarloStream(trials: Int, test: () ⇒ Boolean) =
    Stream.continually(if (test()) 1.0 else 0.0).take(trials).reduceLeft(_ + _) / trials

  def monteCarloParallel(trials: Int, test: () ⇒ Boolean) = {
    val pr = (1 to trials).par
    val s = pr.aggregate(0.0d)((a, _) ⇒ a + (if (test()) 1.0d else 0.0d), _ + _)
    s / trials
  }

  val pointInUnitSquare = Rng.choosedouble(0.0, 1.0) zip Rng.choosedouble(0.0, 1.0)

  val insideCircle = pointInUnitSquare.map { case (x, y) ⇒ x * x + y * y <= 1 }

  def monteCarloPurelyFunctional(trials: Int): Rng[Double] =
    EphemeralStream.range(0, trials).foldLeftM(0) {
      case (acc, _) ⇒ insideCircle.map(_.fold(1, 0) + acc)
    }.map(_ / trials.toDouble)

  (3 to 7).foreach { i ⇒
    val runs = math.pow(10, i).toInt
    println()
    println("Method | Runs | Result | Milliseconds")
    println("-------|------|--------|------------")
    println(row(runs, profile({ monteCarloMap(runs, circleTest) }, msg = "Map over range")))
    println(row(runs, profile({ monteCarloFold(runs, circleTest) }, msg = "Fold over range")))
    println(row(runs, profile({ monteCarloIterator(runs, circleTest) }, msg = "Fold over iterator")))
    println(row(runs, profile({ monteCarloStream(runs, circleTest) }, msg = "Fold over stream")))
    println(row(runs, profile({ monteCarloParallel(runs, circleTest) }, msg = "Aggregate over parallel")))
    println(row(runs, profile({ monteCarloRecurse(runs, circleTest) }, msg = "Recurse")))
    println(row(runs, profile({ monteCarloPurelyFunctional(runs).run.unsafePerformIO }, msg = "Purely Functional")))

  }
}

Here are the results for 10^7 runs:

Method	Runs	Result	Milliseconds
Map over range	10000000	0.7854594	964
Fold over range	10000000	0.7855624	519
Fold over iterator	10000000	0.7851361	666
Fold over stream	10000000	0.785246	1024
Aggregate over parallel	10000000	0.7852919	4990
Recurse	10000000	0.7854042	420
Purely Functional	10000000	0.7852436	9463