Implement polarization correction

[SimonHeybrock]

This adds function to apply the inverse of the polarizer and analyzer matrices. I think this should work for both supermirror and He3, I kept this part generic. There are small helper function which serve as glue (connectors) that can be used to, e.g., specify that the TransmissionFunction[Analyzer] should use He3TransmissionFunction[Analyzer] (which will then result in calling the entire workflow).

Tests are basic, I expect more changes. We will also need to look into performance and memory consumption, there are definitely ways to improve this somewhat.

~Draft for now, since there are more TODOs than anticipated initially.~

Fixes #57. Fixes #24.

[jokasimr]

My feeling is that it becomes quite complicated to split everything into different steps like this. Are the intermediate results here really interesting?

The alternative would be to do it in a single step, something like (pseudo code):

def _elements(channels, transmission):
    p = np.array([transmission(c, 'plus') for c in channels])           
    m = np.array([transmission(c, 'minus') for c in channels])           
    d = p*2 - m**2                                                                                                                                                                                                                      
    p /= d  
    m /= d  
    m *= -1 
    return p, m                                                                                                                                                                   

def correct_polarization(channels: SpinChannelIntensities, transmission_polarizer: TransmissionFunction[Polarizer], transmission_analyzer: TransmissionFunction[Analyzer]) -> PolarizationCorrectedIntensities:

    pp, pm = _elements(channels, transmission_polarizer)                 
    ap, am = _elements(channels, transmission_analyzer)
    # pp : polarizer spin 'plus'
    # pm : polarizer spin 'minus'
    # etc               

    mat = [                                                             
        [ap[0]*pp[0], am[1]*pp[1], ap[2]*pm[2], am[3]*pm[3]],
        [am[0]*pp[0], ap[1]*pp[1], am[2]*pm[2], ap[3]*pm[3]],
        [ap[0]*pm[0], am[1]*pm[1], ap[2]*pp[2], am[3]*pp[3]],
        [am[0]*pm[0], ap[1]*pm[1], am[2]*pp[2], ap[3]*pp[3]],
    ]

    res = zeros_like(channels)
    for row in mat:
        for col in mat:
            res[row] += channels[col] * mat[row, col]

    return res

Are there reasons for not doing it in that way?

[SimonHeybrock]

Are there reasons for not doing it in that way?

• Memory consumption.
• Exploiting parallelism.
• Computing results if not all channels have been measured yet.
• Avoiding code duplication and loops which may result in more bugs?
• Ability to histogram each channel before combining, for less memory consumption.

[SimonHeybrock]

My feeling is that it becomes quite complicated to split everything into different steps like this.

Note also that it is further simplified in the follow-up PR (despite adding more matrices into the mix).

[jokasimr]

Memory consumption.

This doesn't seem right to me. Why would it consume more memory? If anything computing the components separately and summing later seems like it should consume 4x more memory (assuming the channels being combined have not been reduced before the correction is applied).

But it doesn't matter. I saw that it was already simplified in the follow up PR you mentioned.

[SimonHeybrock]

Memory consumption.

This doesn't seem right to me. Why would it consume more memory? If anything computing the components separately and summing later seems like it should consume 4x more memory (assuming the channels being combined have not been reduced before the correction is applied).

Several reasons. First and foremost:

• By performing polarization correction, we get 4x more events.
• Concatenating event list temporarily consumes more memory (need all inputs as was as the output), i.e., there is a higher peak consumption.
• If performed naively as written above, also also duplicate all the events coords, even though they are the same. Writing it in a way that avoids this requires significantly more code.

But it doesn't matter. I saw that it was already simplified in the follow up PR you mentioned.

It does matter, since it is important that everyone can understand how, e.g. event data works!