From original slack post:

Hi all,
  I spent some time this weekend hacking on ts_wep so it would run on my Mac (does anyone here use it on a Mac already?). I got it mostly working, with the one caveat that a handful of the tests only seem to work when multiprocessing is turned off. Should I make a PR? Jira ticket? Just fork and keep this to myself?

 
While hacking I noticed that the only c++ functionality in ts_wep is essentially zernike and polynomial gradients/evaluation. GalSim has all the equivalent functionality and I think is implicitly a dependency through the stack conda environment. I got a bit nerd-sniped into doing some accuracy and speed comparisons. In speed, GalSim’s functions are typically factors of a few faster than ts_wep. (Mike Jarvis spent a fair amount of time optimizing the GalSim polynomial evaluators; there’s some really esoteric stuff in there now: (https://github.com/GalSim-developers/GalSim/blob/releases/2.3/src/math/Horner.cpp)). No idea yet if this is actually an interesting efficiency boost to the overall wavefront estimation.
 
As for accuracy, everything seems to match quite well (part per billion or better). (Caveat: this is GalSim v2.3. Accuracy is only part per thousand in GalSim v2.2, which is the current version in conda Rubin-env.) ts_wep is very likely formally more accurate in most cases since it hard-codes coefficients instead of computing them numerically like GalSim does. For example, the Hessian determinant (jacobian with order=“2nd”) of circular Z6 comes out to -24.0 in ts_wep and -23.99999998976111 in GalSim. But I doubt the difference matters in practice.
 
The one glaring exception is ZernikeAnnularJacobian order=“2nd” for Noll indices 7&8 (coma), where differences as large as 30% exist. (Still get ~part per billion agreement or better on all other indices). Since the GalSim implementation has no hard-coded numbers in it, this makes me pretty suspicious of the hard-coded coefficients in lines 397-400 of https://github.com/lsst-ts/ts_wep/blob/develop/src/mathcwfs.cc, or maybe line 308.

So anyway, I’d recommend investigating the above and hopefully swapping the ts_wep c++ code for GalSim (and pressuring Rubin-env to upgrade to GalSim v2.3), which would also carry the benefit of making ts_wep pure-python.

I'm attaching a python script and associated plots supporting the above.  I think the biggest question is if there really is a mistake in ZernikeAnnularJacobian with order="2nd".  This depends on my understanding that by "2nd order Jacobian", what is meant is the Hessian determinant.  This appears to be the case from the comparison with Zernikes other than 7 & 8, but I didn't see any better description in the ts_wep or cwfs code itself.  Similarly, I'm assuming that "1st order Jacobian" means Laplacian.

Some other updates:, the speed difference apparently depends a bit on the input array size.  For instance, with 10_000 points, ts_wep ZernikeAnnularJacobian order="1st" is a bit faster (on my mac mini) than GalSim's equivalent (0.0098 seconds for 10 repeats vs. 0.0066 seconds).  But with 100_000 points, GalSim edges out ts_wep 0.0541 to 0.0621.  Other functions are generally 1.5x to 3x faster in GalSim. 

An interesting exception is ZernikeAnnularEval, which is ~15x faster in GalSim.  I took a brief look at the ts_wep code here and discovered it's the only zernike function with trig calls.  This isn't actually necessary, since zernikes always reduce to 2d xy polynomials (despite often being defined in terms of trig functions).  So we should at the very least rewrite that function.

I also compared and timed mathcwfs.poly10_2D and mathcwfs.poly10Grad.  Story here is that agreement is generally floating-point precision limited, and GalSim is 2x to 3x faster.

Finally, I'll note that the GalSim functions all accept arbitrary orders, so could handle a 12th-order 2d polynomial, or Zernike coefficients up to Z66, for example.