use new methods in test suite; added a few new test problems

Author: JohannesBuchner

.gitignore

@@ -3,3 +3,4 @@
 *.pdf
 [_]build/
 *.egg-info/
+*.so

setup.py

@@ -7,7 +7,7 @@
 	long_description = f.read()
 
 setup(name='nested_sampling',
-	version='0.1',
+	version='0.2',
 	author='Johannes Buchner',
 	url='https://bitbucket.org/JohannesBuchner/nested-sampling',
 	author_email='[email protected]',

testsuite/algorithms/__init__.py

@@ -4,13 +4,12 @@
 import multinest
 import cuba
 import nest
-#import runnestle
+import runnestle
 
 algorithms = \
 	multinest.configs + \
+	runnestle.configs + \
 	cuba.configs + \
 	nest.configs
-	# + \
-	#runnestle.configs
 
 

testsuite/algorithms/multinest.py

@@ -25,6 +25,10 @@ def myloglike(cube, ndim, nparams):
 			sys.exit(-127)
 		return l
 	nlive_points = config['nlive_points']
+	if config.get('unlimited_sampling', False):
+		max_samples = 0
+	else:
+		max_samples = 2000000
 	mn_args = dict(
 		importance_nested_sampling = config['importance_nested_sampling'],
 		outputfiles_basename = output_basename + 'out_',
@@ -36,7 +40,7 @@ def myloglike(cube, ndim, nparams):
 		const_efficiency_mode = False,
 		evidence_tolerance = 0.5,
 		seed = config['seed'],
-		max_iter = 2000000,
+		max_iter = max_samples,
 	)
 	starttime = time.time()
 	pymultinest.run(myloglike, myprior, mn_args['n_params'], **mn_args)

testsuite/algorithms/nest.py

@@ -5,7 +5,7 @@
 import nestle
 import numpy
 from numpy import log, exp
-
+import sys
 import time
 
 def run_nestle(**config):
@@ -14,21 +14,38 @@ def run_nestle(**config):
 	def priortransform(u):
 		assert len(u) == ndim, u
 		return u
+	def dump_callback(info):
+		sys.stderr.write("\r%d|%d|logz=%.4f|eff=%f%%    " % (info['it'], info['ncall'], info['logz'], info['it']*100./info['ncall']))
+		
+		#if info['it'] % 50 != 0: return
+		#print "Replacements: %d" % (info['ncall'])
+		#print "Samples: %d" % (info['it'])
+		#print "Efficiency: %f" % (info['ncall']/info['it'])
+		#print "Nested Sampling ln(Z): %f" % (info['logz'])
 	if 'seed' in config:
 		numpy.random.seed(config['seed'])
 
 	# can use directly
 	loglikelihood = config['loglikelihood']
 	nlive_points = config['nlive_points']
 	method = config['method']
-	
+	if config.get('unlimited_sampling', False):
+		max_samples = None
+	else:
+		max_samples = 2000000
+	print
 	print 'running nestle ...'
+	options = dict()
+	#if 'enlarge' in config:
+	#	options['enlarge'] = config['enlarge']
 	starttime = time.time()
 	result = nestle.sample(loglikelihood=loglikelihood, prior_transform=priortransform, ndim=ndim, npoints=nlive_points,
-		method=method, update_interval=None, maxcall=2000000, dlogz=0.5, rstate=numpy.random)
+		method=method, update_interval=None, maxcall=max_samples, dlogz=0.5, rstate=numpy.random, callback=dump_callback,
+		**options)
 	endtime = time.time()
 	output_basename = config['output_basename']
-	print 'nestle done lnZ = %(logz).1f +- %(logzerr).1f' % (result)
+	print
+        print 'nestle done lnZ = %(logz).1f +- %(logzerr).1f' % (result)
 
 	if config.get('seed', 0) == 0:
 		import matplotlib.pyplot as plt
@@ -64,18 +81,26 @@ def priortransform(u):
 	return dict(
 		Z_computed = float(result['logz']),
 		Z_computed_err = float(result['logzerr']),
-		niterations = result['ncall'],
+		niterations = result['niter'],
 		duration = endtime - starttime,
 	)
 
 configs = [
 	[
-		#dict(nlive_points=100),
+		dict(nlive_points=100),
 		dict(nlive_points=400),
 		dict(nlive_points=1000),
 	], [
+		dict(method='classic'),
 		dict(method='single'),
 		dict(method='multi'),
+		#dict(method='single-robust'),
+		#dict(method='multi-robust'),
+		#dict(method='single-veryrobust'),
+		#dict(method='multi-veryrobust'),
+		#dict(method='multi-limitedrobust'),
+		#dict(method='multi-simplelimitedrobust'),
+		dict(method='multi-rememberingrobust'),
 	]
 ]
 configs = [dict([[k, v] for d in config for k, v in d.iteritems()]) for config in itertools.product(*configs)]

testsuite/algorithms/runnestle.py

@@ -3,43 +3,77 @@
 """
 import gauss
 import loggamma
-import real.timeseries
+#import real.timeseries
+import real.bexvar
 
 problems = [
 	#gauss.create_problem_gauss(ndim=1, problem_name='gauss1d'),
 	gauss.create_problem_gauss(ndim=3, problem_name='gauss3d'),
 	gauss.create_problem_gauss(ndim=10, problem_name='gauss10d'),
-	gauss.create_problem_vargauss(ndim=10, problem_name='vargauss10d'),
-	gauss.create_problem_vargauss(ndim=20, problem_name='vargauss20d'),
+	#gauss.create_problem_gauss_i(ndim=5, problem_name='gauss_i_5d'),
+	#gauss.create_problem_gauss_i(ndim=10, problem_name='gauss_i_10d'),
+	#gauss.create_problem_vargauss(ndim=5, problem_name='vargauss5d'),
+	#gauss.create_problem_vargauss(ndim=10, problem_name='vargauss10d'),
+	#gauss.create_problem_corrgauss(ndim=5, problem_name='corrgauss5d'),
+	gauss.create_problem_corrgauss(ndim=10, problem_name='corrgauss10d'),
+	#gauss.create_problem_corrgauss(ndim=10, problem_name='corrgauss10d-2', difficulty=2),
+	#loggamma.create_problem_bananas(ndim=6, problem_name='bananas6d'),
+	#loggamma.create_problem_bananas(ndim=10, problem_name='bananas10d'),
+	#loggamma.create_problem_poisson_counts(ndim=3, problem_name='poissoncounts3d'),
+	#loggamma.create_problem_poisson_counts(ndim=10, problem_name='poissoncounts10d'),
+	#gauss.create_problem_vargauss(ndim=20, problem_name='vargauss20d'),
 	#gauss.create_problem_gauss(ndim=50, problem_name='vargauss50d'),
+	gauss.create_problem_gauss_i(ndim=2, problem_name='gauss_i_2d'),
+	gauss.create_problem_gauss_i(ndim=3, problem_name='gauss_i_3d'),
+	gauss.create_problem_gauss_i(ndim=4, problem_name='gauss_i_4d'),
 	gauss.create_problem_shell(ndim=2, problem_name='shell2d'),
 	#gauss.create_problem_shell(ndim=2, problem_name='shell2dtilt10', tilt=10),
 	gauss.create_problem_shell(ndim=10, problem_name='shell10d'),
+	#gauss.create_problem_shell(ndim=10, problem_name='shell10d-thin', width=[0.001/12]*2),
 	gauss.create_problem_gen_gauss_sequence(ndim=7, problem_name='norm_sequence'),
 	gauss.create_problem_complement_gen_gauss_sequence(ndim=7, problem_name='norm_sequence_complement'),
 	gauss.create_problem_eggbox(problem_name='eggbox'),
 	loggamma.create_problem_rosenbrock(ndim=2, problem_name='rosenbrock2d'),
 	loggamma.create_problem_loggamma(ndim=2, problem_name='loggamma2d'),
-	loggamma.create_problem_loggamma(ndim=2, problem_name='loggamma10d'),
-	loggamma.create_problem_funnel(ndim=2, problem_name='funnel2d'),
-	loggamma.create_problem_funnel(ndim=5, problem_name='funnel5d'),
+	loggamma.create_problem_loggamma(ndim=10, problem_name='loggamma10d'),
+	#loggamma.create_problem_funnel(ndim=2, problem_name='funnel2d'),
+	#loggamma.create_problem_funnel(ndim=2, difficulty=0, problem_name='funnel2d-0'),
+	#loggamma.create_problem_funnel(ndim=2, difficulty=1, problem_name='funnel2d-1'),
+	#loggamma.create_problem_funnel(ndim=2, difficulty=2, problem_name='funnel2d-2'),
+	#loggamma.create_problem_funnel(ndim=5, problem_name='funnel5d'),
+	#loggamma.create_problem_funnel(ndim=5, difficulty=0, problem_name='funnel5d-0'),
+	#loggamma.create_problem_funnel(ndim=5, difficulty=1, problem_name='funnel5d-1'),
+	#loggamma.create_problem_funnel(ndim=5, difficulty=2, problem_name='funnel5d-2'),
+	#loggamma.create_problem_funnel(ndim=10, problem_name='funnel10d'),
+	#loggamma.create_problem_funnel(ndim=20, problem_name='funnel20d'),
+	loggamma.create_problem_ffunnel(ndim=2, problem_name='ffunnel2d-3', difficulty=3),
+	loggamma.create_problem_ffunnel(ndim=2, problem_name='ffunnel2d-5', difficulty=5),
+	loggamma.create_problem_ffunnel(ndim=5, problem_name='ffunnel5d-3', difficulty=3),
+	loggamma.create_problem_ffunnel(ndim=5, problem_name='ffunnel5d-5', difficulty=5),
+	loggamma.create_problem_ffunnel(ndim=10, problem_name='ffunnel10d-5', difficulty=5),
+	#real.exvar.create_problem_exvar(ndim=5, problem_name='exvar5d'),
+	real.bexvar.create_problem_bexvar(ndim=10, problem_name='bexvar10d', variance=1),
+	#real.bexvar.create_problem_bexvar(ndim=10, problem_name='bexvar10d-1', variance=-1),
+	#real.exvar.create_problem_exvar(ndim=20, problem_name='exvar20d'),
+	#loggamma.create_problem_ffunnel(ndim=20, problem_name='ffunnel20d', difficulty=5),
 	loggamma.create_problem_spikeslab(problem_name='spikeslab2d_difficulty3', ndim=2, difficulty=3),
 	loggamma.create_problem_spikeslab(problem_name='spikeslab2d_difficulty4', ndim=2, difficulty=4),
 	#loggamma.create_problem_spikeslab(problem_name='spikeslab2d_difficulty6', ndim=2, difficulty=6),
 	#loggamma.create_problem_spikeslab(problem_name='spikeslab5d_difficulty2', ndim=5, difficulty=2),
 	#loggamma.create_problem_spikeslab(problem_name='spikeslab5d_difficulty6', ndim=5, difficulty=6),
 	# the slowest last
-	loggamma.create_problem_loggamma_multimodal(ndim=2, problem_name='loggamma_multimodal2d'),
-	loggamma.create_problem_loggamma_multimodal(ndim=5, problem_name='loggamma_multimodal5d'),
-	loggamma.create_problem_loggamma_multimodal(ndim=10, problem_name='loggamma_multimodal10d'),
-	loggamma.create_problem_eyes(ndim=5, hardness=1, problem_name='eyes5d'),
-	loggamma.create_problem_eyes(ndim=10, hardness=1, problem_name='eyes10d'),
-	loggamma.create_problem_eyes(ndim=5, hardness=5, problem_name='eyes5d-5'),
-	loggamma.create_problem_eyes(ndim=10, hardness=5, problem_name='eyes10d-5'),
-	real.timeseries.create_problem_RVexoplanet(nplanets=0, problem_name='RVexoplanet-0'),
-	real.timeseries.create_problem_RVexoplanet(nplanets=1, problem_name='RVexoplanet-1'),
+	loggamma.create_problem_loggammaI_multimodal(ndim=2, problem_name='loggammaI_multimodal2d'),
+	#loggamma.create_problem_loggammaI_multimodal(ndim=5, problem_name='loggammaI_multimodal5d'),
+	loggamma.create_problem_loggammaI_multimodal(ndim=10, problem_name='loggammaI_multimodal10d'),
+	#loggamma.create_problem_loggammaI_multimodal(ndim=20, problem_name='loggammaI_multimodal20d'),
+	#loggamma.create_problem_eyes(ndim=5, hardness=1, problem_name='eyes5d'),
+	#loggamma.create_problem_eyes(ndim=10, hardness=1, problem_name='eyes10d'),
+	#loggamma.create_problem_eyes(ndim=5, hardness=5, problem_name='eyes5d-5'),
+	#loggamma.create_problem_eyes(ndim=10, hardness=5, problem_name='eyes10d-5'),
+	#real.timeseries.create_problem_RVexoplanet(nplanets=0, problem_name='RVexoplanet-0'),
+	#real.timeseries.create_problem_RVexoplanet(nplanets=1, problem_name='RVexoplanet-1'),
 	#loggamma.create_problem_eyes(ndim=20, hardness=5, problem_name='eyes20d-5'),
-	real.timeseries.create_problem_RVexoplanet(nplanets=2, problem_name='RVexoplanet-2'),
+	#real.timeseries.create_problem_RVexoplanet(nplanets=2, problem_name='RVexoplanet-2'),
 ]
 
 

testsuite/problems/__init__.py

@@ -74,7 +74,7 @@ def draw_constrained(x, Lmin):
 
 def create_problem_gauss_multimodal(**config):
 	ndim = config.get('ndim', 2)
-	assert ndim >= 2
+	assert ndim >= 2, ('ndim must be at least 2')
 
 	def loglikelihood(y):
 		x = numpy.asarray(y)
@@ -109,25 +109,98 @@ def volume_computation(u=None, x=None, L=None):
 
 def create_problem_vargauss(**config):
 	ndim = config.get('ndim', 10)
-	i = numpy.arange(ndim) + 3.0
+	assert ndim > 3, ('ndim must be at least 3')
+	i = numpy.arange(ndim - 3)
 	mu = 0.54321
-	s = 6./((3.*i)**1.4)
+	s = 10**(-((i + 2))*0.5)
 	def loglikelihood(x):
-		z = numpy.asarray(x)
+		z = numpy.asarray(x[3:])
 		return -0.5 * (((z - mu)/s)**2).sum()
 
 	Z_analytic = 0.5 * log(2*pi * s**2).sum()
 	config['loglikelihood'] = loglikelihood
 	config['Z_analytic'] = Z_analytic
 	config['description'] = """Sequence of independent gaussians which become
-	narrower and narrower with the dimension number. Symmetric, %d dimensions.
+	narrower and narrower with the dimension number. 
+	The first three dimensions are unconstrained.
+	Difficult because the parameter sizes are orders of magnitudes different,
+	and region of interest is tiny.
+	Symmetric, %d dimensions. 
+	""" % ndim
+	return config
+
+
+def create_problem_corrgauss(**config):
+	ndim = config.get('ndim', 10)
+	difficulty = config.get('difficulty', 1)
+
+	logmatrix = numpy.zeros((ndim, ndim), dtype=int)
+	matrix = numpy.zeros((ndim, ndim))
+
+	for i in range(ndim):
+		for j in range(i+1):
+			logmatrix[i,j] = (((j*i) % (5+i)) // 2) % 5
+			logmatrix[j,i] = logmatrix[i,j]
+	
+	for i in range(ndim):
+		for j in range(i+1):
+			scalei = numpy.exp(-logmatrix[i,i] - 2)
+			scalej = numpy.exp(-logmatrix[j,j] - 2)
+			if i == j:
+				matrix[i,i] = scalei*scalej
+			else:
+				matrix[i,j] = (1 - numpy.exp(-logmatrix[i,j]*difficulty))*scalei*scalej
+				matrix[j,i] = (1 - numpy.exp(-logmatrix[j,i]*difficulty))*scalei*scalej
+	# ensure positive semi-definite
+	_, matrix = scipy.linalg.polar(matrix)
+	matrix = numpy.matrix(matrix)
+	invmatrix = numpy.linalg.inv(matrix)
+	mu = 0.5 + numpy.zeros(ndim)
+	assert numpy.linalg.det(matrix) > 0, numpy.linalg.det(matrix)
+	norm = -0.5 * log(numpy.linalg.det(2 * pi * matrix))
+	def loglikelihood(x):
+		z = numpy.asarray(x) - mu
+		return -0.5 * numpy.dot(numpy.dot(z, invmatrix), z.T) + norm
+	
+	Z_analytic = 0
+	config['loglikelihood'] = loglikelihood
+	config['Z_analytic'] = Z_analytic
+	config['description'] = """Correlated gaussian likelihood in %d dimensions. 
+	-log(1 - correlation matrix) is: <pre>%s</pre>
+	Centred at 0.5.
+	""" % (ndim, logmatrix*difficulty)
+	return config
+
+
+
+def create_problem_gauss_i(**config):
+	ndim = config.get('ndim', 10)
+	mu1 = 0.654321
+	mu2 = 0.123456
+	s1 = 0.1
+	s2 = 0.01
+	norm1 = -0.5 * log(2*pi*s1**2) * ndim
+	norm2 = -0.5 * log(2*pi*s2**2) * ndim + log(100)
+	def loglikelihood(x):
+		z = numpy.asarray(x)
+		return numpy.logaddexp(norm1 - 0.5 * ((((z - mu1)/s1)**2)).sum(), 
+			norm2 - 0.5 * (((z - mu2)/s2)**2).sum())
+	
+	Z_analytic = log(101)
+	config['loglikelihood'] = loglikelihood
+	config['Z_analytic'] = Z_analytic
+	config['description'] = """Two gaussian solutions, one large 0.654321+-0.1, one small 0.123456+-0.01, in %d dimensions.
+	
+	<p>This function may be difficult if the smaller solution is lost,
+	because by area it is small; but later it turns out to be the important one.
 	""" % ndim
 	return config
 
 
+
 def create_problem_halfgauss(**config):
 	ndim = config.get('ndim', 2)
-	assert ndim >= 2
+	assert ndim >= 2, ('ndim must be at least 2')
 
 	def loglikelihood(y):
 		x = numpy.asarray(y)
@@ -295,7 +368,7 @@ def loglikelihood(x):
 
 def create_problem_pyramid(**config):
 	ndim = config.get('ndim', 20)
-	assert ndim >= 2
+	assert ndim >= 2, ('ndim must be at least 2')
 
 	def loglikelihood(x):
 		x = numpy.asarray(x)