def f3d_to_fcylavg(f, x, y, z, bins=None, log=False):
    """Cylindrically average a function defined on a 3D grid

    Parameters
    ----------
    f :
    x :
    y :
    z :
    bins : 
        Number of bins, in the form (ks, kpar)

    Returns
    -------
    
    """
    if bins == None:
        # Count the maximum number of gridpts outward from origin in x, y, or z, then subtract 1
        bins_prp = max([ max( np.count_nonzero(a>0), np.count_nonzero(a<0) ) for a in (x,y) ]) - 1
        bins_par = max( np.count_nonzero(z>0), np.count_nonzero(z<0) ) - 1.
        bins = (bins_prp, bins_par)

    # Get (r,z) cylindrical bin edges from (x,y,z) cartesian grid
    rcylbins    = xyz_to_rcylbins(x,y,z, bins=bins, log=log)

    xx,yy,zz = np.meshgrid(x, y, z, indexing='ij')
    rr  = np.sqrt( xx**2 + yy**2 )                                  # 3D grid of r (cylindrical radius from origin)
    nn  = np.histogram2d(rr.ravel(), zz.ravel(), bins=rcylbins)[0]  # Number of cells used to compute average
    favg= np.histogram2d(rr.ravel(), zz.ravel(), bins=rcylbins, weights=f.ravel())[0] / nn

    rmid_prp    = _bin_midpoints(rcylbins[0])
    rmid_par    = _bin_midpoints(rcylbins[1])

    return rmid_prp, rmid_par, favg

 def test_binparameter_combination(self):
     x = array(
         [0, 0.09207008,  0.64575234,  0.12875982,  0.47390599,  
          0.59944483, 1])
     y = array(
         [0, 0.14344267,  0.48988575,  0.30558665,  0.44700682,  
          0.15886423, 1])
     edges = (0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1)
     H, xe, ye = histogram2d(x, y, (edges, 4))
     answer = array(
         [[ 2.,  0.,  0.,  0.],
          [ 0.,  1.,  0.,  0.],
          [ 0.,  0.,  0.,  0.],
          [ 0.,  0.,  0.,  0.],
          [ 0.,  1.,  0.,  0.],
          [ 1.,  0.,  0.,  0.],
          [ 0.,  1.,  0.,  0.],
          [ 0.,  0.,  0.,  0.],
          [ 0.,  0.,  0.,  0.],
          [ 0.,  0.,  0.,  1.]])
     assert_array_equal(H, answer)
     assert_array_equal(ye, array([0., 0.25, 0.5, 0.75, 1]))
     H, xe, ye = histogram2d(x, y, (4, edges))
     answer = array(
         [[ 1.,  1.,  0.,  1.,  0.,  0.,  0.,  0.,  0.,  0.],
          [ 0.,  0.,  0.,  0.,  1.,  0.,  0.,  0.,  0.,  0.],
          [ 0.,  1.,  0.,  0.,  1.,  0.,  0.,  0.,  0.,  0.],
          [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  1.]])
     assert_array_equal(H, answer)
     assert_array_equal(xe, array([0., 0.25, 0.5, 0.75, 1]))

 def test_simple(self):
     x = array(
         [0.41702200, 0.72032449, 1.1437481e-4, 0.302332573, 0.146755891])
     y = array(
         [0.09233859, 0.18626021, 0.34556073, 0.39676747, 0.53881673])
     xedges = np.linspace(0, 1, 10)
     yedges = np.linspace(0, 1, 10)
     H = histogram2d(x, y, (xedges, yedges))[0]
     answer = array(
         [[0, 0, 0, 1, 0, 0, 0, 0, 0],
          [0, 0, 0, 0, 0, 0, 1, 0, 0],
          [0, 0, 0, 0, 0, 0, 0, 0, 0],
          [1, 0, 1, 0, 0, 0, 0, 0, 0],
          [0, 1, 0, 0, 0, 0, 0, 0, 0],
          [0, 0, 0, 0, 0, 0, 0, 0, 0],
          [0, 0, 0, 0, 0, 0, 0, 0, 0],
          [0, 0, 0, 0, 0, 0, 0, 0, 0],
          [0, 0, 0, 0, 0, 0, 0, 0, 0]])
     assert_array_equal(H.T, answer)
     H = histogram2d(x, y, xedges)[0]
     assert_array_equal(H.T, answer)
     H, xedges, yedges = histogram2d(list(range(10)), list(range(10)))
     assert_array_equal(H, eye(10, 10))
     assert_array_equal(xedges, np.linspace(0, 9, 11))
     assert_array_equal(yedges, np.linspace(0, 9, 11))

def plot_predicted_vs_observed_pks(obs_score, pred_score, ofname):
    plt.figure()
    heatmap, xedges, yedges = numpy.histogram2d(
        rankdata(-obs_score, method='ordinal'), 
        rankdata(pred_score, method='ordinal'), 
        bins=20)
    heatmap, xedges, yedges = numpy.histogram2d(
        numpy.clip(-numpy.log(1+obs_score), -0.1, 0),
        numpy.clip(numpy.log(1+pred_score), 0, 0.1), 
        bins=100)
    #heatmap, xedges, yedges = numpy.histogram2d(
    #    numpy.clip(-numpy.log(1+data['ATAC_mean']), -0.1, 0),
    #    numpy.clip(numpy.log(y), 0, 0.1), 
    #    bins=100)

    #heatmap, xedges, yedges = numpy.histogram2d(
    #    rankdata(-data['ATAC_mean'], method='average'), 
    #    rankdata(obs_score, method='average'), 
    #    bins=20)
    extent = [xedges[0], xedges[-1], yedges[0], yedges[-1]]
    plt.clf()
    plt.imshow(heatmap, extent=extent)
    #plt.scatter(rankdata(obs_score, method='ordinal'), 
    #            rankdata(pred_score, method='ordinal'))
    plt.savefig(ofname)
    plt.close()
    return

def radprof3D(field, xx, yy, z, center, nbins):
    """
    compute radial profile from the 'center' point for 3D field
    assumes square, periodic domain, given by x,y.
    
    """
    #nx = len(xx[:,0])
    #ny = len(yy[:,0])
    nz = z.size
    width = xx[0,-1]
    delx = np.abs(xx - center[0])
    dely = np.abs(yy - center[1])
    delx[delx > width/2] = width - delx[delx > width/2]
    dely[dely > width/2] = width - dely[dely > width/2]
    r = np.sqrt(np.power(delx, 2) + np.power(dely, 2))
    
    r = r.ravel()
    field = field.reshape(nz, r.size)
    #r = np.repeat(r[np.newaxis,:], z.size, axis=0)
    rr, zz = np.meshgrid(r, z)
  
    fieldsums, redges, zedges = np.histogram2d(rr.ravel(), zz.ravel(), bins=nbins, weights=field.ravel())
    nr, redges, zedges = np.histogram2d(rr.ravel(), zz.ravel(), bins=nbins)
    
    fieldmeans = fieldsums/nr
    
    return (redges, zedges, fieldmeans)

def bootstrap_estimate(RA, DEC, gamma, RAbins, DECbins, Nbootstrap=1000):
    N = len(gamma)
    
    gamma_sum = 0
    gamma2_sum = 0

    for i in range(Nbootstrap):
        ind = np.random.randint(N, size=N)

        RAi = RA[ind]
        DECi = DEC[ind]
        gammai = gamma[ind]

        count = np.histogram2d(RAi, DECi, (RAbins, DECbins))
        shear = np.histogram2d(RAi, DECi, (RAbins, DECbins), weights=gammai)

        mask = (count == 0)
        count[mask] = 1

        shear /= count

        gamma_sum += shear
        gamma2_sum += abs(shear) ** 2

    gamma_sum /= Nbootstrap
    gamma2_sum /= Nbootstrap

    d2gamma = gamma2_sum - abs(gamma_sum ** 2)

  def cached(data):
    cols = datatable.get_cols(*markers)
    if not range:
      fixed_range = [
          min(cols[0]),
          min(cols[1]),
          max(cols[0]),
          max(cols[1])]    
    else:
      fixed_range = range
      
    hist, data.x_edges, data.y_edges = np.histogram2d(
        cols[0], 
        cols[1], 
        [
            np.r_[fixed_range[0]:fixed_range[2]:no_bins], 
            np.r_[fixed_range[1]:fixed_range[3]:no_bins]])
    data.final_hist = np.sign(np.subtract(
        np.clip(np.abs(hist), min_cells_per_bin, np.inf),
        min_cells_per_bin))

    if color_marker:
      data.is_colored = True
      weights = datatable.get_cols(color_marker)[0]     
      weighted_hist, x_edges, y_edges = np.histogram2d(
          cols[0], 
          cols[1], 
          [
              np.r_[fixed_range[0]:fixed_range[2]:no_bins], 
              np.r_[fixed_range[1]:fixed_range[3]:no_bins]], None, False, weights)
      data.colored_hist = np.multiply(
          np.true_divide(weighted_hist, hist), data.final_hist)
    else:
      data.is_colored = False

def conn_tests(solver, save_format):
    print '>>>', datetime.now(), 'Begin conn tests'
    solver.net.set_phase_test()
    solver.test_nets[0].share_with(solver.net)
    net = solver.test_nets[0]
    save_dir = save_format.format(solver.iter)
    os.mkdir(save_dir)
    hist = np.zeros((2, 2))
    for fname in fnames_val:
        net.forward()
        im = Image.fromarray(net.blobs['left-conn'].data[0].argmax(0)
                .astype(np.uint8), mode='P')
        im.save(os.path.join(save_dir, fname + '-L.png'))
        im = Image.fromarray(net.blobs['top-conn'].data[0].argmax(0)
                .astype(np.uint8), mode='P')
        im.save(os.path.join(save_dir, fname + '-T.png'))
        h, _ , _ = np.histogram2d(net.blobs['left-gt'].data[0, 0].flatten(),
                net.blobs['left-conn'].data[0].argmax(0).flatten(),
                bins = 2, range=[[0, 2], [0, 2]])
        hist += h
        h, _ , _ = np.histogram2d(net.blobs['top-gt'].data[0, 0].flatten(),
                net.blobs['top-conn'].data[0].argmax(0).flatten(),
                bins = 2, range=[[0, 2], [0, 2]])
        hist += h
    print '>>>', datetime.now(), 'Iteration', solver.iter, 'overall accuracy', \
        np.diag(hist).sum() / hist.sum()
    print '>>>', datetime.now(), 'Iteration', solver.iter, 'total IU', \
            hist[1, 1] / (hist[0, 1] + hist[1, 0] + hist[1, 1])
    solver.net.set_phase_train()

def contourplotgeneric(vx,vy,wgt=None,ax=None,wgtopt=1,\
               nxbins=25,nybins=25,flist=[0.683,0.954],
               linestyle='-',color='k',linewidths=3):
  """
  Set wgtopt = 1 to weight the points according to the field 'weight'
  in the chain.  Set wgtopt = 0 to weight the points equally.
  flist is fraction of chain you want to enclose.
  """
  if wgt is None:
    wgtopt = 0
  if wgtopt == 1:  
    H, xedges, yedges = np.histogram2d(vx,vy,[nxbins,nybins],weights=wgt)
  else:
    H, xedges, yedges = np.histogram2d(vx,vy,[nxbins,nybins])

  Hflip = np.zeros([nybins,nxbins])  # for some idiotic reason, contour progam bins are flipped.
  xcen = (xedges[1:] + xedges[:-1])*0.5
  ycen = (yedges[1:] + yedges[:-1])*0.5
  for x in range(nxbins):
    for y in range(nybins):
      Hflip[y,x] = H[x,y]

  clist = getclevels(Hflip,flist)
  if ax is None:
    f = plt.figure()
    ax = f.add_subplot(111)
    ax.contour(xcen,ycen,Hflip,levels=clist,linestyles=linestyle,colors=color,linewidths=linewidths)
    return f, ax
  else:
    ax.contour(xcen,ycen,Hflip,levels=clist,linestyles=linestyle,colors=color,linewidths=linewidths)

def dualHistogram2D(X1,Y1,X2,Y2,res=16,title1='Pairs',title2='Isolated',xLabel= "Proj. Dis (Mpc) binned",yLabel= "dv (km/s) binned"):
    
    hist1,xedges1, yedges1= numpy.histogram2d(X1,Y1,bins=res)
    hist2,xedges2, yedges2= numpy.histogram2d(X2,Y2,bins=res)
    #print xedges1," ",yedges1
    #create bar graphs   
    fig = plt.figure()
    
    #set ranges
    #extent1 = [yedges1[0], yedges1[-1], xedges1[-1], xedges1[0]]
    plot1 = fig.add_subplot(121)
    im1 = plot1.imshow(hist1,interpolation='nearest')
    plot1.set_ylabel(yLabel)
    plot1.set_xlabel(xLabel)
    plot1.set_title(title1)
    fig.colorbar(im1) 
    
    #set ranges
    extent2 = [yedges2[0], yedges2[-1], xedges2[-1], xedges2[0]]
    plot2 = fig.add_subplot(122)
    im2 =  plot2.imshow(hist2,interpolation='nearest')
    plot2.set_ylabel(yLabel)
    plot2.set_xlabel(xLabel)
    plot2.set_title(title2)
    fig.colorbar(im2)
    
    
    plt.show()

def get_corr_map(coo1, coo2, skypos, skyrange):
  imsz = imagetools.deg2pix(skypos, skyrange, 0.001)
  count = np.zeros(imsz)
  print(imsz)
  co_rel = np.array([[0,0]])
  len1 = coo1.shape[0]
  len2 = coo2.shape[0]
  print(len1,len2)
  wcs = imagetools.define_wcs(skypos,skyrange,width=False,height=False,verbose=0,pixsz=0.001)
  if len2>len1:
    for i in range(len2):
      print(i)
      co_rel = np.concatenate((co_rel, np.roll(coo2, i, axis=0)[0:len1,:]-coo1), axis = 0)
      if (i+1)%200 == 0:
        foc = wcs.sip_pix2foc(wcs.wcs_world2pix(co_rel[1:],1),1)
        H,xedges,yedges=np.histogram2d(foc[:,1]-0.5, foc[:,0]-0.5,\
                                   bins=imsz, range=([ [0,imsz[0]],[0,imsz[1]] ]))
        count += H
        co_rel = np.array([[0, 0]])
  else:
    for i in range(len1):
      print(i)
      co_rel = np.concatenate((co_rel, coo2-np.roll(coo1, i, axis=0)[0:len2,:]), axis = 0)
      if (i+1)%200 == 0:
        foc = wcs.sip_pix2foc(wcs.wcs_world2pix(co_rel[1:],1),1)
        H,xedges,yedges=np.histogram2d(foc[:,1]-0.5, foc[:,0]-0.5,\
                                   bins=imsz, range=([ [0,imsz[0]],[0,imsz[1]] ]))
        count += H
        co_rel = np.array([[0, 0]])
  if co_rel.shape[0]>1:
    foc = wcs.sip_pix2foc(wcs.wcs_world2pix(co_rel[1:],1),1)
    H,xedges,yedges=np.histogram2d(foc[:,1]-0.5, foc[:,0]-0.5,\
                               bins=imsz, range=([ [0,imsz[0]],[0,imsz[1]] ]))
    count += H
  return count

def hist2d(x,y,nBinsX,nBinsY,rngX=None,rngY=None):
    if rngX == None and rngY == None:
        h2d, xp, yp = np.histogram2d(y,x,bins=(nBinsY,nBinsX), normed = True)
    else:
        h2d, xp, yp = np.histogram2d(y,x,bins=(nBinsY,nBinsX),range=[rngY,rngX], normed = True)
    extent = [yp[0],yp[-1],xp[0],xp[-1]]
    return h2d, extent

    def plot2D(self,i,j):
      
        fig, axes = plt.subplots(1, 1)
     
        array = []
        labels=[]
        x = []
        y = []
        
        self.LastChain.setSample(0)
        while self.LastChain.getNextPoint():
	    x.append( self.Parameters.getVariable(j).getValue() )
	    y.append( self.Parameters.getVariable(i).getValue() )
     
        nbins = 50
        H, xedges, yedges = np.histogram2d(x,y,bins=nbins)
        while(H.max()>100):
            nbins = nbins + 1 
            H, xedges, yedges = np.histogram2d(x,y,bins=nbins)
  
        H = np.rot90(H)
        H = np.flipud(H)
        Hmasked = np.ma.masked_where(H==0,H) 
        axes.pcolormesh(xedges,yedges,Hmasked)
        axes.locator_params(nbins = 3)
        
        axes.set_xlabel("x")
        axes.set_ylabel("y")
        fig.show()

    def computeDensityFieldForHaloNumber(self, path_to_DF, path_to_RS, outputFile, gridSize=2048, subgridSize = 256 ) :
		"""
        Extracts the distribution of quantity 'name' out of all snapshots of the Multidark simulation.        
        :param path_to_DF: path to the density field file
        :param path_to_RS: path to the rockstar halo catalog file
        :param outputFile: where the histogram iswritten
        :param bins: binning scheme to compute the historgram.
		:param gridSize: grid size from the density field
		:param subgridSize: grid size to compute histograms on and write outputs.
        """		
        #In each cell, average the number of counts in the RS file : Ncen Nsat
		dL = 1000./gridSize
		NperBatch = subgridSize**3
		Ntotal = gridSize**3
		Nratio = Ntotal / NperBatch

		hf=open(path_to_RS,'r')
		DEFINE x, y, z, c_o_s, 
		REWRITE halos of interest ?
		
		f=open(path_to_DF,'r')
		out = n.empty( (subgridSize**3, 3) )
		count = 0
		countOut = 0
		for kk in n.arange(gridSize):
			for jj in n.arange(gridSize):
				for ii in n.arange(gridSize):
					sel =( x > dL*(0.5 + ii) ) & ( x < dL*(0.5 + ii + 1) ) & ( y > dL*(0.5 + jj) ) & ( y < dL*(0.5 + jj + 1) ) & ( x > dL*(0.5 + kz) ) & ( z < dL*(0.5 + kk + 1) )
					Nhalos = len(sel.nonzero()[0])
					selCen =  (sel) & (CONDITION_CEN)
					NCentrals = len(selCen.nonzero()[0])
					deltaValue = n.fromfile(f,dtype="float64",count=1)
					out[count] = n.array([deltaValue, Nhalos, Ncentrals])
					if count == subgridSize**3 :
						dataAB_tot = n.histogram2d(n.log10(out.T[0]), n.log10(out.T[1]) ,bins = [binsA,binsB])
						dataAB_cen = n.histogram2d(n.log10(out.T[0]), n.log10(out.T[2]) ,bins = [binsA,binsB])
						f = open(outputFile + "_" +str(countOut)+".pkl" ,'w')
						cPickle.dump([binsA,binsB,dataAB_tot, dataAB_cen],f)
						f.close()

						out = n.empty( (subgridSize**3, 3) )
						countOut +=1 
					
					count += 1
					
		#GATHER RESULTS
		fileList = glob.glob(outputFile + "_*.pkl")
		out_all = n.empty( (Nratio, len(binsA)-1, len(binsB)-1) )
		out_cen = n.empty( (Nratio, len(binsA)-1, len(binsB)-1) )
		for ii, el in enumerate(fileList):
			f = open(el ,'r')
			binsA,binsB,dataAB_tot, dataAB_cen = cPickle.dump([binsA,binsB,dataAB_tot, dataAB_cen],f)
			f.close()
			out_all[ii] = dataAB_tot
			out_cen[ii] = dataAB_cen

		
		f = open(outputFile + "_all.pkl" ,'w')
		cPickle.dump([binsA,binsB,n.sum(out_all, xis=0), n.sum(out_cen, xis=0)],f)
		f.close()

def basic_imshow(ax,xfield,yfield1,yfield2,plot_type='holding',xkey='temp'):
	if plot_type=='emis':
		yfield1,yfield2,xfield = yfield1.flat,yfield2.flat,xfield.flat
		idx1 = np.where(yfield1 > 1.)[0]
		yfield1,yfield2,xfield = yfield1[idx1],yfield2[idx1],xfield[idx1]
		idx2 = np.where(yfield2 > 1.)[0]
		yfield1,yfield2,xfield = yfield1[idx2],yfield2[idx2],xfield[idx2]
		ion_frac = np.log10(yfield1/yfield2)
	if plot_type=='coldens':
		yfield1,yfield2,xfield = yfield1.flat,yfield2.flat,xfield.flat
                idx1 = np.where(np.log10(yfield1) > 12.)[0]
                yfield1,yfield2,xfield = yfield1[idx1],yfield2[idx1],xfield[idx1]
                idx2 = np.where(np.log10(yfield2) > 12.)[0]
                yfield1,yfield2,xfield = yfield1[idx2],yfield2[idx2],xfield[idx2]
                ion_frac = np.log10(yfield1/yfield2)
	else:
		ion_frac = yfield1/yfield2
		xfield = xfield.flat
		ion_frac = np.log10(ion_frac.flat)
	if xkey == 'temp':
		hist,xedges,yedges = np.histogram2d(np.log10(xfield),ion_frac,range=[[4,6.5],[-3,3]],bins=100)
		extent = [yedges[0], yedges[-1], xedges[0], xedges[-1]]
        	ax.imshow(np.log10(hist.T),origin='lower',extent=[4,6.5,-3,3],interpolation='nearest',cmap='jet',vmin=-3.,vmax=3.) 
	if xkey == 'dens':
		hist,xedges,yedges = np.histogram2d(np.log10(xfield),ion_frac,range=[[-5,-1],[-4,2]],bins=100)
		extent = [yedges[0], yedges[-1], xedges[0], xedges[-1]]
        	ax.imshow(np.log10(hist.T),origin='lower',extent=[-5,-1,-4,2],interpolation='nearest',cmap='jet',vmin=-3.,vmax=3.)
	if xkey == 'NH':
		hist,xedges,yedges = np.histogram2d(np.log10(xfield),ion_frac,range=[[18,22],[-6,3]],bins=100)
		extent = [yedges[0], yedges[-1], xedges[0], xedges[-1]]
		ax.imshow(np.log10(hist.T),origin='lower',extent=[18,22,-6,3],interpolation='nearest',cmap='jet',vmin=-3.,vmax=3.)
	ax.set_aspect('auto')
	return

def contourplot(x,y,wgt=None,ax=None,\
                 nxbins=25,nybins=25,flist=[0.683,0.954],
                 linestyle='-',color='k',linewidths=3):
  """
  Input optional weight to make a weighted 2d histogram.
  flist is fraction of points you want to enclose.
  """
  if wgt is None:
    H, xedges, yedges = np.histogram2d(x,y,[nxbins,nybins])
  else:
    H, xedges, yedges = np.histogram2d(x,y,[nxbins,nybins],weights=wgt)
  Hflip = np.zeros([nybins,nxbins])  # for some idiotic reason, contour progam bins are flipped.
  xcen = (xedges[1:] + xedges[:-1])*0.5
  ycen = (yedges[1:] + yedges[:-1])*0.5
  for x in range(nxbins):
    for y in range(nybins):
      Hflip[y,x] = H[x,y]

  clist = getclevels(Hflip,flist)
  if ax is None:
    f = plt.figure()
    ax = f.add_subplot(111)
    ax.contour(xcen,ycen,Hflip,levels=clist,linestyles=linestyle,colors=color,linewidths=linewidths)
    return f, ax
  else:
    ax.contour(xcen,ycen,Hflip,levels=clist,linestyles=linestyle,colors=color,linewidths=linewidths)    

def Get_TProfile(ar1, ar2, nbBin1, we2):
	'''
		Create a plot similar to a ROOT TProfile
	'''

	number, axisX, axisY = numpy.histogram2d(ar1, ar2, (nbBin1,1))
	weight, axisX, axisY = numpy.histogram2d(ar1, ar2, (nbBin1,1), weights=we2)
        mean,   axisX, axisY = numpy.histogram2d(ar1, ar2, (nbBin1,1), weights=we2*ar2)
        err,    axisX, axisY = numpy.histogram2d(ar1, ar2, (nbBin1,1), weights=we2*(ar2**2.))

	### find the axis X
        axisX = numpy.array([ axisX[i]+(axisX[i+1]-axisX[i])/2. for i in range(0,axisX.size-1) ])
	
	### Get only not empty bins
	bool_number = (number[:,0]>1)
	
	axisX  = axisX[bool_number]
	number = number[:,0][bool_number]
	weight = weight[:,0][bool_number]
	mean   = mean[:,0][bool_number]
	err    = err[:,0][bool_number]

	mean  = mean/weight
	err   = numpy.sqrt((err/weight-mean**2.)/number)

	return axisX, mean, err, number

    def convert_3dps_to_2dps(self, k_edges_p, k_edges_v):
        
        k_bin_x = 2. * np.pi * np.fft.fftshift(np.fft.fftfreq(self.boxshape[0],
                                                              self.boxunit))
        k_bin_y = 2. * np.pi * np.fft.fftshift(np.fft.fftfreq(self.boxshape[1],
                                                              self.boxunit))
        k_bin_z = 2. * np.pi * np.fft.fftshift(np.fft.fftfreq(self.boxshape[2],
                                                              self.boxunit))

        k_bin_p = np.abs(k_bin_x)
        k_bin_v = np.sqrt( (k_bin_y**2)[:,None] + (k_bin_z**2)[None,:] )

        k_bin_2d = np.zeros(shape=(2,)+ k_bin_p.shape + k_bin_v.shape)
        k_bin_2d[0] = k_bin_p[:, None, None]
        k_bin_2d[1] = k_bin_v[None, :, :]

        ones_weight = np.ones_like(self.ps_3d)
        kn_2d, xedges, yedges = np.histogram2d(k_bin_2d[0].flatten(), 
                                               k_bin_2d[1].flatten(),
                                               (k_edges_p, k_edges_v),
                                               weights=ones_weight.flatten())

        ps_2d, xedges, yedges = np.histogram2d(k_bin_2d[0].flatten(), 
                                               k_bin_2d[1].flatten(), 
                                               (k_edges_p, k_edges_v),
                                               weights=self.ps_3d.flatten())

        kn_2d[kn_2d==0] = np.inf
        ps_2d /= kn_2d
        kn_2d[kn_2d==np.inf] = 0 

        self.kn_2d = kn_2d
        self.ps_2d = ps_2d

def create_contingency_table ( data, dico, x, y, z ):
    # détermination de la taille de z
    size_z = 1
    offset_z = np.zeros ( len ( z ) )
    j = 0
    for i in z:
        offset_z[j] = size_z      
        size_z *= len ( dico[i] )
        j += 1

    # création du tableau de contingence
    res = np.zeros ( size_z, dtype = object )

    # remplissage du tableau de contingence
    if size_z != 1:
        z_values = np.apply_along_axis ( lambda val_z : val_z.dot ( offset_z ),
                                         1, data[z,:].T )
        i = 0
        while i < size_z:
            indices, = np.where ( z_values == i )
            a,b,c = np.histogram2d ( data[x,indices], data[y,indices],
                                     bins = [ len ( dico[x] ), len (dico[y] ) ] )
            res[i] = ( indices.size, a )
            i += 1
    else:
        a,b,c = np.histogram2d ( data[x,:], data[y,:],
                                 bins = [ len ( dico[x] ), len (dico[y] ) ] )
        res[0] = ( data.shape[1], a )
    return res

def sample_agency_magnitude_pairs(data, xbins, ybins, number_samples=1):
    """

    """
    keys = data.keys()
    n_data = len(data[keys[0]])
    if not number_samples or (number_samples == 1):
        # Only one sample, return simple histogram
        #print xbins, ybins
        return np.histogram2d(np.around(data[keys[0]], 2),
                              np.around(data[keys[2]], 2),
                              bins=[xbins, ybins])[0]
    elif (np.max(data[keys[1]]) < 1E-15) and (np.max(data[keys[3]]) < 1E-15):
        # No uncertainty on magnitudes
        return np.histogram2d(np.around(data[keys[0]], 2),
                              np.around(data[keys[2]], 2),
                              bins=[xbins, ybins])[0]
    else:
        counter = np.zeros([len(xbins) - 1, len(ybins) - 1])
        for i in xrange(number_samples):
            # Sample data sets
            data_x = data[keys[0]] + data[keys[1]] * np.random.normal(0., 1.,
                                                                      n_data)
            data_y = data[keys[2]] + data[keys[3]] * np.random.normal(0., 1.,
                                                                      n_data)
            counter += np.histogram2d(data_x, data_y, bins=[xbins, ybins])[0]
        
        return counter / float(number_samples)