commit b06114622ff92eadb9f774f78fd753e5988df232
Author: Guillaume Martres <smarter@ubuntu.com>
Date:   Tue Feb 19 20:44:23 2019 +0100

    Add barneshut assignment

diff --git a/.gitignore b/.gitignore
new file mode 100644
index 0000000..349d2e8
--- /dev/null
+++ b/.gitignore
@@ -0,0 +1,16 @@
+# General
+*.DS_Store
+*.swp
+*~
+
+# Dotty
+*.class
+*.tasty
+*.hasTasty
+
+# sbt
+target/
+
+# Dotty IDE
+/.dotty-ide-artifact
+/.dotty-ide.json
diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
new file mode 100644
index 0000000..0ded59d
--- /dev/null
+++ b/.gitlab-ci.yml
@@ -0,0 +1,36 @@
+# DO NOT EDIT THIS FILE
+
+stages:
+  - build
+  - grade
+
+compile:
+  stage: build
+  image: lampepfl/moocs:dotty-2020-02-12
+  except:
+    - tags
+  tags:
+    - cs206
+  script:
+    - sbt packageSubmission
+  artifacts:
+    expire_in: 1 day
+    paths:
+      - submission.jar
+
+grade:
+  stage: grade
+  except:
+    - tags
+  tags:
+    - cs206
+  image:
+    name: smarter3/moocs:parprog1-barneshut-2020-03-09
+    entrypoint: [""]
+  allow_failure: true
+  before_script:
+    - mkdir -p /shared/submission/
+    - cp submission.jar /shared/submission/submission.jar
+  script:
+    - cd /grader
+    - /grader/grade | /grader/feedback-printer
diff --git a/.vscode/settings.json b/.vscode/settings.json
new file mode 100644
index 0000000..a35362b
--- /dev/null
+++ b/.vscode/settings.json
@@ -0,0 +1,8 @@
+{
+  "dotty": {
+    "trace": {
+      "remoteTracingUrl": "wss://lamppc36.epfl.ch/dotty-remote-tracer/upload/lsp.log",
+      "server": { "format": "JSON", "verbosity": "verbose" }
+    }
+  }
+}
diff --git a/assignment.sbt b/assignment.sbt
new file mode 100644
index 0000000..f7d3ae1
--- /dev/null
+++ b/assignment.sbt
@@ -0,0 +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"
+)
diff --git a/build.sbt b/build.sbt
new file mode 100644
index 0000000..21e5745
--- /dev/null
+++ b/build.sbt
@@ -0,0 +1,13 @@
+course := "parprog1"
+assignment := "barneshut"
+
+scalaVersion := "0.23.0-bin-20200211-5b006fb-NIGHTLY"
+scalacOptions ++= Seq("-language:implicitConversions", "-deprecation")
+libraryDependencies ++= Seq(
+  "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))
+
+testOptions in Test += Tests.Argument(TestFrameworks.JUnit, "-a", "-v", "-s")
+testSuite := "barneshut.BarnesHutSuite"
diff --git a/grading-tests.jar b/grading-tests.jar
new file mode 100644
index 0000000..b13b869
Binary files /dev/null and b/grading-tests.jar differ
diff --git a/project/MOOCSettings.scala b/project/MOOCSettings.scala
new file mode 100644
index 0000000..171244f
--- /dev/null
+++ b/project/MOOCSettings.scala
@@ -0,0 +1,46 @@
+package ch.epfl.lamp
+
+import sbt._
+import sbt.Keys._
+
+/**
+ * Coursera uses two versions of each assignment. They both have the same assignment key and part id but have
+ * different item ids.
+ *
+ * @param key Assignment key
+ * @param partId Assignment partId
+ * @param itemId Item id of the non premium version
+ * @param premiumItemId Item id of the premium version (`None` if the assignment is optional)
+ */
+case class CourseraId(key: String, partId: String, itemId: String, premiumItemId: Option[String])
+
+/**
+  * Settings shared by all assignments, reused in various tasks.
+  */
+object MOOCSettings extends AutoPlugin {
+
+  object autoImport {
+    val course = SettingKey[String]("course")
+    val assignment = SettingKey[String]("assignment")
+    val options = SettingKey[Map[String, Map[String, String]]]("options")
+    val courseraId = settingKey[CourseraId]("Coursera-specific information identifying the assignment")
+    val testSuite = settingKey[String]("Fully qualified name of the test suite of this assignment")
+    // Convenient alias
+    type CourseraId = ch.epfl.lamp.CourseraId
+    val CourseraId = ch.epfl.lamp.CourseraId
+  }
+
+  import autoImport._
+
+  override val globalSettings: Seq[Def.Setting[_]] = Seq(
+    // supershell is verbose, buggy and useless.
+    useSuperShell := false
+  )
+
+  override val projectSettings: Seq[Def.Setting[_]] = Seq(
+    parallelExecution in Test := false,
+    // Report test result after each test instead of waiting for every test to finish
+    logBuffered in Test := false,
+    name := s"${course.value}-${assignment.value}"
+  )
+}
diff --git a/project/StudentTasks.scala b/project/StudentTasks.scala
new file mode 100644
index 0000000..7604830
--- /dev/null
+++ b/project/StudentTasks.scala
@@ -0,0 +1,318 @@
+package ch.epfl.lamp
+
+import sbt._
+import Keys._
+
+// import scalaj.http._
+import java.io.{File, FileInputStream, IOException}
+import org.apache.commons.codec.binary.Base64
+// import play.api.libs.json.{Json, JsObject, JsPath}
+import scala.util.{Failure, Success, Try}
+
+/**
+  * Provides tasks for submitting the assignment
+  */
+object StudentTasks extends AutoPlugin {
+
+  override def requires = super.requires && MOOCSettings
+
+  object autoImport {
+    val packageSourcesOnly = TaskKey[File]("packageSourcesOnly", "Package the sources of the project")
+    val packageBinWithoutResources = TaskKey[File]("packageBinWithoutResources", "Like packageBin, but without the resources")
+    val packageSubmissionZip = TaskKey[File]("packageSubmissionZip")
+    val packageSubmission = inputKey[Unit]("package solution as an archive file")
+    val runGradingTests = taskKey[Unit]("run black-box tests used for final grading")
+  }
+
+
+  import autoImport._
+  import MOOCSettings.autoImport._
+
+  override lazy val projectSettings = Seq(
+    packageSubmissionSetting,
+    // submitSetting,
+    runGradingTestsSettings,
+
+    fork := true,
+    connectInput in run := true,
+    outputStrategy := Some(StdoutOutput),
+  ) ++ packageSubmissionZipSettings
+
+  lazy val runGradingTestsSettings = runGradingTests := {
+    val testSuiteJar = "grading-tests.jar"
+    if (!new File(testSuiteJar).exists) {
+      throw new MessageOnlyException(s"Could not find tests JarFile: $testSuiteJar")
+    }
+
+    val classPath = s"${(Test / dependencyClasspath).value.map(_.data).mkString(File.pathSeparator)}${File.pathSeparator}$testSuiteJar"
+    val junitProcess =
+      Fork.java.fork(
+        ForkOptions(),
+        "-cp" :: classPath ::
+        "org.junit.runner.JUnitCore" ::
+        (Test / testSuite).value ::
+        Nil
+      )
+
+    // Wait for tests to complete.
+    junitProcess.exitValue()
+  }
+
+
+  /** **********************************************************
+    * SUBMITTING A SOLUTION TO COURSERA
+    */
+
+  val packageSubmissionZipSettings = Seq(
+    packageSubmissionZip := {
+      val submission = crossTarget.value / "submission.zip"
+      val sources = (packageSourcesOnly in Compile).value
+      val binaries = (packageBinWithoutResources in Compile).value
+      IO.zip(Seq(sources -> "sources.zip", binaries -> "binaries.jar"), submission)
+      submission
+    },
+    artifactClassifier in packageSourcesOnly := Some("sources"),
+    artifact in (Compile, packageBinWithoutResources) ~= (art => art.withName(art.name + "-without-resources"))
+  ) ++
+  inConfig(Compile)(
+    Defaults.packageTaskSettings(packageSourcesOnly, Defaults.sourceMappings) ++
+    Defaults.packageTaskSettings(packageBinWithoutResources, Def.task {
+      val relativePaths =
+        (unmanagedResources in Compile).value.flatMap(Path.relativeTo((unmanagedResourceDirectories in Compile).value)(_))
+      (mappings in (Compile, packageBin)).value.filterNot { case (_, path) => relativePaths.contains(path) }
+    })
+  )
+
+  val maxSubmitFileSize = {
+    val mb = 1024 * 1024
+    10 * mb
+  }
+
+  /** Check that the jar exists, isn't empty, isn't crazy big, and can be read
+    * If so, encode jar as base64 so we can send it to Coursera
+    */
+  def prepareJar(jar: File, s: TaskStreams): String = {
+    val errPrefix = "Error submitting assignment jar: "
+    val fileLength = jar.length()
+    if (!jar.exists()) {
+      s.log.error(errPrefix + "jar archive does not exist\n" + jar.getAbsolutePath)
+      failSubmit()
+    } else if (fileLength == 0L) {
+      s.log.error(errPrefix + "jar archive is empty\n" + jar.getAbsolutePath)
+      failSubmit()
+    } else if (fileLength > maxSubmitFileSize) {
+      s.log.error(errPrefix + "jar archive is too big. Allowed size: " +
+        maxSubmitFileSize + " bytes, found " + fileLength + " bytes.\n" +
+        jar.getAbsolutePath)
+      failSubmit()
+    } else {
+      val bytes = new Array[Byte](fileLength.toInt)
+      val sizeRead = try {
+        val is = new FileInputStream(jar)
+        val read = is.read(bytes)
+        is.close()
+        read
+      } catch {
+        case ex: IOException =>
+          s.log.error(errPrefix + "failed to read sources jar archive\n" + ex.toString)
+          failSubmit()
+      }
+      if (sizeRead != bytes.length) {
+        s.log.error(errPrefix + "failed to read the sources jar archive, size read: " + sizeRead)
+        failSubmit()
+      } else encodeBase64(bytes)
+    }
+  }
+
+  /** Task to package solution to a given file path */
+  lazy val packageSubmissionSetting = packageSubmission := {
+    val args: Seq[String] = Def.spaceDelimited("[path]").parsed
+    val s: TaskStreams = streams.value // for logging
+    val jar = (packageSubmissionZip in Compile).value
+
+    val base64Jar = prepareJar(jar, s)
+
+    val path = args.headOption.getOrElse((baseDirectory.value / "submission.jar").absolutePath)
+    scala.tools.nsc.io.File(path).writeAll(base64Jar)
+  }
+
+/*
+  /** Task to submit a solution to coursera */
+  val submit = inputKey[Unit]("submit solution to Coursera")
+  lazy val submitSetting = submit := {
+    // Fail if scalafix linting does not pass.
+    scalafixLinting.value
+
+    val args: Seq[String] = Def.spaceDelimited("<arg>").parsed
+    val s: TaskStreams = streams.value // for logging
+    val jar = (packageSubmissionZip in Compile).value
+
+    val assignmentDetails =
+      courseraId.?.value.getOrElse(throw new MessageOnlyException("This assignment can not be submitted to Coursera because the `courseraId` setting is undefined"))
+    val assignmentKey = assignmentDetails.key
+    val courseName =
+      course.value match {
+        case "capstone" => "scala-capstone"
+        case "bigdata"  => "scala-spark-big-data"
+        case other      => other
+      }
+
+    val partId = assignmentDetails.partId
+    val itemId = assignmentDetails.itemId
+    val premiumItemId = assignmentDetails.premiumItemId
+
+    val (email, secret) = args match {
+      case email :: secret :: Nil =>
+        (email, secret)
+      case _ =>
+        val inputErr =
+          s"""|Invalid input to `submit`. The required syntax for `submit` is:
+              |submit <email-address> <submit-token>
+              |
+              |The submit token is NOT YOUR LOGIN PASSWORD.
+              |It can be obtained from the assignment page:
+              |https://www.coursera.org/learn/$courseName/programming/$itemId
+              |${
+                premiumItemId.fold("") { id =>
+                  s"""or (for premium learners):
+                     |https://www.coursera.org/learn/$courseName/programming/$id
+                   """.stripMargin
+                }
+              }
+          """.stripMargin
+        s.log.error(inputErr)
+        failSubmit()
+    }
+
+    val base64Jar = prepareJar(jar, s)
+    val json =
+      s"""|{
+          |   "assignmentKey":"$assignmentKey",
+          |   "submitterEmail":"$email",
+          |   "secret":"$secret",
+          |   "parts":{
+          |      "$partId":{
+          |         "output":"$base64Jar"
+          |      }
+          |   }
+          |}""".stripMargin
+
+    def postSubmission[T](data: String): Try[HttpResponse[String]] = {
+      val http = Http("https://www.coursera.org/api/onDemandProgrammingScriptSubmissions.v1")
+      val hs = List(
+        ("Cache-Control", "no-cache"),
+        ("Content-Type", "application/json")
+      )
+      s.log.info("Connecting to Coursera...")
+      val response = Try(http.postData(data)
+                         .headers(hs)
+                         .option(HttpOptions.connTimeout(10000)) // scalaj default timeout is only 100ms, changing that to 10s
+                         .asString) // kick off HTTP POST
+      response
+    }
+
+    val connectMsg =
+      s"""|Attempting to submit "${assignment.value}" assignment in "$courseName" course
+          |Using:
+          |- email: $email
+          |- submit token: $secret""".stripMargin
+    s.log.info(connectMsg)
+
+    def reportCourseraResponse(response: HttpResponse[String]): Unit = {
+      val code = response.code
+      val respBody = response.body
+
+       /* Sample JSON response from Coursera
+      {
+        "message": "Invalid email or token.",
+        "details": {
+          "learnerMessage": "Invalid email or token."
+        }
+      }
+      */
+
+      // Success, Coursera responds with 2xx HTTP status code
+      if (response.is2xx) {
+        val successfulSubmitMsg =
+          s"""|Successfully connected to Coursera. (Status $code)
+              |
+                |Assignment submitted successfully!
+              |
+                |You can see how you scored by going to:
+              |https://www.coursera.org/learn/$courseName/programming/$itemId/
+              |${
+            premiumItemId.fold("") { id =>
+              s"""or (for premium learners):
+                 |https://www.coursera.org/learn/$courseName/programming/$id
+                       """.stripMargin
+            }
+          }
+              |and clicking on "My Submission".""".stripMargin
+        s.log.info(successfulSubmitMsg)
+      }
+
+      // Failure, Coursera responds with 4xx HTTP status code (client-side failure)
+      else if (response.is4xx) {
+        val result = Try(Json.parse(respBody)).toOption
+        val learnerMsg = result match {
+          case Some(resp: JsObject) =>
+            (JsPath \ "details" \ "learnerMessage").read[String].reads(resp).get
+          case Some(x) => // shouldn't happen
+            "Could not parse Coursera's response:\n" + x
+          case None =>
+            "Could not parse Coursera's response:\n" + respBody
+        }
+        val failedSubmitMsg =
+          s"""|Submission failed.
+              |There was something wrong while attempting to submit.
+              |Coursera says:
+              |$learnerMsg (Status $code)""".stripMargin
+        s.log.error(failedSubmitMsg)
+      }
+
+      // Failure, Coursera responds with 5xx HTTP status code (server-side failure)
+      else if (response.is5xx) {
+        val failedSubmitMsg =
+          s"""|Submission failed.
+              |Coursera seems to be unavailable at the moment (Status $code)
+              |Check https://status.coursera.org/ and try again in a few minutes.
+           """.stripMargin
+        s.log.error(failedSubmitMsg)
+      }
+
+      // Failure, Coursera repsonds with an unexpected status code
+      else {
+        val failedSubmitMsg =
+          s"""|Submission failed.
+              |Coursera replied with an unexpected code (Status $code)
+           """.stripMargin
+        s.log.error(failedSubmitMsg)
+      }
+    }
+
+    // kick it all off, actually make request
+    postSubmission(json) match {
+      case Success(resp) => reportCourseraResponse(resp)
+      case Failure(e) =>
+        val failedConnectMsg =
+          s"""|Connection to Coursera failed.
+              |There was something wrong while attempting to connect to Coursera.
+              |Check your internet connection.
+              |${e.toString}""".stripMargin
+        s.log.error(failedConnectMsg)
+    }
+
+   }
+*/
+
+  def failSubmit(): Nothing = {
+    sys.error("Submission failed")
+  }
+
+  /**
+    * *****************
+    * DEALING WITH JARS
+    */
+  def encodeBase64(bytes: Array[Byte]): String =
+    new String(Base64.encodeBase64(bytes))
+}
diff --git a/project/build.properties b/project/build.properties
new file mode 100644
index 0000000..a919a9b
--- /dev/null
+++ b/project/build.properties
@@ -0,0 +1 @@
+sbt.version=1.3.8
diff --git a/project/buildSettings.sbt b/project/buildSettings.sbt
new file mode 100644
index 0000000..8fac702
--- /dev/null
+++ b/project/buildSettings.sbt
@@ -0,0 +1,5 @@
+// Used for Coursera submission (StudentPlugin)
+// libraryDependencies += "org.scalaj" %% "scalaj-http" % "2.4.2"
+// libraryDependencies += "com.typesafe.play" %% "play-json" % "2.7.4"
+// Used for Base64 (StudentPlugin)
+libraryDependencies += "commons-codec" % "commons-codec" % "1.10"
diff --git a/project/plugins.sbt b/project/plugins.sbt
new file mode 100644
index 0000000..017735d
--- /dev/null
+++ b/project/plugins.sbt
@@ -0,0 +1,2 @@
+addSbtPlugin("org.scala-js" % "sbt-scalajs" % "0.6.28")
+addSbtPlugin("ch.epfl.lamp" % "sbt-dotty" % "0.4.0")
diff --git a/src/main/scala/barneshut/BarnesHut.scala b/src/main/scala/barneshut/BarnesHut.scala
new file mode 100644
index 0000000..d84f337
--- /dev/null
+++ b/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()
+  }
+
+}
diff --git a/src/main/scala/barneshut/Interfaces.scala b/src/main/scala/barneshut/Interfaces.scala
new file mode 100644
index 0000000..9d06659
--- /dev/null
+++ b/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
+}
diff --git a/src/main/scala/barneshut/SimulationCanvas.scala b/src/main/scala/barneshut/SimulationCanvas.scala
new file mode 100644
index 0000000..24f7d4c
--- /dev/null
+++ b/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()
+    }
+  })
+
+}
diff --git a/src/main/scala/barneshut/SimulationModel.scala b/src/main/scala/barneshut/SimulationModel.scala
new file mode 100644
index 0000000..eb27d56
--- /dev/null
+++ b/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
+  }
+
+}
diff --git a/src/main/scala/barneshut/Simulator.scala b/src/main/scala/barneshut/Simulator.scala
new file mode 100644
index 0000000..7fdafb2
--- /dev/null
+++ b/src/main/scala/barneshut/Simulator.scala
@@ -0,0 +1,103 @@
+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 = {
+    ???
+  }
+
+  def mergeBoundaries(a: Boundaries, b: Boundaries): Boundaries = {
+    ???
+  }
+
+  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
+    ???
+  }
+
+  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
+    ???
+  }
+
+  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)
+  }
+
+}
diff --git a/src/main/scala/barneshut/conctrees/Conc.scala b/src/main/scala/barneshut/conctrees/Conc.scala
new file mode 100644
index 0000000..6afb50e
--- /dev/null
+++ b/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)
+  }
+
+}
diff --git a/src/main/scala/barneshut/conctrees/ConcBuffer.scala b/src/main/scala/barneshut/conctrees/ConcBuffer.scala
new file mode 100644
index 0000000..37f451d
--- /dev/null
+++ b/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)
+  }
+
+}
diff --git a/src/main/scala/barneshut/conctrees/package.scala b/src/main/scala/barneshut/conctrees/package.scala
new file mode 100644
index 0000000..627d07c
--- /dev/null
+++ b/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]
+  }
+}
diff --git a/src/main/scala/barneshut/package.scala b/src/main/scala/barneshut/package.scala
new file mode 100644
index 0000000..20b05b4
--- /dev/null
+++ b/src/main/scala/barneshut/package.scala
@@ -0,0 +1,273 @@
+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 = ???
+    def massY: Float = ???
+    def mass: Float = ???
+    def total: Int = ???
+    def insert(b: Body): Quad = ???
+  }
+
+  case class Fork(
+    nw: Quad, ne: Quad, sw: Quad, se: Quad
+  ) extends Quad {
+    val centerX: Float = ???
+    val centerY: Float = ???
+    val size: Float = ???
+    val mass: Float = ???
+    val massX: Float = ???
+    val massY: Float = ???
+    val total: Int = ???
+
+    def insert(b: Body): Fork = {
+      ???
+    }
+  }
+
+  case class Leaf(centerX: Float, centerY: Float, size: Float, bodies: coll.Seq[Body])
+  extends Quad {
+    val (mass, massX, massY) = (??? : Float, ??? : Float, ??? : Float)
+    val total: Int = ???
+    def insert(b: Body): Quad = ???
+  }
+
+  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) =>
+          // add force contribution of each body by calling addForce
+        case Fork(nw, ne, sw, 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 = {
+      ???
+      this
+    }
+
+    def apply(x: Int, y: Int) = matrix(y * sectorPrecision + x)
+
+    def combine(that: SectorMatrix): SectorMatrix = {
+      ???
+    }
+
+    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)
+  }
+}
diff --git a/src/test/scala/barneshut/BarnesHutSuite.scala b/src/test/scala/barneshut/BarnesHutSuite.scala
new file mode 100644
index 0000000..bf947c1
--- /dev/null
+++ b/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
+        }
+  }
+}
+