Commit a90792a3 authored by Johann Dreo's avatar Johann Dreo
Browse files

first import

parents
numpy
matplotlib
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 )
########################################################################
# 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.
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
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
"""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
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
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
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
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()
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}"),