From e31c8f27e8e1ac8452d4ea6d42e2cb619811e4b3 Mon Sep 17 00:00:00 2001 From: Nicholas Kern Date: Fri, 24 May 2019 18:22:46 -0700 Subject: [PATCH 1/5] added time-dependent PN in generate_noise modified: hera_pspec/grouping.py modified: hera_pspec/tests/test_grouping.py modified: hera_pspec/tests/test_uvpspec.py modified: hera_pspec/uvpspec.py --- hera_pspec/grouping.py | 1 - hera_pspec/tests/test_grouping.py | 20 ++++++++++++-------- hera_pspec/tests/test_uvpspec.py | 2 ++ hera_pspec/uvpspec.py | 13 ++++++++++--- 4 files changed, 24 insertions(+), 12 deletions(-) diff --git a/hera_pspec/grouping.py b/hera_pspec/grouping.py index 2350296b..4b3ee14a 100644 --- a/hera_pspec/grouping.py +++ b/hera_pspec/grouping.py @@ -334,7 +334,6 @@ def average_spectra(uvp_in, blpair_groups=None, time_avg=False, if store_cov: bpg_cov.append(cov * w[:, None]) - # Take integration-weighted averages, with clipping to deal # with zeros bpg_data = np.sum(bpg_data, axis=0) \ diff --git a/hera_pspec/tests/test_grouping.py b/hera_pspec/tests/test_grouping.py index 72f4b46d..a4838648 100644 --- a/hera_pspec/tests/test_grouping.py +++ b/hera_pspec/tests/test_grouping.py @@ -215,16 +215,20 @@ def test_validate_bootstrap_errorbar(): taper='none', sampling=False, little_h=False, spw_ranges=[(0, 50)], verbose=False) # bootstrap resample - seed = 0 - Nsamples = 200 - uvp_avg, uvp_boots, uvp_wgts = grouping.bootstrap_resampled_error(uvp, time_avg=False, Nsamples=Nsamples, - seed=seed, normal_std=True, - blpair_groups=[uvp.get_blpairs()]) + # uvp.data_array[0] = (stats.norm.rvs(0,1,450*50)+1j*stats.norm.rvs(0,1,450*50)).reshape(450,50,1) + bs_std_zscr = [] + for seed in range(40): + Nsamples = 50 + uvp_avg, uvp_boots, uvp_wgts = grouping.bootstrap_resampled_error(uvp, time_avg=False, Nsamples=Nsamples, + seed=seed, normal_std=True, + blpair_groups=[uvp.get_blpairs()]) + + # assert z-score has std of ~1.0 along time ax to within 1/sqrt(Nsamples) + bs_std_zscr_real = np.sqrt(np.mean(uvp_avg.data_array[0].real**2)) / np.mean(uvp_avg.stats_array['bs_std'][0].real) + bs_std_zscr_imag = np.sqrt(np.mean(uvp_avg.data_array[0].imag**2)) / np.mean(uvp_avg.stats_array['bs_std'][0].imag) + bs_std_zscr.append(bs_std_zscr_real + 1j * bs_std_zscr_imag) - # assert z-score has std of ~1.0 along time ax to within 1/sqrt(Nsamples) - bs_std_zscr_real = np.std(uvp_avg.data_array[0].real) / np.mean(uvp_avg.stats_array['bs_std'][0].real) nt.assert_true(np.abs(1.0 - bs_std_zscr_real) < 1/np.sqrt(Nsamples)) - bs_std_zscr_imag = np.std(uvp_avg.data_array[0].imag) / np.mean(uvp_avg.stats_array['bs_std'][0].imag) nt.assert_true(np.abs(1.0 - bs_std_zscr_imag) < 1/np.sqrt(Nsamples)) diff --git a/hera_pspec/tests/test_uvpspec.py b/hera_pspec/tests/test_uvpspec.py index 7397af1f..fb6d84b2 100644 --- a/hera_pspec/tests/test_uvpspec.py +++ b/hera_pspec/tests/test_uvpspec.py @@ -354,6 +354,8 @@ def test_sense(self): # test a blpair selection P_N = uvp.generate_noise_spectra(0, 1515, 500, form='Pk', component='real') + # test as a dictionary + def test_average_spectra(self): uvp = copy.deepcopy(self.uvp) # test blpair averaging diff --git a/hera_pspec/uvpspec.py b/hera_pspec/uvpspec.py index c6660dd7..b77622ad 100644 --- a/hera_pspec/uvpspec.py +++ b/hera_pspec/uvpspec.py @@ -1608,8 +1608,9 @@ def generate_noise_spectra(self, spw, polpair, Tsys, blpairs=None, polpair int. Strings are expanded into polarization pairs, e.g. 'XX' becomes ('XX,'XX'). - Tsys : float - System temperature in Kelvin. + Tsys : dictionary, float or array + System temperature in Kelvin for each blpair. Key is blpair-integer, + value is Tsys float or ndarray. If fed as an ndarray, shape=(Ntimes,) blpairs : list List of unique blair tuples or i12 integers to calculate noise @@ -1667,6 +1668,12 @@ def generate_noise_spectra(self, spw, polpair, Tsys, blpairs=None, # Get delays dlys = self.get_dlys(spw) + # handle Tsys + if not isinstance(Tsys, (dict, odict)): + if not isinstance(Tsys, np.ndarray): + Tsys = np.ones(self.Ntimes) * Tsys + Tsys = dict([(blp, Tsys) for blp in blpairs]) + # Iterate over blpairs to get P_N P_N = odict() for i, blp in enumerate(blpairs): @@ -1687,7 +1694,7 @@ def generate_noise_spectra(self, spw, polpair, Tsys, blpairs=None, k = None # Get noise power spectrum - pn = noise.calc_P_N(scalar, Tsys, t_int, k=k, + pn = noise.calc_P_N(scalar, Tsys[blp][j], t_int, k=k, Nincoherent=n_samp, form=form, component=component) # Put into appropriate form From 643cc07bcf7a0efc0f369ec2fc29db07bd6ac7cb Mon Sep 17 00:00:00 2001 From: Nicholas Kern Date: Fri, 24 May 2019 18:41:43 -0700 Subject: [PATCH 2/5] added uvpspec PN test --- hera_pspec/tests/test_uvpspec.py | 21 ++++++++++++++------- 1 file changed, 14 insertions(+), 7 deletions(-) diff --git a/hera_pspec/tests/test_uvpspec.py b/hera_pspec/tests/test_uvpspec.py index fb6d84b2..0eca35cb 100644 --- a/hera_pspec/tests/test_uvpspec.py +++ b/hera_pspec/tests/test_uvpspec.py @@ -335,26 +335,33 @@ def test_sense(self): uvp = copy.deepcopy(self.uvp) # test generate noise spectra - P_N = uvp.generate_noise_spectra(0, 1515, 500, form='Pk', component='real') + polpair = ('xx', 'xx') + P_N = uvp.generate_noise_spectra(0, polpair, 500, form='Pk', component='real') nt.assert_equal(P_N[101102101102].shape, (10, 30)) # test smaller system temp - P_N2 = uvp.generate_noise_spectra(0, 1515, 400, form='Pk', component='real') + P_N2 = uvp.generate_noise_spectra(0, polpair, 400, form='Pk', component='real') nt.assert_true((P_N[101102101102] > P_N2[101102101102]).all()) # test complex - P_N2 = uvp.generate_noise_spectra(0, 1515, 500, form='Pk', component='abs') + P_N2 = uvp.generate_noise_spectra(0, polpair, 500, form='Pk', component='abs') nt.assert_true((P_N[101102101102] < P_N2[101102101102]).all()) # test Dsq - Dsq = uvp.generate_noise_spectra(0, 1515, 500, form='DelSq', component='real') + Dsq = uvp.generate_noise_spectra(0, polpair, 500, form='DelSq', component='real') nt.assert_equal(Dsq[101102101102].shape, (10, 30)) nt.assert_true(Dsq[101102101102][0, 1] < P_N[101102101102][0, 1]) - # test a blpair selection - P_N = uvp.generate_noise_spectra(0, 1515, 500, form='Pk', component='real') + # test a blpair selection and int polpair + blpairs = uvp.get_blpairs()[:1] + P_N = uvp.generate_noise_spectra(0, 1515, 500, form='Pk', blpairs=blpairs, component='real') + nt.assert_equal(P_N[101102101102].shape, (10, 30)) - # test as a dictionary + # test as a dictionary of arrays + Tsys = dict([(uvp.antnums_to_blpair(k), 500 * np.ones((uvp.Ntimes, uvp.Ndlys)) * np.linspace(1, 2, uvp.Ntimes)[:, None]) for k in uvp.get_blpairs()]) + P_N = uvp.generate_noise_spectra(0, 1515, Tsys, form='Pk', blpairs=blpairs, component='real') + # assert time gradient is captured: 2 * Tsys results in 4 * P_N + nt.assert_true(np.isclose(P_N[101102101102][0, 0] * 4, P_N[101102101102][-1, 0])) def test_average_spectra(self): uvp = copy.deepcopy(self.uvp) From d4e7a6437d0a1fd19d209cbf6ad95699be2990da Mon Sep 17 00:00:00 2001 From: Nicholas Kern Date: Mon, 27 May 2019 11:29:35 -0700 Subject: [PATCH 3/5] modified: hera_pspec/uvpspec.py --- hera_pspec/uvpspec.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/hera_pspec/uvpspec.py b/hera_pspec/uvpspec.py index b77622ad..fefeced5 100644 --- a/hera_pspec/uvpspec.py +++ b/hera_pspec/uvpspec.py @@ -1679,7 +1679,7 @@ def generate_noise_spectra(self, spw, polpair, Tsys, blpairs=None, for i, blp in enumerate(blpairs): # get indices inds = self.blpair_to_indices(blp) - + assert isinstance(Tsys[blp], (float, np.float, int, np.int)) or Tsys[blp].shape[0] == self.Ntimes, "Tsys must be a float or an ndarray with shape[0] == Ntimes" P_blp = [] # iterate over time axis for j, ind in enumerate(inds): From 14755ec4282ad4312992d5f0b8a04f59d6dae03e Mon Sep 17 00:00:00 2001 From: Nicholas Kern Date: Wed, 29 May 2019 13:57:20 -0700 Subject: [PATCH 4/5] modified: hera_pspec/noise.py --- hera_pspec/noise.py | 7 ++----- 1 file changed, 2 insertions(+), 5 deletions(-) diff --git a/hera_pspec/noise.py b/hera_pspec/noise.py index 03825fbc..27fb9195 100644 --- a/hera_pspec/noise.py +++ b/hera_pspec/noise.py @@ -44,11 +44,8 @@ def calc_P_N(scalar, Tsys, t_int, Ncoherent=1, Nincoherent=None, form='Pk', k=No assert form in ('Pk', 'DelSq'), "form must be either 'Pk' or 'DelSq' for P(k) or Delta^2(k) respectively" assert component in ['abs', 'real', 'imag'], "component must be one of 'real', 'imag', 'abs'" - # convert to mK - Tsys *= 1e3 - - # construct prefactor - P_N = scalar * Tsys**2 + # construct prefactor in mK^2 + P_N = scalar * (Tsys * 1e3)**2 # Multiply in effective integration time P_N /= (t_int * Ncoherent) From de7da81252e751cacd93397493853a914d621406 Mon Sep 17 00:00:00 2001 From: Nicholas Kern Date: Sun, 9 Jun 2019 14:22:56 -0700 Subject: [PATCH 5/5] update zscr in grouping bs validation test --- hera_pspec/tests/test_grouping.py | 25 ++++++++++--------------- hera_pspec/uvpspec.py | 6 +++--- 2 files changed, 13 insertions(+), 18 deletions(-) diff --git a/hera_pspec/tests/test_grouping.py b/hera_pspec/tests/test_grouping.py index a4838648..9caaf681 100644 --- a/hera_pspec/tests/test_grouping.py +++ b/hera_pspec/tests/test_grouping.py @@ -215,21 +215,16 @@ def test_validate_bootstrap_errorbar(): taper='none', sampling=False, little_h=False, spw_ranges=[(0, 50)], verbose=False) # bootstrap resample - # uvp.data_array[0] = (stats.norm.rvs(0,1,450*50)+1j*stats.norm.rvs(0,1,450*50)).reshape(450,50,1) - bs_std_zscr = [] - for seed in range(40): - Nsamples = 50 - uvp_avg, uvp_boots, uvp_wgts = grouping.bootstrap_resampled_error(uvp, time_avg=False, Nsamples=Nsamples, - seed=seed, normal_std=True, - blpair_groups=[uvp.get_blpairs()]) - - # assert z-score has std of ~1.0 along time ax to within 1/sqrt(Nsamples) - bs_std_zscr_real = np.sqrt(np.mean(uvp_avg.data_array[0].real**2)) / np.mean(uvp_avg.stats_array['bs_std'][0].real) - bs_std_zscr_imag = np.sqrt(np.mean(uvp_avg.data_array[0].imag**2)) / np.mean(uvp_avg.stats_array['bs_std'][0].imag) - bs_std_zscr.append(bs_std_zscr_real + 1j * bs_std_zscr_imag) - - nt.assert_true(np.abs(1.0 - bs_std_zscr_real) < 1/np.sqrt(Nsamples)) - nt.assert_true(np.abs(1.0 - bs_std_zscr_imag) < 1/np.sqrt(Nsamples)) + Nsamples = 1000 + seed = 0 + uvp_avg, uvp_boots, uvp_wgts = grouping.bootstrap_resampled_error(uvp, time_avg=False, Nsamples=Nsamples, + seed=seed, normal_std=True, + blpair_groups=[uvp.get_blpairs()]) + # assert z-score has std of ~1.0 along time ax to within 1/sqrt(Nsamples) + zscr_real = np.std(uvp_avg.data_array[0].real / uvp_avg.stats_array['bs_std'][0].real) + zscr_imag = np.std(uvp_avg.data_array[0].imag / uvp_avg.stats_array['bs_std'][0].imag) + nt.assert_true(np.abs(1.0 - zscr_real) < 1/np.sqrt(Nsamples)) + nt.assert_true(np.abs(1.0 - zscr_imag) < 1/np.sqrt(Nsamples)) def test_bootstrap_run(): diff --git a/hera_pspec/uvpspec.py b/hera_pspec/uvpspec.py index fefeced5..6131cf96 100644 --- a/hera_pspec/uvpspec.py +++ b/hera_pspec/uvpspec.py @@ -1572,8 +1572,8 @@ def generate_noise_spectra(self, spw, polpair, Tsys, blpairs=None, little_h=True, form='Pk', num_steps=2000, component='real'): """ - Generate the expected 1-sigma noise power spectrum given a selection of - spectral window, system temp., and polarization. This estimate is + Generate the expected RMS noise power spectrum given a selection of + spectral window, system temp. [K], and polarization. This estimate is constructed as: P_N = scalar * (Tsys * 1e3)^2 / (integration_time) / sqrt(Nincoherent) @@ -1588,7 +1588,7 @@ def generate_noise_spectra(self, spw, polpair, Tsys, blpairs=None, If the polarizations specified are pseudo Stokes pol (I, Q, U or V) then an extra factor of 2 is divided. - If form == 'DelSq' then a factor of k^3 / (2pi^2) is multiplied. + If form == 'DelSq' then a factor of |k|^3 / (2pi^2) is multiplied. If real is True, a factor of sqrt(2) is divided to account for discarding imaginary noise component.