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

snp.py

@@ -48,6 +48,9 @@ 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")
 
@@ -165,26 +168,28 @@ if __name__=="__main__":
     # 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()