HERA-Team / hera_pspec

HERA power spectrum estimation code and data formats
http://hera-pspec.readthedocs.io/en/latest/
BSD 3-Clause "New" or "Revised" License
5 stars 3 forks source link

Allow incoherent averaging in time chunks, rather than over all times. #395

Open jsdillon opened 6 months ago

jsdillon commented 6 months ago

Right now, if I have 6 UVPSpec objects with similar times (because they come from interleaves, see #379 ) and I want to incoherently average them together, I have to do something like the following:

# combine all interleaves into UVPSpec object
interleaved_uvp = recursive_add(uvps)

# select each individual time and average interleaves together incoherently
to_recombine = []
for times in interleaved_uvp.time_avg_array.reshape(-1, len(uvps)):
    to_recombine.append(interleaved_uvp.select(times=times, inplace=False))
    to_recombine[-1].average_spectra(time_avg=True, error_weights='P_N', error_field=['P_N', 'P_SN'])

# combine all single-integration UVPSpec objects    
interleaved_uvp = recursive_add(to_recombine)

where

def recursive_add(objects):
    '''Generalized method for faster combination of e.g. UVPSpec objects.'''
    if len(objects) == 0:
        raise ValueError('Cannot run recursive_add on length-0 objects.')
    if len(objects) == 1:
        # Base case: only one object left, return it
        return objects[0]
    elif len(objects) == 2:
        # Base case: two objects, add them together
        return objects[0] + objects[1]
    else:
        # Recursive case: split the list in half and add each half
        midpoint = len(objects) // 2
        left_sum = recursive_add(objects[:midpoint])
        right_sum = recursive_add(objects[midpoint:])
        return left_sum + right_sum

This is awkward and the final concatenation is quite slow, even with recursive adding. It'd be nice to have a way to perform incoherent time averaging that doesn't require you to average over all times in the UVPSpec object.

I had a similar problem in H1C where I wanted to do averages over fields. I had to do a select first, then an average.