sobol sequence solver

bin/examples/python/parabola.py

@@ -11,7 +11,7 @@
 
 # we run experiments a number of times to estimate each solver's variance
 particle_details = None
-for i in range(200):
+for i in range(100):
     xoff = random.random()
     yoff = random.random()
     def f(x, y):
@@ -27,8 +27,8 @@ def f(x, y):
 
 # plot results
 print('plotting results')
-colors = ['r', 'g', 'b', 'y', 'k', 'y']
-markers = ['x', '+', 'o', 's', 'p']
+colors =  ['r', 'g', 'b', 'y', 'k', 'y', 'r', 'g']
+markers = ['x', '+', 'o', 's', 'p', 'x', '+', 'o']
 
 delta = 0.025
 x = np.arange(-5.0, 5.0, delta)
@@ -42,6 +42,7 @@ def f(x, y):
 CS = plt.contour(X, Y, Z)
 plt.clabel(CS, inline=1, fontsize=10, alpha=0.5)
 for i, solver in enumerate(solvers):
+    print(solver)
     plt.scatter(logs[solver][0,:], logs[solver][1,:], c=colors[i], marker=markers[i], alpha=0.80)
 
 plt.xlim([-5, 5])
@@ -67,16 +68,4 @@ def f(x, y):
 plt.yticks(np.arange(len(means))+0.4, list(means.keys()))
 plt.tight_layout()
 plt.show()
-#plt.savefig('parabola_solver_precision.png', transparant=True)
-
-
-
-
-num_evals = len(particle_details.call_log['values'])
-best_fitness = [min(particle_details.call_log['values'][:x]) for x in range(1, 1 + num_evals)]
-avg_fitness = [sum(particle_details.call_log['values'][:x]) / x for x in range(1, 1 + num_evals)]
-
-plt.figure(3)
-plt.plot(range(num_evals), best_fitness, 'b')
-plt.plot(range(num_evals), avg_fitness, 'r')
-plt.show()
+plt.savefig('parabola_solver_precision.png', transparant=True)

bin/examples/python/sobol_vs_uniform.py

@@ -0,0 +1,32 @@
+# In this example we will show the difference between a 2-d Sobol sequence
+# and sampling uniformly at random in 2 dimensions.
+# The Sobol sequence has far lower discrepancy, i.e., the generated samples
+# are spread out better in the sampling space.
+#
+# This example requires matplotlib to generate figures.
+
+import matplotlib.pyplot as plt
+import optunity
+import random
+
+num_pts = 200   # the number of points to generate
+skip = 5000     # the number of initial points of the Sobol sequence to skip
+
+# generate Sobol sequence
+res = optunity.solvers.Sobol.i4_sobol_generate(2, num_pts, skip)
+x1_sobol, x2_sobol = zip(*res)
+
+# generate uniform points
+x1_random = [random.random() for _ in range(num_pts)]
+x2_random = [random.random() for _ in range(num_pts)]
+
+# plot results
+plt.figure(1)
+plt.plot(x1_sobol, x2_sobol, 'o')
+plt.title('Sobol sequence')
+plt.draw()
+
+plt.figure(2)
+plt.plot(x1_random, x2_random, 'ro')
+plt.title('Uniform random samples')
+plt.show()

docs/user/solvers/sobol.rst

@@ -8,6 +8,7 @@ hypercube with lower discrepancy than completely random sampling (e.g. |randomse
 sequences in up to 40 dimensions (e.g. 40 hyperparameters). 
 
 The figures below show the differences between a Sobol sequence and sampling uniformly at random.
+These figures can be recreated using the code in `bin/examples/python/sobol_vs_random.py`.
 
 .. figure:: sobol.png
     :alt: 200 points sampled in 2D with a Sobol sequence.

optunity/solvers/ParticleSwarm.py

@@ -39,6 +39,7 @@
 from .solver_registry import register_solver
 from .util import Solver, _copydoc, uniform_in_bounds
 from . import util
+from .Sobol import Sobol
 
 @register_solver('particle swarm',
                  'particle swarm optimization',
@@ -128,6 +129,8 @@ def __init__(self, num_particles, num_generations, max_speed=None, phi1=1.5, phi
         self._num_particles = num_particles
         self._num_generations = num_generations
 
+        self._sobolseed = random.randint(100,2000)
+
         if max_speed is None:
             max_speed = 0.7 / num_generations
 #            max_speed = 0.2 / math.sqrt(num_generations)
@@ -147,6 +150,12 @@ def phi1(self):
     def phi2(self):
         return self._phi2
 
+    @property
+    def sobolseed(self): return self._sobolseed
+
+    @sobolseed.setter
+    def sobolseed(self, value): self._sobolseed = value
+
     @staticmethod
     def suggest_from_box(num_evals, **kwargs):
         """Create a configuration for a ParticleSwarm solver.
@@ -210,7 +219,12 @@ def bounds(self):
 
     def generate(self):
         """Generate a new Particle."""
-        part = ParticleSwarm.Particle(position=array.array('d', uniform_in_bounds(self.bounds)),
+        if len(self.bounds) < Sobol.maxdim():
+            sobol_vector, self.sobolseed = Sobol.i4_sobol(len(self.bounds), self.sobolseed)
+            vector = util.scale_unit_to_bounds(sobol_vector, self.bounds.values())
+        else: vector = uniform_in_bounds(self.bounds)
+
+        part = ParticleSwarm.Particle(position=array.array('d', vector),
                                       speed=array.array('d', map(random.uniform,
                                                                  self.smin, self.smax)),
                                       best=None, fitness=None, best_fitness=None)