How to use the treecorr.GGCorrelation function in TreeCorr

config['ng_statistic'] = 'simple'
    config['nn_statistic'] = 'simple'  # For the later Norm tests
    treecorr.corr2(config)
    corr2_output = np.genfromtxt(os.path.join('output','ng_nmap.out'), names=True)
    np.testing.assert_allclose(corr2_output['NMap'], nmap, rtol=1.e-3)
    np.testing.assert_allclose(corr2_output['NMx'], nmx, atol=1.e-3)
    np.testing.assert_allclose(corr2_output['sig_nmap'], np.sqrt(varnmap), rtol=1.e-3)

    # And check the norm output file, which adds a few columns
    dd = treecorr.NNCorrelation(bin_size=0.1, min_sep=0.5, max_sep=40., sep_units='arcmin',
                                verbose=1)
    dd.process(lens_cat)
    rr = treecorr.NNCorrelation(bin_size=0.1, min_sep=0.5, max_sep=40., sep_units='arcmin',
                                verbose=1)
    rr.process(rand_cat)
    gg = treecorr.GGCorrelation(bin_size=0.1, min_sep=0.5, max_sep=40., sep_units='arcmin',
                                verbose=1)
    gg.process(source_cat)
    napsq, varnap = dd.calculateNapSq(rr=rr)
    mapsq, mapsq_im, mxsq, mxsq_im, varmap = gg.calculateMapSq(m2_uform='Crittenden')
    nmap_norm = nmap**2 / napsq / mapsq
    napsq_mapsq = napsq / mapsq
    corr2_output = np.genfromtxt(os.path.join('output','ng_norm.out'), names=True)
    np.testing.assert_allclose(corr2_output['NMap'], nmap, rtol=1.e-3)
    np.testing.assert_allclose(corr2_output['NMx'], nmx, atol=1.e-3)
    np.testing.assert_allclose(corr2_output['sig_nmap'], np.sqrt(varnmap), rtol=1.e-3)
    np.testing.assert_allclose(corr2_output['Napsq'], napsq, rtol=1.e-3)
    np.testing.assert_allclose(corr2_output['sig_napsq'], np.sqrt(varnap), rtol=1.e-3)
    np.testing.assert_allclose(corr2_output['Mapsq'], mapsq, rtol=1.e-3)
    np.testing.assert_allclose(corr2_output['sig_mapsq'], np.sqrt(varmap), rtol=1.e-3)
    np.testing.assert_allclose(corr2_output['NMap_norm'], nmap_norm, rtol=1.e-3)
    np.testing.assert_allclose(corr2_output['Nsq_Mapsq'], napsq_mapsq, rtol=1.e-3)

ra = np.random.uniform(0, 20, N)
    dec = np.random.uniform(0, 20, N)

    g1 = np.random.uniform(-0.05, 0.05, N)
    g2 = np.random.uniform(-0.05, 0.05, N)
    k = np.random.uniform(-0.05, 0.05, N)

    config = {
        'min_sep': 0.5,
        'max_sep': 100.0,
        'nbins': 10,
        'bin_slop':0.,
        'sep_units':'arcmin',
    }

    gg1 = treecorr.GGCorrelation(config)
    cat1 = treecorr.Catalog(ra=ra, dec=dec, g1=g1, g2=g2, k=k, ra_units='deg', dec_units='deg')
    gg1.process(cat1)

    # Now add some Nan's to the end of the data
    # TreeCorr's message warns about this but says it's ignoring the rows
    ra = np.concatenate((ra, [0.0]))
    dec = np.concatenate((dec, [0.0]))
    g1 = np.concatenate((g1, [np.nan]))
    g2 = np.concatenate((g2, [np.nan]))
    k = np.concatenate((k, [np.nan]))

    gg2 = treecorr.GGCorrelation(config)
    cat2 = treecorr.Catalog(ra=ra, dec=dec, g1=g1, g2=g2, k=k, ra_units='deg', dec_units='deg',
                            keep_zero_weight=True)

    assert cat2.nobj == cat1.nobj  # same number of non-zero weight

g1_col='g1', g2_col='g2', k_col='k',
                               patch_centers=patch_centers)
    gk_cat2 = treecorr.Catalog(file_name,
                               ra_col='ra', dec_col='dec', ra_units='rad', dec_units='rad',
                               g1_col='g1', g2_col='g2', k_col='k',
                               patch_centers=patch_centers, save_patch_dir='output')

    gg1 = treecorr.GGCorrelation(bin_size=0.5, min_sep=1., max_sep=30., sep_units='arcmin')
    t0 = time.time()
    s0 = hp.heap().size if hp else 0
    gg1.process(gk_cat1)
    t1 = time.time()
    s1 = hp.heap().size if hp else 0
    print('GG1: ',s1, t1-t0, s1-s0)
    gk_cat1.unload()
    gg2 = treecorr.GGCorrelation(bin_size=0.5, min_sep=1., max_sep=30., sep_units='arcmin')
    t0 = time.time()
    s0 = hp.heap().size if hp else 0
    gg2.process(gk_cat2, low_mem=True)
    t1 = time.time()
    s1 = hp.heap().size if hp else 0
    print('GG2: ',s1, t1-t0, s1-s0)
    gk_cat2.unload()
    np.testing.assert_allclose(gg1.xip, gg2.xip)
    np.testing.assert_allclose(gg1.xim, gg2.xim)
    np.testing.assert_allclose(gg1.weight, gg2.weight)

    # NG, still with the same file, but cross correlation.
    ng1 = treecorr.NGCorrelation(bin_size=0.5, min_sep=1., max_sep=30., sep_units='arcmin')
    t0 = time.time()
    s0 = hp.heap().size if hp else 0
    ng1.process(gk_cat1, gk_cat1)

gg.write(out_file_name)
    data = fitsio.read(out_file_name)
    np.testing.assert_allclose(data['r_nom'], np.exp(gg.logr))
    np.testing.assert_allclose(data['meanr'], gg.meanr)
    np.testing.assert_allclose(data['meanlogr'], gg.meanlogr)
    np.testing.assert_allclose(data['xip'], gg.xip)
    np.testing.assert_allclose(data['xim'], gg.xim)
    np.testing.assert_allclose(data['xip_im'], gg.xip_im)
    np.testing.assert_allclose(data['xim_im'], gg.xim_im)
    np.testing.assert_allclose(data['sigma_xip'], np.sqrt(gg.varxip))
    np.testing.assert_allclose(data['sigma_xim'], np.sqrt(gg.varxim))
    np.testing.assert_allclose(data['weight'], gg.weight)
    np.testing.assert_allclose(data['npairs'], gg.npairs)

    # Check the read function
    gg2 = treecorr.GGCorrelation(bin_size=0.1, min_sep=1., max_sep=100., sep_units='arcmin')
    gg2.read(out_file_name)
    np.testing.assert_allclose(gg2.logr, gg.logr)
    np.testing.assert_allclose(gg2.meanr, gg.meanr)
    np.testing.assert_allclose(gg2.meanlogr, gg.meanlogr)
    np.testing.assert_allclose(gg2.xip, gg.xip)
    np.testing.assert_allclose(gg2.xim, gg.xim)
    np.testing.assert_allclose(gg2.xip_im, gg.xip_im)
    np.testing.assert_allclose(gg2.xim_im, gg.xim_im)
    np.testing.assert_allclose(gg2.varxip, gg.varxip)
    np.testing.assert_allclose(gg2.varxim, gg.varxim)
    np.testing.assert_allclose(gg2.weight, gg.weight)
    np.testing.assert_allclose(gg2.npairs, gg.npairs)
    assert gg2.coords == gg.coords
    assert gg2.metric == gg.metric
    assert gg2.sep_units == gg.sep_units
    assert gg2.bin_type == gg.bin_type

print('gg.xim = ',gg1.xim)
    print('theory_gammat = ',theory_gQ)
    print('ratio = ',gg1.xim / theory_gQ)
    print('diff = ',gg1.xim - theory_gQ)
    print('max diff = ',max(abs(gg1.xim - theory_gQ)))
    assert max(abs(gg1.xim - theory_gQ)) < 4.e-5

    # Can also do brute force just on one cat or the other.
    # In this case, just cat2 is equivalent to full brute force, since nlens << nsource.
    gg1a = treecorr.GGCorrelation(bin_size=bin_size, min_sep=min_sep, max_sep=max_sep, verbose=1,
                                  metric='Rlens', bin_slop=0, brute=1)
    gg1a.process(lens_cat, source_cat)
    assert lens_cat.field.brute == True
    assert source_cat.field.brute == False
    gg1b = treecorr.GGCorrelation(bin_size=bin_size, min_sep=min_sep, max_sep=max_sep, verbose=1,
                                  metric='Rlens', bin_slop=0, brute=2)
    gg1b.process(lens_cat, source_cat)
    assert lens_cat.field.brute == False
    assert source_cat.field.brute == True

    assert max(abs(gg0.xim - true_gQ)) < 2.e-6
    assert max(abs(gg1.xim - true_gQ)) < 2.e-5
    assert max(abs(gg1a.xim - true_gQ)) < 2.e-5
    assert max(abs(gg1b.xim - true_gQ)) < 2.e-6
    assert max(abs(gg0.xip - true_gCr)) < 2.e-6
    assert max(abs(gg1.xip - true_gCr)) < 2.e-5
    assert max(abs(gg1a.xip - true_gCr)) < 2.e-5
    assert max(abs(gg1b.xip - true_gCr)) < 2.e-6

    # Now use a more normal value for bin_slop.
    gg2 = treecorr.GGCorrelation(bin_size=bin_size, min_sep=min_sep, max_sep=max_sep, verbose=1,

# Test the singleBin function for TwoD bintype.

    # N random points in 2 dimensions
    rng = np.random.RandomState(8675309)
    N = 5000
    x = rng.uniform(-20, 20, N)
    y = rng.uniform(-20, 20, N)
    g1 = rng.uniform(-0.2, 0.2, N)
    g2 = rng.uniform(-0.2, 0.2, N)
    cat = treecorr.Catalog(x=x, y=y, g1=g1, g2=g2, k=10.*np.ones_like(x))

    max_sep = 21.
    nbins = 5  # Use very chunky bins, so more pairs of non-leaf cells can fall in single bin.

    # First use brute force for reference
    gg0 = treecorr.GGCorrelation(max_sep=max_sep, nbins=nbins, bin_type='TwoD', brute=True)
    t0 = time.time()
    gg0.process(cat)
    t1 = time.time()
    print('t for bs=0 = ',t1-t0)

    # Now do bin_slop = 0
    gg1 = treecorr.GGCorrelation(max_sep=max_sep, nbins=nbins, bin_type='TwoD', bin_slop=0)
    t0 = time.time()
    gg1.process(cat)
    t1 = time.time()
    print('t for bs=1.e-10 = ',t1-t0)
    print('max abs diff xip = ',np.max(np.abs(gg1.xip - gg0.xip)))
    print('max abs diff xim = ',np.max(np.abs(gg1.xim - gg0.xim)))
    np.testing.assert_array_equal(gg1.npairs, gg0.npairs)
    np.testing.assert_allclose(gg1.xip, gg0.xip, atol=1.e-10)
    np.testing.assert_allclose(gg1.xim, gg0.xim, atol=1.e-5)

print('theory_gammat = ',theory_gQ)
    print('ratio = ',gg0s.xim / theory_gQ)
    print('diff = ',gg0s.xim - theory_gQ)
    print('max diff = ',max(abs(gg0s.xim - theory_gQ)))
    assert max(abs(gg0s.xim - theory_gQ)) < 4.e-5

    # This should be identical to the 3d version, since going all the way to leaves.
    # (The next test will be different, since tree creation is different.)
    assert max(abs(gg0s.xim - gg0.xim)) < 1.e-7
    assert max(abs(gg0s.xip - gg0.xip)) < 1.e-7
    assert max(abs(gg0s.xim_im - gg0.xim_im)) < 1.e-7
    assert max(abs(gg0s.xip_im - gg0.xip_im)) < 1.e-7
    assert max(abs(gg0s.npairs - gg0.npairs)) < 1.e-7

    # Now use a more normal value for bin_slop.
    ggs2 = treecorr.GGCorrelation(bin_size=bin_size, min_sep=min_sep, max_sep=max_sep, verbose=1,
                                  metric='Rlens', bin_slop=0.3)
    ggs2.process(lens_cat, source_cat)
    Rlens = ggs2.meanr
    theory_gQ = gamma0 * np.exp(-0.5*Rlens**2/R0**2)

    print('Results with bin_slop = 0.3')
    print('gg.npairs = ',ggs2.npairs)
    print('gg.xim = ',ggs2.xim)
    print('theory_gammat = ',theory_gQ)
    print('ratio = ',ggs2.xim / theory_gQ)
    print('diff = ',ggs2.xim - theory_gQ)
    print('max diff = ',max(abs(ggs2.xim - theory_gQ)))
    # Not quite as accurate as above, since the cells that get used tend to be larger, so more
    # slop happens in the binning.
    assert max(abs(ggs2.xim - theory_gQ)) < 5.e-5
    print('gg.xim_im = ',ggs2.xim_im)

print('ratio = ',corr2_output['xim']/gg2.xim)
    print('diff = ',corr2_output['xim']-gg2.xim)
    np.testing.assert_allclose(corr2_output['xim'], gg2.xim, rtol=1.e-3)
    np.testing.assert_allclose(corr2_output['xim_im'], gg2.xim_im, rtol=1.e-3)
    np.testing.assert_allclose(corr2_output['xip'], gg2.xip, rtol=1.e-3)
    np.testing.assert_allclose(corr2_output['xip_im'], gg2.xip_im, rtol=1.e-3)

    # Repeat with the sources being given as RA/Dec only.
    ral, decl = coord.CelestialCoord.xyz_to_radec(xl,yl,zl)
    ras, decs = coord.CelestialCoord.xyz_to_radec(xs,ys,zs)
    lens_cat = treecorr.Catalog(ra=ral, dec=decl, ra_units='radians', dec_units='radians', r=rl,
                                g1=gl.real, g2=gl.imag)
    source_cat = treecorr.Catalog(ra=ras, dec=decs, ra_units='radians', dec_units='radians',
                                  g1=g1, g2=g2)

    gg0s = treecorr.GGCorrelation(bin_size=bin_size, min_sep=min_sep, max_sep=max_sep, verbose=1,
                                  metric='Rlens', brute=True)
    gg0s.process(lens_cat, source_cat)

    Rlens = gg0s.meanr
    theory_gQ = gamma0 * np.exp(-0.5*Rlens**2/R0**2)

    print('Results with brute=True')
    print('gg.npairs = ',gg0s.npairs)
    print('true_npairs = ',true_npairs)
    np.testing.assert_array_equal(gg0s.npairs, true_npairs)
    print('gg.xim = ',gg0s.xim)
    print('true_gQ = ',true_gQ)
    print('ratio = ',gg0s.xim / true_gQ)
    print('diff = ',gg0s.xim - true_gQ)
    print('max diff = ',max(abs(gg0s.xim - true_gQ)))
    assert max(abs(gg0s.xim - true_gQ)) < 2.e-6

assert len(i1) == 100
    assert len(i2) == 100
    assert len(sep) == 100
    actual_sep = ((x1[i1]-x2[i2])**2 + (y1[i1]-y2[i2])**2)**0.5
    np.testing.assert_allclose(sep, actual_sep, rtol=0.1)
    np.testing.assert_array_less(sep, nn.right_edges[b])
    np.testing.assert_array_less(nn.left_edges[b], sep)

    # To exercise case 3, we need to go to larger separations, so the recursion
    # more often stops before getting to the leaves.
    # Also switch to 3d coordinates.

    cat1 = treecorr.Catalog(x=x1, y=y1, z=z1, w=w1, g1=g11, g2=g21, k=k1, keep_zero_weight=True)
    cat2 = treecorr.Catalog(x=x2, y=y2, z=z2, w=w2, g1=g12, g2=g22, k=k2, keep_zero_weight=True)

    gg = treecorr.GGCorrelation(min_sep=0.4, nbins=10, bin_size=0.1, max_top=0)
    gg.process(cat1, cat2)
    print('rnom = ',gg.rnom)
    print('npairs = ',gg.npairs.astype(int))
    for b in [0,5]:
        i1, i2, sep = gg.sample_pairs(100, cat1, cat2,
                                      min_sep=gg.left_edges[b], max_sep=gg.right_edges[b])

        print('len(npairs) = ',len(gg.npairs))
        print('npairs = ',gg.npairs)
        print('i1 = ',i1)
        print('i2 = ',i2)
        print('sep = ',sep)
        assert len(i1) == 100
        assert len(i2) == 100
        assert len(sep) == 100
        actual_sep = ((x1[i1]-x2[i2])**2 + (y1[i1]-y2[i2])**2 + (z1[i1]-z2[i2])**2)**0.5

mean_varxim = data['mean_varxim']

    print('mean_xip = ',mean_xip)
    print('mean_xim = ',mean_xim)
    print('mean_varxip = ',mean_varxip)
    print('mean_varxim = ',mean_varxim)
    print('var_xip = ',var_xip)
    print('ratio = ',var_xip / mean_varxip)
    print('var_xim = ',var_xim)
    print('ratio = ',var_xim / mean_varxim)

    np.random.seed(1234)
    # First run with the normal variance estimate, which is too small.
    x, y, g1, g2 = generate_shear_field(nside)
    cat = treecorr.Catalog(x=x, y=y, g1=g1, g2=g2)
    gg1 = treecorr.GGCorrelation(bin_size=0.3, min_sep=10., max_sep=50.)
    t0 = time.time()
    gg1.process(cat)
    t1 = time.time()
    print('Time for non-patch processing = ',t1-t0)

    print('npairs = ',gg1.npairs)
    print('xip = ',gg1.xip)
    print('xim = ',gg1.xim)
    print('varxip = ',gg1.varxip)
    print('varxim = ',gg1.varxim)
    print('pullsq for xip = ',(gg1.xip-mean_xip)**2/var_xip)
    print('pullsq for xim = ',(gg1.xim-mean_xim)**2/var_xim)
    print('max pull for xip = ',np.sqrt(np.max((gg1.xip-mean_xip)**2/var_xip)))
    print('max pull for xim = ',np.sqrt(np.max((gg1.xim-mean_xim)**2/var_xim)))
    np.testing.assert_array_less((gg1.xip - mean_xip)**2/var_xip, 25) # within 5 sigma
    np.testing.assert_array_less((gg1.xim - mean_xim)**2/var_xim, 25)