From 226a5cf41d17cec1b84157d4e14e950329295810 Mon Sep 17 00:00:00 2001
From: Manodeep Sinha <manodeep@gmail.com>
Date: Wed, 7 Jun 2017 05:56:33 +1000
Subject: [PATCH 01/15] For issue #763

---
 .../cpairs/npairs_s_mu_engine.pyx             | 60 +++++++++++--------
 1 file changed, 34 insertions(+), 26 deletions(-)

diff --git a/halotools/mock_observables/pair_counters/cpairs/npairs_s_mu_engine.pyx b/halotools/mock_observables/pair_counters/cpairs/npairs_s_mu_engine.pyx
index 0049b6cb1..1fcc1bdb7 100644
--- a/halotools/mock_observables/pair_counters/cpairs/npairs_s_mu_engine.pyx
+++ b/halotools/mock_observables/pair_counters/cpairs/npairs_s_mu_engine.pyx
@@ -8,7 +8,7 @@ cimport cython
 from libc.math cimport ceil
 from libc.math cimport sqrt
 
-__author__ = ('Andrew Hearin', 'Duncan Campbell')
+__author__ = ('Andrew Hearin', 'Duncan Campbell', 'Manodeep Sinha')
 __all__ = ('npairs_s_mu_engine', )
 
 @cython.boundscheck(False)
@@ -50,8 +50,9 @@ def npairs_s_mu_engine(double_mesh, x1in, y1in, z1in, x2in, y2in, z2in,
         separated by a distance less than the corresponding entry of ``s_bins``.
 
     """
-    cdef cnp.float64_t[:] s_bins = s_bins_in
-    cdef cnp.float64_t[:] mu_bins = mu_bins_in
+    cdef cnp.float64_t[:] sqr_s_bins = s_bins_in * s_bins_in
+    cdef cnp.float64_t[:] sqr_mu_bins = mu_bins_in * mu_bins_in
+
     cdef cnp.float64_t xperiod = double_mesh.xperiod
     cdef cnp.float64_t yperiod = double_mesh.yperiod
     cdef cnp.float64_t zperiod = double_mesh.zperiod
@@ -60,8 +61,8 @@ def npairs_s_mu_engine(double_mesh, x1in, y1in, z1in, x2in, y2in, z2in,
     cdef int PBCs = double_mesh._PBCs
 
     cdef int Ncell1 = double_mesh.mesh1.ncells
-    cdef int num_s_bins = len(s_bins)
-    cdef int num_mu_bins = len(mu_bins)
+    cdef int num_s_bins = len(sqr_s_bins)
+    cdef int num_mu_bins = len(sqr_mu_bins)
     cdef cnp.int64_t[:,:] counts = np.zeros((num_s_bins, num_mu_bins), dtype=np.int64)
     cdef cnp.int64_t[:,:] counts_sum = np.zeros((num_s_bins, num_mu_bins), dtype=np.int64)
 
@@ -105,7 +106,8 @@ def npairs_s_mu_engine(double_mesh, x1in, y1in, z1in, x2in, y2in, z2in,
     cdef cnp.float64_t x2shift, y2shift, z2shift, dx, dy, dz, dxy_sq, dz_sq
     cdef cnp.float64_t x1tmp, y1tmp, z1tmp, s, mu
     cdef int Ni, Nj, i, j, k, l, g, max_k
-    cdef cnp.float64_t s_max = np.max(s_bins_in), mu_max = np.max(mu_bins_in)
+    # cdef cnp.float64_t s_max = np.max(s_bins_in), mu_max = np.max(mu_bins_in)
+    cdef cnp.float64_t sqr_s_max = np.max(sqr_s_bins), sqr_mu_max = np.max(sqr_mu_bins)
 
     cdef cnp.float64_t[:] x_icell1, x_icell2
     cdef cnp.float64_t[:] y_icell1, y_icell2
@@ -188,26 +190,32 @@ def npairs_s_mu_engine(double_mesh, x1in, y1in, z1in, x2in, y2in, z2in,
                                     dz_sq = dz*dz
 
                                     #transform to s and mu
-                                    s = sqrt(dz_sq + dxy_sq)
-                                    if s!=0:
-                                        mu = sqrt(dxy_sq)/s
-                                    else:
-                                        mu=0.0
-
-                                    if (s <= s_max) & (mu <= mu_max):
-
-                                        k = num_s_bins-2
-                                        while k!=-1:
-                                            if s > s_bins[k]: break
-                                            k=k-1
-
-                                        g = num_mu_bins-2
-                                        while g!=-1:
-                                            if mu > mu_bins[g]: break
-                                            g=g-1
-
-                                        # Only counts pairs in that bin.
-                                        counts[k+1,g+1] += 1
+                                    sqr_s = dz_sq + dxy_sq
+
+                                    if sqr_s > sqr_s_max:
+                                        continue
+
+                                    sqr_mu = sqr_s > 0 ? dxy_sq/sqr_s : 0.0
+                                    if sqr_mu > sqr_mu_max:
+                                        continue
+
+                                    # MS: I have intentionally split
+                                    # up the loop. Since division
+                                    # is slow, computing mu is a bottle-neck.
+                                    # Computing the 's' bin however can proceed
+                                    # in the meantime.
+                                    k = num_s_bins-2
+                                    while k!=-1:
+                                        if sqr_s > sqr_s_bins[k]: break
+                                        k=k-1
+
+                                    g = num_mu_bins-2
+                                    while g!=-1:
+                                        if sqr_mu > sqr_mu_bins[g]: break
+                                        g=g-1
+
+                                    # Only counts pairs in that bin.
+                                    counts[k+1,g+1] += 1
 
     # Adds counts for all bins where s < s_bin and mu < mu_bin.
     for k in range(num_s_bins):

From 89ab24adeb3e3bb46ab73620c0e7c8e5ca63fcfc Mon Sep 17 00:00:00 2001
From: Manodeep Sinha <manodeep@gmail.com>
Date: Wed, 7 Jun 2017 06:03:54 +1000
Subject: [PATCH 02/15] Removed commented out line

---
 .../mock_observables/pair_counters/cpairs/npairs_s_mu_engine.pyx | 1 -
 1 file changed, 1 deletion(-)

diff --git a/halotools/mock_observables/pair_counters/cpairs/npairs_s_mu_engine.pyx b/halotools/mock_observables/pair_counters/cpairs/npairs_s_mu_engine.pyx
index 1fcc1bdb7..6e8f2add7 100644
--- a/halotools/mock_observables/pair_counters/cpairs/npairs_s_mu_engine.pyx
+++ b/halotools/mock_observables/pair_counters/cpairs/npairs_s_mu_engine.pyx
@@ -106,7 +106,6 @@ def npairs_s_mu_engine(double_mesh, x1in, y1in, z1in, x2in, y2in, z2in,
     cdef cnp.float64_t x2shift, y2shift, z2shift, dx, dy, dz, dxy_sq, dz_sq
     cdef cnp.float64_t x1tmp, y1tmp, z1tmp, s, mu
     cdef int Ni, Nj, i, j, k, l, g, max_k
-    # cdef cnp.float64_t s_max = np.max(s_bins_in), mu_max = np.max(mu_bins_in)
     cdef cnp.float64_t sqr_s_max = np.max(sqr_s_bins), sqr_mu_max = np.max(sqr_mu_bins)
 
     cdef cnp.float64_t[:] x_icell1, x_icell2

From 63e8914a5378a8b946181a13f85756d3741e664a Mon Sep 17 00:00:00 2001
From: Manodeep Sinha <manodeep@gmail.com>
Date: Wed, 7 Jun 2017 06:32:45 +1000
Subject: [PATCH 03/15] Ternary operators in Cython are the python kind, not
 the C kind

---
 .../pair_counters/cpairs/npairs_s_mu_engine.pyx                 | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/halotools/mock_observables/pair_counters/cpairs/npairs_s_mu_engine.pyx b/halotools/mock_observables/pair_counters/cpairs/npairs_s_mu_engine.pyx
index 6e8f2add7..646c769ce 100644
--- a/halotools/mock_observables/pair_counters/cpairs/npairs_s_mu_engine.pyx
+++ b/halotools/mock_observables/pair_counters/cpairs/npairs_s_mu_engine.pyx
@@ -194,7 +194,7 @@ def npairs_s_mu_engine(double_mesh, x1in, y1in, z1in, x2in, y2in, z2in,
                                     if sqr_s > sqr_s_max:
                                         continue
 
-                                    sqr_mu = sqr_s > 0 ? dxy_sq/sqr_s : 0.0
+                                    sqr_mu = (dxy_sq/sqr_s if sqr_s > 0.0 else 0.0)
                                     if sqr_mu > sqr_mu_max:
                                         continue
 

From dde85da3f338db281c186a8528c49ba998edf014 Mon Sep 17 00:00:00 2001
From: Duncan Campbell <duncan.alan.campbell@gmail.com>
Date: Fri, 9 Jun 2017 15:56:55 -0400
Subject: [PATCH 04/15] added brute force s and mu pair counter

---
 .../mock_observables/pair_counters/pairs.py   | 145 +++++++++++++++++-
 1 file changed, 144 insertions(+), 1 deletion(-)

diff --git a/halotools/mock_observables/pair_counters/pairs.py b/halotools/mock_observables/pair_counters/pairs.py
index a69e945a3..544add5c4 100755
--- a/halotools/mock_observables/pair_counters/pairs.py
+++ b/halotools/mock_observables/pair_counters/pairs.py
@@ -352,6 +352,96 @@ def xy_z_wnpairs(sample1, sample2, rp_bins, pi_bins, period=None, weights1=None,
     return n
 
 
+def s_mu_npairs(sample1, sample2, s_bins, mu_bins, period=None):
+    """
+    Calculate the number of pairs with 3D radial separations less than or equal to
+    :math:`s`, and angular seperations along the LOS, :math:`\mu=\cos(\theta_{\rm LOS})`.
+
+    Assumes the first N-1 dimensions are perpendicular to the line-of-sight (LOS), and
+    the final dimension is parallel to the LOS.
+
+    Parameters
+    ----------
+    sample1 : array_like
+        N by k numpy array of k-dimensional positions. Should be between zero and
+        period
+
+    sample2 : array_like
+        N by k numpy array of k-dimensional positions. Should be between zero and
+        period
+
+    s_bins : array_like
+        numpy array of shape (num_s_bin_edges, ) storing the :math:`s`
+        boundaries defining the bins in which pairs are counted.
+
+    mu_bins : array_like
+        numpy array of shape (num_mu_bin_edges, ) storing the
+        :math:`\cos(\theta_{\rm LOS})` boundaries defining the bins in
+        which pairs are counted. All values must be between [0,1].
+
+    period : array_like, optional
+        length k array defining  periodic boundary conditions. If only
+        one number, Lbox, is specified, period is assumed to be np.array([Lbox]*k).
+        If none, PBCs are set to infinity.
+
+    Returns
+    -------
+    N_pairs : ndarray of shape (num_s_bin_edges, num_mu_bin_edges) storing the 
+        number counts of pairs with seperations less than ``s_bins`` and ``mu_bins``
+    
+    Notes
+    -----
+    Along the first dimension of ``N_pairs``, :math:`s` (the radial seperation) increases.
+    Along the second dimension,  :math:`\mu` (the cosine of :math:`\theta_{\rm LOS}`) 
+    decreases, i.e. :math:`\theta_{\rm LOS}` increases.
+    """
+
+    sample1 = np.atleast_2d(sample1)
+    sample2 = np.atleast_2d(sample2)
+    s_bins = np.atleast_1d(s_bins)
+    mu_bins = np.atleast_1d(mu_bins)
+
+    # Check to make sure both data sets have the same dimension. Otherwise, throw an error!
+    if np.shape(sample1)[-1] != np.shape(sample2)[-1]:
+        raise ValueError("sample1 and sample2 inputs do not have the same dimension.")
+        return None
+
+    # Process period entry and check for consistency.
+    if period is None:
+        period = np.array([np.inf]*np.shape(sample1)[-1])
+    else:
+        period = np.asarray(period).astype("float64")
+        if np.shape(period) == ():
+            period = np.array([period]*np.shape(sample1)[-1])
+        elif np.shape(period)[0] != np.shape(sample1)[-1]:
+            raise ValueError("period should have len == dimension of points")
+            return None
+    
+    # create N1 x N2 x 2 array to store **all** pair separation distances
+    # note that this array can be very large for large N1 and N2
+    N1 = len(sample1)
+    N2 = len(sample2)
+    dd = np.zeros((N1*N2, 2))
+    
+    # calculate distance between every point and every other point
+    for i in range(0, N1):
+        x1 = sample1[i, :]
+        x2 = sample2
+        dd[i*N2:i*N2+N2, 0] = distance(x1, x2, period)
+        dd[i*N2:i*N2+N2, 1] = np.cos(theta_LOS(x1, x2, period))
+    
+    # put mu bins in increasing theta_LOS order
+    mu_bins = np.sort(mu_bins)[::-1]
+    
+    # bin distances in s and mu bins
+    n = np.zeros((s_bins.size, mu_bins.size), dtype=np.int)
+    for i in range(s_bins.size):
+        for j in range(mu_bins.size):
+            n[i, j] = np.sum((dd[:, 0] <= s_bins[i]) & (dd[:, 1] >= mu_bins[j]))
+    
+    return n
+
+
 def distance(x1, x2, period=None):
     """
     Find the Euclidean distance between x1 & x2, accounting for box periodicity.
@@ -371,7 +461,6 @@ def distance(x1, x2, period=None):
     Returns
     -------
     distance : array
-
     """
 
     x1 = np.atleast_2d(x1)
@@ -477,3 +566,57 @@ def perpendicular_distance(x1, x2, period=None):
     distance = np.sqrt(np.sum(m*m, axis=len(np.shape(m))-1))
 
     return distance
+
+
+def theta_LOS(x1, x2, period=None):
+    """
+    Find the seperation angle from the LOS between x1 & x2, accounting for box periodicity.
+
+    Assumes the first N-1 dimensions are perpendicular to the line-of-sight (LOS).
+
+    Parameters
+    ----------
+    x1 : array_like
+        N by k numpy array of k-dimensional positions. Should be between zero and period
+
+    x2 : array_like
+        N by k numpy array of k-dimensional positions. Should be between zero and period.
+
+    period : array_like
+        Size of the simulation box along each dimension. Defines periodic boundary
+        conditioning.  Must be axis aligned.
+
+    Returns
+    -------
+    theta_LOS : array
+        angle from LOS in radians
+        
+    Notes
+    -----
+    theta_LOS is set to 0.0 if the distance between points is 0.0
+    """
+
+    x1 = np.atleast_2d(x1)
+    x2 = np.atleast_2d(x2)
+    if period is None:
+        period = np.array([np.inf]*np.shape(x1)[-1])
+
+    # check for consistency
+    if np.shape(x1)[-1] != np.shape(x2)[-1]:
+        raise ValueError("x1 and x2 list of points must have same dimension k.")
+    else:
+        k = np.shape(x1)[-1]
+    if np.shape(period)[0] != np.shape(x1)[-1]:
+        raise ValueError("period must have length equal to the dimension of x1 and x2.")
+    
+    r_perp = perpendicular_distance(x1, x2, period=period)
+    r_parallel = parallel_distance(x1, x2, period=period)
+    
+    # deal with zero seperation
+    r = np.sqrt(r_perp**2 + r_parallel**2)
+    mask = (r>0.0)
+    
+    theta = np.zeros(len(r)) # set to zero if r==0
+    theta[mask] = np.pi/2.0 - np.arctan2(r_parallel[mask],r_perp[mask])
+
+    return theta

From 7c6e388a513e18514557bdf411bb89a32ddde998 Mon Sep 17 00:00:00 2001
From: Duncan Campbell <duncan.alan.campbell@gmail.com>
Date: Fri, 9 Jun 2017 16:04:36 -0400
Subject: [PATCH 05/15] added non-tirival s mu pair tests

---
 .../test_pair_counters/test_npairs_s_mu.py    | 98 ++++++++++++++++---
 1 file changed, 83 insertions(+), 15 deletions(-)

diff --git a/halotools/mock_observables/pair_counters/test_pair_counters/test_npairs_s_mu.py b/halotools/mock_observables/pair_counters/test_pair_counters/test_npairs_s_mu.py
index ebf31d028..61c2f9387 100755
--- a/halotools/mock_observables/pair_counters/test_pair_counters/test_npairs_s_mu.py
+++ b/halotools/mock_observables/pair_counters/test_pair_counters/test_npairs_s_mu.py
@@ -8,6 +8,7 @@
 from ..npairs_s_mu import npairs_s_mu
 from ....mock_observables import npairs_3d
 # load comparison simple pair counters
+from ..pairs import s_mu_npairs as pure_python_brute_force_npairs_s_mu
 
 from astropy.tests.helper import pytest
 from astropy.utils.misc import NumpyRNGContext
@@ -15,7 +16,8 @@
 slow = pytest.mark.slow
 fixed_seed = 43
 
-__all__ = ('test_npairs_s_mu_periodic', 'test_npairs_s_mu_nonperiodic')
+__all__ = ('test_npairs_s_mu_periodic', 'test_npairs_s_mu_nonperiodic',
+           'test_npairs_s_mu_point_surrounded_by_sphere')
 
 # set up random points to test pair counters
 Npts = 1000
@@ -55,26 +57,40 @@ def test_npairs_s_mu_periodic():
     """
     test npairs_s_mu with periodic boundary conditions.
     """
-
+    
+    # define bins
     s_bins = np.array([0.0, 0.1, 0.2, 0.3])
     N_mu_bins = 100
     mu_bins = np.linspace(0, 1.0, N_mu_bins)
     Npts = len(random_sample)
-
+    
+    # count pairs using optimized double tree pair counter
     result = npairs_s_mu(random_sample, random_sample, s_bins, mu_bins,
         period=period, num_threads=num_threads)
-
+    
     msg = 'The returned result is an unexpected shape.'
     assert np.shape(result) == (len(s_bins), N_mu_bins), msg
-
-    result = np.diff(result, axis=1)
-    result = np.sum(result, axis=1) + Npts
-
+    
+    # count pairs using 3D pair counter
     test_result = npairs_3d(random_sample, random_sample, s_bins,
         period=period, num_threads=num_threads)
-
+    
+    # summing pairs counts along mu axis should match the 3D radial pair result
+    result = np.diff(result, axis=1)
+    result = np.sum(result, axis=1) + Npts
+    
     msg = "The double tree's result(s) are not equivalent to simple pair counter's."
     assert np.all(result == test_result), msg
+    
+    # compare to brute force s and mu pair counter on a set of random points
+    result = npairs_s_mu(random_sample, random_sample, s_bins, mu_bins,
+        period=period, num_threads=num_threads)
+    
+    test_result = pure_python_brute_force_npairs_s_mu(random_sample, random_sample,
+        s_bins, mu_bins, period=period)
+    
+    msg = "The double tree's result(s) are not equivalent to simple pair counter's."
+    assert np.all(test_result == result), msg
 
 
 def test_npairs_s_mu_nonperiodic():
@@ -82,20 +98,72 @@ def test_npairs_s_mu_nonperiodic():
     test npairs_s_mu without periodic boundary conditions.
     """
 
+    # define bins
     s_bins = np.array([0.0, 0.1, 0.2, 0.3])
     N_mu_bins = 100
     mu_bins = np.linspace(0, 1.0, N_mu_bins)
-
+    Npts = len(random_sample)
+    
+    # count pairs using optimized double tree pair counter
     result = npairs_s_mu(random_sample, random_sample, s_bins, mu_bins,
         num_threads=num_threads)
-
+    
     msg = 'The returned result is an unexpected shape.'
     assert np.shape(result) == (len(s_bins), N_mu_bins), msg
-
+    
+    # count pairs using 3D pair counter
+    test_result = npairs_3d(random_sample, random_sample, s_bins,
+        num_threads=num_threads)
+    
+    # summing pairs counts along mu axis should match the 3D radial pair result
     result = np.diff(result, axis=1)
     result = np.sum(result, axis=1) + Npts
-
-    test_result = npairs_3d(random_sample, random_sample, s_bins, num_threads=num_threads)
-
+    
     msg = "The double tree's result(s) are not equivalent to simple pair counter's."
     assert np.all(result == test_result), msg
+    
+    # compare to brute force s and mu pair counter on a set of random points
+    result = npairs_s_mu(random_sample, random_sample, s_bins, mu_bins,
+        num_threads=num_threads)
+    
+    test_result = pure_python_brute_force_npairs_s_mu(random_sample, random_sample,
+        s_bins, mu_bins)
+    
+    msg = "The double tree's result(s) are not equivalent to simple pair counter's."
+    assert np.all(test_result == result), msg
+
+
+def test_npairs_s_mu_point_surrounded_by_sphere():
+    """
+    test npairs_s_mu without periodic boundary conditions on a point surrounded by an
+    evenly distribued circle of points perpdindicular to the LOS.
+    """
+    
+    # put one point per degree in the circle
+    theta = np.linspace(0,2*np.pi-2*np.pi/360,360)+(2.0*np.pi)/(360*2)
+    r = 1.5
+    x = r*np.cos(theta)
+    y = r*np.zeros(360)
+    z = r*np.sin(theta)
+    sample2 = np.vstack((x,y,z)).T
+    
+    # find pairs between the circle and the origin
+    sample1 = np.array([[0,0,0]])
+    
+    # define bins
+    theta_bins = np.linspace(0.0,np.pi/2.0,91) # make 1 degree mu bins
+    mu_bins = np.sort(np.cos(theta_bins))
+    s_bins = np.array([0,1,2])
+    
+    # count pairs using optimized double tree pair counter
+    result = npairs_s_mu(sample1, sample2, s_bins, mu_bins, num_threads=num_threads)
+    pairs = np.diff(np.diff(result, axis=0),axis=1)
+    
+    # note that there should be 4 pairs per mu bin
+    # since each quadrant of te circle counts once for each mu bin
+    # and there is one point per degree
+    
+    msg = "The number of pairs between the origin and a circle with one point per degree is incorrect."
+    assert np.all(pairs[1,:] == 4), msg
+
+

From 1b0ea51a930d388924f33824918a0674893c6605 Mon Sep 17 00:00:00 2001
From: Duncan Campbell <duncan.alan.campbell@gmail.com>
Date: Fri, 9 Jun 2017 16:05:04 -0400
Subject: [PATCH 06/15] added non-tirival s mu pair tests

---
 .../pair_counters/test_pair_counters/test_npairs_s_mu.py      | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/halotools/mock_observables/pair_counters/test_pair_counters/test_npairs_s_mu.py b/halotools/mock_observables/pair_counters/test_pair_counters/test_npairs_s_mu.py
index 61c2f9387..7bf1e694e 100755
--- a/halotools/mock_observables/pair_counters/test_pair_counters/test_npairs_s_mu.py
+++ b/halotools/mock_observables/pair_counters/test_pair_counters/test_npairs_s_mu.py
@@ -17,7 +17,7 @@
 fixed_seed = 43
 
 __all__ = ('test_npairs_s_mu_periodic', 'test_npairs_s_mu_nonperiodic',
-           'test_npairs_s_mu_point_surrounded_by_sphere')
+           'test_npairs_s_mu_point_surrounded_by_circle')
 
 # set up random points to test pair counters
 Npts = 1000
@@ -133,7 +133,7 @@ def test_npairs_s_mu_nonperiodic():
     assert np.all(test_result == result), msg
 
 
-def test_npairs_s_mu_point_surrounded_by_sphere():
+def test_npairs_s_mu_point_surrounded_by_circle():
     """
     test npairs_s_mu without periodic boundary conditions on a point surrounded by an
     evenly distribued circle of points perpdindicular to the LOS.

From 46b920ee342ee0dd381b4c3cefb03d82c06de773 Mon Sep 17 00:00:00 2001
From: Duncan Campbell <duncan.alan.campbell@gmail.com>
Date: Fri, 9 Jun 2017 16:16:08 -0400
Subject: [PATCH 07/15] updated docstring and transformed mu_bins internally to
 sine convention which is used by the engine.

---
 .../pair_counters/npairs_s_mu.py              | 22 ++++++++++---------
 1 file changed, 12 insertions(+), 10 deletions(-)

diff --git a/halotools/mock_observables/pair_counters/npairs_s_mu.py b/halotools/mock_observables/pair_counters/npairs_s_mu.py
index 6cc184c2a..685c6a9fd 100644
--- a/halotools/mock_observables/pair_counters/npairs_s_mu.py
+++ b/halotools/mock_observables/pair_counters/npairs_s_mu.py
@@ -21,7 +21,7 @@ def npairs_s_mu(sample1, sample2, s_bins, mu_bins, period=None,
         verbose=False, num_threads=1, approx_cell1_size=None, approx_cell2_size=None):
     r"""
     Function counts the number of pairs of points separated by less than
-    radial separation, *s,* and :math:`\mu\equiv\sin(\theta_{\rm los})`,
+    radial separation, *s,* and :math:`\mu\equiv\cos(\theta_{\rm los})`,
     where :math:`\theta_{\rm los}` is the line-of-sight angle
     between points and :math:`s^2 = r_{\rm parallel}^2 + r_{\rm perp}^2`.
 
@@ -99,15 +99,13 @@ def npairs_s_mu(sample1, sample2, s_bins, mu_bins, period=None,
         number of pairs separated by less than (s, mu)
 
     Notes
-    ------
-    The quantity :math:`\mu` is defined as the :math:`\sin(\theta_{\rm los})`
-    and not the conventional :math:`\cos(\theta_{\rm los})`. This is
-    because the pair counter has been optimized under the assumption that its
-    separation variable (in this case, :math:`\mu`) *increases*
-    as :math:`\theta_{\rm los})` increases.
-
+    -----
+    Along the first dimension of ``num_pairs``, :math:`s` (the radial seperation) increases.
+    Along the second dimension,  :math:`\mu` (the cosine of :math:`\theta_{\rm LOS}`) 
+    decreases, i.e. :math:`\theta_{\rm LOS}` increases.
+    
     One final point of clarification concerning double-counting may be in order.
-    Suppose sample1==sample2 and rbins[0]==0. Then the returned value for this bin
+    Suppose sample1==sample2 and s_bins[0]==0. Then the returned value for this bin
     will be len(sample1), since each sample1 point has distance 0 from itself.
 
     Examples
@@ -158,6 +156,10 @@ def npairs_s_mu(sample1, sample2, s_bins, mu_bins, period=None,
         msg = ("\n Input `mu_bins` must be a monotonically increasing \n"
                "1D array with at least two entries")
         raise ValueError(msg)
+    # convert to mu=sin(theta_los) binning used by the cython engine.
+    mu_bins_prime = np.sin(np.arccos(mu_bins))
+    mu_bins_prime = np.sort(mu_bins_prime)
+    # increasing mu_prime now corresponds to increasing theta_LOS
 
     search_xlength, search_ylength, search_zlength = rmax, rmax, rmax
 
@@ -176,7 +178,7 @@ def npairs_s_mu(sample1, sample2, s_bins, mu_bins, period=None,
 
     # # Create a function object that has a single argument, for parallelization purposes
     engine = partial(npairs_s_mu_engine,
-        double_mesh, x1in, y1in, z1in, x2in, y2in, z2in, s_bins, mu_bins)
+        double_mesh, x1in, y1in, z1in, x2in, y2in, z2in, s_bins, mu_bins_prime)
 
     # # Calculate the cell1 indices that will be looped over by the engine
     num_threads, cell1_tuples = _cell1_parallelization_indices(

From 6f400388ee3707eb674e0d40ac4eda6944c78f2d Mon Sep 17 00:00:00 2001
From: Duncan Campbell <duncan.alan.campbell@gmail.com>
Date: Fri, 9 Jun 2017 16:26:52 -0400
Subject: [PATCH 08/15] removed conversion to mu sine convention and updated
 doc string

---
 .../two_point_clustering/s_mu_tpcf.py         | 39 +++++++------------
 1 file changed, 15 insertions(+), 24 deletions(-)

diff --git a/halotools/mock_observables/two_point_clustering/s_mu_tpcf.py b/halotools/mock_observables/two_point_clustering/s_mu_tpcf.py
index d4bdae9d0..f7f9748a6 100644
--- a/halotools/mock_observables/two_point_clustering/s_mu_tpcf.py
+++ b/halotools/mock_observables/two_point_clustering/s_mu_tpcf.py
@@ -61,9 +61,6 @@ def s_mu_tpcf(sample1, s_bins, mu_bins, sample2=None, randoms=None,
         :math:`\cos(\theta_{\rm LOS})` boundaries defining the bins in
         which pairs are counted. All values must be between [0,1].
 
-        Note that using the sine function is not common convention for
-        calculating the two point correlation function (see notes).
-
     sample2 : array_like, optional
         Npts2 x 3 array containing 3-D positions of points.
         Passing ``sample2`` as an input permits the calculation of
@@ -214,7 +211,7 @@ def s_mu_tpcf(sample1, s_bins, mu_bins, sample2=None, randoms=None,
         do_auto, do_cross, estimator, num_threads, max_sample_size,
         approx_cell1_size, approx_cell2_size, approx_cellran_size, seed)
 
-    sample1, s_bins, mu_prime_bins, sample2, randoms, period, do_auto, do_cross, num_threads,\
+    sample1, s_bins, mu_bins, sample2, randoms, period, do_auto, do_cross, num_threads,\
         _sample1_is_sample2, PBCs = _s_mu_tpcf_process_args(*function_args)
 
     # what needs to be done?
@@ -230,11 +227,11 @@ def s_mu_tpcf(sample1, s_bins, mu_bins, sample2=None, randoms=None,
         # this is arbitrarily set, but must remain consistent!
         NR = N1
 
-    D1D1, D1D2, D2D2 = pair_counts(sample1, sample2, s_bins, mu_prime_bins, period,
+    D1D1, D1D2, D2D2 = pair_counts(sample1, sample2, s_bins, mu_bins, period,
         num_threads, do_auto, do_cross, _sample1_is_sample2,
         approx_cell1_size, approx_cell2_size)
 
-    D1R, D2R, RR = random_counts(sample1, sample2, randoms, s_bins, mu_prime_bins,
+    D1R, D2R, RR = random_counts(sample1, sample2, randoms, s_bins, mu_bins,
         period, PBCs, num_threads, do_RR, do_DR, _sample1_is_sample2,
         approx_cell1_size, approx_cell2_size, approx_cellran_size)
 
@@ -258,7 +255,7 @@ def s_mu_tpcf(sample1, s_bins, mu_bins, sample2=None, randoms=None,
             return xi_11, xi_22
 
 
-def spherical_sector_volume(s, mu_prime):
+def spherical_sector_volume(s, mu):
     """
     This function is used to calculate analytical randoms.
 
@@ -267,12 +264,12 @@ def spherical_sector_volume(s, mu_prime):
 
     Note that the extra factor of 2 is to get the reflection.
     """
-    theta = np.arcsin(mu_prime)
+    theta = np.arccos(mu)
     vol = (2.0*np.pi/3.0) * np.outer((s**3.0), (1.0-np.cos(theta)))*2.0
     return vol
 
 
-def random_counts(sample1, sample2, randoms, s_bins, mu_prime_bins,
+def random_counts(sample1, sample2, randoms, s_bins, mu_bins,
         period, PBCs, num_threads, do_RR, do_DR, _sample1_is_sample2,
         approx_cell1_size, approx_cell2_size, approx_cellran_size):
     r"""
@@ -290,7 +287,7 @@ def random_counts(sample1, sample2, randoms, s_bins, mu_prime_bins,
     # PBCs and randoms.
     if randoms is not None:
         if do_RR is True:
-            RR = npairs_s_mu(randoms, randoms, s_bins, mu_prime_bins, period=period,
+            RR = npairs_s_mu(randoms, randoms, s_bins, mu_bins, period=period,
                              num_threads=num_threads,
                              approx_cell1_size=approx_cellran_size,
                              approx_cell2_size=approx_cellran_size)
@@ -298,7 +295,7 @@ def random_counts(sample1, sample2, randoms, s_bins, mu_prime_bins,
         else:
             RR = None
         if do_DR is True:
-            D1R = npairs_s_mu(sample1, randoms, s_bins, mu_prime_bins, period=period,
+            D1R = npairs_s_mu(sample1, randoms, s_bins, mu_bins, period=period,
                               num_threads=num_threads,
                               approx_cell1_size=approx_cell1_size,
                               approx_cell2_size=approx_cellran_size)
@@ -309,7 +306,7 @@ def random_counts(sample1, sample2, randoms, s_bins, mu_prime_bins,
             D2R = None
         else:
             if do_DR is True:
-                D2R = npairs_s_mu(sample2, randoms, s_bins, mu_prime_bins, period=period,
+                D2R = npairs_s_mu(sample2, randoms, s_bins, mu_bins, period=period,
                                   num_threads=num_threads,
                                   approx_cell1_size=approx_cell2_size,
                                   approx_cell2_size=approx_cellran_size)
@@ -325,7 +322,7 @@ def random_counts(sample1, sample2, randoms, s_bins, mu_prime_bins,
         NR = len(sample1)
 
         # do volume calculations
-        dv = spherical_sector_volume(s_bins, mu_prime_bins)
+        dv = spherical_sector_volume(s_bins, mu_bins)
         dv = np.diff(dv, axis=1)  # volume of wedges
         dv = np.diff(dv, axis=0)  # volume of wedge 'pieces'
         global_volume = period.prod()
@@ -348,14 +345,14 @@ def random_counts(sample1, sample2, randoms, s_bins, mu_prime_bins,
         raise ValueError('Un-supported combination of PBCs and randoms provided.')
 
 
-def pair_counts(sample1, sample2, s_bins, mu_prime_bins, period,
+def pair_counts(sample1, sample2, s_bins, mu_bins, period,
         num_threads, do_auto, do_cross, _sample1_is_sample2,
         approx_cell1_size, approx_cell2_size):
     """
     Count data pairs.
     """
     if do_auto is True:
-        D1D1 = npairs_s_mu(sample1, sample1, s_bins, mu_prime_bins, period=period,
+        D1D1 = npairs_s_mu(sample1, sample1, s_bins, mu_bins, period=period,
             num_threads=num_threads,
             approx_cell1_size=approx_cell1_size,
             approx_cell2_size=approx_cell1_size)
@@ -369,7 +366,7 @@ def pair_counts(sample1, sample2, s_bins, mu_prime_bins, period,
         D2D2 = D1D1
     else:
         if do_cross is True:
-            D1D2 = npairs_s_mu(sample1, sample2, s_bins, mu_prime_bins,
+            D1D2 = npairs_s_mu(sample1, sample2, s_bins, mu_bins,
                 period=period, num_threads=num_threads,
                 approx_cell1_size=approx_cell1_size,
                 approx_cell2_size=approx_cell2_size)
@@ -377,7 +374,7 @@ def pair_counts(sample1, sample2, s_bins, mu_prime_bins, period,
         else:
             D1D2 = None
         if do_auto is True:
-            D2D2 = npairs_s_mu(sample2, sample2, s_bins, mu_prime_bins, period=period,
+            D2D2 = npairs_s_mu(sample2, sample2, s_bins, mu_bins, period=period,
                 num_threads=num_threads,
                 approx_cell1_size=approx_cell2_size,
                 approx_cell2_size=approx_cell2_size)
@@ -414,12 +411,6 @@ def _s_mu_tpcf_process_args(sample1, s_bins, mu_bins, sample2, randoms,
     # process angular bins
     mu_bins = get_line_of_sight_bins_array(mu_bins)
 
-    # work with the sine of the angle between s and the LOS.  Only using cosine as the
-    # input because of convention.  sin(theta_los) increases as theta_los increases, which
-    # is required in order to get the pair counter to work.
-    theta = np.arccos(mu_bins)
-    mu_prime_bins = np.sin(theta)[::-1]  # must be increasing, remember to reverse result.
-
     if (np.min(mu_bins) < 0.0) | (np.max(mu_bins) > 1.0):
         msg = "`mu_bins` must be in the range [0,1]."
         raise ValueError(msg)
@@ -443,5 +434,5 @@ def _s_mu_tpcf_process_args(sample1, s_bins, mu_bins, sample2, randoms,
 
     verify_tpcf_estimator(estimator)
 
-    return sample1, s_bins, mu_prime_bins, sample2, randoms, period,\
+    return sample1, s_bins, mu_bins, sample2, randoms, period,\
         do_auto, do_cross, num_threads, _sample1_is_sample2, PBCs

From 0d8e88899f9f6ac020dd3c583b634ccc113d53e3 Mon Sep 17 00:00:00 2001
From: Duncan Campbell <duncan.alan.campbell@gmail.com>
Date: Fri, 9 Jun 2017 16:53:01 -0400
Subject: [PATCH 09/15] log API change npairs_s_mu

---
 CHANGES.rst | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/CHANGES.rst b/CHANGES.rst
index a060bc284..fea78431d 100644
--- a/CHANGES.rst
+++ b/CHANGES.rst
@@ -1,7 +1,7 @@
 0.6 (unreleased)
 ----------------
 
-- No changes yet
+- Changed the API for mock_observables.pair_counters.n_pairs_s_mu which now requires ``mu_bins`` to be in the conventional mu=cos(theta_LOS) format instead of mu=sin(theta_LOS).
 
 
 0.5 (2017-05-31)

From 86c93c8af6d3aa6364ea62b23b54a9a269e0cde7 Mon Sep 17 00:00:00 2001
From: Duncan Campbell <duncan.alan.campbell@gmail.com>
Date: Mon, 12 Jun 2017 13:34:50 -0400
Subject: [PATCH 10/15] fixed analytic randoms to work with new mu bin
 convention

---
 halotools/mock_observables/two_point_clustering/s_mu_tpcf.py | 4 +++-
 1 file changed, 3 insertions(+), 1 deletion(-)

diff --git a/halotools/mock_observables/two_point_clustering/s_mu_tpcf.py b/halotools/mock_observables/two_point_clustering/s_mu_tpcf.py
index f7f9748a6..f8ccb47d1 100644
--- a/halotools/mock_observables/two_point_clustering/s_mu_tpcf.py
+++ b/halotools/mock_observables/two_point_clustering/s_mu_tpcf.py
@@ -265,6 +265,7 @@ def spherical_sector_volume(s, mu):
     Note that the extra factor of 2 is to get the reflection.
     """
     theta = np.arccos(mu)
+    
     vol = (2.0*np.pi/3.0) * np.outer((s**3.0), (1.0-np.cos(theta)))*2.0
     return vol
 
@@ -322,7 +323,8 @@ def random_counts(sample1, sample2, randoms, s_bins, mu_bins,
         NR = len(sample1)
 
         # do volume calculations
-        dv = spherical_sector_volume(s_bins, mu_bins)
+        mu_bins_reverse_sorted = np.sort(mu_bins)[::-1]
+        dv = spherical_sector_volume(s_bins, mu_bins_reverse_sorted)
         dv = np.diff(dv, axis=1)  # volume of wedges
         dv = np.diff(dv, axis=0)  # volume of wedge 'pieces'
         global_volume = period.prod()

From d274e2afef22168a673b83cc0c67f98cefcdb98d Mon Sep 17 00:00:00 2001
From: Andrew Hearin <andrew.hearin@yale.edu>
Date: Mon, 12 Jun 2017 11:21:52 -0700
Subject: [PATCH 11/15] Added a couple of forgotten declarations

---
 .../pair_counters/cpairs/npairs_s_mu_engine.pyx        | 10 ++++++++--
 1 file changed, 8 insertions(+), 2 deletions(-)

diff --git a/halotools/mock_observables/pair_counters/cpairs/npairs_s_mu_engine.pyx b/halotools/mock_observables/pair_counters/cpairs/npairs_s_mu_engine.pyx
index 646c769ce..8c73ddd19 100644
--- a/halotools/mock_observables/pair_counters/cpairs/npairs_s_mu_engine.pyx
+++ b/halotools/mock_observables/pair_counters/cpairs/npairs_s_mu_engine.pyx
@@ -106,7 +106,9 @@ def npairs_s_mu_engine(double_mesh, x1in, y1in, z1in, x2in, y2in, z2in,
     cdef cnp.float64_t x2shift, y2shift, z2shift, dx, dy, dz, dxy_sq, dz_sq
     cdef cnp.float64_t x1tmp, y1tmp, z1tmp, s, mu
     cdef int Ni, Nj, i, j, k, l, g, max_k
-    cdef cnp.float64_t sqr_s_max = np.max(sqr_s_bins), sqr_mu_max = np.max(sqr_mu_bins)
+    cdef cnp.float64_t sqr_s_max = np.max(sqr_s_bins)
+    cdef cnp.float64_t sqr_mu_max = np.max(sqr_mu_bins)
+    cdef cnp.float64_t sqr_s, sqr_mu
 
     cdef cnp.float64_t[:] x_icell1, x_icell2
     cdef cnp.float64_t[:] y_icell1, y_icell2
@@ -194,7 +196,11 @@ def npairs_s_mu_engine(double_mesh, x1in, y1in, z1in, x2in, y2in, z2in,
                                     if sqr_s > sqr_s_max:
                                         continue
 
-                                    sqr_mu = (dxy_sq/sqr_s if sqr_s > 0.0 else 0.0)
+                                    if sqr_s > 0.0:
+                                        sqr_mu = dxy_sq/sqr_s
+                                    else:
+                                        sqr_mu = 0.0
+
                                     if sqr_mu > sqr_mu_max:
                                         continue
 

From 9adaa55a612ca14197589e69ad3c2f8b4b57d215 Mon Sep 17 00:00:00 2001
From: Duncan Campbell <duncan.alan.campbell@gmail.com>
Date: Mon, 12 Jun 2017 16:10:16 -0400
Subject: [PATCH 12/15] edited doc string

---
 .../pair_counters/npairs_s_mu.py              | 50 ++++++++++++-------
 1 file changed, 33 insertions(+), 17 deletions(-)

diff --git a/halotools/mock_observables/pair_counters/npairs_s_mu.py b/halotools/mock_observables/pair_counters/npairs_s_mu.py
index 685c6a9fd..d461ace43 100644
--- a/halotools/mock_observables/pair_counters/npairs_s_mu.py
+++ b/halotools/mock_observables/pair_counters/npairs_s_mu.py
@@ -21,19 +21,20 @@ def npairs_s_mu(sample1, sample2, s_bins, mu_bins, period=None,
         verbose=False, num_threads=1, approx_cell1_size=None, approx_cell2_size=None):
     r"""
     Function counts the number of pairs of points separated by less than
-    radial separation, *s,* and :math:`\mu\equiv\cos(\theta_{\rm los})`,
-    where :math:`\theta_{\rm los}` is the line-of-sight angle
-    between points and :math:`s^2 = r_{\rm parallel}^2 + r_{\rm perp}^2`.
-
-    Note that if sample1 == sample2 that the
-    `~halotools.mock_observables.npairs_s_mu` function double-counts pairs.
-    If your science application requires sample1==sample2 inputs and also pairs
-    to not be double-counted, simply divide the final counts by 2.
+    radial separation, :math:`s`, given by ``s_bins`` and 
+    angular distance, :math:`\mu\equiv\cos(\theta_{\rm los})`, given by ``mu_bins``,
+    where :math:`\theta_{\rm los}` is the angle between :math:`\vec{s}` and 
+    the line-of-sight (LOS).
+    
+    The first two dimensions (x, y) define the plane for perpendicular distances.
+    The third dimension (z) defines the LOS.  i.e. x,y positions are on
+    the plane of the sky, and z is the radial distance coordinate.  This is the 'distant
+    observer' approximation.
 
     A common variation of pair-counting calculations is to count pairs with
-    separations *between* two different distances *r1* and *r2*. You can retrieve
-    this information from the `~halotools.mock_observables.npairs_s_mu`
-    by taking `numpy.diff` of the returned array.
+    separations *between* two different distances, e.g. [*s1*,*s2*] and [*mu1*,*mu2*]. 
+    You can retrieve this information from `~halotools.mock_observables.npairs_s_mu`
+    by taking `numpy.diff` of the returned array along each axis.
 
     See Notes section for further clarification.
 
@@ -60,9 +61,6 @@ def npairs_s_mu(sample1, sample2, s_bins, mu_bins, period=None,
         :math:`\cos(\theta_{\rm LOS})` boundaries defining the bins in
         which pairs are counted. All values must be between [0,1].
 
-        Note that using the sine function is not common convention for
-        calculating the two point correlation function (see notes).
-
     period : array_like, optional
         Length-3 sequence defining the periodic boundary conditions
         in each dimension. If you instead provide a single scalar, Lbox,
@@ -100,9 +98,27 @@ def npairs_s_mu(sample1, sample2, s_bins, mu_bins, period=None,
 
     Notes
     -----
-    Along the first dimension of ``num_pairs``, :math:`s` (the radial seperation) increases.
-    Along the second dimension,  :math:`\mu` (the cosine of :math:`\theta_{\rm LOS}`) 
-    decreases, i.e. :math:`\theta_{\rm LOS}` increases.
+    Let :math:`\vec{s}` be the radial vector connnecting two points.
+    The magnitude, :math:`s`, is:
+
+    .. math::
+        s = \sqrt{r_{\parallel}^2+r_{\perp}^2},
+
+    where :math:`r_{\parallel}` is the separation parallel to the LOS
+    and :math:`r_{\perp}` is the separation perpednicular to the LOS.  :math:`\mu` is
+    the cosine of the angle, :math:`\theta_{\rm LOS}`, between the LOS
+    and :math:`\vec{s}`:
+
+    .. math::
+        \mu = \cos(\theta_{\rm LOS}) \equiv r_{\parallel}/s.
+    
+    Along the first dimension of ``num_pairs``, :math:`s` increases.
+    Along the second dimension,  :math:`\mu` decreases, 
+    i.e. :math:`\theta_{\rm LOS}` increases.
+    
+    If sample1 == sample2 that the `~halotools.mock_observables.npairs_s_mu` function 
+    double-counts pairs. If your science application requires sample1==sample2 inputs 
+    and also pairs to not be double-counted, simply divide the final counts by 2.
     
     One final point of clarification concerning double-counting may be in order.
     Suppose sample1==sample2 and s_bins[0]==0. Then the returned value for this bin

From f6ef4fe7fadc31678c8b0697fdf83cb5e5ffdc58 Mon Sep 17 00:00:00 2001
From: Duncan Campbell <duncan.alan.campbell@gmail.com>
Date: Mon, 12 Jun 2017 16:14:50 -0400
Subject: [PATCH 13/15] edited doc string

---
 halotools/mock_observables/pair_counters/npairs_s_mu.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/halotools/mock_observables/pair_counters/npairs_s_mu.py b/halotools/mock_observables/pair_counters/npairs_s_mu.py
index d461ace43..78e241ddf 100644
--- a/halotools/mock_observables/pair_counters/npairs_s_mu.py
+++ b/halotools/mock_observables/pair_counters/npairs_s_mu.py
@@ -32,7 +32,7 @@ def npairs_s_mu(sample1, sample2, s_bins, mu_bins, period=None,
     observer' approximation.
 
     A common variation of pair-counting calculations is to count pairs with
-    separations *between* two different distances, e.g. [*s1*,*s2*] and [*mu1*,*mu2*]. 
+    separations *between* two different distances, e.g. [s1 ,s2] and [mu1, mu2]. 
     You can retrieve this information from `~halotools.mock_observables.npairs_s_mu`
     by taking `numpy.diff` of the returned array along each axis.
 

From db2dcd2880c3934f67926a0da090ef602f298a90 Mon Sep 17 00:00:00 2001
From: Duncan Campbell <duncan.alan.campbell@gmail.com>
Date: Mon, 12 Jun 2017 17:06:03 -0400
Subject: [PATCH 14/15] fixed typos in doc strings and comments

---
 .../cpairs/npairs_s_mu_engine.pyx             | 23 +++++++++++--------
 .../pair_counters/npairs_s_mu.py              |  4 ++--
 .../mock_observables/pair_counters/pairs.py   | 12 +++++-----
 .../test_pair_counters/test_npairs_s_mu.py    |  8 +++----
 .../two_point_clustering/s_mu_tpcf.py         |  4 ++--
 5 files changed, 28 insertions(+), 23 deletions(-)

diff --git a/halotools/mock_observables/pair_counters/cpairs/npairs_s_mu_engine.pyx b/halotools/mock_observables/pair_counters/cpairs/npairs_s_mu_engine.pyx
index 8c73ddd19..1d2ae32e6 100644
--- a/halotools/mock_observables/pair_counters/cpairs/npairs_s_mu_engine.pyx
+++ b/halotools/mock_observables/pair_counters/cpairs/npairs_s_mu_engine.pyx
@@ -16,7 +16,8 @@ __all__ = ('npairs_s_mu_engine', )
 @cython.nonecheck(False)
 def npairs_s_mu_engine(double_mesh, x1in, y1in, z1in, x2in, y2in, z2in,
     s_bins_in, mu_bins_in, cell1_tuple):
-    r""" Cython engine for counting pairs of points as a function of projected separation.
+    r""" Cython engine for counting pairs of points as a function of radial separation, s,
+    and the angle between the line-of-sight (LOS) and s.
 
     Parameters
     ------------
@@ -48,7 +49,11 @@ def npairs_s_mu_engine(double_mesh, x1in, y1in, z1in, x2in, y2in, z2in,
     counts : array
         Integer array of length len(s_bins) giving the number of pairs
         separated by a distance less than the corresponding entry of ``s_bins``.
-
+    
+    Notes
+    -----
+    mu is defined as the sin(theta_LOS) so that as theta_LOS increases, mu increases.
+    
     """
     cdef cnp.float64_t[:] sqr_s_bins = s_bins_in * s_bins_in
     cdef cnp.float64_t[:] sqr_mu_bins = mu_bins_in * mu_bins_in
@@ -175,22 +180,22 @@ def npairs_s_mu_engine(double_mesh, x1in, y1in, z1in, x2in, y2in, z2in,
                         z_icell2 = z2[ifirst2:ilast2]
 
                         Nj = ilast2 - ifirst2
-                        #loop over points in cell1 points
+                        # loop over points in cell1 points
                         if Nj > 0:
                             for i in range(0,Ni):
                                 x1tmp = x_icell1[i] - x2shift
                                 y1tmp = y_icell1[i] - y2shift
                                 z1tmp = z_icell1[i] - z2shift
-                                #loop over points in cell2 points
+                                # loop over points in cell2 points
                                 for j in range(0,Nj):
-                                    #calculate the square distance
+                                    # calculate the square distance
                                     dx = x1tmp - x_icell2[j]
                                     dy = y1tmp - y_icell2[j]
                                     dz = z1tmp - z_icell2[j]
                                     dxy_sq = dx*dx + dy*dy
                                     dz_sq = dz*dz
 
-                                    #transform to s and mu
+                                    # transform to s and mu
                                     sqr_s = dz_sq + dxy_sq
 
                                     if sqr_s > sqr_s_max:
@@ -204,9 +209,9 @@ def npairs_s_mu_engine(double_mesh, x1in, y1in, z1in, x2in, y2in, z2in,
                                     if sqr_mu > sqr_mu_max:
                                         continue
 
-                                    # MS: I have intentionally split
-                                    # up the loop. Since division
-                                    # is slow, computing mu is a bottle-neck.
+                                    # The loop has been intentionally split up
+                                    # Since division is slow,
+                                    # computing mu is a bottle-neck.
                                     # Computing the 's' bin however can proceed
                                     # in the meantime.
                                     k = num_s_bins-2
diff --git a/halotools/mock_observables/pair_counters/npairs_s_mu.py b/halotools/mock_observables/pair_counters/npairs_s_mu.py
index 78e241ddf..4ab64e153 100644
--- a/halotools/mock_observables/pair_counters/npairs_s_mu.py
+++ b/halotools/mock_observables/pair_counters/npairs_s_mu.py
@@ -192,11 +192,11 @@ def npairs_s_mu(sample1, sample2, s_bins, mu_bins, period=None,
         approx_x2cell_size, approx_y2cell_size, approx_z2cell_size,
         search_xlength, search_ylength, search_zlength, xperiod, yperiod, zperiod, PBCs)
 
-    # # Create a function object that has a single argument, for parallelization purposes
+    # Create a function object that has a single argument, for parallelization purposes
     engine = partial(npairs_s_mu_engine,
         double_mesh, x1in, y1in, z1in, x2in, y2in, z2in, s_bins, mu_bins_prime)
 
-    # # Calculate the cell1 indices that will be looped over by the engine
+    # Calculate the cell1 indices that will be looped over by the engine
     num_threads, cell1_tuples = _cell1_parallelization_indices(
         double_mesh.mesh1.ncells, num_threads)
 
diff --git a/halotools/mock_observables/pair_counters/pairs.py b/halotools/mock_observables/pair_counters/pairs.py
index 544add5c4..c2ac53bbb 100755
--- a/halotools/mock_observables/pair_counters/pairs.py
+++ b/halotools/mock_observables/pair_counters/pairs.py
@@ -7,7 +7,7 @@
 from __future__ import absolute_import, division, print_function, unicode_literals
 import numpy as np
 
-__all__ = ['npairs', 'wnpairs', 'xy_z_npairs', 'xy_z_wnpairs']
+__all__ = ['npairs', 'wnpairs', 'xy_z_npairs', 'xy_z_wnpairs', 's_mu_npairs']
 __author__ = ['Duncan Campbell']
 
 
@@ -355,7 +355,7 @@ def xy_z_wnpairs(sample1, sample2, rp_bins, pi_bins, period=None, weights1=None,
 def s_mu_npairs(sample1, sample2, s_bins, mu_bins, period=None):
     """
     Calculate the number of pairs with 3D radial separations less than or equal to
-    :math:`s`, and angular seperations along the LOS, :math:`\mu=\cos(\theta_{\rm LOS})`.
+    :math:`s`, and angular separations along the LOS, :math:`\mu=\cos(\theta_{\rm LOS})`.
 
     Assumes the first N-1 dimensions are perpendicular to the line-of-sight (LOS), and
     the final dimension is parallel to the LOS.
@@ -387,11 +387,11 @@ def s_mu_npairs(sample1, sample2, s_bins, mu_bins, period=None):
     Returns
     -------
     N_pairs : ndarray of shape (num_s_bin_edges, num_mu_bin_edges) storing the 
-        number counts of pairs with seperations less than ``s_bins`` and ``mu_bins``
+        number counts of pairs with separations less than ``s_bins`` and ``mu_bins``
     
     Notes
     -----
-    Along the first dimension of ``N_pairs``, :math:`s` (the radial seperation) increases.
+    Along the first dimension of ``N_pairs``, :math:`s` (the radial separation) increases.
     Along the second dimension,  :math:`\mu` (the cosine of :math:`\theta_{\rm LOS}`) 
     decreases, i.e. :math:`\theta_{\rm LOS}` increases.
     """
@@ -570,7 +570,7 @@ def perpendicular_distance(x1, x2, period=None):
 
 def theta_LOS(x1, x2, period=None):
     """
-    Find the seperation angle from the LOS between x1 & x2, accounting for box periodicity.
+    Find the separation angle from the LOS between x1 & x2, accounting for box periodicity.
 
     Assumes the first N-1 dimensions are perpendicular to the line-of-sight (LOS).
 
@@ -612,7 +612,7 @@ def theta_LOS(x1, x2, period=None):
     r_perp = perpendicular_distance(x1, x2, period=period)
     r_parallel = parallel_distance(x1, x2, period=period)
     
-    # deal with zero seperation
+    # deal with zero separation
     r = np.sqrt(r_perp**2 + r_parallel**2)
     mask = (r>0.0)
     
diff --git a/halotools/mock_observables/pair_counters/test_pair_counters/test_npairs_s_mu.py b/halotools/mock_observables/pair_counters/test_pair_counters/test_npairs_s_mu.py
index 7bf1e694e..a0524f8f5 100755
--- a/halotools/mock_observables/pair_counters/test_pair_counters/test_npairs_s_mu.py
+++ b/halotools/mock_observables/pair_counters/test_pair_counters/test_npairs_s_mu.py
@@ -55,7 +55,7 @@
 
 def test_npairs_s_mu_periodic():
     """
-    test npairs_s_mu with periodic boundary conditions.
+    test npairs_s_mu with periodic boundary conditions on a random set of points.
     """
     
     # define bins
@@ -95,7 +95,7 @@ def test_npairs_s_mu_periodic():
 
 def test_npairs_s_mu_nonperiodic():
     """
-    test npairs_s_mu without periodic boundary conditions.
+    test npairs_s_mu without periodic boundary conditions on a random set of points.
     """
 
     # define bins
@@ -136,7 +136,7 @@ def test_npairs_s_mu_nonperiodic():
 def test_npairs_s_mu_point_surrounded_by_circle():
     """
     test npairs_s_mu without periodic boundary conditions on a point surrounded by an
-    evenly distribued circle of points perpdindicular to the LOS.
+    evenly distributed circle of points aligned with the LOS.
     """
     
     # put one point per degree in the circle
@@ -160,7 +160,7 @@ def test_npairs_s_mu_point_surrounded_by_circle():
     pairs = np.diff(np.diff(result, axis=0),axis=1)
     
     # note that there should be 4 pairs per mu bin
-    # since each quadrant of te circle counts once for each mu bin
+    # since each quadrant of the circle counts once for each mu bin
     # and there is one point per degree
     
     msg = "The number of pairs between the origin and a circle with one point per degree is incorrect."
diff --git a/halotools/mock_observables/two_point_clustering/s_mu_tpcf.py b/halotools/mock_observables/two_point_clustering/s_mu_tpcf.py
index f8ccb47d1..5ed3ea6aa 100644
--- a/halotools/mock_observables/two_point_clustering/s_mu_tpcf.py
+++ b/halotools/mock_observables/two_point_clustering/s_mu_tpcf.py
@@ -235,8 +235,8 @@ def s_mu_tpcf(sample1, s_bins, mu_bins, sample2=None, randoms=None,
         period, PBCs, num_threads, do_RR, do_DR, _sample1_is_sample2,
         approx_cell1_size, approx_cell2_size, approx_cellran_size)
 
-    # return results.  remember to reverse the final result because we used
-    # mu_prime = sin(theta_los) bins instead of the user passed in mu = cos(theta_los).
+    # return results.  remember to reverse the final result since
+    # the pair counts are done in order of increasing theta_LOS (i.e. decreasing mu)
     if _sample1_is_sample2:
         xi_11 = _TP_estimator(D1D1, D1R, RR, N1, N1, NR, NR, estimator)[:, ::-1]
         return xi_11

From fefd22f5c1fb2509cf63230015b2d13172d1ae0f Mon Sep 17 00:00:00 2001
From: Duncan Campbell <duncan.alan.campbell@gmail.com>
Date: Mon, 12 Jun 2017 18:58:02 -0400
Subject: [PATCH 15/15] fixed issue with doc strings that affects python 3
 builds

---
 .../mock_observables/pair_counters/pairs.py   | 26 +++++++------------
 .../test_pair_counters/test_npairs_s_mu.py    |  9 ++++---
 2 files changed, 15 insertions(+), 20 deletions(-)

diff --git a/halotools/mock_observables/pair_counters/pairs.py b/halotools/mock_observables/pair_counters/pairs.py
index c2ac53bbb..034deddd7 100755
--- a/halotools/mock_observables/pair_counters/pairs.py
+++ b/halotools/mock_observables/pair_counters/pairs.py
@@ -1,4 +1,4 @@
-"""
+r"""
 simple python brute force pair counting functions.  The primary purpose of these functions
 is as a sanity check on more complex pair counting techniques.  These functions should not
 be used on large data sets, as memory usage is very large, and runtimes can be very slow.
@@ -12,7 +12,7 @@
 
 
 def npairs(sample1, sample2, rbins, period=None):
-    """
+    r"""
     Calculate the number of pairs with separations less than or equal to rbins[i].
 
     Parameters
@@ -38,7 +38,6 @@ def npairs(sample1, sample2, rbins, period=None):
     -------
     N_pairs : array of length len(rbins)
         number counts of pairs
-
     """
 
     sample1 = np.atleast_2d(sample1)
@@ -83,7 +82,7 @@ def npairs(sample1, sample2, rbins, period=None):
 
 
 def xy_z_npairs(sample1, sample2, rp_bins, pi_bins, period=None):
-    """
+    r"""
     Calculate the number of pairs with parellal separations less than or equal to
     pi_bins[i], and perpendicular separations less than or equal to rp_bins[i].
 
@@ -116,7 +115,6 @@ def xy_z_npairs(sample1, sample2, rp_bins, pi_bins, period=None):
     -------
     N_pairs : ndarray of shape (len(rp_bins),len(pi_bins))
         number counts of pairs
-
     """
 
     sample1 = np.atleast_2d(sample1)
@@ -159,7 +157,7 @@ def xy_z_npairs(sample1, sample2, rp_bins, pi_bins, period=None):
 
 
 def wnpairs(sample1, sample2, r, period=None, weights1=None, weights2=None):
-    """
+    r"""
     Calculate the weighted number of pairs with separations less than or equal to rbins[i].
 
     Parameters
@@ -191,7 +189,6 @@ def wnpairs(sample1, sample2, r, period=None, weights1=None, weights2=None):
     -------
     wN_pairs : array of length len(rbins)
         weighted number counts of pairs
-
     """
 
     sample1 = np.atleast_2d(sample1)
@@ -251,7 +248,7 @@ def wnpairs(sample1, sample2, r, period=None, weights1=None, weights2=None):
 
 
 def xy_z_wnpairs(sample1, sample2, rp_bins, pi_bins, period=None, weights1=None, weights2=None):
-    """
+    r"""
     Calculate the number of weighted pairs with parellal separations less than or equal to
     pi_bins[i], and perpendicular separations less than or equal to rp_bins[i].
 
@@ -291,7 +288,6 @@ def xy_z_wnpairs(sample1, sample2, rp_bins, pi_bins, period=None, weights1=None,
     -------
     wN_pairs : ndarray of shape (len(rp_bins),len(pi_bins))
         weighted number counts of pairs
-
     """
 
     sample1 = np.atleast_2d(sample1)
@@ -353,7 +349,7 @@ def xy_z_wnpairs(sample1, sample2, rp_bins, pi_bins, period=None, weights1=None,
 
 
 def s_mu_npairs(sample1, sample2, s_bins, mu_bins, period=None):
-    """
+    r"""
     Calculate the number of pairs with 3D radial separations less than or equal to
     :math:`s`, and angular separations along the LOS, :math:`\mu=\cos(\theta_{\rm LOS})`.
 
@@ -443,7 +439,7 @@ def s_mu_npairs(sample1, sample2, s_bins, mu_bins, period=None):
 
 
 def distance(x1, x2, period=None):
-    """
+    r"""
     Find the Euclidean distance between x1 & x2, accounting for box periodicity.
 
     Parameters
@@ -483,7 +479,7 @@ def distance(x1, x2, period=None):
 
 
 def parallel_distance(x1, x2, period=None):
-    """
+    r"""
     Find the parallel distance between x1 & x2, accounting for box periodicity.
 
     Assumes the last dimension is the line-of-sight.
@@ -503,7 +499,6 @@ def parallel_distance(x1, x2, period=None):
     Returns
     -------
     distance : array
-
     """
 
     x1 = np.atleast_2d(x1)
@@ -526,7 +521,7 @@ def parallel_distance(x1, x2, period=None):
 
 
 def perpendicular_distance(x1, x2, period=None):
-    """
+    r"""
     Find the perpendicular distance between x1 & x2, accounting for box periodicity.
 
     Assumes the first N-1 dimensions are perpendicular to the line-of-sight.
@@ -546,7 +541,6 @@ def perpendicular_distance(x1, x2, period=None):
     Returns
     -------
     distance : array
-
     """
 
     x1 = np.atleast_2d(x1)
@@ -569,7 +563,7 @@ def perpendicular_distance(x1, x2, period=None):
 
 
 def theta_LOS(x1, x2, period=None):
-    """
+    r"""
     Find the separation angle from the LOS between x1 & x2, accounting for box periodicity.
 
     Assumes the first N-1 dimensions are perpendicular to the line-of-sight (LOS).
diff --git a/halotools/mock_observables/pair_counters/test_pair_counters/test_npairs_s_mu.py b/halotools/mock_observables/pair_counters/test_pair_counters/test_npairs_s_mu.py
index a0524f8f5..fc1db2753 100755
--- a/halotools/mock_observables/pair_counters/test_pair_counters/test_npairs_s_mu.py
+++ b/halotools/mock_observables/pair_counters/test_pair_counters/test_npairs_s_mu.py
@@ -1,4 +1,5 @@
-"""
+r"""
+test module for pair counts in s and mu bins
 """
 from __future__ import absolute_import, division, print_function, unicode_literals
 
@@ -54,7 +55,7 @@
 
 
 def test_npairs_s_mu_periodic():
-    """
+    r"""
     test npairs_s_mu with periodic boundary conditions on a random set of points.
     """
     
@@ -94,7 +95,7 @@ def test_npairs_s_mu_periodic():
 
 
 def test_npairs_s_mu_nonperiodic():
-    """
+    r"""
     test npairs_s_mu without periodic boundary conditions on a random set of points.
     """
 
@@ -134,7 +135,7 @@ def test_npairs_s_mu_nonperiodic():
 
 
 def test_npairs_s_mu_point_surrounded_by_circle():
-    """
+    r"""
     test npairs_s_mu without periodic boundary conditions on a point surrounded by an
     evenly distributed circle of points aligned with the LOS.
     """