Question: the origin of SMALLER_EDGE2OPEDGE_DIST table

[hombit]

At LINCC Frameworks, we are working on spatial algorithms using HEALPix grids, most over on cross-matching. The minimum distance between opposite edges and other geometrical properties of the HEALPix tiles are very important to us.

Could you please tell me the origin of this array with minimum distances between opposite edges of a pixel? https://github.com/cds-astro/cds-healpix-rust/blob/d209a61b46969db35081034806d0c3bed9496d33/src/lib.rs#L93

Is there any paper about that? The comment section also lists other properties. Where did they come from?

[fxpineau]

It comes from my own work:

• although I briefly mention those HEALPix developments in this ADASS 2018 poster and associated proceedings, I have not published more than this;
• FWIW, we tested spark-based cross-matches, see e.g. this IVOA 2016 presentation, before I implemented a pure Rust demonstrator (ExXmatch) presented at ADASS 2019: see the presentation and the proceedings. For the radec-tree (and HEALPIx indexed binary file, by blocks of columns), see this ADASS2010 proceeding.

[hombit]

A bit more context: We are working on a Healpix geometry analysis paper and have some preliminary results regarding the minimum distances between edges. We checked the numbers today and found that, at least for some Healpix orders, we can find slightly smaller distances than specified in the code.

The analysis and results are still preliminary, but they look practically useful for this project. We are happy to schedule a meeting or discuss it in any other way you find convenient. Please let me know what you think about it here or by email: malanchev@cmu.edu

[hombit]

@fxpineau thank you! I was writing my last comment just when you posted your one, I'll check all the references!

[fxpineau]

Ok, I will write to you tomorrow. The values in the Rust code come from the Java library. I may dig into my (possibly unpublished?) code to find how I computed it. (I was off for a few days and will be back at work tomorrow).

[fxpineau]

@hombit, I was back from a long trip yesterday, my mind is more clear this morning.

First of all, thank you for your message and for the double-check.

Having too large values in SMALLER_EDGE2OPEDGE_DIST is problematic since it may lead to miss some associations if used in a cross-match. This array is not used inside the Healpix library, it may contain wrong values that have not been spotted yet.

I have two questions:

• how large the discrepancies between the values you computed and the ones in SMALLER_EDGE2OPEDGE_DIST are?
• how did you computed those values?

I dug into my (private) Java code, here an excerpt:

private static double computeCellHashSmallerEdgeToOppositeEdgeDistance(int depth) {
    final double nside = (double) Healpix.nside(depth);
    final F f = new DistInNorthPolarCapNegativeSlope(0d, TRANSITION_LATITUDE, 2 / nside, 0.25 * PI / nside);
    final double xMin = TRANSITION_LATITUDE;
    final double xMax = TRANSITION_LATITUDE + 1 / nside;
    final int nIterMax = 100;
    final double xRelErrTol = 1e-16;
    final double yRelErrTol = 1e-16;
    final double latRad = Minimization.xOfMinimum(f, xMin,  xMax, nIterMax, xRelErrTol, yRelErrTol);
    final double squareOfsinOfhalfDist = f.f(latRad);
    final double dist = 2 * asin(sqrt(squareOfsinOfhalfDist));
    return dist;
  }

  private static final class DistInNorthPolarCapNegativeSlope implements F {

    private final AngularDistanceComputer dComputer;

    private final double cLon;
    private final double cLat;
    private final double cLatCos;
    private final double b;

    private DistInNorthPolarCapNegativeSlope(final double centerLonRad, final double centerLatRad,
        final double b, final double angleScaleInRad) {
      this.dComputer = AngularDistanceComputer.getComputer(angleScaleInRad);
      this.cLon = centerLonRad;
      this.cLat = centerLatRad;
      this.cLatCos = cos(this.cLat);
      this.b = b;
    }
    @Override
    public double f(double x) {
      final double lat = x;
      final double z = sin(lat);
      final double bigZ = sqrt(3 * (1 - z));
      final double y2d = 1 - bigZ;
      // X = Z * (4/pi * lon - 1) + 1
      // y2D = -x2d + b
      final double x2d = this.b - y2d;
      final double lon = ((x2d - 1) / bigZ + 1) * PI * 0.25;
      return dComputer.squareOfsinOfhalfDistInRad2(lon - this.cLon, lat - this.cLat,
          this.cLatCos, cos(lat));
    }
  };

private static final class DistInNorthPolarCapNegativeSlope implements F {

    private final AngularDistanceComputer dComputer;

    private final double cLon;
    private final double cLat;
    private final double cLatCos;
    private final double b;

    private DistInNorthPolarCapNegativeSlope(final double centerLonRad, final double centerLatRad,
        final double b, final double angleScaleInRad) {
      this.dComputer = AngularDistanceComputer.getComputer(angleScaleInRad);
      this.cLon = centerLonRad;
      this.cLat = centerLatRad;
      this.cLatCos = cos(this.cLat);
      this.b = b;
    }
    @Override
    public double f(double x) {
      final double lat = x;
      final double z = sin(lat);
      final double bigZ = sqrt(3 * (1 - z));
      final double y2d = 1 - bigZ;
      // X = Z * (4/pi * lon - 1) + 1
      // y2D = -x2d + b
      final double x2d = this.b - y2d;
      final double lon = ((x2d - 1) / bigZ + 1) * PI * 0.25;
      return dComputer.squareOfsinOfhalfDistInRad2(lon - this.cLon, lat - this.cLat,
          this.cLatCos, cos(lat));
    }
  };

private static final void genSMALLER_EDGE2OPEDGE_DIST() {
    System.out.println("double[] SMALLER_EDGE2OPEDGE_DIST = new double[]{");
    for (int d = 0; d <= Healpix.DEPTH_MAX; d++) {
      double dist = computeCellHashSmallerEdgeToOppositeEdgeDistance(d);
      System.out.print(dist);
      if (d != Healpix.DEPTH_MAX) {
        System.out.println(", // depth = " + d);
      }
    }
    System.out.println("};");
    // A partir d'un momment, faire la flat sky approximation ??
    // center = (0, TRANS_Z)
    // Strait line defined by:
    // P1: x1 = lon1 = 0, y1 = lat(Y=1/nside) => y1 = lat1 = asin(1 - (1 - 1/nside)^2/3 - TRANSITION_LATITUDE);
    // P2: x2 = PI/4 * 1/nside, y2 = lat2 = 0 = TRANSITION_LATITUDE - TRANSITION_LATITUDE
    // => y = (y2-y1)/(x2 - x1) x + y1
    // => trouver le perpendiculaire a cette droite, passant par l'origine
  }

In which Minimization uses the Golden Section Search in One dimansion of the Numerecial Recipes 3, section 10.2.

So:

• if the discrepancies are very small, they may come from numerical inaccuracies (or parameters) in the minimization algorithm
• if they are large, they may come:
  • from a bug in my code
  • from a miss-identification of the location on the sky in which the distance is the smallest one (maybe the location depends on the depth)

You may have find a more robust way to estimate those values (e.g. using path_along_cell_edge at order 29).

Let's continue by email (or zoom, or ...) as you suggested: francois-xavier.pineau@astro.unistra.fr Thanks again.

[fxpineau]

So, the location of the minimum distance depends on the depth. The location is the same for depth 0 and 1, but then it changes. I made a Rust code to compute new values from new assumptions (intuitive assumptions that must be checked more carefully). The starting point of the new assumptions are the values (distance + center of the cell) provided by @hombit for depth from 0 to 14. The new values and the code to compute them is available in commit eb82866a8cadc926ea73c9b37f05e5c7e676e42a.

[hombit]

Thank you, @fxpineau. I propose to keep this issue open until I and @smcguire-cmu finish our analysis and report the results back. The rough time estimate is a month

[fxpineau]

I agree. FYI, I pushed on crates.io an new version of the library with refined values using 100_000_000 sub-segments at orders 1 to 9 (inclusive), 10_000_000 sub-segments at order 10 and 1_000_000 sub-segments from 11 to 15 (inclusive). See here. I am waiting for your analysis to possibly further improve those values (and to confirm that the assumptions are correct). Thank you again.

[hombit]

Sounds great!