From f71278d807eaf360e3195739215aa9b6fbd96e72 Mon Sep 17 00:00:00 2001 From: "Paul T. Baker" Date: Thu, 24 Jun 2021 14:22:37 -0400 Subject: [PATCH 1/4] add options to return efac, equad, ecorr vectors --- libstempo/toasim.py | 18 +++++++++++++++--- 1 file changed, 15 insertions(+), 3 deletions(-) diff --git a/libstempo/toasim.py b/libstempo/toasim.py index 66f43ed..fe9d5fb 100644 --- a/libstempo/toasim.py +++ b/libstempo/toasim.py @@ -159,7 +159,8 @@ def make_ideal(psr): psr.fit() -def add_efac(psr, efac=1.0, flagid=None, flags=None, seed=None): +def add_efac(psr, efac=1.0, flagid=None, flags=None, + return_vec=False, seed=None): """Add nominal TOA errors, multiplied by `efac` factor. Optionally take a pseudorandom-number-generator seed.""" @@ -184,8 +185,12 @@ def add_efac(psr, efac=1.0, flagid=None, flags=None, seed=None): psr.stoas[:] += efacvec * psr.toaerrs * (1e-6 / day) * N.random.randn(psr.nobs) + if return_vec: + return efacvec -def add_equad(psr, equad, flagid=None, flags=None, seed=None): + +def add_equad(psr, equad, flagid=None, flags=None, + return_vec=False, seed=None): """Add quadrature noise of rms `equad` [s]. Optionally take a pseudorandom-number-generator seed.""" @@ -210,6 +215,9 @@ def add_equad(psr, equad, flagid=None, flags=None, seed=None): psr.stoas[:] += (equadvec / day) * N.random.randn(psr.nobs) + if return_vec: + return equadvec + def quantize(times, dt=1): bins = N.arange(N.min(times), N.max(times) + dt, dt) @@ -264,7 +272,8 @@ def quantize_fast(times, flags=None, dt=1.0): # print N.sum((t - t2)**2), N.all(U == U2) -def add_jitter(psr, ecorr, flagid=None, flags=None, coarsegrain=0.1, seed=None): +def add_jitter(psr, ecorr, flagid=None, flags=None, coarsegrain=0.1, + return_vec=False, seed=None): """Add correlated quadrature noise of rms `ecorr` [s], with coarse-graining time `coarsegrain` [days]. Optionally take a pseudorandom-number-generator seed.""" @@ -295,6 +304,9 @@ def add_jitter(psr, ecorr, flagid=None, flags=None, coarsegrain=0.1, seed=None): psr.stoas[:] += (1 / day) * N.dot(U * ecorrvec, N.random.randn(U.shape[1])) + if return_vec: + return np.dot(U * ecorrvec, np.ones(U.shape[1])) + def add_rednoise(psr, A, gamma, components=10, seed=None): """Add red noise with P(f) = A^2 / (12 pi^2) (f year)^-gamma, From fb1e3467b0124fa2cd26645344409c850aef1e38 Mon Sep 17 00:00:00 2001 From: "Paul T. Baker" Date: Thu, 24 Jun 2021 14:29:57 -0400 Subject: [PATCH 2/4] add epoch averaging function --- libstempo/toasim.py | 50 +++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 50 insertions(+) diff --git a/libstempo/toasim.py b/libstempo/toasim.py index fe9d5fb..ca3d17a 100644 --- a/libstempo/toasim.py +++ b/libstempo/toasim.py @@ -937,3 +937,53 @@ def computeORFMatrix(psr): ORF[ll, kk] = 2.0 return ORF + +def compute_epoch_ave(times, res, err, ecorr=None, dt=1.0, flags=None): + """Compute epoch averaged TOAs including uncertainty + + :param times: TOAs in seconds + :param res: Residuals in seconds + :param err: Scaled (by EFAC and EQUAD) error bars in seconds + :param ecorr: (optional) ECORR value for each point in s^2 [default None] + :param dt: (optional) Time bins for averaging [default 1 s] + + :returns: average TOAs, error bars, residuals in seconds + :returns: (optional) average flags + + from PAL2 (Ellis & van Haasteren 2017) + https://doi.org/10.5281/zenodo.251456 + """ + isort = np.argsort(times) + + bucket_ref = [times[isort[0]]] + bucket_ind = [[isort[0]]] + + for i in isort[1:]: + if times[i] - bucket_ref[-1] < dt: + bucket_ind[-1].append(i) + else: + bucket_ref.append(times[i]) + bucket_ind.append([i]) + + avetoas = np.array([np.mean(times[l]) for l in bucket_ind],'d') + + if flags is not None: + aveflags = np.array([flags[l[0]] for l in bucket_ind]) + + aveerr = np.zeros(len(bucket_ind)) + averes = np.zeros(len(bucket_ind)) + + for i,l in enumerate(bucket_ind): + M = np.ones(len(l)) + C = np.diag(err[l]**2) + if ecorr is not None: + C += np.ones((len(l), len(l))) * ecorr[l[0]] + + avr = 1/np.dot(M, np.dot(np.linalg.inv(C), M)) + aveerr[i] = np.sqrt(avr) + averes[i] = avr * np.dot(M, np.dot(np.linalg.inv(C), res[l])) + + if flags is not None: + return avetoas, aveerr, averes, aveflags + else: + return avetoas, aveerr, averes From 6c734a19aa10d53817db18507afb24f8366e00c1 Mon Sep 17 00:00:00 2001 From: "Paul T. Baker" Date: Thu, 24 Jun 2021 14:35:06 -0400 Subject: [PATCH 3/4] np -> N --- libstempo/toasim.py | 24 ++++++++++++------------ 1 file changed, 12 insertions(+), 12 deletions(-) diff --git a/libstempo/toasim.py b/libstempo/toasim.py index ca3d17a..0cd1f31 100644 --- a/libstempo/toasim.py +++ b/libstempo/toasim.py @@ -305,7 +305,7 @@ def add_jitter(psr, ecorr, flagid=None, flags=None, coarsegrain=0.1, psr.stoas[:] += (1 / day) * N.dot(U * ecorrvec, N.random.randn(U.shape[1])) if return_vec: - return np.dot(U * ecorrvec, np.ones(U.shape[1])) + return N.dot(U * ecorrvec, N.ones(U.shape[1])) def add_rednoise(psr, A, gamma, components=10, seed=None): @@ -953,7 +953,7 @@ def compute_epoch_ave(times, res, err, ecorr=None, dt=1.0, flags=None): from PAL2 (Ellis & van Haasteren 2017) https://doi.org/10.5281/zenodo.251456 """ - isort = np.argsort(times) + isort = N.argsort(times) bucket_ref = [times[isort[0]]] bucket_ind = [[isort[0]]] @@ -965,23 +965,23 @@ def compute_epoch_ave(times, res, err, ecorr=None, dt=1.0, flags=None): bucket_ref.append(times[i]) bucket_ind.append([i]) - avetoas = np.array([np.mean(times[l]) for l in bucket_ind],'d') + avetoas = N.array([N.mean(times[l]) for l in bucket_ind],'d') if flags is not None: - aveflags = np.array([flags[l[0]] for l in bucket_ind]) + aveflags = N.array([flags[l[0]] for l in bucket_ind]) - aveerr = np.zeros(len(bucket_ind)) - averes = np.zeros(len(bucket_ind)) + aveerr = N.zeros(len(bucket_ind)) + averes = N.zeros(len(bucket_ind)) for i,l in enumerate(bucket_ind): - M = np.ones(len(l)) - C = np.diag(err[l]**2) + M = N.ones(len(l)) + C = N.diag(err[l]**2) if ecorr is not None: - C += np.ones((len(l), len(l))) * ecorr[l[0]] + C += N.ones((len(l), len(l))) * ecorr[l[0]] - avr = 1/np.dot(M, np.dot(np.linalg.inv(C), M)) - aveerr[i] = np.sqrt(avr) - averes[i] = avr * np.dot(M, np.dot(np.linalg.inv(C), res[l])) + avr = 1/N.dot(M, N.dot(N.linalg.inv(C), M)) + aveerr[i] = N.sqrt(avr) + averes[i] = avr * N.dot(M, N.dot(N.linalg.inv(C), res[l])) if flags is not None: return avetoas, aveerr, averes, aveflags From b0ee7f2133b099f4538073dcf5073262b6e97fc7 Mon Sep 17 00:00:00 2001 From: "Paul T. Baker" Date: Thu, 24 Jun 2021 14:40:42 -0400 Subject: [PATCH 4/4] fix linting --- libstempo/toasim.py | 18 ++++++++---------- 1 file changed, 8 insertions(+), 10 deletions(-) diff --git a/libstempo/toasim.py b/libstempo/toasim.py index 0cd1f31..811d1ad 100644 --- a/libstempo/toasim.py +++ b/libstempo/toasim.py @@ -159,8 +159,7 @@ def make_ideal(psr): psr.fit() -def add_efac(psr, efac=1.0, flagid=None, flags=None, - return_vec=False, seed=None): +def add_efac(psr, efac=1.0, flagid=None, flags=None, return_vec=False, seed=None): """Add nominal TOA errors, multiplied by `efac` factor. Optionally take a pseudorandom-number-generator seed.""" @@ -189,8 +188,7 @@ def add_efac(psr, efac=1.0, flagid=None, flags=None, return efacvec -def add_equad(psr, equad, flagid=None, flags=None, - return_vec=False, seed=None): +def add_equad(psr, equad, flagid=None, flags=None, return_vec=False, seed=None): """Add quadrature noise of rms `equad` [s]. Optionally take a pseudorandom-number-generator seed.""" @@ -272,8 +270,7 @@ def quantize_fast(times, flags=None, dt=1.0): # print N.sum((t - t2)**2), N.all(U == U2) -def add_jitter(psr, ecorr, flagid=None, flags=None, coarsegrain=0.1, - return_vec=False, seed=None): +def add_jitter(psr, ecorr, flagid=None, flags=None, coarsegrain=0.1, return_vec=False, seed=None): """Add correlated quadrature noise of rms `ecorr` [s], with coarse-graining time `coarsegrain` [days]. Optionally take a pseudorandom-number-generator seed.""" @@ -938,6 +935,7 @@ def computeORFMatrix(psr): return ORF + def compute_epoch_ave(times, res, err, ecorr=None, dt=1.0, flags=None): """Compute epoch averaged TOAs including uncertainty @@ -965,7 +963,7 @@ def compute_epoch_ave(times, res, err, ecorr=None, dt=1.0, flags=None): bucket_ref.append(times[i]) bucket_ind.append([i]) - avetoas = N.array([N.mean(times[l]) for l in bucket_ind],'d') + avetoas = N.array([N.mean(times[l]) for l in bucket_ind], "d") if flags is not None: aveflags = N.array([flags[l[0]] for l in bucket_ind]) @@ -973,13 +971,13 @@ def compute_epoch_ave(times, res, err, ecorr=None, dt=1.0, flags=None): aveerr = N.zeros(len(bucket_ind)) averes = N.zeros(len(bucket_ind)) - for i,l in enumerate(bucket_ind): + for i, l in enumerate(bucket_ind): M = N.ones(len(l)) - C = N.diag(err[l]**2) + C = N.diag(err[l] ** 2) if ecorr is not None: C += N.ones((len(l), len(l))) * ecorr[l[0]] - avr = 1/N.dot(M, N.dot(N.linalg.inv(C), M)) + avr = 1 / N.dot(M, N.dot(N.linalg.inv(C), M)) aveerr[i] = N.sqrt(avr) averes[i] = avr * N.dot(M, N.dot(N.linalg.inv(C), res[l]))