Handling of "dr"

[MikhailRyazanov]

All transform methods have the dr parameter, but they treat it inconsistently:
basex has it in both the high-level basex_transform() and the low-level basex_core_transform(), the scaling is done in basex_core_transform().
dasch has it in _dasch_transform() (aliased by all the high-level <method>_transform()) but not in the low-level dasch_transform(), the scaling is done in _dasch_transform().
direct has only direct_transform(), hopefully working correctly for both forward and inverse.
hansenlaw has only hansenlaw_transform(), hopefully working correctly for both forward and inverse.
linbasex has dr as a "dummy variable for call compatibility with the other methods" and does not use it at all.
onion_bordas has only onion_bordas_transform() and only inverse transform, apparently correct.

This question arose during implementation of the forward transform in BASEX, since it should scale intensities by multiplying by dr instead of dividing. Therefore all functions that take dr must also take direction.

In addition to implementing proper dr handling in linbasex (why not?), I see two general approaches to make it consistent:

1. All transform functions (high- and low-level) should take dr (and direction, if implemented), and the scaling is done at the low level.
2. Methods that pre-compute matrices should include the dr effect in these matrices. Then high-level methods need to pass dr to their get_bs_cached(), and low-level methods simply apply given matrices and always produce correct results without the need to know dr and the transform direction.

The second approach seems to be cleaner and should work slightly faster for repetitive transforms. For default parameters (dr=1), the results will not change in either case.

[stggh]

The role of dr is purely for Jacobian scaling. I think in the past I may have been confused by this, despite @DanHickstein's gentle guidance.

proper dr handling in linbasex (why not?)

The linbasex code was adapted from Thomas Gerber's jupyter python note book, just before his retirement. It is a little bit of a black box. You are right dr should scale the inverse Abel transform in the same way as the other methods.

1. All PyAbel methods do take the direction parameter.

2. dr effect in these matrices.

   I am not sure this is true. The PyAbel methods operate on pixels, dr fixes any intensity scaling(?). There should be no change to the basis arrays(?)

[MikhailRyazanov]

@stggh, I mean

1. The low-level dasch does not take dr (and direction, which is not needed). If somebody uses get_bs_cached() and dasch_transform(), there is no way to take dr into account except manually dividing the dasch_transform() output.
2. The basis matrices are calculated assuming that r_i = i and x_j = j. If we scale these variables by dr, this dr will appear in the Jacobian for the corresponding basis integrals. In any case, get_bs_cached() functions for basex and dasch actually return the transform matrices, and the transform should include the Jacobian regardless of whether it comes from computing the basis in the scaled coordinates or from scaling the problem to the "standard" basis.

[stggh]

1. Ok, so either dr should be passed to dasch_transform() or the function be hidden. This function is only present for the benchmark timing.
2. There is no change to the basis arrays, only the final transform intensity is scaled by dr, right?

[MikhailRyazanov]

1. Another possibility is to pass dr to get_bs_cached(), then dasch_transform() would not need it.
   If you say that the low-level transform functions were needed only for benchmarks — I suppose that was done to avoid the basis-loading overhead? But since now all methods use memory caching, the high-level transform functions should have negligible overhead and can be benchmarked directly, so the low-level transform functions can be eliminated completely, moving that single matrix multiplication to the high-level functions?
2. The basis arrays are defined in a "canonical form" at integer points, like D in equation (1) and W in (10) in the Dasch article, so yes, they do not depend on dr. The transform matrices are (1/Δr) D in (1) and Δr W in (10), so it makes sense to include dr there.
   Another consideration is that BASEX includes the σ factor (which determines the spacing between basis functions, so is somewhat similar to Δr) in its basis matrices. (The article also defines Δ as the pixel pitch, but does not use it anywhere, implicitly assuming that Δ = 1.)

[stggh]

1. avoid the basis-loading overhead

   Yes. Although, I did not write the abel.benchmark code, only hacked it occasionally. The cleverness is due to @rth and @DanHickstein.

   Coming from hansenlaw I consider it to be almost cheating, to partition the calculation (basis/non-basis) for timing purposes. However, for the time consuming basis methods, I "understand".

   high-level transform functions should have negligible overhead and can be benchmarked directly,

   For the benchmarks, the basis calculation time information would still need to be determined.

2. pass dr to get_bs_cached() ... it makes sense to include dr there.

   A good argument. I agree!

[MikhailRyazanov]

For the benchmarks, we still have the low-level _bs_... functions that do basis calculation without any caching. But I need to think about BASEX benchmarks, since there are basis calculations (very slow) and transform calculations (comparable to applying these transforms to images).

[DanHickstein]

Sorry that I missed this discussion earlier. dasch_transform() is basically not a real transform method, and so I think that it's fine if it doesn't obey the full requirements for a standard transform method. It probably should start with an underscore to indicate that it's intended as a private method for the dasch.py namespace. We can still access it from benchmark.py if it has an underscore. A single underscore is mostly a suggestion that the method is not something that most users should be calling.

[MikhailRyazanov]

There is _dasch_transform() already, which is actually a higher-level method than dasch_transform(). Maybe it would make sense to rename dasch_transform() to dasch_core_transform() (like basex_core_transform() in abel.basex, or both of them to just core_transform() within each module) and then _dasch_transform() to dasch_transform()?

In any case, passing dr to get_bs_cached() is probably reasonable. I've done this in basex (when implementing the forward transform), and it works fine there. It also avoids dividing each image by dr, but might complicate the caching logic a little bit.

[MikhailRyazanov]

So #240 has moved dr handling in BASEX from basex_core_transform() to get_bs_cached().

I do not know whether the Dasch methods would benefit from the same approach.

And linbasex still needs some dr handling. By the way, the total scaling for velocity distributions should be dr², I think.