Johannes Buchner

If you're new here, you may want to subscribe to my RSS feed. Thanks for visiting!

In many cases, the parameter space (e.g. simulation parameters) is a n-dimensional box/cube, scalable to [0:1]^n, where n donates the number of dimensions = number of parameters.

Since simulations for one parameter set take time, and one wants to exhaust the parameter space efficiently (visit every subspace uniformly), one has to come up with some method of generating or suggesting points to visit.

Often, a Monte-Carlo approach is used by just using random points (uniformly, to be exact).

Here I want to present a simple enumeration scheme that eventually visits all points, by ordering the  rational numbers in an suitable way.

First, lets look at one dimension. The strategy is to visit: 0, 1, 1/2, 1/4, 3/4. Then continue even deeper with 1/8, 3/8, 5/8, 7/8 and so on. Using fractions, a pattern is clearly visible. Here is a visualization:

The y-axis shows the one dimensional parameter space, the x-axis just helps you to see the ordering of the enumeration. Each animation step goes one step deeper.

The generating scheme up to a certain deepness is:

{ (j*2+1) / 2^deepness | j in [0..(2^deepness) / 2] }

Higher dimensions

So far, so easy. Lets generalize this to higher dimensions. We want to keep the property that the enumeration sort of “goes deeper”, becomes more accurate by filling the voids between the previously visited points.

Here, a two-dimensional parameter space is used. First, the edges (0,0),(1,0),(0,1),(1,1) are visited. Then the space in between and so on and so on. Each animation step goes one step deeper.

The easiest method for generating a enumeration scheme for arbitrary dimensions is

0. no visited points
1. Q1 <- generate enumerations of one dimension, up to deepness o
2. build all permutations [a_1,a_2,...,a_3] where a_i in Q1
3. remove already visited points
4. increase deepness, go to 1

Use: The more points you generate, the more you can exhaust the parameter space and the better the results get. At some point you have to stop your calculations, and this enumeration makes sure you have visited all subspaces homogeneously.

Appendix

Appended is a python script (numbering.py) I used to generate points in this enumerated way. the parameter space can easily be scaled by piping into awk (e.g. | awk '{print $1*3,$2+1}').

If you don’t want to or can’t run the script, the output for 1, 2 and 3 dimensional parameter spaces can be found in http://johannes.jakeapp.com/files/enumerate_rational/. These contains the first million points (deepness 22, 11 and 8).

import math
import sys 
 
def getSegment(i):
	if i == 0:
		return 0.
	#if i == 1:
	#	return 1
	o = int(math.ceil(math.log(i, 2)))
	j = i - (2**o / 2) - 1
	#print "i=", i, "o=", o, "j=", j
	return float(j*2+1) / 2**o
 
def getSegmentsInDeepnessRange(a, b):
	if a == -1:
		first = 0
	else:
		first = 2**a + 1
	last  = 2**b + 1
	# print "o=", o, "first=", first, "last=", last
	return [getSegment(j) for j in range(first, last)]
 
#crossproduct = lambda ss,row=[],level=0: len(ss)>1 \
#	and reduce(lambda x,y:x+y,[crossproduct(ss[1:],row+[i],level+1) for i in ss[0]]) \
#	or [row+[i] for i in ss[0]]
def crossproduct(ss, row=[], level=0):
	if len(ss)>1:
		subcrossproduct = [crossproduct(ss[1:],row+[i],level+1) for i in ss[0]]
		return reduce(lambda x,y:x+y,subcrossproduct)
	else:
		return [row+[i] for i in ss[0]]
 
if len(sys.argv) != 4:
	print "SYNOPSIS: dim min_deepness max_deepness"
	print
	print "\tdim\tdimensions of the parameter space"
	print "\tdeepness\thow many stages to go deep (set min=0)"
	print
	print "Q: how many numbers will be produced?"
	print "A: when min=0: (2^(max-1)+1)^dim; as a table:"
	print "  max   |    dim=1           dim=2     ->"
	for i in range(20):
		print "%6d" % i,
		for dim in range(1,6):
			if ((2**(i-1)+1)**dim) > 10000000:
				print "\t    ...    ",
			else:
				print "\t%10d" % (2**(i-1)+1)**dim,
		print
	exit()
 
dim = int(sys.argv[1])
 
i = 0
 
mindeepness = int(sys.argv[2])
maxdeepness = int(sys.argv[3])
 
if dim > 1:
	# o = deepness
	for o in range(mindeepness, maxdeepness):
		p    = getSegmentsInDeepnessRange(-1, o)
		pnew = getSegmentsInDeepnessRange(o-1,  o)
		if dim == 1:
			newvals = [[nv] for nv in p]
		else:
			newvals = crossproduct([p]*dim)
		# print only the new ones:
		for nv in newvals:
			for pn in pnew:
				if pn in nv:
					for v in nv:
						print v,
					print
					i = i + 1
					break
else: # dim = 1:
#for i in range(20):
#	print getSegment(i)
 
# one dimension 
 
	i = 0
	print i
	i = 1
	print i
	for o in range(mindeepness, maxdeepness):
		for j in range((2**o) / 2):
			print ("%.30f" % ((j*2+1) * pow(2,-o))).rstrip('0')