Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

fix: use correct antpairs when finding bad pairs #53

Merged
merged 1 commit into from
Aug 26, 2024
Merged

fix: use correct antpairs when finding bad pairs #53

Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
141 changes: 109 additions & 32 deletions hera_notebook_templates/notebooks/lststack.ipynb
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -132,7 +132,7 @@
},
"outputs": [],
"source": [
"fileconf: str = \"/lustre/aoc/projects/hera/h6c-analysis/IDR2/lstbin-outputs/redavg-smoothcal-inpaint/file-config.h5\"\n",
"fileconf: str = \"/lustre/aoc/projects/hera/h6c-analysis/IDR2/lstbin-outputs/redavg-smoothcal-inpaint-500ns-lstcal/file-config.h5\"\n",
"fileidx: int = 380\n",
"blchunk: int = 0\n",
"\n",
Expand Down Expand Up @@ -1081,7 +1081,76 @@
" color='k', lw=2\n",
" )\n",
" \n",
" ax[0,0].legend(ncols=3)"
" ax[0,0].legend(ncols=3)\n",
"\n",
" for axx in ax[-1]:\n",
" axx.set_xlabel(\"Frequency [MHz]\")\n",
"\n",
" \n",
"def make_auto_plot_multi_autos(auto_stacks: list[UVData], lstbin: list[dict]):\n",
" \n",
" fig, ax = plt.subplots(\n",
" len(stackconf.pols), len(auto_stacks), \n",
" sharex=True, squeeze=False, constrained_layout=True,\n",
" figsize=(12, 2 + len(stackconf.pols)*1.5)\n",
" )\n",
"\n",
" for i, (stack, avg) in enumerate(zip(auto_stacks, lstbin)):\n",
" for p, pol in enumerate(stackconf.pols):\n",
" axx = ax[p, i]\n",
" \n",
" for k, t in enumerate(stack.times):\n",
" flg = stack.flags[:, k, :, p]\n",
" d = np.where(flg, np.nan, np.abs(stack.data[:, k, :, p]))\n",
" percentiles = np.nanpercentile(d, [1, 99], axis=0)\n",
" mean = np.nanmean(d, axis=0)\n",
" \n",
" axx.fill_between(\n",
" stack.freq_array / 1e6,\n",
" percentiles[0],\n",
" percentiles[1],\n",
" alpha=0.4,\n",
" **styles[int(t)]\n",
" )\n",
" axx.plot(\n",
" stack.freq_array / 1e6,\n",
" mean,\n",
" label=f\"{int(t)}\" if not p else None,\n",
" **styles[int(t)]\n",
" )\n",
" \n",
" dd = np.where(np.isnan(d), 1000, d)\n",
" p0 = np.where(np.isnan(percentiles[0]), 0, percentiles[0])\n",
" outliers = np.where(np.any(dd < p0, axis=1))[0].tolist()\n",
" \n",
" dd = np.where(np.isnan(d), -1000, d)\n",
" p0 = np.where(np.isnan(percentiles[1]), np.inf, percentiles[1])\n",
" outliers += np.where(np.any(dd > p0, axis=1))[0].tolist()\n",
" outliers = np.unique(outliers)\n",
" \n",
" for idx in outliers:\n",
" plt.plot(\n",
" stack.freq_array / 1e6,\n",
" d[idx],\n",
" lw=1,\n",
" **styles[int(t)]\n",
" )\n",
" \n",
" axx.set_yscale('log')\n",
" axx.set_title(f\"pol={pol}, LST {stackconf.lst_grid[i]*12/np.pi:.3f} hr\")\n",
"\n",
" # plot the mean (over all autos)\n",
" axx.plot(\n",
" stack.freq_array / 1e6,\n",
" np.nanmean(np.where(avg['flags'][:, :, p], np.nan, np.abs(avg['data'][:, :, p])), axis=0),\n",
" label='mean',\n",
" lw=2,\n",
" color='k'\n",
" )\n",
" \n",
" ax[0,0].legend(ncols=3)\n",
" for axx in ax[-1]:\n",
" axx.set_xlabel(\"Frequency [MHz]\")"
]
},
{
Expand All @@ -1094,7 +1163,10 @@
"outputs": [],
"source": [
"if make_plots:\n",
" make_auto_plot(auto_stacks, autos_lstavg)"
" if len(stackconf.autopairs)>1:\n",
" make_auto_plot_multi_autos(auto_stacks, autos_lstavg)\n",
" else:\n",
" make_auto_plot(auto_stacks, autos_lstavg)"
]
},
{
Expand Down Expand Up @@ -1379,16 +1451,15 @@
" subband = slice(max(worst_idx - 50, 0), min(worst_idx+50, len(diff)))\n",
"\n",
" fq = cross_stacks[0].freq_array[band][subband] / 1e6\n",
" ax[kk].plot(\n",
" fq, \n",
" np.abs(\n",
" np.where(\n",
" cross_lstavg[0]['flags'][apidx, band, polidx], np.nan, original_data_mean[0][apidx, band, polidx]\n",
" )[subband]\n",
" ), \n",
" lw=3, ls='-', color='k', label='flagged mean'\n",
" flagged_mean = np.abs(\n",
" np.where(\n",
" np.all(cross_stacks[0].flags[:, apidx, band, polidx], axis=0), np.nan, original_data_mean[0][apidx, band, polidx]\n",
" )[subband]\n",
" )\n",
" ax[kk].plot(fq, np.abs(cross_lstavg[0]['data'][apidx, band, polidx][subband]), lw=3, ls='--', color='red', label='simul. inpaint')\n",
" ax[kk].plot(fq, flagged_mean, lw=3, ls='-', color='k', label='flagged mean')\n",
" inpmean =np.abs(cross_lstavg[0]['data'][apidx, band, polidx][subband])\n",
" ax[kk].plot(fq, inpmean, lw=3, ls='--', color='red', label='simul. inpaint')\n",
"\n",
"\n",
" #ax[kk].plot(fq, np.abs(inpaint_dpss_models[0][blpol][band][subband]), lw=2, color='purple', ls=':', label='model')\n",
"\n",
Expand All @@ -1405,7 +1476,13 @@
" ax[kk].scatter(fq[flg], d[flg], marker='o', edgecolors=styles[jd]['color'], facecolors='none', alpha=0.4, label='flagged datum' if not gotlabel else None)\n",
" gotlabel = True\n",
" ax[kk].legend(ncols=2)\n",
" ax[kk].set_title(str(blpol))\n",
" ax[kk].set_title(f\"({int(blpol[0])}, {int(blpol[1])}, {blpol[2]})\")\n",
"\n",
" plotmin = min(np.nanmin(flagged_mean), np.nanmin(inpmean))\n",
" plotmax = max(np.nanmax(flagged_mean), np.nanmax(inpmean))\n",
" rng = plotmax - plotmin\n",
" ax[kk].set_ylim(plotmin - rng/5, plotmax + rng/5)\n",
"\n",
"\n",
" kk += 1\n",
" fig.suptitle(\"Examples of Biggest Differences in Inpainted Solutions\")"
Expand Down Expand Up @@ -2534,10 +2611,9 @@
" for line in range(0, 1536, 200):\n",
" axx.axvline(meta.freq_array[line]/1e6, color='gray', alpha=0.4)\n",
" axx.set_xlim(freqs[0], freqs[-1])\n",
" \n",
" all_figs.append(fig)\n",
" \n",
" return all_figs "
"\n",
" yield fig\n",
" plt.close(fig)\n"
]
},
{
Expand All @@ -2552,8 +2628,8 @@
"def get_worst_mean_over_each_band(zscores, n=1, autopols_only=True):\n",
" bad_fellas = {}\n",
"\n",
" nights0 = [data_jd_ints.index(jd) for jd in zscores[0].time_array[::zscores[0].Nbls].astype(int)]\n",
" nights1 = [data_jd_ints.index(jd) for jd in zscores[1].time_array[::zscores[1].Nbls].astype(int)]\n",
" nights0 = [data_jd_ints.index(jd) for jd in zscores[0].times.astype(int)]\n",
" nights1 = [data_jd_ints.index(jd) for jd in zscores[1].times.astype(int)]\n",
" \n",
" if autopols_only:\n",
" polidx = [i for i, p in enumerate(stackconf.pols) if len(set(p))==1]\n",
Expand All @@ -2562,15 +2638,15 @@
" \n",
" npols = len(polidx)\n",
" \n",
" newmeans = {band: np.ones((len(zscores), len(data_jd_ints), len(stackconf.antpairs), npols))*np.nan for band in metrics['band_reduced_mean']}\n",
" newmeans = {band: np.ones((len(zscores), len(data_jd_ints), len(zscores[0].antpairs), npols))*np.nan for band in metrics['band_reduced_mean']}\n",
" \n",
" for band, zsqs in metrics['band_reduced_mean'].items():\n",
" # zsqs is length(lstbins), where each is an array of shape (nights, antpairs, pols)\n",
" # however, the number of nights for each lstbin could be different, so make them the same here....\n",
" newmeans[band][0, nights0] = zsqs[0][..., polidx]\n",
" newmeans[band][1, nights1] = zsqs[1][..., polidx]\n",
" \n",
" lst_night_bl_pols = [(lst, jd, bl + (pol,)) for lst in range(len(zscores)) for jd in data_jd_ints for bl in stackconf.antpairs for j, pol in enumerate(stackconf.pols) if j in polidx]\n",
" lst_night_bl_pols = [(lst, jd, bl + (pol,)) for lst in range(len(zscores)) for jd in data_jd_ints for bl in zscores[0].antpairs for j, pol in enumerate(stackconf.pols) if j in polidx]\n",
" \n",
" for band, zsq in newmeans.items():\n",
" zsq = np.where(np.isnan(zsq.flatten()), -1, zsq.flatten())\n",
Expand All @@ -2587,7 +2663,7 @@
" \n",
" bad_fellas[(lst, bl)].append((jd, z, fr\"Worst $Z^2$ in band {band[0]}-{band[1]}\", band))\n",
"\n",
" return bad_fellas\n"
" return bad_fellas"
]
},
{
Expand Down Expand Up @@ -2654,8 +2730,8 @@
"source": [
"def get_worst_continuous_bad_zscore(zscores, n=1, autopols_only=True):\n",
" bad_fellas = {}\n",
" nights0 = [data_jd_ints.index(jd) for jd in zscores[0].time_array[::zscores[0].Nbls].astype(int)]\n",
" nights1 = [data_jd_ints.index(jd) for jd in zscores[1].time_array[::zscores[1].Nbls].astype(int)]\n",
" nights0 = [data_jd_ints.index(jd) for jd in zscores[0].times.astype(int)]\n",
" nights1 = [data_jd_ints.index(jd) for jd in zscores[1].times.astype(int)]\n",
" \n",
" smallbands = [b for b in bands_considered if b[1] - b[0] <= 200]\n",
" \n",
Expand All @@ -2667,7 +2743,7 @@
" npols = len(polidx)\n",
"\n",
" newmeans = np.ones(\n",
" (len(smallbands), len(zscores), len(data_jd_ints), len(stackconf.antpairs), npols)\n",
" (len(smallbands), len(zscores), len(data_jd_ints), len(zscores[0].antpairs), npols)\n",
" )*np.nan\n",
" \n",
" for i, band in enumerate(smallbands):\n",
Expand All @@ -2677,7 +2753,7 @@
" newmeans[i, 0, nights0] = zsqs[0][..., polidx]\n",
" newmeans[i, 1, nights1] = zsqs[1][..., polidx]\n",
"\n",
" lst_night_bl_pols = [(lst, jd, bl + (pol,)) for lst in range(len(zscores)) for jd in data_jd_ints for bl in stackconf.antpairs for j, pol in enumerate(stackconf.pols) if j in polidx]\n",
" lst_night_bl_pols = [(lst, jd, bl + (pol,)) for lst in range(len(zscores)) for jd in data_jd_ints for bl in zscores[0].antpairs for j, pol in enumerate(stackconf.pols) if j in polidx]\n",
"\n",
" zsq = np.nanmin(newmeans, axis=0)\n",
" \n",
Expand Down Expand Up @@ -2807,16 +2883,17 @@
" freq_ranges = [sum((vv[-1] for vv in v), start=()) for v in badstuff.values()]\n",
" freq_ranges = [(min(v), max(v)) for v in freq_ranges]\n",
"\n",
" plot_visibilities_per_type(\n",
" lstbin_blpols= list(badstuff.keys()), \n",
" for fig in plot_visibilities_per_type(\n",
" lstbin_blpols= list(badstuff.keys()),\n",
" stacks= cross_stacks,\n",
" stats= cross_stats,\n",
" comments=[\"\\n\".join([f\"{vv[-2]}: {vv[0]}\" for vv in v]) for v in badstuff.values()],\n",
" freq_range=freq_ranges,\n",
" alpha=0.5,\n",
" zscores=zsquare,\n",
" auto_stats=auto_stats\n",
" );"
" ):\n",
" plt.show()"
]
},
{
Expand Down Expand Up @@ -2906,9 +2983,9 @@
],
"metadata": {
"kernelspec": {
"display_name": "Python 3 (ipykernel)",
"display_name": "h6c",
"language": "python",
"name": "python3"
"name": "h6c"
},
"language_info": {
"codemirror_mode": {
Expand All @@ -2920,7 +2997,7 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.9.6"
"version": "3.10.6"
},
"toc": {
"base_numbering": 1,
Expand Down
Loading