From 606e38454c6d7e6dd9d6d6657af3593ed8f8ae00 Mon Sep 17 00:00:00 2001 From: Paul Chichura Date: Tue, 7 Aug 2018 06:29:28 -0400 Subject: [PATCH 1/9] added delay_wedge() to plot.py --- hera_pspec/plot.py | 175 +++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 175 insertions(+) diff --git a/hera_pspec/plot.py b/hera_pspec/plot.py index 13bdacfa..9af6a743 100644 --- a/hera_pspec/plot.py +++ b/hera_pspec/plot.py @@ -5,6 +5,8 @@ import matplotlib.pyplot as plt import copy from collections import OrderedDict as odict +import astropy.units as u +import astropy.constants as c def delay_spectrum(uvp, blpairs, spw, pol, average_blpairs=False, @@ -494,3 +496,176 @@ def delay_waterfall(uvp, blpairs, spw, pol, component='real', average_blpairs=Fa if new_plot: return fig +def delay_wedge(uvps, blpairs, fold=False, cosmo=True, center_line=True, horizon_lines=True, cosmo_cbar=True, suptitle=''): + """ + Plot a 2D delay spectrum (or spectra) for a group of baselines. + + Parameters + ---------- + uvps : UVPspec + List of UVPSpec objects, containing delay spectra for a set of baseline-pairs, + times, polarizations, and spectral windows. + + blpairs : list of tuples + List of baseline-pair tuples. + + fold : bool, optional + Whether to fold the power spectrum in :math:`|k_\parallel|`. + Default: False. + + cosmo : bool, optional + Whether to convert axes to cosmological units :math:`|k_\parallel|` and + :math:`|k_\perpendicular|`. + Default: True. + + center_line : bool, optional + Whether to plot a dotted line at :math:`|k_\parallel|` =0. + Default: True. + + horizon_lines : bool, optional + Whether to plot dotted lines along the horizon. + Default: True. + + cosmo_cbar : bool, optional + Whether power has been converted to appropriate units through calibration. + Default: True. + + suptitle : string, optional + Suptitle for the plot. If not provided, no suptitle will be plotted. + Default: ''. + + Returns + ------- + f : matplotlib.pyplot.Figure + Matplotlib Figure instance. + """ + + if (type(uvps) != list) and (type(uvps) != np.ndarray): + raise AttributeError("The uvps paramater should be a list of UVPSpec objects.") + + #Initialize axes objects + ncols = len(uvps) + f, axes = plt.subplots( + ncols=ncols, + nrows=1, + sharex=True, + sharey=True, + squeeze=False, + figsize=(20, 10)) + + #Plot each uvp + pols = [] + plt.subplots_adjust(wspace=0, hspace=0.1) + for uvp, ax in zip(uvps, axes.flatten()): + #Find redundant-length baseline groups and their baseline lengths + BIN_WIDTH = 0.3 #redundancy tolerance in meters + NORM_BINS = np.arange(0.0, 10000.0, BIN_WIDTH) + sorted_blpairs = {} + bllens = [] + blpair_seps = uvp.get_blpair_seps() + blpair_array = uvp.blpair_array + for antnums in blpairs: + blpair = uvp.antnums_to_blpair(antnums) + bllen = blpair_seps[np.where(blpair_array==blpair)[0][0]] + bllen = np.round(np.digitize(bllen, NORM_BINS) * BIN_WIDTH, 1) + if bllen in bllens: + sorted_blpairs[bllen].append(antnums) + else: + bllens.append(bllen) + sorted_blpairs[bllen] = [antnums] + bllens = sorted(bllens) + + #Average the spectra along redundant baseline groups and time + uvp.average_spectra(blpair_groups=sorted_blpairs.values(), time_avg=True) + + #Grab polarization string for naming the plot later on + pol_int_to_str = {1: 'pI', 2: 'pQ', 3: 'pU', 4: 'pV', -5: 'XX', -6: 'YY', -7: 'XY', -8: 'YX'} + pol = pol_int_to_str[uvp.pol_array[0]] + pols.append(pol) + + #Format data array + data = uvp.data_array[0][:, :, 0] + if fold: + uvp.fold_spectra() + data = uvp.data_array[0][:, data.shape[1] // 2:, 0] + + #Format k_para axis + x_axis = (uvp.get_dlys(0)*1e9).tolist() + if cosmo: + x_axis = uvp.get_kparas(0).tolist() + if fold: + x_axis.insert(0, 0) + ax.set_xlim((x_axis[0] / 2, x_axis[-1] / 2)) + + #Format k_perp axis + #Calculate kpr indices, and stretch wedgeslices to length of baseline norms + bllen_indices = [int(bllen * 10) for bllen in bllens] + stretched_data = np.zeros((bllen_indices[-1], data.shape[-1]), dtype=np.float64) + j = 0 + for i in range(len(bllen_indices)): + stretched_data[j:bllen_indices[i]] = data[i] + j = bllen_indices[i] + data = stretched_data[...] + + #Find kpr mid-indicies for the tickmarks + bllen_midindices = [] + for i in range(len(bllen_indices)): + if i == 0: + bllen_midindices.append(bllen_indices[i] / 2.) + else: + bllen_midindices.append((bllen_indices[i] + bllen_indices[i - 1]) / 2.) + + #Setting y axis ticks and labels + ax.set_yticks(bllen_midindices) + ax.set_yticklabels(bllens, fontsize=10) + + #Plot the data + im = ax.imshow( + np.log10(np.abs(data)), + aspect='auto', + interpolation='nearest', + extent=[x_axis[0], x_axis[-1], bllen_indices[-1], 0]) + + #Plot the horizon lines if requested + if horizon_lines: + horizons = [(bllen*u.m/c.c).to(u.ns).value for bllen in bllens] + #if cosmo: + #XXX still need to convert horizons to appropriate kpara values + j = 0 + for i, (horizon, bllen) in enumerate(zip(horizons, bllens)): + x1, y1 = [horizon, horizon], [j, bllen_indices[i]] + x2, y2 = [-horizon, -horizon], [j, bllen_indices[i]] + ax.plot(x1, y1, x2, y2, color='#ffffff', linestyle='--', linewidth=1) + j = bllen_indices[i] + + #plot center line at k_para=0 if requested + if center_line: + ax.axvline(x=0, color='#000000', ls='--', lw=1) + + ax.tick_params(axis='both', direction='inout') + ax.set_title(pol, fontsize=15) + + #add colorbar + cbar_ax = f.add_axes([0.9125, 0.25, 0.025, 0.5]) + cbar = f.colorbar(im, cax=cbar_ax) + + #add axis labels + if cosmo: + f.text(0.5, 0.05, r'$k_{\parallel}\ [h\ Mpc^{-1}]$', fontsize=15, ha='center') + f.text(0.07, 0.5, r'$k_{\perp}\ [\rm h\ Mpc^{-1}]$', fontsize=15, va='center', rotation='vertical') + else: + f.text(0.5, 0.05, r'$\tau\ [ns]$', fontsize=15, ha='center') + f.text(0.07, 0.5, r'Baseline Group $[m]$', fontsize=15, va='center', rotation='vertical') + + #add colorbar label + if cosmo_cbar: + cbar.set_label(r"$P(k_{\parallel})\ [mK^2h^{-3}Mpc^3]$", fontsize=15, ha='center') + else: + cbar.set_label(r'$P(\tau)\ [mK^2]$', fontsize=15, ha='center') + + #Add suptitle if requested + if len(suptitle) > 0: + f.suptitle(suptitle, fontsize=15) + + #Return Figure: the axis is an attribute of figure + return f \ No newline at end of file From 06c3ee184ad3d652033c9b41f5483277530b2ad2 Mon Sep 17 00:00:00 2001 From: Paul Chichura Date: Wed, 5 Sep 2018 07:21:47 -0400 Subject: [PATCH 2/9] made revisions to plot.py to match comments --- hera_pspec/plot.py | 192 ++++++++++++++++++++++++++++++--------------- 1 file changed, 130 insertions(+), 62 deletions(-) diff --git a/hera_pspec/plot.py b/hera_pspec/plot.py index 9af6a743..f6a669e3 100644 --- a/hera_pspec/plot.py +++ b/hera_pspec/plot.py @@ -1,6 +1,6 @@ import numpy as np import pyuvdata -from hera_pspec import conversions +from hera_pspec import conversions, uvpspec, utils import matplotlib import matplotlib.pyplot as plt import copy @@ -496,67 +496,123 @@ def delay_waterfall(uvp, blpairs, spw, pol, component='real', average_blpairs=Fa if new_plot: return fig -def delay_wedge(uvps, blpairs, fold=False, cosmo=True, center_line=True, horizon_lines=True, cosmo_cbar=True, suptitle=''): +def delay_wedge(uvp, blpairs, spw, pol, fold=False, delay=False, + center_line=True, horizon_lines=True, suptitle='', subtitle=None, ax=None): + """ Plot a 2D delay spectrum (or spectra) for a group of baselines. Parameters ---------- - uvps : UVPspec - List of UVPSpec objects, containing delay spectra for a set of baseline-pairs, - times, polarizations, and spectral windows. + uvp : UVPspec + UVPSpec object or list of UVPSpec objects which will be plotted side + by side, containing delay spectra for a set of baseline-pairs, times, + polarizations, and spectral windows. blpairs : list of tuples List of baseline-pair tuples. + + spw : int or str, or list of ints or strs + Which spectral window to plot for each uvp object to be plotted. If + plotting multiple UVPSpec objects, this should be a list of the same + length as the number of objects to be plotted. + + pol : str or list of strs + Which pol (or pols) to plot. If plotting multiple UVPSpec objects, + this should be a list of of the same length as the number of + objects to be plotted. + get_reds_from : + fold : bool, optional - Whether to fold the power spectrum in :math:`|k_\parallel|`. + Whether to fold the power spectrum in k_parallel. Default: False. - cosmo : bool, optional - Whether to convert axes to cosmological units :math:`|k_\parallel|` and - :math:`|k_\perpendicular|`. - Default: True. + delay : bool, optional + Whether to plot the axes in tau. If False, axes will be plotted in k. + Default: False. center_line : bool, optional - Whether to plot a dotted line at :math:`|k_\parallel|` =0. + Whether to plot a dotted line at k_parallel=0. Default: True. horizon_lines : bool, optional Whether to plot dotted lines along the horizon. Default: True. - cosmo_cbar : bool, optional - Whether power has been converted to appropriate units through calibration. - Default: True. - suptitle : string, optional Suptitle for the plot. If not provided, no suptitle will be plotted. Default: ''. + + subtitle : string, list, optional + Title for each subplot. If plotting multiple UVPSpec objects, this + should be a list of strings of the same length as the number of + objects to be plotted. If not provided, no subtitles will be plotted. + Default: None. + + ax : matplotlib.axes, optional + Use this to pass in an existing Axes object, which the power spectra + will be added to. (Warning: Labels and legends will not be altered in + this case, even if the existing plot has completely different axis + labels etc.) If None, a new Axes object will be created. + Default: None. Returns ------- - f : matplotlib.pyplot.Figure + fig : matplotlib.pyplot.Figure Matplotlib Figure instance. """ - - if (type(uvps) != list) and (type(uvps) != np.ndarray): - raise AttributeError("The uvps paramater should be a list of UVPSpec objects.") - - #Initialize axes objects - ncols = len(uvps) - f, axes = plt.subplots( - ncols=ncols, - nrows=1, - sharex=True, - sharey=True, - squeeze=False, - figsize=(20, 10)) + + #Make sure uvp is a UVPSpec or list of UVPSpec + if isinstance(uvp, uvpspec.UVPSpec): + uvp = [uvp] + if isinstance(uvp, list): + for uv in uvp: + if not isinstance(uv, uvpspec.UVPSpec): + raise AttributeError('The uvp parameter should be a UVPSpec ' + 'object or list of UVPSpec objects.') + else: + raise AttributeError('The uvp parameter should be a UVPSpec ' + 'object or list of UVPSpec objects.') + + #check to make sure spw, pol, and subtitle are appropriately formatted + if isinstance(spw, (int,str)): + spw = [spw] + if not isinstance(spw, list) or not len(spw)==len(uvp): + raise AttributeError('The same number of arguments should be provided ' + 'for spw and uvp.') + if isinstance(pol, str): + pol = [pol] + if not isinstance(pol, list) or not len(pol)==len(uvp): + raise AttributeError('The same number of arguments should be provided ' + 'for pol and uvp.') + if subtitle is not None: + if isinstance(subtitle, str): + subtitle = [subtitle] + if not isinstance(subtitle, list) or not len(subtitle)==len(uvp): + raise AttributeError('The same number of arguments should be ' + 'provided for subtitle and uvp.') + else: + subtitle = ['' for i in range(len(uvp))] + + + #Create new Axes if none specified + new_plot = False + if ax is None: + ncols = len(uvp) + fig, axes = plt.subplots( + ncols=ncols, + nrows=1, + sharex=True, + sharey=True, + squeeze=False, + figsize=(20, 10)) + new_plot = True #Plot each uvp pols = [] plt.subplots_adjust(wspace=0, hspace=0.1) - for uvp, ax in zip(uvps, axes.flatten()): + for uvp, ax, spw, pol, subtitle in zip(uvp, axes.flatten(), spw, pol, subtitle): #Find redundant-length baseline groups and their baseline lengths BIN_WIDTH = 0.3 #redundancy tolerance in meters NORM_BINS = np.arange(0.0, 10000.0, BIN_WIDTH) @@ -576,22 +632,20 @@ def delay_wedge(uvps, blpairs, fold=False, cosmo=True, center_line=True, horizon bllens = sorted(bllens) #Average the spectra along redundant baseline groups and time - uvp.average_spectra(blpair_groups=sorted_blpairs.values(), time_avg=True) - - #Grab polarization string for naming the plot later on - pol_int_to_str = {1: 'pI', 2: 'pQ', 3: 'pU', 4: 'pV', -5: 'XX', -6: 'YY', -7: 'XY', -8: 'YX'} - pol = pol_int_to_str[uvp.pol_array[0]] - pols.append(pol) + uvp = uvp.average_spectra(blpair_groups=sorted_blpairs.values(), time_avg=True, inplace=False) #Format data array - data = uvp.data_array[0][:, :, 0] if fold: uvp.fold_spectra() - data = uvp.data_array[0][:, data.shape[1] // 2:, 0] - + data = [] + for blpair in uvp.get_blpairs(): + data.append(uvp.get_data((spw, blpair, pol))[0]) + data = np.array(data) + #Format k_para axis - x_axis = (uvp.get_dlys(0)*1e9).tolist() - if cosmo: + if delay: + x_axis = (uvp.get_dlys(0)*1e9).tolist() + else: x_axis = uvp.get_kparas(0).tolist() if fold: x_axis.insert(0, 0) @@ -617,7 +671,13 @@ def delay_wedge(uvps, blpairs, fold=False, cosmo=True, center_line=True, horizon #Setting y axis ticks and labels ax.set_yticks(bllen_midindices) - ax.set_yticklabels(bllens, fontsize=10) + if delay: + ax.set_yticklabels(bllens, fontsize=10) + else: + avg_z = uvp.cosmo.f2z(np.mean(uvp.freq_array[uvp.spw_to_freq_indices(spw)])) + perpfactor = uvp.cosmo.bl_to_kperp(avg_z, little_h=True) + bl_lens = [np.round(bllen * perpfactor, 3) for bllen in bllens] + ax.set_yticklabels(bl_lens, fontsize=10) #Plot the data im = ax.imshow( @@ -629,8 +689,10 @@ def delay_wedge(uvps, blpairs, fold=False, cosmo=True, center_line=True, horizon #Plot the horizon lines if requested if horizon_lines: horizons = [(bllen*u.m/c.c).to(u.ns).value for bllen in bllens] - #if cosmo: - #XXX still need to convert horizons to appropriate kpara values + if not delay: + avg_z = uvp.cosmo.f2z(np.mean(uvp.freq_array[uvp.spw_to_freq_indices(spw)])) + parafactor = uvp.cosmo.tau_to_kpara(avg_z, little_h=True) + horizons = [horizon * parafactor * 10**-9 for horizon in horizons] j = 0 for i, (horizon, bllen) in enumerate(zip(horizons, bllens)): x1, y1 = [horizon, horizon], [j, bllen_indices[i]] @@ -643,29 +705,35 @@ def delay_wedge(uvps, blpairs, fold=False, cosmo=True, center_line=True, horizon ax.axvline(x=0, color='#000000', ls='--', lw=1) ax.tick_params(axis='both', direction='inout') - ax.set_title(pol, fontsize=15) + ax.set_title(subtitle, fontsize=15) #add colorbar - cbar_ax = f.add_axes([0.9125, 0.25, 0.025, 0.5]) - cbar = f.colorbar(im, cax=cbar_ax) + if new_plot: + cbar_ax = fig.add_axes([0.9125, 0.25, 0.025, 0.5]) + cbar = fig.colorbar(im, cax=cbar_ax) + + #add units + psunits = uvp.units + if "h^-1" in psunits: psunits = psunits.replace("h^-1", "h^{-1}") + if "h^-3" in psunits: psunits = psunits.replace("h^-3", "h^{-3}") + if "Mpc" in psunits and "\\rm" not in psunits: + psunits = psunits.replace("Mpc", r"{\rm Mpc}") + if "pi" in psunits and "\\pi" not in psunits: + psunits = psunits.replace("pi", r"\pi") + cbar.set_label("P [$%s$]" % psunits, fontsize=15, ha='center') #add axis labels - if cosmo: - f.text(0.5, 0.05, r'$k_{\parallel}\ [h\ Mpc^{-1}]$', fontsize=15, ha='center') - f.text(0.07, 0.5, r'$k_{\perp}\ [\rm h\ Mpc^{-1}]$', fontsize=15, va='center', rotation='vertical') - else: - f.text(0.5, 0.05, r'$\tau\ [ns]$', fontsize=15, ha='center') - f.text(0.07, 0.5, r'Baseline Group $[m]$', fontsize=15, va='center', rotation='vertical') - - #add colorbar label - if cosmo_cbar: - cbar.set_label(r"$P(k_{\parallel})\ [mK^2h^{-3}Mpc^3]$", fontsize=15, ha='center') - else: - cbar.set_label(r'$P(\tau)\ [mK^2]$', fontsize=15, ha='center') + if new_plot and delay: + fig.text(0.5, 0.05, r'$\tau\ [ns]$', fontsize=15, ha='center') + fig.text(0.07, 0.5, r'Baseline Group $[m]$', fontsize=15, va='center', rotation='vertical') + elif new_plot: + fig.text(0.5, 0.05, r'$k_{\parallel}\ [h\ Mpc^{-1}]$', fontsize=15, ha='center') + fig.text(0.07, 0.5, r'$k_{\perp}\ [\rm h\ Mpc^{-1}]$', fontsize=15, va='center', rotation='vertical') #Add suptitle if requested - if len(suptitle) > 0: - f.suptitle(suptitle, fontsize=15) + if len(suptitle) > 0 and new_plot: + fig.suptitle(suptitle, fontsize=15) #Return Figure: the axis is an attribute of figure - return f \ No newline at end of file + if new_plot: + return fig \ No newline at end of file From 19d32a3be047aae39ff5673f74319983c7e245e8 Mon Sep 17 00:00:00 2001 From: Nicholas Kern Date: Thu, 6 Sep 2018 21:17:39 +0800 Subject: [PATCH 3/9] minor edit to utils.get_blvec_reds and plot.delay_wedge --- hera_pspec/plot.py | 70 ++++++++++++++++++++------------------------- hera_pspec/utils.py | 8 +++--- 2 files changed, 35 insertions(+), 43 deletions(-) diff --git a/hera_pspec/plot.py b/hera_pspec/plot.py index f6a669e3..5c89e7a6 100644 --- a/hera_pspec/plot.py +++ b/hera_pspec/plot.py @@ -495,51 +495,45 @@ def delay_waterfall(uvp, blpairs, spw, pol, component='real', average_blpairs=Fa # Return Axes if new_plot: return fig - -def delay_wedge(uvp, blpairs, spw, pol, fold=False, delay=False, - center_line=True, horizon_lines=True, suptitle='', subtitle=None, ax=None): + +def delay_wedge(uvp, spw, pol, blpairs=None, fold=False, delay=True, + center_line=True, horizon_lines=True, suptitle='', + subtitle=None, ax=None): """ Plot a 2D delay spectrum (or spectra) for a group of baselines. Parameters ---------- - uvp : UVPspec - UVPSpec object or list of UVPSpec objects which will be plotted side - by side, containing delay spectra for a set of baseline-pairs, times, - polarizations, and spectral windows. - - blpairs : list of tuples - List of baseline-pair tuples. - - spw : int or str, or list of ints or strs - Which spectral window to plot for each uvp object to be plotted. If - plotting multiple UVPSpec objects, this should be a list of the same - length as the number of objects to be plotted. - - pol : str or list of strs - Which pol (or pols) to plot. If plotting multiple UVPSpec objects, - this should be a list of of the same length as the number of - objects to be plotted. - - get_reds_from : + uvp : UVPSpec + UVPSpec object containing delay spectra to plot. + + spw : integer + Which spectral window to plot. + + pol : int or str + Polarization integer or string t + + blpairs : list of tuples, optional + List of baseline-pair tuples to use in plotting. fold : bool, optional Whether to fold the power spectrum in k_parallel. Default: False. - + delay : bool, optional - Whether to plot the axes in tau. If False, axes will be plotted in k. - Default: False. - + Whether to plot the axes in tau (ns). If False, axes will + be plotted in cosmological units. + Default: True. + center_line : bool, optional Whether to plot a dotted line at k_parallel=0. Default: True. - + horizon_lines : bool, optional Whether to plot dotted lines along the horizon. Default: True. - + suptitle : string, optional Suptitle for the plot. If not provided, no suptitle will be plotted. Default: ''. @@ -556,14 +550,13 @@ def delay_wedge(uvp, blpairs, spw, pol, fold=False, delay=False, this case, even if the existing plot has completely different axis labels etc.) If None, a new Axes object will be created. Default: None. - + Returns ------- fig : matplotlib.pyplot.Figure Matplotlib Figure instance. """ - - #Make sure uvp is a UVPSpec or list of UVPSpec + # Make sure uvp is a UVPSpec or list of UVPSpec if isinstance(uvp, uvpspec.UVPSpec): uvp = [uvp] if isinstance(uvp, list): @@ -574,8 +567,8 @@ def delay_wedge(uvp, blpairs, spw, pol, fold=False, delay=False, else: raise AttributeError('The uvp parameter should be a UVPSpec ' 'object or list of UVPSpec objects.') - - #check to make sure spw, pol, and subtitle are appropriately formatted + + # Check to make sure spw, pol, and subtitle are appropriately formatted if isinstance(spw, (int,str)): spw = [spw] if not isinstance(spw, list) or not len(spw)==len(uvp): @@ -594,9 +587,8 @@ def delay_wedge(uvp, blpairs, spw, pol, fold=False, delay=False, 'provided for subtitle and uvp.') else: subtitle = ['' for i in range(len(uvp))] - - #Create new Axes if none specified + # Create new Axes if none specified new_plot = False if ax is None: ncols = len(uvp) @@ -608,12 +600,12 @@ def delay_wedge(uvp, blpairs, spw, pol, fold=False, delay=False, squeeze=False, figsize=(20, 10)) new_plot = True - - #Plot each uvp + + # Plot each uvp pols = [] plt.subplots_adjust(wspace=0, hspace=0.1) for uvp, ax, spw, pol, subtitle in zip(uvp, axes.flatten(), spw, pol, subtitle): - #Find redundant-length baseline groups and their baseline lengths + # Find redundant-length baseline groups and their baseline lengths BIN_WIDTH = 0.3 #redundancy tolerance in meters NORM_BINS = np.arange(0.0, 10000.0, BIN_WIDTH) sorted_blpairs = {} @@ -631,7 +623,7 @@ def delay_wedge(uvp, blpairs, spw, pol, fold=False, delay=False, sorted_blpairs[bllen] = [antnums] bllens = sorted(bllens) - #Average the spectra along redundant baseline groups and time + # Average the spectra along redundant baseline groups and time uvp = uvp.average_spectra(blpair_groups=sorted_blpairs.values(), time_avg=True, inplace=False) #Format data array diff --git a/hera_pspec/utils.py b/hera_pspec/utils.py index 84a933d0..f092b230 100644 --- a/hera_pspec/utils.py +++ b/hera_pspec/utils.py @@ -701,15 +701,15 @@ def config_pspec_blpairs(uv_templates, pol_pairs, group_pairs, exclude_auto_bls= def get_blvec_reds(blvecs, bl_error_tol=1.0): """ - Given a blvecs dictionary, form groups of baseline-pair objects based on + Given a blvecs dictionary, form groups of baseline-pairs based on redundancy in ENU coordinates. Note: this only uses the East-North components of the baseline vectors to calculate redundancy. Parameters: ----------- blvecs : dictionary (or UVPSpec object) - A dictionary with baseline vectors as values. Alternatively, this - can be a UVPSpec object with baseline-pairs and baseline vectors. + A dictionary with len-2 or 3 ndarray baseline vectors as values. + Alternatively, this can be a UVPSpec object. bl_error_tol : int Redundancy tolerance of baseline vector in meters. @@ -764,7 +764,7 @@ def get_blvec_reds(blvecs, bl_error_tol=1.0): bl_tag = "{:03.0f}_{:03.0f}".format(bl_len, bl_ang) # append to list if unique within tolerance - match = [np.all(np.isclose(blv, bl_vec, bl_error_tol)) for blv in red_bl_vec] + match = [np.all(np.isclose(blv, bl_vec, atol=bl_error_tol)) for blv in red_bl_vec] if np.any(match): match_id = np.where(match)[0][0] red_bl_grp[match_id].append(bl) From 3f1f41dccc4a1a73ec81319e4c321f57f348ca07 Mon Sep 17 00:00:00 2001 From: Nicholas Kern Date: Mon, 10 Sep 2018 16:54:38 +0800 Subject: [PATCH 4/9] plot.delay_wedge consistency overhaul + bug fixes --- hera_pspec/plot.py | 441 ++++++++++++++++++++++++++------------------- 1 file changed, 252 insertions(+), 189 deletions(-) diff --git a/hera_pspec/plot.py b/hera_pspec/plot.py index ad33784d..a8a5d4fb 100644 --- a/hera_pspec/plot.py +++ b/hera_pspec/plot.py @@ -214,9 +214,11 @@ def delay_spectrum(uvp, blpairs, spw, pol, average_blpairs=False, cax = None if lines: if logscale: - cax, = ax.plot(x, np.abs(y), marker='None', label=label, **kwargs) + _y = np.abs(y) else: - cax, = ax.plot(x, y, marker='None', label=label, **kwargs) + _y = y + + cax, = ax.plot(x, _y, marker='None', label=label, **kwargs) if markers: if cax is None: @@ -245,7 +247,11 @@ def delay_spectrum(uvp, blpairs, spw, pol, average_blpairs=False, c = None else: c = cax.get_color() - cax = ax.errorbar(x, np.abs(y), fmt='none', ecolor=c, yerr=cast(errs[i]), **kwargs) + if logscale: + _y = np.abs(yp) + else: + _y = y + cax = ax.errorbar(x, _y, fmt='none', ecolor=c, yerr=cast(errs[i]), **kwargs) else: raise KeyError("Error variable '%s' not found in stats_array of UVPSpec object." % error) @@ -377,7 +383,7 @@ def delay_waterfall(uvp, blpairs, spw, pol, component='real', average_blpairs=Fa blpairs.append([blpairs_in[i],]) else: blpairs.append(blpairs_in[i]) - + # iterate through and make sure they are blpair integers _blpairs = [] for blpgrp in blpairs: @@ -582,11 +588,16 @@ def delay_waterfall(uvp, blpairs, spw, pol, component='real', average_blpairs=Fa return fig -def delay_wedge(uvp, spw, pol, blpairs=None, fold=False, delay=True, - center_line=True, horizon_lines=True, suptitle='', - subtitle=None, ax=None): +def delay_wedge(uvp, spw, pol, blpairs=None, times=None, fold=False, delay=True, + rotate=False, component='real', log10=True, loglog=False, + red_tol=1.0, center_line=False, horizon_lines=False, + title=None, ax=None, cmap='viridis', figsize=(8, 6), + deltasq=False, colorbar=False, cbax=None, vmin=None, vmax=None, + edgecolor='none', flip_xax=False, flip_yax=False, lw=2, **kwargs): """ - Plot a 2D delay spectrum (or spectra) for a group of baselines. + Plot a 2D delay spectrum (or spectra) from a UVPSpec object. Note that + all integrations and redundant baselines are averaged (unless specifying times) + before plotting. Parameters ---------- @@ -597,11 +608,15 @@ def delay_wedge(uvp, spw, pol, blpairs=None, fold=False, delay=True, Which spectral window to plot. pol : int or str - Polarization integer or string t + Polarization integer or string blpairs : list of tuples, optional List of baseline-pair tuples to use in plotting. + times : list, optional + An ndarray or list of times from uvp.time_avg_array to + select on before plotting. Default: None. + fold : bool, optional Whether to fold the power spectrum in k_parallel. Default: False. @@ -611,207 +626,255 @@ def delay_wedge(uvp, spw, pol, blpairs=None, fold=False, delay=True, be plotted in cosmological units. Default: True. + rotate : bool, optional + If False, use baseline-type as x-axis and delay as y-axis, + else use baseline-type as y-axis and delay as x-axis. + Default: False + + component : str, optional + Component of complex spectra to plot. Options=['real', 'imag', 'abs'] + Default: 'real'. + + log10 : bool, optional + If True, take log10 of data before plotting. Default: True + + loglog : bool, optional + If True, turn x-axis and y-axis into log-log scale. Default: False + + red_tol : float, optional + Redundancy tolerance when solving for redundant groups in meters. + Default: 1.0 + center_line : bool, optional - Whether to plot a dotted line at k_parallel=0. - Default: True. + Whether to plot a dotted line at k_perp = 0. + Default: False. horizon_lines : bool, optional Whether to plot dotted lines along the horizon. - Default: True. - - suptitle : string, optional - Suptitle for the plot. If not provided, no suptitle will be plotted. - Default: ''. + Default: False. - subtitle : string, list, optional - Title for each subplot. If plotting multiple UVPSpec objects, this - should be a list of strings of the same length as the number of - objects to be plotted. If not provided, no subtitles will be plotted. - Default: None. + title : string, optional + Title for subplot. Default: None. ax : matplotlib.axes, optional - Use this to pass in an existing Axes object, which the power spectra - will be added to. (Warning: Labels and legends will not be altered in - this case, even if the existing plot has completely different axis - labels etc.) If None, a new Axes object will be created. - Default: None. + If not None, use this axes as a subplot for delay wedge. + + cmap : str, optional + Colormap of wedge plot. Default: 'viridis' + + figsize : len-2 integer tuple, optional + If ax is None, this is the new figure size. + + deltasq : bool, optional + Convert to Delta^2 before plotting. Default: False + + colorbar : bool, optional + Add a colorbar to the plot. Default: False + + cbax : matplotlib.axes, optional + Axis object for adding colorbar if True. Default: None + + vmin : float, optional + Minimum range of colorscale. Default: None + + vmax : float, optional + Maximum range of colorscale. Default: None + + edgecolor : str, optional + Edgecolor of bins in pcolormesh. Default: 'none' + + flip_xax : bool, optional + Flip xaxis if True. Default: False + + flip_yax : bool, optional + Flip yaxis if True. Default: False + + lw : int, optional + Line-width of horizon and center lines if plotted. Default: 2. + + kwargs : dictionary + Additional keyword arguments to pass to pcolormesh() call. Returns ------- fig : matplotlib.pyplot.Figure - Matplotlib Figure instance. + Matplotlib Figure instance if ax is None. """ - # Make sure uvp is a UVPSpec or list of UVPSpec - if isinstance(uvp, uvpspec.UVPSpec): - uvp = [uvp] - if isinstance(uvp, list): - for uv in uvp: - if not isinstance(uv, uvpspec.UVPSpec): - raise AttributeError('The uvp parameter should be a UVPSpec ' - 'object or list of UVPSpec objects.') - else: - raise AttributeError('The uvp parameter should be a UVPSpec ' - 'object or list of UVPSpec objects.') - - # Check to make sure spw, pol, and subtitle are appropriately formatted - if isinstance(spw, (int,str)): - spw = [spw] - if not isinstance(spw, list) or not len(spw)==len(uvp): - raise AttributeError('The same number of arguments should be provided ' - 'for spw and uvp.') - if isinstance(pol, str): - pol = [pol] - if not isinstance(pol, list) or not len(pol)==len(uvp): - raise AttributeError('The same number of arguments should be provided ' - 'for pol and uvp.') - if subtitle is not None: - if isinstance(subtitle, str): - subtitle = [subtitle] - if not isinstance(subtitle, list) or not len(subtitle)==len(uvp): - raise AttributeError('The same number of arguments should be ' - 'provided for subtitle and uvp.') - else: - subtitle = ['' for i in range(len(uvp))] + # type checking + uvp = copy.deepcopy(uvp) + assert isinstance(uvp, uvpspec.UVPSpec), "input uvp must be a UVPSpec object" + assert isinstance(spw, (int, np.integer)) + assert isinstance(pol, (int, str, np.integer, np.str)) - # Create new Axes if none specified + # Create new ax if none specified new_plot = False if ax is None: - ncols = len(uvp) - fig, axes = plt.subplots( - ncols=ncols, - nrows=1, - sharex=True, - sharey=True, - squeeze=False, - figsize=(20, 10)) + fig, ax = plt.subplots(figsize=figsize) new_plot = True + else: + fig = ax.get_figure() - # Plot each uvp - pols = [] - plt.subplots_adjust(wspace=0, hspace=0.1) - for uvp, ax, spw, pol, subtitle in zip(uvp, axes.flatten(), spw, pol, subtitle): - # Find redundant-length baseline groups and their baseline lengths - BIN_WIDTH = 0.3 #redundancy tolerance in meters - NORM_BINS = np.arange(0.0, 10000.0, BIN_WIDTH) - sorted_blpairs = {} - bllens = [] - blpair_seps = uvp.get_blpair_seps() - blpair_array = uvp.blpair_array - for antnums in blpairs: - blpair = uvp.antnums_to_blpair(antnums) - bllen = blpair_seps[np.where(blpair_array==blpair)[0][0]] - bllen = np.round(np.digitize(bllen, NORM_BINS) * BIN_WIDTH, 1) - if bllen in bllens: - sorted_blpairs[bllen].append(antnums) - else: - bllens.append(bllen) - sorted_blpairs[bllen] = [antnums] - bllens = sorted(bllens) - - # Average the spectra along redundant baseline groups and time - uvp = uvp.average_spectra(blpair_groups=sorted_blpairs.values(), time_avg=True, inplace=False) - - #Format data array - if fold: - uvp.fold_spectra() - data = [] - for blpair in uvp.get_blpairs(): - data.append(uvp.get_data((spw, blpair, pol))[0]) - data = np.array(data) - - #Format k_para axis + # Select out times if provided + if times is not None: + uvp.select(times=times, inplace=True) + + # Average across redundant groups and time + # this also ensures blpairs are ordered from short_bl --> long_bl + blp_grps, lens, angs, tags = utils.get_blvec_reds(uvp, bl_error_tol=red_tol) + uvp.average_spectra(blpair_groups=blp_grps, time_avg=True, inplace=True) + + # Convert to DeltaSq + if deltasq: + uvp.convert_to_deltasq(inplace=True) + + # Fold array + if fold: + uvp.fold_spectra() + + # Format ticks + if delay: + x_axis = uvp.get_dlys(spw) * 1e9 + y_axis = uvp.get_blpair_seps() + else: + x_axis = uvp.get_kparas(spw, little_h=little_h) + y_axis = uvp.get_kperps(spw, little_h=little_h) + if rotate: + _x_axis = y_axis + y_axis = x_axis + x_axis = _x_axis + + # Conigure Units + psunits = "({})^2\ {}".format(uvp.vis_units, uvp.norm_units) + if "h^-1" in psunits: psunits = psunits.replace("h^-1", "h^{-1}\ ") + if "h^-3" in psunits: psunits = psunits.replace("h^-3", "h^{-3}\ ") + if "Hz" in psunits: psunits = psunits.replace("Hz", r"{\rm Hz}\ ") + if "str" in psunits: psunits = psunits.replace("str", r"\,{\rm str}\ ") + if "Mpc" in psunits and "\\rm" not in psunits: + psunits = psunits.replace("Mpc", r"{\rm Mpc}") + if "pi" in psunits and "\\pi" not in psunits: + psunits = psunits.replace("pi", r"\pi") + if "beam normalization not specified" in psunits: + psunits = psunits.replace("beam normalization not specified", + r"{\rm unnormed}") + + # get data casting + if component == 'real': + cast = np.real + elif component == 'imag': + cast = np.imag + elif component == 'abs': + cast = np.abs + else: + raise ValueError("Did not understand component {}".format(component)) + + # get data with shape (Nblpairs, Ndlys) + data = cast([uvp.get_data((spw, blp, pol)).squeeze() for blp in uvp.get_blpairs()]) + + # take log10 + if log10: + data = np.log10(np.abs(data)) + + # rotate + if rotate: + data = np.rot90(data[:, ::-1], k=1) + + # Get bin edges + xdiff = np.diff(x_axis) + x_edges = np.array([x_axis[0]-xdiff[0]/2.0] + list(x_axis[:-1]+xdiff/2.0) + [x_axis[-1]+xdiff[-1]/2.0]) + ydiff = np.diff(y_axis) + y_edges = np.array([y_axis[0]-ydiff[0]/2.0] + list(y_axis[:-1]+ydiff/2.0) + [y_axis[-1]+ydiff[-1]/2.0]) + X, Y = np.meshgrid(x_edges, y_edges) + + # plot + cax = ax.pcolormesh(X, Y, data, cmap=cmap, edgecolor=edgecolor, lw=lw, + vmin=vmin, vmax=vmax, **kwargs) + + # loglog + if loglog: + ax.set_xscale('log') + ax.set_yscale('log') + + # Configure ticks + if rotate: + ax.set_xticks(np.around(x_axis, 4)) + else: + ax.set_yticks(np.around(y_axis, 4)) + + # Add colorbar + if colorbar: + if cbax is None: + cbax = ax + cbar = fig.colorbar(cax, ax=cbax) if delay: - x_axis = (uvp.get_dlys(0)*1e9).tolist() + p = "P({},\ {})".format(r'\tau', r'|\vec{b}|') else: - x_axis = uvp.get_kparas(0).tolist() - if fold: - x_axis.insert(0, 0) - ax.set_xlim((x_axis[0] / 2, x_axis[-1] / 2)) - - #Format k_perp axis - #Calculate kpr indices, and stretch wedgeslices to length of baseline norms - bllen_indices = [int(bllen * 10) for bllen in bllens] - stretched_data = np.zeros((bllen_indices[-1], data.shape[-1]), dtype=np.float64) - j = 0 - for i in range(len(bllen_indices)): - stretched_data[j:bllen_indices[i]] = data[i] - j = bllen_indices[i] - data = stretched_data[...] - - #Find kpr mid-indicies for the tickmarks - bllen_midindices = [] - for i in range(len(bllen_indices)): - if i == 0: - bllen_midindices.append(bllen_indices[i] / 2.) - else: - bllen_midindices.append((bllen_indices[i] + bllen_indices[i - 1]) / 2.) - - #Setting y axis ticks and labels - ax.set_yticks(bllen_midindices) - if delay: - ax.set_yticklabels(bllens, fontsize=10) + p = "P({},\ {})".format(r'k_\parallel', r'k_\perp') + if log10: + psunits = r"$\log_{{10}}\ {}\ [{}]$".format(p, psunits) + else: + psunits = r"${}\ [{}]$".format(p, psunits) + cbar.set_label(psunits, fontsize=14) + + # Configure tick labels + if delay: + xlabel = r"$\tau$ $[{\rm ns}]$" + ylabel = r"$|\vec{b}|$ $[{\rm m}]$" + else: + xlabel = r"$k_{\parallel}\ [h\ \rm Mpc^{-1}]$" + ylabel = r"$k_{\perp}\ [h\ \rm Mpc^{-1}]$" + if rotate: + _xlabel = ylabel + ylabel = xlabel + xlabel = _xlabel + if ax.get_xlabel() == '': + ax.set_xlabel(xlabel, fontsize=16) + if ax.get_ylabel() == '': + ax.set_ylabel(ylabel, fontsize=16) + + # Configure center line + if center_line: + if rotate: + ax.axhline(y=0, color='#000000', ls='--', lw=lw) else: + ax.axvline(x=0, color='#000000', ls='--', lw=lw) + + # Plot horizons + if horizon_lines: + # get horizon in ns + horizons = uvp.get_blpair_seps() / conversions.units.c * 1e9 + + # convert to cosmological wave vector + if not delay: + # Get average redshift of spw avg_z = uvp.cosmo.f2z(np.mean(uvp.freq_array[uvp.spw_to_freq_indices(spw)])) - perpfactor = uvp.cosmo.bl_to_kperp(avg_z, little_h=True) - bl_lens = [np.round(bllen * perpfactor, 3) for bllen in bllens] - ax.set_yticklabels(bl_lens, fontsize=10) - - #Plot the data - im = ax.imshow( - np.log10(np.abs(data)), - aspect='auto', - interpolation='nearest', - extent=[x_axis[0], x_axis[-1], bllen_indices[-1], 0]) - - #Plot the horizon lines if requested - if horizon_lines: - horizons = [(bllen*u.m/c.c).to(u.ns).value for bllen in bllens] - if not delay: - avg_z = uvp.cosmo.f2z(np.mean(uvp.freq_array[uvp.spw_to_freq_indices(spw)])) - parafactor = uvp.cosmo.tau_to_kpara(avg_z, little_h=True) - horizons = [horizon * parafactor * 10**-9 for horizon in horizons] - j = 0 - for i, (horizon, bllen) in enumerate(zip(horizons, bllens)): - x1, y1 = [horizon, horizon], [j, bllen_indices[i]] - x2, y2 = [-horizon, -horizon], [j, bllen_indices[i]] - ax.plot(x1, y1, x2, y2, color='#ffffff', linestyle='--', linewidth=1) - j = bllen_indices[i] - - #plot center line at k_para=0 if requested - if center_line: - ax.axvline(x=0, color='#000000', ls='--', lw=1) + horizons *= uvp.cosmo.tau_to_kpara(avg_z, little_h=little_h) - ax.tick_params(axis='both', direction='inout') - ax.set_title(subtitle, fontsize=15) - - #add colorbar - if new_plot: - cbar_ax = fig.add_axes([0.9125, 0.25, 0.025, 0.5]) - cbar = fig.colorbar(im, cax=cbar_ax) - - #add units - psunits = uvp.units - if "h^-1" in psunits: psunits = psunits.replace("h^-1", "h^{-1}") - if "h^-3" in psunits: psunits = psunits.replace("h^-3", "h^{-3}") - if "Mpc" in psunits and "\\rm" not in psunits: - psunits = psunits.replace("Mpc", r"{\rm Mpc}") - if "pi" in psunits and "\\pi" not in psunits: - psunits = psunits.replace("pi", r"\pi") - cbar.set_label("P [$%s$]" % psunits, fontsize=15, ha='center') - - #add axis labels - if new_plot and delay: - fig.text(0.5, 0.05, r'$\tau\ [ns]$', fontsize=15, ha='center') - fig.text(0.07, 0.5, r'Baseline Group $[m]$', fontsize=15, va='center', rotation='vertical') - elif new_plot: - fig.text(0.5, 0.05, r'$k_{\parallel}\ [h\ Mpc^{-1}]$', fontsize=15, ha='center') - fig.text(0.07, 0.5, r'$k_{\perp}\ [\rm h\ Mpc^{-1}]$', fontsize=15, va='center', rotation='vertical') - - #Add suptitle if requested - if len(suptitle) > 0 and new_plot: - fig.suptitle(suptitle, fontsize=15) - - #Return Figure: the axis is an attribute of figure + # iterate over bins and plot lines + if rotate: + bin_edges = x_edges + else: + bin_edges = y_edges + for i, hor in enumerate(horizons): + if rotate: + ax.plot(bin_edges[i:i+2], [hor, hor], color='#ffffff', ls='--', lw=lw) + if not uvp.folded: + ax.plot(bin_edges[i:i+2], [-hor, -hor], color='#ffffff', ls='--', lw=lw) + else: + ax.plot([hor, hor], bin_edges[i:i+2], color='#ffffff', ls='--', lw=lw) + if not uvp.folded: + ax.plot([-hor, -hor], bin_edges[i:i+2], color='#ffffff', ls='--', lw=lw) + + # flip axes + if flip_xax: + plt.gca().invert_xaxis() + if flip_yax: + plt.gca().invert_yaxis() + + # add title + if title is not None: + ax.set_title(title, fontsize=12) + + # return figure if new_plot: return fig @@ -829,7 +892,7 @@ def plot_uvdata_waterfalls(uvd, basename, data='data', plot_mode='log', Input data object. Waterfalls will be stored for all baselines and polarizations within the object; use uvd.select() to remove unwanted information. - + basename : str Base filename for the output plots. This must have two placeholders for baseline ID ('bl') and polarization ('pol'), From 94113c02474a2587a0608fa28b376f35b0934931 Mon Sep 17 00:00:00 2001 From: Nicholas Kern Date: Mon, 10 Sep 2018 21:56:32 +0800 Subject: [PATCH 5/9] updated plot.delay_wedge, added tests and added kwarg to utils.get_blvec_reds --- hera_pspec/plot.py | 21 ++++++++---- hera_pspec/tests/test_plot.py | 60 ++++++++++++++++++++++++++++++++++ hera_pspec/tests/test_utils.py | 6 ++++ hera_pspec/utils.py | 17 +++++++--- 4 files changed, 93 insertions(+), 11 deletions(-) diff --git a/hera_pspec/plot.py b/hera_pspec/plot.py index a8a5d4fb..249c1c9f 100644 --- a/hera_pspec/plot.py +++ b/hera_pspec/plot.py @@ -706,6 +706,9 @@ def delay_wedge(uvp, spw, pol, blpairs=None, times=None, fold=False, delay=True, assert isinstance(spw, (int, np.integer)) assert isinstance(pol, (int, str, np.integer, np.str)) + # check pspec units for little h + little_h = 'h^-3' in uvp.norm_units + # Create new ax if none specified new_plot = False if ax is None: @@ -720,7 +723,7 @@ def delay_wedge(uvp, spw, pol, blpairs=None, times=None, fold=False, delay=True, # Average across redundant groups and time # this also ensures blpairs are ordered from short_bl --> long_bl - blp_grps, lens, angs, tags = utils.get_blvec_reds(uvp, bl_error_tol=red_tol) + blp_grps, lens, angs, tags = utils.get_blvec_reds(uvp, bl_error_tol=red_tol, match_bl_lens=True) uvp.average_spectra(blpair_groups=blp_grps, time_avg=True, inplace=True) # Convert to DeltaSq @@ -774,6 +777,15 @@ def delay_wedge(uvp, spw, pol, blpairs=None, times=None, fold=False, delay=True, if log10: data = np.log10(np.abs(data)) + # loglog + if loglog: + ax.set_xscale('log') + ax.set_yscale('log') + ax.yaxis.set_major_formatter(matplotlib.ticker.LogFormatterSciNotation()) + ax.yaxis.set_minor_formatter(matplotlib.ticker.NullFormatter()) + ax.xaxis.set_major_formatter(matplotlib.ticker.LogFormatterSciNotation()) + ax.xaxis.set_minor_formatter(matplotlib.ticker.NullFormatter()) + # rotate if rotate: data = np.rot90(data[:, ::-1], k=1) @@ -786,14 +798,9 @@ def delay_wedge(uvp, spw, pol, blpairs=None, times=None, fold=False, delay=True, X, Y = np.meshgrid(x_edges, y_edges) # plot - cax = ax.pcolormesh(X, Y, data, cmap=cmap, edgecolor=edgecolor, lw=lw, + cax = ax.pcolormesh(X, Y, data, cmap=cmap, edgecolor=edgecolor, lw=0.01, vmin=vmin, vmax=vmax, **kwargs) - # loglog - if loglog: - ax.set_xscale('log') - ax.set_yscale('log') - # Configure ticks if rotate: ax.set_xticks(np.around(x_axis, 4)) diff --git a/hera_pspec/tests/test_plot.py b/hera_pspec/tests/test_plot.py index bd1e4060..8ac9bcf5 100644 --- a/hera_pspec/tests/test_plot.py +++ b/hera_pspec/tests/test_plot.py @@ -291,5 +291,65 @@ def test_uvdata_waterfalls(self): for f in figfiles: os.remove(f) + def test_delay_wedge(self): + """ Tests for plot.delay_wedge """ + # construct new uvp + reds, lens, angs = utils.get_reds(self.ds.dsets[0], pick_data_ants=True) + bls1, bls2, blps, _, _ = utils.calc_blpair_reds(self.ds.dsets[0], self.ds.dsets[1], + exclude_auto_bls=False, exclude_permutations=True) + uvp = self.ds.pspec(bls1, bls2, (0, 1), ('xx','xx'), + spw_ranges=[(300, 350)], + input_data_weight='identity', norm='I', + taper='blackman-harris', verbose=False) + + # test basic delay_wedge call + f1 = plot.delay_wedge(uvp, 0, 'xx', blpairs=None, times=None, fold=False, delay=True, + rotate=False, component='real', log10=False, loglog=False, red_tol=1.0, + center_line=False, horizon_lines=False, title=None, ax=None, cmap='viridis', + figsize=(8, 6), deltasq=False, colorbar=False, cbax=None, vmin=None, vmax=None, + edgecolor='none', flip_xax=False, flip_yax=False, lw=2) + plt.close() + + # specify blpairs and times + f2 = plot.delay_wedge(uvp, 0, 'xx', blpairs=uvp.get_blpairs()[-5:], times=uvp.time_avg_array[:1], + fold=False, delay=True, component='imag', + rotate=False, log10=False, loglog=False, red_tol=1.0, + center_line=False, horizon_lines=False, title=None, ax=None, cmap='viridis', + figsize=(8, 6), deltasq=False, colorbar=False, cbax=None, vmin=None, vmax=None, + edgecolor='none', flip_xax=False, flip_yax=False, lw=2) + plt.close() + + # fold, deltasq, cosmo and log10, loglog + f3 = plot.delay_wedge(uvp, 0, 'xx', blpairs=None, times=None, fold=True, delay=False, component='abs', + rotate=False, log10=True, loglog=True, red_tol=1.0, + center_line=False, horizon_lines=False, title='hello', ax=None, cmap='viridis', + figsize=(8, 6), deltasq=True, colorbar=False, cbax=None, vmin=None, vmax=None, + edgecolor='none', flip_xax=False, flip_yax=False, lw=2) + plt.close() + + # colorbar, vranges, flip_axes, edgecolors, lines + f4 = plot.delay_wedge(uvp, 0, 'xx', blpairs=None, times=None, fold=False, delay=False, component='abs', + rotate=False, log10=True, loglog=False, red_tol=1.0, + center_line=True, horizon_lines=True, title='hello', ax=None, cmap='viridis', + figsize=(8, 6), deltasq=False, colorbar=True, cbax=None, vmin=6, vmax=15, + edgecolor='grey', flip_xax=True, flip_yax=True, lw=2) + plt.close() + + # feed axes, red_tol + fig, ax = plt.subplots() + cbax = fig.add_axes([0.85, 0.1, 0.05, 0.9]) + cbax.axis('off') + plot.delay_wedge(uvp, 0, 'xx', blpairs=None, times=None, fold=False, delay=True, component='abs', + rotate=True, log10=True, loglog=False, red_tol=10.0, + center_line=False, horizon_lines=False, ax=ax, cmap='viridis', + figsize=(8, 6), deltasq=False, colorbar=True, cbax=cbax, vmin=None, vmax=None, + edgecolor='none', flip_xax=False, flip_yax=False, lw=2) + plt.close() + + + + + + if __name__ == "__main__": unittest.main() diff --git a/hera_pspec/tests/test_utils.py b/hera_pspec/tests/test_utils.py index 1312b1a2..f7d2a615 100644 --- a/hera_pspec/tests/test_utils.py +++ b/hera_pspec/tests/test_utils.py @@ -314,6 +314,12 @@ def test_get_blvec_reds(): red_bl_tag) = utils.get_blvec_reds(uvp, bl_error_tol=0.0) nt.assert_equal(len(red_bl_grp), uvp.Nblpairs) + # test combine angles + uvp = testing.uvpspec_from_data(fname, reds[:3], spw_ranges=[(10, 40)]) + (red_bl_grp, red_bl_len, red_bl_ang, + red_bl_tag) = utils.get_blvec_reds(uvp, bl_error_tol=1.0, match_bl_lens=True) + nt.assert_equal(len(red_bl_grp), 1) + def test_job_monitor(): # open empty files diff --git a/hera_pspec/utils.py b/hera_pspec/utils.py index 1c023132..37388495 100644 --- a/hera_pspec/utils.py +++ b/hera_pspec/utils.py @@ -680,7 +680,7 @@ def config_pspec_blpairs(uv_templates, pol_pairs, group_pairs, exclude_auto_bls= return groupings -def get_blvec_reds(blvecs, bl_error_tol=1.0): +def get_blvec_reds(blvecs, bl_error_tol=1.0, match_bl_lens=False): """ Given a blvecs dictionary, form groups of baseline-pairs based on redundancy in ENU coordinates. Note: this only uses the East-North components @@ -692,8 +692,12 @@ def get_blvec_reds(blvecs, bl_error_tol=1.0): A dictionary with len-2 or 3 ndarray baseline vectors as values. Alternatively, this can be a UVPSpec object. - bl_error_tol : int - Redundancy tolerance of baseline vector in meters. + bl_error_tol : float, optional + Redundancy tolerance of baseline vector in meters. Default: 1.0 + + match_bl_lens : bool, optional + Combine baseline groups of identical baseline length but + differing angle (using bl_error_tol). Default: False Returns: -------- @@ -745,7 +749,12 @@ def get_blvec_reds(blvecs, bl_error_tol=1.0): bl_tag = "{:03.0f}_{:03.0f}".format(bl_len, bl_ang) # append to list if unique within tolerance - match = [np.all(np.isclose(blv, bl_vec, atol=bl_error_tol)) for blv in red_bl_vec] + if match_bl_lens: + # match only on bl length + match = [np.all(np.isclose(bll, bl_len, atol=bl_error_tol)) for bll in red_bl_len] + else: + # match on full bl vector + match = [np.all(np.isclose(blv, bl_vec, atol=bl_error_tol)) for blv in red_bl_vec] if np.any(match): match_id = np.where(match)[0][0] red_bl_grp[match_id].append(bl) From d33192681b92bff9f9326f6eff588fd2f12a5791 Mon Sep 17 00:00:00 2001 From: Nicholas Kern Date: Mon, 10 Sep 2018 22:11:03 +0800 Subject: [PATCH 6/9] updated travis install order --- .travis.yml | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/.travis.yml b/.travis.yml index e6bd969b..95c01bfe 100644 --- a/.travis.yml +++ b/.travis.yml @@ -36,17 +36,16 @@ install: # create environment and install dependencies - conda create -q -n test-environment python=$TRAVIS_PYTHON_VERSION numpy scipy nose pip matplotlib coverage - source activate test-environment - - conda install -c conda-forge aipy + - conda install -c conda-forge healpy aipy - pip install coveralls + - pip install h5py - pip install git+https://github.com/HERA-Team/pyuvdata.git$PYUVDATA_VERSION - pip install git+https://github.com/HERA-Team/omnical.git - pip install git+https://github.com/HERA-Team/linsolve.git - pip install git+https://github.com/HERA-Team/hera_qm.git - pip install git+https://github.com/HERA-Team/uvtools.git - pip install git+https://github.com/HERA-Team/hera_cal.git - - pip install healpy - pip install scikit-learn - - pip install h5py - pip install pyyaml - pip install multiprocess - python setup.py install From f7c2d2b4534d266d5fa5825835012c0b111721b4 Mon Sep 17 00:00:00 2001 From: Nicholas Kern Date: Mon, 10 Sep 2018 22:52:27 +0800 Subject: [PATCH 7/9] bumbed pyuvdata install commit to HEAD 9/10/2018 --- .travis.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.travis.yml b/.travis.yml index 95c01bfe..b5080c9f 100644 --- a/.travis.yml +++ b/.travis.yml @@ -8,7 +8,7 @@ cache: pip: true env: - - PYUVDATA_VERSION="@9c94be0366eb41692f5501caac99981026a7ffe6" + - PYUVDATA_VERSION="@fae18f1a97e6c521b822b2ffc3c1afef4e2ac3b3" - PYUVDATA_VERSION="" allow_failures: From f51e76fda3df07dd49b706435a094cd58198e4de Mon Sep 17 00:00:00 2001 From: Nicholas Kern Date: Mon, 10 Sep 2018 22:57:22 +0800 Subject: [PATCH 8/9] minor bug / test fixes... --- hera_pspec/plot.py | 12 +- hera_pspec/tests/test_pipes.py | 187 ----- hera_pspec/tests/test_plot.py | 8 +- hera_pspec/tests/test_pspecdata.py | 8 +- hera_pspec/tests/test_uvpspec.py | 6 +- hera_pspec/utils.py | 4 +- .../idr2_preprocessing/preprocess_data.py | 737 ------------------ .../idr2_preprocessing/preprocess_params.yaml | 122 --- pipelines/pspec_pipeline/pspec_batch.sh | 26 - pipelines/pspec_pipeline/pspec_pipe.py | 247 ------ pipelines/pspec_pipeline/pspec_pipe.yaml | 106 --- setup.py | 2 - 12 files changed, 21 insertions(+), 1444 deletions(-) delete mode 100644 hera_pspec/tests/test_pipes.py delete mode 100755 pipelines/idr2_preprocessing/preprocess_data.py delete mode 100644 pipelines/idr2_preprocessing/preprocess_params.yaml delete mode 100644 pipelines/pspec_pipeline/pspec_batch.sh delete mode 100755 pipelines/pspec_pipeline/pspec_pipe.py delete mode 100644 pipelines/pspec_pipeline/pspec_pipe.yaml diff --git a/hera_pspec/plot.py b/hera_pspec/plot.py index 249c1c9f..69e0c865 100644 --- a/hera_pspec/plot.py +++ b/hera_pspec/plot.py @@ -248,7 +248,7 @@ def delay_spectrum(uvp, blpairs, spw, pol, average_blpairs=False, else: c = cax.get_color() if logscale: - _y = np.abs(yp) + _y = np.abs(y) else: _y = y cax = ax.errorbar(x, _y, fmt='none', ecolor=c, yerr=cast(errs[i]), **kwargs) @@ -812,10 +812,14 @@ def delay_wedge(uvp, spw, pol, blpairs=None, times=None, fold=False, delay=True, if cbax is None: cbax = ax cbar = fig.colorbar(cax, ax=cbax) + if deltasq: + p = "\Delta^2" + else: + p = "P" if delay: - p = "P({},\ {})".format(r'\tau', r'|\vec{b}|') + p = "{}({},\ {})".format(p, r'\tau', r'|\vec{b}|') else: - p = "P({},\ {})".format(r'k_\parallel', r'k_\perp') + p = "{}({},\ {})".format(p, r'k_\parallel', r'k_\perp') if log10: psunits = r"$\log_{{10}}\ {}\ [{}]$".format(p, psunits) else: @@ -854,7 +858,7 @@ def delay_wedge(uvp, spw, pol, blpairs=None, times=None, fold=False, delay=True, if not delay: # Get average redshift of spw avg_z = uvp.cosmo.f2z(np.mean(uvp.freq_array[uvp.spw_to_freq_indices(spw)])) - horizons *= uvp.cosmo.tau_to_kpara(avg_z, little_h=little_h) + horizons *= uvp.cosmo.tau_to_kpara(avg_z, little_h=little_h) / 1e9 # iterate over bins and plot lines if rotate: diff --git a/hera_pspec/tests/test_pipes.py b/hera_pspec/tests/test_pipes.py deleted file mode 100644 index 635ba8c0..00000000 --- a/hera_pspec/tests/test_pipes.py +++ /dev/null @@ -1,187 +0,0 @@ -""" -Testing file for scripts in pipelines/ directory -""" -import unittest -import nose.tools as nt -import numpy as np -import os, copy, sys -from hera_pspec.data import DATA_PATH -from pyuvdata import UVData -import subprocess -import hera_pspec as hp -import yaml -import shutil -import fnmatch -import glob - -## Note! In order to debug this script, you should open -## an ipython session, and copy-paste `test_preprocess_run` -## and `test_pspec_run` into the session, and check -## the log.out and err.out files generated... - -class Test_PreProcess(unittest.TestCase): - """ Testing class for preprocess_data.py pipeline """ - - def setUp(self): - pass - - def tearDown(self): - files = ["preproc.yaml", "log.out", "err.out"] - for f in files: - if os.path.exists(f): - try: - os.remove(f) - except: - shutil.rmtree(f) - - def test_preprocess_run(self): - # load default config - cf = hp.utils.load_config(os.path.join(DATA_PATH, "../../pipelines/idr2_preprocessing/preprocess_params.yaml")) - - # edit parameters - cf['analysis']['bl_reformat'] = False - cf['analysis']['rfi_flag'] = False - cf['analysis']['timeavg_sub'] = True - cf['analysis']['time_avg'] = True - cf['analysis']['form_pstokes'] = True - cf['analysis']['fg_filt'] = True - cf['analysis']['multiproc'] = False - cf['analysis']['maxiter'] = 1 - cf['io']['joinlog'] = True - cf['io']['verbose'] = False - cf['io']['work_dir'] = "./" - cf['io']['out_dir'] = "./" - cf['io']['logfile'] = "log.out" - cf['io']['errfile'] = 'err.out' - cf['io']['overwrite'] = True - cf['io']['plot'] = False - cf['data']['data_template'] = os.path.join(DATA_PATH, cf['data']['data_template']) - - # write new config - config = "./preproc.yaml" - with open(config, "w") as f: - f.write(yaml.dump(cf)) - - # clean space - exts = ['T', 'TF', 'TFP', 'TFD'] - pols = ['xx', 'yy'] - stokes = ['pI', 'pQ'] - files = [] - for e in exts: - files += glob.glob("zen.all.??.*uvA{}".format(e)) - files += glob.glob("zen.all.??.tavg.uvA") - for f in files: - if os.path.exists(f): - try: - os.remove(f) - except: - shutil.rmtree(f) - - # run preprocess pipe - out = subprocess.call(["preprocess_data.py", "preproc.yaml"]) - nt.assert_equal(out, 0) - - # check time average (xtalk) subtraction step - for p in pols: - # check files exist - nt.assert_true(os.path.exists("zen.all.{}.tavg.uvA".format(p))) - nt.assert_true(os.path.exists("zen.all.{}.LST.1.06964.uvAT".format(p))) - uvd = UVData() - uvd.read_miriad("zen.all.{}.tavg.uvA".format(p)) - nt.assert_equal(uvd.Ntimes, 1) - nt.assert_equal(uvd.Nbls, 3) - - # check time averaging (FRF) & pstokes steps - for p in pols + stokes: - # check files exist: glob is used b/c starting LST changed slightly - nt.assert_true(len(glob.glob("zen.all.{}.LST.1.0696?.uvATF".format(p)))==1) - - # check FG filtering step - for p in pols + stokes: - # check files exist: glob is used b/c starting LST changed slightly - nt.assert_true(len(glob.glob("zen.all.{}.LST.1.0696?.uvATFD".format(p)))==1) - nt.assert_true(len(glob.glob("zen.all.{}.LST.1.0696?.uvATFP".format(p)))==1) - - # clean space - files = [] - for e in exts: - files += glob.glob("zen.all.??.*uvA{}".format(e)) - files += glob.glob("zen.all.??.tavg.uvA") - for f in files: - if os.path.exists(f): - try: - os.remove(f) - except: - shutil.rmtree(f) - os.remove("preproc.yaml") - os.remove("log.out") - - - -class Test_PspecPipe(unittest.TestCase): - """ Testing class for pspec_pipe.py pipeline """ - - def setUp(self): - pass - - def tearDown(self): - files = ["pspec.h5", "pspec.yaml", "log.out", "err.out"] - for f in files: - if os.path.exists(f): - try: - os.remove(f) - except: - shutil.rmtree(f) - - def test_pspec_run(self): - # load default config - cf = hp.utils.load_config(os.path.join(DATA_PATH, "../../pipelines/pspec_pipeline/pspec_pipe.yaml")) - - # edit parameters - cf['analysis']['run_diff'] = False - cf['analysis']['run_pspec'] = True - cf['analysis']['run_bootstrap'] = True - cf['analysis']['multiproc'] = False - cf['data']['data_template'] = os.path.join(DATA_PATH, cf['data']['data_template']) - cf['io']['joinlog'] = True - cf['io']['verbose'] = False - cf['io']['work_dir'] = "./" - cf['io']['out_dir'] = "./" - cf['io']['logfile'] = "log.out" - cf['io']['errfile'] = 'err.out' - cf['io']['overwrite'] = True - cf['algorithm']['pspec']['beam'] = os.path.join(DATA_PATH, cf['algorithm']['pspec']['beam']) - - # write new config - config = "./pspec.yaml" - with open(config, "w") as f: - f.write(yaml.dump(cf)) - - # clean space - if os.path.exists(cf['algorithm']['pspec']['outfname']): - os.remove(cf['algorithm']['pspec']['outfname']) - - # run preprocess pipe - out = subprocess.call(["pspec_pipe.py", "pspec.yaml"]) - nt.assert_equal(out, 0) - - # check pspec and bootstrap steps - nt.assert_true(os.path.exists(cf['algorithm']['pspec']['outfname'])) - psc = hp.PSpecContainer(cf['algorithm']['pspec']['outfname'], 'r') - groups = psc.groups() - nt.assert_true(len(groups) > 0) - for grp in groups: - spectra = psc.spectra(grp) - nt.assert_true(len(spectra), 2) - nt.assert_true("_avg" not in spectra[0]) - nt.assert_true("_avg" in spectra[1]) - - # clean space - if os.path.exists(cf['algorithm']['pspec']['outfname']): - os.remove(cf['algorithm']['pspec']['outfname']) - os.remove("pspec.yaml") - os.remove("log.out") - - -if __name__ == "__main__": - unittest.main() diff --git a/hera_pspec/tests/test_plot.py b/hera_pspec/tests/test_plot.py index 8ac9bcf5..abcad905 100644 --- a/hera_pspec/tests/test_plot.py +++ b/hera_pspec/tests/test_plot.py @@ -346,10 +346,10 @@ def test_delay_wedge(self): edgecolor='none', flip_xax=False, flip_yax=False, lw=2) plt.close() - - - - + # test exceptions + nt.assert_raises(ValueError, plot.delay_wedge, uvp, 0, 'xx', component='foo') + plt.close() + if __name__ == "__main__": unittest.main() diff --git a/hera_pspec/tests/test_pspecdata.py b/hera_pspec/tests/test_pspecdata.py index 97333125..6c968969 100644 --- a/hera_pspec/tests/test_pspecdata.py +++ b/hera_pspec/tests/test_pspecdata.py @@ -857,7 +857,7 @@ def test_rephase_to_dset(self): # rephase and get pspec ds.rephase_to_dset(0) uvp2 = ds.pspec(bls, bls, (0, 1), pols=('xx','xx'), verbose=False) - blp = (0, ((37,39),(37,39)), 'XX') + blp = (0, ((37,39),(37,39)), 'xx') nt.assert_true(np.isclose(np.abs(uvp2.get_data(blp)/uvp1.get_data(blp)), 1.0).min()) def test_Jy_to_mK(self): @@ -1063,7 +1063,7 @@ def test_pspec(self): spw_ranges=[(20, 30)]) nt.assert_equal(len(ds._identity_Y), len(ds._identity_G), len(ds._identity_H)) nt.assert_equal(len(ds._identity_Y), 1) - nt.assert_equal(ds._identity_Y.keys()[0], ((0, 24, 25, 'XX'), (1, 24, 25, 'XX'))) + nt.assert_equal(ds._identity_Y.keys()[0], ((0, 24, 25, 'xx'), (1, 24, 25, 'xx'))) # assert caching is not used when inappropriate ds.dsets[0].flag_array[ds.dsets[0].antpair2ind(37, 38, ordered=False), :, 25, :] = True uvp = ds.pspec([(24, 25), (37, 38)], [(24, 25), (37, 38)], (0, 1), ('xx', 'xx'), @@ -1071,8 +1071,8 @@ def test_pspec(self): spw_ranges=[(20, 30)]) nt.assert_equal(len(ds._identity_Y), len(ds._identity_G), len(ds._identity_H)) nt.assert_equal(len(ds._identity_Y), 2) - nt.assert_true(((0, 24, 25, 'XX'), (1, 24, 25, 'XX')) in ds._identity_Y.keys()) - nt.assert_true(((0, 37, 38, 'XX'), (1, 37, 38, 'XX')) in ds._identity_Y.keys()) + nt.assert_true(((0, 24, 25, 'xx'), (1, 24, 25, 'xx')) in ds._identity_Y.keys()) + nt.assert_true(((0, 37, 38, 'xx'), (1, 37, 38, 'xx')) in ds._identity_Y.keys()) def test_normalization(self): # Test Normalization of pspec() compared to PAPER legacy techniques diff --git a/hera_pspec/tests/test_uvpspec.py b/hera_pspec/tests/test_uvpspec.py index f70dfd5d..2737a7e4 100644 --- a/hera_pspec/tests/test_uvpspec.py +++ b/hera_pspec/tests/test_uvpspec.py @@ -82,8 +82,8 @@ def test_get_funcs(self): nt.assert_equal(blps, [((1, 2), (1, 2)), ((1, 3), (1, 3)), ((2, 3), (2, 3))]) # test get all keys keys = self.uvp.get_all_keys() - nt.assert_equal(keys, [(0, ((1, 2), (1, 2)), 'XX'), (0, ((1, 3), (1, 3)), 'XX'), - (0, ((2, 3), (2, 3)), 'XX')]) + nt.assert_equal(keys, [(0, ((1, 2), (1, 2)), 'xx'), (0, ((1, 3), (1, 3)), 'xx'), + (0, ((2, 3), (2, 3)), 'xx')]) # test omit_flags self.uvp.integration_array[0][self.uvp.blpair_to_indices(((1, 2), (1, 2)))[:2]] = 0.0 nt.assert_equal(self.uvp.get_integrations((0, ((1, 2), (1, 2)), 'xx'), omit_flags=True).shape, (8,)) @@ -448,7 +448,7 @@ def test_combine_uvpspec(self): out = uvpspec.combine_uvpspec([uvp1, uvp2], verbose=False) nt.assert_equal(out.Npols, 2) nt.assert_true(len(set(out.pol_array) ^ set([-5, -6])) == 0) - key = (0, ((37, 38), (38, 39)), 'XX') + key = (0, ((37, 38), (38, 39)), 'xx') nt.assert_true(np.all(np.isclose(out.get_nsamples(key), np.ones(10, dtype=np.float64)))) nt.assert_true(np.all(np.isclose(out.get_integrations(key), 190 * np.ones(10, dtype=np.float64), atol=5, rtol=2))) diff --git a/hera_pspec/utils.py b/hera_pspec/utils.py index 37388495..3626fa6a 100644 --- a/hera_pspec/utils.py +++ b/hera_pspec/utils.py @@ -751,10 +751,10 @@ def get_blvec_reds(blvecs, bl_error_tol=1.0, match_bl_lens=False): # append to list if unique within tolerance if match_bl_lens: # match only on bl length - match = [np.all(np.isclose(bll, bl_len, atol=bl_error_tol)) for bll in red_bl_len] + match = [np.all(np.isclose(bll, bl_len, rtol=0.0, atol=bl_error_tol)) for bll in red_bl_len] else: # match on full bl vector - match = [np.all(np.isclose(blv, bl_vec, atol=bl_error_tol)) for blv in red_bl_vec] + match = [np.all(np.isclose(blv, bl_vec, rtol=0.0, atol=bl_error_tol)) for blv in red_bl_vec] if np.any(match): match_id = np.where(match)[0][0] red_bl_grp[match_id].append(bl) diff --git a/pipelines/idr2_preprocessing/preprocess_data.py b/pipelines/idr2_preprocessing/preprocess_data.py deleted file mode 100755 index bf8b038f..00000000 --- a/pipelines/idr2_preprocessing/preprocess_data.py +++ /dev/null @@ -1,737 +0,0 @@ -#!/usr/bin/env python2 -""" -preprocess_data.py ------------------------------------------ -Copyright (c) 2018 The HERA Collaboration - -This script is used in the IDR2.1 power -spectrum pipeline as a pre-processing step after -calibration, RFI-flagging and LSTbinning. This -additional processing includes RFI-flagging, timeavg subtraction, -fringe-rate filtering, pseudo-stokes visibility formation, -and foreground filtering. - -See preprocess_params.yaml for relevant parameter selections. -""" -import multiprocess -import numpy as np -import hera_cal as hc -import hera_pspec as hp -import hera_qm as hq -import pyuvdata -from pyuvdata import UVData -import pyuvdata.utils as uvutils -import os -import sys -import glob -import yaml -from datetime import datetime -import uvtools as uvt -import json -import itertools -import aipy -import shutil -from collections import OrderedDict as odict - - -#------------------------------------------------------------------------------- -# Parse YAML Configuration File -#------------------------------------------------------------------------------- -# Get config and load dictionary -config = sys.argv[1] -cf = hp.utils.load_config(config) - -# Consolidate IO, data and analysis parameter dictionaries -params = odict(cf['io'].items() + cf['data'].items() + cf['analysis'].items()) -assert len(params) == len(cf['io']) + len(cf['data']) + len(cf['analysis']), ""\ - "Repeated parameters found within the scope of io, data and analysis dicts" -algs = cf['algorithm'] - -# Extract certain parameters used across the script -verbose = params['verbose'] -overwrite = params['overwrite'] -pols = params['pols'] -data_template = os.path.join(params['data_root'], params['data_template']) -data_suffix = os.path.splitext(data_template)[1][1:] - -# Open logfile -logfile = os.path.join(params['out_dir'], params['logfile']) -if os.path.exists(logfile) and params['overwrite'] == False: - raise IOError("logfile {} exists and overwrite == False, quitting pipeline...".format(logfile)) -lf = open(logfile, "w") -if params['joinlog']: - ef = lf -else: - ef = open(os.path.join(params['out_dir'], params['errfile']), "w") - -time = datetime.utcnow() -hp.utils.log("Starting preprocess pipeline on {}\n{}\n".format(time, '-'*60), - f=lf, verbose=verbose) -hp.utils.log(json.dumps(cf, indent=1) + '\n', f=lf, verbose=verbose) - -# Change to working dir -os.chdir(params['work_dir']) - -# out_dir should be cleared before each run: issue a warning if not the case -outdir = os.path.join(params['work_dir'], params['out_dir']) -oldfiles = glob.glob(outdir+'/*') -if len(oldfiles) > 0: - hp.utils.log("\n{}\nWARNING: out_dir should be cleaned before each new run to " \ - "ensure proper functionality.\nIt seems like some files currently " \ - "exist in {}\n{}\n".format('-'*50, outdir, '-'*50), f=lf, verbose=verbose) - -# Define history prepend function -def prepend_history(action, param_dict): - """ create a history string to prepend to data files """ - dict_str = '\n'.join(["{} : {}".format(*_d) for _d in param_dict.items()]) - time = datetime.utcnow() - hist = "\nRan preprocess_data.py {} step at\nUTC {} with \nhera_pspec [{}], "\ - "hera_cal [{}],\nhera_qm [{}] and pyuvdata [{}]\nwith {} algorithm "\ - "attrs:\n{}\n{}\n".format(action, time, - hp.version.git_hash[:10], - hc.version.git_hash[:10], - hq.version.git_hash[:10], - pyuvdata.version.git_hash[:10], - action, '-'*50, - dict_str) - return hist - -# Assign iterator function -if params['multiproc']: - pool = multiprocess.Pool(params['nproc']) - M = pool.map -else: - M = map - -#------------------------------------------------------------------------------- -# Reformat Data by Baseline Type -#------------------------------------------------------------------------------- -if params['reformat']: - - # Start block - time = datetime.utcnow() - hp.utils.log("\n{}\nstarting baseline reformatting: {}\n".format("-"*60, time), - f=lf, verbose=verbose) - - # Get datafiles - datafiles, datapols = uvt.utils.search_data( - data_template.format(group=params['groupname']), - pols, matched_pols=False, - reverse_nesting=False, - flatten=False ) - - # Get redundant groups as a good split for parallelization - reds, lens, angs = hp.utils.get_reds(datafiles[0][0], bl_error_tol=1.0, - add_autos=True, - bl_len_range=params['bl_len_range'], - bl_deg_range=params['bl_deg_range']) - - # Iterate over polarization group - for i, dfs in enumerate(datafiles): - - # Setup bl reformat function - def bl_reformat(j, i=i, datapols=datapols, dfs=dfs, lens=lens, - angs=angs, reds=reds, data_suffix=data_suffix, - p=algs['reformat'], params=params): - try: - if not p['bl_len_range'][0] < lens[j] < p['bl_len_range'][1]: - return 0 - outname = p['reformat_outfile'].format(len=int(round(lens[j])), - deg=int(round(angs[j])), - pol=datapols[i][0], - suffix=data_suffix) - outname = os.path.join(params['out_dir'], outname) - if os.path.exists(outname) and overwrite == False: - return 1 - uvd = UVData() - uvd.read_miriad(dfs, bls=reds[j]) - uvd.write_miriad(outname, clobber=True) - uvd.history = "{}{}".format(prepend_history("BL REFORMAT", p), - uvd.history) - if params['plot']: - hp.plot.plot_uvdata_waterfalls(uvd, - outname + ".{pol}.{bl}", - data='data', - plot_mode='log', - format='png') - except: - hp.utils.log("\njob {} threw exception:".format(j), - f=ef, tb=sys.exc_info(), verbose=verbose) - return 1 - return 0 - - # Launch jobs - failures = hp.utils.job_monitor(bl_reformat, range(len(reds)), \ - "BL REFORMAT: pol {}".format(pol), - M=M, lf=lf, maxiter=params['maxiter'], - verbose=verbose) - - # Edit data template - data_template = os.path.join( - params['out_dir'], - algs['reformat']['new_data_template'].format( - pol='{pol}', - suffix=data_suffix)) - time = datetime.utcnow() - hp.utils.log("\nfinished baseline reformatting: {}\n{}".format(time, "-"*60), - f=lf, verbose=verbose) - -#------------------------------------------------------------------------------- -# RFI-Flag -#------------------------------------------------------------------------------- -if params['rfi_flag']: - - # Start block - time = datetime.utcnow() - hp.utils.log("\n{}\nstarting RFI flagging: {}\n".format("-"*60, time), - f=lf, verbose=verbose) - - # Get datafiles - datafiles, datapols = uvt.utils.search_data( - data_template.format(group=params['groupname']), - pols, matched_pols=False, - reverse_nesting=False, - flatten=True ) - - # Setup RFI function - def run_xrfi(i, datafiles=datafiles, p=cf['algorithm']['xrfi'], params=params): - try: - # Setup delay filter class as container - df = datafiles[i] - F = hc.delay_filter.Delay_Filter() - F.load_data(df) # load data - - # RFI flag if desired - if run_xrfi: - for k in F.data.keys(): - new_f = hq.xrfi.xrfi(F.data[k], - f=F.flags[k], - **p['xrfi_params']) - F.flags[k] += new_f - - # Write to file - add_to_history = prepend_history("XRFI", p) - outname = os.path.join(params['out_dir'], - os.path.basename(df) + p['file_ext']) - hc.io.update_vis(df, outname, - filetype_in='miriad', - filetype_out='miriad', - data=F.data, - flags=F.flags, - add_to_history=add_to_history, - clobber=overwrite) - - except: - hp.utils.log("\njob {} threw exception:".format(i), - f=ef, tb=sys.exc_info(), verbose=verbose) - return 1 - return 0 - - # Launch jobs - failures = hp.utils.job_monitor(run_xrfi, range(len(datafiles)), - "XRFI", M=M, lf=lf, - maxiter=params['maxiter'], - verbose=verbose) - - # Update template - data_template = os.path.join(params['out_dir'], - os.path.basename(data_template) \ - + algs['rfi_flag']['file_ext']) - data_suffix += algs['rfi_flag']['file_ext'] - - time = datetime.utcnow() - hp.utils.log("\nfinished RFI flagging: {}\n{}".format(time, "-"*60), - f=lf, verbose=verbose) - -#------------------------------------------------------------------------------- -# Time Average Subtraction -#------------------------------------------------------------------------------- -if params['timeavg_sub']: - - # Start block - time = datetime.utcnow() - hp.utils.log("\n{}\nstarting full time-average spectra and subtraction: {}\n".format("-"*60, time), - f=lf, verbose=verbose) - - # get datafiles - datafiles, datapols = uvt.utils.search_data( - data_template.format(group=params['groupname'], pol='{pol}'), - pols, matched_pols=False, - reverse_nesting=False, - flatten=False) - - # get redundant groups as a good split for parallelization - reds, lens, angs = hp.utils.get_reds(datafiles[0][0], - bl_error_tol=1.0, - add_autos=True, - bl_len_range=params['bl_len_range'], - bl_deg_range=params['bl_deg_range']) - - # iterate over pols - for i, dfs in enumerate(datafiles): - pol = datapols[i][0] - - # write full tavg function - def full_tavg(j, pol=pol, lens=lens, angs=angs, reds=reds, dfs=dfs, - data_suffix=data_suffix, p=algs['timeavg_sub'], - params=params): - try: - # load data into uvdata - uvd = UVData() - # read data, catch ValueError - try: - uvd.read_miriad(dfs, bls=reds[j], polarizations=[pol]) - except ValueError: - hp.utils.log("job {} failed b/c no data is present given bls and/or pol selection".format(j)) - return 0 - - # Perform full time-average to get spectrum - F = hc.frf.FRFilter() # instantiate FRF object - F.load_data(uvd) # load data - F.timeavg_data(1e10, rephase=False, verbose=verbose) - - # Delay Filter if desired - if p['dly_filter']: - - # RFI Flag - for k in F.avg_data.keys(): - # RFI Flag - new_f = hq.xrfi.xrfi(F.avg_data[k], - f=F.avg_flags[k], - **p['rfi_params']) - F.avg_flags[k] += new_f - - # Delay Filter - DF = hc.delay_filter.Delay_Filter() - DF.load_data(F.input_data) - DF.data = F.avg_data - DF.flags = F.avg_flags - DF.run_filter(**p['dly_params']) - - # Replace timeavg spectrum with CLEAN model - F.avg_data = DF.CLEAN_models - - # Unflag all frequencies - for k in F.avg_flags.keys(): - # only unflag spectra from non-xants - if np.min(F.avg_flags[k]) == False: - F.avg_flags[k][:] = False - - # Write timeavg spectrum - _len = lens[j] - _deg = angs[j] - _tavg = "zen.{group}.{pol}.{len:03d}_{deg:03d}.{tavg_tag}.{suffix}" - tavg_file = _tavg.format( group=params['groupname'], - pol=pol, len=int(_len), deg=int(_deg), - tavg_tag=p['tavg_tag'], - suffix=data_suffix ) - tavg_file = os.path.join(params['out_dir'], tavg_file) - add_to_history = prepend_history("FULL TIME AVG", p) - F.write_data(tavg_file, write_avg=True, overwrite=overwrite, - add_to_history=add_to_history) - except: - hp.utils.log("\njob {} threw exception:".format(j), - f=ef, tb=sys.exc_info(), verbose=verbose) - return 1 - return 0 - - # Launch jobs - failures = hp.utils.job_monitor(full_tavg, range(len(reds)), - "FULL TAVG: pol {}".format(pol), - M=M, lf=lf, maxiter=params['maxiter'], - verbose=verbose) - - # Collate tavg spectra into a single file - tavgfiles = sorted(glob.glob( - os.path.join( - params['out_dir'], - "zen.{group}.{pol}.*.{tavg_tag}.{suffix}".format( - group=params['groupname'], - pol=pol, - tavg_tag=algs['timeavg_sub']['tavg_tag'], - suffix=data_suffix)) )) - uvd = UVData() - uvd.read_miriad(tavgfiles[::-1]) - tavg_out = os.path.join( params['out_dir'], - "zen.{group}.{pol}.{tavg_tag}.{suffix}".format( - group=params['groupname'], - pol=pol, - tavg_tag=algs['timeavg_sub']['tavg_tag'], - suffix=data_suffix)) - - # Write tavg data file - uvd.write_miriad(tavg_out, clobber=overwrite) - - # Plot waterfalls if requested - if params['plot']: - hp.plot.plot_uvdata_waterfalls(uvd, - tavg_out + ".{pol}.{bl}", - data='data', - plot_mode='log', - format='png') - for tf in tavgfiles: - if os.path.exists(tf): - shutil.rmtree(tf) - - # Write tavg subtraction function - def tavg_sub(j, dfs=dfs, pol=pol, tavg_file=tavg_out, - p=cf['algorithm']['timeavg_sub'], params=params): - try: - # Load data file - uvd = UVData() - df = dfs[j] - uvd.read_miriad(df) - - # Get full tavg spectra - tavg = UVData() - tavg.read_miriad(tavg_file) - - # Subtract timeavg spectrum from data - for bl in np.unique(uvd.baseline_array): - bl_inds = np.where(uvd.baseline_array == bl)[0] - polnum = uvutils.polstr2num(pol) - pol_ind = np.where(uvd.polarization_array == polnum)[0] - if bl in tavg.baseline_array: - uvd.data_array[bl_inds, :, :, pol_ind] -= tavg.get_data(bl)[None] - else: - uvd.flag_array[bl_inds, :, :, pol_ind] = True - if verbose: - print "baseline {} not found in time-averaged spectrum".format(bl) - - # Put uniq_bls in if it doesn't exist - if not uvd.extra_keywords.has_key('uniq_bls'): - uvd.extra_keywords['uniq_bls'] = json.dumps( - np.unique(uvd.baseline_array).tolist() ) - - # Write tavg-subtracted data - out_df = os.path.join(params['out_dir'], - os.path.basename(df) + p['file_ext']) - uvd.history = "{}{}".format(prepend_history("TAVG SUB", p), - uvd.history) - uvd.write_miriad(out_df, clobber=overwrite) - - - # Plot waterfalls if requested - if params['plot']: - hp.plot.plot_uvdata_waterfalls(uvd, - out_df + ".{pol}.{bl}", - data='data', - plot_mode='log', - format='png') - - except: - hp.utils.log("\njob {} threw exception:".format(i), - f=ef, tb=sys.exc_info(), verbose=verbose) - return 1 - return 0 - - # Launch jobs - failures = hp.utils.job_monitor( tavg_sub, - range(len(dfs)), - "TAVG SUB: pol {}".format(pol), - M=M, lf=lf, maxiter=params['maxiter'], - verbose=verbose ) - - time = datetime.utcnow() - hp.utils.log("\nfinished full time-average spectra and subtraction: {}\n{}".format(time, "-"*60), - f=lf, verbose=verbose) - - data_template = os.path.join( params['out_dir'], - os.path.basename(data_template) \ - + algs['timeavg_sub']['file_ext'] ) - data_suffix += algs['timeavg_sub']['file_ext'] - -#------------------------------------------------------------------------------- -# Time Averaging (i.e. Fringe Rate Filtering) -#------------------------------------------------------------------------------- -if params['time_avg']: - # Start block - time = datetime.utcnow() - hp.utils.log("\n{}\nstarting time averaging: {}\n".format("-"*60, time), - f=lf, verbose=verbose) - - # Get datafiles - datafiles, datapols = uvt.utils.search_data( - data_template.format(group=params['groupname'], - pol='{pol}'), - pols, matched_pols=False, - reverse_nesting=False, - flatten=False) - - # Get redundant groups as a good split for parallelization - reds, lens, angs = hp.utils.get_reds( datafiles[0][0], - bl_error_tol=1.0, - add_autos=True, - bl_len_range=params['bl_len_range'], - bl_deg_range=params['bl_deg_range'] ) - - # Iterate over pol groups - for i, dfs in enumerate(datafiles): - pol = datapols[i][0] - - def time_average(j, dfs=dfs, pol=pol, lens=lens, angs=angs, reds=reds, - data_suffix=data_suffix, p=cf['algorithm']['tavg'], - params=params): - try: - # Load data into uvdata - uvd = UVData() - - # Read data, catch ValueError - try: - uvd.read_miriad(dfs, bls=reds[j], polarizations=[pol]) - except ValueError: - hp.utils.log("job {} failed w/ ValueError, probably b/c no data is present given bls and/or pol selection".format(j)) - return 0 - - # Perform time-average - F = hc.frf.FRFilter() # instantiate FRF object - F.load_data(uvd) # load data - F.timeavg_data( p['t_window'], - rephase=p['tavg_rephase'], - verbose=verbose ) - - # Write timeavg spectrum - _len = lens[j] - _deg = angs[j] - tavg_file = "zen.{group}.{pol}.{len:03d}_{deg:03d}.{suffix}".format( - group=params['groupname'], - pol=pol, - len=int(_len), - deg=int(_deg), - suffix=data_suffix) - tavg_file = os.path.join(params['out_dir'], - tavg_file + p['file_ext']) - add_to_history = prepend_history("TIME AVERAGE", p) - F.write_data(tavg_file, write_avg=True, overwrite=overwrite, - add_to_history=add_to_history) - except: - hp.utils.log("\njob {} threw exception:".format(j), - f=ef, tb=sys.exc_info(), verbose=verbose) - return 1 - return 0 - - # Launch jobs - failures = hp.utils.job_monitor( time_average, - range(len(reds)), - "TIME AVERAGE: pol {}".format(pol), - M=M, lf=lf, - maxiter=params['maxiter'], - verbose=verbose ) - - # Collate averaged data into time chunks - tavg_files = os.path.join( params['out_dir'], - "zen.{group}.{pol}.*.{suffix}".format( - group=params['groupname'], - pol=pol, - suffix=data_suffix \ - + algs['tavg']['file_ext'])) - tavg_files = sorted(glob.glob(tavg_files)) - assert len(tavg_files) > 0, "len(tavg_files) == 0" - - # Pick one file to get full time information from - uvd = UVData() - uvd.read_miriad(tavg_files[-1]) - times = np.unique(uvd.time_array) - Ntimes = len(times) - - # break into subfiles - Nfiles = int(np.ceil(Ntimes / float(algs['tavg']['file_Ntimes']))) - times = [ times[ i*algs['tavg']['file_Ntimes'] : - (i+1)*algs['tavg']['file_Ntimes'] ] - for i in range(Nfiles) ] - - def reformat_files(j, pol=pol, tavg_files=tavg_files, times=times, - data_suffix=data_suffix, p=cf['algorithm']['tavg'], - params=params): - try: - uvd = UVData() - uvd.read_miriad(tavg_files[::-1], - time_range=[times[j].min()-1e-8, - times[j].max()+1e-8], - polarizations=[pol]) - lst = uvd.lst_array[0] - np.median(uvd.integration_time) / 2. \ - * 2 * np.pi / (3600. * 24) - outfile = os.path.join( - params['out_dir'], - "zen.{group}.{pol}.LST.{LST:.5f}.{suffix}".format( - group=params['groupname'], - pol=pol, LST=lst, - suffix=data_suffix + p['file_ext'])) - uvd.write_miriad(outfile, clobber=overwrite) - - # Plot waterfalls if requested - if params['plot']: - hp.plot.plot_uvdata_waterfalls(uvd, - outfile + ".{pol}.{bl}", - data='data', - plot_mode='log', - format='png') - - except: - hp.utils.log("\njob {} threw exception:".format(j), - f=ef, tb=sys.exc_info(), verbose=verbose) - return 1 - - return 0 - - # Launch jobs - failures = hp.utils.job_monitor(reformat_files, - range(len(times)), - "TAVG REFORMAT: pol {}".format(pol), - M=M, lf=lf, maxiter=params['maxiter'], - verbose=verbose) - - # Clean up time averaged files - for f in tavg_files: - if os.path.exists(f): - shutil.rmtree(f) - - data_template = os.path.join(params['out_dir'], - os.path.basename(data_template) \ - + algs['tavg']['file_ext']) - data_suffix += algs['tavg']['file_ext'] - - time = datetime.utcnow() - hp.utils.log("\nfinished time averaging: {}\n{}".format(time, "-"*60), - f=lf, verbose=verbose) - -#------------------------------------------------------------------------------- -# Form Pseudo-Stokes Visibilities -#------------------------------------------------------------------------------- -if params['form_pstokes']: - # Start block - time = datetime.utcnow() - hp.utils.log("\n{}\nstarting pseudo-stokes: {}\n".format("-"*60, time), f=lf, verbose=verbose) - - # Get datafiles with reversed nesting - datafiles, datapols = uvt.utils.search_data( - data_template.format(group=params['groupname'], - pol='{pol}'), - pols, reverse_nesting=True) - - # Write pseudo-Stokes function - def make_pstokes(i, datafiles=datafiles, datapols=datapols, - p=cf['algorithm']['pstokes'], params=params): - try: - # Get all pol files for unique datafile - dfs = datafiles[i] - dps = datapols[i] - - # Load data - dsets = [] - for df in dfs: - uvd = UVData() - uvd.read_miriad(df) - dsets.append(uvd) - - # Iterate over outstokes - for pstokes in p['outstokes']: - try: - ds = hp.pstokes.filter_dset_on_stokes_pol(dsets, pstokes) - ps = hp.pstokes.construct_pstokes(ds[0], ds[1], pstokes=pstokes) - outfile = os.path.basename(dfs[0]).replace( - ".{}.".format(dps[0]), - ".{}.".format(pstokes)) - outfile = os.path.join(params['out_dir'], outfile) - ps.history = "{}{}".format(prepend_history("FORM PSTOKES", p), - ps.history) - ps.write_miriad(outfile, clobber=overwrite) - - # Plot waterfalls if requested - if params['plot']: - hp.plot.plot_uvdata_waterfalls(ps, - outfile + ".{pol}.{bl}", - data='data', - plot_mode='log', - format='png') - - except AssertionError: - hp.utils.log("failed to make pstokes {} for job {}:".format(pstokes, i), - f=ef, tb=sys.exc_info(), verbose=verbose) - - except: - hp.utils.log("job {} threw exception:".format(i), - f=ef, tb=sys.exc_info(), verbose=verbose) - return 1 - - return 0 - - # Launch jobs - failures = hp.utils.job_monitor( make_pstokes, - range(len(datafiles)), - "PSTOKES", - M=M, lf=lf, maxiter=params['maxiter'], - verbose=verbose) - - # Add pstokes pols to pol list for downstream calculations - pols += algs['pstokes']['outstokes'] - - time = datetime.utcnow() - hp.utils.log("\nFinished pseudo-stokes: {}\n{}".format(time, "-"*60), - f=lf, verbose=verbose) - -#------------------------------------------------------------------------------- -# Foreground Filtering (and data in-painting) -#------------------------------------------------------------------------------- -if params['fg_filt']: - - # Start block - time = datetime.utcnow() - hp.utils.log("\n{}\nstarting foreground filtering: {}\n".format("-"*60, time), - f=lf, verbose=verbose) - - # Get flattened datafiles - datafiles, datapols = uvt.utils.search_data( - data_template.format(group=params['groupname'], - pol='{pol}'), - pols, matched_pols=False, - reverse_nesting=False, - flatten=True) - - # Write fgfilt function - def fg_filter(i, datafiles=datafiles, p=cf['algorithm']['fg_filt'], - params=params): - try: - # Get datafile and start FGFilter - df = datafiles[i] - DF = hc.delay_filter.Delay_Filter() - DF.load_data(df, filetype='miriad') - - # Run filter - DF.run_filter(**p['filt_params']) - - # Write filtered term - outfile = os.path.join(params['out_dir'], - os.path.basename(df) + p['filt_file_ext']) - add_to_history = prepend_history("FG FILTER", p) - DF.write_filtered_data( outfile, - filetype='miriad', - clobber=overwrite, - add_to_history=add_to_history ) - - # Write original data with in-paint - outfile = os.path.join(params['out_dir'], - os.path.basename(df) + p['inpaint_file_ext']) - add_to_history = prepend_history("DATA INPAINT", p) - DF.write_filtered_data( outfile, - filetype='miriad', - clobber=overwrite, - write_filled_data=True, - add_to_history=add_to_history ) - except: - hp.utils.log("job {} threw exception:".format(i), - f=ef, tb=sys.exc_info(), verbose=verbose) - return 1 - return 0 - - # Launch jobs - failures = hp.utils.job_monitor( fg_filter, - range(len(datafiles)), - "FG FILTER", - M=M, lf=lf, maxiter=params['maxiter'], - verbose=verbose ) - - time = datetime.utcnow() - hp.utils.log("\nfinished foreground-filtering: {}\n{}".format(time, "-"*60), - f=lf, verbose=verbose) - diff --git a/pipelines/idr2_preprocessing/preprocess_params.yaml b/pipelines/idr2_preprocessing/preprocess_params.yaml deleted file mode 100644 index 78da2fc2..00000000 --- a/pipelines/idr2_preprocessing/preprocess_params.yaml +++ /dev/null @@ -1,122 +0,0 @@ -# preprocess_params.yaml -# -# hera_pspec IDR2.1 preprocessing pipeline -# configuration file - -#--------------------------------------------------------------- -# IO Parameters -#--------------------------------------------------------------- -io : - work_dir : './' # directory to work in - out_dir : './' # directory to put output: root is work_dir - logfile : 'preprocess_out.log' # logfile: root is work_dir - errfile : 'preprocess_err.log' # error file: root is work_dir - joinlog : False # redirect error output into logfile - overwrite : True # overwrite - plot : True - verbose : True # report feedback to standard output - -#--------------------------------------------------------------- -# Analysis Parameters -#--------------------------------------------------------------- -analysis : - reformat : False # reformat data by baseline type - rfi_flag : False # RFI-Flag data - timeavg_sub : True # full time-avg subtraction - time_avg : True # time average data - form_pstokes : True # form pseudo-Stokes visibilities - fg_filt : True # foreground filter (and/or data-inpaint) - - multiproc : False # use multiprocess module - nproc : 8 # number of processes to spawn - maxiter : 1 - -#--------------------------------------------------------------- -# Data Parameters -#--------------------------------------------------------------- -data : - # list of linear polarizations of data to format data template with - pols : - - 'xx' - - 'yy' - - # a groupname to format data template with - groupname : "all" - - # a single path to a (glob-parseable) template with {pol} and {group} fields - data_root : "./" # data root dir - data_template : 'zen.{group}.{pol}.LST.1.?????.uvA' # data basename - - # data baseline-type length range [meters] and angle range [degrees] - bl_len_range : [14, 16] - bl_deg_range : [0, 1] - -#--------------------------------------------------------------- -# Algorithm Parameters -#--------------------------------------------------------------- -algorithm : - - # Data Reformatting - reformat : - # output file template basename: root is out_dir - reformat_outfile : "zen.{group}.{len:03d}_{deg:03d}.{pol}.{suffix}" - # a new data template to adopt w/ {group}, {pol} and {suffix} fields: root is out_dir - new_data_template : "zen.{group}.*.{pol}.{suffix}" - - # RFI Flagging - rfi_flag : - file_ext : 'R' # file extension - run_xrfi : True # run xrfi if True - xrfi_params: # if run_xrfi, set xrfi parameters - Kt : 10 - Kf : 15 - sig_init : 6 - sig_adj : 2 - - # Time-Avg Subtraction - timeavg_sub : - file_ext : "T" # file extension of timeavg-subtracted data - tavg_tag : "tavg" # filename tag of timeavg-spectrum - dly_filter : True # RFI flag and delay filter timeavg spectrum - rfi_params : # if dly_inpaint, these are the xrfi params - sig_init : 2 - sig_adj : 1 - Kf : 20 - Kt : 1 - dly_params : # if dly_filter, these are the delay filter params - min_dly : 1500.0 - horizon : 0.0 - tol : 0.0001 - maxiter : 100 - window : 'blackman-harris' - flag_nchan_low : 25 - flag_nchan_high : 25 - gain : 0.1 - - # Time Averaging - tavg : - file_ext : "F" # file extension of time-averaged data - t_window : 200.0 # width of averaging window in seconds - file_Ntimes : 10 # output file duration in seconds: Ntimes ~ t_file / t_window - tavg_rephase : True # rephase data to window-center LST before average - - # Forming pStokes - pstokes : - outstokes : ['pI', 'pQ'] - - # Foreground Filtering - fg_filt : - filt_file_ext : "D" - inpaint_file_ext : "P" - filt_params : # kwarg parameters in hera_cal.delay_filer.run_filter - standoff : 50.0 - horizon : 1.0 - min_dly : 0.0 - tol : 0.0001 - maxiter : 100 - window : 'tukey' - skip_wgt : 0.5 - flag_nchan_low : 50 - flag_nchan_high : 50 - gain : 0.1 - diff --git a/pipelines/pspec_pipeline/pspec_batch.sh b/pipelines/pspec_pipeline/pspec_batch.sh deleted file mode 100644 index d8889282..00000000 --- a/pipelines/pspec_pipeline/pspec_batch.sh +++ /dev/null @@ -1,26 +0,0 @@ -#!/bin/bash -#PBS -q hera -#PBS -V -#PBS -j oe -#PBS -o pspec_batch.out -#PBS -l nodes=1:ppn=8 -#PBS -l walltime=24:00:00 -#PBS -l vmem=128gb, mem=128gb - -# start script -echo "starting power spectrum pipeline: $(date)" - -# copy parameter files to more discrete names -# and put a readonly lock on them -cp pspec_pipe.yaml .pspec_pipe.yaml -chmod 444 .pspec_pipe.yaml - -# run scripts -pspec_pipe.py .pspec_pipe.yaml - -# unlock and clean-up parameter files -chmod 775 .pspec_pipe.yaml -rm .pspec_pipe.yaml - -# end script -echo "ending power spectrum pipeline: $(date)" diff --git a/pipelines/pspec_pipeline/pspec_pipe.py b/pipelines/pspec_pipeline/pspec_pipe.py deleted file mode 100755 index 03015076..00000000 --- a/pipelines/pspec_pipeline/pspec_pipe.py +++ /dev/null @@ -1,247 +0,0 @@ -#!/usr/bin/env python2 -""" -pspec_pipe.py ------------------------------------------ -Copyright (c) 2018 The HERA Collaboration - -This script is used as the IDR2.1 power -spectrum pipeline. - -See pspec_pipe.yaml for relevant parameter selections. -""" -import multiprocess -import numpy as np -import hera_pspec as hp -from pyuvdata import UVData -import pyuvdata.utils as uvutils -import os -import sys -import glob -import yaml -from datetime import datetime -import uvtools as uvt -import json -import itertools -import aipy -import shutil -from collections import OrderedDict as odict - - -#------------------------------------------------------------------------------- -# Parse YAML Configuration File -#------------------------------------------------------------------------------- -# get config and load dictionary -config = sys.argv[1] -cf = hp.utils.load_config(config) - -# consolidate IO, data and analysis paramemeter dictionaries -params = odict(cf['io'].items() + cf['data'].items() + cf['analysis'].items()) -assert len(params) == len(cf['io']) + len(cf['data']) + len(cf['analysis']), ""\ - "Repeated parameters found within the scope of io, data and analysis dicts" -algs = cf['algorithm'] -params['data_template'] = os.path.join(params['data_root'], params['data_template']) -if params['std_template'] is not None: - params['std_template'] = os.path.join(params['data_root'], params['std_template']) - -# open logfile -logfile = os.path.join(params['out_dir'], params['logfile']) -if os.path.exists(logfile) and params['overwrite'] == False: - raise IOError("logfile {} exists and overwrite == False, quitting pipeline...".format(logfile)) -lf = open(logfile, "w") -if params['joinlog']: - ef = lf -else: - ef = open(os.path.join(params['out_dir'], params['errfile']), "w") -time = datetime.utcnow() -hp.utils.log("Starting pspec pipeline on {}\n{}\n".format(time, '-'*60), f=lf, verbose=params['verbose']) -hp.utils.log(json.dumps(cf, indent=1) + '\n', f=lf, verbose=params['verbose']) - -# define history prepend function -def prepend_history(action, param_dict): - """ create a history string to prepend to data files """ - dict_str = '\n'.join(["{} : {}".format(*_d) for _d in param_dict.items()]) - time = datetime.utcnow() - hist = "\nRan pspec_pipe.py {} step at\nUTC {} with \nhera_pspec [{}], "\ - "and pyuvdata [{}]\nwith {} algorithm "\ - "attrs:\n{}\n{}\n".format(action, time, hp.version.git_hash[:10], - pyuvdata.version.git_hash[:10], action, '-'*50, dict_str) - return hist - -# Create multiprocesses -if params['multiproc']: - pool = multiprocess.Pool(params['nproc']) - M = pool.map -else: - M = map - -# change to working dir -os.chdir(params['work_dir']) - -# out_dir should be cleared before each run: issue a warning if not the case -outdir = os.path.join(params['work_dir'], params['out_dir']) -oldfiles = glob.glob(outdir+"/*") -if len(oldfiles) > 0: - hp.utils.log("\n{}\nWARNING: out_dir should be cleaned before each new run to " \ - "ensure proper functionality.\nIt seems like some files currently " \ - "exist in {}\n{}\n".format('-'*50, outdir, '-'*50), f=lf, verbose=params['verbose']) - -#------------------------------------------------------------------------------- -# Run Visibility Data Difference -#------------------------------------------------------------------------------- -if params['run_diff']: - # start block - time = datetime.utcnow() - hp.utils.log("\n{}\nstarting {} visibility data difference: {}\n".format("-"*60, algs['diff']['diff_type'], time), f=lf, verbose=params['verbose']) - - raise NotImplementedError - - -#------------------------------------------------------------------------------- -# Run OQE Pipeline -#------------------------------------------------------------------------------- -if params['run_pspec']: - # start block - time = datetime.utcnow() - hp.utils.log("\n{}\nstarting pspec OQE: {}\n".format("-"*60, time), f=lf, verbose=params['verbose']) - - # configure dataset groupings - groupings = hp.utils.config_pspec_blpairs(params['data_template'], params['pol_pairs'], params['group_pairs'], exclude_auto_bls=algs['pspec']['exclude_auto_bls'], - exclude_permutations=algs['pspec']['exclude_permutations'], bl_len_range=params['bl_len_range'], - bl_deg_range=params['bl_deg_range'], xants=params['xants']) - - # create dictionary of individual jobs to launch - jobs = odict() - for pair in groupings: - # make Nblps / Nblps_per_job jobs for this pair - blps = groupings[pair] - Nblps = len(blps) - Njobs = Nblps // algs['pspec']['Nblps_per_job'] + 1 - job_blps = [blps[i*algs['pspec']['Nblps_per_job']:(i+1)*algs['pspec']['Nblps_per_job']] for i in range(Njobs)] - for i, _blps in enumerate(job_blps): - key = tuple(hp.utils.flatten(pair) + [i]) - jobs[key] = _blps - - # setup output - outfname = os.path.join(params['out_dir'], algs['pspec']['outfname']) - - # create pspec worker function - def pspec(i, outfname=outfname, jobs=jobs, params=params, alg=algs['pspec'], ef=ef): - try: - # get key - key = jobs.keys()[i] - # parse dsets - dsets = [params['data_template'].format(group=key[0], pol=key[2]), params['data_template'].format(group=key[1], pol=key[3])] - dset_labels = [os.path.basename(d) for d in dsets] - if params['std_template'] is not None: - dsets_std = [params['std_template'].format(group=key[0], pol=key[2]), params['std_template'].format(group=key[1], pol=key[3])] - else: - dsets_std = None - - # pspec_run - hp.pspecdata.pspec_run(dsets, outfname, dsets_std=dsets_std, dset_labels=dset_labels, - dset_pairs=[(0, 1)], spw_ranges=alg['spw_ranges'], n_dlys=alg['n_dlys'], - pol_pairs=params['pol_pairs'], blpairs=jobs[key], input_data_weight=alg['input_data_weight'], - norm=alg['norm'], taper=alg['taper'], beam=alg['beam'], cosmo=alg['cosmo'], - rephase_to_dset=alg['rephase_to_dset'], trim_dset_lsts=alg['trim_dset_lsts'], - broadcast_dset_flags=alg['broadcast_dset_flags'], time_thresh=alg['time_thresh'], - Jy2mK=alg['Jy2mK'], overwrite=params['overwrite'], psname_ext="_{}".format(key[4]), - verbose=params['verbose']) - except: - hp.utils.log("\nPSPEC job {} errored with:".format(i), f=ef, tb=sys.exc_info(), verbose=params['verbose']) - return 1 - - return 0 - - # launch pspec jobs - failures = hp.utils.job_monitor(pspec, range(len(jobs)), "PSPEC", lf=lf, maxiter=algs['pspec']['maxiter'], verbose=params['verbose']) - - # print failures if they exist - if len(failures) > 0: - hp.utils.log("\nSome PSPEC jobs failed after {} tries:\n{}".format(algs['pspec']['maxiter'], '\n'.join(["job {}: {}".format(i, str(jobs.keys()[i])) for i in failures])), f=lf, verbose=params['verbose']) - - # Merge power spectrum files from separate jobs - hp.utils.log("\nStarting power spectrum file merge: {}\n".format(time), f=lf, verbose=params['verbose']) - - # Get all groups - psc = hp.PSpecContainer(outfname, 'r') - groups = psc.groups() - del psc - - # Define merge function - def merge(i, groups=groups, filename=outfname, ef=ef, params=params): - try: - psc = hp.PSpecContainer(filename, mode='rw') - grp = groups[i] - hp.container.combine_psc_spectra(psc, groups=[grp], overwrite=params['overwrite']) - except: - hp.utils.log("\nPSPEC MERGE job {} errored with:".format(i), f=ef, tb=sys.exc_info(), verbose=params['verbose']) - return 1 - - return 0 - - # launch pspec merge jobs - failures = hp.utils.job_monitor(merge, range(len(groups)), "PSPEC MERGE", lf=lf, maxiter=algs['pspec']['maxiter'], verbose=params['verbose']) - - # print failures if they exist - if len(failures) > 0: - hp.utils.log("\nSome PSPEC MERGE jobs failed after {} tries:\n{}".format(algs['pspec']['maxiter'], '\n'.join(["group {}: {}".format(i, str(groups[i])) for i in failures])), f=lf, verbose=params['verbose']) - - # print to log - time = datetime.utcnow() - hp.utils.log("\nFinished PSPEC pipeline: {}\n{}".format(time, "-"*60), f=lf, verbose=params['verbose']) - - -#------------------------------------------------------------------------------- -# Run Bootstrap Pipeline -#------------------------------------------------------------------------------- -if params['run_bootstrap']: - # start block - time = datetime.utcnow() - hp.utils.log("\n{}\nStarting BOOTSTRAP resampling pipeline: {}\n".format("-"*60, time), f=lf, verbose=params['verbose']) - - # ensure outfname is same as pspec - if params['run_pspec'] and (algs['pspec']['outfname'] != algs['bootstrap']['psc_name']): - raise ValueError("bootstrap psc_name {} doesn't equal pspec outfname {}".format(algs['bootstrap']['psc_name'], algs['pspec']['outfname'])) - - # open container - psc_name = os.path.join(params['out_dir'], algs['bootstrap']['psc_name']) - psc = hp.PSpecContainer(psc_name, mode='r') - - # get groups - groups = psc.groups() - del psc - - # define bootstrap function - def bootstrap(i, groups=groups, ef=ef, alg=algs['bootstrap'], params=params, psc_name=psc_name): - try: - # get container - psc = hp.PSpecContainer(psc_name, mode='rw') - - # get spectra - group = groups[i] - spectra = [os.path.join(group, sp) for sp in psc.spectra(group)] - - # run bootstrap - hp.grouping.bootstrap_run(psc, spectra=spectra, time_avg=alg['time_avg'], Nsamples=alg['Nsamples'], - seed=alg['seed'], normal_std=alg['normal_std'], robust_std=alg['robust_std'], - cintervals=alg['cintervals'], keep_samples=alg['keep_samples'], - bl_error_tol=alg['bl_error_tol'], overwrite=params['overwrite'], verbose=params['verbose']) - - - except: - hp.utils.log("\nBOOTSTRAP job {} errored with:".format(i), f=ef, tb=sys.exc_info(), verbose=params['verbose']) - return 1 - - return 0 - - # launch bootstrap jobs - failures = hp.utils.job_monitor(bootstrap, range(len(groups)), "BOOTSTRAP", lf=lf, maxiter=algs['bootstrap']['maxiter'], verbose=params['verbose']) - - # print failures if they exist - if len(failures) > 0: - hp.utils.log("\nSome BOOTSTRAP jobs failed after {} tries:\n{}".format(algs['bootstrap']['maxiter'], '\n'.join(["group {}: {}".format(i, str(groups[i])) for i in failures])), f=lf, verbose=params['verbose']) - - # print to log - time = datetime.utcnow() - hp.utils.log("\nFinished BOOTSTRAP pipeline: {}\n{}".format(time, "-"*60), f=lf, verbose=params['verbose']) - diff --git a/pipelines/pspec_pipeline/pspec_pipe.yaml b/pipelines/pspec_pipeline/pspec_pipe.yaml deleted file mode 100644 index 488a62b2..00000000 --- a/pipelines/pspec_pipeline/pspec_pipe.yaml +++ /dev/null @@ -1,106 +0,0 @@ -# pspec_pipe.yaml -# -# hera_pspec IDR2.1 pipeline -# configuration file -# -# Note: only strings, booleans, integers, floats, -# lists, and lists of lists may be specified here. -# When PSpecData parameters have required type None -# use 'None' here, and for list of tuples, -# use list of lists here. - -#--------------------------------------------------------------- -# IO Parameters -#--------------------------------------------------------------- -io : - work_dir : './' # directory to work in - out_dir : './' # directory for output: root is work_dir - logfile : 'pspipe_out.log' # logfile: root is work_dir - errfile : 'pspipe_err.log' # error file: root is work_dir - joinlog : False # redirect error output into logfile - overwrite : True # overwrite - verbose : True # report feedback to standard output - -#--------------------------------------------------------------- -# Analysis Parameters -#--------------------------------------------------------------- -analysis : - run_diff : False # Run visibility data difference - run_pspec : True # Run power spectrum estimation - run_bootstrap : True # Run bootstrapping for errorbars - multiproc : False # use multiprocess module - nproc : 8 # number of processes to spawn - -#--------------------------------------------------------------- -# Data Parameters -#--------------------------------------------------------------- -data : - # Datafile group pairs to search for - group_pairs : - - ['2458042', '2458042'] - - # Datafile pol pairs to search for - pol_pairs : - - ['xx', 'xx'] - - # data root dir: absolute path - data_root : './' - # data template basename with {pol} and {group} fields - data_template : 'zen.{group}.?????.{pol}.HH.uvXA' - # std template basename with {pol} and {group} fields - std_template : None - - # baseline type restriction in length (meters) and angle (degrees) - bl_len_range : [14.0, 16.0] - bl_deg_range : [0, 10] - - # denote bad antennas - xants : [0, 1, 50, 98] - -#--------------------------------------------------------------- -# Algorithm Parameters -#--------------------------------------------------------------- -algorithm : - - # jacknife data difference - diff : - split_type : "antenna" - - # pspec pipeline - pspec : - outfname : "pspec.h5" - input_data_weight : "identity" - norm : "I" - taper : "none" - beam : "HERA_NF_dipole_power.beamfits" # work_dir is root - cosmo : "None" - spw_ranges : - - [10, 20] - n_dlys : - - 5 - Jy2mK : True - rephase_to_dset : None - trim_dset_lsts : False - broadcast_dset_flags : True - time_thresh : 0.2 - Nblps_per_job : 200 - exclude_auto_bls : True - exclude_permutations : True - bl_len_range : [0, 20] - bl_deg_range : [40, 80] - maxiter : 1 - - # bootstrap pipeline - bootstrap : - psc_name : "pspec.h5" - time_avg : True - Nsamples : 100 - seed : None - normal_std : True - robust_std : True - cintervals : None - keep_samples : False - bl_error_tol : 1.0 - maxiter : 1 - - diff --git a/setup.py b/setup.py index 7df503be..6af90132 100644 --- a/setup.py +++ b/setup.py @@ -33,8 +33,6 @@ def package_files(package_dir, subdirectory): 'install_requires': ['numpy>=1.14', 'scipy','matplotlib>=2.2'], 'include_package_data': True, 'scripts': ['scripts/pspec_run.py', 'scripts/pspec_red.py', - 'pipelines/idr2_preprocessing/preprocess_data.py', - 'pipelines/pspec_pipeline/pspec_pipe.py', 'scripts/bootstrap_run.py'], 'zip_safe': False, } From 96fe18f26e95a912108984c79933ec767c7251de Mon Sep 17 00:00:00 2001 From: Nicholas Kern Date: Tue, 11 Sep 2018 00:33:57 +0800 Subject: [PATCH 9/9] ignore deltasq if delay is True in plot.delay_wedge modified: plot.py --- hera_pspec/plot.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/hera_pspec/plot.py b/hera_pspec/plot.py index 69e0c865..1dfeb1b5 100644 --- a/hera_pspec/plot.py +++ b/hera_pspec/plot.py @@ -666,7 +666,8 @@ def delay_wedge(uvp, spw, pol, blpairs=None, times=None, fold=False, delay=True, If ax is None, this is the new figure size. deltasq : bool, optional - Convert to Delta^2 before plotting. Default: False + Convert to Delta^2 before plotting. This is ignored if delay=True. + Default: False colorbar : bool, optional Add a colorbar to the plot. Default: False @@ -727,7 +728,7 @@ def delay_wedge(uvp, spw, pol, blpairs=None, times=None, fold=False, delay=True, uvp.average_spectra(blpair_groups=blp_grps, time_avg=True, inplace=True) # Convert to DeltaSq - if deltasq: + if deltasq and not delay: uvp.convert_to_deltasq(inplace=True) # Fold array