commit a90792a323a4c4ab85dde709bff9db4a7826f59f
Author: Johann Dreo <johann.dreo@thalesgroup.com>
Date:   Mon Dec 9 11:23:26 2019 +0100

    first import

diff --git a/requirement.txt b/requirement.txt
new file mode 100644
index 0000000..aa094d9
--- /dev/null
+++ b/requirement.txt
@@ -0,0 +1,2 @@
+numpy
+matplotlib
diff --git a/sho/__init__.py b/sho/__init__.py
new file mode 100644
index 0000000..ff06111
--- /dev/null
+++ b/sho/__init__.py
@@ -0,0 +1,41 @@
+import numpy as np
+
+__all__ = [
+        'x', 'y', 'distance',
+       'algo',
+       'make',
+       'iters',
+       'num',
+       'bit',
+       'plot',
+       'pb',
+   ]
+
+########################################################################
+# Utilities
+########################################################################
+
+def x(a):
+    """Return the first element of a 2-tuple.
+    >>> x([1,2])
+    1
+    """
+    return a[0]
+
+
+def y(a):
+    """Return the second element of a 2-tuple.
+    >>> y([1,2])
+    2
+    """
+    return a[1]
+
+
+def distance(a,b):
+    """Euclidean distance (in pixels).
+
+    >>> distance( (1,1),(2,2) ) == math.sqrt(2)
+    True
+    """
+    return np.sqrt( (x(a)-x(b))**2 + (y(a)-y(b))**2 )
+
diff --git a/sho/algo.py b/sho/algo.py
new file mode 100644
index 0000000..c7188ec
--- /dev/null
+++ b/sho/algo.py
@@ -0,0 +1,41 @@
+########################################################################
+# Algorithms
+########################################################################
+import numpy as np
+
+def random(func, init, again):
+    """Iterative random search template."""
+    best_sol = init()
+    best_val = func(best_sol)
+    val,sol = best_val,best_sol
+    i = 0
+    while again(i, best_val, best_sol):
+        sol = init()
+        val = func(sol)
+        if val >= best_val:
+            best_val = val
+            best_sol = sol
+        i += 1
+    return best_val, best_sol
+
+
+def greedy(func, init, neighb, again):
+    """Iterative randomized greedy heuristic template."""
+    best_sol = init()
+    best_val = func(best_sol)
+    val,sol = best_val,best_sol
+    i = 1
+    while again(i, best_val, best_sol):
+        sol = neighb(best_sol)
+        val = func(sol)
+        # Use >= and not >, so as to avoid random walk on plateus.
+        if val >= best_val:
+            best_val = val
+            best_sol = sol
+        i += 1
+    return best_val, best_sol
+
+# TODO add a simulated-annealing template.
+
+# TODO add a population-based stochastic heuristic template.
+
diff --git a/sho/bit.py b/sho/bit.py
new file mode 100644
index 0000000..e6a38f5
--- /dev/null
+++ b/sho/bit.py
@@ -0,0 +1,68 @@
+import numpy as np
+import copy
+
+from . import x,y,pb
+
+########################################################################
+# Objective functions
+########################################################################
+
+def cover_sum(sol, domain_width, sensor_range):
+    """Compute the coverage quality of the given array of bits."""
+    domain = np.zeros((domain_width,domain_width))
+    sensors = to_sensors(sol)
+    return np.sum(pb.coverage(domain, sensors, sensor_range))
+
+
+def to_sensors(sol):
+    """Convert an square array of d lines/columns containing n ones
+    to an array of n 2-tuples with related coordinates.
+
+    >>> to_sensors([[1,0],[1,0]])
+    [(0, 0), (0, 1)]
+    """
+    sensors = []
+    for i in range(len(sol)):
+        for j in range(len(sol[i])):
+            if sol[i][j] == 1:
+                sensors.append( (j,i) )
+    return sensors
+
+
+########################################################################
+# Initialization
+########################################################################
+
+def rand(domain_width, nb_sensors):
+    """"Draw a random domain containing nb_sensors ones."""
+    domain = np.zeros( (domain_width,domain_width) )
+    for x,y in np.random.randint(0, domain_width, (nb_sensors, 2)):
+        domain[y][x] = 1
+    return domain
+
+
+########################################################################
+# Neighborhood
+########################################################################
+
+def neighb_square(sol, scale, domain_width):
+    """Draw a random array by moving ones to adjacent cells."""
+    # Copy, because Python pass by reference
+    # and we may not want to alter the original solution.
+    new = copy.copy(sol)
+    for py in range(len(sol)):
+        for px in range(len(sol[py])):
+            if sol[py][px] == 1:
+                new[py][px] = 0 # Remove original position.
+                d = np.random.randint(-scale//2,scale//2,2)
+                if py+y(d) < 0 :
+                    d[1] = np.random.randint(-py,scale//2)
+                if py+y(d) >= domain_width :
+                    d[1] = np.random.randint(-scale//2,domain_width-py)
+                if px+y(d) < 0 :
+                    d[0] = np.random.randint(0,scale//2)
+                if px+x(d) >= domain_width :
+                    d[0] = np.random.randint(-scale//2,domain_width-px)
+                new[py+y(d)][px+x(d)] = 1
+    return new
+
diff --git a/sho/iters.py b/sho/iters.py
new file mode 100644
index 0000000..33b325c
--- /dev/null
+++ b/sho/iters.py
@@ -0,0 +1,75 @@
+import sys
+from collections import deque
+
+########################################################################
+# Stopping criterions
+########################################################################
+
+def max(i, val, sol, nb_it):
+    """Stop after reaching nb_it iterations."""
+    if i < nb_it:
+        return True
+    else:
+        return False
+
+
+def target(i, val, sol, target):
+    """Stop after reaching target value."""
+    if val < target:
+        return True
+    else:
+        return False
+
+
+class steady:
+    """Stop if improvement is lesser than epsilon, in the last delta iterations."""
+
+    def __init__(self, delta, epsilon = 0):
+        self.epsilon = epsilon
+        self.delta = delta
+        self.delta_vals = deque()
+
+    def __call__(self, i, val, sol):
+        if i < self.delta: # Always wait the first delta iterations.
+            self.delta_vals.append(val)
+            return True
+
+        else:
+            #FILO stack.
+            self.delta_vals.popleft()
+            self.delta_vals.append(val)
+
+            if val - self.delta_vals[0] <= self.epsilon:
+                return False # Stop here.
+            else:
+                return True
+
+
+# Stopping criterions that are actually just checkpoints.
+
+def several(i, val, sol, agains):
+    """several stopping criterions in one."""
+    over = []
+    for again in agains:
+        over.append( again(i, val, sol) )
+    return all(over)
+
+
+def save(i, val, sol, filename="run.csv", fmt="{it} ; {val} ; {sol}\n"):
+    """Save all iterations to a file."""
+    # Append a line at the end of the file.
+    with open(filename.format(it=i), 'a') as fd:
+        fd.write( fmt.format(it=i, val=val, sol=sol) )
+    return True # No incidence on termination.
+
+
+def history(i, val, sol, history):
+    history.append((val,sol))
+    return True
+
+
+def log(i, val, sol, fmt="{it} {val}\n"):
+    """Print progress on stderr."""
+    sys.stderr.write( fmt.format(it=i, val=val) )
+    return True
+
diff --git a/sho/make.py b/sho/make.py
new file mode 100644
index 0000000..4da141b
--- /dev/null
+++ b/sho/make.py
@@ -0,0 +1,35 @@
+"""Wrappers that captures parameters of a function
+and returns an operator with a given interface."""
+
+def func(cover, **kwargs):
+    """Make an objective function from the given function.
+    An objective function takes a solution and returns a scalar."""
+    def f(sol):
+        return cover(sol,**kwargs)
+    return f
+
+
+def init(init, **kwargs):
+    """Make an initialization operator from the given function.
+    An init. op. returns a solution."""
+    def f():
+        return init(**kwargs)
+    return f
+
+
+def neig(neighb, **kwargs):
+    """Make an neighborhood operator from the given function.
+    A neighb. op. takes a solution and returns another one."""
+    def f(sol):
+        return neighb(sol, **kwargs)
+    return f
+
+
+def iter(iters, **kwargs):
+    """Make an iterations operator from the given function.
+    A iter. op. takes a value and a solution and returns
+    the current number of iterations."""
+    def f(i, val, sol):
+        return iters(i, val, sol, **kwargs)
+    return f
+
diff --git a/sho/num.py b/sho/num.py
new file mode 100644
index 0000000..c989954
--- /dev/null
+++ b/sho/num.py
@@ -0,0 +1,48 @@
+import numpy as np
+
+from . import pb
+
+########################################################################
+# Objective functions
+########################################################################
+
+# Decoupled from objective functions, so as to be used in display.
+def to_sensors(sol):
+    """Convert a vector of n*2 dimension to an array of n 2-tuples.
+
+    >>> to_sensors([0,1,2,3])
+    [(0, 1), (2, 3)]
+    """
+    sensors = []
+    for i in range(0,len(sol),2):
+        sensors.append( ( int(round(sol[i])), int(round(sol[i+1])) ) )
+    return sensors
+
+
+def cover_sum(sol, domain_width, sensor_range):
+    """Compute the coverage quality of the given vector."""
+    domain = np.zeros((domain_width,domain_width))
+    sensors = to_sensors(sol)
+    return np.sum(pb.coverage(domain, sensors, sensor_range))
+
+
+########################################################################
+# Initialization
+########################################################################
+
+def rand(dim, scale):
+    """Draw a random vector in [0,scale]**dim."""
+    return np.random.random(dim) * scale
+
+
+########################################################################
+# Neighborhood
+########################################################################
+
+def neighb_square(sol, scale, domain_width):
+    """Draw a random vector in a square of witdh `scale`
+    around the given one."""
+    # TODO handle constraints
+    new = sol + (np.random.random(len(sol)) * scale - scale/2)
+    return new
+
diff --git a/sho/pb.py b/sho/pb.py
new file mode 100644
index 0000000..3bec753
--- /dev/null
+++ b/sho/pb.py
@@ -0,0 +1,67 @@
+from . import distance
+
+########################################################################
+# Objective functions
+########################################################################
+
+def coverage(domain, sensors, sensor_range):
+    """Set a given domain's cells to on if they are visible
+    from one of the given sensors at the given sensor_range.
+
+    >>> coverage(np.zeros((5,5)),[(2,2)],2)
+    array([[ 0.,  0.,  0.,  0.,  0.],
+           [ 0.,  1.,  1.,  1.,  0.],
+           [ 0.,  1.,  1.,  1.,  0.],
+           [ 0.,  1.,  1.,  1.,  0.],
+           [ 0.,  0.,  0.,  0.,  0.]])
+    """
+    for py in range(len(domain)):
+        for px in range(len(domain[py])):
+            p = (px,py)
+            for x in sensors:
+                if distance(x,p) < sensor_range:
+                    domain[py][px] = 1
+                    break
+    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
+
diff --git a/sho/plot.py b/sho/plot.py
new file mode 100644
index 0000000..611345a
--- /dev/null
+++ b/sho/plot.py
@@ -0,0 +1,63 @@
+import numpy as np
+from matplotlib import cm
+import matplotlib.pyplot as plt
+import itertools
+from mpl_toolkits.mplot3d import Axes3D
+
+from . import x,y,distance
+
+def sphere(x,offset=0.5):
+    """Computes the square of a multi-dimensional vector x."""
+    f = 0
+    for i in range(len(x)):
+        f += (x[i]-offset)**2
+    return -1 * f
+
+
+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) )
+
+    X = np.arange(0,shape[0],1)
+    Y = np.arange(0,shape[1],1)
+    X,Y = np.meshgrid(X,Y)
+    #ax.plot_surface(X, Y, Z, cmap=cm.viridis)
+    ax.plot_surface(X, Y, Z)
+
+
+def path(ax, shape, history):
+    def pairwise(iterable):
+        a, b = itertools.tee(iterable)
+        next(b, None)
+        return zip(a, b)
+
+    k=0
+    for i,j in pairwise(range(len(history)-1)):
+        xi = history[i][1][0]
+        yi = history[i][1][1]
+        zi = history[i][0]
+        xj = history[j][1][0]
+        yj = history[j][1][1]
+        zj = history[j][0]
+        x = [xi, xj]
+        y = [yi, yj]
+        z = [zi, zj]
+        ax.plot(x,y,z, color=cm.RdYlBu(k))
+        k+=1
+
+
+def highlight_sensors(domain, sensors, val=2):
+    """Add twos to the given domain, in the cells where the given
+    sensors are located.
+
+    >>> highlight_sensors( [[0,0],[1,1]], [(0,0),(1,1)] )
+    [[2, 0], [1, 2]]
+    """
+    for s in sensors:
+        # `coverage` fills the domain with ones,
+        # adding twos will be visible in an image.
+        domain[y(s)][x(s)] = val
+    return domain
+
diff --git a/snp.py b/snp.py
new file mode 100644
index 0000000..4c060c0
--- /dev/null
+++ b/snp.py
@@ -0,0 +1,143 @@
+import numpy as np
+import matplotlib.pyplot as plt
+
+from sho import make, algo, iters, plot, num, bit, pb
+
+########################################################################
+# Interface
+########################################################################
+
+if __name__=="__main__":
+    import argparse
+
+    # Dimension of the search space.
+    d = 2
+
+    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")
+
+    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"]
+    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("-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")
+
+
+    the = can.parse_args()
+
+    # Minimum checks.
+    assert(0 < the.nb_sensors)
+    assert(0 < the.sensor_range <= 1)
+    assert(0 < the.domain_width)
+    assert(0 < the.iters)
+
+    # Do not forget the seed option,
+    # in case you would start "runs" in parallel.
+    np.random.seed(the.seed)
+
+    # Weird numpy way to ensure single line print of array.
+    np.set_printoptions(linewidth = np.inf)
+
+
+    # Common termination and checkpointing.
+    history = []
+    iters = make.iter(
+                iters.several,
+                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.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)
+                ]
+            )
+
+    # Erase the previous file.
+    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(
+                make.func(num.cover_sum,
+                    domain_width = the.domain_width,
+                    sensor_range = the.sensor_range * the.domain_width),
+                make.init(num.rand,
+                    dim = d * the.nb_sensors,
+                    scale = the.domain_width),
+                make.neig(num.neighb_square,
+                    scale = the.domain_width/10,
+                    domain_width = the.domain_width),
+                iters
+            )
+        sensors = num.to_sensors(sol)
+
+    elif the.solver == "bit_greedy":
+        val,sol = algo.greedy(
+                make.func(bit.cover_sum,
+                    domain_width = the.domain_width,
+                    sensor_range = the.sensor_range),
+                make.init(bit.rand,
+                    domain_width = the.domain_width,
+                    nb_sensors = the.nb_sensors),
+                make.neig(bit.neighb_square,
+                    scale = the.domain_width/10,
+                    domain_width = the.domain_width),
+                iters
+            )
+        sensors = bit.to_sensors(sol)
+
+    # Fancy output.
+    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)
+    ax2.imshow(domain)
+
+    plt.show()
diff --git a/snp_landscapes.py b/snp_landscapes.py
new file mode 100644
index 0000000..5339911
--- /dev/null
+++ b/snp_landscapes.py
@@ -0,0 +1,68 @@
+import numpy as np
+import matplotlib.pyplot as plt
+
+from sho import *
+
+def yonly_cover_sum(sol, domain_width, sensor_range, fixed_x = (10,30)):
+    """Compute the coverage quality of the given vector."""
+    domain = np.zeros((domain_width,domain_width))
+    sensors = [ (fixed_x[0], sol[0]), (fixed_x[1], sol[1]) ]
+    return np.sum(pb.coverage(domain, sensors, sensor_range))
+
+
+if __name__ == "__main__":
+
+    d = 2
+    w = 40
+    n = 2
+    r = 0.3 * w
+
+    # Common termination and checkpointing.
+    history = []
+    iters = make.iter(
+                iters.several,
+                agains = [
+                    make.iter(iters.max,
+                        nb_it = 100),
+                    make.iter(iters.log,
+                        fmt="\r{it} {val}"),
+                    make.iter(iters.history,
+                        history = history)
+                ]
+            )
+
+    x0,x1 = 0.25*w, 0.75*w
+
+    val,sol = algo.greedy(
+            make.func(yonly_cover_sum,
+                domain_width = w,
+                sensor_range = r,
+                fixed_x = (x0,x1) ),
+            make.init(num.rand,
+                dim = 2, # Two sensors moving along y axis.
+                scale = w),
+            make.neig(num.neighb_square,
+                scale = 0.1 * w),
+            iters
+        )
+    sensors = [ (int(x0), int(round(sol[0]))), (int(x1), int(round(sol[1]))) ]
+
+    print("\n{} : {}".format(val,sensors))
+
+    shape=(w,w)
+    fig = plt.figure()
+    ax1 = fig.add_subplot(121, projection='3d')
+    ax2 = fig.add_subplot(122)
+
+    f = make.func(yonly_cover_sum,
+                    domain_width = w,
+                    sensor_range = r)
+    plot.surface(ax1, shape, f)
+    plot.path(ax1, shape, history)
+
+    domain = np.zeros(shape)
+    domain = pb.coverage(domain, sensors, r)
+    domain = plot.highlight_sensors(domain, sensors)
+    ax2.imshow(domain)
+
+    plt.show()