Commits (11)
#encoding: utf-8
from os import system, remove, path
import matplotlib.pyplot as plt
import numpy as np
########################################################################
# Interface
########################################################################
if __name__=="__main__":
import argparse
can = argparse.ArgumentParser()
can.add_argument("-N", "--nb-run", metavar="NB", default=100, type=int,
help="Number of runs")
can.add_argument("-v", "--for-value", metavar="NB", type=int,
help="Print ERT with probabilities to obtain this value")
arg = can.parse_known_args()
the = arg[0]
print(the)
# Minimum checks.
assert(0 < the.nb_run)
if path.isfile('run.csv'):
remove('run.csv')
for i in range(the.nb_run):
print("\rRun {}".format(i), end = '')
system("python3 snp.py -o run -no "+" ".join(arg[1]))
print("\r")
# Data extraction
data = []
attempt = {}
current_best = 0
temps_max = 0
value_max = 0
with open('run.csv', 'r') as fd:
for line in fd:
line = line.split(' ; ')
line[0] = int(line[0])
line[1] = float(line[1])
if line[0] > temps_max:
temps_max = line[0]
if line[1] > value_max:
value_max = line[1]
if line[0] == 1:
if len(attempt) > 0:
data.append(attempt)
current_best = float("inf")
attempt = {}
if line[1] < current_best:
attempt[line[0]] = line[1]
current_best = line[1]
data.append(attempt)
def probability(data,time,value):
num = 0
for attempt in data:
good = False
for t,v in attempt.items():
if t < time and v < value:
good = True
break
if good:
num += 1
return num/len(data)
probability = np.vectorize(probability, excluded=[0])
# Plot
if the.for_value is not None:
temps = np.linspace(0,temps_max,100)
P = probability(data,temps,the.for_value)
plt.plot(temps,P)
plt.show()
else:
temps = np.linspace(0,temps_max,100)
values = np.linspace(0,value_max,100)
X, Y = np.meshgrid(temps, values)
Z = probability(data,X,Y)
fig, ax = plt.subplots()
cs = ax.contourf(X, Y, Z, 20)
cbar = fig.colorbar(cs)
plt.show()
remove('run.csv')
...@@ -35,7 +35,25 @@ def greedy(func, init, neighb, again): ...@@ -35,7 +35,25 @@ def greedy(func, init, neighb, again):
i += 1 i += 1
return best_val, best_sol return best_val, best_sol
# TODO add a simulated-annealing template.
def simulated_annealing(func, init, neighb, temp, again):
"""Iterative randomized simulated-annealing 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
elif np.random.random() < np.exp(-(best_val-val)/temp(i, best_val, best_sol)):
best_val = val
best_sol = sol
i += 1
return best_val, best_sol
# TODO add a population-based stochastic heuristic template. # TODO add a population-based stochastic heuristic template.
########################################################################
# Temperature programme
########################################################################
class continus:
"""Continus decrease temperature from T_0 and with T_{i+1} = \alpha x T_i"""
def __init__(self, T_0, alpha):
self.T = T_0
self.alpha = alpha
def __call__(self, i, val, sol):
T = self.T
self.T *= self.alpha
return T
class level:
"""Level decrease temperature from T_0 and with T_{i+1} = T_i - \alpha after \beta changes"""
def __init__(self, T_0, alpha, beta):
self.T = T_0
self.alpha = alpha
self.beta = beta
self.val = None
self.count = 0
def __call__(self, i, val, sol):
T = self.T
if val != self.val:
self.count += 1
self.val = val
if self.count > self.beta:
self.count = 0
self.T -= self.alpha
print()
print("T = "+str(T))
return T
\ No newline at end of file
...@@ -3,7 +3,7 @@ import math ...@@ -3,7 +3,7 @@ import math
import numpy as np import numpy as np
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
from sho import make, algo, iters, plot, num, bit, pb from sho import make, algo, iters, plot, num, bit, pb, temp
######################################################################## ########################################################################
# Interface # Interface
...@@ -32,7 +32,7 @@ if __name__=="__main__": ...@@ -32,7 +32,7 @@ if __name__=="__main__":
can.add_argument("-s", "--seed", metavar="VAL", default=None, type=int, can.add_argument("-s", "--seed", metavar="VAL", default=None, type=int,
help="Random pseudo-generator seed (none for current epoch)") help="Random pseudo-generator seed (none for current epoch)")
solvers = ["num_greedy","bit_greedy"] solvers = ["num_greedy","bit_greedy","num_simulated_annealing","bit_simulated_annealing"]
can.add_argument("-m", "--solver", metavar="NAME", choices=solvers, default="num_greedy", can.add_argument("-m", "--solver", metavar="NAME", choices=solvers, default="num_greedy",
help="Solver to use, among: "+", ".join(solvers)) help="Solver to use, among: "+", ".join(solvers))
...@@ -48,6 +48,12 @@ if __name__=="__main__": ...@@ -48,6 +48,12 @@ if __name__=="__main__":
can.add_argument("-a", "--variation-scale", metavar="RATIO", default=0.3, type=float, can.add_argument("-a", "--variation-scale", metavar="RATIO", default=0.3, type=float,
help="Scale of the variation operators (as a ration of the domain width)") help="Scale of the variation operators (as a ration of the domain width)")
can.add_argument("-no", "--no-graphical-output", action='store_true', default=False,
help="Disable graphical output")
can.add_argument("-o", "--output-file", metavar="NAME", default=None, type=str,
help="Output file where log are put")
the = can.parse_args() the = can.parse_args()
# Minimum checks. # Minimum checks.
...@@ -63,6 +69,10 @@ if __name__=="__main__": ...@@ -63,6 +69,10 @@ if __name__=="__main__":
# Weird numpy way to ensure single line print of array. # Weird numpy way to ensure single line print of array.
np.set_printoptions(linewidth = np.inf) np.set_printoptions(linewidth = np.inf)
# Default filename
if the.output_file == None:
the.output_file = the.solver
# Common termination and checkpointing. # Common termination and checkpointing.
history = [] history = []
...@@ -72,7 +82,7 @@ if __name__=="__main__": ...@@ -72,7 +82,7 @@ if __name__=="__main__":
make.iter(iters.max, make.iter(iters.max,
nb_it = the.iters), nb_it = the.iters),
make.iter(iters.save, make.iter(iters.save,
filename = the.solver+".csv", filename = the.output_file+".csv",
fmt = "{it} ; {val} ; {sol}\n"), fmt = "{it} ; {val} ; {sol}\n"),
make.iter(iters.log, make.iter(iters.log,
fmt="\r{it} {val}"), fmt="\r{it} {val}"),
...@@ -121,29 +131,65 @@ if __name__=="__main__": ...@@ -121,29 +131,65 @@ if __name__=="__main__":
) )
sensors = bit.to_sensors(sol) sensors = bit.to_sensors(sol)
elif the.solver == "num_simulated_annealing":
val,sol = algo.simulated_annealing(
make.func(num.cover_sum,
domain_width = the.domain_width,
sensor_range = the.sensor_range,
dim = d * the.nb_sensors),
make.init(num.rand,
dim = d * the.nb_sensors,
scale = the.domain_width),
make.neig(num.neighb_square,
scale = the.variation_scale,
domain_width = the.domain_width),
temp.continus(100,0.99),
iters
)
sensors = num.to_sensors(sol)
elif the.solver == "bit_simulated_annealing":
val,sol = algo.simulated_annealing(
make.func(bit.cover_sum,
domain_width = the.domain_width,
sensor_range = the.sensor_range,
dim = d * the.nb_sensors),
make.init(bit.rand,
domain_width = the.domain_width,
nb_sensors = the.nb_sensors),
make.neig(bit.neighb_square,
scale = the.variation_scale,
domain_width = the.domain_width),
temp.continus(100,0.99),
iters
)
sensors = bit.to_sensors(sol)
# Fancy output. # Fancy output.
print("\n{} : {}".format(val,sensors)) print("\n{} : {}".format(val,sensors))
shape=(the.domain_width, the.domain_width) if the.no_graphical_output is False:
shape=(the.domain_width, the.domain_width)
fig = plt.figure() fig = plt.figure()
if the.nb_sensors ==1 and the.domain_width <= 50: if the.nb_sensors ==1 and the.domain_width <= 50:
ax1 = fig.add_subplot(121, projection='3d') ax1 = fig.add_subplot(121, projection='3d')
ax2 = fig.add_subplot(122) ax2 = fig.add_subplot(122)
f = make.func(num.cover_sum, f = make.func(num.cover_sum,
domain_width = the.domain_width, domain_width = the.domain_width,
sensor_range = the.sensor_range * the.domain_width) sensor_range = the.sensor_range * the.domain_width)
plot.surface(ax1, shape, f) plot.surface(ax1, shape, f)
plot.path(ax1, shape, history) plot.path(ax1, shape, history)
else: else:
ax2=fig.add_subplot(111) ax2=fig.add_subplot(111)
domain = np.zeros(shape) domain = np.zeros(shape)
domain = pb.coverage(domain, sensors, domain = pb.coverage(domain, sensors,
the.sensor_range * the.domain_width) the.sensor_range * the.domain_width)
domain = plot.highlight_sensors(domain, sensors) domain = plot.highlight_sensors(domain, sensors)
ax2.imshow(domain) ax2.imshow(domain)
plt.show() plt.show()