From ff38f225fc4954a1dca17ab7e961d95fe67f74e8 Mon Sep 17 00:00:00 2001 From: Josh Dillon Date: Fri, 18 Oct 2024 14:01:58 -0700 Subject: [PATCH] Various clean-ups, fixes, and updates based on trying to run the notebook Some related to Bobby's signal loss stuff, some just prettying things up, a few bug fixes --- ...le_baseline_postprocessing_and_pspec.ipynb | 180 +++++++----------- 1 file changed, 67 insertions(+), 113 deletions(-) diff --git a/notebooks/single_baseline_postprocessing_and_pspec.ipynb b/notebooks/single_baseline_postprocessing_and_pspec.ipynb index fad387c..3d926a9 100644 --- a/notebooks/single_baseline_postprocessing_and_pspec.ipynb +++ b/notebooks/single_baseline_postprocessing_and_pspec.ipynb @@ -7,7 +7,7 @@ "source": [ "# Single Baseline Filtering and Power Spectrum Estimation\n", "\n", - "**by Josh Dillon, Bobby Pascua, and Mike Wilensky**, last updated October 16, 2024\n", + "**by Josh Dillon, Bobby Pascua, and Mike Wilensky**, last updated October 18, 2024\n", "\n", "This notebook is designed to take a single redundantly-averaged unique baseline (typically after LST-binning) and push it through all the way to the power spectrum. It operates on single files that contain a single baseline for all LSTs and both `'ee'` and `'nn'` polarizations. It then can:\n", "* Throw out highly flagged times and/or channels\n", @@ -32,8 +32,8 @@ "# [• Figure 3: Waterfalls After Delay Filtering and/or Inpainting](#Figure-3:-Waterfalls-After-Delay-Filtering-and/or-Inpainting)\n", "# [• Figure 4: First Set of De-Interleaved Waterfalls after Cross-Talk Filtering](#Figure-4:-First-Set-of-De-Interleaved-Waterfalls-after-Cross-Talk-Filtering)\n", "# [• Figure 5: First Set of De-Interleaved Waterfalls after Main-Beam Fringe-Rate Filtering](#Figure-5:-First-Set-of-De-Interleaved-Waterfalls-after-Main-Beam-Fringe-Rate-Filtering)\n", - "# [• Figure 6: First Set of De-Interleaved Waterfalls after Forming Pseudo-Stokes I](#Figure-6:-First-Set-of-De-Interleaved-Waterfalls-after-Forming-Pseudo-Stokes-I)\n", - "# [• Figure 7: First Set of De-Interleaved Waterfalls after Coherent Time Averaging](#Figure-7:-First-Set-of-De-Interleaved-Waterfalls-after-Coherent-Time-Averaging)\n", + "# [• Figure 6: First Set of De-Interleaved Waterfalls after Coherent Time Averaging](#Figure-6:-First-Set-of-De-Interleaved-Waterfalls-after-Coherent-Time-Averaging)\n", + "# [• Figure 7: First Set of De-Interleaved Waterfalls after Forming Pseudo-Stokes I](#Figure-7:-First-Set-of-De-Interleaved-Waterfalls-after-Forming-Pseudo-Stokes-I)\n", "# [• Figure 8: Interleave-Averaged Power Spectra (Pseudo-Stokes I, Q, U, & V) vs. LST](#Figure-8:-Interleave-Averaged-Power-Spectra-(Pseudo-Stokes-I,-Q,-U,-&-V)-vs.-LST)\n", "# [• Figure 9: Interleave-Averaged Power Spectrum SNR vs. LST (Real and Imaginary for pI)](#Figure-9:-Interleave-Averaged-Power-Spectrum-SNR-vs.-LST-(Real-and-Imaginary-for-pI))\n", "# [• Figure 10: High Delay Power Spectrum SNR Histograms Before and After Incoherent Averaging](#Figure-10:-High-Delay-Power-Spectrum-SNR-Histograms-Before-and-After-Incoherent-Averaging)\n", @@ -74,6 +74,7 @@ "import warnings\n", "from astropy import units\n", "from scipy import constants, interpolate, special\n", + "from scipy.signal.windows import dpss\n", "from scipy.stats import norm\n", "from pyuvdata import UVFlag, UVData, UVCal, UVBeam\n", "from pyuvdata import utils as uvutils\n", @@ -114,8 +115,8 @@ "OUT_PSPEC_FILE = os.environ.get(\"OUT_PSPEC_FILE\", SINGLE_BL_FILE.replace('.uvh5', '.pspec.h5'))\n", "\n", "# Band settings\n", - "BAND_STR = os.environ.get(\"BAND_STR\", '50.1~62.2,62.7~73.8,74.6~87.4,108.0~124.5,125.3~136.2,138.3~148.2,'\n", - " '148.5~159.2,159.3~175.2,175.3~189.2,191.5~208.5,208.7~222.9,223.6~231.1') # in MHz\n", + "BAND_STR = os.environ.get(\"BAND_STR\", '50.1~62.2,63.3~73.5,74.6~85.4,108.0~116.1,117.3~124.4,125.3~136.2,138.3~148.2,'\n", + " '150.1~159.2,159.3~169.9,171.9~181.1,181.4~196.4,198.5~208.4,212.3~220.6,224.4~231.1') # in MHz\n", "\n", "# Inpainting settings \n", "ALREADY_INPAINTED = os.environ.get(\"ALREADY_INPAINTED\", \"FALSE\").upper() == \"TRUE\"\n", @@ -144,14 +145,14 @@ " PIXEL_FLAG_CUT = 0.0\n", " \n", "# FRF and time-averaging settings\n", - "NINTERLEAVE = int(os.environ.get(\"EIGENVAL_CUTOFF\", 4)) # number of interleaves to independently FRF\n", + "NINTERLEAVE = int(os.environ.get(\"NINTERLEAVE\", 4)) # number of interleaves to independently FRF\n", "XTALK_FR = float(os.environ.get(\"XTALK_FR\", 0.01)) # Fringe rate half-width in Hz used for fringe rate filtering crosstalk\n", "FR_SPECTRA_FILE = os.environ.get(\"FR_SPECTRA_FILE\", \"/lustre/aoc/projects/hera/zmartino/hera_frf/spectra_cache/spectra_cache_hera_core.h5\")\n", "FR_QUANTILE_LOW = float(os.environ.get(\"FR_QUANTILE_LOW\", 0.05))\n", "FR_QUANTILE_HIGH = float(os.environ.get(\"FR_QUANTILE_HIGH\", 0.95))\n", "FR_EIGENVAL_CUTOFF = float(os.environ.get(\"FR_EIGENVAL_CUTOFF\", 1e-12))\n", - "TARGET_AVERAGING_TIME = 300 # coherent integration time in seconds. Actual time might be less to so that all interleaves have the same number of samples averaged\n", - "USE_CORR_MATRIX = True # Whether or not to use the correlation matrix to do the noise statistics correction factor\n", + "TARGET_AVERAGING_TIME = float(os.environ.get(\"TARGET_AVERAGING_TIME\", 300.0)) # coherent integration time in seconds. Actual time might be less to so that all interleaves have the same number of samples averaged\n", + "USE_CORR_MATRIX = os.environ.get(\"USE_CORR_MATRIX\", \"TRUE\").upper() == \"TRUE\" # Whether or not to use the correlation matrix to do the noise statistics correction factor\n", "\n", "# Power spectrum settings\n", "EFIELD_HEALPIX_BEAM_FILE = os.environ.get(\"EFIELD_HEALPIX_BEAM_FILE\", \"/lustre/aoc/projects/hera/h6c-analysis/IDR2/beams/NF_HERA_Vivaldi_efield_beam_healpix.fits\")\n", @@ -168,14 +169,14 @@ "outputs": [], "source": [ "# Example set of settings\n", - "SINGLE_BL_FILE = '/lustre/aoc/projects/hera/h6c-analysis/IDR2/lstbin-outputs/redavg-smoothcal-inpaint-500ns-lstcal/inpaint/single_baseline_files/zen.LST.baseline.0_4.sum.uvh5' \n", - "OUT_PSPEC_FILE = '/lustre/aoc/projects/hera/mwilensk/H6C/PSPEC/pspec_out/zen.LST.baseline.0_4.sum.pspec.h5'\n", - "ALREADY_INPAINTED = True\n", - "PERFORM_INPAINT = False\n", - "PERFORM_DLY_FILT = False\n", - "USE_BAND_AVG_NSAMPLES = True \n", - "CHANNEL_FLAG_CUT = 0.0\n", - "SAVE_RESULTS = False" + "# SINGLE_BL_FILE = '/lustre/aoc/projects/hera/h6c-analysis/IDR2/lstbin-outputs/redavg-smoothcal-inpaint-500ns-lstcal/inpaint/single_baseline_files/zen.LST.baseline.0_4.sum.uvh5' \n", + "# OUT_PSPEC_FILE = '/lustre/aoc/projects/hera/jsdillon/H6C/PSPEC/pspec_out/zen.LST.baseline.0_4.sum.pspec.h5'\n", + "# ALREADY_INPAINTED = True\n", + "# PERFORM_INPAINT = False\n", + "# PERFORM_DLY_FILT = False\n", + "# USE_BAND_AVG_NSAMPLES = True \n", + "# CHANNEL_FLAG_CUT = 0.0\n", + "# SAVE_RESULTS = False" ] }, { @@ -192,7 +193,7 @@ " 'ALREADY_INPAINTED', 'PERFORM_INPAINT', 'INPAINT_MIN_DLY', 'INPAINT_HORIZON', 'INPAINT_STANDOFF', 'INPAINT_EIGENVAL_CUTOFF',\n", " 'PERFORM_DLY_FILT', 'DLY_FILT_MIN_DLY', 'DLY_FILT_HORIZON', 'DLY_FILT_STANDOFF', 'DLY_FILT_EIGENVAL_CUTOFF',\n", " 'USE_BAND_AVG_NSAMPLES', 'FM_CUT_FREQ', 'PIXEL_FLAG_CUT', 'INTEGRATION_FLAG_CUT', 'CHANNEL_FLAG_CUT', \n", - " 'NINTERLEAVE', 'XTALK_FR', 'FR_SPECTRA_FILE', 'FR_QUANTILE_LOW', 'FR_QUANTILE_HIGH', 'FR_EIGENVAL_CUTOFF', 'TARGET_AVERAGING_TIME',\n", + " 'NINTERLEAVE', 'XTALK_FR', 'FR_SPECTRA_FILE', 'FR_QUANTILE_LOW', 'FR_QUANTILE_HIGH', 'FR_EIGENVAL_CUTOFF', 'TARGET_AVERAGING_TIME', 'USE_CORR_MATRIX',\n", " 'EFIELD_HEALPIX_BEAM_FILE', 'TAPER', 'INCLUDE_INTERLEAVE_AUTO_PS', 'STORE_WINDOW_FUNCTIONS']:\n", " if issubclass(type(eval(setting)), str):\n", " print(f'{setting} = \"{eval(setting)}\"')\n", @@ -748,7 +749,7 @@ " max_contiguous_edge_flags=len(data.times), cache=cache, filter_dims=0)\n", " frf_data[bl] *= 0\n", " frf_data[bl][tslice, :] = np.where(wgts[bl][tslice, :] == 0, 0, d_mdl[tslice, :])\n", - " \n", + "\n", " return frf_data" ] }, @@ -1227,9 +1228,8 @@ " modes = modes[:cut].T * frf_phasor[:,None] # Shape (n_times, n_dpss)\n", " \n", " # Compute the filter matrix according to a weighted least-squares fit.\n", - " return wgts[:,None] * modes @ np.linalg.inv(\n", - " modes.T.conj() @ (wgts[:,None] * modes)\n", - " ) @ (wgts[:,None] * modes).T.conj()" + " modes = wgts[:,None] * modes\n", + " return modes @ np.linalg.inv(modes.T.conj() @ modes) @ modes.T.conj()" ] }, { @@ -1345,19 +1345,17 @@ " deint_xtalk_filt_data.append(d)\n", " deint_frf_data.append(d)\n", " \n", - "\n", - " \n", " # Coherent time-averaging, rephasing to a set of lsts that's consistent across interleaves\n", " dlst = [target_lsts[i] - l for i in range(n_avg_int) for l in np.unwrap(d.lsts)[i * Navg:(i+1) * Navg]]\n", " rephased_frf_d = copy.deepcopy(deint_frf_data[-1])\n", " utils.lst_rephase(rephased_frf_d, bl_vec, d.freqs, dlst, lat=hd.telescope_location_lat_lon_alt_degrees[0], inplace=True)\n", "\n", " pstokes_pols = sorted([pol for pol in rephased_frf_d.pols() if utils.polstr2num(pol, x_orientation=hd.x_orientation) > 0])\n", - " avg_d, avg_f, avg_n = timeavg_data(rephased_frf_d, f, n, Navg=int(AVERAGING_TIME // (dt * NINTERLEAVE)), \n", - " pols=[\"ee\", \"nn\", \"ne\", \"en\"], rephase=False)\n", + " avg_d, avg_f, avg_n = timeavg_data(rephased_frf_d, f, n, Navg=int(AVERAGING_TIME // (dt * NINTERLEAVE)), rephase=False)\n", " \n", - " # Form pseudo-Stokes I and Q from ee and nn\n", + " # Form pseudo-Stokes I, Q, U, and V from ee and nn\n", " form_pseudostokes(avg_d, avg_f, avg_n)\n", + "\n", " # Store Results\n", " deint_avg_data.append(avg_d)\n", " deint_avg_flags.append(avg_f)\n", @@ -1369,7 +1367,9 @@ "id": "a560e664", "metadata": {}, "source": [ - "# Calculate FRF noise covariances (set USE_CORR_MATRIX=True in globals cell to run)" + "## Calculate FRF noise covariances (set USE_CORR_MATRIX=True in globals cell to run)\n", + "\n", + "### TODO: Mike, what's the minimal set of the following we actually need in this notebook?" ] }, { @@ -1377,7 +1377,7 @@ "id": "6fae5075", "metadata": {}, "source": [ - "## First check that the alternate FRF operator doesn't do something insane by comparing a few streams at a given spectral window to the pipeline results" + "### First check that the alternate FRF operator doesn't do something insane by comparing a few streams at a given spectral window to the pipeline results" ] }, { @@ -1395,7 +1395,7 @@ "def get_frop_wrapper(pol=\"nn\", stream_ind=0, band_ind=0, t_avg=AVERAGING_TIME, \n", " rephase=True, wgt_tavg_by_nsample=True, nsamples=None, bl_vec=None,\n", " dlst=None, coherent_avg=True):\n", - " \n", + " '''TODO: document'''\n", " bl=(ANTPAIR[0], ANTPAIR[1], pol)\n", " band = bands[band_ind]\n", " band_slice = band_slices[band_ind]\n", @@ -1419,7 +1419,7 @@ "def get_alt_frf_dat(pol=\"nn\", stream_ind=0, band_ind=0, t_avg=AVERAGING_TIME, \n", " rephase=True, wgt_tavg_by_nsample=True, nsamples=None, bl_vec=None,\n", " dlst=None, coherent_avg=True):\n", - " \n", + " '''TODO: document'''\n", " bl=(ANTPAIR[0], ANTPAIR[1], pol)\n", " band = bands[band_ind]\n", " band_slice = band_slices[band_ind]\n", @@ -1435,13 +1435,11 @@ " return deint_frf_data_test, alt_frf_data\n", "\n", "def plot_streams(test_data, alt_data, band_ind=0, amp=True, linthresh=1e-6):\n", - "\n", + " '''TODO: document'''\n", " band_slice = band_slices[band_ind]\n", "\n", " fig, ax = plt.subplots(nrows=1, ncols=4, figsize=[12, 8])\n", " \n", - " \n", - "\n", " for stream_ind in range(4):\n", " if amp:\n", " dat_plot = np.abs(alt_data[stream_ind] / test_data[stream_ind]) - 1\n", @@ -1477,7 +1475,7 @@ " return\n", "\n", "def get_alt_coavg_dat(stream_ind=0, band_ind=0):\n", - "\n", + " '''TODO: document'''\n", " dlst = [target_lsts[i] - l for i in range(n_avg_int) for l in np.unwrap(deint_filt_data[stream_ind].lsts)[i * Navg:(i+1) * Navg]]\n", " dlst = np.array(dlst).flatten()\n", "\n", @@ -1492,8 +1490,7 @@ " elif hd.pol_convention == \"avg\":\n", " alt_frf_pI = 0.25 * (alt_frf_nn + alt_frf_ee)\n", " \n", - " return alt_frf_pI\n", - " " + " return alt_frf_pI" ] }, { @@ -1531,13 +1528,10 @@ "cell_type": "code", "execution_count": null, "id": "0562df6d", - "metadata": { - "scrolled": false - }, + "metadata": {}, "outputs": [], "source": [ "if USE_CORR_MATRIX:\n", - "\n", " plot_streams(deint_frf_data_test, alt_frf_data, amp=True, linthresh=1e-7)" ] }, @@ -1553,9 +1547,7 @@ "cell_type": "code", "execution_count": null, "id": "dafdcc8a", - "metadata": { - "scrolled": false - }, + "metadata": {}, "outputs": [], "source": [ "if USE_CORR_MATRIX:\n", @@ -1578,7 +1570,7 @@ "id": "bffaaeda", "metadata": {}, "source": [ - "# Fractional error after FRF + coherent average + convert to pI" + "## Fractional error after FRF + coherent average + convert to pI" ] }, { @@ -1589,13 +1581,9 @@ "outputs": [], "source": [ "if USE_CORR_MATRIX:\n", - "\n", - "\n", - "\n", " dat_plot = np.abs(alt_frf_pI/deint_avg_data[0][(0, 4, \"pI\")][:, band_slices[0]] - 1)\n", " coherent_lsts = lst_rad_to_hours(deint_avg_data[0].lsts)\n", "\n", - "\n", " extent = [bands[0][0], bands[0][1], coherent_lsts[-1], coherent_lsts[0]]\n", " im = plt.matshow(dat_plot, norm=matplotlib.colors.LogNorm(), extent=extent)\n", " plt.gca().xaxis.set_label_position('top')\n", @@ -1614,7 +1602,7 @@ " axHist.set_xlim([10**amin, 10**amax])\n", " axHist.axvline(np.median(dat_plot), linestyle=\"--\", color=\"black\")\n", " axHist.axvline(np.quantile(dat_plot, norm.cdf(-2)), linestyle=\":\", color=\"black\")\n", - " axHist.axvline(np.quantile(dat_plot, norm.cdf(2)), linestyle=\":\", color=\"black\")\n" + " axHist.axvline(np.quantile(dat_plot, norm.cdf(2)), linestyle=\":\", color=\"black\")" ] }, { @@ -1622,7 +1610,7 @@ "id": "79980689", "metadata": {}, "source": [ - "# Close enough for now. Let's press on to covariance calculation." + "## Close enough for now. Let's press on to covariance calculation." ] }, { @@ -1659,7 +1647,6 @@ " nsamples=deint_nsamples[stream_ind][cross_antpairpol][:, freq_slice],\n", " dlst=dlst, bl_vec=bl_vec[cross_antpairpol])\n", " \n", - " \n", " cov_dict[band][pol] = frf.get_FRF_cov(frop_dict[band][pol], var_dict[band][pol])\n", " cov_dict[band][\"pI\"] = cov_dict[band][\"ee\"] + cov_dict[band][\"nn\"]\n", " if hd.pol_convention == \"avg\":\n", @@ -1669,18 +1656,6 @@ " deint_covs.append(cov_dict)" ] }, - { - "cell_type": "code", - "execution_count": null, - "id": "8776d437", - "metadata": {}, - "outputs": [], - "source": [ - "if USE_CORR_MATRIX:\n", - "\n", - " example_cov = deint_covs[0][bands[0]][\"pI\"]" - ] - }, { "cell_type": "code", "execution_count": null, @@ -1689,10 +1664,9 @@ "outputs": [], "source": [ "if USE_CORR_MATRIX:\n", - "\n", " co_lst_extent = [coherent_lsts[0], coherent_lsts[-1], coherent_lsts[-1], coherent_lsts[0]]\n", " cov_freq_ind = 50\n", - "\n", + " example_cov = deint_covs[0][bands[0]][\"pI\"]\n", " plt.matshow(np.abs(example_cov[cov_freq_ind, :, :]), norm=matplotlib.colors.LogNorm(), \n", " extent=co_lst_extent)\n", " plt.colorbar(label=r\"abs(cov) (Jy$^2$)\")\n", @@ -1709,10 +1683,7 @@ "outputs": [], "source": [ "if USE_CORR_MATRIX:\n", - "\n", " corr = frf.get_corr(example_cov)\n", - "\n", - "\n", " plt.matshow(np.abs(corr[cov_freq_ind]), norm=matplotlib.colors.LogNorm(), extent=co_lst_extent)\n", " plt.colorbar(label=\"abs(Corr. Coeff.)\")\n", " plt.gca().xaxis.set_label_position('top')\n", @@ -1728,10 +1699,7 @@ "outputs": [], "source": [ "if USE_CORR_MATRIX:\n", - "\n", " lst_ind = 68\n", - "\n", - "\n", " plt.title(f\"Correlation at %.1f MHz with LST=%.2f Hours\" % (data.freqs[band_slices[0]][cov_freq_ind] * 1e-6, coherent_lsts[lst_ind]))\n", " plt.plot(coherent_lsts, np.abs(corr[cov_freq_ind, lst_ind, :]))\n", " plt.yscale(\"log\")\n", @@ -1761,9 +1729,9 @@ "metadata": {}, "outputs": [], "source": [ - "corr_factor = frf.get_correction_factor_from_cov(example_cov)\n", - "\n", - "print(f\"Correction factor for example covariance (including edge times): {corr_factor}\")" + "if USE_CORR_MATRIX:\n", + " corr_factor = frf.get_correction_factor_from_cov(example_cov)\n", + " print(f\"Correction factor for example covariance (including edge times): {corr_factor}\")" ] }, { @@ -1771,15 +1739,11 @@ "id": "64bcc708", "metadata": {}, "source": [ - "# Calculate the covariance matrices" - ] - }, - { - "cell_type": "markdown", - "id": "7046d482", - "metadata": {}, - "source": [ - "# For the PS covariance (TODO: implement this if desired):\n", + "## Calculate the covariance matrices\n", + "\n", + "### For the PS covariance (TODO: implement this if desired):\n", + "\n", + "### TODO: Mike, can we move this out of this notebook?\n", "\n", "$p^{ij}_{\\tau t t'} = x^{i}_{\\tau t} (x^{j}_{\\tau t'})^*$\n", "\n", @@ -1915,39 +1879,39 @@ }, { "cell_type": "markdown", - "id": "52232ecf", + "id": "c89de03e", "metadata": {}, "source": [ - "# *Figure 6: First Set of De-Interleaved Waterfalls after Forming Pseudo-Stokes I*" + "# *Figure 6: First Set of De-Interleaved Waterfalls after Coherent Time Averaging*" ] }, { "cell_type": "code", "execution_count": null, - "id": "db91c7fc", + "id": "34aa4f63", "metadata": {}, "outputs": [], "source": [ "# TODO: The deint_frf_data no longer has pI since we go to pI after coherent average\n", "\n", "if PLOT:\n", - " plot_waterfall(deint_frf_data[0], bl=(ANTPAIR + ('ee',)), flags=deint_flags[0], nsamples=deint_nsamples[0], tslice=None)\n", - " plot_real_delay_vs_lst(deint_frf_data[0], bl=(ANTPAIR + ('ee',)), flags=deint_flags[0], xlim=[-3999, 3999], clim=[-2e2, 2e2], linthresh=1, tslice=None)\n", - " plot_dly_vs_fr(deint_frf_data[0], bl=(ANTPAIR + ('ee',)), xlim=[-1999, 1999], clim=[1e0, 1e5], tslice=None)" + " plot_waterfall(deint_avg_data[0], bl=(ANTPAIR + ('ee',)), flags=deint_avg_flags[0], nsamples=deint_avg_nsamples[0], tslice=None)\n", + " plot_real_delay_vs_lst(deint_avg_data[0], bl=(ANTPAIR + ('ee',)), flags=deint_avg_flags[0], xlim=[-3999, 3999], clim=[-2e2, 2e2], linthresh=1, tslice=None)\n", + " plot_dly_vs_fr(deint_avg_data[0], bl=(ANTPAIR + ('ee',)), xlim=[-1999, 1999], clim=[1e0, 1e5], tslice=None)" ] }, { "cell_type": "markdown", - "id": "c89de03e", + "id": "52232ecf", "metadata": {}, "source": [ - "# *Figure 7: First Set of De-Interleaved Waterfalls after Coherent Time Averaging*" + "# *Figure 7: First Set of De-Interleaved Waterfalls after Forming Pseudo-Stokes I*" ] }, { "cell_type": "code", "execution_count": null, - "id": "34aa4f63", + "id": "db91c7fc", "metadata": {}, "outputs": [], "source": [ @@ -1991,20 +1955,11 @@ " avg_hd.polarization_array = np.array([utils.polstr2num(pol) for pol in pstokes_pols])\n", " avg_hd._determine_pol_indexing()\n", " \n", - " avg_data_copy = copy.deepcopy(avg_data)\n", - " avg_flags_copy = copy.deepcopy(avg_flags)\n", - " avg_nsamples_copy = copy.deepcopy(avg_nsamples)\n", - " nonpstokes_pols = [\"ee\", \"nn\", \"ne\", \"en\"]\n", - " for ob in [avg_data_copy,avg_flags_copy,avg_nsamples_copy]:\n", - " bad_keys = []\n", - " for key in ob:\n", - " if key[2] in nonpstokes_pols:\n", - " bad_keys.append(key)\n", - " for key in bad_keys:\n", - " del ob[key]\n", - " # put data into avg_hd\n", - " avg_hd.update(data=avg_data_copy, flags=avg_flags_copy, nsamples=avg_nsamples_copy)\n", - "\n", + " # update pstokes only data in avg_hd\n", + " avg_data_for_update = datacontainer.DataContainer({bl: avg_data[bl] for bl in avg_data if bl[2] in pstokes_pols})\n", + " avg_flags_for_update = datacontainer.DataContainer({bl: avg_flags[bl] for bl in avg_flags if bl[2] in pstokes_pols})\n", + " avg_nsamples_for_update = datacontainer.DataContainer({bl: avg_nsamples[bl] for bl in avg_nsamples if bl[2] in pstokes_pols})\n", + " avg_hd.update(data=avg_data_for_update, flags=avg_flags_for_update, nsamples=avg_nsamples_for_update)\n", " \n", " # add in autocorrelations copies for all antennas for use in noise calculations\n", " auto_aps = [ap for ap in avg_hd.get_antpairs() if ap[0] == ap[1]]\n", @@ -2078,7 +2033,7 @@ "metadata": {}, "outputs": [], "source": [ - "# TODO: this function should probably be graduated into hera_pspec\n", + "# TODO: this function should probably be graduated into hera_pspec or hera_cal\n", "\n", "def dpss_coherent_avg_correction(spw):\n", " '''This function computes an approximate correction to the noise calculation after fringe-rate filtering. It assumes the that number of integrations\n", @@ -2148,7 +2103,7 @@ "for spw, band in enumerate(bands):\n", " if USE_CORR_MATRIX:\n", " # Tiny variation over frequency and streams, just average to get one number\n", - " corr_factor_streams = [np.mean(get_correction_factor_from_cov(deint_covs[stream_ind][band]))\n", + " corr_factor_streams = [np.mean(frf.get_correction_factor_from_cov(deint_covs[stream_ind][band]['pI']))\n", " for stream_ind in range(NINTERLEAVE)]\n", " coherent_avg_correction_factor = np.mean(corr_factor_streams)\n", " else:\n", @@ -2186,7 +2141,6 @@ "uvps_time_avg = [uvp.average_spectra(time_avg=True, error_weights='P_N', error_field=['P_N', 'P_SN'], inplace=False) for uvp in uvps]\n", "\n", "for spw, band in enumerate(bands):\n", - " #coherent_avg_correction_factor = dpss_coherent_avg_correction(spw)\n", " coherent_avg_correction_factor = coherent_avg_correction_factors[spw]\n", " for i, uvp in enumerate(uvps_time_avg): \n", " # loop over all pols, but only this spw\n", @@ -2257,7 +2211,7 @@ " SNRs.append(np.ravel((uvp.get_data(key) / uvp.get_stats('P_N', key).real)[:, high_dlys]))\n", " for i, uvp in enumerate(uvps_time_avg):\n", " high_dlys = np.abs(uvp.get_dlys(key[0]) * 1e9) > 1000\n", - " tavg_SNRs.append(np.ravel((uvp.get_data(key) / uvp.get_stats('P_N', key).real)[:, high_dlys]))\n" + " tavg_SNRs.append(np.ravel((uvp.get_data(key) / uvp.get_stats('P_N', key).real)[:, high_dlys]))" ] }, { @@ -2640,9 +2594,9 @@ ], "metadata": { "kernelspec": { - "display_name": "canary", + "display_name": "Python 3 (ipykernel)", "language": "python", - "name": "canary" + "name": "python3" }, "language_info": { "codemirror_mode": { @@ -2654,7 +2608,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.11.10" + "version": "3.12.0" }, "toc": { "base_numbering": 1,