Ac/power spectrum map #245
Conversation
Estimate power spectrum from ringmap
| @@ -2534,6 +2534,229 @@ def freq(self): | |||
| return self.attrs["freq"] | |||
|
|
|||
|
|
|||
| class SpatialDelayCube(ContainerBase): | |||
There was a problem hiding this comment.
The containers with a delay axis should just inherit from DelayContainer instead of re-implementing the delay axis
| return self.index_map["pol"] | ||
|
|
||
| @property | ||
| def k_parallel(self): |
There was a problem hiding this comment.
I think that some of these index_map entries should actually be datasets, including u, v, uv_mask, and k_parallel. These datasets would then have axes kx, ky, kx,ky`, etc...
My reasoning for this is that index maps should be used for fixed axis coordinates, whereas these are effectively data points at each kx, ky, etc... which depend on the cosmology being used.
| """Get the pol axis.""" | ||
| return self.index_map["pol"] | ||
|
|
||
| @property |
There was a problem hiding this comment.
Same comment as above, I think a few properties can be removed
| return self.attrs["freq_center"] | ||
|
|
||
|
|
||
| class Powerspec3D(ContainerBase): |
There was a problem hiding this comment.
I think the naming of these should be a bit more formal, like PowerSpectrum3D, PowerSpectrum2D, etc...
| _dataset_spec: ClassVar = { | ||
| "ps3D": { | ||
| "axes": ["pol", "delay", "kx", "ky"], | ||
| "dtype": float, |
There was a problem hiding this comment.
| "dtype": float, | |
| "dtype": np.float64, |
| vis_cube.attrs["volume"] = vol_cube | ||
|
|
||
| # Create index map for (u,v) coordinates and k_parallel | ||
| fourier_axes = ["u", "v", "k_parallel", "uv_mask"] |
There was a problem hiding this comment.
I mention this in the containers file as well, but I think that these should actually be datasets in the container rather than index maps, since two containers with the same kx ky axes could have different u and v depending on the cosmology used (I think).
It's also cleaner since index_map entries are generally expected to be used by at least one dataset, so having some extra axes doesn't fit into the container code design as well.
|
|
||
| # Save the spatial window as an attribute | ||
| if self.apply_spatial_window: | ||
| vis_cube.attrs["window_spatial"] = self.spatial_window |
There was a problem hiding this comment.
I think this can also be a dataset instead of an attribute.
| for lde, de in vis_cube.data_tau_uv[:].enumerate(axis=1): | ||
| slc = (pp, lde, slice(None), slice(None)) | ||
| data = np.ascontiguousarray(data_view.local_array[slc]) | ||
| data_uv, NEB_ra, NEB_dec = image_to_uv( |
There was a problem hiding this comment.
Like mentioned above, this could probably be sped up using caput.fftw
| ps_norm = volume_cube * NEB | ||
|
|
||
| # Initialise the 3D power spectrum container | ||
| ps_cube = containers.Powerspec3D( |
There was a problem hiding this comment.
This could just be axes_from=data_1
|
|
||
| # Create index map for (u,v) coordinates,k_parallel axis | ||
| # and uv-mask | ||
| fourier_axes = ["u", "v", "k_parallel", "uv_mask"] |
There was a problem hiding this comment.
Same thing as mentioned before that these should probably be datasets
This will estimate power spectrum (3D,2D,1D) from a ringmap. The short descriptions of tasks are mentioned below:
This task is for check. The 1D estimate from 3D to 1D or 2D to 1D should match. I checked that and it is consistent.
A example run on HI signal is stored here : /project/rpp-chime/arnab92/work/21cm/halofit_new/ringmap/masked_map/ps_HI
Follow the config file there. It will take ~3mins to run. One can also just run inside a login node by caput-pipeline run *.yaml, for time saving.