From ace5ff9d916d64ac1c52894029c218d70ba764e3 Mon Sep 17 00:00:00 2001 From: Carolin Hofer Date: Thu, 2 Jul 2020 11:30:28 -0700 Subject: [PATCH 1/4] feat(kltransform): allow F/S diagonalisation order (cherry picked from commit 3af5152218f19c271b0e26d7e1d8e6523f6f8772) --- drift/core/doublekl.py | 19 ++++++++++-- drift/core/kltransform.py | 62 +++++++++++++++++++++++++++++++-------- 2 files changed, 65 insertions(+), 16 deletions(-) diff --git a/drift/core/doublekl.py b/drift/core/doublekl.py index b3dc958d..c98a2abb 100644 --- a/drift/core/doublekl.py +++ b/drift/core/doublekl.py @@ -28,6 +28,7 @@ class DoubleKL(kltransform.KLTransform): """ foreground_threshold = config.Property(proptype=float, default=100.0) + _dkl = True def _transform_m(self, mi): @@ -49,11 +50,18 @@ def _transform_m(self, mi): cs, cn = [cv.reshape(nside, nside) for cv in self.sn_covariance(mi)] # Find joint eigenbasis and transformation matrix - evals, evecs2, ac = kltransform.eigh_gen(cs, cn) + if self.diagonalisation_order == "sf": + evals, evecs2, ac = kltransform.eigh_gen(cs, cn) + else: + evals, evecs2, ac = kltransform.eigh_gen(cn, cs) + evecs = evecs2.T.conj() # Get the indices that extract the high S/F ratio modes - ind = np.where(evals > self.foreground_threshold) + if self.diagonalisation_order == "sf": + ind = np.where(evals > self.foreground_threshold) + else: + ind = np.where(evals < self.foreground_threshold) # Construct evextra dictionary (holding foreground ratio) evextra = {"ac": ac, "f_evals": evals.copy()} @@ -90,7 +98,12 @@ def _ev_save_hook(self, f, evextra): kltransform.KLTransform._ev_save_hook(self, f, evextra) # Save out S/F ratios - f.create_dataset("f_evals", data=evextra["f_evals"]) + if self.diagonalisation_order == "sf": + f_evals = evextra["f_evals"] + else: + f_evals = evextra["f_evals"][::-1] + + f.create_dataset("f_evals", data=f_evals) def _collect(self): def evfunc(mi): diff --git a/drift/core/kltransform.py b/drift/core/kltransform.py index 147f8dd0..72a9e869 100644 --- a/drift/core/kltransform.py +++ b/drift/core/kltransform.py @@ -173,6 +173,8 @@ class KLTransform(config.Reader): threshold = config.Property(proptype=float, default=0.1, key="threshold") + diagonalisation_order = config.enum(["sf", "fs"], default="sf") + _foreground_regulariser = config.Property( proptype=float, default=1e-14, key="regulariser" ) @@ -187,6 +189,7 @@ class KLTransform(config.Reader): _cvfg = None _cvsg = None + _dkl = None @property def _evfile(self): @@ -346,7 +349,12 @@ def _transform_m(self, mi): # Perform the generalised eigenvalue problem to get the KL-modes. st = time.time() - evals, evecs, ac = eigh_gen(cvb_sr, cvb_nr) + + if self.diagonalisation_order == "sf": + evals, evecs, ac = eigh_gen(cvb_sr, cvb_nr) + else: + evals, evecs, ac = eigh_gen(cvb_nr, cvb_sr) + et = time.time() print("Time =", (et - st)) @@ -396,18 +404,26 @@ def transform_save(self, mi): evalsf = np.zeros(nside, dtype=np.float64) if evals.size != 0: evalsf[(-evals.size) :] = evals + f.create_dataset("evals_full", data=evalsf) # Discard eigenmodes with S/N below threshold if requested. if self.subset: i_ev = np.searchsorted(evals, self.threshold) - - evals = evals[i_ev:] - evecs = evecs[i_ev:] - print( - "Modes with S/N > %f: %i of %i" - % (self.threshold, evals.size, evalsf.size) - ) + if self.diagonalisation_order == "sf" or self._dkl: + evals = evals[i_ev:] + evecs = evecs[i_ev:] + print( + "Modes with S/N > %f: %i of %i" + % (self.threshold, evals.size, evalsf.size) + ) + else: + evals = evals[:i_ev] + evecs = evecs[:i_ev] + print( + "Modes with N/S < %f: %i of %i" + % (self.threshold, evals.size, evalsf.size) + ) # Write out potentially reduced eigen spectrum. f.create_dataset("evals", data=evals) @@ -416,7 +432,10 @@ def transform_save(self, mi): if self.inverse: if self.subset: - inv = inv[i_ev:] + if self.diagonalisation_order == "sf" or self._dkl: + inv = inv[i_ev:] + else: + inv = inv[:i_ev] f.create_dataset("evinv", data=inv) @@ -466,6 +485,7 @@ def evfunc(mi): if f["evals_full"].shape[0] > 0: ev = f["evals_full"][:] evf[-ev.size :] = ev + f.close() return evf @@ -500,6 +520,9 @@ def generate(self, regen=False): if mpiutil.rank0: st = time.time() print("======== Starting KL calculation ========") + print( + "======= Diagonalisation order: %s =======" % self.diagonalisation_order + ) # Iterate list over MPI processes. for mi in mpiutil.mpirange(self.telescope.mmax + 1): @@ -571,7 +594,10 @@ def modes_m(self, mi, threshold=None): if startind == evals.size: modes = None, None else: - modes = (evals[startind:], f["evecs"][startind:]) + if self.diagonalisation_order == "sf" or self._dkl: + modes = (evals[startind:], f["evecs"][startind:]) + else: + modes = (evals[:startind], f["evecs"][:startind]) # If old data file perform complex conjugate modes = ( @@ -628,7 +654,10 @@ def evals_m(self, mi, threshold=None): if startind == evals.size: modes = None else: - modes = evals[startind:] + if self.diagonalisation_order == "sf" or self._dkl: + modes = evals[startind:] + else: + mode = evals[:startind] f.close() @@ -661,7 +690,10 @@ def invmodes_m(self, mi, threshold=None): if threshold != None: nevals = evals.size - inv = inv[(-nevals):] + if self.diagonalisation_order == "sf" or self._dkl: + inv = inv[(-nevals):] + else: + inv = inv[:nevals] return inv.T @@ -896,7 +928,11 @@ def project_sky(self, sky, mlist=None, threshold=None, harmonic=False): def _proj(mi): p1 = self.project_sky_vector_forward(mi, alm[:, :, mi], threshold) p2 = np.zeros(nmodes, dtype=np.complex128) - p2[-p1.size :] = p1 + if self.diagonalisation_order == "sf" or self._dkl: + p2[-p1.size :] = p1 + else: + p2[: p1.size] = p1 + return p2 # Map over list of m's and project sky onto eigenbasis From b7733b173cc3c93fb9bab18a8f5ebf42adec34c0 Mon Sep 17 00:00:00 2001 From: Carolin Hofer Date: Mon, 1 Feb 2021 19:27:04 -0800 Subject: [PATCH 2/4] fix(doublekl): black --- drift/core/doublekl.py | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/drift/core/doublekl.py b/drift/core/doublekl.py index fd6cf81a..e4102346 100644 --- a/drift/core/doublekl.py +++ b/drift/core/doublekl.py @@ -44,9 +44,13 @@ def _transform_m(self, mi): # Find joint eigenbasis and transformation matrix if self.diagonalisation_order == "sf": - evals, evecs2, ac = kltransform.eigh_gen(cs, cn, message="m = %d; KL step 1" % mi) + evals, evecs2, ac = kltransform.eigh_gen( + cs, cn, message="m = %d; KL step 1" % mi + ) else: - evals, evecs2, ac = kltransform.eigh_gen(cn, cs, message="m = %d; KL step 1" % mi) + evals, evecs2, ac = kltransform.eigh_gen( + cn, cs, message="m = %d; KL step 1" % mi + ) evecs = evecs2.T.conj() From 853e72fa44b00f011508dfa0cbfd7368c7c7d110 Mon Sep 17 00:00:00 2001 From: Carolin Hofer Date: Thu, 25 Mar 2021 21:50:25 -0700 Subject: [PATCH 3/4] docs(KLTransform): explain diagonalisation order and _dkl --- drift/core/doublekl.py | 2 +- drift/core/kltransform.py | 6 +++++- 2 files changed, 6 insertions(+), 2 deletions(-) diff --git a/drift/core/doublekl.py b/drift/core/doublekl.py index e4102346..435d98c9 100644 --- a/drift/core/doublekl.py +++ b/drift/core/doublekl.py @@ -100,7 +100,7 @@ def _ev_save_hook(self, f, evextra): if self.diagonalisation_order == "sf": f_evals = evextra["f_evals"] else: - f_evals = evextra["f_evals"][::-1] + f_evals = evextra["f_evals"] f.create_dataset("f_evals", data=f_evals) diff --git a/drift/core/kltransform.py b/drift/core/kltransform.py index b7955e4e..a6fa4470 100644 --- a/drift/core/kltransform.py +++ b/drift/core/kltransform.py @@ -154,6 +154,9 @@ class KLTransform(config.Reader): If True, throw away modes below a S/N `threshold`. threshold : scalar S/N threshold to cut modes at. + diagonalisation_order : string, optional + Options are 'sf' for eigen value decomposition with S/F eigen values and + 'fs' for F/S eigen values. inverse : boolean If True construct and cache inverse transformation. use_thermal, use_foregrounds : boolean @@ -405,6 +408,7 @@ def transform_save(self, mi): # Discard eigenmodes with S/N below threshold if requested. if self.subset: i_ev = np.searchsorted(evals, self.threshold) + # The second step in the double KL will be eigen values of S/N regardless of first step KL step. We save the double KL values in S/N order. if self.diagonalisation_order == "sf" or self._dkl: evals = evals[i_ev:] evecs = evecs[i_ev:] @@ -652,7 +656,7 @@ def evals_m(self, mi, threshold=None): if self.diagonalisation_order == "sf" or self._dkl: modes = evals[startind:] else: - mode = evals[:startind] + modes = evals[:startind] f.close() From 14b571cdeaddde568b69ab008390352e47351121 Mon Sep 17 00:00:00 2001 From: Carolin Hofer Date: Sun, 28 Mar 2021 17:04:34 -0700 Subject: [PATCH 4/4] fix(doublekl): remove unnecessary if/else statement --- drift/core/doublekl.py | 7 +------ 1 file changed, 1 insertion(+), 6 deletions(-) diff --git a/drift/core/doublekl.py b/drift/core/doublekl.py index 435d98c9..af34f0aa 100644 --- a/drift/core/doublekl.py +++ b/drift/core/doublekl.py @@ -97,12 +97,7 @@ def _ev_save_hook(self, f, evextra): kltransform.KLTransform._ev_save_hook(self, f, evextra) # Save out S/F ratios - if self.diagonalisation_order == "sf": - f_evals = evextra["f_evals"] - else: - f_evals = evextra["f_evals"] - - f.create_dataset("f_evals", data=f_evals) + f.create_dataset("f_evals", data=evextra["f_evals"]) def _collect(self): def evfunc(mi):