Repository: manodeep/Corrfunc (language: python)

Issue:
Implement per-axis periodicity and box size
This PR allows the user to specify a periodic box length in each dimension, rather than assuming a cube.  Internally, Corrfunc already supports this, so the only necessary change was to add support for passing a Python tuple as the box size instead of a scalar.

This change is API and ABI backward-compatible (EDIT: no longer ABI backwards compatible, see https://github.com/manodeep/Corrfunc/pull/276#discussion_r882974028).  Currently, I've only implemented it for `theory.DD`, but all the other modules continue to work without modification.

@johannesulf, would you be willing to test this PR? I've tested that  `boxsize=L` and `boxsize=(L,L,L)` get the same answer, but I haven't actually verified the code gets the right answer for `Lx!=Ly!=Lz`!

Will close #247.
BUG DDrppi buggy for cross-corr in periodic spaces
<!--
Fill in the information below before opening an issue.
If you have a question, please consider emailing the mailing list:
  https://groups.google.com/forum/#!forum/corrfunc
-->

## General information

- Corrfunc version: 2.3.2
- platform: osx
- installation method (pip/source/other?): conda

## Issue description 

I am computing the pairs by hand and not getting what corrfunc reports.

### Expected behavior 

It should match the brute force computation.

### Actual behavior

It doesn't match it.

### What have you tried so far?

A bunch of stuff to try and see if it will miss one particle. It does not, so the bug is more subtle. It does not appear to be a bin edge thing.

## Minimal failing example

<!-- In this section, you should include or link to a code snippet that demonstrates the issue. -->

```python
import numpy as np
import Corrfunc

lbox = 120
zmax = 60
n_pi = int(zmax)
rmin = 40
rmax = 45
n_r = 2
rbins = np.array([0, rmin, rmax])

lbox_2 = lbox / 2

npts = 11
rng = np.random.RandomState(seed=42)
x1 = rng.uniform(low=0, high=lbox, size=npts)
y1 = rng.uniform(low=0, high=lbox, size=npts)
z1 = rng.uniform(low=0, high=lbox, size=npts)
w1 = rng.uniform(size=npts) * 0 + 1

nparts = 100
px1 = rng.uniform(low=0, high=lbox, size=nparts)
py1 = rng.uniform(low=0, high=lbox, size=nparts)
pz1 = rng.uniform(low=0, high=lbox, size=nparts)

print("\n")

dd_cf = Corrfunc.theory.DDrppi(
    False,
    1,
    zmax,
    rbins,
    x1,
    y1,
    z1,
    weights1=w1,
    X2=px1,
    Y2=py1,
    Z2=pz1,
    weights2=np.ones_like(px1),
    periodic=True,
    boxsize=lbox,
    weight_type="pair_product",
    verbose=True,
    isa="fallback",
)
dd_cf = dd_cf["weightavg"].reshape((n_r, n_pi)) * dd_cf["npairs"].reshape(
    (n_r, n_pi)
)

print("\n")
print("testing particles one by one...")

dd = np.zeros(n_r)
for hind, (hx, hy, hz) in enumerate(zip(x1, y1, z1)):

    for px, py, pz in zip(px1, py1, pz1):
        dx = np.abs(hx - px)
        dy = np.abs(hy - py)
        dz = np.abs(hz - pz)

        if dx > lbox_2:
            dx = lbox - dx
        if dy > lbox_2:
            dy = lbox - dy
        if dz > lbox_2:
            dz = lbox - dz

        if dz <= zmax:
            r = np.sqrt(dx * dx + dy * dy)

            if r < rmin:
                dd[0] += w1[hind]
                if dd_cf[0, int(dz)] == 0:
                    print(
                        "    def missed part:",
                        r,
                        dz,
                        int(dz),
                        dd_cf[1, int(dz)],
                    )

            elif r >= rmin and r < rmax:
                dd[1] += w1[hind]
                if dd_cf[1, int(dz)] == 0:
                    print(
                        "    def missed part:",
                        r,
                        dz,
                        int(dz),
                        dd_cf[1, int(dz)],
                    )

print("dd corrfunc:", np.sum(dd_cf, axis=1))
print("dd brute   :", dd)

dd_cf_sum = np.sum(dd_cf, axis=1)
assert np.array_equal(dd_cf_sum, dd)
```

This script outputs

```bash
(anl) clarence:Desktop beckermr$ python corrfunc_test.py 


[Warning] The CPU supports AVX2 but the compiler does not.  Can you try another compiler?
ND1 =           11 [xmin,ymin,zmin] = [2.470139,16.739263,5.574050], [xmax,ymax,zmax] = [114.085717,116.389182,94.221115]
ND2 =          100 [xmin,ymin,zmin] = [0.662654,0.607390,1.727219], [xmax,ymax,zmax] = [118.426432,118.278054,118.806462]
Running with points in [xmin,xmax] = 0.662654,118.426432 with periodic wrapping = 120.000000
Running with points in [ymin,ymax] = 0.607390,118.278054 with periodic wrapping = 120.000000
Running with points in [zmin,zmax] = 1.727219,118.806462 with periodic wrapping = 120.000000
In gridlink_double> Running with [nmesh_x, nmesh_y, nmesh_z]  = 5,5,2.  Time taken =   0.000 sec
countpairs_rp_pi_double> gridlink seems inefficient. nmesh = (5, 5, 2); avg_np = 0.22. Boosting bin refine factor - should lead to better performance
xmin = 0.662654 xmax=118.426432 rpmax = 45.000000
In gridlink_double> Running with [nmesh_x, nmesh_y, nmesh_z]  = 8,8,2.  Time taken =   0.000 sec
In gridlink_double> Running with [nmesh_x, nmesh_y, nmesh_z]  = 8,8,2.  Time taken =   0.000 sec
Using fallback kernel
0%.........10%.........20%.........30%.........40%.........50%.........60%.........70%.........80%.........90%.........100% done. Time taken =  0.001 secs


testing particles one by one...
    def missed part: 41.087356477185 48.94941175529149 48 0.0
dd corrfunc: [362.  93.]
dd brute   : [387.  96.]
Traceback (most recent call last):
  File "corrfunc_test.py", line 99, in <module>
    assert np.array_equal(dd_cf_sum, dd)
AssertionError
```


Discussion / hints:
Hello @lgarrison! Thanks for updating this PR. We checked the lines you've touched for [PEP 8](https://www.python.org/dev/peps/pep-0008) issues, and found:













There are currently no PEP 8 issues detected in this Pull Request. Cheers! :beers: 

##### Comment last updated at 2022-07-19 15:10:00 UTC
Fantastic! I'll check this PR against halotools in the next few days and report back.
I currently don't find agreement with halotools. Take the following code that counts pairs in spheres with radius of 10 in a uniformly sampled, non-cubic box.

```import numpy as np
import Corrfunc
from Corrfunc.theory.DD import DD
from halotools.mock_observables import tpcf

N = 100000
boxsize = (100.0, 100.0, 50.0)
nthreads = 4
autocorr = 1
seed = 42
np.random.seed(seed)
X = np.random.uniform(0, boxsize[0], N)
Y = np.random.uniform(0, boxsize[1], N)
Z = np.random.uniform(0, boxsize[2], N)
r_bins = np.array([1e-10, 10])
weights = np.ones_like(X)
results = DD(autocorr, nthreads, r_bins, X, Y, Z, weights1=weights,
             weight_type='pair_product', boxsize=boxsize, periodic=True)

pos = np.vstack([X, Y, Z]).T
xi = tpcf(pos, r_bins, period=boxsize)
print('Number of pairs...')
print('halotools: {:.0f}'.format(
    (xi + 1) * 4 * np.pi / 3 * r_bins[-1]**3 / np.prod(boxsize) * N**2))
print('Corrfunc: {:.0f}'.format(results['npairs'][0]))
print('Expectation: {:.0f}'.format(
    N**2 * 4 * np.pi / 3 * r_bins[-1]**3 / np.prod(boxsize)))
```

I get agreement with halotools and the expectation value (within error) for cubic boundary conditions. For the above code with a non-cubic box, Corrfunc currently does not agree with halotools. However, if I only consider points with Z<40, there is agreement with halotools. So it seems that the periodicity along the z-axis might not be taken into account by Corrfunc in the above code.
Awesome! Just tested the updated version and everything looks good.
The current PR doesn't implement non-cubic boundary conditions for Corrfunc.theory functions other than DD, right? Is it possible to also add it for the other pair counters, especially DDsmu?
Great! And yes, I haven't ported this feature to the other modules yet. It's not too bad, just involves echoing the `_countpairs.c` change to each module in that file (and updating the docstrings), and then updating each `theory/*/*impl.c.src` with the snippet from `theory/DD/countpairs_impl.c.src`. The last bit is updating the `Corrfunc/theory/*.py` interfaces (and docstrings).

@manodeep, can you review the change to DD before we proceed to the other modules?
I took one quick look and noted some comments, but I need to look more closely. Different boxsizes per dimension might impact the min-separation + periodic features, particularly in the calculation of the cell-pairs. 
I agree, if you could check those areas of the code, that would be great. Fortunately, there were few places in the code that used `boxsize`—most functions take `xdiff`, `ydiff`, and `zdiff` separately. And Corrfunc has always supported non-cubic periodic boxes whenever `boxsize=0` was used, although that code path was probably less tested.  It's also reassuring that non-cubic DD agrees with halotools exactly!
I've implemented the semantics from https://github.com/manodeep/Corrfunc/pull/276#discussion_r882985310. Specifically, a boxsize of -1 along any dimension now means that dimension is not periodic.  The surface area of the change was pretty minimal; it just affected the generation of the cell pairs.

I've tested that turning off periodicity per-dimension works by comparing to the periodic answer where the boxsize along that dimension has been extended beyond the particle distribution, such that it is effectively non-periodic. The answers agree exactly.

@manodeep, can you take a look?
Thanks for the review; I think we've uncovered an unintended behavior in the binsize calculation predating this PR. Can you take a look at https://github.com/manodeep/Corrfunc/pull/276#discussion_r888974272 and see if you agree?
Hey @manodeep, can you take a look at this when you get the chance?
I had looked at the comment and it seems fine. (Realised just now that I never "submitted" the comment - so it was pending ...)
Great! The code now uses the spatial extent instead of the boxsize, and all the tests are passing, so I think we're ready to port to the other modules. Will work on that next week.
@manodeep: xi and wp take a `double boxsize` argument in their C API, so they actually ignore `options->boxsize`. Do you think we should leave xi and wp alone in this PR, or should they also support non-cubic (periodic) boxes?  If so, do you think we should add `double boxsize_y, double boxsize_z` to the C function, or remove `double boxsize` and only use `options->boxsize`?
We could leave the C command line as is but change the shared library to use options->boxsize. This does result in a divergence of supported options from the command-line vs the shared library. 

Either way, if we support non-cubic boxes, then the normalization calculation within wp and xi will have to be updated. 


The other route would be to remove the boxsize argument from the command-line. In which case, perhaps it would call for a major release rather than a minor. 

(We have accumulated quite a few upgrades - so a major release could be in order)

Not sure if this answers your question - happy to discuss further to refine
I think I'd be okay with declaring that xi and wp only support cubic boxes. They are already restricted to periodic boxes, and we can point people to DD and DDrppi if they need the fully flexible, non-cubic and/or non-periodic options. It seems like overkill to change multiple levels of the API and add and test this new feature when there are good alternatives available!
I am fine with that - cubic periodic boxes should be the most common and the intended target for wp and xi anyway. 

Just so I understand properly - we will keep wp and xi as is and only add the non-cubic feature to the other pair counters


Yes, exactly!
@johannesulf This PR should now support per-axis boxsize and periodicity for the other modules, including `DDsmu`. Would you like to check that `DDsmu` works as expected before we merge?
> @johannesulf This PR should now support per-axis boxsize and periodicity for the other modules, including `DDsmu`. Would you like to check that `DDsmu` works as expected before we merge?

Yes, wonderful! Let me do a few checks and report back to you.
I performed some initial tests of `DDsmu` against halotools using per-axis periodicity and the results look good.
Excellent, thanks. @manodeep, good to merge?
I did some more in-depth tests and now I get some errors but only sometimes. Let me see whether I can find out what exactly triggers it.
@lgarrison I get errors when the number of tracers is very small. For example, if I take my code above and set `N=2`, it'll crash for me. I'm not sure whether this is supposed to happen. In principle, with a maximum search radius of 10 and box dimensions of 100, 100 and 50, no tracer should be counted twice.

`../../utils/gridlink_utils_double.c> ERROR:  nlattice = 2 is so small that with periodic wrapping the same cells will be counted twice ....exiting`
`../../utils/gridlink_utils_double.c> Please reduce Rmax = 10.000000 to be a smaller fraction of the particle distribution region = 13.333546`
`../../utils/gridlink_utils_double.c> ERROR:  nlattice = 1 is so small that with periodic wrapping the same cells will be counted twice ....exiting`
`../../utils/gridlink_utils_double.c> Please reduce Rmax = 10.000000 to be a smaller fraction of the particle distribution region = 0.001206`
`Received xstatus = 0 ystatus = 1 zstatus = 1. Error`
@johannesulf Corrfunc requires a minimum of 3 cells per dimension and so will crash based on the calculated number of cells -- [check is here](https://github.com/manodeep/Corrfunc/blob/master/utils/gridlink_utils.c.src#L37-L44). It looks like the actual particle extent is significantly smaller than 100 (or 50) - e.g., the two tracers are separated by ~13 Mpc in Y (-> only 2 cells) and ~1e-3 in Z (-> only 1 cell). 

You can avoid this crash by explicitly passing the boxsize. 


I just checked, the crash still happens if the boxsize is passed (as Johannes' above example does). But this makes sense, because we changed the gridding to only span the particle extent, whereas before it was using the boxsize. So now, regardless of the boxsize, we may have few cells if the particle extent is narrow.

I think the fix is actually pretty simple: only apply the 3 cell minimum if any pairs can be counted across the wrap; that is, if `xdiff + Rmax >= xwrap`. This is actually a little conservative, but I think will cover almost all practical cases.
Works much better. But I'm getting an out-of-bounds error when calculating an auto-correlation function with e.g. `Corrfunc.theory.DD ` using just 1 particle, even if the particle is within the bounds.
Wonderful! Everything looks good now.
Great! @manodeep do you want to check my last two commits?
Thanks @manodeep, I've added the warning and the parentheses as you suggested.

I'm not sure I follow the suggestion about the inverse xdiff. Wouldn't using the boxsize when xdiff is zero just create a lot of empty cells to process? And if the user requested automatic detection of the particle extent, then we don't even have a sensible non-zero value of the boxsize to use.

I think the code handles the planar case fine: once all the particles are gridded, it no longer matters that they were in a plane, or even what the size of their containing cell is.  The binsize calculation will still raise an error if the user tries to apply periodicity in that dimension that could cause repeat counts.
> 



> Thanks @manodeep, I've added the warning and the parentheses as you suggested.
> 
> I'm not sure I follow the suggestion about the inverse xdiff. Wouldn't using the boxsize when xdiff is zero just create a lot of empty cells to process? And if the user requested automatic detection of the particle extent, then we don't even have a sensible non-zero value of the boxsize to use.
> 

Agree that we need a single cell (i.e. cell index == 0) and the easiest way to do that is by setting inv_xyzdiff to 0. However, that might be confusing to read - may be add some lines of comments to show why inv_xyzdiff is being set to 0. That way, other users and later developers will not be confused. 

Yup - agree that when the particles are co-planar (or worse), we have a meaningful boxsize only in the periodic cases and not in the non-periodic ones. 

> I think the code handles the planar case fine: once all the particles are gridded, it no longer matters that they were in a plane, or even what the size of their containing cell is. The binsize calculation will still raise an error if the user tries to apply periodicity in that dimension that could cause repeat counts.

If there is only a single cell, then I am not sure what conditions needs to be checked - (xdiff + rmax) >= xwrap or (xdiff + rmax) >= 0.5*xwrap. Did you have a simple unit-test added that can check?


> If there is only a single cell, then I am not sure what conditions needs to be checked - (xdiff + rmax) >= xwrap or (xdiff + rmax) >= 0.5*xwrap. Did you have a simple unit-test added that can check?

I'm thinking of the condition as: "is this problem periodic?" (this previous condition was literally `options->periodic`). Which means that one edge of the particle extent can touch the other across the wrap. For that to happen, they need to be within `Rmax` of each other, or `xwrap - xdiff >= rmax`, which is equivalent to your first condition.

There's a unit test added here to make sure the code gives the right answer in the single-cell case (I've updated it slightly in the last commit): https://github.com/manodeep/Corrfunc/pull/276/files#diff-c7d4836bf94aa48cd9d4e4410ddab869be6ec08d83dfc0999ee5eac26cd8d30cR100

If we revisit #210, we might be able to relax this condition to allow a greater fraction of the boxsize. But I think this version is the most literal translation of the existing code to this PR.
@aphearin @manodeep I have been banging my head against this for a week. Here is my most minimal example of something I don't understand. Hopefully, it is my code!
Thanks for the bug report @beckermr - let me check whether I can reproduce it. I did have a couple of questions though - i) Does the bug still appear if you reduce zmax to a value smaller than Lbox/2, say 55 ? and ii) Can you reproduce the bug with Corrfunc installed from source/pip? (Related, I thought that the conda installed Corrfunc was called `corrfunc` - is that not the case?)
Just confirming that I can run your script and get the same output; thanks
for making the minimal reproducer!  Haven't had a chance to track down the
bug yet.

On Tue, Jan 28, 2020 at 3:04 PM Manodeep Sinha <notifications@github.com>
wrote:

> Thanks for the bug report @beckermr <https://github.com/beckermr> - let
> me check whether I can reproduce it. I did have a couple of questions
> though - i) Does the bug still appear if you reduce zmax to a value smaller
> than Lbox/2, say 55 ? and ii) Can you reproduce the bug with Corrfunc
> installed from source/pip? (Related, I thought that the conda installed
> Corrfunc was called corrfunc - is that not the case?)
>
> —
> You are receiving this because you are subscribed to this thread.
> Reply to this email directly, view it on GitHub
> <https://github.com/manodeep/Corrfunc/issues/210?email_source=notifications&email_token=ABLA7S6STYSQOINX2HARYGTRACFWLA5CNFSM4KMVBGQKYY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOEKEWPHQ#issuecomment-579430302>,
> or unsubscribe
> <https://github.com/notifications/unsubscribe-auth/ABLA7S3JT5UZ7TAWWHTCM4TRACFWLANCNFSM4KMVBGQA>
> .
>


-- 
Lehman Garrison <lgarrison@flatironinstitute.org>
Flatiron Research Fellow, Cosmology X Data Science Group
Center for Computational Astrophysics, Flatiron Institute
lgarrison.github.io

@manodeep yes I can reproduce the bug in both of those cases.
I suspect this is connected to the relative size of rmax and zmax to Lbox. For instance, if I change (xbin, ybin, zbin) refine-factors to `(1, 1, 1)`, then we have many mis-matches. And if I increase all the bin-refine-factors to >= 3, then the bug disappears. Since allowing such large relative separations is a new feature introduced with #189 (via commit https://github.com/manodeep/Corrfunc/commit/e53bce3050139337dae5723a3601e707f1f37920), may be a bug was introduced. 

@lgarrison At minimum, we will need a separate test case (from a brute force output) to test these sorts of large separations. 
If I disable the duplicate cell-checking in [here](https://github.com/manodeep/Corrfunc/blob/master/utils/gridlink_impl.c.src#L460-L462), then the bug disappears. However, we need to work through the logic and make sure we are not double-counting under any circumstances...

@manodeep is it possible that we we are missing particles when we do duplicate cell elimination because we are only considering one "wrapping" direction?  I think we need to consider both, the question is how to do so without double counting pairs.  Probably the answer lies in some restriction llike Rmax < CellSize/2...
Once this is resolved, do y'all have a brute force pair-counter test in the unit-testing to protect against regression? I find it's really hard to fool those kind of tests. Let me know in case that would be useful for me to pass along your way. 
@lgarrison That is definitely what is happening here. Consider zmax = Lbox/2, and we have 2 bins. In that case, the top half of `iz=0` can be within `zmax` of the bottom half of the `iz=1`, and the bottom half of the `iz=0` can be within the top half of `iz=1` cell *after periodic wrapping*. However, we will only associate (iz=0, iz=1) once and without periodic wrapping. Essentially the duplicate check has to needs to make sure both directions of the cell-pairs are included. For this specific scenario, no double-counting should occur but I am not sure about all the possible combinations of bin-refs, min-sep-opt, and various zmax/Lbox ratios ranging between [0.3-0.5] (i.e, between `2 <= nzbins <= 3`. I am reasonably confident that once `nzbins > 3`, Corrfunc works correctly. We would need to run the tests with the COMPREHENSIVE_INTEGRATION_TESTS` option enabled. 

@aphearin Thanks for the offer - we will absolutely need to add tests for large Rmax/zmax for **all** the pair-counters. 

@lgarrison In the meantime, should we revert the previous change (introduced in v2.3.2) for allowing large ratios of (Rmax, zmax)/Lboox?
@manodeep Yes, I think since we happen to be making a release right now anyways we should revert this change while we consider the best way to fix this.  Let's double-check with @beckermr that it fixes his problem before releasing!
@lgarrison If we revert the change, then Corrfunc will report an error for `(nx, ny, nzbins < 3)` and @beckermr will not be able to run this example

Fix patch:
diff --git a/CHANGES.rst b/CHANGES.rst
index fd0f774d..16f2f853 100644
--- a/CHANGES.rst
+++ b/CHANGES.rst
@@ -17,6 +17,7 @@ New features
 Enhancements
 ------------
 - Allow user to specify periodicity and box size per dimension [#276]
+- Allow larger Rmax (up to half the boxsize) [#277]
 
 Changes
 -------
@@ -28,6 +29,7 @@ Fixes
 -----
 - Add additional check to tell if it's safe to redirect stdout/err [#270]
 - Check and fix ``z`` vs ``cz`` in ``DDrppi_mocks`` and ``DDsmu_mocks`` only if comoving distance flag is not set [#275]
+- Update GNU assembler bug detection [#278]
 
 
 2.4.0 (2021-09-30)
diff --git a/common.mk b/common.mk
index ef128dd0..f098ba2b 100644
--- a/common.mk
+++ b/common.mk
@@ -347,12 +347,10 @@ ifeq ($(DO_CHECKS), 1)
     CLINK += -lm
   endif #not icc
 
-  # The GNU Assembler (GAS) has an AVX-512 bug in versions 2.30 to 2.31.1
-  # So we turn off AVX-512 if the compiler reports it is using one of these versions.
-  # This works for gcc and icc.
-  # clang typically uses its own assembler, but if it is using the system assembler, this will also detect that.
+  # The GNU Assembler (GAS) has an AVX-512 bug in some versions (originally 2.30 to 2.31.1)
+  # So we compile a test program and check the assembly for the correct output.
   # See: https://github.com/manodeep/Corrfunc/issues/193
-  GAS_BUG_DISABLE_AVX512 := $(shell $(CC) $(CFLAGS) -xc -Wa,-v -c /dev/null -o /dev/null 2>&1 | \grep -Ecm1 'GNU assembler version (2\.30|2\.31|2\.31\.1)')
+  GAS_BUG_DISABLE_AVX512 := $(shell echo 'vmovaps 64(,%rax), %zmm0' | $(CC) $(CFLAGS) -xassembler -c - -oa.out 2>/dev/null && objdump -dw a.out | \grep -q 'vmovaps 0x40(,%rax,1),%zmm0'; RET=$$?; rm -f a.out; echo $$RET)
 
   ifeq ($(GAS_BUG_DISABLE_AVX512),1)
     # Did the compiler support AVX-512 in the first place? Otherwise no need to disable it!
diff --git a/utils/gridlink_impl.c.src b/utils/gridlink_impl.c.src
index 92739de4..c8232b2a 100644
--- a/utils/gridlink_impl.c.src
+++ b/utils/gridlink_impl.c.src
@@ -460,7 +460,7 @@ struct cell_pair_DOUBLE * generate_cell_pairs_DOUBLE(struct cellarray_DOUBLE *la
         return NULL;
     }
     const int64_t num_self_pairs = totncells;
-    const int64_t num_nonself_pairs = totncells * max_ngb_cells / (1 + autocorr);
+    const int64_t num_nonself_pairs = totncells * max_ngb_cells;
 
     const int64_t max_num_cell_pairs = num_self_pairs + num_nonself_pairs;
     int64_t num_cell_pairs = 0;
@@ -529,7 +529,7 @@ struct cell_pair_DOUBLE * generate_cell_pairs_DOUBLE(struct cellarray_DOUBLE *la
                     //Check if the ngb-cell has already been added - can happen under periodic boundary
                     //conditions, with small value of nmesh_x/y/z (ie large Rmax relative to BoxSize)
                     if(check_for_duplicate_ngb_cells) {
-                        CHECK_AND_CONTINUE_FOR_DUPLICATE_NGB_CELLS_DOUBLE(icell, icell2, num_cell_pairs, num_ngb_this_cell, all_cell_pairs);
+                        CHECK_AND_CONTINUE_FOR_DUPLICATE_NGB_CELLS_DOUBLE(icell, icell2, off_xwrap, off_ywrap, off_zwrap, num_cell_pairs, num_ngb_this_cell, all_cell_pairs);
                     }
 
                     struct cellarray_DOUBLE *second = &(lattice2[icell2]);
diff --git a/utils/gridlink_mocks_impl.c.src b/utils/gridlink_mocks_impl.c.src
index fd4ae0ec..e09ab727 100644
--- a/utils/gridlink_mocks_impl.c.src
+++ b/utils/gridlink_mocks_impl.c.src
@@ -1585,7 +1585,12 @@ struct cell_pair_DOUBLE * generate_cell_pairs_mocks_theta_ra_dec_DOUBLE(cellarra
                         continue;
                     }
 
-                    CHECK_AND_CONTINUE_FOR_DUPLICATE_NGB_CELLS_DOUBLE(icell, icell2, num_cell_pairs, num_ngb_this_cell, all_cell_pairs);
+                    // No periodicity with mocks; wrap value is always zero.
+                    const DOUBLE xwrap = ZERO;
+                    const DOUBLE ywrap = ZERO;
+                    const DOUBLE zwrap = ZERO;
+
+                    CHECK_AND_CONTINUE_FOR_DUPLICATE_NGB_CELLS_DOUBLE(icell, icell2, xwrap, ywrap, zwrap, num_cell_pairs, num_ngb_this_cell, all_cell_pairs);
 
                     XRETURN(num_cell_pairs < max_num_cell_pairs, NULL,
                             "Error: Assigning this existing cell-pair would require accessing invalid memory.\n"
@@ -1626,6 +1631,10 @@ struct cell_pair_DOUBLE * generate_cell_pairs_mocks_theta_ra_dec_DOUBLE(cellarra
                     this_cell_pair->closest_y1 = closest_y1;
                     this_cell_pair->closest_z1 = closest_z1;
 
+                    this_cell_pair->xwrap = xwrap;
+                    this_cell_pair->ywrap = ywrap;
+                    this_cell_pair->zwrap = zwrap;
+
                     this_cell_pair->same_cell = (autocorr == 1 && icell2 == icell) ? 1:0;
 
                     num_cell_pairs++;
diff --git a/utils/gridlink_utils.c.src b/utils/gridlink_utils.c.src
index 24f68f40..593b414d 100644
--- a/utils/gridlink_utils.c.src
+++ b/utils/gridlink_utils.c.src
@@ -33,16 +33,24 @@ int get_binsize_DOUBLE(const DOUBLE xdiff, const DOUBLE xwrap, const DOUBLE rmax
 {
     int nmesh=(int)(refine_factor*xdiff/rmax);
     nmesh = nmesh < 1 ? 1:nmesh;
-    if( xwrap > 0 && ((xdiff + rmax) >= xwrap) ) {
+    
+    if(xwrap > 0 && rmax >= xwrap/2){
+        fprintf(stderr,"%s> ERROR: rmax=%f must be less than half of periodic boxize=%f to avoid double-counting particles\n",
+            __FILE__,rmax,xwrap);
+        return EXIT_FAILURE;
+    }
+
+    /*if( xwrap > 0 && ((xdiff + rmax) >= xwrap) ) {
         if (nmesh<(2*refine_factor+1))  {
             fprintf(stderr,"%s> ERROR:  nlattice = %d is so small that with periodic wrapping the same cells will be counted twice ....exiting\n",__FILE__,nmesh) ;
             fprintf(stderr,"%s> Please reduce Rmax = %"REAL_FORMAT" to be a smaller fraction of the particle distribution region = %"REAL_FORMAT"\n",
                     __FILE__,rmax, xdiff);
             return EXIT_FAILURE;
         }
-    }
+    }*/
 
     if (nmesh>max_ncells)  nmesh=max_ncells;
+    if (nmesh<2)  nmesh=2;  // to avoid forming duplicate cell pairs
     *xbinsize = xdiff/nmesh;
     *nlattice = nmesh;
 
diff --git a/utils/gridlink_utils.h.src b/utils/gridlink_utils.h.src
index d73ed36e..48645367 100644
--- a/utils/gridlink_utils.h.src
+++ b/utils/gridlink_utils.h.src
@@ -43,7 +43,7 @@ extern "C" {
                                                                    DOUBLE *sqr_sep_min, DOUBLE *sqr_sep_max);
 
 
-#define CHECK_AND_CONTINUE_FOR_DUPLICATE_NGB_CELLS_DOUBLE(icell, icell2, num_cell_pairs, num_ngb_this_cell, all_cell_pairs) { \
+#define CHECK_AND_CONTINUE_FOR_DUPLICATE_NGB_CELLS_DOUBLE(icell, icell2, off_xwrap, off_ywrap, off_zwrap, num_cell_pairs, num_ngb_this_cell, all_cell_pairs) { \
         int duplicate_flag = 0;                                         \
         XRETURN(num_cell_pairs - num_ngb_this_cell >= 0, NULL,          \
                 "Error: While working on detecting (potential) duplicate cell-pairs on primary cell = %"PRId64"\n" \
@@ -57,10 +57,16 @@ extern "C" {
                     "For cell-pair # %"PRId64", the primary cellindex (within cell-pair) = %"PRId64" should be *exactly* " \
                     "equal to current primary cellindex = %"PRId64". Num_cell_pairs = %"PRId64" num_ngb_this_cell = %"PRId64"\n", \
                     icell, num_cell_pairs - jj, this_cell_pair->cellindex1, icell, num_cell_pairs, num_ngb_this_cell); \
-            if (this_cell_pair->cellindex2 == icell2) {                 \
+            _Pragma("GCC diagnostic push")                              \
+            _Pragma("GCC diagnostic ignored \"-Wfloat-equal\"")         \
+            if (this_cell_pair->cellindex2 == icell2 &&                 \
+                this_cell_pair->xwrap == off_xwrap   &&                 \
+                this_cell_pair->ywrap == off_ywrap   &&                 \
+                this_cell_pair->zwrap == off_zwrap) {                   \
                 duplicate_flag = 1;                                     \
                 break;                                                  \
             }                                                           \
+            _Pragma("GCC diagnostic pop")                               \
         }                                                               \
         if(duplicate_flag == 1) continue;                               \
     }


Test patch:
diff --git a/Corrfunc/tests/test_theory.py b/Corrfunc/tests/test_theory.py
index 777796e6..6b94d60e 100644
--- a/Corrfunc/tests/test_theory.py
+++ b/Corrfunc/tests/test_theory.py
@@ -131,6 +131,104 @@ def test_narrow_extent(N, isa='fastest', nthreads=maxthreads()):
         assert np.all(results['npairs'] == [0, 0])
 
 
+@pytest.mark.parametrize('autocorr', [0, 1], ids=['cross', 'auto'])
+@pytest.mark.parametrize('binref', [1, 2, 3])
+@pytest.mark.parametrize('min_sep_opt', [False, True])
+@pytest.mark.parametrize('maxcells', [1, 2, 3])
+def test_duplicate_cellpairs(autocorr, binref, min_sep_opt, maxcells,
+                             isa='fastest', nthreads=maxthreads()):
+    '''A test to implement Manodeep's example from
+    https://github.com/manodeep/Corrfunc/pull/277#issuecomment-1190921894
+    where a particle pair straddles the wrap.
+
+    Also tests the large Rmax case, where Rmax is almost as large as L/2.
+    '''
+    from Corrfunc.theory import DD
+    boxsize = 432.
+    r_bins = np.array([0.01, 0.4]) * boxsize
+    pos = np.array([[0.02, 0.98],
+                    [0., 0.],
+                    [0., 0.0],
+                    ]) * boxsize
+
+    results = DD(autocorr, nthreads, r_bins, pos[0], pos[1], pos[2],
+                 X2=pos[0], Y2=pos[1], Z2=pos[2],
+                 boxsize=boxsize, periodic=True, isa=isa,
+                 verbose=True, max_cells_per_dim=maxcells,
+                 xbin_refine_factor=binref, ybin_refine_factor=binref,
+                 zbin_refine_factor=binref, enable_min_sep_opt=min_sep_opt,
+                 )
+
+    assert np.all(results['npairs'] == [2])
+
+    r_bins = np.array([0.2, 0.3, 0.49]) * boxsize
+    pos = np.array([[0., 0.],
+                    [0., 0.],
+                    [0., 0.48],
+                    ]) * boxsize
+
+    results = DD(autocorr, nthreads, r_bins, pos[0], pos[1], pos[2],
+                 X2=pos[0], Y2=pos[1], Z2=pos[2],
+                 boxsize=boxsize, periodic=True, isa=isa,
+                 verbose=True, max_cells_per_dim=maxcells,
+                 xbin_refine_factor=binref, ybin_refine_factor=binref,
+                 zbin_refine_factor=binref, enable_min_sep_opt=min_sep_opt,
+                 )
+
+    assert np.all(results['npairs'] == [0, 2])
+
+
+@pytest.mark.parametrize('autocorr', [0, 1], ids=['cross', 'auto'])
+@pytest.mark.parametrize('binref', [1, 2, 3])
+@pytest.mark.parametrize('min_sep_opt', [False, True])
+@pytest.mark.parametrize('maxcells', [1, 2, 3])
+def test_rmax_against_brute(autocorr, binref, min_sep_opt, maxcells,
+                            isa='fastest', nthreads=maxthreads()):
+    '''Generate two small point clouds near particles near (0,0,0)
+    and (L/2,L/2,L/2) and compare against the brute-force answer.
+
+    Use close to the max allowable Rmax, 0.49*Lbox.
+    '''
+    np.random.seed(1234)
+    npts = 100
+    eps = 0.2  # as a fraction of boxsize
+    boxsize = 123.
+    bins = np.linspace(0.01, 0.49*boxsize, 20)
+
+    # two clouds of width eps*boxsize
+    pos = np.random.uniform(low=-eps, high=eps,
+                            size=(3, npts))*boxsize
+    pos[:, npts//2:] += boxsize/2.  # second cloud is in center of box
+    pos %= boxsize
+
+    # Compute the pairwise distance between particles with broadcasting.
+    # Broadcasting (3,1,npts) against (3,npts,1) yields (3,npts,npts),
+    # which is the array of all npts^2 (x,y,z) differences.
+    pdiff = np.abs(pos[:, np.newaxis] - pos[:, :, np.newaxis])
+    pdiff[pdiff >= boxsize/2] -= boxsize
+
+    # Compute dist^2 = x^2 + y^2 + z^2, flattening because we don't
+    # care which dist came from which particle pair
+    sqr_pdiff = (pdiff**2).sum(axis=0).reshape(-1)
+    brutecounts, _ = np.histogram(sqr_pdiff, bins=bins**2)
+
+    # spot-check that we have non-zero counts
+    assert np.any(brutecounts > 0)
+
+    from Corrfunc.theory import DD
+
+    results = DD(autocorr, nthreads, bins, pos[0], pos[1], pos[2],
+                 X2=pos[0], Y2=pos[1], Z2=pos[2],
+                 boxsize=boxsize, periodic=True, isa=isa,
+                 verbose=True, max_cells_per_dim=maxcells,
+                 xbin_refine_factor=binref, ybin_refine_factor=binref,
+                 zbin_refine_factor=binref, enable_min_sep_opt=min_sep_opt,
+                 )
+
+    assert np.all(results['npairs'] == brutecounts)
+
+
+
 @pytest.mark.parametrize('isa,nthreads', generate_isa_and_nthreads_combos())
 def test_DDrppi(gals_Mr19, isa, nthreads):
     from Corrfunc.theory import DDrppi