diff --git a/ecdf.py b/ecdf.py
new file mode 100644
index 0000000..872ffb9
--- /dev/null
+++ b/ecdf.py
@@ -0,0 +1,137 @@
+import sys
+import csv
+import argparse
+import numpy as np
+import matplotlib.pyplot as plt
+from difflib import SequenceMatcher
+
+
+def guess_number_evals(filenames):
+    """Guess the number of evals from first file."""
+    with open(filenames[0], 'r') as fd:
+        nevals = len(fd.readlines())
+    return nevals
+
+
+def along_runtime(filenames, data):
+    for fid,filename in enumerate(filenames):
+        with open(filename, 'r') as fd:
+            strdata = csv.reader(fd, delimiter=';')
+            for i,row in enumerate(strdata):
+                evals = int(row[0])
+                val = float(row[1])
+                data[evals,fid] = val
+    return data
+
+
+def cumul(data, delta, optim = None, do_min = False):
+    # Keep only best values along columns.
+    for i in range(1,len(data)):
+        for j in range(len(data[i])):
+            data[i,j] = max( data[i,j], data[i-1,j] )
+
+    if not optim:
+        optim = data.max()
+
+    # Normalize.
+    norm = data/optim
+
+    # Threshold.
+    if do_min:
+        ecdf = (norm < delta)
+    else:
+        ecdf = (norm > delta)
+
+    # Sum across rows.
+    return ecdf.sum(axis=1)/data.shape[1]
+
+
+def parse(filenames, delta, nb_rows = None, optim = None, do_min = False):
+    if not nb_rows:
+        nb_rows = guess_number_evals(filenames)
+
+    data = np.zeros( (nb_rows+1, len(filenames)) )
+    data = along_runtime(filenames,data)
+    ert = cumul(data, delta, optim, do_min)
+    return ert
+
+
+def make_name(names, delta, erts, name_strip = [], do_min = False):
+    common = names[0]
+    for run in names:
+        match = SequenceMatcher(None, common, run).find_longest_match(0, len(common), 0, len(run))
+        common = common[match.a: match.a + match.size]
+
+    for strp in name_strip:
+        common = common.replace(strp,"")
+
+    name = u"{} $\Delta={}$".format(common,delta)
+
+    if name in erts:
+        i += 1
+        name += " ({})".format(i)
+
+    return name
+
+
+if __name__ == "__main__":
+
+
+    can = argparse.ArgumentParser()
+
+    can.add_argument("-e", "--evals", metavar="NB", default=None, type=int,
+            help="Max number of evaluations to consider")
+
+    # can.add_argument("-q", "--quality", action='store_true',
+    #         help="Produce Expected Quality ECDF, instead of Expected Runtime ECDF.")
+
+    can.add_argument("-m", "--min", action='store_true',
+            help="Minimization problem, instead of maximization.")
+
+    can.add_argument("-o", "--optimum", metavar="VAL", default=None, type=float,
+            help="Best value used for normalization (else, default to the max in the data).")
+
+    can.add_argument("-s", "--name-strip", metavar="STR", default=[],
+            type=str, action='append',
+            help="Remove this string from the labels.")
+
+    can.add_argument("-d", "--delta", metavar="PERC",
+            action='append', type=float, required=True,
+            help="Target(s), as a percentage of values normalized against optimum.")
+
+    can.add_argument("-r", "--runs", metavar="FILES", nargs='*', required=True, action='append')
+
+    the = can.parse_args()
+
+    print(the.name_strip)
+
+    erts = {}
+    names = []
+    i = 0
+    for runs in the.runs:
+        for delta in the.delta:
+            ert = parse(
+                    runs, delta,
+                    nb_rows = the.evals, optim = the.optimum, do_min = the.min
+                )
+
+            name = make_name(runs, delta, erts, the.name_strip, the.min)
+            erts[name] = ert
+
+    fig = plt.figure()
+    for name in erts:
+        plt.plot(erts[name], label=name)
+
+    plt.ylim([0,1])
+
+    if the.min:
+        comp = "<"
+    else:
+        comp=">"
+    # plt.ylabel(r"$P\left(f\left(\hat{x})\right)/"+str(the.optimum)+comp+r"\Delta\right)$")
+    plt.ylabel(r"$P\left(1/"+str(the.optimum)+r"\cdot f\left(\hat{x})\right)"+comp+r"\Delta\right)$")
+    plt.xlabel("Time (#function evals)")
+    plt.title("Expected RunTime Empirical Cumulative Density Function")
+    plt.legend()
+    plt.show()
+
diff --git a/expe.py b/expe.py
new file mode 100644
index 0000000..1f8a2d2
--- /dev/null
+++ b/expe.py
@@ -0,0 +1,46 @@
+
+if __name__ == "__main__":
+    import os
+    import subprocess
+
+    # can = argparse.ArgumentParser()
+    #
+    # can.add_argument("-n", "--nb-sensors", metavar="NB", default=3, type=int,
+    #         help="Number of sensors")
+    #
+    # can.add_argument("-r", "--sensor-range", metavar="RATIO", default=0.3, type=float,
+    #         help="Sensors' range (as a fraction of domain width)")
+    #
+    # can.add_argument("-w", "--domain-width", metavar="NB", default=30, type=int,
+    #         help="Domain width (a number of cells)")
+    #
+    # can.add_argument("-i", "--iters", metavar="NB", default=100, type=int,
+    #         help="Maximum number of iterations")
+    #
+    # the = can.parse_args()
+
+    const_args=" --nb-sensors 5 --sensor-range 0.2 --domain-width 50 --iters 10000"
+    solvers = ["num_greedy","bit_greedy","num_rand","bit_rand"]
+    nbruns = 100
+    outdir = "results"
+
+    if not os.path.exists(outdir):
+        os.mkdir(outdir)
+
+    for seed in range(nbruns):
+        procs = []
+        for solver in solvers:
+            print(seed,solver)
+            p = subprocess.Popen(
+                    "python3 snp.py "
+                    + const_args
+                    + " --no-plot --dir {} --seed {} --solver {}"
+                    .format(outdir,seed,solver),
+                    shell=True
+                )
+            procs.append(p)
+
+        for proc in procs:
+            proc.wait()
+
+
diff --git a/sho/func.py b/sho/func.py
new file mode 100644
index 0000000..be713ac
--- /dev/null
+++ b/sho/func.py
@@ -0,0 +1,26 @@
+
+########################################################################
+# Wrappers around objective functions
+########################################################################
+
+class Dump:
+    """A wrapper around an objective function that
+    dumps a line in a file every time the objective function is called."""
+
+    def __init__(self, func, filename="run.csv", fmt="{it} ; {val} ; {sol}\n", sepsol=" , "):
+        self.func = func
+        self.filename = filename
+        self.fmt = fmt
+        self.sepsol = sepsol
+        self.counter = 0
+        # Erase previous file.
+        with open(self.filename, 'w') as fd:
+            fd.write("")
+
+    def __call__(self, sol):
+        val = self.func(sol)
+        self.counter += 1
+        with open(self.filename, 'a') as fd:
+            fmtsol = self.sepsol.join([str(i) for i in sol])
+            fd.write( self.fmt.format(it=self.counter, val=val, sol=fmtsol) )
+        return val
diff --git a/snp.py b/snp.py
index ef55f00..ae97448 100644
--- a/snp.py
+++ b/snp.py
@@ -1,7 +1,8 @@
+import os
 import numpy as np
 import matplotlib.pyplot as plt
 
-from sho import *
+from sho import algo, bit, func, iters, make, num, pb, plot
 
 ########################################################################
 # Interface
@@ -30,18 +31,23 @@ if __name__=="__main__":
     can.add_argument("-s", "--seed", metavar="VAL", default=None, type=int,
             help="Random pseudo-generator seed (none for current epoch)")
 
-    solvers = ["num_greedy","bit_greedy"]
+    solvers = ["num_greedy","bit_greedy","num_rand","bit_rand"]
     can.add_argument("-m", "--solver", metavar="NAME", choices=solvers, default="num_greedy",
             help="Solver to use, among: "+", ".join(solvers))
 
-    can.add_argument("-t", "--target", metavar="VAL", default=30*30, type=float,
-            help="Objective function value target")
+    # can.add_argument("-t", "--target", metavar="VAL", default=30*30, type=float,
+    #         help="Objective function value target")
+    #
+    # can.add_argument("-y", "--steady-delta", metavar="NB", default=50, type=float,
+    #         help="Stop if no improvement after NB iterations")
+    # can.add_argument("-e", "--steady-epsilon", metavar="DVAL", default=0, type=float,
+    #         help="Stop if the improvement of the objective function value is lesser than DVAL")
 
-    can.add_argument("-y", "--steady-delta", metavar="NB", default=50, type=float,
-            help="Stop if no improvement after NB iterations")
-    can.add_argument("-e", "--steady-epsilon", metavar="DVAL", default=0, type=float,
-            help="Stop if the improvement of the objective function value is lesser than DVAL")
+    can.add_argument("-p", "--no-plot", action='store_true',
+            help="Do not display plots.")
 
+    can.add_argument("-d", "--dir", metavar="DIR", default="", type=str,
+            help="Directory to which output written files.")
 
     the = can.parse_args()
 
@@ -66,29 +72,34 @@ if __name__=="__main__":
                 agains = [
                     make.iter(iters.max,
                         nb_it = the.iters),
-                    make.iter(iters.save,
-                        filename = the.solver+".csv",
-                        fmt = "{it} ; {val} ; {sol}\n"),
+                    # make.iter(iters.save,
+                    #     filename = os.path.join(the.dir,the.solver+".csv"),
+                    #     fmt = "{it} ; {val} ; {sol}\n"),
                     make.iter(iters.log,
                         fmt="\r{it} {val}"),
                     make.iter(iters.history,
                         history = history),
-                    make.iter(iters.target,
-                        target = the.target),
-                    iters.steady(the.steady_delta, the.steady_epsilon)
+                    # make.iter(iters.target,
+                    #     target = the.target),
+                    # iters.steady(the.steady_delta, the.steady_epsilon)
                 ]
             )
 
     # Erase the previous file.
-    with open(the.solver+".csv", 'w') as fd:
-        fd.write("# {} {}\n".format(the.solver,the.domain_width))
+    # with open(the.solver+".csv", 'w') as fd:
+    #     fd.write("# {} {}\n".format(the.solver,the.domain_width))
 
     val,sol,sensors = None,None,None
     if the.solver == "num_greedy":
-        val,sol = algo.greedy(
+        fdump = func.Dump(
                 make.func(num.cover_sum,
                     domain_width = the.domain_width,
                     sensor_range = the.sensor_range * the.domain_width),
+                filename = os.path.join(the.dir,"{s}_run_{i}.csv".format(s=the.solver, i=the.seed)),
+                fmt = "{it} ; {val} ; {sol}\n"
+            )
+        val,sol = algo.greedy(
+                fdump,
                 make.init(num.rand,
                     dim = d * the.nb_sensors,
                     scale = the.domain_width),
@@ -98,11 +109,34 @@ if __name__=="__main__":
             )
         sensors = num.to_sensors(sol)
 
+    if the.solver == "num_rand":
+        fdump = func.Dump(
+                make.func(num.cover_sum,
+                    domain_width = the.domain_width,
+                    sensor_range = the.sensor_range * the.domain_width),
+                filename = os.path.join(the.dir,"{s}_run_{i}.csv".format(s=the.solver, i=the.seed)),
+                fmt = "{it} ; {val} ; {sol}\n"
+            )
+        val,sol = algo.random(
+                fdump,
+                make.init(num.rand,
+                    dim = d * the.nb_sensors,
+                    scale = the.domain_width),
+                iters
+            )
+        sensors = num.to_sensors(sol)
+
+
     elif the.solver == "bit_greedy":
-        val,sol = algo.greedy(
+        fdump = func.Dump(
                 make.func(bit.cover_sum,
                     domain_width = the.domain_width,
                     sensor_range = the.sensor_range),
+                filename = os.path.join(the.dir,"{s}_run_{i}.csv".format(s=the.solver, i=the.seed)),
+                fmt = "{it} ; {val} ; {sol}\n"
+            )
+        val,sol = algo.greedy(
+                fdump,
                 make.init(bit.rand,
                     domain_width = the.domain_width,
                     nb_sensors = the.nb_sensors),
@@ -112,30 +146,49 @@ if __name__=="__main__":
             )
         sensors = bit.to_sensors(sol)
 
+    elif the.solver == "bit_rand":
+        fdump = func.Dump(
+                make.func(bit.cover_sum,
+                    domain_width = the.domain_width,
+                    sensor_range = the.sensor_range),
+                filename = os.path.join(the.dir,"{s}_run_{i}.csv".format(s=the.solver, i=the.seed)),
+                fmt = "{it} ; {val} ; {sol}\n"
+            )
+        val,sol = algo.random(
+                fdump,
+                make.init(bit.rand,
+                    domain_width = the.domain_width,
+                    nb_sensors = the.nb_sensors),
+                iters
+            )
+        sensors = bit.to_sensors(sol)
+
+
 
     # Fancy output.
     print("\n{} : {}".format(val,sensors))
 
-    shape=(the.domain_width, the.domain_width)
+    if not the.no_plot:
+        shape=(the.domain_width, the.domain_width)
 
-    fig = plt.figure()
+        fig = plt.figure()
 
-    if the.nb_sensors ==1 and the.domain_width <= 50:
-        ax1 = fig.add_subplot(121, projection='3d')
-        ax2 = fig.add_subplot(122)
+        if the.nb_sensors ==1 and the.domain_width <= 50:
+            ax1 = fig.add_subplot(121, projection='3d')
+            ax2 = fig.add_subplot(122)
 
-        f = make.func(num.cover_sum,
-                        domain_width = the.domain_width,
-                        sensor_range = the.sensor_range * the.domain_width)
-        plot.surface(ax1, shape, f)
-        plot.path(ax1, shape, history)
-    else:
-        ax2=fig.add_subplot(111)
+            f = make.func(num.cover_sum,
+                            domain_width = the.domain_width,
+                            sensor_range = the.sensor_range * the.domain_width)
+            plot.surface(ax1, shape, f)
+            plot.path(ax1, shape, history)
+        else:
+            ax2=fig.add_subplot(111)
 
-    domain = np.zeros(shape)
-    domain = pb.coverage(domain, sensors,
-            the.sensor_range * the.domain_width)
-    domain = plot.highlight_sensors(domain, sensors)
-    ax2.imshow(domain)
+        domain = np.zeros(shape)
+        domain = pb.coverage(domain, sensors,
+                the.sensor_range * the.domain_width)
+        domain = plot.highlight_sensors(domain, sensors)
+        ax2.imshow(domain)
 
-    plt.show()
+        plt.show()