Should ABD objects optimize for when some fields are omitted?

[duetosymmetry]

We may be in the situation where we use and abd object to store only strain, or just strain and \psi_i through \psi_4 with some smallest value of i. In this case, there are a lot of calculations in e.g. AsymptoticBondiData.transform() that could be skipped. For example, in the lines

scri/scri/asymptotic_bondi_data/transformations.py

Lines 325 to 336 in 6db96ad

	# ψ0(u, θ', ϕ') exp(2iλ)
	ψ0 = sf.Grid(self.psi0.evaluate(distorted_grid_rotors), spin_weight=2)
	# ψ1(u, θ', ϕ') exp(iλ)
	ψ1 = sf.Grid(self.psi1.evaluate(distorted_grid_rotors), spin_weight=1)
	# ψ2(u, θ', ϕ')
	ψ2 = sf.Grid(self.psi2.evaluate(distorted_grid_rotors), spin_weight=0)
	# ψ3(u, θ', ϕ') exp(-1iλ)
	ψ3 = sf.Grid(self.psi3.evaluate(distorted_grid_rotors), spin_weight=-1)
	# ψ4(u, θ', ϕ') exp(-2iλ)
	ψ4 = sf.Grid(self.psi4.evaluate(distorted_grid_rotors), spin_weight=-2)
	# σ(u, θ', ϕ') exp(2iλ)
	σ = sf.Grid(self.sigma.evaluate(distorted_grid_rotors), spin_weight=2)

we can simply skip all the fields which are not stored; and similarly in the lines

scri/scri/asymptotic_bondi_data/transformations.py

Lines 342 to 387 in 6db96ad

	fprime_of_timenaught_directionprime = np.empty((6, self.n_times, n_theta, n_phi), dtype=complex)
	# ψ0'(u, θ', ϕ')
	fprime_temp = ψ4.copy()
	fprime_temp *= ðuprime_over_k
	fprime_temp += -4 * ψ3
	fprime_temp *= ðuprime_over_k
	fprime_temp += 6 * ψ2
	fprime_temp *= ðuprime_over_k
	fprime_temp += -4 * ψ1
	fprime_temp *= ðuprime_over_k
	fprime_temp += ψ0
	fprime_temp *= one_over_k_cubed
	fprime_of_timenaught_directionprime[0] = fprime_temp
	# ψ1'(u, θ', ϕ')
	fprime_temp = -ψ4
	fprime_temp *= ðuprime_over_k
	fprime_temp += 3 * ψ3
	fprime_temp *= ðuprime_over_k
	fprime_temp += -3 * ψ2
	fprime_temp *= ðuprime_over_k
	fprime_temp += ψ1
	fprime_temp *= one_over_k_cubed
	fprime_of_timenaught_directionprime[1] = fprime_temp
	# ψ2'(u, θ', ϕ')
	fprime_temp = ψ4.copy()
	fprime_temp *= ðuprime_over_k
	fprime_temp += -2 * ψ3
	fprime_temp *= ðuprime_over_k
	fprime_temp += ψ2
	fprime_temp *= one_over_k_cubed
	fprime_of_timenaught_directionprime[2] = fprime_temp
	# ψ3'(u, θ', ϕ')
	fprime_temp = -ψ4
	fprime_temp *= ðuprime_over_k
	fprime_temp += ψ3
	fprime_temp *= one_over_k_cubed
	fprime_of_timenaught_directionprime[3] = fprime_temp
	# ψ4'(u, θ', ϕ')
	fprime_temp = ψ4.copy()
	fprime_temp *= one_over_k_cubed
	fprime_of_timenaught_directionprime[4] = fprime_temp
	# σ'(u, θ', ϕ')
	fprime_temp = σ.copy()
	fprime_temp -= ððα
	fprime_temp *= one_over_k
	fprime_of_timenaught_directionprime[5] = fprime_temp

we can skip the calculations for fields \psi_j where j<i, i representing the smallest subscript of Weyl scalar we keep, as above.

If I remember correctly, @akhadse has implemented this as a speedup. Was that right, Akshay?

@moble do you think this is a feature worth having for ABD objects?

[akhadse]

@duetosymmetry - Yes that is right. I check there if self.fields_present[i]==True and skip the calculations for the corresponding field if they are not present

[moble]

I think the objective is good, but it would pay to be thoughtful about the implementation. In particular, it's dangerous to start putting in flags and branches, because it increases code complexity, and we would need to make sure that those flags are supported in every piece of code. Plus, it makes it harder to ensure that what we write reflects the mathematics. I think there's a much easier solution that should be more generic and elegant.

I put some thought into the design of the base AsymptoticBondiData class, specifically with this use case in mind. I decided that it would be best for safety if the base class used the naive full-storage data format (so that we could design and test with minimal complexity), but then abstracted away access to the parts of that data. So obviously, this line may create far more storage than necessary when, e.g., psi0 is essentially 0. But then I define self._psi0 as a ModesTimeSeries that wraps a slice of that big array, and further hide that slice behind setter/getter properties for self.psi0. And only self.psi0 should be used in any code that uses an ABD.

So a cleaner and simpler approach would be to subclass AsymptoticBondiData with a different constructor, where that big array isn't created, then just define self._psi0 as essentially 0, and let numpy handle broadcasting. (Actually, numpy will frequently short-circuit operations involving 0, in which case there's literally no cost.) Now, because we need to be able to do things like self.psi0.evaluate, we still need self._psi0 to be a ModesTimeSeries object. But I wrote that class with this possibility in mind; you can actually create an MTS from a literal 0, even with a nontrivial time array. (Though you'll probably want to specify ell_max.)

The big advantage of this approach is that code implementing math can be generic, free of branches, and actually look like the equations. It won't be precisely as efficient as if you special-case for 0 Newman-Penrose components, but it won't be anywhere close to a bottleneck, in either speed or memory. Obviously, I haven't actually tested these things. And in particular, it looks like ModesTimeSeries.interpolate needs to handle the case where there's just one value. But otherwise, I think it should just work.

[moble]

@duetosymmetry I just realized you may want more efficiency not just on the input to this transform function that you linked to above, but also on the output. That is you may want to just not calculate, e.g., the output psi0 for some cases, regardless of what any input may be. I can see why that could be useful when you only need h, for example. But it strikes me as very dangerous to do that at the generic object level, because people could try to use the output in applications where they actually need psi0. Essentially, if some of the data is intentionally incorrect, I don't think we should consider it an AsymptoticBondiData object.

Still, it might be worth implementing this special behavior just in the transform function, but ensure that the output is a sort of dead-end object that doesn't get to participate in all the rest of the ABD goodness. If this is what you mean, then I think a better approach would be to take the minimum i (for psi_i) as an argument to transform rather than as a property of the input object, and emit a different subclass of ABD, which just defines self._psi0 to raise an error (or maybe just as None, if you want to be able to check for this condition). Then, if anyone tries to actually use the psi0 field for math, they would correctly get an error — but, e.g., sigma would still be accessible. (In fact, such an object could even be passed to transform without any harm as long as that new argument is big enough.)