Detectors as indices of arrays instead of objects

[dpole]

It would be extremely good (and timely) to take an explicit decision about the data layout.

Current structure

Every simulation owns a list of observations. Every observation owns a detector from which you can get a timestream array, a quaternion, a PSD array, a white noise level etc. Firstly, I would reduce the number of indirections from simulation.observations[3].detector.get_tod() to simulation.observations[3].tod. Anyway, even in the latter case, the paradigm is

sim.observations[3].tod  # 1D np.array
sim.observations[3].wn_level  # float
print(sim.observations[3])

The advantage is that the code reflects the logical hierarchy and organization of the data we are storing. On the practical side, it allows to easily plug in and remove detectors, print or dump to disk all the attributes of the detector. This model is well suited for random operations on random detectors.

Proposal

I propose the following alternative paradigm. The simulation directly owns an array of quaternions, a 2d array of TODs, a 2d array of PSDs, a 1d array of white noise levels etc. These arrays stack the data of the individual detectors.

sim.tod  # 2D np.ndarray: first index for detectors, second for time
sim.tod[3]  # same result as before
sim.wn_level[3]  # same result as before
sim.print_detector(3)  # same result as before

Pros

Operations that involve more than one detector are likely to be easier to write. This is sure if the operation involves all the detectors.
If these operations are done by compiled codes, the interface is much easier (pass one big array, instead of many small values/arrays)
All numpy functions integrate seamlessly with ndarrays
...many other technical/philosophical things that probably shouldn’t drive the decision...

Examples

Standard deviations of the TOD of each detector: sim.tod.std(1)
TOD of the standard deviation across the focal plane: sim.tod.std(0)
TOD of the mean signal across the focal plane: sim.tod.mean(0)
Add white noise to the TOD:
sim.tod += np.random.normal(size=sim.tod.shape) * sim.wn_level[:, np.newaxis]
Add cross-talks given by a crosstalk matrix C: sim.tod = C @ sim.tod

Note

If we think it is necessary, we could still implement our current API with a caveat: sim.observations[3].tod has to be a view of sim.tod[3] (same for other attributes). In practice, it means that sim.observations[3].tod = … is not allowed (you have to do sim.observations[3].tod[:] = …) and sim.observations[3].wn_level is a 1-element array.

You got to the end, you're a hero. One last effort: comment!

[dpole]

I've copied this proposal in a google doc https://docs.google.com/document/d/1aqTuBUrUPZGJx-fGgjH3nfB7kvJOG6Z27H4ybRrgLmQ/edit?usp=sharing
I think it is way more practical for small comments/questions on specific sentences. The caveat, of course, is that only the issue keeps permanent records.

[mreineck]

As I mentioned in the google doc, I don't think that the 2D arrays with (detector, time) indexing will work, because they won't necessarily have the same length ... or is sampling frequency identical for all LiteBIRD detectors?
If my assumption is correct, and 2D arrays cannot be used, we could gain some more flexibility and make Simulation hold a dict of Observations instead of a list.

[dpole]

This is a good point. We have to be able to handle tod at different sampling rates (even if not envisaged at this point). We can gather all the detectors with the same sampling frequencies and time span into one Observation and organize then in arrays.
So, the current organization is

Simulation
|->Observation1   (det 1)
|     |-> sampling frequency
|     |-> tod    (1D array)
|     '-> wn_level
|
|-> Observation2  (det 2)
|     |-> sampling frequency       (possibly the same)
|     |-> tod    (1D array)
|     '-> wn_level
:

I propose

Simulation
|->Observation1 (possibly, but not necessarily, all dets with same sampling freq)
|     |-> sampling frequency
|     |-> tod    (2D array)
|     '-> wn_level (1D array)
|
|-> Observation2  
|     |-> sampling frequency       (possibly different)
|     |-> tod    (2D array)
|     '-> wn_level    (1D array)
:

Even in this case, the TODs that are likely to "communicate" are those at the same sampling frequency: a regular data structure would make this easier and more efficient. Please point out if I'm missing major shortcomings.

[ziotom78]

Mmmm, the more I read your proposal, the more I like it. @dpole , may you please submit a PR with the proposed changes? I think it would be beneficial to have some real code to examine.

[dpole]

Sure, happy to do it. It will take a bit of time though. In the meanwhile all suggestions and comments are helpful.

[dpole]

#42 implemented this