Skip to content

Commit

Permalink
eliminated uvpspec.generate_noise pstokes correction
Browse files Browse the repository at this point in the history
and added noise_validation test to test_noise
  • Loading branch information
nkern committed Aug 18, 2018
1 parent e4fc532 commit 45a0e2d
Show file tree
Hide file tree
Showing 3 changed files with 51 additions and 10 deletions.
4 changes: 2 additions & 2 deletions hera_pspec/testing.py
Original file line number Diff line number Diff line change
Expand Up @@ -201,8 +201,8 @@ def uvpspec_from_data(data, bl_grps, data_std=None, spw_ranges=None, beam=None,
return uvp


def noise_sim(data, Tsys, beam, Nextend=0, seed=None, inplace=False, whiten=False,
run_check=True):
def noise_sim(data, Tsys, beam, Nextend=0, seed=None, inplace=False,
whiten=False, run_check=True):
"""
Generate a simulated Gaussian noise realization.
Expand Down
48 changes: 47 additions & 1 deletion hera_pspec/tests/test_noise.py
Original file line number Diff line number Diff line change
Expand Up @@ -4,7 +4,7 @@
import os
import sys
from hera_pspec.data import DATA_PATH
from hera_pspec import uvpspec, conversions, parameter, pspecbeam, noise
from hera_pspec import uvpspec, conversions, pspecdata, pspecbeam, noise, testing, utils
import copy
import h5py
from collections import OrderedDict as odict
Expand Down Expand Up @@ -64,4 +64,50 @@ def test_calc_P_N(self):
nt.assert_true(Dsq[1] < P_N)


def test_noise_validation():
"""
make sure that the noise.py code produces
correct noise 1-sigma amplitude using a
noise simulation.
"""
# get simulated noise in Jy
bfile = os.path.join(DATA_PATH, 'HERA_NF_dipole_power.beamfits')
beam = pspecbeam.PSpecBeamUV(bfile)
uvfile = os.path.join(DATA_PATH, "zen.even.xx.LST.1.28828.uvOCRSA")
Tsys = 300.0 # Kelvin

# generate noise
seed = 0
uvd = testing.noise_sim(uvfile, Tsys, beam, seed=seed, whiten=True, inplace=False, Nextend=9)

# get redundant baseline group
reds, lens, angs = utils.get_reds(uvd, pick_data_ants=True, bl_len_range=(10, 20),
bl_deg_range=(0, 1))
bls1, bls2, blps = utils.construct_blpairs(reds[0], exclude_auto_bls=True, exclude_permutations=True)

# setup PSpecData
ds = pspecdata.PSpecData(dsets=[copy.deepcopy(uvd), copy.deepcopy(uvd)], wgts=[None, None], beam=beam)
ds.Jy_to_mK()

# get pspec
uvp = ds.pspec(bls1, bls2, (0, 1), [('xx', 'xx')], input_data_weight='identity', norm='I',
taper='none', sampling=False, little_h=True, spw_ranges=[(0, 50)], verbose=False)

# get noise spectra from one of the blpairs
P_N = uvp.generate_noise_spectra(0, 'xx', Tsys, blpairs=uvp.get_blpairs()[:1], num_steps=2000).values()[0][0, 0]

# get P_std of real spectra for each baseline across time axis
P_stds = np.array([np.std(uvp.get_data((0, bl, 'xx')).real, axis=1) for bl in uvp.get_blpairs()])

# get average P_std_avg and its standard error
P_std_avg = np.mean(P_stds)

# assert close to P_N: 2%
# This should be updated to be within standard error on P_std_avg
# when the spw_range-variable pspec amplitude bug is resolved
nt.assert_true(np.abs(P_std_avg - P_N) / P_N < 0.02)





9 changes: 2 additions & 7 deletions hera_pspec/uvpspec.py
Original file line number Diff line number Diff line change
Expand Up @@ -1570,13 +1570,8 @@ def generate_noise_spectra(self, spw, pol, Tsys, blpairs=None, little_h=True,
if form == 'Pk':
pn = np.ones(len(dlys), np.float) * pn

# If pseudo stokes pol (as opposed to linear or circular pol),
# divide by extra factor of 2
if isinstance(pol, (np.str, str)):
pol = uvutils.polstr2num(pol)

if pol in (1, 2, 3, 4): pn /= 2.0 # pseudo stokes pol
if real: pn /= np.sqrt(2) # if real divide by sqrt(2)
if real:
pn /= np.sqrt(2) # if real divide by sqrt(2)

# append to P_blp
P_blp.append(pn)
Expand Down

0 comments on commit 45a0e2d

Please sign in to comment.