diff --git a/nuclockutils/barycorr.py b/nuclockutils/barycorr.py index d443d85..521540e 100644 --- a/nuclockutils/barycorr.py +++ b/nuclockutils/barycorr.py @@ -165,7 +165,6 @@ def apply_clock_correction( else: hdu.header[keyname] = corrected_time - hdu.header['CREATOR'] = f'NuSTAR Clock Utils - v. {version}' hdu.header['DATE'] = Time.now().fits hdu.header['PLEPHEM'] = f'JPL-{ephem}' diff --git a/nuclockutils/clean_clock/app.py b/nuclockutils/clean_clock/app.py index eec6124..731e347 100644 --- a/nuclockutils/clean_clock/app.py +++ b/nuclockutils/clean_clock/app.py @@ -59,8 +59,8 @@ def recalc(outfile='save_all.pickle'): load_and_flag_clock_table(clockfile=CLOCKFILE, shift_non_malindi=True) table_times = temptable_raw['met'] - met_start = table_times[0] - met_stop = table_times[-1] + met_start = clock_offset_table['met'][0] + met_stop = clock_offset_table['met'][-1] + 30 clock_jump_times = np.array([78708320, 79657575, 81043985, 82055671, 293346772]) clock_jump_times += 30 # Sum 30 seconds to avoid to exclude these points @@ -74,7 +74,8 @@ def recalc(outfile='save_all.pickle'): time_resolution=10, craig_fit=False, hdf_dump_file='dump.hdf5') table_new = eliminate_trends_in_residuals( - table_new, clock_offset_table_corr, gtis) + table_new, clock_offset_table_corr, gtis, + fixed_control_points=np.arange(291e6, 295e6, 86400)) mets = np.array(table_new['met']) start = mets[0] @@ -82,6 +83,7 @@ def recalc(outfile='save_all.pickle'): good_mets = clock_offset_table['met'] < stop clock_offset_table = clock_offset_table[good_mets] + # print(clock_offset_table['met'][-10:], table_new['met'][-1]) clock_mets = clock_offset_table['met'] clock_mjds = clock_offset_table['mjd'] diff --git a/nuclockutils/diagnostics/bary_and_fold_all.py b/nuclockutils/diagnostics/bary_and_fold_all.py index 289fb65..25db873 100644 --- a/nuclockutils/diagnostics/bary_and_fold_all.py +++ b/nuclockutils/diagnostics/bary_and_fold_all.py @@ -9,6 +9,7 @@ import astropy from astropy import log from astropy.io import fits +import pint from .get_crab_ephemeris import get_crab_ephemeris from .fold_to_ephemeris import get_events_from_fits, \ @@ -96,6 +97,13 @@ def main(args=None): parser.add_argument('--plot-phaseogram', default=False, action='store_true', help='Plot the phaseogram (requires PINT)') + parser.add_argument("--expocorr", default=False, + action='store_true', + help="Calculate the exposure from (NuSTAR only, the" + "event file must be unfiltered)") + parser.add_argument("--use-standard-barycorr", default=False, + action='store_true', + help="Use standard barycorr instead of nuclockutils") args = parser.parse_args(args) @@ -135,10 +143,20 @@ def main(args=None): mkdir(outdir) if not os.path.exists(bary_file): - cmd = f'{fname} {attorb_file} -p {parfile} ' \ - f'-c {args.clockfile} -o {bary_file}' - - nubarycorr(cmd.split()) + if args.use_standard_barycorr: + model = pint.models.get_model(parfile) + ra = model.RAJ.to("degree").value + dec = model.DECJ.to("degree").value + cmd = f'barycorr {fname} {bary_file} ' \ + f'orbitfiles={attorb_file} ' \ + f'clockfile={args.clockfile} ' \ + f'refframe=ICRS ra={ra} dec={dec}' + sp.check_call(cmd.split()) + else: + cmd = f'{fname} {attorb_file} -p {parfile} ' \ + f'-c {args.clockfile} -o {bary_file}' + + nubarycorr(cmd.split()) if args.plot_phaseogram: cmd = f'photonphase {bary_file} {parfile} ' \ diff --git a/nuclockutils/diagnostics/compare_pulses.py b/nuclockutils/diagnostics/compare_pulses.py index dcfaa66..86e914d 100644 --- a/nuclockutils/diagnostics/compare_pulses.py +++ b/nuclockutils/diagnostics/compare_pulses.py @@ -78,13 +78,14 @@ def main(args=None): mjds = [] plt.figure(figsize=(6, 9)) ref_profile = None + phase_time_plot_template, prof_plot_template = None, None if args.template is not None and os.path.exists(args.template): mjd, phase, prof, _, _, period_ms = \ format_profile_and_get_phase(args.template, template=None) - phase_time_plot, prof_plot = format_for_plotting(phase, prof, period_ms) + phase_time_plot_template, prof_plot_template = format_for_plotting(phase, prof, period_ms) local_max = phase[np.argmax(prof)] * period_ms - plt.plot(phase_time_plot, prof_plot, + plt.plot(phase_time_plot_template, prof_plot_template, drawstyle='steps-mid', label=args.template, color='k') ref_profile = prof ref_max = local_max @@ -100,13 +101,18 @@ def main(args=None): local_max = phase[np.argmax(prof)] * period_ms - phases.append(phase_res * period_ms) + local_phase_time = phase_res * period_ms + phases.append(local_phase_time) phase_errs.append(phase_res_err * period_ms) mjds.append(mjd) maxs.append(local_max) plt.plot(phase_time_plot, (i + 1) * 0.2 + prof_plot, drawstyle='steps-mid', label=f, alpha=0.5, color='grey') + if phase_time_plot_template is not None: + plt.plot(phase_time_plot_template + local_phase_time, + (i + 1) * 0.2 + prof_plot_template, + drawstyle='steps-mid', label=f, alpha=0.5, color='grey') if len(files) < 2: continue for plot_shift in [0, period_ms, 2 * period_ms]: diff --git a/nuclockutils/diagnostics/fftfit.py b/nuclockutils/diagnostics/fftfit.py index 274381d..598f4d1 100644 --- a/nuclockutils/diagnostics/fftfit.py +++ b/nuclockutils/diagnostics/fftfit.py @@ -1,166 +1,204 @@ -import numpy as np -from scipy.optimize import minimize, brentq - - -def find_delay_with_ccf(amp, pha): - nh = 32 - nprof = nh * 2 - CCF = np.zeros(64, dtype=np.complex) - CCF[:nh] = amp[:nh] * np.cos(pha[:nh]) + 1.0j * amp[:nh] * np.sin(pha[:nh]) - CCF[nprof: nprof - nh: -1] = np.conj(CCF[nprof: nprof - nh: -1]) - CCF[nh // 2: nh] = 0 - CCF[nprof - nh // 2: nprof - nh: -1] = 0 - ccf = np.fft.ifft(CCF) - - imax = np.argmax(ccf.real) - cmax = ccf[imax] - shift = normalize_phase_0d5(imax / nprof) - - # plt.figure() - # plt.plot(ccf.real) - # plt.show() - # fb=np.real(cmax) - # ia=imax-1 - # if(ia == -1): ia=nprof-1 - # fa=np.real(ccf[ia]) - # ic=imax+1 - # if(ic == nprof): ic=0 - # fc=np.real(ccf[ic]) - # if ((2*fb-fc-fa) != 0): - # shift=imax+0.5*(fa-fc)/(2*fb-fc-fa) - # shift = normalize_phase_0d5(shift / nprof) - return shift - - -def best_phase_func(tau, amp, pha, ngood=20): - # tau = params['tau'] - # good = slice(1, idx.size // 2 + 1) - good = slice(1, ngood + 1) - idx = np.arange(1, ngood + 1, dtype=int) - res = np.sum(idx * amp[good] * np.sin(-pha[good] + TWOPI * idx * tau)) - # print(tau, res) - return res - - -TWOPI = 2 * np.pi - -def chi_sq(b, tau, P, S, theta, phi, ngood=20): - # tau = params['tau'] - # good = slice(1, idx.size // 2 + 1) - good = slice(1, ngood + 1) - idx = np.arange(1, ngood + 1, dtype=int) - angle_diff = phi[good] - theta[good] + TWOPI * idx * tau - exp_term = np.exp(1.0j * angle_diff) - - to_square = P[good] - b * S[good] * exp_term - res = np.sum((to_square * to_square.conj())) - - return res.real - - -def chi_sq_alt(b, tau, P, S, theta, phi, ngood=20): - # tau = params['tau'] - # good = slice(1, idx.size // 2 + 1) - good = slice(1, ngood + 1) - idx = np.arange(1, ngood + 1, dtype=int) - angle_diff = phi[good] - theta[good] + TWOPI * idx * tau - chisq_1 = P[good] ** 2 + b**2 * S[good] ** 2 - chisq_2 = -2 * b * P[good] * S[good] * np.cos(angle_diff) - res = np.sum(chisq_1 + chisq_2) - - return res - - -def fftfit(prof, template): - """Align a template to a pulse profile. - - Parameters - ---------- - prof : array - The pulse profile - template : array, default None - The template of the pulse used to perform the TOA calculation. If None, - a simple sinusoid is used - - Returns - ------- - mean_amp, std_amp : floats - Mean and standard deviation of the amplitude - mean_phase, std_phase : floats - Mean and standard deviation of the phase - """ - prof = prof - np.mean(prof) - - nbin = len(prof) - - template = template - np.mean(template) +"""Use FFT techniques to align a template with a pulse profile. - temp_ft = np.fft.fft(template) - prof_ft = np.fft.fft(prof) - freq = np.fft.fftfreq(prof.size) - good = freq == freq +This should be accompanied by test_fftfit.py, which uses hypothesis to check +its reliability. If not, don't edit this, you don't have the git version. - P = np.abs(prof_ft[good]) - theta = np.angle(prof_ft[good]) - S = np.abs(temp_ft[good]) - phi = np.angle(temp_ft[good]) +Anne Archibald, 2020 - assert np.allclose(temp_ft[good], S * np.exp(1.0j * phi)) - assert np.allclose(prof_ft[good], P * np.exp(1.0j * theta)) +""" +from __future__ import division - amp = P * S - pha = theta - phi - - mean = np.mean(amp) - ngood = np.count_nonzero(amp >= mean) - - dph_ccf = find_delay_with_ccf(amp, pha) +import numpy as np +import scipy.optimize +import scipy.stats - idx = np.arange(0, len(P), dtype=int) - sigma = np.std(prof_ft[good]) - def func_to_minimize(tau): - return best_phase_func(-tau, amp, pha, ngood=ngood) +class FFTFITResult: + pass - start_val = dph_ccf - start_sign = np.sign(func_to_minimize(start_val)) - count_down = 0 - count_up = 0 - trial_val_up = start_val - trial_val_down = start_val - while True: - if np.sign(func_to_minimize(trial_val_up)) != start_sign: - best_dph = trial_val_up - break - if np.sign(func_to_minimize(trial_val_down)) != start_sign: - best_dph = trial_val_down - break - trial_val_down -= 1/nbin - count_down += 1 - trial_val_up += 1/nbin - count_up += 1 +def wrap(a): + """Wrap a floating-point number or array to the range -0.5 to 0.5.""" + return (a + 0.5) % 1 - 0.5 - a, b = best_dph - 2 / nbin, best_dph + 2 / nbin - shift, res = brentq(func_to_minimize, a, b, full_output=True) +def zap_nyquist(profile): + if len(profile) % 2: + return profile + else: + c = np.fft.rfft(profile) + c[-1] = 0 + return np.fft.irfft(c) - nmax = ngood - good = slice(1, nmax) - big_sum = np.sum( - idx[good] ** 2 * - amp[good] * np.cos(-pha[good] + 2 * np.pi * idx[good] * -shift) +def vonmises_profile(kappa, n, phase=0): + """Generate a profile based on a von Mises distribution. + The von Mises distribution is a cyclic analogue of a Gaussian distribution. The width is + specified by the parameter kappa, which for large kappa is approximately 1/(2*pi*sigma**2). + """ + return np.diff( + scipy.stats.vonmises(kappa).cdf( + np.linspace(-2 * np.pi * phase, 2 * np.pi * (1 - phase), n + 1) ) + ) - b = np.sum(amp[good] * np.cos(-pha[good] + 2 * np.pi * idx[good] * -shift) - ) / np.sum(S[good]**2) - eshift = sigma**2 / (2 * b * big_sum) +def upsample(profile, factor): + """Produce an up-sampled version of a pulse profile. + This uses a Fourier algorithm, with zero in the new Fourier coefficients. + """ + output_len = len(profile) * factor + if output_len % 2: + raise ValueError("Cannot cope with odd output profile lengths") + c = np.fft.rfft(profile) + output_c = np.zeros(output_len // 2 + 1, dtype=complex) + output_c[: len(c)] = c * factor + output = np.fft.irfft(output_c) + assert len(output) == output_len + return output + + +def shift(profile, phase): + """Shift a profile in phase. + This is a shift towards later phases - if your profile has a 1 in bin zero + and apply a phase shift of 1/4, the 1 will now be in bin n/4. If the + profile has even length, do not modify the Nyquist component. + """ + c = np.fft.rfft(profile) + if len(profile) % 2: + c *= np.exp(-2.0j * np.pi * phase * np.arange(len(c))) + else: + c[:-1] *= np.exp(-2.0j * np.pi * phase * np.arange(len(c) - 1)) + return np.fft.irfft(c, len(profile)) + + +def irfft_value(c, phase, n=None): + """Evaluate the inverse real FFT at a particular position. + If the phase is one of the usual grid points the result will agree with + the results of `np.fft.irfft` there. + No promises if n is small enough to imply truncation. + """ + natural_n = (len(c) - 1) * 2 + if n is None: + n = natural_n + phase = np.asarray(phase) + s = phase.shape + phase = np.atleast_1d(phase) + c = np.array(c) + c[0] /= 2 + if n == natural_n: + c[-1] /= 2 + return ( + ( + c[:, None] + * np.exp(2.0j * np.pi * phase[None, :] * np.arange(len(c))[:, None]) + ) + .sum(axis=0) + .real + * 2 + / n + ).reshape(s) + + +def fftfit_full( + template, profile, compute_scale=True, compute_uncertainty=True, std=None +): + # We will upsample the cross-correlation function to ensure + # that the highest peak is not missed + upsample = 8 + + t_c = np.fft.rfft(template) + if len(template) % 2 == 0: + t_c[-1] = 0 + p_c = np.fft.rfft(profile) + if len(profile) % 2 == 0: + p_c[-1] = 0 + n_c = min(len(t_c), len(p_c)) + t_c = t_c[:n_c] + p_c = p_c[:n_c] + + ccf_c = np.conj(t_c).copy() + ccf_c *= p_c + ccf_c[0] = 0 + n_long = 2 ** int(np.ceil(np.log2(2 * (n_c - 1) * upsample))) + ccf = np.fft.irfft(ccf_c, n_long) + i = np.argmax(ccf) + assert ccf[i] >= ccf[(i - 1) % len(ccf)] + assert ccf[i] >= ccf[(i + 1) % len(ccf)] + x = i / len(ccf) + l, r = x - 1 / len(ccf), x + 1 / len(ccf) + + def gof(x): + return -irfft_value(ccf_c, x, n_long) + + res = scipy.optimize.minimize_scalar( + gof, bounds=(l, r), method="Bounded", options=dict(xatol=1e-5 / n_c) + ) + if not res.success: + raise ValueError("FFTFIT failed: %s" % res.message) + # assert gof(res.x) <= gof(x) + r = FFTFITResult() + r.shift = res.x % 1 + + if compute_scale or compute_uncertainty: + # shifted template corefficients + s_c = t_c * np.exp(-2j * np.pi * np.arange(len(t_c)) * r.shift) + assert len(s_c) == len(p_c) + n_data = 2 * len(s_c) - 1 + a = np.zeros((n_data, 2)) + b = np.zeros(n_data) + a[0, 1] = len(template) + a[0, 0] = s_c[0].real + b[0] = p_c[0].real + b[1 : len(p_c)] = p_c[1:].real + b[len(p_c) :] = p_c[1:].imag + a[1 : len(s_c), 0] = s_c[1:].real + a[len(s_c) :, 0] = s_c[1:].imag + + lin_x, res, rk, s = scipy.linalg.lstsq(a, b) + assert lin_x.shape == (2,) + + r.scale = lin_x[0] + r.offset = lin_x[1] + + if compute_uncertainty: + if std is None: + resid = ( + r.scale * shift(template, r.shift) + r.offset - profile + ) + std = np.sqrt(np.mean(resid ** 2)) + + J = np.zeros((2 * len(s_c) - 2, 2)) + J[: len(s_c) - 1, 0] = ( + -r.scale * 2 * np.pi * s_c[1:].imag * np.arange(1, len(s_c)) + ) + J[len(s_c) - 1 :, 0] = ( + r.scale * 2 * np.pi * s_c[1:].real * np.arange(1, len(s_c)) + ) + J[: len(s_c) - 1, 1] = s_c[1:].real + J[len(s_c) - 1 :, 1] = s_c[1:].imag + cov = scipy.linalg.inv(np.dot(J.T, J)) + assert cov.shape == (2, 2) + # FIXME: std is per data point, not per real or imaginary + # entry in s_c; check conversion + r.uncertainty = std * np.sqrt(len(profile) * cov[0, 0] / 2) + r.cov = cov + + return r + + +def fftfit_basic(template, profile): + """Compute the phase shift between template and profile + + We should have fftfit_basic(template, shift(template, s)) == s + """ + r = fftfit_full(template, profile, compute_scale=False, compute_uncertainty=False) + return r.shift - eb = sigma**2 / (2 * np.sum(S[good]**2)) - return b, np.sqrt(eb), normalize_phase_0d5(shift), np.sqrt(eshift) +def fftfit(profile, template): + res = fftfit_full(template, profile) + return res.scale, 0, normalize_phase_0d5(res.shift), res.uncertainty def normalize_phase_0d5(phase): diff --git a/nuclockutils/diagnostics/fftfit_save.py b/nuclockutils/diagnostics/fftfit_save.py new file mode 100644 index 0000000..274381d --- /dev/null +++ b/nuclockutils/diagnostics/fftfit_save.py @@ -0,0 +1,184 @@ +import numpy as np +from scipy.optimize import minimize, brentq + + +def find_delay_with_ccf(amp, pha): + nh = 32 + nprof = nh * 2 + CCF = np.zeros(64, dtype=np.complex) + CCF[:nh] = amp[:nh] * np.cos(pha[:nh]) + 1.0j * amp[:nh] * np.sin(pha[:nh]) + CCF[nprof: nprof - nh: -1] = np.conj(CCF[nprof: nprof - nh: -1]) + CCF[nh // 2: nh] = 0 + CCF[nprof - nh // 2: nprof - nh: -1] = 0 + ccf = np.fft.ifft(CCF) + + imax = np.argmax(ccf.real) + cmax = ccf[imax] + shift = normalize_phase_0d5(imax / nprof) + + # plt.figure() + # plt.plot(ccf.real) + # plt.show() + # fb=np.real(cmax) + # ia=imax-1 + # if(ia == -1): ia=nprof-1 + # fa=np.real(ccf[ia]) + # ic=imax+1 + # if(ic == nprof): ic=0 + # fc=np.real(ccf[ic]) + # if ((2*fb-fc-fa) != 0): + # shift=imax+0.5*(fa-fc)/(2*fb-fc-fa) + # shift = normalize_phase_0d5(shift / nprof) + return shift + + +def best_phase_func(tau, amp, pha, ngood=20): + # tau = params['tau'] + # good = slice(1, idx.size // 2 + 1) + good = slice(1, ngood + 1) + idx = np.arange(1, ngood + 1, dtype=int) + res = np.sum(idx * amp[good] * np.sin(-pha[good] + TWOPI * idx * tau)) + # print(tau, res) + return res + + +TWOPI = 2 * np.pi + +def chi_sq(b, tau, P, S, theta, phi, ngood=20): + # tau = params['tau'] + # good = slice(1, idx.size // 2 + 1) + good = slice(1, ngood + 1) + idx = np.arange(1, ngood + 1, dtype=int) + angle_diff = phi[good] - theta[good] + TWOPI * idx * tau + exp_term = np.exp(1.0j * angle_diff) + + to_square = P[good] - b * S[good] * exp_term + res = np.sum((to_square * to_square.conj())) + + return res.real + + +def chi_sq_alt(b, tau, P, S, theta, phi, ngood=20): + # tau = params['tau'] + # good = slice(1, idx.size // 2 + 1) + good = slice(1, ngood + 1) + idx = np.arange(1, ngood + 1, dtype=int) + angle_diff = phi[good] - theta[good] + TWOPI * idx * tau + chisq_1 = P[good] ** 2 + b**2 * S[good] ** 2 + chisq_2 = -2 * b * P[good] * S[good] * np.cos(angle_diff) + res = np.sum(chisq_1 + chisq_2) + + return res + + +def fftfit(prof, template): + """Align a template to a pulse profile. + + Parameters + ---------- + prof : array + The pulse profile + template : array, default None + The template of the pulse used to perform the TOA calculation. If None, + a simple sinusoid is used + + Returns + ------- + mean_amp, std_amp : floats + Mean and standard deviation of the amplitude + mean_phase, std_phase : floats + Mean and standard deviation of the phase + """ + prof = prof - np.mean(prof) + + nbin = len(prof) + + template = template - np.mean(template) + + temp_ft = np.fft.fft(template) + prof_ft = np.fft.fft(prof) + freq = np.fft.fftfreq(prof.size) + good = freq == freq + + P = np.abs(prof_ft[good]) + theta = np.angle(prof_ft[good]) + S = np.abs(temp_ft[good]) + phi = np.angle(temp_ft[good]) + + assert np.allclose(temp_ft[good], S * np.exp(1.0j * phi)) + assert np.allclose(prof_ft[good], P * np.exp(1.0j * theta)) + + amp = P * S + pha = theta - phi + + mean = np.mean(amp) + ngood = np.count_nonzero(amp >= mean) + + dph_ccf = find_delay_with_ccf(amp, pha) + + idx = np.arange(0, len(P), dtype=int) + sigma = np.std(prof_ft[good]) + + def func_to_minimize(tau): + return best_phase_func(-tau, amp, pha, ngood=ngood) + + start_val = dph_ccf + start_sign = np.sign(func_to_minimize(start_val)) + + count_down = 0 + count_up = 0 + trial_val_up = start_val + trial_val_down = start_val + while True: + if np.sign(func_to_minimize(trial_val_up)) != start_sign: + best_dph = trial_val_up + break + if np.sign(func_to_minimize(trial_val_down)) != start_sign: + best_dph = trial_val_down + break + trial_val_down -= 1/nbin + count_down += 1 + trial_val_up += 1/nbin + count_up += 1 + + a, b = best_dph - 2 / nbin, best_dph + 2 / nbin + + shift, res = brentq(func_to_minimize, a, b, full_output=True) + + nmax = ngood + good = slice(1, nmax) + + big_sum = np.sum( + idx[good] ** 2 * + amp[good] * np.cos(-pha[good] + 2 * np.pi * idx[good] * -shift) + ) + + b = np.sum(amp[good] * np.cos(-pha[good] + 2 * np.pi * idx[good] * -shift) + ) / np.sum(S[good]**2) + + eshift = sigma**2 / (2 * b * big_sum) + + eb = sigma**2 / (2 * np.sum(S[good]**2)) + + return b, np.sqrt(eb), normalize_phase_0d5(shift), np.sqrt(eshift) + + +def normalize_phase_0d5(phase): + """Normalize phase between -0.5 and 0.5 + + Examples + -------- + >>> normalize_phase_0d5(0.5) + 0.5 + >>> normalize_phase_0d5(-0.5) + 0.5 + >>> normalize_phase_0d5(4.25) + 0.25 + >>> normalize_phase_0d5(-3.25) + -0.25 + """ + while phase > 0.5: + phase -= 1 + while phase <= -0.5: + phase += 1 + return phase diff --git a/nuclockutils/diagnostics/fold_to_ephemeris.py b/nuclockutils/diagnostics/fold_to_ephemeris.py index 2956263..ec65ea9 100644 --- a/nuclockutils/diagnostics/fold_to_ephemeris.py +++ b/nuclockutils/diagnostics/fold_to_ephemeris.py @@ -11,6 +11,7 @@ from astropy import log from astropy.table import Table +from astropy.io import fits import matplotlib.pyplot as plt from pint import models @@ -51,7 +52,8 @@ def calculate_profile_from_phase(phase, nbin=256, expo=None): if expo is not None: prof_corr = prof / expo - t = Table({'phase': np.linspace(0, 1, nbin + 1)[:-1], 'profile': prof_corr, 'profile_raw': prof}) + t = Table({'phase': np.linspace(0, 1, nbin + 1)[:-1] + 0.5 / nbin, + 'profile': prof_corr, 'profile_raw': prof}) if expo is not None: t['expo'] = expo @@ -80,6 +82,113 @@ def prepare_TOAs(mjds, ephem): return toalist +# def get_expo_kristin(evfile, nbin, n_phot_max=1000000): +# with fits.open(evfile) as hdul: +# times = np.array(hdul[1].data['TIME'][:n_phot_max]) +# priors = np.array(hdul[1].data['PRIOR'][:n_phot_max]) +# mjdref = hdul[1].header['MJDREFF'] + hdul[1].header['MJDREFI'] +# +# for i in range(times.size): +# istart = np.floor((lcDT[i] - evt[i].prior) / (nbin * p0)) +# iend = np.floor(lcDT[i] / (nbin * p0)) +# if (istart lt 0) then begin +# livetime[0:iend] += 1 +# livetime[200 + istart: *] += 1 +# +# +# endif +# if (istart ge 0) then begin +# if (iend eq 200) then begin +# iend = 199 +# livetime[0] += 1 +# endif +# livetime[istart:iend] += 1 +# endif +# endfor +def get_expo_corr(evfile, parfile, nbin=128, n_phot_max=None): + cache_file = os.path.basename( + os.path.splitext(evfile)[0]) + f'_expo_{nbin}.pickle' + if cache_file is not None and os.path.exists(cache_file): + log.info("Using cached exposure information") + return pickle.load(open(cache_file, 'rb')) + if n_phot_max is None: + n_phot_max = 500 * nbin + + # ephem = get_model(parfile) + + with fits.open(evfile) as hdul: + times = np.array(hdul[1].data['TIME'][:n_phot_max]) + priors = np.array(hdul[1].data['PRIOR'][:n_phot_max]) + mjdref = hdul[1].header['MJDREFF'] + hdul[1].header['MJDREFI'] + + hdul.close() + + expo = get_exposure_per_bin(times, priors, parfile, + nbin=nbin, cache_file=cache_file, mjdref=mjdref) + return expo + + +def get_exposure_per_bin(event_times, event_priors, parfile, + nbin=16, cache_file=None, mjdref=0): + """ + Examples + -------- + # >>> event_times = np.arange(1, 10, 0.01) + # >>> event_priors = np.zeros_like(event_times) + # >>> event_priors[0] = 1 + # >>> prof = get_exposure_per_bin(event_times, event_priors, 1/10, nbin=10) + # >>> prof[0] + # 0 + # >>> np.allclose(prof[1:], 1) + # True + """ + log.info("Calculating exposure") + + if cache_file is not None and os.path.exists(cache_file): + log.info("Using cached exposure information") + return pickle.load(open(cache_file, 'rb')) + + start_live_time = event_times - event_priors + + m = get_model(parfile) + + sampling_time = 1 / m.F0.value / nbin / 7.1234514351515132414 + chunk_size = np.min((sampling_time * 50000000, 1000)) + mjdstart = (event_times[0] - 10) / 86400 + mjdref + mjdstop = (event_times[-1] + 10) / 86400 + mjdref + + phase_correction_fun = get_phase_from_ephemeris_file(mjdstart, mjdstop, + parfile) + + prof_all = 0 + prof_dead = 0 + for start_time in tqdm.tqdm(np.arange(start_live_time[0], event_times[-1], chunk_size)): + sample_times = np.arange(start_time, start_time + chunk_size, sampling_time) + + idxs_live = np.searchsorted(start_live_time, sample_times, side='right') + idxs_dead = np.searchsorted(event_times, sample_times, side='right') + + dead = idxs_live == idxs_dead + + sample_times = sample_times / 86400 + mjdref + # phases_all = _fast_phase_fddot(sample_times, freq, fdot, fddot) + # phases_dead = _fast_phase_fddot(sample_times[dead], freq, fdot, fddot) + phases_all = phase_correction_fun(sample_times) + # phases_dead = phase_correction_fun(sample_times[dead]) + phases_all = phases_all - np.floor(phases_all) + phases_dead = phases_all[dead] + + prof_all += histogram(phases_all.astype(np.double), bins=nbin, ranges=[0, 1]) + prof_dead += histogram(phases_dead.astype(np.double), bins=nbin, ranges=[0, 1]) + + expo = (prof_all - prof_dead) / prof_all + + if cache_file is not None: + pickle.dump(expo, open(cache_file, 'wb')) + + return expo + + def get_phase_from_ephemeris_file(mjdstart, mjdstop, parfile, ntimes=1000, ephem="DE405", return_pint_model=False): @@ -134,6 +243,9 @@ def main(args=None): help="Maximum energy of photons, in keV") parser.add_argument('--n-phot-max', default=None, type=int, help="Maximum number of photons") + parser.add_argument("--expocorr", default=None, type=str, + help="Unfiltered file to calculate the exposure from " + "(NuSTAR only)") args = parser.parse_args(args) @@ -148,7 +260,22 @@ def main(args=None): n_phot_max = n_phot_per_bin * nbin ephem = get_ephemeris_from_parfile(args.parfile) - + expo = None + if args.expocorr is not None: + if args.expocorr.endswith('pickle'): + expo = pickle.load(open(args.expocorr, 'rb')) + if expo.size < nbin: + log.warning("Interpolating exposure. Use at your own risk.") + else: + expo = get_expo_corr( + args.expocorr, args.parfile, + nbin=nbin * 4, n_phot_max=300 * 4 * nbin) + if expo.size != nbin: + expo_fun = interp1d( + np.linspace(0, 1, expo.size + 1) + 0.5 / expo.size, + np.concatenate((expo, [expo[0]]))) + expo = expo_fun( + np.linspace(0, 1, nbin + 1)[:nbin] + 0.5 / nbin) log.info(f"Loading {evfile}") events = get_events_from_fits(evfile) @@ -166,7 +293,8 @@ def main(args=None): mjdstop = events_time[-1] + 0.01 log.info("Calculating correction function...") - correction_fun = get_phase_from_ephemeris_file(mjdstart, mjdstop, args.parfile) + correction_fun = \ + get_phase_from_ephemeris_file(mjdstart, mjdstop, args.parfile) del events @@ -179,11 +307,13 @@ def main(args=None): phase = correction_fun(local_events) - t = calculate_profile_from_phase(phase, nbin=nbin) + t = calculate_profile_from_phase(phase, nbin=nbin, expo=expo) label = f'{event_start}-{event_start + n_phot_max}' if emin is not None or emax is not None: label += f'_{emin}-{emax}keV' + if args.expocorr: + label += '_deadcorr' for attr in ['F0', 'F1', 'F2']: t.meta[attr] = getattr(ephem, attr).value @@ -193,9 +323,18 @@ def main(args=None): t.write(os.path.splitext(evfile)[0] + f'_{label}.ecsv', overwrite=True) fig = plt.figure() - plt.plot(t['phase'], t['profile'], label='corrected') - - plt.legend() + plt.title(f"MJD {t.meta['mjd']}") + phase = np.concatenate((t['phase'], t['phase'] + 1)) + profile = np.concatenate((t['profile'], t['profile'])) + plt.plot(phase, profile, label='corrected') + if args.expocorr: + profile_raw = np.concatenate((t['profile_raw'], t['profile_raw'])) + expo = np.concatenate((t['expo'], t['expo'])) + + plt.plot(phase, profile_raw, alpha=0.5, label='raw') + plt.plot(phase, expo * np.max(t['profile']), + alpha=0.5, label='expo') + plt.legend() plt.savefig(os.path.splitext(evfile)[0] + f'_{label}.png', overwrite=True) plt.close(fig) diff --git a/nuclockutils/diagnostics/get_crab_ephemeris.py b/nuclockutils/diagnostics/get_crab_ephemeris.py index 250b769..86a9e2d 100644 --- a/nuclockutils/diagnostics/get_crab_ephemeris.py +++ b/nuclockutils/diagnostics/get_crab_ephemeris.py @@ -52,7 +52,7 @@ def get_best_ephemeris(MJD): def get_crab_ephemeris(MJD, fname=None): - log.info("Getting correct ephemeris") + log.info(f"Getting correct ephemeris at MJD {MJD}") ephem = get_best_ephemeris(MJD) if fname is None: @@ -72,22 +72,20 @@ def get_crab_ephemeris(MJD, fname=None): return ephem -def main(evfile): - with fits.open(evfile) as hdul: - header = hdul[1].header - t0 = high_precision_keyword_read(header, 'TSTART') - t1 = high_precision_keyword_read(header, 'TSTOP') - mjdref = high_precision_keyword_read(header, 'MJDREF') - MJD = (t1 + t0) / 2 / 86400 + mjdref +def main(args=None): + import argparse + parser = argparse.ArgumentParser(description="") - get_crab_ephemeris( - MJD, fname=os.path.splitext(os.path.basename(evfile))[0] + '.par') + parser.add_argument("fnames", help="Input file names", nargs='+') + args = parser.parse_args(args) + for evfile in args.fnames: + with fits.open(evfile, memmap=False) as hdul: + header = hdul[1].header + t0 = high_precision_keyword_read(header, 'TSTART') + t1 = high_precision_keyword_read(header, 'TSTOP') + mjdref = high_precision_keyword_read(header, 'MJDREF') + MJD = (t1 + t0) / 2 / 86400 + mjdref -if __name__ == '__main__': - for fname in sys.argv[1:]: - print(f"Looking for crab ephemeris close to {fname}") - try: - main(fname) - except FileNotFoundError: - print("File not found") + get_crab_ephemeris( + MJD, fname=os.path.splitext(os.path.basename(evfile))[0] + '.par') diff --git a/nuclockutils/nustarclock.py b/nuclockutils/nustarclock.py index 57f8218..0b3a2f5 100644 --- a/nuclockutils/nustarclock.py +++ b/nuclockutils/nustarclock.py @@ -117,12 +117,17 @@ def load_and_flag_clock_table(clockfile="latest_clock.dat", shift_non_malindi=Fa return clock_offset_table -def spline_detrending(clock_offset_table, temptable, outlier_cuts=None): +def spline_detrending(clock_offset_table, temptable, outlier_cuts=None, + fixed_control_points=None): tempcorr_idx = np.searchsorted(temptable['met'], clock_offset_table['met']) + tempcorr_idx[tempcorr_idx >= temptable['met'].size] = \ + temptable['met'].size - 1 clock_residuals = \ np.array(clock_offset_table['offset'] - temptable['temp_corr'][tempcorr_idx]) + clock_mets = clock_offset_table['met'] + if outlier_cuts is not None: log.info("Cutting outliers...") better_points = np.array(clock_residuals == clock_residuals, @@ -134,11 +139,17 @@ def spline_detrending(clock_offset_table, temptable, outlier_cuts=None): i]) | ((clock_residuals[better_points] - mm[better_points]) < outlier_cuts[0]) better_points[better_points] = ~wh + # Eliminate too recent flags, in the last month of solution. + one_month = 86400 * 30 + do_not_flag = clock_mets > clock_mets.max() - one_month + better_points[do_not_flag] = True + clock_offset_table = clock_offset_table[better_points] clock_residuals = clock_residuals[better_points] - detrend_fun = spline_through_data(clock_offset_table['met'], - clock_residuals, downsample=20) + detrend_fun = spline_through_data( + clock_offset_table['met'], clock_residuals, downsample=20, + fixed_control_points=fixed_control_points) r_std = residual_roll_std( clock_residuals - detrend_fun(clock_offset_table['met'])) @@ -156,18 +167,25 @@ def spline_detrending(clock_offset_table, temptable, outlier_cuts=None): def eliminate_trends_in_residuals(temp_table, clock_offset_table, - gtis, debug=False): + gtis, debug=False, + fixed_control_points=None): - good = clock_offset_table['met'] < np.max(temp_table['met']) - clock_offset_table = clock_offset_table[good] + # good = clock_offset_table['met'] < np.max(temp_table['met']) + # clock_offset_table = clock_offset_table[good] temp_table['temp_corr_raw'] = temp_table['temp_corr'] tempcorr_idx = np.searchsorted(temp_table['met'], clock_offset_table['met']) + tempcorr_idx[tempcorr_idx == temp_table['met'].size] = \ + temp_table['met'].size - 1 clock_residuals = \ clock_offset_table['offset'] - temp_table['temp_corr'][tempcorr_idx] + # Only use for interpolation Malindi points; however, during the Malindi + # problem in 2013, use the other data for interpolation but subtracting + # half a millisecond + use_for_interpol, bad_malindi_time = \ get_malindi_data_except_when_out(clock_offset_table) @@ -241,8 +259,9 @@ def eliminate_trends_in_residuals(temp_table, clock_offset_table, # print(f'df/f = {(p(stop) - p(start)) / (stop - start)}') - btis = np.array( - [[g0, g1] for g0, g1 in zip(gtis[:-1, 1], gtis[1:, 0])]) + bti_list = [[g0, g1] for g0, g1 in zip(gtis[:-1, 1], gtis[1:, 0])] + bti_list += [[gtis[-1, 1], clock_offset_table['met'][-1] + 10]] + btis = np.array(bti_list) # Interpolate the solution along bad time intervals for g in btis: @@ -251,16 +270,33 @@ def eliminate_trends_in_residuals(temp_table, clock_offset_table, temp_idx_start, temp_idx_end = \ np.searchsorted(temp_table['met'], g) - if temp_idx_end - temp_idx_start == 0: + if temp_idx_end - temp_idx_start == 0 and \ + temp_idx_end < len(temp_table): continue table_new = temp_table[temp_idx_start:temp_idx_end] cl_idx_start, cl_idx_end = \ np.searchsorted(clock_offset_table['met'], g) local_clockoff = clock_offset_table[cl_idx_start - 1:cl_idx_end + 1] + clock_off = local_clockoff['offset'] + clock_tim = local_clockoff['met'] + last_good_tempcorr = temp_table['temp_corr'][temp_idx_start - 1] - next_good_tempcorr = temp_table['temp_corr'][temp_idx_end + 1] last_good_time = temp_table['met'][temp_idx_start - 1] - next_good_time = temp_table['met'][temp_idx_end + 1] + if temp_idx_end < temp_table['temp_corr'].size: + next_good_tempcorr = temp_table['temp_corr'][temp_idx_end + 1] + next_good_time = temp_table['met'][temp_idx_end + 1] + clock_off = np.concatenate( + ([last_good_tempcorr], clock_off, [next_good_tempcorr])) + clock_tim = np.concatenate( + ([last_good_time], clock_tim, [next_good_time])) + else: + clock_off = np.concatenate( + ([last_good_tempcorr], clock_off)) + clock_tim = np.concatenate( + ([last_good_time], clock_tim)) + + next_good_tempcorr = clock_off[-1] + next_good_time = clock_tim[-1] if cl_idx_end - cl_idx_start < 2: log.info("Not enough good clock measurements. Interpolating") @@ -271,23 +307,19 @@ def eliminate_trends_in_residuals(temp_table, clock_offset_table, q + (table_new['met'] - last_good_time) * m continue - clock_off = local_clockoff['offset'] - clock_tim = local_clockoff['met'] - - clock_off = np.concatenate( - ([last_good_tempcorr], clock_off, [next_good_tempcorr])) - clock_tim = np.concatenate( - ([last_good_time], clock_tim, [next_good_time])) order = np.argsort(clock_tim) clock_off_fun = interp1d( - clock_tim[order], clock_off[order], kind='linear', assume_sorted=True) + clock_tim[order], clock_off[order], kind='linear', + assume_sorted=True) table_new['temp_corr'][:] = clock_off_fun(table_new['met']) log.info("Final detrending...") + table_new = spline_detrending( clock_offset_table, temp_table, - outlier_cuts=[-0.002, -0.001, -0.0006, -0.0004]) + outlier_cuts=[-0.002, -0.001], + fixed_control_points=fixed_control_points) return table_new @@ -780,8 +812,14 @@ def __init__(self, temperature_file, mjdstart=None, mjdstop=None, else: mjdstart = mjdstart - additional_days / 2 + self.clock_offset_file = clock_offset_file + self.clock_offset_table = \ + read_clock_offset_table(self.clock_offset_file, + shift_non_malindi=True) if mjdstop is None: - mjdstop = sec_to_mjd(self.temptable['met'].max()) + last_met = max(self.temptable['met'].max(), + self.clock_offset_table['met'].max()) + mjdstop = sec_to_mjd(last_met) mjdstop = mjdstop + additional_days / 2 @@ -799,23 +837,30 @@ def __init__(self, temperature_file, mjdstart=None, mjdstop=None, self.hdf_dump_file = hdf_dump_file self.plot_file = label + "_clock_adjustment.png" - self.clock_offset_file = clock_offset_file - self.clock_offset_table = \ - read_clock_offset_table(self.clock_offset_file, - shift_non_malindi=True) self.clock_jump_times = \ np.array([78708320, 79657575, 81043985, 82055671, 293346772]) + self.fixed_control_points = np.arange(291e6, 295e6, 86400) + # Sum 30 seconds to avoid to exclude these points + # from previous interval self.gtis = find_good_time_intervals( - self.temptable, self.clock_jump_times) + self.temptable, self.clock_jump_times + 30) + + # table_new = temperature_correction_table( + # met_start, met_stop, temptable=temptable_raw, + # freqchange_file=FREQFILE, + # time_resolution=10, craig_fit=False, hdf_dump_file='dump.hdf5') + # + # table_new = eliminate_trends_in_residuals( + # table_new, clock_offset_table_corr, gtis) self.temperature_correction_data = \ temperature_correction_table( self.met_start, self.met_stop, - force_divisor=self.force_divisor, + # force_divisor=self.force_divisor, time_resolution=10, temptable = self.temptable, - hdf_dump_file=self.hdf_dump_file, + # hdf_dump_file=self.hdf_dump_file, freqchange_file=self.freqchange_file) if adjust_absolute_timing: @@ -854,7 +899,8 @@ def adjust_temperature_correction(self): self.temperature_correction_data, load_clock_offset_table( self.clock_offset_file, shift_non_malindi=True), self.gtis, - debug=False) + debug=False, + fixed_control_points=self.fixed_control_points) table_new['temp_corr_nodetrend'] = table_new['temp_corr'] table_new['temp_corr'] = table_new['temp_corr_detrend'] table_new.remove_column('temp_corr_detrend') @@ -863,7 +909,7 @@ def adjust_temperature_correction(self): self.temperature_correction_data = table_new def write_clock_file(self, filename=None, save_nodetrend=False, - shift_times=0.): + shift_times=0., highres=False): from astropy.io.fits import Header, HDUList, BinTableHDU, PrimaryHDU from astropy.time import Time @@ -876,7 +922,8 @@ def write_clock_file(self, filename=None, save_nodetrend=False, tempcorr_idx = np.searchsorted(table_new['met'], clock_offset_table['met']) - tempcorr_idx[tempcorr_idx >= table_new['met'].size] = table_new['met'].size -1 + tempcorr_idx[tempcorr_idx >= table_new['met'].size] = \ + table_new['met'].size -1 clock_residuals_detrend = clock_offset_table['offset'] - \ table_new['temp_corr'][tempcorr_idx] @@ -891,7 +938,8 @@ def write_clock_file(self, filename=None, save_nodetrend=False, clockerr = clock_err_fun(table_new['met']) - btis = np.array(list(zip(self.gtis[:-1, 1], self.gtis[1:, 0]))) + bti_list = list(zip(self.gtis[:-1, 1], self.gtis[1:, 0])) + btis = np.array(bti_list) for g in btis: start, stop = g @@ -902,7 +950,7 @@ def write_clock_file(self, filename=None, save_nodetrend=False, if temp_idx_end - temp_idx_start == 0: continue clockerr[temp_idx_start:temp_idx_end] = 0.001 - + clockerr[table_new['met'] > self.gtis[-1, 1]] = 0.001 new_clock_table = Table( {'TIME': table_new['met'], 'CLOCK_OFF_CORR': -table_new['temp_corr'] + shift_times, @@ -911,8 +959,22 @@ def write_clock_file(self, filename=None, save_nodetrend=False, table_new['met'], edge_order=2), 'CLOCK_ERR_CORR': clockerr}) - new_clock_table_subsample = new_clock_table[::100] + allmets = new_clock_table['TIME'] + + if not highres: + good_for_clockfile = np.zeros(allmets.size, dtype=bool) + good_for_clockfile[::100] = True + twodays = 86400 * 2 + for jumptime in self.clock_jump_times: + idx0, idx1 = np.searchsorted( + allmets, [jumptime - twodays, jumptime + twodays]) + # print(idx0, idx1) + good_for_clockfile[idx0:idx1:10] = True + new_clock_table_subsample = new_clock_table[good_for_clockfile] + else: + new_clock_table_subsample = new_clock_table del new_clock_table + if save_nodetrend: new_clock_table_nodetrend = Table( {'TIME': table_new['met'], @@ -922,8 +984,13 @@ def write_clock_file(self, filename=None, save_nodetrend=False, table_new['met'], edge_order=2), 'CLOCK_ERR_CORR': clock_err_fun( table_new['met'])}) - new_clock_table_subsample_nodetrend = \ - new_clock_table_nodetrend[::100] + if not highres: + new_clock_table_subsample_nodetrend = \ + new_clock_table_nodetrend[good_for_clockfile] + else: + new_clock_table_subsample_nodetrend = \ + new_clock_table_nodetrend + del new_clock_table_nodetrend t = Time.now() @@ -1084,6 +1151,15 @@ def plot_scatter(new_clock_table, clock_offset_table, shift_times=0): from bokeh.models import HoverTool yint, good_mets = interpolate_clock_function(new_clock_table, clock_offset_table['met']) + import matplotlib.pyplot as plt + plt.figure() + plt.scatter(clock_offset_table[good_mets]['met'], + clock_offset_table[good_mets]['offset'] - shift_times) + plt.scatter( + clock_offset_table[good_mets]['met'], + clock_offset_table[good_mets]['offset'] - shift_times + yint) + plt.plot(new_clock_table['TIME'], -new_clock_table['CLOCK_OFF_CORR']) + plt.show() yint = - yint clock_offset_table = clock_offset_table[good_mets] @@ -1538,6 +1614,10 @@ def main_create_clockfile(args=None): help="Save un-detrended correction in separate FITS " "extension", action='store_true', default=False) + parser.add_argument("--high-resolution", + help="Create a high-resolution file " + "(100 times larger)", + action='store_true', default=False) args = parser.parse_args(args) @@ -1548,7 +1628,8 @@ def main_create_clockfile(args=None): freqchange_file=args.frequency_changes) clockcorr.write_clock_file(args.outfile, save_nodetrend=args.save_nodetrend, - shift_times=args.shift_times) + shift_times=args.shift_times, + highres=args.high_resolution) clock_offset_table = read_clock_offset_table(args.offsets) plot = plot_scatter(Table.read(args.outfile, hdu="NU_FINE_CLOCK"), diff --git a/nuclockutils/utils.py b/nuclockutils/utils.py index 82de500..6d90c8f 100644 --- a/nuclockutils/utils.py +++ b/nuclockutils/utils.py @@ -243,7 +243,7 @@ def rolling_std(a, window, pad='center'): def spline_through_data(x, y, k=2, grace_intv=1000., smoothing_factor=0.001, - downsample=10): + downsample=10, fixed_control_points=None): """Pass a spline through the data Examples @@ -259,12 +259,14 @@ def spline_through_data(x, y, k=2, grace_intv=1000., smoothing_factor=0.001, control_points = \ np.linspace(lo_lim + 2 * grace_intv, hi_lim - 2 * grace_intv, x.size // downsample) + if fixed_control_points is not None and len(fixed_control_points) > 0: + print(f'Adding control points: {fixed_control_points}') + control_points = np.sort( + np.concatenate((control_points, fixed_control_points))) detrend_fun = LSQUnivariateSpline( x, y, t=control_points, k=k, - bbox=[lo_lim - grace_intv, hi_lim + grace_intv]) - - detrend_fun.set_smoothing_factor(smoothing_factor) + bbox=[lo_lim, hi_lim]) return detrend_fun diff --git a/setup.cfg b/setup.cfg index e335fec..37012f9 100644 --- a/setup.cfg +++ b/setup.cfg @@ -57,6 +57,8 @@ nustar_plot_diagnostics = nuclockutils.nustarclock:main_plot_diagnostics nustar_test_new_clockfile = nuclockutils.diagnostics.bary_and_fold_all:main nustar_compare_pulses = nuclockutils.diagnostics.compare_pulses:main nustar_compare_clock_files = nuclockutils.diagnostics.compare_clock_files:main +nustar_get_crab_ephemeris = nuclockutils.diagnostics.get_crab_ephemeris:main +nustar_fold_to_ephemeris = nuclockutils.diagnostics.fold_to_ephemeris:main [config.logging_helper] # Threshold for the logging messages. Logging messages that are less severe