Switch default 3pt to use LogSAS and ordered=True

CHANGELOG.rst

@@ -25,9 +25,15 @@ API Changes
   normally used directly by users, so it shouldn't be noticeable in user code. (#155)
 - Removed all deprecations from the 4.x series. (#156)
 - Removed support for reading back in output files from the 3.x series. (#165)
-- Removed the 3pt CrossCorrelation classes.  See the new ``ordered=True`` option to the 3pt
-  ``process`` method instead to get correlations where the order of the three catalogs is fixed.
-  This is simpler and more intuitive, and for many use cases it is more efficient. (#165)
+- Removed the 3pt CrossCorrelation classes, which used to be the way to get ordered three-point
+  correlations.  But they were rather unwieldy and not very intuitive.  The new ``ordered``
+  option to the three-point ``process`` methods is much simpler and more efficient for the common
+  case of only wanting a single order for the catalogs. (#165)
+- Switched the default behavior of 3pt cross-correlations to respect the order of the catalogs
+  in the triangle definitions.  That is, points from cat1 will be at P1 in the triangle,
+  points from cat2 at P2, and points from cat3 at P3.  To recover the old behavior, you may
+  use the new ``ordered=False`` option. (#166)
+- Switched the default binning for three-point correlations to LogSAS, rather than LogRUV. (#166)
 
 
 Performance improvements

include/BinType.h

@@ -417,9 +417,9 @@ struct BinTypeHelper<LogRUV>
         double minu, double minusq, double maxu, double maxusq,
         double minv, double minvsq, double maxv, double maxvsq)
     {
-        dbg<<"Stop111: "<<std::sqrt(d1sq)<<"  "<<std::sqrt(d2sq)<<"  "<<std::sqrt(d3sq)<<std::endl;
-        dbg<<"sizes = "<<s1<<"  "<<s2<<"  "<<s3<<std::endl;
-        dbg<<"sep range = "<<minsep<<"  "<<maxsep<<std::endl;
+        xdbg<<"Stop111: "<<std::sqrt(d1sq)<<"  "<<std::sqrt(d2sq)<<"  "<<std::sqrt(d3sq)<<std::endl;
+        xdbg<<"sizes = "<<s1<<"  "<<s2<<"  "<<s3<<std::endl;
+        xdbg<<"sep range = "<<minsep<<"  "<<maxsep<<std::endl;
 
         if (!O) {
             Assert(d1sq >= d2sq);
@@ -457,19 +457,19 @@ struct BinTypeHelper<LogRUV>
 
             // (d2 + s1+s3) < (d3 - s1-s2)
             if (O == 3 && d3sq > SQR(d2 + 2*s1+s2+s3)) {
-                dbg<<"d2 cannot be larger than d3\n";
+                xdbg<<"d2 cannot be larger than d3\n";
                 return true;
             }
             // (d1 + s2+s3) < (d2 - s1-s3)
             double ss = s1+s2+2*s3;
             if (ss < d2 && d1sq < SQR(d2 - ss)) {
-                dbg<<"d1 cannot be larger than d2\n";
+                xdbg<<"d1 cannot be larger than d2\n";
                 return true;
             }
             d1 = sqrt(d1sq);
             // (d1 + s2+s3) < (d3 - s1-s2)
             if (d3sq > SQR(d1 + s1+2*s2+s3)) {
-                dbg<<"d1 cannot be larger than d3\n";
+                xdbg<<"d1 cannot be larger than d3\n";
                 return true;
             }
         }
@@ -679,25 +679,27 @@ struct BinTypeHelper<LogRUV>
     static bool isTriangleInRange(const BaseCell<C>& c1, const BaseCell<C>& c2,
                                   const BaseCell<C>& c3,
                                   const MetricHelper<M,0>& metric,
+                                  double d1sq, double d2sq, double d3sq,
                                   double d1, double d2, double d3, double& u, double& v,
                                   double logminsep,
-                                  double minsep, double maxsep, double binsize, double nbins,
-                                  double minu, double maxu, double ubinsize, double nubins,
-                                  double minv, double maxv, double vbinsize, double nvbins,
+                                  double minsep, double maxsep, double binsize, int nbins,
+                                  double minu, double maxu, double ubinsize, int nubins,
+                                  double minv, double maxv, double vbinsize, int nvbins,
                                   double& logd1, double& logd2, double& logd3,
                                   int ntot, int& index)
     {
         // Make sure all the quantities we thought should be set have been.
         Assert(d1 > 0.);
         Assert(d3 > 0.);
         Assert(u > 0.);
-        Assert(v >= 0.);  // v can potentially == 0.
 
         if (O && !(d1 >= d2 && d2 >= d3)) {
             xdbg<<"Sides are not in correct size ordering d1 >= d2 >= d3\n";
             return false;
         }
 
+        Assert(v >= 0.);  // v can potentially == 0.
+
         if (d2 < minsep || d2 >= maxsep) {
             xdbg<<"d2 not in minsep .. maxsep\n";
             return false;
@@ -722,7 +724,7 @@ struct BinTypeHelper<LogRUV>
         Assert(kr >= 0);
         Assert(kr <= nbins);
         if (kr == nbins) --kr;  // This is rare, but can happen with numerical differences
-                                 // between the math for log and for non-log checks.
+                                // between the math for log and for non-log checks.
         Assert(kr < nbins);
 
         int ku = int(floor((u-minu)/ubinsize));
@@ -747,8 +749,7 @@ struct BinTypeHelper<LogRUV>
         Assert(kv < nvbins);
 
         // Now account for negative v
-        if (!metric.CCW(c1.getData().getPos(), c2.getData().getPos(),
-                        c3.getData().getPos())) {
+        if (!metric.CCW(c1.getPos(), c2.getPos(), c3.getPos())) {
             v = -v;
             kv = nvbins - kv - 1;
         } else {
@@ -783,7 +784,7 @@ struct BinTypeHelper<LogSAS>
 {
     enum { sort_d123 = false };
 
-    static int calculateNTot(int nbins, int nphibins, int nvbins)
+    static int calculateNTot(int nbins, int nphibins, int )
     { return nbins * nbins * nphibins; }
 
     static bool tooSmallS2(double s2, double halfminsep, double minphi, double )
@@ -900,7 +901,7 @@ struct BinTypeHelper<LogSAS>
 
         // If we are not swapping 2,3, stop if orientation cannot be counter-clockwise.
         if (O > 1 &&
-            !metric.CCW(c1.getData().getPos(), c3.getData().getPos(), c2.getData().getPos())) {
+            !metric.CCW(c1.getPos(), c3.getPos(), c2.getPos())) {
             // For skinny triangles, be careful that the points can't flip to the other side.
             // This is similar to the calculation below.  We effecively check that cosphi can't
             // increase to 1.
@@ -1149,11 +1150,12 @@ struct BinTypeHelper<LogSAS>
     static bool isTriangleInRange(const BaseCell<C>& c1, const BaseCell<C>& c2,
                                   const BaseCell<C>& c3,
                                   const MetricHelper<M,0>& metric,
+                                  double d1sq, double d2sq, double d3sq,
                                   double d1, double d2, double d3, double& phi, double& cosphi,
                                   double logminsep,
-                                  double minsep, double maxsep, double binsize, double nbins,
-                                  double minphi, double maxphi, double phibinsize, double nphibins,
-                                  double , double , double , double ,
+                                  double minsep, double maxsep, double binsize, int nbins,
+                                  double minphi, double maxphi, double phibinsize, int nphibins,
+                                  double , double , double , int ,
                                   double& logd1, double& logd2, double& logd3,
                                   int ntot, int& index)
     {
@@ -1181,12 +1183,12 @@ struct BinTypeHelper<LogSAS>
         }
 
         if (O > 1 &&
-            !metric.CCW(c1.getData().getPos(), c3.getData().getPos(), c2.getData().getPos())) {
+            !metric.CCW(c1.getPos(), c3.getPos(), c2.getPos())) {
             xdbg<<"Triangle is not CCW.\n";
             return false;
         }
 
-        XAssert(metric.CCW(c1.getData().getPos(), c3.getData().getPos(), c2.getData().getPos()));
+        XAssert(metric.CCW(c1.getPos(), c3.getPos(), c2.getPos()));
 
         if (phi < minphi || phi >= maxphi) {
             xdbg<<"phi not in minphi .. maxphi\n";
@@ -1223,7 +1225,7 @@ struct BinTypeHelper<LogSAS>
         Assert(kphi < nphibins);
 
         xdbg<<"d1,d2,d3,phi = "<<d1<<", "<<d2<<", "<<d3<<",  "<<phi<<std::endl;
-        index = (kr2 * nphibins + kphi) * nbins + kr3;
+        index = (kr2 * nbins + kr3) * nphibins + kphi;
         xdbg<<"kr2,kr3,kphi = "<<kr2<<", "<<kr3<<", "<<kphi<<":  "<<index<<std::endl;
         Assert(index >= 0);
         Assert(index < ntot);

include/Corr3.h

@@ -41,6 +41,7 @@ class BaseCorr3
     virtual ~BaseCorr3() {}
 
     virtual std::shared_ptr<BaseCorr3> duplicate() =0;
+    virtual void writeZeta(std::ostream& os, int n) const = 0;
 
     virtual void addData(const BaseCorr3& rhs) =0;
 
@@ -162,6 +163,9 @@ class Corr3 : public BaseCorr3
 
     void clear();  // Set all data to 0.
 
+    void writeZeta(std::ostream& os, int n=1) const
+    { _zeta.write_full(os, n); }
+
     template <int C>
     void finishProcess(
         const BaseCell<C>& c1, const BaseCell<C>& c2, const BaseCell<C>& c3,

include/Field.h

@@ -51,7 +51,14 @@ class BaseField
     Position<C> getCenter() const { return _center; }
     double getSize() const { return std::sqrt(_sizesq); }
     long getNTopLevel() const { BuildCells(); return long(_cells.size()); }
-    const std::vector<BaseCell<C>*>& getCells() const { BuildCells(); return _cells; }
+    const std::vector<const BaseCell<C>*>& getCells() const
+    {
+        BuildCells();
+        // const_cast is insufficient to turn this into a vector of const BaseCell*.
+        // cf. https://stackoverflow.com/questions/19122858/why-is-a-vector-of-pointers-not-castable-to-a-const-vector-of-const-pointers
+        // But reinterpret_cast here is safe.
+        return reinterpret_cast<std::vector<const BaseCell<C>*>&>(_cells);
+    }
     long countNear(double x, double y, double z, double sep) const;
     void getNear(double x, double y, double z, double sep, long* indices, long n) const;
 

include/Metric.h

@@ -712,9 +712,9 @@ struct MetricHelper<Arc, P>
         // Rather than do trig though, recompute the chord lengths so we can compute
         // cos(phi) with just regular arithmetic and a sqrt.
         //
-        double d1sq = (c2.getData().getPos() - c3.getData().getPos()).normSq();
-        double d2sq = (c1.getData().getPos() - c3.getData().getPos()).normSq();
-        double d3sq = (c1.getData().getPos() - c2.getData().getPos()).normSq();
+        double d1sq = (c2.getPos() - c3.getPos()).normSq();
+        double d2sq = (c1.getPos() - c3.getPos()).normSq();
+        double d3sq = (c1.getPos() - c2.getPos()).normSq();
 
         double cosphi = (d2sq + d3sq - 0.5*d2sq*d3sq - d1sq);
         cosphi /= 2. * std::sqrt( d2sq * d3sq * (1.-0.25*d2sq) * (1.-0.25*d3sq) );

include/ProjectHelper.h

@@ -72,7 +72,7 @@ struct ProjectHelper<Flat>
         const Cell<D1,Flat>& c1, const Cell<D,Flat>& c2, std::complex<double>& z2)
     {
         // Project given spin-s quantity to the line connecting them.
-        std::complex<double> r(c2.getData().getPos() - c1.getData().getPos());
+        std::complex<double> r(c2.getPos() - c1.getPos());
         z2 *= _expmsialpha<D>(r);
     }
 
@@ -82,7 +82,7 @@ struct ProjectHelper<Flat>
         std::complex<double>& z1, std::complex<double>& z2)
     {
         // Project given spin-s quantities to the line connecting them.
-        std::complex<double> r(c2.getData().getPos() - c1.getData().getPos());
+        std::complex<double> r(c2.getPos() - c1.getPos());
         std::complex<double> expmsialpha = _expmsialpha<D>(r);
         z1 *= expmsialpha;
         z2 *= expmsialpha;
@@ -94,9 +94,9 @@ struct ProjectHelper<Flat>
         std::complex<double>& z1, std::complex<double>& z2, std::complex<double>& z3)
     {
         // Project given spin-s quantities to the line connecting each to the centroid.
-        const Position<Flat>& p1 = c1.getData().getPos();
-        const Position<Flat>& p2 = c2.getData().getPos();
-        const Position<Flat>& p3 = c3.getData().getPos();
+        const Position<Flat>& p1 = c1.getPos();
+        const Position<Flat>& p2 = c2.getPos();
+        const Position<Flat>& p3 = c3.getPos();
         Position<Flat> cen = (p1 + p2 + p3)/3.;
         std::complex<double> r1(cen - p1);
         std::complex<double> r2(cen - p2);
@@ -189,8 +189,8 @@ struct ProjectHelper<Sphere>
     static void Project(
         const Cell<D1,Sphere>& c1, const Cell<D,Sphere>& c2, std::complex<double>& z2)
     {
-        const Position<Sphere>& p1 = c1.getData().getPos();
-        const Position<Sphere>& p2 = c2.getData().getPos();
+        const Position<Sphere>& p1 = c1.getPos();
+        const Position<Sphere>& p2 = c2.getPos();
         std::complex<double> r = calculate_direction(p1,p2);
         z2 *= _expmsialpha<D>(r);
     }
@@ -200,8 +200,8 @@ struct ProjectHelper<Sphere>
         const Cell<D,Sphere>& c1, const Cell<D,Sphere>& c2,
         std::complex<double>& z1, std::complex<double>& z2)
     {
-        const Position<Sphere>& p1 = c1.getData().getPos();
-        const Position<Sphere>& p2 = c2.getData().getPos();
+        const Position<Sphere>& p1 = c1.getPos();
+        const Position<Sphere>& p2 = c2.getPos();
         std::complex<double> r12 = calculate_direction(p1,p2);
         std::complex<double> expmsialpha = _expmsialpha<D>(r12);
         z2 *= expmsialpha;
@@ -220,9 +220,9 @@ struct ProjectHelper<Sphere>
         const Cell<D,Sphere>& c1, const Cell<D,Sphere>& c2, const Cell<D,Sphere>& c3,
         std::complex<double>& z1, std::complex<double>& z2, std::complex<double>& z3)
     {
-        const Position<Sphere>& p1 = c1.getData().getPos();
-        const Position<Sphere>& p2 = c2.getData().getPos();
-        const Position<Sphere>& p3 = c3.getData().getPos();
+        const Position<Sphere>& p1 = c1.getPos();
+        const Position<Sphere>& p2 = c2.getPos();
+        const Position<Sphere>& p3 = c3.getPos();
         Position<Sphere> cen((p1 + p2 + p3)/3.);
         z1 *= _expmsialpha<D>(calculate_direction(cen,p1));
         z2 *= _expmsialpha<D>(calculate_direction(cen,p2));
@@ -240,8 +240,8 @@ struct ProjectHelper<ThreeD>
     static void Project(
         const Cell<D1,ThreeD>& c1, const Cell<D,ThreeD>& c2, std::complex<double>& z2)
     {
-        const Position<ThreeD>& p1 = c1.getData().getPos();
-        const Position<ThreeD>& p2 = c2.getData().getPos();
+        const Position<ThreeD>& p1 = c1.getPos();
+        const Position<ThreeD>& p2 = c2.getPos();
         Position<Sphere> sp1(p1);
         Position<Sphere> sp2(p2);
         std::complex<double> r = ProjectHelper<Sphere>::calculate_direction(sp1,sp2);
@@ -253,8 +253,8 @@ struct ProjectHelper<ThreeD>
         const Cell<D,ThreeD>& c1, const Cell<D,ThreeD>& c2,
         std::complex<double>& z1, std::complex<double>& z2)
     {
-        const Position<ThreeD>& p1 = c1.getData().getPos();
-        const Position<ThreeD>& p2 = c2.getData().getPos();
+        const Position<ThreeD>& p1 = c1.getPos();
+        const Position<ThreeD>& p2 = c2.getPos();
         Position<Sphere> sp1(p1);
         Position<Sphere> sp2(p2);
         std::complex<double> r12 = ProjectHelper<Sphere>::calculate_direction(sp1,sp2);
@@ -271,9 +271,9 @@ struct ProjectHelper<ThreeD>
         const Cell<D,ThreeD>& c1, const Cell<D,ThreeD>& c2, const Cell<D,ThreeD>& c3,
         std::complex<double>& z1, std::complex<double>& z2, std::complex<double>& z3)
     {
-        const Position<ThreeD>& p1 = c1.getData().getPos();
-        const Position<ThreeD>& p2 = c2.getData().getPos();
-        const Position<ThreeD>& p3 = c3.getData().getPos();
+        const Position<ThreeD>& p1 = c1.getPos();
+        const Position<ThreeD>& p2 = c2.getPos();
+        const Position<ThreeD>& p3 = c3.getPos();
         Position<Sphere> sp1(p1);
         Position<Sphere> sp2(p2);
         Position<Sphere> sp3(p3);

include/Split.h

@@ -126,11 +126,11 @@ inline bool Check(
     // Checks that d1,d2,d3 are correct for the three Cells given.
     // Used as a debugging check.
     bool ok=true;
-    if (Dist(c3.getData().getPos(),c2.getData().getPos())-d1 > 0.0001)
+    if (Dist(c3.getPos(),c2.getPos())-d1 > 0.0001)
     { std::cerr<<"d1\n"; ok = false; }
-    if (Dist(c1.getData().getPos(),c3.getData().getPos())-d2 > 0.0001)
+    if (Dist(c1.getPos(),c3.getPos())-d2 > 0.0001)
     { std::cerr<<"d2\n"; ok = false; }
-    if (Dist(c2.getData().getPos(),c1.getData().getPos())-d3 > 0.0001)
+    if (Dist(c2.getPos(),c1.getPos())-d3 > 0.0001)
     { std::cerr<<"d3\n"; ok = false; }
     if (d1 > d2+d3+0.0001) { std::cerr<<"sum d1\n"; ok = false; }
     if (d2 > d1+d3+0.0001) { std::cerr<<"sum d2\n"; ok = false; }

src/Corr2.cpp

@@ -150,6 +150,8 @@ void BaseCorr2::process(const BaseField<C>& field, bool dots)
     dbg<<"field has "<<n1<<" top level nodes\n";
     Assert(n1 > 0);
 
+    const std::vector<const BaseCell<C>*>& cells = field.getCells();
+
 #ifdef _OPENMP
 #pragma omp parallel
     {
@@ -176,10 +178,10 @@ void BaseCorr2::process(const BaseField<C>& field, bool dots)
 #endif
                 if (dots) std::cout<<'.'<<std::flush;
             }
-            const BaseCell<C>& c1 = *field.getCells()[i];
+            const BaseCell<C>& c1 = *cells[i];
             bc2.template process2<B,M,P>(c1, metric);
             for (long j=i+1;j<n1;++j) {
-                const BaseCell<C>& c2 = *field.getCells()[j];
+                const BaseCell<C>& c2 = *cells[j];
                 bc2.process11<B,M,P,BinTypeHelper<B>::do_reverse>(c1, c2, metric);
             }
         }
@@ -226,6 +228,9 @@ void BaseCorr2::process(const BaseField<C>& field1, const BaseField<C>& field2,
     Assert(n1 > 0);
     Assert(n2 > 0);
 
+    const std::vector<const BaseCell<C>*>& c1list = field1.getCells();
+    const std::vector<const BaseCell<C>*>& c2list = field2.getCells();
+
 #ifdef _OPENMP
 #pragma omp parallel
     {
@@ -251,9 +256,9 @@ void BaseCorr2::process(const BaseField<C>& field1, const BaseField<C>& field2,
 #endif
                 if (dots) std::cout<<'.'<<std::flush;
             }
-            const BaseCell<C>& c1 = *field1.getCells()[i];
+            const BaseCell<C>& c1 = *c1list[i];
             for (long j=0;j<n2;++j) {
-                const BaseCell<C>& c2 = *field2.getCells()[j];
+                const BaseCell<C>& c2 = *c2list[j];
                 bc2.process11<B,M,P,false>(c1, c2, metric);
             }
         }