draw_ert.py

+#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')

sho/algo.py

@@ -35,7 +35,25 @@ def greedy(func, init, neighb, again):
         i += 1
     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.
 

sho/temp.py

+
+########################################################################
+# 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

snp.py

@@ -3,7 +3,7 @@ import math
 import numpy as np
 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
@@ -32,7 +32,7 @@ if __name__=="__main__":
     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"]
+    solvers = ["num_greedy","bit_greedy","num_simulated_annealing","bit_simulated_annealing"]
     can.add_argument("-m", "--solver", metavar="NAME", choices=solvers, default="num_greedy",
             help="Solver to use, among: "+", ".join(solvers))
 
@@ -48,6 +48,12 @@ if __name__=="__main__":
     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)")
 
+    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()
 
     # Minimum checks.
@@ -63,6 +69,10 @@ if __name__=="__main__":
     # Weird numpy way to ensure single line print of array.
     np.set_printoptions(linewidth = np.inf)
 
+    # Default filename
+    if the.output_file == None:
+        the.output_file = the.solver
+
 
     # Common termination and checkpointing.
     history = []
@@ -72,7 +82,7 @@ if __name__=="__main__":
                     make.iter(iters.max,
                         nb_it = the.iters),
                     make.iter(iters.save,
-                        filename = the.solver+".csv",
+                        filename = the.output_file+".csv",
                         fmt = "{it} ; {val} ; {sol}\n"),
                     make.iter(iters.log,
                         fmt="\r{it} {val}"),
@@ -121,29 +131,65 @@ if __name__=="__main__":
             )
         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.
     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:
-        ax1 = fig.add_subplot(121, projection='3d')
-        ax2 = fig.add_subplot(122)
+        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)
+            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)
+        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()
+        plt.show()