Merge e31c8f2 into 3132ad5

HERA-Team · May 25, 2019 · bdd451c · bdd451c
2 parents 3132ad5 + e31c8f2
commit bdd451c
Show file tree

Hide file tree

Showing 4 changed files with 24 additions and 12 deletions.
diff --git 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
@@ -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
@@ -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
@@ -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