.gitlab-ci.yml

@@ -25,7 +25,7 @@ grade:
   tags:
     - cs206
   image:
-    name: smarter3/moocs:reactive-actorbintree-2020-04-15
+    name: smarter3/moocs:parprog1-barneshut-2020-03-09
     entrypoint: [""]
   allow_failure: true
   before_script:

assignment.sbt

@@ -1,4 +1,9 @@
 // Student tasks (i.e. submit, packageSubmission)
 enablePlugins(StudentTasks)
 
-
+courseraId := ch.epfl.lamp.CourseraId(
+  key = "itfW99oJEeWXuxJgUJEB-Q",
+  itemId = "z1ugn",
+  premiumItemId = Some("xGkV0"),
+  partId = "ep95q"
+)

build.sbt

@@ -1,24 +1,13 @@
-course := "reactive"
-assignment := "actorbintree"
-
-testOptions in Test += Tests.Argument(TestFrameworks.JUnit, "-a", "-v", "-s")
-parallelExecution in Test := false
-
-val akkaVersion = "2.6.0"
+course := "parprog1"
+assignment := "barneshut"
 
 scalaVersion := "0.23.0-bin-20200211-5b006fb-NIGHTLY"
-
-scalacOptions ++= Seq(
-  "-feature",
-  "-deprecation",
-  "-encoding", "UTF-8",
-  "-unchecked",
-  "-language:implicitConversions"
-)
-
+scalacOptions ++= Seq("-language:implicitConversions", "-deprecation")
 libraryDependencies ++= Seq(
-  "com.typesafe.akka"        %% "akka-actor"         % akkaVersion,
-  "com.typesafe.akka"        %% "akka-testkit"       % akkaVersion % Test,
-  "com.novocode"             % "junit-interface"     % "0.11"      % Test
+  "com.storm-enroute" %% "scalameter-core" % "0.19",
+  "org.scala-lang.modules" %% "scala-parallel-collections" % "0.2.0",
+  "com.novocode" % "junit-interface" % "0.11" % Test
 ).map(_.withDottyCompat(scalaVersion.value))
-testSuite := "actorbintree.BinaryTreeSuite"
+
+testOptions in Test += Tests.Argument(TestFrameworks.JUnit, "-a", "-v", "-s")
+testSuite := "barneshut.BarnesHutSuite"

src/main/scala/actorbintree/BinaryTreeSet.scala

@@ -1,189 +0,0 @@
-/**
- * Copyright (C) 2009-2013 Typesafe Inc. <http://www.typesafe.com>
- */
-package actorbintree
-
-import akka.actor._
-import scala.collection.immutable.Queue
-
-object BinaryTreeSet {
-
-  trait Operation {
-    def requester: ActorRef
-    def id: Int
-    def elem: Int
-  }
-
-  trait OperationReply {
-    def id: Int
-  }
-
-  /** Request with identifier `id` to insert an element `elem` into the tree.
-    * The actor at reference `requester` should be notified when this operation
-    * is completed.
-    */
-  case class Insert(requester: ActorRef, id: Int, elem: Int) extends Operation
-
-  /** Request with identifier `id` to check whether an element `elem` is present
-    * in the tree. The actor at reference `requester` should be notified when
-    * this operation is completed.
-    */
-  case class Contains(requester: ActorRef, id: Int, elem: Int) extends Operation
-
-  /** Request with identifier `id` to remove the element `elem` from the tree.
-    * The actor at reference `requester` should be notified when this operation
-    * is completed.
-    */
-  case class Remove(requester: ActorRef, id: Int, elem: Int) extends Operation
-
-  /** Request to perform garbage collection */
-  case object GC
-
-  /** Holds the answer to the Contains request with identifier `id`.
-    * `result` is true if and only if the element is present in the tree.
-    */
-  case class ContainsResult(id: Int, result: Boolean) extends OperationReply
-
-  /** Message to signal successful completion of an insert or remove operation. */
-  case class OperationFinished(id: Int) extends OperationReply
-
-}
-
-
-class BinaryTreeSet extends Actor {
-  import BinaryTreeSet._
-  import BinaryTreeNode._
-
-  def createRoot: ActorRef = context.actorOf(BinaryTreeNode.props(0, initiallyRemoved = true))
-
-  var root = createRoot
-
-  // optional
-  var pendingQueue = Queue.empty[Operation]
-
-  // optional
-  def receive = normal
-
-  // optional
-  /** Accepts `Operation` and `GC` messages. */
-  val normal: Receive = {
-    case op:Operation => root ! op
-    case GC => {
-      val newRoot = createRoot;
-      root ! CopyTo(newRoot)
-      context.become(garbageCollecting(newRoot))
-    }
-   }
-
-  // optional
-  /** Handles messages while garbage collection is performed.
-    * `newRoot` is the root of the new binary tree where we want to copy
-    * all non-removed elements into.
-    */
-  def garbageCollecting(newRoot: ActorRef): Receive = {
-    case op:Operation => pendingQueue = pendingQueue.enqueue(op)
-    case CopyFinished =>
-      pendingQueue.foreach(newRoot ! _) //foreach preserves order of a queue (same as dequeueing)
-      root ! PoisonPill //Will also stop all of its children
-      pendingQueue = Queue.empty
-      root = newRoot;
-      context.become(normal)
-    //Ignore GC messages here
-  }
-}
-
-object BinaryTreeNode {
-  trait Position
-
-  case object Left extends Position
-  case object Right extends Position
-
-  case class CopyTo(treeNode: ActorRef)
-  case object CopyFinished
-
-  def props(elem: Int, initiallyRemoved: Boolean) = Props(classOf[BinaryTreeNode],  elem, initiallyRemoved)
-}
-
-class BinaryTreeNode(val elem: Int, initiallyRemoved: Boolean) extends Actor {
-  import BinaryTreeNode._
-  import BinaryTreeSet._
-
-  var subtrees = Map[Position, ActorRef]()
-  var removed = initiallyRemoved
-
-  // optional
-  def receive = normal
-
-  def goDownTo(elem : Int) : Position = if(elem < this.elem) Left else Right
-  // optional
-  /** Handles `Operation` messages and `CopyTo` requests. */
-  val normal: Receive = {
-    case Insert  (requester, id, elem) =>
-      if(elem == this.elem && !removed){
-        requester ! OperationFinished(id)
-      }else{
-        val nextPos = goDownTo(elem)
-        
-        subtrees get nextPos match{
-          case Some(node) => node ! Insert(requester, id, elem)
-          case None => {
-              val newActorSubtree = (nextPos, context.actorOf(BinaryTreeNode.props(elem, false)))
-              subtrees = subtrees + newActorSubtree
-              requester ! OperationFinished(id);
-          }
-        }
-      }
-    case Contains(requester, id, elem) =>
-      if(elem == this.elem && !removed)
-        requester ! ContainsResult(id, true)
-      else{
-        //Need to search subtrees
-        subtrees get goDownTo(elem) match{
-          case Some(node) => node ! Contains(requester, id, elem)
-          case None => requester ! ContainsResult(id, false)
-        }
-      }
-        
-    case Remove  (requester, id, elem) => 
-      if(elem == this.elem && !removed){
-        removed = true
-        requester ! OperationFinished(id)
-      }else{
-        subtrees get goDownTo(elem) match{
-          case Some(node) => node ! Remove(requester, id, elem)
-          case None => requester ! OperationFinished(id) // (elem isn't in the tree)
-        }
-      }
-
-    case CopyTo(newRoot) =>
-      //We are already done, nothing to do
-      if(removed && subtrees.isEmpty) context.parent ! CopyFinished
-      else{
-        if(!removed) newRoot ! Insert(self, elem, elem)
-        subtrees.values foreach(_ ! CopyTo(newRoot)) //Copy subtrees elems
-        //val insertConfirmed = if(removed) true else false, hence we can simply pass removed
-        context.become(copying(subtrees.values.toSet, removed))
-      }
-  }
-
-  // optional
-  /** `expected` is the set of ActorRefs whose replies we are waiting for,
-    * `insertConfirmed` tracks whether the copy of this node to the new tree has been confirmed.
-    */
-  def copying(expected: Set[ActorRef], insertConfirmed: Boolean): Receive = {
-    //To catch the insert of this node into the new tree beeing finished
-    case OperationFinished(_) => {
-      if(expected.isEmpty) context.parent ! CopyFinished
-      else context.become(copying(expected, true))
-    }
-    
-    case CopyFinished => {
-      val newExp = expected-sender
-      if(insertConfirmed && newExp.isEmpty){
-        context.parent ! CopyFinished
-      }else{
-        context.become(copying(newExp, insertConfirmed))
-      }
-    }
-  }
-}

src/main/scala/barneshut/BarnesHut.scala

@@ -0,0 +1,151 @@
+package barneshut
+
+import java.awt._
+import java.awt.event._
+import javax.swing._
+import javax.swing.event._
+import scala.collection.parallel._
+import scala.collection.mutable.ArrayBuffer
+import scala.reflect.ClassTag
+
+object BarnesHut {
+
+  val model = new SimulationModel
+
+  var simulator: Simulator = _
+
+  def initialize(parallelismLevel: Int, pattern: String, nbodies: Int): Unit = {
+    model.initialize(parallelismLevel, pattern, nbodies)
+    model.timeStats.clear()
+    simulator = new Simulator(model.taskSupport, model.timeStats)
+  }
+
+  class BarnesHutFrame extends JFrame("Barnes-Hut") {
+    setDefaultCloseOperation(WindowConstants.EXIT_ON_CLOSE)
+    setSize(1024, 600)
+    setLayout(new BorderLayout)
+
+    val rightpanel = new JPanel
+    rightpanel.setBorder(BorderFactory.createEtchedBorder(border.EtchedBorder.LOWERED))
+    rightpanel.setLayout(new BorderLayout)
+    add(rightpanel, BorderLayout.EAST)
+
+    val controls = new JPanel
+    controls.setLayout(new GridLayout(0, 2))
+    rightpanel.add(controls, BorderLayout.NORTH)
+
+    val parallelismLabel = new JLabel("Parallelism")
+    controls.add(parallelismLabel)
+
+    val items = (1 to Runtime.getRuntime.availableProcessors).map(_.toString).toArray
+    val parcombo = new JComboBox[String](items)
+    parcombo.setSelectedIndex(items.length - 1)
+    parcombo.addActionListener(new ActionListener {
+      def actionPerformed(e: ActionEvent) = {
+        initialize(getParallelism, "two-galaxies", getTotalBodies)
+        canvas.repaint()
+      }
+    })
+    controls.add(parcombo)
+
+    val bodiesLabel = new JLabel("Total bodies")
+    controls.add(bodiesLabel)
+
+    val bodiesSpinner = new JSpinner(new SpinnerNumberModel(25000, 32, 1000000, 1000))
+    bodiesSpinner.addChangeListener(new ChangeListener {
+      def stateChanged(e: ChangeEvent) = {
+        if (frame != null) {
+          initialize(getParallelism, "two-galaxies", getTotalBodies)
+          canvas.repaint()
+        }
+      }
+    })
+    controls.add(bodiesSpinner)
+
+    val stepbutton = new JButton("Step")
+    stepbutton.addActionListener(new ActionListener {
+      def actionPerformed(e: ActionEvent): Unit = {
+        stepThroughSimulation()
+      }
+    })
+    controls.add(stepbutton)
+
+    val startButton = new JToggleButton("Start/Pause")
+    val startTimer = new javax.swing.Timer(0, new ActionListener {
+      def actionPerformed(e: ActionEvent): Unit = {
+        stepThroughSimulation()
+      }
+    })
+    startButton.addActionListener(new ActionListener {
+      def actionPerformed(e: ActionEvent): Unit = {
+        if (startButton.isSelected) startTimer.start()
+        else startTimer.stop()
+      }
+    })
+    controls.add(startButton)
+
+    val quadcheckbox = new JToggleButton("Show quad")
+    quadcheckbox.addActionListener(new ActionListener {
+      def actionPerformed(e: ActionEvent): Unit = {
+        model.shouldRenderQuad = quadcheckbox.isSelected
+        repaint()
+      }
+    })
+    controls.add(quadcheckbox)
+
+    val clearButton = new JButton("Restart")
+    clearButton.addActionListener(new ActionListener {
+      def actionPerformed(e: ActionEvent): Unit = {
+        initialize(getParallelism, "two-galaxies", getTotalBodies)
+      }
+    })
+    controls.add(clearButton)
+
+    val info = new JTextArea("   ")
+    info.setBorder(BorderFactory.createLoweredBevelBorder)
+    rightpanel.add(info, BorderLayout.SOUTH)
+
+    val canvas = new SimulationCanvas(model)
+    add(canvas, BorderLayout.CENTER)
+    setVisible(true)
+
+    def updateInformationBox(): Unit = {
+      val text = model.timeStats.toString
+      frame.info.setText("--- Statistics: ---\n" + text)
+    }
+
+    def stepThroughSimulation(): Unit = {
+      SwingUtilities.invokeLater(new Runnable {
+        def run() = {
+          val (bodies, quad) = simulator.step(model.bodies)
+          model.bodies = bodies
+          model.quad = quad
+          updateInformationBox()
+          repaint()
+        }
+      })
+    }
+
+    def getParallelism = {
+      val selidx = parcombo.getSelectedIndex
+      parcombo.getItemAt(selidx).toInt
+    }
+
+    def getTotalBodies = bodiesSpinner.getValue.asInstanceOf[Int]
+
+    initialize(getParallelism, "two-galaxies", getTotalBodies)
+  }
+
+  try {
+    UIManager.setLookAndFeel(UIManager.getSystemLookAndFeelClassName())
+  } catch {
+    case _: Exception => println("Cannot set look and feel, using the default one.")
+  }
+
+  val frame = new BarnesHutFrame
+
+  def main(args: Array[String]): Unit = {
+    frame.repaint()
+  }
+
+}

src/main/scala/barneshut/Interfaces.scala

@@ -0,0 +1,23 @@
+package barneshut
+
+import conctrees.ConcBuffer
+
+// Interfaces used by the grading infrastructure. Do not change signatures
+// or your submission will fail with a NoSuchMethodError.
+
+trait SectorMatrixInterface {
+  def +=(b: Body): SectorMatrix
+  def combine(that: SectorMatrix): SectorMatrix
+  def apply(x: Int, y: Int): ConcBuffer[Body]
+}
+
+trait QuadInterface {
+  def massX: Float
+  def massY: Float
+  def mass: Float
+  def centerX: Float
+  def centerY: Float
+  def size: Float
+  def total: Int
+  def insert(b: Body): Quad
+}

src/main/scala/barneshut/SimulationCanvas.scala

@@ -0,0 +1,146 @@
+package barneshut
+
+import java.awt._
+import java.awt.event._
+import javax.swing._
+import javax.swing.event._
+
+class SimulationCanvas(val model: SimulationModel) extends JComponent {
+
+  val MAX_RES = 3000
+
+  val pixels = new Array[Int](MAX_RES * MAX_RES)
+
+  override def paintComponent(gcan: Graphics) = {
+    super.paintComponent(gcan)
+
+    val width = getWidth
+    val height = getHeight
+    val img = new image.BufferedImage(width, height, image.BufferedImage.TYPE_INT_ARGB)
+
+    // clear canvas pixels
+    for (x <- 0 until MAX_RES; y <- 0 until MAX_RES) pixels(y * width + x) = 0
+
+    // count number of bodies in each pixel
+    for (b <- model.bodies) {
+      val px = ((b.x - model.screen.minX) / model.screen.width * width).toInt
+      val py = ((b.y - model.screen.minY) / model.screen.height * height).toInt
+      if (px >= 0 && px < width && py >= 0 && py < height) pixels(py * width + px) += 1
+    }
+
+    // set image intensity depending on the number of bodies in the pixel
+    for (y <- 0 until height; x <- 0 until width) {
+      val count = pixels(y * width + x)
+      val intensity = if (count > 0) math.min(255, 70 + count * 50) else 0
+      val color = (255 << 24) | (intensity << 16) | (intensity << 8) | intensity
+      img.setRGB(x, y, color)
+    }
+
+    // for debugging purposes, if the number of bodies is small, output their locations
+    val g = img.getGraphics.asInstanceOf[Graphics2D]
+    g.setColor(Color.GRAY)
+    if (model.bodies.length < 350) for (b <- model.bodies) {
+      def round(x: Float) = (x * 100).toInt / 100.0f
+      val px = ((b.x - model.screen.minX) / model.screen.width * width).toInt
+      val py = ((b.y - model.screen.minY) / model.screen.height * height).toInt
+      if (px >= 0 && px < width && py >= 0 && py < height) {
+        g.drawString(s"${round(b.x)}, ${round(b.y)}", px, py)
+      }
+    }
+
+    // render quad if necessary
+    if (model.shouldRenderQuad) {
+      g.setRenderingHint(RenderingHints.KEY_ANTIALIASING, RenderingHints.VALUE_ANTIALIAS_ON)
+      val green = new Color(0, 225, 80, 150)
+      val red = new Color(200, 0, 0, 150)
+      g.setColor(green)
+      def drawQuad(depth: Int, quad: Quad): Unit = {
+        def drawRect(fx: Float, fy: Float, fsz: Float, q: Quad, fill: Boolean = false): Unit = {
+          val x = ((fx - model.screen.minX) / model.screen.width * width).toInt
+          val y = ((fy - model.screen.minY) / model.screen.height * height).toInt
+          val w = ((fx + fsz - model.screen.minX) / model.screen.width * width).toInt - x
+          val h = ((fy + fsz - model.screen.minY) / model.screen.height * height).toInt - y
+          g.drawRect(x, y, w, h)
+          if (fill) g.fillRect(x, y, w, h)
+          if (depth <= 5) g.drawString("#:" + q.total, x + w / 2, y + h / 2)
+        }
+        quad match {
+          case Fork(nw, ne, sw, se) =>
+            val cx = quad.centerX
+            val cy = quad.centerY
+            val sz = quad.size
+            drawRect(cx - sz / 2, cy - sz / 2, sz / 2, nw)
+            drawRect(cx - sz / 2, cy, sz / 2, sw)
+            drawRect(cx, cy - sz / 2, sz / 2, ne)
+            drawRect(cx, cy, sz / 2, se)
+            drawQuad(depth + 1, nw)
+            drawQuad(depth + 1, ne)
+            drawQuad(depth + 1, sw)
+            drawQuad(depth + 1, se)
+          case Empty(_, _, _) | Leaf(_, _, _, _) =>
+            // done
+        }
+      }
+      drawQuad(0, model.quad)
+
+    }
+    gcan.drawImage(img, 0, 0, null)
+  }
+
+  // zoom on mouse rotation
+  addMouseWheelListener(new MouseAdapter {
+    override def mouseWheelMoved(e: MouseWheelEvent): Unit = {
+      val rot = e.getWheelRotation
+      val cx = model.screen.centerX
+      val cy = model.screen.centerY
+      val w = model.screen.width
+      val h = model.screen.height
+      val factor = {
+        if (rot > 0) 0.52f
+        else if (rot < 0) 0.48f
+        else 0.5f
+      }
+      model.screen.minX = cx - w * factor
+      model.screen.minY = cy - h * factor
+      model.screen.maxX = cx + w * factor
+      model.screen.maxY = cy + h * factor
+      repaint()
+    }
+  })
+
+  // reset the last known mouse drag position on mouse press
+  var xlast = Int.MinValue
+  var ylast = Int.MinValue
+  addMouseListener(new MouseAdapter {
+    override def mousePressed(e: MouseEvent): Unit = {
+      xlast = Int.MinValue
+      ylast = Int.MinValue
+    }
+  })
+
+  // update the last known mouse drag position on mouse drag,
+  // update the boundaries of the visible area
+  addMouseMotionListener(new MouseMotionAdapter {
+    override def mouseDragged(e: MouseEvent): Unit = {
+      val xcurr = e.getX
+      val ycurr = e.getY
+      if (xlast != Int.MinValue) {
+        val xd = xcurr - xlast
+        val yd = ycurr - ylast
+        val w = model.screen.width
+        val h = model.screen.height
+        val cx = model.screen.centerX - xd * w / 1000
+        val cy = model.screen.centerY - yd * h / 1000
+        model.screen.minX = cx - w / 2
+        model.screen.minY = cy - h / 2
+        model.screen.maxX = cx + w / 2
+        model.screen.maxY = cy + h / 2
+        println(model.screen)
+      }
+      xlast = xcurr
+      ylast = ycurr
+      repaint()
+    }
+  })
+
+}

src/main/scala/barneshut/SimulationModel.scala

@@ -0,0 +1,73 @@
+package barneshut
+
+import java.awt._
+import java.awt.event._
+import javax.swing._
+import javax.swing.event._
+import scala.collection.parallel.{TaskSupport, defaultTaskSupport}
+import scala.{collection => coll}
+
+class SimulationModel {
+
+  var screen = new Boundaries
+
+  var bodies: coll.Seq[Body] = Nil
+
+  var quad: Quad = Empty(screen.centerX, screen.centerY, Float.MaxValue)
+
+  var shouldRenderQuad = false
+
+  var timeStats = new TimeStatistics
+
+  var taskSupport: TaskSupport = defaultTaskSupport
+
+  def initialize(parallelismLevel: Int, pattern: String, totalBodies: Int): Unit = {
+    taskSupport = new collection.parallel.ForkJoinTaskSupport(
+      new java.util.concurrent.ForkJoinPool(parallelismLevel))
+
+    pattern match {
+      case "two-galaxies" => init2Galaxies(totalBodies)
+      case _ => sys.error(s"no such initial pattern: $pattern")
+    }
+  }
+
+  def init2Galaxies(totalBodies: Int): Unit = {
+    val bodyArray = new Array[Body](totalBodies)
+    val random = new scala.util.Random(213L)
+
+    def galaxy(from: Int, num: Int, maxradius: Float, cx: Float, cy: Float, sx: Float, sy: Float): Unit = {
+      val totalM = 1.5f * num
+      val blackHoleM = 1.0f * num
+      val cubmaxradius = maxradius * maxradius * maxradius
+      for (i <- from until (from + num)) {
+        val b = if (i == from) {
+          new Body(blackHoleM, cx, cy, sx, sy)
+        } else {
+          val angle = random.nextFloat * 2 * math.Pi
+          val radius = 25 + maxradius * random.nextFloat
+          val starx = cx + radius * math.sin(angle).toFloat
+          val stary = cy + radius * math.cos(angle).toFloat
+          val speed = math.sqrt(gee * blackHoleM / radius + gee * totalM * radius * radius / cubmaxradius)
+          val starspeedx = sx + (speed * math.sin(angle + math.Pi / 2)).toFloat
+          val starspeedy = sy + (speed * math.cos(angle + math.Pi / 2)).toFloat
+          val starmass = 1.0f + 1.0f * random.nextFloat
+          new Body(starmass, starx, stary, starspeedx, starspeedy)
+        }
+        bodyArray(i) = b
+      }
+    }
+
+    galaxy(0, bodyArray.length / 8, 300.0f, 0.0f, 0.0f, 0.0f, 0.0f)
+    galaxy(bodyArray.length / 8, bodyArray.length / 8 * 7, 350.0f, -1800.0f, -1200.0f, 0.0f, 0.0f)
+
+    bodies = bodyArray.toSeq
+
+    // compute center and boundaries
+    screen = new Boundaries
+    screen.minX = -2200.0f
+    screen.minY = -1600.0f
+    screen.maxX = 350.0f
+    screen.maxY = 350.0f
+  }
+
+}

src/main/scala/barneshut/Simulator.scala

@@ -0,0 +1,112 @@
+package barneshut
+
+import java.awt._
+import java.awt.event._
+import javax.swing._
+import javax.swing.event._
+import scala.{collection => coll}
+import scala.collection.parallel.{TaskSupport, Combiner}
+import scala.collection.parallel.mutable.ParHashSet
+import scala.collection.parallel.CollectionConverters._
+
+class Simulator(val taskSupport: TaskSupport, val timeStats: TimeStatistics) {
+
+  def updateBoundaries(boundaries: Boundaries, body: Body): Boundaries = {
+    boundaries.minX = Math.min(boundaries.minX, body.x);
+    boundaries.minY = Math.min(boundaries.minY, body.y);
+    boundaries.maxX = Math.max(boundaries.maxX, body.x);
+    boundaries.maxY = Math.max(boundaries.maxY, body.y);
+    boundaries
+  }
+
+  def mergeBoundaries(a: Boundaries, b: Boundaries): Boundaries = {
+    val bnd = new Boundaries
+    bnd.minX = Math.min(a.minX, b.minX)
+    bnd.minY = Math.min(a.minY, b.minY)
+    bnd.maxX = Math.max(a.maxX, b.maxX)
+    bnd.maxY = Math.max(a.maxY, b.maxY)
+    bnd
+  }
+
+  def computeBoundaries(bodies: coll.Seq[Body]): Boundaries = timeStats.timed("boundaries") {
+    val parBodies = bodies.par
+    parBodies.tasksupport = taskSupport
+    parBodies.aggregate(new Boundaries)(updateBoundaries, mergeBoundaries)
+  }
+
+  def computeSectorMatrix(bodies: coll.Seq[Body], boundaries: Boundaries): SectorMatrix = timeStats.timed("matrix") {
+    val parBodies = bodies.par
+    parBodies.tasksupport = taskSupport
+    parBodies.aggregate(new SectorMatrix(boundaries, SECTOR_PRECISION))((accSM, b) => accSM += b, (sm1, sm2) => sm1.combine(sm2))
+  }
+
+  def computeQuad(sectorMatrix: SectorMatrix): Quad = timeStats.timed("quad") {
+    sectorMatrix.toQuad(taskSupport.parallelismLevel)
+  }
+
+  def updateBodies(bodies: coll.Seq[Body], quad: Quad): coll.Seq[Body] = timeStats.timed("update") {
+    val parBodies = bodies.par
+    parBodies.tasksupport = taskSupport
+    parBodies.map(_.updated(quad)).seq
+  }
+
+  def eliminateOutliers(bodies: coll.Seq[Body], sectorMatrix: SectorMatrix, quad: Quad): coll.Seq[Body] = timeStats.timed("eliminate") {
+    def isOutlier(b: Body): Boolean = {
+      val dx = quad.massX - b.x
+      val dy = quad.massY - b.y
+      val d = math.sqrt(dx * dx + dy * dy)
+      // object is far away from the center of the mass
+      if (d > eliminationThreshold * sectorMatrix.boundaries.size) {
+        val nx = dx / d
+        val ny = dy / d
+        val relativeSpeed = b.xspeed * nx + b.yspeed * ny
+        // object is moving away from the center of the mass
+        if (relativeSpeed < 0) {
+          val escapeSpeed = math.sqrt(2 * gee * quad.mass / d)
+          // object has the espace velocity
+          -relativeSpeed > 2 * escapeSpeed
+        } else false
+      } else false
+    }
+
+    def outliersInSector(x: Int, y: Int): Combiner[Body, ParHashSet[Body]] = {
+      val combiner = ParHashSet.newCombiner[Body]
+      combiner ++= sectorMatrix(x, y).filter(isOutlier)
+      combiner
+    }
+
+    val sectorPrecision = sectorMatrix.sectorPrecision
+    val horizontalBorder = for (x <- 0 until sectorPrecision; y <- Seq(0, sectorPrecision - 1)) yield (x, y)
+    val verticalBorder = for (y <- 1 until sectorPrecision - 1; x <- Seq(0, sectorPrecision - 1)) yield (x, y)
+    val borderSectors = horizontalBorder ++ verticalBorder
+
+    // compute the set of outliers
+    val parBorderSectors = borderSectors.par
+    parBorderSectors.tasksupport = taskSupport
+    val outliers = parBorderSectors.map({ case (x, y) => outliersInSector(x, y) }).reduce(_ combine _).result
+
+    // filter the bodies that are outliers
+    val parBodies = bodies.par
+    parBodies.filter(!outliers(_)).seq
+  }
+
+  def step(bodies: coll.Seq[Body]): (coll.Seq[Body], Quad) = {
+    // 1. compute boundaries
+    val boundaries = computeBoundaries(bodies)
+
+    // 2. compute sector matrix
+    val sectorMatrix = computeSectorMatrix(bodies, boundaries)
+
+    // 3. compute quad tree
+    val quad = computeQuad(sectorMatrix)
+
+    // 4. eliminate outliers
+    val filteredBodies = eliminateOutliers(bodies, sectorMatrix, quad)
+
+    // 5. update body velocities and positions
+    val newBodies = updateBodies(filteredBodies, quad)
+
+    (newBodies, quad)
+  }
+
+}

src/main/scala/barneshut/conctrees/Conc.scala

@@ -0,0 +1,148 @@
+package barneshut
+package conctrees
+
+import scala.annotation.tailrec
+
+sealed trait Conc[@specialized(Int, Long, Float, Double) +T] {
+  def level: Int
+  def size: Int
+  def left: Conc[T]
+  def right: Conc[T]
+  def normalized = this
+}
+
+object Conc {
+
+  case class <>[+T](left: Conc[T], right: Conc[T]) extends Conc[T] {
+    val level = 1 + math.max(left.level, right.level)
+    val size = left.size + right.size
+  }
+
+  sealed trait Leaf[T] extends Conc[T] {
+    def left = sys.error("Leaves do not have children.")
+    def right = sys.error("Leaves do not have children.")
+  }
+
+  case object Empty extends Leaf[Nothing] {
+    def level = 0
+    def size = 0
+  }
+
+  class Single[@specialized(Int, Long, Float, Double) T](val x: T) extends Leaf[T] {
+    def level = 0
+    def size = 1
+    override def toString = s"Single($x)"
+  }
+
+  class Chunk[@specialized(Int, Long, Float, Double) T](val array: Array[T], val size: Int, val k: Int)
+  extends Leaf[T] {
+    def level = 0
+    override def toString = s"Chunk(${array.mkString("", ", ", "")}; $size; $k)"
+  }
+
+  case class Append[+T](left: Conc[T], right: Conc[T]) extends Conc[T] {
+    val level = 1 + math.max(left.level, right.level)
+    val size = left.size + right.size
+    override def normalized = {
+      def wrap[T](xs: Conc[T], ys: Conc[T]): Conc[T] = (xs: @unchecked) match {
+        case Append(ws, zs) => wrap(ws, zs <> ys)
+        case xs => xs <> ys
+      }
+      wrap(left, right)
+    }
+  }
+
+  def concatTop[T](xs: Conc[T], ys: Conc[T]) = {
+    if (xs == Empty) ys
+    else if (ys == Empty) xs
+    else concat(xs, ys)
+  }
+
+  private def concat[T](xs: Conc[T], ys: Conc[T]): Conc[T] = {
+    val diff = ys.level - xs.level
+    if (diff >= -1 && diff <= 1) new <>(xs, ys)
+    else if (diff < -1) {
+      if (xs.left.level >= xs.right.level) {
+        val nr = concat(xs.right, ys)
+        new <>(xs.left, nr)
+      } else {
+        val nrr = concat(xs.right.right, ys)
+        if (nrr.level == xs.level - 3) {
+          val nl = xs.left
+          val nr = new <>(xs.right.left, nrr)
+          new <>(nl, nr)
+        } else {
+          val nl = new <>(xs.left, xs.right.left)
+          val nr = nrr
+          new <>(nl, nr)
+        }
+      }
+    } else {
+      if (ys.right.level >= ys.left.level) {
+        val nl = concat(xs, ys.left)
+        new <>(nl, ys.right)
+      } else {
+        val nll = concat(xs, ys.left.left)
+        if (nll.level == ys.level - 3) {
+          val nl = new <>(nll, ys.left.right)
+          val nr = ys.right
+          new <>(nl, nr)
+        } else {
+          val nl = nll
+          val nr = new <>(ys.left.right, ys.right)
+          new <>(nl, nr)
+        }
+      }
+    }
+  }
+
+  def appendTop[T](xs: Conc[T], ys: Leaf[T]): Conc[T] = (xs: @unchecked) match {
+    case xs: Append[T] => append(xs, ys)
+    case _ <> _ => new Append(xs, ys)
+    case Empty => ys
+    case xs: Leaf[T] => new <>(xs, ys)
+  }
+  @tailrec private def append[T](xs: Append[T], ys: Conc[T]): Conc[T] = {
+    if (xs.right.level > ys.level) new Append(xs, ys)
+    else {
+      val zs = new <>(xs.right, ys)
+      xs.left match {
+        case ws @ Append(_, _) => append(ws, zs)
+        case ws if ws.level <= zs.level => ws <> zs
+        case ws => new Append(ws, zs)
+      }
+    }
+  }
+
+  def traverse[@specialized(Int, Long, Float, Double) T, @specialized(Int, Long, Float, Double) U](xs: Conc[T], f: T => U): Unit = (xs: @unchecked) match {
+    case left <> right =>
+      traverse(left, f)
+      traverse(right, f)
+    case s: Single[T] =>
+      f(s.x)
+    case c: Chunk[T] =>
+      val a = c.array
+      val sz = c.size
+      var i = 0
+      while (i < sz) {
+        f(a(i))
+        i += 1
+      }
+    case Empty =>
+    case Append(left, right) =>
+      traverse(left, f)
+      traverse(right, f)
+    case _ =>
+      sys.error("All cases should have been covered: " + xs + ", " + xs.getClass)
+  }
+
+  def iterator[@specialized(Int, Long, Float, Double) T](xs: Conc[T]): Iterator[T] = (xs: @unchecked) match {
+    case left <> right       => iterator(left) ++ iterator(right)
+    case s: Single[T]        => Iterator.single(s.x)
+    case c: Chunk[T]         => c.array.iterator.take(c.size)
+    case Empty               => Iterator.empty
+    case Append(left, right) => iterator(left) ++ iterator(right)
+    case _                   => sys.error("All cases should have been covered: " + xs + ", " + xs.getClass)
+  }
+
+}

src/main/scala/barneshut/conctrees/ConcBuffer.scala

@@ -0,0 +1,86 @@
+package barneshut
+package conctrees
+
+import scala.reflect.ClassTag
+import scala.collection.parallel.CollectionConverters._
+import org.scalameter._
+
+class ConcBuffer[@specialized(Byte, Char, Int, Long, Float, Double) T: ClassTag](
+  val k: Int, private var conc: Conc[T]
+) extends Iterable[T] {
+  require(k > 0)
+
+  def this() = this(128, Conc.Empty)
+
+  private var chunk: Array[T] = new Array(k)
+  private var lastSize: Int = 0
+
+  def iterator: Iterator[T] = Conc.iterator(conc) ++ chunk.iterator.take(lastSize)
+
+  final def +=(elem: T): this.type = {
+    if (lastSize >= k) expand()
+    chunk(lastSize) = elem
+    lastSize += 1
+    this
+  }
+
+  final def combine(that: ConcBuffer[T]): ConcBuffer[T] = {
+    val combinedConc = this.result <> that.result
+    this.clear()
+    that.clear()
+    new ConcBuffer(k, combinedConc)
+  }
+
+  private def pack(): Unit = {
+    conc = Conc.appendTop(conc, new Conc.Chunk(chunk, lastSize, k))
+  }
+
+  private def expand(): Unit = {
+    pack()
+    chunk = new Array(k)
+    lastSize = 0
+  }
+
+  def clear(): Unit = {
+    conc = Conc.Empty
+    chunk = new Array(k)
+    lastSize = 0
+  }
+
+  def result: Conc[T] = {
+    pack()
+    conc
+  }
+}
+
+object ConcBufferRunner {
+
+  val standardConfig = config(
+    Key.exec.minWarmupRuns -> 20,
+    Key.exec.maxWarmupRuns -> 40,
+    Key.exec.benchRuns -> 60,
+    Key.verbose -> false
+  ).withWarmer(new Warmer.Default)
+
+  def main(args: Array[String]): Unit = {
+    val size = 1000000
+
+    def run(p: Int): Unit = {
+      val taskSupport = new collection.parallel.ForkJoinTaskSupport(
+        new java.util.concurrent.ForkJoinPool(p))
+      val strings = (0 until size).map(_.toString)
+      val time = standardConfig measure {
+        val parallelized = strings.par
+        parallelized.tasksupport = taskSupport
+        parallelized.aggregate(new ConcBuffer[String])(_ += _, _ combine _).result
+      }
+      println(s"p = $p, time = ${time.value}")
+    }
+
+    run(1)
+    run(2)
+    run(4)
+    run(8)
+  }
+
+}

src/main/scala/barneshut/conctrees/package.scala

@@ -0,0 +1,16 @@
+package barneshut
+
+import org.scalameter._
+
+package object conctrees {
+
+  implicit class ConcOps[T](val self: Conc[T]) extends AnyVal {
+    def foreach[U](f: T => U) = Conc.traverse(self, f)
+    def <>(that: Conc[T]) = Conc.concatTop(self.normalized, that.normalized)
+  }
+
+  // Workaround Dotty's handling of the existential type KeyValue
+  implicit def keyValueCoerce[T](kv: (Key[T], T)): KeyValue = {
+    kv.asInstanceOf[KeyValue]
+  }
+}

src/main/scala/barneshut/package.scala

@@ -0,0 +1,334 @@
+import java.util.concurrent._
+import scala.{collection => coll}
+import scala.util.DynamicVariable
+import barneshut.conctrees._
+
+package object barneshut {
+
+  class Boundaries {
+    var minX = Float.MaxValue
+
+    var minY = Float.MaxValue
+
+    var maxX = Float.MinValue
+
+    var maxY = Float.MinValue
+
+    def width = maxX - minX
+
+    def height = maxY - minY
+
+    def size = math.max(width, height)
+
+    def centerX = minX + width / 2
+
+    def centerY = minY + height / 2
+
+    override def toString = s"Boundaries($minX, $minY, $maxX, $maxY)"
+  }
+
+  sealed abstract class Quad extends QuadInterface {
+    def massX: Float
+
+    def massY: Float
+
+    def mass: Float
+
+    def centerX: Float
+
+    def centerY: Float
+
+    def size: Float
+
+    def total: Int
+
+    def insert(b: Body): Quad
+  }
+
+  case class Empty(centerX: Float, centerY: Float, size: Float) extends Quad {
+    def massX: Float = centerX
+    def massY: Float = centerY
+    def mass: Float = 0
+    def total: Int = 0
+    def insert(b: Body): Quad = Leaf(centerX, centerY, size, Seq(b))
+  }
+
+  case class Fork(
+    nw: Quad, ne: Quad, sw: Quad, se: Quad
+  ) extends Quad {
+
+    val centerX: Float = (nw.centerX + ne.centerX)/2
+    val centerY: Float = (nw.centerY + sw.centerY)/2
+    val size: Float = nw.size + ne.size
+    val mass: Float = Seq(nw, ne, sw, se).map(_.mass).sum
+    val massX: Float = if(mass == 0) centerX else Seq(nw, ne, sw, se).map(q => q.mass * q.massX).sum/mass
+    val massY: Float = if(mass == 0) centerY else Seq(nw, ne, sw, se).map(q => q.mass * q.massY).sum/mass
+    val total: Int = Seq(nw, ne, sw, se).map(_.total).sum;
+
+    def insert(b: Body): Fork = {
+      if(b.x <= centerX){ //W
+        if(b.y <= centerY){ // NW
+          Fork(nw.insert(b), ne, sw, se)
+        }
+        else{ //SW
+          Fork(nw, ne, sw.insert(b), se)
+        }
+      }
+      else{ //EAST
+        if(b.y <= centerY){ //NE
+          Fork(nw, ne.insert(b), sw, se)
+        }
+        else{ //SE
+          Fork(nw, ne, sw, se.insert(b))
+        }
+      }
+    }
+  }
+
+  case class Leaf(centerX: Float, centerY: Float, size: Float, bodies: coll.Seq[Body])
+  extends Quad {
+    val mass = bodies.map(_.mass).sum
+    val massX = bodies.map(b => b.mass * b.x).sum / mass
+    val massY = bodies.map(b => b.mass * b.y).sum / mass
+    val total: Int = bodies.length
+    def insert(b: Body): Quad = {
+      if(size > minimumSize){
+        val wCorr = centerX - size/4
+        val eCorr = centerX + size/4
+        val nCorr = centerY - size/4
+        val sCorr = centerY + size/4
+
+        val newSize = size/2
+
+        val fork = Fork(
+          Empty(wCorr, nCorr, newSize),
+          Empty(eCorr, nCorr, newSize),
+          Empty(wCorr, sCorr, newSize),
+          Empty(eCorr, sCorr, newSize))
+        (bodies :+ b).foldLeft(fork)((f, b) => f.insert(b))
+      }
+      else {
+        Leaf(centerX, centerY, size, bodies :+ b)
+      }
+    }
+  }
+
+  def minimumSize = 0.00001f
+
+  def gee: Float = 100.0f
+
+  def delta: Float = 0.01f
+
+  def theta = 0.5f
+
+  def eliminationThreshold = 0.5f
+
+  def force(m1: Float, m2: Float, dist: Float): Float = gee * m1 * m2 / (dist * dist)
+
+  def distance(x0: Float, y0: Float, x1: Float, y1: Float): Float = {
+    math.sqrt((x1 - x0) * (x1 - x0) + (y1 - y0) * (y1 - y0)).toFloat
+  }
+
+  class Body(val mass: Float, val x: Float, val y: Float, val xspeed: Float, val yspeed: Float) {
+
+    def updated(quad: Quad): Body = {
+      var netforcex = 0.0f
+      var netforcey = 0.0f
+
+      def addForce(thatMass: Float, thatMassX: Float, thatMassY: Float): Unit = {
+        val dist = distance(thatMassX, thatMassY, x, y)
+        /* If the distance is smaller than 1f, we enter the realm of close
+         * body interactions. Since we do not model them in this simplistic
+         * implementation, bodies at extreme proximities get a huge acceleration,
+         * and are catapulted from each other's gravitational pull at extreme
+         * velocities (something like this:
+         * http://en.wikipedia.org/wiki/Interplanetary_spaceflight#Gravitational_slingshot).
+         * To decrease the effect of this gravitational slingshot, as a very
+         * simple approximation, we ignore gravity at extreme proximities.
+         */
+        if (dist > 1f) {
+          val dforce = force(mass, thatMass, dist)
+          val xn = (thatMassX - x) / dist
+          val yn = (thatMassY - y) / dist
+          val dforcex = dforce * xn
+          val dforcey = dforce * yn
+          netforcex += dforcex
+          netforcey += dforcey
+        }
+      }
+
+      def traverse(quad: Quad): Unit = (quad: Quad) match {
+        case Empty(_, _, _) =>
+          // no force
+        case Leaf(_, _, _, bodies) => bodies.foreach(b => addForce(b.mass, b.x, b.y))
+          // add force contribution of each body by calling addForce
+        case Fork(nw, ne, sw, se) => if(quad.size / distance(x, y, quad.centerX, quad.centerY) < theta)
+        addForce(quad.mass, quad.massX, quad.massY)
+      else {
+        traverse(nw)
+        traverse(ne)
+        traverse(sw)
+        traverse(se)
+      }
+          // see if node is far enough from the body,
+          // or recursion is needed
+      }
+
+      traverse(quad)
+
+      val nx = x + xspeed * delta
+      val ny = y + yspeed * delta
+      val nxspeed = xspeed + netforcex / mass * delta
+      val nyspeed = yspeed + netforcey / mass * delta
+
+      new Body(mass, nx, ny, nxspeed, nyspeed)
+    }
+
+  }
+
+  val SECTOR_PRECISION = 8
+
+  class SectorMatrix(val boundaries: Boundaries, val sectorPrecision: Int) extends SectorMatrixInterface {
+    val sectorSize = boundaries.size / sectorPrecision
+    val matrix = new Array[ConcBuffer[Body]](sectorPrecision * sectorPrecision)
+    for (i <- 0 until matrix.length) matrix(i) = new ConcBuffer
+
+    def +=(b: Body): SectorMatrix = {
+      val x = Math.max(Math.min(b.x, boundaries.maxX), boundaries.minX)
+      val y = Math.max(Math.min(b.y, boundaries.maxY), boundaries.minY)
+
+      //Distance from top-left corner
+      val dx = x - boundaries.minX;
+      val dy = y - boundaries.minY;
+
+      //Corresponding sector
+      val xSect = (dx / sectorSize).toInt
+      val ySect = (dy / sectorSize).toInt
+
+      apply(xSect, ySect) += b
+
+      this
+    }
+
+    def apply(x: Int, y: Int) = matrix(y * sectorPrecision + x)
+
+    def combine(that: SectorMatrix): SectorMatrix = {
+      var i = 0;
+      while(i < matrix.length){
+        matrix(i) = matrix(i).combine(that.matrix(i));
+        i += 1
+      }
+      this
+    }
+
+    def toQuad(parallelism: Int): Quad = {
+      def BALANCING_FACTOR = 4
+      def quad(x: Int, y: Int, span: Int, achievedParallelism: Int): Quad = {
+        if (span == 1) {
+          val sectorSize = boundaries.size / sectorPrecision
+          val centerX = boundaries.minX + x * sectorSize + sectorSize / 2
+          val centerY = boundaries.minY + y * sectorSize + sectorSize / 2
+          var emptyQuad: Quad = Empty(centerX, centerY, sectorSize)
+          val sectorBodies = this(x, y)
+          sectorBodies.foldLeft(emptyQuad)(_ insert _)
+        } else {
+          val nspan = span / 2
+          val nAchievedParallelism = achievedParallelism * 4
+          val (nw, ne, sw, se) =
+            if (parallelism > 1 && achievedParallelism < parallelism * BALANCING_FACTOR) parallel(
+              quad(x, y, nspan, nAchievedParallelism),
+              quad(x + nspan, y, nspan, nAchievedParallelism),
+              quad(x, y + nspan, nspan, nAchievedParallelism),
+              quad(x + nspan, y + nspan, nspan, nAchievedParallelism)
+            ) else (
+              quad(x, y, nspan, nAchievedParallelism),
+              quad(x + nspan, y, nspan, nAchievedParallelism),
+              quad(x, y + nspan, nspan, nAchievedParallelism),
+              quad(x + nspan, y + nspan, nspan, nAchievedParallelism)
+            )
+          Fork(nw, ne, sw, se)
+        }
+      }
+
+      quad(0, 0, sectorPrecision, 1)
+    }
+
+    override def toString = s"SectorMatrix(#bodies: ${matrix.map(_.size).sum})"
+  }
+
+  class TimeStatistics {
+    private val timeMap = collection.mutable.Map[String, (Double, Int)]()
+
+    def clear() = timeMap.clear()
+
+    def timed[T](title: String)(body: =>T): T = {
+      var res: T = null.asInstanceOf[T]
+      val totalTime = /*measure*/ {
+        val startTime = System.currentTimeMillis()
+        res = body
+        (System.currentTimeMillis() - startTime)
+      }
+
+      timeMap.get(title) match {
+        case Some((total, num)) => timeMap(title) = (total + totalTime, num + 1)
+        case None => timeMap(title) = (0.0, 0)
+      }
+
+      println(s"$title: ${totalTime} ms; avg: ${timeMap(title)._1 / timeMap(title)._2}")
+      res
+    }
+
+    override def toString = {
+      timeMap map {
+        case (k, (total, num)) => k + ": " + (total / num * 100).toInt / 100.0 + " ms"
+      } mkString("\n")
+    }
+  }
+
+  val forkJoinPool = new ForkJoinPool
+
+  abstract class TaskScheduler {
+    def schedule[T](body: => T): ForkJoinTask[T]
+    def parallel[A, B](taskA: => A, taskB: => B): (A, B) = {
+      val right = task {
+        taskB
+      }
+      val left = taskA
+      (left, right.join())
+    }
+  }
+
+  class DefaultTaskScheduler extends TaskScheduler {
+    def schedule[T](body: => T): ForkJoinTask[T] = {
+      val t = new RecursiveTask[T] {
+        def compute = body
+      }
+      Thread.currentThread match {
+        case wt: ForkJoinWorkerThread =>
+          t.fork()
+        case _ =>
+          forkJoinPool.execute(t)
+      }
+      t
+    }
+  }
+
+  val scheduler =
+    new DynamicVariable[TaskScheduler](new DefaultTaskScheduler)
+
+  def task[T](body: => T): ForkJoinTask[T] = {
+    scheduler.value.schedule(body)
+  }
+
+  def parallel[A, B](taskA: => A, taskB: => B): (A, B) = {
+    scheduler.value.parallel(taskA, taskB)
+  }
+
+  def parallel[A, B, C, D](taskA: => A, taskB: => B, taskC: => C, taskD: => D): (A, B, C, D) = {
+    val ta = task { taskA }
+    val tb = task { taskB }
+    val tc = task { taskC }
+    val td = taskD
+    (ta.join(), tb.join(), tc.join(), td)
+  }
+}

src/test/scala/actorbintree/BinaryTreeSuite.scala

@@ -1,126 +0,0 @@
-/**
- * Copyright (C) 2009-2015 Typesafe Inc. <http://www.typesafe.com>
- */
-package actorbintree
-
-import akka.actor.{ActorRef, ActorSystem, Props, actorRef2Scala, scala2ActorRef}
-import akka.testkit.{ImplicitSender, TestKit, TestProbe}
-import org.junit.Test
-import org.junit.Assert._
-
-import scala.util.Random
-import scala.concurrent.duration._
-
-class BinaryTreeSuite extends TestKit(ActorSystem("BinaryTreeSuite")) with ImplicitSender {
-
-  import actorbintree.BinaryTreeSet._
-
-  def receiveN(requester: TestProbe, ops: Seq[Operation], expectedReplies: Seq[OperationReply]): Unit =
-    requester.within(5.seconds) {
-      val repliesUnsorted = for (i <- 1 to ops.size) yield try {
-        requester.expectMsgType[OperationReply]
-      } catch {
-        case ex: Throwable if ops.size > 10 => sys.error(s"failure to receive confirmation $i/${ops.size}\n$ex")
-        case ex: Throwable                  => sys.error(s"failure to receive confirmation $i/${ops.size}\nRequests:" + ops.mkString("\n    ", "\n     ", "") + s"\n$ex")
-      }
-      val replies = repliesUnsorted.sortBy(_.id)
-      if (replies != expectedReplies) {
-        val pairs = (replies zip expectedReplies).zipWithIndex filter (x => x._1._1 != x._1._2)
-        fail("unexpected replies:" + pairs.map(x => s"at index ${x._2}: got ${x._1._1}, expected ${x._1._2}").mkString("\n    ", "\n    ", ""))
-      }
-    }
-
-  def verify(probe: TestProbe, ops: Seq[Operation], expected: Seq[OperationReply]): Unit = {
-    val topNode = system.actorOf(Props[BinaryTreeSet])
-
-    ops foreach { op =>
-      topNode ! op
-    }
-
-    receiveN(probe, ops, expected)
-    // the grader also verifies that enough actors are created
-  }
-
-  @Test def `proper inserts and lookups (5pts)`(): Unit = {
-    val topNode = system.actorOf(Props[BinaryTreeSet])
-
-    topNode ! Contains(testActor, id = 1, 1)
-    expectMsg(ContainsResult(1, false))
-
-    topNode ! Insert(testActor, id = 2, 1)
-    topNode ! Contains(testActor, id = 3, 1)
-
-    expectMsg(OperationFinished(2))
-    expectMsg(ContainsResult(3, true))
-    ()
-  }
-
-  @Test def `instruction example (5pts)`(): Unit = {
-    val requester = TestProbe()
-    val requesterRef = requester.ref
-    val ops = List(
-      Insert(requesterRef, id=100, 1),
-      Contains(requesterRef, id=50, 2),
-      Remove(requesterRef, id=10, 1),
-      Insert(requesterRef, id=20, 2),
-      Contains(requesterRef, id=80, 1),
-      Contains(requesterRef, id=70, 2)
-    )
-
-    val expectedReplies = List(
-      OperationFinished(id=10),
-      OperationFinished(id=20),
-      ContainsResult(id=50, false),
-      ContainsResult(id=70, true),
-      ContainsResult(id=80, false),
-      OperationFinished(id=100)
-    )
-
-    verify(requester, ops, expectedReplies)
-  }
-
-
-  @Test def `behave identically to built-in set (includes GC) (40pts)`(): Unit = {
-    val rnd = new Random()
-    def randomOperations(requester: ActorRef, count: Int): Seq[Operation] = {
-      def randomElement: Int = rnd.nextInt(100)
-      def randomOperation(requester: ActorRef, id: Int): Operation = rnd.nextInt(4) match {
-        case 0 => Insert(requester, id, randomElement)
-        case 1 => Insert(requester, id, randomElement)
-        case 2 => Contains(requester, id, randomElement)
-        case 3 => Remove(requester, id, randomElement)
-      }
-
-      for (seq <- 0 until count) yield randomOperation(requester, seq)
-    }
-
-    def referenceReplies(operations: Seq[Operation]): Seq[OperationReply] = {
-      var referenceSet = Set.empty[Int]
-      def replyFor(op: Operation): OperationReply = op match {
-        case Insert(_, seq, elem) =>
-          referenceSet = referenceSet + elem
-          OperationFinished(seq)
-        case Remove(_, seq, elem) =>
-          referenceSet = referenceSet - elem
-          OperationFinished(seq)
-        case Contains(_, seq, elem) =>
-          ContainsResult(seq, referenceSet(elem))
-      }
-
-      for (op <- operations) yield replyFor(op)
-    }
-
-    val requester = TestProbe()
-    val topNode = system.actorOf(Props[BinaryTreeSet])
-    val count = 1000
-
-    val ops = randomOperations(requester.ref, count)
-    val expectedReplies = referenceReplies(ops)
-
-    ops foreach { op =>
-      topNode ! op
-      if (rnd.nextDouble() < 0.1) topNode ! GC
-    }
-    receiveN(requester, ops, expectedReplies)
-  }
-}

src/test/scala/barneshut/BarnesHutSuite.scala

@@ -0,0 +1,138 @@
+package barneshut
+
+import java.util.concurrent._
+import scala.collection._
+import scala.math._
+import scala.collection.parallel._
+import barneshut.conctrees.ConcBuffer
+import org.junit._
+import org.junit.Assert.{assertEquals, fail}
+
+class BarnesHutSuite {
+  // test cases for quad tree
+
+import FloatOps._
+  @Test def `Empty: center of mass should be the center of the cell`: Unit = {
+    val quad = Empty(51f, 46.3f, 5f)
+    assert(quad.massX == 51f, s"${quad.massX} should be 51f")
+    assert(quad.massY == 46.3f, s"${quad.massY} should be 46.3f")
+  }
+
+  @Test def `Empty: mass should be 0`: Unit = {
+    val quad = Empty(51f, 46.3f, 5f)
+    assert(quad.mass == 0f, s"${quad.mass} should be 0f")
+  }
+
+  @Test def `Empty: total should be 0`: Unit = {
+    val quad = Empty(51f, 46.3f, 5f)
+    assert(quad.total == 0, s"${quad.total} should be 0")
+  }
+
+  @Test def `Leaf with 1 body`: Unit = {
+    val b = new Body(123f, 18f, 26f, 0f, 0f)
+    val quad = Leaf(17.5f, 27.5f, 5f, Seq(b))
+
+    assert(quad.mass ~= 123f, s"${quad.mass} should be 123f")
+    assert(quad.massX ~= 18f, s"${quad.massX} should be 18f")
+    assert(quad.massY ~= 26f, s"${quad.massY} should be 26f")
+    assert(quad.total == 1, s"${quad.total} should be 1")
+  }
+
+
+  @Test def `Fork with 3 empty quadrants and 1 leaf (nw)`: Unit = {
+    val b = new Body(123f, 18f, 26f, 0f, 0f)
+    val nw = Leaf(17.5f, 27.5f, 5f, Seq(b))
+    val ne = Empty(22.5f, 27.5f, 5f)
+    val sw = Empty(17.5f, 32.5f, 5f)
+    val se = Empty(22.5f, 32.5f, 5f)
+    val quad = Fork(nw, ne, sw, se)
+
+    assert(quad.centerX == 20f, s"${quad.centerX} should be 20f")
+    assert(quad.centerY == 30f, s"${quad.centerY} should be 30f")
+    assert(quad.mass ~= 123f, s"${quad.mass} should be 123f")
+    assert(quad.massX ~= 18f, s"${quad.massX} should be 18f")
+    assert(quad.massY ~= 26f, s"${quad.massY} should be 26f")
+    assert(quad.total == 1, s"${quad.total} should be 1")
+  }
+
+  @Test def `Empty.insert(b) should return a Leaf with only that body (2pts)`: Unit = {
+    val quad = Empty(51f, 46.3f, 5f)
+    val b = new Body(3f, 54f, 46f, 0f, 0f)
+    val inserted = quad.insert(b)
+    inserted match {
+      case Leaf(centerX, centerY, size, bodies) =>
+        assert(centerX == 51f, s"$centerX should be 51f")
+        assert(centerY == 46.3f, s"$centerY should be 46.3f")
+        assert(size == 5f, s"$size should be 5f")
+        assert(bodies == Seq(b), s"$bodies should contain only the inserted body")
+      case _ =>
+        fail("Empty.insert() should have returned a Leaf, was $inserted")
+    }
+  }
+
+  // test cases for Body
+
+  @Test def `Body.updated should do nothing for Empty quad trees`: Unit = {
+    val b1 = new Body(123f, 18f, 26f, 0f, 0f)
+    val body = b1.updated(Empty(50f, 60f, 5f))
+
+    assertEquals(0f, body.xspeed, precisionThreshold)
+    assertEquals(0f, body.yspeed, precisionThreshold)
+  }
+
+  @Test def `Body.updated should take bodies in a Leaf into account (2pts)`: Unit = {
+    val b1 = new Body(123f, 18f, 26f, 0f, 0f)
+    val b2 = new Body(524.5f, 24.5f, 25.5f, 0f, 0f)
+    val b3 = new Body(245f, 22.4f, 41f, 0f, 0f)
+
+    val quad = Leaf(15f, 30f, 20f, Seq(b2, b3))
+
+    val body = b1.updated(quad)
+
+    assert(body.xspeed ~= 12.587037f)
+    assert(body.yspeed ~= 0.015557117f)
+  }
+
+  // test cases for sector matrix
+
+  @Test def `'SectorMatrix.+=' should add a body at (25,47) to the correct bucket of a sector matrix of size 96 (2pts)`: Unit = {
+    val body = new Body(5, 25, 47, 0.1f, 0.1f)
+    val boundaries = new Boundaries()
+    boundaries.minX = 1
+    boundaries.minY = 1
+    boundaries.maxX = 97
+    boundaries.maxY = 97
+    val sm = new SectorMatrix(boundaries, SECTOR_PRECISION)
+    sm += body
+    val res = sm(2, 3).size == 1 && sm(2, 3).find(_ == body).isDefined
+    assert(res, s"Body not found in the right sector")
+  }
+
+  @Rule def individualTestTimeout = new org.junit.rules.Timeout(10 * 1000)
+}
+
+object FloatOps {
+  val precisionThreshold = 1e-4
+
+  /** Floating comparison: assert(float ~= 1.7f). */
+  implicit class FloatOps(val self: Float) extends AnyVal {
+    def ~=(that: Float): Boolean =
+      abs(self - that) < precisionThreshold
+  }
+
+  /** Long floating comparison: assert(double ~= 1.7). */
+  implicit class DoubleOps(val self: Double) extends AnyVal {
+    def ~=(that: Double): Boolean =
+      abs(self - that) < precisionThreshold
+  }
+
+  /** Floating sequences comparison: assert(floatSeq ~= Seq(0.5f, 1.7f). */
+  implicit class FloatSequenceOps(val self: Seq[Float]) extends AnyVal {
+    def ~=(that: Seq[Float]): Boolean =
+      self.size == that.size &&
+        self.zip(that).forall { case (a, b) =>
+          abs(a - b) < precisionThreshold
+        }
+  }
+}
+