Repository: rmjarvis/TreeCorr (language: python)

Issue:
multiply-occurring objects bias results low
I got unexpected shear-shear correlation results with a catalog including galaxies that can appear more than once. To test it I just compared results from a "normal" catalog of galaxies with one where I put in every galaxy twice and I got the following example result:

![check_treecorr](https://user-images.githubusercontent.com/44603577/134537605-97bb2291-11c2-4053-a03a-7d6129e34ee0.png)

The signal of the doubled catalog is consistently suppressed by ~30% for low theta while for large theta they roughly agree. Unfortunately, I could not really pinpoint where the calculation goes wrong. It is, e.g., not just the number counts per bin that are wrong.

My current work-around is to shift the position of every single duplicated galaxy by <0.5 arc seconds in a random way. Then the original and duplicated catalog agree perfectly.

I thought, however, you might find an easy fix for this or at least mention it in the documentation (maybe I just missed it though ;)

Cheers, Sandra Unruh


Discussion / hints:
Hi Sandra.  I couldn't reproduce this with a straightforward test that duplicated all the galaxies.  I'm getting answers that are the same to machine precision, so there must be something in your script that I'm not including in my test.  

Can you please send me a self-contained script that shows the effect?  Either post here, or send me by email or Slack.
Dear Mike,

I put together a minimal script that reproduces my error. It can be found here 

https://www.dropbox.com/sh/8ko0en0xz7edp89/AABJKY_IpdhcXsTzyRwVF4hQa?dl=0 <https://www.dropbox.com/sh/8ko0en0xz7edp89/AABJKY_IpdhcXsTzyRwVF4hQa?dl=0>

with a 26MB catalog of input data (which shouldn’t matter but still). Please change the data directory (data_dir) to yours, the it should run through.

I use these versions of Python and TreeCorr:
TreeCorr.__version__ = ‚4.2.3'
Python 3.8.10 | packaged by conda-forge

 I really hope it is no stupid mistake from my side! Thank you for your time, I enjoy working with TreeCorr a lot.
Cheers, Sandra

> Am 23.09.2021 um 18:31 schrieb Mike Jarvis ***@***.***>:
> 
> 
> Hi Sandra. I couldn't reproduce this with a straightforward test that duplicated all the galaxies. I'm getting answers that are the same to machine precision, so there must be something in your script that I'm not including in my test.
> 
> Can you please send me a self-contained script that shows the effect? Either post here, or send me by email or Slack.
> 
> —
> You are receiving this because you authored the thread.
> Reply to this email directly, view it on GitHub <https://github.com/rmjarvis/TreeCorr/issues/134#issuecomment-925972766>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AKUJROO3TJPSDHRKNGI7O43UDNI7TANCNFSM5EUAQKFA>.
> Triage notifications on the go with GitHub Mobile for iOS <https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675> or Android <https://play.google.com/store/apps/details?id=com.github.android&referrer=utm_campaign%3Dnotification-email%26utm_medium%3Demail%26utm_source%3Dgithub>. 
> 



Fix patch:
diff --git a/src/Cell.cpp b/src/Cell.cpp
index 4246ea4e..cbc089df 100644
--- a/src/Cell.cpp
+++ b/src/Cell.cpp
@@ -108,7 +108,7 @@ double Cell<D,C>::calculateInertia() const
             double w = leaves[k]->getW();
             inertia2 += w * (p-cen).normSq();
         }
-        if (abs(inertia2 - inertia) > 1.e-8 * inertia) {
+        if (std::abs(inertia2 - inertia) > 1.e-8 * inertia) {
             dbg<<"Cell = "<<*this<<std::endl;
             dbg<<"cen = "<<cen<<std::endl;
             dbg<<"w = "<<getW()<<std::endl;
@@ -223,7 +223,7 @@ std::complex<double> ParallelTransportShift(
         double normBsq = sinB*sinB + cosB*cosB;
         xxdbg<<"A = atan("<<sinA<<"/"<<cosA<<") = "<<atan2(sinA,cosA)*180./M_PI<<std::endl;
         xxdbg<<"B = atan("<<sinB<<"/"<<cosB<<") = "<<atan2(sinB,cosB)*180./M_PI<<std::endl;
-        if (normAsq == 0. || normBsq == 0.) {
+        if (normAsq < 1.e-12 && normBsq < 1.e-12) {
             // Then this point is at the center, no need to project.
             dwg += vdata[i].first->getWG();
         } else {


Test patch:
diff --git a/tests/Tutorial.ipynb b/tests/Tutorial.ipynb
index 8eed961d..3d480f5b 100644
--- a/tests/Tutorial.ipynb
+++ b/tests/Tutorial.ipynb
@@ -254,14 +254,14 @@
      "text": [
       "## {'coords': 'spherical', 'metric': 'Euclidean', 'sep_units': 'arcmin', 'bin_type': 'Log'}\n",
       "#   r_nom       meanr       meanlogr       xip          xim         xip_im       xim_im     sigma_xip    sigma_xim      weight       npairs   \n",
-      "  1.0304e+00   1.0311e+00   3.0480e-02   7.5773e-06   1.6730e-07   7.0287e-08  -3.5157e-08   1.8992e-08   1.8992e-08   2.0181e+05   2.0181e+05\n",
+      "  1.0304e+00   1.0311e+00   3.0480e-02   7.5773e-06   1.6730e-07   7.0288e-08  -3.5156e-08   1.8992e-08   1.8992e-08   2.0181e+05   2.0181e+05\n",
       "  1.0940e+00   1.0948e+00   9.0389e-02   7.3854e-06   3.5900e-07  -2.0456e-08   1.5109e-08   1.7948e-08   1.7948e-08   2.2597e+05   2.2597e+05\n",
-      "  1.1616e+00   1.1623e+00   1.5029e-01   7.1456e-06   2.6500e-07  -1.3577e-08  -5.0994e-08   1.6930e-08   1.6930e-08   2.5396e+05   2.5396e+05\n",
-      "  1.2333e+00   1.2342e+00   2.1029e-01   7.0705e-06   2.1738e-07  -4.6802e-08  -8.8672e-08   1.6049e-08   1.6049e-08   2.8261e+05   2.8261e+05\n",
-      "  1.3095e+00   1.3104e+00   2.7017e-01   6.8516e-06   3.4167e-07   6.9888e-08  -8.9430e-08   1.5142e-08   1.5142e-08   3.1746e+05   3.1746e+05\n",
-      "  1.3903e+00   1.3914e+00   3.3014e-01   6.6381e-06   2.9299e-07  -3.1255e-08  -1.4265e-08   1.4349e-08   1.4349e-08   3.5354e+05   3.5354e+05\n",
-      "  1.4762e+00   1.4772e+00   3.8999e-01   6.6007e-06   2.7230e-07  -4.1687e-08  -2.6884e-08   1.3533e-08   1.3533e-08   3.9745e+05   3.9745e+05\n",
-      "  1.5673e+00   1.5685e+00   4.4994e-01   6.4226e-06   4.2499e-07   3.8318e-08  -1.7131e-09   1.2785e-08   1.2785e-08   4.4530e+05   4.4530e+05\n",
+      "  1.1616e+00   1.1623e+00   1.5029e-01   7.1456e-06   2.6500e-07  -1.3585e-08  -5.0992e-08   1.6930e-08   1.6930e-08   2.5396e+05   2.5396e+05\n",
+      "  1.2333e+00   1.2342e+00   2.1029e-01   7.0705e-06   2.1739e-07  -4.6803e-08  -8.8672e-08   1.6049e-08   1.6049e-08   2.8261e+05   2.8261e+05\n",
+      "  1.3095e+00   1.3104e+00   2.7017e-01   6.8516e-06   3.4166e-07   6.9887e-08  -8.9427e-08   1.5142e-08   1.5142e-08   3.1746e+05   3.1746e+05\n",
+      "  1.3903e+00   1.3914e+00   3.3014e-01   6.6381e-06   2.9299e-07  -3.1259e-08  -1.4264e-08   1.4349e-08   1.4349e-08   3.5354e+05   3.5354e+05\n",
+      "  1.4762e+00   1.4772e+00   3.8999e-01   6.6007e-06   2.7230e-07  -4.1688e-08  -2.6885e-08   1.3533e-08   1.3533e-08   3.9745e+05   3.9745e+05\n",
+      "  1.5673e+00   1.5685e+00   4.4994e-01   6.4226e-06   4.2499e-07   3.8314e-08  -1.7152e-09   1.2785e-08   1.2785e-08   4.4530e+05   4.4530e+05\n",
       "\n"
      ]
     }
diff --git a/tests/test_gg.py b/tests/test_gg.py
index abbb5e19..c063efa3 100644
--- a/tests/test_gg.py
+++ b/tests/test_gg.py
@@ -1353,6 +1353,65 @@ def test_varxi():
     np.testing.assert_allclose(gg.varxip, var_xip, rtol=0.2)
     np.testing.assert_allclose(gg.varxim, var_xim, rtol=0.2)
 
+@timer
+def test_double():
+    # Test in response to issue #134.  Check that duplicating all the lenses gives the same answer.
+
+    gamma0 = 0.05
+    r0 = 20.
+    ngal = 10000
+    rng = np.random.RandomState(8675309)
+    x = rng.normal(0,r0, (ngal,) )
+    y = rng.normal(0,r0, (ngal,) )
+    z = rng.normal(0,r0, (ngal,) ) + 200
+    r2 = (x**2 + y**2)/r0**2
+    print('min/max r2 = ',np.min(r2),np.max(r2))
+    g1 = -gamma0 * np.exp(-r2/2.) * (x**2-y**2)/r0**2
+    g2 = -gamma0 * np.exp(-r2/2.) * (2.*x*y)/r0**2
+    print('max g = ',np.max(g1**2+g2**2)**0.5)
+    w = rng.random_sample(ngal) * 10
+    #w = np.ones_like(x)
+    ra, dec = coord.CelestialCoord.xyz_to_radec(x,y,z)
+
+    print('Normal case (all lenses different)')
+    cat = treecorr.Catalog(ra=ra, dec=dec, g1=g1, g2=g2, w=w, ra_units='rad', dec_units='rad')
+    gg = treecorr.GGCorrelation(bin_size=0.1, min_sep=80., max_sep=100., sep_units='arcmin',
+                                verbose=1)
+    gg.process(cat)
+
+    print('xip = ',gg.xip)
+    print('npairs = ',gg.npairs)
+    print('weight = ',gg.weight)
+    #print('xip_im = ',gg.xip_im)
+    #print('xim = ',gg.xim)
+    #print('xim_im = ',gg.xim_im)
+
+    print('Duplicated galaxies')
+    ra_2 = np.concatenate([ra,ra])
+    dec_2 = np.concatenate([dec,dec])
+    g1_2 = np.concatenate([g1,g1])
+    g2_2 = np.concatenate([g2,g2])
+    w_2 = np.concatenate([w,w])
+    cat_2 = treecorr.Catalog(ra=ra_2, dec=dec_2, g1=g1_2, g2=g2_2, w=w_2, ra_units='rad', dec_units='rad')
+    gg_2 = treecorr.GGCorrelation(bin_size=0.1, min_sep=80., max_sep=100., sep_units='arcmin',
+                                  verbose=1)
+    gg_2.process(cat_2)
+
+    print('xip = ',gg_2.xip)
+    print('npairs = ',gg_2.npairs)
+    print('weight/4 = ',gg_2.weight/4)
+    #print('xip_im = ',gg_2.xip_im)
+    #print('xim = ',gg_2.xim)
+    #print('xim_im = ',gg_2.xim_im)
+    print('max diff xip = ',np.max(np.abs(gg_2.xip-gg.xip)))
+    print('max diff xip_im = ',np.max(np.abs(gg_2.xip_im-gg.xip_im)))
+    print('max diff xim = ',np.max(np.abs(gg_2.xim-gg.xim)))
+    print('max diff xim_im = ',np.max(np.abs(gg_2.xim_im-gg.xim_im)))
+    np.testing.assert_allclose(gg_2.xip, gg.xip, rtol=1.e-7)
+    np.testing.assert_allclose(gg_2.xip_im, gg.xip_im, rtol=1.e-7)
+    np.testing.assert_allclose(gg_2.xim, gg.xim, rtol=1.e-7)
+    np.testing.assert_allclose(gg_2.xim_im, gg.xim_im, rtol=1.e-7)
+
 
 if __name__ == '__main__':
     test_direct()
@@ -1364,4 +1423,5 @@ def test_varxi():
     test_aardvark()
     test_shuffle()
     test_haloellip()
-    test_varxi
+    test_varxi()
+    test_double()
diff --git a/tests/test_ng.py b/tests/test_ng.py
index 9e8158b6..5e98962b 100644
--- a/tests/test_ng.py
+++ b/tests/test_ng.py
@@ -1378,6 +1378,62 @@ def test_varxi():
     np.testing.assert_allclose(varxi, var_xi_2, rtol=0.5)
 
 
+@timer
+def test_double():
+    # Test in response to issue #134.  That was about GG correlations, but the same bug could
+    # sho up here, so check that duplicating all the lenses gives the same answer.
+    # Use same signal as in test_ng.
+
+    nlens = 1000
+    nsource = 100000
+    gamma0 = 0.05
+    r0 = 10.
+    L = 50. * r0
+    rng = np.random.RandomState(8675309)
+    xl1 = (rng.random_sample(nlens)-0.5) * L
+    yl1 = (rng.random_sample(nlens)-0.5) * L
+    xs = (rng.random_sample(nsource)-0.5) * L
+    ys = (rng.random_sample(nsource)-0.5) * L
+    g1 = np.zeros( (nsource,) )
+    g2 = np.zeros( (nsource,) )
+    for x,y in zip(xl1,yl1):
+        dx = xs-x
+        dy = ys-y
+        r2 = dx**2 + dy**2
+        gammat = gamma0 * np.exp(-0.5*r2/r0**2)
+        g1 += -gammat * (dx**2-dy**2)/r2
+        g2 += -gammat * (2.*dx*dy)/r2
+
+    # Double the lenses
+    xl2 = np.concatenate([xl1,xl1])
+    yl2 = np.concatenate([yl1,yl1])
+
+    lens_cat1 = treecorr.Catalog(ra=xl1, dec=yl1, ra_units='arcmin', dec_units='arcmin')
+    lens_cat2 = treecorr.Catalog(ra=xl2, dec=yl2, ra_units='arcmin', dec_units='arcmin')
+    source_cat = treecorr.Catalog(ra=xs, dec=ys, g1=g1, g2=g2, ra_units='arcmin', dec_units='arcmin')
+    ng1 = treecorr.NGCorrelation(bin_size=0.1, min_sep=1., max_sep=20., sep_units='arcmin',
+                                 verbose=1)
+    ng2 = treecorr.NGCorrelation(bin_size=0.1, min_sep=1., max_sep=20., sep_units='arcmin',
+                                 verbose=1)
+
+    print('Normal case (all lenses different)')
+    ng1.process(lens_cat1, source_cat)
+    print('ng.xi = ',ng1.xi)
+    print('ng.xi_im = ',ng1.xi_im)
+
+    print('Duplicated lenses')
+    ng2.process(lens_cat2, source_cat)
+    print('ng.xi = ',ng2.xi)
+    print('ng.xi_im = ',ng2.xi_im)
+
+    print('Difference')
+    print('delta ng.xi = ',ng2.xi-ng1.xi)
+    print('delta ng.xi_im = ',ng2.xi_im-ng1.xi_im)
+    print('max diff = ',max(abs(ng2.xi - ng1.xi)))
+    np.testing.assert_allclose(ng2.xi, ng1.xi, rtol=1.e-3)
+    np.testing.assert_allclose(ng2.xi_im, ng1.xi_im, rtol=1.e-3, atol=1.e-6)
+
+
 if __name__ == '__main__':
     test_direct()
     test_direct_spherical()
@@ -1390,3 +1446,4 @@ def test_varxi():
     test_pieces()
     test_haloellip()
     test_varxi()
+    test_double()
diff --git a/tests/test_periodic.py b/tests/test_periodic.py
index 1c0b97ef..06a2fa35 100644
--- a/tests/test_periodic.py
+++ b/tests/test_periodic.py
@@ -310,7 +310,7 @@ def test_halotools():
 
     # TreeCorr is only a little faster than halotools for this, which makes sense, since
     # this is the exact calculation (bin_slop=0).  But with just a little slop, we get
-    # almost the same result bute faster.  This matters even more for larger N.
+    # almost the same result but faster.  This matters even more for larger N.
     corr = treecorr.NNCorrelation(min_sep=rbins[0], max_sep=rbins[-1], nbins=len(rbins)-1,
                                   period=Lbox, bin_slop=0.2)
     t4 = time.time()