Px data #2

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Draft

andreufont wants to merge 12 commits into main from px_data

cupix/px_data/data_healpix.py

-Original file line number
+Diff line change
@@ -0,0 +1,160 @@
+    import numpy as np
+    import h5py
+    from cupix.px_data import read_healpix_px, px_window
+    class HealpixPxArchive(object):
+        '''Collection of PX measurements from healpixels'''
+        def __init__(self, fname, list_hp=None, list_px=None,
+                        compute_window=False):
+            '''Setup collection, from file or by hand'''
+            if fname is not None:
+                assert (list_hp is None) and (list_px is None)
+                self.fname = fname
+                print('read healpix info from', fname)
+                # read information from file (Sindhu's format)
+                reader = read_healpix_px.HealpixPxReader(fname)
+                # use information to construct list_px
+                self.list_hp, self.list_px = reader.get_list_px(compute_window)
+            else:
+                self.fname = None
+                assert len(list_hp) == len(list_px)
+                self.list_hp = list_hp
+                self.list_px = list_px
+            print(f'Got Px for {len(self.list_hp)} healpixels')
+            # get binning information
+            self.z_bins = self.list_px[0].z_bins
+            self.t_bins = self.list_px[0].t_bins
+            self.k_bins = self.list_px[0].k_bins
+            # length in Angstroms
+            self.L_A = self.list_px[0].L_A
+            # mean pixel resolution in Angstroms (could be z dependent eventually)
+            self.sig_A = self.list_px[0].sig_A
+            for px in self.list_px:
+                assert len(self.z_bins) == len(px.z_bins)
+                assert len(self.t_bins) == len(px.t_bins)
+                assert len(self.k_bins) == len(px.k_bins)
+                assert self.L_A == px.L_A
+                assert self.sig_A == px.sig_A
+            return
+        def get_mean_px(self):
+            list_px_z = []
+            for iz, z_bin in enumerate(self.z_bins):
+                list_px_zt = []
+                for it, t_bin in enumerate(self.t_bins):
+                    # this should be properly computed, not hard-coded
+                    Nk = len(self.k_bins)
+                    # add the contribution from each healpixel
+                    F_m = np.zeros(Nk)
+                    W_m = np.zeros(Nk)
+                    for px in self.list_px:
+                        px_zt = px.list_px_z[iz].list_px_zt[it]
+                        F_m += px_zt.F_m
+                        W_m += px_zt.W_m
+                    px_zt = px_window.Px_zt_w.from_unnormalized(
+                                    z_bin, t_bin, self.k_bins,
+                                    F_m=F_m, W_m=W_m,
+                                    L_A=self.L_A, sig_A=self.sig_A,
+                                    compute_window=False)
+                    list_px_zt.append(px_zt)
+                px_z = px_window.Px_z_w(self.t_bins, list_px_zt)
+                list_px_z.append(px_z)
+            mean_px = px_window.Px_w(self.z_bins, list_px_z)
+            return mean_px
+        def get_mean_and_cov(self, compute_window=False):
+            list_px_z = []
+            for iz, z_bin in enumerate(self.z_bins):
+                list_px_zt = []
+                for it, t_bin in enumerate(self.t_bins):
+                    # collect PX from each healpixel into (N_hp, Nk) array
+                    Nk = len(self.k_bins)
+                    N_hp = len(self.list_hp)
+                    all_P_m = np.empty((N_hp,Nk))
+                    all_V_m = np.empty((N_hp,Nk))
+                    for ipix in range(N_hp):
+                        px_zt = self.list_px[ipix].list_px_z[iz].list_px_zt[it]
+                        all_P_m[ipix] = px_zt.P_m
+                        all_V_m[ipix] = px_zt.V_m
+                    # compute weighted mean and covariance (from Picca)
+                    mean_P_m = (all_P_m * all_V_m).sum(axis=0)
+                    sum_V_m = all_V_m.sum(axis=0)
+                    w = sum_V_m > 0.
+                    mean_P_m[w] /= sum_V_m[w]
+                    meanless_P_V_m = all_V_m * (all_P_m - mean_P_m)
+                    C_mn = meanless_P_V_m.T.dot(meanless_P_V_m)
+                    sum_V2 = sum_V_m * sum_V_m[:, None]
+                    w = sum_V2 > 0.
+                    C_mn[w] /= sum_V2[w]
+                    # check if we need to compute the mean window matrix
+                    if compute_window:
+                        mean_U_mn = self.get_mean_window_zt(iz, it, all_V_m)
+                    else:
+                        mean_U_mn = None
+                    # setup object (normalized, with covariance)
+                    px_zt = px_window.Px_zt_w(z_bin, t_bin, self.k_bins,
+                                    P_m=mean_P_m, V_m=sum_V_m, C_mn=C_mn,
+                                    F_m=None, W_m=None, U_mn=mean_U_mn)
+                    list_px_zt.append(px_zt)
+                px_z = px_window.Px_z_w(self.t_bins, list_px_zt)
+                list_px_z.append(px_z)
+            mean_px = px_window.Px_w(self.z_bins, list_px_z)
+            return mean_px
+        def get_mean_window_zt(self, iz, it):
+            '''Compute mean window matrix, for a given (z,t) bin'''
+            # collect PX from each healpixel into (N_hp, Nk) array
+            Nk = len(self.k_bins)
+            N_hp = len(self.list_hp)
+            all_P_m = np.empty((N_hp, Nk))
+            all_V_m = np.empty((N_hp, Nk))
+            # collect individual window matrices (per healpixel)
+            all_U_mn = np.empty((N_hp, Nk, Nk))
+            for ipix in range(N_hp):
+                px_zt = self.list_px[ipix].list_px_z[iz].list_px_zt[it]
+                all_P_m[ipix] = px_zt.P_m
+                all_V_m[ipix] = px_zt.V_m
+                if px_zt.U_mn is None:
+                    # compute U_mn for each healpixel
+                    R2_m = px_window.compute_R2_m(px_zt.k_bins, px_zt.sig_A)
+                    U_mn = px_window.compute_U_mn(px_zt.W_m, R2_m, px_zt.L_A)
+                    all_U_mn[ipix] = U_mn
+                else:
+                    all_U_mn[ipix] = px_zt.U_mn
+            # compute weighted mean
+            mean_P_m = (all_P_m * all_V_m).sum(axis=0)
+            sum_V_m = all_V_m.sum(axis=0)
+            w = sum_V_m > 0.
+            mean_P_m[w] /= sum_V_m[w]
+            # weighted sum of all the windows
+            sum_U_mn = np.einsum('im,imn->mn', all_V_m, all_U_mn)
+            mean_U_mn = sum_U_mn / sum_V_m[:,None]
+            return sum_V_m, mean_P_m, mean_U_mn
+        def rebin(self, rebin_t_factor, rebin_k_factor, include_k_0=True):
+            '''Return a new Px archive, after rebinning in theta (t) and k'''
+            new_list_px = []
+            for px in self.list_px:
+                new_px = px.rebin(rebin_t_factor, rebin_k_factor, include_k_0)
+                new_list_px.append(new_px)
+            new_archive = HealpixPxArchive(fname=None,
+                                list_hp=self.list_hp, list_px=new_list_px)
+            return new_archive

cupix/px_data/px_binning.py

-Original file line number
+Diff line change
@@ -0,0 +1,161 @@
+    import numpy as np
+    class Bin_z(object):
+        '''Description of a particular redshift bin'''
+        def __init__(self, min_z, max_z, label=None):
+            '''Construct by providing edges, and plotting label'''
+            self.min_z = min_z
+            self.max_z = max_z
+            if label is None:
+                self.label = f'{min_z} < z < {max_z}'
+            else:
+                self.label = label
+        def mean(self):
+            return 0.5*(self.min_z+self.max_z)
+    class Bin_t(object):
+        '''Description of a particular theta bin'''
+        def __init__(self, min_t, max_t, label=None):
+            '''Construct by providing edges, and plotting label'''
+            self.min_t = min_t
+            self.max_t = max_t
+            if label is None:
+                self.label = f'{min_t} < theta [arcmin] < {max_t}'
+            else:
+                self.label = label
+        def mean(self):
+            return 0.5*(self.min_t+self.max_t)
+        def B_t(self, t):
+            '''Contribution of theta value to this bin'''
+            # this could be made more efficient by allowing arrays, if needed
+            if (t >= self.min_t) and (t < self.max_t):
+                return 1.0
+            else:
+                return 0.0
+        def B_t_bins(self, t_bins):
+            '''Same than B_t for a list of narrow theta bins'''
+            t_a = np.array( [ t_bin.mean() for t_bin in t_bins ])
+            B_a = np.zeros_like(t_a)
+            in_bin = [ (t >= self.min_t) and (t < self.max_t) for t in t_a ]
+            B_a[in_bin] = 1.0
+            return B_a
+    class Discrete_k(object):
+        '''Description of a particular wavenumber k (no binning)'''
+        def __init__(self, k, label=None):
+            '''Construct by providing value and plotting label'''
+            self.k = k
+            if label is None:
+                self.label = f'k = {k} 1/A'
+            else:
+                self.label = label
+    class Bin_k(object):
+        '''Description of a particular k bin'''
+        def __init__(self, min_k, max_k, label=None):
+            '''Construct by providing edges, and plotting label'''
+            self.min_k = min_k
+            self.max_k = max_k
+            if label is None:
+                self.label = f'{min_k} < k [1/A] < {max_k}'
+            else:
+                self.label = label
+        def mean(self):
+            return 0.5*(self.min_k+self.max_k)
+        def B_k(self, k):
+            '''Contribution of k value to this k-bin'''
+            # this could be made more efficient by allowing arrays, if needed
+            if (k >= self.min_k) and (k < self.max_k):
+                return 1.0
+            else:
+                return 0.0
+        def B_k_bins(self, k_bins):
+            '''Same than B_k, for list of narrow bins'''
+            k_m = np.array( [ k_bin.k for k_bin in k_bins ])
+            B_m = np.zeros_like(k_m)
+            in_bin = [ (k >= self.min_k) and (k < self.max_k) for k in k_m ]
+            B_m[in_bin] = 1.0
+            return B_m
+    def get_coarser_k_bins(in_k_bins, rebin_factor, include_k_0=True):
+        '''Return coarser bins, rebin_factor coarser than in_k_bins.
+           If include_k_0=True, include k=0 instead of k_Ny.'''
+        # this function for now should be called on Discrete_k bins
+        assert isinstance(in_k_bins[0], Discrete_k)
+        in_k = [ k_bin.k for k_bin in in_k_bins ]
+        in_dk = in_k[1]-in_k[0]
+        in_Nk = len(in_k)
+        in_k_max = np.max(np.abs(in_k))
+        #print(f'input Nk = {in_Nk}, dk = {in_dk:.4f}, k_max = {in_k_max:.4f}')
+        # extra factor of 2 here since in_k_bins include negative frequencies
+        out_Nk = int(in_Nk / 2 / rebin_factor)
+        if include_k_0:
+            # include k=0 in first bin
+            k_shift = -0.5*in_dk
+        else:
+            # include k_Ny in last bin
+            k_shift = +0.5*in_dk
+        out_k_edges = np.linspace(k_shift, in_k_max+k_shift, num = out_Nk+1, endpoint=True)
+        out_dk = out_k_edges[1] - out_k_edges[0]
+        out_k_max = np.max(np.abs(out_k_edges))
+        #print(f'out Nk = {out_Nk}, dk = {out_dk:.4f}, k_max = {out_k_max:.4f}')
+        out_k_bins = []
+        for ik in range(out_Nk):
+            k_bin = Bin_k(out_k_edges[ik], out_k_edges[ik+1])
+            out_k_bins.append(k_bin)
+        return out_k_bins
+    def get_coarser_t_bins(in_t_bins, rebin_factor):
+        '''Return coarser bins, rebin_factor coarser than in_t_bins'''
+        in_Nt = len(in_t_bins)
+        assert in_Nt % rebin_factor == 0, 'size mismatch'
+        in_min_t = [ t_bin.min_t for t_bin in in_t_bins ]
+        in_max_t = [ t_bin.max_t for t_bin in in_t_bins ]
+        out_min_t = in_min_t[::rebin_factor]
+        out_max_t = in_max_t[rebin_factor-1::rebin_factor]
+        out_Nt = len(out_min_t)
+        out_t_bins = []
+        for it in range(out_Nt):
+            #print(f'{it}, {out_min_t[it]} < t < {out_max_t[it]}')
+            t_bin = Bin_t(out_min_t[it], out_max_t[it])
+            out_t_bins.append(t_bin)
+        return out_t_bins

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Px data #2

Uh oh!

Diff view

Diff view

There are no files selected for viewing

Uh oh!

Px data #2

Are you sure you want to change the base?

Uh oh!

Px data #2

Uh oh!

Uh oh!

Diff view

Diff view

There are no files selected for viewing

Uh oh!