Open Fil opened 4 years ago
@toja I suspect there may be faster possible implementations than either of these. I haven’t looked very closely but do you know of https://dlmf.nist.gov/19.36 ?
@jrus Ah no, I wasn't aware of it. But while I was working on the topic I also implemented Carlson’s elliptic integral of the first kind, RF(x, y, z) after Numerical Recipes, p. 312 for comparison, but it was slower.
I also came across T. Fukushima, Fast computation of incomplete elliptic integral of first kind by half argument transformation. Numer. Math. 116 (4), pp. 687–719, that is mentioned at https://dlmf.nist.gov/19.36. I want to give it a try at some point, but it requires a lot more code than Bulirsch's algorithm.
I’m pretty confused about the code:
var k_ = (sqrt2 - 1) / (sqrt2 + 1), k = sqrt(1 - k_ * k_), K = ellipticF(halfPi, k * k),
Why do we take the square root and then immediately square it, never using the square rooted version again? Why aren’t these constants (k_
, k*k
, K
) all precomputed?
Here’s one paper, but I’m unconvinced it actually has the fastest possible methods for this. https://doi.org/10.1080/2151237X.2009.10129277
The conformal sphere (with 4 slits from one pole to the equator) -> square projection can probably be most cleanly thought of starting from the stereographic projection. From there it is a map from the quarter of the complex plane with positive real and imaginary parts to the square where the 4 corners come from 0, 1, ∞, i. [With the other 4 quarters of the projection found by copying the original signs.]
Alternately we can think of it as the conformal map from unit disk to to a square with vertices at [1, i, –1, –i ], or from a quarter disk to a triangle with vertices at [0, 1, i], with parts outside the disk filled in by symmetry [invert the magnitude of any points outside the disk, then reflect across the line x + y = 1 at the end]
I should experiment a bit to figure out precisely which of these maps are fastest to approximate. I think we will be able to make a complex rational approximation without needing to call any additional special functions at runtime.
For example, using the technology from
https://arxiv.org/pdf/1911.03696
https://arxiv.org/pdf/1804.08127
Is there any standard in d3-geo-projection for how accurate we want projections to be? i.e. how many digits of precision we are looking for? I expect if we can get away with something accurate to 8 digits, that should be somewhat faster to compute than if we need 13 or 15 digits.
we take the square root and then immediately square it
Good catch!
Is there any standard in d3-geo-projection for how accurate we want projections to be?
I'm not aware of any, and we might want to do a survey. I think it varies a lot between projections, for good and bad reasons. I don't think accuracy is very important in "fancy" projections, which are never going to be used in life or death situations, and speed matters more (esp. for animations/interactions). There is value in having a certain regularity of the derivatives—an approximation which "jumps" might create some Moiré or jagged images.
@jrus If you want to revisit the Guyou / Peirce quincuncial you might want to take a look here: https://observablehq.com/d/fc1d2db223f7ed90
I've started to implement them following Snyder's Album of map projections. Since these version use a static k2 = 0.5 it's possible to use a Chebyshev series to compute the elliptic integral. But there're some problems at the poles and also the inverse needs a lot optimization
I haven’t read this and am headed to sleep, but it looks promising:
http://www.acsel-lab.com/arithmetic/arith22/papers/8664a050.pdf https://scholar.google.com/scholar?cluster=17190421556177330303
I have been trying to make rational-function approximations to the quincuncial projection using the AAA algorithm (https://people.maths.ox.ac.uk/trefethen/AAAfinal.pdf), with pretty good success. If I reduce the domain to half an octant of the (stereographically projected) globe, and then approximate a map returning (1 - p(z))^2, where p is the quincuncial projection, a degree 10 rational approximation is accurate to about 14.5 digits throughout the domain. It takes a degree 9 approximant for the inverse map accurate to about 14.5 digits. [I think these could be made accurate to another digit or two, but would have to be constructed using higher-precision arithmetic, which I can’t easily do from Matlab.]
So these functions are going to require overall a bit of logic and possibly a couple conditional branches for domain reduction, then about 12–15 divisions, and 3 square roots (to perform the single complex square root; if this is a bottleneck it might be possible to do something more optimized than the naive implementation), plus maybe 40 additions and multiplications.
I think they should turn out quite snappy, but I still have to implement some barycentric complex rational function evaluation code in Javascript before I can show you something concrete.
Basically, the idea is if we map between these two shapes, there are no square root singularities nearby, making it much easier for a rational function to handle. The black dots are the nodes used for the barycentric formulation of the approximation, and the red dots are the poles of the approximation.
https://i.imgur.com/4cicDa9.png
(I rotated the output shape by 180° in this picture so it’s easier to see which vertex goes to which.)
The maximum errors around the boundary are ~3e-15 in the forward direction and ~2e-15 in the inverse direction; should be a comparable level of error in the interior.
Still working on this (e.g. I'm going to switch to nodes that are slightly outside the reduced domain so that there’s no accidental chance of a NaN from evaluating Infinity / Infinity), but feel free to play with it:
https://observablehq.com/d/89028e12f98e34aa
I was a bit optimistic with my estimate of the number of multiplications/additions. More like 300 of them, plus the 12–14 divisions and 3 square roots, per projected point. (If Javascript could do SIMD, we could use all the SIMD lanes and cut that down by a factor of 4 or whatever, since most of the arithmetic is independent and can be done in parallel.) Still should be pretty quick though I expect.
I’ll try to test the speed out when I get a chance. Edit: it can handle about 9 million projected points per second on my desktop computer.
@toja here’s a fork of your notebook where I only replaced the forward direction but not the inverse
https://observablehq.com/d/47ceb9539492ff76
also ping @Fil
@toja @Fil okay, here’s the inverse as well https://observablehq.com/@jrus/fast-peirce-quincuncial
I'm testing @jrus's implementation in d3-geo-projection.
The (trivial) difference is that the current projection has a default rotate of [-90, -90, 45], which creates a weird default aspect (shown in the difference (black&white) map below). For a nicer aspect, one has to configure it as, e.g. d3.geoPeirceQuincuncial().rotate([-90 + 25, -90, 45]).
With Jacob's implementation the aspect given by the suggested rotate([25, 0, 0]) is perfect, as in this yellow map:
Other than that, the difference shown below (with the new implementation at [-90, 0, 0] to match the original aspect) is minimal—I suppose it comes from the new projection being more precise :)
Notebook: https://observablehq.com/d/6bbd20a7462e68f5
One question if we want to use the new implementation, is to know if we can replace the Guyou projection too, or if we're going to have a diverging code base?
Added fast PeirceQuincuncial in #199
Can definitely use this to implement Guyou as well. I might not have the bandwidth to look at it in the next few days though.
Implementation by @toja — https://observablehq.com/d/9c104687d45ef00e