diff --git a/.gitignore b/.gitignore
index c6254490..2eb0cc02 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1,6 +1,9 @@
 *.swp
+*.iml
+.idea
 .DS_Store
 sbt
 project/project
 project/target
 target
+logs
diff --git a/project/Build.scala b/project/Build.scala
index b2423aa8..46830ad2 100644
--- a/project/Build.scala
+++ b/project/Build.scala
@@ -15,7 +15,14 @@ object ApplicationBuild extends Build {
 
     showSuccess := false,
       
-    resolvers += "sonatype-public" at "https://oss.sonatype.org/content/groups/public"
+    resolvers ++= Seq(
+      "sonatype-public" at 
+        "https://oss.sonatype.org/content/groups/public",
+      "Sonatype Snapshots" at 
+        "https://oss.sonatype.org/content/repositories/snapshots/",
+      "Sonatype Releases" at 
+        "https://oss.sonatype.org/content/repositories/releases/"
+    )
   )
 
   val coreSettings = commonSettings ++ Seq(
@@ -25,7 +32,8 @@ object ApplicationBuild extends Build {
 
     libraryDependencies ++= Seq(
       "org.scalatest" % "scalatest_2.11" % "2.2.4" % "test",
-      "com.github.scopt" %% "scopt" % "3.3.0"
+      "com.github.scopt" %% "scopt" % "3.3.0",
+      "org.scalanlp" %% "breeze" % "0.12"
     ),
 
     mainClass in (Compile, run) := Some("tps.Main"),
diff --git a/scripts/run-scalability-analysis b/scripts/run-scalability-analysis
new file mode 100755
index 00000000..0eabd6b9
--- /dev/null
+++ b/scripts/run-scalability-analysis
@@ -0,0 +1,7 @@
+#!/bin/bash
+
+sbt "project tps-core" "run-main tps.evaluation.ScalabilityAnalysis \
+  --network data/networks/directed-pin-with-resource-edges.tsv \
+  --source EGF_HUMAN \
+  --threshold 0.01
+"
diff --git a/tps-core/src/main/scala/tps/ArgHandling.scala b/tps-core/src/main/scala/tps/ArgHandling.scala
index b7ee04aa..f59df7fd 100644
--- a/tps-core/src/main/scala/tps/ArgHandling.scala
+++ b/tps-core/src/main/scala/tps/ArgHandling.scala
@@ -2,8 +2,6 @@ package tps
 
 import tps.util.LogUtils
 
-import java.io.File
-
 object ArgHandling {
   def parseOptions(args: Array[String]) = {
     val opts = new Options()
diff --git a/tps-core/src/main/scala/tps/Graphs.scala b/tps-core/src/main/scala/tps/Graphs.scala
index 0395ccc7..f9c01720 100644
--- a/tps-core/src/main/scala/tps/Graphs.scala
+++ b/tps-core/src/main/scala/tps/Graphs.scala
@@ -1,7 +1,5 @@
 package tps
 
-import tps.util.LogUtils._
-
 object Graphs {
 
   case class Vertex(id: String) extends Serializable {
@@ -33,23 +31,24 @@ object Graphs {
       assert(v1.id <= v2.id)
     }
 
+    private var neighborMap = Map[Vertex, Set[Vertex]]().withDefaultValue(
+      Set.empty[Vertex])
+    for (e <- E) {
+      neighborMap += e.v1 -> (neighborMap(e.v1) + e.v2)
+      neighborMap += e.v2 -> (neighborMap(e.v2) + e.v1)
+    }
+
     override def toString = {
       V.mkString("V = {", ", ", "}") + "\n" +
       E.mkString("E = {", ", ", "}") + "\n" +
       sources.mkString("SRC = {", ", ", "}")
     }
 
-    def neighbors(v: Vertex): Set[Vertex] = {
-      E collect {
-        case Edge(v1, v2) if v1 == v => v2
-        case Edge(v1, v2) if v2 == v => v1
-      }
-    }
+    def neighbors(v: Vertex): Set[Vertex] = neighborMap(v)
 
     def incidentEdges(v: Vertex): Set[Edge] = {
-      E filter {
-        case Edge(v1, v2) => v1 == v || v2 == v
-      }
+      val ns = neighbors(v)
+      ns.map(n => lexicographicEdge(n, v))
     }
 
     def contains(e: Edge): Boolean = {
@@ -88,20 +87,18 @@ object Graphs {
         distances += v -> 0
       }
 
-      while (!toProcess.isEmpty) {
+      while (toProcess.nonEmpty) {
         distance += 1
-        val nextStep = toProcess.flatMap{ v =>
+        val nextStep = toProcess.flatMap { v =>
           this.neighbors(v)
-          }.filter{ v =>
-            !alreadySeen.contains(v)
-          }
+        }.diff(alreadySeen)
 
-          for (v <- nextStep) {
-            distances += v -> distance
-          }
+        for (v <- nextStep) {
+          distances += v -> distance
+        }
 
-          alreadySeen ++= nextStep
-          toProcess = nextStep
+        alreadySeen ++= nextStep
+        toProcess = nextStep
       }
 
       distances
@@ -145,6 +142,9 @@ object Graphs {
     }
   }
 
+  def lexicographicEdge(v1: Vertex, v2: Vertex): Edge =
+    GraphParsing.lexicographicEdge(v1.id, v2.id)
+
   private def reverseDirection(d: EdgeDirection): EdgeDirection = d match {
     case Forward => Backward
     case Backward => Forward
diff --git a/tps-core/src/main/scala/tps/PINParser.scala b/tps-core/src/main/scala/tps/PINParser.scala
index 36b5453e..f3c2c5a1 100644
--- a/tps-core/src/main/scala/tps/PINParser.scala
+++ b/tps-core/src/main/scala/tps/PINParser.scala
@@ -1,7 +1,6 @@
 package tps
 
 import Graphs._
-import GraphParsing._
 
 object PINParser {
   def run(f: java.io.File): DirectedGraph = {
@@ -10,11 +9,13 @@ object PINParser {
       val Seq(id1, id2, weight, orientation) = tuple
       orientation match {
         case "U" => 
-          lexicographicEdge(id1, id2) -> Set[EdgeDirection](Forward, Backward)
+          GraphParsing.lexicographicEdge(id1, id2) -> Set[EdgeDirection](
+            Forward, Backward)
         case "D" => 
-          lexicographicEdge(id1, id2) -> Set[EdgeDirection](lexicographicForwardDirection(id1, id2))
+          GraphParsing.lexicographicEdge(id1, id2) -> Set[EdgeDirection](
+            GraphParsing.lexicographicForwardDirection(id1, id2))
       }
     }
-    aggregateLabels(pairs)
+    GraphParsing.aggregateLabels(pairs)
   }
 }
diff --git a/tps-core/src/main/scala/tps/ReferenceParser.scala b/tps-core/src/main/scala/tps/ReferenceParser.scala
index bf1c6af1..b85dd133 100644
--- a/tps-core/src/main/scala/tps/ReferenceParser.scala
+++ b/tps-core/src/main/scala/tps/ReferenceParser.scala
@@ -1,7 +1,6 @@
 package tps
 
 import Graphs._
-import GraphParsing._
 
 object ReferenceParser {
   def run(f: java.io.File): (SignedDirectedGraph, Map[Edge, String]) = {
@@ -12,7 +11,7 @@ object ReferenceParser {
 
     val tuples = data.tuples.collect{ case tuple if tuple.size >= 7 =>
       val Seq(src, tgt, lra, lri, rla, rli, rest @ _*) = tuple
-      val edge = lexicographicEdge(src, tgt)
+      val edge = GraphParsing.lexicographicEdge(src, tgt)
 
       assert(!evidencePerEdge.isDefinedAt(edge))
       evidencePerEdge += edge -> rest.mkString(", ")
@@ -28,11 +27,12 @@ object ReferenceParser {
       if (labelValue(rla)) originalLabels += ActiveEdge(Backward, Activating)
       if (labelValue(rli)) originalLabels += ActiveEdge(Backward, Inhibiting)
 
-      val lexicOrientedLabels = originalLabels map { l => lexicographicLabel(src, tgt, l) }
+      val lexicOrientedLabels = originalLabels map { l =>
+        GraphParsing.lexicographicLabel(src, tgt, l) }
 
       (edge, lexicOrientedLabels)
     }
 
-    (aggregateLabels(tuples), evidencePerEdge)
+    (GraphParsing.aggregateLabels(tuples), evidencePerEdge)
   }
 }
diff --git a/tps-core/src/main/scala/tps/SignedDirectedGraphParser.scala b/tps-core/src/main/scala/tps/SignedDirectedGraphParser.scala
index 9acb2361..ccc2c9bf 100644
--- a/tps-core/src/main/scala/tps/SignedDirectedGraphParser.scala
+++ b/tps-core/src/main/scala/tps/SignedDirectedGraphParser.scala
@@ -8,7 +8,7 @@ object SignedDirectedGraphParser {
     val data = new TSVSource(f, noHeaders = true).data
     val pairs = data.tuples map { tuple =>
       val Seq(id1, tpe, id2) = tuple
-      val edge = lexicographicEdge(id1, id2)
+      val edge = GraphParsing.lexicographicEdge(id1, id2)
       tpe match {
         case "N" => edge -> Set(
           lexicographicActivation(id1, id2), 
diff --git a/tps-core/src/main/scala/tps/UndirectedGraphOps.scala b/tps-core/src/main/scala/tps/UndirectedGraphOps.scala
index 1a0eff9c..30c777ba 100644
--- a/tps-core/src/main/scala/tps/UndirectedGraphOps.scala
+++ b/tps-core/src/main/scala/tps/UndirectedGraphOps.scala
@@ -26,4 +26,13 @@ object UndirectedGraphOps {
   }
 
   def emptyGraph: UndirectedGraph = UndirectedGraph(Set.empty, Set.empty, Set.empty)
+
+  def fromDirectedGraph(g: DirectedGraph): UndirectedGraph = {
+    var V: Set[Vertex] = Set.empty
+    val E = g.keySet
+    for (e <- E) {
+      V ++= Set(e.v1, e.v2)
+    }
+    UndirectedGraph(V, E, Set.empty)
+  }
 }
diff --git a/tps-core/src/main/scala/tps/evaluation/GraphStats.scala b/tps-core/src/main/scala/tps/evaluation/GraphStats.scala
new file mode 100644
index 00000000..9c1a011c
--- /dev/null
+++ b/tps-core/src/main/scala/tps/evaluation/GraphStats.scala
@@ -0,0 +1,74 @@
+package tps.evaluation
+
+import java.io.File
+
+import tps.Graphs.UndirectedGraph
+import tps.synthesis.Synthesis
+import tps.{TimeSeries, UndirectedGraphParser}
+import tps.util.MathUtils
+
+/**
+  * Computes basic graph statistics
+  */
+object GraphStats {
+  def computeGraphStats(g: UndirectedGraph): Unit = {
+    // # nodes, # edges, avg degree, median degree
+    println(s"# vertices: ${g.V.size}")
+    println(s"# edges   : ${g.E.size}")
+
+    val vertexDegrees = g.V.map(v => g.neighbors(v).size.toDouble)
+    println(s"Average vertex degree: ${MathUtils.mean(vertexDegrees)}")
+    println(s"Median vertex degree : ${MathUtils.median(vertexDegrees)}")
+  }
+
+  def computeDataCoverageStats(g: UndirectedGraph, ts: TimeSeries): Unit = {
+    val profileIds = ts.profiles.map(_.id).toSet
+    val verticesWithData = g.V.map(_.id).intersect(profileIds)
+    val coverageRatio = verticesWithData.size.toDouble / g.V.size
+    println(s"Ratio of vertices with data: ${coverageRatio}")
+  }
+
+  /**
+    * Prints the mean number of phosphosites per protein in the given graph.
+    *
+    * The graph and time series data are not expanded.
+    */
+  def printMeanNbPhosphosites(
+    g: UndirectedGraph, ts: TimeSeries, ppm: Map[String, Set[String]]
+  ): Unit = {
+    val phosphositeCardinalities = g.V.toSeq map { v =>
+      // get peptides that map to the protein
+      val matchingPeptides = ppm.filter{ case (pep, prots) =>
+        prots contains v .id }.keySet
+      val matchingProfiles = ts.profiles.filter{ p =>
+        matchingPeptides contains p.id
+      }
+      matchingProfiles.size
+    }
+
+    val mean = MathUtils.mean(phosphositeCardinalities.map(_.toDouble))
+    println(s"Mean number of phosphosites: $mean")
+  }
+
+  def computeProfileStats(
+    ts: TimeSeries,
+    firstScores: Map[String, Seq[Double]],
+    prevScores: Map[String, Seq[Double]],
+    threshold: Double
+  ): Unit
+  = {
+    val nbSigMeasurements = ts.profiles map { p =>
+      Synthesis.nbSignificantMeasurements(p, firstScores, prevScores, threshold)
+    }
+    val med = MathUtils.median(nbSigMeasurements.map(_.toDouble))
+    println(s"Median number of significant time points: $med")
+    val mean = MathUtils.mean(nbSigMeasurements.map(_.toDouble))
+    println(s"Mean number of significant time points: $mean")
+  }
+
+  def main(args: Array[String]): Unit = {
+    val graphName = args(0)
+    val graph = UndirectedGraphParser.run(new File(graphName))
+    computeGraphStats(graph)
+  }
+}
diff --git a/tps-core/src/main/scala/tps/evaluation/ScalabilityAnalysis.scala b/tps-core/src/main/scala/tps/evaluation/ScalabilityAnalysis.scala
new file mode 100644
index 00000000..2d349b14
--- /dev/null
+++ b/tps-core/src/main/scala/tps/evaluation/ScalabilityAnalysis.scala
@@ -0,0 +1,50 @@
+package tps.evaluation
+
+import tps._
+import tps.simulation.{RandomGraphGenerator, RandomTimeSeriesGenerator}
+import tps.synthesis.{Synthesis, SynthesisOptions}
+import tps.util.{Stopwatch, TimingUtil}
+
+/**
+  * Evaluates solver running time on randomly simulated data.
+  */
+object ScalabilityAnalysis {
+
+  private val MIN_GRAPH_SIZE = 1000
+  private val MAX_GRAPH_SIZE = 128000
+
+  private val SIGNIFICANCE_THRESHOLD = 0.01
+
+  def main(args: Array[String]): Unit = {
+    val resultReporter = new NoopReporter()
+
+    var size = MIN_GRAPH_SIZE
+    do {
+      println(s"Evaluating with size $size")
+
+      TimingUtil.timeReplicates(s"Scalability analysis for $size", 3) {
+        val g = RandomGraphGenerator.generateRandomGraph(size)
+        val ppm = RandomTimeSeriesGenerator.generateRandomPeptideProteinMap(g)
+        val ts = RandomTimeSeriesGenerator.generateRandomTimeSeries(ppm.keySet)
+        val scores = RandomTimeSeriesGenerator.generateSignificanceScores(ts)
+
+        Synthesis.run(
+          g,
+          ts,
+          scores,
+          scores,
+          Set.empty,
+          ppm,
+          g.sources map (_.id),
+          SIGNIFICANCE_THRESHOLD,
+          SynthesisOptions(),
+          resultReporter
+        )
+
+        println("Nb. phosphosites: " + ppm.keySet.size)
+      }
+
+      size = size * 2
+    } while (size <= MAX_GRAPH_SIZE)
+  }
+}
diff --git a/tps-core/src/main/scala/tps/simulation/RandomGraphGenerator.scala b/tps-core/src/main/scala/tps/simulation/RandomGraphGenerator.scala
new file mode 100644
index 00000000..965a4092
--- /dev/null
+++ b/tps-core/src/main/scala/tps/simulation/RandomGraphGenerator.scala
@@ -0,0 +1,44 @@
+package tps.simulation
+
+import tps.Graphs
+import tps.Graphs.{Edge, UndirectedGraph, Vertex}
+
+/**
+  * Generates a random graph of given size, according to a predetermined
+  * likelihood of choosing between adding a new node and adding a new edge to
+  * it.
+  */
+object RandomGraphGenerator {
+
+  private val RANDOM_SEED = 131161511
+  private val random = new scala.util.Random(RANDOM_SEED)
+
+  private val NODE_CREATION_PROBABILITY = 0.5
+
+  def generateRandomGraph(nbEdges: Int): UndirectedGraph = {
+    val src = Vertex("src")
+    var V = Set[Vertex](src)
+    var E = Set[Edge]()
+
+    var i = 0
+    while (E.size < nbEdges) {
+      val srcVertex = V.toIndexedSeq(random.nextInt(V.size))
+      val tgtVertex = if (random.nextDouble() < NODE_CREATION_PROBABILITY) {
+        // add a new node
+        i += 1
+        Vertex(s"v_$i")
+      } else {
+        // add an edge to an existing node
+        V.toIndexedSeq(random.nextInt(V.size))
+      }
+      // do not add self edges
+      if (srcVertex != tgtVertex) {
+        V += tgtVertex
+        E += Graphs.lexicographicEdge(srcVertex, tgtVertex)
+      }
+    }
+
+    UndirectedGraph(V, E, Set(src))
+  }
+
+}
diff --git a/tps-core/src/main/scala/tps/simulation/RandomTimeSeriesGenerator.scala b/tps-core/src/main/scala/tps/simulation/RandomTimeSeriesGenerator.scala
new file mode 100644
index 00000000..db2993af
--- /dev/null
+++ b/tps-core/src/main/scala/tps/simulation/RandomTimeSeriesGenerator.scala
@@ -0,0 +1,82 @@
+package tps.simulation
+
+import breeze.stats.distributions.Poisson
+import tps.Graphs.UndirectedGraph
+import tps.PeptideExpansion.PeptideProteinMap
+import tps.{Profile, TimeSeries}
+
+/**
+  * Generates random time series data for a given [[tps.Graphs.UndirectedGraph]]
+  */
+object RandomTimeSeriesGenerator {
+
+  private val RANDOM_SEED = 2483967
+  private val random = new scala.util.Random(RANDOM_SEED)
+
+  private val MAX_PROFILE_VALUE = 10
+
+  private val NB_TIME_POINTS = 8
+
+  // Poisson parameter for drawing number of phosphosites
+  private val NB_SITES_POISSON_PARAMETER = 1.7
+  private val nbSitesDistribution = Poisson.distribution(
+    NB_SITES_POISSON_PARAMETER)
+
+  private val NB_SIGNIFICANT_TIME_POINTS_POISSON_PARAMETER = 0.55
+  private val nbSignificantTimePointsDistribution = Poisson.distribution(
+    NB_SIGNIFICANT_TIME_POINTS_POISSON_PARAMETER)
+
+  def generateRandomPeptideProteinMap(
+    graph: UndirectedGraph
+  ): PeptideProteinMap = {
+    var mapping: PeptideProteinMap = Map.empty
+    for (v <- graph.V) {
+      val nbSites = nbSitesDistribution.draw()
+      for (i <- 0 until nbSites) {
+        val siteName = s"${v.id}#site${i}"
+        assert(!mapping.isDefinedAt(siteName))
+        mapping += siteName -> Set(v.id)
+      }
+    }
+    mapping
+  }
+
+  def generateRandomTimeSeries(pepIDs: Set[String]): TimeSeries = {
+    val profiles = pepIDs.toSeq map { id =>
+      generateProfile(id)
+    }
+    TimeSeries(generateLabels(), profiles)
+  }
+
+  /**
+    * Generates random significance scores for the given time series data.
+    */
+  def generateSignificanceScores(
+    timeSeries: TimeSeries
+  ): Map[String, Seq[Double]] = {
+    val pairs = for (p <- timeSeries.profiles) yield {
+      // decide how many points should be significant
+      val nbSignificantTimePoints = nbSignificantTimePointsDistribution.draw()
+      val significanceProb = nbSignificantTimePoints.toDouble / NB_TIME_POINTS
+
+      // compute scores for all time points except the first
+      val sigScores = p.values.tail.map { v =>
+        if (random.nextDouble() < significanceProb) 0.0 else 1.0
+      }
+      p.id -> sigScores
+    }
+    pairs.toMap
+  }
+
+  private def generateLabels(): Seq[String] = {
+    (1 to NB_TIME_POINTS).map(i => s"t$i")
+  }
+
+  private def generateProfile(id: String): Profile = {
+    val values = (0 until NB_TIME_POINTS) map { i =>
+      Some(random.nextDouble() * MAX_PROFILE_VALUE)
+    }
+    Profile(id, values)
+  }
+
+}
diff --git a/tps-core/src/main/scala/tps/synthesis/Synthesis.scala b/tps-core/src/main/scala/tps/synthesis/Synthesis.scala
index 3f1de4d8..d8bd1dbf 100644
--- a/tps-core/src/main/scala/tps/synthesis/Synthesis.scala
+++ b/tps-core/src/main/scala/tps/synthesis/Synthesis.scala
@@ -3,6 +3,8 @@ package tps.synthesis
 import tps.Graphs._
 import tps.SignedDirectedGraphOps._
 import tps._
+import tps.evaluation.GraphStats
+import tps.util.TimingUtil
 
 object Synthesis {
   def run(
@@ -71,26 +73,44 @@ object Synthesis {
       )
     }
 
+    // debug and stats printing
     printCollapsedInterpretation()
+    println("Non-expanded graph stats:")
+    println("Unfiltered data:")
+    GraphStats.printMeanNbPhosphosites(networkWithSources, timeSeries,
+      peptideProteinMap)
+    GraphStats.computeProfileStats(timeSeries, firstScores, prevScores,
+      significanceThreshold)
+    println("Filtered, significant data:")
+    GraphStats.printMeanNbPhosphosites(networkWithSources,
+      significantTimeSeries, ppmWithData)
+    GraphStats.computeProfileStats(significantTimeSeries, firstScores,
+      prevScores, significanceThreshold)
+    println("Ratio of significant profiles: " +
+      (significantTimeSeries.profiles.size.toDouble / timeSeries.profiles.size))
+
+    println("Expanded graph stats:")
+    GraphStats.computeGraphStats(expandedNetwork)
+    GraphStats.computeDataCoverageStats(expandedNetwork, expandedTimeSeries)
 
     // dispatch solver
     val solver = opts.solver match {
-      case "dataflow" => 
+      case "dataflow" =>
         new DataflowSolver(
-          expandedNetwork, 
-          expandedPartialModel, 
-          opts, 
+          expandedNetwork,
+          expandedPartialModel,
+          opts,
           interpretation,
           resultReporter
         )
-      case "naive" => 
+      case "naive" =>
         new NaiveSymbolicSolver(
           expandedNetwork,
           expandedPartialModel,
           opts,
           interpretation
         )
-      case "bilateral" => 
+      case "bilateral" =>
         new BilateralSolver(
           expandedNetwork,
           expandedPartialModel,
@@ -99,7 +119,9 @@ object Synthesis {
         )
     }
 
-    val expandedSol = solver.summary()
+    val expandedSol = TimingUtil.time("solver") {
+      solver.summary()
+    }
     collapseSolution(expandedSol, ppmWithData)
   }
 
@@ -109,13 +131,23 @@ object Synthesis {
     prevScores: Map[String, Seq[Double]],
     threshold: Double
   ): Boolean = {
+    nbSignificantMeasurements(p, firstScores, prevScores, threshold) > 0
+  }
+
+  // TODO move
+  def nbSignificantMeasurements(
+    p: Profile,
+    firstScores: Map[String, Seq[Double]],
+    prevScores: Map[String, Seq[Double]],
+    threshold: Double
+  ): Int = {
     val fs = firstScores(p.id)
     val ps = prevScores(p.id)
     val toEval = p.values.tail
     val filteredValues = for ((v, (f, p)) <- toEval zip (fs zip ps)) yield {
       if (f < threshold || p < threshold) v else None
     }
-    filteredValues.exists(_.isDefined)
+    filteredValues.filter(_.isDefined).size
   }
 
 }
diff --git a/tps-core/src/main/scala/tps/util/MathUtils.scala b/tps-core/src/main/scala/tps/util/MathUtils.scala
index ddb9ec56..4909ed8c 100644
--- a/tps-core/src/main/scala/tps/util/MathUtils.scala
+++ b/tps-core/src/main/scala/tps/util/MathUtils.scala
@@ -8,6 +8,22 @@ object MathUtils {
   def max(xs: Seq[Double]): Double = xs.reduceLeft(Math.max)
   def min(xs: Seq[Double]): Double = xs.reduceLeft(Math.min)
 
+  def mean(vs: Iterable[Double]): Double = {
+    assert(vs.nonEmpty)
+    vs.sum / vs.size
+  }
+
+  def median(vs: Iterable[Double]): Double = {
+    assert(vs.nonEmpty)
+    val med = vs.size / 2
+    val sorted = vs.toIndexedSeq.sorted
+    if (vs.size % 2 == 0) {
+      (sorted(med - 1) + sorted(med)) / 2.0
+    } else {
+      sorted(med)
+    }
+  }
+
   def log2(x: Double) = scala.math.log(x) / scala.math.log(2)
 
   def foldChanges(vs: Seq[Double]): Seq[Double] = {
diff --git a/tps-core/src/main/scala/tps/util/TimingUtil.scala b/tps-core/src/main/scala/tps/util/TimingUtil.scala
new file mode 100644
index 00000000..be423ab6
--- /dev/null
+++ b/tps-core/src/main/scala/tps/util/TimingUtil.scala
@@ -0,0 +1,34 @@
+package tps.util
+
+/**
+  * Provides utility functions to measure running time.
+  */
+object TimingUtil {
+  def time[R](label: String)(block: => R): R = {
+    val (result, time) = computeAndTime(block _)
+    println(s"[${label}] Elapsed time: ${time}s")
+    result
+  }
+
+  def timeReplicates[R](label: String, replicates: Int)(block: => R): Unit = {
+    var times = Set[Double]()
+    for (i <- 1 to replicates) {
+      val (_, time) = computeAndTime(block _)
+      times += time
+    }
+    val min = times.min
+    val max = times.max
+    val mean = MathUtils.mean(times)
+    val median = MathUtils.median(times)
+    println(s"[$label] Elapsed: min = $min, max = $max, mean = $mean, median " +
+      s"= $median")
+  }
+
+  private def computeAndTime[R](block: () => R): (R, Double) = {
+    val t0 = System.nanoTime()
+    val result = block()    // call-by-name
+    val t1 = System.nanoTime()
+    val s = (t1 - t0) / 1000000000.0
+    (result, s)
+  }
+}
diff --git a/tps-core/src/test/scala/tps/SolverOutputTest.scala b/tps-core/src/test/scala/tps/SolverOutputTest.scala
index 2a2cd13d..5bef0394 100644
--- a/tps-core/src/test/scala/tps/SolverOutputTest.scala
+++ b/tps-core/src/test/scala/tps/SolverOutputTest.scala
@@ -1,5 +1,4 @@
 package tps
-package test
 
 import org.scalatest.FunSuite
 import org.scalatest.Matchers
@@ -7,16 +6,10 @@ import org.scalatest.Matchers
 import Graphs._
 import tps.synthesis.SynthesisOptions
 
-import java.io.File
+import TestResourceUtil.testFile
 
 class SolverOutputTest extends FunSuite with Matchers {
 
-  private val resourceBaseFolder = "/simple-example"
-  private def testFile(name: String): File = {
-    val url = getClass.getResource(s"$resourceBaseFolder/$name")
-    new File(url.getFile())
-  }
-
   test("output of solvers match") {
     val network = UndirectedGraphParser.run(testFile("network.tsv"))
     val timeSeries = TimeSeriesParser.run(testFile("time-series.tsv"))
diff --git a/tps-core/src/test/scala/tps/TestResourceUtil.scala b/tps-core/src/test/scala/tps/TestResourceUtil.scala
new file mode 100644
index 00000000..8f0559de
--- /dev/null
+++ b/tps-core/src/test/scala/tps/TestResourceUtil.scala
@@ -0,0 +1,12 @@
+package tps
+
+import java.io.File
+
+object TestResourceUtil {
+  private val RESOURCE_BASE_FOLDER = "/simple-example"
+
+  def testFile(name: String): File = {
+    val url = getClass.getResource(s"$RESOURCE_BASE_FOLDER/$name")
+    new File(url.getFile())
+  }
+}
diff --git a/tps-core/src/test/scala/tps/simulation/RandomGraphGeneratorTest.scala b/tps-core/src/test/scala/tps/simulation/RandomGraphGeneratorTest.scala
new file mode 100644
index 00000000..0b54ed3a
--- /dev/null
+++ b/tps-core/src/test/scala/tps/simulation/RandomGraphGeneratorTest.scala
@@ -0,0 +1,21 @@
+package tps.simulation
+
+import org.scalatest.FunSuite
+import org.scalatest.Matchers
+import tps.Graphs.Vertex
+import tps.UndirectedGraphParser
+import tps.TestResourceUtil.testFile
+
+class RandomGraphGeneratorTest extends FunSuite with Matchers {
+
+  test("completely random graph has one source and right number of edges") {
+    val generated = RandomGraphGenerator.generateRandomGraph(5)
+
+    // Check that there is only one source, and it is in the vertex set.
+    generated.sources.size should equal (1)
+    generated.V should contain (generated.sources.head)
+
+    generated.E.size should equal (5)
+  }
+
+}
diff --git a/tps-core/src/test/scala/tps/util/MathUtilsTest.scala b/tps-core/src/test/scala/tps/util/MathUtilsTest.scala
index 9954b55a..7d1d346e 100644
--- a/tps-core/src/test/scala/tps/util/MathUtilsTest.scala
+++ b/tps-core/src/test/scala/tps/util/MathUtilsTest.scala
@@ -27,4 +27,20 @@ class MathUtilsTest extends FunSuite with Matchers {
     val c = MathUtils.combination(list, 0)
     c should equal (List(Nil))
   }
+
+  test("mean of list") {
+    val l1 = List(1.0)
+    MathUtils.mean(l1) should equal (1.0)
+
+    val l2 = List(8.0, 2.0, 6.0, 4.0)
+    MathUtils.mean(l2) should equal (5.0)
+  }
+
+  test("median of list") {
+    val l1 = List(1.0)
+    MathUtils.median(l1) should equal (1.0)
+
+    val l2 = List(8.0, 2.0, 6.0, 4.0)
+    MathUtils.median(l2) should equal (5.0)
+  }
 }