refactor as a package
This commit is contained in:
parent
dcf9b798dc
commit
650d93585b
12 changed files with 477 additions and 520 deletions
41
sho/__init__.py
Normal file
41
sho/__init__.py
Normal file
|
|
@ -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 )
|
||||||
|
|
||||||
40
sho/algo.py
Normal file
40
sho/algo.py
Normal file
|
|
@ -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.
|
||||||
|
|
||||||
|
|
||||||
70
sho/api.py
70
sho/api.py
|
|
@ -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()
|
|
||||||
|
|
@ -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()
|
|
||||||
61
sho/bit.py
Normal file
61
sho/bit.py
Normal file
|
|
@ -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
|
||||||
|
|
||||||
40
sho/iters.py
Normal file
40
sho/iters.py
Normal file
|
|
@ -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
|
||||||
|
|
||||||
35
sho/make.py
Normal file
35
sho/make.py
Normal file
|
|
@ -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
|
||||||
|
|
||||||
48
sho/num.py
Normal file
48
sho/num.py
Normal file
|
|
@ -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
|
||||||
|
|
||||||
27
sho/pb.py
Normal file
27
sho/pb.py
Normal file
|
|
@ -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
|
||||||
|
|
||||||
|
|
||||||
75
sho/plot.py
75
sho/plot.py
|
|
@ -1,17 +1,31 @@
|
||||||
import numpy as np
|
import numpy as np
|
||||||
from matplotlib import cm
|
from matplotlib import cm
|
||||||
|
import matplotlib.pyplot as plt
|
||||||
import itertools
|
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):
|
def surface(ax, shape, f):
|
||||||
Z = np.zeros( shape )
|
Z = np.zeros( shape )
|
||||||
for y in range(shape[0]):
|
for y in range(shape[0]):
|
||||||
for x in range(shape[1]):
|
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)
|
X = np.arange(0,shape[0],1)
|
||||||
Y = np.arange(0,shape[1],1)
|
Y = np.arange(0,shape[1],1)
|
||||||
X,Y = np.meshgrid(X,Y)
|
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 path(ax, shape, history):
|
||||||
def pairwise(iterable):
|
def pairwise(iterable):
|
||||||
|
|
@ -21,11 +35,11 @@ def path(ax, shape, history):
|
||||||
|
|
||||||
k=0
|
k=0
|
||||||
for i,j in pairwise(range(len(history)-1)):
|
for i,j in pairwise(range(len(history)-1)):
|
||||||
xi = history[i][1][0]*shape[0]
|
xi = history[i][1][0]
|
||||||
yi = history[i][1][1]*shape[1]
|
yi = history[i][1][1]
|
||||||
zi = history[i][0]
|
zi = history[i][0]
|
||||||
xj = history[j][1][0]*shape[0]
|
xj = history[j][1][0]
|
||||||
yj = history[j][1][1]*shape[1]
|
yj = history[j][1][1]
|
||||||
zj = history[j][0]
|
zj = history[j][0]
|
||||||
x = [xi, xj]
|
x = [xi, xj]
|
||||||
y = [yi, yj]
|
y = [yi, yj]
|
||||||
|
|
@ -33,3 +47,52 @@ def path(ax, shape, history):
|
||||||
ax.plot(x,y,z, color=cm.RdYlBu(k))
|
ax.plot(x,y,z, color=cm.RdYlBu(k))
|
||||||
k+=1
|
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()
|
||||||
|
|
|
||||||
374
sho/snp.py
374
sho/snp.py
|
|
@ -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()
|
|
||||||
|
|
||||||
116
snp.py
Normal file
116
snp.py
Normal file
|
|
@ -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()
|
||||||
|
|
||||||
Loading…
Add table
Add a link
Reference in a new issue