How to use the treecorr.NNNCorrelation function in TreeCorr

rng = np.random.RandomState(8675309)
    x = rng.normal(0,s, (ngal,) )
    y = rng.normal(0,s, (ngal,) )
    cat = treecorr.Catalog(x=x, y=y)

    min_sep = 1.
    max_sep = 50.
    nbins = 50
    min_u = 0.13
    max_u = 0.89
    nubins = 10
    min_v = 0.13
    max_v = 0.59
    nvbins = 10

    ddd = treecorr.NNNCorrelation(min_sep=min_sep, max_sep=max_sep, nbins=nbins,
                                  min_u=min_u, max_u=max_u, nubins=nubins,
                                  min_v=min_v, max_v=max_v, nvbins=nvbins,
                                  brute=True, verbose=1)
    ddd.process(cat)

    log_min_sep = np.log(min_sep)
    log_max_sep = np.log(max_sep)
    true_ntri = np.zeros( (nbins, nubins, 2*nvbins) )
    bin_size = (log_max_sep - log_min_sep) / nbins
    ubin_size = (max_u-min_u) / nubins
    vbin_size = (max_v-min_v) / nvbins
    for i in range(ngal):
        for j in range(i+1,ngal):
            for k in range(j+1,ngal):
                dij = np.sqrt((x[i]-x[j])**2 + (y[i]-y[j])**2)
                dik = np.sqrt((x[i]-x[k])**2 + (y[i]-y[k])**2)

print('max rel diff = ',np.max(np.abs(((ddd.meand1-ddd.meand2)/ddd.meand3-np.abs(ddd.meanv))/ddd.meanv)))
    np.testing.assert_allclose(ddd.meanlogd1, np.log(ddd.meand1), rtol=1.e-3)
    np.testing.assert_allclose(ddd.meanlogd2, np.log(ddd.meand2), rtol=1.e-3)
    np.testing.assert_allclose(ddd.meanlogd3, np.log(ddd.meand3), rtol=1.e-3)
    np.testing.assert_allclose(ddd.meand3/ddd.meand2, ddd.meanu, rtol=1.e-5 * tol_factor)
    np.testing.assert_allclose((ddd.meand1-ddd.meand2)/ddd.meand3, np.abs(ddd.meanv),
                                  rtol=1.e-5 * tol_factor, atol=1.e-5 * tol_factor)
    np.testing.assert_allclose(ddd.meanlogd3-ddd.meanlogd2, np.log(ddd.meanu),
                                  atol=1.e-3 * tol_factor)
    np.testing.assert_allclose(np.log(ddd.meand1-ddd.meand2)-ddd.meanlogd3,
                                  np.log(np.abs(ddd.meanv)), atol=2.e-3 * tol_factor)

    rx = (rng.random_sample(nrand)-0.5) * L
    ry = (rng.random_sample(nrand)-0.5) * L
    rand = treecorr.Catalog(x=rx,y=ry, x_units='arcmin', y_units='arcmin')
    rrr = treecorr.NNNCorrelation(min_sep=min_sep, max_sep=max_sep, nbins=nbins,
                                  min_u=min_u, max_u=max_u, min_v=min_v, max_v=max_v,
                                  nubins=nubins, nvbins=nvbins,
                                  sep_units='arcmin', verbose=1)
    rrr.process(rand)
    #print('rrr.ntri = ',rrr.ntri)

    d1 = ddd.meand1
    d2 = ddd.meand2
    d3 = ddd.meand3
    #print('rnom = ',np.exp(ddd.logr))
    #print('unom = ',ddd.u)
    #print('vnom = ',ddd.v)
    #print('d1 = ',d1)
    #print('d2 = ',d2)
    #print('d3 = ',d3)
    true_zeta = (1./(12.*np.pi**2)) * (L/s)**4 * np.exp(-(d1**2+d2**2+d3**2)/(6.*s**2)) - 1.

nubins=nubins, nvbins=nvbins, bin_slop=0.1, verbose=1)
    ddd.process(data_cats)
    print('From multiple catalogs: ddd.ntri = ',ddd.ntri)

    # Now do the same thing with one big catalog
    dddx = treecorr.NNNCorrelation(min_sep=min_sep, max_sep=max_sep, nbins=nbins,
                                   min_u=min_u, max_u=max_u, min_v=min_v, max_v=max_v,
                                   nubins=nubins, nvbins=nvbins, bin_slop=0.1, verbose=1)
    data_catx = treecorr.Catalog(x=x.reshape( (ngal*ncats,) ), y=y.reshape( (ngal*ncats,) ))
    dddx.process(data_catx)
    print('From single catalog: dddx.ntri = ',dddx.ntri)
    # Only test to rtol=0.1, since there are now differences between the auto and cross related
    # to how they characterize triangles especially when d1 ~= d2 or d2 ~= d3.
    np.testing.assert_allclose(ddd.ntri, dddx.ntri, rtol=0.1)

    rrr = treecorr.NNNCorrelation(min_sep=min_sep, max_sep=max_sep, nbins=nbins,
                                  min_u=min_u, max_u=max_u, min_v=min_v, max_v=max_v,
                                  nubins=nubins, nvbins=nvbins, bin_slop=0.1, verbose=1)
    rrr.process(rand_cats)
    print('rrr.ntri = ',rrr.ntri)

    rrrx = treecorr.NNNCorrelation(min_sep=min_sep, max_sep=max_sep, nbins=nbins,
                                   min_u=min_u, max_u=max_u, min_v=min_v, max_v=max_v,
                                   nubins=nubins, nvbins=nvbins, bin_slop=0.1, verbose=1)
    rand_catx = treecorr.Catalog(x=rx.reshape( (nrand*ncats,) ), y=ry.reshape( (nrand*ncats,) ))
    rrrx.process(rand_catx)
    print('rrrx.ntri = ',rrrx.ntri)
    np.testing.assert_allclose(ddd.ntri, dddx.ntri, rtol=0.1)

    zeta, varzeta = ddd.calculateZeta(rrr)
    zetax, varzetax = dddx.calculateZeta(rrrx)
    print('zeta = ',zeta)

ddd6 = treecorr.NNNCorrelation(min_sep=min_sep, max_sep=max_sep, nbins=nbins*2,
                                   min_u=min_u, max_u=max_u, nubins=nubins,
                                   min_v=min_v, max_v=max_v, nvbins=nvbins)
    with assert_raises(ValueError):
        ddd2 += ddd6
    ddd7 = treecorr.NNNCorrelation(min_sep=min_sep, max_sep=max_sep, nbins=nbins,
                                   min_u=min_u-0.1, max_u=max_u, nubins=nubins,
                                   min_v=min_v, max_v=max_v, nvbins=nvbins)
    with assert_raises(ValueError):
        ddd2 += ddd7
    ddd8 = treecorr.NNNCorrelation(min_sep=min_sep, max_sep=max_sep, nbins=nbins,
                                   min_u=min_u, max_u=max_u+0.1, nubins=nubins,
                                   min_v=min_v, max_v=max_v, nvbins=nvbins)
    with assert_raises(ValueError):
        ddd2 += ddd8
    ddd9 = treecorr.NNNCorrelation(min_sep=min_sep, max_sep=max_sep, nbins=nbins,
                                   min_u=min_u, max_u=max_u, nubins=nubins*2,
                                   min_v=min_v, max_v=max_v, nvbins=nvbins)
    with assert_raises(ValueError):
        ddd2 += ddd9
    ddd10 = treecorr.NNNCorrelation(min_sep=min_sep, max_sep=max_sep, nbins=nbins,
                                    min_u=min_u, max_u=max_u, nubins=nubins,
                                    min_v=min_v-0.1, max_v=max_v, nvbins=nvbins)
    with assert_raises(ValueError):
        ddd2 += ddd10
    ddd11 = treecorr.NNNCorrelation(min_sep=min_sep, max_sep=max_sep, nbins=nbins,
                                    min_u=min_u, max_u=max_u, nubins=nubins,
                                    min_v=min_v, max_v=max_v+0.1, nvbins=nvbins)
    with assert_raises(ValueError):
        ddd2 += ddd11
    ddd12 = treecorr.NNNCorrelation(min_sep=min_sep, max_sep=max_sep, nbins=nbins,
                                    min_u=min_u, max_u=max_u, nubins=nubins,

np.testing.assert_almost_equal(nnn.v1d[nnn.nvbins-1], -nnn.min_v - 0.5*nnn.vbin_size)
        np.testing.assert_equal(nnn.v.shape, (nnn.nbins, nnn.nubins, 2*nnn.nvbins) )
        np.testing.assert_almost_equal(nnn.v[0,0,:], nnn.v1d)
        np.testing.assert_almost_equal(nnn.v[-1,-1,:], nnn.v1d)

    def check_defaultuv(nnn):
        assert nnn.min_u == 0.
        assert nnn.max_u == 1.
        assert nnn.nubins == np.ceil(1./nnn.ubin_size)
        assert nnn.min_v == 0.
        assert nnn.max_v == 1.
        assert nnn.nvbins == np.ceil(1./nnn.vbin_size)

    # Check the different ways to set up the binning:
    # Omit bin_size
    nnn = treecorr.NNNCorrelation(min_sep=5, max_sep=20, nbins=20, bin_type='LogRUV')
    #print(nnn.min_sep,nnn.max_sep,nnn.bin_size,nnn.nbins)
    #print(nnn.min_u,nnn.max_u,nnn.ubin_size,nnn.nubins)
    #print(nnn.min_v,nnn.max_v,nnn.vbin_size,nnn.nvbins)
    assert nnn.min_sep == 5.
    assert nnn.max_sep == 20.
    assert nnn.nbins == 20
    check_defaultuv(nnn)
    check_arrays(nnn)

    # Specify min, max, n for u,v too.
    nnn = treecorr.NNNCorrelation(min_sep=5, max_sep=20, nbins=20,
                                  min_u=0.2, max_u=0.9, nubins=12,
                                  min_v=0., max_v=0.2, nvbins=2)
    #print(nnn.min_sep,nnn.max_sep,nnn.bin_size,nnn.nbins)
    #print(nnn.min_u,nnn.max_u,nnn.ubin_size,nnn.nubins)
    #print(nnn.min_v,nnn.max_v,nnn.vbin_size,nnn.nvbins)

rng = np.random.RandomState(8675309)
    x = (rng.random_sample(ngal)-0.5) * Lx
    y = (rng.random_sample(ngal)-0.5) * Ly
    cat = treecorr.Catalog(x=x, y=y)

    min_sep = 1.
    max_sep = 40.  # This only really makes sense if max_sep < L/4 for all L.
    nbins = 50
    min_u = 0.13
    max_u = 0.89
    nubins = 10
    min_v = 0.13
    max_v = 0.59
    nvbins = 10

    ddd = treecorr.NNNCorrelation(min_sep=min_sep, max_sep=max_sep, nbins=nbins,
                                  min_u=min_u, max_u=max_u, nubins=nubins,
                                  min_v=min_v, max_v=max_v, nvbins=nvbins,
                                  bin_slop=0, xperiod=Lx, yperiod=Ly, brute=True)
    ddd.process(cat, metric='Periodic', num_threads=1)
    #print('ddd.ntri = ',ddd.ntri)

    log_min_sep = np.log(min_sep)
    log_max_sep = np.log(max_sep)
    true_ntri = np.zeros( (nbins, nubins, 2*nvbins) )
    bin_size = (log_max_sep - log_min_sep) / nbins
    ubin_size = (max_u-min_u) / nubins
    vbin_size = (max_v-min_v) / nvbins
    for i in range(ngal):
        for j in range(i+1,ngal):
            for k in range(j+1,ngal):
                xi = x[i]

ku = int(np.floor( (u-min_u) / ubin_size ))
                if v > 0:
                    kv = int(np.floor( (v-min_v) / vbin_size )) + nvbins
                else:
                    kv = int(np.floor( (v-(-max_v)) / vbin_size ))
                assert 0 <= kr < nbins
                assert 0 <= ku < nubins
                assert 0 <= kv < 2*nvbins
                true_ntri[kr,ku,kv] += 1

    #print('true_ntri = ',true_ntri)
    #print('diff = ',ddd.ntri - true_ntri)
    np.testing.assert_array_equal(ddd.ntri, true_ntri)

    # Repeat with binslop = 0
    ddd = treecorr.NNNCorrelation(min_sep=min_sep, max_sep=max_sep, nbins=nbins,
                                  min_u=min_u, max_u=max_u, nubins=nubins,
                                  min_v=min_v, max_v=max_v, nvbins=nvbins,
                                  bin_slop=0, verbose=1)
    ddd.process(cat1, cat2, cat3)
    #print('binslop > 0: ddd.ntri = ',ddd.ntri)
    #print('diff = ',ddd.ntri - true_ntri)
    np.testing.assert_array_equal(ddd.ntri, true_ntri)

    # And again with no top-level recursion
    ddd = treecorr.NNNCorrelation(min_sep=min_sep, max_sep=max_sep, nbins=nbins,
                                  min_u=min_u, max_u=max_u, nubins=nubins,
                                  min_v=min_v, max_v=max_v, nvbins=nvbins,
                                  bin_slop=0, verbose=1, max_top=0)
    ddd.process(cat1, cat2, cat3)
    #print('max_top = 0: ddd.ntri = ',ddd.ntri)
    #print('true_ntri = ',true_ntri)

assert nnn.nvbins == np.ceil(1./nnn.vbin_size)

    # Check the different ways to set up the binning:
    # Omit bin_size
    nnn = treecorr.NNNCorrelation(min_sep=5, max_sep=20, nbins=20, bin_type='LogRUV')
    #print(nnn.min_sep,nnn.max_sep,nnn.bin_size,nnn.nbins)
    #print(nnn.min_u,nnn.max_u,nnn.ubin_size,nnn.nubins)
    #print(nnn.min_v,nnn.max_v,nnn.vbin_size,nnn.nvbins)
    assert nnn.min_sep == 5.
    assert nnn.max_sep == 20.
    assert nnn.nbins == 20
    check_defaultuv(nnn)
    check_arrays(nnn)

    # Specify min, max, n for u,v too.
    nnn = treecorr.NNNCorrelation(min_sep=5, max_sep=20, nbins=20,
                                  min_u=0.2, max_u=0.9, nubins=12,
                                  min_v=0., max_v=0.2, nvbins=2)
    #print(nnn.min_sep,nnn.max_sep,nnn.bin_size,nnn.nbins)
    #print(nnn.min_u,nnn.max_u,nnn.ubin_size,nnn.nubins)
    #print(nnn.min_v,nnn.max_v,nnn.vbin_size,nnn.nvbins)
    assert nnn.min_sep == 5.
    assert nnn.max_sep == 20.
    assert nnn.nbins == 20
    assert nnn.min_u == 0.2
    assert nnn.max_u == 0.9
    assert nnn.nubins == 12
    assert nnn.min_v == 0.
    assert nnn.max_v == 0.2
    assert nnn.nvbins == 2
    check_arrays(nnn)

#print('true_ntri => ',true_ntri)
    #print('diff = ',ddd.ntri - true_ntri)
    np.testing.assert_array_equal(ddd.ntri, true_ntri)

    # Repeat with binslop = 0
    ddd = treecorr.NNNCorrelation(min_sep=min_sep, max_sep=max_sep, nbins=nbins,
                                  min_u=min_u, max_u=max_u, nubins=nubins,
                                  min_v=min_v, max_v=max_v, nvbins=nvbins,
                                  bin_slop=0, verbose=1)
    ddd.process(cat)
    #print('ddd.ntri = ',ddd.ntri)
    #print('diff = ',ddd.ntri - true_ntri)
    np.testing.assert_array_equal(ddd.ntri, true_ntri)

    # And again with no top-level recursion
    ddd = treecorr.NNNCorrelation(min_sep=min_sep, max_sep=max_sep, nbins=nbins,
                                  min_u=min_u, max_u=max_u, nubins=nubins,
                                  min_v=min_v, max_v=max_v, nvbins=nvbins,
                                  bin_slop=0, verbose=1, max_top=0)
    ddd.process(cat)
    #print('ddd.ntri = ',ddd.ntri)
    #print('true_ntri => ',true_ntri)
    #print('diff = ',ddd.ntri - true_ntri)
    np.testing.assert_array_equal(ddd.ntri, true_ntri)

    # And compare to the cross correlation
    ddd.clear()
    ddd.process(cat,cat,cat)
    #print('ddd.ntri = ',ddd.ntri)
    #print('true_ntri => ',true_ntri)
    #print('diff = ',ddd.ntri - true_ntri)
    np.testing.assert_array_equal(ddd.ntri, true_ntri)

ddr.process(cat1,cat2,rand3)
            logger.info("Done DDR calculations.")
            logger.warning("Performing RDR calculations...")
            rdr = treecorr.NNNCorrelation(config,logger)
            rdr.process(rand1,cat2,rand3)
            logger.info("Done RDR calculations.")
            logger.warning("Performing RRD calculations...")
            rrd = treecorr.NNNCorrelation(config,logger)
            rrd.process(rand1,rand2,cat3)
            logger.info("Done RRD calculations.")
            logger.warning("Performing DRD calculations...")
            drd = treecorr.NNNCorrelation(config,logger)
            drd.process(cat1,rand2,cat3)
            logger.info("Done DRD calculations.")
            logger.warning("Performing RDD calculations...")
            rdd = treecorr.NNNCorrelation(config,logger)
            rdd.process(rand1,cat2,cat3)
            logger.info("Done RDD calculations.")
        ddd.write(config['nnn_file_name'],rrr,drr,rdr,rrd,ddr,drd,rdd)
        logger.warning("Wrote NNN correlation to %s",config['nnn_file_name'])

    # Do KKK correlation function if necessary
    if 'kkk_file_name' in config:
        logger.warning("Performing KKK calculations...")
        kkk = treecorr.KKKCorrelation(config,logger)
        kkk.process(cat1,cat2,cat3)
        logger.info("Done KKK calculations.")
        kkk.write(config['kkk_file_name'])
        logger.warning("Wrote KKK correlation to %s",config['kkk_file_name'])