def bcesp(y1,y1err,y2,y2err,cerr,nsim=10000):
	"""
Parallel implementation of the BCES with bootstrapping.
Divide the bootstraps equally among the threads (cores) of
the machine. It will automatically detect the number of
cores available.

Usage:

>>> a,b,aerr,berr,covab=bcesp(x,xerr,y,yerr,cov,nsim)

:param x,y: data
:param xerr,yerr: measurement errors affecting x and y
:param cov: covariance between the measurement errors (all are arrays)
:param nsim: number of Monte Carlo simulations (bootstraps)

:returns: a,b - best-fit parameters a,b of the linear regression 
:returns: aerr,berr - the standard deviations in a,b
:returns: covab - the covariance between a and b (e.g. for plotting confidence bands)

.. seealso:: Check out ~/work/projects/playground/parallel python/bcesp.py for the original, testing, code. I deleted some line from there to make the "production" version.

* v1 Mar 2012: serial version ported from bces_regress.f. Added covariance output.
* v2 May 3rd 2012: parallel version ported from nemmen.bcesboot.

.. codeauthor: Rodrigo Nemmen, http://goo.gl/8S1Oo
	"""	
	import time	# for benchmarking
	import multiprocessing
	
	print "BCES,", nsim,"trials... ",
	tic=time.time()
	
	# Find out number of cores available
	ncores=multiprocessing.cpu_count()
	# We will divide the processing into how many parts?
	n=2*ncores
	
	"""
	Must create lists that will be distributed among the many
	cores with structure 
	core1 <- [y1,y1err,y2,y2err,cerr,nsim/n]
	core2 <- [y1,y1err,y2,y2err,cerr,nsim/n]
	etc...
	"""
	pargs=[]	# this is a list of lists!
	for i in range(n):
		pargs.append([y1,y1err,y2,y2err,cerr,nsim/n])
	
	# Initializes the parallel engine
	pool = multiprocessing.Pool(processes=ncores)	# multiprocessing package
			
	"""
	Each core processes ab(input)
		return matrixes Am,Bm with the results of nsim/n
		presult[i][0] = Am with nsim/n lines
		presult[i][1] = Bm with nsim/n lines
	"""
	presult=pool.map(ab, pargs)	# multiprocessing
	pool.close()	# close the parallel engine
	
	# vstack the matrixes processed from all cores
	i=0
	for m in presult:
		if i==0:
			# Initialize the matrixes
			am,bm=m[0].copy(),m[1].copy()
		else: 
			am=numpy.vstack((am,m[0]))
			bm=numpy.vstack((bm,m[1]))
		i=i+1
	
	# Computes the bootstrapping results on the stacked matrixes
	a=numpy.array([ am[:,0].mean(),am[:,1].mean(),am[:,2].mean(),am[:,3].mean() ])
	b=numpy.array([ bm[:,0].mean(),bm[:,1].mean(),bm[:,2].mean(),bm[:,3].mean() ])

	# Error from unbiased sample variances
	erra,errb,covab=numpy.zeros(4),numpy.zeros(4),numpy.zeros(4)
	for i in range(4):
		erra[i]=numpy.sqrt( 1./(nsim-1) * ( numpy.sum(am[:,i]**2)-nsim*(am[:,i].mean())**2 ))
		errb[i]=numpy.sqrt( 1./(nsim-1) * ( numpy.sum(bm[:,i]**2)-nsim*(bm[:,i].mean())**2 ))
		covab[i]=1./(nsim-1) * ( numpy.sum(am[:,i]*bm[:,i])-nsim*am[:,i].mean()*bm[:,i].mean() )
	
	print "%f s" % (time.time() - tic)	
	
	return a,b,erra,errb,covab
Parallel implementation of the BCES with bootstrapping.
Divide the bootstraps equally among the threads (cores) of
the machine. It will automatically detect the number of
cores available.

Usage:

>>> a,b,aerr,berr,covab=bcesp(x,xerr,y,yerr,cov,nsim)

:param x,y: data
:param xerr,yerr: measurement errors affecting x and y
:param cov: covariance between the measurement errors (all are arrays)
:param nsim: number of Monte Carlo simulations (bootstraps)

:returns: a,b - best-fit parameters a,b of the linear regression 
:returns: aerr,berr - the standard deviations in a,b
:returns: covab - the covariance between a and b (e.g. for plotting confidence bands)

.. seealso:: Check out ~/work/projects/playground/parallel python/bcesp.py for the original, testing, code. I deleted some line from there to make the "production" version.

* v1 Mar 2012: serial version ported from bces_regress.f. Added covariance output.
* v2 May 3rd 2012: parallel version ported from nemmen.bcesboot.

.. codeauthor: Rodrigo Nemmen, http://goo.gl/8S1Oo