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..3ff1faf --- /dev/null +++ b/sho/algo.py @@ -0,0 +1,40 @@ + +######################################################################## +# Algorithms +######################################################################## + +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, val, 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) + if val > best_val: + best_val = val + best_sol = sol + i += 1 + return best_val, best_sol + +# TODO add a simulated annealing solver. +# TODO add a population-based stochastic heuristic template. + + diff --git a/sho/api.py b/sho/api.py deleted file mode 100644 index aa43245..0000000 --- a/sho/api.py +++ /dev/null @@ -1,70 +0,0 @@ -import numpy as np - -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 f - -def square(sol,scale=1): - """Gnerate a random vector close at thegiven scale to the given sol.""" - return sol + np.random.random(len(sol))*scale - -def greedy(objective_function, dimension, iterations, target=1e-3, neighborhood=square, scale=1/100, history=None): - """Search the given objective_function of the given dimension, - during the given number of iterations, generating solution - with the given neighborhood. - Returns the best value of the function and the best solution.""" - best_sol = np.random.random(dimension) - best_val = objective_function(best_sol) - for i in range(iterations): - sol = neighborhood(best_sol,scale) - val = objective_function(sol) - if val < best_val: - best_val = val - best_sol = sol - if history is not None: - history.append((val,sol)) - if val < target: # Assume the optimum is zero - break - return best_val, best_sol - - -if __name__=="__main__": - import matplotlib.pyplot as plt - from mpl_toolkits.mplot3d import Axes3D - import argparse - import plot - - parser = argparse.ArgumentParser() - - parser.add_argument("-d", "--dim", metavar="NB", default=2, type=int, - help="Number of dimensions") - - functions = {"sphere":sphere} - parser.add_argument("-f", "--func", metavar="NAME", choices=functions, default="sphere", - help="Objective function") - - parser.add_argument("-i", "--iter", metavar="NB", default=1000, type=int, - help="Maximum number of iterations") - - parser.add_argument("-t", "--target", metavar="VAL", default=1e-3, type=float, - help="Function value target delta") - - parser.add_argument("-s", "--seed", metavar="VAL", default=0, type=int, - help="Random pseudo-generator seed (0 for epoch)") - - asked = parser.parse_args() - - np.random.seed(asked.seed) - - history = [] - val,sol = greedy(functions[asked.func], asked.dim, asked.iter, asked.target, square, 0.03, history) - - fig = plt.figure() - ax = fig.gca(projection='3d') - shape = (20,20) - plot.surface(ax, shape, sphere) - plot.path(ax, shape, history) - plt.show() diff --git a/sho/basics.py b/sho/basics.py deleted file mode 100644 index f60a762..0000000 --- a/sho/basics.py +++ /dev/null @@ -1,70 +0,0 @@ -import numpy as np - -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 f - -def onemax(x): - """Sum the given bitstring.""" - s = 0 - for i in x: - s += i - return s - -def numerical_random(d): - """Draw a random multi-dimensional vector in [0,1]**d""" - return np.random.random(d) - -def bitstring_random(d): - """Draw a random bistring of size d, with P(1)=0.5.""" - return [int(round(i)) for i in np.random.random(d)] - -def search(objective_function, dimension, iterations, generator, history=None): - """Search the given objective_function of the given dimension, - during the given number of iterations, generating random solution - with the given generator. - Returns the best value of the function and the best solution.""" - best_val = float("inf") - best_sol = None - for i in range(iterations): - sol = generator(dimension) - val = objective_function(sol) - if val < best_val: - best_val = val - best_sol = sol - if history is not None: - history.append((val,sol)) - return best_val, best_sol - - -if __name__=="__main__": - import matplotlib.pyplot as plt - from mpl_toolkits.mplot3d import Axes3D - import plot - - print("Random search over 10-OneMax") - print("After 10 iterations:") - val,sol = search(onemax, 10, 10, bitstring_random) - print("\t",val,sol) - print("After 1000 iterations:") - val,sol = search(onemax, 10, 1000, bitstring_random) - print("\t",val,sol) - - print("Random search over 2-Sphere") - print("After 10 iterations:") - val,sol = search(sphere, 2, 10, numerical_random) - print("\t",val,sol) - print("After 50 iterations:") - history = [] - val,sol = search(sphere, 2, 50, numerical_random, history) - print("\t",val,sol) - - fig = plt.figure() - ax = fig.gca(projection='3d') - shape = (20,20) - plot.surface(ax, shape, sphere) - plot.path(ax, shape, history) - plt.show() diff --git a/sho/bit.py b/sho/bit.py new file mode 100644 index 0000000..866f54e --- /dev/null +++ b/sho/bit.py @@ -0,0 +1,61 @@ +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): + """Draw a random array by moving ones to adjacent cells.""" + # Copy, because Python pass by reference + # and we may not the 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. + # TODO handle constraints + d = np.random.randint(-scale//2,scale//2,2) + 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..7744c0d --- /dev/null +++ b/sho/iters.py @@ -0,0 +1,40 @@ +import sys + +######################################################################## +# Stopping criterions +######################################################################## + +def max(i, val, sol, nb_it): + if i < nb_it: + return True + else: + return False + +# Stopping criterions that are actually just checkpoints. + +def several(i, val, sol, agains): + """several 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..1e7bbe6 --- /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): + """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..89c5c45 --- /dev/null +++ b/sho/pb.py @@ -0,0 +1,27 @@ +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 + + diff --git a/sho/plot.py b/sho/plot.py index dd5a7ae..251dca4 100644 --- a/sho/plot.py +++ b/sho/plot.py @@ -1,17 +1,31 @@ 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/shape[0],y/shape[1]), 0.5 ) + Z[y][x] = f( (x,y), shape[0]/2 ) 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, cmap=cm.viridis) + ax.plot_surface(X, Y, Z) + def path(ax, shape, history): def pairwise(iterable): @@ -21,11 +35,11 @@ def path(ax, shape, history): k=0 for i,j in pairwise(range(len(history)-1)): - xi = history[i][1][0]*shape[0] - yi = history[i][1][1]*shape[1] + xi = history[i][1][0] + yi = history[i][1][1] zi = history[i][0] - xj = history[j][1][0]*shape[0] - yj = history[j][1][1]*shape[1] + xj = history[j][1][0] + yj = history[j][1][1] zj = history[j][0] x = [xi, xj] y = [yi, yj] @@ -33,3 +47,52 @@ def path(ax, shape, history): 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 + + +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/sho/snp.py b/sho/snp.py deleted file mode 100644 index 0914842..0000000 --- a/sho/snp.py +++ /dev/null @@ -1,374 +0,0 @@ -import sys -import numpy as np -import matplotlib.pyplot as plt -import copy - - -######################################################################## -# 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 ) - - -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 - - -######################################################################## -# 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. - - >>> snp.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 - - -# Decoupled from objective functions, so as to be used in display. -def num_to_sensors(sol): - """Convert a vector of n*2 dimension to an array of n 2-tuples. - - >>> num_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 bit_to_sensors(sol): - """Convert an square array of d lines/columns containing n ones - to an array of n 2-tuples with related coordinates. - - >>> bit_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 - - -def bit_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 = bit_to_sensors(sol) - return np.sum(coverage(domain, sensors, sensor_range)) - - -def num_cover_sum(sol, domain_width, sensor_range): - """Compute the coverage quality of the given vector.""" - domain = np.zeros((domain_width,domain_width)) - sensors = num_to_sensors(sol) - return np.sum(coverage(domain, sensors, sensor_range)) - - -def make_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 - - -######################################################################## -# Initialization -######################################################################## - -def num_rand(dim, scale): - """Draw a random vector in [0,scale]**dim.""" - return np.random.random(dim) * scale - - -def bit_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 - - -def make_init(init, **kwargs): - """Make an initialization operator from the given function. - An init. op. returns a solution.""" - def f(): - return init(**kwargs) - return f - - -######################################################################## -# Neighborhood -######################################################################## - -def num_neighb_square(sol, scale): - """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 - - -def bit_neighb_square(sol, scale): - """Draw a random array by moving ones to adjacent cells.""" - # Copy, because Python pass by reference - # and we may not the 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. - # TODO handle constraints - d = np.random.randint(-scale//2,scale//2,2) - new[py+y(d)][px+x(d)] = 1 - return new - - -def make_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 - - -######################################################################## -# Stopping criterions -######################################################################## - -def iter_max(i, val, sol, nb_it): - if i < nb_it: - return True - else: - return False - -def make_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 - - -# Stopping criterions that are actually just checkpoints. - -def combine(i, val, sol, agains): - """Combine several stopping criterions in one.""" - res = True - for again in agains: - res = res and again(i, val, sol) - return res - - -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 iter_log(i, val, sol, fmt="{it} {val}\n"): - """Print progress on stderr.""" - sys.stderr.write( fmt.format(it=i, val=val) ) - return True - - -######################################################################## -# Algorithms -######################################################################## - -def random(func, init, again): - """Iterative random search template.""" - best_sol = None - best_val = - np.inf - val,sol = best_val,best_sol - i = 0 - while again(i, val, 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, val, sol): - sol = neighb(best_sol) - val = func(sol) - if val > best_val: - best_val = val - best_sol = sol - i += 1 - return best_val, best_sol - -# TODO add a simulated annealing solver. -# TODO add a population-based stochastic heuristic template. - - -######################################################################## -# 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=100, 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)) - - # TODO add the corresponding stopping criterion. - can.add_argument("-t", "--target", metavar="VAL", default=1e-3, type=float, - help="Function value target delta") - - 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) - - domain = np.zeros((the.domain_width, the.domain_width)) - - # Common termination and checkpointing. - iters = make_iter( - combine, - agains = [ - make_iter(iter_max, - nb_it = the.iters), - make_iter(save, - filename = the.solver+".csv", - fmt = "{it} ; {val} ; {sol}\n"), - make_iter(iter_log, - fmt="\r{it} {val}") - ] - ) - - # Erase the previous file. - with open(the.solver+".csv", 'w') as fd: - fd.write("# {} {}\n".format(the.solver,the.domain_width)) - - if the.solver == "num_greedy": - val,sol = 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), # TODO think of an alternative. - iters - ) - sensors = num_to_sensors(sol) - - elif the.solver == "bit_greedy": - val,sol = 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), - iters - ) - sensors = bit_to_sensors(sol) - - - # Fancy output. - print("\n",val,":",sensors) - - domain = coverage(domain, sensors, - the.sensor_range * the.domain_width) - domain = highlight_sensors(domain, sensors) - plt.imshow(domain) - plt.show() - diff --git a/snp.py b/snp.py new file mode 100644 index 0000000..2732543 --- /dev/null +++ b/snp.py @@ -0,0 +1,116 @@ +import sys +import numpy as np +import matplotlib.pyplot as plt +import copy + +from sho import * + +######################################################################## +# 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=100, 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)) + + # TODO add the corresponding stopping criterion. + can.add_argument("-t", "--target", metavar="VAL", default=1e-3, type=float, + help="Function value target delta") + + 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) + + domain = np.zeros((the.domain_width, the.domain_width)) + + # Common termination and checkpointing. + 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}") + ] + ) + + # 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), # TODO think of an alternative. + 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), + iters + ) + sensors = bit.to_sensors(sol) + + + # Fancy output. + print("\n",val,":",sensors) + + domain = pb.coverage(domain, sensors, + the.sensor_range * the.domain_width) + domain = plot.highlight_sensors(domain, sensors) + plt.imshow(domain) + plt.show() +