diff --git a/.gitignore b/.gitignore new file mode 100644 index 00000000..73eaaff4 --- /dev/null +++ b/.gitignore @@ -0,0 +1,12 @@ +# compiled code +*.pyc +*__pycache__/ + +# build +build/ +env/ +dist/ +eggs/ +*.eggs/ + + diff --git a/.travis.yml b/.travis.yml new file mode 100644 index 00000000..c83402d5 --- /dev/null +++ b/.travis.yml @@ -0,0 +1,43 @@ +language: python +python: + # We don't actually use the Travis Python, but this keeps it organized. + - "2.7" +install: + # ensure that we have the full tag information available for version.py + - git fetch --tags + - sudo apt-get update + # We do this conditionally because it saves us some downloading if the + # version is the same. + - if [[ "$TRAVIS_PYTHON_VERSION" == "2.7" ]]; then + wget https://repo.continuum.io/miniconda/Miniconda2-latest-Linux-x86_64.sh -O miniconda.sh; + else + wget https://repo.continuum.io/miniconda/Miniconda3-latest-Linux-x86_64.sh -O miniconda.sh; + fi + - bash miniconda.sh -b -p $HOME/miniconda + - export PATH="$HOME/miniconda/bin:$PATH" + - hash -r + - conda config --set always_yes yes --set changeps1 no + - conda install -q conda=4.3.25 + # Useful for debugging any issues with conda + - conda info -a + + # 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 + - pip install coveralls + - pip install git+https://github.com/HERA-Team/pyuvdata.git + - pip install git+https://github.com/HERA-Team/omnical.git + - pip install git+https://github.com/HERA-Team/hera_cal.git + - pip install git+https://github.com/HERA-Team/linsolve.git + - python setup.py install + +before_script: + - "export DISPLAY=:99.0" + - "sh -e /etc/init.d/xvfb start" + - sleep 3 + +script: nosetests hera_pspec --with-coverage --cover-package=hera_pspec + +after_success: + - coveralls diff --git a/LICENSE b/LICENSE new file mode 100644 index 00000000..c1f3eec7 --- /dev/null +++ b/LICENSE @@ -0,0 +1,29 @@ +BSD 3-Clause License + +Copyright (c) 2017, Hydrogen Epoch of Reionization Array (HERA) +All rights reserved. + +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions are met: + +* Redistributions of source code must retain the above copyright notice, this + list of conditions and the following disclaimer. + +* Redistributions in binary form must reproduce the above copyright notice, + this list of conditions and the following disclaimer in the documentation + and/or other materials provided with the distribution. + +* Neither the name of the copyright holder nor the names of its + contributors may be used to endorse or promote products derived from + this software without specific prior written permission. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE +DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE +FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL +DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR +SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER +CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, +OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE +OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. diff --git a/MANIFEST.in b/MANIFEST.in new file mode 100644 index 00000000..02cd31e4 --- /dev/null +++ b/MANIFEST.in @@ -0,0 +1,5 @@ +include *.md +include LICENSE +include hera_pspec/VERSION +include hera_pspec/GIT_INFO + diff --git a/README.md b/README.md index f1458715..890a1dc1 100644 --- a/README.md +++ b/README.md @@ -1 +1,28 @@ -# hera_pspec \ No newline at end of file +# hera_pspec + +[![Build Status](https://travis-ci.org/HERA-Team/hera_pspec.svg?branch=master)](https://travis-ci.org/HERA-Team/hera_pspec) + +**hera_pspec description** + +## Installation + +### Code Dependencies + +* numpy >= 1.10 +* pyuvdata (`pip install pyuvdata` or use https://github.com/HERA-Team/pyuvdata.git) +* aipy (```conda install -c conda-forge aipy```) +* scipy >= 0.19 +* astropy >= 2.0 + +For anaconda users, we suggest using conda to install astropy, numpy and scipy. + +### Installing hera_pspec +Clone the repo using +`git clone https://github.com/HERA-Team/hera_pspec.git` + +Navigate into the directory and run `python setup.py install`. + +## Package Details and Usage + + + diff --git a/hera_pspec/VERSION b/hera_pspec/VERSION new file mode 100644 index 00000000..77d6f4ca --- /dev/null +++ b/hera_pspec/VERSION @@ -0,0 +1 @@ +0.0.0 diff --git a/hera_pspec/__init__.py b/hera_pspec/__init__.py index 8acc066f..2e5b950f 100644 --- a/hera_pspec/__init__.py +++ b/hera_pspec/__init__.py @@ -1 +1,9 @@ -import pspec +""" +__init__.py file for hera_pspec +""" +import version +import conversions +from dataset import DataSet +import pspec as legacy # XXX: This will eventually be deprecated + +__version__ = version.version diff --git a/hera_pspec/conversions.py b/hera_pspec/conversions.py new file mode 100644 index 00000000..9ddfdf92 --- /dev/null +++ b/hera_pspec/conversions.py @@ -0,0 +1,256 @@ +""" +conversions.py +-------- + +Cosmological and instrumental +conversion functions for hera_pspec + + +""" +import numpy as np +from scipy import integrate +try: + from astropy.cosmology import LambdaCDM + _astropy = True +except: + print("could not import astropy") + _astropy = False + + +class units: + """ + fundamental constants and conversion constants + in ** SI ** units + + c : speed of light m s-1 + ckm : speed of light in km s-1 + kb : boltzmann constant 2 kg s-2 K-1 + hp : planck constant m2 kg s-1 + sb : stefan-boltzmann constant W m-2 K-4 + f21 : 21cm frequency in Hz + w21 : 21cm wavelength in meters + """ + c = 2.99792458e8 # speed of light in m s-1 + ckm = c / 1e3 # speed of light in km s-1 + G = 6.67408e-11 # Newton constant m3 kg-1 s-2 + kb = 1.38064852e-23 # boltzmann constant m2 kg s-2 K-1 + hp = 6.62607004e-34 # planck constant m2 kg s-1 + sb = 5.670367e-8 # stefan boltzmann constant W m-2 K-4 + f21 = 1.420405751e9 # frequency of 21cm transition in Hz + w21 = 0.211061140542 # 21cm wavelength in meters + H0_to_SI = 3.24078e-20 # km s-1 Mpc-1 to s-1 + +class Cosmo_Conversions(object): + """ + Cosmo_Conversions class for mathematical conversion functions, + some of which require a cosmological model, others of which + do not require a predefined cosmology. + + Default parameter values are Planck 2015 TT,TE,EE+lowP. + (Table 4 of https://doi.org/10.1051/0004-6361/201525830) + + Distance measures come from: + Hogg 1999 (astro-ph/9905116) + Furlanetto 2006 (2006PhR...433..181F) + + Note that all distance measures are by default in h-1 Mpc, + because the default H0 = 100 km / sec / Mpc. If the user + changes the input H0, distance measures are in Mpc. + """ + # astropy load attribute + _astropy = _astropy + + def __init__(self, Om_L=0.68440, Om_b=0.04911, Om_c=0.26442, H0=100.0, + Om_M=None, Om_k=None): + """ + + + Default parameter values are Planck 2015 TT,TE,EE+lowP. + (Table 4 of https://doi.org/10.1051/0004-6361/201525830) + + Parameters: + ----------- + Om_L : float, Omega Lambda naught, default=0.68440 + cosmological constant energy density fraction at z=0 + + Om_b : float, Omega baryon naught, default=0.04911 + baryon energy density fraction at z=0 + + Om_c : float, Omega cold dark matter naught, default=0.26442 + cdm energy density fraction at z=0 + + H0 : float, Hubble constant, default=67.31 km s-1 Mpc-1 + + Om_M : float, Omega matter naught, default=None [optional] + matter energy density faction at z=0 + + Om_k : float, Omega curvature naught, default=None [optional] + curvature energy density at z=0 + + Notes: + ------ + Note that all distance measures are by default in h-1 Mpc, + because the default H0 = 100 km / sec / Mpc. If the user + changes the input H0, the distance measures are in Mpc. + """ + # Setup parameters + if Om_M is not None: + if np.isclose(Om_b + Om_c, Om_M, atol=1e-5) is False: + Om_b = Om_M * 0.156635 + Om_c = Om_M * 0.843364 + else: + Om_M = Om_b + Om_c + + if Om_k is None: + Om_k = 1 - Om_L - Om_M + + ### TODO add radiation component to class and to distance functions + self.Om_L = Om_L + self.Om_b = Om_b + self.Om_c = Om_c + self.Om_M = Om_M + self.Om_k = Om_k + self.H0 = H0 + self.h = self.H0 / 100.0 + + if _astropy: + self.lcdm = LambdaCDM(H0=self.H0, Om0=self.Om_M, Ode0=self.Om_L) + + def f2z(self, freq, ghz=False): + """ + convert frequency to redshift for 21cm line + + Parameters: + ----------- + freq : frequency in Hz, type=float + + ghz : boolean, if True: assume freq is GHz + + Output: + ------- + z : float + redshift + """ + if ghz: + freq = freq * 1e9 + + return (units.f21 / freq - 1) + + def z2f(self, z, ghz=False): + """ + convert redshift to frequency in Hz for 21cm line + + Parameters: + ----------- + z : redshift, type=float + + ghz: boolean, if True: convert to GHz + + Output: + ------- + freq : float + frequency in Hz + """ + freq = units.f21 / (z + 1) + if ghz: + freq /= 1e9 + + return freq + + def E(self, z): + """ + ratio of hubble parameters: H(z) / H(z=0) + Hogg99 Eqn. 14 + + Parameters: + ----------- + z : redshift, type=float + """ + return np.sqrt(self.Om_M*(1+z)**3 + self.Om_k*(1+z)**2 + self.Om_L) + + def DC(self, z): + """ + line-of-sight comoving distance in Mpc + Hogg99 Eqn. 15 + + Parameters: + ----------- + z : redshift, type=float + """ + return integrate.quad(lambda z: 1/self.E(z), 0, z)[0] * units.ckm / self.H0 + + def DM(self, z): + """ + transverse comoving distance in Mpc + Hogg99 Eqn. 16 + + Parameters: + ----------- + z : redshift, type=float + """ + DH = units.ckm / self.H0 + if self.Om_k > 0: + DM = DH * np.sinh(np.sqrt(self.Om_k) * self.DC(z) / DH) / np.sqrt(self.Om_k) + elif self.Om_k < 0: + DM = DH * np.sin(np.sqrt(np.abs(self.Om_k)) * self.DC(z) / DH) / np.sqrt(np.abs(self.Om_k)) + else: + DM = self.DC(z) + + return DM + + def DA(self, z): + """ + angular diameter (proper) distance in Mpc + Hogg99 Eqn. 18 + + Parameters: + ----------- + z : redshift, type=float + """ + return self.DM(z) / (1 + z) + + def dRperp_dtheta(self, z): + """ + conversion factor from angular size (radian) to transverse + comoving distance (Mpc) at a specific redshift: [Mpc / radians] + + Parameters: + ----------- + z : float, redshift + """ + return self.DM(z) + + def dRpara_df(self, z, ghz=False): + """ + conversion from frequency bandwidth to radial + comoving distance at a specific redshift: [Mpc / Hz] + + Parameters: + ----------- + z : float, redshift + ghz : convert output to [Mpc / GHz] + """ + y = (1 + z)**2.0 / self.E(z) * units.ckm / self.H0 / units.f21 + if ghz: + return y * 1e9 + else: + return y + + def X2Y(self, z): + """ + Conversion from radians^2 Hz -> Mpc^3 + at a specific redshift. + + Parameters: + ----------- + z : float, redshift + + Notes: + ------ + Calls Cosmo_Conversions.dRperp_dtheta() and Cosmo_Conversions.dRpara_df(). + """ + return self.dRperp_dtheta(z)**2 * self.dRpara_df(z) + + + + diff --git a/hera_pspec/data/__init__.py b/hera_pspec/data/__init__.py new file mode 100644 index 00000000..22c9eda2 --- /dev/null +++ b/hera_pspec/data/__init__.py @@ -0,0 +1,2 @@ +"""init file for hera_pspec/data/""" +DATA_PATH = __path__[0] diff --git a/hera_pspec/dataset.py b/hera_pspec/dataset.py new file mode 100644 index 00000000..abdc31eb --- /dev/null +++ b/hera_pspec/dataset.py @@ -0,0 +1,771 @@ +import numpy as np +import aipy +from .utils import hash, cov + + +class DataSet(object): + + def __init__(self, dsets=None, wgts=None, lsts=None, conj=None, + lmin=None, lmode=None): + """ + Object to store a set of visibilities and perform operations such as + power spectrum estimation on them. + + Parameters + ---------- + dsets : TODO, optional + TODO. Default: None. + wgts : TODO, optional + TODO. Default: None. + lsts : TODO, optional + TODO. Default: None. + conj : TODO, optional + TODO. Default: None. + lmin : TODO, optional + TODO. Default: None. + lmode : TODO, optional + TODO. Default: None. + """ + self.x, self.w = {}, {} + self.clear_cache() + self.lmin = lmin + self.lmode = lmode + #if not npzfile is None: self.from_npz(npzfile) + if dsets is not None: + self.add_data(dsets, wgts=wgts, conj=conj, replace=True) + + def flatten_data(self, data): + """ + TODO. + + Parameters + ---------- + data : TODO + TODO. + + Returns + ------- + d : TODO + TODO. + """ + if data is None: return None + d = {} + for k in data: + for bl in data[k]: + for pol in data[k][bl]: + key = (k, bl, pol) + d[key] = data[k][bl][pol] + return d + + def set_data(self, dsets, wgts=None, conj=None): + """ + Deprecated. Use add_data() instead. + """ + DeprecationWarning("This function will be removed. Use add_data() instead.") + """ + if type(dsets.values()[0]) == dict: + dsets,wgts = self.flatten_data(dsets), self.flatten_data(wgts) + self.x, self.w = {}, {} + for k in dsets: + self.x[k] = dsets[k].T + try: self.w[k] = wgts[k].T + except(TypeError): self.w[k] = np.ones_like(self.x[k]) + try: + if conj[k[1]]: self.x[k] = np.conj(self.x[k]) + except(TypeError,KeyError): pass + """ + + def add_data(self, dsets, wgts=None, conj=None, replace=False): + """ + TODO. + + Parameters + ---------- + dsets : TODO + TODO. + + wgts : TODO, optional + TODO. Default: None. + + conj : TODO, optional + TODO. Default: None. + + replace : bool, optional + Whether to completely replace the x, w data, or just add to it. + """ + if type(dsets.values()[0]) == dict: + dsets,wgts = self.flatten_data(dsets), self.flatten_data(wgts) + if replace: + self.x, self.w = {}, {} + + for k in dsets: + self.x[k] = dsets[k].T + try: self.w[k] = wgts[k].T + except(TypeError): self.w[k] = np.ones_like(self.x[k]) + try: + if conj[k[1]]: self.x[k] = np.conj(self.x[k]) + except(TypeError, KeyError): pass + + def clear_cache(self, keys=None): + """ + Clear stored covariance data (or some subset of it). + + Parameters + ---------- + keys : TODO, optional + TODO. Default: None. + """ + # XXX right now clear_cache munges I if I != ones on the diagnal. rethink I or clear_cache of I + if keys is None: + self._C, self._Cempirical, self._I, self._iC = {}, {}, {}, {} + self._iCt = {} + else: + for k in keys: + try: del(self._C[k]) + except(KeyError): pass + try: del(self._Cempirical[k]) + except(KeyError): pass + try: del(self._I[k]) + except(KeyError): pass + try: del(self._iC[k]) + except(KeyError): pass + + def C(self, k, t=None): + """ + TODO. + + Parameters + ---------- + k : TODO + TODO. + t : TODO, optional + TODO. Default: None. + """ + if not self._C.has_key(k): # defaults to true covariance matrix + self.set_C({k:cov(self.x[k], self.w[k])}) + self._Cempirical[k] = self._C[k] # save computing this later, if we are making it now + if t is None: return self._C[k] + # If t is provided, Calculate C for the provided time index, including flagging + w = self.w[k][:,t:t+1] + return self._C[k] * (w * w.T) + + def set_C(self, d): + """ + Set covariance to user-provided values. + + Parameters + ---------- + d : dict + Dictionary containing new covariance values for each key. + """ + self.clear_cache(d.keys()) + for k in d: self._C[k] = d[k] + + def C_empirical(self, k): + """ + Calculate empirical covariance from the data (with appropriate + weighting). + + Parameters + ---------- + k : TODO + Key. TODO. + + Returns + ------- + C_empirical : TODO + Empirical covariance for the specified key. + """ + if not self._Cempirical.has_key(k): + # Must be the actual empirical covariance; no overwriting + self._Cempirical[k] = cov(self.x[k], self.w[k]) + return self._Cempirical[k] + + def I(self, k): + """ + TODO. + + Parameters + ---------- + k : TODO + TODO. + + Returns + ------- + I : TODO + TODO. + """ + if not self._I.has_key(k): + nchan = self.x[k].shape[0] + self.set_I({k:np.identity(nchan)}) + return self._I[k] + + def set_I(self, d): + """ + Set the values along the diagonal of the 'identity' covariance matrix + to user-specified values. + + Parameters + ---------- + d : dict + Dictionary of values that should be inserted along the diagonal of + the identity covariance, I. + """ + for k in d: self._I[k] = d[k] + + def iC(self, k, t=None, rcond=1e-12): + """ + TODO. + + Parameters + ---------- + k : TODO + TODO. + t : TODO, optional + TODO. + rcond : TODO + TODO. Default: 1e-12 + + Returns + ------- + iCt : TODO + TODO. + """ + if not self._iC.has_key(k): + C = self.C(k) + U,S,V = np.linalg.svd(C.conj()) # conj in advance of next step + if self.lmin != None: S += self.lmin # ensure invertibility + if self.lmode != None: S += S[self.lmode-1] + self.set_iC({k:np.einsum('ij,j,jk', V.T, 1./S, U.T)}) + + if t is None: + return self._iC[k] + + # If t is provided, calculate iC for the provided time index, including flagging + # XXX this does not respect manual setting of iC with ds.set_iC + UserWarning("This does not respect manual setting of iC with ds.set_iC") + w = self.w[k][:,t:t+1] + m = hash(w) + if not self._iCt.has_key(k): self._iCt[k] = {} + if not self._iCt[k].has_key(m): + self._iCt[k][m] = np.linalg.pinv(self.C(k,t), rcond=rcond) + return self._iCt[k][m] + + def set_iC(self, d): + """ + TODO. + + Parameters + ---------- + d : TODO + TODO. + """ + for k in d: self._iC[k] = d[k] + + def to_npz(self, filename): + """ + This function will be removed. + """ + DeprecationWarning("This function will be removed.") + data = {} + for k in self.x: + sk = str(k) + data['x '+sk] = self.x[k] + try: data['C '+sk] = self.C[k] + except(AttributeError): pass + try: data['iC '+sk] = self.iC[k] + except(AttributeError): pass + for k in self.lsts: + data['lst '+k] = self.lsts[k] + np.savez(filename, **data) + + def from_npz(self, filename): + """ + This function will be removed. + """ + DeprecationWarning("This function will be removed.") + npz = np.load(filename) + self.x = {} + for k in [f for f in npz.files if f.startswith('x')]: + self.x[eval(k[2:])] = npz[k] + self.C = {} + for k in [f for f in npz.files if f.startswith('C')]: + self.C[eval(k[2:])] = npz[k] + self.iC = {} + for k in [f for f in npz.files if f.startswith('iC')]: + self.iC[eval(k[3:])] = npz[k] + self.lsts = {} + for k in [f for f in npz.files if f.startswith('lst')]: + self.lsts[k[4:]] = npz[k] + dsets, bls, pols = {}, {}, {} + for k in self.x: + dsets[k[0]] = None + bls[k[1]] = None + pols[k[2]] = None + self.dsets = dsets.keys() + self.bls = bls.keys() + self.pols = pols.keys() + + def gen_bl_boots(self, nboots, ngroups=5, ngps=None): + """ + TODO. This function will be moved elsewhere. + + Parameters + ---------- + nboots : int + Number of bootstraps to perform. + ngroups : int, optional + TODO. Default: 5. + """ + if ngps is not None: + DeprecationWarning("The 'ngps' kwarg has been renamed to 'ngroups'.") + ngroups = ngps + DeprecationWarning("This function will be moved elsewhere.") + _bls = {} + for k in self.x: _bls[k[1]] = None + for boot in xrange(nboots): + bls = _bls.keys()[:] + random.shuffle(bls) + gps = [bls[i::ngps] for i in range(ngps)] + gps = [[random.choice(gp) for bl in gp] for gp in gps] + yield gps + return + + def gen_gps(self, bls, ngroups=5, ngps=None): + """ + TODO. This function will be moved elsewhere. + + Parameters + ---------- + bls : TODO + TODO. + ngroups : int, optional + TODO. Default: 5. + + Returns + ------- + gps : TODO + TODO. + """ + # Remove this check when the code freezes + if ngps is not None: + DeprecationWarning("The 'ngps' kwarg has been renamed to 'ngroups'.") + ngroups = ngps + DeprecationWarning("This function will be moved elsewhere.") + + random.shuffle(bls) + gps = [bls[i::ngps] for i in range(ngps)] + #gps = [[random.choice(gp) for bl in gp] for gp in gps] # Sample w/replacement inside each group + # All independent baselines within a group, except the last one which + # is sampled randomly + gps = [ + [ gps[gp][i] for i in np.append(np.arange(0,len(gps[gp])-1), + random.choice(np.arange(0,len(gps[gp]))))] + for gp, G in enumerate(gps) ] + return gps + + def group_data(self, keys, gps, use_cov=True): + """ + TODO. + + Parameters + ---------- + keys : tuple + TODO. Keys have format (k, bl, POL). + + gps : TODO + TODO. + use_cov : bool, optional + If True, include the original empirical covariance, summed into the + specified groups, in the returned DataSet object. If False, set + the covariance to the identity. Default: True. + + Returns + ------- + newkeys : TODO + TODO. + ds : DataSet + New DataSet object containing TODO. + """ + # XXX avoid duplicate code for use_cov=True vs False + # (i.e. no separate dsC & dsI) + sets = np.unique([key[0] for key in keys]) + POL = keys[0][2] + nchan = self.x[keys[0]].shape[0] + + newkeys = [] + dsC_data, dsI_data = {}, {} + iCsum, iCxsum, Ixsum, Isum = {}, {}, {}, {} + + for s in sets: + # Sum up data for each group and make new keys + for gp in range(len(gps)): + newkey = (s, gp) + newkeys.append(newkey) + + # XXX: Replace with np.sum? + iCsum[newkey] = sum([self.iC((s,bl,POL)) for bl in gps[gp]]) + iCxsum[newkey] = sum([np.dot(self.iC((s,bl,POL)), + self.x[(s,bl,POL)]) for bl in gps[gp]]) + Isum[newkey] = sum([np.identity(nchan) for bl in gps[gp]]) + Ixsum[newkey] = sum([self.x[(s,bl,POL)] for bl in gps[gp]]) + + # Find effective summed up x based on iCsum and iCxsum + dsC_data[newkey] = np.dot(np.linalg.inv(iCsum[newkey]), + iCxsum[newkey]).T + # Find effective summed up x based on Isum and Ixsum + dsI_data[newkey] = np.dot(np.linalg.inv(Isum[newkey]), + Ixsum[newkey]).T + + # Create new DataSet objects to store result of grouping operation + # (I has to be separate because it has different x's populated into it) + dsC = DataSet(); dsC.add_data(dsets=dsC_data, replace=True) + dsI = DataSet(); dsI.add_data(dsets=dsI_data, replace=True) + + # Override default empirical covariance, since they're incorrect if + # computed from x + dsC.set_iC(iCsum) + dsI.set_I(Isum) + + if use_cov: + return newkeys, dsC + else: + return newkeys, dsI + + def q_hat(self, k1, k2, use_cov=True, use_fft=True, cov_flagging=False): + """ + TODO. + + Parameters + ---------- + k1, k2 : TODO + TODO. + use_cov : bool, optional + TODO. Default: True. + use_fft : bool, optional + TODO. Default: True. + cov_flagging : bool, optional + TODO. Default: False. + """ + try: + if self.x.has_key(k1): k1 = [k1] + except(TypeError): pass + try: + if self.x.has_key(k2): k2 = [k2] + except(TypeError): pass + nchan = self.x[k1[0]].shape[0] + if use_cov: + if not cov_flagging: + #iC1,iC2 = self.iC(k1), self.iC(k2) + #iC1x, iC2x = np.dot(iC1, self.x[k1]), np.dot(iC2, self.x[k2]) + iC1x, iC2x = 0, 0 + for k1i in k1: iC1x += np.dot(self.iC(k1i), self.x[k1i]) + for k2i in k2: iC2x += np.dot(self.iC(k2i), self.x[k2i]) + else: + raise NotImplementedError("cov_flagging not implemented yet.") + """ + # XXX make this work with k1,k2 being lists + iCx = {} + for k in (k1,k2): + iCx[k] = np.empty_like(self.x[k]) + inds = {} + w = self.w[k] + ms = [hash(w[:,i]) for i in xrange(w.shape[1])] + for i,m in enumerate(ms): inds[m] = inds.get(m,[]) + [i] + iCxs = {} + for m in inds: + x = self.x[k][:,inds[m]] + #x = self.x[k].take(inds[m], axis=1) + iC = self.iC(k,inds[m][0]) + iCx[k][:,inds[m]] = np.dot(iC,x) + #iCx[k].put(inds[m], np.dot(iC,x), axis=1) + iC1x,iC2x = iCx[k1], iCx[k2] + """ + else: + # XXX make this work with k1,k2 being lists + #iC1x, iC2x = self.x[k1].copy(), self.x[k2].copy() + iC1x, iC2x = 0, 0 + for k1i in k1: iC1x += np.dot(self.I(k1i), self.x[k1i]) + for k2i in k2: iC2x += np.dot(self.I(k2i), self.x[k2i]) + #iC1x, iC2x = np.dot(self.I(k1), self.x[k1]), np.dot(self.I(k2), self.x[k2]) + # Whether to use FFT or slow direct method + if use_fft: + _iC1x = np.fft.fft(iC1x.conj(), axis=0) + _iC2x = np.fft.fft(iC2x.conj(), axis=0) + + # Conjugated because inconsistent with pspec_cov_v003 otherwise + return np.conj( np.fft.fftshift(_iC1x,axes=0).conj() + * np.fft.fftshift(_iC2x,axes=0) ) + else: + # Slow method, used to explicitly cross-check fft code + # XXX make this work with k1,k2 being lists + q = [] + for i in xrange(nchan): + Q = self.get_Q(i, nchan) + iCQiC = np.einsum('ab,bc,cd', iC1.T.conj(), Q, iC2) # C^-1 Q C^-1 + qi = np.sum(self.x[k1].conj() * np.dot(iCQiC,self.x[k2]), axis=0) + q.append(qi) + return np.array(q) + + def get_F(self, k1, k2, use_cov=True, cov_flagging=False): + """ + TODO. + + Parameters + ---------- + k1, k2 : + TODO. + use_cov : bool, optional + TODO. Default: True. + cov_flagging : bool, optional + TODO. Default: False. + + Returns + ------- + F : TODO + Fisher matrix. + """ + try: + if self.x.has_key(k1): k1 = [k1] + except(TypeError): pass + try: + if self.x.has_key(k2): k2 = [k2] + except(TypeError): pass + + nchan = self.x[k1[0]].shape[0] + if use_cov: + + if not cov_flagging: + F = np.zeros((nchan,nchan), dtype=np.complex) + #iC1,iC2 = self.iC(k1), self.iC(k2) + iC1, iC2 = 0, 0 + for k1i in k1: iC1 += self.iC(k1i) + for k2i in k2: iC2 += self.iC(k2i) + #if False: Cemp1, Cemp2 = self.C_empirical(k1), self.C_empirical(k2) + else: + raise NotImplementedError("cov_flagging does not currently work.") + """ + # XXX make this work with k1,k2 being lists + # This is for the "effective" matrix s.t. W=MF and p=Mq + F = {} + w1,w2 = self.w[k1], self.w[k2] + m1s = [hash(w1[:,i]) for i in xrange(w1.shape[1])] + m2s = [hash(w2[:,i]) for i in xrange(w2.shape[1])] + for m1,m2 in zip(m1s,m2s): F[(k1,m1,k2,m2)] = None + for k1,m1,k2,m2 in F.keys(): + #for m1 in self._iCt[k1]: # XXX not all m1/m2 pairs may exist in data + # for m2 in self._iCt[k2]: + F[(k1,m1,k2,m2)] = np.zeros((nchan,nchan), dtype=np.complex) + iCQ1,iCQ2 = {}, {} + for ch in xrange(nchan): # this loop is nchan^3 + Q = self.get_Q(ch,nchan) + iCQ1[ch] = np.dot(self._iCt[k1][m1],Q) #C^-1 Q # If ERROR: Compute q_hat first + iCQ2[ch] = np.dot(self._iCt[k2][m2],Q) #C^-1 Q + for i in xrange(nchan): # this loop goes as nchan^4 + for j in xrange(nchan): + F[(k1,m1,k2,m2)][i,j] += np.einsum('ij,ji', iCQ1[i], iCQ2[j]) #C^-1 Q C^-1 Q + return F + """ + else: + # Use identity matrix instead of default covariance + # XXX make this work with k1,k2 being lists + #iC1 = np.linalg.inv(self.C(k1) * np.identity(nchan)) + #iC2 = np.linalg.inv(self.C(k2) * np.identity(nchan)) + iC1, iC2 = 0, 0 + for k1i in k1: iC1 += self.I(k1i) + for k2i in k2: iC2 += self.I(k2i) + + # Do this to get the effective F (see below) + if False: Cemp1, Cemp2 = self.I(k1), self.I(k2) + F = np.zeros((nchan,nchan), dtype=np.complex) + #Cemp1, Cemp2 = self.C_empirical(k1), self.C_empirical(k2) + + if False: + # This is for the "true" Fisher matrix + CE1, CE2 = {}, {} + + for ch in xrange(nchan): + Q = self.get_Q(ch,nchan) + # C1 Cbar1^-1 Q Cbar2^-1, C2 Cbar2^-1 Q Cbar1^-1 + CE1[ch] = np.dot(Cemp1, np.dot(iC1, np.dot(Q, iC2))) + CE2[ch] = np.dot(Cemp2, np.dot(iC2, np.dot(Q, iC1))) + + for i in xrange(nchan): + for j in xrange(nchan): + F[i,j] += np.einsum('ij,ji', CE1[i], CE2[j]) # C E C E + else: + # This is for the "effective" matrix s.t. W=MF and p=Mq + iCQ1, iCQ2 = {}, {} + + for ch in xrange(nchan): # this loop is nchan^3 + Q = self.get_Q(ch, nchan) + iCQ1[ch] = np.dot(iC1, Q) #C^-1 Q + iCQ2[ch] = np.dot(iC2, Q) #C^-1 Q + + for i in xrange(nchan): # this loop goes as nchan^4 + for j in xrange(nchan): + F[i,j] += np.einsum('ij,ji', iCQ1[i], iCQ2[j]) #C^-1 Q C^-1 Q + return F + + def get_MW(self, F, mode='F^-1'): + """ + TODO + + Parameters + ---------- + F : TODO + TODO. + mode : str, optional + TODO. Default: 'F^-1' + + Returns + ------- + M : TODO + TODO. + + W : TODO + Weight matrix. + """ + if type(F) is dict: # recursive case for many F's at once + M,W = {}, {} + for key in F: M[key],W[key] = self.get_MW(F[key], mode=mode) + return M,W + modes = ['F^-1', 'F^-1/2', 'I', 'L^-1']; assert(mode in modes) + if mode == 'F^-1': + M = np.linalg.pinv(F, rcond=1e-12) + #U,S,V = np.linalg.svd(F) + #M = np.einsum('ij,j,jk', V.T, 1./S, U.T) + elif mode == 'F^-1/2': + U,S,V = np.linalg.svd(F) + M = np.einsum('ij,j,jk', V.T, 1./np.sqrt(S), U.T) + elif mode == 'I': + M = np.identity(F.shape[0], dtype=F.dtype) + else: + # Cholesky decomposition to get M (XXX: Needs generalizing) + order = np.array([10, 11, 9, 12, 8, 20, 0, + 13, 7, 14, 6, 15, 5, 16, + 4, 17, 3, 18, 2, 19, 1]) + iorder = np.argsort(order) + F_o = np.take(np.take(F,order, axis=0), order, axis=1) + L_o = np.linalg.cholesky(F_o) + U,S,V = np.linalg.svd(L_o.conj()) + M_o = np.dot(np.transpose(V), np.dot(np.diag(1./S), np.transpose(U))) + M = np.take(np.take(M_o, iorder, axis=0), iorder, axis=1) + W = np.dot(M, F) + norm = W.sum(axis=-1); norm.shape += (1,) + M /= norm; W = np.dot(M, F) + return M,W + + def get_Q(self, mode, n_k, window='none', delay=False): + """ + Response of the covariance to a given bandpower, dC / dp_alpha. This + currently assumes that Q will operate on a visibility vector in frequency + space. + + Parameters + ---------- + mode : float + Central wavenumber of the bandpower. + n_k : int + Number of k bins that will be . + window : str, optional + TODO. Default: 'none'. + delay : bool, optional + TODO. Default: False. + + Returns + ------- + Q : array_like + Matrix operator that encodes the response of the data covariance C to a + change in bandpower p_alpha. + """ + if delay: + # XXX need to have this depend on window + if window is not 'none': + raise NotImplementedError("Window function not yet supported in delay mode.") + Q = np.zeros_like(C) + Q[mode,mode] = 1 + else: + _m = np.zeros((n_k,), dtype=np.complex) + _m[mode] = 1. # delta function at specified delay mode + + # FFT to convert to frequency domain + m = np.fft.fft(np.fft.ifftshift(_m)) * aipy.dsp.gen_window(n_k, window) + Q = np.einsum('i,j', m, m.conj()) # dot it with its conjugate + return Q + + + def p_hat(self, M, q, scalar=1.): + """ + TODO + + Parameters + ---------- + M : TODO + TODO. + q : TODO + TODO. + scalar : int, optional + TODO. Default: 1. + + Returns + ------- + TODO + """ + if type(M) is dict: # we have different M's for different times + (k1,m1,k2,m2) = M.keys()[0] + w1,w2 = self.w[k1], self.w[k2] + m1s = [hash(w1[:,i]) for i in xrange(w1.shape[1])] + m2s = [hash(w2[:,i]) for i in xrange(w2.shape[1])] + inds = {} + for i,(m1,m2) in enumerate(zip(m1s,m2s)): + inds[(k1,m1,k2,m2)] = inds.get((k1,m1,k2,m2),[]) + [i] + p = np.zeros_like(q) + for key in inds: + qi = q[:,inds[key]] + #qi = q.take(inds[key], axis=1) + p[:,inds[key]] = np.dot(M[key], qi) * scalar + #p.put(inds[key], np.dot(M[key], qi) * scalar, axis=1) + return p + else: + return np.dot(M, q) * scalar + + def pspec(self, keys, weights): + """ + Calculate power spectrum for this DataSet using the optimal quadratic + estimator (OQE) from (CITE PAPER). + + Parameters + ---------- + keys : dict + xxx + weights : TODO + TODO. + + Returns + ------- + pspec : list + Optimal quadratic estimate of the power spectrum for the baselines + specified in 'keys'. + """ + pvs = [] + for k, key1 in enumerate(keys): + if k == 1 and len(keys) == 2: + # NGPS = 1 (skip 'odd' with 'even' if we already did 'even' + # with 'odd') + continue + + for key2 in keys[k:]: + if len(keys) > 2 and (key1[0] == key2[0] or key1[1] == key2[1]): + # NGPS > 1 + continue + if key1[0] == key2[0]: # don't do 'even' with 'even', for example + continue + else: + Fv = self.get_F(key1, key2) + qv = self.q_hat(key1, key2) + Mv, Wv = self.get_MW(Fv, mode=weights) + pv = self.p_hat(Mv, qv, scalar=scalar) + pvs.append(pv) + return pvs + diff --git a/hera_pspec/pspec.py b/hera_pspec/pspec.py index 4ff2830a..a9b1d221 100644 --- a/hera_pspec/pspec.py +++ b/hera_pspec/pspec.py @@ -1,4 +1,6 @@ '''Units in mK, GHz, Mpc unless stated otherwise''' +#XXX: This code is provided for legacy purposes only. + import numpy as n, aipy as a #from binning import LST_RES, UV_RES, uv2bin, bin2uv, rebin_log diff --git a/hera_pspec/tests/test_conversions.py b/hera_pspec/tests/test_conversions.py new file mode 100644 index 00000000..e7b67359 --- /dev/null +++ b/hera_pspec/tests/test_conversions.py @@ -0,0 +1,56 @@ +import unittest +import nose.tools as nt +import numpy as np +import os +import sys +from hera_pspec.data import DATA_PATH +from hera_pspec import conversions + + +class Test_Cosmo(unittest.TestCase): + + def setUp(self): + self.C = conversions.Cosmo_Conversions(Om_L=0.68440, Om_b=0.04911, Om_c=0.26442, H0=100.0, + Om_M=None, Om_k=None) + + def tearDown(self): + pass + + def runTest(self): + pass + + def test_init(self): + # test empty instance + C = conversions.Cosmo_Conversions() + # test parameters exist + C.Om_b + C.H0 + # test parameters get fed to class + C = conversions.Cosmo_Conversions(H0=25.5) + self.assertAlmostEqual(C.H0, 25.5) + + def test_distances(self): + self.assertAlmostEqual(self.C.f2z(100e6), 13.20405751) + self.assertAlmostEqual(self.C.z2f(10.0), 129127795.54545455) + self.assertAlmostEqual(self.C.E(10.0), 20.450997530682947) + self.assertAlmostEqual(self.C.DC(10.0), 6499.708111027144) + self.assertAlmostEqual(self.C.DM(10.0), 6510.2536925709637) + self.assertAlmostEqual(self.C.DA(10.0), 591.84124477917851) + self.assertAlmostEqual(self.C.dRperp_dtheta(10.0), 6510.2536925709637) + self.assertAlmostEqual(self.C.dRpara_df(10.0), 1.2487605057418872e-05) + self.assertAlmostEqual(self.C.dRpara_df(10.0, ghz=True), 12487.605057418872) + self.assertAlmostEqual(self.C.X2Y(10.0), 529.26719942209002) + + + + + + + + + + + + +if __name__ == "__main__": + unittest.main() diff --git a/hera_pspec/tests/test_version.py b/hera_pspec/tests/test_version.py new file mode 100644 index 00000000..623f1473 --- /dev/null +++ b/hera_pspec/tests/test_version.py @@ -0,0 +1,46 @@ +"""Tests for version.py.""" +import nose.tools as nt +import sys +import os +from StringIO import StringIO +import subprocess +import hera_pspec + +def test_main(): + version_info = hera_pspec.version.construct_version_info() + + saved_stdout = sys.stdout + try: + out = StringIO() + sys.stdout = out + hera_pspec.version.main() + output = out.getvalue() + nt.assert_equal(output, 'Version = {v}\ngit origin = {o}\n' + 'git branch = {b}\ngit description = {d}\n' + .format(v=version_info['version'], + o=version_info['git_origin'], + b=version_info['git_branch'], + d=version_info['git_description'])) + finally: + sys.stdout = saved_stdout + + +def test_main(): + version_info = hera_pspec.version.construct_version_info() + + saved_stdout = sys.stdout + try: + out = StringIO() + sys.stdout = out + hera_pspec.version.main() + output = out.getvalue() + nt.assert_equal(output, 'Version = {v}\ngit origin = {o}\n' + 'git branch = {b}\ngit description = {d}\n' + .format(v=version_info['version'], + o=version_info['git_origin'], + b=version_info['git_branch'], + d=version_info['git_description'])) + finally: + sys.stdout = saved_stdout + + diff --git a/hera_pspec/utils.py b/hera_pspec/utils.py new file mode 100644 index 00000000..838059b7 --- /dev/null +++ b/hera_pspec/utils.py @@ -0,0 +1,133 @@ + +import numpy as np +import aipy +import md5 + +DELAY = False + +def hash(w): + """ + Return an MD5 hash of a set of weights. + """ + return md5.md5(w.copy(order='C')).digest() + + +def cov(d1, w1, d2=None, w2=None): + """ + Calculate empirical covariance of data vectors with weights. + + Parameters + ---------- + d1, w1 : array_like + Data vector and weight vector, which should be the same shape. If only + these are specified, the covariance of d1 with itself (conjugated) will + be calculated. + + d2, w2 : array_like, optional + Second data vector and weight vector. If specified, the covariance of + d1 with d2 will be calculated. + + Returns + ------- + cov : array_like + Empirical covariance matrix of the weighted data vector(s). + """ + if d2 is None: d2, w2 = d1.conj(), w1 + + # Weighted data vectors + d1sum, d1wgt = (w1*d1).sum(axis=1), w1.sum(axis=1) + d2sum, d2wgt = (w2*d2).sum(axis=1), w2.sum(axis=1) + + # Weighted means + x1 = d1sum / np.where(d1wgt > 0, d1wgt, 1) + x2 = d2sum / np.where(d2wgt > 0, d2wgt, 1) + x1.shape = (-1,1); x2.shape = (-1,1) + d1x = d1 - x1 + d2x = d2 - x2 + + # Calculate covariance and weight matrices + C = np.dot(w1*d1x, (w2*d2x).T) + W = np.dot(w1, w2.T) + return C / np.where(W > 1, W-1, 1) + +""" +def noise(size): + sig = 1./np.sqrt(2) + return np.random.normal(scale=sig, size=size) + 1j*np.random.normal(scale=sig, size=size) + +def lst_grid(lsts, data, wgts=None, lstbins=6300, wgtfunc=lambda dt,res: np.exp(-dt**2/(2*res**2))): + lstgrid = np.linspace(0, 2*np.pi, lstbins) + lstres = lstgrid[1]-lstgrid[0] + if wgts is None: wgts = np.where(np.abs(data) == 0, 0, 1.) + sumgrid,wgtgrid = 0, 0 + for lst,d,w in zip(lsts,data,wgts): + dt = lstgrid - lst + wf = wgtfunc(dt,lstres); wf.shape = (-1,) + (1,)*(data.ndim-1) + d.shape, w.shape = (1,-1), (1,-1) + wgtgrid += w * wf + sumgrid += d * w * wf + datgrid = np.where(wgtgrid > 1e-10, sumgrid/wgtgrid, 0) + return lstgrid, datgrid, wgtgrid + +def lst_grid_cheap(lsts, data, wgts=None, lstbins=6300, wgtfunc=lambda dt,res: np.exp(-dt**2/(2*res**2))): + lstgrid = np.linspace(0, 2*np.pi, lstbins) + lstres = lstgrid[1]-lstgrid[0] + if wgts is None: wgts = np.where(np.abs(data) == 0, 0, 1.) + sumgrid = np.zeros((lstbins,)+data.shape[1:], dtype=data.dtype) + wgtgrid = np.zeros(sumgrid.shape, dtype=wgts.dtype) + for lst,d,w in zip(lsts,data,wgts): + i,j = int(np.floor(lst/lstres)), int(np.ceil(lst/lstres)) + wi,wj = wgtfunc(lst-lstgrid[i],lstres), wgtfunc(lst-lstgrid[j],lstres) + sumgrid[i] += d * w * wi; wgtgrid[i] += w * wi + sumgrid[j] += d * w * wj; wgtgrid[j] += w * wj + datgrid = np.where(wgtgrid > 1e-10, sumgrid/wgtgrid, 0) + return lstgrid, datgrid, wgtgrid + +def lst_align(lsts, lstres=.001, interpolation='none'): + lstgrid = np.arange(0, 2*np.pi, lstres) + lstr, order = {}, {} + for k in lsts: #orders LSTs to find overlap + order[k] = np.argsort(lsts[k]) + lstr[k] = np.around(lsts[k][order[k]] / lstres) * lstres + lsts_final = None + for i,k1 in enumerate(lstr.keys()): + for k2 in lstr.keys()[i:]: + if lsts_final is None: lsts_final = np.intersect1d(lstr[k1],lstr[k2]) #XXX LSTs much match exactly + else: lsts_final = np.intersect1d(lsts_final,lstr[k2]) + inds = {} + for k in lstr: #selects correct LSTs from data + inds[k] = order[k].take(lstr[k].searchsorted(lsts_final)) + return inds + +def lst_align_data(inds, dsets, wgts=None, lsts=None): + for k in dsets: + k0 = k[0] + dsets[k] = dsets[k][inds[k0]] + if not wgts is None: wgts[k] = wgts[k][inds[k0]] + if not lsts is None: + for k0 in lsts: lsts[k0] = lsts[k0][inds[k0]] + return [d for d in [dsets,wgts,lsts] if not d is None] + +def boot_waterfall(p, axis=1, nsamples=None, nboots=1000, usemedian=True, verbose=False): + boots = [] + dim = p.shape[axis] + if nsamples is None: nsamples = dim + for i in xrange(nboots): + if verbose and i % 10 == 0: print ' ', i, '/', nboots + inds = np.random.randint(0,dim,size=nsamples) + if usemedian: pi = np.median(p.take(inds,axis=axis), axis=1) + else: pi = np.average(p.take(inds,axis=axis), axis=1) + boots.append(pi) + # XXX deal with folding + boots = np.array(boots) + if verbose: print 'Sorting bootstraps...' + pk = np.average(boots, axis=0) #average over all boots + #this is excluding imag component in noise estimate ` + boots = np.sort(boots.real, axis=0) #dropping imag component here + up_thresh = int(np.around(0.975 * boots.shape[0])) #2 sigma, single tail + dn_thresh = int(np.around(0.025 * boots.shape[0])) #2 sigma, single tail + #important to only include real component in estimation of error + err_up = (boots[up_thresh] - pk.real) / 2 #effective "1 sigma" derived from actual 2 sigma + err_dn = -(boots[dn_thresh] - pk.real) / 2 #effective "1 sigma" derived from actual 2 sigma + return pk, (err_up,err_dn), boots +""" diff --git a/hera_pspec/version.py b/hera_pspec/version.py new file mode 100644 index 00000000..bda01c40 --- /dev/null +++ b/hera_pspec/version.py @@ -0,0 +1,72 @@ +import os +import subprocess +import json + +hera_pspec_dir = os.path.dirname(os.path.realpath(__file__)) +version_file = os.path.join(hera_pspec_dir, 'VERSION') +version = open(version_file).read().strip() + + +def construct_version_info(): + hera_pspec_dir = os.path.dirname(os.path.realpath(__file__)) + version_file = os.path.join(hera_pspec_dir, 'VERSION') + version = open(version_file).read().strip() + + try: + git_origin = subprocess.check_output(['git', '-C', hera_pspec_dir, 'config', + '--get', 'remote.origin.url'], + stderr=subprocess.STDOUT).strip() + git_hash = subprocess.check_output(['git', '-C', hera_pspec_dir, 'rev-parse', 'HEAD'], + stderr=subprocess.STDOUT).strip() + git_description = subprocess.check_output(['git', '-C', hera_pspec_dir, + 'describe', '--dirty', '--tag', '--always']).strip() + git_branch = subprocess.check_output(['git', '-C', hera_pspec_dir, 'rev-parse', + '--abbrev-ref', 'HEAD'], + stderr=subprocess.STDOUT).strip() + git_version = subprocess.check_output(['git', '-C', hera_pspec_dir, 'describe', '--always', + '--tags', '--abbrev=0']).strip() + except: # pragma: no cover - can't figure out how to test exception. + try: + # Check if a GIT_INFO file was created when installing package + git_file = os.path.join(hera_pspec_dir, 'GIT_INFO') + with open(git_file) as data_file: + data = [x.encode('UTF8') for x in json.loads(data_file.read().strip())] + git_origin = data[0] + git_hash = data[1] + git_description = data[2] + git_branch = data[3] + except: + git_origin = '' + git_hash = '' + git_description = '' + git_branch = '' + + version_info = {'version': version, 'git_origin': git_origin, + 'git_hash': git_hash, 'git_description': git_description, + 'git_branch': git_branch} + return version_info + +version_info = construct_version_info() +version = version_info['version'] +git_origin = version_info['git_origin'] +git_hash = version_info['git_hash'] +git_description = version_info['git_description'] +git_branch = version_info['git_branch'] + +# String to add to history of any files written with this version of pyuvdata +hera_pspec_version_str = ('hera_pspec version: ' + version + '.') +if git_hash is not '': + hera_pspec_version_str += (' Git origin: ' + git_origin + + '. Git hash: ' + git_hash + + '. Git branch: ' + git_branch + + '. Git description: ' + git_description + '.') + + +def main(): + print('Version = {0}'.format(version)) + print('git origin = {0}'.format(git_origin)) + print('git branch = {0}'.format(git_branch)) + print('git description = {0}'.format(git_description)) + +if __name__ == '__main__': + main() diff --git a/setup.py b/setup.py index 0658cc85..3c910e9e 100644 --- a/setup.py +++ b/setup.py @@ -1,22 +1,38 @@ from setuptools import setup -import glob -import os import sys +import os +from hera_pspec import version import json +data = [version.git_origin, version.git_hash, version.git_description, version.git_branch] +with open(os.path.join('hera_pspec', 'GIT_INFO'), 'w') as outfile: + json.dump(data, outfile) + +def package_files(package_dir, subdirectory): + # walk the input package_dir/subdirectory + # return a package_data list + paths = [] + directory = os.path.join(package_dir, subdirectory) + for (path, directories, filenames) in os.walk(directory): + for filename in filenames: + path = path.replace(package_dir + '/', '') + paths.append(os.path.join(path, filename)) + return paths +data_files = package_files('hera_pspec', 'data') + setup_args = { - 'name': 'hera_pspec', - 'author': 'HERA Team', - 'url': 'https://github.com/HERA-Team/hera_pspec', - 'license': 'BSD', - 'description': 'collection of calibration routines to run on the HERA instrument.', - 'package_dir': {'hera_pspec': 'hera_pspec'}, - 'packages': ['hera_pspec'], + 'name': 'hera_pspec', + 'author': 'HERA Team', + 'url': 'https://github.com/HERA-Team/hera_pspec', + 'license': 'BSD', + 'version': version.version, + 'description': 'HERA Power Spectrum Estimator Code.', + 'packages': ['hera_pspec'], + 'package_dir': {'hera_pspec': 'hera_pspec'}, + 'package_data': {'hera_pspec': data_files}, 'include_package_data': True, - 'scripts': [], - 'zip_safe': False, + 'zip_safe': False, } - if __name__ == '__main__': apply(setup, (), setup_args) diff --git a/transiting_code/.gitignore b/transiting_code/.gitignore new file mode 100644 index 00000000..0d20b648 --- /dev/null +++ b/transiting_code/.gitignore @@ -0,0 +1 @@ +*.pyc diff --git a/transiting_code/file_format/parameter.py b/transiting_code/file_format/parameter.py new file mode 100644 index 00000000..c2ef409b --- /dev/null +++ b/transiting_code/file_format/parameter.py @@ -0,0 +1,25 @@ +import numpy as np + +class psparam(object): + """ + psp == "power spectrum parameter" + """ + + def __init__(self, name, required=True, value=None, form=(), units=None, + description='', expected_type=None, tols=(1e-5, 1e-08)): + """Init power spectrum parameter object.""" + self.name = name + self.required = required + self.value = value + self.description = description + self.units = units + self.form = form + if self.form == 'str': + self.expected_type = str + else: + self.expected_type = expected_type + if np.size(tols) == 1: + # Only one tolerance given, assume absolute, set relative to zero + self.tols = (0, tols) + else: + self.tols = tols # relative and absolute tolerances to be used in np.isclose \ No newline at end of file diff --git a/transiting_code/file_format/pspec.py b/transiting_code/file_format/pspec.py new file mode 100644 index 00000000..184f4e27 --- /dev/null +++ b/transiting_code/file_format/pspec.py @@ -0,0 +1,327 @@ +import numpy as np +import parameter as prm +from pspecbase import pspecBase + +class pspec(pspecBase): + """ + A class for storing power spectrum data. + """ + + def __init__(self, Ndims, length_units=None, temp_units=None,\ + freq_units=None): + """ + Container class for power spectrum data. + + Parameters + ---------- + Ndims : int + Number of dimensions of the power spectrum. + + length_units, temp_units, freq_units : str + String specifying the units that have been assumed for lengths + (e.g. Mpc or Mpc/h), temperatures (e.g. K or mK), and frequencies + (e.g. MHz or GHz). These arguments are compulsory. + """ + # Test whether units have been specified + if length_units is None: + raise KeyError('Must specify length_units e.g. Mpc/h') + if temp_units is None: + raise KeyError('Must specify temp_units e.g. mK') + if freq_units is None: + raise KeyError('Must specify freq_units e.g. GHz') + self.check_units(length_units,temp_units) + + if Ndims not in (1, 2, 3): + raise ValueError("Power spectrum cannot have %d dimensions." \ + % Ndims) + + # Units used + desc = ('Units used in the power spectrum object' + 'Acceptable length units: Mpc, Mpc/h, Gpc, Gpc/h' + 'Acceptable temperature units: K, mK, uK, muK, microK') + self._lengthunits = prm.psparam('lengthunits', description=desc,\ + expected_type=str) + self._tempunits = prm.psparam('tempunits', description=desc,\ + expected_type=str) + + # Get dimensionality of power spectrum + desc = ('Type of pspec: 1 = spherically averaged P(k); ' + '2 = cylindrically averaged P(kperp, kpara); ' + '3 = full 3D cube P(kx,ky,kz)') + self._Ndims = prm.psparam('Ndims', description=desc, expected_type=int) + + # Define binning of power spectrum in wavenumber, k + desc = ('Number of bins along each k direction') + self._Nks = prm.psparam('Nks', description=desc, expected_type=int, + form=('Ndims',)) + desc = ('Boundaries of k bins. List (length Ndims) of arrays, ' + 'each of which is of shape (Nks[i],2), where each row ' + 'of the array gives the bottom edge and top edge of a ' + 'k bin.') + self._kbounds = prm.psparam('kbounds', description=desc, + expected_type=np.float, form=('Ndims',), + units='(%s)^-1' % length_units) + + # Define array for storing the power spectrum data + desc = ('Power spectrum data. An array that gives the power spectrum ' + 'values in each kbin. 1D, 2D, or 3D array depending on Ndims ' + 'parameter.') + self._pspec = prm.psparam('pspec', description=desc, + expected_type=np.float, + units='(%s)^2 (%s)^3' + % (temp_units, length_units)) + + # Specify central redshift of power spectrum + desc = ('Central redshift of power spectrum') + self._z0 = prm.psparam('z0', description=desc, expected_type=np.float) + + # Specify frequency bins that were used to form the power spectrum + desc = ('Frequencies used in forming the power spectrum') + self._freqs = prm.psparam('freqs', description=desc, + expected_type=np.float, + units='%s' % freq_units) + + # Specify window function array + desc = ('Window functions. An array where each row corresponds to ' + 'the linear combination of all the other k bins of the ' + 'true power spectrum that is probed when a power spectrum ' + 'estimate is made at a particular k bin.') + self._window = prm.psparam('window', description=desc, + expected_type=np.float) + + # Initialize base class with the set of parameters specified above + super(pspec, self).__init__() + + # Set the number of dimensions of the power spectrum + # E.g. Ndim = 1 is P(k), while Ndim = 3 is P(\vec{k}) + self._Ndims = Ndims + + self._lengthunits = length_units + self._tempunits = temp_units + + + def set_pspec(self, pk_data,kbin_edges=None): + """ + Set the power spectrum data to a new set of values. If kbins are not + provided, it is assumed that the input power spectrum is replacing + an existing power spectrum and has the same binning. The input data are + assumed to be in the same units that were used to initialize this + object. + + Parameters + ---------- + pkdata : array_like + New power spectrum data, assumed to have the same dimensionality, + binning, and units as the current data that it will replace. + kbin_edges : array_like + Edges of the k-bins. List (length Ndims) of arrays. + """ + if pk_data.ndim != self._Ndims: + raise ValueError("Input power spectrum has %s dimensions, but\ + the power spectrum container has %s dimensions." \ + % (pk_data.ndim,self._Ndims)) + + if kbin_edges is not None: + self._kbounds = kbin_edges + self._Nks = np.zeros(self._Ndims, dtype=int) + for i in range(self._Ndims): + self._Nks[i] = kbin_edges[i].shape[0] + + + + + # elif pk_data.ndim != len(kbin_edges): + # raise ValueError("Input power spectrum has dimensions %s, so\ + # we expect %s separate arrays for the kbin edges, but instead\ + # got %s arrays" \ + # % (pk_data.ndim, pk_data.ndim, len(kbin_edges))) + + # for i in range(self._Ndims): + # if pk_data. + + self._pspec = pk_data + + def bin_centers(self): + """ + Return the central k values for the power spectrum bins. + + Returns + ------- + k_centers : array_like + Central k values for each bin, in units of (length_units)^-1. + """ + k_centers = [] + for i in range(self._Ndims): + k_centers.append(np.mean(self._kbounds[i], axis=1)) + return k_centers + + def DeltaSq(self): + """ + Return the 'dimensionless' power spectrum, related to the usual power + spectrum by Delta^2 = k^3 P(k) / (2\pi^2). The units are actually + (temp_units)^2. + + Returns + ------- + DeltaSq : array_like + 'Dimensionless' power spectrum, in (temp_units)^2. + """ + + pspec = np.atleast_3d(self._pspec) + if self._Ndims == 1: + pspec = np.rollaxis(pspec,0,-1) + k_mag = np.zeros_like(self._pspec) + k_mids = self.bin_centers() + for i in range(self._Nks[0]): + for j in range(self._Nks[1]): + for k in range(self._Nks[2]): + k_mag[i,j,k] = k_mids[i]**2 + k_mids[j]**2 + k_mids[k]**2 + k_mag[i,j,k] = np.sqrt(k_mag[i,j,k]) + + return np.squeeze(pspec * k_mag**3 / (2. * np.pi)) + + def check_units(self, length_units, temp_units): + """ + Checks a set of units to see if they are understood by our object + + + Parameters + ---------- + length_units, temp_units : str + Units that the power spectrum should be returned in. + Acceptable length units: 'Mpc', 'Mpc/h', 'Gpc', 'Gpc/h' + Acceptable temperature units: 'K', 'mK', 'uK', 'muK', 'microK' + + """ + + # Define acceptable units (i.e. that the code knows how to deal with) + acceptable_length = ['Mpc', 'Mpc/h', 'h^-1 Mpc', 'h^-1Mpc', + 'Gpc', 'Gpc/h'] + acceptable_temp = ['K', 'mK', 'uK', 'muK', 'microK'] + + # Sanity checks on units + if length_units.lower() not in [l.lower() for l in acceptable_length]: + raise ValueError("Length units '%s' not recognized." \ + % length_units) + if temp_units.lower() not in [l.lower() for l in acceptable_temp]: + raise ValueError("Temperature units '%s' not recognized." \ + % temp_units) + + + def in_units(self, length_units, temp_units): + """ + Return the power spectrum in a different set of units. + + Parameters + ---------- + length_units, temp_units : str + Units that the power spectrum should be returned in. + Acceptable length units: 'Mpc', 'Mpc/h', 'Gpc', 'Gpc/h' + Acceptable temperature units: 'K', 'mK', 'uK', 'muK', 'microK' + + Returns + ------- + pspec : array_like + Array containing the power spectrum converted to the new set of + units (the stored power spectrum is not modified). + """ + + pspec = self._pspec.copy() + # Convert old and new units to a standard set of units (Mpc, K) + # Variables length_units and temp_units refer to the new units + # while self._lengthunits and self._tempunits are the old units + if 'h' in self._lengthunits and 'h' not in length_units: + pspec /= (cosmo_units.Ho)**3 + elif 'h' not in self._lengthunits and 'h' in length_units: + pspec *= (cosmo_units.Ho)**3 + + micro_old = False + if ('m' in self._lengthunits) or ('u' in self._lengthunits): + micro_old = True + micro_new = False + if ('m' in length_units) or ('u' in length_units): + micro_new = True + + if micro_old and not micro_new: + pspec /= 1000000. + elif micro_new and not micro_old: + pspec *= 1000000. + + # QUESTION: MAYBE WE SHOULD RETURN K VALUES AS WELL? + + return pspec + + def rebin(self, kbounds, method=None): + """ + Rebin the power spectrum. + + Parameters + ---------- + kbounds : array_like + Boundaries of k bins. List (length Ndims) of arrays, each of which + is of shape (Nks[i],2), where each row of the array gives the + bottom edge and top edge of a k bin. + + method : str + Rebinning method to use. Options are: + 'simple': Combine bins using a simple average. + 'Ninv': Combine bins using inverse noise weights. + + Returns + ------- + pspec : array_like + Power spectrum array with the new k binning. + """ + raise NotImplementedError() + # rebins the pspec without changing dimensionality, produces a new object + + def cylindrical(self): + """ + Return the cylindrically-averaged (i.e. 2D) power spectrum. This method + will only work if the stored power spectrum is at least 2D. + + If the stored power spectrum is 3D, it will perform a uniformly- + weighted cylindrical average of the data. + + Returns + ------- + pspec : pspec object + New 2D power spectrum object in cylindrical (k_par, k_perp) bins. + """ + raise NotImplementedError() + # returns a P(kperp, kpara) if available or does a uniform summing. + + def spherical(self): + """ + Return the spherically-averaged (i.e. 1D) power spectrum. This method + will only work if the stored power spectrum is at least 1D. + + If the stored power spectrum is 2D or 3D, it will perform a uniformly- + weighted spherical average of the data. + + Returns + ------- + pspec : pspec object + New 1D power spectrum object in spherical (|k|) bins. + """ + raise NotImplementedError() + # same as cylindrical + + def reweight(self): + """ + Apply a new set of weights to the stored power spectrum. + + Parameters + ---------- + w : array_like + Matrix operator for the new weights. + + Returns + ------- + pspec : pspec object + New power spectrum object, with power spectrum P_new = w^T P_old + and covariance C_new = w^T C_old w. + """ + raise NotImplementedError() + # applying a new M matrix, returns a new object + diff --git a/transiting_code/file_format/pspecbase.py b/transiting_code/file_format/pspecbase.py new file mode 100644 index 00000000..06267a21 --- /dev/null +++ b/transiting_code/file_format/pspecbase.py @@ -0,0 +1,41 @@ +import numpy as np +import parameter as prm + +class pspecBase(object): + def __init__(self): + """Create properties from pspec attributes.""" + for p in self: + print p + this_param = getattr(self, p) + attr_name = this_param.name + setattr(self.__class__, attr_name, property(self.prop_fget(p),self.prop_fset(p))) + + + + + def __iter__(self): + """Iterator for all UVParameter attributes.""" + attribute_list = [a for a in dir(self) if not a.startswith('__') and + not callable(getattr(self, a))] + param_list = [] + for a in attribute_list: + attr = getattr(self, a) + if isinstance(attr, prm.psparam): + param_list.append(a) + for a in param_list: + yield a + + def prop_fget(self, param_name): + """Getter method for PspecParameter properties.""" + def fget(self): + this_param = getattr(self, param_name) + return this_param.value + return fget + + def prop_fset(self, param_name): + """Setter method for PspecParameter properties.""" + def fset(self, value): + this_param = getattr(self, param_name) + this_param.value = value + setattr(self, param_name, this_param) + return fset \ No newline at end of file diff --git a/transiting_code/file_format/sandbox.ipynb b/transiting_code/file_format/sandbox.ipynb new file mode 100644 index 00000000..1fac431d --- /dev/null +++ b/transiting_code/file_format/sandbox.ipynb @@ -0,0 +1,365 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "import numpy as np, parameter as prm\n", + "from pspecbase import pspecBase" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "class pspec(pspecBase):\n", + " \"\"\"\n", + " A class for storing power spectrum data\n", + " \"\"\"\n", + " \n", + " def __init__(self, Ndims, length_units=None, temp_units=None, freq_units=None):\n", + " \n", + " if length_units is None:\n", + " raise KeyError('Must specify length_units e.g. Mpc/h')\n", + " if temp_units is None:\n", + " raise KeyError('Must specify temp_units e.g. mK')\n", + " if freq_units is None:\n", + " raise KeyError('Must specify freq_units e.g. GHz')\n", + " \n", + " desc = ('Type of pspec: 1 = spherically averaged P(k); '\n", + " '2 = cylindrically averaged P(kperp, kpara); '\n", + " '3 = full 3D cube P(kx,ky,kz)')\n", + " self._Ndims = prm.psparam('Ndims', description=desc, expected_type=int)\n", + " \n", + " desc = ('Number of bins along each k direction')\n", + " self._Nks = prm.psparam('Nks', description=desc, expected_type=int, form=('Ndims',))\n", + " \n", + " desc = ('Boundaries of k bins. List (length Ndims) of arrays, '\n", + " 'each of which is of shape (Nks[i],2), where each row '\n", + " 'of the array gives the bottom edge and top edge of a '\n", + " 'k bin.')\n", + " self._kbounds = prm.psparam('kbounds', description=desc, expected_type=np.float, form=('Ndims',),\\\n", + " units='(%s)^-1' % length_units)\n", + " \n", + " desc = ('Power spectrum data. An array that gives the power spectrum '\n", + " 'values in each kbin. 1D, 2D, or 3D array depending on Ndims '\n", + " 'parameter.')\n", + " self._pspec = prm.psparam('pspec', description=desc, expected_type=np.float, \\\n", + " units='(%s)^2 (%s)^3' % (temp_units,length_units))\n", + " \n", + " desc = ('Central redshift of power spectrum')\n", + " self._z0 = prm.psparam('z0', description=desc, expected_type=np.float)\n", + " \n", + " desc = ('Frequencies used in forming the power spectrum')\n", + " self._freqs = prm.psparam('freqs', description=desc, expected_type=np.float, \\\n", + " units='%s' % freq_units)\n", + " \n", + " desc = ('Window functions. An array where each row corresponds to '\n", + " 'the linear combination of all the other k bins of the '\n", + " 'true power spectrum that is probed when a power spectrum '\n", + " 'estimate is made at a particular k bin.')\n", + " self._window = prm.psparam('window', description=desc, expected_type=np.float)\n", + " \n", + " super(pspec, self).__init__()\n", + " \n", + " def set_pspec():\n", + " \n", + " \n", + " def bin_centers(self):\n", + " k_centers = []\n", + " for i in range(self._Ndims):\n", + " k_centers.append(np.mean(self._kbounds[i], axis=1))\n", + " return k_centers\n", + " \n", + " def DeltaSq(self):\n", + " k_mag = np.zeros_like(self._pspec)\n", + " pspec = np.atleast_3d(self._pspec)\n", + " if self._Ndims == 1:\n", + " pspec = np.rollaxis(pspec,0,-1)\n", + " \n", + " raise NotImplementedError()\n", + " # returns Delta^2(k)\n", + " \n", + " def in_units(self):\n", + " raise NotImplementedError()\n", + " # returns a power spectrum in different units\n", + " \n", + " def rebin(self):\n", + " raise NotImplementedError()\n", + " # rebins the pspec without changing dimensionality, produces a new object\n", + " \n", + " def cylindrical(self):\n", + " raise NotImplementedError()\n", + " # returns a P(kperp, kpara) if available or does a uniform summing.\n", + " \n", + " def spherical(self):\n", + " raise NotImplementedError()\n", + " # same as cylindrical\n", + " \n", + " def reweight(self):\n", + " raise NotImplementedError()\n", + " # applying a new M matrix, returns a new object\n", + " \n", + " \n", + " \n", + " \n", + "# def check windows" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "_Ndims\n", + "_Nks\n", + "_freqs\n", + "_kbounds\n", + "_pspec\n", + "_window\n", + "_z0\n" + ] + } + ], + "source": [ + "gulu = pspec(3, length_units='Mpc/h', temp_units='mK', freq_units='GHz')" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "ename": "NotImplementedError", + "evalue": "", + "output_type": "error", + "traceback": [ + "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", + "\u001b[0;31mNotImplementedError\u001b[0m Traceback (most recent call last)", + "\u001b[0;32m\u001b[0m in \u001b[0;36m\u001b[0;34m()\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mgulu\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mDeltaSq\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m", + "\u001b[0;32m\u001b[0m in \u001b[0;36mDeltaSq\u001b[0;34m(self)\u001b[0m\n\u001b[1;32m 47\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 48\u001b[0m \u001b[0;32mdef\u001b[0m \u001b[0mDeltaSq\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 49\u001b[0;31m \u001b[0;32mraise\u001b[0m \u001b[0mNotImplementedError\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 50\u001b[0m \u001b[0;31m# returns Delta^2(k)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 51\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n", + "\u001b[0;31mNotImplementedError\u001b[0m: " + ] + } + ], + "source": [ + "gulu.DeltaSq()" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "_Ndims\n", + "_Nks\n", + "_freqs\n", + "_kbounds\n", + "_pspec\n", + "_wind\n" + ] + } + ], + "source": [ + "for hey in gulu:\n", + " print hey" + ] + }, + { + "cell_type": "code", + "execution_count": 16, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 16, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "gulu._Ndims" + ] + }, + { + "cell_type": "code", + "execution_count": 16, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "class C(object):\n", + " def __init__(self):\n", + " self._x = None\n", + " setattr(self.__class__, attr_name, property(self.prop_fget(p),self.prop_fset(p)))\n", + "\n", + " def getx(self):\n", + " return self._x\n", + "\n", + " def setx(self, value):\n", + " self._x = value\n", + "\n", + " def delx(self):\n", + " del self._x\n", + " \n", + "\n", + "\n", + "# x = property(getx, setx, delx, \"I'm the 'x' property.\")\n", + " \n", + " def prop_fget(self, param_name):\n", + " \"\"\"Getter method for UVParameter properties.\"\"\"\n", + " def fget(self):\n", + " this_param = getattr(self, param_name)\n", + " return this_param.value\n", + " return fget\n", + "\n", + " def prop_fset(self, param_name):\n", + " \"\"\"Setter method for UVParameter properties.\"\"\"\n", + " def fset(self, value):\n", + " this_param = getattr(self, param_name)\n", + " this_param.value = value\n", + " setattr(self, param_name, this_param)\n", + " return fset" + ] + }, + { + "cell_type": "code", + "execution_count": 25, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "ename": "AttributeError", + "evalue": "'str' object has no attribute 'kbounds'", + "output_type": "error", + "traceback": [ + "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", + "\u001b[0;31mAttributeError\u001b[0m Traceback (most recent call last)", + "\u001b[0;32m\u001b[0m in \u001b[0;36m\u001b[0;34m()\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mgulu\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mprop_fget\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'kbounds'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'kbounds'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m", + "\u001b[0;32m/Users/acliu/Research/pspecPipelineHERA/hera_pspec/transiting_code/file_format/pspecbase.py\u001b[0m in \u001b[0;36mfget\u001b[0;34m(self)\u001b[0m\n\u001b[1;32m 28\u001b[0m \u001b[0;34m\"\"\"Getter method for UVParameter properties.\"\"\"\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 29\u001b[0m \u001b[0;32mdef\u001b[0m \u001b[0mfget\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 30\u001b[0;31m \u001b[0mthis_param\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mgetattr\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mparam_name\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 31\u001b[0m \u001b[0;32mreturn\u001b[0m \u001b[0mthis_param\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mvalue\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 32\u001b[0m \u001b[0;32mreturn\u001b[0m \u001b[0mfget\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", + "\u001b[0;31mAttributeError\u001b[0m: 'str' object has no attribute 'kbounds'" + ] + } + ], + "source": [ + "gulu.prop_fget('kbounds')('kbounds')" + ] + }, + { + "cell_type": "code", + "execution_count": 17, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "property?" + ] + }, + { + "cell_type": "code", + "execution_count": 20, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "(4, 1, 1)" + ] + }, + "execution_count": 20, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "np.rollaxis(np.atleast_3d(np.arange(4)),0,-1).shape" + ] + }, + { + "cell_type": "code", + "execution_count": 19, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "(2, 3, 1)" + ] + }, + "execution_count": 19, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "np.atleast_3d(np.array([[1,2,8],[3,4,2]])).shape" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 2", + "language": "python", + "name": "python2" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 2 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython2", + "version": "2.7.13" + } + }, + "nbformat": 4, + "nbformat_minor": 0 +} diff --git a/transiting_code/file_format/tests/__init__.py b/transiting_code/file_format/tests/__init__.py new file mode 100644 index 00000000..e69de29b diff --git a/transiting_code/file_format/tests/test_pspec.py b/transiting_code/file_format/tests/test_pspec.py new file mode 100644 index 00000000..b957600f --- /dev/null +++ b/transiting_code/file_format/tests/test_pspec.py @@ -0,0 +1,50 @@ +""" +Tests for pspec power spectrum object. +""" +import nose.tools as nt +import numpy as np +import pspec + +def setup_package(): + """ + Construct simple power spectrum and window function data structures. + """ + # Trivial binning and power spectrum data + kbins = np.logspace(-2., 0., 20) + pspec3d = np.ones((kbins.size-1, kbins.size-1, kbins.size-1)) + pspec2d = np.ones((kbins.size-1, kbins.size-1)) + pspec1d = np.ones((kbins.size-1,)) + + +def test_units(): + """ + Make sure that valid units are specified for power spectrum data. + """ + # Errors should be raised if full set of units is not specified + nt.assert_raises(KeyError, pspec.pspec, 3) + nt.assert_raises(KeyError, pspec.pspec, 3, 'Mpc') + nt.assert_raises(KeyError, pspec.pspec, 3, 'Mpc', 'mK') + nt.assert_raises(KeyError, pspec.pspec, 3, length_units='Mpc') + nt.assert_raises(KeyError, pspec.pspec, 3, None, None, None) + nt.assert_raises(KeyError, pspec.pspec, 3, None, 'mK', 'GHz') + nt.assert_raises(KeyError, pspec.pspec, 3, '', 'mK', 'GHz') + nt.assert_raises(KeyError, pspec.pspec, 3, 'Mpc', '', 'GHz') + nt.assert_raises(KeyError, pspec.pspec, 3, 'Mpc', 'mK', '') + nt.assert_raises(KeyError, pspec.pspec, 3, ' ', 'mK', 'GHz') + nt.assert_raises(KeyError, pspec.pspec, 3, ' ', 'mK', 'GHz') + nt.assert_raises(KeyError, pspec.pspec, 3, ['Mpc',], 'mK', 'GHz') + nt.assert_raises(KeyError, pspec.pspec, 3, ('Mpc',), 'mK', 'GHz') + + # Make sure valid instantiations of the pspec class work as expected + nt.assert_true(isinstance(pspec.pspec(3, 'Mpc', 'mK', 'GHz'), pspec.pspec)) + nt.assert_true(isinstance(pspec.pspec(3, 'Mpc/h', 'mK', 'GHz'), pspec.pspec)) + nt.assert_true(isinstance(pspec.pspec(3, 'h^-1 Mpc', 'mK', 'GHz'), pspec.pspec)) + nt.assert_true(isinstance(pspec.pspec(3, 'Mpc', 'mK', 'ghz'), pspec.pspec)) + + # Check that power spectrum units are correct + units = pspec.pspec(3, 'Mpc', 'mK', 'GHz')._pspec.units + nt.assert_true('(mK)^2' in units) + nt.assert_true('(Mpc)^3' in units) + +if __name__ == '__main__': + test_units() diff --git a/transiting_code/oqe.py b/transiting_code/oqe.py index e677be42..28c97174 100644 --- a/transiting_code/oqe.py +++ b/transiting_code/oqe.py @@ -248,7 +248,9 @@ def gen_bl_boots(self, nboots, ngps=5): def gen_gps(self, bls, ngps=5): random.shuffle(bls) gps = [bls[i::ngps] for i in range(ngps)] - gps = [[random.choice(gp) for bl in gp] for gp in gps] #sample w/replacement inside each group + #gps = [[random.choice(gp) for bl in gp] for gp in gps] #sample w/replacement inside each group + # All independent baselines within a group, except the last one which is sampled randomly + gps = [[gps[gp][i] for i in np.append(np.arange(0,len(gps[gp])-1), random.choice(np.arange(0,len(gps[gp]))))] for gp,G in enumerate(gps)] return gps def group_data(self, keys, gps, use_cov=True): #XXX keys have format (k,bl,POL) # XXX avoid duplicate code for use_cov=True vs False (i.e. no separate dsC & dsI) @@ -262,13 +264,26 @@ def group_data(self, keys, gps, use_cov=True): #XXX keys have format (k,bl,POL) for gp in range(len(gps)): newkey = (s,gp) newkeys.append(newkey) + """ # Jack-Knife Test: change sign of baselines going into each group + gps[gp] = gps[gp][:(len(gps[gp])/2)*2] # ensure even number of baselines per group + factors = [] + for b in range(len(gps[gp])): + if b < len(gps[gp])/2: #for half of the baselines + #self.add_data(dsets={(s,gps[gp][b],POL): -1*self.x[(s,gps[gp][b],POL)].T}) + factors.append(-1) + else: factors.append(1) + iCsum[newkey] = sum([self.iC((s,bl,POL)) for bl in gps[gp]]) + iCxsum[newkey] = sum([np.dot(self.iC((s,bl,POL)),factors[i]*self.x[(s,bl,POL)]) for i,bl in enumerate(gps[gp])]) + Isum[newkey] = sum([np.identity(nchan) for bl in gps[gp]]) + Ixsum[newkey] = sum([factors[i]*self.x[(s,bl,POL)] for i,bl in enumerate(gps[gp])]) + """ iCsum[newkey] = sum([self.iC((s,bl,POL)) for bl in gps[gp]]) iCxsum[newkey] = sum([np.dot(self.iC((s,bl,POL)),self.x[(s,bl,POL)]) for bl in gps[gp]]) Isum[newkey] = sum([np.identity(nchan) for bl in gps[gp]]) Ixsum[newkey] = sum([self.x[(s,bl,POL)] for bl in gps[gp]]) dsC_data[newkey] = np.dot(np.linalg.inv(iCsum[newkey]),iCxsum[newkey]).T #finding effective summed up x based on iCsum and iCxsum dsI_data[newkey] = np.dot(np.linalg.inv(Isum[newkey]),Ixsum[newkey]).T #finding effective summed up x based on Isum and Ixsum - dsC = DataSet(); dsC.add_data(dsets=dsC_data) + dsC = DataSet(); dsC.set_data(dsets=dsC_data) dsI = DataSet(); dsI.set_data(dsets=dsI_data) #I has to be a separate dataset because it has different x's populated into it dsC.set_iC(iCsum) #override since if they're computed from x, they're incorrect dsI.set_I(Isum)