diff --git a/sho/algo.py b/sho/algo.py
index 3ff1faf..d9e2a70 100644
--- a/sho/algo.py
+++ b/sho/algo.py
@@ -12,7 +12,7 @@ def random(func, init, again):
     while again(i, val, sol):
         sol = init()
         val = func(sol)
-        if val > best_val:
+        if val >= best_val:
             best_val = val
             best_sol = sol
         i += 1
@@ -28,7 +28,8 @@ def greedy(func, init, neighb, again):
     while again(i, best_val, best_sol):
         sol = neighb(best_sol)
         val = func(sol)
-        if val > best_val:
+        # Use >= and not >, so as to fallback to random walk on plateus.
+        if val >= best_val:
             best_val = val
             best_sol = sol
         i += 1
diff --git a/sho/iters.py b/sho/iters.py
index 7744c0d..de3c4fa 100644
--- a/sho/iters.py
+++ b/sho/iters.py
@@ -10,6 +10,7 @@ def max(i, val, sol, nb_it):
     else:
         return False
 
+
 # Stopping criterions that are actually just checkpoints.
 
 def several(i, val, sol, agains):
diff --git a/sho/pb.py b/sho/pb.py
index 89c5c45..a09e9f3 100644
--- a/sho/pb.py
+++ b/sho/pb.py
@@ -25,3 +25,45 @@ def coverage(domain, sensors, sensor_range):
     return domain
 
 
+def line(x0, y0, x1, y1):
+    """Compute the set of pixels (integer coordinates) of the line
+    between the given line (x0,y0) -> (x1,y1).
+    Use the Bresenham's algorithm.
+    This make a generator that yield the start and the end points.
+    """
+    dx = x1 - x0
+    dy = y1 - y0
+
+    if dx > 0:
+        xs = 1
+    else:
+        xs = -1
+
+    if dy > 0:
+        ys = 1
+    else:
+        xs = -1
+
+    dx = abs(dx)
+    dy = abs(dy)
+
+    if dx > dy:
+        ax, xy, yx, ay = xs, 0, 0, ys
+    else:
+        dx, dy = dy, dx
+        ax, xy, yx, ay = 0, ys, xs, 0
+
+    D = 2 * dy - dx
+    y = 0
+
+    for x in range(dx + 1):
+        yield x0 + x*ax + y*yx , y0 + x*xy + y*ay
+
+        if D >= 0:
+            y += 1
+            D -= 2 * dx
+
+    D += 2 * dy
+
+
+#TODO maximin trajectories
diff --git a/sho/plot.py b/sho/plot.py
index 251dca4..611345a 100644
--- a/sho/plot.py
+++ b/sho/plot.py
@@ -18,7 +18,7 @@ def surface(ax, shape, f):
     Z = np.zeros( shape )
     for y in range(shape[0]):
         for x in range(shape[1]):
-            Z[y][x] = f( (x,y), shape[0]/2 )
+            Z[y][x] = f( (x,y) )
 
     X = np.arange(0,shape[0],1)
     Y = np.arange(0,shape[1],1)
@@ -61,38 +61,3 @@ def highlight_sensors(domain, sensors, val=2):
         domain[y(s)][x(s)] = val
     return domain
 
-
-if __name__=="__main__":
-    import snp
-
-    w = 100
-    shape = (w,w)
-    history = []
-
-    val,sol = snp.greedy(
-            snp.make_func(sphere,
-                offset = w/2),
-            snp.make_init(snp.num_rand,
-                dim = 2 * 1,
-                scale = w),
-            snp.make_neig(snp.num_neighb_square,
-                scale = w/10),
-            snp.make_iter(
-                    snp.several,
-                    agains = [
-                        snp.make_iter(snp.iter_max,
-                            nb_it = 100),
-                        snp.make_iter(snp.history,
-                            history = history)
-                    ]
-                )
-        )
-    sensors = snp.num_to_sensors(sol)
-
-    #print("\n".join([str(i) for i in history]))
-
-    fig = plt.figure()
-    ax = fig.gca(projection='3d')
-    surface(ax, shape, sphere)
-    path(ax, shape, history)
-    plt.show()
diff --git a/snp.py b/snp.py
index 2732543..2847bea 100644
--- a/snp.py
+++ b/snp.py
@@ -1,7 +1,5 @@
-import sys
 import numpy as np
 import matplotlib.pyplot as plt
-import copy
 
 from sho import *
 
@@ -23,7 +21,7 @@ if __name__=="__main__":
     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=100, type=int,
+    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,
@@ -55,9 +53,9 @@ if __name__=="__main__":
     # Weird numpy way to ensure single line print of array.
     np.set_printoptions(linewidth = np.inf)
 
-    domain = np.zeros((the.domain_width, the.domain_width))
 
     # Common termination and checkpointing.
+    history = []
     iters = make.iter(
                 iters.several,
                 agains = [
@@ -67,7 +65,9 @@ if __name__=="__main__":
                         filename = the.solver+".csv",
                         fmt = "{it} ; {val} ; {sol}\n"),
                     make.iter(iters.log,
-                        fmt="\r{it} {val}")
+                        fmt="\r{it} {val}"),
+                    make.iter(iters.history,
+                        history = history)
                 ]
             )
 
@@ -106,11 +106,28 @@ if __name__=="__main__":
 
 
     # Fancy output.
-    print("\n",val,":",sensors)
+    print("\n{} : {}".format(val,sensors))
 
+    shape=(the.domain_width, the.domain_width)
+
+    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)
+
+        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)
-    plt.imshow(domain)
-    plt.show()
+    ax2.imshow(domain)
 
+    plt.show()